寧建國(guó) 原新鵬 馬天寶李 健

(北京理工大學(xué)爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京100081)

-生物、工程及交叉力學(xué)

計(jì)算動(dòng)力學(xué)中的偽弧長(zhǎng)方法研究1)

寧建國(guó) 原新鵬 馬天寶2)李 健

(北京理工大學(xué)爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京100081)

動(dòng)力學(xué)問題通常采用微分方程來描繪,但由于工程實(shí)際問題的復(fù)雜性,微分方程模型常伴隨著解的不連續(xù)性、剛性或激波間斷奇異性特點(diǎn),傳統(tǒng)方法很難求解,奇異性問題是計(jì)算動(dòng)力學(xué)難點(diǎn),同時(shí)也是國(guó)內(nèi)外學(xué)者研究的熱點(diǎn).偽弧長(zhǎng)數(shù)值算法是針對(duì)計(jì)算動(dòng)力學(xué)中的奇異性問題所提出的,其基本思想為通過在解曲線上引入偽弧長(zhǎng)參數(shù),并增加一個(gè)約束方程,在偽弧長(zhǎng)參數(shù)作用下,使得原始離散單元發(fā)生扭曲形變,從而達(dá)到消除或減弱奇異性的目的.本文首先介紹偽弧長(zhǎng)方法求解定常對(duì)流--擴(kuò)散方程的奇異性問題,并提出針對(duì)雙曲守恒定律的局部偽弧長(zhǎng)算法,其思想在于首先通過間斷解的梯度變換來確定強(qiáng)間斷所處位置,進(jìn)而通過局部網(wǎng)格點(diǎn)重構(gòu)以及數(shù)值修正來達(dá)到強(qiáng)間斷處奇異性消除與降低的目的.針對(duì)高維問題,提出全局偽弧長(zhǎng)方法,通過對(duì)整個(gè)計(jì)算區(qū)域內(nèi)的網(wǎng)格點(diǎn)進(jìn)行重構(gòu),使得所有網(wǎng)格點(diǎn)向奇異間斷點(diǎn)處移動(dòng),從而降低間斷點(diǎn)的影響域,達(dá)到降低奇異性的目的.重點(diǎn)討論了三維全局偽弧長(zhǎng)算法問題的計(jì)算難點(diǎn),即三維空間網(wǎng)格扭曲大變形導(dǎo)致的數(shù)值算法不收斂，并提出在算法設(shè)計(jì)過程中采用分塊重構(gòu)與整體計(jì)算相結(jié)合的策略，實(shí)現(xiàn)了三維空間中的偽弧長(zhǎng)數(shù)值算法,最后通過數(shù)值實(shí)驗(yàn)來驗(yàn)證偽弧長(zhǎng)算法對(duì)于奇異性問題的有效性.

計(jì)算動(dòng)力學(xué),偽弧長(zhǎng),自適應(yīng),爆炸沖擊波,三維

引言

計(jì)算動(dòng)力學(xué)是利用科學(xué)計(jì)算方法來研究動(dòng)力學(xué)現(xiàn)象與特性的科學(xué)[1].但由于動(dòng)力學(xué)問題較復(fù)雜,許多實(shí)驗(yàn)的開展面臨著周期長(zhǎng)、代價(jià)大等困難,通過理論分析又難以獲得大多數(shù)現(xiàn)實(shí)工程問題的解析解,數(shù)值方法就成為一種有效的替代手段.特別是在第二次世界大戰(zhàn)以后,電子計(jì)算機(jī)的出現(xiàn)以及隨后數(shù)十年中大型計(jì)算機(jī)的迅速普及和數(shù)值計(jì)算方法的迅猛發(fā)展[2],為計(jì)算科學(xué)的發(fā)展提供了物質(zhì)與理論基礎(chǔ).因?yàn)閯?dòng)力學(xué)問題常常伴隨著時(shí)間和空間的演化,所以計(jì)算動(dòng)力學(xué)問題通常采用基于連續(xù)介質(zhì)力學(xué)問題的微分方程模型來描述.而且由于工程實(shí)際問題的復(fù)雜性,動(dòng)力系統(tǒng)微分方程模型常伴隨著解的不連續(xù)性、剛性或激波間斷等奇異性特點(diǎn).奇異性問題是計(jì)算力學(xué)中的難點(diǎn),同時(shí)也是國(guó)內(nèi)外學(xué)者研究的熱點(diǎn)[35].1979年,Riks[6]首次應(yīng)用弧長(zhǎng)法求解了帶有參數(shù)的非線性方程中的極值問題.Boor證明了利用參數(shù)變換求解微分方程的可行性,并且提出高階逼近函數(shù)的方法[7].之后Crisfiel[8]通過變分法將非線性幾何問題的本構(gòu)方程變換為非線性代數(shù)方程,進(jìn)而引入弧長(zhǎng)參數(shù),建立了弧長(zhǎng)算法的基本理論.此外Chan[9]在1984年通過引入偽弧長(zhǎng)控制參數(shù)將其引申為偽弧長(zhǎng)方法,并將其應(yīng)用于求解非線性方程的奇異性問題中.在此基礎(chǔ)之上,忻孝康等[10]、陳國(guó)謙等[11]將偽弧長(zhǎng)方法推廣，用于求解常微分方程動(dòng)力系統(tǒng)的奇異性問題.受此啟發(fā),武際可等[12]將常微分方程動(dòng)力系統(tǒng)的偽弧長(zhǎng)算法推廣,用來求解微分動(dòng)力系統(tǒng)中的雙曲偏微分方程的Burgers奇異性問題,總結(jié)出偽弧長(zhǎng)算法的基本思想.通過在解曲線上引入偽弧長(zhǎng)參數(shù),并增加一個(gè)約束方程,在偽弧長(zhǎng)參數(shù)作用下,使得原始離散單元發(fā)生扭曲形變,達(dá)到消除或減弱奇異性的目的.此時(shí)雙曲偏微分方程的偽弧長(zhǎng)算法的框架思想已經(jīng)基本形成.與此同時(shí),出現(xiàn)了一種以網(wǎng)格自適應(yīng)為目的,通過網(wǎng)格點(diǎn)的重新分布與數(shù)值解的重新插值來求解奇異雙曲問題的數(shù)值算法——移動(dòng)網(wǎng)格方法.該方法在1983年由Harten等[13]第一次提出.Ren等[14]對(duì)該方法進(jìn)行了改進(jìn),提出利用弧長(zhǎng)參數(shù)來控制網(wǎng)格點(diǎn)的自適應(yīng)移動(dòng).在此基礎(chǔ)之上,Huang等[15]根據(jù)均分原則[16],提出了結(jié)合偽弧長(zhǎng)控制參數(shù)的移動(dòng)網(wǎng)格控制方程.此后,大量學(xué)者利用移動(dòng)網(wǎng)格控制方程與雙曲控制方程構(gòu)成的耦合系統(tǒng)來研究移動(dòng)網(wǎng)格方法[15,1718].對(duì)比分析偽弧長(zhǎng)方法與移動(dòng)網(wǎng)格方法的特點(diǎn),可以發(fā)現(xiàn),兩種數(shù)值方法具有一定聯(lián)系,移動(dòng)網(wǎng)格法就是一種具有偽弧長(zhǎng)性質(zhì)的雙曲方程數(shù)值算法,屬于計(jì)算動(dòng)力學(xué)中的偽弧長(zhǎng)數(shù)值算法的一部分.對(duì)于雙曲偏微分方程來說,雖然偽弧長(zhǎng)方法與移動(dòng)網(wǎng)格方法的出發(fā)點(diǎn)不一樣,但是移動(dòng)網(wǎng)格方法具有偽弧長(zhǎng)數(shù)值算法的特點(diǎn),符合偽弧長(zhǎng)算法的定義,可以被稱作是偽弧長(zhǎng)移動(dòng)網(wǎng)格算法.此外，由于移動(dòng)網(wǎng)格方法的出發(fā)點(diǎn)是網(wǎng)格的移動(dòng)重構(gòu),缺乏網(wǎng)格的整體性與直觀性,特別是在三維問題的計(jì)算過程中,認(rèn)為六面體網(wǎng)格在移動(dòng)后,某個(gè)面發(fā)生畸變,變形為非六面體,從而導(dǎo)致計(jì)算失敗.當(dāng)前對(duì)三維問題都是采用工程的簡(jiǎn)化計(jì)算方法來完成的,所獲得的計(jì)算結(jié)果不可避免地會(huì)失去原有工程問題的物理本質(zhì).再者大多數(shù)工程實(shí)際問題都是三維問題，所以針對(duì)三維問題的大規(guī)模高效算法的研究是迫切與必要的.

近年來,寧建國(guó)等將偽弧長(zhǎng)算法用于求解爆炸沖擊波問題,先后發(fā)展出局部偽弧長(zhǎng)算法與全局偽弧長(zhǎng)移動(dòng)網(wǎng)格算法[1924],通過偽弧長(zhǎng)變換來捕捉爆炸沖擊波陣面,建立了爆炸沖擊波問題的偽弧長(zhǎng)算法的基礎(chǔ)理論體系.先引入一般形式的Runge-Kutta法,得到雙曲守恒控制方程與網(wǎng)格映射控制方程所構(gòu)成的限制微分動(dòng)力系統(tǒng)[2526],進(jìn)一步利用牛頓迭代法,將非線性問題轉(zhuǎn)化為類線性問題來處理,繼而對(duì)類線性問題的系數(shù)矩陣進(jìn)行正則化歸約分析,建立了偽弧長(zhǎng)數(shù)值算法的穩(wěn)定性理論,得到了偽弧長(zhǎng)算法的收斂性結(jié)論[24].針對(duì)反應(yīng)氣體動(dòng)力學(xué)模型源項(xiàng)的剛性特點(diǎn)以及反應(yīng)區(qū)域的沖擊波與稀疏波的復(fù)雜作用導(dǎo)致計(jì)算過程密度、壓力等物理量出現(xiàn)負(fù)值的問題[2728],提出了偽弧長(zhǎng)數(shù)值算法的保正性理論[23].

本文在研究計(jì)算動(dòng)力學(xué)問題中的三類偽弧長(zhǎng)算法的基礎(chǔ)之上,從偽弧長(zhǎng)算法降低奇異性的特點(diǎn)出發(fā),提出了三維空間中六面體網(wǎng)格單元分割與整體結(jié)合的思想,實(shí)現(xiàn)了三維空間中爆炸沖擊問題的偽弧長(zhǎng)法的數(shù)值模擬.進(jìn)而針對(duì)二維、三維爆轟波對(duì)于板型障礙物的爆炸沖擊問題進(jìn)行數(shù)值模擬,實(shí)現(xiàn)了網(wǎng)格對(duì)于爆轟波陣面的捕捉,對(duì)障礙物后方的多個(gè)標(biāo)定點(diǎn)處的壓力大小進(jìn)行分析,得到障礙物后方的安全位置.

1 定常對(duì)流--擴(kuò)散問題中的偽弧長(zhǎng)方法

定常對(duì)流--擴(kuò)散問題廣泛存在于化學(xué)工程、傳熱傳質(zhì)、大氣和水質(zhì)污染等領(lǐng)域中[2930].其模型通常采用常微分方程描繪,在數(shù)值計(jì)算中常常會(huì)出現(xiàn)偽振蕩、非物理的負(fù)濃度解、激波層或邊界層被拉寬等奇異性現(xiàn)象.采用偽弧長(zhǎng)算法求解此類問題,可以有效地避免與降低常微分方程的奇異性問題.定常對(duì)流--擴(kuò)散方程的一維無量綱形式為

其中，u為對(duì)流速度,q為源匯項(xiàng),c為污染物的濃度,Pe稱為Peclet數(shù),定義為

式中，U為特征對(duì)流速度,L為特征長(zhǎng)度,α為擴(kuò)散系數(shù).

在大Peclet數(shù)情況下,即Pe?1,可令

方程(1)可以寫成奇異攝動(dòng)型的二階常微分方程式

其中，y=c,f,g,h均可是x和y(c)的函數(shù),α,β為常數(shù).引入解曲線的弧長(zhǎng)s,則有

令v=dy/dx,則式(4)可轉(zhuǎn)化為如下一階常微分方程組

假設(shè)x=a端的弧長(zhǎng)為零,x=b端的弧長(zhǎng)為Smax,于是式(5)化為如下的邊界條件

因此,奇異攝動(dòng)兩點(diǎn)邊值問題式(4)和式(5)就轉(zhuǎn)化為一階常微分方程的邊值問題式(7)和式(8).對(duì)于方程組(7)通常采用Rosenbrock格式求解[31].對(duì)于邊值問題式(8)可采用打靶法進(jìn)行求解[32].

2 非定常對(duì)流問題的偽弧長(zhǎng)方法

非定常對(duì)流問題通常采用雙曲型偏微分方程來刻畫[33].這類問題通常伴隨著激波間斷性而導(dǎo)致解出現(xiàn)強(qiáng)間斷奇異性[34].數(shù)值求解這類奇異間斷性問題的核心在于對(duì)間斷處的精確捕捉與計(jì)算.對(duì)于這類問題的求解,偽弧長(zhǎng)方法包括局部偽弧長(zhǎng)方法和全局偽弧長(zhǎng)方法.局部偽弧長(zhǎng)方法為偽弧長(zhǎng)方法的早期形式,其核心在于首先通過間斷解的梯度變換來確定強(qiáng)間斷所處位置,進(jìn)而通過局部網(wǎng)格點(diǎn)重構(gòu)以及數(shù)值修正來達(dá)到強(qiáng)間斷處奇異性消除與降低的目的,如圖1(a)所示.而全局偽弧長(zhǎng)方法(偽弧長(zhǎng)移動(dòng)網(wǎng)格法)則是通過對(duì)整個(gè)計(jì)算區(qū)域內(nèi)的網(wǎng)格點(diǎn)進(jìn)行重構(gòu),使得所有網(wǎng)格點(diǎn)向奇異間斷點(diǎn)處進(jìn)行移動(dòng),從而降低間斷點(diǎn)的影響域,進(jìn)而達(dá)到降低奇異性的目的,如圖1(b)所示.

2.1 局部偽弧長(zhǎng)方法

考慮一維雙曲守恒系統(tǒng)

圖1 雙曲微分方程偽弧長(zhǎng)方法示意圖Fig.1 Schematic diagram for the hyperbolic di ff erential equations pseudo arc-length method

利用雅克比矩陣A(U)=?UF,守恒型方程可以轉(zhuǎn)換到如下形式

引入弧長(zhǎng)參量ξ,其滿足關(guān)系式

其中u是U的分量形式,進(jìn)而由上式可得

進(jìn)而聯(lián)立控制方程(10)與偽弧長(zhǎng)限制方程(13)可以得到

中,采用下式計(jì)算

其中ε是一個(gè)小的正數(shù).進(jìn)而對(duì)式(14)進(jìn)行空間離散,得

在計(jì)算空間中可采用均勻網(wǎng)格,ui=u(ξi,t),Δξ=ξi+1-ξi為計(jì)算空間網(wǎng)格步長(zhǎng),物理空間網(wǎng)格步長(zhǎng)為Δxi=x(ξi,τ+ Δτ)-x(ξi,τ).進(jìn)一步可對(duì)式 (16)利用Runge-Kutta[3536]進(jìn)行時(shí)間離散求解.

為提高上述過程計(jì)算效率,令參數(shù)

其中,Δu=u(xi+1,t)-u(xi,t),Δx為物理空間網(wǎng)格步長(zhǎng),在求解過程中,對(duì)于每個(gè)離散單元,檢驗(yàn)Φ值,當(dāng)Φ≤Φ0(其中Φ0為一個(gè)小的正數(shù)),采用引入弧長(zhǎng)變量ξ的方程,而其他的單元仍采用原有的方程.

2.2 全局偽弧長(zhǎng)方法

由于在高維空間中,間斷點(diǎn)的跟蹤與控制方程在間斷處的變形修正都是十分困難的,因此相比于局部偽弧長(zhǎng)算法,全局偽弧長(zhǎng)方法更適合處理二維、三維問題.在二維空間中,網(wǎng)格變形如圖2所示,在計(jì)算過程中需要保證每個(gè)單元網(wǎng)格面積大于零,避免網(wǎng)格面積為零或?yàn)樨?fù)值的出現(xiàn),相比于二維問題,三維問題要復(fù)雜很多,由于在三維空間中六面體單元畸變以后,其外表面各點(diǎn)會(huì)出現(xiàn)不共面的情況,如圖3所示,這樣會(huì)導(dǎo)致計(jì)算過程中的物理量重構(gòu)不守恒,以及扭曲變形后網(wǎng)格單元的控制方程演化失敗.為解決這一難題,本文采用分塊重構(gòu)與整體計(jì)算相結(jié)合的策略,即在物理量重構(gòu)階段,將六面體單元每個(gè)空間外表面拆分成兩個(gè)面(如圖4過程),則每個(gè)六面體在X,Y,Z每個(gè)方向上被拆分成兩個(gè)三棱柱分別進(jìn)行處理(如圖5所示).在網(wǎng)格移動(dòng)與控制方程的演化計(jì)算階段再回歸整體單元區(qū)域進(jìn)行計(jì)算.分塊重構(gòu)過程可以避免網(wǎng)格的重構(gòu)過程中由于多個(gè)點(diǎn)不在一個(gè)平面導(dǎo)致物理量重構(gòu)不守恒問題.整體移動(dòng)與計(jì)算可以避免多個(gè)塊體分別計(jì)算導(dǎo)致計(jì)算量過大以及多個(gè)塊體的同時(shí)移動(dòng)導(dǎo)致網(wǎng)格錯(cuò)位以及交叉網(wǎng)格的出現(xiàn).

圖2 二維空間網(wǎng)格變形示意圖Fig.2 Two-dimensional spatial grid deformation

圖3 三維空間網(wǎng)格變形示意圖Fig.3 Three-dimensional grid deformation

圖4 三維曲面計(jì)算示意圖Fig.4 Computational diagram of three-dimensional curved surface

對(duì)于三維非定常、無黏、無熱傳導(dǎo)反應(yīng)流體動(dòng)力學(xué)問題,考慮三維空間中一步反應(yīng)流體歐拉方程組

其中,x為三維物理空間向量；三維空間矩陣函數(shù)

物理量 w=(ρ,ρu,ρv,ρw,E,ρR)T；化學(xué)反應(yīng)源項(xiàng)函數(shù)為s(w)=(0,0,0,0,0,ω)T.對(duì)于一步Arrhenius化學(xué)反應(yīng)模型ω

一步化學(xué)反應(yīng)狀態(tài)方程為

引入弧長(zhǎng)參數(shù)ξ=ξ(x,t),其滿足弧長(zhǎng)表達(dá)式

進(jìn)一步,我們可以定義監(jiān)視函數(shù)

通過引入偽時(shí)間來得到多維空間中的網(wǎng)格移動(dòng)速度的偏微分方程

對(duì)上式可以采用Gauss-Seidel循環(huán)迭代求解

其中

對(duì)方程(18)在每個(gè)網(wǎng)格單元上K積分得到并化簡(jiǎn)得

其中

進(jìn)而可對(duì)系統(tǒng)式 (22)和式 (23)采用 Runge-Kutta[3536]耦合求解.耦合求解過程中,要保證網(wǎng)格物理量值插值重構(gòu)的守恒性

其中,wK,K為重構(gòu)前后的物理量值,AK,K為重構(gòu)前后的網(wǎng)格體積.令Dl為式(24)一次循環(huán)求解后,外表面所掃過的體積.所以得到對(duì)于單個(gè)網(wǎng)格的守恒插值策略

其中,Dl為外表面掃過區(qū)域產(chǎn)生的體積改變量,這里以計(jì)算D1為例,D1為圖6所示的三棱柱塊體.D1可以分解為三個(gè)四面體進(jìn)行求解,因此可得

通過四面體體積公式求解出每個(gè)四面體有向體積VOABC，即當(dāng)ABC為逆時(shí)針排列時(shí),其值為正,順時(shí)針排列時(shí),其值為負(fù).

圖6 單個(gè)網(wǎng)格面掃過的體積Fig.6 Swept volume of the typical mesh surface

除此之外,對(duì)于化學(xué)反應(yīng)流問題,在計(jì)算過程中要應(yīng)用保正性策略,保證計(jì)算過程中的壓力、密度、化學(xué)反應(yīng)質(zhì)量分?jǐn)?shù)等為正.筆者在文獻(xiàn)[23]已經(jīng)討論了二維空間中保正性問題,建立了全局偽弧長(zhǎng)保正性算法的體系結(jié)構(gòu).全局偽弧長(zhǎng)算法的保正性理論包括兩個(gè)方面：(1)控制方程離散的保正性;(2)物理量自適應(yīng)重構(gòu)的保正性.其實(shí)現(xiàn)步驟也包括兩個(gè)方面：(1)通過添加時(shí)間步長(zhǎng)和網(wǎng)格移動(dòng)步長(zhǎng)的約束條件實(shí)現(xiàn)網(wǎng)格均值的保正性;(2)通過補(bǔ)充算法實(shí)現(xiàn)網(wǎng)格的每個(gè)點(diǎn)物理量值的保正性.

3 數(shù)值算例

算例1 首先通過具有解析解的一維Sod’s問題來說明全局偽弧長(zhǎng)算法能夠降低激波間斷奇異性.對(duì)于Euler方程(18)的一維初值問題

計(jì)算域取[0,1],最終計(jì)算時(shí)間為T=0.2.計(jì)算過程中取監(jiān)視函數(shù)為

表1 中對(duì)比給出偽弧長(zhǎng)法與固定網(wǎng)格(fi mesh)方法的計(jì)算效果.從表 1中可以看出,利用偽弧長(zhǎng)(pseudo arc-length)法50個(gè)網(wǎng)格點(diǎn)就可以獲得比固定網(wǎng)格方法100個(gè)網(wǎng)格點(diǎn)更好的計(jì)算效果,而且計(jì)算時(shí)間并沒有增加.利用100個(gè)網(wǎng)格點(diǎn)的不同的弧長(zhǎng)參數(shù)的偽弧長(zhǎng)法可以獲得比固定網(wǎng)格方法150,200,300個(gè)更好的計(jì)算效果.特別是L2與L∞誤差范數(shù)顯著降低,說明偽弧長(zhǎng)算法對(duì)于激波間斷奇異性問題有很好的計(jì)算效果.

表1 偽弧長(zhǎng)法與固定網(wǎng)格方法計(jì)算效果對(duì)比Table 1 Computational comparison for pseudo arc-length method and fi ed grid method

算例2 考慮二維板型障礙物爆轟衍射問題.計(jì)算域如圖7所示,紅色為障礙物區(qū)域,其空間位置為[1.6,1.9]×[0,1.0].初始反應(yīng)區(qū)選定為X方向的小于1的區(qū)域.計(jì)算過程中取未反應(yīng)區(qū)R為1,完全反應(yīng)區(qū)R為0.指前因子?K為2566.4,活化能?T為50.一步化學(xué)反應(yīng)狀態(tài)方程為

其中,熱釋放率q為50,反應(yīng)氣體常數(shù)γ為1.2.以Zeldovich Neuman-Doring(ZND)模型[3738]解析解作為反應(yīng)初值.初始計(jì)算域如圖7所示,障礙物為剛性反射邊界,地面為剛性反射邊界,計(jì)算域前后及上方均為無限邊界,初始網(wǎng)格步長(zhǎng)為Δx=Δy=圖8給出了當(dāng)t=0.35時(shí)的模擬結(jié)果.結(jié)果表明,網(wǎng)格隨著爆轟波陣面而進(jìn)行自適應(yīng)變化,實(shí)現(xiàn)了網(wǎng)格對(duì)爆炸沖擊波陣面的捕捉.選定板型障礙物前后12個(gè)位置點(diǎn),根據(jù)位置點(diǎn)與障礙物的位置關(guān)系分成4組，考察其各點(diǎn)處的壓力變化,12點(diǎn)的位置如圖7中所示,其壓力變化如圖9所示.選定一般研究認(rèn)為的壓力安全線3.0為參考?jí)毫€.從圖9中可以看出,在障礙物前方區(qū)域壓力較高,為危險(xiǎn)區(qū)域.在障礙物后方區(qū)域中,位置較高點(diǎn)的爆轟波壓力較大,位置較低的點(diǎn)壓力相對(duì)較小,且越靠近障礙物的后方區(qū)域,其安全性越高.在本文所選取的12個(gè)點(diǎn)中,安全性最高的是e(2.0,0.5)位置的點(diǎn).f(2.0,0.1)位置由于地面反射波的原因,稍有危險(xiǎn).所以在爆炸發(fā)生時(shí),最安全的區(qū)域是障礙物后方的中部位置.

圖7 二維爆轟模擬初態(tài)示意圖Fig.7 Initial state of two-dimensional detonation

圖8 二維板型障礙物繞射問題Fig.8 Two-dimensional di ff raction problem for plate-shape obstacle

圖9 二維板型障礙物繞射各位置點(diǎn)壓力變化Fig.9 Pressure history of typical locations in two-dimensional plate-shape obstacle di ff raction problem

圖9 二維板型障礙物繞射各位置點(diǎn)壓力變化(續(xù))Fig.9 Pressure history of typical locations in two-dimensional plate-shape obstacle di ff raction problem(continued)

算例3 下面分析考察三維空間中的板型障礙物繞射問題.計(jì)算域如圖10(a)所示,紅色為障礙物區(qū)域,其空間位置為[1.6,1.9]×[0,1.0]×[0,1.0].初始反應(yīng)區(qū)選定為X方向的小于1的區(qū)域,其余為未反應(yīng)區(qū).計(jì)算過程中取未反應(yīng)區(qū)R為1,完全反應(yīng)區(qū)R為0.計(jì)算參數(shù)與算例2相同,并且同樣以Zeldovich Neuman-Doring(ZND)模型解析解作為反應(yīng)初值.障礙物為剛性反射邊界,地面為剛性反射邊界,計(jì)算域前后左右及上方均為無限邊界.計(jì)算中采用100萬網(wǎng)格.圖11中給出的是三維空間中壓力與密度云圖,計(jì)算結(jié)果表明障礙物后方存在一個(gè)低密度區(qū).圖12中給出障礙物處切片的網(wǎng)格分布圖,特別是圖12(b)中給出圖12(a)中A部分的網(wǎng)格變形前后的對(duì)比圖,黑色為網(wǎng)格畸變前的,紅色表示網(wǎng)格畸變后的.圖13中給出非障礙物處切片的網(wǎng)格分布圖,可以明顯看出三維空間中網(wǎng)格的移動(dòng)變形.計(jì)算結(jié)果表明,偽弧長(zhǎng)算法能夠?qū)崿F(xiàn)三維空間中網(wǎng)格的自適應(yīng)變化.圖14中,給出圖11中點(diǎn)C處的密度計(jì)算效果對(duì)比,其中參照解是采用2000萬網(wǎng)格的計(jì)算效果,從圖中可以看出,偽弧長(zhǎng)算法的計(jì)算效果與參照解的計(jì)算效果基本一致,且明顯優(yōu)于固定網(wǎng)格算法的100萬網(wǎng)格的計(jì)算效果.

圖10 三維爆轟模擬初態(tài)示意圖Fig.10 Initial state diagram of three-dimensional detonation

圖11 三維板型障礙物繞射問題Fig.11 Three-dimensional di ff raction problem for plate-shape obstacle

圖12 三維板型障礙物繞射問題障礙物處切片F(xiàn)ig.12 Slices at obstacle for three-dimensional plate-shape obstacle di ff raction problem

圖13 三維板型障礙物繞射問題非障礙物處切片F(xiàn)ig.13 Slices at non-obstacle for three-dimensional plate-shape obstacle di ff raction problem

圖14 點(diǎn)C處的密度計(jì)算效果對(duì)比Fig.14 Comparison of densities at point C

進(jìn)而選取障礙物前后X方向的a,b,c,d四個(gè)橫切面,如圖10(a)所示,繼而在每個(gè)橫切面上選取5個(gè)位置點(diǎn),各點(diǎn)的選取如圖10(b)所示.根據(jù)數(shù)值模擬結(jié)果測(cè)定各個(gè)位置點(diǎn)處的壓力變化,如圖15所示.圖15(a)中給出的是障礙物前方a位置處各點(diǎn)的壓力變化曲線.從圖中可以看出,板前區(qū)域的壓力普遍較大,均遠(yuǎn)高于壓力安全線3.0.圖15(b)中給出的是障礙物后方臨近障礙物b位置處的若干點(diǎn)的壓力變化情況.從圖中可以看出位置b4處的壓力是處于安全線以下的.圖15(c)和圖15(d)給出的是障礙物后方c位置和d位置各點(diǎn)的壓力變化曲線.從圖中可以看出,其各點(diǎn)均處于危險(xiǎn)區(qū)域中,但是d位置處的最大壓力峰值相對(duì)滯后.圖16中給出了更為直觀的各個(gè)位置點(diǎn)處的最大壓力分布圖,從圖中可以明顯看出第一組的峰值壓力高于其他各組,最小的峰值壓力在第二組中.

圖15 各位置點(diǎn)處壓力變化Fig.15 Pressure history of typical locations

圖16 各位置點(diǎn)處最大壓力值Fig.16 Maximum pressure at typical locations

4 總結(jié)與討論

本文針對(duì)于計(jì)算動(dòng)力學(xué)問題,分析討論了定常對(duì)流--擴(kuò)散問題的偽弧長(zhǎng)法以及求解非定常對(duì)流問題的偽弧長(zhǎng)法.其中非定常對(duì)流問題的偽弧長(zhǎng)法包括局部偽弧長(zhǎng)法和全局偽弧長(zhǎng)法.針對(duì)三維爆炸沖擊奇異性問題的計(jì)算難點(diǎn),提出了爆炸沖擊問題的三維全局偽弧長(zhǎng)算法,與前人研究的移動(dòng)網(wǎng)格方法不同,由于移動(dòng)網(wǎng)格方法的出發(fā)點(diǎn)是網(wǎng)格的移動(dòng)重構(gòu),缺乏網(wǎng)格的整體性與直觀性,認(rèn)為六面體網(wǎng)格移動(dòng)后,某個(gè)面發(fā)生畸變,變形為非六面體,從而導(dǎo)致計(jì)算失敗.而偽弧長(zhǎng)算法的計(jì)算目標(biāo)不是網(wǎng)格移動(dòng),而是減少奇異性.算例表明,本文提出的六面體網(wǎng)格單元分割與整體結(jié)合的計(jì)算方法成功實(shí)現(xiàn)了三維空間中的爆炸與沖擊問題的數(shù)值模擬.針對(duì)二維,三維爆轟波對(duì)于板型障礙物的爆炸沖擊問題,本文得到了以下結(jié)論：

(1)偽弧長(zhǎng)算法可以有效實(shí)現(xiàn)網(wǎng)格對(duì)于爆轟波陣面的捕捉,從而提高計(jì)算精度,其計(jì)算精度要明顯優(yōu)于固定網(wǎng)格方法.

(2)障礙物前方區(qū)域均為危險(xiǎn)區(qū)域.

(3)障礙物后方臨近障礙物的區(qū)域會(huì)存在某些安全區(qū)域.

(4)在障礙物后方臨近障礙物的區(qū)域中,遠(yuǎn)離障礙物邊緣與地面的地方為最安全區(qū)域.在二維空間中,最安全的區(qū)域是障礙物后方的中部位置.在三維空間中,最安全的區(qū)域是障礙物后方離地面與邊緣最遠(yuǎn)的位置.

1馬天寶,任會(huì)蘭,李健等.爆炸與沖擊問題的大規(guī)模高精度計(jì)算.力學(xué)學(xué)報(bào),2016,48(3)：599-608(Ma Tianbao,Ren Huilan,Li Jian,et al.Large scale high precision computation for explosion and impactproblems.ChineseJournalofTheoreticalandAppliedMechanics,2016,48(3)：599-608(in Chinese))

2楊露斯,黎煉.論計(jì)算機(jī)發(fā)展史及展望.信息與電腦,2010,6(1)：188-188(Yang Lusi,Li Lian.Theory of computer development and prospects.China Computer&Communication,2010,6(1)：188-188(in Chinese))

3 Luo ST,Tran HV,Yu Y.Some inverse problems in periodic homogenization of hamilton-jacobi equations.Archive for Rational Mechanics and Analysis,2016,221(3)：1585-1617

4 Deleuze Y,Chiang CY,Thiriet M,et al. Numerical study of plume patterns in a chemotaxis–di ff usion–convection coupling system.Computers&Fluids,2016,126(1)：58-70

5 Tian DS.Multiple positive periodic solutions for second-order differential equations with a singularity.Acta Applicandae Mathematicae,2016,144(1)：1-10

6 Riks E.An incremental approach to the solution of snapping and buckling problems.International Journal of Solids&Structures,1979,15(7)：529-551

7 Boor CD.On local spline approximation by moments.Journal of Mathematics and Mechanics,1968,17(1)：729-735

8 Crisfiel MA.An arc-length method including line searches and accelerations.International Journal for Numerical Methods in Engineering,1983,19(1)：1268-1289

9 Chan TF.Newton-like pseudo-arclength methods for computing simple turning points.Siam Journal on Scientifi&Statistical Computing,1984,5(1)：135-148

10忻孝康,唐登海.一維定常對(duì)流–擴(kuò)散方程的偽弧長(zhǎng)參數(shù)法.應(yīng)用力學(xué)學(xué)報(bào),1988,5(1)：45-54(Xin Xiaokang,Tang Denghai.An arc-length parameter technique for steady di ff usion-convection equation.Chinese Journal of Applied Mechanics,1988,5(1)：45-54(in Chinese))

11陳國(guó)謙,高智.對(duì)流擴(kuò)散方程的迎風(fēng)變換及相應(yīng)有限差分方法.力學(xué)學(xué)報(bào),1991,23(4)：418-425(Chen Guoqian,Gao Zhi.Transformation of the convective di ff usion equation with corresponding finit di ff erence method.Chinese Journal of Theoretical and Applied Mechanics,1991,23(4)：418-425(in Chinese))

12武際可,許為厚,丁紅麗.求解微分方程初值問題的一種弧長(zhǎng)法.應(yīng)用數(shù)學(xué)和力學(xué),1999,20(8)：875-880(Wu Jike,Xu Weihou,Ding Hongli.Arc-lengthmethodfordi ff erentialequations.AppliedMathematics and Mechanics,1999,20(8)：875-880(in Chinese))

13 Harten A,Hyman JM.Self adjusting grid methods for onedimensional hyperbolic conservation laws.Journal of Computational Physics,1983,50(2)：235-269

14 Ren YH,Russell RD.Moving mesh techniques based upon equidistribution,and their stability.Siam Journal on Scientifi&Statistical Computing,1992,13(6)：1265-1286

15 Huang WZ,Ren YH,Russell RD.Moving mesh partial di ff erential equations(MMPDEs)based upon the equidistribution principle.Siam Journal on Numerical Analysis,1994,31(3)：709-730

16 White AB.On selection of equidistributing meshes for two-point boundary-value problems.Siam Journal on Numerical Analysis,1979,16(3)：472-502

17 Stockie JM,Mackenzie JA,Russell RD.A moving mesh method for one-dimensional hyperbolic conservation laws.Siam Journal on Scientifi Computing,2001,22(5)：1791-1813

18 Tang HZ,Tang T.Adaptive mesh methods for one-and twodimensional hyperbolic conservation laws.Siam Journal on Numerical Analysis,2003,41(2)：487-515

19 Ning JG,Wang X,Ma TB,et al.Numerical simulation of shock wave interaction with a deformable particle based on the pseudo arclength method.Science China Technological Sciences,2015,58(5)：848-857

20王星,馬天寶,寧建國(guó).雙曲偏微分方程的局部偽弧長(zhǎng)方法研究.計(jì)算力學(xué)學(xué)報(bào),2014,31(3)：384-389(Wang Xing,Ma Tianbao,Ning Jianguo.Local pseudo arc-length method for hyperbolic partial di ff erential equation.Chinese Journal of Computational Mechanics,2014,31(3)：384-389(in Chinese))

21 Wang X,Ma TB,Ning JG.A pseudo arc-length method for numerical simulation of shock waves.Chinese Physics Letter,2014,31(3)：030201-030201

22 Wang X,Ma TB,Ren HL,et al.A local pseudo arc-length method for hyperbolic conservation laws.Acta Mechanica Sinica,2014,30(6)：956-965

23 Ning JG,Yuan XP,Ma TB,et al.Positivity-preserving moving mesh scheme for two-step reaction model in two dimensions.Computers&Fluids,2015,123(1)：72-86

24 Yuan XP,Ning JG,Ma TB,et al.Stability of newton tvd runge–kutta scheme for one-dimensional Euler equations with adaptive mesh.Applied Mathematics&Computation,2016,282(1)：1-16

25 Sun LP,Cong YH,Kuang JX.Asymptotic behavior of nonlinear delay di ff erential–algebraic equations and implicit Euler methods.Applied Mathematics&Computation,2014,228(1)：395-403

26 Harwood SM,Kai H,Barton PI.Efficient solution of ordinary differential equations with a parametric lexicographic linear program embedded.Numerische Mathematik,2016,133(4)：623-653

27 Perthame B,Shu CW.On positivity preserving finitvolume schemes for Euler equations.Numerische Mathematik,1996,73(1)：119-130

28 Wang C,Zhang X,Shu CW,et al.Robust high order discontinuous galerkin schemes for two-dimensional gaseous detonations.Journal of Computational Physics,2012,231(2)：653-665

29魏濤,許明田,汪引.求解對(duì)流擴(kuò)散問題的積分方程法.化工學(xué)報(bào),2015,66(10)：3888-3894(Wei Tao,Xu Mingtian,Wang Yin.Integral equation approach to convection-di ff usion problems.Ciesc Journal,2015,66(10)：3888-3894(in Chinese))

30 Zhou XF,Guo J,Li F.Stability,accuracy and numerical di ff usion analysis of nodal expansion method for steady convection di ff usion equation.Nuclear Engineering&Design,2015,295(1)：567-575

31 Sandu A.Rosenbrock methods with an explicit firs stage.International Journal of Computer Mathematics,2015,93(6)：1-18

32馮帆,王自發(fā),唐曉.一個(gè)基于打靶法的大氣污染源反演自適應(yīng)算法.大氣科學(xué),2016,40(4)：719-729(Feng Fan,Wang Zifa,Tang Xiao.Development of an adaptive algorithm based on the shooting method and its application in the problem of estimating air pollutant emissions.Chinese Journal of Atmospheric Sciences,2016,40(4)：719-729(in Chinese))

33 Chen Z,Huang H,Yan J.Third order maximum-principle-satisfying direct discontinuous galerkin methods for time dependent convection di ff usion equations on unstructured triangular meshes.Journal of Computational Physics,2015,308(1)：198-217

34劉朝陽,王振國(guó),孫明波,等.一種求解激波問題的中心差分--WENO混合方法研究.推進(jìn)技術(shù),2016,37(8)：1522-1528(Liu Chaoyang,Wang Zhenguo,Sun Mingbo,et al.Investigation of a hybrid central di ff erence-WENO method for shock wave problems.Journal of Propulsion Technology,2016,37(8)：1522-1528(in Chinese))

35 Shu CW,Osher S.Efficient implementation of essentially nonoscillatory shock-capturing schemes. Journal of Computational Physics,1989,83(1)：32-78

36 Gotlieb S,Shu CW.Total variation diminishing Runge-Kutta schemes.Mathematics of Computation,1996,67(221)：73-85

37 Doring W.On detonation processes in gases.Annals of Physics,1943,43(9)：421-436

38 Fickett W,Wood WW.Flow calculations for pulsating onedimensional detonations.The Physics of Fluids,1966,9(5)：903-916

PSEUDO ARC-LENGTH NUMERICAL ALGORITHM FOR COMPUTATIONAL DYNAMICS1)

Ning Jianguo Yuan Xinpeng Ma Tianbao2)Li Jian
(State Key Laboratory of Explosion Science and Technology，Beijing Institute of Technology，Beijing 100081，China)

Di ff erential equations are often used to describe the dynamic problems.Classical approaches are always hard to solve it in engineering practice due to its characteristics of strong discontinuity,rigidity and shock singularity,among which singularity problem is one of the most difficult and hot issues among scholars.Pseudo arc-length numerical algorithm is proposed for singularity problems in computational dynamics,whose basic idea is to introduce a pseudo arc-length parameter in the original governing equations so that a constraint equation is added.Under the action of a pseudo arc-length parameter,the original discrete elements are distorted to achieve the goal of eliminating or weakening singularity.Firstly,this paper gave an introduction about the pseudo arc-length method for solving the singularity problem in steady di ff usion-convection equations.Then the local pseudo arc-length algorithm is proposed for hyperbolic conservation laws.This algorithm has two steps.The firs step is to determine the location of the strong discontinuity and the second one aims to eliminate or reduce the singularity by reconstructing the local mesh.Secondly,a global pseudoarc-length algorithm is put forward for high dimensional problems,which can reconstruct the mesh in whole area.Since the reconstructed mesh can adaptively catch the singularity points,the singularity is greatly reduced.Thirdly,the threedimensional global pseudo arc-length algorithm and its computational difficulties in numerical algorithm convergence caused by large grid distortion in three-dimensional space are presented.Then the combination of block reconstruction and integral calculation strategy is adopted in the algorithm design process to achieve the pseudo arc-length numerical algorithm in three-dimensional space.Finally,numerical examples were employed to verify the validity of the pseudo arc-length algorithm.

computational dynamics,pseudo arc-length,adaptive,explosion shock wave,three-dimensional

O302,O175

：A

10.6052/0459-1879-16-385

2016–12–19 收稿，2017–03–14 錄用,2017–03–14 網(wǎng)絡(luò)版發(fā)表.

1)國(guó)家自然科學(xué)基金資助項(xiàng)目(11390363,11532012).

2)馬天寶,副教授,主要研究方向：計(jì)算爆炸力學(xué).E-mail：madabal@bit.edu.cn

寧建國(guó),原新鵬,馬天寶,李健.計(jì)算動(dòng)力學(xué)中的偽弧長(zhǎng)方法研究.力學(xué)學(xué)報(bào),2017,49(3)：703-715

Ning Jianguo,Yuan Xinpeng,Ma Tianbao,Li Jian.Pseudo arc-length numerical algorithm for computational dynamics.Chinese Journal of Theoretical and Applied Mechanics,2017,49(3)：703-715