Diagnostic quantities

Post-treatment routines to compute physical quantities (Nusselt number, Q criterion, kinetic energy, etc.). More...

Namespaces
module	mod_compute_ib_surface_diagnostics
	Compute immersed boundary surface diagnostics (surface, pressure drag, shear drag, nusselt_number, etc.).
module	enum_diagnostic_quantities
	Diagnostic quantity enumerators.
module	fields_diagnostic_quantities
	Diagnostic quantity fields.
module	variables_diagnostic_quantities
	Diagnostic quantity switches.

Functions
subroutine	mod_compute_lagrangian_acceleration::compute_lagrangian_acceleration (acceleration)
	Compute the lagrangian acceleration field on cells.
subroutine	mod_compute_dns_scales::compute_batchelor_scale (batchelor_scale, tke_dissipation_rate, viscosity, density, conductivity, specific_heat_cp, is_first_stat_iteration)
	Compute batchelor scale.
subroutine	mod_compute_dns_scales::compute_kolmogorov_scale (kolmogorov_scale, tke_dissipation_rate, viscosity, density, is_first_stat_iteration)
	Compute kolmogorov scale.
subroutine	mod_compute_dns_scales::compute_kolmogorov_timescale (kolmogorov_timescale, tke_dissipation_rate, viscosity, density, is_first_stat_iteration)
	Compute Kolmogorov timescale.
double precision function	mod_compute_domain_volume::compute_domain_volume (grid_volume)
	Compute the domain volume, as used by other physical tools.
subroutine	mod_compute_energy_budget::compute_energy_budget (internal_energy_domain, energy_from_bc, energy_loss)
	Compute energy budget.
double precision function	mod_compute_stefan_problem_front_position_melting::compute_stefan_problem_front_position_melting ()
	Compute front position.
double precision function	mod_compute_stefan_problem_front_position_solidification::compute_stefan_problem_front_position_solidification ()
	Compute front position.
double precision function, public	mod_compute_ib_surface_diagnostics::compute_ib_surface (cut_cell_id)
	Compute the surface of the immersed boundary.
double precision function, dimension(3), public	mod_compute_ib_surface_diagnostics::compute_ib_pressure_drag_lift (cut_cell_id)
	Compute the drag force due to pressure variations at the surface.
double precision function, dimension(3), public	mod_compute_ib_surface_diagnostics::compute_ib_shear_drag_lift (cut_cell_id)
	Compute the drag force due to shear at the surface.
subroutine	mod_compute_instantaneous_tke_dissipation_rate::compute_instantaneous_tke_dissipation_rate_incompressible (instant_tke_dissipation_rate, viscosity, velocity_fluctuation)
	Compute turbulent kinetic energy dissipation rate for incompressible flow.
subroutine	mod_compute_instantaneous_tke_dissipation_rate::compute_instantaneous_tke_dissipation_rate_compressible (instant_tke_dissipation_rate, viscosity, velocity_instantaneous, velocity_fluctuation)
	Compute turbulent kinetic energy dissipation rate for compressible flow.
subroutine	mod_compute_k_plus::compute_k_plus ()
	Compute k+.
subroutine	mod_compute_kinetic_energy::compute_kinetic_energy (kinetic_energy, velocity)
	Compute the kinetic energy field on cells.
subroutine	mod_compute_kinetic_energy_per_mass_unit::compute_kinetic_energy_per_mass_unit (kinetic_energy, velocity)
	Compute the kinetic energy field on cells.
subroutine	compute_local_wall_nusselt_number::compute_local_wall_nusselt (temperature, energy_bc)
	Compute local wall skin friction coef (physical and immersed boundaries)
subroutine	mod_compute_local_wall_skin_friction_coef::compute_local_wall_skin_friction_coef ()
	Compute local wall skin friction coef (physical and immersed boundaries)
double precision function	mod_compute_max_field::compute_max_field (field)
	Compute the max value of a field restricted to the inner domain.
double precision function, dimension(4)	mod_compute_max_velocity::compute_max_velocity (velocity)
	Compute the maximum velocity on the whole domain, for each component and for the velocity magnitude.
double precision function	mod_compute_mean_kinetic_energy::compute_mean_kinetic_energy (velocity)
	Compute the mean kinetic energy on the whole domain.
double precision function	mod_compute_mean_kinetic_energy::compute_kinetic_energy_integral (velocity)
	Compute the mean kinetic energy on the whole domain.
double precision function	mod_compute_mean_longitudinal_pressure_difference::compute_mean_longitudinal_pressure_difference (pressure)
	Compute the mean pressure difference between left and right boundaries.
double precision function	mod_compute_min_field::compute_min_field (field)
	Compute the min value of a field restricted to the inner domain.
subroutine	mod_compute_min_max_mean_velocity_magnitude::compute_min_max_mean_velocity_magnitude (velocity, min_velocity, max_velocity, mean_velocity)
	Compute the min, the max, and the min velocity magnitude on the whole domain.
double precision function, dimension(spatial_dimension)	mod_compute_momentum_integral::compute_momentum_integral (density_face, velocity)
	Compute the momentum \( \rho u \) integral.
subroutine	mod_compute_physical_diagnostics::compute_initial_physical_diagnostics ()
	Compute physical diagnostics before time loop.
subroutine	mod_compute_physical_diagnostics::compute_physical_diagnostics ()
	Execute enabled physical tools that returns a scalar value (Nusselt, mean quantities, etc.).
subroutine	mod_compute_q_criterion::compute_q_criterion (q_criterion, velocity)
	Compute Q-criterion of velocity field.
double precision function	mod_compute_spatial_averaged_field::compute_spatial_averaged_field (field, grid_volume)
	Compute the spatial averaged value of a field restricted to the inner domain.
subroutine	mod_compute_strain_rate_magnitude::compute_strain_rate_magnitude (strain_rate_magnitude, velocity)
	Compute strain rate magnitude of the velocity field.
subroutine, public	mod_compute_symmetric_difference_area::compute_symmetric_difference_area (phases1, phases2, area)
	Compute the area of symmetric difference between two phases.
subroutine, public	mod_compute_symmetric_difference_area::compute_field_symmetric_difference_area (phases1, phases2, area)
	Compute the area of symmetric difference between two phases. Returns a field.
subroutine	mod_compute_turbulent_kinetic_energy::compute_turbulent_kinetic_energy (turbulent_kinetic_energy, velocity_rms)
	Compute the turbulent_kinetic energy field on cells.
subroutine	mod_compute_u_plus::compute_u_plus ()
	Compute u+ (physical and immersed boundaries)
subroutine	mod_compute_viscous_dissipation_rate::compute_viscous_dissipation_rate (viscous_dissipation_rate, velocity, viscosity)
	Compute viscous dissipation rate per unit time and volume.
subroutine	mod_compute_vorticity_magnitude::compute_vorticity_magnitude (vorticity_magnitude, velocity)
	Compute the vorticity magnitude:
subroutine	mod_compute_wall_nusselt_sherwood_number::compute_wall_nusselt_sherwood_number (temperature, energy_bc, diagnostic_delta, left_nusselt, right_nusselt, bottom_nusselt, top_nusselt, back_nusselt, front_nusselt, ib_wall_nusselt_left, ib_wall_nusselt_right, ib_wall_nusselt_bottom, ib_wall_nusselt_top, ib_wall_nusselt_back, ib_wall_nusselt_front, ib_wall_nusselt)
	Compute wall Nusselt or Sherwood numbers.
subroutine	mod_compute_wall_shear_stress::compute_wall_shear_stress (velocity, wall_shear_stress_bottom, wall_shear_stress_top)
	Computes the wall shear stres on top and bottom boundaries (2D/3D)
subroutine	mod_compute_wall_shear_velocity::compute_wall_shear_velocity (wall_shear_stress_bottom, wall_shear_stress_top, wall_shear_velocity_bottom, wall_shear_velocity_top)
	Computes the wall skin friction coefficient.
subroutine	mod_compute_wall_skin_friction_coef::compute_wall_skin_friction_coef (velocity, wall_shear_stress_bottom, wall_shear_stress_top, wall_skin_friction_coef_bottom, wall_skin_friction_coef_top)
	Computes the wall skin friction coefficient.
subroutine	mod_compute_y_plus::compute_y_plus ()
	Compute y+ (physical and immersed boundaries)
subroutine	mod_drag_output::prepare_drag_output_file ()
	Print physical diagnostics.

Detailed Description

This directory contains routines that compute physical quantities:

• the kinetic energy defined on cells (compute_kinetic_energy.f90),
• the strain rate magnitude defined on cells (compute_strain_rate_magnitude.f90),
• the Q-criterion defined on cells (compute_q_criterion.f90).
• the lagrangian acceleration defined on cells (compute_acceleration.f90).

Physical diagnostics are scalar physical quantities. Routines that compute physical diagnostics are:

• the mean kinetic energy (compute_mean_kinetic_energy.f90),
• the mean longitudinal pressure difference (compute_mean_longitudinal_pressure_difference.f90),
• the mean kinetic temperature (compute_mean_temperature.f90),
• the mean velocity magnitude (compute_mean_velocity_magnitude.f90),
• the wall Nusselt and Sherwood numbers (compute_wall_nusselt_sherwood_number.f90).

More physical diagnostics can be computed via lower-level routines:

• the mean of a field defined on cells (compute_spatial_averaged_field.f90).

The list of physical diagnostics is in the module variables_diagnostic_quantities. This module contains switches that enable or disable their computation. Note the user interface is a convenient way the activate the switches. The routine 'compute_physical_diagnostics.f90`, called in the time loop, actually computes enabled physical diagnostics.

Function Documentation

compute_batchelor_scale()

subroutine mod_compute_dns_scales::compute_batchelor_scale	(	double precision, dimension(:,:,:), intent(inout)	batchelor_scale,
		double precision, dimension(:,:,:), intent(in)	tke_dissipation_rate,
		double precision, dimension(:,:,:), intent(in)	viscosity,
		double precision, dimension(:,:,:), intent(in)	density,
		double precision, dimension(:,:,:), intent(in)	conductivity,
		double precision, dimension(:,:,:), intent(in)	specific_heat_cp,
		logical, intent(in)	is_first_stat_iteration )

Compute batchelor scale is computed as \( \eta_B \).

\[ \eta_B = \eta_K / (\mathrm{Pr})^{1/2} = \Big( \frac{\nu \alpha^2}{\varepsilon_u}\Big)^{1/4} \]

with \( \overline{\nu} = \overline{mu} / \overline{\rho} \) the kinematic viscosity, \( \overline{\alpha}= \overline{\lambda} / (\overline{\rho} \overline{c_p}) \) the thermal diffusivity with \( \lambda \) \( \rho \) and \( c_p \) and \( \varepsilon_u \) the turbulent kinetic energy dissipation rate. The relation trought the Prandlt (or Schmidt for diffusion of species) is also given. From with scalar field, one can compute the time averaged Kolmogorov timescale \( \overline{\eta_B} \)

compute_initial_physical_diagnostics()

subroutine mod_compute_physical_diagnostics::compute_initial_physical_diagnostics

This routine is called before the time loop.

compute_instantaneous_tke_dissipation_rate_compressible()

subroutine mod_compute_instantaneous_tke_dissipation_rate::compute_instantaneous_tke_dissipation_rate_compressible	(	double precision, dimension(:,:,:), intent(inout)	instant_tke_dissipation_rate,
		double precision, dimension(:,:,:), intent(in)	viscosity,
		type(t_face_field), intent(in)	velocity_instantaneous,
		type(t_face_field), intent(in)	velocity_fluctuation )

Compute turbulent kinetic energy dissipation rate \( \varepsilon_u \) for compressible flow [2, eq (5.42)]:

\[ \overline{\rho} \varepsilon_u = \overline{\tau_{ij} S_{ij}^{''} } \]

with \( \tau_{ij} = 2 \mu S_{ij} - \frac{2\mu}{3} S_{ll} \delta_{ij} \) the viscous stress tensor of instantaneous velocity and \( S_{ij}^{''} \) the symmetric strain rate tensor based on velocity Favre fluctuations. Note that its subroutine only compute the instantaneous scalar field \( \tau_{ij} S_{ij}^{''} \) while the averaging is done in the "compute_statistics" subroutine.

Reference : [2] Turbulence Modeling for CFD, D.C. Wilcox, 2006, Chapter 5 Favre Averaging.

Terms from the diagonal

Term \(\text{I}\)

                   ┏━━━━━┓
    y              ┃     ┃
    │              >  o  >
    o──x           ┃     ┃
   z               ┗━━━━━┛
                   i  i  ip

\[ \Big( \frac {\partial u} {\partial x} \Big)_{\mathtt{i,j}} = \mathtt{ \frac {u_{ip,j} - u_{i,j} } {\Delta x_{i}} } \]

Terms \(\text{II}\) and \(\text{III}\) are computed likewise.

Terms from extra diag

Term \(\text{IV}\)

                ┌─────┐
                │     │
             jp >     >
                │     │
                ┢━━━━━┪     jp ┌──∧──┲━━∧━━┱──∧──┐
    y           ┃     ┃        │     ┃     ┃     │
    │         j >  o  >      j │     ┃  o  ┃     │
    o──x        ┃     ┃        │     ┃     ┃     │
   z            ┡━━━━━┩      j └──∧──┺━━∧━━┹──∧──┘
                │     │           im    i     ip
             jm >     >
                │     │
                └─────┘
                i  i  ip

\[ \Big( \frac {\partial u} {\partial y} \Big)_{\mathtt{i,j}} = \frac {1} {4} \Big( \mathtt{ \frac {u_{ip,jp} - u_{ip,j} } {\Delta yv_{jp}} } + \mathtt{ \frac {u_{i,jp} - u_{i,j} } {\Delta yv_{jp}} } + \mathtt{ \frac {u_{ip,j} - u_{ip,jm}} {\Delta yv_{j}} } + \mathtt{ \frac {u_{i,j} - u_{i,jm} } {\Delta yv_{j}} } \Big) \\ \Big( \frac {\partial v} {\partial x} \Big)_{\mathtt{i,j}} = \frac {1} {4} \Big( \mathtt{ \frac {v_{ip,jp} - v_{i,jp}} {\Delta xu_{ip}} } + \mathtt{ \frac {v_{ip,j} - v_{i,j} } {\Delta xu_{ip}} } + \mathtt{ \frac {v_{i,jp} - v_{im,jm}} {\Delta xu_{i}} } + \mathtt{ \frac {v_{i,j} - v_{im,j} } {\Delta xu_{i}} } \Big) \]

compute_instantaneous_tke_dissipation_rate_incompressible()

subroutine mod_compute_instantaneous_tke_dissipation_rate::compute_instantaneous_tke_dissipation_rate_incompressible	(	double precision, dimension(:,:,:), intent(inout)	instant_tke_dissipation_rate,
		double precision, dimension(:,:,:), intent(in)	viscosity,
		type(t_face_field), intent(in)	velocity_fluctuation )

Compute turbulent kinetic energy dissipation rate \( \varepsilon_u \) for incompressible flow [1, eq (4.6) so called 'true' dissipation rate]:

\[ \varepsilon_u = \overline{2\nu S_{ij}^{'} S_{ij}^{'}} \]

with \( \nu = \frac{\mu}{\rho} \), \( S_{ij}^{'} \) the symmetric strain rate tensor based on velocity fluctuations (Reynolds averaged). Note that its subroutine only compute the instantaneous scalar field \( 2\nu S_{ij}^{'} S_{ij}^{'} \) while the averaging is done in the "compute_statistics" subroutine.

Reference: [1] Turbulence Modeling for CFD, D.C. Wilcox, 2006, Chapter 4 One-Equation and Two-Equation Models.

compute_kinetic_energy()

subroutine mod_compute_kinetic_energy::compute_kinetic_energy	(	double precision, dimension(:,:,:), intent(out)	kinetic_energy,
		type(t_face_field), intent(in)	velocity )

The kinetic energy field is defined as:

\[ \frac {1} {2} \rho (u^2 + v^2 + w^2) \]

where the velocity components are interpolated on cell centers.

compute_kinetic_energy_integral()

double precision function mod_compute_mean_kinetic_energy::compute_kinetic_energy_integral

(

type(t_face_field), intent(in)

velocity

)

The mean velocity magnitude is defined as:

\[ \int_{V}^{} 1/2 \rho (u^2 + v^2 + w^2) \, \mathrm{d}v \]

where the velocity components are interpolated on cell centers.

compute_kinetic_energy_per_mass_unit()

subroutine mod_compute_kinetic_energy_per_mass_unit::compute_kinetic_energy_per_mass_unit	(	double precision, dimension(:,:,:), intent(out)	kinetic_energy,
		type(t_face_field), intent(in)	velocity )

The kinetic energy per mass unit field is defined as:

\[ \frac {1} {2} (u^2 + v^2 + w^2) \]

where the velocity components are interpolated on cell centers.

compute_kolmogorov_scale()

subroutine mod_compute_dns_scales::compute_kolmogorov_scale	(	double precision, dimension(:,:,:), intent(inout)	kolmogorov_scale,
		double precision, dimension(:,:,:), intent(in)	tke_dissipation_rate,
		double precision, dimension(:,:,:), intent(in)	viscosity,
		double precision, dimension(:,:,:), intent(in)	density,
		logical, intent(in)	is_first_stat_iteration )

Compute kolmogorov scale is computed as \( \eta_K \).

\[ \eta_K = \Big( \frac{\nu^3}{\varepsilon_u}\Big)^{1/4} \]

with \( \nu \) the kinematic viscosity and \( \varepsilon_u \) the turbulent kinetic energy dissipation rate. From with scalar field, one can compute the time averaged Kolmogorov timescale \( \overline{\eta_K} \)

compute_kolmogorov_timescale()

subroutine mod_compute_dns_scales::compute_kolmogorov_timescale	(	double precision, dimension(:,:,:), intent(inout)	kolmogorov_timescale,
		double precision, dimension(:,:,:), intent(in)	tke_dissipation_rate,
		double precision, dimension(:,:,:), intent(in)	viscosity,
		double precision, dimension(:,:,:), intent(in)	density,
		logical, intent(in)	is_first_stat_iteration )

Compute Kolmogorov timescale \( \tau_K \).

\[ \tau_K = \Big( \frac{\nu}{\varepsilon_u}\Big)^{1/2} \]

with \( \nu \) the kinematic viscosity and \( \varepsilon_u \) the turbulent kinetic energy dissipation rate. From with scalar field, one can compute the time averaged Kolmogorov timescale \( \overline{\tau_K} \)

compute_lagrangian_acceleration()

subroutine mod_compute_lagrangian_acceleration::compute_lagrangian_acceleration

(

double precision, dimension(:,:,:,:), intent(out)

acceleration

)

The lagrangian acceleration field is defined as:

\[ a_x = \frac {\partial u_x} {\partial t} + u_x \frac {\partial u} {\partial x} + v_y \frac {\partial u_x} {\partial y} + w_z \frac {\partial u_x} {\partial z} \]

\[ a_y = \frac {\partial v_y} {\partial t} + u_x \frac {\partial v_y} {\partial x} + v_y \frac {\partial v} {\partial y} + w_z \frac {\partial v_y} {\partial z} \]

\[ a_z = \frac {\partial w_z} {\partial t} + u_x \frac {\partial w_z} {\partial x} + v_y \frac {\partial w_z} {\partial y} + w_z \frac {\partial w} {\partial z} \]

where the velocity components ( \(u, v, w\)) are interpolated on cell centers ( \(u_x, v_y, w_z\)).

Parameters

[out]

acceleration

Lagrangian acceleration on cell centers.

compute_max_field()

double precision function mod_compute_max_field::compute_max_field

(

double precision, dimension(:,:,:), intent(in)

field

)

The argument field is defined on cells or faces.

compute_max_velocity()

double precision function, dimension(4) mod_compute_max_velocity::compute_max_velocity

(

type(t_face_field), intent(in)

velocity

)

The result is a vector:

• \((u_max, v_max, 0, |U| max)\) in 2D
• \((u_max, v_max, w_max, |U| max)\) in 3D

where the velocity components are interpolated on cell centers for computing \(|U| max\).

Parameters

[in]

velocity

velocity field

compute_mean_kinetic_energy()

double precision function mod_compute_mean_kinetic_energy::compute_mean_kinetic_energy

(

type(t_face_field), intent(in)

velocity

)

The mean velocity magnitude is defined as:

\[ \frac {1} {2V} \int_{V}^{} (u^2 + v^2 + w^2) \, \mathrm{d}v \]

where the velocity components are interpolated on cell centers.

compute_mean_longitudinal_pressure_difference()

double precision function mod_compute_mean_longitudinal_pressure_difference::compute_mean_longitudinal_pressure_difference

(

double precision, dimension(nx,ny,nz)

pressure

)

The mean mean pressure difference is defined as:

\[ \frac {1} {S_{left}} \int_{S_{left}}^{} p \, \mathrm{d}s - \frac {1} {S_{right}} \int_{S_{right}}^{} p \, \mathrm{d}s \]

compute_min_field()

double precision function mod_compute_min_field::compute_min_field

(

double precision, dimension(:,:,:), intent(in)

field

)

The argument field is defined on cells or faces.

compute_min_max_mean_velocity_magnitude()

subroutine mod_compute_min_max_mean_velocity_magnitude::compute_min_max_mean_velocity_magnitude	(	type(t_face_field), intent(in)	velocity,
		double precision, intent(out), optional	min_velocity,
		double precision, intent(out), optional	max_velocity,
		double precision, intent(out), optional	mean_velocity )

The mean velocity magnitude is defined as:

\[ \frac {1} {V} \int_{V}^{} \sqrt{(u^2 + v^2 + w^2)} \, \mathrm{d}v \]

where the velocity components are interpolated on cell centers.

compute_physical_diagnostics()

subroutine mod_compute_physical_diagnostics::compute_physical_diagnostics

This routine is called in the time loop.

compute_q_criterion()

subroutine mod_compute_q_criterion::compute_q_criterion	(	double precision, dimension(:,:,:), intent(out)	q_criterion,
		type(t_face_field), intent(in)	velocity )

The Q-criterion is given the second invariant of the velocity gradient tensor \( \mathbf{L} = \mathbf{\nabla} \mathbf{u} \). The definition follows :

\[ Q = \frac{1}{2}\Big( \mathrm{Trace}(\mathbf{L})^2 - \mathrm{Trace}(\mathbf{L}\cdot\mathbf{L}) \Big) = \frac{1}{2}\Big( (\mathbf{\nabla} \cdot \mathbf{u})^2 - \mathbf{\nabla} \mathbf{u} : \mathbf{\nabla} \mathbf{u}^T \Big) = \frac{1}{2}\Big( (\mathbf{\nabla} \cdot \mathbf{u})^2 + ||\mathbf{\Omega}||_F^2 - ||\mathbf{S}||_F^2 \Big) \]

where \(||.||_F\) denotes the Frobenius norm, \(\mathbf{S}\) is the symmetric strain rate tensor, and \(\mathbf{\Omega}\) is the skew-symmetric strain rate (or vorticity) tensor. The \(||\mathbf{S}||_F\) term is computed as described in the [compute_strain_rate_magnitude] routine, while the vorticity tensor norm develops as follows:

\[ ||\mathbf{\Omega}||_F^2 = \frac{1}{2} \Big( \Big( \underbrace{ \frac {\partial u} {\partial y}}_{\mathrm{I}} - \underbrace{ \frac {\partial v} {\partial x}}_{\mathrm{II}} \Big)^2 + \Big( \frac {\partial u} {\partial z} - \frac {\partial w} {\partial x} \Big)^2 + \Big( \frac {\partial v} {\partial z} - \frac {\partial w} {\partial y} \Big)^2 \Big) \]

The terms \(\mathrm{I}\) and \(\mathrm{II}\) are discretized using second-order finites differences as follows:

                   ┌─────┐
                   │     │
                   >     >  jp
                   │     │
                   ┢━━━━━┪        ┌──∧──┲━━∧━━┱──∧──┐  jp
    y              ┃     ┃        │     ┃     ┃     │
    │              >  o  >  j     │     ┃  o  ┃     │  j
    o──x           ┃     ┃        │     ┃     ┃     │
   z               ┡━━━━━┩        └──∧──┺━━∧━━┹──∧──┘  j
                   │     │
                   >     >  jm       im    i     ip
                   │     │
                   └─────┘
                   i  i  ip

\[ \Big( \frac {\partial u} {\partial y} \Big)_{\mathtt{i,j}} = \frac {1} {4} \Big( \mathtt{ \frac {u_{ip,jp} - u_{ip,j} } {\Delta yv_{jp}} } + \mathtt{ \frac {u_{i,jp} - u_{i,j} } {\Delta yv_{jp}} } + \mathtt{ \frac {u_{ip,j} - u_{ip,jm}} {\Delta yv_{j}} } + \mathtt{ \frac {u_{i,j} - u_{i,jm} } {\Delta yv_{j}} } \Big) \\ \Big( \frac {\partial v} {\partial x} \Big)_{\mathtt{i,j}} = \frac {1} {4} \Big( \mathtt{ \frac {v_{ip,jp} - v_{i,jp}} {\Delta xu_{ip}} } + \mathtt{ \frac {v_{ip,j} - v_{i,j} } {\Delta xu_{ip}} } + \mathtt{ \frac {v_{i,jp} - v_{im,jm}} {\Delta xu_{i}} } + \mathtt{ \frac {v_{i,j} - v_{im,j} } {\Delta xu_{i}} } \Big) \, . \]

The computations of the other partial derivatives terms are given in the [compute_strain_rate_magnitude] routine documentation.

compute_spatial_averaged_field()

double precision function mod_compute_spatial_averaged_field::compute_spatial_averaged_field	(	double precision, dimension(:,:,:), intent(in)	field,
		double precision, dimension(:,:,:), intent(in)	grid_volume )

The argument field can be defined on cells or faces.

compute_stefan_problem_front_position_melting()

double precision function mod_compute_stefan_problem_front_position_melting::compute_stefan_problem_front_position_melting

The front position of the Stefan problem is computed.

compute_stefan_problem_front_position_solidification()

double precision function mod_compute_stefan_problem_front_position_solidification::compute_stefan_problem_front_position_solidification

The front position of the Stefan problem is computed.

compute_strain_rate_magnitude()

subroutine mod_compute_strain_rate_magnitude::compute_strain_rate_magnitude	(	double precision, dimension(:,:,:), intent(out)	strain_rate_magnitude,
		type(t_face_field), intent(in)	velocity )

The strain rate magnitude is the Frobenius norm of the strain rate tensor \( || \mathbf{S} ||_F = \sqrt{\mathbf{S}:\mathbf{S}} \), which develops as:

\[ || \mathbf{S} ||_F ^2 = \Big( \underbrace{ \frac {\partial u} {\partial x}}_{\text{I}} \Big)^2 + \Big( \underbrace{ \frac {\partial v} {\partial y}}_{\text{II}} \Big)^2 + \Big( \underbrace{ \frac {\partial w} {\partial z}}_{\text{III}} \Big)^2 + \frac{1}{2}\Big( \underbrace{ \frac {\partial u} {\partial y} + \frac {\partial v} {\partial x}}_{\text{IV}} \Big)^2 + \frac{1}{2}\Big( \underbrace{ \frac {\partial u} {\partial z} + \frac {\partial w} {\partial x}}_{\text{V}} \Big)^2 + \frac{1}{2}\Big( \underbrace{ \frac {\partial v} {\partial z} + \frac {\partial w} {\partial y}}_{\text{VI}} \Big)^2 \]

where the derivatives are discretized using second-order finite differences as presented below. We use the following definition of the symmetric strain rate tensor \( \mathbf{S} = (\mathbf{\nabla} \mathbf{u} +\mathbf{\nabla} \mathbf{u}^T) / 2 \)

Term \(\text{I}\)

                   ┏━━━━━┓
    y              ┃     ┃
    │              >  o  >
    o──x           ┃     ┃
   z               ┗━━━━━┛
                   i  i  ip

\[ \Big( \frac {\partial u} {\partial x} \Big)_{\mathtt{i,j}} = \mathtt{ \frac {u_{ip,j} - u_{i,j} } {\Delta x_{i}} } \]

Terms \(\text{II}\) and \(\text{III}\) are computed likewise.

Term \(\text{IV}\)

                ┌─────┐
                │     │
             jp >     >
                │     │
                ┢━━━━━┪     jp ┌──∧──┲━━∧━━┱──∧──┐
    y           ┃     ┃        │     ┃     ┃     │
    │         j >  o  >      j │     ┃  o  ┃     │
    o──x        ┃     ┃        │     ┃     ┃     │
   z            ┡━━━━━┩      j └──∧──┺━━∧━━┹──∧──┘
                │     │           im    i     ip
             jm >     >
                │     │
                └─────┘
                i  i  ip

\[ \Big( \frac {\partial u} {\partial y} \Big)_{\mathtt{i,j}} = \frac {1} {4} \Big( \mathtt{ \frac {u_{ip,jp} - u_{ip,j} } {\Delta yv_{jp}} } + \mathtt{ \frac {u_{i,jp} - u_{i,j} } {\Delta yv_{jp}} } + \mathtt{ \frac {u_{ip,j} - u_{ip,jm}} {\Delta yv_{j}} } + \mathtt{ \frac {u_{i,j} - u_{i,jm} } {\Delta yv_{j}} } \Big) \\ \Big( \frac {\partial v} {\partial x} \Big)_{\mathtt{i,j}} = \frac {1} {4} \Big( \mathtt{ \frac {v_{ip,jp} - v_{i,jp}} {\Delta xu_{ip}} } + \mathtt{ \frac {v_{ip,j} - v_{i,j} } {\Delta xu_{ip}} } + \mathtt{ \frac {v_{i,jp} - v_{im,jm}} {\Delta xu_{i}} } + \mathtt{ \frac {v_{i,j} - v_{im,j} } {\Delta xu_{i}} } \Big) \]

compute_turbulent_kinetic_energy()

subroutine mod_compute_turbulent_kinetic_energy::compute_turbulent_kinetic_energy	(	double precision, dimension(:,:,:), intent(out)	turbulent_kinetic_energy,
		type(t_face_field), intent(in)	velocity_rms )

The turbulent_kinetic energy field for incompressible flows is defined as:

\[ \frac {1} {2} (\overline{u^{'}u^{'}} + \overline{v^{'}v^{'}} + \overline{w^{'}w^{'}}) \]

or for compressible flows:

\[ \frac {1} {2} (\widetilde{u^{'}u^{'}} + \widetilde{v^{'}v^{'}} + \widetilde{w^{'}w^{'}}) \]

where the velocity components are interpolated on cell centers. \( \widetilde{u^{'}u^{'}} \) is the Favre averaging and \( \overline{u^{'}u^{'}} \) is the Reynolds averaging

compute_viscous_dissipation_rate()

subroutine mod_compute_viscous_dissipation_rate::compute_viscous_dissipation_rate	(	double precision, dimension(:,:,:), intent(out)	viscous_dissipation_rate,
		type(t_face_field), intent(in)	velocity,
		double precision, dimension(:,:,:), intent(in)	viscosity )

Compute viscous energy dissipation rate per unit time and volume \( \Phi_d \)

\[ \Phi_d = 2 \mu S_{ij} S_{ij} - \frac{2\mu}{3} \Big(\frac{\partial u_k}{\partial x_k}\Big)^2 \, , \]

with \( S_{ij} \) the strain rate tensor. The developed formula implemented is

\[ \Phi_d = 2\mu \Big( \underbrace{ \frac {\partial u} {\partial x}}_{\text{I}} \Big)^2 + 2\mu \Big( \underbrace{ \frac {\partial v} {\partial y}}_{\text{II}} \Big)^2 + 2\mu \Big( \underbrace{ \frac {\partial w} {\partial z}}_{\text{III}} \Big)^2 +\mu \Big( \underbrace{ \frac {\partial u} {\partial y} + \frac {\partial v} {\partial x}}_{\text{IV}} \Big)^2 +\mu \Big( \underbrace{ \frac {\partial u} {\partial z} + \frac {\partial w} {\partial x}}_{\text{V}} \Big)^2 +\mu \Big( \underbrace{ \frac {\partial v} {\partial z} + \frac {\partial w} {\partial y}}_{\text{VI}} \Big)^2 -\frac {2\mu} {3} \Big( \underbrace{ \frac {\partial u} {\partial x}}_{\text{I}} + \underbrace{ \frac {\partial v} {\partial y}}_{\text{II}} + \underbrace{ \frac {\partial w} {\partial z}}_{\text{III}} \Big)^2 \]

where the derivatives are discretized using second-order finite differences as presented below. The routine call 'compute_strain_rate_magnitude'.

Term \(\text{I}\)

                   ┏━━━━━┓
    y              ┃     ┃
    │              >  o  >
    o──x           ┃     ┃
   z               ┗━━━━━┛
                   i  i  ip

\[ \Big( \frac {\partial u} {\partial x} \Big)_{\mathtt{i,j}} = \mathtt{ \frac {u_{ip,j} - u_{i,j} } {\Delta x_{i}} } \]

Terms \(\text{II}\) and \(\text{III}\) are computed likewise.

compute_vorticity_magnitude()

subroutine mod_compute_vorticity_magnitude::compute_vorticity_magnitude	(	double precision, dimension(:,:,:), intent(out)	vorticity_magnitude,
		type(t_face_field), intent(in)	velocity )

\[ \mathbf{\Omega}^2 = + \Big( \underbrace{ \frac {\partial u} {\partial y} - \frac {\partial v} {\partial x}}_{\text{IV}} \Big)^2 + \Big( \underbrace{ \frac {\partial u} {\partial z} - \frac {\partial w} {\partial x}}_{\text{V}} \Big)^2 + \Big( \underbrace{ \frac {\partial v} {\partial z} - \frac {\partial w} {\partial y}}_{\text{VI}} \Big)^2 \]

Term \(\text{IV}\)

                ┌─────┐
                │     │
             jp >     >
                │     │
                ┢━━━━━┪     jp ┌──∧──┲━━∧━━┱──∧──┐
    y           ┃     ┃        │     ┃     ┃     │
    │         j >  o  >      j │     ┃  o  ┃     │
    o──x        ┃     ┃        │     ┃     ┃     │
   z            ┡━━━━━┩      j └──∧──┺━━∧━━┹──∧──┘
                │     │           im    i     ip
             jm >     >
                │     │
                └─────┘
                i  i  ip

\[ \Big( \frac {\partial u} {\partial y} \Big)_{\mathtt{i,j}} = \frac {1} {4} \Big( \mathtt{ \frac {u_{ip,jp} - u_{ip,j} } {\Delta yv_{jp}} } + \mathtt{ \frac {u_{i,jp} - u_{i,j} } {\Delta yv_{jp}} } + \mathtt{ \frac {u_{ip,j} - u_{ip,jm}} {\Delta yv_{j}} } + \mathtt{ \frac {u_{i,j} - u_{i,jm} } {\Delta yv_{j}} } \Big) \\ \Big( \frac {\partial v} {\partial x} \Big)_{\mathtt{i,j}} = \frac {1} {4} \Big( \mathtt{ \frac {v_{ip,jp} - v_{i,jp}} {\Delta xu_{ip}} } + \mathtt{ \frac {v_{ip,j} - v_{i,j} } {\Delta xu_{ip}} } + \mathtt{ \frac {v_{i,jp} - v_{im,jm}} {\Delta xu_{i}} } + \mathtt{ \frac {v_{i,j} - v_{im,j} } {\Delta xu_{i}} } \Big) \]

compute_wall_nusselt_sherwood_number()

subroutine mod_compute_wall_nusselt_sherwood_number::compute_wall_nusselt_sherwood_number	(	double precision, dimension(:,:,:), intent(in)	temperature,
		type(t_boundary_condition), intent(in)	energy_bc,
		double precision	diagnostic_delta,
		double precision, intent(out)	left_nusselt,
		double precision, intent(out)	right_nusselt,
		double precision, intent(out)	bottom_nusselt,
		double precision, intent(out)	top_nusselt,
		double precision, intent(out)	back_nusselt,
		double precision, intent(out)	front_nusselt,
		double precision, intent(out), optional	ib_wall_nusselt_left,
		double precision, intent(out), optional	ib_wall_nusselt_right,
		double precision, intent(out), optional	ib_wall_nusselt_bottom,
		double precision, intent(out), optional	ib_wall_nusselt_top,
		double precision, intent(out), optional	ib_wall_nusselt_back,
		double precision, intent(out), optional	ib_wall_nusselt_front,
		double precision, intent(out), optional	ib_wall_nusselt )

This routine compute the wall Nusselt number for each boundary where a Dirichlet boundary condition is set on temperature.

If H is the length of the boundary, L the distance between two oppositve boundaries, \( \Delta T \) the temperature difference between these two walls, the Nusselt number is given by:

\[ Nu = \frac {1} {\lambda_0 H} \int_{0}^{H} \frac {\lambda\frac {\partial T} {\partial n}} {\frac {\Delta T} {L}} \, \mathrm{d}x \]

The normal derivative is computed at second order in space thanks to the compute_uncentered_boundary_cell_gradient routine.

The mass transfer analog of the Nusselt number is the Sherwood number so this routine can be used for both computation (see compute_physical_diagnostics.f90)