楔形體入波浪水面數(shù)值模擬

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

南京航空航天大學(xué) 航空學(xué)院 飛行器先進設(shè)計技術(shù)國防重點學(xué)科實驗室，南京 210016

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

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

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

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

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

1 數(shù)值計算方法

1.1 流體方程及求解

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

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

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

圖1 波浪參數(shù)及波形位置

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

(1)

式中：H為波高；k為波數(shù)；ω為圓頻率。

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

(2)

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

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

(3)

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

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

(4)

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

1.3 整體動網(wǎng)格法

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

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

圖2 空氣體積分?jǐn)?shù)

2 算例驗證

2.1 數(shù)值造波及波形驗證

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

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

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

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

表1 網(wǎng)格參數(shù)

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

圖4 規(guī)則波波形對比

圖5 不規(guī)則波波形對比

圖6 波形及流場分布

2.2 楔形體靜水面入水

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

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

圖7 楔形體橫截面[26]

圖8 網(wǎng)格計算域及初始位置

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

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

圖9 楔形體入水速度變化

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

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

3 計算結(jié)果與分析

3.1 二維楔形體入波浪水面數(shù)值模擬

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

圖12 楔形體入不規(guī)則波波面著水位置

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

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

圖13 楔形體入規(guī)則波波浪水面的載荷及運動響應(yīng)結(jié)果

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

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

圖16 楔形體入不規(guī)則波波浪水面的載荷及運動響應(yīng)結(jié)果

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

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

3.2 三維楔形體入規(guī)則波波浪水面數(shù)值模擬

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

圖19 三維楔形體模型

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

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

圖20 三維楔形體入規(guī)則波波浪水面的載荷及運動響應(yīng)結(jié)果

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

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

4 結(jié) 論

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

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

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

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

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