二维机织复合材料力学分析中的周期性边界条件研究

薛亚红， 陈继刚， 闫世程， 骆俊廷

(1. 燕山大学 机械工程学院， 河北 秦皇岛 066004； 2. 燕山大学 自润滑关节轴承共性技术航空科技重点实验室，河北 秦皇岛 066004； 3. 燕山大学 先进锻压成形技术与科学教育部重点实验室， 河北 秦皇岛 066004)

二维机织复合材料力学分析中的周期性边界条件研究

薛亚红1，2， 陈继刚1，2， 闫世程1，2， 骆俊廷1,3

(1. 燕山大学 机械工程学院， 河北 秦皇岛 066004； 2. 燕山大学 自润滑关节轴承共性技术航空科技重点实验室，河北 秦皇岛 066004； 3. 燕山大学 先进锻压成形技术与科学教育部重点实验室， 河北 秦皇岛 066004)

为了精确地进行二维机织复合材料力学性能的数值分析,需建立单胞模型的准确边界条件。基于周期边界条件理论，提出了简便通用的二维机织复合材料周期边界方程，并给出了周期边界条件下各弹性常数在有限元分析中的求解方法；为验证周期边界条件的正确性，建立了9个单胞构成的九宫格结构，取中央单胞作为参考单胞，对不同边界条件下独立单胞的变形和应力分布与参考单胞进行对比。研究结果表明：即使在单向拉伸载荷下，单胞各个边界面也不保持平面状态，而是出现凹凸翘曲变形，即存在边界周期性；通过边界周期性条件，可正确地获得二维机织织物的工程弹性常数。

二维机织复合材料； 周期边界条件； 单胞； 有限元模型； 弹性常数

二维机织复合材料是一种新型轻质高效结构材料，具有质量轻，抗冲击性好，比强度和比刚度高等优点，被广泛应用于航空航天、国防、电力、医药和纺织等领域[1]。对于二维机织复合材料综合力学性能的预测是其用于结构设计的重要前提。织物复合材料在结构上一般呈现良好的周期性，通常基于单胞计算分析其宏观力学性能[2]，主要包括理论解析法和有限元计算法。理论解析法基于等应力或应变假设，按所占的体积分数加权平均组分材料刚度矩阵来得到织物总体的等效刚度。它能很好地预测织物的宏观力学性能，但不能得到材料细观层次上的应力分布。有限元计算法通过建立织物的单胞模型，并施加一定的边界条件，即可得到材料细观层次的应力应变分布及宏观力学性能，因此受到众多学者的青睐。

基于单胞的织物有限元分析中，合理的边界条件是保证计算准确的关键。目前普遍假设边界面保持一致[3]，即约束一个面的法向，对其对立面施加荷载，而其他面进行法向耦合使其保持平面。这种边界条件在有限元计算中很方便施加，但形成了过约束，破坏了相邻单胞边界面上应力的连续性，仅适用于单向载荷状态具有反对称结构的单胞[4-5]。Hori等[6]研究结果表明：等应变假设的边界条件将得到织物弹性常数的上限，应力连续性难以满足；等应力假设的边界条件将得到织物弹性常数的下限，不能满足相邻单胞变形的连续。为此，文献[7]通过建立多个单胞来减小边界效应，虽然能获得较好的结果，但大大降低了计算效率。

Whitcomb[8]、Xia[9]提出了周期边界条件理论，较精确地获得了织物的力学性能。Li[10]、Xu[11]、Carvalho[12]根据所研究织物的具体结构，给出了各种单胞下周期边界的表达形式，虽然较好地预测了织物的力学性能，但这些表达式均需事先给定单胞的平均应变，通过有限元模拟得到相应的应力。而工程上通常在已知载荷的情况下，用简单拉伸及剪切实验测量试件的位移，求得工程弹性常数。另外，周期边界方程依赖于单胞的选取、织物的结构及受载情况。文献[13]利用对称和反对称性，将单胞简化为1/16单胞模型，在降低计算量的同时，极大地增加了周期性边界条件的施加难度；张超等[14]介绍了三维织物复合材料的周期边界方程及弹性常数求解方法，但并不适用于只具有2对周期边界的二维机织物。文献[15]将三维周期性边界表达式应用在具有二维机织复合材料上，将非周期性边界面施加周期边界方程，与实际并不相符。

为此，本文给出了一种对于二维机织复合材料通用而简便的周期边界方程，通过建立多个单胞组成的九宫格结构，取中间单胞以减小边界效应的方法，验证周期边界方程的准确性。

1 周期边界条件理论

织物复合材料在结构上一般具有周期性，只需进行平移而不需要旋转就能形成整个织物结构的最小重复单元，称为单胞。图1示出具有周期性结构的织物胞体二维示意图。

图1 具有周期性结构的织物胞体二维示意图

胞体内任意点的位移场[9]可表示为

(1)

(2)

(3)

式中：k+和k-分别表示法向沿Xk正负2个方向的周期边界面。

(4)

(5)

从式(5)可以看出，方程满足边界上位移的周期性与连续性，并且平行的周期边界面上对应节点的位移差保持为常数。Xia等[9]还进一步证明施加式(5)后相邻单胞边界处能够同时满足应力的连续性，因此可采用式(5)对单胞进行分析。

2 周期边界条件有限元分析模型

2.1 单胞分析模型

斜纹织物衬垫为典型的具有二维周期性结构的机织复合材料。图2示出一上二下左斜纹织物的细观结构三视图。图中虚线部分为其单胞，沿着2个编织方向，单胞结构上具有周期性。为分析方便，建立了织物单胞坐标位置图，如图3所示。单胞在坐标系中位置为0≤x≤2a，0≤y≤2b，0≤z≤2c，坐标原点建立在点D，2对周期边界面分别沿x轴和z轴，织物厚度方向沿着y轴，结构上不具有周期性。

图2 斜纹织物的细观结构三视图

图3 织物单胞坐标位置图

2.2 有限元力学分析的周期边界条件设定

2.2.1 刚体位移约束

为了消除单胞刚体位移，选择顶点A、D、C、H，约束式(6)中的6个自由度。

(6)

式中：U、V、W分别表示沿x、y、z方向的位移；下标字

母表示顶点处的节点；下标数字表示棱边上的节点；下标为等式则表示对应面上的节点。

2.2.2 面上节点约束

在2对平行的周期边界面中，分别选择x=0与z=0的面为主平面，与主平面相对的面为从平面，主平面上的节点为主节点，从平面上的节点为从节点。选择节点对(D，C)，(D，A)作为参考，由式(5)可知，对于法向沿x方向的对应主节点与从节点，其位移差即为节点C与D的位移差。由于节点D固定，位移差即为节点C的位移。同理，法向沿z方向的主节点与从节点的位移差，即节点A的位移。

法向沿x轴的面内节点约束：

(7)

法向沿z轴的面内节点约束：

(8)

连续体沿y轴方向的结构不具有周期性，因此，法向沿y轴的面不需要施加周期边界约束。

2.2.3 棱边节点约束

对于自由的4个顶点(B、E、F、G)和4个棱边(9，10，11，12)，处于两周期边界面的交界处，其节点同时满足式(7)、(8)，如果直接把这些约束都施加到有限元分析中，就会导致有些节点既为主节点，又为从节点，计算将无法进行，因此，需将这些棱边和顶点单独施加相互独立的方程。

以棱9为参考，对于棱10，只施加方程(7)，其节点约束为

(9)

对于棱12，只施加方程(8)，其节点约束为

(10)

对于棱11，由式(7)、(8)可得其节点约束为

(11)

(12)

2.2.4 顶点约束

以顶点D为参考，对于顶点G，只施加方程(7)，其节点约束为

(13)

对于顶点E，只施加方程(8)，其节点约束为

(14)

对于顶点F，由式(7)、(8)可得其节点约束为

(15)

以顶点D为参考，对于顶点B，其节点约束为

(16)

式(6)～(16)，可应用有限元分析软件中的边界设定，利用其多点约束功能并结合编写的相应程序加以实现。

3 周期边界下织物弹性常数求解方法

3.1 宏观弹性常数求解方法

在利用细观有限元法研究二维机织复合材料宏观力学性能时，施加周期边界条件的目的，仅仅保证

可以用1个单胞代替整个织物的建模，如通过单胞研究织物的宏观弹性常数进行强度预测等。欲求得织物的宏观力学性能，需要在周期边界条件的基础上，施加一定的载荷。二维机织复合材料为各项异性材料，其等效本构方程可记为

(17)

(18)

式中：Ex、Ey、Ez分别为二维机织复合材料在x、y、z方向的弹性模量；Gxy、Gyz、Gxz分别为xy、yz、xz面内的剪切模量；μxy、μyz、μxz、μyx、μzy、μzx分别为xy、yz、xz、yx、zy、zx面内的泊松比。若要得到织物的9个工程弹性常数，可以施加6组线性不相关的宏观应力场。周期边界条件下的6组简单宏观应力场如表1所示。

表1 周期边界条件下的6组简单宏观应力场Tab.1 Six groups of simple macroscopic stress fieldparameters under periodic boundary condition

施加表1所示的1种应力场，可得到1个或多个二维机织复合材料的宏观弹性常数。施加全部6种应力场，即可得到织物全部9个宏观弹性常数。如施加k=1的应力场，可求得Ex、μxy、μxz。

(19)

3.2 周期边界条件下的加载方法

对单胞施加周期边界方程后，对于表1中的6组宏观应力场，只需分别对6个节点自由度施加相应的集中力。

(20)

式中Fx|C为在顶点C施加x轴方向的集中力。

通过有限元后处理，由节点位移求得连续体宏观平均应变值。

(21)

周期边界条件下，二维机织复合材料宏观弹性常数的具体求解方法为：首先在单胞模型的周期边界面施加周期边界方程(6)～(16)，然后给定表1所示的其中1种或几种宏观应力场，按照式(20)的载荷施加方式对节点施加集中力，计算后通过式(21)获得宏观平均应变值，最后由式(19)获得织物宏观弹性常数。

由于单胞y=0及y=2b的2面为非周期边界面，未对节点进行关联约束，因此，对于k=2及k=5的应力场，为避免应力集中，在施加集中力前，先将2面内的节点进行耦合处理。对于k=2的应力场：需将面y=2b内部(不包括4个棱边和4个顶点)的所有节点与顶点H进行y方向的耦合，将面y=0内部的所有节点与顶点D进行y方向的耦合。对于k=5的应力场：需要将面y=2b内部的所有节点与顶点H进行z方向的耦合；将面y=0内部的所有节点与顶点D进行z方向的耦合。

4 实例验证

4.1 周期边界效应存在性验证

为验证周期边界效应的存在，进行了斜纹织物衬垫压缩实验。斜纹织物衬垫由美国杜邦公司提供，其由径向芳纶和纬向PTFE纤维编织而成，基体为酚醛树脂。试件尺寸为100 mm×50 mm×0.4 mm，压缩实验在TOX电子压机上进行，压缩方向为垂直纤维编织方向，载荷为75 MPa。图4示出显微镜下斜纹织物衬垫压缩前后细观结构图。可以看出，斜纹织物衬垫经压缩后，由于斜纹织物的细观组织结构的不对衬，单胞的边界面变形不一致，明显看出径向纤维束由直线状经过压缩出现屈曲变形，且这种变形形态存在重复性和周期性。

图4 二维织物衬垫的受压前与受压后对比图

4.2 周期边界条件设置的正确性验证

为验证周期边界条件设置的正确性，以斜纹织物为对象，进行力学性能的有限元模拟及实验。其试件材料及尺寸同衬垫压缩实验，力学性能实验在微机控制电子式万能试验机上进行。图5示出斜纹织物的单胞的几何模型。建立了如图6所示的由9个单胞组成的九宫格结构，其中间单胞边界效应较小，作为参考单胞。

以施加x方向的单轴拉伸为例，对于不同的胞体载荷施加方式为：1) 按照边界保持一致的假设对九宫格结构施加约束。即固定ADEH面，使其法向位移为0；对CDHG、ABFE这2个边界面施加法向耦合约束，使其法向变形保持平面；在BCGF面上施加沿x方向、大小为150 N的单向拉力。2) 按照本文方法对单胞施加周期边界方程后，在单胞的节点C施加大小为50 N的集中力。3) 按照边界保持一致的假设对单胞施加约束后，在BCGF面上施加沿x方向大小为50 N的单向拉力。图7示出不同边界下胞体的变形形态。

图5 斜纹织物单个胞体的几何模型

图6 九宫格结构中参考单胞分析示意图

图7 不同边界条件下胞体的等效应力图

从图7可看出，由于边界效应的影响，九宫格结构中沿x轴方向排布的3个单胞应力分布各不相同。周期边界条件下的单胞变形形态与参考单胞的基本一致。此外，由于结构的不对称性，即使在单向应力状态下的边界面也不再保持平面，而是出现了翘曲，但平行的周期边界面保持相同的形状变形，其位移差相等，从而保证了相邻单胞位移的连续性。周期边界下的单胞与参考单胞，应力分布有较小偏差，这是由于对九宫格结构施加了简化的边界条件，参考单胞受到一定的边界效应的影响，并不能模拟连续体真实的周期性状态。从图5可看出，在平面保持一致的假设下单胞的等效应力整体偏大，这是由于边界面保持一致，相当于施加了附加力，一定的载荷下变形减弱，使得单胞整体刚度增加。

表2给出了不同边界条件下有限元模拟与实验获得的部分弹性常数值。可以看出，周期边界条件下的弹性常数值与实验结果较为接近，其预测的斜纹织物弹性模量误差能控制在6%以内；由于单胞边界面在变形过程中被迫保持为平面状态，使得边界保持一致假设下单胞的弹性常数预测值偏大，这与前面的分析一致。因此，基于本文给出的周期边界条件的设置方法对预测二维织物复合材料的弹性常数有足够的精度，及相对于传统边界保持一致的边界条件设置方法具有绝对的优势。

表2 不同边界条件下有限元模拟值与实验结果比较

5 结 论

本文针对二维机织复合材料的有限元分析模型问题，开展了单胞模型周期边界条件的研究，得出了具有二维周期性结构的机织复合材料单胞的周期边界方程，给出了在有限元中织物宏观弹性常数求解方法。通过斜纹织物压缩实验，证实了织物周期边界效应的存在。将九宫格结构的中央单胞作为参考单胞，对比不同边界条件下单胞变形和应力分布，并进行斜纹织物弹性常数预测的模拟和实验，证明了本文给出的单胞模型周期边界条件的正确性。研究表明，织物复合材料的单胞即使在单向拉伸载荷下各个边界面也不保持平面状态，而是出现凹凸翘曲变形，即存在边界周期性；本文给出的边界周期性条件，相对于传统有限元方法预测二维机织复合材料弹性常数有更高的精度。

FZXB

[1] QIU M, MIAO Y W, LI Y C, et al. Film-forming mechanis ms for self-lubricating radial spherical plain bearings with hybrid PTFE/aramid fabric liners modified by ultrasonic [J]. Tribology International, 2015(87): 132-138.

[2] BARBERO E J, TROVILLION J, MAYUGO J A, et al. Finite element modeling of plain weave fabrics from photomicrograph measurements [J]. Composite Structures, 2006, 73(1): 41-52.

[3] 张赋, 李旭东. 织物增强复合材料弹性常数的有限元法预测 [J]. 机械工程材料, 2014, (4)38: 82-86. ZHANG Fu, LI Xudong. Finite element prediction of elastic constant of fabric reinforced composites [J]. Materials for Mechanical Engineering, 2014, 38(4): 82-86.

[4] NEEDLEMAN A, TVERGAARD V. Comparison of crystal plasticity and isotropic hardening predictions for metal-matrix composites [J]. ASME Journal of Applied Mechanics, 1993,60(1): 70-76.

[5] SUN C T, VAIDYA R S. Prediction of composite properties from a representative volume element [J]. Composites Science and Technology, 1996, 56(2):171-179.

[6] HORI M, NEMAT-NASSER S. On two micromechanics theories for determining micro-macro relations in heterogeneous solids[J]. Mechanics of Materials, 1999, 31(10): 667-682.

[7] FANG Z, YAN C, SUN W, et al. Homogenization of heterogeneous tissue scaffold: a comparison of mechanics, asymptotic homogenization, and finite element approach [J]. Appl Bionics Biomech, 2005, 2(1): 17-29.

[8] WHITCOMB J D, CHAPMAN C D, TANG X D. Derivation of boundary conditions for micromechanics analysis of plain and satin weave composites [J]. Journal of Composite Materials, 2000, 34(9): 724-747.

[9] XIA Z H, ZHOU C W, YONG Q L, et al. On selection of repeated unit cell model and application of unified periodic boundary conditions in micromechanical analysis of composites [J]. International Journal of Solids and Structures, 2006, 43(2): 266-278.

[10] LI S G, WONGSTO A. Unit cells for micromechanical analyses of particle-reinforced composites [J]. Mechanics of Materials, 2004,36(7): 543-572.

[11] XU K, XU X W. Finite element analysis of mechanical properties of 3D five-directional braided compo-sites [J]. Materials Science and Engineering A, 2008(487): 499-509.

[12] CARVALHO N, PINHO S, ROBINSON P. Reducing the domain in the mechanical analysis of periodic structures, with application to woven composites [J]. Composites Science and Technology, 2011, 71(7): 969-979.

[13] LI S, ZHOU C, YU H, et al. Formulation of a unit cell of a reduced size for plain weave textile composites [J]. Computational Materials Science, 2011(50):1770-1780.

[14] 张超, 徐希武, 严雪. 纺织复合材料细观力学分析的一般性周期性边界条件及其有限元实现 [J]. 航空学报, 2013,34(7):1636-1645. ZHANG Chao, XU Xiwu, YAN Xue. General periodic boundary conditions and their application to micromech-anical finite element analysis of textile composite [J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(7): 1636-1645.

[15] BAKAR I A A. KRAMER O, BORDAS S, et al. Optimization of elastic properties and weaving patterns of woven composites [J]. Composite Structures, 2013, 100(5): 575-591.

Periodic boundary conditions for mechanical property analysis of 2-D woven fabric composite

XUE Yahong1，2, CHEN Jigang1，2, YAN Shicheng1，2, LUO Junting1,3

(1.SchoolofMechanicalEngineering,YanshanUniversity,Qinhuangdao,Hebei066004,China; 2.AviationKeyLaboratoryofScienceandTechnologyonGenericTechnologyofSelf-LubricatingSphericalPlainBearing,YanshanUniversity,Qinhuangdao,Hebei066004,China; 3.KeyLaboratoryofAdvancedForging&StampingTechnologyandScience,MinistryofEducation,YanshanUniversity,Qinhuangdao,Hebei066004,China)

In order to accurately implement the numerical analysis for the mechanics performance of two-dimensional (2-D) woven fabric composite, the accurate boundary condition of unit cell model should be established. Based on the periodic boundary condition theory, a set of simple and universal periodic boundary equations was proposed for the 2-D woven fabric composite, and the solving method of elastic constants for engineering in the finite element analysis under the periodic boundary conditions was given. In order to verify the correctness of periodic boundary conditions, the nine-block-box structure including 9 unite cells of 2-D fabric composite were established. Taking the central cell as the reference cell, the deformation and stress distribution of the single cells under different boundary conditions were compared with the reference cell. The results indicate that all boundary surfaces of the fabric composite do not keep planar state, but present the concave and convex buckling deformation under uniaxial tensile load. In other words, the periodic property of the unite cell boundary faces is demonstrated. Furthermore, the engineering constants of the 2-D woven fabric composite can be obtained properly under the periodic boundary conditions.

2-D woven fabric composite; periodic boundary condition; unit cell; finite element model; elastic constant

10.13475/j.fzxb.20150800808

2015-08-05

2016-03-10

国家科技支撑计划资助项目(2014BAF08B03)

薛亚红(1989—)，女，硕士生。研究方向为自润滑材料、精密成形。陈继刚，通信作者，E-mail:24000082@qq.com。

TB 332

A