热物理过程的数值模拟-计算传热学3

四、非线笥问题迭代式解法的收敛性

每一层次上满足迭代法求解的收敛条件+相邻次间代数方程的系数变化不太大（亦即未知量的变化不太大←多数情形下非线性问题迭代式解法是可以收敛的）。

使相邻两层次间未知量变化不太大的措施： 1、欠松弛迭代 常用逐次欠弛线迭法（SLUR ）：一组临时系数下逐线迭代求解+对所得的解施以欠松弛，再用欠松弛后的解去计算新的系数，常数，以进入下一层次的迭代。

实施：常把欠松弛处理纳入迭代过程，而不是在一个层次迭代完成后再行欠松弛。

t (p n +1) =t (p n ) +ω(

a p

∑a n bt n b (n )

-t p ) a p

a p

t (p n )

(

ω

) t (p n +1) =∑a n b +t n b +b +(1-ω)

ω

a ' p t (p n +1) =∑a n bt n b +b '

a ' p =a p ω, b ' =b +(1-ω)(a p ω) t (p n ) ，用交替方向线迭代法求解这一方程，就实现了SLUR

的迭代求解。为一般化起见，上式中t n b 上没有标以迭代层次的符号（J ，GS 时不相同）。 2、采用拟非稳态法

前面已指出，稳态问题的迭代解法与非稳态问题的步进法十分相似。对于非线性稳态问题，从代数方程的一组临时系数进入到另一组临时系数亦好象非稳态问题前进了一个时间层，非稳态问题的物理特性：系数热惯性越大（a p =ρc ∆v /∆τ↑），温度变化越慢，仿此，对稳态非线性问题，可在离散方程中加入拟非稳态项，以减小未知量托两个层次间的变化，即

由

(n +1) o (n ) (∑a n b -S p ∆V ) t (p n +1) =∑a n bt n b +b ⇒(∑a n b -S p ∆V +a o ) t =∑a bt b +b +a p p n n p t p

(n )

∑a n bt n b +b +a o p t p

o

t (p n +1) =

∑a n b -S p ∆V +a o p

一直进行到t p , t n b 收敛，虚拟时间步∆τ的大小通过计算实践确定。

3、采用Jacobi 点迭代法

中止迭代的判据（该层次迭代）除前述变化率判据外，还可以规定迭代的轮数，例如规定进行4-6次ADI 线迭代就结束该层次上的计算。此时，用收敛速度低的丁迭代也就起到了欠松弛的作用。

五、迭代法的收敛速度 1、收敛速度

对给定的代数方程组（包括是临时系数的情形），采用不同的迭代方法求解时，使一定的初始误差缩小成α倍所需要的迭代轮数K 是不相的。α

e (k ) =k e (0) ≤G k e (0)

2

因为G 是对称的，所以G 所以

=ρ(G T G ) =ρ(G 2) =ρ(G )

2

e (k )

2

≤G k /e (0)

2

2

e (0) =ρ(G )) k e (0)

2

α=e (k )

2

≤(ρ(G )) k , l n α≤kl n ρ(G )

即 k ≥l n α/l n ρ(G ) =-l n α/(-l n ρ(G )) α, ρ(G ) 均

R =-l n ρ(G )

表示迭代法的收敛速度，则

k ≥-l n α/R or R ≥-l n α/k

即所需迭代轮数与收敛速度成反比，收敛速度又与谱半径成反比，收敛速度愈快，迭代轮数愈少。

注意：不同的迭代方法每进行一轮迭代所需的运算次数不同，最终所需的计算时间的多少取决于迭代轮数及每一轮迭代所需的时间。

2、收敛性的定性分析

为什么不同的迭代方法的收敛速度不同，亦即为达到满足一定精度要求所需的迭代轮数不同？以二维常物性、无内热源、稳态导热问题来进行讨论。

∂2t ∂' t

+=0 ∂x 2∂y 2

1B 、C a p t p =a E t E +a w t w +a N t N +a S t S

TB

迭代法需要假定一个初场，例如假定一个均场，从微分方程为看，均匀是其一个解，但却不是所研究问题的解，为什么？因为它虽然满足内部节点上的离散方程，却不满足与边界有关的节点的离散方程（图中红点），即不满足边界条件。所以迭代法的实质是要通过迭代，尽快建立起与边界条件相适应的φ变量场，

j t wB t RB

关键：必须使B 、C 的影响迅速传入计算区域内部，以改进节点φ变量值，尽快与B 、C 相应，B 、C 的影响传入愈快，逼近真解就愈快，收敛就越快！

B 、C 的影响传入计算区域内部的快慢与哪些因素有关？ （1）与迭代方法有关

J 迭代：节点温度的更新均用上轮迭代所得的“旧”值来计算，所以完成一轮迭代后，B 、C

的影响只能传入与边界相邻的一批节点上，即仅可传入一个网格，且扫描方向与收敛快慢无关。要在以后各轮迭代中，B 、C 的影响才由这些节点逐步向内渗透，所以收敛慢。

GS 迭代：假设从左向右扫描，则每做完一轮迭代，左边界和下边界的影响传遍全区域，而右边界的影响只能传入一个网格，且收敛速度受迭代扫描方向的影响。

线迭代：GS 线迭代。自左→右扫描，完成一轮迭代不仅左边界的影响逐步传入，而且在每一列的直接求解中，上、下边界的影响全部传入到该列的各节点上，即一轮迭代使左、上、下边界影响传入全区域，但右边界影响仍仅传入一个网格。

ADI ：一轮迭代包括一次逐行、一次逐列的扫描；所以在每一轮迭代后所有边界的影响均传入计算区域内部，从而加快了收敛速度。

（2）与边界条件的性质有关

t

t

tf

t

1B 、C 规定了边界节点的温度，影响直接传入计算区域内部； 3B 、C 规定了环境温度及定向点位置，

λ∂t

l ∂x

=t -t ∞，对边界温度的限定程度比1B 、C 时

弱，所以对内部的影响也较弱；或将t f 视为外部温度，其对计算区域内部的影响被外部换热热阻削弱，而1B 、C 可视为α→∞或外部热阻→0的极限情况，故3B 、C 的影响比1B 、C 时弱。

2B 、C 仅规定了壁面的钭率，壁温完全不确定，对内部节点温度值的改进提供的信息最少，收敛最慢。

可见，为了提高代数方程迭代解法的收敛速度，应力求使边界条件的影响迅速传入计算区域内部，措施：

①增加迭代解法中直接解法的成分，从点迭代→线迭代→ADI ；

②适当选择扫描的始边，多以1B 、C 或3B 、C 的边界为始边，少以2B 、C （尤其是绝热边界）为扫描始边。

5-4 不规则区域的处理—网格生成技术

如何对不规则区域进行有效的处理，以便于进行传热与流动过程的数值模拟，是近年来计算传热研究中的一个重要课题。

以上讨论的传热过程大都发生在规则而简单的区域中，但许多实际的热传递现象是在不规则

的区域中进行的，例如：

①套片管中肋片的传热 ②渐扩通道中的流动与换热

③环形空间中的自然对流 ④流体外掠管束

以上四种情形中的流动与换热不是直角坐标、圆柱轴对称坐标或极坐标所能方便地予以描述。 虽然有限元法在处理不规则边界方面显示了极大的优越性，但就流动与传热而言，在计算技巧与方法方面，有限差分法都比有限元法成熟。

用有限差分法处理这类问题的方法可以归纳为以下几种。 1、采用阶梯形边界（网格）

用阶梯形边界近似代替四分之一圆弧边界。阶梯形边界（网格）是采用有限差分法计算不规则区域的最普通的方法。

缺点：程序缺少通用性；曲面边界上网格必须划得比较细密（否

则会引起较大误差）。

2、采用区域扩充法

当计算区域的边界不规则程度不很严重时，可以采用区域扩充法，把计算区域扩充为直角坐标，圆柱坐标等常规正交坐标系中易于描述的形状。条件：保证原计算区域的情况不变！’

α,t f

② ①

q

v=实际值

v=1030

流动：把计算区域扩充到图中虚线所示的整个圆形通道，从而可以应用圆柱轴对称坐标系中的控制方程加以描述。

如何保证原来的情况不变？孔板区视为粘性无限大的“流体”，而其余区域的流体粘性值就等于真实值，边界上流速赋为零，计算中零边值将迅速传到孔板区域内，有效地模拟了孔板的存在。

传热：对三种不同的边界条件，具体的处理方式不同 （1）均匀壁温边界条件

令扩充区域中的导热系数为无限大，而扩充后的区域边界温度则等于已知值。 （2）绝热的边界条件 q 令扩充区域中的材料导热系数为零即可实现此条件 （3）均匀热流边界条件

可应用附加源项法来实现，真实边界上均匀热流可以附加源项的形式置于与真实边界相邻的控制容积中去，而扩充区域则处于绝热状态。

周期性二维渐扩、渐缩通道中的换热，倾斜的边界上实际边界 作用有均匀的热流。采用左下所示阶梯形网格，并把计算d

p 区域扩充到一个长方形，以便利用直角坐标系求解。 δ e

a

看一个网格单元的放大后的情形，P 控制容积的附加计算 λ扩≡0 q 边界 源项为

q

c b

S ad =qL of /∆V abcd

L ef —实际边界与控制容积P 的两条边界相交部分的长度；

∆V abcd —控制容积P 的体积

扩充区域λ扩=0，则控容P 中的附加源项S ad 不会向扩充区域传递，从而实现了实附边界上的均匀热流加热条件。

（4）外部对流换热边界条件，l -t f

根据附加源项法，此时P 控制容积的两个附加源项为

S c , ad =

t f

1/α+δ/λ∆V

⋅

L of

S p , ad =-

L 1

⋅e

1/α+δ/λ∆V

δ—网格节点P 到实际边界的距离λ扩=0

区域扩充法的优点：可以用按规则区域编制的通用程序来计算非规则区域的问题；易于实施。缺点：浪费一些计算机内存及计算时间。 5

3、采用三角形网格 4 e f 外节点法对于不规则区域中的导热问题，采用三角形网格可以得

d 到比较满意的结果。

从不规则区域的三角形网格中划出围绕节点P 的多边形来分析。

3 c

确定节点P 所代表的控制容积：①三角形外心作为控容的顶点，要求：b 三角形为锐角三角形，以保证外心一定在三角形内；

2

②以三角形重心作为控容的顶点，对三角形形状无限制；外心为控制

6

1

容积顶点，P 控制容积如图所示，各部分平均导热系数分别为λa , λb , λc , λf ，用控制容积能量平衡法建立离散方程。

P-1之间的热导c p 1=λa /L a 1' /L p 1+λb L b 1' /L p 1

11

L p 1cot β1L b 1' =L p 1cot β2 22111

所以C p 1=λa cot β1+λb cot β2=(λa cot β1+λb cot β2)

222

1

常物性时，C p 1=λ(cotβ1+cot β2)

2

L a 1' =

其它各节点（2-6）与P 节点之间的热导亦按上式计算，而只需把β1, β2改换成相应顶角即可，P 控制容积的热容：

变物性时，图示各部阴影部分的热容量之和由控制容积P 的热平衡可得

∑(ρc ∆V

i

p , i i

)

∑(ρc ∆V )

i

t p -t o p ∆τ

=∑C pi (t i -t p ) +S (∆V p )

i

常物性时，上式变为

ρc ∆V p

t p -t o p ∆τ

=∑C pi (t i -t p ) +S (∆V p )

i

以上二式右端温度值的时层未标出，它取决于所采用的格式。 边界条件的处理：

第一类边界条件，无需专门处理；

第二、三类B 、C 边界上的节点，需要仿上述方法建立补充方程和图第二、三类边界上的节点P 的控制容积

b

2

1 β2 a

c

d

β4

3

P 1,t f

C p 1=C p 2

11

λb cot β1, C p 3=λc cot β3 221

=(λb cot β2+λc cot β4) 2

边界加入热量则计入到P 控制容积的源项中，对第三类边界条件，传入热量α(t f -t p ) L ad ，其中t p 的时层取决于所采用的格式。 Winslow 指出，当导热系数为常数时，三角形网格所建立的离散方程的系数矩阵具有正定性，可用SOR 求解；导热系数与温度有关时，为求迭代收敛，不宜采用SOR ，且两次迭代之间的C pi 可用欠松弛。

根据实际问题的需要，可以采用三角形网格与矩形网格的组合结构，例如二维纵向肋片管。

三角形网格的缺点：节点位置的确定、编号，节点间距的计算比较费时程序比较复杂。 内节点法：以每个三角形的外心为节点，其余同上。 4、组合坐标法

二维平直一弯曲段组成的通道，可采用坐标系组合的方法，例如直角坐标与极坐标的组合，把该计算系统分成三个区域，在区域I 、III 中采用直角坐标系，在区域II 中采用极坐标系。

三个区域中垂直于主流方向的坐标增量是一致的，dy =dr ，但在主流方向上的坐标不一致，为求得一致，在I 、III 区域中定义~x =x /R i ，则I 、

y

x

Ri θ r

x 也是弧度，III 区域中的~与II 区中一致，从而可以

x 坐标。 从入口截面开始就统一地连续计量~

在各区域中分别按相应的控制方程建立其离散

方程，然后在整个区域内统一地求解代数方程。

三个区域的两个分界线作为主控制容积的界面线，所以当采用交错网格时，主流方向速度的控制容积将跨越两个区域，如图所示。该速控制容积的a E , a w 分别按级坐标和直角坐标系计算，而a N , a S 可按并联处理：先分别按两区域中的半控制容积计算，后相加，源项亦仿此。

5、采用特殊的正交坐标系

三维正交曲线坐标系中，稳态对流一扩散问题的通用控制方程为

速度u 的控制容积 P u E

1∂1∂1∂

(ρu 1h 2h 3φ) +(ρu 2h 1h 3φ) +(ρu 3h 1h 2φ)

h 1h 2h 3∂x 1h 1h 2h 3∂x 2h 1h 2h 3∂x 3

=

h h ∂φh 1h 2∂φ1∂1∂

(Γ132) +(Γ) +S

h 1h 2h 3∂x 1h 2∂x h 1h 2h 3∂x 3h 3∂x 3

控制方程写成这种形式时，一切由曲线坐标，滞流模型引起的附加项都包括在S 中去了，因

此，S 的表达式一般比较复杂。

曲面边界平面通道中层流充分发展对流换热，设x 3主为主流方向，且为直线坐标，再利用充分发展条件，可得

h ∂φh ∂φ∂∂(Γ2) +(Γ1) =-Sh 1h 2=C ∂x 1h 1∂x 1∂x 2h 2∂x 2

式中，Γ, C 例于表中

在该流道截面的二维正交曲线坐标系中取出一个控制容积P 来分析，x 3方向控制容积厚度为

∆x 3，将上面的控制方程在该控制容积上作积分，等式左端第一项为

I =∆x 3⎰

n

s

h 2∂φ∂(Γ⎰ω∂x 1h 1∂x 1) dx 1dx 2

e

d

N

ω W

=∆x 3

⎰

n

s

[(Γ

h 2∂φh ∂φ) e -(Γ2) ω]dx 2 h 1∂x 1h 1∂x 1

n a

e P s

E b x 2

x 1

c

利用度规系数的性质，上式进一步变为

S

I ≅∆x 3[Γe

φE -φp (δs 1) PE

⎰

n

s

h 2

dx 2e

-Γω

φp -φωn

h 2dx 2] ⎰s (δs 1) wp ω

=∆x 3[Γe (φE -φP )

(∆s 2) e (∆s 1) e

-Γω(φp -φw ) ]

(δs 1) PE (δs 2) SP (∆s 1) n (∆s 1) s

-Γs (φp -φS ) ]

(δs 2) NP (δs 2) SP

类似地，左端第二项积分为：

II ≈∆x 3[Γn (φN -φp )

等号右端的积分为：C =-(S c +S p φp ) h 1h 2

III =∆x 3⎰

n

s

⎰ω

e

Cdx 1dx 2=-(S c +S p φp ⎰

n

s

⎰ωh h dx dx

12

1

e

2

=-(S c +S p φp ) ∆V p

其中控容体积∆V p 为：

⎡(∆s ) +(∆s 2) ω⎤⎡(∆s 1) n +(∆s 1) s ⎤

∆V p ≅⎢2e ∆x 3 ⎥⎢⎥22⎣⎦⎣⎦

将I 、II 、III 代入积分式，整理得到

a p φp =a E φE +a w φw +a N φN +a S φS +b

其中

a E =Γe (

(∆s 2) e (∆s 2) ω(∆s 1) n (∆s 1) s

), a w =, a N =Γn , a S =Γs

(δs 1) PE (δs 1) wp (δs 2) NP (δs 2) SP

a p =a E +a w +a N +a S -S p ∆V p , b =S C ∆V p

前面导出的三种典型坐标系中的二维导热的离散方程只不过是上式的特殊形式，从上式出发，可以编制正交坐标系中二维导热的通用程序，不同坐标系中只是规定h 1, h 2的不同计算式即

可。

（对流动问题，采用原始变量法时，亦需应用交错网格）

对于扩散（导热）问题，当计算区域的边界正好与正交曲线坐标系的坐标相吻合时，采用相当导热体的概念可以使这类问题的计算得到简化，无内热源时，能量方程为

或

h h h 3∂t h 1h 2h 3∂t h 1h 2h 3∂t ∂∂(λ12) +(λ)(λ) =0 222∂x 1∂x 2h 1∂x 1h 2∂x 3h 3∂x 3h 1h 2h 3∂t ∂

(λ) =0 ∑2∂x ∂x h i =1i i i

3

如果把x 1, x 2及x 3视为一新直角坐标系中的三个坐标，把λ

h 1h 2h 3

作为相当导热系数，记着2

h 1

λ*i ，则上式化为

∂*∂t

(λi ) =0 ∑∂x ∂x i =1i i

3

这样，在原直角坐标系中的一个曲面物体就被转换成了新直角坐标系中的规则长方体，导热系数由λ变成了各向异性的相当导热系数λ*原直角坐标系所代表的空间称为物理空间，而新直i ，角坐标系所代表的则称为计算空间，如图所示：

Z x 1 x 1 x

x 3

x 2

y

x 1

计算空间

x 2

物理空间

于是，该导热问题可先在计算空间中求解，由于它在计算空间中是规则区域，数值计算易于

进行，然后把计算空间中的计算结果传送到物理空间中去，两个空间中点的对应关系由曲线坐标的定义所规定。

在由物理空间向计算空间转换时，边界条件亦应作相应的转换而成为计算空间长方体的边界条件。

第一类边界条件：把给定温度赋给计算空间中相应的点即可。

x 2=const

把物理空间中的曲线坐标系当作计算空x

间中的直解坐标系，从而把物理空间中不规则y ,t f 形状的求解区域变换成计算空间中规则的长方x *

a 2=h 1h 3α 体（三维问题）或矩形（二维问题）

，这种做不

仅适用于导热问题，也可用于对流换热问题，

h h h *

2=123q n =h 1h 3q n 于是，有限差分法的一个主要弱点—不易处理

x 2 h 2不规则边界问题，可以说原则上已根本解决了。 x 1

现在的问题是，现成的曲线坐标系的数目

毕竟是有限的，而不规则形状的物体则是千变万化的，不可能满足上述变换的需要。要使上述变换思想付诸实现，必须发展一套方法，对于那些没有现成曲线坐标系可以利用的复杂形状，可以利用这套方法来建立一个与该物体相适应的坐标系，这就是适体坐标系的思想。

2

5-5 适体坐标系（Body-Fitted coordinates）

一、适体坐标系的基本概念

如上所述，最理想的坐标系是其各坐标轴与所计算物体的边界一一相符合的坐标系，称为适体坐标系，又称贴体坐标系，附体坐标系，无现成的坐标系可资利用时，希望通过计算来构造这样的坐标系，如图所示。

y

η η2 D

ξ1

A

η1

η

C

ξ2

ξ B

x

η1 η2 Δη A

ξ1

Δξ

B ξ ξ2 C

x-y 坐标系中的不规则区域ABCD ，

把其相交的两个边界作为曲线坐标系的两个轴，记为ξ, η，且：

（1）一条边界上，只能一个坐标单值地发生为化，另一个坐标保持为常数；

（2）在两条对应边界上，另一组曲线坐标的最大值与最小值对应相等，以便在计算平面上得出矩形区域。

问题：物理平面上，点（x ，y ）与计算平面上（ξ, η）的对应关系？求解结果的传送？

如果把ξ, η视为物理平面上的两个未知函数，那么上述确定ξ, η的问题就是物理平面上的一个边值问题，因此，从物理平面上来说，所谓要生成一个适体坐标系，实际上相当于要求解物理平面上的一个边值问题。

反之，首先把物理平面上的区域ABCD 按已规定的边界上的ξ, η值，画成计算平面上的ξ, η坐标系中对应的矩形，然后以均匀网格离散化，于是问题变为：已知计算区域边界上各节点（ξ, η）相应的（ξ, η），问在计算区域内部任意一点（ξ, η）对应的（ξ, η）值是多少？这样，如把x 、y 视为计算平面的未知函数，则生成适体坐标系的问题即为计算平面中的一个边值问题。

从数值计算的观点，对生成的适体坐标系有以下几点要求：

（1）物理平面上的点与计算平面上的点一一对应，同一簇中的曲线不能相变，不同簇中的两曲线只能相交一次；

（2）在适体坐标系中，每一个节点应当是一系列曲线坐标轴的交点，而不是一群三角形元素的顶点或无序的点群，以便设计有效，经济的算法及程序，采用矩形网格即可；

（3）物理平面求解区域内网格疏密程度易于控制；

（4）在物理平面的适体坐标的边界上，网格线最好与边界线正交或接近正交，以便于边界条件的离散化。

二、生成适体坐标系的方法分类 大致可分为三类 1、复变函数法

利用复变函数的映射，可以把相当一批二维不规则区域变换成矩形区域，而且可以得出解极或部分解析的变换关系式。

θ2π

r 12r

η1

ξ=(r -r 1) /(r 2-t 1)

η=0/2π

相应地：

x (ξη)

=[r 1+(r 2-r 1) ξ]cos(2πη)

y (ξ, η)

=[r 1+(r 2-r 1) ξ]xin (2πη)

η 2n 1n 2

复变函数映射

W =l n z =l n (re i θ) =l n r +i θ

ξ+i η ⎧ξ=l n r

⎩η=θ

所以⎨

∂2t ∂t 1∂2t

+2=0 控制方程的变化：2+2∂r r ∂θ∂r

因为η=θ，∴∂2t /∂θ2=∂2t /∂y 2

∂t ∂ξ∂2t

ξ=l nr ，∴∂t /∂r =+=0„仍然是各向同性！

∂ξ∂r ∂η2

代入原控制方程，得

∂2t ∂2t

（条件r ≠0）+2=0„仍然是各向同性。 2

∂ξ∂y

可见，正交坐标系法与复变函数法并不相同！

2、代数变换法—利有一些代数关系式来进行区域变换，具体实施方法颇多，最简单的是边界规范化的方法。

3、解微分方程的方法—通过求解边值问题的微分方程建立物理平面与计算平面上各点间的对应关系。

对该边值问题的控制方程的类型，物理问题本身无任何限定，可以比较自由地按照对所生成网格的要求来选择控制方程。

三、计算平面上求解区域形状的选择

1、物理平面上的区域由四条两两相交曲线构成的单连通域→计算平面上取为正方形或矩形； 2、物理平面的L 型区域：有两种选择

6

5 4

1

x

2

η

1

ξ 6

5 4

① 3 2

② 5 η

6 ξ

1

2

4

3

y 3

3、物理平面上的型区域：两种选择

特别注意：节点间一一对应关系不要受破坏，计算平面中求解节点方程时，对重迭线上的节点，每一个变量应有两套数组，以分别存放从左边和右边计算而得之值。

①

8 7 注意：计算平面上的尖角点未必对应于8 7

物理平面上的尖角点，反之亦然，例如物理

3 4 y 3 4 平面上一条连续封闭的曲线所围成的区域可η

以转换成计算平面上的一个正方形。 1 2 5 6 1 2 5 6

x ξ 选择计算平面上求解区域形状的一条

8 7 8 7 件：生成的网格要适用于计算的问题，通常，

9 在计算平面上都采用正方形均匀网格。 3 4

四、用适体坐标系求解流动与换热问题3 9 4 y η 的总体步骤

1 2 5 6 1 2 5 6

ξ x 1、网格生成：在计算平面上选择与物

理平面上复杂区域相应的求解区域，并找出两个区域内部节点之间的对应关系，

x =x (ξ, η), y =y (ξ, η) 或其反函数。这些关系可以是解析的，也可以是数值的。

2、控制方程的生成与离散，把物理平面上的控制方程及边界条件转换成计算平面上的相应形式，并利用控制容积积人法律立离散方程。

3、离散方程的求解及解的传递。在计算平面上获得收敛的解后，根据节点间的对应关系，传送到物理平面上去。

5-6 控制方程的转换及离散化

一、方程的转换

讨论如何把物理平面上的稳态对流一扩散方程的通用形式转换成一般曲线坐标系ξ-η中的形式，物理平面上的通用形式为

∂(ρu φ) ∂(ρu φ) ∂∂φ∂∂φ

+=(Γ) +(Γ) +S (x , y ) ∂x ∂y ∂x ∂x ∂y ∂y ∂φ∂φ∂ξ∂φ∂η=+

∂x ∂ξ∂x ∂η∂x

∂φ∂φ∂ξ∂φ∂η=+ ∂y ∂ξ∂y ∂η∂y

y =y (ξ, η) 要求找出反函数ξ, η的上述导数（用已知函数

设已有函数x =x (ξ, η)

x =x (ξ, η), y =y (ξ, η) 的导数来表示），这涉及到多元函数的微分关系：

∂ξ1∂y 1∂η-1∂y 1

==y η, ==-y ξ, ∂x J ∂ηJ ∂x J ∂ξJ ∂ξ1∂x 1∂η1∂x 1

=-=-x η, ==x ξ ∂y J ∂ηJ ∂y J ∂ξJ

式中，J 代有Jacobi 行列式：J =

∂x ∂y ∂x ∂y

-=x ξy η-x ηy ξ

∂ξ∂η∂y ∂ξ

∂φ1∂φ∂y ∂φ∂y 11∂(φy η) ∂(φy ξ) ∴=[-]=(φξy η-φηy ξ) =[-] ∂x J ∂ξ∂η∂η∂ξJ J ∂ξ∂η

=

1

[(φy η) ξ-(φy ξ) η] J

∂φ11

=(-φξx η+φηx ξ) =[-(-φx η) ξ+(φx ξ) η] ∂y J J

同理，对函数ρu φ, ρu φ可以写出：

∂(ρu φ) 1∂(ρu φy η) ∂(ρu φy ξ) =[-] ∂x J ∂ξ∂η∂(ρu φ) 1∂(ρu φx ξ) ∂(ρu φx η)

=[-] ∂y J ∂η∂ξ

扩散项的转换： 由：有

∂φ1∂φ∂φ=(y η-y ξ) ∂x J ∂ξ∂η

：

⎫∂∂φ∂11⎧∂Γ∂Γ(Γ) =[Γ(φξy η-φηy ξ)]=⎨[(φξy η-φηy ξ) y η]-[(φξy η-φηy ξ) y ξ]⎬∂x ∂x ∂x J J ⎩∂ξJ ∂ηJ ⎭

=

⎫1⎧∂Γ∂Γ22

[(φy -φy y )]-[(φx -φx x )]⎨ξηηξηηξξηξ⎬ J ⎩∂ξJ ∂ηJ ⎭

同理有：

⎫∂∂φ1⎧∂Γ∂Γ2(Γ) =-⎨[(φηx ξx η-φξx η)]-[(φηx ξ2-φξx ηx ξ)]⎬ ∂y ∂y J ⎩∂ξJ ∂ηJ ⎭

源项：S (x , y ) =S [x (ξ, η), y (ξ, η)]=S (ξ, η)

1∂(ρu φy η) ∂(ρu φy ξ) 1∂(ρu φx ξ) ∂(ρu φx η)

两对流项相加：得[-]+[-]

J ∂ξ∂ηJ ∂η∂ξ

=

1∂[ρφ(u η-ux η)]∂[ρφ(ux ξ-uy ξ)]1∂(ρu φ) ∂(ρv φ) {+=[+] J ∂ξ∂ηJ ∂ξ∂η

式中，U =uy η-ux η, V =ux ξ-uy ξ，可分别视为计算平面上ξ及η方向上的速度分量。 两扩散项相加：得

1∂Γ2∂Γ1∂Γ∂Γ2{[(x η+y η) φξ]-[(y ξy η+x ξx η) φη]}+[(x ξ2+y ξ2) φη]-[(y ξy η+x ξx η) φξ]}J ∂ξJ ∂ξJ J ∂ηJ ∂ηJ

=

1∂Γ1∂Γ

[(l φξ-βφη)]+[(-βφξ+x φη)] J ∂ξJ J ∂ηJ

2

2

式中α=x η+y η, β=x ξx η+y ξy η=计算平面上的守恒方程为：

∂x ∂x ∂y ∂y

+, r =x ξ2+y ξ2

∂ξ∂η∂ξ∂η

1∂1∂1∂Γ1Γ

(ρu φ) +(ρv φ) =[(αφξ-βφη)]+[(-βφξ+r φη)]+S (ξ, η) J ∂ξJ ∂ηJ ∂ξJ J J

可见，式中各项仍然保持了直角坐标中相应各项的意义，其中源项S 完全由直角坐标中的源项转换而来，其它各项在转换过程中并不给源项增添新的成分。

生成的曲线坐标在物理平面上是否正交？

β=0的点处，曲线坐标正交，且α和r 即为相应度规系数的平方；

β≠0的点处，曲线坐标斜交，β偏离零的多少反应了局部曲线坐标偏离正交的程度。

从物理平面到计算平面区域的简化（规则化）是以控制方程的复杂化为代价的。 二、方程的离散化