駱寒冰，董德龍，季紅葉，趙中榜，高 原

（天津大學 a.水利工程仿真與安全國家重點實驗室；b.船舶與海洋工程系，天津 300072）

基于ALE算法的無轉角和有轉角鋁制加筋板楔形體水彈性砰擊的數值模擬研究

駱寒冰a，b，董德龍b，季紅葉b，趙中榜b，高 原b

（天津大學 a.水利工程仿真與安全國家重點實驗室；b.船舶與海洋工程系，天津 300072）

文章采用基于任意拉格朗日—歐拉（ALE）算法的顯式有限元技術研究水彈性砰擊現象，針對已開展的鋁制加筋板楔形體結構入水砰擊模型實驗，開展了數值模擬比較工作。該楔形體底部斜升角為20度，底部兩側是包含三根縱骨和兩根橫梁的加筋板結構，兩側結構剛度不同。預報了模型無轉角和有轉角典型工況的砰擊入水過程，得到的入水加速度、底部加筋板結構縱骨應力和橫梁響應與模型實驗結果吻合較好。研究表明該ALE算法具備模擬船舶局部結構的水彈性砰擊流固耦合問題的能力。

ALE算法；鋁制加筋板結構；水彈性砰擊；轉角；數值計算

0 引 言

在波浪中高速航行船舶的首部容易遭受嚴重的水動力砰擊作用[1]。準確預報砰擊載荷并合理確定其設計載荷是艏部結構設計的一個難題。砰擊載荷在復雜船首外表面的時間和空間分布很復雜。船首局部加筋板結構設計時候，通常把這復雜的水動力載荷簡化為準靜態(tài)的、均勻分布的載荷處理，這與實際的時空分布相比差異很大。而且，當局部結構的斜升角較小時，或者砰擊載荷的持續(xù)時間和局部結構的固有周期接近時，Faltinsen[2]指出，需要考慮局部結構的水彈性效應，砰擊問題應當作水動力學和結構力學的聯合體來研究，即需要同時分析砰擊載荷和結構動響應的時間和空間分布特性。

由于砰擊問題的復雜性，采用解析或半解析方法研究水彈性砰擊機理時候，通常把復雜結構簡化成二維的楔形體、圓柱體等，分析結構落體到平靜水面的典型入水問題。通常流體域采用勢流理論，結構域采用有限元（FEM）方法，來考慮結構和流體兩者之間耦合作用。如盧熾華等[3]結合邊界元方法（BEM）和FEM法，分析二維彈性楔形體砰擊，Sun和Faltinsen[4]結合BEM和結構解析方法，研究圓柱殼體水彈性砰擊。Korobkin等[5]采用Wagner[6]理論和梁有限元方法，分析二維彈性楔形體的水彈性效應。除了典型的回轉體這類簡單結構之外，目前對三維水彈性砰擊研究尚不多見。如Peseux等[7]采用三維Wagner理論和FEM法，分析了圓錐體砰擊問題。

近年來，基于顯式有限元技術的數值方法開始應用于砰擊問題。Stenius等[8]應用顯式有限元LSDYNA軟件的ALE算法研究二維楔形體剛體入水砰擊問題，Stenius等[9]還探討了二維彈性體的水彈性效應。Aquelet等[10]討論了LS-DYNA軟件預報砰擊時罰函數耦合算法對計算結果的影響。Luo等[11]針對二維剛體砰擊討論了網格大小對數值模擬結果精度的影響。Luo等[12]選取船舶最基本的結構—加筋板結構，應用ALE方法數值分析了某加筋板楔形體入水砰擊問題，預報的加速度和結構應力響應與實驗吻合較好。由于需要描述局部流體砰擊射流，LS-DYNA數值模擬時候要求流體的網格很小，網格才不會畸變，導致時間步長變得很小，模型的網格數量大，在微機上計算時間很長，駱寒冰等[13]研究了并行計算技術來提高計算效率。上述研究表明，基于ALE算法的顯式有限元技術在砰擊預報上有很好的應用前景，具備預報結構水彈性砰擊的潛力，但是，如果要推廣應用到船舶上，還需要繼續(xù)做數值驗證工作。

目前，針對有轉角結構水彈性砰擊的研究還較少見，而船舶在首斜浪航行時候，會有明顯的橫搖運動，此時船首入水砰擊的斜升角會比頂浪時候的要小，角度小時候水彈性效應可能更明顯，因此，研究有轉角工況具有重要的工程應用價值。本文根據最近在天津大學進行的鋁制加筋板楔形體入水模型實驗[14]，應用基于ALE算法的顯式有限元技術，開展數值預報工作，分析在無轉角和有轉角情況下模型結構的入水砰擊過程，預報砰擊加速度和結構應力響應，與模型實驗結果進行對比分析。通過本文的研究，驗證算法在上述問題上的適用性，為今后應用本方法預報實際船舶結構水彈性砰擊做技術儲備。

1 基于ALE算法的顯式有限元技術

ALE算法最早出現于數值模擬流體動力學問題的有限差分方法中。該算法同時具備Lagrange方法和 Euler方法兩者的特長[15-16]。引進了Lagrange方法的特點來處理結構邊界運動，能夠有效地跟蹤物質結構邊界的運動；吸收了Euler方法的長處來劃分內部網格，使內部網格單元獨立于物質實體而存在，不同于Euler方法的是，ALE方法網格可以根據定義的參數在求解過程中適當地調整位置，使得網格不致出現嚴重的畸變。這種方法在處理整個物體有空間的大位移并且本身有大變形的問題時非常有優(yōu)勢，可以克服單元嚴重畸變引起的數值計算困難，并可以實現流體—固體耦合的動態(tài)分析。

在軟件LS-DYNA[17]中水和空氣采用空白材料模型與狀態(tài)方程聯用表示，多物質歐拉網格，自由表面的生成利用VOF方法，結構采用彈性材料模型，拉格朗日網格。在歐拉流體和拉格朗日固體接觸面，采用罰函數耦合算法，可以模擬流體與固體之間的耦合效應。罰函數偶合系數追蹤拉格朗日節(jié)點（結構，從物質）和歐拉流體（主物質）物質位置間的相對位移d，檢查每一個主從節(jié)點對主物質表面的貫穿，如果從節(jié)點不出現貫穿，就不進行任何操作；如果發(fā)生從節(jié)點對主物質表面的貫穿，界面力F就會分布到歐拉流體的節(jié)點上。界面力的大小與發(fā)生的貫穿數量成正比：

式中：ki表示基于主從節(jié)點質量模型特性的剛度系數。

由于采用了上述任意的拉格朗日—歐拉（ALE）算法和罰函數耦合算法，LS-DYNA具備模擬水彈性砰擊流固耦合問題的能力。

2 鋁制加筋板楔形體入水砰擊實驗介紹

設計了一個鋁制三維加筋板楔形體模型，其主尺度長×寬×高為：1.440m×1.354m×0.746m。圖1是模型典型剖面結構示意圖。模型底部斜升角為20°，模型結構質量180 kg。底部兩側是包含三根‘1’字形縱骨和兩根‘T’字形橫梁的加筋板結構，是實驗測試區(qū)域。底部左右兩側的外板厚度分別為4 mm、3mm，兩側的縱骨和橫梁的尺寸及厚度不同，通過這樣的設計，在一次落體砰擊實驗中可以測試得到兩組不同剛度結構的應變響應結果。模型上部結構剛度設計得比較大，采用的鋁板厚度為8mm，目的是給底部加筋板結構提供一個剛性的四周邊界條件。模型前后端的擋板向楔形體下方做了適當延伸，目的是減少落體模型入水砰擊的三維效應。

圖1 模型典型橫剖面結構示意圖Fig.1 Mid-section structure of themodel

圖2 實驗模型測點布置示意圖Fig.2 Distribution ofmeasuring points on themodel

圖3 無轉角入水砰擊過程典型照片（高度1.20m）Fig.3 Photos forwater entry without roll angle(height1.20m)

模型實驗主要測量入水砰擊過程的落體加速度、底部外板的砰擊壓力，底部加筋板結構的應變動響應等參數。測點布置見圖2。模型縱桁前后布置2個加速度測點（以A為首字母），測量垂向入水加速度；底部外板板格的中心位置共布置了8個砰擊壓力測點（以P為首字母），測量垂直于外板的砰擊壓力；在縱骨上表面、縱向中間的位置粘貼18個應變片，在橫梁面板上表面、橫向粘貼8個應變片，合計26個應變測點（以S為首字母）。實驗總共36個測點。

在天津大學建造了一個實驗水池，其主尺度長×寬×高為：5.00m×4.40m×1.30m。實驗時候水池水深為1.20m。在實驗水池上方設計安裝了一個門字架，在門字架內側垂直安裝了兩個不銹鋼圓柱導軌。模型的前后檔板上部分別固定兩個滑動軸承，滑動軸承套在圓柱導軌上，滑動軸承的摩擦阻力可以忽略不計，于是利用該導軌可以確保模型是垂直入水的。在門字架一側安裝了一套手動葫蘆和滑輪組，用于人工升降模型。設計了一套釋放/抱緊裝置，一端與模型連接，另一端與鋼絲繩連接。通過人工啟動釋放/抱緊裝置的釋放開關，可以實現瞬間模型自由下落入水砰擊。通過調整安裝在模型前后檔板上滑動軸承的位置，模型還可以旋轉一定的角度，進行不同角度的入水砰擊實驗。

對于無轉角砰擊實驗工況，落體高度分別為0.30 m、0.50 m、0.75 m、1.00 m、1.20 m和 1.40 m。有轉角實驗工況包括左轉5°和右轉5°兩大類，落體高度分別為0.30 m、0.50 m和0.75 m。每組工況進行不少于3次的重復實驗。實驗數據采樣頻率10 kHz。圖3顯示了1.20m高度時候無轉角入水砰擊實驗的典型照片，觀測到了明顯的射流。圖4顯示了0.50m高度、右轉動5°時候模型實驗照片。

圖4 模型右轉5°時落體實驗照片（高度0.50m）Fig.4 Photos forwater entrymodelswith roll angle 5 degrees starboard side(height0.50m)

3 無轉角入水砰擊結果比較

對于上述鋁制三維加筋板楔形體的入水實驗，建立計算模型，采用LS-DYNA軟件中的ALE算法和罰函數耦合算法，數值模擬水彈性砰擊過程。由于模型四周和上部結構的板比底部板厚，為提高流固耦合計算效率，只考慮加筋板結構為彈性，分析其與流體的耦合問題，模型四周和上部結構簡化處理為剛體。為了減少計算時間，初始入水高度設為0.055m，可以充分減少計算時間，根據公式 ，可以確定不同高度下落時對應的模型的初始速度。數值模擬實際實驗的時間長度為0.15 s。

圖5 無轉角入水砰擊LS-DYNA計算3D模型 Fig.5 3Dmodelwithout roll angle forwater entry in LS-DYNA

圖6 無轉角入水砰擊流體域網格劃分Fig.6 Meshes in air and water domain ofmodel without roll angle

為了提高計算精度，建立了左右兩個模型，即根據對稱面X=0，針對楔形體左、右側底部加筋板結構，分別建立了一半模型。對稱面采用了對稱條件，無反射邊界條件施加于除對稱面以外的外部邊界。整個落體模型質量為190 kg，包括結構質量180 kg，軸承和螺栓共6 kg，測量傳感器及導線估算4 kg，那么一半模型的質量為95 kg。圖5顯示了無轉角入水砰擊的左側的LS-DYNA計算模型。

流體區(qū)域采用非均勻網格。距離楔形體近的網格為均勻布置，其它區(qū)域，距離楔形體越遠，網格尺寸越大。為能夠合理反映結構應變響應，縱骨在垂向上劃分了4個單元。水和空氣的流體域單元數量總共489 280個，底部外板厚度為3 mm的右側模型結構單元數量7 264個，底部外板厚度為4 mm的左側模型結構單元數量8 120個。圖6顯示了左側計算模型的流體域網格劃分情況。

選取落體高度為1.20m典型工況，對數值計算和模型實驗的入水加速度和結構動響應結果進行比較。圖7顯示了加速度時域曲線比較情況，其中，實驗加速度開始突變時刻34.084 5 s，作為結構開始入水砰擊的時間，“A1-Test”表示實驗A1測點的結果，而“A-Num.”表示加速度數值模擬結果。圖8顯示了結構應力時域曲線比較情況。S13、S12和S11是左側縱骨上縱向應力測點，分別對應右側測點S14、S15和S16，S17和S18是左側橫梁面板上橫向應力測點，對應右側測點為S20、S19。圖8中，“S13-Test”表示實驗S13測點的結果，“S13-Num.”表示S13數值模擬結果，其它標注含義以此類推。從圖7、8比較中可以看出：

圖7 砰擊加速度實驗與數值計算比較（高度1.20m）Fig.7 Comparison of the impactacceleration(height1.20m)

（1）模型實驗和數值模擬的加速度和和結構應力響應在時域上吻合很好，正確反映了整個入水過程，首先在時間34.084 5 s加速度突變，模型開始入水，接著應變測點信號從模型底部向上陸續(xù)出現響應；

（2）模型前后加速度測點A1、A2在變化趨勢和極值響應上吻合好，峰值出現的時間分別是34.095 2s、34.095 4 s，大小分別是135.8 m/s2、129.4m/s2，可以推斷出模型前后是同時入水的；數值預報的加速度峰值為139.3m/s2，略大于實驗測試結果，出現的時間是34.094 2 s，略早于實驗時間；

圖8 砰擊結構應力實驗與數值計算比較（高度1.20m）Fig.8 Comparison of stress responses(height1.20m)

（3）左側結構應力比右側的要小，這是由于左側底部加筋板結構的剛度比右側的大，在受到相同的外界壓力載荷情況下，左側縱骨和橫梁上的應力響應就小一些；縱骨應力是我們關注的重點，比橫梁上應力要大；

（4）靠近模型底部的縱骨應力最大。對于左側底部結構，縱骨測點S13應力響應峰值最大，實驗和數值模擬分別為63.7 MPa、59.7 MPa，出現的時間分別為34.096 1 s、34.096 7 s；對于右側底部結構，縱骨測點S14應力響應峰值最大，實驗和數值模擬分別為102.3 MPa、98.8 MPa，出現的時間分別為34.097 2 s、34.097 7 s；

（5）左右兩側結構對應測點位置的應力響應出現峰值、增大和減少變化趨勢相同，反映實驗入水過程保持左右對稱垂直入水，左右斜升角一樣；

（6）實驗縱骨測點的應力在峰值后觀測到一定的高頻振蕩，這是水彈性砰擊的表現之一，不過，數值預報的振蕩不太明顯，在后續(xù)的數值研究中將繼續(xù)探討該問題。

4 有轉角入水砰擊結果比較

對于有轉角入水問題，選取了落體高度0.50m、右轉5°工況開展比較分析。建立了整個落體結構的計算模型g，模型水和空氣的流體域單元數量總共544 640個，模型結構單元數量16 931個。流體域單元總數與前面無轉角模型的相當，由于無轉角模型只取了一半的流體域，所以本部分的模型網格密度比前部分的略稀疏些。外部邊界施加無反射邊界條件。圖9顯示了該數值計算模型的典型橫剖面。結構質量設為190 k。初始入水高度設為0.055 m，數值模擬實際實驗的時間長度為0.15 s。

圖10、11分別顯示了加速度時域曲線、結構應力時域曲線比較情況。從前部分無轉角入水砰擊結構應力結果來看，橫梁上面板應力較小，本部分就不討論橫梁的應力情況。通過比較可以看出：

（1）對有轉角情況，實驗和數值模擬的加速度和和結構應力響應在時域上吻合很好。

（2）模型前后加速度測點A1、A2峰值出現的時間分別是10.848 3 s、10.849 7 s，大小分別是53.5m/s2、56.8m/s2，表明實驗模型前后是同時入水的；數值預報的加速度峰值為61.4m/s2，出現的時間是10.849 6 s；

圖9 右轉5°入水砰擊LS-DYNA計算模型Fig.9 Modelwith roll angle 5 degrees starboard for water entry in LS-DYNA

圖10 砰擊加速度實驗與數值計算比較（高度0.50m，右轉5°）Fig.10 Comparison of the impactacceleration(height 0.50m，roll angle 5 degrees starboard side)

（3）由于該模型右轉5°，相當于左側底部加筋板斜升角為25°，右側的斜升角為15°，于是，入水過程中右側的砰擊壓力大于左側的，斜升角小一側的濕長度相對大一些。對于左右對稱的縱骨結構應力測點，右側S14測點比左側S13的先產生應力動響應，S14應力出現峰值的時間也比S13的早，同樣適用于右側S15與左側S12、右側S16與左側S11之間的關系；

（4）對于左側底部結構，縱骨測點S13應力響應峰值最大，實驗和數值模擬分別為17.9MPa、18.7 MPa，出現的時間分別為10.855 4 s、10.855 0 s；對于右側底部結構，縱骨測點S14應力響應峰值最大，實驗和數值模擬分別為61.7 MPa、65.5MPa，出現的時間分別為10.850 0 s、10.849 3 s；

（5）圖10數值模擬的加速度毛刺相對圖7的多，入水砰擊之前毛刺比較多，入水之后少一些，可能是加速度對模型流體單元網格密度要求高，而本部分模型的密度略小一些；相比較，無轉角和有轉角模型的結構應力預報的結果都很穩(wěn)定。

圖11 砰擊結構應力實驗與數值計算比較（高度0.50m，右轉5°）Fig.11 Comparison of stress responses(height 0.50m,roll angle 5 degrees starboard)

5 結 論

根據近期開展的鋁制加筋板楔形體結構入水砰擊模型實驗的結果，本文開展了數值模擬比較工作。采用基于ALE算法和罰函數耦合算法的顯式有限元技術，針對無轉角和有轉角入水砰擊典型工況，數值模擬了入水過程中的加速度、縱骨上端縱向應力和橫梁上面板橫向應力響應，并與實驗結果進行了比較，取得了以下結論：

（1）為了提高精度，減少計算時間，建立計算模型的時候，模型的底部加筋板結構采用了彈性材料，均勻網格，其它結構為剛性材料；流體區(qū)域采用了非均勻網格技術，距離楔形體底部近的，采用均勻布置，其它區(qū)域，距離越遠，網格密度越?。?/p>

（2）對于無轉角入水砰擊工況，分別建立了左右兩個模型進行分析。左右側底部加筋板結構對稱測點位置的應力響應在時間上的變化趨勢相同，由于右側的剛度比左側的要小，因此右側的應力大于左側的；

（3）對于右轉5°入水砰擊工況，建立了一個整體模型進行分析。由于右側底部斜升角小于左側的，那么右側結構對應縱骨測點應力響應峰值出現的時間早于左側的。

（4）模型實驗縱骨應力在峰值后觀測到一定的高頻振蕩，而數值模擬的振蕩現象不明顯，今后需要繼續(xù)深入研究該現象。

（5）數值模擬的加速度和應力響應與模型實驗的結果在時域上吻合較好，這表明本文提出的水彈性砰擊計算技術的合理性和可靠性，能夠預報無轉角和有轉角入水砰擊問題，具備推廣到預報船首結構水彈性砰擊的潛力。

參 考 文 獻：

[1]駱寒冰,徐 慧,余建星,萬正權.艦船砰擊載荷及結構動響應研究綜述[J].船舶力學,2010,14(4):439-450. Luo Hanbing,XüHui,YüJianxing,Wan Zhengquan.Review of the state of the art of dynamic responses induced by slamming loads on ship structures[J].Journal of Ship Mechanics,2010,14(4):439-450.(in Chinese)

[2]Faltinsen OM.Hydroelastic slamming[J].Journal ofMarine Science and Technology,2000,5(2):49-65.

[3]盧熾華,何友聲.二維彈性結構入水沖擊過程中的流固耦合效應[J].力學學報,2000,32(2):129-140. Lu Chihua,He Yousheng.Coupled analysis of nonlinear interaction between fluid and structure during impact[J].Journal of Theoretical and Applied Mechanics,2000,32(2):129-140.(in Chinese)

[4]Sun H,Faltinsen O M.Water impact of horizontal circular cylinders and cylindrical shells[J].Applied Ocean Research, 2006,28(5):299-311.

[5]Korobkin A,Gueret R,Malenica?.Hydroelastic coupling of beam finite elementmodelwith Wagner theory ofwater impact [J].Journal of Fluids and Structures,2006,22(4):493-504.

[6]Wagner H.über Sto?-und Gleitvorg?ngean der Oberflüche von Flüssigkeiten[J].ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik,1932,12(4):193-215.

[7]Peseux B,Gornet L,Donguy B.Hydrodynamic impact:Numerical and experimental investigations[J].Journal of Fluids and Structures,2005,21(3):277-303.

[8]Stenius I,Rosn A,Kuttenkeuler J.Explicit FE-modeling of fluid-structure interaction in hull-water impacts[J].International Shipbuilding Progress,2006,53:1031-121.

[9]Stenius I,Rosén A,Kuttenkeuler J.Explicit FE-modelling of hydroelasticity in panel-water impacts[J].International Shipbuilding Progress,2007,54(2):111-127.

[10]Aquelet N,Souli M,Olovsson L.Euler-Lagrange coupling with damping effects:Application to slamming problems[J]. Computer Methods in Applied Mechanics and Engineering,2006,195(1):110-132.

[11]Luo H B,Wang S,Guedes Soares C.Numerical prediction of slamming loads on rigid wedge for water entry problem by an explicit finite elementmethod[M].Advances in Marine Structures,GuedesSoares C,and FrickeW(Eds),Taylor&Francis, UK,2011.

[12]Luo H,Wang H,Soares C G.Comparative study of hydroelastic impact for one free-drop wedge with stiffened panels by experimental and explicit finite elementmethods[C]//ASME 2011 30th International Conference on Ocean,Offshore and Arctic Engineering.American Society ofMechanical Engineers,2011:119-127.

[13]駱寒冰,吳景健,王 珊,徐慧.基于顯式有限元方法的二維楔形剛體入水砰擊載荷并行計算預報[J].船舶力學, 2012,16(8):907-914. Luo Hanbing,Wu Jingjian,Wang Shan,XüHui.Parallel computing simulation ofwater entry of a 2D rigid wedge using an explicit finite elementmethod[J].Journal of Ship Mechanics,2012,16(8):907-914.(in Chinese)

[14]駱寒冰,劉 鑫,董德龍,趙中榜,季紅葉,林維學.鋁制加筋板楔形體入水砰擊模型試驗研究[J].水動力學研究與進展,2014,29(4):460-468. Luo Hanbing,Liu Xin,Dong Delong,Zhao Zhongbang,Ji Hongye,Lin Weixue.Experimental investigation ofwater impact of one free-drop wedge with stiffened aluminum panels[J].Journal of Hydrodynamics,2014,29(4):460-468.(in Chinese)

[15]李裕春,時黨勇,趙 遠.ANSYS 11.0、LS-DYNA基礎理論與工程實踐[M].中國水利水電出版社,2008. Li Yuchun,Shi Dangyong,Zhao Yuan.ANSYS 11.0,LS-DYNA basic theory and engineering practice[M].China Water-Power Press,2008.(in Chinese)

[16]LS-DYNA,Keyword User’s Manual,Livermore Software Technology Corporation[M].Version 971/Rev5,May 2001.

Numerical investigation of hydroelastic impact on one free-drop wedgew ith alum inum stiffened panels considering roll angles using ALE algorithm

LUO Han-binga,b,DONG De-longb,JIHong-yeb,ZHAO Zhong-bangb,GAO Yuanb
(a.State Key Laboratory of Hydraulic Engineering Simulation and Safety;b.Departmentof Naval Architecture and Ocean Engineering,Tianjin University,Tianjin 300072,China)

The explicit finite elementmethod using the Arbitrary Lagrangian-Eulerian(ALE)algorithm is applied to investigate the hydroelastic impact of one free-drop wedge with aluminum stiffened panels.The dead-rise angle is 20 degrees.The stiffened panel on each side of the wedge ismade up of 3 longitudinal stiffeners and 2 transverse frameswith different sizes of cross sections.The acceleration motion,stress responses on longitudinal stiffeners and transversewebs of the stiffened panels are obtained.The numerical results are validated with the new experimental results from drop testswith and without roll angles.This solution shows the capability of ALE algorithms and its potential to predict the hydroelastic slamming of local structures on ships.

ALE algorithm;aluminum stiffened panels;hydroelastic impact;roll angle;numerical simulation

U661

：Adoi:10.3969/j.issn.1007-7294.2017.02.007

2016-05-28

國家自然科學基金資助項目（51109159）

駱寒冰（1976-），男，副教授，E-mail:luohanbing@tju.edu.cn；董德龍（1988-），男，工程師。

1007-7294（2017）02-0175-09