sandbox/easystab/stab2014/NavierStokesStationnaryCylindric

    clear all; clf; format compact; %{

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

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

    This is a little change from /NavierStokesStationnary.m, 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 /NavierStokesStationnary.mfor 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.m 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