曲线坐标下三维水动力学生态综合模型的建立.pdf

水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 1 文章编号 0559 9350 2005 08 0891 09 曲线坐标下三维水动力学生态综合模型的建立 倪浩清 李福田 中国水利水电科学研究院 北京 100044 摘要摘要 为了获得高分辨率 高精度及快速的数值模拟计算成果 以便解决复杂边界的工程实际问题 本文研究了通用 形式表示的湍流代数应力模型流动基本方程组在拟合坐标系中的转换 通过该模型方程组的转换 并结合具有生物化 学作用过程的生态子模型 建立了曲线坐标下三维水动力学的生态综合模型 便于复杂边界工程问题的方程组求解及 对未来情况的预测和评价 关键词关键词 湍流基本方组转换 贴体曲线坐标 水动力生态综合模型 中图分类号中图分类号 X171 文献标识码文献标识码 A 本文中所提出的水动力学生态综合模型实质上是污染物质在水环境中因物理 化学 生物作用及影响因素 之间相互关系的数学描述 这方面比较有代表性的模型有Grenney 1 在QUAL I综合模型基础上发展的QUAL 综 合水质模型 倪浩清等 2 提出的 三维湍流代数应力 热流模型 简化的统一二阶矩湍流模型 沈永明 3 发表 的 近海环境准三维k 湍流流动和水质的精细模拟 雒文生等 4 发表的 垂向二维湍流与水温水质耦合模型 李清雪 5 提出的 海湾浮游生物及氮营养盐生态水动力学模拟 等 这些模型对水环境数学模型发展均取得积 极的推进作用 但在这些模型中没有综合考虑海湾的弯曲边界 各向异性和生态动力作用对湍流数值模拟计算 的影响 因此 考虑到湍流 温度 生态及流动边界相互耦合影响的关系而建立一个完整的曲线坐标下三维水 动力学生态综合模型是十分必要的 1 控制方程在拟合坐标系中的守恒方程 1 1 模型的选择与拟合坐标系1 1 模型的选择与拟合坐标系 由于江河湖海的岸边界弯弯曲曲极不规则 其中的流体质点受到离心力及 各向异性湍流雷诺应力场的作用而导致产生二次流 它与弯道轴线方向的主流合成各向异性的螺旋流 这种二 次流能导致主流流场 温度场及污染物浓度场再分布 采用湍流代数应力 热流模型 则能模拟出离心力和各 向异性湍流应力场所引起的二次流 因此 选用具有考虑生化反应影响的各向异性湍流代数应力 热流模型作 为水环境数学模型是符合天然水域生态水动力学特性的 同时 针对天然河流 海湾边界弯曲不规则性 采用 拟合坐标系变换办法 则是解决数值模拟计算问题中重要途径之一 1 2 变换坐标系中的守恒方程式1 2 变换坐标系中的守恒方程式 时均化的N S方程及质量输运方程可以写成如下的笛卡儿坐标向量微分 方程的守恒形式 Q z H y G x F t M 1 收稿日期 2004 04 02 作者简介 倪浩清 1931 男 江苏无锡人 教授级高级工程师 主要从事水利水电工程 海洋工程和火 核电冷却水工 程中的水动力学 波动力学及水环境动力学研究 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 2 z Dw P w wv wu w H y Dv vw P v vu v G x Du uw uv P u u F S f f f O Q w v u P M zz zz zy zx yy yz yy yx xx xz xy xx z y x 2 2 2 2 其中 jiijij v v 3 这里 j i i j ij x u x u 为层流剪切应力 jiv v 为雷诺应力 假设边界拟合坐标变换将物理域 笛卡儿坐标系 x y z t 变换到计算域 非正交曲线坐标系 即 tzyxtzyxtzyxt 4 同时将式 1 的笛卡儿坐标向量微分方程的守恒形式转换到曲线坐标系中微分方程的守恒形式中去 文献 6 通过上述边界拟合坐标变换及其方程转换 便得到了Lagrange型曲线坐标系中的微分方程的守恒形式 本 文作者将Lagrange型中的函数行列式D的倒数转化为Eulerian型中Jacobin行列式并考虑了Lagrange型中变换 坐标函数的导数 即 ij i j a Jx 1 5 式中 Jacobin行列式为 zyx zyx zyx J 6 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 3 aij为下列伴随矩阵A 的元素 度量系数 yxyxyxyxyxyx zxzxzxzxzxzx zyzyzyzyzyzy A 7 这样便可得到下列Eulerian型曲线坐标系中的转换控制方程守恒形式 QHGFM 8 其中 332313322212312111 HaaGFaHHaGaFaGHaGaFaFJQQJMM 9 2 代数应力控制方程转换及其离散 2 1 基本方程2 1 基本方程 在天然水域中 由于温度 盐度或含沙量的不同 经常出现水体分层现象 倪浩清 周力行 等于1987年4月首次提出了改进的双方程k 模型 7 首次运用于分层流的模拟中 并取得了一定的成功 但 由于改进的k 模型属于经验性关系 求解缺乏系统性 稍后倪浩清 周力行等于1987年7月提出了改进的代 数应力模型 8 取得了满意的结果 这就充分发挥了代数应力湍流模型固有的各向异性的特征 从而克服了双 方程中k 模型经验性的弱点 从此 对分层流的模拟取得了较大进展 但对于不规则边界的天然河流及海湾 等问题 其使用受到了很大限制 为此 采用坐标变换是解决该问题的途径之一 现将笛卡儿坐标的三维湍流 代数应力 热流模型的通用形式列表1 以便该模型方程组的转换 表1 表1 S xx SS xx v xt iiii j j 10 S 方程 S S 连续 1 0 0 0 动量 vi e Tg x v xx P i i j e ji ji ji j j i j vv xx v x vk c x 2 能量 T T e H j jjTj v xx Tk c x 2 浓度 C C e z w SQ C sCC cv xx Ck c x j jjcj 2 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 4 紊动能 k k e G B BPGP x kk c xx k vv k c x bk jjl ijk j 2 紊动能 耗散率 e k cB k cG k c l 2 23 B k cP k cG k cP k c x kk c xx vv k c x bK jjl ij j 331 1 2 1 3 2 11 132 32 1 cPcPcPcPc c k vv bkijbijijji m i mTiTTi T i x v vcgcP c k v 23 1 1 k k x T v k c 1 1 m i mcicCi c j x v cvcccgcP c k cv 23 1 1 应力 热流 浓度 代数式 k k c x C cv k c cc 1 1 表中 SSS 为紊动应力代数模型总的源项 其中 S 为k 模型原有的湍流黏性 各向同性 的源项 k S 为考虑紊流应力 热流各向异性影响的附加源项 H为热交换项 包括太阳短波辐射 蒸发和热传导 c为 不同生态变量 状态变量 浓度 Q C为生化反应项 亦即是状态变量之间生物化学相互作用项 又称内源 SC是 源汇项 外源的输入及输出 如内陆河道径流的输入 潮流量向内陆输入 C s w 为悬浮碎屑和浮游植物的颗粒 物质沉降速度 通用常数为 c1 2 2 c2 0 55 c3 0 55 ck 0 24 c 0 15 c 1 1 44 c 2 1 92 c 3 0 8 5 0 5 0 2 3 8 0 3 1 0 1 0 1 5 0 5 0 2 3 8 0 33113211 iib iT e i j i i j j i ete ccccTkTTT vgP x T gB x v x v x v G ccccCccc ijjibij l i jl l j liij j i jik vgvgP x v vv x v vvP x v vvP 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 5 k i k k kiCi k i k k kiTi x v vc x C vvP x v v x T vvP 2 2 方程转换2 2 方程转换 主要根据变换坐标系中的控制方程的守恒形式 8 式 9 进行转换 1 对流项 i i v x 转换 将F u G v H w 代入控制方程的守恒式 可得 332313322212312111 awavauawavauawavau iji i jij i i i i i i i avavawavau 321 令此式中 iiij Uva 即 i j iij i i i U J va J v x 11 11 如果F G H 类似上法推导得 j ij i ij ji a Jx a Jx 11 或 12 2 扩散项 ii xx 的转换 将 j j a Jx F 1 j j j j a Jz Ha Jy G 32 代入控制方程的守恒式 可得 vi 直角坐标系中的速度分量 Ui 曲线坐标系中逆变速度分量 如果在计算空间中求解时 可采用选择逆变速度分量作为 求解变量 9 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 6 j njni i j ij j ij j ij ii j j j j j j j j j j j i j j j j j aa J aa J aa J aa J aa J aa J aa J aa J aa J aa J aa J aa J aa J 332211333232131 323222121313212111 即 j njni iii aa JJxx 1 13 3 动量源 y w zy v yy u x S x w zx v yx u x S vu z w zz v yz u x S w 转换 首先转换S u 将 j j u a Jx u F 1 j j j j w a Jx w H v a Jx v G 11 代入控制方程的守恒式 可得 n j nji i j n nij ij j j j j j j j j j j j j j j v j j j v aa J v aa J w aa J v aa J u aa J w aa J v aa J u aa J w aa J aa J u aa J 1 1331231131 321221121 311211111 14 S v Sw的转换 同上法 可得 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 7 n j nji i v aa J 2 15 n j nji i v aa J 3 16 所以动量源项转换式为 n j mnji ij i i v aa JJx v x 1 17 4 产生项 j i i j j i e x v x v x v G 转换 应用转换公式 j ij i a Jx 1 得 222 2 2 1 in n j jn n i in n j jn n ie j i i j j e a v a v a v a v Jx v x v x v G 18 5 压力项 i x p 转换 应用转换式 12 可得 pa Jx p ji ji 1 19 6 浮力项 y T gB i T e 转换 应用转换公式 12 可得 Ta J gB ij j i T e 1 20 应用拟合坐标转换式 11 20 转换后取得了代数应力湍流模型拟合坐标的基本方程组的通用形式 列于 表2中 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 8 表 2 表 2 JSaa J SJSJaa J U t J l nlnl ll nlnl l i i 21 S 方程 S S 连续 1 0 0 0 动量 Ui e TgS v aa J pa i n j mnjl e l ij i i jijk k n j mnjl t ll k nlnl t l vva v aa J v aa J 能量 T T e i SH jjk ki nini t i va T aa J 浓度 C C e i SCwaSQ C j j CC 3 cva C aa J jjk ki nini t i 紊动能 k k e i S JBJG JBJPJGJP k aa J ka J vva k c bk i nini k i lm m ijjik i 1 紊动能 耗散率 e i S k c JB k cJG k c 2 2 31 JB k cJP k c JG k cJP k caa J a J vva k c b k i nini i lm m ijji i 33 1 1 1 1 3 2 11 132 32 1 cPcPcPcPc c k vv bkijbijijji imn n i mTiTTi T i va v J vcgcP c k v 1 1 33 1 应力 热流 浓度 代数式 Ta J v k c J kn n k 11 1 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 9 imn n mciccij va J cvcccgcP c kJ cv 1 1 23 1 Ca C J cv k c ccJ kn n k 11 1 表中 H Q C SC c s w 的定义同表1下注明 Ui即U V W是计算空间域上 的速度分量 称为逆变速度 SSSvaU jjni 9 为紊动应力代数模型总的源项 其中 S 是k 模型原有的紊流黏性 各向同性 的 源项 S 为考虑紊流应力 热流及浓度各向异性影响的附加源项和k 模型各向同性的扩散项 2 321 3 33 32 22 21 11 1 JJJ S JJJ aa J aa J aa J aa J i nnnnnn i nini i 其中 1 2 3为相应坐标轴 与 与 与 交角大小的3个参数 aij为Jacobian矩阵J的伴 随矩阵A的元素 度量系数 三维Jacobian矩阵J及其伴随矩阵A见式 6 7 2 kckc k tte 11 11 11 1 1 2 1 222 2 ikj j knnkj j kiCi ikj j kkj j kiTi ijjibijilk k lijlk k liij ilj j likiibij jT e i in n j jn n i in n j jn n ie va J vcCa J vvP va J vTa J vvP vgvgPva J vvva J vvP va J vvPvgPTa J gB a v a v a v a v J G n 2 3 边界条件转换2 3 边界条件转换 物理域上的控制方程转换到计算域上后 边界条件亦应作相应的转换 物理空间 域上的第一类边界条件 变换后其值保持不变 但含有导数项的边界条件 则应作变换处理 梯度变换表达式 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 10 llj j l m lm ka J ka J 11 22 标量函数的法向量变换表达式 2 3 2 2 2 1 ln 2 3 2 2 2 1 332211 iii l mmm mmm i i aaa ka aaa kakaka n i 23 综合式 22 23 得第二类边界条件的标量函数的法向导数变换表达式为 j mjmk inni aa aaJ n n i i 1 24 式中 i为自由指标 2 4 方程离散2 4 方程离散 任意曲线坐标系下的代数应力通用输运方程式 21 展成 JS JJ a J SJSJ JJ a J WVU t J 25 式中 为通用变量 如u v w k 和T等 分别为 方向的通用扩散系数 S 为源项 不同的基本方程式 对应 取不同的值或表达式 在计算空间域中的交错 网格上来离散通用控制方程式 25 将压强p 湍动能k 耗散率 温度T或浓度C等变量布置在主网格的节点 上 将速度分量上U V W放在主网格中间的交界面 在时间间隔 t t t 内任意曲线坐标系中 应用有限体 积法将式 25 对计算域中的控制容积P积分有 tt t e w n s h l tt t e w n s h l tt t e w n s h l tt t e w n s h l dtdddS J dtddd JJJ dtdddWVUdtddd tJ 1 1 26 经整理得如下的三维迎风差分格式离散方程式 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 11 bAAAAAAA LLHHSSNNWWEEPP 27 运用SIMPLEC算法即可对此方程进行求解 上式中 P i i P PLHSNWE P P S J AAS J AAAAAAAA 11 28 PP P P P AS J b tJ A 1 1 29 系数AE等取值取决于所采用格式 如当界面上的函数及其导数均采用分段线性的线型时 采用迎风格式 对流 项在界面上取的分布函数为 0 0 0 0 0 0 0 0 0 0 0 0 W W W W V V V V U U U U H P h P L l P S s N P n P W w E P e 30 系数有 llLhhHssS nnNwwWeeE FDAFDAFDA FDAFDAFDA 0 0 0 0 0 0 31 这里符号 表示取各量之中的最大值 F与D是界面上的流量与扩散系数 其计算式为 llhhss nnwwee WFWFVF VFUFUF 32 l l h h s s n n w W e e J D J D J D J D J D J D 33 式 29 33 中各量定义 为界面间距离 为节点间距离 3 生态子模型的引入 富营养是由于水体内营养成分的增加 导致生物 特别是浮游生物即藻类 的生产能力异常增加的过程 从 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 12 理论上讲 富营养化在任一水生生态系统内都有可能发生 但实际上主要出现在湖泊 水库 河口 海湾等较 封闭的水域 富营养化最显著的特征是水面藻类 主要是蓝藻 绿藻 异常增殖 成片成团地覆盖在水面上 其 中出现在湖面上的称为 水华 或 湖绽 出现在海湾水面上称为赤潮 水质富营养化被列入水污染范畴还 是近二三十年的事情 但由于发展快 危害大 处理难 恢复慢 富营养化已被列为全球性的重点水污染问题 加以研究 引起富营养化的物质 主要是浮游生物生长所必须的碳 磷 氮 硫 镁 钾等20余种元素 其中特别是 碳 磷 氮更为重要 文献 5 通过对海湾的调查研究 提出了浮游生物生态系统子模型 本文作者将该子模 型引入前面所建立的湍流模型 即为曲线坐标下的三维水动力学生态综合模型 3 1 海湾浮游生物生态系与氮营养盐生态子模型描述3 1 海湾浮游生物生态系与氮营养盐生态子模型描述 该模型将海湾生态系分为浮游植物 浮游动物 营养盐 悬浮碎屑四个功能群 包括NO3 N1 NH 24 N 浮游植物 P 浮游动物 Z 悬浮碎屑 D 五个状态 变量 该生态系统模型框图如图1所示 其中浮游植物代表初级生产者 浮游动物代表二次生产者 氮营养盐 考虑了氨氮和硝酸盐氮 浮游植物吸收氨氮和硝酸盐氮后增长 二次生长者摄食初级生产者后增殖 浮游植物 和浮游动物死亡后成为悬浮碎屑 碎屑通过水解转化为溶解性有机氮然后进一步矿化为氨氮 但在模型中没有 考虑溶解性有机氮这个状态变量 而是认为碎屑直接矿化为氨氮 悬浮碎屑和浮游植物沉降到底层形成底栖碎屑 底层 不考虑硝化作用只考虑反硝化作用 也可通过再悬浮 变成悬浮碎屑 水层的化学过程只考虑了硝化过程 并直接由氨氮转化为硝酸盐氮 没有考虑亚硝酸盐这 一中间变量 模型中没有考虑细菌这一状态变量 细 菌在氮循环中的作用隐含在营养盐再生和碎屑的矿 化参数中 模型中状态变量单位为mmol Nm 3 3 2 生态子模型中生物化学相互作用项的数学 描述 3 2 生态子模型中生物化学相互作用项的数学 描述 1 浮游植物方程中生物的相互作用项 考虑光 营养盐 10 及温度影响的浮游植物比增长速率 浮游动 物的摄食速率 浮游植物的死亡率及浮游植物的呼吸 率这三种浮游植物损失因素 浮游植物对硝酸盐氮和 P为浮游植物 Z为浮游动物 N1为硝酸盐氮 N2为氨氮 D为悬浮碎屑 图1 海湾生态系统模型中的氮循环 氨氮总吸收 从而建立了浮游植物方程中生物的相互作用项表示为 P NK N e NK N e I I eQ a N n I I sat zTa p P sat z 2 2 1 1 1 max1 21 0 1 34 式中 1为浮游植物呼吸和排泄比 pmax为浮游植物0 时的最大比增长速率 a1为温度影响系数 a1 0 061 T为温度 Isat为饱和光强 Isat 112000Lux 33 4W m 2 I z为透光层的平均光强 光线沿水深的变化可以根据比尔 定理确定 即Iz 0 5I0 H e kzdz 其中k为衰减系数 它由三部分组成 k k w kc kd kw m 1 代表纯水引起的 衰减 kc m 1 是考虑浮游植物自遮由叶绿素cha1 m 1 引起的衰减 k c schal s m 1mg chal m3 为浮游植物自 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 13 遮系数 kd m 1 是由水体中无机悬浮物质引起的光的衰减 I 0为海面入射光的强度 海面入射光强随时日 地 理纬度 天空状况及海面反射率而变化 根据Dobon et al 11 分析结果 由云量估计太阳辐射通量公式 I0 t 1 R Q0S Ai BiS 其中 Q0 1368Wm 2为太阳辐射常数 R为海面反射率 A i和Bi是有关云量的常数 根 据Brock 12 的公式 S有如下表达式 S sin sin cos cos cos 0 2618t 0 40927 sin 2 284 n 365 其中 是地理纬度 是入射倾角 t是小时变化 n为一年中日数 N1为硝酸盐氮 N2为氨氮 Kn和Ka为硝酸 盐氮和氨氮吸收的半饱和常数 它反映了浮游植物硝酸盐和氨氮吸收的容易程度 为经验常数 它表示氨氮 存在对浮游植物硝酸盐氮吸收的抑制 z为最大摄食速率 为Ivlev摄食参数 一般取 0 33 2 0 mmol N 1 z为浮游动物浓度 P为浮游植物浓度 P 为浮游动物停止摄食的下限阈值 pm为浮游植物自 然死亡率 2 浮游动物方程中的生物相互作用项 考虑了浮游动物对浮游植物和悬浮碎屑的摄食同化速率 浮游动 物排泄速度 浮游动物死亡率及高级动物对浮游动物的摄食速率从而建立了浮游动物方程中的生物相互作用项 可用下式表示 Q z z DZ 1 e D z PZ 1 e p p zeZ zmZ zgZ 35 式中 为浮游动物的动化系数 D为碎屑浓度 ze为浮游动物的颗率粪便排泄速度 zm为浮游动物的死亡 率 zg为高级动物对浮游动物的摄食速率 3 悬浮碎屑方程中的生物相互作用项 考虑了浮游植物的比死亡速率 浮游动物的比死亡速率 浮游动 物对浮游植物净摄食同化速率 浮游动物对碎屑摄食同化速率及碎屑砂化速率从而建立了悬浮碎屑方程中的生 物相互作用项可用下式表示 Q D pmP zmZ 1 z 1 e p p zDZ 1 e D D 36 式中 z为浮游动物的最大比增长速率 为碎屑的矿化速率 4 氨氮方程中的生物化学相互作用项 氨氮包括浮游植物的吸收 悬浮碎屑在微生物作用下的矿化 浮 游植物与动物的呼吸与排泄 还有在硝化细菌作用下将氨氮氧化为硝酸盐氮硝化作用 在本模型中没有考虑反 硝化作用 硝化过程的速率与环境中溶解氧的浓度和温度有关 在此没有考虑溶解氧的影响 因此 氨氮方程 中的生物化学相互作用项可表示为 DZPNkP NK N e I I eQ ep T an a I I sat zTa p N sat z 12 20 2 2 1 max 2 20 0 1 37 式中 kan为水温在20 时硝酸盐氮的比硝化速度 为温度修正系数 p浮游植物的比增长速率 1为浮游 植物的呼吸与排泄比率 ze为浮游动物比排泄速率 5 硝酸盐氮方程中的生物化学相互作用项 硝酸盐氮方程中的生物化学相互作用项只考虑了氨氮的消化 和浮游植物的吸收 具体表达式如下 水水 利利 学学 报报 2005 年 8 月 SHUILI XUEBAO 第 36 卷 第 8 期 14 2 20 1 1 1 max 1 20 0 21 NkPe NK N e I I eQ T an N n I I sat zTa p N sat z 38 5 结语 综观本文所述 对湍流模型的选择 转换及水动力学的生态综合模型的建立作以下几点结论 1 选择各 向异性的湍流代数应力 热流模型是较为合适的 它能预报旋流 分层流及环流等效果较好 2 关于不规则复 杂边界问题 应用拟合坐标变换办法对控制方程进行转换 这是非常必要的 3 根据复杂的流动和复杂的边 界条件 确定了转换后的湍流代数应力 热流模型并在其浓度的对流扩散方程中加上考虑生化过程的源项 这 样就形成了较为完整的三维水动力学生态综合模型 此种模型可同时模拟几十个状态变量及其相互作用 并统 一考虑了各种物理 化学和生物过程的综合作用 其功能齐全 理论先进 能进行实时预报 能起到环保作用 4 代数应力模型湍流方程的简化 源项的改进及方程的求解和实例分析 将另文论述 参参 考考 文文 献 献 1 Greney QUAL 综合水质模型 A 环境水力学进展 C 武汉 武汉水利电力大学出版社 1998 2 倪浩清 等 紊流数学模型在环境预测中研究及其展望 J 水利学报 1986 5 24 38 3 Shen Yongming Three dimensional modeling of hydrodynamics and water quality in coastal waters R ALCS Report HR Walling ford 1994 4 雒文生 河流多维综合水质模型 垂向二维湍流与水温水质耦合模型 环境水力学进展 C 武汉 武汉水利电力大学 出版社 1998 505 512 5 李清雪 海湾浮游生物及氮营养盐生态水动力学模型 R 天津 天津大学力学与工程测试系 2000 6 吴江航 韩庆书 计算流体力学的理论 方法及应用 M 北京 科学出版社 1988 7 倪浩清 周力行 等 明渠温差异重流中湍浮力回流的数学模型 J 中国科学 A辑 1987 4 8 Ni Haoqing Zhou Lixing et al Numerical simulation of Nonisotropic Turbulent Buoyant Recirculating Flows By an Algebraic Stress Model A Proceedings of th International Cmference on Fluid Mechanics C Beijing 1987 9