sandbox/easystab/free_surface_gravity_particles.m

    Free-surface: displaying the eigenmodes using particles

    In this code, we take the code that computes the eigenmodes for free surface waves free_surface_gravity.m and we visualize the flow motion of the first eigenmode by expanding the solution to physical space and advecting tracer particles, like we did in a simple case in particles.m.

    clear all; clf;
    
    n=100;      % number of gridpoints
    alpha=1;    % wavenumber in x
    L=2;        % Fluid height in y
    rho=1;      % fluid density
    mu=0.0001;    % fuid viscosity
    g=1;        % gravity
    
    % differentiation matrices
    scale=-2/L;
    [y,DM] = chebdif(n,2); 
    D=DM(:,:,1)*scale;    
    DD=DM(:,:,2)*scale^2;    
    y=(y-1)/scale; 
    I=eye(n); Z=zeros(n,n);
    
    % renaming the matrices
    dy=D; dyy=DD;
    dx=i*alpha*I; dxx=-alpha^2*I;
    Delta=dxx+dyy;
    
    % System matrices
    A=[mu*Delta, Z, -dx, Z(:,1); ...
       Z, mu*Delta, -dy, Z(:,1); ...
       dx, dy, Z, Z(:,1); ...
       Z(1,:),I(n,:),Z(1,:),0];
    
    E=blkdiag(rho*I,rho*I,Z,1);
    
    % boundary conditions
    loc=[1,n,n+1,2*n];
    C=[I(1,:),Z(1,:),Z(1,:),0; ... 
       Z(1,:),I(1,:),Z(1,:),0; ...
       Z(1,:),Z(1,:),-I(n,:),rho*g; ...
       dy(n,:),dx(n,:),Z(1,:),0]; 
    
    E(loc,:)=0;  
    A(loc,:)=C;
    
    % compute eigenmodes
    [U,S]=eig(A,E);
    s=diag(S);  [t,o]=sort(-real(s)); s=s(o); U=U(:,o);
    rem=abs(s)>1000; s(rem)=[]; U(:,rem)=[];

    The animation

    We show the position of the free surface stored in q(eta), and the particles are initialy set as a mesh, and stretched in y to fit with the initial position of the free surface for the top prticles.

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % particle animation of the eigenmodes
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % parameters
    nper=1      % number of periods of oscillation
    nt=30       % number of time steps per period
    nx=40       % number of points in x
    modesel=2;  % which more do animate
    
    % select the eigenmode
    u=1:n; v=u+n; p=v+n; eta=3*n+1;
    q=U(:,modesel); 
    lambda=s(modesel);
    
    % The time and x extent
    tvec=linspace(0,nper*2*pi/abs(lambda),nper*nt); 
    dt=tvec(2)-tvec(1);
    Lx=2*pi/alpha; 
    x=linspace(0,Lx,nx);
    
    % scale mode amplitude 
    q=0.05*q/abs(q(eta)); 
    
    % initialize tracer particles
    [px,py]=meshgrid(linspace(0,Lx,60),linspace(0,L,30)); 
    py=py.*(1+2*real(exp(lambda*tvec(1))*q(eta)*exp(i*alpha*px))/L);
    px=px(:);py=py(:);  
    
    % time loop
    for ind=1:nper*nt
    

    Expand the solution to physical space

    We have assumed a wavelike behaviour of the flow solution in the x direction, thus we have writen for all the flow variables, for instance for the u velocity \displaystyle u(x,y,t)=\hat{u(y)}\exp(\alpha x+s t)+CC where CC is the complex conjugate of the other term. This way we have transformed a variable u depending on three coordinates, to a variable \hat{u} depending on one single coordinate, thus transforming everything into a 1D problem. What we do now is to transform back the 1D problem into a full original problem by using the exponential.

    A complex quantity plus its complex conjugate gives twice the real part of it, this is what we code here.

    In the code, we do the expansion to hysical space for all the variables (u,v,p,\eta) at once in a single command.

        % expand mode to physical space 
        qq=2*real(exp(lambda*tvec(ind))*q*exp(i*alpha*x));
            
        % plot pressure
        surf(x,y,qq(p,:)-10,'facealpha',0.3); view(2); shading interp; hold on
        
        % plot free surface
        plot(x,L+qq(eta,:),'k-',x,0*x+L,'k--'); hold on
        
        %%%% plot the particles
        plot(mod(px,Lx),py,'k.');
    

    Interpolation of the particle positions

    The particles are not on a cartesian grid, they can be wherever the field has brought them, the only thing that matter here for the interpolation is where they are on the y grid, and then we can do the expansion in x as a second step, this is less computations.

        % compute velocity at position of particles
        pu=interp1(y,q(u),py);
        pv=interp1(y,q(v),py);
    

    The Taylor expansion

    When we wrote the boundary conditionat the free surface for setting up the linear model, we have used Taylor expansions to express the velocity and pressure at the position of the moving interface as a function of the values at the fixed position L. here we do again the same thing for the particles that are above L.

        % For particles above L, use Taylor expansion for velocity
        pu(py>L)=q(u(n))+q(eta)*dy(n,:)*q(u);
        pv(py>L)=q(v(n))+q(eta)*dy(n,:)*q(v);
        
        % expand to physical space
        puu=2*real(exp(lambda*tvec(ind))*pu.*exp(i*alpha*px));
        pvv=2*real(exp(lambda*tvec(ind))*pv.*exp(i*alpha*px));
     
        % advect particles
        px=px+puu*dt;
        py=py+pvv*dt;
        
        xlabel('x');    
        ylabel('y'); 
        axis equal; axis([0,Lx,0,1.3*L]); 
        grid off
        hold off
        drawnow
    end
    
    set(gcf,'paperpositionmode','auto');
    print('-dpng','-r100','free_surface_gravity_particles.png');
    
    The figure at final iteration of the time loop, shwing the free surface at the top, the particles in black, and the disturbance pressure field in color

    The figure at final iteration of the time loop, shwing the free surface at the top, the particles in black, and the disturbance pressure field in color

    Exercices/Contributions

    • Please make a movie instead of just an image to display on the page –> free_surface_gravity_particles_gif.m
    • Please play with the parameters of the particles, and the amplitude of the eigenmode
    • Please check that the first and second eigenmodes correspond to waves propagating in opposite directions –>free_surface_gravity_first_modes.m
    • Please use this animation to check the theory for the wave speed of gravity waves
    • Please plot also the vector field using quiver in a way that confirms the motion of the particles that we observe –>free_surface_gravity_vector_field.m
    • Please investigate the motion of the other eigenmodes (which correspond to damped stationnary waves, for now, there is a bug that makes things strange for the non propagative modes…)
    • Please consider the unstable case with g negative and find the nice scaling of the mode amplitude and the time of the animation so that we nicely see the flow in this instability
    • Please investigate how the Taylor expansino is not accurate when the surface deformation is large (you can see this with the particles that are close to the interface)
    • Please do the animation showing the total pressure instead of only the pressure disturbance (that is, add the hydrostatic pressure field as well to see what this changes)free_surface_gravity_total_pressure
    • Please do the same visualisation with particle for the other flow instabilities