A Variational Multiscale Method with Multifractal Subgrid-scale Modeling for Large-eddy Simulation of Turbulent Flow

A variational multiscale method with multifractal subgrid-scale modeling is proposed for large-eddy simulation of turbulent flow. In the multifractal subgrid-scale model-ing approach, the subgrid-scale velocity is evaluated from a multifractal description of the subgrid-scale vorticity, which is based on the multifractal scale similarity of gradient fields in turbulent flow. The multifractal subgrid-scale modeling approach is integrated into a vari-ational multiscale formulation, demonstrating a new field of application of the variational multiscale concept. In addition, the application of the multifractal subgrid-scale modeling approach to wall-bounded turbulent flow is considered in this study. For this purpose, a near-wall limit of the multifractal subgrid-scale modeling approach is developed. The novel computational approach of multifractal subgrid-scale modeling within a variational multiscale formulation is then applied to turbulent flow over a backward-facing step. The results confirm a very good performance of the proposed method, and improved results are obtained compared to a dynamic Smagorinsky model and a residual-based variational multiscale method. Keywords: Large-eddy simulation, Multifractal subgrid-scale modeling, Variational multi-scale method, Wall-bounded turbulent flow, Backward-facing step.


INTRODUCTION
Turbulence is driven by vorticity and its self-induced velocity field acting on the vortical structures by stretching and folding them.The repeated stretching and folding of the vorticity field as well as the strain rate represents a multiplicative process causing multifractal scale similarity in gradient fields in turbulent flows such as kinetic energy dissipation and enstrophy, see, e.g., [16].This observation enables a novel approach to modeling turbulent flow.Multifractal subgrid-scale (MFS) modeling to estimate unresolved scales in Large-Eddy Simulation (LES, see, e.g., [15] for an overview) was presented in comprehensive form in [2] and tested for homogeneous isotropic turbulence in the accompanying study [3].It was proposed by those authors to estimate the subgrid-scale velocity field from the subgrid-scale vorticity field via the law of Biot-Savart, where the integrated subgrid-scale vorticity field is reconstructed by a multiplicative process.
The variational multiscale approach to LES was originally introduced in [9].Up to now, the variational multiscale method (VMM) has either been used based on a three-scale separation including a subgrid-viscosity term acting only on the smaller resolved scales as proposed in [9] and review in [5] or as a residual-based two-scale version using the residual together with an appropriate parameter to approximate the unresolved-scale quantities (see, e.g., [1]).An important aspect of the variational multiscale concept is that scales are separated by variational projection rather than by filtering.Particularly this feature offers a formulation that enables a straightforward introduction of structural subgrid-scale models into the projected Navier-Stokes equations.
In this study, we propose multifractal subgrid-scale modeling within a variational multiscale method, which represents the first approach to integrate multifractal subgrid-scale modeling into the variational multiscale framework; see [14].By introducing such a structural turbulence model into the variational multiscale framework, the variational multiscale concept for LES is broadened, and a new field of application is demonstrated.Important features of the method are scale separation by plain aggregation algebraic multigrid methods to further decompose the resolved scales, as proposed in [7] and used, e.g., in [6], and the development of an appropriate near-wall limit of the multifractal subgrid-scale modeling approach.The reader is referred to [14] for full methodical details, more elaborate discussion of the approach and further numerical examples.
The present study is organized as follows.In section 2, we present the variational multiscale formulation for LES and address the potential introduction of structural velocityestimation models.Section 3 is then devoted to the subgrid-scale modeling.The multifractal subgrid-scale modeling approach is derived, and its near-wall limit is addressed.The estimated subgrid-scale velocity is introduced into a residual-based form of the variational multiscale formulation.Multifractal subgrid-scale modeling within a variational multiscale method is then tested for the widely-used example of turbulent flow over a backward-facing step in section 4. Conclusions are drawn in section 5.

Variational formulation of incompressible Navier-Stokes equations
Fluid motion in the domain Ω described by the incompressible Navier-Stokes equations, given in convective form, where u denotes the velocity, ε (u) = 1 2 ∇u + (∇u) T the rate-of-deformation tensor, p the pressure, ν the kinematic viscosity and f a given body force vector, is considered.Appropriate initial and boundary conditions close the system.At t = 0, a divergence-free velocity field u 0 is prescribed.On the boundary Γ, Dirichlet and Neumann boundary conditions are given as where σ = −pI + 2νε(u) denotes the Cauchy stress tensor and n the outer unit normal vector on the boundary Γ.Moreover, it is assumed that Assuming appropriate solution function spaces S u for u and S p for p as well as weighting function spaces V u for the velocity weighting function v and V p for the pressure weighting function q, the system of equations ( 1)-( 2) is multiplied by v ∈ V u and q ∈ V p and integrated over the domain Ω.Viscous and pressure term are integrated by parts, with boundary conditions (3) and (4) applied to the resulting boundary integrals.The variational formulation of the incompressible Navier-Stokes equations is given as follows: ( The form on the left-hand side, comprising momentum and continuity part, is defined as The linear form (v), including the Neumann boundary condition, is given as Here, (•, •) = (•, •) Ω and (•, •) Γ N denote the usual L 2 -inner product on Ω and Γ N , respectively.

Variational multiscale formulation for large-eddy simulation
For the variational multiscale formulation of the Navier-Stokes equations, velocity and pressure are decomposed into resolved and unresolved (or subgrid) components as where resolved scales are identified by a spatial discretization of characteristic element length h.Subgrid scales are denoted by (•).According to the decomposition of the solution functions, direct sum decompositions of the underlying function spaces in the form S u = S h u ⊕ Ŝu and S p = S h p ⊕ Ŝp , respectively, are assumed.Based on the variational multiscale concept, we assume a variational projection for separating resolved and unresolved scales.A Galerkin finite element method may be interpreted as projection; see, e.g., [8] for a discussion of this issue.Therefore, we introduce a direct sum decomposition of the weighting function spaces Inserting the decomposition of velocity and pressure ( 8)-( 9) into the weak form (5), weighting separately by the resolved and the subgridscale part of the decomposed weighting functions and omitting the equation projected onto the space of unresolved scales, the variational multiscale formulation is obtained as: (10) where is the projection of the cross-stress terms and the projection of the Reynolds-stress term onto the space of resolved scales.The form contains the remaining linear terms in the unresolved-scale quantities.
The variational multiscale formulation (10) constitutes a mathematical framework for turbulence modeling in LES.In particular, it allows for directly modeling the crossand Reynolds-stress terms by structural velocity-estimation models such as the multifractal subgrid-scale modeling approach.

Multifractal subgrid-scale modeling
In the multifractal-subgrid scale modeling approach, the subgrid-scale velocity û is explicitly evaluated via the law of Biot-Savart Therein, the associated subgrid-scale vorticity field ω = ω e ω is recovered by a multifractal process consisting of two steps, which are outlined in this section.For an exhaustive presentation of the derivation, the reader is referred to [2] and, in the particular form for wall-bounded turbulent flow, to [14].First, the magnitude of the subgrid-scale vorticity field ω is derived by a multiplicative cascade distributing the subgrid-scale enstrophy Q = ω 2 within each element.Therefore, the subgrid-scale enstrophy Q is determined as a function of the enstrophy at the smaller resolved scales δQ h = δω h 2 , i.e., a scale range between the element length h and a larger scale αh, which is assumed to be located in the inertial range.The enstrophy spectrum in the inertial range scales as ∼ k 1 3 , where k is the wave number.Integrating the enstrophy spectrum both from the wave number k h to the viscous wave number k ν and from the smaller wave number k αh to k h enables a formulation for the subgrid-scale enstrophy depending on the enstrophy of the smaller resolved scales: A three-dimensional stochastic multiplicative cascade distributes the subgrid-scale enstrophy over each element.After a sufficient number of cascade steps, the resulting field becomes highly intermittent and displays multifractal scale-similarity; see, e.g., [4].The repeated application of a scale-invariant distribution of multipliers M n , mapping the subgridscale enstrophy in consecutive cascade steps to successively smaller subelements, leads to the following expression for the magnitude of the subgrid-scale vorticity: where N denotes the number of cascade steps.Second, the orientation of the subgrid-scale vorticity êω is derived from an additive decorrelation cascade.The orientation in the subgrid-scale vorticity field ω decorrelates at successively smaller scales from the orientation of the smaller resolved scales δe h ω as where f n denotes the decorrelation increments.As the strongest vortical structures maintain a preferred alignment with the local strain rate tensor over a relatively large range of scales, an intermittency factor I is defined from a correlation between ω and δω h as The orientation of the subgrid-scale vorticity can then be reformulated: where f * n are modified decorrelation increments.
After combining both cascades and introducing the resulting subgrid-scale vorticity ω, which is assumed approximately equal to its expectation value, into the law of Biot-Savart ( 14), the subgrid-scale velocity can be calculated as The ratio of the element length h to the viscous length scale λ ν determines the number of cascade steps N via The ratio of h to λ ν scales as ∼ Re h , where Re h is the local element Reynolds number.The local element Reynolds number as proposed in [14] is given by The necessary proportionality constant is set to 0.1, based on the experimental value given in [13].The subgrid-scale velocity û should become independent of Re h for Re h → ∞ and, consequently, N → ∞.This implies an expression for the intermittency factor I subject to N .Finally, the subgrid-scale velocity û reads as where The proportionality constant C sgs arises from the high Reynolds number limit of I and is set to 0.25 in this study.It is referred to [14] for a detailed discussion of C sgs .The required smaller resolved scales δu h are identified by further explicitly decomposing the resolved scales u h via level-transfer operators from plain algebraic multigrid methods, as already used for scale separation in, e.g., [6], where α = 3 was chosen.Eventually, the modeled forms of the crossand Reynolds-stress terms, ( 11) and ( 12), respectively, read For wall-bounded turbulent flow, a near-wall limit of the proposed multifractal subgridscale modeling approach was derived in [14].In the near-wall region, among other things, the vorticity field becomes highly anisotropic.The resulting stronger correlations in the orientation of the subgrid-scale vorticity cause an increase of the intermittency factor I. Higher intermittency factors are associated with an increase of the proportionality constant C sgs .Based on these considerations, C sgs is enhanced by an anisotropy factor C ai for wall-bounded turbulent flow: Assuming that (i) C nw sgs tends to a finite value as Re h → ∞ and C nw sgs → 0 as Re h → 1, (ii) the intermittency factor I is bounded as 0 ≤ I ≤ 1 and (iii) the strain rate represents a measure for anisotropy, the anisotropy factor C ai is defined as where is the element Reynolds number based on the norm of the strain rate tensor.

Residual-based multiscale modeling
As an introduction of artificial dissipation is not intended by the multifractal subgridscale modeling approach, its is incorporated into a residual-based form of the variational multiscale formulation to account for potential destabilizing effects induced by the numerical scheme.Appropriate stabilization terms are obtained by a residual-based approximation of subgrid-scale velocity and pressure, as shown, e.g., in [10].The final modeled variational multiscale formulation is then given as with the multifractal subgrid-scale modeling terms in the first line (second and third term) and the residual-based multiscale modeling terms in the second line (first three terms), which constitute, from left to right, a Streamline/Upwind Petrov-Galerkin (SUPG), a Pressure Stabilizing Petrov-Galerkin (PSPG) and a grad-div term, where R h M and R h C denote the discrete residuals of momentum and continuity equation, respectively.The stabilization parameters τ M and τ C as proposed in [17] and [18] are given by where ∂x j utilizes the coordinate system ξ of the element domain.The timestep size of the temporal discretization is denoted by ∆t, and C I is a positive constant, which is chosen to be 36.0for (tri-)linearly interpolated finite elements as used in this work.

TURBULENT FLOW OVER A BACKWARD-FACING STEP
The proposed method is investigated for turbulent flow over a backward-facing step.Based on the step height H and the mean centerline velocity U c at the inlet, the step Reynolds number Re H = U c H/ν is defined.The expansion ratio ER = 1.5 , which is the ratio of the channel height downstream and upstream of the step, is chosen according to the geometry of the experiment reported in [11].The step height H = 0.041m is also chosen similar to that study.The Reynolds number examined in that study was Re H = 5540.A Direct Numerical Simulation (DNS) of flow over a backward-facing step was presented in [12], where a step Reynolds number of 5100, which is close to the one in [11], and ER = 1.2 were used.Fig. 1 depicts the geometry of the step as well as the channel used to generate a turbulent inflow for the backward-facing step.
At the inlet of the backward-facing step, the velocity profile at a selected plane of the channel is prescribed as Dirichlet condition.No-slip boundary conditions are prescribed at the top and bottom wall, including the vertical step wall.Periodic boundary conditions are assumed in spanwise x 3 -direction.At the outlet, a zero-traction Neumann boundary condition (i.e., h = 0) is prescribed.In addition, a convective term is applied to ensure that potentially arising eddies at the outflow boundary are appropriately convected out of the domain; see, e.g., [14] for the respective formulation.At the inflow channel, no-slip boundary conditions are prescribed at the top and bottom wall.In homogeneous streamwise and spanwise directions, periodic boundary conditions are assumed.A constant pressure gradient in streamwise direction drives the flow through the channel.In [11], a friction Reynolds number Re τ = 290 was evaluated upstream of the step.Based on the step height, which equals the channel halfwidth, and the kinematic viscosity ν = 1.5268 • 10 −5 m 2 /s, the prescribed pressure gradient dp/dx 1 in the inflow section is 0.2844518N/m 3 .
A generalized-α time-integration scheme is applied.The parameter ρ ∞ is chosen to be 0.5, such that second-order accuracy in time is achieved.A constant time-step size ∆t = 0.0008s is used.After a sufficient number of time steps for turbulence to fully develop, statistics are collected during another 5000 time steps, representing approximately six flowthrough times.Two spatial discretizations are considered.The coarser discretizations contains 18 × 32 × 48 elements upstream of the step, 68 × 48 × 48 downstream and 32 × 32 × 48 in the channel in streamwise, wall-normal and spanwise direction, respectively.The finer discretization is obtained by doubling the number of elements.While the elements are uniformly distributed in spanwise direction, they are refined towards the step in streamwise direction and, in wall-normal direction, towards a horizontal line defined by the upper corner of the step as well as towards the upper and lower wall.
The results obtained with the proposed method are compared to results obtained with a dynamic Smagorinsky model (DSM) and the residual-based variational multiscale method (RBVMM) according to [1].Fig. 2 depicts the mean streamwise velocity at various locations behind the step.As differences between the results obtained with the coarser and the finer discretization are only of small amount, the mean streamwise velocity profiles of the simulations using the finer discretization are merely depicted.Differences between the models are hardly observable.All models provide results that are in good agreement with the experimental results.Small discrepancies between experiment and simulations occur around x 1 /H = 8 and further downstream.These deviations are attributed to the relative coarse mesh used in this region also for the finer discretization.
Figure 2. Backward-facing step: normalized mean streamwise velocity u 1 U c at various locations x 1 H using the finer discretization.
The mean reattachment length is determined by the location of zero wall-shear stress τ W or skin-friction coefficient C f = 2τ W /U c , respectively.Fig. 3 shows the skin-friction coefficient along the bottom wall.Results obtained with the coarser and finer discretization are marked by "(C)" and "(F)", respectively.DNS data from [12], named "DNS LMK97', are also included.Using the coarser discretization, MFS, RBVMM and DSM yield very similar results.Compared to DNS, they show a higher negative peak value.This peak is also somewhat shifted to the left.For the finer discretization, MFS captures the peak quite well, while RBVMM and DSM produce slightly too negative values.The predicted mean reattachment lengths X r /H obtained with the finer discretization are summarized in Tab. 1, where results from other studies are also included.MFS provides a value relatively close to the one taken from [12].Smaller reattachment lengths are obtained for DSM and RBVMM.present results numerical results experimental results MFS (F) RBVMM (F) DSM (F) DNS LKM95 [12] Exp KM95 [11] X r H 6.18 5.72 5.81 6.28 6.51 Differences between the results obtained with the methods considered here are observable for the more sensitive Reynolds-stress values, for instance u 1 u 2 , shown in Fig. 4. In the recirculation zone as well as behind the reattachment region, MFS provides the best approximation for both discretizations.RBVMM and DSM provide similar approximations.Moreover, convergence to the experimental data is observed for all methods.Evaluations of computing times reveal notably reduced computational cost compared to DSM and only marginally enhanced cost compared to RBVMM.

CONCLUSIONS
A variational multiscale method with multifractal subgrid-scale modeling has been proposed for large-eddy simulation of turbulent flow.In the multifractal subgrid-scale modeling approach, the subgrid-scale velocity is evaluated based on a multifractal reconstruction of the subgrid-scale vorticity and calculated via the law of Biot-Savart.The required smaller resolved scales have been identified by further separating the resolved scales via level-transfer  operators from plain algebraic multigrid methods.Moreover, a near-wall limit of the multifractal subgrid-scale modeling approach has been included, which allows for particularly taking into account the near-wall effects of turbulent flow.The proposed multifractal subgridscale modeling approach has been embedded into a residual-based variational multiscale formulation.This novel approach has shown the integration of structural subgrid-scale models into the variational multiscale framework.
Turbulent flow over a backward-facing step has been examined to evaluate the proposed method including the near-wall limit of the multifractal subgrid-scale modeling approach.The results have been compared to results obtained with the widely-used dynamic Smagorinsky model and the residual-based variational multiscale method.Experimental and numerical results have also been included for comparison.Only slight differences between the different methods have been observed for the mean velocity.Notable differences between the respective methods have occurred for the more sensitive values.For the component of the Reynolds-stress tensor, the proposed form of the multifractal subgrid-scale modeling approach has revealed notably better approximations than the two other methods.For the finer discretization, the reattachment length has also been approximated better by the proposed method than by the two other methods.
In our future work, we intend to extend the proposed method to scalar transport in turbulent flow, in particular for variable-density flow at low Mach number and combustion.

Figure 1 .
Figure 1.Geometry of backward-facing step and channel used to generate a turbulent inflow velocity profile.

Table 1 .
bottom wall of backward-facing step: red lines: coarser discretization; blue lines: finer discretization.Mean reattachment lengths Xr H of backward-facing step flow from present simulations and other numerical and experimental results.