计算流体力学一维稳态导热编程作业

The Finite Volume Method for One-Dimensional
Diffusion Problems
I. Problem of 1-dimensional Steady-State Source-free Heat Conduction
Consider the problem of source-free heat conduction in an insulated rod whose ends are maintained at constant temperatures of 100℃ and 500℃ respectively. The one-dimensional problem sketched in the Figure 1 is governed by
0=⎪⎭
⎫ ⎝⎛dx dT k dx d Calculate the steady state temperature distribution in the rod. Thermal conductivity k equals 1000W/m/K, cross-sectional area A is 10-2m 2
.
Fig. 1 Physical Model
1网格划分
条件
2
210500100,//1000,
1.05/,5.0m A T T K m W k m L x m L B A -======∆=℃℃，
2方程离散
0.5 m
T B =500
T A =100
Area(A) B
A
求解域内共有5个节点，节点2、3、4的离散方程：
W w WP w E e PE e W P w WP w e PE
T A x k T A x k T T A x k A x k ⋅⎪⎪⎭⎫ ⎝⎛+⋅⎪⎪⎭⎫
⎝⎛+⋅=⋅⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛δδδδ0e 由于k k k w e ==，x x x WP Pe δδδ==,A A A W e ==均为常数，因此对节点2、3、4有离散方程：
E E W W P P T a T a T a +=
式中A x k a W δ=
,A x
k a E δ=，E W P a a a +=, 节点1的离散方程：0=---AP
A
P w w PE P E e
e x T T A k x T T A k δδ
A w AP w E e PE e W P w AP w e PE
e T A x k T A x k T T A x k A x k ⋅⎪⎪⎭⎫
⎝⎛+⋅⎪⎪⎭⎫ ⎝⎛+⋅=⋅⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛δδδδ0 可写为：u E E W W P P S T a T a T a ++= 其中e PE
e
E A x k a δ=
, 0=W a , P E W P S a a a -+=, w AP
w P A x k
S δ2-=, A w AP w u
T A x k S ⎪⎪⎭
⎫ ⎝⎛=δ2 同理，节点5的离散方程0=---WP
W P w w PB P
B e
e x T T A k x T T A k δδ
B e PB e W w WP w E P e PB e w WP
w T A x k T A x k T T A x k A x k ⋅⎪⎪⎭⎫
⎝⎛+⋅⎪⎪⎭⎫ ⎝⎛+⋅=⋅⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛δδδδ0 可写为：u E E W W P P S T a T a T a ++= 其中0=E a , w WP
w
W A x k a δ=
, P E W P S a a a -+=, PB
e
P x k S δ2-
=, B e PB e u T A x k S ⋅⎪⎪⎭
⎫ ⎝⎛=δ2 所以得到各节点的a W , a E , a P , S c , S p 的值
各节点离散方程系数
值 节点
a W a E
S p
S u
a P= a W +a E -S p
1 0
A x k δ A x
k δ2
- A AT x
k
δ2
A x k δ3
2 A x k δ A x k δ 0 0 A x k δ2
3 A x
k δ A x
k δ 0 0 A x k δ2 4 A x
k δ A x k δ 0
A x k δ2 5
A x
k δ 0
A x
k δ2
- B AT x
k
δ2
A x k δ3 即：
值 节点
a W a E S p S u
a P= a W +a E -S p
1 0 100 -200 A T 200
300 2 100 100 0 0 200 3 100 100 0 0 200 4 100 100 0 0
200 5
100
-200
B T 200
300
从而得下述代数方程组 ⎪⎪⎪⎩⎪
⎪⎪⎨⎧+=+=+=+=+=B A T T T T
T T T T T T
T T T T T 200100300100100200100100200100100200200100300455
3442331221 写成矩阵形式有：
⎥⎥⎥⎥⎥⎥⎦⎤⎢
⎢
⎢
⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--------B A T T T T T T T 200000200300100000100200100000100200100000100200100000100
30054321
将500,100
==B A T T 代入，解得此方程组为：
⎥
⎥
⎥
⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡46038030022014054321T T T T T 3程序
C 程序内容如下：
#include "stdafx.h" #include <stdio.h> #define N 5 int main() {
int i;
double GL,TA,TB,dx,k,Area,Lenth;
double AW[N+1],AE[N+1],AP[N+1],SP[N+1],SU[N+1]; double
a[N+1],b[N+1],c[N+1],f[N+1],M[N+1],L[N+1],U[N+1],Y[N+1],T[N+1];
Lenth=0.5; //长度
TA=100; //西边界温度 TB=500; //东边界温度 dx=Lenth/N; //网格宽度 k=1000; //导热系数 Area=0.01; //截面积
//第一点计算 AW[1]=0.0;
AE[1]=k*Area/dx; SP[1]=-2*k*Area/dx; SU[1]=2*k*Area/dx*TA; // 最后一点计算 AW[N]=k*Area/dx; AE[N]=0.0;
SP[N]=-2*k*Area/dx; SU[N]=2*k*Area/dx*TB; // 中间点计算
for(i=2;i<=N-1;i++) {
AW[i]=k*Area/dx; AE[i]=k*Area/dx; SP[i]=0.0; SU[i]=0.0;
}
//计算Ap
for(i=1;i<=N;i++) AP[i]=AW[i]+AE[i]-SP[i];
//三对角方程组追赶法求解
for(i=1;i<=N;i++)
{
a[i+1]=-AE[i];
b[i]=AP[i];
c[i]=-AW[i+1];
f[i]=SU[i];
}
for(i=1;i<=N-1;i++)
U[i]=c[i];
L[1]=b[1];
for(i=2;i<=N;i++)
{M[i]=a[i]/L[i-1];
L[i]=b[i]-M[i]*U[i-1];}
Y[1]=f[1];
for(i=2;i<=N;i++)
Y[i]=f[i]-M[i]*Y[i-1];
T[N]=Y[N]/L[N];
for(i=N-1;i>=1;i--)
T[i]=(Y[i]-U[i]*T[i+1])/L[i];
printf("方程组ax=b的解为:\n");
for(i=1;i<=N;i++)
printf(" T[%d]=%.3f\n",i,T[i]);
return 0;
}。