基于MLPG法的船舶結(jié)構(gòu)無網(wǎng)格自適應(yīng)分析技術(shù)

陳建平，唐文勇，徐曼平

(1.上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海200240;2.廣州航海學(xué)院 船舶工程學(xué)院，廣東 廣州510725)

0 引 言

當(dāng)前船舶結(jié)構(gòu)分析最常用也最為有效的計(jì)算工具之一就是有限元方法，有限元方法在分析、處理和模擬船舶結(jié)構(gòu)變形、應(yīng)力和疲勞分析等都有著廣泛和成熟的應(yīng)用［1］。但有限元法因其自身的特點(diǎn)，在處理分析船舶結(jié)構(gòu)場(chǎng)量(位移場(chǎng)和應(yīng)力場(chǎng)等)變化劇烈的高梯度區(qū)域時(shí)，會(huì)出現(xiàn)計(jì)算精度降低甚至計(jì)算中斷現(xiàn)象。為了解決這個(gè)問題，有限元法通常采用對(duì)該區(qū)域的網(wǎng)格進(jìn)行加密(細(xì)分)或者采用高階單元，這樣就要求方法具有較強(qiáng)的自適應(yīng)分析能力，其結(jié)果是加大了有限元法前后處理的工作量，從而降低了計(jì)算效率，事實(shí)上這種方法也并不能從根本上消除問題產(chǎn)生的根源。

作為與有限元法相對(duì)應(yīng)的另一種數(shù)值分析方法——無網(wǎng)格法，在近20年中得到了很大發(fā)展［2-6］。無網(wǎng)格法是建立在系列獨(dú)立離散點(diǎn)的基礎(chǔ)上，通過構(gòu)造點(diǎn)的近似函數(shù)來求解問題。與傳統(tǒng)有限元法相比，無網(wǎng)格法無需網(wǎng)格背景，在計(jì)算過程中可以根據(jù)需要任意增減節(jié)點(diǎn)，而不需要處理節(jié)點(diǎn)之間的拓?fù)湫畔ⅲ貏e適合用來進(jìn)行自適應(yīng)分析計(jì)算。另外由于無網(wǎng)格法沒有必要?jiǎng)?chuàng)建網(wǎng)格，節(jié)點(diǎn)可以由計(jì)算機(jī)以自主方式進(jìn)行創(chuàng)建，可以節(jié)省花費(fèi)在創(chuàng)建和處理網(wǎng)格的時(shí)間和計(jì)算資源。無網(wǎng)格法目前在航空材料、高速碰撞、動(dòng)態(tài)裂紋擴(kuò)展、加工成型、節(jié)理巖體分析等諸多領(lǐng)域都得到了較為廣泛的應(yīng)用［7-10］。

本文鑒于上述現(xiàn)狀和背景，提出基于局部無網(wǎng)格Petrov-Galerkin法的船舶板結(jié)構(gòu)無網(wǎng)格自適應(yīng)分析技術(shù)。首先將需要分析的船舶板結(jié)構(gòu)簡(jiǎn)化為平面板結(jié)構(gòu)，并將此結(jié)構(gòu)定義為問題域;再根據(jù)二維彈性理論，采用在離散節(jié)點(diǎn)上建立緊支試函數(shù)，運(yùn)用加權(quán)余量法建立系統(tǒng)的離散方程;然后由移動(dòng)最小二乘法(MLS)構(gòu)造離散節(jié)點(diǎn)的形函數(shù)，通過離散積分方程的等效弱形式來進(jìn)行求解。運(yùn)用Delaunay 三角化細(xì)分方案來執(zhí)行節(jié)點(diǎn)加密自適應(yīng)計(jì)算方案。本文的本質(zhì)邊界條件采用全轉(zhuǎn)化法的方法來處理。最后選取了2 種典型船舶結(jié)構(gòu)進(jìn)行應(yīng)力分析計(jì)算，并與有限元法的計(jì)算結(jié)果進(jìn)行比較來驗(yàn)證本文方法的有效性。

1 船舶結(jié)構(gòu)MLPG法表達(dá)式

由于船體板的板厚遠(yuǎn)小于其另外2個(gè)方向(長(zhǎng)度和寬度)的幾何尺寸，在分析其受力時(shí)，可以運(yùn)用Kirchhoff-Love板殼理論和Mindlin-Reissner板殼理論。本文基于板殼的Mindlin-Reissner的基本理論，運(yùn)用局部無網(wǎng)格Petrov-Galerkin法(MLPG)，來分析研究無網(wǎng)格法下的船舶板結(jié)構(gòu)應(yīng)力。

1.1 節(jié)點(diǎn)控制方程

板結(jié)構(gòu)經(jīng)離散后二維彈性理論為基礎(chǔ)的節(jié)點(diǎn)系統(tǒng)方程［11］

式中:i和j 為板的2個(gè)維度;σij，i為應(yīng)力σij在i 方向的分量;bi為板在i 方向的體積力;Ω 為問題域;Γ 為Ω的邊界;ui=;Γu為問題的本質(zhì)邊界條件;σijnj=;Γt為問題的自然邊界條件;nj為自然邊界上單位外法向矢量的第j個(gè)分量。

利用加權(quán)余量法［12］，可以得到節(jié)點(diǎn)I 處系統(tǒng)方程的微分方程強(qiáng)形式為

對(duì)節(jié)點(diǎn)I 依式(2)在其積分域上進(jìn)行積分，可以建立起它的系統(tǒng)方程。這樣對(duì)每個(gè)節(jié)點(diǎn)都采用式(2)在其積分域上積分，可以得到所有離散節(jié)點(diǎn)的系統(tǒng)方程，將這些離散節(jié)點(diǎn)的系統(tǒng)方程組裝起來就能夠獲得問題域的整體系統(tǒng)方程。

利用移動(dòng)最小二乘法(MLS)得到節(jié)點(diǎn)的積分點(diǎn)支撐域內(nèi)的形函數(shù)，從而獲得問題域的位移逼近函數(shù)［12］。

式中:ΦT(X)為根據(jù)MLS所得的形函數(shù)矩陣;uI為離散節(jié)點(diǎn)I的節(jié)點(diǎn)值;N 為積分點(diǎn)的支撐域Ωs中的節(jié)點(diǎn)數(shù)。

根據(jù)彈性力學(xué)應(yīng)力-應(yīng)變關(guān)系有

國(guó)內(nèi)外已有使用激光進(jìn)行高速目標(biāo)測(cè)速的相關(guān)研究，常用方法大致可分為三類：一是基于直接光譜技術(shù)，利用法布里-珀羅干涉儀直接測(cè)量譜線的頻移[5]；二是雙激光拍頻技術(shù)，利用兩個(gè)線偏振光同時(shí)傳感物體的速度，可大大提高最高可測(cè)量速度[5-10]；三是利用外調(diào)制技術(shù)，使用微波作為模擬調(diào)制信號(hào)對(duì)激光進(jìn)行強(qiáng)度調(diào)制，通過檢測(cè)微波調(diào)制信號(hào)的多普勒頻率，來實(shí)現(xiàn)相對(duì)運(yùn)動(dòng)速度的測(cè)量[11-13]。相比于雙激光拍頻技術(shù)，采用外調(diào)制技術(shù)降低了對(duì)激光器線寬和穩(wěn)定度的要求，更容易實(shí)現(xiàn)高精度和高動(dòng)態(tài)的測(cè)量。2016年，郝文澤等人采取強(qiáng)度調(diào)制和直接探測(cè)的方式，成功獲得了可靠的測(cè)速實(shí)驗(yàn)數(shù)據(jù)[13]。

式中:D 為材料彈性矩陣;B 為幾何矩陣，且Bj=

式中:ΓQt為積分域ΩQ和自然邊界的重合部分;N為自然邊界上的單位外法向矢量矩陣。

1.2 剛度方程

把節(jié)點(diǎn)系統(tǒng)方程式(3)和式(4)代入節(jié)點(diǎn)控制方程式(5)，可以得到

式(6)可簡(jiǎn)記為矩陣形式

其中

這樣，由每個(gè)離散節(jié)點(diǎn)的控制方程總裝成整個(gè)結(jié)構(gòu)系統(tǒng)的控制方程表達(dá)為

由于在MLS 中，采用的節(jié)點(diǎn)支撐域都是緊支的，系統(tǒng)的總體剛度矩陣K 是帶狀稀疏矩陣，可以減少計(jì)算負(fù)荷。根據(jù)文獻(xiàn)［13］對(duì)有限元法中罰因子的分析，建議罰因子α 在范圍內(nèi)105～12× E 取值，E為楊氏彈性模量。

2 自適應(yīng)細(xì)分計(jì)算

圖1 細(xì)分方案示意圖Fig.1 Refinement scheme by local Delaunay algorithm

算法流程為:

1)影響域三角化:在確定不符合精度要求的節(jié)點(diǎn)i 后，在其影響域內(nèi)搜索相鄰節(jié)點(diǎn)并利用這些節(jié)點(diǎn)進(jìn)行三角化;

2)插入新節(jié)點(diǎn):遍歷環(huán)繞節(jié)點(diǎn)i的所有三角形，在三角形各條邊的中點(diǎn)上加入新的節(jié)點(diǎn);

3)刪除冗余節(jié)點(diǎn):刪除重復(fù)節(jié)點(diǎn)，使得每個(gè)位置只有1個(gè)節(jié)點(diǎn)。若任意2個(gè)節(jié)點(diǎn)允許的最小距離為dm，判斷兩新節(jié)點(diǎn)間距離小于dm，則需刪除其中一個(gè)節(jié)點(diǎn)。

4)刪除三角形:在局部區(qū)域內(nèi)刪除三角形的邊。

自適應(yīng)計(jì)算誤差確定采用如下的相對(duì)誤差估計(jì)方法［14］:

當(dāng)相對(duì)誤差值小于給定的計(jì)算誤差值，該節(jié)點(diǎn)計(jì)算自動(dòng)終止，轉(zhuǎn)入下一節(jié)點(diǎn)的計(jì)算，并執(zhí)行下一節(jié)點(diǎn)的自適應(yīng)計(jì)算方案。

3 本質(zhì)邊界條件施加

因無網(wǎng)格法形函數(shù)ΦT(X)不滿足克羅內(nèi)克條件，對(duì)于剛度方程式(7)中的節(jié)點(diǎn)I的參數(shù)u 并非其真實(shí)位移，所以本質(zhì)邊界條件不能夠像有限元法那樣直接施加。Chen J S 等提出了全轉(zhuǎn)換法(Full transformationmethod)［15］來處理其邊界條件。本文按照全轉(zhuǎn)換法的原理對(duì)剛度方程(控制方程)進(jìn)行修正得到

4 算 例

為了驗(yàn)證文章方法的正確性，本節(jié)給出船舶結(jié)構(gòu)中“實(shí)肋板”和“裂紋板”2 種典型結(jié)構(gòu)的分析計(jì)算。本文算例中材料的楊氏彈性模量為E=2.1 ×105MPa，泊松比μ=0.3。

4.1 實(shí)肋板

中間開有減輕孔的實(shí)肋板結(jié)構(gòu)是船舶普通結(jié)構(gòu)之一。圖2 簡(jiǎn)化為上下兩端自由、另兩端承受均勻拉力的船體實(shí)肋板。板側(cè)均布拉力q=1 MPa。

圖2 典型實(shí)肋板示意圖Fig.2 The typical solid floor

根據(jù)本文算法，采用不規(guī)則初始離散節(jié)點(diǎn)245個(gè)。經(jīng)過6 步自適應(yīng)計(jì)算，每個(gè)自適應(yīng)步離散圖和計(jì)算出的應(yīng)力云圖如圖3所示。

各自適應(yīng)步在實(shí)肋板沿Y 方向?qū)ΨQ軸上的應(yīng)力比較(考慮X 方向也對(duì)稱，只取一半比較)如圖4所示。從圖中可看出，曲線變化趨勢(shì)存在局部反彈，但總體趨勢(shì)是在收斂，且收斂方向趨于一致。

圖3 自適應(yīng)離散圖和應(yīng)力云圖Fig.3 The discrezition figure and fringe firgure

圖4 Y 各自適應(yīng)步方向軸線上應(yīng)力變化圖Fig.4 The curve of Mises stress along Y-Direction inmidsection

圖5 為本文方法計(jì)算結(jié)果與有限元軟件Ansys和Nastran 以及用解析法計(jì)算的結(jié)果比較圖。從圖示各圖線來看，本文方法比在靠近板邊緣和孔的邊緣附近處，比有限元法有著更好的精度。

圖5 各種方法計(jì)算應(yīng)力結(jié)果比較圖Fig.5 The curve of Mises stress with different ways

各自適應(yīng)步在實(shí)肋板沿Y 方向?qū)ΨQ軸上的應(yīng)力相對(duì)誤差比較(考慮X 方向也對(duì)稱，只取一半比較)如圖6所示。由圖中可看出，相對(duì)誤差雖然存在局部震蕩情況，但每一步變化趨勢(shì)都是在減小，說明計(jì)算精度在逐步提高，說明文章方法具有很好的計(jì)算精度。

圖6 相對(duì)誤差變化圖Fig.6 The curve of the relative error

4.2 裂紋板

一長(zhǎng)邊為1 m，短邊為0.5 m，在長(zhǎng)邊中部邊緣沿Y 方向有一裂紋長(zhǎng)為100 mm。板上下兩端為自由邊，側(cè)邊受均布拉力p=10 N/m。

根據(jù)本文算法，采用不規(guī)則初始離散節(jié)點(diǎn)137個(gè)。經(jīng)過5 步自適應(yīng)計(jì)算，每個(gè)自適應(yīng)步離散圖和計(jì)算出的應(yīng)力云圖如圖7所示。

圖7 自適應(yīng)離散圖和應(yīng)力云圖Fig.7 The discrezition figure and fringe firgure

各自適應(yīng)步從邊緣到裂紋尖端沿Y 方向的應(yīng)力比較如圖8所示。從圖中可看出，曲線在兩端變化較大，但總體趨勢(shì)是在收斂，且收斂方向趨于一致。

各自適應(yīng)步沿裂紋方向的應(yīng)力相對(duì)誤差比較如圖9所示。由圖可知，相對(duì)誤差雖然存在局部震蕩情況，但每一步變化趨勢(shì)都在減小，說明計(jì)算精度在逐步提高，說明本方法具有良好的計(jì)算精度。

圖8 沿裂紋方向應(yīng)力變化曲線Fig.8 The curve of the Mises stress along the crack line

圖9 相對(duì)誤差變化圖Fig.9 The curve of the relative error

5 結(jié) 語(yǔ)

本文基于無網(wǎng)格局部MLPG法，提出了船舶結(jié)構(gòu)的無網(wǎng)格分析方法，并采用Delaunay 三角化方法對(duì)節(jié)點(diǎn)進(jìn)行自動(dòng)加密計(jì)算。最后運(yùn)用所提方法對(duì)典型的船舶結(jié)構(gòu)——船體實(shí)肋板和裂紋板的分析計(jì)算，并通過與有限元法的計(jì)算結(jié)果進(jìn)行比較，可以看出本文提出的無網(wǎng)格自適應(yīng)分析方法對(duì)于求解船舶結(jié)構(gòu)的變形(位移)和應(yīng)力可行，并且具有良好的精度，驗(yàn)證了本文方法的有效性和準(zhǔn)確性。

［1］孫麗萍，李力波.船舶結(jié)構(gòu)有限元分析［M］.哈爾濱:哈爾濱工程大學(xué)出版社，2013.

［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:3-47.

［3］LIU W K，HAO S，BELYTSCHKO T，et al.Multiple scale meshless methods for damage fracture and localization［J］.Comput.Mater.Sci.，1999，16:197-206.

［4］ATLURI S N，SHEN S.The basis of meshless domain discretization: the meshless local Petrov- Galerkin(MLPG) method ［J］.Advances in Computational Mathematics，2005，23:73-93.

［5］ODEN J T，DUARTE C A，ZIENKIEWICZ O C.A new could-based hp finite element method［J］.Int.J.Num.Meth.Engng.，1998，50:160-170.

［6］LIU G R，GU Y T.Meshless local Petrov-Galerkinmethod in combination with finite element and boundary element approaches［J］.Computational Mechanics，2000，26:536-646.

［7］何沛祥，李子然，吳長(zhǎng)春.無網(wǎng)格與有限元的耦合在動(dòng)態(tài)斷裂研究中的應(yīng)用［J］.應(yīng)用力學(xué)學(xué)報(bào)，2006，23(2):195-198.HE Pei-xiang，LI Zi-ran，WU Chang-chun.Coupled finite element-element-free galerkinmethod for dynamic fracture［J］.Chinese Journal of Applied Mechanics，2006，23(2):195-198.

［8］段念，王文珊，于怡青，等.基于FEM 與SPH 耦合算法的單顆磨粒切削玻璃的動(dòng)態(tài)過程仿真［J］.中國(guó)機(jī)械工程，2013，24(20):2716-2721.DUAN Nian，WANG Wen-shan，YU Yi-qing，et al.Dynamic simulation of single grain cutting of glass by coupling FEM and SPH ［J］.Chinese Mechanical Engineering，2013，24(20):2716-2721.

［9］JOHNSON R G，STRYK A R，BEISSEL R S，et al.An algorithm to automatically convert distorted finite element into meshless particles during dynamic deformation［J］.International Journal of Impact Engineering，2002，27:997-1013.

［10］胡德安，韓旭，肖毅華，等.光滑粒子法及其與有限元耦合算法的研究進(jìn)展［J］.力學(xué)學(xué)報(bào)，2013，45(5):639-652.HU De-an，HAN Xu，XIAO Yi-hua，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):39-652.

［11］LIU G R.Meshfree methods moving beyond the finite element method［M］.CRC Press，2003.

［12］ATLURI S N，ZHU T.A new meshless local Petrov-Galerkin(MLPG)approach in computational mechanics［J］.Computational Mechanics，1998，22:117-179.

［13］ZIENKIEWICZ O C.The finite element method［C］//4th ed.，McGraw-Hill，London，1989.

［14］GAVETE L，F(xiàn)ALCON S，RUIZ A.An error indicator for the element free Galerkinmethod ［J］.European Journal Mechanics-A/Solids，2001，20(3):327-341.

［15］周培德.計(jì)算幾何［M］.北京:清華大學(xué)出版社，2000.