Frederic Dehoux

Yannick Lecocq

Sofiane Benhamadouche
email: sofiane.benhamadouche@edf.fr

Remi Manceau

Laurent-Emmanuel Brizzi

Laurent-Emmanuel Brizzi

L.-E Brizzi

Algebraic modeling of the turbulent heat fluxes using the elliptic blending approach. Application to forced and mixed convection regimes

Keywords: Turbulent heat fluxes, Elliptic blending, Thermal length scale, Forced and mixed convection

Introduction

Many industrial applications, in particular in the field of energy production, are still treated with RANS eddy-viscosity models and the Simple Gradient Diffusion Hypothesis (SGDH), based on a turbulent Prandtl number, to model the turbulent heat fluxes in the Reynolds-averaged temperature equation. These simple models have several limitations, either for the prediction of the dynamic of the flow or for the estimation of the Nusselt number distribution. This can be due to the inherent complexity of the dynamics of the flow, e.g., in impingement or swirling regions, as well as to complex coupling mechanisms, such as buoyancy effects.

The present article focuses on the so-called low-Reynolds number modeling for both the Reynolds stresses and the turbulent heat fluxes. In order to avoid the use of damping functions to replicate the effect of the presence of solid boundaries, based on a non-dimensional distance to the wall and/or a local turbulent Reynolds number, Durbin [START_REF] Durbin | Near-wall turbulence closure modeling without "damping functions[END_REF] introduced the Elliptic Relaxation approach. A simplified version, the Elliptic Blending approach, was proposed by Manceau and Hanjalić [START_REF] Manceau | An improved version of the elliptic blending model. Application to non-rotating and rotating channel flows[END_REF][START_REF] Manceau | Elliptic blending model: a new near-wall Reynolds-stress turbulence closure[END_REF], leading to the so-called Elliptic Blending Reynolds-Stress Model (EB-RSM), which was successfully applied to several configurations, in particular in isothermal and forced convection flows, in the context of RANS [START_REF] Borello | Prediction of cascade flows with innovative second-moment closures[END_REF][START_REF] Fadai-Ghotbi | Revisiting URANS computations of the backwardfacing step flow using second moment closures. Influence of the numerics[END_REF][START_REF] Han | Turbulent heat transfer of supercritical carbon dioxide in square cross-sectional duct flow[END_REF][START_REF] Ho | Compressibility factor effect on the turbulence heat transfer of super-critical carbon dioxide by an elliptic-blending second moment closure[END_REF][START_REF] Shin | Modeling of turbulent heat transfer in an axially rotating pipe flow[END_REF][START_REF] Shin | Numerical analysis of turbulent flow and heat transfer in a square sectioned U-bend duct by elliptic-blending second moment closure[END_REF][START_REF] Shin | Refinement of a second moment closure by the elliptic blending equation and its application to turbulent rotating channel flows[END_REF][START_REF] Thielen | Predictions of flow and heat transfer in multiple impinging jets with an elliptic-blending second-moment closure[END_REF][START_REF] Törnblom | A Reynolds stress closure description of separation control with vortex generators in a plane asymmetric diffuser[END_REF][START_REF] Viti | Numerical study of stress-transport turbulence models: Implementation and validation issues[END_REF], as well as, recently, hybrid RANS/LES [START_REF] Fadai-Ghotbi | A seamless hybrid RANS-LES model based on transport equations for the subgrid stresses and elliptic blending[END_REF][START_REF] Fadai-Ghotbi | Temporal filtering: a consistent formalism for seamless hybrid RANS-LES modeling in inhomogeneous turbulence[END_REF] (see also the review articles [START_REF] Durbin | Limiters and wall treatments in applied turbulence modeling[END_REF][START_REF] Gatski | Current trends in modeling research for turbulent aerodynamic flows[END_REF]). This model is used in the present article to model the Reynolds stresses.

The modeling of the turbulent heat fluxes has attracted less efforts from the research community. As a consequence, several industrial computations are carried out with sophisticated Reynolds Stress Models but still use a Simple Gradient Diffusion Hypothesis for the heat fluxes, which is not satisfactory in mixed and natural convection regimes (e.g., see Kenjereš et al. [START_REF] Kenjereš | Contribution to elliptic relaxation modelling of turbulent natural and mixed convection[END_REF]). The Generalized Gradient Diffusion Hypothesis (GGDH) [START_REF] Ince | On the computation of buoyancy-driven turbulent flows in rectangular enclosures[END_REF] accounts for the influence of the anisotropy of the Reynolds stresses on turbulent heat fluxes, but is not sufficiently general to correctly represent buoyancy effects [START_REF] Hanjalić | One-point closure models for buoyancy-driven turbulent flows[END_REF]. To overcome this problem, Algebraic Flux Models (AFMs) [START_REF] Dol | A comparative assessment in the second-moment differential and algebraic models in turbulent natural convection[END_REF][START_REF] Hanjalić | Natural convection in partitioned two-dimensional enclosures at higher Rayleigh numbers[END_REF][START_REF] Kenjereš | Prediction of turbulent thermal convection in concentric and eccentric horizontal annuli[END_REF][START_REF] Launder | On the effects of a gravitational field on the turbulent transport of heat and momentum[END_REF] preserves the different production mechanisms arising in the exact transport equation for the turbulent heat fluxes.

However, in the existing AFMs, the effects of the walls are not taken into account. Similarly to the case of the Reynolds stresses mentioned above, in order to avoid the use of damping functions, the present paper aims at introducing the elliptic blending approach in algebraic flux models. Such models can be derived from the model recently proposed by Shin et al. [START_REF] Shin | Elliptic relaxation second moment closure for the turbulent heat fluxes[END_REF], based on differential equations for the turbulent heat fluxes (EB-DFM, Elliptic Blending-Differential Flux Model). Therefore, the present paper aims at -providing an in-depth justification for using elliptic relaxation/blending for the turbulent heat fluxes, based on a theoretical derivation using Green's formalism and an analysis of DNS budgets, -investigating the validity of the simplifying hypothesis used in Shin et al. [START_REF] Shin | Elliptic relaxation second moment closure for the turbulent heat fluxes[END_REF], which consists in using the same length scales in the elliptic blending equations for the Reynolds stresses and the turbulent heat fluxes, -deriving Elliptic Blending-Algebraic Flux Models from the Elliptic Blending-Differential Flux Model and validating them in forced and mixed convection cases.

The Elliptic Blending Approach

In this section, the main ideas leading to the elliptic blending approach for the velocity-pressure gradient correlation in the Reynolds stress transport equation are first recalled and its validity in cases with buoyancy is investigated. Finally, the extension of the approach to the modeling of the temperature pressure-gradient correlation term that enters the turbulent heat fluxes equations is presented.

Modeling the effect of the wall on the Reynolds-stresses

In a Reynolds Stress Model (RSM), in order to replicate the blocking of the wall-normal fluctuations, it can be shown, using a simple model for the two-point correlations and Green's formalism [START_REF] Durbin | A Reynolds stress model for near-wall turbulence[END_REF][START_REF] Manceau | Inhomogeneity and anisotropy effects on the redistribution term in Reynolds-averaged Navier-Stokes modelling[END_REF], that the elliptic relaxation model

φ * ij -ε ij -L 2 ∇ 2 φ * ij -ε ij = (φ h ij -ε h ij ) (1)
can be used for the difference between the velocity-pressure-gradient correlation term φ * ij and dissipation term ε ij . The right hand side, φ h ij -ε h ij , can be modeled using any available high-Reynolds number model. The main feature of this model is the satisfaction of the asymptotic, near-wall balance

φ * ij -ε ij = -D ν ij , (2)
where D ν ij denotes the molecular diffusion tensors. In order to avoid the resolution of six, strongly coupled and numerically unstable differential equations, a simplified version, the Elliptic Blending-RSM, was proposed by Manceau and Hanjalić [START_REF] Manceau | Elliptic blending model: a new near-wall Reynolds-stress turbulence closure[END_REF] and improved by Manceau [START_REF] Manceau | An improved version of the elliptic blending model. Application to non-rotating and rotating channel flows[END_REF], based on a blending of any standard model for φ h ij -ε h ij , valid far from the wall, and a near-wall model for φ w ij -ε w ij that satisfies the asymptotic balance (Eq. 2) in the near-wall region,

φ * ij -ε ij = (1 -α 3 )(φ w ij -ε w ij ) + α 3 (φ h ij -ε h ij ). ( 3
)
In order to preserve the non-local character of the blocking effect, similarly to Durbin's elliptic relaxation model, the blending function α is obtained from an elliptic relaxation equation

α -L 2 ∇ 2 α = 1 with α| w = 0 , (4)
and goes from 0 at the wall to 1 far from the wall. If Rotta's model [START_REF] Rotta | Statistische theorie nichthomogener turbulenz[END_REF] (ε w ij = u i u j /k ε) is used, the asymptotic analysis in the vicinity of the wall shows that the following model for φ w ij is consistent with the wall limits

φ w ij = -5 ε k u i u k n j n k + u j u k n i n k - 1 2 u k u l n k n l n i n j + δ ij ,
where n is a unit vector

n = ∇α/ ∇α , ( 5
)
providing a generalized wall-normal direction. Far from the wall, the SSG model (Speziale et al. [START_REF] Speziale | Modeling the pressure-strain correlation of turbulence: an invariant dynamical system approach[END_REF]) is used for the redistribution term and an isotropic model for the dissipation tensor (ε h ij = 2/3 ε δ ij ). In cases where buoyancy is present, the issue of the validity of the EB-RSM, which is based on the satisfaction of the near-wall balance (Eq. 2), arises. Close to a wall located at y = 0, Taylor-series expansions lead to the behaviors given in Table 1 for the Reynolds stress budgets when the temperature fluctuation is zero at the wall, which is the case in the present article. It can be seen that the production

D ν ij φ * ij ε ij P ij D T ij G ij u 2 O(1) O(y) O(1) O(y 3 ) O(y 3 ) O(y 2 ) v 2 O(y 2 ) O(y 2 ) O(y 2 ) O(y 5 ) O(y 5 ) O(y 3 ) w 2 O(1) O(y) O(1) O(y 3 ) O(y 3 ) O(y 2 ) u v O(y) O(y) O(y) O(y 4 ) O(y 4 ) O(y 2 ) u w O(1) O(y) O(1) O(y 3 ) O(y 3 ) O(y 2 ) v w O(y) O(y) O(y) O(y 4 ) O(y 4 ) O(y 2 )
due to buoyancy G ij = -β u i θ g j + u j θ g i does not affect the balance given by Eq. 2, and, in particular, that G ij /φ * ij = O(y) for all i and j. It can be concluded that the EB-RSM does not require any further modification while introducing buoyancy terms. However, a modification of the transition from the near-wall behavior to the far-from-the-wall behavior can be expected, which would question the validity of Eq. 4. Nevertheless, using the DNS data of Kasagi and Nishimura [START_REF] Kasagi | Direct numerical simulation of combined forced and natural turbulent convection in a vertical plane channel[END_REF], in the mixed convection case of a vertical channel with a pressure gradient and a temperature difference between the two walls (a more detailed description of this case is provided in Section 4), it can be seen in Fig. 1 that the blending function α is only marginally modified by buoyancy. It is worth recalling here that the model is formulated to reproduce the wall-blocking effect, i.e., to specifically damp the redistribution of energy toward the wall-normal fluctuations. Therefore, the evaluation of the blending parameter α and the length scale L to be reproduced by the model can be performed by using the components φ * 22 and ε 22 of the DNS databases (a priori evaluation), by computing

α 3 = (φ * 22 -ε 22 ) -(φ w 22 -ε w 22 ) (φ h 22 -ε h 22 ) -(φ w 22 -ε w 22 ) , ( 6
)
Similar to the isothermal case, it can be seen that α varies between 0 and 1, and even though it becomes slightly asymmetrical, the slope in the near-wall region (i.e., the extent of the influence of the wall) is only weakly modified. This modification and the asymmetry can be explained (and thus reproduced) by the influence of buoyancy on turbulence, and thus on the length scale L = C L max k 3/2 /ε, C η ν 3/4 /ε 1/4 . Therefore, it is considered that the EB-RSM does not require any modification when buoyancy effects are present.

2.2 Extension to the modeling of the effect of the wall on the heat fluxes As it is the case for the redistribution term in the Reynolds stress transport equations, a model for the scrambling term in the transport equation for the turbulent heat fluxes can be derived from the Poisson equation of the fluctuating pressure.

The momentum and continuity equations are written in the RANS form with a Boussinesq approximation

∂u j ∂ x j = 0 ( 7
)
∂u i ∂t + ∂ u i u j ∂ x j = - 1 ρ 0 ∂ p ∂ x i + ν ∂ 2 u i ∂ x j ∂ x j - ∂ R ij ∂ x j + ρ ρ 0 g i ( 8
)
where u i are the velocity components, p the pressure, R ij the Reynolds stresses, ν the kinematic viscosity, ρ the (temperature-dependent) density and ρ 0 a reference density.

The mean temperature equation is given by

∂θ ∂t + ∂ u j θ ∂ x j = ∂ ∂ x j ν Pr ∂θ ∂ x j -u j θ (9)
A Poisson equation can be obtained for the gradient of the fluctuating pressure

∇ 2 ∂ p ∂ x k = -ρ 0 ∂ ∂ x k 2 ∂u i ∂ x j ∂u j ∂ x i + ∂u i ∂ x j ∂u j ∂ x i - ∂u i ∂ x j ∂u j ∂ x i - 1 ρ 0 g l ∂ρ ∂ x l ( 10
)
Using the Green function of the domain and assuming that ∂ p ∂ x k satisfies a Neuman boundary condition (see Manceau et al. [START_REF] Manceau | Inhomogeneity and anisotropy effects on the redistribution term in Reynolds-averaged Navier-Stokes modelling[END_REF]), the solution of Eq. 10 can be obtained in the integral form

∂ p ∂ x k (x) = ∇ 2 ∂ p ∂ x k (y)G(x, y)dy
The scrambling term (temperature pressure-gradient one-point correlation) that arises in the turbulent heat flux transport equation is

φ * θi = - 1 ρ 0 θ ∂ p ∂ x i , such that ρ 0 φ * θi = θi (x, y)G(x, y)dy
where θi is the two-point correlation between the temperature and the Laplacian of the pressure gradient:

θi (x, y) = -θ (x)∇ 2 ∂ p ∂ x k (y)
Similarly to what was proposed by Durbin [START_REF] Durbin | Near-wall turbulence closure modeling without "damping functions[END_REF] and justified by Manceau et al. [START_REF] Manceau | Inhomogeneity and anisotropy effects on the redistribution term in Reynolds-averaged Navier-Stokes modelling[END_REF] for the velocity-pressure gradient two-point correlation, an exponential decay is assumed

θi (x, y) = θi (y, y)e -r L θ , (11)
where r = xy , and the correlation length scale L θ is to be determined. Using this assumption, the scrambling term reads

ρ 0 φ * θi = 1 4π θi (y, y) exp (-r/L θ ) r dy . Since 1 4π exp(-r/Lθ ) r
is the Green's function associated to the operator

-∇ 2 + 1 L 2 θ Id, an elliptic relaxation equation is obtained for φ * θi φ * θi -L 2 θ ∇ 2 φ * θi = - 1 ρ 0 L 2 θ θ ∂ p ∂ x i .
In a homogeneous flow, this equation reduces to

φ * θi = - 1 ρ 0 L 2 θ θ ∂ p ∂ x i
Consequently, the right hand side of the equation can be modeled using any quasihomogeneous model, such that the elliptic relaxation provides the appropriate nearwall treatment for standard high Reynolds models

φ * θi -L 2 θ ∇ 2 φ * θi = φ h θi . ( 12
)
The above derivation justifies the use of the elliptic relaxation approach to the turbulent heat fluxes, as proposed by Shin et al. [START_REF] Shin | Elliptic relaxation second moment closure for the turbulent heat fluxes[END_REF] based on a simple analogy with the Reynolds stresses. Instead of solving three equations for individual components of the heat fluxes and following Manceau and Hanjalić [START_REF] Manceau | Elliptic blending model: a new near-wall Reynolds-stress turbulence closure[END_REF], Shin et al. [START_REF] Shin | Elliptic relaxation second moment closure for the turbulent heat fluxes[END_REF] extended the elliptic blending approach to the heat fluxes, leading to the Elliptic Blending-Differential Flux Model (EB-DFM)

φ * θi -ε θi = (1 -α 3 θ ) φ w θi -ε w θi + α 3 θ φ h θi -ε h θi . ( 13
)
The parameter α θ , used to blend the near-wall and far-from-the-wall forms of the heat fluxes model is obtained from the elliptic relaxation equation

α θ -L 2 θ ∇ 2 α θ = 1 (14)
with the boundary condition α θ = 0 at the wall. This equation ensures that α θ → 1 far from the wall. Contrary to Shin et al. [START_REF] Shin | Elliptic relaxation second moment closure for the turbulent heat fluxes[END_REF], it is not assumed here that the length scale L θ is equal to the dynamic length scale L used in Eq. 4, i.e., for the difference φ * ij -ε ij in the Reynolds-stress transport equations, which is, following Durbin [START_REF] Durbin | Near-wall turbulence closure modeling without "damping functions[END_REF], modeled by

L = C L max k 3/2 ε , C η ν 3/4 ε 1/4 (15)
with C L = 0.125 and C η = 80. Table 2 gives the asymptotic behavior of the different terms on the r.h.s. of the transport equations for the turbulent heat fluxes for zero temperature fluctuations at the wall (see Eq. 24). Similar to the case of the Reynolds stresses, the near-wall budget of the wall normal heat flux reduces, at leading order, to a balance between scrambling, dissipation and molecular diffusion,

φ * θ 2 -ε θ 2 + D ν θ 2 = 0.
In order to ensure the correct asymptotic behavior at the wall of φ w θi -ε w θi , Shin et al. [START_REF] Shin | Elliptic relaxation second moment closure for the turbulent heat fluxes[END_REF] showed that the near-wall terms can be modeled as

φ w θi = -1 + 1 2 1 + 1 Pr ε k u j θ n i n j ( 16
)
and

ε w θi = 1 2 1 + 1 Pr ε k u i θ . ( 17
)
Note that this model does not satisfy the individual asymptotic behaviors of φ w θi and ε w θi , but only of their difference. This problem can be simply overcome by using

φ w θi = - ε k u j θ n i n j ( 18
)
and

ε w θi = 1 2 1 + 1 Pr ε k u i θ + u j θ n i n j , ( 19
)
which is, in practice, strictly equivalent to the previous model, since the difference φ w θi and ε w θi is the same. Far from the wall, the model 13 tend to the high-Reynolds number model φ h θi -ε h θi . Shin et al. [START_REF] Shin | Elliptic relaxation second moment closure for the turbulent heat fluxes[END_REF] proposed to use the standard linear model [START_REF] Launder | On the computation of convective heat transfer in complex turbulent flows[END_REF],

φ h θi = -C θ 1 ε k u i θ + C θ 2 u j θ ∂u i ∂ x j + C θ 2 u i u j ∂θ ∂ x j + C θ 3 g i βθ 2 (20)
and

ε h θi = 0 . ( 21
)
Table 2 Asymptotic behaviors of the terms of the heat fluxes budget (Shin et al. [START_REF] Shin | Elliptic relaxation second moment closure for the turbulent heat fluxes[END_REF])

P U θi + P T θi G θi φ * θi D ν θi -ε θi D t θi u θ O(y 3 ) O(y 2 ) O(y) O(1) O(1) O(y 3 ) v θ O(y 4 ) O(y 2 ) O(y) O(y) O(y) O(y 4 ) w θ O(y 3 ) O(y 2 ) O(y) O(1) O(1) O(y 3 )
y denotes the distance to the wall In order to write equations as general as possible, the high-Reynolds number model φ h θi for the pressure scrambling term used here is the nonlinear model of Kenjereš et al. [START_REF] Kenjereš | Contribution to elliptic relaxation modelling of turbulent natural and mixed convection[END_REF], consisting of the rapid part

φ r θi = C θ 2 u j θ ∂u i ∂ x j + C θ 2 u i u j ∂θ ∂ x j + C θ 3 g i βθ 2 (22)
and the slow part

φ s θi = -C θ 1 ε k u i θ -C θ 1 a ij u j θ ( 23
)
where a ij is the stress-anisotropy tensor,

a ij = R ij /k -2/3 δ ij .
The linear high-Reynolds number model (Eq. 20) is obviously a particular case of this model.

Algebraic Modeling of the Heat Fluxes

Although differential flux modeling is the appropriate level to represent the physical mechanisms driving the evolution of the turbulent heat fluxes, industrial application still rely on much simpler models. The aim of the present section is to make available an intermediate model based on the elliptic blending method to represent the nearwall effects, i.e., an algebraic flux model derived from the EB-DFM presented above.

After recalling the method used to derive the standard AFM, this method is applied to the EB-DFM, leading to the Elliptic Blending-Algebraic Flux Model (EB-AFM). Similar to the case of GGDH, which can be written as a simplification of the AFM, a simple model, the EB-GGDH (Elliptic Blending-Generalized Gradient Diffusion Hypothesis), is readily obtained.

The standard AFM

Using the Boussinesq approximation, the transport equations for the turbulent heat fluxes can be written

Du i θ Dt Cθi = -u j θ ∂u i ∂ x j P U θi -u i u j ∂θ ∂ x j P T θi -g i βθ 2 G θi - θ ρ ∂ p ∂ x i φ * θi -(λ + ν) ∂θ ∂ x j ∂u i ∂ x j ε θi + ∂ ∂ x j -θ u i u j D t θi + ∂ ∂ x j λ ∂θ ∂ x j u i + νθ ∂u i ∂ x j D ν θi ( 24
)
where C θi , P U θi , P T θi , G θi , φ * θi , ε θi , D t θi and D ν θi denote the material derivative, the production by the mean velocity gradient, by the mean temperature gradient, by buoyancy, the temperature-pressure gradient correlation, the dissipation, the turbulent and molecular diffusions, respectively.

Algebraic forms of the turbulent heat flux transport equation can be derived, similarly to the Reynolds stress transport equations (see, for example, [START_REF] Oceni | Introduction of wall effects in explicit algebraic stress models through elliptic blending[END_REF]).

As detailed by Hanjalić [START_REF] Hanjalić | One-point closure models for buoyancy-driven turbulent flows[END_REF], algebraic models for the heat fluxes rely on the weak equilibrium assumptions for convection

du i θ dt = 1 2 1 k dk dt + 1 θ 2 dθ 2 dt u i θ , ( 25
)
and diffusion

D θi = 1 2 1 k D k + 1 θ 2 D θ 2 u i θ , ( 26
)
where D θi , D k and D θ 2 denote the total diffusion of u i θ , k and θ 2 , respectively, which yields the turbulent heat flux algebraic equation

P θi + φ * θi -ε θi - u i θ 2k (P k + G k -ε) - u i θ 2θ 2 P θ 2 -ε θ 2 = 0 ,
where P θi = P U θi + P T θi + G θi . If the production and the dissipation terms of both k and θ 2 are assumed to be locally in balance, i.e., turbulence is assumed to be in local equilibrium,

P k + G k = ε ( 27
)
and

P θ 2 = ε θ 2 , ( 28
)
then Eq. 27 reduces to

P θi + φ * θi -ε θi = 0 . ( 29
)
Introducing the nonlinear model (Eqs. 22 and 23) for the scrambling term in this relation yields

C θ1 ε k u i θ = (1 -C θ 2 ) P U θi + 1 -C θ 2 P T θi + (1 -C θ 3 ) G θi + C θ 1 C θ 1 ε k a ij u j θ -ε h θi ( 30
)
For Prandtl numbers not small compared to unity, temperature fluctuations at small scales can be considered isotropic [START_REF] Tennekes | A First course in Turbulence[END_REF], such that the dissipation term is assumed to be zero far from the wall (ε h θi = 0). Equation 24 thus reduces to the Non-Linear Algebraic Flux Model (NL-AFM), Kenjereš et al. [START_REF] Kenjereš | Contribution to elliptic relaxation modelling of turbulent natural and mixed convection[END_REF],

u i θ = -C θ k ε ζ u i u j ∂θ ∂ x j + ξ u j θ ∂u i ∂ x j + ηβg i θ 2 -χ ε k a ij u j θ ( 31
)
where

C θ = C θ /C θ1 , ζ = 1 -C θ 2 , ξ = 1 -C θ 2 , η = 1 -C θ 3 and χ = C θ 1 C θ1 .
Two simplifications of this model will be utilized in the following: the linear Algebraic Flux Model (AFM, χ = 0) and the Generalized Gradient Diffusion Hypothesis (GGDH, ξ = η = χ = 0). Throughout the present paper, the coefficients that are generally recommended [START_REF] Hanjalić | One-point closure models for buoyancy-driven turbulent flows[END_REF], C θ 1 = 3.0, C θ 2 = 0, C θ 2 = C θ 3 = 0.55, are used . Note that the coefficient C θ , which should be equal to 1/C θ 1 , is modified by the introduction of the coefficient C θ ; indeed, the assumptions used in the algebraic methodology make necessary the introduction of this recalibration coefficient.

Elliptic blending-algebraic flux model

The model for the scrambling term given by Eqs. 22 and 23, and consequently the AFM, does not account for the effect of the wall on turbulence. In order to overcome this limitation, the main proposal of the present paper is to apply the elliptic blending strategy to the heat fluxes in order to derive the Elliptic Blending-Algebraic Flux Model (EB-AFM).

The same derivation as the one done in the Section 3.1 for the high-Reynolds number model but now applied to the EB-DFM described in Section 2.2 then leads to

P U θi + P T θi + G θi - u i θ 2k (P k + G k ) - u i θ 2θ 2 P θ 2 -(1 -α 3 θ ) ε k u j θ n i n j -α 3 θ C θ 1 ε k u i θ + α 3 θ C θ 1 C θ 1 ε k a ij u j θ -α 3 θ C θ 2 P U θi -α 3 θ C θ2 P T θi -α 3 θ C θ3 G i -(1 -α 3 θ )C ε ε k u i θ + u j θ n i n j -α 3 θ ε h θi + u i θ ε 2k + ε θ 2 2θ 2 = 0 , (32)
such that the algebraic model now reads, assuming equilibrium (Eqs. 27 and 28),

u i θ = ξ P U θi + ζ P T θi + ηG θi -γ ε k u j θ n i n j + χ ε k a ij u j θ -α 3 θ ε h iθ α 3 θ C θ 1 + (1 -α 3 θ )C ε ε k ( 33
)
where

ξ = (1 -α 3 θ C θ 2 ), ζ = (1 -α 3 θ C θ 2 ), χ = α 3 θ C θ 1 C θ1 , η = (1 -α 3 θ C θ 3 ), γ = (1 - α 3 θ ) (1 + C ε ), C ε = 1 2 (1 + 1/Pr).
Assuming that ε h iθ = 0, this relation reduces to

u i θ = -C θ k ε ζ u i u j ∂θ ∂ x j + ξ u j θ ∂u i ∂ x j + ηβg i θ 2 -χ ε k a ij u j θ + γ ε k u j θ n i n j , ( 34
)
where

C θ = C θ α 3 θ C θ 1 + (1 -α 3 θ )C ε . ( 35
)
This model only differs from Eq. 31 by the additional term γ ε k u j θ n i n j , which sensitizes the model to the orientation of the wall, and, above all, by the fact that the coefficients ξ , ζ , χ , η and C θ are now dependent on the blending function α θ . It is easy to see that in regions far from the wall, where α θ → 1, the influence of the elliptic blending method vanishes and the model (Eq. 31) is recovered.

In the following, the usual constants of the linear high-Reynolds number model are used: C θ1 = 3, C θ 2 = 0.55, C θ 2 = 0, C θ 3 = 0.55 and C θ1 = 0. Note that the elliptic blending method does not introduce any new coefficient.

Moreover, the model can be further simplified by using ζ = 1 and ξ = η = χ = 0, which yields the EB-GGDH model

u i θ = -C θ k ε u i u j ∂θ ∂ x j + γ ε k u j θ n i n j . ( 36
)

Equation for the temperature variance

The transport equation for the temperature variance θ 2 , involved in Eq. 34, is written as

∂θ 2 ∂t + ∂ u j θ 2 ∂ x j = ∂ ∂ x k ν Pr δ kl + C θθ ρu k u l T ∂θ 2 ∂ x l + P θθ -ε θθ ( 37
)
In this equation, the production term, equal to

P θθ = -2 u i θ ∂θ ∂ x i
does not require modeling, and the dissipation rate ε θθ is obtained via the ratio R of thermal to mechanical time scales

R = θ 2 ε θθ ε k .
For large turbulent Peclet numbers, which can be considered valid far from the wall, for fluids with not too small Prandtl numbers, assuming a constant thermal-tomechanical time-scale ratio R = R h , where the superscript h denotes homogeneous or quasi-homogeneous, is a usual practice [START_REF] Hanjalić | One-point closure models for buoyancy-driven turbulent flows[END_REF]. In the present case, where the Prandtl number is close to unity, the widely admitted value R h = 0.5 is used. However, in the near-wall region, this assumption fails and, in particular, R approaches the Prandtl number at the wall. In order to account for this limit, in the framework of the elliptic blending approach, R can be modeled as

R = α 3 θ R h + (1 -α 3 θ )Pr (38)
The DNS database of Abe and Kawamura [START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1,020 with Pr = 0.025 and 0.71[END_REF], in which the ratio R is available, enables the evaluation of the performance of this simple model in forced convection regime. The flow between two parallel flat plates is driven by an imposed pressure gradient (see Fig. 3), and a constant wall heat flux qw is imposed. Buoyancy terms are neglected, such that temperature is a passive scalar. The two parameters of the flow are the friction Reynolds number Re τ and the Prandtl number Pr = 0.71.

Figure 2 shows the time-scale ratio R obtained with Eq. 38, for several value of the length scale used in the equation of α 3 θ (see Section 4 for the discussion about this length scale). In this a priori procedure, in order to obtain α θ , Eq. 14 is first solved, with the mechanical length scale L evaluated using DNS data and L θ taken proportional to L. It can be seen that using the simple model (Eq. 38) yields a significant improvement of the prediction of R in the near-wall region, although the complex shape below y + = 50 is not reproduced. Therefore, Eq. 38 will be used in the computations presented in Section 5 for the EB-AFM and EB-GGDH models. It is found in Fig. 2 that the optimal value for L θ /L is 1.75. Next section is devoted to the investigation of the role of this length scale in the modeling of φ * θi -ε θi .

Influence of the Modeling of the Thermal Length Scale

As mentioned in Section 2.2, Shin et al. [START_REF] Shin | Elliptic relaxation second moment closure for the turbulent heat fluxes[END_REF] assumed that the blending function α θ used in Eq. 13 is equal to the blending function α used in Eq. 3, which is a strong assumption. Since α and α θ are obtained from similar elliptic Eqs. 4 and 14, this is equivalent to the assumption L θ = L. Actually, L θ originates from the approximation of the two-point correlations by exponential functions in Eq. 11, and the correlation length scale has no reason to be the same for temperature-pressure correlations as for velocity-pressure correlations. At large turbulent Peclet numbers Pe t = Re t Pr, the large-scale temperature fluctuations are simply convected, such that L θ can be considered entirely determined by the velocity field, and assuming a constant L θ /L is reasonable, but this does not imply that the constant is exactly one. Therefore, the present section focuses on the modeling of the thermal length scale L θ , and its influence on the predictions of φ * θi and ε θi , using a priori tests in forced and mixed convection regimes.

In order to cover a wide range of regimes, it is usual to assume that the time-scale used in heat fluxes models is a function k, ε, k θ and ε θ , i.e., a hybrid of mechanical and thermal time-scales (see, e.g., Abe et al. [START_REF] Abe | A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows-II. Thermal field calculations[END_REF]). This hypothesis can be extended to L θ , and a dimensional analysis shows that the ratio

L θ /L is then a function of R L θ L = f (R) , ( 39
)
where L = k 3/2 /ε. In order to account for the variations of R in the near-wall region in a simple manner, Eq. 39 can be approximated by a linear relation f (R)

L θ = C L R L . ( 40
)
Moreover, as mentioned above, assuming a constant thermal-to-mechanical timescale ratio R is usual for large Peclet numbers. This assumption implies a constant ratio L θ /L. In order to assess a priori the validity of the simple relation 40 and to evaluate the coefficient C L , two DNS databases of fully developed turbulent channel flows are used, in forced and mixed convection regimes. The a priori procedure consists in first solving Eq. 14, with the modeled L θ evaluated using DNS data, and then computing the model 13, also from DNS data. The standard coefficients [START_REF] Launder | On the computation of convective heat transfer in complex turbulent flows[END_REF] of the model for the scrambling term Eqs. 22 and 23 are used (C θ 1 = 3, C θ 2 = 0.55, C θ2 = 0, C θ 3 = 0.55, C θ1 = 0), which yields for the EB-AFM model (Eq. 34)

ξ = (1 -0.55 α 3 θ ) , ζ = 1 , χ = 0 , η = (1 -0.55 α 3 θ ) , γ = (1 -α 3 θ ) (1 + C ε ) , C ε = 1 2 1 + 1 Pr ,
and

C θ = C θ 3 α 3 θ + (1 -α 3 θ ) C ε ,
where C θ = 0.68 is used.

In addition to case of Abe and Kawamura [START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1,020 with Pr = 0.025 and 0.71[END_REF] described above and in Fig. 3, the DNS database of Kasagi and Nishimura [START_REF] Kasagi | Direct numerical simulation of combined forced and natural turbulent convection in a vertical plane channel[END_REF] for the mixed convection regime is used. Buoyancy is accounted for in a vertical channel flow driven by a pressure gradient, with an imposed wall-temperature difference (see Fig. 4). The Reynolds number based on the mean friction velocity is Re τ = 150. The Prandtl number is Pr = 0.71 and the Grashof number based on the temperature difference Gr = 9.6 × 10 5 .

Figures 5 and6 show φ * θ 1 -ε θ 1 (streamwise direction) and φ * θ 2 -ε θ 2 (wall-normal direction), respectively, obtained with different values of L θ . Figures 7 and8 show the same results for the mixed convection case. The difference of sensitivity to a modification of the length of these two components is striking. Contrary to φ * θ 1 -ε θ 1 , for which doubling the length scale yields a significant modification of the results, φ * θ2 -ε θ2 remains nearly unaffected. The reason for this behavior lies in the term γ ε k u k θ n i n k arising when the elliptic blending approach is introduced (see Eq. 34). Indeed, the normal-wise component writes Fig. 4 Channel flow in mixed convection regime, Kasagi and Nishimura [START_REF] Kasagi | Direct numerical simulation of combined forced and natural turbulent convection in a vertical plane channel[END_REF] which can be recast as

v θ = -C θ k ε v 2 ∂θ ∂ y + γ ε k v θ , ( 41
)
v θ = -C * θ k ε v 2 ∂θ ∂ y , ( 42
)
where

C * θ = C θ C θ 1 α 3 θ + (1 -α 3 θ ) C ε + (1 + C ε ) C θ C θ 1 -1
.

Far from the wall, C * θ → C θ /C θ 1 , and Eq. 42 tends to the form given by standard models (GGDH or AFM).

For Pr = 0.71, C ε ≈ 1.204, such that C * θ goes from 0.23 far from the wall to 0.25 at the wall. Therefore, whatever the shape of α θ , i.e., whatever the model used for L θ , v θ remains nearly the same. This remark leads to an important conclusion: the reason why standard models (GGDH and AFM) successfully reproduce the mean temperature profile in channel flows for fluids with a Prandtl number close to unity, despite the fact that they do not account for near-wall effects on the heat fluxes, is a coincidence. Indeed, in such flows, the streamwise heat flux is not active, and the wall-normal heat flux is correctly reproduced because, for such Prandlt numbers, C * θ is almost a constant. Of course, this feature relies on the correct prediction of the wall-normal component of the Reynolds stress, for which accounting for near-wall effects (which is the case with the EB-RSM) is crucial. However, the situation is different in more complex geometries, where a correct prediction of all the component of the turbulent heat fluxes is important, as well as for Prandtl numbers far from unity: Fig. 9 illustrates the distribution of the coefficient C * θ in the domain for different values of the Prandtl number. It is clearly seen that the variation of C * θ cannot be neglected close to the wall for large or small Prandtl numbers, and, consequently, accounting for wall effects in the model for the turbulent heat fluxes is necessary to correctly reproduce temperature profiles. It is worth pointing out that, for such values of the Prandtl number, other modeling issues arise, such as the modeling of ε h θi , which cannot be neglected for small Prandtl numbers, since the small-scale temperature fluctuations can be considered isotropic only at the limit of large Peclet numbers [START_REF] Kenjereš | Numerical simulation of magnetic control of heat transfer in thermal convection[END_REF][START_REF] Tennekes | A First course in Turbulence[END_REF].

On the contrary, for both cases, the streamwise component φ * θ 1 -ε θ 1 is significantly improved by increasing the ratio L θ /L (albeit not shown here, the improvement is particularly pronounced for φ * θ 1 ). Considering the present a priori test and the conclusion of Section 3.3, the optimal value is L θ /L = 1.75. In the forced convection case, since the time-scale ratio R is provided by the database, the linear hypothesis (Eq. 40) can be also evaluated with a variable R. It can be observed that the results obtained with L θ = 3 R L and L θ = 1.75 L are almost identical, such that using a constant ratio L θ /L is considered sufficient.

Computational Results

EDF in-house open source (http://www.code-saturne.org) CFD tools Code_Saturne is used for the present computations. Code_Saturne is an unstructured, collocated finite volume solver for cells of any shape. Spatial discretization is second order accurate. The Navier-Stokes equations are solved for turbulent incompressible flows using a SIMPLEC algorithm for pressure-velocity coupling (see Archambeau et al. [START_REF] Archambeau | Code Saturne: a finite volume code for the computation of turbulent incompressible flows-industrial applications[END_REF] for additional details about the code). It includes several RANS models with first and second moment closures such as the standard k -ε, k -ω-SST, stabilized v 2f (φmodel), LRR and SSG models. The code can also perform Large-Eddy Simulation using standard or dynamic Smagorinsky models. The elliptic blending approach has been introduced into the standard SSG model available in Code_Saturne in order to implement the EB-RSM, which is used for all the computations. For all the solutions presented in the present paper, grid convergence has been carefully ensured. The first computational point for both test cases are located around y + = 0.1.

The different turbulent heat fluxes models used in the present work, and summarized below, have been implemented:

-GGDH (Generalized Gradient Diffusion Hypothesis)

u i θ = -C θ k ε u i u j ∂θ ∂ x j with C θ = 0.235. -AFM (Algebraic Flux Model) u i θ = -C θ k ε u i u j ∂θ ∂ x j + ξ u j θ ∂u i ∂ x j + ηβg i θ 2
with C θ = 0.235, ξ = 0.45 and η = 0.45.

-EB-GGDH (Elliptic Blending-Generalized Gradient Diffusion Hypothesis)

u i θ = -C θ k ε u i u j ∂θ ∂ x j + γ ε k u k θ n i n k (43) with C θ = 0.68 3 α 3 θ + (1 -α 3 θ )C ε , γ = 1 -α 3 θ 1 + C ε and C ε = 1 2 1 + 1 Pr
-The EB-AFM (Elliptic Blending-Algebraic Flux Model)

u i θ = -C θ k ε u i u j ∂θ ∂ x j + ξ u j θ ∂ u i ∂ x j + ηβg i θ 2 + γ ε k u j θ n i n j
with ξ = 1 -0.55α 3 θ and η = 1 -0.55α 3 θ (these coefficients tend to the coefficients of the AFM far from the wall). The temperature variance θ 2 is obtain with the Eqs. 37 and 38.

It is worth pointing out that the coefficients used in the GGDH, AFM, EB-GGDH and EB-AFM all originate from the same calibration [START_REF] Launder | On the computation of convective heat transfer in complex turbulent flows[END_REF] of the underlying model for the scrambling term (Eqs. Computations were performed in order to evaluate the advantage of the introduction of the elliptic blending approach combined either with the GGDH or the AFM model on the prediction of the turbulent heat fluxes and mean velocity/temperature profiles, as well as to confirm the necessity of using a thermal length scale different from the mechanical length scale. The two DNS test-cases used above for a priori tests are investigated.

Forced convection regime In a fully developed turbulent channel flow, the EB-RSM model gives very satisfactory results for the mean velocity and the Reynolds stresses, as shown in Fig. 10 and in Fig. 11 , but this is of course independent of the turbulent heat fluxes model. Figure 12 shows the temperature profile for the different turbulent heat fluxes models. Note that, in forced convection, the AFM and EB-AFM results give exactly the same results as the GGDH or the EB-GGDH, respectively. As mentioned in Section 4, the accounting for the near-wall effects only weakly modifies the prediction of the mean temperature profile, since the wall-normal heat flux is only marginally affected by the elliptic blending approach. This is confirmed by Fig. 13, in θ which it can be seen that the predictions given by all the heat fluxes models are not distinguishable.

For the sake of clarity, the predictions for the streamwise heat flux by the GGDH and AFM models have been split in the two Figs. 14 and 15. These two figures clearly show the improvement in the near-wall region introduced by elliptic blending, in particular in the case of the AFM. Moreover, the benefit of using a thermal length scale different from the mechanical length scale is confirmed, although it does not constitute a major breakthrough. It is worth pointing out that differentiating these two length scales has a computational cost, since it makes necessary the resolution of an additional elliptic equation for α θ (otherwise α θ = α), but the cost of the resolution of this equation is very low (typically less than 1% of the total computational time).

Mixed convection For the mixed convection regime, the EB-RSM gives satisfactory results for the mean velocity (Fig. 16) and the Reynolds Tensor (Fig. 17) whatever the heat fluxes model. Contrary to the case of forced convection, it can be seen in Figs. 18 and 19 that the predictions of the wall-normal heat flux and, consequently, the mean temperature, are slightly influenced by the introduction of the near-wall effects. The most significant improvement is again for the streamwise component (Fig. 20), but also for the temperature variance, as seen in Fig. 21. Overall, the value L θ /L = 1.75 is an appropriate choice, which confirms the conclusion of Section 3.3.

Conclusions

The main contribution of the present paper is the proposal of a new algebraic model for the turbulent heat fluxes, the EB-AFM, as well as a simplified form, the EB-GGDH. These models are similar to usual AFM and GGDH, but reproduce the influence of the wall on the heat fluxes. They are algebraic versions of the EB-DFM [START_REF] Shin | Elliptic relaxation second moment closure for the turbulent heat fluxes[END_REF], based on transport equations for the turbulent heat fluxes in which the blocking effect of the wall on the scrambling term is reproduced by the elliptic blending approach (hence the EB prefix).

Moreover, the analytical justification of the use of the elliptic relaxation approach for the turbulent heat fluxes was introduced, which forms the basis for the EB-DFM. In addition, it was proposed to refine the modeling of the thermal-to-mechanical time-scale ratio R, used to evaluate the dissipation of the temperature variance, in order to impose the correct wall limit, as well as to distinguish the mechanical and thermal length scales in the elliptic blending approach for the Reynolds stresses and the heat fluxes, respectively.

The validation of these proposals was carried out using both a priori tests and numerical simulations in forced and mixed convection regimes. The predictions of the streamwise turbulent heat flux and the temperature variance are significantly improved when the near-wall effects are accounted for by elliptic blending. Using a thermal length scale larger than the mechanical length scale proved also beneficial, although the improvement is moderate. On the contrary, in channel flow cases, for fluids with Prandtl numbers close to unity, the wall-normal turbulent heat flux is only moderately affected by the introduction of the effects of the wall, which provides an explanation for the success of standard models (AFM and GGDH) in channel flow configurations, although they derive from a model not valid in the near-wall region. This is to be traced to the fact that, fortuitously, the two limiting values (at the wall and far from the wall) of the variable coefficient driving this component of the heat fluxes are almost equal for Pr = 0.71. Accounting for the near-wall effects might have an important influence in cases where the Prandtl number is far from unity.

Fig. 1

 1 Fig.1Influence of buoyancy on α evaluated a priori by Eq. 6. DNS data from Kasagi and Nishimura[START_REF] Kasagi | Direct numerical simulation of combined forced and natural turbulent convection in a vertical plane channel[END_REF]

Fig. 2

 2 Fig.2Thermal-to-mechanical time-scale ratio R obtained from Eq. 38. A priori test in forced convection, Re τ = 640)[START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1,020 with Pr = 0.025 and 0.71[END_REF]

Fig. 3

 3 Fig.3Channel flow in forced convection regime, Abe and Kawamura[START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1,020 with Pr = 0.025 and 0.71[END_REF]

Fig. 5 AFig. 6 2 -Fig. 7 A 5 DNSFig. 8 * θ 2 -ε θ 2 0with

 56275822 Fig.5A priori test of the model for φ * θ 1 -ε θ 1 in forced convection regime[START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1,020 with Pr = 0.025 and 0.71[END_REF] (Re τ = 640)

  22 and 23):C θ 1 = 3, C θ 2 = 0.55, C θ2 = 0, C θ 3 = 0.55, C θ1 = 0. The recalibration coefficient (see the discussion in Section 3.1) C θ , involved in C θ = C θ /C θ1is chosen for each model in order to optimize the reproduction of the mean temperature profile in the forced convection case. The corresponding values are C θ = 0.705 for GGDH and AFM, and C θ = 0.68 for EB-GGDH and EB-AFM.

Fig. 10

 10 Fig.[START_REF] Fadai-Ghotbi | Temporal filtering: a consistent formalism for seamless hybrid RANS-LES modeling in inhomogeneous turbulence[END_REF] The velocity profile (computational results for the forced convection regime, Re τ = 640)

Fig. 11 Fig. 12 Fig. 13 Fig. 14 Fig. 15 Fig. 16 Fig. 18 Fig. 19 Fig. 20 Fig. 21 Temperature variance θ 2 (

 1112131415161819202 Fig. 11 The Reynolds tensor profile. (Symbols: DNS of Kawamura Re τ = 640 [1]) (computational results for the forced convection regime)

Table 1

 1 Asymptotic behaviors of the terms of the Reynolds stress budget in buoyancy affected flows


Claudio Vergari
email: c.vergari@gmail.com

Laurent Gajny

Isabelle Courtois

Eric Ebermeyer

Kariman Abelin-Genevois

Youngwoo Kim

Tristan Langlais

Raphaël Vialle

Ayman Assi

Ismat Ghanem

Eric Ebemeyer

Raphael Vialle

Jean Dubousset

Wafa Skalli

Quasi-automatic early detection of progressive idiopathic scoliosis from biplanar radiography: a preliminary validation

Keywords: adolescent idiopathic scoliosis, 3D reconstruction, reliability, feature extraction, severity index

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Adolescent idiopathic scoliosis (AIS) is a three-dimensional deformity of the spine which, if left untreated in its progressive form, can lead to respiratory or cardiovascular impairment and more in general a decrease in quality of life [START_REF] Negrini | SOSORT guidelines: orthopaedic and rehabilitation treatment of idiopathic scoliosis during growth[END_REF]. Corrective surgery is invasive, and it is accompanied by non-negligible risks of complications. Bracing can be used to efficiently slow down or stop the progression, and patients often prefer it to observation [START_REF] Weinstein | Effects of Bracing in Adolescents with Idiopathic Scoliosis[END_REF][START_REF] Negrini | Braces for idiopathic scoliosis in adolescents[END_REF]. Early detection of progressive curves offers the opportunity to undergo effective, non-operative treatments to the patients, thus reducing the need of surgery. This concept lead several associations to endorse screening for scoliosis [START_REF] Roussouly | Pre-and post-operative sagittal balance in idiopathic scoliosis: a comparison over the ages of two cohorts of 132 adolescents and 52 adults[END_REF][START_REF] Grivas | The pendulum swings back to scoliosis screening: screening policies for early detection and treatment of idiopathic scoliosiscurrent concepts and recommendations[END_REF], in order to detect the pathology as early as possible. Nevertheless, such screening remains controversial [START_REF] Hresko | Early Detection of Idiopathic Scoliosis in Adolescents[END_REF], mainly because of cost and potentially unnecessary radiation exposure, although today low-dose radiography is available. Moreover, the benefits of detection of early scoliosis cases are hindered by the lack of robust methods to recognized which of these cases are at risk of progression.

Several factors are associated with curve progression, the most common being the skeletal maturity and menarcheal age, Cobb angle, and scoliosis in the family. Other geometrical factors could be associated with progression, such as the rib-vertebra angle as proposed by Mehta [START_REF] Mehta | The rib-vertebra angle in the early diagnosis between resolving and progressive infantile scoliosis[END_REF], the back surface shape [START_REF] Korvin G De | Detection of progressive idiopathic scoliosis during growth using back surface topography: A prospective study of 100 patients[END_REF], spinal torsion [START_REF] Poncet | Geometric torsion in idiopathic scoliosis: threedimensional analysis and proposal for a new classification[END_REF], as well as patient characteristics such as alterations of the vestibular system or the spinal cord [START_REF] Kong | Variation in Anisotropy and Diffusivity along the Medulla Oblongata and the Whole Spinal Cord in Adolescent Idiopathic Scoliosis: A Pilot Study Using Diffusion Tensor Imaging[END_REF]. In particular, it was suggested that the threedimensional shape of the scoliotic curve could be indicative of progression risk [START_REF] Perdriolle | A study of scoliotic curve. The importance of extension and vertebral rotation (author's transl)[END_REF][START_REF] Dubousset | Three-dimensional Analysis of the Scoliotic Deformity[END_REF]. Indeed, a 3D phenotype of progressive curves was recently defined that can discriminate progressive and stable scoliosis with an overall accuracy of 82%, and 89% of correctly detected progressive curves [START_REF] Skalli | Early detection of progressive adolescent idiopathic scoliosis: A severity index[END_REF]. However, a limitation of the previous study was that a fully experienced user was needed to perform a thorough 3D reconstruction of the spine from biplanar radiography, which might be difficult to implement in clinical routine.

The aim of this work was to develop a fast and semi-automatic method to calculate a robust severity index (s-index), to predict progression in untreated IS during growth.

Methods

Subjects

Data were collected prospectively from four hospitals, after approval of the ethical committee. Inclusion criteria were: diagnosis of AIS, Cobb angle between 10° and 25°, European Risser sign < 3, age > 10 years, no previous treatment for scoliosis. All curve topologies were included (thoracic, thoracolumbar and lumbar). Exclusion criteria were: supernumerary vertebrae, transition anomaly, non-standard position in biplanar radiography. All patients underwent a low-dose biplanar radiography (EOS system, EOS Imaging, Paris, France) at inclusion, which corresponded to their diagnosis of AIS, in the standardized free-standing position [START_REF] Faro | Evaluation of a functional position for lateral radiograph acquisition in adolescent idiopathic scoliosis[END_REF]. Patients who did not respect this position were excluded (arms held too high, too low or asymmetrically, obvious leaning forwards or backwards, etc.).

Part of this cohort (thirty patients) was already used in the previous validation work [START_REF] Skalli | Early detection of progressive adolescent idiopathic scoliosis: A severity index[END_REF], but some patients were imaged with older stereoradiographic devices with low image quality; these patients were excluded from the current study. Cohort size for this sensitivity/specificity study was defined based on the hypothesis of 40 % prevalence for progressive AIS [START_REF] Danielsson | A prospective study of brace treatment versus observation alone in adolescent idiopathic scoliosis: A follow-up mean of 16 years after maturity[END_REF][START_REF] Bunnell | The Natural History of Idiopathic Scoliosis Before Skeletal Maturity[END_REF][START_REF] Felice | The Natural History of Idiopathic Scoliosis during Growth: A Meta-Analysis[END_REF] and a statistical power of 0.8 [START_REF] Bujang | Requirements for Minimum Sample Size for Sensitivity and Specificity Analysis[END_REF], leading to a minimal cohort of 50 patients with at least 20 progressive cases.

Patients were followed until prescription of brace (considered as progressive scoliosis) or skeletal maturity without bracing (Risser sign ≥ 3 and Cobb angle < 25°, considered as stable scoliosis). Follow-up was also performed with EOS biplanar x-rays. Patients were considered progressive when showing Risser sign ≤ 2 and a Cobb angle higher than 25°, consistently with SOSORT guidelines [START_REF] Negrini | SOSORT guidelines: orthopaedic and rehabilitation treatment of idiopathic scoliosis during growth[END_REF], but also if they showed an increase of more than 5° Cobb angle and/or vertebral axial rotation in 6 months [START_REF] Lonstein | The prediction of curve progression in untreated idiopathic scoliosis during growth[END_REF], together with worsening of the clinical appearance with hump increase. The progressive character of each patient was confirmed by consensus of the prescribing clinician and an expert surgeon (J.D.), who analyzed the biplanar radiography of these patients at decision of treatment while blinded to the value of the severity index. Patients for whom there was no consensus were excluded.

Quasi-automatic 3D reconstruction

A quasi-automatic 3D reconstruction of the spine was performed in three steps, which will be detailed below: the operator first selected a few anatomical landmarks on the frontal and lateral radiographs (Figure 1A), then an automatic algorithm provided an initial solution of 3D reconstruction, on which the operator could perform fast manual adjustment of some key vertebrae (the apex, junctions and the adjacent to the junctions, Figure 1B) which were projected over the original radiography (Figure 2). More in detail, the following landmarks were digitized in both views (Figure 1): the sacral plate, left and right acetabula and the spinal midline through the centre of all vertebral bodies, from the odontoid apophysis of C2 to L5 lower endplate. The upper endplate of C7 and lower endplate of T12 were also selected in the lateral view, while the two endplates delimiting the main scoliotic curve were selected in the frontal view. The junctional levels of the main curve Fig. 1 (A) Manual initialization of the 3D reconstruction. The sacral plate, acetabula and the spinal midline from the odontoid apophysis of C2 vertebra to L5 lower endplate were selected in both views. The upper endplate of C7 and lower endplate of T12 were selected in the lateral view, while the main scoliotic curve was selected in the frontal view. (B) 3D reconstruction of the 7 key vertebrae that were manually adjusted: the apex, the junctional vertebrae and the adjacent to the junctional vertebrae.

were defined to obtain the maximization of Cobb angle, and they were often characterized by a local discontinuity of vertebral axial rotation (i.e., a sudden change of vertebral axial orientation).

The initial solution was automatically calculated through statistical inference [START_REF] Humbert | 3D reconstruction of the spine from biplanar X-rays using parametric models based on transversal and longitudinal inferences[END_REF][START_REF] Gajny | Quasi-automatic 3D reconstruction of the full spine from low-dose biplanar X-rays based on statistical inferences and image analysis[END_REF] and image analysis [START_REF] Ebrahimi | Lumbar spine posterior corner detection in X-rays using Haar-based features[END_REF][START_REF] Ebrahimi | Vertebral Corners Detection on Sagittal X-rays based on Shape Modelling, Random Forest Classifiers and Dedicated Visual Features[END_REF]. The 3D reconstructed vertebrae thus obtained were retro-projected (Figure 2) on the frontal radiography. The operator was only required to adjust horizontal and vertical position and axial rotation of 7 key vertebrae: the curve apex, two junctional vertebrae and the adjacent to the junctional vertebrae. The aim of this step was to better fit the projection of the 3D model to the frontal radiography; the sagittal radiography was not used at this stage since the initial solution was sufficiently accurate in this plane.

Severity index

The patient's severity index was automatically calculated, as previously described [START_REF] Skalli | Early detection of progressive adolescent idiopathic scoliosis: A severity index[END_REF], from the following parameters relative to the main scoliotic curve: Cobb angle, vertebral axial rotation at the apical level, intervertebral axial rotation at the upper and lower junctional levels, torsion index of the curve and hypokyphosis index of the apical vertebra. All parameters were calculated in the patient's reference frame. Predictive discriminant analysis was used to compare each patient's parameters to a previously built database of 53 healthy subjects, 17 prebracing and 31 severe AIS patients. The comparison was summarized in a severity index: an index lower than 0.4 was associated with a stable scoliosis while an index higher than 0.6 with a progressive one. Values in between were considered doubtful (non-classified).

Validation

The relevance of the severity index calculated with this quasi-automatic method was evaluated by comparing the index's classification with the actual patients' outcome (progressive or stable). A confusion matrix was built with the number of true negatives, false negatives, true positives and false negatives. Sensitivity (true positives rate) and specificity (true negatives rate) were also calculated.

A reproducibility study was performed to determine reliability of the severity index: three operators (two researchers with two and five-years' experience in 3D reconstruction, Fig. 2 Position of the vertebral model for apical vertebra (L2) before (panel A) and after manual adjustment. In panel B, the contour does not perfectly match the vertebra because only fast rigid transformations (translation and axial rotation) were required to compute the severity index. In this example, the vertebra was displaced by 3 mm (corresponding to 16 pixels) and it was axially rotated by less than 1°.

respectively, and one surgeon with one day experience) repeated the quasi-automatic reconstruction process twice on a subset of 30 patients. Cohen's kappa and percentage agreement were calculated to assess intra-operator repeatability while Intraclass Correlation Coefficient (ICC) was calculated to determine inter-operator reproducibility. All calculations were performed in Matlab 2016b.

Results

Patients

Fifty-five patients were included, forty-six girls (12 ± 1 years old, European Risser sign ≤ 2, 16 ± 4° Cobb angle) and nine boys (13 ± 1 years old, Risser sign ≤ 2,16 ± 3° Cobb angle). Table 1 reports patient demographical data and their geometrical parameters related to severity index. No patient dropped out of the study. Thirty-six patients were progressive while 19 were stable.

Severity index

Table 2 reports the confusion matrix of the severity index classification based on the fast 3D reconstruction of the spine. Overall, 85% of the patients were correctly classified, with a sensitivity of 92% and specificity of 74%. Three patients (8 %) were false negative, i.e., their progressive character was not recognized, but they all had a Cobb angle lower than 12°. Four patients were false positive (stable curves detected as progressive).

Table 1. Clinical parameters of the included cohort (average ± standard deviation). Curve topologies: T = Thoracic (apex in T1-T11), TL = thoraco-lumbar (apex in T12-L1), Lumbar (apex in L2-L5). Age, and European Risser sign were expressed as median [1st-3rd quartile]. Reliability 180 reconstructions were performed (30 patients x 3 operators x 2 reconstructions). "Substantial" intra-operator agreement [START_REF] Landis | The Measurement of Observer Agreement for Categorical Data[END_REF] was observed for all three operators (Cohen's kappa between 0.6 and 0.7), with percentage agreement between 81 and 85%, while ICC showed good inter-operator agreement (0.63). Fourteen patients in this cohort were progressive, corresponding to 84 reconstructions (14 x 3 operators x 2 repetitions): 80% of these reconstructions were correctly classified as progressive. Figure 3 reports a receiver operating characteristic curve for the detection of progressive scoliosis.

Cohort

Quasi-automatic 3D reconstruction

The initial manual identification of the spinal curve and specific endplates lasted 2.5 minutes in average, while the final manual adjustment of seven key vertebrae lasted less than one minute. Fig. 2 shows an example of rigid manual adjustment.

Overall, the apexes were displaced by 0.1 ± 0.2 mm (max 4.2 mm) and rotated by 0.1 ± 0.1° (max 1.8°).

Discussion

A semi-automatic method to calculate a predictive index of progressive scoliosis from biplanar radiography was described in the present work. The performance of the s-index confirmed the results of the previous validation, which was based on manual 3D reconstruction: at least 80% of the progressive patients were detected in both studies, even when inter-operator variability was introduced. However, operator-time time was reduced from about 10-15 minutes for an expert operator to less than 4 minutes for a newly trained one. Operator-time did not depend on the operator' experience; however, all operators had knowledge of spinal anatomy and of the radiological appearance of the important landmarks. Moreover, this implementation of the severity index showed good to substantial intra-and inter-operator reliability.

Three patients were false negatives, i.e. they were progressive, but they were classified as stable scoliosis. However, these patients were at a very early stage, with Cobb angle lower than 12°. Therefore, a low severity index for such an early stage of scoliosis should not be interpreted as a zero risk of progression. Rather, the patient is not at an immediate risk of progression, and the next follow-up exam and s-index evaluation could be delayed. Four false positives were stable scoliosis which were classified as progressive. These patients had Cobb angle between 16 and 25°, and they all presented relatively high apical axial rotation (between 4° and 8°). Further refinement of the s-index through larger databases and additional clinical parameters could improve these rates.

Rate of scoliosis progression was 65 % in this study, since 36 progressive patients out of 55. A recent meta-analysis by Di Felice et al. [START_REF] Felice | The Natural History of Idiopathic Scoliosis during Growth: A Meta-Analysis[END_REF] reported a progression rate of 49%, with a large 95% confidence interval of 19%-79%. The present study falls within the higher spectrum of this confidence interval, while the study by Danielsson et al [START_REF] Danielsson | A prospective study of brace treatment versus observation alone in adolescent idiopathic scoliosis: A follow-up mean of 16 years after maturity[END_REF] falls on the lower part of the spectrum, with a 40% progression rate. Lonstein and Carlson [START_REF] Lonstein | The prediction of curve progression in untreated idiopathic scoliosis during growth[END_REF] reported an even lower rate of 23 % progression. However, progression rates should be compared with care, because of differences of inclusion and progression criteria. For instance, Danielsson et al. included patients with 23-39° Cobb angle (some of which would be considered as already progressed in this study), while Lonstein and Carlson included 12% of patients with Cobb angles lower than 10°, which would not be considered as scoliotic according to SRS and SOSORT criteria.

This work has some limitations. First, fifty-five patients still represent a relatively small cohort of patients; data collection is slow because of the long follow-up required to ascertain patient's outcome. Part of the cohort of patients included in this study was already published in the first validation of the severity index [START_REF] Skalli | Early detection of progressive adolescent idiopathic scoliosis: A severity index[END_REF]. In the present work there are actually less total patients because the previous cohort included several radiographies acquired using stereoradiographic devices [START_REF] Pomero | Fast accurate stereoradiographic 3Dreconstruction of the spine using a combined geometric and statistic model[END_REF] with lower image quality. These were not adapted to the image analysis techniques used for the initialization of the 3D reconstruction in the present method, where all acquisitions were performed with an EOS system.

Indeed, the proposed severity index requires the availability of calibrated biplanar radiographies of the patient, which are not available in every clinical center and could decrease the diffusion of the method. Moreover, inclusions were limited to AIS: non-idiopathic scoliosis were excluded, as well as early-onset scoliosis. Nevertheless, this study confirmed the previously obtained index's discriminatory power on adolescent patients [START_REF] Skalli | Early detection of progressive adolescent idiopathic scoliosis: A severity index[END_REF]. Moreover, this study allowed to validate a novel user-friendly method to calculate the severity index, based on the digitization of well-known anatomical landmarks (vertebral endplates, pelvic acetabula) and on simple horizontal translation and axial rotation of only 7 vertebrae. A multicentric international data collection is under way to include more patients at first exam and follow them to skeletal maturity or brace treatment, in order to form a larger cohort for extensive validation.

The second limitation is that criteria for brace decision vary between clinicians and hospitals, and it is difficult to enforce quantitative standards such as the SRS recommendations. The main criteria for brace decision were "an increase of 5° Cobb angle and/or vertebral axial rotation in 6 months, together with worsening of the clinical appearance with rib hump increase". However, the "worsening of the clinical appearance" could be subjective, and the complete clinical exam was taken into account by the clinician when deciding for brace. To overcome this limitation, radiographies of all progressive patient were reviewed by an expert senior surgeon, who was blinded to the value of the severity index, in order to establish collegiate clinical agreement on the progressive character of those patients. No more than a couple of patients were excluded from the study because of disagreement between the expert and the prescribing clinician, i.e., patients that the expert would not have treated by brace; these cases were not further documented. This collegiate review took place before the calculation of the severity index, so it is completely independent from it, allowing to confirm the progressive character of all included patients, as was done in our previous work [START_REF] Skalli | Early detection of progressive adolescent idiopathic scoliosis: A severity index[END_REF].

Finally, patients' follow-up was not complete for this study, since it stopped at the decision of bracing (progressive patients) or at European Risser sign ≥ 3 with Cobb angle < 25° (stable patients). On the one hand, progression after bracing is largely influenced by the patient's compliance and response to bracing, so these aspects become important factors that should be included in the estimation of progression risk towards need for surgery. On the other hand, rate of progression decreases at Risser 3 for patients presenting low Cobb angle, so estimation of risk is less crucial. Still, progression risk is not zero event at this stage, so patients can continue their radiological examinations, albeit less frequently [START_REF] Kohashi | A New Method Using Top Views of the Spine to Predict the Progression of Curves in Idiopathic Scoliosis During Growth[END_REF].

The proposed method to calculate the severity index requires manual intervention to initialize and finely adjust a 3D reconstruction of the main scoliotic curve. Although the whole process takes less than 4 minutes of operator-time, work is in progress towards a fully automatic 3D reconstruction, based on larger database of patients and image processing and analysis, to obtain a more accurate and robust calculation of the severity index.

In the present work, a semi-automatic method was proposed to calculate the severity index from biplanar radiography. Its main advantage is that it is fast and user-friendly. The resulting severity index was able to correctly detect 80% of the progressive patients. Once extensively validated, such index could allow very early initiation of conservative treatment for progressive patients, thus increasing treatment efficacy and therefore reducing the need for surgery.

Fig. 3

 3 Fig. 3 Receiver operating characteristic curve for the detection of progressive scoliosis. Area under curve (AUC) was 0.75.

Table 2 .

 2 Confusion matrix of progressive vs stable scoliosis detection.

		S-index ≤ 0.4	0.4 < S-index < 0.6 S-index ≥ 0.6
	Stable scoliosis (N = 19)	14 (74 %)	1	4 (21 %)
	Progressive scoliosis (N = 36)	3 (8 %)	0	33 (92 %)

Acknowledgements

The authors are grateful to the ParisTech BiomecAM chair program on subject-specific musculoskeletal modelling (with the support of ParisTech and Yves Cotrel Foundations, Société Générale, Proteor and Covea).