2. Unsteady flow due to sudden movement of the plate

#### 2.1. Governing equations

An infinite length plate is set in a stationary fluid as an initial condition. Let us consider the situation that the infinite length plate suddenly moves along its parallel direction at a constant speed uw. This problem was first solved by Stokes [2] in his famous treatment of the pendulum. Since Lord Rayleigh [3] also treated this flow, it is often called the Rayleigh problem in the literature. One takes that x is the plate movement direction and y is distance from the plate. Since the velocity component perpendicular to the plate v is zero, the momentum equation is simplified and is shown as a diffusion equation

$$\frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial y^2} \tag{1}$$

Here, u is the velocity component parallel to the plate direction, t is the time, and ν is the kinematic viscosity. The boundary conditions for this partial differential equation are as follows:

$$\begin{cases} y = 0: & \quad \mu = \mu\_w \\ y \to \ast : & \quad \mu \to 0 \end{cases} \tag{2}$$

In order to reduce the partial differential equation to an ordinary equation, the following dimensionless velocity U and the similar variable η are introduced

$$
\mu = \mu\_w \mathcal{U}(\eta), \quad \eta = \frac{y}{2\sqrt{\nu t}} \tag{3}
$$

Then, the following ordinary differential equation can be obtained

$$\frac{d^2\mathcal{U}}{d\eta^2} + 2\eta \frac{d\mathcal{U}}{d\eta} = 0\tag{4}$$

The boundary condition for the ordinary differential equation is as follows using the similar variable η instead of y:

$$\begin{cases} \eta = 0 : \quad \mathcal{U} = 1\\ \eta \to \stackrel{\sim}{\sim} : \quad \mathcal{U} \to 0 \end{cases} \tag{5}$$

As a consequence, one needs to solve this boundary value problem. The theoretical solution can be easily obtained and expressed using the error function

$$\mathcal{U} = \frac{\mu}{\mu\_w} = 1 - \text{erf}(\eta) = 1 - \frac{2}{\sqrt{\pi}} \int\_0^{\eta} \exp\left(-\xi^2\right) d\xi \tag{6}$$

The velocity profile is shown in Figure 1.

be handled and it is valid when the flow velocity is much slower than the sound velocity and/ or the temperature difference in the fluid is small enough to consider the thermal expansion coefficient is independent to the temperature. The former situation is known the low Mach

Another simplification on the incompressible flows is the reduction of dimension due to the characteristic of similarity and periodicity. For the boundary layer flows such as the Blasius flow, the stagnation-point flow, and the von Kármán rotating disk flow have the similar solution where the flow transition from laminar to turbulence does not occur. In those cases, a combined dimensionless variable (similar variable) η is introduced and the velocity distribution can be only a function of η. While for the onset of instability such as the Rayleigh-Bénard convection, the Bénard-Marangoni convection, and the Taylor-Couette flow, the periodic characteristic of flow structure is observed. At the stage of onset of instability, the non-linear term is negligible and therefore the function of flow field is separated into the amplitude part and periodic part, respectively. This makes the effort on numerical analysis to reduce significantly

This chapter consists of three main bodies. First, a numerical technique for solving the boundary value problem called the first Stokes problem or the Rayleigh problem [1] is introduced. The differential equation is transferred into an ordinary equation and it is solved by a finite difference method using the Jacobi method. Second, similar solution of natural convection heat transfer heated from a vertical plate with uniform heat flux is introduced together with the method how to obtain the system of ordinary differential equations. The obtained Nusselt numbers are compared with some previous studies. Third, for example, of the linear stability analysis, one shows that the HSMAC method can be applied to obtain the critical values for the onset of secondary flow such

as the Taylor-Couette flow. The Eigen functions of flow and pressure fields are visualized.

An infinite length plate is set in a stationary fluid as an initial condition. Let us consider the situation that the infinite length plate suddenly moves along its parallel direction at a constant speed uw. This problem was first solved by Stokes [2] in his famous treatment of the pendulum. Since Lord Rayleigh [3] also treated this flow, it is often called the Rayleigh problem in the literature. One takes that x is the plate movement direction and y is distance from the plate. Since the velocity component perpendicular to the plate v is zero, the momentum equation is

> ∂u <sup>∂</sup><sup>t</sup> <sup>¼</sup> <sup>ν</sup>

∂2 u

Here, u is the velocity component parallel to the plate direction, t is the time, and ν is the kinematic viscosity. The boundary conditions for this partial differential equation are as

<sup>∂</sup>y<sup>2</sup> (1)

number approximation, while the latter one the Boussinesq approximation.

244 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

and also to contribute the augmentation of accuracy of the results.

2. Unsteady flow due to sudden movement of the plate

2.1. Governing equations

follows:

simplified and is shown as a diffusion equation

Figure 1. Velocity profile.

### 2.2. Numerical method for solving the ordinary differential equation using finite difference method

For numerical solution, it is necessary to define the range of η, as recognized from Figure 1, η = 4 is enough. Hence, the boundary condition shown below is used instead of Eq. (5)

$$\begin{cases} \eta = 0 : & \mathcal{U} = 1 \\ \eta = 4 : & \mathcal{U} = 0 \end{cases} \tag{7}$$

η ¼ 0 : U<sup>1</sup> ¼ 1 η ¼ 4 : UN ¼ 0

In the following, the case of N = 7 is considered, for example. By substituting i = 2–6 into Eq. (8),

1

0

U<sup>2</sup> U<sup>3</sup> U<sup>4</sup> U<sup>5</sup> U<sup>6</sup> 1

�α2U<sup>1</sup>

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer

http://dx.doi.org/10.5772/intechopen.72263

1

CCCCCCCCCA

1

0

U<sup>2</sup> U<sup>3</sup> U<sup>4</sup> U<sup>5</sup> U<sup>6</sup> 1

CCCCCCCCCA

(11)

CCCCCCCCCA

BBBBBBBBB@

�β6U<sup>7</sup>

0 0 0

0

BBBBBBBBB@

1

CCCCCCCCCA

CCCCCCCCCA ¼

CCCCCCCCCA

This kind of tridiagonal matrix is often seen and can be solved by a direct numerical method, such as Tomas method. However, the rank of the matrix is usually extremely large and one

In general, the rank of the matrix appearing in computational fluid dynamics (CFD) is large and iterative methods such as Jacobi, Gauss-Seidel, or successive over relaxation (SOR) methodare employed. In this subsection, the Jacobi method is explained. The matrix can be divided into

0

BBBBBBBBB@

�α2U<sup>1</sup>

�β<sup>N</sup>�<sup>1</sup>UN

0 0 0

0

BBBBBBBBB@

In the Jacobi method, only the diagonal part is put in the left-hand side (n + 1 step), while the

BBBBBBBBB@

(9)

247

(10)

(

�2 β<sup>2</sup> 0 00 α<sup>3</sup> � 2 β<sup>3</sup> 0 0 0 α<sup>4</sup> � 2 β<sup>4</sup> 0 0 0 α<sup>5</sup> � 2 β<sup>5</sup> 000 α<sup>6</sup> � 2

introduces an iterative method for solving the king-size matrix.

1

0

U<sup>2</sup> U<sup>3</sup> U<sup>4</sup> U<sup>5</sup> U<sup>6</sup> 1

CCCCCCCCCA þ

CCCCCCCCCA

1

0

U<sup>2</sup> U<sup>3</sup> U<sup>4</sup> U<sup>5</sup> U<sup>6</sup> 1

CCCCCCCCCA ¼

CCCCCCCCCA

lower and upper parts are moved to the right-hand side (n step)

BBBBBBBBB@

BBBBBBBBB@

the following simultaneous equation is obtained:

0

BBBBBBBBB@

2.3. Iterative method for matrix solver

0

BBBBBBBBB@

0

þ

BBBBBBBBB@

three parts of lower, diagonal, and upper as follows:

As illustrated in Figure 2, in which vertical and horizontal axes are exchanged from Figure 1, one needs to obtain each value of dimensionless velocity numerically. The approximated velocity profile is expressed by connecting these values smoothly. For simplicity, the intervals between neighboring two points are the same and it is noted as Δη. When the second-order central difference method is used, Eq. (4) is as follows:

$$\frac{\mathcal{U}\_{i+1} - 2\mathcal{U}\_i + \mathcal{U}\_{i-1}}{\left(\Delta \eta\right)^2} + 2\eta\_i \frac{\mathcal{U}\_{i+1} - \mathcal{U}\_{i-1}}{2\left(\Delta \eta\right)} = 0, \quad (i = 2, 3, \dots \text{ $N - 1$ })$$

Here, N is total number of grids and in this chapter, the first grid point starts from 1 as its definition. The above equation becomes

$$\underbrace{\left\{1-\eta\_{i}(\Delta\eta)\right\}}\_{a\_{i}}\mathcal{U}\_{i-1}-2\mathcal{U}\_{i}+\underbrace{\left\{1+\eta\_{i}(\Delta\eta)\right\}}\_{f\_{i}}\mathcal{U}\_{i+1}=0,\quad(i=2,3,\cdots N-1).\tag{8}$$

Here, <sup>η</sup><sup>i</sup> <sup>¼</sup> ð Þ <sup>i</sup> � <sup>1</sup> Δη � � and <sup>α</sup><sup>i</sup> and <sup>β</sup><sup>i</sup> are coefficients determined by the number of grids. The boundary condition (7) is modified

Figure 2. Equidistant grids discretized.

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer http://dx.doi.org/10.5772/intechopen.72263 247

$$\begin{cases} \eta = 0: & \mathcal{U}\_1 = 1 \\ \eta = 4: & \mathcal{U}\_N = 0 \end{cases} \tag{9}$$

In the following, the case of N = 7 is considered, for example. By substituting i = 2–6 into Eq. (8), the following simultaneous equation is obtained:

$$
\begin{pmatrix} -2 & \beta\_2 & 0 & 0 & 0 \\ \alpha\_3 & -2 & \beta\_3 & 0 & 0 \\ 0 & \alpha\_4 & -2 & \beta\_4 & 0 \\ 0 & 0 & \alpha\_5 & -2 & \beta\_5 \\ 0 & 0 & 0 & \alpha\_6 & -2 \end{pmatrix} \begin{pmatrix} U\_2 \\ U\_3 \\ U\_4 \\ U\_5 \\ U\_6 \end{pmatrix} = \begin{pmatrix} -\alpha\_2 \boldsymbol{\mu} \boldsymbol{I}\_1 \\ 0 \\ 0 \\ 0 \\ 0 \\ -\beta\_6 \boldsymbol{\mu} \boldsymbol{I}\_7 \end{pmatrix} \tag{10}
$$

This kind of tridiagonal matrix is often seen and can be solved by a direct numerical method, such as Tomas method. However, the rank of the matrix is usually extremely large and one introduces an iterative method for solving the king-size matrix.

#### 2.3. Iterative method for matrix solver

2.2. Numerical method for solving the ordinary differential equation using finite

η = 4 is enough. Hence, the boundary condition shown below is used instead of Eq. (5)

�

central difference method is used, Eq. (4) is as follows:

246 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

<sup>2</sup> þ 2η<sup>i</sup>

Ui�<sup>1</sup> � 2Ui þ 1 þ η<sup>i</sup>

Uiþ<sup>1</sup> � 2Ui þ Ui�<sup>1</sup> ð Þ Δη

definition. The above equation becomes

1 � η<sup>i</sup> ð Þ <sup>Δ</sup><sup>η</sup> � � |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} <sup>α</sup><sup>i</sup>

boundary condition (7) is modified

Figure 2. Equidistant grids discretized.

Here, η<sup>i</sup> ¼ ð Þ i � 1 Δη

For numerical solution, it is necessary to define the range of η, as recognized from Figure 1,

η ¼ 0 : U ¼ 1 η ¼ 4 : U ¼ 0

As illustrated in Figure 2, in which vertical and horizontal axes are exchanged from Figure 1, one needs to obtain each value of dimensionless velocity numerically. The approximated velocity profile is expressed by connecting these values smoothly. For simplicity, the intervals between neighboring two points are the same and it is noted as Δη. When the second-order

Uiþ<sup>1</sup> � Ui�<sup>1</sup>

Here, N is total number of grids and in this chapter, the first grid point starts from 1 as its

ð Þ <sup>Δ</sup><sup>η</sup> � � |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} βi

<sup>2</sup>ð Þ <sup>Δ</sup><sup>η</sup> <sup>¼</sup> <sup>0</sup>, ið Þ <sup>¼</sup> <sup>2</sup>; <sup>3</sup>; <sup>⋯</sup><sup>N</sup> � <sup>1</sup>

� � and α<sup>i</sup> and β<sup>i</sup> are coefficients determined by the number of grids. The

Uiþ<sup>1</sup> ¼ 0, ið Þ ¼ 2; 3; ⋯N � 1 : (8)

(7)

difference method

In general, the rank of the matrix appearing in computational fluid dynamics (CFD) is large and iterative methods such as Jacobi, Gauss-Seidel, or successive over relaxation (SOR) methodare employed. In this subsection, the Jacobi method is explained. The matrix can be divided into three parts of lower, diagonal, and upper as follows:

$$
\begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ \alpha\_3 & 0 & 0 & 0 & 0 \\ 0 & \alpha\_4 & 0 & 0 & 0 \\ 0 & 0 & \alpha\_5 & 0 & 0 \\ 0 & 0 & 0 & \alpha\_6 & 0 \end{pmatrix} \begin{pmatrix} U\_2 \\ U\_3 \\ U\_4 \\ U\_5 \\ U\_6 \\ U\_6 \end{pmatrix} + \begin{pmatrix} -2 & 0 & 0 & 0 & 0 \\ 0 & -2 & 0 & 0 & 0 \\ 0 & 0 & -2 & 0 & 0 \\ 0 & 0 & 0 & -2 & 0 \\ 0 & 0 & 0 & 0 & -2 \end{pmatrix} \begin{pmatrix} U\_2 \\ U\_3 \\ U\_4 \\ U\_5 \\ U\_6 \end{pmatrix} \\
$$

$$
+ \begin{pmatrix} 0 & \beta\_2 & 0 & 0 & 0 \\ 0 & 0 & \beta\_3 & 0 & 0 \\ 0 & 0 & 0 & \beta\_4 & 0 \\ 0 & 0 & 0 & 0 & \beta\_5 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} U\_2 \\ U\_3 \\ U\_4 \\ U\_5 \\ U\_6 \\ U\_7 \end{pmatrix} = \begin{pmatrix} -a\_2 \, \boldsymbol{I}\_1 \\ 0 \\ 0 \\ 0 \\ 0 \\ -\beta\_{N-1} \, \boldsymbol{I}\_N \end{pmatrix} \tag{11}
$$

In the Jacobi method, only the diagonal part is put in the left-hand side (n + 1 step), while the lower and upper parts are moved to the right-hand side (n step)

�20 0 0 0 0 � 20 00 0 0 � 200 000 � 2 0 0000 � 2 0 BBBBBB@ 1 CCCCCCA |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} D U<sup>2</sup> U<sup>3</sup> U<sup>4</sup> U<sup>5</sup> U<sup>6</sup> 0 BBBBBB@ 1 CCCCCCA nþ1 |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} U ! nþ1 ¼ �α2U<sup>1</sup> 0 0 0 �β<sup>N</sup>�<sup>1</sup>UN 0 BBBBBB@ 1 CCCCCCA |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} b ! � 00000 α<sup>3</sup> 00 0 0 0 α<sup>4</sup> 00 0 0 0 α<sup>5</sup> 0 0 000 α<sup>6</sup> 0 0 BBBBBB@ 1 CCCCCCA |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} L U<sup>2</sup> U<sup>3</sup> U<sup>4</sup> U<sup>5</sup> U<sup>6</sup> 0 BBBBBB@ 1 CCCCCCA n |fflfflfflfflffl{zfflfflfflfflffl} U ! n � 0 β<sup>2</sup> 000 0 0 β<sup>3</sup> 0 0 000 β<sup>4</sup> 0 0000 β<sup>5</sup> 00000 0 BBBBBB@ 1 CCCCCCA |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} U U<sup>2</sup> U<sup>3</sup> U<sup>4</sup> U<sup>5</sup> U<sup>6</sup> 0 BBBBBB@ 1 CCCCCCA n |fflfflfflfflffl{zfflfflfflfflffl} U ! n (12)

Here, n is the old iteration step and n + 1 is the new iteration step. Hence, the following equation is repeatedly used:

$$
\overrightarrow{\boldsymbol{U}}^{n+1} = \boldsymbol{D}^{-1} \left[ \overrightarrow{\boldsymbol{b}} - (\boldsymbol{L} + \boldsymbol{U}) \overrightarrow{\boldsymbol{U}}^{n} \right] \tag{13}
$$

shown in Eqs. (15)–(17) together with the boundary condition (18). Here, one defines that x axis is in the vertical direction and its velocity component is u, and y axis is in the direction

> ∂2 u

y ¼ 0 : u ¼ v ¼ 0, q ¼ �kð Þ ∂T=∂y

Here, β is the thermal expansion coefficient, g is the acceleration due to gravity, α is the thermal

First, dimensionless variables, such as velocity and temperature, are set as follows using the

vaxa yaua |ffl{zffl} ½ � 1

> ∂2 U

∂V <sup>∂</sup><sup>Y</sup> <sup>¼</sup> <sup>0</sup>

, V <sup>¼</sup> <sup>v</sup> va

<sup>∂</sup>Y<sup>2</sup> <sup>þ</sup> <sup>g</sup>βð Þ Tb � <sup>T</sup><sup>∞</sup> xa ua 2 |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} ½ � 3

, <sup>θ</sup> <sup>¼</sup> <sup>T</sup> � Tb Ta

> <sup>þ</sup> <sup>g</sup>βTaxa ua 2 |fflfflffl{zfflfflffl} ½ � 4

θ

, U <sup>¼</sup> <sup>u</sup> ua

∂U ∂X þ

y ! ∞ : u ! 0, T ! T<sup>∞</sup>

∂2 T

<sup>∂</sup><sup>y</sup> <sup>¼</sup> <sup>0</sup> (15)

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer

http://dx.doi.org/10.5772/intechopen.72263

<sup>∂</sup>y<sup>2</sup> <sup>þ</sup> <sup>g</sup>βð Þ <sup>T</sup> � <sup>T</sup><sup>∞</sup> (16)

<sup>∂</sup>y<sup>2</sup> (17)

(18)

249

(19)

∂u ∂x þ ∂v

perpendicular to the vertical plate and its velocity component is v.

u ∂u ∂x þ v ∂u <sup>∂</sup><sup>y</sup> <sup>¼</sup> <sup>ν</sup>

�

unknown reference value denoted with subscripts a and b:

Equation (19) is substituted into Eqs. (15)–(18), and one gets

vaxa yaua |ffl{zffl} ½ � 1

<sup>V</sup> <sup>∂</sup><sup>U</sup> <sup>∂</sup><sup>Y</sup> <sup>¼</sup> <sup>ν</sup>xa ya <sup>2</sup>ua |ffl{zffl} ½ � 2

, Y <sup>¼</sup> <sup>y</sup> ya

<sup>X</sup> <sup>¼</sup> <sup>x</sup> xa

<sup>U</sup> <sup>∂</sup><sup>U</sup> ∂X þ

diffusivity, k is the thermal conductivity, and T is the temperature.

u ∂T ∂x þ v ∂T <sup>∂</sup><sup>y</sup> <sup>¼</sup> <sup>α</sup>

Continuity of mass

Momentum equation

Energy equation

Boundary equation

3.3. Non-dimensionalization

This is equivalent to the following equation:

$$\mathcal{U}\mathcal{U}\_{i}^{n+1} = \frac{1}{2} (\alpha\_i \mathcal{U}\_{i-1}^n + \beta\_i \mathcal{U}\_{i+1}^n), \quad (i = 2, 3, 4, 5, 6) \tag{14}$$

By using Eq. (9), Eq. (14) is computed repeatedly and then the value of each grid gradually converges to a certain solution. The Gauss-Seidel and SOR methods are known as the faster convergence method.

### 3. Similarity solution for natural convection heated from a vertical plate

#### 3.1. Introduction

In this section, let us consider the natural convection heat transfer for a vertical plate heated with uniform heat flux in the wide range of Prandtl number from zero to infinity. In order to explain the numerical method as how to solve the governing equations, one assumes that the flow and temperature fields formed in the vicinity of the heated plate have a similarity and then one introduces the finite difference method to obtain numerical results.

#### 3.2. Governing equations

One assumes that the flow is incompressible laminar and boundary layer equations are used in this analysis. The governing equations with presuming the Boussinesq approximation are shown in Eqs. (15)–(17) together with the boundary condition (18). Here, one defines that x axis is in the vertical direction and its velocity component is u, and y axis is in the direction perpendicular to the vertical plate and its velocity component is v.

Continuity of mass

�20 0 0 0 0 � 20 00 0 0 � 200 000 � 2 0 0000 � 2


0 β<sup>2</sup> 000 0 0 β<sup>3</sup> 0 0 000 β<sup>4</sup> 0 0000 β<sup>5</sup> 00000


equation is repeatedly used:

convergence method.

3.1. Introduction

3.2. Governing equations

1

0

248 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

BBBBBB@

U<sup>2</sup> U<sup>3</sup> U<sup>4</sup> U<sup>5</sup> U<sup>6</sup>


U<sup>2</sup> U<sup>3</sup> U<sup>4</sup> U<sup>5</sup> U<sup>6</sup> 1

nþ1

¼

0

BBBBBB@

1

CCCCCCA

�

0

BBBBBB@

00000 α<sup>3</sup> 00 0 0 0 α<sup>4</sup> 00 0 0 0 α<sup>5</sup> 0 0 000 α<sup>6</sup> 0


1

0

BBBBBB@

U<sup>2</sup> U<sup>3</sup> U<sup>4</sup> U<sup>5</sup> U<sup>6</sup>


1

n

CCCCCCA

(13)

CCCCCCA

�β<sup>N</sup>�<sup>1</sup>UN


Here, n is the old iteration step and n + 1 is the new iteration step. Hence, the following

Un iþ1

By using Eq. (9), Eq. (14) is computed repeatedly and then the value of each grid gradually converges to a certain solution. The Gauss-Seidel and SOR methods are known as the faster

3. Similarity solution for natural convection heated from a vertical plate

In this section, let us consider the natural convection heat transfer for a vertical plate heated with uniform heat flux in the wide range of Prandtl number from zero to infinity. In order to explain the numerical method as how to solve the governing equations, one assumes that the flow and temperature fields formed in the vicinity of the heated plate have a similarity and

One assumes that the flow is incompressible laminar and boundary layer equations are used in this analysis. The governing equations with presuming the Boussinesq approximation are

�ð Þ L þ U U !<sup>n</sup> � �

� �, ið Þ <sup>¼</sup> <sup>2</sup>; <sup>3</sup>; <sup>4</sup>; <sup>5</sup>; <sup>6</sup> (14)

<sup>¼</sup> <sup>D</sup>�<sup>1</sup> <sup>b</sup> !

<sup>i</sup>�<sup>1</sup> <sup>þ</sup> <sup>β</sup><sup>i</sup>

CCCCCCA


n

1

CCCCCCA

U !nþ<sup>1</sup>

<sup>2</sup> <sup>α</sup>iU<sup>n</sup>

then one introduces the finite difference method to obtain numerical results.

(12)

CCCCCCA

1

0

BBBBBB@

CCCCCCA

This is equivalent to the following equation:

U<sup>n</sup>þ<sup>1</sup> <sup>i</sup> <sup>¼</sup> <sup>1</sup>

0

BBBBBB@

0

BBBBBB@

�

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \tag{15}$$

Momentum equation

$$
\mu \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \nu \frac{\partial^2 u}{\partial y^2} + \text{g}\beta (T - T\_{\text{\textquotedblleft}}) \tag{16}
$$

Energy equation

$$
\mu \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} = \alpha \frac{\partial^2 T}{\partial y^2} \tag{17}
$$

Boundary equation

$$\begin{cases} y = 0: \quad u = v = 0, \quad q = -k(\partial T/\partial y) \\ y \to \rightsquigarrow \quad u \to 0, \quad T \to T\_{\rightsquigarrow} \end{cases} \tag{18}$$

Here, β is the thermal expansion coefficient, g is the acceleration due to gravity, α is the thermal diffusivity, k is the thermal conductivity, and T is the temperature.

#### 3.3. Non-dimensionalization

First, dimensionless variables, such as velocity and temperature, are set as follows using the unknown reference value denoted with subscripts a and b:

$$X = \frac{x}{x\_a}, \quad Y = \frac{y}{y\_a}, \quad \mathcal{U} = \frac{u}{u\_a}, \quad V = \frac{v}{v\_a}, \quad \Theta = \frac{T - T\_b}{T\_a} \tag{19}$$

Equation (19) is substituted into Eqs. (15)–(18), and one gets

$$\frac{\partial \mathcal{U}}{\partial \mathcal{X}} + \underbrace{\frac{\upsilon\_a \mathbf{x}\_a}{\underbrace{y\_a u\_a}} \frac{\partial \mathcal{V}}{\partial \mathcal{Y}}}\_{[1]} = \mathbf{0}$$

$$\mathcal{U}\frac{\partial \mathcal{U}}{\partial \mathcal{X}} + \underbrace{\frac{\upsilon\_a \mathbf{x}\_a}{y\_a u\_a} \mathcal{V}\frac{\partial \mathcal{U}}{\partial \mathcal{Y}}}\_{[1]} = \underbrace{\frac{\nu \mathbf{x}\_a}{y\_a^2 u\_a} \frac{\partial^2 \mathcal{U}}{\partial \mathcal{Y}^2}}\_{[2]} + \underbrace{\frac{g\beta (T\_b - T\_\Leftrightarrow) \mathbf{x}\_a}{u\_a^2}}\_{[3]} + \underbrace{\frac{g\beta T\_a \mathbf{x}\_a}{u\_a^2} \Theta}\_{[4]}$$

$$\begin{aligned} \underbrace{\mathcal{U}\frac{\partial\Theta}{\partial\mathcal{X}} + \underbrace{\frac{v\_{a}x\_{a}}{y\_{a}u\_{a}}}\_{[1]}}\_{[1]}\underbrace{\delta\Theta}\_{\partial Y} &= \underbrace{\frac{\alpha x\_{a}}{y\_{a}^{2}u\_{a}}}\_{[5]}\frac{\partial^{2}\Theta}{\partial Y^{2}} \\\\ \left\{ \begin{aligned} Y &= 0 &: \quad \mathcal{U} = V = 0, \quad & \underbrace{(qy\_{a})/(kT\_{a})}\_{[6]} = -\partial\theta/\partial Y \\\\ Y &\to \infty &: \quad \mathcal{U} = 0, \quad & \theta = \underbrace{(T\_{\infty} - T\_{b})/T\_{a}}\_{[7]} \end{aligned} \right. \end{aligned}$$

Y ¼ 0 : U ¼ V ¼ 0, ∂θ=∂Y ¼ �1

xa Ra<sup>∗</sup> ð Þ Pr �1=<sup>5</sup> , U <sup>¼</sup> <sup>u</sup>

<sup>k</sup> Ra<sup>∗</sup> ð Þ Pr �1=<sup>5</sup>

qxa

Furthermore, one assumes that the velocity and temperature fields has a similarity along the direction of vertical plate, so one puts X = 1. These equations are useful for analyzing low

> <sup>þ</sup> Pr <sup>d</sup><sup>2</sup> U

� <sup>V</sup> <sup>d</sup><sup>θ</sup> dη þ d2 θ

η ¼ 0 : U ¼ V ¼ 0, dθ=dη ¼ �1

<sup>∗</sup> <sup>¼</sup> <sup>g</sup>βqx<sup>4</sup>

αν<sup>k</sup> , Pr <sup>¼</sup> <sup>ν</sup>

, V <sup>¼</sup> <sup>v</sup> α x Rax <sup>∗</sup> ð Þ Pr <sup>1</sup> 5 ,

α

α xa Ra<sup>∗</sup> ð Þ Pr <sup>2</sup>=<sup>5</sup> ,

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer

http://dx.doi.org/10.5772/intechopen.72263

(29)

<sup>d</sup>η<sup>2</sup> <sup>þ</sup> <sup>θ</sup> <sup>¼</sup> <sup>0</sup> (30)

<sup>d</sup>η<sup>2</sup> <sup>¼</sup> <sup>0</sup> (31)

(27)

251

(28)

(32)

(33)

Y ! ∞ : U ¼ 0, θ ¼ 0

Ra<sup>∗</sup> ð Þ Pr <sup>1</sup>=<sup>5</sup> , <sup>θ</sup> <sup>¼</sup> <sup>T</sup> � <sup>T</sup><sup>∞</sup>

dV <sup>d</sup><sup>η</sup> <sup>¼</sup> <sup>1</sup> <sup>5</sup> <sup>η</sup> dU <sup>d</sup><sup>η</sup> � <sup>3</sup><sup>U</sup>

<sup>d</sup><sup>η</sup> � <sup>V</sup> dU dη

η ! ∞ : U ¼ θ ¼ 0

The dimensionless variables and non-dimensional numbers are defined as follows:

, U <sup>¼</sup> <sup>u</sup> α x Rax <sup>∗</sup> ð Þ Pr <sup>2</sup> 5

<sup>U</sup> dV

U <sup>5</sup> <sup>η</sup> dθ <sup>d</sup><sup>η</sup> � <sup>θ</sup> 

5

5 , Rax

The local Nusselt number can be obtained by the following derivation:

<sup>η</sup> <sup>¼</sup> <sup>y</sup> x Rax <sup>∗</sup> ð Þ Pr �<sup>1</sup>

<sup>θ</sup> <sup>¼</sup> <sup>T</sup> � <sup>T</sup><sup>∞</sup> qx <sup>k</sup> Rax <sup>∗</sup> ð Þ Pr �<sup>1</sup>

The dimensionless variables are summarized as follows:

, Y <sup>¼</sup> <sup>y</sup>

<sup>X</sup> <sup>¼</sup> <sup>x</sup> xa

<sup>V</sup> <sup>¼</sup> <sup>v</sup> α xa

Prandtl number cases and summarized as follows:

Low Prandtl number

Continuity of mass

Momentum equation

Energy equation

Boundary conditions

At the moment stage, xa is recognized as the height of the vertical plate. Putting [3] = 0, and one obtains Tb ¼ T∞.Hence [7] becomes θ ¼ 0.

Putting [6] = 1, and one gets qya kTa <sup>¼</sup> <sup>1</sup> ) Ta <sup>¼</sup> qya k Putting [5] = 1, and one gets <sup>α</sup>xa ya <sup>2</sup>ua <sup>¼</sup> <sup>1</sup> ) ya <sup>¼</sup> <sup>α</sup>xa ua � �<sup>1</sup>=<sup>2</sup> Putting [1] = 1, vaxa yaua <sup>¼</sup> <sup>1</sup> ) va <sup>¼</sup> yaua xa

$$\text{Putting } [4] = 1, \frac{\mathfrak{g}\mathcal{T}\_{\mathfrak{s}}\mathfrak{x}\_{\mathfrak{s}}}{u\_{\mathfrak{s}}^{2}} = 1 \quad \Rightarrow \quad \mathfrak{u}\_{\mathfrak{a}} = \left(\mathfrak{g}\mathcal{T}\_{\mathfrak{a}}\mathbf{x}\_{\mathfrak{a}}\right)^{1/2} = \left(\mathfrak{g}\mathcal{T}\_{\mathfrak{k}}^{\mathfrak{g}\mathfrak{z}}\mathbf{x}\_{\mathfrak{a}}\right)^{1/2} = \left(\mathfrak{g}\mathcal{T}\_{\mathfrak{k}}^{\mathfrak{g}}\left(\frac{\alpha\mathbf{x}\_{\mathfrak{a}}}{u\_{\mathfrak{a}}}\right)^{1/2}\mathbf{x}\_{\mathfrak{a}}\right)^{1/2}$$

$$\mathbf{u} \cdot \mathbf{u}\_d = \left\{ \left( \mathbf{g} \beta \frac{\eta}{k} \right)^2 \alpha \mathbf{x}\_d^{\ 3} \right\}^{1/5} = \left( \frac{\mathbf{g} \beta \mathbf{q} \mathbf{x}\_d^4}{k a^2} \right)^{2/5} \frac{\alpha}{\mathbf{x}\_d} = (\mathbf{R} \mathbf{a}^\* \mathbf{P} \mathbf{r})^{2/5} \frac{\alpha}{\mathbf{x}\_d}, \quad \because \mathbf{R} \mathbf{a}^\* = \frac{\mathbf{g} \beta \mathbf{q} \mathbf{x}\_d^{\ 4}}{k \alpha \mathbf{v}} \tag{20}$$

$$\mathbf{y}\_d = \left(\frac{a\mathbf{x}\_d}{\mathbf{u}\_d}\right)^{1/2} = \left(\frac{a\mathbf{x}\_d}{\left(Ra^\*Pr\right)^{2/5}\frac{a}{\mathbf{x}\_d}}\right)^{1/2} = \left(\frac{\mathbf{x}\_d^{\;2}}{\left(Ra^\*Pr\right)^{2/5}}\right)^{1/2} = \mathbf{x}\_d \left(Ra^\*Pr\right)^{-1/5} \tag{21}$$

$$T\_a = \frac{qy\_a}{k} = \frac{q}{k} \frac{x\_a}{(Ra^\*Pr)^{1/5}} = \frac{q\chi\_a}{k} \left(Ra^\*Pr\right)^{-1/5} \tag{22}$$

$$w\_a = \frac{y\_a u\_a}{\chi\_a} = \frac{\frac{\chi\_a}{(Ra^\*Pr)^{1/5}} (Ra^\*Pr)^{2/5} \frac{\alpha}{\chi\_a}}{\chi\_a} = (Ra^\*Pr)^{1/5} \frac{\alpha}{\chi\_a} \tag{23}$$

In the above process, finally one obtains the dimensionless equations as follows:

$$
\frac{\partial \mathcal{U}}{\partial X} + \frac{\partial V}{\partial Y} = 0 \tag{24}
$$

$$
\Delta U \frac{\partial \mathcal{U}}{\partial X} + V \frac{\partial \mathcal{U}}{\partial Y} = Pr \frac{\partial^2 \mathcal{U}}{\partial Y^2} + \theta \tag{25}
$$

$$U\frac{\partial\theta}{\partial X} + V\frac{\partial\theta}{\partial Y} = \frac{\partial^2\theta}{\partial Y^2} \tag{26}$$

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer http://dx.doi.org/10.5772/intechopen.72263 251

$$\begin{cases} Y = 0: \quad \mathcal{U} = V = 0, \quad \partial\theta/\partial Y = -1\\ Y \to \rightsquigarrow \quad \mathcal{U} = 0, \quad \theta = 0 \end{cases} \tag{27}$$

The dimensionless variables are summarized as follows:

$$\begin{aligned} X &= \frac{\mathbf{x}}{\mathbf{x}\_d}, \quad Y = \frac{y}{\mathbf{x}\_d (Ra^\*Pr)^{-1/5}}, \quad \mathcal{U} = \frac{u}{\frac{\alpha}{\mathbf{x}\_d} (Ra^\*Pr)^{2/5}},\\ V &= \frac{v}{\frac{\alpha}{\mathbf{x}\_d} (Ra^\*Pr)^{1/5}}, \quad \Theta = \frac{T - T\_{\infty}}{\frac{q\mathbf{x}\_d}{k} (Ra^\*Pr)^{-1/5}} \end{aligned} \tag{28}$$

Furthermore, one assumes that the velocity and temperature fields has a similarity along the direction of vertical plate, so one puts X = 1. These equations are useful for analyzing low Prandtl number cases and summarized as follows:

#### Low Prandtl number

Continuity of mass

<sup>U</sup> <sup>∂</sup><sup>θ</sup> ∂X þ

250 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

8 >>>><

>>>>:

Putting [6] = 1, and one gets qya

Putting [5] = 1, and one gets <sup>α</sup>xa

ua

∴ ua ¼ gβ

ya <sup>¼</sup> <sup>α</sup>xa ua � �<sup>1</sup>=<sup>2</sup>

Putting [1] = 1, vaxa

Putting [4] = 1, <sup>g</sup>βTaxa

vaxa yaua |ffl{zffl} ½ � 1

<sup>Y</sup> ! <sup>∞</sup> : <sup>U</sup> <sup>¼</sup> <sup>0</sup>, <sup>θ</sup> <sup>¼</sup> ð Þ <sup>T</sup><sup>∞</sup> � Tb <sup>=</sup>Ta |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}

Y ¼ 0 : U ¼ V ¼ 0, qya

At the moment stage, xa is recognized as the height of the vertical plate.

kTa <sup>¼</sup> <sup>1</sup> ) Ta <sup>¼</sup> qya

xa

<sup>2</sup>ua <sup>¼</sup> <sup>1</sup> ) ya <sup>¼</sup> <sup>α</sup>xa

<sup>¼</sup> <sup>g</sup>βqxa

xa

xa Ra<sup>∗</sup> ð Þ Pr <sup>1</sup>=<sup>5</sup> <sup>¼</sup> qxa

Ra<sup>∗</sup> ð Þ Pr <sup>1</sup>=<sup>5</sup> Ra<sup>∗</sup> ð Þ Pr <sup>2</sup>=<sup>5</sup> <sup>α</sup>

xa

∂V

<sup>∂</sup><sup>Y</sup> <sup>¼</sup> Pr <sup>∂</sup><sup>2</sup>

<sup>∂</sup><sup>Y</sup> <sup>¼</sup> <sup>∂</sup><sup>2</sup>

Putting [3] = 0, and one obtains Tb ¼ T∞.Hence [7] becomes θ ¼ 0.

ya

<sup>2</sup> ¼ 1 ) ua ¼ gβTaxa

<sup>¼</sup> <sup>α</sup>xa

Ta <sup>¼</sup> qya

va <sup>¼</sup> yaua xa ¼

Ra<sup>∗</sup> ð Þ Pr <sup>2</sup>=<sup>5</sup> <sup>α</sup>

<sup>k</sup> <sup>¼</sup> <sup>q</sup> k

xa

In the above process, finally one obtains the dimensionless equations as follows:

<sup>U</sup> <sup>∂</sup><sup>U</sup>

∂U ∂X þ

<sup>∂</sup><sup>X</sup> <sup>þ</sup> <sup>V</sup> <sup>∂</sup><sup>U</sup>

<sup>U</sup> <sup>∂</sup><sup>θ</sup>

<sup>∂</sup><sup>X</sup> <sup>þ</sup> <sup>V</sup> <sup>∂</sup><sup>θ</sup>

!<sup>1</sup>=<sup>2</sup>

yaua <sup>¼</sup> <sup>1</sup> ) va <sup>¼</sup> yaua

αxa 3 � �<sup>1</sup>=<sup>5</sup>

q k � �<sup>2</sup> <sup>V</sup> <sup>∂</sup><sup>θ</sup> <sup>∂</sup><sup>Y</sup> <sup>¼</sup> <sup>α</sup>xa ya <sup>2</sup>ua |ffl{zffl} ½ � 5

∂2 θ ∂Y<sup>2</sup>

¼ �∂θ=∂Y

� �=ð Þ kTa |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} ½ � 6

½ � 7

k

� �<sup>1</sup>=<sup>2</sup> <sup>¼</sup> <sup>g</sup><sup>β</sup> qya

4 kα<sup>2</sup> � �<sup>2</sup>=<sup>5</sup>

ua � �<sup>1</sup>=<sup>2</sup>

> α xa

<sup>¼</sup> xa

xa

U

θ

<sup>k</sup> xa � �<sup>1</sup>=<sup>2</sup> <sup>¼</sup> <sup>g</sup><sup>β</sup> <sup>q</sup>

2 Ra<sup>∗</sup> ð Þ Pr <sup>2</sup>=<sup>5</sup>

!<sup>1</sup>=<sup>2</sup>

<sup>¼</sup> Ra<sup>∗</sup> ð Þ Pr <sup>2</sup>=<sup>5</sup> <sup>α</sup>

<sup>¼</sup> Ra<sup>∗</sup> ð Þ Pr <sup>1</sup>=<sup>5</sup> <sup>α</sup>

k αxa ua � �<sup>1</sup>=<sup>2</sup>

xa

xa

<sup>¼</sup> xa Ra<sup>∗</sup> ð Þ Pr �1=<sup>5</sup> (21)

4 <sup>k</sup>αν (20)

(23)

� �<sup>1</sup>=<sup>2</sup>

, <sup>∵</sup>Ra<sup>∗</sup> <sup>¼</sup> <sup>g</sup>βqxa

<sup>k</sup> Ra<sup>∗</sup> ð Þ Pr �1=<sup>5</sup> (22)

xa

<sup>∂</sup><sup>Y</sup> <sup>¼</sup> <sup>0</sup> (24)

<sup>∂</sup>Y<sup>2</sup> <sup>þ</sup> <sup>θ</sup> (25)

<sup>∂</sup>Y<sup>2</sup> (26)

$$\frac{dV}{d\eta} = \frac{1}{5} \left( \eta \frac{d\mathcal{U}}{d\eta} - \mathfrak{G}\mathcal{U} \right) \tag{29}$$

Momentum equation

$$
\mathcal{U}\frac{dV}{d\eta} - V\frac{d\mathcal{U}}{d\eta} + \mathcal{P}r\frac{d^2\mathcal{U}}{d\eta^2} + \theta = 0\tag{30}
$$

Energy equation

$$\frac{dI}{5}\left(\eta \frac{d\theta}{d\eta} - \theta\right) - V\frac{d\theta}{d\eta} + \frac{d^2\theta}{d\eta^2} = 0\tag{31}$$

Boundary conditions

$$\begin{cases} \eta = 0: \quad \mathcal{U} = V = 0, \quad d\theta/d\eta = -1\\ \eta \to \stackrel{\sim}{\sim}: \quad \mathcal{U} = \theta = 0 \end{cases} \tag{32}$$

The dimensionless variables and non-dimensional numbers are defined as follows:

$$\begin{aligned} \eta &= \frac{y}{\text{x}(Ra\_{\text{x}}^{\*}Pr)^{-\frac{1}{5}}}, \quad U = \frac{u}{\frac{\alpha}{\text{x}}(Ra\_{\text{x}}^{\*}Pr)^{\frac{2}{5}}}, \quad V = \frac{v}{\frac{\alpha}{\text{x}}(Ra\_{\text{x}}^{\*}Pr)^{\frac{1}{5}}},\\ \Theta &= \frac{T - T\_{\text{w}}}{\frac{qX}{k}(Ra\_{\text{x}}^{\*}Pr)^{-\frac{1}{5}}}, \quad Ra\_{\text{x}}^{\*} = \frac{g\beta qx^{4}}{\alpha vk}, \quad Pr = \frac{\nu}{\alpha} \end{aligned} \tag{33}$$

The local Nusselt number can be obtained by the following derivation:

$$Nu\_x = \frac{h\_x x}{k} = \frac{qx}{(T\_w - T\_\circ)k} = \frac{qx}{T\_a \theta\_w k} = \underbrace{\frac{\left(\frac{g\beta q}{\alpha^2 k}\right)^{1/5}}{qx^{1/5}}}\_{T\_x - 1} \frac{k}{\theta\_w k} \tag{34}$$

$$= \left(\frac{g\beta q}{\alpha^2 k}\right)^{1/5} \frac{x^{4/5}}{\theta\_w} = \left(\frac{g\beta qx^4}{\alpha^2 k}\right)^{1/5} \frac{1}{\theta\_w} = (Ra\_x^\* Pr)^{1/5} \frac{1}{\theta\_w}$$

$$Nu\_x = (Ra\_x^\* Pr)^{\frac{1}{5}} (\theta|\_{\eta=0})^{-1} = (Ra\_x Nu\_x Pr)^{\frac{1}{5}} (\theta|\_{\eta=0})^{-1} \tag{35}$$

Therefore, the local Nusselt number can be obtained just from the dimensionless temperature at the wall using Eq. (36)

$$\frac{N\mu\_{\chi}}{(R a\_{\mathfrak{x}} Pr)^{\frac{1}{\mathfrak{x}}}} = \left(\left.\Theta\right|\_{\mathfrak{y}=0}\right)^{-\frac{5}{4}}\tag{36}$$

Nux ð Þ Rax 1 4 ¼ θj η¼0 �<sup>5</sup> 4

the cases of Pr ≥ 1, while the right-hand side ones the cases of Pr ≤ 1

of discrepancy exists. In this study, the boundary condition for Pr ! 0

η ! ∞ : U ¼ θ ¼ 0

Figure 3 shows the numerical result for the various Prandtl number cases. The upper figures indicate the vertical velocity and lower ones the temperature. The left-hand side figures show

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer

http://dx.doi.org/10.5772/intechopen.72263

Table 1 shows the summary of the local Nusselt number for various Prandtl number cases together with the reference of Churchill and Ozoe for comparison [4]. The agreement is quite good except for the extreme cases such as Pr ! 0 and ∞. In such extreme cases, a small amount

η ¼ 0 : dU=dη ¼ V ¼ 0, dθ=dη ¼ �1

Figure 3. Vertical velocity and temperature distributions for various Prandtl numbers. The left-hand side indicates high

Prandtl number cases while the right-hand side low Prandtl number cases.

3.4. Numerical results

and that for Pr ! ∞

(42)

253

#### High Prandtl number

If the Prandtl number is higher than unity, the following equations are useful:

Continuity of mass

$$\frac{dV}{d\eta} = \frac{1}{5} \left( \eta \frac{d\mathcal{U}}{d\eta} - 3\mathcal{U} \right) \tag{37}$$

Momentum equation

$$
\rho \mathcal{U} \frac{dV}{d\eta} - V \frac{d\mathcal{U}}{d\eta} + Pr \frac{d^2 \mathcal{U}}{d\eta^2} + Pr \theta = 0 \tag{38}
$$

Energy equation

$$\frac{dI}{5}\left(\eta\frac{d\theta}{d\eta}-\theta\right)-V\frac{d\theta}{d\eta}+\frac{d^2\theta}{d\eta^2}=0\tag{39}$$

Boundary conditions

$$\begin{cases} \eta = 0: \quad \mathcal{U} = V = 0, \quad d\theta/d\eta = -1\\ \eta \to \stackrel{\sim}{\sim}: \quad \mathcal{U} = \theta = 0 \end{cases} \tag{40}$$

The dimensionless variables and non-dimensional numbers are defined as follows:

$$\begin{aligned} \eta &= \frac{y}{\mathfrak{x}(Ra\_{\mathfrak{x}}^{\*})^{-\frac{1}{5}}}, \quad \mathcal{U} = \frac{u}{\frac{\alpha}{\mathfrak{x}}(Ra\_{\mathfrak{x}}^{\*})^{\frac{2}{5}}}, \quad \mathcal{V} = \frac{v}{\frac{\alpha}{\mathfrak{x}}(Ra\_{\mathfrak{x}}^{\*})^{\frac{1}{5}}},\\ \Theta &= \frac{T - T\_{\infty}}{\frac{q\chi}{k}(Ra\_{\mathfrak{x}}^{\*})^{-\frac{1}{5}}}, \quad Ra\_{\mathfrak{x}}^{\*} = \frac{g\beta q \mathfrak{x}^{4}}{\alpha vk}, \quad Pr = \frac{\nu}{\alpha} \end{aligned} \tag{41}$$

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer http://dx.doi.org/10.5772/intechopen.72263 253

$$\frac{N\mu\_{\rm x}}{(Ra\_{x})^{\frac{1}{4}}} = \left(\theta|\_{\eta=0}\right)^{-\frac{1}{4}}\tag{42}$$

#### 3.4. Numerical results

Nux <sup>¼</sup> hxx

252 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

<sup>¼</sup> <sup>g</sup>β<sup>q</sup> α2k � �<sup>1</sup>=<sup>5</sup> x<sup>4</sup>=<sup>5</sup>

Nux ¼ Rax

at the wall using Eq. (36)

High Prandtl number

Continuity of mass

Momentum equation

Energy equation

Boundary conditions

<sup>k</sup> <sup>¼</sup> qx

θ<sup>w</sup>

<sup>5</sup> θj η¼0 � ��<sup>1</sup>

Nux ð Þ RaxPr <sup>1</sup> 4 ¼ θj η¼0 � ��<sup>5</sup> 4

If the Prandtl number is higher than unity, the following equations are useful:

dV <sup>d</sup><sup>η</sup> <sup>¼</sup> <sup>1</sup> <sup>5</sup> <sup>η</sup> dU <sup>d</sup><sup>η</sup> � <sup>3</sup><sup>U</sup> � �

<sup>d</sup><sup>η</sup> � <sup>V</sup> dU dη

<sup>þ</sup> Pr <sup>d</sup><sup>2</sup> U

> � <sup>V</sup> <sup>d</sup><sup>θ</sup> dη þ d2 θ

η ¼ 0 : U ¼ V ¼ 0, dθ=dη ¼ �1

<sup>∗</sup> <sup>¼</sup> <sup>g</sup>βqx<sup>4</sup>

, V <sup>¼</sup> <sup>v</sup> α x Rax <sup>∗</sup> ð Þ<sup>1</sup> 5 ,

α

αν<sup>k</sup> , Pr <sup>¼</sup> <sup>ν</sup>

η ! ∞ : U ¼ θ ¼ 0

The dimensionless variables and non-dimensional numbers are defined as follows:

, U <sup>¼</sup> <sup>u</sup> α x Rax <sup>∗</sup> ð Þ<sup>2</sup> 5

<sup>U</sup> dV

U <sup>5</sup> <sup>η</sup> dθ <sup>d</sup><sup>η</sup> � <sup>θ</sup> � �

�

<sup>η</sup> <sup>¼</sup> <sup>y</sup> x Rax <sup>∗</sup> ð Þ�<sup>1</sup> 5

<sup>θ</sup> <sup>¼</sup> <sup>T</sup> � <sup>T</sup><sup>∞</sup> qx <sup>k</sup> Rax <sup>∗</sup> ð Þ�<sup>1</sup> 5 , Rax

<sup>∗</sup> ð Þ Pr <sup>1</sup>

ð Þ Tw � <sup>T</sup><sup>∞</sup> <sup>k</sup> <sup>¼</sup> qx

<sup>¼</sup> <sup>g</sup>βqx<sup>4</sup> α<sup>2</sup>k � �<sup>1</sup>=<sup>5</sup> 1

Therefore, the local Nusselt number can be obtained just from the dimensionless temperature

Taθwk <sup>¼</sup>

θ<sup>w</sup>

<sup>¼</sup> ð Þ RaxNuxPr <sup>1</sup>

gβq α2k � �<sup>1</sup>=<sup>5</sup> k

<sup>¼</sup> Ra<sup>∗</sup>

qx<sup>1</sup>=<sup>5</sup> |fflfflfflfflffl{zfflfflfflfflffl} Ta �1

<sup>x</sup>Pr � �<sup>1</sup>=<sup>5</sup> <sup>1</sup>

<sup>5</sup> θj η¼0 � ��<sup>1</sup>

qx θwk

θ<sup>w</sup>

<sup>d</sup>η<sup>2</sup> <sup>þ</sup> Pr<sup>θ</sup> <sup>¼</sup> <sup>0</sup> (38)

<sup>d</sup>η<sup>2</sup> <sup>¼</sup> <sup>0</sup> (39)

(34)

(35)

(36)

(37)

(40)

(41)

Figure 3 shows the numerical result for the various Prandtl number cases. The upper figures indicate the vertical velocity and lower ones the temperature. The left-hand side figures show the cases of Pr ≥ 1, while the right-hand side ones the cases of Pr ≤ 1

Table 1 shows the summary of the local Nusselt number for various Prandtl number cases together with the reference of Churchill and Ozoe for comparison [4]. The agreement is quite good except for the extreme cases such as Pr ! 0 and ∞. In such extreme cases, a small amount of discrepancy exists. In this study, the boundary condition for Pr ! 0

$$\begin{cases} \eta = 0: \quad d\mathcal{U}/d\eta = V = 0, \quad d\theta/d\eta = -1\\ \eta \to \stackrel{\sim}{\sim}: \quad \mathcal{U} = \theta = 0 \end{cases}$$

and that for Pr ! ∞

Figure 3. Vertical velocity and temperature distributions for various Prandtl numbers. The left-hand side indicates high Prandtl number cases while the right-hand side low Prandtl number cases.


Table 1. Local Nusselt number for various values of Prandtl number (the upper: present results, the lower: Churchill and Ozoe [4]).

$$\begin{cases} \eta = 0: \quad \mathcal{U} = V = 0, \quad d\theta/d\eta = -1\\ \eta \to \stackrel{\sim}{\sim}: \quad d\mathcal{U}/d\eta = \theta = 0. \end{cases}$$

are used. Owing to this kind of special treatments for the boundary condition of such extreme cases, one can obtain accurate numerical results for the system of ordinary equations. The results between the solution of the present method and that of Le Fevre [5] for the case of constant temperature of heated wall are identical to each other. The value for Pr ! ∞ is 0.5027 and that for Pr = 0 is 0.6004.

### 4. Linear stability of Taylor-Couette flow

#### 4.1. Governing equations

In the text book of Chandrasekar [6], various examples of the linear stability analysis such as the Rayleigh-Bénard convection, the Taylor-Couette flow, and the Rayleigh-Taylor instability were studied extensively. More recently, Koschmieder [7] described the research focusing on the Bénard cells and the Taylor vortices. In this section, only the Taylor-Coette flow is considered. Figure 4 shows the schematic model considered for the Taylor-Couette flow. In this section, the fluid flow inside of the co-axial double cylindrical enclosure is assumed to be incompressible Newtonian, isothermal and axisymmetric. The gray part represents the computational domain. It is known that the stationary secondary flow is generated at a certain condition under the influence of centrifugal force due to the rotation of primary basic flow which is in azimuthal direction. The continuity of mass and momentum equations are shown in the cylindrical coordinate system as follows:

$$\frac{\partial u\_r}{\partial r} + \frac{u\_r}{r} + \frac{\partial u\_z}{\partial z} = 0 \tag{43}$$

∂uz <sup>∂</sup><sup>t</sup> <sup>þ</sup> ur

Figure 4. Schematic model for the Taylor-Couette flow.

4.2. Basic state and linearization

Azimuthal velocity

Pressure

follows:

∂uz ∂r þ uz ∂uz <sup>∂</sup><sup>z</sup> ¼ � <sup>1</sup> r ∂p ∂z þ ν

uθð Þ¼� r

p rð Þ¼ ; z

r2 <sup>2</sup>Ω<sup>2</sup> � <sup>r</sup><sup>2</sup>

r2 <sup>1</sup> � <sup>r</sup><sup>2</sup> 2 ∂2 uz ∂r<sup>2</sup> þ

Here, it is indicated that r is the radial, θ is the azimuthal, and z is the axial components.

the basic states for the azimuthal component of velocity and pressure are as follows:

The cylindrical enclosure is long enough to neglect the top and bottom ends. In that situation,

<sup>1</sup>Ω<sup>1</sup>

<sup>ð</sup> <sup>r</sup>f g <sup>u</sup>θð Þ<sup>r</sup>

r

Here, Ω<sup>1</sup> is the angular velocity at the inner cylinder, Ω<sup>2</sup> is the angular velocity at the outer cylinder, p is the pressure, r is the density, and g is the acceleration due to gravity. In order to derive disturbance equations for the linear stability, the three components of velocity and pressure are represented as a summation of basic state and infinitesimal disturbance as

r þ r2 1r2

2

1 r ∂uz ∂r þ ∂2 uz ∂z<sup>2</sup>

� �

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer

http://dx.doi.org/10.5772/intechopen.72263

255

<sup>2</sup>ð Þ Ω<sup>2</sup> � Ω<sup>1</sup> r2 <sup>1</sup> � <sup>r</sup><sup>2</sup> 2

1

dr � rgz þ p<sup>0</sup> (48)

� g (46)

<sup>r</sup> (47)

$$\frac{\partial u\_r}{\partial t} + u\_r \frac{\partial u\_r}{\partial r} + u\_z \frac{\partial u\_r}{\partial z} - \frac{u\_\theta}{r} = -\frac{1}{\rho} \frac{\partial p}{\partial r} + \nu \left( \frac{\partial^2 u\_r}{\partial r^2} + \frac{1}{r} \frac{\partial u\_r}{\partial r} - \frac{u\_r}{r^2} + \frac{\partial^2 u\_r}{\partial z^2} \right) \tag{44}$$

$$\frac{\partial u\_{\theta}}{\partial t} + u\_r \frac{\partial u\_{\theta}}{\partial r} + u\_z \frac{\partial u\_{\theta}}{\partial z} + \frac{u\_r u\_{\theta}}{r} = \nu \left( \frac{\partial^2 u\_{\theta}}{\partial r^2} + \frac{1}{r} \frac{\partial u\_{\theta}}{\partial r} - \frac{u\_{\theta}}{r^2} + \frac{\partial^2 u\_{\theta}}{\partial z^2} \right) \tag{45}$$

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer http://dx.doi.org/10.5772/intechopen.72263 255

Figure 4. Schematic model for the Taylor-Couette flow.

$$\frac{\partial u\_z}{\partial t} + u\_r \frac{\partial u\_z}{\partial r} + u\_z \frac{\partial u\_z}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial z} + \nu \left( \frac{\partial^2 u\_z}{\partial r^2} + \frac{1}{r} \frac{\partial u\_z}{\partial r} + \frac{\partial^2 u\_z}{\partial z^2} \right) - g \tag{46}$$

Here, it is indicated that r is the radial, θ is the azimuthal, and z is the axial components.

#### 4.2. Basic state and linearization

The cylindrical enclosure is long enough to neglect the top and bottom ends. In that situation, the basic states for the azimuthal component of velocity and pressure are as follows:

Azimuthal velocity

$$\overline{u}\_{\theta}(r) = -\frac{r\_2^2 \Omega\_2 - r\_1^2 \Omega\_1}{r\_1^2 - r\_2^2} r + \frac{r\_1^2 r\_2^2 (\Omega\_2 - \Omega\_1)}{r\_1^2 - r\_2^2} \frac{1}{r} \tag{47}$$

Pressure

η ¼ 0 : U ¼ V ¼ 0, dθ=dη ¼ �1

Table 1. Local Nusselt number for various values of Prandtl number (the upper: present results, the lower: Churchill and

0.456

0.4564 0.456

0.5234 0.524

N/A

0.5495 0.550

0.5631 0.5627

are used. Owing to this kind of special treatments for the boundary condition of such extreme cases, one can obtain accurate numerical results for the system of ordinary equations. The results between the solution of the present method and that of Le Fevre [5] for the case of constant temperature of heated wall are identical to each other. The value for Pr ! ∞ is 0.5027

In the text book of Chandrasekar [6], various examples of the linear stability analysis such as the Rayleigh-Bénard convection, the Taylor-Couette flow, and the Rayleigh-Taylor instability were studied extensively. More recently, Koschmieder [7] described the research focusing on the Bénard cells and the Taylor vortices. In this section, only the Taylor-Coette flow is considered. Figure 4 shows the schematic model considered for the Taylor-Couette flow. In this section, the fluid flow inside of the co-axial double cylindrical enclosure is assumed to be incompressible Newtonian, isothermal and axisymmetric. The gray part represents the computational domain. It is known that the stationary secondary flow is generated at a certain condition under the influence of centrifugal force due to the rotation of primary basic flow which is in azimuthal direction. The continuity of mass and momentum equations are shown in the

<sup>∂</sup><sup>z</sup> <sup>¼</sup> <sup>0</sup> (43)

(44)

(45)

∂ur ∂r þ ur r þ ∂uz

2 <sup>r</sup> ¼ � <sup>1</sup> r ∂p ∂r þ ν ∂2 ur ∂r<sup>2</sup> þ

> 1 r ∂u<sup>θ</sup> <sup>∂</sup><sup>r</sup> � <sup>u</sup><sup>θ</sup> r<sup>2</sup> þ ∂2 uθ ∂z<sup>2</sup>

∂2 uθ ∂r<sup>2</sup> þ 1 r ∂ur <sup>∂</sup><sup>r</sup> � ur r<sup>2</sup> þ ∂2 ur ∂z<sup>2</sup>

η ! ∞ : dU=dη ¼ θ ¼ 0:

Pr 0 0.01 0.1 1 10 100 ∞

0.5970 0.597

0.6694 0.670

N/A 0.4564

254 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

4. Linear stability of Taylor-Couette flow

cylindrical coordinate system as follows:

∂ur ∂r þ uz ∂ur <sup>∂</sup><sup>z</sup> � <sup>u</sup><sup>θ</sup>

> ∂u<sup>θ</sup> ∂r þ uz ∂u<sup>θ</sup> ∂z þ uru<sup>θ</sup> <sup>r</sup> <sup>¼</sup> <sup>ν</sup>

∂ur <sup>∂</sup><sup>t</sup> <sup>þ</sup> ur

> ∂u<sup>θ</sup> <sup>∂</sup><sup>t</sup> <sup>þ</sup> ur

and that for Pr = 0 is 0.6004.

0.7107 0.6922

Nux ð Þ Rax <sup>1</sup> = 4

Nux ð Þ RaxPr <sup>1</sup><sup>=</sup> 4

Ozoe [4]).

4.1. Governing equations

$$\overline{p}(r,z) = \int \frac{\rho \{\overline{u}\_{\theta}(r)\}^2}{r} \, dr - \rho \mathbf{g}z + p\_0 \tag{48}$$

Here, Ω<sup>1</sup> is the angular velocity at the inner cylinder, Ω<sup>2</sup> is the angular velocity at the outer cylinder, p is the pressure, r is the density, and g is the acceleration due to gravity. In order to derive disturbance equations for the linear stability, the three components of velocity and pressure are represented as a summation of basic state and infinitesimal disturbance as follows:

$$\begin{aligned} u\_{\theta}(r,z,t) &= \overline{u}\_{\theta}(r) + \overline{v}'(r,z,t), & u\_{r}(r,z,t) &= \overline{u}'(r,z,t), & u\_{z}(r,z,t) &= \overline{w}'(r,z,t), \\ p(r,z,t) &= \overline{p}(r,z) + \overline{p}'(r,z,t) \end{aligned} \tag{49}$$

After neglecting the second-order disturbance, the following linearized equations are obtained:

$$\frac{\partial \mu'}{\partial r} + \frac{\mu'}{r} + \frac{\partial w'}{\partial z} = 0 \tag{50}$$

<sup>R</sup> <sup>¼</sup> <sup>r</sup> r2

Re<sup>Ω</sup> <sup>¼</sup> <sup>Ω</sup>1r<sup>2</sup>

, <sup>U</sup><sup>θ</sup> <sup>¼</sup> <sup>u</sup><sup>θ</sup>

2 ν

the computational results:

Ta <sup>¼</sup> <sup>4</sup>Ω<sup>1</sup>

2 r4 1 ν2

Ω1r<sup>2</sup>

The boundary conditions are as follows:

, k <sup>¼</sup> <sup>r</sup>2a, <sup>η</sup> <sup>¼</sup> <sup>r</sup><sup>1</sup>

(

, <sup>U</sup><sup>~</sup> ; <sup>V</sup><sup>~</sup> ; <sup>W</sup><sup>~</sup> � � <sup>¼</sup> ð Þ <sup>u</sup>~; <sup>~</sup>v; <sup>w</sup><sup>~</sup>

, <sup>μ</sup> <sup>¼</sup> <sup>Ω</sup><sup>2</sup> Ω<sup>1</sup>

r2

<sup>1</sup> � <sup>μ</sup> � � <sup>1</sup> � <sup>4</sup><sup>μ</sup> � �

found in the recent papers published by the present author [9, 10]

Figure 5. The staggered grids in the radius direction together with the points of each variable definition.

<sup>1</sup> � <sup>η</sup><sup>2</sup> ð Þ<sup>2</sup> <sup>¼</sup> <sup>4</sup>Re <sup>2</sup>

Ω1r<sup>2</sup>

, <sup>P</sup><sup>~</sup> <sup>¼</sup> <sup>~</sup><sup>p</sup>

, S <sup>¼</sup> <sup>s</sup> ν=r<sup>2</sup>

<sup>R</sup> <sup>¼</sup> <sup>η</sup> : <sup>U</sup><sup>~</sup> <sup>¼</sup> <sup>V</sup><sup>~</sup> <sup>¼</sup> <sup>W</sup><sup>~</sup> <sup>¼</sup> 0 Inner wall ð Þ <sup>R</sup> <sup>¼</sup> <sup>1</sup> : <sup>U</sup><sup>~</sup> <sup>¼</sup> <sup>V</sup><sup>~</sup> <sup>¼</sup> <sup>W</sup><sup>~</sup> <sup>¼</sup> 0 Outer wall ð Þ:

After Chandrasekar [6], the following two non-dimensional numbers are introduced to verify

Ω

In this section, it is assumed that S = 0. This indicates that the secondary flow caused by the centrifugal instability is stationary and it contains toroidal vortices. To deal with the simultaneous ordinary differential equations for the boundary value problem, a one-dimensional staggered grid system is employed as shown in Figure 5. All the equations are discretized by the fourth order central difference method with a given wavenumber k using the HSMAC method [8] during which Re<sup>Ω</sup> is obtained by the Newton method. The following equations are used for correction of the pressure and velocity simultaneously. Here, the subscript i indicates grid location, while the superscripts m and n indicate the iteration of the corrections for the convergence of Eq. (56) and the time step, respectively. The more detailed explanation can be

rνΩ<sup>1</sup> ,

<sup>η</sup><sup>4</sup> <sup>1</sup> � <sup>μ</sup> � � <sup>1</sup> � <sup>4</sup><sup>μ</sup> � �

<sup>1</sup> � <sup>η</sup><sup>2</sup> ð Þ<sup>2</sup> , <sup>κ</sup> <sup>¼</sup> <sup>1</sup> � <sup>μ</sup>=η<sup>2</sup>

<sup>2</sup> , D <sup>¼</sup> <sup>d</sup>

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer

dR , D<sup>∗</sup> <sup>¼</sup> <sup>d</sup>

http://dx.doi.org/10.5772/intechopen.72263

dR þ 1 R (60)

257

(61)

<sup>1</sup> � <sup>μ</sup> (62)

$$\frac{\partial \overline{u'} }{\partial t} = -\frac{1}{\rho} \frac{\partial p'}{\partial r} + \nu \left( \frac{\partial^2 \underline{u'}}{\partial r^2} + \frac{1}{r} \frac{\partial \underline{u'}}{\partial r} - \frac{\underline{u'}}{r^2} + \frac{\partial^2 \underline{u'}}{\partial z^2} \right) + \frac{2 \overline{u}\_\theta \underline{v'}}{r} \tag{51}$$

$$\frac{\partial \overline{v}'}{\partial t} = \nu \left( \frac{\partial^2 \overline{v}'}{\partial r^2} + \frac{1}{r} \frac{\partial \overline{v}'}{\partial r} - \frac{\overline{v}'}{r^2} + \frac{\partial^2 \overline{v}'}{\partial z^2} \right) - \left( \frac{d \overline{u}\_\theta}{dr} + \frac{\overline{u}\_\theta}{r} \right) \underline{u}' \tag{52}$$

$$\frac{\partial w^{\prime}}{\partial t} = -\frac{1}{\rho} \frac{\partial p^{\prime}}{\partial z} + \nu \left( \frac{\partial^2 w^{\prime}}{\partial r^2} + \frac{1}{r} \frac{\partial w^{\prime}}{\partial r} + \frac{\partial^2 w^{\prime}}{\partial z^2} \right) \tag{53}$$

By considering the periodicity of the secondary flow which could be happened, each component of infinitesimal disturbance is assumed to be given in the following form. Here, a is the axial wavenumber (real number) and s is angular frequency (complex number)

$$\frac{\stackrel{\text{u}'}{\tilde{u}(r)}}{\tilde{u}(r)} = \frac{\stackrel{\text{v}'}{\tilde{v}(r)}}{\tilde{v}(r)} = \frac{\stackrel{\text{w}'}{\tilde{w}(r)}}{\tilde{p}(r)} = \frac{\stackrel{\text{p}'}{\tilde{p}(r)}}{\tilde{p}(r)} = \exp\left(iaz + st\right) \tag{54}$$

#### 4.3. Linear stability analysis

The dimensionless simultaneous ordinary equations are summarized as follows:

Basic velocity

$$\overline{\mathcal{U}}\mathcal{U}\_{\theta}(\mathcal{R}) = \frac{\mu - \eta^2}{1 - \eta^2} \mathcal{R} + \frac{\eta^2 \left(1 - \mu\right)}{1 - \eta^2} \frac{1}{\mathcal{R}} \tag{55}$$

Disturbance equations for amplitude functions

$$D\_\*\tilde{\mathcal{U}} + i\mathcal{k}\tilde{\mathcal{W}} = 0\tag{56}$$

$$
\Delta \tilde{s} \tilde{\boldsymbol{L}} = -D \tilde{\boldsymbol{P}} + \left( \boldsymbol{D} \boldsymbol{D}\_{\ast} - \boldsymbol{k}^{2} \right) \tilde{\boldsymbol{\mathcal{U}}} + \mathrm{Re}\_{\Omega} \frac{2 \mathcal{U}\_{\theta}}{R} \tilde{\boldsymbol{V}} \tag{57}
$$

$$
\tilde{SV} = \left(\mathbf{D}\mathbf{D}\_\* - k^2\right)\tilde{V} - \mathrm{Re}\_\Omega \left(\mathbf{D}\_\* \overline{\mathbf{U}}\_\theta\right)\tilde{\mathbf{U}}\tag{58}
$$

$$\mathcal{S}\tilde{\mathcal{W}} = -\mathrm{i}k\tilde{P} + \left(D\_{\ast}D - k^{2}\right)\tilde{\mathcal{W}}\tag{59}$$

Here, the dimensionless variables and non-dimensional numbers are as follows. The outer radius r<sup>2</sup> is taken as the characteristic length

$$\begin{aligned} R &= \frac{r}{r\_2}, \quad \overline{\mathcal{U}}\_\theta = \frac{\overline{u}\_\theta}{\Omega\_1 r\_2}, \quad \left( \check{\mathcal{U}}, \check{V}, \check{W} \right) = \frac{(\check{u}, \check{v}, \check{w})}{\Omega\_1 r\_2}, \quad \check{P} = \frac{\check{p}}{\rho v \Omega\_1}, \\\ Re\_\Omega &= \frac{\Omega\_1 r\_2^2}{\nu}, \quad k = r\_2 a, \quad \eta = \frac{r\_1}{r\_2}, \quad \mu = \frac{\Omega\_2}{\Omega\_1}, \quad S = \frac{s}{\nu/r\_2^2}, \quad D = \frac{d}{d\mathcal{R}}, \quad D\_\* = \frac{d}{d\mathcal{R}} + \frac{1}{R} \end{aligned} \tag{60}$$

The boundary conditions are as follows:

uθð Þ¼ r; z; t uθð Þþ r v

p rð Þ¼ ; z; t p rð Þþ ; z p

∂u 0 <sup>∂</sup><sup>t</sup> ¼ � <sup>1</sup> r ∂p 0 ∂r þ ν

∂v 0 <sup>∂</sup><sup>t</sup> <sup>¼</sup> <sup>ν</sup> 0

256 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

0

∂2 v 0 ∂r<sup>2</sup> þ 1 r ∂v 0 <sup>∂</sup><sup>r</sup> � <sup>v</sup> 0 r<sup>2</sup> þ ∂2 v 0 ∂z<sup>2</sup>

∂w<sup>0</sup> <sup>∂</sup><sup>t</sup> ¼ � <sup>1</sup> r ∂p 0 ∂z þ ν

u 0 u r <sup>~</sup>ð Þ <sup>¼</sup> <sup>v</sup>

Disturbance equations for amplitude functions

radius r<sup>2</sup> is taken as the characteristic length

4.3. Linear stability analysis

Basic velocity

ð Þ r; z; t , urð Þ¼ r; z; t u

∂u 0 ∂r þ u 0 r þ ∂w<sup>0</sup>

> ∂2 u 0 ∂r<sup>2</sup> þ 1 r ∂u 0 <sup>∂</sup><sup>r</sup> � <sup>u</sup> 0 r<sup>2</sup> þ ∂2 u 0 ∂z<sup>2</sup>

axial wavenumber (real number) and s is angular frequency (complex number)

The dimensionless simultaneous ordinary equations are summarized as follows:

<sup>U</sup>θð Þ¼ <sup>R</sup> <sup>μ</sup> � <sup>η</sup><sup>2</sup>

SU<sup>~</sup> ¼ �DP<sup>~</sup> <sup>þ</sup> DD<sup>∗</sup> � <sup>k</sup>

SV<sup>~</sup> <sup>¼</sup> DD<sup>∗</sup> � <sup>k</sup>

0 <sup>~</sup>v rð Þ <sup>¼</sup> <sup>w</sup><sup>0</sup>

After neglecting the second-order disturbance, the following linearized equations are obtained:

∂2 w0 ∂r<sup>2</sup> þ

By considering the periodicity of the secondary flow which could be happened, each component of infinitesimal disturbance is assumed to be given in the following form. Here, a is the

0

<sup>1</sup> � <sup>η</sup><sup>2</sup> <sup>R</sup> <sup>þ</sup> <sup>η</sup><sup>2</sup> <sup>1</sup> � <sup>μ</sup>

<sup>2</sup> <sup>U</sup><sup>~</sup> <sup>þ</sup> Re<sup>Ω</sup>

<sup>2</sup> <sup>V</sup><sup>~</sup> � Re<sup>Ω</sup> <sup>D</sup>∗U<sup>θ</sup>

Here, the dimensionless variables and non-dimensional numbers are as follows. The outer

1 � η<sup>2</sup>

1

<sup>D</sup>∗U<sup>~</sup> <sup>þ</sup> ikW<sup>~</sup> <sup>¼</sup> <sup>0</sup> (56)

2U<sup>θ</sup>

SW<sup>~</sup> ¼ �ikP<sup>~</sup> <sup>þ</sup> <sup>D</sup>∗<sup>D</sup> � <sup>k</sup><sup>2</sup> <sup>W</sup><sup>~</sup> (59)

w r <sup>~</sup> ð Þ <sup>¼</sup> <sup>p</sup>

1 r ∂w<sup>0</sup> ∂r þ ∂2 w0 ∂z<sup>2</sup>

0

ð Þ r; z; t , uzð Þ¼ r; z; t w<sup>0</sup>

þ 2uθv 0

uθ r 

u 0

<sup>~</sup>p rð Þ <sup>¼</sup> exp ð Þ iaz <sup>þ</sup> st (54)

<sup>R</sup> (55)

<sup>R</sup> <sup>V</sup><sup>~</sup> (57)

U~ (58)

� du<sup>θ</sup> dr þ

<sup>∂</sup><sup>z</sup> <sup>¼</sup> <sup>0</sup> (50)

ð Þ <sup>r</sup>; <sup>z</sup>; <sup>t</sup> (49)

ð Þ r; z; t ,

<sup>r</sup> (51)

(52)

(53)

$$\begin{cases} R = \eta: & \tilde{\mathcal{U}} = \tilde{V} = \tilde{W} = 0 \quad \text{(Inner wall)}\\ R = 1: & \tilde{\mathcal{U}} = \tilde{V} = \tilde{W} = 0 \quad \text{(Outer wall)}. \end{cases} \tag{61}$$

After Chandrasekar [6], the following two non-dimensional numbers are introduced to verify the computational results:

$$Ta = \frac{4\Omega\_1^2 r\_1^4}{\nu^2} \frac{(1-\mu)\left(1-4\mu\right)}{\left(1-\eta^2\right)^2} = 4\text{Re}\_\Omega^2 \frac{\eta^4 \left(1-\mu\right)\left(1-4\mu\right)}{\left(1-\eta^2\right)^2}, \quad \kappa = \frac{1-\mu/\eta^2}{1-\mu} \tag{62}$$

In this section, it is assumed that S = 0. This indicates that the secondary flow caused by the centrifugal instability is stationary and it contains toroidal vortices. To deal with the simultaneous ordinary differential equations for the boundary value problem, a one-dimensional staggered grid system is employed as shown in Figure 5. All the equations are discretized by the fourth order central difference method with a given wavenumber k using the HSMAC method [8] during which Re<sup>Ω</sup> is obtained by the Newton method. The following equations are used for correction of the pressure and velocity simultaneously. Here, the subscript i indicates grid location, while the superscripts m and n indicate the iteration of the corrections for the convergence of Eq. (56) and the time step, respectively. The more detailed explanation can be found in the recent papers published by the present author [9, 10]

Figure 5. The staggered grids in the radius direction together with the points of each variable definition.

$$\begin{split} \, ^{m+1}\hat{P}\_{i}^{n+1} &= ^{m}\hat{P}\_{i}^{n+1} + ^{m}(\delta\hat{P})\_{i}^{n+1} \\ &= ^{m}\hat{\boldsymbol{U}}\_{i}^{n+1} - \frac{\, ^{m}\hat{\boldsymbol{U}}\_{i+1}^{n+1} + 2\mathcal{T}^{m}\hat{\boldsymbol{U}}\_{i}^{n+1} - 2\mathcal{T}^{m}\hat{\boldsymbol{U}}\_{i-1}^{n+1} + \, ^{m}\hat{\boldsymbol{U}}\_{i-2}^{n+1}}{24(\Delta\boldsymbol{R})} + \frac{\, ^{m}\hat{\boldsymbol{U}}\_{i+1}^{n+1} + \theta^{m}\hat{\boldsymbol{U}}\_{i}^{n+1} + \theta^{m}\hat{\boldsymbol{U}}\_{i-1}^{n+1} - \, ^{m}\hat{\boldsymbol{U}}\_{i-2}^{n+1}}{16R\_{\rm R}} + \, \text{ik} \cdot ^{m}\hat{\boldsymbol{W}}\_{i}^{n+1} \end{split} \tag{63}$$

$$\mathbf{U}^{m+1} \tilde{\mathbf{U}}\_i^{n+1} = \prescript{m}{}{\mathbf{U}}\_i^{n+1} + \frac{\Delta \mathbf{r}}{\Delta \mathbf{R}} \cdot \prescript{m}{}{\left(\delta \tilde{\mathbf{P}}\right)^{n+1}\_i} \quad \overset{m+1}{\ }{\mathbf{U}}\_{i-1}^{n+1} = \prescript{m}{}{\tilde{\mathbf{U}}}\_{i-1}^{n+1} - \frac{\Delta \mathbf{r}}{\Delta \mathbf{R}} \cdot \prescript{m}{}{\left(\delta \tilde{\mathbf{P}}\right)^{n+1}\_i} \tag{64}$$

$$\boldsymbol{^{m+1}i\tilde{W}^{n+1}\_{i}} = {}^{m}\tilde{W}^{n+1}\_{i} - i\mathbf{k}(\Delta\boldsymbol{\pi}) \cdot {}^{m} \left(\delta\tilde{P}\right)^{n+1}\_{i} \tag{65}$$

Table 2 shows the computational results for various rotation speeds at η = 0.5. When μ > 0.25, the basic flow is always stable due to the Rayleigh's criterion. The present results exhibit slightly smaller values of Taylor number than those of Chandrasekar. Figures 6 and 7 show the amplitude functions and Eigen functions, respectively, for the case of μ = 0 (the outer cylinder is stationary), and Figures 8 and 9 show the case of μ = �0.5 (the outer cylinder rotates with half angular velocity in opposite direction to the inside rotation).

The simultaneous ordinary equations from (56) to (59) were divided into the real and imaginary parts. However, only four equations among the eight equations are necessary to solve in this problem because of the symmetricity and anti-symmetricity of the complex variables. In Figures 6 and 8, the real part of U, ~ V , ~ P~ and the imaginary part of W~ are shown. For the visualization shown in Figures 7 and 9, the Stokes stream function Ψ is defined as follows:


<sup>U</sup><sup>~</sup> <sup>ℜ</sup> � cos ð Þ¼ kZ <sup>1</sup>

the similar manner using the trigonometric functions.

function, azimuthal velocity, and pressure.

Figure 6. Amplitude functions (η = 0.5, μ = 0, k = 6.325).

R ∂Ψ

Here, the subscripts ℜ and ℑ represent the real part and the imaginary part, respectively. The visualization of other variables, such as the azimuthal velocity and the pressure, are treated in

Figure 7. Visualization of Eigen functions for two wavelengths (η = 0.5, μ = 0, k = 6.325). From left to right, Stokes stream

<sup>∂</sup><sup>Z</sup> , <sup>W</sup><sup>~</sup> <sup>ℑ</sup> � sin ð Þ¼ kZ <sup>1</sup>

R ∂Ψ

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer

http://dx.doi.org/10.5772/intechopen.72263

259

<sup>∂</sup><sup>R</sup> (66)

Table 2. Computational results and comparison with Chandrasekar (η = 0.5).

Figure 6. Amplitude functions (η = 0.5, μ = 0, k = 6.325).

mþ1 <sup>P</sup><sup>~</sup> <sup>n</sup>þ<sup>1</sup> <sup>i</sup> <sup>¼</sup> mP<sup>~</sup>

> <sup>¼</sup> mP<sup>~</sup> <sup>n</sup>þ<sup>1</sup> <sup>i</sup> �

nþ1 <sup>i</sup> <sup>þ</sup> <sup>m</sup>ðδP~<sup>Þ</sup>

mþ1 <sup>U</sup><sup>~</sup> <sup>n</sup>þ<sup>1</sup>

defined as follows:

�mU<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

nþ1 i

258 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

<sup>i</sup>þ<sup>1</sup> <sup>þ</sup> <sup>27</sup>mU<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

<sup>i</sup> <sup>¼</sup> mU<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

<sup>i</sup> þ

mþ1 <sup>W</sup><sup>~</sup> <sup>n</sup>þ<sup>1</sup>

variables. In Figures 6 and 8, the real part of U,

Table 2. Computational results and comparison with Chandrasekar (η = 0.5).

<sup>i</sup> � <sup>27</sup>mU<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

Δτ

with half angular velocity in opposite direction to the inside rotation).

<sup>Δ</sup><sup>R</sup> � <sup>m</sup> <sup>δ</sup>P<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

<sup>i</sup> <sup>¼</sup> mW<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

<sup>i</sup>�<sup>1</sup> <sup>þ</sup> mU<sup>~</sup> <sup>n</sup>þ<sup>1</sup> i�2 <sup>24</sup>ðΔR<sup>Þ</sup> <sup>þ</sup> �mU<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

Δτf2=ðΔRÞ

<sup>i</sup> , <sup>m</sup>þ<sup>1</sup>

Table 2 shows the computational results for various rotation speeds at η = 0.5. When μ > 0.25, the basic flow is always stable due to the Rayleigh's criterion. The present results exhibit slightly smaller values of Taylor number than those of Chandrasekar. Figures 6 and 7 show the amplitude functions and Eigen functions, respectively, for the case of μ = 0 (the outer cylinder is stationary), and Figures 8 and 9 show the case of μ = �0.5 (the outer cylinder rotates

The simultaneous ordinary equations from (56) to (59) were divided into the real and imaginary parts. However, only four equations among the eight equations are necessary to solve in this problem because of the symmetricity and anti-symmetricity of the complex

shown. For the visualization shown in Figures 7 and 9, the Stokes stream function Ψ is

κ μ Critical wave number Critical Ta number Wavenumber Ta number 0 1/4 6.286 15316 6.4 15332 0.4 1/6 6.293 19518 6.4 19542 0.6 2/17 6.299 22617 6.4 22644 1.0 0 6.325 33062 6.4 33100 4/3 �1/8 6.403 53210 6.4 53280 1.6 �1/4 6.715 98520 6.4 99072 1.8 �4/11 7.819 197715 7.8 199540 1.9 �9/21 8.733 288761 8.6 293630 2.0 �1/2 9.602 417734 9.6 428650

~ V ,

Present (201 grids) Chandrasekar [6]

<sup>i</sup>þ<sup>1</sup> <sup>þ</sup> <sup>9</sup>mU<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

<sup>i</sup>�<sup>1</sup> <sup>¼</sup> mU<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

<sup>2</sup> <sup>þ</sup> <sup>k</sup> 2 g

<sup>U</sup><sup>~</sup> <sup>n</sup>þ<sup>1</sup>

<sup>i</sup> � ikð Þ� Δτ <sup>m</sup> <sup>δ</sup>P<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

<sup>i</sup> <sup>þ</sup> <sup>9</sup>mU<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

16RPi

<sup>i</sup>�<sup>1</sup> � Δτ

<sup>i</sup>�<sup>1</sup> � mU<sup>~</sup> <sup>n</sup>þ<sup>1</sup> i�2

<sup>Δ</sup><sup>R</sup> � <sup>m</sup> <sup>δ</sup>P<sup>~</sup> <sup>n</sup>þ<sup>1</sup>

~ P~ and the imaginary part of W~ are

<sup>i</sup> (65)

<sup>þ</sup> ik � mW<sup>~</sup> <sup>n</sup>þ<sup>1</sup> i

<sup>i</sup> (64)

(63)

Figure 7. Visualization of Eigen functions for two wavelengths (η = 0.5, μ = 0, k = 6.325). From left to right, Stokes stream function, azimuthal velocity, and pressure.

$$
\tilde{M}\_{\mathfrak{R}} \cdot \cos(kZ) = \frac{1}{R} \frac{\partial \Psi}{\partial Z}, \quad \tilde{W}\_{\mathfrak{F}} \cdot \sin(kZ) = \frac{1}{R} \frac{\partial \Psi}{\partial R} \tag{66}
$$

Here, the subscripts ℜ and ℑ represent the real part and the imaginary part, respectively. The visualization of other variables, such as the azimuthal velocity and the pressure, are treated in the similar manner using the trigonometric functions.

References

Springer; 1999

zine. 1911;21:697-711

78(794):1680-1695

Japan; August 2014

Journal of Heat Transfer. 1973;95(4):540-541

Laboratory of University of California. 1975. p. LA-5852

Papers III

1961

[1] Schlichting H, Gersten K. Boundary Layer Theory. 8th Revised and Enlarged ed. Berlin:

Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer

http://dx.doi.org/10.5772/intechopen.72263

261

[2] Stokes GG. On the effect of the internal friction of fluids on the motion of pendulums. Transactions of the Cambridge Philosophical Society. 1856;9(Part II):8-106 or Collected

[3] Rayleigh L. On the motion of solid bodies through viscous liquids. Philosophical Maga-

[4] Churchill SW, Ozoe H. A correlation for laminar free convection from a vertical plate.

[5] Le Fevre EJ. Laminar free convection from a vertical plane surface. In: Proceedings of the 9th International Congress of Applied Mechanics; Brussels; Vol. 4; 1956. pp. 168-174 [6] Chandrasekar S. Hydrodynamic and Hydromagnetic Stability. Oxford: Clarendon Press;

[7] Koschmieder EL. Bénard cells and Taylor vortices. In: Cambridge Monographs on Mechanics and Applied Mathematics. United State of America: Cambridge University Press; 1993

[8] Hirt CW, Nichols BD, Romero NC. SOLA: A Numerical Solution Algorithm for Transient Fluid Flows. Los Alamos, New Mexico, United State of America: Los Alamos Scientific

[9] Tagawa T, Egashira R. Fluid flow of a liquid metal in a cylinder driven by a rotating magnetic field. Transactions of the Japan Society of Mechanical Engineers, Part B. 2012;

[10] Tagawa T. Numerical investigation of Bénard-Marangoni convection of paramagnetic liquid in annular layers. In: The 15th International Heat Transfer Conference; Kyoto,

Figure 8. Amplitude functions (η = 0.5, μ = 0.5, k = 9.602).

Figure 9. Visualization of Eigen functions for two wavelengths (η = 0.5, μ = 0.5, k = 9.602). From left to right, Stokes stream function, azimuthal velocity, and pressure.

### Author details

Toshio Tagawa

Address all correspondence to: tagawa-toshio@tmu.ac.jp

Department of Aerospace Engineering, Tokyo Metropolitan University, Japan
