基于修正任意拉格朗日歐拉算法的結構體入水仿真

曹力飛,李炳彥,劉 艷

(1.西安機電信息技術研究所,陜西 西安 710065;2.中國航發成都發動機有限公司,四川 成都 610503)

0 引言

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

針對結構體入水的研究主要有CFD仿真方法和ALE流固耦合方法:任意拉格朗日-歐拉(ALE)方法分析結構體對于流體的流動特性,該方法可以分析高速入水與沖擊過程,可以確定結構體入水的整個過程的力學特性,但是計算量過大。CFD算法關于流體對于結構體的粘性摩擦分布特性的研究,使用范圍主要集中在連續流體領域和低速流體領域,對于流-固高速沖擊問題無法計算,并且對于計算的網格質量要求較高,否則容易產生發散。

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

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

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

結構體的入水過程是一個非線性過程,流體域包含了氣液多相流流動,ALE算法主要包含質量、動量和能量守恒方程。

質量守恒方程:

動量守恒方程:

式中,X為ALE坐標;x為空間坐標;ci為ALE描述的對流速度,ci=vi-wi;vi為流體指點的物質速度;wi為網格速度;ρw為流體速度;bi為流體體積;σij為應力漲量。

能量守恒方程:

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

(3)

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

2 C-ALE算法的結構體入水計算方法

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

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

具體計算過程中,檢查每個計算步長的網格節點是否有物質穿透,沒有則對該網格不做處理。如果發生物質穿透,則通過界面力F參與計算。計算過程中界面力與貫穿量C成正比:

F=kt·C,

(4)

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

(5)

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

通過貫穿量C的引入,加快了在耦合界面的能量傳遞參數的計算速度,尤其在高速流體沖擊過程中,對于結構體的沖擊應力的計算。但是,對于低速入水的結構體入水過程影響不大。

2.1 算例模型

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

圖2 結構體入水仿真模型Fig.2 Structure water entry simulation model

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

圖3 流體域網格劃分Fig.3 Fluid domain meshing

通過粘性應力本構關系,建立水的Gruneisen狀態方程材料模型,模型參數如表1所示。

(6)

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

圖4 固體域網格劃分Fig.4 Solid domain meshing

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

入水結構體選用鋁合金2A12,材料模型選取Johnson-Cook模型,模型參數如表2所示。Johnson-Cook模型本質是將應變、應變率和溫度三個變量進行分離,利用乘積關系處理三者對動態屈服應力的影響,具有形式簡單、各項物理意義明確的優點。

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

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

2.2 流體域計算準確性

通過空泡理論對仿真的準確性進行驗證。入水速度取350 m/s。對于結構體入水難以準確測量,只能通過空泡的外形進行對比驗證。基于空泡理論可以相對準確的計算入水過程的空泡外形[9]。

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

根據采用C-ALE算法進行對比可以得到仿真空泡輪廓,彈丸入水過程體積分數變化如圖5所示。

圖5 彈丸入水過程體積分數變化Fig.5 Variation of volume fraction of projectile during water

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

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

2.3 固體域計算準確性

以350 m/s的速度進行基于沖擊修正的任意拉格朗日-歐拉(C-ALE)方法的結構體入水數值分析。入水過程結構體速度和加速度變化曲線如圖7、圖8所示。

仿真結果顯示:彈丸入水過程結構體速度穩定下降,原因是由于空泡的產生只有入水結構體的部分接觸水面并產生阻力;至6 ms后下降速度加快,分析原因為錐形段部分開始浸水導致摩擦阻力增加。

圖7 入水過程結構體速度變化Fig.7 Variation of structure velocity during water entry

入水結構體的變形主要發生在入水結構體接觸水面的瞬間,此時過載雖然沒有形成最大過載,但是此時接觸面積較小。

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

在入水結構體中形成超過極限屈服載荷的壓力,入水結構體產生形變。最大形變為3.520 6 mm,時間發生在0.122 ms,與入水過載峰值時間基本保持一致。隨后形變發生一定回彈,并保持穩定,最終形變為3.231 1 mm。入水結構體的變形階段主要發生在入水結構體與水接觸的初期。當入水結構體侵入水面之后,隨著錐形彈丸與水的接觸面積的增加,變形不再變化。

根據圖9入水前后結構體形變可知:仿真結果顯示形變主要發生在結構體的頭部,在錐形段部分沒有發生相關形變。

外場試驗時,測試彈丸以350 m/s的速度進行入水沖擊試驗,試驗前后入水結構體變形情況如圖10所示。

試驗結果表明:1) 結構體入水的形變主要發生在結構體的頭部,這與仿真結果保持一致;2) 外場試驗實測產品產生3.10 mm的變形,仿真計算入水結構體產生3.231 1 mm變形,仿真計算與外場試驗結果基本一致。

圖9 入水后結構體形變云圖Fig.9 Deformation of structure after water entry

圖10 入水前后結構體變形(試驗結果)Fig.10 Variation of projectile acceleration during water entry with structure

驗證了沖擊修正的任意拉格朗日-歐拉(C-ALE)方法對錐形殼體入水沖擊計算準確,可以用于指導結構體高速入水的仿真。

3 適應性與數據分析

C-ALE算法與ALE算法采用同樣的網格模型,因此可以通過ALE算法進一步驗證C-ALE算法對結構體入水過程仿真的適用性。

3.1 C-ALE算法與結構體形變收斂特性

由圖11和表3可知,采用同一網格模型,在相同輸入條件下,通過ALE算法進行計算。ALE算法計算收斂時間約為8 976 s,收斂后,入水結構體的結構形變為3.76 mm。C-ALE算法計算收斂時間約為3 342 s,入水結構體變形3.230 9 mm。相對誤差從21.2%減小到4.2%,絕對誤差由0.66 mm減小到0.13 mm,收斂時間8 976 s減小到3 342 s。可以得出:針對結構體入水計算過程,相同輸入條件下,C-ALE算法比ALE算法具有更小的計算誤差和更快的收斂速度。

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

表3 不同計算方法的結構體入水計算結果Tab.3 Calculation results of water inflow of structure with different calculation

3.2 入水深度與計算時間的關系

為了進一步驗證C-ALE算法對于結構體入水深度的計算時間,在350 m/s入水速度條件下,采用同樣的網格和材料模型,對比ALE算法對結構體入水50 mm深度的仿真過程進行監測,獲得C-ALE和ALE算法入水深度與計算時間關系如圖12所示。

圖12 C-ALE和ALE算法入水深度與計算時間關系Fig.12 C-ALE和ALE algorithm relationship between water depth and calculation time

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

3.3 入水速度與計算時間的關系

為確定不同入水速度對于算法仿真速度的影響,入水速度分別取50、100、150、200、250、300、350 m/s,采用同樣的網格和材料模型,對比ALE算法對結構體入水50 mm深度的仿真過程進行監測,獲得不同入水速度C-ALE和ALE下入水深度與計算時間如圖13所示。

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

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

表4 50 mm深度不同入水速度下結構體形變與計算時間Tab.4 Structure deformation and calculation time with different water entry velocity at 50mm depth

根據表4可知,入水速度在150 m/s時候,結構體開始產生形變,隨著入水速度的增加,結構體的形變增加,C-ALE算法中引入貫穿量C對于仿真的影響增大,對比ALE算法,進一步縮短了仿真時間。

3.4 網格數量與收斂性

進一步研究不同網格數量對C-ALE算法在入水仿真過程的收斂特性的影響,分別對如下網格特性進行仿真計算。A網格特性:流體域網格數量675,固體域網格數量373;B網格特性:流體域網格數量2 325,固體域網格數量1 562;C網格特性:流體域網格數量6 566,固體域網格數量2 987。不同網格特性的收斂過程如圖14所示。仿真條件為350 m/s入水沖擊過程,實際形變為沖擊試驗形變。C-ALE算法下不同網格特性下結構體形變收斂時間如表5所示。

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

表5 不同網格特性下結構體形變與收斂時間Tab.5 Structure deformation and calculation time with different mesh methods

由圖14和表5可知,不同流體域與固體域的網格數量僅影響仿真的收斂速度,對仿真的準確度沒有影響。

4 結論

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