分片光滑物態方程的能量方程非線性迭代解法

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

（1.北京應用物理與計算數學研究所，北京市海淀區豐豪東路2號，100094；2.計算物理實驗室，北京市海淀區花園路6號，100088）

文章編號：1001?246X（2015）05?0505?09

分片光滑物態方程的能量方程非線性迭代解法

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

（1.北京應用物理與計算數學研究所，北京市海淀區豐豪東路2號，100094；2.計算物理實驗室，北京市海淀區花園路6號，100088）

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

間斷狀態方程；能量方程；非線性迭代算法

0 引言

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

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

向下跳段是一種非物理的現象.但在高能量密度情形，實際采用的物態方程中能量函數在跳段處的間斷量遠小于該處能量函數的值，物理上認為該跳段對整體物理過程的影響較小.因此，帶跳段的物態方程仍廣泛應用于實際問題的數值模擬中.

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

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

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

為了解決物態方程跳段引起的數值計算收斂慢甚至不收斂的問題，并在計算中對于跳段帶來的能量盈虧進行精確的統計，本文提出一種適用于計算跳段過程的非線性迭代算法.對不出現跳段的光滑問題，該算法退化到經典的Newton迭代；對跳段問題，與傳統算法相比，該算法迭代收斂更快，并能精確地統計跳段過程中的能量盈虧.

1 能量跳段問題的傳統迭代算法

在輻射流體力學計算中，內能形式的能量方程為其中內能e＝e（ρ，T）和比熱cV＝（?e（ρ，T）/?T）ρ都是溫度和密度的非線性函數，p是物質壓強，u是物質速度，κ是熱傳導系數，Q是源項.在求解輻射流體力學問題時，通常先計算流體部分，確定速度和密度，在能量方程計算時密度保持不變.因此，只需考慮由于溫度變化引起的物態方程跳段問題，以及該問題對能量方程數值求解的影響.

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

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

其中s表示非線性迭代次數．注意到，物理上cV（T）＞0，該方法可保證當系統能量加源為正時溫度會上升，這符合物理規律的要求．不足之處是需要計算T?p/?T，并且易出現迭代不收斂的現象：考慮圖1的情況，在跳段處左邊的比熱小，右邊的比熱大，則在適當的壓力作功和加源情形，溫度T從小于T0（此時采用跳段左邊的比熱）剛好上升越過跳段到右邊（T＞T0），然后，下一步迭代時由于采用右邊的比熱，使得溫度返回到跳段前（T＜T0），即出現溫度振蕩，從而導致迭代難以收斂甚至不收斂.

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

為討論簡單起見，忽略流體運動.此時，對于網格j，考慮由于求解（4）導致的j網格上的能量變化

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

其中Vj是網格j的體積.而通過網格的溫度變化計算出的能量變化為

這兩種計算方式得到的能量相差

其中Dv，j（）是能量函數關于溫度的二階導數，是介于，之間的某個溫度值.由此可知，如果一個網格上的溫度發生變化，采用上面的方法求解導致的能量的守恒誤差與溫度差的平方和能量關于溫度的二階導數的乘積成正比，從而無法保證系統的能量嚴格守恒.為了得到所需的守恒性，需要縮短時間步長使得溫度變化非常小.

這一傳統的處理方法使得能量函數在一個時間步內連續可微，迭代過程中不發生來回振蕩，使得非線性迭代收斂.然而，該方法設計中只在跳段處保持比熱連續，延拓后相應的比熱與原來的比熱沒有直接關系，這導致如下兩個問題：①如上所述，算法不能保證能量嚴格守恒；②如果除了能量不連續外，跳段處左右兩邊的比熱也不連續，則溫度下降經過跳段與溫度上升經過跳段處導致的能量盈虧的大小不相等（見圖2）.它們都可以與能量跳段值不相等，此時，數值解的能量守恒誤差難以精確統計，不利于準確評估跳段現象對整體能量守恒的影響.

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

總之，上述兩種處理方法都不能很好地描述物理問題，也不能解決迭代收斂的問題.本文針對上述物態方程的跳段引起的迭代求解困難的問題，引入輔助能量函數，將問題轉化為一個新的非線性問題，形成相應的迭代算法.該算法可以結合流體計算進行，為簡單起見以下部分忽略流體運動對于能量計算的影響.

2 適應跳段的迭代算法設計

其中m表示質量．對e（T）采用Newton線性化（上式右端項仍采用簡單的Picard線性化），得到

其中s表示非線性迭代次數.在不會引起混淆的情形，下面對非線性迭代有時將忽略時間上標.

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

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

對應的線性化方程（右端項仍簡單線性化）為

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

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

圖3 跳段導致非物理解Fig.3 Unphysical solution due to discontinuity in EOS

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

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

由式（12）解出Tn＋1（s＋1）后，判斷是否發生跳段．如未發生跳段，取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））．

這里算法的關鍵是判斷溫度是否經過跳段.本文設計如下的跳段判斷：當小于0，或者大于某個大于1的常數（通常取為5），則判斷為經過跳段.這兩種情況分別對應于向下跳段和向上跳段，如圖5所示.

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

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

適用于狀態方程跳段問題的非線性迭代算法

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

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

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

4）計算e（Tn＋1（s＋1））和跳段條件（14），自動判斷是否越過跳段：

① 未經跳段時，令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））；

5）判斷迭代是否收斂，如未收斂，轉3）；

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

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

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

迭代能量函數f（T）是分片光滑的連續函數，相比原間斷函數，迭代求解顯然更穩定.對于光滑問題，在適當條件下Newton迭代是局部二階收斂的，在某些特殊情形，Newton迭代全局收斂，參見文獻［4－5］.在文［4］中，當非線性函數的一階導數為Lipschitz連續時給出了收斂性結果.

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

3 數值算例

算例1 考慮一維熱傳導方程

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

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

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

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

整個過程計算到時間t＝1，時間步長為Δt＝10－4．

采用原方法（10）的計算結果如圖8所示.

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

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

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

采用本文提出的新算法的計算結果如圖9所示.

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

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

算例2 方程及初邊界條件與算例1基本相同，不同之處：計算區域為［0，0.1］，采用100個網格等分.離子能量為Ei（T）＝0.1T，電子能量的比熱cVe＝0.9，采用跳段的物態方程，在T0＝0.6處向下跳段－0.5，如圖10所示.

采用原方法（10）的計算結果如圖11所示.

由于計算區域較小，此問題的部分網格的溫度很快通過跳段點，直到時間0.4以后，所有的網格溫度都通過了跳段點.在時間0.4以前迭代次數較多，主要原因是：此時段部分網格的溫度T超過0.6而發生跳段，非線性迭代不易收斂.實際上，在溫度經過0.6的區域，可以看到由于跳段導致溫度出現不光滑的分布（見圖11（c）的底部）.圖11（b）表明整個系統的能量守恒且沒有漏失.

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

采用本文提出的新算法的計算結果如圖12所示.

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

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

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

4 小結

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

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

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

［2］ 徐錫申，等.實用物態方程理論導引［M］.北京：科學出版社，1986.

［3］ 常鐵強，張鈞，等.激光等離子體相互作用與激光聚變［M］.長沙：湖南科學技術出版社，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

國家自然科學基金（11171036，11001023）和中物院科學技術發展基金（2014A0202009）資助項目

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

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