# sandbox/easystab/free_surface_navier_stokes.m

This is a viscous flow described with the Navier-Stokes equation. The top is a free surface that can deform accordingly to the pressure jump through the interface due to surface tension, and the normal and tangential stress in the flow. At the inflow we impose a profile that looks like a boundary layer.

clear all; clf; format compact

%paramvec=linspace(5,1,7);
paramvec=linspace(1,0.85,10);
%paramvec=0.95

for tre=1:length(paramvec);

disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
param=paramvec(tre)

% parameters
Re=500; % Reynolds number
Nx=100;  % number of grid nodes in z
Ny=50;  %number of grid nodes in r
Lx=10;  % domain length
Ly=1;   % domain height
sig=5;  % surface tension
V=Lx*Ly*param; % desired volume
umax=1; % max inflow velocity
mu=1/Re;

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

% useful
l.top=[l.ctl; l.top; l.ctr]; % include the corners to the top
l.u=(1:NN)'; l.v=l.u+NN; l.p=l.v+NN; l.e=l.p(end)+(1:Nx)'; l.p0=l.e(end)+1; % location vectors
l.ue=[l.u;l.e]; l.ve=[l.v;l.e]; l.pe=[l.p;l.e]; % compound location vectors

p0loc=2*Ny; % where to impose pressure=p0
t=l.top; % a shortcut because its often used

n=3*NN+Nx+1; % numbers of degrees of freedom
ZZ=spalloc(n,n,0);    zx=spalloc(Nx,Nx,0);
II=speye(n);        ix=speye(Nx);
Iu=II(l.u,:); Iv=II(l.v,:); Ip=II(l.p,:); Ie=II(l.e,:); Ip0=II(l.p0,:);
T=kron(speye(Nx,Nx),ones(Ny,1)); % transformation operator from 1D x to 2D grid
yy=Y(:);

% initial guess
v=zeros(NN,1);
%u=umax*Y(:).*(2-Y(:));
u=umax*tanh(Y/0.1);
p=-mu*X*umax; p=p-p(p0loc);
e=Ly+0*sin(pi*x);
p0=0;
q0=[u(:);v(:);p(:);e;p0];

% Newton iterations
disp('Newton loop'); tlu=0;
if tre==1; q=q0; end
quit=0;count=0;
while ~quit
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the present solution and its computational-space derivatives
e=q(l.e); ex=d.x*e; exx=d.xx*e; E=T*e; Ex=T*ex; Exx=T*exx; a=1+ex.^2;
u=q(l.u); v=q(l.v); p=q(l.p); p0=q(l.p0);
ux=Dc.x*u; uy=Dc.y*u; uxy=Dc.y*ux; uyy=Dc.yy*u;
vx=Dc.x*v; vy=Dc.y*v; vxy=Dc.y*vx; vyy=Dc.yy*v;
px=Dc.x*p; py=Dc.y*p;


# Physical space differentiation matrices

This is just like we did in stretching_formula.m.

        % physical-space differentation matrices
D.lap=Dc.xx ...
+spd(-2*yy.*E.^-1.*Ex)*(Dc.x*Dc.y) ...
+spd(yy.^2.*E.^-2.*Ex.^2+E.^-2)*Dc.yy ...
+spd(yy.*(2*E.^-2.*Ex.^2-E.^-1.*Exx))*Dc.y;
D.x=Dc.x+spd(-yy.*E.^-1.*Ex)*Dc.y;
D.y=spd(E.^-1)*Dc.y;


# Jacobians of the differentiation matrices

Since the domain deformation is one of the unknown of th system, the differentiation is no longer a linear operation (because it depends nonlinearly on the domain stretching \eta). We do just like we did in the test case jacobian_formula.m.

        % jacobians of the physical-space differentiation operators
AA=-2*spd(yy)*T*(diag(e.^-1)*d.x-diag(e.^-2.*ex));
BB=2*(spd(yy.^2)*T*(diag(e.^-2.*ex)*d.x-diag(e.^-3.*ex.^2))-T*diag(e.^-3));
CC=spd(yy)*T*(diag(4*e.^-2.*ex)*d.x-diag(e.^-1)*d.xx+diag(-4*e.^-3.*ex.^2+e.^-2.*exx));
DD=T*(diag(-e.^-1)*d.x+diag(e.^-2.*ex));
EE=T*diag(-e.^-2);

D.xue=spd(yy.*uy)*DD;  D.yue=spd(uy)*EE;  D.lapue=spd(uxy)*AA+spd(uyy)*BB+spd(uy)*CC;
D.xve=spd(yy.*vy)*DD;  D.yve=spd(vy)*EE;  D.lapve=spd(vxy)*AA+spd(vyy)*BB+spd(vy)*CC;
D.xpe=spd(yy.*py)*DD;  D.ype=spd(py)*EE;

Iux=[D.x,D.xue]*II(l.ue,:);  Iuy=[D.y,D.yue]*II(l.ue,:);  Iulap=[D.lap,D.lapue]*II(l.ue,:);
Ivx=[D.x,D.xve]*II(l.ve,:);  Ivy=[D.y,D.yve]*II(l.ve,:);  Ivlap=[D.lap,D.lapve]*II(l.ve,:);
Ipx=[D.x,D.xpe]*II(l.pe,:);  Ipy=[D.y,D.ype]*II(l.pe,:);

% physical-space derivatives
ux=D.x*u; uy=D.y*u;
vx=D.x*v; vy=D.y*v;
px=D.x*p; py=D.y*p;


# The nonlinear function

The first thing we do is to compute the expression of the normal and tangential stresses at the free surface. The stress tensor in a viscous Newtonian fluid is \displaystyle \sigma= \begin{pmatrix} 2\mu u_x-p & \mu(u_y+v_x) \\ \mu(u_y+v_x) & 2\mu v_y-p \end{pmatrix} We are interested at the components of the stress along our free-surface. We first need to find the expresion of its normal and tangent vectors. We can say that the normal is along the gradient of a function for which the free-surface is an iso-line. For instance take \displaystyle F(x,y)=y-\eta(x) then its gradient, normalized to unit one is the normal to the interface \displaystyle \overrightarrow{n}= \begin{pmatrix} -\frac{\eta_x}{\sqrt{1+\eta_x^2}}\\ \frac{1}{\sqrt{1+\eta_x^2}}\\ \end{pmatrix} and the tangent vector is \displaystyle \overrightarrow{t}= \begin{pmatrix} \frac{1}{\sqrt{1+\eta_x^2}}\\ \frac{\eta_x}{\sqrt{1+\eta_x^2}}\\ \end{pmatrix} Then, the normal stress is \displaystyle \sigma_n=\overrightarrow{n}\sigma \overrightarrow{n} =\frac{\eta_x^2(2\mu u_x-p)-2\mu\eta_x(u_y+v_x)+2\mu v_y-p}{1+\eta_x^2} and the tangential stress is \displaystyle \sigma_t=\overrightarrow{t}\sigma \overrightarrow{n} =\mu\frac{(1-\eta_x^2)(u_y+v_x)+2\mu\eta_x(v_y-u_x)}{1+\eta_x^2} The conditions at the interface are that the normal stress has a jump equal to the pressure jump through a curved interface with surface tension, and the tangential stress is zero.

We have then Navier-Stokes in the flow domain. And the last equation imposes the value of the fluid volume.

        % nonlinear function
stressn=-sig*(exx.*a.^-1.5)+(ex.^2.*(2*mu*ux(t)-p(t))-2*mu*ex.*(uy(t)+vx(t))+2*mu*vy(t)-p(t))./a; % normal stress
stresst=mu*((1-ex.^2).*(uy(t)+vx(t))+2*ex.*(vy(t)-ux(t)))./a; % tangential stress

f=[-(u.*ux+v.*uy)+D.lap*u/Re-px; ...
-(u.*vx+v.*vy)+D.lap*v/Re-py; ...
ux+vy; ...
stressn; ... % jump in normal stress
d.wx*Ly*e-V]; % volume

# The Jacobian

We build here the expression of the Jacobian of f.

        % analytical Jacobian
Astressn=-sig*(spd(a.^-1.5)*d.xx-spd(3*ex.*exx.*a.^-2.5)*d.x)*Ie ...
+diag(2*mu*ex.^2./a)*Iux(t,:)-diag(2*mu*ex./a)*Iuy(t,:) ...
-diag(2*mu*ex./a)*Ivx(t,:)+diag(2*mu./a)*Ivy(t,:) ...
-Ip(t,:);

Astresst=mu*(diag((1-ex.^2)./a)*Iuy(t,:)-diag(2*ex./a)*Iux(t,:) ...
+diag((1-ex.^2)./a)*Ivx(t,:)+diag(2*ex./a)*Ivy(t,:));

A=[-(spd(u)*Iux+spd(ux)*Iu+spd(v)*Iuy+spd(uy)*Iv)+Iulap/Re-Ipx; ...
-(spd(u)*Ivx+spd(vx)*Iu+spd(v)*Ivy+spd(vy)*Iv)+Ivlap/Re-Ipy; ...
Iux+Ivy; ...
Astressn; ...
d.wx*Ly*Ie];

# Boundary conditions

They are:

• Dirichlet for u and v at inflow (imposed inflow profile set by the initial guess) and bottom (no-slip)
• Neumann for u and v at outflow (x derivative is zero)
• The pressure must equal the reference pressure p_0 close to the bottom left corner
• The Neumann condition on pressure to avoid spurious pressure modes
• The height of the interface is imposed at x=0 and x=L_x.
• The free surface is a streamline, or said differently: no penetration through the interface: \displaystyle v-u\eta_x=0
• The tangential stress is zero
        % Boundary conditions
dir=[l.u([l.left;l.bot]); l.v([l.left;l.bot])]; % where to put Dirichlet on u and v
neuploc=[l.ctl;l.ctr;l.ctr-Ny];  % where to impose the neumann condition on the pressure

loc=[dir; ...
l.u(l.right); ...
l.v(l.right); ...
l.p(p0loc); ...
l.p(neuploc); ...
l.e([1,Nx]); ...
l.u(t);
l.v(t)];

f(loc)=[q(dir)-q0(dir); ... % all Dirichlet
ux(l.right); ... % Homogeneous Neuman on u at outflow
vx(l.right); ...
p(p0loc)-p0; % pressure equal to p0 at p0loc
D.lap(neuploc,:)*v/Re-px(neuploc); ... % neumann constraint on pressure
e(1)-Ly; ... % free-surface pinned at both ends
e(Nx)-Ly; ...
v(t)-u(t).*ex; ... % free-surface is a streamline
stresst];

C=[II(dir,:); ...     % all Dirichlet
Iux(l.right,:); ...
Ivx(l.right,:); ...
Ip(p0loc,:)-Ip0;  % pressure equal to p0 at p0loc
Ivlap(neuploc,:)/Re-Ipx(neuploc,:); ... % Neumann constraint on pressure
ix(1,:)*Ie; ... % free-surface pinned at both ends
ix(Nx,:)*Ie; ...
Iv(t,:)-diag(ex)*Iu(t,:)-spd(u(t))*d.x*Ie; ... % free-surface is a streamline
Astresst];

A(loc,:)=C;
tjac=toc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% plotting
Yp=Y.*reshape(E,Ny,Nx); % physical mesh
ww=reshape(u,Ny,Nx); vv=reshape(v,Ny,Nx);
vo=vx-uy; % vorticity
surf(X,Yp,reshape(vo,Ny,Nx)-10,'facealpha',1); view(2); shading interp; hold on
sely=1:2:Ny; selx=round(linspace(1,Nx,10));
plot(x,Ly+0*x,'b--',Lx+y*0,y,'b--');
axis([0,1.1*Lx,0,1.2*Ly]);
xlabel('x'); ylabel('y'); title('Vorticity and velocity field'); grid off;hold off
%caxis([-30,0]-10)
drawnow

% convergence test
res=norm(f);
disp([num2str(count) '  ' sprintf('%9.3g',res) '  tjac:' sprintf('%7.3g',tjac) '  tlu:' sprintf('%7.3g',tlu)]);
if count>50||res>1e5; disp('no convergence');break; end
if res<1e-5; quit=1; disp('converged'); continue; end

% Newton step
tic; q=q-A\f; tlu=toc;
count=count+1;
end

end
set(gcf,'paperpositionmode','auto')
print('-dpng','-r80','free_surface_navier_stokes.png')


Here is the final figure: