Answers to Exercises Computational Fluid Dynamics

(1)

Answers to

Exercises Computational Fluid Dynamics

Exercise 1 - Artificial diffusion

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

k=0.1000

upwind computations

x k=0.0100

k=0.0010 k=0.0010 exact

For k = 0.1, 0.01 and 0.001, and for N = 20, the discrete upwind solutions outside the boundary layer approach, the analytic solution for k = 0.025. Explanation: the effective viscosity of the upwind solution is given by

k_{ef f} = k + uh

2 = k + 1 40, and this approaches 0.025 when k → 0.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

k=0.1000

upwind computations

x k=0.0100

k=0.0010 k=0.0250 exact

(2)

Exercise 2 - Various discretization methods

0 0.5 1

upwind

0 0.5 1

smart upwind

0 0.5 1

B3

0 0.5 1

QUICK

0 0.5 1

central

0 0.5 1

upwind nonuniform

0 0.5 1

central nonuniform A

0 0.5 1

central nonuniform B N=20

k=0.01 discrete exact

Solutions for k = 0.01

0 0.5 1

upwind

0 0.5 1

smart upwind

0 0.5 1

B3

0 0.5 1

QUICK

0 0.5 1

central

0 0.5 1

upwind nonuniform

0 0.5 1

Solutions for k = 0.001

(3)

Question 4: Stretched grid

With Method B, the first time one of the eigenvalues of the coefficient matrix crosses the imaginary axis is for k ≈ 0.00425. For a slightly smaller value of k, viz. k ≈ 0.00342576, an eigenvalue passes through zero (within MATLAB accuracy) and the matrix becomes singular.

Question 5: Monotonicity of λ-schemes

The λ-schemes applied to our convection-diffusion equation yield a difference scheme given by λ

hφ_i−2− (k

h² + 6λ + 1

2h ) φ_i−1+ (2k h² +3λ

h ) φ_i− ( k

h² +2λ − 1

2h ) φ_i+1= RHS.

This scheme is of the type

a φ_i−2+ b φ_i−1+ c φ_i+ d φ_i+1= RHS.

The characteristic roots r of this equation satisfy d r³ + c r² + b r + a = 0. Since r = 1 is a solution of this equation (why?), the factor r − 1 can be divided. There remains

d r²+ (c + d) r − a = 0,

where the coefficient of r, given by c + d ≡ _h^k2 + ^4λ+1_2h , is always positive. The product of both roots is given by −a/d, whereas their sum is given by −(c + d)/d. With a always positive, it is observed that for d < 0 both the sum and the product of the roots are positive, hence both roots are positive. But for d > 0 the product is negative, hence one of the roots is negative.

In the latter situation (d > 0) with a negative root, oscillations are possible. More precise, when applied to the difference equation for the λ-scheme, oscillations can occur when

−d ≡ k

h² + 2λ − 1 2h < 0.

• For the B3 scheme, with λ = 1/2, this can never happen.

• For the QUICK scheme, with λ = 1/8, the oscillation boundary P = 8/3 results, where P is the mesh-P´eclet number. By determining the discrete solution for a number of values of k close to the boundary (for N = 20 it is given by k = 0.01875), it is observed that also in practice this boundary separates the monotone QUICK solutions from the oscillating ones.

Remark Note that the coefficient of φ_i−2is always positive. Hence the discretization will never lead to a positive operator, and the theory from Appendix A.1.3 is not applicable.

Question 6: One interior grid point

With one interior grid point located at x = 1 − 2k the discretization for Methods A becomes (φ₊− φ−) − k(1 − 2k)φ+− φ0+ 2kφ₋

k(1 − 2k) = 0.

and that of Method B

(1 − 2k)²φ₊+ ((2k)²− (1 − 2k)²)φ₀− (2k)²φ₋

2k(1 − 2k) − k(1 − 2k)φ+− φ0+ 2kφ₋ k(1 − 2k) = 0.

(4)

0 0.5 1 0

0.5 1

upwind

0 0.5 1

smart upwind

0 0.5 1

B3

0 0.5 1

QUICK

0 0.5 1

central

0 0.5 1

upwind nonuniform

0 0.5 1

Solutions for k = 0.00342576 With φ₋ = 0 and φ₊= 1, the solution in the interior

grid point becomes

Method A: φ₀ = 0, and

Method B: φ₀ = (1 − 2k)(1 − 4k) 1 − 6k .

For small k the latter value is close to 1. The analytical solution in x = 1 − 2k is given by

φ0 = 1 − e^(1−2k)/k 1 − e^1/k , which tends to e⁻² ≈ 0.135 for small k.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.2 0.4 0.6 0.8 1

exact Lagrange symmetric

Remark The studied problem is a slightly simplified version of the problem with which the QUICK method was introduced (B.P. Leonard: A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comp. Meth. Appl. Mech. Engng. 19, pp. 59–98, 1979). In this paper the B3 method was not entered as a comparison. Many years later Leonard did make this comparison (Int. J. Num. Meth. Engng. 30, pp. 729–766, 1990) and found, just like we did, that B3 is a serious competitor for QUICK.

(5)

Exercise 3 - JOR and SOR

1. Theoretically, Gauss–Seidel converges twice as fast as Jacobi.

uniform Jacobi method Gauss–Seidel method grid convergent? # iterations convergent? # iterations

upwind yes 26 yes 10

central no no

2. Empirically obtained optimum and maximum values for the relaxation factor ω:

uniform JOR method SOR method

grid ωopt # it. ωmax # it. ωopt # it. ωmax # it.

upwind 1.000 26 1.31 11411 1.09 6 1.999 15499

central 0.0485 308 0.0880 > 50000 0.3452 26 0.3534 > 20000 The relevant Jacobi eigenvalues are 0.52572 (upwind) and 4.6592i (central).

3. Stretched grids: Jacobi versus Gauss–Seidel.

stretched Jacobi method Gauss-Seidel method grid convergent? # iterations convergent? # iterations

upwind yes 83 yes 41

central A no no

central B no no

Again, when converging, Gauss–Seidel is twice as fast as Jacobi (upwind eigenvalue 0.8918).

The central methods are in trouble; observe the position of the eigenvalues with respect to the unit circle. Central A: 7.6843i; Central B: 7.6423i and 1.0015.

4. Stretched grids with tuned relaxation factors:

stretched JOR method SOR method

grid ωopt # it. ωmax # it. ωopt # it. ωmax # it.

upwind 0.983 79 1.057 23108 1.382 15 1.999 12256

central A 0.0247 1173 0.0333 > 100000 0.2294 153 0.2303 > 100000

central B 0.00707 1719^† – – 0.00703 1722^† – –

†) The iterations reach the convergence criterion, but will eventually diverge.

Note that the upwind method and the central method A can be made convergent with suitably chosen relaxation factors. Central method B cannot be made convergent, because of the Jacobi eigenvalue 1.0015.

5. The upwind matrix reads

diag

−u h − k

h², u h +2k

h², − k h²

, or with some scaling

diag[−P − 1, P + 2, −1],

where P is the mesh-P´eclet number. The Jacobi matrix then becomes diag

−P + 1

P + 2, 0, − 1 P + 2

.

(6)

Its eigenvalues µ_J,upw are given by µ_J,upw = 2

P + 2

√P + 1 cos(iπ

N), i = 1, . . . , N − 1, (1) where N is the number of grid cells. These are obviously real for all P ≥ 0 and satisfy

−1 ≤ µJ ≤ 1.

As already given, for the central discretization the Jacobi eigenvalues are given by µ_J,cent = ¹₂^p4 − P² cos(iπ

N), i = 1, . . . , N − 1, (2) 6. Let the eigenvalues of the Jacobi matrix be either real µ_J = µ_r or purely imaginary

µ_J = iµ_I.

– In case of real eigenvalues the optimum and maximum relaxation factors for JOR and SOR are given by

JOR: ωJOR,opt= 1, ωJOR,max= 2

1 + µ_R, (3)

(derive these yourself), and

SOR: ω_SOR,opt= 2 1 +^q1 − µ²R

, ω_SOR,max= 2 (4)

(see Section A.3.2).

– For purely imaginary Jacobi eigenvalues, the optimum and maximum relaxation factors are given by

JOR: ω_JOR,opt= 1

1 + µ²_I, ω_JOR,max = 2

1 + µ²_I, (5)

and

SOR: ωSOR,opt= 2 1 +^q1 + µ²_I

, ωSOR,max = 2

1 + µ_I, (6)

respectively.

7. When k = 0.01 and N = 10, the mesh P´eclet number is given by P = 10. On a uniform grid the Jacobi eigenvalues can be predicted theoretically. From (1) we find that the relevant Jacobi eigenvalue for upwind discretization becomes µJ,upw = 0.5257. Similarly, from (2), for central discretization the relevant Jacbi eigenvalue is µ_J,cent = 4.6592. Now the optimum and maximum relaxation factors can be estimated with the formulas from Question 6:

uniform equations JOR method SOR method grid for ω ωopt ωmax ωopt ωmax

upwind (3), (4) 1.000 1.3109 1.0807 2.0 central (5), (6) 0.0440 0.0880 0.3469 0.3534 Notice the excellent agreement with the findings from Question 2.

(7)

8. On a stretched grid, the Jacobi eigenvalues cannot be predicted theoretically. Instead we have to use the computed Jacobi eigenvalues found in Question 3 to produce estimates for optimum and maximum relaxation factors. As in the previous question, the optimum and maximum relaxation factors for JOR and SOR can be obtained from the above formulas.

Note that for the Central A the numbers below are based on the largest imaginary eigenvalue, although also real eigenvalues exist which may interfere. Central B never converges because of the Jacobi eigenvalue 1.0015.

stretched equations JOR method SOR method grid for ω ω_opt ω_max ω_opt ω_max upwind (3), (4) 1.000 1.0572 1.3770 2.0 central A (5), (6) 0.0167 0.0334 0.2286 0.2303

central B – – – –

(8)

Exercise 4 - Time integration

The discretization with the generalized Crank-Nicolson method with parameter ω reads:

T_iⁿ⁺¹− Tiⁿ− ωd

2(T_i+1ⁿ⁺¹− 2Tiⁿ⁺¹+ T_i−1ⁿ⁺¹) − (1 − ω)d

2(T_i+1ⁿ − 2Tiⁿ+ T_i−1ⁿ ) = 0, where d = 2 δt/h². A positive operator results if and only if 0 < d < 1/(1 − ω).

a) Explicit

0 0.5 1

0 1 2 3 4

t

T expliciet

d = 1

-1e+07 -5e+06 0 5e+06 1e+07

0 1 2 3 4

t

T expliciet

d = 1.02

For ω = 0 the Fourier amplification factor is given by g(θ) = 1 − d(1 − cos θ).

Hence, the explicit method is Fourier stable if and only if 0 ≤ d ≤ 1. For the chosen value of h this limit corresponds with δt = 0.005. This coincides with the limit of being a positive operator.

b) Crank-Nicolson

0 0.5 1

0 1 2 t 3 4

T ω = 0.5

d = 10

0 0.5 1

0 1 2 t 3 4

T ω = 0.5

d = 100

For general ω the amplification factor reads

g(θ) = 1 − (1 − ω) d (1 − cos θ) 1 + ω d (1 − cos θ) ,

hence the method is stable if and only if 0 ≤ d ≤ 1/(1 − 2ω) in case ω < 1/2, and unconditionally stable when ω ≥ 1/2. As a particular case, the second-order Crank-Nicolson method is unconditionally stable.

(9)

For large time steps δt, hence large d, the amplification factor behaves as g(θ) ∼ (ω − 1)/ω. For ω = 1/2 we have g(θ) ∼ −1 as d → ∞. For a slightly larger value of ω the oscillations drastically decrease in amplitude, as shown below for ω = 0.6 and δt = 0.5 (d = 100).

0 0.5 1

0 1 2 t 3 4

T ω = 0.6

d = 100

c) Fully implicit

0 0.5 1

0 1 2 t 3 4

T ω = 1.0

d = 10

0 0.5 1

0 1 2 t 3 4

T ω = 1.0

d = 100

Above it has been shown already that the fully implicit method (ω = 1) is unconditionally stable. It always gives rise to a positive operator as the sign of all neighbouring coefficients is OK. An even easier way is to look at the sign of the amplification factor, which is positive for all frequencies θ, hence no wiggles can occur (these would require g < 0).

(10)

Exercise 5 - Navier–Stokes Part 1: Treatment of div u

ⁿ

Question 0 – ’exact’ solution TFin = 25, δt = 0.025

Time Div U Iter U_bottom U_mid U_top

T= 0.25000E+02 div= 0.61845E-06 1 -0.01981 -0.20055 0.50161

t = 2.500e+01 : streamfunction

−0.09

−0.08

−0.07

−0.06

−0.05

−0.04

−0.03

−0.02

−0.01 0

Questions without div uⁿ

1) ε = 10⁻¹, initial solution p = 0:

a) TFin=25, δt = 0.025

T= 0.20000E+02 div= 0.20544E+02 1 0.00615 -0.00257 0.34755 T= 0.21000E+02 div= 0.20724E+02 2 0.00654 0.01535 0.37700 T= 0.22000E+02 div= 0.20927E+02 1 0.00704 0.03234 0.41041 T= 0.23000E+02 div= 0.21111E+02 1 0.00741 0.04475 0.44161 T= 0.24000E+02 div= 0.21271E+02 2 0.00758 0.05372 0.46994 T= 0.25000E+02 div= 0.21412E+02 2 0.00754 0.06027 0.49501

TFin=100, δt = 0.025

T= 0.95003E+02 div= 0.19522E+02 2 -0.01519 0.09989 0.52061 T= 0.96004E+02 div= 0.19515E+02 2 -0.01506 0.10205 0.51959 T= 0.97004E+02 div= 0.19508E+02 2 -0.01496 0.10409 0.51851 T= 0.98004E+02 div= 0.19500E+02 2 -0.01486 0.10602 0.51739 T= 0.99004E+02 div= 0.19492E+02 2 -0.01479 0.10782 0.51623 T= 0.10000E+03 div= 0.19484E+02 2 -0.01474 0.10948 0.51502

The magnitude of the divergence div u keeps changing linearly with time: div uⁿ = div u⁰− ε n δt ≈ div u⁰− εt. The solution is completely wrong.

(11)

−0.06

−0.04

−0.02 0 0.02 0.04 0.06 0.08 0.1 0.12

−0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

b) δt = 0.0025

T= 0.20000E+02 div= 0.20547E+02 2 0.00681 -0.00953 0.34392 T= 0.21000E+02 div= 0.20692E+02 2 0.00667 0.00436 0.36397 T= 0.22000E+02 div= 0.20867E+02 2 0.00694 0.02018 0.39315 T= 0.23000E+02 div= 0.21044E+02 2 0.00722 0.03355 0.42277 T= 0.24001E+02 div= 0.21207E+02 2 0.00734 0.04410 0.45054 T= 0.25001E+02 div= 0.21354E+02 1 0.00728 0.05215 0.47605

Note that the time step has hardly any influence on div u ∼ −εt, independent of the time step. Also, the solution at 25 seconds does not change that much, although in time it is rapidly changing (compare the solution after 25 seconds with that after 100 seconds).

2) ε = 10⁻³, initial solution p = 0.

Tfin=25, δt = 0.025

T= 0.20000E+02 div= 0.47469E+00 10 -0.01815 -0.19822 0.49768 T= 0.21000E+02 div= 0.49685E+00 10 -0.01806 -0.19805 0.49745 T= 0.22000E+02 div= 0.51888E+00 10 -0.01798 -0.19787 0.49722 T= 0.23000E+02 div= 0.54078E+00 10 -0.01789 -0.19770 0.49698 T= 0.24000E+02 div= 0.56256E+00 10 -0.01781 -0.19752 0.49674 T= 0.25000E+02 div= 0.58420E+00 10 -0.01772 -0.19735 0.49650

The level of div u has dropped two orders of magnitude, proportional to the decrease in ε:

div uⁿ∼ −εt. The solution begins to look like the correct solution.

−0.09

−0.08

−0.07

−0.06

−0.05

−0.04

−0.03

−0.02

−0.01 0

(12)

3) ε = 10⁻³, initial solution p = pold δt = 0.025

T= 0.20000E+02 div= 0.37338E-01 1 -0.02009 -0.20106 0.50118 T= 0.21000E+02 div= 0.37330E-01 1 -0.02009 -0.20106 0.50118 T= 0.22000E+02 div= 0.37329E-01 1 -0.02009 -0.20106 0.50118 T= 0.23000E+02 div= 0.37326E-01 1 -0.02009 -0.20106 0.50118 T= 0.24000E+02 div= 0.37326E-01 1 -0.02009 -0.20106 0.50118 T= 0.25000E+02 div= 0.37325E-01 1 -0.02009 -0.20106 0.50118

The divergence div u is not yet really small.

4) ε = 10⁻¹, initial solution p = pold δt = 0.025

With larger ε the solution becomes worse (as might be expected).

Questions with div uⁿ

5) ε = 10⁻¹, initial solution p = 0 a) δt = 0.025

b) δt = 0.0025

The level of div u has dropped an order of magnitude, proportional to the decrease in δt: div uⁿ∼ −ε δt.

(13)

6) ε = 10⁻³, initial solution p = 0 δt = 0.025

The reduction of ε has led to a proportional reduction of div u.

7) ε = 10⁻³, initial solution p = pold δt = 0.025

The divergence becomes quite small (in fact machine accuracy).

8) ε = 10⁻¹, initial solution p = pold δt = 0.025

Even at higher ε, the divergence div u is at a low level (in fact machine accuracy, hence the irregularities), because each time step at least one Poisson iteration is made (in fact:

has to be made before convergence can be checked), although the convergence criterion has been met already.

9) By comparing Question 1 and Question 8, we notice that even with a crude convergence criterion a good solution can be obtained, provided the term div uⁿ is kept in the right- hand side of the Poisson equation.

(14)

Part 2: Influence of boundary conditions

Questions with prescribed outflow 10) ν = 0.01

−21 0 2 4 6 8 10 12

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

X

horizontal velocity

Y= 0.5

T=15 T= 5 T= 4 T= 3 T= 2 T= 1

Figure 1: Horizontal velocity along center line for ν = 0.01 at outflow condition u = 1.

11) ν = 0.005

−21 0 2 4 6 8 10 12

1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5

X

horizontal velocity

Y= 0.5

T=15 T= 5 T= 4 T= 3 T= 2 T= 1

12) ν = 0.001.

−21 0 2 4 6 8 10 12

1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4

X

horizontal velocity

Y= 0.5

T=15 T= 5 T= 4 T= 3 T= 2 T= 1

(15)

t = 1.500e+01 : velocity

0.2 0.4 0.6 0.8 1 1.2

−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4

−0.2 0 0.2 0.4 0.6 0.8 1 1.2

horizontal velocity

Y

T= 14.9998 X= 9.5312

Figure 4: Left: Absolute velocity and velocity vectors for x ∈ [9, 10]. Right: Horizontal velocity profile at x = 9.53 (i = 32) .

Questions with free outflow 13) ν = 0.01

−21 0 2 4 6 8 10 12

1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5

X

horizontal velocity

Y= 0.5

T=15 T= 5 T= 4 T= 3 T= 2 T= 1

Figure 5: Horizontal velocity at center line for ν = 0.01 with free outflow.

0 0.5 1 1.5

−0.2 0 0.2 0.4 0.6 0.8 1 1.2

horizontal velocity

Y

T= 14.9998 X= 9.5312

Figure 6: Horizontal velocity profile at x = 9.53 (i = 32) for ν = 0.01. It has a more or less parabolic shape.

(16)

14) ν = 0.005

−21 0 2 4 6 8 10 12

1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5

X

horizontal velocity

Y= 0.5

T=15 T= 5 T= 4 T= 3 T= 2 T= 1

15) ν = 0.001

−21 0 2 4 6 8 10 12

1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4

X

horizontal velocity

Y= 0.5

T=15 T= 5 T= 4 T= 3 T= 2 T= 1

0 0.2 0.4 0.6 0.8 1 1.2 1.4

−0.2 0 0.2 0.4 0.6 0.8 1 1.2

horizontal velocity

Y

T= 14.9999 X= 9.5312

Figure 9: Horizontal velocity profile at x = 9.53 (i = 32) for ν = 0.001. The profile does not really look like a parabola.

(17)

−5 0 5 10 15 20 25 1

1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45

X

horizontal velocity

T= 29.9997 Y= 0.5

−101 0 10 20 30 40 50 60

1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5

X

horizontal velocity

T= 49.9995 Y= 0.5

Figure 10: The horizontal velocity at the center line for ν = 0.001 and a free outflow condition is prescribed, at channel length 20 (left) and 50 (right), respectively.

t = 5.000e+01 : velocity

0.2 0.4 0.6 0.8 1 1.2 1.4

0 0.5 1 1.5

−0.2 0 0.2 0.4 0.6 0.8 1 1.2

horizontal velocity

Y

T= 50.0007 X= 49.2188

Figure 11: Left: Absolute velocity and velocity vectors for x ∈ [49, 50]. Right: Horizontal velocity profile near exit. Viscosity ν = 0.001; the parabola has developed now.

Conclusion: Smaller viscosity requires more length for a flow to become fully developed.