CN117390778A - 一种考虑薄板应变强化效应的塑性安定上限载荷计算方法 - Google Patents

Abstract

本发明提出一种考虑薄板应变强化效应的塑性安定上限载荷计算方法，属于物理技术领域，包括如下步骤：ST1、准备计算数据，得到薄板各个三角形子域对应的光滑广义应变‑位移速度关系矩阵；ST2、计算广义弹性应力场，分别得到角点载荷和恒定载荷作用时的光滑广义弹性应力场；ST3、在初始迭代时，假设整个薄板处于非屈服状态，求解线性方程组；ST4、在第h(h≥1)次迭代时，根据第h‑1次迭代的计算结果，集成和求解相应的线性方程组，最终得到第h次迭代时的安定上限载荷乘子，判断是否终止迭代；ST5、计算结果后处理。本发明解决了考虑薄板应变强化效应安定上限分析极小化迭代格式的建立与线性化求解问题，实现了考虑薄板应变强化效应安定上限载荷的精确高效稳定求解。

Description

一种考虑薄板应变强化效应的塑性安定上限载荷计算方法

技术领域

本发明属于塑性力学技术领域，具体涉及一种考虑薄板应变强化效应的塑性安定上限载荷计算方法。

背景技术

薄板是航天、航空、核能、船舶等领域广泛使用的一种结构形式。精确可靠地计算出薄板结构在循环载荷作用下的塑性极限承载能力，对板结构的优化设计和安全评定有着非常重要的意义。安定分析是一种研究结构在循环变化载荷作用下塑性失效的直接有效方法，不涉及到弹塑性增量分析所需要的载荷变化历史，具有很强的实用性和可操作性。与常规的弹性分析相比，安定性分析所获得的安定载荷可以作为反映结构塑性失效的重要参数，更能反映结构的实际承载潜力和安全程度，已经在许多工程结构的设计规范和法规中得到了广泛应用。

目前，安定分析的理论研究和工程应用已经有了长足的发展，但是薄板结构安定分析的研究工作还相对较少。薄板安定分析的理论和实验研究受到结构几何、边界和作用载荷等条件的限制，目前仅有极少量的研究见诸报导。数值计算能准确求解相对复杂薄板的极限承载力，因而一直是该领域的研究热点之一。许多学者将成熟的数值分析技术与安定分析的数学规划理论相结合，重点研究薄板安定分析的高效优化算法和数值求解程序，并且取得了一些令人关注的研究成果。由于许多金属或合金具有明显的应变强化性质，且材料的应变强化效应对结构的安定性影响很大，忽略应变强化效应得到的安定载荷偏于保守。然而，目前这些有限的薄板安定分析的研究成果中，都没有考虑应变强化效应对薄板安定性的影响。因此，非常有必要使用精度好、效率高的数值计算方法开发出精确高效的求解程序，来研究材料应变强化效应对薄板结构安定失效的影响。

发明内容

本发明提供一种考虑薄板应变强化效应的塑性安定上限载荷计算方法，基于具有计算精度和效率高、数值稳定好、后处理方便等优良性能的C¹节点自然单元法(C¹-nodal-NEM)，目的是解决考虑薄板应变强化效应安定上限分析极小化格式的建立与线性化求解问题，实现考虑薄板应变强化效应安定上限载荷的精确、高效和稳定求解。

本发明的目的是通过如下技术方案实现的：

一种考虑薄板应变强化效应的塑性安定上限载荷计算方法，包括如下步骤：

ST1、准备计算数据，得到薄板结构各个三角形子域对应的光滑广义应变-位移速度关系矩阵；

ST2、计算广义弹性应力场，分别得到角点载荷和恒定载荷作用时薄板结构各个三角形子域对应的光滑广义弹性应力场；

ST3、在初始迭代时(h＝0)，假设整个薄板结构处于非屈服状态，求解线性方程组，依次得到初始迭代时的拉格朗日乘子、残余位移增量、每个角点载荷作用时各个三角形子域对应的光滑广义塑性应变、考虑薄板应变强化效应的安定上限载荷乘子；

ST4、在第h(h≥1)次迭代时，根据第h-1次迭代的计算结果，集成和求解相应的线性方程组，依次得到第h次迭代时的拉格朗日乘子、残余位移增量、每个角点载荷作用时各个三角形子域对应的光滑广义塑性应变、考虑薄板应变强化效应的安定上限载荷乘子，根据迭代收敛条件判断是否终止迭代计算；

ST5、计算结果后处理，得到薄板结构各个节点的塑性耗散功率；

作为进一步优化，准备计算数据步骤具体包括：

ST1.1、准备薄板结构的离散节点、Delaunay三角形、位移边界节点、载荷(恒定载荷、变化载荷)、几何参数、材料参数(杨氏模量E、泊松比v、屈服应力Y、强化参数)、偏移系数γ、误差容限vol1和vol2；

ST1.2、根据薄板的节点、Delaunay三角形和位移边界节点的信息，采用C¹节点自然单元法依次确定围绕每个节点x_r(r＝1～NP，NP为离散节点的总数)Voronoi子域S_r的顶点及其逆时针排列顺序，Voronoi结构顶点微小偏移后的坐标，每个三角形子域S_rs(RS为子域S_r的三角形子域总数)的面积A_rs、每条边的切向矢量n₁和法向矢量n₂、每条边上每个积分点的自然临近节点I(I＝1～np，np为采用空圆准则确定的每个积分点的自然临近节点总数)及其C¹自然临近形函数的一阶导数/> x_rs为对应于三角形子域S_rs的虚拟节点，上标w、θ_x、θ_y分别是薄板的挠度、x方向的旋度、y方向的旋度，下标x、y分别表示对x方向和y方向求一阶导数。进而计算每个积分点自然临近节点对应的光滑广义应变-位移速度关系矩阵/>

式中，Γ_rs为子域S_rs的边界；

ST1.3、“对号入座”集成各个三角形子域S_rs全部积分点自然临近节点对应的矩阵(/> 为合并全部积分点相同自然临近节点后的节点总数)，得到各个三角形子域S_rs对应的光滑广义应变-位移速度关系矩阵/>

作为进一步优化，计算广义弹性应力场步骤具体包括：

ST2.1、根据各个三角形子域对应的光滑广义应变-位移速度关系矩阵和弹性本构关系，计算和集成薄板的整体弹性刚度矩阵K^e：

式中，D_b为弹性关系矩阵；

为薄板的弯曲刚度；

h₀为板厚；

ST2.2、根据作用在薄板上每组独立变化载荷的范围和组数确定角点载荷和其总数/>对于每个角点载荷中的均布力，按照Voronoi子域S_r的面积/>将其等效为作用在节点上的等效载荷列阵/>对于每个角点载荷中的集中力，将其直接施加在所作用的节点上。由各个节点的等效载荷列阵/>“对号入座”集成每个角点载荷作用时薄板的整体弹性载荷列阵/>

ST2.3、引入位移边界条件修改薄板的整体弹性刚度矩阵K^e和每个角点载荷作用时的整体弹性载荷列阵求解薄板弹性问题的线性控制方程/>得到每个角点载荷作用时薄板的弹性位移场/>

ST2.4、根据弹性本构关系和弹性位移场，计算每个角点载荷作用时各个三角形子域S_rs对应的光滑广义弹性应力

ST2.5、仿照步骤ST2.2～2.4，计算恒定载荷作用时各个三角形子域S_rs对应的光滑广义弹性应力

作为进一步优化，在初始迭代时，假设整个薄板结构处于非屈服状态，求解线性方程组步骤具体包括：

ST3.1、假设整个薄板结构处于非屈服状态，取和/>依次计算出初始迭代时各个三角形子域S_rs对应的中间变量/>

式中，M_P＝Yh₀ ²/4为塑性极限弯矩；

Q是正定对称的常数矩阵，且Q^-1＝2D/3，D为处理薄板塑性不可压条件引入的常数矩阵：

ST3.2、计算和集成整体恒定载荷列阵F₀、初始迭代时考虑薄板应变强化效应的整体刚度矩阵K⁰和整体变化载荷列阵F₁ ⁰：

ST3.3、引入位移边界条件修改整体恒定载荷列阵F₀、初始迭代时考虑薄板应变强化效应的整体刚度矩阵K⁰和整体变化载荷列阵F₁ ⁰，分别求解线性方程K⁰(Δa₀)⁰＝F₀和K⁰(Δa₁)⁰＝F₁ ⁰，得到初始迭代时残余位移增量的中间变量(Δa₀)⁰和(Δa₁)⁰；

ST3.4、求解初始迭代时的拉格朗日乘子λ⁰：

式中，中间变量和/>分别为：

ST3.5、求解初始迭代时的残余位移增量(Δa)⁰和在每个角点载荷作用下各个三角形子域S_rs对应的光滑广义塑性应变

(Δa)⁰＝(Δa₀)⁰+λ⁰(Δa₁)⁰

ST3.6、根据得到的残余位移增量和光滑广义塑性应变，依次计算出初始迭代时各个三角形子域S_rs对应的变量和μ¹(x_rs)、/>

式中，β₁和β₂均为远小于1的正数，且分别设定为：

进而得到初始迭代时考虑薄板应变强化效应的安定上限载荷乘子s⁰：

作为进一步优化，在第h(h≥1)次迭代时，根据第h-1次迭代的计算结果，集成和求解相应的线性方程组步骤具体包括：

ST4.1、根据第h-1次迭代的计算结果确定出变量和/>的值，依次计算出第h次迭代时各个三角形子域S_rs对应的中间变量/>

ST4.2、计算和集成第h次迭代时考虑薄板应变强化效应的整体刚度矩阵K^h和整体变化载荷列阵F₁ ^h：

ST4.3、引入位移边界条件修改整体恒定载荷列阵F₀、第h次迭代时的整体刚度矩阵K^h和整体变化载荷列阵F₁ ^h，分别求解线性方程K^h(Δa₀)^h＝F₀和K^h(Δa₁)^h＝F₁ ^h，得到第h次迭代时残余位移增量的中间变量(Δa₀)^h和(Δa₁)^h；

ST4.4、求解第h次迭代时的拉格朗日乘子λ^h：

式中，中间变量和/>分别为：

ST4.5、求解第h次迭代时的残余位移增量(Δa)^h和在每个角点载荷作用下各个三角形子域S_rs对应的光滑广义塑性应变

(Δa)^h＝(Δa₀)^h+λ^h(Δa₁)^h

ST4.6、根据得到的残余位移增量和光滑广义塑性应变，依次计算出第h次迭代时各个三角形子域S_rs对应的变量和μ^h+1(x_rs)、/>

从而计算出第h次迭代时考虑薄板应变强化效应的安定上限载荷乘子s^h：

ST4.7、根据设定的误差容限vol1、vol2和下列收敛条件，判断是否终止迭代：

||(Δa)^h-(Δa)^h-1||/||(Δa)^h-1||≤vol1，|s^h-s^h-1|/s^h-1≤vol2。

作为进一步优化，计算结果后处理步骤具体包括：

ST5.1、假设第h次迭代终止时满足收敛条件，计算终止迭代时各个三角形子域S_rs对应的光滑塑性耗散功率D^h(x_rs)：

ST5.2、根据各个三角形子域S_rs与Voronoi子域S_r的面积关系，得到薄板结构各个节点的塑性耗散功率D^h(x_r)：

本发明所取得的有益技术效果是：

1、本发明发展的计算方法具有格式较为简单、通用性强、易于程序实现、计算精度和效率高、数值稳定性好等优点，能够准确获取不同强化模型下薄板结构的安定上限载荷，弥补了目前薄板结构安定分析相对较少、且尚未考虑应变强化效应的不足，是对薄板结构塑性安定分析的进一步丰富和发展。

2、本发明采用的C¹节点自然单元法(C¹-nodal-NEM)兼具C¹自然临近插值和三角形子域稳定相容积分的优点，只需要计算C¹自然临近形函数的一阶导数、而不需要计算其二阶导数，非常有利于构造薄板安定上限分析的高精度光滑广义应变场，并且数值计算结果的后处理方便。

3、通过本发明获得的薄板不同材料模型的安定载荷和能真实反映结构安定失效模式的塑性耗散功率云图，既能为薄板结构的工程设计和安全评估提供了借鉴，也能为进一步深入研究应变强化效应对薄板结构安定失效的影响提供参考。

附图说明

图1是本发明的实施流程框图；

图2是本发明的三角形子域划分示意图；

图3是本发明具体实施例中进行安定上限分析的夹支圆形薄板受均布力q作用示意图；

图4是图3中的夹支圆形薄板801个离散节点示意图；

图5是图3中的夹支圆形薄板1536个Delaunay三角形分布示意图；

图6是图3中的夹支圆形薄板安定上限载荷乘子的迭代收敛曲线图；

图7是图3中的夹支圆形薄板在材料模型工况0≤q≤q^max时的塑性耗散功率云图(10⁶N·m)；

图8是图4中的夹支圆形薄板在材料模型工况0≤q≤q^max时的塑性耗散功率云图(10⁶N·m)；

图9是图5中的夹支圆形薄板在材料模型和/>工况-q^max≤q≤q^max时的塑性耗散功率云图(10⁶N·m)。

具体实施方式

下面结合附图和具体实施方式对本发明的技术方案做进一步详细说明。显然，所描述的实施例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明要求保护的范围。

一种考虑薄板应变强化效应的塑性安定上限载荷计算方法具体实施例，根据安定分析的Koiter定理和两面屈服准则，建立了考虑薄板应变强化效应的塑性安定上限分析数学规划格式，并且采用一种具有计算精度和效率高、数值稳定好、后处理方便等优良性能的C¹节点自然单元法(C¹-nodal-NEM)，构造薄板的光滑广义塑性应变增量，对该含有时间积分、带等式约束的非线性极小化规划问题进行离散。采用理论处理载荷域内各个角点载荷对应的广义塑性应变增量，来消除规划格式中的时间积分；建立一套直接迭代格式来线性化和求解该非线性问题。在每次迭代求解中，根据上一次的计算结果对目标函数和约束条件进行相应的修改，实现了非线性目标函数的线性化，将每次迭代计算转化为一组线性方程组的求解。最终得到考虑薄板应变强化效应安定上限分析的光滑广义塑性应变场、安定上限载荷和塑性耗散功率。

如图1所示，本具体实施例中考虑薄板应变强化效应的塑性安定上限载荷计算方法包括如下步骤：

ST1、准备计算数据

ST1.1、准备薄板结构的离散节点、Delaunay三角形、位移边界节点、载荷(恒定载荷、变化载荷)、几何参数、材料参数(杨氏模量E、泊松比v、屈服应力Y、强化参数)、偏移系数γ、误差容限vol1和vol2。

ST1.2、根据薄板的节点、Delaunay三角形和位移边界节点的信息，采用C¹节点自然单元法依次确定围绕每个节点x_r(如图2所示，r＝1～NP，NP为离散节点的总数)Voronoi子域S_r的顶点及其逆时针排列顺序，Voronoi结构顶点微小偏移后的坐标，每个三角形子域S_rs(RS为子域S_r的三角形子域总数)的面积A_rs、每条边的切向矢量n₁和法向矢量n₂、每条边上每个积分点的自然临近节点I(I＝1～np，np为采用空圆准则确定的每个积分点的自然临近节点总数)及其C¹自然临近形函数的一阶导数/> x_rs为对应于三角形子域S_rs的虚拟节点，上标w、θ_x、θ_y分别是薄板的挠度、x方向的旋度、y方向的旋度，下标x、y分别表示对x方向和y方向求一阶导数。进而计算每个积分点自然临近节点对应的光滑广义应变-位移速度关系矩阵

式(1)中Γ_rs为子域S_rs的边界。

ST1.3、“对号入座”集成各个三角形子域S_rs全部积分点自然临近节点对应的矩阵(/> 为合并全部积分点相同自然临近节点后的节点总数)，得到各个三角形子域S_rs对应的光滑广义应变-位移速度关系矩阵/>

ST2、计算广义弹性应力场

ST2.1、根据各个三角形子域对应的光滑广义应变-位移速度关系矩阵和弹性本构关系，计算和集成薄板的整体弹性刚度矩阵K^e：

式(4)中，D_b为弹性关系矩阵，为薄板的弯曲刚度，h₀为板厚。

ST2.2、根据作用在薄板上每组独立变化载荷的范围和组数确定角点载荷和其总数/>对于每个角点载荷中的均布力，按照Voronoi子域S_r的面积/>将其等效为作用在节点上的等效载荷列阵/>对于每个角点载荷中的集中力，将其直接施加在所作用的节点上。由各个节点的等效载荷列阵/>“对号入座”集成每个角点载荷作用时薄板的整体弹性载荷列阵/>

ST2.3、引入位移边界条件修改薄板的整体弹性刚度矩阵K^e和每个角点载荷作用时的整体弹性载荷列阵求解薄板弹性问题的线性控制方程/>得到每个角点载荷作用时薄板的弹性位移场/>

ST2.4、根据弹性本构关系和弹性位移场，计算每个角点载荷作用时各个三角形子域S_rs对应的光滑广义弹性应力

ST2.5、仿照步骤ST2.2～2.4，计算恒定载荷作用时各个三角形子域S_rs对应的光滑广义弹性应力

ST3、初始迭代(h＝0)

ST3.1、假设整个薄板结构处于非屈服状态，取和/>依次计算出初始迭代时各个三角形子域S_rs对应的中间变量/>

式(7)和式(10)中，M_P＝Yh₀ ²/4为塑性极限弯矩；Q是正定对称的常数矩阵，且Q^-1＝2D/3，D为处理薄板塑性不可压条件引入的常数矩阵：

ST3.2、计算和集成整体恒定载荷列阵F₀、初始迭代时考虑薄板应变强化效应的整体刚度矩阵K⁰和整体变化载荷列阵F₁ ⁰：

ST3.3、引入位移边界条件修改整体恒定载荷列阵F₀、初始迭代时考虑薄板应变强化效应的整体刚度矩阵K⁰和整体变化载荷列阵F₁ ⁰，分别求解线性方程K⁰(Δa₀)⁰＝F₀和K⁰(Δa₁)⁰＝F₁ ⁰，得到初始迭代时残余位移增量的中间变量(Δa₀)⁰和(Δa₁)⁰。

ST3.4、求解初始迭代时的拉格朗日乘子λ⁰：

式(15)中的中间变量和/>分别为：

ST3.5、求解初始迭代时的残余位移增量(Δa)⁰和在每个角点载荷作用下各个三角形子域S_rs对应的光滑广义塑性应变

(Δa)⁰＝(Δa₀)⁰+λ⁰(Δa₁)⁰ (18)

ST3.6、根据得到的残余位移增量和光滑广义塑性应变，依次计算出初始迭代时各个三角形子域S_rs对应的变量和μ¹(x_rs)、/>

/>

式(21)中的β₁和式(23)中的β₂均为远小于1的正数，且分别设定为：

进而得到初始迭代时考虑薄板应变强化效应的安定上限载荷乘子s⁰：

ST4、第h(h≥1)次迭代

ST4.1、根据第h-1次迭代的计算结果确定出变量和/>的值，依次计算出第h次迭代时各个三角形子域S_rs对应的中间变量/>

ST4.2、计算和集成第h次迭代时考虑薄板应变强化效应的整体刚度矩阵K^h和整体变化载荷列阵F₁ ^h：

ST4.3、引入位移边界条件修改整体恒定载荷列阵F₀、第h次迭代时的整体刚度矩阵K^h和整体变化载荷列阵F₁ ^h，分别求解线性方程K^h(Δa₀)^h＝F₀和K^h(Δa₁)^h＝F₁ ^h，得到第h次迭代时残余位移增量的中间变量(Δa₀)^h和(Δa₁)^h。

ST4.4、求解第h次迭代时的拉格朗日乘子λ^h：

式(33)中的中间变量和/>分别为：/>

ST4.5、求解第h次迭代时的残余位移增量(Δa)^h和在每个角点载荷作用下各个三角形子域S_rs对应的光滑广义塑性应变

(Δa)^h＝(Δa₀)^h+λ^h(Δa₁)^h (36)

ST4.6、根据得到的残余位移增量和光滑广义塑性应变，依次计算出第h次迭代时各个三角形子域S_rs对应的变量和μ^h+1(x_rs)、/>

从而计算出第h次迭代时考虑薄板应变强化效应的安定上限载荷乘子s^h：

ST4.7、根据设定的误差容限vol1、vol2和下列收敛条件，判断是否终止迭代：

||(Δa)^h-(Δa)^h-1||/||(Δa)^h-1||≤vol1，|s^h-s^h-1|/s^h-1≤vol2 (43)

ST5、计算结果后处理

ST5.1、假设第h次迭代终止时满足收敛条件，计算终止迭代时各个三角形子域S_rs对应的光滑塑性耗散功率D^h(x_rs)：

ST5.2、根据各个三角形子域S_rs与Voronoi子域S_r的面积关系，得到薄板结构各个节点的塑性耗散功率D^h(x_r)：

下面以一组具体的数据为例，对上述具体实施例作出进一步地说明：

如图3所示，以考虑应变强化效应夹支圆形薄板分别在两种工况(考虑变化均布力的两种工况：0≤q≤q^max、-q^max≤q≤q^max，不考虑恒定载荷)、两种材料模型(理想弹塑性/>有限随动强化/>)下的安定上限分析为例，实施方式如下：/>

步骤1：准备薄板的离散节点(如图4所示，总数为NP＝801)、Delaunay三角形(如图5所示，总数为1536个)、位移边界节点(图4中圆弧上的节点，总数为64个)、变化载荷的基准值q＝1.0N、角点载荷数目l＝2、板厚h₀＝0.01m、半径R＝1.0m、杨氏模量E＝210GPa、泊松比v＝0.3、屈服应力Y＝200MPa、强化参数和/>偏移系数γ＝1.0×10^-5、误差容限vol1＝vol2＝1.0×10^-4等信息。

步骤2：计算并形成每个积分点自然临近节点对应的光滑广义应变-位移速度关系矩阵，并且“对号入座”集成全部三角形子域对应的光滑广义应变-位移速度关系矩阵。

步骤3：计算和集成薄板的整体弹性刚度矩阵和每个角点载荷作用时的整体弹性载荷列阵。

步骤4：引入位移边界条件修改薄板的整体弹性刚度矩阵和每个角点载荷作用时的整体弹性载荷列阵，求解得到每个角点载荷作用时薄板的弹性位移场和此时各个三角形子域对应的光滑广义弹性应力场。

步骤5：依次计算出初始迭代时各个三角形子域对应的相关中间变量，并且计算和集成初始迭代时考虑薄板应变强化效应的整体刚度矩阵和整体变化载荷列阵。

步骤6：引入位移边界条件修改初始迭代时考虑薄板应变强化效应的整体刚度矩阵和整体变化载荷列阵，求解相应的线性方程，依次得到初始迭代时残余位移增量的中间变量、拉格朗日乘子、残余位移增量、每个角点载荷作用下各个三角形子域对应的光滑广义塑性应变、考虑薄板应变强化效应的安定上限载荷乘子。

步骤7：根据第h-1次迭代的计算结果计算出第h次迭代时各个三角形子域对应的相关中间变量，并且计算和集成第h次迭代时考虑薄板应变强化效应的整体刚度矩阵和整体变化载荷列阵。

步骤8：引入位移边界条件修改第h次迭代时的整体刚度矩阵和整体变化载荷列阵，求解相应的线性方程，依次得到第h次迭代时残余位移增量的中间变量、拉格朗日乘子、残余位移增量、每个角点载荷作用下各个三角形子域对应的光滑广义塑性应变、考虑薄板应变强化效应的安定上限载荷乘子。

步骤9：根据设定的误差容限和收敛条件，判断是否终止迭代。

步骤10：计算终止迭代时各个三角形子域对应的塑性耗散功率，并且根据各个三角形子域与Voronoi子域的面积关系，得到薄板结构各个节点的塑性耗散功率。

表1给出了采用本具体实施例(基于C¹节点自然单元法，C¹-nodal-NEM，采用三角形子域稳定相容节点积分)求解获得的安定上限载荷乘子以及计算时间，也给出了基于统一强度理论得到的理论解[1]、基于C¹自然单元法(C¹-NEM，采用传统的背景网格积分)求解获得的数值解[2]及计算时间。

表1夹支圆形薄板的安定上限载荷乘子和计算时间

从表1可以看出：

1、在理想弹塑性模型时时，使用本具体实施例得到的安定上限数值解，均分别与基于统一强度理论的理论解[1]、基于C¹-NEM的数值解吻合良好，这说明本具体实施例具有良好的计算精度。

2、在有限随动强化模型时，夹支圆形薄板在-q^max≤q≤q^max工况为交变塑性失效，使用本具体实施例和C¹-NEM得到的安定载荷乘子分别与理想弹塑性时的相同；在0≤q≤q^max工况为累积塑性破坏，使用本具体实施例和C¹-NEM得到的安定载荷乘子比理想弹塑性时的分别增加了22.371％和19.861％，这验证了忽略材料的强化效应会使得到的安定载荷偏于保守的结论。

3、如图6所示，整体而言(除个别情况外，如材料模型工况0≤q≤q^max时)，采用本具体实施例完成计算的迭代步数通常较少，其迭代计算时间大约只是采用C¹-NEM所花计算时间的70％左右。这说明本具体实施例采用的基于三角形子域稳定相容节点积分的C¹-nodal-NEM，相对于采用传统背景网格积分的C¹-NEM，通常具有更高的迭代计算效率。

图6给出了获得表1中夹支圆形薄板安定上限载荷乘子的迭代收敛曲线图，可见本具体实施例经过11～78次迭代，就能单调、稳定地收敛到所需要的计算精度。

图7～9示出了采用本具体实施例后处理得到的夹支圆形薄板安定极限状态的塑性耗散功率云图，真实、直观地反映了夹支圆形薄板的塑性安定失效模式。

综上，本发明采用性能优良的C¹节点自然单元法，解决了考虑薄板应变强化效应安定上限分析极小化迭代格式的建立与线性化求解问题，能够精确、高效和稳定地求解考虑薄板应变强化效应的安定上限载荷，具有格式较为简单、通用性强、易于程序实现、计算精度和效率高、数值稳定性好、后处理方便的技术效果。

Claims (6)

1.一种考虑薄板应变强化效应的塑性安定上限载荷计算方法，其特征在于，包括如下步骤：

ST1、准备计算数据，得到薄板结构各个三角形子域对应的光滑广义应变-位移速度关系矩阵；

ST2、计算广义弹性应力场，分别得到角点载荷和恒定载荷作用时薄板结构各个三角形子域对应的光滑广义弹性应力场；

ST3、在初始迭代时(h＝0)，假设整个薄板结构处于非屈服状态，求解线性方程组，依次得到初始迭代时的拉格朗日乘子、残余位移增量、每个角点载荷作用时各个三角形子域对应的光滑广义塑性应变、考虑薄板应变强化效应的安定上限载荷乘子；

ST4、在第h(h≥1)次迭代时，根据第h-1次迭代的计算结果，集成和求解相应的线性方程组，依次得到第h次迭代时的拉格朗日乘子、残余位移增量、每个角点载荷作用时各个三角形子域对应的光滑广义塑性应变、考虑薄板应变强化效应的安定上限载荷乘子，根据迭代收敛条件判断是否终止迭代计算；

ST5、计算结果后处理，得到薄板结构各个节点的塑性耗散功率。

2.根据权利要求1所述的塑性安定上限载荷计算方法，其特征在于，所述准备计算数据步骤具体包括：

ST1.1、准备薄板结构的离散节点、Delaunay三角形、位移边界节点、载荷(恒定载荷、变化载荷)、几何参数、材料参数(杨氏模量E、泊松比v、屈服应力Y、强化参数)、偏移系数γ、误差容限vol1和vol2；

ST1.2、根据薄板的节点、Delaunay三角形和位移边界节点的信息，采用C¹节点自然单元法依次确定围绕每个节点x_r(r＝1NP，NP为离散节点的总数)Voronoi子域S_r的顶点及其逆时针排列顺序，Voronoi结构顶点微小偏移后的坐标，每个三角形子域S_rs(RS为子域S_r的三角形子域总数)的面积A_rs、每条边的切向矢量n₁和法向矢量n₂、每条边上每个积分点的自然临近节点I(I＝1np，np为采用空圆准则确定的每个积分点的自然临近节点总数)及其C¹自然临近形函数的一阶导数/> x_rs为对应于三角形子域S_rs的虚拟节点，上标w、θ_x、θ_y分别是薄板的挠度、x方向的旋度、y方向的旋度，下标x、y分别表示对x方向和y方向求一阶导数，进而计算每个积分点自然临近节点对应的光滑广义应变-位移速度关系矩阵/>

式中，Γ_rs为子域S_rs的边界；

ST1.3、“对号入座”集成各个三角形子域S_rs全部积分点自然临近节点对应的矩阵(J＝1/>为合并全部积分点相同自然临近节点后的节点总数)，得到各个三角形子域S_rs对应的光滑广义应变-位移速度关系矩阵/>

3.根据权利要求2所述的塑性安定上限载荷计算方法，其特征在于，所述计算广义弹性应力场步骤具体包括：

ST2.1、根据各个三角形子域对应的光滑广义应变-位移速度关系矩阵和弹性本构关系，计算和集成薄板的整体弹性刚度矩阵K^e：

式中，D_b为弹性关系矩阵；

为薄板的弯曲刚度；

h₀为板厚；

ST2.2、根据作用在薄板上每组独立变化载荷的范围和组数确定角点载荷和其总数对于每个角点载荷中的均布力，按照Voronoi子域S_r的面积/>将其等效为作用在节点上的等效载荷列阵/>对于每个角点载荷中的集中力，将其直接施加在所作用的节点上；由各个节点的等效载荷列阵/>“对号入座”集成每个角点载荷作用时薄板的整体弹性载荷列阵/>

ST2.3、引入位移边界条件修改薄板的整体弹性刚度矩阵K^e和每个角点载荷作用时的整体弹性载荷列阵求解薄板弹性问题的线性控制方程/>得到每个角点载荷作用时薄板的弹性位移场/>

ST2.4、根据弹性本构关系和弹性位移场，计算每个角点载荷作用时各个三角形子域S_rs对应的光滑广义弹性应力

ST2.5、计算恒定载荷作用时各个三角形子域S_rs对应的光滑广义弹性应力

4.根据权利要求3所述的塑性安定上限载荷计算方法，其特征在于，所述在初始迭代时，假设整个薄板结构处于非屈服状态，求解线性方程组步骤具体包括：

ST3.1、假设整个薄板结构处于非屈服状态，取和/>依次计算出初始迭代时各个三角形子域S_rs对应的中间变量/>

式中，M_P＝Yh₀ ²/4为塑性极限弯矩；

Q是正定对称的常数矩阵，且Q^-1＝2D/3，D为处理薄板塑性不可压条件引入的常数矩阵：

ST3.2、计算和集成整体恒定载荷列阵F₀、初始迭代时考虑薄板应变强化效应的整体刚度矩阵K⁰和整体变化载荷列阵F₁ ⁰：

ST3.3、引入位移边界条件修改整体恒定载荷列阵F₀、初始迭代时考虑薄板应变强化效应的整体刚度矩阵K⁰和整体变化载荷列阵F₁ ⁰，分别求解线性方程K⁰(Δa₀)⁰＝F₀和K⁰(Δa₁)⁰＝F₁ ⁰，得到初始迭代时残余位移增量的中间变量(Δa₀)⁰和(Δa₁)⁰；

ST3.4、求解初始迭代时的拉格朗日乘子λ⁰：

式中，中间变量和/>分别为：

ST3.5、求解初始迭代时的残余位移增量(Δa)⁰和在每个角点载荷作用下各个三角形子域S_rs对应的光滑广义塑性应变

(Δa)⁰＝(Δa₀)⁰+λ⁰(Δa₁)⁰

ST3.6、根据得到的残余位移增量和光滑广义塑性应变，依次计算出初始迭代时各个三角形子域S_rs对应的变量和μ¹(x_rs)、/>

式中，β₁和β₂均为远小于1的正数，且分别设定为：

进而得到初始迭代时考虑薄板应变强化效应的安定上限载荷乘子s⁰：

5.根据权利要求4所述的塑性安定上限载荷计算方法，其特征在于，所述在第h(h≥1)次迭代时，根据第h-1次迭代的计算结果，集成和求解相应的线性方程组步骤具体包括：

ST4.1、根据第h-1次迭代的计算结果确定出变量和/>的值，依次计算出第h次迭代时各个三角形子域S_rs对应的中间变量/>

ST4.2、计算和集成第h次迭代时考虑薄板应变强化效应的整体刚度矩阵K^h和整体变化载荷列阵F₁ ^h：

ST4.3、引入位移边界条件修改整体恒定载荷列阵F₀、第h次迭代时的整体刚度矩阵K^h和整体变化载荷列阵F₁ ^h，分别求解线性方程K^h(Δa₀)^h＝F₀和K^h(Δa₁)^h＝F₁ ^h，得到第h次迭代时残余位移增量的中间变量(Δa₀)^h和(Δa₁)^h；

ST4.4、求解第h次迭代时的拉格朗日乘子λ^h：

式中，中间变量和/>分别为：

ST4.5、求解第h次迭代时的残余位移增量(Δa)^h和在每个角点载荷作用下各个三角形子域S_rs对应的光滑广义塑性应变

(Δa)^h＝(Δa₀)^h+λ^h(Δa₁)^h

ST4.6、根据得到的残余位移增量和光滑广义塑性应变，依次计算出第h次迭代时各个三角形子域S_rs对应的变量和μ^h+1(x_rs)、/>

从而计算出第h次迭代时考虑薄板应变强化效应的安定上限载荷乘子s^h：

ST4.7、根据设定的误差容限vol1、vol2和下列收敛条件，判断是否终止迭代：

||(Δa)^h-(Δa)^h-1||/||(Δa)^h-1||≤vol1，|s^h-s^h-1|/s^h-1≤vol2。

6.根据权利要求5所述的塑性安定上限载荷计算方法，其特征在于，所述计算结果后处理步骤具体包括：

ST5.1、假设第h次迭代终止时满足收敛条件，计算终止迭代时各个三角形子域S_rs对应的光滑塑性耗散功率D^h(x_rs)：

ST5.2、根据各个三角形子域S_rs与Voronoi子域S_r的面积关系，得到薄板结构各个节点的塑性耗散功率D^h(x_r)：