空中點源爆炸波陣面參數的誤差控制求解方法

張 鏖,李 鋒,孫云厚,蘇榮華

(軍事科學院 國防工程研究院,北京 100850)

0 引言

在人防工程設計過程中確定沖擊荷載是進行結構安全性核算的重要一環,而確定空中點源爆炸波陣面參數是后續確定沖擊波傳播和結構荷載的初始計算步,其計算的精確性嚴重影響后續沖擊波傳播及工程結構爆炸荷載的核算,因此能夠獲取更加精確的自模擬解對工程防護結構抗力優化意義重大。空氣中球面對稱裝藥強爆炸存在移動波陣面邊界,該波陣面邊界和其自身控制方程相關聯,該問題的求解通常稱為自模擬解。Taylor[1],Von Neumann[2],Sedov[3]及國內科研人員[4]分別給出了空氣中強爆炸的自模擬解。盡管形式上有所不同,但基本原理和過程是一致的。首先構建常微分形式的無量綱的沖擊波壓力、質點速度、密度與比例距離的關系式,該常微分的繼續求解則存在多種可能的形式:精確解析解、近似解析解、數值解。其中精確解析解最終并沒有給出顯式的最終解[3],近似解近似過程則仍然存在誤差且難以求解[1],目前通常的求解辦法是采用數值解通過經典四階Runge-Kutta方法(RK4)進行求解[5,6],該求解方法可以將波陣面參數的計算誤差控制在四階。現在數值計算軟件如Autodyn中將球對稱自模擬解計算結果映射到三維網格上作為沖擊波初始參數并傳播到結構產生后續荷載[7],這種算法對自模擬解精度提出了更高要求。本文中改進原有算法,采用誤差控制的自適應步長Runge-Kutta法,以實現了對空中點源爆炸自模擬過程更精確求解。

1 空中球面對稱裝藥強爆炸及其控制方程

本文中以空中點源爆炸波陣面參數為主要研究對象,該過程可由自模擬解求得沖擊波陣面上的壓力,質點速度和密度3個重要指標。空中點源爆炸的示意圖如圖1所示,這一過程主要基于幾點假設:① 爆炸物相對于求解范圍足夠小(如核爆炸);② 空氣為理想氣體;③ 在這一過程中與外界不發生熱交換。

圖1 空中點源爆炸示意圖

在以上假設條件下可建立空中點源爆炸自模擬過程的基本方程,在球對稱條件下原來的三維球坐標問題則可以縮減為一維問題,另外根據自模擬過程的定義,知道所有變量均與一個變量相關,因此假定一個基本參量——波陣面半徑rs,而波陣面質點速度可寫為

(1)

式中:D為波陣面質點速度;rs為波陣面半徑;t為時間。根據Rankine-Hugoniot條件[8-11],沖擊波陣面上存在動量守恒、質量守恒和能量守恒

(2)

式中:ρs為波陣面上的空氣密度;Ps為波陣面上的沖擊波壓力;us為波陣面上的質點速度;γ0空氣的沖擊絕熱系數;ρ0為常溫常壓下空氣密度。

球對稱坐標下,理想氣體的運動方程滿足動量守恒、質量守恒和能量守恒方程[12]

(3)

式中:u為質點密度;p為壓力;ρ為密度;r為球坐標位置。式(2)、式(3)即是求解球面對稱裝藥強爆炸自模擬解的控制方程,通過式(2)、式(3)求解對應的波陣面上質點速度、壓力、密度等參數的過程即為自模擬解的求解過程。為求解這些待求未知量,需要對式(2)、式(3)進行無量綱化,并應用邊界條件,然后進行差分求解。

2 控制方程的無量綱化、微分形式的自模擬解及其邊界條件

根據量綱分析中的Π定理[13],比例距離和主要變量之間存在以下聯系

(4)

式中:R為無量綱化的比例距離;E0為爆心的初始能量;ρ0為爆心的初始密度。能量、密度、時間表示成可量測的變量質量M,時間T以及長度L

L∝MαL2αT-2αTβMδL-3δ

(5)

式中:α、β、δ為常量。兩邊量綱相同,因此:

α+δ=0,2α-3δ=1,-2α+β=0

(6)

求得:

α=1/5,β=2/5,δ=-1/5

(7)

根據式(4)中的比例距離和主要變量關系以及式(5)—式(7)計算得到的量綱系數,構建得到自模擬解的獨立變量無量綱的比例距離η,表示為

(8)

以無量綱的比例距離η為基本變量,根據Rankine-Hugoniot條件,待求解變量沖擊波壓力p、氣體密度u、密度ρ,分別滿足如下關系:

(9)

u=AR-3/2φ(η)

(10)

ρ=ρ0ψ(η)

(11)

將式(9)—式(11)代入式(3)進行變換(此處省略簡化公式的詳細步驟),得到待求解無量綱的壓力f,質點速度φ及密度ψ微分形式的表達式[2,6,7]

(12)

(13)

(14)

其中;γ為絕熱系數,取1.4[2]。波陣面邊界上存在邊界條件:

(15)

另外在球心對稱點處為靜止狀態滿足質點速度為0,因此φ(0)=0。

常微分方程組(式(12)—式(14))通常采用經典四階Runge-Kutta方法(RK4)進行求解[5-6],該求解方法具備4階精度。

3 自適應Runge-Kutta方法

為進一步提高計算精度,滿足更高精度的工程計算要求,本研究中采用可控精度的5階自適應Runge-Kutta方法進行求解,在滿足5階精度的情況下,可以有效減少迭代步數,同時因為誤差控制的引入,可以顯式的控制結果的誤差,從而實現人為可控的高精度求解的要求。對于任意常微分方程存在:

g′(ξ)=G(ξ,g(ξ),g(ξ0))=g0

(16)

任意函數經泰勒展開后可表示為

(17)

式中:bi為Runge-Kutta方法中的增量步的權重系數;ki為增量步i上的導數值,Runge-Kutta方法可以給出不同的增量步數和權重值對常微分方程進行求解,其中比較典型的Runge-Kutta求解方法有前向歐拉、后向歐拉、Newmark方法以及四階Runge-Kutta方法[14]。不同迭代方法的當前步第i步的導數ki均可表示為

(18)

其中:

(19)

泰勒展開式誤差可以用不同階數的展開形式求差獲得,因此對于自適應Runge-Kutta方法的誤差通過比較2種不同階數p、p-1得到[15]

(20)

為進一步提高原有四階Runge-Kutta方法沖擊波自模擬解的計算精度,采用四階-五階自適應的Runge-Kutta-Fehlberg方法(簡稱RKF方法)[16],p=5階和p-1=4階2種情況下的Butcher tableau(Runge-Kutta方法中的系數矩陣a通常通過Butcher tableau的形式給出[17]),RKF的Butcher tableau具體形式如式(22)。該方法因具有五階精度且能夠根據當前誤差自動控制步長,目前在微分方程和偏微分方程的求解中得到了廣泛應用[18-19]。

(21)

(22)

該方法的計算步長可通過如下公式獲得:

(23)

RKF方法與原有的四階Runge-Kutta方法求解自模擬解程序實現方法上的主要不同[5-6]之處在于:① Butcher tableau數據不同;② 計算完各增量點出的導數k之后,需要采用不同的權重和分別計算最終函數值,其中以五階權重對應的結果為準;③ 計算完成進行誤差的校核;④ 評估誤差是否滿足精度要求;⑤ 如不滿足修正增量步長重新計算,如滿足以當前步長作為輸入值,進入下一步計算。通過Matlab程序實現了該算法,計算程序設計流程見圖2所示。

圖2 RKF算法計算程序流程

4 計算結果與分析

沖擊波的自模擬問題沒有顯式的精確解析解,為進行對照分析,本研究中同時采用四階Runge-Kutta方法,取極小步長Δη=0.000 1對方程進行求解,作為對照分析的基準結果。然后采用固定步長RKF算法對自模擬方程進行求解,并分析了各主要爆炸波待求參量f、φ、ψ求解誤差隨比例距離η變化的分布情況。最后利用自適應步長的RKF算法對自模擬方程進行求解,得到了誤差控制的自模擬解。

4.1 RKF算法固定步長條件下的誤差分布

分析等步長計算時的誤差分布可為后續驗證RKF方法是否能在誤差較大區域自動縮減步長提供驗證性依據。因此本文中首先采用較長的步長Δη=0.005進行等步長的RKF算法計算,分析固定步長條件下的誤差分布。圖3為固定步長求解結果,從其中可以看出,3個無量綱參數均從初始值η=1.0開始迭代,隨著η的遞減,誤差逐步降低,接近數值畸點0時,誤差開始陡增。從誤差數值上看,無量綱壓力f的誤差最小,其次是無量綱密度ψ,無量綱速度φ再次之。然而數量上無量綱速度φ接近畸點是陡增至0.8,為解決該問題則需要大幅度增加計算步數。在其他數值精度均已相對可靠的情況下,因部分數值的誤差需要計算系統的額外開銷。因此,總結圖3所示的規律可以發現,誤差主要來自于自模擬過程的無量綱速度φ,為解決接近畸點處的精度問題,需要額外的計算步數。這一結果符合對算法在強非線性段加密步長的預期。

圖3 固定步長求解結果

4.2 自適應步長RKF算法求解精度分析

為解決4.1小節中無量綱速度φ接近畸點處誤差較大的問題,通過采用自適應步長的RKF算法進行求解,理論上能夠在畸點附近開始自動調節計算補充。計算中取αsafety=0.8,最大增加幅度αmax=0.4,最小縮減幅度αmin=0.125。RKF算法為五階精度,因此計算步長是p=5。盡管RKF算法能夠實現誤差1.0×10-12以下的高精度計算,但是為能夠在圖上清楚的顯示步長的變化,設置了3個較大的誤差控制參數1.0×10-5、1.0×10-6、1.0×10-7,以便誤差控制過程現實的更加明顯,計算結果如圖4所示。圖3顯示了兩頭的誤差較大,而圖4中兩頭的數據點比較密集,說明自適應步長RKF算法實現了根據誤差對計算步長的自動控制。

4.3 不同算法的計算精度和效率分析

自適應步長RKF方法除了對誤差顯式控制上存在優勢,在計算精度和計算效率方面是否依然存在優勢還需要進一步的驗證,因此本文中分別用RK4、固定步長RKF以及自動步長RKF三種方法進行了計算,其中自動步長RKF的誤差控制參數同第4.2節。計算中采用的計算硬件平臺為i7-7700k和64 G內存,軟件平臺為Windows 10 和Matlab 2019,不同算法的計算控制參數、計算時間及最大絕對誤差如表1所示。

表1 不同算法的計算精度與計算時間

從表1中可以看出RK4用時最省,為0.03 s,自動步長RKF所需時間最長為0.08 s,RKF算法在計算效率上略低于RK4算法,但是所有算法的計算時長均在0.1 s以內,因此RKF的計算效率的降低仍在可接受范圍之內。對比固定步長RKF和自動步長RKF的計算效率可以發現,即使是自動步長RKF方法在最大誤差小于固定步長RKF方法時,用時仍然小于固定步長RKF方法,因此自動步長RKF方法對比固定步長RKF方法計算效率上具有明顯優勢。分析不同計算方法的最大絕對誤差,可以看出同計算步長的RKF方法的計算誤差明顯小于RK4方法,而自動步長RKF方法則可以很好地將誤差控制在設置范圍內。

5 結論

空中點源爆炸波陣面參數的誤差會在累計到后續的裝備毀傷計算中,因此提高其精度對提升毀傷分析的可靠性有重要意義。本研究中改進了空中點源爆炸波陣面參數的求解方法,將原常用的四階Runge-Kutta方法改進為五階精度的自適應步長Runge-Kutta方法,并對這2種方法的結果進行了對比。結果表明,這一改進確能將計算精度提高了一階,且解決了原先的誤差不能顯式控制的問題。改進后的方法在計算效率上略有降低,但考慮到總體計算時長較短,總體效率上可以接受。沖擊波陣面參數求解中引入自適應步長RKF算法可為沖擊波參數的精密計算提供了理論方法和技術支持,為后續裝備毀傷評估的可靠性提供有力保證。