薄壳在包括土木建筑、航天航空、机械电子等在内的工程科技领域有着广泛应用，如穹顶、航天器、舰艇以及工业设备与产品的外壳等. 壳体造型多为复杂曲面，受力后将同时发生面内位移和横向位移，整体的薄膜状态和弯曲状态相互耦合，需要同时分析，这使得薄壳的力学求解通常比较复杂. 早期均通过建立微分方程来获得其解析解^[1]，但面对非线性问题时就非常困难，因而目前主要是根据虚功或能量变分原理建立弱形式的等效积分方程，并借助数值方法来求出一组离散解. 但面对实际工程中不断出现的一些复杂动力问题，尤其是同时涉及到壳体结构的大变形、大转动、屈曲、碰撞、破坏等高度非线性的力学行为时，许多传统方法常由于计算量和数值稳定等方面的局限性而求解困难^[2-3]. 因此，亟需发展力学描述与计算列式更为简洁且能够有效处理薄壳动力强非线性问题的数值分析方法.

有限单元法是当前薄壳非线性分析中最常用的数值方法. 依据不同的几何描述、壳体理论、插值函数及各异的处理技巧，已开发了诸多有限元列式，大致可将其分为两类^[4]：一类是直接建立于壳体理论平衡方程的等效积分弱形式，如DKT单元^[5]、BWC单元^[6]、BST单元^[7]、ANCF单元^[8]、IGAS单元^[9]等；另一类是通过在连续体单元上强制引入壳体理论的运动假设建立壳单元（即基于连续体（CB）的壳），如MITC4单元^[10]、HL单元^[11]等. 对于前者，当涉及到非线性问题时，由于理论上的限制，要建立与壳单元变分弱形式等效的控制方程会有本质上的困难；而后者通过CB方法建立壳单元离散方程，比较直观，适用于任意壳体的大变形问题，但由于采用的是等效R-M运动假设，可能出现剪切和薄膜自锁现象，需要引入设置缩减积分或假设应变场等特殊技巧才能够完成分析^[12].

有别于传统分析力学架构下的计算方法，有限质点法是基于向量式力学理论提出的面向工程结构复杂行为的新型数值分析方法^[13]. 该方法采取多质点离散的物理分析模式，根据牛顿定律通过追踪质点的运动轨迹来获得结构整体的运动和变形情况，动力和几何非线性描述是该方法的本质特征；各质点运动方程相互独立，不涉及刚度矩阵的集成和求逆，也不需要迭代求解，避免了方程病态或刚度矩阵奇异而造成的数值求解困难，在结构复杂非线性问题分析中具有优势和灵活性. 目前该方法已被应用于机构运动、接触碰撞、结构屈曲、倒塌破坏、固体大变形、膜结构找形、结构多尺度精细化等结构复杂行为的分析^[14-15]，均取得了较好的效果. 同时，有许多学者同样基于向量式力学和数值计算相结合的概念发展了向量式有限元方法^[16]，开发了包括DKT壳单元在内的众多类型单元，对薄壳的弹塑性本构、碰撞和穿透等问题进行了初步研究^[17-18].

本文将通过理论推导、程序编制及算例，进一步发展通用的有限质点法薄壳动力非线性行为的分析框架，完善此类方法的计算理论体系. 首先介绍薄壳的有限质点法基本理论和控制方程，推导有限质点法处理薄壳纯变形和质点内力（矩）计算公式；然后根据有限质点法的计算特点对薄壳质点的转动惯量、时间步长和阻尼等关键参数取值及质点控制方程的数值求解提出合理可行的处理策略；同时利用经典的J₂流动理论处理薄壳平面应力弹塑性问题，建立算法的具体实现步骤；最后通过算例的模拟和分析，检验方法的正确性和有效性.

有限质点法采用标准化向量式力学分析架构，可同时求解结构体系受力后的位置变化与几何变形，其基本概念包括：1）点值描述；2）途径单元；3）刚体运动和纯变形分离的物理模式处理；4）基于牛顿第二定律的强形式广义质点运动控制方程. 关于有限质点法的基本思想详见文献[13]，这里结合薄壳理论仅简单介绍其要点.

薄壳属于退化三维固体模型，应同时考虑薄膜刚度和弯曲刚度的影响. 有限质点法将薄壳离散为一组具有空间平动与转动自由度的质点集合，质点间由一系列具备描述面内外变形能力的薄壳元相连. 薄壳中所有物理量均由各个质点的点值描述. 薄壳元没有质量，仅用于表示质点间相互作用，随质点运动发生变形，由此产生的内力反作用于相连质点.

如图1所示，假设一薄壳构件经历一个连续变形和运动过程，从初始时刻t₀的构形S₀经中间时刻（ ${t_a} = {t_k} - \tau $）的构形S_a至当前时刻t_k的构形S_k. 该时段内（ ${t_0} - {t_k}$）位于薄壳离散模型中面上的任一质点j的位置从 ${{{x}}_j}_0$ 经 ${{{x}}_{ja}}$ 运动至 ${{{x}}_{jk}}$，其轨迹被一组时间点 $\left( {{t_0},{t_1},{t_2}, \cdots ,{t_a},{t_b}, \cdots ,{t_k}, \cdots } \right)$ 分割成若干个独立的时间片段，每个时段内的运动轨迹即为1个途径单元. 通过追踪质点轨迹可获得薄壳整体的运动变形情况，从而将问题转化为求解壳面上的质点集合在途径单元内的坐标、转角、速度等物理量的变化情况.

有限质点法从物理角度考虑，认为所有质点均处于动平衡状态，在每个途径单元（如t_a≤t≤t_b）内的运动都遵循牛顿第二定律，且控制方程独立. 薄壳中的质点运动是广义的，对于壳面上任一质点j，其运动变量包括沿坐标轴方向的3个平动自由度和绕轴旋转的3个转动自由度，对应3个坐标方向的质点力和力矩，因此薄壳质点运动控制方程式可分为点的平动公式和点的转角公式. 前者用于求解中面内变形引起的质点位置改变（可直接利用薄膜质点空间运动方程^[19]进行求解），后者用于求解面外弯曲变形引起的质点截面转动.

根据经典薄壳理论，质点截面转动只产生弯曲变形而没有壳面内扭转，因此需要转换到质点的切平面坐标系下建立其转动控制方程. 如图2所示，令坐标平面 $\tilde x\tilde o\tilde y$（坐标轴单位向量为 $\tilde {{e}}_x^j$ 和 $\tilde {{e}}_y^j$）位于过壳面上任意一点j的切平面内， $\tilde {\rm{z}} $轴单位向量 $\tilde {{e}}_z^j =\left(\tilde {{e}}_x^j \times \tilde {{e}}_y^j\right)\Big/\Big| \tilde {{e}}_x^j \times \tilde {{e}}_y^j \Big|$. 切平面坐标系（ $\tilde x,\tilde y,\tilde z$）与全域坐标系（ $x,x,z$）的转换矩阵为 $\tilde {{\varOmega }}$，则切平面坐标下的质点转动向量为 $ {\tilde {{\beta }}^j} = $ ${\left[ {\begin{array}{*{20}{c}} {\tilde \beta _x^j},\ {\tilde \beta _y^j},\ {\tilde \beta _z^j} \end{array}} \right]^{\rm{T}}} = \tilde {{\varOmega }}\;{{{\beta }}^j}$，其中独立转动自由度只有2个绕 $\tilde x$ 和 $\tilde y$ 轴方向的，可按下列运动方程式计算：

(1) $ {\tilde {{I}}^*}_j {\ddot {\tilde {{\beta }}}}_j^* = {\tilde {{M}}}_j^*,\;\;{\tilde {{I}}^*}_j = {\left[\!\!\! {\begin{array}{*{20}{c}} {{{\tilde I}_{xx}}}&{{{\tilde I}_{xy}}} \\ {{{\tilde I}_{yx}}}&{{{\tilde I}_{yy}}} \end{array}} \!\!\!\right]_j},\;\;\ddot {\tilde {{\beta }}}_j^* = {\left\{\!\!\! {\begin{array}{*{20}{c}} {{{\ddot {\tilde \beta} }_x}} \\ {{{\ddot {\tilde \beta} }_y}} \end{array}} \!\!\!\right\}_j},\;\;\tilde {{M}}_j^* = {\left\{\!\!\! {\begin{array}{*{20}{c}} {{{\tilde M}_x}} \\ {{{\tilde M}_y}} \end{array}} \!\!\!\right\}_j}. $

式中： ${\tilde I_{xx}}$、 ${\tilde I_{xy}}$、 ${\tilde I_{yy}}$ 为质点j切平面坐标系下质量惯性矩阵元素； ${\tilde M_x}$、 ${\tilde M_y}$ 为切平面坐标下质点j的合力矩分量，可由下式求得：

${\tilde {{M}}_j} = {\tilde {{\varOmega }}}{{{M}}_j} = {\tilde {{\varOmega }}}\left( {{{M}}_j^{\rm{ext}} + {{M}}_j^{\rm{int}} + {{M}}_j^{\rm{dmp}}} \right).$

其中， ${{{M}}_j}$ 为全域坐标下质点j的合力矩， ${{M}}_j^{\rm{ext}}$、 ${{M}}_j^{\rm{int}}$、 ${{M}}_j^{\rm{dmp}}$ 分别为质点j的外力矩向量、内力矩向量、阻尼力矩向量. 绕 $\tilde x$ 和 $\tilde y$ 轴的转角是耦合的，因而2个方向的转动方程需要联立求解.

质点内力矩来自于质点转动引起的壳元弯曲变形，可通过其相连壳元的节点等效内力矩 ${{m}}_\alpha ^{\rm{int}}$ 组集得到，即

${{M}}_j^{\rm{int}} = - \sum\nolimits_{\alpha = 1}^{{n_{\rm c}}} \ {{{m}}_\alpha ^{\rm{int}}}, $

式中：n_c为质点j相连薄壳元个数.

质点外力矩包括作用在质点上的集中力矩 ${{m}}_j^{\rm{ext}}$ 以及由相连壳元上分布荷载转化而来的等效力矩 ${{m}}_\alpha ^{\rm{ext}}$，即

${{M}}_j^{\rm{ext}} = {{m}}_j^{\rm{ext}} + \sum\nolimits_{\alpha = 1}^{{n_{\rm c}}}\ {{{m}}_\alpha ^{\rm{ext}}} .$

阻尼力矩可写成瑞利比例阻尼的形式根据壳体黏滞特性计算：

${{M}}_j^{\rm{dmp}} = {{m}}_j^{\rm{mass}} + {{m}}_j^{\rm{stiff}}.$

其中，质量阻尼力矩

${{m}}_j^{\rm{mass}} = - {a_{\rm{m}}}{{{I}}_j}{\dot {{\beta }}_j},$

I_j为质点j在全域坐标系下的惯性矩阵；刚度阻尼力矩

${{m}}_j^{\rm{stiff}} = - {a_{\rm{s}}}\sum\nolimits_{\alpha = 1}^{{n_{\rm c}}}\ {{{m}}_\alpha ^{\rm{stiff}}} \left( {{{\dot {{\beta }}}_j}} \right),$

其中，a_m和a_s为阻尼参数， ${\dot {{\beta }}_j}$ 为质点j的截面转速， ${{m}}_\alpha ^{\rm{stiff}}\left( {{{\dot {{\beta }}}_j}} \right)$ 为通过壳元内黏性应力等效积分求得的节点阻尼力矩.

绕 $\tilde z$ 轴旋转角 $\tilde \beta _{jz}$ 被视作刚体转动，通过质点j与相邻质点间的相对位移进行估算，即

(2) $ \tilde \beta _{jz}^{n + 1} = \tilde \beta _{jz}^n + \frac{1}{{{N_j}}}\sum\nolimits_{I = 1}^{{N_j}}\ {\Delta {{{\phi }}_I} \cdot \tilde {{e}}_z^j} , $

(3) $ \Delta {{{\phi }}_I} = {\sin ^{ - 1}}\left[ {\displaystyle\frac{{\left| {{ x}_n^{jI} \times { x}_{n + 1}^{jI}} \right|}}{{\left| {{ x}_n^{jI}} \right| \times \left| {{ x}_{n + 1}^{jI}} \right|}}} \right]\displaystyle\frac{{{ x}_n^{jI} \times { x}_{n + 1}^{jI}}}{{\left| {{ x}_n^{jI} \times { x}_{n + 1}^{jI}} \right|}}. $

式中：N_j为质点j相邻质点数； ${{x}}_n^{jI}$ 和 ${{x}}_{n + 1}^{jI}$ 分别为t_n和t_n+1时刻点j与相邻点I的相对位置向量；法向量 $\tilde {{e}}_z^j$ 可通过相邻壳元法向量的面积A_α加权平均得到，即

$\tilde { e}_z^j{\rm{ = }}{{\sum\nolimits_{\alpha {\rm{ = }}1}^{{N_j}} {{A_\alpha }{{n}}_j^\alpha } } \Bigg/ {\sum\nolimits_{\alpha {\rm{ = }}1}^{{N_j}} {{A_\alpha }} }}.$

式中： $n_j^\alpha $为与质点j相连壳元α的单位法向量.

对应薄壳质点的平动和转动分量，壳元的变形、内力也可分为两部分，分别与面内薄膜刚度和面外弯曲刚度相关. 构造一种通用的三角形薄壳元来模拟壳面质点间的相互作用，其内力由薄壳中面内的薄膜应力和沿厚度变化的面外弯扭应力累加构成，壳元内部自身不考虑2种变形之间的耦合. 前者可借助T3薄膜元直接求解，相关公式参见文献[19]，这里不再赘述；后者可通过引入BCIZ形函数^[20]构建一种弯曲元进行计算，下面给出主要推导过程.

如图3所示，三节点薄壳元（厚度为h）的空间位置和几何构形由位于中面内的3个节点（编号i=1，2，3）确定，每个节点包含3个平动自由度和3个转动自由度. 节点（i=1，2，3）与质点（I=p，q，r）刚性联结，t_a和t时刻节点位置向量和转动向量分别记作 $({{x}}_i^{\rm a},{{\beta }}_i^{\rm a}) = ({{x}}_I^{\rm a},{{\beta }}_I^{\rm a})$ 和 $({{{x}}_i},{{{\beta }}_i}) = ({{{x}}_I},{{{\beta }}_I})$. 以t_a时刻薄壳元构形 ${V_{\rm a}}({1_{\rm a}} - {2_{\rm a}} - {3_{\rm a}})$ 为基础架构，t_a~t时段内节点线位移和转角增量分别为

$\Delta {{{u}}_i} = {{{x}}_i} - {{x}}_i^{\rm a},\,\, \Delta {{{\beta }}_i} = {{{\beta }}_i} - {{\beta }}_i^{\rm a},\;\;i = 1,2,3.$

壳面质点内力（矩）仅与相连壳元的纯变形相关，所涉及到的应变、应力物理量均通过节点的纯变形量计算，因此需要对Δu_i和 $\Delta {{{\beta }} _i}$ 中的刚体运动分量进行扣除，这是有限质点法求解大变形问题的关键. 参照T3薄膜元中基于“虚拟运动”的物理模式估算方法，求解薄壳元节点纯变形量主要步骤如下：1）估算t_a~t时段内薄壳元的刚体平移Δu₁和刚体转动 ${{\theta }}$（分成面外转动向量 ${{{\theta }}_{\rm{op}}} = {\theta _{\rm{op}}}{\bar {{n}}_{\rm{op}}}$ 和面内转动向量 ${{{\theta }}_{\rm{ip}}} = {\theta _{\rm{ip}}}{\bar {{n}}_{\rm{ip}}}$）；2）以虚拟逆向平移（−Δu₁）和逆向转动（− ${{\theta }}$）方式扣除刚体运动，估算t_a−t时段内节点变形位移 ${{\eta }}_i^{\rm{d}}$ 和变形转角 ${{\beta }}_i^{\rm{d}}$ 为

(4) $ {{\eta }}_1^{\rm{d}} = 0; $

(5) $ {{\eta }}_i^{\rm{d}} = {{{\eta }}_i} + {{\eta }}_i^{\rm{r}} = \left( {\Delta {{{u}}_i} - \Delta {{{u}}_1}} \right) + \left[ {{{R}}\left( { - {{\theta }}} \right) - {{I}}} \right]\left( {{{{x}}_i} - {{{x}}_1}} \right); $

(6) $ {{\beta }}_i^{\rm{d}} = \Delta {{{\beta }}_i} + \left( { - {{{\theta }}_{\rm{op}}}} \right) + \left( { - {{{\theta }}_{\rm{ip}}}} \right). $

式中：η_i和 ${{\eta }}_i^{\rm{r}}$ 分别为节点i相对于参考点（节点1）的全位移和逆向刚体转动位移；R为转动矩阵.

与传统有限元法不同，式（4）~（6）求出的节点变量是变形位移而非全位移. 内插函数包含刚体模式，将其用于描述薄壳元内的变形分布时，需要将对应于刚体运动模式自由度的3个多余节点变量舍去. 薄壳元经历逆向刚体运动后，处于虚拟位置上的各个节点与t_a时刻基础架构位置共面，因此与弯曲元刚体模式相关的3个面外位移分量已被完全扣除，剩下的节点变形转角分量均为独立变量. 为了与膜元部分合算内力，这里仍参照T3薄膜元的处理方式，在壳元基础架构上定义一组变形坐标 $\hat {{x}} = {\left[ {\hat x,\;\hat y,\;\hat z} \right]^{\rm{T}}}$，基底向量为 $({\hat {{e}}_x},{\hat {{e}}_y},{\hat {{e}}_z})$，坐标原点与坐标轴方向设置方式均与薄膜元相同.

令 $\hat {{Q}} = {\left[ {\hat {{e}}_x ,\; \hat {{e}}_y ,\; \hat {{e}}_z} \right]^{\rm{T}}}$ 为变形坐标与全域坐标间的转换矩阵，则节点线位移和转角位移用变形坐标分量可表示为

$\hat {{u}}_i = {\left[ {{{\hat u}_i}{\rm{ }},\;{{\hat v}_i}{\rm{ }},\;{{\hat w}_i}} \right]^{\rm{T}}} = \hat {{Q}}{{\eta }}_i^{\rm{d}},$

${\hat {{\varphi }}_i} = {\left[ {\begin{aligned} {\hat \varphi _{xi}},\;{\hat \varphi _{yi}},\;{\hat \varphi _{zi}} \end{aligned}} \right]^{\rm{T}}} = \hat {{Q}}{{\beta }}_i^{\rm{d}}.$

对于 ${{\beta }}_i^{\rm{d}}$，因刚体转动θ_ip和θ_op满足 ${\bar {{n}}_{\rm{ip}}} = {\hat {{e}}_z} = {{{n}}_a}$ 及 ${{{\theta }}_{{\rm{op}}}} \bot {{{n}}_a}$，有

${\hat \theta _z} = {{\theta }} \cdot {\hat {{e}}_z} = {\theta _{\rm{op}}}\left( {{{\bar {{n}}}_{\rm{op}}} \cdot {{\hat {{e}}}_z}} \right) + {\theta _{\rm{ip}}} = {\theta _{\rm{ip}}},$

另外2个变形坐标方向上的刚体转动分量则由θ_op转换而来，即

$ {\hat \theta _x} = {{\theta }} \cdot {\hat {{e}}_x} = {\theta _{\rm{op}}}\left( {{{\bar {{n}}}_{\rm{op}}} \cdot {{\hat {{e}}}_x}} \right), \; {\hat \theta _y} = {{\theta }} \cdot {\hat e_y} = {\theta _{\rm{op}}}\left( {{{\bar {{n}}}_{\rm{op}}} \cdot {{\hat {{e}}}_y}} \right). $

于是，式（7）中节点变形转角分量可具体表示如下：

(7) $ {{\hat{ \varphi }}_i} = \left\{ \begin{array}{l} \mathop {\hat \varphi _{xi}}\limits \\ \hat \varphi _{yi}\\ \mathop {\hat \varphi _{zi}}\limits \end{array} \right\} = \left\{ \begin{aligned} & \Delta \hat \beta _{xi} + \left( { - {{\hat \theta }_x}} \right)\\ & \Delta \beta _{yi} + \left( { - {{\hat \theta }_y}} \right)\\ & \Delta \beta _{zi} + \left( { - {{\hat \theta }_z}} \right) \end{aligned} \right\}. $

遵循经典板壳理论不考虑面内扭转刚度的基本假设，计算中可忽略 $\hat \varphi _{zi}$. 再根据壳元面内与面外变形特点，将剩余非零独立节点变量相应分成两部分.

1) $\hat x$ 和 $\hat y$ 轴方向线位移（计算薄膜内力）：

$\hat {{u}}_{\rm{m}}^* = {\left[ {\begin{array}{*{20}{c}} {{{\hat u}_{2x}}},\;{{{\hat u}_{3x}}},\;{{{\hat u}_{3y}}} \end{array}} \right]^{\rm{T}}};$

2) 绕 $\hat x$ 和 $\hat y$ 轴转角（计算弯扭内力（矩））：

$\hat {{u}}_{\rm b}^* = {\left[ {\begin{array}{*{20}{c}} {\hat \varphi _{x1}},\;{\hat \varphi _{y1}},\;{\hat \varphi _{x2}},\;{\hat \varphi _{y2}},\;{\hat \varphi _{x3}},\;{\hat \varphi _{y3}} \end{array}} \right]^{\rm{T}}}.$

为了描述壳元内的变形和应力，引入传统有限元中的BCIZ平板元形函数作为内插函数，但节点变量仅考虑6个独立转角量 $\hat {{u}}_{{\rm b}i}^* = {\left[ {\begin{array}{*{20}{c}} {\hat \varphi _{xi}},\;{\hat \varphi _{yi}} \end{array}} \right]^{\rm{T}}}$，则t_a~t时段内壳元中面内的变形挠度分布模式经修正后可表示为

(8) $ \hat w\left( {\hat x,\hat y} \right) = \hat {{N}}_1^*\hat {{u}}_{{\rm{b}}1}^* + \hat {{N}}_2^*\hat {{u}}_{{\rm{b}}2}^* + \hat {{N}}_3^*\hat {{u}}_{{\rm{b}}3}^* = {\hat {{N}}^*}\hat {{u}}_{\rm{b}}^*. $

式中： ${\hat {{N}}^*} = \left[ {\begin{aligned} {\hat {{N}}_1^*},\;{\hat {{N}}_2^*},\;{\hat {{N}}_3^*} \end{aligned}} \right] $，是以三角形面积坐标 $({\hat L_i}\left( {\hat x,\hat y} \right))$表示的修正后的形函数矩阵，由原形函数矩阵删去对应于法向位移 ${\hat w_i}$的元素后得到，即 $\hat {{N}}_i^* = \left[ {\begin{aligned} {\hat N_{xi}},\;{\hat N_{yi}} \end{aligned}} \right]$.

根据途径单元内小变形假设，应力和应变可采用工程应力和微应变形式，无需引用复杂应力率和应变率公式，本构关系可由材料试验直接给出. 利用 $\hat w\left( {\hat x,\hat y} \right)$ 可得t_a~t时段内的弯曲变形应变增量为

(9) $ \Delta \hat {{\varepsilon }}_{\rm{r}}^{\rm{b}} = \hat z\;{{\kappa }} = \hat z\hat {{B}}_{\rm b}^*\hat {{u}}_{\rm b}^* = \hat z\left[ {\begin{aligned} {\hat {{B}}_{{\rm b}1}^*},\;{\hat {{B}}_{{\rm b}2}^*},\;{\hat {{B}}_{{\rm b}3}^*} \end{aligned}} \right]\hat {{u}}_{\rm b}^*. $

式中： $\Delta \hat {{\varepsilon }}_{\rm{r}}^{\rm{b}} = \left[ {\begin{aligned} {\Delta \hat \varepsilon _x^{\rm{b}}},\;{\Delta \hat \varepsilon _y^{\rm{b}}},\;{\Delta \hat \gamma _{xy}^{\rm{b}}} \end{aligned}} \right]^{\rm{T}}; $ ${{\kappa }} = - \Big[ {\begin{aligned} {\hat w{,_{xx}}}\,\,\,\,{\hat w{,_{yy}}}\end{aligned}} \, $ $\,\,\,\,{2\hat w{,_{xy}}}\Big]^{\rm{T}} , $为曲率向量； $\hat {{B}}_{\rm b}^*$ 是将BCIZ形函数求导后的变形矩阵考虑 ${\hat w_1} = {\hat w_2} = {\hat w_3} = 0$ 后修正得到的. 此外，薄膜变形引起的应变增量为 $\Delta \hat {{\varepsilon }}_{\rm{r}}^{\rm{m}} = \hat {{B}}_{\rm{m}}^*\hat {{u}}_{\rm{m}}^*$. 于是，以t_a时刻构形为基础，当前任一点的总应变增量与合应力为

(10) $ \Delta \hat {{\varepsilon }}_{\rm{r}} = \Delta \hat {{\varepsilon }}_{\rm{r}}^{\rm{m}} + \Delta \hat {{\varepsilon }}_{\rm{r}}^{\rm{b}}; $

(11) $ \hat {{\sigma }}_{\rm{r}} = \hat {{\sigma }}_{\rm a} + {\hat {{D}}_{\rm a}}\Delta \hat {{\varepsilon }}_{\rm{r}} = \hat {{\sigma }}_{\rm a} + \Delta \hat {{\sigma }}_{\rm{r}}^{\rm{m}} + \Delta \hat {{\sigma }}_{\rm{r}}^{\rm{b}}. $

式中： ${\hat {{D}}_{\rm a}}$ 为应力等于 $\hat {{\sigma }}_{\rm a}$ 时沿变形坐标轴 $\hat x$ 和 $\hat y$ 方向的材料本构矩阵. 若材料进入塑性，则需要根据平面应力弹塑性增量理论对 $\hat {{\sigma }}_{\rm{r}}$ 进行修正.

不考虑 $\Delta \hat {{\varepsilon }}_{\rm{r}}^{\rm{m}}$ 与 $\Delta \hat {{\varepsilon }}_{\rm{r}}^{\rm{b}}$ 之间的耦合作用，由变形体的虚功原理可得t时刻的等效节点内力和弯矩分别为

(12) $ \hat {{f}}^* = \int_{{A_{\rm{a}}}} {\int_{{{ - {h_{\rm{a}}}} / 2}}^{{{{h_{\rm{a}}}} / 2}} {{{\left( {\hat {{B}}_{\rm{m}}^*} \right)}^{\rm{T}}}} } {\hat {{\sigma }}_{\rm{r}}}{\rm d}\hat z{\rm d}\hat A, $

(13) $ {{\hat{ m}}^*} = \int_{{A_{\rm{a}}}} {\int_{{{ - {h_{\rm{a}}}} / 2}}^{{{{h_{\rm{a}}}} / 2}} {\hat z{{\left( {{\hat{ B}}_{\rm{b}}^*} \right)}^{\rm T}}} } {{\hat{ \sigma }}_{\rm{r}}}{\rm d}\hat z{\rm d}\hat A. $

式中：A_a为t_a时刻的薄壳元面积，h_a为t_a时刻的薄壳元厚度.

虽然这里内力求解过程中不生成刚度矩阵而无法直接计算刚度阻尼，但考虑到刚度阻尼本质上源于材料黏性效应^[21]，因此阻尼力（矩）也可作为附加内力（矩）通过等效积分来计算得到，即

(14) $ \hat {{f}}^{\rm{stiff}} = \frac{1}{{\Delta {t_{\rm{r}}}}}\int_{{A_{\rm{a}}}} {\int_{{{ - {h_{\rm{a}}}} / 2}}^{{{{h_{\rm{a}}}} / 2}} {{{\left( {\hat {{B}}_{\rm{m}}^*} \right)}^{\rm{T}}}\left( {{{\hat \sigma }_{\rm{r}}} - {{\hat \sigma }_{\rm{a}}}} \right)} } \;{\rm d}\hat z{\rm d}\hat A; $

(15) $ \hat {{m}}^{\rm{stiff}} = \frac{1}{{\Delta {t_{\rm{r}}}}}\int_{{A_{\rm{a}}}} {\int_{{{ - {h_{\rm{a}}}} / 2}}^{{{{h_{\rm{a}}}} / 2}} {\hat z{{\left( {\hat {{B}}_{\rm{b}}^*} \right)}^{\rm{T}}}} } \left( {{{\hat \sigma }_{\rm{r}}} - {{\hat \sigma }_{\rm{a}}}} \right){\rm d}\hat z{\rm d}\hat A. $

式（14）、（15）中应力、应变均随位置变化，直接积分较为不便，可借助数值方法进行求解. 对于弹性和弹塑性问题，可分别采用沿厚度方向的2点Gauss积分方案和5点Labatto积分方案，中面内则分别采用三角形3点和7点高斯积分方案.

式（12）和（13）只能求出部分节点力分量 $\left( {{{\hat f}_{2x}},{{\hat f}_{3x}},{{\hat f}_{3y}}} \right)$ 和弯矩分量 $\left( {{{\hat m}_{1x}},{{\hat m}_{1y}},{{\hat m}_{2x}},{{\hat m}_{2y}},{{\hat m}_{3x}},{{\hat m}_{3y}}} \right)$，其余6个分量对应于被扣除的刚体自由度，需通过壳元在平面内、外的静力平衡条件来确定，即

(16) $\left. \begin{aligned} & \displaystyle\sum {{M_{\hat z}}} = 0,\;\displaystyle\sum {{F_{\hat x}}} = 0,\;\displaystyle\sum {{F_{\hat y}}} = 0 \Rightarrow {{\hat f}_{1x}},{{\hat f}_{1y}},{{\hat f}_{2y}},\\ & \displaystyle\sum {{M_{\hat x}}} = 0,\;\displaystyle\sum {{M_{\hat y}}} = 0,\;\displaystyle\sum {{F_{\hat z}}} = 0 \Rightarrow {{\hat f}_{1z}},{{\hat f}_{2z}},{{\hat f}_{3z}}. \end{aligned}\right\} $

上述薄壳面内、外变形内力的等效节点力 ${\hat {{f}}_i} = $ ${\left[ {\begin{aligned} {{{\hat f}_{ix}}},\;{{{\hat f}_{iy}}},\;{{{\hat f}_{iz}}} \end{aligned}} \right]^{\rm{T}}}$ 和节点弯矩 ${\hat {{m}}_i} = {\left[ {\begin{aligned} {{{\hat m}_{ix}}},\;{{{\hat m}_{iy}}},\;{{{\hat m}_{iz}}} \end{aligned}} \right]^{\rm{T}}}$ 都是在虚拟位置的变形坐标系下推导的，还需要通过正向运动与坐标变换转到当前 $t$ 时刻的全局坐标系下，即

$ {{{f}}_i} = {{R}}\left( {{\theta }} \right){\hat {{Q}}^{\rm{T}}}\hat {{f}}_i, \; {{{m}}_i} = {{R}}\left( {{\theta }} \right){\hat {{Q}}^{\rm{T}}}\hat {{m}}_i, $

最后，按照薄壳元节点与质点的对应关系，将上述f_i和m_i反组集到节点i所连接的质点上.

如上文所述，薄壳所有质量均由质点承担，要求离散质点的质量分布收敛于实际情况. 与变形相对应，薄壳质点的广义质量有两部分，即点的平动质量M_j和转动惯量I_j. 前者的计算与薄膜质点质量的计算方法相同^[19]；后者把质点视为具有一定形状的均布质量截面，通过组集质点及其相连壳元均匀分配的等效转动惯量得到，即

${{I}}_j = {{I}}_j^0 + \sum\limits_{\alpha = 1}^{{n_{\rm c}}} \ {{{I}}_\alpha ^{\rm{e}}}. $

等效转动惯量 ${{I}}_\alpha ^{\rm{e}}$ 可根据d’Alembert原理通过计算惯性力矩的转动虚功得到. 以质点j相连节点1为例（见图4），节点等效质量为 $\overline m_1^{\rm{e}} = {{{m_{\rm{e}}}} / 3}$，m_e为薄壳元的质量，则截面上任意一点的转动位移可表示为

$ \hat {{u}}_{1}^{\rm{b}} = \left[\hat u_{1x}^{\rm{b}} ,\hat u_{1y}^{\rm{b}}\right]^{\rm T}= - z{\hat \varphi _{1}}. $

惯性力矩的转动虚功可表示为

$\delta W_1^I = -{\int_{A_{\rm a}^1} {\left( {\delta \hat {{u}}_1^{\rm{b}}} \right)} ^{\rm{T}}} {\ddot{ \hat {{u}}}_1^{\rm{b}}} d \overline m_1^{\rm{e}} = - {\left( {\delta {\ddot {\hat {{\varphi }}}_1}} \right)^{\rm{T}}}\hat {{I}}_1^*\hat {{\varphi }}_1^*$

式中： $ {\hat { I}} $为该点与中性面的距离， $ {\hat {{\varphi}} _1} = {\left[ {{{\hat \varphi }_{1x}},{{\hat \varphi }_{1y}}} \right]^{\rm{T}}} $为节点1的截面转动量。

由此可得局部坐标系下节点1的截面等效转动惯量矩阵为

${\hat {{I}}_1} = {\rm{diag}}\left[ {\begin{aligned} {{{\hat I}_{1x}}},\;{{{\hat I}_{1y}}},\;{{{\hat I}_{1z}}} \end{aligned}} \right] = \displaystyle\frac{1}{{12}}{\rm{diag}}\left[ {\begin{aligned} {\overline m_1^{\rm{e}}{{\left( {h_{\rm{a}}^{\rm{e}}} \right)}^2}},\;{\overline m_1^{\rm{e}}{{\left( {h_{\rm{a}}^{\rm{e}}} \right)}^2}},\;0 \end{aligned}} \right],$

式中： $ h_{\rm{a}}^{\rm e}$为薄壳元在t_a时刻的厚度. 再将其转换到全域坐标系下，可得 ${{I}}_1^{\rm{e}} = {\hat {{\varOmega }}^{\rm{T}}}{\hat {{I}}_1}\hat {{\varOmega }}$.

有限质点法采用显式直接积分法求解运动方程，以避免复杂迭代和收敛问题. 考虑数值稳定与物理精度两方面约束以及冲击振动中材料变形的影响，本文采用变时间步长法^[4]，其临界步长满足

(17) $ \Delta {t_{\rm{cr}}} = \frac{2}{{{\omega _{\max }}}}\left( {\sqrt {{\xi ^2} + 1} - \xi } \right) \leqslant \mathop {\min }\limits_{\rm shell}\; \frac{{l_{\rm{c}}^{\rm s}}}{{{c^{\rm s}}}}\left( {\sqrt {{\xi ^2} + 1} - \xi } \right). $

式中： $ \omega _{\max }$为系统的最大固有振动频率； $\xi $为系统最高振型对应的阻尼比系数；c^s为壳元内当前波速；

${c^{\rm s}} = \sqrt {{E / {\left[ {\rho {{\left( {1 - \nu } \right)}^2}} \right]}}} ,$

其中，E、ν和ρ分别为材料的弹性模量、泊松比和密度； $l_{\rm c}^{\rm s}$ 为壳元的最小特征长度，

$l_{\rm c}^{\rm s} = {{2{A_{\rm{s}}}} / {\max \left( {{L_1},{L_2},{L_3}} \right)}},$

其中，A_s为壳元面积，L_i（i=1，2，3）为边长.

利用中央差分格式，可将质点的速度和加速度写成位移和时间的表达式，代入式（1）后可得各时间步质点j切平面内沿独立转动自由度方向的转角和转速分别为（为表示方便，省略下标j）

(18) $ \tilde {\dot {{\beta }}}_{n + 1}^* = {c_2}\tilde {\dot {{\beta }}}_n^* + {c_1}\Delta {t_n}{\left( {\tilde {{I}}_n^*} \right)^{ - 1}}\left( {\tilde {{M}}_n^{*{\rm ext}} + \tilde {{M}}_n^{{\rm{*}}{\rm int}}} \right); $

(19) $ \begin{split} &\tilde {{\beta }}_{n + 1}^* = \frac{1}{{1 + {c_2}}} \times \\ &\left[ {2{c_1}\tilde {{\beta }}_n^* + 2{c_2}\Delta {t_n}{\tilde {\dot {{\beta }}}}_n^* + {c_1}\Delta t_n^2{{\left( {\tilde {{I}}_n^*} \right)}^{ - 1}}\left( {\tilde {{M}}_n^{*{\rm ext}} + \tilde {{M}}_n^{{\rm{*int}}}} \right)} \right]. \end{split}$

式中： ${c_1} = {\left( {1 + {{{\alpha _{\rm{m}}} \cdot \Delta {t_n}} / 2}} \right)^{{\rm{ - 1}}}}$, ${c_2} = {c_1}\left( {1 - {{{\alpha _{\rm{m}}} \cdot \Delta {t_n}} / 2}} \right),\;$ $ \Delta {t_n} =$ $ {t_n}_{ + 1} - {t_n}\,$.

实际计算中为保证一定余量，将初始步长Δt₀在临界步长基础上折减，即 $\Delta {t_0} = \kappa \Delta t_{\rm{cr}}^0$，这里 $\kappa $ 为考虑非线性不稳定因素的折减系数，通常可取 $\kappa $=0.9. 另外，后续第n步的临界步长若满足 $ \Delta t_{\rm{cr}}^n \geqslant$ $ \zeta \Delta {t_{n - 1}}$（ $ \zeta $为可靠度系数，可取1.05~1.10），则步长不变，即 $\Delta {t_n} = \Delta t_{n - 1}$，否则需重新设定时间步长，取 $\Delta t_n = \alpha \Delta t_{\rm{cr}}^n$.

根据瑞利比例阻尼的定义知阻尼比ξ_i与圆频率ω_i的关系^[4]：

(20) $ {{{a_{\rm{m}}}} / {\left( {2{\omega _i}} \right)}} + {{\left( {{a_{\rm{s}}}{\omega _i}} \right)} / 2} = {\xi _i}. $

对于固定的阻尼比值，a_m和a_s是随频率变化的. 作为一种近似，通常可由所关心频率范围上、下限的阻尼比（一般由实验方法或经验确定），根据式（21）得到联立方程定出a_m和a_s，作为该频率段的系数使用. 需要指出，ω_i的选取与主导振频有关，对于多数低频响应占主导的动力学问题，计算时可只取最低几阶频率范围，而对于冲击等问题则必须考虑激发的高频响应，具体取值可通过结构自由振动数值模拟试验确定. 另外，刚度阻尼主要对高频起作用，而且若a_s过大，临界步长会大大缩短，因此对于低频动力问题也可忽略刚度阻尼而通过适当提高a_m方式来考虑总阻尼力或阻尼力矩中a_s的影响.

基于经典J₂塑性理论^[22]，讨论平面应力回退映射算法及其在本文方法中的实现步骤.

为满足壳面内平面应力状态的限制条件，特别构建应力和应变二维子空间，相应的应力、应变矢量可表示为

$ \begin{split} &\hat {{\sigma}} = {\left[ {\begin{aligned} {{{\hat \sigma }_{11}}},\;{{{\hat \sigma }_{22}}},\;{{{\hat \tau }_{12}}} \end{aligned}} \right]^{\rm{T}}}, \; \hat {\varepsilon} = {\left[ {\begin{aligned} {{{\hat {{\varepsilon }}}_{11}}},\;{{{\hat \varepsilon }_{22}}},\;{{{\hat \gamma }_{12}}} \end{aligned}} \right]^{\rm{T}}},\\ & \hat {\alpha } = {\left[ {\begin{aligned} {{{\hat {{\alpha }}}_{11}}},\;{{{\hat \alpha }_{22}}},\;{{{\hat \alpha }_{12}}} \end{aligned}} \right]^{\rm{T}}}. \end{split}$

根据J₂理论流动准则和硬化律，利用向后Euler差分法可得t时刻塑性回退映射各状态量如下：

(21) $ {}^t{{\hat{ \varepsilon }}^{\rm{p}}} = {\hat{ \varepsilon }}_a^{\rm{p}} + \Delta {{\hat{ \varepsilon }}^{\rm{p}}} = {\hat{ \varepsilon }}_a^{\rm{p}} + \Delta \lambda {{P}}{}^t{\hat{ \xi }}, $

(22) $ {}^t{\bar \varepsilon ^{\rm{p}}} = \bar \varepsilon _a^{\rm{p}} + \Delta {\bar \varepsilon ^{\rm{p}}} = \bar \varepsilon _a^{\rm{p}} + \frac{2}{3}\Delta \lambda {}^t{\bar \sigma _{\rm{e}}}, $

(23) $ {}^t\hat {{\alpha }} = {\hat {{\alpha }}_a} + \Delta \hat {{\alpha }} = {\hat {{\alpha }} _a} + \frac{2}{3}\Delta \lambda \left( {1 - \beta } \right){H'_{\rm{p}}}{}^t\hat { \xi} . $

式中： $\Delta \lambda $ 为塑性流动因子； ${\hat {{\varepsilon }}^{\rm{p}}}$ 为塑性应变； ${\bar \varepsilon ^{\rm{p}}}$ 为等效塑性应变； $ {\bar \sigma _{\rm{e}}} = {\left( {1.5{{\hat {{\xi }}}^{\rm{T}}}{{P}}\hat {{\xi }}} \right)^{{\rm{1/2}}}} $为等效应力； $\hat {{\xi }} = \hat {{\sigma }} - \hat {{\alpha }}$ 为相对应力； $\hat {{\alpha }}$ 为背应力，代表子空间屈服面中心； $\beta $ 为混合硬化参数， $\beta = 0$ 和 $\beta = 1$时对应的硬化率分别简化为各向同性强化与随动强化； ${H'_{\rm{p}}}$ 为材料塑性模量，对应σ_T- $\varepsilon _{\ln }^{\rm{p}}$试验曲线斜率（σ_T为真应力， $\varepsilon _{\ln }^{\rm{p}}$为对数塑性应变），即 ${H'_{\rm{p}}} = {{\rm d}{{\bar \sigma }_{\rm{e}}}} /$ $ {{\rm d}{{\bar \varepsilon }^{\rm{p}}}} = {{{\rm d}\sigma _{\rm{T}}} / {{\rm d}\varepsilon _{\ln }^{\rm{p}}}}$；P为偏量投影矩阵.

将式（22）~（24）代入本构方程

${}^t\hat {{\sigma }} = {\hat {{\sigma }}_a} + {\hat {{D}}_{\rm{e}}}\left( {\Delta \hat {{\varepsilon }} - \Delta {{\hat {{\varepsilon }}}^{\rm{p}}}} \right) = {}^t{\hat {{\sigma }}^*} - {\hat {{D}}_{\rm{e}}}\Delta {\hat {{\varepsilon }}^{\rm{p}}},$

可得到如下参量更新公式：

(24) $ {}^t{\hat {{\xi }}^*} = {}^t{\hat {{\sigma }}^*} - {\hat {{\alpha }}_a}, $

(25) $ {}^t\hat {{\xi }} = {\left[ {1 + \frac{{\rm{2}}}{{\rm{3}}}\Delta \lambda \left( {1 - \beta } \right){{H'}_{\rm{p}}}} \right]^{{\rm{ - 1}}}}\hat {{D}}_{\rm{e}}^{\rm{cor}}\left( {\Delta \lambda } \right)\hat {{D}}_{\rm{e}}^{ - 1}{}^t{\hat {{\xi }}^*}. $

式中： ${}^t{\hat {{\sigma }}^*}$、 ${}^t\hat {{\xi }}$ 分别为弹性试探应力与相对应力； ${\hat {{D}}_{\rm{e}}}$ 为平面应力问题的弹性本构矩阵； $\hat {{D}}_{\rm{e}}^{\rm{cor}}$ 相当于对 ${\hat {{D}}_{\rm{e}}}$ 的修正，

$\hat {{D}}_{\rm{e}}^{\rm{cor}} = {\left[ {\hat {{D}}_{\rm{e}}^{ - 1} + {{\Delta \lambda {{P}}} \Bigg/ {\left( {1 + \frac{2}{3}\left( {1 - \beta } \right){{H'}_{\rm{p}}}\Delta \lambda } \right)}}} \right]^{ - 1}}.$

由式（22）~（26）可知，各状态量的塑性修正值均取决于 $\Delta \lambda $，可以由t时刻的一致性条件（J₂材料服从Mises屈服准则）来确定，即

(26) $ {}^t\varphi \left( {\Delta \lambda } \right) = {}^t\bar \sigma _{\rm{e}}^2 - \sigma _{\rm{s}}^2\left( {{}^t\bar \varepsilon ^{\rm{p}}} \right) \equiv 0. $

式中： ${}^t\bar \sigma _{\rm{e}}$ 为t时刻等效应力； ${\sigma _{\rm{s}}}$ 为屈服应力，可由材料应力强化曲线按当前等效塑性应变值确定.

由于 ${}^t\hat \xi $ 是 $\Delta \lambda $ 的一个非线性函数，因此式（27）构成了一个关于 $\Delta \lambda $ 的非线性代数方程. 对于弹性各向同性材料，利用矩阵 $\hat {{D}}_{\rm{e}}$ 和P的结构特征条件：

$ {{P}} = {{U}}{{{\varLambda }}_{\rm{P}}}{{{U}}^{\rm{T}}}, \; {\hat {{D}}_{\rm{e}}} = {{U}}{{{\varLambda }}_{\rm{D}}}{{{U}}^{\rm{T}}}. $

式中：U为正交矩阵， ${{{\varLambda }}_{\rm{P}}}$、 ${{{\varLambda }}_{\rm{D}}}$ 分别为 ${{P}}$ 与 ${\hat {{D}}_{\rm{e}}}$ 的对角矩阵. 可以将一致性条件改写成如下简单形式：

(27) $ \left. {\begin{array}{*{20}{l}} {{}^t\bar \sigma _{\rm{e}}^2\left( {\Delta \lambda } \right) = \displaystyle\frac{1}{2}\left( {\displaystyle\frac{{{\varGamma _1}}}{{{\varXi _1}}} + \displaystyle\frac{{{\varGamma _2}}}{{{\varXi _2}}}} \right)},\\ {\sigma _{\rm{s}}^2\left( {\Delta \lambda } \right) = \sigma _{\rm{s}}^2\left[ {{}^t\bar \varepsilon _a^{\rm{p}}\left( {\Delta \lambda } \right),\;\;{}^t\!\!\;\;\dot {\bar \varepsilon }_a^{\rm{p}}\left( {\Delta \lambda } \right)} \right]},\\ {{}^t\varphi \left( {\Delta \lambda } \right) = {}^t\bar \sigma _{\rm{e}}^2\left( {\Delta \lambda } \right) - \sigma _{\rm{s}}^2\left( {\Delta \lambda } \right) \equiv 0}. \end{array}} \right\} $

式中：

(28) $ \left. \begin{aligned} & {\varGamma _1} = {\left( {{}^t\hat \eta _{{\rm{u}}\_11}^*} \right)^2},\\ & {\varXi _1} = {\left[ {{\rm{1 + }}\left( {\frac{{{E}}}{{{\rm{3}}\left( {{\rm{1}} - {{v}}} \right)}} + \frac{2}{3}\left( {1 - \beta } \right){{H'}_{\rm{p}}}} \right)\Delta \lambda } \right]^2},\\ & {\varGamma _2} = 3{\left( {{}^t\hat \eta _{{\rm{u}}\_22}^*} \right)^2} + 6{\left( {{}^t\hat \eta _{{\rm{u}}\_12}^*} \right)^2},\\ & {\varXi _2} = {\left[ {{\rm{1 + }}\left( {2G + \frac{2}{3}\left( {1 - \beta } \right){{H'}_{\rm{p}}}} \right)\Delta \lambda } \right]^2},\\ & {\left[ {\begin{array}{*{20}{c}} {\hat \eta _{{\rm{u}}\_11}^*},\;{\hat \eta _{{\rm{u}}\_22}^*},\;{\hat \eta _{{\rm{u}}\_12}^*} \end{array}} \right]^{\rm{T}}} = \hat {{\eta }}_u^* = {{{U}}^{\rm{T}}}\hat {{\eta }}^* . \end{aligned} \right\} $

式（28）可利用Newton法局部迭代求解得到 $\Delta \lambda $，再代入式（22）~（26），就可以得到新的状态量.

上述基于平面应力J₂流动理论的塑性应力更新算法在本文方法中的具体实现步骤如下.

1）由各积分点t_a−t时段内的平面应变增量 $\Delta \hat {{\varepsilon }}$（假设为弹性应变）计算弹性试探应力 ${}^t{\hat {{\sigma }}^*} = {\hat {{\sigma }}_{\rm{a}}} + $ $ {\hat {{D}}_{\rm{e}}}\Delta \hat {{\varepsilon }}$.

2）将t时刻的应力 ${}^t{\hat {{\sigma }}^*}$ 及t_a时刻的塑性内变量 $\bar \varepsilon _{\rm{a}}^{\rm{p}}$ 和 ${\hat {{\alpha }}_{\rm{a}}}$ 代入屈服准则（式（27））计算弹性试探值：

$\varphi \left( {{}^t{{\hat {{\sigma }}}^*},{{\hat {{\alpha }}}_{\rm{a}}},\bar \varepsilon _{\rm{a}}^{\rm{p}}} \right) = \bar \sigma _{\rm{e}}^2\left( {{}^t{{\hat {{\sigma }}}^*} - {{\hat {{\alpha }}}_{\rm{a}}},} \right) - \sigma _{\rm{s}}^2\left( {\bar \varepsilon _{\rm{a}}^{\rm{p}}} \right)$.

若 $\varphi \leqslant 0$，说明为弹性步，有 $\Delta \hat {{\sigma }} = \Delta {\hat {{\sigma }}^*}$， $\Delta \bar \varepsilon ^{\rm{p}} = \Delta \hat {{\varepsilon }}^{\rm{p}} = $ $\Delta \hat {{\alpha }} = 0 $；若 $\varphi > 0$，说明为塑性加载步，应进行下面的塑性修正.

3）利用回退映射算法的应力修正过程如下.

a）设迭代初始值：i=0， ${}^t\bar \varepsilon ^{{\rm p}\left( 0 \right)} = \bar \varepsilon _{\rm{a}}^{\rm{p}}$， $\Delta {\lambda ^{\left( 0 \right)}} = 0$， ${}^t{\hat {{\sigma }}^{\left( 0 \right)}} = {\hat {{\sigma }}_{\rm{a}}}$， ${}^t{\hat {{\varepsilon }}^{{\rm p}\left( 0 \right)}} = \hat {{\varepsilon }}_{\rm{a}}^{\rm{p}}$， ${}^t{\hat {{\alpha }}^{\left( 0 \right)}} = {\hat {{\alpha }}_{\rm{a}}}$；

b）在第i次迭代时由式（27）检查屈服条件，如果函数值 ${}^t{\varphi ^{\left( i \right)}} < {\rm TOL}$（其中TOL是容差限值），说明迭代已收敛，可停止迭代；

c）利用式(28)由

$ {}^t{\varphi ^{\left( {i + 1} \right)}} = {}^t{\varphi ^{\left( i \right)}} + {\left( {{{{\rm d}{}^t\varphi } / {{\rm d}\Delta \lambda }}} \right)^{\left( i \right)}}$ $\delta {\lambda ^{\left( i \right)}} = 0$

计算塑性流动因子 $\Delta {\lambda ^{\left( i \right)}}$；

d）将 $\Delta {\lambda ^{\left( i \right)}}$ 代入式(22)~(26)，计算应力和内变量，依次得到 ${}^t{\hat {{\xi }}^{\left( {i + 1} \right)}}$、 ${}^t{\hat {{\alpha }}^{\left( {i + 1} \right)}}$、 ${}^t{\bar {{\varepsilon }}^{{\rm p}\left( {i + 1} \right)}}$、 ${}^t{\hat {{\varepsilon }}^{{\rm p}\left( {i + 1} \right)}}$，同时令i=i+1，转至2)继续迭代，直至收敛；

4）更新当前时刻的状态量：

${}^t{\bar \varepsilon ^{\rm{p}}} = \bar \varepsilon _{\rm{a}}^{\rm{p}} + \Delta {\bar \varepsilon ^{\rm{p}}},\;{}^t{\hat {{\varepsilon }} ^{\rm{p}}} = \hat {{\varepsilon }}_{\rm{a}}^{\rm{p}} + \Delta \hat {{\varepsilon }}^{\rm{p}},$

${}^t\hat {{\alpha }} = {\hat {{\alpha }}_{\rm{a}}} + \Delta \hat {{\alpha }} ,$

${}^t\hat {{\sigma }} = {\hat {{\sigma }}_{\rm{a}}} + {\hat {{D}}_{\rm{e}}}\left( {\Delta \hat \varepsilon - \Delta {{\hat \varepsilon }^{\rm{p}}}} \right),$

${}^t{\hat \varepsilon _z} = - {\nu \Big/ { {E\left( {{}^t{{\hat \sigma }_x} + {}^t{{\hat \sigma }_y}} \right)} }} - \left( {{}^t\varepsilon _x^{\rm{p}} + {}^t\varepsilon _y^{\rm{p}}} \right).$

本文方法的质点控制方程内力项( ${}^t{{M}}^{\rm{int}}$)基于当前步末各壳元内的应力状态( ${}^t\hat {{\sigma }},{}^t\hat {{\alpha }} ,{}^t\bar \varepsilon ^{\rm{p}}$)，而这些状态量可通过当前时刻已知的位移增量直接求出，因此对方法本身不需再作额外修正.

基于前文理论计算公式推导及相关参数的取值策略，同时引入弹塑性应力修正的处理，采用C++语言在Visual Studio平台实现计算分析程序的编写和运行.

整个有限质点法分析程序的实现流程分成前处理、计算求解和后处理. 前处理包括结构模型信息读入（如位置、材料、几何等）和初始分析参数设置（如步长、约束等）. 后处理主要是输出计算结果. 求解部分由2层循环构成：外层按时间步循环，内层包含点值循环和单元循环. 点值循环调用该步的内、外力（矩）和阻尼力（矩）等信息采用显式积分求解运动方程获得并更新质点位移（转角）和速度（转速）；单元循环通过虚拟运动计算壳元纯变形并获得内力（矩），然后集成到相对应的质点上，并把结果返回给点值循环供下一步质点运动方程积分计算使用. 材料非线性等特殊情况处理可在内层循环中作为单独的计算模块接入，在每次循环中对单元或质点状态作判断，然后调用相应模块. 下面通过若干典型算例分析验证理论推导和所编制程序的有效性和正确性.

受持续均布压力作用的扁球壳结构是传统薄壳非线性分析中的经典算例^[23]，常用来验证算法的正确性. 扁球壳周边固定约束，从t=0时刻起，上表面施加恒定均布荷载p=6.137 MPa（为方便与文献对比，这里采用英制单位），作用时间为1 ms，不考虑阻尼效应. 几何尺寸、约束条件、材料性质及受力情况如图5所示. 图中，R为扁球壳半径，α为弧度，h为厚度，E为材料弹性模量，v为泊松比，ρ为密度，σ_y为屈服应力，E_t为切线模量。利用对称性条件，只取1/4部分进行计算. 壳体材料模型考虑线弹性和双线性等向强化弹塑性2种情况. 采用质点数分别为127、469和631的3种离散模型（记作模型Ⅰ ~ Ⅲ，见图6）进行计算，检验算法的收敛性.

扁球壳在突加恒定面荷载作用下产生振动. 图7（a）给出由3种离散模型计算得到的弹性球壳顶点竖向位移w的时程曲线，与Oñate等^[24-25]采用非线性有限元法的求解结果进行对比. 可以看出，采用不同质点密度的离散模型得到的顶点位移时程曲线的变化趋势与文献结果均比较接近，并且随着模型细化，模拟结果趋于吻合，体现了本文方法良好的收敛性. 另外，图7（b）还给出了由模型Ⅲ计算得到的顶点弹塑性竖向位移时程曲线，与文献[23]和[24]结果基本一致. 对比图7（a）和（b），可以看到t=0.16 ms时扁球壳内有部分区域开始进入塑性阶段，之后受刚度软化和塑性变形的影响，与弹性模型相比结构振幅有所减小.

表1列出了不同方法计算得到的若干中间时刻的顶点竖向位移值. 其中，本文方法与Bathe等^[23-25]给出结果的最大偏差仅分别为17.1%、4.6%和5.2%，进一步证明了本文方法具有较高的计算精度. 通过本算例初步验证了有限质点法处理薄壳结构动力非线性问题的有效性.

如图8所示，本算例分析周边固定的120°柱面曲壳受到局部瞬时冲击作用后的动力大变形过程. Witmer等^[26]曾对此做过试验研究，其结果被许多学者用来验证各类数值计算理论的可靠性^{[24, 27]}. 柱面壳的几何尺寸、材料性质及初始冲击作用如图8所示. 在图中25.9 cm×7.8 cm条形阴影区域内对质点施加143.5 m/s初始法向速度，用以模拟爆炸冲击作用，分析柱面壳受冲击后1.0 ms内的振动响应. 利用对称性条件，只取一半模型计算，共布置13×33个离散质点. 采用理想弹塑性模型进行分析，计算中不考虑阻尼效应.

图9给出了该柱面壳受冲击后的变形过程，其中沿x向和y向对称面（即y=15.95 cm和x=0截面）的变形后截面形状（t=1 ms）分别如图10（a）和（b）所示. 图11详细记录了纵轴线上y=23.93 cm和y=15.95 cm两点处的竖向位移w随时间的变化情况. 图10和11中的模拟结果与试验值^[26]的对比表明，即使在高频模态被激发且振动响应非线性程度较高的情况下，本文计算结果与试验值仍能较好吻合，说明该方法可以用于求解薄壳弹塑性冲击动力问题.

薄壁金属方管具备良好的耗能特性，在结构抗震和汽车防撞吸能盒设计中都有广泛应用. 以如图12所示几何尺寸和材料的薄壁方管为例，将其一端固定，另一端用一大质量刚体块沿轴向撞击，刚块质量m₀=800 kg，撞击速度v₀=9.5 m/s. 由于惯性，撞击会导致管壁产生连续动力屈曲，并伴随波浪式卷曲折叠的传播与扩展. 该过程包括了材料大变形、动力屈曲、碰撞接触等多种强非线性行为，也是检验动力学算法可靠性的典型算例^[28-29]. 为减少计算量，只取方管的1/4对称部分分析，各侧面布置129×31个离散质点，模拟时长共18 ms. 对于管壁皱曲过程中的自接触现象，采用显式拉格朗日乘子法^[4]求出接触力并反向施加到相应范围的质点上. 管体材料的应力-应变关系符合Johnson-Cook模型，计入撞击后的应力硬化和应变率效应影响，满足混合强化法则并取参数β=0.5.

采用以上模型计算得到刚块撞击反力R随时间变化情况及管端部一点A处的纵向加速度a_z时程曲线分别如图13（a）和（b）所示，图中的一系列峰值点对应着管体受撞击后发生的连续多次屈曲现象. 图14给出了撞击后动力屈曲演化过程中t=3.5、9.6、12.8、17.4 ms方管的变形情况. 可以看到，管端被撞击后承受的冲击力迅速增大，使其靠近刚块一端进入塑性状态，并导致约t=3.5 ms第1层波状皱曲出现并形成塑性铰；之后，由于管壁皱曲，能量被吸收，同时方管刚度减小，刚块撞击反力也有所降低；但随着皱曲扩展，撞击应力波回弹，同时波浪式翻卷的管壁之间发生挤压，接触力逐渐增大，皱曲部位继续产生塑性应变强化，因此刚块反力又呈现出一定的上升趋势；最终，撞击能量通过管壁连续屈曲形式被完全吸收，在各个侧面上形成4层明显的波状皱曲. 图15还给出了第4层皱曲发生后（约t=13.4 ms）管壁内外两侧的等效塑性应变云图，最大值分别为1.348和1.502，表明此时已产生较大塑性变形.

算例结果表明，尽管管壁材料已产生极大变形，但本文方法仍然实现了对撞击后全过程的有效模拟. 整个计算中没有额外设置特殊技巧或进行人工参数的修正，也没有出现因强非线性行为而导致的数值不稳定现象，反映出本文方法的良好性能.

本文将有限质点法应用于求解薄壳结构的动力非线性问题. 构造了具有大转动大变形分析能力的薄壳离散计算模型，建立了质点和壳元的基本控制方程与计算公式，并对质点转动惯量、临界时间步长等关键参数取值问题提出了合理可行的处理策略. 针对薄壳质点位移及壳元变形与内力均按照面内拉压和面外弯扭两部分进行拆分和组合，并着重推导了与面外弯扭刚度相对应的质点内力（矩）计算公式，给出了质点切平面外转角变步长显式时间积分式. 与直接构建高阶的等参曲面壳元相比，这种组合内力分块计算方式简单，运算量小，也避免了因剪切自锁和零能模态等问题而带来的求解困难. 此外，还在质点内力求解模块引入了基于J₂塑性理论的平面应力通用本构积分算法，建立了算法的具体实现步骤，可用于分析理想弹塑性、混合应力强化、率相关塑性等多种材料模型. 通过多个典型数值算例对本文方法进行了测试，并将模拟结果与相关文献的数值和试验结果进行了对比，表明该方法从理论到程序均能有效求解考虑薄壳几何与材料大变形的爆破冲击、撞击等动力非线性问题，且具有较高精度.

有限质点法不对结构的线性和非线性运动变形加以区别，也不存在非线性迭代的概念，而采取统一的分析模式. 另外，质点在内力和外力作用下，其运动状态是永恒的，质点控制方程的解本身就是结构的动力响应，其动力分析过程更加真实结构的行为过程. 这些特点使得非线性和动力分析在有限质点法中成为一种自然的过程. 因而该方法在薄壳结构动荷载下的强非线性行为计算中有很大的优势和灵活性，给薄壳研究中的这类问题求解提供了新的思路和强有力的解决方法.