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
- Chapter 1 : Oscillation modes for the case of a flat surface
- CHAPTER 2 : Eigenvalue computation for a flat surface, free condition, INVISCID
- Chapter 2b : plot the results
- CHAPTER 3 : eigenvalues for flat surface, pined condition
- Chapter 4 : meniscus (theta = 45?), free conditions
- Chapter 5 : meniscus (theta = 45?), VISCOUS, m=1
- APPENDIX B : meniscus (theta = 45?), free conditions m=0
- APPENDIX C : meniscus (theta = 45?), free conditions m=2
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 the azimuthal wavenumber and
radial wavenumber (number of nodal circles), we have
with
(
being the nth zero of the bessel funcion
)
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 :
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+'));