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）边界条件
边界条件指流体在其流动边界上应该满足的条件。铸液充型流动过程中经常遇到的边界条件有：
①自由表面速度 自由表面是指铸液流动前沿的非约束表面（即不与铸型壁接触的表面）。在处理铸液流动的自由表面时，动量守恒方程依然可用，但连续方程因流动域的改变而不再适用，因此必须仔细设置其速度边界条件。如果铸液前沿液/气界面互不渗透，而且又满足不发生分离的连续性条件，则界面处的法向流速相等。
②自由表面压力 自由表面的边界压力一般由两部分组成：一是因铸型排气不畅产生的阻碍铸液前沿流动的空气压力（背压）；二是铸液前沿自由表面的张力。通常情况下，如果型腔内的气体压力已知（例如背压P0），则自由表面压力的法向分量Pn和切向分量Pt满足Pn=-P0，Pt=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法求解速度场和压力场的基本思路如下。
①粗算速度值 以初始速度和初始压力u0、v0、w0、p0或当前时刻的速度和压力u、v、w、p为基础，利用差分方程（3-33），粗算下一时刻的速度值u'、v'、w'。
②压力、速度的修正及连续性校验 以粗算速度值u'、v'、w'为基础，利用式（3-32）校正压力值，并将该校正值回代到式（3-33）中修正速度值。校正压力值的目的是迫使铸液充型流动的速度场满足连续性条件［即质量守恒方程式（3-32）］，也就是通过修正压力和修正速度将非零值的散度拉回到零或使之趋近于零。修正压力和修正速度的计算过程是一个循环迭代过程（见图3-4）。其中，每一步迭代计算获得的速度逼近值，同时也是下一步计算修正压力的速度校验值；而每一步迭代计算获得的压力修正值又是下一步计算修正速度的初始值；如此循环，直到计算速度场满足连续性条件或接近连续性条件为止。
图3-4 SOLA-VOF法计算流程
由于压力值校正过程不能破坏流体流动的动量守恒关系，因此，必须在N-S方程的约束下进行压力值校正。计算压力值校正的方法简述如下。
分别用u'+δu、px+δp取代方程（3-33a）中的u'和px，再减去未取代前的原方程，得
对于微元x+Δx，采用类似步骤也有
将δu改写成迭代形式，并注意到δun=un+1-un，于是可得上述两式的迭代过程表达式：
（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=Vflow/V；
Vflow——微元内的流体体积；
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）