剛性球體三維入水沖擊數(shù)值仿真

孟 龍， 萬(wàn)美慧， 肖京平

(中國(guó)空氣動(dòng)力研究與發(fā)展中心， 四川 綿陽(yáng) 621000)

0 引 言

水上飛機(jī)和宇宙飛船的水上降落、太空艙濺落至海洋中等情況均屬于結(jié)構(gòu)物入水沖擊問(wèn)題，結(jié)構(gòu)物的入水沖擊過(guò)程持續(xù)時(shí)間非常短，瞬時(shí)砰擊壓力大，流體呈現(xiàn)較強(qiáng)的彈性效應(yīng)，并產(chǎn)生振動(dòng)[1]，流固耦合作用復(fù)雜，有些流動(dòng)機(jī)理尚待研究。

早在1897年WORTHINGTON等[2]就開(kāi)始入水問(wèn)題的研究，采用拍攝技術(shù)對(duì)外形為圓形的結(jié)構(gòu)物進(jìn)行入水界面處的拍攝，對(duì)水氣界面處的液滴飛濺和液面閉合情況進(jìn)行細(xì)致的研究。VON KARMAN[3]對(duì)入水問(wèn)題進(jìn)行系統(tǒng)的理論研究，但是并沒(méi)有考慮流體的堆積效應(yīng)。WAGNER[4]基于勢(shì)流理論提出近似平板理論，該理論適用于小斜升角并且考慮結(jié)構(gòu)物入水時(shí)引起的堆積效應(yīng)，奠定了后續(xù)入水問(wèn)題理論研究和數(shù)值模型的基礎(chǔ)，然而該理論仍需要引進(jìn)修正系數(shù)進(jìn)行精確計(jì)算[5]。國(guó)內(nèi)眾多學(xué)者[6-7]也開(kāi)展大量的結(jié)構(gòu)物入水?dāng)?shù)值研究，馬林[8]主要采用光滑粒子流體動(dòng)力學(xué)(Smoothed Particle Hydrodynamics，SPH)方法開(kāi)展結(jié)構(gòu)物入水問(wèn)題的研究，但是其建立的SPH模型對(duì)于三維入水問(wèn)題的使用性較差，對(duì)計(jì)算資源要求較高，計(jì)算時(shí)間較長(zhǎng)。

隨著數(shù)值方法(如有限差分法[9-10]、邊界元法[11-12]、有限元法[13-14]等)的發(fā)展及其在流體流動(dòng)問(wèn)題上的應(yīng)用，對(duì)于結(jié)構(gòu)物入水沖擊載荷的試驗(yàn)研究[15]也迅速發(fā)展起來(lái)：張?jiān)狼嗟萚16]開(kāi)展楔形結(jié)構(gòu)物入水問(wèn)題的研究，主要研究結(jié)構(gòu)物在入水瞬間砰擊時(shí)的運(yùn)動(dòng)響應(yīng)，對(duì)于其他入水砰擊問(wèn)題具有一定的指導(dǎo)作用；ALAOUIA等[17]與張?jiān)狼嗟妊芯肯嗨?，但是兩者研究的結(jié)構(gòu)物外形不同，ALAOUIA等主要研究金字塔形狀結(jié)構(gòu)物入水瞬時(shí)的砰擊問(wèn)題；SHAMS等[18]對(duì)入水結(jié)構(gòu)體在外形不對(duì)稱情況下開(kāi)展入水砰擊研究，并且進(jìn)行入水參數(shù)對(duì)載荷影響分析；CHUANG等[19-24]開(kāi)展平板、楔形體、圓錐體和船模等的入水試驗(yàn)，得出砰擊壓力的回歸經(jīng)驗(yàn)公式；佘文軒等[25]開(kāi)展垂直入水砰擊試驗(yàn)，主要觀察宏觀運(yùn)動(dòng)響應(yīng)和細(xì)節(jié)流場(chǎng)結(jié)構(gòu)；朱貴勛[26]運(yùn)用多相流模型模擬兩相入水砰擊問(wèn)題，但是其研究主要應(yīng)用于二維，對(duì)于三維情況還需要進(jìn)行改進(jìn)。

對(duì)于結(jié)構(gòu)體入水問(wèn)題，大部分研究?jī)H進(jìn)行定性分析，并未給出定量的關(guān)系式。本文對(duì)不同尺度模型的數(shù)值仿真進(jìn)行系統(tǒng)規(guī)劃與研究，基于CFD方法開(kāi)展球體入水沖擊載荷特性研究。在不同入水速度、入水角度等參數(shù)下進(jìn)行結(jié)構(gòu)體入水沖擊研究，分析模型在不同參數(shù)條件下的入水沖擊載荷特性問(wèn)題，進(jìn)一步揭示流固耦合問(wèn)題的復(fù)雜機(jī)理。

1 結(jié)構(gòu)體-流體數(shù)值仿真模型計(jì)算工況及模型驗(yàn)證

1.1 數(shù)值仿真模型

在三維結(jié)構(gòu)物入水沖擊數(shù)值仿真中，采用笛卡爾坐標(biāo)系，x軸正方向水平向右，y軸正方向指向圖1內(nèi)部并與x軸形成右手螺旋定則，z軸正方向垂直向上，坐標(biāo)原點(diǎn)選在平均水線面處，結(jié)構(gòu)體入水沖擊示例如圖1所示，其中：H為剛性球體距離水平面的高度；h為剛性球體底部距離上方水平面的高度；R為剛性球體半徑；r為剛性球體入水某時(shí)刻觸水半徑。剛性球體以一定的初速度和角度入水，入水角度為結(jié)構(gòu)體與水平面之間的銳角夾角，圖1中入水角度為90°。

采用商業(yè)CFD軟件和重疊網(wǎng)格技術(shù)開(kāi)展數(shù)值仿真研究，建立多個(gè)數(shù)值仿真模型，如表1所示，質(zhì)量為2 370～8 000 g，球體半徑為0.06～0.30 m。球體數(shù)值仿真模型簡(jiǎn)化為剛體，具有六自由度運(yùn)動(dòng)特性，可在x、y、z等3個(gè)方向上進(jìn)行平動(dòng)和轉(zhuǎn)動(dòng)。流體域分為上部空氣域和下部水域，在整個(gè)計(jì)算域中有4種邊界條件：速度進(jìn)口邊界、壓力出口邊界、對(duì)稱邊界和壁面邊界。剛性球體在入水砰擊過(guò)程中會(huì)導(dǎo)致自由液面產(chǎn)生較大變化，自由液面的變化同時(shí)也會(huì)對(duì)剛性球體產(chǎn)生重要影響，為了更加精確地得到模型入水時(shí)刻的沖擊載荷特性，在自由液面附近進(jìn)行網(wǎng)格加密，同時(shí)也對(duì)仿真模型進(jìn)行加密處理。

表1 數(shù)值仿真模型

為消除網(wǎng)格密度對(duì)仿真結(jié)果的影響，采用網(wǎng)格密度分別為505萬(wàn)個(gè)、755萬(wàn)個(gè)和913萬(wàn)個(gè)網(wǎng)格開(kāi)展網(wǎng)格無(wú)關(guān)性驗(yàn)證。隨著網(wǎng)格密度的增大，計(jì)算結(jié)果呈減小趨勢(shì)，后兩者網(wǎng)格仿真結(jié)果的誤差為3.78%，但是計(jì)算量顯著增大。綜合考慮，選取網(wǎng)格數(shù)量為755萬(wàn)個(gè)的網(wǎng)格開(kāi)展數(shù)值仿真，采用重疊網(wǎng)格技術(shù)實(shí)現(xiàn)剛體在計(jì)算域中的運(yùn)動(dòng)，仿真網(wǎng)格情況如圖2所示。

圖2 仿真網(wǎng)格情況

以某工況為例，模型質(zhì)量選取4 kg，入水速度為15 m/s，入水角度為90°，仿真獲得的球體入水時(shí)水面處水體積分?jǐn)?shù)變化如圖3所示。

圖3 球體入水時(shí)水面處水體積分?jǐn)?shù)變化

1.2 控制方程和湍流模型

基于雷諾時(shí)均的三維黏性不可壓Navier-Stokes方程的連續(xù)性方程和動(dòng)量方程分別為

(1)

(2)

采用剪應(yīng)力傳輸模型(Shear Stress Transfer，SST)k-ω湍流模型，該模型結(jié)合k-ω模型在近壁區(qū)域計(jì)算的魯棒性和精確性等優(yōu)點(diǎn)與k-ε模型適用于遠(yuǎn)場(chǎng)計(jì)算的優(yōu)點(diǎn)，其輸送方程可寫為

(3)

(4)

1.3 重疊網(wǎng)格技術(shù)

重疊網(wǎng)格是指在計(jì)算域中網(wǎng)格相互重疊在一起進(jìn)行計(jì)算的一種網(wǎng)格方案。激活單元、接收單元和失效單元組成重疊網(wǎng)格技術(shù)的3種基本單元，激活單元與接收單元之間的插值如圖4所示，背景網(wǎng)格與重疊網(wǎng)格共同組成計(jì)算域網(wǎng)格[27]。重疊網(wǎng)格與背景網(wǎng)格重疊邊界處選取相同大小的網(wǎng)格，采用線性插值方式進(jìn)行數(shù)據(jù)傳遞。在數(shù)值仿真過(guò)程中，先進(jìn)行穩(wěn)態(tài)仿真計(jì)算，待運(yùn)行穩(wěn)定后，在穩(wěn)態(tài)計(jì)算結(jié)果基礎(chǔ)之上開(kāi)展非穩(wěn)態(tài)數(shù)值仿真研究，其中仿真時(shí)間步長(zhǎng)選取1×10-6s，采用二階臨時(shí)離散。

圖4 激活單元與接收單元之間的插值

1.4 計(jì)算工況和CFD仿真模型驗(yàn)證

采用建立的數(shù)值仿真模型，考慮入水速度、入水角度等參數(shù)的變化，進(jìn)行不同工況的數(shù)值仿真研究，具體工況如表2所示，其中入水角度為結(jié)構(gòu)體速度方向與平均水線面之間的夾角。根據(jù)設(shè)計(jì)的數(shù)值仿真工況，數(shù)值仿真主要研究?jī)?nèi)容如下：

表2 數(shù)值仿真工況

(1) 研究結(jié)構(gòu)體入水時(shí)不同體積情況下的入水沖擊載荷特性。保持結(jié)構(gòu)體質(zhì)量相同、入水速度相同、入水角度相同，采用5種不同體積的模型，計(jì)算工況為L(zhǎng)C1～LC5。由于研究球體入水，因此在保證質(zhì)量相同但是體積不同時(shí)，采用空心球模型，通過(guò)控制空心球的壁厚來(lái)控制球體體積。

(2) 研究結(jié)構(gòu)體入水時(shí)不同入水速度情況下的入水沖擊載荷特性。采用同一個(gè)M1試驗(yàn)?zāi)Ｐ?，保持結(jié)構(gòu)體入水角度相同，改變結(jié)構(gòu)體入水速度，計(jì)算工況為L(zhǎng)C1、LC6和LC9。

(3) 研究結(jié)構(gòu)體入水時(shí)不同入水角度情況下的入水沖擊載荷特性。采用同一個(gè)M1試驗(yàn)?zāi)Ｐ?，保持結(jié)構(gòu)體入水速度相同，采用多種不同入水角度，計(jì)算工況為L(zhǎng)C10～LC15。

(4) 研究縮尺模型的縮尺效應(yīng)對(duì)結(jié)構(gòu)體入水沖擊載荷的影響。保持結(jié)構(gòu)體入水角度相同，質(zhì)量、頭部直徑和入水速度按照縮尺比進(jìn)行試驗(yàn)，計(jì)算工況為L(zhǎng)C7和LC8。此時(shí)采用的縮尺比為傳統(tǒng)海洋工程使用的縮尺法則，具體推導(dǎo)公式可參照文獻(xiàn)[28]。

采用CFD方法和文獻(xiàn)[29]中的理論計(jì)算方法進(jìn)行球體半徑分別為0.10 m、0.15 m和0.30 m的入水過(guò)程比較，其中球體質(zhì)量為4 kg、入水角度為90°。2種方法得出的沖擊載荷峰值如表3所示。

表3 數(shù)值仿真與理論計(jì)算結(jié)果比較

隨著球體半徑的增大，數(shù)值仿真結(jié)果較理論計(jì)算值偏大。在理論計(jì)算中，理論模型未考慮結(jié)構(gòu)體入水沖擊過(guò)程中產(chǎn)生的水堆，水堆的出現(xiàn)使球體周圍等效水面抬升，導(dǎo)致入水沖擊載荷峰值減小。對(duì)于小體積球體，數(shù)值仿真與理論計(jì)算誤差較小，說(shuō)明所建立的結(jié)構(gòu)體-流體數(shù)值仿真模型計(jì)算結(jié)果具有一定的可信度，可用于后續(xù)研究。

2 參數(shù)影響分析

采用控制變量法進(jìn)行參數(shù)影響分析，主要分析不同剛性球體體積、不同入水速度、不同入水角度等參數(shù)對(duì)剛性球體入水沖擊載荷特性的影響。對(duì)于變體積工況，在保證質(zhì)量相同的情況下采用空心球模型。

2.1 體積參數(shù)影響

在球體質(zhì)量、入水速度和入水角度相同的情況下，采用中空球體，改變球體壁厚以此改變球體半徑，得到?jīng)_擊載荷隨著球體體積變化的規(guī)律如圖5所示。

圖5 沖擊載荷峰值隨球體體積變化的規(guī)律

球體體積越大，所受到的沖擊載荷越大，沖擊載荷與球體半徑的關(guān)系近似為

F=ARα+BRβ+CRγ+D

(5)

式中：F為球體入水沖擊載荷峰值；R為球體半徑，取值范圍為0.06～0.30 m。對(duì)模擬數(shù)據(jù)進(jìn)行數(shù)據(jù)擬合，得到如表4所示的參數(shù)。

表4 體積參數(shù)影響下的數(shù)據(jù)擬合參數(shù)

根據(jù)數(shù)值仿真結(jié)果，當(dāng)2個(gè)體積不同的模型進(jìn)行入水沖擊試驗(yàn)時(shí)，若球體半徑滿足式(6)，則2個(gè)模型入水沖擊載荷峰值滿足式(7)：

(6)

(7)

式(6)和式(7)中：λR為幾何縮尺比。

2.2 入水速度參數(shù)影響

在球體質(zhì)量、球體體積和入水角度相同的情況下，改變?nèi)胨俣?，得到?jīng)_擊載荷隨入水速度變化的規(guī)律如圖6所示，不同入水速度結(jié)構(gòu)體沖擊特性曲線對(duì)比如圖7所示，入水速度越大，所受到的沖擊載荷越大，但是入水沖擊載荷的峰值持續(xù)時(shí)間會(huì)縮短。對(duì)模擬數(shù)據(jù)進(jìn)行擬合，沖擊載荷與入水速度的關(guān)系近似為

圖6 沖擊載荷峰值隨球體入水速度變化的規(guī)律

圖7 不同入水速度結(jié)構(gòu)體沖擊特性曲線對(duì)比

F=11.621v2+27.851v+1 486.5

(8)

式中：v為球體入水速度，取值范圍為15～30 m/s。

根據(jù)數(shù)值仿真結(jié)果，當(dāng)2個(gè)入水速度不同的模型進(jìn)行入水沖擊試驗(yàn)時(shí)，若入水速度滿足式(9)，則2個(gè)模型入水沖擊載荷峰值滿足式(10)：

(9)

(10)

2.3 入水角度參數(shù)影響

在球體質(zhì)量、球體體積和入水速度相同的情況下，改變?nèi)胨嵌龋玫經(jīng)_擊載荷隨入水角度變化的規(guī)律如圖8所示，不同入水角度結(jié)構(gòu)體沖擊特性曲線對(duì)比如圖9所示。入水角度越大，所受到的沖擊載荷越大。對(duì)模擬數(shù)據(jù)進(jìn)行擬合，沖擊載荷與入水角度的關(guān)系近似為

圖8 沖擊載荷峰值隨球體入水角度變化的規(guī)律

圖9 不同入水角度結(jié)構(gòu)體沖擊特性曲線對(duì)比

F=0.000 3θ5-0.049 1θ4+3.413 6θ3-107.07θ2+1 547.8θ-4 791.3

(11)

式中：θ為球體入水角度，取值范圍為10°～60°。

根據(jù)數(shù)值仿真結(jié)果，當(dāng)2個(gè)入水角度不同的模型進(jìn)行入水沖擊試驗(yàn)時(shí)，若入水角度滿足式(12)，則2個(gè)模型入水沖擊載荷峰值滿足式(13)：

(12)

(13)

2.4 質(zhì)量速度耦合影響

采用中空的剛性球體，在球體體積不變的情況下，改變球體壁厚以控制球體的質(zhì)量。在球體體積、入水角度相同但球體的質(zhì)量和入水速度不相同的情況下，對(duì)模擬數(shù)據(jù)進(jìn)行擬合，球體所受到的沖擊載荷與球體質(zhì)量和入水速度的關(guān)系近似為

(14)

式中：λm和λv分別為質(zhì)量縮尺比和速度縮尺比。

2.5 縮尺比影響

在保證球體入水角度相同的情況下，球體質(zhì)量、體積、入水速度按照縮尺法則進(jìn)行參數(shù)設(shè)計(jì)，幾何縮尺比為λL，如式(15)和式(16)所示，根據(jù)數(shù)值仿真結(jié)果，球體所受到的沖擊載荷滿足式(17)：

(15)

(16)

(17)

3 結(jié) 論

基于計(jì)算流體力學(xué)方法建立結(jié)構(gòu)體高速入水仿真模型，并且與參考文獻(xiàn)中理論計(jì)算方法進(jìn)行對(duì)比，驗(yàn)證該仿真方法的有效性。在該有效的仿真模型中，比較不同體積、不同入水速度和不同入水角度等參數(shù)，得到不同參數(shù)對(duì)結(jié)構(gòu)物入水沖擊載荷特性的影響規(guī)律，主要結(jié)論如下：

(1) 將所建立的計(jì)算流體力學(xué)仿真模型與參考文獻(xiàn)中理論模型進(jìn)行比較，入水沖擊載荷峰值誤差隨著球體半徑的增大呈現(xiàn)增大的趨勢(shì)，但是對(duì)于小體積結(jié)構(gòu)物兩者吻合性較好，建立的仿真模型具有可信度，可用于后續(xù)研究。

(2) 入水沖擊載荷峰值在結(jié)構(gòu)體入水初期達(dá)到最大值，呈現(xiàn)出先增大后減小的趨勢(shì)。

(3) 采用控制變量法，得到剛性結(jié)構(gòu)體入水沖擊載荷峰值在一定范圍內(nèi)與結(jié)構(gòu)體入水速度、入水角度和結(jié)構(gòu)體體積等參數(shù)之間的定量關(guān)系，得出計(jì)算公式并且給出計(jì)算公式適用的范圍，該計(jì)算公式是從模型數(shù)據(jù)計(jì)算出實(shí)型數(shù)據(jù)的計(jì)算法則，有利于結(jié)構(gòu)體的初始設(shè)計(jì)。

建立一種計(jì)算結(jié)構(gòu)體高速入水載荷的CFD仿真方法，并為物理試驗(yàn)提供指導(dǎo)，縮短開(kāi)發(fā)周期，具有一定的理論意義和工程應(yīng)用價(jià)值。