{'fraction_non_alphanumeric': 0.02709378920953576, 'fraction_numerical': 0.006579360100376411, 'mean_word_length': 4.7084471104346655, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Active model split hybrid RANS/LES

Reliably predictive simulation of complex flows requires a level of model sophistication and robustness exceeding the capabilities of current Reynolds-averaged Navier-Stokes (RANS) models. The necessary capability can often be provided by well-resolved large eddy simulation (LES), but, for many flows of interest, such simulations are too computationally intensive to be performed routinely. In principle, hybrid RANS/LES (HRL) models capable of transitioning through arbitrary levels of modeled and resolved turbulence would ameliorate both RANS deficiencies and LES expense. However, these HRL approaches have led to a host of unique complications, in addition to those already present in RANS and LES. This work proposes a modeling approach aimed at overcoming such challenges. The approach presented here relies on splitting the turbulence model into three distinct components: two responsible for the standard subgrid model roles of either providing the unresolved stress or dissipation and a third which reduces the model length scale by creating resolved turbulence. This formulation renders blending functions unnecessary in HRL. Further, the split-model approach both reduces the physics-approximation burden on simple eddy-viscosity-based models and provides convenient flexibility in model selection. In regions where the resolution is adequate to support additional turbulence, fluctuations are generated at the smallest locally resolved scales of motion. This active forcing drives the system towards a balance between RANS and grid-resolved LES for any combination of resolution and flow while the split-model formulation prevents local disruption to the total stress. The model is demonstrated on fully-developed, incompressible channel flow [1] and the periodic hill [2], in which it is shown to produce accurate results and avoid common HRL shortcomings, such as model stress depletion.


I. INTRODUCTION
It has long been recognized that current Reynolds-averaged Navier-Stokes models (RANS) are inadequate for the prediction of complex turbulent flows. On the other hand, sufficiently well-resolved large-eddy simulation (LES) models provide an accurate representation of turbulence in many circumstances, but are commonly too expensive to apply in practice. Turbulence modeling methods that allow for a flexible balance between resolving some of the turbulent fluctuations, as in LES, and modeling their effect, as in RANS, offer an attractive compromise between these approaches. In particular, by resolving turbulence fluctuations only where RANS models are deficient, significant improvements in mean flow predictions are possible with minimal sacrifice of computational efficiency. For instance, RANS models are well-known to be inadequate in regions of flow separation, reattachment, or three-dimensionality in the mean [3][4][5][6][7]. Enabling LES to be active in such regions will effectively avoid the deficiencies of the RANS models, while also avoiding the cost of using LES resolution where RANS models are sufficient.
Due to the combination of limited computational resources and the need for high fidelity simulations, methods that partially resolve turbulent fluctuations can be expected to remain a necessity for engineering applications until the turn of the next century [8]. In response, a myriad of such turbulence modeling techniques have been proposed and developed with varying degrees of success. The techniques fall into two main categories: wall-modeled large eddy simulation (WMLES) [9][10][11][12][13], which focuses specifically on relieving original work, a transport-based subgrid model was used for the fluctuating portion while a wall-sensitive velocity-difference model was used for the mean. The main advantage of this structure was to allow for a more standard Reynolds stress closures to take over near walls were grid scales become large in comparison to some characteristic mixing length scale. Thus, this may actually be considered the first hybrid method. Until recently, mean and fluctuating stress decomposition methods have only seen pure-LES application [45,46].
The TVHRL [40] method blends distinct RANS and LES models (category 2 above), and builds on the innovation of [44] to have each model act on a different portion of the resolved strain rate tensor. In particular, the RANS eddy viscosity acts only on the mean strain rate while the LES eddy viscosity acts only on the resolved fluctuating strain rate. Additionally, all quantities entering the RANS model are mean values while a standard Smagorinsky model acts on the fluctuating portion of the strain. This spitting of the model was designed to allow the resolved turbulent stress to develop independently from the RANS viscosity. The authors further noted that the portion of the model acting on the mean is intended to provide the entirety of the subgrid contribution to the mean stress while the fluctuating portion is to contribute only to the rate of transfer of energy from the resolved turbulent motions to the unresolved scales. This assumes, however, that the fluctuating eddy viscosity is uncorrelated with the resolved strain rate magnitude, which is not the case for the Smagorinsky model as used in that work. Time averaging over many eddy turnovers was used to establish the mean velocity. The model still made use of an ad-hoc hyperbolic tangent lengthscale comparison blending function to switch between RANS and LES models, which can result in drastic model shifts depending on the local ratio of model-to-grid length scales. Nonetheless, as discussed further in Sec. III A and Sec. IV, the two velocity model splitting approach is significant and may indeed be a necessity for any LES based on eddy viscosity models (EVM), when a non-negligible portion of the mean turbulent stress must be provided by the model.
The RSC-SGS model [41] eliminated the ad-hoc blending used in [40] by directly reduced the mean-strain contribution, as calculated with either the SA model [47] or the van Driest mixing-length approximation, to the total stress by subtracting an averaged resolved stress. The dynamic Smagorinsky model (DSM) [48] was used on the fluctuating strain. However, the fluctuating portion of the model will again contribute to the mean Reynolds stress due the correlation between the Smagorinsky model and resolved strain. Without this contribution, the approach will yield the unaltered RANS stress in expectation. The method also relied on specification of the wall-normal location where the model transitioned from the "constrained" mixed model to a standard DSM. Thus, from the hybrid perspective, this would be considered a zonal approach.
The DHRL method [42,43] extends the two-velocity approach by blending between RANS, acting only on the mean, and an implicit SGS model (monotonically integrated LES [49]) using the ratio of resolved turbulent production to the difference between the production resulting from use of only the RANS model and the production from use of only the SGS model. Thus, the state of the resolved field is considered in the model blending. Further, both models are active throughout the entire domain. The DMHRL [18] approach represents the extreme of interfaced methods by coupling distinct LES and RANS simulations on separate grids. Additional forcing terms are added to each set of governing equations which enforce consistency of the two simulations in expectation. In the LES simulation, regions are designated as either RANS or LES regions based on a measure of the local resolution. In designated LES regions, the difference between the RANS simulation velocity and time-average LES simulation velocity is used along with a time scale to construct an artificial forcing acceleration that is added to the RANS simulation. In regions designated as under-resolved in the LES simulation, the difference between the total mean stress, as determined by the RANS model and the subgrid scale model plus the resolved fluctuations, is added to the LES simulation to effectively enforce the RANS stress.
Despite extensive modeling developments, hybrid models based on blending between RANS and LES have not led to generally reliable HRL methods capable of traversing through arbitrary levels of resolved turbulence. There are two key reasons for this shortcoming. First, HRL models often rely on passive generation of resolved turbulence in the presence of reduced modeled stress. This approach necessarily leads to inconsistency between the modeled, resolved, and total turbulent stresses, as exhibited in regions of modeled-stress-depletion. This issue is described further in Sec. II A. Second, most HRL approaches rely on scalar eddy viscosity models to simultaneously represent both the unresolved portion of the mean turbulent stress as well as the transfer of energy from the resolved fluctuations to the subgrid, with the notable exception of the TVHRL method and related approaches described previously. However, a single eddy viscosity is insufficient for this task, as discussed in Sec. II B. Additional problems arise due to ad-hoc blending approaches, discussed in Sec. II C, and inappropriate application of transport equations formulated to govern mean turbulence properties in RANS closure models to fluctuating quantities, as described in Sec. II D.
Here, we introduce an HRL method that directly addresses these common issues, which we will refer to as the "active model split" (AMS) method. The remainder of the paper is organized as follows. Section II describes in more detail the common HRL deficiencies that AMS is designed to address. Section III discusses general modeling issues that arise for partial turbulence resolving models, and details of AMS are described in Section IV. Results of applying AMS in channel flow and a periodic hill case [2] are evaluated in Section V and in Section VI conclusions and opportunities for further developments are discussed.

II. HYBRID MODELING ISSUES
In this section, we recall several common issues with HRL and propose strategies to either correct or circumvent them. As this work is primarily focused on a framework for predictive HRL, those issues inherent to either RANS or LES are not considered here to allow focus on those unique to hybrid methods. However, as discussed in Sec. III, formulating a robust HRL has led to modeling concepts applicable to LES in general. The particular hybrid modeling issues addressed here are those of (1) maintaining consistency between resolved and modeled turbulence, (2) over-reliance on simple eddy viscosity models, (3) ad-hoc RANS to LES blending, and (4) misuse of the RANS transport models. Many of these issues were first identified in [50] with the exception of the first, which has been reported on extensively elsewhere [32,35,51,52].

A. Consistency of Resolved and Modeled Turbulence
For any scale resolving turbulence model to be predictive, the resolved and modeled contributions to the mean turbulent stress (the Reynolds stress) must be consistent. That is, they should sum to the correct total stress. As more turbulence is resolved, the resolved contribution to the Reynolds stress should increase and the modeled Reynolds stress should decrease by the same amount, until ultimately all of the Reynolds stress is being carried by resolved turbulence fluctuations as in a DNS.
In typical HRL methods, the ability of the simulation to resolve turbulence fluctuations in some region is signaled by a reduction in the modeled Reynolds stress. In response, assuming the mean shear is large enough, natural flow instabilities will then lead to resolved fluctuations. By construction, such HRL methods must exhibit regions of modeled-stress-depletion (MSD) [29,34,51] as the instabilities cannot develop without the total stress being depleted first. In fact, this is the best case scenario, as it is possible that turbulent regions of flow exist where reducing the model stress to zero will not result in the development of resolved fluctuations. Simulations by a variety of researchers with disparate models have shown that some type of active forcing to generate resolved fluctuations is needed in simulations of mixed levels of resolved turbulence [15,52,53]. It appears that explicitly introducing resolved fluctuations, at a rate consistent with the reduction of modeled stress is the only way to prevent regions of MSD. Otherwise, common consequences of MSD such as log-layer mismatch, reduced body forces, premature flow separation, and delayed flow reattachment will persist in HRL.
Existing forcing methods typically address this issue by generating fluctuations only at a prescribed LES inlet [38,[54][55][56]. In addition to requiring specification of the distinct LES region and the entire spectrum of locally resolved synthetic turbulent scales and intensities, the LES inlet location must be sufficiently upstream of the flow features of interest so that the artificial inflow condition can "heal" to a realistic state. An alternative is to use spatially distributed body forcing to introduce synthetic turbulence slowly so that the fluctuations maintain a realistic turbulent state through the entire forcing region. In this way, none of the simulated domain and associated computational cost is sacrificed to the healing process, which should result in an accurate solution in the entire simulated domain. Further, the simulation domain size can be reduced or turbulence-resolved regions limited, since no healing regions are needed, While it is clear that some forcing is required, the method and form it should take is not. The problem of introducing turbulent structures in a hybrid simulation requires these structures to be artificially constructed. While one can make principled estimates of the appropriate strength of the body forcing based on the smallest resolved scales and the length and time scales of the modeled turbulence, the structure of the forcing must be prescribed.
Formally, body forcing can be considered a model for the commutator of the filter operator, denoted ·, defining the resolved scales of turbulence and the substantial derivative in the momentum equation, C e = Du i /Dt − Du i /Dt, where here the filter operator is the projection of the turbulent field onto the discrete representation of the numerical solution (a so-called implicit filter, see [57]). From this perspective, the only way to avoid "making up" fluctuations would be to have the full DNS solution as a function of time, which would, of course, obviate the LES. However, formulating the forcing to be as representative of the commutator as possible should result in more realistic forced fluctuations which can be introduced more rapidly while maintaining a realistic resolved turbulent state. Improving the prescribed forcing structure will allow the use of shorter regions of LES resolution upstream of flow features of interest, such as flow separation.
There is also evidence that even unrealistic forcing structures may be sufficient to reduce inconsistencies between modeled and resolved turbulence. For instance, in channel flow, Piomelli et.al. [52] found that introducing a region of broad spectrum, tuned stochastic forcing in a thin region along the RANS/LES interface region in an HRL eliminated log-layer mismatch in fully developed channel flow. While this forcing is problem dependent and difficult to implement in general, this numerical test showed that forcing with random fluctuations in regions of transition from RANS to LES can address problems related to MSD. The notion of adding user-specified regional body forcing was further explored by Menter et.al. in the SAS-F model [15]. In their work, a broad spectrum of velocity fluctuations were produced using the random flow generator method (RFG) [58] and used to construct an acceleration term for each discrete timestep. The unique aspect of their construction is that the basic SAS model form allows for the model to detect the added fluctuations and naturally respond by lowering the internal model length scale and corresponding eddy viscosity. Even so, channel flow results still exhibited log layer mismatch.
In sum, inconsistency in the resolved and modeled turbulence arise from combinations of reliance on passive self-generation of resolved turbulence in the presence of reduced modeled stress and the overarching hybridization strategy that only manipulates the modeled turbulence. In addition to representing the unresolved contribution to the turbulent stress, a hybrid method should ensure accurate representation of the total turbulence stress. We propose that body forcing is a necessary component of an HRL framework and that modeled stress must only be reduced in response to the presence of locally resolved fluctuations. As presented in Sec. IV D, a body forcing method is formulated with no explicit time step dependance and a structure that only introduces fluctuations at the scale of the smallest locally resolved turbulence, thereby reducing the potential corruption of existing, and presumably realistic, resolved turbulence. Further, no user specification of forcing-regions is necessary.

B. Simple Eddy Viscosities
The second issue stems from a limitation of simple eddy viscosity models that has been known in the context of LES for at least two decades [59], but has been ignored in hybrid model formulations. HRL requires the model to both make a non-negligible contribution to the unresolved mean turbulence stress (Reynolds stress) and represent the transfer of energy from the resolved scales of motion to the modeled subgrid scales. In LES, simple eddy viscosities only represent the latter while significantly under-predicting the former. Jimenéz and Moser [59] demonstrated that the dynamic Smagorinsky model produces subgrid stress values that are nearly five times smaller than the true filtered stress for fully-developed channel flow. In these cases, the LES is "well-resolved" in that the contribution of the subgrid turbulence to the mean turbulent shear stress is small, so that the large relative errors in the mean subgrid shear stress have minimal impact on the mean velocity. The ability of subgrid models to accurately predict the dissipation but not the momentum transport has also been observed in experimental measurements of a turbulent plume [60]. The canonical SGS model test case of homogeneous isotropic turbulence (HIT) further demonstrates this point as the divergence of the mean stress is zero by homogeneity and any contribution to local stress from the model is irrelevant. In either of the above cases, the function of the subgrid "stress" model is only to model the transfer of energy from the resolved to unresolved scales, with no concern for the mean stress. On the other hand, even for flows where an eddy viscosity model provides a good approximation of the total mean stress for RANS, there is no reason to expect the same model to provide the correct energy transfer for any resolved fluctuations.
The difficulty arises because the statistical characteristics of the mean subgrid stress and the mean energy transfer are different. The latter is governed by the correlation coefficient between the subgrid stress and the resolved strain rate tensor, as opposed to the former which is just the mean of the stress. In filtered turbulence, the correlation coefficient is small (of order 20% or less [59,61]), but with an eddy viscosity model, the correlation coefficient is one or nearly so. To get the energy transfer rate right, it is necessary to reduce the eddy viscosity compared to that needed to get the magnitude, and the mean, of the subgrid stresses right. Therefore, a simple eddy viscosity model is not equipped to model both stress and energy transfer to the subgrid scales, as commonly demanded in HRL. It seems natural then to divorce these two model functions in HRL. That is, rather than use a single SGS model, introduce a model primarily responsible for the mean stress, operating only on the mean gradients, and a model primarily responsible for energy transfer, operating only on resolved fluctuating gradients. This "model-split" approach is presented in Sec. IV.

C. Model Blending
Third, the methods used to reduce the model stress in HRL often involve ad-hoc blending functions and simplified resolution indicators. Typically, a blending function is specified to internally modify a model's behavior, transition between two separate models, or directly reduce the modeled Reynolds stress. In all cases, it has a direct effect on the physics of the flow being simulated, and unintended consequences may arise.
In some methods [25,36], the complexity of determining the correct blending provided the grid, domain, model, and flow is placed entirely on the user with ab intitio specification of the level of resolved turbulence. The optimal balance of resolved and modeled turbulence will generally be a function of space and time making explicit specification by the user difficult.
To avoid this explicit specification, it is more common to construct the blending function based on a resolution indicator that is intended to determine when and where LES fluctuations can be resolved based on the local grid and flow. However, commonly used indicators use simplified measures of the generally anisotropic resolution and resolved turbulence scales, if they make use of information from the resolved fluctuations at all. For instance, the most commonly used grid length scale, the cube-root of the cell volume [30,62], tends to the smallest cell length scale as the cell aspect ratio increases. Therefore, it implicitly assumes structures of the smallest grid scale are resolved in every direction, which over-estimates the effective resolution of the grid. From this perspective, grid measures based on the cell diagonal are preferable. Certainly, with highly anisotropic resolutions common to near-wall meshes, different definitions will provide different results for any particular combination of grid and flow. In determining the turbulence length scale, the state of resolved field is often neglected in favor of turbulence length scales derived from the model. This practice implicitly assumes isotropy of the resolved scales, and consistency between the resolved and unresolved turbulence. However, the resolved turbulence is generally not isotropic both because the large scales are generally not isotropic and because anisotropic resolution will necessarily yield anisotropy. Anisotropy of both the resolution and the resolved turbulence should clearly be accounted for in a resolution indicator.
Ideally, HRL methods should be constructed that eliminate blending entirely. We do so here with the aforementioned model-split approach, which yields two separate models, both active everywhere in space and time so that no transition functions need be specified. Additionally, ambiguous grid measures are eliminated by introducing a resolution tensor that represents the anisotropy of the resolution. Finally, the indicator metric depends on both modeled and resolved quantities everywhere in the flow.

D. Transport Models
The final issue is the use of transport equations formulated to govern mean turbulence properties for RANS models. These models can dramatically misbehave when used with resolved fluctuations. This problem is specific to single-model unified HRL approaches based on RANS models, a class including many common methods (DES [34], PANS [24], PITM [25]). However, the fundamental issue has largely been ignored or unrecognized because its symptoms have been ameliorated with numerical dissipation, coefficient tuning for a particular numerical method and code, or the averaging of selected terms.
In these HRL methods it is often tacitly assumed that nonlinear source terms in the turbulence model transport equations, which have been constructed and calibrated to operate on mean quantities in a RANS model, will perform similarly when applied to fluctuating quantities. Of course, these terms are themselves models for actual terms in transport equations, and in general, their use with fluctuating quantities will result in spurious correlations with non-negligible effects on the mean behavior. For instance, consider the model for the destruction of dissipation in the ε transport equation [63] which is used in all k-ε based RANS models. The model has little physical justification but after coefficient tuning, it is functional for RANS in conjunction with models of other unclosed terms. The impact of applying this model to fluctuations can be seen by expanding 1/k in a Taylor series about 1/ k and retaining up to quadratic terms in k . The result is Thus, the mean destruction of dissipation becomes a direct, and complicated, function of the amount of resolved turbulence. Even the standard RANS production term, P k = 2ν t S 2 , is a model and when used with fluctuating quantities will exhibit spurious correlations, depending on the particular RANS model used to determine ν t . This difficulty would appear to be inconsistent with the fact that transport-based k-SGS models have been used successfully for some time [62]. However, in those applications, the model length scale is taken as the grid scale, which results in production and dissipation of kinetic energy that scale like k 1/2 and k 3/2 , respectively. The stronger dependence of dissipation on k effectively damps the k fluctuations, minimizing the effects of this problem by effectively adjusting the mean of each modeled term so that the length scale will remain as prescribed. In a hybrid context, this approach is not possible due to the lack of an intrinsic model length scale, unless the local simulation is fully LES or RANS. Among other issues, such an approach implies production scaling with k 1/2 which is entirely incorrect when the SGS length scale is between the integral and grid length scales.
To fully address the problem of using RANS transport models in HRL without resorting to problem-and numerics-specific coefficient tuning for all non-linear model terms, such models should be used to represent mean quantities as a function of only mean quantities, as they were designed to do. This is not to say that they have no place in HRL, nor that they cannot function as subgrid models, only that what goes in to, and what comes out of, RANS-based transport models must be expected values.

E. An Alternate Hybrid Modeling Formulation
Motivated by the HRL difficulties discussed in Sec. II A-II D, an alternate approach has been developed. The proposed HRL formulation includes: 1) a generalizations of the resolution adequacy indicator used in DES [17], 2) different models for use in RANS and LES similar to TVHRL and DHRL though without blending, and 3) forcing to generate resolved fluctuations as in SAS-F [15] though with a different structure. Referred to as the active model-split (AMS) formulation, the critical features are: 1) the use of distinct turbulence models to act on different portions of the resolved velocity gradient tensor, and 2) directly incorporating forcing into the model formulation so that the modeled stress can be reduced only in response to the local level of resolved turbulence. Resolved turbulence is continuously generated where indicated by a resolution adequacy parameter which accounts for the anisotropy in the resolution and resolved velocity field along with the model quantities. Before presenting the new formulation in Sec. IV, we present a simple argument which motivates a significant change to LES in general and directly leads to a merger of RANS and LES modes of operation in HRL.

III. GENERAL CONSIDERATIONS FOR TURBULENCE-RESOLVING METHODS
To begin, we argue that the model-split approach offers an attractive paradigm for constructing general turbulence resolving methods. Toward this end, consider the filtered Navier-Stokes equations: where u i denotes the filtered velocity, p is the filtered pressure, τ visc ij is the filtered viscous stress, and τ sgs ij = −(u i u j − u i u j ) represents the subgrid stress. It is common to represent the deviatoric part of the subgrid stress with an eddy viscosity based model: where 2k sgs is the trace of −τ sgs ij , S ij is the filtered strain rate tensor, ν sgs is the eddy viscosity, and the filter operation has been assumed to commute with differentiation. Of course, this model form cannot be correct in general, simply because it relates two tensors through a scalar. Further, it has been shown by a priori tests [59,61] that common eddy-viscosity-based subgrid models (e.g., Smagorinsky) do not lead to high correlation between the modeled and actual subgrid stesses.
Fortunately, accurate prediction of the mean flow does not require that (3) hold in an instantaneous sense. Instead, it is only necessary to match important statistics involving the subgrid scales. At a minimum, for general flows, it is necessary that two conditions are satisfied. First, the mean of the subgrid contribution to the total mean stress anisotropy must be well-represented, to ensure that errors are not introduced into the mean momentum balance. This implies that Second, the contribution of the subgrid stress to the mean rate of transfer of energy from the resolved to the unresolved scales must also be well-represented, such that the resolved fluctuations have the correct energy, or Other statistics may also be important. For example, one may also wish to correctly capture the anisotropy of the transfer from the resolved to unresolved turbulence. However, it is sufficient for the present purposes to consider just (4) and (5), because they are enough to reveal both the deficiencies of common approaches and one of the major benefits of model-splitting. If a model existed such that (3) were instantaneously true or nearly so, then (4) and (5) would be satisfied automatically. However, as discussed in Sec. II B, it is known that existing RANS and LES models do not do a good job at simultaneously predicting both the mean subgrid stress and the energy transfer rate. Indeed, they are generally not even designed with this goal in mind. RANS models are built entirely to capture (4). In RANS, since the only "resolved" component is the mean, an accurately predicted mean stress is sufficient to accurately predict the transfer from resolved to unresolved, which is just the usual production of TKE. However, once the state is fluctuating, there is no longer any reason to expect a RANS-like eddy viscosity to give a good prediction of the energy transfer, even if it is adequate for the mean subgrid stress.
Alternatively, LES models are commonly built for and tested on either homogenous flows, where the mean subgrid stress does not vary in space and is therefore dynamically insignificant, or for well-resolved simulations, where the mean subgrid stress is negligible. Thus, common LES models are calibrated to capture the energy transfer rate, which leads to poor predictions of the mean subgrid stress, as discussed in Sec. II B. Advanced subgrid models [48,[64][65][66], developed for transitional and wall-bounded flows, are constructed so the eddy viscosity vanishes where both the resolved and modeled stress become dynamically insignificant. Again, the model construction does not speak to the general scenario of both the modeled and total stress being dynamically significant. Successful application of such models to well-resolved LES has no bearing on HRL.

A. Model Splitting
In general turbulence resolving simulations, one cannot expect to have good LES resolution everywhere, so at least in some regions of the flow, both (4) and (5) are important. This is especially true in a HRL as the absence of such a region would mean the simulation is just a regular well-resolved LES. Further, in the context of HRL methods specifically, there is no reason to expect an ad hoc blending of models designed to work in the limits of RANS or well-resolved LES will succeed in simultaneously capturing (4) and (5) in the intermediate regime. Instead, recognizing these requirements motivates splitting the model into two parts: one aimed at (4) and the second at (5). In general, these two models need not be based on eddy viscosities. For instance, upwinding of the convective term is known to introduce numerical diffusivity [67] and could conceivably be used as the energy transfer portion. But, since one of the aims of this work is to develop an HRL which works with existing RANS models, we proceed in the vein of eddy viscosity models. In this situation, the goal is to develop two eddy-viscosity-based models, one denoted ν s to represent the mean subgrid stress, and a second denoted ν e to capture the energy transfer from the resolved fluctuations to unresolved scales where the greater than superscript, (·) > , denotes the resolved fluctuation of the quantity about its mean. Note that the component of production of subgrid energy due to the mean, which is given by 2ν s S ij ∂ j u i , is not part of (7), but if (6) is satisfied, the corresponding modeled production will also be automatically be correct.
The above relations (6) and (7) suggest the following instantaneous model: In general, due to correlations between the eddy viscosities and the velocity fluctuations, this form would result in a the contribution from ν s to the energy transfer, 2 ν s S ij ∂ j u > i , and a contribution from ν e to the stress, 2 ν e S > ij . To eliminate these contributions, we require that both ν s and ν e not fluctuate, i.e. they must only depend on mean quantities or be explicitly averaged after construction. While not disqualifying, such correlations would complicate the roles, and hence the construction, of the two desired models.
As with (3), there is no reason to believe that (8) represents a good approximation of the instantaneous subgrid stress, but it has two substantial advantages. First, for flows with a one-dimensional mean, it gives sufficient flexibility to simultaneously satisfy the two requirements of predicting the mean subgrid stress and the mean energy transfer rate. For complex three-dimensional flows, at a minimum, the mean production ( u i u j S ij ) and subgrid energy transfer constraints can be simultaneously satisfied. Second, since ν s is responsible for the mean subgrid stress and ν e is responsible for the portion of the energy transfer not associated with the mean, these components map naturally to existing RANS and LES models, offering a formulation in which the models can function together simultaneously rather than requiring ad hoc blending.

B. Mean Stress Scaling
The model-split form provides a natural segregation of the model in a partial turbulence-resolving simulation into "RANS-like" and "LES-like" components that can function simultaneously. However, it does not give any insight into how these components should scale with the level of resolved turbulence, which is important to the performance of any practical model. To deduce this scaling, we apply basic eddy viscosity arguments to a decomposition of the subgrid stress tensor.
Begin by decomposing the velocity field into its mean and fluctuating parts: u i = u i + u i . Further, apply the standard triple decomposition: where u > i denotes the resolved fluctuation, and u < i denotes the subgrid fluctuation, so u i = u > i + u < i . In this notation, the filtered field is given by In the following, it is assumed that u > i = 0, with the result that u i = u i . While this does not hold for general filters applied to general fields, one can always construct a filter that is consistent with this assumption by applying any standard filter to just the fluctuating field and then adding the mean. Applying this expansion to the true subgrid stress term, we have where B is a generic "breaking operator" given by B(a, c) = ac − ac. The first three terms are those most commonly considered as constituting the subgrid stress tensor. The four terms on the second line written in terms of B arise because of non-idealities that may occur when the mean is not well resolved in the LES and due to the non-linear interaction of the resolved fluctuations. The standard Leonard stress is formed from [68]. Stress contributions left in B-form are akin to aliasing errors in the numerical discretization of nonlinear terms, and are often neglected or explicitly discarded with dealiased numerical approximations. The term B(u > i , u > j ) may be of modeling interest to implicitly filtered LES. However, we currently neglect all terms on the second line of (9) to obtain To ensure that a turbulence-resolving simulation produces the correct mean flow, it is necessary that the expected value of the filtered Navier-Stokes equations including the modeled subgrid stress, is consistent with the RANS equations. In terms of the total convective flux, this requires which, using the triple decomposition and properties of the filter introduced above, implies that A turbulence-resolving model is an improvement over RANS due to the direct representation of the resolved Reynolds stress, (u > i u > j ), but we are left with the difficulty of modeling the expected value of the subgrid stress τ sgs ij , that is In the model-split paradigm, the mean subgrid stress is modeled separately, which, assuming an eddyviscosity form for the model is The objective is then to model ν s , and we do so by considering the theoretical basis for a linear eddy viscosity model.
Gradient-diffusion models represent the expected convective transport of some conserved fluctuating quantity, φ, by velocity fluctuations through a linear approximation of their covariance [69][70][71], where χ j is the stochastic displacement of a fluid particle when φ fluctuations become decorrelated with themselves. This decorrelation occurs over a time T φ , which is defined such that φ (0)φ (T φ ) / φ 2 1. Alternatively, T φ may be taken as the the integral timescale of the Lagrangian two-point correlation between the i th and j th velocity components and ensemble averaged over the Lagrangian particles [71]. However, such a timescale does not ensure the validity of modeling the fluctuating convective transport of some quantity through its gradient. It seems more natural for the timescale to be determined from the conserved quantity leading to the substitution u i χ j = u i u j T φ in (12). Assuming isotropy then yields, which is a standard eddy viscosity model for the transport of φ, with the eddy viscosity given by ν φ = Ck c T φ . Here the subscript "c" indicates that k c is the kinetic energy associated with the convecting velocity fluctuations and the subscript "φ" that T φ is the decorrelation time for the fluctuations of the conserved quantity φ. In LES, we essentially extend this approximation to filtered quantities and write However, u i φ is now itself a fluctuating quantity, and there is no reason to expect that the fluctuation of the model on the right hand side of (14), which was formulated to represent the expected value, will match those of u i φ . Therefore, a subgrid model of the form (14) takes on an implicit dual role. The equality in (14) is while the model fluctuations, if correlated with fluctuations in the gradient of φ, are responsible for the dissipation φ 2 when φ is the momentum. We will proceed with the approximation of (14) and return to its expected value in due course.
The purpose of this seemingly banal recounting is to highlight a subtlety that is often overlooked when constructing subgrid stress models. In LES modeling, it is important to consider how the fluctuations being modeled are defined; that is, whether they are the standard fluctuations about the mean, as in the Reynolds decomposition (Reynolds fluctuations), or subgrid fluctuations defined relative to the resolved turbulence in the LES. When the fluctuations of u k and φ are defined in the same way (i.e. both Reynolds or both subgrid), there are natural choices for k c , T φ , and the gradient. For instance, if we were interested in the product of the total fluctuations, e.g. u i φ , the approximation is proportional to k tot T tot ∂ i φ where the "tot" subscript indicates values derived from the total fluctuations. On the other hand, were we concerned with the product of the subgrid fluctuations, e.g. u < i φ < , the approximation is proportional to k sgs T sgs ∂ i φ where the "sgs" subscript indicates quantities derived from subgrid scale fluctuations. Alternatively, for "mixed" terms, such as u i φ < , the three quantities forming the model, k c , T φ , and the gradient of φ, must also be mixed. In particular, k c is associated with the definition of the convecting velocity fluctuations, T φ is related to the definition of fluctuations in the transported quantity (φ), and the gradient is taken of the complement of the φ fluctuations.
For mixed scale terms, we therefore have and Let the level of resolved turbulence be measured by α = k sgs /k tot , assume T tot = k tot /ε and T sgs = k sgs /ε = αT tot , and define a total eddy viscosity as ν tot = Ck tot T tot . Then the possible combinations of filtered products would be approximated as and To apply these scaling arguments, we rearrange the unresolved stress contributions in (10) as The first two terms are transposes of each other so their sum is symmetric, as is of course the third term. However, a straight-forward application (17)(18)(19) with φ = u j does not yield a symmetric result for τ sgs . Instead, we take the forms (17)(18)(19) as guidance and symmetrize the velocity gradients that appear, Gathering like terms we obtain where ν s = α(2 − α)ν tot and ν e = α(1 − α)ν tot . Thus, from applying simple eddy viscosity arguments to the total subgrid stress tensor, we see that the model-split form arises naturally. Further, this analysis indicates the contribution to the subgrid stress from the mean and fluctuating gradients do not scale in the same way, i.e. the form τ sgs ij − 1 3 τ sgs kk δ ij = 2ν t S ij is not generally valid. In expectation, (23) is simply and the contribution of the model (4) to the Reynolds stress depends only on the mean. On the other hand, the energy transfer provided by the model (5) arises entirely from the correlation of S > ij S > ij , scaled by 2ν e . It is quite apparent why this mixed scaling issue has been overlooked. For RANS, it is not an issue as all scales are the same. Subgrid stress models were patterned off RANS and adopted identical forms using length scales from the filter width, timescales from the inverse of the filtered velocity gradient, and the filtered velocity gradient itself. Therefore, such models are only directly relevant to the subgrid-subgrid term of (III B), i.e. u < i u < j . Further, for filters and projections which minimize the L 2 error of ||u > −u ||, the mixed terms are zero as u > i and u < i live in orthogonal subspaces of the full space of Navier-Stokes solutions. Many numerical method do not satisfy this requirement. For instance, the second-order finite volume method used in this work (more details in Sec. V) most closely approximates box filtering in physical space. With such common methods, the decomposed velocity fluctuations are not uncorrelated and the above arguments are relevant. This raises another practical issue. The above arguments made use of α = u < i u < i / u j u j . In a simulation, we only have access to the information of the resolved field and an approximation of u j u j provided by the RANS model. Therefore, we directly have an approximation of only β = 1 − u > i u > i / u j u j or τ kk / u j u j . When the mixed terms are non-zero, α and β are not equivalent. We must then find a relationship between α and β for the particular numerics of the code in use.
In isotropic turbulence at infinite Reynolds number, the Kolmogorov inertial range expression for the two-point correlation leads to α ∝ β, with the proportionality constant depending on the details of the filter. However, we need to use this relationship for wall-bounded shear flows, and so here we estimate α as a function of β by applying box filters of varying size to channel flow fields at Re τ = 1000 [1], excluding the region y + < 300 since it is wall and viscous dominated. Averages were computed over the homogeneous spatial directions and six snapshots in time separated by two channel flow throughs. The values of α and β determined at various y locations and for six different filters sizes were aggregated to determine a single relationship between α and β. In Fig. 1, to reduce statistical noise, the data were binned by α value to obtain average β's. We find a simple relation of α = β 1.7 as shown in Figure 1(a). Clearly, given the difference between the results for isotropic turbulence and channel flow, the relationship between α and β depends on and (b) ratio of total subgrid shear stress (including mixed terms) and total shear stress, overlaid with α(2 − α). characteristics of the turbulence. It will be useful to generalize this relationship to reflect that, but for now, we have used the relationship inferred from the channel flow DNS.
Of interest to the general eddy viscosity arguments presented here, we can also examine the total subgrid stress scaling. Figure 1(b) shows this relation for the box filtered channel shear stress data. The α(2 − α) scaling determined from the arguments above is consistent with the data and is used in the models discussed here, though α 9/8 (2 − α) fits the data slightly better.

IV. ACTIVE MODEL-SPLIT HRL
The hybridization strategy proposed here is motivated by deficiencies of previous methods, as described in Sec. II. To overcome these problems, as described in Sec. III, the total model stress is split into two separately modeled components, and an appropriate scaling of the mixed stress terms is incorporated. Further, a stirring force is used to introduce resolved fluctuations where they are capable of being resolved. In this "active model-split" (AMS) hybridization, the resolved momentum equation takes the following form: where u i is the resolved velocity, p is the resolved pressure, τ s ij is the "mean" subgrid stress, τ e ij is the "fluctuating" subgrid stress which results in energy transfer from the resolved to unresolved scales, and F i is the forcing. The forcing term subsumes all filter-width commutation errors.The models for τ s ij and τ e ij will be formulated in terms of "pseudo-mean" quantities (thus the quotes above), which are described in the following subsection before the detailed description of the hybrid formulation.

A. The Pseudo-Mean
The arguments presented in Sec. III have tacitly assumed knowledge of the expected value for all required quantities. In particular, u i is used to define u > i and then α and β through the resolved TKE as 2k res = u > i u > i . Additional averaged terms will appear in Sec. IV B and Sec. IV D. Yet, this information is not readily available in LES or HRL, necessitating a method for obtaining mean quantities. The natural assumption would be that the method for approximating the mean should be as close to the true expected value as possible, perhaps even warranting solution of a separate set of RANS mean velocity equations as in [18]. However, it was determined emperically that using the true mean in the AMS model formulation yields poor simulation results. The difficulty arises because of the use of an eddy viscosity model to represent the mixed term u < i u j (23), which represents the transport of the unresolved fluctuations, u < i , by the total fluctuations, u j . The total fluctuations include resolved fluctuations with much larger length and time scales than the unresolved scales. Interactions between these very large resolved scales and the unresolved scales are poorly described by dissipative processes, as the dominant effect is the large-scale advection of the unresolved scales. That is, these slowly varying and relatively energetic fluctuations behave more like an unsteady mean in the context of an eddy viscosity model (EVM). For this reason, using a pseudo-mean that filters out all but the largest and slowest scales of motion when defining the resolved fluctuations leads to a more accurate representation of the mean stress and a more consistent definition of α for the purpose of stress modeling.
Anticipating use of the AMS formulation in complex geometries in which the turbulence is inhomogeneous, time averaging is used to define the pseudo-mean, rather than spatial averaging. For stationary flows and flows with slowly evolving means a causal time average with an exponentially decaying kernel should be sufficient, and this approach is used here. In this case, the pseudo-mean {φ} of a quantity φ evolves according to where the time constant is proportional to the large-scale turbulent time scale, or T ave = C φ k tot /ε. For all averaged quantities C φ = 1, with the exception of P k , the production of turbulent kinetic energy, due to special circumstances as described in Sec. IV B. The performance of the model appears somewhat insensitive to the precise value of C φ , but no attempt has been made to optimize its value or to determine whether it needs to vary depending on flow characteristics. Using (25) with C φ of order one, eliminates all but the lowest frequency fluctuations from the definition of {φ}.

B. Subgrid Mean Stress Model
Following the argument presented in Sec. III B, the deviatoric subgrid contribution to the unresolved Reynolds stress tensor is assumed to scale with τ s ij ≈ α(2 − α)τ R ij where α is computed from the psuedo-mean resolved turbulent kinetic energy as α = k sgs /k tot = β 1.
denotes the Reynolds stress determined by a RANS closure evaluated using only pseudo-mean quantities, to be consistent with the arguments outlined in Sec. II D.
Due to the presence of α in the subgrid stress model described above, it is necessary that k tot be known. This can be obtained from the k-equation of any two-equation RANS model, provided the resulting k is intended to be representative of k tot . For example "k" in the SST model [72] does not represent the turbulent kinetic energy in attached boundary layers, but rather is more akin to 3 2 v 2 where v 2 is the wallnormal component of the Reynolds stress tensor. Use of this "k" would be inconsistent with the current evaluation of the resolved TKE and the AMS approach would perform poorly near walls. Recall that the velocity fluctuations u i = u i − {u i } is defined in terms of {·}, which is essentially a low-pass temporal filter rather than the expected value. Therefore, u i excludes the low frequency fluctuations u i = {u i } − u i . For consistency in the definition of α, k tot is defined as u i u i /2, i.e. excluding u i contributions. In determining a model equation for k tot , we assume that cross terms that would be zero if {·} were the true expected value are negligible and arrive at the usual k equation, with the exception of the convection and production terms (C k and P k respectively), which are Thus, the modeled k-equation from a RANS model can be used, with the exception of the production and convection terms. In the production term, we use the model for τ s ij to account for the unresolved portion of u i u j and approximate the expected value by the pseudo mean (25) with C φ = 4 (signified as [·]). This yields For simplicity u i in the convection term is approximated as {u i }, though approximating it as [u i ] would be more consistent. This simplification did not appear to make a significant difference.

C. Energy Transfer Model
As discussed in Sec. III A, when τ s ij is not fluctuating, the transfer of energy from resolved to unresolved turbulence is entirely due to τ e ij . Traditional LES SGS models primarily function to model this transfer of energy with the majority of the Reynolds stress assumed to be resolved [59]. Accordingly, one can use typical LES SGS model forms for τ e ij in addition to the scaling argument-based model presented in Sec. III B. In principle, any standard SGS model can be applied. However, for typical eddy viscosity model forms (e.g. Smagorinsky and variants), τ e ij will contribute to the mean stress due to non-vanishing correlations between fluctuations in the model viscosity and the velocity gradient tensor. In this case, these contributions to the mean stress must be accounted for in the subgrid Reynolds stress formulation described in Sec. III B.
On the other hand, if the eddy viscosity does not fluctuate or if its fluctuations are uncorrelated with the fluctuating velocity gradient, then τ e ij = 0 so that τ s ij and τ e ij are the mean and fluctuating parts of the model stress. The AMS formulation can thus be considered a technique to segregate the mean and fluctuating model stresses so they can be modeled differently, consistent with their differing roles of representing mean momentum transport and the transfer of energy from resolved to unresolved turbulence.
Here, we use the tensor eddy viscosity M43 model, which is formulated to account for anisotropy of the LES resolution [73]. Its coefficients were determined for consistency with second-order finite volume numerics as described in Appendix A of [73], with the resulting coefficients provided in Appendix B here. Details of the model formulation are provided in [73] but a modification is required here. As derived, M43 is an algebraic subgrid model which dissipates at the anisotropic resolution scales. It is formulated in terms of the mean dissipation rate, which in [73] was prescribed a priori from knowledge of the forcing energy injection rate in simple forced homogeneous isotropic turbulence. Here, the mean dissipation is determined from the dissipation equation in the RANS model. In this case, the M43 model becomes a two-equation (k and ε) tensor eddy viscosity model. Where a non-negligible portion of the total dissipation is resolved, is not true, the dissipation used in the M43 model should be reduced by the resolved dissipation. This modification would correctly observe DNS limits. However, when resolutions are relatively coarse, as is the case with the simulations performed here, the resolved dissipation may be neglected.

D. Active Forcing
As discussed in Sec. II A, resolved turbulence can be passively generated only by creating regions of modeled-stress depletion. Thus, relying on natural self-generation of turbulence fluctuations inherently requires reducing the fidelity of the model. When the resolved flow includes turbulent fluctuations the energy cascade will naturally populate additional finer scales of motion as they can be resolved, if the modeled dissipation is also reduced. However, the near-RANS regions will not self-generate finer turbulent scales without explicitly disrupting the mean stress. A high fidelity hybrid formulation therefore requires a mechanism for generating turbulent fluctuations, and here we propose a fluctuating body force (active forcing) driven by a resolution metric.
An active forcing formulation requires three ingredients: 1) identification of regions where more turbulence can be resolved, 2) determination of the rate at which resolved fluctuations should be introduced, and 3) specification of the structure of the generated velocity fluctuations. Where finer velocity fluctuations can be resolved, they should be introduced with the length scale of the smallest resolved turbulence so as not to disrupt existing larger, and presumably accurate, structures. Conversely, in under-resolved regions, τ e ij will naturally dissipate resolved energy. However, due to grid and flow inhomogeneity, τ e ij may not affect this transfer fast enough. Here, the forcing scheme is only formulated to add energy to regions that can resolve more turbulence. A forcing formulation capable of actively removing energy in under-resolved regions would also be useful. Currently, a modification to τ e ij , discussed in Appendix B, specifically (B4), partially alleviates the issue of under-resolution. Each of the three necessary forcing components is presented next.

Resolution evaluation
Identifying regions where forcing is needed requires evaluating the grid's ability to resolve more, or any, of the local turbulent fluctuations. Historically, scalar measures of grid resolution (e.g., cell diagonal or volume cube-root) have been compared to a scalar turbulent length scale to make this assessment. With anisotropic grids and/or turbulence, such measures are incomplete indicators and generally insufficient. Instead, the evaluation should be based on the locally least-resolved orientation to ensure adequate resolution in all directions. Naturally, this measure must depend on both the grid and the flow state.
First, the resolution capacity of the grid is described by the resolution tensor, M ij , [74] which characterizes resolution anisotropy. The eigenvalues, λ M i , of M ij are the resolution scale in i th principal directions, while its eigenvectors define those directions. Common grid measures are invariants or eigenvalues of M ij (e.g. δ vol = (det(M ij )) 1/3 ). Defining a resolution metric based on M ij provides a tensorially consistent mechanism to account for resolution anisotropy.
Second, the size and anisotropy of the largest unresolved turbulence scales, which are also the smallest resolved turbulence scales, can be characterized by the production of unresolved turbulence and unresolved kinetic energy. This is analogous to the expression of the integral scale as proportional to k 3/2 /ε in stationary, homogeneous, isotropic turbulence, since then the dissipation, ε, is also the rate of production. To represent the anisotropy of the unresolved length scale, we adopt the ansatz that the tensor representing the inverse unresolved length scale L −1 is given by In (31), the non-deviatoric part, 2 3 βk tot δ ij , is included because it affects the scale anisotropy of the production, despite the fact that it does not contribute to the production of kinetic energy, and is therefore often ignored. Note that in (28), the indices of R −3/2 are contracted with the velocity component indices of P, so that only length scale anisotropy is represented in L.
The ratio ρ(e) of the grid resolution to the smallest resolved scale in the direction defined by the unit vector e is then given by Clearly adequate representation of the resolved fluctuations requires that the maximum of this ratio over all directions be less than some constant, here taken to be one. Also, if the maximum ratio is significantly less than one, then smaller scale fluctuations can be resolved, indicating that forcing can be applied to generate such fluctuations.
However, the anisotropy of the unresolved fluctuation covariance R ij is commonly not available. In many RANS models, such anisotropy is only accounted for near the wall, often through the use of wall functions. As discussed by Durbin [75], such wall functions represent the effects of near wall anisotropy, and in the context of the v 2 -f model can be seen to scale with the ratio of the wall normal velocity variance (v 2 ) to the turbulent kinetic energy. Regardless of the two-equation RANS model being used, we can thus take the v 2 -f eddy viscosity relation as an expression for v 2 and determine the anisotropy measure ζ as: where C ζ = 7.5, as determined from the coefficients of the v 2 -f model, and T = k/ε is a turbulent time scale. Since v 2 is an estimate of the minimum eigenvalue of the Reynolds stress, and because the unresolved velocity covariance is expected to be less anisotropic than the Reynolds stress, the following inequality should hold: where the maximum on the right hand side refers to the maximum eigenvalue of its tensor argument. The coefficient C r is related to how many grid spacings are required to resolve the structures produced at the model cutoff. With the exception of the study described in Sec. V C we take this coefficient to be unity. The quantity r M is easily determined, and (34) implies that using it as a resolution indicator is conservative in the sense that it will indicate that scales resolvable on the grid are larger than it can actually support. This is the resolution indicator used here. In particular, when r M < 1, the grid is locally capable of resolving smaller scale fluctuations. When the grid resolution is isotropic and the turbulence is isotropic and in equilibrium, r M reduces to a length scale comparison similar to that used by DES [17]. Thus, it can be viewed as an anisotropic generalization of the DES resolution metric.

Forcing structure
In regions where r M < 1, the goal of forcing is to gradually introduce resolved fluctuation, thereby allowing the forced resolved structures to evolve into actual turbulence without corrupting the mean. Uniform forcing over some region where r M < 1 is clearly not desirable, as it would necessarily result in forcing the mean. Instead, we must prescribe the spatial structure of the acceleration field. Ideally, this artificial structure would retain turbulence characteristics, such as spatial correlation and intensity, of the largest of the unresolved turbulence. As more turbulence is added to the resolved field, the characteristics of the largest of the unresolved scales necessarily changes. The construction of this field as used here is ad hoc. It is presented here as an example, which was implemented to assess the utility of the greater AMS HRL framework. Improvements are certainly possible and desirable.
An artificial turbulence-like vortex field is defined based on the structure of a Taylor-Green (TG) vortex field with variable length scale: , where the magnitudes are somewhat arbitrary and selected such that h i h i ≤ 1 with A = 1, B = −1/3, and C = −2/3. The desired vortex sizes should closely mimic the local length scale at the implicit cutoff and be resolvable everywhere in the simulation domain. Thus, we select the vortex scale to be = min(N L L sgs , δ wall ), where L sgs = (βk tot ) 3/2 /ε, δ wall is the distance to the wall, and N L is a constant empirically set to 8. Applying to (35), we have where D (i) is the domain size in periodic directions, "nint" is the nearest integer, and subcripts in parenthesis indicate no summation is implied. This form is designed to ensure periodicity of the h's in periodic directions. The TG vortex coordinates, x p i , are specified to translate with the mean flow as to provide a degree of temporal correlation in addition to the desired length scale. With the prescribed structure in hand, we must now determine the appropriate magnitude for the forcing acceleration.

Forcing magnitude
As previously discussed, the magnitude of the forcing acceleration should be based on the largest of the unresolved fluctuations. Additionally, the rate of injection of manufactured fluctuations should be specified by the eddy turnover time of the smallest resolved structures so that the added fluctuations can be "healed" into realistic turbulence, as they are added, by the resolved turbulence. Making use of near-wall anisotropy, the target forcing is then defined as where C F is currently an empirical constant set to 8. Note the √ β in the denominator of (38) results from both the velocity and timescale corresponding to the subgrid scales. With the target acceleration magnitude (38) and structure (35) in hand, additional modifications are necessary to respect RANS and DNS limits. The RANS limit will be naturally detected by r M while identifying near-DNS conditions requires an approximation of the local Kolmogorov length scale. With the information provided by the RANS model, this limit is approximated as Limiters are incorporated in a scaling coefficient, η, as where D lim enforces the DNS limit and F r responds to the resolution evaluation and can thus be considered both a grid-resolved LES and RANS limiter. In the vicinity of either of these three limits, the scaling of h i is attenuated to prevent sharp transitions between full and no forcing. Gradual attenuation is prescribed in the following ad hoc functional forms.
Overall forcing behavior is controlled through F r as When r M < 1, forcing is activated with F r → 1. In the vicinity of r M ≈ 1, F r smoothly goes to zero. Finally, in RANS and under-resolved regions, r M > 1 and so F r = 0. The DNS limit is enforced with Specifically, as DNS resolution is approached, β → β kol ,β → 0, and D lim → F r . Thus, η also goes to zero, and the forcing effectively "shuts off" in the limit of DNS resolution. Alternatively, when β β kol , D lim ≈ 0 and the limiter does not affect the forcing. Note that (38) is poorly behaved in the limit of β → 0. To avoid this, one could replace β in (38) with max(β, β kol ). In the RANS limit, the structure of h i may be applied everywhere in the domain to affect turbulence generation. However, with existing turbulent fluctuations, applying a generic acceleration structure, such as h i , will on average result in both the addition and removal of energy. So, the magnitude of h i must be further modified. This issue is addressed by testing the resolved production due to h i as which will be used to clip (35). With these considerations, the forcing acceleration vector is The specification of F i along with the two model terms (Sec. IV B and Sec. IV C) and the pseudo-mean (Sec. IV A) closes the AMS approach. While there are many advantages of the AMS hybrid formulation, the model-split form introduces new challenges as well. Before presenting several tests of the new approach, we first briefly discuss these issues.

E. Discussion of AMS Formulation
The components of the AMS hybrid RANS/LES approach include: a RANS mean stress model, a subgrid scale dissipation model, under-utilized resolution evaluation, active fluctuation forcing, and computation of the resolved turbulent stress with a causal time average. The resolved mean stress is used throughout AMS, most notably to scale the RANS stress to the appropriate subgrid scale stress. All of these components are dynamically and continuously coupled in a single numerical simulation, on a single grid. There are two critical features of the AMS formulation that are necessary for good model performance. First, is the explicit separation of the mean stress and dissipation roles of the subgrid scale turbulence model. This splitting allows the eddy-viscosity models to be tailored to each role, enabling the overall AMS formulation to remain valid regardless of resolution. It is essentially an LES that by construction respects both RANS and DNS limits. In particular, in the limit of DNS, all the model terms (τ s , τ e and F ) go to zero, leaving only the numerical approximation of the Navier-Stokes equations.
The second critical feature is that the model hybridization responds only to the turbulence fluctuations that are resolved, rather than the fluctuations that could be resolved. This ensures that the modeled mean stress is always consistent with the resolved mean stress, which is important to avoid common problems in HRL, such a model stress depletion. Many HRL techniques promote the development of fluctuations where they can be resolved by reducing the modeled stress and thereby promoting instabilities, but while also disrupting the mean. In AMS, by contrast, resolved fluctuations are introduced explicitly through forcing when needed. The details of the forcing as described in Section IV D are necessarily ad hoc, as is the promotion of instabilities used in other approaches, since fluctuation information must be created. Using explicit forcing, however, provides an opportunity to design the forcing that will allow fluctuations to be introduced as rapidly as possible while ensuring that they realistically represent their contributions to the mean stress, as discussed below.
The model-split hybridization approach has numerous advantages over traditional HRL methods. First, following the argument in Sec. IV, blending, which is responsible for many of the difficulties in HRL, is obviated. Second, it uses the RANS model as designed. The turbulence model state variables represent mean features of the turbulence, and the governing PDEs depend only on mean, or psuedo-mean, quantities. Thus, pathological behaviors of the RANS transport models due to fluctuating state variables are avoided. Third, the RANS eddy viscosity appearing in τ s ij makes no contribution to the dissipation of resolved fluctuations, allowing turbulence at all resolvable scales to form naturally without being overly dissipated. The AMS formulation will also naturally be inactive in laminar flow regions. The RANS model that is part of the AMS formulation governs this behavior. When the RANS model predicts zero or negligible TKE as in a laminar region, the forcing will be zero and the eddy viscosity will also go to zero. Finally, nearly any combination of base mean and fluctuating models can be used. Because of this flexibility, advanced models are easily incorporated and used to treat complex flow and domain features such as the effects of resolution anisotropy with the M43 model. Since stress anisotropies are mostly carried by the largest of turbulent structures, it is reasonable to expect that even modestly resolving turbulent kinetic energy will result in significant improvement over RANS to the representation of total stress anisotropy. Nonetheless, the AMS approach does present the opportunity to include treatments for stress anisotropy such as explicit algebraic Reynolds stress models (EARSM) [76].
However, there are several issues demanding future attention beyond the scope of this paper. Perhaps the most glaring is the specifics of the psuedo-mean and precisely what large scale fluctuations should be excluded from the α scaling in the model. This concept, that some of the interactions with large-scale turbulence should not be modeled with gradient diffusion, is new to both HRL and LES in general. The degree of exclusion of low frequency turbulence from the model is effectively controlled by the coefficient used in the psuedo-mean. Another way to approach the modeled Reynolds stress would be to adopt the form of [41] and simply subtract the deviatoric portion of the mean resolved stress from the RANS model stress. This method would sidestep some of the ambiguity in RANS scalar terms but, as formulated, leads to no improvement over basic RANS behavior in expectation. Nonetheless, this should be explored. Further, the use of the causal time average results in a lagged response to large-scale unsteadiness and is strictly only valid for stationary flows. However, similar to the assumptions of unsteady RANS, if the unsteady timescales are much larger than the largest turbulence timescales and the averaging time, the averaging will produce an appropriate pseudo-mean for use in the AMS formulation.
The second main issue is the rudimentary form of the forcing. As discussed in Sec. IV D the prescribed structure of the synthetic forcing is necessarily ad hoc. The particular method outlined here has potential drawbacks. Foremost, the resulting F i is not divergence free. This is a result of TG fields only being divergence-free when constructed using a uniform overall scaling and vortex length scale and without clipping. In an incompressible solver, the dilatational portion is projected out, which alters the structure of the effective forcing. The result is that the forcing can remove resolved turbulence at some points in space and time. However, on average, the approach does result in net production of resolved fluctuations where desired. In a compressible solver, the method may result in spurious acoustic sources. Further, it should be possible to force more strongly, i.e. increase the coefficient C F (currently set to 8), if the forcing structure produced fluctuations that were more consistent with unresolved turbulence in equilibrium with the current u i . The forcing design described in Section IV D is certainly not the best that can be done in this regard. Improving forcing structure so that forcing strength can be increased would allow the size of the LES-resolved region upstream of an area of interest to be minimized. However, while the forcing used here is ad hoc, the AMS approach is not sensitive to the ad hoc characteristics of the forcing. That is, AMS responds to the presence of resolved turbulence and not to how the turbulence is introduced. As a result, any new or improved forcing technique can easily be substituted into AMS. Indeed, we regard the current forcing as a "place holder" to close the overall AMS modeling framework.
The RANS models used for τ s are improved relative to their use in a purely RANS computation through the inclusion of the resolved stress in the production term. However, for the sake of simplicity, we have not analogously included an explicitly computed resolved contribution ∂ j u > i u > i u > j to the turbulent transport of k (i.e. ∂ j u i u i u j ), relying instead on the standard eddy viscosity model to represent the entire term. This simplification can lead to local inconsistency between the resolved and modeled TKE, which we have observed as small regions of negative β where the total TKE is small. Such negative values were clipped to the Kolmogorov microscale value of β, given by β kol (39). The results of Section V suggest that, at least for the flows considered, this treatment is adequate. Yet we are essentially discarding useful information by not computing ∂ j u > i u > i u > j explicitly, and so we expect that doing so would improve the veracity and robustness of the model.
Finally, while we have touted model-flexibility as a benefit, we must also caution that different RANS models may require different stress scalings and a different resolution-adequacy coefficient. This is an unfortunate side-effect of how RANS models function through a series of error cancellations to arrive at a reasonable eddy viscosity for a given mean flow. That is, the scalar terms which adopt the names of "turbulent kinetic energy", "turbulent dissipation", etc. are not, strictly speaking, models for those terms. In particular the k from one model may be significantly different from the k from another.

V. AMS PERFORMANCE TESTS
The proposed AMS hybrid formulation is evaluated using an implementation in a branch of the finite volume incompressible Navier Stokes solver CDP v2.4 [77,78], developed at the Stanford Center for Turbulence Research. The solver is 2 nd order accurate in time and space with no upwinding used for convective fluxes in the momentum equation and time advancement performed with Crank-Nicolson. The convective term is linearized with an Adams-Bashforth prediction of the convective velocity. Two base RANS models are considered: Chien's k-ε [79] and the "code-friendly" version [80] of the v 2 -f RANS model [75]. Details of these models are provided in appendix A. The energy transfer model is the M43 model [73], which is specified in appendix B.
The results presented here include AMS simulations with varying resolutions. However, these should not be confused with convergence studies in the usual sense of numerical analysis. Indeed, even the finest resolutions reported here are very coarse by LES standards, with significant contributions of unresolved scales to the mean stresses, and there are always RANS regions near the walls. This is the resolution regime in which hybrid methods are intended to operate and standard SGS models are known to fail. Within this regime, we consider a sequence of AMS models, with different resolutions, to investigate sensitivities to the resolution. Models with finer resolution will resolve more of the turbulent fluctuations, relying less on RANS, and so one would generally expect them to yield more accurate simulations. However, this may not always be true for all quantities. We are most concerned that the mean velocity be accurately predicted, since its prediction is the most common objective for HRL in practice, with the expectation that HRL will improve on RANS predictions. For a HRL formulation to be useful it is important that its results not exhibit extreme sensitivities to the resolution, and that relevant features of the solution, particularly the mean, improve or at least not degrade as resolution increases. It is in this context that we examine AMS model results at different resolutions in the following subsections.

A. Periodic Channel
To demonstrate the potential of the proposed hybrid modeling formulation, fully-developed, incompressible, turbulent channel flow at Re τ ≈ 5200 is simulated. DNS data is available for this case [1] allowing detailed evaluation of the results. It is common for RANS models to perform very well for channel flow. Therefore, successful hybrid simulations for this case would simply not degrade RANS mean velocity profiles. This may seem like a rather modest goal, however, the relatively coarse levels of LES resolution used here place the simulations firmly in a regime where existing SGS models will fail.
Results are presented in the following three subsections with focus on 1) hybrid steady-state with varying resolutions, 2) hybrid state evolution in time, and 3) hybrid state evolution in space. The term "hybrid state" is used to indicate the level of resolved turbulence which, in general, will not be at a grid-resolved LES levels and will vary in both time and space. We emphasize that evaluation of HRL in the transition from RANS to the statistically stationary hybrid state, is absolutely critical for an HRL. While all but one of the flow scenarios considered in this paper employ periodic boundary conditions, in practice, a HRL will be required to transition in space from a RANS, even if only at inlet boundaries, to a hybrid LES-RANS state. If this transition behavior is corrupted by the HRL formulation, the downstream solution may be corrupted.

Stationary channel with varying resolution
The domain for the hybrid simulation is 8πδ × 3πδ × 2δ, where δ is the channel half-width; i.e. identical to the DNS domain in [1]. The base RANS model used is Chien's k-ε [79] in which u τ is specified a priori for use in the wall functions based on DNS. A spatially uniform streamwise body force (mean pressure gradient) was applied throughout the channel, which varied in time to maintain a constant bulk velocity. Multiple grid resolutions are considered as shown in Table I. For all resolutions, the wall-normal grid spacing is fixed with N y = 110, ∆ + y (wall) ≈ 1, and ∆ + y (center) ≈ 345 while the spanwise and streamwise spacing ranges from approximately 700 − 1500 wall units. Note that the finest resolution considered here is still some 70  times too coarse to approach the DNS limit for even Re τ = 180 [1]. Results are obtained using τ s ij scaling of α(2 − α) as described in Sec. III B. The hybrid simulation is initialized from a steady-state RANS solution. Snapshots of the resolved streamwise velocity fields are shown in Fig. 2 to illustrate the difference in resolved turbulence for each resolution case. Statistically stationary mean velocity and resolution level (β) profiles  Table 1.
are shown in Fig. 3a. As expected, the basic k-ε RANS model does a good job of reproducing the mean velocity profile with only a slight under-prediction in the buffer layer. For all resolutions levels, the AMS hybrid formulation also adequately predicts the mean flow with a very small deviation from the RANS and DNS solution in the center of the channel. As qualitatively illustrated in Fig. 2, and quantitatively shown by β in Fig. 3, the AMS formulation allows the hybrid state to evolve to an appropriate level of resolved turbulence, in the sense that the mean is correctly predicted, for the different grids in each simulation. With each increase in resolution in Table 1, turbulence is resolved further into the log-layer while introducing no log-layer mismatch. Naturally, the most resolved turbulence occurs in the fine case with β as low as 0.4 at the channel center. For typical LES, β 0.1 is expected over the entire domain for the model to perform well [59]. The model-splitting formulation presented here retains good performance well past this threshold, enabling true "coarse" LES. The nearly identical mean profiles for such varied levels of resolved turbulence is an indication of the robustness of the formulation.
The response of the underlying RANS model to the modified production (Eq. 27) is shown Fig. 3b. Here, k determined from the RANS transport equation in the hybrid simulations remains virtually unchanged from that in the pure RANS model, despite the modified production terms and the use of the pseudo-mean velocity in the RANS equations. This confirms that the formulation is successful in preserving the RANS solution characteristics in the presence of resolved fluctuations. Ostensibly, it may seem that the RANS k should move towards the DNS value with increasing resolution; however, RANS models function through a series or error cancelling. That is, they are tuned to produce an appropriate eddy viscosity and not the correct individual scalar quantities whose name they bare (k, ε, etc.). Thus, changes to the RANS k might actually break the underlying RANS behavior. Both the RANS k and the turbulent kinetic energy including the resolved fluctuations ( 1 2 u > i u > i + βk) remain virtually unchanged. With increasing resolution, the two measures of TKE do increase towards the DNS values for y + > 300 due to increased resolved kinetic energy and resolved production. The fact that these two definitions of TKE are nearly indistinguishable indicates consistency between the resolved a modeled turbulence. However, resolved TKE must ultimately approach the DNS as α → 0, that is as all the fluctuations are resolved. Therefore, the modeled TKE and resolved TKE may separate with reducing α, which begins to be visible in the finest resolution case in Fig 3b. Such a separation is due to deficiencies in the underlying RANS model.
By changing the spanwise and streamwise grid spacing, we are also varying the resolution anisotropy, from 2:1 cell aspect ratios at the center of the channel in the fine case to 20:1 at y + = 400 in the extra coarse case. The converged steady-state mean profiles for these different degrees of resolution anisotropy indicate that the resolution adequacy parameter proposed in Sec. IV D successfully quantifies the turbulence that an anisotropic grid is capable of resolving in the presence of inhomogeneous mean shear. Further, the use of the M43 model for the energy transfer portion of the model, and its scaling with the pseudo-mean ε obtained from the RANS transport equations, appears to be valid.

Channel with temporally evolving hybrid-state
The results in the previous section are for stationary hybrid states. In this section, we examine the temporal development of the hybrid state for the case of the "coarse" resolution. The forcing method presented in Sec. IV D was formulated with the goal that that hybrid solutions remain valid through transition from RANS to an arbitrary LES state. For stationary flows, the transient from an initial field to the final solution is of no consequence. However, for unsteady flows, transition through varying turbulence states in time is of prime concern and may continue throughout the duration of a simulation. The case examined here represents an extreme example of this situation, i.e. transition from no resolved turbulence to a gridresolved LES. Such a drastic transition would not actually occur in most unsteady simulations. The opposite transition from less to more modeled turbulence is also of interest, but is not considered here.
As previously discussed, the hybridization is driven by forcing. If no forcing is applied, the simulation would remain a RANS simulation. The process of evaluating the resolution and introducing resolved turbulence is summarized through the resolution adequacy parameter (34), forcing field (45), and resolved velocity for the coarse case in Fig. 4. An example of the wall-parallel forcing field is shown in Fig. 6. Note that by construction, a variable length scale and clipped Taylor-Green field is not divergence-free. For an incompressible solver, the divergence of the forcing field is projected out in the pressure solution step so overall behavior of the hybrid method is not affected by forcing divergence. Thus, for incompressible applications, an additional pressure-Poisson-like solve is not necessary however, it may be important for compressible solvers. At t f = 0, forcing is activated, and the state evolves in time from the RANS initial condition to a Numbers in the legend refer to x/δ, the down-stream distance normalized by channel half-width.
statistically stationary hybrid state as shown in the previous section. During the first few steps of forcing (Fig.4, t f = 0.0 + ) a large region with excess resolution is present, indicated by r M < 1, and shown as blue in the figure. Note that r M is much larger than 1.5 near the wall in regions which are incapable of resolving turbulence and remain RANS. Coherent, large fluctuations are excited as the scale of the forcing is the large scale turbulent lengthscale. Some clipping of the forcing (44) still occurs due to the approximation of means through time-averaging. The local magnitude of the forcing follows the RANS k profile and the hybrid solution is still nearly identical to the RANS state at this point. The forcing field is characterized by very long streamwise, but thin wall-normal, structures. Shortly after (t f = 0.4), some fluctuations become evident and the forcing field is drastically altered with much shorter streamwise structures. This is due to reduction of the subgrid turbulence length scale (L sgs ∼ (βk tot ) 3/2 /ε) and as a result, the forcing lengthscale, and to more local clipping where the prescribed Taylor-Green field would remove energy. By t f = 0.8, distinct turbulent structures are visible and the over-resolved core has been reduced and flanked by bands of r M ≈ 1 extending towards the center of the channel. At t f = 1.6, we begin to see a shortcoming in the current formulation. A grid-resolved LES region extends through the entire channel (green regions of r M ≈ 1) as intended, however, there are still small region of slight active forcing and bands of under-resolved regions (red regions of r M > 1). This is another artifact of time averaging being used as a surrogate for the true expected value of r M . That is, small fluctuations in the puesdo-mean of r M result in local forcing activation when r M is only slightly less than unity whereas there is no counterpart when r M is slightly greater than one. The result, on average, is that excess energy is being continually added to the small scales. An ad hoc modification to the M43 model coefficient to remove this excess energy is described in Appendix B (see Eq. B4), but, it appears to do so too slowly. The effects of this under-resolution may contribute to the temporary disruption of the mean velocity profiles in Fig. 5a. Keeping in mind that the steady-state β for the coarse mesh only achieves a minimum of about 0.6 ( Fig. 3a), we see the resolved field over-shoots this value to 0.4 at t f = 1.6. Thus, as indicated by the red r M > 1 bands in Fig. 4, more resolved turbulence has been added than can be properly resolved. The blunting of the mean velocity profile in the center of the channel may be a result. After this over-shooting, the field gradually heals with the M43 modification removing excess energy. By t f = 8, the mean velocity has nearly attained its steady-state Sec. V A 1. tot , at multiple wall-parallel planes a t = 0.04 after initiating forcing from the steady RANS solution. Entire spanwise domain is displayed. Early forcing experiences no clipping yielding large coherent acceleration with the size of the structures decreasing with proximity to the wall. Note that ktot is dependent on wall distance.
Though we have not performed a rigorous spectral analysis of the resolved velocity field, it appears that the largest scales of turbulence are not excited until later in the simulation. For instance, compare the streamwise oscillations at t = 0.4 and t = 8.0 snapshots of Fig. 4. Contrary to the intended forcing behavior, the largest turbulence scale should be excited at the start of the simulation, suggesting a potential deficiency in the current forcing structure. This was alluded to in the previous remark about the wall-normal forcing structure being too small. More interestingly, because the forcing is apparently not exciting such large scale motions directly, a transfer of energy to these large scales of motion is necessary. While some of the energy may be coming from the excited turbulence scales, it may also be transferred directly from the mean. A drain of energy from the mean would also be consistent with the temporary blunting of the velocity profiles at early time. The mean then recovers through the action of the applied mean pressure gradient. The structure of the forcing at early time is plotted in several wall-parallel planes in Fig. 6 in addition to the wall-normal plane shown in Fig. 5. The prescribed forcing field does result in large structures away from the wall, but only in the wall-parallel planes. Such granularity in the wall-normal direction forcing certainly contributes to the largest turbulence structures not being excited. Correcting the prescribed early forcing structure would likely help reduce the small distortion of the mean during the forcing.
A potential remedy to the discussed under-resolution would be to allow the forcing field to also remove energy in regions satisfying the conditions α < 1, r M > 1, and P test F < 1 (44). Alternatively, and certainly more simply, the coefficient C F (38) could be reduced so that the forcing is more gradual and there is no healing time necessary to correct the overshoot. Drastic transitions in time from RANS to grid-resolved LES would not actually occur in a simulation other than during the initial startup or in some highly unsteady problem. Therefore, such a coefficient reduction may not be generally necessary for practical application of the AMS method. This viewpoint is supported by the results of the next section.

Channel with spatially evolving hybrid-state
Next we consider the case of a hybrid-state developing in space. Unlike the previous case of transition from RANS to LES in time, transition in space will occur in many practical HRL simulations. That is, one would like to construct a mesh with LES-quality resolution only where needed, specify only mean inlet velocity and scalar profiles with no resolved turbulence, and rely on the HRL framework to transition from RANS to LES without any loss of simulation fidelity. Poor transitional behavior may corrupt the entire solution downstream of the transition region. Therefore, spatial transition is important to the stationary state and needs to be tested.
In this section, rather than using periodicity in the streamwise direction, a RANS inlet profile is specified. Again using the resolution of the coarse case, the streamwise length of the domain is extended to 16π while maintaining the same span and height. The inlet and outlet sections of the extended channel are shown in Fig. 7 and mean velocity profiles, averaged over specific streamwise slices, are presented in Fig. 3b. Contrary to the temporally developing case, no β overshoot and no LES regions of r M > 1 are present. By the end of the channel, the mean velocity and β have nearly reached the coarse grid stationary hybridstate. There is still a small over-resolved core-region at the outlet of the channel. Further, the mid channel velocity is slightly lower than the DNS and periodic hybrid simulation. Forcing is most active just at the inlet (x = 0) but it is sustained at low amplitude over the entire channel length. This early region, x < 1, produces the largest relative drop in β. Such a rapid addition in resolved turbulence without disruption to the mean is highly desirable for hybrid applications as it indicates that only small regions of LES resolution upstream of a flow feature of interest are necessary.
Unfortunately, the relatively slow drop in β after x = π/2 indicate the current forcing formulation is not able to rapidly fill in all resolvable scales. This behavior seems to be particular to the spatial development case. For instance, after t = 1.6 in the time-developing channel, β has dropped below 0.5. Based on the bulk velocity, the fluid has traveled just over 12π at this time. However, after traveling 12π from the RANS inlet in the spatial developing case, β has only dropped to just under 0.7. Thus, it seems there are subtleties to the current forcing structure that relies on periodicity in the streamwise direction to be most effective. Nonetheless, the AMS approach is successful in transitioning a hybrid state from RANS to LES in space.
This also indicates that the current C F value is not excessive for practical applications of AMS.

High Reynolds number channel
To conclude our study of basic AMS in channel flow, we consider a much higher Reynolds number, Re τ = 20, 000. We have also exercised AMS in lower Reynolds number cases (Re τ = 180 and 2000, not shown) with results not substantively different than the Re τ = 5200 cases discussed in the previous sections. However, at higher Reynolds number, we can expect the outer kinetic energy peak to increase due to additional large scale turbulent structures away from the wall [81].In general, RANS models do not represent this outer peak well, yet they can still produce excellent mean velocity profiles with a well-defined log-layer. This is a result of the previously mentioned error-cancellation in RANS models. As indicated in Fig. 3, the AMS hybrid formulation improves the prediction of the outer k-peak, and one might suspect that this improvement could disrupt the error cancellation in the RANS representation of the subgrid Reynolds stress used in AMS and thereby affect the mean. At Re τ = 5K and the resolutions considered, this improvement is small and the overall performance of AMS is still excellent. The largest structures that are responsible for the increase in the outer kinetic energy peak with increasing with Reynolds number, should be easily resolved in AMS. So, at higher Reynolds number, AMS will produce larger improvements in the outer kinetic energy peak, making a higher Reynolds number channel simulation a simple test of whether improving the representation of k in AMS can degrade the performance of the RANS-based subgrid Reynolds stress model. To this end, an AMS simulation was performed with the fine resolution of Table I at Re τ = 20K but with the first wall-normal grid point shifted to maintain y + = 1 for the higher Re τ . Therefore, the wall-normal resolution in the center of the channel is more coarse than the Re τ = 5200 fine resolution case. Chien's k-ε model was found to perform poorly at this Reynolds number in the near-wall to buffer region. Since we do not expect AMS to improve the simulation this close to the wall with such a grid, it was necessary to use a different base RANS model for this AMS test. While there are no DNS or detailed experimental results available at Re τ = 20K, the v 2 -f model yields expected mean velocity and is used with AMS here. Mean velocity and k profiles from the Re τ = 20K simulation are shown in Fig. 8. There is a depature from the log profile from y + = 200 to y + = 2000. In comparison with Fig. 3, β is lower in the center of the channel and values below one extend to lower values of y for the same resolution. This is a direct result of large turbulence scales which are not present at lower Reynolds number and are resolvable at the current resolution. The largest deviation from the baseline RANS profiles is in the outer-layer peak in k. This is so for both k obtained directly from the hybrid RANS k equation and for k determined from the resolved fluctuations plus α times the hybrid RANS k. While we cannot directly confirm how well the AMS outer peak corresponds to the true peak, the existence of the peak is already an improvement over the basic RANS model. As suspected, the increase in k relative to RANS affects the RANS model for the subgrid Reynolds stress, resulting in a distortion of the mean. Note that the departure starts where k begins to increase relative to baseline RANS and ends at about the location of the outer peak. The problem with the RANS model for the Reynolds stress is that the large scales responsible for the outer layer peak in k and for the majority of the Reynolds stress away from the wall also extend closer to the wall, where they do not carry significant Reynolds stress [81]. As a consequence, in this intermediate region where the AMS mean velocity departs from expectations, k is not a good predictor of the Reynolds stress. In baseline RANS this is mitigated by the compensating error that the outer layer peak in k is absent. In AMS, the improved representation of k will have two consequences: the hybrid RANS estimate of the turbulent viscosity will be inconsistent; and β will not adequately represent the contribution of the resolved turbulence to mean momentum transfer.
These results highlight the sensitivities of the AMS formulation to the fidelity of the hybridized RANS model. To improve this fidelity, it may be fruitful to further modify the hybrid RANS model to make more use of knowledge of the resolved scales. Another possibility is to reformulate the hybrid RANS model to represent only the subgrid turbulence or to enrich the representation of the subgrid Reynolds stress and its dependence on β. These are out of scope for the current paper, but are important directions for future development.

B. Periodic Hill
In the previous channel cases, RANS models perform well, and the simulations demonstrate that the AMS hybrid framework does not disrupt this performance as resolution varies and hybrid states evolve. In this section, we move to a case where RANS is known to fail and evaluate whether the AMS formulation can improve over the RANS results. Flow separation and smooth-wall reattachment are of great interest in many engineering applications. However, RANS models tend to incorrectly predict the recirculation region, often delaying reattachment. This leads to erroneous prediction of critical quantities of interest such as lift and drag. With localized regions of model deficiency, such flow scenarios are precisely where HRL should be of greatest benefit. The periodic hill test case geometry consists of a channel with a 2D-hill which produces a separation region of approximately three hill heights. Having been studied both experimentally [2,82] and numerically [2], this case offers a rare combination of "truth" data as well as established resolution requirements for accurate simulation of the flow with existing LES models. The periodic hill is therefore an ideal test case to assess the utility of the AMS framework. In what follows, we therefore present periodic hill simulation results using both a baseline RANS model and an AMS hybrid formulation.
Flow over a periodic hill geometry at Re h = 10, 600 (where h is the hill height), case UFR 3-30 (available: http://qnet-ercoftac.cfms.org.uk/w/index.php/UFR_3-30_Test_Case) in the ERCOFTAC database, is simulated. The computational domain consists of a single peak-to-peak period of the periodic hill in the streamwise direction of dimensions 9.5h × 3.035h × 4.5h. Periodic boundary conditions are applied in both the streamwise and spanwise directions. The domain is discretized as shown in Table II with plus units based on the friction velocity at the peak of the hill, or at x = 0. Simulation grids consist of approximately 0.9M, 0.5M, and 0.2M cells in total. The wall-resolved LES used for comparison [2] uses over 13M cells. As in the channel, a streamwise body force (mean pressure gradient) is required to maintain a constant bulk velocity. In this case it is applied only in the upper portion of the domain (y > 1). This avoids a non-physical interaction of the forcing with the inclined hill surfaces. The base RANS model used for the periodic hill is a slightly modified version of the "code-friendly" [80] v 2 -f model [75] (appendix A, Eq. A3). Elliptic  x is an average value as streamwise clustering in the recirculation region is used on all grids. Grid size reductions are reported in the last column as the ratio of the WRLES [2] grid size to that of the AMS simulations.
relaxation RANS models have the advantage of not requiring wall distance or wall shear stress, making them attractive in simulations of complex domains with unstructured grids. The low-Re k-ω model [83] would also be a potential candidate for AMS application to complex geometries as it satisfies the requirements of providing both a turbulent length and time scale and the model k captures near-wall and log-layer behavior. Contrary to the channel results, where profiles were simply averaged over planes at a fixed y-location, profiles  Table II. Instantaneous AMS snapshots are shown here to illustrate the degree of resolved turbulence, mean velocity profiles for AMS are shown in Fig. 10.
presented here are obtained by averaging quantities (Sec. IV) using 50 samples over 20 flow-throughs. The pseudo-mean was found to be slowly fluctuating even after the simulation was brought to a quasi-steady state after 20 flow-throughs. Further, the resolved stress as defined through the psuedo-mean excludes lowfrequency structures (see Sec. IV B). To enable calculation of the total Reynolds stress to compare with experiments, the resolved stress is calculated directly through u i u j − u i u j . Again, the simulation was initialized from the RANS solution. Figure 9 shows streamwise velocity contours for both the baseline RANS v 2 -f model and an instantaneous field from hybrid simulations for all grids, illustrating the scales of resolved structures in the hybrid simulation. Comparisons of experimental, LES, RANS, and AMS mean velocity profiles at various streamwise locations are displayed in Fig. 10. Clearly, AMS improves RANS predictions, matching experiments and WRLES quite well for all resolutions considered. Of primary importance to aerodynamic applications is the  Table II, along with data from the experiments of [82] and WRLES of [2]. Note that the vertical velocity is multiplied by a factor of four to enhance visibility.
reattachment point at x ≈ 4. A gradual improvement of the reattachment length is obtained with increasing resolution for AMS while the RANS reattachment is delayed downstream to approximately x = 7. There is little difference in reattachment length between the medium and fine resolutions, indicating consistency of the hybrid solutions once sufficiently high resolution is achieved. In addition to the reattachment, the hybrid simulations also improves the streamwise velocity profile nearly everywhere in the domain. In particular, a drastic improvement is seen in the near-wall region at the top of the hill and everywhere after the reattachment point on the bottom wall. The fine resolution considered here yields further improvements near the wall and just past the reattachment point. Similar improvements to the mean vertical velocity profiles with increasing resolution are also observed. Though the vertical velocity is improved over nearly all the streamwise locations, at x = 0.5 and 1.0, the hybrid results are nearly identical to RANS. The likely reason for this is apparent from examination of the hybrid state in this region. Spatial distributions of the fraction of unresolved turbulence energy and the resolution adequacy param- eter are shown in Fig. 11 for the fine grid case. As expected, β = 1.0 in the RANS regions near the walls.
In the large hybrid regions throughout most of the domain β ≈ 0.5 with a minimum of about 0.3 near the top of the channel. In the recirculation region, β increases to about 0.7 and very nearly unity in the actual separation shear layer. This is a clear indication that the fine resolution used here is not sufficient to resolve the fine near-wall structures that propagate into the shear layer. This is the likely cause of the poor vertical velocity predictions in this region. However, the resolution used here is sufficient to improve representation of all other flow characteristic, in particular the separation, indicating resolving fairly large structures is all that is generally necessary for good prediction of this flow. That is, representing the mixing of high-momentum upper channel fluid with the recirculation region requires only coarse resolutions.
As expected, r M is very high in RANS regions near the wall. Interestingly, r M varies significantly in hybrid regions from just below unity in the recirculation regions, to approximately unity above the hill, to well above unity (r M ≈ 4) in the center of the channel. As discussed in Sec. IV D, hybrid regions of r M > 1 indicate locally under-resolved turbulent structures. Perhaps a combination of the resolved structures in this region being more three-dimensional than those around the hill, coupled with the fact that the spanwise grid spacing is coarser than the streamwise and wall-normal resolution, lead to this region of under-resolution. However, this does indicate that the modification to the M43 model (see Appendix B, Eq. B4) is not generally sufficient to remove resolved turbulence as it is convected into under-resolved regions. This suggests that introducing additional terms to the resolution adequacy parameter that are sensitive to the convective gradient of β or r M , may be necessary. Fortunately, the mean flow seems to be tolerant of these under-resolved regions, as we see no local disagreement between the hybrid results and experiments in the mean velocity profiles. Perhaps, this insensitivity is indicative of C r in (34) being overly restrictive. We will evaluate this hypothesis in the next section. The local over-resolution in the recirculation region leads to continuous local forcing of resolved structures. The observed distribution of r M indicates that, for complex geometries and flows, the convection of resolved turbulence leads to a continuous need to locally remove excess resolved turbulence in some regions while also having to add resolved fluctuations in others.  Table II, along with data from the experiments of [82] and WRLES of [2]. variance (top) vertical velocity variance (middle) and Reynolds shear stress (bottom). The prime here indicates the total fluctuations, that is, using fluctuations relative to the true mean. Note amplification by a factor of 10 for visibility. Despite the minor discrepancies described above, the AMS formulation yields nearly identical results to the wall-resolved LES of [2] at a resolution reduction of over an order-of-magnitude and computational costs reduced by a factor of 30 or more due to larger timesteps. These improvements are a direct result of the model-split formulation being capable of providing improved modeled stress predictions for arbitrary levels of resolved turbulence. This claim is supported by examining components of the Reynolds stress tensor in Fig. 12. AMS stress profiles generally show large improvements over RANS. However, some hybrid stress components become excessive. For instance, the fine resolution results in excesses of u u and v v in the recirculation region. The aforementioned sustained local addition of artificial fluctuations in the recirculation region may contribute to these disagreements. Another possibility is the improved production used in the RANS model (27) which enhances production of TKE in the separated shear layer from x = 0 to x = 1. While correct, in the context of the RANS system and its numerous modeling limitations this may cause an excess in TKE and the observed excesses in both normal stress components in the vicinity of the shear layer at x = 2, which then diffuses downstream. This is another example of the limitations imposed by the RANS model on the fidelity of the TKE. Shear stress, u v , is uniformly improved with the fine resolution only deviating from the experimental values in the separated shear layer just after separation at x = 1.
Of course, the periodic hill has been simulated successfully by several other hybrid and bridging frameworks. Successful simulations have been performed using DDES [84], PANS [85], and PITM [86] all at approximately 1M cell counts using codes with second-order numerics. AMS and these other hybrid methods are clearly more cost-effective than WRLES (Table II). Note that at worst only a factor of 2 more equations are solved in AMS, due to the inclusion of RANS, than in the dynamic Smagorinsky WRLES of [2], while obviating the dynamic procedure, making AMS on even the finest grid at least an order of magnitude less costly. But, it is appropriate to ask how the cost of AMS would compare to a WMLES of the same problem. Assuming that the wall-model is valid in the separation, recirculation, and reattachment regions and that it does not require a separate highly-refined grid [12], we can estimate the cost for a WMLES based on resolution requirements for boundary layers [10,87] in which the streamwise, spanwise and wall normal grid sizes should be no larger than 8%, 2% and 5% of the boundary layer thickness, respectively. The most restrictive point is the top of the hill, where δ 99 ≈ 0.08h, and for a structured grid, this will set the number of grid cells in the wall-normal and spanwise directions throughout the domain. Streamwise spacing, however, can adjust to the local boundary layer thickness, which is so small for only about 0.5h in the streamwise direction. Assuming that δ 99 ≈ 1h in the remainder of the domain yields a total required grid of approximately 12M cells, which is only slightly smaller than that required in WRLES [2]. The primary driver of this large cell count is the spanwise resolution requirement at the top of the hill which, due to the structured grid, is used throughout the domain. If one could use a grid in which the spanwise and wall-normal resolution varies in the streamwise direction, then the WMLES grid could be reduced to about 1.6M cells, which is still larger than the fine AMS grid in Table II. However, this analysis greatly understates the cost reductions of HRL relative to WMLES for the types of external flows in which HRL are likely to be applied (e.g. an airfoil or even aircraft). When applied to such flows, a HRL can represent the turbulent boundary layer over most of the body with RANS, requiring very coarse LES resolution only in critical regions such as near separation. In a WMLES, however, LES resolution is needed for all turbulent boundary layers, which will be orders of magnitude more expensive.
Caveat the slightly excessive normal stress components just after the the separation, the reduced grid sizes necessary with AMS are encouraging. In addition to revealing an avenue for formulation improvement, this test case has shown the ability of AMS models to accurately predict complex flow features at reasonable computational costs. The tolerance of under-resolution (Fig. 11) indicates that perhaps even more coarse resolutions could be used by relaxing the r M requirements. We examine this possibility next.

C. Sensitivity to resolution parameter
The assumption that C r = 1 in the evaluation of the available resolution in Eq. 34 has not been strongly justified to this point. This coefficient essentially indicates how many grid lengths are necessary to resolve some convecting turbulent structure. Clearly this coefficient should change depending on the characteristics of the spatial discretization, including the order of accuracy and dispersion relation [88]. For increasing order of accuracy and decreasing dispersion errors, small structures approaching the grid scale are effectively better resolved, and hence a lower C r should be possible. Determining exact values of C r for different order methods is beyond the scope this paper. Motivated by the observed tolerance in the under-resolved regions of periodic hill problem, we now examine the effects of reducing C r with the second-order finite volume numerics of CDP. First, consider the coarse resolution channel with C r = 0.5, 0.2, and 0.1 (Fig. 13). Lowering C r from unity to 0.5 results in a drop in the center-channel β from just above 0.5 to 0.2 with resolved turbulence moving towards the wall from y + ≈ 500 to 100. Only a slight change in the log-layer slope and middle channel mean velocity are observed. Both measures of the total TKE (modeled and resolved plus model times β) move toward the DNS value with a small excess TKE at the center of the channel. However, reducing C r below 0.5 results in both distortions to the mean velocity and large errors in TKE. Errors introduced into the TKE are so pronounced that the outer peak is predicted to be larger than the inner with C r = 0.1. Thus, it seems the small errors incurred with C r = 0.5 may be a tolerable exchange for the increase in resolved turbulence on a given grid, but moving below C r = 0.5 does not appear to be acceptable.  Table II using different Cr values, along with data from the experiments of [82] and WRLES of [2]. Shown are the streamwise velocity variance (top) vertical velocity variance (middle) and Reynolds shear stress (bottom). Note amplification by a factor of 10 for visibility.
Next, we examine the effects of reduced C r in the more complex case of the periodic hill. Figure 14 shows how reducing C r has the same basic effect as increasing the resolution with C r = 1. For C r = 0.5, the reattachment location is well predicted though profiles at x = 6 and 7 still deviate slightly from the experimental data. Contrary to the channel, C r = 0.2 actually gives the best prediction of the mean streamwise velocity with values nearly identical to the experiment and the fine resolution hybrid case with C r = 1. However, the Reynolds stress components (Fig. 15) are over-predicted as observed in the channel. This is most pronounced in the prediction of the streamwise variance where the hybrid profiles are significantly higher than the experimental value at x = 2, 3, and 4 in the recirculation region. The good mean velocity results for C r = 0.2 appear to be a result of error canceling. Again consistent with the channel results, lowering C r to 0.5 appears to introduce only small errors in exchange for the most resolved turbulence possible for a given resolution. So, it would seem that C r = 1 is too conservative and the value of 0.5 may be more appropriate for the numerical schemes used here.

D. Alternative energy transfer models
To this point, we have used the tensor-diffusivity M43 model for the energy transfer portion of the modelsplit formulation. Though the model has been shown to perform well, many existing RANS-based codes do not currently support a tensor diffusivity, making implementation and adoption of the AMS method more difficult. It is therefore useful to investigate the performance of AMS using more easily implemented models for τ e ij . The three additional models discussed here are formulated in terms of expected values so that they do not contribute to mean stress portion of AMS. The first such model (labeled α-based) was introduced in Sec. III B. Based on expanding the total sugbrid term in light of arguments leading to eddy viscosity formulations, it was suggested that the fluctuating gradient contribution to the model should scale with α as α(1 − α)ν tot S > ij (see Eq. 23), where ν tot is the RANS eddy viscosity. This α dependence has the correct RANS and DNS limits, since it vanishes where there are no resolved fluctuations (α = 1) and where all turbulence is resolved (α = 0). The second simple energy transfer model (iM43 for isotropic M43) is an isotropic version of the M43 model, with the length scale defined as the cell-diagonal, i.e. ν e = C(M)ε 1/3 δ 4/3 diag where δ diag = (M ij M ji ) 1/2 . Finally, the third model (iSGET for implicit subgrid energy transfer) is not formulated in terms of an explicit eddy viscosity. It is well-known that first-order upwinding of the convection term results in numerical dissipation. Motivated by the streamwise-upwind Petrov-Galerkin (SUPG) method [67], we introduce upwinding depending on a cell-Reynolds number, but in this case the Reynolds number is based on the subgrid turbulence intensity, as determined from βk tot (see Appendix C for details).
Results for these alternative models are presented in Fig. 16 for the fully-developed channel. Here, the fine resolution defined in Table I is used because the energy transfer models are most important with higher levels of resolved turbulence. The previously presented M43 results are also included for comparison. While the anisotropic M43 model performs the best, both the α-based and iSGET models also perform well. The α-based model produces only a slight deviation from the log law centered around y + ≈ 300, while iSGET slightly distorts the mean velocity profile above y + ≈ 2000. Interestingly, the α-based model results in the least resolved turbulence and iSGET resolves significantly more turbulence above y + ≈ 2000. The likely reason for the latter is that we have not counted the numerical dissipation as part of the production of subgrid energy in the definition of r M in Eq. 34. This is an inconsistency in the implementation of iSGET into AMS, which would need to be corrected by calculating the mean rate of numerical dissipation introduced by the upwinding. Many CFD codes, including CDP used here, are not instrumented for this. The relatively good results obtained using this incomplete integration of iSGET into AMS suggest that implementing the numerical dissipation diagnostics required to complete the AMS integration may be worthwhile. Perhaps correcting this would also improve the slight mean velocity distortion. The good results with iSGET also confirm that the AMS formulation can be used even with implicit LES models for the energy transfer. However, the iM43 model performs poorly, with a significant shift in the log-layer due to over-prediction of the wall shear stress. The reason for this is clear from the turbulent kinetic energy. While the other models produce similar k tot profiles, the iM43 results in k tot values that are far too high for y + 1000. The higher turbulence levels enhance momentum transfer towards the wall, increasing the wall shear stress. Thus, it appears the M43 model does require the anisotropic form to perform well.
From the perspective of both flexibility in model selection and implementation of AMS, these results are encouraging. The α-based model piggy-backs on the mean stress portion of the split-model while iSGET simply exploits commonly used numerical approximations. Neither model requires a tensor-diffusivity. As mentioned above, both these models are formulated in terms of expected values of eddy viscosity or cell Reynolds number so they do not contribute to the mean stress. In AMS, using energy transfer models with fluctuating eddy viscosities, e.g. Smagorinsky, Vreman, WALE, etc., or a fluctuating cell Re in implicit methods, would require that the contribution of the energy transfer term to the mean stress be computed, and the mean stress model adjusted accordingly.

VI. CONCLUSIONS
The deficiencies of RANS models in representing complex flow features like separation, and the expense of using LES to represent, for example, broad expanses of a boundary layer, make a hybridization of these two model paradigms (HRL) compelling. Unfortunately, active development of such HRL for more than 20 years has not resulted in widely applicable robust and accurate HRL models. As detailed in Sec. II, there are four interrelated deficiencies common to most HRL formulations that appear responsible for this state of affairs: first is the use of a single eddy viscosity to represent both the mean subgrid stress and the transfer of energy to the small scales; second is the blending of RANS and LES models to hybridize them; third is the application of RANS-based transport models to fluctuating quantities; and last is reliance on passive self-generation of resolved turbulence fluctuations when such fluctuations need to be introduced. These deficiencies lead to common HRL failure modes, such as log-layer miss-match and model stress depletion.
The Active Model Split (AMS) hybrid formulation described in Section IV was developed specifically to eliminate these deficiencies. It does so by "splitting" the model into separate models to represent the mean subgrid stress and the transfer of energy to the unresolved scales, and by actively stirring where necessary to produce resolved fluctuations. The model splitting eliminates the overloading of a single eddy viscosity with two roles, since separate models are used for mean stress and energy transfer. It also eliminates blending because the model is essentially LES everywhere, with a consistent model for the mean subgrid stress. This allows LES to be used with much coarser resolution than would otherwise be required, since it is not necessary for the subgrid contribution to the mean stress to be negligible, which is important since in HRL there will generally be regions of very coarse LES. Further, model splitting allows RANS models to be used as designed, only for the mean stress and acting only on average quantities. Finally, active stirring explicitly eliminates the need for fluctuations to develop due to natural instabilities of the mean.
The model split formulation does not just address challenges in hybrid RANS LES modeling, it also addresses similar challenges in LES. In particular, by eliminating the need for the resolved fluctuations to carry the majority of the Reynolds stress, model splitting reduces the resolution requirements for LES of inhomogeneous turbulent flows. Further, by reverting to RANS if the LES resolution is not sufficient to represent near-wall fluctuations, the model splitting formulation provides a natural wall representation for LES, as demonstrated in the channel results reported in Section V A. One can argue then that model splitting like that described here is how LES should generally be performed.
The model tests in channel flow and the periodic hill indicate that the AMS formulation has indeed addressed the primary shortcomings of HRL models. In the channel flow, it produces generally consistent results for horizontal resolution ranging over a factor of at least 2.2, showing both resolution independence and a lack of log-layer miss-match. Further, the AMS solution remains consistent as the model transitions from RANS to LES in either time or space. In the more complex periodic hill case, in which RANS performs poorly due to flow separation, the AMS formulation produces mean velocity in good agreement with both experiments and wall-resolved LES, with a grid that is up to 60 times smaller than the LES. These test results suggest that the AMS approach is a solution to the common hybrid modeling issues identified above and successfully eliminates hybridization artifacts such as modeled-stress depletion. Evaluations of the AMS formulation in more complex flow scenarios are clearly warranted. These should include common aerodynamic test cases with varying degrees of smooth wall separation and reattachment, and flows that are unsteady in the mean.
While the model test results are quite good, the AMS solutions are not perfect. Further, there are several details that must be specified to complete the definition of a particular AMS implementation. For the model and tests presented here, these details were described in Section IV and are summarized in Appendix A and B, but further investigation and refinements would clearly be useful. The most significant of these are listed here.
1. Forcing Formulation: The forcing formulation described in Section IV D is an ad hoc place holder and has shortcomings. A formulation that is more realistic in structure, does not need to be clipped, has a more controllable energy injection rate and can extract resolved energy when needed would be a great improvement. Such improvements could allow more rapid transitions from RANS to LES.

Resolution Inhomogeneity:
As is well known, when LES resolution is inhomogeneous, the filter operator that defines the resolved scales and the spatial derivative operator do not commute. When there is a mean flow through the inhomogeneous grid, the commutator represents the transfer of energy between resolved and unresolved fluctuations, as needed as the resolution changes. This generally goes unmodeled, and its effects are particularly acute with HRL because hybrid simulations are commonly done with strongly inhomogeous grids, indeed that is the objective. Thus, the commutator needs to be modeled, and the forcing formulating may be useful in this regard.
or the averaging time scale have been carefully investigated, so refinements are likely to be appropriate. For example, in the presence of a mean velocity, the convective turbulent time scale k 3/2 /(ε| u |) is also relevant. The appropriate averaging time scale is of some importance because, as discussed in Section IV B, the interactions of the largest turbulent scales with the subgrid are poorly described by gradient transport models, so these largest scale fluctuations are best included in the pseudomean. Finally, the interaction of the pseudo-mean time averaging with mean unsteadiness needs to be investigated.
These opportunities for further refinement are essentially to improve the LES component of the AMS formulation. Despite the need for some refinements, by addressing the primary shortcomings of most HRL formulations, AMS provides a framework for robust, predictive, and cost-effective simulation complex turbulent flows.  where the turbulent Reynolds number is Re T = k 2 /(νε). Finally, the following coefficients close Chien's model C ε1 = 1.35, C ε2 = 1.8.
The length scale used in the f -equation for v 2 -f is given by L = max min k 3/2 ε , k 3/2 √ 6C µ v 2 S , C η ν 3/4 ε 1/4 (A2) and an additional time scale modification is made here with which was observed to not affect basic RANS behavior but to be necessary for use with AMS. Without this modification, v 2 becomes excessive in hybrid simulations. The general applicability of this modification deserves further study. The model is complete with the following coefficients    This modification is motivated by the fact that r M is a length scale ratio. When r M > 1 while resolved fluctuations are not zero, the simulation is locally under-resolved by an amount indicated by r M . Thus, the length scales in the M43 eddy viscosity are modified by this length scale ratio. However, it is essentially ad-hoc. A more principled formulation is needed.