周 蕊 李 理 田保林 ,?,

* (北京應(yīng)用物理與計算數(shù)學研究所,北京 100094)

? (北京大學應(yīng)用物理與技術(shù)研究中心,北京 100091)

引言

爆轟波是包含快速化學反應(yīng)的強間斷面.凝聚炸藥爆轟過程釋放的巨大壓力或能量可以使周圍介質(zhì)變形、屈服,甚至產(chǎn)生巨大的損傷效應(yīng)[1-2].因此,凝聚炸藥爆轟在惰性材料約束下的傳播發(fā)展過程,及其驅(qū)動效應(yīng)一直是工程應(yīng)用中被廣泛研究的關(guān)鍵問題之一[3-5].目前,對于惰性金屬材料與凝聚炸藥爆轟的相互作用機制認識還存在不足,尤其對于不同約束條件下非理想爆轟過程的規(guī)律性認識還有待系統(tǒng)深入地研究[6].

考慮到實驗成本和觀測手段的局限性,多年來數(shù)值模擬一直是研究爆轟驅(qū)動多介質(zhì)問題的重要手段.這一問題是典型的爆轟彈塑性多介質(zhì)大變形問題,模擬該問題的數(shù)值方法大概可以分為兩類.一類是歐拉框架下的高精度數(shù)值方法[7-12],這類方法的優(yōu)點是網(wǎng)格固定不動,對大變形問題的適應(yīng)性強,但存在無法清晰捕捉物質(zhì)界面,界面附近存在數(shù)值耗散等不足.另一類是Lagrange 框架下的數(shù)值方法[6,13-14],這類方法因網(wǎng)格跟隨流體自適應(yīng)運動,可以自然地捕捉物質(zhì)界面或自由面.特別地,Lagrange 框架下適應(yīng)多分區(qū)的考慮接觸滑移的交錯型數(shù)值方法在模擬彈塑性、多介質(zhì)、薄殼和大空腔等復(fù)雜問題時具有獨特的優(yōu)勢[15-16],成為工程應(yīng)用中爆轟驅(qū)動多介質(zhì)問題不可或缺的數(shù)值模擬手段[17-19].

為了揭示炸藥爆轟和惰性金屬介質(zhì)相互作用時內(nèi)部的細致流場結(jié)構(gòu)和演化過程,需要在相當密的網(wǎng)格下開展數(shù)值模擬研究.然而,Lagrange 方法在較小的網(wǎng)格尺度下對多介質(zhì)大變形問題進行數(shù)值模擬時,存在因網(wǎng)格畸變導致計算不健壯的問題.在上述問題中,我們主要關(guān)注的是爆轟反應(yīng)區(qū)附近及金屬材料中沖擊波附近區(qū)域的流場結(jié)構(gòu)和相互作用規(guī)律,為了避免整體密網(wǎng)格Lagrange 計算帶來的效率和健壯性困難,發(fā)展適應(yīng)凝聚炸藥爆轟驅(qū)動多介質(zhì)問題的Lagrange 網(wǎng)格自適應(yīng)方法(adaptive mesh refinement,AMR)是十分必要的.這可以在保障我們關(guān)心區(qū)域網(wǎng)格分辨率的同時,粗化爆轟波或激波后膨脹區(qū)的網(wǎng)格尺度,提升波后大變形區(qū)域Lagrange計算的健壯性.

已有AMR 方法多是基于歐拉框架[20],采用分層的結(jié)構(gòu)網(wǎng)格塊自適應(yīng)加密所關(guān)心的流場區(qū)域.在歐拉框架下因網(wǎng)格固定不動,不同層網(wǎng)格分別計算后,一般通過時間自適應(yīng)處理不同層網(wǎng)格邊界之間的銜接,這方面已發(fā)展了相對成熟的數(shù)值算法.在Lagrange 框架下的AMR 方法研究并不多見,Anderson等[21]曾將歐拉框架下結(jié)構(gòu)網(wǎng)格AMR 方法拓展到Lagrange 框架下,發(fā)展了相適應(yīng)的層間耦合算法.Wang[22]發(fā)展了適應(yīng)復(fù)雜材料響應(yīng)過程的ALE-AMR方法.王瑞利等[23]基于多介質(zhì)流體力學程序LAD2D發(fā)展了非結(jié)構(gòu)網(wǎng)格Lagrange 自適應(yīng)方法.然而,以上工作都未見在復(fù)雜幾何邊界、多分區(qū)和多介質(zhì)問題上的應(yīng)用.Morrell 等[24-25]基于CORVUS 程序發(fā)展了基于cell 的非結(jié)構(gòu)網(wǎng)格自適應(yīng)方法,實現(xiàn)了跨不同物質(zhì)區(qū)的自適應(yīng)模擬.然而,他們采用基于“連通數(shù)組”的數(shù)據(jù)結(jié)構(gòu),這類數(shù)據(jù)結(jié)構(gòu)與程序緊耦合,不易推廣,且隨著加密層數(shù)的增多,程序?qū)崿F(xiàn)的復(fù)雜程度越來越大.

本文針對凝聚炸藥爆轟驅(qū)動多介質(zhì)問題精密化模擬需求,提出了一種新的Lagrange 非結(jié)構(gòu)網(wǎng)格多層自適應(yīng)方法,及其與接觸滑移的自適應(yīng)耦合算法,實現(xiàn)了彈塑性、多介質(zhì)和多分區(qū)問題的Lagrange 自適應(yīng)數(shù)值模擬能力.該方法既發(fā)揮了非結(jié)構(gòu)網(wǎng)格分層數(shù)據(jù)結(jié)構(gòu)的靈活性優(yōu)勢,也通過將有效網(wǎng)格壓至一層進行Lagrange 計算,避免了傳統(tǒng)AMR 中不同層級網(wǎng)格分別計算帶來的層間耦合困難.采用該方法開展了拐角爆轟、多層炸藥隔爆和有限尺寸彎道爆轟等復(fù)雜問題的數(shù)值模擬研究.計算結(jié)果顯示,在保障爆轟波、激波附近網(wǎng)格分辨率的前提下,該模擬方法可大幅縮減可計算網(wǎng)格規(guī)模,提升Lagrange計算的健壯性,獲得了豐富的凝聚炸藥爆轟與惰性金屬材料相互作用的流場信息,可為揭示其中的相互作用規(guī)律提供豐富的定量數(shù)據(jù).

1 適應(yīng)多分區(qū)的Lagrange 非結(jié)構(gòu)網(wǎng)格多層自適應(yīng)方法

本文基于二維爆轟彈塑性流體力學Lagrange 程序發(fā)展非結(jié)構(gòu)網(wǎng)格多層自適應(yīng)方法,基礎(chǔ)算法的具體細節(jié)可參考文獻[13,26].該程序的主要控制方程是考慮彈塑性[27]和復(fù)雜狀態(tài)方程的流體力學守恒方程組,具體形式如下

其中 β=1 為平面問題,β=0 為以z方向為對稱軸的軸對稱問題.ρ 和e分別代表密度和比內(nèi)能.vr和vz分別是r方向和z方向的速度.σrr,σzz,σθθ和σrz是Cauchy 應(yīng)力張量.τ 是偏應(yīng)力張量,應(yīng)力與偏應(yīng)力滿足如下關(guān)系

式中 τ 由以下本夠方程求得,壓力p由狀態(tài)方程求得.

為了使方程系統(tǒng)封閉,需要引入不同材料的狀態(tài)方程,本文凝聚炸藥采用JWL 狀態(tài)方程[28].其中未反應(yīng)炸藥的狀態(tài)方程形式如下

炸藥爆轟反應(yīng)產(chǎn)物的狀態(tài)方程形式如下

表1 炸藥狀態(tài)方程參數(shù)取值Table 1 The equation of state parameters for high-explosive

凝聚炸藥爆轟反應(yīng)過程用Lee-Tarver 三項式反應(yīng)率模型,具體形式如下

其中 λ 為炸藥的反應(yīng)份額,λ=0 時代表未反應(yīng)炸藥,λ=1時代表爆轟產(chǎn)物.式中的I,b,a,x,G1,c,d,y,G2,e,g,z為反應(yīng)率模型參數(shù).本文計算所用的炸藥反應(yīng)率模型參數(shù)取值如表2 所示.

表2 炸藥Lee-Tarver 三項式反應(yīng)率模型參數(shù)取值Table 2 The Lee-Tarver parameters for high-explosive

其中e為比內(nèi)能,a0,s,Γ0為狀態(tài)方程參數(shù).本文數(shù)值算例涉及的金屬材料包括鋼和鋁,它們的狀態(tài)方程參數(shù)取值如表3 所示.

表3 金屬材料 M ie-Grneisen 狀態(tài)方程參數(shù)取值Table 3 The M ie-Grneisen equation of state for steel and aluminium

表3 金屬材料 M ie-Grneisen 狀態(tài)方程參數(shù)取值Table 3 The M ie-Grneisen equation of state for steel and aluminium

上述控制方程的離散基于交錯型Lagrange 格式,詳細的離散過程可參考文獻[26],篇幅原因,本文不再贅述.所有的熱力學量定義在單元中心上,例如 ρ,τ,p,e等;速度v定義在網(wǎng)格節(jié)點上.在Lagrange計算中,每一個單元的質(zhì)量是守恒的,單元密度由初始質(zhì)量除以當前體積求得.動量方程采用基于面格式的有限元方法離散,內(nèi)能方程基于有限體積格式離散.在上述方法體系下,需要引入人為黏性捕捉激波間斷,使得人為黏性耗散產(chǎn)生的熵增等于激波引起的物理熵增,以此來達到正確模擬激波波前和波后狀態(tài)關(guān)系的目的.本文采用文獻[26]中人為黏性形式,具體形式如下

其中,二維問題中l(wèi)=為網(wǎng)格的特征尺度,A為網(wǎng)格面積,為應(yīng)變率張量,a為聲速.C0和C1為無量綱參數(shù),在本文的算例中取值分別為1.5 和0.06.

為了適應(yīng)多分區(qū)多介質(zhì)問題,采用接觸滑移算法處理不同物質(zhì)之間的相互作用.本文采用基于分配參數(shù)思想的接觸滑移算法[29-30],它可以處理接觸后兩接觸面具有相對切向滑移、由分開到接觸、由接觸到分開等復(fù)雜情況的問題.這類方法在我們關(guān)心的多介質(zhì)大變形相關(guān)工程問題中被廣泛應(yīng)用.

自適應(yīng)數(shù)據(jù)結(jié)構(gòu)的構(gòu)建和層間耦合策略是自適應(yīng)方法的核心.此外,考慮到我們關(guān)心的問題是多介質(zhì)、多分區(qū)的,網(wǎng)格自適應(yīng)與接觸滑移的耦合算法也是本研究需要解決的關(guān)鍵問題之一.下邊將分別介紹這些內(nèi)容.

1.1 非結(jié)構(gòu)數(shù)據(jù)結(jié)構(gòu)及加密準則

非結(jié)構(gòu)多層數(shù)據(jù)結(jié)構(gòu)是實現(xiàn)Lagrange 多層AMR 的重要基礎(chǔ).在這個數(shù)據(jù)結(jié)構(gòu)中,我們設(shè)計了“單元類”、“邊類”、“節(jié)點類”等自定義數(shù)據(jù)類型.

(1) “單元類”中提供單元相關(guān)的幾何信息、物理量信息、自適應(yīng)信息,如表4 所示.具體如下.

表4 “單元類”中定義的變量及其物理含義Table 4 The variables and their meanings in “cell type”

①單元幾何信息:單元在當前網(wǎng)格層上的非結(jié)構(gòu)標識號cell_ID;組成單元的節(jié)點在當前網(wǎng)格層上的標識號cell_point(1:N_Point),其中N_point 為組成單元的節(jié)點數(shù),本文取值為4;組成單元的邊在當前網(wǎng)格層上的標識號cell_edge(1:N_edge),其中N_edge 為組成單元的邊數(shù),本文取值為4;共邊單元在當前網(wǎng)格層上的標識號sameedge_cell(1:N_edge).

②單元物理量信息:材料標識mat;自適應(yīng)加密準則用到的單元量,例如壓力pressure、密度density、壓力梯度P_grad、密度梯度ρ_grad、爆轟反應(yīng)份額fac 等.

③單元自適應(yīng)信息:在當前網(wǎng)格層上單元即將被加密的層數(shù)level;父單元在其所屬網(wǎng)格層上的局部單元標識號father_cell;子單元在下一網(wǎng)格層上的局部單元標識號son_cell(1:4),如果當前網(wǎng)格層上單元被加密,那么相應(yīng)單元將一分為4,形成4 個小單元加入到下一網(wǎng)格層中.

(2) “邊類”中提供與邊相關(guān)的幾何信息和自適應(yīng)信息,如表5 所示.具體如下.

表5 “邊類”中定義的變量及其物理含義Table 5 The variables and their meanings in “edge type”

①邊幾何信息:邊在當前網(wǎng)格層上的非結(jié)構(gòu)標識號edge_ID;組成邊的節(jié)點在當前網(wǎng)格層上的標識號edge_point(1:N_point),其中N_point 為組成邊的節(jié)點數(shù),本文取值為2;邊兩側(cè)單元在當前網(wǎng)格層上的標識號edge_cell(1:N_cell),其中N_cell 為邊兩側(cè)的單元數(shù),如果邊是計算的邊界,那么N_cell 取值為1,否則取值為2.

②邊自適應(yīng)信息:父邊在其所屬網(wǎng)格層上的局部標識號father_edge;如果邊所在的單元被加密,子邊在下一層網(wǎng)格上的局部標識號son_edge(1:2).

(3) “節(jié)點類”中提供與網(wǎng)格節(jié)點相關(guān)的幾何信息、物理量信息和自適應(yīng)信息,如表6 所示.具體如下.

表6 “節(jié)點類”中定義的變量及其物理含義Table 6 The variables and their meanings in “node type”

①節(jié)點幾何信息:節(jié)點在當前網(wǎng)格層上的非結(jié)構(gòu)標識號node_ID;節(jié)點周圍的單元在當前網(wǎng)格層上的標識號node_cell(1:N_cell),其中N_cell 為節(jié)點周圍的單元數(shù)量,取值在1～4 之間;節(jié)點相連的邊在當前網(wǎng)格層上的標識號node_edge(1:N_edge),其中N_edge 為與節(jié)點相連的邊的個數(shù),取值為2 或4.

②節(jié)點物理量信息:節(jié)點位置node_r,node_z;節(jié)點邊界條件標識node_boundarycondition,存在3 種類型,0 代表自由邊界,1 代表r方向約束,2 代表z方向約束,3 代表固定不動,4 代表滑移邊界.

③節(jié)點自適應(yīng)信息:對應(yīng)上一網(wǎng)格層上的節(jié)點局部標識號father_node;對應(yīng)下一網(wǎng)格層上的節(jié)點局部標識號son_node.如果當前層節(jié)點在上一層或下一層網(wǎng)格上沒有對應(yīng)的節(jié)點,那么對應(yīng)節(jié)點標識號為0.

以上數(shù)據(jù)結(jié)構(gòu)是普適且獨立的,它不依賴于程序底層和數(shù)值算法.根據(jù)物理需求,我們首先確定基礎(chǔ)層網(wǎng)格,它可以是整個計算域的網(wǎng)格,也可以是其中的一部分.將基礎(chǔ)層網(wǎng)格信息裝配到上述非結(jié)構(gòu)數(shù)據(jù)結(jié)構(gòu)中,根據(jù)物理過程特點,確定每個基礎(chǔ)層單元的加密層數(shù).對于凝聚炸藥區(qū),我們一般采用反應(yīng)份額作為網(wǎng)格加密判據(jù);對于金屬材料區(qū),我們采用密度梯度或壓力梯度作為網(wǎng)格加密判據(jù).首先,我們設(shè)置最高加密層數(shù)為L,通過加密判據(jù)先確定基礎(chǔ)層網(wǎng)格上加密層數(shù)為L的單元.舉例來說,對于炸藥區(qū),我們將反應(yīng)份額大于0.01,小于0.99 的單元加密層級標識為L.隨后,通過非結(jié)構(gòu)數(shù)據(jù)結(jié)構(gòu)中的“單元類”中的共邊單元信息,確定加密L-1 層的單元.以此類推,直到加密層級為1 的單元確定完成為止.因此,所有基礎(chǔ)層網(wǎng)格上的單元加密與否,以及加密層數(shù)被確定.在上述加密策略中,相鄰基礎(chǔ)層單元的加密層級最多相差一層,這個約束條件對于后文中的層間耦合是非常重要的.

為了便于理解,我們給出一個簡單的示例來更清楚地解釋基礎(chǔ)層網(wǎng)格單元加密層級確定的過程.如圖1 所示,在一個方形平面計算域中,初始時刻在原點附近給定一個高壓,形成散心爆轟波向外傳播.在某時刻,基礎(chǔ)層網(wǎng)格上的反應(yīng)份額分布如圖1 所示,我們將反應(yīng)份額在0.01～0.99 之間的單元加密層級確定為3 層,由此確定了如圖2 左圖所示的紅色單元的加密層級為3.根據(jù)基礎(chǔ)層網(wǎng)格上紅色單元的共邊單元信息,確定了如圖2 中間圖所示的綠色單元加密層級為2.進一步地,由綠色單元的共邊單元信息,確定了如圖2 右圖所示的藍色單元加密層級為1.至此,確定了基礎(chǔ)層網(wǎng)格上所有單元的加密層級.

圖1 基礎(chǔ)層網(wǎng)格及爆轟反應(yīng)份額分布(僅顯示了反應(yīng)份額大于0.01 的單元)Fig.1 Basic coarse mesh and its reactive variable distribution

圖2 基礎(chǔ)層網(wǎng)格上單元加密層級Fig.2 Refinement levels of the basic coarse cells

在非結(jié)構(gòu)多層數(shù)據(jù)結(jié)構(gòu)中,第0 層網(wǎng)格信息與基礎(chǔ)層網(wǎng)格信息一致.按照前述過程確定了第0 層網(wǎng)格上每個單元的加密層級.如圖3 網(wǎng)格細分過程示意圖所示,假設(shè)第0 層網(wǎng)格上的4 個單元C0_1,C0_2,C0_3,C0_4 加密層級分別為3,1,2,0.對于單元C0_4 不存在子單元,不會出現(xiàn)在下一層網(wǎng)格信息中.第一層網(wǎng)格數(shù)據(jù)信息的構(gòu)建以第0 層數(shù)據(jù)信息為基礎(chǔ),在第0 層網(wǎng)格上,如果單元在當前網(wǎng)格層上的加密層級不小于1,那么這個單元被細分為4 個小單元,這4 個小單元加入到第1 層網(wǎng)格數(shù)據(jù)信息中.在第1 層網(wǎng)格上,所有被0 層網(wǎng)格細分后的小單元都有一個局部標識號,如圖3 中第一層網(wǎng)格上單元C1_1,···,C1_12.這12 個單元在當前網(wǎng)格層上的加密層級等于父單元的加密層級減1,即C1_1,···,C1_4 這4 個單元在第一層網(wǎng)格上的局部加密層級為2,C1_5,···,C1_8 這4 個單元在第一層網(wǎng)格上的局部加密層級為0,C1_9,···,C1_12 這4 個單元在第一層網(wǎng)格上的局部加密層級為1.按照這個思路,就完成了第一層網(wǎng)格非結(jié)構(gòu)數(shù)據(jù)信息的構(gòu)建.同樣地,基于第1 層網(wǎng)格的非結(jié)構(gòu)數(shù)據(jù)信息,如果單元當前加密層級不小于1,那么這個單元被細分為4 個小單元,這4 個小單元加入到第2 層數(shù)據(jù)信息中,如圖3中第2 層網(wǎng)格上的單元C2_1,···,C2_32.這32 個單元在當前網(wǎng)格層上的加密層級等于父單元的加密層級減1,即C2_1,···,C2_16 這16 個單元在第2 層網(wǎng)格上的局部加密層級為1,C2_16,···,C2_32 這16 個單元在第二層網(wǎng)格上的局部加密層級為0.這樣就構(gòu)建了第二層網(wǎng)格上的非結(jié)構(gòu)數(shù)據(jù)信息.以此類推,基于第2 層網(wǎng)格的非結(jié)構(gòu)數(shù)據(jù)信息,如果單元當前加密層級不小于1,那么這個單元被細分為4 個小單元,這4 個小單元加入到第3 層數(shù)據(jù)信息中.如圖3中第3 層網(wǎng)格上的單元C3_1,···,C3_64.

圖3 網(wǎng)格細分過程示意圖Fig.3 Schematic diagram of mesh subdivision progress

按照上述網(wǎng)格細分過程,在每一層網(wǎng)格上都形成了當前層的“單元類”、“邊類”、“節(jié)點類”信息,因此,當前層上局部單元-邊(cell_edge)、單元-節(jié)點(cell_point)、共邊單元(sameedge_cell);邊-單元(edge_cell)、邊-節(jié)點(edge_point);節(jié)點-單元(node_cell)、節(jié)點-邊(node_edge)等幾何信息被構(gòu)建.此外,相鄰網(wǎng)格層之間的連通信息也被構(gòu)建,例如單元-子單元(son_cell(1:4))、單元-父單元(father_cell)、邊-子邊(son_edge(1:2))、邊-父邊(father_edge)、點-上一層節(jié)點(father_node)、點-下一層節(jié)點(son_node)等.圖1 所示算例按照上述過程,形成了非結(jié)構(gòu)分層網(wǎng)格信息如圖4 所示.

圖4 非結(jié)構(gòu)分層網(wǎng)格信息Fig.4 Unstructured hierarchical mesh information

1.2 層間耦合方法

在傳統(tǒng)的AMR 方法中,在每一層網(wǎng)格上分別進行求解,通過相鄰層之間的邊界條件處理不同層級網(wǎng)格之間的耦合.這類方法在網(wǎng)格固定不動的歐拉框架下是相對成熟和完善的[20].在Lagrange 框架下,網(wǎng)格跟隨流體運動,不同層級網(wǎng)格如果分別進行計算,在層間邊界上會涉及復(fù)雜的時間和空間雙重自適應(yīng),實現(xiàn)起來是相當困難的.本文通過提取不同層級網(wǎng)格之間的連通信息,實現(xiàn)不同層級上有效網(wǎng)格壓至一層進行Lagrange 計算的AMR 策略,避免了Lagrange 分層計算層間耦合的困難.

具體來說,在第0 層網(wǎng)格上(圖3 第1 列),如果單元的當前加密層級為0,那么該單元加入到有效計算單元隊列中,即C0_4 單元加入到有效計算網(wǎng)格中.在第1 層網(wǎng)格上(圖3 第2 列),如果單元在當前層上的局部加密層級為0,那么該單元被加入到有效計算單元隊列中,即C1_5,···,C1_8 單元加入到有效計算網(wǎng)格中.以此類推,圖3 中第3 列的C2_17,···,C2_32 單元,以及第4 列中的C3_1,···,C3_64 單元加入到有效計算網(wǎng)格中.可見,加入擬計算非結(jié)構(gòu)網(wǎng)格中的有效單元都沒有子單元,在當前層有子單元的網(wǎng)格會被下一層子單元所覆蓋.按照上述思路,圖4所示的多層網(wǎng)格中,有效計算網(wǎng)格如圖5(a)所示,不同顏色代表不同層網(wǎng)格提取的有效網(wǎng)格單元,最終用于Lagrange 計算的整體最密非結(jié)構(gòu)網(wǎng)格構(gòu)建完成,如圖5(b)所示.通過上述過程,實現(xiàn)了多層非結(jié)構(gòu)網(wǎng)格分層存儲,有效網(wǎng)格壓至一層進行Lagrange計算的AMR 策略.

圖5 不同層級上有效網(wǎng)格被壓至一層,形成Lagrange 計算的整體非結(jié)構(gòu)網(wǎng)格Fig.5 Multi-level meshes are flattened onto the final global unstructured mesh

當不同層級網(wǎng)格被壓至一層時,懸點會出現(xiàn)在層間邊界上.在相鄰單元加密層最多相差一層的約束下,網(wǎng)格邊上最多只有1 個懸點,如圖6 所示.本文借鑒文獻[25]中的方法對懸點施加約束.假設(shè)懸點為無質(zhì)量和動量的虛擬點,小單元計算得到的懸點上的節(jié)點質(zhì)量和節(jié)點動量按比例分配給相連非懸點節(jié)點,懸點的速度和Lagrange 運動位置由相連非懸點節(jié)點的速度和位置插值得到.

圖6 懸點周圍網(wǎng)格分布示意圖Fig.6 Schematic diagram of mesh distribution around the hang point

具體來說,在圖6 中,節(jié)點nh為懸點,節(jié)點n1和n2為與nh相連的兩個規(guī)則節(jié)點.l1為節(jié)點n1與節(jié)點nh的距離,l為兩個規(guī)則節(jié)點n1與n2之間的距離.定義比例因子f由下式計算

那么,懸點nh的速度和坐標位置由規(guī)則節(jié)點n1和n2的速度和坐標位置計算得到,具體形式如下

兩個規(guī)則節(jié)點的質(zhì)量和節(jié)點力按如下公式進行修正

以上懸點計算方法在文獻[24-25]中被成功應(yīng)用,該約束方法解決了層間耦合的問題.然而,對于懸點計算方法的深入系統(tǒng)分析還有待后續(xù)進一步開展.

1.3 自適應(yīng)接觸滑移耦合方法

本文采用的接觸滑移算法需要人為事先給定主線和從線,由此得到對應(yīng)的主點、主單元、從點和從單元.具體計算流程包括:判斷主從點接觸狀態(tài)、對符合條件的主點施加接觸約束、計算相互作用后主點的加速度、計算產(chǎn)生碰撞后的主點速度、計算從點的速度和加速度、調(diào)整下一步時間步長.在Lagrange 自適應(yīng)計算的過程中,如果主單元或從單元被自適應(yīng)加密或粗化,那么滑移線的幾何拓撲勢必會發(fā)生改變.本文提出了接觸滑移與Lagrange 自適應(yīng)算法的“緊耦合”策略.具體來說,在1.1 節(jié)中構(gòu)建的分層非結(jié)構(gòu)數(shù)據(jù)結(jié)構(gòu)中,在每一層網(wǎng)格上增加滑移線的幾何信息,通過已經(jīng)構(gòu)建的不同層級網(wǎng)格之間的連通信息,即能提取出滑移線的單元、節(jié)點在相鄰層級上的連通信息.通過提取這些連通信息,可以在形成整體Lagrange 計算非結(jié)構(gòu)密網(wǎng)格的同時,自適應(yīng)形成新的滑移線幾何拓撲.

圖7 給出了一個示例,計算域中包含兩塊炸藥區(qū),初始時刻在原點附近給定一個高壓,散心爆轟波逐漸向外傳播,從下邊的炸藥區(qū)傳至上邊的炸藥區(qū).上邊的炸藥區(qū)定義為滑移線的主側(cè),下邊定義為從側(cè).計算到某一時刻,自適應(yīng)網(wǎng)格分布如圖7(a)所示,圖中藍色線為滑移線所在位置.我們將紅色矩形圍起的區(qū)域放大,如圖7(b)和圖7(c)所示.在基礎(chǔ)層網(wǎng)格上,從單元的全局單元編號依次為211,212,213,···,對應(yīng)主單元的全局單元編號為191,192,193,···,如圖7(b)所示.通過提取和搜索相鄰網(wǎng)格層上的滑移線幾何連通信息,形成新的整體非結(jié)構(gòu)密網(wǎng)格上的滑移線幾何拓撲關(guān)系.如圖7(c)所示,新的從單元全局編號依次為184,328,331,540,543,552,555,···,對應(yīng)的主單元的全局編號為173,325,326,529,530,533,534.

圖7 接觸滑移算法與Lagrange AMR 方法耦合,形成新的接觸滑移幾何拓撲Fig.7 New geometric topology of the sliding interface is rebuilt in the adaptive coupling strategy

1.4 Lagrange 自適應(yīng)計算流程

完成上述AMR 方法設(shè)計的基礎(chǔ)上,按照如圖8所示的流程實現(xiàn)上述AMR 策略與Lagrange 主體算法的耦合.首先,提取用于自適應(yīng)計算的基礎(chǔ)層網(wǎng)格,將其裝配至非結(jié)構(gòu)自適應(yīng)數(shù)據(jù)結(jié)構(gòu)中.然后,根據(jù)所關(guān)心物理問題的特點確定自適應(yīng)加密準則,基于此,確定基礎(chǔ)層網(wǎng)格上每個單元的加密層級.隨后,分層構(gòu)建非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)信息,得到每一層網(wǎng)格上所有的幾何信息、自適應(yīng)連通信息以及必要的物理量信息.最后,提取出用于Lagrange 計算的有效非結(jié)構(gòu)網(wǎng)格信息,以及懸點和滑移線信息,將它們傳回給Lagrange 主體算法開展Lagrange 計算.在Lagrange計算的過程中,需要實時監(jiān)測自適應(yīng)網(wǎng)格分布是否滿足物理過程特點,根據(jù)監(jiān)測判據(jù)開展新的自適應(yīng)信息的構(gòu)建.

圖8 AMR 與Lagrange 主體算法耦合流程Fig.8 Coupling low between the AMR and Lagrange algorithm

以前文中提到的簡單算例為例,給出自適應(yīng)網(wǎng)格更新的步驟如下.

(1) 假設(shè)自適應(yīng)計算至某一時間步,Lagrange 計算網(wǎng)格分布如圖9(a)所示,自適應(yīng)非結(jié)構(gòu)數(shù)據(jù)結(jié)構(gòu)中存儲了完備的用于計算的非結(jié)構(gòu)網(wǎng)格幾何信息和物理量信息.基于這個密網(wǎng)格,得到對應(yīng)基礎(chǔ)層網(wǎng)格上的必要物理量.基礎(chǔ)層網(wǎng)格的節(jié)點位置由自適應(yīng)網(wǎng)格上對應(yīng)節(jié)點位置直接賦值得到,物理量由相交重映計算得出.需要強調(diào)的是,這個過程中只重映與加密判據(jù)相關(guān)的物理量即可.例如,對于爆轟問題來說,只重映反應(yīng)份額即可.基礎(chǔ)層網(wǎng)格上的反應(yīng)份額分布如圖9(b)所示.

圖9 自適應(yīng)網(wǎng)格更新流程Fig.9 Update process of adaptive mesh

(2) 根據(jù)圖2 所示流程,確定基礎(chǔ)層網(wǎng)格上每個單元是否加密,以及加密層級,如圖9(c)所示.圖中紅色單元加密層級為3,綠色單元加密層級為2,藍色單元加密層級為1.

(3) 按照圖3 所示流程,分層構(gòu)建非結(jié)構(gòu)數(shù)據(jù)信息,如圖9(d)所示.由此得到每一網(wǎng)格層上局部幾何拓撲關(guān)系,以及相鄰層級之間的自適應(yīng)連通信息.

(4) 通過提取不同層級之間的連通信息,獲得有效計算非結(jié)構(gòu)密網(wǎng)格的數(shù)據(jù)信息,新的壓至一層的自適應(yīng)非結(jié)構(gòu)網(wǎng)格分布如圖9(e)所示.

(5) 通過以上步驟,得到舊自適應(yīng)網(wǎng)格(圖9(a))的幾何信息和物理量信息,以及新自適應(yīng)網(wǎng)格(圖9(e))的幾何信息.采用非結(jié)構(gòu)相交重映方法獲得新自適應(yīng)網(wǎng)格上的物理量信息.

(6) 圖9(e)所示的自適應(yīng)網(wǎng)格幾何信息和物理量信息傳回給主體程序進行后續(xù)的Lagrange 計算.隨后,根據(jù)加密判據(jù),返回步驟(1).

2 算例驗證

考慮到采用的爆轟彈塑性流體力學Lagrange程序已經(jīng)經(jīng)過多種復(fù)雜算例的驗證與確認,對凝聚炸藥爆轟過程的模擬有較高的置信度.例如對于典型PBX9404 炸藥的一維爆轟傳播過程,數(shù)值計算得到的穩(wěn)定爆轟波傳播速度為 0.8817 cm/μs,相同條件下,根據(jù)Rayleigh-Hugoniot 關(guān)系推導的理論爆轟波傳播速度為 0.8808 cm/μs,可見數(shù)值模擬結(jié)果與理論解的偏差僅為0.1%.因此,在本節(jié)的數(shù)值驗證中,一般將一致均勻密網(wǎng)格的數(shù)值模擬結(jié)果作為基準解,將Lagrange 自適應(yīng)模擬結(jié)果與其定量比較,來說明所提出自適應(yīng)方法模擬結(jié)果的正確性.

我們采用一維爆轟傳播問題考核所提出Lagrange非結(jié)構(gòu)網(wǎng)格多層自適應(yīng)方法的正確性.圖10 為一維爆轟傳播問題的數(shù)值模擬結(jié)果.初始時刻,在炸藥左端設(shè)置高壓區(qū)起爆爆轟波,在爆轟波傳播的過程中,自適應(yīng)密網(wǎng)格一直跟隨爆轟波陣面移動,使得爆轟反應(yīng)區(qū)附近一直保持預(yù)設(shè)的網(wǎng)格分辨率.將自適應(yīng)加密3 層的模擬結(jié)果與同程度均勻密網(wǎng)格模擬結(jié)果進行定量比較,如圖10 所示,相同時刻爆轟反應(yīng)區(qū)附近的網(wǎng)格分布及壓力分布完全一致.兩種方式計算得到的最大壓力均為36.98 GPa,爆轟波位置均在z=12.5567 cm 處,驗證了Lagrange 自適應(yīng)模擬方法的正確性.自適應(yīng)模擬在保證爆轟反應(yīng)區(qū)網(wǎng)格分辨率的前提下,相比于均勻一致加密,節(jié)省95%的網(wǎng)格規(guī)模.

圖10 一維爆轟傳播問題 (上:AMR(加密3 層)計算結(jié)果,下:一致密網(wǎng)格計算結(jié)果)Fig.10 1D detonation propagation (up:AMR calculation,down:uniformly fine mesh calculation)

為了定量考核Lagrange 多層自適應(yīng)方法與接觸滑移耦合算法的正確性,我們設(shè)計了如圖11 所示的一維爆轟傳播算例.這個算例包括3 個物質(zhì)區(qū).2 條滑移線、2 種材料,同時還存在滑移線相交的情況.該問題本質(zhì)上是一維問題,初始時刻,在左側(cè)金屬材料鋼區(qū)與右側(cè)兩塊炸藥之間存在空隙,需采用開穴滑移算法計算.兩塊炸藥之間的滑移線不考慮相對滑動,采用束縛滑移算法計算.初始時刻,在炸藥區(qū)左側(cè)設(shè)置高壓區(qū)起爆爆轟波,爆轟產(chǎn)物膨脹會使金屬與炸藥之間的空隙閉合.如圖11 所示,Lagrange自適應(yīng)方法能正確模擬含復(fù)雜滑移線(含相交)情況的算例.圖11(b)和圖11(c)分別為采用自適應(yīng)方法模擬和同程度均勻密網(wǎng)格模擬的壓力分布,兩種方式模擬在相同時刻的爆轟波位置和爆轟波附近壓力分布定量一致,最大壓力均為36.9548 GPa,爆轟波面位置在z=3.8349 cm 處,驗證了Lagrange 自適應(yīng)方法與接觸滑移耦合算法的正確性.

圖11 含多條滑移線的一維爆轟傳播問題,藍色實線為滑移線所在位置Fig.11 1D detonation propagation with slide line,blue solid line is slide line

更進一步地,我們設(shè)計了如圖12 所示的二維爆轟傳播與對碰問題,來考察Lagrange 多層自適應(yīng)方法對于多點起爆、爆轟波對碰等復(fù)雜波系結(jié)構(gòu)的適應(yīng)性.我們分別在兩個角點處設(shè)置局部高壓起爆爆轟波,兩道爆轟波相向傳播直至對碰.如圖13 和圖14 所示,在爆轟波傳播與對碰的過程中,自適應(yīng)密網(wǎng)格有效地跟隨激波運動.隨著兩道爆轟波靠近、對碰,兩塊非結(jié)構(gòu)自適應(yīng)密網(wǎng)格自動銜接.同時,在相同時刻自適應(yīng)加密3 層的模擬結(jié)果與同程序均勻一致密網(wǎng)格模擬的爆轟波陣面位置和波系結(jié)構(gòu)一致,體現(xiàn)了所發(fā)展Lagrange 多層自適應(yīng)方法對復(fù)雜問題的適應(yīng)性.

圖12 二維爆轟傳播與對碰問題、壓力分布 (t=3 μs)Fig.12 2D detonation propagation and collision,pressure distribution (t=3 μs)

3 爆轟驅(qū)動多介質(zhì)問題的Lagrange 自適應(yīng)模擬

本節(jié)開展拐角爆轟、多層炸藥隔爆、有限尺寸彎道爆轟等復(fù)雜爆轟驅(qū)動多介質(zhì)問題的Lagrange自適應(yīng)數(shù)值模擬研究.

3.1 拐角爆轟問題

本文所模擬的拐角爆轟問題,初始計算域為(0.0,20.0 cm)×(0.0,10.0 cm),其中金屬材料鋼的初始幾何尺寸為(0.0,8.0 cm)×(0.0,6.0 cm),剩下的區(qū)域為凝聚炸藥.基礎(chǔ)層網(wǎng)格劃分100×50 個單元,對應(yīng)網(wǎng)格尺寸為0.2 cm.在這個軸對稱問題中,初始時刻在炸藥區(qū)左側(cè)設(shè)置高壓區(qū)域,形成爆轟波遇金屬材料拐角形成復(fù)雜的爆轟波繞爆過程.在炸藥中,采用反應(yīng)份額作為加密判據(jù),當單元反應(yīng)份額介于0.01～0.99 之間時,單元加密3 層.在金屬材料中,采用密度梯度作為加密判據(jù),當密度梯度大于0.05 時,相應(yīng)單元加密3 層.自適應(yīng)加密后,最小網(wǎng)格尺寸為0.25 mm.如圖15 和圖16 所示,自適應(yīng)計算的密網(wǎng)格始終跟隨爆轟波陣面及金屬中沖擊波移動,有效模擬了爆轟波沿金屬邊界傳播,及遇金屬拐角繞爆的過程.一致密網(wǎng)格模擬得到的爆轟波傳播和繞爆過程與自適應(yīng)模擬結(jié)果一致.自適應(yīng)計算拐角處節(jié)點起跳時間為8.8143 μs,一致密網(wǎng)格計算拐角處節(jié)點起跳時間為8.8066 μs,相差小于0.1%.自適應(yīng)計算中爆轟波附近加密3 層實際計算網(wǎng)格總數(shù)不超過3 萬,而同程度一致密網(wǎng)格計算網(wǎng)格總數(shù)為32 萬.

圖16 軸對稱多介質(zhì)拐角爆轟問題模擬結(jié)果 (t=10 μs)Fig.16 Axisymmetric multi-material corner detonation problem (t=10 μs)

3.2 多層炸藥隔爆問題

在實際工程中,除了經(jīng)常出現(xiàn)爆轟繞爆這一復(fù)雜現(xiàn)象,爆轟波過金屬隔爆也是常見現(xiàn)象.為此,設(shè)計了如圖17 所示的多分區(qū)、多介質(zhì)問題,這個算例包含4 個物質(zhì)區(qū)、3 條滑移線、3 種材料,存在爆轟波傳播、過金屬隔爆等復(fù)雜物理過程.炸藥與金屬之間允許切向移動,接觸滑移計算采用單純滑移算法.兩塊炸藥的初始位置分別為(0.0,20.0 cm)×(0.0,2.0 cm) 和 (0.0,20.0 cm)×(2.4 cm,4.4 cm);兩塊炸藥之間的金屬為鋁,其初始幾何位置為(0.0,20.0 cm)×(2.0 cm,2.4 cm);炸藥外邊界處的金屬為鋼,其初始幾何位置為(0.0,20.0 cm)×(4.4 cm,4.8 cm).初始時刻,在坐標原點附近設(shè)置高壓點起爆爆轟波,爆轟波在下方第一塊炸藥中形成散心爆轟向外傳播,驅(qū)動中間層金屬材料,壓縮引爆上方第2 塊炸藥,產(chǎn)生隔爆.隨后在兩塊炸藥中形成兩個爆轟波同步向右側(cè)傳播.如圖17 和圖18 所示,我們僅在炸藥區(qū)采用反應(yīng)份額作為自適應(yīng)加密判據(jù),將爆轟波陣面附近網(wǎng)格自適應(yīng)加密3 層.在整個爆轟傳播發(fā)展的過程中,自適應(yīng)密網(wǎng)格都能有效地捕捉爆轟波所在位置,實現(xiàn)爆轟波陣面附近的高分辨率模擬.自適應(yīng)加密3 層與同程度均勻密網(wǎng)格模擬的爆轟波陣面附近壓力分布是一致的,驗證了該Lagrange 多層自適應(yīng)方法對復(fù)雜問題的適應(yīng)性.在這個問題中,自適應(yīng)模擬最大網(wǎng)格規(guī)模為1.3 萬,同程度均勻密網(wǎng)格模擬網(wǎng)格規(guī)模為12.84 萬,可見自適應(yīng)模擬節(jié)省了90%以上的計算量.

圖17 多分區(qū)炸藥爆轟隔爆問題模擬結(jié)果 (t=4 μs)Fig.17 Multi-block multi-material detonation problem (t=4 μs)

圖18 多分區(qū)炸藥爆轟隔爆問題模擬結(jié)果 (t=10 μs)Fig.18 Multi-block multi-material detonation problem (t=10 μs)

3.3 有限尺寸彎道爆轟問題

有限尺寸彎道中的爆轟問題廣泛存在于軍用裝備起爆系統(tǒng)的設(shè)計中,因炸藥截面尺寸較小,一般在毫米量級,若采用歐拉高精度方法,在彎道與炸藥之間的界面處會存在明顯的數(shù)值耗散,使得無法準確獲得爆轟虧損效應(yīng)等定量結(jié)果.針對有限尺寸彎道爆轟問題,我們提出的基于爆轟彈塑性多介質(zhì)流體力學程序的Lagrange 非結(jié)構(gòu)網(wǎng)格多層自適應(yīng)方法會比較適用.我們設(shè)計了如圖19 所示的算例,直管段長度為1 cm,彎道炸藥內(nèi)半徑為1 cm,炸藥截面厚度為0.2 cm,彎道炸藥內(nèi)、外各放置一層金屬鋼,厚度均為0.2 cm.爆轟波在直管端部起爆后,形成一維爆轟進入彎道溝槽中,因金屬材料邊界曲率的影響,凝聚炸藥爆轟與惰性金屬材料之間存在復(fù)雜的相互作用.初始時刻沿圓周方向炸藥區(qū)設(shè)置210 個網(wǎng)格,半徑方向設(shè)置10 個網(wǎng)格,內(nèi)側(cè)鋼和外側(cè)鋼網(wǎng)格規(guī)模與炸藥區(qū)一致.初始徑向網(wǎng)格尺寸為0.2 mm,計算中設(shè)置自適應(yīng)加密層數(shù)為3,因此最小網(wǎng)格尺寸為0.025 mm.在炸藥區(qū),自適應(yīng)加密判據(jù)為0.01<λ<0.99,其中 λ 為反應(yīng)份額;在金屬材料鋼中,自適應(yīng)加密判據(jù)為 ?p>3.5,p為網(wǎng)格單元的壓力.采用接觸滑移算法處理炸藥與金屬鋼之間的相互作用,在這個算例中,采用切向不做約束的單純滑移算法,實際計算中炸藥與鋼之間會存在摩擦,這方面的影響在日后的研究中會深入開展.

圖19 有限尺寸彎道爆轟問題初始幾何 (綠色:高能炸藥,藍色:金屬鋼,紅色:CJ 狀態(tài)爆轟產(chǎn)物)Fig.19 Initial geometry of detonation problem in finite-size curved pipeline (green:high-explosive,blue:steel,red:reactive product of CJ state)

圖20 給出了彎道爆轟不同時刻的流場分布.如圖20(a)所示,在直管中,形成穩(wěn)態(tài)爆轟波自下而上傳播,流場和自適應(yīng)加密網(wǎng)格沿直管中線左右對稱.需要注意的是,因側(cè)向約束的影響,爆轟波陣面在炸藥與金屬交界處略有彎曲,這個側(cè)向影響會使爆轟波傳播能量產(chǎn)生一定虧損.隨后,直管中爆轟波傳入彎道中,如圖20(b)～圖20(d)所示,靠近內(nèi)界面的爆轟波被稀疏,靠近外界面的爆轟波被壓縮,出現(xiàn)了復(fù)雜的非理想爆轟現(xiàn)象.爆轟波在彎道中傳播的過程中,密網(wǎng)格始終跟隨爆轟波陣面和金屬鋼中激波自適應(yīng)加密,波后膨脹區(qū)網(wǎng)格自適應(yīng)稀疏.在這個復(fù)雜幾何構(gòu)型中,Lagrange 自適應(yīng)模擬即能清晰地保持物質(zhì)界面,保障激波間斷處網(wǎng)格分辨率的同時,緩解均勻密網(wǎng)格計算該問題時波后流場因網(wǎng)格大變形、畸變導致的計算不健壯的問題.在這個問題中,自適應(yīng)模擬最大網(wǎng)格規(guī)模約4 萬,同程度均勻密網(wǎng)格模擬網(wǎng)格規(guī)模為40.32 萬,可見自適應(yīng)模擬節(jié)省了約90%的計算量.

圖20 有限尺寸彎道爆轟傳播問題不同時刻壓力和網(wǎng)格分布 (左:整體流場壓力分布,中:局部放大壓力分布,右:局部放大網(wǎng)格分布)Fig.20 Simulation results for detonation problem in finite-size curved pipeline (left:pressure distribution,middle:local enlarged pressure distribution,right:local enlarged mesh distribution)

4 結(jié)論

本文針對爆轟驅(qū)動多介質(zhì)問題的高分辨率健壯數(shù)值模擬需求,提出了適應(yīng)于爆轟彈塑性多介質(zhì)Lagrange 流體力學程序的非結(jié)構(gòu)網(wǎng)格多層自適應(yīng)加密與稀疏技術(shù),在發(fā)揮非結(jié)構(gòu)多層自適應(yīng)數(shù)據(jù)結(jié)構(gòu)靈活性優(yōu)勢的基礎(chǔ)上,通過將有效網(wǎng)格壓至一層進行Lagrange 計算,避免了Lagrange 框架下多層網(wǎng)格分別計算帶來的層間耦合困難.更進一步地,實現(xiàn)了Lagrange 多層AMR 與接觸滑移算法的自適應(yīng)耦合,使其對包含復(fù)雜幾何邊界的多分區(qū)、多介質(zhì)問題有很好的適應(yīng)性.采用上述AMR 方法,實現(xiàn)了拐角爆轟、多層炸藥隔爆、有限尺寸彎道爆轟等多個復(fù)雜多介質(zhì)、多分區(qū)問題的Lagrange 多層AMR 模擬,獲得了豐富的凝聚炸藥爆轟與彈塑性介質(zhì)相互作用的流場圖像.

后續(xù)工作將針對凝聚炸藥爆轟與惰性金屬材料相互作用規(guī)律開展深入系統(tǒng)的研究,包括本文涉及的復(fù)雜爆轟驅(qū)動多介質(zhì)問題的計算條件和物理模型參數(shù)的定量標定,在此基礎(chǔ)上獲得不同炸藥、不同約束材料和不同幾何構(gòu)型下爆轟約束效應(yīng)及其非理想行為的規(guī)律性認識.