[数据结构与算法]板翅式换热器（锯齿形）matlab的

?试设计某实验平台所需要使用的一台气水换热器。

设计校核计算通过matlab实现，请求斧正！

clear,clc
%申明所有压力单位统一为Pa
%% 预备计算
%%预备计算
%%假设换热器效率n=0.8，任务书中要求的是大于等于0.7，x假设为为冷测水的进口流量
%初始值设定
t1_1=90;t2_1=25;n=0.8
x=890;
t1_11 = t1_1 - n*(t1_1 - t2_1);%%计算出口气温
%%由查表得出的空气的物性参数
cp1 = 1003+0.02*43+4*0.00001*43^2;
%%由查表得出的水的物性参数
cp2 = 4184.4-0.6964*21+1.06964*21+1.036*10^(-2)*21^2;
qm1=(1000+20*35)/3600;%36为学号
qm2=x/3600;
W1=qm1*cp1;
W2=qm2*cp2;
C=W1/W2;%计算C
%由冷测进口流量和C*的关系可判断，C*在x=1600以内是永远大于0.5的
t2_11=t2_1+W1/W2*(t1_1-t1_11);
%Wmax侧（冷侧）取算数平均温度
tm2=(t2_1+t2_11)/2
%两流体的平均温差为
deta_t1m=((t1_1-tm2)-(t1_11-tm2))/log((t1_1-tm2)/(t1_11-tm2));
%热侧平均温度为
tm1=tm2+deta_t1m;
%根据平均温度求物性参数
mu1=1.50619*10^(-6)*(tm1+273)^1.5;%空气的动力黏度miu
lam1=2.456*10^(-4)*(tm1+273)^0.823;%空气导热系数λ1
Pr1=mu1*cp1/lam1;%计算热侧的普朗特数Pr1
mu2=10^(230.298/(tm2+126.203)-4.5668);%水的动力黏度
lam2=0.598+1.373*10^(-3)*tm2-5.333*10^(-6)*(tm2^2);%水的导热系数
Pr2=mu2*cp2/lam2%计算冷侧的普朗特数Pr2
mixsum=[x;mu1;lam1;Pr1;mu2;lam2;Pr2]';
%% 产品结构规划及其计算
%参数说明     冷端          热端
%间断长度     ls1            ls2
%翅片间距     Pf1            Pf2
%板间距       s1             s2 
%翅片厚度     Deltaf1        Deltaf2
%翅片高度     h1             h2
%翅片层数     m1             m2
%翅片厚度           Deltap        
%侧板厚度           Deltas       
%封条宽度     bs1            bs2
%热侧进口宽度     l2
%冷侧出口         l1
%热侧出口         l3

    ls1=240       ;     ls2=60     ;   %间断长度 
     Pf1=50     ;      Pf2=40       %翅片间距
     s1= 150       ;     s2= 50    ;    %板间距 
     Deltaf1=5.5    ;    Deltaf2=3.35   ; % 翅片厚度     
     h1=75.5       ;     h2=65       ;   %翅片高度
    m1=370       ;     m2=400     ;      %翅片层数 
   Deltap=13   ;  %   Deltap2=       %翅片厚度  
    Deltas=40   ; %    Deltas2=     % 侧板厚度 
    bs1=150   ;         bs2= 90    ;     %封条宽度 
   l2=  1000       % 热侧进口宽度  
     l1=  1000       %  冷侧出口
      %l3=        %热侧出口 
l3 = m1*s1+m2*s2+2*m2*Deltap;    %计算非流动方向长度
de1=2*(Pf1-2*Deltaf1)*(s1-Deltaf1)/((Pf1-Deltaf1)+(s1-Deltaf1)+(s1-Deltaf1)*Deltaf1/ls1);
de2=2*(Pf2-2*Deltaf2)*(s2-Deltaf2)/((Pf2-Deltaf2)+(s2-Deltaf2)+(s2-Deltaf2)*Deltaf2/ls2);
%分别计算当量直径
phi1=(s1-Deltaf1)/(Pf1+2*Deltaf1+s1);
phi2=(s2-Deltaf2)/(Pf2+2*Deltaf2+s2);%计算翅片面积比
Ay1=(l2-2*bs1)*(l3-2*Deltas);
Ay2=(l1-1.5*bs2)*(l3-2*Deltas);%计算迎风面积
Vp1=(l1-1.5*bs2)*(l2-2*bs1)*m1*s1;
Vp2=(l1-1.5*bs2)*(l2-2*bs1)*m2*s2;%计算板间体积
%传热面积密度为
%以各板间体积为基准
A1=2*((Pf1-Deltaf1)+(s1-Deltaf1));
A2=2*((Pf2-Deltaf2)+(s2-Deltaf2));
beta1=A1/Vp1;
beta2=A2/Vp2;
%总传热面积为
A1=beta1*Vp1;%这里的单位是平方毫米
A2=beta2*Vp1;
%最小流通面积为
Ac1=de1*A1/(4*l1);
Ac2=de2*A2/(4*l2);
%每侧孔度为
sig1=Ac1/Ay1;
sig2=Ac2/Ay2;
Ap=(2*m1+2)*(l1-1.5*bs2)*(l2-2*bs1);%一次传热面积（隔板导热面积）(㎡)
%% 计算质量流速、雷诺数、对流表面传热系数、翅片效率和表面效率
%质量流速为
gm1=qm1/Ac1;
gm2=qm2/Ac2;
%计算雷诺数Re
Re1=gm1*de1/mu1;
Re2=gm2*de2/mu2;
%对于空气或气体工质，层流区Re<=1000时，适用Weiting拟合关系式（3-26）
%温流区Re>=2000时，适用Weiting拟合关系式（3-27）
%对于过渡区需要根据下式确定Re^*_j,
%%%%%%%%%%！！！！！！！！！！！！
%设置判断语句
%判断Re1,根据雷诺数爬满断面热侧空气的传热系数j1
if Re1<=1000
    j1=0.483*(ls1/de1)^0.162 * ((Pf1-Deltaf1)/(s1-Deltaf1))^(-0.184)*Re1^(-0.536);
else if Re1>=2000
        j1=0.242*(ls1/de1)^(-0.322)*(Deltaf1/de1)^0.089*Re1^(-0.368);
    else
        Re_j=61.9*((ls1/de1)^0.952)*((Pf1-Deltaf1)/(s1-Deltaf1))^(-1.1)*Re1^(-0.53);
        if Re1>=Re_j
        j1=0.242*(ls1/de1)^(-0.322)*((Pf1-Deltaf1)/(s1-Deltaf1))^(0.089)*Re1^(-0.368);
        else j1=0.483*(ls1/de1)^0.384 * ((Pf1-Deltaf1)/(s1-Deltaf1))^(-0.184)*Re1^(-0.368);
        end
    end
end
%Re_j=61.9*((ls1/de1)^0.952)*((pf1-Deltaf1)/(s1-Deltaf1))^(-1.1)*Re^(-0.53);
%根据题意，Re1>Re_j，因此热侧空气传热系数可用Weiting拟合公式（3-26）确定j
%否则用(3-27)确定
%j1=0.242*(ls1/de1)^(-0.322)*((pf1-Deltaf1)/(s1-Deltaf1))^(0.089)*Re^(-0.368);
%对于冷测的传热系数采用下式计算
Pr=2;%Pr暂定为2，需核实，！！！！！！！！！！！！！！！
j2=0.287*Re1^(-0.322)*Pr^(0.167);%%这里的Pr不知道取冷端还是热端。还是平均数；!!!!
%对流传热面系数
alp1=j1*gm1*cp1/Pr1^(2/3);
alp2=j2*gm2*cp2/Pr2^(2/3);
%翅片参数
m1=sqrt(2*alp1*(1+Deltaf1/ls1)/(lam1*Deltaf1));
m2=sqrt(2*alp2*(1+Deltaf2/ls2)/(lam2*Deltaf2));
nf1=tanh(m1*h1)/(m1*h1);
nf2=tanh(m2*h2)/(m2*h2);
%两侧翅片表面效率为
n01=1-phi1*(1-nf1);%%书上例题公式有Af1，不知其所指？？？？？已解决，使用等式后边的公式进行计算
n02=1-phi2*(1-nf2);
%% 计算壁面热阻和传热系数
lamw=169;
Rew=Deltap/(lamw*Ap);%这里的λw是指什么
%对于气-水换热器，暂不考虑污垢热阻，故有
KA1=1/(n01*alp1*A1)+Rew+1/(n02*alp2*A2);%KA1是指KA的倒数
KA=1/KA1;
%因为KA=K1A1=K2A2,当以热流体侧的总传热面积A1为基准时，
%对应的换热系数为
K1=KA/A1;

%% 计算传热单元数、换热器效率和出口流体温度
W1=qm1*cp1;
W2=qm2*cp2;
C=W1/W2;%C*=Wmin/Wmax;
NTU=KA/W1;%Wmin=W1;
%两流体各自非混流的叉流换热器效率，此处按德雷克近似关系式计算单个芯体的效率，有
ni=1-exp(((NTU^0.22)/C)*exp((-C*NTU^0.78)-1));
%换热器总效率为  n1
n1=(((1-C*ni)/(1-ni))^2-1/(((1-C*ni)/(1-ni))^2-C));%n1为设计换热效率。n为假设效率，第一次计算n=0。74
%此处需要判断是否符合当初假定的效率！！！！！！
%换热器的传热热流量为
Phi=n1*W1*(t1_1-t2_1);
%流体出口温度为
t1_11=t1_1-Phi/W1;
t2_11=t2_1-Phi/W2;
%与假设相对误差为
detan=abs((n1-n)/n);
%补充数组换热量与效率的关系。mixsum2=[];
%% 计算阻力
detap1=100;%预设的detap1
%比体积计算
p1_1=200*1000;
p2_1=100*1000;%任务书所给的初识进口压力；
p1_11=p1_1-detap1;
p2_11=p2_1-detap1;
R = 287;%是一个常数，定值
v1_1=R*t1_1/p1_1;
v1_11=R*t1_11/p1_11;
vm1=(v1_1+v1_11)/2;%计算热侧的参数，均为国际单位。
v2_1=R*t2_1/p2_1;
v2_11=R*t2_11/p2_11;
vm2=(v2_1+v2_11)/2;%计算冷测的参数；
%计算热侧阻力
%由于书中图3-12中查询K_1和K_11较繁琐，故将其曲线用线性关系近似代替
K1_1=(-3/7)*sig1+1.342857143;
K1_11=(-13/8)*sig1+0.825;
K2_1=(-3/7)*sig2+1.342857143;
K2_11=(-13/8)*sig2+0.825;
%设置判断语句
%判断Re1,根据雷诺数爬满断面热侧空气的f
if Re1<=1000
    f1=7.661*(ls1/de1)^0.384 * ((Pf1-Deltaf1)/(s1-Deltaf1))^(-0.092)*Re1^(-0.712);
else if Re1>=2000
        f1=1.136*(ls1/de1)^(-0.781)*(Deltaf1/de1)^0.534*Re1^(-0.198);
    else
        Re_f=41*((ls1/de1)^0.772)*((Pf1-Deltaf1)/(s1-Deltaf1))^(-0.179)*Re1^(-1.04);
        if Re1>=Re_j
        f1=7.661*(ls1/de1)^0.384 * ((Pf1-Deltaf1)/(s1-Deltaf1))^(-0.092)*Re1^(-0.712);
        else f1=1.136*(ls1/de1)^(-0.781)*(Deltaf1/de1)^0.534*Re1^(-0.198);
        end
    end
end
%此时已经求得f1；计算实际的detap1，实际的命名为detap1_1;
%A_1=1-sig1^2+K1_1;
%B_1=2*(v1_11/v1_1-1);
%C_1=(4*f1*2*l1*vm1)/(de1*v1_1);
%D_1=(1-sig1^2+K1_11)*(v1_11/v1_1);
%E_1=gm1^2*v1_1/2;
%detap1_1=E_1*(A_1+B_1+C_1-D_1);
detap1_1=((1-sig1^2+K1_1)+2*(v1_11/v1_1-1)+(4*f1*2*l1*vm1)/(de1*v1_1)-((1-sig1^2+K1_11)*v1_11/v1_1))*(gm1^2*v1_1/2);
%计算冷测阻力
%冷测为双流程，对于单流程预设detap20=80pa，采用相同的算法
%设置判断语句
%判断Re1,根据雷诺数爬满断面热侧空气的f
if Re2<=1000
    f2=7.661*(ls2/de2)^0.384 * ((Pf2-Deltaf2)/(s2-Deltaf2))^(-0.092)*Re2^(-0.712);
else if Re2>=2000
        f2=1.136*(ls2/de2)^(-0.781)*(Deltaf2/de2)^0.534*Re2^(-0.198);
    else
        Re_f=41*((ls2/de2)^0.772)*((Pf2-Deltaf2)/(s2-Deltaf2))^(-0.179)*Re2^(-1.04);
        if Re2>=Re_j
        f2=7.661*(ls2/de2)^0.384 * ((Pf2-Deltaf2)/(s2-Deltaf2))^(-0.092)*Re2^(-0.712);
        else f2=1.136*(ls2/de2)^(-0.781)*(Deltaf2/de2)^0.534*Re2^(-0.198);
        end
    end
end
%此时已经求得f2；计算实际的detap2，实际的命名为detap2_1;
%A_2=1-sig2^2+K2_1;
%B_2=2*(v2_11/v2_1-1);
%C_2=(4*f2*2*l2*vm1)/(de2*v2_1);
%D_2=(1-sig2^2+K2_11)*(v2_11/v2_1);
%E_2=gm1^2*v2_1/2;
%detap2_1=E_2*(A_2+B_2+C_2-D_2);
detap2_1=(((1-sig2^2+K2_1)+2*(v2_11/v2_1-1)+(4*f2*2*l2*vm2)/(de2*v2_1)-(((1-sig2^2+K2_11)*v2_11)/v2_1))*(gm2^2*v2_1/2));
detap2_2=detap2_1*2;%双流程压降为detap2_2
%
%% 强度校核
 %换热器热侧空气的入口温度为90°，由于铝材料随着温度升高强度指标有所下降，故选取100°时的强度进行强度校核。
 %防锈铝LF21-M的需应力
 %sig_Al = 75*1000*1000;
 sig_Al=23.75*1000*1000;
 %计算翅片厚度，校核结果的厚度用Deltaf1_1表示
 Deltaf1_1=((p1_1*2.5)/(sig_Al*1))+0.01;%这里乘1是指phi，phi是一个常数，取1.
 Deltaf2_1=((p2_1*2.5)/(sig_Al*1))+0.01;%算出来的东西单位是毫米mm。
 %计算隔板厚度，校核结果用Deltaf1_2表示
 Deltaf1_2=2.5*sqrt((6*p1_1)/(8*sig_Al))+0.05;
 Deltaf2_2=2*sqrt((6*p1_1)/(8*sig_Al))+0.05;