下载:
-
在船舶结构设计中，考虑到为减轻结构自身重量、方便设备和管线穿过以及人员通行等因素，在结构强度允许的前提下，大量使用开孔梁和开孔板结构。在船舶制造和修理等施工过程中，有时也需要在结构上开设大开口，例如船舶内部结构拆换、大型设备出入舱等。由于这些孔洞的存在，必然削弱结构自身强度，造成局部的应力集中甚至结构失稳。因此有必要对开孔结构重新进行强度校核。图1所示为船舶典型双层底结构。
图 1 船舶典型双层底结构
Figure 1. The typical opening hole
除产生强度削弱问题外，开孔结构还有一个突出问题就是容易产生局部屈曲效应。如在高剪切应力作用下，可能会产生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法（mesh-free local petrov-galerkin method，MLPG），提出开孔结构局部屈曲问题的无网格数值分析方法。MLPG法首先定义结构的问题域并进行离散，然后在离散节点上建立其紧支函数，采用加权余量法建立系统位移场（应力场）的离散方程，最后运用移动最小二乘法（moving least square method，MLS）来逼近离散节点的形函数。本文通过算例来评估提出方法对开孔梁板弹性屈曲问题分析的有效性和正确性。通过分析开孔板梁的无网格法数值分析特性，为分析各种形状的蜂窝状板梁结构提供了一种新的思路。
-
在目前的研究中，分析板屈曲的一般计算公式都是在切线刚度矩阵（
${{{K}}_{\rm{T}}}$ ）奇异性的基础上发展起来的，一般由问题域材料刚度矩阵（${{{K}}_{\rm{E}}}$ ）和几何刚度矩阵（${{{K}}_{\rm{G}}}$ ）组成。而KE只是单纯基于Kirchhoff板理论，KG的计算仅仅需要一阶运动学方程的弹簧旋转类比（rotational spring analogy，RSA）原理类推简化而来[20]。在预先给定的确定平面应力的线性分析中，这些都较容易获得。作为大多数弹性结构问题，KE可被认为是独立的内部力量。这种方法虽考虑了材料的型线响应，但在发生变形屈曲点之前，忽略了小变形，而这是一个很重要的假定，正是基于这个假定，从而把屈曲问题作为一个线性特征值问题，直接得出其简化计算公式。本文采用无网格局部Petrov-Galerkin法来离散问题域，然后使用最小二乘法来构建问题域位移场近似函数。MLPG法是一种纯粹的无网格法，由于它在真正意义上没有网格，与其他无网格法相比，它更适合于自适应计算。
在建立开孔梁的严格屈曲模式过程中，文章采用假定模式和互补模式结合的方法，从而组成一个降低特征值求解的二阶问题，然后通过迭代的方法[21]使原来多自由度问题能够得到有效解答。在迭代运算中，局部域的应用在求解的精度和计算要求方面都非常有效[24~26]。
-
不同于有限元法所采用的分段函数（位移函数的形函数）是预定义在所有单元子域上的，MLPG法位移场函数是通过最小二乘法MLS来获得的。这种方法的主要特征是基于某一问题域
${\bf{\Omega }}$ 内一系列给定节点的曲线拟合来获得一个连续逼近函数。通过采用多项式组合来完成，如$${{u}}^{\rm{h}}\left( x \right) = {p}^{\rm{T}}\left( x \right){a}\left( x \right)$$ (1) 式中：
${u}^{\rm{h}}\left( x \right)$ 为求解区域场函数近似函数；${p}^{\rm{T}}\left( x \right)$ 为多项式基函数；${a}\left( x \right)$ 为待定系数。对于$x y$ 平面内的问题，${{{X}}} = {[x\;y ]^{\rm{T}}}$ 为空间坐标系内的向量。${{\bf{p}}^T}\left( x \right) = [ 1\;x\;y\;{{x^2}}\;{xy}\;{{y^2}} ]$ 可以看做为一个二次基函数。参见文献[21]，MLS方法规定的节点参数值${u_I}$ 来构造近似函数，${\bf{a}}\left( x \right)$ 可以通过加权L2范数最小化来建立，可以表示为$$\frac{{\partial {\bf{\Delta }}}}{{\partial {\bf{a}}}} = 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\bf{\Delta }} = \sum\limits_{I = 1}^N {w\left( {{\bf{x}} - {{\bf{x}}_I}} \right){{\left[ {{{\bf{u}}^{\rm{h}}}\left( {{x_I}} \right) - {{\bf{u}}_I}} \right]}^2}} $$ (2) 其中
$w\left( {{\bf{x}} - {{\bf{x}}_I}} \right)$ 为权函数。把式（1）代入式（2），有
$${\bf{\Delta }} = \sum\limits_{I = 1}^N {w\left( {{\bf{x}} - {{\bf{x}}_I}} \right){{\left[ {{{\bf{u}}^{\rm{h}}}\left( {{x_I}} \right) - {{\bf{u}}_I}} \right]}^2}} $$ (3) 对式（3）泛函
${\bf{\Delta }}$ 中取极小值，可得$${\bf{A}}\left( x \right){\bf{a}}\left( x \right) = {\bf{B}}\left( x \right){{\bf{u}}_I}$$ (4) 其中
$$\left\{ { \begin{aligned} & {{\bf{A}}\left( x \right) = \sum\limits_{I = 1}^N {w\left( {{\bf{x}} - {{\bf{x}}_I}} \right){{\bf{p}}^{\rm{T}}}\left( {{{\bf{x}}_I}} \right){\bf{p}}\left( {{{\bf{x}}_I}} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } } \\ & {{\bf{B}}\left( x \right) =[ {\begin{array}{*{20}{c}} {{{\bf{b}}_1}}& \cdots &{{{\bf{b}}_N}} \end{array}} ],{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\bf{b}}_I} = w\left( {{\bf{x}} - {{\bf{x}}_I}} \right){\bf{p}}\left( {{{\bf{x}}_I}} \right)} \end{aligned}} \right.$$ (5) 定义
${\bf{\Phi }}\left( x \right)$ 为由节点参数得到的场逼近函数的形函数矩阵，设${\bf{\Phi }}\left( x \right)$ 为$${\bf{\Phi }}\left( x \right) = {{\bf{p}}^{\rm{T}}}\left( x \right)\left[ {{\bf{A}}{{\left( x \right)}^{ - 1}}{\bf{B}}\left( x \right)} \right]$$ (6) 这样可以把式（1）进一步简化表达为
$${{\bf{u}}^{\rm{h}}}\left( x \right) = {\bf{\Phi }}\left( x \right){{\bf{u}}_I}$$ (7) 对于表达公式中的权函数在最小二乘法中具有非常重要的作用，它可以用来控制影响域的位置，一般要求其具有紧支性。文章采用三次样条函数
$$w\left( r \right) = \left\{ { \begin{aligned} & {{2 / 3} - 4{r^2} + 4{r^3}}&{r \leqslant {1 / 2}}\qquad\qquad\qquad\;\;\\& {{4 / 3} - 4r + 4{r^2} - {{4{r^3}} / 3}}&{{1 / 2} < r \leqslant 1,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} r = \left| {x - {x_I}} \right|} \\& 0&{r > 1} \qquad\qquad\qquad\quad\; \end{aligned}} \right.$$ (8) -
为了求得承受平面载荷梁的外平面几何刚度矩阵
${{\bf{K}}_{\rm{G}}}$ ，需要得知与平面响应相关的应力分布。在使用MLPG法离散整个问题域时，简化离散过程是必要也是可能的，这样对于那些具有相同几何形状的单元来说，可以避免一些重复不必要的计算。文章提出方法是通过将开孔梁的整体响应进行分割成梁单元来进行简化，其中每个单元构成一个被降低了自由度的超级单元。图2为一个简化单元，共有四个质心超级节点，位置对应质心，每个节点有三个平面自由度$\left( {{u_i},{v_i},{\theta _i}} \right)$ 。在不规则孔洞板梁情况下，需要使用更多的单元模型来表示每个不同开口。因此，基于离散法得到的单元对于那些开有不规则孔洞的问题更为有益。
图 2 简化单元力学模型
Figure 2. simplified element mechanical model
单个单元（即不包括相同的单元）通过MLPG法离散建立起它的特征单元响应，然后再利用一个标准的离散装配过程[23]把所有单元组装成整体的梁响应。特征单元的响应是通过考虑单元力交变载荷的情况下获得的，在节点1处，当单元在外载荷作用下能够达到平衡，那么对于节点2，3，4，每一个单元对应于一个超级单元自由度，如图2所示。每个负载的情况都与MLPG近似一致的一组节点位移有关，按照MLPG方法逼近，这为转化4个节点超级单元平面刚度去替代单个单元响应提供了一个灵活的方法。
由于是基于单元的全梁分析，得到的应力仅在局部单元域是连续的，其产生的相关误差对于典型的梁结构来说是可以忽略不计的。此外，一个连续平面应力场可以在整体梁内，通过MLS来近似实现。下面给出单个单元平面应力场的计算公式
$${\bf{\sigma }} = \left\{ {\begin{array}{*{20}{c}} {{{\rm{\sigma}} _x}} \\ {{{\rm{\sigma}} _y}} \\ {{{\rm{\sigma}} _{\textit{z}}}} \end{array}} \right\} = {\bf{D}}\varepsilon = {\bf{D}}{{\bf{B}}_{xy}}\overbrace {\left[ {{{\bf{T}}_{xy}}\left\{ {\begin{array}{*{20}{c}} {\bf{F}} \\ \cdots \\ w \end{array}} \right\}} \right]}^{{u_I}}$$ (9) 式中
${\bf{D}}$ 为平面应力弹性矩阵，${{\bf{B}}_{xy}}$ 为应变位移矩阵，${{\bf{T}}_{xy}}$ 是在MLPG单元域内节点参数对应超级单元离散载荷和分布载荷的转换矩阵。式（9）封装了平面离散分析的两个层次，即MLPG和超级单元两个层次。 -
由文献[20]可知，平面外响应的切线刚度矩阵
${{\bf{K}}_T}$ 是由${{\bf{K}}_{\rm{E}}}$ 和${{\bf{K}}_{\rm{G}}}$ 这两个部分组成。 -
根据Kirchhoff薄板理论[22]，运用MLPG法，
${{\bf{K}}_{\rm{E}}}$ 可以分解为两个分量一个是弯曲刚度矩阵（${\bf{K}}$ ），另一个是施加本质边界条件得到的罚函数刚度矩阵（${{\bf{K}}_\alpha }$ ）$${{\bf{K}}_{\rm{E}}} = {\bf{K}} + {{\bf{K}}_\alpha }$$ (10) 其中
$${\bf{K}} = \int_\Omega {{\bf{B}}_{xy}^{\rm{T}}{\bf{D}}{{\bf{B}}_{xy}}{\rm{d}}\Omega } $$ (11) 其中
$${\bf{D}} = \frac{{{\bf{E}}t_w^3}}{{12\left( {1 - {{\rm{\mu}} ^2}} \right)}}\left[ {\begin{array}{*{20}{c}} 1&{\rm{\mu}} &0 \\ {\rm{\mu}} &1&0 \\ 0&0&{\dfrac{{\left( {1 - {{\rm{\mu}} ^2}} \right)}}{2}} \end{array}} \right]$$ (12) 式中
${\rm{E}}$ 为杨氏模量，$\mu $ 为泊松比，${t_w}$ 为梁腹板厚度。${{\bf{B}}_{xy}}$ 是应变位移矩阵，由下式确定，$${{\bf{B}}_{xy}} = \left[ {\begin{array}{*{20}{c}} {{\bf{\Phi }}{{\left( x \right)}_{,xx}}} \\ {{\bf{\Phi }}{{\left( x \right)}_{,yy}}} \\ {2{\bf{\Phi }}{{\left( x \right)}_{,{\textit{z}}{\textit{z}}}}} \end{array}} \right]$$ (13) 式（13）中由MLS得出的形函数矩阵
${\bf{\Phi }}\left( x \right)$ 下标号表示为对空间下标的导数。进一步，罚函数刚度矩阵（${{\bf{K}}_{\rm{\alpha}} }$ ）可以定义为，$${{\bf{K}}_\alpha } = \int_{{\Gamma _u}} {{\bf{R}}_{xy}^{\rm{T}}\left[ {\begin{array}{*{20}{c}} {{\alpha _ }}&0 \\ 0&{{\alpha _\theta }} \end{array}} \right]} {{\bf{R}}_{xy}}{\rm{d}}\Gamma $$ (14) 式中
${\Gamma _u}$ 是规定的位移边界，${{\bf{R}}_{xy}}$ 是连接节点位移和节点旋转自由度的一个转换矩阵。罚因子${\alpha _z}$ 和${\alpha _\theta }$ 为常数，通常其范围为刚度矩阵最大对角元素值的倍数，位数取值范围为${10^{4 - 13}}$ [23]。从施加约束出发，罚因子应是无穷的，但在实际计算中是不可能实现的。实际计算中如罚因子取值过小，可能会造成约束不能精确施加，如果取值过大，则计算可能会溢出或者得到病态解，为了简化计算，本文罚因子取值为弯曲刚度矩阵${\bf{K}}$ 最大对角线元素值的倍数，且使得${{\rm{\alpha}} _z} = {{\rm{\alpha}} _{\rm{\theta}} }$ 。 -
几何刚度
${{\bf{K}}_{\rm{G}}}$ 是根据RSA方法从平面应力场${{\rm{\sigma}} _x}$ 、${{\rm{\sigma}} _y}$ 和${{\rm{\tau}} _{xy}}$ 得到的[20]。使用相同的离散化方法对式（10）中的几何刚度矩阵${{\bf{K}}_{\rm{G}}}$ 进行求解，得到$${{\bf{K}}_{\rm{G}}} = \int_\Omega {{\bf{T}}_{xy}^{\rm{T}}\left[ {\begin{array}{*{20}{c}} {{{\rm{\sigma}} _x}}&{{{\rm{\tau}} _{xy}}} \\ {{{\rm{\tau}} _{yx}}}&{{{\rm{\sigma}} _y}} \end{array}} \right]} {{\bf{T}}_{xy}}{t_w}{\rm{d}}\Omega $$ (15) 式中，
${{\bf{T}}_{xy}}$ 是节点自由度做旋转等效变换的一个几何变换矩阵，可以通过形函数的导数来确定，如下式$${{\bf{T}}_{xy}} = \left[ {\begin{array}{*{20}{c}} {{\bf{\Phi }}{{\left( x \right)}_{,x}}} \\ {{\bf{\Phi }}{{\left( x \right)}_{,y}}} \end{array}} \right]$$ (16) 式中，
${\bf{\Phi }}{\left( x \right)_{,x}}$ 和${\bf{\Phi }}{\left( x \right)_{,y}}$ 为由MLS法求得的形函数的一阶偏导数。 -
折边对开孔梁的局部屈曲响应的影响可以通过一种简化模型来实现。对于顶部和底部都有折边的组合工字梁来说，如果这些折边能够沿着连接腹板内边缘起到完全约束作用的话，那么通常可以只考虑它们的材料和几何刚度。如果是有上下板约束的话，则可以把板结构进行等效处理，从而简化为梁结构。相对于三维问题的全方位建模（因为折边在横向平面，即
$x - {\textit{z}}$ 平面）来说，上面提到的简化方法仅需要采用一维模型即可达到简化折边响应。同时这种方法保留了二维平面（即$x - y$ 平面）模型的适用性。考虑一个一维折边
${\Gamma _f}$ ，其尺寸为${t_f} \times {b_f}$ ，其中${t_f}$ 表示为折边厚度，${b_f}$ 表示为折边宽度，在忽略次要变形的影响下，那么折边材料刚度贡献来自其扭转刚度${\bf{GJ}}$ ，则折边的材料刚度矩阵为$${{\bf{K}}_{{\rm{E}}f}} = \int_{{\Gamma _f}} {{\bf{\Phi }}_{,xy}^{\rm{T}}{\bf{GJ}}{{\bf{\Phi }}_{,xy}}} {\rm{d}}\Gamma $$ (17) 式中
${{\bf{\Phi }}_{,xy}}$ 为关于梁的横向位移形函数的二阶偏导数，对于薄板的扭转刚度由下式确定$${\bf{GJ}} = \frac{{\bf{E}}}{{2\left( {1 + {\rm{\mu}} } \right)}}\left[ {\frac{{{b_f}t_f^3}}{3}} \right]$$ (18) 另一方面，为了建立相应的几何刚度矩阵，假定沿着折边宽度和厚度方向的应力
${{\rm{\sigma}} _x}$ 为一个常数，在沿着$x$ 轴方向的内部压缩或者拉伸应力变量${{\rm{\sigma}} _x}$ 作用下，可以认为该梁为一个一维等效梁。因此，分布在横截面上每单位面积的等效转动几何刚度可以取为${{\bf{K}}_{\rm{\theta}} } = {{\bf{\sigma }}_x}$ ，等效弯曲的侧向旋转角由下式确定[21]$${\bf{\theta }} = \sqrt {{y^2} + {{\textit{z}}^2}} \left( {{{\bf{\Phi }}_{,xy}}{{\bf{u}}_{\textit{z}}}} \right) \Rightarrow {{\bf{T}}_{\rm{\theta}} } = \sqrt {{y^2} + {{\textit{z}}^2}} {{\bf{\Phi }}_{,xy}}$$ (19) 式中
${{\bf{u}}_z}$ 代表腹板平面平移自由度。因此，根据RSA法，得到折边的几何刚度$$ \begin{split} & \qquad{{\bf{K}}_{{\rm{G}}f}} = \int_{{\Gamma _f}} {\int_\Omega {{\bf{T}}_{\rm{\theta}} ^{\rm{T}}{\rm{diag}}\left( {{{\bf{K}}_{\rm{\theta}} }} \right){{\bf{T}}_{\rm{\theta}} }{\rm{d}}\Omega {\rm{d}}\Gamma } } = \\& \int_{{\Gamma _f}} {{\bf{\Phi }}_{,xy}^{\rm{T}}{\rm{diag}}\left( {{{\bf{\sigma }}_x}} \right){{\bf{\Phi }}_{,xy}}\left( {\int_\Omega {\left( {{y^2} + {{\textit{z}}^2}} \right){\rm{d}}\Omega } } \right){\rm{d}}\Gamma } = \\& \qquad\qquad{I_o}\int_{{\Gamma _f}} {{\bf{\Phi }}_{,xy}^T\left( {{{\bf{\sigma }}_x}} \right){{\bf{\Phi }}_{,xy}}d\Gamma } \end{split} $$ (20) 式中，
${I_O}$ 表示面积二次极矩。 -
目前对开孔梁的失稳评估模型一般是在仅取梁的某一适当部分的基础上建立起来的。文章采用局部模型，沿梁的长度方向取至少四个单元，这样具有更好的分析效果。文章提出的方法只需要在无网格MLPG区域内计算材料和几何刚度矩阵，大大减少了刚度矩阵的规模，从而提高了计算效率。在梁单元局部域内，按照比例施加载荷情况下的梁的临界载荷因子（
${{\rm{\lambda}} _{{\rm{cr}}}}$ ）可以由相关的${{\bf{K}}_{\rm{E}}}$ 和${{\bf{K}}_{\rm{G}}}$ 的特征值求得，然而由于相关矩阵的大小直接导致求解其特征值的工作量庞大，从而大大增加了计算成本。考虑到RSA方法中[20]，当梁在屈曲模态时，${{\rm{\lambda}} _{{\rm{cr}}}}$ 也可以作为内部能量吸收和转动等效转化的比例系数。文章提出的方法提供了一个实际屈曲载荷因子的近似模态代替了精确模态的束缚，典型的模态可以近似地通过对线性结构执行屈曲模态来获得[20, 21]。
图 3 问题域离散节点图
Figure 3. Problem domain discrete node graph
考虑一个近似平面外模态（
${{\bf{u}}_z}$ ）由对应于外部的平面力（${{\bf{f}}_z}$ ）获得和，忽略屈曲前的平面内变形的影响，则临界屈曲载荷因子${\lambda _{{\rm{cr}}}}$ 可以由下式确定[22 ]，$${{\rm{\lambda}} _{{\rm{cr}}}} = - \frac{{{\bf{u}}_{\textit{z}}^{\rm{T}}{{\bf{f}}_{\textit{z}}}}}{{{\bf{u}}_{\textit{z}}^{\rm{T}}{{\bf{K}}_{\rm{G}}}{{\bf{u}}_{\textit{z}}}}} = - \frac{{{\bf{u}}_{\textit{z}}^{\rm{T}}{{\bf{K}}_{\rm{E}}}{{\bf{u}}_{\textit{z}}}}}{{{\bf{u}}_{\textit{z}}^{\rm{T}}{{\bf{K}}_{\rm{G}}}{{\bf{u}}_{\textit{z}}}}}$$ (21) 基于假定的屈曲模态，对预定的模式通过进行迭代运算来进一步改进，式（21）可以用来获得一个初始解的一个近似临界屈曲载荷因子
${{\rm{\lambda}} _{{\rm{cr}}}}$ 。文章方法中，通过局部域的手段来降低问题的维数，从而使得假定的模态的阶数降低，这样特征值的求解变得容易[24~26]。根据这种方法，相应于一个任意的平面外载荷模式（${{\bf{f}}_{{\textit{z}}{\rm{A}}}}$ ），最初的近似模态（${{\bf{u}}_{{\textit{z}}{\rm{A}}}}$ ）在一个新的负载模式（${{\bf{f}}_{{\textit{z}}{\rm{B}}}}$ ）作用下，可以建立一个互补的模态（${{\bf{u}}_{{\textit{z}}{\rm{B}}}}$ ）。${{\bf{f}}_{{\textit{z}}{\rm{B}}}}$ 定义如下，$${{\bf{f}}_{{\textit{z}}{\rm{B}}}} = - \left[ {{{\bf{f}}_{{\textit{z}}{\rm{A}}}} + {\lambda _{cr\left( {\rm{A}} \right)}}\left( {{{\bf{K}}_{\rm{G}}}{{\bf{u}}_{{\textit{z}}{\rm{A}}}}} \right)} \right] \Rightarrow {{\bf{u}}_{{\textit{z}}{\rm{B}}}} = {\bf{K}}_{\rm{G}}^{ - 1}{{\bf{f}}_{{\textit{z}}{\rm{B}}}}$$ (22) 因此，两个互补模态
${{\bf{u}}_{{\textit{z}}{\rm{A}}}}$ 和${{\bf{u}}_{{\textit{z}}{\rm{B}}}}$ 用于计算和解答经降阶了的模态特征值，结果经过每一次迭代计算直至收敛。 -
算例采用一连续多孔板，在板沿长度方向等距开有半径
$r = 250\;{\rm{mm}}$ 的系列圆形孔，对称半长$L = 40\;000\;{\rm{mm}}$ ，宽度$B = 2\;000\;{\rm{mm}}$ ，板厚$t = 5\;{\rm{mm}}$ 。这种结构在船舶结构中非常普遍，例如肋板、船底纵桁等，这些结构的一般特征就是开有一定数量的减轻孔或者人孔。图4为该连续多孔板结构的靠近两端的带4.5个圆形孔的部分简化形式。算例中计算用到的计算常量有，材料杨氏弹性模量${\rm{E}} = 3.0 \times {10^5}\;{\rm{MPa}}$ ，泊松比$\mu = 0.3$ ，密度$\rho = 7.5 \times {10^3}\;{\rm{Kg/}}{{\rm{m}}^{\rm{3}}}$ 。板的的上下端为自由边，左侧边界刚性固定，右侧平直边部分受到与板中性面方向一致的均布压力$q = 2 \times {10^4}\;{\rm{N}}$ 。
图 4 四个孔单元结构图（单位：mm）
Figure 4. Structural drawing of four hole units (unit: mm)
算例运用MLPG无网格算法，初始节点42108个（如图5a)所示），采用高斯积分法，积分背景初始网格为83194个，经过五步自适应计算，得到图5b)计算结果。
图 5 无网格MLPG法计算结果（单位：Pa）
Figure 5. calculation results of meshless MLPG method (Unit: Pa)
为了进行比较分析，算例运用有限元法进行验证性计算（使用NASTRAN 2016），图6a)为有限元离散网格，共4900个节点，4583个网格，计算得到图6b)计算结果。
图 6 有限元法（MSC. NASTRAN 2016）计算结果（单位：Pa）
Figure 6. calculation results of finite element method (MSC. NASTRAN 2016) (Unit: Pa)
从两种方法的计算结果来看，计算的应力分布高度相似，并且应力集中的区域几乎都落在相同的区域。从无网格方法的计算结果来看，由于其采用的是连续的应力场函数，不仅可以清晰的得到每一点的应力分布，甚至在那些边缘应力出现间断都可以获得。这就是无网格法相对于有限元的重要优势。
为了进一步分析文章所提方法对结构局部屈曲载荷因子分析的有效性，下面采用传统方法（有限元法）与之分析结果进行比较，两种方法计算结果如图7所示。从图7所示图线，可以得出结论，有限元法基本上对于整个结构的屈曲载荷因子不敏感，对于结构当中的孔洞所带来的的局部结构屈曲并没有明显的响应。而运用文章所提出的方法不难得出结论，在结构的两个端部屈曲载荷因子较大，结构的中部屈曲载荷因子降低，该结果与原结构的值吻合得较好。
图 7 运用MLPG法和有限元法计算出屈曲载荷因子的变化图线
Figure 7. Calculate the change diagram of buckling load factor with MLPG method and FEM
为了更好地分析文章所提出的方法能够更好地分析临界屈曲载荷因子的大小随结构变化的敏感性，即分析结构屈曲载荷因子随结构变化而变化情况，将原结构中两侧第二和第三个圆形孔按照以半径作为边长修改为正方形，如图8所示。
图 8 修改后的结构示意图（单位：mm）
Figure 8. Modified structure diagram (unit: mm)
运用MLPG无网格算法，初始节点44 400个（如图9a)所示），采用高斯积分法，积分背景初始网格为87 923个，经过五步自适应计算，得到图9b)计算结果。
图 9 无网格MLPG法计算结果（单位：Pa）
Figure 9. calculation results of meshless MLPG method (Unit: Pa)
同样采用有限元法进行比较分析（使用NASTRAN 2016），图10a)为有限元离散网格，共3585个节点，3127个网格，计算得到图10b)计算结果。
图 10 有限元法（MSC. NASTRAN 2016）计算结果（单位：Pa）
Figure 10. Calculation results of finite element method (MSC. NASTRAN 2016) (Unit: Pa)
将经过修改的结构，再次运用两种方法来计算，从计算应力云图结果来看，应力分布依然高度相似，并且应力集中趋势也基本一致。再来看看重新计算出的屈曲载荷因子，如图11所示，有限元法计算对结构修改后的载荷因子与原结构载荷因子几乎相同，说明其对于孔洞的细微修改并不不敏感，而文章所提出的方法则非常敏感，两侧第二和第三孔洞的扩大（由半径
$r = $ 250 mm扩大为边长$a = $ 250 mm的正方形）明显降低了其所在位置的屈曲载荷因子，同时也对其两侧的结构也产生了一定的影响，而再于其外的结构基本不产生影响，这个计算结果也与原结构有着更高的契合度。
图 11 运用MLPG法和有限元法计算出屈曲载荷因子的变化图线
Figure 11. Change diagram of buckling load factor calculated by MLPG method and finite element method
-
文章采用无网格局部Petrov-Galerkin法计算和评估开孔梁和开孔板弹性屈曲问题。在建立开孔梁的严格屈曲模式过程中，采用假定模式和互补模式结合的方法，进行降阶求解模态特征值。在迭代运算中，局部域的应用在求解的精度和计算要求方面都是非常有效的。
算例运用有限元法和文章提出方法进行计算，通过对应力计算结果分析和比较，得到两种方法的计算应力分布高度相似，应力集中区域基本相同。文章方法由于其采用的是连续应力场函数，可以清晰的得到场点的应力分布。
通过算例结构局部屈曲载荷因子的计算和比较，得出文章所提出的方法能够较好地仿真结构两个端部屈曲载荷因子，该结果与原结构有着较好的吻合度。
算例还测试了文章方法和有限元法临界屈曲载荷因子的大小随结构变化的敏感性，传统方法（文章采用有限元法）对于板梁孔洞的细微修改不敏感，而文章所提出的方法则非常敏感，这个计算结果与原结构有着更好的契合度。
MLPG Analysis Method for Local Buckling of Ship Opening Beams
-
摘要:
目的 文章针对船舶开孔梁板易产生局部屈曲效应、受到较高的集中载荷作用易产生明显结构响应和屈曲破坏力的问题，提出采用无网格数值方法对该类问题进行数值分析。 方法 该方法以问题域离散节点为基础，建立其紧支函数，运用加权余量法建立系统位移场离散方程，然后由移动最小二乘法来构造离散节点的形函数，从而获得位移函数的逼近函数。最后结合经简化的屈曲评估法方法，有效地将屈曲方程特征值求解进行降阶处理，最后得到更为精确的屈曲载荷因子。 结果 文章对所提出方法和传统方法计算连续开孔梁的屈曲载荷因子进行比较，验证研究方法的有效性和优点。通过算例的演示和比较，验证了文章方法的有效性和优点。 结论 传统方法对于整个结构的屈曲载荷因子不敏感，对于结构当中的孔洞所带来的的局部结构屈曲并没有明显的响应，文章所提方法在结构的两个端部屈曲载荷因子可以得到放大，计算结果与原结构结果吻合得较好；对于各个孔洞的屈曲载荷因子，传统方法对于孔洞的细微修改不敏感，而文章所提方法明显降低了其所在位置的屈曲载荷因子，同时也对其两侧的结构也产生了一定的影响，这个计算结果与原结构也有着更高的契合度。 Abstract:Objective in order to solve the problems of local buckling effect and obvious structural response and buckling failure force caused by high concentrated load, a meshless numerical method is proposed to simulate the problem. Methods Based on the discrete nodes in the problem domain, this method establishes its compact support function, uses the weighted residual method to establish the discrete equation of the system displacement field, and then uses the moving least square method to construct the shape function of the discrete nodes to finally obtain the approximation function of the displacement function. Finally, combined with the simplified buckling evaluation method, the eigenvalues of the buckling equation are effectively solved to reduce the order, and finally a more accurate buckling load factor is obtained. results in this paper, the proposed method and the traditional finite element method are used to calculate the buckling load factors of continuous beams with openings, and the results are compared to verify the effectiveness and advantages of the research method. Through the demonstration and comparison of examples, the effectiveness and advantages of the research method are verified. Conclusion ① the classical method is not sensitive to the buckling load factor of the whole structure, and has no obvious response to the local structure buckling caused by the holes in the structure. By using the method proposed in this paper, the two end buckling load factors of the structure can be amplified, and the calculated results are in good agreement with the experimental results. ② for the buckling load factors of each hole, the traditional method is not sensitive to the minor modification of the hole, while the method proposed in this paper significantly reduces the buckling load factor of its location, but also has a certain impact on the structure on both sides of it. This calculation result has a higher agreement with the experimental results. -
Key words:
- Ship structure /
- slab girder with holes /
- buckling analysis /
- load factor /
- meshless analysis
-
图 5 无网格MLPG法计算结果（单位：Pa）
Figure 5. calculation results of meshless MLPG method (Unit: Pa)
图 6 有限元法（MSC. NASTRAN 2016）计算结果（单位：Pa）
Figure 6. calculation results of finite element method (MSC. NASTRAN 2016) (Unit: Pa)
图 7 运用MLPG法和有限元法计算出屈曲载荷因子的变化图线
Figure 7. Calculate the change diagram of buckling load factor with MLPG method and FEM
图 9 无网格MLPG法计算结果（单位：Pa）
Figure 9. calculation results of meshless MLPG method (Unit: Pa)
图 10 有限元法（MSC. NASTRAN 2016）计算结果（单位：Pa）
Figure 10. Calculation results of finite element method (MSC. NASTRAN 2016) (Unit: Pa)
-
KERDAL D, NETHERCOT D A. Failure modes for castellated beams[J]. Journal of Constructional Steel Research, 1984, 4(4): 295–315. doi: 10.1016/0143-974X(84)90004-X LAWSON R W, LIM J, HICKS S J, et al. Design of composite asymmetric cellular beams and beams with large web openings[J]. Journal of Constructional Steel Research, 2006, 62(6): 614–629. doi: 10.1016/j.jcsr.2005.09.012 张文福, 邓云, 赵文艳, 等. 基于板梁理论的工字形钢—混组合双跨连续梁弯扭屈曲分析[J]. 东北石油大学学报, 2017, 41(4): 107–115. doi: 10.3969/j.issn.2095-4107.2017.04.012 ZHANG W F, DENG Y, ZHAO W Y, et al. Flexural-torsional buckling analysis on I-shaped steel-concrete composite double-span continuous beams based on plate-beam theory[J]. Journal of Northeast Petroleum University, 2017, 41(4): 107–115 (in Chinese). doi: 10.3969/j.issn.2095-4107.2017.04.012 LIU T C H, CHUNG K F. Steel beams with large web openings of various shapes and sizes: finite element investigation[J]. Journal of Constructional Steel Research, 2003, 59(9): 1159–1176. doi: 10.1016/S0143-974X(03)00030-0 陈彦廷, 于昌利, 桂洪斌. 船体板和加筋板的屈曲及极限强度研究综述[J]. 中国舰船研究, 2017, 12(1): 54–62. doi: 10.3969/j.issn.1673-3185.2017.01.009 CHEN Y T, YU C L, GUI H B. Research development of buckling and ultimate strength of hull plate and stiffened panel[J]. Chinese Journal of Ship Research, 2017, 12(1): 54–62 (in Chinese). doi: 10.3969/j.issn.1673-3185.2017.01.009 British Standards Institution. Structural use of steelwork in building——Part 1: Code of practice for design—rolled and welded sections: BS 5950-1: 2000[S]. London, UK: British Standards Institution, 2001. British Standards Institution. Eurocode 3: Design of steel structures——Part 1-1: General rules and rules for buildings: BS EN 1993-1-1: 2005[S]. London, UK: British Standards Institution, 2005. 殷毓基, 梁珂, 孙秦. 基于非线性有限元降阶方法的网格加筋筒壳结构屈曲承载特性研究[J]. 宇航总体技术, 2019, 3(1): 17–22. YIN Y J, LIANG K, SUN Q. Buckling load-carrying capability analysis of gridstiffened cylinders using a nonlinear finite element reduced-order method[J]. Astronautical Systems Engineering Technology, 2019, 3(1): 17–22 (in Chinese). 张旭东, 邱志平, 李琦. 模糊变分原理在求解结构屈曲临界载荷中的应用[J]. 工程力学, 2015, 32(8): 29–35. ZHANG X D, QIU Z P, LI Q. Application of the fuzzy variational principle in solving the critical buckling load of structures[J]. Engineering Mechanics, 2015, 32(8): 29–35 (in Chinese). 沃周华, 周上然. 船体板梁结构屈曲强度分析[J]. 中国设备工程, 2019(5): 165–166. WO Z H, ZHOU S R. Buckling strength analysis of hull plate girder structure[J]. China Plant Engineering, 2019(5): 165–166 (in Chinese). 程昊, 秦朝红, 孔凡金, 等. 壁板结构热屈曲后模态特性试验[J]. 振动、测试与诊断, 2019, 39(2): 306–310, 443. CHENG H, QIN Z H, KONG F J, et al. Experimental research on modal characteristics of a thermal post-buckling panel[J]. Journal of Vibration, Measurement & Diagnosis, 2019, 39(2): 306–310, 443 (in Chinese). ODEN J T, DUARTE C A M, ZIENKIEWICZ O C. A new cloud-based hp finite element method[J]. Computer Methods in Applied Mechanics and Engineering, 1998, 153(1–2): 117–126. doi: 10.1016/S0045-7825(97)00039-X LIU G R, GU Y T. Meshless local Petrov–Galerkin (MLPG) method in combination with finite element and boundary element approaches[J]. Computational Mechanics, 2000, 26(6): 536–546. doi: 10.1007/s004660000203 何沛祥, 李子然, 吴长春. 无网格与有限元的耦合在动态断裂研究中的应用[J]. 应用力学学报, 2006, 23(2): 195–198. doi: 10.3969/j.issn.1000-4939.2006.02.006 HE P X, LI Z R, WU C C. Coupled finite element-element-free Galerkin method for dynamic fracture[J]. Chinese Journal of Applied Mechanics, 2006, 23(2): 195–198 (in Chinese). doi: 10.3969/j.issn.1000-4939.2006.02.006 段念, 王文珊, 于怡青, 等. 基于FEM与SPH耦合算法的单颗磨粒切削玻璃的动态过程仿真[J]. 中国机械工程, 2013, 24(20): 2716–2721. doi: 10.3969/j.issn.1004-132X.2013.20.004 DUAN N, WANG W S, YU Y Q, et al. Dynamic simulation of single grain cutting of glass by coupling FEM and SPH[J]. China Mechanical Engineering, 2013, 24(20): 2716–2721 (in Chinese). doi: 10.3969/j.issn.1004-132X.2013.20.004 JOHNSON G R, STRYK R A, BEISSEL S R, et al. An algorithm to automatically convert distorted finite elements into meshless particles during dynamic deformation[J]. International Journal of Impact Engineering, 2002, 27(10): 997–1013. doi: 10.1016/S0734-743X(02)00030-1 高欣. 一致性高阶无单元伽辽金法及裂纹扩展分析[D]. 大连: 大连理工大学, 2018. GAO X. Consistent high order element-free Galerkin method and crack propagation analysis[D]. Dalian: Dalian University of Technology, 2018 (in Chinese). LIU G R. Mesh Free Methods: Moving Beyond the Finite Element Method[M]. Boca Raton, FL, USA: CRC Press, 2003. 王毅刚. 子域无网格伽辽金法及其在固体力学中的应用[D]. 长沙: 湖南大学, 2018. WANG Y G. A sub-domain element free Galerkin method for solid mechanics problems[D]. Changsha: Hunan University, 2018 (in Chinese). IZZUDDIN B A. Rotational spring analogy for buckling analysis[J]. Journal of Structural Engineering, 2007, 133(5): 739–751. doi: 10.1061/(ASCE)0733-9445(2007)133:5(739) IZZUDDIN B A. Simplified buckling analysis of skeletal structures[J]. Proceedings of the Institution of Civil Engineers-Structures and Buildings, 2006, 159(4): 217–228. doi: 10.1680/stbu.2006.159.4.217 郭瑶. 基于无网格伽辽金法的结构动力分析[D]. 西安: 西安理工大学, 2017. GUO Y. Dynamic analysis of a structure based on meshless method[D]. Xi’an: Xi’an University of Technology, 2017 (in Chinese). LANCASTER F D. Meshless methods and their application to cellular beams[R]. South Kensington, London: Imperial College London, 2010. 李冬睿, 许统德. 自适应邻域选择的数据可分性降维方法[J]. 计算机应用, 2012, 32(8): 2253–2257. LI D R, XU T D. Dimensionality reduction method with data separability based on adaptive neighborhood selection[J]. Journal of Computer Applications, 2012, 32(8): 2253–2257 (in Chinese). 张延良, 卢冰, 洪晓鹏, 等. 微表情识别的局部区域方法[J]. 计算机应用, 2019. doi: 10.11772/j.issn.1001-9081.2018102090 ZHANG Y L, LU B, HONG X P, et al. Micro-expression recognition based on local region method[J]. Journal of Computer Applications, 2019 (in Chinese). doi: 10.11772/j.issn.1001-9081.2018102090 蒲玲. 自适应局部线性降维方法[J]. 计算机应用与软件, 2013, 30(4): 255–257. doi: 10.3969/j.issn.1000-386x.2013.04.073 PU L. Adaptive local linear dimensional reduction method[J]. Computer Applications and Software, 2013, 30(4): 255–257 (in Chinese). doi: 10.3969/j.issn.1000-386x.2013.04.073 -
点击查看大图
图(11)
计量
- 文章访问数: 32
- HTML全文浏览量: 4
- PDF下载量: 1142
- 被引次数: 0
