圓柱體垂直入水三維數(shù)值模擬及影響因素研究

鄒田春， 高 飛， 魏家威， 解 江

(中國民航大學(xué) 適航學(xué)院，天津 300300)

結(jié)構(gòu)物入水沖擊涉及氣-固-液三者之間的相互作用，是典型的流固耦合過程，在結(jié)構(gòu)入水過程中，結(jié)構(gòu)與自由液面沖擊產(chǎn)生的較大沖擊力極易導(dǎo)致外部結(jié)構(gòu)破損，內(nèi)部構(gòu)件的失效，返回艙入水、魚雷入水、水上飛機的起降等都屬于結(jié)構(gòu)入水的范疇。因此結(jié)構(gòu)物入水的研究在軍事、船舶、海洋工程、航空航天等領(lǐng)域都有重要的工程應(yīng)用背景。

近十年來研究人員針對結(jié)構(gòu)入水沖擊問題的分析研究大多以楔形體[1-4]、圓柱體[5-8]、球體[9-10]等作為研究對象，采取試驗和仿真方法開展研究。Shah進(jìn)行了剛性楔形體入水沖擊試驗，利用高速攝像成功捕捉到了整個二維流體流動動態(tài)過程，并運用光滑粒子流體動力學(xué)(smoothed particle hydrodynamics，SPH)方法建立相應(yīng)模型，對SPH方法模擬結(jié)構(gòu)入水的準(zhǔn)確性進(jìn)行了驗證。夏維學(xué)等以圓柱體為研究對象開展入水試驗研究，對入水空泡壁面運動特性與空泡演化過程之間的關(guān)系進(jìn)行分析，并基于空泡壁面運動特性劃分空泡演化各階段，分析了空泡演化在各階段的流動機制。Luo等[11]對三維鋼質(zhì)楔形體進(jìn)行了一系列自由落體入水沖擊試驗，研究楔形體入水沖擊加速度、壓力及應(yīng)力響應(yīng)，并建立了楔形體試驗件的有限元模型，通過在模型上施加試驗測得的入水沖擊壓力來預(yù)測結(jié)構(gòu)的瞬態(tài)響應(yīng)，并對水彈性效應(yīng)進(jìn)行討論。施書文等[12]針對入水抨擊問題提出了一種弱可壓光滑粒子水動力與光滑點插值流固耦合法，運用該方法計算了彈性斜板低速入水和楔形體高速入水問題，計算結(jié)果與試驗值和半解析解吻合良好，證明了該方法對模擬彈性體入水問題的有效性。結(jié)構(gòu)入水仿真方法主要有有限體積法(finite volume method，F(xiàn)VM)[13-17]、SPH方法[18-20]以及任意拉格朗日歐拉(arbitrary Lagrangian-Eulerian，ALE)方法[21-26]等。ALE方法兼具拉格朗日法和歐拉法的優(yōu)點，在結(jié)構(gòu)入水問題研究領(lǐng)域采用較多且較為成熟。在傳統(tǒng)ALE理論和算法不變的情況下，S-ALE(structured-arbitrary Lagrange-Euler)方法則利用網(wǎng)格的規(guī)律性，將ALE模型規(guī)則網(wǎng)格通過輸入三維點位坐標(biāo)方式由程序本身自行創(chuàng)建網(wǎng)格，極大降低數(shù)值模擬的計算成本。

針對結(jié)構(gòu)入水問題的研究，在國防、軍事上大多以楔形體為研究對象，當(dāng)前我國正在大力發(fā)展民用客機且目前民用大型客機都是典型筒段機身結(jié)構(gòu)，考慮到民機機身水上迫降的工程背景，本文選擇圓柱體為研究對象?；趫A柱體垂直入水沖擊試驗研究結(jié)果，本文建立了圓柱體三維數(shù)值模型，采用S-ALE算法對圓柱體垂直入水沖擊過程進(jìn)行了數(shù)值模擬，分析其入水沖擊響應(yīng)特性，并通過對比圓柱體入水位移及自由液面動態(tài)變化來驗證模型有效性；基于驗證的數(shù)值模型，進(jìn)一步研究圓柱體入水速度、質(zhì)量和長度對其垂直入水沖擊的影響規(guī)律。

1 基本理論

1.1 S-ALE基礎(chǔ)方程

S-ALE在實現(xiàn)上是獨立于ALE而全新開發(fā)的，程序簡潔高效，但在理論方面S-ALE求解器與ALE完全相同。在ALE方法的描述中，引入了除Lagrange和Euler坐標(biāo)系之外的第三個任意參照坐標(biāo)系，與參照坐標(biāo)相關(guān)的材料微商可以采用式(1)描述

(1)

式中：Xi為Lagrange坐標(biāo)；xi為Euler坐標(biāo)；wi為相對速度。

假設(shè)物質(zhì)速度為v，網(wǎng)格速度為u，引入相對速度w=v-u。則ALE方法的控制方程可以由如下守恒方程給出

(1)質(zhì)量守恒方程

(2)

(2)動量守恒方程

控制固定域內(nèi)的牛頓流體流動問題的增強形式由控制方程和對應(yīng)的初始及邊界條件組成，控制流體問題的方程是Navier-Stokes方程的ALE描述

(3)

應(yīng)力張量為

σij=-pδij+μ(vi,j+vj,i)

(4)

式中：p為壓強；bi為單位質(zhì)量的體積力；δij為Kronecker δ函數(shù)；μ為動力黏性系數(shù)。該方程可與下列邊界條件和初始條件聯(lián)立求解

(5)

σijnj=0,Γ2邊界

(6)

Γ1∪Γ2=Γ,Γ1∩Γ2=0

(7)

式中：Г為計算域的完整邊界，分為Г1和Г2兩個部分，其中，Г1為固體邊界，Г2為自由液面邊界[27]；上標(biāo)0為該參數(shù)的初始指定值；nj為邊界的外法線單位向量。假設(shè)整個計算域在t=0時刻的速度場是靜止的，即

vi(xi,0)=0

(8)

(3)能量守恒方程

(9)

采用算子分離法求解ALE方程，該方法將每個時間步分為兩個階段。首先執(zhí)行拉格朗日過程，此時網(wǎng)格隨物質(zhì)一起運動。該過程計算速度、內(nèi)外力引起的內(nèi)能改變量，平衡方程為

(10)

(11)

單元邊界沒有物質(zhì)流過，所以滿足質(zhì)量守恒。第二個計算階段為對流階段，在對流階段之前采用VOF(volume of fluid)方法捕捉不同材料之間的界面位置，然后對穿過單元邊界的質(zhì)量運輸、內(nèi)能以及動量進(jìn)行計算，即將拉格朗日階段的移位網(wǎng)格重映射回初始位置或者任意位置。

對于每個節(jié)點，速度和位移按下列等式進(jìn)行更新

(12)

xn+1=xn+1+Δtun+1/2

(13)

外力矢量又分為體積力、邊界力、非反射邊界條件及接觸力。本文采用標(biāo)準(zhǔn)罰函數(shù)算法計算流體與結(jié)構(gòu)之間的接觸力。罰函數(shù)算法中允許從節(jié)點穿透到主段，所以在每一個時間步，系統(tǒng)首先檢查從節(jié)點是否穿透主面，未穿透前，不做任何處理，當(dāng)發(fā)生穿透，則引入一個與穿透距離d成比例的界面接觸力F并將其分配到流體單元節(jié)點上，計算公式為

F=ki·d

(14)

式中,ki為基于主從節(jié)點質(zhì)量模型特性的剛度系數(shù)。

1.2 水和空氣狀態(tài)方程

圓柱體垂直入水模型中，水和空氣采用空材料模型，通過一個截斷壓力值來定義壓力的更小邊界值，式(15)為任意方向黏性應(yīng)力的表達(dá)式。通過狀態(tài)方程來定義壓力，其中水采用Gruneisen狀態(tài)方程，空氣采用Polynomial狀態(tài)方程。

(15)

Gruneisen狀態(tài)方程對于壓縮材料由式(16)計算，膨脹材料由式(17)計算。

(γ0+au)E

(16)

P=ρ0C2u+(γ0+au)E

(17)

式中：C為沖擊波速度；E為流體的比內(nèi)能；ρ0為流體密度；γ0為Gruneisen伽馬；a為對γ0的一階體積修正；S1，S2，S3為曲線斜率的系數(shù)；由于材料不受初始載荷，初始能量項E和修正因子a均設(shè)為1。

線性多項式狀態(tài)方程

P=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E

(18)

式中：C0，C1，C2，C3，C4，C5，C6為自定義的常數(shù)；E為單位初始體積的內(nèi)能。

水的Gruneisen狀態(tài)方程參數(shù)，如表1所示?？諝獾腜olynomial狀態(tài)方程參數(shù)，如表2所示。

表1 水的Gruneisen狀態(tài)方程參數(shù)Tab.1 Parameters of Gruneisen state equation for water

表2 空氣的Polynomial狀態(tài)方程參數(shù)Tab.2 Parameters of Polynomial state equation for air

2 有限元建模

圓柱體垂直入水三維模型，如圖1所示。圓柱體直徑D=160 mm，長L=295 mm；水和空氣域長L=810 mm，寬W=660 mm，高H=700 mm其中水深h=480 mm。模型采用實體單元建模，圓柱體采用Lagrange單元，空氣和水采用Euler單元。通過定義ALE_STRUCTURED_MESH_CONTROL_POINTS和ALE_STRUCTURED_MESH兩個關(guān)鍵字在計算過程中由程序自動生成水和空氣的網(wǎng)格。

為防止圓柱體在試驗過程中發(fā)生變形，將PVC圓柱體內(nèi)部用泡沫材料進(jìn)行實心填充，所以圓柱體采用剛體進(jìn)行建模。水和空氣采用NULL材料模型和流體狀態(tài)方程來模擬，并對模型施加向下的重力加速度g，可以產(chǎn)生水的壓力梯度。

約束流體域底面單元節(jié)點的所有自由度來模擬水池底面邊界條件；在流體域四周設(shè)置無反射邊界條件，模擬無限水域，可減小流體域大小，降低計算成本，還可防止邊界產(chǎn)生的反射波與結(jié)構(gòu)發(fā)生流固耦合作用。運用S-ALE罰函數(shù)耦合算法計算圓柱體與流體之間的耦合作用。

基于Russo等開展的圓柱體入水試驗結(jié)果，選取質(zhì)量為m0(2.214 kg)試驗件提升1 m高度的入水試驗工況，如圖2所示。作為本文圓柱體入水沖擊仿真的初始工況。初始時刻圓柱體距液面高度為零，圓柱體初始入水速度為v0，通過運用GetData Graph Digitizer確定初始入水速度v0為3.6 m/s。將圓柱體入水沖擊仿真結(jié)果與試驗結(jié)果進(jìn)行對比驗證。

圖1 數(shù)值模擬模型Fig.1 Numerical simulation model

圖2 不同質(zhì)量試驗件提升高度與入水時刻速度關(guān)系圖Fig.2 Relation diagram of lifting height and water entry velocities of test pieces with different masses

3 有限元仿真及驗證

3.1 S-ALE與ALE計算對比

采用10 mm網(wǎng)格分別建立ALE模型和S-ALE模型，圓柱體與空氣、水網(wǎng)格比例為1 ∶1，設(shè)置計算總時長為0.25 s，使用LS-DYNA 971求解器在主頻為2.8 GHz的工作站上采用8核進(jìn)行計算。兩種模型在相同設(shè)置和工況下計算時間與最大位移值，如表3所示。S-ALE模型計算時間相對于ALE模型減少了20.2%。兩種模型計算得到的圓柱體位移-時間曲線對比，如圖3所示。從圖3可知，變化趨勢完全一致。其中，ALE模型得到最大位移值為0.231 6 m，S-ALE模型得到最大位移值為0.231 8 m，兩者相差可以忽略不計。因此，S-ALE模型在獲得相同仿真精度條件下可以大大減少計算時間。

表3 計算時間及最大位移值Tab.3 Calculation time and maximum displacement

圖3 S-ALE模型與ALE模型仿真結(jié)果對比Fig.3 Comparison of simulation results between S-ALE model and ALE model

3.2 模型驗證

圓柱體與空氣、水的網(wǎng)格尺寸相同，分別選取7.5 mm、10 mm、15 mm和20 mm進(jìn)行計算。采用S-ALE模型獲得的圓柱體入水位移-時間曲線與試驗結(jié)果對比，如圖4所示。

隨著網(wǎng)格的細(xì)化，仿真獲得的圓柱體最大位移值增加，且更接近試驗最大位移值(0.241 m)；網(wǎng)格細(xì)化到10 mm時，仿真獲得的圓柱體最大位移值趨于穩(wěn)定。當(dāng)網(wǎng)格尺寸為10 mm時，圓柱體最大位移值為0.231 m，與試驗值相差4.15%；網(wǎng)格尺寸為7.5 mm時，圓柱體最大位移值為0.235 m，與試驗值相差2.49%。因此，選用7.5 mm網(wǎng)格尺寸，模型總網(wǎng)格數(shù)為898 372，圓柱體不同入水時刻獲得的仿真位移值與試驗值對比情況，如表4所示。位移誤差均在3%以內(nèi)，因此S-ALE模型能夠精確模擬圓柱體入水位移時間歷程。

圖4 不同網(wǎng)格大小計算結(jié)果與試驗結(jié)果對比Fig.4 Comparison of calculation results of different grid sizes and test results

表4 不同時刻圓柱體仿真位移值與試驗結(jié)果對比Tab.4 The simulated displacement value of cylinder at different times compared with the experimental results

圓柱體入水仿真結(jié)果與試驗結(jié)果對比,如圖5所示。當(dāng)圓柱體接觸到自由液面時，圓柱體動能開始轉(zhuǎn)移為水的動能，圓柱體所擠占空間處的水在接觸處被噴射出來形成射流；圓柱體在下沉過程中高速的自由液面射流不斷增長與擴大(T<τ/5)(τ為圓柱體達(dá)到最大向下位移的時刻)；在(τ/5)

圖5 仿真結(jié)果和試驗結(jié)果對比Fig.5 Comparison of simulation and experimental results

S-ALE模型計算得到的液面運動情況與試驗之間有微小的差異，對于射流的飛濺效果，模擬結(jié)果不如試驗明顯，但仿真模擬基本上再現(xiàn)了圓柱體撞擊水面、自由液面射流形成、空氣腔打開、空氣腔閉合的典型現(xiàn)象。整體來看自由液面動態(tài)仿真模擬結(jié)果與試驗基本一致。

3.3 圓柱體入水速度分析

7.5 mm尺寸網(wǎng)格模型仿真得到的圓柱體入水速度變化曲線,如圖6所示。從圖6可知，圓柱體入水整個運動過程可分為撞擊水面、減速下移和上浮三個階段。圓柱體在撞擊水面時受到較大沖擊載荷，速度迅速下降，由于圓柱體密度小于水，在沖擊波衰減過后，其受水的作用力仍大于重力，因此圓柱體在撞擊水面后做加速度相對較小的減速運動，并在0.198 s時速度減小到0，此時圓柱體在水中處于懸浮狀態(tài)，隨后在浮力的作用下圓柱體開始上浮。

圖6 圓柱體入水速度-時間曲線Fig.6 Velocity-time history curve of cylinder water entry

4 不同因素影響研究

4.1 圓柱體入水速度影響

圓柱體(D=160 mm，L=295 mm)不同入水速度下的加速度-時間歷程曲線，如圖7所示。從圖7可知，三種速度下加速度曲線都有一較大峰值，此時段為圓柱體沖擊水面時段，在較大沖擊載荷作用下會對結(jié)構(gòu)產(chǎn)生較大破壞，也是結(jié)構(gòu)物入水研究關(guān)注的重點，對比加速度曲線可知，入水速度對圓柱體加速度峰值影響非常明顯，隨著速度的增大加速度峰值增大，而對加速度脈沖持續(xù)時間幾乎沒有影響。

圖7 圓柱體不同入水速度加速度-時間曲線Fig.7 Acceleration-time history curves of cylinder with different water entry velocities

圓柱體(D=160 mm，L=295 mm)在不同速度下的入水位移-時間歷程曲線，如圖8所示。從圖8可知，三種速度下圓柱體位移趨勢一致，隨著入水速度增大，圓柱體向下運動的越快且最大位移增大，但達(dá)到最大位移的時刻基本保持不變。

圖8 圓柱體不同入水速度位移-時間曲線Fig.8 Displacement-time history curves of cylinder with different water entry velocities

4.2 圓柱體質(zhì)量影響

通過在圓柱體上方增加集中質(zhì)量的方式改變圓柱體質(zhì)量從而獲得2.214 kg，3.214 kg，4.214 kg，5.214 kg圓柱體入水沖擊加速度。圖9為圓柱體(D=160 mm，L=295 mm)入水速度3.6 m/s，不同質(zhì)量工況下加速度曲線，如圖9所示。從圖9可知，隨著質(zhì)量增大，加速度峰值有減小趨勢，但質(zhì)量增加未改變加速度脈沖持續(xù)時間。圓柱體受力示意圖如圖10所示。由牛頓第二定律式(19)可以得到圓柱體入水沖擊在垂直方向上所受水的合外力F(包括浮力、動量傳輸引起的附加力、形狀阻力)，不同質(zhì)量圓柱體入水所受阻力曲線，如圖11所示。由圖11可知，在圓柱體外形不變的條件下不同質(zhì)量圓柱體入水沖擊所受水的合外力基本不變，則由式(20)可知，外形相同的圓柱體入水沖擊時，質(zhì)量m越大，加速度a越小。

F-mg=ma

(19)

a=F/m-g

(20)

式中：m為圓柱體質(zhì)量；g為重力加速度；a為圓柱體加速度。

圖9 不同質(zhì)量圓柱體入水加速度-時間曲線Fig.9 The water entry acceleration-time history curves of cylinders with different masses

不同質(zhì)量圓柱體(D=160 mm，L=295 mm)以2.2 m/s速度入水時位移-時間歷程曲線，如圖12所示。從圖12可知，由于加速度大小與質(zhì)量呈負(fù)相關(guān)，則質(zhì)量越大圓柱體入水減速越緩，向下運動的越快；隨著質(zhì)量增加，圓柱體入水最大位移及其到達(dá)時刻都增大；當(dāng)質(zhì)量增加到圓柱體重力大于浮力時，圓柱體入水過程最終將會達(dá)到穩(wěn)定狀態(tài)，保持勻速一直下落。

圖10 圓柱體入水沖擊時的受力圖Fig.10 The force diagram of cylinder impact water

圖11 不同質(zhì)量圓柱體入水所受阻力曲線Fig.11 The water entry drag force curves of the cylinders with different masses

圖12 不同質(zhì)量圓柱體入水位移-時間歷程曲線Fig.12 The water entry displacement-time history curves of cylinders with different masses

4.3 圓柱體長度影響

為研究圓柱體入水三維效應(yīng)，針對給定直徑D的圓柱體，保持其密度ρ不變，只改變長度L，研究長度對圓柱體入水特性的影響。

選取D=160 mm圓柱體，不同長度圓柱體入水加速度-時間歷程曲線，如圖13所示。不同長度(D=160 mm)圓柱體加速度峰值大小，如圖14所示。從圖14可知，加速度峰值隨著長度的增大而增大，且當(dāng)L/D>1時，增加趨勢逐漸變緩，當(dāng)L/D>3時，L/D的增加對加速度峰值幾乎沒有影響。L/D對加速度脈沖持續(xù)時間幾乎沒有影響。相對于速度，長度改變所造成的加速度變化幅度明顯要小。

圖13 不同長度(D=160 mm)圓柱體入水加速度-時間曲線Fig.13 The water entry acceleration-time history curves of cylinders(D=160 mm)of different lengths

圖14 不同長度(D=160 mm)圓柱體入水加速度峰值Fig.14 The water entry peak accelerations of cylinders (D=160 mm) with different lengths

不同長度(D=80 mm)圓柱體入水位移-時間歷程曲線，如圖15所示。從圖15可知，隨著L/D增大，圓柱體最大入水位移與其達(dá)到最大位移的時刻都增大；當(dāng)L/D>2時，L/D的增大對最大位移的影響逐漸降低；當(dāng)L/D>4時，L/D的增加不再影響最大位移。

圖15 不同長度(D=80 mm)圓柱體入水位移-時間曲線Fig.15 The water entry displacement-time history curves of cylinders(D=80 mm)with different lengths

不同長度(D=120 mm)圓柱體入水位移-時間歷程曲線，如圖16所示。當(dāng)L/D<3時，隨著L/D的增加，圓柱體最大入水位移與其達(dá)到最大位移的時刻都增大；當(dāng)L/D>3時，隨著L/D的增加，最大入水位移開始減小。

圖16 不同長度(D=120 mm)圓柱體入水位移-時間曲線Fig.16 The water entry displacement-time curves of cylinders(D=120 mm)with different lengths

圖17給出了不同長度(D=160 mm)圓柱體入水位移-時間歷程曲線，當(dāng)L/D<2時，圓柱體最大入水位移與其達(dá)到最大位移的時刻都增大；當(dāng)L/D>2時，隨著L/D的增加，最大入水位移逐漸減小。

圖17 不同長度(D=160 mm)圓柱體入水位移-時間曲線Fig.17 The water entry displacement-time curves of cylinders(D=160 mm)with different lengths

L/D對圓柱體入水的影響在開始時間階段較小，且L/D越小向下運動的越快，隨著向下運動的進(jìn)行，L/D對圓柱體位移的影響變的更加明顯，且因直徑的不同而呈現(xiàn)不同的規(guī)律。

圓柱體最大入水位移與不同L/D的關(guān)系曲線，如圖18所示。從圖18可知，隨著直徑的增大，最大位移增大；三種直徑圓柱體，在L/D<2時，L/D對位移的影響較大；當(dāng)L/D>2時，影響逐漸變緩，且不同直徑呈現(xiàn)不同的影響規(guī)律。

圖18 不同長度、不同直徑圓柱體入水位移曲線Fig.18 The water entry displacement curves of cylinders with different lengths and diameters

5 結(jié) 論

(1)相對ALE模型，S-ALE模型在保證計算精度的前提下，可大大提升計算效率。S-ALE模型可以完整的再現(xiàn)圓柱體沖擊水面、自由液面射流形成、空氣腔打開、空氣腔閉合的整個過程。通過圓柱體入水過程典型時刻自由液面動態(tài)以及位移的數(shù)值計算結(jié)果與試驗結(jié)果對比，驗證了本文所建立模型的有效性和準(zhǔn)確性。

(2)隨著圓柱體入水沖擊速度的增大，圓柱體加速度峰值及最大位移值都增大，但對加速度脈沖持續(xù)時間和達(dá)到最大位移時刻沒有影響。

(3)加速度大小與質(zhì)量呈負(fù)相關(guān)，但圓柱體質(zhì)量改變，加速度脈沖持續(xù)時間基本不變，質(zhì)量增大最大位移及其到達(dá)時刻增大。

(4)長度增加，圓柱體入水沖擊加速度增加，但對加速度脈沖持續(xù)時間沒有影響，且相對于速度，長度改變引起的加速度變化幅度明顯要?。蝗胨灰圃陂_始時間段受長度的影響較小，且L/D越小圓柱體向下運動的越快，隨著圓柱體向下運動的進(jìn)行，L/D對圓柱體位移的影響變的更加明顯，且因直徑的不同而呈現(xiàn)不同的規(guī)律。