预成形体渗透率预测及剪切变形的影响

杨 波，金天国，毕凤阳，2，李建广

（1哈尔滨工业大学 机电工程学院，哈尔滨150001；2黑龙江工程学院 机电工程学院，哈尔滨150050）

复合材料液体模塑成型工艺（Liquid Composite Molding，LCM）如今正广泛应用于制造复杂结构的复合材料构件，然而以试验和经验为主要手段的模具设计、工艺参数的确定方法成本高、效率低，很大程度上推升了该工艺制件的价格，限制了LCM工艺的应用推广。因此采用数值模拟代替实验来预测模具填充过程中可能产生的缺陷，进而优化模具设计和工艺参数的方法越来越受到重视。

树脂在预成型体中的流动可以看做是不可压缩黏性流体在多孔介质中的渗流问题，通常采用达西定律描述，联立求解达西定律与连续性方程即可获得流场参数，进而预测工艺过程中可能产生的缺陷。在此过程中，渗透率是进行填充仿真过程中最重要的输入参数，对仿真结果精度有极大的影响。

渗透率的研究是为树脂充模仿真服务的，然而仅仅研究单胞在经纬纱线正交状态下的渗透率是不足的，因为织物铺覆在模具上时会在约束下产生变形，这样渗透率会发生较大变化。织物铺覆过程中可能发生的变形有剪切、矫直、起皱、拉伸和滑移等，其中剪切是织物变形的主要形式。通常忽略其他变形模式，认为织物铺覆在细观上的体现就是单胞剪切变形［1，2］。以某雷达天线屏蔽器的LCM工艺模具填充仿真为例，其仿真过程如图1所示。首先，建立模具型腔的几何模型，并进行离散；然后利用铺覆算法计算出织物在模具内受力变形导致的局部剪切角，根据渗透率与剪切角的关系获得预成形体的渗透率分布；最后根据渗透率分布求解达西定律获得填充过程参数。可以看出，研究预成形体在剪切变形后的渗透率变化有重要意义。

图1 填充仿真过程Fig.1 Illustration of filling simulation for LCM

目前主要通过实验测定渗透率，然而实验成本高、时间长、对环境敏感，尤其是剪切渗透率很难准确测定。因此，通过数值方法模拟流体流过预成形体单胞，获得流场参数后回带入达西定律预测复杂预成型体渗透率的方法正逐渐显出其优势。

Verleye等［3，4］针对无皱褶织物和平纹织物建立了细观几何模型，采用有限差分法求解Stokes方程预测单胞渗透率，并且研究了剪切的影响，最后通过实验进行了验证。Naoki Takano等［5］采用渐进均匀理论研究了织物细观和微观渗透率，采用有限元法求解流动方程，并对单胞剪切渗透率进行了深入探讨。Chen等［6］采用FLOTRAN CFD软件仿真了平纹织物单胞中的细观流动，预测了细观渗透率，并通过大量仿真研究了几何参数对渗透率的影响。国内研究者目前大都采用实验法测定渗透率［7－10］，对数值模拟方法预测渗透率的研究很少。其中，戴福洪等［11］采用均匀化法预测了单胞渗透率，并通过与文献中的实验结果比较，验证了其方法的合理性。倪爱清等［12］采用均匀化理论分析了流体在多孔介质中的流动问题，推导出广义达西定律，分别对单向纤维织物和二维平面织物的渗透率进行了验证，考察了单胞微观结构对渗透率的影响。吴炎等［13］采用有限元法针对双向缝合毡单胞进行了流动分析，预测了渗透率，并通过一组实验进行了验证。

目前的研究主要有以下局限：（1）计算大都将纱线视为不可渗透的，忽略了纱线内部的树脂流动，这样的预测结果不能体现预成形体单胞的双尺度特性；（2）受几何建模或计算方法的限制，预成形体单胞的几何模型被大大简化：纱线截面简化为矩形，纱线的卷曲路径也采用分段直线来代替。这样的简化损害了渗透率预测的精度；（3）研究大都针对正交单胞，对剪切变形单胞的渗透率求解方法研究较少，渗透率预测结果对实际应用的帮助有限。本工作考虑纱线的内部流动建立树脂流动模型，建立逼近真实的预成形体正交单胞几何模型；基于Chorin投影法和Adams－Bashforth差分格式构造求解树脂流动控制方程的高分辨率TVD数值格式，预报预成形体渗透率；在此基础上，针对剪切单胞，采用贴体坐标法实现其渗透率的预测，研究单胞剪切变形对其渗透率的影响。

1 流动数学建模

根据Darcy定律，预成形体渗透率可以用下式表示：

其中，K表示多孔介质的渗透率张量；u＝（u，v，w）表示树脂流动速度，上划线表示体积平均；μ表示树脂的黏度；p表示压力。流场速度和压力分布为未知量，需要通过数值求解预成形体中流场控制方程获得。

宏观预成形体由纱线编织成纤维布，再通过铺覆过程获得，结构复杂度非常高，不可能一次性仿真构件预成形体内的流动，但是根据其周期性可知，预成形体是由大量重复排列的编织结构成，称之为细观单胞，以下简称单胞，图2所示为平纹织物单胞。考虑纱线的可渗透性，单胞内部的树脂流动由两个尺度的流动组成，一是纱线间区域的树脂流动，称之为细观流动；二是树脂在纱线内部纤维单丝间的流动，称之为微观流动。

图2 平纹织物单胞内的树脂流动Fig.2 The resin flow of plain fabric unit－cell

将树脂看做不可压缩黏性牛顿流体，填充过程视为层流、恒温的流动过程，这样，纱线间的树脂细观流动可以用不可压缩Navier－Stokes方程表示：

其中，f表示液体的单位质量力，ρ表示流体密度，Δ表示Laplacian算子。

由于树脂流动速度很慢，忽略非线性对流项，另外，树脂的密度对渗透率的预测没有影响，因此树脂在纱线间的细观流动可以用下式表达：

对于纱线内部纤维间的流动，可以将纱线看做多孔介质，采用Brinkman方程描述流动：

其中，Kyarn表示纱线内部的微观渗透率张量，用来表达纱线内部纤维间孔隙的渗透性能，受到纤维分布状态的影响，Kyarn通常通过经验公式（如Gebart模型［14］）获得。对比式（3）和式（4）可以发现，Brinkman方程仅细观流动控制方程多一个余项μu／Kyarn，因此，可以将单胞内的整个流动全部采用Brinkman方程来描述，纱线外部的细观流动可看做是K－1yarn＝0的Brinkman流动。

联立求解式（4）和连续性方程，即可获得单胞内流场的压力分布和速度分布，回带入式（1），即可求得织物单胞的渗透率值K。

2 正交单胞渗透率预测

2.1 几何建模及边界条件

预成形体单胞几何模型通常采用通用CAD（Computer Aided Design，计算机辅助设计）软件（如Pro／E、CATIA等）建立，受其限制，很难建立逼真的单胞几何模型，因此需要将纱线截面近似为多边形，纱线中心线路经简化为分段直线，这样的简化模型与单胞的真实结构差距很大，虽然方便了建模和计算，但是损害了渗透率预报精度。本工作利用ACIS几何建模内核，采用曲线曲面建模方法建立该逼真的织物正交单胞几何模型如图3所示。

图3 平纹织物正交单胞几何模型Fig.3 Geometry model of plain fabric orthogonal unit－cell

预测单胞细观X向（Y向类似）渗透率时的求解区域及边界条件设置如图4所示，在入口Γ3和出口Γ6施加恒压边界条件分别为P0和0；由于预成形体是由多层织物堆叠而成，而且单胞具有周期性，在其余边界施加周期性边界条件，即Γ1和Γ4，Γ2和Γ5分别成对施加周期边界条件。

图4 求解区域及边界条件设置Fig.4 Calculation region and boundary conditions

基于有限差分法求解单胞中树脂流动控制方程。对流动区域几何模型采用矩形网格离散，为了避免计算产生没有物理意义的锯齿形流动速度或压力分布，本工作采用交错网格布置方式。以二维网格布置为例，交错网格布置方式如图5所示。压力在由实心点表示的网格点上计算，水平和垂直流动速度分别在由灰色实心点和空心点表示的网格点上计算。这样压力和速度就在不同的网格点上计算，可以消除产生锯齿形流速和压力分布的可能性。三维交错网格与二维相似，不再赘述。

图5 交错网格布置Fig.5 The staggered grid

2.2 细观流动数值求解

由于本工作的流动计算考虑了纱线的可渗透性，而且纱线的渗透率较低，纱线内部的流动阻力会在纱线和孔隙界面处造成明显的速度不连续现象，这样在采用普通差分格式计算流动时，可能会呈现出震荡现象，影响计算的稳定性，这种震荡完全是因为数值的原因。为了避免上述现象的发生，应该选择一种满足总变差减小（TVD）条件的高分辨率时间离散格式，本工作采用满足该条件的Adams－Bashforth离散格式来进行时间离散。

Adams－Bashforth格式用于对流方程时是一个一步，三层和时间前差的差分格式，它的截断误差为O（Δt2，Δx2），它可以被看成是对二阶时间导数有限差分近似。Adams－Bashforth的基本形式为：

其中，Γ表示所计算的物理量，n表示时间步，Δt表示时间步长，i，j，k表示节点序号。

采用Chorin投影法解耦单胞内流动的动量守恒方程（4）和连续性方程，它是一种预估校正法。第一步利用动量守恒方程求解获得预测速度场，这个预测速度场不满足连续性方程。第二步通过连续性方程推导出用于求解压力场的泊松类方程，通过求解该方程获得下一时间步的压力分布。最后利用第二步求得的压力场将预测速度场映射到无散度的速度场中，该速度场即为下一时间步的速度场分布。将Brinkman方程的时间项利用上述Adams－Bashforth差分格式离散，并应用到该映射法中，则整个计算过程可以表示为：

（1）已知第n时间步的速度分布un，求解下式获得预测速度场u＊：

（2）联立：

与连续性方程

获得泊松类方程：

求解式（9）获得n＋1时间步的压力分布pn＋1。

（3）将第一步求得的u＊和第二步求得的pn＋1带入式（7），获得n＋1时间步的速度分布。

上述公式中Δt表示时间步长，对于显式计算的时间推进法，时间步长有如下稳定性条件：

其中，Δd表示最小空间网格步长（包含三个方向）。

对于压力项和扩散项的离散，均采用具有二阶精度的中心差分格式：

在以上算法中，每一步迭代均需要且仅需要使用前两个时间步的流场信息，因此启动该迭代需要提供起始两个时间步的流场信息，这两个时间步的流场信息可以通过很多方法获得，本工作前两步采用隐式算法获得。上述显式差分计算算法流程如图6所示。

图6 算法流程Fig.6 The algorithm flowchart

2.3 细观流动求解结果分析

为了验证以上算法的正确性，本节预测Syncoglass R420织物（表1）的渗透率，其正交单胞几何建模结果如图7所示。

表1 Syncoglass R420织物几何参数Table1 Geometry parameters of Syncoglass R420fabric

图7 Syncoglass R420织物单胞模型Fig.7 Meso－geometry model of Syncoglass R420fabric

根据经验公式求得Syncoglass R420织物纱线轴向和径向渗透率分别为kalong＝2.71×10－12m2，ktrans＝5.88×10－13m2。采用本文所述算法计算该织物单胞内的流场参数，边界条件为入口压力P0＝1Pa，流体黏度μ＝1Pa·s。X方向流动仿真结果的压力分布和速度分布分别如图8，9所示。

图8 X向流动压力分布Fig.8 Pressure distribution of X－direction flow

图9 X向流动速度分布Fig.9 Velocity distribution of X－direction flow

将计算所得压力和速度分布回带入式（1）获得单胞渗透率，本工作算法预报的渗透率为Kxx＝1.657×10－10m2，Kyy＝1.873×10－10m2。戴福洪等［11］针对同一织物，采用均匀化方法，不考虑纱线的可渗透性而预报的结果为Kxx＝1.29×10－10m2，Kyy＝1.228×10－10m2。而实验测定的该织物渗透率为Kxx＝1.43×10－10m2，Kyy＝1.79×10－10m2。比较这几个结果可以看出，在忽略纱线可渗透性情况下预测的渗透率值小于实验值，而本工作的渗透率预测值略大于实验值。考虑到实验过程中模具对预成形体的挤压作用，本工作的预测结果较为准确，证明了该算法的可行性。

3 剪切变形对渗透率的影响

3.1 剪切单胞渗透率预测

预成型体在发生剪切变形时，不可避免地引起厚度以及纱线内部纤维排布和状态的变化，但是由于这些变形量很小，而且相对于剪切角度，这些因素对单胞内流动的影响非常小，因此本工作忽略这些因素，仅考虑剪切角对细观渗透率的影响。单胞剪切锁合角度（最大剪切角度）一般为15～35°，很少超过40°，与织物拓扑和材料相关。对正交单胞模型进行简单处理即可获得剪切单胞模型，剪切角为θ的平纹织物单胞如图10所示。

图10 剪切角为θ的平纹织物单胞模型Fig.10 Plain fabric unit－cell with shear angleθ

对剪切单胞树脂流动区域做贴体坐标系如图11所示（仅剪切平面），将有剪切的物理域流动求解问题转换到无剪切的计算域进行求解。（x，y）和（ε，η）分别为物理域和计算域的自变量。

对于三维剪切单胞的贴体坐标的物理域和计算域自变量有如下关系：

其中，z和ζ分别为单胞厚度方向在物理域和计算域的自变量；θ为剪切角。

可以看出，求解区域从物理域到计算域的转换可以大大简化数值计算的难度，但是树脂流动控制方程式（4）仅适用于物理域，需要将其向计算域转化，可以说求解区域的简化是以控制方程的复杂化为代价的。下面给出转换到计算域中控制方程的具体形式。

图11 贴体坐标系转换示意图（a）物理域；（b）计算域Fig.11 Coordinate transformation of body－fitted coordinates（a）physical domain；（b）computational domain

对于Brinkman方程式（4），在计算域上的表达式为：其中，S＝K－1yarn。连续性方程式（8）在计算域上的具体形式为：

即：

由于单胞剪切造成的纱线微观渗透率主方向的与坐标系方向产生夹角，在计算之前需要进行渗透率方向变换。

设渗透率主方向与坐标系方向夹角角度为θ，为了进行数值计算，需将该处纱线渗透率Kyarn转换到坐标系方向K′上，如图12所示。根据张量理论，有如下转换关系：

式中，kalong和ktrans分别为纱线的微观轴向和径向渗透率；且有k21＝k12。

图12 渗透率方向转换Fig.12 The transform of permeability direction

将式（15）和（16）作为控制方程，利用上一节所述显式计算方法进行剪切单胞下的流动参数的求解，即可获得流场参数，进而预测渗透率。

3.2 渗透率与剪切角的关系

单胞剪切角为20°时，以Y向为压力差方向的流场仿真结果如图13，14，15所示，分别表示压力、Y向速度和X向速度分布。

图13 剪切单胞压力分布Fig.13 Pressure distribution of sheared unit－cell

图14 剪切单胞Y向流动速度分布Fig.14 Y－direction velocity distribution of sheared unit－cell

图15 剪切单胞X向流动速度分布Fig.15 X－direction velocity distribution of sheared unit－cell

从上述数值计算结果可以看出，单胞发生剪切变形后，Y方向的压力差会产生较大的X 方向的树脂流动，因此坐标轴X，Y向不再是单胞渗透率的主方向。根据仿真结果可以分别计算Kxx，Kyy，Kxy，进而通过张量的标准型变换可获得剪切单胞主渗透率值K11和K22。剪切角每隔5°进行流动仿真，获得剪切单胞主渗透率比值与剪切角的关系如图16所示。通过与文献［3］，［8］比较发现，该关系曲线与实验和理论计算所得主渗透率值－剪切角关系相吻合，证明了以上剪切单胞渗透率预测算法的正确性。

图16 单胞主渗透率比值与剪切角的关系Fig.16 Relationship between the ratio of principal permeability and shear angle

4 结论

（1）考虑纱线的可渗透性，利用Brinkman方程建立了表达织物单胞内部树脂流动的控制方程，达到了将纱线内部和外部流动统一表达的目的。

（2）建立逼近真实的正交单胞几何模型，基于Adams－Bashforth差分格式和Chorin投影法构造了求解树脂流动控制方程的高分辨率TVD数值格式，利用达西定律预测了单胞的渗透率，算例表明该算法预测值与实验值更接近，验证了算法的准确性。

（3）采用贴体坐标法完成了剪切变形后流动控制方程从物理域到计算域的转换，进而实现了剪切单胞渗透率的数值预测，考察了剪切单胞主渗透率比与剪切角的关系，比较表明与文献中的实验和理论计算结果相吻合，证明了该剪切单胞渗透率预测算法的合理性。

［1］LOUIS M，HUBER U.Investigation of shearing effects on the permeability of woven fabrics and implementation into LCM simulation［J］.Composites Science and Technology，2003，63：2081－2088.

［2］HAMILA N，BOISSE P.Simulations of textile composite reinforcement draping using a new semi－discrete three node finite element［J］.Composites：Part B，2008，39：999－1010.

［3］VERLEYE B，CROCE R，GRIEBEL M.Permeability of textile reinforcements：simulation，influence of shear and validation［J］.Composites Science and Technology，2008，68：2804－2810.

［4］VERLEYE B，LOMOV S V，LONG A.Permeability prediction for the meso－macro coupling in the simulation of the impregnation stage of resin transfer moulding［J］.Composites：Part A，2010，41：29－35.

［5］TAKANO N，ZAKO M，OKAZAKI T，et al.Microstructurebased evaluation of the influence of woven architecture on permeability by asymptotic homogenization theory［J］.Composites Science and Technology，2002，62：1347－1356.

［6］CHEN Z R，YE L，LU M.Permeability predictions for woven fabric preforms［J］.Journal of Composite Materials，2010，44（13）：1569－1586.

［7］李海晨，张明福，王彪.RTM工艺增强纤维渗透率测量方法研究［J］.航空材料学报，2001，21（1）：51－54.LI H C，ZHANG M F，WANG B.Method on measuring fibre permeabilities in resin transfer molding［J］.Journal of Aeronautical Materials，2001，21（1）：51－54.

［8］田正刚，祝颖丹，张垣，等.剪切效应对纤维增强材料渗透率的影响［J］.武汉理工大学学报，2005，27（2）：4－12.TIAN Z G，ZHU Y D，ZHANG Y，et al.Influence of shearing effects on the permeability of fiber reinforcement［J］.Journal of Wuhan University of Technology，2005，27（2）：4－12.

［9］祝君军，段跃新，陈吉平，等.碳纤维经编织物定型参数及渗透特性［J］.复合材料学报，2012，29（3）：42－48.ZHU J J，DUAN Y X，CHEN J P，et al.Packifier parameters and permeability characteristics of non－crimp stitched carbon fabric［J］.Acta Materiae Compositae Sinica，2012，29（3）：42－48.

［10］罗云烽，段跃新，王伟，等.VIP工艺中L型构件渗透规律的实验研究［J］.材料工程，2009，（增刊2）：15－19.LUO Y F，DUAN Y X，WANG W，et al.Experimental study on permeability of L－type preform for vacuum infusion process［J］.Journal of Materials Engineering，2009，（Suppl 2）：15－19.

［11］戴福洪，张博明，杜善义.用均匀化方法预报平纹织物的渗透率［J］.复合材料学报，2009，26（2）：90－93.DAI F H，ZHANG B M，DU S Y.Permeability prediction of fabric preform using homogenization method［J］.Acta Materiae Compositae Sinica，2009，26（2）：90－93.

［12］倪爱清，王继辉，朱以文.复合材料液体模塑成型工艺中预成型体渗透率张量的数值预测［J］.复合材料学报，2007，24（6）：50－56.NI A Q，WANG J H，ZHU Y W.Numerical prediction of saturated permeability tensor of a woven fabric for use in the fluid simulation of liquid composite molding［J］.Acta Materiae Compositae Sinica，2007，24（6）：50－56.

［13］吴炎，晏石林，谭华.具有空隙尺寸双尺度特性的纤维织物渗透率的预测［J］.固体力学学报，2008，29（S）：80－84.WU Y，YAN S L，TAN H.Predication of permeability in a glass－fiber mat with dual－scale pore－size［J］.Chinese Journal of Solid Mechanics，2008，29（S）：80－84.

［14］GEBART B A.Permeability of unidirectional reinforcements for RTM［J］.Journal of Composite Materials，1992，26（8），1100－1133.