用于有限元模拟弹塑性分析的统一能量优化

《Computers & Structures》：Unified energy optimization for finite elements simulating elasto-plastic analysis

　　

摘要 基于增量迭代方法的弹塑性有限元结构行为分析并非易事，且通常涉及复杂的数值计算，无法保证收敛性。替代方案是应用强大且稳定的凸优化求解器来最小化余能或势能，这可能为非线性超弹性材料在给定单调加载情况下的解提供一步完成的方法，从而模拟弹塑性材料行为。通过平衡和协调的弱形式表述，本文展示了两种最小化方法的新型表述互为对偶，产生相同的优化解，且对偶参数的解释是显而易见的。这种统一的混合解决方案能得出应力、应变和位移的准确结果，通常优于通过严格余能或势能方法获得的结果。提供了钢和钢筋混凝土结构的验证及引言平面结构示例，显示了快速收敛。此外，提供了一个钢筋混凝土梁建模的综合结构示例，显示了令人信服的结果。

引言 考虑非线性材料行为的结构有限元分析通常基于增量方法，结合虚功原理和弱形式表述。一般来说，增量求解策略考虑更新的线性化控制方程和平衡迭代，以便详细映射加载路径，参见例如 [1], [2]。因此，结构行为的非线性有限元分析并非易事，且通常涉及复杂的数值计算，这些计算可能难以控制，可能导致收敛问题。

在有限元的早期阶段，Maier [3] 建议利用二次规划来建模桁架结构的弹塑性行为。可能由于当时求解器能力非常有限，这种建模方法未被进一步追求。后来 [4], [5] 和 [6] 将二次规划应用于增量求解方法中的弹塑性步骤建模。最近，Vestergaard 等人 [7] 和 Ferradi 等人 [8] 表明，非线性增量结构分析的替代方案可以是一步模拟的弹塑性分析。他们的方法基于平衡单元和最小化余能，应用二阶锥规划。

对于非线性弹性材料，基于最小化余能的解析解与基于最小化势能的解是相同的，这是一个众所周知的事实。对于有限元解，这通常不成立。在此，提出了一种新的有限元解决方案，它包含平衡和协调要求的弱形式表述，以及用于本构计算的特定积分描述。这种有限元建模统一了基于余能和势能原理的模型，产生相同的解。

基于最小化余能或势能的方法的一个后果是，非线性材料行为被建模为超弹性，这排除了建模滞回的可能性。因此，该方法的局限性在于仅适用于单调加载，对于此类情况，建模非常高效。模拟弹塑性分析 (SEPA) 的一个主要优势是它可以在一步内提供给定载荷水平的解。对于需要弹塑性分析的设计情况，SEPA 将是一个稳定且有效的选择。特别是对于钢筋混凝土结构的设计，由于组合相互作用使其成为一项困难的任务，SEPA 是一个显而易见的选择。

SEPA 高效的一个主要原因是凸优化求解器的稳定性和速度，这主要归功于 Karmarkar [9] 引入的内点法。如今，凸优化领域已非常成熟，参见例如 [10]，并且有现代稳定的求解器如 [11] 可用。Vestergaard 等人 [7] 提出的方法基于最小化余能，应用带有应力参数的严格平衡单元 [12]。然而，在本文中，两种不同方法（分别最小化余能或势能）被独立表述，但方式使得它们可以证明互为对偶。为实现这一目标，应用了平衡和协调要求的弱形式表述。这种表述意味着任何一种方法都会得出完全相同的结果。此外，它允许对对偶参数进行清晰且显而易见的解释。这种表述规避了 Vestergaard 等人 [7] 应用的严格平衡单元所看到的困难，例如过度约束的边界条件和运动学解释中虚假模式的出现，参见 [13]。此外，基于弱形式的三角形单元通常显示出比严格平衡单元更快的收敛性，参见 [14]。最近，Larsen 等人 [15] 有效地应用最小势能原理来优化钢筋混凝土面板。一般来说，目前的统一混合解决方案得出的结果优于通过严格余能或势能方法获得的结果。

这里展示的工作借鉴了 Poulsen 和 Olesen [14], [16] 提出的有限元极限分析 (FELA) 之前的工作，其中考虑了刚性/理想塑性材料行为。这些平行点包括平衡的弱形式、一些协调要求以及屈服准则的二阶锥表述。

本文讨论平面结构问题，即处于平面应力状态的面板。原则上，材料行为被假设为多线性的，即由连续的线性弹性阶段组成。然而，为了理论部分的清晰，这里仅考虑两个阶段。另一方面，这将允许建模不同的有用情况，并与有限元极限分析 (FELA) 的结果进行直接比较。

从一个阶段到另一个阶段的过渡由作为屈服准则的凸关系给出。在本表述中，对于平面结构问题，屈服准则可以有效地使用二阶锥规划来表述。这里，两阶段材料模型应用于具有 von Mises 屈服准则 [17] 的钢结构分析，以及具有 Mohr-Coulomb 破坏准则 [13] 和钢筋简单线性屈服准则的钢筋混凝土结构分析。

基于 von Mises 的钢模型被实现，以说明在单调加载情况下的准确性和有效性。由于 von Mises 准则的简单性，它也被用于验证有限元实现和理论方面，特别是原始参数和对偶参数之间的独特耦合。

由于材料的组合结构，钢筋混凝土结构在非线性的开裂阶段难以分析。这导致了诸如撑杆 - 拉杆法 [18] 等简化方法，以及具有清晰物理意义并适用于弹塑性分析的应力场法 [19]。在这方面，本表述最大的潜力可能在于钢筋混凝土结构领域。

为了说明本方法的优势，分析了一个带有孔洞的钢筋混凝土梁的综合示例。基于结果，该示例进一步分析了后处理中的裂缝形成。尽管假设在全局应力分布中抗拉强度为零，但在考虑实际抗拉强度和恒定粘结应力 [20] 的子模型中，可以解决裂缝分布和裂缝宽度。这种方法给出了现实的裂缝宽度，避免了局部重新分布期间的迭代；然而，它这样做并没有包括受拉刚度效应。

本文首先陈述能量优化原理，随后展示所应用的非线性本构模型，形式为连续的线性阶段。然后发展有限元的余能方法，随后展示有限元的势能方法。这种新的统一能量优化从这些计算能量方法的一对一对偶性中展示出来，随后是一种有效的计算表述。论文的其余部分是数值模型的验证和应用说明，该模型基于六节点三角形有限元。

章节片段 能量优化原理
余能和势能原理都可以形成优化表述的基础，在最优时提供满足特定加载下所有场和边界条件的解。它们可以表述如下：
余能方法 – 任何满足平衡方程和应力边界条件，并使余能达到最小值的应力状态，将确保协调条件的满足。
本构模型
应力 - 应变关系可以是多线性的，如图 1 所示的三线性情况，具有三个应力等级区间。应力状态和应变状态均表示为线性贡献的总和，σ=∑σi, ε=∑εi，对于区间 i，应力 - 应变关系由下式给出
εi=(ci?ci?1)σi
其中 c0=0
其中柔度矩阵 ci 由表征每个等级 i 的刚度矩阵 di 的逆给出，即 ci=di?1。
应力等级之间的过渡 有限元插值
应用有限元方法以建立离散的最小余能和势能优化问题。这些表述是通用的，这里它们应用于平面问题。
单元中的位移场和应力场，u 和 σ，分别从节点应力值，v 和 S 插值：
u=Nuv
σ=NσS
其中 Nu 和 Nσ 是插值矩阵。
应变插值矩阵 B 通过位移场经线性应变定义获得，
弱形式平衡
平衡以弱形式强制执行，类似于 [14]，要求对于任何运动学上允许的虚位移场 δv，虚内部功和虚外部功相等。虚内部功和虚外部功的单元贡献是：
δwint=∫ΩeδεTσdΩe=δvT∫ΩeBTNσdΩeS
δwext=∫ΩeδuTpdΩe=δvT∫ΩeNuTpdΩe
其中 p 是单元上的载荷。(13) 中的积分涉及单元平衡矩阵 h 的定义，(14) 中的积分产生单元节点载荷
弱形式协调
协调要求分为两部分：首先位移必须与应变协调，其次塑性应变必须与转换的塑性应变协调，后者进入锥限制。
位移和应变之间的协调以弱形式强制执行，要求独立的应变场 ε 产生与由 (12) 给出的位移场导出的应变场相同数量的虚内部功。在
能量方法的对偶性
基于附录 A，并预先识别对偶参数 yT=[?V,E?p]T，(36) 的对偶可以写为：
minimize12STCeS+12(S?Se)TCep(S?Se)+CTE?p?RTV
s.t.[CeS+Cep(S?Se)?Cep(S?Se)]+[HT00TT][?VE?p]=[00]
E?p?C0
利用
CeS=MEe
和
Cep(S?Se)=CepSp=MEp
以下解释是有效的：
[CeS+Cep(S?Se)?Cep(S?Se)]=[MEe+MEp?MEp]=[ME?MEp]
其中引入了总节点应变向量，E=Ee+Ep。
交换 (62) 中目标函数里余能的表达式
计算表述
当二次最小化问题包含锥限制时，以锥限制参数集的形式重写二次目标函数在计算上可能更有利。因此，在实现余能方法时，(36) 中的目标函数根据以下式子重写：
Wc=12STCeS+12(S?Se)TCep(S?Se)=[SSe]T(12[Ce+Cep?Cep?CepCep])[SSe]=[SSe]TTcTTc[SSe]
其中用 Cholesky 分解替换二次对称矩阵，使能量
三角形有限元
应用厚度恒定的六节点三角形单元进行数值建模，参见图 2。应用广义应力，形式为应力乘以单元厚度。三个顶点节点处的广义应力分量线性插值，所有六个节点处的位移分量二次插值。
单元的实现描述在 Poulsen 和 Olesen [14] 中，其中单元矩阵的插值和积分是
钢 - von Mises
为了验证模拟弹塑性材料行为的实现，建立了一个简单示例，即 4m×4m 的正方形板，在面内单轴拉伸 p=400MPa 下加载。该板用 2 个三角形单元建模，应用以下参数值：t=1m, Ee=210GPa, ν=0, fy=300MPa, Ep=Ee/2, νh=0。为了获得一个良好条件的优化问题，所有数据均以 SI 单位应用，其中长度乘以 100 (cm)，力除以 1000 (kN)。这
深梁 - 钢
在本节中，第一个示例考虑图 6 左侧所示的深梁。

该材料由钢材制成，具有双线性应力 - 应变关系，其中第一阶段以弹性模量 Ee=E 为特征，第二阶段以 Ep=10?6?E 为特征，称为硬化模量。这模拟了线弹性 - 理想塑性材料，从而能够与 FELA 解 [14] 进行比较。该面板采用以下参数值进行建模：h=12?=1m，t=0.04m，E=210GPa，ν=结论。

一种替代增量迭代方法，用于通过有限元分析解决非线性结构行为，被提出。本方法基于互补能或势能的最小化，应用强大的凸优化稳定求解器，有可能针对给定的载荷情况在一步内提供模拟的弹塑性解。作为一种创新，这两种最小化方法的统一公式被证明恰好是彼此的。

CRediT 作者贡献声明
Peter Noe Poulsen：撰写 – 初稿，可视化，验证，软件，方法论，调查，形式分析。John Forbes Olesen：撰写 – 初稿，可视化，验证，软件，方法论，调查，形式分析。

利益冲突声明
作者声明，他们没有已知的竞争性财务利益或可能影响本文报告工作的个人关系。

Peter Noe Poulsen|John Forbes Olesen