sandbox/easystab/NavierStokesStationnaryCylindric

    clear all; clf; format compact; %{

    Cartesian ,Stationnary, Uncompressible, Navier-Stokes with Jacobian resolution

    Coded by Stab2014/Lus.m and Stab2014/PaulValcke.m, as a mix between ventury.m an NavierStokesStationnary code.

    This is a little change from NavierStokesStationnary, here you have an axiscylindrical flow instead of a rectangular one. The derivatives on Y are differents in consequences. We just changed the expression of few elements. We also changed the written structure.

    See NavierStokesStationnaryfor more details about the structure of the code.

    Instead of just doing a classic poiseuille, we added a geometry shape of ventury, see ventury.m an geometry for more informations

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %PHYSICAL PARAMETERS Re=1000; % reynolds number, it change the impact of the diffusion term compared to pression %GEOMETRIC PARAMETERS Lx = 30; Ly=2; amp=0.8; % radius at venturi xpos=8; % position of the neck %NUMERICAL PARAMETERS Nx=50; % number of grid nodes in x Ny=10; %number of grid nodes in y alpha=0.7; % Under-relaxation maxIter = 200; % Maximum number of iteration for Newton method. resSTOP = 1e-5; % Precision required to stop convergence %FIXED PARAMETERS OR INTERMEDIATE FORMULAS Ymin=-Ly/2; pts=5; % number of points in finite difference stencils %%%%%%%%%%%%%%%%%%%%%%%%%%%%%DIFFERENTIATION MATRIXES%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DIFFERENTIATION [d.x,d.xx,d.wx,x]=dif1D(‘fd’,0,Lx,Nx,pts); [d.y,d.yy,d.wy,y]=dif1D(‘cheb’,Ymin,Ly,Ny); [D,l,X,Y,Z,I,NN]=dif2D(d,x,y); % PLAN TRANSFORMATION etay=1-(1-amp)exp(-((x-xpos)/3).^2); Y=Y.repmat(etay’,Ny,1);
    D=map2D(X,Y,D); ym1=spd(1./Y); % FORMULAS ARE EASIER TO WRITE WITH THESES ym2=spd(1./Y.^2); D.lap=D.yy+D.xx+ym1D.y; %WE ADD A TERM COMPARED TO CARTESIAN COORDINATES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%BOUNDARIES GESTION%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% POSITIONS OF EVERY VECTORS u=(1:NN)‘; v=u+NN; p=v+NN; II=speye(3NN); % EXAMPLE OF DIRICHLET ENTRY dirInlet = [l.left]; % PRESSURE CONDITIONS neuploc=[l.ctl;l.ctr;l.ctr-Ny]; % where to impose the neumann condition on the pressure p0loc=2Ny; % where to impose zero pressure % BOUNDARY POSITIONS dir=[l.ctl;l.ctr;l.left;l.top;l.bot;]; % where to put Dirichlet on u and v % Loc will be used to apply boundary condition to the function loc=[u(dir); v(dir);p(p0loc); … u(l.right); … v(l.right); … p(neuploc)]; %MATRIX OF CONSTRAINTS
    C=[II([u(dir);v(dir);p(p0loc)],:); % … % Dirichlet on u,v,p D.x(l.right,:), Z(l.right,:), Z(l.right,:); … % Neuman on u at outflow Z(l.right,:), D.x(l.right,:),Z(l.right,:); … -Z(neuploc,:), D.lap(neuploc,:)/Re, -D.x(neuploc,:)]; % Neumann on v at outflow %%%%%%%%%%%%%%%%%%%%%%%%%%%%%INITIAL GUESS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Initial condition on U U=zeros(NN,1); U=(1-2(Y(:,1)/Lyones(1,Nx)).^2); %Initial condition on V V=zeros(NN,1); %Initial condition on P P=-X/Re; %Constant gradient of pressure P=P-P(p0loc); %Change of scale %Initial Solution sol0=[U(:);V(:);P(:)]; %%%%%%%%%%%%%%%%%RESOLUTION ALGORITHM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % NEWTON ALGORITHM sol=sol0; solM=sol; quit=0;count=0; while ~quit
    % RELAXATION METHOD % DEFINITION OF THE UNKNOWS USED U=alpha.sol(u) + (1-alpha).solM(u); V=alpha.sol(v) + (1-alpha).solM(v); P=alpha.sol(p) + (1-alpha).solM(p); % DEFINITION OF THE DERIVATES
    Ux=D.xU; Uy=D.yU; Vx=D.xV; Vy=D.yV; Px=D.xP; Py=D.yP; %DEFINITION OF THE FUNCTIONS YOU WANT TO RESOLVE f=[-(U.Ux+V.Uy)+D.lapU/Re-Px; … -(U.Vx+V.Vy)+(D.lapV-ym2V)/Re-Py; … D.xU+(ym1+D.y)*V]; %JACOBIAN OF THE EQUATIONS A=[-( spd(Ux) + spd(U)D.x + spd(V)D.y )+(D.lap)/Re, -spd(Uy), -D.x; … -spd(Vx), -(spd(Vy) + spd(V)D.y + spd(U)D.x )+(D.lap-ym2)/Re, -D.y; … D.x, ym1+D.y, Z]; %BOUNDARY CONDITIONS f(loc)=C(sol-sol0); A(loc,:)=C; %CONVERGENCE TEST res=norm(f); disp([num2str(count) ’ ’ num2str(res)]); if (count>maxIter) || (res>1e8); disp(’no convergence’);break; end if res<resSTOP; quit=1; disp(’converged’); continue; end %NEWTON STEP FOR RESOLUTION solM = sol; sol=solM - A;%-gmres( A , f, 1000, 1e-7);
    count=count+1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%PLOTTING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Place here figures you want at each step of the iterations %Plot of the Module figure(1) subplot(2,1,1); Module = sqrt((U.2+V.2)); surf(X,Y,reshape(Module,Ny,Nx)); view(2); shading interp; hold on axis ([0,Lx,Ymin,-Ymin]) quiver3(X,Y,ones(Ny,Nx).
    max(Module),reshape(U,Ny,Nx),reshape(V,Ny,Nx),… zeros(Ny,Nx),’k’);
    xlabel(‘x’); ylabel(‘r’); title(‘Module of the speed’); grid off;hold off %colorbar; subplot(2,1,2); surf(X,Y,reshape(Px,Ny,Nx)); view(2); shading interp; hold on axis ([0,Lx,Ymin,-Ymin])
    xlabel(‘x’); ylabel(‘r’); title(‘Pressure variation Px’); grid off;hold off %colorbar; drawnow set(gcf,‘paperpositionmode’,‘auto’) print(‘-dpng’,‘-r80’,‘NSRSUmodpre.png’)
    figure(2) subplot(2,1,1); surf(X,Y,reshape(U-1,Ny,Nx)); view(2); shading interp; hold on sely=1:Ny; selx=1:6:Nx; ww=reshape(V,Ny,Nx); uu=reshape(U,Ny,Nx); quiver3(X,Y,ones(Ny,Nx).
    max(Module),reshape(U,Ny,Nx),reshape(V,Ny,Nx),… zeros(Ny,Nx),‘k’);
    axis([0,Lx,-1,1]); xlabel(‘z’); ylabel(‘r’); title(‘radial velocity U’); grid off;hold off subplot(2,1,2); surf(X,Y,reshape(V,Ny,Nx)); view(2); shading interp; xlabel(‘x’); ylabel(‘r’); title(‘axial velocity W’); grid off drawnow set(gcf,‘paperpositionmode’,‘auto’) print(‘-dpng’,‘-r80’,‘NSRSUvelVel.png’) end %Place here figures you only want once at the end

    Figures

    Here for the parameters written in the code

    alt text alt text