High-order conservative formulation of viscous terms for variable viscosity flows

The work presents a general strategy to design high-order conservative co-located finite-difference approximations of viscous/diffusion terms for flows featuring extreme variations of diffusive properties. The proposed scheme becomes equivalent to central finite-difference derivatives with corresponding order in the case of uniform flow properties, while in variable viscosity/diffusion conditions it grants a strong preservation and a proper telescoping of viscous/diffusion terms. Presented tests show that standard co-located discretisation of the viscous terms is not able to describe the flow when the viscosity field experiences substantial variations, while the proposed method always reproduces the correct behaviour. Thus, the process is recommended for such flows whose viscosity field highly varies, in both laminar and turbulent conditions, relying on a more robust approximation of diffuse terms in any situation. Hence, the proposed discretisation should be used in all these cases and, for example, in large eddy simulations of turbulent wall flows where the eddy viscosity abruptly changes in the near-wall region.

accounts for the cell-bound fluxes. In the field of computational fluid dynamics (CFD), the problem is not particularly felt when direct numerical simulations (DNS) are dealt with, since the resolution often mitigates it. Conversely, it is dominant in the case of LES. In particular, the LES approach, employing the same high-order discretisation techniques of the DNS, wants the computational mesh to be sufficiently coarse compared to the DNS setup. The method represents a superior technology if compared to a more standard Reynolds averaged Navier-Stokes (RANS) approach since the flow-dependent/energy-dominant large eddies are resolved rather than being modelled [15,16]. On the other hand, coarsening the grid and avoiding to inject artificial viscosity impacts on the correct behaviour of the local diffusion and quickly accumulates aliasing errors which result in the blow-up of the simulation. Since the majority of the stability issues are linked to the Navier-Stokes equations' nonlinear terms, many contributions have been proposed to better treat these components. In the last decades, the formulation of high-order shock-capturing numerical methods such the Essentially Non-Oscillatory (ENO) and the Weighted-Essentially Non-Oscillatory (WENO) schemes [17][18][19][20][21][22][23][24], today allows to evolve complex interactions between shocks and discontinuities stably. The extension of these schemes to a hybrid-compact/central formulations [25][26][27][28] combined with shock-detection strategies [29,30] improved the numerical resolution and reduced the computational efforts in the discretisation process. In particular, avoiding any artificial diffusion and dealing with high-order central-FD discretisation, the idea of splitting the convective derivatives and replicate the energy-conservation properties of the Euler equations represents one of the most valuable attempts in the field of accurate numerical simulation of compressible flows [28]. In this path, from the pioneering work of [31], many results have been proposed in the literature [32][33][34][35], involving more and more properties to the numerical discretisation. Indeed, the split formulations of the convective terms proposed by [27] represent the state of art in this field. Even if a lot of effort has been spent accounting for a more accurate and stable description of the Navier-Stokes nonlinear contributions, only a few words have been written concerning a better way to treat the viscous terms especially in a high-order FD framework and, in this field, no more than standard co-located approximations are employed.
In the present work, a high-order conservative FD discretisation of viscous fluxes for variable viscosity flows is presented. The method is proposed for compressible flow simulations but results in a generally conservative strategy in any flow conditions. In particular, the authors present a numerical strategy which conservatively discretised the incompressible contribution of the viscous terms, while the rest of the diffusive contributions, which represents a high-order series expansion of incompressible phenomena, still be computed in the co-located nodes. The discretisation process grants the proper telescoping of the viscous terms and strong conservation properties of the stress components in the limit of incompressible or weakly compressible flows, relying on a more robust approximation, especially in the under-resolved flow regions. The proposed procedure is then coupled with the energy-preserving scheme developed by Kennedy, Grüber, and Pirozzoli [27] in the treatment of the nonlinear convective terms. Thus, the numerical recipe results in a thoroughly conservative approximation of the Navier-Stokes system of equations in the limit of incompressible or weakly compressible flows. Since the variability of eddy-viscosity is inherent in the field of LES of wall-bounded flows, the method is a promising and a superior strategy for such kind of computations. Moreover, even if the conservation properties of the proposed method are recovered just for weakly compressible flows, a high-order conservative discretisation of viscous terms is of considerable help in any flow conditions, also in multiphase or stratified flows or in shocked states, if a high-order discretisation framework wants to be taken into account.
The present paper is organised as follows: In Sect. 2, the mathematical model is described, while, in Sect. 3, the numerical method employed in the discretisation of the Navier-Stokes equations is presented. In Sect. 4, the method is applied to both laminar and turbulent shockless flow conditions, providing the accuracy of the proposal for a wide variety of low Mach number test cases. Finally, in Sect. 5 concluding remarks are given.

Governing equations
The present study is carried out using URANOS (Unsteady Robust All-around Navier-StOkes Solver) [36,37], a general-purpose DNS/LES solver recently developed at the University of Padova. The solver deals with the compressible Navier-Stokes system of equations in a conservative formulation. In particular, the filtered version of the equations is employed. Being φ a general scalar quantity, both the Reynolds (φ φ + φ ) and the Favre (φ φ + φ ,φ ρφ/ρ) decompositions are used for the model description. Thus, the obtained Navier-Stokes system is: whereρ is the filtered density,ũ i is the Favre velocity component in the i-th direction,p is the filtered thermodynamic pressure,Ẽ ẽ + u i u i /2 is the Favre total energy per unit mass,ẽ is the Favre internal energy per unit mass, andJ j is the j-th component of the filtered molecular heat flux. The thermodynamic variables are correlated by the filtered pressure field expressed asp ρ RT whereT denotes the Favre temperature and R is the specific gas constant. The molecular viscous stress tensorτ i j and the j-th component of the filtered molecular heat flux,J j , are expressed as: whereS i j 1 2 g i j + g ji is the resolved strain-rate tensor, g i j ∂ũ i /∂ x j is the resolved velocity gradient, and μ(T ) is the molecular viscosity. Apart where different stated, the latter is assumed to obey the Sutherland's law. λ(T ) denotes the molecular thermal diffusivity, expressed as λ(T ) γ R/(γ − 1)μ(T )/Pr, while γ c p /c v and Pr are the specific heat ratio and the molecular Prandtl number, respectively. The latter parameters are set equal to 1.4 and 0.71. f j denotes the external forcing term components whose corresponding power (i.e. f jũ j ) is also included in the right-hand side (RHS) of the total energy equation.
Concerning the SGS stress tensor, T SG S i j ρu i u j −ρũ iũ j , the latter is modelled via the canonical Boussinesq's hypothesis. Thus, its deviator part is expressed as: where μ SG S is the SGS eddy-viscosity and T SG S kk the isotropic contribution. In the present model, T SG S kk is set to zero since it is expected to assume large values near shocks and discontinues (see for more details [16]). Finally, the SGS term of the filtered energy equation is expanded as: The SGS heat flux Q j is modelled as Q j −λ SG S ∂T /∂ x j where λ SG S μ SG S γ R/(γ − 1)/Pr SG S , while ψ ii j , the triple velocity correlation tensor, is neglected. The SGS Prandtl number Pr SG S is set equal to 0.9. The model is made fully non-dimensional since the groups γ, Re ∞ , Ma ∞ , and Pr univocally determine the flow conditions. The proposed version of the filtered Navier-Stokes equations relies on an efficient formulation since it allows to switch between a DNS approach to a LES model easily. In particular, if the grid resolution lies in the range of a DNS application, the SGS viscosity is set to zero, and the system recovers the Navier-Stokes unfiltered version completely. Otherwise, when the grid resolution is expected to induce some SGS turbulence effects, a specification for μ SG S is demanding, and the filtered version of the equations is automatically employed.
Concerning the SGS viscosity, the Wall-Adaptive Large-Eddy (WALE) viscosity model by [38] is employed. Among the variety of SGS models, the latter has peculiar characteristics in the field of wall turbulence since the SGS viscosity automatically vanishes between the bulk of the flow and the boundary layer without prescribing any damping function or artificial transition. In particular, the WALE model defines the SGS viscosity as where is the traceless symmetric part of the square of the resolved velocity gradient tensor. [38] proved that the model automatically detects both the laminar and the near-wall regions of the flow and provides the correct μ SG S /μ ∼ O(y + ) 3 asymptotic behaviour at the wall locations. C W 0.325 accurately fits a wide variety of flow conditions, while¯ ( x y z) 1/3 is the local mesh size.

Spatial discretisation of viscous terms
The entire model is discretised via a high-order central-FD approach. In particular, the solver allows both uniform and non-uniform structured Cartesian grids discretisation. As said, the peculiarity of the discretisation process stands in the numerical treatment of the viscous terms, which provides strong preservation of viscous stresses in the limit of incompressible or weakly compressible flows. In the following, the authors provide a detailed description of the numerical scheme related to the viscous contributions. Thanks to a simple manipulation of the viscous momentum fluxes of the Navier-Stokes equations, a formulation can be obtained in which the incompressible and the compressible contributions are highlighted.
In particular, the former acts in both compressible and incompressible flow conditions, while the latter is expected to take large values just near shocks and discontinuities or in highly expanding/deleting regions. μ tot μ + μ SG S denotes the overall viscosity. A numerical discretisation which accounts for the variability of the overall viscosity usually takes advantage of the standard Laplacian formulation: The method expands all the mixed derivatives and accounts for the viscosity and the velocity gradients separately in a co-located fashion way (see, e.g. [28,39]). Here α (1) l and α (2) l maximise the formal order of accuracy of a central approximation of 2L-size stencils for the first and the second derivative, respectively. Supposing the grid resolution in the range of the scheme convergence and the overall viscosity vary smoothly, the technique allows to resolve the viscous stresses accurately. However, especially in the case of poorly resolved flows or where the total viscosity experiences highly local variations, the method can provide erroneous behaviours of the diffusive contributions resulting in a wrong prediction of the entire flow field. In these situations, a numerical approximation which strongly conserves the shear stress components represents a more suitable description of incompressible viscous contributions. Let us describe the method for uniform Cartesian grids. The extension to non-uniform meshes is straightforward, and the interested reader is addressed to look at Appendix 6 for the details and the implementation issues. Thus, recalling the definition of a conservative method [27], the incompressible viscous stress components can be cast as: Hereτ i+1/2 is a high-order representation of the viscous stresses at the cell interface. The latter can be expressed as a linear combination which accounts for both the viscosity and the velocity gradient interpolations at the i + 1/2 node in a high-order path:τ The interpolation coefficients {β l , γ l } m l n depend on the desired order of accuracy, and the authors derived their values as listed in Table 1. In particular, the {β l } m l n coefficients provide an (m − n + 1)-order explicit interpolation for μ i+1/2 at the cell bound, while the {γ l } m l n allow for the interpolation of the velocity gradient in the intermediate node maximising the formal order of accuracy of a standard 2L-size stencils second derivative approximation, therefore respecting the following constraint: Here it is not worthless to be mentioned that, apart from 2nd-order of accuracy, the proposed {γ l } m l n coefficients do not represent the only possible choice for the interpolation of cell-bound gradients. In particular, according to [40], these coefficients do not maximise the formal order of accuracy of a first cell-bound derivative (d f/dx) i+1/2 . However, if the (d f/dx) i+1/2 is interpolated according to [40], the process leads to a cell-centred second derivative (d 2 f /dx 2 ) i with 2nd-order of accuracy if this is resembled in a conservative formulation, according to Eq. (10). Thus, the proposed {γ l } m l n coefficients are the only which provide a (d f/dx) i+1/2 interpolation able to fall back to standard co-located FD second derivatives, with corresponding order of accuracy, if uniform viscosity/diffusivity distributions are taken into account. In the following, the method is addressed as conservative discretisation. As said, the method recovers the standard Laplacian formulation (i.e. Eq. (9.1)), with corresponding order of accuracy, in the case of uniform viscosity fields. Besides, due to the inherent characteristic of conservative schemes, in the limit of incompressible or weakly compressible flows, the conservative approach provides the high-order conservation of the viscous terms independently to the resolution and grid stretching (see Appendix 6 for details). Thus, the method relies on a more robust strategy of diffusive terms in any situation, even if the flow embedded shocks.
To conclude, a note is demanded concerning the heat flux components in the total energy equation. Even if the conservative method is explained looking at the momentum viscous components, nothing forbids in applying it for the discrete treatment of the heat flux divergence observing that is structurally similar to the one given in Eq. (10). If employed, the process relies in the conservation also of the heat flux components. Hereλ tot λ + λ SG S is the total diffusion. In the present work, apart from where differently stated, the viscous terms are treated with 6th-order formulas.

Spatial discretisation of convective terms
As far as the convective terms of the Navier-Stokes system of equations are concerned, the latter are discretised using the central 6th-order energy-preserving method by Kennedy, Grüber, and Pirozzoli (KGP) [27]. As shown by [27], the scheme guarantees that the total kinetic energy is discretely conserved in the limit case of inviscid incompressible flows. The property is also effective on strongly non-uniform Cartesian grids [41,42]; thus, combined with the proposed discretisation form of the viscous components, the entire Navier-Stokes system of equations results in a central strong-conservative formulation without the addition of artificial viscosity.

Temporal integration
Finally, the solution is advanced in time with the 3rd-order low-storage Total Variation Diminishing (TVD) Runge-Kutta scheme proposed by [43]. The method provides stability up to a maximum Courant-Friedrich-Lewy equal to C F L max 1. In the authors' numerical computations, a CFL equal to 0.5 is assumed apart where differently stated. The value represents a good compromise between stability and computational efficiency.

Results and discussion
The present Section aims at validating and discussing the results obtained with the proposed numerical method. The latter is firstly validated concerning its theoretical error decay. Thereafter, a wide variety of flows in laminar conditions is presented as benchmark for severe-varying viscosity conditions. Finally, the numerical assessment of a channel configuration is analysed, providing the effectiveness of the method also in complex three-dimensional turbulent conditions.

Kolmogorov flow
To quantify the error scaling of the proposed interpolation method, a doubly periodic Kolmogorov flow (KF) in laminar and steady-state conditions is solved. The problem consists in forcing the two-dimensional Navier-Stokes system of equations along with a single Cartesian direction with a stationary monochromatic harmonic, Here ω k 2πk/L y is the reduced frequency, k is the wave number, α is the force amplitude, and L y is the length of the domain along with the y-coordinate. Employing a doubly periodic configuration and keeping the system in low-speed and laminar conditions, the entire Navier-Stokes system reduces to the following equilibrium equation: In the case of uniform viscosity, Eq. (15) admits the following analytical solution: The KF is of considerable help in verifying the error scaling of a numerical method since the numerical assessment could rely on a fully periodic setup so that no uncertainties can be ascribed to the numerical treatment of the boundary conditions. To reproduce the problem numerically, the Navier-Stokes system of equations is solved in a square box of size L x × L y L 0 × L 0 . As far as the external force is concerned, the amplitude α is set equal to 0.1u ∞ , and the fourth wave number (k 4) is employed. To keep the system in laminar conditions, μ ∞ and u ∞ are set in such a way that reference Mach and Reynolds numbers of The flow is initialised with zero velocity components, and the non-dimensional temperature T /T 0 and the non-dimensional pressure p/ p 0 are set equal to one. The computation is carried out until the steady-state convergence of the system. To avoid any are retained as metrics of the error scaling, and in particular the p 1 and p 2 definitions are used. Here u(y j ) and u j denote the analytical and the numerical solutions, respectively. Whatever the type of norm is used, the error should asymptotically decay as L M i where C is a constant value and q denotes the error order of accuracy. Thus, comparing the ratio of the error decay, the accuracy order results as with r 2 representing the grid ratio of the computations. Figure 1 reports the error scaling of the present numerical method for a wide variety of formulas accuracy, while the corresponding data are listed in Table 2. From the obtained results, it can be claimed that the current conservative form of the incompressible viscous terms provides the expected high order of accuracy, validating the correctness of the obtained interpolation coefficients.   Due to its inherent bi-periodical arrangement, the KF numerical solution provides a valid benchmark also in studying the convergency of the proposed method even if severe-varying viscosity flows are accounted for. In particular, assuming the viscosity to vary as this profile grants an analytical solution to Eq. (15) which holds as: The problem is treated numerically in a similar manner to the previous one, employing the same Reynolds and Mach regime. Therefore, a uniform grid is used along with the y-coordinate, and the meshes are increasingly refined from a minimum of N y 5 up to a maximum of N y 80 computational nodes. The grid ratio of r is kept equal to 2. Again, the fourth wave number is employed in the forcing term. Since the analytical solution is not differentiable, whatever type of (i) norm used as a metric of the numerical error and (ii) theoretical order of accuracy employed in the discretisation scheme, the latter can scale at most first order (if convergent). Moreover, according to [44], if even just one discontinuity is present in a dynamical system, a non-conservative scheme does not converge to any consistent solution. Thus, the non-conservative approach is expected to fail in reproducing the physical behaviour of the flow in the case of step-gradient viscosity. Figure 2 shows the results of the computation. In particular, a comparison between the standard and the proposed conservative method is presented. The obtained results reveal that a standard co-located discretisation fails in the flow description. As expected, the standard strategy provides an erroneous prediction of the system dynamics, and the method is not able to converge to the analytical reference. On the other hand, the conservative scheme converges to the analytical solution proving an asymptotic accuracy order around unity for all the interpolation formulas. The data related to the error scaling are reported in Table 3.

Poiseuille flow with non-trivial viscosity distributions
To highlight the effectiveness of the present method, a Poiseuille flow (PF) with non-trivial viscosity distributions is considered as the first benchmark in wall-bounded conditions. The computation is performed in a two-dimensional rectangular box with a size equal to L x × L y 4h × 2h along the x-and y-coordinate, respectively, where h is the channel half-height. A uniform points distribution is employed in both the streamwise and wall-normal directions, and in particular, the computations featured an increasingly wall-normal resolution. The number of points, along with the streamwise direction, is kept constant and equal to N x 20. The computational grid is built up in such a way that the first and the last cell interfaces coincide with the walls. Periodic boundary conditions are set along the streamwise direction since the flow infinitely recycles in the channel and is dynamically sustained thought as a constant pressure drop dp/dx. The latter, added as a forcing term to the RHS of the Navier-Stokes system, is defined as where Ma u ∞ /c w is the Mach number of the flow and Re u ∞ h/ν ∞ is the Reynolds number. Here, u ∞ denotes the reference speed, c w is the wall speed of sound, and ν ∞ the reference kinematic viscosity. To avoid any compressible effects and to keep the system in laminar conditions, the latter parameters are set in such a way that a Mach and a Reynolds number equal to 0.1 and 10 are recovered, respectively. The flow is initialised with ambient conditions setting to zero all the velocity components and imposing an initial non-dimensional pressure field p/ p 0 1. The lateral boundaries are modelled with a no-slip isothermal wall discretely enforcing to zero all the velocity components at the locations of the walls and setting to T w /T 0 1 the wall non-dimensional temperature. With these assumptions, the entire Navier-Stokes system of equations reduces to the steady-state equilibrium expressed by the following equation: Equation (23) states that the shear stress τ μdu/dy must be balanced by the pressure drop, independently from the viscosity distribution, resulting in that τ (y) is a linear function of the y-coordinate only. Independently from the employed numerical discretisation, the shear stress behaviour must be granted to accurately resolved the entire flow field. In the present analysis, the Navier-Stokes system of equations is solved until its stationarity. Thus, dynamically converging to the steady-state equilibrium is expressed by Eq. (23).
To convince the reader that a standard co-located discretisation quickly fails the analytical behaviour expressed by Eq. (23), two different viscosity laws are employed. Then the results with two numerical discretisations are compared. The first viscosity law expresses a smooth highly varying profile in the form of where a 0.03, y 2 −y 1 0.5, while the second law provides an intense step-gradient distribution which reads as In Fig. 3, the results related to the first setup (μ 1 ) are given and compared in terms of velocity and stress distributions. Figure 4a and c shows the results obtained with the conservative method, while Fig. 4b and d reports the solution obtained with a standard discretisation as a function of the grid resolution. Compared to the higher-resolved case (solid black line), the classical method can capture the steady-state system equilibrium both in terms of velocity and shear stress. Nevertheless, to capture the analytical distribution accurately, the method requires the grid to place a sufficient number of points inside the intense viscosity-variation region. This Velocity profile and shear stress distributions for a Poiseuille flow with a smooth highly varying viscosity law. a reports the velocity profiles for a Poiseuille flow with sinusoidal viscosity law obtained with the current conservative approach, while b reports the same information obtained with a standard co-located discretisation. c and d draw the shear stress distribution obtained with the two methods as a function of the grid resolution issue is not detectable in the case of the proposed conservative discretisation for which the stress distribution is fully respected regardless of the grid spacing (see Fig. 3). The test can be considered as additional proof on how highly varying diffusive properties act as discontinuities if the grid resolution is not sufficient, thus according to the Hou & Le Floch theorem [44], that leads to non-physical solutions.
Even if the previous results already show the superiority of the conservative approach compared to a standard co-located description, Fig. 4 depicts the behaviour of the same system if a step-gradient (μ 2 ) viscosity law is employed. Again, Fig. 4a and b reports the velocity distributions, while Fig. 4d and c shows the stress profiles. In this case not only a conservative approach can guarantee the correct decaying of the shear stress along with the channel, but even a co-located approach provides an entirely erroneous behaviour of the system producing a non-physical solution.
This result allows for making further comments. Although the authors' method can strictly conserve the viscous stresses in the limit of incompressible or weakly compressible flows, this approach can be very helpful also in shocked embedded conditions. In these cases, upstream of a discontinuity the viscosity field experiences strong gradients related to the high local temperature variations. For this reason, a numerical approximation of viscous terms such as the proposed one is superior also in these situations.

First Stokes problem with different viscosity ratios
To verify the effectiveness of the present method in a time-dependent condition using a non-uniform grid, the so-called First Stokes Problem (FSP) is selected as the benchmark. The problem consists of a y-infinite region Velocity profile and shear stress distributions for a Poiseuille flow with a step-gradient viscosity law. a and b report the velocity distributions for a Poiseuille flow with step-gradient viscosity obtained with the present conservative method and with a standard discretisation, respectively, as a function of the grid resolution. c and d report the comparison of the two methods applied to the shear stress distribution occupied by two different layers of fluids, initially at rest. Suddenly, one layer (the driver) gains a steady speed of u ∞ and starts to perturb the system. Due to the fluid viscosity, a velocity perturbation propagates in the system and diffuses the speed of the driving region. According to this hypothesis, in the limit of incompressible flows, the Navier-Stokes system reduces to the following equation: As a result, the flow speed is a function of the y-coordinate and the time only. Equation (26) states that spatial variation of the shear stress τ μ∂u/∂ y is dynamically balanced by the speed rate ∂u/∂t changing. Thus, a wrong prediction of the shear stress distribution can result in an erroneous behaviour of the entire system. The problem is a classic textbook benchmark since Eq. (26) admits an analytical solution if the two layers are perfectly homogeneous. In this case, the velocity profile is expressed as where erf is the Gauss' error function expressed as and μ ∞ is the uniform reference viscosity. A detailed description of the problem can be found in [45,46]. In the present work, an ad hoc version of the original problem is considered. In particular, as sketched in Fig. 5, two layers of fluid with different viscosity are employed. This setup can be considered as a prototype for shear layer dynamics in laminar conditions and represents a crucial test for the proposed method accounting for both temporal dependency and non-uniform grids. To reproduce the problem numerically, the flow is simulated in a pseudo-infinite configuration employing a twodimensional box of sizes L x × L y L 0 × 40L 0 . Periodic boundary conditions are enforced along with the x-direction, while the field variables are extrapolated at the top and the bottom sides. To avoid any compressible effect, the driver layer is made to flow at a constant speed u ∞ which induces a reference Mach number Ma ∞ u ∞ /c ∞ equal to 0.1. Moreover, the flow is kept in laminar conditions setting the reference kinematic viscosity ν ∞ in such a way that a Reynolds number Re ∞ u ∞ L 0 /ν ∞ of 10 is granted. The computational grid features N x × N y 10 × 120 grid points. A hyperbolic tangent function is employed to cluster the nodes in the interface region. The flow is initialised with zero speed and imposing the non-dimensional temperature T /T 0 and pressure p/ p 0 equal to one, respectively.
Even if an analytical description of the problem is not available for arbitrary viscosity ratios (i.e. R μ 2 /μ 1 ), any consistent physical behaviour would see the system dynamic to change gradually, while the R parameter increases. To verify this statement, the R parameter is varied in the range [1; 5]. Consequently, the dynamics of the system is observed at a reduced time t/t 0 equal to 100. The obtained results are reported in Fig. 6. In particular, Fig. 6a and b shows the velocity distribution as a function of the y-coordinate obtained with the conservative and the standard methods, respectively. Figure 6c and d reports the shear stress distributions. The reader can easily recognise that the standard co-located approach fails in transferring the information from the driver layer to the one at rest, and the driver layer seems to be insensitive to the fluid at its top, especially while the viscosity ratio increases. This behaviour leads to an unstable dynamics for the interface while approaching the highest viscosity ratio (R 5). Thus, the shear stress suddenly changes its sign. On the other hand, the conservative approach grants the proper role to the diffusion, and the interface region continuously smears, dynamically sharing the stress between the two fluid areas. Thus, the presented results show the effectiveness of the conservative method, even in the case of time-dependent problems with non-uniform grids. Therefore, the effectiveness of the proposed method compared to a standard co-located description of the incompressible viscous terms is also confirmed in this test case.

Large eddy simulation of turbulent channel flow
The proposed numerical method is of real help in the field of LES. All of the eddy-viscosity LES models provide an SGS contribution which is always a function of the velocity field and the thermodynamic quantities. It automatically follows that the latter can inherit all the intrinsic nonlinearities to the mechanics of compressible flow, experiencing large local variations and intense gradients. Moreover, dealing with LES, the grid resolution is often far from the scheme convergence, providing under-resolved flow regions of conspicuous extension. For these reasons, the present Section reports how the proposed numerical method better performs in accurately  Figure 6a and b reports the velocity distributions of the problem at a reduced time t/t 0 100 as a function of the viscosity ratio. Figure 6c and d reports the shear distributions at the same time level. The solution obtained with conservative method (left panels) is compared with the one drawn with a standard strategy (right panels). The line colours shade from blue to red increasing the density ratio parameter heading a poorly resolved LES, since the turbulent stresses generated by the SGS model can better propagate in the domain. Thus, the strategy results in a higher accuracy compared to DNS references.
To provide this effectiveness, the numerical assessment of a turbulent channel flow is retained as a benchmark for the proposed approach. In particular, an LES computation is performed in a nearly incompressible condition. Anyhow, independently of the Mach regime, the channel flow configuration does not exhibit compressible features like shocks or discontinuities relying on a proper setup to provide the effectiveness of the proposed method. An LES model with M b u b /c w 0.1 is run and compared with the DNS database of [47]. It is not worthless to be said that dealing with such low Mach conditions in wall-bounded configurations represents an extremely challenging target for a compressible solver, since the acoustical phenomena play a crucial role concerning the aliasing of the pressure field. The problem often results in the instability and the blow-up of the simulations if an accurate numerical description of the flow is not employed. The Reynolds number of the flow is set equal to Re b 2ρ b u b h/μ w 21579. Here, ρ b 1/V V ρdV is the bulk channel density, u b 1/(ρ b V ) V ρudV is the bulk channel speed, while μ w and c w are the wall viscosity and the wall speed of sound, respectively. The configuration corresponds to a friction Reynolds number Re τ ρ w u τ h/μ w of 581 for the conservative method and 593 for the standard discretisation. u τ √ τ w /ρ w is the friction velocity, and τ w is the wall shear stress.
The simulation is carried out in a three-dimensional rectangular box with a size of L x × L y × L z 2π h × 2h × π h. As in the Poiseuille flow, h is the channel half-height. To stress the method, an under-resolved LES is performed. In particular, the mesh features N x × N y × N z 48 × 48 × 48 computational nodes resulting in an internal spacing of x + × y + × z + 77 × (3 ÷ 40) × 38. A uniform grid is used along with the wall parallel directions, while an error mapping function, stretched near the wall regions, is employed in the wall-normal direction. The latter, in particular, is prescribed as follows: Here η j ( j − 1/2)/N y is the computational coordinate, α 2.8 is the stretching parameter, while erf denotes the Gauss' error function as in Eq. (28). Such resolution is not recommended for LES of wall-bounded flows since the numerical scheme hardly captures most of the nonlinear wall turbulent features. However, the combination with high-order and strong-conservative methods has made the computation not excessively impacted by the coarseness of the grid, and the obtained results feature the DNS wall-normal statistics quite accurately.
Periodic boundary conditions are employed along with the homogeneous directions, while a no-slip isothermal wall condition is imposed at the channel walls. The lateral boundaries are modelled with a no-slip isothermal wall enforcing to zero the velocity components and setting to T w /T 0 1 the wall non-dimensional temperature. Like in the Poiseuille flow setup, the mesh is staggered in such a way that the nominal wall location coincides with a cell interface. The process grants strong conservation of all the conservative variables and the viscous stresses. The flow initialisation takes advantage of the [48] procedure. Therefore, to the analytical solution of the Poiseuille flow, a vortex pair is superimposed promoting an early transition to turbulence.
Concerning the forcing term f in the RHS of the Navier-Stokes equations, the latter is evaluated at each time step to enforce a constant mass flow rate discretely. The corresponding power, fu, is added to the RHS of the total energy equation. This arrangement ensures the mass flow rate to be strictly conserved during the time; hence, the bulk Reynolds number Re b is constant during the whole simulation. The simulation is carried out until the flow variables are statistically steady. Figure 7 shows the results of the computation and the comparison between a standard method and the proposed conservative one. In particular, Fig. 7a reports the mean velocity profiles; Fig. 7b shows the Favrefiltered turbulent velocity fluctuations scaled by the friction velocity (i.e. ( u i u i /u 2 τ ) 1/2 ); Fig. 7c depicts the Favre-filtered Reynolds' stresses (i.e. − u v /u 2 τ ) as a function of y + y Re τ coordinate, and finally 7d provides the stress budget for the channel flow. From the obtained results it can be noticed that the proposed numerical method, in combination with wall-consistent LES algebraic model, can accurately predict the DNS reference and the entire highly nonlinear behaviour related to wall-turbulence, even if the resolution is very poor.
In particular, the scaled mean velocity profile u + u/u τ is remarkably accurate, and the velocity fluctuations are in a good agreement with the DNS reference even if the resolution is very poor. However, the most exciting result is the stress budget (Fig. 7d). Since the turbulent viscosity changes by two orders of magnitude in the two first two grid points, the present computation relies on a perfect benchmark for the proposed numerical method. Thus, as shown in Fig. 7d, the conservative form better conserves the turbulent shear stress components along with the channel height, and the overall stress budget shows the scheme much more aligned with the theoretical datum compared with the standard procedure. The result definitively proves the superiority of the present method even in modelling turbulent conditions using a varying eddy viscosity, resulting in a superior way compared to a standard co-located discretisation.

Conclusions
A general strategy to design high-order conservative finite-difference (FD) approximations of viscous/diffusion terms for variable viscosity/diffusion flows is proposed. Results show that a standard co-located discretisation of viscous contributions in laminar conditions fails in the description of flows in the case of highly varying viscosity conditions, even in laminar flows. The proposed method instead grants the high-order proper telescoping of the shear stresses independently of the grid stretching and resolution, resulting in a superior numerical technique even if step-gradient viscosity behaviours are provided. Numerical tests have been performed for laminar and turbulent test cases and nearly incompressible Mach regimes. In particular, from steady laminar flows to unsteady problems, up to complex three-dimensional turbulent flows, the proposed method is always found superior compared to a standard viscous term discretisation when substantial variations of the viscosity are considered. 1 and Re τ 590 with different viscous terms discretisation. a reports the mean streamwise velocity u + u/u τ as a function of inner wall distance y + y Re τ . b shows the scaled velocity fluctuation ( u i u i /u 2 τ ) 1/2 as a function of y + . c reports the scaled Reynolds' stresse − u v /u 2 τ , and finally d provides the comparison of the stress budgets for the standard (S) and the conservative (C) methods. Present data are compared to incompressible DNS results of [47] (black circles) The method seems very promising for large eddy simulations (LES) of turbulent wall-bounded flows where strong eddy viscosity variations occur in the near-wall region, relying on a thoroughly conservative approximation of the Navier-Stokes system of equations in the limit of incompressible or weakly compressible flows. In these concerns, the present scheme for viscous terms, in combination with conservative treatment of convective fluxes [27], allows an accurate representation of the total stress balance.
Since the present formulation provides strong preservation properties of viscous stresses, even in poorly resolved conditions and with highly varying overall viscosity, the method is a promising strategy for LES of complex three-dimensional compressible flows even in case of shocked conditions. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
where {β l } 2 p i 1 are the polynomial coefficients to be determined. The interpolation task is easily translated in finding the coefficients' values which fit the nodal values ofμ tot over a selected stencil S {i − p + 1, . . . , i, . . . , i + p}. Enforcing Eq. (34) to pass through the nodal values ofμ tot we get: This system of equations can be resembled in a matrix formulation as: A mlβl μ tot i− p+m m, l 1, . . . , 2 p.
Hereμ tot i− p+m denotes the array of the nodal values ofμ tot over the stencil S, while A ml is the matrix that collects the grid's metrics. The latter reads as: A ml x i− p+l − x i+1/2 l−1 m, l 1, . . . , 2 p .
With these positions, the interpolation coefficients can be computed as In particular, since Eq. (34) states thatμ tot (x i+1/2 ) β 1 , we obtain concluding thatβ m A −1 m,1 represents the set of the interpolation coefficients which yield to a 2p-order of accuracy for theμ tot i+1/2 value on an arbitrary non-uniform grid. It is not worthless to be said that theβ m β m ∀m if the grid is uniform.
To conclude, since theγ m and theβ m coefficients equal γ m and β m in a uniform grid case, the proposed coefficients fall back to standard FD coefficients if a uniform mesh is employed. Conversely, if non-uniform grids are considered, the method provides the proper cell-bound interpolated values taking into account the local grid variance.