·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 期