求解自然对流换热问题的高效方法
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1.2 雅格比坐标变换与直接求解法的比较 对于紧致有限差分格式应用到非结构化 网格上求解 Navier-Stokes 方程主要有两种方 法：一种是雅克比坐标变换法（JT）。 这种方法 被用在大多数的稳态流动中， 但是不经常用 在非稳态流动的直接数值模拟中。 另一种方 法是直接求解法 （FIM）， 就是直接在非等距 网格上进行数值计算。 ① 雅克比变换法：只需要构造在等距网 格上的有限差分方法， 通过物理空间到计算 空间的变换将其应用到非等距网格中， 一阶 导数和二阶导数的变换过程为：
［2］Y. MA，D. FU， Numerical solution of the incompressible Navier -Stokes equations with an upwind compact difference scheme ［J］， Numer Meth Fiulds， 30:509-521， 1992.
+
坠2ν 坠y2
）-（u
坠ν 坠x
+ν
坠ν 坠y
）
+fy
（2.6b）
为了保证边界处的压力满足连续性方
程，设
ω= 坠ν - 坠u 坠x 坠y
（2.7）
由 连 续 性 方 程 （2.2）、 （2.6a） 和 （2.6b） 可 得：
坠p 坠x
=-
坠u 坠t
-
1 Re
坠ω 坠y
-（u
坠u 坠x
+ν
坠u 坠y
当求解压力时， 压力边界条件采用的是
Neumann 边界条件，由（2.1）分 别 展 开 为 x 方
向和 y 方向上的动量方程得：
坠p 坠x
=- 坠u 坠t
+
1 Re
（
坠2u 坠x2
+
坠2u 坠y2
）-（u
坠u 坠x
+ν 坠u 坠y
）
+fx
（2.6a）
坠p 坠x
=-
坠ν 坠t
+
1 Re
（
坠2ν 坠x2
β=-
（h2i+1+hi+1hi-h2i）hi （h3i+h3i+1+4hi+1h2i+4hi h2i+1）
方程（1.4）中系数：
A1=
6 h22-h32
，
B1=-
6 （h2-h3）h3
C1=-
6 h2h3（h22-h32）
，
β1=
（2h2+h3） （h2-h3）
参考文献
［1］S. K. Lele， Compact finite difference schemes with spectral -like resolution ［J］， Journal of Computational Physics， 103: 16-42， 1992.
式中的压力梯度项 p 可以通过求解一阶偏 导数的四阶紧致差分格式求得，（2.3） 式中的 扩散项 2ν軋 可以通过求解二阶偏导数的三阶 紧致差分格式求得， 压力泊松方程采用二阶 中心差分格式离散求解。 由（2.3）式得到先验
速度軋ν* 不满足连续性方程，在（2.5）式中被修
正。
2.2 边界条件的处理
2．求解 N-S 方程的高精度算法 2.1 投影法
无量纲化后的二维纳维斯托克斯方程的 守恒形式和连续性方程为：
坠ν軋+（ν軋· ）ν軋+ p= 1 2ν軋+ f軆
坠t
Re
（2.1）
·ν軋=0
（2.2）
其 中 ν軋=（u，ν）为 速 度 矢 量 ，p 为 压 力 ， f軆= （fx，fy）为源项，Re 为雷诺数。
点，并比较了在三种不同的非等距网格下，直
接数值求解法与雅克比变换法的区别与联
系，为方法的选择提供了理论依据。 将紧致差
分格式与投影法相结合，构造了求解 Navier-
Stokes 方程的新方法。 通过数值实验表明该
方法比经典有限差分具有更高的精度， 是一
种求解 Navier-Stokes 方程的有效方法。
附录
方程（1.1）中系数：
Ai=
12hi+1 （h3i+h3i+1+4hi+1h2i+4hih2i+1）
Bi=
12 （h2i+1+3hi+1hi+h2i）
Ci=百度文库
12hi+1 （h3i+h3i+1+4hi+1h2i+4hih2i+1）
α=-
hi+1（-hi+1hi-h2i+h2i+1） （h3i+h3i+1+4hi+1h2i+4hih2i+1）
意义。
②直接求解法： 直接构造基于非等距网
格的差分方法， 网格的选取比 JT 方法广泛，
可以选择不光滑的网格，也能保证有效的解。
下面以模型方程 f=cos （ωx+准）（其中 ω=
4，准=1）为例来说明在三种不同的非等距网格
SINH2-BASED MESH，TANH2-BASED MESH
和 RANDOM-BASED MESH［7］下 两 种 方 法 的
［3］Y. MA，D. FU， Super compact finite difference method with uniform and non -uniform grid system ［J］， Proc Sixth Intern Symp on CFD， Lake Tahoe， Nevada， 1435-1439， 1995.
图二
- 42 -
网 格 点 数 （k） 8×8
16×16 32×32 42×42
表一
紧 致 差 分 方 法 （r.m.s） 1.606e-3 4.859e-5 4.151e-6 3.875e-6
区域供热 2010.1 期
经 典 差 分 方 法 （r.m.s） 2.895e-3 9.510e-5 5.531e-6 4.076e-6
式中未知系数应满足如下关系：
Ai+Bi+Ci=0
-Ai hi+Ci hi+1=0 Ai hi2+Ci h2i+1=2！（αi +1+βi） -Ai hi3+Ci h3i+1=3！（αi hi+βi hi+1）
Ai hi4+Ci
h4i+1=
4！ 2！
（αi
h2i+βi h2i+1）
（1.2）
截断误差为：
1 基于非等距网格紧致差分格式
1.1 紧致差分格式构造及截断误差
计 算 网 格 为 ｛xk｝，网 格 间 距 为 hk=xk-xk-1。 非等距网格下二阶偏导数的紧致差分格式为：
αi fi-1″+ fi″+βi fi+1″=Ai fi-1+Bi fi+Ci fi+1 （1.1） 根据泰勒展式系数匹配的方法 ［5］，（1.1）
区域供热 2010.1 期
坠f 坠x
=
坠f/坠ξ 坠x/坠ξ
，
坠2f 坠x2
=
坠2f 坠ξ2
-
坠f 坠ξ
坠2x 坠ξ2
（ 坠x ）2
坠ξ
（1.6）
公式中包 含 了 网 格 的 转 换 过 程 ξ=ξ（x），
公式右边的每个导数都用等距的有限差分方
法进行计算。 这种方法的要求是所选取的网
格必须充分光滑才可以保证 坠x/坠ξ，坠2x/坠ξ2 有
ε=［-Ai
hi5+Ci
h5i+1+
5！ 3！
（αi
h3i-βi
h3i+1）］
h （5） i
5！
（1.3）
以上构造的差分格式适用于内结点，对
于边界点， 可采用降低精度的紧致差分逼近
格式为：
f1″+β1 f2″=A1 f1+B1 f2+C1 f3
（1.4）
截断误差为：
ε=［Bi
h4i+1+Ci
（hi+1+hi+2）4-
（3.1）
在初始时刻，速度和压力均为零，任意时
刻的压力边界条件被给。 驱动力为：
fx = - costsin2πxsinπycosπy + πsin2tsin3πx cosπxsin2πy + πsin tcosπy - π2 sint （ 6sin2πx -
Re
2cos2πx）sinπycos πy
姨 其中：r.m.s=
1 k
∑k（unumerical-uexact）2
4．结论
紧致差分格式是一种基架点少， 精度高
的差分格式， 目前在复杂的多尺度数值模拟
中得到了较广泛的应用。 本文构造了基于
非等距网格的紧致差分格式， 对计算区域进
行直接数值求解， 该方法克服了传统的雅克
比变换法在网格剧烈变化时降低精度的缺
例：瞬时的周期性流动
考虑具有固定边界的单位方腔内短暂流
动。 流动是由于周期性的体力驱动的，该问题
的 精 确 解 为 ［9］：
u（x，y，t）=-sin tsin2 πxsinπxcosπy
ν（x，y，t）=sin tsin πxcosπxsin2πy
p（x，y，t）=sin tsin πxcosπy
fy = costsinπxcosπxsin2πy + πsin2tsin2πx sin3πycosπy -πsin tsinπxsinπy + π2 sintsinπx
Re
cos πx（6sin2πy-2cos2πy）
（3.2）
从图二和表一可以看出， 本文所构造的
紧致差分方法是一种高效求解 Navier-Stokes
4！ 2！
βh2i+1］
f （4）
i
4！
（1.5）
- 40 -
同理， 可以得出非等距网格下一阶偏导 数 的 紧 致 差 分 格 式 ［6］。 格 式 中 的 系 数 只 与 计 算 网格有关，一旦网格给定，这些系数就固定下 来。 在数值计算中，如果计算网格不随时间变 化，则这些系数只需计算一次，因此这种方法 并不比等距网格下的紧致格式增加计算量。
通过显式投影法可以得到如下半离散形 式 ［8］：
- 41 -
区域供热 2010.1 期
軋ν*=軋νn+Δt［-（ν軋· ）軋νn+ 1 Re
2ν→n+ f軆］
（2.3）
2pn+1= 1 ·軋ν* Δt
（2.4）
軋νn+1=軋νn-Δt pn+1
（2.5）
其中（2.3）式中的对流项（ν軋· ）軋νn 与（2.5）
区域供热 2010.1 期
求解自然对流换热问题的高效方法
吉林市热力有限公司 黄建利
【摘 要】 通过泰勒展式系数匹配的方法构造了基于非等距网格的紧致差分格 式，并将其与投影法相结合，得出了求解自然对流换热问题的新算法，该算法用基于 交错网格的压力泊松方程求解压力项，满足相容性条件。 数值实验表明，紧致差分格 式比经典差分格式能够得到更高精度。
【关键词】紧致差分格式 非等距网格 投影法
0 引言 近 年 来 高 精 度 的 紧 致 差 分 格 式 （CFDS） 引起了国内外普遍重视， 目的是更准确地模 拟复杂流场的流动。 与高精度的经典差分 （FDS） 格 式 相 比 ， 在 相 同 精 度 下 ， 紧 致 差 分 格 式有着风格基架点少和分辨率高的优点。 20 世纪 90 年代以来，紧致差分格式的研究在国 内外都得到了长足的进展。 Lele 给出不限于 三 点 的 对 称 型 紧 致 格 式 ［1］，傅 德 蕙 和 马 延 文 于 1992 年提出迎风紧致格式 ［2］， 又于 1995 年 提 出 了 超 紧 致 格 式 ［3］。 张 涵 信 和 沈 孟 育 于 2003 年将上述紧致格式进一步推广，提出了 广义紧致格式， 前面提到的三种紧致差分格 式 都 是 它 的 特 例 ［4］。 然而，到目前为至大多数紧致差分格式是 基于等距网格的；对于非均匀网格，通常是通 过雅克比变换将其变为计算域的均匀网格，这 在物理域网格变换剧烈的地方会产生较大的 误差。 本文构造出直接针对非等距网格的紧致 差分格式，该格式不但具有较高的精度，还可 适应空间网格的剧烈变化。 将其与投影方法结 合，得出求解二维不可压 Navier-Stokes 方程的 高效算法。 数值计算的结果表明，该方法是求 解 Navier-Stokes 方程的有效方法。
） +Su
（2.8a）
坠p 坠y
=- 坠ν 坠t
+1 Re
坠ω 坠x
-（u
坠ν 坠x
+ν 坠ν 坠y
） +Sν
（2.8b）
压 力 泊 松 方 程 求 解 的 关 键 在 于 其 Neu-
mann 边界条件与源项间要满足相容性条件，
采用交错网格解决了边界处的残余散度问
题，使相容性条件得到满足。
3.数值模拟
方程的方法，与经典的差分格式相比，它具有
更高的精度。 图二：当雷诺数 Re=104，无量纲时间 t=0.8
时，紧致差分方法计算结果，（a）压力 p 分布，
（b）横向速度 u 分布，（c）纵向速度 v 分布。 表一： 当雷诺数 Re=104， 无量纲时间 t=
0.8 时，紧致差分方法与经典差分方法对比
区别与联系。 二阶导数精确解与数值解的误
差为 ε2。
从图一可以看出： 对于变化剧烈的网格
图一 雅克比变换法和直接求解法在三种不同网格下求解方程二阶偏导数的对比
RANDOM-BASED MESH，直接求解方法要明 显优于雅克比变换法， 后者甚至不能够得到 有 效 的 解 ； 对 于 光 滑 网 格 TANH2 -BASED MESH， 直接求解方法略好于雅克比变换法； 对于网格 SINH2-BASED MESH，直接求解方 法和雅克比坐标变换法得到的精度基本一 致。 所以大多数情况下，选择直接求解法更为 有效。
［2］Y. MA，D. FU， Numerical solution of the incompressible Navier -Stokes equations with an upwind compact difference scheme ［J］， Numer Meth Fiulds， 30:509-521， 1992.
+
坠2ν 坠y2
）-（u
坠ν 坠x
+ν
坠ν 坠y
）
+fy
（2.6b）
为了保证边界处的压力满足连续性方
程，设
ω= 坠ν - 坠u 坠x 坠y
（2.7）
由 连 续 性 方 程 （2.2）、 （2.6a） 和 （2.6b） 可 得：
坠p 坠x
=-
坠u 坠t
-
1 Re
坠ω 坠y
-（u
坠u 坠x
+ν
坠u 坠y
当求解压力时， 压力边界条件采用的是
Neumann 边界条件，由（2.1）分 别 展 开 为 x 方
向和 y 方向上的动量方程得：
坠p 坠x
=- 坠u 坠t
+
1 Re
（
坠2u 坠x2
+
坠2u 坠y2
）-（u
坠u 坠x
+ν 坠u 坠y
）
+fx
（2.6a）
坠p 坠x
=-
坠ν 坠t
+
1 Re
（
坠2ν 坠x2
β=-
（h2i+1+hi+1hi-h2i）hi （h3i+h3i+1+4hi+1h2i+4hi h2i+1）
方程（1.4）中系数：
A1=
6 h22-h32
，
B1=-
6 （h2-h3）h3
C1=-
6 h2h3（h22-h32）
，
β1=
（2h2+h3） （h2-h3）
参考文献
［1］S. K. Lele， Compact finite difference schemes with spectral -like resolution ［J］， Journal of Computational Physics， 103: 16-42， 1992.
式中的压力梯度项 p 可以通过求解一阶偏 导数的四阶紧致差分格式求得，（2.3） 式中的 扩散项 2ν軋 可以通过求解二阶偏导数的三阶 紧致差分格式求得， 压力泊松方程采用二阶 中心差分格式离散求解。 由（2.3）式得到先验
速度軋ν* 不满足连续性方程，在（2.5）式中被修
正。
2.2 边界条件的处理
2．求解 N-S 方程的高精度算法 2.1 投影法
无量纲化后的二维纳维斯托克斯方程的 守恒形式和连续性方程为：
坠ν軋+（ν軋· ）ν軋+ p= 1 2ν軋+ f軆
坠t
Re
（2.1）
·ν軋=0
（2.2）
其 中 ν軋=（u，ν）为 速 度 矢 量 ，p 为 压 力 ， f軆= （fx，fy）为源项，Re 为雷诺数。
点，并比较了在三种不同的非等距网格下，直
接数值求解法与雅克比变换法的区别与联
系，为方法的选择提供了理论依据。 将紧致差
分格式与投影法相结合，构造了求解 Navier-
Stokes 方程的新方法。 通过数值实验表明该
方法比经典有限差分具有更高的精度， 是一
种求解 Navier-Stokes 方程的有效方法。
附录
方程（1.1）中系数：
Ai=
12hi+1 （h3i+h3i+1+4hi+1h2i+4hih2i+1）
Bi=
12 （h2i+1+3hi+1hi+h2i）
Ci=百度文库
12hi+1 （h3i+h3i+1+4hi+1h2i+4hih2i+1）
α=-
hi+1（-hi+1hi-h2i+h2i+1） （h3i+h3i+1+4hi+1h2i+4hih2i+1）
意义。
②直接求解法： 直接构造基于非等距网
格的差分方法， 网格的选取比 JT 方法广泛，
可以选择不光滑的网格，也能保证有效的解。
下面以模型方程 f=cos （ωx+准）（其中 ω=
4，准=1）为例来说明在三种不同的非等距网格
SINH2-BASED MESH，TANH2-BASED MESH
和 RANDOM-BASED MESH［7］下 两 种 方 法 的
［3］Y. MA，D. FU， Super compact finite difference method with uniform and non -uniform grid system ［J］， Proc Sixth Intern Symp on CFD， Lake Tahoe， Nevada， 1435-1439， 1995.
图二
- 42 -
网 格 点 数 （k） 8×8
16×16 32×32 42×42
表一
紧 致 差 分 方 法 （r.m.s） 1.606e-3 4.859e-5 4.151e-6 3.875e-6
区域供热 2010.1 期
经 典 差 分 方 法 （r.m.s） 2.895e-3 9.510e-5 5.531e-6 4.076e-6
式中未知系数应满足如下关系：
Ai+Bi+Ci=0
-Ai hi+Ci hi+1=0 Ai hi2+Ci h2i+1=2！（αi +1+βi） -Ai hi3+Ci h3i+1=3！（αi hi+βi hi+1）
Ai hi4+Ci
h4i+1=
4！ 2！
（αi
h2i+βi h2i+1）
（1.2）
截断误差为：
1 基于非等距网格紧致差分格式
1.1 紧致差分格式构造及截断误差
计 算 网 格 为 ｛xk｝，网 格 间 距 为 hk=xk-xk-1。 非等距网格下二阶偏导数的紧致差分格式为：
αi fi-1″+ fi″+βi fi+1″=Ai fi-1+Bi fi+Ci fi+1 （1.1） 根据泰勒展式系数匹配的方法 ［5］，（1.1）
区域供热 2010.1 期
坠f 坠x
=
坠f/坠ξ 坠x/坠ξ
，
坠2f 坠x2
=
坠2f 坠ξ2
-
坠f 坠ξ
坠2x 坠ξ2
（ 坠x ）2
坠ξ
（1.6）
公式中包 含 了 网 格 的 转 换 过 程 ξ=ξ（x），
公式右边的每个导数都用等距的有限差分方
法进行计算。 这种方法的要求是所选取的网
格必须充分光滑才可以保证 坠x/坠ξ，坠2x/坠ξ2 有
ε=［-Ai
hi5+Ci
h5i+1+
5！ 3！
（αi
h3i-βi
h3i+1）］
h （5） i
5！
（1.3）
以上构造的差分格式适用于内结点，对
于边界点， 可采用降低精度的紧致差分逼近
格式为：
f1″+β1 f2″=A1 f1+B1 f2+C1 f3
（1.4）
截断误差为：
ε=［Bi
h4i+1+Ci
（hi+1+hi+2）4-
（3.1）
在初始时刻，速度和压力均为零，任意时
刻的压力边界条件被给。 驱动力为：
fx = - costsin2πxsinπycosπy + πsin2tsin3πx cosπxsin2πy + πsin tcosπy - π2 sint （ 6sin2πx -
Re
2cos2πx）sinπycos πy
姨 其中：r.m.s=
1 k
∑k（unumerical-uexact）2
4．结论
紧致差分格式是一种基架点少， 精度高
的差分格式， 目前在复杂的多尺度数值模拟
中得到了较广泛的应用。 本文构造了基于
非等距网格的紧致差分格式， 对计算区域进
行直接数值求解， 该方法克服了传统的雅克
比变换法在网格剧烈变化时降低精度的缺
例：瞬时的周期性流动
考虑具有固定边界的单位方腔内短暂流
动。 流动是由于周期性的体力驱动的，该问题
的 精 确 解 为 ［9］：
u（x，y，t）=-sin tsin2 πxsinπxcosπy
ν（x，y，t）=sin tsin πxcosπxsin2πy
p（x，y，t）=sin tsin πxcosπy
fy = costsinπxcosπxsin2πy + πsin2tsin2πx sin3πycosπy -πsin tsinπxsinπy + π2 sintsinπx
Re
cos πx（6sin2πy-2cos2πy）
（3.2）
从图二和表一可以看出， 本文所构造的
紧致差分方法是一种高效求解 Navier-Stokes
4！ 2！
βh2i+1］
f （4）
i
4！
（1.5）
- 40 -
同理， 可以得出非等距网格下一阶偏导 数 的 紧 致 差 分 格 式 ［6］。 格 式 中 的 系 数 只 与 计 算 网格有关，一旦网格给定，这些系数就固定下 来。 在数值计算中，如果计算网格不随时间变 化，则这些系数只需计算一次，因此这种方法 并不比等距网格下的紧致格式增加计算量。
通过显式投影法可以得到如下半离散形 式 ［8］：
- 41 -
区域供热 2010.1 期
軋ν*=軋νn+Δt［-（ν軋· ）軋νn+ 1 Re
2ν→n+ f軆］
（2.3）
2pn+1= 1 ·軋ν* Δt
（2.4）
軋νn+1=軋νn-Δt pn+1
（2.5）
其中（2.3）式中的对流项（ν軋· ）軋νn 与（2.5）
区域供热 2010.1 期
求解自然对流换热问题的高效方法
吉林市热力有限公司 黄建利
【摘 要】 通过泰勒展式系数匹配的方法构造了基于非等距网格的紧致差分格 式，并将其与投影法相结合，得出了求解自然对流换热问题的新算法，该算法用基于 交错网格的压力泊松方程求解压力项，满足相容性条件。 数值实验表明，紧致差分格 式比经典差分格式能够得到更高精度。
【关键词】紧致差分格式 非等距网格 投影法
0 引言 近 年 来 高 精 度 的 紧 致 差 分 格 式 （CFDS） 引起了国内外普遍重视， 目的是更准确地模 拟复杂流场的流动。 与高精度的经典差分 （FDS） 格 式 相 比 ， 在 相 同 精 度 下 ， 紧 致 差 分 格 式有着风格基架点少和分辨率高的优点。 20 世纪 90 年代以来，紧致差分格式的研究在国 内外都得到了长足的进展。 Lele 给出不限于 三 点 的 对 称 型 紧 致 格 式 ［1］，傅 德 蕙 和 马 延 文 于 1992 年提出迎风紧致格式 ［2］， 又于 1995 年 提 出 了 超 紧 致 格 式 ［3］。 张 涵 信 和 沈 孟 育 于 2003 年将上述紧致格式进一步推广，提出了 广义紧致格式， 前面提到的三种紧致差分格 式 都 是 它 的 特 例 ［4］。 然而，到目前为至大多数紧致差分格式是 基于等距网格的；对于非均匀网格，通常是通 过雅克比变换将其变为计算域的均匀网格，这 在物理域网格变换剧烈的地方会产生较大的 误差。 本文构造出直接针对非等距网格的紧致 差分格式，该格式不但具有较高的精度，还可 适应空间网格的剧烈变化。 将其与投影方法结 合，得出求解二维不可压 Navier-Stokes 方程的 高效算法。 数值计算的结果表明，该方法是求 解 Navier-Stokes 方程的有效方法。
） +Su
（2.8a）
坠p 坠y
=- 坠ν 坠t
+1 Re
坠ω 坠x
-（u
坠ν 坠x
+ν 坠ν 坠y
） +Sν
（2.8b）
压 力 泊 松 方 程 求 解 的 关 键 在 于 其 Neu-
mann 边界条件与源项间要满足相容性条件，
采用交错网格解决了边界处的残余散度问
题，使相容性条件得到满足。
3.数值模拟
方程的方法，与经典的差分格式相比，它具有
更高的精度。 图二：当雷诺数 Re=104，无量纲时间 t=0.8
时，紧致差分方法计算结果，（a）压力 p 分布，
（b）横向速度 u 分布，（c）纵向速度 v 分布。 表一： 当雷诺数 Re=104， 无量纲时间 t=
0.8 时，紧致差分方法与经典差分方法对比
区别与联系。 二阶导数精确解与数值解的误
差为 ε2。
从图一可以看出： 对于变化剧烈的网格
图一 雅克比变换法和直接求解法在三种不同网格下求解方程二阶偏导数的对比
RANDOM-BASED MESH，直接求解方法要明 显优于雅克比变换法， 后者甚至不能够得到 有 效 的 解 ； 对 于 光 滑 网 格 TANH2 -BASED MESH， 直接求解方法略好于雅克比变换法； 对于网格 SINH2-BASED MESH，直接求解方 法和雅克比坐标变换法得到的精度基本一 致。 所以大多数情况下，选择直接求解法更为 有效。