软化材料的弹塑性有限元分析
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
［ ３ ］
１ 弹塑性本构积分
塑性屈服面在应力空间中可记为 Ｆ （ ， ｋ ） ＝０ σ 考虑等向强化（ 软化） 的情况， 此时式（ １ ） 可记为 Ｆ （ ， ｋ ） ＝ｆ （ ）－ｋ＝０ σ σ （ ２ ） ｆ （ ） 为有效应力； ｋ 为强化（ 软化） 参数， 强化时 ｋ 增 式中： σ 加， 软化时 ｋ 减少。通常 ｋ 是塑性内变量如塑性功、 等效塑 性应变或塑性体应变的函数。由于这 ３者均与塑性应变有
Ι Π ｕ ｕ ｕ δ δ ｉ ＝δ ｉ ＋λ ｉ ｉ － １ Ｔ Ι ｕ Ｂ ｄ ｖ δ σ ｉ ＝ －Ｋ ｉ － １ ｉ － １ － １ Π ｕ δ ｉ ＝Ｋ ｉ － １ ｑ
式（ ７ ） 求出 Ｈ ′ ， 并由式（ ８ ） 求出： Ｆ ： ｄ ′ σ ′ σ Δ λ＝ Ｆ ： Ｄ ： Ｆ ′ ′ ′ σ σ ＋Ｈ ′ －Δ Ｄ ： Ｆ σ ＝σ λ ′ σ ′ Ｆ ε λ Ｐ ＝ε Ｐ ＋Δ ′ σ ｋ＝ｋ （ ） ε ｐ ｋ ， ′ 转入 ５ ） ， 直到 Ｆ （ ， ｋ ）＜ Ｔ Ｏ Ｌ 。 ε σ ｐ＝ ε ｐ ， ２ 配合 Ｔ ｈ ｏ ｍ ａ ｓ 加速的初应力平衡迭代 如前所述， 在结构未达到极限荷载前， 迭代是通过初 应力法配合 Ｔ ｈ ｏ ｍ ａ ｓ 加速来控制的。关于初应力法本文不 再详述， 下面将给出 Ｔ ｈ ｏ ｍ ａ ｓ 加速算法。 带加速因子的初应力迭代可表示为 Δｕ ＝ Ｋ Ｒ
∫
（ ３ １ ）
２ ５ ） 可改写为 这样式（ Ｋ ｕ δ ｉ － １ ｉ ＝－ｇ ｉ － １ 当系统达到平衡时， 系统总势能最小， 即 （ ）＝ φ ρ Δｗ Ｔ ＝δ ｕ （ ｕ ） Δｕ δ ｉｇ ｉ － １ ＋ρ ｉ ρ （ ３ ３ ） 图２ 强化（ 软化） 参数 ｋ 与塑性体应变 θ 满足分段线性关系
１ ］ 得到了工程界的广泛应用 ［ ； 相比较而言， 由于显示法较
每个增量步中， 能通过弧长约束条件使迭代在位移—荷载 空间沿一弧形路径收敛至平衡位置， 且该平衡位置位于上 一增量步平衡位置附近， 这样就可稳定地越过极值点， 从 而避免计算发散的危险。 ｕ ｌ ｅ ｒ 针对上述问题， 本文中根据计算经验选用了显式 Ｅ 法作为本构积分方案。在结构未达到极限荷载前， 采用初 应力法配合 Ｔ ｈ ｏ ｍ ａ ｓ 加速方案作为平衡迭代的方法； 而当结 构接近极限荷载时， 当迭代转换为弧长法控制， 越过极值 点进入软化区直至破坏。
直接， 也易于实现， 并且还可通过误差控制的自适应本构 积分 来 提 高 计 算 精 度 和 效 率， 故近年也得到了新的
２ ］ 发展 ［ 。
平衡迭代是弹塑性非线性有限元分析的又一重要问 题。平衡迭代就是将不平衡力施加于系统得出位移增量， 重复这个过程直至系统离平衡状态的差距可以忽略为止。 根据所采用线性方程组刚度矩阵的不同， 平衡迭代法可分 ｅ ｗ ｔ ｏ ｎ Ｒ ａ ｐ ｈ ｓ ｏ ｎ 法、 拟Ｎ ｅ ｗ ｔ ｏ ｎ － Ｒ ａ ｐ ｈ ｓ ｏ ｎ 法和初应力法。 为Ｎ 关于各种方法及其比较有很多文献可供参考， 这里不再一 一列举。由于初应力法具有极强的稳定性和广泛的适用 性， 且对于非线性较强的行为（ 如应变局部化等） 更具有优 势，因此该方法一直在弹塑性有限元分析中占据着重要地 位。尽管如此， 该法却有一个明显缺点， 那就是收敛速度 缓慢。为了改进初应力法的收敛特性， 研究者已经给出了 Ｓ ｌ ｏ ａ ｎ 等也对此作了详细评述。 很多种加速方案， 此外， 对于进入软化阶段的平衡迭代还需解决 ２个问
ｐ １ ）给出初始应力 σ ， 初始塑性应变 ε ０ ０ 和本荷载步应
（ ２ １ ）
（ ７ ）
这样式（ １ ７ ） 可表示为
ｉ ｉ － １ ｉ － １ ｉ ｉ ｕ ＝ｕ ＋ａ Δｕ 珘 ｅ ＋Δ ｕ ｅ
源自文库
（ ２ ２ ）
（ ８ ）
令 ａ＝ １ ， 利用式（ ２ １ ） 、 式（ ２ ２ ） 便可完成迭代。 关于 Ｔ ｈ ｏ ｍ ａ ｓ 加速算法还有其它形式， 根据经验， 本文 中的这种形式对于结构内塑性区很大的情况具有很强的 适应性， 并且这种形式更符合传统的有限元格式。
Ι ｕ Δｌ －Δｕ Δｕ ｉ － １ （ ｉ － １ ＋δ ｉ） Ｔ Π ｕ Δｕ ｉ － １ δ ｉ ２ Ｔ
（ ２ ９ ）
式（ ２ ９ ） 即为弧长法荷载增量计算公式。根据不同的简化 可根椐需要选用。 条件， λ ｉ的表示方法也有好几种， 线性搜索在进入软化的有限元分析中十分重要， 将式
（ １ ６ ） （ １ ７ ）
ｐ 为塑性应变 ε 的函数 关， 本文中记 ｋ
（ １ ）
随着塑性变形的发生， 塑性屈服面随之改动。这里仅
收稿日期： ２ ０ １ ０－ ０ １－ ０ ７ 作者简介： 夏怀孝（ １ ９ ７ ４ —） ， 男， 博士， 讲师， 主要从事自适应有限元研究。
夏怀孝： 软化材料的弹塑性有限元分析 １ ３ １
ｐ
（ ３ ２ ）
式中： ｗ为系统总势能。 φ （ ） 为一高度非线性且不可显 Δ ρ 示的函数， 一般利用线性搜索确定 ρ 的大致范围， 即： （ ０ ） φ ρ＝－ （ １ ）－φ （ ０ ） φ
Ｔ Ｔ （ ０ ） ＝δ ｕ （ ） ＝ －δ ｕ ｕ φ Δｕ δ ｉｇ ｉ － １ ｉＫ ｉ － １ ｉ
ｐ ５ ）令 σ ′＝ σ ， ｋ ′＝ ｋ ， ′ ｄ ′＝ σ ， 由 ε σ σ ｅ ０ ｐ＝ ε ０， ｅ－ ０
３ 弧长法控制的迭代
在非线性有限元分析中， 荷载是通过若干小的增量步
{
＜０ 卸载 ９ ） （
进行施加的。在每个小的增量步中， 通过若干次迭代达到 为某增量步的第 ｉ 次迭代， ｕ 平衡。设 ｉ δ ｉ为该迭代步位移 增量， 则该荷载增量步第 ｉ 次迭代的位移可记为 （ １ ０ ） ｕ ｕ ｕ Δ ｉ ＝Δ ｉ － １ ＋δ ｉ 弧长法的约束条件记为 Δｕ Δｕ ｉ － １ ｉ ＝ Δｌ ｌ 保持不变， ｕ Δ δ ｉ可表示为
Ｐ ： （ Δ σ ＝Ｄ Δ ε－Δ ε）
Ｔ ｈ ｏ ｍ ａ ｓ 假定
ｉ ｉ ｉ － １ ｉ ｉ ａ Δｕ 珘 ｅ ≈ ａ Δｕ ｅ ＋Δ ｕ ｅ
（ １ ８ ） １ ９ ） （ （ ２ ０ ）
（ ４ ）
其中
ｉ － １ ｉ － １ ｉ － １ ｉ } {Ｒ ｕ Δ 珘 （ ｕ ＋ａ Δｕ ｅ ｅ ＝Ｋ ｅ）
ｉ ｉ － １ ｉ ｉ ｅ － １ ｉ － １ ｅ ｉ ｅ
（ １ ２ ） （ １ ３ ） （ １ ４ ） （ １ ５ ）
（ ２ ６ ） （ ２ ７ ） （ ２ ８ ）
∫
如Ｆ （ ， ｋ ）＜ Ｔ Ｏ Ｌ ， 转入下一荷载步； 否则令 σ ′＝ σ ， ｋ ′＝ σ
将式（ ２ ６ ） 代入式（ ２ ４ ） 可得 λ ｉ＝
（ ３ ４ ） （ ３ ５ ） （ ３ ６ ）
对于严格的凸问题， 线性搜索要求 （ ０ ）≤ ０ ， φ （ １ ） ＞φ （ ０ ） φ ３ ５ ） 、 式（ ３ ６ ） 可知， 当Ｋ 为非正定（ 结构进入软化） ， 由式（ ｉ － １
第３ １卷 第 ３期 四 川 兵 工 学 报 ２ ０ １ ０年 ３月 【 其他研究】
软化材料的弹塑性有限元分析
夏怀孝
（ 三峡大学 土木水电学院， 湖北 宜昌 ４ ４ ３ ０ ０ ２ ）
摘要： 给出了软化材料弹塑性有限元分析的基本算法：选用显式 Ｅ ｕ ｌ ｅ ｒ 法作为本构积分方案；在结构未达到极限 荷载前采用初应力法配合 Ｔ ｈ ｏ ｍ ａ ｓ 加速方案作为平衡迭代方法； 而当结构接近极限荷载时， 程序转换为弧长法控 制迭代， 使结构越过极值点进入软化区直至破坏。最后， 讨论了该领域面临的重要问题— — —由于材料软化而导 致的应变局部化现象及其有限元模拟。 关键词： 有限元分析； 弹塑性； 软化； 应变局部化 中图分类号： Ｏ ３ ４ ４ ． ３ 文献标识码： Ａ 文章编号： １ ０ ０ ６－ ０ ７ ０ ７ （ ２ ０ １ ０ ） ０ ３－ ０ １ ３ ０－ ０ ３ ｌ ｉ ｍ ｉ ｔ 题， 即如何求 得 极 限 荷 载， 以及如何越过该极值点（ ｐ ｏ ｉ ｎ ｔ ） 使受力—位移过程得以继续发展。这 ２个问题一直 是非线性有限元分析的难点， 研究者也提出了许多方法进 行解决， 如虚拟弹簧法、 位移控制法、 强制迭代法等， 关于
ｐ ｋ＝ｋ （ ε）
（ ３ ）
ｉ － １ ｉ Ｒ 为上一迭代末的不平衡力； ａ 为加速因子。 式中：
假定塑 性 势 与 塑 性 屈 服 面 重 合， 那么根据正交流动法 则， 有 Ｆ ｐ Δ ε ＝Δ λ σ 假设塑性屈服面连续变化， 对式（ １ ） 进行微分可得 Ｆ ｋ ： Ｆ Δ σ ＋Ｆ λ ＝０ ｐ ｋ σ： ε σΔ 由增量弹塑性应力应变关系可得
０
变增量 Δ ， 根据式（ ３ ） 求出 ｋ 。 ε ０ ２ ）判断该荷载步为加载或卸载。由于应力空间中加 载、 卸载准则仅适用于强化阶段， 故采用应变空间中加载 卸载准则： Ｆ ： ｄ ε ＝ ＝０ 中性变载 ε ＞０ 加载 Ｆ Ｆ Ｆ ： ｄ Ｄ ： ｄ ε＝ ： ε＝ ： Δ σ ｅ ε σ σ σ σ ｅ ＝σ ０ ＋Δ ｅ ３ ）若材料处于卸载， ， 转入下一荷载步。 σ＝ σ ｅ ４ ）若材料处于加载， 如Ｆ （ ， ｋ ）＜ ０ ， 则 σ＝ ， 转入 σ σ ｅ ０ ｅ 下一荷载步。
－ １ Ｔ ｕ ｑ－ Ｂ ｄ ｖ ） δ λ σ ｉ ＝Ｋ ｉ － １ （ ｉ ｉ － １ Ｔ ２
（ ２ ３ ） （ ２ ４ ）
式中： Δ σ ｅ为弹性应力增量。相应定义弹性试探应力为 （ １ １ ）
其中： ｌ 为每一增量步的弧长。在该增量步的所有迭代中 Δ
∫
（ ２ ５ ）
其中： Ｋ 为该迭代步的刚度矩阵； ｑ 是单位荷载增量； λ ｉ － １ ｉ 为本迭代步荷载。在弧长法中， ｕ δ ｉ的表示形式有好几种， 这里的表示方法可以较好地与后文将提到的线性搜索相 配合。 ｕ 部分： 将δ ｉ分为 ２
令 （ ５ ） （ ６ ）
ｉ ｉ － １ ａ ＝ａ ＋ ｉ ｉ ｉ ｉ － １ ｉ ε ＝ ａΔｕ 珘 ｅ －ａ Δ ｕ ｅ －Δ ｕ ｅ
利用最小二乘法求出 ａ， 使得 ε 最小
ｉ { Δｕｉ }Ｔ{ Δｕ 珘} ｅ ｅ ｅ ｅ ｉ { Δｕｉ }Ｔ{ Δｕ }
ｉ
定义塑性模量 Ｈ为 Ｈ ＝－Ｆ ｋ ： Ｆ ｐ ｋ ε σ 则由式（ ４ ）～ （ ６ ） 可得 Ｆ Ｄ ： Δ ε σ： Δ λ＝ Ｆ ： Ｄ ： Ｆ ＋Ｈ σ σ ｕ ｌ ｅ ｒ 本构积分的算法。 下面具体给出 Ｅ
［ ４ ］ ｅ ｒ ｇ ａ ｎ 也有详细评述。众多研究表明， Ｒ ｉ ｓ ｋ ｓ 和 这个问题 Ｂ ５ － ６ ］ Ｃ ｒ ｉ ｓ ｆ ｉ ｅ ｌ ｄ 等［ 提出的弧长法是较为有效的方法。该法在
近 ３ ０年来， 弹塑性有限元分析一直是国际数值计算领 域研究的重点。随着工程规模的扩大和相应研究的深入， 结构的极限荷载及其越过极限荷载后的破坏情况等逐渐 引起了工程界的高度重视。因此， 含软化阶段的弹塑性有 限元分析也随之被引入。 本构积分是弹塑性分析最重要的环节之一， 它对计算 结果的精度和计算的收敛性均有着重要影响。弹塑性本 构积分方案大致有显式法和隐式法 ２种。隐式法由于能给 出一致切线刚度矩阵（ ｃ ｏ ｎ ｓ ｉ ｓ ｔ ｅ ｎ ｔ ｔ ａ ｎ ｇ ｅ ｎ ｔ ｓ ｔ ｉ ｆ ｆ ｎ ｅ ｓ ｓ ｍ ａ ｔ ｒ ｉ ｘ ） ， 且 对于 Ｎ ｅ ｗ ｔ ｏ ｎ Ｒ ａ ｐ ｈ ｓ ｏ ｎ法可达到二阶收敛速度， 因而该方法
（ ２ ３ ） 改写为 ｕ ｕ ｕ Δ δ ｉ ＝Δ ｉ － １ ＋ρ ｉ （ ３ ０ ）
ｕ ＝ ｕ ＋ａΔｕ
四川兵工学报 １ ３ ２
式中： 为加速因子。记不平衡力为 ｇ ， 有 ρ
Ｔ ｇ ｄ ｖ －λ ｑ ｉ － １ ＝ Ｂ σ ｉ － １ ｉ
１ 弹塑性本构积分
塑性屈服面在应力空间中可记为 Ｆ （ ， ｋ ） ＝０ σ 考虑等向强化（ 软化） 的情况， 此时式（ １ ） 可记为 Ｆ （ ， ｋ ） ＝ｆ （ ）－ｋ＝０ σ σ （ ２ ） ｆ （ ） 为有效应力； ｋ 为强化（ 软化） 参数， 强化时 ｋ 增 式中： σ 加， 软化时 ｋ 减少。通常 ｋ 是塑性内变量如塑性功、 等效塑 性应变或塑性体应变的函数。由于这 ３者均与塑性应变有
Ι Π ｕ ｕ ｕ δ δ ｉ ＝δ ｉ ＋λ ｉ ｉ － １ Ｔ Ι ｕ Ｂ ｄ ｖ δ σ ｉ ＝ －Ｋ ｉ － １ ｉ － １ － １ Π ｕ δ ｉ ＝Ｋ ｉ － １ ｑ
式（ ７ ） 求出 Ｈ ′ ， 并由式（ ８ ） 求出： Ｆ ： ｄ ′ σ ′ σ Δ λ＝ Ｆ ： Ｄ ： Ｆ ′ ′ ′ σ σ ＋Ｈ ′ －Δ Ｄ ： Ｆ σ ＝σ λ ′ σ ′ Ｆ ε λ Ｐ ＝ε Ｐ ＋Δ ′ σ ｋ＝ｋ （ ） ε ｐ ｋ ， ′ 转入 ５ ） ， 直到 Ｆ （ ， ｋ ）＜ Ｔ Ｏ Ｌ 。 ε σ ｐ＝ ε ｐ ， ２ 配合 Ｔ ｈ ｏ ｍ ａ ｓ 加速的初应力平衡迭代 如前所述， 在结构未达到极限荷载前， 迭代是通过初 应力法配合 Ｔ ｈ ｏ ｍ ａ ｓ 加速来控制的。关于初应力法本文不 再详述， 下面将给出 Ｔ ｈ ｏ ｍ ａ ｓ 加速算法。 带加速因子的初应力迭代可表示为 Δｕ ＝ Ｋ Ｒ
∫
（ ３ １ ）
２ ５ ） 可改写为 这样式（ Ｋ ｕ δ ｉ － １ ｉ ＝－ｇ ｉ － １ 当系统达到平衡时， 系统总势能最小， 即 （ ）＝ φ ρ Δｗ Ｔ ＝δ ｕ （ ｕ ） Δｕ δ ｉｇ ｉ － １ ＋ρ ｉ ρ （ ３ ３ ） 图２ 强化（ 软化） 参数 ｋ 与塑性体应变 θ 满足分段线性关系
１ ］ 得到了工程界的广泛应用 ［ ； 相比较而言， 由于显示法较
每个增量步中， 能通过弧长约束条件使迭代在位移—荷载 空间沿一弧形路径收敛至平衡位置， 且该平衡位置位于上 一增量步平衡位置附近， 这样就可稳定地越过极值点， 从 而避免计算发散的危险。 ｕ ｌ ｅ ｒ 针对上述问题， 本文中根据计算经验选用了显式 Ｅ 法作为本构积分方案。在结构未达到极限荷载前， 采用初 应力法配合 Ｔ ｈ ｏ ｍ ａ ｓ 加速方案作为平衡迭代的方法； 而当结 构接近极限荷载时， 当迭代转换为弧长法控制， 越过极值 点进入软化区直至破坏。
直接， 也易于实现， 并且还可通过误差控制的自适应本构 积分 来 提 高 计 算 精 度 和 效 率， 故近年也得到了新的
２ ］ 发展 ［ 。
平衡迭代是弹塑性非线性有限元分析的又一重要问 题。平衡迭代就是将不平衡力施加于系统得出位移增量， 重复这个过程直至系统离平衡状态的差距可以忽略为止。 根据所采用线性方程组刚度矩阵的不同， 平衡迭代法可分 ｅ ｗ ｔ ｏ ｎ Ｒ ａ ｐ ｈ ｓ ｏ ｎ 法、 拟Ｎ ｅ ｗ ｔ ｏ ｎ － Ｒ ａ ｐ ｈ ｓ ｏ ｎ 法和初应力法。 为Ｎ 关于各种方法及其比较有很多文献可供参考， 这里不再一 一列举。由于初应力法具有极强的稳定性和广泛的适用 性， 且对于非线性较强的行为（ 如应变局部化等） 更具有优 势，因此该方法一直在弹塑性有限元分析中占据着重要地 位。尽管如此， 该法却有一个明显缺点， 那就是收敛速度 缓慢。为了改进初应力法的收敛特性， 研究者已经给出了 Ｓ ｌ ｏ ａ ｎ 等也对此作了详细评述。 很多种加速方案， 此外， 对于进入软化阶段的平衡迭代还需解决 ２个问
ｐ １ ）给出初始应力 σ ， 初始塑性应变 ε ０ ０ 和本荷载步应
（ ２ １ ）
（ ７ ）
这样式（ １ ７ ） 可表示为
ｉ ｉ － １ ｉ － １ ｉ ｉ ｕ ＝ｕ ＋ａ Δｕ 珘 ｅ ＋Δ ｕ ｅ
源自文库
（ ２ ２ ）
（ ８ ）
令 ａ＝ １ ， 利用式（ ２ １ ） 、 式（ ２ ２ ） 便可完成迭代。 关于 Ｔ ｈ ｏ ｍ ａ ｓ 加速算法还有其它形式， 根据经验， 本文 中的这种形式对于结构内塑性区很大的情况具有很强的 适应性， 并且这种形式更符合传统的有限元格式。
Ι ｕ Δｌ －Δｕ Δｕ ｉ － １ （ ｉ － １ ＋δ ｉ） Ｔ Π ｕ Δｕ ｉ － １ δ ｉ ２ Ｔ
（ ２ ９ ）
式（ ２ ９ ） 即为弧长法荷载增量计算公式。根据不同的简化 可根椐需要选用。 条件， λ ｉ的表示方法也有好几种， 线性搜索在进入软化的有限元分析中十分重要， 将式
（ １ ６ ） （ １ ７ ）
ｐ 为塑性应变 ε 的函数 关， 本文中记 ｋ
（ １ ）
随着塑性变形的发生， 塑性屈服面随之改动。这里仅
收稿日期： ２ ０ １ ０－ ０ １－ ０ ７ 作者简介： 夏怀孝（ １ ９ ７ ４ —） ， 男， 博士， 讲师， 主要从事自适应有限元研究。
夏怀孝： 软化材料的弹塑性有限元分析 １ ３ １
ｐ
（ ３ ２ ）
式中： ｗ为系统总势能。 φ （ ） 为一高度非线性且不可显 Δ ρ 示的函数， 一般利用线性搜索确定 ρ 的大致范围， 即： （ ０ ） φ ρ＝－ （ １ ）－φ （ ０ ） φ
Ｔ Ｔ （ ０ ） ＝δ ｕ （ ） ＝ －δ ｕ ｕ φ Δｕ δ ｉｇ ｉ － １ ｉＫ ｉ － １ ｉ
ｐ ５ ）令 σ ′＝ σ ， ｋ ′＝ ｋ ， ′ ｄ ′＝ σ ， 由 ε σ σ ｅ ０ ｐ＝ ε ０， ｅ－ ０
３ 弧长法控制的迭代
在非线性有限元分析中， 荷载是通过若干小的增量步
{
＜０ 卸载 ９ ） （
进行施加的。在每个小的增量步中， 通过若干次迭代达到 为某增量步的第 ｉ 次迭代， ｕ 平衡。设 ｉ δ ｉ为该迭代步位移 增量， 则该荷载增量步第 ｉ 次迭代的位移可记为 （ １ ０ ） ｕ ｕ ｕ Δ ｉ ＝Δ ｉ － １ ＋δ ｉ 弧长法的约束条件记为 Δｕ Δｕ ｉ － １ ｉ ＝ Δｌ ｌ 保持不变， ｕ Δ δ ｉ可表示为
Ｐ ： （ Δ σ ＝Ｄ Δ ε－Δ ε）
Ｔ ｈ ｏ ｍ ａ ｓ 假定
ｉ ｉ ｉ － １ ｉ ｉ ａ Δｕ 珘 ｅ ≈ ａ Δｕ ｅ ＋Δ ｕ ｅ
（ １ ８ ） １ ９ ） （ （ ２ ０ ）
（ ４ ）
其中
ｉ － １ ｉ － １ ｉ － １ ｉ } {Ｒ ｕ Δ 珘 （ ｕ ＋ａ Δｕ ｅ ｅ ＝Ｋ ｅ）
ｉ ｉ － １ ｉ ｉ ｅ － １ ｉ － １ ｅ ｉ ｅ
（ １ ２ ） （ １ ３ ） （ １ ４ ） （ １ ５ ）
（ ２ ６ ） （ ２ ７ ） （ ２ ８ ）
∫
如Ｆ （ ， ｋ ）＜ Ｔ Ｏ Ｌ ， 转入下一荷载步； 否则令 σ ′＝ σ ， ｋ ′＝ σ
将式（ ２ ６ ） 代入式（ ２ ４ ） 可得 λ ｉ＝
（ ３ ４ ） （ ３ ５ ） （ ３ ６ ）
对于严格的凸问题， 线性搜索要求 （ ０ ）≤ ０ ， φ （ １ ） ＞φ （ ０ ） φ ３ ５ ） 、 式（ ３ ６ ） 可知， 当Ｋ 为非正定（ 结构进入软化） ， 由式（ ｉ － １
第３ １卷 第 ３期 四 川 兵 工 学 报 ２ ０ １ ０年 ３月 【 其他研究】
软化材料的弹塑性有限元分析
夏怀孝
（ 三峡大学 土木水电学院， 湖北 宜昌 ４ ４ ３ ０ ０ ２ ）
摘要： 给出了软化材料弹塑性有限元分析的基本算法：选用显式 Ｅ ｕ ｌ ｅ ｒ 法作为本构积分方案；在结构未达到极限 荷载前采用初应力法配合 Ｔ ｈ ｏ ｍ ａ ｓ 加速方案作为平衡迭代方法； 而当结构接近极限荷载时， 程序转换为弧长法控 制迭代， 使结构越过极值点进入软化区直至破坏。最后， 讨论了该领域面临的重要问题— — —由于材料软化而导 致的应变局部化现象及其有限元模拟。 关键词： 有限元分析； 弹塑性； 软化； 应变局部化 中图分类号： Ｏ ３ ４ ４ ． ３ 文献标识码： Ａ 文章编号： １ ０ ０ ６－ ０ ７ ０ ７ （ ２ ０ １ ０ ） ０ ３－ ０ １ ３ ０－ ０ ３ ｌ ｉ ｍ ｉ ｔ 题， 即如何求 得 极 限 荷 载， 以及如何越过该极值点（ ｐ ｏ ｉ ｎ ｔ ） 使受力—位移过程得以继续发展。这 ２个问题一直 是非线性有限元分析的难点， 研究者也提出了许多方法进 行解决， 如虚拟弹簧法、 位移控制法、 强制迭代法等， 关于
ｐ ｋ＝ｋ （ ε）
（ ３ ）
ｉ － １ ｉ Ｒ 为上一迭代末的不平衡力； ａ 为加速因子。 式中：
假定塑 性 势 与 塑 性 屈 服 面 重 合， 那么根据正交流动法 则， 有 Ｆ ｐ Δ ε ＝Δ λ σ 假设塑性屈服面连续变化， 对式（ １ ） 进行微分可得 Ｆ ｋ ： Ｆ Δ σ ＋Ｆ λ ＝０ ｐ ｋ σ： ε σΔ 由增量弹塑性应力应变关系可得
０
变增量 Δ ， 根据式（ ３ ） 求出 ｋ 。 ε ０ ２ ）判断该荷载步为加载或卸载。由于应力空间中加 载、 卸载准则仅适用于强化阶段， 故采用应变空间中加载 卸载准则： Ｆ ： ｄ ε ＝ ＝０ 中性变载 ε ＞０ 加载 Ｆ Ｆ Ｆ ： ｄ Ｄ ： ｄ ε＝ ： ε＝ ： Δ σ ｅ ε σ σ σ σ ｅ ＝σ ０ ＋Δ ｅ ３ ）若材料处于卸载， ， 转入下一荷载步。 σ＝ σ ｅ ４ ）若材料处于加载， 如Ｆ （ ， ｋ ）＜ ０ ， 则 σ＝ ， 转入 σ σ ｅ ０ ｅ 下一荷载步。
－ １ Ｔ ｕ ｑ－ Ｂ ｄ ｖ ） δ λ σ ｉ ＝Ｋ ｉ － １ （ ｉ ｉ － １ Ｔ ２
（ ２ ３ ） （ ２ ４ ）
式中： Δ σ ｅ为弹性应力增量。相应定义弹性试探应力为 （ １ １ ）
其中： ｌ 为每一增量步的弧长。在该增量步的所有迭代中 Δ
∫
（ ２ ５ ）
其中： Ｋ 为该迭代步的刚度矩阵； ｑ 是单位荷载增量； λ ｉ － １ ｉ 为本迭代步荷载。在弧长法中， ｕ δ ｉ的表示形式有好几种， 这里的表示方法可以较好地与后文将提到的线性搜索相 配合。 ｕ 部分： 将δ ｉ分为 ２
令 （ ５ ） （ ６ ）
ｉ ｉ － １ ａ ＝ａ ＋ ｉ ｉ ｉ ｉ － １ ｉ ε ＝ ａΔｕ 珘 ｅ －ａ Δ ｕ ｅ －Δ ｕ ｅ
利用最小二乘法求出 ａ， 使得 ε 最小
ｉ { Δｕｉ }Ｔ{ Δｕ 珘} ｅ ｅ ｅ ｅ ｉ { Δｕｉ }Ｔ{ Δｕ }
ｉ
定义塑性模量 Ｈ为 Ｈ ＝－Ｆ ｋ ： Ｆ ｐ ｋ ε σ 则由式（ ４ ）～ （ ６ ） 可得 Ｆ Ｄ ： Δ ε σ： Δ λ＝ Ｆ ： Ｄ ： Ｆ ＋Ｈ σ σ ｕ ｌ ｅ ｒ 本构积分的算法。 下面具体给出 Ｅ
［ ４ ］ ｅ ｒ ｇ ａ ｎ 也有详细评述。众多研究表明， Ｒ ｉ ｓ ｋ ｓ 和 这个问题 Ｂ ５ － ６ ］ Ｃ ｒ ｉ ｓ ｆ ｉ ｅ ｌ ｄ 等［ 提出的弧长法是较为有效的方法。该法在
近 ３ ０年来， 弹塑性有限元分析一直是国际数值计算领 域研究的重点。随着工程规模的扩大和相应研究的深入， 结构的极限荷载及其越过极限荷载后的破坏情况等逐渐 引起了工程界的高度重视。因此， 含软化阶段的弹塑性有 限元分析也随之被引入。 本构积分是弹塑性分析最重要的环节之一， 它对计算 结果的精度和计算的收敛性均有着重要影响。弹塑性本 构积分方案大致有显式法和隐式法 ２种。隐式法由于能给 出一致切线刚度矩阵（ ｃ ｏ ｎ ｓ ｉ ｓ ｔ ｅ ｎ ｔ ｔ ａ ｎ ｇ ｅ ｎ ｔ ｓ ｔ ｉ ｆ ｆ ｎ ｅ ｓ ｓ ｍ ａ ｔ ｒ ｉ ｘ ） ， 且 对于 Ｎ ｅ ｗ ｔ ｏ ｎ Ｒ ａ ｐ ｈ ｓ ｏ ｎ法可达到二阶收敛速度， 因而该方法
（ ２ ３ ） 改写为 ｕ ｕ ｕ Δ δ ｉ ＝Δ ｉ － １ ＋ρ ｉ （ ３ ０ ）
ｕ ＝ ｕ ＋ａΔｕ
四川兵工学报 １ ３ ２
式中： 为加速因子。记不平衡力为 ｇ ， 有 ρ
Ｔ ｇ ｄ ｖ －λ ｑ ｉ － １ ＝ Ｂ σ ｉ － １ ｉ