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
