陳建平，唐文勇，徐曼平

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

0 引 言

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

作為與有限元法相對應(yīng)的另一種數(shù)值分析方法——無網(wǎng)格法，在近20年中得到了很大發(fā)展［2-6］。無網(wǎng)格法是建立在系列獨立離散點的基礎(chǔ)上，通過構(gòu)造點的近似函數(shù)來求解問題。與傳統(tǒng)有限元法相比，無網(wǎng)格法無需網(wǎng)格背景，在計算過程中可以根據(jù)需要任意增減節(jié)點，而不需要處理節(jié)點之間的拓?fù)湫畔ⅲ貏e適合用來進(jìn)行自適應(yīng)分析計算。另外由于無網(wǎng)格法沒有必要創(chuàng)建網(wǎng)格，節(jié)點可以由計算機以自主方式進(jìn)行創(chuàng)建，可以節(jié)省花費在創(chuàng)建和處理網(wǎng)格的時間和計算資源。無網(wǎng)格法目前在航空材料、高速碰撞、動態(tài)裂紋擴展、加工成型、節(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é)構(gòu)，并將此結(jié)構(gòu)定義為問題域;再根據(jù)二維彈性理論，采用在離散節(jié)點上建立緊支試函數(shù)，運用加權(quán)余量法建立系統(tǒng)的離散方程;然后由移動最小二乘法(MLS)構(gòu)造離散節(jié)點的形函數(shù)，通過離散積分方程的等效弱形式來進(jìn)行求解。運用Delaunay 三角化細(xì)分方案來執(zhí)行節(jié)點加密自適應(yīng)計算方案。本文的本質(zhì)邊界條件采用全轉(zhuǎn)化法的方法來處理。最后選取了2 種典型船舶結(jié)構(gòu)進(jìn)行應(yīng)力分析計算，并與有限元法的計算結(jié)果進(jìn)行比較來驗證本文方法的有效性。

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

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

1.1 節(jié)點控制方程

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

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

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

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

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

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

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

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

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

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

1.2 剛度方程

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

式(6)可簡記為矩陣形式

其中

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

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

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

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

算法流程為:

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

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

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

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

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

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

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

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

4 算 例

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

4.1 實肋板

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

圖2 典型實肋板示意圖Fig.2 The typical solid floor

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

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

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

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

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

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

4.2 裂紋板

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

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

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

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

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

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

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

5 結(jié) 語

本文基于無網(wǎng)格局部MLPG法，提出了船舶結(jié)構(gòu)的無網(wǎng)格分析方法，并采用Delaunay 三角化方法對節(jié)點進(jìn)行自動加密計算。最后運用所提方法對典型的船舶結(jié)構(gòu)——船體實肋板和裂紋板的分析計算，并通過與有限元法的計算結(jié)果進(jìn)行比較，可以看出本文提出的無網(wǎng)格自適應(yīng)分析方法對于求解船舶結(jié)構(gòu)的變形(位移)和應(yīng)力可行，并且具有良好的精度，驗證了本文方法的有效性和準(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］何沛祥，李子然，吳長春.無網(wǎng)格與有限元的耦合在動態(tài)斷裂研究中的應(yīng)用［J］.應(yīng)用力學(xué)學(xué)報，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 耦合算法的單顆磨粒切削玻璃的動態(tài)過程仿真［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é)報，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］周培德.計算幾何［M］.北京:清華大學(xué)出版社，2000.