Nouveau-dossier-(5).rar_instability_kelvin

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Résolution des équations de Navier-Stokes incompressibles %%% instationnaires : %%% - 2D, différences finies centrées ordre 2, grille cartésienne, colocation %%% - schéma en temps basé sur la formulation PPE (article H. Johnston -J.-G. Liu, JCP 2004) %%% - terme de diffusion explicite en temps %%% - cas test : instablités de Kelvin-Helmoltz, CL périodiques (tiré du livre "Introduction au calcul scientifique par la pratique", de Danaila et al., Dunod) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % parametres Lx=2; Ly=1; imax=65; jmax=65; %imax=10; jmax=10; dx=Lx/(imax-1); dy=Ly/(jmax-1); x=0:dx:Lx; y=0:dy:Ly; rey=1000; nu=1/rey; % CI = jet U0=1; Rj=Ly/4; Pj=20; Ax=0.25; lamx=0.5*Lx; [xx,yy]=meshgrid(x,y);xx=xx';yy=yy'; rr=yy-Ly/2; ux=U0*0.5d0*(1.d0+tanh(0.5*Pj*(1-abs(rr)/Rj))); ux=ux+Ax*ux.*sin(2*pi/lamx*xx); uy=0*ux; p=0*ux; %--- marqueur m=2.d0+tanh(0.5*Pj*(1-abs(rr)/Rj)); %--- indices i=1:imax; ip=i+1; ip(imax)=1; ip2=i+2; ip2(imax-1)=1; ip2(imax)=2; im=i-1; im(1)=imax; im2=i-2; im2(1)=imax-1; im2(2)=imax; j=1:jmax; jp=j+1; jp(jmax)=1; jp2=j+2; jp2(jmax-1)=1; jp2(jmax)=2; jm=j-1; jm(1)=jmax; jm2=j-2; jm2(1)=jmax-1; jm2(2)=jmax; %--- projection du champ de vitesse initial %-- eq pression : 2nd membre divu=(ux(ip,j)-ux(im,j))/(2*dx)+(uy(i,jp)-uy(i,jm))/(2*dy); rhs=-divu; %-- eq pression : mise sous forme vecteur du 2nd membre rhs=reshape(rhs',imax*jmax,1); %-- eq pression : matrice pmat=sparse(zeros(imax*jmax)); for k=1:imax for l=1:jmax ind=(i(k)-1)*jmax+j(l); indip1=(ip(k)-1)*jmax+j(l); indim1=(im(k)-1)*jmax+j(l); indjp1=(i(k)-1)*jmax+jp(l); indjm1=(i(k)-1)*jmax+jm(l); pmat(ind,indim1)=1/dx^2; pmat(ind,indjm1)=1/dy^2; pmat(ind,ind)=-2/dx^2-2/dy^2; pmat(ind,indjp1)=1/dy^2; pmat(ind,indip1)=1/dx^2; if k==1 & l==1 pmat(ind,ind)=0; end end end %-- r�solution eq poisson p=pmat\rhs; p=reshape(p,jmax,imax)'; %-- projection ux=ux+(p(ip,j)-p(im,j))/(2*dx); uy=uy+(p(i,jp)-p(i,jm))/(2*dy); divu=(ux(ip,j)-ux(im,j))/(2*dx)+(uy(i,jp)-uy(i,jm))/(2*dy); divumax(1)=max(max(abs(divu))); %--- boucle en temps t=0; for n=1:210 %-- cfl dt=0.9/max(max(abs(ux)/dx+abs(uy)/dy)); t=t+dt %-- eq pression : 2nd membre rhs=-1/(2*dx)*( ux(ip,j).*(ux(ip2,j)-ux(i,j))/(2*dx) - ux(im,j).*(ux(i,j)-ux(im2,j))/(2*dx) ) ... -1/(2*dx)*( uy(ip,j).*(ux(ip,jp)-ux(ip,jm))/(2*dy) - uy(im,j).*(ux(im,jp)-ux(im,jm))/(2*dy)) ... -1/(2*dy)*( ux(i,jp).*(uy(ip,jp)-uy(im,jp))/(2*dx) - ux(i,jm).*(uy(ip,jm)-uy(im,jm))/(2*dx)) ... -1/(2*dy)*( uy(i,jp).*(uy(i,jp2)-uy(i,j))/(2*dy) - uy(i,jm).*(uy(i,j)-uy(i,jm2))/(2*dy)) ... ; %-- eq pression : mise sous forme vecteur du 2nd membre rhs=reshape(rhs',imax*jmax,1); %-- eq pression : matrice pmat=sparse(zeros(imax*jmax)); for k=1:imax for l=1:jmax ind=(i(k)-1)*jmax+j(l); indip1=(ip(k)-1)*jmax+j(l); indim1=(im(k)-1)*jmax+j(l); indjp1=(i(k)-1)*jmax+jp(l); indjm1=(i(k)-1)*jmax+jm(l); pmat(ind,indim1)=1/dx^2; pmat(ind,indjm1)=1/dy^2; pmat(ind,ind)=-2/dx^2-2/dy^2; pmat(ind,indjp1)=1/dy^2; pmat(ind,indip1)=1/dx^2; if k==1 & l==1 pmat(ind,ind)=0; end end end %-- résolution eq poisson p=pmat\rhs; p=reshape(p,jmax,imax)'; %-- eq vitesse (schema RK4 pour la stabilite, donne de bien meilleurs % resultats que EE) [dux1,duy1]=SMu(ux,uy,p,dx,dy,i,j,ip,jp,im,jm,nu); [dux2,duy2]=SMu(ux+dt/2*dux1,uy+dt/2*duy1,p,dx,dy,i,j,ip,jp,im,jm,nu); [dux3,duy3]=SMu(ux+dt/2*dux2,uy+dt/2*duy2,p,dx,dy,i,j,ip,jp,im,jm,nu); [dux4,duy4]=SMu(ux+dt*dux3,uy+dt*duy3,p,dx,dy,i,j,ip,jp,im,jm,nu); ux=ux+dt/6*(dux1+2*dux2+2*dux3+dux4); uy=uy+dt/6*(duy1+2*duy2+2*duy3+duy4); %-- calcul et affichage vorticite w=(uy(ip,j)-uy(im,j))/(2*dx)-(ux(i,jp)-ux(i,jm))/(2*dy); %- affichage matlab figure(1); pcolor(x,y,w'); axis equal;shading('interp');drawnow %- sortie visit % fich=strcat('SORTIES/resultats.',num2str(n),'.vtk'); % fid = fopen(fich,'w'); % fprintf(fid,'# vtk DataFile Version 2.0\n'); % fprintf(fid,'convection 2D\n'); % fprintf(fid,'ASCII\n'); % fprintf(fid,'DATASET RECTILINEAR_GRID\n'); % fprintf(fid,'%s %d %d %d\n','DIMENSIONS ',imax,jmax,1); % fprintf(fid,'%s %d %s\n','X_COORDINATES',imax,' double'); % fprintf(fid,'%12.8f\n',x); % fprintf(fid,'%s %d %s\n','Y_COORDINATES',jmax,' double'); % fprintf(fid,'%12.8f\n',y); % fprintf(fid,'%s %d %s\n','Z_COORDINATES',1,' double'); % fprintf(fid,'%12.8f\n',0); % fprintf(fid,'%s %d\n','POINT_DATA',imax*jmax); % fprintf(fid,'%s\n','SCALARS w double'); % fprintf(fid,'%s\n','LOOKUP_TABLE default'); % for l=1:jmax % for k=1:imax % fprintf(fid,'%12.8f\n',w(k,l)); % end % end %-- calcul divergence divu=(ux(ip,j)-ux(im,j))/(2*dx)+(uy(i,jp)-uy(i,jm))/(2*dy); divumax(n+1)=max(max(abs(divu))); %-- evolution et affichage du marqueur % (schema decentre EE) m=m-dt* ( ((ux(ip,j)-abs(ux(ip,j)))/2.*m(ip,j)+(ux(ip,j)+abs(ux(ip,j)))/2.*m(i,j))/dx ... -((ux(i,j)-abs(ux(i,j)))/2.*m(i,j)+(ux(i,j)+abs(ux(i,j)))/2.*m(im,j))/dx) ... -dt* ( ((uy(i,jp)-abs(uy(i,jp)))/2.*m(i,jp)+(uy(i,jp)+abs(uy(i,jp)))/2.*m(i,j))/dy ... -((uy(i,j)-abs(uy(i,j)))/2.*m(i,j)+(uy(i,j)+abs(uy(i,j)))/2.*m(i,jm))/dy); %-affichage matlab figure(3) pcolor(x,y,m'); axis equal;shading('interp');drawnow %- sortie matlab % fprintf(fid,'%s\n','SCALARS m double'); % fprintf(fid,'%s\n','LOOKUP_TABLE default'); % for l=1:jmax % for k=1:imax % fprintf(fid,'%12.8f\n',m(k,l)); % end % end % fclose(fid); end %--- affichage divergence figure(2) plot(1:n+1,divumax,'o-');