8.3
general documentation
cs_vof.cpp File Reference
#include "cs_defs.h"
#include "bft_mem.h"
#include "cs_array.h"
#include "cs_base.h"
#include "cs_blas.h"
#include "cs_boundary_conditions.h"
#include "cs_convection_diffusion.h"
#include "cs_convection_diffusion_priv.h"
#include "cs_divergence.h"
#include "cs_equation.h"
#include "cs_equation_iterative_solve.h"
#include "cs_face_viscosity.h"
#include "cs_field.h"
#include "cs_field_default.h"
#include "cs_field_operator.h"
#include "cs_field_pointer.h"
#include "cs_log.h"
#include "cs_log_iteration.h"
#include "cs_math.h"
#include "cs_mesh.h"
#include "cs_mesh_location.h"
#include "cs_mesh_quantities.h"
#include "cs_parall.h"
#include "cs_physical_constants.h"
#include "cs_prototypes.h"
#include "cs_rotation.h"
#include "cs_sles_default.h"
#include "cs_turbomachinery.h"
#include "cs_turbulence_model.h"
#include "cs_vof.h"
+ Include dependency graph for cs_vof.cpp:

Functions

cs_vof_parameters_tcs_get_glob_vof_parameters (void)
 
void cs_vof_field_create (void)
 Create VoF fields. More...
 
void cs_vof_log_setup (void)
 Log setup of VoF model. More...
 
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). More...
 
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). More...
 
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: More...
 
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). More...
 
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. More...
 
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...
 
void cs_vof_solve_void_fraction (int iterns)
 Solve the void fraction $ \alpha $ for the Volume of Fluid method (and hence for cavitating flows). More...
 
cs_cavitation_parameters_tcs_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: More...
 

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_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