分片光滑物態(tài)方程的能量方程非線性迭代解法

杭旭登，李雙貴，楊 容，袁光偉，2，?

（1.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京市海淀區(qū)豐豪東路2號(hào)，100094；2.計(jì)算物理實(shí)驗(yàn)室，北京市海淀區(qū)花園路6號(hào)，100088）

文章編號(hào)：1001?246X（2015）05?0505?09

分片光滑物態(tài)方程的能量方程非線性迭代解法

杭旭登1，李雙貴1，楊 容1，袁光偉1，2，?

（1.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京市海淀區(qū)豐豪東路2號(hào)，100094；2.計(jì)算物理實(shí)驗(yàn)室，北京市海淀區(qū)花園路6號(hào)，100088）

實(shí)際應(yīng)用中的物態(tài)方程由分片光滑曲面拼接而成，拼接處存在間斷.隱式求解相應(yīng)的能量方程時(shí)，經(jīng)常出現(xiàn)迭代收斂慢的情況和非物理解.本文通過(guò)構(gòu)造對(duì)應(yīng)的新的非線性問(wèn)題，提出一種非線性迭代算法.該算法適用于求解有間斷的分片光滑物態(tài)方程的非線性能量方程，其中引入一個(gè)度量能量變化的參數(shù)用于自動(dòng)判斷跳段是否發(fā)生，在求解時(shí)無(wú)需事先知道物態(tài)方程間斷的位置，且能精確計(jì)算物態(tài)方程間斷帶來(lái)的能量盈虧，用于評(píng)估物態(tài)方程間斷對(duì)能量的影響.典型算例驗(yàn)證了新算法具有穩(wěn)定的收斂性，并給出符合物理規(guī)律的解.

間斷狀態(tài)方程；能量方程；非線性迭代算法

0 引言

在慣性約束聚變（ICF）等實(shí)際問(wèn)題中，系統(tǒng)的物質(zhì)狀態(tài)（密度、溫度等）在很小的時(shí)空區(qū)域中發(fā)生很大的變化.物態(tài)方程e＝e（ρ，T）描述物質(zhì)內(nèi)能e與密度ρ和溫度T等熱力學(xué)量之間的依賴關(guān)系.當(dāng)物態(tài)經(jīng)歷固體、液體、氣體和等離子體等不同狀態(tài)時(shí)，物態(tài)方程由不同的函數(shù)形式表示；另外，在不同的溫度段，由參數(shù)庫(kù)擬合的物態(tài)方程也是采用不同的函數(shù)表示［1－2］.在不同狀態(tài)的連接處，這些函數(shù)表達(dá)式給出的曲面（或曲線）無(wú)法光滑拼接，導(dǎo)致能量作為溫度的函數(shù)是間斷的，此現(xiàn)象稱為物態(tài)方程的跳段.根據(jù)跳段處兩邊能量值的大小關(guān)系，可以分為向上跳段和向下跳段兩種情況，典型的向下跳段的情況如圖1所示，其中e＝e（ρ，T）在T＝T0處出現(xiàn)跳段.

圖1 向下跳段示意圖Fig.1 A“jump down”discontinuity in EOS

向下跳段是一種非物理的現(xiàn)象.但在高能量密度情形，實(shí)際采用的物態(tài)方程中能量函數(shù)在跳段處的間斷量遠(yuǎn)小于該處能量函數(shù)的值，物理上認(rèn)為該跳段對(duì)整體物理過(guò)程的影響較小.因此，帶跳段的物態(tài)方程仍廣泛應(yīng)用于實(shí)際問(wèn)題的數(shù)值模擬中.

在ICF研究中，求解非平衡非線性能量方程

常用的非線性迭代方法是Newton?Picard迭代算法［6－7，11］.Newton?Picard迭代算法對(duì)非平衡能量方程中的耦合項(xiàng)采用Newton展開，而對(duì)于擴(kuò)散項(xiàng)則采用Picard迭代.Knoll等［8，10］提出了無(wú) Jacobi的Newton?Krylov（JFNK）算法，它對(duì)于整個(gè)方程組采用Newton迭代，其中對(duì)Jacobi矩陣乘以向量，采用非線性函數(shù)的有限差分來(lái)逼近.該算法對(duì)于光滑問(wèn)題收斂速度較快，但由于偏微分方程離散格式導(dǎo)致的非線性函數(shù)的計(jì)算非常費(fèi)時(shí)，并且還受限于Newton算法的局部收斂性，在實(shí)際應(yīng)用中其計(jì)算效率并不高，穩(wěn)定性也較差.袁光偉等［9］提出了Picard?Newton迭代算法，主要思想是對(duì)于擴(kuò)散項(xiàng)在連續(xù)形式下先線性化，再對(duì)線性化的方程中不同性質(zhì)的項(xiàng)進(jìn)行有針對(duì)性的離散.由于選取了合適的離散格式，該算法不僅收斂快、效率較高，而且提高了求解的穩(wěn)定性.在這些研究中，對(duì)于擴(kuò)散項(xiàng)和耦合項(xiàng)的研究較多，而對(duì)于能量函數(shù)的非線性迭代往往進(jìn)行簡(jiǎn)化而忽略了，對(duì)于間斷的能量函數(shù)的研究更少.

物態(tài)方程中的能量函數(shù)（電子能量、離子能量）的間斷對(duì)能量方程的迭代求解帶來(lái)了很大的困難.在迭代過(guò)程中，當(dāng)溫度越過(guò)兩段物態(tài)方程的連接處（即跳段處）時(shí)，由于能量和比熱等參數(shù)在跳段處間斷，極易導(dǎo)致迭代不收斂，這在實(shí)際問(wèn)題的求解中非常普遍.常用的處理方法是將每一段物態(tài)方程在跳段處進(jìn)行外插，在迭代結(jié)束之后再進(jìn)入下一段物態(tài)方程曲線.這種做法容易導(dǎo)致迭代溫度在跳段處在不同時(shí)間步之間來(lái)回變化，導(dǎo)致不同時(shí)間步之間的解不斷振蕩，且不能精確計(jì)算在跳段處帶來(lái)的能量盈虧，難以考察跳段對(duì)整體計(jì)算結(jié)果的影響，不利于評(píng)估整個(gè)模擬過(guò)程的置信度.

為了解決物態(tài)方程跳段引起的數(shù)值計(jì)算收斂慢甚至不收斂的問(wèn)題，并在計(jì)算中對(duì)于跳段帶來(lái)的能量盈虧進(jìn)行精確的統(tǒng)計(jì)，本文提出一種適用于計(jì)算跳段過(guò)程的非線性迭代算法.對(duì)不出現(xiàn)跳段的光滑問(wèn)題，該算法退化到經(jīng)典的Newton迭代；對(duì)跳段問(wèn)題，與傳統(tǒng)算法相比，該算法迭代收斂更快，并能精確地統(tǒng)計(jì)跳段過(guò)程中的能量盈虧.

1 能量跳段問(wèn)題的傳統(tǒng)迭代算法

在輻射流體力學(xué)計(jì)算中，內(nèi)能形式的能量方程為其中內(nèi)能e＝e（ρ，T）和比熱cV＝（?e（ρ，T）/?T）ρ都是溫度和密度的非線性函數(shù)，p是物質(zhì)壓強(qiáng)，u是物質(zhì)速度，κ是熱傳導(dǎo)系數(shù)，Q是源項(xiàng).在求解輻射流體力學(xué)問(wèn)題時(shí)，通常先計(jì)算流體部分，確定速度和密度，在能量方程計(jì)算時(shí)密度保持不變.因此，只需考慮由于溫度變化引起的物態(tài)方程跳段問(wèn)題，以及該問(wèn)題對(duì)能量方程數(shù)值求解的影響.

方程（3）含比熱系數(shù)，不直接含能量函數(shù).若n記為時(shí)間層，則常用的時(shí)間離散為

由此，基于能量方程（3）相應(yīng)的Picard非線性迭代為

其中s表示非線性迭代次數(shù)．注意到，物理上cV（T）＞0，該方法可保證當(dāng)系統(tǒng)能量加源為正時(shí)溫度會(huì)上升，這符合物理規(guī)律的要求．不足之處是需要計(jì)算T?p/?T，并且易出現(xiàn)迭代不收斂的現(xiàn)象：考慮圖1的情況，在跳段處左邊的比熱小，右邊的比熱大，則在適當(dāng)?shù)膲毫ψ鞴图釉辞樾?，溫度T從小于T0（此時(shí)采用跳段左邊的比熱）剛好上升越過(guò)跳段到右邊（T＞T0），然后，下一步迭代時(shí)由于采用右邊的比熱，使得溫度返回到跳段前（T＜T0），即出現(xiàn)溫度振蕩，從而導(dǎo)致迭代難以收斂甚至不收斂.

為解決上述迭代不收斂的問(wèn)題，通常的處理方法是將跳段一側(cè)的物態(tài)方程向另一側(cè)光滑延拓（保持一階導(dǎo)數(shù)在跳段處連續(xù)），如圖2所示.當(dāng)初始溫度在跳段左邊、下一迭代步溫度進(jìn)入跳段右邊時(shí)，采用跳段左邊的比熱進(jìn)行能量函數(shù)的延拓（如圖2中向右上虛箭頭所示）；而當(dāng)初始溫度在右邊、下一迭代步溫度進(jìn)入跳段左邊時(shí)，采用跳段右邊的比熱進(jìn)行計(jì)算（如圖2中向左下的虛箭頭所示）.

為討論簡(jiǎn)單起見，忽略流體運(yùn)動(dòng).此時(shí)，對(duì)于網(wǎng)格j，考慮由于求解（4）導(dǎo)致的j網(wǎng)格上的能量變化

圖2 傳統(tǒng)的能量跳段處理方式Fig.2 Traditional approach for discontinuity in EOS

其中Vj是網(wǎng)格j的體積.而通過(guò)網(wǎng)格的溫度變化計(jì)算出的能量變化為

這兩種計(jì)算方式得到的能量相差

其中Dv，j（）是能量函數(shù)關(guān)于溫度的二階導(dǎo)數(shù)，是介于，之間的某個(gè)溫度值.由此可知，如果一個(gè)網(wǎng)格上的溫度發(fā)生變化，采用上面的方法求解導(dǎo)致的能量的守恒誤差與溫度差的平方和能量關(guān)于溫度的二階導(dǎo)數(shù)的乘積成正比，從而無(wú)法保證系統(tǒng)的能量嚴(yán)格守恒.為了得到所需的守恒性，需要縮短時(shí)間步長(zhǎng)使得溫度變化非常小.

這一傳統(tǒng)的處理方法使得能量函數(shù)在一個(gè)時(shí)間步內(nèi)連續(xù)可微，迭代過(guò)程中不發(fā)生來(lái)回振蕩，使得非線性迭代收斂.然而，該方法設(shè)計(jì)中只在跳段處保持比熱連續(xù)，延拓后相應(yīng)的比熱與原來(lái)的比熱沒(méi)有直接關(guān)系，這導(dǎo)致如下兩個(gè)問(wèn)題：①如上所述，算法不能保證能量嚴(yán)格守恒；②如果除了能量不連續(xù)外，跳段處左右兩邊的比熱也不連續(xù)，則溫度下降經(jīng)過(guò)跳段與溫度上升經(jīng)過(guò)跳段處導(dǎo)致的能量盈虧的大小不相等（見圖2）.它們都可以與能量跳段值不相等，此時(shí)，數(shù)值解的能量守恒誤差難以精確統(tǒng)計(jì)，不利于準(zhǔn)確評(píng)估跳段現(xiàn)象對(duì)整體能量守恒的影響.

針對(duì)物態(tài)方程的跳段，如果采用將間斷函數(shù)進(jìn)行人為的磨光連接處理的方法，也將導(dǎo)致如下困難：一方面，人為的磨光函數(shù)不能保證物態(tài)方程物理上的自洽；另一方面，也不能直接解決迭代求解的問(wèn)題.對(duì)向上跳段情形，會(huì)導(dǎo)致局部的比熱太大，在迭代過(guò)程中溫度變化即使很小，也會(huì)使得能量增加很多；對(duì)于向下跳段，會(huì)導(dǎo)致出現(xiàn)負(fù)的比熱，產(chǎn)生非物理的現(xiàn)象：溫度降低能量反而增加了.

總之，上述兩種處理方法都不能很好地描述物理問(wèn)題，也不能解決迭代收斂的問(wèn)題.本文針對(duì)上述物態(tài)方程的跳段引起的迭代求解困難的問(wèn)題，引入輔助能量函數(shù)，將問(wèn)題轉(zhuǎn)化為一個(gè)新的非線性問(wèn)題，形成相應(yīng)的迭代算法.該算法可以結(jié)合流體計(jì)算進(jìn)行，為簡(jiǎn)單起見以下部分忽略流體運(yùn)動(dòng)對(duì)于能量計(jì)算的影響.

2 適應(yīng)跳段的迭代算法設(shè)計(jì)

其中m表示質(zhì)量．對(duì)e（T）采用Newton線性化（上式右端項(xiàng)仍采用簡(jiǎn)單的Picard線性化），得到

其中s表示非線性迭代次數(shù).在不會(huì)引起混淆的情形，下面對(duì)非線性迭代有時(shí)將忽略時(shí)間上標(biāo).

在流體采用Lagrange描述時(shí)mn＋1＝mn，非線性迭代收斂時(shí)此方程保持能量守恒.在出現(xiàn)能量和比熱跳段時(shí)，與非線性迭代（3）類似，傳統(tǒng)的Newton迭代法（10）依然迭代不收斂，或者在跳段局部找不到符合物理的解.如圖3所示，在求解方程（10）時(shí)，設(shè)前一迭代步溫度T（s）在T0左側(cè)附近，如果方程（10）右端為正使得溫度越過(guò)跳段T0，那么由能量的守恒性，圖3中大于跳段左側(cè)的能量只能在遠(yuǎn)離跳段T0的右側(cè)（T（s＋1）?T0）才能找到，該解遠(yuǎn)大于T（s），它是非物理的解.

下面我們構(gòu)造新的算法，給出迭代收斂且符合物理的解.為此，由式（9）構(gòu)造如下的非線性問(wèn)題

對(duì)應(yīng)的線性化方程（右端項(xiàng)仍簡(jiǎn)單線性化）為

其中f＝f（T）＝f＋（T）或f＝f－（T），它們分別定義如下

構(gòu)造f（x）的基本思想是：將能量函數(shù)在跳段兩側(cè)上下平移使得兩段曲線連續(xù)拼接（f＋和f－如圖4所示），該函數(shù)用于非線性迭代求解過(guò)程中，稱其為迭代能量函數(shù).

圖3 跳段導(dǎo)致非物理解Fig.3 Unphysical solution due to discontinuity in EOS

圖4 新的迭代能量函數(shù)fFig.4 New iterative function of energy

在線性化迭代方程（12）中，當(dāng)Tn＋1（s）≠T0，（s≥0）時(shí)，cV（Tn＋1（s）） ＝?e（Tn＋1（s））/?T；如果迭代值Tn＋1（s）＝T0（s＞1），則cV（Tn＋1（s））取為與Tn＋1（s－1）同一側(cè)能量函數(shù)的單側(cè)導(dǎo)數(shù).

由式（12）解出Tn＋1（s＋1）后，判斷是否發(fā)生跳段．如未發(fā)生跳段，取f（Tn＋1（s＋1））＝e（Tn＋1（s＋1））；否則取f（Tn＋1（s＋1））＝f（Tn＋1（s））＋cV（Tn＋1（s））（Tn＋1（s＋1）－Tn＋1（s））．

這里算法的關(guān)鍵是判斷溫度是否經(jīng)過(guò)跳段.本文設(shè)計(jì)如下的跳段判斷：當(dāng)小于0，或者大于某個(gè)大于1的常數(shù)（通常取為5），則判斷為經(jīng)過(guò)跳段.這兩種情況分別對(duì)應(yīng)于向下跳段和向上跳段，如圖5所示.

圖5 向下跳段和向上跳段的判斷Fig.5 Judgement of jump up and jump down

下面給出能量方程新的非線性迭代算法（僅限于與能量函數(shù)相關(guān)的部分）.

適用于狀態(tài)方程跳段問(wèn)題的非線性迭代算法

1）取Tn＋1（0）＝Tn和f（Tn＋1（0）） ＝e（Tn＋1（0）），確定采用f＋還是f－，計(jì)算cV（Tn＋1（0））；

2）非線性迭代開始，s＝0，1，…；

3）計(jì)算cV（Tn＋1（s）），求解方程（12），解出Tn＋1（s＋1）；

4）計(jì)算e（Tn＋1（s＋1））和跳段條件（14），自動(dòng)判斷是否越過(guò)跳段：

① 未經(jīng)跳段時(shí)，令f（Tn＋1（s＋1）） ＝e（Tn＋1（s＋1））；

② 如經(jīng)過(guò)跳段，令f（Tn＋1（s＋1）） ＝ f（Tn＋1（s）） ＋cV（Tn＋1（s））（Tn＋1（s＋1）－Tn＋1（s））；

5）判斷迭代是否收斂，如未收斂，轉(zhuǎn)3）；

6）如迭代收斂，計(jì)算能量盈虧為Δen＋1（s＋1）＝e（Tn＋1（s＋1））－f（Tn＋1（s＋1））．

由本算法可知其具有如下性質(zhì)：加源為正則溫度升高，加源為負(fù)則溫度降低，迭代能量函數(shù)連續(xù)，迭代解不發(fā)生大的跳躍，符合物理規(guī)律.對(duì)于迭代能量函數(shù)，本算法只涉及單個(gè)網(wǎng)格本身的量，在每個(gè)時(shí)間層不同網(wǎng)格采用的迭代能量函數(shù)相互獨(dú)立.在任一迭代步都可以計(jì)算兩個(gè)能量函數(shù)e（T）和f（T）的差

如果求解越過(guò)跳段，則此值剛好是跳段處的能量盈虧±δ；如果求解越過(guò)跳段后又折回，則該值為0；由此可見，該算法可以精確統(tǒng)計(jì)跳段的能量，適用于分片光滑的物態(tài)方程跳段問(wèn)題，也適用于光滑問(wèn)題，并且跳段判斷和迭代結(jié)合在一起，求解時(shí)無(wú)需事先知道物態(tài)方程跳段的具體位置，能自動(dòng)判斷；另外，該算法能直接推廣到多個(gè)跳段的情形.

迭代能量函數(shù)f（T）是分片光滑的連續(xù)函數(shù)，相比原間斷函數(shù)，迭代求解顯然更穩(wěn)定.對(duì)于光滑問(wèn)題，在適當(dāng)條件下Newton迭代是局部二階收斂的，在某些特殊情形，Newton迭代全局收斂，參見文獻(xiàn)［4－5］.在文［4］中，當(dāng)非線性函數(shù)的一階導(dǎo)數(shù)為L(zhǎng)ipschitz連續(xù)時(shí)給出了收斂性結(jié)果.

本文中迭代能量函數(shù)f形式上將兩個(gè)光滑函數(shù)連在一起，如果兩段函數(shù)分別滿足上述收斂性定理的條件，則迭代過(guò)程也收斂.圖6給出了三個(gè)例子，說(shuō)明對(duì)于迭代能量函數(shù)f迭代是收斂的.

3 數(shù)值算例

算例1 考慮一維熱傳導(dǎo)方程

圖6 新算法的迭代收斂軌跡Fig.6 Convergence trajectory of the new algorithem

其中E（T）＝Ee（T）＋Ei（T）＋aT4，Ee，Ei分別為電子內(nèi)能和離子內(nèi)能，離子能量Ei（T）＝0.01T，電子能量函數(shù)如圖7所示，其中cVe＝0.09，在T0＝0.6處有一個(gè)向上的跳段，間斷大小為0.5．a(chǎn)T4是輻射能量，a＝5.5656×10－3．計(jì)算區(qū)域?yàn)椋?，50］，采用100個(gè)網(wǎng)格均勻剖分．初值條件為T｜t＝0＝3×10－4；區(qū)域右邊界為絕熱條件Tx＝0，左邊界為凈流－κTx＝acTb4/4，其中c＝2.997 924 58×104為光速，Tb≡1，該能量源使得系統(tǒng)的溫度上升并經(jīng)過(guò)跳段．熱傳導(dǎo)系數(shù)為κ＝0.011＋400acT3/3．

分別采用原方法（10）和新方法（12）進(jìn)行計(jì)算，線性方程采用追趕法進(jìn)行求解，非線性迭代的收斂準(zhǔn)則是溫度的相對(duì)變差滿足

圖7 算例1電子能量的跳段Fig.7 Discontinuity of electron energy in Example 1

整個(gè)過(guò)程計(jì)算到時(shí)間t＝1，時(shí)間步長(zhǎng)為Δt＝10－4．

采用原方法（10）的計(jì)算結(jié)果如圖8所示.

圖8 原算法的非線性迭代次數(shù)（a）、能量守恒與能量漏失（b）、溫度分布（c）Fig.8 Nonlinear iterations，conservation error，energy leakage and temperature profiles by the old algorithm

從圖8（a）可以看出：迭代次數(shù)在時(shí)間0.3以后有一個(gè)較大的增長(zhǎng)；此時(shí)溫度T超過(guò)0.6而發(fā)生跳段，非線性迭代不易收斂．迭代收斂的解已不符合能量非負(fù)的物理要求．圖8（c）給出了時(shí)間t＝0.01，i×0.05，（i ＝1，…7，8，10，…，20）共14個(gè)時(shí)刻的溫度分布，從i＝7開始，數(shù)值解都是負(fù)值.定義能量漏失為

圖8（b）給出能量守恒和能量漏失的統(tǒng)計(jì)結(jié)果，表明整個(gè)系統(tǒng)是能量守恒的，沒(méi)有能量漏失.當(dāng)溫度經(jīng)過(guò)跳段點(diǎn)后，由于系統(tǒng)的總能量突然升高，為了保持能量守恒，由方程（10）可知，系統(tǒng)的溫度只能降低.注意到在總內(nèi)能中輻射能量aT4恒為正，因此，只有將電子和離子溫度降低到負(fù)值，才能達(dá)到系統(tǒng)能量平衡.本算例表明，采用原來(lái)的方法進(jìn)行迭代求解，一方面難以收斂，另一方面即使收斂，也是收斂到非物理的解.

采用本文提出的新算法的計(jì)算結(jié)果如圖9所示.

圖9 新算法的非線性迭代次數(shù)（a）、能量守恒與能量漏失（b）、溫度分布（c）Fig.9 Nonlinear iterations，conservation error，energy leakage and temperature profiles by the new algorithm

圖9（a）表明，新算法的非線性迭代收斂較快，整個(gè)過(guò)程計(jì)算穩(wěn)定.圖9（b）表明能量的盈虧剛好是100×0.5×0.5＝25.在扣除了這些能量漏失后得到的能量守恒誤差非常小，保證了計(jì)算精度，驗(yàn)證了能量漏失統(tǒng)計(jì)的精確性.圖9（c）給出了溫度隨時(shí)間變化的分布，表明新算法避免了原算法存在的能量出負(fù)現(xiàn)象，溫度逐步升高的軌跡平滑，符合能量增加溫度升高的物理規(guī)律.

算例2 方程及初邊界條件與算例1基本相同，不同之處：計(jì)算區(qū)域?yàn)椋?，0.1］，采用100個(gè)網(wǎng)格等分.離子能量為Ei（T）＝0.1T，電子能量的比熱cVe＝0.9，采用跳段的物態(tài)方程，在T0＝0.6處向下跳段－0.5，如圖10所示.

采用原方法（10）的計(jì)算結(jié)果如圖11所示.

由于計(jì)算區(qū)域較小，此問(wèn)題的部分網(wǎng)格的溫度很快通過(guò)跳段點(diǎn)，直到時(shí)間0.4以后，所有的網(wǎng)格溫度都通過(guò)了跳段點(diǎn).在時(shí)間0.4以前迭代次數(shù)較多，主要原因是：此時(shí)段部分網(wǎng)格的溫度T超過(guò)0.6而發(fā)生跳段，非線性迭代不易收斂.實(shí)際上，在溫度經(jīng)過(guò)0.6的區(qū)域，可以看到由于跳段導(dǎo)致溫度出現(xiàn)不光滑的分布（見圖11（c）的底部）.圖11（b）表明整個(gè)系統(tǒng)的能量守恒且沒(méi)有漏失.

圖10 算例2的電子能量跳段Fig.10 Discontinuity of electron energy in Example 2

采用本文提出的新算法的計(jì)算結(jié)果如圖12所示.

圖11 原算法的非線性迭代次數(shù)（a）、能量守恒與能量漏失（b）、溫度分布（c）Fig.11 Nonlinear iterations，conservation error，energy leakage and temperature profiles by the old algorithm

圖12（a）表明，新算法的非線性迭代收斂較快，整個(gè)過(guò)程計(jì)算穩(wěn)定.圖12（b）表明能量漏失剛好是0.1×（－0.5）＝－0.05.在扣除了這些能量漏失后得到的能量守恒誤差非常小，確實(shí)保證了計(jì)算精度，也驗(yàn)證了整個(gè)能量漏失統(tǒng)計(jì)的精確性.圖12（c）給出了溫度隨時(shí)間變化的分布，溫度逐步升高的軌跡平滑，不存在原算法中出現(xiàn)的溫度分布抖動(dòng)的現(xiàn)象，這符合物理規(guī)律.

圖12 新算法的非線性迭代次數(shù)（a）、能量守恒與能量漏失（b）、溫度分布（c）Fig.12 Nonlinear iterations，conservation error，energy leakage and temperature profiles by the new algorithm

4 小結(jié)

一些應(yīng)用問(wèn)題的數(shù)值模擬中經(jīng)常出現(xiàn)物態(tài)方程跳段問(wèn)題，以往的算法導(dǎo)致迭代收斂慢或不收斂.針對(duì)該問(wèn)題，本文簡(jiǎn)要介紹和分析了已有的處理方法，提出了一個(gè)新非線性問(wèn)題和相應(yīng)的非線性迭代算法.該算法設(shè)計(jì)適用于求解分片光滑的物態(tài)方程的問(wèn)題，因此有廣泛的適應(yīng)性.本文給出了數(shù)值算例，分析和比較了新算法與原算法的計(jì)算結(jié)果.計(jì)算表明，對(duì)于物態(tài)方程向上跳段和向下跳段，新算法都比原算法高效，并得到合理的解.新算法可以應(yīng)用于更復(fù)雜的情況，例如慣性約束聚變數(shù)值模擬中，物態(tài)方程由實(shí)際氣體物態(tài)方程向理想氣體物態(tài)方程過(guò)渡的跳段問(wèn)題.

致謝：感謝審稿人提出的修改建議.

［1］ 湯文輝，張若棋.物態(tài)方程理論及計(jì)算概論［M］.第二版.北京：高等教育出版社，2008.

［2］ 徐錫申，等.實(shí)用物態(tài)方程理論導(dǎo)引［M］.北京：科學(xué)出版社，1986.

［3］ 常鐵強(qiáng)，張鈞，等.激光等離子體相互作用與激光聚變［M］.長(zhǎng)沙：湖南科學(xué)技術(shù)出版社，1991.

［4］ Dennis JE，Schnabel R B.Numericalmethods for unconstrained optimization and nonlinear equations［M］.Prentice?Hall，Englewood Cliffs，NJ，1983.

［5］ Deuflhard P.Newton methods for nonlinear problems，affine invariance and adaptive algorithms［M］.Springer，2004.

［6］ RiderW J，Knoll D A.Amultigrid Newton?Krylovmethod formultimaterial equilibrium radiation diffusion［J］.JComp Phys，1999，152（1）：164－191.

［7］ Knoll D A，Chacon L.On balanced approximations for time integration of multiple time scale systems［J］.JComp Phys，2003，185（2）：583.

［8］ Mousseau V A，Knoll D A.Physics?based preconditioning and the Newton?Krylov method for non?equilibrium radiation diffusion［J］.JComp Phys，2000，160（2）：743－765.

［9］ Yuan GW，Hang X D.Accelerating nonlinear iteration for nonlinear parabolic equations［J］.JCM 2006，24（3）：412－424. ［10］ Knoll D A，Keyes D E.Jacobian?free Newton?Krylov methods：A survey of approaches and applications［J］.JComp Phys，2004，193（1）：357－397.

［11］ Brown P N，Shumaker D E.Fully implicit solution of large?scale non?equilibrium radiation diffusion with high order time integration［J］.JComp Phys，2005，204（2）：760－783.

A Nonlinear Iterative M ethod for Energy Equations w ith Piecew ise Smooth EOS

HANG Xudeng1，LIShuanggui1，YANG Rong1，YUAN Guangwei1，2
（1.Institute of Applied Physics and Computational Physics，No.2 Fenghao Donglu，Haidian，Beijing 100094，China；2.Laboratory ofComputational Physics，No.6 Huayuan Road，Haidian，Beijing 100088，China）

In practical applications，equation of states（EOS）consists of several piecewise smooth surfaces，which leads to discontinuity at interface.As a traditional nonlinear iterative algorithm is applied to an energy equation with discontinuous EOS，itmay lead to slow convergence and unphysical solutions.To overcome the difficulties，a nonlinear problem is designed，and a nonlinear iterative algorithm is proposed to solve the problem.The algorithm is fit for energy equations with discontinuous EOS of piecewise smooth functions.A parameter of energy change is defined in the algorithm so that it is unnecessary to know discontinuity position in advance.The algorithm calculates precisely net gain or leakage of energy，which can be used to assess influence of discontinuity in EOS.Typical numerical experiments verify that the algorithm converges stably，and gives physical solutions.

discontinuous EOS；energy equation；nonlinear iterative algorithm

O242

A

2014－08－04；

2014－11－29

國(guó)家自然科學(xué)基金（11171036，11001023）和中物院科學(xué)技術(shù)發(fā)展基金（2014A0202009）資助項(xiàng)目

杭旭登（1976－），男，江蘇無(wú)錫，博士，研究員，從事輻射輸運(yùn)數(shù)值模擬研究，E?mail：hang_xudeng＠iapcm.ac.cn?通訊作者：袁光偉，E?mail：yuan_guangwei＠iapcm.ac.cn

Received date：2014－08－04；Revised date： 2014－11－29