sandbox/easystab/free_surface_2D.m

    Free surface oscillations in 2D (with surface tension)

    We have a perfect fluid on a rectangular box with the free surface at the top. Here we chose to not have the Navier-Stokes equations for the fluid but just describe its flow with a velocity potential (and the Bernoulli equation for the relation between the flow and the pressure). Thus, all the dynamics is just the dynamics of the free surface, and the flow is just a constraint for the evolution of this free surface. On the contrary to free_surface_gravity.m, here the oscillations of the free surface are not gravity waves, they are just capillary waves, so the theory is a little different from (and complementary to) free_surface_gravity.m#validation.

    See free_surface_2D_functions.m for the same thing but writen with the functions dif1D.m and dif2D.m.

    Dependency:

    • chebdif.m for the Chebychev differentiation matrices
    clear all; clf

% parameters
Nx=20; % gridpoints in x
Ny=20; % gridpoints in y 
Lx=1; % domain size in x
Ly=1; % domain size in y
pts=5; % number of points in finite difference stencils
sigma=1; % surface tension
rho=1; % fluid density

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

    Location vectors

    As it is for venturi.m, and maybe here a little bit more, one of the central aspects of building the system matrices and imposing the boundary conditions has a lot to knowing where the different variables are in the state vector and in the matrices. This is done using the “location vectors”, which are nothing else than a memory of the indices that correspond to the different parts of the state variable q. Here we have the velocity potential \phi(x,y,t) defined on a 2D grid, and we have \eta the position of the free surface, defined just along the top of the rectagular box (it is a 1D variable). So one nice aspect of the present example is that we have to combine variables of different dimensionalities.

    Here the top of the box are the grid points at the top and excluding the corners. The corners are instead included to the right and left location vectors. This is good because we impose boundary conditions for the free surface at the corners, so we can use top just to access the inside points of the top of the mesh.

    Implicitely, the top, bottom, left and right location vectors pertain to the velocity potential \phi which thus occupies the first Nx*Ny elements of the state vector. We then have the values of \eta following.

    % vectors for coordinate selections 
NN=Nx*Ny;
l.left=[l.cbl; l.left; l.ctl];
l.right=[l.cbr; l.right; l.ctr];
l.phi=(1:NN)';
l.eta=((NN+1):(NN+Nx))';

% useful matrices
Z=zeros(NN,NN); ZZ=zeros(NN+Nx,NN+Nx);    zx=zeros(Nx,Nx);
I=eye(NN);      II=eye(NN+Nx);            ix=eye(Nx);

    System matrices

    The continuity equation for \phi is nothing else than its Laplacian equal zero \displaystyle \Delta \phi=0 the velocity field is recovered from \phi by derivations \displaystyle u=\phi_x, v=\phi_y

    This equation for \phi is not an evolution equation, there is no time derivative. The ingredient of evolution comes for the dynamics of the free surface. Just as for free_surface_gravity.m#boundary-conditions, the interface is locally advected vertically by v thus \displaystyle \eta_t=\phi_y(y=L)

    To build the dynamic matrix A we stack the Laplacian of the 2D grid and a matrix of zeros of the size of the \eta grid-that is-with Nx points. And then we add the block that connects the vertical velocity at the top to \eta. For this we use the location vectors. The state vector is thus \displaystyle q=\left(\begin{array}{c} \phi\\ \eta \end{array}\right)

    The left hand-side matrix E is full of zeros for the part of \phi (its equation has no time-derivative: it is just that \Delta \phi=0). The \eta part of E is the identity matrix.

    % system matrices
A=blkdiag(D.xx+D.yy,zx);
A(l.eta(2:end-1),l.phi)=D.y(l.top,:);
E=blkdiag(Z,ix);

    Boundary conditions

    We use no-penetration (slip) of the fluid (since it is not viscous, we cannot impose no-slip). This means that on the lateral boundaries we have u=\phi_x=0, this is the first equation in c

    Dx([left;right],:)*II(phi,:)

    At the right, there is the “big” identity matrix II for which we extract only the rows corresponding to \phi such that

    II(phi,:)*q=q(phi)

    and then we multiply this with the x differentiation matrix Dx but with just the rows that correspond to the left and right gridpoints. Thus this line of the constraint matrix tells that all this must be zero whatever happens.

    At the bottom boundary we have v=\phi_y=0, this is the second equation in c

    Dy(bot,:)*II(phi,:)

    We are now done with the boundary conditions for \phi and we can continue with the boundary conditions for \eta. We tell that it should have a zeros derivative at x=0 and x=Lx. This is good here because we can write by hand in #validation a theory for the validation of this case. This is writen in the third row of c where we use the 1D differentiation matrix in x dx, and selecting only the first and last rows of it

    dx([1,Nx],:)*II(eta,:)

    We then need to impose that the interface induces no variation of volume compared to the square domain-that is-its integral should be zero. This is just a scalar constraint. For this we use the integration vector ix built at the time of building the 1D matrices. This is the mast row of c

    d.wx*II(eta,:)

    Now we need yet one last constraint. We have writen in the dynamic matrix how the velocity field affects the interface (vertical advection), but we did not write how the interface affects the fluid. Across a free surface tensed by surface tension, there is the “Laplace pressure jump” proportional to the interface curvature times the value \sigma of the surface tension \displaystyle p_{in}-p_{out}=\sigma/R where R is the radius of curvature of the interface. We assume that the outside pressure is zero, and if the interface is close to flat, the expression of the curvature is simple (linear) \displaystyle p_{in}=-\sigma\eta_{xx}. If the second derivative of \eta is positive, the concave side of the interface is below, so the pressure should be less in the fluid than outside, this is why there is the minus sign.

    Now we need to relate \phi to the pressure inside the domain, for this we use the instationnary Bernoulli equation linearized about a zero base flow \displaystyle \rho \phi_t+p=0 thus we have \displaystyle \rho \phi_t=\sigma \eta_{xx} which is our last boundary condition. Note that this is not just a simple boundary conditions as we had until now, it is a little more subtle because there is a time derivative of \phi in this boundary condition. This means that in the numerical expression of this constraint, we will need to have a term in the left-hand-side of \displaystyle Eq_t=Aq so there will not only be zeros in E. This is coded in Ca and Ce. Ca is the part of the constraint that goes into A and Ce is the part of it that goes into E.

    % boundary conditions
loc=[l.top;l.left;l.bot;l.right;l.eta([1,2,Nx])];

c=[D.x([l.left;l.right],:)*II(l.phi,:); ...
   D.y(l.bot,:)*II(l.phi,:); ...
   d.x([1,Nx],:)*II(l.eta,:); ...
   d.wx*II(l.eta,:)];

Ca=[c; ...
    sigma*d.xx(2:end-1,:)*II(l.eta,:)];

Ce=[0*c; ...
    rho*II(l.top,:)]; 

A(loc,:)=Ca;
E(loc,:)=Ce;

    The structure of these matrices is quite sophisticated so it is interesting to see what they look like

    % show the matrices
figure(1); 
subplot(1,2,1); spy(E); title('E');
subplot(1,2,2); spy(A); title('A');
set(gcf,'paperpositionmode','auto');
print('-dpng','-r100','free_surface_2D_spy.png');
    spy of the dynamic matrices

    spy of the dynamic matrices

    Eigenmodes

    There is no dissipation and no instability in this flow, so all the modes should be neutral (the real part of the eigenvalues s should be 0), but there are oscillations so the imaginary parts should not be zero). The spectrum should be right-left symetric because the system itself is, so we remove the eigenmodes with negative imaginary parts.

    % compute eigenmodes
disp('computing eigenmodes');
[U,S]=eig(full(A),full(E));
s=diag(S);  [t,o]=sort(abs(s)); s=s(o); U=U(:,o);
rem=abs(s)>1000; s(rem)=[]; U(:,rem)=[];
rem=imag(s)<0; s(rem)=[]; U(:,rem)=[];

    Validation

    Just as for the example in wave_like.m#validation, we have chosen boundary conditions on \eta so that by chance our numerical solution should be the same as a wavelike solution periodic in x. We do here a theory similar to that of free_surface_gravity.m#validation but with surface tension (pressure jump through the interface) instead of gravity (pressure gradient in the fluid). And also as a difference we use velocity potential instead of stream-function (for a change…).

    We have the wavelike assumption \displaystyle \phi(x,y,t)=\hat{\phi}(y)\exp(i\alpha x+st)+C.C. so the Laplacian of \phi becomes \displaystyle -\alpha^2\hat{\phi}+\hat{\phi}_{yy}=0 and the solution should be like \displaystyle \hat{\phi}=A\cosh(\alpha y)+B\sinh(\alpha y) The bottom boundary condition is \hat{\phi}_y(y=0)=0 thus A=0.

    The advection of the interface gives \displaystyle s\hat{\eta}=\hat{\phi}(y=L)=\alpha B \sinh(\alpha L) and the presure jump through the interface gives \displaystyle s\hat{\phi}=-\alpha^2\sigma \hat{\eta} thus \displaystyle \hat{\eta}=-\frac{s}{\alpha^2\sigma}B\cosh(\alpha L) combining the two equations we get the dispertion relation for this system \displaystyle s^2=-\alpha^3\tanh(\alpha L)\sigma thus \displaystyle s=\pm i\sqrt{\sigma \alpha^3\tanh(\alpha L)} So we see that the eigenvalues are purely imaginary (neutral modes) and oscillate faster and faster as the wavenumber \alpha increases (the wavelength \lambda decreases). When the domain is very deep, L is large compared to 2\pi/\alpha so \alpha L becomes large, we have the limit of deep water waves \displaystyle s=\pm i \sqrt{\sigma \alpha ^3}

    To validate our computations, we plot the spectrum (the eigenvalues plotted in the complex plane) and we add this theory. For this we have to chose the values of \alpha of the wavelike theory that are relevant for our rectangular box. We say that thanks to the boundary conditions on \eta, there can be half a wavelength in Lx, one full wavelength and so on, that is \displaystyle \lambda=L_x/(k/2), k=1,2,3 ... We show this in the left plot of the figure below, showing a good agreement between the theory and the computation. To get even beter agreement, please increase the grid resolution (the computation is a little longer and there will be too much vectors on the velocity field plot).

    % validation
subplot(1,4,1)
lambda=Lx./([1:15]/2); 
alpha=2*pi./lambda;

stheo=i*sqrt(sigma*alpha.^3.*tanh(alpha*Ly))';
plot(real(s),imag(s),'b.',real(stheo),imag(stheo),'ro');
axis([-1,1,-10,100]);
xlabel('real part of eigenvalue'); ylabel('imaginary part of eigenvalue');
title('spectra'); legend('numeric','theory','location','north')
grid on

    Velocity fields

    It is nice to see what the flow looks like for the different modes of oscillation of our fluid/interface system. We need to get the velocity field and also the connection with the interface. We plot both the real part (in blue) and the imaginary part (in red).

    You can recognize directly by eye that these configurations corresponds to stable oscillations because at the top of the domain, the velocity field is in the direction opposite to the deflection of the free surface: the flow pushes back the surface toward rest (but it will overshoot and go to the other side, so the fluid will push back again and so on). If this was unstable, the flow would instead push the interface away from its rest position.

    % show velocity field and free surface
sel=[2,3,4];
for ind=1:length(sel);
    subplot(1,4,ind+1)
    
    % select eigenvector
    q=U(:,sel(ind)); 
    
    % extract u and v and reshape
    u=reshape(D.x*q(l.phi),Ny,Nx);
    v=reshape(D.y*q(l.phi),Ny,Nx);
    quiver(x,y,real(u),real(v),'b'); hold on 
    quiver(x,y,imag(u),imag(v),'r');
    
    % free surface
    e=q(l.eta);
    plot(x,real(e)+Ly,'b-',x,imag(e)+Ly,'r-')
    xlabel('x'); ylabel('y'); title(['mode ' num2str(sel(ind))]);
    axis equal; axis([0,Lx,0,2*Ly]);
    grid on


end
hold off

% print figure
set(gcf,'paperpositionmode','auto');
print('-dpng','-r100','free_surface_2D.png');
    The spectrum and the velocity field of the first three eigenmodes

    The spectrum and the velocity field of the first three eigenmodes

    You can see on this figure that we have a numerical eigenvalue s=0, why so? it is just because \phi comes into the equations only through its derivatives, so it is immaterial to the system if we add a constant to \phi. This means that if we build a state q where \phi is constant and \eta is 0, then \displaystyle Aq=0 so s=0 is the associated “eigenvalue”. To remove this uninteresting (so-called “spurious”) mode, we sould add a constraint on the state, just like we do in peristalsis.m or venturi.m for the pressure, saying that the value of \phi somewhere in the domain should be zero.

    Exercices/Contribution

    • Please change the boundary conditions on \eta to have fixed ends of the interface \eta(x=0)=0, \eta(x=Lx)=0 and compare how it changes the eigenmodes of oscillations.
    • Please do like in free_surface_gravity_particles.m, show the time evolution of the eigenmodes using the exponential in time. An then if you have the courage, advect also particles using this time dependent velocity field to give a realistic impression of the fluid/interface motion.
    • Please code all this using Navier-Stokes instead of velocity potential (for Masters! … master students?)
    • Please code gravity waves instead of capillary waves
    • Please using the gravity waves as above, look at the stable/unstable eigenmodes and re-evaluate what was said here about the flow pushing the interface back to rest (stable) or moving it away from rest (unstable).
    • Please add a constraint to remove the spurious mode with s=0