% ************************************************* % resolution of the cavity flow problem * % schema explicite centre d'ordre 2 * % pour la vorticite * % resolution directe du laplacien * % pour la fonction courant. * % ************************************************* % * M. Garbey 00 * % ************************************************* N=Ncell+1; % space grid: if cas_test==1, Lx=2.*pi; Ly=2.*pi; else Lx=1.; Ly=1.; end hx=Lx/(N-1);hy=Ly/(N-1); % time step: deltaT_Diffusion=0.9*1/4*min(hx,hy)^2*Re; deltaT=deltaT_Diffusion; for i=1:N, x(i)=(i-1)*hx; y(i)=(i-1)*hy; end % initialisation: newomega=omega; % boundary condition on the sliding wall: if cas_test==1, % exact analytical solution for j=1:N, g(j)=sin(y(j)); end end if cas_test==2, % bell function for sliding wall velocity for j=1:N, g(j)=-16*x(j)^2*(1.-x(j))^2; end end if cas_test==3, % uniform sliding wall velocity for j=1:N, g(j)=-1.; end end % getting the right hand side corresponding to the BC: if cas_test==1, % exact analytical solution. for i=1:N, for j=1:N, rhsomega(i,j)=4./Re*sin(x(j))*sin(y(i)); end end end if cas_test>1, % numerical solution for sliding wall. rhsomega=zeros(N,N); end if cas_test==1, for i=1:N, for j=1:N, omega_exact(i,j)=2.*sin(x(j))*sin(y(i)); end end psi_exact=0.5.*omega_exact; end temps=0.; % operators for the Poisson problem: M=N-2; M2=(N-2)*(N-2); DC=spalloc(M2,M2,5*M2); for i=1:M2, DC(i,i)=(-2./hx/hx-2./hy/hy); if i>1, DC(i,i-1)=1./hy/hy; end; if i 0, DC(i,i-M)=1./hx/hx; end; end; for i=1:M-1, DC(i*M,i*M+1)=0; DC(i*M+1,i*M)=0; end; if solveur==1, [LP,UP,PP]=lu(DC); end if solveur==2, [LP,UP]=luinc(DC,1e-6); end cpt=0; % ATTENTION test=1e-9; flops(0); test_iteration=tol+1; % initialisation of test. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% kt=0; while test_iteration>tol, % test on steady solution. temps=temps+deltaT; % update time kt=kt+1; % Boundary conditions: if cas_test==1, if CL==2, omega(1,2:N-1)=1./(2.*hy^2)*(psi(3,2:N-1)-8.*psi(2,2:N-1))+3./hy*g(2:N-1); omega(N,2:N-1)=1./(2.*hy^2)*(psi(N-2,2:N-1)-8.*psi(N-1,2:N-1))-3./hy*g(2:N-1); omega(2:N-1,1)=1./(2.*hx^2)*(psi(2:N-1,3)-8.*psi(2:N-1,2))+3./hx*g(2:N-1)'; omega(2:N-1,N)=1./(2.*hx^2)*(psi(2:N-1,N-2)-8.*psi(2:N-1,N-1))-3./hx*g(2:N-1)'; end if CL==1, omega(1,2:N-1)=-2./hy^2*psi(2,2:N-1)+2./hy*g(2:N-1); omega(N,2:N-1)=-2./hy^2*psi(N-1,2:N-1)-2./hy*g(2:N-1); omega(2:N-1,1)=-2./hx^2*psi(2:N-1,2)+2./hy*g(2:N-1)'; omega(2:N-1,N)=-2./hx^2*psi(2:N-1,N-1)-2./hy*g(2:N-1)'; end if CL==3, omega(1,2:N-1)=-3./hy^2*psi(2,2:N-1)-0.5*omega(2,2:N-1)+3./hy*g(2:N-1)+hy/2.*g(2:N-1); omega(N,2:N-1)=-3./hy^2*psi(N-1,2:N-1)-0.5*omega(N-1,2:N-1)-3./hy*g(2:N-1)-hy/2.*g(2:N-1); omega(2:N-1,1)=-3./hx^2*psi(2:N-1,2)-0.5*omega(2:N-1,2)+3./hy*g(2:N-1)'+hy/2.*g(2:N-1)'; omega(2:N-1,N)=-3./hx^2*psi(2:N-1,N-1)-0.5*omega(2:N-1,N-1)-3./hy*g(2:N-1)'-hy/2.*g(2:N-1)'; end end if cas_test>1, if CL==1, omega(1,2:N-1)=-2./hy^2*psi(2,2:N-1)+2./hy*g(2:N-1); omega(N,2:N-1)=-2./hy^2*psi(N-1,2:N-1); omega(2:N-1,1)=-2./hx^2*psi(2:N-1,2); omega(2:N-1,N)=-2./hx^2*psi(2:N-1,N-1); end if CL==2, omega(1,2:N-1)=1./(2.*hy^2)*(psi(3,2:N-1)-8.*psi(2,2:N-1))+3./hy*g(2:N-1); omega(N,2:N-1)=1./(2.*hy^2)*(psi(N-2,2:N-1)-8.*psi(N-1,2:N-1)); omega(2:N-1,1)=1./(2.*hy^2)*(psi(2:N-1,3)-8.*psi(2:N-1,2)); omega(2:N-1,N)=1./(2.*hy^2)*(psi(2:N-1,N-2)-8.*psi(2:N-1,N-1)); end end % finite difference computation of derivatives: for i=2:N-1, for j=2:N-1, psiy(i,j)=(psi(i+1,j)-psi(i-1,j))/(2*hy); psix(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*hx); omegay(i,j)=(omega(i+1,j)-omega(i-1,j))/(2*hy); omegax(i,j)=(omega(i,j+1)-omega(i,j-1))/(2*hx); % check time step for stability condition: CFL=max(max(max(abs(psix))),max(max(abs(psiy)))); if CFL<=0. deltaT=0.01*min(hx,hy); else deltaT=0.25*min(hx,hy)/CFL; end if deltaT>1/4*min(hx,hy)^2*Re, deltaT=deltaT_Diffusion; end convect(i,j)=psiy(i,j)*omegax(i,j)-psix(i,j)*omegay(i,j); omegayy(i,j)=(omega(i+1,j)+omega(i-1,j)-2*omega(i,j))/hy^2; omegaxx(i,j)=(omega(i,j+1)+omega(i,j-1)-2*omega(i,j))/hx^2; diffu(i,j)=1./Re*(omegaxx(i,j)+omegayy(i,j)); end end % Time stepping for vorticity: if verif ==1, newomega(2:N-1,2:N-1)=omega(2:N-1,2:N-1)-deltaT*convect(2:N-1,2:N-1)+deltaT*(diffu(2:N-1,2:N-1)+rhsomega(2:N-1,2:N-1)); else omega(2:N-1,2:N-1)=omega(2:N-1,2:N-1)-deltaT*convect(2:N-1,2:N-1)+deltaT*(diffu(2:N-1,2:N-1)+rhsomega(2:N-1,2:N-1)); end % verification of the vorticity solution: if verif ==1, for i=2:N-1, for j=2:N-1, omegayy(i,j)=(omega(i+1,j)+omega(i-1,j)-2*omega(i,j))/hy^2; omegaxx(i,j)=(omega(i,j+1)+omega(i,j-1)-2*omega(i,j))/hx^2; argu=(newomega(i,j)-omega(i,j))/deltaT; argu=argu-(omegaxx(i,j)+omegayy(i,j))/Re-rhsomega(i,j); argu=argu+psiy(i,j)*omegax(i,j)-psix(i,j)*omegay(i,j); residus1(i,j)=argu; end end omega(2:N-1,2:N-1)=newomega(2:N-1,2:N-1); end histoire(kt,1)=max(max(abs(convect(2:N-1,2:N-1)-diffu(2:N-1,2:N-1)-rhsomega(2:N-1,2:N-1)))); test_iteration=histoire(kt,1); histoire(kt,2)=deltaT; histoire(kt,1) if cas_test==1, histoire(kt,3)=max(max(max(abs(omega-omega_exact))),max(max(abs(psi-psi_exact)))); end % resolution of the Poisson problem: RHSC=reshape(omega(2:N-1,2:N-1),M2,1); % choice is for direct solver or iterative solver: if solveur==1, % LU decomposition solver. % newpsi=DC\RHSC; newpsi=UP\(LP\(PP*RHSC)); end if solveur==2, % preconjugate gradient solver. if choix==0, test=1e-9; else test=min(histoire(kt,1)*min(hx,hy)^(1.5),1e-6); end newpsi=pcg(DC,RHSC,test,10,LP,UP); end psi=zeros(N); psi(2:N-1,2:N-1)=-reshape(newpsi,M,M); % verification of the Poisson solution: if verif ==1, for i=2:N-1, for j=2:N-1, psiyy(i,j)=(psi(i+1,j)+psi(i-1,j)-2*psi(i,j))/hy^2; psixx(i,j)=(psi(i,j+1)+psi(i,j-1)-2*psi(i,j))/hx^2; residus(i,j)=psixx(i,j)+psiyy(i,j)+omega(i,j); end end end end % end time stepping. cpt=cpt+flops; cpt=round(cpt/10^6) figure(1) subplot(2,2,1) contour(omega,10); title('vorticity'); subplot(2,2,2) contour(psi,10); title('courant'); subplot(2,2,3); plot(log10(histoire)); title('steady solution?'); if verif==1, subplot(2,2,4); mesh(residus); title('verification?'); pause subplot(2,2,4); mesh(residus1); title('verification?'); else subplot(2,2,4); mesh(convect-diffu); title('stationary solution?'); end % computation of the numerical error when the analytical solution is available: if cas_test==1, AA=max(max(abs(omega-omega_exact))) BB=max(max(abs(psi-psi_exact))) end