F Nikfarjam

Y Cheny

O Botella

The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries

Keywords: Computational fluid dynamics (CFD), Non-Newtonian fluids, Immersed boundary methods, Cut-cell methods, High-performance computing (HPC)

The LS-STAG method is an immersed boundary/cut-cell method for viscous incompressible flows based on the staggered MAC arrangement for Cartesian grids, where the irregular boundary is sharply represented by its level-set function, results in a significant gain in computer resources (wall time, memory usage) compared to commercial body-fitted CFD codes. The 2D version of LS-STAG method is now well-established (Y. Cheny & O. Botella, J. Comput. Phys. Vol. 229, 1043-1076, 2010), and this paper presents its extension to 3D geometries with translational symmetry in the z direction (hereinafter called 3D extruded configurations). This intermediate step towards the fully 3D implementation can be applied to a wide variety of canonical flows and will be regarded as the keystone for the full 3D solver, since both discretization and implementation issues on distributed memory machines are tackled at this stage of development. The LS-STAG method is then applied to various Newtonian and non-Newtonian flows in 3D extruded geometries (axisymmetric pipe, circular cylinder, duct with an abrupt expansion) for which benchmark results and experimental data are available. The purpose of these investigations are (a) to investigate the formal order of accuracy of the LS-STAG method, (b) to asses the versatility of method for flow applications at various regimes (Newtonian and shear-thinning fluids, steady and unsteady laminar to turbulent flows) (c) to compare its performance with well-established numerical methods (body-fitted and immersed boundary methods).

Introduction

Advanced computational fluid dynamics (CFD) applications concern flows in highly irregular three-dimensional geometries, whose process of grid generation, from the early stage of CAD design to final volume meshing, is very time-consuming and requires highly trained users. Typical unstructured body-conformal meshes are made of elements (cells of hexahedral, tetrahedral shape, . . . ) with irregular connectivity, which makes the solution algorithms extremely memory and CPU intensive. Moreover, it is difficult to assess the quality of the grid (smooth, undistorted cells properly clustered at the wall, . . . ), which may negatively impact the convergence and accuracy of the numerical solution. In contrast, cut-cell/immersed boundary CFD codes (see [START_REF] Mittal | Immersed boundary methods[END_REF] for a review) are based on Cartesian methods which are extremely robust and optimized, and thus requires only a fraction of the computer resources of an unstructured code. Cartesian grid generation is greatly simplified compared to the body-conformal approach, and the grid regularity and quality is not significantly impacted by the complexity of the geometry. Moreover, the computation of flows with moving boundaries can be performed on fixed grids, without the domain remeshing which is required by body-conformal methods.

The major difficulty of traditional IB methods is the treatment of the cut-cells (cells of irregular shape which are formed by the intersection of the Cartesian mesh with the immersed boundary). Classical IB methods such as the momentum forcing method [START_REF] Mittal | Immersed boundary methods[END_REF] discard the discretization of flow equations in the cut-cells, and use special interpolations instead. Thus, strict conservation of mass is not observed, which leads to non-divergence free velocities or non-physical pressure oscillations in the vicinity of the immersed boundary [START_REF] Kang | Accurate and efficient immersed-boundary interpolations for viscous flows[END_REF]. On the other hand, cut-cell methods aim for discretizing the flow equations in cut-cells by ad hoc treatments, whose major drawback is the loss of the 5-points structure (in 2D) of Cartesian methods. Moreover, these methods rely on cell-merging techniques for which 3D extension is known to be non-trivial [START_REF] Mittal | Immersed boundary methods[END_REF]. Therefore, the development of an accurate and cost-effective 3D solver based upon the IB/cut-cell approach is still a challenging task.

In Ref. [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF] we have developed an innovative cut-cell discretization for 2D incompressible flows called "LS-STAG" (Level-Set, STAGgered), which is based on the classical MAC method for Cartesian staggered grids, where the irregular boundary is sharply represented by its level-set function. In three basic types of 2D cut-cells (pentagonal, trapezoidal and triangular, see [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF]), the discretization is designed such as the global invariants of the flow (total mass, momentum and kinetic energy) are preserved at the discrete level. The LS-STAG discretization in the cut-cells is consistent with the MAC discretization used in Cartesian fluid cells, and has the ability to preserve the 5-point Cartesian structure of the stencil, resulting in a highly computationally efficient method. The LS-STAG method has been successfully applied to Newtonian flows at moderate Reynolds number in fixed and moving geometries [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF], pseudoplastic flows [START_REF] Botella | The LS-STAG immersed boundary method for non-Newtonian flows in irregular geometries : Flow of shear-thinning liquids between eccentric rotating cylinders[END_REF], and viscoelastic flows [START_REF] Botella | Application of the LS-STAG immersed boundary/cut-cell method to viscoelastic flow computations[END_REF].

The purpose of this paper is to present a progress report towards a fully 3D LS-STAG code for Newtonian and pseudoplastic liquids, totally optimized for high-performance computing. One of the major difficulty in cut-cell discretization for 3D geometries is the large number of cases to consider: we have enumerated 16 different type of cut-cells in 3D whereas there are only 3 in 2D. In order to mitigate this issue we first consider building the LS-STAG discretization for 3D configurations with translational symmetry in the z direction, that will be subsequently called 3D extruded geometries. This intermediate step towards the fully 3D implementation can be applied to a wide variety of canonical flows (flow in Taylor-Couette reactors, flow past cylinder of arbitrary cross section shape, abrupt axisymmetric expansion/contraction flow, . . . ) and will be regarded as the keystone for the full 3D solver, since both discretization and implementation issues on distributed memory machines are tackled at this stage of development. In these geometries, 3D cut-cells are only extruded extensions of pentagonal, trapezoidal and triangular cut-cells in 2D, which enables us to extend the principles of the energy-conserving discretization of Ref. [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF] to 3D : the discretization of the continuity equation, pressure gradient, convective fluxes, normal stresses are direct extension of the 2D LS-STAG formulas. The discretization of velocity gradient components such as ∂w/∂x and ∂w/∂y, which were absent from the 2D cases, needs further attention. These terms are tentatively discretized with the formula proposed by van der Heiden et al. [START_REF] Van Der Plas | Efficiently simulating viscous flow effects by means of regularization turbulence modeling and local grid refinement[END_REF].

The LS-STAG method is then applied to various Newtonian and non-Newtonian flows in 3D extruded geometries (Taylor-Couette reactor, axisymmetric pipe, circular cylinder, duct with sudden expansion) for which benchmark results and experimental data are available. The purpose of these investigations are (a) to evaluate the formal order of accuracy of the LS-STAG method, (b) to asses the versatility of the method for flow applications at various regimes (Newtonian and shear-thinning fluids, steady and unsteady laminar to turbulent flows) (c) to compare its performance with well-established numerical method (unstructured method, IBM method). First, a grid convergence study is conducted on the Hagen-Poiseuille flow for power-law liquids. A thorough comparison on the accuracy and computational costs (memory usage, CPU time) with respect to a commercial CFD package on unstructured grid (ANSYS FLUENT) is performed on the DFG benchmark for cylinder flow in unsteady laminar regime [START_REF] Bayraktar | Benchmark computations of 3D laminar flow around a cylinder with CFX, OpenFOAM and FeatFlow[END_REF]. The ability to perform computation of high Reynolds number flow is addressed on the (weakly) turbulent cylinder flow at Re = 300 [START_REF] Mittal | On the inclusion of three-dimensional effects in simulations of twodimensional bluff body wake flows[END_REF][START_REF] Kravchenko | B-spline method and zonal grids for simulations of complex turbulent flows[END_REF], for which comparison with the results of a well-established IB method [START_REF] Balaras | Modeling complex boundaries using an external force field on fixed cartesian grids in large-eddy simulations[END_REF] is provided. The paper is concluded by studying the flow of shear-thinning fluids in an axisymmetric sudden expansion geometry, for which rheological characterization, experimental flow measurements and FLUENT c simulations have been performed in [START_REF]Comportement de fluides complexes sous écoulement : Approche expérimentale par résonance magnétique nucléaire et techniques optiques et simulations numeriques[END_REF].

Governing equations for non-Newtonian flows

Let Ω be a 3D computational box domain and Γ its surface. The governing equations are the incompressible Navier-Stokes equations for non-Newtonian fluids written in conservative form. The continuity equation reads :

Γ v • n dS = 0, (1)
where v = (u, v, w) is the velocity, and the momentum equations in the x, y and z directions are respectively :

ρ d dt Ω u dV + ρ Γ (v • n) u dS + Γ p e x • ndS - Γ η( γ)( γxx e x + γxy e y + γxz e z ) • ndS = 0, (2a) ρ d dt Ω v dV + ρ Γ (v • n) v dS + Γ p e y • ndS - Γ η( γ)( γyx e x + γyy e y + γyz e z ) • ndS = 0, (2b) ρ d dt Ω w dV + ρ Γ (v • n) w dS + Γ p e z • ndS - Γ η( γ)( γzx e x + γzy e y + γzz e z ) • ndS = 0, ( 2c
)
where ρ is the fluid density, p is the pressure and γ = ∇v + ∇v T is the rate-of-strain tensor.

In the following two rheological models are employed to describe the shear-dependent dynamic viscosity η( γ), where the scalar shear-rate γ is defined as :

γ = 1 2 γ : γ = [ γxx ] 2 + [ γyy ] 2 + [ γzz ] 2 + [ γxy ] 2 + [ γxz ] 2 + [ γyz ] 2 . ( 3
)
The simplest one is the power-law (Ostwald-de Waele) model :

η( γ) = k γn-1 , (4)
where k and n are the consistency and power-law index respectively. Based on the value of the power-law index, the fluid exhibits 3 distinct non-Newtonian behaviors: shear-thinning or pseudoplastic (n < 1), Newtonian (n = 1) and shear-thickening or dilatant (n > 1). The powerlaw model is widely used in engineering applications and its simplicity allows to derive analytical solutions of the Navier-Stokes equations for canonical flows [START_REF] Bird | Dynamics of Polymeric liquids[END_REF]. This rheological model will be used for evaluating the accuracy of the LS-STAG method in Section 6.1. However the powerlaw model describes only approximatively the behavior of a real non-Newtonian fluid and most notably fails to predict the Newtonian plateaus that reaches the viscosity at low (η ∼ = η 0 ) and high (η ∼ = η ∞ ) shear-rates. These limitations are overcomed by the Cross model :

η( γ) = η ∞ + (η 0 -η ∞ ) 1 + (K γ) m , ( 5
)
where K and m are 2 fitting parameters. The Cross parameters can be correlated to the power-law parameters (4) by observing that k = η 0 K -m and n = 1 -m in the limiting case where η 0 η ∞ , K γ 1 and η ∞ is small. The Cross model will be used for evaluating the LS-STAG method against experimental results in Sect. 6.3.

The LS-STAG mesh for 3D extruded geometries

We consider a 3D irregular solid domain Ω ib which is embedded in the computational domain Ω, such that Ω f = Ω \ Ω ib represents the fluid domain where the flow equations are to be discretized. To keep track of the irregular boundary Γ ib , we employ a signed distance function φ(x) (i.e. , the level-set function [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF]) such that φ(x) is negative in the fluid region Ω f , φ(x) is positive in the solid region Ω ib , and such that the boundary Γ ib corresponds to the zero level-set of this function. This leads to the extension of the well-known MAC mesh (see e.g. [START_REF] Verstappen | Symmetry-preserving discretization of turbulent flow[END_REF]) that will be subsequently referred to as the LS-STAG mesh. In the case of irregular geometries extruded in the z direction, the computational cells Ω i,j,k are either Cartesian cells of size ∆x i × ∆y j × ∆z k or cut-cells, whose generic types (pentagonal, trapezoidal and triangular) are depicted in Figure 1. For all cell types, the velocities unknowns are located at face centers and the pressure p i,j,k is at cell centroids. In the cut-cell shown in Figure 2, the immersed boundary is represented by the rectangular plane Γ ib i,j,k whose extremities are defined by linear interpolation of the level-set function φ i,j,k = φ(x i , y j , z k ) at cell corners. This leads to the following decomposition of the surface Γ i,j,k of a cut-cell :

Γ i,j,k = Γ w i,j,k ∪ Γ e i,j,k ∪ Γ s i,j,k ∪ Γ n i,j,k ∪ Γ b i,j,k ∪ Γ t i,j,k ∪ Γ ib i,j,k , (6)
by using the usual compass notations for Cartesian discretizations [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. In [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF], the level-set function was employed to obtain analytical expressions of the geometry parameters of 2D cut-cells, such as their face areas and volume. These definitions are now extended to the 3D case, by calculating first the fluid portion of the cut-cell edges. For the cut-cell of Figure 2, the length of Edge E y i,j,k = {(x i , y, z k ), y ∈ [y j-1 , y ib ]} is y ib -y j-1 = α y i,j,k ∆y j , where the cell-edge fraction ratio α y i,j,k is computed by linear interpolation of φ(x i , y, z k ) as :

α y i,j,k = φ i,j-1,k φ i,j-1,k -φ i,j,k , (7)
since φ(x i , y ib , z k ) = 0. Note that cell-edge fraction ratios take values in [0, 1] by construction, and due to extrusion α x i,j,k-1 = α x i,j,k , α y i,j,k-1 = α y i,j,k and α z i,j,k is either 0 or 1. The cell-edge fraction ratios are subsequently used to calculate the fluid area of cut-cell faces. For example, the fluid area of the top face Γ t i,j,k in Figure 2 is θ w i,j,k ∆x i ∆y j where the cell-face fraction ratio is :

θ w i,j,k = 1 - 1 2 (1 -α x i,j,k )(1 -α y i,j,k ), (8)
for a pentagonal face. Analytical expressions for trapezoidal and triangular faces (which correspond to volumes in a 2D mesh) are given in [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF]. Finally, the volume V i,j,k of a cut-cell is simply :

V i,j,k = θ w i,j,k ∆x i ∆y j ∆z k . ( 9
)
since each cut-cell is a prism of base Γ t i,j,k due to extrusion. These geometry parameters will extensively be used to discretize the flow equations in the next sections.

Discretization of the continuity equation

The starting point of the discretization concerns the continuity equation [START_REF] Mittal | Immersed boundary methods[END_REF], which is written as the net mass flux through each elementary face of control volume (CV) Ω i,j,k according to (6) :

ṁi,j,k ≡ -u i-1,j,k + u i,j,k -v i,j-1,k + v i,j,k -w i,j,k-1 + w i,j,k + U ib i,j,k = 0. ( 10
)
In the following, we consider homogeneous boundary conditions at the immersed boundary Γ ib , so that the mass flux through the immersed boundary

U ib i,j,k ≡ Γ ib i,j,k v ib • n ib i,j
,k dS is zero. The mass fluxes through fluid faces (denoted with a bar) are approximated with midpoint quadrature, yielding for example at the top face :

w i,j,k ≡ Γ t i,j,k v • e z dS ∼ = θ w i,j,k ∆x i ∆y j w i,j,k . ( 11
)
After similar discretizations at the other faces, the discrete continuity equation reads :

ṁi,j,k ∼ = θ u i,j,k u i,j,k -θ u i-1,j,k u i-1,j,k ∆y j ∆z k + θ v i,j,k v i,j,k -θ v i,j-1,k v i,j-1,k ∆x i ∆z k + θ w i,j,k w i,j,k -θ w i,j,k-1 w i,j,k-1 ∆x i ∆y j = 0, (12)
and is valid for any type of fluid cells.

Momentum discretization

Momentum discretization will be illustrated for the x-momentum equation (2a) in the staggered CV Ω u i,j,k represented in Figure 3. This CV is centered around u i,j,k and is formed by halves of neighboring cut-cells Ω i,j,k and Ω i+1,j,k , such as its volume is 1 2 V i,j,k + 1 2 V i+1,j,k and its surface is decomposed as :

Γ u i,j,k = Γ u,w i,j,k ∪ Γ u,e i,j,k ∪ Γ u,s i,j,k ∪ Γ u,ib i,j,k ∪ Γ u,b i,j,k ∪ Γ u,t i,j,k . ( 13
)
Figure 3: Staggered arrangement of the variables for two adjacent cut-cells on the LS-STAG mesh. CV Ω u i,j,k ∈ Ω i,j,k ∪ Ω i+1,j,k is indicated in red dashed line (--).

Basics of the LS-STAG discretization

As in 2D [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF][START_REF] Botella | The LS-STAG immersed boundary method for non-Newtonian flows in irregular geometries : Flow of shear-thinning liquids between eccentric rotating cylinders[END_REF], the key idea of the LS-STAG discretization is to enforce the strict conservation of total mass Ω f ∇ • v dV , total momentum P(t) = ρ Ω f v dV and total kinetic energy E c (t) = 1 2 Ω f |v| 2 dV (when viscosity η becomes negligible) in each of the half CVs of Ω u i,j,k independently, such that any combination of half CVs yields a consistent discretization with the aforementioned global conservation properties. This results to a discrete pressure gradient which is dual to the discrete divergence [START_REF] Balaras | Modeling complex boundaries using an external force field on fixed cartesian grids in large-eddy simulations[END_REF], a skew-symmetric discretization for the convective terms, and a viscous stress discretization which is positive definite. Some of these terms are direct extensions of the 2D formulas: the discrete divergence (10), the discrete pressure gradient which reads in the x-direction :

Γ u i,j,k p e x • n dS ∼ = [G x P ] i,j,k = θ u i,j,k ∆y j ∆z k (p i+1,j,k -p i,j,k ) , ( 14
)
and the skew-symmetric discretization of the convective terms, which amounts to using centered approximations with constant weights irrespective of the cell size [START_REF] Verstappen | Symmetry-preserving discretization of turbulent flow[END_REF], for example at top face Γ u,t i,j,k :

Γ u,t i,j,k u (v • e z ) dS ∼ = u i,j,k + u i,j,k+1 2
w i,j,k + w i+1,j,k 2 . ( 15
)
These formulas are valid for any type of fluid cells, and in the particular case of Cartesian fluid cells (such that the cell-face fraction ratios are equal to 1), the MAC discretization on non-uniform Cartesian grids of Verstappen & Veldman [START_REF] Verstappen | Symmetry-preserving discretization of turbulent flow[END_REF] is recovered. In contrast, the discretization of the viscous stresses in extruded cut-cells introduces additional difficulties compared to the 2D case, which will be discussed in the next sections.

Discretization of the velocity gradient

We first consider the discretization of the velocity gradient :

∇v =         ∂u ∂x ∂v ∂x ∂w ∂x ∂u ∂y ∂v ∂y ∂w ∂y ∂u ∂z ∂v ∂z ∂w ∂z         , ( 16
)
which plays an important role in the numerical modeling of the viscous fluxes presented in the next section since it defines the strain-rate tensor γ and the shear-rate γ from which the local viscosity is evaluated. Due to the staggering of the velocity unknowns, the diagonal and offdiagonal components of ∇v are respectively located at centroids and edges of the computational cells (see Figure 1), leading thus to distinct discretization formulas.

The computation of the diagonal components of ∇v at cell centroids is performed as in the 2D case by enforcing the Green theorem Ω ∇ • v dV = Γ v • n dS at the discrete level [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF]. This leads to a unique formula in all types of cut-cells, e.g. for ∂u/∂x :

∂u ∂x i,j,k ∼ = θ u i,j,k u i,j,k -θ u i-1,j,k u i-1,j,k V i,j,k /(∆z k ∆y j ) , (17)
which reduces to a standard two-point formula in Cartesian cells.

An issue related to the 3D extension of the LS-STAG method concerns the computation of the off-diagonal components of ∇v at the edge midpoints where the orthogonality property of the MAC mesh is lost [START_REF] Moukalled | The finite volume method in computational fluid dynamics[END_REF]. In the case of extruded geometries, it only concerns terms ∂w/∂x and ∂w/∂y in the cut-cells (see Figure 3, where Edge E y i,j,k and the line connecting w i,j,k and w i+1,j,k is not orthogonal), where the use of a standard two-point finite-difference formula would be inaccurate unless a non-orthogonal correction is taken into account [START_REF] Moukalled | The finite volume method in computational fluid dynamics[END_REF]. The design of such a scheme, which requires data interpolation and leads to an enlarged stencil near the cut-cells is left to future work. Instead, we employ the discretization of Van der Plas et al. [START_REF] Van Der Plas | Efficiently simulating viscous flow effects by means of regularization turbulence modeling and local grid refinement[END_REF] who suggest to compute these terms by analogy with the pressure gradient in 2D :

∂w ∂x i,j,k ∼ = α y i,j,k ∆y j (w i+1,j,k -w i,j,k ) 1 2 ( θ w i+1,j,k ∆x i+1 ∆y j + θ w i,j,k ∆x i ∆y j ) , ( 18
)
and ∂w/∂y is computed similarly. The impact of this ad hoc two-point formula on the overall accuracy of the method will be investigated in details in Section 6.1.

For other terms standard 2-point formulas are used, such as for ∂u/∂y at the midpoint of Edge E z i,j-1,k in Figure 3 :

∂u ∂y i,j-1,k ∼ = u i,j,k -u i,j-1,k ∆y , ∆y = 1 2 α y i,j,k ∆y j + 1 2 ∆y j-1 , (19)
and at the midpoint of solid edge E z i,j,k :

∂u ∂y i,j,k ∼ = 0 -u i,j,k ∆y ib , ∆y ib = 1 2 α y i,j,k ∆y j . ( 20
)
for the case of homogeneous boundary conditions. However, for terms such as ∂v/∂x at the solid boundary, similar formulas cannot be derived since θ v i,j,k and θ v i+1,j,k are both null. To address this issue, we have derived a discretization of ∂v/∂x| i,j,k by enforcing that the Schwarz theorem (symmetry of second derivatives) :

∂ ∂y ∂v ∂x = ∂ ∂x ∂v ∂y , ( 21
)
holds discretely on cell face Γ e i,j,k . Firstly, the LHS of ( 21) is readily discretized as :

∂ ∂y ∂v ∂x i,j,k = 1 2∆y ib ∂v ∂x i,j,k - ∂v ∂x i,j-1,k , ( 22
)
where ∂v/∂x| i,j-1,k has an expression similar to [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF]. Secondly, the RHS of ( 21) is discretized by analogy with the pressure gradient ( 14) since ∂v/∂y is also located at cell centroids, leading to :

∂ ∂x ∂v ∂y i,j,k = θ u i,j,k ∆y j ∆z k 1 2 V i,j,k + 1 2 V i+1,j,k ∂v ∂y i+1,j,k - ∂v ∂y i,j,k . ( 23
)
Finally, the identification of ( 22) with [START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF] gives the expression of ∂v/∂x| i,j,k :

∂v ∂x i,j,k = ∂v ∂x i,j-1,k + θ u i,j,k ∆y j ∆z k ∆y ib 1 4 V i,j,k + 1 4 V i+1,j,k ∂v ∂y i+1,j,k - ∂v ∂y i,j,k , ( 24
)
which is completed with discretization formulas for ∂v/∂y| i,j,k and ∂v/∂y| i+1,j,k analogous to [START_REF] Gropp | Using MPI: Portable Parallel Programming with the Message-Passing Interface[END_REF].

Discretization of the viscous fluxes

The discrete velocity gradients are now used for discretizing the viscous fluxes of the xmomentum equation (2a) :

F v u = Γ u i,j,k η( γ) γxx e x • n dS + Γ u i,j,k η( γ) γxy e y • n dS + Γ u i,j,k η( γ) γxz e z • n dS. ( 25
)
As in our previous works [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF][START_REF] Botella | The LS-STAG immersed boundary method for non-Newtonian flows in irregular geometries : Flow of shear-thinning liquids between eccentric rotating cylinders[END_REF], the normal stress flux is computed by analogy with the pressure gradient ( 14) :

Γ u i,j,k η( γ) γxx e x • n dS ∼ = θ u i,j,k ∆y j ∆z k 2η( γe ) ∂u ∂x i+1,j,k -2η( γw ) ∂u ∂x i,j,k , ( 26
)
and the discretization is completed with Eq. ( 17).

For the shear stress fluxes, midpoint quadrature is applied to each of the faces of Γ u i,j,k (given by ( 13)), leading to:

Γ u,s i,j,k ∪Γ u,ib i,j,k η( γ) γxy e y • n dS ∼ = ∆S ib η( γn ) ∂u ∂y + ∂v ∂x i,j,k -∆S s η( γs ) ∂u ∂y + ∂v ∂x i,j-1,k , (27a) Γ u,t i,j,k ∪Γ u,b i,j,k η( γ) γxz e z • n dS ∼ = ∆S t η( γt ) ∂u ∂z + ∂w ∂x i,j,k -∆S b η( γb ) ∂u ∂z + ∂w ∂x i,j,k-1 . (27b)
In these formulas, ∆S ib = 1 2 (∆x i + ∆x i+1 ) is the projected area of the solid face along the y axis, and the fluid face areas are evaluated by mean of the cell-face fraction ratios, for instance the surface area of Γ u,t in Equation (27b) reads :

∆S t = 1 2 ( θ w i,j,k ∆x i + θ w i+1,j,k ∆x i+1 )∆y j , ( 28
)
which ensures local flux conservativity across fluid faces.

For evaluating the non-Newtonian viscosity in ( 26) and ( 27), the shear-rate γ (3) has to be computed at cell centroids for the normal stresses (e.g. γw in ( 26)) and at edge midpoints for the shear stresses (e.g. γs in [START_REF] Rajani | Numerical simulation of laminar flow past a circular cylinder[END_REF]). This is achieved by first computing γ at cell centroids, for example:

γw = γi,j,k ∼ = γxx | 2 i,j,k + γyy | 2 i,j,k + γzz | 2 i,j,k + γxy 2 i,j,k + γxz 2 i,j,k + γyz 2 i,j,k , ( 29
)
where the off-diagonal terms are interpolated from values at the surrounding edges by arithmetic average [ • ]. Then, the shear-rate at cell edges is obtained from its neighboring cell-centroid values by using the volume-weighted interpolation first employed in [START_REF] Botella | Application of the LS-STAG immersed boundary/cut-cell method to viscoelastic flow computations[END_REF] in the framework of viscoelastic flow simulations. For example, the value of γs at Edge E z i,j-1,k is interpolated as :

γs = [β V γ] i,j,k + [β V γ] i-1,j,k + [β V γ] i-1,j-1,k + [β V γ] i,j-1,k / V i,j,k , (30)
where V i,j,k is the weighted sum of the neighboring cell volumes :

V i,j,k = β i,j,k V i,j,k + β i-1,j,k V i-1,j,k + β i-1,j-1,k V i-1,j-1,k + β i,j-1,k V i,j-1,k , (31)
and the coefficients β i,j,k depend on the type of cell Ω i,j,k :

β i,j,k =    0 Solid cells, 1/3 Triangular cells, 1/4 Otherwise, (32)
These coefficients have been defined in [START_REF] Botella | Application of the LS-STAG immersed boundary/cut-cell method to viscoelastic flow computations[END_REF] for deriving a conservative discretization of the elastic shear stresses, and then successfully employed in [START_REF] Botella | The LS-STAG immersed boundary method for non-Newtonian flows in irregular geometries : Flow of shear-thinning liquids between eccentric rotating cylinders[END_REF] for 2D non-Newtonian computations.

Solution strategy and parallel implementation

The time-advancement of the Navier-Stokes equations is based on the IMEX fractional-step method described in [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF] for Newtonian flows and extended in [START_REF] Botella | The LS-STAG immersed boundary method for non-Newtonian flows in irregular geometries : Flow of shear-thinning liquids between eccentric rotating cylinders[END_REF] to non-Newtonian flows. For the steady computations reported in the next section, first-order time stepping with Euler schemes is used until the solution reaches steady-state. Unsteady computations of Newtonian flows are performed with second-order methods (BDF2 time-stepping with Adams-Bashforth extrapolations).

The LS-STAG flow solver is a Message Passing Interface (MPI) [START_REF] Gropp | Using MPI: Portable Parallel Programming with the Message-Passing Interface[END_REF] code written in FORTRAN-90 where parallelism is achieved by performing a Cartesian decomposition of the domain. Each cuboid-shaped sub-domain of the decomposition is identified with a coordinate triplet (I, J, K) and is surrounded by one overlapping layer of buffer cells (see Figure 4), which allows to compute the discrete fluxes at the boundary of the sub-domain. The data in the layer of buffer cells are updated from the neighboring sub-domains (I ± 1, J ± 1, K ± 1) by using border exchange communication MPI routines (MPI_SENDRECV) at the end of each key step of the fractional-step algorithm, as summarized in Figure 5. The resulting sequence of operations is the same as for any standard Cartesian method, which notably shows that the LS-STAG cut-cell method introduces no specific issues related to parallelism.

The HYPRE [START_REF] Falgout | Center for Applied Scientific Computing (CASC)[END_REF] library (V2.10) is employed for the iterative solution of linear systems, most notably the Poisson pressure equation. All the simulations reported in the next section were performed on the computer cluster of Institut Jean Lamour (Université de Lorraine) with an infiniband network and an Intel R Xeon E5-2660 2.20 Ghz processor. A comparison of the LS-STAG solver performance with the commercial CFD code Ansys-FLUENT c will be given in Section 6.2.1.

Results and discussion

The LS-STAG flow solver is now evaluated on various benchmark flows of Newtonian and non-Newtonian fluids in 3D extruded geometries, such as Hagen-Poiseuille flow, flows past a circular cylinder and through a duct with an abrupt expansion. For all computations, the Reynolds number is defined as

Re D = ρU 2-n ref D n /K ( 33
)
where D and U ref are reference length and velocity, n > 0 and K are the power-law parameters (Sect. 2).

Accuracy study of the LS-STAG method

The spatial accuracy of the LS-STAG method for 2D non-Newtonian fluids has thoroughly been studied on Taylor-Couette flows in [START_REF] Botella | The LS-STAG immersed boundary method for non-Newtonian flows in irregular geometries : Flow of shear-thinning liquids between eccentric rotating cylinders[END_REF], where second-order accuracy has been observed in Cartesian cells away from the complex geometry and super-linear convergence in the cut-cells. Although a Taylor-Couette cell is an extruded geometry, the independence of exact solution in the axial direction makes impossible to evaluate the accuracy of LS-STAG discretization in that direction. Hence, a more suitable test-case is the Hagen-Poiseuille flow, where a constant pressure drop drives the flow through a circular cylinder of radius R with no-slip boundaries. For power-law fluids, the analytical solution reads in polar coordinates:

w(r) = w c 1 - 2r D n+1 n , ( 34
)
where w c is the streamwise velocity at the centerline. This test-case has already been considered in [START_REF] Van Der Plas | Efficiently simulating viscous flow effects by means of regularization turbulence modeling and local grid refinement[END_REF] for a Newtonian fluid and, following their recommendations, the pipe dimensions are 2R × L z = 1×2 with periodic boundary conditions at the streamwise boundaries. The constant pressure drop is modeled as a source term in the Navier-Stokes equations. Computations are performed at Re D = 2000, where D and w c are the reference length and velocity, on a series of refined meshes with square cells of size h. In order to evaluate the LS-STAG cut-cell treatment, we have also performed computation with the staircase method which corresponds to a stepwise approximation of the irregular boundaries with Cartesian cells. As a result the staircase treatment only affects the discretization in cut-cells, while it remains unchanged in Cartesian cells. Figure 6 compares the errors on the streamwise velocity measured in the maximum norm on the whole computational domain and in Cartesian cells away from the irregular boundary for both Newtonian (n = 1) and shear-thinning (n = 0.5) fluids. As in the 2D case [START_REF] Botella | The LS-STAG immersed boundary method for non-Newtonian flows in irregular geometries : Flow of shear-thinning liquids between eccentric rotating cylinders[END_REF], the LS-STAG errors are more than a decade lower than those of the first-order staircase method, which highlights the benefits of our cut-cell treatment. Second order accuracy is obtained in Cartesian cells away from the irregular boundary for both fluid types. The magnitude of errors are slightly higher for the power-law fluids due to additional errors introduced by the computation of the non-Newtonian viscosity. However, compared to the 2D case, the error on the whole domain (which actually measures the error in cut-cells) reduces from superlinear to linear only. This behavior, also observed in the cut-cell method of Ref. [START_REF] Van Der Plas | Efficiently simulating viscous flow effects by means of regularization turbulence modeling and local grid refinement[END_REF], is certainly caused by the ad hoc computation of the off-diagonal gradients ∂w/∂x and ∂w/∂y (see Sect. 5.2) which were absent from the 2D version. Further insight on the structure of discretization error is given in Figure 7, which displays the pointwise error on mesh points along cross-section z = 1 for h = 5.5 × 10 -3 . First, it can be observed that the largest errors are localized near the complex boundary and, for the LS-STAG method, the errors drop to several orders of magnitude away from the cut-cells. In contrast, the crude treatment of the staircase method pollutes the solution in the whole fluid domain. For non-Newtonian fluids, the difference between the staircase and LS-STAG method is less pronounced. This trend was also observed in the 2D version [START_REF] Botella | The LS-STAG immersed boundary method for non-Newtonian flows in irregular geometries : Flow of shear-thinning liquids between eccentric rotating cylinders[END_REF], and is attributed to the additional errors introduced by the evaluation of the non-Newtonian viscosity in the viscous fluxes.

Flows past circular cylinder

In this part, we consider the flow past circular cylinder at two different flow regimes: an unsteady laminar flow with time-varying Reynolds number and a weakly turbulent flow. The computational domain is depicted in Figure 8 ; its dimensions and conditions at the boundaries will be given below for each configuration.

Comparison with unstructured body-fitted solver Ansys-FLUENT c

The purpose of this part is to address 2 crucial questions regarding flow computations in complex geometries which, to the best of our knowledge, have seldomly been addressed in the literature: how accurate is the solution given by a cut-cell/IBM code compared to a body-fitted unstructured solver of same formal order of accuracy? How do the computational resources of the cut-cell/IBM code (memory usage, CPU time) compare with those of the body-fitted solver?

For that purpose, we have conducted a series of computations with the LS-STAG code and the commercial CFD package FLUENT c (V14.57) on the benchmark 3D-3Z initially proposed in [19], which concerns the unsteady laminar flow past a circular cylinder in a confined square duct. The computational domain is depicted in Figure 8, with spanwise and cross-stream dimensions H = 0.41 m and L z = 0.41 m, respectively. The upstream and downstream distances of the cylinder (of diameter D = 0.1 m) are set to L u = 0.45 m and L z = 1.95 m with an 1.0 cm offset in the cross-stream direction (H d = 0.15 m and H u = 0.16 m). At the spanwise and cross-stream boundaries, no-slip conditions are applied, while a parabolic stream-wise velocity profile is set at the inlet:

u inlet = 16U max yz sin( πt 8 ) (H -y)(H -z) H 4 , U max = 2.25 m/s. ( 35
)
At the outlet, a convective outflow condition [START_REF] Orlanski | A simple boundary condition for unbounded hyperbolic flows[END_REF] is applied. The fluid viscosity is ν = 10 -3 m 2 /s, yielding a time-varying Reynolds number Re(t) ∈ [0, 100] based on the cylinder diameter D and the mean flow velocity U ref = 4U max /9 at the inlet. For this configuration, Bayraktaret al. [START_REF] Bayraktar | Benchmark computations of 3D laminar flow around a cylinder with CFX, OpenFOAM and FeatFlow[END_REF] recently conducted a thorough comparison of the performances of 3 unstructured CFD packages (OpenFOAM, Ansys-CFX and FeatFlow) on a sequence of body-fitted curvilinear grids and provided highly accurate benchmark results. The same grids have been used for our FLUENT c computations: a xy view of the coarsest mesh we considered, labeled L2 in [START_REF] Bayraktar | Benchmark computations of 3D laminar flow around a cylinder with CFX, OpenFOAM and FeatFlow[END_REF], is shown in Figure 9(a). The refined grid layer around the cylinder surface consists in near-uniform hexahedral cells of size ∆ 2 /D 0.31 (based on the main diagonal) at the cylinder surface. The 3D mesh is obtained by extrusion in the z direction with a uniform distribution of cells, which gives a number of 49152 3D cells. Finer meshes of cell size ∆ 2 /2 and ∆ 2 /4 are obtained by successive refinement, yielding meshes L3 and L4 with 393216 and 3145728 cells respectively. Corresponding block-structured Cartesian grids M2-M4 have been constructed for the LS-STAG solver with the following constraints: the total number of fluid cells and the grid resolution near the complex boundary should be as close as possible to the curvilinear meshes L2-L4, in spite of the inferior flexibility of Cartesian grids. The xy view of grid M2 is shown in Figure 9(b) : grid clustering around the cylinder is obtained by defining a rectangular block of L x /D × L y /D = 4.0 × 2.0 with uniform cells of size h x /D × h y /D = 0.118 × 0.105, while the grid size is relaxed upstream and downstream of this block. In the spanwise direction, a uniform cell distribution of size ∆z/D = 0.256 is used. The resulting 3D mesh has 49136 fluid cells of average size ∆ 2 /D 0.30 near the cylinder body, which is similar to the L2 mesh. The main characteristics of the M2-M4 meshes are given in Table 1, along with the time-step used for the LS-STAG and FLUENT c computations.

The FLUENT c computations have been performed with the pressure-based solver, with second-order spatial discretizations [21]: QUICK scheme for convection discretization, Green-Gauss node-based interpolation for gradient calculation, and PRESTO! for pressure interpolation. Second-order implicit time-stepping (BDF2) has been selected. For pressure-velocity coupling, 2 variants of the segregated PISO time-advancement have been evaluated: the iterative (default) and non-iterative (NITA) versions. The default PISO algorithm performs several inner iterations per time-step to solve the momentum and pressure equations while the NITA variant solves these equations only once per time-step, yielding thus a time-advancement which resembles the most to our AB/BDF2 fractional step method. All other numerical parameters have been left to their default settings.

As in [START_REF] Bayraktar | Benchmark computations of 3D laminar flow around a cylinder with CFX, OpenFOAM and FeatFlow[END_REF][START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF], the LS-STAG and FLUENT c simulations are conducted for the time interval [0, 8s] corresponding to a half period of the inflow condition (35), and are evaluated against a highly accurate FeatFlow solution obtained on a fine L5 mesh [22]. Figure 10 displays the time evolution of the force coefficients obtained on the various meshes. First, it can be seen that both LS-STAG and FLUENT c solvers fail to capture the flow dynamics on the M2/L2 mesh due to insufficient grid resolution. Accurate results are obtained on the 2 finer meshes although, as shown in the close-ups of Figure 10, the LS-STAG results appear to be slightly inferior to FLUENT c . A quantitative comparison is provided by Table 2 where, as proposed in Ref. [START_REF] Bayraktar | Benchmark computations of 3D laminar flow around a cylinder with CFX, OpenFOAM and FeatFlow[END_REF], Finally, Table 3 presents the run-time and memory usage for the LS-STAG and FLUENT c computations. Due to the limited availability of FLUENT c licenses in our institution, the comparisons have been performed on 2 processors for the coarse mesh computation (M2 and L2) and 4 processors for the two finer meshes. The FLUENT c computations with the default PISO algorithm proved to be exceedingly time consuming (more than 22 times slower than LS-STAG) due to the cost of inner iterations and have not been pursued on the finest mesh. Compared to the NITA variant of the PISO algorithm, the LS-STAG code is still substantially faster (approximately 3 to 4 times faster than FLUENT c ). The memory savings of the LS-STAG methods are also sub- stantial (4.5 times less memory needed than FLUENT c on M2/L2 and superior to 3 times for the finer meshes) due to the on-the-fly computation of cell geometry parameters of the LS-STAG code. Finally, we mention that these comparisons are made with the current implementation of LS-STAG code that uses the HYPRE library for solving the linear systems. Since these linear systems typically take up to 60% of the computational cost of a simulation, the optimization of LS-STAG code will always be limited by the performances of its algebraic solver package.

Turbulent flow at Re D = 300

The purpose of this section is to assess the ability of LS-STAG method to perform accurate direct numerical simulation (DNS) of turbulent flows in complex geometries. The selected testcase is the (weakly) turbulent flow at Re = 300, at beginning of the wake turbulent shedding regime [START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF], for which accurate treatment of the fine scale structures in the near-cylinder wake is needed to properly capture the flow dynamics. It has been selected due to the availability of experimental results [START_REF] Norberg | Flow around a circular cylinder: Apects of fluctuating lift[END_REF] and accurate body-conforming DNS (the spectral simulations of Mittal & Balachandar [START_REF] Mittal | On the inclusion of three-dimensional effects in simulations of twodimensional bluff body wake flows[END_REF] and the zonal B-splines computations of Kravchenko et al. [START_REF] Kravchenko | B-spline method and zonal grids for simulations of complex turbulent flows[END_REF]). These studies will serve as reference for our simulations. In addition, this flow has also been computed by a state-of-the-art immersed boundary method [START_REF] Balaras | Modeling complex boundaries using an external force field on fixed cartesian grids in large-eddy simulations[END_REF]. In order to address the comparison of LS-STAG method with this IBM computation, we have used the same domain definition (L = 50D, H = 30D and L z = 6D in Figure 8) and built Cartesian grids with similar resolution. The coarse mesh M1 (correspond to the Case 1 mesh in [START_REF] Balaras | Modeling complex boundaries using an external force field on fixed cartesian grids in large-eddy simulations[END_REF]) has 243 × 133 × 64 cells and a grid distribution similar to the mesh displayed in Figure 9(b): the grid distribution is uniform in the spanwise direction, and in the xy plane a uniform block of size 4D × 4D is used around the circular cylinder with grid resolution h/D = 3.88 × 10 -2 . A finer mesh M2 of 422 × 226 × 64 cells with grid resolution h/D = 1.94 × 10 -2 near the cylinder (similar to Case 2 in [START_REF] Balaras | Modeling complex boundaries using an external force field on fixed cartesian grids in large-eddy simulations[END_REF]) has also been considered.

The boundary conditions are the following: uniform velocity u = U ref at the inlet boundary, periodicity at the spanwise boundaries, traction-free conditions at cross-stream boundaries, and convective outflow condition at the outlet boundary. The time step is set as ∆tU ref /D = 9 × 10 -3 and 4.5×10 -3 on mesh M1 and M2 respectively, yielding a maximal cell CFL number of 0.41 during one shedding period. The initial condition is defined as the rest with random perturbations in the spanwise and crosswise direction to trigger the flow instability. The simulations were advanced in time up to a statistically steady flow pattern. Then, statistics were calculated by averaging in space in the spanwise direction and in time over 15 shedding cycles (T ∼ 80D/U ref ). In addition to LS-STAG computations on M1 and M2 meshes, a staircase simulation (see Section 6.1) on M1 mesh has also been performed in order to evaluate the benefits of our cut-cell treatment. For both methods, the time history of the drag and lift coefficients is displayed in Figure 11. Although both methods give force signals with comparable amplitudes, the crude staircase treatment of the complex geometry introduces a spurious lag in the force signals which yield incorrect prediction of the shedding frequency, as we shall see next. The flow at Re D = 300 is characterized by the presence of streamwise-oriented coherent structures of spanwise length order of D [START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF] in the cylinder wake. These large vortical structures are best visualized by the Q-criterion [START_REF] Hunt | Eddies, streams, and convergence zones in turbulent flows[END_REF], which identifies vortical flow regions to positive values of the second invariant of the velocity gradient Q. Figure 12 displays isosurfaces of Q in the near-wake of circular cylinder at an instantaneous time. It clearly shows the presence of large counter-rotating streamwise vortices in the near cylinder wake, whose number decreases further downstream. A similar qualitative observation of the wake structure has been performed in [START_REF] Kravchenko | B-spline method and zonal grids for simulations of complex turbulent flows[END_REF] by visualizing the instantaneous streamwise vorticity. [START_REF] Mittal | On the inclusion of three-dimensional effects in simulations of twodimensional bluff body wake flows[END_REF][START_REF] Kravchenko | B-spline method and zonal grids for simulations of complex turbulent flows[END_REF], experimental measurements [START_REF] Wieselsberger | New data on the law of hydro and aerodynamic resistance[END_REF][START_REF] Norberg | Flow around a circular cylinder: Apects of fluctuating lift[END_REF], IBM computations [START_REF] Balaras | Modeling complex boundaries using an external force field on fixed cartesian grids in large-eddy simulations[END_REF] and the recent DNS of Rajani et al. [START_REF] Rajani | Numerical simulation of laminar flow past a circular cylinder[END_REF] with a finite-volume solver on multiblock polar grids. It can be clearly seen that the shedding frequency is the easiest quantity to predict: only the IBM and staircase method on the coarsest mesh M1 and Ref. [START_REF] Rajani | Numerical simulation of laminar flow past a circular cylinder[END_REF] misestimate its value. For the drag coefficient, accurate reference computations bracket its value in the range [1.26, 1.28], for which the LS-STAG computations on both coarse and fine meshes give a very close prediction, C D = 1.29. The rms lift coefficient proves to be the most sensitive value to predict, with reference values spread in the range [0.380, 0.435], with an upper value given by the experiments [START_REF] Norberg | Flow around a circular cylinder: Apects of fluctuating lift[END_REF]. Both coarse and fine mesh LS-STAG simulation fall in this range, and in particular the M2 prediction is the closest computational result to the experimental value. The mean velocity profiles, time-averaged velocity fluctuations and Reynolds shear stress at three different downstream locations in the wake of circular cylinder are displayed in Figures 13 to [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF], respectively. The LS-STAG and staircase computations are compared to IBM simulations [START_REF] Balaras | Modeling complex boundaries using an external force field on fixed cartesian grids in large-eddy simulations[END_REF] and reference spectral results [START_REF] Mittal | On the inclusion of three-dimensional effects in simulations of twodimensional bluff body wake flows[END_REF] (The B-spline solution of Ref. [START_REF] Kravchenko | B-spline method and zonal grids for simulations of complex turbulent flows[END_REF] are virtually indistinguishable from the spectral data and not shown in these figures). Firstly, it can be observed that the LS-STAG results on both coarse and fine meshes show excellent agreement with the reference results, which indicates that the M1 grid resolution is sufficiently fine to capture the wake dynamics at this flow regime. For this grid resolution, the IBM computation exhibits notable discrepancies on the cross-flow velocities (Figure 13) and velocity fluctuations (Figure 14) in the cylinder wake. It was conjectured that these discrepancies were caused by insufficient grid resolution, see [START_REF] Balaras | Modeling complex boundaries using an external force field on fixed cartesian grids in large-eddy simulations[END_REF]. However, on the light of LS-STAG computations, these discrepancies are more likely due to the inferior treatment of complex geometries given by the IBM compared to a cut-cell method of formally same order of accuracy. It is also interesting to observe that the staircase method, which is formally first-order accurate according to Section 6.1, gives superior results to the IBM method at similar grid resolution, see Table 4 and Figures 13 and 14.

Flow of shear-thinning liquids in a 1:2 axisymmetric sudden expansion

Finally, the LS-STAG method is applied to the flow of Newtonian and shear-thinning fluids in an axisymmetric sudden expansion geometry of ratio ER = 2, for which rheological characterization, experimental flow measurements (nuclear magnetic resonance, NMR) and FLUENT c simulations have recently been performed by Rigal [START_REF]Comportement de fluides complexes sous écoulement : Approche expérimentale par résonance magnétique nucléaire et techniques optiques et simulations numeriques[END_REF]. Following this work, the geometry

w(r) = 3n + 1 n + 1 U ref 1 - 2r D 1 n+1 n , ( 36
)
where U ref is the inlet mean velocity and n is the power-law index. 3 types of fluids are considered: water (Newtonian) and two Xanthan solutions (of 0.1% and 0.3% mass concentration) whose shearthinning behavior has been fitted with the power-law (4) and Cross (5) models in Table 5. The power-law parameters are used to calculate the flow Reynolds number 33 and inlet velocity 36, while the Cross model is used to model the shear-thinning viscosity in the LS-STAG computations. First, a grid convergence study is conducted on 4 refined meshes whose properties are reported in Table 6. In the axial direction, the grid size h z is kept constant in the vicinity of the singular expansion (z ∈ [-5 mm, 15 mm]), and is relaxed upstream and downstream. In the other 2 directions, the grid is uniform of size h. For the Xanthan 0.1% flow at Re D1 = 56, Figure 17 displays axial velocity profiles at various cross sections (defined in Figure 16(right)) given by the LS-STAG computations and the experiments. Negative velocities at z = 1 mm and z = 9 mm indicate the presence of a recirculation zone which vanishes downstream (z = 21 mm). An excellent agreement is observed with the experimental results measured by NMR velocimetry, although a close inspection of the profile at z = 1 mm next to the abrupt expansion indicates a lack of grid resolution for the 2 coarsest meshes only.

- The main feature of expansion flows is the presence of a recirculating vortex caused by the separation of fluid flow at the singularity and its reattachment at the wall further downstream. For non-Newtonian fluids, the vortex characteristics in axisymmetric geometries have been studied experimentally in [START_REF] Halmos | The behavior of a power-law fluid flowing through a sudden expansion. Part ii. experimental verification[END_REF][START_REF] Pak | Separation and reattachment of non-newtonian fluid flows in a sudden expansion pipe[END_REF] and numerically in [START_REF] Halmos | The behavior of a power-law fluid flowing through a sudden expansion: Part i. a numerical solution[END_REF][START_REF] Nag | Variation of the recirculation length of newtonian and non-newtonian power-law fluids in laminar flow through a suddenly expanded axisymmetric geometry[END_REF]. The quantities of interest are the vortex length L v , its intensity ∆ψ characterized by the maximal value of the stream function in the recirculation zone, and the establishment length L e which is the downstream distance to the expansion plane from which the axial velocity becomes constant. These quantities, non-dimensionalized with the step-height d = (D 2 -D 1 )/2, are reported in Table 7 at Re D1 = 56 for meshes M1-M4. This table shows that nearly grid independent results are obtained on grid M3, with discrepancies below 1% for L v and L e , and below 3% for ∆ψ compared to the M4 results. It can be observed that the higher is the shear-thinning index, the weaker is the intensity of vortex. In effect, since the drop in velocity gradients is more pronounced in the recirculation zone for high shear-thinning fluids, the bulk flow is able to drive the secondary flow at a lower rate. Furthermore, we observe that the establishment length increases significantly with the shear-thinning index, while the vortex size tends to shorten. A similar trend has also been observed in the numerical computations of [START_REF] Halmos | The behavior of a power-law fluid flowing through a sudden expansion: Part i. a numerical solution[END_REF] for power-law fluids and the experiments of [START_REF] Pak | Separation and reattachment of non-newtonian fluid flows in a sudden expansion pipe[END_REF] for purely viscous non-Newtonian fluids in the laminar regime.

2 -1 0 1 2 r/R 1 0 0.5 1 1.5 w/U ref z = -9 mm z = 1 mm z = 9 mm z = 21 mm
More insight in the vortex characteristics is given by varying the Reynolds number, with computations performed at Re D1 = 9, 22, 56, 77.5, 100, 160 and 217 on the M3 mesh. Figure 18 shows the vortex length predicted by the LS-STAG code and the axisymmetric FLUENT c computations of Rigal [START_REF]Comportement de fluides complexes sous écoulement : Approche expérimentale par résonance magnétique nucléaire et techniques optiques et simulations numeriques[END_REF] with respect to the generalized Reynolds number [START_REF] Nag | Variation of the recirculation length of newtonian and non-newtonian power-law fluids in laminar flow through a suddenly expanded axisymmetric geometry[END_REF] Re gen = 2 3-n ( n 3n+1 ) n Re D1 , where n is the power-law index given in Table 5. As first observed by Halmos et al. [START_REF] Halmos | The behavior of a power-law fluid flowing through a sudden expansion: Part i. a numerical solution[END_REF][START_REF] Halmos | The behavior of a power-law fluid flowing through a sudden expansion. Part ii. experimental verification[END_REF], the length of recirculation zone is a linear function of Re gen for all fluid types. Deviation from linearity for the FLUENT c computations of the Xanthan 0.3% flows may be attributed to grid under-resolution [START_REF]Comportement de fluides complexes sous écoulement : Approche expérimentale par résonance magnétique nucléaire et techniques optiques et simulations numeriques[END_REF]. In addition, the empirical curve fitting of Nag & Datta [START_REF] Nag | Variation of the recirculation length of newtonian and non-newtonian power-law fluids in laminar flow through a suddenly expanded axisymmetric geometry[END_REF] have been added to this figure. For expansion ratio ER = 2, it predicts a linear growth of the recirculation zone as: Figure 19 displays the evolution of establishment length L e and secondary vortex intensity ∆ψ with respect to the Reynolds number for the various fluids. The monotonic growth of L e with the Reynolds number is directly related to the size of vortex. In effect, a longer recirculation zone allows the flow to decelerate more gradually, retarding thus flow establishment as the Reynolds number increases. Vortex intensity is a monotically increasing function of the Reynolds number that reaches a plateau at high Reynolds number. In [START_REF] Halmos | The behavior of a power-law fluid flowing through a sudden expansion: Part i. a numerical solution[END_REF], this behaviour is attributed to the increase in rate of momentum transfer from the bulk flow to the recirculation zone, until the momentum dissipated in the wall region reaches an asymptotical value. Similar trends were observed in [START_REF] Scott | A finite element analysis of laminar flows through planar and axisymmetric abrupt expansions[END_REF][START_REF] Badekas | Eddy correlations for laminar axisymmetric sudden expansion flows[END_REF] for Newtonian flows with various expansion ratios.

L v d =

Concluding remarks

This paper has presented the extension of the LS-STAG cut-cell method to 3D extruded geometries in a parallel computing environment. A series of benchmark Newtonian and non-Newtonian flows have been conducted for assessing the performance of the LS-STAG solver (accuracy, computer resources) compared to well established flow solvers (Ansys-FLUENT c on body-fitted grids, immersed boundary method). To our knowledge, these issues have seldomly been addressed in the literature to such extent.

The results of these benchmark cases are very promising for cut-cell methods. The computations of the confined cylinder benchmark (Sect. 6.2.1) showed that a cut-cell solver was able to provide solution nearly as accurate as a commercial body-fitted CFD package, at a fraction of the computational cost. The DNS of turbulent cylinder flow (Sect. 6.2.2) highlighted the superior accuracy of the LS-STAG treatment of complex geometries compared to more traditional immersed boundary methods. Finally, the computation of Xanthan flows in axisymmetric sudden expansion illustrates the versatility of the method.

However, a close inspection on the accuracy of the LS-STAG method shows that the convergence rate is not as smooth as for methods on body-fitted grids. As stated in Sect. 6.1 and 6.2.1, this behaviour is certainly caused by the ad-hoc discretization of the off-diagonal velocity gradients ∂w/∂x and ∂w/∂y in cut-cells, which is not yet as mature as the discretization of other terms in the Navier-Stokes equations. Future work will be dedicated to the improvement of the discretization of these terms, by adding non-orthogonal corrections in the fashion of the diamond scheme for unstructured grids [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection-diffusion problem[END_REF][START_REF] Moukalled | The finite volume method in computational fluid dynamics[END_REF].

Figure 1 :

 1 Figure 1: Four basic types of fluid cells present in the LS-STAG mesh for extruded geometries. Velocity unknowns are located at face centroids, pressure and diagonal components of ∇v at cell centroids (•), off-diagonal components of ∇v at edge midpoints : ∂u ∂y , ∂v ∂x ( ), ∂w ∂x , ∂u ∂z ( ) and ∂w ∂y , ∂v ∂z ( ).

Figure 2 :

 2 Figure 2: Geometry parameters of pentagonal cut-cell Ω i,j,k . The top elementary fluid face Γ t i,j,k is colored in green, the immersed boundary face Γ ib i,j,k is outlined with red dashed line (--) and the level-set values φ i,j,k are stored at cell vertices ( ).

Figure 4 :

 4 Figure 4: Information interchange between two adjacent sud-domains (I, J, K) and (I + 1, J, K). Data in the buffer cells Ω Nx+1,j,k and Ω 0,j,k (in red) are duplicated respectively from the appropriate values in the computational cells Ω 1,j,k and Ω Nx,j,k (in blue).

Figure 5 :

 5 Figure 5: Flowchart of the parallel fractional-step algorithm, the corresponding sequential algorithm is recovered by discarding steps in shaded boxes. The steps marked in italic are required for simulations of non-Newtonian flows only.

Figure 6 :

 6 Figure 6: Maximal error on the streamwise velocity versus grid size h for the Hagen-Poiseuille flow. At left: on 90% of the fluid domain; at right: on the whole fluid domain.

Figure 7 :

 7 Figure 7: Pointwise error profile for h = 5.5 × 10 -3 on cross-section z = 1 of Hagen-Poiseuille flow.

Figure 8 :

 8 Figure 8: Computational domain for the cylinder flows of Section 6.2.

Figure 9 :

 9 Figure 9: 2D view of the coarsest meshes for the circular cylinder benchmark: (a) FLUENT c L2 mesh and (b) LS-STAG M2 mesh.

Figure 10 :

 10 Figure 10: Time evolution of the force coefficients computed by LS-STAG, FLUENT c , and reference solution [22] at the cylinder. At left: on the total time interval; at right: close up around t = 4 s.

Figure 11 :

 11 Figure 11: Time evolution of the force coefficients for the cylinder flow at Re D = 300 on the M1 mesh.

Figure 12 :

 12 Figure 12: Vortex visualization of the cylinder flow at Re D = 300 by the Q-criterion on mesh M2. As in [10], the isosurfaces displayed correspond to the value Q = 0.1 and are colored with streamwise vorticity values: red is ωxD/U ref = 0.25 and blue is ωxD/U ref = -0.25.

Figure 13 :

 13 Figure 13: Mean velocity profiles at different locations in the wake of circular cylinder at Re D = 300 : stream-wise (left) and cross-flow (right) components.

Figure 14 :Figure 15 :

 1415 Figure 14: time-averaged velocity fluctuations at different locations in the wake of circular cylinder at Re D = 300. At left: stream-wise component and right: cross-flow component.

Figure 16 :

 16 Figure 16: Sketch of the axisymmetric sudden expansion geometry: At left: isometric view and at right: 2D view. The sudden expansion plane is placed at z = 0.

Figure 17 :

 17 Figure17: Grid convergence of axial velocity profiles at various cross sections before and after the expansion plane for the Xanthan 0.1% flow at Re D 1 = 56. Circles (•) correspond to experimental results[START_REF]Comportement de fluides complexes sous écoulement : Approche expérimentale par résonance magnétique nucléaire et techniques optiques et simulations numeriques[END_REF] while the dash-dotted (-•-), dotted (••), dashed (--) and solid (-) lines represent the LS-STAG solutions on meshes M1 to M4 respectively.

Figure 18 :

 18 Figure 18: Variation of the vortex length with respect to the generalized Reynolds number.

Figure 19 :

 19 Figure 19: Evolution of vortex intensity (left) and establishment length (right) with the Reynolds number: (--) LS-STAG and (-• * -•) FLUENT c [11] results.

  Divergence-free velocity field U n and pressure field P n available Compute γn at cell centers Communicate the values γn between neighboring processors

		U n , P n , γn
		Prediction Step :
	P n	Solve equation for the provisional velocity field	Ũ
	U n , P n+1	Communicate the values Ũ between neighboring processors Ũ , P n
		Projection step I :
	U n+1	Solve Poisson equation for P n+1
		Communicate the values P n+1 between neighboring processors
		Ũ , P n+1
		Projection step II :
		Update velocity field
		Communicate the values U n+1 between neighboring processors

Table 1 :

 1

				Cartesian Cut-cell Solid
	M2	0.016	83 × 38 × 16	97.4%	1.1%	1.5%	49136	0.118 × 0.105
	M3	0.008	159 × 77 × 32	97.6%	0.6%	1.8%	391776	0.060 × 0.051
	M4	0.004	317 × 153 × 64	97.8%	0.3%	1.9%	3104064	0.030 × 0.026

Salient properties of the LS-STAG meshes used for the circular cylinder benchmark. Mesh ∆tU ref /D Nx × Ny × Nz Cell Types # fluid cells hx/D × hy/D

Table 2 :

 2 Relative errors (in percent) of the force coefficients with respect to the reference solution[22] (C Dmax = 3.2978, C Lmax = 0.00276 and C L min = -0.0110) and L 2 errors for t ∈ [0, 8s]. grid convergence as smooth as FLUENT c . As already observed in Sect. 6.1, this behavior is certainly due to the ad hoc discretization of the off-diagonal gradients ∂w/∂x and ∂w/∂y in the cut-cells.

	Solver	Mesh	% Error C Dmax C Lmax C Lmin	% L 2 Error C D C L
		M2	3.72	11.05 36.61 3.74 48.06
	LS-STAG	M3	2.23	2.96	18.07 2.42 16.94
		M4	2.96	0.49	10.02 3.05 10.22
		L2	2.35	4.40	16.29 2.11 17.03
	FLUENT c	L3	0.96	2.16	10.96 0.81 10.01
		L4	0.25	0.08	2.89	0.21 2.67

errors on the extrema of force coefficients and L 2 errors for t ∈ [0, 8s] are compared. The LS-STAG results display a monotonic convergence for the lift coefficient, C Lmin , C Lmax and C L . The errors are comparable to FLUENT c on the first 2 meshes, and only slightly inferior to FLUENT c on M4/L4. In addition, results on the drag coefficient, C Dmax and C D , although satisfactory, do not display a

Table 3 :

 3 Clock time and memory usage of the LS-STAG and FLUENT c computations on the IJL-LEMTA parallel machine (np refers to the number of processors).

			Wall Time (min)	Memory Usage (mb)
	Mesh	np	LS-STAG	FLUENT c PISO NITA	LS-STAG	FLUENT c PISO NITA
	M2/L2	2	16	363	55	90	406	400
	M3/L3	4	132	2943	488	411	1368	1388
	M4/L4	4	2164	-	7841	2417	-	7927

Table 4 :

 4 Mean flow parameters for the cylinder flow at Re D = 300.

	Case	C D C Lrms	St
	Measurement [26]	1.22	-	-
	Measurement [24]	-	0.435 0.203
	Spectral simulation [8]	1.26 0.380 0.203
	B-spline simulation [9]	1.28	0.40	0.203
	LS-STAG, M1	1.29 0.420 0.203
	LS-STAG, M2	1.29 0.437 0.203
	Staircase, M1	1.24 0.346 0.217
	IBM Case 1 [10]	1.41	0.44	0.20
	IBM Case 2 [10]	1.27	0.42	0.21
	Body-fitted finite-volume [27] 1.28 0.499 0.195

Table 4

 4 compares the main flow characteristics, namely the mean drag coefficient, C D , root mean square average of the lift coefficient, C Lrms , and the Strouhal number St = f D/U ref , where f is the shedding frequency, with reference body-fitted simulations

Table 5 :

 5 Parameters of the power-law and Cross model reported byRigal [11] for the 2 Xanthan solutions.

	Fluid type	n	k (Pa.s n ) η 0 (Pa.s) η ∞ (Pa.s)	m	K (s)
	Xanthane 0.1% 0.543	0.128	0.376	0.00193	0.631	1.64
	Xanthane 0.3% 0.322	1.11	9.84	0.0048	0.764	13.1

Table 6 :

 6 Main properties of the grids for axisymmetric expansion flow.

	Mesh	N x × N y × N z	Cell Types Cartesian Cut-cell	Solid	h/D 1	h z /D 1
	M1	50 × 50 × 185	254844	32336	175320	4.08 × 10 -2 16.7 × 10 -2
	M2	100 × 100 × 370	2147480	130200	1422320 2.02 × 10 -2 8.33 × 10 -2
	M3	150 × 150 × 555	7321320	294780	4871400 1.35 × 10 -2 5.56 × 10 -2
	M4	200 × 200 × 740 17428564	523916 11647520 1.01 × 10 -2 4.17 × 10 -2

Table 7 :

 7 Grid convergence study of the dimensionless vortex length (Lv/d), vortex intensity ∆ψ * = ∆ψ/(πd 2 U ref ) and establishment length (Le/d) at Re D 1 = 56. The establishment length is defined as the location where the axial and outlet velocity difference becomes lower than 2%. type Mesh L v /d ∆ψ * × 10 3 L e /d of our computational results is excellent for the lower Reynolds number, while a deviation from this fit is observed as the Reynolds number is increased, most notably for the Xanthan 0.3% computations at Re D1 ≥ 56 (which corresponds to Re gen 200.1).

	Fluid water	M1 M2	5.05 5.00	15.86 15.32	14.33 14.17
		M3	4.97	15.23	14.11
		M4	4.96	15.11	14.08
		M1	4.25	6.24	15.33
	Xanthan 0.1%	M2	4.10	5.68	15.33
		M3	4.04	5.55	15.56
		M4	4.00	5.43	15.58
		M1	4.37	2.98	21.00
	Xanthan 0.3%	M2	4.16	2.44	21.00
		M3	4.08	2.27	21.22
		M4	4.03	2.16	21.33
	why the agreement

Preprint submitted to Journal of L A T E X TemplatesDecember 20, 2017


B Portelenelle

Olivier Botella

Y Cheny

ACCURATE DISCRETIZATION OF DIFFUSION IN THE LS-STAG CUT-CELL METHOD USING DIAMOND CELL TECHNIQUES

Keywords: Computational fluid dynamics (CFD), Incompressible Flows, Heat Transfer, Cut-Cell Method, Gradient Discretisation, Data Interpolation

HAL is

INTRODUCTION

The LS-STAG method [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF] is a Cartesian method for incompressible flow computations in irregular geometries which aims at discretizing accurately the flow equations in the cut-cells, i.e. cells of complex polyhedral shape formed by the intersection of the Cartesian mesh with the immersed boundary. In the three basic types of 2D cut-cells (pentagonal, trapezoidal and triangular, see Fig. 1(a)), the discretization is designed such as the global invariants of the flow (total mass, momentum and kinetic energy) are preserved at the discrete level. The LS-STAG discretization in the cut-cells is consistent with the MAC discretization used in Cartesian fluid cells, and has the ability to preserve the 5-point Cartesian structure of the stencil, resulting in a highly computationally efficient method. We have successfully applied the LS-STAG method to Newtonian flows at moderate Reynolds number in fixed and moving geometries [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF], pseudoplastic flows [START_REF] Botella | The LS-STAG immersed boundary method for non-Newtonian flows in irregular geometries : Flow of shear-thinning liquids between eccentric rotating cylinders[END_REF], and viscoelastic flows [START_REF] Botella | Application of the LS-STAG immersed boundary/cut-cell method to viscoelastic flow computations[END_REF]. In a recent paper [START_REF] Nikfarjam | The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries[END_REF], we have extended the LS-STAG discretization for 3D configurations with translational symmetry in the z direction (subsequently called 3D extruded configurations), where are only present the 4 types of cells depicted in Fig. 1(b). In such geometries, only extruded 2D cut-cells are present, which enables us to easily extend the principles of the energy-conserving discretization of Ref. [START_REF] Cheny | The LS-STAG method : A new immersed boundary / level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF] to the cutcells of Fig. 1(b) : the discretization of the continuity equation, pressure gradient, convective fluxes, normal stresses are straightforward extensions of the 2D LS-STAG formulas.

However, the discretization of cell-face fluxes involving shear stresses such as ∂w/∂x and ∂w/∂y, which were absent from the 2D case, needed further attention. In effect, due to the non-orthogonality of the cut-cells, the use of 2-point formulas for computing these fluxes proves to be inaccurate. A way to improve the accuracy is to compute the whole velocity gradient at the cut-cell faces, thus decomposing the flux as an orthogonal contribution (using a standard 2-point formula) and non-orthogonal correction (using data at cell vertices). The data at cutcell vertices are then interpolated from face-normal velocities and boundary conditions. This gradient reconstruction technique is commonly denominated "secondary gradients" [START_REF] Mathur | A pressure-based method for unstructured meshes[END_REF] in the CFD community and "diamond cell method" [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection-diffusion problem[END_REF] in the applied mathematics community.

This issue is also found for computing heat transfer problems with the 2D version of the LS-STAG method. Hence, for the ease of discussion the diamond cell method will first be described and analyzed for the diffusion fluxes of the heat equation, using various interpolation For the 3D extruded mesh, the velocities are discretized at cell faces, pressure and normal stresses at cell center (•), shear stresses at edge centers : ∂u ∂y , ∂v ∂x ( ), ∂u ∂z , ∂w ∂x ( ) and ∂v ∂z , ∂w ∂y ( ). For heat transfer problems, the temperature is located at cell centers.

schemes (inverse distance weighting, least-squares, barycentric interpolation on triangulation). The accuracy of the discretization will firmly be assessed on a series of benchmark problems (Taylor-Couette flow, natural convection from a cylinder in an enclosure [START_REF] Demirdžić | Fluid flow and heat transfer test problems for non-orthogonal grids: Bench-mark solutions[END_REF]) by inspecting the formal order of accuracy and the heat flux distribution at the immersed boundary. Finally, the diamond cell technique is employed for enhancing the accuracy of the shear stresses in 3D extruded geometries. Comparison with body-fitted CFD codes will be provided in terms of accuracy for benchmark flow past circular cylinder [START_REF] Bayraktar | Benchmark computations of 3D laminar flow around a cylinder with CFX, OpenFOAM and FeatFlow[END_REF].

DISCRETIZATION OF DIFFUSION WITH THE DIAMOND CELL TECH-NIQUES

The various discretization of diffusion that we consider are best detailed for the 2D temperature equation :

∂T ∂t + ∇ • (vT ) -κ∇ 2 T = 0, ( 1
)
where T is the temperature, v = (u, v) is the convection velocity and κ is the thermal diffusivity. In cut-cell Ω i,j depicted in Fig. 2, whose volume is denoted V i,j , the temperature T i,j is located at the centroid x G i,j = (x G i,j , y G i,j ), and the faces of the cut-cell are decomposed as

Γ i,j = Γ w i,j ∪ Γ e i,j ∪ Γ s i,j ∪ Γ n i,j ∪ Γ ib i,j
, where Γ ib i,j denotes the solid part of the cut-cell, while the fluid faces are denoted with the standard compass notations. The diffusive flux at the east face F d e = Γ e i,j κ∇T • n e dS is discretized with midpoint quadrature, yielding :

F d e ∼ = θ u i,j ∆y j κ ∂T ∂x e , ( 2
)
where θ u i,j ∈ [0, 1] is the fluid part of the face, such as θ u i,j ∆y j = x N -x S in Fig. 2, and ∂T /∂x| e = ∇T | e • n e is the face-normal temperature gradient.

Several alternate formulas can be considered for computing ∂T /∂x| e . The simplest is the following two-point approximation :

∂T ∂x e ∼ = T i+1,j -T i,j d , ( 3
)
where d = x G i+1,j -x G i,j is the vector joining the cell centroids. In the case where d and the face normal vector n e = e x are colinear, formulae (2)-(3) recover the usual second-order Cartesian discretization. However, for the cut-cell depicted in Fig. 2, the colinearity is lost and formula (3) loses its accuracy as the angle α = mes(e x , d) becomes large. This issue, which is analogous to the non-orthogonality of curvilinear and unstructured grids [START_REF] Moukalled | The Finite Volume Method in Computational Fluid Dynamics[END_REF], has first been observed in [START_REF] Van Der Plas | Efficiently simulating viscous flow effects by means of regularization turbulence modeling and local grid refinement[END_REF] for the LS-STAG method in 3D extruded geometries. In this reference, they proposed to compute off-diagonal components of the velocity gradients of the Navier-Stokes equations with a formula consistent with the pressure gradient discretization. For the case of the temperature equation, this ad hoc formula amounts to :

∂T ∂x e ∼ = T i+1,j -T i,j ( 1 2 V i+1,j + 1 2 V i,j )/θ u i,j ∆y j . ( 4
)
However, the Navier-Stokes computations in [START_REF] Nikfarjam | The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries[END_REF] reported that the use of the above two-point formulas diminishes the accuracy of the LS-STAG method, and that only a superlinear order of convergence was obtained. The reason is that, when α = 0, the face-normal gradient ∂T /∂x| e cannot be written as a function of T i,j and T i+1,j only, it also has a component tangential to the face also. It is thus necessary to compute the whole temperature gradient, which reads in the local basis of the east face :

∇T | e = ∂T ∂x e e x + ∂T ∂y e e y , (5)
since n e = e x in Fig. 2. An elegant way to proceed is to approximate the gradient in the d direction as :

∇T | e • d d = cos α ∂T ∂x e + sin α ∂T ∂y e ∼ = T i+1,j -T i,j d , ( 6
)
while the gradient tangential to the face is computed as : where T N and T S are located at the face endpoints x N and x S . Then, substitution of ( 7) in ( 6) gives the following approximation of the face-normal gradient :

∂T ∂y e ∼ = T N -T S θ u i,j ∆y j , ( 7
)
∂T ∂x e ∼ = T i+1,j -T i,j d cos α - T N -T S θ u i,j ∆y j tan α, (8a)
with

d cos α = x G i+1,j -x G i,j , tan α = y G i+1,j -y G i,j x G i+1,j -x G i,j . ( 8b
)
The first quotient in the RHS of (8a) is the orthogonal contribution of the gradient discretization, while the second quotient is the non-orthogonal contribution that vanishes when d and n e are colinear (α = 0), yielding a two-point approximation identical to (3) and (4). A similar procedure is employed to approximate diffusion at the immersed boundary. In the case of homogeneous boundary condtions T = T ib , the non-orthogonal contribution vanishes and the face normal gradient to Γ ib i,j is :

∂T ∂n ib ∼ = T ib -T i,j ∆ , ( 9
)
where ∆ is the perpendicular distance from point x G i,j to segment Γ ib i,j . The gradient formula (8) can also be derived by applying the Green-Gauss theorem in the diamond-shaped control volume formed by joining the cell centers x G i,j , x G i+1,j and the face endpoints x N , x S , and is known as the diamond cell scheme. The diamond cell scheme has been introduced in [START_REF] Martinelli | Calculations of viscous flows with a multigrid method[END_REF][START_REF] Martinelli | Validation of a multigrid method for the Reynolds averaged equations[END_REF][START_REF] Jacon | A Navier-Stokes algorithm for turbulent flows using an unstructured grid and flux difference splitting[END_REF] for viscous compressible computations on unstructured meshes. Mathematical analysis has been performed in [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection-diffusion problem[END_REF][START_REF] Bertolazzi | On vertex reconstructions for cell-centered finite volume approximations of 2D anisotropic diffusion problems[END_REF] for diffusion equations. Furthermore, formula [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection-diffusion problem[END_REF] is equivalent to the over-relaxed approach for diffusion discretization presented in [START_REF] Moukalled | The Finite Volume Method in Computational Fluid Dynamics[END_REF].

PRACTICAL IMPLEMENTATION OF THE DIAMOND CELL SCHEME

In formula (8a), vertex value at x N is obtained from the boundary conditions (T N = T ib ), while the value at x S need to be reconstructed from neighboring cell-centered values :

T S ∼ = 1 k,l=0 ω i+k,j-1+l T i+k,j-1+l . ( 10
)
Bertolazzi & Manzini [START_REF] Bertolazzi | On vertex reconstructions for cell-centered finite volume approximations of 2D anisotropic diffusion problems[END_REF] reviewed various possibilities for defining the reconstruction weights ω i,j . Basic choices includes the inverse-volume interpolation, with weights based on the cell volumes :

ω i,j = V i,j 1 k,l=0 V i+k,j-1+l , ( 11
)
and the inverse-distance interpolation, with weights :

ω i,j = x G i,j -x S -p 1 k,l=0 x G i+k,j-1+l -x S -p , ( 12
)
where p ≥ 1 is a given exponent (p = 1 corresponds to the Euclidian distance). These two choices of weights give exact reconstruction for constant functions only, reducing thus the formal accuracy of the DCM formula to first-order only. In [START_REF] Bertolazzi | On vertex reconstructions for cell-centered finite volume approximations of 2D anisotropic diffusion problems[END_REF], an algorithm is discussed to reach second-order accuracy, by altering the values of the weights in order to fulfill the condition for exact reconstruction of linear functions. Note however that this algorithm involves centroid selection, solution of linear systems, and the resulting set of weights is not unique.

Second-order accuracy can also be reached by using least-square weights that minimize the reconstruction of linear functions. This least-square reconstruction is described in [START_REF] Jacon | A Navier-Stokes algorithm for turbulent flows using an unstructured grid and flux difference splitting[END_REF] for the solution of the compressible Navier-Stokes equations on unstructured grids, and has been analyzed in [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection-diffusion problem[END_REF][START_REF] Bertolazzi | On vertex reconstructions for cell-centered finite volume approximations of 2D anisotropic diffusion problems[END_REF] for model diffusion problems. However, our numerical experiments (see next section) showed that negative weights may appear for some cut-cells, giving rise to oscillatory solutions. For this reason, we have implemented a reconstruction based on triangulation [START_REF] Amidror | Scattered data interpolation methods for electronic imaging systems: A survey[END_REF] that ensures positivity of the weights. Out of the 4 triangles formed by the neighboring centroids of a given vertex, we select the triangle of largest minimal angle that contains the vertex. Value at the vertex is then obtained by barycentric interpolation.

In the LS-STAG code, diffusion is treated with implicit time-integration method, so the use of the DCM with vertex reconstruction has the effect to increase the stencil of the linear systems. A shown in Fig. 3, the DCM discretization needs to be applied to the four fluid faces of cut-cell Ω i,j , yielding a 9-point stencil instead of the standard 5-point stencil for 2-point discretization. Note that DCM needs only to be applied to neighboring cells of the immersed boundary. In the vast majority of the computational domain, mesh is orthogonal and standard 5-point stencil is recovered. The linear systems are solved with the HYPRE library (V2.10) using the semi-structured-grid system interface [START_REF] Falgout | Center for Applied Scientific Computing (CASC)[END_REF].

NUMERICAL RESULTS

The diamond cell method (8) (labeled herafter DCM) is implemented in the 2D and 3Dextruded LS-STAG code for the discretization of face-normal gradients in the temperature and Navier-Stokes equations. First, results are presented for the 2D temperature equation with no coupling with the Navier-Stokes equations. For interpolation at cell vertex, the following reconstruction schemes are compared :

-LS-STAG/DCM-ID1: first-order inverse distance weights [START_REF] Martinelli | Calculations of viscous flows with a multigrid method[END_REF] with p = 2, -LS-STAG/DCM-ID2: second-order inverse distance weights (with p = 2) using the algorithm discussed in [START_REF] Bertolazzi | On vertex reconstructions for cell-centered finite volume approximations of 2D anisotropic diffusion problems[END_REF],

-LS-STAG/DCM-LSQ: least-square weights [START_REF] Jacon | A Navier-Stokes algorithm for turbulent flows using an unstructured grid and flux difference splitting[END_REF],

-LS-STAG/DCM-TRI: barycentric interpolation on triangulation.

The performance of DCM is also compared with the 2-point gradient (4), which is labeled LS-STAG/2PT. For all methods, Dirichlet boundary conditions are imposed with Eq. ( 9). Subsequently, the temperature equation is coupled with the Navier-Stokes equations for computing 2D natural convection problems. The purpose is to assess the ability to compute accurately local heat flux at immersed boundaries, which is a sensible quantity since it involves the temperature gradient at the sold face of cut-cells. Once the DCM has been assessed on 2D flows, the DCM is applied to the Navier-Stokes equations in 3D-extruded geometries to enhance the discretization of isothermal incompressible flows computed with the LS-STAG method with 2-point discretization reported in [START_REF] Nikfarjam | The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries[END_REF]. First the spatial accuracy of the various variants of DCM is assessed on the test case of pure thermal conduction (Eq (1) with v = 0) between concentric cylinders of diameters r = R 1 and r = R 2 , where a temperature difference ∆T = T 2 -T 1 is applied. The steady exact solution is purely radial :

TEMPERATURE CONDUCTION IN TAYLOR-COUETTE CELL

T ex (r) = ln (r/R 2 ) ln(R 1 /R 2 ) (T 1 -T 2 ) + T 2 . ( 13
)
For the numerical experiments, we set ∆T = 1K and employ the same geometry (R 1 = 1 m and R 2 = 4 m) and uniform meshes (of size ∆x = ∆y = h) considered in [START_REF] Botella | The LS-STAG immersed boundary method for non-Newtonian flows in irregular geometries : Flow of shear-thinning liquids between eccentric rotating cylinders[END_REF] for non-Newtonian computations with the LS-STAG method.

For the various discretization of the face-normal temperature gradient, Fig. 4 compares the maximal error :

E ∞ h (T ) = max CVs Ω i,j T i,j -T ex (x G i,j ) max CVs Ω i,j T ex (x G i,j ) (14)
obtained on the whole computational domain and away from the immersed boundary. As expected, the 2-point method gives the worst accuracy, measured as first-order only. The errors given by all DCM variants are more than one decade lower, highlighting the benefits of this gradient computation. Note however that the errors of DCM-ID1 quickly saturates to a first-order rate when the grid size is reduced, due to the low accuracy of the vertex reconstruction scheme. All other reconstruction schemes (DCM-ID2, DCM-LSQ, DCM-TRI) gives comparable second-order accuracy. For these methods, it is worthwhile to note that the errors on 90% and 100% of the domain are of same magnitude, showing that the discretization in cut-and Cartesian cells has a comparable accuracy. The second test-case concerns a fluid flow and heat transfer problem in complex geometries, for which Demirdžić et al [START_REF] Demirdžić | Fluid flow and heat transfer test problems for non-orthogonal grids: Bench-mark solutions[END_REF] published accurate benchmark solutions. This test case concerns B. Portelenelle, O. Botella and Y. Cheny natural convection of a heated cylinder placed eccentrically in an inclosure (test C4 in [START_REF] Demirdžić | Fluid flow and heat transfer test problems for non-orthogonal grids: Bench-mark solutions[END_REF]), where the heat equation ( 1) is coupled with the Navier-Stokes equations in the Boussinesq approximation. The computational domain and boundary conditions are shown in Fig. 5 (left), with gravity acting in the downwards direction: since the solution possesses a vertical symmetry at x = 0, only half of the geometry is considered. As in [START_REF] Demirdžić | Fluid flow and heat transfer test problems for non-orthogonal grids: Bench-mark solutions[END_REF], the dimensions and fluid properties are: enclosure length L = 1 m, density ρ = 1 kg • m -3 , gravity constant g = 1 m • s -2 , thermal expansion coefficient β = 0.1K -1 , thermal diffusivity κ = 10 -4 m 2 • s -1 , dynamic viscosity µ = 10 -3 kg • m -1 • s -1 and temperature difference T H -T C = 1 K, that set the Rayleigh and Prandtl number respectively as:

NATURAL CONVECTION FROM A CYLINDER IN AN ENCLOSURE

Ra ≡ βρgL 3 µκ (T H -T C ) = 10 6 , ( 15
)
Pr ≡ µ ρκ = 10. ( 16
)
The purpose of these computations is to assess the ability of the LS-STAG/DCM method to compute accurately gradients at the immersed boundary, in the present case the average Nusselt number on the half-cylinder of length |Γ ib | :

Nu = 1 |Γ ib | (T H -T C )/L Γ ib ∇T • n dS, (17)
and the distribution of the local Nusselt number along the half cylinder, which is calculated at the center of the solid boundary of the cut-cells using formula [START_REF] Demirdžić | Fluid flow and heat transfer test problems for non-orthogonal grids: Bench-mark solutions[END_REF]. The LS-STAG computations are performed on a series of refined meshes M1-M5 of square uniform Cartesian cells, where the coarsest mesh M1 represented in Fig. 5 (right) has 32 × 64 cells and the finer meshes are obtained by doubling the number of cells in each direction, such as the finest mesh M5 has 512 × 1024 cells. For comparison purpose, we have also computed solutions with the commercial CFD package FLUENT c (V15.07), on a series of curvilinear meshes L1-L5 (see Fig. 5 (middle)) of same number of cells, with cell clustering near the heated cylinder. The FLUENT c computations have been performed with the pressure-based solver with the steady SIMPLE algorithm, with second-order spatial discretizations: QUICK scheme for convection discretization, Green-Gauss node-based interpolation for gradient calculation, and body-force weighted pressure interpolation. Fig. 6 shows the convergence of the average Nusselt number as the grid is refined. Irrespective of the gradient discretization scheme, all LS-STAG variants converge to the same grid-independent value at a similar rate. It is worthwhile to observe that on the 2 coarsest mesh, the FLUENT c computations are more accurate than LS-STAG. This is most certainly due to the high orthogonal quality of the curvilinear meshes near the circular cylinder. Fig. 7 displays the distribution of the local Nusselt number along the immersed cylinder obtained for the LS-STAG and FLUENT c computations on the M4/L4 mesh, along with the reference results [START_REF] Demirdžić | Fluid flow and heat transfer test problems for non-orthogonal grids: Bench-mark solutions[END_REF]. In order to quantify the influence of the cut-cells shape on the accuracy of the results, the LS-STAG profiles are colorized with the volume-change ratio of cut-cell Ω i,j :

RV i,j = min (k,l)∈N b V i,j V k,l , ( 18
)
where N b refers to the neighboring cells that share a face with Ω i,j (solid cells are excluded).

For each cut-cell RV i,j ∈ ]0, 1] since at least one neighboring cell is totally fluid, and a value close to 0 indicates a cut-cell of small size compared to the Cartesian cells. First it can be observed that the LS-STAG method with 2-point gradients gives a highly oscillatory profile, confirming the inaccuracy of this approximation of diffusion. The use of the diamond scheme reduces the level of oscillations for the first-order method DCM-ID1 and for DCM-LSQ. For both methods, the oscillations are clearly located in cut-cells with small volume-change ratio RV i,j , indicating inaccurate approximation in small cut-cells. In contrast, the profiles obtained with the other second-order methods (DCM-ID2 and TRI) are virtually free of any oscillations, and perfectly match the reference results.

As a conclusion, these computations show that accurate solutions can be achieved with the DCM discretization. Nonetheless, the accuracy of the DCM approximation is highly dependent on the interpolation scheme that is used, otherwise inaccurate results are encountered in the smaller cut-cells. These results show that the second-order methods DCM-ID2 and DCM-TRI stand as the best discretization methods, yielding indistinguishable results on the local Nusselt distribution. Owing to the simplicity of the triangular interpolation compared to the inverse distance interpolation of second order, the DCM-TRI variant is selected for implementation in the 3D-extruded code.

ISOTHERMAL FLOWS IN 3D EXTRUDED GEOMETRIES

The last test-case concerns the unsteady laminar flow past a circular cylinder in a confined square duct, for which accurate benchmark results with unstructured solvers are reported in [START_REF] Bayraktar | Benchmark computations of 3D laminar flow around a cylinder with CFX, OpenFOAM and FeatFlow[END_REF]. The flow unsteadiness is induced by a time-varying inflow boundary condition, yielding a timeperiodic Reynolds number Re ∈ [0, 100]. This benchmark flow has previously been performed in [START_REF] Nikfarjam | The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries[END_REF] with the LS-STAG method for 3D extruded geometries, using 2-point formulas for the diffusive fluxes of the Navier-Stokes equations.

A sketch of the velocity staggering in extruded cut-cells is shown in Fig. 8 : it is clearly seen that diamond cell discretization is needed for enhancing the discretization of for the off-diagonal components ∂w/∂x and ∂w/∂y of the diffusive fluxes of the Navier-Stokes equations. Due to the extrusion in the z direction, the other velocity gradients are accurately computed with 2-point formulas [START_REF] Nikfarjam | The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries[END_REF]. Note also that due to the staggering of the discrete velocities, the discretization of ∂w/∂x and ∂w/∂y involves co-planar unknowns (on planes of constant index k), such that only the 2D DCM formula ( 8) is needed.

Table 1 compares the results obtained on a series of refined meshes (500,000 to 3M cells) with an unstructured CFD solver (FLUENT c ), the LS-STAG solver with 2-point discretization [START_REF] Nikfarjam | The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries[END_REF] and the present DCM discretization. We refer to [START_REF] Nikfarjam | The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries[END_REF] for a complete description of the geometry, boundary conditions, meshes and discretizations used by both solvers. As observed in [START_REF] Nikfarjam | The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries[END_REF], the 2-point discretization of the LS-STAG solver gives only monotonic convergence for the lift quantities, while the drag quantities fails to convergence. In contrast, the DCM discretization show a dramatic increase in the accuracy of the LS-STAG discretization, yielding a smooth convergence for all force coefficients, with an accuracy comparable to the FLUENT c solver on all 3 grids.

Figure 1 :

 1 Figure 1: Cartesian cell and 3 basic types of cut-cells present in the LS-STAG mesh in 2D (a) and 3D extruded (b) geometries.For the 3D extruded mesh, the velocities are discretized at cell faces, pressure and normal stresses at cell center (•), shear stresses at edge centers : ∂u ∂y , ∂v ∂x ( ), ∂u ∂z , ∂w ∂x ( ) and ∂v ∂z , ∂w ∂y ( ). For heat transfer problems, the temperature is located at cell centers.

Figure 2 :

 2 Figure 2: Relevant notations used for the discretization of the diffusion flux at the east face of cut-cell Ω i,j .

Figure 3 :

 3 Figure 3: Stencil used for the discretization of diffusion in cut-cell Ω i,j . In black: 5-point stencil for 2-point discretization. For the DCM discretization with reconstruction at vertices ( ), unknowns at the green centroids are added to the stencil.

  /2PT LS-STAG/DCM-TRI LS-STAG/DCM-ID2 LS-STAG/DCM-ID1 LS-STAG/DCM-LSQ 1st order slope 2nd order slope

Figure 4 :

 4 Figure 4: Maximal error on the temperature versus grid size h for the pure conduction test case. At left : on 90% of the domain, at right : on the whole domain.

Figure 5 :

 5 Figure 5: Sketch of the computational domain for the natural convection test-case (with L = 1 m, T H = 1 K, T C = 0 K), curvilinear mesh L1 and Cartesian mesh M1 of size 32 × 64 used respectively for the FLUENT c and LS-STAG computations.

Figure 6 :

 6 Figure 6: Convergence of the average Nusselt number with respect to the number of fluid cells.

Figure 7 :

 7 Figure 7: Profiles of the Nusselt number along the half cylinder for the various discretization of diffusion in the LS-STAG method, and comparison with the FLUENT c computation and reference results of Demirdžić et al [9]. The profiles are colorized with respect to the volume-change ratio RV i,j (18) of the cut-cells.

Figure 8 :

 8 Figure 8: Staggered mesh for 3D extruded geometries in the z direction, and location of the off-diagonal components of ∇w.

Table 1 :

 1 Error C Dmax C Lmax C L min C DC Relative errors (in percent) of the force coefficients with respect to the reference solution[2] and L 2 errors for one-half period of the simulation (t ∈ [0, 8s]).

	Solver	Mesh		% Error	% L 2 06
		M3	2.23	2.96 18.07 2.42 16.94
		M4	2.96	0.49 10.02 3.05 10.22
		M2	1.31	6.99 34.00 1.48 41.91
	LS-STAG/DCM-TRI	M3	0.63	0.72 10.78 0.65 10.06
		M4	0.24	0.06	5.42 0.22 5.46
		L2	2.35	4.40 16.29 2.11 17.03
	FLUENT c	L3	0.96	2.16 10.96 0.81 10.01
		L4	0.25	0.08	2.89 0.21 2.67

L LS-STAG/2PT [16] M2 3.72 11.05 36.61 3.74 48.

Concluding remarks

This paper has presented a progress report on the application of the diamond cell method (DCM) for improving the accuracy of diffusion in cut-cell methods. The DCM is employed for cell-centred data (temperature in 2D) and face-centred data (Navier-Stokes equations in 3Dextruded geometries), which illustrates the versatility of the approach, that can be applied to other cut-cell methods than the LS-STAG method. One of the major findings is the dependence of the accuracy on the reconstruction scheme used for data interpolation. The use of barycentric interpolation allowed us to obtain local accuracy in cut-cells comparable to the one reached in Cartesian cells, irrespective of the size and shape of the cut-cell. Progress is on the way to employ DCM for full 3D computations with the LS-STAG method.