主程序：

while(1)
L=input('请输入网格单元的长度L（mm）(L<10mm)：');
if  (mod(10,L)==0)
    break;
end
end
M=10/L+1;
Temp=zeros(M,1);
Th=input('请输入热边温度（K）(Th>300K)：');
Tc=300;
P=input('请输入环境压力（Pa）：');
Temp(1,1)=Th;%热壁面温度赋值
Temp(M,1)=Tc;%冷壁面温度赋值
t=zeros(M,1);%验证矩阵t

Coe=zeros(M,1);
Coe(1,1)=Ouhe(Th,P);%热壁面温度对应k值求解
Coe(M,1)=Ouhe(Tc,P);%冷壁面温度对应k值求解
c=zeros(M,1);%验证矩阵c

while(1)
    T=Temp;%临时存放数值矩阵T
    C=Coe;%临时存放数值矩阵C
for l=2:M-1
   
Temp(l,1)=(Coe(l+1,1)*(Coe(l-1,1)+Coe(l,1))*Temp(l+1,1)+Coe(l-1,1)*(Coe(l+1,1)+Coe(l,1))*Temp(l-1,1))/(Coe(l,1)*(Coe(l-1,1)+Coe(l+1,1))+2*Coe(l-1,1)*Coe(l+1,1));
Coe(l,1)=Ouhe(Temp(l,1),P);
        
end
for l=1:M%判断结果是否已经收敛
if (abs(Temp(l,1)-T(l,1))<1e-1)&&(abs(Coe(l,1)-C(l,1))<1e-4)
            t(l,1)=1;
            c(l,1)=1;%判断赋值矩阵每一项是否都为1，即每项收敛
end
end
if  (t==ones(M,1))&(c==ones(M,1))
   break;
end
end
k = harmmean(C);
sprintf('该温度下的有效导热系数为 %f w/m.C',k)

调用函数：

function K=Ouhe(T,P)
syms f r Cp Pr B  %定义变量：固体比 比热率 比热 普朗特数 消光系数
syms Ks Kr Kg K %定义变量：固体导热 辐射导热 气体导热 耦合导热
syms Lc s p %定义变量：特征长度 分子平均自由程 密度
p=24;%p先假设密度24kg/m3
s=1.38*10^(-23)*T/(sqrt(2)*pi*(3.798*10^(-10))^2*P);
Lc=pi*2.4*10^(-6)/(4*p/2600);
B=5728.624+0.188*T;
r=1040.05-0.0116361*T+4.639*10^(-4)*T^2-2.82*10^(-7)*T^3;
Pr=0.734577-1.4446*10^(-4)*T+1.16345*10^(-6)*T^2+3.207*10^(-10)*T^3;
Kr=16*(5.67*10^(-8))*T^3/(3*B);
Ks=(p/2600)^2*(0.653+1.49*10^(-3)*T);
Kg=(0.0241+7.407*10^(-5)*T+2.462*10^(-8)*T^2)/(1+4*r*s/(Pr*Lc*(r+1)));
K=0.8*(p/2600*(Ks)+(1-p/2600)*Kg)+0.2*Ks*Kg/((1-p/2600)*Ks+p/2600*Kg)+Kr;
end

我把步长设置为1即L为10是可以出结果的，但是L为1时就计算不出，求解决办法 返回小木虫查看更多