基于FFT方法的复合材料有效性能预报①

田晓晓，孟松鹤，矫利闯，易法军，解维华

(1.哈尔滨工业大学 复合材料与结构研究所，哈尔滨 150080;2.巴黎北大学，LPMTM，巴黎 93430)

0 引言

复合材料是由2种或多种不同性质的材料用物理和化学方法在宏观尺度上组成的具有新性能的材料［1］，具有比强度高、比刚度大、抗疲劳性能好、减震性能好和材料性能可设计等优点［2］，应用范围非常广。由于复合材料细观结构非常复杂，有效性能就比较难以获得。很多国内外学者从实验［3－4］和数值［5－6］角度进行了复合材料有效性能预报。这些数值方法往往只能对简单结构进行力学性能预报，而不能应用到复杂结构中去。

快速Fourier变换是计算离散Fourier变换的一种快速算法(简称FFT)，主要用于数字信号处理，计算大整数乘法，求解偏微分方程等。Moulinec和 Suquet［7－8］首次使用 FFT 方法提出一种数值迭代方法，用来计算脆性或准脆性材料的全局响应问题。Michel［9］等随后又用该方法给出了线性和非线性均匀复合材料的全局响应。本文从复合材料的平衡方程出发，将复合材料问题转化为含极化应力的各向同性材料问题，结合Lippmann-Schwinger方程实现其形式解，并通过FFT方法实现其数值计算，最终实现复合材料有效性能的预报，并进一步讨论了分辨率和尺寸效应对FFT方法的影响。

1 复合材料平衡方程解

考虑一个体积为Ω边界为∂Ω的周期性规则实体，可以看作是其单胞构造，也就是代表性体积单元，简写为RVE。用2个不同的空间坐标来描述:宏观尺度参数x用来描述整个实体的全局应力场;微观尺度参数y可用来描述局部应力场。其关系可定义为y=(x －)/ε，其中是单胞的坐标几何中心。经过多尺度渐进分析后，RVE内含应变梯度的平衡方程可转化为求解如下方程［10］:

式中 cijkl为复合材料弹性常数;fi为体积力;σij为应力;为剪应力;εkl为应变;#代表为周期性的;－#为反周期的。设为参考材料模量，为各向同性材料。令τij=(cijkl－)εkl+，则方程可化简为

这里可考虑更细致一点的微观结构，如考虑结构的非均匀性、界面层、纤维束的强度分布、缺陷分布等因素，可将这些性质转化为含极化应力的均匀参考材料进行研究，对研究复合材料的破坏机理奠定了良好的基础。

式(2)的解可通过周期Green算子Γ0在Fourier空间表示［8］:

其衍生问题可解决x处受指定应变E时非均匀弹性复合材料的问题:

其实空间和Fourier空间的解分别为

2 FFT运算法则

初始值:

经过多次实验，当参考材料满足下式时收敛速度最快［8］:

3 实例分析

通过单向玻璃纤维增强环氧树脂基复合材料算例来验证FFT方法的有效性，并讨论了分辨率和尺寸效应对FFT方法的影响。

单向玻璃纤维增强环氧树脂基复合材料二维单胞结构如图1所示;其组分性能［12］见表1。

表1 单向玻璃纤维增强环氧树脂基复合材料组分材料性能Table 1 Properties of components of unidirectional glass fibre/epoxy resin composite

唐邵锋等［12］曾经用多尺度法和Mori-Tanaka法计算过该材料的有效杨氏模量，并和文献［13］中的实验值进行了对比。本文将FFT的计算结果与多尺度法和实验结果进行对比。为了分析尺寸效应和分辨率对FFT方法的影响，分别建立2个模型并进行计算对比，模型及对比结果如下:

3.1 分辨率对FFT方法的影响

取1个单胞(图1)，采用FFT方法分别计算分辨率为100×100和400×400下材料的有效模量，并与其他方法和实验结果进行对比，如图2所示。

由图2可见，当纤维百分含量不超过40%时，FFT方法和多尺度方法数值计算结果几乎吻合;当纤维百分含量超过40%时，FFT和多尺度数值计算方法结果开始微有不同。通过与文献［13］中的实验值相比较，可看出FFT方法计算结果比多尺度方法更接近实验值;另外，从图2中可看出，2种分辨率计算结果随纤维百分含量的增加差别开始变大，说明分辨率对FFT方法的影响随着纤维体积分数的增大而增大。但与多尺度等其他方法相比，FFT计算方法可进一步考虑更细致的结构，如加入材料本身缺陷、界面、具有统计概率分布的纤维束强度等因素，引入这些参数对FFT计算方法来说是非常容易的，其计算结果更接近真实的实验值。由于材料固有缺陷、界面性质和纤维束的强度分布也是较难确定的，故在本算例中暂不考虑，只研究FFT方法的可行性及有效性。从计算结果可看出，该方法是可行的和有效的，且后期计算非常简洁，易于操作。

3.2 尺寸效应对FFT方法的影响

为了验证FFT方法的尺寸效应，选取一个包含若干单胞的结构，如图3所示。选取的是9单胞结构，呈3×3排列，中间单胞完整包含在整体结构内。采用FFT方法分别计算分辨率为100×100和400×400的2种情况下的有效模量，并与多尺度方法计算结果和实验结果进行对比，如图4所示，与单胞时规律相同。

为了验证FFT方法的尺寸效应，将单胞和9单胞结构计算结果进行对比。如图5所示，可看出在高分辨率下，单胞结构和9单胞结构计算结果相差很小，几乎吻合;低分辨率下，单胞结构和9单胞结构计算结果有差别，但相差很微小，几乎可忽略。这是因为此结构为周期性结构，没有考虑材料的非均匀性，从理论上来看，单胞结构和9单胞组合结构的计算结果应是完全一样的，这里的微弱差别是由于数值计算原因造成的。

4 结论

(1)将复合材料平衡方程求解问题转化为含极化应力的各向同性材料平衡方程求解问题，并结合FFT方法实现其数值解。

(2)FFT方法预报复合材料有效性能与实验结果吻合较好，与其他细观力学方法相比，FFT方法精度较高，且FFT方法可实现程序化的运算，后期计算非常简洁，易于操作。

(3)FFT数值计算时，分辨率对计算结果有影响。分辨率对FFT方法的影响随着纤维体积分数的增大而增大。

(4)FFT数值计算时，分别采用1单胞和9单胞结构进行计算，两者计算结果基本一致，说明了FFT方法不受尺寸效应影响。因此，可用单胞代替整体结构进行计算。

［1］沈观林，胡更开.复合材料力学［M］.北京:清华大学出版社，2006.

［2］罗祖道，王震鸣.复合材料力学进展［M］.北京:北京大学出版社，1992.

［3］杨成鹏，矫桂琼，王波.2D-C/SiC复合材料力学性能的实验研究［J］.复合材料:创新与可持续发展，2010(8):729.

［4］Macander A B，Crane R M.Fabrication and mechanical properties of multi-dimensionally braided composite materials［C］//Whitney J M.Composite Materials:Testing and Design(7th Conf.)，STPS893，ASTM，1986:422-443.

［5］Naik N K，Ganesh V K.An analytical method of plain weave fabric composites［J］.Composites，1995，26(4):281-289.

［6］顾震隆，吕文涛，孟兰，等.多向炭/炭复合材料弹性常数的理论预报［J］.复合材料学报，1988，5(4).

［7］Moulinec H，Suquet P.A fast numerical method for computing the linear and the nonlinear mechanical properties of composites［J］.CR Acad Sci.Paris Ser II，1994;318:1417-1423.

［8］Moulinec H，Suquet P.A numerical method for computing the overall response of nonlinear composites with complex microstructure［J］.Comput.Methods Appl.Mech.Engrg.，1998，157:69-94.

［9］Michel J C，Moulinec H，Suquet P A.Computational method based on augmented lagrangians and fast Fourier transforms for composites with high contrast［J］.Comput.Modelling Engrg.Sci.，2000，1(2):79-88.

［10］Tian Xiao-xiao，Li J，Meng Song-he，et al.A homogenized method including strain gradients based FFT prediction of the mechanical properties of composites［J］.Polymers ＆Polymer Composites，2011，19(4＆5).

［11］Kroner E.Statistical continuum mechanics［M］.Wien-New York，Springer-Verlag，1972:110.

［12］唐绍锋.复合材料热/力学性能的双尺度渐进分析［D］.哈尔滨工业大学，2009.

［13］Tsai S W.Structural behavior of composite materials［R］.NASA CR-71，1964.