Oscillation modes of a partially filled cylindrical container

The aim of this tutorial program is to demonstrate how to use StabFem for capillary oscillations of static free-surface problem.

The investigated situation corresponds to a partially filled cylindrical container with radius R. we note L the height of liquid along the vertical wall (including the meniscus), gamma the surface tension, theta the contact angle and rhog the product between gravity acceleration and liquid density.

[[PUBLISH]] -> this is a tag for automatic publication as html

NB case updated sept 2020 using new solvers with lineic meshes... see old case in folder "OBSOLETE_CASES"

Contents

Initialization

close all;
addpath('../../SOURCES_MATLAB')
SF_Start('verbosity',4,'ffdatadir','./WORK/');
SF_core_setopt('eigensolver','SLEPC'); % this set of programs uses arpack ; to be rewritten someday with Slepc
SF_core_arborescence('cleanall');

system('mkdir FIGURES');
figureformat = 'png';

NOTICE  - Initialization already done
WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead

Chapter 1 : Oscillation modes for the case of a flat surface

We first create a mesh with rectangular shape (flat surface). See the Freefem++ program < https://gitlab.com/stabfem/StabFem/blob/master/STABLE_CASES/SLOSHING/MeshInit_Vessel.edp MeshInit_Vessel.edp > to see how it works.

% Dimensions of the initial mesh
R = 1; % radius of the vessel
L = 2; % Height of the vessel
density=50; % this is the mesh density in the initial mesh

% Physical parameters for this flat surface
gamma = 0.002;
rhog = 1;

freesurf = SF_Init('Mesh_Vessel.edp','Params',[L density],'problemtype','AxiFSStatic');
freesurf = SF_Deform(freesurf,'P',0,'gamma',gamma,'rhog',rhog,'typestart','pined','typeend','axis');

CHAPTER 2 : Eigenvalue computation for a flat surface, free condition, INVISCID

In this case there is a classical theoretical solution. Noting $m$ the azimuthal wavenumber and $n$ radial wavenumber (number of nodal circles), we have $\omega_{m,n} = \sqrt{ k_{m,n} g + \gamma / \rho k_{m,n}^3}$ with $k_{m,n} R = j'_{m,n}$ ($j'_{m,n}$ being the nth zero of the bessel funcion $J_m'$)

First we compute this theoretical solution for m=1 and n =(1,2,3) :

j1p(1) = 1.8412; j1p(2) = 5.3314; j1p(3) = 8.5363; % zeros of j1'
k = j1p(:)/R;
g= 1;gamma=0.002;rho=1;H=2;

disp( 'frequencies of m = 1 modes : ');
evTheory = sqrt( (g*k+gamma*k.^3).*tanh(k*H))

frequencies of m = 1 modes : 

evTheory =

    1.3606
    2.3737
    3.1274

We now compute the eigenvalues using the generic driver SF_Stability.m which invokes the FreeFem++ solver StabAxi_FS_Potential.edp

[evm1,emm1] =  SF_Stability(freesurf,'nev',10,'m',1,'shift',2.1i,'typestart','freeV','typeend','axis');
evm1FLATFREE = evm1;

NB alternative syntax

%[evm1FLATFREE,emm1] =  SF_Stability(freesurf,'solver','StabAxi_FS_Potential.edp','Options',{'nev',10,'m',1,'shift',2.1i,'typestart','freeV','typeend','axis'})
%evm1 = evm1FLATFREE

Note that we have used a shift centered on the (m,n) = (1,2) mode. with this choice the modes (m,n) = (1,1), (1,2) and (1,3) will be found with index (2,1,3)

Here we compare the results with the theory :

error1 = abs(imag(evm1FLATFREE(2))/evTheory(1)-1)+abs(imag(evm1FLATFREE(1))/evTheory(2)-1)+abs(imag(evm1FLATFREE(3))/evTheory(3)-1)

error1 =

   2.0417e-04

Chapter 2b : plot the results

Secondly plot the structure of the modes (we plot the imaginary part of the velocity potential along with the free surface deformation).

This is the first mode :

figure;SF_Plot(emm1(2),'phi.im','title',{'Mode (m,n)= (1,1)',['freq = ',num2str(imag(evm1(1)))]},'symmetry','YA','xlim',[-1 1]);
hold on;SF_Plot_ETA(emm1(2),'symmetry','YA');

ploting the first three modes

figure(2);
%suptitle('m=1 sloshing modes : Horizontal surface, free condition');% does not work with octave
hold on;
subplot(1,3,1);
SF_Plot(emm1(2),'phi.im','title',{'Mode (m,n)= (1,1)',['freq = ',num2str(imag(evm1(2)))]},'symmetry','YA','xlim',[-1 1]);
hold on;SF_Plot_ETA(emm1(2),'symmetry','YA');
subplot(1,3,2);
SF_Plot(emm1(1),'phi.im','title',{'Mode (m,n)= (1,2)',['freq = ',num2str(imag(evm1(1)))]},'symmetry','YA','xlim',[-1 1]);
hold on;SF_Plot_ETA(emm1(1),'symmetry','YA');
subplot(1,3,3);
SF_Plot(emm1(3),'phi.im','title',{'Mode (m,n)= (1,3)',['freq = ',num2str(imag(evm1(3)))]},'symmetry','YA','xlim',[-1 1]);
hold on;SF_Plot_ETA(emm1(3),'symmetry','YA');
pos = get(gcf,'Position'); pos(3)=pos(4)*2.6;set(gcf,'Position',pos); % resize aspect ratio

Note that the option ('symmetry','YA') allows to plot the structure in the symmetrized range r = [-1,1] taking into account the antisymmetry of phi (for m=1).

pause(1);

CHAPTER 3 : eigenvalues for flat surface, pined condition

disp('### Second test : flat surface, pined condition')

[evm1FLATPINED,emm1] = SF_Stability(freesurf,'solver','StabAxi_FS_Potential.edp','Options',{'nev',10,'m',1,'shift',2.1i,'typestart','pined','typeend','axis'});

%OLD SYNTAX : [evm1FLATPINED,emm1] =  SF_Stability(freesurf,'nev',10,'m',1,'shift',2.1i,'typestart','pined','typeend','axis');

evTheory = [1.4444    2.4786    3.2671]; % There is a paper by Miles with analytical solution...
error2 = abs(imag(evm1FLATPINED(2))/evTheory(1)-1)+abs(imag(evm1FLATPINED(1))/evTheory(2)-1)+abs(imag(evm1FLATPINED(3))/evTheory(3)-1)

### Second test : flat surface, pined condition

error2 =

    1.1336

The figures are very similar to the previous case, uncomment the next lines to see them !

%figure(1);
%plot(imag(evm1),0*real(evm1),'b+');hold on;
%legend('flat, pined');
%figure(3);suptitle('m=1 sloshing modes : Horizontal surface, H/R = 2, Bo = 500, pined condition');hold on;
%subplot(1,3,1);
%SF_Plot(emm1(2),'phi.im','title',{'Mode (m,n)= (1,1)',['freq = ',num2str(imag(evm1(2)))]},'symmetry','YA');hold on;SF_Plot_ETA(emm1(2),'symmetry','YA');
%subplot(1,3,2);
%SF_Plot(emm1(1),'phi.im','title',{'Mode (m,n)= (1,2)',['freq = ',num2str(imag(evm1(1)))]},'symmetry','YA');hold on;SF_Plot_ETA(emm1(1),'symmetry','YA');
%subplot(1,3,3);
%SF_Plot(emm1(3),'phi.im','title',{'Mode (m,n)= (1,3)',['freq = ',num2str(imag(evm1(3)))]},'symmetry','YA');hold on;SF_Plot_ETA(emm1(3),'symmetry','YA');
%pos = get(gcf,'Position'); pos(3)=pos(4)*2.6;set(gcf,'Position',pos); % resize aspect ratio
%pause(1);

Chapter 4 : meniscus (theta = 45?), free conditions

This case is the same as in Viola, Gallaire & Brun

We first compute equilibrium shape and corresponding mesh. The equilibrium shape is the solution of the Young-Laplace equation :

$p + \rho g z + \gamma K = 0$ along the free surface

We solve this equation using an iterative method implying a mesh deformation at each step. The method is described in this (unfinished and unpublished) paper : JFM_Menisci.pdf

A simpler matlab implementation of the same algorithm can be found here: meniscus.m

gamma = 0.002;
thetaE = pi/4; % this is the contact angle used in Viola et al.
hmeniscus = sqrt(2*gamma*(1-sin(thetaE))); % this is the expected height of the meniscus (valid for large Bond)
P = -freesurf.rhog*hmeniscus; % pressure in the liquid at z=0 (altitude of the contact line) ;
% the result will be to lower the interface by this ammount

Starting with the previously computed "flat surface" mesh as an initial "guess", we use the driver SF_Mesh_Deform.m to compute a mesh with a free surface satisfying the Young-Laplace solution. The driver launches the FreeFem++ solver Newton_Axi_FreeSurface_Static.edp .

freesurf = SF_Init('Mesh_Vessel.edp','Params',[L, density],'problemtype','AxiFSStatic'); % this line should not be necessary

freesurf = SF_Deform(freesurf,'P',P,'gamma',gamma,'rhog',rhog,'typestart','pined','typeend','axis');

We plot this mesh :

figure;SF_Plot(freesurf,'mesh','symmetry','ys');

Specific plots of the free surface geometry and curvature component

figure(12);plot(freesurf.mesh.meshlin.rsurf,freesurf.mesh.meshlin.zsurf,'-+');xlabel('r');ylabel('z');
figure(13);plot(freesurf.mesh.meshlin.rsurf,freesurf.mesh.meshlin.K0a,'-+b',freesurf.mesh.meshlin.rsurf,freesurf.mesh.meshlin.K0b,'-xr');
xlabel('r');ylabel('K0a,K0b');

A few statistics about the free surface shape :

zR = freesurf.mesh.meshlin.zsurf(1) % altitude of the contact line
z0 = freesurf.mesh.meshlin.zsurf(end) % altitude of the center
%Vol0 = freesurf.Vol % volume
%alphastart = freesurf.alpha(1)*180/pi % this should be 225 degrees (angle with respect to vertical = 45 degrees)

zR =

   8.0290e-06


z0 =

   -0.0342

We now compute m=1 eigenvalues using this mesh :

[evm1,emm1] =  SF_Stability(freesurf,'nev',10,'m',1,'shift',2.1i,'typestart','freeV','typeend','axis');
evm1MENISCUSFREE= evm1

evTheory = [1.3587    2.3630    3.1118]; % the values are compatible with results from viola et al.
error3 = abs(imag(evm1MENISCUSFREE(2))/evTheory(1)-1)+abs(imag(evm1MENISCUSFREE(1))/evTheory(2)-1)+abs(imag(evm1MENISCUSFREE(3))/evTheory(3)-1)

evm1MENISCUSFREE =

   0.0000 + 2.3563i
  -0.0000 + 1.3552i
  -0.0000 + 3.1005i
  -0.0000 + 3.8290i
  -0.0000 + 4.5932i
   0.0000 + 5.4141i
   0.0000 - 1.3552i
   0.0000 + 6.3011i
  -0.0000 - 2.3563i
  -0.0000 + 7.2587i


error3 =

    0.0091

We plot these modes

figure(4);
%suptitle('m=1 sloshing modes : Meniscus (45?), H/R = 2, Bo = 500, free condition');hold on;
subplot(1,3,1);
SF_Plot(emm1(2),'phi.im','title',{'Mode (m,n)= (1,1)',['freq = ',num2str(imag(evm1(2)))]},'symmetry','YA','xlim',[-1 1]);
hold on;SF_Plot_ETA(emm1(2),'Amp',0.05,'symmetry','YA');
subplot(1,3,2);
SF_Plot(emm1(1),'phi.im','title',{'Mode (m,n)= (1,2)',['freq = ',num2str(imag(evm1(1)))]},'symmetry','YA','xlim',[-1 1]);
hold on;SF_Plot_ETA(emm1(1),'Amp',0.05,'symmetry','YA');
subplot(1,3,3);
SF_Plot(emm1(3),'phi.im','title',{'Mode (m,n)= (1,3)',['freq = ',num2str(imag(evm1(3)))]},'symmetry','YA','xlim',[-1 1]);
hold on;SF_Plot_ETA(emm1(3),'Amp',0.05,'symmetry','YA');
pos = get(gcf,'Position'); pos(3)=pos(4)*2.6;set(gcf,'Position',pos); % resize aspect ratio

check if boundary condition is correctly verified

%figure(51);title('eta (plain), - d eta /ds + K0a cot(alpha) eta (dashed), limit (dotted)');
%DetaDs = diff(emm1(1).eta)./diff(freesurf.S0);
%plot(freesurf.xsurf,real(emm1(1).eta),'-'); hold on;
%plot((freesurf.xsurf(1:end-1)+freesurf.xsurf(2:end))/2,DetaDs,'--');
%plot(freesurf.xsurf,-freesurf.K0a.*cot(freesurf.alpha).*(abs(cot(freesurf.alpha))<1e2).*emm1(1).eta);

Chapter 5 : meniscus (theta = 45?), VISCOUS, m=1

NB in the viscous case, the driver SF_Stability.m invokes the FreeFem++ solver StabAxi_FreeSurface_Viscous.edp

nu = 1e-4;
m=1;sym =  'YA';% YS if m is even, YA if m is odd
[evm1VISCOUS,emm1] =  SF_Stability(freesurf,'nev',10,'m',1,'nu',nu,'shift',2.1i,'typestart','freeV','typeend','axis','plotspectrum',true);

At this stage it is useful to adapt the mesh to capture the thin boundary layer

 freesurf = SF_Adapt(freesurf,emm1(3),'Hmax',.1) ;
 figure(7);SF_Plot(freesurf,'mesh');

Now recompute the modes

[evm1V,emm1V] =  SF_Stability(freesurf,'nev',10,'m',1,'nu',nu,'shift',2.1i,'typestart','freeV','typeend','axis');

WARNING -  baseflow.iter = 0 ! it seems your baseflow was projected after mesh adaptation but not recomputed !

Plot the results

figure(5);sym='YA';symUZ ='YS';hold off;
%suptitle(['Sloshing modes : Meniscus (45?), H/R = 2, Bo = ' num2str(1/gamma) '; Oh = ' num2str(nu)  '; m = ' num2str(m) ]);hold on;
subplot(1,3,1);
SF_Plot(emm1(2),'uz.im','title',{'Mode (m,n)= (1,1)',['\omega_r = ',num2str(imag(evm1V(1))),', \omega_i = ',num2str(real(evm1V(1)))]},'symmetry',sym,'xlim',[-1 1]);hold on;
SF_Plot_ETA(emm1V(2),'Amp',0.05,'symmetry',sym);xlim([-1 1]);ylim([-2,.5]);
subplot(1,3,2);
SF_Plot(emm1(1),'uz.im','title',{'Mode (m,n)= (1,2)',['\omega_r = ',num2str(imag(evm1V(2))),', \omega_i = ',num2str(real(evm1V(2)))]},'symmetry',sym,'xlim',[-1 1]);hold on;
SF_Plot_ETA(emm1V(1),'Amp',0.05,'symmetry',sym);xlim([-1 1]);ylim([-2,.5]);
subplot(1,3,3);
SF_Plot(emm1V(3),'uz.im','title',{'Mode (m,n)= (1,3)',['\omega_r = ',num2str(imag(evm1V(3))),', \omega_i = ',num2str(real(evm1V(3)))]},'symmetry',sym,'xlim',[-1 1]);hold on;
SF_Plot_ETA(emm1V(3),'Amp',0.05,'symmetry',sym);xlim([-1 1]);ylim([-2,.5]);
pos = get(gcf,'Position'); pos(3)=pos(4)*2.6;set(gcf,'Position',pos); % resize aspect ratio

figure(6);
SF_Plot(emm1V(2),'uz.im','title',{'Mode (m,n)= (1,1)',['\omega_r = ',num2str(imag(evm1V(1))),', \omega_i = ',num2str(real(evm1V(1)))]},'symmetry',sym,'xlim',[-1 1]);hold on;
SF_Plot_ETA(emm1V(2),'Amp',0.05,'symmetry',sym);xlim([-1 1]);ylim([-2,.5]);
title('Sloshing mode in a cylindrical container (viscous problem)')




figure(7);
plot(imag(evm1FLATFREE),real(evm1FLATFREE),'g+');hold on;
plot(imag(evm1FLATPINED),real(evm1FLATPINED),'ro');hold on;
plot(imag(evm1MENISCUSFREE),real(evm1MENISCUSFREE),'bx');hold on;
plot(imag(evm1VISCOUS),real(evm1VISCOUS),'k*');hold on;
xlabel('omega_r');ylabel('omega_i');
legend('flat, free', 'flat, pined','meniscus, free','meniscus, pined, viscous' );
xlim([0 6]);ylim([-.1 .05]);

APPENDIX B : meniscus (theta = 45?), free conditions m=0

(uncomment the following to see how it works...)

%his case is the same as in Viola, Gallaire & Brun
%
% nu = 1e-3;
% m=0;sym =  'YS';% YS if m is even, YA if m is odd
% [evm1,emm1] =  SF_Stability(freesurf,'nev',10,'m',m,'nu',nu,'shift',2.1i,'typestart','freeV','typeend','axis');
%
%
% %evTheory = [1.3587    2.3630    3.1118];
% %error4 = abs(imag(evm1(2))/evTheory(1)-1)+abs(imag(evm1(1))/evTheory(2)-1)+abs(imag(evm1(3))/evTheory(3)-1)
%
% figure(6);hold off;
% suptitle(['Sloshing modes : Meniscus (45?), H/R = 2, Bo = ' num2str(1/gamma) '; Oh = ' num2str(nu)  '; m = ' num2str(m) ]);hold on;
% subplot(1,3,1);
% SF_Plot(emm1(1),'uz1.im','title',{'Mode (m,n)= (0,1)',['\omega_r = ',num2str(imag(evm1(1))),', \omega_i = ',num2str(real(evm1(1)))]},'symmetry',sym);
% hold on;SF_Plot_ETA(emm1(1),'Amp',0.05,'symmetry','YS');xlim([-1 1]);ylim([-2,.5]);
% subplot(1,3,2);
% SF_Plot(emm1(2),'uz1.im','title',{'Mode (m,n)= (0,2)',['\omega_r = ',num2str(imag(evm1(2))),', \omega_i = ',num2str(real(evm1(2)))]},'symmetry',sym);
% hold on;SF_Plot_ETA(emm1(2),'Amp',0.05,'symmetry','YS');xlim([-1 1]);ylim([-2,.5]);
% subplot(1,3,3);
% SF_Plot(emm1(3),'uz1.im','title',{'Mode (m,n)= (0,3)',['\omega_r = ',num2str(imag(evm1(3))),', \omega_i = ',num2str(real(evm1(3)))]},'symmetry',sym);
% hold on;SF_Plot_ETA(emm1(3),'Amp',0.05,'symmetry','YS');xlim([-1 1]);ylim([-2,.5]);
% pos = get(gcf,'Position'); pos(3)=pos(4)*2.6;set(gcf,'Position',pos); % resize aspect ratio
%
%
%
%
%

APPENDIX C : meniscus (theta = 45?), free conditions m=2

% This case is the same as in Viola, Gallaire & Brun
% nu = 1e-3;
% [evm1,emm1] =  SF_Stability(freesurf,'nev',10,'m',2,'nu',nu,'shift',2.1i,'typestart','freeV','typeend','axis');
% m=2;sym =  'YS';% YS if m is even, YA if m is odd
%
% %evTheory = [1.3587    2.3630    3.1118];
% %error4 = abs(imag(evm1(2))/evTheory(1)-1)+abs(imag(evm1(1))/evTheory(2)-1)+abs(imag(evm1(3))/evTheory(3)-1)
%
% figure(7);hold off;
% suptitle(['Sloshing modes : Meniscus (45?), H/R = 2, Bo = ' num2str(1/gamma) '; Oh = ' num2str(nu)  '; m = ' num2str(m) ]);hold on;
% subplot(1,3,1);
% SF_Plot(emm1(1),'uz1.im','title',{'Mode (m,n)= (2,1)',['\omega_r = ',num2str(imag(evm1(1))),', \omega_i = ',num2str(real(evm1(1)))]},'symmetry',sym);
% hold on;SF_Plot_ETA(emm1(1),'Amp',0.05,'symmetry','YS');xlim([-1 1]);ylim([-2,.5]);
% subplot(1,3,2);
% SF_Plot(emm1(2),'uz1.im','title',{'Mode (m,n)= (2,2)',['\omega_r = ',num2str(imag(evm1(2))),', \omega_i = ',num2str(real(evm1(2)))]},'symmetry',sym);
% hold on;SF_Plot_ETA(emm1(2),'Amp',0.05,'symmetry','YS');xlim([-1 1]);ylim([-2,.5]);
% subplot(1,3,3);
% SF_Plot(emm1(3),'uz1.im','title',{'Mode (m,n)= (2,3)',['\omega_r = ',num2str(imag(evm1(3))),', \omega_i = ',num2str(real(evm1(3)))]},'symmetry',sym);
% hold on;SF_Plot_ETA(emm1(3),'Amp',0.05,'symmetry','YS');xlim([-1 1]);ylim([-2,.5]);
% pos = get(gcf,'Position'); pos(3)=pos(4)*2.6;set(gcf,'Position',pos); % resize aspect ratio










%
fclose(fopen('SCRIPT_sloshing.success','w+'));


Published with MATLAB® R2020a