1
OpenCFD‐EC 理论手册
李新亮
中国科学院力学研究所
2
版权说明
OpenCFD‐EC (Open Computational Fluid Dynamic code for Engineering
Computation) 是可压缩 Navier‐Stokes 方程多块结构网格有限体积求解器。主要
用于计算工程问题。该软件是作者 OpenCFD
计划
项目进度计划表范例计划下载计划下载计划下载课程教学计划下载
的一部分。
作者的 OpenCFD 计划包含本软件（OpenCFD‐EC）及另外一个软件
OpenCFD‐SC (Open Computational Fluid Dynamic code for Scientific Computation，即
原先的 Hoam‐OpenCFD)。OpenCFD‐SC 是可压缩 N‐S 方程高精度差分求解器，主
要用来计算科学问题。
OpenCFD‐EC 及 OpenCFD‐SC 是由作者李新亮开发的计算流体力学（CFD）
软件。目前该软件尚未开源，未获作者许可，获得该软件代码的用户请不要泄露
给第三方。
如在科研或工程计算中使用 OpenCFD‐EC及 OpenCFD‐SC，请务必在论文及
报告
软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载
中进行标注。即在文中致谢并引用作者的论文（作者的论文可使用 SCI搜索
作者=“LIXL AND FUDX”即可）。
参考文献：
J. Blazek: Computational Fluid Dynamics: Principles and Applications, Elsevier 2005
傅德薰 主编 《计算空气动力学》
任玉新等 《计算流体力学基础》
李新亮 《计算流体力学》 课件
（下载地址： www.cfluid.com (流体中文网） -> “流体论坛” ->“ CFD基础理”；
http://cid-1cc0dcbff560c149.skydrive.live.com/browse.aspx/.Public ）
朱自强等： 《应用计算流体力学》
E. F. Toro: Riemann Solvers and numerical methods for fluid dynamics
第一章 二维求解器
§1.1 基本方程
本章以二维问题为例，说明本软件求解器的理论方法。
控制方程为二维守恒型 Navier‐Stokes方程
v v
t x y x y
1 21 2 F (U) F (U)F (U) F (U)U （1.1.1）
本软件采用无量纲形式计算。密度、速度、温度和压力分别采用特征密度 ，
特征速度U ，特征温度T 及动压 2U 无量纲化。
3
无量纲 N‐S方程的具体形式见傅德薰《计算空气动力学》。
在以(I,J)点为中心的控制体上，对 N‐S方程进行积分，得到如下方程
1 1 0IJ v
IJ IJ
ds ds
U F n F nt （1.1.2）
左端第 1项为时间导数项，第 2项为无粘项，第 3项为粘性项。
控制体的选取主要有网格中心型及节点中心型两种方式（见图 1）。 相对而
言，网格中心型控制体在处理角点、壁面及交接面时具有一定的优势（见图 2），
因而本程序采用网格中心型控制体。
在网格中心型控制体中，物理量（密度、速度、压力等）存储在控制体中心。
本手册约定：下标（i,j）表示网格节点，下标（I,J）表示控制体中心。
图 1.1 网格及控制体示意图
图 1.2 网格中心型及节点中心型控制体积的对比
其中，控制体中心定义为：
1
cr rd
4
即该四边形的重心。
具体计算时，四边形 ABCD的中心可以用三角形 ABD及 BCD的中心加权平
均得到。即：
, ,ABD c ABD BCD c BCD
c
ABD BCD
r r
r
其中，三角形中心的坐标为三个顶点坐标的代数平均。即：
,
1 ( )
3c ABD A B C
r r r r
三角形的面积
公式
小学单位换算公式大全免费下载公式下载行测公式大全下载excel公式下载逻辑回归公式下载
为：
1 1 [( )( ) ( )( )]
2 2ABD AB AD B A D A B A D A
r r x x y y y y x x
图 1.3 示意图： 四边形 ABCD的中心
§1.2 无粘项的离散
1.2.1 积分方程
（1.1.2）式左侧第 2项为无粘项，其表达式如下：
1/2, 1/2, , 1/2 , 1/2
1 1 ( )I J I J I J I J
IJ IJ
ds F n H H H H （1.2.1）
其中 IJ 是以包围中心点（I,J）的控制体（图 1中红色虚线围成四边形）的面积。
具体计算时，可用该四边形对角线的向量进行叉乘获得。
5
图 1.4 示意图：四边形的面积= 1
2
AC DB
图 1.5 控制体示意图
穿过控制体右边界面（图 1.5中箭头所指的界面）的无粘通量的表达式为
1/2, 1 1 2 2 1/2,1/2,I J I JI JH n F n F s （1.2.2）
其中：
2
1
( )
u
u p
uv
E p u
F 2 2
( )
v
uv
v p
E p v
F （1.2.3）
1 2( , )
Tn nn 为控制体右边界面的法方向， 1/2,I Js 为该界面的长度（在三维情况
下为面积）。控制体其他边界面通量的计算方法与右边界的计算方法相同。
无粘通量计算的关键是计算(I+1/2,J)点处的通量 1/2,I JF 。如采用中心型格式，
该值可右左、右两侧值平均获得，即：
1/2, , 1,( ) / 2I J I J I J F F F （1.2.4）
或
1/2, , 1,( ) / 2I J I J I J F F U U （1.2.5）
直接使用(1.2.4)或（1.2.5）会带来不稳定。为了消除数值振荡，通常在（1.2.4）
或（1.2.5）式右端加入人工粘性项，以便于稳定程序。
即使添加了人工粘性，采用中心格式计算的Mach数仍不会太高（例如不超
过 3）。本程序采用迎风格式计算无粘项。迎风格式需要把无粘通量按照特征方
向进行分解（以便于确定“迎风”方向）。根据分解形式不同，可分为流通矢量
分裂(FVS)及通量差分分裂(FDS)两类。
常用的流通矢量分裂(FVS)包括 Lax‐Friedichs 分裂，Steger‐Warming 分裂，
Van Leer分裂 及 Liou‐Steffen分裂等。 Lax‐Friedichs分裂最为简单，且能保证函
数的光滑性，但其耗散最大。Steger‐Warming分裂耗散较小，但无法保证函数音
6
速点的光滑性。Van Leer分裂克服了 Steger‐Warming分裂的缺点，但其耗散略大。
Liou‐Steffen分裂即 AUSM方法，是在 Van Leer分裂基础上，将压力项单独处理
的一种方法，在目前得到了广泛的应用。
通量差分分裂(FDS)通过 Riemann 解来计算数值通量，因为更好地利用了双
曲方程的特征方向，因而其激波捕捉能力更强，数值振荡更小。但其计算量要大
于流通矢量分裂。常用的 FDS 方法包括精确 Riemann 解(Godnov 方法），Roe 近
似 Riemann解及 HLL/HLLC近似 Riemann解等。
本软件提供 Steger‐Warming FVS，HLL/HLLC 及 Roe FDS 四种通量分裂方法。
1.2.2 通量差分分裂(FDS)方法
先重构出控制界面(I+1/2,J)上的值。在迎风型方法中，利用（偏）左侧及（偏）
右侧网格点，重构出该点的左值 1/2LiU 及右值 1/2RiU 。可利用原始变量、守恒变量
或特征变量进行重构。 采用特征变量进行重构数值振荡最小，但计算量最大。
本软件采用守恒变量进行重构。
重构方法与有限差分法构造差分格式的形式完全相同。本软件采用 2阶或 3
阶 3阶迎风偏斜重构（见 2.3节）。
重构出 1/2LIU 及 1/2RIU 以后，求解 Riemann问题，即一维 Euler方程及如下初
值条件：
1/2 1/2
1/2
L
I I
R
I
U if x x
U
U otherwise
即可得到穿过(i+1/2,j)界面的通量 1/2,i jH 。
需要注意，由于本问题属于二维问题，可通过坐标旋转，转化成为垂直界面
方向的（扩展）一维 Riemann问题，然后求解。具体公式见 任玉新 《计算流体
力学基础》 （p164‐168).
本软件目前使用 HLL/HLLC 积分型近似 Riemann 求解器及 Roe 微分型近似
Riemann 求解器计算通量。具体公式见 Toro: Riemann Solvers and numerical
methods for fluid dynamics。
1.2.3 流通矢量分裂(FVS)方法
与 FDS方法相同，首先重构出界面（I+1/2,J）处的函数值： 1/2LIU 及 1/2RIU 。
利用 Steger‐Warming流通矢量分裂，将控制界面上的通量分解为正通量及
负通量：
1/2, 1/2, 1/2,I J I J I JH H H
具体公式见 傅德薰《计算空气动力学》158‐159页 （请留意书中的印刷错误）。
7
（2.5）式中的系数为：
其中正通量用 1/2LIU 计算出，负通量用 1/2RIU 计算出，即：
1/2, 1/2, 1/2, 1/2, 1/2,( ) ( )
L R
I J I J I J I J I JH H U H U
与 FDS方法相同，对于二维问题，需要进行坐标旋转，转化为“扩展一维问题”。
1.2.4 重构
重构是利用节点上的控制体平均值 iU 计算出控制体界面（半点）上的值
1/2iU 。应当注意重构与插值的区别（见作者的《计算流体力学》课件）。对于不
超过 2阶精度的重构，由于节点值与控制体平均值相同，因此重构形式与插值相
同。但对于高阶方法二者不同。 当网格充分光滑时，有限体积法与有限差分法
的重构形式相同， 有限差分法中的公式可直接应用。
本软件使用的重构方法有如下四种：
1）2阶 NND格式 (NND2)
1/2 1 1
1/2 1 1 2 1
1/ 2 min mod( , )
1/ 2 min mod( , )
L
I I I I I I
R
I I I I I I
U U U U U U
U U U U U U
（1.2.6）
2）3阶迎风偏斜格式(UCD3)
1/2 1 1
1/2 1 2
( 5 2 ) / 6
(2 5 ) / 6
L
I I I I
R
I I I I
U U U U
U U U U
（1.2.7）
3） 3阶WENO格式 (WENO3)
2 2
1 1 2 1
2
1 2
1 2
1/2 1 1 2 1
( ) , ( )
/ ( ) , / ( )
1/ 3, 2 / 3
( / 2 3 / 2 ) ( / 2 / 2)
I I I I
k k k k k
L
I I I I I
IS U U IS U U
C IS
C C
U U U U U
（1.2.8）
4） 2阶MUSCL格式 (MUSCL2)
1/2 1 11/ 2min mod( ,2( ))
L
I I I I I IU U U U U U （1.2.9）
1.2.5 原始变量、守恒变量与特征变量的重构
本软件支持使用原始变量、守恒变量或特征变量来重构出半点 I+1/2处的值。
所谓原始变量，指密度、速度及压力。本软件中为网格内的平均密度、平均
速度及平均压力。在二阶精度范围内，可认为是网格中心处（I点）的值。因为
8
这些量可通过实验直接观测，是最直观的变量，因此称为“原始变量”。原始变
量为：
I
I
u
v
p
q
所谓守恒变量，指的是密度（质量密度）、动量密度及能量密度。本软件仍然
使用的是网格中心的值（更严格的讲上是网格内的平均值）。因为这些量满足质
量守恒、动量守恒及能量守恒关系式（1.1），因此成为“守恒变量”。守恒变量
为：
I
I
u
v
E
U
其中能量密度的表达式为：
2 21 ( )
1 2
pE u v
为单位体积内的总能（内能+动能）。
特征变量是根据双曲方程特征理论，通过对守恒变量进行线性变换得到的一
组新变量。即：
1/2I I IV = S U （1.2.10）
其中矩阵 S为 Jocabi矩阵
FA
U
的左特征矩阵。其具体表达式见傅德薰《计算
空气动力学》第 4.7节 （158‐162页， 需留意其中的印刷错误，公式以本软件
为准）。
对于常系数偏微方程组，根据特性理论，采用特征变量可以把方程组解耦，
变换成若干独立的单变量微分方程求解。对于像 N‐S 方程这样复杂的偏微分方
程，仍可采用冻结系数法，利用 I+1/2点的特征矩阵 1/2IS 将周围几个点（本程序
中使用 I-1, I, I+1三个点）上的守恒变量通过（2.5）式变换为特征变量 IV 。 使
用特征变量的好处是各物理量可尽量避免相互干扰，减少数值振荡。
本软件支持使用原始变量、守恒变量及特征变量来重构 I+1/2点上的物理量。
对于每种方法，用户均可从 NND,WENO,UCD,MUSCL四种格式中选择一种进行重
构（见公式(2.1)‐(2.4))。
例如，用户选择 2阶 NND格式重构 I+1/2点左值的公式为：
1/2 1 11 / 2 min mod( , )LI I I I I Iq q q q q q （原始变量重构）
1/2 1 11/ 2min mod( , )LI I I I I IU U U U U U （守恒变量重构）
9
1/2 1 11/ 2min mod( , )
L
I I I I I IV V V V V V （特征变量重构）
然后，可利用 U与 q的关系式或 U与 V的关系式（ 1U S V ）计算出 I+1/2点的
U。
§1.3 粘性项的离散
3.1 粘性通量的表达式
方程（1.2）最后一项为粘性项，其表达式如下：
1/2, 1/2, , 1/2 , 1/2
1 1 ( )v v v vv I J I J I J I J
ij ij
ds F n H H H H （1.3.1）
穿过控制体右表面粘性通量的表达式为
1/2, 1 1 2 2 1/2,1/2,
v v v
i j I JI J
H n F n F s （1.3.2）
11
1 12
11 12
0
vF
Tk u v
x
12
2 22
12 22
0
vF
Tk u v
y
（1.3.3）
2
Re 3
ji k
ij ij
j i k
uu u
x x x
（1.3.4）
Pr Re
pck
从表达式中可以看出，只需计算出(I+1/2,J)点处的导数 ,i
j j
u T
x x
，就可计算出穿过
右表面的粘性通量。穿过其他三个表面的粘性通量同样计算。
3.2 导数的计算方法
在（3.4）‐（3.5）中，需要计算物理量的导数
x
和 y
（为速度或温度）。
计算方法有两种选择。
10
方法 1： 利用控制体上的 Gauss公式
图 1.6 计算导数的控制体（红色虚线）
利用以（I+1/2,J）为中心的控制体（图 1.6中红色虚线围成）的高斯积分公
式，可得到该导数的近似值。根据高斯公式：
1/2, 1/2,I J I J
ds d
n
假设控制体内梯度均匀分布，有：
1/2,1/2,
1
I JI J
ds
n
利用中点公式将积分展开，即可得(I+1/2,J)点的导数值。
方法 2： 利用 Jocabian变换
图 1.7 示意图：利用周围点信息计算(I+1/2,J)点的导数值
在有限差分法中，通常利用 Jocabian变换计算曲线坐标系下的导数值，有限
体积法也可同样计算。计算公式为：
11
x xx
（1.3.5）
y yy
（1.3.6）
其中 Jocabian系数的计算方法为：
x
y
x
y
Jy
Jx
Jy
Jx
1
1( , ) ( )
( , )
x yJ x y x y
（1.3.7）
利用（I+1/2,J）周围 4个点（图 3.2中红色的点）的信息，可以计算出上式
中的量：
1, ,
1, 1 1,
I J I J
i j i j
（1.3.8）
其中 可以是函数（速度或温度），也可以是坐标 x或 y.
（i+1,j）及（i+1,j+1）点上的坐标值是已知的，而这些点上的函数值可用周
围格心点上的函数值插值获得。
利用（3.6）‐（3.9）即可计算(I+1/2,J)点上的导数值。通过对比发现，该方
法得到的公式与方法 1完全相同。
§4 网格的多块对接
为了处理更为复杂的几何形状，本软件支持多块对接网格。该网格是由多块
网格对接而成，每块网格均为（相对光滑）的结构网格。
12
图 1.8 对接网格示意图
对接网格如图 1.8所示，对接线上的计算与流体内点的计算方式相同。需
要指出的是，由于对接线两侧的网格光滑性通常无法保证，因此对接线上通常
采用 2阶格式计算。 如采用高阶格式，需考虑网格的非光滑性，形式比较复杂。
在程序中，采用“虚网格点”法，实现不同区域之间的数据交换。设流场
的网格点为 i=1,2.....Nx, j=1,2....Ny; 定义物理量在 i=0,...Nx+1, j=0,.....,Ny+1。 其中
i=0及 Nx+1; j=0及 Ny+1的点为“虚网格点”，其物理采用相邻区域内的信息。
虚网格点上的信息根据网格连接文件（bc2d_1.in)的指示，从相邻 Block中获
取。
网格连接信息文件 bc2d_1.in 示例如下：
(算例： 30P30N三段翼型绕流）
2D BC file for multi‐block mesh: 30P30N
# Blocks
60 总块数
Block 1 第 1块
Subface 该块的子面数
4
f_no, face, istart,iend, jstart, jend, neighb, subface orient
1 1 1 1 1 65 ‐20 0 0
2 4 1 129 65 65 ‐20 0 0
3 2 1 129 1 1 17 3 2
4 3 129 129 1 65 2 2 2
Block 2
Subface
4
f_no, face, istart,iend, jstart, jend, neighb, subface orient
1 4 1 33 65 65 ‐20 0 0
2 1 1 1 1 65 1 4 2
13
3 2 1 33 1 1 18 3 2
4 3 33 33 1 65 3 3 2
Block 3
Subface
4
f_no, face, istart,iend, jstart, jend, neighb, subface orient
1 4 1 33 65 65 ‐20 0 0
2 2 1 33 1 1 19 2 2
3 1 1 1 1 65 2 4 2
4 3 33 33 1 65 4 5 2
Block 4
Subface
5
f_no, face, istart,iend, jstart, jend, neighb, subface orient
1 4 17 41 65 65 ‐20 0 0
2 4 1 17 65 65 ‐20 0 0
3 2 1 41 1 1 20 3 2
4 3 41 41 1 65 5 2 2
5 1 1 1 1 65 3 4 2
......
f_no : 子边界号
face: 子边界类型：1 i‐; 2 j‐ ; 3 i+ ;4 j+
istar, iend: 起始、终止的 i值；
jstar, jend起始、终止的 j值;
neighb: 相邻的 Block号
subface: 相邻的子边界号
ori: 连接方向 （2 顺序连接； 1，3 逆序连接）
14
§4 二维程序测试算例
4.1在随机网格上计算导数
测试算例：
令
( , ) sin( )cos( )u x y x y ，
计算网格公式如下：
, 1
, 2
1 0.4
1 1
1 0.4
1 1
i j
i j
ix r
N N
jy r
N N
其中 r1和 r2是（‐1,1）之间的随机数，网格点 101*101.
利用§2.2中的方法计算偏导数 ,u u
x y
，与精确解相比计算其误差。
计算网格及局部放大图
最大误差： x方向导数： 0.077; y方向导数 0.065
15
y方向导数的精确解（左）与数值解（右）
4.2 RAE2822 翼型绕流
该算例是考核程序的
标准
excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载
算例，计算工况及网格见 NASA网站：
http://www.grc.nasa.gov/WWW/wind/valid/raetaf/raetaf.html
x/C
-C
p
0 0.2 0.4 0.6 0.8 1
-1.5
-1
-0.5
0
0.5
1
1.5
Experiment
2th finite volume method
WENO5 finite difference method
2th FVM (conservative variables)
3rd FVM
RAE2822 Airfoil
表面系数分布： OpenCFD‐EC v0.98 的计算结果
16
4.3 30P30N三段翼型
计算网格
压力分布及流线
17
x/L
-C
p
0 0.2 0.4 0.6 0.8 1 1.2
-1
0
1
2
3
4
5
OpenCFD-EC 0.98: NND+HLLC+BL
Experimential data
翼型表面的压力分布
18
x/L
-C
p
0 0.2 0.4 0.6 0.8 1 1.2
-1
0
1
2
3
4
5
OpenCFD-EC 0.98: NND+HLLC+BL
Experimential data
OpenCFD-EC 0.98: 3rd upwind+HLLC+BL
翼型表面的压力分布（三阶迎风 vs 二阶 NND）
4.4 前台阶问题 （Ma=3）
计算域及区域分解
19
前台阶问题的压力分布
20
第二章 三维求解器
§2.1 基本方程
控制方程为三维守恒型 Navier‐Stokes方程
3 3v v v
t x y z x y z
1 21 2 F (U) F (U) F (U) F (U)F (U) F (U)U （2.1.1）
在以控制体上，对 N‐S方程进行积分得到如下方程
1 1 0IJK v
IJK IJK
ds ds
U F n F nt （2.1.2）
左端第 1项为时间导数项，第 2项为无粘项，第 3项为粘性项。
本程序采用网格中心型控制体。
图 2.1 以（IJK）点为中心的控制体
1.1 控制体体积、表面积、表面法方向及中心点的计算公式
（Blazek et al. : Computational Fluid Dynamics, Principles and Applications, section
5.1.2)
侧面积及法向量计算公式为：
1/2, , 1, 1, 1, , 1 1, 1, 1 1, ,
1 ( ) ( )
2I J K i j k i j k i j k i j k
S r r r r
S Sn
控制体积的计算公式为：
1 ( )
3 c mm
r S
其中下标 m表示该六面体的侧面。
21
控制体中心点的计算公式为：
3 ( )
4 ( )
c m cm
m
c
c m
m
r S r
r
r S
§2.2 无粘项的离散
积分方程（2.1.2）中，无粘项为：
1 1
m m m
mIJK IJK
ds n S F n F
因此，关键是计算各侧面的流通量。
以(I+1/2,J,K)面的通量为例，说明无粘通量的计算步骤：
step 1: 重构，计算出（I+1/2,J,K）点处的流体变量的左、右值： 1/2 1/2,L RI IU U 。
本软件使用的数值方法（NND,3阶迎风，MUSCL,WENO3）均为三点格式。左
值 1/2LIU 利用（I-1,J,K），（I,J,K），（I+1,J,K）三个点上的值构造； 右值 1/2RIU 利用
（I,J,K），（I+1,J,K），（I+2,J,K）三个点上的值构造。具体公式见（1.2.6）-（1.2.9）。
用户可选择使用原始变量、守恒变量或特征变量进行重构。详见第 1.2.5节。
Step 2: 利用 FVS或 FDS方法，根据 1/2 1/2,L RI IU U 计算出穿过(I+1/2,J,K)点所在界
面的通量。
FVS 方法的计算公式为： 1/2 1/2( ) ( )L RI IF F U F U 。 本软件使用
Steger‐Warming FVS方法，具体公式见 傅德薰《计算空气动力学》4.7.2节。
FDS是通过求解 Riemann问题获知穿过界面的通量。由于使用精确 Riemann
解（Godnov 方法）计算量较大，本软件使用近似 Riemann 解。 包括 HLL/HLLC
及 Roe 近似 Riemann 解。 具体公式见 Toro : Riemann solvers and numerical
methods for fluid dynamics, 第 10章及第 11章。也可参考作者的“计算流体力学”
课件。
利用坐标旋转，将三维 Riemann 问题转化成界面法向的“扩展一维 Riemann
问题”（见作者的课件），以简化计算。以二维问题为例，示意图如下。
图 2.3 示意图： 通过坐标旋转，将二维 Riemann问题转化为“扩展”一维问题
22
具体步骤如下：
图 2.4 控制体及界面的法、切方向
a) 计算出(I+1/2,J,K)所在界面的法方向和切方向（见图 2.1）：
1, 1, 1 1, , 1, , 1 1, 1,( ) ( )i j k i j k i j k i j kn r r r r
1 1, 1, 1 1, ,i j k i j kr r
2 1n
b) 通过坐标变换，将通过 Step 1得到的 1/2 1/2,L RI IU U 转化到（ 1 2, ,n ）坐标系
下。守恒变量U 从原始坐标系到（ 1 2, ,n ）坐标系的变换公式如下：
u
U v
w
E
（2.1.3）
其中：
1 1 2 2, ,nu u u n v u u w u u （2.1.4）
分别为速度在新坐标系三个坐标轴方向的分量。 （2.1.3）显示，由于密度 和
总能量（密度）E是标量，在坐标旋转过程中保持不变。动量（密度） u 是矢
量，坐标变换过程中遵循旋转公式(2.1.4)。
将 1/2 1/2,L RI IU U 带入（2.1.3）‐（2.1.4）式，即可得到新坐标系中的守恒变量
1/2 1/2,
L R
I IU U 。
23
c) 在界面法方向n求解扩展一维 Riemann 问题，即可得到流通量（ 1 2, ,n 坐
标系下的值） 1/2IF 。这里所说的“扩展”一维问题指，问题本身是一维的（所
以变量仅在 n方向变化，在其他两个方向均匀分布），但物理量是“扩展的”。因
为纯一维问题的基本物理量是 , , Tu E ，而本问题的物理量是
, , , , Tu v w E ，多了两个切向动量分量。在求解 Riemann问题中，“扩展”
一维问题与原一维 Riemann问题的解法相同，切向分量的性质与被动标量相同，
不会给求解带来任何不便。具体公式见 Toro书的第 10及 11章。
d) 得到流通量后，为了便于和（I,J,K）点周围其他几个界面的流通量累加以求
出进入该控制体积的总流通量，需通过坐标旋转，将 1/2IF 变回原先的(x,y,z)坐标
系。令
1
2
3
4
5
F
F
F
F
F
F ，
1
2
3
4
5
F
F
F
F
F
F
分别表示（x,y,z）及（ 1 2, ,n ）坐标系下的流通量。
则：
1 1 5 5,F F F F
2 2 3 1 4 2x x xF F n F F
3 2 3 1 4 2y y yF F n F F
4 2 3 1 4 2z z zF F n F F
上式显示，质量通量 1F及能量通量 5F 是标量，不随坐标旋转而变化。而动量通
量 2 3 4, ,F F F 是矢量的分量，随坐标旋转而变化。
Step 3: 将计算出的（I+1/2,J,K）点的流通量 1/2IF 乘以该界面面积，即得到穿过该
界面的流通量： 1/2 1/2I IF S 。
采用同样方法，可求出(I,J+1/2,K)及（I,J,K+1/2)所在界面的流通量。这样，计算无
粘流通量就完成了。
24
§2.3粘性项的计算
1） 控制方程
1 1 ( )v x x y y z z m m
mIJK IJK
ds n n n S v v vF n F F F （2.3.1a)
其中：
0
xx
vx xy
xz
xx xy xz
F
Tu v w k
x
0
yx
vy yy
yz
yx yy yz
F
Tu v w k
y
0
zx
vz zy
zz
zx zy zz
F
Tu v w k
z
（2.3.1b）
其中 / Prk Cp 为（无量纲）导热系数。Cp为定压比热（无量纲），在本程序
的无量纲约定下，该系数为： 2
1
( 1)
Cp
Ma 。Pr为 Prantdl数，对于理想空气，
通常取 Pr=0.77。
（无量纲）粘性系数的计算可采用如下 Sutherland公式：
3/2
1 288.15 110.4
Re 288.15 110.4
T
T
（2.3.2）
其中T TT 为有量纲温度。
25
4 2 ( )
3 3
4 2 ( )
3 3
4 2 ( )
3 3
u v w u v u w
x y z y x z x
u v v u w v w
y x y x z z y
u w v w w u v
z x z y z x y
τ （2.3.3）
为粘性应力张量。
2） 偏导数的计算
从上段的公式可以看出，计算粘性项的关键在于偏导数 i
j
u
x
及 j
T
x
。计算出
导数后，就可求出粘性应力张量及粘性通量。
与二维程序相同，可使用 Gauss 公式，通过将体积分转化为面积分来计算
导数。也可通过 Jacobian 变换来计算导数。二者得到的公式完全相同。写成
Jocabian 变换的形式更为简洁。Jocabian 变换是将物理空间（x,y,z）变换到均匀
网格的计算空间( , , )。由于计算空间采用均匀网格，不妨假设计算空间网格
间距为 1，这样 ( , , ) ( , , )I J K 。及计算空间的坐标值就是网格的索引（下标）
值。
在网格点（I+1/2,J,K）上计算导数的具体公式如下：
1/2, , 1/2, ,
x x x
I J K I J Kx
1/2, ,1/2, ,
y y y
I J KI J Ky
1/2, , 1/2, ,
z z z
I J K I J Kz
（2.3.4）
其中：
1, , , ,
1/2, ,
I J K I J K
I J K
1/2, 1/2, 1/2, 1/2,
1/2, ,
I J K I J K
I J K
1/2, , 1/2 1/2, , 1/2
1/2, ,
I J K I J K
I J K
(2.3.5)
26
这里为任意函数，可以是速度 , ,u v w、温度T，也可以是坐标 , ,x y z。
物理量（速度、温度）定义在格心，及整数（I,J,K）点。半点处的值可以用
整点处的值插值得到（本程序使用）。例如：
1/2, 1/2, , , 1, , , 1, , 1, 1
1 ( )
4I J K I J K I J K I J K I J K
（2.3.6）
这种插值方法当网格接近平行四边形时具有二阶精度，当网格严重扭曲时为一阶
精度。当网格扭曲严重时，如仍需二阶精度，需用更复杂的双线性插值（任玉新
《计算流体力学基础》5.5.2节）。
（2.3.4）中变换系数的计算公式为（见傅德薰《计算空气动力学》 28页）：
( ), ( ), ( )
( ), ( ), ( )
( ), ( ), ( )
x y z
x y z
x y z
J y z z y J z x x z J x y y x
J y z z y J z x x z J x y y x
J y z z y J z x x z J x y y x
（2.3.7）
其中 Jocabian行列式为：
1
1
( , , )
( , , )
x x x
x y zJ y y y
z z z
（2.3.8）
物理坐标（x,y,z）对计算坐标（即下标）的导数同样采用（2.3.5）式计算。例：
1, , , ,
1/2, ,
I J K I J K
I J K
x x x
1/2, 1/2, 1/2, 1/2,
1/2, ,
I J K I J K
I J K
x x x
实际计算中，需注意格心坐标与格点坐标的区别。本手册中，下标为大写字母的，
表示格心值。下标为小写字母的，表示格点值 （见图 2.4）。
浙公网安备 33021202002483