sandbox/easystab/NavierStokesStationnary

    clear all; clf; format compact; %{

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

    Coded by stab2014/Luis.m and stab2014/PaulValcke.m, on the initial venturi.m code.

    We took the Ventury code, initially created for cylindrical pipe flow. We changed the expression of the jacobian to adapt it to cartesian coordinates.

    This version explain all the structure of the code, there are alternates version of the code here :

    We present here a Poiseuille flow, we developped the code also for :

    1. Poiseuille flow which is how a uniform Flow become a Poiseuille flow

    2. backwardfacing step which is a Poiseuille Flow in a part of the left boundary

    3. Cavity Lid Driven Which is 4 dirichlet condition for evey unknown, all homogenous but one (horizontal speed at the top). A smooth approach for easier convergence at high Reynolds is also presented.

    4. Ventury Which is a reformulation of Heopffner’s (../Ventury.m) code. (in cylindric)

    Used Codes : you may find theses codes in diffpack :

    How does this code works ?

    I you are not totally familiar with the structures of codes in this website, you should read first pedagogy.m

    • First, we define the physical and numerical parameters of our simulation. We then calculate the differentiation matrixes : they contains all the information of the schemes you can use. They replace continous operators like derivates with Matrixes. The type of matrixe you use define the kind of interpolation polynoms you use (Tchebitchev, Fourier, Finite differences…). See diffmat.m and pedagogy.m for more informations
    • We do geometric transformation of our domain, transforming rectangular initial domain to the form we want. See Geometry.m (how to use it) and [map2D.m] (what is the tool) for more informations. New differentiations matrixes are calculated according to plan transformation.
    • We define our Boundary conditions : we select all the index of boundaries in “dir” and this kind of parameters. pedagogy.m will explain you in details the theory
    • We then define our initial guess. All our unknowns are on the same vectors, one after the other one (for Navier-Stokes, we have U,V,P). We stock all of them in the same vector Sol, u,v,p contain only the indexes of U,V,P in Sol. (again, pedagogy.m for more informations).
    • We are now on the convergence loop We define the functions we solve, in f, one after the other one (Here, momentum conservation on x-axis, y-axis, mass conservation). We then define the jacobi of the functions, this matrix contains all the informations about the derivates. It’s a kind of ultimate tool for linearisation, as it contains every slope you may ever need on your problem. We apply our boundary conditions of the Jacobi and the function We then do Newton iterations to find the solution, (see wikipedia page here : ). As the jacobi is the general notion for the derivate, it works the same. We use relaxation method, wich imply more iteration to converge but it is more stable.

    Comments

    We take on purpose a very bad initial guess, just to show that it still converge easily

    We have a constant input, wich imply a zone where the flow is not established.

    We need to publish the figure, to check the form of the velocity, and to check how long (in meters) it takes to correspond to poiseuille flow.

    The code

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %PHYSICAL PARAMETERS Re=100; % reynolds number, it change the impact of the diffusion term compared to pression Uentry=0.1; %GEOMETRIC PARAMETERS Lx = 0.1; Ly = 0.01; %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 %Uinlet = Re1.79e-5/(1.225DiamInlet); pts=5; % number of points in finite difference stencils %%%%%%%%%%%%%%%%%%%%%%%%%%%%%DIFFERENTIATION MATRIXES%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DIFFERENTIATION % We create here our matrixes of differentiatons of our system, on a % rectangle domain. [d.x,d.xx,d.wx,x]=dif1D(‘fd’,0,Lx,Nx,pts); [d.y,d.yy,d.wy,y]=dif1D(‘fd’,0,Ly,Ny,pts); [D,l,X,Y,Z,I,NN]=dif2D(d,x,y); % PLAN TRANSFORMATION %You can transform X and Y as done in geometry % We use Map2D to change our differentiations matrixes according to the % transformation of X and Y D=map2D(X,Y,D); D.lap=D.xx+D.yy; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%BOUNDARIES GESTION%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %We always proceed on the same way : we localise all the coordinates of the %boundaries we are looking for, then we put them in the matrix C. %For Neumann, use D.vector %For DIRICHLET, use vector %%%% POSITIONS OF EVERY VECTORS %u,v,p contain only the position u,v,p will have in SOL. u=(1:NN)‘; v=u+NN; p=v+NN; II=speye(3NN); %%%% EXAMPLE OF DIRICHLET ENTRY dirInlet = [l.left]; % PRESSURE CONDITIONS % Actually, we don’t have any condition on Pressure and it works. % On a poiseuille Flow, it would be smart to define D.xxP instead of % neumann condition, as Dp/Dx=cte in theorical solution. % Neumann condition will change the behaviours close to boundaries, so I % consider it’s better to put no conditions in % BOUNDARY POSITIONS %l.ctl correspond to the point at corner top left, cbr corner bottom right… %l.bot cotains all the indexes of point of the bottom but corner dir=[l.ctl;l.ctr;l.cbl;l.cbr;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); u(l.right); … v(l.right)]; %MATRIX OF CONSTRAINTS
    C=[II([u(dir);v(dir)],:); % … % Dirichlet on u,v 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,:)]; % … % Neumann on v at outflow %%%%%%%%%%%%%%%%%%%%%%%%%%%%%INITIAL GUESS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %It needs to satisfy the boundary conditions, otherwise it will be very %hard to converge to the good results. The best your initial guess is, the %fastest your algorithm will work. If you have no convergence, try to do a %better initial guess, this will help. %Initial condition on U U=(1-((Y(:,1)-Ly/2)/Lyones(1,Nx)).^2); U(dirInlet)=1; %Initial condition on V V=zeros(NN,1); V(dirInlet) = 0; %Initial condition on P P=-X/Re; %Constant gradient of pressure P=P-P(1,1); %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 % We take a mix between the two last iteration on convergence into do a % relaxation. The solution converge slower but more easily. % 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.x
    U; Uy=D.yU; Vx=D.xV; Vy=D.yV; Px=D.xP; Py=D.yP; %DEFINITION OF THE FUNCTIONS YOU WANT TO RESOLVE %if you have some coupled equations, just add them with ; %If they aren’t coupled, copy and paste this section for each equations f=[-(U.Ux+V.Uy)+(D.lapU)/Re-Px; … -(U.Vx+V.Vy)+(D.lapV)/Re-Py; … D.xU+D.yV]; %JACOBIAN OF THE EQUATIONS %separate with , for every unknow of the equation, ; for every 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)/Re, -D.y; … D.x, 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 %Module and pressure variation 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,0,Ly]) quiver3(X,Y,ones(Ny,Nx).max(Module),reshape(U,Ny,Nx),reshape(V,Ny,Nx),… zeros(Ny,Nx),‘k’);
    xlabel(‘x’); ylabel(‘y’); 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,0,Ly])
    xlabel(‘x’); ylabel(‘y’); title(‘Pressure variation Px’); grid off;hold off %colorbar; drawnow %Vertical and horizontal velocity
    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’); % colorbar; axis([0,Lx,0,Ly]); xlabel(‘z’); ylabel(‘y’); title(‘Horizontal velocity U’); grid off;hold off subplot(2,1,2); surf(X,Y,reshape(V,Ny,Nx)); view(2); shading interp; xlabel(‘x’); ylabel(‘y’); title(‘Vertical velocity W’); grid off %colorbar drawnow end %Place here figures you only want once at the end

    Figures

    Here is the result for Poiseuille flow with uniform entry :

    I have problems with colorbar, they reduce the vertical size

    alt text alt text

    Contribute !

    • Please do a mesh convergence
    • Please compare the flow to theory
    • Please compare the length of non-established zone to boundary layer theory
    • Please do some experiments with geometry