周书涛¹, 马斌捷¹, 侯传涛¹, 童军¹, 巨亚堂¹, 刘应华²
1. 北京强度环境研究所, 北京 100076;
2. 清华大学 航天航空学院, 工程力学系, 北京 100084
收稿日期：2020-05-25
基金项目：中国博士后科学基金（2013M540934）
作者简介：周书涛(1983-), 男, 高级工程师。E-mail: zhoushutao-1983@163.com

摘要：采用C¹自然单元法研究了不同工况下圆形、菱形、等边多边形薄板的极限承载力。根据薄板极限上限分析的迭代求解格式，构造出了满足平衡方程和边界条件的广义应力场，并由极限下限定理和得到的广义应力场，建立了求解薄板结构极限下限载荷乘子的迭代格式。提出的数值方法克服了极限下限定理中约束条件的强非线性，降低了下限分析的计算规模，具有易于程序实现的优点。该数值方法与极限上限分析方法相结合可以有效估算出薄板结构极限载荷的范围。数值算例表明，提出的求解薄板结构上、下限载荷的方法是有效的，具有较高的计算精度和较快的收敛性。
关键词：极限分析上限定理下限定理C¹自然单元法直接迭代算法
C¹ natural element method for plastic limit analyses of thin plates
ZHOU Shutao¹, MA Binjie¹, HOU Chuantao¹, TONG Jun¹, JU Yatang¹, LIU Yinghua²
1. Beijing Institute of Structure & Environment Engineering, Beijing 100076, China;
2. Department of Engineering Mechanics, School of Aerospace, Tsinghua University, Beijing 100084, China

Abstract: The C¹ natural element method (C¹ NEM) was used to study the limiting loads of circular, rhombic, and equilateral polygon thin plates subjected to various loading conditions. An iterative solution for the upper load limits of the thin plates made the generalized stress fields satisfy the equilibrium equations and the boundary conditions. Iterative solutions were also used to calculate the lower limits of the load multipliers of thin plates using the lower bound theorem to obtain the generalized stress fields. This numerical method overcomes the difficulties introduced by the strong nonlinearity of the constraint condition in the lower bound theorem and reduces the calculations for the lower bound analysis in an easily implemented algorithm. This numerical approach can also be incorporated into upper bound analyses to estimate the limiting loads of thin plates. Numerical examples show that this numerical method can accurately and quickly predict the upper and lower load limits of thin plates.
Key words: limit analysisupper bound theoremlower bound theoremC¹ natural element method (C¹ NEM)direct iteration algorithm
板结构是航空、航天、船舶等工程领域广泛使用的一种结构形式。准确、可靠获得板结构的极限承载能力，对其工程设计和安全评估有着非常重要的意义。相对于传统的线弹性设计而言，将塑性力学的相关理论应用到板结构的设计和评估中，除了能在充分利用材料塑性性能的基础上得到能真实反映结构安全裕度的重要参数之外，还在降低成本和节约资源方面具有重要的理论研究意义和实际应用价值。塑性极限分析是塑性力学的一个重要分支，能够直接获得结构在比例载荷作用下的极限载荷，是研究板结构极限承载能力的一种直接而有效的方式。
塑性极限分析理论建立于Drucker和Hill等^[1]提出的极限上、下限定理。该理论认为，只有当结构整体进入屈服状态，且在外载荷作用下变成几何可变结构时，结构才会因承载能力的丧失而破坏。极限上、下限定理适用于不同的领域，上限定理适用于塑性压力加工和岩土工程，下限定理适用于结构工程设计和安全评定^[1]。国内外的众多****在应用塑性极限分析理论研究板结构的极限承载能力方面开展了大量系统而深入的工作。一部分****使用Mises、Tresca、统一强度等准则和理论，通过理论推导获得板结构极限载荷的解析解。国内外相关系统性的研究成果可见文[2-7]。另一部分****借助于目前广泛使用的有限元、无网格等数值方法，研究和发展了一系列获得板结构极限载荷的高效优化算法。例如：Hodge和Belytschko^[8]采用有限单元法获得了方板和矩形板的极限载荷; Capsoni和Corradi^[9]采用有限单元法对理想刚塑性Kirchhoff和Mindlin板结构进行了极限上限分析; Turco和Caracciolo^[10]采用混合三角形(HS)有限单元和Newton-Raphson算法对圆形和方形薄板进行了极限分析; Le等^[11]基于二阶凸规划(SOCP)理论和无单元法(EFG)，得到了Kirchhoff薄板的极限上限载荷; 同样基于SOCP求解算法，Bleyer和Buhan^[12]分别采用T6、T6b、H3有限单元，获得了基于Johansen准则的简支和夹支方板，基于Mises准则的夹支方板、简支和夹支圆板、简支L形板的极限上限载荷; Capsoni和Vicente da Silva^[13]通过使用混合b-bar应变法和Mindlin板单元，数值研究了圆板和斜板的极限上限解; Baghani和Fereidoonnezhad^[14]采用变分迭代法和Mises准则对由梯度功能材料组成的简支圆板进行了极限分析，获得了极限载荷的解析解和数值解; Makrodimopoulos^[15]使用6节点三角形单元和Nielsen、Mises准则，对夹支方薄板、不规则薄板进行了极限上限分析; Nguyen-Thoi等^[16]采用CS-DSG3、DSG3、MIN3等3种有限板单元，数值研究了基于Mises准则和Mindlin假设的夹支方板、简支矩形板、夹支菱形板、夹支圆形板、夹支等边三角形板的极限上限载荷。Zhou等^[17-18]基于极限上限分析理论发展了一套直接迭代算法，分别采用非协调矩形薄板单元ACM和C¹自然单元法(C¹ NEM)，对多种不同形状的薄板进行了极限上限分析。由于上述数值求解算法能够较为准确地获得各种复杂几何、载荷和边界条件下板结构的极限承载力，因而更加受广大****的关注。
工程设计与安全评定的技术人员往往更加关心下限解。相对于极限下限分析的不等式约束，极限上限分析的等式约束条件更容易处理。如何在极限上限分析的基础上开展极限下限分析，也是相关****的研究重点之一。陈浩峰^[19]和Zhou等^[20]分别利用三维、二维结构极限上限迭代格式中得到的中间变量，相应地构造出了结构下限分析，并验证了采用该应力场进行下限分析的有效性和准确性。鉴于此，本文采用新颖的C¹自然单元法(C¹ NEM)，充分发挥其计算精度高、便于施加位移边界条件、使用灵活等优势，在Zhou等^[18]薄板极限上限迭代格式的基础上，构造出了满足平衡方程和边界条件、适合薄板结构下限分析的广义应力场，并由极限下限定理建立了求解薄板结构极限下限载荷的迭代求解格式。该格式有效克服了极限下限定理中的强非线性约束条件，降低了传统下限分析的计算规模，只需在极限上限分析的计算程序中增加若干求解代码，易于程序实现。据此开发出了相应的求解程序，对不同工况下圆形、菱形、等边多边形薄板结构进行了极限分析，得到了相应的极限上、下限载荷乘子。同时，采用移动最小二乘技术，得到了能够直观反映结构塑性失效模式的Mises广义等效应力云图。
1 薄板结构极限分析理论1.1 薄板结构极限上限分析理论极限分析的上限定理可表述为^[1]：满足几何约束条件，能构成几何可变机构，且外载荷在其上所做的功率不小于结构内部塑性耗散功率的结构位移速率场，称为机动容许速率场，与机动容许速率场对应的载荷是极限载荷的上限。这些上限解中的最小者就是极限载荷。根据该定理可将薄板结构极限上限分析的数学规划格式写为^[17-18]：

${v^ + } = \min \int_S D \left( {{{\dot \kappa }_{ij}}(\mathit{\boldsymbol{x}})} \right){\rm{d}}S,$

(1)

${\rm{s}}{\rm{. t}}.{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \int_{{\mathit{\Gamma }_t}} {{{\dot a}_i}} {t_i}{\rm{d}}\mathit{\Gamma } = 1,$

(2)

${{\dot a}_{i,i}} = 0.$

(3)

其中：v⁺为极限上限载荷乘子; $D\left(\dot{\kappa}_{i j}\right)$是对应于机动容许广义应变率$\dot{\kappa}_{i j}(\boldsymbol{x})$的塑性耗散功率，利用Mises屈服准则可以写为$D\left(\dot{\kappa}_{i j}\right)=\sqrt{2 / 3} M_{\mathrm{p}} \sqrt{\dot{\kappa}_{i j} \dot{\kappa}_{i j}}$, M_p=σ_sh²/4是屈服应力为σ_s、厚度为h板结构的塑性极限弯矩; S和Γ_t分别为薄板结构的中面和应力边界; $\dot{a}_{i}$和t_i分别是机动容许的位移速率场和表面力。式(2)和(3)分别表示归一化条件和塑性不可压条件。
式(1)—(3)所示的上限规划格式为一个带等式约束的非线性规划问题，其求解的主要困难在于非线性目标函数的线性化和塑性不可压条件的处理。
1.2 薄板结构极限下限分析理论极限分析的下限定理可以表述为^[1]：如果结构的某个应力场满足平衡条件，且处处满足材料的屈服条件，则称为静力容许应力场。与该应力场对应的载荷是极限载荷的下限解，极限载荷是这些下限解中的最大者。根据该定理可将薄板结构极限下限分析的数学规划格式写为

$v^{-}=\max \beta, $

(4)

$f\left[\beta M_{i j}(\boldsymbol{x})\right]-M_{\mathrm{p}}^{2} \leqslant 0, \quad \forall \boldsymbol{x} \in S .$

(5)

其中: β是极限下限载荷乘子，v^-为这些乘子中的最大者，M_ij(x)为满足平衡方程和边界条件的广义应力场。式(5)为由屈服函数f[·]所表示的屈服条件。
式(4)和(5)所示的下限规划格式为一个需要同时满足等式约束(平衡方程)和不等式约束(屈服条件)的非线性规划问题，其求解的主要困难在于约束条件的强非线性。虽然将屈服条件近似线性化后可以较为方便地将问题进行处理，但是往往会极大地增加求解规模^[19]。
塑性极限载荷位于上、下限解之间。在上述塑性极限上、下限分析中采用的基本假设包括^[1]：理想塑性材料、小变形、比例加载、在达到极限载荷前结构不失稳。
2 C¹自然单元法(C¹ NEM)C¹自然单元法是Sukumar和Moran^[21]提出的一种基于离散点的Voronoi结构和Delaunay三角化、适合求解涉及到高阶偏微分方程问题的新型无网格法。该方法采用伽辽金法来离散和生成系统方程，采用C¹自然临近插值来构造形函数。C¹自然临近插值是将Sibson自然临近坐标插入到Bernstein Bézier曲面的单纯形中，构造出满足C¹连续性要求、能同时对节点值和节点梯度值进行插值的数值方法。
假设空间中的某点x具有n个自然临近节点，这n个自然临近节点的Sibson插值形函数分别表示为Φ₁(x), Φ₂(x), …, Φ_n(x)。Sibson插值形函数的计算草图见图 1(图中n=5，节点1、2、3、6、7为基于空圆准则确定的、点x的自然临近节点)，其具体计算方法见Sukumar等^[22]介绍的Watson算法。将$\mathit{\tilde \Phi }$=(Φ₁(x), Φ₂(x), …, Φ_n(x))表示成点x的自然临近坐标，并将其代入到Bernstein- Bézier曲面单纯形的表达式中，可构造得到如下的C¹自然临近插值形式^{[18, 21]}：

图 1 Sibson插值形函数计算草图

图选项

$w^{h}(\widetilde{\varPhi})=\boldsymbol{B}^{\mathrm{T}}(\widetilde{\varPhi}) \boldsymbol{b}=\boldsymbol{B}^{\mathrm{T}}(\widetilde{\varPhi}) \boldsymbol{Ta} =\varPsi^{\mathrm{T}}(\widetilde{\varPhi}) \boldsymbol{a}.$

(6)

$\begin{aligned}\varPsi^{\mathrm{T}}(\widetilde{\varPhi})=& \boldsymbol{B}^{\mathrm{T}}(\widetilde{\varPhi}) \boldsymbol{T}=\left[\varPsi_{1}(\widetilde{\varPhi}) \varPsi_{2}(\widetilde{\varPhi}) \varPsi_{3}(\widetilde{\varPhi}) \cdots\right.\\&\left.\varPsi_{3 n-2}(\widetilde{\varPhi}) \varPsi_{3 n-1}(\widetilde{\varPhi}) \varPsi_{3 n}(\widetilde{\varPhi})\right], \end{aligned}$

(7)

$\boldsymbol{a}=\left[\begin{array}{lllllll}w_{1} & \theta_{1 x} & \theta_{1 y} & \cdots & w_{n} & \theta_{n x} & \theta_{n y}\end{array}\right]^{\mathrm{T}}.$

(8)

其中：B($\mathit{\tilde \Phi }$)和b分别为Bernstein- Bézier基函数和Bézier坐标的列向量; T为转换矩阵; Ψ^T($\mathit{\tilde \Phi }$)为C¹自然临近形函数的矩阵形式; a是节点自由度矩阵，为3n维列向量。对于薄板问题而言，w_n=w(x_n)表示节点n的挠度，θ_nx=w_{, x}(x_n)和θ_ny=w_{, y}(x_n)则分别表示该点的旋度值。

$\varPsi_{3 I-2}\left(\boldsymbol{x}_{J}\right)=\delta_{I J}, \varPsi_{3 I-1}\left(\boldsymbol{x}_{J}\right)=0, \varPsi_{3 I}\left(\boldsymbol{x}_{J}\right)=0;$

(9)

$\varPsi_{3 I-2, x}\left(\boldsymbol{x}_{J}\right)=0, \varPsi_{3 I-1, x}\left(\boldsymbol{x}_{J}\right)=\delta_{I J}, \varPsi_{3 I, x}\left(\boldsymbol{x}_{J}\right)=0;$

(10)

$\varPsi_{3 I-2, y}\left(\boldsymbol{x}_{J}\right)=0, {\mathtt{ψ}}_{3 I-1, y}\left(\boldsymbol{x}_{J}\right)=0, {\mathtt{ψ}}_{3 I, y}\left(\boldsymbol{x}_{J}\right)=\delta_{I J}.$

(11)

C¹自然临近形函数Ψ($\mathit{\tilde \Phi }$)具有非负性、单位分解性、二次相容性等许多重要性质，其中最重要的性质是其对节点函数和节点梯度值的Kronecker delta插值特性(式(9)—(11))。该特性使其能够像有限元一样精确地施加位移边界条件。图 2b~2d给出了在图 2a所示支撑域上，节点I的C¹自然临近插值形函数云图。

图 2 (网络版彩图)节点支撑域和C1自然临近插值形函数云图

图选项

3 基于C¹ NEM的薄板极限分析3.1 薄板极限上限分析的迭代格式根据Mises屈服准则和C¹ NEM建立的薄板极限上限分析(式(1)—(3))的离散格式为^[18]：

$v^{+}=\min \sqrt{2 / 3} M_{\mathrm{p}} \sum\limits_{i \in \mathrm{IG}} \rho_{i}|\boldsymbol{J}|_{i} \sqrt{\boldsymbol{\dot a}^{\mathrm{T}} \boldsymbol{K}_{i} \boldsymbol{\dot d}} ,$

(12)

$\text { s. t. } \ \boldsymbol{F}^{\mathrm{T}} \dot{\boldsymbol{a}}=1.$

(13)

其中：IG为薄板结构离散后的全部积分点数; ρ_i和|J|_i分别是积分点i的积分权和Jacobi行列式; $\dot{\boldsymbol{a}}$为薄板结构离散后的整体节点位移速率矢量，是由采用空圆准则确定的、积分点i的n个自然临近节点的位移向量$\mathit{\boldsymbol{\dot {\tilde a}}}=\left[\begin{array}{ccc}\dot{w}_{1} \; \dot{\theta}_{1 x}\; \dot{\theta}_{1 y} \; \cdots \; \dot{w}_{n} \; \dot{\theta}_{n x} \; \dot{\theta}_{n y}\end{array}\right]$扩展组装而来，$\dot{\omega}_{j}、\dot{\theta}_{j x}$和$\dot{\theta}$_jy(j=1, 2, …, n)分别表示节点j处的挠度率、绕x轴的转角率和绕y轴的转角率; K_i=B_i^TDB_i为与薄板结构全部节点对应的矩阵，B_i为应变位移关系矩阵，D为一常数矩阵，D矩阵的引入使得薄板结构的塑性不可压条件(式(3))在目标函数中得到了满足; F为薄板结构的等效节点载荷向量。
目标函数式(12)为带根号的非线性项，为了直接对其进行求解，Zhou等^[18]构造出了在第k次迭代时只含节点位移速度向量$\mathit{\boldsymbol{\dot a}}$_k和拉格朗日乘子λ_k 2个未知变量的线性化迭代格式：

$\sum\limits_{i \in \operatorname{IG}} \boldsymbol{H}_{i}^{k} \boldsymbol{a}_{k} =\lambda_{k} \boldsymbol{F} ,$

(14)

$\boldsymbol{F}^{\mathrm{T}} \dot{\boldsymbol{a}}_{k} =1$

(15)

其中，

$\boldsymbol{H}_{i}^{k}=\sqrt{2 / 3} M_{\mathrm{P}} \rho_{i}\ |\boldsymbol{J}|_{i} \boldsymbol{K}_{i} / \hat{X}_{i}^{k}.$

(16)

在第k(k≥1)次迭代时，$\hat{X}_{i}^{k}$为已知变量。根据第(k-1)次迭代的结果，将全部积分点集IG按照如下的准则划分成塑性区子集P_k和刚性区子集R_k：

$\mathrm{IG}=P_{k} \cup R_{k}, $

(17)

$P_{k}=\left\{i \in \mathrm{IG}, X_{i}^{k}=\sqrt{\dot{\boldsymbol{a}}_{k-1}^{\mathrm{T}} \boldsymbol{K}_{i} \dot{\boldsymbol{a}}_{k-1}}>0\right\}, $

(18)

$R_{k}=\left\{i \in \mathrm{IG}, {X}_{i}^{k}=\sqrt{\dot{\boldsymbol{a}}_{k-1}^{\mathrm{T}} \boldsymbol{K}_{i} \dot{\boldsymbol{a}}_{k-1}}=0\right\} .$

(19)

考虑到实际计算中存在的数值误差，当某个积分点的X_i^k值太小而趋近于一个远小于1的非负小数γ(通常设置为γ=(10^－4~10^－8)×$\sum\limits_{i=1}^{\mathrm{IG}} \sqrt{\dot{\boldsymbol{a}}_{0}^{\mathrm{T}}\boldsymbol{K}_{i} \dot{\boldsymbol { a }}_{0}}$/IG))时，就认为该点进入刚性状态，此时取$\hat{X}_{i}^{k}$=γ。反之，如果积分点i的X_i^k值大于γ，就认为该点处于塑性状态，此时取$\hat{X}_{i}^{k}=X_{i}^{k}=\sqrt{\dot{\boldsymbol a}_{k-1}^{\mathrm{T}}\boldsymbol{K}_{i} \dot{\boldsymbol{a}}_{k-1}}$。按照这一划分准则，目标函数在每次迭代中都是有意义的，刚性积分点随着迭代的进行也逐渐被识别出来。
式(14)和(15)所示的线性化迭代格式对应如下带一个等式约束的二次极小化问题：

$\nu_{k}^{+} =\min \sum\limits_{i \in \operatorname{IG}} \dot{\boldsymbol{a}}^{\mathrm{T}} \boldsymbol{H}_{i}^{k} \dot{\boldsymbol{a}}, $

(20)

$\text { s. } \mathrm{t.} \ \ \boldsymbol{F}^{\mathrm{T}} \dot{\boldsymbol{a}}=1.$

(21)

使用Lagrange乘子λ将式(21)表示的归一化条件代入到式(20)表示的目标函数中去，再根据极小化条件得到如下关于未知变量$\dot{\boldsymbol { a }}$和λ的线性方程组：

$\sum\limits_{i \in \operatorname{IG}} \boldsymbol{H}_{i}^{k} \dot{\boldsymbol{a}}=\lambda \boldsymbol{F}, $

(22)

$\boldsymbol{F}^{\mathrm{T}} \dot{\boldsymbol{a}}=1.$

(23)

定义：

$\dot{\boldsymbol{a}}=\lambda \boldsymbol{\delta}.$

(24)

分别将式(24)代入到式(22)和式(23)中去，组装矩阵$\boldsymbol{G}_{k}=\sum\limits_{i=1}^{\mathrm{IG}} \boldsymbol{H}_{i}^{k}$并引入位移边界条件对G_k和F进行修改，可以相应地计算出δ=G_k^－1F和λ=1/(F^Tδ)。将得到的δ和λ代回到式(24)中，可以获得相应的节点位移速度向量$\dot{\boldsymbol { a }}$。从而可得第k次迭代时薄板极限上限载荷乘子为

$\nu_{k}^{+}=\sqrt{2 / 3} M_{\mathrm{P}} \sum\limits_{i \in \mathrm{IG}} \rho_{i}\ |\boldsymbol{J}|_{i} \ \sqrt{\dot{\boldsymbol{a}}^{\mathrm{T}} \boldsymbol{K}_{i} \dot{\boldsymbol{a}}_{k}} .$

(25)

3.2 薄板极限下限分析的迭代格式根据薄板结构极限分析的下限定理(式(4)—(5))，采用C¹ NEM离散后的极限下限载荷乘子可以写为

$v^{-} =\min \limits_{i \in \operatorname{IG}} \beta_{i},$

(26)

$\text { s. t. }\ f\left(\beta_{i} M_{i}\right)=0, \quad i=1,2, \cdots, \text { IG. }$

(27)

对于薄板结构而言，根据von Mises屈服条件可得

$f(\boldsymbol{M})=\boldsymbol{M}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{M}-1=0.$

(28)

其中，

$\boldsymbol{M}=\left[\begin{array}{llll}M_{11} & M_{22} & M_{12}\end{array}\right]^{\mathrm{T}},$

(29)

$\boldsymbol{A}=\frac{1}{\boldsymbol{M}_{\mathrm{P}}^{2}}\left[\begin{array}{ccc}1 & -0.5 & 0 \\-0.5 & 1 & 0 \\0 & 0 & 3\end{array}\right] .$

(30)

由式(14)—(15)所表示的、上限分析的第k次迭代格式可以写为

$\left(\sqrt{2 / 3} M_{\mathrm{P}} \sum\limits_{i \in \mathrm{IG}} \rho_{i}|\boldsymbol{J}|_{i} \boldsymbol{B}_{i}^{\mathrm{T}} \boldsymbol{D} \boldsymbol{B}_{i} / \hat{X}_{i}^{k}\right) \dot{\boldsymbol{a}}_{k}=\lambda_{k} \boldsymbol{F} ,$

(31)

$\boldsymbol{F}^{\mathrm{T}} \dot{\boldsymbol{a}}_{k}=1.$

(32)

根据式(24)可知，在第k次迭代时有$\dot{\boldsymbol { a }}$_k=λδ_k。将其代入到式(31)中可得

$\sum\limits_{i \in \mathrm{IG}} \rho_{i}|\boldsymbol{J}|_{i} \boldsymbol{B}_{i}^{\mathrm{T}} \boldsymbol{M}_{i k}=\boldsymbol{F}.$

(33)

其中，

$\boldsymbol{M}_{i k}=\sqrt{2 / 3} M_{\mathrm{P}} \boldsymbol{D} \boldsymbol{B}_{i} \boldsymbol{\delta}_{k} / \hat{X}_{i}^{k}.$

(34)

由式(33)确定的第k次迭代的广义应力场M_ik显然是与所施加的比例载荷F平衡。同时由式(34)可知，变量δ_k和$\hat{X}_{i}^{k}$均是在满足位移边界条件下计算得到的，因此广义应力场M_ik显然也是满足位移边界条件的。
将每个积分点上的广义应力场β_{i, k}M_ik代入式(28)表示的屈服条件中进行检查，即：

$\left(\beta_{i, k} \boldsymbol{M}_{i k}\right)^{\mathrm{T}} \boldsymbol{A}\left(\beta_{i, k} \boldsymbol{M}_{i k}\right)=1 .$

(35)

进一步可得：

$\beta_{i, k}^{2}\left(\sqrt{2 / 3} \sigma_{\mathrm{s}}\right)^{2} \boldsymbol{\delta}_{k}^{\mathrm{T}} \boldsymbol{B}_{i}^{\mathrm{T}} \boldsymbol{D}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{D} \boldsymbol{B}_{i} \boldsymbol{\delta}_{k}=\left(\hat{X}_{i}^{k}\right)^{2} .$

(36)

由于β_{i, k}>0，故：

$\beta_{i, k}=\hat{X}_{i}^{k} /\left(\sqrt{2 / 3} M_{\mathrm{P}} \sqrt{\boldsymbol{\delta}_{k}^{\mathrm{T}}} \boldsymbol{B}_{i}^{\mathrm{T}} \boldsymbol{D}^{\mathrm{T}} \boldsymbol{A D B}_{i} \boldsymbol{\delta}_{k}\right).$

(37)

取每个积分点上得到的β_{i, k}的最小者，即为第k次迭代时薄板极限下限载荷乘子：

$v_{k}^{-}=\min \limits_{i \in {\rm{IC}}} \beta_{i, k}.$

(38)

3.3 迭代求解步骤根据上述薄板极限上、下限分析的求解过程，构造如下的迭代求解步骤。
第0次迭代：在迭代开始时认为整个薄板结构处于理想塑性状态，令$\hat{X}_{i}^{0}=1$。由式(16)得H_i⁰=$\sqrt{2 / 3} M_{\mathrm{P}} \rho_{i}|\boldsymbol{J}|{ }_{i} \boldsymbol{K}_{i}$，求解如下初始的二次极小化数学规划问题：

$v=\min \sum\limits_{i \in \operatorname{IG}} \dot{\boldsymbol{a}}^{\mathrm{T}} \boldsymbol{H}_{i}^{0} \dot{\boldsymbol{a}},$

(39)

$\text { s. } \mathrm{t.} \ \ \boldsymbol{F}^{\mathrm{T}} \dot{\boldsymbol{a}}=1.$

(40)

按照式(22)—(24)所示的处理过程，可依次得到δ₀和$\dot{\boldsymbol { a }}$₀，再分别由式(25)和(38)，可以计算出初始迭代时薄板结构极限上、下限载荷乘子：

$v_{0}^{+}=\sqrt{2 / 3} M_{\mathrm{P}} \sum\limits_{i \in \operatorname{IG}} \rho_{i}|\boldsymbol{J}|_{i} \sqrt{\dot{\boldsymbol{a}}_{0}^{\mathrm{T}} \boldsymbol{K}_{i} \dot{\boldsymbol{a}}_{0}},$

(41)

$v_{0}^{-}=\min \limits_{i \in \operatorname{IG}} \beta_{i, 0}.$

(42)

第k(k≥1)次迭代：根据第(k-1)次迭代的计算结果$\dot{\boldsymbol { a }}$_k－1，先由式(17)—(19)所给的划分准则确定出刚性区子集R_k和塑性区子集P_k，并由此计算出每个积分点上$\hat{X}_{i}^{k}$和H_i^k的值，再求解如下的二次极小化数学规划问题：

$\nu=\min \sum\limits_{i \in \operatorname{IG}} \dot{\boldsymbol{a}}^{\mathrm{T}} \boldsymbol{H}_{i}^{k} \dot{\boldsymbol{a}}, $

(43)

$\text { s. t. } \ \boldsymbol{F}^{\mathrm{T}} \dot{\boldsymbol{a}}=1.$

(44)

类似于上述求解过程计算得到δ_k和$\dot{\boldsymbol { a }}$_k，相应地可得到本次迭代薄板的极限上、下限载荷乘子：

$v_{k}^{+}=\sqrt{2 / 3} M_{\mathrm{P}} \sum\limits_{i \in \operatorname{IG}} \rho_{i}\ |\boldsymbol{J}|_{i}\ \sqrt{\dot{\boldsymbol{a}}_{k}^{\mathrm{T}} \boldsymbol{K}_{i} \dot{\boldsymbol{a}}_{k}}, $

(45)

$v_{k}^{-}=\min \limits_{i \in \operatorname{IG}} \beta_{i, k}.$

(46)

当上述迭代过程满足以下收敛条件时，迭代终止：

$\begin{array}{c}\left\|\dot{\boldsymbol{a}}_{k}-\dot{\boldsymbol{a}}_{k-1}\right\| /\left\|\dot{\boldsymbol{a}}_{k-1}\right\| \leqslant \mathrm{vol} 1, \\\left|v_{k}^{+}-v_{k-1}^{+}\right| / v_{k-1}^{+} \leqslant \mathrm{vol} 2.\end{array}$

(47)

这里vol1和vol2分别是根据计算所需精度确定的误差容限。
4 数值算例按照上述薄板结构极限上、下限分析的迭代求解步骤，本文在对Zhou等^[18]的C¹ NEM极限上限分析程序改进优化之后，进一步采用C¹ NEM编制了求解圆形、菱形、等边多边形薄板极限上、下限载荷的计算程序。这些程序中所取的各项参数如下：屈服应力σ_s=200 MPa，Young模量E=2.1×10⁵ MPa，Poisson比v=0.3，误差容限vol1=vol2=1×10^－4。根据Zhou等^[18]在采用C¹ NEM进行薄板结构极限上限分析时关于积分方案选择的详细讨论和分析情况，本文在每个Delaunay三角形上采用3点积分法则来进行C¹自然单元法的数值积分。采用移动最小二乘技术^[17-18]，由积分点上的Mises广义等效应力拟合出节点上的Mises广义等效应力，直观显示出薄板结构在塑性极限状态下的失效云图。
4.1 圆形薄板本算例对受横向集中力P、周边简支和均布压力q、周边夹支作用的圆板进行研究。该圆板受集中力P作用、周边夹支时的示意图如图 3所示。计算中取圆板的半径为R=1.0 m、厚度为h=0.01 m，分别采用Zhou等^[18]中的457和801个节点对该圆板进行离散，见图 4。

图 3 受集中力作用的夹支圆形薄板

图选项

图 4 圆形薄板的节点布置[18]

图选项

表 1给出了采用本文求解程序得到的数值解，同时也列出了徐秉业和刘信声^[2]的理论解、以及张丕辛^[23]、Zhou等^[18]的上限数值解。可见：1) 本文在集中力、周边简支和Zhou等^[18]在均布压力、周边夹支时的上限解均略小于张丕辛^[23]的上限解，本文获得的下限解均与相应的上限解非常接近; 2) 整体而言，本文的数值解与徐秉业和刘信声^[2]的理论解吻合较好。在集中力、周边简支时，采用2种节点获得的上、下限解均略大于理论解; 在均布力、周边夹支时，采用2种节点获得的下限解均略小于理论解。上述结果的微小数值差异主要与节点布置和数值计算误差有关。
表 1 圆形薄板的极限载荷乘子

工况	徐秉业和刘信声[2]	张丕辛[23]上限	本文(上限/下限)
			457个节点	801个节点
集中力，简支(Mp/P)	6.283	6.971	6.673/6.550	6.612/6.453
均布力，夹支(Mp/qR2)	12.5	12.56	12.541[18]/12.317	12.486[18]/12.043

表选项






上述2种工况下圆形薄板极限载荷乘子的迭代收敛曲线(801个节点)如图 5所示。可以看出：上限载荷乘子单调下降; 本文的计算程序只需8~9次迭代，所得到的极限上、下限载荷乘子就能趋于稳定并且收敛到预设的误差容限。

图 5 圆形板极限载荷乘子的迭代收敛曲线(801个节点)

图选项

4.2 菱形薄板本算例研究在简支和固支边界约束、受均布压力载荷q作用的菱形薄板。该板的几何参数和节点布置如图 6所示。在本文的数值计算中，板的厚度取为h=0.02 m，图 6a中所示的半径取为R=0.5m，变化倾斜角α分别取为30°、45°、60°、75°、90°。针对该菱形板，Capsoni和Vicente da Silva^[13]和Nguyen-Thoi等^[16]分别采用Mises屈服准则和不同的Mindlin板单元，获得了该菱形板的上限解，相应的数值结果如表 2所示。从该表可以看出，本文得到的下限解比文[13, 16, 18]中的上限解略小、均与其非常接近。

图 6 菱形薄板

图选项

表 2 受均布载荷作用菱形薄板的极限载荷乘子(M_P/(qL²))

α/(°)	简支				夹支
	文[13]上限	文[18]上限	本文下限		文[13]上限	文[16]上限	文[18]上限	本文下限
30	5.140	4.808	4.743		9.575	9.831	8.901	8.700
45	5.609	5.447	5.393		10.596	11.021	9.939	9.751
60	5.966	5.916	5.853		11.394	11.898	10.781	10.550
75	6.197	6.186	6.108		11.893	12.438	11.065	10.800
90	6.278	6.267	6.198		12.062	12.636	11.296	11.034

表选项






在倾斜角α取90°时，上述菱形板就变为正方形板。许多****^{[6, 8-11, 13, 16]}采用不同的方法研究了受均布载荷作用方形薄板的极限承载能力。表 3给出了上述****获得的方形薄板的极限载荷乘子。为了便于比较，表 3对不同****获得的极限载荷乘子进行了统一转换，例如6.198M_P/(qR²)≈24.79M_P/(qL²)。同时，为了比较不同数值方法的计算精度，本文也采用非协调矩形薄板单元ACM，对该方形板进行了极限分析，计算结果也列于表 3中。从该表可以看出，本文采用C¹ NEM得到极限上、下限解均与不同****的参考解非常吻合，且与采用ACM获得的数值解之间的差异非常小。
表 3 受均布载荷作用方形板极限载荷乘子的比较(M_P/(qL²))

作者	方法	简支	夹支
Le等[11]	上限	25.01	45.07
Capsoni和Corradi[9]	上限	25.02	45.29
Nguyen-Thoi等[16]	上限	25.06	50.55
本文(C1 NEM)	上限/下限	25.07[18]/24.79	45.18[18]/44.14
本文(ACM)	上限/下限	25.00/24.73	46.35/45.64
Capsoni和Vicente da Silva[13]	上限	25.11	48.25
Hodge和Belytschko[8]	上限	26.54	49.25
Lubliner[6]	上限	27.71	52.01

表选项






数值算法的迭代次数和计算时间是反映其计算效率的重要指标。迭代次数一般与算法本身性能、节点数目等因素有关，除了上述因素外，计算时间还与计算机性能、程序编制技巧等因素有关。表 4给出了采用不同方法对方形板进行极限分析的计算效率比较。可见：1)本文采用C¹ NEM和ACM对简支方板进行极限分析的迭代次数均为19次; 采用C¹ NEM对夹支方板进行极限分析的迭代次数为8次，介于Le等^[11]采用EFG程序的6次和本文ACM程序的11次之间。2)统计自同一台计算机(Intel(R) Core (TM) i5-3470，主频3.2 GHz，内存3.47 GB)的计算时间表明，基于C¹ NEM极限分析程序的计算时间比基于ACM的略少。需要指出的是^[18]：C¹ NEM比ACM使用更加灵活，不仅可以用来分析矩形板，也可用来分析其他形状的板; 与Le等^[11]使用的EFG程序相比，本文的C¹ NEM极限分析程序可以像ACM和其他有限板单元一样容易、精确地施加位移边界条件。
表 4 方形板极限分析的计算效率比较

作者	简支			夹支
	迭代次数	计算时间/s		迭代次数	计算时间/s
Le等[11]	—	—		6	—
本文(C1 NEM)	19	202.27		8	91.02
本文(ACM)	19	204.92		11	124.20

表选项






基于极限上限分析得到的塑性耗散功率云图和下限分析得到的Mises广义等效应力云图是对结构塑性失效模式的直观反映。图 7给出了该菱形薄板在极限状态下倾斜角为α=60°和α=90°、边界条件为简支和夹支时的Mises广义等效应力云图。这些云图与Capsoni和Vicente da Silva^[13]、Nguyen-Thoi等^[16]和Zhou等^[18]提供的、对应的塑性耗散功率分布云图具有良好的一致性。这说明了本文方法的有效性和合理性。

图 7 (网络版彩图)菱形薄板极限状态下的Mises广义等效应力(MPa)

图选项

图 8和9分别给出了本文和Zhou等^[18]在简支和夹支约束、采用5种不同倾斜角时得到的该菱形薄板极限载荷乘子的迭代收敛曲线。可以看出：在简支和夹支约束时，分别只需19~22次和8~12次迭代，得到的极限上、下限载荷乘子就能分别收敛于一系列稳定的极值。

图 8 简支菱形薄板极限载荷乘子的迭代收敛曲线

图选项

图 9 夹支菱形薄板极限载荷乘子的迭代收敛曲线

图选项

4.3 等边多边形薄板本算例研究受均布载荷q作用的等边多边形薄板，如图 10a和10b所示。计算中取板的外接圆半径为R=0.5 m、厚度为h=0.02 m，采用图 10c和10d所示的节点布置对该组等边多边形薄板进行离散。

图 10 等边多边形板

图选项

Wood^[24]也对该组受均布载荷作用的等边多边形薄板进行了上限分析。表 5给出了本文的数值解和参考文献中的上限解。从表中可以看出，本文得到的数值解与参考解很吻合，夹支时的上、下限解介于Wood^[24]的上限解之间。综合表 3和表 5可以看出，在相同载荷作用下，板结构在夹支边界时的极限载荷略小于简支边界时的2倍。
表 5 等边多边形薄板的极限载荷乘子(M_p/(qR²))

等边多边形板	简支		夹支
	本文上限/下限		本文上限/下限	Wood[24]上限
等边三边形	5.545[18]/5.479		10.039[18]/9.821	9.999
等边五边形	6.180/6.056		11.566/11.231	11.358

表选项






图 11给出了简支和夹支等边三边形板在极限状态下的Mises广义等效应力云图。该云图是对称和合理的，与Zhou等^[18]得到的塑性耗散功率分布云图吻合良好。

图 11 (网络版彩图)等边三边形薄板在极限状态下的Mises广义等效应力(MPa)

图选项

等边多边方形薄板极限载荷乘子的迭代收敛曲线如图 12所示。综合图 8、图 9和图 12的迭代收敛曲线可见，当误差容限相同时，本文的迭代求解算法在板结构夹支时所需的迭代步数比简支时的少。

图 12 等边多边方形薄板极限载荷乘子的迭代收敛曲线

图选项

5 结论本文根据薄板塑性极限上限分析的迭代求解格式，构造出了满足平衡方程和边界条件的广义应力场，并将该广义应力场代入极限下限定理，得到了与极限上限载荷乘子相对应的、易于程序实现的下限载荷乘子的迭代计算格式。充分发挥C¹ NEM计算精度高、便于施加位移边界条件、使用灵活等优势，基于该方法开发了求解程序，获得了不同工况下圆形、菱形、等边多边形薄板的极限上、下限载荷乘子。同时，采用移动最小二乘技术，拟合得到了能够直观反映薄板结构塑性失效模式的Mises广义等效应力云图。计算结果表明：
1) 得到的极限上、下限载荷乘子与参考文献中的解析解和数值解吻合良好，具有较高的计算精度。
2) 得到的Mises广义等效应力云图是对称合理的，与相同工况下参考文献得到的塑性耗散功率云图具有很好的一致性，验证了本文构造出的薄板结构极限下限迭代格式的有效性。
3) 迭代算法具有较快的收敛速度，简支时只需20次左右、夹支时只需10次左右迭代就能在指定误差容限内获得本文薄板的极限上、下限载荷乘子。
鉴于本文提出的板结构上、下限迭代求解方法的有效性、精确性和适用性，可将其与其他适用于中厚板结构的单元或数值方法相结合，进一步研究任意形状(圆板、斜板、多边形板、不规则板)中厚板结构的极限承载能力。

参考文献

[1]	陈钢, 刘应华. 结构塑性极限与安定分析理论及工程方法[M]. 北京: 科学出版社, 2006. CHEN G, LIU Y H. Numerical theories and engineering methods for structural limit and shakedown analyses[M]. Beijing: Science Press, 2006. (in Chinese)
[2]	徐秉业, 刘信声. 结构塑性极限分析[M]. 北京: 中国建筑工业出版社, 1985. XU B Y, LIU X S. Plastic limit analysis of structures[M]. Beijing: China Architecture & Building Press, 1985. (in Chinese)
[3]	YU M H, MA G W, LI J C. Structural plasticity: Limit, shakedown and dynamic plastic analyses of structures[M]. Berlin: Springer Press, 2009.
[4]	HODGE P G. Limit analysis of rotationally symmetric plates and shells[M]. New Jersey: Prentice-Hall, 1963.
[5]	SAVE M A, MASSONNET C C. Plastic analysis and design of plates, shells, and disks[M]. Amsterdam: Elsevier Science Ltd., 1972.
[6]	LUBLINER J. Plasticity theory[M]. New York: Macmillan, 1990.
[7]	FOX E N. Limit analysis for plates: The exact solution for a clamped square plate of isotropic homogeneous material obeying the square yield criteron and loaded by uniform pressure[J]. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1974, 277(1265): 121-155.
[8]	HODGE P G Jr, BELYTSCHKO T. Numerical methods for the limit analysis of plates[J]. Journal of Applied Mechanics, 1968, 35(4): 795-802.
[9]	CAPSONI A, CORRADI L. Limit analysis of plates: A finite element formulation[J]. Structural Engineering and Mechanics, 1999, 8(4): 325-341. DOI:10.12989/sem.1999.8.4.325
[10]	TURCO E, CARACCIOLO P. Elasto-plastic analysis of Kirchhoff plates by high simplicity finite elements[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 190(5-7): 691-706. DOI:10.1016/S0045-7825(99)00438-7
[11]	LE C V, GILBERT M, ASKES H. Limit analysis of plates using the EFG method and second-order cone programming[J]. International Journal for Numerical Methods in Engineering, 2009, 78(13): 1532-1552. DOI:10.1002/nme.2535
[12]	BLEYER J, DE BUHAN P. On the performance of non-conforming finite elements for the upper bound limit analysis of plates[J]. International Journal for Numerical Methods in Engineering, 2013, 94(3): 308-330. DOI:10.1002/nme.4460
[13]	CAPSONI A, DA SILVA M V. A finite element formulation of Mindlin plates for limit analysis[J]. International Journal for Numerical Methods in Biomedical Engineering, 2011, 27(1): 143-156. DOI:10.1002/cnm.1300
[14]	BAGHANI M, FEREIDOONNEZHAD B. Limit analysis of FGM circular plates subjected to arbitrary rotational symmetric loads using von-Mises yield criterion[J]. Acta Mechanica, 2013, 224(8): 1601-1608. DOI:10.1007/s00707-013-0828-z
[15]	MAKRODIMOPOULOS A. An upper bound limit analysis formulation for thin plates[J]. Engineering and Computational Mechanics, 2015, 168(4): 133-143.
[16]	NGUYEN-THOI T, PHUNG-VAN P, NGUYEN-THOI M H, et al. An upper-bound limit analysis of Mindlin plates using CS-DSG3 method and second-order cone programming[J]. Journal of Computational and Applied Mathematics, 2015, 281: 32-48. DOI:10.1016/j.cam.2014.12.006
[17]	周书涛, 刘应华, 陈莘莘. 基于非协调矩形弯曲单元的薄板极限上限分析方法[J]. 清华大学学报(自然科学版), 2011, 51(12): 1887-1893. ZHOU S T, LIU Y H, CHEN S S. Upper-bound limit analysis method of thin plates based on the nonconforming rectangular bending element[J]. Journal of Tsinghua University (Science and Technology), 2011, 51(12): 1887-1893. (in Chinese)
[18]	ZHOU S T, LIU Y H, CHEN S S. Upper bound limit analysis of plates utilizing the C1 natural element method[J]. Computational Mechanics, 2012, 50(5): 543-561. DOI:10.1007/s00466-012-0688-8
[19]	陈浩峰. 多组载荷下结构极限分析和参考应力确定的数值方法及工程应用[D]. 北京: 清华大学, 1998. CHEN H F. Numerical methods for limit analysis and reference stress determination of structures under multi-loading systems and their engineering applications[D]. Beijing: Tsinghua University, 1998. (in Chinese)
[20]	ZHOU S T, LI Y, LIU Y H, et al. A numerical estimate method for limit analysis of plane structures by adopting the natural element method[C]//The 22nd Conference on Structural Mechanics in Reactor Technology. California, USA, 2013: 1279-1287.
[21]	SUKUMAR N, MORAN B. C1 natural neighbour interpolant for partial differential equations[J]. Numerical Methods for Partial Differential Equations, 1999, 15(4): 417-447. DOI:10.1002/(SICI)1098-2426(199907)15:4<417::AID-NUM2>3.0.CO;2-S
[22]	SUKUMAR N, MORAN B, BELYTSCHKO T. The natural element method in solid mechanics[J]. International Journal for Numerical Methods in Engineering, 1998, 43(5): 839-887. DOI:10.1002/(SICI)1097-0207(19981115)43:5<839::AID-NME423>3.0.CO;2-R
[23]	张丕辛. 极限分析的无搜索数学规划算法及其应用[D]. 北京: 清华大学, 1989. ZHANG P X. No-search mathematical programming algorithm for limit analysis and its engineering applications[D]. Beijing: Tsinghua University, 1989. (in Chinese)
[24]	MANSFIELD E H. Studies in collapse analysis of rigid-plastic plates with a square yield diagram[J]. Proceedings of the Royal Society A-Mathematical, and Physical and Engineering Sciences, 1957, 241(1226): 311-338.