传热学上机C程序源答案之一维稳态导热的数值计算

一维稳态导热的数值计算

1.1物理问题

一个等截面直肋，处于温度t ∞=80

的流体中。肋表面与流体之间的对流换热系数为

h =45W/(m 2?℃),肋基处温度t w =300℃，肋端绝热。肋片由铝合金制成，其导热系数为λ=110W/(m ?℃)，肋片厚度为δ=0.01m ，高度为H=0.1m 。试计算肋内的温度分布及肋的总换热量。

1.2数学描述及其解析解

引入无量纲过余温度θ

=

t?t ∞t w ?t ∞

，则无量纲温度描述的肋片导热微分方程及其边界条件：

22

20d m dx

θθ-= x=0,θ=θw =1 x=H,

0x

θ?=?

其中m =

上述数学模型的解析解为：[()]

()()

w ch m x H t t t t ch mH ∞∞--=-?

()()w hp

t t th mH m

∞?=

-

1.3数值离散

1.3.1区域离散

计算区域总节点数取N 。

1.3.2微分方程的离散

对任一借点i 有：22

2

0i d m dx

θ

θ??-= ???

用θ在节点i 的二阶差分代替θ在节点i 的二阶导数，得：211

2

20i i i i m x θθθθ+--+-=

整理成迭代形式：()1122

1

2i i i m x θθθ+-=++ (i=2,3……,N-1)

1.3.3边界条件离散

补充方程为：11w θθ==

右边界为第二类边界条件，边界节点N 的向后差分得：1

0N N x

θθ--= ，将此式整理为

迭代形式，得：N 1N θθ-=

1.3.4最终离散格式

11w θθ== ()1122

1

2i i i m x θθθ+-=

++ (i=2,3……,N-1)

N 1N θθ-=

1.3.5代数方程组的求解及其程序

假定一个温度场的初始发布，给出各节点的温度初值：01θ，02θ，….，0

N θ。将这些初值代

入离散格式方程组进行迭代计算，直至收敛。假设第K 步迭代完成，则K+1次迭代计算式为：

K 11w θθ+=

()

11

11

2212i i K K K i m x

θθθ+-++=

++ (i=2,3……,N-1) 1

11N K K N θθ-++=

#include #include #define N 11 main() { inti;

float cha;/*cha 含义下面用到时会提到*/

float t[N],a[N],b[N];

float h,t1,t0,r,D,H,x,m,A,p; /*r 代表λ,x 代表Δx ，D 代表δ*/ printf("\t\t\t 一维稳态导热问题\t\t"); printf("\n\t\t\t\t\t\t----何鹏举\n");

printf("\n 题目：补充材料练习题一\n");

printf("已知：h=45，t1=80, t0=200, r=110, D=0.01, H=0.1 (ISO)\n"); /*下面根据题目赋值*/

h=45.0; t1=80.0; t0=300.0; r=110.0; D=0.01; H=0.1;

x=H/N; A=3.1415926*D*D/4; p=3.1415926*D; m=sqrt((h*p)/(r*A)); /*x 代表步长，p 代表周长，A 代表面积*/

printf("\n 请首先假定一个温度场的初始分布，即给出各节点的温度初值：\n");

{

scanf("%f",&t[i]);

a[i]=(t[i]-t1)/(t0-t1);

b[i]=a[i];/*这里b[i]用记录一下a[i]，后面迭代条件及二阶采用温度初场要用到*/ }

/*采用一阶精度的向后差分法数值离散*/

cha=1;

while(cha>0.0001)

{

a[0]=1;

for(i=1;i

a[i]=(a[i+1]+a[i-1])/(2+m*m*x*x);

a[N-1]=a[N-2];

cha=0;

for(i=0;i

cha=cha+a[i]-b[i];

cha=cha/N;/*cha代表每次迭代后与上次迭代各点温度差值的平均值*/

}

for(i=0;i

t[i]=a[i]*(t0-t1)+t1;

printf("\n\n经数值离散(一阶精度的向后差分法)计算得肋片的温度分布为:\n");

for(i=0;i

printf("%4.2f\t",t[i]);

printf("\n\n");

getchar();

/*采用二阶精度的元体平衡法数值离散(温度初值还用设定的初场，便于比较)*/ for(i=0;i

a[i]=b[i];

cha=1;

while(cha>0.0001)

{

a[0]=1;

for(i=1;i

a[i]=(a[i+1]+a[i-1])/(2+m*m*x*x);

a[N-1]=a[N-2]/(1+0.5*m*m*x*x);

cha=0;

for(i=0;i

cha=cha+a[i]-b[i];

cha=cha/N;

}

for(i=0;i

t[i]=a[i]*(t0-t1)+t1;

printf("\n\n经数值离散(二阶精度的元体平衡法)计算得肋片的温度分布为:\n");

printf("%4.2f\t",t[i]); printf("\n\n");

getchar();

}