相對論氣體動力論之數值方法（2/2）
計畫編號：NSC 89-2612-E-002-006
執行期限：89 年 8 月 1 日至 90 年 7 月 31 日
主持人：楊照彥 國立臺灣大學應用力學所
執行人：蔡尚熹 國立臺灣大學應用力學所
一、中文摘要 本文分為兩部分：第一部份根據熱力學第二定 律，推導在無外力場作用下達平衡時，滿足具相對 論效應之波茲曼方程式之微觀分布函數。第二部份 係利用由第一部份推導得到的平衡態微觀分布函 數，求其對應之巨觀物理量及所滿足之守衡律方程 式，並以「通量分離法(kinetic flux vector splitting method, KFVS)」求其數值解。Abstract
The aim of this work is two fold. First is to derive the equilibrium distribution function (EDF,fE)
for the relativistic Boltzmann equation in the absence of external field by the maximization of entropy and the principle of variation. The second purpose is to calculate the corresponding macroscopic quantities and to analyze numerically the governing equations in conservation law form by the kinetic flux vector splitting method.
二、緣由與目的
相對論氣體動力學主要的應用領域包括天文 學、電漿物理(plasma physics)與高能物理等方面， 如 超 新 星 爆 炸 所 產 生 的 強 震 波 、 電 磁 音 震 波 (magnetoacoustic shock wave)、重離子碰撞、高能 分子束碰撞等問題，由於其速度均極大，因此必須 考慮相對論效應。在無外力場作用下可以僅考慮平 坦 的 時 空 結 構(flat or Minkowski space-time structure)，以狹義相對論(special theory of relativity) 分析物理系統的相對論效應，亦即度量張量gµν為 (1) ( − − − = 1, 1, 1, 1 gµν diag ) 此時流體力學的基本方程式仍與古典力學的控制 方程式(如 Navier-Stokes 方程式)十分一致，均保持 雙曲型守衡律方程組形式，但方程式間的關係較為 複雜，且各物理量間又有所限制，因此計算較為困 難。而在外力場存在下，如重力場、電磁場等，則 需考慮度量張量因外力場造成的變化，亦即所謂 「時空曲率」或空間彎曲的問題。此時除了流體 力子間的交互作用，即使流速不大，仍需考慮粒 子與外力場間的交互作用及場方程式，方能正確 而完整的描述流體在外力場中的運動模式。 本文延伸古典氣體動力學中的微觀模型，假 設每一氣體分子均遵守相噸論的運動原則，應用 波茲曼方程式探討在平衡態下氣體分子的運動行 為，根據在平衡態下熵函數達最大值之熱力學定 律，推導平衡態之微觀分布函數。 數值方法部分則推導含相對論效應之通量分 離法，並比較不同初使條件所造成的相對論效應 及改善方式。 三、理論與數值方法 (一) 無外力場作用下之微觀分布函數 無外力場作用達平衡時粒子密度(n)、能量密 度( )、壓力( )與熵函數(S)均為定值，且可由 分布函數 ε p E f 求出： =
∫
E 3 n f d p (2) =c∫
f p d pE 0 3 ε (3) =∫
30 ∆ 3 E c d p p p p µ ν µν p f (4) ( ) =∫
−1 E 3 S lnf f d p ) (5) 式(3) 中 為 光 速 , 為” 動 量 - 能 量 四 維 向 量 ” (momentum four-vector) ， 為 投 影 算 子 (projection tensor)，定義為 c pµ ∆µν (6) ( − − − ∆µν =diag 0, 1, 1, 1 根據熱力學第二定律,達平衡時系統之熵函 數 達 最 大 值 ， 亦 即 可 以(2) 與 (3) 為 限 制 條 件 (constraints)重新定義熵函數如下： 1( )
(
(
)
= − − − − −∫
∫
)
∫
3 3 0 3 1 E E E E S d p lnf f a n f d p b ε c p f d p (7) G 則對(7)取變分(variation)應為零： ( ) =∫
3 E + + 0 E =0 S d p lnf a cbp f δ δ 由此可得fE為 ( = ⋅ − 0 E ) f A exp bcp (8) 此處常數 可由直接由(2)求得。若考慮理想氣體 (ideal gas)則其壓力與密度滿足 A (9) = B p nk T 其中 為波茲曼常數，T 為溫度。因此 可由(4)與 (9)求得。由此可得無外力場作用下相對於靜止座標 系滿足具有相對論效應之波茲曼方程式的微觀分 布函數 B k b E f : ( )(
− = 3 0 2 2 4 E B n b cp k T f exp z z K π)
(10) 式中b= 1kBT， = 2 B c k T z m 。觀察(10)式中cp 項為粒子相對於靜止座標系之總能量，為一不變量 (invariant) ， 可 將 其 改 寫 為 相 對 於 以 四 維 速 度 0(
)
= 1, Uµ γ UGc ， =( )− − 2 1 2 1 U γ G 為羅侖茲因子 (Lorentz factor)，移動之慣性座標(MFU)之總能量， 其形式如下： (11) = 0 cp p Uµ µ 上式中使用加法原則(summation convention)。除了 該項外其餘均為與座標系無關之物理量。因此可將 (10)改寫為相對於 MFU 之形式如下： ( )(
)
− = 2 3 2 4 E B n b p U k T f exp z z K µ µ π (12) 無外力場作用達平衡態時，具相對論效應之波 茲曼方程式為： ∂ = ∂ 0 E E df f p dτ µ xµ = (13) 此處 為靜止座標之”時間”(proper time)， x 為四 維時空座標(space-time coordinates)。上式若對「碰 撞不變量」(collisional invariants)積分亦成立： τ µ ( ) ∂ = ∂∫
3 0 E 0 d p p p f x p µ µ ψ (14) 若令ψ 可分別得到粒子數守衡、動 量守衡與能量守衡方程式，可表為相對於 MFU (U U (p)=1, ,p pi 0 ( ) = , ,V W )之守衡律形式如下： (15) ∂tQ+ ∂1F+ ∂2G+ ∂3H=0 ∂ ∂ = ∂ t t, ∂ ∂ = ∂ k xk , k = 1,2, 3 其中Q 為巨觀之守衡量，F、G與 分別為三個 空間方向之通量向量(flux vector)。分別定義如下： H ( ) ( ) ( ) ( ) = = = =∫
∫
∫
∫
3 1 3 0 2 3 0 3 3 0 E E E E p f d p p p d p f p p p d p f p p p d p f p ψ ψ ψ ψ Q F G H (16) 上式中若對應的積分區域分別取為 ， − − − −∞ ≤ ≤ −∞ ≤ ≤ −∞ ≤ ≤ 1 2 3 : 0 : 0 : 0 p p p F G H + + + ≤ ≤ ∞ ≤ ≤ ∞ ≤ ≤ ∞ 1 2 3 : 0 : 0 : 0 p p p F G H則得三維之分離通量(split flux vector) 、G 與 。 ± F ± ± H (二) 數值方法 本文採用通量分離法解一維具相對論效應之 氣體動力學問題。令U =(U, 0, 0 並分別求 、G 與 H 之分離通量則(16)可簡化為： F ) G (17) ± ∂tQ+ ∂1F =0 其中 ( ) ( ) + = + − 2 2 n e p U e p p γ γ γ Q (18) ( ) ( ) ( ) =±∞ ± = ± ± ± = ± + + + = + + +
∫
∫
11 1 2 1 0 0 1 0 2 2 2 01 0 2 2 10 0 2 2 2 2 p E p p p d p dp f p nU N e p U p UT e p U T ψ γ γ γ γ γ γ F K K K (19) 2上式中 ( ) ( ( ) ) ± = ± + − 2 0 2 1 2 k z exp z N cn x K z z K K (20) ( ) ( ( ) ) ± = ± 2 + + − 0 2 0 2 3 3 2 k n x z z exp z T K z z β K K (21) modified Bessel function of the second
kind ( )= 2 z K 四、結果與討論 本文第一部份由熱力學第二定律出發，以變分 法求平衡態下氣體的微觀分布函數，因此必滿足最 大熵值的條件，較傳統以H-theorem 推導微觀分布 函數的方式為簡潔，並對推導有外力場存在之分布 函數提供一個值得參考的方式。 第二部分則以通量分離法分析黎曼問題，在數 值研究中發現相對論效應主要來自羅侖茲因子與 Bessel 函數：低速運動時羅侖茲因子約與 Bessel 函 數可以下式近似 ( ) − ≈ ≈ 1 2 z n K z e z γ π (22) 相對論效應不明顯，計算結果與古典氣體動力學之 結果相當一致；僅在高速運動下相對論效應方才顯 現，此時羅侖茲因子仍可直接計算，然而 Bessel 函數則需特別處理。本文以四點求積分法求Bessel 函數之近似值，並與(22)之低階近似比較，如附圖。 圖(1)與(2)為測試本研究數值方法之兩算例，均 假設為理想氣體，其氣體之gamma value=1.4，總 計算時間= 0.48，空間範圍 ，圓圈符號 為以(22)近似之 Bessel 函數，紅色線條為四點求積 分法修正近似之 Bessel 函數。圖(1)為低速流場， 相對論效應較不明顯，計算條件如下：左右初始狀 態分為 ≤ ≤ 0 x 1 左：(n v pL, ,L L)=(1,0,1) 右：(n v pR, R, R)=(0.125, 0, 0.1) 總格點數=100，圓圈之 CFL number=0.1，紅色線 條之 CFL number=0.6。圖(2)為高速流場，相對論 效應較明顯，計算條件如下：左右初始狀態分為 左：(n v pL, , )=(10,0,40/3)L L 右：( , , )=(1.0,0.0,2/3*10 )-6 R R R n v p 總格點數=200， CFL number 均為 0.1。由計算結 果可知在相同的格點數條件下，以高階近似的 Bessel 函數除了較佳之震波解析結果，也可使用較 大之CFL number。 五、致謝 本研究由國科會 NSC 89-2612-E-002-002 部 分經費支持，謹此致謝。 六、參考文獻
1. Ray D’Inverno, “Introducing Einstein’s
Relativity”, 1998.
2. S. R. de Groot, W. A. van Leeuwen, and Ch.
G. van Weert, ”Relativistic Kinetic Theory: Principles and Applications”, 1980.
3. W. M. Eppard, B. Grossman, “A
multi-dimensional kinetic-based upwind solver for the Euler equations”, AIAA-93-3303-CP.
4. G. A. Bird; “Molecular Gas Dynamics and the
Direct Simulation of Gas Flows”, 1994.
5. S. Chapman, T.G. Cowlings, “The
mathematical Theory of Non-Uniform Gases”, 1970.
6. Z. Kopal, ”Numerical analysis”, 1961
7. 陳旻宏,”相對論氣體動力學之數值研究”,1995
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Number D ensit y 圖(1.a).密度分布 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.1 0.2 0.3 0.4 0.5 Ve lo ci ty 圖(1.b).速度分布 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Pr es su re 圖(1.c).壓力分布 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.8 1.6 2.4 3.2 4.0 En ergy 圖(1.d).能量分布 0.0 0.2 0.4 0.6 0.8 1.0 0.0 2.0 4.0 6.0 8.0 10.0 Num be r Den si ty 圖(2.a).密度分布 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 Ve lo ci ty 圖(2.b).速度分布 0.0 0.2 0.4 0.6 0.8 1.0 0.0 4.0 8.0 12.0 16.0 Pr es su re 圖(2.c).壓力分布 0.0 0.2 0.4 0.6 0.8 1.0 0.0 10.0 20.0 30.0 40.0 50.0 En ergy 圖(2.d).能量分布 4
