#include "fractions.h"Go to the source code of this file.
Functions | |
| for (int _d=0;_d< 2;_d++) static double vof_concentration_gradient_x(Point point | |
The gradient of a VOF-concentration t is computed using a standard three-point scheme if we are far enough from the interface (as controlled by cmin), otherwise a two-point scheme biased away from the interface is used. | |
| if (t.inverse) cl | |
| if (cc >=cmin &&t.gradient !=zero) | |
| static void | vof_concentration_refine (Point point, scalar s) |
| On trees, VOF concentrations need to be refined properly i.e. | |
| void | event_defaults (void) |
| On trees, we need to setup the appropriate prolongation and refinement functions for the volume fraction fields. | |
| void | event_stability (void) |
| We need to make sure that the CFL is smaller than 0.5 to ensure stability of the VOF scheme. | |
| if (tracers) | |
| reconstruction (c, n, alpha) | |
| We reconstruct the interface normal \(\mathbf{n}\) and the intercept \(\alpha\) for each cell. | |
| delete (gfl) | |
| free (gfl) | |
| if (cfl > 0.5+1e-6) fprintf(ferr | |
| We warn the user if the CFL condition has been violated. | |
| src vof | fflush (ferr) |
| delete (tfluxl) | |
| free (tfluxl) | |
| void | vof_advection (scalar *interfaces, int i) |
| void | event_vof (void) vof_advection(interfaces |
| Event: vof (i++) | |
Variables | |
| attribute | |
| bool | inverse |
| scalar * | interfaces |
| We will need basic functions for volume fraction computations. | |
| vector | uf |
| We allocate the (face) velocity field. | |
| double | dt |
| These come from the multilayer solver. | |
| scalar | c |
| scalar scalar | t |
| double | cl = c[-1] |
| double | cc = c[] |
| double | cr = c[1] |
| return | |
| scalar scalar * | tcl |
| scalar | alpha [] |
| scalar | flux [] |
| double | cfl = 0. |
| scalar * | tracers = c.tracers |
| If we are also transporting tracers associated with \(c\), we need to compute their gradient i.e. | |
| scalar * | gfl = NULL |
| scalar * | tfluxl = NULL |
| src vof | h |
| src vof | LINENO |
| void | i |
Function Documentation
◆ delete() [1/2]
| delete | ( | gfl | ) |
◆ delete() [2/2]
| delete | ( | tfluxl | ) |
◆ event_defaults()
On trees, we need to setup the appropriate prolongation and refinement functions for the volume fraction fields.
Boundary conditions for VOF-advected tracers usually depend on boundary conditions for the VOF field.
Event: defaults (i = 0)
The default density field is set to unity (times the metric).
We reset the multigrid parameters to their 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.
We switch to a pure minmod limiter by default for increased robustness.
With the MUSCL scheme we use the CFL depends on the dimension of the problem.
On trees we need to replace the default bilinear refinement/prolongation with linear so that reconstructed values also use slope limiting.
The restriction/refine attributes of the charge density are those of a tracer otherwise the conservation is not guaranteed.
By default the permittivity is unity and other quantities are zero.
The (velocity) CFL is limited by the unsplit advection scheme, so is dependent on the dimension. The (gravity wave) CFL is set to 1/2 (if not already set by the user).
The gradient and prolongation/restriction functions are set for all tracer fields.
We setup the default display.
Boundary conditions
By default we set a zero Neumann boundary condition for all the components except if the bottom is an axis of symmetry.
We use (strict) minmod slope limiting for all components.
We reset the multigrid parameters to their default values.
The pressures are never dumped.
The default density field is set to unity (times the metric and the solid factors).
On trees, refinement of the face-centered velocity field needs to preserve the divergence-free condition.
When using embedded boundaries, the restriction and prolongation operators need to take the boundary into account.
We set the dimensions of the velocity field.
On trees, the refinement and restriction functions above rely on the volume fraction field f being refined/restricted before the components of velocity. To ensure this, we move f to the front of the field list (all).
We then set the refinement and restriction functions for the components of the velocity field. The boundary conditions on \(\mathbf{u}\) now depend on those on \(f\).
The default display.
We set limiting for q, h.
As well as default values.
On trees, we ensure that limiting is also applied to prolongation.
We setup the default display.
By default, all the layers have the same relative thickness.
We overload the default 'advance' and 'update' functions of the predictor-corrector scheme and setup the prolongation and restriction methods on trees.
On trees we make sure that slope-limiting is also used for refinement and prolongation. The prolongation/restriction functions for \(\eta\) are set and they depend on boundary conditions on \(z_b\) and \(h\).
We setup the default display.
If the viscosity is non-zero, we need to allocate the face-centered viscosity field.
We add the interface to the default display.
Definition at line 110 of file vof.h.
References c, fraction_refine(), restriction_volume_average(), set_prolongation(), set_restriction(), t, tracers, vof_concentration_refine(), and x.
◆ event_stability()
We need to make sure that the CFL is smaller than 0.5 to ensure stability of the VOF scheme.
We need to overload the stability event so that the CFL is taken into account (because we set stokes to true).
Event: stability (i++)
We first compute the minimum and maximum values of \(\alpha/f_m = 1/\rho\), as well as \(\Delta_{min}\).
The maximum timestep is set using the sum of surface tension coefficients.
We first compute the minimum and maximum values of \(\alpha/f_m = 1/\rho\), as well as \(\Delta_{min}\).
The maximum timestep is set using the sum of surface tension coefficients.
Definition at line 145 of file vof.h.
References CFL.
◆ event_vof()
◆ fflush()
Referenced by apply_bc(), balance(), cartesian_debug(), check_snd_rcv_matrix(), compdir(), display_command(), event__progress(), event_perfs(), gdb(), mg_solve(), mpi_recv_check(), msolve(), multigrid_debug(), output_fluxes(), output_gauges(), output_vtk(), rcv_pid_receive(), rcv_pid_send(), save(), solve(), stencil_val(), stencil_val_a(), and timer_print().
◆ for()
| for | ( | ) |
The gradient of a VOF-concentration t is computed using a standard three-point scheme if we are far enough from the interface (as controlled by cmin), otherwise a two-point scheme biased away from the interface is used.
Once we have computed the fluxes on all faces, we can update the volume fraction field according to the one-dimensional advection equation.
One-dimensional advection
The simplest way to implement a multi-dimensional VOF advection scheme is to use dimension-splitting i.e. advect the field along each dimension successively using a one-dimensional scheme.
We implement the one-dimensional scheme along the x-dimension and use the for (int _d = 0; _d < dimension; _d++) operator to automatically derive the corresponding functions along the other dimensions.
\[ \partial_tc = -\nabla_x\cdot(\mathbf{u}_f c) + c\nabla_x\cdot\mathbf{u}_f \]
The first term is computed using the fluxes. The second term – which is non-zero for the one-dimensional velocity field – is approximated using a centered volume fraction field cc which will be defined below.
For tracers, the corresponding one-dimensional update is
\[ \partial_tt_j = -\nabla_x\cdot(\mathbf{u}_f t_j) + t_j\nabla_x\cdot\mathbf{u}_f \]
The compressive second-term can be removed by defining the NO_1D_COMPRESSION macro.
To compute the volume fraction flux, we check the sign of the velocity component normal to the face and compute the index i of the corresponding upwind cell (either 0 or -1).
We also check that we are not violating the CFL condition.
If we assume that un is negative i.e. s is -1 and i is 0, the volume fraction flux through the face of the cell is given by the dark area in the figure below. The corresponding volume fraction can be computed using the rectangle_fraction() function.
When the upwind cell is entirely full or empty we can avoid this computation.
Once we have the upwind volume fraction cf, the volume fraction flux through the face is simply:
If we are transporting tracers, we compute their flux using the upwind volume fraction cf and a tracer value upwinded using the Bell–Collela–Glaz scheme and the gradient computed above.
Definition at line 200 of file vof.h.
References alpha, c, cfl, cm, cs, dt, flux, fm, gfl, i, min, n, rectangle_fraction(), s, SEPS, sign(), t, tfluxl, tracers, uf, un, vector::x, and x.
◆ free() [1/2]
| free | ( | gfl | ) |
◆ free() [2/2]
| free | ( | tfluxl | ) |
◆ if() [1/4]
| if | ( | cc >=cmin &&t.gradient ! | = zero | ) |
◆ if() [2/4]
◆ if() [3/4]
| if | ( | t. | inverse | ) |
◆ if() [4/4]
| if | ( | tracers | ) |
◆ reconstruction()
We reconstruct the interface normal \(\mathbf{n}\) and the intercept \(\alpha\) for each cell.
Then we go through each (vertical) face of the grid.
◆ vof_advection()
Multi-dimensional advection
The multi-dimensional advection is performed by the event below.
We first define the volume fraction field used to compute the divergent term in the one-dimensional advection equation above. We follow Weymouth & Yue, 2010 and use a step function which guarantees exact mass conservation for the multi-dimensional advection scheme (provided the advection velocity field is exactly non-divergent).
We then apply the one-dimensional advection scheme along each dimension. To try to minimise phase errors, we alternate dimensions according to the parity of the iteration index i.
Definition at line 330 of file vof.h.
References _i, c, cc, d, dimension, free(), i, list_append(), restriction_volume_average(), set_prolongation(), set_restriction(), t, tcl, tracers, vof_concentration_refine(), and x.
Referenced by event_vof().
◆ vof_concentration_refine()
On trees, VOF concentrations need to be refined properly i.e.
using volume-fraction-weighted linear interpolation of the concentration.
Definition at line 85 of file vof.h.
References cm, dimension, f, g, point, s, vector::x, and x.
Referenced by event_defaults(), and vof_advection().
Variable Documentation
◆ alpha
◆ attribute
| attribute |
Volume-Of-Fluid advection
We want to approximate the solution of the advection equations
\[ \partial_tc_i + \mathbf{u}_f\cdot\nabla c_i = 0 \]
where \(c_i\) are volume fraction fields describing sharp interfaces.
This can be done using a conservative, non-diffusive geometric VOF scheme.
We also add the option to transport diffusive tracers (aka "VOF concentrations") confined to one side of the interface i.e. solve the equations
\[ \partial_tt_{i,j} + \nabla\cdot(\mathbf{u}_ft_{i,j}) = 0 \]
with \(t_{i,j} = c_if_j\) (or \(t_{i,j} = (1 - c_i)f_j\)) and \(f_j\) is a volumetric tracer concentration (see Lopez-Herrera et al, 2015).
The list of tracers associated with the volume fraction is stored in the tracers attribute. For each tracer, the "side" of the interface (i.e. either \(c\) or \(1 - c\)) is controlled by the inverse attribute).
◆ c
| scalar c |
Definition at line 57 of file vof.h.
Referenced by add_cexpr(), args(), bubbles_are_close(), cache_append(), cache_level_append(), cache_level_init(), cache_level_shrink(), cache_shrink(), centroids_curvature_fit(), cfilter(), coeffs1(), coeffs2(), colormap_color(), comment(), comment(), curvature(), default_stencil(), display_send(), echo_c(), eigenvalues(), event_defaults(), event_face_fields(), event_init(), event_stability(), event_tracer_advection(), event_tracer_diffusion(), event_vof(), facet_normal(), fclone(), flux(), for(), foreach_tree(), fraction_refine(), fractions_cleanup(), free_cexpr(), get_cexpr(), getput(), getput(), givens(), h_relax(), half_column(), height_curvature(), height_curvature_fit(), height_myc_normal(), height_normal(), heights(), if(), if(), inclex(), input(), input(), input0(), input_gfs(), input_xy(), interface_area(), interface_normal(), interface_normal(), interfacial(), line_alpha(), list_clone(), main(), matrix(), matrix1(), merge(), mycs(), next_char(), next_string(), no_coalescence(), ohmic_flux(), periodic(), pmfunc_alloc(), pmfunc_free(), pmfunc_trace(), PointTriangleDistance(), postlex(), quadratic(), receive_tree(), reconstruction(), relax(), relax_diffusion(), relax_nh(), relax_nh1(), replace_(), residual(), residual_diffusion(), residual_nh(), residual_nh3(), right_value(), solve_hessenberg(), tracer_is_close(), update_distance(), update_green_naghdi(), util_matrix(), vertical_diffusion(), vertical_viscosity(), vof_advection(), volumez(), youngs_normal(), yyunput(), yyunput(), and yyunput().
◆ cc
Definition at line 60 of file vof.h.
Referenced by embed_fraction_refine(), event_pressure(), fraction_refine(), if(), main(), util_findp_bisection(), and vof_advection().
◆ cfl
◆ cl
◆ cr
◆ dt
|
extern |
These come from the multilayer solver.
Ohmic conduction
This function approximates implicitly the EHD equation set given by the electric Poisson equation and the ohmic conduction term of the charge density conservation (the advection term is computed elsewhere using a tracer advection scheme),
\[ \nabla \cdot( \epsilon \nabla \phi^{n+1}) = -\rho_e^{n+1} \quad \mathrm{and} \quad (\rho_e^{n+1}-\rho_e^n) = \Delta t\nabla \cdot (K\nabla \phi^{n+1}) \]
where \(\rho_e\) is the charge density, and \(\phi\) is the electric potential. \(K\) and \(\epsilon\) are the conductivity and permittivity respectively.
Substituting the Poisson equation into the conservation equation gives the following time-implicit scheme,
\[ \nabla \cdot [(K \, \Delta t + \epsilon) \nabla \phi^{n+1}] = -\rho_e^{n} \]
We need the Poisson solver and the timestep dt.
Generic time loop
The run() function below implements a generic time loop which executes events until termination.
The timestep dt can be accessed as a global variable.
Definition at line 18 of file predictor-corrector.h.
Referenced by advance_generic(), for(), run(), and run().
◆ flux
| scalar flux[] |
Definition at line 166 of file vof.h.
Referenced by advect(), advection(), boundary_internal(), boundary_stencil(), event_viscous_term(), for(), if(), output_fluxes(), riemann(), segment_flux(), tracer_fluxes(), update_tracer(), and vertical_fluxes().
◆ gfl
◆ h
◆ i
| void i |
Definition at line 387 of file vof.h.
Referenced by for(), and vof_advection().
◆ interfaces
|
extern |
We will need basic functions for volume fraction computations.
The list of volume fraction fields interfaces, will be provided by the user.
The face velocity field uf will be defined by a solver as well as the timestep.
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().
◆ inverse
◆ LINENO
◆ return
◆ t
Definition at line 57 of file vof.h.
Referenced by event_defaults(), for(), if(), if(), and vof_advection().
◆ tcl
Definition at line 163 of file vof.h.
Referenced by vof_advection().
◆ tfluxl
◆ tracers
| scalar* tracers = c.tracers |
If we are also transporting tracers associated with \(c\), we need to compute their gradient i.e.
Here we set the gradient functions for each tracer (as defined in the user-provided tracers list).
\(\partial_xf_j = \partial_x(t_j/c)\) or \(\partial_xf_j = \partial_x(t_j/(1 - c))\) (for higher-order upwinding) and we need to store the computed fluxes. We first allocate the corresponding lists.
Definition at line 176 of file vof.h.
Referenced by advection(), event_defaults(), event_half_advection(), event_remap(), event_tracer_advection(), event_vof(), for(), if(), vertical_remapping(), and vof_advection().
◆ uf
|
extern |
We allocate the (face) velocity field.
An advection solver
We wish to solve the advection equations
\[ \partial_tf_i+\mathbf{u}\cdot\nabla f_i=0 \]
where \(\mathbf{u}\) is the velocity field and \(f_i\) are a list of passive tracers. This can be done with a flux-based advection scheme such as the 2nd-order, unsplit, upwind scheme of Bell-Collela-Glaz, 1989.
The main time loop is defined in [run.h](). A stable timestep needs to respect the CFL condition. For compatibility with the other solvers, we allocate it as uf and define an alias. The gradient function is used to set the type of slope-limiting required. The default is to not use any limiting (i.e. a purely centered slope estimation).
Definition at line 29 of file advection.h.
Referenced by event_acceleration(), event_adapt(), event_advection_term(), event_defaults(), event_end_timestep(), event_face_fields(), event_init(), event_pressure(), event_projection(), event_stability(), event_tracer_advection(), for(), prediction(), project(), tracer_fluxes(), and update_tracer().
