连铸有限差分法传热数值模拟

杨铮...
已于 2024-03-07 10:22:27 修改
阅读量940
点赞数 4
分类专栏：建模文章标签：连铸计算传热学
于 2022-12-04 22:36:07 首次发布
本文链接：https://blog.csdn.net/weixin_37928884/article/details/128178574
版权
建模专栏收录该内容
42 篇文章
订阅专栏
在这里插入图片描述
有限差分法建立方坯连铸机传热方程，计算中对各冷却段采用不同的冷却方式。
参考：
大圆坯连铸传热数值模拟*
连铸坯凝固传热的数值模拟**
大方钢坯连铸二冷区传热与凝固过程数值模拟***
以连铸板坯为对象，横断面建立三维传热模型，计算分析连铸过程中的铸坯特征点温度、坯壳厚度，其结果对结晶器、二冷配水及轻压下区间优化提供指导，对改善板坯内裂纹与中心偏析提供帮助。
在这里插入图片描述
clc;clear;close all;
Tim_Mic=0.2;                                                       %离散时间间隔
X_Mic=0.004;Y_Mic=0.004;                                           %X和Y轴的离散间隔 
S_X=0.18/2;S_Y=0.18/2;                                             %铸坯尺寸180mm × 180mm 大方坯（结晶器断面尺寸）180指宽度b，180指高度h
num_x=round(S_X/X_Mic);                                          %X轴上的离散点个数
num_y=round(S_Y/Y_Mic);                                          %Y轴上的离散点个数 
%---------------------------------------工艺参数——————————————————————————————% 
Speed=1.3/60;                                                    %拉速1.3m/min 拉速m/s
Tem_C=1510;                                                      %浇注温度
Tem_W=25;                                                        %冷却水温度
Tem_S=20;                                                        %环境温度
Tem_M=200;                                                       %结晶器铜板温度
ETAM=1.1;                                                        %调整系数                                        
%---------------------------------------设备参数——————————————————————————————%
Lm=0.85;                                                         %结晶器有效长度（入口出口之间的距离）-12页
L_(1)=0.24;L_(2)=2.07;L_(3)=2.23;L_(4)=2.23;                     %各冷却区长度
L_left=16;                                                        %二冷段出口到矫直机的距离
time_m=round(Lm/(Speed*Tim_Mic));                                  %微元体从弯月面到结晶器出口经过的离散时间间隔数
time2=round((Lm+L_(1))/(Speed*Tim_Mic));                            %微元体从弯月面到足辊段出口经过的离散时间间隔数
time3=round((Lm+L_(1)+L_(2))/(Speed*Tim_Mic));                        %微元体从弯月面到一冷段出口经过的离散时间间隔数
time4=round((Lm+L_(1)+L_(2)+L_(3))/(Speed*Tim_Mic));                    %微元体从弯月面到二冷段出口经过的离散时间间隔数
time5=round((Lm+L_(1)+L_(2)+L_(3)+L_(4))/(Speed*Tim_Mic));                %微元体从弯月面到二冷段出口经过的离散时间间隔数
time_air=round((Lm+L_(1)+L_(2)+L_(3)+L_(4)+L_left)/(Speed*Tim_Mic));       %微元体从弯月面到拉娇机入口经过的离散时间间隔数
time_n(1)=time_m;
time_n(2)=time2;
time_n(3)=time3;
time_n(4)=time4;
time_n(5)=time5;
time_n(6)=time_air;
%---------------------------------------热物理参数—————————————————————————%
% 7.40 g/cm^3=7.40 t/m^3=7400 kg/m^3
S_Ds=7400;S_Dl=7000;S_Dsl=7200;                                   %钢坯的密度-24页
Ts=1378;Tl=1478;                                                  %钢坯的固液相线温度-24页
Lf=272000;                                                        %凝固潜热-20页
Cs=710;%钢坯在固相区比热-23页
Cl=810;%钢坯在液相区比热-23页 
Csl=(Cs+Cl)/2+Lf/(Tl-Ts);%两相区等效比热-23页
Ks=28.2;%钢坯在固相区导热系数-22页
Ksl=28.2;%两相区等效比热导热系数-22页
Kl=2*28.2;%钢坯在液相区导热系数-22页
%---------------------------------------计算初始化—————————————————————————%
Tem_Mic=ones(num_x,num_y,time_air);                                 %定义三维数组记录铸坯所有微元体在各时间点的温度
State_Mic=zeros(num_x,num_y);                                       %定义二维数组记录铸坯所有微元体的状态，0表示液态，1表示固液两相，2表示固相
Tem_Mic(:,:,1)=Tem_C*ones(num_x,num_y);                             %钢坯温度初始化
Tem_Center=zeros(time_air);Tem_Center(1)=Tem_C;                        %存储中心温度
Tem_Surface=zeros(time_air);Tem_Surface(1)=Tem_C;                      %存储表面温度
Thick=zeros(time_air);Thick(1)=0;                                      %存储坯壳厚度
%---------------------------------------计算热边界条件—————————————————————————%
h(1)=800;h(2)=800;h(3)=800; h(4)=800;                                     %二冷传热系数（换热系数）-14页
T_target(1)=1080;T_target(2)=1150;T_target(3)=1100;T_target(4)=1050;      %目标温度

%*******************************************计算结晶器传热******************************************
for i=2:time_m
    h_o=837.2-(0.0279*(Speed*100*(i-1)*Tim_Mic)-3.2209)*Speed*60*100;%经验公式
    h_m=h_o*ETAM;
    heat_transfer(i)=h_m;%传热系数（换热系数）
    for p=1:num_x
        for q=1:num_y
            if(Tem_Mic(p,q,i-1)>Tl)
                K(p,q,i-1)=Kl;C(p,q,i-1)=Cl;D(p,q,i-1)=S_Dl;
            else if(Tem_Mic(p,q,i-1)>Ts)
                K(p,q,i-1)=Ksl;C(p,q,i-1)=Csl;D(p,q,i-1)=S_Dsl;
                else
                K(p,q,i-1)=Ks;C(p,q,i-1)=Cs;D(p,q,i-1)=S_Ds;
                end
            end
        end
    end
     %中心点O点
     Tem_Mic(num_x,num_y,i)=Tem_Mic(num_x,num_y,i-1)+2*Tim_Mic/(D(num_x,num_y,i-1)*C(num_x,num_y,i-1)*X_Mic^2)*K(num_x,num_y,i-1)*(Tem_Mic(num_x-1,num_y,i-1)+Tem_Mic(num_x,num_y-1,i-1)-2*Tem_Mic(num_x,num_y,i-1)) ;
     Tem_Center(i)=Tem_Mic(num_x,num_y,i); 
     %边界点A点
     Tem_Mic(num_x,1,i)=Tem_Mic(num_x,1,i-1)+2*Tim_Mic/(D(num_x,1,i-1)*C(num_x,1,i-1)*X_Mic^2)*(K(num_x,1,i-1)*(Tem_Mic(num_x-1,1,i-1)+Tem_Mic(num_x,2,i-1)-2*Tem_Mic(num_x,1,i-1))+h_m*(Tem_M-Tem_Mic(num_x,1,i-1))*X_Mic);
     %边界点B点
     Tem_Mic(1,num_y,i)=Tem_Mic(1,num_y,i-1)+2*Tim_Mic/(D(1,num_y,i-1)*C(1,num_y,i-1)*X_Mic^2)*(K(1,num_y,i-1)*(Tem_Mic(2,num_y,i-1)+Tem_Mic(1,num_y-1,i-1)-2*Tem_Mic(1,num_y,i-1))+h_m*(Tem_M-Tem_Mic(1,num_y,i-1))*X_Mic);
     Tem_Surface(i)=Tem_Mic(1,num_y,i);
     %钢坯外部角部C点
     Tem_Mic(1,1,i)=Tem_Mic(1,1,i-1)+4*Tim_Mic/(D(1,1,i-1)*C(1,1,i-1)*X_Mic^2)*(1/2*K(1,1,i-1)*(Tem_Mic(2,1,i-1)+Tem_Mic(1,2,i-1)-2*Tem_Mic(1,1,i-1))+h_m*(Tem_M-Tem_Mic(1,1,i-1))*X_Mic);  
     for j=2:num_y-1
         %OA边
         Tem_Mic(num_x,j,i)=Tem_Mic(num_x,j,i-1)+Tim_Mic/(D(num_x,j,i-1)*C(num_x,j,i-1)*X_Mic^2)*K(num_x,j,i-1)*(2*Tem_Mic(num_x-1,j,i-1)+Tem_Mic(num_x,j+1,i-1)+Tem_Mic(num_x,j-1,i-1)-4*Tem_Mic(num_x,j,i-1));
         %钢坯外部CB边
         Tem_Mic(1,j,i)=Tem_Mic(1,j,i-1)+2*Tim_Mic/(D(1,j,i-1)*C(1,j,i-1)*X_Mic^2)*(1/2*K(1,j,i-1)*(2*Tem_Mic(2,j,i-1)+Tem_Mic(1,j+1,i-1)+Tem_Mic(1,j-1,i-1)-4*Tem_Mic(1,j,i-1))+h_m*(Tem_M-Tem_Mic(1,j,i-1))*X_Mic);
     end
     for j=2:num_x-1
         %OB边
         Tem_Mic(j,num_y,i)=Tem_Mic(j,num_y,i-1)+Tim_Mic/(D(j,num_y,i-1)*C(j,num_y,i-1)*X_Mic^2)*K(j,num_y,i-1)*(Tem_Mic(j+1,num_y,i-1)+Tem_Mic(j-1,num_y,i-1)+2*Tem_Mic(j,num_y-1,i-1)-4*Tem_Mic(j,num_y,i-1)) ;
         %钢坯外部CA边
         Tem_Mic(j,1,i)=Tem_Mic(j,1,i-1)+2*Tim_Mic/(D(j,1,i-1)*C(j,1,i-1)*X_Mic^2)*(1/2*K(j,1,i-1)*(Tem_Mic(j+1,1,i-1)+Tem_Mic(j-1,1,i-1)+2*Tem_Mic(j,2,i-1)-4*Tem_Mic(j,1,i-1))+h_m*(Tem_M-Tem_Mic(j,1,i-1))*X_Mic);
     end
     for j=2:num_x-1
         for k=2:num_y-1
         %内部节点
         Tem_Mic(j,k,i)=Tem_Mic(j,k,i-1)+Tim_Mic/(D(j,k,i-1)*C(j,k,i-1)*X_Mic^2)*(K(j,k,i-1)*(Tem_Mic(j+1,k,i-1)+Tem_Mic(j-1,k,i-1)+Tem_Mic(j,k+1,i-1)+Tem_Mic(j,k-1,i-1)-4*Tem_Mic(j,k,i-1)));
         end
     end
end
 %*******************************************计算二冷传热************************************************            
for n=1:4 %冷却区
    counter=0;
    while(abs(Tem_Mic(1,num_y,time_n(n+1))-T_target(n))>5 && counter<100)    
        for i=(time_n(n)+1):time_n(n+1)
            heat_transfer(i)=h(n);
    for p=1:num_x
        for q=1:num_y
            if(Tem_Mic(p,q,i-1)>Tl)
                K(p,q,i-1)=Kl;C(p,q,i-1)=Cl;D(p,q,i-1)=S_Dl;
            else if(Tem_Mic(p,q,i-1)>Ts)
                K(p,q,i-1)=Ksl;C(p,q,i-1)=Csl;D(p,q,i-1)=S_Dsl;
                else
                K(p,q,i-1)=Ks;C(p,q,i-1)=Cs;D(p,q,i-1)=S_Ds;
                end
            end
        end
    end
            %中心点O点
            Tem_Mic(num_x,num_y,i)=Tem_Mic(num_x,num_y,i-1)+2*Tim_Mic/(D(num_x,num_y,i-1)*C(num_x,num_y,i-1)*X_Mic^2)*K(num_x,num_y,i-1)*(Tem_Mic(num_x-1,num_y,i-1)+Tem_Mic(num_x,num_y-1,i-1)-2*Tem_Mic(num_x,num_y,i-1)) ;
            Tem_Center(i)=Tem_Mic(num_x,num_y,i); 
            %边界点A点
            Tem_Mic(num_x,1,i)=Tem_Mic(num_x,1,i-1)+2*Tim_Mic/(D(num_x,1,i-1)*C(num_x,1,i-1)*X_Mic^2)*(K(num_x,1,i-1)*(Tem_Mic(num_x-1,1,i-1)+Tem_Mic(num_x,2,i-1)-2*Tem_Mic(num_x,1,i-1))+h(n)*(Tem_W-Tem_Mic(num_x,1,i-1))*X_Mic);
            %边界点B点
            Tem_Mic(1,num_y,i)=Tem_Mic(1,num_y,i-1)+2*Tim_Mic/(D(1,num_y,i-1)*C(1,num_y,i-1)*X_Mic^2)*(K(1,num_y,i-1)*(Tem_Mic(2,num_y,i-1)+Tem_Mic(1,num_y-1,i-1)-2*Tem_Mic(1,num_y,i-1))+h(n)*(Tem_W-Tem_Mic(1,num_y,i-1))*X_Mic);
            Tem_Surface(i)=Tem_Mic(1,num_y,i);
            %钢坯外部角部C点
            Tem_Mic(1,1,i)=Tem_Mic(1,1,i-1)+4*Tim_Mic/(D(1,1,i-1)*C(1,1,i-1)*X_Mic^2)*(1/2*K(1,1,i-1)*(Tem_Mic(2,1,i-1)+Tem_Mic(1,2,i-1)-2*Tem_Mic(1,1,i-1))+h(n)*(Tem_W-Tem_Mic(1,1,i-1))*X_Mic);  
            for j=2:num_y-1
                %OA边
                Tem_Mic(num_x,j,i)=Tem_Mic(num_x,j,i-1)+Tim_Mic/(D(num_x,j,i-1)*C(num_x,j,i-1)*X_Mic^2)*K(num_x,j,i-1)*(2*Tem_Mic(num_x-1,j,i-1)+Tem_Mic(num_x,j+1,i-1)+Tem_Mic(num_x,j-1,i-1)-4*Tem_Mic(num_x,j,i-1));
                %钢坯外部CB边
                Tem_Mic(1,j,i)=Tem_Mic(1,j,i-1)+2*Tim_Mic/(D(1,j,i-1)*C(1,j,i-1)*X_Mic^2)*(1/2*K(1,j,i-1)*(2*Tem_Mic(2,j,i-1)+Tem_Mic(1,j+1,i-1)+Tem_Mic(1,j-1,i-1)-4*Tem_Mic(1,j,i-1))+h(n)*(Tem_W-Tem_Mic(1,j,i-1))*X_Mic);
            end
            for j=2:num_x-1
                %OB边
                Tem_Mic(j,num_y,i)=Tem_Mic(j,num_y,i-1)+Tim_Mic/(D(j,num_y,i-1)*C(j,num_y,i-1)*X_Mic^2)*K(j,num_y,i-1)*(Tem_Mic(j+1,num_y,i-1)+Tem_Mic(j-1,num_y,i-1)+2*Tem_Mic(j,num_y-1,i-1)-4*Tem_Mic(j,num_y,i-1)) ;
                %钢坯外部CA边
                Tem_Mic(j,1,i)=Tem_Mic(j,1,i-1)+2*Tim_Mic/(D(j,1,i-1)*C(j,1,i-1)*X_Mic^2)*(1/2*K(j,1,i-1)*(Tem_Mic(j+1,1,i-1)+Tem_Mic(j-1,1,i-1)+2*Tem_Mic(j,2,i-1)-4*Tem_Mic(j,1,i-1))+h(n)*(Tem_W-Tem_Mic(j,1,i-1))*X_Mic);
            end
            for j=2:num_x-1
                for k=2:num_y-1
                    %内部节点
                    Tem_Mic(j,k,i)=Tem_Mic(j,k,i-1)+Tim_Mic/(D(j,k,i-1)*C(j,k,i-1)*X_Mic^2)*(K(j,k,i-1)*(Tem_Mic(j+1,k,i-1)+Tem_Mic(j-1,k,i-1)+Tem_Mic(j,k+1,i-1)+Tem_Mic(j,k-1,i-1)-4*Tem_Mic(j,k,i-1)));
                end
            end
        end 
        if ((Tem_Mic(1,num_y,time_n(n+1))-T_target(n))>5)
            h(n)=h(n)+20;
        else if ((Tem_Mic(1,num_y,time_n(n+1))-T_target(n))<-5)
                 h(n)=h(n)-20;
            else
                 h(n)=h(n);
             end
        end 
        counter=counter+1;                      
    end %end while
end
 %*******************************************计算空冷传热******************************************
epsilon = 0.8;%黑度系数（发射率）
sigama = 5.67*10^(-8);%波兹曼常数
 for i=(time5+1):time_air
    Q_air=epsilon*sigama;                                           
    heat_transfer(i)=Q_air*((Tem_Mic(1,num_y,i-1)+273)^2+(Tem_S+273)^2)*((Tem_Mic(1,num_y,i-1)+273)+(Tem_S+273));%辐射传热系数（换热系数）
    for p=1:num_x
        for q=1:num_y
            if(Tem_Mic(p,q,i-1)>Tl)
                K(p,q,i-1)=Kl;C(p,q,i-1)=Cl;D(p,q,i-1)=S_Dl;
            else if(Tem_Mic(p,q,i-1)>Ts)
                K(p,q,i-1)=Ksl;C(p,q,i-1)=Csl;D(p,q,i-1)=S_Dsl;
                else
                K(p,q,i-1)=Ks;C(p,q,i-1)=Cs;D(p,q,i-1)=S_Ds;
                end
            end
        end
    end
    %中心点O点
    Tem_Mic(num_x,num_y,i)=Tem_Mic(num_x,num_y,i-1)+2*Tim_Mic/(D(num_x,num_y,i-1)*C(num_x,num_y,i-1)*X_Mic^2)*K(num_x,num_y,i-1)*(Tem_Mic(num_x-1,num_y,i-1)+Tem_Mic(num_x,num_y-1,i-1)-2*Tem_Mic(num_x,num_y,i-1)) ;
    Tem_Center(i)=Tem_Mic(num_x,num_y,i); 
    %边界点A点
    Tem_Mic(num_x,1,i)=Tem_Mic(num_x,1,i-1)+2*Tim_Mic/(D(num_x,1,i-1)*C(num_x,1,i-1)*X_Mic^2)*(K(num_x,1,i-1)*(Tem_Mic(num_x-1,1,i-1)+Tem_Mic(num_x,2,i-1)-2*Tem_Mic(num_x,1,i-1))+Q_air*((Tem_Mic(num_x,1,i-1)+273)^2+(Tem_S+273)^2)*((Tem_Mic(num_x,1,i-1)+273)+(Tem_S+273))*(Tem_S-Tem_Mic(num_x,1,i-1))*X_Mic);
    %边界点B点
    Tem_Mic(1,num_y,i)=Tem_Mic(1,num_y,i-1)+2*Tim_Mic/(D(1,num_y,i-1)*C(1,num_y,i-1)*X_Mic^2)*(K(1,num_y,i-1)*(Tem_Mic(2,num_y,i-1)+Tem_Mic(1,num_y-1,i-1)-2*Tem_Mic(1,num_y,i-1))+Q_air*((Tem_Mic(1,num_y,i-1)+273)^2+(Tem_S+273)^2)*((Tem_Mic(1,num_y,i-1)+273)+(Tem_S+273))*(Tem_S-Tem_Mic(1,num_y,i-1))*X_Mic);
    Tem_Surface(i)=Tem_Mic(1,num_y,i);
    %钢坯外部角部C点
    Tem_Mic(1,1,i)=Tem_Mic(1,1,i-1)+4*Tim_Mic/(D(1,1,i-1)*C(1,1,i-1)*X_Mic^2)*(1/2*K(1,1,i-1)*(Tem_Mic(2,1,i-1)+Tem_Mic(1,2,i-1)-2*Tem_Mic(1,1,i-1))+Q_air*((Tem_Mic(1,1,i-1)+273)^2+(Tem_S+273)^2)*((Tem_Mic(1,1,i-1)+273)+(Tem_S+273))*(Tem_S-Tem_Mic(1,1,i-1))*X_Mic); 
    for j=2:num_y-1
        %OA边
        Tem_Mic(num_x,j,i)=Tem_Mic(num_x,j,i-1)+Tim_Mic/(D(num_x,j,i-1)*C(num_x,j,i-1)*X_Mic^2)*K(num_x,j,i-1)*(2*Tem_Mic(num_x-1,j,i-1)+Tem_Mic(num_x,j+1,i-1)+Tem_Mic(num_x,j-1,i-1)-4*Tem_Mic(num_x,j,i-1));
        %钢坯外部CB边
        Tem_Mic(1,j,i)=Tem_Mic(1,j,i-1)+2*Tim_Mic/(D(1,j,i-1)*C(1,j,i-1)*X_Mic^2)*(1/2*K(1,j,i-1)*(2*Tem_Mic(2,j,i-1)+Tem_Mic(1,j+1,i-1)+Tem_Mic(1,j-1,i-1)-4*Tem_Mic(1,j,i-1))+Q_air*((Tem_Mic(1,j,i-1)+273)^2+(Tem_S+273)^2)*((Tem_Mic(1,j,i-1)+273)+(Tem_S+273))*(Tem_S-Tem_Mic(1,j,i-1))*X_Mic);
    end
    for j=2:num_x-1
        %OB边
        Tem_Mic(j,num_y,i)=Tem_Mic(j,num_y,i-1)+Tim_Mic/(D(j,num_y,i-1)*C(j,num_y,i-1)*X_Mic^2)*K(j,num_y,i-1)*(Tem_Mic(j+1,num_y,i-1)+Tem_Mic(j-1,num_y,i-1)+2*Tem_Mic(j,num_y-1,i-1)-4*Tem_Mic(j,num_y,i-1)) ;
        %钢坯外部CA边
        Tem_Mic(j,1,i)=Tem_Mic(j,1,i-1)+2*Tim_Mic/(D(j,1,i-1)*C(j,1,i-1)*X_Mic^2)*(1/2*K(j,1,i-1)*(Tem_Mic(j+1,1,i-1)+Tem_Mic(j-1,1,i-1)+2*Tem_Mic(j,2,i-1)-4*Tem_Mic(j,1,i-1))+Q_air*((Tem_Mic(j,1,i-1)+273)^2+(Tem_S+273)^2)*((Tem_Mic(j,1,i-1)+273)+(Tem_S+273))*(Tem_S-Tem_Mic(j,1,i-1))*X_Mic);
    end
    for j=2:num_x-1
        for k=2:num_y-1
            %内部节点
            Tem_Mic(j,k,i)=Tem_Mic(j,k,i-1)+Tim_Mic/(D(j,k,i-1)*C(j,k,i-1)*X_Mic^2)*(K(j,k,i-1)*(Tem_Mic(j+1,k,i-1)+Tem_Mic(j-1,k,i-1)+Tem_Mic(j,k+1,i-1)+Tem_Mic(j,k-1,i-1)-4*Tem_Mic(j,k,i-1)));
        end
    end
end
 %*******************************************计算结果******************************************
 
 %---------------------------------------表面中心角部点温度——————————————————————%
t=0:Tim_Mic*Speed:(time_air-1)*Tim_Mic*Speed;
for i=1:length(Tem_Mic(1,1,:))
    Tcorner(i) = Tem_Mic(1,1,i);
end

figure(10);
str0 = '中心温度';
str1 = '表面温度';
str2 = '角部温度';

plot(t,Tem_Surface,'r');axis([0 18 600 1600]);text(9,Tem_Surface(1100),str0)
hold on
plot(t,Tem_Center,'g');axis([0 18 600 1600]);text(9,Tem_Center(1100),str1)
hold on
plot(t,Tcorner,'b');axis([0 18 600 1600]);text(9,Tcorner(1100),str2)
xlabel('距离弯月面 m');ylabel('温度 ℃');
title('铸坯表面温度、中心温度、角部温度图')
 %---------------------------------------固相分数—————————————————————————%
for i=1:time_air
    for j=1:num_y
        if (Tem_Mic(num_x,j,i)>Ts&&Tem_Mic(num_x,j,i)<Tl)
            Solid_frac(i,j)=1-(Tem_Mic(num_x,j,i)-Ts)/(Tl-Ts);
        else if (Tem_Mic(num_x,j,i)>=Tl)
            Solid_frac(i,j)=0;
            else
             Solid_frac(i,j)=1;   
            end
        end
    end
end
x=0:Tim_Mic*Speed:(time_air-1)*Tim_Mic*Speed;
y=0:Y_Mic:(num_y-1)*Y_Mic;
%figure(2);pcolor(x,y,Solid_frac');shading interp;xlabel('距离弯月面');ylabel('固相分率');colorbar;
 %---------------------------------------坯壳厚度—————————————————————————%
for i=1:time_air
for j=1:num_y
    if (Tem_Mic(num_x,j,i)>Ts)
        Thick(i)=(j-1)*X_Mic;break;
    else 
            Thick(i)=j*X_Mic;
        
    end
end
%Thick(i)=(num_y-1)*X_Mic;
end
% figure(3);plot(t,Thick);xlabel('距离弯月面');ylabel('坯壳厚度 mm'); %凝固坯壳

figure(11);
plot(t,heat_transfer)
xlabel('距离弯月面 m');ylabel('传热系数 W/（m^2*K）');