9.0
general documentation
Loading...
Searching...
No Matches
cs_vof.cpp File Reference
#include "base/cs_defs.h"
#include "bft/bft_mem.h"
#include "base/cs_array.h"
#include "base/cs_base.h"
#include "alge/cs_blas.h"
#include "base/cs_boundary_conditions.h"
#include "alge/cs_convection_diffusion.h"
#include "alge/cs_convection_diffusion_priv.h"
#include "cs_dispatch.h"
#include "alge/cs_divergence.h"
#include "cdo/cs_equation.h"
#include "base/cs_equation_iterative_solve.h"
#include "alge/cs_face_viscosity.h"
#include "base/cs_field.h"
#include "base/cs_field_default.h"
#include "base/cs_field_operator.h"
#include "base/cs_field_pointer.h"
#include "base/cs_log.h"
#include "base/cs_log_iteration.h"
#include "base/cs_math.h"
#include "mesh/cs_mesh.h"
#include "mesh/cs_mesh_location.h"
#include "mesh/cs_mesh_quantities.h"
#include "base/cs_parall.h"
#include "base/cs_physical_constants.h"
#include "base/cs_prototypes.h"
#include "base/cs_rotation.h"
#include "alge/cs_sles_default.h"
#include "base/cs_turbomachinery.h"
#include "turb/cs_turbulence_model.h"
#include "base/cs_vof.h"
Include dependency graph for cs_vof.cpp:

Functions

cs_vof_parameters_t * cs_get_glob_vof_parameters (void)
void cs_vof_field_create (void)
 Create VoF fields.
void cs_vof_log_setup (void)
 Log setup of VoF model.
void cs_vof_compute_linear_rho_mu (const cs_mesh_t *m)
 Compute the mixture density, mixture dynamic viscosity given fluid volume fractions and the reference density and dynamic viscosity $ \rho_l, \mu_l $ (liquid), $ \rho_v, \mu_v $ (gas).
void cs_vof_update_phys_prop (const cs_mesh_t *m)
 Compute the mixture density, mixture dynamic viscosity and mixture mass flux given the volumetric flux, the volume fraction and the reference density and dynamic viscosity $ \rho_l, \mu_l $ (liquid), $ \rho_v, \mu_v $ (gas).
void cs_vof_log_mass_budget (const cs_mesh_t *m, const cs_mesh_quantities_t *mq)
 Write in main log the global mixture mass budget:
void cs_vof_surface_tension (const cs_mesh_t *m, const cs_mesh_quantities_t *mq, cs_real_3_t stf[])
 Compute the surface tension momentum source term following the CSF model of Brackbill et al. (1992).
void cs_vof_deshpande_drift_flux (const cs_mesh_t *m, const cs_mesh_quantities_t *mq)
 Compute a relative velocity $ \vect u _d $ directly at internal faces (drift flux), following the approach described by Suraj S. Deshpande et al 2012 Comput. Sci. Disc. 5 014016. It is activated with the option idrift = 1.
void cs_vof_drift_term (int imrgra, int nswrgp, int imligp, int iwarnp, cs_real_t epsrgp, cs_real_t climgp, cs_real_t *restrict pvar, const cs_real_t *restrict pvara, cs_real_t *restrict rhs)
 Add the divergence of the drift velocity term in the volume fraction equation.
void cs_vof_solve_void_fraction (int iterns)
 Solve the void fraction $ \alpha $ for the Volume of Fluid method (and hence for cavitating flows).
cs_cavitation_parameters_t * cs_get_glob_cavitation_parameters (void)
void cs_cavitation_compute_source_term (const cs_real_t pressure[], const cs_real_t voidf[])
 Compute the vaporization source term $ \Gamma_V \left(\alpha, p\right) = m^+ + m^- $ using the Merkle model:
void cs_vof_contact_angle_set (const cs_vof_contact_angle_t choice)
 Set the contact angle mode.

Detailed Description

VOF model data.

Function Documentation

◆ cs_cavitation_compute_source_term()

void cs_cavitation_compute_source_term ( const cs_real_t pressure[],
const cs_real_t voidf[] )

Compute the vaporization source term $ \Gamma_V \left(\alpha, p\right) = m^+ + m^- $ using the Merkle model:

\[m^+ = -\dfrac{C_{prod} \rho_l \min \left( p-p_V,0 \right)\alpha(1-\alpha)} {0.5\rho_lu_\infty^2t_\infty}, \]

\[m^- = -\dfrac{C_{dest} \rho_v \max \left( p-p_V,0 \right)\alpha(1-\alpha)} {0.5\rho_lu_\infty^2t_\infty}, \]

with $ C_{prod}, C_{dest} $ empirical constants, $ t_\infty=l_\infty/u_\infty $ a reference time scale and $p_V$ the reference saturation pressure. $ l_\infty $, $ u_\infty $ and $p_V$ may be provided by the user (user function).

Note that the r.h.s. of the void fraction transport equation is $ \Gamma_V/\rho_v $.

Parameters
[in]pressurePressure array
[in]voidfVoid fraction array

◆ cs_get_glob_cavitation_parameters()

cs_cavitation_parameters_t * cs_get_glob_cavitation_parameters ( void )

◆ cs_get_glob_vof_parameters()

cs_vof_parameters_t * cs_get_glob_vof_parameters ( void )

◆ cs_vof_compute_linear_rho_mu()

void cs_vof_compute_linear_rho_mu ( const cs_mesh_t * m)

Compute the mixture density, mixture dynamic viscosity given fluid volume fractions and the reference density and dynamic viscosity $ \rho_l, \mu_l $ (liquid), $ \rho_v, \mu_v $ (gas).

Computation is done as follows on cells:

\[\rho_\celli = \alpha_\celli \rho_v + (1-\alpha_\celli) \rho_l, \]

\[\mu_\celli = \alpha_\celli \mu_v + (1-\alpha_\celli) \mu_l, \]

A similar linear formula is followed on boundary using fluid volume fraction value on the boundary.

Parameters
[in]mpointer to mesh structure

◆ cs_vof_contact_angle_set()

void cs_vof_contact_angle_set ( const cs_vof_contact_angle_t choice)

Set the contact angle mode.

Parameters
[in]choiceContact angle model to set.

◆ cs_vof_deshpande_drift_flux()

void cs_vof_deshpande_drift_flux ( const cs_mesh_t * m,
const cs_mesh_quantities_t * mq )

Compute a relative velocity $ \vect u _d $ directly at internal faces (drift flux), following the approach described by Suraj S. Deshpande et al 2012 Comput. Sci. Disc. 5 014016. It is activated with the option idrift = 1.

Compute the flux of the drift velocity $ \vect u _d $, by using the flux of the standard velocity $ \vect u $ following the approach described by Suraj S Deshpande et al 2012 Comput. Sci. Disc. 5 014016. It is activated with the option idrift = 1.

Using the notation:

\[\begin{cases} \left ( \vect u ^{n+1} . \vect S \right ) _{\face} = \Dot{m}_{\face}\\ \left ( \vect u _d^{n+1} . \vect S \right ) _{\face} = \Dot{m^d}_{\face} \end{cases} \]

The drift flux is computed as:

\[\Dot{m^d}_{\face} = min \left ( C_{\gamma} \dfrac{\Dot{m}_{\face}} {\vect S_{\face}}, \underset{\face'}{max} \left [ \dfrac{\Dot{m}_{\face'}} {\vect S_{\face'}} \right ] \right ) \left ( \vect n \cdot \vect S \right ) _{\face} \]

Where $ C_{\gamma} $ is the drift flux factor defined with the variable cdrift, $ \vect n _{\face} $ the normal vector to the interface. The gradient is computed using a centered scheme:

\[\vect {n} _{\face} = \dfrac{\left ( \grad \alpha \right ) _{\face}} {\norm {\left ( \grad \alpha \right ) _{\face} + \delta}}, \text{ with: } \left ( \grad \alpha \right ) _{\face _{\celli \cellj}} = \dfrac{\left ( \grad \alpha \right ) _\celli + \left ( \grad \alpha \right ) _\cellj}{2}, \text{ and: } \delta = 10^{-8} / \overline{\vol \celli} ^{1/3} \]

Parameters
[in]mpointer to mesh structure
[in]mqpointer to mesh quantities structure

◆ cs_vof_drift_term()

void cs_vof_drift_term ( int imrgra,
int nswrgp,
int imligp,
int iwarnp,
cs_real_t epsrgp,
cs_real_t climgp,
cs_real_t *restrict pvar,
const cs_real_t *restrict pvara,
cs_real_t *restrict rhs )

Add the divergence of the drift velocity term in the volume fraction equation.

More precisely, the right hand side $ Rhs $ is updated as follows:

\[Rhs = Rhs - \sum_{\fij \in \Facei{\celli}} \left( \alpha_\celli^{n+1} \left( 1 - \alpha_\cellj^{n+1} \right) \left( \dot{m}_\fij^{d} \right)^{+} + \alpha_\cellj^{n+1} \left( 1 - \alpha_\celli^{n+1} \right) \left( \dot{m}_\fij^{d} \right)^{-} \right) \]

Parameters
[in]imrgraindicator
  • 0 iterative gradient
  • 1 least squares gradient
[in]nswrgpnumber of reconstruction sweeps for the gradients
[in]imligpclipping gradient method
  • < 0 no clipping
  • = 0 by neighboring gradients
  • = 1 by the mean gradient
[in]iwarnpverbosity
[in]epsrgprelative precision for the gradient reconstruction
[in]climgpclipping coefficient for the computation of the gradient
[in]pvarsolved variable (current time step)
[in]pvarasolved variable (previous time step)
[in,out]rhsright hand side $ \vect{Rhs} $

◆ cs_vof_field_create()

void cs_vof_field_create ( void )

Create VoF fields.

◆ cs_vof_log_mass_budget()

void cs_vof_log_mass_budget ( const cs_mesh_t * m,
const cs_mesh_quantities_t * mq )

Write in main log the global mixture mass budget:

\[\sum_\celli\left( |\Omega_\celli|\dfrac{\rho_\celli^n - \rho_\celli^{n-1}}{\Delta t} + \sum_{\cellj\in\Face{\celli}}\left(\rho\vect{u}\vect{S}\right)_{\ij}^n \right). \]

Parameters
[in]mpointer to mesh structure
[in]mqpointer to mesh quantities structure

◆ cs_vof_log_setup()

void cs_vof_log_setup ( void )

Log setup of VoF model.

◆ cs_vof_solve_void_fraction()

void cs_vof_solve_void_fraction ( int iterns)

Solve the void fraction $ \alpha $ for the Volume of Fluid method (and hence for cavitating flows).

This function solves:

\[\dfrac{\alpha^n - \alpha^{n-1}}{\Delta t} + \divs \left( \alpha^n \vect{u}^n \right) + \divs \left( \left[ \alpha^n \left( 1 - \alpha^{n} \right) \right] \vect{u^r}^n \right) = \dfrac{\Gamma_V \left( \alpha^{n-1}, p^n \right)}{\rho_v} \]

with $ \Gamma_V $ the eventual vaporization source term (Merkle model) in case the cavitation model is enabled, $ \rho_v $ the reference gas density and $ \vect{u^r} $ the drift velocity for the compressed interface.

Parameters
[in]iternsNavier-Stokes iteration number

◆ cs_vof_surface_tension()

void cs_vof_surface_tension ( const cs_mesh_t * m,
const cs_mesh_quantities_t * mq,
cs_real_3_t stf[] )

Compute the surface tension momentum source term following the CSF model of Brackbill et al. (1992).

Parameters
[in]mpointer to mesh structure
[in]mqpointer to mesh quantities structure
[out]stfsurface tension momentum source term

◆ cs_vof_update_phys_prop()

void cs_vof_update_phys_prop ( const cs_mesh_t * m)

Compute the mixture density, mixture dynamic viscosity and mixture mass flux given the volumetric flux, the volume fraction and the reference density and dynamic viscosity $ \rho_l, \mu_l $ (liquid), $ \rho_v, \mu_v $ (gas).

For the computation of mixture density, mixture dynamic viscosity, see cs_vof_compute_linear_rho_mu.

Computation of mass flux is as follows:

\[\left( \rho\vect{u}\cdot\vect{S} \right)_\ij = \\ \left\lbrace \begin{array}{ll} \rho_\celli (\vect{u}\cdot\vect{S})_\ij &\text{ if } (\vect{u}\cdot\vect{S})_\ij>0, \\ \rho_\cellj (\vect{u}\cdot\vect{S})_\ij &\text{ otherwise }, \end{array} \right. \]

\[\left( \rho\vect{u}\cdot\vect{S} \right)_\ib = \\ \left\lbrace \begin{array}{ll} \rho_\celli (\vect{u}\cdot\vect{S})_\ib &\text{ if } (\vect{u}\cdot\vect{S})_\ib>0, \\ \rho_b (\vect{u}\cdot\vect{S})_\ib &\text{ otherwise }. \end{array} \right. \]

Parameters
[in]mpointer to mesh structure