计算流体—流过平板的超声速流动

kerbalXspace

已于 2024-04-23 10:53:04 修改

阅读量443

点赞数 5

文章标签：学习

于 2024-04-23 10:52:40 首次发布

本文链接：https://blog.csdn.net/kerbalXspace/article/details/138115714

版权

计算流体力学基础及其应用

有没有人能看出这是哪的错误？我算的结果第二步温度就成复数了，一直看不出来那有错

感觉出问题的就是结果对数据太敏感，只有变了一点点，就直接偏离结果。

欢迎交流。

前三个是算一些参数的函数

function k=THERMC(Pr,miu,Cp)
    k=miu*Cp/Pr;
return

function tao=TAUXY(u,v,T,miu,k,i,j,n)
    delta_x=1.449e-07;
    delta_y=1.19e-07;
    lamuda=-2/3;
    if i==1
        partial_u_y=(u(i+1,j)-u(i,j))/delta_y;
        partial_v_y=(v(i+1,j)-v(i,j))/delta_y;
        partial_T_y=(T(i+1,j)-T(i,j))/delta_y;
    elseif i==71
        partial_u_y=(u(i,j)-u(i-1,j))/delta_y;
        partial_v_y=(v(i,j)-v(i-1,j))/delta_y;
        partial_T_y=(T(i,j)-T(i-1,j))/delta_y;
    else
        partial_u_y=(u(i+1,j)-2*u(i,j)+u(i-1,j))/2/delta_y;
        partial_v_y=(v(i+1,j)-2*v(i,j)+v(i-1,j))/2/delta_y;
        partial_T_y=(T(i+1,j)-2*T(i,j)+T(i-1,j))/2/delta_y;
    end
    if j==1
        partial_u_x=(u(i,j+1)-u(i,j))/delta_x;
        partial_v_x=(v(i,j+1)-v(i,j))/delta_x;
        partial_T_x=(T(i,j+1)-T(i,j))/delta_x;
    else
        partial_u_x=(u(i,j)-u(i,j-1))/delta_x;
        partial_v_x=(v(i,j)-v(i,j-1))/delta_x;
        partial_T_x=(T(i,j)-T(i,j-1))/delta_x;
    end
        switch n
            case 1
                tao=lamuda*(partial_u_x+partial_v_y)+2*miu(i,j)*partial_u_x;
            case 2
                tao=miu(i,j)*(partial_v_x+partial_u_y);
            case 3
                tao=lamuda*(partial_u_x+partial_v_y)+2*miu(i,j)*partial_v_y;
            case 4
                tao=-k(i,j)*partial_T_x;
            case 5
                tao=-k(i,j)*partial_T_y;  
        end
return

function miu=DYNVIS(miu0,T0,T)
    miu=miu0*((T0+110)/(T+110))*(T/T0)^(3/2);
return

这里才是主要程序，代码水平很差，没怎么学过，只会用matlab。

clc
clear all
%%初始条件
i_max=71;
j_max=71;
maxit=zeros(i_max,j_max);
Ma=4;
T=288.16;%来流温度
a=340.28;
LHORI=0.00001;
P=101325;
Tw=T;%壁面温度
gamma=1.4;
Pr=0.71;
miu0=1.789e-5;
R=287;
T0=288.16;%初始温度
miu=DYNVIS(miu0,T0,T);
Cv=R/(gamma-1);
Cp=gamma*Cv;
rho=P/(R*T);
ReL=rho*Ma*a*LHORI/miu;
e=Cv*T;
k=THERMC(Pr,miu,Cp);
delta=5*LHORI/sqrt(ReL);
LVERT=5*delta;
delta_x=LHORI/(j_max-1);
delta_y=LVERT/(i_max-1);
K=0.6;
for time =1:100
    
    if time==1
        u=ones(size(maxit))*Ma*a;
        v=zeros(size(maxit));
        rho=ones(size(maxit))*rho;
        p=ones(size(maxit))*P;
        e=ones(size(maxit))*e;
        T=ones(size(maxit))*T;
        miu=ones(size(maxit))*miu;
        k=ones(size(maxit))*k;
        a=ones(size(maxit))*a;
    end
    %% TSTEP
    vv=4*gamma*miu./(3*rho*Pr);
    delta_t=1./(abs(u)/delta_x+abs(v)/delta_y+a*sqrt(1/delta_x^2+1/delta_y^2)+2*vv*(1/delta_x^2+1/delta_y^2));
    delta_tmin=min(min(K*delta_t));
    
    %% MAC
    U1=rho;
    U2=rho.*u;
    U3=rho.*v;
    U4=rho.*(e+(u.^2+v.^2)/2);
    for i=1:71
        for j=1:71
            E1(i,j)=rho(i,j)*u(i,j);
            E2(i,j)=rho(i,j)*u(i,j)*u(i,j)+p(i,j)-TAUXY(u,v,T,miu,k,i,j,1);
            E3(i,j)=rho(i,j)*u(i,j)*v(i,j)-TAUXY(u,v,T,miu,k,i,j,2);
            E4(i,j)=(U4(i,j)+p(i,j))*u(i,j)-u(i,j)*TAUXY(u,v,T,miu,k,i,j,1)-v(i,j)*TAUXY(u,v,T,miu,k,i,j,2)+TAUXY(u,v,T,miu,k,i,j,4);
            F1(i,j)=rho(i,j)*v(i,j);
            F2(i,j)=rho(i,j)*u(i,j)*v(i,j)-TAUXY(u,v,T,miu,k,i,j,2);
            F3(i,j)=rho(i,j)*v(i,j)*v(i,j)+p(i,j)-TAUXY(u,v,T,miu,k,i,j,2);
            F4(i,j)=(U4(i,j)+p(i,j))*v(i,j)-u(i,j)*TAUXY(u,v,T,miu,k,i,j,2)-v(i,j)*TAUXY(u,v,T,miu,k,i,j,3)+TAUXY(u,v,T,miu,k,i,j,5);
        end
    end
    U_conbined=cat(3,U1,U2,U3,U4);
    E_conbined=cat(3,E1,E2,E3,E4);
    F_conbined=cat(3,F1,F2,F3,F4);  
    partial_U_t=zeros(size(U_conbined));
    partial_U_t(2:70,2:70,1:4)=-(E_conbined(2:70,3:71,1:4)-E_conbined(2:70,2:70,1:4))/delta_x-(F_conbined(3:71,2:70,1:4)-F_conbined(2:70,2:70,1:4))/delta_y;
    U_m=U_conbined+delta_tmin*partial_U_t;
    for i=2:70
        for j=2:70
         rho_m(i,j)=U_m(i,j,1);
         u_m(i,j)=U_m(i,j,2)/U_m(i,j,1);
         v_m(i,j)=U_m(i,j,3)/U_m(i,j,1);
         e_m(i,j)=U_m(i,j,4)/U_m(i,j,1)-(u_m(i,j)^2+v_m(i,j)^2)/2;
         T_m(i,j)=e_m(i,j)/Cv;
         p_m(i,j)=rho_m(i,j)*R*T_m(i,j);
        end
    end 
   
    %边界条件
    %出口边界条件，不包括上下边界的交点
    p_m(2:70,71)=2*p_m(2:70,70)-p_m(2:70,69);
    u_m(2:70,71)=2*u_m(2:70,70)-u_m(2:70,69);    
    v_m(2:70,71)=2*v_m(2:70,70)-v_m(2:70,69);   
    T_m(2:70,71)=2*T_m(2:70,70)-T_m(2:70,69);   
    %物面边界条件
    u_m(1,:)=0;
    v_m(1,:)=0;
    p_m(1,:)=2*p_m(2,:)-p_m(3,:);
    T_m(1,:)=288.16;
    %进口边界条件
    u_m(:,1)=340.28*4;
    v_m(:,1)=0;
    p_m(:,1)=101325;
    T_m(:,1)=288.16;
    %上边界条件
    u_m(71,:)=340.28*4;
    v_m(71,:)=0;
    p_m(71,:)=101325;
    T_m(71,:)=288.16;
    %前缘点边界条件
    u_m(1,1)=0;
    v_m(1,1)=0;
    
    rho_mm=p_m./T_m/R;
    e_m=Cv*T_m;
    a_m=sqrt(R*gamma*T_m);
  for i=1:71
       for j=1:71
           miu_m(i,j)=DYNVIS(miu0,T0,T_m(i,j));
           k_m(i,j)=THERMC(miu_m(i,j),Cp,Pr);
        end
   end
      for i=1:71
        for j=1:71
            nE1(i,j)=rho_mm(i,j)*u_m(i,j);
            nE2(i,j)=rho_mm(i,j)*u_m(i,j)*u_m(i,j)+p_m(i,j)-TAUXY(u_m,v_m,T_m,miu_m,k_m,i,j,1);
            nE3(i,j)=rho_mm(i,j)*u_m(i,j)*v_m(i,j)-TAUXY(u_m,v_m,T_m,miu_m,k_m,i,j,2);
            nE4(i,j)=(rho_mm(i,j)*(e_m(i,j)+(u_m(i,j)^2+v_m(i,j)^2)/2)+p_m(i,j))*u_m(i,j)-u_m(i,j)*TAUXY(u_m,v_m,T_m,miu_m,k_m,i,j,1)-v(i,j)*TAUXY(u_m,v_m,T_m,miu_m,k_m,i,j,2)+TAUXY(u_m,v_m,T_m,miu_m,k_m,i,j,4);
            nF1(i,j)=rho_mm(i,j)*v_m(i,j);
            nF2(i,j)=rho_mm(i,j)*u_m(i,j)*v_m(i,j)-TAUXY(u_m,v_m,T_m,miu_m,k_m,i,j,2);
            nF3(i,j)=rho_mm(i,j)*v_m(i,j)*v_m(i,j)+p_m(i,j)-TAUXY(u_m,v_m,T_m,miu_m,k_m,i,j,2);
            nF4(i,j)=(rho_mm(i,j)*(e_m(i,j)+(u_m(i,j)^2+v_m(i,j)^2)/2)+p_m(i,j))*v_m(i,j)-u_m(i,j)*TAUXY(u_m,v_m,T_m,miu_m,k_m,i,j,2)-v_m(i,j)*TAUXY(u_m,v_m,T_m,miu_m,k_m,i,j,3)+TAUXY(u_m,v_m,T_m,miu_m,k_m,i,j,5);
        end
      end
    E_Ncombined = cat(3, nE1, nE2, nE3, nE4);
    F_Ncombined = cat(3, nF1, nF2, nF3, nF4);
    partial_U_tnew=zeros(size(U_conbined));
    partial_U_tnew(2:70,2:70,1:4)=-(E_Ncombined(2:70,2:70,1:4)-E_Ncombined(2:70,1:69,1:4))/delta_x-(F_Ncombined(2:70,2:70,1:4)-F_Ncombined(1:69,2:70,1:4))/delta_y;
    U=U_conbined+0.5*(partial_U_tnew+partial_U_t)*delta_tmin;
    for i=2:70
        for j=2:70
         rho_n(i,j)=U(i,j,1);
         u_n(i,j)=U(i,j,2)/U(i,j,1);
         v_n(i,j)=U(i,j,3)/U(i,j,1);
         e_n(i,j)=U(i,j,4)/U(i,j,1)-(u_n(i,j)^2+v_n(i,j)^2)/2;
       %  T_n(i,j)=e_n(i,j)/Cv;
       %  p_n(i,j)=rho_n(i,j)*R*T_n(i,j);
        end
    end 
    %边界条件
    %出口边界条件，不包括上下边界的交点
    p_n(2:70,71)=2*p_n(2:70,70)-p_n(2:70,69);
    u_n(2:70,71)=2*u_n(2:70,70)-u_n(2:70,69);    
    v_n(2:70,71)=2*v_n(2:70,70)-v_n(2:70,69);   
    T_n(2:70,71)=2*T_n(2:70,70)-T_n(2:70,69);   
    %物面边界条件
    u_n(1,:)=0;
    v_n(1,:)=0;
    p_n(1,:)=2*p_n(2,:)-p_n(3,:);
    T_n(1,:)=288.16;
    %进口边界条件
    u_n(:,1)=340.28*4;
    v_n(:,1)=0;
    p_n(:,1)=101325;
    T_n(:,1)=288.16;
    %上边界条件
    u_n(71,:)=340.28*4;
    v_n(71,:)=0;
    p_n(71,:)=101325;
    T_n(71,:)=288.16;
    %前缘点边界条件
    u_n(1,1)=0;
    v_n(1,1)=0;
    
    rho_n=p_n./T_n/R;
    e_n=Cv*T_n;
    a_n=sqrt(R*gamma*T_n);
        u=u_n;
        v=v_n;
        rho=rho_n;
        p=p_n;
        e=e_n;
        T=T_n;
        a=a_n;
  for i=1:71
       for j=1:71
           miu(i,j)=DYNVIS(miu0,T0,T(i,j));
           k(i,j)=THERMC(miu(i,j),Cp,Pr);
       end
   end
  %  if rem(time, 1) == 0
     %    disp(['已迭代到', num2str(time), '次']);
     %   [X,Y]=meshgrid(1:71,1:71);
      %  mesh(X,Y,rho);
  %  end
    
    
end