0.   引　言

在船舶结构设计中，为减轻结构自身重量，便于设备、管线和人员通行，在结构强度允许的前提下，会使用大量开孔梁和开孔板结构。在进行船舶修理时，如需拆换船舶内部设备，在大型设备需要出入舱的过程中，有时也需要在结构上开设大开口。这些孔洞的存在，会削弱结构自身强度，造成局部的应力集中甚至是结构失稳。因此，有必要对开孔结构重新进行强度校核。

船舶结构上的孔洞除会削弱结构强度外，还容易产生局部屈曲效应。如在高剪切应力作用下，大部分开有不规则孔洞的板和梁（如腹板和面板）会因孔洞的原因而受到较高的集中载荷，可能会产生Web-Post屈曲、Vierendeel Bending和破坏^[1-2]，从而导致明显不同的结构载荷响应和屈曲破坏力^[3-5]。对船舶开孔结构的校核，各大船级社都有相关规定。针对这些开孔板和开孔梁的设计及评估方法依旧局限于BS5950：Part 1^[6]和Eurocode 3^[7]，这些规则与标准考虑到工程安全，大部分是通过逐步校核梁、板的稳性及横截面临界能力来进行评判的，方法通常偏于保守。针对具有一定几何形状（包括布局和尺寸范围）的有限域的问题，很大程度上是基于简化的分析方法^{[2, 8-9]}，而这些模型计算强度一般较大，在实践中较难推广。

在船舶结构分析的传统方法中，有限元法是最常用也是最为有效的计算方法，其在处理船舶结构变形、应力和疲劳分析等方面都有着成熟和广泛的应用^[10-12]，但在分析船舶结构场量（如位移场和应力场等）剧烈变化的高梯度区域时，会出现计算精度降低甚至计算中断等现象。通常，采用加密该区域网格或者采用高阶单元来解决此问题。但是，这给有限元法前、后处理带来了巨大的工作量，导致计算效率降低，同时，这种方法并不能从根本上消除问题的根源。

与有限元法原理相对应的另一种数值分析方法是无网格法，该方法在近20年中也得到了很大的发展^[13-16]。无网格法是建立在系列独立的离散节点上，通过构造这些离散点的近似函数来对问题域进行求解。无网格法在计算过程中可以根据需要任意增加节点，不需要处理节点之间的拓扑信息，特别适合自适应分析计算。同时，由于无网格法的节点可以由计算机以自主的方式进行创建，故可节省花费在创建和处理网格上的时间和计算资源。无网格法目前在航空材料、高速碰撞、动态裂纹扩展、加工成型等诸多领域得到了较为成功的应用^[17-19]。

鉴于上述背景，基于无网格局部Petrov-Galerkin（MLPG）法，拟提出开孔结构局部屈曲问题的无网格数值分析方法。首先，基于MLPG法定义结构的问题域并进行离散，在离散节点上建立其紧支函数，采用加权余量法建立系统位移场（应力场）的离散方程；然后，运用移动最小二乘法（MLS）逼近离散节点的形函数；最后，通过算例评估该方法对开孔梁板弹性屈曲问题分析的有效性和正确性。

1.   基于MLPG的梁屈曲问题

在目前的研究中，分析板屈曲问题的一般计算公式是在切线刚度矩阵（${{{\mathit{\boldsymbol{K}}}}_{\rm{T}}}$）奇异性的基础上发展起来的，一般由问题域材料刚度矩阵（${{{\mathit{\boldsymbol{K}}}}_{\rm{E}}}$）和几何刚度矩阵（${{{\mathit{\boldsymbol{K}}}}_{\rm{G}}}$）组成。K_E只是单纯基于Kirchhoff板理论，K_G的计算仅仅需要一阶运动学方程的弹簧旋转类比（RSA）原理类推简化而来^[20]。在预先给定的确定平面应力的线性分析中，容易获得这些参数。作为大多数弹性结构问题，K_E可以被认为是独立的内部力量。这种方法虽然考虑了材料的型线响应，但在发生变形屈曲之前，通常会忽略小的变形，从而把屈曲问题作为一个线性特征值问题，直接得出其简化计算公式。

本文采用MLPG法来离散问题域，然后使用MLS构建问题域位移场的近似函数。在建立开孔梁的严格屈曲模式过程中，文章采用假定模式与互补模式结合的方法，组成一个降低特征值求解的二阶问题，然后通过迭代的方法^[21]使多自由度问题能够得到有效解决。在迭代运算中，局部域的应用在求解精度和计算要求方面对两者有很大的帮助^[22-24]。

1.1   位移场逼近函数

不同于有限元法所采用的分段函数（位移函数的形函数）预定义在所有单元子域上，MLPG法位移场函数是通过MLS获得的。这种方法的主要特征是，基于某一问题域${{\mathit{\boldsymbol{\varOmega }}}}$内一系列给定节点的曲线拟合来获得一个连续逼近函数，采用多项式组合来完成：

式中：${u}^{\rm{h}}\left( x \right)$为求解区域场函数的近似函数；${{{\mathit{\boldsymbol{p}}}}}^{\rm{T}}\left( x \right)$为多项式基函数；${{{\mathit{\boldsymbol{a}}}}}\left( x \right)$为待定系数。对于$x y$平面内的问题，${{{\mathit{\boldsymbol{p}}}}^{\rm{T}}}\left( x \right) = [ 1\;x\;y\;{{x^2}}\;{xy}\;{{y^2}} ]$，可以看做为一个二次基函数。参考文献[21]，采用MLS方法规定的节点参数，即节点位置x_I和节点位移u_I (I=1, 2,···, N)来构造近似函数；${{{\mathit{\boldsymbol{a}}}}}\left( x \right)$可以通过加权L2范数最小化来建立，可以表示为

式中：$w\left( {{{\mathit{\boldsymbol{x}}}} - {{{\mathit{\boldsymbol{x}}}}_I}} \right)$为权函数，${{{\mathit{\boldsymbol{x}}}}} = {[x\;y ]^{\rm{T}}}$为空间坐标系内的向量；${{\mathit{\boldsymbol{\varDelta }}}}$为节点参数值的加权残值。

泛函${{\mathit{\boldsymbol{\varDelta }}}}$中取极小值，可得

其中

定义${{\mathit{\boldsymbol{\varPhi }}}}\left( x \right)$为由节点参数得到的场逼近函数形函数矩阵，设${{\mathit{\boldsymbol{\varPhi }}}}\left( x \right)$为

由此，式（1）可进一步简化为

权函数在MLS中具有非常重要的作用，它可以用来控制和调节影响域的位置，一般要求其具有紧支性。本文采用三次样条函数W(r)来表达权函数 ：

式中：$r = \left| {{{\mathit{\boldsymbol{x}}}} - {{{\mathit{\boldsymbol{x}}}}_I}} \right|$。

1.2   梁平面响应

为了求得载荷梁的外平面几何刚度矩阵${{{\mathit{\boldsymbol{K}}}}_{\rm{G}}}$，需要得知与平面响应相关的应力分布。在使用MLPG法离散整个问题域时，有必要简化离散过程，这样对于那些具有相同几何形状的单元来说，可以避免重复计算。本文提出将开孔梁分割成梁单元以进行简化，其中每个单元便构成一个被降低了自由度的超级单元。图1所示为一个简化单元，共有4个质心超级节点，节点位置如图所示，每个节点有3个平面自由度$\left( {{u}, {v}, {\theta }} \right)$，每个单元受到集中载荷F和力矩M的作用。在不规则开口情况下，需要使用更多的单元模型来表示每个不同的开口。

图  1  简化单元力学模型

Figure  1.  Simplified element mechanical model

下载: 全尺寸图片 幻灯片

通过MLPG法对单个单元离散建立起其特征单元响应，然后再利用一个标准的离散装配过程^[23]，把所有单元组装成整体的梁响应。特征单元的响应通过考虑单元力交变载荷的情况获得。在节点1处，当单元在外载荷作用下达到平衡时，节点2，3，4每一个单元对应于一个超级单元自由度，如图1所示。按照MLPG方法逼近，这为转化4个节点超级单元平面刚度去替代单个单元响应提供了一种灵活的方法。

由于是基于单元的全梁分析，得到的应力仅在局部单元域连续，故其产生的相关误差对于典型的梁结构来说可以忽略不计。此外，一个连续平面应力场可以在整体梁内通过MLS来近似实现。单个单元平面应力场的计算公式如下：

式中：${{\mathit{\boldsymbol{D}}}}$为平面应力弹性矩阵；ε 为应变；${{{\mathit{\boldsymbol{B}}}}_{xy}}$为应变位移矩阵；${{{\mathit{\boldsymbol{T}}}}_{xy}}$为在MLPG单元域内节点参数对应超级单元离散载荷和分布载荷的转换矩阵；F为集中载荷 ；ω为y方向的均布载荷。式（9）封装了平面离散分析的2个层次，即MLPG和超级单元。

1.3   平面外响应

由文献[20]可知，平面外响应的切线刚度矩阵${{{\mathit{\boldsymbol{K}}}}_{\rm{T}}}$是由${{{\mathit{\boldsymbol{K}}}}_{\rm{E}}}$和${{{\mathit{\boldsymbol{K}}}}_{\rm{G}}}$这2个部分组成的。

1.3.1   材料刚度矩阵

根据Kirchhoff薄板理论^[22]，运用MLPG法，${{{\mathit{\boldsymbol{K}}}}_{\rm{E}}}$可以分解为弯曲刚度矩阵${{\mathit{\boldsymbol{K}}}}$和施加本质边界条件后得到的罚函数刚度矩阵${{{\mathit{\boldsymbol{K}}}}_\alpha }$：

式中：${E}$为杨氏模量；$\mu$为泊松比；${t_{\rm{w}}}$为梁腹板厚度；形函数矩阵${{\mathit{\boldsymbol{\varPhi }}}}\left( x \right)$的下标表示对不同空间的导数。

进一步地，罚函数刚度矩阵${{{\mathit{\boldsymbol{K}}}}_{\rm{\alpha}} }$可定义为

式中：${\varGamma _u}$为规定的位移边界；${{{\mathit{\boldsymbol{R}}}}_{xy}}$为连接节点位移和节点旋转自由度的转换矩阵；罚因子${\alpha _{\textit{z}}}$和${\alpha _\theta }$为常数，通常其范围为刚度矩阵最大对角元素值的倍数，位数取值范围为10⁴~10^13[23]。从施加约束出发，罚因子应为无穷大，但在实际计算中不可能实现。在实际计算中，如果罚因子取值过小，可能会造成约束不能精确施加；如果取值过大，则计算可能会溢出或者得到病态解。为简化计算，本文罚因子的取值为弯曲刚度矩阵${{\mathit{\boldsymbol{K}}}}$最大对角线元素值的倍数，且使得${{\rm{\alpha}} _{\textit{z}}} = {{\rm{\alpha}} _{\rm{\theta}} }$。

1.3.2   几何刚度矩阵

几何刚度${{{\mathit{\boldsymbol{K}}}}_{\rm{G}}}$是根据RSA方法从平面应力场${{\rm{\sigma}} _x}$，${{\rm{\sigma}} _y}$和${{\rm{\tau}} _{xy}}$得到^[20]。使用相同的离散化方法对${{{\mathit{\boldsymbol{K}}}}_{\rm{G}}}$进行求解，得

式中：${{\mathit{\boldsymbol{\varPhi }}}}{\left( x \right)_{,x}}$和${{\mathit{\boldsymbol{\varPhi }}}}{\left( x \right)_{,y}}$为由MLS法求得的形函数的一阶偏导数；${{\rm{\sigma}} _x}$，${{\rm{\sigma}} _y}$为x，y方向的应力；${{\rm{\tau}} _{xy}}$为xy平面的剪应力。

1.3.3   折边刚度矩阵

折边对开孔梁的局部屈曲响应的影响可以通过一种简化模型来实现。对于顶部和底部都有折边的组合工字梁来说，如果这些折边能够沿着连接腹板内边缘起到完全约束作用的话，那么通常可以只考虑它们的材料和几何刚度；如果有上、下板约束，则可以把板结构进行等效处理从而简化为梁结构。相对于建立三维问题模型来说，上面提到的简化方法仅需采用一维模型即可达到简化折边响应的目的。同时，这种方法还保留了二维平面模型的适用性。

考虑一个一维折边${\varGamma _{\rm{f}}}$，其尺寸为${t_{\rm{f}}} \times {b_{\rm{f}}}$，其中${t_{\rm{f}}}$为折边厚度，${b_{\rm{f}}}$为折边宽度。在忽略次要变形的影响下，折边材料刚度贡献来自其扭转刚度${{\mathit{\boldsymbol{GJ}}}}$，则折边的材料刚度矩阵为

式中，${{{\mathit{\boldsymbol{\varPhi }}}}_{,xy}}$为关于梁的横向位移形函数的二阶偏导数。

另一方面，为了建立相应的几何刚度矩阵，假定沿着折边宽度和厚度方向的应力${{\rm{\sigma}} _x}$为一个常数，在沿着$x$轴方向的内部压缩或拉伸应力变量${{{\sigma}} _x}$的作用下，可以认为该梁为一个一维等效梁。因此，分布在横截面上每单位面积的等效转动几何刚度可以取为${{{\mathit{\boldsymbol{K}}}}_{\rm{\theta}} } = {{{\sigma }}_x}$，等效弯曲的侧向旋转角T_θ可以通过以下公式求得^[21]。

式中：${{{\mathit{\boldsymbol{u}}}}_{\textit{z}}}$为腹板平面平移自由度。因此，根据RSA法，可得折边的几何刚度为

式中，${I_O}$为面积二次极矩。

1.4   简化后的屈曲分析

目前，针对开孔梁的失稳评估模型一般是在仅取梁的某一适当部分的基础上建立起来的。本文采用局部模型，沿梁的长度方向取至少4个单元，这样具有更好的分析效果。本文所提方法只需要在无网格MLPG区域内计算材料和几何刚度矩阵即可，大幅减少了刚度矩阵的规模，提高了计算效率。在梁单元局部域内，在按照比例施加载荷情况下梁的临界屈曲载荷因子${{\rm{\lambda}} _{{\rm{cr}}}}$可以由相关的${{{\mathit{\boldsymbol{K}}}}_{\rm{E}}}$和${{{\mathit{\boldsymbol{K}}}}_{\rm{G}}}$的特征值求得，然而相关矩阵的大小导致求解其特征值的工作量庞大，大幅增加了计算成本。考虑到RSA方法^[20]，当梁在屈曲模态时，${{\rm{\lambda}} _{{\rm{cr}}}}$也可以作为内部能量吸收和转动等效转化的比例系数。本文所提方法采用实际屈曲载荷因子的近似模态代替精确模态的束缚，典型的模态可以近似地通过对线性结构执行屈曲模态来获得^[20-21]。

考虑一个近似平面外模态${{\mathit{\boldsymbol{u}}}}_{\textit{z}}$由对应于外部的平面力${{\mathit{\boldsymbol{f}}}}_{\textit{z}}$获得，忽略屈曲前平面内变形的影响，则临界屈曲载荷因子${\lambda _{{\rm{cr}}}}$可由下式确定。^{[22 ]}

基于假定的屈曲模态，对预定的模式通过进行迭代运算来进一步改进，式（21）可以用来获得一个初始解的一个近似临界屈曲载荷因子${{\rm{\lambda}} _{{\rm{cr}}}}$。本文方法是通过局部域的手段来降低问题的维数，从而使得假定的模态阶数降低，这样，特征值的求解就变得容易^[24-26]。根据这种方法，相应于一个任意的平面外载荷模式${{{{{\mathit{\boldsymbol{f}}}}}}_{{\textit{z}}{\rm{A}}}}$，最初的近似模态${{{{{\mathit{\boldsymbol{u}}}}}}_{{\textit{z}}{\rm{A}}}}$在一个新的负载模式${{{{{\mathit{\boldsymbol{f}}}}}}_{{\textit{z}}{\rm{B}}}}$作用下可以建立一个互补的模态${{{{{\mathit{\boldsymbol{u}}}}}}_{{\textit{z}}{\rm{B}}}}$。可以定义为：

因此，2个互补模态${{{{{\mathit{\boldsymbol{u}}}}}}_{{\textit{z}}{\rm{A}}}}$和${{{{{\mathit{\boldsymbol{u}}}}}}_{{\textit{z}}{\rm{B}}}}$用于计算和解答经降阶了的模态特征值，结果经过每一次迭代计算直至收敛。

2.   数值模拟算例

算例采用一连续多孔板。在板沿长度方向，等距开有半径$R = 250\;{\rm{mm}}$的系列圆形孔，对称半长$L = 40\;000\;{\rm{mm}}$，宽度$B = 2\;000\;{\rm{mm}}$，板厚$t = 5\;{\rm{mm}}$。这种结构在船舶结构中非常普遍，例如肋板、船底纵桁等，这些结构的一般特征就是开有一定数量的减轻孔或者人孔。图2所示为该连续多孔板结构靠近两端的带有4.5个圆形孔的部分简化形式，图中数值的单位为mm。算例中用到的计算常量包括：材料的杨氏弹性模量${{E}} =$$3.0 \times {10^5}\;{\rm{MPa}}$，泊松比$\mu = 0.3$，密度$\rho =$$7.5 \times {10^3}\;{\rm{kg/}}{{\rm{m}}^{\rm{3}}}$。板的上、下端为自由边，左侧边界刚性固定，右侧平直边部分受到与板中性面方向一致的均布压力$\omega = 2 \times {10^4}\;{\rm{N}}$。

图  2  4个孔的单元结构图

Figure  2.  Structural drawing of four hole units

下载: 全尺寸图片 幻灯片

算例采用MLPG无网格算法，初始节点42 108个，如图3（a）所示；采用高斯积分法，积分背景初始网格83 194个，经过五步自适应计算，得到如图3（b）所示的计算结果。

为了进行比较分析，算例采用有限元法（NASTRAN 2016）进行了验证性计算。图4（a）所示为有限元离散网格，共4 900个节点，4 583个网格，计算得到如图4（b）所示的计算结果。

图  3  无网格MLPG法计算结果

Figure  3.  Calculation results of meshless MLPG method

下载: 全尺寸图片 幻灯片

图  4  有限元法计算结果

Figure  4.  Calculation results of finite element method

下载: 全尺寸图片 幻灯片

从2种方法的计算结果来看，计算的应力分布高度相似，并且应力集中区域几乎都落在相同的区域。从无网格方法的计算结果来看，由于采用的是连续应力场函数，不仅可以清晰地得到每一点的应力分布，甚至在那些边缘应力出现间断时也都可以获得。这就是无网格法相对于有限元的重要优势。

为进一步分析文章所提方法对结构局部屈曲载荷因子分析的有效性，对有限元法与本文所提方法的计算结果进行了比较，结果如图5所示。从图中可以看出，有限元法基本上对于整个结构的屈曲载荷因子都不敏感，对于结构当中由孔洞所带来的的局部结构屈曲并没有明显的响应；而本文所提的运用MLPG方法计算得出在结构2个端部的屈曲载荷因子较大，结构中部屈曲载荷因子降低，该结果与原结构的值吻合较好。

图  5  运用MLPG法和有限元法计算的屈曲载荷因子变化图

Figure  5.  Change diagram of buckling load factor calculated by MLPG method and finite element method

下载: 全尺寸图片 幻灯片

为了更好地确认本文所提方法能够更好地分析临界屈曲载荷因子大小随结构变化的敏感性，即分析结构屈曲载荷因子随结构变化而变化的情况，将原结构两侧第2和第3个圆形孔按照以直径作为边长修改为了正方形（称为改型方案），如图6所示。图中，数值单位为mm。

图  6  改型方案结构示意图

Figure  6.  Structural diagram of modification scheme

下载: 全尺寸图片 幻灯片

运用MLPG无网格算法进行计算，初始节点44 400个，如图7（a）所示；采用高斯积分法，积分背景初始网格为87 923个，经过五步自适应计算，得到如图7（b）所示的计算结果。

图  7  改型算例无网格MLPG法计算结果

Figure  7.  Calculation results of modified example by meshless MLPG method

下载: 全尺寸图片 幻灯片

同样，采用有限元法（使用NASTRAN 2016）进行比较分析。图8（a）所示为有限元离散网格，共3 585个节点，3 127个网格，得到如图8（b）所示的计算结果。

图  8  改型算例有限元法计算结果

Figure  8.  Calculation results of modified example by finite element method

下载: 全尺寸图片 幻灯片

将经过修改的结构再次运用2种方法予以计算。从计算应力云图结果来看，应力分布依然高度相似，并且应力集中趋势也基本一致。重新计算出的屈曲载荷因子如图9所示。有限元法计算结果显示结构修改后的载荷因子与原结构载荷因子几乎相同，说明其对于孔洞的细微修改并不敏感，而本文所提出的方法则非常敏感，两侧第2和第3个孔洞的扩大（由半径$R =$250 mm扩大为边长$a =$500 mm的正方形）明显降低了其所在位置的屈曲载荷因子，同时也对其两侧结构产生了一定的影响，而对于其外的结构基本不产生影响，这个计算结果也与原结构有着更高的契合度。

图  9  运用MLPG法和有限元法计算的改型算例屈曲载荷因子变化图

Figure  9.  Change diagram of modified example buckling load factor calculated by MLPG method and finite element method

下载: 全尺寸图片 幻灯片

3.   结　论

本文采用MLPG法对开孔梁和开孔板的弹性屈曲问题予以了计算和评估。在建立开孔梁的严格屈曲模式过程中，采用假定模式与互补模式相结合的方法，降阶求解了模态特征值。在迭代运算中，局部域的应用在求解精度和计算要求方面均非常有效。

运用有限元法和本文所提方法对算例进行了计算，通过对应力计算结果的分析和比较，得到2种方法计算的应力分布高度相似，应力集中区域基本相同。本文方法由于采用的是连续应力场函数，所以可以清晰地得到场点的应力分布。

通过对算例中结构局部屈曲载荷因子的计算和比较，得出本文所提方法能够较好地仿真结构2个端部的屈曲载荷因子，其结果与原结构吻合较好。

算例还测试了本文方法和有限元法临界屈曲载荷因子的大小随结构变化的敏感性，结果显示传统方法（有限元法）对板梁孔洞的细微修改不敏感，而本文所提方法则非常敏感，该计算结果与原结构的契合度更好。