# Title: Shape oscillation of an inviscid droplet # # Description: # # A two-dimensional elliptical droplet (density ratio 1/1000) is # released in a large domain. Under the effect of surface-tension the # shape of the droplet oscillates around its (circular) equilibrium # shape. The fluids inside and outside the droplet are inviscid so # ideally no damping of the oscillations should occur. As illustrated # on Figure \ref{kinetic} some damping occurs in the simulation due to # numerical dissipation. # # This simulation is also a stringent test case of the accuracy of the # surface tension representation as no explicit viscosity can damp # eventual parasitic currents. # # \begin{figure}[htbp] # \caption{\label{kinetic}Evolution of the kinetic energy as a function # of time for the spatial resolutions indicated in the legend. The # black lines are fitted decreasing exponential functions.} # \begin{center} # \includegraphics[width=0.9\hsize]{k.eps} # \end{center} # \end{figure} # # The initial shape of the droplet is given by: # $$r(\theta) = r_0 + \alpha\cos(n\theta)$$ # The oscillation frequency is then \cite{torres00}: # $$\omega_n^2={(n^3-n)\sigma\over (\rho_d+\rho_e)r_0^3}$$ # # A comparison between the theoretical and numerical values of the # frequency is given in Figure \ref{frequency}. # # \begin{figure}[htbp] # \caption{\label{frequency}Relative error in the oscillation # frequency as a function of resolution.} # \begin{center} # \includegraphics[width=0.8\hsize]{frequency.eps} # \end{center} # \end{figure} # # The amount of numerical damping can be estimated by computing an # equivalent viscosity. With viscosity, kinetic energy is expected to # decrease as: # $$\exp(-C\nu/D^2t)$$ # where $C$ is a constant, $\nu$ the viscosity and $D$ the droplet # diameter. Using curve fitting the damping coefficient $b=C\nu/D^2$ # can be estimated (black curves on Figure \ref{kinetic}). An # equivalent Laplace number can then be computed as: # $$La={\sigma D\over \rho\nu^2}={\sigma C^2 \over \rho b^2 D^3}$$ # The equivalent Laplace number depends on spatial resolution as # illustrated in Figure \ref{laplace}. # # \begin{figure}[htbp] # \caption{\label{laplace}Equivalent Laplace number estimated from the # numerical damping of kinetic energy.} # \begin{center} # \includegraphics[width=0.8\hsize]{laplace.eps} # \end{center} # \end{figure} # # Author: St\'ephane Popinet # Command: sh oscillation.sh # Version: 1.2.0 # Required files: oscillation.sh fit.ref # Running time: 3 minutes # Generated files: frequency.eps k.eps laplace.eps # Define DIAMETER 0.2 Define EPSILON 0.05 Define VAR(T,min,max) (min + CLAMP(T,0,1)*(max - min)) Define RHO(T) VAR(T, 1., 1e-3) Define RADIUS(x,y) (DIAMETER/2.*(1. + EPSILON*cos (2.*atan2 (y, x)))) 1 0 GfsSimulation GfsBox GfsGEdge {} { Time { end = 1 } Refine LEVEL VariableTracerVOF T VariableFiltered T1 T 1 InitFraction T ({ x += 0.5; y += 0.5; return x*x + y*y - RADIUS(x,y)*RADIUS(x,y); }) PhysicalParams { alpha = 1./RHO(T1) } VariableCurvature K T SourceTension T 1. K # fixme: A small amount of viscosity seems to be necessary to # obtain a non-increasing kinetic energy at high-resolution (8 # levels). SourceViscosityExplicit 7e-5*RHO(T1) AdaptFunction { istep = 1 } { cmax = 0.01 maxlevel = LEVEL } (T > 0 && T < 1 ? 1. : fabs (Vorticity)*ftt_cell_size (cell)) OutputScalarSum { istep = 1 } k-LEVEL { v = RHO(T1)*Velocity2 } EventScript { start = end } { cat <&1 | awk '{if ($1 == "result:") print LEVEL,$2,$3,$4;}' k(t)=a*exp(-b*t)*(cos(c*t+3.14159265359)+1.) a = 3e-4 b = 1.5 c = 153 fit k(x) 'k-LEVEL' u 3:5 via a,b,c print "result: ", a, b, c EOF rm -f fit.log } } GfsBox { left = Boundary bottom = Boundary }