%{
# 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](http://en.wikipedia.org/wiki/Laminar-turbulent_transition)
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:
$$ \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
%{
# Parameters
The domain consists on a rectangular mesh, where you can adjust its
dimensions through $$Lx$$ and $$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.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);
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(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;
II=speye(3*NN);
% Condition at the top
dir=[l.ctl;l.left;l.top;l.bot;l.cbl]; % where to put Dirichlet on u and v
dirInlet = y_inlet_ind;
loc=[u(dir); v(dir);
u(l.right);
v(l.right)];
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);
V=zeros(NN,1);
P=ones(NN,1);
q0=[U(:);V(:);P(:)];
%{
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.
%}
% TIME MARCHING
% Newton iterations
disp('Newton loop')
ReDer = dt*(ReMax-Re)/ReAccelTime; % Calculates the increase rate of Reynolds
q=q0;
qNm1 = q0;
qM=q;
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; ...
D.x*U+D.y*V];
% 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
f(loc)=C*(q-q0);
A(loc,:)=C;
% convergence test
res=norm(f);
disp([num2str(count) ' ' num2str(res)]);
if (count>maxIter) || (res>1e5); NoConv=1;disp('no convergence');break;else NoConv=0; end
if res