编辑: Cerise银子 | 2019-07-06 |
5 5 t ρ uj +
5 5 xi (ρ ui uj) = -
5 p
5 xj +
5 5 xi μ(
5 uj
5 xi +ρ gj) (1) 式中 ρ = α g ρ g +α lρ l (2) μ = α gμg +α lμ l (3) gj 为j向重力加速度;
p 为压力. 1.
2 质量守恒方程 气液两相间界面位置的确定须通过求解由连续 性方程导出的相体积分数方程完成 ,即55t(α k ρ k) + ui
5 5 xi (α k ρ k) =
0 (4) 对不可压缩流体 ,应用连续性方程 ,可得散度形式连 续性方程
5 5 t α k +
5 5 xi (α k ui) =
0 (5) 式中 :α k 代表 k 相的体积分数 ( k 为g时表示气相 , k 为l时表示液相) ;
ui 代表 i 方向速度. 式(5) 的求 解只对分散相进行 ,主相体积分数按照下式求解 α g +α l =
1 (6) 1.
3 压力方程 将连续性方程和动量方程构成的求解方程组在 非均分(δx ≠ δ y) 交错网格上进行差分 ,对于非稳态 加速项 ,时间应用一次向前差分 ,速度计算应用显式 近似 ,即将前一次迭代的计算结果作为本次计算迭 代的初值. 对流加速项应用 QU ICK 差分格式 ,粘性 应力项应用中心差分. 如此对连续性方程和动量方 程进行差分分解 ,得到的差分求解方程并不封闭 ,还 需要与之对应的压力方程. 本文对于速度与压力的 耦合采用 SIMPL EC 算法. 由给定压力 ,按次序求解 u 、 v 动量方程 , 得到的速度场未必满足连续性方 程 ,因此须对压力做出修正. 由修正的压力去改进速 度 ,得出在这一迭代层次上满足连续性方程的解 ,然 后以此新的速度场开始下一层迭代. 由连续性方程 推导得到的压力修正方程如下式所示 ,其详细推导 过程参见文献[6 ]. Δpi , j aP - Δpi -
1 , j aE - Δpi +1 , j aw - Δpi , j -
1 aS - Δpi , j +1 aN = b (7) 式中 : aP = aE + aW + aS + aN ;
b 为源相;
压力修正 Δp 对压力 p 按照 pn +
1 = pn +Δp 进行修正 , n 表示 时间层. 当界面移动时 ,应用式 (5) 确定界面位置 ,对液 体控制容积 ,应用上述压力修正方程求解动量方程 , 对气体控制容积 ,不求解动量方程 ,认为Δp = 0. 对 于气液界面应用垂直界面边界条件代替 ,即5τn
5 n =
0 (8) 此处 n 为界面法向 ,τ n 为法向应力 ,在界面上忽略 粘性 ,则ps = pg + K σ 式中 : ps 为界面上液体侧的压力;
pg 为气泡压力;
K 为界面曲率;
σ为表面张力. 为得 到界面控制容积的压力,应用线性插值[7 ] ,如图
1 所示 ,得pn+1 i , j = η i , j ps + (1 - η i , j) pl (9) 由于压力方程为Δp 的线性方程 ,故式 (9) 可表 示为 Δpi , j + (η i , j - 1)Δpl = η i , j ps + (1 - η i , j) p n l - p n i , j (10) 图1表面压力插值方法 式中 : pl 为界面液体 侧全控制容 积的压力;
η i , j = d/ dc. 由式 (8) 、 式(10) 可得到压 力修正 Δp. 由压力 修正可得到速度修正 值 ,平移界面 ,改进压 力 ,开始下一时层迭 代. 求解过程为防止 发散 ,本文对界面插值应用施主Ο 受主格式进行离 散 ,如图
2 所示. 在界面附近 ,这种格式识别流体从 控制单元界面流出的相单元为施主单元 ,流体流入 另一相单元为受主单元 ,用有效扩散系数来确定离 散方程的系数. 这一方法的具体实施过程见文献 [4 ].
2 计算方法 2.........