求解温盐双扩散系统的一种高精度方法

2 U k a Z+ V k Z Z+ e k Z+ f k a) + Pr Rt yZ ( T - RdS ) ( T - RdS ) - ya /J ( 2) Z a 1 J j T - j T = 1 (T Ta a2 Za a Z J ( 3)
T + t
2 U F a Z+ V TZ Z+ e TZ+ f Ta) S + t 1 J j S- j S = 1 (T Sa a2 Z a aZ Le J
1, j , k - 8
两个条件没有本质上的区别 ,只是流态或温盐场的 图形倒过来了 。 其余边界的边界条件与文 [ 2] 的相 同 。 式中 , A = 5为腔体的长宽比。 计算中 , 选取 的网格节点为 81× 51, Rt = 8 × 105 , Pr = 1. 5,
- 5 Le = 1, 时间步长为 Δt = 10 , 由于温度 Raylei gh
2 2 2 2 2 2
354
计
算
力
学
学
报
第 19 卷
度 控 制 函数 , j为 流 函 数 ,k为 涡 量 , T 为 温 度 , S 为盐度 , Pr 为 Prandt l数 , Rt 为温度 Rayleig h 数 , Rd为浮力比 , Le 为 Lewis 数。边界条件为 在固壁边界上 u = v = j= 0, ( 5) T 在绝热边界上 = 0, ( 6) n S 在非穿透边界上 n = 0, ( 7) 如 果将顶部边界看成是海面 , 则温度 、流函数 和涡量可给定 [ 3 ]: j= k= 0, T = Ts ( x ) 对于盐度 ,可给定盐度分布或盐度通量 ,即 S = Ss ( x ) , ( 9a ) S = W Qs ( x ) ( 9b) z 其中 ,参数 W 表示盐度通量的强度 , Qs ( x ) 盐度通量 的空间分布结构。 ( 8)
在 时 间 离 散 方 面 , 采 用 下 列 的 Adams
第 3期
詹杰民 , 等 : 求解温盐双扩散系统的一种高精度方法 x 1 A 2 6 π sin p pπ 2
- 4
355
Qs ( x ) = 3cos pπ
( 16)
S′ = - 5× 10 | sin( 2 π ( 1000t - x ) )| 0. 07≤ t < 0. 08
一正一负先后输入扰动能使扰动总量保持平衡 。 在 盐度通量条件下 , 当对称状态受到扰动时 , 质点将获得加速 。 这种连续的加速可能会引起不稳 定性 , 当 W 较小时 , 对称状态可能不会被 破坏 (图 2)。 随着 W 增加 , 不对称结构和反转结构将会出现 , 甚 至变成反转 对称结构 (图 3-7)。数值结 果与文 [ 13 ] 和文 [ 2] 的非常一致 。
Z+ V Z+ e 2 U Sa SZ SZ+ f Sa)
( 4)
而 T= xZ+ yZ, U= x axZ+ yayZ, V= xa+ ya, J = xayZ - yaxZ, e= J Q 和 f= J P 式中 t 为时间 , a ,Z 为边界拟合坐标 , P , Q 为网格密
收稿日期 : 2001-01-10; 修改稿收到日期 : 2001-08-10. 基金项目 : 国家自然科学基金 ( 19802023)及国家教育部留学 回国人员科研启动基金 。 作者简介 : 詹杰民 , ( 1963-) ,男 ,教授 ,博士生导师 .
界层的流动问题时 , 有限元法较有限差分法更灵 活 , 但当引入边界拟合坐标系后 , 有限差分法除了 保持其编程简单的优点外 , 还具有有限元法灵活布 置网格的优点。
1 詹杰民 1 , 李毓湘 2 , 郑王 君婷 ( 1. 中山大学 应用力学与工程系 , 广州 510275; 2. 香港理工大学 土木及结构工程系 )
摘 要 : 引入边界拟 合坐标系来研究温 盐双扩散系统 。 为了提高 求解的精确性 ,对流 项采用四阶 精度的迎风 格 式 , 扩散项和涡量方程的浮力项则采用四阶精度的中心差 , 因此本文的方法是高精度的方法 。首先针对温度占优 ( Rd= 0. 32 ) 和 盐度占优 ( Rd = 1. 68 ) 的情形进行了验证性计算 ,得到了与前人一 致的结果 。 进一步 , 本文系统给 出了不同的盐度通量强度下的流动形态 ,包括对称结构 , 不对称结构 ,反转结构等 , 结果与前人的吻合 。 关键词 : 温盐双扩散 ; 高精度格式 ; 多平衡态 中图分类号 : O315. 2 文献标识码 : A
参考文献 ( References ):
[1 ] L ee J W , Hyun J M . Doublediffusiv e conv ectio n in a rectang le with opposing ho rizo ntal tempe ratur e gr adient and co ncentr atio n gr adients [ J ]. Int . J . : 1619-1632. Heat Mass Transf er , 1990, 33 [2 ] Q uo n C, Ghil M. so luta l conv ection M ultiple equilibria in thermo due to saltflux bounda ry
2 控制方程和边界条件
本 文考虑的是温度 、 盐度并存的双扩散系统。 控制方程为结合了 Boussi nesq 假设的二维非定常 Navi er-Stokes 方程。在边界拟合坐标系下 [ 9-11 ] , 无 量纲形式的流函数、涡量、 温度和盐度方程为 1 j a a- 2 U j a Z+ V j Z Z+ e j Z+ f j a) = k ( 1) 2 (T J k + t Pr 1 j k jk = 2 (T k a aJ Za a Z J
第 19 卷第 3期 2 00 2年 8月
计 算 力 学 学 报 Chinese Journal of Computational Mechanics
V ol . 19 , N o. 3 Aug ust 2002
文章编号: 1007-4708( 2002) 03-0353-06
求解温盐双扩散系统的一种高精度方法
Bashf or th 方法: f n+ 1 = f n + Δ t 3 2 f t
n
-
1 2
f t
n- 1
( 14)
4 算例与讨论
4. 1 方法的验证 首先针对温度占优 ( Rd = 0. 32) 和盐度占优 ( Rd = 1. 68) 的情形进行了验证性计算 , 为了便于 比较 , 物理方程采用了文 [ 2] 的形式 , 由于采用不 同的特征尺度来定义无量纲量 , 文 [ 2] 采用的方程 与 ( 1) -( 4) 对应的物理平面下的形式略有不同。边 界拟合坐标系如图 1 所示 , 局部加密网格有两个目 的 , 一是适应物理量的剧烈变化 , 二是在边界附近 差分格式的精度会有所下降 ,加密网格有助于提高 精度 , 弥补这一缺陷 。为了便于将模型比拟为两端 为南北极 , 中央为赤道的海洋模型 , 本文不是采用 文 [ 2] 底部形式的边界条件 , 而是采用文 [ 13 ] 等海 面形式的边界条件 ,即 : Ss ( x ) = Ts ( x ) = x 1 1 cos 2 π A 2 2 + 1 ( 15)
2, j
u1 k = ( u1 )i , j [ - 6(k)i+ a i, j 60 (k)i+
1, j
+
+ 40(k ) i , j - 120(k) i - 1, j +
30 (k)i - 2, j - 4(k) i - 3, j ] / ( 120 Δa ) ( 12) 如果 u 1 < 0, 则 u1 k = ( u 1 ) i , j [ 4 (k)i+ a i, j
3
8
f
a
2=
16 ( f i+ 1, j , k + f i - 1, j , k ) - 30f i , j , k - ( f i+ 2, j , k + f i - 2, j , k ) 12 Δa
2
( 11)
wk.baidu.com
而扩散项则采用四阶精度的迎风差 如 , 对于
[ 14 ]
来近似。例
1 j 1 jk , 如果 u1 = > 0, 则 J Z J Za
- fi 2 Δa
1, j , k
-
1 3
f i+
2, j , k
- f i4 Δa
2, j , k
( 10)
2
数相当于文 [ 2] 中的 Ra = A Rt = 10 , 因此 , 这里 的条件与文 [ 2]的一致 。 计算结果如图 2, 3所示 , 与 文 [ 2, 14 ] 的一致 , 说明方法能正确反映温盐双扩 散对流的运动情况 , 由于温度和盐度在 Le = 1下方 程和边界条件完全相同 ,因此温度场和盐度场也完 全相同 。 为了定量验证方法的准确性 , 我们曾采取 与文 [ 2] 完全相同的方程形式与边界条件针对盐 度占优 ( Rd= 1. 68) 的情形进行数值实验 , 结果是: 两者涡心位置流线重合 , 其它位置流线也相当吻 合。 4. 2 温盐驱动流中的对称破坏和翻转形态 由于该算例与上例几何形状相同 ,因此仍采用 了图 1 所示边界拟合坐标系。温度边界条件仍取 ( 15) 式 , 盐度边界条件为通量边界条件 , 能够将南
1, j 3, j
- 30 (k)i+
2, j
+
120 (k)i+
- 40 (k)i , j - 60(k)i - 1, j + ( 13)
6 (k)i - 2, j ] / ( 120 Δa )
北极冰川融化 , 赤道海水蒸发的现象反映出来。具 体形式采用文 [ 13 ] 的形式:
S′ = 5× 10 |sin( 2π( 1000t - x ) )| 0. 03≤ t < 0. 04和
- 4
356
计
算
力
学
学
报
第 19 卷
第 3期
詹杰民 , 等 : 求解温盐双扩散系统的一种高精度方法
357
5 结 论
本文提出的边界拟合坐标系下求解双扩散问 题的方法是一种高精度的方法。由于对流项采用的 是 4阶精度的迎风格式 , 该方法具有较好的数值稳 定性。 同时 , 计算表明 , 方法能够将对称结构 ,不对 称结构 , 反转结构等复杂的多平衡状态准确地模拟 出来。 因此 ,适于直接模拟双扩散问题 , 并能适应复 杂的边界条件和几何形状。
3 高精度差分格式
在 本文 , 方程 ( 1) 用超松驰方法求解 , 计算精 度 控 制 为 前 后 两 次 迭 代值 之 差 的 绝 对 值 小 于 10 。 在 空间离散方面 , 除了非线性的对流项之外 , 方程 ( 1) -( 4) 中的一阶和二阶 导数采用四阶精度 的中心差:
= a f 4 3 f i+
其中 p= 2. 6, 其余边界的边界条件与文 [ 3] 的相同。 计 算 中 Rt = Rs = 4 × 104 , P r = 2. 25, Rd = 1, Le = 1, Δt = 10- 5。要得到对称破坏的流 动形态 , 需要引入扰动。本文在距腔体顶部的相邻 坐标线上的左半部分 , 引入如下形式的扰动 :
1 引 言
由 于诸如干燥过程、 晶体生长技术 、 室内空气 调节 , 核反应堆冷却 , 太阳能吸收器和温盐海洋环 流等的自然界和工业流动系统的输运过程中存在 着温盐对流现象 , 因此 , 这样的双对流扩散系统正 越来越受到重视。 目前 ,国内外已有大量关于双对流扩散系统问 题的研究 , 文 [ 1 ] 对在一个封闭腔体内的双扩散对 流系统进行研究 , 发现了当温盐浮力比 λ 是中等大 小时 , 流场内出现了多层的流动结构。文 [ 2 ] 对温 盐溶液的二维对流系统进行研究 , 指出了瑞利数和 盐度通量因子是不稳定的关键参数 ; 文 [ 3 ] 还对具 有较小高度 -长度比的二维腔体中的温盐溶液 , 研 究了单调或非单调温度 、盐度边界条件对流体运动 的影响 , 发现当盐通量强度足够大时 , 对两种边界 条件系统都会发生 2倍分叉行为 —— Hopf分叉。 求 解双对流扩散系统的方法有很多 , 例如 , 文 [ 4 ] 用 Galerki n 方法针对垂直和倾斜的腔体研究 了非线性现象 ; 文 [ 5 ] 采用有限元法进行了研究 , 并 与文 [ 6 ] 等的实验结果进行了比较 ; 文 [ 7 ] 用差 分 法 模 拟 了 Chen 等 的 实 验 结 果 ; 文 [8 ] 则 用 Cheby shev 配置法研究了温盐分层系统中层与层 之间的合并现象。在求解边界形状复杂 ,并带有边