Higher-order realizable algebraic Reynolds stress modeling based on the square root tensor

TL;DR

This paper develops a higher-order algebraic Reynolds stress model based on the square root tensor, ensuring realizability and improving predictions for complex three-dimensional turbulent flows.

Contribution

It introduces a quartic-nonlinear eddy-viscosity model that maintains realizability even with higher-order nonlinearities, extending applicability to complex turbulent flows.

Findings

Model performs well in turbulent channel flow

Accurately predicts mean swirl flow in rotating pipe

Ensures realizability with higher-order nonlinear terms

Abstract

In this study, realizable algebraic Reynolds stress modeling based on the square root tensor [Phys. Rev. E \textbf{92}, 053010 (2015)] is further developed for extending its applicability to more complex flows. In conventional methods, it was difficult to construct an algebraic Reynolds stress model satisfying the realizability conditions when the model involves higher-order nonlinear terms on the mean velocity gradient. Such higher-order nonlinear terms are required to predict turbulent flows with three-dimensional mean velocity. The present modeling based on the square root tensor enables us to make the model always satisfy the realizability conditions, even when it involves higher-order nonlinearity. To construct a realizable algebraic Reynolds stress model applicable to turbulent flows with three-dimensional mean velocity, a quartic-nonlinear eddy-viscosity model is proposed. The performance of the model is numerically verified in a turbulent channel flow, a homogeneous turbulent shear flow, and an axially rotating turbulent pipe flow. The present model gives a good result in each turbulent flow. Note that the mean swirl flow in an axially rotating turbulent pipe flow is reproduced because the present model involves cubic nonlinearity. Such a higher-order realizable algebraic Reynolds stress model, involving quartic nonlinearity on the mean velocity, is expected to be useful in numerically stable predictions of turbulent flows with three-dimensional mean velocity.

Tables

Table 1. Table 1: Model constants.

Model	$C_1$	$C_2$	$C_3$	$C_{\nu }$	$C_{\epsilon 1}$	$C_{\epsilon 2}$	${\sigma }_K$	${\sigma }_{\epsilon }$	$a_1$	$a_2$	$a_3$	$a_{\epsilon 1}$	$a_{\epsilon 2}$	$a_{\epsilon 3}$
present	0.13	0.021	0.018	0.09	1.5	1.9	1.4	1.4	14	5	200	3.1	0.3	6.5
linear Abe et al. (1994)	-	-	-	0.09	1.5	1.9	1.4	1.4	14	5	200	3.1	0.3	6.5

Equations

$$R_{\alpha \alpha }\ge 0,$$

$$R_{\alpha \alpha }{R}_{\beta \beta }\ge R_{\alpha \beta }^2,$$

$$\det \left[R_{ij}\right]\ge 0,$$

$${\frac{\mathrm{d}F}{\mathrm{d}t}|}_{F=0}=0,{\frac{{\mathrm{d}}^{2}F}{\mathrm{d}{t}^2}|}_{F=0}>0,$$

$${\frac{\mathrm{d}F}{\mathrm{d}t}|}_{F=0}>0.$$

$$R_{ij}=\frac{2}{3}K{\delta }_{ij}-{\nu }_{\mathrm{T}}{S}_{ij}.$$

$$S_{ij}=\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}.$$

$${\nu }_{\mathrm{T}}=C_{\nu }\frac{K^2}{\epsilon }$$

$$R_{xx}=\frac{2}{3}K-{\nu }_{\mathrm{T}}\frac{\partial U_x}{\partial x}=\frac{2}{3}K\left(1-\frac{3}{2}{C}_{\nu }\frac{K}{\epsilon }\frac{\partial U_x}{\partial x}\right).$$

$$R_{ij}$$

$$+{\zeta }_{SS}\left(S_{ia}{S}_{aj}-\frac{1}{3}{S}^2{\delta }_{ij}\right)+{\zeta }_{SW}\left(S_{ia}{W}_{aj}+S_{ja}{W}_{ai}\right)+{\zeta }_{WW}\left(W_{ia}{W}_{aj}+\frac{1}{3}{W}^2{\delta }_{ij}\right)$$

$$+{\zeta }_{SSW}\left(S_{ia}{S}_{ab}{W}_{bj}+S_{ja}{S}_{ab}{W}_{bi}\right)+{\zeta }_{SWW}\left(S_{ia}{W}_{ab}{W}_{bj}+S_{ja}{W}_{ab}{W}_{bi}-\frac{2}{3}{S}_{ab}{W}_{bc}{W}_{ca}{\delta }_{ij}\right)$$

$$+\cdots ,$$

$$W_{ij}=\frac{\partial U_i}{\partial x_j}-\frac{\partial U_j}{\partial x_i}.$$

$$R_{ij}={\sqrt{R}}_{ia}{\sqrt{R}}_{ja}.$$

$${\sqrt{R}}_{ij}={\gamma }_0{\delta }_{ij}-{\gamma }_S{S}_{ij}.$$

$$R_{ij}={\gamma }_0^2{\delta }_{ij}-2{\gamma }_0{\gamma }_S{S}_{ij}+{\gamma }_S^2{S}_{ia}{S}_{ja}.$$

$${\sqrt{R}}_{ij}={\gamma }_0{\delta }_{ij}+\sum _{\mathrm{m}}{\gamma }_{\mathrm{m}}{T}_{ij}^{\mathrm{m}},$$

$$2K=R_{ii}={\sqrt{R}}_{ij}{\sqrt{R}}_{ij}=3{\gamma }_0^2+2\sum _{\mathrm{m}}{\gamma }_0{\gamma }_{\mathrm{m}}{T}_{ii}^{\mathrm{m}}+\sum _{\mathrm{m},\mathrm{n}}{\gamma }_{\mathrm{m}}{\gamma }_{\mathrm{n}}{T}_{ij}^{\mathrm{m}}{T}_{ij}^{\mathrm{n}}.$$

$$2K=3{\gamma }_0^2\left(1+\frac{2}{3}\sum _{\mathrm{m}}{f}_{\mathrm{m}}{T}_{ii}^{\mathrm{m}}+\frac{1}{3}\sum _{\mathrm{m},\mathrm{n}}{f}_{\mathrm{m}}{f}_{\mathrm{n}}{T}_{ij}^{\mathrm{m}}{T}_{ij}^{\mathrm{n}}\right).$$

$${\gamma }_0=\sqrt{\frac{2K}{3D}},D=1+\frac{2}{3}\sum _{\mathrm{m}}{f}_{\mathrm{m}}{T}_{ii}^{\mathrm{m}}+\frac{1}{3}\sum _{\mathrm{m},\mathrm{n}}{f}_{\mathrm{m}}{f}_{\mathrm{n}}{T}_{ij}^{\mathrm{m}}{T}_{ij}^{\mathrm{n}}.$$

$${\sqrt{R}}_{ij}={\gamma }_0{\delta }_{ij}-{\gamma }_S{S}_{ij}-{\gamma }_W{W}_{ij}-{\gamma }_N{N}_{ij},$$

$$R_{ij}=\frac{2K}{3D}$$

$$+f_N{f}_S\left(N_{ia}{S}_{aj}+N_{ja}{S}_{ai}\right)-f_N{f}_W\left(N_{ia}{W}_{aj}+N_{ja}{W}_{ai}\right)+f_N^2{N}_{ia}{N}_{ja}],$$

$$D=1+\frac{1}{3}{f}_S^2{S}^2+\frac{1}{3}{f}_W^2{W}^2-\frac{2}{3}{f}_N{N}_{ii}+\frac{2}{3}{f}_S{f}_N{S}_{ij}{N}_{ij}+\frac{2}{3}{f}_W{f}_N{W}_{ij}{N}_{ij}+\frac{1}{3}{f}_N^2{N}^2,$$

$${\nu }_{\mathrm{T}}=\frac{4K{f}_S}{3D}=\frac{4C_S}{3D}\frac{K^2}{\epsilon }=\frac{4C_S}{3}\frac{1}{1+C_S^2{\hat{S}}^2/3+C_W^2{\hat{W}}^2/3+\cdots }\frac{K^2}{\epsilon },$$

$${\sqrt{R}}_{ij}={\gamma }_0{\delta }_{ij}-{\gamma }_S{S}_{ij}-{\gamma }_W{W}_{ij}-{\gamma }_C{C}_{ij},$$

$$C_{ij}=S_{ia}{W}_{aj}+S_{ja}{W}_{ai}.$$

$$R_{ij}=\frac{2K}{3D}$$

$$+f_C{f}_S\left(C_{ia}{S}_{aj}+C_{ja}{S}_{ai}\right)-f_C{f}_W\left(C_{ia}{W}_{aj}+C_{ja}{W}_{ai}\right)+f_C^2{C}_{ia}{C}_{aj}],$$

$$D=1+\frac{1}{3}{f}_S^2{S}^2+\frac{1}{3}{f}_W^2{W}^2+\frac{1}{3}{f}_C^2{C}^2,$$

$$R_{xx}$$

$$R_{yy}$$

$$R_{zz}$$

$$R_{xy}$$

$$D$$

Full text

I Introduction

Turbulent flows are ubiquitous in real-world fluid flows. Despite remarkable advancements in modern computational technology, most of real-world turbulent flows of high-Reynolds number are still out of reach of direct numerical simulation (DNS). In this sense, the Reynolds averaged Navier-Stokes (RANS) modeling, which treats the statistically averaged quantities of turbulent flows, is a useful tool in predicting high-Reynolds-number turbulent flows.

In a variety of RANS modeling, Reynolds stress transport modeling (RSM, also referred to as second-order closure modeling) is expected to have a potential to give a good prediction for several complex turbulent flows. This is because RSM explicitly calculates the production, dissipation, and diffusion of the Reynolds stress. However, RSM has a problem of realizability; namely, it sometimes gives unphysical predictions such as negative normal stresses Hanjalić and Launder (2011); Lumley (1978). Schumann Schumann (1977) showed three inequality conditions, referred to as the realizability conditions for the Reynolds stress $R_{ij}=\langle u_{i}^{\prime}u_{j}^{\prime}\rangle$, where $u_{i}^{\prime}$ denotes the velocity fluctuation. They are written as follows:

[TABLE]

where $\alpha,\beta=1,2,3$, $\mathrm{det}[R_{ij}]=\epsilon_{ij\ell}\epsilon_{abc}R_{ia}R_{jb}R_{\ell c}/6$, and $\epsilon_{ij\ell}$ is the alternating tensor. Note that summations are not taken for Greek indices. The above conditions, Eqs. (1a)–(1c), are mathematically rigorous and provide a physical importance. For example, the first condition, Eq. (1a), prevents an excessively anisotropic condition caused by the negative normal stresses $R_{\alpha\alpha}<0$ and the second condition, Eq. (1b), gives the upper bound of the Reynolds shear stress. These lead to the numerical stability of the model. Hence, the realizability conditions given by Eqs. (1a)–(1c) form a fundamental guideline for the development of physically consistent and numerically stable RANS modeling. However, many simple RANS models, even for the standard eddy-viscosity model, do not necessarily satisfy the realizability conditions.

For RSM, realizability is often discussed in terms of the time development of the Lumley’s flatness $F\{=\mathrm{det}[R_{ij}/(R_{\ell\ell}/3)]\}$ around $F=0$. Here, $F=0$ indicates that turbulence field involves only two-component velocity fluctuation, which is referred to as two-component limit (TCL) condition. In this context, two types of realizability conditions are argued Pope (2000); one is the strong realizability,

[TABLE]

and the other is the weak realizability,

[TABLE]

Durbin and Speziale Durbin and Speziale (1994) gave a stochastic analysis on the weak realizability of the conventional RSM in terms of Langevin equations. Craft and Launder Craft and Launder (1996); Craft and Laudner (2001) proposed highly sophisticated RSM by considering TCL condition in the vicinity of the solid wall. On further development of the realizability conditions for RSM, Girimaji Girimaji (2004) discussed and proposed a new constraint on the pressure–strain correlation, and its utility was demonstrated by Mishra and Girimaji Mishra and Girimaji (2014).

Although RSM has an advantage that it involves production, dissipation, and diffusion process of the Reynolds stress, its numerical cost is still large because RSM, by nature, demands to solve six individual transport equations for the Reynolds stress. To reduce the numerical cost while retaining the advantages of RSM, the algebraic Reynolds stress modeling (ARSM, or the so-called eddy-viscosity type modeling) is a useful method. ARSM gives an algebraic expression for the Reynolds stress instead of numerically integrating its transport equation. In case of explicit algebraic Reynolds stress modeling (EARSM) which is constructed based on RSM, the realizability of the based RSM is tightly connected to the property of the proposed model. However, when we give an algebraic expression independent of RSM, we should concern Eqs. (1a)–(1c) instead of Eqs. (2) or (3). In this paper, we discuss the realizability conditions for ARSM independent of a specific RSM. Note that the third condition (1c) is equivalent to the positive semidefiniteness of the Lumley’s flatness, $F\geq 0$. Therefore, ARSM satisfying Eqs. (1a)–(1c) ensures the weak realizability at least.

To date, many studies on ARSM have been proposed (see, e.g., Ref. Pope (1975); Yoshizawa (1984); Taulbee (1992); Gatski and Speziale (1993); Shih et al. (1993, 1995); Girimaji (1996); Abe et al. (1997); Wallin and Johanson (2000); Hamba (2001); Abe et al. (2003); Ariki (2019)). However, only a few studies on ARSM consider the realizability conditions. A representative study of realizable ARSM was conducted by Shih et al. Shih et al. (1993, 1995). They constructed ARSM incorporated with quadratic nonlinearity on the mean velocity gradient. Note that the model proposed by Shih et al. Shih et al. (1993, 1995) is a quadratic-nonlinear model, and hence, is adequate for turbulent flows with two-dimensional mean velocity Pope (1975). Such flows denote the flows in which the mean velocity field is settled in two dimensions, such as turbulent flows in a plane channel or a straight pipe. Hereafter, we simply refer to turbulent flows with two- or three-dimensional mean velocity as two- or three-dimensional flows, respectively. Some other ARSM which consider the realizability Abe et al. (1997); Hamba (2001); Abe et al. (2003) are also quadratic-nonlinear models. However, it was pointed out that cubic nonlinearity is essential for predicting the mean swirl or tangential velocity in an axially rotating turbulent pipe flow, which is an example of the three-dimensional flow Speziale et al. (2000). The realizable ARSM involving cubic or higher-order nonlinearity is required to make numerically stable predictions of three-dimensional flows of high-Reynolds number.

Generally, the conventional realizable ARSM can hardly be extended to the models incorporated with higher-order nonlinearity on the mean velocity gradient. Recently, Ariki Ariki (2015) proposed an alternative approach to realizable ARSM by introducing the square root tensor of the Reynolds stress. In other words, all of eigenvalues of the Reynolds stress constructed by this modeling must be positive semidefinite. Although this approach is a sufficient condition that the model is realizable, it provides a rigorous realizability in a systematic procedure. In the present study, we further develop the square root modeling method, aiming at its application to more complex turbulent flows. Considering the application to three-dimensional flows, a model involving quartic nonlinearity on the mean velocity gradient is proposed and its performance is numerically verified in a turbulent channel flow, a homogeneous turbulent shear flow as basic two-dimensional flow, and an axially rotating turbulent pipe flow as an example of the three-dimensional flow.

The rest of this paper is organized as follows. In Sec. II, we review the conventional realizable turbulence modeling. Here, we also give an introduction and further development of the turbulence modeling based on the square root tensor Ariki (2015). In Sec. III, we construct a quartic-nonlinear model of the Reynolds stress and examine the basic property of the model for the unidirectional shear flow. In Sec. IV, numerical verifications of the proposed model are given. A summary is given and conclusions are discussed in Sec V.

II Modeling the Reynolds stress based on the square root tensor

II.1 Realizability and conventional turbulence modeling

The simplest and frequently used model for the Reynolds stress is the linear eddy-viscosity model:

[TABLE]

where $K(=\langle u_{i}^{\prime}u_{i}^{\prime}\rangle/2))$ denotes the turbulent energy, $\nu_{\mathrm{T}}$ denotes the eddy viscosity, and $S_{ij}$ denotes the mean strain tensor, which is written in an inertial frame of the Cartesian coordinates as

[TABLE]

Here, $U_{i}$ denotes the mean velocity. In case of the standard $K$-$\varepsilon$ model, $\nu_{\mathrm{T}}$ is written as

[TABLE]

where $\varepsilon=\nu\langle(\partial u_{i}^{\prime}/\partial x_{j})^{2}\rangle$ denotes the dissipation rate of the turbulent energy and $C_{\nu}$ is the model constant often optimized as $C_{\nu}=0.09$ in a turbulent channel flow Yoshizawa (1998); Pope (2000). In this model, $R_{xx}$ is written as

[TABLE]

As seen from Eq. (7), $R_{xx}$ can take a negative value when $(3C_{\nu}/2)(K/\varepsilon)\partial U_{x}/\partial x>1$, which violates Eq. (1a). Such an unphysical behavior of the Reynolds stress can cause numerical instability. An easy approach to rectify this problem is to extend the constant $C_{\nu}$ to a functional of the mean strain rate, e.g., $C_{\nu}\to C_{\nu}/(1+C_{\nu}^{\prime}K^{2}S^{2}/\varepsilon^{2})$, where $S^{2}=S_{ij}S_{ij}$. Such a functionalization of the coefficient is a primitive method for the realizable modeling. It is also known that the eddy-viscosity model with $C_{\nu}=0.09$ often overestimates the turbulent energy production rate in homogeneous turbulent shear flows Yoshizawa (1998). In this sense, the eddy-viscosity coefficient $\nu_{\mathrm{T}}$ should be generally modeled as a functional form in terms of $S^{2}$ or other scalar variables in addition to $K$ and $\varepsilon$. In this course, a systematic way of functionalization of the eddy-viscosity coefficient is a matter of modeling.

For the generalization of ARSM, the Reynolds stress is often expanded by the mean velocity gradient Pope (1975), which leads to the following nonlinear eddy-viscosity model:

[TABLE]

where $\zeta$’s are dimensional coefficients, $W^{2}=W_{ij}W_{ij}$, and $W_{ij}$ denotes the mean absolute vorticity tensor, which is written in an inertial frame of the Cartesian coordinates as

[TABLE]

For a representative of realizable ARSM, Shih et al. Shih et al. (1993, 1995) constructed a simple expression for the Reynolds stress in which $\nu_{\mathrm{T}}$ and $\zeta_{SW}$ are retained but other $\zeta$’s are neglected. In their modeling, the functional form of $\nu_{\mathrm{T}}$ and $\zeta_{SW}$ were determined to satisfy the realizability conditions for simple shear or rotating turbulent flows. Hamba Hamba (2001) discussed a realizable ARSM through the statistical theory for inhomogeneous turbulence referred to as the two-scale direct-interaction approximation (TSDIA) Yoshizawa (1984). Abe et al. Abe et al. (1997) modified the ARSM proposed by Gatski and Speziale Gatski and Speziale (1993) to satisfy the realizability conditions for simple flows and Abe et al. Abe et al. (2003) further developed its model with careful consideration on the near wall TCL condition. Note that all of them are quadratic-nonlinear eddy-viscosity models. Although these models can capture the anisotropic property of the Reynolds stress in turbulent shear flows owing to the quadratic-nonlinear terms, it is almost impossible to extend these models to higher-order nonlinear ones because the freedom for the dimensional coefficients $\zeta$’s is too large to analytically determine their form. However, Speziale et al. Speziale et al. (2000) pointed out that cubic nonlinearity is essential for predicting the mean swirl or tangential velocity in an axially rotating turbulent pipe flow, which is an example of the three-dimensional flow. Furthermore, the Reynolds stress should be affected not only by the mean velocity gradient but also the turbulent helicity Yokoi and Yoshizawa (1993); Yokoi and Brandenburg (2016); Inagaki et al. (2017) or the time history effect of the mean strain Speziale (1987); Hamilington and Dahm (2008); Ariki (2019). The conventional realizable modeling method is hardly extended to such a more general turbulent flow.

II.2 Modeling based on the square root tensor

Ariki Ariki (2015) presented another possibility of realizable ARSM by introducing the square root tensor of the Reynolds stress $\sqrt{R}_{ij}$, which is defined as

[TABLE]

Note that $\sqrt{R}_{ij}$ does not signify the square root of each component of the Reynolds stress $R_{ij}$, i.e., $\sqrt{R}_{ij}\neq\sqrt{R_{ij}}$. A physical interpretation of the square root tensor is discussed in Appendix A. In this modeling, we choose the square root tensor $\sqrt{R}_{ij}$ as a target of modeling instead of $R_{ij}$ itself. It should be emphasized that the ARSM proposed by this methodology always satisfies the realizability conditions given by Eqs. (1a)–(1c). In other words, all of eigenvalues of the Reynolds stress are always positive semidefinite. A primitive model for $\sqrt{R}_{ij}$ is an eddy-viscosity type model:

[TABLE]

The model given by Eqs. (10) and (11) yields the following model of the Reynolds stress:

[TABLE]

Thus, even a linear model for $\sqrt{R}_{ij}$ corresponds to a quadratic-nonlinear expression for $R_{ij}$ on the mean velocity gradient The linear eddy-viscosity model can be derived when we add $-\gamma_{SS}S_{ia}S_{aj}$ to Eq. (11) and choose $\gamma_{SS}$ to vanish $S_{ia}S_{aj}$ term in $R_{ij}$. Let us write a general expression for the model of $\sqrt{R}_{ij}$ as

[TABLE]

where $\{T_{ij}^{\mathrm{m}}\}$ indicate arbitrary algebraic model terms for $\sqrt{R}_{ij}$, such as $S_{ij}$, $W_{ij}$, or higher-order tensors composed by them. In a previous study by Ariki Ariki (2015), $\{T_{ij}^{\mathrm{m}}\}$ were restricted to symmetric real tensors, but this was not a necessary condition. In this study, $\{T_{ij}^{\mathrm{m}}\}$ are allowed to include antisymmetric real tensors.

II.3 Further development of square root modeling

As a further development of square root modeling Ariki (2015), we propose a systematical approach to determine the functional form of $\gamma_{\mathrm{m}}$’s following Hamba Hamba (2001). As the trace of the Reynolds stress must correspond to the turbulent energy $K=\langle u_{i}^{\prime}u_{i}^{\prime}\rangle/2$ as $R_{ii}=2K$, the model given by Eqs. (10) and (13) should satisfy

[TABLE]

Here, we introduce a fraction $f_{\mathrm{m}}(=\gamma_{\mathrm{m}}/\gamma_{0})$ by assuming $\gamma_{0}\neq 0$. Then, Eq. (14) reads

[TABLE]

Thus, $\gamma_{0}$ can be defined as

[TABLE]

Using this definition of $\gamma_{0}$, the model given by Eqs. (10) and (13) always satisfies Eq. (14). Note that $D>0$ always holds because it is expressed by a square value, namely $D=(\delta_{ij}+\sum_{\mathrm{m}}f_{\mathrm{m}}T^{\mathrm{m}}_{ij})^{2}/3$. Hence, the present modeling gives regular dimensional coefficients corresponding to $\nu_{\mathrm{T}}$ and $\zeta$’s, in contrast to the singular behavior of coefficients in EARSM Gatski and Speziale (1993); Girimaji (1996); Abe et al. (1997); Wallin and Johanson (2000). Although the present way to determine $\gamma_{\mathrm{m}}$’s is not a unique method for making the model satisfy Eq. (14), it suggests a physically interesting feature of the model. To see this point clearly, we consider the following model:

[TABLE]

where $N_{ij}$ represents an additional model term. To construct a higher-order nonlinear eddy-viscosity model, we substitute $S_{ia}S_{ja}$ or $S_{ia}W_{aj}+S_{ja}W_{ai}$ into $N_{ij}$. To construct a more general model for the Reynolds stress, we substitute the helicity term Yokoi and Yoshizawa (1993); Yokoi and Brandenburg (2016); Inagaki et al. (2017) or the time history effect of the mean strain tensor Speziale (1987); Hamilington and Dahm (2008); Ariki (2019) into $N_{ij}$. It should be emphasized that the square root modeling makes the Reynolds stress always realizable, even when the complex model terms stated above are adopted to $N_{ij}$. The Reynolds stress given by Eqs. (10), (16), and (17) is written as

[TABLE]

where

[TABLE]

and $N^{2}=N_{ij}N_{ij}$. To observe the physical property of the model more clearly, we adopt a simple coefficient form to $f_{S}$ and $f_{W}$; namely, we give $f_{S}=C_{S}K/\varepsilon$ and $f_{W}=C_{W}K/\varepsilon$, where $C_{S}$ and $C_{W}$ are constant. In such a case, the dimensional coefficient corresponds to the eddy viscosity $\nu_{\mathrm{T}}$ [see Eq. (8)] in the model given by Eqs. (18) and (19) is written as follows:

[TABLE]

where $\hat{S}^{2}=K^{2}S^{2}/\varepsilon^{2}$ and $\hat{W}^{2}=K^{2}W^{2}/\varepsilon^{2}$. Namely, the eddy viscosity is not simply expressed by Eq. (6) but is incorporated with the effect of the mean strain rate and vorticity, $\hat{S}^{2}$ and $\hat{W}^{2}$, through its denominator. Equation (20) suggests that the eddy viscosity depends on the mean strain or rotation rate. The mean strain rate- and vorticity-dependent expression of eddy viscosity is physically natural because $\nu_{\mathrm{T}}$ given by Eq. (6) with a constant $C_{\nu}$ is not universal for several turbulent shear flows, as mentioned in Sec. II.1. This functionalization of the dimensional coefficients leads to the regular behavior of the Reynolds stress for strongly strained or rotating turbulent flows. Moreover, the effect of the additional term $N_{ij}$ deductively enters the denominator of each term on the right-hand side of the Reynolds stress, as shown in Eqs. (18) and (19). The present modeling provides a systematic way to determine a functional form of the dimensional coefficients for more general cases such as one that involves higher-order nonlinearity on the mean velocity gradient.

The dependence of the eddy viscosity on the mean strain or rotation rate was analytically suggested by Okamoto Okamoto (1994) through TSDIA Yoshizawa (1984). In their study, such an effect of the mean velocity gradient on the eddy viscosity appears as part of the cubic-nonlinear term. In the present modeling, it rather comes from the constraint given by Eqs. (15) and (16). The mean strain rate- and vorticity-dependent form of the eddy viscosity was also discussed as the synthesized or composite time scale Yoshizawa et al. (2006). Hamba Hamba (2001) analytically obtained a quadratic-nonlinear realizable model through TSDIA Yoshizawa (1984). In contrast to the theoretical approach by Hamba Hamba (2001), we have to choose a set of terms for $\{T^{\mathrm{m}}_{ij}\}$, such as $S_{ij}$ and $W_{ij}$, in the present modeling. However, this freedom of choice of $\{T^{\mathrm{m}}_{ij}\}$ enables us to easily extend the realizable modeling to more general turbulent flows involving the turbulent helicity effect Yokoi and Yoshizawa (1993); Yokoi and Brandenburg (2016); Inagaki et al. (2017) or the time-history effect of the mean strain rate Speziale (1987); Hamilington and Dahm (2008); Ariki (2019).

III Construction of the model based on the square root technique

III.1 Quartic-nonlinear model of the Reynolds stress

To construct a realizable model applicable to three-dimensional flows, we consider a simple quadratic-nonlinear form of the square root tensor on the velocity gradient:

[TABLE]

Using Eq. (10), the Reynolds stress is written as

[TABLE]

where

[TABLE]

and $C^{2}=C_{ij}C_{ij}$. Note that $S_{ij}C_{ij}=W_{ij}C_{ij}=0$. Because $C_{ij}$ is a quadratic-nonlinear term on the mean velocity gradient, the model given by Eq. (23) corresponds to a quartic-nonlinear eddy-viscosity model. Although Eq. (23) can be put into irreducible form through Cayley-Hamlitom’s theorem as suggested by Pope Pope (1975), we retain Eq. (23) as the present form to simply express coefficients only by $f_{S}$, $f_{W}$, and $f_{C}$. The model given by Eq. (23) comprises seven independent tensor bases except for $\delta_{ij}$. In the conventional turbulence modeling, there is no systematic way to determine the seven-dimensional coefficients, which makes the Reynolds stress satisfy the realizability. The square root modeling reduces the freedom of the dimensional coefficients by imposing the realizability as a guiding constraint. Consequently, we only need to determine three coefficients, $f_{S}$, $f_{W}$, and $f_{C}$, in the present modeling.

To see some basic properties of the model given by Eq. (23), we consider a unidirectional shear flow defined by $\bm{U}=(U_{x},U_{y},U_{z})^{\mathrm{t}}=(U_{x}(y),0,0)^{\mathrm{t}}$. In this flow, the Reynolds stress given by Eq. (23) reads

[TABLE]

where $G=\partial U_{x}/\partial y$ and $R_{xz}=R_{yz}=0$. There are two points that should be noted. The first point is the anisotropy of the Reynolds stress. From many DNS and experiments of shear flows, such as a turbulent channel flow or a homogeneous turbulent shear flow, it is well-known that the diagonal components of the Reynolds stress show the inequality $R_{xx}>R_{zz}>R_{yy}$. However, the model given by Eqs. (25a)–(25c) results in $R_{zz}\leq R_{yy}$ when $f_{C}=0$. Hence, we demand $f_{C}>0$. As seen from Eqs. (25a) and (25b), $C_{ij}$ term plays the role of redistribution of the turbulence intensity between $R_{xx}$ and $R_{yy}$. When $f_{C}>0$, the model given by Eqs. (25a)–(25c) can accurately predict the inequality for turbulent shear flows, $R_{xx}>R_{zz}>R_{yy}$.

The second point is the sign of the shear stress $R_{xy}$. If the shear rate increases and $2f_{W}f_{C}G^{2}$ becomes larger than $f_{S}$, $f_{S}<2f_{W}f_{C}G^{2}$, the sign of $R_{xy}$ becomes the same as the mean velocity gradient $G$. This implies that the effective viscosity can be negative. The negative effective viscosity not only causes numerical instability but also gives a physically invalid condition for the unidirectional shear flow of the high-shear rate condition. Hence, we demand $f_{S}>2f_{W}f_{C}G^{2}$ for an arbitrary shear rate $G$ and adopt an appropriate expression for $f_{W}$ or $f_{C}$.

Considering the above discussion, we propose a $K$-$\varepsilon$ model expression for the present model. Here, the following expressions for the coefficients $f_{S}$, $f_{W}$, and $f_{C}$ are adopted:

[TABLE]

where $f_{\nu}$ is a damping function introduced for the connection of the model to the solid wall and $\hat{S}^{2}$ and $\hat{W}^{2}$ are already defined in connection to Eq. (20). $C_{1}$, $C_{2}$, and $C_{3}$ are positive model constants. Because $f_{C}$ depends on the mean strain and rotation rate, it can be calibrated to satisfy $f_{S}>2f_{W}f_{C}G^{2}$ for a high-shear-rate condition of the unidirectional shear flow. Of course, Eqs. (26a) and (26b) are not a unique choice of $f_{S}$, $f_{W}$, and $f_{C}$, but are a simple form that avoids the irregular condition for the shear stress $R_{xy}$. In addition, we do not pay much attention to the behavior of the model around TCL condition in the present study. To consider TCL condition of ARSM in the vicinity of the solid wall, we may have to put additional terms expressing the near wall effect Abe et al. (2003). Such a modeling of TCL physics is a future study. In the present study, we rather focus on the bulk behavior away from the solid wall.

III.2 Anisotropy tensor in unidirectional shear flow

We examine the anisotropic behavior of the model given by Eqs. (25a)–(25e), (26a), and (26b) for a high-shear-rate case in the unidirectional shear flow. The anisotropy tensor of the Reynolds stress $b_{ij}$ is defined as

[TABLE]

The anisotropy tensor for the model given by Eqs. (25a)–(25e), (26a), and (26b) is written as

[TABLE]

where $\hat{G}=KG/\varepsilon$ and $b_{xz}=b_{yz}=0$. Note that the damping function $f_{\nu}$ is not concerned here; namely, $f_{\nu}=1$. In this case, $b_{ij}$ is determined only by the non-dimensional shear rate $\hat{G}$. It is clearly seen that the condition $b_{zz}>b_{yy}$ is satisfied for finite $\hat{G}$ when $C_{2}/(2C_{3})<1$. Moreover, the sign of $b_{xy}$ is always opposite to that of $\hat{G}$ as long as $C_{2}/(2C_{3})<1$. Figure 1 shows the profile of the anisotropy tensor against the non-dimensional shear rate $\hat{G}$. Here, $C_{1}=0.13$, $C_{2}=0.021$, and $C_{3}=0.018$ are adopted. The first equation of the realizability condition given by Eq. (1a) is rewritten as

[TABLE]

It is confirmed that the present model exactly satisfies Eq. (29) and reproduces the physical features in shear flows, namely $b_{xx}>b_{zz}>b_{yy}$ and $b_{xy}<0$, even when the non-dimensional shear rate $\hat{G}$ is large.

III.3 Shear stress in rotating shear flow

For another simple shear flow case, we consider a spanwise rotating shear flow Bardina et al. (1983). We consider the same unidirectional shear flow as previous sections and set the rotation axis in the $z$ direction, $\bm{\Omega}^{\mathrm{F}}=(0,0,\Omega^{\mathrm{F}})^{\mathrm{t}}$ where $\bm{\Omega}^{\mathrm{F}}$ denotes the angular velocity of frame rotation. In this condition, $b_{xy}$ is written as

[TABLE]

Here, a new parameter, $\Omega^{\mathrm{F}}/G$, appears compared with Eqs. (28d) and (28e). Of course, Eqs. (30a) and (30b) reduce to Eqs. (28d) and (28e) when $\Omega^{\mathrm{F}}/G=0$. Figure 2 shows the profiles of $b_{xy}$ against the non-dimensional shear rate $\hat{G}$ for $\Omega^{\mathrm{F}}/G=0,0.25$, and $-0.25$. When the rotation is against the shear, $\Omega^{\mathrm{F}}/G=0.25$, the shear stress $b_{xy}$ decreases, while it increases when the rotation is with the shear, $\Omega^{\mathrm{F}}/G=-0.25$. This result is consistent with the large-eddy simulation performed by Bardina et al. Bardina et al. (1983).

IV Numerical verification

In this section, we numerically verify the performance of the present model in fundamental turbulent shear flows. First, we perform a turbulent channel flow as a basic wall shear flow. Second, we perform a homogeneous turbulent shear flow as an example where the linear eddy-viscosity model given by Eqs. (4) and (6) fails to predict the evolution of the turbulent kinetic energy. Third, we perform an axially rotating turbulent pipe flow as an example of the three-dimensional flow. In this flow, no quadratic-nonlinear eddy-viscosity model can predict the mean tangential or swirl velocity Speziale et al. (2000). In contrast to the conventional realizable model, such as that given by Shih et al. Shih et al. (1993, 1995), the present model involves cubic nonlinearity on the velocity gradient, as seen in Eq. (23), and thus, is able to predict the mean tangential or swirl velocity.

IV.1 Turbulent channel flow

In the turbulent channel flow, the governing equations for the $K$-$\varepsilon$ model are given by

[TABLE]

where $x$ and $y$ denote the streamwise and wall-normal directions, respectively. The length and velocity are normalized by the channel half width $h$ and the wall friction velocity $u_{\tau}(=\sqrt{\nu|\partial U_{x}/\partial y|_{\text{wall}}|})$, respectively. $\sigma_{K}$, $C_{\varepsilon 1}$, $C_{\varepsilon 2}$, and $\sigma_{\varepsilon}$ are the model constants. $f_{\varepsilon}$ denotes the damping function for the destruction rate of the dissipation rate $\varepsilon$. The performance of the present model is compared with that of the linear eddy-viscosity model, in which the Reynolds stress is given by Eq. (4). For the linear eddy-viscosity model, which is applicable to the non-slip boundary condition, we use that given by Abe et al. Abe et al. (1994) (hereafter the AKN model). The results are compared with those of the DNS performed by Moser et al. Moser et al. (1999) at $\mathrm{Re}_{\tau}(=u_{\tau}h/\nu)=590$. We adopt the same damping functions for the present model as those for the AKN model, which are described by the distance from the wall normalized by the Kolmogorov length scale $\eta[=(\nu^{3}/\varepsilon)^{1/4}]$ and the turbulent Reynolds number $\mathrm{Re}_{\mathrm{T}}(=K^{2}/\nu\varepsilon)$:

[TABLE]

An advantage of using the Kolmogorov length scale for the damping functions is that they can have a non-zero value even when the mean velocity gradient at the wall is zero, $\partial U_{x}/\partial y|_{\text{wall}}=0$, while the conventional damping functions described by the distance from the wall normalized by $\nu/u_{\tau}$ gives $f_{\nu}=0$. The eddy viscosity $\nu_{\mathrm{T}}$ in Eqs. (32) and (33) is given by

[TABLE]

For the linear eddy-viscosity model given by Eq. (4), we use the same expression as Eq. (36). In the present model, the same model constants as those in the AKN model are adopted and these values are given in Table 1. For the unidirectional shear flows, the only difference between the present and AKN models is the model of the Reynolds shear stress $R_{xy}$; in the present model, $R_{xy}$ is given by Eqs. (25d), (26a), and (26b), while in the AKN model, it is given by

[TABLE]

Owing to the damping function $f_{\nu}$, both the present and AKN models achieve the near wall asymptotic condition for the Reynolds shear stress, $R_{xy}\sim y^{3}$.

The profile of the mean velocity is shown in Fig. 3(a). Both the present and linear models give a good result compared to DNS. Figure 3(b) and (c) show the profiles of the turbulent energy and its dissipation rate, respectively. The prediction of the present model is almost similar to that of the linear eddy-viscosity model. In the near wall region at $y^{+}<20$, both the present and linear models slightly underestimate the turbulent energy $K$ and overestimate the dissipation rate $\varepsilon$. In this sense, some modification is required for predicting the near wall behavior. Overall, however, the present model gives a good prediction for the basic properties in the turbulent channel flow. Note that the damping function works only at $y^{+}<100$, so that the behavior of outer region, $y^{+}>100$, is independent of the damping function.

In the RANS model of the unidirectional shear flows, the diagonal components of the Reynolds stress do not appear in the governing equations. In a turbulent flow in a square duct, however, the anisotropic property of the Reynolds stress is essential for predicting the secondary flow Speziale (1987). In this sense, it is worth examining the anisotropic property of the model for a wall-bounded turbulent flow. It is well-known that the linear eddy-viscosity model given by Eq. (4) cannot predict the anisotropy of the Reynolds stress in the unidirectional shear flow; namely, the linear eddy-viscosity model leads to

[TABLE]

while the experiments and DNS show $R_{xx}>R_{zz}>R_{yy}$ in the turbulent channel flow. As discussed in Sec. III, the present model can predict this inequality. The profile of each component of the Reynolds stress for the present model is shown in Fig. 4(a). The profile of the anisotropy tensor of the Reynolds stress $b_{ij}$ is also shown in Fig. 4(b). It is seen that the present model well-predicts the anisotropic property at $100<y^{+}<400$. In the near wall region at $y^{+}<50$, the present model underestimates the streamwise component of the Reynolds stress, as seen in Fig. 4(a) and (b). This is mainly caused by the failure of prediction of the non-dimensional shear rate $\hat{G}$. The profile of $\hat{G}$ is shown in Fig. 5. It is seen that $\hat{G}$ at $y^{+}\sim 10$ is about two thirds the DNS for the present model. This tendency is the same as that of the linear model. This result comes from the underestimation of $K$ and the overestimation of $\varepsilon$ at $y^{+}\sim 10$, because the mean velocity is well-predicted by both the present and linear models. Hence, we need to modify the transport equations for $K$ and $\varepsilon$ at the near wall region to improve the accuracy of the model. However, we do not treat this point further.

IV.2 Homogeneous turbulent shear flow

As discussed in Sec. II.1, the standard type of the linear eddy-viscosity model given by Eqs. (4) and (6) overestimates the time evolution of the turbulent energy in a homogeneous turbulent shear flow. As the conventional ARSM overcomes this shortfall with the aid of the mean strain or rotation rate dependence of the dimensional coefficients Taulbee (1992); Gatski and Speziale (1993); Girimaji (1996); Abe et al. (1997); Shih et al. (1993, 1995); Hamba (2001); Yoshizawa et al. (2006), the present model also does. The governing equations for the $K$-$\varepsilon$ model are given by Eqs. (32) and (33), where the diffusion terms, the third term in them, vanish and $f_{\nu}=f_{\varepsilon}=1$ in the homogeneous turbulent shear flow. The reference DNS was performed by Hamba Hamba (2001), where the mean shear rate was set to $G=7.67$, the initial value of the Taylor microscale Reynolds number $\mathrm{Re}_{\lambda}(=\sqrt{20K^{2}/3\nu\varepsilon})$ was set to $\mathrm{Re}_{\lambda}=24$, the initial value of the turbulent energy was set to $K=1/2$, and the initial energy spectrum was proportional to $k^{4}\exp[-2k^{2}/k_{\mathrm{p}}^{2}]$, where $k_{\mathrm{p}}=9$. As the initial conditions for the $K$-$\varepsilon$ model, the DNS values of $K$ and $\varepsilon$ at time $Gt=6$ are adopted, where the energy spectrum is fully developed in the high-wavenumber region.

The time evolution of the turbulent energy is shown in Fig. 6(a), which also plots the result of the linear eddy-viscosity model for reference, where the set of the model constants is the same as that for the AKN model Abe et al. (1994) given in Table 1. As seen in Fig. 6(a), the present model well-predicts the time evolution of the turbulent energy, while the linear model overestimates it. The time evolution of the anisotropy tensor of the Reynolds stress is shown in Fig, 6(b). The good agreement of $b_{xy}$ with the DNS value provides the good prediction of the turbulent energy. However, the diagonal components of the anisotropy tensor are somewhat overestimated. This result indicates that the present model tends to overestimate the anisotropy for high-shear-rate turbulent flows.

In the present study, we construct a model for $\sqrt{R}_{ij}$ in terms of $S_{ij}$, $W_{ij}$, and $C_{ij}$ for simplicity of the model expression. Hence, there is still enough room for further development, e.g., introduction of additional terms such as $S_{ia}S_{aj}$, $W_{ia}W_{aj}$, or the time history effect of the mean strain tensor Speziale (1987); Hamilington and Dahm (2008); Ariki (2019) to $\sqrt{R}_{ij}$, or modifications of $f_{S}$, $f_{W}$, and $f_{C}$. Namely, the overestimation of anisotropy in the homogeneous turbulent shear flow can be improved by more sophisticated modeling.

IV.3 Axially rotating turbulent pipe flow

A prominent feature of the present model is that it involves cubic nonlinearity on the mean velocity gradient, which is an essential for predicting the mean swirl flow in an axially rotating turbulent pipe flow Speziale et al. (2000). This point is a critical difference between the present and the conventional realizable models Shih et al. (1993, 1995); Abe et al. (1997); Hamba (2001). Assuming the homogeneity of the turbulence field in the axial and azimuthal directions, the mean velocity can be written as $\bm{U}=(U_{r},U_{\theta},U_{z})=(0,U_{\theta}(r),U_{z}(r))$ in the cylindrical coordinates. The governing equations for the $K$-$\varepsilon$ model in an inertial frame of the cylindrical coordinates are written as follows:

[TABLE]

Here, the velocity and length scale are normalized by the bulk mean velocity $U_{\mathrm{m}}[=(2/R^{2})\int_{0}^{R}\mathrm{d}r\ rU_{z}]$ and the pipe radius $R$, respectively. $f^{\mathrm{ex}}$ denotes the external forcing that keeps the bulk Reynolds number $\mathrm{Re}_{\mathrm{m}}(=U_{\mathrm{m}}R/\nu)$ constant. In the cylindrical coordinates, we can write a second-rank tensor $\bm{\mathsf{T}}$ as

[TABLE]

In an axially rotating turbulent pipe flow, the mean strain rate and mean absolute vorticity tensors are respectively written as follows:

[TABLE]

In this matrix form, the tensor $\bm{\mathsf{C}}$ defined by Eq. (22) can be calculated as $\bm{\mathsf{C}}=\bm{\mathsf{S}}\bm{\mathsf{W}}+(\bm{\mathsf{S}}\bm{\mathsf{W}})^{\mathrm{t}}(=\bm{\mathsf{S}}\bm{\mathsf{W}}-\bm{\mathsf{W}}\bm{\mathsf{S}})$. Although the model expression of the Reynolds stress seems to be complicated, as seen in Eq. (23), it can be easily calculated from the square root tensor as $\bm{\mathsf{R}}=\sqrt{\bm{\mathsf{R}}}\sqrt{\bm{\mathsf{R}}}^{\mathrm{t}}$ with Eq. (21). This is another advantage of the usage of the square root tensor. We compare the results of the present model with the experimental data provided by Imao et al. Imao et al. (1996), where the bulk Reynolds number is set to $\mathrm{Re}_{\mathrm{m}}=10000$. The rotation rate $N$ is defined as $N=U_{\theta,\text{wall}}/U_{\mathrm{m}}$ and experiments of three parameters, $N=0$, $0.5$, and $1$, were performed.

Figure 7(a) shows the profile of the mean axial velocity $U_{z}$ in a steady state for each rotation parameter. The result of the linear eddy-viscosity model is also plotted for reference. It is seen that the present model predicts the rotation-dependent property of the mean axial velocity; namely, the mean velocity gradient becomes steep as the rotation parameter increases. Figure 7(b) shows the profile of the mean swirl velocity $U_{\theta}$ in a steady state for two rotating cases. Moreover, neither the linear eddy-viscosity model nor quadratic-nonlinear models can predict the mean swirl velocity Speziale et al. (2000); they just predict the solid body rotation solution, $U_{\theta}=rU_{\theta,\text{wall}}/R$. As seen in Fig. 7(b), the present model predicts the rotation-dependent mean swirl velocity. Although both trends of the mean axial velocity and mean swirl velocity are reproduced by the present model, there is still room for improvement in the quantitative sense. This point should be further discussed in the future. However, it should be rather emphasized that the present model is a realizable quartic-nonlinear ARSM which has never been proposed. It can reproduce a feature of three-dimensional flow owing to the quartic nonlinearity on the mean velocity gradient with ensuring the realizability conditions.

V Conclusions

In this study, the realizable ARSM based on the square root tensor of the Reynolds stress Ariki (2015) was further developed. By imposing the realizability conditions, we could develop a physically consistent and numerically stable turbulence model. In contrast to the conventional method, the square root modeling provides a systematic and simpler construction of the realizable model, even when it involves higher-order nonlinearity on the mean velocity gradient. The cubic or higher-order nonlinearity on the mean velocity gradient is required to predict three-dimensional flows, which denote the turbulent flows with three-dimensional mean velocity. In fact, the cubic nonlinearity is required to predict the mean swirl flow in an axially rotating turbulent pipe flow Speziale et al. (2000).

For the development from the previous study by Ariki Ariki (2015), we proposed a systematic way to determine the dimensional coefficients of the model. Consequently, a dependence of the eddy viscosity on the mean strain and rotation rate is deductively introduced. Moreover, the dimensional coefficients depend on the other additional terms chosen for an expansion basis of the square root tensor. In this sense, the present modeling gives a generic way to determine a functional form of the dimensional coefficients.

As a simple model applicable to three-dimensional flows, a quartic-nonlinear ARSM was proposed through the square root tensor. The performance of the proposed model was numerically verified in a turbulent channel flow, a homogeneous turbulent shear flow, and an axially rotating turbulent pipe flow. The present model gave reasonable predictions for both mean velocity profile in the turbulent channel flow and time evolution of the turbulent energy in the homogeneous turbulent shear flow. Moreover, the rotation rate-dependent property of the mean axial and swirl velocity in the axially rotating turbulent pipe flow was predicted, which cannot be predicted by any quadratic-nonlinear eddy-viscosity model, including the conventional realizable models Shih et al. (1993, 1995); Abe et al. (1997); Hamba (2001); Abe et al. (2003).

The present model overestimated the anisotropy of the Reynolds stress in the homogeneous turbulent shear flow and underestimated the mean axial and swirl velocity in the axially rotating turbulent pipe flow; namely, there remains enough room for further improvement in the quantitative sense. For simplicity of the model, we did not consider the effect of the turbulent helicity Yokoi and Yoshizawa (1993); Yokoi and Brandenburg (2016); Inagaki et al. (2017) and the time history effect of the mean strain rate Speziale (1987); Hamilington and Dahm (2008); Ariki (2019) in this study. Hence, the model can be further developed. However, note that the present model is a quartic-nonlinear ARSM that always satisfies the realizability conditions, which has never been proposed. Such a higher-order realizable turbulence model, involving quartic nonlinearity on the mean velocity gradient, is expected to be useful in numerically stable predictions of turbulent flows with three-dimensional mean velocity.

Appendix A Physical interpretation of the square root tensor

Let us introduce an isotropic stochastic vector $\xi_{i}$, which satisfies $\langle\xi_{i}\rangle=0$ and $\langle\xi_{i}\xi_{j}\rangle=\delta_{ij}$. We consider the anisotropic velocity fluctuation $u_{i}^{\prime}$ expressed by the following mapping form:

[TABLE]

where $\sqrt{R}_{ij}$ denotes the anisotropic mapping operator. Note that we consider a non-fluctuating $\sqrt{R}_{ij}$; namely, $\langle\sqrt{R}_{ij}\rangle=\sqrt{R}_{ij}$ and $\langle\sqrt{R}_{ij}\xi_{j}\rangle=\sqrt{R}_{ij}\langle\xi_{j}\rangle=0$. Then, the Reynolds stress is calculated by its definition, $R_{ij}=\langle u_{i}^{\prime}u_{j}^{\prime}\rangle$, as

[TABLE]

This is just the definition of the square root tensor of the Reynolds stress given by Eq. (10). Hence, $\sqrt{R}_{ij}$ can be interpreted as an envelope representing the anisotropy of the velocity fluctuation.

Here, we discuss an expression for the velocity fluctuation in anisotropic turbulence in terms of the mean velocity gradient. For theoretical convenience, we expand the velocity fluctuation around the isotropic part as $u_{i}^{\prime}=u^{(0)}_{i}+u^{(1)}_{i}$, where $u^{(0)}_{i}$ denotes the isotropic velocity fluctuation and $u^{(1)}_{i}$ denotes the anisotropic velocity fluctuation around $u^{(0)}_{i}$. The perturbation equation for $u^{(1)}_{i}$ from isotropic turbulence can be written as Yoshizawa (1984, 1998, 1995)

[TABLE]

where

[TABLE]

Equation (A3) can be formally integrated as

[TABLE]

This inverse operation $L_{u}^{-1}$ denotes the time integration with a characteristic time scale of turbulence and we approximate this integration as Yoshizawa (1995, 1998)

[TABLE]

where $\tau$ represents a characteristic time scale of turbulence. When the isotropic velocity fluctuation is written as $u^{(0)}_{i}=\sqrt{2K_{0}/3}\xi_{i}$, the velocity fluctuation $u_{i}^{\prime}$ reads

[TABLE]

Comparing Eqs. (A1) and (A7), the square root tensor $\sqrt{R}_{ij}$ can be expressed as

[TABLE]

This expression of $\sqrt{R}_{ij}$ corresponds to Eq. (17), where $\gamma_{0}=\sqrt{2K_{0}/3}$, $\gamma_{S}=\gamma_{0}f_{S}$, $\gamma_{W}=\gamma_{0}f_{W}$, $\gamma_{N}=0$, and $f_{S}=f_{W}=\tau/2$. The higher-order nonlinear term can be obtained by considering the higher-order perturbation. Consequently, the modeling based on the square root tensor is interpreted as the modeling based on the generation mechanism of the anisotropy of the velocity fluctuation.