cfd2_lc.m

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

function  CFD2()

N=100;
x=linspace(-1,1,N);   
dx=x(2)-x(1);         %空间步长
k=1.4;
r=0.8;                %步长比

U=zeros(3,N) %%%%%%U是变量%%%%%%

out=zeros(3,N);%%%%%%%输出结果%%%%%%
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);
A=zeros(1,N);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%黎曼问题的几种初值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

flag=4;
    if flag==1                % 两道激波
       Density_left=1.7;
       Density_right=1.2;
       Velocity_left=1;
       Velocity_right=-1;
       P_left=1.3;
       P_right=1.5;
    elseif flag==2           %  一膨胀波一激波
       Density_left=2.0;
       Density_right=1.5;
       Velocity_left=0;
       Velocity_right=0;
       P_left=1.5;
       P_right=0.6;
   elseif flag==3           %  一激波一膨胀波
       Density_left=1.1;
       Density_right=2;
       Velocity_left=0.0;
       Velocity_right=0.0;
       P_left=0.8;
       P_right=1.5;
   elseif flag==4           %  两膨胀波
       Density_left=2.2;
       Density_right=1.6;
       Velocity_left=-0.5;
       Velocity_right=0.5;
       P_left=2.3;
       P_right=1.8;
   end


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

%  间断左边初值数据
   for mark=1:N/2
       U(1,mark)=Density_left;
       U(2,mark)=Density_left*Velocity_left;
       U(3,mark)=P_left/(k-1)+Density_left*Velocity_left*Velocity_left/2;
   end


% 间断右边初值数据
   for mark=1+N/2:N
       U(1,mark)=Density_right;
       U(2,mark)=Density_right*Velocity_right;
       U(3,mark)=P_right/(k-1)+Density_right*Velocity_right*Velocity_right/2;
   end
   


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





for num=1:50                % 迭代步数

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

     rou=U(1,:);
     u=U(2,:)./U(1,:);
     E=U(3,:)./U(1,:);
     e=E-0.5*u.*u;
     p=(k-1)*e.*rou;
     a=sqrt(k*p./rou);
     M=u./a;
     H=E+p./rou;

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

     for j=1:N-1

         uL(j)=u(j);                %%%%将速度u分裂%%%%%%%
         uR(j)=u(j+1);

         rouL(j)=rou(j);            %%%%将密度rou分裂%%%%%%%
         rouR(j)=rou(j+1);

        HL(j)=H(j);                 %%%%将比焓分裂%%%%%%%
        HR(j)=H(j+1); 
     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));  % 此步计算出各个节点上新的变量值，但是注意计算域边界上的值在每次迭代中都没有变化%
  end
        U(:,1)=U(:,2);     %%%%%%%%%%neumann边界条件%%%%%%%%%%%%%
        U(:,N)=U(:,N-1);    %%%%%%%%%%neumann边界条件%%%%%%%%%%%%%

        out(1,:)=rou;
        out(2,:)=u;
        out(3,:)=p;
        plot(x,out,'-p');

        axis([-1 1 -1.5 5])
        legend('Density','velocity','pressure');
        drawnow;


end