#include "all-mach.h"
#include "vof.h"

Include dependency graph for two-phase.h:

This graph shows which files directly or indirectly include this file:

Macros
#define	NO_1D_COMPRESSION 1
	We transport VOF tracers without the one-dimensional compressive term.

#define	mu(f) (mu1mu2/(clamp(f,0,1)(mu2 - mu1) + mu1))
	By default the Harmonic mean is used to compute the phase-averaged dynamic viscosity.

#define	lambdav(f) (clamp(f,0,1)*(lambdav1 - lambdav2) + lambdav2)
	The volumetric viscosity uses arithmetic average by default.

Functions
double	sound_speed (Point point)
	These functions are provided by the Equation Of State.

double	average_pressure (Point point)

double	bulk_compressibility (Point point)

double	internal_energy (Point point, double fc)

void	event_stability (void)
	Event: stability (i++)

void	fE_refine (Point point, scalar fE)

void	event_defaults (void)
	Event: defaults (i = 0)

void	event_init (void)
	Event: init (i = 0)

void	event_vof (void)
	Event: vof (i++)

void	event_tracer_advection (void)
	We set the list of interfaces back to its default value.

void	event_properties (void)
	Event: properties (i++)

void	event_end_timestep (void)
	Event: end_timestep (i++)

void	event_cleanup (void)
	At the end of the simulation we clean the tracer list.

Variables
scalar	f []
	The primary fields are:

scalar *	interfaces = {f}
	The height functions are stored in the vector field associated with each VOF tracer.

scalar	frho1 []

scalar	frho2 []

scalar	fE1 []

scalar	fE2 []

double	CFLac = HUGE
	The timestep can be limited by a CFL based on the speed of sound.

double	mu1 = 0.
	The dynamic viscosities for each phase, as well as the volumetric viscosity coefficients.

double	mu2 = 0.

double	lambdav1 = 0.

double	lambdav2 = 0.

scalar	rhoc2v []

vector	alphav []

scalar	rhov []

const vector	lambdav0 [] = {0,0,0}

const vector	lambdav = lambdav0

static scalar *	interfaces1 = NULL

Macro Definition Documentation

◆ lambdav

#define lambdav ( f ) (clamp(f,0,1)*(lambdav1 - lambdav2) + lambdav2)

The volumetric viscosity uses arithmetic average by default.

Definition at line 94 of file two-phase.h.

◆ mu

#define mu ( f ) (mu1*mu2/(clamp(f,0,1)*(mu2 - mu1) + mu1))

By default the Harmonic mean is used to compute the phase-averaged dynamic viscosity.

Definition at line 85 of file two-phase.h.

◆ NO_1D_COMPRESSION

#define NO_1D_COMPRESSION 1

We transport VOF tracers without the one-dimensional compressive term.

A compressible two-phase flow solver

This solves the two-phase compressible Navier-Stokes equations including the total energy equation.

\[ \frac{\partial (f \rho_i)}{\partial t } + \nabla \cdot (f \rho_i \mathbf{u}) = 0 \]

\[ \frac{\partial (\rho_i \mathbf{u})}{\partial t } + \nabla \cdot ( \rho_i \mathbf{u} \mathbf{u}) = -\nabla p + \nabla \cdot \tau_i' \]

\[ \frac{\partial [\rho_i (e_i + \mathbf{u}^2/2)]}{\partial t } + \nabla \cdot [ \rho_i \mathbf{u} (e_i + \mathbf{u}^2/2)] = -\nabla \cdot (\mathbf{u} p_i) + \nabla \cdot \left( \tau'_i \mathbf{u} \right) \]

an advection equation for the volume fraction \(f\)

\[ \frac{\partial f}{\partial t} + \mathbf{u} \cdot \nabla f = 0 \]

and an Equation Of State (EOS) that needs to be defined.

The method is described in detail in Fuster & Popinet, 2018 and relies on a time splitting technique were we first solve the advection step using the VOF method for \(f\), \(f_i \rho_i\), \(f_i \rho_i \mathbf{u}\) \(f_i \rho_i e_{tot,i}\) and then perform the projection step using the all-Mach solver.

Definition at line 40 of file two-phase.h.

Function Documentation

◆ average_pressure()

double average_pressure ( Point point )

extern

Definition at line 59 of file Mie-Gruneisen.h.

References clamp(), dimension, f, fE1, fE2, frho1, frho2, PIGAMMA, q, sq(), vector::x, and x.

Referenced by event_properties().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ bulk_compressibility()

double bulk_compressibility ( Point point )

extern

Bulk compressibility of the mixture

i.e. \(\rho c^2\).

Bulk compressibility of the mixture

i.e. \(\rho c^2\).

Definition at line 75 of file Mie-Gruneisen.h.

References b1, b2, clamp(), f, frho1, frho2, gamma1, gamma2, p, PI1, PI2, PIGAMMA, and x.

Referenced by event_properties().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ event_cleanup()

void event_cleanup ( void )

At the end of the simulation we clean the tracer list.

Event: cleanup (i = end)

Definition at line 554 of file two-phase.h.

References f, free(), and x.

Here is the call graph for this function:

◆ event_defaults()

void event_defaults ( void )

Event: defaults (i = 0)

Initialisation

We set the default values.

If the viscosity is non-zero, we need to allocate the face-centered viscosity field.

We also initialize the list of tracers to be advected with the VOF function \(f\) (or its complementary function).

We set limiting.

On trees, we ensure that limiting is also applied to prolongation and refinement.

We add the interface and the density to the default display.

Definition at line 150 of file two-phase.h.

References _i, alpha, alphav, display(), f, fE1, fE2, frho1, frho2, lambdav, list_copy(), minmod2(), mu, mu1, mu2, p, refine_linear(), rho, rhoc2, rhoc2v, rhov, s, and x.

Here is the call graph for this function:

◆ event_end_timestep()

void event_end_timestep ( void )

Event: end_timestep (i++)

Update of the total energy

After projection we update the values of the total energy adding the term missing from the projection step.

For a Newtonian fluid, the stress tensor comprises the pressure term, the viscous uncompressible term and the viscous compressible term,

\[ \tau = -p \mathbf{I} + \mu (\nabla \mathbf{u} + \nabla \mathbf{u}^T) + \lambda_v (\nabla \cdot \mathbf{u}) \mathbf{I} \]

where \(\mathbf{I}\) is the unity tensor and \(\lambda_v = \mu_v - 2/3 \mu\). The volumetric viscosity \(\mu_v\) is negligible in most cases.

The contribution of the pressure to the energy of each phase is lacking.

This is the contribution of the incompressible viscous term.

The velocity \(\mathbf{u}\) is first obtained from the updated momentum \(\mathbf{q}\).

We add the incompressible contribution of the \(\nabla \cdot (u_i \tau'_i)\) term. Note that the compressible contribution, which is typically small, is neglected.

Finally we recover momentum.

Definition at line 404 of file two-phase.h.

References _i, clamp(), cm, dimension, dt, energy(), f, fE1, fE2, frho1, frho2, lambdav, left, momentum(), mu1, mu2, p, q, right, u, uf, vector::x, x, and vector::y.

Here is the call graph for this function:

◆ event_init()

void event_init ( void )

Event: init (i = 0)

We update the fluid properties.

We set the initial timestep (this is useful only when restoring from a previous run).

For the associated tracers we use the gradient defined by f.gradient.

Definition at line 202 of file two-phase.h.

References _i, DT, dtmax, event(), f, fm, frho1, frho2, q, s, uf, vector::x, and x.

Here is the call graph for this function:

◆ event_properties()

void event_properties ( void )

Event: properties (i++)

Pressure and density

During the projection step we compute the provisional pressure from the EOS using the values after advection.

We also compute \(\rho c^2\).

If viscosity is present we obtain the averaged viscosity at the cell faces.

We also compute the averaged density at the cell faces.

The all-Mach solver needs rho*cm.

Definition at line 344 of file two-phase.h.

References _i, alphav, average_pressure(), bulk_compressibility(), cm, f, fm, frho1, frho2, lambdav, mu, mu1, mu2, muv, point, ps, rhoc2v, rhov, vector::x, and x.

Here is the call graph for this function:

◆ event_stability()

void event_stability ( void )

Event: stability (i++)

We need to overload the stability event so that the CFL is taken into account (because we set stokes to true).

Time step restriction based on the speed of sound

Definition at line 114 of file two-phase.h.

References c, CFLac, dtmax, HUGE, min, point, sound_speed(), and x.

Here is the call graph for this function:

◆ event_tracer_advection()

void event_tracer_advection ( void )

We set the list of interfaces back to its default value.

Event: tracer_advection (i++)

Definition at line 333 of file two-phase.h.

References interfaces, and interfaces1.

◆ event_vof()

void event_vof ( void )

Event: vof (i++)

We split the total momentum \(q\) into its two components \(q1\) and \(q2\) associated with \(f\) and \(1 - f\) respectively.

Momentum \(q2\) is associated with \(1 - f\), so we set the inverse attribute to true. We use the same limiting for q1 and q2.

The refinement function is modified by vof_advection(). To be able to restore it, we store its value.

We associate the transport of \(q1\) and \(q2\) with \(f\) and transport all fields consistently using the VOF scheme.

We recover the total momentum.

We avoid negative densities and energies which may have been caused by round-off during VOF advection.

We restore the refinement function for the total momentum.

We set the list of interfaces to NULL so that the default vof() event does nothing (otherwise we would transport \(f\) twice).

Definition at line 238 of file two-phase.h.

References _i, dimension, f, fE1, fE2, fE_refine(), free(), frho1, frho2, i, interfaces, interfaces1, list_concat(), q, q1, q2, refine(), s, set_prolongation(), tracers, u, vof_advection(), vector::x, and x.

Here is the call graph for this function:

◆ fE_refine()

void fE_refine	(	Point	point,
		scalar	fE
	)

Energy refinement function

The energy is refined from the refined pressures, momentum and densities, using the equation of state.

Definition at line 130 of file two-phase.h.

References dimension, f, frho1, frho2, internal_energy(), point, q, sq(), vector::x, and x.

Referenced by event_vof().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ internal_energy()

double internal_energy	(	Point	point,
		double	fc
	)

extern

Internal energy

Definition at line 85 of file Mie-Gruneisen.h.

References p, PIGAMMA, and x.

Referenced by fE_refine().

Here is the caller graph for this function:

◆ sound_speed()

double sound_speed ( Point point )

extern

These functions are provided by the Equation Of State.

See for example the Mie–Gruneisen Equation of State.

Sound speed

In mixture cells, this function returns the maximum between the speeds in both phases.

Sound speed

In mixture cells, this function returns the maximum between the speeds in both phases.

Definition at line 28 of file Mie-Gruneisen.h.

References clamp(), dimension, f, fE1, fE2, frho1, frho2, gamma1, gamma2, max, p, PI1, PI2, q, sq(), vector::x, and x.

Referenced by event_stability().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

◆ alphav

vector alphav[]

Definition at line 105 of file two-phase.h.

Referenced by event_defaults(), and event_properties().

◆ CFLac

double CFLac = HUGE

The timestep can be limited by a CFL based on the speed of sound.

Definition at line 62 of file two-phase.h.

Referenced by event_stability().

◆ f

scalar f[]

The primary fields are:

\[ \begin{aligned} \text{frho1} & = f \rho_1, \\ \text{frho2} & = (1-f) \rho_2, \\ \text{q} & = f \rho_1 \mathbf{u} + (1-f) \rho_2 \mathbf{u}, \\ \text{fE1} & = f \rho_1 (e_1 + \mathbf{u}^2/2), \\ \text{fE2} & = (1-f) \rho_2 (e_2 + \mathbf{u}^2/2) \end{aligned} \]

Definition at line 56 of file two-phase.h.

Referenced by advance(), advect(), advection(), advection_centered(), advection_upwind(), average_pressure(), average_temperature(), bulk_compressibility(), delete(), event_acceleration(), event_cleanup(), event_defaults(), event_end_timestep(), event_tracer_advection(), event_tracer_diffusion(), event_velocity(), event_vof(), face_fraction(), fclone(), flux(), fraction(), free_solver(), gradients(), height_position(), if(), if(), kdt_query_sum(), kpp_bottom_layer(), kpp_do_kpp(), kpp_surface_layer(), momentum(), momentum_refine(), momentum_restriction(), no_coalescence(), normf(), ohmic_flux(), okada(), output_fluxes(), periodic_function(), pmfunc_alloc(), pmfunc_free(), pmfunc_trace(), push(), query_sum(), remap_central(), remap_neumann(), remap_robin(), remove_droplets(), riemann(), right_value(), run(), sound_speed(), statsf(), thermal_expansion(), tracer_fluxes(), update_conservation(), update_tracer(), vof_concentration_refine(), which(), and zarea().

◆ fE1

scalar fE1[]

Definition at line 57 of file two-phase.h.

Referenced by average_pressure(), event_defaults(), event_end_timestep(), event_vof(), and sound_speed().

◆ fE2

scalar fE2[]

Definition at line 57 of file two-phase.h.

Referenced by average_pressure(), event_defaults(), event_end_timestep(), event_vof(), and sound_speed().

◆ frho1

scalar frho1[]

Definition at line 57 of file two-phase.h.

Referenced by average_pressure(), average_temperature(), bulk_compressibility(), event_acceleration(), event_defaults(), event_end_timestep(), event_init(), event_properties(), event_vof(), fE_refine(), and sound_speed().

◆ frho2

scalar frho2[]

Definition at line 57 of file two-phase.h.

Referenced by average_pressure(), average_temperature(), bulk_compressibility(), event_acceleration(), event_defaults(), event_end_timestep(), event_init(), event_properties(), event_vof(), fE_refine(), and sound_speed().

◆ interfaces

scalar * interfaces = {f}

The height functions are stored in the vector field associated with each VOF tracer.

We will need basic functions for volume fraction computations.

They need to be updated every time the VOF field changes. For the centered Navier-Stokes solver, this means after initialisation and after VOF advection.

Note that strictly speaking this should be done for each sweep of the direction-split VOF advection, which we do not do here i.e. we use the normal at the beginning of the timestep and assume it is constant during each sweep. This seems to work fine.

Definition at line 56 of file two-phase.h.

Referenced by event_cleanup(), event_defaults(), event_tracer_advection(), event_vof(), and no_coalescence().

◆ interfaces1

scalar* interfaces1 = NULL

static

VOF advection of momentum

We overload the vof() event to transport consistently the volume fraction and the momentum of each phase.

Definition at line 235 of file two-phase.h.

Referenced by event_tracer_advection(), and event_vof().

◆ lambdav

const vector lambdav = lambdav0

Definition at line 108 of file two-phase.h.

◆ lambdav0

const vector lambdav0[] = {0,0,0}

Definition at line 107 of file two-phase.h.

◆ lambdav1

double lambdav1 = 0.

Definition at line 69 of file two-phase.h.

◆ lambdav2

double lambdav2 = 0.

Definition at line 69 of file two-phase.h.

◆ mu1

double mu1 = 0.

The dynamic viscosities for each phase, as well as the volumetric viscosity coefficients.

Definition at line 68 of file two-phase.h.

Referenced by event_defaults(), event_end_timestep(), and event_properties().

◆ mu2

double mu2 = 0.

Definition at line 68 of file two-phase.h.

Referenced by event_defaults(), event_end_timestep(), and event_properties().

◆ rhoc2v

scalar rhoc2v[]

Auxilliary fields

Auxilliary fields need to be allocated. The quantity \(\rho c^2\), the average density rhov and its inverse alphav.

Definition at line 104 of file two-phase.h.

Referenced by event_defaults(), and event_properties().

◆ rhov

scalar rhov[]

Definition at line 106 of file two-phase.h.

Referenced by event_defaults(), and event_properties().

Macros

Functions

Variables

Macro Definition Documentation

◆ lambdav

◆ mu

◆ NO_1D_COMPRESSION

A compressible two-phase flow solver

Function Documentation

◆ average_pressure()

◆ bulk_compressibility()

Bulk compressibility of the mixture

Bulk compressibility of the mixture

◆ event_cleanup()

◆ event_defaults()

Initialisation

◆ event_end_timestep()

Update of the total energy

◆ event_init()

◆ event_properties()

Pressure and density

◆ event_stability()

Time step restriction based on the speed of sound

◆ event_tracer_advection()

◆ event_vof()

◆ fE_refine()

Energy refinement function

◆ internal_energy()

Internal energy

Internal energy

◆ sound_speed()

Sound speed

Sound speed

Variable Documentation

◆ alphav

◆ CFLac

◆ f

◆ fE1

◆ fE2

◆ frho1

◆ frho2

◆ interfaces

◆ interfaces1

VOF advection of momentum

◆ lambdav

◆ lambdav0

◆ lambdav1

◆ lambdav2

◆ mu1

◆ mu2

◆ rhoc2v

Auxilliary fields

◆ rhov