二維楔形體入水問題的數值和實驗研究

陳光茂，鄭小波，毋曉妮，張懷新，李 曄

(1.上海交通大學 海洋工程國家重點實驗室，上海 200240；2.上海交通大學 船舶海洋與建筑工程學院，上海 200240)

0 引 言

當船舶遭遇惡劣海況時，船體和波浪會有很強的相互作用，因此導致船首與水面產生瞬時劇烈的砰擊。這種結構物入水砰擊過程是一個涉及到氣液固三相的強流固耦合問題，入水瞬間產生的砰擊力很大，極易危害船舶安全，對船舶結構甚至是船上人員造成極大的危險。在救生船下水、導彈入水等過程中同樣涉及到類似的砰擊過程，因此研究此過程中結構物和自由液面的相互作用意義重大。

入水砰擊問題的研究開始于von Karman[1]和Wagner[2]的研究，von Karman[1]基于二維剛性楔形體模型和附加質量理論，首次給出了入水問題的理論解法，Wagner[2]基于線性勢流理論并考慮自由面抬升對入水問題解法提出了有效的改進；之后，Dobrovol’ skaya[3]基于Wagner[2]的理論首次給出了二維楔形體入水的相似解；Greenhow和Lin[4]通過一系列實驗拍攝了不同速度和形狀的圓球和楔形體入水過程中水面的變化情況，但是他們的結果并不清晰；Zhao和Faltinsen[5]運用完全非線性的勢流理論求解了楔形體入水過程中速度與壓力的變化并通過實驗驗證了此方法在壓力分布計算上有一定的準確性。隨著計算流體力學（CFD）理論和計算機能力的發展，基于Navier-Stokes 方程的數值方法也逐漸被運用于楔形體入水過程的計算中。Kleefsman等[6]采用有限體積法（FVM）并使用流體體積比函數（VOF）方法進行自由面捕捉，并將垂向速度和砰擊力與實驗結果進行對比，證明此方法有很好的準確性；王加夏等[7]將VOF法和有限元法耦合，研究三維楔形體結構入水砰擊時結構動態響應特性及板厚對結構變形效應的影響規律；Zhang等[8]使用結合隱式界面追蹤法（level set）的浸沒邊界法（Immersed Boundary Method）來研究楔形體入水問題，但所得到的楔形體上的壓力分布結果并不理想。無網格方法由于其便于處理自由表面和大變形的特點，在楔形體入水問題中也得到了很好的應用。Koshizuka等[9]使用移動粒子半隱式方法（MPS）模擬了破碎波與浮體的相互作用，并由此驗證了無網格方法在處理大變形問題上的能力；陳翔等[10]將此方法應用在楔形體入水問題中；Oger等[11]使用經典的弱可壓縮SPH方法模擬了楔形體入水問題并與Zhao和Faltinsen[5]的實驗進行了對比，結果吻合良好，驗證了SPH方法在解決入水砰擊問題上的適用性；龔凱等[12]使用改進的固壁處理方法得到了更準確的計算結果。

SPH方法是一種純拉格朗日的無網格方法，它首先被應用于天體物理的模擬中[13-14]，后來被用于多種流體力學的研究中，如造波[15]、多相流[16]和流固耦合[17]等。該方法將整個流場離散為一系列相互作用的粒子，并通過插值將流場中的物理量在離散的粒子上進行表達。SPH方法的主要優勢在于無需定義空間網格，從而避免了自由表面大變形問題中網格的扭曲變形導致的計算精度下降和計算時間增加。隨著高性能計算的發展，SPH方法計算效率低的缺點也逐漸被克服，它在處理各種問題，特別是傳統網格方法無法很好解決的大變形問題上有很大的潛力。

本文主要工作為使用SPH方法研究二維楔形體入水砰擊問題，分析楔形體斜升角和入水速度對入水作用力的影響。首先介紹SPH方法的基本理論，隨后從2個方面驗證此方法的準確性：一是驗證本方法計算入水過程中楔形體運動、受力和壓力分布的準確性；二是進行楔形體入水實驗，驗證本方法計算自由表面變形的準確性。最后利用SPH數值模擬得到不同斜升角和入水速度的楔形體受力情況，分析這2種因素的影響。

1 SPH 方法介紹

1.1 基本理論

SPH方法的理論基礎是積分插值理論[18]?？紤]粘性作用，可以得到如下質量守恒方程和動量守恒方程：

其中 Πij是用于減少不穩定性的人工粘性項，其形式如下：

經典的SPH方法認為流體是弱可壓縮的。因此可以使用狀態方程來計算流體壓力而不需要引入泊松方程，這使得計算更有效率并且適合于并行計算。Monaghan[19]提出了如下的狀態方程：

1.2 關鍵設置

1.2.1 光滑函數

無網格方法的核心問題之一是如何基于以任意方式分散的一組節點有效地執行函數逼近而不使用提供節點連接性的預定義網格。SPH方法使用光滑函數來解決這個問題。光滑函數的選擇至關重要，因為它決定了插值模式并且定義了粒子的影響區域。選擇合適的光滑函數，可以提高精度和穩定性。

光滑函數是無量綱的粒子間距q的函數，q=r/h，其中r是任意2個粒子的間距，h是光滑長度。本文使用的五次光滑函數[20]如下：

其 中 αD在 二 維 問 題 中 為 7 /4πh2， 在 三 維 問 題 中 為21/16πh3。

1.2.2 邊界設置

SPH方法中邊界的設置并不像網格方法中那樣簡單，邊界的設置影響著計算的穩定性和準確性。原因在于邊界附近粒子的支持域是截斷的，缺少了一部分支持域的粒子參數計算會出現很大誤差，邊界條件的設置需要盡可能減少這種誤差。本文使用的是動態邊界條件（DBC）[21]，邊界固體粒子滿足與流體粒子相同的方程，但是它們不會在流體的作用力下運動，而是保持固定或者根據預先設置的運動方程運動。當流體粒子靠近壁面時，壁面粒子的密度會增大，壓力也隨之增大，對流體粒子形成的排斥力可以阻止流體穿過邊界。

1.2.3 時間積分格式

多種時間積分方法都可用于SPH方程的積分求解。本文使用的Verlet方法[22]計算耗時少且結果的準確性較好。其計算方法如下：

2 楔形體入水過程中的運動和受力分析

楔形體入水問題中最受關注的是楔形體的運動狀態和受力情況。本節對二維楔形體自由入水過程進行模擬，分析此過程中楔形體運動速度的變化，砰擊力的變化以及不同瞬間的壓力分布，驗證SPH方法計算楔形體入水砰擊時運動狀態和受力情況的準確性。

2.1 數值模型

針對二維楔形體入水問題建立數值模型，其中水池長2.4 m，高1.3 m，以此保證入水砰擊的過程中壁面反射不會對楔形體受力產生干擾。楔形體形狀為等腰三角形，頂邊長度為0.5 m，斜升角（即斜邊與水平方向的夾角）為30°。在楔形體入水前的瞬間，它的垂向速度為6.15 m/s，隨后自由入水，如圖1所示。

圖1 楔形體入水問題圖示Fig.1 Sketch of wedge water entry

計算時流體中的聲速c0為67.7 m/s，是入水速度的10倍以上，以此確保流體壓縮率小于1%。計算中使用的二維楔形體為實心剛體，重量為241 kg，在二維平面內自由運動。粒子初始間距設置為0.001 2 m，粒子總數近300萬，計算的總下落時間為0.024 s。

2.2 結果分析

首先分析自由入水后的垂向速度變化，本次計算結果與Zhao和Faltinsen[5]的實驗結果進行對比（見圖2）。首先可以看出兩者在變化趨勢上有高度的一致性，下落加速度在計算時間內呈現出相同的先增大后減小的特征。由于模型實驗存在三維效應，數值和實驗結果會有一定程度的差異。在入水后的0.003 s內，本文結果與實驗非常吻合，其后SPH方法速度衰減開始快于實驗，在0.015 s時SPH方法速度比實驗結果小4.1%，最終速度的差異為6.7%。

圖2 入水垂向速度隨時間變化Fig.2 Comparison of velocity time history

在砰擊入水的過程中，流體對楔形體的垂向瞬時砰擊力變化是研究中關注的重點問題之一。本文將SPH方法的砰擊力計算結果與Zhao和Faltinsen[5]的實驗結果和解析解對比，結果如圖3所示。可以看出，三者的變化趨勢一致且結果十分接近。入水后的0.06 s內，解析法與實驗結果吻合度良好，SPH方法結果略微偏大；0.06 s之后解析法和SPH結果都大于實驗值，但SPH方法結果更接近實驗：在砰擊力峰值處，SPH方法的結果相比實驗大10%，而解析解相比實驗大30%；入水后0.02 s，SPH方法和實驗相差達到最大值30%，此時解析解和實驗相差為50%。說明SPH方法能更準確地計算楔形體入水砰擊過程中的垂向砰擊力變化。同時可以看出，砰擊力的計算結果和入水速度的變化趨勢相符，SPH計算的砰擊力大于實驗值，因此入水減速更快。

圖3 入水垂向砰擊力隨時間變化Fig.3 Vertical slamming force time history

通過以上分析可知，SPH方法在砰擊力的計算上有很好的準確性。接下來通過研究不同瞬間楔形體表面壓力分布的計算結果，更深入地驗證SPH方法在分析楔形體所受作用力上的準確性。本文選取了2個特定的時間點（t=0.004 4 s 和t=0.015 8 s）來分析沿著楔形體表面的壓力分布，將其與Zhao和Faltinsen[5]的實驗結果和解析解以及Oger等[11]的SPH結果進行對比，結果如圖4所示。由于SPH方法存在壓力不穩定的問題，特別是在邊界上問題更嚴重，本文在SPH方法中使用了delta-SPH[23]來減少這種不穩定的影響。delta-SPH即在連續性方程中加入擴散項，使其變為：

其中 δΦ的值取為0.1。由于SPH方法在邊界上的壓力計算不準確，本文參考Oger等[11]的做法，使用流體粒子壓力的平均值代替邊界上的壓力值，具體做法為以楔形體表面邊界點為圓心，取半徑為6 h的范圍內所有流體壓力的平均值作為此邊界點壓力值。圖中橫坐標為楔形體與自由表面接觸點的無量綱縱坐標，縱坐標為無量綱壓力。其中P為壓力，V為楔形體運動速度，豎直向下為正方向，Z為楔形體與自由表面接觸點的縱坐標，原點在靜水面，豎直向上為正方向。

圖4左圖為2個時刻壓力分布對比結果，右圖為對應時刻SPH方法計算的壓力分布圖。從右圖可以看出，delta-SPH的使用有效保證了SPH結果中壓力分布的穩定性，2個時刻的液體內部壓力分布都有比較合理的梯度變化。從左圖可以看出，對于2個時刻的壓力分布，3種非實驗方法得到的結果比實驗大，其中同樣有實驗三維效應的原因。圖4（a）中本文的SPH方法很好地捕捉到了壓力峰值的位置，并且與解析解的結果整體一致；圖4（b）中本文所用SPH方法的壓力分布結果在壓力峰值點前都與解析解十分吻合，并且壓力峰值點也能準確捕捉。圖4（b）的結果在楔形體上緣迅速下降但是并不為0，除了SPH方法壓力計算存在的問題外，一部分原因是提取壓力時采用的平均算法導致邊緣壓力變大。對比Oger等[11]的SPH結果，本次研究的結果在t=0.004 4 s時整體更接近解析解，較好地捕獲到了壓力分布的峰值；在t=0.015 8 s時大部分區域都更接近解析解，只有靠近邊緣的區域由于平均算法的問題誤差相對較大。

圖4 壓力分布對比，入水時間分別為 0.004 4 s和 (b) 0.015 8 sFig.4 Comparisons of pressure distribution at three instants

以上分析證明SPH方法在處理楔形體入水砰擊問題時，能夠準確地計算速度變化、砰擊力變化和不同時刻壓力分布，說明SPH方法在這一問題上有很好的適用性。

3 自由表面變形驗證實驗

楔形體砰擊入水是楔形體和液體表面互相作用的過程，其中既有楔形體的受力和運動，也有自由表面由于楔形體砰擊和擠壓產生的變形和射流，同時自由表面的變化也將在入水的后續過程中影響楔形體的受力和運動，如氣穴的產生和封閉等將對楔形體產生二次作用力。過去的大多數研究都只關注了楔形體受力而沒有考慮自由表面變形，為從這一角度驗證SPH方法在楔形體入水問題上的適用性，進行相應的實驗研究。

3.1 實驗設置

本次實驗使用自行設計和搭建的實驗裝置完成，主要包括以下設備：上表面長0.3 m，寬0.1 m，斜升角 30°，總重量為 4 kg 的楔形體；長 1 m，寬 0.35 m，高 1 m 的有機玻璃缸；型號為 Vision Miro eX4，幀數1 000 fps，像素 800×600 的高速相機；波長為 532 nm的Nd:YAG激光器。實驗設置如圖5所示。

本次實驗采用的是自由落體的入水方式。為保證實驗的二維性，將水缸和楔形體的寬度設置地盡量接近，同時可以減少正面濺射的水花，使拍攝結果更加清晰。使用電磁鐵作為釋放裝置，用2根光滑導軌限制楔形體釋放之后的自由度，使其在下落過程中不會翻轉和偏移。為防止外界光源和其他物體的干擾，保證實驗拍攝效果，本次實驗在黑暗環境中進行，使用激光器發出的片光作為光源。

圖5 實驗設置圖Fig.5 Experimental setup

實驗設置2種不同的入水高度，分別為0.05 m和0.1 m。實驗開始時，首先將楔形體固定在指定高度，開啟高速相機和激光器，將相機和激光器分別調整至光軸正對楔形體寬度方向和長度方向的中軸線，相機成像清晰。隨后啟動相機錄制，將電磁鐵斷電，釋放楔形體，利用高速相機捕捉楔形體入水過程中的自由表面變形和射流飛濺現象。

3.2 結果對比

利用SPH方法建立與本次實驗相同的容器和楔形體的二維數值模型，計算相同高度下楔形體入水砰擊過程，可以得到數值模擬的自由表面變形結果。下面將SPH方法的結果和實驗結果對比以驗證此方法的準確性。

圖6為入水高度0.05 m時的對比，其中左側為實驗結果，右側為SPH方法的結果，可以看出實驗結果中水和空氣以及水和楔形體的交界面都被很好地捕捉到了，射流的輪廓也清晰可見。圖6（a）入水時間為0.05 s，自由表面由于楔形體的作用向兩側產生變形，并產生高速射流。可以看出SPH方法對此時的自由表面變形有比較準確的計算，與實驗結果的形狀十分吻合，在射流的模擬上相對較弱，其中部分原因是數值方法的粒子尺度較大，同時未考慮表面張力的作用。圖6（b）入水時間為0.136 s，自由表面變形增大，射流也更加明顯，可以看出此時SPH方法不僅對變形有很好的準確性，同時也能夠有效模擬出射流飛濺的效果，射流分離位置和角度吻合良好。此時隨著楔形體的下落，上方的自由表面開始重新向中間靠攏，即將產生閉合，這一時刻楔形體上方的自由表面接近豎直，本方法對此也有準確的計算。

圖6 入水過程自由表面變形對比，入水高度 0.05 mFig.6 Comparisons of free surface deformation，height=0.05 m

圖7為入水高度0.1 m時的對比，相對于0.05 m的入水高度，高度增加之后入水速度增大，因此將拍攝范圍增大以捕捉到相同時長的入水過程。圖7（a）入水時間為0.05 s，可以看出SPH方法對此時的自由表面變形同樣有比較準確的計算，由于入水速度增加，砰擊力也增大，因此自由表面變形更大，擠壓變形的趨勢在大砰擊力的作用下明顯更向豎直方向傾斜。圖7（b）入水時間為0.136 s，SPH方法同樣很好地模擬出了自由表面變形和射流飛濺的現象，也可看出楔形體上方自由表面即將向中間靠攏的傾向。在入水高度為0.05 m時，這一時刻被楔形體擠壓而抬升的自由表面已經有明顯的回落趨勢，而入水高度增大時，由于砰擊力變大，自由表面上液體上升速度增大，因此速度降為0并開始回落的時間更晚。SPH計算結果中同樣可以看到不同高度的這一區別：相對于入水高度0.05 m的情況，入水高度0.1 m時，楔形體附近自由表面靠攏趨勢更不明顯，靜水面以上的自由表面尚未有明顯往下塌陷的趨勢。

圖7 入水過程自由表面變形對比，入水高度 0.1 mFig.7 Comparisons of free surface deformation，height=0.1 m

從以上分析可以看出，2種入水高度下，SPH計算結果和實驗結果都有很好的吻合度。楔形體入水時本方法對自由表面形狀變化的模擬十分接近實驗結果，射流的情況也有一定程度的相似性，說明SPH方法可以很好地模擬入水砰擊時的自由表面變化情況，也從另一個側面表明SPH對楔形體入水砰擊過程的計算是十分準確的。

4 楔形體斜升角和入水速度的影響

4.1 楔形體斜升角的影響

楔形體的斜升角會對入水砰擊力產生影響，其中需要重點關注的是砰擊作用時間和最大砰擊力。本文選取 20°，30°，40°，50°和 60°五種斜升角的楔形體，計算其入水砰擊力隨時間的變化，計算過程中保持每個楔形體頂邊長度為0.1 m，質量為20 kg，入水速度為3 m/s，限制楔形體只能在豎直方向做直線運動。

每個楔形體的計算結果如圖8（a）所示。從圖中可以看出，斜升角越小，入水砰擊力在越短的時間內達到最大值，且下降速度更快?？梢钥闯觯S著斜升角的增大，最大入水砰擊力的減小非常明顯，其中斜升角為 20°時最大砰擊力為 2 909.4 N，60°時最大砰擊力為516.0 N，前者是后者的5.64倍，說明改變斜升角對砰擊力的影響很大。圖8（b）給出了不同斜升角下的無量綱的最大入水砰擊力及其3次擬合曲線，曲線方程為：

圖8 不同斜升角的楔形體入水砰擊力對比Fig.8 Comparison of slamming force between wedges with different rise angles

4.2 入水速度的影響

除了斜升角以外，楔形體入水的速度也是需要考慮的重點因素之一，本文選擇從1 m/s到6 m/s共6種入水速度，分析入水速度變化對砰擊力的影響。計算中采用的楔形體頂邊長度為0.1 m，斜升角為30°，質量為20 kg，限制其只能在豎直方向作直線運動。

每個楔形體的計算結果如圖9（a）所示。從圖中可以看出，入水速度越大，砰擊力在越短的時間內達到最大值，隨著速度的增大，最大入水砰擊力明顯增大，速度從1 m/s增加到6 m/s時，最大砰擊力增大了21.23倍。圖9（b）為6種速度下的最大入水砰擊力及其二次擬合曲線，曲線方程為：

5 總結和展望

圖9 不同速度的楔形體入水砰擊力對比Fig.9 Comparison of slamming force between wedges with different speeds

本文使用SPH方法研究了二維楔形體入水砰擊問題，分析楔形體斜升角和入水速度對入水砰擊力的影響。首先通過分析楔形體入水過程中的速度變化、砰擊力變化和不同時刻壓力分布，驗證SPH方法在處理楔形體入水問題時，能夠準確地計算楔形體受力和運動。隨后進行楔形體入水實驗并成功拍攝到了十分清晰的自由表面變形結果，將上述結果與SPH方法計算結果作對比，證明2種入水高度下，SPH模擬的自由表面變化十分接近實驗結果，產生的射流也有一定程度的相似性，從而驗證本方法計算自由表面變形的準確性。通過以上2個方面，驗證了對于楔形體入水這種瞬時強相互作用過程，SPH方法有很好的適用性和準確性。

在上述研究的基礎上，從最大砰擊力和入水砰擊時間上分析楔形體斜升角和入水速度對砰擊過程的影響。從結果可以看出，改變這2種參數都會明顯改變砰擊過程，從擬合結果可以推斷，最大砰擊力與楔形體斜升角的三次方成正比，與入水速度的二次方成正比。

本文的主要創新點，一是使用合理的實驗設置和拍攝手段，成功拍攝了楔形體入水過程中清晰的自由表面變化情況；二是給出了最大砰擊力與斜升角和入水速度的定量關系。之后的研究將重點繼續這兩方面的工作，第一是實驗覆蓋更多的入水速度和斜升角并拍攝入水后期的自由表面變化，第二是計算更多不同入水速度和斜升角的情況，提高定量關系的準確性。