传热学环肋肋效率matlab程序,利用matlab进行环肋肋效率数值计算

2016-zylzylzylzylzylzylzylzylzyl

(写在前面：Matlab相较于C来说界面更友好，操作起来也更方便，还可以直接做出好看的图，工科学生掌握一下可以避免跟C“苦苦纠缠”；迭代方法用的高斯赛德尔，其实还有其他一些数值计算方法，我还保存着一部分，会在后面文章中写一下；程序有注释，希望不要纯抄，看看注释很容易明白的；在编写过程中，我感觉节点微元怎么取很关键，取的不恰当结果会差很多，程序中的微元面积应该是扇形面积，请注意观察。)

正文：根据《传热学》中对导热方程的离散化，对离散方程进行进一步化简，可得：

(此处有图发不出来，其实是把离散方程同项合并一下，方便后面写成矩阵相乘的形式)

显然，可将离散方程表示为AΘ=b的形式，其中A为一个N×N的系数矩阵，Θ为一个1×N的解矩阵，b为一个1×N的右端项矩阵。这样做的目的是，利用矩阵的形式进行数值计算，较课本上的计算方法更有序，更容易从数值计算的角度理解。

以高斯-赛德尔迭代法为核心设计程序，程序设计框图如下：

(图发不出来)

高斯赛德尔迭代函数gaussseidel.m：

function[

x ] = gaussseidel( A,b,x0 )

format;%四位精度

%至少输入方程组的系数矩阵，方程组右端项，以及初始解向量。

D=diag(diag(A));

L=-tril(A,-1);

U=-triu(A,1);

G=(D-L)\U;

f=(D-L)\b;

x=G*x0+f;

n=1;

while(n<100001)%此循环决定迭代次数，此处统一迭代了100000次，方便满足取各节点值时的精度，也可根据需要修改循环，以达到某允许误差值。

x0=x;

x=G*x0+f;

n=n+1;

End

主程序heat.m:

clear

clc

format;

disp('环肋肋效率matlab程序.');

k=input('请输入k=r2/r1的值.');

k

N=input('请输入节点数N的值.');

N

m=input('请输入参数m的值.');

m

DeltaR = 1/(N-1)%定义步长

R=[]%定义无量纲半径矩阵

fori=1:1:N

R(i)=(1/(k-1))+((i-1)/(N-1));%为无量纲半径赋值

end

b=zeros(N,1);%定义一个N行1列的零矩阵，作为右端项。

b(1)=1%定义b的第1个值为1，以便表达出Θ1=1这个条件。

theta0=ones(N,1)%设置Θ的迭代初值。

A=zeros(N,N)%定义一个系数矩阵A。因为其中大部分元素为0，先使它成为一个零矩阵，再对特殊元素赋值。

A(1,1)=1;%定义A(1,1)的值为1，以便表达出Θ1=1这个条件。这个条件的表达必须将A(1，1)设为1，否则在调用迭代函数时会因有一项对角元为0而出错。

forj=2:1:N-1

A(j,j-1)=(1/(DeltaR^2))-(1/(2*DeltaR*R(j)));

A(j,j)=-(2/(DeltaR^2))-2*m^2;

A(j,j+1)=(1/(DeltaR^2))+(1/(2*DeltaR*R(j)));

end%此循环为系数矩阵A赋值。

A(N,N-1)=-1;

A(N,N)=1;%这两个赋值是为了表达出ΘN-ΘN-1=0,即ΘN=ΘN-1这个顶端绝热条件。

theta=gaussseidel(A,b,theta0)%调用gauss-seidel函数。

%以下开始计算肋效率

s1

= 0;

s2

= 0;

fori=2:1:N-1

s1

= s1+(R(i)-DeltaR/2)*DeltaR*theta(i);

s2

= s2+(R(i)-DeltaR/2)*DeltaR;

end

s1=s1+R(1)*DeltaR/2*theta(1)+(R(N)-DeltaR/2)*DeltaR/2*theta(N);

s2=s2+R(1)*DeltaR/2+(R(N)-DeltaR/2)*DeltaR/2;

Fi

= s1/s2;

disp('肋效率为：');

Fi

两个程序放在同一个文件夹内，运行主程序即可，分成两部分的原因在于迭代程序换成其他的迭代方法也未尝不可，故提供了再编辑的可能。

根据所设计的程序，可任意输入节点数，获得相应的肋效率。表一列出了r2/r1=2,m=2,节点数N取16、20、36、64、100、150、200、250、300时肋效率的值。

表一 节点数对肋效率的影响(r2/r1=2,m=2)

N

16

20

36

64

100

150

200

250

300

η

0.2773

0.2760

0.2741

0.2733

0.2730

0.2728

0.2728

0.2727

0.2727

由表一可知，在不同节点数下，肋效率可保证两位有效数字不再改变，特别是当节点数N取到150以上时，可使在三位有效数字下的解稳定。

根据网格无关性检验的结果，可取节点数N=150，分别计算当r2/r1=2、3、4，m取0.1、0.5、1.0、1.5、2.0、2.5时肋效率η的值，列于表二中。

表二 环肋肋效率随r2/r1及m的变化

mr2/r1

0.1

0.5

1

1.5

2

2.5

2

0.9907

0.8141

0.5428

0.3710

0.2728

0.2135

3

0.9887

0.7816

0.4909

0.3222

0.2307

0.1773

4

0.9871

0.7573

0.4553

0.2904

0.2040

0.1547

(1)表一中数据与《传热学》教材表4-2进行对比分析，两表的最终结果都在三位有效数字要求下稳定到0.272；将表二中数据与《传热学》教材表4-3进行对比分析，两表中各对应数据之差不超过0.001。

(2)根据表二作出η-m曲线图(见附录)，其中系列1为r2/r1=2的曲线,系列2为r2/r1=3的曲线，将其与《传热学》教材图2-20进行对比，两组曲线具有很高的拟合度。(图发不出来)实际上，《传热学》教材图2-20中环肋片效率曲线的绘制考虑了顶端换热的情况，会造成些许的误差，不过两者的高拟合度依然可以证明程序计算结果的可靠性。

(3)增大节点数N至200，重新计算表二中的数据，将计算结果列于表三中。与表二进行对比，各数据在三位有效数字下的解稳定，可认为获得了与网格无关的解，节点数对计算结果的影响可忽略不计。

表三 环肋肋效率随r2/r1及m的变化

mr2/r1

0.1

0.5

1

1.5

2

2.5

2

0.9908

0.8144

0.5426

0.3709

0.2728

0.2134

3

0.9888

0.7826

0.4907

0.322

0.2307

0.1773

4

0.9874

0.7591

0.4552

0.2903

0.2041

0.1547

首先对肋效率随m和r2/r1的变化进行定性分析：肋效率的物理意义为实际散热量与假设整个肋表面处于肋基温度下的散热量之比。由它的物理意义可以看出，肋片表面的平均温度越接近肋基温度，则肋效率就会越大。而要使肋片表面的平均温度接近肋基温度，则要求肋片导热系数尽量大以及表面传热系数h尽量小(参数m的值小)、肋片尺寸尽量小(r2/r1小)。由表二数据以及η-m曲线可以得出，环肋肋效率η随着m的增大而减小、随着r2/r1的增大而减小。这与定性分析的结果是相符的。