材料热力学论文教授推荐8篇之第五篇：热-力耦合作用下复合材料的跨尺度分析
摘要：由于复合材料内部纤维与树脂的热膨胀系数差异很大，尤其是树脂性能对温度载荷较为敏感，服役时复合材料环境的高低温变化将使其热力学性能与常温状态产生较大差异。采用Maxwell本构模型，探讨了温度变化对树脂材料本构关系的影响。假设纤维为稳定材料，即其性能不随温度变化，依据复合材料细观力学理论选择六边形代表体积元为分析对象，建立了复合材料在温度载荷下热力学的本构模型。并分别讨论了温度载荷下复合材料内部纤维体分比和纤维排列方式变化对其热力学性能的影响，实现了热-力耦合作用下复合材料的跨尺度分析。
关键词：复合材料； 本构； 跨尺度； 细观力学有限元；
Abstract:
As the temperature of severing environment variating severely,there would be a great alteration in the thermodynamics properties of composites. Due to the great difference of thermal expansion coefficient between internal fiber and resin in the composites,the resin performance is more sensitive to temperature environment. In this paper,based on the Maxwell constitutive model,the relationship between resin constitutive and temperature is studied. Assuming that the fiber performance does not vary with temperature,according to the mesoscopic composite theory,a hexagon representative volume element is employed,and the finite element method（ FEM） is adopted to establish the thermodynamics constitutive model of composites under temperature. The influence of fiber arrangement and fiber content on the composites thermodynamics constitutive is discussed,respectively. The process would be applied for a multi-scale solution for composite under thermodynamics coupling loads.
Keyword:
composite; constitutive; multi-scale; micromechanics FEM;
1 前言
具有轻质高强、耐腐蚀、抗疲劳等优良特性的复合材料日益受到低温工程的青睐，随着应用范围的扩大，对其温度载荷下的热力学性能响应的研究也逐渐深入和广泛。有文献表明纤维的热力学性能较为稳定，基本不随温度变化[1-3],而高分子树脂的热力学性能可随温度产生较大变化。这是因为在低温下高分子树脂材料的分子链会发生"收缩"现象，使得分子间结合得更加致密，从而导致其性能对温度较为敏感[4,5].目前对于高性能复合材料，大多采用环氧树脂，因此温度载荷对环氧树脂热力学性能的影响开展较多。Ueki等[6]利用实验测试方法对具有不同化学结构的几种环氧树脂热力学性能进行分析，发现具有二维网状线性结构的聚合物无论在常温还是低温下都具有更好的性能。Walsh和Reed[5]分别对氰酸酯树脂和环氧树脂在295 K、76 K和4 K三种温度载荷下的力学性能进行了研究，发现其室温时的拉伸强度比低温时的低，并且模量更高的树脂体系倾向于具有更低强度。王嵘等[7]通过研究D-400/DETD/EP树脂体系的力学性能、成型工艺性，并通过增强材料的筛选，得到了环氧树脂体系在低温和室温环境下的热力学变化规律。
实验测试虽然能够获得相对准确的试验数据，但是温度载荷试验环境的创造需要耗费大量物力、财力，且试验操作和观测均充满挑战；为了获取足够的测试数据，通常需要对不同的材料、不同的温度点进行大量的试验，耗费大量人力、物力。而数值方法利用有限温度点的材料性能测量数据，通过数值分析获得相应温度域内的材料性能。目前已发展出多种数值分析方法，其中一个最有效的方法是细观力学有限元法。
基于细观力学理论的细观力学有限元法是一种实现从宏观到微观的跨尺度结构"场分析"的方法。其优点在于通过获得完整的细观物理场以反映宏观物理场的响应，实现细观结构对宏观性能影响的定量分析。随着数值分析方法的改进，对于理论分析法具有很大难度的高纤维体分比和复杂微观结构等问题，该方法都能很好地处理，并逐步发展出了多尺度有限元[8].Hou等[8]提出了一种基于双尺度渐进分析的有限元方法，用于求解局部区域内含有周期性裂缝的复合材料力学问题。范建华等[9]提出了一种通过在复合材料细观模型上施加多组特定形式的均匀边界条件，实现刚度系数求解的方法，并将其应用于求解计算复合材料模量的三维有限元模型中。吴世平等[10]利用细观力学有限元法探讨了复合材料内部纤维的体分比、截面形状以及微观排列方式等因素对其宏观等效性能的影响。崔俊芝等[11]通过对周期性复合材料的热力耦合问题采用双尺度有限元法进行分析，得到了包含物理、力学参数及热力耦合解的双尺度解析表达式，并据此发展了相应的多尺度有限元。
本文采用细观力学有限元法，分析、预测了单向复合材料IM7/977-3的宏观力学性能，并分析其各项热力学性能在温度载荷下的变化规律。
2 理论基础
通常高分子树脂材料具有黏弹性质，要研究温度载荷对其热力学性能的影响，须获得其力学性能参数与温度之间的关系。研究表明，Maxwell本构关系式可较好地反映这两者间的关系，即树脂材料弹性模量和温度载荷T之间满足如式（1）所示关系式[12]:
式中，T0为参考点温度；T为当前温度；Em0为T0点树脂的弹性模量；Em为T点树脂的弹性模量；α、β为待求的材料参数，可依据实验数据采用最小二乘法求得。
根据文献[13]研究，温度载荷下高性能环氧树脂的泊松比几乎不变。而对温度变化较为敏感的热膨胀系数，可采用最简单二次多项式拟合函数进行拟合，其表达式如下：
同理，采用最小二乘法求得系数a、b、c,从而获得环氧树脂材料在温度载荷影响下的热膨胀系数变化。
依据复合材料力学理论，温度载荷作用下的复合材料本构关系可表示为：
其中，εi、γij分别为复合材料的拉伸应变和剪切应变分量；ΔT为温差；σi、τij分别为复合材料的拉伸应力分量和剪切应力分量；αij（T）为复合材料的热膨胀系数；Cij（T）为复合材料的等效刚度阵。Cij（T）和αij（T）表示随温度变化的函数。建立随温度变化的复合材料跨尺度本构关系，就是要依据微观尺度的结构方式和组分材料获得温度载荷下的本构关系，以求得宏细观尺度下复合材料Cij（T）和αij（T）的值[14,15,16].这些参数可采用细观力学有限元法，通过对计算模型施加合理的载荷边界条件并计算其平均应力响应获得，进而得到其等效工程弹性常数。
3 结果与讨论
3.1 算例分析
依据文献[13]数据，环氧树脂977-3弹性模量随温度变化的测试结果如图1所示，其热膨胀系数随温度变化的测试结果如图2所示。
图1 温度载荷下树脂977-3弹性模量的测试及拟合结果Fig.1 Experiment and numerical results for temperature dependent elasticity modulus of epoxy 977-3
图2 随温度变化树脂977-3热膨胀系数的测试及拟合结果Fig.2 Experimental and numerical datas for temperature dependent CTE of the epoxy 977-3
依据式（1），采用最小二乘法对测试结果进行拟合求解，可求得α、β分别为2、0.718.即环氧树脂977-3随温度变化弹性模量具体的函数表达式为：
式中，T的单位为开尔文，K;Em的单位为GPa.图1为依据式（4）描绘的拟合曲线与试验曲线的对比图。
同理，依据式（2），采用最小二乘法对测试结果进行拟合求解，求得a、b、c参数分别为-0.159×10-3、0.233、1.673.最终得到环氧树脂977-3随温度变化热膨胀系数具体的函数表达式为：
图2为依据式（5）描绘的拟合曲线与测试数据的对比图。由图2可知，温度载荷下环氧树脂977-3热膨胀系数的变化规律可用式（5）较为准确地表示；对比测试结果与拟合结果表明，温度载荷下环氧树脂977-3的热膨胀系数随温度降低而降低，呈近似线性关系。当温差为-400 K,即从455 K降低到50 K时，其数值将由74.33×10-6/K降低至12.93×10-6/K,变化幅度达到-82.6%.
与树脂相比，纤维的热力学性能较为稳定，几乎不随温度变化。依据文献数据，获得石墨纤维IM7的相关热力学性能参数如表1所示。
表1 石墨纤维IM7的热力学参数Table 1 Thermodynamics properties for the graphite fiber IM7
在常规有限元中，建模、加载和求解过程都必须施加在具体的结构上，而由大量纤维与树脂通过物理方式复合而成的复合材料具有非均匀的微观结构，因而无法对复合材料结构通过微观层次上的建模而实现精确分析。基于细观力学理论的细观力学有限元采用一个代表体积元（Representative Volume Element,简称"RVE"）来近似代替典型复合材料，在细观力学有限元中，通过控制代表纤维尺寸与代表单胞尺寸比例的实现其纤维体分比的表达，即代表体积元的实际尺寸不具有直接的物理意义，只需通过比例参数保证所采用RVE的纤维体分比满足要求即可。
用于数值分析的RVE模型可通过观测复合材料的微观结构确定。本文以IM7/977-3复合材料为例，依据其SEM照片获得微观结构，发现采用纤维呈六边形排列的RVE可近似满足其实际纤维分布，如图3所示。
图3 IM7/977-3复合材料SEM结构及其六边形代表体积元Fig.3 SEM structure and hexagonal RVE of IM7/977-3 composite
依据图3所示的六边形RVE结构，建立相应的有限元分析模型，如图4所示。
图4 六边形代表体积元有限元模型Fig.4 FE model of hexagonal RVE
通常作为结构件用的复合材料具有较高的纤维体分比，为60%左右，故以具有该数值纤维体分比的复合材料为研究对象。基体为环氧树脂977-3,增强纤维为IM7,其具体的热力学性能分别如上节所述。
为了研究温度载荷下IM/977-3复合材料热力学性能的变化规律，以50 K为间隔分别对100 K、150 K、200 K、250 K、300 K、350 K和400 K这7个温度点的热力学性能采用细观力学有限元进行计算，依据计算结果研究其变化规律。其中，边界条件为周期性对称边界条件，单元为六面体结构化单元。相应的热力学性能计算结果如表2所示，为便于对比亦给出了理论分析的结果。
表2 不同温度下IM7/977-3复合材料的热力学性能Table 2 Thermodynamics properties for IM7/977-3 at different temperatures
从表2可知，温度载荷对复合材料轴向拉伸模量E1和泊松比μ12、μ13、μ23的影响较小，其变化趋势为随着温度降低而增加，但幅度极小，温度从400 K降低至100 K时其轴向拉伸模量和泊松比变化分别为0.9 GPa和0.002,增加幅度不到0.2%.出现该变化规律主要是由于在复合材料轴向拉伸载荷下，纤维为主要承载者，而树脂仅为载荷传递者，此外，从数值上看，纤维轴向的拉伸模量比树脂的要高两个数量级。所以即使温差很大（-300 K），且树脂拉伸模量增加了96.2%,但由于纤维热力学性能几乎不变，所以复合材料的轴向模量几乎不变。
复合材料的横向模量E2、E3主要受树脂性能影响，可采用与树脂拉伸模量相同的表达式对其进行表示。同理，采用最小二乘法解得温度载荷下E2、E3的函数表达式为：
从表2的数据可知，温度载荷下纤维横向拉伸模量的变化趋势为随着温度的降低而增加，且其变化速率随着温度降低而逐渐放缓，当温差达到-300K时，其变化的幅度为53.8%.产生该现象是由于从横向上看，复合材料内的增强纤维和树脂基体相当于是串联结构。依据复合材料细观力学理论[13],若非纤维体分比特别高，否则增强纤维对复合材料横向性能的提高起不了多大作用，起决定作用的主要是基体的性能，所以温度载荷下复合材料的横向性能E2、E3会产生如式（6）的变化。
同理，对树脂性能较为敏感的复合材料剪切性能G12、G13、G23变化规律可采用与树脂拉伸模量相同的分析方法，最终获得温度载荷下各剪切性能参数的函数关系表达式为：
由表2数据可知，温度载荷下IM7/977-3复合材料的剪切模量G12、G13、G23均具有与树脂弹性模量相类似的变化趋势，即随着温度的降低而增加。其中，面内剪切模量G12、G13的变化规律近似服从函数关系式（7），横向剪切模量G23的变化规律近似服从函数关系式（8），温差-300 K载荷下，其性能分别增加77%和46%.依据复合材料细观力学理论[13]可知，常规纤维体分比环境下，增强纤维的性能并不会对复合材料的剪切性能产生明显的影响。虽然树脂基体在复合材料中主要承担载荷传递作用，但对比式（7）和式（8）可知，基体性能对面内剪切模量G12、G13的影响要比对横向剪切模量G23的影响大。
由表2数据可看出，温度载荷下复合材料轴向热膨胀系数α11的变化呈抛物线变化趋势，即随着温度降低先升高而后降低，其极大值出现在室温300K附近，温差-300 K载荷下最大幅度可达到67.6%.分析其原因，是因为增强纤维IM7的纵向热膨胀系数为负值而树脂的为正值，同时温度载荷作用下树脂性能出现了较为明显的变化，使复合材料内部出现复杂的应力场分布。即在低温环境下，由于树脂模量增大，纤维与树脂间在横向上出现更大的约束能力，在树脂收缩时由于泊松效应，纤维将阻碍其纵向收缩，使得从低温升至常温时复合材料α11逐渐增大并趋近于0.而在较高温度环境下，则出现与之相反的变化趋势。
同理，对温度载荷下复合材料横向热膨胀系数α22、α33的变化趋势采用与式（2）相似的表达式进行表征，最终获得的函数表达式为：
由表2的数据可知，温度载荷下复合材料横向热膨胀系数α22、α33的变化规律为单调递增，符合式（9）的变化规律，即随着温度的增加而增加，温差变化到-300K时，降低的幅度为-60.3%.
为进一步验证本文细观力学有限元分析结果的准确性，对计算数据分别采用已发表文献中的解析模型进行验证[17].研究结果表明，对于复合材料两个面内泊松比μ12、μ13和轴向模量E1,若采用混合律模型，可获得与实验数据误差较少的预测结果，而对于复合材料横向模量E2、E3,剪切模量G12、G13、G23以及μ23,Halpin-Tsai公式的预测结果较为准确。此外，对于热膨胀系数α11、α22、α33,可采用Schapery公式进行理论预测。这些理论模型的分析结果亦列于表2中，以与细观力学有限元的分析结果进行对比。
由对比结果可知，在不同温度载荷下，两种方法的预测结果均十分接近。其中，采用混合律模型计算的轴向模量E1与有限元法分析结果几乎完全一致，其最大误差仅为0.06%;类似地，μ12、μ13的理论分析结果与有限元法的结果亦十分接近，其最大误差仅为3.08%.而采用Halpin-Tsai公式预测的温度载荷下复合材料横向模量E2、E3的剪切模量G12、G13、G23与有限元法分析结果相比，其最大误差分别为6.15%和3.24%.但是对于μ23,两种方法的预测结果出现较大差异，其最大误差达到了20.3%.分析出现该误差的原因主要是因为在Halpin-Tsai公式中考虑了纤维和树脂μ23比值的影响，而在本文中认为纤维和树脂μ23相同，为0.35,且不随温度变化。所以本文中采用Halpin-Tsai公式无法得到随温度变化的μ23.对于温度载荷下复合材料轴向热膨胀系数α11,两种方法的预测结果最大误差为12.3%,且都是Schapery结果小于有限元结果，推测其原因可能是由于分析时RVE的边界仍然偏"刚硬",因而内应力偏大；而对于横向热膨胀系数α22、α33,两种方法的预测结果最大误差为10.38%,推测其原因可能是因为Schapery公式中忽略复合材料泊松比的耦合效应造成的。
上述分析结果表明，细观力学有限元法能够通过单个RVE模型有效预测出复合材料的宏观力学性能且具有较高的精度，从而建立宏细观一体化本构分析模型，运用这一分析模型能够快速、准确地获得不同温度载荷下复合材料的热力学性能参数。
3.2 纤维排列的影响
在实际应用中，复合材料内部增强纤维的排列是随机分布的，所以可用于细观力学有限元法分析的RVE形式并不唯一。理论上，只要所选取的RVE单元满足周期性分布规律，均可应用细观力学有限元法。为分析简便，本文选取了结构最为简单的四边形RVE作为对比，分析该RVE代表复合材料在温度载荷下热力学性能的变化规律。为对比方便，纤维体分比同样选取为60%.根据周期对称边界条件建立如图5所示的四方单胞有限元模型。
图5 纤维排列及相应四方代表体积元有限元模型Fig.5 SEM and FE model for square RVE
采用与六边形RVE相同的建模和计算方法，对纤维呈正方形排列的IM7/977-3复合材料采用四边形RVE进行分析，其中纤维体分比仍然为60%,且各组分材料属性和温度载荷均与上节相同。计算得到复合材料热力学性能数据和六边形RVE的分析结果均列于表3中，采用符号S表示四边形RVE的结果，符号H表示六边形RVE的结果。为节省篇幅，仅列出数据差异较大的E2、E3,G12、G13、G23和μ23的计算结果。
表3 温度载荷下纤维四方排列复合材料的热力学性能Table 3 Thermodynamics properties of square RVE at different temperatures
对比表3中的数据可知，对于复合材料横向模量E2、E3,四边形RVE计算的结果均大于六边形RVE的计算结果，各个温度点的预测结果大0.8~1GPa.推测产生该现象的原因是由于从横向上看，复合材料的纤维与树脂呈串联模式，即其横向模量主要受各组分性能和组分比共同决定，其中树脂性能起主要作用。对比图3和图5可以看出，纤维体分比相同时，四边形RVE中纤维在横向的密度大于六边形RVE的密度，所以四边RVE的计算结果大于六边形RVE的计算结果。
对比表3中的数据可知，对于面内剪切模量G12、G13,四边形RVE的计算结果均大于六边形RVE的计算结果，且随着温差的增大，其差值逐渐减小，相差的百分比由400 K时的18.0%降至100 K时的13.7%.产生该现象是因为面内剪切模量G12、G13对横向纤维密度变化较为敏感，当纤维体分比相同时，四边形RVE中的横向纤维密度大于六边形RVE中的横向纤维密度，导致四边形RVE的计算数据偏大。此外，当温度载荷增大时，树脂的弹性模量增大，即对纤维的约束能力增大，导致纤维排列对复合材料剪切模量G12、G13的影响减弱。
对比表3中的数据可知，对于面外剪切模量G23,四边形RVE的计算结果均偏小，其差值几乎稳定为0.34 GPa.但随着温度载荷的增大，其差值百分比由400 K时的20.9%降至100 K时的8.4%.推测出现该变化是由于面外剪切模量G23主要受纤维的"环向密度",即纤维周围直接接触相关的纤维数量和树脂的性能决定。对比图3和图5的纤维排列可看出，六边形RVE的每根纤维均有6根纤维分布在其周围，比四边形纤维排列多2根。在横向剪切载荷作用下，中心的纤维将受周围的纤维保护而避免发生剪切，因而四边形RVE的面外剪切模量G23的计算结果小于六边形RVE的。
对比表3中的数据可看出，对于μ23,四边形RVE的计算结果均小于六边形RVE的，此外采用六边形RVE计算的温度载荷下μ23的变化极小，而在四边形RVE中，温度载荷下μ23的计算结果出现了增长，温差300 K时其增幅达到15.3%.分析其原因是由于与树脂的弹性模量相比，纤维的弹性模量很大，即树脂比纤维"软"很多，所以树脂是复合材料泊松效应横向收缩量的主要来源。对比图3和图5可知，四边形RVE中的横向纤维密度小于六边形RVE的，导致四边形RVEμ23的计算结果均小于六边形RVE的。
上述分析表明，六边形RVE分析模型比四边形RVE在预测对于树脂性能较为敏感的横向模量E2、E3和剪切模量G12、G13、G23时具有更高的计算精度。
3.3 纤维体分比的影响
作为承载构件使用的复合材料结构通常按照60%纤维体分比进行设计。但在实际制备生产中，由于工艺参数的波动和技术条件的限制，树脂分布并不十分均匀，局部区域容易出现贫胶或富胶现象，导致该区域复合材料热力学性能发生变化。为进一步研究温度载荷下复合材料纤维体分比变化对其热力学性能的影响，分别对纤维体分比为50%、60%、70%的复合材料进行分析，并将变化较为明显的计算结果分别列于表4中。
表4 温度载荷下纤维体分比分别为50%、60%、70%时的复合材料热力学性能Table 4 Thermodynamics properties of composites with different fiber content at different temperatures
从表4中的数据可以看出，在温度载荷下，纤维体分比不同，复合材料的热力学性能分布变化规律相似。对于复合材料的横向模量E2、E3来说，纤维体分比越高，其横向模量E2、E3数值增加越快。根据复合材料细观力学分析理论可推测出产生该变化的原因，即横向模量E2、E3的表达式为：
若将上式对纤维体分比cf求偏导数，可得：
由于纤维弹性模量Ef远大于树脂弹性模量Em,且在固定温度点下Ef-Em为定值，所以式（11）右边的分子为大于0的正值。表达式中只有一个变量cf,当cf增加时，分母减小，所以?E2/?cf随着cf的增加而增大，且呈二次方变化，即横向模量E2、E3随着纤维体分比的增加而快速增加。
同理，对影响复合材料的面内剪切模量G12、G13进行理论分析。依据复合材料细观力学可知面内剪切模量G12、G13的表达式为：
若将上式对cf求偏导数，可得：
由于树脂剪切模量Gm小于纤维剪切模量Gf,且在固定温度点下Gf-Gm为不变量，则式（13）右边的分子为大于0的正值。式（13）中只有一个变量cf,当cf增加时分母减小，所以?G12/?cf随着cf的增加而增大，且呈二次方变化，即面内剪切模量G12、G13随着纤维体分比的增加而快速增加。
由表4中的数据可以看出，温度载荷下不同纤维体分比的复合材料面外剪切模量G23均随着纤维体分比的增加而增加。产生该现象变化的原因可根据Halpin-Tsai公式推测出，已知G23的表达式为：
其中，ζ、η为与纤维排列方式、几何形状和载荷情况有关的参数，η≤1.同理，将上式对cf求偏导数，可得：
由式（14）可知，面外剪切模量G23随着纤维体分比cf的增加而增大。此外，由于η≤1,即cf增加时，1-ηcf则单调降低。从式（15）可知，cf越大，面外剪切模量G23增加越快，且为二次方关系。
对比表4中的数据可知，温度载荷下纤维体分比不同复合材料μ23的变化趋势一致，即随着纤维体分比的增加而降低，且温度载荷越高，μ23的变化越快。分析产生该现象的原因是由于复合材料的横向泊松效应主要是由树脂的泊松收缩提供的，当纤维体分比增加时，相应地，纤维对树脂的约束将增加，抑制纤维的横向收缩，导致纤维体分比越高时，μ23越低。此外，高温载荷下树脂的弹性模量降低，即由于与纤维间的性能差异增大，凸显了纤维对树脂收缩的抑制作用，从而导致温度载荷越大，μ23随着纤维体分比变化的趋势越快。
4 结论
本文采用Maxwell本构模型获得了低温下环氧树脂977-3的本构关系，进一步依据细观力学建立温度载荷下复合材料热力学本构的分析模型，并以IM7/977-3复合材料为例验证了该分析模型的准确性和有效性。通过多个算例间的数值对比和理论分析可知：
（1）本文采用的计算方法是预测复合材料热力学性能的有效方法；
（2）六边形RVE具有比四边形RVE更高的分析精度，尤其是在预测对树脂性能较为敏感的宏观性能时；
（3）纤维体分比的变化对复合材料宏观性能的影响并不呈线性关系。
参考文献
[1]Reed R P,Golda M.Cryogenic composite supports:a review of strap and strut properties[J].Cryogenics,1997,37（5）：233-250.
[2] 马帅。低温复合材料贮箱跨尺度层次化分析[D].大连：大连理工大学，2014.
[3]Schutz J B.Properties of composite materials for cryogenic applications[J].Cryogenics,1998,38（1）：3-12.
[4]王嵘，郝春功，陈振坤，等。混合型固化剂对环氧树脂室温和低温力学性能的影响[J].复合材料学报，2008（04）：7-12.
[5]尹志娟，王丽雪，姜珊。玻璃纤维增强环氧树脂基复合材料的低温性能研究[J].黑龙江工程学院学报（自然科学版），2010（01）：50-52.
[6]Ueki T,Nishijima S,Izumi Y.Designing of epoxy resin systems for cryogenic use[J].Cryogenics,2005,45（2）：141-148.
[7]王嵘，郝春功，陈振坤，等。碳纤维增强环氧树脂室温低温力学性能的研究[C]//第十七届玻璃钢/复合材料学术年会。中国广东广州，2008.
[8]Hou T Y,Wu X H,Cai Z.Convergence of a multiscale finite element method for elliptic problems with rapidly oscillating coefficients[J].Math Comp,1988,68:913-943.
[9]范建华，许庆余，孙秦，等。一种通用的计算复合材料刚度的有限元方法[J].工程力学，2002（04）：166-170.
[10]吴世平，唐绍锋，梁军，等。周期性复合材料热力耦合性能的多尺度方法[J].哈尔滨工业大学学报，2006（12）：2049-2053.
[11]曹礼群，崔俊芝。复合材料拟周期结构的均匀化方法[J].计算数学，1999,21（3）：331-344.
[12]Yamini S,Young R J.The mechanical properties of epoxy resins.Part I,mechanisms of plastic deformation[J].Journal of Materials Science,1981,15:1814-1822.
[13]Choi S.Micromechanics,fracture mechanics and gas permeability of composite laminates for cryogenic storage systems[D].Florida:University of Florida,2005.
[14]Hui D,Dutta P K.Cryogenic temperature effects on performance of polymer composites[R].2002.
[15]Yu Y H,Kim B G,Lee D G.Cryogenic reliability of composite insulation panels for liquefied natural gas（LNG）ships[J].Composite Structures,2012,94（2）：462-468.
[16]Roy S,Benjamin M.Modeling of permeation and damage in graphite/epoxy laminates for cryogenic fuel storage[J].Composites Science and Technology,2004,64（13-14）：2051-2065.
[17]沈观林，胡更开。复合材料力学[M].北京：清华大学出版社，2006.
