2 The incompressible Navier-Stokes equations

(1)

faculty of mathematics and natural sciences

Modeling turbulent flow

Bachelor Project Applied Mathematics

July 2016 Student: Leo Sok

First supervisor: Prof.dr.ir. R.W.C.P. Verstappen Second supervisor: Dr. M.A. Grzegorczyk

(2)

Abstract

Fluid flow is mathematically modelled by the Navier-Stokes equations.

In this thesis we will focus at turbulent flows. Therefore we will first look at the incompressible Navier-Stokes equations. To simulate turbulent flows with the incompressible Navier-Stokes, it is important that the kinetic energy of the flow is conserved. Therefore, we will derive the conservation of the kinetic energy. Then we will discretize the Navier- Stokes equations using the finite volume method and derive in this case the kinetic energy. At last we will derive a model called the QR model.

Therefore we start with a large eddy simulation(LES), which means we will add a spatial filter to a solution of the incompressible Navier-Stokes equations. From this LES, we will derive an eddy-viscosity model called the QR model.

(3)

1 Introduction

The Navier-Stokes equations provide a model for turbulent flow. In this thesis we will follow the study of Henry Bandringa [1].

In section 2 we will look at the incompressible Navier-Stokes equations. Be- cause the conservation of the kinetic energy is important, we will compute the evolution of the kinetic energy in time and see that the kinetic energy will decrease.

In section 3 we will discretize the incompressible Navier-Stokes equations and then compute the evolution of the kinetic energy in time of these discretized Navier-Stokes equation. We will see that the equation for this evolution of the kinetic energy the discrete version is of the equation we computed in the first section and that for the discrete version also the kinetic energy will decrease.

In section 4 we will look at large-eddy simulation (LES). For LES, we add a filter to the incompressible Navier-Stokes equations. This filter introduces a closure model τ . We will use scale-truncation to get a condition on this τ . Then we will introduce an eddy viscosity model for τ , and using the condition on τ , we will derive a model called the QR-model [3].

2 The incompressible Navier-Stokes equations

In this section we will see that the incompressible Navier-Stokes equations are

∂u

∂t + (u · ∇)u + ∇p − ν∆u = 0

∇ · u = 0

From these equations, we will derive the evolution of the kinetic energy en see that the kinetic energy decreases.

2.1 The incompressible Navier-Stokes equations

In this subsection we will give the incompressible Navier-Stokes equations.

Therefore we introduce the following symbols: u is the velocity vector in R³, ρ the density, p the pressure and µ the dynamic viscosity. For the incompressible Navier-Stokes equations, we assume that the density, the temperature and viscosity are constant.

So the equation for conservation of mass

∂ρ

∂t + ∇ · (ρu) = 0 becomes, using the constancy of the density,

∇ · u = 0 (1)

The conservation of momentum is given by ρ∂u

∂t + ρ∇ · (u ⊗ u) + ∇p − µ∆u = 0

(5)

where ⊗ indicates the outer product, which is defined as u ⊗ v = uv^TWhen we divide this equation by ρ and define ν = ^µ_ρ and ˆp = ^p_ρ and call this again p, we get

∂u

∂t + ∇ · (u ⊗ u) + ∇p − ν∆u = 0 (2) Now we rewrite ∇ · (u ⊗ u):

∇ · (u ⊗ u) =X

ij

∂

∂xi

(uiuj) =X

ij

∂u_i

∂xi

uj+ ui

∂uj

∂xi

using equation 1,P

i

∂u_i

∂xi = ∇ · u = 0, so we get

∇ · (u ⊗ u) =X

ij

ui

∂uj

∂xi

= (u · ∇)u (3)

So we end up with the following incompressible Navier-Stokes equations:

∂u

∂t + (u · ∇)u + ∇p − ν∆u = 0 (4)

∇ · u = 0 (5)

2.2 The evolution of the kinetic energy

In this subsection we will compute the evolution of the kinetic energy 1

2 d

dt(u, u) = −ν Z

Ω

∇u : ∇u dΩ + boundary terms

and see that if we neglect the boundary terms, the kinetic energy decreases.

We define the inner product (·, ·) as (f, g) = R

Ωf · g dΩ. Then the time derivative kinetic energy is given by

1 2

d dt

Z

Ω

u · u dΩ = 1 2

d

dt(u, u) = 1 2

∂

∂tu, u

+1

2

u, ∂

∂tu

Using (f, g) =R

Ωf · g dΩ =R

Ωg · f dΩ = (g, f ) we get 1

2 d

dt(u, u) = ∂u

∂t, u

Now we rewrite equation 4 as ^∂u_∂t = −(u · ∇)u − ∇p + ν∆u and put it into the above equation. Then we get

1 2

d

dt(u, u) = −((u · ∇)u, u) − (∇p, u) + ν(∆u, u) (6) We will work this out with the following product rule

∇ · (f u) = ∇f · u + f (∇ · u) (7) When we replace f by u · u and use conservation of mass (equation 5), we get

∇ · ((u · u)u) = ∇(u · u) · u

(6)

Using ∇(u · u) · u =P

ij

∂

∂xi(u_ju_j)u_i=P

ij2u_i^∂u_∂x^j

iu_j= 2(u · ∇)u · u, we get 1

2∇ · ((u · u)u) = (u · ∇)u · u Now we take the integral over Ω and get

1 2

Z

Ω

∇ · ((u · u)u) dΩ = Z

Ω

(u · ∇)u · u dΩ

= ((u · ∇)u, u) Using Gauss’s divergence theorem, we get

Z

Ω

∇ · ((u · u)u) dΩ = Z

Γ

(u · u)(u · n) dΓ

where Γ is the boundary of Ω and n is the orthonormal vector at Γ. So when we combine the two above equations, we get

((u · ∇)u, u) =1 2

Z

Γ

(u · u)(u · n) dΓ So equation 6 becomes

1 2

d

dt(u, u) = −1 2

Z

Γ

|u|²(u · n) dΓ − (∇p, u) + ν(∆u, u) (8)

where |u|²= u · u.

Now we set in equation 7 f = p and see that (∇p, u) =

Z

Ω

∇ · (pu) dΩ − (p, ∇ · u) (9) When we now use conservation of mass (equation 5) and Gauss’s divergence theorem, we get

(∇p, u) = Z

Γ

(u · n)p dΓ So equation 8 becomes

1 2

d

dt(u, u) = − Z

Γ

(u · n)(1

2|u|²+ p) dΓ + ν(∆u, u) (10) When we replace in equation 7 f by a tensor A, it becomes

∇ · (Au) = (∇ · A) · u + A : ∇u Where : is the Frobenius inner product, which is defined by

A : B =X

i,j

AijBij = tr(A^TB) = tr(B^TA)

When we fill in A = ∇u, we get

∇ · ((∇u)u) = (∆u) · u + ∇u : ∇u (11)

(7)

Thus

ν(∆u, u) = ν Z

Ω

(∇ · ((∇u)u) − ∇u : ∇u) dΩ When we now apply Gauss’s divergence theorem, we get

ν(∆u, u) = ν Z

Γ

((∇u)u) · n dΓ − ν Z

Ω

∇u : ∇u dΩ

When we now use ((∇u)u) · n = u ·^∂u_∂n, we end up with the following equation for the evolution of the kinetic energy in time

1 2

d

dt(u, u) = − Z

Γ

(u · n)(1

2|u|²+ p) dΓ + ν Z

Γ

u ·∂u

∂n dΓ − ν Z

Ω

∇u : ∇u dΩ (12) So if we use the no-slip boundary condition, the first two terms on the right hand side vanishes and we get

1 2

d

dt(u, u) = −ν Z

Ω

∇u : ∇u dΩ

When we now use ν > 0 and ∇u : ∇u = P

ij

∂ui

∂x_j

²

≥ 0, we see that

1 2 d

dt(u, u) ≤ 0 and we get ¹₂_dt^d(u, u) = 0 only if u is constant, in which case u = 0 due to the no-slip boundary condition. So if u 6= 0, the kinetic energy will decrease.

3 Discretization of the Navier-Stokes equation

In this section we will derive the folliwing discrete incompressible Navier-Stokes equations.

∂ui

∂t Ωi+ X

j∈Ni

uijuij· nijAij+ X

j∈Ni

pijnijAij− X

j∈Ni

ν uj− ui

(xj− xi) · nij

Aij = 0 X

j∈Ni

uij· nijAij = 0

Furthermore we will show that for these discrete incompressible Navier-Stokes equations the kinetic energy decreases too.

3.1 Discretization of the Navier-Stokes equation

In this section, we will approximate the Navier-Stokes equations using the finite volume method. Furthermore we need a computational grid. We will use a body conforming unstructured grid. We will use the co-located cell-centered grid arrangement.

The conservation of mass in cell i is discretized by taking the integral of equation 5 over the control volume Ωi:

Z

Ωi

∇ · u dΩi= Z

Γi

u · n dΓi= X

j∈Ni

uij· nijAij = 0 (13)

(8)

Figure 1: A part of an unstructured grid

Where N_iconsist of indices j of cells adjacent to cell i, A_ij the surface between cells i and j and n_ijthe normal vector to this surface A_ijin the direction of cell j as can be seen in figure 3.1. Furthermore Γ_i is the boundary of Ω_i and u_ij is the face velocity which is computed as

uij= γijui+ (1 − γij)uj

where ui is the velocity vector in cell i and γij is the geometric interpolation factor, defined as

γ_ij= (x_j− x_ij) · n_ij (xj− xi) · nij

(14) Integrating the momentum equation 2 over an arbitrary control volume Ωigives

Z

Ω_i

∂u

∂t dΩ_i+ Z

Ω_i

∇ · (uu) dΩi+ Z

Ω_i

∇p dΩi− Z

Ω_i

∇ · (ν∇u) dΩi= 0 (15) Now we will discretize the four parts of this equation. When we discretize the first part, we get

Z

Ωi

∂u

∂t dΩi≈∂ui

∂t Ωi (16)

For the second part, we use first Gauss’s divergence theorem and get Z

Ω_i

∇ · (uu) dΩi= Z

Γ_i

u(u · n) dΓ_i≈ X

j∈N_i

u_iju_ij· nijA_ij (17)

Before discretizing the third part, we again first apply Gauss’s divergence theorem.

Z

Ω_i

∇p dΩ_i= Z

Γ_i

pn dΓ_i≈ X

j∈N_i

p_ijn_ijA_ij (18) The last part is a little bit harder to discretize. When we start fromR

Ωi∇ · (ν∇u) dΩi, we get, using Gauss’s divergence theorem,

Z

Ωi

∇ · (ν∇u) dΩi= Z

Γi

ν∇u · n dΓi≈ X

j∈N_i

ν(∇u)ij· nijAij (19)

(9)

When we now assume that the fluxes are normal to the faces, we get (∇u)_ij·n_ij =

∂u

∂n

ij ≈_(x^u^j^−uⁱ

j−xi)·n_ij, so equation 19 becomes Z

Ω_i

∇ · (ν∇u) dΩi≈ X

j∈N_i

ν(∇u)_ij· nijA_ij≈ X

i∈N_i

ν u_j− u_i (xj− xi) · nij

A_ij (20)

Combining equations 15, 16, 17, 18 and 20, we get the following discrete form of the conservation of momentum equation

∂ui

∂t Ω_i+X

j∈N_i

u_iju_ij· n_ijA_ij+X

j∈N_i

p_ijn_ijA_ij−X

j∈N_i

ν uj− ui

(x_j− xi) · n_ijA_ij= 0 (21)

3.2 The evolution of the kinetic energy

In this section we will compute evolution of the kinetic energy in time for the discretized incompressible Navier-Stokes equation and show that the kinetic energy decreases.

The evolution of the kinetic energy is given by X

i

u_i·∂u_i

∂t Ω_i Using equation 21, this equation becomes

X

i

ui· ∂ui

∂t Ωi = −X

i

ui· X

j∈N_i

uijuij· nijAij−X

i

ui· X

j∈N_i

pijnijAij

+X

i

u_i· X

j∈N_i

ν uj− ui

(x_j− xi) · n_ijA_ij (22) We will now rewrite this equation using the following identity [1]

φi

X

j∈Ni

ψijQij+ ψi

X

j∈Ni

φijQij= X

j∈Ni

φψc_ijQij+ φiψi

X

j∈Ni

Qij (23)

Where φ and ψ are general variables, Qij is known on the cell face, φij = γ_ijφ_i+ (1 − γ_ij)φ_j, ψ_ij = γ_ijψ_i+ (1 − γ_ij)ψ_j, with γ_ij the interpolation function defined in equation 14, and cφψ_ij = ¹₂(φ_iψ_j+ ψ_iφ_j)

Now we take in equation 23 φ = u, ψ = φ and Q_ij = u_ij · n_ijA_ij and interpolate the face velocity as u_ij =¹₂(u_i+ u_j). Then we get

X

i

ui· X

j∈N_i

φijuij· nijAij+X

i

φi· X

j∈N_i

1

2(ui+ uj)uij· nijAij

=X

i

X

j∈N_i

1

2(ui· φj+ φi· uj)uij· nijAij+X

i

uiφi

X

j∈N_i

uij· nijAij (24)

Using equation 13, the last term vanishes. When we split the second term into two partsP

i 1

2φi· uiP

j∈N_iuij· nijAij+P

i

P

j∈N_i 1

2(φi· uj)uij· nijAij we see, using again equation 13, that the first part vanishes and the second part is also on the other hand of the equal sign in equation 24, so equation 24 becomes

X

i

ui· X

j∈N_i

φijuij· nijAij =X

i

X

j∈N_i

1

2(ui· φj)uij· nijAij (25)

(10)

When we now substitute φ = u, ψ = p and Q_ij = n_ijA_ij into equation 23, we get

X

i

u_i· X

j∈N_i

p_ijn_ijA_ij+X

i

p_i X

j∈N_i

u_ij· n_ijA_ij

=X

i

X

j∈N_i

cup_ijn_ijA_ij+X

i

u_ip_i· X

j∈N_i

n_ijA_ij

Using equation 13, we see that the second term vanishes. Using geometric arguments, we see that P

j∈Nin_ijA_ij = 0, so the last term vanishes also. Thus we end up with

X

i

u_i· X

j∈N_i

p_ijn_ijA_ij =X

i

X

j∈N_i

upc_ijn_ijA_ij (26) Putting equations 22, 25 (with φ = u) and 26 together, we get the following equation for the kinetic energy

X

i

ui·∂ui

∂t Ωi= −X

i

X

j∈Ni

1

2(ui· uj)uij· nijAij−X

i

X

j∈Ni

upc_ijnijAij

+ νX

i

ui· X

j∈N_i

uj− ui

(x_j− xi) · n_ijAij (27) So this equation is the discrete version of equation 12. In this equation, we see that the first two terms on the right hand side vanishes at the interior. This is because uij = uji, nij = −nji, Aij = Aji and upc_ij =upc_ji. The only relation here which is not trivial is uij = uji, so we will prove this relation. Therefore we will first show that γij = 1 − γji.

1 − γji=(xi− xj) · nji

(x_i− x_j) · n_ji −(xi− xji) · nji

(x_i− x_j) · n_ji =(xj− xij) · nij

(x_j− x_i) · n_ij = γij

With this, we see that

uij = γijui+ (1 − γij)uj= (1 − γji)ui+ γjiuj = uji

So when we use the no-slip boundary condition, the first two terms of equation 27 vanishes and we end up with

X

i

ui· ∂u_i

∂t Ωi = νX

i

ui· X

j∈Ni

u_j− ui

(xj− xi) · nij

Aij

To see why this is non-positive, we use the fact that if j ∈ Ni, then i ∈ Nj, so for a term νui·_(x^u^j^−uⁱ

j−xi)·nijAij with j ∈ Ni, we get also a term νuj·_(x^uⁱ^−u^j

i−xj)·njiAji

with i ∈ Nj. So we get the sum over pairs X

i

ui·∂ui

∂t Ωi= νX

i

X

j∈Ni

ui· uj− ui

(xj− xi) · nij

Aij

= 1 2νX

i

X

j∈Ni

ui· uj− ui

(xj− xi) · nij

Aij+ uj· ui− uj

(xi− xj) · nji

Aji

(11)

When we now use u_j·_(x^uⁱ^−u^j

i−xj)·njiA_ji= −u_j· _(x^u^j^−uⁱ

j−xi)·nijA_ij we get X

i

u_i·∂u_i

∂t Ω_i=1 2νX

i

X

j∈N_i

(u_i− u_j) · u_j− u_i (xj− xi) · nij

A_ij

=1 2 ν

|{z}

>0

X

i

X

j∈N_i

(u_i− uj) · (u_j− ui)

| {z }

≤0

1 (xj− xi) · nij

| {z }

>0

A_ij

|{z}

>0

≤ 0

From this we see that the kinetic energy equals zero if ui= ujfor all pairs i and j ∈ Ni, which means that u is constant. Because we use the no-slip boundary condition, u = 0 in this case. So if u 6= 0, the kinetic energy will be negative. So the discrete and continuous form of the incompressible Navier-Stokes equations give the same conclusion.

4 QR model

In this section we will derive a model called the QR model. Therefore we will first give some notation and useful lemmas which we will use in the derivation of the OR model. Then we will start from large-eddy simulation (LES) and use scale-truncation to get the QR model

τ −1

3tr(τ )I = −2Cδ

max {r(v), 0}

q(v) S

4.1 notation and useful lemmas

In this subsection, we will give some notations and lemmas which we will use during this section.

4.1.1 notation

For clarity, we will first recall the definition of the Frobenius inner product A : B =X

i,j

AijBij = tr(A^TB) = tr(B^TA)

We will use the following notations

|A|²:= A : A

|A|³:= A²: A

||A||²:=

Z

Ωδ

A : A dx = Z

Ωδ

|A|² dx

||A||³:=

Z

Ωδ

A²: A dx = Z

Ωδ

|A|³ dx

(∇v)_ij:= ∂vj

∂x_i

(12)

−

→e_i is the i^th standard basic vector in R³. We define the Levi-Civita symbol _ijk as

_ijk:=







1 if (i, j, k) is (1, 2, 3), (2, 3, 1) or (3, 1, 2)

−1 if (i, j, k) is (1, 3, 2), (2, 1, 3) or (3, 2, 1) 0 if i = j, i = k or j = k

(28)

For the rest of this section, we will call S and T the symmetric and antisymmetric part of ∇v:

S := 1

2 ∇v + ∇v^T , T := 1

2 ∇v − ∇v^T We define ω as the curl of v.

The QR model uses the invariants of S. We define q(v) and r(v) as minus the second and third invariant of S. Then q(v) and r(v) are given by

q(v) = 1

2|S|², r(v) = −1 3|S|³ We will show this in lemma 4.1.

In this section we will not write down the summation sign (P) for clarity reasons. So one should read for example P

ijAijBij for AijBij. When we integrate over Ωδ, we use periodic boundary conditions, so the boundary terms are equal to zero. So every time we use integration by parts, we will neglect the boundary terms.

4.1.2 useful lemmas

Lemma 4.1. Let v be in R³such that ∇ · v = 0 and let S be the symmetric part of ∇v, then the second and third invariant of S, q(v) and r(v) are given by

q(v) = 1

2|S|², r(v) = −1 3|S|³

Proof. Let v ∈ R³ be such that ∇ · v = 0 and let S be the symmetric part of

∇v. Then the invariants of S are defined using the characteristic polynomial of S:

λ³+ Iλ²+ IIλ + III and are given by

I = tr(S) II =1

2 tr(S)²− tr(S²) III = det(S) = 1

6tr(S)³−1

2tr(S²)tr(S) + 1 3tr(S³) When we now use ∇ · v = 0, these invariants become

I = tr(S) = ∇ · v = 0 II = 1

2 tr(S)²− tr(S²) = −1 2tr(S²) III = 1

6tr(S)³−1

2tr(S²)tr(S) +1

3tr(S³) = 1 3tr(S³)

(13)

When we now call q(v) minus the second invariant and use the definition of the Frobenius inner product, we get

q(v) = 1

2tr(S²) = 1

2tr(S^TS) =1 2|S|²

When we now call r(v) minus the third invariant and use the definition of the Frobenius inner product, we get

r(v) = −1

3tr(S³) = −1

3tr(S^TS²) = −1 3|S|³

Lemma 4.2. If A and B are respectively a symmetric and a skew-symmetric matrix, then A : B = 0

Proof. Let A be a symmetric matrix and B an antisymmetric matrix. Then A : B = AijBij = AjiBji= A^T

ij B^T

ij

= A^T : B^T = A : (−B) = −A : B So A : B = −A : B, thus A : B has to be zero

Corollary 4.3. If S is a symmetric matrix and A a matrix, then S : A = S : B

where B is the symmetric part of A

Lemma 4.4. If v is a vector in R³with ∇ · v = 0, then ∆v = ∇ · ∇v = 2∇ · S, where S is the symmetric part of ∇v

Proof. Let v = (v1, v2, v3) and ∇ · v = ^∂v_∂xⁱ

i = 0 then 2∇ · S = ∇ · (∇v + ∇v^T)

= ∂

∂xi

∂v_j

∂xi

+ ∂v_i

∂xj

−→ej

= ∂²v_j

∂x²_i

−

→e_j+ ∂

∂xj

∂v_i

∂xi

−

→e_j

This last term is equal to zero because ^∂v_∂xⁱ

i = 0, so 2∇ · S =∂²vj

∂x²_i

−

→ej = ∆v

Lemma 4.5. If ω is the curl of a vector v ∈ R³, then we can write ωi as ωi= ijk

∂vk

∂xj

(14)

Proof. Let ω be the curl of v. Then, by definition,

ω =





−

→e₁ −→e₂ −→e₃

∂

∂x₁

∂

∂x₂

∂

∂x₃

v1 v2 v3





= ∂v3

∂x₂ − ∂v2

∂x₃

−→e1+ ∂v1

∂x₃− ∂v3

∂x₁

−→e2+ ∂v2

∂x₁− ∂v1

∂x₂

−→e3

= ijk

∂v_k

∂xj

−

→ei

So ωi= ijk∂v_k

∂x_j

Lemma 4.6. If v is a vector in R³, T the antisymmetric part of ∇v and ω the curl of v, then we can write Tij as

Tij= 1 2ijkωk

Proof. Let v ∈ R³ and let T and ω be respectively the antisymmetric part and the curl of v. Then we get when we write T out

T =1

2 ∇v − (∇v)^T

=1 2







0 ^∂v_∂x²

1 −^∂v_∂x¹

2

∂v₃

∂x₁ −^∂v_∂x¹

∂v₁ 3

∂x₂ −_∂x^∂v²

1 0 _∂x^∂v³

2 −^∂v_∂x²

∂v₁ 3

∂x₃ −_∂x^∂v³

1

∂v₂

∂x₃ −^∂v_∂x³

2 0







When we now apply lemma 4.5, we get

T =1 2





0 ω3 −ω2

−ω3 0 ω1

ω2 −ω1 0



= 1 2





0 123ω3 132ω2

213ω3 0 231ω1

132ω2 321ω1 0





So

Tij= 1 2ijkωk

Lemma 4.7. If a is a vector in R³ with ∇ · a = 0, then

∇ × (∇ × a) = −∆a

Proof. Let a ∈ R³ be such that ∇ · a = 0, then we use 4.5 to see that the i^th element of a is given by

(∇ × (∇ × a))_i= _ijk ∂

∂x_j(∇ × a)_k = _ijk ∂

∂x_j

_klm∂am

∂x_l

Now we use _klm= _lmk to get

(∇ × (∇ × a))_i= _ijk_lmk ∂

∂x_l

∂am

∂x_j

(15)

The product _ijk_lmk is not zero if i 6= j 6= k and l 6= m 6= k, in which cases i = l and j = m or i = m and j = l. Thus

(∇ × (∇ × a))_i= _ijk_ijk ∂

∂a_i

∂vj

∂x_j + _ijk_jik ∂

∂x_j

∂ai

∂x_j

= γij

∂

∂x_i

∂aj

∂x_j − γij

∂²ai

∂x²_j where γij = 0 if i = j

1 if i 6= j. When we now use ∇ · a = _∂x^∂a^j

j = 0, the first term becomes

γij

∂

∂xi

∂a_j

∂xj

= − ∂

∂xi

∂a_i

∂xi

= −∂²a_i

∂x²_i so we end up with

(∇ × (∇ × a)) = −∂²ai

∂x²_i − γ_ij∂²ai

∂x²_j

!

−

→e_i = −∂²ai

∂x²_j

−

→e_i = −∆a

Lemma 4.8. If a and b are vectors in R³, then, if we neglect boundary terms, Z

(∇ × a) · b dx = Z

a · (∇ × b) dx Proof. Let a, b ∈ R³. Then

Z

(∇ × a) · b dx = Z

_ijk∂ak

∂x_jb_i dx Now we use integration by parts to get

Z

(∇ × a) · b dx = − Z

ijkak

∂bi

∂xj

dx

where we neglect the boundary terms. Now we use ijk = −kji to get Z

(∇ × a) · b dx = Z

a_k_kji∂b_i

∂xj

dx = Z

a_k(∇ × b)_k dx = Z

a · (∇ × b) dx

Corollary 4.9. When we combine lemmas 4.7 and 4.8, we get Z

(∇ × a) · (∇ × b) dx = Z

a · (∇ × (∇ × b)) dx − Z

a · ∆b dx Lemma 4.10. If a is a vector in R³ with ∇ · a = 0, then

∇ · (∆a) = 0 Proof. Let a ∈ R³such that ∇ · a = _∂x^∂aⁱ

i = 0, then

∇ · (∆a) = ∂

∂xi

∂²ai

∂x²_j = ∂²

∂x²_j

∂ai

∂xi

= 0

(16)

Lemma 4.11. If a is a vector in R³, S the symmetric part of ∇a and ˜S the symmetric part of ∇∆a, then

S = ∆S˜

Proof. Let a ∈ R³, S be the symmetric part of ∇a and ˜S be the symmetric part of ∇∆a. Then we get when we write ˜Sij =¹₂ (∇∆a)ij+ ((∇∆a)^T)ij out

S˜_ij = 1 2

∂(∆a)_j

∂xi

+∂(∆a)_i

∂xj

=1 2

∂

∂xi

∂²a_j

∂x²_k + ∂

∂xj

∂²a_i

∂x²_k

When we now reorder terms we get S˜ij = ∂²

∂x²_k

1 2

∂aj

∂x_i +∂aj

∂x_i

= ∂²

∂x²_kSij= (∆S)ij

Lemma 4.12. If v is a vector in R³ with ∇ · v = 0 and we use the definition of S stated at the start of this subsection, then

Z

(v · ∇)v · ∆v dx = − Z

∇v^T∇v : S dx

Proof. Let v ∈ R³ be such that ∇ · v = 0 and let S be the symmetric part of

∇v.

Writing outR (v · ∇)v · ∆v dx gives Z

(v · ∇)v · ∆v dx = Z

vk

∂

∂xk

vj

∂²

∂x²_ivj dx (29) Integration by parts gives

Z vk

∂

∂xk

vj

∂²

∂x²_ivj dx = − Z ∂

∂xi

vk

∂vj

∂xk

∂vj

∂xi

dx

= − Z ∂v_k

∂xi

∂v_j

∂xk

∂v_j

∂xi

dx − Z

vk

∂²v_j

∂xi∂xk

∂v_j

∂xi

dx (30) We will look a little bit closer on this last term. For clarity, we call a = ^∂v_∂x^j

i, so we get

Z vk

∂a

∂x_ka dx =

Z ∂

∂x_k (vka) a dx − Z ∂vk

∂x_ka² dx Using ∇ · v = ^∂v_∂x^k

k = 0 the last term vanishes, so Z

v_k ∂a

∂xk

a dx =

Z ∂

∂xk

(v_ka) a dx When we now use integration by parts we get

Z v_k ∂a

∂xk

a dx = − Z

v_ka ∂a

∂xk

dx = − Z

v_k ∂a

∂xk

a dx So we conclude

Z vk

∂a

∂xk

a dx = 0

(17)

Now, as we go back to equation 30, we get Z

v_k ∂

∂xk

v_j ∂²

∂x²_iv_j dx = − Z ∂v_k

∂xi

∂v_j

∂xk

∂v_j

∂xi

dx (31)

When we work this out, we get

∂vk

∂xi

∂vj

∂xk

∂vj

∂xi

= (∇v)ik(∇v)kj(∇v)ij = (∇v∇v)ij(∇v)ij = (∇v²) : ∇v With the definition of the Frobenius inner product, we rewrite this as

∂vk

∂x_i

∂vj

∂x_k

∂vj

∂x_i = (∇v²) : ∇v = tr((∇v²)^T∇v) = tr(∇v^T(∇v^T∇v)) = (∇v^T∇v) : ∇v So equation 31 becomes, using corollary 4.3,

Z v_k ∂

∂xk

v_j ∂²

∂x²_iv_j dx = Z

(∇v^T∇v) : ∇v dx = Z

(∇v^T∇v) : S dx Thus we get, using equation 31,

Z

(v · ∇)v · ∆v dx = − Z

∇v^T∇v : S dx

Lemma 4.13. If v is a vector in R³ with ∇ · v = 0 and we use the definitions of ω and ijk stated at the start of this subsection, then

∇ × ((v · ∇)v) = −(ω · ∇)v + (v · ∇)ω

Proof. Let v ∈ R³be such that ∇ · v = 0 and let ω be the curl of v and ijk be the Levi-Civita symbol.

Using lemma 4.5, we get

∇ × ((v · ∇)v) = ∇ ×

v_i∂vj

∂x_i

−

→e_j

= _jkl ∂

∂x_k

v_i∂vl

∂x_i

−→e_j

= jkl

∂vi

∂xk

∂vl

∂xi

−

→ej+ jklvi

∂²vl

∂xi∂xk

−

→ej (32)

We will work this two terms out. For jkl∂vi

∂x_k

∂vl

∂xi

−

→ej, the j^th component is given by jkl

∂v_i

∂x_k

∂v_l

∂x_i −^∂v_∂xⁱ

l

∂v_k

∂x_i

, where k and l are chosen such that k 6= l 6= j. Now we write out the sum over i by filling in j, k and l for i and add_∂x^∂v^j

k

∂vj

∂x_l−^∂v_∂x^j

l

∂vj

∂x_k.

jkl

∂v_i

∂xk

∂v_l

∂xi

−∂v_i

∂xl

∂v_k

∂xi

= jkl

∂v_l

∂xj

∂v_j

∂xk

−∂v_k

∂xj

∂v_j

∂xl

+ ∂v_l

∂xk

∂v_k

∂xk

− ∂v_k

∂xk

∂v_k

∂xl

+∂vl

∂xl

∂vl

∂xk

−∂vk

∂xl

∂vl

∂xl

+ ∂vj

∂xk

∂vj

∂xl

−∂vj

∂xl

∂vj

∂xk

= −jkl

∂v_l

∂xk

−∂v_k

∂xl

−∂v_k

∂xk

−∂v_l

∂xl

+ ∂v_j

∂xl

− ∂v_l

∂xj

∂v_j

∂xk

+ ∂v_k

∂xj

− ∂v_j

∂xk

∂v_j

∂xl

(18)

Now we use ^∂v_∂x^j

j = −_∂x^∂v^k

k −^∂v_∂x^l

l to get

_jkl ∂v_i

∂xk

∂v_l

∂xi

− ∂v_i

∂xl

∂v_k

∂xi

= −_jkl ∂v_l

∂xk

−∂v_k

∂xl

∂v_j

∂xj

+ ∂v_j

∂xl

− ∂v_l

∂xj

∂v_j

∂xk

+ ∂v_k

∂xj

− ∂v_j

∂xk

∂v_j

∂xl

When we now use lemma 4.5, we get

_jkl ∂v_i

∂xk

∂v_l

∂xi

−∂v_i

∂xl

∂v_k

∂xi

= −

ω_j∂v_j

∂xj

+ ω_k∂v_j

∂xk

+ ω_l∂v_j

∂xl

= −

ω_i ∂

∂xi

v_j Thus

jkl

∂vi

∂xk

∂vl

∂xi

−

→ej = −

ωi

∂

∂xi

vj−→ej = −(ω · ∇)v (33) For jklvi ∂²vl

∂x_i∂x_k

−

→ej, we get, using again lemma 4.5,

_jklv_i ∂²v_l

∂xi∂xk

−

→e_j = v_i ∂

∂xi

_jkl∂v_l

∂xk

−

→e_j = (v · ∇)ω (34)

Combining equations 32, 33 and 34 gives us

∇ × ((v · ∇)v) = −(ω · ∇)v + (v · ∇)ω

Lemma 4.14. If v is a vector in R³ with ∇ · v = 0 and we use the definition of ω stated at the start of this subsection, then

Z

((v · ∇)ω) · ω dx = 0

Proof. Let v ∈ R³be such that ∇ · v = 0 and let ω be the curl of v.

Writing outR ((v · ∇)ω) · ω dx gives Z

((v · ∇)ω) · ω dx = Z

vi

∂ω_j

∂xi

ωj dx = Z ∂ω_j

∂xi

viωj dx Now we use integration by parts

Z

((v · ∇)ω) · ω dx = − Z

ωj

∂viωj

∂x_i dx

= − Z

ω_j∂v_i

∂xi

ω_j dx − Z

ω_jv_i∂ω_j

∂xi

dx Now we use ∇ · v = ^∂v_∂xⁱ

i = 0 and see that the first term is equal to zero. So we get

Z

((v · ∇)ω) · ω dx = − Z

vi

∂ωj

∂x_iωj dx

= − Z

((v · ∇)ω) · ω dx

(19)

So we conclude

Z

((v · ∇)ω) · ω dx = 0

Lemma 4.15. If v is a vector in R³ with ∇ · v = 0 and we use the definitions of ω, _ijk, S and T stated at the start of this subsection, then

((ω · ∇)v) · ω = S_ijω_iω_j

Proof. Let v ∈ R³ be such that ∇ · v = 0 and use the definitions of ω, ijk, S and T stated at the start of this subsection. Then writing out ((ω · ∇)v) · ω gives

((ω · ∇)v) · ω =

ωi

∂

∂xi

vj

ωj = ωi

∂vj

∂xi

ωj

When we now use ^∂v_∂x^j

i = (∇v)_ij= S_ij+ T_ij we get

((ω · ∇)v) · ω = S_ijω_iω_j+ T_ijω_iω_j (35) Let ω123= ω1ω2ω3, then, using lemma 4.6 andP

ijkijk= 0 we get Tijωiωj =1

2ijkωkωiωj= 1

2ijkω123= 0 When we put this into equation 35, we get

((ω · ∇)v) · ω = Sijωiωj

Lemma 4.16. If v is a vector in R³ with ∇ · v = 0 and we use the definitions of ω, _ijk, S and T stated at the start of this subsection, then

(∇v^T∇v) : S = |S|³−1

4Sjkωjωk

Proof. Let v ∈ R³ be such that ∇ · v = 0 and use the definitions of ω, ijk, S and T stated at the start of this subsection. We rewrite (∇v^T∇v) : S using

∇v = S + T

(∇v^T∇v) : S = (S + T )^T(S + T ) : S

= S^TS + S^TT + T^TS + T^TT : S

Using the definition of the Frobenius inner product, we get

S^TT : S + T^TS : S = tr(T S^TS) + tr(T^TSS^T) = tr(T SS) − tr(T SS) = 0 So

(∇v^T∇v) : S = S²: S + T^TT : S = |S|³+ TikTijSkj (36)

(20)

We compute S_jkT_ikT_ij using lemma 4.6:

SjkTikTij = Sjk

1

2ikmωm

1 2ijlωl

= Sjk

1

4ikmijlωmωl

For the product ikmijl we have two possibilities, namely

_ikm_ijl= 1 if k = j and l = m

−1 if l = k and m = j So we get

S_jkT_ikT_ij= 1

4S_jjω_mω_m−1

4S_jkω_jω_k Using S_jj = ^∂v_∂x^j

j = 0 this equation becomes SjkTikTij= −1

4Sjkωjωk

When we fill this into equation 36, we get (∇v^T∇v) : S = |S|³−1

4S_jkω_jω_k

4.2 LES

In this subsection, we introduce a methodology called large-eddy simulation (LES). For LES, we add a filter to the incompressible Navier-Stokes equations.

In this filter, only the large eddies in a flow are resolved. So we get a filtered solution ¯u which is much easier to compute by a computer.

We start again with the Navier-Stokes equations (equations 2 and 5)

∂u

∂t + ∇ · (u ⊗ u) + ∇p − ν∆u = 0

∇ · u = 0

Now we will use a filter, which we will denote by an overline. Because it will be a linear filter, we can apply it on each component:

∂u

∂t + ∇ · (u ⊗ u) + ∇p − ν∆u = 0 = 0 When we now assume that ∇ · a = ∇ · a, we get

∂u

∂t + ∇ · u ⊗ u + ∇p − ν∆u = 0

When we rewrite this a little bit, we get the filtered incompressible Navier-Stokes equations:

∂ ¯u

∂t + ∇ · (¯u ⊗ ¯u) + ∇¯p − ν∆¯u = −∇ · (u ⊗ u − ¯u ⊗ ¯u)

(21)

∇ · ¯u = 0

with filtered velocity field ¯u. When we now call τ = u ⊗ u − ¯u ⊗ ¯u and use the approximation v ≈ ¯u and p ≈ ¯p, we get, using equation 3

∂v

∂t + (v · ∇)v − ν∆v + ∇p = −∇ · τ (37)

∇ · v = 0 (38)

4.3 Finding the QR model

In this section we will find the QR model by splitting v into an average ¯v and residue v⁰ and find the following condition for which the residual field will be suppressed.

4 Z

Ω_δ

r(v) dx ≤ Z

Ω_δ

τ : ∆S dx We set v = ¯v + v⁰, where

¯ v = 1

|Ωδ| Z

Ωδ

v(x, t) dx (39)

is the average of v on a periodic box Ωδ with diameter δ. The residual field in v is then v⁰. When we want to use this ¯v in our model, the residual field should not become significant. Therefore we bound the kinetic energy from above

1 2

d

dt||v⁰||²≤ 0 (40)

Note that when we satisfy this condition, and we add the initial condition

||v⁰||(t = 0) = 0 and boundary condition v⁰ = 0, we get ||v⁰|| = 0 for all t > 0. In order to satisfy condition 41, we use the Poincar´e inequality [3]

1

2||v⁰||²≤ 1

2Cδ||∇v||²

where C_δ > 0 is the Poincar´e constant. So when we satisfy the condition 1

2 d

dt||∇v||²≤ 0 (41)

we also satisfy our condition 40:

1 2

d

dt||v⁰||²≤ Cδ

1 2

d

dt||∇v||²≤ 0

With this Poincar´e condition, we do not have to compute v⁰ explicitly.

Now we will compute ¹₂_dt^d||∇v||² and see that it is given by 1

2 d

dt||∇v||²= −4

3||S||³− 2ν||∇S||²− Z

Ω_δ

τ : ∆S dx

When we compute ¹₂_dt^d||∇v||², we get, using integration by parts, 1

2 d

dt||∇v||²= Z

Ωδ

∇∂v

∂t · ∇v

dx =

Z

Ωδ

−∂v

∂t · ∆v

dx

(22)

Now we fill in equation 37 and get 1

2 d

dt||∇v||²= Z

Ω_δ

(v · ∇)v · ∆v dx − Z

Ω_δ

ν∆v · ∆v dx +

Z

Ω_δ

∇p · ∆v dx + Z

Ω_δ

(∇ · τ ) · ∆v dx (42) We will first rewrite the first part (R

Ω_δ(v · ∇)v · ∆v dx) of equation 42 and see that it is equal to −⁴₃||S||³.

We will start again with our filtered Navier-Stokes equations 37 and 38

∂v

∂t + (v · ∇)v − ν∆v + ∇p = −∇ · τ (43)

∇ · v = 0 (44)

From this equations, we will compute ¹₂_dt^d||∇v||² and ¹₂_dt^d||ω||²and see that these two quantities are the same and from this will follow that R

Ω_δ(v · ∇)v ·

∆v dx = −⁴₃||S||³

Computing ¹₂_dt^d||∇v||² gives again equation 42:

1 2

d

dt||∇v||²= Z

Ω_δ

(v · ∇)v · ∆vdx + Z

Ω_δ

∇p · ∆v dx

− ν||∆v||²+ Z

Ωδ

(∇ · τ ) · ∆v dx

Using equation 53R

Ωδ∇p · ∆v dx = 0

, this becomes 1

2 d

dt||∇v||²= Z

Ω_δ

(v · ∇)v · ∆vdx − ν||∆v||²+ Z

Ω_δ

(∇ · τ ) · ∆v dx (45) Now we take the curl of equation 43 and use ∇ × ∇p = 0 and get

∂ω

∂t + ∇ × ((v · ∇)v) − ν∇ × ∆v = −∇ × (∇ · τ ) Using lemma 4.13 this equation becomes

∂ω

∂t − (ω · ∇)v + (v · ∇)ω − ν∇ × ∆v = −∇ × (∇ · τ ) We will now compute ¹₂_dt^d||ω||²using the above equation

1 2

d

dt||ω||²= Z

Ωδ

∂ω

∂t · ω dx

= Z

Ω_δ

− ((v · ∇)ω) · ω + ((ω · ∇)v) · ω + ν(∇ × ∆v) · ω − (∇ × (∇ · τ )) · ω dx Using lemma 4.14 this equation becomes

1 2

d

dt||ω||²= Z

Ω_δ

((ω · ∇)v) · ω + ν(∇ × ∆v) · ω − (∇ × (∇ · τ )) · ω dx (46)

2 The incompressible Navier-Stokes equations

Modeling turbulent flow

Bachelor Project Applied Mathematics

Contents

1 Introduction

2 The incompressible Navier-Stokes equations

2.1 The incompressible Navier-Stokes equations

2.2 The evolution of the kinetic energy

3 Discretization of the Navier-Stokes equation

3.1 Discretization of the Navier-Stokes equation

3.2 The evolution of the kinetic energy

4 QR model

4.1 notation and useful lemmas

4.2 LES

4.3 Finding the QR model