3.2.2　鑄液充型過程的數值模擬

鑄液在充型流動過程中會產生氧化、裹氣、散熱，以及壓力、溫度、速度、黏度波動等一系列化學和物理變化，因此，充型流動與鑄件質量密切相關。采用數值模擬實驗，不僅可以仿真鑄液在鑄型和澆冒系統中的流動狀態（包括流速、壓力等物理量的分布與變化），而且還可以仿真鑄液流動過程中的溫度分布及其變化，從而根據流速、壓力、溫度等物理量變化特征或規律優化澆冒系統設計，防止鑄液吸氣和氧化，減輕鑄液對鑄型的沖蝕，控制澆注不足和冷隔等缺陷的產生。鑄液充型過程的數值模擬主要涉及鑄液流動的自由表面處理、非穩態流場中的速度和壓力求解，以及流動與傳熱的耦合計算等多個方面。需要說明的是，目前比較成熟的充型流動數值模擬技術大都基于層流模型或修正的層流模型，這給流動充型的實際仿真帶來一定誤差，因為鑄液充型時的真實流動通常為非完全展開的紊流流動。盡管業界對鑄液充型過程的紊流流動進行了長期深入的研究，而且也取得了一些有益的成果，但是，鑄液充型的紊流問題仍然是當今流體力學和計算力學研究的熱點。

3.2.2.1　鑄液流動充型的數學模型

鑄液充型過程的流動屬于帶有自由表面的不可壓縮黏性非穩態流動。該流動過程包含質量傳遞、動量傳遞和能量傳遞，可用相應的數學模型（基本方程）加以描述。

（1）連續方程（質量守恒方程）

連續方程是質量守恒定律在流體力學中的具體體現。對于帶有自由表面不可壓縮黏性非穩態流體流動，有：

　　（3-29）

式中　　　D——散度；

u、v、w——流體在x、y、z三個方向上的流速分量。

方程（3-29）表明：對于不可壓縮流體的無源流動場而言，在充滿流體的流動域中任何一點，流速的散度等于0，即無源無漏，質量守恒。

（2）動量守恒方程

動量守恒方程是根據牛頓第二定律導出的黏性流體運動方程，該運動方程又稱納維-斯托克斯（Navier-Stokes）方程（簡稱N-S方程）。對于帶有自由表面不可壓縮黏性非穩態流體流動，動量守恒方程的表現形式如下。

　　（3-30a）

　　（3-30b）

　　（3-30c）

式中　g——重力加速度；

　ρ——流體密度；

　p——流體壓強；

　μ——流體運動黏度，μ=η/ρ；

　η——流體動力黏度。

方程（3-30）表明：由微元體內流體重力、微元體表面壓力和流體自身運動的動力（慣性力與黏性力之差）所產生的動量之和為零。其中：等式左邊第一項代表加速力，第二項代表慣性力；等式右邊第一項代表作用在微元體表面的壓力，第二項代表流體重力，第三項代表黏性力。

（3）能量守恒方程

能量守恒方程是熱力學第一定律在流體力學中的具體體現。對于帶有自由表面不可壓縮黏性非穩態流體流動，有

　　（3-31）

式中　c——流體比熱容；

　Q——流體內熱源；

　λ——流體熱導率。

方程（3-31）表明：流體流動引起的溫度變化主要由流體自身導熱和流體對流傳熱造成。其中：等式左邊第一項代表同溫度變化相關的能量，第二項代表同對流傳熱相關的能量；等式右邊第一項代表同流體自身導熱相關的能量，第二項代表同流體內熱源相關的能量。

3.2.2.2　求解流動充型問題的初邊值條件

（1）初始條件

鑄液充型流動的初始條件通常包括鑄液進入鑄型型腔或澆注系統瞬間（t=0）的初始速度和初始壓力。

①初始速度　初始速度一般是指t=0時的澆注系統入口處或鑄型內澆道入口處（如果不考慮澆注系統）的鑄液流動速度，該速度與澆注方式有關。

②初始壓力　初始壓力一般是指t=0時的鑄型內壓。當鑄型排氣良好時，鑄型內壓可設為零，否則應根據具體情況將初始壓力設為背壓。

③其他初始條件　如果在求解速度場和壓力場時需要耦合溫度場計算，則還應給出初始溫度條件。一般情況下，鑄液流動充型的初始溫度取鑄液的澆注溫度或澆注系統入口處的鑄液溫度。

（2）邊界條件

邊界條件指流體在其流動邊界上應該滿足的條件。鑄液充型流動過程中經常遇到的邊界條件有：

①自由表面速度　自由表面是指鑄液流動前沿的非約束表面（即不與鑄型壁接觸的表面）。在處理鑄液流動的自由表面時，動量守恒方程依然可用，但連續方程因流動域的改變而不再適用，因此必須仔細設置其速度邊界條件。如果鑄液前沿液/氣界面互不滲透，而且又滿足不發生分離的連續性條件，則界面處的法向流速相等。

②自由表面壓力　自由表面的邊界壓力一般由兩部分組成：一是因鑄型排氣不暢產生的阻礙鑄液前沿流動的空氣壓力（背壓）；二是鑄液前沿自由表面的張力。通常情況下，如果型腔內的氣體壓力已知（例如背壓P₀），則自由表面壓力的法向分量P_n和切向分量P_t滿足P_n=-P₀，P_t=0。

③約束表面的速度與壓力　約束表面一般是指鑄液與型腔壁接觸的表面（亦稱為鑄液/鑄型界面）。當鑄液沿固定型腔壁流動時，其法向和切向速度分別為零，這就是所謂無滑移邊界條件。當鑄液沿運動型腔壁（例如離心鑄造）流動時，液/固體界面處的鑄液流速等于型腔壁速度。當型腔壁多孔（例如砂型鑄造）且有鑄液穿越壁面時，則切向速度為零，而法向速度等于鑄液穿過壁面的速度。

④溫度邊界條件　是否給出溫度邊界條件一般由鑄液充型流動數值模擬是否需要耦合溫度場計算［即同時計算能量控制方程（3-31）］決定。耦合求解流動場與溫度場能夠更加準確地描述鑄液充型的真實流動狀態。

3.2.2.3　數學模型的差分離散

為了統一本章的差分格式，以方便學習和理解，現采用非交錯差分格式離散相應的連續方程（3-29）、動量方程（3-30）和能量方程（3-31）。

（1）連續方程的離散

　　（3-32a）

或

　　（3-32b）

（2）N-S方程的離散

以式（3-30a）為例，將方程中的各偏微分項用相應的差分項代替

將上述四式代入方程（3-30a），取Δx=Δy=Δz，并化簡成數值迭代形式

（3-33a）

同理，也有

　　（3-33b）

　　（3-33c）

式中　u、v、w ——坐標為（x，y，z）的流體微元當前時刻的流速分量；

　u'、v'、w'——同一微元在下一時刻的流速分量；帶_足標的u、v、w為當前時刻與微元（x，y，z）相鄰的其他6個微元的流速分量；

　　　　 Δt——時間步長。

借助式（3-33），可以利用當前時刻的微元流速分量u，v，w求得同一微元下一時刻的流速分量u'，v'，w'。

（3）能量方程的離散

將能量方程（3-31）等式兩邊各偏微分格式轉換成相應的差分格式

把上述兩式代入能量方程（3-31），令Δx=Δy=Δz，并化簡成迭代形式，有

　　（3-34）

式中　　　　a——熱擴散系數，；

　　T、T'——當前時刻與下一時刻的微元體溫度；

帶_足標的T——當前時刻與微元體（x，y，z）緊鄰的其他6個微元的溫度值（見圖3-2）。

3.2.2.4　SOLA-VOF求解法

求解帶有自由表面非穩態流動問題的關鍵在于確定自由表面的位置，跟蹤自由表面的移動，處理自由表面的邊界條件。鑒于SOLA-VOF法在解決上述三個關鍵環節上的優勢（例如：計算速度快、內存要求低等），目前，鑄液充型過程數值模擬多采用SOLA-VOF法。

SOLA-VOF法由兩部分組成，其中的SOLA部分負責迭代求解流動域內各微元體的速度場和壓力場，VOF部分負責處理流動前沿（自由表面）的推進變化。為了簡化求解過程，暫且不考慮鑄液充型流動中的熱量散失，即僅利用式（3-33）和式（3-32）求解鑄液的恒溫流動。

（1）SOLA法求解速度場和壓力場

SOLA法求解速度場和壓力場的基本思路如下。

①粗算速度值　以初始速度和初始壓力u₀、v₀、w₀、p₀或當前時刻的速度和壓力u、v、w、p為基礎，利用差分方程（3-33），粗算下一時刻的速度值u'、v'、w'。

②壓力、速度的修正及連續性校驗　以粗算速度值u'、v'、w'為基礎，利用式（3-32）校正壓力值，并將該校正值回代到式（3-33）中修正速度值。校正壓力值的目的是迫使鑄液充型流動的速度場滿足連續性條件［即質量守恒方程式（3-32）］，也就是通過修正壓力和修正速度將非零值的散度拉回到零或使之趨近于零。修正壓力和修正速度的計算過程是一個循環迭代過程（見圖3-4）。其中，每一步迭代計算獲得的速度逼近值，同時也是下一步計算修正壓力的速度校驗值；而每一步迭代計算獲得的壓力修正值又是下一步計算修正速度的初始值；如此循環，直到計算速度場滿足連續性條件或接近連續性條件為止。

由于壓力值校正過程不能破壞流體流動的動量守恒關系，因此，必須在N-S方程的約束下進行壓力值校正。計算壓力值校正的方法簡述如下。

分別用u'+δu、p_x+δp取代方程（3-33a）中的u'和p_x，再減去未取代前的原方程，得

對于微元x+Δx，采用類似步驟也有

將δu改寫成迭代形式，并注意到δuⁿ=uⁿ⁺¹-uⁿ，于是可得上述兩式的迭代過程表達式：

　　（3-35a）

式中　n、n+1——第n次迭代計算和第n+1迭代計算。

同理，可推導獲得

　　（3-35b）

　　（3-35c）

將式（3-35）代入連續方程（3-32a），整理，得

或

考慮近似關系，于是有

　　（3-36）

式中　。

由于D是連續性校驗中計算出的散度值，校正壓力的作用是要抵消不等于0的散度，使D回到0，因此，式（3-36）前面有一負號。

當Δx=Δy=Δz時，式（3-36）變成

　　（3-37）

式（3-37）即為立方體微元校正壓力值的計算公式。將該式計算結果代入式（3-33）便可得到新的修正速度。

在實際計算中，常常應用松弛迭代法來提高數值解的收斂速度，即給式（3-35）右邊第二項乘上一個松弛因子ω，使之成為

　　（3-38）

可以根據情況調整ω的取值，以改變迭代收斂的速度。一般，0<ω<2；其中，ω∈（0，1）為低松弛，ω∈（1，2）為超松弛。

（2）VOF法處理流動前沿的推進變化

利用VOF法處理鑄液充型流動前沿推進變化的基本原理是：為流動域中的每一個微元定義一個標志變量F，用以跟蹤流動前沿的推進；其中，流動前沿微元被定義成已填充區域與未填充區域之間的邊界微元，而標志變量F被定義成微元內的流體體積與該微元容積之比，稱為體積函數。由于微元內的流體凈體積由相鄰微元穿過界面的流量（即流速乘以時間）積累獲得，因此，F=0表明該微元為空，處于未填充域；0<F<1表明該微元部分充填，處于流動前沿；F=1表明該微元已充填結束，處于流動域內。

從上述基本原理可知，確定自由表面的移動，需要求解體積函數方程

　　（3-39）

式中　　F——體積函數，F=V_flow/V；

V_flow——微元內的流體體積；

　　V——微元容積。

邊長為Δx的立方體微元在Δt時間段內的體積函數變化可由下式計算獲得

　　（3-40）

實際上，式（3-40）表示由6個相鄰微元供給微元（x，y，z）的流體凈體積。顯然：當ΔF>0時，流入微元（x，y，z）的流體量大于其流出量，意味著該微元內的流體體積增加；ΔF=0，流入量等于流出量，微元內的流體體積不變；ΔF<0，流入量小于流出量，微元內的流體體積減小。若出現ΔF>1，則表明微元（x，y，z）已經填充滿，并由此產生新的邊界微元。

在求解計算流動場之初，令澆注系統入口處或內澆道入口處（若不考慮澆注系統填充）微元層的標志變量F=1，充型區（包括澆注系統與型腔或只含型腔）內其他所有微元的F=0。此外，作為流體填充源的澆注系統入口處或內澆道入口處的微元層還必須賦予非0的初始速度值。一旦利用SOLA法迭代計算獲得當前時刻微元的u、v、w、p值，就將其速度解代入式（3-40）計算體積函數的變化，然后根據變化的積累值判斷流動前沿（自由表面）微元的狀態與位置，并重新處理和設置其邊界條件。

3.2.2.5　流動與傳熱的耦合計算

嚴格說來，高溫鑄液充型流動過程總伴隨有熱量的散失，特別是鑄液在金屬型中流動時（例如金屬型鑄造、高壓鑄造、低壓鑄造等），其熱量損失速度將比較快。如果因熱量損失導致溫度下降過大，則會影響鑄液的充型流動，最終可能在鑄件中形成諸如冷隔、欠澆等質量缺陷。由動量守恒方程（3-30）和能量守恒方程（3-31）可知，鑄液充型流動的溫度變化將影響鑄液的熱熔、密度、熱導率，以及黏度和流速等物理量，進而改變鑄液充型的流動形式與流動狀態。

流動與傳熱的耦合計算實際上就是將能量方程（3-31）納入鑄液充型流動過程一道求解，即利用每一時刻（增量步）的鑄液流速更新鑄液溫度，再利用更新后的鑄液溫度修正動量方程和能量方程中的相關物理量（如鑄液黏度μ、密度ρ、比熱容c和熱導率λ），為下一時刻（下一增量步）迭代計算速度場和壓力場準備參數。

3.2.2.6　充型流動數值計算的穩定性條件

同凝固過程數值計算穩定性條件的處理方式相同，就是怎樣選取充型流動控制方程（3-32）～方程（3-34）中增量計算的時間步長Δt。時間步長的選取一般應考慮以下幾個因素：

①在一個時間步長內，鑄液流動前沿充滿的微元數不超過一個（一維流動）或一層（二維和三維流動），于是有

　　（3-41）

②在一個時間步長內，動量擴散不超過一個或一層微元，由此得

　　（3-42）

③在一個時間步長內，表面張力不得穿過一個以上或一層以上的網格微元，于是有

　　（3-43）

④如果考慮權重因子a，則有

　　（3-44）

最終，可根據下式選取滿足數值計算穩定性條件的最小時間步長Δt

　　（3-45）