nnd.m

NND格式激波管源码,显示出NND格式在激波捕捉方面的优势

%NND格式，时间二阶RUNGE-KUTTA


n=200;
k=0.15/n;
NN=100;
h=1/NN;
alpha=k/h;
gama=1.4;
  for K=1:NN/2
         rou(1,K)=0.445;
         u(1,K)=0.698;
         p(1,K)=3.528;
         E(1,K)=p(1,K)/rou(1,K)/(gama-1)+0.5*u(1,K)*u(1,K);
  end
  for L=NN/2+1:NN
         rou(1,L)=0.5;
         u(1,L)=0;
         p(1,L)=0.571;
         E(1,L)=p(1,L)/rou(1,L)/(gama-1)+0.5*u(1,L)*u(1,L);
  end
  
  for i=1:n  
    for j=1:NN
      U{i,j}=[rou(i,j);rou(i,j)*u(i,j);rou(i,j)*E(i,j)];
    end
    for j=1:NN
      F{i,j}=[rou(i,j)*u(i,j);rou(i,j)*u(i,j)*u(i,j)+p(i,j);(rou(i,j)*E(i,j)+p(i,j))*u(i,j)];
    end
    for j=1:NN
      E(i,j)=p(i,j)/rou(i,j)/(gama-1)+0.5*u(i,j)*u(i,j);
      H(i,j)=E(i,j)+p(i,j)/rou(i,j);
      c(i,j)=sqrt(((gama*p(i,j)/rou(i,j))));    
      R{i,j}=[-(gama-1)/c(i,j)/c(i,j),-0.5/c(i,j),0.5/c(i,j);-(gama-1)*u(i,j)/c(i,j)/c(i,j),-0.5*(u(i,j)-c(i,j))/c(i,j),0.5*(u(i,j)+c(i,j))/c(i,j);-(gama-1)*0.5*u(i,j)*u(i,j)/c(i,j)/c(i,j),-0.5*(H(i,j)-u(i,j)*c(i,j))/c(i,j),0.5*(H(i,j)+u(i,j)*c(i,j))/c(i,j)];
      R_C{i,j}=[0.5*u(i,j)*u(i,j)-c(i,j)*c(i,j)/(gama-1),-u(i,j),1;-u(i,j)-(gama-1)*0.5/c(i,j)*u(i,j)*u(i,j),1+(gama-1)*u(i,j)/c(i,j),-(gama-1)/c(i,j);-u(i,j)+(gama-1)*0.5/c(i,j)*u(i,j)*u(i,j),1-(gama-1)*u(i,j)/c(i,j),(gama-1)/c(i,j)];
      lamda{i,j}=[u(i,j),0,0;0,u(i,j)-c(i,j),0;0,0,u(i,j)+c(i,j)];
      
      A{i,j}=R{i,j}*lamda{i,j}*R_C{i,j};
      A_R{i,j}=[0,1,0;(gama-3)*0.5*u(i,j)*u(i,j),(3-gama)*u(i,j),(gama-1);(gama-1)*u(i,j)*u(i,j)*u(i,j)-gama*E(i,j)*u(i,j),-3*(gama-1)*0.5*u(i,j)*u(i,j)+gama*E(i,j),gama*u(i,j)];
      
      
      %C(i,j)为声速的平方
     % C(i,j)=gama*p(i,j)/rou(i,j);
      %A{i,j}=[0 1 0;0.5*(gama-3)*u(i,j)*u(i,j) (3-gama)*u(i,j) (gama-1);0.5*(gama-2)*(u(i,j))^3-u(i,j)*C(i,j)/(gama-1) C(i,j)/(gama-1)+0.5*(3-2*gama)*(u(i,j))^2 gama*u(i,j)];
      
      
    end
    for j=1:NN
        %[M{i,j},N{i,j}]=eig(A{i,j});
%         for m=1:3
            %if N{i,j}(m,m)>=0          
%                    ALPHA_REAL{i,j}(m,m)=(lamda{i,j}(m,m)+abs(lamda{i,j}(m,m)))/2;
%             %else
%                    ALPHA_FALSE{i,j}(m,m)=(lamda{i,j}(m,m)-abs(lamda{i,j}(m,m)))/2;
% lamda{i,j}
% lamda{i,j}*(lamda{i,j}>zeros(3,3))
% lamda{i,j}*(lamda{i,j}<zeros(3,3))
                    ALPHA_REAL{i,j}=lamda{i,j}*(lamda{i,j}>zeros(3,3));
%             %else
                    ALPHA_FALSE{i,j}=lamda{i,j}*(lamda{i,j}<zeros(3,3));
            %end
%      end
        A_REAL{i,j}=R{i,j}*ALPHA_REAL{i,j}*R_C{i,j};
%         A_FALSE{i,j}=R{i,j}*ALPHA_FALSE{i,j}*R_C{i,j};
        A_FALSE{i,j}=A_R{i,j}-A_REAL{i,j};
        F_REAL{i,j}=A_REAL{i,j}*U{i,j};
        F_FALSE{i,j}=A_FALSE{i,j}*U{i,j};
    end  
%          for j=1:NN
%                DDD(j)=max(max(u));
%          end
  
 
         
%          for j=1:50
%                DDD(j)=max(max(ALPHA_FALSE{i,j}+ALPHA_REAL{i,j}-lamda{i,j}));
%          end
%          disp(i)
%          max(DDD)
         
%     for j=1:200
%         c(i,j)=sqrt(abs(gama*p(i,j)*rou(i,j)));
%         E(i,j)=p(i,j)/rou(i,j)/(gama-1)+0.5*u(i,j)*u(i,j);
%         if u(i,j)>c(i,j)
%             F_REAL{i,j}=[rou(i,j)*u(i,j),p(i,j)+rou(i,j)*u(i,j)*u(i,j),(rou(i,j)*E(i,j)+p(i,j))*u(i,j)]';
%             F_FALSE{i,j}=[0 0 0]';
%         else F_REAL{i,j}=0.5*rou(i,j)/gama*[(2*gama-1)*u(i,j)+c(i,j),2*(gama-1)*u(i,j)*u(i,j)+(u(i,j)+c(i,j))^2,(gama-1)*u(i,j)*u(i,j)*u(i,j)+0.5*(u(i,j)+c(i,j))^3+0.5*(3-gama)*(u(i,j)+c(i,j))*c(i,j)*c(i,j)/(gama-1)]';
%             F_FALSE{i,j}=0.5*rou(i,j)/gama*[u(i,j)-c(i,j),(u(i,j)-c(i,j))*(u(i,j)-c(i,j)),0.5*(u(i,j)-c(i,j))*(u(i,j)-c(i,j))*(u(i,j)-c(i,j))+0.5*(3-gama)*(u(i,j)-c(i,j))*c(i,j)*c(i,j)/(gama-1)]';
%         end
%     end
    for j=1:NN-1
        delta_F_REAL{i,j}=F_REAL{i,j+1}-F_REAL{i,j};
    end
    for j=1:NN-1
        delta_F_FALSE{i,j}=F_FALSE{i,j+1}-F_FALSE{i,j};
    end
    j=NN;delta_F_FALSE{i,j}=[0 0 0]'; delta_F_REAL{i,j}=[0 0 0]';
    for j=3:NN-1
        UU{i+1,j}=U{i,j}-alpha*(F_REAL{i,j}+0.5*minmod(delta_F_REAL{i,j-1},delta_F_REAL{i,j})+F_FALSE{i,j+1}-0.5*minmod(delta_F_FALSE{i,j},delta_F_FALSE{i,j+1})-(F_REAL{i,j-1}+0.5*minmod(delta_F_REAL{i,j-2},delta_F_REAL{i,j-1})+F_FALSE{i,j}-0.5*minmod(delta_F_FALSE{i,j-1},delta_F_FALSE{i,j})));
    end
    j = 1;
    UU{i+1,j} =U{i,j};
    j = 2;
    UU{i+1,j} =U{i,j};
    j = NN;
    UU{i+1,j} =U{i,j};
    
    for j=1:NN   
      rou(i+1,j)=UU{i+1,j}(1);        
      u(i+1,j)=UU{i+1,j}(2)/rou(i+1,j);
      p(i+1,j)=(-0.5*rou(i+1,j)*u(i+1,j)*u(i+1,j)+UU{i+1,j}(3))*(gama-1);
      E(i+1,j)=p(i+1,j)/rou(i+1,j)/(gama-1)+0.5*u(i+1,j)*u(i+1,j);
    end
  end
%   c='.plt';
%   for j=1:n
%      a(j)=j;
%      b=num2str(a(j));
%      d=[b c];
%      fn=fopen(d,'w');
%      for i=1:NN
%          fprintf(fn,' %12.8f\n',u(j,i));
%      end
%      fclose(fn);
%   end
%       MM=moviein(n);
%       XX=linspace(0,1,200);
%       YY=linspace(0,1.005,201);
%       for i=1:n
%           plot(YY,u_MAC(i,:),'r',XX,u(i,:),'b');
%           m(i)=getframe;
%       end
%       movie(m,20)
      
%       
%       MM=moviein(n);
%       XX=linspace(0,1,200);
%       YY=linspace(0,1.005,201);
%       for i=1:n
%           plot(XX,u_MAC(i,:),'r',XX,u_NND(i,:),'b',XX,p_MAC(i,:),'y',XX,p_NND(i,:),'g',XX,rou_MAC(i,:),'c',XX,rou_NND(i,:),'k');
%           m(i)=getframe;
%       end
%       movie(m,20)
%       
%