9 存在相变的自由表面流体问题模拟

第九章 存在相变的自由表面流动问题模拟——三相流

T.L, MARIN

Deparment of Mining Engineering , University of Chile,

Av. Tupper 2069 , Santiago, Chile

E-mail:tmarin@ing.uchile.cl

本章介绍了流体自由表面凝固过程的固定网格数值模型，并采用水平集方法来描述自由表面的流动。V oller 和Prakash 对此方法进行了改进，考虑了液相凝固问题，包括对流和传递过程。在该方法中，液体物性与温度相关，通过设定0到1之间的孔隙率使得体积力依赖于孔隙率和温度的变化，并通过对Navier-Stokes 方程的修正来模拟液相或者固相。此外，在传热方程中使用改进的热容表达式来计算融化潜热。

1．绪论

液体在凝固时主要通过对流来传热，由于涉及到移动界面问题，所以通常难于模拟。此外，对流体有效的控制方程对固体不再适用。而且，如果问题中包含自由表面，那么流动将变得更加复杂，需要考虑相变或者凝固问题，这时都需要跟踪相界面的变化。

目前有很多种处理液体凝固的计算方法，考虑封闭空间中存在导热和自然对流的情况，V oller 和Prakash 对此作了总结[1]。进展之一为采用变形网格来处理液体－固体界面的移动问题。也可用固定网格以及焓变随温度的变化来处理该问题，这种方法的特点是较为简便，但在研究零速度封闭容器中的固体时会存在问题。对于给定的计算单元，可以简单的设置速度为零，或者通过设置粘度为潜热容[2]函数的方法来实现[3]。在这类情况下，当潜热容接近于零时，粘性会很大，这样才能模拟固体物性。V oller 和Prakash[1]研究了在一定温度范围内发生相变的情况，他们将流体描述为多孔介质，通过在Navier-Stokes 方程中引入源项来模拟相变过程，用于研究速度为零的固体情况。

存在相变的自由表面流动问题的研究难点在于其实验和计算验证均很困难。Pasandideh-Fard 等[4]采用流体体积元方法(VOF)对锡液滴和水滴热表面上[5]的凝固进行了实验和计算研究，对自由表面及过程中的焓变进行跟踪。他们的模型引入了液滴接触角以及液滴底部界面传热系数，该系数由实验测试得到并作为模拟的边界条件。

2．控制方程

纯净流体凝固过程由熔点决定。在液体冷却过程中，一旦达到这个温度，在温度继续降低之前液体开始释放相变潜热。但是在多组分体系中，相变存在一个温度范围，从固体开始出现时的液体温度开始，一直到最后一种液体凝固时的温度为止。在这种情况下，融化潜热在温度改变的同时不断释放。

相变期间固体组分(Fs)可以表示为温度的函数：

()()()()()021m m

s m m m T T T T F T T T T T T εεεεεε≥+⎧⎪+-⎪=+>≥-⎨⎪<-⎪⎩ (1)

其中，T 为系统温度，T m 为液相线和固相线的温度平均值，ε为液相线和固相线温度差值的一半。因此，固相和液相温度为：

S m T T ε=-

(2) L m T T ε=+ (3)

前面已经提到，总的体系热容H 由两部分组成，显热h 和潜热△H 。显热通过以下方程计算：

p h c T =⋅ (4)

潜热可以表示为温度的函数，根据之前对固体比例的定义：

()()()()()()10m s m m m L T T H T L F T T T T T εεεε≥+⎧⎪∆=-+>≥-⎨⎪<-⎩

(5)

导热和对流传热方程用系统温度形式表示为：

()0p p T c k T c Tu t

ρρ∂⋅+∇⋅-∇+⋅=∂ (6) 在这种情况下，释放潜热带来的影响可以通过对有效热容方程的重新定义来包含到热容项中：

pT p pH c c c =+

(7) ()()2,pH m c L T L fldc hs T T δε*=⋅≈⋅-

(8) 这里，()T δ*表示平滑delta 方程，在COMSOL Multiphysics 中由fldc2hs 函数建立。注意，这里需要对该函数在整个温度范围内积分，但是只有m T T ε=-到m T T ε=+的温度范围。只有在模拟纯净物质凝固过程的时候，ε才为零，此时()T δ*将变为真正的Dirac ’s Delta 。

通过这种方法，系统总热容可以通过对方程(7)积分得到：

()222111T T T pT p T T T H c dT c dT L T dT δ*=⋅=⋅+⋅⋅⎰⎰⎰ (9)

速度场和压力场的控制方程即Navier-Stokes 方程：

()()

()T u u u u u p F t

ρμρ∂-∇⋅∇+∇+⋅∇+∇=∂ (10) 0u ∇⋅= (11)

体积力F 包括重力和水平集方法处理的表面张力，同时F 也包含依赖于凝固过程的固体比例源项。体积力分量如下：

()x x F S x φσκδφ∂=+∂ (12)

()y y F g S y φσκδφρ∂=++∂ (13)

其中，σ是液体表面张力，κ是交界面曲率，φ是水平集函数，g 是重力加速度，S x 和S y 项表述如下。

为了将固化过程并入体积力中，可以将液相看作是一种多孔介质，孔隙率λ依赖于温度。全液相对应于孔隙率为1的状态，全固相对应于孔隙率为0的状态。孔隙率定义如下：

()()()()()

1120

m m s m m m T T T T F T T T T T T εελεεεε≥+⎧⎪+-⎪=-=+>≥-⎨⎪<-⎪⎩ (14) 源项定义如下： x S Au =-

(15) y S Av =-

(16) 这里A 是孔隙率的函数，根据以下方程定义：

()()2

31C A q λλ-=+ (17)

下面介绍一下这些源项的作用。当温度位于液相线以上时，系统处于全液态，源项取值为零，对Navier-Stokes 方程没有任何改动。多孔区域意味着温度处于液相线和固相线之间，A 值增大并影响瞬态、对流和扩散项，流动方程近似相当于多孔介质中的Darcy 定律。当温度进一步降低，孔隙率接近于0时（固相），该源项决定了所有其它源项，使得速度值趋于零（固相）。方程(17)中的常数C 和q 是任意选取的，依赖于求解的具体问题，通常C 取较大值而q 取较小值，避免当λ变为零时方程被0除。

水平集函数φ定义了两个初态相（例如，气和液）的交界面。但是当固化（或者反问题中融化）时会出现第三相，一般出现在两个初始相的其中之一，所以需要修改源项S x 和S y ，使其严格作用在某一相中。求解该问题非常简单，只需要将这几个源项乘以φ的Heaviside 函数即可，出于我们对相的选择，我们限制了这些源项的作用范围。这样，方程(10)中的体积力最终可以表示为：

()()x x F S H x

φσκδφφ∂=+∂ (18) ()()y y F g S H y φσκδφρφ∂=++∂

(19) 类似水平集方法，Heaviside 函数同样用来定义连续相和非连续相间剧烈的物性变化。这些物性包括密度，粘度和热传导率。