explicit algebraic Reynolds stress models.
By Igor A. Grigoriev
1& Werner M.J. Lazeroms
1,21Linn´e FLOW Centre, KTH Mechanics, SE-100 44 Stockholm, Sweden
2 Institute for Marine and Atmospheric Research Utrecht, Utrecht University, 3508 TA Utrecht, The Netherlands
Methodological note
A direct solution to a tensorial equation which constitutes a basis for explicit algebraic Reynolds stress models is derived. We consider equations linear and quasilinear in the strain tensor and show how the independent tensor groups emerge. Solution of an extended model with a linearly coupled active scalar, governed by a linear in anisotropy tensor equation, is also outlined.
1. Introduction
In Lazeroms and Grigoriev (2015) a solution to the following quasilinear equa- tion has been obtained and analyzed:
( λx = L(x) + ϕ,
λ = α + f (x). (1)
Specifically, if V is a n-dimensional vector space, then L : V → V is a linear operator and ϕ ∈ V is a given vector and f : V → R is a linear functional on V while α ∈ R is a given constant. To find a solution means that both x ∈ V and λ ∈ R are determined.
Due to Cayley-Hamilton theorem some power n of L, in fact equal to the dimension of the operator, can be expressed as a sum of lower powers:
(cnLn+ cn−1Ln−1+ · · · + c0I)(x) = 0, for all x ∈ V , cn= 1, cn−1= −tr (L), · · · c0= (−1)ndet(L), c−1= 0,
where I is the identity operator and all ckdepend only on L. However, a special form of ϕ may imply that x also has a specific structure and on this contracted set an mthpower of L(x) can be expressed through a sum of lower powers with m < n. Consequently, computing the powers of the first of equations in (1) we can show that it is possible to solve the equation for x and form an (n + 1)th
191
(or (m + 1)th) order polynomial equation for λ. Indeed, from
L1(x) = L(x) = λx − ϕ, L2(x) = λ2x − λϕ − L(ϕ), · · ·
· · · Ln(x) = λnx − λn−1ϕ − λn−2L(ϕ) − . . . − Ln−1(ϕ), the general solution can be derived
x =
n−1
X
k=0
βkLk(ϕ), βk= ck+1+ ck+2λ + · · · + cnλn−k−1 c0+ c1λ + · · · + cnλn , λn+1+ dnλn+ dn−1λn−1+ · · · + d0= 0,
(2)
where dk = ck−1− αck−Pn−k−1
j=0 cj+k+1f [Lj(ϕ)].
For our purposes, we may consider the above method as a generalization of an arbitrary system of quasilinear equations for a tensor aij and a vector
˚ζi which treats both quantities on an equal footing through a vector x. Here we have more modest aim and we will show how to directly solve an algebraic tensor equation
N aij= −C1Sij+ aikΩkj− Ωikakj+ + C2
aikSkj+ Sikakj−2
3(akmSmk) δij
, (3) where Sij is a traceless symmetric tensor, Ωij – antisymmetric tensor, while C1, C2and N are constants.
Equation (3) is encountered in explicit algebraic Reynolds stress modeling of so called ’weakly-equilibrium’ states of turbulence (Wallin and Johansson (2000)). Then aij, Sij, Ωij represent the anisotropy tensor, strain tensor and vorticity tensor, respectively, while N should be additionally determined by solving an equation
N = c01− C3tr (aijSjk) (4) with constants c01 and C3. A common method of solving tensorial algebraic equations for the anisotropy tensor, in general non-linear in aij as well as in Sij and Ωij, based on representing aij as a sum of the complete set of tensor groups expressible in terms of algebraic combinations of Sij and Ωij, was proposed in Pope (1975). Essentially this method has been employed in Wallin and Johansson (2000) to solve the equation (3). Grundestam et al.
(2005) used a modified version of this method, basing on an incomplete set of tensor groups, to solve the nonlinear (both in aij and Sij) extensions of (3) using e.g. extensions by Hallb¨ack et al. (1990), Johansson and Hallb¨ack (1994) and Sj¨ogren and Johansson (2000).
A system including a coupling with an active scalar has been considered in Lazeroms et al. (2013) who generalized Pope’s method by building a broader tensor basis for a two-dimensional flow. Then, Grigoriev et al. (2015) have
found a general three-dimensional solution for the coupled system system sim- ilar to Lazeroms et al. (2013) (assuming C2 = 0), without constructing an extended tensor basis but appropriately using the original Pope’s basis.
Observe that if aij is an algebraic combination of Sij, Ωij and N , the equation (4) becomes a polynomial equation for N . Here we do not analyze its solution and the choice of the physical root. Here we will only present a direct solution to the equation (3) (and discuss a possible extension to a coupled active scalar equation). By itself, it will not give any new knowledge but it is instructive to develop a direct method. As it happens, methods often turn out to be more important than the immediate result they give.
The only preliminary knowledge we need are the consequences of Cayley- Hamilton theorem. Namely, we note that for arbitrary symmetric traceless tensor s and arbitrary antisymmetric tensor o (for brevity we will skip both indices and bold font for tensors later in the paper)
(s + o)2= s2+ s o + o s + o2, tr (s + o)2= IIs+ Io, Io= tr o2, IIs= tr s2, (s + o)3= s3+ (s2o + s o s + o s2) + (s o2+ o2s + o s o) + o3,
tr (s + o)3= III + 3 IV, IV = tr (s o2), III = tr s3, V = tr (s2o2). (5) Then,
det (s + o)3= tr (s + o)3
3 , (s + o)3= IIs+ Io
2 (s + o) + III 3 + IV
I3, (6)
where I3 is the three-dimensional unit tensor, and equating the equal powers of s and o we find
s o2+ o2s + o s o =Io
2 s + IV I3, s2o + s o s + o s2= IIs 2 o, o3= Io
2 o, s3= IIs
2 s + III 3 I3.
(7)
2. Equation linear in S
ijWith C2= 0 the equation (3) can be rewritten as simple as
a = −s + (a o − o a), s = C1S/N, o = Ω/N. (8) Multiplying (8) by o2, taking the trace and using an analogy to the first relation of (7) for the symmetric tensor a we arrive at
tr (a o2) = tr
[−s + a o − o a] o2
= −IV, a o2+ o2a + o a o = Io
2 a − IV I3. (9)
The first step recursive employment of (8) aided by (9) leads to
a = −s − (s o − o s) +
(a o − o a) o − o (a o − o a)
=
= −s−(s o−o s)+
a o2+o2a−2 o a o
= −s−(s o−o s)−3 o a o−IV I3+Io 2 a,
1 − Io
2
a = −s − (s o − o s) − IV I3− 3 o a o . (10)
The second recursive usage of (8) transforms the only unknown o a o to o a o = − o s o + o (a o − o a) o = − o s o + o2s o − o s o2+
+
o (o2a + a o2− 2 o a o) o
, (11) and computing the last square bracket using (9) and (7) we obtain an explicit expression for o a o
...
= o (3 a o2+ 3 o2a − Ioa + 2 IV I3) o = 2 IV o2+ 2 Ioo a o,
(1 − 2 Io) o a o = − o s o + o2s o − o s o2+ 2 IV o2 . (12) Multiplying (10) by (1 − 2 Io) we obtain the answer
1−Io
2
1−2 Io
a = −
1−2 Io
s−
1−2 Io
(s o−o s)−
1−2 Io
IV I3+ + 3 o s o − 3 (o2s o − o s o2) − 6 IV o2= −
1 −7
2Io
s −
1 − 2 Io
(s o − o s)−
− 6 IV
o2−Io
3 I3
− 3
s o2+ o2s − 2 3IV I3
− 3 (o2s o − o s o2). (13) It can be transformed to its common form by Wallin and Johansson (2000) substituting s = C1S/N and o = Ω/N .
3. Equation quasilinear in S
ijIf C26= 0 equation (3) can be rewritten as a = −s + (a o − o a) + α
a s + s a −2
3tr (a s) I3
, α = C2/C1. (14) This general case seems to be much less viable to be solved directly than (8) since (14) is quasilinear in s and allegedly all the tensor groups based on s and o will be involved in the solution, not only linear in s. Nevertheless, the steps leading to the solution are essentially the same as (10), (11-12) and (13) where
relation similar to (7) and (9) have to be employed. Putting (14) as
˜
a = −s + ˜a o1− o2˜a, ˜a =
a + α2
3tr (a s) I3
1 + α24 3tr (a s)
,
o1= o + α s, o2= o − α s, (15)
we note that formally the complications are due to the use of mixed (nor an- tisymmetric neither symmetric) tensors o1 (’right’ o) and o2 (’left’ o) instead of single o. The employment of non-traceless pseudo-anisotropy tensor ˜a (com- prising a and unity tensor proportional to α tr (a s)) does not complicate the consideration by itself. Indeed, using (7) we arrive at similar relations involving s, o1 and o2
o31=IIα
2 o1+ α IVαI3, o2s o1+ o22s + s o21= IIα
2 s + IVαI3, o32= IIα
2 o2− α IVαI3, IIα= Io+ α2IIs, IVα= IV + α2III
3 . (16)
We cannot obtain the relation corresponding to (9) just by substituting a to the second equality in (16). Instead, recursively using (15) with the aid of (16) we derive
o2˜a o1= −o2s o1+ o2a o˜ 21− o22a o˜ 1= −o2s o1+ (˜a o1− ˜a − s) o21+ + o22(−o2˜a − ˜a − s), o2a o˜ 1+ ˜a o21+ o22˜a =
= −o2s o1− s o21− o22s +IIα
2 s + IIα
2 + 2 α IVα
˜ a, o2˜a o1+ ˜a o21+ o22˜a = IIα
2 + 2 α IVα
˜
a − IVαI3. (17)
Incorporating the last formula into the one step recursion of (15) we arrive at the exact analogue of (10) (only the factor in square brackets becomes extended)
a = −(s + s o˜ 1− o2s) + ˜a o21− 2 o2˜a o1+ o22˜a,
1 − IIα
2 + 2 α IVα
˜ a = −
s + s o1− o2s
− IVαI3− 3 o2a o˜ 1 . (18)
In its turn, the exact analogue to (11) is obtained by using (15) and looks like
o2˜a o1= o2
− s − s o1+ o2s
o1+
o2
˜
a o21− 2 o2˜a o1+ o22˜a
o1
, (19)
which due to (16-17) leads to
...
= 3 o2
˜
a o21+ o22˜a
o1− 2 o2
˜
a o21+ o2˜a o1+ o22˜a
o1=
= 3 IIαo2˜a o1+ 3 α IVα(o2˜a − ˜a o1) − 2 o2 IIα
2 + 2 α IVα
˜
a − IVαI3
o1=
= −3 α IVα˜a − 3 α IVαs + 2 o2IVαo1+
2 IIα− 4 α IVα
o2˜a o1, (20)
1 −
2 IIα− 4 α IVα
o2˜a o1= o2(−s − s o1+ o2s) o1+ 2 IVαo2o1−
− 3 α IVαs − 3 α IVα˜a.
This expression is inexact analogue of (12) due to two terms proportional to α on the right-hand side. Especially, the last term presents a complication showing that o2˜a o1 also formally depends on ˜a.
The system to be solved consists of equations (18) and (20) which we rewrite here
1 − IIα
2 + 2 α IVα
˜ a = −
s + s o1− o2s
− IVαI3− 3 o2˜a o1,
1 −
2 IIα− 4 α IVα
o2˜a o1= o2
− s − s o1+ o2s
o1+ 2 IVαo2o1−
− 3 α IVαs − 3 α IVα˜a.
(21)
Multiplying the first of the equations by [1 − (2 IIα− 4 α IVα)], substituting o2˜a o1 from the second one and moving the term with ˜a to the left-hand side we arrive at the following result:
1 −
2 IIα− 4 α IVα
1 − IIα
2 + 2 α IVα
− 9 α IVα
˜a =
= −
1 −
2 IIα− 4 α IVα
s + s o1− o2s
+ 3 o2
s + s o1− o2s
o1−
− 6 IVαo2o1+ 9 α IVαs −
1 −
2 IIα− 4 α IVα
IVαI3, (22) which is a generalized analogue of (13).
Taking into account that o s o s + 2 o s2o + s o s o + o2s2+ s2o2= IIso2 (a consequence of the second relation in (7)) we can rewrite the right-hand side of (22) as a sum of the traceless symmetrical tensor groups and unity tensor multiplied by a function of tensor invariants. Then, we split the equation (22) into two parts equating the traceless part of ˜a =
1 + α2 43tr (a s)
−1 a +
α23tr (a s) I3
to the former constituent on the right-hand side and the re- maining part – to the latter. The last equation, after factorization by 23α I3, is
1−
2 IIα−4 α IVα
1− IIα
2 +2 α IVα
−9 α IVα
tr (a s) 1 + α24
3tr (a s)
=
= −
IIs
1 −7
2IIα
+ 6 IVα2+ 6 V + 2 α
III − IIsIVα
+ 3 α2IIs2
, (23) from which we immediately obtain the expression for tr (a s) as
tr (as) = −∆
1 + α24 3∆
,
∆ = IIs
1 −7
2IIα
+ 6 IVα2+ 6 V + 2 α
III − IIsIVα
+ 3 α2IIs2
1 −
2 IIα− 4 α IVα
1 − IIα
2 + 2 α IVα
− 9 α IVα
. (24)
Now, the expression for a is
1 + α24
3∆
1 −
2 IIα−4 α IVα
1 − IIα
2 + 2 α IVα
−9 α IVα
a =
= −
1 + 3 α2 IIs 2 + α2
3III
−
2 IIα+ 5 α IVα
s−
−
1 + 3 α2IIs−
2 IIα+ 2 α IVα
s o − o s
+
+ 3
o s o − IV
3 I3+ o (s o − o s) o
− (6 IVα− 3 α IIs)
o2+ α2s2−IIα
3 I3
−
− 2 α
1 + 3 α2IIs−
2 IIα+ 2 α IVα
s2−IIs 3 I3
−
− 3 α
s2o − o s2+ 2
s2o2+ o2s2−2 3V I3
− 3 α2
s2o s − s o s2
. (25)
Thus, we have obtained all the tensor groups in the course of direct solution.
The proposed method also allowed us to compute tr (a s) automatically which can be used to verify the result. The only independent tensor group, expressible in s and o, which is not in (25) is o s2o2− o2s2o, consistent with Wallin and Johansson (2000).
4. Extension for the case of coupling with active scalar
For the case when an active scalar ˚ζi is coupled to the anisotropy equation in a simple linear way through a forcing factor Υj, we have to rewrite (3) and (4) as
N aij= −C1Sij+ aikΩkj− Ωikakj+ C2
aikSkj+ Sikakj−
−2
3(akmSmk) δij
− CζSij(+), Sij(+)=
˚ζiΥj+ ˚ζjΥi−2 3
˚ζkΥkδij
, (26) with a constant Cζ, and
N = c01− C3
tr (aijSji) + ˚ζkΥk
. (27)
The solution of (26) has to be linear in ˚ζi and Υj, which has been used in Grigoriev et al. (2015) for the case with C2= 0. When C26= 0 the anisotropy tensor is not linear in Sij and we cannot treat the last tensor group in (26) on an equal footing with Sij. Fortunately, for solving the quasilinear equation with the additional tensor group Sij(+)we only need to slightly modify the appoach in section 3. Indeed, it is easy to show that in this case (15) has to be modified as
˜
a = −˜s + ˜a o1− o2˜a, ˜s = s + s(+) 1 + α24
3tr (a s)
, s(+)=CζSij(+)
N ,
o1= o + α s, o2= o − α s, ˜a =
a + α2
3tr (a s) I3 1 + α24
3tr (a s)
. (28)
Following the solution procedure (16)-(22) but with ˜s not s as a free-standing quantity on the right-hand side of (28), we see that the structure of o1- and o2- interaction with ˜a does not change and only additional linear terms proportional to s(+)/[1 + α24/3 tr (a s)] emerge in all transformed versions of the linear equation for ˜a. Thus, the resulting analogue of solution (22) has the described terms only on the right-hand side. Multiplying both sides of the equation by [1 + α24/3 tr (a s)] and taking the part I3we will obtain linear in s(+)equation for tr(a s). Using the result for equating the remaining traceless part, we will arrive at the expression for a which is also linear in s(+).
Having found this formal solution which through s(+) depends on the un- determined yet ˚ζi, we can substitute it to the equation for the active scalar. In its simplest form ˚ζi-equation looks like
Nζδik+ cSSik+ cΩΩik
˚ζk = −
aik+2
3δik
Γk−2
3cΥΥi (29) where Nζ = N
2 C3 + C4 and C4 is a constant. The last equation with aij
linear in the active scalar can be easily solved for ˚ζi after some cumbersome but straightforward transformations. Substituting it back to the expression (26) for aij we arrive at the answer and using (27) formulate a polynomial N -equation. Although in a more or less general case this equation is unlikely to be solved analytically, an iterative procedure will fairly quickly lead to a sufficiently precise result, as has been demonstrated for C2= 0 case in Lazeroms et al. (2015) and Grigoriev et al. (2015).
S.B. Pope, “A more general effective-viscosity hypothesis,” J. Fluid Mech. 72, 331-340 (1975).
S. Wallin and A.V. Johansson, “An explicit algebraic Reynolds stress model for in- compressible and compressible turbulent flows,” Journal of Fluid Mechanics 403, 89 (2000).
W.M.J. Lazeroms, G. Brethouwer, S. Wallin and A.V. Johansson, “An explicit alge- braic Reynolds-stress and scalar-flux model for stably stratified flows,” Journal of Fluid Mechanics 723, 91 (2013).
A. V. Johansson and M. Hallb¨ack, “Modelling of rapid pressure-strain in Reynolds- stress closures,” Journal of Fluid Mechanics 269, 143 (1994).
T. Sj¨ogren and A.V. Johansson, “Development and calibration of algebraic non-linear models for terms in the Reynolds stress transport equations,” Phys. Fluids 12, 1554 (2000).
M. Hallb¨ack, J. Groth and A.V. Johansson, “An algebraic model for nonisotropic turbulent dissipation rate in Reynolds stress closures,” Phys. Fluids A 2, 859 (1990).
O. Grundestam, S. Wallin and A.V. Johansson, “Techniques for deriving explicit algebraic Reynolds stress models based on incomplete set of basis tensors and predictions of fully developed rotating pipe flow,” Phys. Fluids 17, 115103 (2005).
I.A. Grigoriev, S. Wallin, G. Brethouwer and A.V. Johansson, “Capturing turbulent density flux effects in variable density flow by an explicit algebraic model,” Phys.
Fluids 27(4), 105112 (2015).
W.M.J. Lazeroms, G. Brethouwer, S. Wallin and A.V. Johansson, “Efficient treatment of the nonlinear features in algebraic Reynolds-stress and heat-flux models for stratified and convective flows,” Int. Journal of Heat and Fluid Flow 53, 15-28 (2015).
W.M.J. Lazeroms and I.A. Grigoriev, “A generalized method for deriving explicit algebraic turbulence models,”published in in PhD Thesis of Werner M.J. Laze- roms Turbulence modelling applied to the atmospheric boundary layer, ISBN 978 − 91 − 7595 − 603 − 9, to be extended.
