工程科学学报.第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 期