[未答复] 为什么温度场无变化呢?想让温度从最右侧传到左侧

[复制链接]
我要写出来 发表于 2021-8-22 11:21:55
  1. clear,clc;close all
  2. Hm=1600;
  3. Hn=1500;
  4. r_1=190/2e3;   
  5. r_2=200/2e3;        
  6. r_3=280/2e3;      
  7. k_b=1.5;           
  8. T_0=12;            
  9. T_f=10;                  
  10. alpha1=1.1e-6;      
  11. k_g=2.5;            
  12. beta1=1.2;           %等比系数
  13. dw=2.7;              
  14. dz=20;
  15. dt=90;
  16. dx=log(beta1);
  17. N=10;
  18. r=r_3*beta1.^(1:N);  
  19. M2=Hm/dz;
  20. Nj=Hn/dz;
  21. z=0:dz:Hm;
  22. N2=N+1;
  23. T=zeros(M2+1,N+1);
  24. %T表示p时刻的温度

  25. for i=1:M2+1
  26.     T(i,:)=T_0+dw/100*z(i);
  27. end

  28. for tt=1:400            
  29.     %A整个系数矩阵
  30.     A=zeros((M2+1)*N2,(M2+1)*N2);
  31.     b=zeros((M2+1)*N2,1);
  32.     %用k来记录方程总的数量
  33.     k=1;
  34.     %
  35.     for j=2:M2
  36.         for i=2:N2-1
  37.             Br=(alpha1*dt)/(r(i-1)*dx)^2;
  38.             Bz=(alpha1*dt)/(dz)^2;
  39.             A(k,[i-1+(j-1)*N2 i+(j-1)*N2 i+1+(j-1)*N2])=[-Br 2*(1+Br) -Br];
  40.             b(k)=Bz*T(j-1,i)+2*(1-Bz)*T(j,i)+Bz*T(j+1,i);
  41.             k=k+1;
  42.         end
  43.     end
  44.     %
  45.     A0=pi*r_3^2*(beta1-1)*(5-beta1)/3;
  46.     B1=(alpha1*2*pi*dt)/(A0*dx);
  47.     Bz1=(alpha1*dt)/(dz^2);
  48.     for j=2:Nj
  49.         A(k,[1+(j-1)*N2 2+(j-1)*N2])=[2+B1 -B1];
  50.         b(k)=Bz1*T(j-1,1)+2*(1-Bz1)*T(j,1)+Bz1*T(j+1,1);
  51.         k=k+1;
  52.     end
  53.     for j=Nj+1
  54.         A(k,[1+(j-1)*N2 2+(j-1)*N2])=[2+B1 -B1];
  55.         b(k)=Bz1*T(j-1,1)+2*(1-Bz1)*T(j,1)+Bz1*T(j+1,1);
  56.         k=k+1;
  57.     end
  58.     for j=Nj+2:M2
  59.          A(k,[1+(j-1)*N2 2+(j-1)*N2])=[2+B1 -B1];
  60.          b(k)=Bz1*T(j-1,1)+2*(1-Bz1)*T(j,1)+Bz1*T(j+1,1);
  61.          k=k+1;
  62.     end

  63.     %进口节点方程
  64.     for j=1
  65.         A(k,1)=1;
  66.         b(k)=T_f;
  67.         k=k+1;
  68.     end
  69.     %边界
  70.     for j=1:M2+1
  71.         A(k,j*N2)=1;
  72.         b(k)=T_0+dw/100*z(j);
  73.         k=k+1;
  74.     end
  75.     %边界
  76.     for i=1:N
  77.         A(k,i+1)=1;
  78.         b(k)=T_0;
  79.         k=k+1;
  80.     end
  81.     %边界
  82.     for i=1:N
  83.         A(k,i+M2*N2)=1;
  84.         b(k)=T_0+dw/100*z(end);
  85.         k=k+1;
  86.     end

  87.     T2=A\b;
  88.     T3=reshape(T2,N2,M2+1)';


  89.     %第二个交替方程
  90.     %B整个系数矩阵
  91.     B=zeros((M2+1)*N2,(M2+1)*N2);
  92.     b2=zeros((M2+1)*N2,1);
  93.     %用k来记录方程总的数量

  94.     %
  95.     k=1;
  96.     for j=2:M2
  97.         for i=2:N2-1
  98.             Br=(alpha1*dt)/(r(i-1)*dx)^2;
  99.             Bz=(alpha1*dt)/(dz)^2;
  100.             B(k,[i+(j-2)*N2 i+(j-1)*N2 i+j*N2])=[-Bz 2*(1+Bz) -Bz];
  101.             b2(k)=Br*T3(j,i-1)+2*(1-Br)*T3(j,i)+Br*T3(j,i+1);
  102.             k=k+1;
  103.         end
  104.     end
  105.    %
  106.     A0=pi*r_3^2*(beta1-1)*(5-beta1)/3;
  107.     B1=(2*pi*alpha1*dt)/(A0*dx);
  108.     Bz1=(alpha1*dt)/(dz^2);
  109.     for j=2:Nj
  110.         B(k,[1+(j-2)*N2 1+(j-1)*N2 1+j*N2])=[-Bz1 2*(1+Bz1) -Bz1];
  111.         b2(k)=(2-B1)*T3(j,1)+B1*T3(j,2);
  112.         k=k+1;
  113.     end
  114.     for j=Nj+1
  115.         B(k,[1+(j-2)*N2 1+(j-1)*N2 1+j*N2])=[-Bz1 2*(1+Bz1) -Bz1];
  116.         b2(k)=(2-B1)*T3(j,1)+B1*T3(j,2);
  117.         k=k+1;
  118.     end
  119.     for j=Nj+2:M2
  120.         B(k,[1+(j-2)*N2 1+(j-1)*N2 1+j*N2])=[-Bz1 2*(1+Bz1) -Bz1];
  121.         b2(k)=(2-B1)*T3(j,1)+B1*T3(j,2);
  122.         k=k+1;
  123.     end

  124.     %进口节点
  125.     for j=1
  126.         B(k,1)=1;
  127.         b2(k)=T_f;
  128.         k=k+1;
  129.     end
  130.     %边界
  131.     for j=1:M2+1
  132.         B(k,j*N2)=1;
  133.         b2(k)=T_0+dw/100*z(j);
  134.         k=k+1;
  135.     end
  136.     %边界
  137.     for i=1:N
  138.         B(k,i+1)=1;
  139.         b2(k)=T_0;
  140.         k=k+1;
  141.     end
  142.     %边界
  143.     for i=1:N
  144.         B(k,i+M2*N2)=1;
  145.         b2(k)=T_0+dw/100*z(end);
  146.         k=k+1;
  147.     end

  148.     T4=B\b2;
  149.     T5=reshape(T4,N2,M2+1)';
  150.     T=T5;
  151. end
复制代码



您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表