Posted in Uncategorized | Leave a comment

Big_Case1Liquid-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.

Posted in Uncategorized | Leave a comment


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

Posted in Uncategorized | Leave a comment

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.

Posted in Uncategorized | Leave a comment

Video gallery: rotating plumes

Posted in Uncategorized | Leave a comment

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).

Posted in Uncategorized | Leave a comment

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.


Posted in Uncategorized | Leave a comment

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/Nt_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.




Posted in Uncategorized | Leave a comment


Tips and Tricks page has been updated!

Posted in Uncategorized | Leave a comment