%~Mouliné par MaN_auto v.0.27.3 2023-07-06 10:28:35
\documentclass[CRMECA, Unicode, XML, screen,published]{cedram}




\usepackage{cleveref}
\usepackage{subcaption}
\usepackage{siunitx}
\usepackage{mathtools}
\usepackage{bm}





\newcommand{\CS}{{code\_saturne}}

\newcommand\vect[1]{\mathbf{\ensuremath\boldsymbol{#1}}}
\newcommand{\tens}[1]{\vect{\vect{#1}}}
\newcommand{\divs}{\text{div}\,}
\newcommand{\divv}{\vect{\mathbf{div}}\,}
\newcommand{\divt}{\tens{{\mathbf{div}}}\,}
\newcommand{\grad}{\vect{\nabla}}
\newcommand{\gradt}{\tens{\nabla}}
\newcommand{\gradtt}{\vect{\tens{\nabla}}}

\newcommand{\transpose}[1]{#1^T}
\newcommand{\trace}{tr}
\newcommand{\ugrad}{\left(\vect{\rans{u} }\cdot \grad \right)}
\newcommand{\dd}{\mathrm{d}}
\newcommand{\der}[2]{\dfrac{\partial #1}{\partial #2}}
\newcommand{\dert}[1]{\partial_t #1}
\newcommand{\dern}[1]{\partial_n #1}
\newcommand{\rans}[1]{\overline{#1}}
\newcommand{\fluct}[1]{{{#1}^\prime}}

\newcommand{\mean}[1]{\widehat{{#1}}}
\newcommand{\face}{{f}}
\newcommand{\celli}{c}
\newcommand{\vertex}{v}

\newcommand{\bbR}{\mathbb{R}}
\newcommand{\bfW}{\mathbf{W}}
\newcommand{\clF}{\mathcal{F}}
\newcommand{\clN}{\mathcal{N}}

\crefname{prop}{Proposition}{Propositions}
\crefname{propc}{Proposition}{Propositions}
\crefname{theo}{Theorem}{Theorems}
\crefname{defi}{Definition}{Definitions}
\crefname{rema}{Remark}{Remarks}
\crefname{remac}{Remark}{Remarks}
\crefname{lemm}{Lemma}{Lemmas}
\crefname{coro}{Corollary}{Corollaries}
\crefname{nota}{Notation}{Notations}
\crefname{clai}{Claim}{Claims}
\crefname{exem}{Example}{Example}
\crefname{exam}{Example}{Example}
\crefname{conj}{Conjecture}{Conjectures}
\crefname{section}{Section}{Sections}
\crefname{subsection}{Section}{Sections}
\crefname{appendix}{Appendix}{Appendices}
\crefname{figure}{Figure}{Figures}
\crefname{table}{Table}{Tables}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcounter{stepcount}
\newenvironment{step}[1]{\refstepcounter{stepcount}\begin{proof}[\textbf{Step~\thestepcount}]{\bf{#1}}}{\let\qed\relax\end{proof}}


\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

%\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\RubricEF{Spontaneous Articles}{Articles spontanés}


\title{Stable schemes for second-moment turbulent models for incompressible flows}

\author{\firstname{Martin} \lastname{Ferrand}\CDRorcid{0000-0002-3995-1503}\IsCorresp}
\address{CEREA, EDF R\&D, \'Ecole des Ponts, Chatou, France}
\email[M. Ferrand]{martin.ferrand@edf.fr}

\author{\firstname{Jean-Marc} \lastname{Hérard}\CDRorcid{0000-0001-9099-6072}}
\address[1]{CEREA, EDF R\&D, \'Ecole des Ponts, Chatou, France}
\address{I2M, Marseille, France}
\email[J.-M. Hérard]{jean-marc.herard@edf.fr}
\author{\firstname{Thomas} \lastname{Norddine}}
\email[T. Norddine]{thomas.norddine@edf.fr}
\address[1]{CEREA, EDF R\&D, \'Ecole des Ponts, Chatou, France}

\author{\firstname{Simon} \lastname{Ruget}}
\address{École des Ponts, Champs-sur-Marne, France}
\email[S. Ruget]{ruget.simon@gmail.com}



\keywords{second-moment turbulent closure, incompressible turbulent flows, hyperbolic systems, Riemann problem, numerical scheme}

\begin{abstract}
A stable scheme is proposed in this paper in order to obtain approximate solutions of second-moment turbulent models for incompressible flows with or without thermal transport equation. The analysis of the convective terms, which includes the solution of the associated Riemann problem, enables to propose a standard projection scheme, and to get rid of spurious oscillations.
\end{abstract}

%\selectlanguage{english}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\dateposted{2023-09-07}
\begin{document}

\maketitle

\section{Introduction}

Reynolds-averaged Navier--Stokes models are commonly used in the industrial and environmental applications in order to predict incompressible or compressible flows (see e.g.~\cite{pope2000turbulent,manceau2015recent}) due to the high Reynolds numbers or Grashoff numbers encountered. These equations require closure on the Reynolds stress tensor (covariance of velocity fluctuations), but second-order closures, which consist in solving a transport equation on this symmetric tensor, are less used. Temperature or pollutant averaged transport equations also provide turbulent thermal/scalar fluxes to be modelled. All these second-order models contain first-order non-conservative terms (see e.g.~\cite{le1989shock,dal1995definition,castro2008many,dumbser2016new}). When focusing on turbulent compressible flows, a stable numerical strategy has been proposed, which is grounded on the hyperbolic structure of the full-convective system embedded in the model (see~\cite{bhju2000,nasr2014low,audebertthesis,audebert2008structural}). More recently, this methodology has been extended to shallow-water models~\cite{dumbser2016new,bhole2019fluctuation,gavrilyuk2020variational,nkonga2022exact}.

The main purpose of the present work consists in applying similar ideas, while considering incompressible flow models including second-moment closures, in order to get rid of non-physical oscillations which occur in many practical situations (see e.g.~\cite{richards2011appropriate}).

We first recall the set of equations and the realizability constraint on the Reynolds stress tensor, which is mandatory to guarantee the hyperbolicity of the system. Then, we propose a basic time-stepping and a special focus is given on the analysis of the full-convection part of the models. Finally, the strategy is demonstrated using a Rusanov-like solver on a wall bounded flow. This strategy is extended to anisothermal flows.



\subsection{Set of equations}
The governing set of equations involves mass and momentum conservations and Reynolds stress tensor $\tens{R} \coloneqq \rans{ \fluct{\vect{u}} \otimes \fluct{\vect{u}}}$ transport equation:
\begin{subequations}\label{eq:system_dyn}
\begin{align}
\divs\left(\vect{\rans{u}} \right) & = 0, \label{eq:mass}
\\
\dert{\vect{\rans{u}}} +
\ugrad \vect{\rans{u}} + \divv \left(\tens{R} \right) - \divv\left(\rans{\tens{\Sigma}} \right) & =\vect{0}, \label{eq:momentum}
\\
\dert{\tens{R}} +
\ugrad \tens{R} +
\tens{R}\cdot \gradt\vect{\rans{u}} +\transpose{\gradt\vect{\rans{u}}} \cdot \tens{R} + \divt \left(\vect{\tens{G}}_{\tens{R}} \right) - \tens{\Phi} & = \tens{0}, \label{eq:rij}
\end{align}
\end{subequations}
with, for Newtonian fluid and incompressible flow, $\rans{\tens{\Sigma}} \coloneqq - \tfrac{\rans{p}}{\rho_0}\tens{1} + \nu_0 (\gradt\rans{\vect{u}} + \transpose{\gradt\rans{\vect{u}}}) $ the mean stress tensor per mass unit, and $\rans{p}$ denotes the mean mechanical pressure.

As a covariance matrix, the Reynolds stress tensor is symmetric positive and half-definite. This refers to the so-called realizability constraint (see~\cite{lumley_1979, pope2000turbulent}).

Models for $\vect{\tens{G}}_{\tens{R}}\coloneqq \rans{\fluct{\vect{u}} \otimes \fluct{\vect{u}} \otimes \fluct{\vect{u}}}$ usually rely on a classical first gradient assumption $\vect{\tens{G}}_{\tens{R}} = -C_R\nu_t(\tens{R}, \, T) \gradtt \tens{R}$, where $C_{R}$ is a positive constant, $\nu_t$ is the turbulent viscosity, and $T$ is the time scale related to turbulent dissipation.

Moreover $\tens{\Phi} \coloneqq \rans{
\fluct{\vect{u}} \otimes \divv(\fluct{\tens{\Sigma}} ) + \divv(\fluct{\tens{\Sigma}}) \otimes \fluct{\vect{u}} }$ is classically split into two parts:
\begin{equation}
\tens{\Phi} = \tens{\Phi}^r \left(\gradt \vect{\rans{u}}, \ \tens{R} \right) + \tens{\Phi}^s \left(\tens{R}, \, T\right),
\end{equation}
where the so-called rapid contribution verifies: $\tens{\Phi}^r(\tens{0}, \tens{R}) = \tens{0}$ and $\tens{\Phi}^r(\gradt \vect{\rans{u}}, \ \tens{0}) = 0$. Moreover $\tens{\Phi}^r$ is linear with respect to $\gradt \vect{\rans{u}}$ (see e.g.~\cite{pope2000turbulent} for a review of the turbulence models).


\subsection{Total kinetic energy tensor}\label{sec:total_kinetic_energy}
The evolution equation on $\tens{R}$ can be combined with the kinetic energy of the mean velocity field $\rans{\vect{u}} \otimes \rans{\vect{u}}$ to give the total kinetic energy $\rans{\vect{u} \otimes \vect{u}} = \rans{\vect{u}} \otimes \rans{\vect{u}} + \tens{R} $ evolution:
\begin{equation}\label{eq:rij_tot}
\dert{\rans{\vect{u} \otimes \vect{u}}} +
\divt \left(\rans{\vect{u} \otimes \vect{u} \otimes \vect{u}} \right) - \rans{ \vect{u} \otimes \divt \left(\tens{\Sigma} \right) +
\divt \left(\tens{\Sigma} \right)
\otimes \vect{u} } = \tens{0}. 
\end{equation}
Equation~\eqref{eq:rij_tot} can also be deduced directly from the mean of the instantaneous kinetic energy tensor, and emphasise that the so-called non-conservative production term $ - (\tens{R}\cdot \gradt\vect{\rans{u}} +\transpose{\gradt\vect{\rans{u}}} \cdot \tens{R})$ acts as energy transfer between kinetic energy of the mean velocity field and turbulent kinetic energy (note that the transfer may not be in the same direction) but is transparent with respect to the total kinetic energy. This is analogous to what happens to compressible flows with the non conservative term $p \divs \vect{u}$ term which cancels out between internal energy and kinetic energy balance. Therefore, the total kinetic energy tensor will be used in the following to derived jump conditions when dealing with non regular solutions.

\subsection{Realizability constraint}
To examine the realizability of the Reynolds stress tensor, we use the following Lemma:

\begin{lemm}\label{lem:realisability}
Assume that the governing equation of a symmetric second-rank tensor $\tens{R} \in \bbR^{3} \times \bbR^{3}$ is:
\[
\dert{\tens{R}} +
\ugrad \tens{R} +
\tens{R}\cdot \tens{H} +\transpose{\tens{H}} \cdot \tens{R} = \tens{0}.
\]
Let $\delta_{R}$ be the determinant of the latter tensor $\tens{R}$. Then the governing equation of $\delta_{R}$ is:
\[
\dert{\delta_{R}} +
\ugrad \delta_{R} + 2 \trace (\tens{H}) \delta_{R} = 0.
\]
\end{lemm}

A proof can be found in~\cite[Appendix~1]{CRAS96}. $\tens{H}$ is deduced from the terms $\tens{R}\cdot \gradt\vect{\rans{u}} +\transpose{\gradt\vect{\rans{u}}} \cdot \tens{R} $, $\tens{\Phi}^r$ and $\tens{\Phi}^s$, and the boundedness of the trace of $\tens{H}$ must be examined to ensure positivity of $\delta_{R}$. In case of the basic Lumley model examined in Section~\ref{sec:evolution_step} where $\tens{\phi}^r=\tens{0}$ and for incompressible flows, we have $\trace (\tens{H}) = \divs \rans{\vect{u}} = 0$ which implies that $\dert{\delta_{R}}$ remains strictly positive for suitable initial and boundary conditions, thus all the eigenvalues of $\tens{R}$ remain strictly positive. This is also true for the same model in the comrpessible framework if $\divs \rans{\vect{u}}$ is bounded.



\section{A first order time scheme}

Starting with initial conditions $(\vect{\rans{u}}^n, \,\tens{R}^n)$ and using suitable boundary conditions, $(\vect{\rans{u}}^{n+1}, \,\tens{R}^{n+1})$ arises from the following first-order time scheme:
\begin{subequations}
\begin{align}
\divs\left(\vect{\rans{u}}^{n+1} \right) &= 0, \label{eq:mass_stokes}
\\
\dfrac{\vect{\rans{u}}^{n+1} - \vect{\rans{u}}^{n}}{\Delta t} + \grad \dfrac{\rans{p}^{n+1} }{\rho_0} - \divv\left(\nu_0 \gradt\rans{\vect{u}}^{n+1}\right) &= -\left(
\ugrad \vect{\rans{u}} + \divv \left(\tens{R} \right)\right)^n, \label{eq:momentum_stokes}
\\
\dfrac{\tens{R}^{n+1} - \tens{R}^{n}}{\Delta t} - \divt \left(C_R \nu_t \gradtt \tens{R}^{n+1} \right) - \tens{\Phi}^s\left(\tens{R}^{n+1}, T \right) &=\nonumber\\ 
-&\left(
\ugrad \tens{R} +
\tens{R}\cdot \gradt\vect{\rans{u}} +\transpose{\gradt\vect{\rans{u}}} \cdot \tens{R} - \tens{\Phi}^r\left(\gradt \vect{\rans{u}}, \tens{R} \right)\right)^{n}. \label{eq:rij_stokes}
\end{align}
\end{subequations}

A Stokes-like step has been applied in the scheme~\eqref{eq:mass_stokes} and~\eqref{eq:momentum_stokes} (see, among others, \cite{temam2001navier, chorin1968numerical, bell1989second}). The space discretization associated with right-hand-side of~\eqref{eq:momentum_stokes}, \eqref{eq:rij_stokes} is derived from the analysis of the evolution system~\eqref{eq:evol-u-R}, which involves all convective terms:
\begin{subequations}\label{eq:evol-u-R}
\begin{align}
\dert{\vect{\rans{u}}} +
\ugrad \vect{\rans{u}} + \divv \left(\tens{R} \right) & = \vect{0}, \label{eq:evol-u}
\\
\dert{\tens{R}} +
\ugrad \tens{R} +
\tens{R}\cdot \gradt\vect{\rans{u}} +\transpose{\gradt\vect{\rans{u}}} \cdot \tens{R} - \tens{\Phi}^r\left(\gradt \vect{\rans{u}}, \tens{R} \right) & = \tens{0}. \label{eq:evol-R}
\end{align}
\end{subequations}

\section{Analysis of the evolution step}\label{sec:evolution_step}
In the following, we restrict to the simple Lumley's model which does not contain so-called rapid terms $\tens{\Phi^r} = \tens{0}$. We also restrict from now on to the two-dimensional framework to ease the presentation. We first recall that the latter system~\eqref{eq:evol-u-R} is hyperbolic if the Reynolds stress tensor is realizable (see~\cite{herard1994basic}).


Eq.~\eqref{eq:evol-u-R} is invariant under frame rotation, so we can rewrite it in the $(\vect{n}, \, \vect{\tau})$ reference frame (where $\vect{n}, \vect{\tau}$ are some unit orthogonal vectors), defining $x_n \coloneqq \vect{x} \cdot \vect{n}$, $u_n \coloneqq \vect{\rans{u}} \cdot \vect{n}$, $u_{\tau} \coloneqq \vect{\rans{u}} \cdot \vect{\tau}$, $R_{nn} \coloneqq \vect{n}^T \tens{R} \vect{n} $, $R_{n\tau} \coloneqq \vect{n}^T \tens{R} \vect{\tau} $ and $R_{\tau\tau} \coloneqq \vect{\tau}^T \tens{R} \vect{\tau} $. We now define the evolution step in the $\vect{n}-$direction, by neglecting the transverse derivatives, which yields:
\begin{subequations}\label{eq:evol-u-R-ntau}
\begin{align}
\dert{u_n} + u_n \dern{u_n} + \dern{R_{nn}} &= 0, \\
\dert{u_\tau}  + u_n \dern{u_\tau} + \dern{R_{n\tau}} &= 0, \\
\dert{R_{nn}} + u_n \dern{R_{nn}} + 2 R_{nn} \dern{u_n} &= 0, \\
\dert{R_{\tau\tau}} + u_n \dern{R_{\tau\tau}} + 2 R_{n\tau} \dern{u_\tau} &= 0, \\
\dert{R_{n\tau}} + u_n \dern{R_{n\tau}} + R_{nn} \dern{u_\tau} + R_{n\tau} \dern{u_n} &= 0.
\end{align}
\end{subequations}
We set $\vect{W} \coloneqq (u_n,\, u_\tau, R_{nn}, \, R_{\tau \tau}, \, R_{n \tau})$.


\subsection{Hyperbolicity of the evolution system}
We assume from now on that initial and boundary conditions comply with the strictly realizability constraint, which will imply that $\tens{R}$ remains strictly realizable (see~\cite{Schumann77}). As system~\eqref{eq:evol-u-R-ntau} preserves the realizability, it implies that $R_{nn} > 0$ whatever the unit vector $\vect{n}$ is.

\begin{prop}[Hyperbolic evolution step]\label{prop:hyperbolic_evolution_step}
System~\eqref{eq:evol-u-R-ntau} is hyperbolic. It admits five real eigenvalues which are:
\begin{equation}\label{syst-rij-eigenvalues}
\begin{split}
\lambda_1 &= u_n - \sqrt{2 R_{nn}}, 
\lambda_2 = u_n - \sqrt{R_{nn}}, 
\lambda_3 = u_n,\\
\lambda_4 &= u_n + \sqrt{R_{nn}},\;\;\: 
\lambda_5 = u_n + \sqrt{2 R_{nn}},
\end{split}
\end{equation}
and the associated right eigenvectors span $\bbR^5$.
\end{prop}


We refer to~\cite{herard:hal-02007060,herard1994basic} for a proof. This can be extended to the three-dimensional framework.

\subsection{Wave structure}

Right eigenvectors are:
\begin{subequations}
\begin{align}
\vect{r}_{1} & = \left(1, \ \frac{R_{n\tau}}{R_{nn}}, \ - \sqrt{2 R_{nn}}, \ -\frac{\sqrt{2}R_{n\tau }^{2}}{(R_{nn})^{\frac{3}{2}}}, \ - \frac{\sqrt{2}R_{n\tau}}{\sqrt{R_{nn}}} \right) ^{T}, \\
\vect{r}_{2} & = \left(0, \ 1, \ 0, \ -\frac{2R_{n\tau }}{\sqrt{R_{nn}}}, \ - \sqrt{R_{nn}} \right) ^{T}, \\
\vect{r}_{3} & = \left(0, \ 0, \ 0, \ 1, \ 0 \right) ^{T}, \\
\vect{r}_{4} & = \left(0, \ 1, \ 0, \ \frac{2R_{n\tau }}{\sqrt{R_{nn}}}, \ \sqrt{R_{nn}} \right) ^{T}, \\
\vect{r}_{5} & = \left(1, \ \frac{R_{n\tau }}{R_{nn}}, \ \sqrt{2 R_{nn}}, \ \frac{\sqrt{2}R_{n\tau }^{2}}{(R_{nn})^{\frac{3}{2}}}, \ \frac{\sqrt{2}R_{n\tau}}{\sqrt{R_{nn}}} \right) ^{T}.
\end{align}
\end{subequations}


Fields associated with eigenvalues $\lambda_{2,3,4}$ are linearly degenerate (LD). The $1$-field and $5$-field are genuinely non linear (GNL).

\subsection{Riemann invariants}
The Riemann invariants $I^{i}_{R}$ associated with the $i^{\rm th}$ field ($f \in I^i_R$, $\vect{\nabla} f \cdot \vect{r}_{i} = 0$) are:
\begin{subequations}\label{riemanninvariants}
\begin{align}
I^{1}_{R} &= \left\{ u_n +\sqrt{2R_{nn}}, \ \frac{R_{n\tau }}{R_{nn}}, \ R_{\tau \tau } - \frac{R_{n\tau }^2}{R_{nn}}, \ u_\tau + R_{n\tau }\sqrt{\frac{2}{R_{nn}}}\right\},\label{eq:riemanninvariants:a}
\\
I^{2}_{R} &= \left\{ u_n, \ R_{nn}, \ u_\tau+\frac{R_{n\tau }}{\sqrt{R_{nn}}}, \ R_{\tau \tau }R_{nn} - R_{n\tau }^{2} \right\},\label{eq:riemanninvariants:b}
\\
I^{3}_{R} &= \Bigg\{ u_n, \ u_\tau, \ R_{nn}, \ R_{n\tau } \Bigg\},\label{eq:riemanninvariants:c}
\\
I^{4}_{R} &= \left\{ u_n, \ R_{nn}, \ u_\tau-\frac{R_{n\tau }}{\sqrt{R_{nn}}}, \ R_{\tau \tau }R_{nn} - R_{n\tau }^{2} \right\},\label{eq:riemanninvariants:d}
\\
I^{5}_{R} &= \left\{ u_n -\sqrt{2R_{nn}}, \ \frac{R_{n\tau }}{R_{nn}}, R_{\tau \tau } -\frac{R_{n\tau }^2}{R_{nn}}, \ u_\tau - R_{n\tau } \sqrt{\frac{2}{R_{nn}}} \right\}.\label{eq:riemanninvariants:e}
\end{align}
\end{subequations}

\subsection{Shock waves}
Due to non-conservative terms in~\cref{eq:evol-u-R-ntau}, shock relations cannot be derived in a classical way. Following~\cite{dal1995definition, le1989shock} \emph{approximate} shock relations can nevertheless be proposed. These relations are valid in the limit of weak shocks. Assuming a linear path from left to right with respect to the variable $\vect{W} = (u_n,\, u_\tau, R_{nn}, \, R_{\tau \tau}, \, R_{n \tau})$, these relations write:

\begin{subequations}\label{eqn:jump_lefloch}
\begin{align}
\sigma [u_n] &= \mean{u_n}[u_n] + [R_{nn}], \\
\sigma [u_\tau] &= \mean{u_n}[u_\tau] + [R_{n\tau}], \\
\sigma [R_{nn}] &= 2\mean{R_{nn}}[u_n] + \mean{u_n}[R_{nn}], \\
\sigma [R_{\tau\tau}] &= 2\mean{R_{n\tau}}[u_\tau] + \mean{u_n}[R_{\tau\tau}],\\
\sigma [R_{n\tau}] &= \mean{R_{n\tau}}[u_n] + \mean{R_{nn}}[u_\tau] + \mean{u_n}[R_{n\tau}],
\end{align}
\end{subequations}
where $[x] \coloneqq x_r - x_l$ denotes jump and $\mean{x} \coloneqq \frac{x_r + x_l}{2}$ is the arithmetic mean between left and right states of the discontinuity travelling at speed $\sigma$.

Approximate relations~\eqref{eqn:jump_lefloch} can also be obtained from transport equations written in pseudo-conservative variables $\tens{R}+\vect{\rans{u}} \otimes \vect{\rans{u}} $ instead of $ \tens{R}$ introduced in Section~\ref{sec:total_kinetic_energy}, as explained in~\cite[p.~18 and Appendix~4]{herard:hal-02007060}.



\subsection{An approximate solution of the one-dimensional Riemann problem}

\begin{theo}[Existence and uniqueness of the approximate solution]\label{exiunith}
For the approximate shock relations~\eqref{eqn:jump_lefloch}, there exists a unique self-similar \emph{realizable} solution ($\vect{W}(x_n, t) = \vect{w}(\frac{x_n}{t})$, with $\vect{W} = (u_n,\, u_\tau, R_{nn}, \, R_{\tau \tau}, \, R_{n \tau})$) to the Riemann problem associated with system~\eqref{eq:evol-u-R-ntau} supplemented with strictly realizable initial conditions $\vect{W}(x_n, 0) = \vect{W}^L $ for $x_n <0$ and $\vect{W}(x_n, 0) = \vect{W}^R $ for $x_n > 0$, if and only if the following condition holds:
\begin{equation}\label{conditionexiuni}
u_n^R - u_n^L < \left(\sqrt{2R_{nn}^L} + \sqrt{2R_{nn}^R}\right).
\end{equation}
\end{theo}

As in the compressible framework~\cite{audebert2008structural,audebertthesis,bhju2000,nasr2014low}, the Riemann problem structure of system~\eqref{eq:evol-u-R} allows to specify wall, symmetry and inlet/outlet boundary conditions.

\begin{proof}
The proof follows the main guidelines recalled in~\cite{smoller2012shock}. A detailed expression of intermediate states are fully given in the sequel for the anisothermal case in Section~\ref{sec:anisothermal_extension}.

Let us consider two initial realizable states $\vect{W}^R$ and $\vect{W}^L$.
\setcounter{stepcount}{-1}
\begin{step}{Preliminary calculation.}
To connect left and right states $\vect{W}^{l,r}$ across the $1^{st}$ and $5^{th}$ waves we set $z \coloneqq \frac{R_{nn}^r}{R_{nn}^l} >0$. To select the physical solution for shocks~\eqref{eqn:jump_lefloch}, we use $[u_n] = u_n^r -u_n^l \leq 0$ (see~\cite[Appendix~4]{herard:hal-02007060}). We can also write from Section~\ref{sec:total_kinetic_energy}:
\begin{equation*}
-\sigma \left[u_n^2 + R_{nn}\right] + \left[u_n^3 + 3 u_n R_{nn}\right] \leq 0,
\end{equation*}
in shocks, and therefore, using~\eqref{eqn:jump_lefloch}, we get: $ (\mean{u_n^2} + \mean{R_{nn}})[u_n] \leq 0 $, and thus once more $[u_n] \leq 0$.

For the 1-shock connexion ($z>1$):
\begin{equation}\label{1shockconnexion}
\begin{split}
R_{nn}^r &= zR_{nn}^l;\quad u_n^r = u_n^l + \frac{1-z}{\sqrt{z+1}}\sqrt{R_{nn}^l};\quad u_\tau^r = u_\tau^l + \frac{1-z}{\sqrt{z+1}}\frac{R_{n\tau}^l}{\sqrt{R_{nn}^l}};\\
 R_{n\tau}^r &= zR_{n\tau}^l;\quad R_{\tau\tau}^r = R_{\tau\tau}^l + (z-1)\frac{\left(R_{n\tau}^l\right)^2}{R_{nn}^l}.
\end{split}
\end{equation}
For the 5-shock connexion ($0<z<1$):
\begin{equation}\label{5shockconnexion}
\begin{split}
R_{nn}^r &= zR_{nn}^l;\quad u_n^r = u_n^l + \frac{z-1}{\sqrt{z+1}}\sqrt{R_{nn}^l};\quad u_{\tau}^r = u_\tau^l + \frac{z-1}{\sqrt{z+1}}\frac{R_{n\tau}^l}{\sqrt{R_{nn}^l}};\\
 R_{n\tau}^r &= zR_{n\tau}^l;\quad R_{\tau \tau}^r = R_{\tau \tau}^l + (z-1)\frac{\left(R_{n\tau}^l\right)^2}{R_{nn}^l}.
\end{split}
\end{equation}

To connect the left and right states through the 1-rarefaction wave (respectively 5-rarefaction wave) we use the Riemann invariants~\eqref{eq:riemanninvariants:a} (respectively~\eqref{eq:riemanninvariants:e}) still following~\cite{smoller2012shock}.
\end{step}

\begin{step}{Solution in terms of $u_n$ and $R_{nn}$ variables}
A glance at~\eqref{eq:riemanninvariants:b}, \eqref{eq:riemanninvariants:c}, \eqref{eq:riemanninvariants:d} shows that $u_n$ and $R_{nn}$ are constant across the three LD waves. First we can focus only on $u_n$ and $R_{nn}$ variables. We denote their intermediate values by $u_n^{\sharp}$ and $R_{nn}^{\sharp}$, and link these to initial states $\vect{W}^{R,L}$:
\begin{equation}\label{lefttoright}
\begin{cases}
u_n^L &= u_n^{\sharp} - h_1(z_1)\sqrt{2R_{nn}^{L}}, \text{\ with }z_1 \coloneqq \frac{R_{nn}^{\sharp}}{R_{nn}^{L}},\\
u_n^R &= u_n^{\sharp} + h_5(z_5)\sqrt{2R_{nn}^{R}}, \text{\ with }z_5 \coloneqq \frac{R_{nn}^R}{R_{nn}^{\sharp}},
\end{cases}
\end{equation}
$h_{1,5}(z)$ are defined using \cref{1shockconnexion,5shockconnexion,eq:riemanninvariants:a,eq:riemanninvariants:e}, following~\cite{smoller2012shock}. Note that $z_1 z_5 = \frac{R_{nn}^R}{R_{nn}^{L}}$. We set:
\begin{equation}\label{eq:psi}
\Psi(z_5) \coloneqq u_n^R -u_n^L - \sqrt{2R_{nn}^R}h_5(z_5) - \sqrt{2R_{nn}^L}h_1(z_1(z_5)) = 0.
\end{equation}
Solving \cref{eq:psi} gives $z_5$, hence $R_{nn}^{\sharp}$ and $u_n^{\sharp}$. The solution is unique since $\Psi$ is a strictly monotonic function. Its lower and upper bounds enable to conclude that $u_n^{\sharp}$ and $R_{nn}^{\sharp}$ exist if~\eqref{conditionexiuni} holds.
\end{step}

\begin{step}{Solution in terms of $u_\tau$, $R_{n\tau}$, and $R_{\tau\tau}$ variables}
For these variables, the connexion from left to right can be done using $z_1$ and $z_5$, and exploiting the Riemann-invariants across the $2^{\rm nd}$, $3^{\rm rd}$ and $4^{\rm th}$ LD waves. $R_{nn}$ and $R_{nn}R_{\tau\tau} - R_{n\tau}^2$ remain constant across the $2^{\rm nd}$ and $4^{\rm th}$ waves (see~\eqref{eq:riemanninvariants:b},\eqref{eq:riemanninvariants:d}), thus $\tens{R}$ is realizable. Moreover $u_\tau$ does not vary across the $3^{\rm rd}$ wave (see~\eqref{eq:riemanninvariants:c}).\qedhere
\end{step}\let\qed\relax
\end{proof}



\section{Numerical scheme for the evolution system}

The objective of the present work is to obtain a stable numerical scheme (unlike those usually used as presented in~\cite{richards2011appropriate}, which may use partial upwinding which respect to the material mean velocity). Therefore we restrict the presentation to low order schemes such as Rusanov scheme. Other approximate Godunov schemes and MUSL techniques might also be considered.

For the sake of brevity the scheme is presented in a one dimensional framework. The evolution step~\eqref{eq:evol-u-R} can be written (using the incompressibility constraint):
\begin{equation}
\der{\vect{W}}{t} + \der{\left(u_n\vect{W}\right)}{x_n} + \tens{C}(\vect{W})\der{\vect{W}}{x_n} = \vect{0},
\end{equation}
$\vect{W}_i^{n+1}$ will be computed using the following Rusanov-like~\cite{rusanov1961calculation} scheme:
\begin{equation}\label{scheme}
\Delta x_i \left(\vect{W}_i^{n+1} - \vect{W}_i^{n}\right) + \Delta t^n
\left(\vect{\clF}^{n}_{i+\frac{1}{2}} - \vect{\clF}^{n}_{i-\frac{1}{2}}\right) + \Delta t^n \vect{\clN}_{i}^n = \vect{0},
\end{equation}
where the time step $\Delta t^n$ satisfies the CFL condition $\max_{i}(|\lambda_{i+\frac{1}{2}}^n|)\cdot\frac{\Delta t^n}{\Delta x_i} \leq \frac{1}{2}$, $\Delta x_i$ denotes the size of cell $i$. The numerical flux $\vect{\clF}_{i+\frac{1}{2}}^n$ and the non-conservative contribution $\vect{\clN}_{i}^n$ are:
\begin{subequations}\label{eq:scheme_rusanov}
\begin{align}
\vect{\clF}_{i+\frac{1}{2}}^{n} &\coloneqq (u_{n})_{i+\frac{1}{2}}^n\frac{\vect{W}_i^n +\vect{W}_{i+1}^n}{2} - \frac{\lambda^{n}_{i+\frac{1}{2}}}{2}\left(\vect{W}_{i+1}^n - \vect{W}_i^n\right),\nonumber\\
&\mkern150mu\text{\ with }\lambda^{n}_{i+\frac{1}{2}} = \max_{k} \max_{l=i, i+1} \left(\left|\lambda_k\left(u_n\tens{I}+ \tens{C}\right)\left(\vect{W}_l^n\right) \right|\right),
\\
\vect{\clN}_{i}^n &\coloneqq \tens{C}\left(\vect{W}_i^n\right)\left(\frac{\vect{W}_{i+1}^n - \vect{W}_{i-1}^n}{2}\right),
\end{align}
\end{subequations}
where the discrete divergence free condition on $(u_{n})_{i+\frac{1}{2}}^n$ holds.

Unlike scheme~\eqref{scheme} and~\eqref{eq:scheme_rusanov}, standard schemes only consider material convection upwinding for stabilization: $ \lambda^{n}_{i+\frac{1}{2}} = |(u_{n})_{i+\frac{1}{2}}^n|$.

\section{Wall bounded flow test case}

The following test case describes a flow in the vicinity of a wall, while prescribing shear stress. The initial states $\vect{W}^{L,R}$ and the intermediate states ($I$, $II$, $III$, $IV$) arising in the $1$-dimensional Riemann problem associated with~\eqref{eq:evol-u-R-ntau} are given in \cref{analyticssolCase2}. Computations have been performed using \CS{} finite volume platform~\cite{code_saturne}. The CFL parameter is set to $0.5$ and the meshes contain from $250$ to $32000$ cells. The L$1$-error and the behaviour of $\vect{W}$ are shown in \cref{Figure:simulation}, using material-upwind scheme and scheme~\eqref{scheme} and~\eqref{eq:scheme_rusanov} accounting for all convective effects.

\begin{table}[h]
\caption{\label{analyticssolCase2}
Analytical solution for the tangential flow in the vicinity of a wall.}
\vspace*{-10pt}
\begin{center}
\hspace*{-20pt}{
\begin{tabular}{ | c || c | c | c | c | c | c |}
\hline zone & $L$ & $I$ & $II$ & $III$ & $IV$ & $R$\\ \hline

$u_n\left(\si{\meter\per\second}\right)$ & $u_n(=0)$ & $u_n$ & $u_n$ & $u_n$ & $u_n$ & $u_n(=0)$\\ \hline

$u_\tau\left(\si{\meter\per\second}\right)$ & $u_\tau^L(=1)$ & $u_\tau^L$ & $\frac{u_\tau^L + u_\tau^R}{2} + \frac{R_{n\tau}^L - R_{n\tau}^R}{2\sqrt{R_{nn}}}$ & $\frac{u_\tau^L + u_\tau^R}{2} + \frac{R_{n\tau}^L - R_{n\tau}^R}{2\sqrt{R_{nn}}}$ & $u_\tau^R$ & $u_\tau^R(=1)$ \\ \hline

$R_{nn}\left(\si{\square\meter\per\square\second}\right)$ & $R_{nn}(=0.5)$ & $R_{nn}$ & $R_{nn}$ & $R_{nn}$ & $R_{nn}$ & $R_{nn}(=0.5)$ \\ \hline

$R_{n\tau}\left(\si{\square\meter\per\square\second}\right)$ & $R_{n\tau}^L(=0.4)$ & $R_{n\tau}^L$ & $\frac{\left(u_\tau^L -u_\tau^R\right) \sqrt{R_{nn}}}{2} + \frac{R_{n\tau}^L + R_{n\tau}^R}{2}$ & $\frac{\left(u_\tau^L -u_\tau^R\right) \sqrt{R_{nn}}}{2} + \frac{R_{n\tau}^L + R_{n\tau}^R}{2}$ & $R_{n\tau}^R$ & $R_{n\tau}^R(=0.5)$ \\ \hline

$R_{\tau\tau}\left(\si{\square\meter\per\square\second}\right)$ & $R_{\tau\tau}^L(=0.3)$ & $R_{\tau\tau}^L$ & $R_{\tau\tau}^L + \frac{\left(R_{n\tau}^{II}\right)^2 - \left(R_{n\tau}^L\right)^2}{R_{nn}}$ & $R_{\tau\tau}^R + \frac{\left(R_{n\tau}^{III}\right)^2 - \left(R_{n\tau}^R\right)^2}{R_{nn}}$ & $R_{\tau\tau}^R$ & $R_{\tau\tau}^R(=0.4)$ \\ \hline
\end{tabular}}
\end{center}
\vspace*{-15pt}
\end{table}

\section{Extension to anisothermal flows with second-order turbulence models}\label{sec:anisothermal_extension}

System~\eqref{eq:system_dyn} is now supplemented with mean temperature $\rans{\theta}$ conservation \cref{eq:evol-u-R-T:a}, turbulent heat fluxes $\rans{ \theta'\vect{u}'}$ and temperature variance governing \cref{eq:evol-u-R-T:b,eq:evol-u-R-T:c}:
\begin{subequations}\label{eq:evol-u-R-T}
\begin{align}
\dert{\rans{\theta}} +
\ugrad
\rans{\theta} + \divs\left(\rans{ \theta'\vect{u}'}\right) - \divs \left(
\vect{\Sigma}_{\rans{\theta}}
\right)&= 0,\label{eq:evol-u-R-T:a}
\\
\der{\rans{\vect{u}' \theta'}}{t} +
\ugrad {\rans{ \theta'\vect{u}'}} + \rans{ \theta'\vect{u}'}\cdot\gradt{\rans{\vect{u}}} + \tens{R}\cdot\grad \rans{\theta} + \divv \tens{G}_{\rans{ \theta'\vect{u}'}} -\vect{\Phi}^s_{\rans{ \theta'\vect{u'}}}\left(\tens{R}, \rans{ \theta'\vect{u}'}, \rans{\theta'^2}, T, T_{\rans{\theta}}\right) &= \vect{0},\label{eq:evol-u-R-T:b}
\\
\der{\rans{\theta'^2}}{t} + \ugrad
\rans{\theta'^2} + 2\left(\rans{ \theta'\vect{u}'}\right) \cdot \vect{\nabla} \rans{\theta} +\divs{\vect{G}_{\rans{\theta'}^2}} -\Phi^s_{\rans{\theta'}^2}\left(\tens{R}, \rans{ \theta'\vect{u}' }, \rans{\theta'^2}, T, T_{\rans{\theta}}\right)&= 0\label{eq:evol-u-R-T:c}.
\end{align}
\end{subequations}
Focus is still given on first-order gradient closure laws for $\vect{G}_{\rans{\theta'}^2}$ and $\tens{G}_{\rans{\vect{u}'\theta'}}$. We consider the one dimensional Riemann problem involving \cref{eq:evol-u-R} with $\tens{\Phi}^r=\tens{0}$, supplemented with the first-order terms of \cref{eq:evol-u-R-T:a,eq:evol-u-R-T:b,eq:evol-u-R-T:c} (thus getting rid of $\vect{G}_{\rans{\theta'}^2}$, $\tens{G}_{\rans{\vect{u}'\theta'}}$, $\vect{\Phi}^s_{\rans{\vect{u'}\theta'}}$, $\Phi^s_{\rans{\theta'}^2}$ and $\vect{\Sigma}_{\rans{\theta}}$).

Using the invariance under frame rotation, the counterpart of~\eqref{eq:evol-u-R-ntau} can be derived by setting to zero the transverse derivatives, which yields:
\begin{subequations}\label{eq:evol-u-R-ntau_temp}
\begin{align}
\dert{u_n} + u_n \dern{u_n} + \dern{R_{nn}} &= 0, \\
\dert{u_\tau} + u_n \dern{u_\tau} + \dern{R_{n\tau}} &= 0, \\
\dert{R_{nn}} + u_n \dern{R_{nn}} + 2 R_{nn} \dern{u_n} &= 0, \\
\dert{R_{\tau\tau}} + u_n \dern{R_{\tau\tau}} + 2 R_{n\tau} \dern{u_\tau} &= 0, \\
\dert{R_{n\tau}} + u_n \dern{R_{n\tau}} + R_{nn} \dern{u_\tau} + R_{n\tau} \dern{u_n} &= 0, \\
\dert{\rans{u_n'\theta'}} + u_n \dern{\rans{u_n'\theta'}} + \rans{u_n'\theta'} \dern{u_n} + R_{nn} \dern{\rans{\theta}}&= 0, \\
\dert{\rans{u_\tau'\theta'}} + u_n \dern{\rans{u_\tau'\theta'}} + \rans{u_n'\theta'} \dern{u_\tau} + R_{n\tau} \dern{\rans{\theta}}&= 0, \\
\dert{\rans{\theta'^2}} + u_n \dern{\rans{\theta'^2}} + 2 \rans{u_n'\theta'} \dern{\rans{\theta}} &= 0, \\
\dert{\rans{\theta}} + u_n \dern{\rans{\theta}} + \dern{ \rans{u_n'\theta'} } &= 0.
\end{align}
\end{subequations}
We set $\vect{W} \coloneqq (u_n,\, u_\tau, R_{nn}, \, R_{\tau \tau}, \, R_{n \tau}, \, \rans{u_n'\theta'}, \,
\rans{u_\tau'\theta'}, \, \rans{\theta'^2}, \, \rans{\theta})$.


\begin{prop}[Hyperbolic evolution step]\label{prop:hyperbolic_eiolution_step_temp}
System~\eqref{eq:evol-u-R-ntau_temp} is hyperbolic. It admits nine real eigenvalues which are:
\begin{equation}\label{syst-rij_temp-eigenvalues}
\begin{array}{
r @{\,}c@{\,}l @{\hspace{0.1\textwidth}} r @{\,}c@{\,}l @{\hspace{0.1\textwidth}} r @{\,}c@{\,}l }
\lambda_1 &=& u_n - \sqrt{2 R_{nn}}, &
\lambda_{2,3} &=& u_n - \sqrt{R_{nn}}, &
\lambda_{4,5,6} &=& u_n,\\
\lambda_{7,8} &=& u_n + \sqrt{R_{nn}}, &
\lambda_{9} &=& u_n + \sqrt{2 R_{nn}},
\end{array}
\end{equation}
and the associated right eigenvectors span $\bbR^9$.

Fields associated with eigenvalues $\lambda_{2,3,4,5,6,7,8}$ are linearly degenerate (LD). The $1$-field and $9$-field are genuinely non linear (GNL).
\end{prop}


We refer to~\cite{herard:hal-02007060,herard1994basic} for a proof.


\subsection{Wave structure}

Right eigenvectors are:
\begingroup
\allowdisplaybreaks
\begin{subequations}
\begin{align}
\vect{r}_{1} & = \left(1, \
\frac{R_{n\tau}}{R_{nn}}, \ - \sqrt{2 R_{nn}}, \ -\frac{\sqrt{2}R_{n\tau }^{2}}{(R_{nn})^{\frac{3}{2}}}, \ -R_{n\tau } \sqrt{\frac{2}{R_{nn}}},\nonumber\right. \\
&\mkern200mu\left. -\frac{\sqrt{2}(\rans{u_n'\theta'})}{\sqrt{R_{nn}}}, \ -\frac{\sqrt{2}(\rans{u_n'\theta'})R_{n\tau}}{(R_{nn})^{\frac{3}{2}}}, \ -\frac{\sqrt{2}(\rans{u_n'\theta'})^2}{(R_{nn})^{\frac{3}{2}}}, \
\frac{\rans{u'_n \theta'}}{R_{nn}} \right)^{T},
\\
\vect{r}_{2} & = \left(0, \ 0, \ 0, \ 0, \ 0, \ -\sqrt{R_{nn}}, \ -\frac{R_{n\tau}}{\sqrt{R_{nn}}}, \ -\frac{2\rans{u_n'\theta'}}{\sqrt{R_{nn}}}, \ 1
\right)^T, \\
\vect{r}_{3} & = \left(0, \ 1, \ 0, \ -\frac{2R_{n\tau}}{\sqrt{R_{nn}}}, \ -\sqrt{R_{nn}}, \ 0, \ -\frac{\rans{u_n'\theta'}}{\sqrt{R_{nn}}}, \ 0, \ 0 \right)^T,\\
\vect{r}_{4} & = \left(0, \ 0, \ 0, \ 1, \ 0, \ 0, \ 0, \ 0, \ 0 \right)^T,\\
\vect{r}_{5} & = \left(0, \ 0, \ 0, \ 0, \ 0, \ 0, \ 1, \ 0, \ 0 \right)^T,\\
\vect{r}_{6} & = \left(0, \ 0, \ 0, \ 0, \ 0, \ 0, \ 0, \ 1, \ 0\right)^T,\\
\vect{r}_{7} & = \left(0, \ 1, \ 0, \ \frac{2R_{n\tau}}{\sqrt{R_{nn}}}, \ \sqrt{R_{nn}}, \ 0, \ \frac{\rans{u_n'\theta'}}{\sqrt{R_{nn}}}, \ 0, \ 0 \right)^T,\\
\vect{r}_{8} & = \left(0, \ 0, \ 0, \ 0, \ 0, \ \sqrt{R_{nn}}, \ \frac{R_{n\tau}}{\sqrt{R_{nn}}}, \ \frac{2\rans{u_n'\theta'}}{\sqrt{R_{nn}}}, \ 1 \right)^T,\\
\vect{r}_{9} & = \left(1, \
\frac{R_{n\tau}}{R_{nn}}, \
\sqrt{2 R_{nn}}, \
\frac{\sqrt{2}R_{n\tau }^{2}}{(R_{nn})^{\frac{3}{2}}}, \ R_{n\tau } \sqrt{\frac{2}{R_{nn}}}, \
\frac{\sqrt{2}(\rans{u_n'\theta'})}{\sqrt{R_{nn}}}, \
\frac{\sqrt{2}(\rans{u_n'\theta'})R_{n\tau}}{(R_{nn})^{\frac{3}{2}}}, \
\frac{\sqrt{2}(\rans{u_n'\theta'})^2}{(R_{nn})^{\frac{3}{2}}}, \
\frac{\rans{u'_n \theta'}}{R_{nn}}
\right)^{T}. 
\end{align}
\end{subequations}
\endgroup



\subsection{Riemann invariants}
The Riemann invariants $I^{i}_{R}$ associated with the $i^{\rm th}$ field ($f \in I^i_R$, $\vect{\nabla} f \cdot \vect{r}_{i} = 0$) are:
\begingroup
\allowdisplaybreaks
\begin{subequations}\label{riemanninvariants_temp}
\begin{align}
I^{1}_{R} &= \Bigg\{ u_n+\sqrt{2R_{nn}}, \ \frac{R_{n\tau }}{R_{nn}}, \
\frac{R_{nn} R_{\tau \tau } - R_{n\tau }^2}{R_{nn}}, \nonumber\\
&\mkern150mu u_\tau + R_{n\tau }\sqrt{\frac{2}{R_{nn}}}, \
\frac{\delta_3}{R_{nn}}, \
\frac{\rans{u_n'\theta'} }{R_{nn}}, \
\frac{\rans{\theta'^2}R_{nn} - \rans{u_n'\theta'}^2}{ R_{nn}},\
\rans{\theta} + \rans{u_n'\theta'} \sqrt{\frac{2}{{R_{nn}}}}
\Bigg\},\label{eq:riemanninvariants_temp:a}
\\
I^{2,3}_{R} &= \Bigg\{ u_n, \ R_{nn}, \ u_\tau+\frac{R_{n\tau }}{\sqrt{R_{nn}}}, \ R_{\tau \tau }R_{nn} - R_{n\tau }^{2},\
\rans{\theta} + \frac{\rans{u_n'\theta'}}{\sqrt{{R_{nn}}}}, \
\frac{\rans{\theta'^2}R_{nn} - \rans{u_n'\theta'}^2}{ R_{nn}}, \
\Psi_{2,3}
\Bigg\},\label{eq:riemanninvariants_temp:b}
\\
I^{4,5,6}_{R} &= \Bigg\{ u_n, \ u_\tau, \ R_{nn}, \ R_{n\tau }, \
\rans{\theta}, \
\rans{u_n'\theta'}
\Bigg\},\label{eq:riemanninvariants_temp:c}
\\
I^{7,8}_{R} &= \Bigg\{ u_n, \ R_{nn}, \ u_\tau-\frac{R_{n\tau }}{\sqrt{R_{nn}}}, \ R_{\tau \tau }R_{nn} - R_{n\tau }^{2}, \
\rans{\theta} - \frac{\rans{u_n'\theta'}}{\sqrt{{R_{nn}}}}, \
\frac{\rans{\theta'^2}R_{nn} - \rans{u_n'\theta'}^2}{ R_{nn}}, \
\Psi_{7,8}
\Bigg\},\label{eq:riemanninvariants_temp:d}
\\
I^{9}_{R} &= \Bigg\{ u_n -\sqrt{2R_{nn}}, \ \frac{R_{n\tau }}{R_{nn}},
\frac{R_{nn} R_{\tau \tau } - R_{n\tau }^2}{R_{nn}}, \nonumber\\
&\mkern100mu u_\tau  - R_{n\tau } \sqrt{\frac{2}{R_{nn}}},\
\frac{\delta_3}{R_{nn}}, \
\frac{\rans{u_n'\theta'} }{R_{nn}}, \
\frac{\rans{\theta'^2}R_{nn} - \rans{u_n'\theta'}^2}{ R_{nn}},\
\rans{\theta} - \rans{u_n'\theta'} \sqrt{\frac{2}{{R_{nn}}}}
\Bigg\},\label{eq:riemanninvariants_temp:e}
\end{align}
\end{subequations}
\endgroup
\noindent with
\begin{align}
\Psi_{2,3} &\coloneqq \sqrt{R_{nn}}
\left(\rans{u_\tau'\theta'} + u_\tau \rans{\theta} \right) +
\rans{u_n'\theta'} u_\tau + R_{n \tau}
\rans{\theta},
\\
\Psi_{7,8} &\coloneqq - \sqrt{R_{nn}}
\left(\rans{u_\tau'\theta'} + u_\tau \rans{\theta} \right) +
\rans{u_n'\theta'} u_\tau + R_{n \tau}
\rans{\theta},
\\
\delta_3 &\coloneqq \rans{\theta'^2}R_{nn}R_{\tau\tau} + 2\left(\overline{u_n'\theta'}\right)\left(\overline{u'_\tau\theta'}\right)R_{n\tau} - R_{nn}\left(\overline{u'_\tau\theta'}\right)^2 - R_{\tau\tau}\left(\overline{u_n'\theta'}\right)^2 - \rans{\theta'^2} R_{n\tau}^2.
\end{align}



Approximate shock relations~\eqref{eqn:jump_lefloch} are supplemented with:
\begin{subequations}\label{eqn:jump_lefloch_theta}
\begin{align}
\sigma \left[\rans{\theta}\right] &= \mean{u_n}\left[\rans{\theta}\right] + \left[\rans{u_n'\theta'}\right],\\
\sigma \left[\rans{u_n'\theta'}\right] &= \mean{u_n}\left[\rans{u_n'\theta'}\right] + \mean{\rans{u_n'\theta'}}[u_n] + \mean{R_{nn}}\left[\rans{\theta}\right],\\
\sigma \left[\rans{u_\tau'\theta'}\right] &= \mean{u_n}\left[\rans{u_\tau'\theta'}\right] + \mean{\rans{u_n'\theta'}}[u_\tau] + \mean{R_{n\tau}}\left[\rans{\theta}\right],\\
\sigma \left[\rans{\theta'^2}\right] &= \mean{u_n}\left[\rans{\theta'^2}\right] + 2\mean{\rans{u_n'\theta'}}\left[\rans{\theta}\right].
\end{align}
\end{subequations}
These relations can also be obtained from pseudo-conservative variables $\rans{\theta^2} =\rans{\theta}^2 + \rans{\theta'^2} $ instead of $\rans{\theta'^2}$ and $\rans{\theta \vect{u}} =\rans{\theta} \, \rans{ \vect{u}}+ \rans{\theta' \vect{u}'} $ instead of $\rans{\theta' \vect{u}'}$.

\begin{theo}[Extension of \cref{exiunith} to anisothermal flows]\label{exiunith_bis}
Let
\begin{equation*}
\vect{W} \coloneqq \left(u_n,\, u_\tau, R_{nn}, \, R_{\tau \tau}, \, R_{n \tau}, \ \rans{u_n'\theta'}, \rans{u_\tau'\theta'}, \rans{\theta'^2}, \rans{\theta}\right),
\end{equation*}
and consider the {approximate} shock relations \cref{eqn:jump_lefloch,eqn:jump_lefloch_theta}. Assume strictly realizable initial conditions $\vect{W}(x_n, 0) = \vect{W}^L $ for $x_n <0$ and $\vect{W}(x_n, 0) = \vect{W}^R $ for $x_n > 0$, then there exists a unique self-similar \emph{realizable} solution: $\vect{W}(x_n, t) = \vect{w}(\frac{x_n}{t})$ to the Riemann problem associated with the evolution step involved in \cref{eq:evol-u-R,eq:evol-u-R-T} if and only if condition~\eqref{conditionexiuni} holds.
\end{theo}

\begin{proof}
The proof follows the main guidelines recalled in~\cite{smoller2012shock} and is only partially detailed in~\cite{herard:hal-02007060}. We consider two initial realizable states $\vect{W}^R$, $\vect{W}^L$. We look for the four intermediate states labelled $I$, $II$, $III$, $IV$ from left to right separating the 5 distinct waves.
\setcounter{stepcount}{-1}
\begin{step}{Preliminary calculation}
To connect left and right states $\vect{W}^{l,r}$ across the $1^{st}$ and $5^{th}$ waves we set $z \coloneqq \frac{R_{nn}^r}{R_{nn}^l} >0$. To select the physical solution for shocks~\eqref{eqn:jump_lefloch}, we use $\left[u_n\right] = u_n^r -u_n^l \leq 0$.

For the $1$-shock connection~($z>1$):
\begin{equation}\label{1shockconnexion_temp}
\begin{array}{l}
R_{nn}^r = zR_{nn}^l;\quad u_n^r = u_n^l + \frac{1-z}{\sqrt{z+1}}\sqrt{R_{nn}^l};\quad u_\tau^r = u_\tau^l + \frac{1-z}{\sqrt{z+1}}\frac{R_{n\tau}^l}{\sqrt{R_{nn}^l}};\quad
\\
R_{n\tau}^r = zR_{n\tau}^l;\quad R_{\tau\tau}^r = R_{\tau\tau}^l + (z-1)\frac{\left(R_{n\tau}^l\right)^2}{R_{nn}^l}.
\end{array}
\end{equation}

For the $9$-shock connection~($0<z<1$):
\begin{equation}\label{9shockconnexion_temp}
\begin{array}{l}
R_{nn}^r = zR_{nn}^l;\quad u_n^r = u_n^l + \frac{z-1}{\sqrt{z+1}}\sqrt{R_{nn}^l};\quad u_{\tau}^r = u_\tau^l + \frac{z-1}{\sqrt{z+1}}\frac{R_{n\tau}^l}{\sqrt{R_{nn}^l}};\quad
\\
R_{n\tau}^r = zR_{n\tau}^l;\quad R_{\tau \tau}^r = R_{\tau \tau}^l + (z-1)\frac{\left(R_{n\tau}^l\right)^2}{R_{nn}^l}.
\end{array}
\end{equation}

To connect the left and right states through the $1$-rarefaction wave (respectively $9$-rarefaction wave) we use the Riemann invariants~\eqref{eq:riemanninvariants:a} (respectively~\eqref{eq:riemanninvariants:e}).
\end{step}

\begin{step}{Solution in terms of $u_n$ and $R_{nn}$ variables}
A glance at~\eqref{eq:riemanninvariants:b}, \eqref{eq:riemanninvariants:c}, \eqref{eq:riemanninvariants:d} shows that $u_n$ and $R_{nn}$ are constant across the three LD waves, as for the isothermal case. First we can focus only on $u_n$ and $R_{nn}$ variables. We denote their intermediate values by $u_n^{\sharp}=u_n^{I}=u_n^{II}=u_n^{III}=u_n^{IV}$ and $R_{nn}^{\sharp}=R_{nn}^{I}=R_{nn}^{II}=R_{nn}^{III}=R_{nn}^{IV}$, and link these to initial states $\vect{W}^{R,L}$:
\begin{equation}\label{lefttoright_bis}
\left\{
\begin{aligned}
u_n^L &= u_n^{\sharp} - h_1(z_1)\sqrt{2R_{nn}^{L}}, \text{\ with $z_1 \coloneqq \frac{R_{nn}^{\sharp}}{R_{nn}^{L}}$} > 0,\\
u_n^R &= u_n^{\sharp} + h_5(z_9)\sqrt{2R_{nn}^{R}}, \text{\ with $z_9 \coloneqq \frac{R_{nn}^R}{R_{nn}^{\sharp}} > 0$},
\end{aligned}
\right.
\end{equation}
$h_{1,9}(z)$ are defined using~\eqref{eq:riemanninvariants_temp:a}, \eqref{eq:riemanninvariants_temp:e}, \eqref{1shockconnexion_temp}, \eqref{9shockconnexion_temp}, following~\cite{smoller2012shock}. Note that $z_1 z_9 = \frac{R_{nn}^R}{R_{nn}^{L}}$. We set:
\begin{equation}\label{eq:psi_temp}
\Psi(z_9) \coloneqq u_n^R -u_n^L - \sqrt{2R_{nn}^R}h_9(z_9) - \sqrt{2R_{nn}^L}h_1(z_1(z_9)) = 0.
\end{equation}
Solving~\eqref{eq:psi_temp} gives $z_9$, hence $R_{nn}^{\sharp}$ and $u_n^{\sharp}$. The solution is unique since $\Psi$ is a strictly monotonic function. Its lower and upper bounds enable to conclude that $u_n^{\sharp}$ and $R_{nn}^{\sharp}$ exist if~\eqref{conditionexiuni} holds. $z_1$ and $z_9$ enable to compute the remaining components $(u_\tau, \, R_{\tau \tau}, \, R_{n \tau}, \, \rans{u'_n\theta'}, \, \rans{u'_\tau \theta'}, \, \rans{\theta'^2}, \, \rans{\theta})$ of intermediate states $I$ and $IV$.
\end{step}

\begin{step}{Solution in terms of $u_\tau$, $R_{n\tau}$, and $\rans{\theta}$, $\rans{u'_n \theta'}$ variables}
As $u_\tau$, $R_{n\tau}$ are invariants of 4-5-6 wave, we set $u_\tau^\star = u_\tau^{II} = u_\tau^{III}$ and $R_{n\tau}^{\star}= R_{n\tau}^{II}= R_{n\tau}^{III}$, given by:
\begin{equation}\label{lefttoright_3_dyn}
\left\{
\begin{aligned}
u_\tau^\star + \dfrac{R_{n\tau}^{\star}}{\sqrt{R_{nn}^{\sharp}}} &= u_\tau^{I} + \dfrac{R_{n\tau}^{I}}{\sqrt{R_{nn}^{\sharp}}},\\
u_\tau^\star - \dfrac{R_{n\tau}^{\star}}{\sqrt{R_{nn}^{\sharp}}} &= u_\tau^{IV} - \dfrac{R_{n\tau}^{IV}}{\sqrt{R_{nn}^{\sharp}}}.
\end{aligned}
\right.
\end{equation}
Equation~\eqref{lefttoright_3_dyn} is obtained from Riemann invariants $I^{2,3}_R$, $I^{7,8}_R$.

And similarly for the thermal part, using $I^{4,5,6}_R$, we define $\rans{\theta}^\star = \rans{\theta}^{II} = \rans{\theta}^{III}$ and $\rans{u'_n \theta'}^{\star}= \rans{u'_n \theta'}^{II}= \rans{u'_n \theta'}^{III}$, obtained through:
\begin{equation}\label{lefttoright_3_thermal}
\left\{
\begin{aligned}
\rans{\theta}^\star + \dfrac{\rans{u'_n \theta'}^{\star}}{\sqrt{R_{nn}^{\sharp}}} &=
\rans{\theta}^{I} + \dfrac{\rans{u'_n \theta'}^{I}}{\sqrt{R_{nn}^{\sharp}}},\\
\rans{\theta}^\star - \dfrac{\rans{u'_n \theta'}^{\star}}{\sqrt{R_{nn}^{\sharp}}} &=
\rans{\theta}^{IV} - \dfrac{\rans{u'_n \theta'}^{IV}}{\sqrt{R_{nn}^{\sharp}}}.
\end{aligned}
\right.
\end{equation}
Equation~\eqref{lefttoright_3_thermal} is also obtained from Riemann invariants $I^{2,3}_R$, $I^{7,8}_R$.
\end{step}

\begin{step}{Solution in terms of $R_{\tau\tau}$, $\rans{\theta'^2}$ and $\rans{u'_\tau \theta'}$ variables.} Then, we deduce the remaining two variables $R_{\tau\tau}$ and $\rans{\theta'^2}$ for intermediate states $II$ (resp. $III$) from Riemann invariants of the $2-3$ wave (resp. $7-8$ wave):
\begin{equation}\label{lefttoright_Rtautau_theta2_utautheta}
\left\{
\begin{aligned}
R_{nn}^{\sharp} R_{\tau\tau}^{II} - R_{n\tau}^{\star, 2} &=
\left(R_{nn} R_{\tau\tau} - R_{n\tau}^{2} \right)^{I},\\
R_{nn}^{\sharp} R_{\tau\tau}^{III} - R_{n\tau}^{\star, 2} &=
\left(R_{nn} R_{\tau\tau} - R_{n\tau}^{2} \right)^{IV},\\
\rans{\theta'^2}^{II} - \dfrac{\rans{u'_n \theta'}^{\star, 2}}{R_{nn}^{\sharp}} &=
\left(\rans{\theta'^2} - \dfrac{\rans{u'_n \theta'}^{2}}{R_{nn} } \right)^{I},\\
\rans{\theta'^2}^{III} - \dfrac{\rans{u'_n \theta'}^{\star, 2}}{R_{nn}^{\sharp}} &=
\left(\rans{\theta'^2} - \dfrac{\rans{u'_n \theta'}^{2}}{R_{nn} } \right)^{IV}.
\end{aligned}
\right.
\end{equation}
The last component $\rans{u'_\tau \theta'}^{II} $ (resp. $\rans{u'_\tau \theta'}^{III} $) is uniquely deduced from the conservation of Riemann invariant $\Psi_{2,3}^{II}= \Psi_{2,3}^{I}$ (resp. $\Psi_{7,8}^{III}=\Psi_{7,8}^{IV}$).\qedhere
\end{step}\let\qed\relax
\end{proof}

One can note that the extension to 3-D for the isothermal case is analogous to the 2-D anisothermal case treated in this section.


\subsection{Numerical simulations}

\subsubsection{Near wall anisothermal Riemann problem} We present here the solution for the study of the Riemann problem associated to~\eqref{eq:evol-u-R-ntau_temp} including passive scalar transport, with initial realizable states:
\begin{equation*}
\left\{
\begin{aligned}
\bfW^L &=& \left(0, u_\tau^L, R_{nn}^L, R_{\tau\tau}^L, R_{n\tau}^L, \rans{u_n'\theta'}^L, \rans{u_\tau'\theta'}^L, \rans{\theta'^2}^L,
\rans{\theta}^L
\right)^T,\\
\bfW^R &=& \left(0, u_\tau^R, R_{nn}^L, R_{\tau\tau}^R, R_{n\tau}^R, \rans{u_n'\theta'}^R, \rans{u_\tau'\theta'}^R, \rans{\theta'^2}^R, \rans{\theta}^R\right)^T.
\end{aligned}
\right.
\end{equation*}
\begin{table}[h]
\caption{General solution for the tangential flow in the vicinity of a wall considering the transport of a scalar quantity.}
\vspace*{-10pt}
\begin{center}
\hspace*{-8pt}{
\begin{tabular}{|c || c | c |c | c | c | c |}
\hline zone & L & I & II & III & IV & R\\ \hline

$u_n$ & $0$ & $0$ & $0$ & $0$& $0$ & $0$\\ \hline

$u_\tau$ & $u_\tau^L$ & $u_\tau^L$ & $\frac{u_\tau^L + u_\tau^R}{2} + \frac{R_{n\tau}^L - R_{n\tau}^R}{2\sqrt{R_{nn}}}$ & $\frac{u_\tau^L + u_\tau^R}{2} + \frac{R_{n\tau}^L - R_{n\tau}^R}{2\sqrt{R_{nn}}}$ & $u_\tau^R$ & $u_\tau^R$\\ \hline

$R_{nn}$ & $R_{nn}$ & $R_{nn}$ & $R_{nn}$ & $R_{nn}$ & $R_{nn}$ & $R_{nn}$ \\ \hline

$R_{n\tau}$ & $R_{n\tau}^L$ & $R_{n\tau}^L$ & $\frac{\left(u_\tau^L -u_\tau^R\right) \sqrt{R_{nn}}}{2} + \frac{R_{n\tau}^L + R_{n\tau}^R}{2}$ & $\frac{\left(u_\tau^L -u_\tau^R\right) \sqrt{R_{nn}}}{2} + \frac{R_{n\tau}^L + R_{n\tau}^R}{2}$ & $R_{n\tau}^R$ & $R_{n\tau}^R$ \\ \hline

$R_{\tau\tau}$ & $R_{\tau\tau}^L$ & $R_{\tau\tau}^L$ & $R_{\tau\tau}^L + \frac{\left(R_{n\tau}^{II}\right)^2 - \left(R_{n\tau}^L\right)^2}{R_{nn}}$ & $R_{\tau\tau}^R + \frac{(R_{n\tau}^{III})^2 - \left(R_{n\tau}^R\right)^2}{R_{nn}}$ & $R_{\tau\tau}^R$ & $R_{\tau\tau}^R$ \\ \hline

$\rans{u_n'\theta'}$ & $\rans{u_n'\theta'}^L$ & $\rans{u_n'\theta'}^L$ & $\frac{\left(\rans{\theta}^L-\rans{\theta}^R\right)\sqrt{R_{nn}}}{2} + \frac{\rans{u_n'\theta'}^L + \rans{u_n'\theta'}^R}{2}$ & $\frac{\left(\rans{\theta}^L-\rans{\theta}^R\right)\sqrt{R_{nn}}}{2} + \frac{\rans{u_n'\theta'}^L + \rans{u_n'\theta'}^R}{2}$ & $\rans{u_n'\theta'}^R$ & $\rans{u_n'\theta'}^R$\\ \hline

$\rans{u_\tau'\theta'}$ & $\rans{u_\tau'\theta'}^L$ & $\rans{u_\tau'\theta'}^L$ & $\rans{u_\tau'\theta'}^{II}$ & $\rans{u_\tau'\theta'}^{III}$ & $\rans{u_\tau'\theta'}^R$ & $\rans{u_\tau'\theta'}^R$ \\ \hline

$\rans{\theta'^2}$ & $\rans{\theta'^2}^L$ & $\rans{\theta'^2}^L$ & $\rans{\theta'^2}^L + \frac{\left(\rans{u_n'\theta'}^{II}\right)^2 - \left(\rans{u_n'\theta'}^L\right)^2}{R_{nn}}$ & $\rans{\theta'^2}^R + \frac{\left(\rans{u_n'\theta'}^{III}\right)^2 - \left(\rans{u_n'\theta'}^R\right)^2}{R_{nn}}$ & $\rans{\theta'^2}^R$ & $\rans{\theta'^2}^R$\\ \hline

$\rans{\theta}$ & $\rans{\theta}^L$ & $\rans{\theta}^L$ & $\frac{\rans{\theta}^L + \rans{\theta}^R}{2} + \frac{\rans{u_n'\theta'}^L - \rans{u_n'\theta'}^R}{2\sqrt{R_{nn}}}$ & $\frac{\rans{\theta}^L + \rans{\theta}^R}{2} + \frac{\rans{u_n'\theta'}^L - \rans{u_n'\theta'}^R}{2\sqrt{R_{nn}}}$ & $\rans{\theta}^R$ & $\rans{\theta}^R$\\ \hline
\end{tabular}}
\end{center}
\vspace*{-5pt}
\end{table}


The numerical simulation of the current Riemann's problem is performed using \CS{}. The length of computational domain is $L=\SI{100}{m}$ and the total time of integration is $T=\SI{20}{s}$, and we take $\Delta x/L = \num{0.002}$, $\Delta t/T = 0.005$, which satisfy the CFL condition: $ \max (|\lambda_i|) \frac{\Delta t}{\Delta x} \leq \frac{1}{2} $.

The left and right state numerical values considered are: $u_n^L= u_n^R = \SI{0.}{\meter\per\second}$, $u_\tau^L= u_\tau^R = \SI{1.}{\meter\per\second}$, $R_{nn}^L=R_{nn}^R =R_{\tau\tau}^L=R_{\tau\tau}^R =\SI{0.5}{\square\meter\per\square\second}$, $R_{n\tau}^L= - R_{n\tau}^R = -\SI{0.1}{\square\meter\per\square\second}$, $\rans{\theta}^L = \rans{\theta}^R = \SI{1}{[\theta]}$, $\rans{u'_n\theta'}^L= - \rans{u'_n\theta'}^R = -\SI{0.3}{[\theta] \meter\per\second}$, $\rans{u'_\tau\theta'}^L= \rans{u'_\tau\theta'}^R = \SI{0.1}{[\theta] \meter\per\second}$ and $\rans{\theta'^2}^L = \rans{\theta'^2}^R = \SI{0.5}{[\theta]^2}$. Both left and right states are realizable and they satisfy the condition of \cref{exiunith_bis} for existence and uniqueness of the solution.

The results of the simulation where we applied the Rusanov-like scheme are presented in \cref{Figure:TestCase3}. As previously we observe the absence of oscillations due to the application of the Rusanov-like scheme.
\cref{fig:convergence_Riemann_anisothermal} shows a one-half convergence rate, as expected due to the contact discontinuities. One can also notice that the Rusanov-like scheme does not recover the exact constant state for $R_{\tau \tau}$, $\rans{u'_\tau\theta'}$ and $\rans{\theta'^2}$, similarly to what one gets on pressure profiles with Rusanov scheme for Euler equations with perfect gas equation of state. However this artefact disappears with space-time convergence.
\newpage
\begin{figure}[!ht]
\centering
\begin{subfigure}{.5\textwidth}
\centering
\includegraphics[width=70mm]{figures/profile_vz_rus1-crop.pdf}
\caption{}
\end{subfigure}\hfill
\begin{subfigure}{.5\textwidth}
\centering
\includegraphics[width=70mm]{figures/profile_vx_rus1-crop.pdf}
\caption{}
\end{subfigure}\\
\begin{subfigure}{.5\textwidth}
\centering
\includegraphics[width=70mm]{figures/profile_rzz_rus1-crop.pdf}
\caption{}
\end{subfigure}\hfill
\begin{subfigure}{.5\textwidth}
\centering
\includegraphics[width=70mm]{figures/profile_rxx_rus1-crop.pdf}
\caption{}
\end{subfigure}\\
\begin{subfigure}{.5\textwidth}
\centering
\includegraphics[width=70mm]{figures/profile_rxz_rus1-crop.pdf}
\caption{}
\end{subfigure}\hfill
\begin{subfigure}{.5\textwidth}
\centering
\includegraphics[width=70mm]{figures/profile_T_rus1-crop.pdf}
\caption{}
\end{subfigure}\\
\begin{subfigure}{.5\textwidth}
\centering
\includegraphics[width=70mm]{figures/profile_uT_rus1-crop.pdf}
\caption{}
\end{subfigure}\hfill
\begin{subfigure}{.5\textwidth}
\centering
\includegraphics[width=70mm]{figures/profile_vT_rus1-crop.pdf}
\caption{}
\end{subfigure}
\end{figure}
\FloatBarrier
\pagebreak
\begin{figure}[!ht]
\setcounter{figure}{0}
\centering
\begin{subfigure}{.5\textwidth}
\setcounter{subfigure}{8}
\centering
\includegraphics[width=70mm]{figures/profile_Tvar_rus1-crop.pdf}
\caption{}
\end{subfigure}\hfill
\begin{subfigure}{.5\textwidth}
\setcounter{subfigure}{9}
\centering
\includegraphics[width=70mm]{figures/convergence_plot-crop.pdf}
\caption{\label{fig:convergence_Riemann_anisothermal}
}
\end{subfigure}
\vspace*{-10pt}
\caption{\label{Figure:TestCase3}
Profile of the variables $u_n$ (a.),$u_\tau$ (b.), $R_{nn}$ (c.), $R_{\tau\tau}$ (d.), $R_{n\tau}$ (e.), $\theta$ (f.), $\rans{u'_n \theta'}$ (g.), $\rans{u'_\tau \theta'}$ (h.), $\rans{\theta'^2}$ (i.) along with convergence graph (j.) obtained applying the Rusanov-like scheme.}
\end{figure}

\begin{figure}[!h]
\begin{subfigure}{\textwidth}
\centering
\includegraphics[width=\textwidth]{figures/slice_c_legacy.png}
\caption{Material-upwinding scheme.}\label{fig:pollutantcliclac}
\end{subfigure}
\begin{subfigure}{\textwidth}
\centering
\includegraphics[width=\textwidth]{figures/slice_c_rusanov.png}
\caption{Rusanov-like scheme.}\label{fig:pollutantcliclacerased}
\end{subfigure}
\vspace*{-10pt}
\caption{Point-source dispersion of a pollutant using a material-upwinding scheme and the Rusanov-like scheme.\label{fig:point_source}
}
\end{figure}

\subsubsection{Point-source dispersion of a pollutant}
We present here the steady state of a point-source dispersion of a passive pollutant in a steady and homogeneous dynamic. The following hypotheses were made:
\begin{itemize}
\item the pollutant is introduced using a Dirac mass and continuously in time at a concentration value of $1$ (computationnaly, it is introduced in a single cell),
\item molecular diffusion is neglected in front of turbulent diffusion,
\end{itemize}

One can observe in \cref{fig:point_source} the results obtained by applying material-upwind scheme which present spurious spatial steady oscillations, and are suppressed by the Rusanov-like scheme.

\begin{figure}[ht]
\centering
\includegraphics[scale=1.6]{figures/Fig3.png} 
%\begin{subfigure}{.5\textwidth}
%\centering
%\includegraphics[width=70mm]{figures/profile_vx_rus_nx_500-crop.pdf}
%\end{subfigure}\hfill
%\begin{subfigure}{.5\textwidth}
%\centering
%\includegraphics[width=70mm]{figures/profile_vy_rus_nx_500-crop.pdf}
%\end{subfigure}\\
%\begin{subfigure}{.5\textwidth}
%\centering
%\includegraphics[width=70mm]{figures/profile_rxx_rus_nx_500-crop.pdf}
%\end{subfigure}\hfill
%\begin{subfigure}{.5\textwidth}
%\centering
%\includegraphics[width=70mm]{figures/profile_ryy_rus_nx_500-crop.pdf}
%\end{subfigure}\\
%\begin{subfigure}{.5\textwidth}
%\centering
%\includegraphics[width=70mm]{figures/profile_rxy_rus_nx_500-crop.pdf}
%\end{subfigure}\hfill
%\begin{subfigure}{.5\textwidth}
%\centering
%\includegraphics[width=70mm]{figures/convergence_plot-crop.pdf}
%\end{subfigure}
\caption{\label{Figure:simulation}
Profiles of the variables $u_n$, $u_\tau$, $R_{nn}$, $R_{\tau\tau}$, and $R_{n\tau}$ along with convergence plots obtained by applying material-upwind scheme and Rusanov-like scheme ($500$ cells, CFL$=0.5$).}
\end{figure}


\section{Conclusion}
A numerical strategy has been proposed to compute incompressible flows with second-order moment turbulence closures, in order to remove numerical oscillations. It is based on the analysis of the full convective part of the system. As in the compressible framework, the realizability of the Reynolds stress tensor is mandatory to obtain the hyperbolicity of the evolution system. The numerical convergence has been demonstrated on a verification test case using a Rusanov-like scheme. It is also suitable for anisothermal flows or mean species transport equations with second-moment turbulence closures.

Of course, more accurate Riemann solvers, or relaxation solvers, and second order extension can be investigated. This methodology can be applied to projection step methods or Uzawa algorithms.

Eventually, more complex models involving non-zero contribution $\tens{\Phi}^r(\gradt \vect{\rans{u}}, \, \tens{R})$ may be considered, using the same approach (see for instance~\cite[Appendix~2]{herard:hal-02007060}).




\section*{Conflicts of interest}

The authors declare no competing financial interest.

\section*{Dedication}

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.



\bibliographystyle{crunsrt}
\bibliography{crmeca20230408}
\end{document}
