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

第九章 存在相变的自由表面流动问题模拟——三相流T.L, MARINDeparment of Mining Engineering , University of Chile,Av. Tupper 2069 , Santiago, ChileE-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 ms 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的状态。

孔隙率定义如下：()()()()()1120m 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 是孔隙率的函数，根据以下方程定义：()()231C 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 函数同样用来定义连续相和非连续相间剧烈的物性变化。

这些物性包括密度，粘度和热传导率。

3．结果与讨论该模拟过程需要三个物理模型：不可压Navier-Stokes流动，对流和扩散传质过程，对流和热传导的传热过程。

求解该问题时用到了COMSOL Multiphysics 的多物理场功能来同时处理和求解三个应用模式。

3.1 下降液滴的固化本节介绍一个液滴在冷表面的流动和固化过程。

为简单起见，将其简化为二维模型。

整个系统用左下角在x=y=0处的无量纲矩形域来表示，底面边界温度固定，其余边界绝热。

液体和连续相（环境）的初始边界条件用以下初始水平集函数来表示：()00.5φ==(20)t这样就在模型数值域的左下角形成四分之一个圆。

注意这种情况下液相φ值为负，连续相φ值为正。

为简化起见，在该问题中我们不考虑液/气交界面处表面张力的影响，集中处理固化问题。

但是就像在水平集方法中提到的，加入这些影响因素也非常简单。

打开COMSOL Multiphysics模型导航栏。

按照表1中的步骤建立模型。

注意Navier-Stokes模型中边界2（底边）为滑移/对称边界条件。

虽然正常边界条件是无滑移的，但是该假设对处理边界上水平集“phi”函数的对流非常有帮助，更容易模拟液滴和底面边界的接触过程。

同时，由于液滴从底面向上固化，速度自然会计算为零，由于对固相进行求解，所以在这种情况下滑移/对称边界条件并不算太脱离现实。

如果使用无滑移边界条件，就需要用到水平集函数的再次初始化。

根据方程(20)，初始水平集函数“phi0”用空间坐标“x”和“y”来定义，液相为正（或者在液相区域）、气相为负，如图1所示。

Heaviside函数“H phi”用phi的负值来定义，所以在液滴内部为1，其它任何区域均为0。

液相和气相的物性用取决于水平集函数“phi”的Heaviside函数来定义。

根据方程(14)，用逻辑表达式定义孔隙率函数，该变量依赖于温度变量“T”和固相线、液相线的温度。

然后用该函数定义其它依赖于温度的物性。

图1 水平集函数初始条件的表面图该问题的新特性在于使用了映射网格。

它提供了棋盘形分布的四方形网格。

图2给出了不同时刻的解。

表面图只给出了负值范围内（液滴域）的水平集函数，同时等高线给出了液相线(T m+ε)、平均融化温度(T m)和固相线(T m－ε)温度的等高线，箭头图给出了速度矢量。

图2 在0.5，1，1.5和2s时的模型计算结果从图2中可以看出固化模型是如何作用的。

气相速度场和液相区域与预想结果相同。

在液、固相共存的多孔区域，速度分布仍然存在，但是更趋近于固相区，固相区计算出的速度几乎为零。

液滴中的热传递要比气相中的快很多，因为液相的热传导系数更高。

3.2 固化速率分析在整个模拟时间段内，固化过程已经发生，但是并不是整个液滴同时固化。

将fem结构导出，在COMSOL Script或MATLAB中运行“analysis.m”程序。

可以使用“postint”函数对 的Heavyside函数在空间域积分，算出不同时刻液滴总面积和完全固化区域的面积。

t=1:101;A=postint(fem,’Hphi’,’Solnum’,t);As=postint(fem,’Hphi*(1-flc2hs(T-Tm+e,0.001))’,’Solnum’,t);Fs=As./A*100;plot(0.02*(t-1),A,0.02*(t-1),As);figureplot(0.02*(t-1),Fs);在以上代码中，A是液滴总面积，它从一开始就应该保持常数0.196m2（质量守恒）；A s是固相面积，通过f的Heaviside函数和另一个Heaviside函数（flc2hs）在固相线(T m－ε)附近相叉得到；F s是液滴固化的百分比。