e7编织复合材料宏细观力学性能的数值模拟 118页

南京航空航天大学博士学位论文基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟姓名：董纪伟申请学位级别：博士专业：机械设计及理论指导教师：孙良新20070301 南京航空航天大学博士学位论文摘要本文首先从三维编织工艺——四步法入手，分析纱线的空间位置，用曲线模拟纱线的路径以形成较精确的三维四向编织复合材料的宏观预制件模型，然后采用合理的单胞划分的方法，对宏观模型进行剖分得到内胞、面胞和角胞三种细观单胞模型。通过对细观几何模型的分析，得到了三维四向编织复合材料细观结构的分布规律，还可得到常用细观几何参数和工艺参数的计算公式。在对三维编织复合材料刚度性能的分析过程中，以基于小参数渐近展开和多尺度摄动方法的均匀化理论为理论基础，用虚位移原理导出了等效弹性模量的均匀化列式及其有限元求解方程，从而通过细观和宏观有限元的数值分析方法预测了三维四向编织复合材料的等效弹性模量，预测结果与实验结果较为吻合。同时，还深入讨论分析了编织角、纤维体积分数等工艺参数对等效弹性模量的影响，得到了一些十分有益的规律和结论。然后，通过均匀化理论和有限元分析相结合的方法，对三维编织复合材料在单向拉伸和三点弯曲情况下的细观应力进行了数值模拟，应力模拟过程中首先详细讨论了基于均匀化理论的三种单胞细观周期性边界条件的施加方法。采用截面法绘制的细观应力分布图可以直观、清晰地反映单胞区域任一截面上的应力波动情况。应力模拟结果充分体现了编织角对三维编织复合材料细观应力分布状态的影响，通过细观应力分析所得的结论与实验分析结论十分吻合。为了从微观的角度研究三维编织复合材料的拉伸强度问题，本文对其进行了多个拉伸载荷步下的非线性损伤分析，用合理的强度失效判据判断三种单胞中各单元损伤与否和损伤类型，并采用刚度折减法对损伤单元进行处理。损伤分析结果表明了不同编织角复合材料微观损伤模式的不同，以此为依据本文确定了表征编织复合材料破坏机制转变的一个重要参量——临界编织角，对小编织角复合材料和大编织角复合材料分别建立了不同的拉伸强度失效准则，并以这些准则对三维编织复合材料的拉伸强度进行了数值预报，强度预报结果与实验结果较为吻合。此外，论文还编制出一套较完整的关于三维四向编织复合材料力学性能的分析软件，该软件提供了交互式的分析界面，使得人们可以更加方便、有效地对这种新型材料进行力学性能研究及设计开发。综上所述，本文采用目前较为流行的适用于研究具有周期性结构复合材料的均匀化理论，较好地解决了从细观和微观的角度研究三维编织复合材料应力和强度的问题，所取得的研究成果中细观应力的模拟、微观损伤模式的分析及强度准则的建立都是前人的研究工作中未曾涉及的，从而有效地填补了三维编织复合材I 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟料强度分析领域的空白。关键词：三维编织，复合材料，单胞，均匀化理论，等效弹性模量，细观应力，微观损伤模式，强度准则II 南京航空航天大学博士学位论文ABSTRACTInthisthesis,thespatiallocationofyarnswasanalyzedfirstaccordingtothe3Dbraidedtechnique—four-stepprocedureandthenpropercurveswereusedtosimulatetheroutesofyarnssothatmacrostructureofthefabricwaspreciselyformed.Byusingapropermethodtodividethefabricintounitcells,thecompletedmacrostructurewasdividedandthreekindsofmicrounitcell—bodycell、facecellandcornercellweresimulatedprecisely.Byanalyzingthemicrostructure,someimportantlawsaboutthemicrostructureof3D4directionalbraidedcompositeswerefoundandsomegeneralgeometricandtechnicalparameterswereeasilycalculated.Inordertoanalyzethestiffnessof3D4directionalbraidedcomposites,thehomogenizationtheorybasedonthemulti-scaleperturbationtheorywasusedasthebasictheory.Thevirtualdisplacementprinciplewasusedtoeducethehomogenizationformulationofeffectiveelasticmodulusandtheschemestosolvetheproblemusingfiniteelementmethod.Theeffectiveelasticmodulusof3Dbraidedcompositeswerepredictedbythismethodandthecalculationresultsagreequitewellwithexperimentresults.Atthesametime,theinfluencesofsometechnicalparameterssuchasbraidangleandfibervolumefractiononeffectivemoduluswerediscussedindetail,andlotsofvaluableconclusionsareobtained.Thenmicro-stressof3Dbraidedcompositesundertheloadingofuniaxialtensionandthreepointbendingwassimulatedbythemethodofthehomogenizationtheorycombinedwithfiniteelementanalysis.Atthebeginningofmicro-stresssimulation,theproblemofperiodicboundaryconditionswhichwereappliedforthreekindsofunitcellswasdiscussedindetailfirst,thenthesectionanalysismethodwasusedtosimulatethefluctuateofmicro-stressonarandomsectionwithinacertainunitcell.Somesimulationresultsshowthatbraidanglehasagreatinfluenceonthedistributionstatusofmicro-stress.Conclusionsobtainedbymicro-stressanalysisagreequitewellwithexperimentconclusions.Inordertostudytensilestrengthof3Dbraidedcompositesfrommicrocosmicangle,nonlineardamageundermanytensileloadstepswasanalyzed,propercriteriaofstrengthwereusedtojudgedamageconditionanddamagepatternofeachelementinunitcells,andthenstiffnessdegradationmethodwasusedtodealwithdamagedelements.Theresultsofdamageanalysisindicatethatcompositeswithdifferentbraidanglehavedifferentmicrofailurepattern.Accordingtothisconclusion,animportantIII 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟parameterwhichreflectthetransformoffailuremechanismof3Dbraidedcomposites—criticalbraidanglewasdefinedandthendifferenttensilecriterionwasestablishedforsmallanglebraidedcompositesandlargeanglebraidedcomposites.Thesecriteriawereusedtopredictthetensilestrengthof3Dbraidedcomposites.Thecalculatedresultsofthestrengthagreequitewellwiththeexperimentalresults.Moreover,inthisthesis,arelativelycompletedsoftwareforanalyzingthemechanicalpropertiesof3D4directionalbraidedcompositeswasdesigned.Thissoftwareprovidesauserinteractiveinterfacesothattheresearch、designanddevelopmentonthekindofnewmaterialcanbeeasilyrealized.Inaword,inthisthesis,apopularmethod—homogenizationtheorywhichwasoftenappliedincompositeswithperiodicstructurewaseffectivelyusedtosolvetheproblemofstressandstrengthfrommicrocosmicangleof3Dbraidedcomposites.Theresearchachievementsofmicro-stresssimulation、microfailurepatternanalysisandmicrostrengthcriteriaestablishmentareallnotmentionedinpreviousresearchworkofotherssothattheblankfieldinmicrostrengthanalysisof3Dbraidedcompositeshasbeeneffectivelymadeup.Keywords:3Dbraided,composite,unitcell,homogenizationtheory,effectivemodulus,micro-stress,failurepattern,strengthcriterionIV 南京航空航天大学博士学位论文图表清单图1.1三维编织复合材料应用的产品示意图.............................................................................2图1.2三维编织机示意图.............................................................................................................3图1.3一个机器循环中携纱器的移动规律.................................................................................4[12]图1.4RTM成型工艺装置示意图...........................................................................................4图2.13×3编织物的纤维束在横截面的投影............................................................................13图2.2编织物的三维视图...........................................................................................................14图2.33×3编织物的纱线运动的动画仿真................................................................................15图2.4编织物中几种单胞分布示意图.......................................................................................15图2.5单胞划分方案一...............................................................................................................16图2.6单胞划分方案二...............................................................................................................17图2.7三维编织结构的单胞模型...............................................................................................17图2.8纱线在内胞和角胞处的过渡...........................................................................................18图2.9纱线在三种单胞处的过渡...............................................................................................18图2.10三种编织角的关系...........................................................................................................19图2.11纤维体积分数Vf随m、n的变化情况（γ＝30°）..................................................21图2.12纤维体积分数Vf随内部编织角γ的变化曲线（n＝4）...........................................21图3.1复合材料的周期性结构及单胞示意图...........................................................................25图3.2面胞和角胞中纤维束路径的折线段假设.......................................................................31图3.3三种单胞的总体坐标系...................................................................................................32图3.4内胞的局部坐标系...........................................................................................................33图3.5三种单胞的有限元模型示意图.......................................................................................34图3.6均匀化理论同刚度体积平均法Ez预测值的比较.........................................................37图3.7弹性模量和泊松比随编织角的变化曲线.......................................................................39图3.8弹性模量和泊松比随纤维体积分数的变化曲线...........................................................40图4.1内胞三个方向的完全周期性...........................................................................................45图4.2面胞两个方向的周期性...................................................................................................45图4.3角胞一个方向的周期性...................................................................................................46图4.4三种单胞的有限元模型...................................................................................................48kl图4.5内胞的等效位移χ(y)及变形图.................................................................................49kl图4.6面胞的等效位移χ(y)及变形图.................................................................................49kl图4.7角胞的等效位移χ(y)及变形图.................................................................................50图4.819×5编织物中单胞分布示意图......................................................................................51图4.9小编织角复合材料的宏观有限元模型...........................................................................52图4.10小编织角复合材料的宏观拉应力σ分布.....................................................................52z图4.11小编织角复合材料的内胞纤维应力分布（γ=21°）....................................................54图4.12小编织角复合材料的面胞纤维应力分布（γ=21°）...................................................54图4.13小编织角复合材料的角胞纤维应力分布（γ=21°）...................................................54图4.14三种单胞在z=h/2处截面示意图（γ=21°）................................................................55图4.15三种单胞在z＝h/2截面处σ分布图（γ=21°）.........................................................55z图4.16三种单胞在z＝h/2截面处mises等效应力σ分布图（γ=21°）............................55vmVII 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟图4.17大编织角复合材料的宏观有限元模型（γ=48°）.......................................................56图4.18大编织角复合材料的宏观拉应力σ分布（γ=48°）.................................................57z图4.19三种单胞在z＝h/2截面处σ分布图（γ=48°）.........................................................57z图4.20三种单胞在z＝h/2截面处mises等效应力σ分布图（γ=48°）............................57vm图4.21纤维束正应力σ随内部编织角的变化情况.................................................................581图4.22纤维束σ、σ、τ随内部编织角的变化情况.........................................................592313图4.23弯曲试件的形状及尺寸...................................................................................................60图4.24小编织角复合材料的宏观弯曲正应力σ的分布图（γ=21°）.................................60z图4.25弯曲细观应力的分析单胞位置.......................................................................................61图4.26单胞内弯曲正应力σ的截面分布图（γ=21°）.........................................................62z图4.27大编织角复合材料的宏观弯曲正应力σ的分布图（γ=48°）.................................62z图4.28单胞内弯曲正应力σ的截面分布图（γ=48°）.........................................................63z图4.29弯曲挠度随编织角变化曲线...........................................................................................64图4.30纤维束最大弯曲正应力及最大mises等效应力随编织角变化曲线............................64图4.31基体最大弯曲mises等效应力随编织角变化曲线........................................................65图5.1最大主应力理论所对应的等效应力分布图...................................................................73图5.2最大主应变理论所对应的等效应力分布图...................................................................74图5.3小编织角复合材料的拉伸应力－应变曲线...................................................................74图5.4失效单元体积比随工艺参数变化情况示意图...............................................................75图5.5计算曲线与实验曲线的比较示意图...............................................................................76图5.6失效单元的体积百分比随编织角的变化情况...............................................................77图5.7编织角为48°复合材料的拉伸应力－应变曲线..........................................................78图5.8几种大编织角复合材料的拉伸应力－应变计算曲线...................................................78图5.9内胞单元损伤演变情况...................................................................................................82图5.10面胞单元损伤演变情况...................................................................................................83图5.11角胞单元损伤演变情况...................................................................................................85图6.1三维四向编织复合材料刚度预测流程图.......................................................................89图6.2三维四向编织复合材料应力模拟流程图.......................................................................90图6.3三维四向编织复合材料应强度分析流程图...................................................................91图6.4三维四向编织复合材料力学分析的工具条按钮界面...................................................92图6.5输入编织工艺参数对话框...............................................................................................93图6.6输入材料参数对话框.......................................................................................................93图6.7刚度预报输出结果示意图...............................................................................................95图6.8宏观力学分析类型的选择...............................................................................................96图6.9宏观力学分析输入载荷条件的对话框示意图...............................................................96图6.10输入截面位置对话框.......................................................................................................97图6.11强度预报输出结果示意图..............................................................................................97表3.1织物所用碳纤维和树脂基体的弹性常数………………………………………………..30表3.2三维编织复合材料试件的工艺参数.................................................................................30表3.3局部坐标系为x’y’z’与总体坐标系xyz的方向余弦....................................................31表3.4内胞中纤维束坐标系间的方向余弦.................................................................................33表3.5工程弹性常数的计算、实验结果比较.............................................................................36表3.6三种单胞所占比例及拉伸模量.........................................................................................36VIII 南京航空航天大学博士学位论文表3.7纤维和树脂的力学参数.....................................................................................................37表3.8试样的细观结构参数.........................................................................................................38表3.9工程弹性常数结果比较.....................................................................................................38表4.1织物所用碳纤维和树脂基体的弹性常数.........................................................................47表5.1损伤类型及参量.................................................................................................................71表5.2拉伸强度和断裂应变的预测值和实验值比较.................................................................79表5.3材料的主要参数.................................................................................................................80表5.4七种材料的λ值.................................................................................................................86IX 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟注释表符号含义d圆截面假设纤维束直径r圆截面假设纤维束半径h花节高度Wm预制件长度Wn预制件宽度γ内部编织角θ表面编织角β角柱编织角Df纤维束等效直径λ纱线线密度ρ纤维密度ε纱线填充因子Ω纤维束实际截面积a纤维束椭圆长轴长b纤维束椭圆短轴长Vf纤维体积分数σ应力χ等效位移Cijkl刚度矩阵H等效刚度矩阵CijklX 承诺书本人郑重声明：所呈交的学位论文，是本人在导师指导下，独立进行研究工作所取得的成果。尽我所知，除文中已经注明引用的内容外，本学位论文的研究成果不包含任何他人享有著作权的内容。对本论文所涉及的研究工作做出贡献的其他个人和集体，均已在文中以明确方式标明。本人授权南京航空航天大学可以有权保留送交论文的复印件，允许论文被查阅和借阅,可以将学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或其他复制手段保存论文。(保密的学位论文在解密后适用本承诺书)作者签名：日期： 南京航空航天大学博士学位论文第一章绪论1.1三维编织复合材料的发展背景复合材料是由两种或两种以上不同性能、不同形态的组分材料通过复合手段组合而成的一种多相材料，通过组分材料的相互补充和关联，获得组分材料所不具有的新的材料特性。由于复合材料具有重量轻，比强度、比模量高，可设计性能好等优点，它已被广泛地应用于航空、航天、建筑、机械、汽车、桥梁、化工、医药及运动器材等各个部门。传统的复合材料是通过把纤维束按一定的角度和一定的顺序进行铺层或缠绕而制成的，由于形成层状结构，因此也称层合（压）复合材料。层合复合材料已在诸多领域得到深入的研究和广泛地应用，但是由于“层”的存在也带来了诸多力学性能的弱点：如层间强度低而易分层、抗损伤冲击能力较差、厚度方向的[1-6]刚度小、冲击韧性和损伤容限水平低、抗面内剪切强度低等。20世纪80年代发展起来的纺织结构复合材料，利用纺织技术，首先以纤维束织造成所需结构的形状，形成预成型件，然后以预成型件作为增强体进行浸胶固化，直接形成复合材料结构。通过增强纤维束在空间交织成一整体，纺织结构复合材料克服了传统层合复合材料厚度方向力学性能的诸多缺陷，使其具有优良的抗损坏性和良好的断裂韧性。此外，通过选择适当的工艺和工艺参数，可以织造出非常复杂形状的结构，从而使该材料具有很强的可设计性。因此，纺织结构复合材料近年来引起了力学和材料研究者的高度重视，并且已被广泛地应用于航空、航天、汽车和军事等各个领域。按照技术上分类，纺织结构复合材料包括编[7]织（braiding）、机织（weaving）、针织（knitting）和缝合（stitching）复合材料。编织复合材料与另三类复合材料相比，更具有方向的均衡性和高抗冲击任性，而[8]受到高度重视和关注。20世纪70年代初，在缠绕工艺的影响下，二维编织工艺被引入复合材料领域，其后得到了迅速发展。80年代，通过纺织界和复合材料界的合作，编织技术由二维发展到三维，为制造高性能复合材料提供了新途径，三维编织结构复合材料由于增强体为三维整体结构，大大提高了其静态和动态性能，从而倍受材料科学界、力学界和工程界的重视并获得了飞速发展，目前已在现代工程技术的各个领域显示了独特优势，成为理想的高性能结构材料。三维编织复合材料首先在航空航天领域得到应用。美国航空航天局（NASA）于80年代末制定了先进复合材料技术（AdvancedCompositesTechnology）计划，目的是突破先进复合材料1 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟制造和应用的关键技术。ACT的重要内容之一是纺织复合材料结构的开发和设计，其中特别提出了开发三维编织技术、自动化加工和开发热塑性树脂等重要内容。美国波音（Boeing）和洛克希勒（LockheedMartin）公司研究了编织技术在[9]飞机机身部件的制造；美国Brunswick公司的国防公司用编织结构复合材料制作了2000多个导弹弹翼和航天器接头。据报道，美国先进技术轰炸机采用了非圆形的特殊碳纤维与玻璃纤维混杂编织成三向织物作为复合材料的增强材料，由于这种编织物增强的复合材料像钢一样坚硬，同时又具有优良的吸透波的性能，用该材料制造B-2飞机的多层吸透波蒙皮，既能满足飞机气动特性的要求，又能减少雷达波的反射截面。三维编织复合材料还用来制造飞机的加筋蒙皮、机身框、[10]舱窗间壁板、机身龙骨框等结构件，使得结构更加耐久、损伤容限更高、重量更轻和成本更低。此外，编织结构复合材料还在机械、建筑、军事、体育、医疗、交通等高科技领域得到应用。总之，三维编织复合材料作为一种新型材料比很多传统复合材料所具有更优越的力学性能，使得它的广泛应用已取得了可喜的成绩。（a）加筋蒙皮（b）开窗板条（c）周围框（d）龙骨梁架（e）锥体织物(f)高弹性壁球拍图1.1三维编织复合材料应用的产品示意图1.2三维编织工艺简介三维编织复合材料的加工工艺包括编织工艺和成型工艺两个步骤。所谓三维编织工艺，就是由增强体纤维束（又称纱线）形成编织预制件（又称织物）的过程，它是在三维编织机上完成的。目前比较流行的三维编织工艺主要有二步法和四步法两种。在二步法编织中，轴向纱不动，编织纱围绕着轴向纱而运动；四步法编织过程中，所有的纱线改变它们的位置而互相缠绕在一起，为了增强纵向的2 南京航空航天大学博士学位论文刚度和强度，经常在结构中放置一些轴向纱，其中不含轴向纱的为四向编织复合材料，含有轴向纱的有五向、六向等多向编织复合材料。根据织机和产品形状，每种方法又分为矩形编织和圆形编织。矩形编织用来生产具有直角结构的织物，如方块、T形梁、工字梁等。圆形编织则生产具有圆形横截面的旋转体，如圆棒、[11]圆管、圆锥管等。图1.1（e）所示的锥体织物就属于圆形编织物。本论文主要研究编织工艺最为简单的四步法三维四向矩形编织复合材料。图1.2为四步法矩形编织的三维编织机示意图。图1.2三维编织机示意图四步法编织时，先把携纱器移到各自初始位置，然后通过多次移动携纱器和对已交织纤维束的打紧操作织出预制件。一个编织循环中每个携纱器都移动四步，故称四步法。根据行和列之间移动位置的相对关系，四步法编织工艺可进一步分为1×1、1×3、1×5和1×1×1/2F等，其中第一位数字代表携纱器每次沿横（列）向移动的步距，第二位数字则代表携纱器每次沿列（横）向移动的步距，1×1×1/2F编织模式除了按1×1方式编织外，在结构内部还存在着编织纱线总数目一半的轴向纱线。不同的编织模式将得到不同类型的织物实体。由于1×1方式最为简单而且应用最为广泛，因此，通常被用作四步法编织工艺的代表来讨论。对于1×1方式编织，携纱器的运动规律如下：第1步：奇数行左移一列，偶数行右移一列。第2步，奇数列下移一行，偶数列上移一行。第3步，奇数行右移一列，偶数行左移一列。第4步，奇数列上移一行，偶数列下移一行。3 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟图1.3为1×1方式四步法矩形编织的一个机器循环中携纱器移动示意图。图1.3一个机器循环中携纱器的移动规律上图中的虚线框内的矩形区域称为中心区，又称主体区。当引纱器作行移动时，中心区外的两行引纱器不动；当列移动时，中心区外的两列引纱器不动。若中心区阵列的行数为m，列数为n，则称阵列为m×n阵列，而称相应的织物、材料和编织为m×n织物、m×n材料和m×n编织。图1.2则为3×3编织情况。1.3成型工艺简介三维编织的预制件需与基体材料复合成型才能得到性能优良的复合材料，因此成型工艺对三维编织复合材料的性能有着显著的影响。适合于编织结构复合材料的成型工艺主要有真空浸渍法、编织－拉挤法、热塑性数值基体编织复合材料成型法、预浸纱成型技术和RTM工艺等，其中RTM工艺的应用最为广泛，被业界认为最理想最有效的成型方法。RTM（ResinTransferMolding）工艺是将预成型的增强体置于密闭的模具中，在真空和压力的条件下，树脂被注入模腔而固化成型的过程。与其他成型工艺相比，RTM成型工艺的优点在于：首先，该工艺最适合生产尺寸比较大、外形与结构比较复杂的产品；其次可以减少环境污染，提高环保水平；最后它易于实现自动化和计算机控制，减少附加设备的费用，具有良好的综合经济效益。图1.4为RTM成型工艺装置示意图。[12]图1.4RTM成型工艺装置示意图4 南京航空航天大学博士学位论文1.4三维编织复合材料的研究进展概述三维编织复合材料作为一种新型高性能的复合材料，自产生以来，就受到国内外学者的高度重视，人们根据它的制造工艺提出了各种各样的力学分析模型，并开始对它的刚度、强度等力学性能进行预测。但是由于三维编织复合材料的空间拓扑结构非常复杂，而且影响三维编织复合材料力学性能的参数很多，使得力学分析较为复杂。近年来国内外关于三维编织复合材料力学性能的研究主要还是以刚度分析为主，涉及强度分析的文章还比较少。在国内，天津工业大学（原天津纺织工学院）率先对三维编织复合材料进行了研究，在编织技术、编织设备、计算机辅助设计、复合固化、微观结构和力学性能等方面取得了显著的成就。经过十多年的不懈努力和科技攻关，在航空航天部门的大力支持下，天津工业大学复合材料研究所已成为国内最先进的编织复合材料的研究和生产基地，能够满足航空、航天等部门对应用三维整体编织技术来制造高性能复合材料的迫切要求[13]。除了天津工业大学以外，中材科技股份有限公司、北京航空航天大学、南京航空航天大学、哈尔滨工业大学、西北工业大学、中国科技大学等科研单位也对三维编织复合材料的细观结构和力学性能进行了深入细致的实验和理论研究。1.4.1力学分析模型研究三维编织复合材料由于在细观结构上具有一定的分布规律，从细观力学的角度出发从理论上研究其力学性能并设计和优化复合材料已成为目前编织复合材料研究和开发的一个重要手段。从上世纪八十年代以来，人们就开始对材料的细观力学性能进行预测，并且根据编织和成型工艺提出了各种各样的力学分析模型。弄清几何细观结构、建立合理的力学模型已成为目前的研究热点。早期的细观力学模型主要是基于层合板理论的分析模型。在这方面，Yang，[14]Ma和Chou首先提出了“纤维倾角模型”，将复合材料的单胞结构处理成具有相同厚度的四块倾斜层板组成，对每一单向层板进行二维应力－应变分析；周光[15][16]明、孙慧玉等把“纤维倾角模型”发展为“层板模型”，考虑厚度方向的效应，对每一层板进行三维应力－应变分析，对三维编织复合材料弹性性能进行了[17]预测；Byun，Whitney，Du和Chou利用经典层合板和刚度平均法对二步法编[18]织复合材料的刚度特性进行了分析；近年来，陈利等对纤维倾角模型和层板模型进行了改进，考虑了编织纱线的排列方式、编织工艺参数、表面纱线对其力学性能的影响，建立了三维四向编织复合材料的叠层板分析模型。这种基于层合板理论的单胞分析模型能够较为准确地预测材料的刚度特性，但是由于其忽略了纤维束之间和层板之间的互锁和相互作用，忽略纤维在材料表面和角柱处的弯曲，与实际细观结构有较大的出入，不利于进行复合材料的细观应力和强度分析。5 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟近年来，大多数的研究对编织工艺和流程进行分析，找出三维编织复合材料细观有规律的、具有代表性的体积单元（RVE）用以反映出结构整体的宏观性能。该代表性体积单元反映了宏观整体结构的组分信息、编织工艺参数等，被称之为[19]单胞（unitcell）。1982年，Ko最早提出了反映编织成型工艺的六面体模型，该模型总共包含四根直纤维束，分别沿着六面体的四个对角线方向，因此又称为“四根直纤维模型”或“米字形模型”。这个简化模型中四根纤维束彼此相交，[20,21]与实际情况很不相符。KalidinliS.R.等根据Ko的模型提出了一个等效的模[22,23]型，该模型考虑了纱线的弯曲，被称为“四根弯曲纤维模型”。LiChen等对编织工艺中纱线的位置进行了仔细地研究，采用“十二根直纤维模型”作为内胞[8]模型，冯淼林对LiChen的内胞模型进行了修正，提出了“十二根弯曲纤维模型”作为代表性的体积单元。这些单一的单胞模型的建立，可使力学分析大大简化，但是忽略了内部、表面和角柱处纤维束路径的不同，因此与实际模型有一定的出入。经较详细研究而建立的“三胞模型”，由美国NorthCarolinaStateUniv.的[24]LiWei等人于1990年提出。DrexelUniv.的WangYouQi等人于1994~1995[25-27]年作了进一步发展。它指出：材料的内部纤维是直线路径，位于平行于材料长度方向且与材料侧面成45°的两组互相正交的平面上；存在三种单胞即内胞、面胞和角胞；纤维在三种单胞中的编织角不同，但彼此间有确定关系。他们导出了相应的几何参数尺寸和刚度的计算式。1993年，北京宇航系统工程研究所的[28、29]吴德隆和郝兆平用此模型导出刚度公式并研究材料的双模量、弹塑性本构[30][31][32]关系和损伤对材料性能的影响。黄小平等分别在在文献、文献和文献中分别分析了四步法1×1矩形编织预成形织物的内部结构、表面结构和角柱结构，对纤维束的走向曲线、位置关系等几何参数做了细致而深入的研究，建立了较准[33]确的几何模型，并在文献中提出了新的单胞划分方法，这些都为宏细观有限元[22]计算提供相应的几何参数。LiChen等通过对纤维束空间走向分析，利用周期性分布的特点，建立了三种不同的单胞模型，并分析了织物几何结构参数间的数学关系。“三胞模型”符合细观结构的实际情况，考虑到了纤维束路径和倾角在材料内部、表面和角柱处的不同，并去除了层合板理论中种种人为假定，对进行编织复合材料应力和强度分析提供可靠的数据，因此这种模型目前已被广泛地采用。1.4.2刚度性能研究刚度（又称等效弹性性能或有效性能）表征了复合材料宏观力学性能。由于关于复合材料刚度研究的方法相对比较简单，且理论已经比较成熟，人们对三维编织复合材料力学性能的研究首先从刚度开始。近年来，国内外关于三维编织复6 南京航空航天大学博士学位论文合材料的刚度性能研究的文献很多，既有实验研究也有理论研究。刚度的实验研[34-40]究主要有拉伸实验、压缩实验、剪切实验和弯曲实验等。早期的理论研究方法主要有修正基体法、层合板理论分析法、刚度平均法、[41]弹性应变能法等。修正基体法是利用平均化的观点，通过修正基体的材料性能，将三维的复合材料简化为二维，然后计算二维编织复合材料的弹性模量，这种方[42]法仅适用于正交三维连续纤维增强复合材料的情况。Kregers提出了预报空间增强复合材料宏观性能的刚度平均化方法，通过对单向复合材料有效性能的体积平均得到编织复合材料的等效弹性性能，其中刚度平均法得到工程弹性常数的上限，柔度平均法得到工程弹性常数的下限，他还研究了含孔隙编织复合材料的热传导系数，利用平均化方法导出了编织复合材料热膨胀系数的上下限公式。Ma[43]等在基于材料单胞应变能分析基础上建立了预报三维编织复合材料刚度系数的方法，该模型注重单胞内纤维束互相交织的作用失效，计算出单胞内各复合材料棒在交织力作用下的拉伸、压缩及弯曲应变能，然后根据卡式定理，导出轴向弹性模量的泊松比，它们被表示为纤维体积比和纤维束取向函数。长井谦宏、横[44、45]山敦土等以具有三维正交织物结构的工字梁为对象，建立了结构力学的刚架模型，通过与有限元法结合，求得拉伸弹性常数及拉伸强度与压缩强度，该模型将材料单胞内的纤维束和基体分别理想化成纤维束梁单元和基体梁单元，复合材料单胞被看成是这些梁单元的集合体。近年来，随着计算机技术的飞速发展，有限元技术结构力学分析中所占的比例越来越大。以代表性体积单元（RVE）－单胞为分析对象，建立单胞的有限元模型，然后进行数值运算来预测三维编织复合材料力学性能，已成为目前比较流[46、47]行的刚度预测方法。庞宝君等建立了含多相介质单元的有限元分析方法，采用规则的八结点体单元对内胞模型进行剖分，可以得到三种单元：基体单元、纤维单元和混合单元。对于混合单元，如果高斯积分点位于基体处，则采用基体的[23]刚度矩阵，反之则采用纤维的刚度矩阵。ChenL等提出了有限多相单元法，这种方法以三胞模型为分析对象，首先对每一种单胞采用六面体单元进行网格划分，根据单元的内部组成，用变分法推导出每个单元的刚度矩阵，通过有限元法计算平均应力、应变，得到单胞的等效弹性矩阵，然后把各个单胞的性能代入到[48、49]整体划分模型，用同样的方法得到整个材料的等效弹性矩阵。刘振国等采用了“米”字形单胞模型，讨论了相应的边界条件和约束条件并应用有限元方法计[50]算了该材料的纵向、横向弹性模量。蔡敢为等提出了一种接近四步法三维编织复合材料实际结构的细观力学分析模型，并根据应变能等效原理，采用有限元法分析该材料的弹性模量和泊松比。运用有限元技术对单胞模型进行刚度分析便于利用有限元模块，自编相应软件，达到比较方便的对特定结构编织材料的力学性7 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟能进行预报的目的，便于形成工程使用的集材料设计、几何模拟和力学性能预报为一体的工程软件。但是，这些方法在单元刚度矩阵的选择、单胞模型的建立、边界条件的选择上都采取了一定程度的近似，都会影响到刚度的计算精度。如何建立比较符合实际情况的细观力学模型、如何恰当地简化算法而又能保证计算精度，是目前三维编织有限元分析技术急待解决的问题。二十世纪七十年代出现的均匀化理论以基于小参数渐近展开的多尺度摄动方法为数学依据，用含宏观和细观两种参量和小参数的渐近展开式来描述各种场量，可以结合有限元技术解决具有周期性分布复杂结构的刚度问题和微观应力分[51-53]布的问题，从而有利于建立合理的强度准则。中国科技大学冯淼林等通过均匀化理论对三维编织复合材料和三维压电编织复合材料的等效弹性模量进行了预测，模量计算结果具有一定的精度。与其它方法相比较，均匀化理论具有严格的数学依据，可以同时解决材料的刚度和细观应力两大问题，同时在施加单胞的边界条件方面具有明显的优越性，通过编制程序可以解决任意参数、任意工作载荷下的刚度和细观应力分布问题。因此，它是一种对具有周期性结构材料研究进行刚度和强度分析的比较理想的方法。1.4.3强度性能研究早期的三维四向编织材料强度方面的研究工作基本上可分两类：一类是做强度试验，观察损伤失效模式，分析原因机理，了解参数对性能的影响，积累数据[54-61][62-66]；另一类是兼有预估和试验，即假定或推导出一些较简单的预估公式，[67]或从层板公式导得编织材料的预估公式，将计算值和试验值比较。施加载荷种类有：剪切，弯曲，纵向拉伸、压缩，横向压缩，纵横双向压缩等。近年来，国内外的一些学者采用有限元的方法对单胞模型的应力分布情况进[68]行了数值模拟。TaoZeng等采用混合单元法对材料的细观力学性能进行研究，这种方法考虑了细观力学参数的不均匀性，形象地模拟了编织物横截面上的应力分布情况，对于材料的失效分析具有一定的参考价值，然而在建模过程中采用了“米”字型模型，忽略了表面和角柱区域效应的影响,对于细观分析精度具有一[69]定的影响。修英姝等提出了细观单胞元模型和细观单胞元方法，利用这种方法可以模拟出宏观应力分布。目前关于三维编织复合材料强度理论的研究还处于起步阶段，理论还很不成[45][70]熟。庞宝君等、TaoZeng等采用蔡－胡准则和Mises准则分别作为纤维和基[71]体失效判据，用刚度折减法对单胞模型进行损伤分析；卢子兴等采用有限元的方法对典型的四纤维内胞模型进行了数值分析，提出了以单元破坏的体积比作为该宏观强度失效判据，这种方法对于细观破坏模式分析及宏观强度预报都有非常重要的意义。8 南京航空航天大学博士学位论文1.4.4存在的不足1.4.3.1几何建模的不足在几何单胞模型的建立上，前人采用的模型主要存在着以下两点局限性：①大多数模型只考虑内部的代表性单胞，而没有考虑表面和角柱处的单胞，模型过于简化；②大多数文献提及的单元体结构并不是根据织物的实际构造而建立，而是采用了种种人为的假设，从而不能保证纱线的连续性和织物的完整性。1.4.3.2力学性能研究的不足关于三维编织复合材料的力学性能研究，前人已经取得了一些成果，但也存在着一些不足之处，主要体现在：①多数研究都围绕着刚度性能的预报展开，有关细观应力和强度分析的文献很少；②三维编织复合材料力学性能的理论研究主要以拉伸性能的研究为主，其它工况下的理论研究较少；③关于三维编织复合材料刚度性能的预报，大多数研究都采用了刚度平均的方法，而事实上这种方法只能得到刚度值的上限；④关于三维编织复合材料细观单胞的边界条件施加，目前尚无统一的施加方法与各种工况下的实际边界条件相符合；⑤关于三维编织复合材料应力计算，许多研究都是从宏观角度，计算宏观平均应力，而宏观应力无法反映单胞内部的应力波动情况；⑥关于三维编织复合材料强度性能研究的文献很少，目前尚无统一的强度准则的可供参考。1.5本文的主要工作和创新点本文的主要研究目的包括：建立符合三维编织复合材料实际构造的三胞模型；提出一种适合三维编织复合材料单胞边界条件的施加方法；采用均匀化理论与有限元技术相结合的方法对三维编织复合材料的等效弹性性能及拉伸和弯曲作用下的细观应力进行数值模拟；分析不同编织角材料的细观失效模式，建立合适的强度准则；用APDL语言对ANSYS软件进行二次开发，建立基于均匀化理论的三维编织复合材料刚度和应力分析平台；通过拉伸实验预测三维四向和三维五向编织复合材料的刚度。具体章节安排如下：第一章，对三维编织复合材料的发展和研究概况及编织工艺和成型工艺作简9 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟要介绍，着重综述并点评了三维编织复合材料力学分析模型研究、刚度性能研究和强度性能研究的国内外研究现状；第二章，用3DMAX软件建立三维编织复合材料的宏观和细观模型。用动画的形式显示纱线的运动规律，采用对整体结构剖分的方法，建立较为完整的三种单胞模型，还给出了一些重要的工艺参数的计算公式；第三章，从基于小参数渐近展开均匀化理论出发，推导出等效弹性矩阵及细观应力的均匀化列式，利用虚位移原理推导均匀化理论的有限元方程，并用此方法预测了三维编织复合材料的等效弹性模量，最后讨论了等效弹性模量随编织角、纤维体积含量等工艺参数的变化情况；第四章，首先讨论了均匀化理论周期性边界条件的施加方法，特别针对几种不同的单胞，给出了不同形式的边界条件的施加方法，然后用基于均匀化理论的有限元方法对三维编织复合材料的细观应力进行了数值模拟，分别讨论了单向拉伸和三点弯曲作用下材料的细观应力分布情况，并着重讨论了小编织角和大编织角情况下应力分布模式的不同；第五章，在前人实验分析结果的基础上，对三维编织复合材料的微观破坏机制进行深入的理论分析，然后根据前一章细观应力的计算结果，分别建立了小编织角和大编织角复合材料的拉伸强度准则，并采用非线性微观损伤分析法确定了用以区分大、小编织角的临界编织角；第六章，以ANSYS软件为二次开发平台，用APDL语言编制了一套对三维编织复合材料进行几何建模、刚度预报、应力模拟和强度分析的软件，该软件可以方便地实现用均匀化理论分析三维编织复合材料力学性能的全过程；最后针对全文内容做出了总结，并对后续工作提出了几点展望。本文的主要创新点有：①利用细观几何结构的周期性分布，采用对三维编织整体结构剖分的方法得到三种单胞模型，通过这种方法得到的单胞更接近于实物；②在用均匀化理论预测三维编织复合材料的等效弹性模量时，比文献[51]多考虑了面胞和角胞的影响；③提出了满足均匀化理论要求的单胞边界条件的统一施加方法，用该方法施加的边界条件与实际工况条件无关，仅满足周期性，这将保证均匀化理论可以方便地解决各种工况下的细观应力分布问题；④根据均匀化理论求解细观应力的有限元法求解步骤，对三维编织复合材料在拉伸和弯曲作用下的细观应力分布进行了数值模拟。目前国内外关于三维编织复合材料细观应力分布的文章非常少，其中多数研究也是直接对一种单胞模型（如“米”字形模型、弯曲纤维模型、直纤维内胞模型）进行拉伸有限元分析。10 南京航空航天大学博士学位论文而本方法考虑了面胞、角胞对细观应力分布的影响，大大提高了细观应力的计算精度；⑤采用微观损伤分析中刚度折减的方法对三维编织复合材料的非线性力学性能进行数值分析，得到了关于织物微观失效模式的较为直观的有限元模拟结果，从而为建立合理的强度失效准则提供了可靠的依据；⑥着重讨论了编织角这一重要的工艺参数对材料力学性能的影响，根据编织角大小的不同，制定了不同的强度准则，并确定了用以区分大、小编织角复合材料的临界编织角，能够从微观角度解决材料的强度问题；⑦在ANSYS软件平台上，采用人机交互式的界面，编制开发了一套适合解决三维编织复合材料建模问题、刚度预测问题、应力模拟问题及强度分析问题的综合性分析软件，这将为三维编织复合材料的研究、设计、开发和使用提供便捷的途径。11 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟第二章三维四向矩形编织复合材料的宏细观几何建模2.1引言合理地建立三维四向矩形编织复合材料的几何模型是其力学分析的前提。所谓合理地建模，就是模型能够较为真实地反映三维编织复合材料预制件的实际形状。只有建立较真实的模型，才能保证其力学性能的数值分析精度。三维编织复合材料由增强长纤维束在空间相互交织成预成型件，空间拓扑结构非常复杂，这使得如果直接对宏观模型进行有限元分析通常是自由度数目繁多，除了占用大量内存外，计算精度也难以得到保证，有时甚至计算结果面目全非。而事实上，三维编织结构的细观单胞具有一定的周期性分布规律，研究其周期性，利用其周期性可以使得宏观结构的力学分析大为简化。早期的倾斜纤维模型、层板模型、“米”字型模型等都与三维编织结构的实际形状有较大的差距，三胞模型能够较真实地反映编织结构的周期性分布规律，[22][33]因此目前已被众多学者广泛地采用。L.Chen、黄小平等利用三种单胞在宏观预制件结构中的分布情况及纤维束在三种单胞中的位置关系，采用了合理的单胞划分的方法，分别建立了内胞、面胞和角胞模型，并计算出各种细观几何参数，对于三维编织复合材料的细观结构的精确建模及分析作出了重要的贡献，但其不足之处是其细观建模直接从细观出发，只考虑了单胞中的主要纤维束，而忽略了一些穿入、穿出的残余纤维束，从而不能保证纤维束的完整性，纤维体积百分含量的精度也会受到影响。事实上，如果首先建立三维编织结构的宏观预制件模型，然后对宏观模型通过布尔运算进行合理的剖分以形成细观单胞模型，方法既简单又能提高细观模型[72]的精度。孙义林等通过计算机图形建模技术使得编织过程形象化，通过空间点对纱线位置的精确描述来建立三维编织宏观预制件模型，本文采用了这种宏观建模方法，首先从三维编织复合材料的编织工艺—四步法入手，较为精确地获得纱线编织结构的各种几何参数，根据计算机对三维空间点的精确定位，逼真地模拟出三维编织预制件模型，用3DMAX软件还可以动画的形式直观地显示出四步法三维编织过程中纱线的具体运动情况，从而为细观建模奠定了基础；然后采用合理的单胞划分的方法，对宏观模型进行剖分，得到三种细观单胞模型，通过对建立的单胞模型的观察和分析，可以得到各种单胞中纤维束路径、单胞周期性分布等重要的规律，这些规律将对下一步对三维编织复合材料进行力学分析提供重要的参考。12 南京航空航天大学博士学位论文2.2宏观结构的计算机图形建模2.2.1基本假设[73]本文在宏观预制件模型的建模过程中采用了如下假设:（1）纤维束的横截面为圆形，直径为d，且在垂直于中心线的截面内沿径向没有变形。这将保证在建模过程中，用最大编织角产生的预成型实体内各纤维束之间的相对位置应处于极限状态，即经过了打紧工艺以后，彼此间既无相互侵入，又无大的空隙，此时的预成型件为最密实的状态；（2）纤维束有足够的柔韧性。这将保证纤维束在面胞和角胞内为平滑过渡，不产生折角，在下面的模型建立过程中，面胞和角胞内的纤维束路经采用贝塞尔（Bezier）曲线拟合，与实物较为吻合；（3）编织工艺方法足够稳定。这将保证整个编织结构沿编织方向以均匀的拓扑形式扩展。2.2.2工艺参数及预制件的基本几何尺寸根据引纱器每一步的移动位置，可以确定每根纤维束所连的引纱器的轨迹，由于纤维束在牵引线的拉紧作用下取最短路径，故纤维束的轨迹在织物横截面上[74]的投影不同于引纱器的轨迹。以3×3编织物为例，通过绘制引纱器的轨迹图，不难发现以下两个规律：（1）3×3编织物每一根纤维束所连引纱器都是经过20步又回到原来的位置；（2）根据引纱器的路线的不同可以把纤维束分为三组，每组有五根纤维束，这五根纤维束所连的引纱器移动路线相同，只是走步位置依次滞后4步。由于在纵向纤维束每经过一步都移动1/4个花节长度h，通过绘制出纤维束在织物横截面上的投影图，就可以确定纤维束经过每一步后的空间位置了。图2.1为3×3编织物在织物横截面上的投影图。图2.13×3编织物的纤维束在横截面的投影13 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟花节长度h是指一个机器循环中纤维束在纵向生成的长度，其计算公式如下：h=4d/tanγ(2.1)其中γ为内部编织角，d为纤维束的直径。并且，可以得出区域（图2.1中虚线包围的区域）为m行n列的预制件的长[74]度Wm和宽度Wn的表达式如下：W=(2m+1)d(2.2)mW=(2n+1)d(2.3)n其中d为纤维束的直径。2.2.3计算机图形模拟通过对编织工艺的分析，可以按照图2.1确定每一根纤维束在每一步后的空间位置点，利用3DMAX软件对这些点进行曲线拟合，可以得到三维编织预制件的宏观模型图，如图2.2所示。俯视图主视图左视图用户视图图2.2编织物的三维视图2.2.4动画仿真我们已经得到了3×3编织物的宏观模型图，为了使用户能够更清楚直观地了解纱线在编织过程中的走向情况，我们用3DMAX软件以动画的形式显示了编织过程中织物内纱线的具体运动情况。图2.3为四步法3×3编织物纱线运动的动画仿真中的几帧画面（取每4步为一帧）。14 南京航空航天大学博士学位论文（a）用户视图观察（b）俯视图观察图2.33×3编织物的纱线运动的动画仿真从图2.3的俯视图不难看出，每一个机器循环（即4步）后，织物在横截面的投影图完全相同，只不过同一位置的纤维束发生了变化，这同前面的每一组的五根纤维走步位置依次滞后四步的结论是完全吻合的。2.3细观结构的计算机图形建模2.3.1单胞的周期性分布三维编织复合材料的空间拓扑结构非常复杂，然而实际上它是由三种具有周期性结构的单胞所构成的，这三种单胞分别是内胞、面胞和角胞，图2.4为8×6[22]织物中几种单胞的分布情况。图2.4编织物中几种单胞分布示意图15 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟2.3.2单胞划分方法如何合理地划分单胞是细观结构建模的关键。总体上来讲，单胞划分应遵循以下两条原则：①周期性原则。即划分的单胞应尽可能的反映如图2.4所示细观结构的周期性分布规律，因为只有这样才能利用其周期性简化细观力学计算；②完整性原则。即宏观预制件结构应基本上完全由内胞、面胞和角胞三种单胞组成，应尽量减少划分完成后还存在多余的区域既不属于内胞，又不属于面胞，又不属于角胞。关于单胞的划分问题，根据内胞尺寸的不同，目前主要有两种划分方法，其中以内胞尺寸为2d×2d×h最为普遍。按图2.4中方式的划分方法即属于这种情况。通过该方法建立单胞模型能够清楚地显示内部纱线的四个方向，易于观察内胞分布的周期性，但是，对于任意m×n的三维织物，如果m、n都为偶数，该方法能使宏观结构完全由三种单胞组成，而如果m、n中有一个为奇数，该方法就不能使宏观结构完全地包含三种单胞，而是在织物的角处出现部分的多余正方形部分，如图2.5所示。黄小平等在文献[33]中采用d×d×h作为一个内胞的尺寸，[75]如图2.6所示，洪平等对黄小平的模型作了进一步的细化，提出了微胞体的概念，采用d×d×h/4作为一个内胞的尺寸，并仔细研究了微胞体之间的相互关系。这些单胞划分的新方法能保证编织物完全由三种单胞组成，且使得单胞的分布具有完全对称性，从而可提高细观分析的精度，但是比较难以观察内胞的周期性分布规律，从而给三维编织复合材料的细观力学分析带来了很大的困难，特别是其边界条件的施加问题仍有待探讨。由于经过一个机器循环后引纱器又回到了初始位置，且内部纱线的轴线分布[74]在沿材料纵向且与材料表面成45°的纵向平面上，于是可以取一个花节的长度为单胞的高度，并沿与材料表面成45°的纵向平面剖分宏观预制件（如图2.4所示）得到内胞，然后根据三种单胞的位置关系，在表面和角柱处对宏观织物进行剖分，可以得到面胞和角胞。（a）2×3织物（b）3×3织物图2.5单胞划分方案一16 南京航空航天大学博士学位论文（a）2×3织物（b）3×3织物图2.6单胞划分方案二显然，按照方案一划分单胞，在织物的角处会出现不完整的小正方形，但是通常三维编织结构的m和n的值都比较大，此时不完整部分所占的体积百分比非常小，可以忽略不计，因此采用方案一划分单胞可基本保证细观分析精度。而且，由于本文的重点在于研究三维编织复合材料的细观力学性能，因此采用2d×2d×h作为一个内胞的尺寸按照图2.4和图2.5的方式划分单胞。从图2.4可以看出，按照方案一划分单胞，可以得到两种形式的内胞，3×3排列方式是完整地包含两种形式的内胞、面胞及角胞的最小阵列形式。因此，本文采用对3×3织物进行剖分得到单胞模型是既简单又合理的，通过剖分得到的三种单胞模型同样也适用于其它形式的织物。2.3.3单胞的计算机图形显示宏观模型建立好以后，按照图2.5所示的单胞划分方案，通过布尔运算，可以把宏观模型剖分成细观模型，图2.7显示了通过这种方法得到的三种单胞模型。从图上可以清楚地看出内部纱线的4个方向，并且可以验证出内胞不只包含4根纱线，而是8根。纱线在织物表面和角柱处的曲线路径也清晰可见。(a)内胞(b)面胞(c)角胞图2.7三维编织结构的单胞模型纱线在内胞和角胞处的过渡情况见图2.8。17 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟图2.8纱线在内胞和角胞处的过渡从图2.8，我们不难看出以下几点：①纱线在织物内部是直线路径，到了角柱处就变成了曲线；②在织物纵向，细观结构呈现周期性，周期为一个花节的长度；③在横向，两个相邻的内胞的结构不同，应该作为两种单胞处理。同样地，通过剖分可以显示出纱线在三种单胞处的过渡情况，如图2.9所示。图2.9纱线在三种单胞处的过渡从图2.9，除了可以看出上述几点以外，还可以看到：①纱线在面胞里的路径也是曲线；②引纱器移动6步，纱线在织物的角柱处形成完整的一段，引纱器移动3步，纱线在织物的表面处形成完整的一段。2.4几种常用重要工艺参数2.4.1编织角编织角是纤维束走向和编织方向之间的夹角，其值在织物内部、表面和角柱18 南京航空航天大学博士学位论文处大小不同，分别记为内部编织角γ、表面编织角θ和角柱编织角β。内部编织角是通过控制花节的长度h来确定的，二者之间满足关系式(2.1)，00[74]0且其取值范围为[0,54.7]。本文建立的模型是针对编织角为54.7的情况，此时纤维束之间无相互挤压，因此可以采用圆形截面假设。纤维束在表面和角柱处的路径是曲线，本文用曲线弦的方向与编织方向的夹角来定义θ和β，如图2.10所示。通过几何关系的计算，不难推导出三种编织角满足以下关系式：tanγ=32tanθ（2.4）tanγ=6tanβ（2.5）图2.10三种编织角的关系2.4.2等效直径和纱线填充因子在纤维丝形成纤维束的过程中，由于二者的截面形状都为圆形，因此纤维束中必然存在大量的孔隙，以致在注胶过程中这些孔隙中也会渗进基体。此外，在四步法的编织过程中，每一步都要经过打紧工艺，使得纤维束相互挤压，从而使纤维束的截面形状由圆形变成椭圆形。这些因素都使编织结构的建模和力学分析显得较为复杂。前面的建模过程中把纤维束等效为为实心圆柱体，等效的实心圆[22]的直径我们称之为等效直径，记为D，可通过下式计算得到。f4λD=（2.6）fπρ3其中：λ为纱线线密度（单位：g/m）；ρ为纤维的密度（单位：g/cm）；这两个参量都可以通过实验的方法测得。[22]定义纤维束中纤维的体积分数为纱线填充因子，记为ε，其计算公式如下：19 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟22πDπDffε==（2.7）4Ω4ab其中：Ω为纤维束的实际截面积；a、b为纤维束椭圆截面的长轴和短轴长，它[22]们之间满足关系式（2.8）。a=3bcosγ（2.8）2.4.3纤维体积分数纤维体积分数Vf是编织预制件经过成型后，纤维体积占织物总体积的体积分数。它可以通过三种单胞区域的纤维体积分数与其所占织物总体积的体积分数之积的和计算得到，如下式所示：V=VV+VV+VV（2.9）fiifssfccf其中：Vi、Vs和Vc分别是内胞、面胞和角胞占整体结构的百分比；Vif、Vsf和Vcf分别是三种单胞区域的纤维体积含量。对于任意m×n的编织物，Vi、Vs和Vc可由式（2.10）、（2.11）得到：若m、n都为偶数，计算公式为：2(m−1)(n−1)+2V=i2mn+m+n3(m+n−4)V=（2.10）s2mn+m+n8V=c2mn+m+n若m、n中至少有一个是奇数，计算公式为：2(m−1)(n−1)V=i2mn+m+n3(m+n−2)V=（2.11）s2mn+m+n4V=c2mn+m+nVif、Vsf和Vcf可由公式（2.12）得到：3πV=εif83πcosγVsf=1+3ε（2.12）24cosθ33πcosγV=εcf17cosβ20 南京航空航天大学博士学位论文通过式（2.9）～（2.12）即可计算出整体结构的纤维体积分数，不难发现：影响纤维体积分数的参数主要有m、n和编织角γ。纱线填充因子的取值一般随[76]纤维束的原料而定，对于碳纤维，ε的值一般取为0.785。我们通过计算得到了纤维体积分数随其它参数的变化情况。图2.11显示了内部编织角为30°时，纤维体积分数随m、n的变化情况；图2.12显示了n＝4时，不同m材料的纤维体积分数随内部编织角γ的变化情况。图2.11纤维体积分数Vf随m、n的变化情况（γ＝30°）图2.12纤维体积分数Vf随内部编织角γ的变化曲线（n＝4）通过以上两图可以看出：（1）纤维体积分数随着m、n的增大都逐渐减小，但是减幅也逐渐减小，当m、n的值逐渐增大时，纤维体积分数的值也逐渐趋于平稳；（2）纤维体积分数随着编织角的增大也逐渐减小，但减幅也逐渐增大；（3）21 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟当n不变时，随着编织角的逐渐增大，m对纤维体积分数的影响逐渐减小，从图2.12可以明显看出，当内部编织角γ大于40°时，三条曲线已趋于重合。2.5小结1）从四步法编织工艺入手，可以精确获得纱线编织结构的各种几何参数，利用计算机对三维空间点的精确定位，可以精确而又直观地模拟每一根纱线的运动轨迹，从而建立起较为准确的宏观模型，还可以用动画的形式逼真地模拟每一根纱线地运动情况。2）采用合理的单胞划分方法，对已建立好的宏观模型进行剖分，得到三种细观单胞模型，便于直观清晰地观察织物的内部细观结构、纱线的走向、单胞的周期性。通过这种方法得到的细观模型比直接分析细观结构得到的模型有以下优点：①如果直接从细观入手建立内胞模型，只能包括四根不同方向的纱线，而事实上，由于织物结构的整体性，在内胞的边角处，也存在着部分纱线（见图2.7），因此通过剖分建立的细观模型更加精确；②对于面胞和角胞，可以不用考虑怎样与内胞组成平滑的连接以及自身复杂的曲线形式，通过本法剖分以上情况自然满足,对保持宏观结构的完整性大有益处；③由于单胞的结构只与编织角和纱线直径有关，因此对于任意m×n织物，本法建立的单胞模型都是适合的。3）介绍了几种常用的工艺参数，给出了其计算公式，尤其重点讨论了纤维体积分数随其他工艺参数的变化情况。这些公式及计算结果都为三维编织复合材料的宏、细观力学性能分析提供了准确可靠的数据支持。22 南京航空航天大学博士学位论文第三章三维编织复合材料的等效弹性性能预测在对具有周期性结构的非均质材料进行力学分析时，通常以其代表性体积单元－单胞作为分析对象，从两方面进行研究：（1）采用平均化的思想，预测材料的等效弹性性能，从而把非均质的材料化为均质材料，这是属于宏观力学范畴的刚度分析方法。（2）研究单胞内细观应力的分布情况，再采用适当的强度失效准则对材料进行损伤分析和强度分析，这是属于细观力学范畴的强度分析方法。前一章已讨论了三维编织复合材料细观结构的周期性，本章主要研究其宏观力学性能，即预测其等效弹性性能。关于材料的细观应力的计算及强度准则的建立将在后续章节中进行讨论。3.1基本理论与方法概述对于三维编织复合材料的等效弹性性能，前人已经提出了许多种预测方法，目前比较流行的方法主要有以下三类：1）基于定义的有限元方法复合材料的等效弹性模量的定义如下：假定复合材料纤维和基体的细观应力场σ，应变场为ε，按体积平均值：ijij11σ=σdε=εd（3.1）ij∫vijvij∫vijVVV*则：σ=Dεij（3.2）ijijkl*D即为此复合材料的等效弹性模量，它是一个对称的方阵。通过有限元的ijkl方法求出单胞内应力和应变的体积平均值，再根据公式（3.2）得到材料的等效弹性矩阵。L.Chen等在文献[23]中就采用了这种方法对三维编织复合材料三种单胞的等效弹性模量进行了预测，并提出了一种新的方法—多相元法。在用这种方[23]法预测等效弹性模量时，通常对细观模型施加如下均匀应变边界条件：0u(s)=εx（3.3）iijj00其中：ε称为名义均匀应变场,由复合材料细观力学可知ε即为该胞体的应变体ijij积平均值εij。取六组独立的均匀应变边界条件，通过有限元计算可以分别求出对应的六组应力体积平均值，然后根据式（3.2）即可以得到等效刚度性矩阵，对其求逆可得到等效柔度矩阵，即可较容易地确定各工程弹性常数的值。23 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟这种基于定义的有限元方法的优点是可以得到比其他方法更精确的等效弹性模量的值，其缺点是由于对单胞施加的边界条件并不是实际工况下的边界条件，该方法只对刚度预测有效，无法进行细观应力计算和强度分析。[77][78]2）Voigt上限近似理论和Reuss下限近似理论在计算等效弹性模量时，为了避免在等效体内产生非均匀的应力、应变场，通常采用均匀应变场或均匀应力场假设，而这种与实际情况是不符的。采用均匀应变场假设为Voigt上限近似理论，得到的等效弹性性能是真实等效弹性性能的上限；而采用均匀应力场假设为Reuss下限近似理论，得到的等效弹性性能是真实等效弹性性能的下限。通常大多数文献所采用的刚度体积平均法就是由Voigt上限近似理论演绎而来，而柔度体积平均法则是由Reuss下限近似理论演绎而来。Voigt上限近似理论所得的等效弹性常数与真实解较为接近，为了减小这种上限或下限所带来的误差，KalidinliSR等在文献[79]中提出了一种通过加权平均求解等效弹性常数的方法，所得的结果与实验结果较为接近。这种基于Voigt上限近似理论（或Reuss下限近似理论）的刚度（或柔度）体积平均的方法其优点在于算法比较简单，比较容易用程序实现；但是由于得到的结果是真实结果的上限（或下限），因此计算精度难以保证，而且该方法也只对刚度预测有效，无法进行细观应力计算和强度分析。3）基于小参数渐进展开的均匀化理论均匀化理论是二十世纪七十年代由法国科学家提出的一种具有严格理论基础的数学方法，它属于目前比较流行的多尺度摄动方法的范畴。通过定义宏观和细观两个尺度，引入小参数ε，把位移、应力、应变等任一结构场表示成含有宏观和细观两种尺度坐标的渐进展开式，根据连续介质所满足的基本方程和边界条件，可以得到一系列的摄动方程，通过求解这些摄动方程，不仅可以得到材料的宏观等效性能，还能得到细观应力场的分布，这将对于下一步进行强度分析提供可靠的保证。国内外对这种方法的讨论和应用非常广泛。B.Hassani等在文献[80]—[82]中详细阐述了周期性结构的细观均匀化过程、解析和数值解的均匀化方程及均匀化理论解决最优准则下的拓扑优化问题，得到了许多非常有价值的计算公式和结论；高野直树等在文献[83]中用均匀化理论对平面机织和三维机织复合材料的等效弹性常数及单向拉伸情况下的细观应力进行了数值预报；PeterW.Chung等在文献[84]中用均匀化理论对二维平面机织复合材料及纤维束正交分布的复合材料的等效弹性模量及细观应力分布进行了数值模拟，得到了大量的数值分析结果；国内大连理工大学的刘书田等在文献[85]—[88]中对均匀化过程进行了详细的分析，并运用它解决具有周期性分布的多孔板的弯曲问题、复合材料24 南京航空航天大学博士学位论文的热膨胀问题及梯度功能材料的优化设计问题，得到了大量的刚度和强度的分析[89]结果；谢先海等运用均匀化理论对具有小孔隙的材料在空洞中没有外力的条件下的等效弹性模量的计算公式进行了推导；中国科技大学的冯淼林等在文献[51]—[53]中第一次用均匀化理论预测三维编织复合材料的等效弹性模量，得到了与实验结果较为吻合的预测结果，然而在分析过程中只分析了内胞的平均刚度，忽略了面胞和角胞的影响。与前两种方法相比，均匀化理论的优点主要体现在：①它不仅可以解决刚度问题，还可以解决细观应力问题，有利于进行强度分析；②不管实际工况多么复杂，对于单胞只需施加周期性边界条件，因此它可以解决任意工况下的细观应力分布问题。本章采用用均匀化理论预测三维编织复合材料的等效弹性模量，并讨论了等效弹性模量随工艺参数的变化情况。3.2均匀化过程及有限元列式3.2.1基于小参数渐近展开摄动技术的均匀化过程图3.1所示的复合材料结构是由具有周期重复结构的单胞在空间堆积而成的。（a）复合材料的周期性结构（b）单胞图3.1复合材料的周期性结构及单胞示意图因此，复合材料的不均匀性主要表现在细观，即认为单胞是不均匀的。如果单胞的尺度比整个结构的尺度小得多，设为ε量级，0<ε<<1,由于细观的非均匀性使得在体力和面力的作用下，结构场变量（如位移、应变、应力等）在宏观位置x的非常小的邻域ε内也有很大的变化。因此，可以认为所有的量均依赖于宏观和细观两种尺度坐标x和y=x/ε。也就是说，如果用Ψ(x)表示任一结构场变量，有：Ψ(x)=Ψ(x,y)，y=x/ε（3.4）由于材料细观结构的周期性特征，这些量对细观坐标的依赖关系也具有周期性，即25 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟Ψ(x,y)=Ψ(x,y+kγ)（3.5）其中γ为周期函数的周期,k为任一整数。于是，任一结构场函数对坐标变量的偏微分可表示为：∂Ψ(x)∂Ψ(x,y)1∂Ψ(x,y)=+（3.6）∂x∂xε∂yiii首先将位移场变量表示为ε的渐近展式：()0()12()2xu(x)=u(x,y)+εu(x,y)+εu(x,y)+?,y=（3.7）iiiiε应变和位移、应力和应变满足弹性力学的基本方程（3.8）、（3.9）：1∂ui∂uje(x)=+（3.8）ij2∂x∂xjiσ(x)=Ce(x)（3.9）ijijklij于是，应变场和应力场变量的渐进展开式为：−1(−1)(0)(1)e(x)=εe(x,y)+e(x,y)+εe(x,y)+?（3.10）ijijijij−1(−1)(0)(1)σ(x)=εσ(x,y)+σ(x,y)+εσ(x,y)+?（3.11）ijijijij其中：(0)(0)(−1)1∂ui∂uje(x,y)=+（3.12）ij2∂y∂yji(k)∂u(k)(k+1)∂u(k+1)(k)1∂uij∂uije(x,y)=+++,k=0,1,2,?（3.13）ij2∂x∂x∂y∂yjiji(m)(m)σ(x,y)=C(x,y)e(x,y),m=−1,0,1,2,?（3.14）ijijklkl系统满足如下的平衡方程：σ+b=0（3.15）ij,xji其中b表示体积力。将式（3.11）代入式（3.15）中，将方程左边整理为关于iε的渐近展开式，这样可以得到一系列的控制方程：−2(−1)O(ε):σ(x,y)=0（3.16）ij,yj−1(−1)(0)O(ε)：σ+σ=0（3.17）ij,xjij,yj0(0)(1)O(ε)：σ+σ+b=0（3.18）ij,xjij,yjik(k)(k+1)O(ε)：σ+σ=0,k=1，2，3，?（3.19）ij,xjij,yj26 南京航空航天大学博士学位论文(0)方程（3.16）两边同乘以u(x,y)，并在整个单胞区域内积分，可得：i(0)(−1)uσdγ=0（3.20）∫γiij,yj对方程（3.20）左边进行分部积分，可得：(0)(−1)(0)(0)uσndγ−uCudγ=0（3.21）∫∂γiijj∫γi,yjijklk,yl其中∂γ表示单胞的边界。由于边界具有周期性，方程（3.21）左边第一项为零。又由于刚度矩阵Cijkl是正定的，从而有：(0)u(x,y)=0（3.22）i,yj(0)(0)即u(x,y)=u(x)（3.23）ii方程（3.23）表明位移渐近展开式（3.7）中ε的第0阶项与细观坐标y无关。由方程（3.22）、（3.12）和（3.14）可得：(−1)(−1)σ(x,y)=e(x,y)=0（3.24）ijij然后考虑控制方程（3.17），由（3.24）可知左边第一项为零，则第二项也为零。再把方程（3.13）、（3.14）代入方程（3.15）中，可得：(0)(1)(C⋅(e(u)+e(u)))=0（3.25）ijklklxkly,yj1∂∂其中:e=+。klx2∂x∂xkl(1)用分离变量的方法求解方程（3.25），可将u(x,y)表示为：i(1)kl(0)u(x,y)=χ(y)e(u(x))（3.26）iiklxikl式中χ(y)是γ-周期函数，我们称之为等效位移。i将式（3.26）代入式（3.25）中，可得：kl(C+C⋅e(χ))=0（3.27）ijklijklklyi,yj将式（3.26）、式（3.27）代入式（3.13）中，可得：(0)(0)(0)kle(x,y)=e(u(x))+e(u(x))⋅e(χ(y))（3.28）ijijxklxijy将式（3.28）代入式（3.14）中，可得：(0)kl(0)σ(x,y)=(C+C⋅e(χ(y)))⋅e(u(x))（3.29）ijijklijklijyklx(0)σ(x,y)是应力场展开式（3.11）中的第二项，称为细观应力场。由式ij27 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟（3.24）可知式（3.11）的第一项为零。因此，当小参数ε趋于零时，有：(0)σ(x)→σ(x,y)（3.30）ijij(0)从式（3.29）可以看出，细观应力场σ(x,y)表达式中包括宏观参数x和ij细观参数y，因此它反映了应力在单胞尺度内的波动情况。对式（3.29）的左右两边在单胞区域内取平均值，可得：(0)H(0)σ(x,y)=C⋅e(u(x))（3.31）ijijklklx其中:H1klC=C+C⋅e(χ(y))dγ（3.32）ijkl∫γijklijklijyγH式（3.31）反映了单元尺度内的平均应力和平均应变的关系，C被称为ijkl材料的等效弹性模量，通过求解，可将细观非均匀的材料化为均质材料。从式（3.32）可以看出，积分式的第一项为等效弹性模量的刚度体积平均值，它是平均刚度的上限值，这里称之为刚度体积平均项；积分式的第二项是只与细观尺度坐标有关的值，这里我们称之为摄动项。3.2.2等效位移的等参有限元列式从式（3.29）、式（3.32）可以看出，不管是求材料的等效弹性模量或是细kl观应力，都需先求出等效位移χ(y)。下面我们采用虚位移原理推导等效位移klχ(y)的有限元列式。kl将（3.27）式两边同乘以广义虚位移δχ，并在单胞区域内积分，有：ikl∂Cijklkl∂[]kl∫δχ⋅dγ+∫δχ⋅Ce(χ(y))dγ=0（3.33）iiijklklyiγ∂yγ∂yjj将式（3.33）分部积分得：klklkl∂χiδχCn+δχCn−∫∂γiijkljds∫∂γiijkljds∂yj（3.34）klklkl∂δχ∂δχ∂χ(y)iimCd−Cd=0∫γijklγ∫γijmnγ∂y∂y∂yjjnklkl由于χ和δχ在单胞区域边界上满足周期性边界条件，因此式（3.34）前ii两项为零，故有：klklkl∂δχ∂δχ∂χ(y)iimCd+Cd=0（3.35）∫γijklγ∫γijmnγ∂y∂y∂yjjn28 南京航空航天大学博士学位论文式（3.35）可以化为如下势能泛函的一阶变分：klklklkl∂χi1∂χi∂χm(y)Π(χ)=Cd+Cd=0（3.36）pi∫γijklγ∫γijmnγ∂y2∂y∂yjjn由于编织复合材料的几何形状非常复杂，采用直边的空间单元很难保证其解的收敛性，因此这里采用10节点四面体等参单元，由于这种单元的每条边都是二次曲线，因此可以较好地拟合结构的实际形状。与位移相对应的节点等效位移场为：112233233112χχχχχχ111111kl112233233112χ=χχχχχχ（3.37）i222222χ11χ22χ33χ23χ31χ12333333式（3.36）的有限元离散方程是：ekleKq=P（3.38）其中：111eTK=∫∫∫−−−BCBJdξdηdζ（3.39）111111eTP=−∫∫∫−−−BCJdξdηdζ（3.40）111112233233112qqqqqq111111kl112233233112q=qqqqqq（3.41）i222222q11q22q33q23q31q12333333kle式中：q是等效节点位移矩阵，B是应变矩阵，J是雅克比矩阵，K是单元ieklkl刚度矩阵，P是等效节点载荷矩阵。q与等效位移场χ满足如下关系式：iikleklχ=Nq（3.42）ieN为等参单元的形状函数。klkl应注意的是：χ和q均为3×6矩阵形式，而不是列向量，这就是说细ikl观均匀化的有限元求解过程需6次，每次可以得到q的一列。从式（3.39）可以看出，等效节点载荷矩阵可相当于具有初应力{}jjσ=C,j=1,2,?,6作用的等效载荷矩阵。因此，细观均匀化问题可化为具有0初应力作用及周期性边界条件的空间有限元问题来求解。29 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟3.3三维编织复合材料等效弹性性能计算实例及分析3.3.1实验原始数据及基本假设本算例采用文献[34]中的四组实验数据，基体（TDE-85树脂）为各向同性材料，纤维（T300碳纤维）为横观各向同性材料，基本性能见表3.1。各试件的编织工艺参数见表3.2（其中纤维单丝直径为df；纤维束所含的单丝根数为Nf；内部编织角为γ，纤维体积百分含量为Vf）。表3.1织物所用碳纤维和树脂基体的弹性常数E1/GpaE2/GpaG12/Gpaν12ν23T300碳纤维22013.89.00.20.25TDE-85基体4.54.50.340.34表3.2三维编织复合材料试件的工艺参数3编号df/µmNf/×10γVf17.69210.4427.66410.5437.69420.5847.69480.45在计算过程中采用如下几点假设:（1）纤维束的实际界面是椭圆形，这里为了简化计算，假定纤维束为圆柱形实心体，上一章已经介绍其截面等效直径D可由公式（2.6）得到，而f在已知Nf和df的情况下，D也可由下式确定：fπ2π2（3.43）(D)=(d)⋅N⇒D=N⋅dffffff44[34]（2）当内部编织角较大时，三维编织复合材料的本构关系呈现非线性特性，主要原因是大编织角复合材料基体对整体性能的影响比较大。本章只考虑线性阶段的等效弹性模量问题，因此假设基体的本构关系为线性关系。关于大编织角复合材料的非线性力学分析将在第五章中讨论。（3）由于面胞和角胞内纤维束的路径是曲线，对这两种单胞内的纤维束段进行坐标转换时，在总体坐标系下纤维束的等效弹性常数是不断变化的，这样用有限元的方法定义材料属性相对比较复杂。本文在求解过程中，假设面30 南京航空航天大学博士学位论文[22]胞和角胞处纤维束的路径是螺旋线，把面胞和角胞内的纤维束以机器循环中的每一步为一个步长，拆分成如图3.2所示的折线段，对每一段折线段定义一种材料属性。根据第二章对织物几何结构的分析可知：纤维束在织物表面经过3步形成完整的一段，而在角柱处经过6步形成完整的一段，而面胞中也有部分直纤维束，因此图3.2中面胞中不同方向的折线段共有8条，而角胞中则共有6条。用户视图俯视图仰视图(a)面胞用户视图俯视图(b)角胞图3.2面胞和角胞中纤维束路径的折线段假设3.3.2坐标系的定义及纤维束刚度坐标转换3.3.2.1刚度坐标转换公式设总体坐标系为xyz，局部坐标系为x’y’z’，两种坐标系的方向余弦见表3.3：表3.3局部坐标系为x’y’z’与总体坐标系xyz的方向余弦x’y’z’xl1m1n1yl2m2n2zl3m3n331 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟[90]根据复合材料力学的知识，总体坐标系下的刚度矩阵[C’]和局部坐标系下的刚度矩阵[C]满足以下关系式：T[C′]=[T][C][T]（3.44）σσ其中：222lmn2mn2nl2lm111111111222lmn2mn2nl2lm222222222l2m2n22mn2nl2lm333333333[Tσ]=（3.45）l2l3m2m3n2n3m2n3+m3n2n2l3+n3l2l2m3+l3m2llmmnnmn+mnnl+nllm+lm313131311331133113l1l2m1m2n1n2m1n2+m2n1n1l2+n2l1l1m2+l2m1表示应力空间转换矩阵，称为哈密顿张量转换矩阵。3.3.2.2总体坐标系和局部坐标系的定义由于单胞中纤维束方向的不同，在进行等效弹性常数计算时，首先要定义局部和总体两种坐标系。在用均匀化理论对细观单胞进行有限元分析时，总体坐标系为单胞结构的坐标系，其坐标轴的方向分别记为x、y和z，规定z方向为编织方向，即纵向；单胞在横向两相邻直角边的方向分别为x方向和y方向。由于内胞在横向的两直角边与面胞和角胞的直角边成45°，因此内胞的计算总体坐标系的x和y轴也与面胞和角胞的x和y轴成45°角，如图3.3所示，这在计算过程中应特别注意。（a）三种单胞的相对位置（b）内胞的总体坐标系（c）面胞和角胞的总体坐标系图3.3三种单胞的总体坐标系定义局部坐标系分别为1、2和3，其中规定1方向为纤维束的轴线方向，2、3方向为与纤维束垂直的方向。对于内胞，因为纤维束与总体坐标系下的x（或y）轴垂直，为了计算简便，取2方向与x（或y）方向相同，如图3.4所示；对于面胞和角胞，由于纤维束的路径我们采用了折线假设，定义局部坐标系的1轴沿着各个折线的方向，2和3在垂直折线的平面上可任意选定，因此，面胞和角32 南京航空航天大学博士学位论文胞中的局部坐标系不只一个，根据基本假设3，面胞中的局部坐标系共有8个，而角胞中的局部坐标系共有6个。图3.4内胞的局部坐标系3.3.2.3纤维束刚度矩阵的坐标转换通过以上的方法建立了总体坐标系和局部坐标系后，可以对纤维束的刚度矩阵进行空间坐标转换，下面分三种情况进行讨论：1）内胞中纤维束刚度矩阵的坐标转换对于内胞来说，存在四个方向的纤维束，其中两个与x轴垂直，两个与y轴垂直。这四种纤维束的总体和局部坐标系的方向余弦如表3.4所示。表3.4内胞中纤维束坐标系间的方向余弦123123x010x010ysinγ0cosγy−sinγ0cosγzcosγ0−sinγzcosγ0sinγ（Ⅰ）（Ⅱ）123123xsinγ0cosγx−sinγ0cosγy010y010zcosγ0−sinγzcosγ0sinγ（Ⅲ）（Ⅳ）2）面胞中纤维束刚度矩阵的坐标转换33 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟根据图3.2（a），面胞中各个折线段的方向矢量，即局部坐标1轴的方向矢量，可由各折线段的方向矢量来确定。2、3轴与1轴两两垂直，可任意选定，通过几何分析，可推导出面胞中的8根折线段所在的局部坐标系与总体坐标系的方向余弦，各个方向余弦的值都是关于编织角γ的函数，其表达式非常复杂，这里不再一一列出，可以通过Matlab软件可以计算出任意编织角情况下的各方向余弦的值。3）角胞中纤维束刚度矩阵的坐标转换与面胞坐标转换方法相似，根据图3.2（b），通过几何分析和Matlab计算，同样可计算出角胞中的6根折线段所在的局部坐标系与总体坐标系的方向余弦，这里不再一一列出。3.3.3有限元模型的建立在第二章中我们已经讨论了三维四向编织复合材料的三种单胞模型，采用10节点四面体等参单元对每种单胞的纤维束和基体分别进行网格划分，所建立的有限元模型如图3.5所示（为清晰起见，只显示纤维束单元，下同）：需注意的是：图3.5(a)的内胞模型与图2.7的内胞模型不同，采用的是四个相邻的小内胞模型组成的4d×4d×h的大内胞模型。因为均匀化理论需要满足周期性边界条件，而图3.5(a)采用的内胞模型在空间满足三个方向的周期性，比较易于施加周期性边界条件。（a）内胞（b）面胞（c）角胞图3.5三种单胞的有限元模型示意图34 南京航空航天大学博士学位论文3.3.4固定边界条件的施加均匀化理论的应用前提是细观结构必须满足周期性边界条件。根据图2.4的单胞周期性分布图可以看出：图3.5中的内胞模型具有三个方向的周期性，而面胞和角胞分别只具有两个和一个方向的周期性。在施加边界条件时，既要考虑周期性，又要满足三种单胞在相互交界处的协调性。对于内胞来说，周期性边界条件可以施加等位移边界条件，即认为相对的两个面对应节点的等效位移相等，这用有限元技术比较容易实现，但是很难保证它在与面胞或是角胞的交界处位移保持协调性。因此，本文采用文献[91]所采用的固定边界条件，即认为内胞边界上的等效位移为零。对于面胞和角胞，则令其具有周期性的表面和三种单胞两两相互接触的表面的等效位移为零。事实上，对于面胞和角胞的外部自由表面不用施加约束，只有在其余表面才令其等效位移为零。由于文献[91]的计算结果证明采用固定边界条件与采用等位移边界条件的等效弹性模量的计算结果较为接近，因此，本章中采用了这种固定边界条件。而实事上，若要对进行细观应力分析，边界条件对计算结果有很大影响，采用固定边界条件是使周期性边界条件强制性地满足，但通常不能保证细观应力的计算精度，因而，在下一章中将重点讨论均匀化理论细观边界条件的施加问题。3.3.5整体等效弹性模量及工程弹性常数的计算建立好几种单胞的有限元模型后，对其分别施加不同的周期性边界条件及初H应力载荷，通过均匀化理论可以分别得到内胞、面胞和角胞的等效弹性矩阵C、iHHC和C。然后利用刚度体积平均的方法，根据式（3.46）计算出整个结构的等sc效弹性矩阵：HHHHC=C⋅V+C⋅V+C⋅V（3.46）iisscc其中，V、V和V分别是内胞、面胞、角胞所占的体积百分比。isc通常工程上采用工程弹性常数表示复合材料的弹性性能。由于三维编织结构在横向具有对称性，因此它是横观各向同性材料，具有5个独立的工程弹性常数，分别是横向拉模量Ex、纵向拉模量Ez、横向泊松比μxy、纵向泊松比μzx和纵截面剪切模量Gzx。。可以通过公式（3.47）中等效刚度矩阵[C]、等效柔度矩阵[S]及工程弹性常数的关系计算这5个工程弹性常数。35 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟1µxyµzx−−000EEExxzµxy1µzx−−000EEExxzµµ1−zx−zx000[S]=[C]−1=EzEzEz（3.47）2(1+µ)xy00000Ex100000Gzx100000Gzx本文对文献[34]采用的不同编织角的几种编织物的等效刚度矩阵和工程弹性常数分别进行了计算，等效刚度矩阵的计算结果基本满足横观各向同性材料刚度矩阵的形式，对称性很好。其中四个工程弹性常数的计算结果同文献[34]的实验结果进行比较，如表3.5所示：表3.5工程弹性常数的计算、实验结果比较Ex/GPaEz/GPaµxyµzx编号内部编织角γ/°计算计算实验计算计算实验1219.969.367.30.2870.7160.64824113.630.430.00.3150.6520.60634216.128.127.10.3160.5030.47944818.222.922.10.3820.6530.617为了研究内胞、面胞和角胞对等效弹性模量预测数值的贡献情况，表3.6分别给出了四种材料中三种单胞所占百分比及计算所求得的拉伸模量的数值。表3.6三种单胞所占比例及拉伸模量内胞面胞角胞编号百分比/%E/GPa百分比/%E/GPa百分比/%E/GPazzz167.364.730.876.91.9107.6272.126.926.537.81.469.8367.323.730.835.41.968.6467.318.430.830.71.958.9从表3.5可以看出，通过均匀化理论计算的纵向拉伸模量及泊松比与实验结果比较吻合，随着内部编织角不断增大，拉伸模量随之减小。从表3.6可以看出，36 南京航空航天大学博士学位论文计算三种单胞所得到的等效拉伸模量数值相差较大，第2组数据的内胞所占的比例较大，因此内胞的等效弹性模量与整体的等效弹性模量较为接近。如果不考虑面胞，则计算会带来较大的误差。角胞所占比例甚小，因此计算等效弹性模量时，往往可以忽略，但是在求解细观应力场时，角胞和面胞处往往会出现应力最大值，此时必须正确计算角胞的等效位移。3.3.6均匀化理论同刚度体积平均法的比较前文提到计算等效弹性矩阵的公式（3.32）中的第一项为刚度体积平均项，它是通过voigt上限近似理论得到的结果；而第二项则是与细观参数相关的摄动项。那么究竟摄动项对等效弹性模量的影响有多大呢？本文利用这两种种方法分别计算了材料的等效弹性常数，所得的纵向拉伸模量分别与实验结果比较，如图3.6所示。图3.6均匀化理论同刚度体积平均法Ez预测值的比较从上图可以看出，两种方法的预测结果都比实验值偏大，而显然均匀化理论的预测结果更接近实验结果。尤其是编织角较小时，刚度体积平均法的预测值与实验值存在较大的误差，而均匀化理论预测值与实验值的误差始终保持在很小的范围内，这就充分体现了均匀化理论预测刚度的准确性。3.3.7均匀化理论同有限元方法的比较本文还对文献[23]的三种碳纤维环氧树脂基的三维编织复合材料进行计算,其纤维和树脂的材料常数见表3.7,试样的细观结构参数见表3.8。用本文的方法对上述试件进行计算,与实验和文献[23]通过多相有限元元法的计算结果比较见表3.9。表3.7纤维和树脂的力学参数37 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟E11(Gpa)E22(Gpa)G12(Gpa)G23(Gpa)υ12碳纤维230402414.30.25树脂3.50.35表3.8试样的细观结构参数0试样号编织角()纤维体积含量(%)试样尺寸(mm)11946.620×6×25023047.220×6×25033747.120×6×250表3.9工程弹性常数结果比较试样1试样2试样3实测FEM本文实测FEM本文实测FEM本文E1158.4754.1662.1527.624.2429.518.0517.6519.80E228.339.128.209.188.3310.29E338.349.178.339.168.5910.17G128.859.3513.1814.8913.5514.10G138.419.2212.3413.8712.4113.90G233.083.86.095.548.097.05υ120.720.690.650.780.700.680.780.610.65υ130.690.670.641.000.660.670.740.540.67υ230.310.270.300.260.350.26从表3.9中不难看出，本文的预测结果与用有限元法的预测结果基本吻合，均匀化理论的弹性模量预测结果都比实验结果偏高，而有限元计算结果比实验结果偏低。这是因为：本文在计算整体织物的等效弹性模量时，用到了Voigt刚度体积平均化的思想，会得到刚度的一个上限值；而文献[23]的多相有限元方法由于采用了混合单元，对混合单元通过判断其高斯积分点在纤维还是在基体内来定义单元材料属性，因此，单元尺寸对刚度预测精度影响很大，通常结果会比真实值偏低。我们对表3.5和表3.9进行数据统计和误差分析，可以得出：虽然通过本文的方法得到的等效弹性模量的预测值比实验值都偏高，但是，其误差相对较小，最大误差不超过10％，且许多值都与实验结果非常接近，这也验证了用均匀化理论预测三维编织复合材料等效弹性性能的有效性。3.4工艺参数对工程弹性常数的影响三维编织复合材料的工艺参数对其力学性能有着非常重要的影响，合理设计38 南京航空航天大学博士学位论文工艺参数可以降低研制成本，缩短研制周期，实现材料的优化设计。本节中采用均匀化理论分析编织角和纤维体积分数这两个工艺参数对材料工程弹性常数的影响。分析过程中假定其余参数不变，定性讨论当其中一个工艺参数变化时，相应工程弹性常数的变化趋势。3.4.1编织角对弹性常数的影响其余工艺参数：纤维单丝直径d=7.6µm，纤维束所含的单丝根数fN=9000，纤维体积分数V=0.5。此时，弹性模量和泊松比随内部编织角的ff变化情况如图3.7所示。从图3.7可以看出，随着编织角的增加，纵向拉伸模量Ez逐渐减小，横向拉伸模量Ex逐渐增大，纵向剪切模量Gzx先增大再减小，横向剪切模量Gxy基本不发生变化，且数值很小；随着编织角的增加，纵向泊松比μzx先增大再减小，横向泊松比μxy呈现递减趋势，且内部编织角小于30°时μxy变化非常小。(a)(b)图3.7弹性模量和泊松比随编织角的变化曲线3.4.2纤维体积分数对弹性常数的影响其余工艺参数：纤维单丝直径d=7.6µm，纤维束所含的单丝根数fN=9000，内部编织角γ=30°。此时，弹性模量和泊松比随纤维体积分数的变f化情况如图3.8所示。39 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟(a)(b)图3.8弹性模量和泊松比随纤维体积分数的变化曲线从图3.8可以看出，随着纤维体积分数的增加，所有拉伸模量和剪切模量都呈线性递增趋势，纵向拉伸模量Ez的增幅最大，横向剪切模量Gxy的增幅最小，基本不发生变化；随着纤维体积分数的增加，纵向泊松比μzx呈现递增趋势，但增幅逐渐减小，横向泊松比μxy呈现递减趋势，减幅也逐渐减小。3.5小结本章采用基于小参数渐进展开的多尺度均匀化理论对三维四向编织复合材料的等效弹性模量进行预测。首先建立了包含宏观和细观两个尺度参量的等效弹性模量的均匀化列式，并通过虚位移原理推导出均匀化理论的有限元求解过程，同时还给出了细观应力场的数学表达式；然后对三种类型的单胞施加以不同的周期性边界条件分别进行有限元计算，预测了三维编织复合材料的等效弹性模量，预测结果与实验结果较为吻合，而且比刚度体积平均法的预测结果更为准确；最后讨论了编织角和纤维体积分数这两个工艺参数对工程弹性常数的影响情况，得到了一些重要的规律。大量的事实证明，通过均匀化理论预测三维编织复合材料的等效弹性模量是非常有效的。此外，该方法的优点还主要体现在与下面细观应力分析的方法相贯通，使得三维编织复合材料的强度性能研究更为方便。40 南京航空航天大学博士学位论文第四章三维四向编织复合材料的宏细观应力数值模拟4.1引言前文介绍了采用均匀化理论和有限元技术相结合的方法预测三维四向编织复合材料的刚度，从而得到织物中三种单胞的宏观弹性性能和整体织物的宏观弹性性能。这样，若把单胞或是整个织物看作是均质体，可以很方便地通过宏观力学的方法计算出单胞或是织物的宏观应力。而事实上，三维编织复合材料在细观尺度—单胞范围内具有高度的不均匀性，使得材料的应力在细观很小的尺度内却发生很大的变化，且大量的实验现象也证明：在外载荷作用下，三维编织复合材料并不是沿着截面均匀地破坏，而是从局部首先有纤维或基体发生破坏。因此，仅从宏观角度分析三维编织复合材料的应力是不够的，为了得到材料的细观破坏机理，制定合适的强度失效判据，必须首先正确地计算细观应力分布情况。目前关于三维编织复合材料细观应力计算的理论研究还比较少，多数研究都[69]围绕着几何细观结构和刚度预报而开展。修英姝等采用细观单胞元的方法对三维编织复合材料的应力分布进行了数值模拟，其基本思想是：用结构单胞作为数值分析的离散单元－细观单胞元，建立了这种细观单胞元刚度矩阵的求解方法，然后通过有限元的方法模拟宏观应力分布情况，该方法具有计算简易、工作量小的优点，但计算得到的只是单胞的宏观平均应力，无法反映单胞内部应力的波动[47][68,70]情况。庞宝君、TaoZeng都采用了混合单元的有限元模型对三维编织复合材料进行细观应力计算和损伤分析，这种力学模型考虑了细观力学参数的不均匀性，对于材料的失效分析具有一定的参考价值，然而，在建模过程中没有考虑表面和角柱效应的影响，所采用的内胞模型或“米”字型模型都与织物的真实结构有一定的偏差，且混合单元本身也采用了刚度平均化假设，这些都会对细观应力的计算结果带来一定的影响。事实上，正如上一章所述，采用均匀化理论的优点主要体现在它不仅可以预测材料的宏观刚度，还可以进一步得到细观应力场的分布，更为有利的是，这种方法不受宏观约束和载荷的限制，在细观只需施加周期性边界条件，就可得到单胞内任何位置处各个应力分量的数值。本章主要通过均匀化理论和宏、细观有限元技术相结合的方法对单向拉伸和三点弯曲作用下三维编织复合材料的宏、细观应力分布进行数值模拟。宏观有限元分析可以得到单胞内的宏观平均应力，细观有限元分析则可以反映单胞内部的应力波动情况。然后结合某些文献所得的实验数据和现象，对应力模拟结果进行了分析和讨论。着重讨论了小编织角和大编织角情况下应力分布模式的不同，这41 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟为下一章对大、小编织角复合材料分别制定不同的宏、细观强度判据提供了可靠的依据。4.2均匀化理论求解宏、细观应力的步骤根据第三章中细观均匀化理论过程的公式推导，可以得到细观应力的计算公式，如式（4.1）所示：(0)(kl)(0)σ(x,y)=C+C⋅e(χ(y))⋅e(u(x))（4.1）ijijklijklijyklxklkl式中：C为各组分材料的弹性刚度矩阵，χ(y)为细观等效位移，e(χ(y))为ijklijy(0)(0)细观等效应变，u(x)为宏观位移，e(u(x))为宏观应变。klxkl这里，C只与材料本身的特性有关；χ(y)为细观参量，它与宏观坐标无ijkl(0)关，可以通过基于均匀化理论的细观有限元分析得到；u(x)为宏观参量，可以通过宏观有限元分析得到。根据式（4.1），可以把用均匀化理论求解细观应力可基本归纳为“局部—整体—局部”的分析方法，具体来说，可分为以下六个步骤：1建立细观有限元模型采用前面介绍的合理的单胞划分方法对整体划分单胞，对单胞模型进行网格划分，建立细观有限元模型。2求等效位移矩阵[χ]ii对细观有限元模型施加周期性边界条件，再以向量{σ0}=C,i=1,2,⋯,6对各组分材料分别施加不同的初应力载荷，然后进行有限元求解。由于每次施加的初应力向量为材料弹性常数的一列，因此需要进行六次细观有限元过程，每一次可得到等效位移矩阵[χ]的一列。这一步即为细观均匀化过程，它是求解细观应力的最重要的环节。关于均匀化理论的求解过程我们编制了相应的程序，可用ANSYS10.0软件调用此程序完成不同工艺参数下的细观均匀化过程，这一部分内容将在第六章中作详细介绍。H3求等效弹性矩阵CijklH等效弹性矩阵C的计算公式如式（4.2）所示。ijkl42 南京航空航天大学博士学位论文H1∫klγC=C+C⋅e(χ(y))d（4.2）ijklijklijklijyγγ[χ]带入式（4.2）中，可求出单胞的等效弹性模量。第三章中已经详细介绍了等效弹性模量的计算方法，这一步骤是细观应力计算过程中必不可少的，它是正确求解宏观位移的前提。4单胞的组合根据编织物试件的尺寸，组合三种单胞，形成宏观几何模型，宏观模型可清晰地反映三种单胞的分布情况。(0)(0)5求宏观位移u(x)和宏观应力σ(x)把每种单胞作为均质体看待，用步骤3求得的等效弹性模量表示单胞的材料属性，然后对宏观模型划分网格，施加宏观载荷和约束，进行宏观有限元求解，(0)(0)可得到宏观位移向量{u}。同时，可以得到宏观应力σ(x)的分布情况。(0)6求细观应力σ(x,y)把分别用步骤2和步骤5求得的等效位移和宏观位移带入式（4.1）中，即可求得细观应力。综上所述，所谓的“局部—整体—局部”的分析方法可解释为：前三步通过细观建模、细观均匀化分析和刚度计算，这属于局部分析；第四、五步的宏观有限元分析求得宏观位移和宏观应力属于整体分析；第六步通过均匀化理论求解细观应力公式，求得单胞内任一位置的细观应力，这又回到了局部分析。其中，第1、2、4、5步的细观和宏观有限元分析，可以通过ANSYS10.0软件实现，第3和第6步等效弹性矩阵和细观应力的计算可以通过MATLAB7.0编程实现。由于ANSYS10.0的APDL语言可以调用MATLAB7.0的M文件，因此，可以将ANSYS10.0作为统一的计算平台，编制一套完整的软件，实现三维编织复合材料宏、细观应力分析，这部分内容在第六章中作详细介绍。4.3细观边界条件的施加方法4.3.1周期性边界条件的一般施加方法对于具有周期性结构复合材料的力学性能分析，通常是首先取代表性体积单元—单胞进行细观力学分析，因为材料的不均匀性质主要体现在细观。而单胞在不同的工况下具有不同的边界条件，分析起来较为复杂。单胞边界条件的施加问题已经成为目前复合材料细观力学的研究热点。43 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟一般来说，如果直接通过有限元的方法对局部单胞进行应力分析，可以按照XiaZiHui等在文献[92]中提出的复合材料代表性体积单元周期性边界条件的统一施加方法，该方法认为单胞相对的表面上的位移差等于沿该方向的平均应变与两表面之间距离的乘积，即单胞沿该方向的宏观变形量，具体表达式见式（4.3）。j+j−ju(x,y,z)−u(x,y,z)=c（4.3）iiij式中：j+表示沿x方向的正向，j−表示沿x方向的负向，c表示单胞在x方jjij向的平均变形量。这种方法对于完全周期性分布的复合材料来讲非常有效，可以直接用细观有[93]限元的方法对单胞进行细观力学分析，模拟细观应力的分布情况；王新峰等运用该方法对三维机织复合材料在各种工况下的变形及应力分布进行了模拟，并得到了与实验结果较为接近的等效弹性常数的数值。然而，三维编织复合材料的细观结构比具有一般完全周期性结构复合材料的细观结构要复杂得多。由于纤维束在内部、表面和角柱处的路径不同，三维编织结构不能简单地被看成由具有完全周期性分布的单胞堆积而成，除了应满足位移周期性以外，还必须保证内胞、面胞和角胞相互邻接面的位移协调性，这就给三维编织复合材料的细观边界条件的施加增加了一定的难度。虽然，对于三维编织复合材料，在拉伸、压缩和剪切等情况下单胞的宏观变形通常具有一定的规律性，可以通过文献[92]的方法施加细观边界条件，但在弯曲、扭转等情况下很难发现每个单胞宏观变形的规律性，因此很难用文献[92]所介绍的方法施加细观周期性边界条件。4.3.2基于均匀化理论的周期性边界条件的施加方法均匀化理论采用了局部和整体相互结合的分析方法，在局部分析中，需要通klkl过有限元的方法求解等效位移χ(y)，χ(y)实质上是均匀化理论在求解摄动方程中引入的一个与宏观坐标x完全无关的参量，它并不代表实际的位移，且它与kl单胞的宏观变形量完全无关。对于一般完全周期性结构的复合材料，χ(y)只需00满足的细观边界条件就是等位移边界条件，即对于y∈[0,y],y∈[0,y]11220kl,y∈[0,y]的单胞区域γ，在单胞的边界∂γ上，等效位移χ(y)应满足在三个33[84]方向上位移相等，应满足：jkjk0χ(0,y,y)=χ(y,y,y),(4.4)i23i123jkjk0χ(y,0,y)=χ(y,y,y),(4.5)i13i123jkjk0χ(y,y,0)=χ(y,y,y).(4.6)i12i123对于三维编织复合材料，在运用均匀化理论施加细观边界条件时，除应满足44 南京航空航天大学博士学位论文klklχ(y)的周期性，还应满足内胞、面胞和角胞相互邻接面χ(y)的协调性，因此，关于基于均匀化理论的三维编织复合材料周期性边界条件的施加方法，可分为以下四个方面讨论：（1）内胞满足三个方向的周期性根据第二章中对织物单胞分布规律的分析，如果内胞的尺寸取为8b×8b×h，则内胞属于完全周期性结构，如图4.1所示。图4.1内胞三个方向的完全周期性在具有周期性结构的x、y和z三个方向上，相对的表面应满足周期性边界条件，即等位移边界条件。（2）面胞满足两个方向的周期性对于按照图4.2所示排列的面胞，具有x和z两个方向的周期性。图4.2面胞两个方向的周期性在具有周期性结构的x和z方向上，相对的表面应满足周期性边界条件，即等位移边界条件。45 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟（3）角胞满足一个方向的周期性对于按照图4.3所示排列的面胞，只具有z方向的周期性。图4.3角胞一个方向的周期性同样，在具有周期性结构的z方向上，相对的表面应满足周期性边界条件，即等位移边界条件。kl（4）内胞、面胞和角胞边界满足χ(y)的协调性为了提高细观应力的计算精度，在对三种单胞施加等位移边界条件以后，kl还应保证在内胞、面胞和角胞边界满足χ(y)的协调性，即在胞与胞的接触面kl上应满足χ(y)连续。刚度计算可使用固定边界条件代替等位移边界条件，以使协调性的条件得到强制性的满足，这是因为：刚度计算中采用了平均化的思想，刚度计算精度受边界条件的影响较小。但是这种固定边界条件并不能保证klχ(y)计算的精度，从而会对细观应力的数值计算结果带来很大的影响。因此，kl在计算细观应力时，不宜采用固定边界条件来保证χ(y)的协调性。协调性的满足可用具有强大的二维和三维插值功能的MATLAB7.0软件实现。在用ANSYS10.0对复杂的单胞结构进行多次有限元分析中，很难控制网格尺寸使得内胞、面胞和角胞相互接触面处的节点位置和数目完全一致。MATLAB7.0软件具有强大的曲线、曲面点插值的功能，可以通过在ANSYS10.0中调用MATLAB7.0对接触面内节点位移的双三次多项式插值、接触线上节点位移的三次多项式插值，来满足胞与胞之间的位移协调条件，这部分算法的实现过程将在第六章中作详细介绍。4.4三维编织复合材料的拉伸应力数值模拟单向拉伸属于最简单的工况，而且目前有关三维编织复合材料拉伸实验的文46 南京航空航天大学博士学位论文章很多，因此，本节首先讨论在单向拉伸作用下三维编织复合材料的宏观应力和细观应力的分布情况。[34]编织角是影响三维编织复合材料拉伸力学性能的最重要的工艺参数之一，因此，本节中主要分小编织角情况和大编织角情况，通过两个算例对三维四向编织复合材料的拉伸应力进行数值模拟，其中在算例一中对均匀化理论求解宏、细观应力的步骤进行详述。4.4.1算例一：小编织角情况采用文献[34]的实验数据，基体（TDE-85树脂）为各向同性材料，纤维（T300碳纤维）为横观各向同性材料，弹性常数见表4.1。其余参数为：纤维的内部编织角γ=21°，纤维体积分数V=0.44，纤维单丝直径d=7.6µm，每根纤维束中ff所含的纤维丝的根数N=9000。试件的尺寸为20mm×6mm×120mm。试件两端f所受的拉伸均布载荷为q=20MPa。表4.1织物所用碳纤维和树脂基体的弹性常数E1/GpaE2/GpaG12/Gpaν12ν23T300碳纤维22013.89.00.20.25TDE-85基体4.54.50.340.344.4.1.1单胞的有限元建模对于内胞、面胞和角胞的几何建模，仍然采用第二章中的整体剖分法，建模过程中，对于纤维束的横断面采用椭圆形截面假设，从内胞的直纤维束到面胞、角胞的弯曲纤维束，采用贝塞尔（Bezier）曲线光滑过渡。关于三种单胞的网格划分，采用对纤维束和基体分别划分单元的方法，内胞中的纤维束在经过布尔运算以后形状仍比较规则，因此在ANSYS10.0中可以采用映射网格的划分方法，单元选择SOLID185六面体单元；面胞和角胞中的纤维束轴线为曲线，且在经过布尔运算以后形状很不规则，在ANSYS10.0中用映射网格划分法难以实现，因此仍采用前面常用的的自由网格划分法，单元选用SOLID187四面体10节点单元；基体由于含有大量的孔洞，形状比纤维束更加不规则，因此，可在ANSYS采用自适应网格的划分方法，这样可以保证截面形状变化较大处单元的细化，单元类型选用SOLID187四面体10节点单元。用ANSYS10.0软件对三种单胞划分单元，建立的有限元模型如图4.4所示。47 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟（a）内胞（b）面胞（c）角胞图4.4三种单胞的有限元模型klH4.4.1.2第一个局部分析——求解等效位移矩阵χ(y)及等效弹性常数Cijklkl等效位移χ(y)是细观均匀化理论分析中最重要的一个参数，因为它对细观kl应力的计算精度有很大影响。在计算等效位移χ(y)之前，首先要对已建立好的三种单胞的有限元模型施加周期性边界条件、协调性边界条件和初应力载荷，关于边界条件和初应力载荷的施加方法，前文都已作过详细的介绍，这里不再多述。至此，第一个局部分析的前处理工作已经完成。下面就需要通过六个载荷步的有klkl限元求解计算等效位移矩阵χ(y)的值了，每计算一步可得到χ(y)一列的值。ANSYS10.0软件有着庞大的坐标系体系，不仅可以定义总体坐标系、局部坐标系，还可以定义单元坐标系和节点坐标系。这里，我们为了能得到纤维束沿轴线方向的细观应力，定义局部坐标系的z’轴沿着纤维束的轴线方向，且单元坐标系和局部坐标系一致。图4.5、图4.6和图4.7分别为三种单胞在单元坐标系下计算klkl等效位移χ(y)及所对应的细观变形图，从图中可以清楚地看出χ(y)的周期性规律。（a）mode11（b）mode22（c）mode3348 南京航空航天大学博士学位论文（d）mode12（e）mode23（f）mode13kl图4.5内胞的等效位移χ(y)及变形图（a）mode11（b）mode22（c）mode33（d）mode12（e）mode23（f）mode13kl图4.6面胞的等效位移χ(y)及变形图49 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟（a）mode11（b）mode22（c）mode33（d）mode12（e）mode23（f）mode13kl图4.7角胞的等效位移χ(y)及变形图kl得到等效位移χ(y)以后，可以根据式（4.2），分别计算出内胞、面胞和角HHH胞的等效刚度矩阵，分别记为C、C和C，再利用刚度体积平均法，根据式iscH（3.44），可以计算出整个结构的等效刚度矩阵C，计算结果见式（4.7）（单位：GPa）。11.23.68.000010.43.75.30003.611.28.00003.710.45.3000（4.7）CH=8.08.075.1000CH=5.35.387.6000is0003.4000003.40000009000006.30000009000006.350 南京航空航天大学博士学位论文10.13.64.300010.93.77.00003.610.14.30003.710.97.0000CH=4.34.388.2000CH=7.07.079.3000c0003.2000003.40000005.20000080000005.2000008(0)(0)4.4.1.3宏观分析——组合单胞、求宏观位移u(x)和宏观应力σ(x)在运用均匀化理论求得三种单胞及整体结构的等效弹性常数以后，可以把单胞作为均质体，对三维编织整体结构进行宏观力学分析。根据式(2.2)、式(2.3)，可算出m=19,n=5，即该试件为19×5编织物。由此，可通过恰当的排列方式，组合单胞形成整体织物。三种单胞在整体结构中的分布情况示意图如图4.8所示。图4.819×5编织物中单胞分布示意图需要注意的是：当m或n中有一个为奇数时，在进行单胞划分时，两个对角处会分别有很小的一部分划不到单胞中去，由于这部分的体积很小，对三维编织结构的应力分布影响通常可以忽略不计，这里在进行有限元分析时，我们把这一部分作为基体处理。按照图4.8的单胞排列方式，可以建立宏观几何模型，这里为了消除圣维南效应带来的应力集中，纵向长度取为3h，宏观应力分析结果取中间的一部分。对于比较规则的宏观模型，在ANSYS10.0中可以采用映射网格的方法划分单元，这样得到的单元规整性和对称性都比较好，如图4.9所示。在宏观有限元模型两(0)端施加q=20MPa拉伸均布载荷，可计算得到宏观位移u(x)和宏观拉应力σ，z其中宏观应力σ的分布情况如图4.10所示。z51 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟图4.9小编织角复合材料的宏观有限元模型图4.10小编织角复合材料的宏观拉应力σ分布z从图4.10可以看出，面胞和角胞处的宏观拉应力σ比内胞处略高，其中面z胞处的宏观拉应力最大，这也与文献[7]的结论“表面通常先发生破坏”的结论相符。这也表明了在对三维编织复合材料进行拉伸强度分析时，不能忽略边角效应的影响，因为边角处的纤维束最容易首先被拉断。52 南京航空航天大学博士学位论文(0)4.4.1.4第二个局部分析——求解细观应力σ(x,y)(0)从式（4.1）可知：细观应力σ(x,y)与宏观坐标和细观坐标均有关系，因此，只有在进行完宏观分析后，再重新回到局部，计算细观应力。细观应力的计kl算可按照式（4.1）编制相应的程序，程序中需要读入细观等效位移χ(y)和宏(0)观位移u(x)两个分别来自前面的细观和宏观有限元分析的结果文件。为了便于以后进行强度分析，对于纤维我们通常关心的是其正轴应力分布情况，对于基体则常按照mises准则得到其等效应力分布。关于计算结果的后处理，这里我们采用单元分析法和截面分析法来讨论细观应力的分布情况。（1）单元分析法所谓单元分析法，就是首先建立单元列表，对于纤维束列出空间六个应力分量的值，对于基体列出mises等效应力值，ANSYS10.0中通过ETABLE命令定义单元列表，得到的单元应力值是单元每个节点的应力平均值。单元分析法可以反映整个三维空间结构的应力分布情况，易于判断纤维束和基体的应力危险点的位置，有利于下一步进行损伤刚度折减和强度分析。图4.11、图4.12和图4.13分别为三种单胞中纤维束各细观应力分量分布图（小编织角复合材料的承载能力[34]和变形主要由纤维控制，因此，限于篇幅，这里只对纤维束的应力分布进行数值模拟）。（a）σ1（b）σ2（c）σ353 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟（d）τ（e）τ（f）τ231312图4.11小编织角复合材料的内胞纤维应力分布（γ=21°）（a）σ（b）σ（c）σ（d）τ（e）τ（f）τ123231312图4.12小编织角复合材料的面胞纤维应力分布（γ=21°）（a）σ（b）σ（c）σ（d）τ（e）τ（f）τ123231312图4.13小编织角复合材料的角胞纤维应力分布（γ=21°）从图4.11－图4.13可以看出：小编织角复合材料织物中，纤维束所承受的轴向正应力σ远大于其它应力分量的值，且σ的数值面胞处最大，其次是角胞，11再次是内胞。54 南京航空航天大学博士学位论文（2）截面分析法单元分析法绘制的细观应力云图可以较好的反映整个单胞区域的应力分布情况，但却显得不够直观。为了能更直观、清晰地反映织物的某个横截面上细观应力的波动情况，我们可以采用截面应力分析法来模拟细观应力的分布，这种方法更有利于得到细观应力的分布规律。截面分析法以单胞的某个截面上的节点为分析对象，采用三次立方内插值法模拟织物横截面上的应力波动情况。图4.14、图4.15和图4.16清楚地显示了三种单胞在z=h/2（h为花节高度）处截面示意图以及纵向拉应力σ和米赛斯（mises）等效应力σ的分布图。zvm（a）内胞（b）面胞（c）角胞图4.14三种单胞在z=h/2处截面示意图（γ=21°）（a）内胞（b）面胞（c）角胞图4.15三种单胞在z＝h/2截面处σ分布图（γ=21°）z（a）内胞（b）面胞（c）角胞图4.16三种单胞在z＝h/2截面处mises等效应力σ分布图（γ=21°）vm55 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟比较以上三个图,可以看出：图4.14的纤维束区域恰好对应到图4.15和图4.16出现应力波峰的区域，这是由于在应变基本一致的情况下，纤维束的纵向模量远大于基体的弹性模量，导致其纵向拉伸应力远大于基体的纵向拉伸应力的原因；图4.15和图4.16的网格曲面形状基本相同，纵向拉应力σ比mises等z效应力σ略小一些，这也证明了小编织角复合材料纤维束处于较均匀的受拉状vm态，拉应力以外的其它应力分量数值很小，基本可忽略不计。显然，截面分析法模拟单胞内的细观应力分布比单元分析法更为形象、直观。通过编写程序，还可以模拟出绘制织物内部任意截面处的应力分布图，这对于判断危险截面及危险点的位置也是非常方便的。4.4.2算例二：大编织角情况同样采用文献[34]的实验数据，基体（TDE-85树脂）为各向同性材料，纤维（T300碳纤维）为横观各向同性材料，弹性常数见表4.1。其余参数为：纤维的内部编织角γ=48°，纤维体积分数V=0.54，纤维单丝直径d=7.6µm，每根纤ff维束中所含的纤维丝的根数N=9000。试件的尺寸为20mm×6mm×120mm。试f件两端所受的拉伸均布载荷为q=20MPa。同小编织角情况一样，我们仍然采用均匀化理论和宏、细观有限元相结合的方法，采取“局部－整体－局部”的分析步骤，宏、细观应力的具体求解过程这里不再多述，只给出宏、细观应力的模拟结果。4.4.2.1宏观应力模拟结果图4.17和图4.18分别为这种大编织角复合材料的宏观有限元模型图和宏观应力σ的分布图。z图4.17大编织角复合材料的宏观有限元模型（γ=48°）56 南京航空航天大学博士学位论文图4.18大编织角复合材料的宏观拉应力σ分布（γ=48°）z4.4.2.2细观应力模拟结果采用截面应力分析方法，模拟三种单胞在z=h/2处截面示意图以及纵向拉应力σ和米赛斯（mises）等效应力σ的分布图，见图4.19和图4.20。zvm（a）内胞（b）面胞（c）角胞图4.19三种单胞在z＝h/2截面处σ分布图（γ=48°）z（a）内胞（b）面胞（c）角胞图4.20三种单胞在z＝h/2截面处mises等效应力σ分布图（γ=48°）vm57 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟4.4.2.3结果分析算例二的细观应力计算结果与算例一进行比较，可以发现以下规律：（1）在相同拉伸载荷作用下，大编织角复合材料纤维束的纵向拉应力和mises等效应力都比小编织角情况大许多，因此相同载荷作用下，大编织角复合材料的纤维束先达到自身的拉伸破坏极限，这恰好解释了大编织角复合材料的拉伸强度比小编织角复合材料低的原因；（2）就纵向拉应力σ和mises等效应力σ进行比较，显然大编织角复合材料的zvmσ比σ的增幅比小编织角复合材料要大许多，这主要是由于大编织角复vmz合材料处于较复杂的应力状态，除拉伸正应力以外，其它应力分量的数值也较大；（3）从图4.20可以看出，对于大编织角复合材料，内胞中的纤维束受综合应力因素的影响比面胞和角胞中的纤维束要大，这是由于表面编织角和角柱编织角比内部编织角小很多，因此，大编织角复合材料的表面和角柱处的纤维束还是基本处于较为简单的单向受拉状态，而内部的纤维束则是处于较复杂的多向应力状态。4.4.3编织角对纤维束拉伸应力的影响通过前面的细观应力分布图不难发现：内胞中的纤维束的各个应力分量都处于比较均匀的分布状态。因此，为了研究编织角对三维编织复合材料各应力分量的影响，可以首先取内胞中的一根纤维束为分析对象，绘制出其平均应力随编织角的变化曲线(q=20MPa)。由于纤维束的轴向正应力分量比其它应力分量的数值大得多，因此，为清晰可见，我们用两个图分别表示轴向正应力及其它应力分量随编织角的变化情况，如图4.21和图4.22所示。图4.21纤维束正应力σ随内部编织角的变化情况158 南京航空航天大学博士学位论文图4.22纤维束σ、σ、τ随内部编织角的变化情况2313从图4.21、图4.22，不难得到以下几点规律：（1）在相同拉伸载荷作用下，小编织角复合材料的各应力分量均较小，这说明小编织角复合材料能够承担较大的拉伸载荷，拉伸强度较高。（2）各应力分量的数值随着编织角的增大均呈递增趋势，但轴向正应力σ远1大于其他应力分量的值，尤其编织角越小，σ与其他应力分量相差越大，1因此小编织角复合材料在建立拉伸失效准则时，可只考虑σ的影响，其1他应力分量可以不计。（3）大编织角复合材料的细观应力分布情况比较复杂，其破坏情况受到诸多应力分量因素的影响，在建立拉伸失效准则时，应考虑这些应力分量因素的综合影响。4.5三维编织复合材料弯曲应力数值模拟三维编织复合材料在弯曲情况下的应力分布要比拉伸情况复杂得多，同拉伸情况相似，编织角仍然是影响弯曲破坏模式的主要因素。本节主要结合文献[40]提供的实验数据，按照4.2节中的均匀化理论求解宏、细观应力的步骤对三维编织复合材料的宏、细观应力分布进行数值模拟，通过模拟可以发现应力分布的一些规律。4.5.1算例一：小编织角情况卢子兴等在文献[40]的弯曲实验制备了6种三维编织复合材料试件，试件由碳纤维（T300）和环氧树脂（TDE－85）组成，各组分的性能参数同上面表4.1。实验在北京航空航天大学的MTS万能材料试验机上进行，加载速度为1mm/min，59 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟环境温度为15°C，湿度为18％。弯曲实验采用三点弯曲，试件的形状及尺寸如图4.23所示（单位：mm）。图4.23弯曲试件的形状及尺寸本算例采用文献[40]中的第2种材料，材料相关参数为：体积密度3ρ=1.51g/cm，纤维的内部编织角γ=21°，纤维体积分数V=0.42，纤维单丝直f径d=7.6µm，每根纤维束中所含的纤维丝的根数N=9000。ff4.5.1.1宏观应力及挠度计算通过均匀化理论可以得到整个材料的等效弹性模量，然后把整个三维编织结构看作均质体，并施加以宏观边界条件，通过宏观有限元分析，可以得到宏观应力分布图，图4.24显示了P＝1000N时，试件2（γ=21°）沿织物纵向的应力σz的分布图。图4.24小编织角复合材料的宏观弯曲正应力σ的分布图（γ=21°）z图4.24中弯曲试件上表面中心附近会产生应力集中现象，这是由数值分析过程中加集中力造成圣维南效应所引起的，跟实际情况存在一定的偏差。从图中可以看出：在弯曲情况下，编织复合材料在横截面上的弯曲正应力并不是均匀的，60 南京航空航天大学博士学位论文上侧受压，下侧受拉。弯曲正应力的大小从试件的中性层到上、下表面方向逐渐递增，这与材料力学里弯曲正应力的分布规律是相同的。通过宏观有限元分析还可以得到三维编织弯曲试件的挠度。通过计算，可得这种小编织角复合材料的最大挠度为f=0.6mm，这个结果与文献[40]的载荷max－挠度曲线在P＝1000N时所对应的挠度结果基本一致。4.5.1.2细观应力分布从图4.24可以看出，宏观应力分布图只能反映沿织物纵向z和厚度方向y的应力变化情况，而不能反映x方向的应力变化情况，从而不能准确反映任一横截面上的应力变化情况。通过均匀化理论可以模拟纤维的细观应力分布，得到任一横截面上的应力波动情况。对于三点弯曲情况，按照材料力学的知识，沿z方向中截面是最危险面，沿y方向上下表面为危险截面，因此上、下表面中心点为两个危险点。这里仍采用单胞分析的方法，取危险点附近的三种单胞作为分析对象（其中面胞和角胞受压侧和受拉侧各取一个），采用截面分析法模拟三种单胞按图4.24所示坐标系的z＝0截面处的应力分布图。图4.25为计算所取的三种单胞在宏观模型中的位置分布，图4.26为这三种单胞在z方向中截面的正应力σ分布图（P＝1000N，图中的x和y坐标为细观坐标系下的坐标）。z图4.25弯曲细观应力的分析单胞位置（a）内胞（b）面胞Ⅰ（受压侧）（c）面胞Ⅱ（受拉侧）61 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟（c）角胞Ⅱ（受拉侧）（e）角胞Ⅱ（受压侧）图4.26单胞内弯曲正应力σ的截面分布图（γ=21°）z应力截面图可以清楚地反映出材料横截面上的应力波动情况，纤维束处的应力处于波峰和波谷，基体处的应力比纤维束处小的多。从图中可以看出：在三点弯曲的情况下，应力分布的不均匀性比拉伸情况更加明显，内部区域的应力大小从中性轴沿y轴正、负方向都逐渐递增，其中上面受压，下面受拉，表面和角柱区域的应力远高于内部区域的应力，表面处的弯曲正应力最大，这些结论文献[40]通过分析实验现象所得的结论十分吻合。4.5.2算例二：大编织角情况本算例采用文献[40]中的第4种材料，材料相关参数为：体积密度3ρ=1.56g/cm，纤维的内部编织角γ=48°，纤维体积分数V=0.45，纤维单丝直f径d=7.6µm，每根纤维束中所含的纤维丝的根数N=9000。试件的尺寸和加ff载方式都和算例一的相同。通过计算可知，这种大编织角复合材料单胞在织物横界面上的分布情况也与算例一的相同。4.5.2.1宏观应力及挠度计算与算例一的计算方法相同，我们可以模拟出大编织角复合材料的宏观应力及变形图，如图4.27所示。图4.27大编织角复合材料的宏观弯曲正应力σ的分布图（γ=48°）z62 南京航空航天大学博士学位论文从图4.27可以看出，大编织角复合材料在三点弯曲的情况下，宏观弯曲正应力分布规律与小编织角情况基本相同。通过计算，可得这种小编织角复合材料的最大挠度为f=1.44mm，结果比文献[40]的载荷－挠度曲线在P＝1000N时max所对应的挠度结果偏小（实验曲线上的挠度值约为1.6mm）,造成误差的原因主要是这种大编织角复合材料的实验载荷－挠度曲线在P＝1000N时已呈现明显的非线性特征，而本算例中采用线性假设讨论宏观弯曲应力分布问题，因此计算求得的弯曲变形量略小于材料的实际弯曲变形量。4.5.2.2细观应力分布通过均匀化理论可以模拟出这种大编织角复合材料单胞内任意位置的细观应力。我们仍以图4.25所示所取的单胞为分析对象，采用截面应力分布法模拟三种单胞在z向中截面的正应力σ分布图，如图4.28所示（图中的x和y坐标z为细观坐标系下的坐标）。（a）内胞（b）面胞Ⅰ（受压侧）（c）面胞Ⅱ（受拉侧）（c）角胞Ⅱ（受拉侧）（e）角胞Ⅱ（受压侧）图4.28单胞内弯曲正应力σ的截面分布图（γ=48°）z从图4.28中可以看出：同拉伸情况类似，大编织角复合材料单胞内的细观应力分布规律同小编织角弯曲情况相似，但其弯曲正应力的数值比小编织角情况大得多。63 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟4.5.3编织角对弯曲挠度和弯曲应力的影响通过上面的方法还可以计算出其它编织角复合材料在P＝1000N作用的弯曲挠度和弯曲应力情况。图4.29为Vf＝0.45时，弯曲挠度随编织角的变化曲线。图4.29弯曲挠度随编织角变化曲线由于弯曲正应力的最大值在面胞内，我们取面胞Ⅱ（受拉侧）为分析对象（如图4.30所示），分别计算出P＝1000N时在整体坐标系下的弯曲正应力和mises等效应力的最大值，它们随编织角的变化情况如图4.30所示（其中σz为最大弯曲正应力，σrm为最大mises等效应力）。图4.30纤维束最大弯曲正应力及最大mises等效应力随编织角变化曲线64 南京航空航天大学博士学位论文用同样的方法，我们还可以得到基体的mises等效应力的最大值随编织角的变化规律，如图4.31所示。图4.31基体最大弯曲mises等效应力随编织角变化曲线从图4.29—图4.31中，不难得到以下几点结论：1．弯曲挠度随编织角的增大基本呈线性递增变化，这表明编织角越小，材料抵抗变形的能力越强，塑性越好；2．在相同载荷作用下，纤维束的最大弯曲正应力和最大mises等效应力随编织角的增大逐渐增大，这表明编织角越大，材料的弯曲强度越低；3．纤维束的最大mises等效应力与最大弯曲正应力的差值随着编织角的增大也逐渐增大，编织角较小时，二者的数值基本相等，这表明编织角越大，纤维束除轴向正应力以外的其它应力分量对其影响越大，而编织角较小时，基本可以认为织物表面的纤维束处于单一的受拉或受压状态；4．基体的mises等效应力随编织角的增大逐渐递增，且曲线的斜率也在不断增大。这表明编织角越大，基体对三维编织复材的弯曲力学性能的影响越大，这里在计算过程中采用线性假设，图4.31中当编织角大于37°时，基体的最大mises等效应力已超过其屈服极限（60MPa）甚至超过强度极限（80MPa）。事实上，此时部分基体材料已经发生微屈服，从而使得材料的载荷－位移曲线呈现一定的非线性特征，且编织角越大，这种非线性特征越明显。这一结论与文献[40]通过绘制实验应力应变曲线所得的结论完全相符。65 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟4.6小结同其它理论方法相比，具有严格数学依据的均匀化理论在模拟三维编织复合材料宏、细观应力方面具有相对较强的优势，因此，本章中采用均匀化理论对单向拉伸和三点弯曲作用下三维四向编织复合材料的宏观应力和细观应力分布进行了数值模拟，并着重比较了小编织角情况和大编织角情况应力分布模式的差异，分析了编织角对应力分布的影响规律，得到以下六点重要结论：（1）单向拉伸作用下，小编织角复合材料纤维束的轴向正应力远大于它的其它应力分量的值，基体的mises等效应力较小，这说明此时织物的承载能力主要由纤维束控制，且纤维束一直处于较均匀的受拉状态，可有效地发挥纤维束的抗拉性能；（2）单向拉伸作用下，大编织角复合材料的纤维束的mises应力比其轴向正应力大许多，这说明大编织角复合材料纤维束处于较复杂的应力状态，而难以充分发挥它的抗拉性能；（3）纤维束在织物内部、表面和角柱处的细观拉伸或弯曲应力有一定的差异。通常，面胞处的拉应力或弯曲正应力最大，其次是角胞，再次是内胞，这与拉伸或弯曲实验通常材料表面先发生破坏的现象是完全吻合的；（4）三点弯曲作用下，三维编织复合材料一侧受拉，一侧受压。受拉侧和受压侧应力大小基本相同。受拉侧的细观应力分布规律与单向拉伸情况有一定的相似性，尤其受编织角影响情况与拉伸情况基本相同。但是，由于面胞和角胞正好处于离中性面最远的位置，使得面胞和角胞处的应力与内胞处的应力的差值远大于拉伸情况的差值，因此，可以认为三点弯曲作用下三维编织复合材料的强度主要由面胞和角胞控制。（5）不论是拉伸还是弯曲，编织角都是决定纤维束应力状态的最重要的因素，在其它工艺参数相同且载荷条件也相同的情况下，编织角越大，纤维束的轴向正应力越大，其它应力分量值也越大，纤维束就越容易发生破断。因此，三维编织复合材料的拉伸或弯曲强度都随编织角的增大逐渐减小。（6）通过均匀化理论和宏、细观有限元分析相结合的方法对细观应力进行数值模拟，得到的应力分布的一些规律与文献[34]和文献[40]的实验结论较为吻合，从而验证了该方法模拟细观应力的有效性。以上六点重要结论都将为下一步对三维四向编织复合材料建立合理的强度失效准则提供可靠的依据。66 南京航空航天大学博士学位论文第五章三维编织复合材料强度准则的建立和临界编织角的确定5.1引言三维编织复合材料由于细观结构的不均匀性，使其即使在拉伸、压缩、剪切、弯曲等简单工况下纤维束和基体的应力状态都非常复杂，而且三维编织复合材料的力学性能还受到编织角、纤维体积百分含量、孔隙率、制作和成型工艺质量等诸多因素的影响，使其强度分析显得尤为困难。目前，关于三维编织复合材料的强度研究主要以实验研究为主，许多学者通过拉伸、压缩、剪切、弯曲等实验，用电镜扫描得到的试件断口形状照片分析材料的破坏机理，通过绘制应力－应变曲线得到材料的刚度和强度性能，这些实验数据都为三维编织复合材料的强度理论研究提供了可靠的依据和简便的方法。然而，由于三维编织复合材料的制作工艺复杂、价格昂贵，仅通过实验的方法来预测强度必然造成很大的浪费，因此，如何从理论上制定合理的强度准则具有非常重要的意义。第四章中，我们利用均匀化理论和有限元方法相结合，对拉伸和弯曲情况下编织材料的宏、细观应力分布进行了数值模拟，得到了应力分布的一些重要规律，这对建立合理的强度准则具有重要的参考意义。文献[34]和文献[40]的实验结论指出在拉伸或弯曲情况下，有一个临界编织角，当编织角小于这个值时，应力应变呈线性关系，而大于该值时则呈非线性关系，即具有不同的破坏机理，通过上一章的分析也不难发现编织角对编织复合材料力学性能的影响非常明显。因此，本章中将重点讨论：1.从微观角度出发，对大、小编织角复合材料分别建立不同的强度失效准则，并依据这些准则对三维编织复合材料的拉伸强度进行数值预报；2.确定单向拉伸情况下的临界编织角，用以区分小编织角复合材料和大编织角复合材料。对于压缩、剪切、弯曲等其它工况下的强度问题，同样可按本章介绍的方法进行分析，但是由于这些工况下细观应力的影响因素相对比较复杂，强度预报工作也比较困难，因此本文不作讨论，这部分工作还有待进一步研究。5.2纤维束及基体的微观损伤分析关于复合材料层合板的强度研究方法，主要有宏观分析法和细观分析法两种。前者是不区分纤维和基体，将层合板的每一层都看为均质各向异性材料，通过蔡－希尔、蔡胡、霍夫曼等经典层合板的强度理论判断层板的失效情况；后者则是把纤维、基体及两者之间的界面作为单独的分析对象，认为引起复合材料层合板破坏的原因有纤维破坏、基体破坏和界面分离三种原因，对于纤维、基体和67 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟界面分别采用不同的强度失效判据。然而，通过上一章的应力模拟结果不难发现，三维编织复合材料的应力状态要比层合板复杂的多，不仅纤维和基体的应力相差很多，光是纤维束，在内部、表面或角柱处由于路径发生了变化，应力也会出现较大的差异，直接采用经典层合板的强度理论的宏观分析法或细观分析法，都会给强度分析精度带来较大的误差。因此，我们需要从微观的角度出发，采用有限元的方法，判断织物细观单胞中每个单元的失效情况，从而制定出相对较为合适的强度准则。由文献[34]的实验结论可知：小编织角复合材料的拉伸破坏模式主要是纤维束的瞬间轴向脆断，而大编织角复合材料的受力状况比较复杂，总是从局部先开始产生微裂纹和微孔洞，然后逐渐扩展，直至整个材料的破坏。局部的微裂纹和微孔洞的产生就会导致微观损伤，使得整个材料的有效承载面积减小，力学性能逐渐恶化。因此，从损伤力学的角度去研究编织复合材料的失效情况，更有利于透析材料的破坏机理和建立较精确的强度准则。刚度折减法是用有限元方法处理损伤问题的非常有效的方法，它是通过对于在某种强度判据下已损伤的单元乘以某个系数或是矩阵（也称损伤（折减）系数或损伤（折减）张量），使得该部分单元的刚度矩阵发生衰减而使其承载能力减小，然后在进行下一步的有限元分析。下面我们将详细讨论损伤系数或损伤张量的求法。5.2.1损伤张量的引入引起材料损伤的原因是由于材料内部的微孔洞和微裂纹的存在，使得材料的[94]承载能力降低，从而导致材料力学性能的恶化。损伤的演化和损伤量的大小与材料中微裂纹和微孔洞的发展及其大小、形状、密度和分布有关。对于各向异性体来说，这些微孔洞和微裂纹是有方向性的，并非在各个方向具有相同的性能，因此仅用一个标量描述这些各向异性损伤性能是不够的，有必要引入向量或张量来作为损伤变量。[95]S.Murakami等创立的几何损伤理论定义损伤变量为一个二阶张量，为了定义损伤状态量，他提出了“虚拟（参考）无损构型”的概念，即认为：由于损伤的影响导致作用在损伤面元上载荷的降低，等价于无损伤（虚拟）的情况下材料的有效承载面积的减小。几何损伤理论所定义的损伤张量Ω可以描述复合材料各向异性损伤状态，定义张量(I−Ω)为由于损伤而造成的面元dA的缩减，即*dA=(I−Ω)dA（5.1）其中：I为二阶单位张量，二阶张量Ω是代表材料的各向异性损伤状态的内变量，称为损伤张量。68 南京航空航天大学博士学位论文假设张量Ω是对称的，并且总是具有三个正交的主方向n(i=1,2,3)，其主值iΩ为实变量，因此损伤张量可以表示为正则形：iΩ=∑Ωini⊗ni（5.2）将上式代入式（5.1），不难得出：损伤张量的主值Ω表示在损伤的主平面i内空隙或裂纹的面密度。纤维束可视为具有横观各向同性性质的单向纤维增强复合材料，选取坐标系L－T－Z（L沿纤维方向，T、Z为垂直纤维方向）。设损伤主向与材料主轴方向重合，则损伤张量可写成ΩLΩLTΩZL[]ΩLTZ=ΩΩΩ（5.3）LTTTZΩΩΩZLTZZ*由于有效承载面积的减小，作用在物体上的有效应力（净应力）σ和名义应力σ及损伤张量Ω存在如下关系：*1[−1−1]σ=(I−Ω)σ+σ(I−Ω)（5.4）2*其中：I为二阶单位张量，σ和σ分别为有效应力张量和名义应力张量，都是对称张量。将式（5.3）代入式（5.4），并将应力由张量形式改为向量形式，则坐标系L－T－Z下的有效应力矢量、名义应力矢量及该坐标系下的损伤张量间的关系可表示为：{*}[]{}σ=Tσ（5.5）其中：矩阵[]T的分量T(i,j=1,?,6)是L－T－Z坐标系下损伤张量诸分量的函ij数。各向异性材料有效应力、应变关系为：{*}[]{}（5.6）σ=Dε则名义应力与应变之间的关系可写为：{}[*]{}σ=Dε（5.7）*−1其中：[]D=[]T[D]为损伤状态下的等效弹性矩阵。当应力张量和主损伤张量的主方向重合时，式（5.5）中的各主应力分量可表达为：69 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟***σ=σ/(1−Ω)，σ=σ/(1−Ω)，σ=σ/(1−Ω)（5.8）1112223335.2.2非线性分析纤维束单元和基体单元的损伤判据5.2.2.1损伤与否的判定孔隙的存在使得纤维束中也有部分基体，因此纤维束的刚度和强度都要远低于纤维单丝的刚度和强度。由于碳纤维单丝是横观各向同性材料，而基体为各向同性材料，因此可以把纤维束视为具有横观各向同性的单向复合材料，采用蔡胡（Tsai－Wu）准则判断纤维束单元是否损伤，具体形式如下：222222Fσ+Fσ+Fσ+Fτ+Fτ+Fτ+2Fσσ1112223334423551366121212（5.9）+2Fσσ+2Fσσ+Fσ+Fσ+Fσ=1131323231122331111其中：F=，F=F=，F=，F=F=，1122334425566*2XtXcYtYc(S)(S)111111F=F=−FF，F=−FF，F=−，F=F=−1213112223223312322XXYYtctc式中，X，X和Y，Y分别为单向复合材料的轴向拉、压强度和横向拉、压强tctc*度；S和S分别为单向复合材料的面内剪切强度和面外剪切强度。考虑到单向复合材料力学的研究相对比较成熟，已存在一些与实验结果符合较好的经验公[71]式，故在计算中可以直接引用。具体计算公式如下：*X=σV+σVtffmmEY=1−()V−V1−mσtffEmf2−EmEεmu()1/3X=σV+V或1−VcffEmνff（5.10）EY=1−()V−V1−mσ−cffEmf2S=τV+τVffmmGS*=1−()V−V1−mτffGmf12式中：E和ν分别为按照混合律模型计算的复合材料弹性模量和泊松比，其中纤维泊松比和基体泊松比分别取为0.2和0.34；V和V分别代表纤维体积比和树脂fm体积比；E和E分别代表纤维轴向弹性模量和基体弹性模量(取值为220GPafm*和4.5GPa)；σ为纤维拉伸破坏强度(取值为3GPa)；σ为对应纤维断裂应变的基fm70 南京航空航天大学博士学位论文体拉伸应力(取值为70.6MPa)；σ为树脂拉伸断裂强度(取值为80MPa)；E为mf2纤维的横向弹性模量(取值为13.8GPa)；G为纤维的剪切模量(取值为9.0GPa)；f12−−σ为纤维的压缩强度(取值为2.07GPa)；σ为树脂压缩破坏强度(取值为fm79MPa)；ε为树脂的拉伸断裂应变(取值为1.7%)；τ和τ分别为纤维的剪切mufm破坏强度和树脂的剪切破坏强度(取值分别为943MPa和46MPa)。对于各向同性的基体单元，可用式（5.11）所示的米赛斯（mises）准则判断其是否损伤，具体形式如下：222(222)2(σ−σ)+(σ−σ)+(σ−σ)+6τ+τ+τ=2σ（5.11）122331122331m式中：σ，σ，σ和τ，τ，τ代表材料内一点的应力应力状态，σ为基123122331m体材料的强度。5.2.2.2损伤模式的判定纤维束内部损伤模式可用四种状态来描述，各损伤模式的损伤张量主值与应[47]力状态有关，如表5.1所示。其中，L型损伤表示垂直于纤维束横断面的破断，T及LT型损伤分别表示垂直于T向与纤维方向平行横断面的拉断及剪断，Z型损伤及ZL型损伤分别表示垂直于Z向与纤维方向平行横断面的拉断及剪断，而TZ型损伤表示横向剪切破断。表5.1损伤类型及参量损伤类型应力状态ΩΩΩLTZσL＞0100LσL≤0000σT＞0010T或LTσT≤0000σZ＞0001Z或ZLσZ≤0000σZ＞0011σT＞0σZ≤0010TZσZ＞0001σT≤0σZ≤0000纤维束的损伤模式由式（5.12）所定义的各方向应力分量与其强度比值来确2定,比值最大的方向为破断方向。譬如,如果式（5.12）中的σ/XX值为最大,Ltc则其损伤模式为L型损伤。222222σLσTσZτTZτZLτLT，，，，*，*（5.12）XtXcYtYcYtYcSSS式中：X，X和Y，Y分别为单向复合材料的轴向拉、压强度和横向拉、压强tctc71 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟*度；S和S分别为单向复合材料的面内剪切强度和面外剪切强度。基体材料为各向同性材料，对于损伤的基体单元其等效刚度缩减为0，本文在计算程序中采用文献[70]所介绍的方法进行处理，即将刚度矩阵乘以损伤系数0.001，然后进行下一载荷步的有限元计算。5.3三维编织复合材料拉伸强度准则的建立上一章的细观应力模拟结果和文献[34]的拉伸实验结果都表明：影响三维四向编织复合材料的拉伸破坏模式的最主要因素是它的编织角。当编织角较小时，破坏形式主要是纤维束的轴向拉伸断裂，材料的变形和承载能力主要由纤维控制，应力－应变关系基本保持线性；而当编织角较大时，织物的破坏断面既有脆性断裂特征，也有剪切韧断特征，基体材料对编织复合材料整体力学性能的影响较大,从而使得应力－应变关系表现出明显的非线性性质。因此，在建立三维编织复合材料拉伸强度准则时，需要分小编织角情况和大编织角情况分别进行讨论。5.3.1小编织角复合材料拉伸强度准则的建立5.3.1.1微观损伤判据——最大主应力理论和最大主应变理论由于小编织角复合材料的拉伸破坏形式主要是纤维束的脆性拉断，因此在对小编织角复合材料制定强度准则时，可以不对基体的损伤情况进行判断，对于纤维束单元，本文采用脆性材料的最大主应力理论和最大主应变理论作为其损伤判据，式（5.13）、式（5.14）分别给出了这两种理论的表达形式：σ=X（5.13）1tσ−µ(σ+σ)=X（5.14）123t其中：σ、σ、σ分别为沿编织复合材料三个主轴方向的正应力，X为纤维束123t的拉伸强度，其计算公式可参考式（5.10）中的第一个公式。由于小编织角复合材料断裂应变很小，其破坏通常是纤维束的瞬间脆断，应力－应变关系基本保持线性，因此，不需用较复杂的非线性损伤分析过程对其进行损伤分析，分析过程中可以对其施加载荷增量或位移增量，用最大主应力理论或最大主应变理论计算每个载荷步下纤维的等效应力，当有很少一部分纤维束单元发生破坏，就可认为已达到材料的拉伸极限强度。5.3.1.2宏观强度准则——内胞纤维束损伤判定法微观损伤判据只能判断微观单元的损伤情况，要最终较精确地预报材料的拉72 南京航空航天大学博士学位论文伸强度值，必须建立合理的宏观强度准则。事实上，虽然面胞和角胞处的纤维首先发生破坏，但是由于内胞在整个材料中占最大比例，属于最主要的承载区域，而从图4.2和图4.3中还可以看出内胞中同一个方向的纤维束的应力分布基本相同，因此本文认为：在承受均匀拉伸载荷时，内胞中只要有纤维束单元发生损伤，将会导致整个材料的瞬间脆断，以此作为小编织角复合材料的宏观强度准则。5.3.1.3算例及分析本算例采用文献[34]的实验数据，基体（TDE-85树脂）为各向同性材料，纤维（T300碳纤维）为横观各向同性材料，弹性常数见表4.1。其余参数为：纤维的内部编织角γ=21°，纤维体积分数V=0.44，纤维单丝直径d=7.6µm，每根ff纤维束中所含的纤维丝的根数N=9000。试件的尺寸为20mm×6mm×120mm。f计算过程中纤维束采用椭圆形截面假设，由于纤维束中也存在部分基体，因此可以把纤维束看作单向复合材料，其拉伸强度可按式(5.2)中的第一个公式求得，式中的纤维体积分数V即为纱线填充因子，可按照第二章中式（2.7）求得。f通过计算可得:纤维束的拉伸强度X=1568MPa。t用位移增量法对该材料施加拉伸载荷步进行线性有限元分析（这里取位移增量dz=0.001h），用均匀化理论可以模拟单胞内的纤维应力分布情况，然后分别采用最大主应力理论和最大主应变理论判断每个载荷步下三种单胞内纤维单元的损伤情况。开始的几个位移增量步三种单胞中均无损伤单元出现，然后首先是面胞出现单元损伤，接着是角胞，当计算进行到第9个位移增量步(即平均应变ε=0.009)时，内胞内也已经有纤维单元的等效应力超过纤维束的拉伸强度极限z（1568MPa）。图5.1和图5.2分别是第9个位移增量步时，三种单胞内按照最大主应力理论和最大主应变理论求得的纤维等效应力分布图。（a）内胞（b）面胞（c）角胞图5.1最大主应力理论所对应的等效应力分布图73 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟图5.2最大主应变理论所对应的等效应力分布图从图5.1、图5.2不难发现以下两个规律：（1）按照最大主应力理论和最大主应变理论求得的等效应力基本相同，这说明在单向拉伸状态下，纤维基本可以认为处于单向应力状态，最大拉应力是引起纤维破坏的主要因素，这与文献[34]的实验结果相吻合；（2）内胞处只有很少量的纤维束单元的等效应力已超过纤维束的拉伸强度从而发生破坏，而面胞和角胞处发生破坏的单元明显比内胞处多得多，而且最大主应力出现在面胞处，说明材料表面通常先发生破坏，这与文献[74]所得的结论相符。通过有限元分析，分别计算应力、应变的体积平均值，还可以绘制出这种小编织角复合材料的应力应变曲线（实验曲线参照文献[34]），如图5.3所示。图5.3小编织角复合材料的拉伸应力－应变曲线图5.3拉伸应力－应变计算模拟曲线和实验曲线基本吻合，计算结果证明：这种小编织角复合材料的应力－应变关系基本保持线性。通过上面的计算不难发现，文献[34]的拉伸实验测得的这种小编织角复合材74 南京航空航天大学博士学位论文料拉伸断裂应变为0.891%，本算例中当计算到第9个位移增量步内胞中有纤维单元开始破坏时的平均应变为0.9%，而且通过计算可以得出此时的平均应力为653MPa,此值与实验测得的拉伸强度665MPa也比较接近。因此，采用内胞纤维束损伤判定法作为小编织角复合材料的宏观强度准则是非常有效的。5.3.2大编织角复合材料拉伸强度准则的建立5.3.2.1宏观强度准则的建立方法对于大编织角复合材料，需采用非线性微观损伤分析法对其进行强度分析，其中微观损伤判据在5.2.2节中已作了详细介绍，而宏观强度的准则的建立，本文主要采用了以下两种方法：方法一：绘制失效单元体积比曲线的方法文献[71]提出了一种强度失效判据，即：按照极限强度给单胞体加载，用有限元的方法求出单胞内的应力分布，然后可以把每一个单元作为一种均质材料，用适当的强度理论判断其是否破坏，这样便可计算出编织复合材料破坏时失效单元的体积比Jc，若将此体积比作为临界值，则可由有限次实验的结果给出该临界值随工艺参数（编织角、纤维体积分数等）变化的曲线，如图5.4所示，图中纵坐标J表示任意载荷下失效单元的体积比。利用这一曲线，便可建立编织复合材料强度的一种失效判据。事实上，这种细观失效判据是一种统计的方法，它对于解决复杂细观应力状态下的强度失效问题是非常有效的。图5.4失效单元体积比随工艺参数变化情况示意图用该方法预报大编织角复合材料的拉伸强度，其优点是：由于仅计算在拉伸极限强度下的失效单元体积比，计算方法比较简单，且利用曲线的插值可计算出任意编织角的拉伸强度，算法比较简单，易于编写程序；缺点是：该方法实际上75 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟是一种统计的方法，失效单元体积比的计算结果不仅与编织角有关，还与纤维体积百分比、试件本身的缺陷等其它诸多工艺参数有关，而实验时试验件往往难以保证只随编织角一种参数变化，从而导致统计计算的结果会带来一定的随机性，对计算精度会有较大的影响，且如果实验点取得较少，难以精确地拟合出失效单元体积比曲线，而如果实验点取得较多，曲线的插值又显得没有什么意义，还会造成大量的材料浪费。方法二：基于微观损伤分析的寻找应力－应变曲线中应力下降点的方法第一种分析方法只考虑了编织复合材料在拉伸极限强度下的微观失效情况。而事实上，由于基体的影响和材料的不均匀性，大编织角复合材料的破坏并不是瞬间脆断，采用5.2节中的非线性微观损伤分析法更符合实际情况。通过非线性有限元分析，按照前面的损伤判据对每个载荷步过后已损伤的单元进行刚度折减，然后计算每个载荷步后的平均应力和平均应变，可以得到该种材料的应力－应变关系，在应力应变曲线中我们可以找到应力的下降点，即应力的最大值点，以此时的应力值作为该材料的拉伸强度，以此时的应变值作为该材料的断裂应变。这是因为：由拉伸实验可知，在加载过程中随着应变的增加，虽然由于非线性导致曲线的斜率减小，但应力一直呈递增趋势，当试件破坏时应力突降为零，而我们在有限元分析中，当平均应变达到某个临界值时，会出现应力突然下降的趋势（如图5.5所示），因此可以认为当应变超过该值时，材料将不能继续承载。图5.5计算曲线与实验曲线的比较示意图与方法一相比，通过方法二预报大编织角复合材料的拉伸强度，虽然计算方法比方法一复杂，但是可以不需要借助于大量的实验，而且由于考虑了材料的损伤演变过程，预测结果也比第一种方法更为精确，因此笔者更倾向于用第二种方法预测大编织角复合材料的拉伸强度。76 南京航空航天大学博士学位论文5.3.2.2算例及分析1．采用方法一对大编织角复合材料的拉伸强度进行数值预报。天津工业大学复合材料所和文献[34]中提供了编织角为37°、41°、42°和48°的拉伸实验数据，我们利用编织角为37°、41°和48°的实验数据，用三次多项插值的方法来模拟拉伸极限强度下的失效单元的体积百分比曲线，如图5.6所示。图5.6失效单元的体积百分比随编织角的变化情况根据图5.6，可以得到编织角为42°时的失效单元的体积分数为V=24.7%，然后再预估一个这种材料的拉伸强度值，用均匀化理论和有限元ef分析相结合的方法计算在预估拉伸强度作用下的失效单元体积分数，如果这个值小于临界参数V，则再适量增加载荷继续进行计算，反之，则适量减少载荷继ef续进行计算。通过这种方法，我们可以预测得到编织角为42°织物的拉伸强度"X=336MPa，与文献[34]中通过拉伸实验得到的拉伸强度283MPa比较接近。t2．采用方法二对大编织角复合材料的拉伸强度进行数值预报。采用文献[34]的实验数据，基体（TDE-85树脂）为各向同性材料，纤维（T300碳纤维）为横观各向同性材料，弹性常数见表4.1。其余参数为：纤维的内部编织角γ=48°，纤维体积分数V=0.45，纤维单丝直径d=7.6µm，每根纤维束中ff所含的纤维丝的根数N=9000。试件的尺寸为20mm×6mm×120mm。f采用位移增量法进行非线性有限元分析，为了提高非线性有限元的计算精度，计算时把位移增量取得比较小（这里取dz=0.001h），通过计算每个载荷步后的平均应力和平均应变，可以得到该种材料的应力－应变关系，并与文献[34]的实验曲线进行比较，如图5.7所示。77 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟图5.7编织角为48°复合材料的拉伸应力－应变曲线从图5.7可以找到预测曲线的应力下降点，由此预测出该材料的拉伸强度σ=230MPa，断裂应变ε=1.42%，而通过实验得到的该材料的拉伸强度和btbt断裂应变分别为254MPa和1.585%，预测结果和实验结果还是比较吻合的。同样，通过绘制应力－应变曲线的方法，我们对文献[34]中的其余几种大编织角复合材料的拉伸强度和断裂应变进行了预测，预测结果都与实验结果吻合较好，几种材料的应力－应变曲线见图5.8,拉伸强度和断裂应变的预测结果见表5.2。图5.8几种大编织角复合材料的拉伸应力－应变计算曲线78 南京航空航天大学博士学位论文表5.2拉伸强度和断裂应变的预测值和实验值比较内部编拉伸强度/MPa断裂应变/%织角/°预测值实验值预测值实验值373173061.11.07413052661.21.32422872831.31.46482302541.41.59表5.2中编织角为41°时的拉伸强度预测结果与实验结果相差较大，这是因为这种材料中存在大量的孔隙导致其拉伸强度较低，而在理论计算过程中，我们没有把孔隙的影响考虑进去。显然，通过方法二得到的拉伸强度和断裂应变的预测值比方法一更精确。5.4临界编织角的确定前面的应力模拟和强度预报过程中都提到了大编织角复合材料和小编织角复合材料的概念，然而何为大编织角？何为小编织角？文献[34]指出：三维编织复合材料存在一个表征编织复合材料破坏机制转变的一个重要参量—临界编织角，当三维编织物的编织角小于此临界编织角时，称之为小编织角复合材料，反之则为大编织角复合材料。显然，大编织角复合材料的损伤分析是一个非常复杂的非线性过程，而对小编织角复合材料进行强度分析时则可不必采用如此复杂的非线性分析方法。因此，确定编织角的一个临界值用于区分小编织角复合材料和大编织角复合材料对于制定合理的强度判据、简化强度计算都具有十分重要的意义。下面将对几种不同编织角的复合材料进行拉伸损伤分析，并根据计算经验得出用于表征编织复合材料拉伸破坏机制转变的重要参量—临界编织角。5.4.1拉伸非线性损伤分析在非线性有限元分析过程中，本文采用载荷增量法或位移增量法，按照均匀化理论可以计算出单胞模型中的应力、应变场，对于纤维束，可以用式（5.9）所示的蔡胡准则判断其失效情况，按照5.2节中的各向异性材料损伤模式判据分析其损伤模式；对于基体，可以用式（5.11）所示的米赛斯准则判断其失效情况，对损伤的纤维单元或基体单元按照表5.1列出的损伤类型及参量进行相应形式的刚度衰减，再重新计算单胞内的应力、应变，⋯⋯，如此迭代直到不再有损伤继续发生为止，同时记录下单胞内每个单元的损伤状态，然后再进行下一步载荷步79 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟的有限元计算。为了分析编织角对三维编织复合材料拉伸破坏机制的影响，较准确地确定临界编织角的值，我们分别对编织角为21°、25°、29°、33°、37°、41°和48°的情况进行非线性有限元数值模拟，分析它们的主要破坏模式。其中：对于21°、41°和48°的情况参照文献[34]的实验数据，各组分材料的弹性常数见表3.1；25°和37°的情况参照天津工业大学复合材料所于1998年提供的拉伸实验数据，其余的几种材料无实验数据可参考，这里仅列出假设的工艺参数。七种编织角情况算例的材料原始参数如下表所示：表5.3材料的主要参数3编号编织角/°纤维密度g/cm横截面尺寸/mm纤维体积分数/%实验断裂应变/%1211.5220×6450.8912251.7622×3.6500.8573291.5420×6504331.5420×6505371.7623×4601.0626411.5420×6541.3207481.5620×6451.585对这七种不同编织角的复合材料，按照前面的方法计算每个载荷步下单胞内的细观应力分布，采用相应的损伤判据对其进行非线性有限元分析，加载方式采用位移增量的方式（取位移增量dz=0.001h），计算应力的方法采用细观均匀化理论。首先以第一种材料（编织角为21°）、第四种材料（编织角为33°）第七种材料（编织角为48°）为例，图5.9—图5.11分别为内胞、面胞和角胞内这三种材料的纤维束单元和基体单元的损伤演变情况（其中ε为应变的体积平均值，z图中的体积百分比为各种类型单元体积占单胞总体积的百分比，限于篇幅，这里只显示三种种单胞损伤演变过程中变化幅度比较大的三步的损伤演变图）。80 南京航空航天大学博士学位论文(a)材料1（编织角为21°）内胞单元损伤演变(b)材料4（编织角为33°）内胞单元损伤演变81 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟(c)材料7（编织角为48°）内胞单元损伤演变图5.9内胞单元损伤演变情况（a）材料1（编织角为21°）面胞单元损伤演变82 南京航空航天大学博士学位论文(b)材料4（编织角为33°）面胞单元损伤演变(c)材料7（编织角为48°）面胞单元损伤演变图5.10面胞单元损伤演变情况83 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟（a）材料1（编织角为21°）角胞单元损伤演变（b）材料4（编织角为33°）角胞单元损伤演变84 南京航空航天大学博士学位论文(c)材料7（编织角为48°）角胞单元损伤演变图5.11角胞单元损伤演变情况从图5.9到图5.11的不同编织角情况下三种单胞的单元损伤演变情况，不难发现以下几点规律：（1）材料1内部纤维束单元的损伤形式主要是L型损伤，当平均应变ε=0.009接z近实验断裂应变0.00891时，已损伤纤维束单元的损伤形式还基本都是L形损伤，而此时基体单元的损伤比例比纤维束单元的体积小得多，这充分证明编织角较小时，织物的主要破坏形式是纤维束的轴向拉断；（2）材料7内部纤维束单元的损伤形式主要是L型和ZL型损伤，开始时，ZL型损伤单元的体积比略高于L型损伤单元，但是由于L型损伤单元的增长率高于ZL型损伤单元，因此当平均应变ε=0.015接近实验断裂应变0.01585时，L型z损伤单元的体积比略高于ZL型损伤单元。而此时基体单元的损伤比例要略高于纤维束单元，这证明当编织角较大时，织物的主要破坏原因主要有三种：一是基体的破坏，二是纤维束的轴向拉断，三是纤维束的剪切破断；（3）材料4内部纤维束的损伤形式主要是L型损伤，而基体单元的损伤比例略低于纤维单元，这证明中等编织角复合材料的主要破坏原因有二：一是纤维束的轴向拉断，二是基体的破坏；（4）随着编织角逐渐增大，损伤单元随平均应变的体积增长率呈逐渐下降趋势，这说明编织角越大，织物抵抗变形能力越强，韧性就越好；而编织角越大，纤维束的非L型损伤和基体的损伤占整体损伤的比例也越大，这说明编织角越大，织物的抗拉性能越难以得到发挥，织物的拉伸强度越低。（5）与内胞相比，面胞和角胞内的纤维束损伤单元不仅数量更多，且分布更不均匀，这是由于两种单胞内纤维束应力分布的不均匀性所引起的。此外，与内胞不同的是，编织角较大时，面胞和角胞的单元损伤形式也主要以纤维束单元的L型损伤为主，而受其它损伤类型的影响较小，这是由于纤维束到了表面和角柱处变为曲线，使得表面编织角和角柱编织角比内部编织角小得多，从而使表面和角柱处的纤维束相对内部的纤维束能更好的发挥其抗拉的85 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟性能。（6）对于每种材料，从出现损伤的顺序来说，面胞处都先出现单元损伤，其次是角胞，再次是内胞，这证明整体织物的破坏总是先从材料表面开始。以上得出的六点结论都与文献[34]和文献[47]总结的实验结论和计算结论相吻合，从而为下一步确立合理的临界编织角和建立合理的强度准则都提供了可靠的依据。5.4.2临界编织角的确定VLD为了确定临界编织角，本文引入一个参量λ，令λ=，其中V是在加载LDVD到实验断裂应变时L型损伤的纤维束单元体积，V是在加载到实验断裂应变时D已损伤的纤维束和基体单元的总体积。表5.4为表5.3所列几种不同编织角复材的λ值。表5.4七种材料的λ值材料编号编织角/°λ值1210.9572250.9343290.8834330.7575370.6396410.4787480.335从表5.4可以看出，随着编织角的增大，λ值逐渐减小，纤维束的非L型损伤和基体的损伤所占的比例逐渐增大。这里，为了找到区分大编织角和小编织角的临界编织角，我们需要假设λ的一个临界值。从表5.4中可知，前三种材料的λ值变化不大，而从第三到第四种材料λ值变化比较大，且当编织角为33°时，λ值已比较低，此时基体的损伤和纤维束的非L型损伤已经不能忽略，因此，根据我a们的计算经验，可以假定λ=0.85为λ的临界值，γ=30为临界编织角。当三crcr维编织物的编织角小于此临界编织角时，我们称之为小编织角复合材料，反之，则称为大编织角复合材料。86 南京航空航天大学博士学位论文5.5小结三维四向编织复合材料在单向拉伸情况下的微观破坏形式并不都是简单的纤维束轴向脆断，而是与其编织角的大小有很大的关系。对于三维小编织角复合材料，其破坏形式主要是纤维束的轴向脆断，因此在计算过程中可以忽略对于基体应力的分析，直接用脆性材料的最大主应力理论或最大主应变理论来作为其微观损伤判据；对于三维大编织角复合材料，其破坏形式较为复杂，计算过程中要综合考虑纤维束的各个应力分量及基体应力的影响，对纤维束和基体分别采用蔡胡准则和米赛斯准则作为其微观损伤判据。本章中首先介绍了纤维束和基体的微观损伤分析方法，以此建立大编织角复合材料的微观损伤判据，对小编织角复合材料则以最大主应力理论或最大主应变理论作为其微观损伤判据。然后分别建立了小编织角复合材料和大编织角复合材料的宏观强度准则，其中小编织角复合材料采用内胞纤维束损伤判定法建立其宏观强度准则，强度预报结果与实验结果较为吻合；而对大编织角复合材料分别采用绘制失效单元体积比曲线的方法和绘制应力－应变曲线寻找应力下降点的方法两种方法建立其宏观强度准则，两种方法的强度预报结果都与实验结果较为吻合，但是由于后一种方法强度预测精度更高且不要借助于大量的实验，因此笔者更倾向于用后一种方法预测三维大编织角复合材料的拉伸强度。最后，我们对几种不同编织角的三维编织复合材料的拉伸损伤演变过程进行了数值模拟，找到不同编织角复合材料之间微观破坏模式的差异，然后利用这种差异确定了一个重要的参量——临界编织角，用以区分小编织角复合材料和大编织角复合材料，临界编织角的确定对于制定合理的强度判据、简化强度计算都具有十分重要的意义。本章中关于三维编织复合材料拉伸强度的预测方法也适合于对压缩、剪切、弯曲等其它工况下的强度性能研究。87 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟第六章三维四向编织复合材料力学性能研究的软件开发6.1引言通过前面几章的分析可知：三维四向编织复合材料的宏细观几何建模和力学性能的分析过程是一个非常复杂的过程，尤其是其力学分析过程中用到了均匀化理论，使得在用有限元进行数值分析的过程中的前处理和后处理过程都比较复杂，光用大型的有限元分析软件如ANSYS、NASTRAN等无法实现这些前后处理的功能。然而，许多大型的有限元分析软件如ANSYS、ABAQUS都提供了自身的参数化设计语言，可以对其进行二次开发来实现想要完成的功能。在多种参数化设计语言中，APDL（ANSYSParametricDesignLanguage）语[96]言是ANSYS的参数化设计语言，它是一种类似Fortran的解释性语言，不仅提供了一般程序语言的功能，如参数、宏、标量、向量及矩阵运算、分支、循环、重复等功能，还可以直接访问ANSYS的有限元数据库，另外还通过按钮、对话框等方式提供简单界面定制功能，可以实现参数交互输入、消息机制、界面驱动和运行应用程序等。在参数化分析过程中可以简单地修改其中的参数达到反复分析各种尺寸、不同材料参数、不同载荷大小的多种设计方案或者序列性产品，极大地提高分析效率，减少分析成本。APDL语言中应用最多的是扩展名为mac的宏文件，它用来组织管理ANSYS的有限元分析命令，可以实现参数化建模、施加参数化载荷、求解及参数化后处理结果的显示，从而实现参数化有限元分析的全过程，同时这也是ANSYS批处理分析的最高技术。基于APDL语言强大的编程功能，又可以ANSYS为开发平台实现简单的交互式界面定制功能，笔者通过定制按钮、对话框等方式，并以多个宏文件实现消息响应，来完成三维编织复合材料的细观几何建模、有限元网格划分、边界条件的施加、细观均匀化求解、有效弹性常数的预报、宏细观应力模拟及强度预报等功能等几个主要模块的分析过程，从而形成了一套较完整的基于均匀化理论的三维四向编织复合材料力学性能分析的软件体系。这一软件体系的形成使得读者能更方便、有效地研究三维四向编织复合材料的力学性能，在软件开发平台上，读者只要简单地在对话框中输入编织参数、材料参数、载荷参数等简单的参数信息，点击一些按钮，就可获得单胞的几何模型、有限元模型、刚度、应力、强度等想要得到的信息，从而大大缩短了人们对三维四向编织复合材料的研究及设计开发周期，减少了分析成本。88 南京航空航天大学博士学位论文6.2程序流程图在APDL语言编制软件的过程中，首先应明确三维四向编织复合材料的力学分析过程，这在前面几章中都作了详细介绍，由此我们可以绘制出分析软件的程序流程图。首先，我们绘出刚度预报的分析流程图，如图6.1所示。图6.1三维四向编织复合材料刚度预测流程图89 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟然后,根据图6.1的相关计算数据和第四章求解细观应力的步骤，我们可以绘制出细观应力模拟的分析流程图，如图6.2所示。图6.2三维四向编织复合材料应力模拟流程图得到细观应力模拟结果以后，我们可以对小编织角复合材料和大编织角复合材料按照第五章的方法分别建立微观拉伸强度失效准则，并预报出材料的强度，其分析流程图如图6.3所示。90 南京航空航天大学博士学位论文图6.3三维四向编织复合材料应强度分析流程图91 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟6.3软件界面的定制及开发6.3.1软件界面的设计ANSYS的参数化设计语言APDL可以实现简单的用户化图形交互界面的定制，其主要功能有：定义单参数及多参数输入对话框、调用ANSYS程序已有的对话框、宏中实现拾取操作、创建程序运行进度对话框、定制工具条与缩写等等。其中定制工具条与缩写这一功能应用最为广泛。本文在编制软件的过程中，应用定制用户化工具条按钮来实现对三维编织复合材料的参数输入、细观几何建模、网格划分、细观均匀化求解、等效弹性模量预报、宏观力学分析、细观应力计算和强度分析共八大模块功能，其中每一大模块按照其实现算法又包括相应的子模块，如细观均匀化求解模块共包括初应力载荷施加、周期性边界条件施加及细观有限元计算共三个子模块，刚度、应力和强度三大主要力学分析模块可按照图6.1－图6.3的程序流程图设计相应的子模块。图6.4为三维四向编织复合材料力学分析的工具条按钮界面。图6.4三维四向编织复合材料力学分析的工具条按钮界面图中ANSYSToolbar的前四个按钮是其本身自带的用于数据库文件存储、打开、退出ANSYS和设置图形显示模式的功能。其余按钮是分析软件自定义的按钮，它们表示的含义分别如下：BRAIDPAR—输入编织工艺参数；MATPAR—输入材料参数；BODYCELL—内胞几何建模；FACECELL—面胞几何建模；CORNCELL—角胞几何建模；BODYMESH—内胞网格划分；FACEMESH—面胞网格划分；CORNMESH—角胞网格划分；BBOUN—施加内胞边界条件；FBOUN—施加面胞边界条件；CBOUN—施加角胞边界条件；BINISTR—内胞施加初应力载荷；FINISTR—面胞施加初应力载荷；CINISTR—角胞施加初应力载荷；BHOMO—内胞均匀化求解；FHOMO—面胞均匀化求解；CHOMO—角胞均匀化求解；BODYMOD—内胞刚度预报；FACEMOD—面胞刚度预报；CORNMOD—角胞刚度预报；MOD—整体刚度预报；MACMODEL—宏观几何建模（只包含单胞的模型）；92 南京航空航天大学博士学位论文MACMESH—宏观模型网格划分；MACSOLU—宏观有限元计算；MICSTRESS—求任意界面细观应力分布；TFAIL—拉伸单元损伤分析；TSTREN—拉伸强度预报。6.3.2基本模块的功能和实现算法下面，我们对图6.4中各基本模块中按钮的功能及实现算法作简要介绍。（1）参数输入模块：BRAIDPAR和MATPAR按钮分别启用多参数输入对话框让读者输入编织参数和材料参数，图6.5和图6.6分别为这两个对话框的示意图。图6.5输入编织工艺参数对话框图6.6输入材料参数对话框93 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟（2）细观建模模块：在输入了编织工艺参数和材料参数后，按下BODYCELL、FACECELL和CORNCELL按钮可以分别实现内胞、面胞和角胞的几何建模；三个按钮分别对应三个宏文件，其算法按第二章中细观单胞的建模方法，在ANSYS10.0的GUI窗口按照点—线—面—体—体的布尔操作的顺序进行建模。（3）网格划分模块：BODYMESH、FACEMESH和CORNMESH三个按钮分别实现对内胞、面胞和角胞划分单元的功能，从而建立三种单胞的有限元模型。为了方便定义材料属性，有限元模型在局部坐标系下建立，即单元坐标系与局部坐标系一致，而局部坐标系则沿着纤维束路径的切线方向。内胞中纤维束路径为直线，应力分布比较均匀，因此采用映射网格划分法得到形状较规则的单元；面胞、角胞处的纤维束和基体则采用自适应自由网格划分法，这样可以保证截面形状变化较大处单元的细化。（4）细观均匀化分析模块：这一模块是对三维编织复合材料进行细观力学分析的最主要模块，它又分为三个子模块——初应力载荷的施加、周期性边界条件的施加及均匀化细观有限元计算。前面两个子模块属于细观均匀化分析的前处理部分，而后一个子模块则是六步均匀化的求解和后处理部分。BINISTR、FINISTR和CINISTR三个按钮分别生成内胞、面胞和角胞各6个初应力载荷文件，初应力载荷的施加方法已在前文作过详细介绍。BBOUN、FBOUN和CBOUN三个按钮分别是对三种单胞施加周期性边界条件的过程，周期性边界条件按照第四章所介绍的方法施加，细观等效位移既要满足周期性又要满足协调性，周期性用ANSYS中的COUPLING命令实现，协调性即满足内胞、面胞和角胞边界的位移连续，其基本实现算法是：先用均匀化理论计算内胞的细观等效位移，然后输出内胞和面胞边界的等效位移，以此作为位移约束施加到面胞的相应表面上，用同样的方法把内胞与角胞、面胞与角胞接触面的等效位移作为位移约束施加到角胞的相应表面上，然而由于对三种单胞分网过程是独立进行的，不可能保证单胞的相互接触面处的节点在三种单胞的有限元模型中完全一致，因此，接触面节点等效位移的协调条件可用多项式插值函数来满足，可调MATLAB软件来实现这一功能，而ANSYS的APDL语言可通过如下宏命令语句调用MATLAB的M文件。/sys,D:MATLAB65binwin32MATLAB.exe/rnmbj!内胞、面胞两个邻接面处位移插值BHOMO、FHOMO和CHOMO三个按钮分别是对三种单胞进行细观均匀化计算，可得到细观等效位移，均匀化的前处理过程要在宏文件中首先读入前面定义的初应力载荷文件，完成初应力载荷的施加，对于面胞和角胞，则还要读入使接触面满足协调条件的位移约束条件。然后再进行均匀化求解，每种单胞的均匀化求解过程共需六个载荷步，每计算完一个载荷步用一个扩展名为db的数据库文件进94 南京航空航天大学博士学位论文行保存，其语句格式为“SAVE,body1,db,F:3Dbraid”。均匀化的后处理过程要根据后面的力学分析需要输出多个扩展名为dat的结果文件，这些结果文件主要包括：三种单胞相互接触面的节点位移、单元在计算坐标系（单元坐标系）下的等效位移、待进行应力分析的织物横截面的节点等效位移等等。（5）刚度预报模块：按钮BODYMOD、FACEMOD、CORNMOD和MOD分别可以预测内胞、面胞、角胞和整体织物的工程等效弹性常数。由于各单胞和织物基本都属于横观各向同性材料，因此预报的工程等效弹性常数共5个。通过第三章中等效弹性矩阵的均匀化列式可以得到三种单胞等效刚度矩阵和等效柔度矩阵，通过计算可以得到三种单胞的工程弹性常数，再采用刚度体积平均的方法可以得到整体织物的等效刚度矩阵和等效柔度矩阵，通过计算可以得到整体织物的工程弹性常数。由于MATLAB软件处理矩阵运算效率比较高，这几个按钮对应的宏文件都采用对MATLAB中M文件的调用。MATLAB中计算内胞的工程弹性常数和整体等效弹性常数在程序最后都有结果输出文件，ANSYS通过宏命令读入这些数据文件，并用写字板的打开方式显示刚度的预报结果。于是，当我们点击BODYMOD等按钮时，刚度预报值都会在用户窗口直接显示出来。下图为按照图6.5和图6.6输入的参数所得的内胞刚度和整体刚度的预测结果。图6.7刚度预报输出结果示意图（6）宏观分析模块：在得到单胞及整体的等效弹性常数以后，可以根据材料的外形尺寸和编织参数，组合三种单胞形成宏观模型，并对其施加不同载荷和约束条件进行宏观有限元计算，本软件仅讨论单向拉伸和三点弯曲两种工况的宏观有限元分析。由于APDL语言无法实现单选按钮或复选框的功能，这里仍通过单参数输入对话框来选择拉伸或弯曲的分析类型。当点击MACMODEL按钮时，会弹出如图6.8所示的对话框：95 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟图6.8宏观力学分析类型的选择其中变量Analyis_Type默认值为‘tensile’，即单向拉伸，若要分析三点弯曲情况，则要在对话框的blank处输入‘bend’。若Analyis_Type＝‘tensile’，弹出输入拉伸载荷的对话框，这里采用材料两端施加位移约束的方式；若Analyis_Type＝‘bend’，则弹出输入弯曲载荷的对话框，这里采用上表面施加集中力的方式；这两个对话框如图6.9所示。（a）输入单向拉伸位移约束对话框（b）输入三点弯曲集中力载荷对话框图6.9宏观力学分析输入载荷条件的对话框示意图上图中，h表示花节高度。从图中可以看出，拉伸位移的默认值为0.003h。根据第四章的分析方法，拉伸宏观有限元分析模型织物纵向长度取为3h，因此拉伸位移的默认值实际上使得整个结构的平均应变为0.001。弯曲载荷的加载位置参见图4.23，其默认值为1000N。（7）细观应力模拟模块：细观应力分析方法可以分为单元分析法和截面分析法两种，其中截面分析法可以反映同一个截面上的应力波动情况，较单元分析法更为直观。本软件中只通过MICSTRESS一个按钮来实现截面法绘制细观应力分布图96 南京航空航天大学博士学位论文的功能，点击此按钮首先弹出输入截面位置的对话框（如图6.10所示），然后读取通过细观均匀化分析模块生成的数据库文件，从中获取要分析截面上节点的等效位移的值，然后调用MATLAB中模拟细观应力分布的M文件，绘制出该截面上细观应力分布图。其中单向拉伸情况包括内胞、面胞和角胞的截面应力分布图；三点弯曲情况则包括内胞、受拉侧面胞、受压侧面胞、受拉侧角胞和受压侧角胞的截面应力分布图，所取单胞的位置及应力模拟结果参考第四章的算例。图6.10输入截面位置对话框（8）强度分析模块。强度分析模块又分为拉伸单元损伤分析和拉伸强度值预报两个子模块，分别用TFAIL和TSTREN两个按钮实现这两个功能。强度分析的算法参考第五章和图6.3的流程图，是一个非线性分析的过程。单元损伤分析主要完成非线性有限元求解过程，拉伸强度预报则按照寻找应力下降点法得到拉伸强度和拉伸断裂应变的预测结果。图6.11为按照图6.5和图6.6定义材料，点击TSTREN按钮得到的拉伸强度和拉伸断裂应变的预测结果。图6.11强度预报输出结果示意图97 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟6.4小结采用均匀化理论对三维四向编织复合材料的宏观和细观力学性能进行分析是一个非常繁琐和复杂的过程，因为这其中要考虑宏观和细观有限元相结合、单胞边界条件的施加、非线性微观损伤分析等直接用ANSYS等有限元软件很难解决的问题，而如果直接用Fortran、Matlab等语言编程，由于三维编织四向编织复合材料的几何模型非常复杂，其基本前处理和有限元求解功能也很难实现。如何编制出一套统一的分析软件，可集成均匀化理论的各种分析功能，对于缩短研究周期，提高研究效率具有非常重要的意义。本章中设计了一套可用于三维四向编织复合材料宏细观几何建模和力学性能分析的CAE二次开发软件，该软件以ANSYS10.0为平台，以参数化设计语言APDL为设计开发语言，通过按钮、对话框等方式能很好地实现人机交互式的分析。该软件集成了基于均匀化理论的三维四向编织复合材料的刚度分析、细观应力模拟和强度预报等功能，并提供了参数输入的对话框，使人们可以通过输入几何参数和材料参数，得到各种不同的三维四向编织复合材料的力学性能预报结果，从而可以更为方便地对各种工艺参数对三维四向编织复合材料力学性能的影响进行深入、细致的研究。该软件主要对三维四向编织复合材料的拉伸和弯曲性能进行研究，但是由于其具有良好的可扩展性，随着研究的深入，可以通过添加新的模块，实现其它各种工况下三维四向编织复合材料的宏观和细观力学性能的研究。98 南京航空航天大学博士学位论文第七章总结与展望7.1本文研究工作总结本文以基于小参数渐近展开和摄动技术的多尺度均匀化理论为理论基础，对三维四向编织复合材料的刚度、应力及强度等宏观及细观力学性能进行了有限元数值模拟，并分析了编织角等工艺参数对这些力学参数的影响，该方法从理论上有效地解决了从细观角度解决编织复合材料强度分析的问题，主要可以得到以下几点结论：1．采用对宏观预制件模型剖分的方法建立四步法三维编织复合材料的细观单胞模型，比直接从细观入手建模，更能保证纱线和整体结构的完整性；通过细观几何建模，得到了关于三维四向编织复合材料细观结构的一些重要规律，并且通过计算得到一些常用重要的工艺参数。2．根据摄动原理，引入小参数ε，导出了等效弹性模量和细观应力场的均匀化列式，再根据虚位移原理推导出均匀化理论的有限元列式，然后通过细观和宏观有限元分析相结合，计算三维四向编织复合材料的等效弹性模量和细观应力。通过该方法得到的刚度及强度的数值计算结果都与实验结果较为吻合。3．综合考虑三维四向编织复合材料的非完全周期性和协调性条件，给出了单胞周期性边界条件的统一施加方法。这种细观边界条件的施加方法适合于任意工况下三维四向编织复合材料的均匀化分析,并使得有限元分析具有较高的计算精度。4．影响三维四向编织复合材料的工程弹性常数的因素很多，本文主要讨论了编织角和纤维体积分数这两个工艺参数对工程弹性常数的影响，得到了许多有价值的结论和规律。还对通过均匀化理论的方法和通过其它计算方法预报编织复合材料工程弹性常数的结果进行了比较，分析了各方法的优缺点。5．对单向拉伸和三点弯曲作用下三维四向编织复合材料的细观应力分布进行了数值模拟，应力分布图可清晰、直观地反映出材料在细观尺度内应力的变化情况。通过分析得到了关于三维四向编织复合材料细观应力分布的一些重要规律，这些规律与许多文献中通过实验分析所得的结论较为吻合，从而为建立合理的强度准则提供了可靠的依据。6.分析并比较了不同编织角复合材料细观应力分布状态的差异。小编织角复合材料的承载能力主要由纤维束控制，轴向正应力起主要作用，材料的本构关系基本保持线性；大编织角复合材料的承载能力由纤维束和基体共同控制，纤维束处于较复杂的应力状态，基体的影响使得材料的本构关系呈现明显的99 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟非线性特征。7．对不同编织角的三维编织复合材料的拉伸损伤演变过程进行了数值模拟，找到不同编织角复合材料之间微观破坏模式的差异，利用这种差异确定了一个重要的参量—临界编织角，用以区分小编织角复合材料和大编织角复合材料。对于三维小编织角复合材料和大编织角复合材料分别采用不同的微观损伤判据。8．对小编织角复合材料和大编织角复合材料的拉伸强度进行了数值预报，在预报大编织角复合材料的拉伸强度时分别用了绘制失效单元体积比曲线和绘制应力－应变曲线寻找应力下降点两种方法，两种方法的强度预报结果都与实验结果较为吻合，但是由于后一种方法强度预测精度更高且不要借助于大量的实验，因此笔者更倾向于后一种方法。9．以ANSYS软件为二次开发平台，用APDL语言编制了一套对三维编织复合材料进行几何建模、刚度预报、应力模拟和强度分析的软件，该软件采用人机交互式的界面，可随意增加新的模块，具有良好的扩展性。7.2后续工作展望本文对三维编织复合材料力学性能的研究工作还是初步的，今后工作的方向主要包括：1．进一步对三维编织复合材料在其它工作载荷如压缩、剪切、弯曲、冲击等情况下的细观破坏机理和失效模式进行研究，提出相应的强度失效准则并建立合理的强度预报模型；2．综合考虑除编织角以外其它工艺参数对三维编织复合材料力学性能的影响，如讨论纤维体积分数、孔隙率、纤维束截面形状等因素对三维编织复合材料失效模式的影响，以建立更全面的强度失效准则；3．进一步深入探讨复合材料界面对三维编织复合材料力学性能的影响；4．对三维多向（如五向、六向等）编织复合材料及其它截面形状的编织复合材料（如圆形、工字型等）的编织和成型工艺、几何细观结构和力学性能进行深入研究；5．对三维编织复合材料进行拉、压、剪、弯、扭等强度实验，用扫描电镜观察材料的细观破坏机理。三维编织复合材料作为一种新型高级的复合材料，其优良的力学性能使其在航空、航天、民用等领域具有广阔的应用前景。随着对其力学性能研究的不断深入，本世纪三维编织复合材料的潜力必将得到大规模的开发和利用。100 南京航空航天大学博士学位论文参考文献[1]Y.M.Tarnopolskii,A.V.Roze,G.G.Portnov.SomeNegativeCharacteristicsofFiberReinforcedMaterials.PolymerMechanics,1969,5(1):115-123[2]N.P.Zhmud,V.Y.Petrov&V.N.Shalygin.LaminatedRingsofGlass-PlasticswithExtraRadialReinforcementintheformofSteelPins.PolymerMechanics,1978,18(2):180-183[3]T.Ishikawa.Anti-symmetricElasticPropertiesofCompositePlatesofSatinWeaveCloth.FiberScienceTechnology,1981,15:127-128[4]T.Ishikawa,T.W.Chou.ElasticBehaviorofWovenHybridComposites.JournalofCompositeMaterials,1982,16(2):2-19[5]李嘉禄.立体多向编织结构对复合材料性能的影响.复合材料学报,1996,13(3):71-75[6]吴德隆,沈怀荣.纺织结构复合材料的力学性能.长沙,国防科技大学出版社,1998.12[7]ChouTW,KoFK.TextileStructuralComposites,ElsevierSciencepublishers,1989[8]冯淼林.三维编织复合材料均匀化理论宏细观数值研究.[博士学位论文],合肥,中国科学技术大学,2000[9]PoeCCJr,DexterHB,RajuIS.AreviewofNASAtextilecompositesresearch[R].AIAAPaperNo.97-1321,1997[10]徐焜.三维编织复合材料刚度预报与力学性能试验研究.[硕士学位论文],南京,南京航空航天大学,2004[11]杨桂,敖大新,张志勇等.编织结构复合材料制作、工艺及工业实践.北京,科学出版社,1999:60-63[12]孙慧玉.三维编织复合材料RTM工艺和力学性能.材料工程,1998,(5):34-40[13]肖丽华.异形整体结构复合材料预成型件的编织技术.产业用纺织品,1995,12(4):25-27[14]YangJM,MaCL,ChouTW.FiberInclinationModelofThree-DimensionalTextileStructureComposites,JournalofCompositeMaterials,1986,20(5):472-484[15]ZhouGuangming,LiuWeining,QiaoXin.Astudyonthefabricationof3Dbraidedthcomposites,advancesincompositematerials.In:The8NationalCompositeConference.Beijing:AeronauticalindustryPress,1994.174(inChinese)[16]孙慧玉,吴长春,卞恩荣.三维编织复合材料面内刚度和强度性能研究.复合材料学报,1998,15(4):102-106[17]Joon-HyungByun,ThomasJ.Whitney,Guang-WuDuandTsu-WeiChou,AnalyticalCharacterizationofTwo-StepBraidedComposites,JournalofCompositeMaterials,1991,25:1599-1991101 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟[18]陈利,李嘉禄,冯志海等.三维四向编织复合材料的叠层板理论分析.宇航材料工艺,2003,33(2):29-33[19]KoF.K.,Three-dimensionalfabricsforcompositesanintroductiontothemagnaweavestructure,In:ProceedingsoftheICCM-4JapanSocietyCompositesMaterial,Tokyo/Japan,1982[20]KalidinliS.R.,andAbusafiehA.,Longitudinalandtransversemoduliandstrengthsof3-Dbraidedtextilecompositeswithlowbraidangles,J.Compos.Master.,1996,30:885-905[21]KalidinliS.R.,andFrancoE.,Numericalevaluationofisostrainandweighted-averagemodelsforelasticmoduliof3-Dcomposites[J].CompositesScienceandTechnology,1997,57:293-305[22]L.Chen,X.M.Tao,C.L.Choy.Onthemicrostructureofthree-dimensionalbraidedperforms.CompositesScienceandTechnology,1999,59:391-404[23]L.Chen,X.M.Tao,C.L.Choy.Mechanicalanalysisof3-Dbraidedcompositesbythefinitemultiphaseelementmethod.CompositesScienceandTechnology,1999,59:2383-2391[24]LiW,HammadM,EI-ShiekhA.Structuralanalysisof3-Dbraidedperformsforcomposites,PartI:Thefour-stepperforms.JournaloftheTextileInstitute,1990,81(4):491-514[25]WangYQ,WangASD.Onthetopologicalyarnstructureof3-Drectangularandtubularbraidedperforms.CompositesScienceandTechnology,1994,51(4):575-586[26]WangYQ,WangASD.Microstructure/Propertyrelationshipsinthree-dimensionallybraidedfibercomposites.CompositesScienceandTechnology,1995,53(2):213-222[27]WangYQ,WangASD.Geometricmappingofyarnstructuresduetoshapechangein3-Dbraidedcomposites.CompositesScienceandTechnology,1995,54(4):359-370[28]吴德隆,郝兆平.五向编织结构复合材料的分析模型.宇航学报,1993,7(3):40-51[29]D.L.Wu.Three-cellmodeland5Dbraidedstructuralcomposites.CompositesScienceandTechnology,1996,56(3):225-233[30]黄小平,徐孝诚,孙良新.四向矩形编织物内部纱线上接触点的分布.南京林业大学学报,2000,24(5):33-36[31]黄小平,徐宁光,孙良新.复合材料三维四向矩形编织物表面结构研究.宇航材料工艺,2000,30(3),14-19[32]黄小平,孙良新,徐孝诚.复合材料三维四向矩形编织物角柱结构研究.复合材料学报,2001,18(4):14-19[33]黄小平,孙良新,徐孝诚.四向矩形编织复合材料单胞的新划分方法.纤维复合材料,2001,15(3):15-17102 南京航空航天大学博士学位论文[34]卢子兴,冯志海,寇长河等.编织复合材料拉伸力学性能的研究.复合材料学报,1999,16(3):130-134[35]卢子兴,胡奇.三维编织复合材料压缩力学性能的实验研究.复合材料学报,2003,20(6):67-72[36]陈利,刘景艳,马振杰等.三维多向编织复合材料压缩性能的实验研究.固体火箭技术.2006,29(1):63-66[37]庞宝君,杜善义,严勇等.三维四向碳/环氧编织复合材料剪切力学性能实验研究.实验力学.1999,14(2):209-215[38]庞宝君,杜善义,韩杰才等.三维四向编织碳/环氧复合材料实验研究.复合材料学报.1999,16(4):136-141[39]王波,矫桂琼.三维编织C/SiC复合材料剪切和弯曲性能的实验研究.机械强度.2007,29(1):44-47[40]卢子兴,冯志海,寇长河等.三维编织复合材料的弯曲实验及其失效分析,《复合材料——生命、环境与高技术》,天津大学出版社,2002,10:1070-1073[41]梁军,陈晓峰,庞宝君等.多向编织复合材料的力学性能研究.力学进展.1999,29(2):197-210[42]KregersA.F.,MelbardisYu.G.Determinationofthedeformabilityofthree-dimensionallyreinforcedcompositesbythestiffnessaveragingmethod[J].PolymerMechanics,1978,14(1):3-7[43]MaChang-Long,YangJenn-Ming,ChouTsuWei.Elasticstiffnessofthree-dimensionalbraidedtextilestructurecomposites.ASTMSTP893,Compositematerials:TestingandDesign,Philadelphia,1986:404-421[44]長井謙宏,橫山敦土,前川善一郎等.三次元強化織維複合材料解析手法の研究（第1報）。日本機械學會論文集（A編），1992,58(11):2099-2113[45]長井謙宏,橫山敦土,前川善一郎等.三次元強化織維複合材料解析手法の研究（第2報）。日本機械學會論文集（A編），1994,60(2):514-519[46]庞宝君,曾涛,杜善义.三维多向编织复合材料有效弹性模量的细观计算力学分析.计算力学学报.2001,18(2):231-234[47]庞宝君,杜善义,韩杰才等.三维多向编织复合材料非线性本构行为的细观数值模拟.复合材料学报.2000,17(1):98-102[48]刘振国,陆萌,麦汉超等.三维四向编织复合材料弹性模量数值预报.北京航空航天大学学报.2000,26(2):182-185[49]刘振国,卢子兴,陆萌等.三维四向编织复合材料剪切性能的数值预报.复合材料学报.2000,17(2):66-69103 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟[50]蔡敢为,梁洁萍等.四步法三维编织复合材料弹性性能的数值细观力学分析.湘潭矿业学院学报.1999,14(1):58-61[51]冯淼林,吴长春,孙慧玉.三维均匀化理论预测编织复合材料等效弹性模量.材料科学与工程.2001,19(3):34-37[52]冯淼林,吴长春.基于三维均匀化理论的复合材料本构数值模拟.中国科学技术大学学报.2000,30(6):693-699[53]Miao-LinFeng,Chang-ChunWu,Studyon3-dimensional4-stepbraidedpiezo-ceramiccompositesbyhomogenizationmethod,CompositesScienceandTechnology,2001,61(5):1889-1998[54]AleksanderB.Macander,RogerM.Crane,etc.Fabricationandmechanicalpropertiesofmulti-dimensionally(X-D)braidedcompositematerials.CompositeMaterials,TestingandDesign(7thconf)[C].Philadelphia:STP893,ASTM,1986:422-443[55]WeiLi,AlyElShiekh.Theeffectofprocessesandprocessingparameterson3-Dbraidedperformsforcomposites.33thInternationalSAMPESymposium,1988:104-115[56]LiaoD,TanTM,KoFK.Compressivebehaviorof3-Dbraidedcomposites,PART1:ExperimentalObservations.36thInternationalSAMPESymposium.1991:129-140[57]L.V.Smith,S.R.Swanson.Failureofbraidedcarbon/Epoxycompositeunderbiaxialcompression.JournalofCompositeMaterials,1994(12):1158-1177[58]G.J.Filatovs,R.L.Sadler,etc.FractureBehaviorofA3-DBraidGraphite/EpoxyComposite.JournalofCompositeMaterials,1994(6),526-541[59]李嘉禄,肖丽华,董孚允.立体多向编织结构对复合材料性能的影响.复合材料学报,1996,13(3):71-75[60]黎观生,李建有,徐新.三维编织复合材料拉伸性能试验研究.材料工程,1997(6):39-41[61]刘谦,李嘉禄,李学明.三维编织工艺参数对复合材料拉伸性能的影响.宇航材料工艺,2000,30(1):55-58[62]SuryaR.Kalidindi,AbdelAbusafieh.Longitudinalandtransversemoduliandstrengthoflowangle3-Dbraidedcomposites.JournalofCompositeMaterials,1996,30(8):885-905[63]Vishnus.Avva,HaroldG.Allen,etc.Through-the-thicknesstensionstrengthof3-Dbraidedcomposites.JournalofCompositeMaterials,1994(1):51-67[64]C.M.Pastore,F.K.Ko.Modellingoftextilestructuralcomposites,PartI:Processing-sciencemodelforthree-dimensionalbraiding.JournaloftheTextileInstitute,1990(4):480-490[65]TerryC.Emehel.Predictionofcompressionstrengthoflaminatedandtextilecomposites104 南京航空航天大学博士学位论文usingatowcollapsemodel.AIAA-96-1477-CP:1468-1478[66]L.V.Smith,S.R.Swanson.Effectofarchitectureonthestrengthofbraidedtubesunderbiaxialtensionandcompression.JournalofEngineeringMaterialsandTechnology,1996:478-484[67]Hui-YuSunandXinQiao.Predictionofthemechanicalpropertiesofthree-dimensionallybraidedcomposites.CompositesScienceandTechnology,1997(57),623-629[68]ZengT,WuLZ,GuoLC.Mechanicalanalysisof3Dbraidedcomposites:afiniteelementmodel.CompositesStructures,2004,64(3):399-404.[69]修英姝,陈利,李嘉禄.三维编织复合材料制件的细观单胞元分析.复合材料学报,2003,20(4):63-66[70]TaoZeng,Dai-ningFang,LiMa,etc.Predictingthenonlinearresponseandfailureof3Dbraidedcomposites.MaterialLetters,2004,58(26):3237-3241[71]卢子兴,刘振国,麦汉超等.三维编织复合材料强度的数值预报.北京航空航天大学学报,2002,28(5):563-565[72]孙义林.三维编织复合材料中织物结构的计算机仿真及力学分析.[硕士学位论文],南京,南京航空航天大学,2000[73]孙义林,黄小平,徐宁光等.三维编织结构的计算机图形显示.中国图像图形学报,2000,5(5):443-446[74]黄小平.四向矩形编织复合材料几何细观结构和力学性能研究.[博士学位论文],南京,南京航空航天大学,2001[75]洪平,孙良新,董纪伟.基于微胞体划分的三维编织复合材料有效性能分析.宇航材料工艺,2005,36(5):32-36[76]钟智丽.三维编织结构的几何模型.天津纺织科技,1996,34(3):20-25[77]W.Voigt,[J].Wied.Ann.,1989,38:573-578[78]A.Z.Rauss,[J].Angew.Math.Math.,1929,9:49-58[79]KalidinliSR,FrancoE.Numericalevaluationofisostrainandweighted-averagemodelsforelasticmoduliofthree-dimensionalcomposites.CompositesScienceandTechnology,1997,57:293-305[80]B.Hassani,E.Hinton.AreviewofhomogenizationandtopologyoptimizationⅠ——homo-genizationtheoryformediawithperiodicstructure.ComputersandStructures,1998,69:707-717[81]B.Hassani,E.Hinton.AreviewofhomogenizationandtopologyoptimizationⅡ——analy-ticalandnumericalsolutionofhomogenizationequations.ComputersandStructures,1998,69:719-738105 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟[82]B.Hassani,E.Hinton.AreviewofhomogenizationandtopologyoptimizationⅢ——topologyoptimizationusingoptimalitycriteria.ComputersandStructures,1998,69:739-756[83]均質化法によゐ織物複合材料の三次元マィク口構造設計.日本機械学会論文集,1995,61(5):1038-1043[84]PeterW.Chung,KumarK.Tamma,RajuR.Namburu.Asymtoticexpansionhomogenizationforheterogeneousmedia:computationalissuesandapplications.Composites,2001,32:1291-1301[85]刘书田,程耿东.复合材料应力分析的均匀化理论.力学学报,1997,29(3):306-313[86]刘书田,程耿东,顾元宪等.基于均匀化理论的多孔板弯曲问题新解法.固体力学学报,1999,20(3):195-200[87]刘书田,程耿东.基于均匀化理论的复合材料热膨胀系数预测方法.大连理工大学学报,1995,35(5):451-457[88]刘书田,程耿东.基于均匀化理论的梯度功能材料优化设计方法.宇航材料工艺,1995,26(6):21-27[89]谢先海,廖道训.均匀化理论中等效弹性模量的计算.华中科技大学学报,2001,29(4):44-46[90]周履,范赋群.复合材料力学,高等教育出版社,1991[91]庄守兵,吴长春,冯淼林等.基于均匀化理论的多孔材料细观力学特性数值研究.材料科学与工程,2001,19(4):9-13[92]ZihuiXia,YunfaZhang,FernandEllyin.Aunifiedperiodicalboundaryconditionsforrepres-entativevolumeelementsofcompositesandapplications.InternationalJournalsofSolidsandStructrures,2003,40:1907-1921[93]王新峰,周光明,周储伟等.基于周期性边界条件的机织复合材料多尺度分析.南京航空航天大学学报,2005,37(6):730-735[94]李兆霞.损伤力学及其应用.北京,科学出版社,2002:41-43[95]MurakamiS,OhnoN.Acreepdamagetensorformicroscopiccavities.TransactionsoftheJapanSocietyofMechanicalEngineers(A),1980,46:940-946[96]博弈创作室.APDL参数化有限元分析技术及其应用实例.北京，中国水利水电出版社,2004106 南京航空航天大学博士学位论文致谢本课题的研究是在导师孙良新教授的悉心指导下完成的。导师广博的学识、敏锐的洞察力、严谨的治学态度、科学的工作方法、精益求精的敬业精神和与人为善的良好品德，都给我留下深刻的印象，使我受益匪浅。论文从选题到最后的完成，孙老师付出了大量的心血,不仅在学术上给予我很大的帮助，在生活上也给予我无微不至的关怀，使我十分感动。在此谨向导师表示诚挚的敬意和衷心的感谢！在攻读博士学位期间，得益于南京航空航天大学结构强度研究所良好的科研条件、浓厚的学术氛围和众多老师和同学的关心帮助，使我度过难忘的美好时光，成为我一生宝贵的精神财富。在此向给予我大力支持和帮助的王鑫伟教授、周储伟教授、周光明教授、洪平博士、杨加明博士、张建云博士、王新峰博士、徐焜博士等表示衷心的感谢！在此论文完成之际，感谢我的父母和哥哥给予我学业和生活上的支持！二零零七年三月于南京航空航天大学航空宇航学院107 基于均匀化理论的三维编织复合材料宏细观力学性能的数值模拟在学期间的研究成果及发表的学术论文1.董纪伟，孙良新，洪平.三维编织复合材料宏细观结构的计算机图形模拟.宇航材料工艺，2004，35（4）：21～242.董纪伟，孙良新，洪平.基于均匀化理论的三维编织复合材料等效弹性性能预测.宇航学报，2005，26（4）：482～486EI收录号：053993889873．董纪伟，孙良新，洪平.基于均匀化理论的三维编织复合材料细观应力数值模拟.复合材料学报，2005，22（6）：139～143EI收录号：060496639844．董纪伟，孙良新，洪平.三维四向大编织角复合材料的细观强度分析.材料科学与工程学报，2006，24（4）：539～5425.董纪伟，孙良新，洪平.编织角对三维编织复合材料力学性能的影响.纺织学报（已投稿）6．洪平，孙良新，董纪伟.基于微胞体划分的三维编织复合材料有效力学性能分析.宇航材料与工艺，2005，36（5）：32～36108