非结构化网格中热传导的数值计算

非结构化网格中热传导的数值计算

张艳丽1 张敏1 John C. Chai 2

1

南京理工大学动力工程学院，南京(210094) 2南洋理工大学机械与航天学院, 新加坡(639798)

摘 要：用基元有限容积方法求解热传导问题，结构化网格和非结构化网格同时采用，并和精确解相比较。虽然从两种网格形式，都能得到满意的结果，但结构化网格较三角形非结构化网格在规则区域更精确一些。为了提高非结构化网格的精度，二次扩散项是十分重要的。 关键词：结构化网格，非结构化网格，热传导

0. 引言

在过去三四十年中，对复杂扩散问题数值求解的研究有了长足的进步和发展。其中有限容积法和有限元法被研究者广泛地采用。起初，这两种方法主要应用在正交结构化网格，如笛卡尔直角坐标，柱坐标和极坐标系中[1]

。后来，它们被扩展应用到复杂几何形状的适体网格之中[2]

。现在，这些数值计算方法在非结构化网格中得到很好的应用

[3,4]

。

非结构化网格是近年来被广泛应用于数值计算的一种网格结构。它不但对几何边界具有很强的适应性，而且对局部网格的分解和组合，以及程序的扩展都具有较好的灵活性和简练性。本文的目的就是在结构化和非结构化网格中，用有限容积法，对标量扩散方程的数值解进行比较分析，以展现这些方法的优越之处。

1. 热控制方程和边界条件

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

0=+⎟⎟

⎠⎞⎜⎜⎝

⎛∂∂Γ∂∂

φφφS x x i i

(1.1)

其中φS 是单位体积中的净源项，φΓ是对应于变量φ 的扩散系数。在扩散问题中，一般情况下，会遇到三种边界条件。它们是：① 第一类边界条件(Dirichlet Problem)，给定变量φ ；② 第二类边界条件(Neumann Problem)，给定法向通量；③第三类边界条件(Robin Problem)，上两种边界条件的组合。它们的数学表达式分别为，

given

B φφ=，given B q q =，B

c B q q q φφ+=

(1.2)

p D ,

(a)

(b)

图1 (a)扩散项和(b)单位矢量

2. 非结构化网格中的离散方程

对方程(1.1)在非结构化网格中的离散可以写成，

0=+S D

(2.1)

在一个控制容积P 中，积分上式有，

1

=∆+∑=p p nb

i i V S D

(2.2)

其中，nb 为连接控制容积P 的面数，i D 是某面上总扩散项，并可以表示为，

i i i i A D →

→•∇Γ=φ

(2.3)

方程(2.3)中i Γ是交界面上的扩散系数。由此，i D 总扩散项可分成基本扩散项i p D ,和二次扩

散项i s D ,之和(见图1a ），

i s i p i D D D ,,+=

(2.4)

它们的物理意义分别是某一交界面上的法向扩散项和切向扩散项。当计算网格为正交时，二

次扩散项为零。基本扩散项和二次扩散项还可以写成，

i s i i

i i

P E i

i p e

A A A ds D ,,ˆ)(••−Γ=→→

→

φφ

(2.5a)

i t i s i s i i

i i

a b i

i s e e e

A A A A D ,,,,ˆˆˆ)(•••−Γ=→→

→

φφ

(2.5b)

这里下标P 代表计算单元的中心点，下标E 为相邻单元的中心点。图(2b)中给出了单位矢量

i s e ,ˆ和i t e ,ˆ。i ds 是P 点和E 点之间的距离。式(2.5a)中除i Γ之外，所有的项都能够计算出来。面积矢量→

i A ，单位矢量i s e

,ˆ和距离i ds 为几何量。相邻单元值E φ可以从上次迭代或初始猜测中获得。i P E ds )(φφ−是i s e

,ˆ方向上的分量。在图1b 平面i 中，式(2.5a)的另一种表达方式[5]有，

i s i ave i

P E e

ds ,,ˆ)()(•∇=−→

φφφ (2.6)

其中，i ave ,)(φ→

∇表示两个邻近单元的平均值。即P 和E 单元。利用式(2.6)和式(2.5a)可得出，

i s i i

i i s i ave i i p e

A A A e

D ,,,,ˆˆ)(•••∇Γ=→

→

→→

φ

(2.7)

二次扩散项或切向扩散项由式(2.5b)给出。 在计算它们时，我们需要知道垂线两端的变量值。由于计算这些量的困难性, 在求解二次扩散项时，我们可以选择另外一种方法。从方程(2.4)中，用总扩散项和基本扩散项的差，可得出二次扩散项，

i p i i s D D D ,,−=

(2.8)

代入总扩散项式(2.3)和基本扩散项式(2.7)，方程(2.8)变成，

i i s i i i i s i ave i i ave i

i i

s ds e A A A e A ds D ⎥⎥⎥

⎦⎤⎢⎢⎢⎣⎡

•••∇−•∇Γ=→→→→→→,,,,,ˆˆ)()(φφ

(2.9)

关于i i ds Γ和i ave ,)(φ→

∇的计算请看参考文献[5]。因此对于控制元P ，离散方程可写成，

0)(1

,1

=∆++−∑∑==P P nb

i i s P i nb i i V S D B φφ

(2.10)

最后简化式(2.10)得，

b a a i nb

i i P P +=

∑=φφ1

(2.11)

其中，

i s i i

i i

i i i e

A A A ds

B a ,ˆ••Γ=

=→→

→

(2.12a)

P P nb

i i s V S D b ∆+=

∑=1

,

(2.12b)

∑==

nb

i i P B a 1

(2.12c)

3. 算例计算和分析

3.1第一类边界条件（恒壁温）下的热传导

图2给出了这个问题的几何描述。顶墙为已知温度 T 2=50。其余的墙为相同的常温T 1=0。控制方程是：

0=+⎟⎟⎠

⎞

⎜⎜⎝⎛∂∂∂∂+⎟⎠⎞⎜⎝⎛∂∂∂∂S y T k y x T k x (3.1a) 这个问题的精确解(Kakac and Y ener,1993[6])是：

[][][]∑

∞

=−−=

−−=

1

121sinh sinh sin ])1(1[2),(),(n n n n n

H y x n T T T y x T y x λλλπ

θ (3.1b)

其中,

L

n n π

λ=

n = 1, 2, 3, ……，图3给出了三种网格，图4给出在400个控制单元下的温度分布。由于这个问题是最简单基本的，只给出它的数值解（图4）。本文的所有温度分布图中，实线代表精确解，虚线表示数值解。

图 2 问题的几何描述