Simulation of Jet stream

    Coded by stab2014/Luis.m and stab2014/PaulValcke.m. This codes shows the simulation of a flow jet through a pipe. We solved this problem using the Unsteady Navier-Stokes equations. Velocity-pressure non-linear coupling is solved by means of a Newton method, based on the Jacobian matrix. In this code we use a first-order march in time for integrating through time.

    There exist no analytical solution of this flow, so a theoretical validation cannot be performed.

    Preliminary thought about the flow

    This flow is very interesting from the point of view of stability. A jet stream flow through a pipe is naturally an unstable flow. Transition article from Wikipedia shows the main idea, and a similar example: the Osborne Reynold’s experiment. At its initial stage, the flow remains stable, and the mixing layer is horizontal. Then, the jet stream starts oscillating, until its main structure is broken and several vortical structures arise.

    Note that the Unsteady Navier-Stokes equations of this code has been implemented using the following form:

    \displaystyle \partial_tu_i + u_j\partial_ju_i = -\partial_ip + \frac{1}{Re}\partial^2_{jj}u_i

    Therefore, the meaning of the Reynolds number will depend on what the user specifies as dimensionless speed.

    clear all; figure('OuterPosition',[0, 0, 800, 400]); format compact


    The domain consists on a rectangular mesh, where you can adjust its dimensions through \displaystyle Lx and \displaystyle Ly. The jet inlet is placed in the centerline of its length, and its diameter can be modified. Moreover, the inlet flow imposed has the shape of a Poiseuille’s flow.

    Simulating high Reynolds numbers is always tricky, but this code implements a simple way of increasing smoothly (actually linearly) the Reynolds number from its inital state until a user-provided value. The user also specifies the time rate of increase of the Reynolds number.

    %%%% parameters 
    Re=1000; % Initial Reynolds Number.
    ReMax = 6000; % Final Reynolds Number
    ReAccelTime = 5; % Reynolds increase lapse-time.
    Lx = 2; % Length of the pipe
    Ly = 0.6; % Height of the pipe
    DiamInlet = 0.025; % Diameter of jet stream inlet.
    Uinlet = 1; % Flow velocity at inlet. For dimensionless calculation use 1.
    Nx=200; % number of grid nodes in x
    Ny=150; %number of grid nodes in y
    pts=3; % number of points in finite difference stencils
    alpha=1; % Under-relaxation. If the code diverges, reducing alpha may favour convergence.
    maxIter = 200; % Newton's method maximum iterations
    resSTOP = 1e-5; % Newton's residual
    PlCase = 'Vorticity'; % User-specified: plots Vorticity or Velocity
    dt = 0.1; % Time step
    tmax = 30; % Time max
    % differentiation 
    D.lap=D.yy+D.xx; % Laplacian
    % Stores the indices of the inlet and its corresponding values
    y_inlet_ind = find(Y(l.left)<=(Ly+DiamInlet)/2 &...
    y_inlet = Y(y_inlet_ind);

    Boundary conditions

    %%%% preparing boundary conditions
    u=(1:NN)'; v=u+NN; p=v+NN;
    % Condition at the top
    dir=[l.ctl;l.left;;;l.cbl]; % where to put Dirichlet on u and v
    dirInlet = y_inlet_ind;
    loc=[u(dir); v(dir);
    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

    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);
    U(dirInlet) = 1.5.*Uinlet.*(1-((2.*(y_inlet-min(y_inlet)))./DiamInlet - 1).^2);

    Main loop resolution

    Here we perform the march in time. For each time step, a Newton’s method iterative resolution is performed by means of the Jacobian of the Navier-Stokes equations.

    % Newton iterations
    disp('Newton loop')
    ReDer = dt*(ReMax-Re)/ReAccelTime; % Calculates the increase rate of Reynolds
    qNm1 = q0;
    quit=0;count=0;time = 0;
    for t = 0:dt:tmax
    while ~quit     
        % the present solution and its derivatives
        UM = qNm1(u); % Previous time-step solution
        VM = qNm1(v);
        PM = qNm1(p);
        U=alpha.*q(u) + (1-alpha).*qM(u);
        V=alpha.*q(v) + (1-alpha).*qM(v);
        P=alpha.*q(p) + (1-alpha).*qM(p);
        Ux=D.x*U; Uy=D.y*U;
        Vx=D.x*V; Vy=D.y*V; 
        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=[UM-U-(U.*Ux+V.*Uy+Px-(D.lap*U)/Re).*dt; ...
           VM-V-(U.*Vx+V.*Vy+Py-(D.lap*V)/Re).*dt; ...
        % Jacobian 
        A=[-spd(ones(NN,1))+dt.*(-(spd(Ux)+spd(U)*D.x + spd(V)*D.y )+(D.lap)/Re), -spd(Uy).*dt, (-D.x).*dt; ...
            -spd(Vx).*dt,  -spd(ones(NN,1))+dt.*(-(spd(Vy) + spd(V)*D.y + spd(U)*D.x )+(D.lap)/Re), (-D.y).*dt; ...
             D.x, D.y, Z];  
        % Boundary conditions 
        % convergence test  
        disp([num2str(count) '  ' num2str(res)]);
        if (count>maxIter) || (res>1e5); NoConv=1;disp('no convergence');break;else NoConv=0; end
        if res<resSTOP; quit=1; disp('converged'); end
        % Newton step
        qM = q;
    count = 0;
    qNm1 = q;
    quit = 0;

    Plotting the results

    For each time-step, a PNG file is saved of wether the Velocity or Vorticity field (depends on what the user has specified). Then, an animated gif is done.

        % Plot
        switch PlCase
            case 'Vorticity'
                Vorticity = -D.y*U + D.x*V;
                surf(X,Y,reshape(Vorticity,Ny,Nx)); view(2); shading interp; axis equal; colorbar; grid off;
                caxis([-10 10]);
            case 'Velocity'
                Module = sqrt((U.^2+V.^2));
                surf(X,Y,reshape(Module,Ny,Nx)); view(2); shading interp; axis equal; colorbar; grid off;
        if NoConv, PlCase='***NO CONVERGENCE***'; end
        xlabel('x'); ylabel('y'); title(sprintf('Jet. %s field. t=%0.2fs ; Reynolds %0.0f',PlCase,time,Re));
    time =time+dt;
    if Re < ReMax, Re = Re + ReDer; end
    if NoConv, break;end
    % Creates animated gif
    frame = 1;
    for t = 0:dt:tmax
       im = imread(sprintf('Jet_t%0.1fs.png',t));
       [imind,cm] = rgb2ind(im,256);
       if frame ==1
       frame = frame+1;


    The following results are presented for the case Reynolds 6000, \displaystyle Lx= 2 and \displaystyle Ly = 0.6. The resulting velocity and vorticity fields are:

    Jet Stream - Velocity field

    Jet Stream - Velocity field

    Jet Stream - Vorticity field

    Jet Stream - Vorticity field

    There are several remarks that is worth making. First of all, after testing many cases for different Reynolds and pipe diameters, we found that the maximum size of the vortex is approximately equal to the diameter. Small vortices tends to merge between them creating bigger ones. The biggest vortices remains in place longer and are more stable structures than the smaller ones, that tend to disappear or more often to merge.