Open access peer-reviewed chapter

# Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer

By Toshio Tagawa

Submitted: May 12th 2017Reviewed: November 7th 2017Published: December 20th 2017

DOI: 10.5772/intechopen.72263

## Abstract

The present chapter introduces incompressible Newtonian fluid flow and heat transfer by using the finite difference method. Since the solution of the Navier-Stokes equation is not simple because of its unsteady and multi-dimensional characteristic, the present chapter focuses on the simplified flows owing to the similarity or periodicity. As a first section, the first Stoke problem is considered numerically by introducing the finite difference method. In the second section, 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. In the third example, linear stability analysis for the onset of secondary flow during the Taylor-Couette flow is numerically treated using the HSMAC method.

### Keywords

• finite difference method
• similar solution
• boundary layer
• linear stability analysis
• HSMAC method

## 1. Introduction

The governing equation for the fluid flow is known as Navier-Stokes equation, which is however difficult to solve analytically; and therefore, a lot of numerical techniques have been proposed and developed. Nevertheless various complex flow phenomena such as turbulent flow, multi-phase flow, compressible flow, combustion, and phase change encountered in the fields of engineering would have still difficulties to circumvent even using both present computational resources and numerical techniques. The present chapter devotes not to elucidate such complex phenomena, but to introduce rather simplified fluid flow by using the finite difference method.

One focuses on incompressible flows, in which physical properties such as the viscosity, the thermal conductivity, the specific heat are constant and even the fluid density is not a thermodynamic variable. This simplified assumption makes the fluid flow phenomena much easier to 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 number approximation, while the latter one the Boussinesq approximation.

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 and also to contribute the augmentation of accuracy of the results.

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  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.

## 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  in his famous treatment of the pendulum. Since Lord Rayleigh  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

ut=ν2uy2E1

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:

y=0:u=uwy:u0E2

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

u=uwUη,η=y2νtE3

Then, the following ordinary differential equation can be obtained

d2Udη2+2ηdU=0E4

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

η=0:U=1η:U0E5

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

U=uuw=1erfη=12π0ηexpξ2E6

The velocity profile is shown in Figure 1.

### 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)

η=0:U=1η=4:U=0E7

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:

Ui+12Ui+Ui1Δη2+2ηiUi+1Ui12Δη=0,i=23N1

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

1ηiΔηαiUi12Ui+1+ηiΔηβiUi+1=0,i=23N1.E8

Here, ηi=i1Δηand αi and βi are coefficients determined by the number of grids. The boundary condition (7) is modified

η=0:U1=1η=4:UN=0E9

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:

2β2000α32β3000α42β4000α52β5000α62U2U3U4U5U6=α2U1000β6U7E10

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

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:

00000α300000α400000α500000α60U2U3U4U5U6+2000002000002000002000002U2U3U4U5U6+0β200000β300000β400000β500000U2U3U4U5U6=α2U1000βN1UNE11

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)

2000002000002000002000002DU2U3U4U5U6n+1Un+1=α2U1000βN1UNb00000α300000α400000α500000α60LU2U3U4U5U6nUn0β200000β300000β400000β500000UU2U3U4U5U6nUnE12

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

Un+1=D1bL+UUnE13

This is equivalent to the following equation:

Uin+1=12αiUi1n+βiUi+1n,i=23456E14

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

ux+vy=0E15

Momentum equation

uux+vuy=ν2uy2+TTE16

Energy equation

uTx+vTy=α2Ty2E17

Boundary equation

y=0:u=v=0,q=kT/yy:u0,TTE18

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=xxa,Y=yya,U=uua,V=vva,θ=TTbTaE19

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

UX+vaxayaua1VY=0
UUX+vaxayaua1VUY=νxaya2ua22UY2+TbTxaua23+Taxaua24θ
UθX+vaxayaua1VθY=αxaya2ua52θY2
Y=0:U=V=0,qya/kTa6=θ/YY:U=0,θ=TTb/Ta7

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

Putting  = 0, and one obtains Tb=T.Hence  becomes θ=0.

Putting  = 1, and one gets qyakTa=1Ta=qyak

Putting  = 1, and one gets αxaya2ua=1ya=αxaua1/2

Putting  = 1, vaxayaua=1va=yauaxa

Putting  = 1, Taxaua2=1ua=Taxa1/2=qyakxa1/2=qkαxaua1/2xa1/2

ua=qk2αxa31/5=qxa4kα22/5αxa=RaPr2/5αxa,Ra=qxa4kανE20
ya=αxaua1/2=αxaRaPr2/5αxa1/2=xa2RaPr2/51/2=xaRaPr1/5E21
Ta=qyak=qkxaRaPr1/5=qxakRaPr1/5E22
va=yauaxa=xaRaPr1/5RaPr2/5αxaxa=RaPr1/5αxaE23

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

UX+VY=0E24
UUX+VUY=Pr2UY2+θE25
UθX+VθY=2θY2E26
Y=0:U=V=0,θ/Y=1Y:U=0,θ=0E27

The dimensionless variables are summarized as follows:

X=xxa,Y=yxaRaPr1/5,U=uαxaRaPr2/5,V=vαxaRaPr1/5,θ=TTqxakRaPr1/5E28

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

dV=15ηdU3UE29

Momentum equation

UdVVdU+Prd2Udη2+θ=0E30

Energy equation

U5ηθV+d2θdη2=0E31

Boundary conditions

η=0:U=V=0,/=1η:U=θ=0E32

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

η=yxRaxPr15,U=uαxRaxPr25,V=vαxRaxPr15,θ=TTqxkRaxPr15,Rax=qx4ανk,Pr=ναE33

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

Nux=hxxk=qxTwTk=qxTaθwk=gβqα2k1/5kqx1/5Ta1qxθwk=gβqα2k1/5x4/5θw=qx4α2k1/51θw=RaxPr1/51θwE34
Nux=RaxPr15θη=01=RaxNuxPr15θη=01E35

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

NuxRaxPr14=θη=054E36

High Prandtl number

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

Continuity of mass

dV=15ηdU3UE37

Momentum equation

UdVVdU+Prd2Udη2+Prθ=0E38

Energy equation

U5ηθV+d2θdη2=0E39

Boundary conditions

η=0:U=V=0,/=1η:U=θ=0E40

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

η=yxRax15,U=uαxRax25,V=vαxRax15,θ=TTqxkRax15,Rax=qx4ανk,Pr=ναE41

NuxRax14=θη=054E42

### 3.4. Numerical results

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 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 shows the summary of the local Nusselt number for various Prandtl number cases together with the reference of Churchill and Ozoe for comparison . 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

η=0:dU/=V=0,/=1η:U=θ=0

and that for Pr → ∞

η=0:U=V=0,/=1η:dU/=θ=0.

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  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.

Pr00.010.1110100
NuxRax14N/A0.4564
0.456
0.5234
0.524
0.5495
0.550
0.5631
0.5627
NuxRaxPr140.7107
0.6922
0.6694
0.670
0.5970
0.597
0.4564
0.456
N/A

### Table 1.

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

## 4. Linear stability of Taylor-Couette flow

### 4.1. Governing equations

In the text book of Chandrasekar , 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  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:

urr+urr+uzz=0E43
urt+ururr+uzurzuθ2r=1ρpr+ν2urr2+1rurrurr2+2urz2E44
uθt+uruθr+uzuθz+uruθr=ν2uθr2+1ruθruθr2+2uθz2E45
uzt+uruzr+uzuzz=1ρpz+ν2uzr2+1ruzr+2uzz2gE46

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

u¯θr=r22Ω2r12Ω1r12r22r+r12r22Ω2Ω1r12r221rE47

Pressure

p¯rz=ρu¯θr2rdrρgz+p0E48

Here, Ω1 is the angular velocity at the inner cylinder, Ω2 is the angular velocity at the outer cylinder, p is the pressure, ρ 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:

uθrzt=u¯θr+v'rzt,urrzt=u'rzt,uzrzt=w'rzt,przt=p¯rz+p'rztE49

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

u'r+u'r+w'z=0E50
u't=1ρp'r+ν2u'r2+1ru'ru'r2+2u'z2+2u¯θv'rE51
v't=ν2v'r2+1rv'rv'r2+2v'z2du¯θdr+u¯θru'E52
w't=1ρp'z+ν2w'r2+1rw'r+2w'z2E53

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)

u'u˜r=v'v˜r=w'w˜r=p'p˜r=expiaz+stE54

### 4.3. Linear stability analysis

The dimensionless simultaneous ordinary equations are summarized as follows:

Basic velocity

U¯θR=μη21η2R+η21μ1η21RE55

Disturbance equations for amplitude functions

DU˜+ikW˜=0E56
SU˜=DP˜+DDk2U˜+ReΩ2U¯θRV˜E57
SV˜=DDk2V˜ReΩDU¯θU˜E58
SW˜=ikP˜+DDk2W˜E59

Here, the dimensionless variables and non-dimensional numbers are as follows. The outer radius r2is taken as the characteristic length

R=rr2,U¯θ=u¯θΩ1r2,U˜V˜W˜=u˜v˜w˜Ω1r2,P˜=p˜ρνΩ1,ReΩ=Ω1r22ν,k=r2a,η=r1r2,μ=Ω2Ω1,S=sν/r22,D=ddR,D=ddR+1RE60

The boundary conditions are as follows:

R=η:U˜=V˜=W˜=0InnerwallR=1:U˜=V˜=W˜=0OuterwallE61

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

Ta=4Ω12r14ν21μ14μ1η22=4ReΩ2η41μ14μ1η22,κ=1μ/η21μE62

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  during which ReΩ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]

P˜m+1in+1=P˜min+1+(δP˜)min+1=P˜min+1U˜mi+1n+1+27U˜min+127U˜mi1n+1+U˜mi2n+124(ΔR)+U˜mi+1n+1+9U˜min+1+9U˜mi1n+1U˜mi2n+116RPi+ikW˜min+1Δτ{2/(ΔR)2+k2}E63
U˜m+1in+1=U˜min+1+ΔτΔRδP˜min+1,U˜m+1i1n+1=U˜mi1n+1ΔτΔRδP˜min+1E64
W˜m+1in+1=W˜min+1ikΔτδP˜min+1E65 Figure 5.The staggered grids in the radius direction together with the points of each variable definition.

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).

Present (201 grids)Chandrasekar 
κμCritical wave numberCritical Ta numberWavenumberTa number
01/46.286153166.415332
0.41/66.293195186.419542
0.62/176.299226176.422644
1.006.325330626.433100
4/3−1/86.403532106.453280
1.6−1/46.715985206.499072
1.8−4/117.8191977157.8199540
1.9−9/218.7332887618.6293630
2.0−1/29.6024177349.6428650

### Table 2.

Computational results and comparison with Chandrasekar (η = 0.5). Figure 6.Amplitude functions (η = 0.5, μ = 0, k = 6.325). 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. 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.

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:

U˜coskZ=1R∂ΨZ,W˜sinkZ=1R∂ΨRE66

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.

chapter PDF
Citations in RIS format
Citations in bibtex format

## More

© 2017 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## How to cite and reference

### Cite this chapter Copy to clipboard

Toshio Tagawa (December 20th 2017). Numerical Analysis of the Incompressible Fluid Flow and Heat Transfer, Finite Element Method - Simulation, Numerical Analysis and Solution Techniques, Răzvan Păcurar, IntechOpen, DOI: 10.5772/intechopen.72263. Available from:

### Related Content

Next chapter

#### Numerical Simulation of Wave (Shock Profile) Propagation of the Kuramoto-Sivashinsky Equation Using an Adaptive Mesh Method

By Denson Muzadziwa, Stephen T. Sikwila and Stanford Shateyi

First chapter

#### Finite Element and Finite Difference Methods for Elliptic and Parabolic Differential Equations

By Aklilu T. G. Giorges

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.