應用真實流向概念重建通量法則於靜直接模擬法速解尤拉方程式的研究

Abstract

近年來一項以模擬粒子為基礎來解非黏滯尤拉方程式的數值方法稱之為靜態直接模擬法（Quiet Direct Simulation, QDS）[Albright et al., 2002]。此方法利用高司¬-赫麥（Gauss-Hermite）積分集合式來取代馬克斯威爾與波茲曼（Maxwell-Boltzmann）速度分佈的表示。不僅可以輕易地處理流體中急遽變化的變數分佈（如大部分的計算流體力學一樣）更可以準確的仿效真實流的運行方式引入真實方向的守恒場，又因為具有地域性高（不需要大範圍考慮周遭的網格）的特點使得平行化的計算程式得以輕易達成。由於在每個網格內都擁有相同的離散速度分佈，因此對於計算時所要求的記憶體需求可減至最低。但所謂有一利必有一弊，因而衍生出數值擴散性（numerical diffusion）太強等問題，驅使Smith等人[Smith et al., 2009]發展出另一項以近二階QDS方法（QDS-2N）來改善相關問題以達到大幅提升數值的精確度。但是，在改善的過程中仍未真正考慮到流體真實流向的概念而有所遺憾。為解決此問題，本論文將對此數值方法提出改善方法並且用理論分析的角度進行探討。 在第一部份方面，就數值方法的觀點來說，本論文一樣是利用以近二階靜態直接模擬法來求解尤拉方程式，目的是發展一項根據流體的真實流向概念而進行的重建通量法則，我們稱之為QDS-N2法。不同於傳統的有限差分法是把通量計算存取於網格間的界面，相較之下，QDS 則是根據流體的真實流向來決定該網格的通量方式，其間，實際計算通過網格的流量並有效地存取於網格內使其達到完全遵循流體真實流向之計算法。在計算通量之前，我們會在每一個網格中引入數顆QDS 粒子，而此粒子根據瞬間加權擁有守恒場的變數值，並包含了質量、動量與能量等計算值。甚至守恒場的變數值變化情形（可以依據多項式表示）在同一網格中皆可被允許隨空間變化。雖然在一維流流場中，QDS-N2法的重建通量表示式與QDS-2N法的重建通量表示式是一模一樣，不過在流場延伸至多維度時兩種方法卻有所不同。為了呈現兩方法的優異性，在二維流的計算中，我們經由幾項不同的算例結果來驗證此新的流量重建法的優異性，並且在第二部份理論分析方面以解析解的方式探討與QDS-2N方法的差異。文中我們對於二維的研究比對項目包含了震波撞擊低壓近似氣泡問題、尤拉四震波的交互作用問題、尤拉四面波交互作用問題、三馬赫波正面衝擊水平台階和水平渦流擾動問題。並且對於模擬時所耗費的計算時間與更高階計算之效果做進一步的探討。根據結果顯示，在相較於QDS-2N法的情況之下，本文所提出的多維度重建法則在計算的精確度上確實有很大的提升。以得到相同數值精確度的標準而言，屏除因考慮流體真實流向而增加計算量所耗費的時間外，就水平渦流擾動問題的結果來說，相較於QDS-2N法，QDS-N2法在計算上所消耗的時間可以大幅減少。此外，誠如先前所提及QDS乃是一個高地域性的數值方法，非常有利於分散式叢集電腦上進行平行計算。因此，發展三維的QDS數值方法也將一併利用平行方式進行。而文中所有提及平行效率的研究與相關的平行計算，皆在台灣國家實驗研究院所屬國家高速網路與計算中心內所提供的ALPS叢集電腦中完成。就平行效率的研究方面，我們在計算大尺度的問題時使用高達256顆處理器分別完成0.5、2、12.5百萬的計算網格，根據strong scaling所得的結果顯示，得到的平行效率可達75%、68.5%與65.5%。另一項weak scaling的研究，比對理想值為1.0的平行效率，使用本方法所得到的效率可達1.2，其中處理器最高達到128顆，而平均每一顆處理器皆包含計算2萬個計算網格。 在第二部份理論分析方面，我們將所有有關質量、動量與能量的通量分別對QDS-2N和QDS-N2兩種數值方法進行解析解的分析推導。比較方式是針對兩運算法在流場中因不同變數值的改變於彼此之間所產生的通量差異。結果顯示二者通量值於低密度、低溫與高速的流域範圍中較易產生較大差異，而此範圍往往也是擴散波出現的地方。除此之外，本研究觀察到當兩數值方法在處理水平和垂直軸（x軸與y軸）的計算問題時彼此所產生的通量值差異甚小，反之在計算斜角流場時會相差較大。因此，可以根據本研究結果知道在未來處理模擬流場等問題時，得以輕易地判斷合適於問題的QDS-N2法或QDS-2N 方法而有效地取得所預期的結果，高精確度以利分析或是快速地得到流場趨勢。 對於主要的研究結果與未來研究方向的建議將總結於文末。

A particle-based quiet direct simulation (QDS) method [Albright et al., 2002] was invented to solve the inviscid Euler equation, in which the Maxwell-Boltzmann velocity distribution is enforced through the use of Gauss-Hermite quadrature integration without using any random number. It is a very fast Euler equation solver, which is deterministic with large dynamic range of flow properties like most conventional CFD methods, employs true-direction conservative fluxes for faithfully mimicking real flow motion, is highly localized (a small stencil) for easier parallelization and requires very low memory because the discrete velocities can be re-used in each cell. However, it is numerically very diffusive and has been extended to a nearly second-order numerical scheme by Smith et al. [2009] without really considering true-direction flux reconstruction. Thus, we intend to further address this problem from both numerical and theoretical viewpoints in this thesis. In the numerical part, a true-direction flux reconstruction of the second-order quiet direct simulation (QDS) as an equivalent Euler equation solver, called QDS-N2, is presented. Because of the true-directional nature of QDS, where volume-to-volume (true direction) fluxes are computed, as opposed to fluxes at cell interfaces as employed by traditional finite volume schemes, a volumetric reconstruction is required to reach a totally true-direction scheme. The conserved quantities are permitted to vary (according to a polynomial expression) across all simulated dimensions. Prior to the flux computation, QDS particles are introduced using properties based on weighted moments taken over the polynomial reconstruction of the conserved variables such as mass, momentum and energy. The resulting flux expressions are shown to exactly reproduce the existing second-order extension for a one-dimensional flow, while providing a means for true multi-dimension reconstruction.The new reconstruction is demonstrated in several verification studies. These include several two-dimensional test cases such as shock bubble interaction problem, an Euler-four-shock interaction, Euler-four-contact interaction, Mach 3 facing over a forward step, and the advection of a vortical disturbance. These results are presented, and the increased computational time and the effect of higher-order extension are discussed. The results show that the proposed multi-dimensional reconstruction provides a significant increase in the accuracy of the solution as compared to the previously developed QDS-2N method. We show that, despite the increase in the computational expense, the computational speed of the proposed QDS-N2 method is several times higher than that of the previously proposed QDS-2N method for a fixed degree of numerical accuracy, at least, for the test problem of the advection of vertical disturbances. As mentioned earlier, QDS method is intrinsically a highly localized numerical scheme, which makes it highly suitable for parallel computing on distributed-memory cluster machines using domain decomposition. With parallel implementation, an extension to three-dimensional QDS method is also demonstrated. The results show that the parallel efficiency, based on a strong scaling study, for a large-scale problem using 0.5, 2, and 12.5 million cells can reach up to 75%, 68.5%, and 65.5% with 256 processors respectively. In addition, the parallel efficiency, based on a weak scaling study, for a shock bubble interaction is 1.2, which the ideal efficiency is 1.0, up to 49 processors for 20,000 cells per processor. Note all the parallel performance tests were performed at the APLS cluster of National Center for High-Performance Computing, Taiwan. In the theoretical part, we have derived the analytical expressions of all the fluxes related to mass, momentum and energy in the two-dimensional QDS-N2 and QDS-2N methods respectively. Comparisons are made systematically between the corresponding fluxes in the two methods by varying flow properties. Results show that a large discrepancy of fluxes between these two methods occurs in the ranges of low density, low temperature, and high velocity.It is also interesting to learn that this range of gas flow often corresponds to an expansion wave region. Moreover, the fluxes using both methods are similar horizontal and vertical directions (x and y-direction), while large discrepancy is found in the fluxes going to the diagonal direction.With this observation, we can evaluate the accuracy of QDS-2Nmethod as compared to QDS-N2 method in the flow field, which may be important in deciding which method to be used for different flow problems. The major findings of the research with the recommendations for future study are summarized at the end of the thesis.