sandbox/easystab/venturi.m

    The Venturi flow

    This is a flow in a cylindrical pipe with a localized neck. So we have the Navier-Stokes equations and we solve for a steady state with the proper boundary conditions. The neck is described using a mapping of the mesh just like we did in diffmat_mapping.m, and peristalsis.m. The new thing here is that it is done, like for pipe_sym.m and pipe.m, in cylindrical coordinates.

    This flow is interesting, becasue there is a quick accelleration of the fluid through the neck, and then detachment if the Reynolds number is large enough. This is strongly the case with the present parameters.

    Dependency:

    clear all; clf; format compact

%%%% parameters 
Re=1000; % reynolds number
Nx=101; % number of grid nodes in z
Ny=40; %number of grid nodes in r
Lx=30; % length in z of the domain [0,Lz]
pts=5; % number of points in finite difference stencils
amp=0.5; % radius at venturi
xpos=8; % position of the neck

% differentiation 
[d.x,d.xx,d.wx,x]=dif1D('fd',0,Lx,Nx,5);
[d.y,d.yy,d.wy,y]=dif1D('cheb',-1,2,Ny);
[D,l,X,Y,Z,I,NN]=dif2D(d,x,y);

% mapping of the mesh
disp('Mapping of the mesh')
etay=1-(1-amp)*exp(-((x-xpos)/3).^2); 
Y=Y.*repmat(etay',Ny,1);  
D=map2D(X,Y,D);

% cylindrical Laplacian
ym1=spd(1./Y);   ym2=spd(1./Y.^2);
D.lap=D.yy+D.xx+ym1*D.y;

    Boundary conditions

    Since there are many boundary conditions, we prepare the vector loc and the constraint matrix C before the loop. Since the boundary conditions are linear, we can build them once for all before the Newton iterations.

    The boundary conditions are Dirichlet everywhere for u (radial velocity) and w (axial velocity), and homogeneous Neumann on u,w out the outflow. For the pressure we impose its value 0 at the second gridpoint from the left of the top boundary, this is teold in p0loc.

    Then we need like we did for peristalsis.m to impose some additional constraints on the pressure, these are imposed for the gridpoints in neuploc. how to impose these conditions is still a little mysterious to me, I do trial and errors until I get something that work. Please contribute if you have ideas more clear that this and write to me if you want to learn more.

    A new thing in this code is that we use sparse matrices, this is very good to reduce the memory usage for our large 2D differentiation matrices, and reduces greatly the computation cost.

    %%%% preparing boundary conditions
NN=Nx*Ny;
u=(1:NN)'; w=u+NN; p=w+NN;
II=speye(3*NN);

neuploc=[l.ctl;l.ctr;l.ctr-Ny];  % where to impose the neumann condition on the pressure
p0loc=2*Ny; % where to impose zero pressure
dir=[l.ctl;l.ctr;l.left;l.top;l.bot]; % where to put Dirichley on u and w

loc=[u(dir); w(dir); p(p0loc); ...
    u(l.right); ...
    w(l.right); ...
    p(neuploc)];

C=[II([u(dir);w(dir);p(p0loc)],:); ...     % Dirichlet on u,w,and 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,:); ...    % Neumann on w at outflow
   Z(neuploc,:), D.lap(neuploc,:)/Re, -D.x(neuploc,:)]; % neuman constraint on pressure

    We chose an initial guess that statisfies the boundary conditions, this is good for Newton, and it makes it also very easy to impose the nonhomogeneous boundary conditions in the lop.

    % initial guess
U=zeros(NN,1);
W=(1-(Y(:,1)*ones(1,Nx)).^2); % mean velocity 1 on the pipe of diameter 1,
P=-X/Re; P=P-P(p0loc); % pressure zero at p0loc
sol0=[U(:);W(:);P(:)];

% Newton iterations
disp('Newton loop')
sol=sol0;
quit=0;count=0;
while ~quit     
 
    % the present solution and its derivatives
    U=sol(u); W=sol(w); P=sol(p);
    Ux=D.x*U; Uy=D.y*U;
    Wx=D.x*W; Wy=D.y*W; 
    Px=D.x*P; Py=D.y*P;

    This is now the heart of the code: the expression of the nonlinear fonction that should become zero, and just after, the expression of its Jacobian, then the boundary conditions.

        % nonlinear function
    f=[-(U.*Uy+W.*Ux)+(D.lap*U-ym2*U)/Re-Py; ...
       -(U.*Wy+W.*Wx)+D.lap*W/Re-Px; ...
      (ym1+D.y)*U+D.x*W];
    
    % Jacobian 
    A=[-(spd(W)*D.x+spd(U)*D.y+spd(Uy))+(D.lap-ym2)/Re, -spd(Ux), -D.y; ...
         -spd(Wy), -(spd(W)*D.x+spd(Wx)+spd(U)*D.y)+D.lap/Re, -D.x; ...
         ym1+D.y, D.x, Z];
     
    % Boundary conditions 
    f(loc)=C*(sol-sol0);
    A(loc,:)=C;

    Plotting

        % plotting
    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(W,Ny,Nx); uu=reshape(U,Ny,Nx); 
    quiver(X(sely,selx),Y(sely,selx),ww(sely,selx),uu(sely,selx),'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(W,Ny,Nx)); view(2); shading interp; 
    xlabel('x'); ylabel('y'); title('axial velocity W'); grid off
    drawnow
    
    % convergence test  
    res=norm(f);
    disp([num2str(count) '  ' num2str(res)]);
    if count>50|res>1e5; disp('no convergence');break; end
    if res<1e-5; quit=1; disp('converged'); continue; end
    
    % Newton step
    sol=sol-A\f;   
    count=count+1;
end


set(gcf,'paperpositionmode','auto')
print('-dpng','-r80','venturi.png')
    The velocity field

    The velocity field

    Validation

    Here we should compare these results to a Gerris simulation, or to the method by Pierre-Yves Lagree with the boundary layer equations.

    Exercices/Contributions

    • Please add as well a mapping on the z cordinate to be able to have better resolution close to the neck. For the moment, this resolution is the limiting factor for how steep the neck can be. you can do that just before the call to the function mapping2D.m, and your modifications will be take into account in the mapped differentiation matrices.
    • Please simulate only the top half of the domain, with a symmetry boundary conditions along the axis, as is suggested also for pipe_sym.m. This is divide by two the ammount of gridpoints. You can also try to use the symmetry method of the differentiation matrices that is used in pipe_sym.m.
    • Please compare these calculations with results from gerris with the same geometry—-> venturi_gerris.gfs.
    • Comparison between the Navier-Stokes equations and the Boundary-layer equations ——-> venturi_pyl.m.

    %}