%-------------------------------------------------------------------------------

% This file is part of Code_Saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2018 EDF S.A.
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA 02110-1301, USA.
%-------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
\section{Introduction}

\hypertarget{boundary}{}

Boundary conditions are required in at least three main cases:

\begin{itemize}
\item calculation of the convection terms (first order derivative in space) at
the boundary: the code uses a mass flux at the boundary and requires the
value of the convected variable when the flow is entering into
the domain (or more general wave relations in the sense of the characteristic curves of the system entering the domain);
\item calculation of the diffusion terms (second order derivative
in space): the code needs
a method to determine the value of the first order spatial derivatives
at the boundary, these define \emph{e.g.} the stresses or the thermal fluxes at the wall;
\item calculation of the cell  gradients: the variables at the boundary faces
 allow the code to define the gradient inside the cell connected to the boundary (\emph{e.g.} the pressure gradient or
the transpose gradient terms in the stress-strain relation).
\end{itemize}

These considerations only concern the computational field variables
(velocity, pressure, Reynolds tensor,
scalars solution of a advection-diffusion equations \emph{etc.}). For these variables
\footnote{
The other variables
(the physical properties for instance) have a different treatment which will
not be detailed here (for instance, for the density, the user defines
directly the values at the boundary. This information is then stored; one
is referred to \fort{cs\_user\_physical\_properties} (see the \doxygenfile{cs__user__physical__properties_8f90.html}{programmers reference of the dedicated subroutine})
 or \fort{phyvar} for more information).
},
the user has to define the boundary conditions at every boundary face.

The boundary conditions could be of Neumann type (when the flux is imposed)
or Dirichlet type (when the value of the field variable is prescribed), or mixed type, also
called Robin type (when a combination linking the field variable to its derivative flux
is imposed).

The code (see the \doxygenfile{condli_8f90.html}{programmers reference of the dedicated subroutine})
 transforms the boundary conditions provided by the user
 into two internal formats of representation of the boundary
 conditions.

 A particular treatment is performed on walls:
 wall functions are used to model the boundary layer flow in the vicinity of the wall when the mesh is too coarse to correctly capture the sharp variations of the fields with a linear profile in the near wall cell.
 This will be detailed in the next sections.

A particular treatment on symmetry boundaries is also performed for vectors and tensors
whereas a symmetry boundary is equivalent to an homogeneous Neumann condition (zero normal gradient) for scalar fields.

The physics model that the user wishes to apply needs to be translated into pairs of coefficients
entering the linear system of equations that the code will solve. For any variable $\varia$ for every boundary faces $\fib$ these coefficients are:

\begin{itemize}
\item $\left( A^g_\fib , \, B^g_\fib\right)$ used by the gradient operator and by the advection operator.
The value at the boundary face $\fib$ of the variable $\varia$ is
then defined as:
\begin{equation}\label{eq:bndcnd:coef_g_def}
\varia_\fib \equiv A^g_\fib + B^g_\fib \varia_\centip.
\end{equation}
%
\item $\left( A^f_\ib , \, B^f_\ib\right)$ used by the diffusion operator.
The value at the boundary face $\fib$ of the diffusive flux $q_\ib$ of the variable $\varia$
is then defined as (see Equation \eqref{eq:spadis:boundary_flux}):
\begin{equation}\label{eq:bndcnd:coef_f_def}
q_\ib \equiv \dfrac{D_\ib \left( K_\fib , \, \varia \right)}{\norm{\vect{S}}_\fib} = -\left(A^f_\ib + B^f_\ib  \varia_\centip\right).
\end{equation}
Note that the diffusive boundary coefficients are oriented, which means that they are such that if $D_\ib \left( K_\fib , \, \varia \right) $ is positive,
this flux is gained by $\varia_\celli$.
 %
\end{itemize}
The definitions are recalled on \figurename~\ref{fig:bndcnd:boundary_flux}.

\begin{figure}[!htpb]
\centering
\includegraphics[width=0.8\textwidth]{fluxbord}
\caption{Boundary cell.\label{fig:bndcnd:boundary_flux}}
\end{figure}

\begin{example}\label{sec:bndcnd:intro}
The heat flux $q_w$ from a solid wall to a laminar flow is defined and discretised as
\begin{equation}
q_w = \lambda \grad T \cdot n_\ib = - \lambda \dfrac{T_\centip -T_\centf}{ \overline{\centip \centf}} = -h_{int} \left(T_\centip - T_\centf \right),
\end{equation}
where $h_{int} / \overline{\centip \centf}$ is the exchange coefficient between the wall and the fluid point $\centip$.

\begin{itemize}
\item If the wall temperature $T_w$ is known, $T_\centf =T_w$ is set as Dirichlet condition, $q_w$
will be a result of the simulation, and the gradient and diffusive flux coefficients
are as follows:
\begin{equation}
\begin{array}{r c l}
\left\lbrace
\begin{array}{r c l}
A^g_\fib & = & T_w , \\
B^g_\fib & = & 0,
\end{array}
\right.
 & &
\left\lbrace
\begin{array}{r c l}
A^f_\ib & = & -h_{int} T_w, \\
B^f_\ib & = & h_{int}.
\end{array}
\right.
\end{array}
\end{equation}

\item If the wall heat flux $q_w$ is known, this is a Neumann condition and $T_\centf = T_w$
will be a result of the simulation and the gradient and diffusive flux coefficients
are as follows:
\begin{equation}
\begin{array}{r c l}
\left\lbrace
\begin{array}{r c l}
A^g_\fib & = & - \dfrac{q_w}{h_{int}} ,\\
B^g_\fib & = & 1,
\end{array}
\right.
& &
\left\lbrace
\begin{array}{r c l}
A^f_\ib & = & q_w ,\\
B^f_\ib & = & 0 .
\end{array}
\right.
\end{array}
\end{equation}

\end{itemize}
\end{example}

%-------------------------------------------------------------------------------
\section{Standard user boundary conditions}

The user generally gives standard boundary conditions, which are:

\begin{description}
\item[Inlet:] it corresponds to a Dirichlet boundary condition on all the transported variables
(and should therefore be given by the user)
and to a homogeneous Neumann on the pressure field.

\item[Outlet:] it correspond to a homogeneous Neumann boundary condition on all the transported variables.
For the pressure field, a Dirichlet boundary condition which is expected to mimic $\displaystyle\frac{\partial^2 P}{\partial n \partial \tau}=0$ for any vector  $\vect{\tau}$ parallel to the outlet face
(see \S \ref{sec:bndcnd:outlet_pressure} for more details).
This condition means that the pressure
profile does not vary in the normal direction of the outlet. Warning: if the outgoing mass-flux is negative,
\emph{i.e.} if the outlet becomes an inlet, then the mass-flux is clipped to zero. Moreover, since
 the pressure field is defined up to a constant, it is fixed to a reference pressure $P_0$
at an arbitrary chosen outlet boundary face.
The user can choose an other desired face where a Dirichlet on pressure is prescribed.

\item[Free inlet/outlet:] it corresponds to the standard outlet when the flow is outgoing (see \S \ref{sec:bndcnd:outlet_pressure}), but to a free inlet when the flow
is ingoing. The Bernoulli relationship is used to derive a boundary condition on the pressure increment to couple velocity and pressure at
these free inlet faces.
Note that no clipping on the velocity is imposed.
The same boundary conditions as for outlet on the other variables is imposed. For more details please refer to \S \ref{sec:bndcnd:inlet_bernoulli_pressure}.


\item[Walls:] This particular treatment will be detailed in the following sections.
For the velocity, the aim is to transform the Dirichlet boundary condition (the velocity at the wall
is equal to zero, or the velocity of a moving wall) into a Neumann boundary condition where the wall shear stress
is imposed function of the local flow velocity and the turbulence characteristics.
A similar treatment using wall functions is done on every transported variable if this variable is prescribed.
The boundary condition on the pressure field is a homogeneous Neumann by default, or alternatively an extrapolation of the gradient.

\item[Symmetries:] This condition corresponds to a homogeneous Neumann for the scalar fields
(\emph{e.g.} the pressure field or the temperature field). For vectors, such as the velocity, it corresponds
to impose a zero Dirichlet on the component normal to the boundary, and a homogeneous Neumann
on the tangential components. Thus, this condition couples the  vector components if the symmetry faces are not
aligned with the reference frame. The boundary condition for tensors, such as the Reynolds stresses, will be detailed in the following sections.
\end{description}

See the \doxygenfile{cs__user__boundary__conditions_8f90.html}{programmers reference of the dedicated subroutine} for further details.

%-------------------------------------------------------------------------------
\section{Internal coding of the boundary conditions -- Discretization}

As already mentioned, the boundary conditions set by the user for the variable $\varia$
are translated into two pairs of coefficients $\left( A^g_\fib , \, B^g_\fib\right)$ which are used by the gradient operator and by the advection operator, and $\left( A^f_\ib , \, B^f_\ib\right)$ for use by the diffusion operator, this for all the boundary faces $\fib$.

Let us first recall the general form of the transport equation of a variable $\varia$, which could be a scalar,
a vector or a tensor:
\begin{equation}\label{eq:bndcnd:gov_eqn_scalar}
\displaystyle C \rho \frac{\partial \varia}{\partial t} + C \grad \varia \cdot \left( \rho \vect{u} \right) = \dive \left(\displaystyle K \, \grad \varia \right) +ST_\varia .
\end{equation}

In the Equation (\ref{eq:bndcnd:gov_eqn_scalar})
$\rho$ is the density of the fluid, $\left( \rho \vect{u} \right)$ the convective mass flux of the variable $\varia$, $K$ its
conductivity or diffusivity and $S$ any additional source terms.
Note that $K$ is the sum of molecular and turbulent diffusivity in case of RANS modelling with an eddy viscosity model.
The dimension of $K$ for different variables is displayed in \tablename~\ref{tab:bdncnd:diffusivity}.
The value of $C$ is $1$ for all the variables except for the temperature where $C$ is the specific heat $C_p$.
If the variable $\varia$ is the variance of another scalar, then its diffusivity
is deduced from the scalar itself.

\begin{table}
{\scriptsize
\begin{center}
\begin{tabular}{||l|l|l||l|l|l||}
\hline
\multicolumn{3}{||c||}{$\varia$}&\multicolumn{3}{|c||}{$K$}\\
\hline
symbol                            & name                              & units                    &
symbol                            & name                              & units                  \\
\hline
$u_i$                             & velocity                          & $m.s^{-1}$                &
$\mu$ or $\mu+\mu_t$              & dynamic viscosity                 & $kg.m^{-1}.s^{-1}$     \\
$P$                               & pressure                          & $kg.m^{-1}.s^{-2}$       &
$\Delta t$                     & time step                         & $s$                      \\
$T$                               & temperature                       & $K$                        &
$\lambda$                    & thermal conductivity              & $kg.m.s^{-3}.K^{-1}$ \\
                                  &                                   &                            &
                                  &                                   & $=W.m^{-1}.K^{-1}$\\
$h$                               & enthalpy                          & $m^{2}.s^{-2}$&
$\lambda/C_p$              & thermal conductivity over specific heat  & $kg.m^{-1}.s^{-1}$     \\
                                  &                                   & $=J.kg^{-1}$&
                                  &                                   &                          \\
$\varia$                     & variable                          & unit of ($\varia$)               &
$K  $                          & conductivity or diffusivity       & $kg.m^{-1}.s^{-1}$     \\
\hline
\end{tabular}
\end{center}
}
\caption{Definitions and units of $K$ for the most common equations.}\label{tab:bdncnd:diffusivity}
\end{table}

%-------------------------------------------------------------------------------
\subsection{Basic Dirichlet boundary conditions}\label{sec:bndcnd:dirichlet_bc}

Imposing a basic Dirichlet condition $\varia^{imp}_\fib$ on a boundary face $\fib$ is translated into:

\begin{equation}
\begin{array}{r c l}
\left\lbrace
\begin{array}{r c l}
A^g_\fib & = &\varia^{imp}_\fib , \\
B^g_\fib & = & 0,
\end{array}
\right.
 & &
\left\lbrace
\begin{array}{r c l}
A^f_\ib & = & -h_{int} \varia^{imp}_\fib, \\
B^f_\ib & = & h_{int}.
\end{array}
\right.
\end{array}
\end{equation}

The term $h_{int}$ is a internal coding coefficient (automatically provided by the code) similar to an exchange coefficient. Its value for particular variables is
given in \tablename~\ref{tab:bndcnd:hint_phi_condli}.

\begin{remark}
$\varia^{imp}_\fib$ must be specified by the user. The boundary type code is $1$ (see \tablename~\ref{tab:ICODCLadm_condli}).
\end{remark}


%-------------------------------------------------------------------------------
\subsection{Neumann boundary conditions}

Imposing a Neumann condition $D^{imp}_\ib$ on a boundary face $\fib$
 is translated into

\begin{equation}
\begin{array}{r c l}
\left\lbrace
\begin{array}{r c l}
A^g_\fib & = & - \dfrac{D^{imp}_\ib}{h_{int}} ,\\
B^g_\fib & = & 1,
\end{array}
\right.
& &
\left\lbrace
\begin{array}{r c l}
A^f_\ib & = & D^{imp}_\ib ,\\
B^f_\ib & = & 0 .
\end{array}
\right.
\end{array}
\end{equation}

\begin{remark}
$D^{imp}_\ib$ must be specified by the user. The boundary type code is $3$ (see \tablename~\ref{tab:ICODCLadm_condli}).
\end{remark}

%-------------------------------------------------------------------------------
\subsection{Mixed or Robin boundary conditions}
As explained in the introduction \ref{sec:bndcnd:intro}, the simple case of a heat flux $q_w$ from a solid wall reads: $q_w = \lambda \grad T \cdot \vect{n}_\ib = - \lambda \dfrac{T_\centip - T_\centf}{\overline{\centip \centf}} = -h_{int} \left( T_\centip - T_\centf \right)$ where $h_{int} = \dfrac{\lambda}{\overline{\centip \centf}}$ is displayed in \tablename~ \ref{tab:bndcnd:hint_phi_condli}.
As presented in \S~ \ref{sec:bndcnd:dirichlet_bc}, Dirichlet condition $T_\centf = T_w$ simply leads to $A^f = - h_{int} T_w$ and $B^f = h_{int}$.
In some cases, heat transfer outside the flow domain is described by a one-dimensional conduction equation $q = -h_{ext} \left(T_{REF} - T_w \right)$
in which some reference temperature $T_{REF}$ is given and the coefficient $h_{ext}$ is known from an analytical or experimental correlation.
Equaling this to the heat flux computed in the first fluid cell at the boundary gives $ h_{ext} \left(T_{REF} - T_w \right) = h_{int} \left( T_\centip - T_\centf \right)$. This relation between the variable and its gradient at the boundary is called a Robin or mixed boundary condition.

More generally, to impose an external Dirichlet $\varia^{imp, \, ext}$ at some distance outside the boundary face,
this is done by a user prescribed external coefficient $h_{ext}$
(see \figurename~\ref{fig:bndcnd:boundary_flux}), and the boundary condition coefficients then read:

\begin{equation}
\begin{array}{r c l}
\left\lbrace
\begin{array}{r c l}
A^g_\fib & = &\dfrac{h_{ext}}{h_{int}+h_{ext}} \varia^{imp, \, ext}, \\
B^g_\fib & = & \dfrac{h_{int}}{h_{int}+h_{ext}} ,
\end{array}
\right.
& &
\left\lbrace
\begin{array}{r c l}
A^f_\ib & = & -h_{eq} \varia^{imp, \, ext} , \\
B^f_\ib & = & h_{eq},
\end{array}
\right.
\end{array}
\end{equation}
where $h_{eq} $ is defined by $h_{eq}=\dfrac{h_{int} h_{ext}}{ h_{int} + h_{ext}}$.
The harmonic mean (as in \eqref{eq:spadis:harmonic_viscosity}) comes from summing of resistances instead of conductances.
Note that this case reduces to Dirichlet condition if $h_{ext}$ tends to the infinity.

\begin{remark}
Both $\varia^{imp, \, ext} $ and $ h_{ext} $ must be specified by the user. The boundary code is $1$ (see \tablename~\ref{tab:ICODCLadm_condli}). Take care that an outgoing flux is counted positively.
\end{remark}

%-------------------------------------------------------------------------------
\subsection{Convective outlet boundary conditions}\label{sec:bndcnd:convective_outlet}

If the user wishes to impose a convective outlet (also called radiative outlet) condition which reads:
\begin{equation}
\der{ \varia}{ t} + C \der{ \varia}{ n} = 0,
\end{equation}
where $C$ denotes the convective celerity, then the internal coding reads:

\begin{equation}
\begin{array}{r c l}
\left\lbrace
\begin{array}{r c l}
A^g_\fib & = & \dfrac{1}{1+ CFL} \varia^n_\fib ,\\
B^g_\fib & = & \dfrac{CFL}{1+ CFL},
\end{array}
\right.
& &
\left\lbrace
\begin{array}{r c l}
A^f_\ib & = &  -  \dfrac{h_{int}}{1+ CFL} \varia^n_\fib,\\
B^f_\ib & = &  \dfrac{h_{int}}{1+ CFL},
\end{array}
\right.
\end{array}
\end{equation}
where $CFL \equiv \dfrac{C \Delta t}{\overline{\centip \centf} }$.

\begin{remark}
Both $C$ and $\varia^n_\fib$ must be specified by the user, the boundary code is $2$ (see \tablename~\ref{tab:ICODCLadm_condli}).
\end{remark}

\begin{table}
%{\tiny
\begin{center}
\begin{tabular}{||l|l|l||l|l||}
\hline
\multicolumn{3}{||c||}{$\varia$} & \multicolumn{2}{c||}{$h_{int}$}     \\
\hline
symbol                        & name                        & units                    & homogeneous to                        & units           \\
\hline
$\vect{u}$                    & velocity                    & $m.s^{-1}$               &$\dfrac{\mu+\mu_t}{ \centip \centf}$   & $kg.m^{-2}.s^{-1}$       \\
$P$                           & pressure                    & $kg.m^{-1}.s^{-2}$       & $\dfrac{\Delta t}{ \centip \centf}$   & $s.m^{-1}$                \\
$T$                           & temperature                 & $K$                      &$\dfrac{\lambda+C_p\mu_t/\sigma_t}{ \centip \centf}$  &$kg.s^{-3}.K^{-1}$\\
                              &                             &                          &                                       & $W.m^{-2}.K^{-1}$\\
$h$                           & enthalpy                    & $m^{2}.s^{-2}$           &$\dfrac{\lambda+C_p\mu_t/\sigma_t}{ \centip \centf}$ &$kg.m^{-2}.s^{-1}$\\
                              &                             & $J.kg^{-1} $             &                                                            &                                  \\
$\varia$                      & scalar                      & unit of($\varia$)       &$\dfrac{K}{ \centip \centf}$      & $kg.m^{-2}.s^{-1}$       \\
\hline
\end{tabular}
\end{center}

\begin{center}
\begin{tabular}{||l|l|l||l|l||}
\hline
\multicolumn{3}{||c||}{$\varia$} & \multicolumn{2}{c||}{$D^{imp}$} \\
\hline
symbol     & name                & units         & homogeneous to             & units                                     \\
\hline
$\vect{u}$ & velocity            & $m.s^{-1}$         &$\left( (\mu+\mu_t)\,\gradv \vect{u} \right)\cdot \vect{n}$  & $kg.m^{-1}.s^{-2} $       \\
$p$        & pressure            & $kg.m^{-1}.s^{-2}$ &$\left( \Delta t \grad P \right) \cdot \vect{n}$             & $kg.m^{-2}.s^{-1}$                        \\
$T$        & temperature         & $K$                &$\left( (\lambda+C_p\mu_t/\sigma_t)\grad T\right) \cdot \vect{n} $ &$kg.s^{-3}$ \\
                 &                         &                                     &                                  &$W.m^{-2}$ \\
$h$          & enthalpy           & $m^{2}.s^{-2}$          &$\left( \lambda/C_p+\mu_t/\sigma_t) \grad H \right) \cdot \vect{n}$&$kg.s^{-3}$ \\
                 &                        & $J.kg^{-1}$               &                                                                                                               & $W.m^{-2}$  \\
$\varia$   & scalar               & unit of ($\varia$)       &$K \,\grad \varia \cdot \vect{n}$              & $kg.m^{-2}.s^{-1}.$ unit of ($\varia$)       \\
\hline
\end{tabular}
\end{center}
%}
\caption{Definitions and units of $h_{int}$ and $D^{imp}$ for the most common equations.}\label{tab:bndcnd:hint_phi_condli}
\end{table}


%-------------------------------------------------------------------------------
\section{Wall boundary conditions}\label{sec:clptur}

\hypertarget{wallboundary}{}

This section is dedicated to wall boundary conditions
 for the velocity, and any transported scalar (such as the temperature, the enthalpy, \emph{etc}.)
when a wall boundary value is prescribed and a wall function is set by the user. The dedicated boundary conditions of the turbulent variables ($k$, $\varepsilon$, $R_{ij}$) are detailed in \chaptername~ \ref{chapter:turbul}.
 Here is a link to the \doxygenfile{clptur_8f90.html}{programmers reference of the dedicated subroutine}.

Simplified balance analysis proves that the fluid velocity profile or the Temperature profile are
highly non-linear  for turbulent flows in the vicinity of walls.

The aim of wall functions is to model the vicinity of the wall by substituting the Dirichlet condition on the fluid velocity and on the scalars such as the temperature ($\vect{u}_\fib = \vect{v}_{wall}$, $T_\fib = T_{wall}$ where $\vect{v}_{wall}$ is the velocity of the wall and $T_{wall}$ is the temperature of the wall) by a Robin-type boundary condition linking the wall shear stress to the velocity of the fluid, and the
thermal wall flux to the temperature within the first cell.

%-------------------------------------------------------------------------------
\subsection{Velocity boundary condition for smooth walls and rough walls}
We assume it is projected onto the tangent plane to the wall (if it is not, then the code projects it).
The fluid velocity in the reference frame attached to the wall ("relative" velocity) projected to the wall writes
$\vect{u}^r _\centip \equiv \left(\tens{1} - \vect{n} \otimes \vect{n} \right)\cdot \left( \vect{u}_\centip - \vect{v}_{wall} \right)$.

The orthonormal coordinate system attached to the wall writes
$\hat {\mathcal R}=(\vect{t}, \, -\vect{n}, \, \vect{b})$:
%
\begin{itemize}
\item [$\bullet$] $\vect{t} = \displaystyle\frac{ \vect{u}^r_{\centip} }{ \norm{ \vect{u}^r_{\centip} }}$ is the unit vector parallel to the projection of the relative velocity at $\centip$, $\vect{u}^r_{\centip}$, in the plane tangent to the wall
 (\emph{i.e.} orthogonal to $\vect{n}$),
\item [$\bullet$] $-\vect{n}$ is the unit vector orthogonal to the wall and directed towards the interior of the computational domain,
\item[$\bullet$] $\vect{b}$ is the unit vector which completes the positively oriented coordinate system.
\end{itemize}
In this reference-frame the wall distance of the cell centre to the wall $\overline{\centip \centf}$ is denoted by $y$.

The objective of wall functions
is to increase the exchange coefficient
 to reflect the higher level of mixing effectively taking place in the near wall cell due to
turbulence and due to effectively much sharper gradients of the computed variable than the
linear profile assumed inside each cell by Finite Volumes discretisation. The end result (Cf \eqref{eq:bndcnd:boundary_condition_wall_function}) will
be to correct the laminar coefficient, $h_{int} = \lambda / y $, by a dimensionless analytical function (a
"wall function") $u^+$ that represents the realistic non-linear profile. The correction is moreover
proportional to the cell centre to wall dimensionless distance, $y^+$, over which a linear profile
assumption may or may not be appropriate.
These functions $u^+$ and $y^+$ that depends on the level of turbulent kinetic energy,
are made dimensionless by:
\begin{itemize}
\item either the wall shear stress $\tau_{wall}$ only (this is called "one scale friction velocity" wall function),
\item or by the wall shear stress $\tau_{wall}$ and the turbulent kinetic energy (this is called "two friction velocity scales").
\end{itemize}

When the mesh is made fine enough to capture the sharp variations of the variables near the
wall, wall functions are no longer needed, variables are simply given Dirichlet values at the
wall and the laminar exchange coefficient is correct. However the more basic versions of the
turbulence models ("high Reynolds" $k-\varepsilon$ or Reynolds stress transport) will not be able to
automatically "relaminarise" (\emph{i.e.} reduce the turbulent mixing effect to zero at the correct rate
of decay as y tends toward $0$) and erroneous predictions will result (generally too much friction
or heat transfer).
When using a refined near wall mesh, "down to the wall" or "low Reynolds" versions of the
turbulence models must be selected (\emph{e.g.} $v^2-f$ or elliptic blending models). However in this
case one must make sure the mesh is fine enough.

To allow use of standard high Reynolds models whatever the mesh density,
"scalable wall functions" can be activated. This consist in shifting slightly the
wall outside of the mesh if the first cell lies in the viscous sub-layer. "Scalable wall function"
only work with the "two friction velocity scales" version.

The wall functions are usually derived from Eddy Viscosity Models such as $k-\varepsilon$, thus the coming sections
assume that $\tens{R} = \dfrac{2}{3}k \tens{1} - 2\nu_T \deviator{\tens{S}}$.


%-------------------------------------------------------------------------------
\subsubsection{One friction velocity scale}\label{sec:bndcnd:1velocityscales}

The simplified momentum balance in the first boundary cell reads:
\begin{equation}\label{eq:bndcnd:simplified_momentum_evm}
\left( \mu + \mu_T \right) \dfrac{\partial u}{\partial y} = \tau_{wall}.
\end{equation}

Let $u^\star$ be the friction velocity defined by:
\begin{equation}\label{eq:bndcnd:ustar_1scale_def}
u^\star \equiv \sqrt{ \dfrac{\tau_{wall}}{\rho}},
\end{equation}
%
the equation \eqref{eq:bndcnd:simplified_momentum_evm} can be rewritten:
%
\begin{equation}\label{eq:bndcnd:simplified_momentum_evm_adim}
\left( 1 + \dfrac{\mu_T}{ \mu} \right) \dfrac{\partial u^+}{\partial y^+} = 1,
\end{equation}
where $u^+ \equiv \norm{ \vect{u}^r_{\centip} } / u^\star$ and $ y^+ \equiv y u^\star / \nu$.

The mixing length model states $\nu_T = L_m^2 S $, where the stain rate $S$ reduces to $ \frac{\partial u^r_{\centip} }{\partial y} $ in the wall vicinity.
The mixing length is proportional to the wall distance according to the Prandtl model
\begin{equation}\label{eq:bndcnd:prandtl_model}
L_m = \kappa y,
\end{equation}
where $\kappa = 0.42$ is the Karman constant.

Two areas are therefore defined, one where $ \mu_T/\mu  \ll 1$ called the viscous sub-layer where $\mu_T/\mu $ is neglected in \eqref{eq:bndcnd:simplified_momentum_evm_adim} the velocity profile is found to be linear:
%
\begin{equation}\label{eq:bndcnd:linear_profile}
u^+ = y^+,
\end{equation}
%
the other where $ \frac{\mu_T}{\mu}  \gg 1$ and the velocity profile becomes logarithmic:
\begin{equation}\label{eq:bndcnd:log_profile}
u^+ = \dfrac{1}{\kappa} \ln \left( y^+\right) + C_{log},
\end{equation}
where $C_{log}=5.2$.

\eqref{eq:bndcnd:linear_profile} is valid for $y^+<5$ and \eqref{eq:bndcnd:log_profile} for $y^+>30$ so there is a gap (corresponding to the so called buffer
layer) which more sophisticated models can cover, but they are not detailed here. Instead we
introduce a dimensionless
 limit distance which crudely separates the viscous sub-layer from the logarithmic region writes $y^+_{lim}$. Its value is $1/ \kappa$ in general (to ensure the continuity of the velocity gradient) and $10.88$ in LES (to ensure the continuity of the velocity).

The $u^\star$ is computed by iteratively solving\eqref{eq:bndcnd:linear_profile} or \eqref{eq:bndcnd:log_profile}:

\begin{equation}
\begin{array}{r c l c l}
u^\star &=& \sqrt{ \dfrac{\norm{ \vect{u}^r_{\centip} } \nu }{ y } } & \textrm{ if } & \dfrac{\norm{ y \, \vect{u}^r_{\centip} }}{\nu} < y^{+\, 2}_{lim}, \\
%
\left\lbrace
\begin{array}{r c l}
\left( u^\star \right)^{0} &=& \exp \left( - \kappa C_{log}\right) \dfrac{\nu}{y} \\
 \left( u^\star \right)^{q+1} &=& \dfrac{ \kappa \norm{ \vect{u}^r_{\centip} } + \left( u^\star \right) ^{q} }{ \ln \left( \dfrac{y \left( u^\star \right) ^{q} }{ \nu} \right) + \kappa C_{log} +1 }.
\end{array}
\right.
& \textrm{ otherwise}& \textrm{solve iteratively in $q$}
\end{array}
\end{equation}

Therefore, the value of $\frac{y^+}{u^+}$ in the two layers reads:
\begin{equation}\label{eq:bndcnd:ypup_1scale}
\left\lbrace
\begin{array}{r c l l}
\dfrac{y^+}{u^+} &=& 1 & \textrm{if } y^+ < y^+_{lim}, \\
\dfrac{y^+}{u^+} &=& \dfrac{y^+}{ \dfrac{1}{\kappa} \ln \left( y^+\right) + C_{log}} & \textrm{otherwise}.
\end{array}
\right.
\end{equation}

%-------------------------------------------------------------------------------
\subsubsection{Two friction velocity scales}\label{sec:bndcnd:2velocityscales}

The simplified momentum balance \eqref{eq:bndcnd:simplified_momentum_evm} still holds, but is non-dimensioned not only
by the wall shear stress but also the turbulent kinetic energy. More precisely, let $u_k = \sqrt{\sqrt{C_\mu} k }$ be the friction velocity based on
the turbulent kinetic energy in the first cell.

This definition is in fact valid only for high-$y^+$ values and a blending is performed in case of the intensity of turbulence is very low:
\begin{equation}\label{eq:def_uk}
u_k \equiv \sqrt{ g  \dfrac{\nu \norm{ \vect{u}^r_{\centip}}}{y} + \left(1-g\right)\sqrt{C_\mu}k },
\end{equation}
where $g$ is a blending factor defined by $g = \exp \left(- \dfrac{\sqrt{k} y}{11 \nu} \right)$ (see \cite{}).

The friction velocity $u^\star $ is now defined by:
\begin{equation}\label{eq:bndcnd:ustar_2scale2_def}
u^\star \equiv  \dfrac{\tau_{wall}}{\rho u_k}.
\end{equation}

Equation \eqref{eq:bndcnd:simplified_momentum_evm} can now be rewritten:
%
\begin{equation}\label{eq:bndcnd:simplified_momentum_evm_2scales}
\left( 1 + \dfrac{\mu_T}{ \mu} \right) \dfrac{\partial u^+}{\partial y^+_k} = 1,
\end{equation}
where $u^+ \equiv \frac{\norm{ \vect{u}^r_{\centip} }}{u^\star}$ but $ y^+_k \equiv \frac{y u_k}{ \nu}$.

Let us now remark that $\nu_T = C_\mu^{\frac{1}{4}}\sqrt{k}  C_\mu^{\frac{3}{4}} \dfrac{k^{\frac{3}{2}}}{\epsilon} = u_k L_m $ where
$L_m = C_\mu^{\frac{3}{4}} \dfrac{k^{\frac{3}{2}}}{\epsilon}$ is the integral length scale. Again using the Prandtl model \eqref{eq:bndcnd:prandtl_model},
the following equation is obtained:
%
\begin{equation}\label{eq:bndcnd:simplified_momentum_evm_adim_2scales}
\left( 1 + \kappa y_k^+ \right) \dfrac{\partial u^+}{\partial y^+_k} = 1.
\end{equation}

The two areas defined in the previous section \eqref{eq:bndcnd:ypup_1scale} that are the viscous sub-layer and the logarithmic layer still hold, but with $y^+_k$ instead of $y^+$:
%
\begin{equation}\label{eq:bndcnd:ypup_2scale}
\left\lbrace
\begin{array}{r c l l}
\dfrac{y^+_k}{u^+} &=& 1 & \textrm{if } y^+_k < y^+_{lim}, \\
\dfrac{y^+_k}{u^+} &=& \dfrac{y^+_k}{ \dfrac{1}{\kappa} \ln \left( y^+_k\right) + C_{log}} & \textrm{otherwise}.
\end{array}
\right.
\end{equation}

%-------------------------------------------------------------------------------
\subsubsection{Rough wall functions}\label{sec:bndcnd:2velocityscales}
When specifying a rough wall function with $z_0$ as roughness, the value of $\frac{y^+}{u^+}$ reads:
\begin{equation}\label{eq:bndcnd:ypup_rough}
\dfrac{y^+}{u^+} = \dfrac{y^+}{ \dfrac{1}{\kappa} \ln \left( \dfrac{y+z_0}{z_0}\right) + C_{log}}.
\end{equation}



%-------------------------------------------------------------------------------
\subsubsection{Wall shear stress with wall functions}
The previous sections can be summarized as follows:

\begin{equation}
\tens{\tau} \cdot \vect{n} = \dfrac{\mu}{\overline{\centip \centf} } \dfrac{y^+}{u^+} \left(\tens{1} - \vect{n} \otimes \vect{n} \right) \cdot \left(\vect{u}_\centip - \vect{u}_{wall} \right)
\end{equation}
%
where the rescaling factor $\frac{y^+}{u^+}$ depends on the wall function is given in \eqref{eq:bndcnd:ypup_1scale} or \eqref{eq:bndcnd:ypup_2scale}. Note that the wall shear stress is therefore parallel to the wall and fully implicit if the $\theta$-scheme on the velocity is implicit.

The internal pair of boundary condition coefficient for the diffusive term for the velocity are
\begin{equation}\label{eq:bndcnd:boundary_condition_wall_function}
\left\lbrace
\begin{array}{r c l}
\vect{A}^f_\ib & = & -h_{fluid}  \left(\tens{1} - \vect{n} \otimes \vect{n} \right) \cdot  \vect{u}_{wall}
 - h_{int} \left( \vect{n} \cdot \vect{u}_{wall}  \right) \vect{n}
, \\
\tens{B}^f_\ib & = & h_{fluid}  \left(\tens{1} - \vect{n} \otimes \vect{n} \right)
+ h_{int} \vect{n} \otimes \vect{n}
,
\end{array}
\right.
\end{equation}
where
\begin{equation}
h_{fluid} = \dfrac{\mu}{\overline{\centip \centf} } \dfrac{y^+}{u^+}.
\end{equation}

%-------------------------------------------------------------------------------
\subsubsection{Velocity gradient boundary condition with wall functions}
The gradient of the fluid velocity is of great importance for turbulence models
as it appears in the production term.
Moreover the profile of dissipation is even less linear than $u^+$ (it is
hyperbolic) and would also require a wall function. A common correction in open literature is
to directly modify these source term in the near wall cell volume integral. However this deeply
modifies the code structure and the linear system solver. In order to implement this while still
remaining within the standard boundary condition structure, ther current code version follows
the suggestion of \cite{BouckerMattei:2000} which consists in prescribing a slip velocity at the wall so
that the velocity gradient evaluated inside the cell by the
standard Finite Volumes linear profile is more consistent with wall
function.

A dedicated boundary condition for the velocity gradient is therefore derived
in presence of wall function to be consistent with the logarithmic law giving
the ideal tangential velocity profile.

On the first hand, the theoretical gradient is:

\begin{equation}\label{eq:bndcnd:vel_grad_wallf_theo}
\left( \dfrac{\partial u^r}{\partial y}\right)_{I'}^\textrm{theo} =\dfrac{u^\star}{\kappa \, y}
\end{equation}

On the other hand,
the normal finite volumes gradient of the tangential velocity reduces in the case
of a regular orthogonal mesh  (see notations in \figurename~
\ref{fig:bndcnd:boundary_cell_ortho_mesh}):
\begin{equation}\label{eq:bndcnd:vel_grad_wallf_calc}
\left( \dfrac{\partial u^r}{\partial y}\right)_{\centip}^\textrm{calc} = \dfrac{u_{G}^r - u_{F}^r}{2y} =
 \dfrac{\dfrac{u_{\centi}^r+u_{\centj}^r}{2} -  u_{\centf}^r}{2y}
\end{equation}

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.7\textwidth]{bordortho}
\caption{\label{fig:bndcnd:boundary_cell_ortho_mesh}Boundary cell - Orthogonal mesh.}
\end{figure}

We then assume that $u_{\centj}^r$ can be obtained from $u_{\centi}^r$
and from the gradient of $u^r$ calculated in $G$
from the logarithmic law
$u_{\centj}^r = u_{\centi}^r +2y \displaystyle\dfrac{u^\star}{\kappa\, 2y}$ and we thus
obtain equalizing \eqref{eq:bndcnd:vel_grad_wallf_theo} and \eqref{eq:bndcnd:vel_grad_wallf_calc}
 (the formula is extended without any precaution to non-orthogonal meshes):
%
\begin{equation}
u_{\centf}^r = u_{\centip}^r -  \dfrac{3 u^\star}{2\kappa }
\end{equation}
The normal derivative at the wall is prescribed to zero.
If the value obtained for $y^+$ at the wall is lower than $y^+_{lim}$,
a no-slip condition is prescribed. Finally, one can also
make explicit the velocity of the wall in the final expression:
%
The internal pair of boundary condition coefficient for the diffusive term for the velocity are
\begin{equation}
\left\lbrace
\begin{array}{r c l}
\vect{A}^g_\ib & = & \left( 1-cofimp \right) \left(\tens{1} - \vect{n} \otimes \vect{n} \right) \cdot  \vect{u}_{wall} + \left(\vect{u}_{wall} \cdot \vect{n} \right) \vect{n}, \\
\tens{B}^g_\ib & = & cofimp \left(\tens{1} - \vect{n} \otimes \vect{n} \right) ,
\end{array}
\right.
\end{equation}
where $cofimp = 1 - \dfrac{3}{2 \kappa u^+}$ if $y^+ > y^+_{lim}$, and $cofimp = 0$ otherwise.

%-------------------------------------------------------------------------------
\subsection{Scalar boundary condition for smooth walls}

In this section, the wall functions applied to transported scalars $\varia$ (such as the temperature $T$ for instance)
are described.
The reference "Convection Heat Transfer",
Vedat S. Arpaci and Poul S. Larsen, Prentice-Hall, Inc was used.
The approach is really similar to what is performed on the fluid velocity:
\begin{itemize}
\item if the wall value $\varia_{wall}$ is imposed, the Dirichlet boundary condition is
substituted by a Robin-type one where the diffusive flux is function of the cell centre scalar value,
\item if a Neumann value is imposed (that is to say that the user has prescribed the flux of $\varia$), then the wall
function is used to evaluate the corresponding wall value $\varia_{wall}$ which is displayed by the code.
\end{itemize}

 Let us recall that $q_\ib \left( K_\fib , \, \varia \right) $ is the wall diffusive flux for the scalar $\varia$,
 and the simplified scalar balance reads:
\begin{equation}\label{eq:bndcnd:wall_diffuive_flux_balance}
q_\ib \left( K_\fib , \, \varia \right)  =
-\left(K_\fib + C\,\frac{\mu_T}{\sigma_T}\right)
                  \dfrac{\partial \varia}{\partial y}
\end{equation}

\begin{remark}
the flux is positive
if it enters the fluid domain, as shown by the orientation of the
$y$ axis
\end{remark}

The following considerations are presented using the general notations.
In particular, the Prandtl-Schmidt number writes
$\sigma=\displaystyle\dfrac{\mu \,C}{K}$.
When the considered scalar $\varia$ is the temperature $T$,
we have
\begin{itemize}
\item $C=C_p$ (specific heat capacity),
\item $K=\lambda$ (molecular conductivity),
\item hence, $\sigma = \displaystyle\frac{\mu \,C_p}{\lambda} = Pr$
       (Prandtl number),
\item $\sigma_T = Pr_T$ (turbulent Prandtl number),
\item $q_\ib = \dfrac{D_\ib \left( \lambda_\fib , \, T \right)}{\norm{\vect{S}}_\fib} =\left(\lambda+\displaystyle\frac{C_p \mu_T}{\sigma_T}\right)
        \displaystyle \dfrac{\partial T}{\partial y}$ (flux in $W.m^{-2}$).
\end{itemize}


In order to make \eqref{eq:bndcnd:wall_diffuive_flux_balance} dimensionless,
we introduce $\varia^\star$ defined using
the flux at the boundary  $q_\ib$:
\begin{equation}\label{eq:bndcnd:def_varia_star}
\varia^\star \equiv -\displaystyle \dfrac{q_\ib}{\rho\,C\,u_k}
\end{equation}
For the temperature, we thus have:
\begin{equation}\label{eq:bndcnd:def_T_star}
T^\star \equiv -\displaystyle\frac{q_\ib}{\rho\,C_p\,u_k}
\end{equation}

We then divide both sides of equation \eqref{eq:bndcnd:wall_diffuive_flux_balance}
by $q_\ib$ and after some algebrae it comes:
\begin{equation}\label{eq:bndcnd:scalar_flux_adim}
1 =  \left(\displaystyle\frac{1}{\sigma}+
              \displaystyle\frac{1}{\sigma_T}\dfrac{\nu_T}{\nu}\right)
                  \displaystyle\dfrac{\partial \varia^+}{\partial y^+}
\end{equation}
%
where
\begin{equation}
       \nu=\displaystyle \dfrac{\mu}{\rho}
\qquad \nu_T=\displaystyle\dfrac{\mu_T}{\rho}
\qquad y^+=\displaystyle\dfrac{y\,u_k}{\nu}
\qquad \varia^+=\displaystyle \dfrac{\varia-\varia_{wall}}{\varia^\star}
\end{equation}

%-------------------------------------------------------------------------------
\subsection{The three layers wall function of Arpaci and Larsen}

In order to compute the theoretical scalar profile, we integrate
equation \eqref{eq:bndcnd:scalar_flux_adim} to obtain
$\varia^+_{\centip}$.
The only difficulty then consists in prescribing a variation law of
$\mathcal{K}=\displaystyle\frac{1}{\sigma}+ \displaystyle\frac{1}{\sigma_T}\frac{\nu_T}{\nu}$.

In the fully developed turbulent region
(far enough from the wall, for $y^+\geqslant y^+_2$),
a mixing length hypothesis models the variations of
$\nu_T$:
\begin{equation}
\nu_T = \kappa \,y\,u_k
\end{equation}
Additionally, the molecular diffusion of $\varia$
(or the conduction when $\varia$ represents the temperature)
is negligible compared to its turbulent diffusion: therefore
we neglect
$\displaystyle\frac{1}{\sigma}$ compared to
$\displaystyle\frac{1}{\sigma_T}\frac{\nu_T}{\nu}$.
Finally we have:
\begin{equation}
\mathcal{K}= \displaystyle\frac{\kappa \,y^+}{\sigma_T}
\end{equation}

On the contrary, in the near-wall region (for $y^+ < y^+_1$)
the turbulent contribution becomes negligible
compared to the molecular contribution and we thus neglect
$\displaystyle\frac{1}{\sigma_T}\frac{\nu_T}{\nu}$ compared to
$\displaystyle\frac{1}{\sigma}$.

It would be possible to restrict ourselves to these
two regions, but Arpaci and Larsen suggest the model
can be improved by introducing an intermediate
region ($y^+_1 \leqslant y^+ < y^+_2$)
in which the following hypothesis is made:
\begin{equation}
\frac{\nu_T}{\nu} = a_1 \left( y^+ \right)^3
\end{equation}
where $a_1$ is a constant whose value is obtained from
experimental correlations:
\begin{equation}
a_1 =\displaystyle\frac{\sigma_T}{1000}
\end{equation}

Thus the following model is used for $\mathcal{K}$
(see a sketch
in \figurename~ \ref{fig:bndcnd:3layer_temperature_wall_function}):
\begin{equation}
\mathcal{K}=\left\{
\begin{array}{ll}
\displaystyle\frac{1}{\sigma}
             &\text{if } y^+ < y^+_1\\
\displaystyle\frac{1}{\sigma}
+\displaystyle\frac{a_1 (y^+)^3}{\sigma_T}
             &\text{if } y^+_1 \leqslant y^+ < y^+_2\\
\displaystyle\frac{\kappa \,y^+}{\sigma_T}
             &\text{if } y^+_2\leqslant y^+\\
\end{array}
\right.
\end{equation}

\begin{figure}[!htp]
\centering
\includegraphics[width=0.8\textwidth]{clthermique}
\caption{$\frac{1}{\sigma} + \frac{1}{\sigma_T}\frac{\nu_T}{\nu} $ as a function of $y^+$ obtained
 for $\sigma=1$ and $\sigma_T=1$. \label{fig:bndcnd:3layer_temperature_wall_function}}
\end{figure}

The values of $y^+_1$ and $y^+_2$ are obtained by calculating
the intersection points of the variations laws used
for $\mathcal{K}$.

The existence of an intermediate region depends upon the
values of $\sigma$.
Let's first consider the case where $\sigma$ cannot be neglected
compared to $1$. In practise we consider  $\sigma > 0,1$
(this is the common case when scalar $\varia$ represents
the air or the water temperature in normal temperature
and pressure conditions). It is assumed that
$\displaystyle\dfrac{1}{\sigma}$ can be neglected compared to
$\displaystyle\dfrac{a_1 (y^+)^3}{\sigma_T}$ in the
intermediate region.
We thus obtain:
\begin{equation}
  y^+_1 =\left(\displaystyle\frac{1000}{\sigma}\right)^\frac{1}{3} \qquad\qquad
  y^+_2 = \sqrt{\displaystyle\frac{1000\kappa}{\sigma_T}}
\end{equation}
The dimensionless equation \eqref{eq:bndcnd:scalar_flux_adim}
is integrated under the same hypothesis and we obtain the law of $\dfrac{y^+}{\varia^+}$:
\begin{equation}\label{eq:bndcnd:wallfunctions_t_3layers}
\left\{
\begin{array}{ll}
\dfrac{y^+}{\varia^+} = \dfrac{1}{\sigma} & \text{if } y^+ < y^+_1 \\
\dfrac{y^+}{\varia^+} = \dfrac{y^+}{a_2 -\displaystyle\frac{\sigma_T}{2\,a_1\,(y^+)^2}}& \text{if } y_1^+ \leqslant y^+ < y_2^+ \\
\dfrac{y^+}{\varia^+} = \dfrac{y^+}{\dfrac{\sigma_T}{\kappa}\,\ln \left(y^+ \right)+a_3} & \text{if } y^+_2 \leqslant y^+
\end{array}
\right.
\end{equation}
where $a_2$ and $a_3$ are integration constants,
which have been chosen to obtain
a continuous  profile of $\varia^+$:
\begin{equation}
a_2=15\sigma^{\frac{2}{3}}\qquad\qquad
a_3=15\sigma^{\frac{2}{3}}-\displaystyle\frac{\sigma_T}{2\kappa}
\left(1+
\ln \left(\displaystyle\frac{1000\kappa}{\sigma_T}\right)\right)
\end{equation}

Let's now study the case where  $\sigma$ is much smaller than $1$.
In practise it is assumed that $\sigma \leqslant 0.1$ (this is for
instance the case for liquid metals whose thermal conductivity is very
large, and who have Prandtl number of values of the order of $0.01$).
The intermediate region then disappears and the coordinate of the
interface between the law used in the near-wall region and the one
used away from the wall is given by:
\begin{equation}
y^+_0= \displaystyle\frac{\sigma_T}{\kappa\sigma}
\end{equation}

The dimensionless equation \eqref{eq:bndcnd:scalar_flux_adim}
is then integrated under the same hypothesis, and the law of
 $\varia^+$ is obtained:
\begin{equation}\label{eq:bndcnd:wallfunctions_t_2layers}
\left\{
\begin{array}{ll}
\dfrac{y^+}{\varia^+} = \dfrac{1}{\sigma} & \text{if } y^+ \leqslant y^+_0 \\
\dfrac{y^+}{\varia^+} = \dfrac{y^+}{\displaystyle\frac{\sigma_t}{\kappa}\,
        \ln\left(\displaystyle\frac{y^+}{y^+_0}\right)+\sigma \,y^+_0}
                   & \text{if } y^+_0 < y^+
\end{array}
\right.
\end{equation}

%-------------------------------------------------------------------------------
\subsubsection{Wall diffusive flux with wall functions}
The previous sections can be summarized as follows:

\begin{equation}
q_\ib = -\dfrac{C \mu}{\overline{\centip \centf} } \dfrac{y^+}{\varia^+} \left(\varia_\centip - \varia_{wall} \right)
\end{equation}
%
where the rescaling factor $\frac{y^+}{\varia^+}$ depends on the wall function is given in \eqref{eq:bndcnd:wallfunctions_t_3layers} if $\sigma >0.1$ or \eqref{eq:bndcnd:wallfunctions_t_2layers} if $\sigma \leq 0.1$.

The internal pair of boundary condition coefficient for the diffusive term for the variable $\varia $ are
\begin{equation}
\left\lbrace
\begin{array}{r c l}
A^f_\ib & = & -h_{fluid} \varia_{wall}, \\
B^f_\ib & = & h_{fluid},
\end{array}
\right.
\end{equation}
where $h_{fluid} = \dfrac{C \mu}{\overline{\centip \centf} } \dfrac{y^+}{\varia^+} $.

%-------------------------------------------------------------------------------
\subsection{Outlet boundary condition on the pressure}\label{sec:bndcnd:outlet_pressure}

In this section the boundary condition on the pressure at the outlet is detailed. Some
assumptions are made to derive this boundary condition which consists of a Dirichlet
(combined with a homogeneous Neumann on the velocity) based on the pressure field
at the previous time step.
On basic configurations such as a channel or a pipe where the outlet is orthogonal
to the flow, the shape of the pressure profile in a surface parallel to the outlet is
approximately the shape of the pressure profile at the outlet. This hypothesis
is valid for established flows, far from any perturbation. In this configuration
one can write:

\begin{equation*}
\displaystyle\frac{\partial^2 P}{\partial n \partial \tau } = 0,
\end{equation*}
where $\vect{n}$ is the outward normal vector and  $\vect{\tau}$ is any
vector in the boundary face.

Then remark that the value at the boundary face $\fib$ is linked
to the pressure value in $\centip$ by the relationship:
\begin{equation*}
P_\fib = P_{\centip} +  \grad_\celli P \cdot \vect{\centip\centf}.
\end{equation*}
If moreover we assume that the pressure gradient in the normal direction
is uniform, and that all the $\centip$ related to the outlet faces are on
a single plane parallel to the outlet, then the value of $ \grad_\celli P \cdot \vect{\centip\centf} $
is constant for all the faces and denoted $R$.

Furthermore, the pressure is defined up to a constant, so the code
chooses to fix the pressure at $P_0$ to a given outlet boundary face
$\fib^{imp}$. Therefore the pressure field is shifted  by the constant
$R_0=P_0- P_\fib^{imp} = P_0 - \left( P_{\centip}^{imp}+R \right)$.

All that together gives
\begin{equation}
\begin{array}{rcl}
P_\fib
   &=& P_{\centip}+R+R_0,\\
   &=& P_{\centip}+R+P_0- \left( P_{\centip}^{imp}+R \right),\\
   &=&P_{\centip} +\underbrace{P_0-P_{\centip}^{imp}}_{\text{constant denoted by $\widetilde{R}$}} ,\\
   &=&P_{\centip} + \widetilde{R}.
\end{array}
\end{equation}

To conclude, the outlet boundary condition on the pressure is a Dirichlet based on
the value in $\centip$ (at the previous time step) and shifted to impose
the value $P_0$ at a given face $\fib$.

%-------------------------------------------------------------------------------
\subsection{Free inlet/outlet boundary condition on the pressure increment}\label{sec:bndcnd:inlet_bernoulli_pressure}

In this section the boundary condition on the pressure increment at the free inlet is detailed.
For the other variables (even the pressure itself), the treatment is similar to \S \ref{sec:bndcnd:outlet_pressure}.

Let us start with the Bernoulli relationship between a boundary face $\face$ linked to a point on the same
stream line far away from the computational domain:
\begin{equation}\label{eq:bndcnd:benoulli_inlet}
P_\face - \rho_\face \vect{g} \cdot  \left( \vect{x}_\face - \vect{x}_0 \right) + \dfrac{1+K}{2}\rho_\face \vect{u}_\face \cdot \vect{u}_\face  = P_\infty - \rho_\infty \vect{g} \cdot \left(  \vect{x}_\infty
- \vect{x}_0 \right)
+ \dfrac{1}{2}\rho_\infty \vect{u}_\infty \cdot \vect{u}_\infty,
\end{equation}
where $K$ is a possible head loss of the fluid between the infinity and the boundary face entrance
(which the user may play with to model the non-computed domain).

Assuming that the stream line is horizontal, that the density is constant over the stream line and that the fluid velocity is zero
at the infinity, the Equation \eqref{eq:bndcnd:benoulli_inlet} becomes:
\begin{equation}\label{eq:bndcnd:benoulli_inlet_simple}
P_\face^\star =- \dfrac{1+K}{2} \rho_\face \vect{u}_\face \cdot \vect{u}_\face ,
\end{equation}
where we recall that $P^\star$ is the dynamic pressure.

For stability reasons, the Equation \eqref{eq:bndcnd:benoulli_inlet_simple} is implicit in time as follows:
\begin{equation}\label{eq:bndcnd:benoulli_inlet_time}
P_\face^{\star, n+1} =- \dfrac{1+K}{2} \rho_\face \vect{u}^{n}_\face \cdot \vect{u}^{n+1}_\face.
\end{equation}
This implicitation is crucial to make the calculations stable.

The prediction-correction velocity-pressure coupling algorithm requires boundary conditions on the pressure increment
(computed in the correction step), and therefore
relation \eqref{eq:bndcnd:benoulli_inlet_time} is derived to obtain boundary conditions on the pressure increment. Recall that the
correction step reads:
 \begin{equation}\label{eq:bndcnd:correction_step}
\dfrac{\left( \rho \vect{u} \right)^{n+1} - \left( \rho^n \widetilde{\vect{u}} \right)}{\Delta t} = - \grad \delta P,
\end{equation}
the pressure increment $ \delta P = P^{\star, n+1} - P^{\star, n}$ at the face is consequently given by:
 \begin{equation}\label{eq:bndcnd:pressure_increment}
\delta P_\face = -P^{\star, n}_\face - \dfrac{1+K}{2} \rho_\face \vect{u}^{n}_\face \cdot \left( \widetilde{\vect{u}} - \dfrac{\Delta t}{\rho} \grad \delta P \right)_\face.
\end{equation}


Neglecting the tangencial velocity component at the entering faces ($\left(\rho \vect{u} \right)^n_\face \simeq \left(\rho \vect{u} \right)^n_\face \cdot \vect{n} \, \vect{n} = \dfrac{ \dot{m}_\face }{\norm{S}_\face} \vect{n} $) and approximating $P^{\star, n}_\face$ with its cell-center value, Equation \ref{eq:bndcnd:pressure_increment} reads:
%
 \begin{equation}\label{eq:bndcnd:pressure_increment_simple}
\delta P_\face = -P^{\star, n}_\celli - \dfrac{1+K}{2} \dfrac{ \dot{m}_\face }{\norm{S}_\face}
\left( \widetilde{\vect{u}} \cdot \vect{n}
- \dfrac{\Delta t}{\rho} \der{\delta P }{n} \right)_\face.
\end{equation}
Noting that $  \left. \der{\delta P }{n} \right|_{\face} = \dfrac{\delta P_\face -  \delta P_\celli }{\overline{\centip \centf}}$,  \eqref{eq:bndcnd:pressure_increment_simple} is rewritten as:
\begin{equation}\label{eq:bndcnd:pressure_increment_coefs_bernoulli}
\delta P_\face = A^g_{\delta P} + B^g_{\delta P} \delta P_\celli,
\end{equation}
with:
\begin{equation}\label{eq:bndcnd:pressure_increment_coefs_bernoulli_A_B}
\begin{array}{r c l}
 A^g_{\delta P}  & = & \dfrac{1}{1 + CFL} \left( - P^{\star, n}_\celli - \dfrac{1+K}{2} \dfrac{ \dot{m}_\face }{\norm{S}_\face}
\widetilde{\vect{u}}_\face \cdot \vect{n}
  \right), \\
 B^g_{\delta P} & = & \dfrac{CFL}{1 + CFL} ,
\end{array}
\end{equation}
where the $CFL$ is defined as:
\begin{equation}\label{eq:bndcnd:pressure_increment_coefs_bernoulli_A_B}
CFL = - \dfrac{1+K}{2} \dfrac{\Delta t \dot{m}_\face }{ \rho_\face \norm{S}_\face}.
\end{equation}
%
\begin{remark}
$CFL$ is a positive quantity for ingoing flows, which ensures stability.
\end{remark}
\begin{remark}
The formulation \eqref{eq:bndcnd:pressure_increment_coefs_bernoulli_A_B} is nothing else but
a convective outlet on the pressure increment (see \S \ref{sec:bndcnd:convective_outlet}).
\end{remark}

%-------------------------------------------------------------------------------
\subsubsection{Checking step}

 Before computing the pairs of boundary condition coefficients, a step of checking is performed. Basically,
the code checks that the user has given a boundary condition to all boundary faces, and that
the setting between all the variables is compatible (see \tablename~\ref{tab:ICODCLadm_condli}).



\begin{table}
%{\tiny
\begin{center}
\begin{tabular}{||c|c||p{0,6cm}|p{0,6cm}|p{0,6cm}|p{0,6cm}|p{0,6cm}|p{0,6cm}|p{0,6cm}||}
\hline
\multicolumn{2}{||c||}{Variable}
        &\multicolumn{7}{c||}{Admissible Values}\\
\hline
Velocity                                & $\vect{u}$                                                                         &  1& 2& 3& 4& 5& 6& 9 \\
Pressure                                & $P$                                                                                   &  1& 2& 3&  &  & & \\
Scalar turbulent variables     & $k$, $\varepsilon$, $\varphi$, $\bar{f}$, $\omega$        &  1& 2& 3&  & 5& 6& \\
Reynolds stresses                & $R_{ij}$                                                                              &  1& 2& 3& 4& 5& 6& \\
$\varia$ (except variances)  &  $\varia$                                                                           &  1& 2& 3&  & 5& 6& \\
Variance of a variable $\varia$ &                                                                                      &  1& 2& 3&  &  & & \\
\hline
\end{tabular}
\end{center}
%}
\caption{Admissible values of boundary conditions for all variables.}\label{tab:ICODCLadm_condli}
\end{table}