曹力飛,李炳彥,劉 艷

(1.西安機(jī)電信息技術(shù)研究所,陜西 西安 710065;2.中國(guó)航發(fā)成都發(fā)動(dòng)機(jī)有限公司,四川 成都 610503)

0 引言

入水沖擊過程是一種典型的流固耦合問題,針對(duì)流固耦合的研究從20世紀(jì)30年代就已經(jīng)開始,1929年Von Karman對(duì)楔形體和平頭體的入水載荷進(jìn)行初步計(jì)算,之后Wagner對(duì)Von Karman的算法進(jìn)行改進(jìn),用于完成分析水面抬升現(xiàn)象,提出了銷斜面升角理論。發(fā)展至今對(duì)于流固耦合的力學(xué)分析計(jì)算主要有兩種方案:CFD仿真方法和ALE流固耦合方法。CFD仿真方案主要是通過動(dòng)網(wǎng)格方法[1-2],實(shí)現(xiàn)作用過程中流體與固體的實(shí)時(shí)耦合計(jì)算[3],但該方案存在以下問題:1) 固體運(yùn)動(dòng)幅度太大會(huì)導(dǎo)致計(jì)算的域過大,進(jìn)而致使網(wǎng)格數(shù)量太大,大幅增加計(jì)算周期與計(jì)算壓力;2) 流體域的網(wǎng)格劃分需要對(duì)固體部分進(jìn)行掏空與綁定,對(duì)于復(fù)雜固體的網(wǎng)格劃分難度較大[4]。例如:機(jī)翼在流場(chǎng)中的擺動(dòng)、橋梁在風(fēng)場(chǎng)中的震動(dòng)等。ALE流固耦合方法的特點(diǎn)是流體與固體單獨(dú)進(jìn)行網(wǎng)格劃分,計(jì)算過程中網(wǎng)格流體網(wǎng)格不發(fā)生變化,只在耦合時(shí)在耦合界面進(jìn)行物理量的交換[5]。例如:流體在箱體中的運(yùn)動(dòng)、高超聲速入水、流水沖擊固體等案例[7]。

針對(duì)結(jié)構(gòu)體入水的研究主要有CFD仿真方法和ALE流固耦合方法:任意拉格朗日-歐拉(ALE)方法分析結(jié)構(gòu)體對(duì)于流體的流動(dòng)特性,該方法可以分析高速入水與沖擊過程,可以確定結(jié)構(gòu)體入水的整個(gè)過程的力學(xué)特性,但是計(jì)算量過大。CFD算法關(guān)于流體對(duì)于結(jié)構(gòu)體的粘性摩擦分布特性的研究,使用范圍主要集中在連續(xù)流體領(lǐng)域和低速流體領(lǐng)域,對(duì)于流-固高速?zèng)_擊問題無法計(jì)算,并且對(duì)于計(jì)算的網(wǎng)格質(zhì)量要求較高,否則容易產(chǎn)生發(fā)散。

隨著工程系統(tǒng)的復(fù)雜程度越來越高,網(wǎng)格的數(shù)量越來越大。一方面,更多的網(wǎng)格數(shù)量可以得到更加逼近真實(shí)的仿真環(huán)境;另一方面,更多的網(wǎng)格數(shù)量又勢(shì)必會(huì)降低仿真效率。為了提高計(jì)算效率,對(duì)仿真服務(wù)器和計(jì)算機(jī)性能要求越來越高,影響產(chǎn)品研制周期,通過網(wǎng)格模型優(yōu)化的方法已經(jīng)越來越無法滿足工程計(jì)算的需求。對(duì)采用定步長(zhǎng)ALE算法的大變形、大數(shù)據(jù)量交換的結(jié)構(gòu)體高速入水耦合計(jì)算影響更加明顯。針對(duì)上述的不足,本文提出基于沖擊修正的任意拉格朗日-歐拉(C-ALE)算法的結(jié)構(gòu)體入水計(jì)算方法。

1 任意拉格朗日歐拉(ALE)算法

ALE(arbitrary Lagrange-Euler)連續(xù)介質(zhì)力學(xué)有兩種經(jīng)典的運(yùn)動(dòng)描述方法,分別是Lagrange描述和Euler描述。Lagrange描述:網(wǎng)格節(jié)點(diǎn)固定在固定物質(zhì)點(diǎn)并隨之運(yùn)動(dòng),因此在描述運(yùn)動(dòng)邊界或者運(yùn)動(dòng)界面時(shí)非常方便;但當(dāng)物質(zhì)發(fā)生大變形時(shí)常常產(chǎn)生網(wǎng)格糾纏,輕則影響單元近似精度,重則使坐標(biāo)中的Jocobian行列式的值近似于零或負(fù)數(shù),從而使計(jì)算終止,或者引起局部嚴(yán)重誤差。Euler描述:網(wǎng)格節(jié)點(diǎn)固定在空間始終不動(dòng),因此在描述大變形時(shí),不會(huì)產(chǎn)生糾纏問題;但是,存在網(wǎng)格和物質(zhì)的相對(duì)處理對(duì)流效應(yīng)更加困難和無法精確確定運(yùn)動(dòng)邊界的問題。為了克服Lagrange描述和Euler描述各自存在的缺點(diǎn),引入相對(duì)坐標(biāo)任意移動(dòng)的方法即ALE方法。

結(jié)構(gòu)體的入水過程是一個(gè)非線性過程,流體域包含了氣液多相流流動(dòng),ALE算法主要包含質(zhì)量、動(dòng)量和能量守恒方程。

質(zhì)量守恒方程:

動(dòng)量守恒方程:

式中,X為ALE坐標(biāo);x為空間坐標(biāo);ci為ALE描述的對(duì)流速度,ci=vi-wi;vi為流體指點(diǎn)的物質(zhì)速度;wi為網(wǎng)格速度;ρw為流體速度;bi為流體體積;σij為應(yīng)力漲量。

能量守恒方程:

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

(3)

式(3)中,ps為水的靜壓;μd為動(dòng)力粘性系數(shù);δij為克羅內(nèi)常數(shù)。

2 C-ALE算法的結(jié)構(gòu)體入水計(jì)算方法

本文基于C-ALE算法實(shí)現(xiàn)結(jié)構(gòu)體入水計(jì)算,其原理是在任意拉格朗日-歐拉(ALE)方法的基礎(chǔ)上引入貫穿量C。求解過程如圖1所示,采用交替計(jì)算的方法求解流固耦合問題,將流體和固體劃分為兩個(gè)獨(dú)立的求解域,在耦合界面進(jìn)行物理參數(shù)交換求解之后,分別對(duì)流體域和固體域進(jìn)行單獨(dú)計(jì)算。既保證流體域和固體域的計(jì)算可靠性,又保證耦合區(qū)間的物理量傳遞。利用C-ALE方法單獨(dú)求解流體域,耦合界面上的結(jié)點(diǎn)力傳遞給結(jié)構(gòu),作為固體域的力邊界條件,單獨(dú)求解動(dòng)力學(xué)方程。當(dāng)流體域和固體域同時(shí)達(dá)到精度要求時(shí)結(jié)束迭代,根據(jù)耦合界面的運(yùn)動(dòng)情況更新單元網(wǎng)格進(jìn)入下一時(shí)間步的求解。

圖1 C-ALE流固耦合計(jì)算模型Fig.1 C-ALE fluid-solid calculation model

具體計(jì)算過程中,檢查每個(gè)計(jì)算步長(zhǎng)的網(wǎng)格節(jié)點(diǎn)是否有物質(zhì)穿透,沒有則對(duì)該網(wǎng)格不做處理。如果發(fā)生物質(zhì)穿透,則通過界面力F參與計(jì)算。計(jì)算過程中界面力與貫穿量C成正比:

F=kt·C,

(4)

式(4)中,kt物質(zhì)模型剛度系數(shù),C貫穿量。

(5)

式(5)中,ΔEm為動(dòng)能變化量,Fd為沖擊阻力,s水流沖擊距離。

通過貫穿量C的引入,加快了在耦合界面的能量傳遞參數(shù)的計(jì)算速度,尤其在高速流體沖擊過程中,對(duì)于結(jié)構(gòu)體的沖擊應(yīng)力的計(jì)算。但是,對(duì)于低速入水的結(jié)構(gòu)體入水過程影響不大。

2.1 算例模型

入水模型由入水結(jié)構(gòu)體、彈體、空氣和水組成,如圖2所示。其中彈體和入水結(jié)構(gòu)體作為固體域參與計(jì)算,水和空氣作為流體域參與計(jì)算。耦合界面包含水-空氣流體耦合、彈體與入水結(jié)構(gòu)體固體耦合、入水結(jié)構(gòu)體和水組成流固耦合。

圖2 結(jié)構(gòu)體入水仿真模型Fig.2 Structure water entry simulation model

流體域主要采用歐拉網(wǎng)格進(jìn)行描述,而固體域則采用拉格朗日網(wǎng)格描述。材料網(wǎng)格之間可以進(jìn)行材料交換,而不需要產(chǎn)生網(wǎng)格運(yùn)行與生成[8]。對(duì)固體和流體域在建立耦合邊界條件之后,分別進(jìn)行網(wǎng)格劃分。網(wǎng)格劃分如圖3、圖4所示。

圖3 流體域網(wǎng)格劃分Fig.3 Fluid domain meshing

通過粘性應(yīng)力本構(gòu)關(guān)系,建立水的Gruneisen狀態(tài)方程材料模型,模型參數(shù)如表1所示。

(6)

式(6)中,μ為比體積;p為壓力;E為單位體積內(nèi)能;C,s1,s2,s3為材料常數(shù);γ0為水的密度;a為對(duì)γ0的一階體積修正。

圖4 固體域網(wǎng)格劃分Fig.4 Solid domain meshing

表1 水和空氣的Gruneisen狀態(tài)方程材料模型參數(shù)Tab.1 Gruneisen equation of state material model parameters of water and air

入水結(jié)構(gòu)體選用鋁合金2A12,材料模型選取Johnson-Cook模型,模型參數(shù)如表2所示。Johnson-Cook模型本質(zhì)是將應(yīng)變、應(yīng)變率和溫度三個(gè)變量進(jìn)行分離,利用乘積關(guān)系處理三者對(duì)動(dòng)態(tài)屈服應(yīng)力的影響,具有形式簡(jiǎn)單、各項(xiàng)物理意義明確的優(yōu)點(diǎn)。

表2 2A12材料Johnson-Cook模型參數(shù)Tab.2 Johnson-Cook model parameters of 2A12

彈體部分處理為質(zhì)量單元,在不影響計(jì)算準(zhǔn)確性的前提下以增加計(jì)算效率。

2.2 流體域計(jì)算準(zhǔn)確性

通過空泡理論對(duì)仿真的準(zhǔn)確性進(jìn)行驗(yàn)證。入水速度取350 m/s。對(duì)于結(jié)構(gòu)體入水難以準(zhǔn)確測(cè)量,只能通過空泡的外形進(jìn)行對(duì)比驗(yàn)證。基于空泡理論可以相對(duì)準(zhǔn)確的計(jì)算入水過程的空泡外形[9]。

式(7)中,db為彈體直徑;Cd為入水阻力系數(shù),該系數(shù)與入水的速度有關(guān)。

根據(jù)采用C-ALE算法進(jìn)行對(duì)比可以得到仿真空泡輪廓,彈丸入水過程體積分?jǐn)?shù)變化如圖5所示。

圖5 彈丸入水過程體積分?jǐn)?shù)變化Fig.5 Variation of volume fraction of projectile during water

C-ALE算法與空泡理論計(jì)算對(duì)比曲線如圖6所示。可知基于C-ALE算法入水計(jì)算曲線與空泡理論曲線特征相近,空泡深度空泡和開口尺寸基本一致,因此,基于C-ALE算法的結(jié)構(gòu)體入水計(jì)算結(jié)果準(zhǔn)確,算法合理可行。

圖6 理論空泡曲線與計(jì)算空泡曲線對(duì)比Fig.6 Comparison of theoretical and calculated cavitation curves

2.3 固體域計(jì)算準(zhǔn)確性

以350 m/s的速度進(jìn)行基于沖擊修正的任意拉格朗日-歐拉(C-ALE)方法的結(jié)構(gòu)體入水?dāng)?shù)值分析。入水過程結(jié)構(gòu)體速度和加速度變化曲線如圖7、圖8所示。

仿真結(jié)果顯示:彈丸入水過程結(jié)構(gòu)體速度穩(wěn)定下降,原因是由于空泡的產(chǎn)生只有入水結(jié)構(gòu)體的部分接觸水面并產(chǎn)生阻力;至6 ms后下降速度加快,分析原因?yàn)殄F形段部分開始浸水導(dǎo)致摩擦阻力增加。

圖7 入水過程結(jié)構(gòu)體速度變化Fig.7 Variation of structure velocity during water entry

入水結(jié)構(gòu)體的變形主要發(fā)生在入水結(jié)構(gòu)體接觸水面的瞬間,此時(shí)過載雖然沒有形成最大過載,但是此時(shí)接觸面積較小。

圖8 入水過程彈丸加速度變化Fig.8 Variation of projectile acceleration during water entry

在入水結(jié)構(gòu)體中形成超過極限屈服載荷的壓力,入水結(jié)構(gòu)體產(chǎn)生形變。最大形變?yōu)?.520 6 mm,時(shí)間發(fā)生在0.122 ms,與入水過載峰值時(shí)間基本保持一致。隨后形變發(fā)生一定回彈,并保持穩(wěn)定,最終形變?yōu)?.231 1 mm。入水結(jié)構(gòu)體的變形階段主要發(fā)生在入水結(jié)構(gòu)體與水接觸的初期。當(dāng)入水結(jié)構(gòu)體侵入水面之后,隨著錐形彈丸與水的接觸面積的增加,變形不再變化。

根據(jù)圖9入水前后結(jié)構(gòu)體形變可知:仿真結(jié)果顯示形變主要發(fā)生在結(jié)構(gòu)體的頭部,在錐形段部分沒有發(fā)生相關(guān)形變。

外場(chǎng)試驗(yàn)時(shí),測(cè)試彈丸以350 m/s的速度進(jìn)行入水沖擊試驗(yàn),試驗(yàn)前后入水結(jié)構(gòu)體變形情況如圖10所示。

試驗(yàn)結(jié)果表明:1) 結(jié)構(gòu)體入水的形變主要發(fā)生在結(jié)構(gòu)體的頭部,這與仿真結(jié)果保持一致;2) 外場(chǎng)試驗(yàn)實(shí)測(cè)產(chǎn)品產(chǎn)生3.10 mm的變形,仿真計(jì)算入水結(jié)構(gòu)體產(chǎn)生3.231 1 mm變形,仿真計(jì)算與外場(chǎng)試驗(yàn)結(jié)果基本一致。

圖9 入水后結(jié)構(gòu)體形變?cè)茍DFig.9 Deformation of structure after water entry

圖10 入水前后結(jié)構(gòu)體變形(試驗(yàn)結(jié)果)Fig.10 Variation of projectile acceleration during water entry with structure

驗(yàn)證了沖擊修正的任意拉格朗日-歐拉(C-ALE)方法對(duì)錐形殼體入水沖擊計(jì)算準(zhǔn)確,可以用于指導(dǎo)結(jié)構(gòu)體高速入水的仿真。

3 適應(yīng)性與數(shù)據(jù)分析

C-ALE算法與ALE算法采用同樣的網(wǎng)格模型,因此可以通過ALE算法進(jìn)一步驗(yàn)證C-ALE算法對(duì)結(jié)構(gòu)體入水過程仿真的適用性。

3.1 C-ALE算法與結(jié)構(gòu)體形變收斂特性

由圖11和表3可知,采用同一網(wǎng)格模型,在相同輸入條件下,通過ALE算法進(jìn)行計(jì)算。ALE算法計(jì)算收斂時(shí)間約為8 976 s,收斂后,入水結(jié)構(gòu)體的結(jié)構(gòu)形變?yōu)?.76 mm。C-ALE算法計(jì)算收斂時(shí)間約為3 342 s,入水結(jié)構(gòu)體變形3.230 9 mm。相對(duì)誤差從21.2%減小到4.2%,絕對(duì)誤差由0.66 mm減小到0.13 mm,收斂時(shí)間8 976 s減小到3 342 s。可以得出:針對(duì)結(jié)構(gòu)體入水計(jì)算過程,相同輸入條件下,C-ALE算法比ALE算法具有更小的計(jì)算誤差和更快的收斂速度。

圖11 不同計(jì)算方法的收斂過程Fig.11 Convergence process of different calculation methods

表3 不同計(jì)算方法的結(jié)構(gòu)體入水計(jì)算結(jié)果Tab.3 Calculation results of water inflow of structure with different calculation

3.2 入水深度與計(jì)算時(shí)間的關(guān)系

為了進(jìn)一步驗(yàn)證C-ALE算法對(duì)于結(jié)構(gòu)體入水深度的計(jì)算時(shí)間,在350 m/s入水速度條件下,采用同樣的網(wǎng)格和材料模型,對(duì)比ALE算法對(duì)結(jié)構(gòu)體入水50 mm深度的仿真過程進(jìn)行監(jiān)測(cè),獲得C-ALE和ALE算法入水深度與計(jì)算時(shí)間關(guān)系如圖12所示。

圖12 C-ALE和ALE算法入水深度與計(jì)算時(shí)間關(guān)系Fig.12 C-ALE和ALE algorithm relationship between water depth and calculation time

根據(jù)圖12可知,在高速入水過程中,隨著入水深度的增加,C-ALE算法相比ALE算仿真時(shí)間大幅減少。

3.3 入水速度與計(jì)算時(shí)間的關(guān)系

為確定不同入水速度對(duì)于算法仿真速度的影響,入水速度分別取50、100、150、200、250、300、350 m/s,采用同樣的網(wǎng)格和材料模型,對(duì)比ALE算法對(duì)結(jié)構(gòu)體入水50 mm深度的仿真過程進(jìn)行監(jiān)測(cè),獲得不同入水速度C-ALE和ALE下入水深度與計(jì)算時(shí)間如圖13所示。

圖13 不同入水速度C-ALE和ALE入水深度與計(jì)算時(shí)間Fig.13 Different entry velocity C-ALE and ALE entry depth and calculation time

由圖13可知:1) 在相同網(wǎng)格和材料模型下,不同入水速度的C-ALE仿真過程均比ALE算法的仿真時(shí)間少;2) 在相同網(wǎng)格和材料模型下,當(dāng)入水速度大于150 m/s,C-ALE算法與ALE算法的仿真時(shí)間差距增大。50 mm深度不同入水速度下結(jié)構(gòu)體形變與計(jì)算時(shí)間如表4所示。

表4 50 mm深度不同入水速度下結(jié)構(gòu)體形變與計(jì)算時(shí)間Tab.4 Structure deformation and calculation time with different water entry velocity at 50mm depth

根據(jù)表4可知,入水速度在150 m/s時(shí)候,結(jié)構(gòu)體開始產(chǎn)生形變,隨著入水速度的增加,結(jié)構(gòu)體的形變?cè)黾?C-ALE算法中引入貫穿量C對(duì)于仿真的影響增大,對(duì)比ALE算法,進(jìn)一步縮短了仿真時(shí)間。

3.4 網(wǎng)格數(shù)量與收斂性

進(jìn)一步研究不同網(wǎng)格數(shù)量對(duì)C-ALE算法在入水仿真過程的收斂特性的影響,分別對(duì)如下網(wǎng)格特性進(jìn)行仿真計(jì)算。A網(wǎng)格特性:流體域網(wǎng)格數(shù)量675,固體域網(wǎng)格數(shù)量373;B網(wǎng)格特性:流體域網(wǎng)格數(shù)量2 325,固體域網(wǎng)格數(shù)量1 562;C網(wǎng)格特性:流體域網(wǎng)格數(shù)量6 566,固體域網(wǎng)格數(shù)量2 987。不同網(wǎng)格特性的收斂過程如圖14所示。仿真條件為350 m/s入水沖擊過程,實(shí)際形變?yōu)闆_擊試驗(yàn)形變。C-ALE算法下不同網(wǎng)格特性下結(jié)構(gòu)體形變收斂時(shí)間如表5所示。

圖14 不同網(wǎng)格特性的收斂過程Fig.14 Convergence process of different mesh methods

表5 不同網(wǎng)格特性下結(jié)構(gòu)體形變與收斂時(shí)間Tab.5 Structure deformation and calculation time with different mesh methods

由圖14和表5可知,不同流體域與固體域的網(wǎng)格數(shù)量?jī)H影響仿真的收斂速度,對(duì)仿真的準(zhǔn)確度沒有影響。

4 結(jié)論

本文提出了基于沖擊修正的任意拉格朗日-歐拉(C-ALE)算法的結(jié)構(gòu)體入水計(jì)算方法,該方法引入貫穿量C對(duì)任意拉格朗日-歐拉方法進(jìn)行修正,采用流體域與固體域定步長(zhǎng)交替計(jì)算的方法進(jìn)行求解。通過空泡理論和外場(chǎng)試驗(yàn)結(jié)果驗(yàn)證了C-ALE計(jì)算結(jié)果的準(zhǔn)確性;同時(shí),對(duì)比采用相同輸入條件和網(wǎng)格模型的ALE算法可以得出,對(duì)于結(jié)構(gòu)體入水過程的求解,C-ALE算法比ALE算法具有更快的收斂速度和更小的計(jì)算誤差。在相同網(wǎng)格模型和材料模型下,不同的入水速度,隨著結(jié)構(gòu)體形變的增加,C-ALE算法可以進(jìn)一步縮減仿真時(shí)間?？梢缘贸鋈缦陆Y(jié)論:通過貫穿量C的引入,加快了耦合界面的參數(shù)傳遞效率,提高了仿真計(jì)算速度。