楔形體入波浪水面數值模擬

金禹彤，陳吉昌，盧昱錦，肖天航，童明波

南京航空航天大學 航空學院 飛行器先進設計技術國防重點學科實驗室，南京 210016

隨著海洋資源的不斷開發和建設，波浪因其具有較高的能量密度及對飛行器水上迫降、起降和滑行性能有重要影響而備受關注[1]。飛機著水是指飛行器為完成巡邏、反潛和救援等任務時，在正常操作程序中降落于水面的航態。其降落著水過程中，機體承受著較為嚴重的水動沖擊載荷[2]。即使是相對溫和的海況也會對飛行器造成臨界載荷和無法控制的運動。在完成波浪條件降落著水任務的同時，保證機體內部人員的生命安全，對機體結構強度提出了更高要求。因此，進一步掌握飛行器波浪水面降落性能及水體對機身的載荷作用，對飛行器入波浪水面的技術研究和工程應用具有重要意義。

使用物理水池模型試驗對結構物水動力性能進行研究和預報始于20世紀中期，該方法為驗證理論和數值預報方法的有效性提供了檢驗標準，如荷蘭的瓦格寧水池、美國著名的船舶耐波性試驗水池——泰勒水池(David Taylor Model Basin, DTMB)[3-6]和挪威海洋工程技術研究院(MARINTEK)的海洋工程水池[7]等。但物理水池的建造和維護需要大量的資金、設備、場地和人員投入；而且，物理水池模型試驗還受水池大小的限制，模擬試驗結果存在比尺效應，試驗結果還會受到試驗探測儀器精密度、探測方法等諸多因素的影響。考慮物理水池試驗的局限性和經濟性，隨著數值仿真技術的成熟，構造數值水池已經成為研究物體入水沖擊過程的重要方法。

目前構造數值波浪水池的方法主要有源項造波法、推搖板造波法和速度入口造波法。Brorsen和Larsen[8]使用源項造波法完成對二維任意類型波浪的模擬。李勝忠[9]同樣采用源項造波法，完成對線性規則波和Stokes波等波浪的數值模擬。但對于復雜的黏性流動問題，源項表達式不易推導，不能快速有效地完成數值波浪水池的構造。儲慧林[10]、王文華[11]和王平[12]等基于CFD數值方法，結合推板造波和海綿阻尼消波原理實現數值造消波，分別研究魚雷、圓柱和楔形體入規則波波面過程，但該數值造波方法網格量較大，產生過高的計算成本。Longuet和Cokelet[13]最先使用速度入口造波法生成規則波，且結果與Stokes理論解吻合較好。Boo[14]通過對線性規則波和不規則波仿真，驗證了速度入口造波法的有效性，進一步研究線性和非線性不規則波在垂直圓柱上的作用力。從基本物理模型的實用性、計算精度和計算成本等方面考慮，速度入口造波法更為合理可靠，且已在結構物波浪水面滑行問題上[15-16]得到廣泛應用。

陳宇翔等[17-18]采用有限體積法模擬圓柱出/入水過程，用動態鋪層法處理模型與水面的相對運動，對出/入水過程的水面變化捕捉較好。但其弊端是只考慮圓柱的一自由度運動，而忽略圓柱的偏轉運動。Abraham等[19]運用動態網格重構法數值模擬球體入水過程，對球形體附近網格加密，通過網格的變形和重構實現對球形體下落過程的模擬，但該動網格方法產生過高的計算成本。Wick等[20]采用有限體積法求解非定常Navier-Stokes方程，結合流體體積模型(Volume of Fluid, VOF)液面捕捉技術和動網格方法，對某型無人機(Unmanned Air Vehicle, UAV)不同高度的入水沖擊過程進行了數值模擬。Wick采用的動網格方法是將計算域劃分成兩部分，一部分是包裹UAV的球形體(子域)，該部分網格采用六面體結構，并隨著UAV的運動而運動，且子域內的網格不發生變形；另一部分是球形體以外的其他計算域網格(外子域)，該部分網格采用四面體結構，通過網格重構模擬UAV濺落過程。Guo等[21]采用同樣的數值求解方法和動網格技術，研究俯仰角對某型運輸機水上迫降特性的影響。Carrica等[6]采用的動網格方法產生的計算成本過大，且計算精度不高。Qu等[22]運用整體運動網格法處理模型與水體的相互作用，對結構物水上迫降性能進行數值模擬。該動網格方法中，整個計算域(包含網格單元和邊界條件)隨著模型的運動而運動，而無需劃分子域和外子域。結果表明，整體運動網格法可以精確捕捉自由液面，有效避免計算量大、網格質量差等問題。盧昱錦等[23]采用有限體積法和整體動網格法對高速平板著水問題進行數值仿真，研究不同俯仰角對空氣泡現象的影響規律。

基于此，本文以FLUENT軟件為數值模擬求解器，采用結合整體動網格法的速度入口造波方法及VOF液面捕捉技術，在構造線性規則波和線性不規則波數值水池的基礎上，耦合三自由度模型，數值模擬楔形體入波浪水面過程，分析規則波波形位置對楔形體受力特性和運動姿態變化的影響，并初步探索二維楔形體入不規則波及三維楔形體入規則波情況，為飛行器波浪水面起降和滑行問題的研究提供參考和技術支持。

1 數值計算方法

1.1 流體方程及求解

流體控制方程為非定常不可壓縮流動雷諾平均Navier-Stokes(URANS)方程和剪切應力輸運(SST)k-ω湍流模型，水氣交界面采用VOF模型進行捕捉。求解器采用壓力基求解器、非穩態時間格式，壓力-速度項采用PISO(Pressure Implicit with Spltting of Operators)算法進行迭代計算，采用有限體積法離散控制方程：壓力差值格式采用PRESEO格式，其余項均選用二階迎風格式。

1.2 速度入口造波與源項消波

速度入口造波法是給定入口邊界處的波浪位置和水質點速度，隨著時間步的迭代，波浪在入口邊界處生成，并逐漸傳播到計算域內。對于線性規則波，即為常深度二維小振幅推進波，其波面高度η(x,t)相對于坐標x或相對于時間t作周期性的變化，并且波形以一定速度c沿x向傳播。波浪參數及波形位置說明如圖1所示。

圖1 波浪參數及波形位置

設常深度二維小振幅推進波的波面方程η(x,t)為

(1)

式中：H為波高；k為波數；ω為圓頻率。

波浪中水質點(x,y)的速度分量u、v的表達式為

(2)

式中：φ為速度勢函數；T為周期；d為水深。

對于線性不規則波，采用Longuet-Higgins提出的線性疊加海浪模型，將隨機海浪這種平穩隨機過程描述為由無限多不同周期、振幅和不同隨機初相位，并且在xoy平面上與x軸成θ度方向角的斜向余弦波疊加而成的波列，單向線性不規則波中θ=0°。單向線性不規則波波面方程及波浪中水質點(x,y)的速度分量Vx、Vy的表達式為

(3)

式中：Hi為單個組成波的波高；Ti為單個組成波的周期；ki為單個組成波的波數；ωi為單個組成波的圓頻率；ai為單個組成波的振幅；εi為單個組成波的初相位。

阻尼消波法采用動量源項消波方式，即在動量方程中增加一個動量衰減的源項，源項表達式[24]為

(4)

式中：x≠x0表示計算域在造波區；x=x0表示在消波區。

1.3 整體動網格法

整體動網格法是指包括網格單元和邊界在內的整個計算域與模型有相同的運動規律，即隨著模型的運動而運動[22-23]。

整體動網格方法不需考慮彈性變形和網格重構，提高計算效率和精準度的同時，可以有效地處理復雜幾何構型、邊界條件和水氣液面捕捉。在整體動網格法中，采用體積分數邊界條件以保證水平面保持在初始高度。該邊界條件是基于歐拉坐標系網格點的高度設定，空氣體積分數αa=0.5表示該網格位于水氣交界面，αa=1.0表示該網格位于水面以上，αa=0表示該網格位于水面以下，體積分數分布如圖2所示。

圖2 空氣體積分數

2 算例驗證

2.1 數值造波及波形驗證

本節數值驗證速度入口造波法和源項消波法的有效性和準確性。波浪參數[25]為：波高H=0.08 m，水深d=0.5 m，周期T=1.25 s，波長L=2.181 m。

由于自由液面區域附近的物理量梯度變化大，則將自由液面加密區域延伸到0.3 m(± 0.15 m上下平均水面)；此外，為了減小由于垂向和縱向網格尺寸的突變而對數值擴散和阻尼消波的影響，實現網格的平滑過渡，取計算域網格的偏差增長率為1.2。網格參數見表1所示，計算域及邊界條件設置見圖3。用相同的波浪參數、計算域、網格尺寸和加密方式生成并驗證不規則波。

計算輸出規則波波形在x=-6.5，-3，0，3，11 m處的數值解，并與理論解析解對比，如圖4所示。數值解波形向上偏移1%，主要考慮湍流模型的影響；且在x=11 m處數值解波形被吸收，說明阻尼消波法有效，不規則波數值解和理論解析解基本吻合(見圖5)。

圖6給出計算域內規則波和不規則波波形和流場分布。綜上，速度入口造波法和阻尼消波法有效且準確。

表1 網格參數

圖3 波浪水面計算域及邊界條件

圖4 規則波波形對比

圖5 不規則波波形對比

圖6 波形及流場分布

2.2 楔形體靜水面入水

數值模擬楔形體在靜水面的入水過程，并與試驗結果[26]和理論解[27-29]對比，驗證數值方法及整體動網格法模擬物體入水過程的準確性和可行性。楔形體模型[26]的橫截面如圖7所示。

采用結構化網格，整體計算域長和寬均為5 m，其最小網格尺寸在楔形體處為0.005 m，最大網格尺寸在邊界處為0.06 m，計算域網格結構及計算模型初始時刻如圖8所示。初始時刻，楔形體位于靜止水面上方且頂點與水面重合，上方為空氣域(αa=1)，下方為水域(αa=0)。初始速度為-6.15 m/s，楔形體入水按指定規律運動(圖9)。

圖7 楔形體橫截面[26]

圖8 網格計算域及初始位置

圖10為楔形體入水過程中不同時刻自由液面的變形情況，計算域隨楔形體移動，但水面始終保持在初始位置；由于楔形體對水的擠壓作用，楔形體邊緣形成射流。圖11給出楔形體入水所受垂向力的數值解與試驗結果[26]和理論解[27-29]的對比情況。在0.008 s時數值解開始偏離試驗結果，數值解減少10%，在0.011 s時數值解開始偏離理論值，主要是考慮三維效應的影響。但總體來說，數值解與試驗結果和理論解變化趨勢始終相同，誤差較小，在可接受范圍內。

由于目前并沒有對楔形體入波浪水面的有效試驗數據與數值計算參考，在完成楔形體入靜水面驗證的基礎上，運用相同的數值方法研究楔形體入波浪水面情況。

圖9 楔形體入水速度變化

圖10 楔形體入水不同時刻水面變化

圖11 楔形體垂向力隨時間變化

3 計算結果與分析

3.1 二維楔形體入波浪水面數值模擬

本節研究楔形體自由(三自由度)入規則波和不規則波的波浪參數[30]為：波高H=0.12 m，波長L=2.5 m，水深d=2.4 m，周期T=1.25 s(邊界條件如圖3所示)。模型質量為241 kg[28]，通過數值計算得到模型繞重心質點的z向轉動慣量為 9.448 kg·m2。應用整體動網格法并耦合三自由度模型對楔形體在線性規則波的波峰、波谷、平衡位置-上行速度和下行速度4個位置(圖1)、平靜水面和線性不規則波入水情況進行數值模擬。由于不規則波沒有明顯的周期性，所以當不規則波充滿整個水域之后，使模型分別在t=5.0，6.0，7.1，9.0，9.5 s時間點入水，波面波形和著水位置如圖12所示。保證楔形體初始速度均為-6.15 m/s的同時，分別對楔形體的垂向力、垂向速度、橫向速度、橫向位移和滾轉角進行對比分析。

圖12 楔形體入不規則波波面著水位置

楔形體在規則波不同位置入水的對比結果如圖13所示。可知，模型在規則波各位置所受垂向力有先增大后減小的趨勢，峰值均出現在入水初期，且平緩值相近。垂向速度變化趨勢相同。模型在不同位置入橫向波對橫向位移的影響很小，其數值變化均少于0.01。在平衡位置-上行和下行速度位置入水過程中，模型的滾轉角和橫向位移變化明顯，其數值變化約為波峰波谷位置處的8倍和10倍。分析原因：平衡位置處波浪內流速度變化較快，模型所受波浪力沖量較大；楔形體斜邊與波浪的相對作用力也參與到楔形體的入水過程，影響了模型形態的改變。圖14和圖15給出了楔形體在波峰和平衡位置-上行速度處入水過程中不同時刻的水面變化，可以看出，在楔形體斜邊均產生射流現象，且在平衡位置-上行速度處入水過程產生的射流顯示出明顯的不對稱，該現象進一步驗證了上述原因。

楔形體在不同時間點T入不規則波的對比結果如圖16所示。可知，模型在不同時間點入水所受垂向力變化趨勢相同，峰值出現于入水初期，平緩值相近。垂向速度變化趨勢相同。模型在不同位置入橫向波對橫向位移的影響很小，其數值變化均少于0.001。T=5.0，9.5 s工況，楔形體位于類似平衡位置-上行和下行速度位置處，模型的滾轉角和橫向速度變化較為明顯，其數值變化均為T=6.0，9.0 s工況的4倍左右，形成原因同楔形體入規則波情況相同。在T=6.0 s工況中，模型的滾轉角出現由正向負、橫向速度和橫向位移出現由負向正的反向變化，主要考慮楔形體兩側斜邊均與波浪產生相對作用力的影響。圖17給出楔形體在T=7.1 s時間點入不規則波波浪水面過程中不同時刻T的水面變化，可見在楔形體斜邊產生射流現象。

圖13 楔形體入規則波波浪水面的載荷及運動響應結果

圖14 楔形體在波峰處入波浪水面變化

圖15 楔形體在平衡位置-上行速度處入波浪水面變化

圖16 楔形體入不規則波波浪水面的載荷及運動響應結果

圖18給出楔形體在平衡位置-下行速度處和在T=7.1 s入不規則波后0.06 s時的計算域變化。可以看出，包含邊界在內的整個計算域均隨模型的運動而運動，波形保持完整的同時，液面捕捉良好且很好地耦合了三自由度模型。不規則波情況下的計算域運動策略與入規則波相同，可以說整體動網格法很好地適用于速度入口造波法和物體著水運動。

圖17 楔形體在T=7.1 s入不規則波波浪水面不同時刻水面變化

3.2 三維楔形體入規則波波浪水面數值模擬

本節研究三維楔形體自由(六自由度)入規則波情況。波浪水面參數同3.1節。楔形體截面與圖7相同，展長為1 m，模型如圖19所示。楔形體質量采用中國特種飛行器研究所進行水箱試驗的數據，即406 kg。通過數值計算得到模型繞重心質點的轉動慣量分別為Iox=34.303 kg·m2，Ioy=4.699 kg·m2，Ioz=38.063 kg·m2。二維模型和三維模型繞重心質點轉動慣量的不同主要考慮模型自身屬性和三維效應的影響。由3.1節可知，楔形體在規則波波峰和波谷位置處入水數值變化趨勢相似，在平衡位置-上行和下行速度處入水數值變化趨勢相似。在不影響研究目的的同時進行簡化運算，本節選取規則波波峰和平衡位置-上行速度位置，研究三維楔形體入規則波情況。

圖19 三維楔形體模型

計算及對比結果如圖20所示，可知，三維楔形體在規則波波峰和平衡位置-上行速度所受垂向力均出現先增大后減小的趨勢，峰值出現在入水初期，且平緩值相似。垂向速度變化趨勢相同。模型在不同位置入橫向波過程對俯仰角、偏航角和橫向位移的影響很小，其數值變化均少于0.01；而對滾轉角影響較大，在平衡位置-上行速度處的滾轉角數值變化約為波峰的4倍。其形成原因與二維楔形體入規則波相同。

圖21給出楔形體平衡位置-上行速度處入水過程中不同時刻的水面變化。三維楔形體入水后計算域變化如圖22所示。

圖20 三維楔形體入規則波波浪水面的載荷及運動響應結果

圖21 三維楔形體在平衡位置-上行速度處入波浪水面不同時刻水面變化

圖22 三維楔形體在平衡位置-上行速度處入水計算域變化

4 結 論

1) 采用結合整體動網格法的速度入口造波法和VOF液面捕捉技術，成功構造數值波浪水池，線性規則波和線性不規則波數值解與理論解析解基本吻合，偏差為1%。

2) 數值模擬二維楔形體分別在規則波波峰、波谷、平衡位置-上行速度和下行速度4個位置處自由入水，其垂向力峰值均出現在入水初期。楔形體在平衡位置-上行和下行速度處入水對模型的滾轉角、橫向速度和橫向位移變化影響較大，其數值變化約為波峰波谷位置處的8倍、4倍和10倍。主要考慮該處波浪內流速度較快，模型所受波浪力沖量較大，且楔形體斜邊與波浪的相互作用力也參與到模型入水過程。

3) 當不規則波充滿水域后，數值模擬楔形體在5個隨機時間點入不規則波情況。在波浪內流速度較大處，模型的滾轉角、橫向位移和橫向速度均出現明顯變化。在t=6.0 s工況，數值出現由正到負或由負到正的反向變化，主要考慮楔形體兩側斜邊均與波浪產生相互作用力。

4) 數值模擬三維楔形體在規則波波峰和平衡位置-上行速度自由(六自由度)入水。三維楔形體在入水過程中所受垂向力出現先增大后減小的趨勢，垂向力峰值均出現在入水初期，且平緩值相似。垂向速度和俯仰角變化趨勢相同。模型在不同位置入橫向波過程對俯仰角、偏航角和橫向位移的影響很小，其數值變化均少于0.01；而對滾轉角影響較大，且在平衡位置-上行速度處的滾轉角數值變化約為波峰的4倍。其形成原因與二維楔形體入規則波情況相同。

5) 整體動網格法結合速度入口造波法、VOF液面捕捉技術并耦合三/六自由度模型數值模擬楔形體入波浪水面情況較好，且波面保持完整。