2. School of Ship Engineering, Guangzhou Maritime Institute, Guangzhou 510725, China
当前船舶结构分析最常用也是最为有效的计算工具之一就是有限元方法[1]。但有限元法因其自身的特点,在处理分析船舶结构场量(位移场和应力场等)变化剧烈的高梯度区域时,会出现计算精度降低甚至计算中断现象,为了解决这个问题,有限元法通常采用的方法就是对该区域的网格进行加密(细分)或者采用高阶单元,这样就要求方法具有较强的自适应分析能力,其结果就加大了有限元法前后处理的工作量,从而降低了计算效率,事实上这种方法也并不能从根本上消除问题产生的根源。作为与有限元法相对应的另一种数值分析方法——无网格法,在近20年中得到了很大的发展[2, 3, 4, 5, 6]。无网格法是建立在系列独立的离散点的基础上,通过构造点的近似函数来求解问题。与传统的有限元法相比,无网格法无需网格背景,在计算过程中可以根据需要任意增减节点,而不需要处理节点之间的拓扑信息,特别适合用来进行自适应分析计算。无网格法目前在航空材料、高速碰撞、动态裂纹扩展、加工成型、节理岩体分析等诸多领域都得到了较为广泛的应用[7, 8, 9, 10]。无网格法的计算精度同有限元法一样因其自身的特点,受节点密度、基函数的选用及其阶次以及支撑域半径大小等因素的影响[11, 12]。文章用局部无网格Petrov-Galerkin法为基础,利用m阶B样条函数作为小波基函数来构造对船舶结构应力场的逼近函数,并采用两尺度分解技术来分解应力场的高梯度成分和低尺度成分。文章最后选取了两种典型船舶结构进行变形和应力分析,并与有限元法(ANSYS软件)的计算结果进行了比较来验证文章方法的有效性。
1 船舶结构MLPG法基本表达由于船体板的板厚远小于其另外两个方向(长度和宽度)的几何尺寸,在分析其受力时,可以运用Kirchhoff-Love板壳理论和Mindlin-Reissner板壳理论。文章基于板壳的Mindlin-Reissner基本理论,运用局部无网格Petrov-Galerkin法(MLPG)来分析研究无网格法下的船舶板结构应力。
1.1 节点控制方程
板结构经离散后二维弹性理论为基础的节点系统方程[11]:
σij,i+bi=0
(1)
定义Ω为问题域,Γ为Ω的边界;问题的本质边界条件为:ui=${{{\bar{u}}}_{i}}$,Γu;问题的自然边界条件为:σijnj=${{{\bar{t}}}_{i}}$,Γt,nj为自然边界上单位外法向矢量的第j个分量。
利用加权余量法[12],可以得到节点I处系统方程的微分方程强形式为
∫ΩQσij,i+bi${{{\tilde{W}}}_{I}}$dΩ-α∫ΓQu(ui-${{{\bar{u}}}_{i}}$)${{{\tilde{W}}}_{I}}$dΓ=0
(2)
对节点I依式(2)在其积分域上进行积分,可以建立起它的系统方程。这样对每一个节点都采用式(2)在其积分域上积分,可以得到所有离散节点的系统方程,将这些离散节点的系统方程组装起来就能够获得问题域的整体系统方程。
利用移动最小二乘法(MLS)得到节点积分点支撑域内的形函数,从而获得问题域的位移逼近函数:
uhX=ΦTX·u=$\sum\limits_{I=1}^{N}{{}}$φI(
(3)
根据弹性力学应力-应变关系有
σ=Dε=$D\left[ \begin{array}{*{35}{l}}
\frac{\partial }{\partial x} & 0 \\
0 & \frac{\partial }{\partial y} \\
\frac{\partial }{\partial y} & \frac{\partial }{\partial x} \\
\end{array} \right]\sum\limits_{j}^{n}{{{\Phi }_{j}}{{u}_{j}}=D\sum\limits_{j}^{n}{{{B}_{j}}{{u}_{j}}}}$
(4)
${{B}_{j}}=\left[ \begin{array}{*{35}{l}}
{{\phi }_{j,x}} & 0 \\
0 & {{\phi }_{j,y}} \\
{{\phi }_{j,y}} & {{\phi }_{i,x}} \\
\end{array} \right]$
设${{{\tilde{V}}}_{I}}$为权函数
${{{\tilde{W}}}_{I}}$的导数矩阵,这样式(2)可以进一步表达为
∫ΩQ
$\tilde{V}_{I}^{T}$σdΩ+α∫ΓQu${{{\tilde{W}}}_{I}}$udΓ-∫ΓQu${{{\tilde{W}}}_{I}}$NσdΓ=
∫ΓQt${{{\tilde{W}}}_{I}}$${\bar{t}}$dΓ+α∫ΓQu${{{\tilde{W}}}_{I}}$${\bar{u}}$dΓ+∫ΩQ${{{\tilde{W}}}_{I}}$bdΩ
(5)
把节点系统方程式(3)和式(4)代入节点控制方程式(5),可以得到
$\begin{align}
& \int{_{{{\Omega }_{Q}}}\tilde{V}_{I}^{T}\sum\limits_{j}^{n}{{{B}_{j}}{{u}_{j}}d\Omega +\alpha }\int{_{{{\Gamma }_{{{Q}_{u}}}}}}}{{{\tilde{W}}}_{I}}\sum\limits_{j}^{n}{{{\Phi }_{j}}{{u}_{j}}} \\
& \int{_{{{\Gamma }_{{{Q}_{u}}}}}}{{{\tilde{W}}}_{I}}N\sum\limits_{j}^{n}{{{B}_{j}}{{u}_{j}}d\Gamma }= \\
& \int{_{{{\Gamma }_{{{Q}_{u}}}}}}{{{\tilde{W}}}_{I}}\bar{t}d\Gamma \alpha \int{_{{{\Gamma }_{{{Q}_{u}}}}}}{{{\tilde{W}}}_{I}}\bar{u}d\Gamma +\int{_{{{\Omega }_{Q}}}}{{{\tilde{W}}}_{I}}bd\Omega \\
\end{align}$
(6)
KIu=FI
(7)
$\begin{align}
& {{K}_{I}}=\sum\limits_{j}^{n}{{{K}_{lj}}=}\int{_{{{\Omega }_{Q}}}}\tilde{V}_{I}^{T}\sum\limits_{j}^{n}{{{B}_{j}}d\Omega }+\sum\limits_{j}^{n}{{}} \\
& \alpha \int{_{{{\Gamma }_{{{Q}_{u}}}}}}{{{\tilde{W}}}_{I}}\sum\limits_{j}^{n}{{{\Phi }_{j}}d\Omega -}\int{_{{{\Gamma }_{{{Q}_{u}}}}}}{{{\tilde{W}}}_{I}}N\sum\limits_{j}^{n}{{{B}_{j}}d\Gamma } \\
\end{align}$
这样由每个离散节点的控制方程总装成整个结构系统的控制方程表达为
Ku=F
(8)
由于在MLS中,采用的节点支撑域都是紧支的,系统的总体刚度矩阵K是带状稀疏矩阵,可以减少计算工作量。根据文献[13]对有限元法中罚因子的分析,建议罚因子α 105~12×E在范围内取值,E为杨氏弹性模量。
2 B样条小波基在MLPG法中,运用最小二乘法(MLS)是求解位移场逼近函数的关键,从节点控制方程的推导过程来看,权函数的选定及其参数的确定对于计算的精度和结果的收敛至关重要。在无网格法中,一般要求权函数为紧支函数,在求解域内具有连续可导、单调递减和正态紧支等特性[14]。由于小波函数具有在不同分辨率水平上表达函数的变尺度能力,而其相应的基函数(小波基)则在具有震荡性、衰减性的同时还具有紧支性或近似紧支性[15, 16];另一方面B样条函数具有分段光滑、局部可微和线性组合等特征[17],所以B样条小波函数的时频局部特性最接近无网格法中常用作权函数的Gauss函数,并且它的紧支性优于Gauss函数,这样保证了在计算处理过程中的相位失真。所以本文将利用小波函数和B样条函数的各自优点,来构造小波基函数,并以此构造出来的B样条小波基函数作为MLS法中的权函数。
对于m阶B样条函数空间是由节点序列{xkj}k=-m+12j+2m-1(节点总数为2j+2m-1)形成的样条函数构成,节点序列{xkj}k=-m+12j+2m-1由式(9)确定[14]:
$\begin{align} & \left\{ \begin{matrix} x_{-m+1}^{j}=x_{-m+1}^{j}=\cdots =0 \\ x_{k}^{j}=k{{2}^{j}} \\ x_{{{2}^{j}}+1}^{j}=x_{{{2}^{j}}+2}^{j}=\cdots x_{{{2}^{j}}+m-1}^{j}=1 \\ \end{matrix} \right. \\ & k=1,{{2}^{1}},\cdots ,{{2}^{j}} \\ \end{align}$ | (9) |
依据式(9)的节点序列可以进一步构造在0,1区间内的m阶局部支撑嵌套的j尺度逼近空间V0,1J,其基函数则可以表达为
$\left\{ \begin{matrix}
B_{m,k}^{j}\left( x \right)={{N}_{m}}\left( {{2}^{j}}\xi -k \right) \\
\sup pB_{m,k}^{j}=\left[ \xi _{k}^{j},\xi _{k+m}^{j} \right] \\
\end{matrix} \right.$
k=-m+1,-m+2,…,2j-1
(10)
在有限区间[0,1]上,若要小波尺度函数
φjm,kξ=Bjm,kξ,需满足条件2j≥2m-1。则对于任意尺度j的m阶B样条尺度函数φm,kj(ξ)可以表示为
φm,kjξ=$\left\{ \begin{matrix}
\phi _{m,k}^{l}\left( {{2}^{j-l}}\xi \right),k=-m+1,-m+2,\cdots ,-1 \\
\phi _{m,2j-m-k}^{l}\left( 1-{{2}^{j-l}}\xi \right),k={{2}^{j}}-m+1,\cdots ,{{2}^{j}}-1 \\
\phi _{m,0}^{l}\left( {{2}^{j-l}}\xi -{{2}^{-l}}k \right),k=0,1,\cdots ,{{2}^{j}}-m \\
\end{matrix} \right.$
(11)
φm,kj(ξ)对应的小波函数φm,kj
(ξ)可以表示为[18]
$\varphi _{m,k}^{j}\left( \xi \right)\left\{ \begin{matrix}
\varphi _{m,k}^{l}\left( {{2}^{j-l}}\xi \right),k=-m+1,-m+2,\cdots ,-1 \\
\varphi _{m,2j-m-k}^{l}\left( 1-{{2}^{j-l}}\xi \right),k={{2}^{j}}-2m-2,\cdots ,{{2}^{j}}-0 \\
\varphi _{m,0}^{l}\left( {{2}^{j-l}}\xi -{{2}^{-l}}k \right),k=0,1,\cdots ,{{2}^{j}}-m+1 \\
\end{matrix} \right.$
(12)
局部域内的小波紧支撑空间表示为
$\sup p\varphi _{m,k}^{j}\left( \xi \right)=\left\{ \begin{matrix}
\left[ k{{2}^{-j}},\left( 2m-1+k \right){{2}^{-j}} \right] \\
\left[ 0,\left( 2m-1+k \right){{2}^{-j}} \right] \\
\left[ k{{2}^{j}},1 \right] \\
\end{matrix} \right.$
(13)
文章选用三次B样条小波函数,即m=3。三次B样条小波函数的基本表达式为[17]
$\varphi \left( x \right)=\frac{1}{6}\left\{ \begin{array}{*{35}{l}}
{{\left( x+2 \right)}^{3}},x\in \left[ -2,-1 \right] \\
{{\left( x+2 \right)}^{3}}-4{{\left( x+1 \right)}^{3}},x\in \left[ -1,0 \right] \\
{{\left( 2-x \right)}^{3}}-4{{\left( 1-x \right)}^{3}},x\in \left[ 0,1 \right] \\
{{\left( 2-x \right)}^{3}},x\in \left[ 1,2 \right] \\
0,\left| x \right|\ge 2 \\
\end{array} \right.$
(14)
因无网格法形函数ΦT
(X)不满足克罗内克条件,对于刚度方程式(7)中的节点I的参数u并非其真实位移,所以本质边界条件不能够像有限元法那样直接施加。Chen J S等提出了全转换法(Full transformation method)[19]来处理其边界条件。本文按照全转换法的原理对刚度方程(控制方程)进行修正,得到
${{{\bar{K}}}_{I}}\bar{u}={{{\bar{F}}}_{I}}$
(15)
这样就可以得到节点的真实位移${\bar{u}}$,然后可以直接施加本质边界条件。
4 算例为了验证文章方法的正确性,本节给出船舶结构中“船体梁”和“开孔板”两种典型结构的计算分析。本文算例中材料的杨氏弹性模量为E=2.1×105 MPa,泊松比μ=0.3。
4.1 船体梁梁结构是船体基本结构之一。图1为一端固支、一端自由的船体简化梁结构。均布载荷p=10 N/m。
图2(a)是根据本文无网格算法经过4步自适应生成的离散节点图。初始离散为90个均布节点,积分背景网格为351个,经过4步自适应分析后,节点数为387个,背景网格为1 210个。图2(b)为ANSYS生成的有限元网格单元。
采用本文无网格方法,计算出算例梁的应力云图,如图3(a)所示。通过与图3(b)比较,应力集中区域表现都是一致的。
图4和图5给出了在小波基尺度scale=2.5下按照本文无网格法和有限元ANSYS计算结果和解析解比较曲线(图中FEM为有限元法解,ANS为解析解,MLPG为文章算法解,后面图例一致)。为了简单起见,仅对比板Y方向中线截面上正应力分布情况和板X方向中线截面上位移变化图。从图可以得出:1)FEM法和MLPG法的计算相对误差范数都在5%以下;2)文章提出的方法比FEM方法更接近解析解。可见本文方法具有很高的计算精度。
表1为在小波基不同尺度下的计算误差值表(图中Lu为位移相对误差范数,Lσ为应力相对误差范数,后面相同),根据表1可以绘制误差曲线图为图6所示。从表1和图6的变化趋势可以得出结论:1)计算精度受尺度变化影响,也就是支撑域的大小影响到计算精度;2)进一步提高文章方法的精度在于确定恰当的小波基尺度,可以由尺度三次B样条函数的尺度函数来调节。
尺度 | 2.0 | 2.1 | 2.5 | 2.8 | 3.0 | 3.5 | 4.0 |
Lu-B | 4.773 | 3.450 | 0.651 | 3.363 | 3.461 | 1.571 | 2.385 |
Lσ-B | 4.904 | 3.426 | 1.311 | 3.347 | 3.459 | 1.686 | 2.347 |
Lσ-G | 12.607 | 15.567 | 7.161 | 1.079 | 5.547 | 9.728 | 4.773 |
Lu-G | 13.093 | 15.678 | 7.224 | 2.139 | 6.038 | 9.742 | 4.905 |
4.2 开孔板
图7为简化为两端自由、另两端承受均匀拉力的平直船体板,板中间开有圆孔。此类型结构在船舶结构上也非常普遍,板侧均布拉力q=1 MPa。
图8(a)是根据本文无网格算法经过四步自适应生成的离散节点图。初始离散为749个均布节点,积分背景网格为1 390个,经过四步自适应分析后,节点数为1 323个,背景网格为2 515个。图8(b)为ANSYS生成的有限元网格单元。
采用本文无网格方法,计算出算例板的应力云图,如图9(a)所示。通过与有限元软件ANSYS计算的应力云图图9(b)所示比较,应力分布表现都是一致的。
进一步比较本文无网格法与有限元ANSYS计算的结果的相似程度。为了简单起见,仅对比板Y方向中线截面上正应力分布情况和板X方向中线截面上位移变化图。结果如图10和图11所示。FEM法的位移相对误差范数最大为2.5%、应力误差范数最大为2.4%,而文章方法的位移相对误差范数为0.3%、应力误差范数最大为0.28%,可见文章方法的计算精度具有明显的优势。
为了进一步了解文章方法的正确性,下面给出在尺度为2.1沿板中线X方向节点使用FEM法和文章方法计算的位移相对误差范数变化图,如图12所示。由图和前面的分析可以得出结论:1)文章方法具有很好的计算精度;2)尺度选定对文章算法有很大的影响,不同尺度下的计算精度有很大差别(如算例1中图6所示规律)。
5 结论文章提出的基于小波基的船舶结构无网格分析技术由于拥有较好的自适应性,在计算精度上比传统的有限元法具有优势。采用有限元法(ANSYS)和提出的方法对船体梁和开孔板进行进行实例计算,并对计算结果进行比较和分析,可以得出结论:
1)文章方法具有很好的计算精度;
2)尺度选定对文章算法有很大的影响,不同尺度下的计算精度有很大差别;
3)当3阶B样条小波函数作为权函数,结果权重函数比高斯权函数更好;
4)文章所提出的方法可以在小支撑域的情况下,达到良好的精度,即它可以实现更高的拟合。
[1] | 孙丽萍, 李力波. 船舶结构有限元分析[M]. 哈尔滨:哈尔滨工程大学出版社, 2013:1-12. SUN Liping, LI Libo. Finite element analysis of ship structures[M]. Harbin:Harbin Engineering University Press, 2013:1-12. |
[2] | BELYTSCHKO T, KRONGAUZ Y, ORGAN D, et al. Meshless methods:An overview and recent developments[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1-4):3-47. |
[3] | LIU W K, HAO S, BELYTSCHKO T, et al. Multiple scale meshless methods for damage fracture and localization[J]. Computational Materials Science, 1999, 16(1/2/3/4):197-205. |
[4] | ATLURI S N, SHEN Shengping. The basis of meshless domain discretization:the meshless local Petrov-Galerkin (MLPG) method[J]. Advances in Computational Mathematics, 2005, 23(1/2):73-93. |
[5] | 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. |
[6] | 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. |
[7] | 何沛祥, 李子然, 吴长春. 无网格与有限元的耦合在动态断裂研究中的应用[J]. 应用力学学报, 2006, 23(2):195-198. HE Peixiang, LI Ziran, WU Changchun. Coupled finite element-element-free Galerkin method for dynamic fracture[J]. Chinese Journal of Applied Mechanics, 2006, 23(2):195-198. |
[8] | 段念, 王文珊, 于怡青, 等. 基于FEM与SPH耦合算法的单颗磨粒切削玻璃的动态过程仿真[J]. 中国机械工程, 2013, 24(20):2716-2721. DUAN Nian, WANG Wenshan, YU Yiqing, et al. Dynamic simulation of single grain cutting of glass by coupling FEM and SPH[J]. China Mechanical Engineering, 2013, 24(20):2716-2721. |
[9] | 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. |
[10] | 胡德安, 韩旭, 肖毅华, 等. 光滑粒子法及其与有限元耦合算法的研究进展[J]. 力学学报, 2013, 45(5):639-652. HU Dean, HAN Xu, XIAO Yihua, et al. Research developments of smoothed particle hydrodynamics method and its coupling with finite element method[J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(5):639-652. |
[11] | LIU G R. Meshfree methods:moving beyond the finite element method[M]. Boca Raton:CRC Press, 2002:56-85. |
[12] | ATLURI S N, ZHU T. A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics[J]. Computational Mechanics, 1998, 22(2):117-127. |
[13] | ZIENKIEWICZ O C, TAYLOR R L. The finite element method[M]. 4th ed. London:McGraw-Hill, 1989:97-101. |
[14] | 杨玉英, 李晶. 无网格Galerkin方法中权函数的研究[J]. 塑性工程学报, 2005, 12(4):5-9. YANG Yuying, LI Jing. A study of weight function in element-free Galerkin method[J]. Journal of Plasticity Engineering, 2005, 12(4):5-9. |
[15] | LONG Shuyao, HU Dean. A study on the weight function of the moving least square approximation in the local boundary integral equation method[J]. Acta Mechanica Solida Sinica, 2003, 16(3):276-282. |
[16] | MALLAT S G. Multiresolution approximations and wavelet orthonormal bases of L2(R)[J]. Transactions of the American Mathematical Society, 1989, 315(1):69-87. |
[17] | 秦荣. 样条无网格法[M]. 北京:科学出版社, 2012:25-46.QIN Rong. Spline meshless method[M]. Beijing:Science Press, 2012:25-46. |
[18] | CHUI C K, QUAK E. Wavelets on a bounded interval[M]//BRAESS D, SCHUMAKER L L. eds. Numerical Methods of Approximation Theory. Basel:Birkhäuser Basel, 1993, 9:53-57. |
[19] | CHEN J S, PAN Chunhui, WU Chengtang, et al. Reproducing kernel particle methods for large deformation analysis of non-linear structures[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1/2/3/4):195-227. |