计算流体力学作业

计算流体力学课程作业

任课教师：魏文礼

姓名：

学号：

指导老师：

目录

1.写出通用方程，并说明如何代表各类守恒方程。 (1)

2.推导流体运动的质量、动量守恒方程。 (2)

3.简述源项线性化、网格划分问题。 (5)

4.用ddxKeTex+S=0，谈谈边界条件如何处理。 (8)

5.用有限体积法离散ρceTet=eexKeTex，并推广到二维、三维问题，写出过程。 (9)

6.从不同角度对流体运动分类。 (12)

7.谈谈物理模型试验与计算流体力学方法的关系。 (12)

8.讨论离散对流项时离散格式的进化过程。 (13)

9.利用幂函数格式离散二维、三维通用方程的离散方程。 (15)

10.解释交错网格的概念。 (15)

11.简述压力校正法解N-S方程的过程。 (16)

12.思考anbvnb′为什么可以省去。 (17)

1.写出通用方程，并说明如何代表各类守恒方程。

答：(1)写出通用方程。

在Cartesian坐标系下单位体积黏性流动N-S方程组微分形式如下：

{

eρ

et

+??(ρV)=0 (1)

eρu

+??(ρuV)=??(μ?u)+

1

μ[

e

(??V)]?

ep

+F x+S mx(2a)eρv

et

+??(ρvV)=??(μ?v)+

1

3

μ[

e

ey

(??V)]?

ep

ey

+F y+S my(2b)

eρw et +??(ρwV)=??(μ?w)+

1

3

μ[

e

ez

(??V)]?

ep

ez

+F z+S mz(2c)eρe

et

+??(ρeV)=??(k?T)?p(??V)+Φ+Q (3)

上述微分形式黏性流动N-S方程组中，式(1)为连续性方程，式(2a)、

(2b)、(2c)分别为x、y、z方向上的动量方程，式(3)为能量方程。

上述方程组中各个方程具有不同变量，代表不同的守恒定律，但他们的形式都十分相似。若引入一个通用的特征变量?，在不同的方程中?代表不同的变量，就可以把它写为通用变量形式。

非定常通用变量N-S方程为：

e(ρ?)

et

+??(ρ?V)=??(Γ???)+S?

若流场中速度等物理量不随时间变化，则e(ρ?)

et

=0，可得定常通用变量N-S方程为：

??(ρ?V)=??(Γ???)+S?

其中，?为通用变量，可代表u、v、w、T等求解变量；Γ?为扩散和热传导系数，S?为方程组源项。

(2)用通用方程代表各类守恒方程

用通用方程代表各类守恒方程是，通用变量在各守恒方程中的取值如表1所示。

表1 在各守恒方程中通用变量的取值

其中，F x、F y、F z为单位质量流体所受体积力在x、y、z方向上的分量，S mx、S my、S mz为单位质量流体的质量源在x、y、z方向上的分量。

2.推导流体运动的质量、动量守恒方程。

答：(1)推导流体力学基本方程组的基本思路

采用Eulerian法，在Cartesian坐标系下，设在时刻t，流场中任意一点(x,y,z)处，取固定不动的六面体单元为控制体，如下图所示。

控制体边长为δx、δy、δz,设流体密度为ρ，某一流动量为?。在δt时间内，从x=x0的δyδz面上流入的流动量为ρ?uδtδyδz，从而取Taylor级数展开

δx}δtδyδz。一阶式，得x=x0+δx的δyδz面上流出的流动量为{ρ?u+e(ρ?u)

ex

同理,δxδz面和δxδy面上流入和流出的流动量?也可得到类似的表达式。

【注：ρ?uδtδyδz=ρ??s??A??=ρ??τ??=?m??】

在δt时间内，通过控制体各表面的流动量?的净增量(对流增量)为：

{e(ρ?u)

ex

+

e(ρ?v)

ey

+

e(ρ?w)

ez

}δtδxδyδz

同时，在控制体内流动量?的净增量(局部增量)为：

e(ρ?)

et

δtδxδyδz

两者之和就是流动量?在δt时间内，在控制体内流动量?的总增量。对它除以ρδtδxδyδz，就得到单位质量流动量?随时间变化的总增量：

1 {e(ρ?)

+

e(ρ?u)

+

e(ρ?v)

+

e(ρ?w)

}

=1

ρ{e(ρ?)

et

+??(ρ?V)}=D?

Dt

+?(1

ρ

Dρ

Dt

+??V)(1)

(2)流体运动的质量守恒方程

对于该控制体，质量守恒定律可表达为：[单位时间内微元体中流体质量的增加]=[同一时间间隔内流入该微元体的净质量]

把单位质量流体?=1代入式(1)，可得Cartesian坐标系下单位质量流体连续方程：

eρet +

e(ρu)

ex

+

e(ρv)

ey

+

e(ρw)

ez

=0

其矢量形式表达式为：

eρ

et

+??(ρV)=0其张量形式表达式为：

eρet +

e(ρu i)

ex i

=0

其中，ρ为流体密度，V为流动速度矢量，u、v、w是其在x、y、z方向上的分量，x i是空间点的坐标，u i为在t时刻x i点的速度分量，i=1,2,3。

对于不可压缩流体，其流体密度为常数，连续性方程可简化为??V=0.

(3)流体运动的动量守恒方程

对控制体分别在三个坐标方向上应用Newton第二定律ma=∑F在流体流动中的表现形式：[微元体中流体动量的增加律]=[作用在微元体上各种力之

和]，并引入Newton 切应力公式及Stokes 表达式，则单位质量动量守恒方程为：

{

Du Dt =eu et +u eu ex +v eu ey +w eu ez =?1ρep ex +F x +1ρ(eτxx ex +eτxy ey +eτxz ez )+S mx Dv Dt =ev et +u ev ex +v ev ey +w ev ez =?1ρep ey +F y +1ρ(eτyx ex +eτyy ey +eτyz ez )+S my Dw Dt =ew et +u ew ex +v ew ey +w ew ez =?1ρep ez +F z +1ρ(eτzx ex +eτzy ey +eτzz ez )+S mz { τxx =2μeu ex +(μ′?23μ)(eu ex +ev ey +ew ez )，τxy =τyx =μ(ev ex +eu ey )τyy =2μev ey +(μ′?23μ)(eu ex +ev ey +ew ez )，τyz =τzy =μ(ew ey +ev ez )τyy =2μew ez +(μ′?23μ)(eu ex +ev ey +ew ez )，τzx =τxz =μ(eu ez +ew ex

) 它的矢量形式表达式为：

DV Dt =eV et +u eV ex +v eV ey +w eV ez =?1ρ

?p +F +τ+S m τ=[τxx τxy τxz

τyx τyy τyz τzx τzy τyy

] 它的张量形式表达式为：

Du i Dt =eu i et +u j eu i ex j =?1ρep ex i +F i +1ρeτi,j ex i

+S mi τij =μ(eu i ex j +eu j ex i )+(μ′?23μ)δi,j eu k ex k

其中，p 为流体压力，F 为单位质量流体所受的体积力，F x 、F y 、F z 是其在x 、y 、z 方向上的分量，F i 为其在时间t 坐标x i 点上的分量，i =1,2,3，τ为流体的黏性应力，τi,j 为其在(i,j )上的张量分量，i =1,2,3和j =1,2,3。μ为流体的动力黏性系数，μ′为膨胀黏性系数。流体黏性系数μ和μ′的大小是由流体分子的性质和分子间的相互作用决定的，它们是温度的函数。由于流体的μ′值往往要比μ值小得多，一般情况下膨胀黏性系数μ′是可以忽略的。上式中δi,j 为

Kronecker 符号，δi,j

={1，i =j 0，i ≠j

；S m 为流体质量源，S mx 、S my 、S mz 是其在

x、y、z方向上的分量，S mi为其在时间t坐标x i点上的分量，i=1,2,3。

对于牛顿流体，流体黏性系数常常可看做是常数，并可忽略膨胀黏性系数μ′，则动量守恒方程式可写为：

{Du

=?

1ep

+F x+υ(

e2u

2

+

e2u

2

+

e2u

2

)+

υe

(

eu

+

ev

+

ew

)+S mx Dv

Dt

=?

1

ρ

ep

ey

+F y+υ(

e2v

ex2

+

e2v

ey2

+

e2v

ez2

)+

υ

3

e

ey

(

eu

ex

+

ev

ey

+

ew

ez

)+S my

Dw Dt =?

1

ρ

ep

ez

+F z+υ(

e2w

ex2

+

e2w

ey2

+

e2w

ez2

)+

υ

3

e

ez

(

eu

ex

+

ev

ey

+

ew

ez

)+S mz

其中，υ=μ

ρ

为流体的运动黏性系数。它的矢量形式表达式为：

DV Dt =?

1

ρ

?p+F+υ?2V+

υ

3

?(??V)+S m

其中，?2=e2

ex2+e2

ey2

+e2

ez2

。

它的张量形式表达式为：

Du i Dt =

eu i

et

+u j

eu i

ex j

=?

1

ρ

ep

ex i

+F i+ν

e

ex j

eu i

ex j

+

υ

3

e

ex j

eu j

ex j

+S mi

3.简述源项线性化、网格划分问题。

答：(1)扩散基本方程源项的线性化

非稳态的扩散方程或导热方程，对一个标量物理变量T可写成：

ρc p eT

=

e

i

(Γ

eT

i

)+S T

其中，S T是单位体积中的净源项，Γ是对应于变量T的扩散系数。

当源项为未知量的函数时，线性化处理比假定源项为常数更合理，线性化处理又是建立线性代数方程所必需的。

把源项局部线性化，亦即假定在未知量微小的变动范围内，源项S T可以表示为该未知量的线性函数。在控制容积P内，它可以表示为以下形式：

S T=S C+S P T P (1)

其中S C为常数部分，S P为S T随T而变化的曲线在P点的斜率，即图1中切线1的斜率。在有限容积方法中，对于一个控制容积P，式(1)为对标量T P的线性方程。

当源项S T为一个非线性函数时，S C和S P两量也将变成标量T P的函数值，此时在数值计算中，我们不得不通过迭代来更新它们的值。在进行S C和S P线性化中，有许多方法可以选择，这其中有一种最佳的方法形式如下式所示：

S T=S?+(dS

dT )

?

(T p?T p?)(2)

在此，假设S T为T的一阶可微函数、星号*表示当前时刻的估计值，如初始值或猜测值。

比较式(1)和式(2)，可得到：

S C=S??(dS

dT )

?

T p?，S P=(dS

dT

)

?

为了保证代数方程迭代求解的收敛，S P≤0，此处要求S P<0；若S P=0，则整个源项为S C。

由代数方程迭代求解的公式

T P=∑a nb T nb+b ∑a nb?S p?V

可见，S p的大小影响到迭代过程中T P的变化，S p的绝对值越大(S P<0)，好像系统的惯性越大，相邻两次迭代之间T P的变化越小，因而收敛速度下降，但有利于克服迭代过程的发散，在图(1)中，如S P取为曲线3的斜率属于此情形。S P的绝对值越小，可使变化率加快，但易引起发散，在图(1)中曲线2即代表此情形。

(2)网格划分问题

网格划分通常有以下两种方法：

(a)外节点法

外节点法即先定节点位置，后定控制容积。这种网格划分问题，在边界上将出现半控制容积。当网格划分不均匀时，节点位置并不落在控制容积的几何中心位置。

(b)内节点法

内节点法即先定控制容积，后定节点位置。这种方法是在边界上附上一层厚度为零的控制容积，代表这个控制容积的边界节点恰好落在边界上。无论网格如何划分都不会出现半控制容积，并且所有节点都位于控制容积的几何中心。

4.用d

dx (KeT

ex

)+S=0，谈谈边界条件如何处理。

答：(1)离散方程

一维稳态导热问题控制方程d

dx (KeT

ex

)+S=0，其中S=S C+S P T P，相应的

离散化方程为a P T P=a E T E+a W T W+b，式中：a E=k e

(δx)e ， a W=k w

(δx)w

，a P=

a E+a W?S P?x，b=S C?x。

此离散方程适用于稳态导热的任何内部节点，为计算一个具体问题，应把边界条件也用于离散方程表示。因为只有离散化方程的个数与带求节点变量的数目相等时，代数方程组才能封闭。

(2)网格划分

采用外节点法划分网格，即先定节点位置，后定控制容积。

(3)构建边界节点附加方程

(a)第一类边界条件

T B已知，不必额外增加边界节点方程，把T B代入邻近节点的代数方程即可。

(b)第二类边界条件

q B已知，可在边界半控制容积内对微分方程积分建立附加方程，热平衡式为：q B?q i+(S C+S P T B)?x=0，规定以进入计算区域的热量为正。

交界面i处的热流为：q i=k i(T B?T I)

(δx)i

，将q i代入上式得：

q B?q i=k i(T B?T I)

(δx)i

+(S C+S P T B)?x=0

q B已知时，方程可整理成：a B T B=a I T I+b，式中：a I=k i

(δx)i

，a B=a I?S P?x，b=S C?x+q B。

(c)第三类边界条件

已知?、t f,将边界热流q B=?(T f?T B)代入上式的B点方程，形式同上：

a B T B=a I T I+b，式中：a I=k i

(δx)i

，a B=a I?S P?x+?，b=S C?x+?T f。

5.用有限体积法离散ρceT

et =e

ex

(KeT

ex

)，并推广到二维、三维

问题，写出过程。

答：(1)一维情况下的离散

非稳态一维热传导方程为：

ρceT

et =e

ex

(KeT

ex

)(1)

通常，我们将假设ρc为常数，c为比热。