工程科学学报.第41卷.第12期：1558-1566.2019年12月 Chinese Journal of Engineering,Vol.41,No.12:1558-1566,December 2019 D0L:10.13374/.issn2095-9389.2019.01.06.001,http://journals.ustb.edu.cn 基于均匀化理论的复合材料安定性分析方法 秦 方)，张乐乐)区，黄松华)，陈耕 1)北京交通大学机械与电子控制工程学院，北京1000442)德国亚探工业大学机械工程材料应用研究所，亚琛52062 ☒通信作者，E-mail:llzhangl@bjtu.edu.cn 摘要周期性非均质复合材料具有微观结构特征，需要均匀化理论进行宏观和微观的多尺度分析来研究其性能表现.针对 其耐久强度性能，应用塑性极限安定下限定理，特别分析了其在长期交变载荷下的安定状态.结合工程应用目标，提出一种 全新的代表性单元边界条件，结合圆锥二次优化算法进行数值计算，可以从材料微结构和组分性能出发，经过弹性应力场求 解确定位移边界载荷数值，最终由优化求解得到复合材料板材的面内塑性性能容许域.所求得的应力域以单向应力为基，可 根据结构宏观的单向应力状态变化幅值直接进行安定状态与否的判定.通过文中的多个算例，验证了所编写的软件及计算 流程的可行性及数值准确性，展示了该方法在工程模型中的应用场合和工程实践意义 关键词复合材料：极限安定分析：均匀化理论；代表性单元：圆锥二次优化 分类号0344.5 Shakedown analysis method for composites based on homogenization theory QIN Fang,ZHANG Le-le HUANG Song-hua,CHEN Geng 1)School of Mechanical,Electronic and Control Engineering,Beijing Jiaotong University,Beijing 100044,China 2)Institute for Materials Applications in Mechanical Engineering,RWTH Aachen University,Aachen 52062,German Corresponding author,E-mail:llzhangl@bjtu.edu.cn ABSTRACT Direct methods of plastic analysis are widely used in composites analysis to determine material strength for safety assessment or lightweight optimization design.Multi-scale processing of periodic heterogeneous composite material is needed due to its existing of microstructure.The standard method is to determine the macroscopic properties from the calculation results of microcosmic representative volume elements(RVEs)by using the homogenization theory.However,in current practice,there are some disadvantages of transforming the micro strain domain to the macro stress shakedown domain when considering multiple external loads.The domain cannot fully demonstrate the shakedown condition,and it is impossible to evaluate a known loading combination only from the knowledge of whether the load leads to the shakedown state.To overcome this disadvantage,a new comprehensive approach was proposed to enhance endurance limit strength of composites under variable loads for long term.Considering the example of in-plane strength analysis,for microcosmic RVEs,a new set of boundary condition was defined in the form of uniform strain.The boundary condition was derived from the elastic response under unit loads by using Hook's law and stiffness matrix.The resulting elastic stress field was used later for plastic shakedown analysis.Based on the lower bound theorem of plastic mechanics,optimization programming for load factor was performed,and after proper mathematical reformulation,the conic quadratic optimization problem could be solved efficiently.Macro-stress shakedown domain can be obtained after scale-transformation of the RVE results.The bases of this stress domain are unidirectional stress in geometry space.The stress amplitude of a structure can be evaluated by this domain for determining the shakedown state in a simple and practical manner.Further,changes in the boundary condition of RVE do not affect the limit and elastic analysis.Finally,few numerical examples were presented for verification and illustration.This approach can be expanded to three 收稿日期：2019-01-06 基金项目：国家自然科学基金资助项目(51475036)
基于均匀化理论的复合材料安定性分析方法 秦 方1)，张乐乐1) 苣，黄松华1)，陈 耕2) 1) 北京交通大学机械与电子控制工程学院，北京 100044 2) 德国亚琛工业大学机械工程材料应用研究所，亚琛 52062 苣通信作者，E-mail：llzhang1@bjtu.edu.cn 摘 要 周期性非均质复合材料具有微观结构特征，需要均匀化理论进行宏观和微观的多尺度分析来研究其性能表现. 针对 其耐久强度性能，应用塑性极限安定下限定理，特别分析了其在长期交变载荷下的安定状态. 结合工程应用目标，提出一种 全新的代表性单元边界条件，结合圆锥二次优化算法进行数值计算，可以从材料微结构和组分性能出发，经过弹性应力场求 解确定位移边界载荷数值，最终由优化求解得到复合材料板材的面内塑性性能容许域. 所求得的应力域以单向应力为基，可 根据结构宏观的单向应力状态变化幅值直接进行安定状态与否的判定. 通过文中的多个算例，验证了所编写的软件及计算 流程的可行性及数值准确性，展示了该方法在工程模型中的应用场合和工程实践意义. 关键词 复合材料；极限安定分析；均匀化理论；代表性单元；圆锥二次优化 分类号 O344.5 Shakedown analysis method for composites based on homogenization theory QIN Fang1) ，ZHANG Le-le1) 苣 ，HUANG Song-hua1) ，CHEN Geng2) 1) School of Mechanical, Electronic and Control Engineering, Beijing Jiaotong University, Beijing 100044, China 2) Institute for Materials Applications in Mechanical Engineering, RWTH Aachen University, Aachen 52062, German 苣 Corresponding author, E-mail: llzhang1@bjtu.edu.cn ABSTRACT Direct methods of plastic analysis are widely used in composites analysis to determine material strength for safety assessment or lightweight optimization design. Multi-scale processing of periodic heterogeneous composite material is needed due to its existing of microstructure. The standard method is to determine the macroscopic properties from the calculation results of microcosmic representative volume elements (RVEs) by using the homogenization theory. However, in current practice, there are some disadvantages of transforming the micro strain domain to the macro stress shakedown domain when considering multiple external loads. The domain cannot fully demonstrate the shakedown condition, and it is impossible to evaluate a known loading combination only from the knowledge of whether the load leads to the shakedown state. To overcome this disadvantage, a new comprehensive approach was proposed to enhance endurance limit strength of composites under variable loads for long term. Considering the example of in-plane strength analysis, for microcosmic RVEs, a new set of boundary condition was defined in the form of uniform strain. The boundary condition was derived from the elastic response under unit loads by using Hook’s law and stiffness matrix. The resulting elastic stress field was used later for plastic shakedown analysis. Based on the lower bound theorem of plastic mechanics, optimization programming for load factor was performed, and after proper mathematical reformulation, the conic quadratic optimization problem could be solved efficiently. Macro-stress shakedown domain can be obtained after scale-transformation of the RVE results. The bases of this stress domain are unidirectional stress in geometry space. The stress amplitude of a structure can be evaluated by this domain for determining the shakedown state in a simple and practical manner. Further, changes in the boundary condition of RVE do not affect the limit and elastic analysis. Finally, few numerical examples were presented for verification and illustration. This approach can be expanded to three 收稿日期: 2019−01−06 基金项目: 国家自然科学基金资助项目（51475036） 工程科学学报，第 41 卷，第 12 期：1558−1566，2019 年 12 月 Chinese Journal of Engineering, Vol. 41, No. 12: 1558−1566, December 2019 DOI：10.13374/j.issn2095-9389.2019.01.06.001; http://journals.ustb.edu.cn
秦方等：基于均匀化理论的复合材料安定性分析方法 .1559· dimensions and employed for more complex structures. KEY WORDS composite material;limit and shakedown analysis;homogenization theory;representative volume element;conic quadratic optimization 非均质复合材料广泛应用于现代工业的方方 换，针对周期性非均质材料的安定性能，通过对 面面，其构成的结构往往长期承受交变载荷的作 RVE边界条件的合理设定，可以直接求出以坐标 用，直至发生塑性破坏.准确计算它们的强度性能 系各向应力为基的宏观安定许用应力域.通过此 在结构设计和安全评估中有重要的作用.基于极 应力域包络线，可以根据交变外载下的应力响应 限安定理论的塑性计算直接法具有不考虑加载路 范围直接判定结构是否满足安定条件，同时，此边 径的优点，通过求解极限及安定域，直接得到结构 界条件的变化也不会影响弹性、塑性极限的结果 破坏的临界值以及可耐久承受的载荷范围.而对 及应用 其分析时，不得不考虑给复合材料带来优异力学 1均匀化理论 特性的材料微观多相物质及其组成的微结构，因 此对于复合材料的研究往往涉及微观及宏观尺度山. 均匀化理论广泛应用于复合材料等效力学性 微观的材料特性及微结构影响着宏观尺度下的复 质的预测.简单讲就是从整体结构中截取一个代 合材料性能表现，两种尺度相关力学性能参数的 表体元，对其施加指定的单位位移边界条件或者 转换可以由均匀化理论方法来解释和处理.Suquet 力边界条件，令其应变能和一个均匀材料的应变 系统的阐释了均匀化思想和计算方法.综合均匀 能等效，从而求出等效性质 化方法和塑性极限安定分析的研究亦有进展.以 若考察宏观非均质连续体V,具有微结构2. 论文[3-4]开始，Weichert等应用下限定理，计算二 通过均匀化理论，将考虑力学问题的两种尺度，宏 维复合材料的安定域.Ponter与Leckie!s则使用上 观和微观以及建立其连接的两种转换方法.通过 限定理，在论文中计算同时考虑温度载荷下的复 缩小尺度或称局部化可以把宏观材料的微结构提 合材料塑性性能 取出来进行分析；而全局化则相反，是将微观代表 经过十几年的发展，Hachemi等I研究了非均 性单元的性质作为宏观层面的平均等效值赋予回 质复合材料的极限状态，综合介绍了使用极限安 宏观层次.宏观参数应变E、应力Σ和微观量具有 定理论对非均质材料进行宏观性能预测的理论发 如下关系： 展、数值计算流程，包含多个结合均匀化方法的算 E=(8〉= bsdv (1) 例.Chen等在此基础之上，即研究了金属基颗粒 1 增强材料(particle-reinforced metal matrix composite,. =()=gloordv (2) PRMMC)性能预测时如何减少代表性单元(RVE) 式中，ε和σ为微观局部的应变和应力.)运算表 尺寸效应的影响可，也从统计学角度对大量极限安 示在RVE区域内的体积平均.当非均质体承受外 定计算结果进行处理以预测材料宏观性能.Zhang 载产生均匀应力应变时，远离边界条件的代表性 等研究了多孔材料在宏观的屈服准则表现.国 体元将跟随几何周期性而产生周期性的应力应 内方面，李华祥等、张宏涛等叫在该领域也有 变.换言之，微观物理场可以分解为平均值和扰动 所研究建树 量两部分： 以上研究计算结果的目的往往是推导材料的 u(x)=E·x+i(x) (3) 宏观屈服强度数值或者等效屈服面，较少直接针 对工程结构的设计和安全校核.本文课题组针对 s(x)=E+(x) (4) 呈现各向异性屈服特性的复合材料也进行了研 σ(x)=∑+G(x) (5) 究，扩展了极限安定计算方法的屈服准则形式， 式中，、(x)和c(x分别是局部速度场、应变和应 适合直接计算材料所构成结构的承载特性.该分 力场的扰动量.这些扰动量具有周期性.根据式 析思路需要材料各方向等效强度参数，往往并不 (1)和(2)可知在RVE区域内： 易获得而具有应用障碍.而本文沿用其优化计算 (〉=0 (6) 方法，更多的关注材料微观和宏观性能的尺度转 (G)=0 (7)
dimensions and employed for more complex structures. KEY WORDS composite material； limit and shakedown analysis； homogenization theory； representative volume element； conic quadratic optimization 非均质复合材料广泛应用于现代工业的方方 面面，其构成的结构往往长期承受交变载荷的作 用，直至发生塑性破坏. 准确计算它们的强度性能 在结构设计和安全评估中有重要的作用. 基于极 限安定理论的塑性计算直接法具有不考虑加载路 径的优点，通过求解极限及安定域，直接得到结构 破坏的临界值以及可耐久承受的载荷范围. 而对 其分析时，不得不考虑给复合材料带来优异力学 特性的材料微观多相物质及其组成的微结构，因 此对于复合材料的研究往往涉及微观及宏观尺度[1] . 微观的材料特性及微结构影响着宏观尺度下的复 合材料性能表现，两种尺度相关力学性能参数的 转换可以由均匀化理论方法来解释和处理. Suquet[2] 系统的阐释了均匀化思想和计算方法. 综合均匀 化方法和塑性极限安定分析的研究亦有进展. 以 论文 [3−4] 开始，Weichert 等应用下限定理，计算二 维复合材料的安定域. Ponter 与 Leckie[5] 则使用上 限定理，在论文中计算同时考虑温度载荷下的复 合材料塑性性能. 经过十几年的发展，Hachemi 等[6] 研究了非均 质复合材料的极限状态，综合介绍了使用极限安 定理论对非均质材料进行宏观性能预测的理论发 展、数值计算流程，包含多个结合均匀化方法的算 例. Chen 等在此基础之上，即研究了金属基颗粒 增强材料（particle-reinforced metal matrix composite, PRMMC）性能预测时如何减少代表性单元（RVE） 尺寸效应的影响[7] ，也从统计学角度对大量极限安 定计算结果进行处理以预测材料宏观性能[8] . Zhang 等[9] 研究了多孔材料在宏观的屈服准则表现. 国 内方面，李华祥等[10]、张宏涛等[11] 在该领域也有 所研究建树. 以上研究计算结果的目的往往是推导材料的 宏观屈服强度数值或者等效屈服面，较少直接针 对工程结构的设计和安全校核. 本文课题组针对 呈现各向异性屈服特性的复合材料也进行了研 究[12] ，扩展了极限安定计算方法的屈服准则形式， 适合直接计算材料所构成结构的承载特性. 该分 析思路需要材料各方向等效强度参数，往往并不 易获得而具有应用障碍. 而本文沿用其优化计算 方法，更多的关注材料微观和宏观性能的尺度转 换，针对周期性非均质材料的安定性能，通过对 RVE 边界条件的合理设定，可以直接求出以坐标 系各向应力为基的宏观安定许用应力域. 通过此 应力域包络线，可以根据交变外载下的应力响应 范围直接判定结构是否满足安定条件，同时，此边 界条件的变化也不会影响弹性、塑性极限的结果 及应用. 1 均匀化理论 均匀化理论广泛应用于复合材料等效力学性 质的预测. 简单讲就是从整体结构中截取一个代 表体元，对其施加指定的单位位移边界条件或者 力边界条件，令其应变能和一个均匀材料的应变 能等效，从而求出等效性质. 若考察宏观非均质连续体 V，具有微结构 Ω. 通过均匀化理论，将考虑力学问题的两种尺度，宏 观和微观以及建立其连接的两种转换方法. 通过 缩小尺度或称局部化可以把宏观材料的微结构提 取出来进行分析；而全局化则相反，是将微观代表 性单元的性质作为宏观层面的平均等效值赋予回 宏观层次. 宏观参数应变 E、应力 Σ 和微观量具有 如下关系： E = ⟨ε⟩ = 1 Ω ∫Ω εdV （1） Σ = ⟨σ⟩ = 1 Ω ∫Ω σdV （2） 式中，ε 和 σ 为微观局部的应变和应力. ⟨·⟩ 运算表 示在 RVE 区域内的体积平均. 当非均质体承受外 载产生均匀应力应变时，远离边界条件的代表性 体元将跟随几何周期性而产生周期性的应力应 变. 换言之，微观物理场可以分解为平均值和扰动 量两部分： u(x) = E · x+u˜ (x) （3） ε(x) = E +ε˜(x) （4） σ(x) = Σ +σ˜ (x) （5） 式中，u˜、ε˜(x) 和σ˜ (x) 分别是局部速度场、应变和应 力场的扰动量. 这些扰动量具有周期性. 根据式 （1）和（2）可知在 RVE 区域内： ⟨ε˜⟩ = 0 （6） ⟨σ˜⟩ = 0 （7） 秦 方等： 基于均匀化理论的复合材料安定性分析方法 · 1559 ·
·1560 工程科学学报，第41卷，第12期 对于分别满足周期性和反周期性的容许速度 →∑cp,=C=0 (16) 场和应力场，结合经典弹塑性理论，对代表性单元 应用虚功原理： 式中，M为单元数目，NGE为每单元高斯点数目， ∑：E=(o:E) (8) NG为全局高斯点数目.通过所有单元矩阵的总 装，得到整体平衡矩阵C.安定下限求解问题最终 此虚功方程被称为Hill-Mandel宏观同质方 程，把材料局部微结构的非均质特性和宏观性质 转换为非线性优化问题： max o 平均量进行了结合，是进一步开展分析工作的基础 st∑Ck=0 2安定下限定理及数值计算优化方法 F(ao+p≤Y (17) 安定下限定理由Melan提出，又被称为静力安 r=1,…,NG;k∈[l,Ww] 定定理.它的物理表述为：如果存在一个自平衡的 式中，σ为弹性外载产生的弹性应力场，y为屈服 残余应力场，与外界载荷产生的弹性应力场叠加 应力数值，v为载荷域角点数目.当结构较为复 后，结构处处不违反屈服准则，则结构满足安定状 杂，单元数目较多时，优化问题求解需要消耗较多 态.下限计算分析就是寻找满足条件的最大的外 的算力，为了提高计算效率，需对问题进行数学等 载荷系数α.和上限定理相比其结果更为保守，适 效转换.下面以von Mises屈服准则，Nv=l为例进 合用于结构安全评估.结合均匀化的条件，其物理 行阐述.结构应力可以分解为残余应力和弹性应 表述为： 力两部分，即 F(aoE+p)<0 (9) r=p+aofk (18) V,p=0 (10) 引入矩阵T∈Rx6,然后以它定义r: p.n=0 (11) 1/21/21/2 -1/21/21/2 p)=0 (12) -1/2-1/21/2 1/w6 (19) 式中，F为屈服准则，p为时间无关残余应力，G为 外载荷产生的弹性应力场，式(10)(12)在RVE单 1/W6 元的体积2内成立，n为体积边界2的外法向单 1/N6 位向量，式(11)在受载边界面厂成立，且pn在的 Vr=T-or (20) 相对表面上具有反周期性 通过这个变换，等式约束条件可变换为： 本文采用纯数学方法，通过数值描述残余应 (21) 力场，使之满足平衡条件和其他约束条件，求解对 ∑c,-a∑ca}=0 应的最优化问题确定其分布数值及对应的载荷系 分解C,Ty,可得： 数.这些工作的第一步是结合有限元方法对物理 CTy,=∑alC,T+c,T (22) 描述进行离散化.下式对自平衡残余应力场进行 处理.此过程中使用了高斯积分法，结合积分点权 其中，[C,T表示[C,T矩阵的第j列，是y,的第 重"，对单元内所有积分点进行计算.式中J为雅 j个张量分量.现在引入新的向量Z,eR以及 各比矩阵，B为单元矩阵，V。为单个单元的体积， x,∈RNG: 5,1,为有限元计算中标准坐标系的三个基 Z ly6s pdv=fuBTpdv=0 (13) (23) -∑Bnv= Z,= (14) ∑666B'oUIagdrd z ÷∑∑，= x,= (24) (15) ∑∑c=0 定义新矩阵A,∈R3M5,B∈R3NN%,以及向 量weR3k及u,∈R5:
对于分别满足周期性和反周期性的容许速度 场和应力场，结合经典弹塑性理论，对代表性单元 应用虚功原理： Σ : E = ⟨σ : ε⟩ （8） 此虚功方程被称为 Hill-Mandel 宏观同质方 程，把材料局部微结构的非均质特性和宏观性质 平均量进行了结合，是进一步开展分析工作的基础. 2 安定下限定理及数值计算优化方法 安定下限定理由 Melan 提出，又被称为静力安 定定理. 它的物理表述为：如果存在一个自平衡的 残余应力场，与外界载荷产生的弹性应力场叠加 后，结构处处不违反屈服准则，则结构满足安定状 态. 下限计算分析就是寻找满足条件的最大的外 载荷系数 α. 和上限定理相比其结果更为保守，适 合用于结构安全评估. 结合均匀化的条件，其物理 表述为： F ( ασ E +ρ ) ⩽ 0 （9） ∇ · ρ = 0 （10） ρ · n = 0 （11） ⟨ρ⟩ = 0 （12） 式中，F 为屈服准则，ρ 为时间无关残余应力，σ E 为 外载荷产生的弹性应力场，式（10）（12）在 RVE 单 元的体积 Ω 内成立，n 为体积边界∂Ω 的外法向单 位向量，式（11）在受载边界面 Γ 成立，且 ρ·n 在的 相对表面上具有反周期性 本文采用纯数学方法，通过数值描述残余应 力场，使之满足平衡条件和其他约束条件，求解对 应的最优化问题确定其分布数值及对应的载荷系 数. 这些工作的第一步是结合有限元方法对物理 描述进行离散化. 下式对自平衡残余应力场进行 处理. 此过程中使用了高斯积分法，结合积分点权 重 w，对单元内所有积分点进行计算. 式中 J 为雅 各比矩阵，B 为单元矩阵，Ve 为单个单元的体积， ξ,η,ζ 为有限元计算中标准坐标系的三个基. ∫V δε T ρdV = ∫V δu TB T ρdV = 0 （13） ⇒ ∑NE 1 ∫Ve B T ρdV = ∑NE 1 ∫ 1 0 ∫ 1 0 ∫ 1 0 B T ρ|J|dξdηdζ （14） ⇒ ∑NE i=1 ∑NGE j=1 wi j|J| i jB T i jρi j = ∑NE i=1 ∑NGE j=1 Ci jρi j = 0 （15） ⇒ ∑NG i=1 Ciρi = [C]{ρ} = 0 （16） 式中，NE 为单元数目，NGE 为每单元高斯点数目， NG 为全局高斯点数目. 通过所有单元矩阵的总 装，得到整体平衡矩阵 C. 安定下限求解问题最终 转换为非线性优化问题： max α s.t. ∑NG r=1 Crρr,k = 0 F ( ασ E r,k +ρr,k ) ⩽ σY,r r = 1,··· ,NG; k ∈ [1,NV] （17） σ E 式中， k 为弹性外载产生的弹性应力场，σY 为屈服 应力数值，NV 为载荷域角点数目. 当结构较为复 杂，单元数目较多时，优化问题求解需要消耗较多 的算力，为了提高计算效率，需对问题进行数学等 效转换. 下面以 von Mises 屈服准则，NV=1 为例进 行阐述. 结构应力可以分解为残余应力和弹性应 力两部分，即 σr = ρ+ασ E r,k （18） T ∈ R 引入矩阵 6×6 ，然后以它定义 νr： T = 1/2 1/2 1/2 −1/2 1/2 1/2 −1/2 −1/2 1/2 1/ √ 6 1/ √ 6 1/ √ 6 （19） νr = T −1σr （20） 通过这个变换，等式约束条件可变换为： ∑NG r=1 CrTνr −α (∑NG r=1 Crσ E r ) = 0 （21） 分解 CrTνr 可得： CrTνr = ∑6 j=1, j,3 [CrT] j ν j r +[CrT] 3 ν 3 r （22） [CrT] j [CrT] ν j r Zr ∈ R 5 xr ∈ R NG 其中， 表示 矩阵的第 j 列 ， 是 νr 的第 j 个张量分量 . 现在引入新的向量 以 及 ： Zr = Z 1 r Z 2 r Z 3 r Z 4 r Z 5 r = v 1 r v 2 r v 4 r v 5 r v 6 r （23） xr = ν 3 r （24） Ar ∈ R 3NK×5 B ∈ R 3NK×NG w ∈ R 3NK ur ∈ R 5 定义新矩阵 ， ，以及向 量 及 ： · 1560 · 工程科学学报，第 41 卷，第 12 期
秦方等：基于均匀化理论的复合材料安定性分析方法 ·1561 A,=[(C,I)1(C,T2(C,I4(C,I5(C,I6] 时普遍采用均布在边界上的位移载荷，这是一种 L-T/V20Yr 均匀边界条件，符合材料实际承载后的微观变形 (25) 特性，同时其加载计算也较简便.本小节将研究微 N=CT3(C2T3…(CT3 (26) 观尺度计算结果在转换为宏观应力场时的一些问 题和改进方法.下面以研究平面应力状态下复合 w=∑co (27) 材料面内性能为例开展.取该材料边长为1的RVE 单元，在微观层面的边界条件为均匀位移约束，即 ur LTZr/V2oYn (28) 互相独立的、4，图1为加载条件示意图 其中的L定义如下： 「2 1/W23/W2 L= (29) 最终等式约束变为： 了G rArur+Nxr-0w=0 (30) 接下来处理优化问题中的不等式约束，将式 (20)代入(17)中的不等式得： 1 2+(}++}+(}+}2+(}2≤2y2 图1RVE边界条件示意图 (31) Fig.1 Boundary condition of RVE 同样，式中严表示应力张量y,的第m个组件 假设RVE的刚度矩阵为K,对于关心的y向 结合式(23)和(28)，不等式约束(31)转换为欧 应力应变关系有如下关系式： 几里得球约束的形式： (x=k118x+k128y (34) zz=,2≤1 (32) σy=k21ex+k22Ey 其中符号为欧几里得范数.综上，对于一个载荷 所研究的面内性能对应了两向独立载荷，即 点的工况，极限求解的优化问题处理为： 载荷域角点Nv=4的工况.结合均匀化理论，载荷 max o 域的顶点在微观位移域和宏观应力域的对应关 st∑A,+N,-aw=0 (33) 系为表1.其中α为弹性或极限安定计算的载荷 Ilul2≤1，r=1,…,NG 系数. 此规划问题有(6NG+1)个变量和(NG+3Nk)个 表1由微观位移载荷城推导宏观应力域 约束条件(NK为结构节点总数目).相比处理变形 Table 1 Derivation of macro-scale stress domain from micro-scale 之前的问题规模得到了明显的缩减.同时，目标函 displacement boundary 数位于等式约束而非先前的不等式约束之中，作 微观位移载荷域 宏观应力域 为圆锥二次优化问题可以得到较高效的求解效率 点 x方向y方向 x方向 y方向 3尺度转换研究 P 0 0 0 0 P2 0 ak114Ⅱ ak4Ⅱ 对具有微观结构特征的复合材料进行分析 aily auy a(kuux+krzuy)/1 a(k2ux+k2zuy)/l 时，往往需要微观宏观两个尺度水平的研究.即在 P 0 ak24,Ⅱ akzzu/l 微观层面对RVE进行承载分析得到微观特性，再 通过尺度转换到宏观层面用以指导设计和校核的 图2展示了微观尺度呈矩形的位移载荷域及 实践.微观分析时RVE的边界条件有均匀边界条 对应到宏观的，表现为平行四边形的应力域 件和周期性边界条件两种设定形式，均可通过特 不同于极限分析和弹性分析，安定分析的载 定的位移或力载荷来实现. 荷域结果为允许载荷的变化范围，是独立载荷组 对于常见的周期性复合材料，进行微观分析 合构成的多维载荷空间.这些独立的载荷可称之
Ar = [ (CrT) 1 (CrT) 2 (CrT) 4 (CrT) 5 (CrT) 6 ] L −T / √ 2σY,r （25） N = [ (C1T) 3 (C2T) 3 ··· ( CNG T )3 ] （26） w = ∑NG i=1 Crσ E r （27） ur = L TZr/ √ 2σY,r （28） 其中的 L 定义如下: L = √ 2 1/ √ 2 √ 3/ √ 2 1 1 1 （29） 最终等式约束变为： ∑NG r=1 Arur + Nxr −αw = 0 （30） 接下来处理优化问题中的不等式约束，将式 （20）代入（17）中的不等式得： ( ν 1 r )2 + ( ν 2 r )2 + ( ν 1 r +ν 2 r )2 + ( ν 4 r )2 + ( ν 5 r )2 + ( ν 6 r )2 ⩽ 2σY,r 2 （31） v m 同样，式中 r 表示应力张量 νr 的第 m 个组件. 结合式（23）和（28），不等式约束（31）转换为欧 几里得球约束的形式： L TZr 2 = ur 2 ⩽ 1 （32） 其中‖∙‖符号为欧几里得范数. 综上，对于一个载荷 点的工况，极限求解的优化问题处理为： max α s.t. ∑NG r=1 Arur + Nxr −αw = 0 ∥ur∥ 2 ⩽ 1,r = 1,··· ,NG （33） 此规划问题有（6NG+1）个变量和（NG+3NK）个 约束条件（NK 为结构节点总数目）. 相比处理变形 之前的问题规模得到了明显的缩减. 同时，目标函 数位于等式约束而非先前的不等式约束之中，作 为圆锥二次优化问题可以得到较高效的求解效率. 3 尺度转换研究 对具有微观结构特征的复合材料进行分析 时，往往需要微观宏观两个尺度水平的研究. 即在 微观层面对 RVE 进行承载分析得到微观特性，再 通过尺度转换到宏观层面用以指导设计和校核的 实践. 微观分析时 RVE 的边界条件有均匀边界条 件和周期性边界条件两种设定形式，均可通过特 定的位移或力载荷来实现. 对于常见的周期性复合材料，进行微观分析 时普遍采用均布在边界上的位移载荷，这是一种 均匀边界条件，符合材料实际承载后的微观变形 特性，同时其加载计算也较简便. 本小节将研究微 观尺度计算结果在转换为宏观应力场时的一些问 题和改进方法. 下面以研究平面应力状态下复合 材料面内性能为例开展. 取该材料边长为 l 的 RVE 单元，在微观层面的边界条件为均匀位移约束，即 互相独立的 ux、uy，图 1 为加载条件示意图. 假设 RVE 的刚度矩阵为 K，对于关心的 xy 向 应力应变关系有如下关系式： { σx = k11εx +k12εy σy = k21εx +k22εy （34） 所研究的面内性能对应了两向独立载荷，即 载荷域角点 NV=4 的工况. 结合均匀化理论，载荷 域的顶点在微观位移域和宏观应力域的对应关 系为表 1. 其中 α 为弹性或极限安定计算的载荷 系数. 图 2 展示了微观尺度呈矩形的位移载荷域及 对应到宏观的，表现为平行四边形的应力域. 不同于极限分析和弹性分析，安定分析的载 荷域结果为允许载荷的变化范围，是独立载荷组 合构成的多维载荷空间. 这些独立的载荷可称之 表 1 由微观位移载荷域推导宏观应力域 Table 1 Derivation of macro-scale stress domain from micro-scale displacement boundary 点 微观位移载荷域 宏观应力域 x方向 y方向 x方向 y方向 P1 0 0 0 0 P2 αux 0 αk11ux /l αk21ux /l P3 αux αuy α ( k11ux +k12uy ) /l α ( k21ux +k22uy ) /l P4 0 αuy αk12uy /l αk22uy /l uy ux l l 图 1 RVE 边界条件示意图 Fig.1 Boundary condition of RVE 秦 方等： 基于均匀化理论的复合材料安定性分析方法 · 1561 ·
1562 工程科学学报，第41卷，第12期 (a) (b) P 图2微观位移载荷域()和对应宏观应力域(b)对应示意图 Fig.2 Displacement load in micro-scale (a)and the corresponding stress domain in macro-scale (b) 为安定载荷域的基，其对应的应力状态为安定容 标系下为矩形形式的容许应力载荷域 许应力域.由图2各载荷角点的对应关系可知，该 f(ux,uy)=sulux+si2luy 宏观应力域虽然绘制在正应力为基的坐标系上， (35) 但它实际表达含义的基仍然是由PP4点的载荷决 g(ux,uy)=s21lux+s22luy 定，而并非沿坐标轴方向的正应力.对于宏观结构 来说，图中的虚线是可以满足安定条件的应力响 1g4,,) 应曲线，一旦加载路径对应的应力超出平行四边 形区域，则不能确定是否满足安定条件 这一特性是计算结果进入工程应用的障碍 在结构设计以及校核检修之中，即便得到受载后 应力的上下限范围，也无法通过此宏观容许应力 f4,4,) 域来判定结构是否满足安定条件.若使用加载路 P 径来判定，就失去了直接法分析的优势.这根本原 图3以人g函数为基的微观位移载荷域 因是宏观容许应力域的基和实际需求不符合.在 Fig.3 Micro-scale displacement load with function f.g as bases 这种情况之下，若能改变宏观容许应力域的基为 4 数值算例 单向应力，则结构安定与否的判定将清晰明了，计 算结果将具有明确的实践指导价值.根据此目的， 4.1算例1 修改微观尺度的位移边界条件，对应关系见表2， 算例1的目的是数值验证论文所采用的极限 表中S为RVE单元柔度矩阵 安定优化算法及求解工具.该程序采用了上文的 安定列式及转换为二次圆锥优化的数学方法，所 表2由宏观应力域推导对应的微观位移载荷域 结合的商业有限元软件为Ansys..有限元软件的作 Table 2 Derivation of micro-scale displacement boundary from macro-scale stress domain 用是计算模型在单位载荷下的纯弹性应力响应 宏观应力域 微观位移载荷域 弹性分析之后，提取有限元模型的单元、节点、材 x方向y方向 x方向 方向 料、边界条件及应力响应数据之后，所有的计算均 Pr 0 0 0 0 在Matlab环境下的以Gurobi为核心求解器的自制 E 0 SuEl 5212l 软件平台中进行 P3 5+52 S21El什s2zl 算例模型为承受两向拉力的圆孔方板模型， 0 古 s122l S2l 它是极限安定相关论文中最常用的对比模型，其 几何材料参数见表3.模型载荷示意图以及安定计 将表2中的位移单位载荷绘制到以f(ux,y)和 算结果见图4. g(4,4,)为基的坐标系中，见图3.此载荷域符合一 图中附加作对比的相关曲线来自文献3均，它 般极限安定计算的需求，单位载荷f(ux,4)和g(山x,4) 们和本文计算结果数值相近，趋势相同.数值误差 的定义见式(34).此单位载荷即为所需的，可以在 可认为是有限元软件单元及网格不同带来的差 尺度转换之后，得到以各向应力为基的，在应力坐 异.结果曲线的吻合从数值上验证优化求解算法
为安定载荷域的基，其对应的应力状态为安定容 许应力域. 由图 2 各载荷角点的对应关系可知，该 宏观应力域虽然绘制在正应力为基的坐标系上， 但它实际表达含义的基仍然是由 P2P4 点的载荷决 定，而并非沿坐标轴方向的正应力. 对于宏观结构 来说，图中的虚线是可以满足安定条件的应力响 应曲线，一旦加载路径对应的应力超出平行四边 形区域，则不能确定是否满足安定条件. 这一特性是计算结果进入工程应用的障碍. 在结构设计以及校核检修之中，即便得到受载后 应力的上下限范围，也无法通过此宏观容许应力 域来判定结构是否满足安定条件. 若使用加载路 径来判定，就失去了直接法分析的优势. 这根本原 因是宏观容许应力域的基和实际需求不符合. 在 这种情况之下，若能改变宏观容许应力域的基为 单向应力，则结构安定与否的判定将清晰明了，计 算结果将具有明确的实践指导价值. 根据此目的， 修改微观尺度的位移边界条件，对应关系见表 2， 表中 S 为 RVE 单元柔度矩阵. f ( ux,uy ) g ( ux,uy ) f ( ux,uy ) g ( ux,uy ) 将表 2 中的位移单位载荷绘制到以 和 为基的坐标系中，见图 3. 此载荷域符合一 般极限安定计算的需求，单位载荷 和 的定义见式（34）. 此单位载荷即为所需的，可以在 尺度转换之后，得到以各向应力为基的，在应力坐 标系下为矩形形式的容许应力载荷域. f ( ux,uy ) = s11lux + s12luy g ( ux,uy ) = s21lux + s22luy （35） 4 数值算例 4.1 算例 1 算例 1 的目的是数值验证论文所采用的极限 安定优化算法及求解工具. 该程序采用了上文的 安定列式及转换为二次圆锥优化的数学方法，所 结合的商业有限元软件为 Ansys. 有限元软件的作 用是计算模型在单位载荷下的纯弹性应力响应. 弹性分析之后，提取有限元模型的单元、节点、材 料、边界条件及应力响应数据之后，所有的计算均 在 Matlab 环境下的以 Gurobi 为核心求解器的自制 软件平台中进行. 算例模型为承受两向拉力的圆孔方板模型， 它是极限安定相关论文中最常用的对比模型，其 几何材料参数见表 3. 模型载荷示意图以及安定计 算结果见图 4. 图中附加作对比的相关曲线来自文献[13−15] ，它 们和本文计算结果数值相近，趋势相同. 数值误差 可认为是有限元软件单元及网格不同带来的差 异. 结果曲线的吻合从数值上验证优化求解算法 表 2 由宏观应力域推导对应的微观位移载荷域 Table 2 Derivation of micro-scale displacement boundary from macro-scale stress domain 点 宏观应力域 微观位移载荷域 x方向 y方向 x方向 y方向 P1 0 0 0 0 P2 Σ1 0 s11Σ1 l s21Σ1 l P3 Σ1 Σ2 s11Σ1 l+s12Σ2 l s21Σ1 l+s22Σ2 l P4 0 Σ2 s12Σ2 l s22Σ2 l (a) uy ux P4 P3 P1 P2 (b) Σy Σx P4 P3 P1 P2 图 2 微观位移载荷域 (a) 和对应宏观应力域 (b) 对应示意图 Fig.2 Displacement load in micro-scale (a) and the corresponding stress domain in macro-scale (b) g(ux , uy ) f(ux , uy ) P4 P3 P1 P2 图 3 以 f、g 函数为基的微观位移载荷域 Fig.3 Micro-scale displacement load with function f, g as bases · 1562 · 工程科学学报，第 41 卷，第 12 期
秦方等：基于均匀化理论的复合材料安定性分析方法 1563· 表3圆孔方板模型几何和材料参数 及软件平台的准确性，算例运算过程也展现了其 Table 3 Geometry and material parameters of plate with hole 优异的计算效率. 长度.mm 孔径，D/mm 厚度，hlmm DI 4.2算例2 100 汤 0.2 算例2的主要目的是展示上文提出的边界设 材料 弹性模量.EMPa泊松比.D 屈服强度，oy/MPa 定方法的求解过程及结果分析.本算例在前述工 作的基础上，结合了本课题组论文的工作，引 钢材 210 0.3 280 入Hill各向异性屈服准则代替von Mises屈服准 0.9 (a) (b) ·一安定域 -米-Chen等1 -◆-…Simon等4 0.6 +Carvelli等s 0.3 0 0.3 0.6 0.9 Play 图4圆孔方板算例.(a)模型：b)安定计算结果 Fig.4 Plate with hole:(a)model;(b)shakedown results 则，扩展了可以计算的复合材料所用材料模型.算 系以及铺层角度描述方式见图5.所使用的材料工 例模型为多层纤维加强层合板，其承载性能与各 程弹性强度参数叨见表4，其中E为弹性模量， 层的铺层形式有关6，具体模型为20层铺层结 G为剪切模量，4为泊松比，T为拉伸强度，S为剪 构，铺层角度为[0°/（±45）3/(90)3]s,关于结构坐标 切强度.每层的材料厚度为0.1308mm (a) (b) 1层 2层 N-1层 纤维 方向 >0 N层 图5层合板模型.()层合结构示意图：(b)纤维角度定义 Fig.5 Lamination model:(a)laminated structure;(b)orientation angle of fiber 表4模型所用面板单层材料力学参数 取一个大小为10mm×10mm×1.308mm的矩形微 Table 4 Material properties of single layer 元作为材料代表性单元进行建模.材料方向坐标 E/GPa E,/GPa Go/GPa Hoy M=T/MPa T,/MPa S/MPa 系为x-0y平面，x轴为材料主方向.为求得RVE 146.9 11.4 6.180.30.41370.6 66.5 133.8 的面内力学特性，初次计算，对RVE的x、y向施 加均布位移作为单位载荷 其他弹性模量及强度参数可以按下式计算 经过计算，得到RVE的弹性及塑性极限安定 E:=Ey,Gx=Gxy.Gxry=0.5Ey/(1+Myz),Sx=Sxy 解见图6(a).图中不同的线型代表了胞元容许的 Sz=S:在建模中，根据此参数，每层结构均视作 位移载荷域，包括弹性解、安定域和极限解.安定 正交各向异性的均质材料.根据对称及周期性，任 域包络线上的点代表了该微元可以承受对应载荷
及软件平台的准确性，算例运算过程也展现了其 优异的计算效率. 4.2 算例 2 算例 2 的主要目的是展示上文提出的边界设 定方法的求解过程及结果分析. 本算例在前述工 作的基础上，结合了本课题组论文[12] 的工作，引 入 Hill 各向异性屈服准则代替 von Mises 屈服准 [ 0 ◦ /(±45◦ )3/(90◦ )3 ] s 则，扩展了可以计算的复合材料所用材料模型. 算 例模型为多层纤维加强层合板，其承载性能与各 层的铺层形式有关[16] ，具体模型为 20 层铺层结 构，铺层角度为 ，关于结构坐标 系以及铺层角度描述方式见图 5. 所使用的材料工 程弹性强度参数[17] 见表 4，其中 E 为弹性模量， G 为剪切模量，μ 为泊松比，T 为拉伸强度，S 为剪 切强度. 每层的材料厚度为 0.1308 mm. Ez = Ey,Gxz = Gxy,Gxy = 0.5Ey/ ( 1+µyz) S xz = S xy, S yz = S xy 其他弹性模量及强度参数可以按下式计算. ， . 在建模中，根据此参数，每层结构均视作 正交各向异性的均质材料. 根据对称及周期性，任 取一个大小为 10 mm×10 mm×1.308 mm 的矩形微 元作为材料代表性单元进行建模. 材料方向坐标 系为 x-0-y 平面，x 轴为材料主方向. 为求得 RVE 的面内力学特性，初次计算，对 RVE 的 x、y 向施 加均布位移作为单位载荷. 经过计算，得到 RVE 的弹性及塑性极限安定 解见图 6（a）. 图中不同的线型代表了胞元容许的 位移载荷域，包括弹性解、安定域和极限解. 安定 域包络线上的点代表了该微元可以承受对应载荷 表 3 圆孔方板模型几何和材料参数 Table 3 Geometry and material parameters of plate with hole 长度，l/mm 孔径，D/mm 厚度，h/mm D/l 100 20 2 0.2 材料 弹性模量，E/MPa 泊松比，υ 屈服强度，σY/MPa 钢材 210 0.3 280 表 4 模型所用面板单层材料力学参数 Table 4 Material properties of single layer Ex /GPa Ey /GPa Gxy/GPa μxy μyz Tx /MPa Ty /MPa Sxy/MPa 146.9 11.4 6.18 0.3 0.4 1370.6 66.5 133.8 Py (a) l h Px D Py/σY (b) 0.9 0.6 0.3 0 0 0.3 0.6 0.9 安定域 Chen等[13] Simon等[14] Carvelli等[15] Px/σY 图 4 圆孔方板算例. (a) 模型; (b) 安定计算结果 Fig.4 Plate with hole: (a) model; (b) shakedown results z (a) N层 N−1层 2层 1层 …x y y x θ θ θ>0 (b) 2 1 纤维 方向 图 5 层合板模型. (a) 层合结构示意图; (b) 纤维角度定义 Fig.5 Lamination model: (a) laminated structure; (b) orientation angle of fiber 秦 方等： 基于均匀化理论的复合材料安定性分析方法 · 1563 ·
.1564 工程科学学报，第41卷，第12期 0.14 900r …弹性解 (b) 安定域 极限解 750 0.09 600 目004 450 -。一·新边界弹性解 L15 -0.05-0.01-0 005 0.1b 0.15 300 ,,.。对强性解 新边界安定域 4./mm ·一对比安定域 150 一·一新边界极限解 -0.06 一南一一对极阴解 400 600 800 -0.11上 /MPa 图6层合结构分析结果.(a)微观容许位移场：b)宏观容许应力场 Fig.6 Analysis results:(a)micro-scale admissible displacement domain;(b)macro-scale admissible stress domain 区域内的任意加载而不失效，极限解则表明了微 观RVE边界条件的单位载荷，重新开展计算，可 元能承受的极限载荷 以得到另外3条以方块线型标记的容许应力场曲 根据尺度转换方法，将微观位移场结果转换 线，此边界条件具体数值见表5.为了展示清晰， 到宏观尺度应力场坐标系中作为对比解，见图6(b) 宏观应力场结果只绘制了第一象限 中以带有三角形线型标记的曲线.同时，改变微 图6(b)中可以看出，两组不同的边界条件得 表5修改后层合板RVE边界条件 Table 5 Modified boundary condition of lamination RVE mm 单位载荷编号 x-0-y平面 x-0-平面 0-:平面 平面x=10 平面)=10 1 对称约束 对称约束 对称约束 4,=0.033 4,=-0.01 2 对称约束 对称约束 对称约束 4,=-0.01 4,=0.0208 到的材料宏观容许应力场中，代表弹性和极限值 4.3算例3 的曲线是完全重合的.这两条曲线上的点代表了 本算例建立了和论文[6]算例1中的相同的 材料特定的应力状态，是材料的固有属性，不会随 模型，即用于管道的双层多孔金属薄壁，结构及 着计算过程中单位载荷的定义变化而变化.它们 RVE示意图见图7.两层等厚材料分别为钢和铝 良好吻合也一定程度在数值上验证了变换过程的 合金，使用完美弹塑性材料模型，材料参数见表6 正确性 RVE模型尺寸为8mm×6mm×4mm,内孔直径为 而代表安定域的曲线则有明显差异，这是由 3 mm. 于它们代表的含义不同.按照新提出的边界条件 根据对照论文，对RVE模型的边界分别施 计算得的安定应力域包络线在坐标系中为单调递 加两方向的均布位移载荷作为单位载荷，所设定 减形式，且横纵坐标在曲线上有一一对应关系，线 的边界条件见表7.其中坐标系定义、模型对称简 上任意一点对坐标轴的投影都会落在此安定域之 化的方式同算例2. 内，符合以坐标轴单位值作为基的安定域的特征 计算结果见图8中三角形标记曲线.不同于 这些表象都不同于标有三角符号的另一条安定域 对照论文中方形标记的安定域曲线，它具有通过 包络线 两向外载荷响应应力直接判定安定与否的能力. 此包络线上的点具有明确的物理及工程意 除此之外，还设定了一个均布力载荷边界条件的 义，其坐标值就代表了材料容许持续加载保持安 对照组.在力边界条件下，RVE模型在加载之后边 定状态的的x向及y向应力最大幅值（最小应力幅 界产生了明显的非对称特性，全局化的结果并不 值为0).若进行层合板结构面内载荷是否满足安 能准确的反应宏观性能.所以尽管它具有和本文 定条件的判定，只需要得到它两向应力上下限幅 所求解曲线相同的物理含义，但在数值上有明显 值范围，即可直接和此曲线进行对比判定.而旧有 较大的偏差，从数值上表明了非均质材料的RVE 边界条件下安定域转换结果不具有此意义， 模型不适宜采用均布力边界条件
区域内的任意加载而不失效，极限解则表明了微 元能承受的极限载荷. 根据尺度转换方法，将微观位移场结果转换 到宏观尺度应力场坐标系中作为对比解，见图 6（b） 中以带有三角形线型标记的曲线. 同时，改变微 观 RVE 边界条件的单位载荷，重新开展计算，可 以得到另外 3 条以方块线型标记的容许应力场曲 线，此边界条件具体数值见表 5. 为了展示清晰， 宏观应力场结果只绘制了第一象限. 图 6（b）中可以看出，两组不同的边界条件得 到的材料宏观容许应力场中，代表弹性和极限值 的曲线是完全重合的. 这两条曲线上的点代表了 材料特定的应力状态，是材料的固有属性，不会随 着计算过程中单位载荷的定义变化而变化. 它们 良好吻合也一定程度在数值上验证了变换过程的 正确性. 而代表安定域的曲线则有明显差异，这是由 于它们代表的含义不同. 按照新提出的边界条件 计算得的安定应力域包络线在坐标系中为单调递 减形式，且横纵坐标在曲线上有一一对应关系，线 上任意一点对坐标轴的投影都会落在此安定域之 内，符合以坐标轴单位值作为基的安定域的特征. 这些表象都不同于标有三角符号的另一条安定域 包络线. 此包络线上的点具有明确的物理及工程意 义，其坐标值就代表了材料容许持续加载保持安 定状态的的 x 向及 y 向应力最大幅值（最小应力幅 值为 0）. 若进行层合板结构面内载荷是否满足安 定条件的判定，只需要得到它两向应力上下限幅 值范围，即可直接和此曲线进行对比判定. 而旧有 边界条件下安定域转换结果不具有此意义. 4.3 算例 3 本算例建立了和论文 [6] 算例 1 中的相同的 模型，即用于管道的双层多孔金属薄壁，结构及 RVE 示意图见图 7. 两层等厚材料分别为钢和铝 合金，使用完美弹塑性材料模型，材料参数见表 6. RVE 模型尺寸为 8 mm×6 mm×4 mm，内孔直径为 3 mm. 根据对照论文[6] ，对 RVE 模型的边界分别施 加两方向的均布位移载荷作为单位载荷，所设定 的边界条件见表 7. 其中坐标系定义、模型对称简 化的方式同算例 2. 计算结果见图 8 中三角形标记曲线. 不同于 对照论文中方形标记的安定域曲线，它具有通过 两向外载荷响应应力直接判定安定与否的能力. 除此之外，还设定了一个均布力载荷边界条件的 对照组. 在力边界条件下，RVE 模型在加载之后边 界产生了明显的非对称特性，全局化的结果并不 能准确的反应宏观性能. 所以尽管它具有和本文 所求解曲线相同的物理含义，但在数值上有明显 较大的偏差，从数值上表明了非均质材料的 RVE 模型不适宜采用均布力边界条件. 表 5 修改后层合板 RVE 边界条件 Table 5 Modified boundary condition of lamination RVE mm 单位载荷编号 x-0-y平面 x-0-z平面 y-0-z平面 平面x=10 平面y=10 1 对称约束 对称约束 对称约束 ux=0.033 uy=−0.01 2 对称约束 对称约束 对称约束 ux=−0.01 uy=0.0208 −0.11 −0.06 −0.01 0.04 0.09 0.14 −0.15 −0.10 −0.05 0 ux /mm 0.05 0.10 0.15 uy /mm 弹性解 安定域 极限解 (a) 0 150 300 450 600 750 900 0 200 400 600 800 Σx /MPa Σy /MPa 新边界弹性解 对比弹性解 新边界安定域 对比安定域 新边界极限解 对比极限解 (b) 图 6 层合结构分析结果. (a) 微观容许位移场; (b) 宏观容许应力场 Fig.6 Analysis results: (a) micro-scale admissible displacement domain; (b) macro-scale admissible stress domain · 1564 · 工程科学学报，第 41 卷，第 12 期
秦方等：基于均匀化理论的复合材料安定性分析方法 ·1565· (a) (b) ⊙ 图7多孔薄壁材料示意图.(a)结构模型：(b)微观模型：(c)RVE模型 Fig.7 Sketch maps of porous material:(a)structure model;(b)microscopic model;(c)RVE model 表6多孔薄壁材料参数 Table 6 Material parameters of porous materials 材料 弹性模量，E/GPa 泊松比，。 屈服强度，MPa 钢材 200 0.3 360 铝合金 72 0.33 100 表7修改后多孔材料RVE边界条件 Table7 Modified boundary condition of porous material RVE mm 单位载荷编号 x-0-y平面 x-0-:平面 0-:平面 =4平面 =3平面 1 对称约束 对称约束 对称约束 4=0.01 4=-0.001998 2 对称约束 对称约束 对称约束 4=-0.003732 4,=0.01 200 接为工程应用所用，包络线上的点代表了满足安 定条件的单向应力上限值.通过多组算例进行了 150 验证和新方法意义的展示说明.未来此方法体系 可扩展到三维问题，为复合材料的应用设计提供 行 更广泛的帮助 50 ·一新边界安定域 ◆一力边界安定域 ·一传统边界安定域 参考文献 50 100 150 200250 [1]Zheng XX,Zheng X T,Gou L H.The research progress on MPa multiscale method for the mechanical analysis of composites.Adv 图8所计算多孔材料宏观容许应力域结果对比 Mech,2010,40(141 Fig.8 Comparison of macro-scale admissible stress domain of porous (郑晓霞，郑锡涛，猴林虎.多尺度方法在复合材料力学分析中 material 的研究进展.力学进展，2010,40(1)：41) [2] Suquet P M.Discontinuities and plasticity//Nonsmooth Mechanics 5结论 and Applications.Vienna:Springer,1988:279 以数值方法研究了周期性非均质复合材料的 [3]Weichert D,Hachemi A,Schwabe F.Application of shakedown analysis to the plastic design of composites.Arch Appl Mech, 安定域计算及工程安全评估应用问题.通过均匀 1999,69(9-10):623 化理论方法建立微观结构和宏观性能的关系，结 [4]Weichert D,Hachemi A,Schwabe F.Shakedown analysis of 合极限安定下限理论对复合材料代表性单元进行 composites.Mech Res Commun,1999,26(3):309 弹塑性分析.分析了现有尺度转换方法存在的问 [5]Ponter A R S,Leckie F A.On the behaviour of metal matrix 题，即无法得到以单向应力为基的安定容许应力 composites subjected to cyclic thermal loading.J Mech Plnys 域.根据微观位移域和宏观应力域的对应关系，给 Solids,1998,46(11):2183 [6] Hachemi A,Chen M,Chen G,et al.Limit state of structures made 出一种方法利用胡克定律计算定义出一组新的微 of heterogeneous materials.Int/Plast,2014,63:124 观边界条件，使之全局化后可以得到以宏观坐标 [] Chen G,Bezold A,Broeckmann C,et al.On the size of the 轴方向应力为基的宏观安定应力域.该域可以直 representative volume element used for the strength prediction:A
5 结论 以数值方法研究了周期性非均质复合材料的 安定域计算及工程安全评估应用问题. 通过均匀 化理论方法建立微观结构和宏观性能的关系，结 合极限安定下限理论对复合材料代表性单元进行 弹塑性分析. 分析了现有尺度转换方法存在的问 题，即无法得到以单向应力为基的安定容许应力 域. 根据微观位移域和宏观应力域的对应关系，给 出一种方法利用胡克定律计算定义出一组新的微 观边界条件，使之全局化后可以得到以宏观坐标 轴方向应力为基的宏观安定应力域. 该域可以直 接为工程应用所用，包络线上的点代表了满足安 定条件的单向应力上限值. 通过多组算例进行了 验证和新方法意义的展示说明. 未来此方法体系 可扩展到三维问题，为复合材料的应用设计提供 更广泛的帮助. 参 考 文 献 Zheng X X, Zheng X T, Gou L H. The research progress on multiscale method for the mechanical analysis of composites. Adv Mech, 2010, 40（1）: 41 （郑晓霞, 郑锡涛, 缑林虎. 多尺度方法在复合材料力学分析中 的研究进展. 力学进展, 2010, 40（1）：41 ） [1] Suquet P M. Discontinuities and plasticity//Nonsmooth Mechanics and Applications. Vienna: Springer, 1988: 279 [2] Weichert D, Hachemi A, Schwabe F. Application of shakedown analysis to the plastic design of composites. Arch Appl Mech, 1999, 69（9-10）: 623 [3] Weichert D, Hachemi A, Schwabe F. Shakedown analysis of composites. Mech Res Commun, 1999, 26（3）: 309 [4] Ponter A R S, Leckie F A. On the behaviour of metal matrix composites subjected to cyclic thermal loading. J Mech Phys Solids, 1998, 46（11）: 2183 [5] Hachemi A, Chen M, Chen G, et al. Limit state of structures made of heterogeneous materials. Int J Plast, 2014, 63: 124 [6] Chen G, Bezold A, Broeckmann C, et al. On the size of the representative volume element used for the strength prediction: A [7] 表 6 多孔薄壁材料参数 Table 6 Material parameters of porous materials 材料 弹性模量，E/GPa 泊松比，υ 屈服强度，σY/MPa 钢材 200 0.3 360 铝合金 72 0.33 100 表 7 修改后多孔材料 RVE 边界条件 Table 7 Modified boundary condition of porous material RVE mm 单位载荷编号 x-0-y平面 x-0-z平面 y-0-z平面 x=4平面 y=3平面 1 对称约束 对称约束 对称约束 ux=0.01 uy=−0.001998 2 对称约束 对称约束 对称约束 ux=−0.003732 uy=0.01 (a) (b) (c) 图 7 多孔薄壁材料示意图. (a) 结构模型；(b) 微观模型；(c)RVE 模型 Fig.7 Sketch maps of porous material: (a) structure model; (b) microscopic model; (c) RVE model Σx /MPa Σy /MPa 0 50 100 150 200 0 50 100 150 200 250 新边界安定域 力边界安定域 传统边界安定域 图 8 所计算多孔材料宏观容许应力域结果对比 Fig.8 Comparison of macro-scale admissible stress domain of porous material 秦 方等： 基于均匀化理论的复合材料安定性分析方法 · 1565 ·
·1566 工程科学学报，第41卷，第12期 statistical survey applied to the particulate reinforce metal matrix (秦方，张乐乐，陈敏，等.正交各向异性材料塑性极限与安定的 composites(PRMMCsV/Advances in Direct Methods for Materials 下限分析.清华大学学报：自然科学版，2018,58(11)：966) and Structures.Cham:Springer,2018:51 [13]Chen S S,Liu Y H,Cen ZZ.Lower bound shakedown analysis by [8]Chen G,Bezold A,Broeckmann C,et al.On the statistical using the element free Galerkin method and non-linear determination of strength of random heterogeneous materials. programming.Comput Methods Appl Mech Eng,2008,197(45. Compos Struci1,2016,149:220 48):3911 [9] Zhang J,Oueslati A,Shen W Q,et al.Shakedown of porous [14]Simon J W,Weichert D.Numerical lower bound shakedown material with Drucker-Prager dilatant matrix under general cyclic analysis of engineering structures.Comput Methods Appl Mech loadings.Composite Structures,2019,220:566 Eg,2011,200(41-44:2828 [10]Li H X,Liu Y H,Feng X Q,et al.Plastic limit analysis of ductile [15]Carvelli V,Cen Z Z,Liu Y,et al.Shakedown analysis of defective composites based on homogenization theory.Acta Mech Sin,2002. pressure vessels by a kinematic approach.Arch Appl Mech,1999, 34(4):550 69(9-10):751 (李华祥，刘应华，冯西桥，等基于均匀化理论韧性复合材料塑 [16]Huang C F,Xiao J Y,Huang Z H,et al.Buckling of composite 性极限分析.力学学报，2002,34(4)：550) cylindrical shells fabricated using thin-ply under axial [11]Zhang H T,Liu Y H,Xu B Y.Lower bound shakedown analysis compression.Chin J Eng,2018,40(07):857 of periodic ductile composites.J Tsinghua Univ Sci Technol,2005. (黄春芳，肖加余，黄展鸿，等.薄铺层复合材料薄壁管轴压屈曲 45(2):267 行为研究.工程科学学报，2018,40(07)：857) (张宏涛，刘应华，徐秉业.周期性韧性复合材料的安定下限分 [17]Wang YQ,Tong M B,Zhu S H.3D nonlinear progressive damage 析.清华大学学报：自然科学版，2005,45(2)：267) analysis model for composite laminates.Acta Mater Compos Sin, [12]Qin F,Zhang LL,Chen M,et al.Lower bound analysis of plastic 2009,26(5:159 limit and shakedown state of orthotropic materials.J Tsinghua (王跃全，童明波，朱书华.三维复合材料层合板渐进损伤非线 Uniy Sci Technol,2018,58(11):966 性分析模型.复合材料学报，2009,26(5)：159)
statistical survey applied to the particulate reinforce metal matrix composites (PRMMCs)//Advances in Direct Methods for Materials and Structures. Cham: Springer, 2018: 51 Chen G, Bezold A, Broeckmann C, et al. On the statistical determination of strength of random heterogeneous materials. Compos Struct, 2016, 149: 220 [8] Zhang J, Oueslati A, Shen W Q, et al. Shakedown of porous material with Drucker-Prager dilatant matrix under general cyclic loadings. Composite Structures, 2019, 220: 566 [9] Li H X, Liu Y H, Feng X Q, et al. Plastic limit analysis of ductile composites based on homogenization theory. Acta Mech Sin, 2002, 34（4）: 550 （李华祥, 刘应华, 冯西桥, 等. 基于均匀化理论韧性复合材料塑 性极限分析. 力学学报, 2002, 34（4）：550 ） [10] Zhang H T, Liu Y H, Xu B Y. Lower bound shakedown analysis of periodic ductile composites. J Tsinghua Univ Sci Technol, 2005, 45（2）: 267 （张宏涛, 刘应华, 徐秉业. 周期性韧性复合材料的安定下限分 析. 清华大学学报: 自然科学版, 2005, 45（2）：267 ） [11] Qin F, Zhang L L, Chen M, et al. Lower bound analysis of plastic limit and shakedown state of orthotropic materials. J Tsinghua Univ Sci Technol, 2018, 58（11）: 966 [12] （秦方, 张乐乐, 陈敏, 等. 正交各向异性材料塑性极限与安定的 下限分析. 清华大学学报: 自然科学版, 2018, 58（11）：966 ） Chen S S, Liu Y H, Cen Z Z. Lower bound shakedown analysis by using the element free Galerkin method and non-linear programming. Comput Methods Appl Mech Eng, 2008, 197（45- 48）: 3911 [13] Simon J W, Weichert D. Numerical lower bound shakedown analysis of engineering structures. Comput Methods Appl Mech Eng, 2011, 200（41-44）: 2828 [14] Carvelli V, Cen Z Z, Liu Y, et al. Shakedown analysis of defective pressure vessels by a kinematic approach. Arch Appl Mech, 1999, 69（9-10）: 751 [15] Huang C F, Xiao J Y, Huang Z H, et al. Buckling of composite cylindrical shells fabricated using thin-ply under axial compression. Chin J Eng, 2018, 40（07）: 857 （黄春芳, 肖加余, 黄展鸿, 等. 薄铺层复合材料薄壁管轴压屈曲 行为研究. 工程科学学报, 2018, 40（07）：857 ） [16] Wang Y Q, Tong M B, Zhu S H. 3D nonlinear progressive damage analysis model for composite laminates. Acta Mater Compos Sin, 2009, 26（5）: 159 （王跃全, 童明波, 朱书华. 三维复合材料层合板渐进损伤非线 性分析模型. 复合材料学报, 2009, 26（5）：159 ） [17] · 1566 · 工程科学学报，第 41 卷，第 12 期