## Abstract

A dissipative scale-similarity subgrid model was recently proposed in which only the dissipative part of the subgrid stresses was added to the momentum equations. This was achieved by adding the gradient of a subgrid stress only when its sign agreed with that of the corresponding viscous term. In the present work, this idea is used the other way around as forcing in hybrid large eddy simulation–Reynolds-averaged Navier–Stokes: only the part of a subgrid stress term that corresponds to back scatter is added to the momentum equations. The forcing triggers resolved turbulence in the transition region between the unsteady Reynolds-averaged Navier–Stokes and large eddy simulation regions. The new approach is evaluated for fully developed channel flow at *Re*_{τ}=4000. It is found that the forcing indeed does increase the resolved turbulence in the transition region. The magnitude of the production (i.e. back scatter) due to forcing in the equation for resolved kinetic energy is of the order of that due to the usual strain-rate production term. The present approach of using back scatter from a scale-similarity model can also probably be useful for triggering transition.

## 1. Introduction

In hybrid large eddy simulation–Reynolds-averaged Navier–Stokes (LES–RANS), forcing is often used in the transition region between the unsteady RANS (URANS) region and the LES region. The resolved turbulence in the URANS region is fairly small because of the large dissipation by the high turbulent viscosity; forcing is used to stimulate the growth of resolved turbulence when a fluid particle leaves the URANS region and enters the LES region. The use of forcing in hybrid LES–RANS can be compared with the use of turbulent fluctuating inlet boundary conditions in LES. The object in both cases is the same: to create resolved turbulence.

Scale-similarity subgrid models are based on the idea that the largest unresolved scales are *similar* to the smallest resolved scales, hence the name *scale similar*. When the first scale-similarity model was proposed, it was found that it is not sufficiently dissipative (Bardina *et al.* 1980). An eddy-viscosity model must be added to make the model sufficiently dissipative; these models are called *mixed* models. A dissipative scale-similarity model was recently presented (Davidson 2008) in which only the dissipative part of the subgrid stress term, ∂*τ*_{ik}/∂*x*_{k}, is added to the momentum equations.

This idea is also employed in the present work, but here the approach is used the other way around: only the part of the subgrid stress term that corresponds to back scatter is added to the momentum equations. Hence, the subgrid stress term acts as a *counter-gradient* diffusion term in the momentum equations.

## 2. Scale-similarity model

### (a) Dissipative model

The filtered Navier–Stokes read (⨪ denotes filtering)(2.1)where D/D*t* and *τ*_{ik} denote the material derivative and the SGS stress tensor, respectively. In the scale-similarity model, the SGS tensor is given by (Speziale 1985)(2.2)This model is not sufficiently dissipative. Let us take a closer look at the equation for the resolved, turbulent kinetic energy, , which reads (〈.〉 denotes averaging in time)(2.3)The first term on the third line is the non-isotropic (i.e. the true) viscous dissipation, *ϵ*^{non}; this is predominately negative. The first term on the right-hand side of the last line is the viscous diffusion term and the second term, *ϵ*, is the viscous dissipation term, which is always positive. The last term, *ϵ*_{SGS}, is a source term arising from the SGS stress tensor, which can be positive or negative. When it is positive, forward scattering takes place (i.e. it acts as a dissipation term); when it is negative, back scattering occurs.

Figure 1 presents SGS dissipation, *ϵ*_{SGS} in equation (2.3), computed from filtered direct numerical simulation (DNS) data. The forward scatter, , and back scatter, , SGS dissipation are defined as the sum of all instants when *ϵ*_{SGS} is positive and negative, respectively. is computed with (see equation (2.7)). As can be seen, the scale-similarity model is slightly dissipative (i.e. *ϵ*_{SGS}>0), but the magnitudes of the forward and back scatter dissipation are both much larger than *ϵ*_{SGS}.

One way to make the SGS stress tensor strictly dissipative is to set the back scatter to zero, i.e. max(*ϵ*_{SGS}, 0). This could be achieved by setting ∂*τ*_{ik}/∂*x*_{k}=0 when its sign is different from that of (see the last term in equation (2.3)). This would work if we were solving for *K*. Usually, we do not, and the equations that we do solve (the filtered Navier–Stokes equations) are not directly affected by the dissipation term, *ϵ*_{SGS}.

Instead, we have to modify the SGS stress tensor as it appears in the filtered Navier–Stokes equations (equation (2.1)). However, since it is the resolved *turbulent fluctuations*, i.e. *K* in equation (2.3), which we want to dissipate, we must consider the filtered Navier–Stokes equations for the fluctuating velocity, . It is the diffusion term in this equation that appears in the first term on the second line in equation (2.3). The viscous diffusion term is dissipative. To ensure that *ϵ*_{SGS}>0, we set −∂*τ*_{ik}/∂*x*_{k} to zero when its sign is different from that of the viscous diffusion term (cf. the first two terms on the third line in equation (2.3)). This is achieved by defining a sign function(2.4)where *M*_{ik}=±1. The problem is that we do not know until the simulations have been carried out. Fortunately, the sign of the second derivative of the resolved velocity fluctuation, , is mostly the same as that of the resolved velocity, . Figure 2*a* presents a comparison of the two second derivatives using DNS data. As can be seen, the r.m.s. of the second derivative of *u*′ is larger—or much larger—than that of . Figure 2*b* shows the correlation of the signs of the two second derivatives. It can be seen that the correlation is larger than 95 per cent for *y*^{+}>40. Hence, equation (2.4) can be replaced by(2.5)Each component of the divergence of SGS stress tensor in equation (2.1) is then simply multiplied by(2.6)so that(2.7)In order to avoid the sign function changing sign between two iterations within a time step, the second derivatives of in equation (2.5) are evaluated using velocities at the old time step. It should be noted that, since the limiter, , operates on each cell rather than on each face, the SGS diffusive fluxes, , are not conservative. However, this is unavoidable since we need to control the *net* force per unit volume, , rather than the stresses at the face, . It can further be noted that, by using equation (2.5) rather than equation (2.4), the model is no longer strictly dissipative in the equation. It is now only 95 per cent dissipative (see figure 3). However, the model is indeed—assuming that the diffusion term, is negligible—strictly dissipative in the equation.

By using the limiter, , we omit the back scatter caused by the SGS stresses. Figure 1 includes , which is computed with .

### (b) Back scatter model

LES is an appropriate method for flows that include turbulent fluctuations with large spatial scales. For this type of flow, a fairly coarse mesh can be used to resolve the large fluctuations. RANS in this type of flow is often a poor alternative, since the large fluctuations are difficult to accurately model with a RANS turbulence model. The problems for LES start when walls are introduced. The reason is that length scales of the largest fluctuations near walls are not large; hence, a fine mesh must be used to resolve these fluctuations.

Hybrid LES–RANS was invented to get around this problem. In this method, RANS is used near walls and LES is used in regions away from the walls. However, difficulties appear in the transition region between RANS and LES. When the flow goes from a RANS region to a LES region, the resolved turbulence is too small. To stimulate growth of resolved turbulence, forcing is often used (Piomelli *et al.* 2003; Batten *et al.* 2004; Davidson & Billson 2006). The forcing in the transition region is related to turbulent fluctuating inlet boundary conditions in LES: the object in both cases is to create resolved turbulence.

A dissipative scale-similarity model was presented above. Via the sign function, only the dissipative part of the subgrid stress term, −∂*τ*_{ik}/∂*x*_{k}, was included in the momentum equations. The subgrid stress term can also be used as a forcing term. In this case, only the back scatter part is included. This is easily accomplished by replacing equations (2.5)–(2.7) by(2.8)where superscript ‘B’ indicates back scatter. In this way, the subgrid tensor term is included whenever the sign of −∂*τ*_{ik}/∂*x*_{k} is *opposite* to that of the viscous diffusion. This means that the subgrid stress term, , in the momentum equations acts as a *counter-gradient* diffusion term.

In order to increase the strength of the forcing, the stresses in equation (2.2) are computed at the test level, i.e. at a grid half as fine as the computational grid. The stresses are filtered using the formula ( denotes here or and *ℓ* is the cell index)(2.9)This filtering is applied consecutively in the three coordinate directions.

## 3. Hybrid LES–RANS

The momentum equations with an added SGS/RANS viscosity read(3.1)where the last term is the forcing term, which is zero outside the forcing region. A one-equation model is employed in both the URANS region and the LES region and readswhere . The location at which the switch is made from URANS to LES is called the interface and is located at *y*_{ml} from each wall. In the inner region (*y*≤*y*_{ml}), *k*_{T} corresponds to RANS turbulent kinetic energy, *k*; in the outer region (*y*>*y*_{ml}), it corresponds to subgrid-scale kinetic turbulent energy, *k*_{SGS}. The coefficients are different in the two regions (see table 1). No special treatment is applied in the equations at the matching plane except that the form of the turbulent viscosity and the turbulent length scale is different in the two regions. *k*_{T}=0 at the walls.

## 4. The numerical method

An incompressible, finite-volume code is used (Davidson & Peng 2003). For space discretization, central differencing is used for all terms except for the convection term in the *k*_{T} equation, for which the hybrid central/upwind scheme is employed. The Crank–Nicolson scheme is used for time discretization of all equations. The numerical procedure is based on an implicit, fractional step technique with a multigrid pressure Poisson solver (Emvin 1997) and a non-staggered grid arrangement.

## 5. Results

In Davidson (2008), the dissipative scale-similarity model was applied to decaying grid turbulence, where it was found to be as dissipative as the Smagorinsky model. Below, the dissipative scale-similarity model is first applied to channel flow at *Re*_{τ}=500. In the second part of this section, the back scatter from the scale-similarity model is used for forcing in hybrid LES–RANS.

### (a) Channel flow at *Re*_{τ}=500. Dissipative mode

The Reynolds number is 500 based on half the channel height and the friction velocity. The mesh has 64×80×64 cells. The extent of the computational domain is 3.2 and 1.6 in the streamwise and spanwise directions, respectively. A stretching of 12 per cent is used in the wall-normal direction.

Figure 3 presents the velocity profiles obtained with the dissipative scale-similarity model, the dynamic model and with no model. No converged results could be obtained with the standard scale-similarity model. As can be seen, no model gives perfect agreement with DNS and the log law. Hence, this flow is not a good test case for evaluating the accuracy of SGS models. Here, it is used to analyse the dissipative scale-similarity model. The dynamic model gives slightly better agreement with DNS than the dissipative scale-similarity model.

Figure 4*a* presents the momentum diffusion terms close to the wall. It can be seen that the SGS diffusion term evaluated using the standard scale-similarity model is of opposite sign to that of the viscous diffusion. When introducing the sign function in equations (2.5)–(2.7), it can be seen that the SGS diffusion term takes the same sign as the viscous diffusion term for *y*^{+}>10. The fact that the two terms have opposite signs for *y*^{+}<10 simply means that the viscous diffusion is very large at instants when the SGS diffusion term is set to zero. The diffusion due to the resolved shear stress is included in the figure. It is, as can be seen, much larger (more than five times) than the SGS term.

Figure 4*b* compares the SGS dissipation from the scale-similarity model with that from the dissipative scale-similarity model (recall that the simulations were carried out with the latter model). As can be seen, the SGS dissipation is indeed much larger with the dissipative model than with the standard model. For comparison, the SGS dissipation where the SGS viscosity is computed with the Smagorinsky–Lilly (*C*_{s}=0.1) model is also included. For *y*^{+}<20, the SGS dissipation from the Smagorinsky is larger than that from the dissipative scale-similarity model, but the situation is reversed further away from the wall.

### (b) Channel flow at *Re*_{τ}=4000. Back scatter mode

The Reynolds number is 4000 based on half the channel height and the friction velocity. Pure LES is not suitable at such a high Reynolds number. Instead, hybrid LES–RANS is used, in which the inner region is modelled by a one-equation RANS model (*k*) and the outer region is modelled by a one-equation LES model (*k*_{SGS}). The location of the interface, *y*_{ml}, has been chosen along a fixed grid line, either at (17 cells in the inner region at each wall) or at (22 cells in the inner region). The last value corresponds to the location of the matching line in DES, i.e. where 0.65*Δ*=*κy* (*Δ*=max{Δ*x*, Δ*y*, Δ*x*}). A 64×80×64 (*x*, *y*, *z*) mesh was used. The size of the domain is 6.4×2×3.2; 15 per cent stretching is used in the *y* direction. The forcing term, , is added in the region for which *y*_{ml}<*y*<2*y*_{ml} (six cells). Outside this region, the forcing term is set to zero.

Figure 5 presents the velocity profiles and the streamwise resolved normal stresses. As can be seen, the forcing increases the resolved streamwise fluctuations. The resolved shear stress is also increased, which increases the wall-normal diffusion, giving a fuller velocity profile (a reduced centreline to friction velocity, 〈*u*_{C}〉/*u*_{τ}). The forcing has a much greater effect on the resolved streamwise fluctuations when the matching line is located at than when it is at . The explanation is found in figure 6*a*, which shows the diffusion due to the forcing term and that due to the resolved shear stress. The forcing diffusion is much larger for than for , presumably because the velocity gradient is larger at the former location. In both cases ( and 270), it can be seen that—except at the two innermost cells—the diffusion due to forcing is (almost) as large as that due to the resolved shear stress.

Figure 6*b* presents some terms in the *K* equation in the forcing region. Production term , including the standard scale-similarity stresses, is small and—on average—dissipative (this term has been computed only for post-processing reasons; it is not used in the simulations). However, production term is, as expected, positive and hence it contributes to the generation of resolved turbulence. For comparison, the production term due to the resolved shear stress, 〈*u*′*v*′〉, is included.

The sign function, , in equation (2.8) governs when the forcing term is non-zero. How often does it happen that and how often does it switch between 0 and 1? Figure 7*a* shows component , which affects the wall-normal SGS diffusion in the streamwise momentum equation. As can be seen, it is larger than 0.5 for both cases except for the innermost node, where it attains a value of 0.2. The reason for this small value is that the viscous diffusion of becomes large at the interface between the URANS and the LES regions because of the large gradients occurring there. Hence, the non-isotropic dissipation, *ϵ*_{non}, and the isotropic one, *ϵ*, differ a great deal, cf. equation (2.3) (at this location, is actually negative). Since we are forcing the SGS term, , to have the same sign as *ϵ*_{non}, the SGS term and *ϵ* may take opposite signs if the viscous diffusion is large.

The fact that is zero almost half the time raises the question: how often does change sign? Does it flip-flop and change sign every second time step? Figure 7*b* shows that—fortunately—this is not the case. The sign function stays positive, on average, for 10 consecutive time steps, which corresponds to 50 viscous time units.

Above, the extent of the forcing region was set to six cells. Figure 8 presents the results when the forcing region is reduced to four cells and extended to eight cells, respectively. As can be seen, the velocity profile (figure 8*a*) for is in excellent agreement with the log-law profile when the forcing region is extended, and it becomes somewhat closer to the velocity profile without any forcing when the extent of the forcing region is reduced (cf. with figure 5*a*). The resolved streamwise fluctuations, *u*_{r.m.s.}, increase both when the extent of the forcing region is increased (eight cells) and when it is decreased (four cells); see figures 5*b* and 8*b*. In the former case, *u*_{r.m.s.} increases simply because the integrated effect of the forcing is increased. In the latter case, the increase in the peak of *u*_{r.m.s.} can probably be explained by a sharper velocity gradient outside the interface (cf. solid lines in figures 8*a* and 5*a*). For , a reduction in the extent of the forcing region to four cells improves the predicted velocity profile (see figures 8*a* and 5*a*). However, when the extent of the forcing region is extended to eight cells the method fails completely, as can be seen in figure 8*a*. This failure can probably be blamed on the positive feedback enhanced by the streamwise periodic boundary conditions. In ongoing work on channel flow with inlet and outlet, forcing can be used in the entire channel without any such problems.

## 6. Concluding comments

The paper presents a new approach to extract either the forward scatter (dissipation) or the back scatter given by the scale-similarity model. When the SGS term in the momentum equation has the same sign as the viscous diffusion term, it is forward scatter and, when the two terms have opposite signs, it is back scatter. In forward scatter mode, the model is used as a usual dissipative SGS model. In back scatter mode, the model is used as forcing in hybrid LES–RANS, stimulating the generation of resolved turbulence.

The idea of locally either dampening the resolved fluctuations (forward scatter) or stimulating the growth of resolved turbulence (back scatter) can be useful in URANS and DES. In these methods, the momentum equations are triggered through instabilities to go unsteady in regions where the grid is fine enough. It can then be useful to add back scatter from the scale-similarity model to promote this transition. On the other hand, if the grid is coarse, it may be a good idea to ensure that the momentum equation operates in fully steady mode. In this case, some extra dissipation from the dissipative scale-similarity model may be needed to damp any unwanted resolved fluctuations.

## Acknowledgments

Financial support by Swedish National Infrastructure for Computing (SNIC) for computer time at C3SE (Chalmers Center for Computational Science and Engineering) is gratefully acknowledged.

## Footnotes

One contribution of 16 to a Discussion Meeting Issue ‘Applied large eddy simulation’.

- © 2009 The Royal Society