Biringen S., Chow C.-Y. An Introduction to Computational Fluid Mechanics by Example

Подождите немного. Документ загружается.

202 VISCOUS FLUID FLOWS

the disturbance in a reference frame moving with the wave speed. This enables us

to deploy all available mesh resolution to capture the nonlinear evolution of one

single wavelength of the disturbance. It is often sufﬁcient to assume the spanwise

z direction to be periodic as well. The assumption of periodicity enables the use

of discrete Fourier transforms, which can be very efﬁciently computed by fast

Fourier transform (FFT) algorithms.

Efﬁcient FFT algorithms require that the periodic function, u(x),isdeﬁnedon

N equidistant grid points, where N is an even number given as N = 2

,where

l is an integer. Generally, to resolve one period of a disturbance wave with

minimum acceptable accuracy, l = 4orN = 16. The discrete Fourier transform

of u(x) can be written as

u(x

) =

n=N /2−1



n=−N /2

ˆu(k

(3.8.65)

Making use of the discrete orthogonality of the Fourier series (e.g., Moin, 2001,

p. 159), we can write the Fourier coefﬁcients in (3.8.65) as

ˆu(k

) =

m=N −1



m=0

u(x

−Ik

(3.8.66)

The discrete Fourier transform is computed at N equidistant points along the

x direction over a computational distance L, separated by a constant grid size

of x,andL = N .x. We can then impose periodicity of the function u(x) by

making u(x

) = u(x

) so that the period ends at x

N −1

and the next period starts

at x

. Hence, we have

= mx (m = 0, 1, ..., N − 1)

2πn

(wave number) n =−N /2, −N /2 +1, ...,0,..., N /2 −1

(3.8.67)

Derivatives along the x direction can also be easily calculated by the pseudo-

spectral method (Orszag, 1972; Fox and Orszag, 1973) according to the following

transforms:

∂u

∂x



n=N /2−1



n=−N /2

ˆu(k

where

∂ ˆu

∂x



= Ik

ˆu(k

)

(3.8.68)

It follows that according to (3.8.68), the ﬁrst derivative of u(x) can be calculated

by Fourier transforming u(x) in the x direction, multiplying it by Ik

, and taking

FLOW STABILITY AND PSEUDO-SPECTRAL METHODS 203

the inverse transform. Similarly, the second derivative can be calculated from

∂

∂x



n=N /2−1



n=−N /2

−k

ˆu(k

where

∂

ˆu

∂x



=−k

ˆu(k

)

(3.8.69)

In computing derivatives using (3.8.68) and (3.8.69), we note that for real u(x)

in the physical space, the Fourier coefﬁcient of the derivative for n =−N/2must

be set to zero to avoid a the calculation of a “noisy” derivative.

As an example, let us consider the ordinary differential equation (Moin, 2001,

p. 191),



−u



−2u = 2 +6sinx − 38 cos 6x,

u(0) = u(2π) = 0

(3.8.70)

We will solve this two-point boundary-value problem by the Fourier transform

method. Writing the derivatives in the Fourier space for wave number k

,we

obtain

−k

ˆu(k

) − Ik

ˆu(k

) − 2ˆu(k

) =

b(k

)

ˆu(k

) =

b(k

)

−k

−Ik

−2

, n = 1, ..., N

(3.8.71)

We note that

b(k

) is the Fourier coefﬁcient of the right-hand-side of (3.8.70),

and the solution will be obtained simply by back transforming (3.8.71). In

Program 3.9, a MATLAB script is given that performs this integration, and the

solution is compared with the exact solution,

u(x) = cos(6x) − 1 (3.8.72)

If this program is executed for different values of N , it becomes apparent that

the number of grid points on which the function is evaluated has a very strong

effect on the solution accuracy, and to have an acceptable solution one would

have at least N = 32 grid points. Figure 3.8.7 presents the solution to (3.8.70)

using 64 collocation points and compares it with the exact solution (3.8.72). As

expected from the periodic nature of u(x), the numerical solution is exact for

this problem.

A more interesting application of the Fourier method can be illustrated by its

application to the two-dimensional convection-diffusion equation,

∂T

∂t

∂T

∂x

∂T

∂y

Pr Re



∂

∂x

∂

∂y



u =−sin(πx) cos(2πy), v = cos(πx) cos(πy) (given velocity ﬁeld) (3.8.73)

x = 0, T = 0; x = 1, T = 1

204 VISCOUS FLUID FLOWS

0 1 2 3 4 5 6 7

−2

−1.8

−1.6

−1.4

−1.2

−1

−0.8

−0.6

−0.4

−0.2

u(x)

FIGURE 3.8.7 Solution of equation (3.8.70) by the Fourier transform method.

According to (3.8.73), temperature, T(x, y) is a periodic function in y with a

period, L = 2; also, Re is the Reynolds number and Pr is the Prandtl number. In

the solution process, ﬁrst, we advance the solution in time using the Euler explicit

scheme and second-order central differences along the nonperiodic x direction.

With

n representing the current time level, the resulting ﬁnite difference equation

reads

n+1

i, j

= T

i, j

−u

i, j

t

i+1,j

−T

i−1,j

2x

−v

i, j

t



∂T

∂y



i, j

t

Pr Re

i+1,j

−T

i, j

i−1,j

x

t

Pr Re



∂

∂y



i, j

(3.8.74)

The derivatives along the y direction will be calculated by the pseudo-spectral

Fourier transform method. Program 3.10 presents a MATLAB script that inte-

grates (3.8.74) until it meets the convergence criterion below:



n+1

i, j

−T

i, j



max

≤ e

−5

(3.8.75)

Because the integration is explicit in time, the time step has to comply with both

the diffusive and convective stability criteria. For Re = 1, and Pr = 0.71 (air),

FLOW STABILITY AND PSEUDO-SPECTRAL METHODS 205

0 0.5 1

0.2

0.4

0.6

0.8

1.2

1.4

1.6

1.8

FIGURE 3.8.8 Velocity vectors given by u and v in (3.8.73).

t = 0.0001 was required to ensure numerical stability, whereas for Re = 1

and Pr = 25 (molten glass), t = 0.001 was sufﬁcient to meet both stability

criteria.

Figure 3.8.8 shows the velocity ﬁeld, which is “frozen” in time, whereas in

Fig. 3.8.9 temperature contours for Pr = 0.71 are presented. Comparing these

with the temperature contours for Pr = 25 (Fig. 3.8.10) clearly indicates that

with increasing Pr, temperature diffusion is much stronger, and high temperatures

generated at the right boundary penetrate much farther toward the left (cold)

boundary.

Problem 3.16 Solve equation (3.8.73) using the Euler scheme (forward in

time) with explicit time advancement for the convective terms and implicit time

advancement for the diffusion terms with ﬁrst-order time advancement (forward

in time). Use second-order central differences in space along the x direction. The

problem is periodic in the y direction so that Fourier transforms can be used

along y (see Program 3.10). Because the velocity components u and v are given

(they are not solution dependent), the convective terms are linear. You can use

any solver from Chapter 2, e.g.,

TRID from Program 2.1, for the solution of the

tridiagonal system of equations resulting from the ﬁnite differences along the

nonperiodic x direction. You can also use the MATLAB backslash operator to

solve the resulting tridiagonal matrix equation.

206 VISCOUS FLUID FLOWS

0 0.5 1

0.2

0.4

0.6

0.8

1.2

1.4

1.6

1.8

FIGURE 3.8.9 Temperature contours for Pr = 0.71. Right wall (x = 1) is ‘‘hot’’ and left wall

is ‘‘cold.’’

0 0.5 1

0.2

0.4

0.6

0.8

1.2

1.4

1.6

1.8

FIGURE 3.8.10 Temperature contours for Pr = 25. Right wall (x = 1) is ‘‘hot’’ and left wall

is ‘‘cold.’’

APPENDIX 207

APPENDIX

TABLE 3.A.1 Tabulated Data for Program 3.1 (Number of

Iterations = 31)

ETA F F

∗

∗∗

0003.3206e-01

2.0000e-01 6.6410e-03 6.6408e-02 3.3198e-01

4.0000e-01 2.6560e-02 1.3276e-01 3.3147e-01

6.0000e-01 5.9735e-02 1.9894e-01 3.3008e-01

8.0000e-01 1.0611e-01 2.6471e-01 3.2739e-01

1.0000e+00 1.6557e-01 3.2978e-01 3.2301e-01

1.2000e+00 2.3795e-01 3.9378e-01 3.1659e-01

1.4000e+00 3.2298e-01 4.5626e-01 3.0787e-01

1.6000e+00 4.2032e-01 5.1676e-01 2.9666e-01

1.8000e+00 5.2952e-01 5.7476e-01 2.8293e-01

2.0000e+00 6.5002e-01 6

.2977e-01 2.6675e-01

2.2000e+00 7.8119e-01 6.8131e-01 2.4835e-01

2.4000e+00 9.2229e-01 7.2898e-01 2.2809e-01

2.6000e+00 1.0725e+00 7.7245e-01 2.0645e-01

2.8000e+00 1.2310e+00 8.1151e-01 1.8401e-01

3.0000e+00 1.3968e+00 8.4604e-01 1.6136e-01

3.2000e+00 1.5691e+00 8.7608e-01 1.3913e-01

3.4000e+00 1.7469e+00 9.0176e-01 1.1788e-01

3.6000e+00 1.9295e+00 9.

2333e-01 9.8086e-02

3.8000e+00 2.1160e+00 9.4112e-01 8.0126e-02

4.0000e+00 2.3057e+00 9.5552e-01 6.4234e-02

4.2000e+00 2.4980e+00 9.6696e-01 5.0520e-02

4.4000e+00 2.6924e+00 9.7587e-01 3.8973e-02

4.6000e+00 2.8882e+00 9.8268e-01 2.9484e-02

4.8000e+00 3.0853e+00 9.8779e-01 2.1871e-02

5.0000e+00 3.2833e+00 9.9154e-01 1.5907e-02

5.2000e+00 3.4819e

+00 9.9425e-01 1.1342e-02

5.4000e+00 3.6809e+00 9.9615e-01 7.9278e-03

5.6000e+00 3.8803e+00 9.9748e-01 5.4321e-03

5.8000e+00 4.0799e+00 9.9838e-01 3.6485e-03

6.0000e+00 4.2796e+00 9.9897e-01 2.4021e-03

6.2000e+00 4.4795e+00 9.9936e-01 1.5502e-03

6.4000e+00 4.6794e+00 9.9961e-01 9.8067e-04

6.6000e+00 4.8793e+00 9.9977e-01 6.0809e-04

6.8000e+

00 5.0793e+00 9.9986e-01 3.6960e-04

7.0000e+00 5.2792e+00 9.9992e-01 2.2019e-04

7.2000e+00 5.4792e+00 9.9996e-01 1.2859e-04

7.4000e+00 5.6792e+00 9.9998e-01 7.3606e-05

7.6000e+00 5.8792e+00 9.9999e-01 4.1300e-05

7.8000e+00 6.0792e+00 9.9999e-01 2.2714e-05

8.0000e+00 6.2792e+00 1.0000e+00 1.2245e-05

8.2000e+00 6.4792e+00 1.0000e+

00 6.4706e-06

(continues)

208 VISCOUS FLUID FLOWS

TABLE 3.A.1 (continued)

ETA F F

∗

∗∗

8.4000e+00 6.6792e+00 1.0000e+00 3.3516e-06

8.6000e+00 6.8792e+00 1.0000e+00 1.7017e-06

8.8000e+00 7.0792e+00 1.0000e+00 8.4688e-07

9.0000e+00 7.2792e+00 1.0000e+00 4.1314e-07

9.2000e+00 7.4792e+00 1.0000e+00 1.9755e-07

9.4000e+00 7.6792e+00 1.0000e+00 9.2598e-08

9.6000e+00 7.8792e+00 1.0000e

+00 4.2545e-08

9.8000e+00 8.0792e+00 1.0000e+00 1.9161e-08

1.0000e+01 8.2792e+00 1.0000e+00 8.4589e-09

TABLE 3.A.2 Tabulated Data for Program 3.2

ETA THETA THETA THETA

(WATER) (AIR) (MERCURY)

02.4877e+00 8.4428e-01 1.5900e-01

1.0000e-01 2.4802e+00 8.4349e-01 1.5895e-01

2.0000e-01 2.4579e+00 8.4113e-01 1.5880e-01

3.0000e-01 2.4208e+00 8.3719e-01 1.5856e-01

4.0000e-01 2.3690e+00 8.3169e-01 1.5822e-01

5.0000e-01 2.3029e+00 8.2462e-01 1.5779e-01

6.0000e-01 2.2231e+00 8.1600e-01 1.5725e-01

7.0000e-01 2.1303e+00 8.0585e-01 1.5663e-01

8.0000e-01 2.0256e+00 7.9418e-01 1.5590e-01

9.0000e-01 1.

9104e+00 7.8103e-01 1.5508e-01

1.0000e+00 1.7862e+00 7.6643e-01 1.5417e-01

1.1000e+00 1.6552e+00 7.5043e-01 1.5317e-01

1.2000e+00 1.5193e+00 7.3309e-01 1.5208e-01

1.3000e+00 1.3809e+00 7.1447e-01 1.5090e-01

1.4000e+00 1.2424e+00 6.9465e-01 1.4964e-01

1.5000e+00 1.1060e+00 6.7371e-01 1.4829e-01

1.6000e+00 9.7404e-01 6.5174e-01 1.4686e-01

1.7000e+00 8

.4847e-01 6.2885e-01 1.4536e-01

1.8000e+00 7.3095e-01 6.0516e-01 1.4379e-01

1.9000e+00 6.2280e-01 5.8078e-01 1.4214e-01

2.0000e+00 5.2488e-01 5.5583e-01 1.4044e-01

2.1000e+00 4.3767e-01 5.3045e-01 1.3867e-01

2.2000e+00 3.6121e-01 5.0478e-01 1.3684e-01

2.3000e+00 2.9521e-01 4.7894e-01 1.3497e-01

2.4000e+00 2.3906e-01 4.5307e-01 1.3305e-01

2.5000e+00 1.9196e-01 4.2732e-01 1.3109e-01

2.6000e+00 1.5296e-01 4

.0179e-01 1.2909e-01

2.7000e+00 1.2104e-01 3.7663e-01 1.2705e-01

2.8000e+00 9.5183e-02 3.5195e-01 1.2499e-01

(continues)

APPENDIX 209

TABLE 3.A.2 (continued)

ETA THETA THETA THETA

(WATER) (AIR) (MERCURY)

2.9000e+00 7.4440e-02 3.2786e-01 1.2291e-01

3.0000e+00 5.7926e-02 3.0446e-01 1.2080e-01

3.1000e+00 4.4868e-02 2.8183e-01 1.1868e-01

3.2000e+00 3.4600e-02 2.6005e-01 1.1654e-01

3.3000e+00 2.6566e-02 2.3919e-01 1.1440e-01

3.4000e+00 2.0305e-02 2.1931e-01 1.1225e-01

3.5000e+00 1.5447e-02 2.0043e-01 1.1009e-01

3.6000e+00 1.1692e-02 1.8259e-01 1.0793e-01

3.7000e+00 8.8015e-03 1.6580e-01 1.0578e-01

8000e+00 6.5873e-03 1.5008e-01 1.0362e-01

3.9000e+00 4.8996e-03 1.3541e-01 1.0148e-01

4.0000e+00 3.6205e-03 1.2179e-01 9.9332e-02

4.1000e+00 2.6570e-03 1.0918e-01 9.7196e-02

4.2000e+00 1.9360e-03 9.7572e-02 9.5068e-02

4.3000e+00 1.4002e-03 8.6916e-02 9.2950e-02

4.4000e+00 1.0050e-03 7.7177e-02 9.0842e-02

4.5000e+00 7.1570e-04 6.8310e-02 8.8745e-02

4.6000e+00 5.0565e-04 6.0268e-02 8.6660e-02

4.7000e+

00 3.5436e-04 5.3004e-02 8.4588e-02

4.8000e+00 2.4629e-04 4.6466e-02 8.2528e-02

4.9000e+00 1.6976e-04 4.0603e-02 8.0482e-02

5.0000e+00 1.1603e-04 3.5367e-02 7.8450e-02

5.1000e+00 7.8628e-05 3.0707e-02 7.6433e-02

5.2000e+00 5.2826e-05 2.6576e-02 7.4430e-02

5.3000e+00 3.5184e-05 2.2926e-02 7.2442e-02

5.4000e+00 2.3229e-05 1.9713e-02 7.0470e-02

5.5000e+00 1.5202e-05 1.6896e-02 6.8514e-02

5.6000e+00 9.

8606e-06 1.4435e-02 6.6573e-02

5.7000e+00 6.3391e-06 1.2292e-02 6.4650e-02

5.8000e+00 4.0388e-06 1.0433e-02 6.2743e-02

5.9000e+00 2.5501e-06 8.8259e-03 6.0853e-02

6.0000e+00 1.5956e-06 7.4421e-03 5.8980e-02

6.1000e+00 9.8933e-07 6.2546e-03 5.7125e-02

6.2000e+00 6.0782e-07 5.2393e-03 5.5288e-02

6.3000e+00 3.7002e-07 4.3743e-03 5.3469e-02

6.4000e+00 2.2319e-07 3.6400e-03 5.1668e-02

6.5000e+00 1.3338e-07 3.

0189e-03 4.9885e-02

(continues)

210 VISCOUS FLUID FLOWS

TABLE 3.A.2 (continued)

ETA THETA THETA THETA

(WATER) (AIR) (MERCURY)

6.6000e+00 7.8976e-08 2.4954e-03 4.8121e-02

6.7000e+00 4.6329e-08 2.0558e-03 4.6376e-02

6.8000e+00 2.6925e-08 1.6880e-03 4.4650e-02

6.9000e+00 1.5503e-08 1.3813e-03 4.2943e-02

7.0000e+00 8.8425e-09 1.1265e-03 4.1256e-02

7.1000e+00 4.9964e-09 9.1565e-04 3.9588e-02

7.2000e+00 2.7967e-09 7.4171e-04 3.7939e-02

7.3000e+00 1.5508e-09 5.9876e-04 3.6310e-02

7.4000e+00 8.5179e-10 4.8171e-04 3.4702e-02

5000e+00 4.6344e-10 3.8620e-04 3.3113e-02

7.6000e+00 2.4977e-10 3.0855e-04 3.1544e-02

7.7000e+00 1.3334e-10 2.4566e-04 2.9995e-02

7.8000e+00 7.0506e-11 1.9489e-04 2.8467e-02

7.9000e+00 3.6928e-11 1.5407e-04 2.6959e-02

8.0000e+00 1.9157e-11 1.2135e-04 2.5471e-02

8.1000e+00 9.8438e-12 9.5233e-05 2.4004e-02

8.2000e+00 5.0098e-12 7.4453e-05 2.2557e-02

8.3000e+00 2.5253e-12 5.7982e-05 2.1131e-02

8.4000e+

00 1.2608e-12 4.4972e-05 1.9725e-02

8.5000e+00 6.2343e-13 3.4735e-05 1.8339e-02

8.6000e+00 3.0531e-13 2.6707e-05 1.6974e-02

8.7000e+00 1.4809e-13 2.0435e-05 1.5630e-02

8.8000e+00 7.1137e-14 1.5553e-05 1.4306e-02

8.9000e+00 3.3842e-14 1.1767e-05 1.3003e-02

9.0000e+00 1.5943e-14 8.8403e-06 1.1720e-02

9.1000e+00 7.4374e-15 6.5869e-06 1.0457e-02

9.2000e+00 3.4346e-15 4.8580e-06 9.2149e-03

9.3000e+00 1.

5693e-15 3.5364e-06 7.9929e-03

9.4000e+00 7.0857e-16 2.5298e-06 6.7911e-03

9.5000e+00 3.1530e-16 1.7659e-06 5.6094e-03

9.6000e+00 1.3737e-16 1.1883e-06 4.4478e-03

9.7000e+00 5.7655e-17 7.5308e-07 3.3062e-03

9.8000e+00 2.2291e-17 4.2640e-07 2.1844e-03

9.9000e+00 6.7570e-18 1.8207e-07 1.0824e-03

1.0000e+01 0 0 0

TABLE 3.A.3 Tabulated Data for Program 3.3 (Velocity Distribution across the Channel Based on 60 Iterations, Volume Flow

Rate = 0.056892

04.1527e-02 7.3396e-02 9.5915e-02 1.0932e-01 1.1378e-01 1.0933e-0 9.5925e-02 7.3407e-02 4.1534e-02 0

04.1353e-02 7.3067e-02 9.5463e-02 1.0879e-01 1.1322e-01 1.0880e-01 9.5474e-02 7.3079e-02 4.1362e-02 0

04.0819e-02 7.2053e-02 9.4069e-02 1.0716e-01 1.1150e-01 1.0716e-01 9.4080e-02 7.2066e-02 4.0828e-02 0

03.9870e-02 7.0251e-02 9.1593e-02 1.0425e-01 1.0845e-01 1.0426e-01 9.1606e-02 7.0265e-02 3.9880e-02 0

03.8410e-02 6.7481e-02 8.7792e-02 9.9793e-02 1.0376e-01 9.9800e-02 8.7805e-02 6.7495e-02 3.8421e-02 0

03.6287e-02 6.3459e-02 8.

2285e-02 9.3345e-02 9.6995e-02 9.3353e-02 8.2298e-02 6.3473e-02 3.6297e-02 0

03.3269e-02 5.7765e-02 7.4521e-02 8.4283e-02 8.7492e-02 8.4290e-0 7.4533e-02 5.7779e-02 3.3279e-02 0

02.9007e-02 4.9781e-02 6.3718e-02 7.1743e-02 7.4367e-02 7.1750e-02 6.3729e-02 4.9793e-02 2.9015e-02 0

02.2930e-02 3.8579e-02 4.8782e-02 5.4568e-02 5.6447e-02 5.4573e-02 4.8790e-02 3.8589e-02 2.2936e-02 0

01.4007e-02 2.2730e-02 2.8205e-02 3.1261e-02 3.2248e-02 3.1264e-02 2.8210e-02 2.2735e-02 1.4011e-02 0

0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0

211