Diagnostic quantities

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. More...
double precision function	mod_compute_stefan_problem_front_position_melting::compute_stefan_problem_front_position_melting ()
	Compute front position. More...
double precision function	mod_compute_stefan_problem_front_position_solidification::compute_stefan_problem_front_position_solidification ()
	Compute front position. More...
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_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. More...
subroutine	mod_compute_kinetic_energy_per_mass_unit::compute_kinetic_energy_per_mass_unit (kinetic_energy, velocity)
	Compute the kinetic energy field on cells. More...
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. More...
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. More...
double precision function	mod_compute_mean_kinetic_energy::compute_mean_kinetic_energy (velocity)
	Compute the mean kinetic energy on the whole domain. More...
double precision function	mod_compute_mean_kinetic_energy::compute_kinetic_energy_integral (velocity)
	Compute the mean kinetic energy on the whole domain. More...
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. More...
double precision function	mod_compute_min_field::compute_min_field (field)
	Compute the min value of a field restricted to the inner domain. More...
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. More...
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. More...
subroutine	mod_compute_physical_diagnostics::compute_physical_diagnostics ()
	Execute enabled physical tools that returns a scalar value (Nusselt, mean quantities, etc.). More...
subroutine	mod_compute_q_criterion::compute_q_criterion (q_criterion, velocity)
	Compute Q-criterion of velocity field. More...
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. More...
subroutine	mod_compute_strain_rate_magnitude::compute_strain_rate_magnitude (strain_rate_magnitude, velocity)
	Compute strain rate magnitude of the velocity field. More...
subroutine	mod_compute_u_plus::compute_u_plus ()
	Compute u+ (physical and immersed boundaries)
subroutine	mod_compute_viscous_energy_dissipation_rate::compute_viscous_energy_dissipation_rate (viscous_energy_dissipation_rate, velocity, viscosity)
	Compute viscous energy dissipation rate per unit time and volume \( \Phi_d \). More...
subroutine	mod_compute_vorticity_magnitude::compute_vorticity_magnitude (vorticity_magnitude, velocity)
	Compute the vorticity magnitude: More...
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. More...
subroutine	mod_compute_wall_shear_stress::compute_wall_shear_stress (velocity, wall_shear_stress_bottom, wall_shear_stress_top)
	Computes the wall skin friction coefficient 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_initial_physical_diagnostics()

subroutine mod_compute_physical_diagnostics::compute_initial_physical_diagnostics

This routine is called before the time loop.

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

\[ Q = \frac{1}{2} \Big( |\mathbf{\Omega}|^2 - |\mathbf{S}|^2 \Big) \]

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

\[ |\mathbf{\Omega}|^2 = \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 \]

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 other terms are discretized likewise.

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} | \), which develops as:

\[ \mathbf{S}^2 = 2 \Big( \underbrace{ \frac {\partial u} {\partial x}}_{\text{I}} \Big)^2 + 2 \Big( \underbrace{ \frac {\partial v} {\partial y}}_{\text{II}} \Big)^2 + 2 \Big( \underbrace{ \frac {\partial w} {\partial z}}_{\text{III}} \Big)^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 \]

where the derivatives are discretized using second-order finite differences as presented below.

Term @f$\text{I}@f$

                   ┏━━━━━┓
    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 @f$\text{IV}@f$

                ┌─────┐
                │     │
             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_viscous_energy_dissipation_rate()

subroutine mod_compute_viscous_energy_dissipation_rate::compute_viscous_energy_dissipation_rate	(	double precision, dimension(:,:,:), intent(out)	viscous_energy_dissipation_rate,
		type(t_face_field), intent(in)	velocity,
		double precision, dimension(:,:,:), intent(in)	viscosity
	)

\[ \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.

Term @f$\text{I}@f$

                   ┏━━━━━┓
    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 @f$\text{IV}@f$

                ┌─────┐
                │     │
             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_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 @f$\text{IV}@f$

                ┌─────┐
                │     │
             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)