Liquid-liquid breakup simulation in a turbulent flow. Interface in yellow. Slices show the  vorticity magnitude (in color), the velocity magnitude (in contours) and the grid.

Gerris3D

Open source flow solver with oct-tree grid refinement used to study oil jets in water environments.Above: Lambda2=-0.005 isosurfaces colored by velocity magnitude. Above: Oil interface colored by streamwise velocity.

Spectral element method: a simple example

Consider the 1D advection problem of a scalar $\phi$ with constant velocity $U=1$:

$\frac{\partial \phi}{\partial t} = - U \frac{\partial \phi}{\partial x}, \qquad -1 \le x \le 1, \qquad t > 0$,

$\phi(x,0) = e^{-ln(2)(x+1)^2/\sigma^2}$,

$\phi(-1,t) = e^{-ln(2)(1-t)^2/\sigma^2}$,

where $\sigma=0.2$.

Fig.1 compares the convergence (norm 2 of the difference between the numerical and exact solution) of Discontinuous Galerkin Method and Finite Difference for the problem above.

Fig.1: Convergence of Discontinuos Galerkin Method (DG) and Finite difference (FD) for a 1D advection problem.

Image gallery: rotating plumes

Figure 1. From left to right: temporally and azimuthally averaged passive scalar volume fraction field for plumes at Ro=$\infty$, 40, 10 and 1. These results show how increasing the rotation rate from Ro=$\infty$ to 10, results in the plume axis deflection and reduces the lateral intrusion height. At lower Rossby numbers, the lateral transport is completely suppressed leading to a change of regime (see Ro=1 case) characterized by the accumulation of passive scalar in the domain. These results are in agreement with the experiments of Helfrich and Battisti [1].

Figure 2. Similar to above. The mean circulation is visualized using streamlines colored by azimuthal velocity magnitude. The plume axis deflection accompanied by the blocking of the radial transport increases and the system rotation rate grows. The complete suppression of the radial transport is evident at Ro=1 characterized by the lack of lateral intrusion signature and the presence of smaller eddies that explain the well mixed passive scalar field shown in Fig.1.

[1] Helfrich, K. R. and Battisti, T. M., “Experiments on baroclinic vortex shedding from hydrothermal plumes“, Journal of Geophysical Research, volume 96, pages 12511-12518 (1991).

Image gallery: hybrid plumes

Figure 1. Partially averaged velocity field for a single phase plume. Left: two isosurfaces of radial velocity at $\langle u_r \rangle=0.1$ (red) and $\langle u_r \rangle = -0.07$ (blue). Note the entrainment associated to negative values of $\langle u_r \rangle$ along the plume stem below the trapping height and the formation of the lateral intrusion associated to positive values of radial velocity. Right: two isosurfaces of vertical velocity at $\langle w \rangle=0.1$ (red) and $\langle w \rangle=-0.15$. Note the plume top surrounded by a ring of downdrafts (in blue) and the formation of gravity waves (outer ring of positive $\langle w \rangle$).

Figure 2. Similar to above for a bubble plume. The presence of gas, characterized by the slip velocity, leads to a destratification across the entire water column and the formation of peeling regions above the main lateral intrusion.

Numerical model for bubble plumes

Here we briefly summarize the derivation of our Eulerian-Eulerian model for multiphase (gas and water) flows used to investigate plumes generated by subsurface blowouts in stratified environments. The buoyancy defect at the source is assumed to be due to thermal effects and the presence of a gas phase. The main assumptions are:

1. Boussinesq flow conditions (density variations only accounted in the buoyancy term).
2. Gas bubbles carry a negligible amount of momentum compared to the liquid phase.
3. Constant slip velocity of the gas phase.
4. Linear background stratification.

The starting point for the model derivation, is the mass conservation equation for a phase $k$:

$\frac{\partial}{\partial \tilde{t}} \left( \alpha_k \rho_k \right)+ \nabla \cdot \left( \alpha_k \rho_k \mathbf{\tilde{u}}_k \right) =0$

where $\rho_k$, $\alpha_k$ and $\mathbf{\tilde{u}}_k=\left( \tilde{u},\tilde{v},\tilde{w} \right)_k$ are density, volume fraction and Cartesian velocity field for the phase $k$ respectively, $\tilde{t}$ is time and the symbol $\tilde{\phantom{.}}$ denotes dimensional quantities.

Similarly, the momentum conservation can be written as:

$\frac{D}{D \tilde{t}}\left( \alpha_k \rho_k \mathbf{\tilde{u}}_k \right)=-\alpha_k\nabla \tilde{p}+\nabla \cdot \alpha_k \widetilde{\boldsymbol{\tau}}_k+\alpha_k \rho_k g \vec{k} - 2\mathbf{\Omega} \times \mathbf{\tilde{u}}_k+\widetilde{\mathbf{M}}_k$

where $\tilde{p}$ is pressure (taken to be equal in each phase [1]), $\nabla \cdot \alpha_k \widetilde{\boldsymbol{\tau}}_k$ are the molecular and turbulent diffusive terms for phase $k$, $g=-9.8 \, m/s^2$ is the gravity acceleration acting in the vertical direction aligned with coordinate $z$, $\mathbf{\Omega}$ is the rotation vector and $D/Dt=\partial/\partial t +\mathbf{\tilde{u}}_k \cdot \nabla$ is the material derivative.
Momentum transfer between phases is accounted in the $\widetilde{\mathbf{M}}_k$ term and we require:

$\sum_k \widetilde{\mathbf{M}}_k=0,$

i.e. phase interactions do not affect the bulk momentum.

As detailed in [1], under the assumptions presented above, the system of transport equations can be simplified as:

$\nabla \cdot \mathbf{u}=0$

$\frac{D \mathbf{u}}{D t} =-\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u} + \left(\theta + Ri \alpha_b\right) \hat{\textbf{k}} - \frac{1}{Ro} \left(\hat{\textbf{k}} \times \mathbf{u}\right)$

$\frac{D \theta}{D t} = \frac{1}{Pe_T} \nabla^2 \theta - \mathbf{u} \cdot \hat{\textbf{k}}$

$\frac{D \alpha_b}{D t} = \frac{1}{Pe_b} \nabla^2 \alpha_b - U_N \frac{\partial \alpha_b}{\partial z}$

$\frac{D \beta}{D t} = \frac{1}{Pe_p} \nabla^2 \beta$

where $\mathbf{u}=(u,v,w)$ is the liquid phase velocity, $\alpha_b$ is the gas volume fraction, $\beta$ is a passive scalar volume fraction, $\theta = T - \zeta z$ is the temperature $T$ perturbation and $\zeta$ is the stratification slope. The system of equations has been non-dimensionalized using the classical velocity, length, time and temperature scales for stratified flows.
Using the inlet buoyancy flux $B_0$ and the buoyancy frequency $N$, these scales are defined as $U_S = (B_0 \, N)^{1/4}$$L_S = U_S/N$$t_S = 1/N$ and $T_S = \zeta L_S$ respectively. $U_N = w_s / U_S$ is the non-dimensional slip velocity and $N$ is the buoyancy frequency defined as:

$N=\sqrt{-\frac{g}{\rho_0}\frac{\partial \rho_e}{\partial z}} = \sqrt{g \gamma \zeta}$

where $\rho_0$ is the reference density value (taken at $z=0$ for an unperturbed environment), $\rho_e = \rho_0 (1-\gamma \zeta z)$ is the environment density profile and $\gamma$ is the thermal expansion coefficient for the liquid phase.

The non-dimensional groups are the Reynolds number $Re=\frac{U_S L_S}{\nu}$, the Richardson number $Ri=\frac{g \, L_S}{U_S^2}$, the Rossby number $Ro=\frac{U_S}{f \, L_S}$ and the Péclet number for a scalar $m$ $Pe_m =\frac{U_s L_S}{\mathcal{D}_m}$. $\nu$ is the liquid phase kinematic viscosity, $\mathcal{D}_m$ is the diffusivity coefficient for the scalar $m$ and $f=2\Omega$ is the Coriolis frequency (the rotation vector $\mathbf{\Omega}$ is here aligned with $z$).

Detailed derivation of the model can be found in:

[1] Fabregat Tomàs, A., Dewar, W. K., Özgökmen, T. M., Poje, A. C. and Wienders, N., “Numerical Simulations of Turbulent Thermal, Bubble and Hybrid Plumes”, Ocean Modelling, Volume 90, pages 16-28, 2015.