气泡在液体中运动过程的数值模拟

2 期 卢作伟等:
气泡在液体中运动过程的数值模拟 129
〔 12〕 式 ( 12) ( 13) ( 14) ( 15) 进行。 ( 3) 为满足连续方程 ( 10) , 按 S I M PL E 方法进行压力修正, 收敛 的条件是全流场的压力最大误差不得超过某一个值。 ( 4) 距离场按方程 ( 11) 进行推进。 ( 5) 距 ( 6) 重复步骤 ( 12) 到 ( 5) , 进行不定常时间的下一个时间步。 离场的重新初始化, 计算方程 ( 16) 。
( 14)
128
计 算 力 学 学 报
14 卷
由于在界面两侧物性参数的变化比较大, 所以为了能够顺利的进行计算, 我们将把物性参 数在界面层中进行平滑的过渡, 以保证控制方程组解的存在。 在一般情况下界面层都 是非常 薄的一层, 因而这种过渡只要比较平滑即可, 在这里我们选用正弦函数来作为过渡函数变化的 基础, 这样可以保证物性参数的一阶导数和它本身在全流场内都是连续的。 具体的计算公式如 下: 1 Υ> Α ) = Xb Xl ( 15) X (Υ Υ< - Α Υ ≤Α ) X + ∃X sin ( Π Υ2 Α X = (X b + X l ) 2 X l ∃X = (X l - X b ) 2 X l 上式中下标 b 表示气相的参数, l 表示液相的参数, X 表示任一流体物性参数, 如密度或粘性系 数等。 当 Υ场是严格的距离场时, 而且它的等零线代表界面, 那么界面的曲率将可以从 Υ场得到 严格的计算公式如下: ) = ( Υ ( 16) k (Υ Υ) = Υ 214 距离函数的重新初始化 我们要求在每一时间积分步 Υ都是严格的距离场, 但是从对流方程 ( 3) 计算所得的 Υ场显 然不能始终保持这一特性。 所以我们从距离场和 Υ = 1 互为充要条件出发, 引入下面的 2 方程 : H am ilton J acob i ) ( 1. 0 Υ Υ) Σ = sign ( Υ ( 17) → → Υ( x , 0) = Υ 0 (x ) → ) 是 Υ的符号函数, Υ 式中: sign ( Υ 在每一完 0 ( x ) 是每一时间步由对流方程 ( 3 ) 得到的 Υ的分布。 整积分步后, 用方程 ( 17) 重新调整 Υ场, 这样得到的 Υ场满足 Υ = 1, 并作为下一个时间步 的距离场。 215 算例的参数与边界条件 在实际计算中气泡的初始形状是球形或其他轴对称形, 因此我们在柱坐标系下进行数值 求解, 并假定气泡在半径足够大的圆柱筒中运动。 采用交错网格来离散求解区域, 对流项用二 阶迎风格式, 时间方向是一阶精度的。 在圆柱边界面上采用无滑移条件, 在垂直于气泡运动的 上下计算域界面上采用第二类边界条件。 我们计算了两类问题, 第一类是文献 〔 11〕 实验的气泡在糖溶液中运动, 溶液的物性参数的 范围是粘性系数从 01082 ( Pa ・ s) 到 218 ( Pa ・ s) , 密度从 131410 ( kg m 3 ) 到 1390 ( kg m 3 ) , 表 面 张力系数从 010769 (N m ) 到 0108 (N m ) 。 另一类就是水中的气泡运动, 粘性系数为 01001137 ( Pa ・ s) , 密度是 100010 ( kg m 3 ) , 表面张力系数是 010728 (N m ) 。 气泡的物性参数 都按空气的参数计算, 粘性系数是 010000178 ( Pa ・ s) , 密度是 11226 ( kg m 3 ) 。 显然在我们的 算例中, 界面处物性参数的变化最小是几百倍, 最大可达几十万倍, 因而有效的处理界面对计 算的顺利进行是至关重要的。 计算的步骤如下: ( 1) 给定初始场 ( 在我们的计算中初始场是静止场) , 特别要给出初始的 距离函数场, 要求它严格符合距离场的性质。 ( 2) 求解动量方程 ( 9) , 其中有关参数的计算按公
P +
→
) + F T + Ρk n ∆Α( Υ
→ →
→
→
( 9) ( 10) ( 11) ( 12) ( 13)
→ →
U= 0
Υ t + U
→ →
→
Υ= 0
→ →
T = 2Λ D
) = ∆Α( Υ
→ →
( 1. 0 + co s ( Π )) 2 Α ΥΑ Υ < Α 0 其它
→ → →
式中: U 是流体速度, Θ是流体的密度, Λ 是流体的动力粘性系数, T 是粘性偏应力张量, D 是应 变率张量, F 是体积力, Υ是距离函数, Ρ 是表面张力系数, k 是界面的曲率, n 是界面的法向矢 ) 是选择函数, P 是压力。 量, ∆Α( Υ 这是 Α值的选取将关系到界面捕捉的精度和计算的顺利进行, 我们按下式估算: Α= m ax ( ( 1 - 2) ∃ y , 式中 ∃ y 表示空间网格的大小。 213 物性参数的过渡与界面曲率的计算 ΧΘ g)
Α
- Α
∫
Υ →±Α
∫
- Α
DU Θ Dt
→ →
Α
n d Υ=
∫
(- Α
P + F+
→
) Ρk n ) T + ∆Α( Υ
Α →
→ →Biblioteka →→n dΥ
或
Α
∫
(- Α
P +
) Ρk n ) T + ∆Α( Υ
→ →
→
→
n d Υ=
∫
- Α
(Θ
DU Dt
F)
→
→
n dΥ
当 Α→ 0 时, 上式右端等于零, 左端经积分后得: ∃ P + ∃Σn + Ρk = 0 ( ) 显然这与前面界面的动力学条件 6 相吻合。 界面运动条件 ( 2) ( 3) 和动力学条件 ( 4) ( 5) 由模型 的连续性可以得到保证。 利用界面层模型, 全流场可用如下的统一控制方程组: Θ DU D t = →
摘 要 本文用数值方法预测气泡在液体中的非定常运动。 运用位标函数进行界面的 隐含跟踪, 并且与有限体积法相结合, 构成一种可行的计算方法。 本文方法允许在界面 处存在很大的物性差, 而且较容易将表面张力引入控制方程。 我们对气液两相流中单个 气泡的运动进行了数值计算, 得到了与实验结果符合很好的数值结果。 关键词 气泡运动; 位标函数; 有限体积法 分类号 TV 312; O 353
2 数值方法的描述
211 基本方程
我们假定气液两相之间没有化学或物理化学作用, 同时都可以作为不可压缩牛顿型流体。 对于有界面的流体运动, 我们可以写出它们的控制方程和边界条件如下: Θ DU D t = → →
P + U= 0
→ → → → →
T + F
→ →
→
( 1. a ) ( 11b )
→
式中 Θ是密度, U 是速度矢量, P 是流体压强, T = 2Λ D 是牛顿型流体的偏应力张量, F 是体积 → 力。 除了壁面无滑移条件外, 在气泡的界面上提供边界条件。 设 Υ( x , t) = 0 是流动界面, 则以下 三个条件都应当得到满足: ( 1) 界面的守恒方程 → ( 2) Υ Υ= 0 t + U ( ) 2 界面上的运动学条件 → → ( 3) V1= V2 ( 3) 界面上的动力学条件 ( 4) Σn1 = Σn2 ( 5) Σt1 = Σt2 ( 6) ∃ P + ∃Σn + Ρk = 0 ∃ P 是界面两侧压力差, ∃Σn 是界面上的法向偏应力, Ρ 是表面张力系数, k 是界面上的曲率, Σt1、 Σt2、 Σn1 和 Σn2 分别是界面两侧两个方向上的切应力, 三式分别表示界面上切向和法向应力平衡。 212 距离函数和界面层 如果我们用以上的基本方程和边界条件进行数值计算, 在求解不定常问题时将不可避免 地遇到运动边界问题, 特别是当涉及到界面破碎和弥合时, 处理起来将非常困难。 因而对界面 进行隐含的追踪处理是较为可行的选择。 → 从方程 ( 2) 可以看到, 如果我们把界面方程 Υ( x , t) = 0 扩展为一个在全流场存在的标量场 → Υ( x , t) , 让它的零值面代表界面, 其它位置上的值是该点到界面的距离, 这样就形成一个距离 函数场。 为使 Υ场是距离函数场, 须满足 Υ = 1, 在后文我们将叙述如何在计算过程中满足 ) , 该层内表面张力 这一条件。 同时为了抑制数值扩散, 我们假设界面是一厚度很小的薄层 ( 2Α 和物性参数是分布的, 我们称这一层为界面层。 如图 1 所示。 界面层内的流动满足下面的方程: Θ DU D t = →
收稿日期: 1996204220; 修改稿收到日期: 1996209226
126
计 算 力 学 学 报
14 卷
果一致, 但严格来说他们都不是数值模拟真实的气泡运动, 计算中的气液密度比或粘性比都没 达到真实的物性参数值。 我们数值模拟了空气泡在水中和糖溶液中的运动, 并以 B haga 和 〔 11〕 W eber 1981 年气泡在糖溶液中运动的实验结果作为对比资料。
第 14 卷第 2 期 1997 年 5 月
计 算 力 学 学 报
CH I N ESE JOU RNAL O F COM PU TA TLONAL M ECHAN ICS
V ol. 14 N o. 2
M ay 1997
气泡在液体中运动过程的数值模拟
卢作伟 崔桂香 张兆顺
( 清华大学工程力学系, 北京, 100084)
P +
→
) + F T + Ρk n ∆Α( Υ
→ →
→
→
( 7. a ) ( 7. b )
U= 0
) 为确定表面张力的分布函数。 方程右端第三项是表面张力在界面层中的分布项, 式中 ∆Α( Υ 它
2 期 卢作伟等:
气泡在液体中运动过程的数值模拟 127
图 1 界面层模型的示意图
Α
) d Υ= 1 和 li m ∆Α( Υ ) = 0。 满足 ∆Α( Υ 计算中的具体表达式为: ) = ( 1. 0 + co s ( Π ) ) 2Α ( 8) ∆Α( Υ ΥΑ 可以证明界面层方程 ( 7) 是渐近收敛于界面的动力学条件。 即当界面层的厚度 Α→ 0 时, 界面层模型的控制方程 ( 7) 应当趋于边界条件 ( 2) ( 3) ( 4) ( 5) ( 6) 。 为此我们对方程组 ( 71a ) 沿 界面法向积分。
1 引 言
气泡运动是气液两相流动中的基本问题, 它是一种具有气液界面的非定常流动。 如果在数 值计算中不对界面附近的物性参数变化进行有效的处理, 而用通常的差分格式进行数值计算, 如差分格式精度低, 过度的数值扩散往往会把界面淹没掉; 但如应用精度较高的差分格式计算 时, 稳定性比较差, 尤其当存在物性参数的急剧变化时, 通常会导致数值振荡或得不到收敛的 解。 因此发展一种适合于气泡运动的计算方法是非常必要的。 准确跟踪气液界面是气泡运动数值模拟的关键。 目前跟踪界面的方法主要有以下几大类, 第一类是面跟踪的方法, 它是把一些标志点放在界面上, 通过标志点的运动来模拟界面的运 动〔1〕 , 第二类是体跟踪的方法, 它是通过计算各离散网格中各成分所占体积百分比来计算界面 的准确位置〔2〕 , 第三类是动网格的方法, 就是把界面作为一离散的网格线或面来处理〔3〕 , 第四 类是位标函数方法, 就是用位标函数来隐含界面, 通过位标函数场的变化来模拟界面的运 〔 4〕 动〔10〕 。 上述方法除位标函数方法外都是显式跟踪界面的思想, 因而它们的计算量一般都比 较大, 而且在界面附近的处理也很复杂。 位标函数是用隐式的方法进行界面的跟踪, 它能够更 容易地与动量方程的求解合在一起, 从而简化求解过程。 本文采用流场中任意点到界面的距离 作为位标函数进行界面的捕捉。 在现有的气泡运动数值研究中, 都对气体的物性作了简化。 例如, 1976 年 Young ren 和 〔 5〕 对气泡运动的整个过程进行数值模拟, 最早是 R ysk in 和 A crivo s 研究了势流中的气泡运动 。 〔 6〕 L ea l 用有限差分的方法与动边界相结合的思想研究了气泡平衡时的形状, 他们只考虑液体 相, 忽略了气相物质的性质。 在以后的研究中人们开始考虑气相的存在, 也就是说把它作为两 相系统进行计算, D andy 和 L ea l〔7〕在 1989 年, A nderson〔8〕在 1985 年, B aker 和 M oo re〔9〕在 1985 年, 分别对气液两相有有限密度比与零密度比的问题进行了研究; U nverd i 和 T ryggva 2 〔 10〕 son 于 1992 年也发表了对气泡运动的数值模拟结果。虽然以上的数值模拟定性地和实验结