cfd2_lc_muscl.m

有限体积法求解一维欧拉方程

function  CFD2()

N=100
x=linspace(0,200,N);   
dx=x(2)-x(1);        %空间步长
k=1.4;
r=0.8;                %步长比
error=1e-5;            %小量    

S=zeros(3,N);
U=zeros(3,N);  %%%%%%U是变量%%%%%%
out=zeros(3,N-1);%%%%%%%输出结果%%%%%%
D=zeros(1,N-1);
u_average=zeros(1,N-1);
h_average=zeros(1,N-1);
a_average=zeros(1,N-1);
d=zeros(N-1,3);
FR=zeros(3,N-1);
FL=zeros(3,N-1);
LAMDA=zeros(3,3,N-1);
R=zeros(3,3,N-1);
qu=zeros(1,N);
qr=zeros(1,N);
qH=zeros(1,N);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%初值

       Density=1.7;
       Velocity=1;
       P=0.8;
%左边界
       Density1=1.7;
       Velocity1=1;
       P1=0.8;
       
%右边界(亚音流出口）

       P2=0.9;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%将初值赋给给变量U%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 初始化  
   for mark=1:N
       U(1,mark)=Density;
       U(2,mark)=Density*Velocity;
       U(3,mark)=P/(k-1)+Density*Velocity^2/2;
   end
%边界条件

       U(1,1)=Density1;
       U(2,1)=Density1*Velocity1;
       U(3,1)=P1/(k-1)+Density1*Velocity1^2/2; %%%%%%%%%%超音流入口%%%%%%%%%%%%%
       
       
       %U(1,N)=Density;
       %U(2,N)=Density*Velocity;
       U(3,N)=P2/(k-1)+Density*Velocity^2/2;    %%%%%%%%%%亚音流出口%%%%%%%%%%%%
       
       
%%%%%%%%%%%%%虚边界，右边亚音流出口%%%%%%%%%%%%%%%%%%%%%
      % U(:,N+1)=U(:,N);
       
       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





for num=1:2000                %迭代步数

%%%%%%%%%%%%%%%%具体的各个节点处气流参数值%%%%%%%%%%%%%%%%%%%%

     rou=U(1,:);
     u=U(2,:)./U(1,:);
     
     U(3,N)=P2/(k-1)+rou(N)*u(N)^2/2;%%%%%%%%%%%%注意每次迭代都需重新计算%%%%%%%%%%%%%%%%%%%%%%
     E=U(3,:)./U(1,:);
     
     e=E-0.5*u.*u;
     p=(k-1)*e.*rou;
  
     p(N)=P2;
     a=sqrt(k*p./rou);
     M=u./a;
     H=E+p./rou;
     
     
       for j=1:N
         A(j)=(0.0014*x(j)-0.3049)/(-0.0007*x(j)^2+0.3049*x(j)+23.57);
       end
       
      for mark=1:N
       S(1,mark)=A(mark)*rou(mark)*u(mark);
       S(2,mark)=A(mark)*rou(mark)*u(mark)^2;
       S(3,mark)=A(mark)*(k*p(mark)*u(mark)/(k-1)+rou(mark)*u(mark)^3/2);
   end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%估计时间步长%%%%%%%%%%%%%%%%%%%%%%%%
     kmax=0;
     for j=1:N
           kmax=max(kmax,abs(u(j))+a(j));
     end
           dt=r*dx/kmax;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   


         uL(1)=u(1); 
         uR(1)=u(2);             %%%%将速度u分裂%%%%%%%
         uR(N-1)=u(N);
         uL(N-1)=u(N-1);
         
         rouL(1)=rou(1);            %%%%将密度rou分裂%%%%%%%
         rouR(1)=rou(2);
         rouL(N-1)=rou(N-1);
         rouR(N-1)=rou(N);
         
        HL(1)=H(1);                 %%%%将比焓分裂%%%%%%%
        HR(1)=H(2); 
        HL(N-1)=H(N-1);
        HR(N-1)=H(N);


        for j=2:N-1
            
        qu(j)=(u(j+1)-u(j))*(u(j)-u(j-1))/((u(j+1)-u(j))^2+(u(j)-u(j-1))^2+error);
        qr(j)=(rou(j+1)-rou(j))*(rou(j)-rou(j-1))/((rou(j+1)-rou(j))^2+(rou(j)-rou(j-1))^2+error);
        qH(j)=(H(j+1)-H(j))*(H(j)-H(j-1))/((H(j+1)-H(j))^2+(H(j)-H(j-1))^2+error);
        end
        qr(N)=0;
        qu(N)=0;
        qH(N)=0;
       
        
        for j=2:N-1
        uL(j)=u(j)+qu(j)*(u(j+1)-u(j-1))*dx/2;                      %%%%将速度u分裂，插值%%%%%%%
        uR(j)=u(j+1)-qu(j+1)*(u(j+1)-u(j-1))*dx/2;
                
         rouL(j)=rou(j)+qr(j)*(rou(j+1)-rou(j-1))*dx/2;            %%%%将密度rou分裂，插值%%%%%%%
         rouR(j)=rou(j+1)-qr(j+1)*(rou(j+1)-rou(j-1))*dx/2;

        HL(j)=H(j)+qH(j)*(H(j+1)-H(j-1))*dx/2;                 %%%%将比焓分裂，插值%%%%%%%
        HR(j)=H(j+1)-qH(j+1)*(H(j+1)-H(j-1))*dx/2;
   
     end
%%%%%%%%%%%%%%%%%%%%%%Roe平均各参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   for j=1:N-1

         D(j)=sqrt(rouR(j)/rouL(j));
         u_average(j)=(uL(j)+D*uR(j))/(1+D);
         h_average(j)=(HL(j)+D*HR(j))/(1+D);
         a_average(j)=sqrt((k-1)*(h_average(j)-0.5*u_average(j)*u_average(j)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%LAMDA对角阵元素为雅各比系数矩阵A的特征值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        d(j,:)=[abs(u_average(j)),abs(u_average(j)-a_average(j)),abs(u_average(j)+a_average(j))];
        LAMDA(:,:,j)=diag(d(j,:));

%%%%%%%%%%R矩阵为将雅各比系数矩阵A变换为对角阵的矩阵%%%%%%%%%%%%%%
     R(:,:,j)=[(1-k)/(a_average(j)*a_average(j))               -0.5/a_average(j)                                     0.5/a_average(j);
        (1-k)*u_average(j)/(a_average(j)*a_average(j))     (a_average(j)-u_average(j))/2/a_average(j)               (a_average(j)+u_average(j))/2/a_average(j);
        (1-k)*u_average(j)^2/2/(a_average(j)*a_average(j)) (u_average(j)*a_average(j)-h_average(j))/2/a_average(j)  (h_average(j)+u_average(j)*a_average(j))/2/a_average(j)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%             
        UL(:,j)=U(:,j);
        UR(:,j)=U(:,j+1);
        FR(:,j)=[UR(2,j)   (3-k)*UR(2,j)^2/2/UR(1,j)+(k-1)*UR(3,j)    k*UR(2,j)*UR(3,j)/UR(1,j)-(k-1)*UR(2,j)^3/2/UR(1,j)^2]';
        FL(:,j)=[UL(2,j)   (3-k)*UL(2,j)^2/2/UL(1,j)+(k-1)*UL(3,j)    k*UL(2,j)*UL(3,j)/UL(1,j)-(k-1)*UL(2,j)^3/2/UL(1,j)^2]';
        
       f(:,j)=0.5*(FL(:,j)+FR(:,j))+0.5*R(:,:,j)*LAMDA(:,:,j)*inv(R(:,:,j))*(UL(:,j)-UR(:,j));
  end


   for j=2:N-1
        U(:,j)=U(:,j)-dt/dx*(f(:,j)-f(:,j-1))+S(:,j);  % 此步计算出各个节点上新的变量值，但是注意计算域zuo边界上的值在每次迭代中都没有变化%
   end
        U(:,N)=U(:,N-1);

   for j=1:N
        out(1,j)=rou(j);
        out(2,j)=u(j)/((k*p(j)/rou(j))^0.5);
        out(3,j)=p(j);
    end
        plot(x,out,'-p');

        axis([1 200 -0.2 4])
        legend('Density','Ma','pressure');
        xlabel('x');
        ylabel('Ma');
        drawnow;


end