Introduction to Finite Element Method

(Elemental contributions and formation of Global Matrix)

Prof. Gautam Biswas

Indian Institute of Technology Kanpur
 Finite Element Methods
Let us recall Eqn. (11) which we obtained via integration by parts:

(11)

Where the terms have been numbered as (i) to (v) for our reference

e
Finite Element Method

(31)
Finite Element Methods

(32)
 Finite Element Method
The convective heat loss term (iii) can be obtained as:

(33)

(34)
 Finite Element Method

3
The shape functions for the vertices of the elemental domain:

(35)
 Finite Element Method

(36)

together with H(e)2 and

We can account for the Convective contribution.

All the above mentioned elemental matrices and vectors have been
obtained from the integration of the expressions for N1, N2 and N3 , as
defined in the equation (11)
 Finite Element Method
All the above mentioned elemental matrices and vectors have been obtained from the
integration of the expressions for N1 , N2 and N3 , as defined in Eq. (11). Let us, for the
present, consider a problem with no radiating boundary. Thus, the nodal equation for
any of the nodes of the FEM can be constructed by appropriately summing up the
contributions from the (i) heat conduction term (ii) heat generation term (iii) convective
boundary terms and (ii) the heat flux boundary. If the point lies on a known temperature
boundary, the nodal equation is very simple, it takes the form Ti = Tb for the node i
where Tb is the prescribed boundary temperature.

To find the properties of the overall system, we must assemble all the individual
element equations, that is, to combine the matrix equations of each element in an
appropriate way such that the resulting matrix represents the behaviour of the entire
solution region of the problem. The boundary conditions must be incorporated after
the assemblage of the individual element contributions that is

{K }{T } ={G}
Where {K } is the global LHS matrix, which is the assemblage of the individual
element LHS matrices, {G} is the global load vector which is assemblage of the
individual element load vectors and {T } is the global unknown vector
 Finite Element Method
Global matrix equation in terms of nodal temperatures:

(37)

Where {K} is conduction matrix and Gi is the heat load vector.

The radiation heat transfer process can be modeled as:

(38)

(39)
 Finite Element Method

Although hr varies over the length of the line element L(e) due to
variation in the value of T* , it may not be worth the trouble if the
element sizes are not too large. In such situations:

(40)

where i and j are the nodes which make up the concerned radiation
boundary element.
 Finite Element Method
The Unsteady Problems:
Governing equation for the problem is:

(41)

(42)
 Finite Element Method
The resulting global matrix equation derived from equation
(42) takes the form:

(43)

Explicit formulation scheme:

{ Tn}K (44)

Implicit formulation scheme:

{ Tn}K+1 (45)

Semi- implicit scheme:

(46)
{Tn}K+1 {Tn}K
 Finite Element Method

For the implicit scheme, substituting for {Tn }k+1 from equation (43), we
get:

(47)

Solving the matrix equation (47), the temperature vector at the (k+1)th
time level can be obtained, knowing the kth level temperature vector.
 Finite Element Method

Elements of Higher Orders:

4-noded quadrilateral interpolation for temperature can be written as:

(48)

4 3

1 2
 Finite Element Method
Elements of Higher Orders::

(49)
 Finite Element Method

(50)

(51)

The size of the elemental matrices will also be different for higher-
order elements.
 Finite Element Method
Data Structure for Assembly
1. IEN array, the connectivity array, creates the connection between Local
node and Global node
IEN (1:Nen, 1: Nel)
Nen is the number of local nodes
Nel is the number of elements
2. ID array, the destination array, establishes the relation between Global
node number and the equation number
ID (1: Nnp)
Nnp is number of (Global) points
ID (A) = { P, P is any integer, denotes Eqn. number
{ 0, if A ∈ ng
Here ng is the Dirichlet boundary. The boundary value is specified. In
such cases, no computation is needed. There will be no equation
number.
 Finite Element Methods
Data Structure for Assembly

3. LM array, the location matrix array, establishes connection between

element number to equation number.
LM (1:Nen, 1: Nel)
LM array uses ID and IEN arrays. It determines the position, where
each element of elemental matrix will be placed in the Global matrix.

Global node
ID
IEN

Equation Element
LM
 Finite Element Methods

Example
Temp.
Specified 2
2 3 4 3

4 3 1
3
3 5 3
3
5 5 2 3
4 2 5

1
1 1
1 2
4 2
1 2 1 2 1

Temp.
Specified 4 5
denotes a local
3 node number
Here Nen=3, Nel=4, Nnp=5. denotes
1 2 element
number
 Finite Element Methods

Example

IEN Array: IEN (3,4)

ID Array: ID (5)
Nel
1 2 3 4
1 1 2 3 4 Nnp
2 2 3 4 1 1 2 3 4 5
3 5 5 5 5 0 0 0 0 1
Nen

LM Array: LM (3,4)
Nel
1 2 3 4
1 0 0 0 0
2 0 0 0 0
Nen 3 1 1 1 1
 Finite Element Method

Fundamentals of the Finite Element Method for Heat and Fluid Flow,
Roland Lewis, Perumal Nithiarasu and K.N. Seetharamu, John Wiley and
Sons, 2004, Chichester, England

Introduction to Finite Elements in Engineering, T. R. Chandrupatla and A.D.

Belegundu, Prentice Hall of India, 2002, New Delhi

Thank You
Chapter 5

Vorticity-Stream Function
Approach

5.1 Introduction
The vorticity-stream function method is one of the most popular methods for
solving 2-D incompressible Navier -Stokes equations. The governing equations
are
 
∂u ∂u ∂p
ρ u +v =− + µ ∇2 u (5.1)
∂x ∂y ∂x
 
∂v ∂v ∂p
ρ u +v =− + µ ∇2 v (5.2)
∂x ∂y ∂y

and
∂u ∂v
+ =0 (5.3)
∂x ∂y

As it was discussed earlier, it is difficult to deal with equation (5.1 and 5.2)
due to the lack of presence of a separate equation for pressure p. We introduce
stream function (ψ) and vorticity (ω) as

∂ψ ∂ψ
u= , v=− (5.4)
∂y ∂x

∂v ∂u
ω= − (5.5)
∂x ∂y

We can readily see that existence of (5.4) automatically satisfies continuity equa-
tion (5.3).

1
5.2 Computational Fluid Dynamics

5.2 Numerical Methodology

If we substitute the dependent variable with stream function, we shall not be
concerned with equation (5.3) any more. Invoking equation (5.4) into (5.5) we
obtain Poisson equation
∂2ψ ∂2ψ
+ = −ω (5.6)
∂x2 ∂y 2
Now we differentiate equation (5.1) with respect to y and equation (5.2) with
respect to x. If we subtract differentiated equation (5.2) from differentiated
equation (5.1) and rearrange the resulting equation, we shall obtain
 2
∂ ω ∂2ω

∂ω ∂ω
u +v =ν + (5.7)
∂x ∂y ∂x2 ∂y 2
This equation is the vorticity transport equation. Let us express the equation
∂2ψ ∂2ψ
(5.6) in terms of finite difference quotient for , and ωi, j to
∂x2 i, j ∂y 2 i, j
get
ψi+1, j − 2ψi, j + ψi−1, j ψi, j+1 − 2ψi, j + ψi, j−1
+ = −ωi, j
(∆x)2 (∆y)2

ψi+1, j + ψi−1, j + λ2 (ψi, − 2 1 + λ2 ψi, = −ωi, j (∆x)2

or, j+1 + ψi, j−1 ) j

2
ψi+1, j + ψi−1, j + (ψi, j+1 + ψi, j−1 ) λ + ωi, j (∆x)2
or, ψi, j =
2 (1 + λ2 )
where λ is the mesh aspect ratio (∆x)/(∆y). Putting ∆x = ∆y, we get
1
∆x2

ψi, ψi+1, j + ψi−1, j + ψi,
j = j+1 + ψi, j−1 + ωi, j (5.8)
4
Now for equation (5.7) we can write

ωi+1, j − ωi−1, j ωi, − ωi, j−1

j+1
ui, j + vi, j
2∆x 2∆y

ωi+1, j − 2ωi, j + ωi−1, ωi, j+1 − 2ωi, j + ωi,
 (5.9)
j j−1
=ν +
(∆x)2 (∆y)2
If we put ∆x = ∆y, we shall obtain
   
1 1 ui, j ∆x 1 ui, j ∆x
ωi, j = 1− ωi+1, j + 1 + ωi−1, j
4 2 ν 2 ν
    
1 vi, j ∆x 1 vi, j ∆x
+ 1− ωi, j+1 + 1 + ωi, j−1 (5.10)
2 ν 2 ν
And we rewrite equation (5.4) as
ψi, j+1 − ψi, j−1
ui, j = (5.11)
2∆y
 Vorticity-Stream Function Approach 5.3

ψi+1, j − ψi−1, j
vi, j =− (5.12)
2∆x
Thus we have now a system of simultaneous equation, (5.8), (5.10), (5.11) and
(5.12) to be solved for ψi, j , ωi, j , ui, j and vi, j . Let us discuss the solution
procedure

Figure 5.1: Two-dimensional Computational Domain (Flow over a

backward facing step)

5.3 The Algorithm

1. Divide the physical domain by a mesh system where ∆x = ∆y

2. Set the known boundary conditions for u, v, ψ and ω.

3. Choose initial values of ψi, j and ωi, j at the interior grid points. Taking
initial values for vorticity as zero is usually acceptable. Initial values for
ψi, j at each i column of points can be calculated from the axial velocity
profile at that location. However, ui, j for interior points may be taken as
U∞ and vi, j = 0.

4. Calculate everywhere ψi, j using equation (5.8). Gauss-Seidel or overre-

laxation type calculation is done, for example:
′ ′  ′ 
ψi,k+1
j = ψ k
i, j + F ψ k+1
i, j − ψ k
i, j

′
where, ψi,k j is the value from previous calculation, ψi,k+1
j is the most recent
′
value, F is the overrelaxation factor and ψi,k+1
j is newly adjusted better
guess.
5.4 Computational Fluid Dynamics

5. Calculate ui, j and vi, j at all the internal grid points using equation (5.11
and 5.12).
6. In a subsequent step, calculate ωi, j at all interior mesh points using equa-
tion (5.10).
7. Apply appropriate boundary conditions (which has been discussed, in de-
tails, in a following section).
8. Go back to step (4) and calculate ψi,k+1j with the help of current value of
k
ωi, j and ψi, j values, where k denotes the previous level. After, evaluating
′
ψi,k+1 k+1
j at all points, find out ψi, j at all points which are indeed improved
values. Start repeating steps 4 to 8 untill the desired degree of convergence
is achieved.

5.3.1 Bottom Wall Boundary

Now let us discuss about the boundary conditions. Consider Fig. 5.1, we
shall call B1 and B3 as bottom wall. Similar kind of boundary conditions are
applicable on B1 and B3. At the nodal points which are coinciding with the
solid wall we can directly put ui, j = 0 and vi, j = 0. Since the line B1−B2−B3
is a streamline, any constant value of ψ on it is acceptable. The usual choice is
ψi, j = 0. The wall vorticity is an extremely important evaluation. At no- slip
boundaries, ωi, j is produced. It is the diffusion and subsequent advection of
the wall produced vorticity which governs the physics. Using boundary B1 as
an example, we expand ψi, js+1 by a Taylor series as
∂ψ 1 ∂2ψ 1 ∂3ψ
ψi, js+1 = ψi, js + ∆y + (∆y)2 + (∆y)3 + .....
∂y i, js 2 ∂y 2 i, js 6 ∂y 3 i, js
(5.13)
∂ψ ∂2ψ ∂u
But = ui, js = 0 by no-slip condition and =
∂y i, js ∂y 2 i, js ∂y i, js
∂v ∂u
Again, ω = −
∂x ∂y
∂v
Along the wall, = 0 [because v = constant = 0].
∂x i, js
∂u
Thus, ωi, js = −
∂y i, js
Substituting this into (5.13) and solving for ωi, js with ψi, js = 0 gives
−2ψi, js+1
ωi, js =
(∆y)2
More general form, regardless of the wall orientation or value of ψ at the bound-
ary, it can be written as
−2(ψw+1 − ψw )
ω|at w = (5.14)
(∆n)2
 Vorticity-Stream Function Approach 5.5

where ∆n is the distance from (w + 1) to (w) in the normal direction [w denotes

at the wall].

5.3.2 Upper Boundary

The upper boundary B5 in Fig 5.1 is having the usual no-slip and impervious
conditions for velocity components , i.e., ui, j = 0, vi, j = 0 . For vorticity
(ωi, j ), (5.14) will apply.
But how to evaluate ψi, j at the upper wall?
The value of ψi, j at the upper wall is constant and may be evaluated by inte-
grating the u velocity profile at the inlet. Integration may be performed through
Simpson’s rule to get
Z y at j=JMAX
ψi, j |w = u(y)inlet dy (5.15)
y at j= js

If we want to model the condition of no boundary at B5, or, in other words, in y-

direction, fluid at infinite extent is assumed, the problem is little more difficult.
However, Thoman and Szewczyk (1966) used a treatment which specifies this
∂ψ ∂ψ
far-field condition of ωi, j = 0 with = u = U∞ and = v = 0. Thus
∂y ∂x
u = U∞ was applied through a Neumann condition at the boundary along B5
as
ψi, ju = ψi, ju−1 + U∞ ∆y (5.16)
where B5 is considered at j = ju

5.3.3 Inlet Boundary

Inlet boundary in Fig. 5.1 cannot have any unique prescription. It will depend
on the physical situation. For the axial velocity (u), uniform or parabolic or
any possible profile can be taken. Most widely used conditions are:

u1, j = U∞ for j = js to JM AX (5.17)

or, "  2 #
ym − y
u1, j = 1.5 Uav 1 − (5.18)
ym
j=js to JMAX

For normal velocity (v), Fromm and Harlow (1963)set

v1, j =0 for j = js to JM AX (5.19)

The stream function ψ1, j can be obtained from the axial velocity profile at the
inlet as Z y at j
ψ1, j = u(y) dy (5.20)
y at js
inlet
5.6 Computational Fluid Dynamics

Vorticity (ωi, j ) also depends on inlet velocity profile. Pao and Daugherty
(1969) used uniform axial velocity profile, v = 0, and then specified ω1, j =
0. Greenspan (1969) fixed up ψ1, j from axial velocity profile and assumed
∂v/∂x = 0, which results in

∂2ψ ∂u
ω1, j =− 2
=− (5.21)
∂y ∂y

5.3.4 Outflow Boundary

B4 is the outflow boundary (Fig. 5.1). If the outflow boundary conditions are
known beforehand, then why are we computing? They are not known explicitly,
but we can prescribe or set some gradients at the outlet which are physically
meaningful. We can imagine about continuative outflow conditions which will
ensure smooth transition through the outlet boundary. For axial and normal
velocities, we can impose less restrictive type conditions, which are

∂u ∂v
= =0 (5.22)
∂x ∂x
Thoman and Szewczyk (1966) developed outflow boundary conditions through
setting
∂ω
=0 (5.23)
∂x
∂v ∂2ψ
Then, from = 0, they derived = 0. For constant ∆x, at i = IM AX,
∂x ∂x2
this gives:
ωIMAX, j = ωIMAX−1, j (5.24)
ψIMAX, j = 2 ψIMAX−1, j − ψIMAX−2, j (5.25)

5.4 Some Important Observations

Some difficulties were experienced if the coefficients of ω’s on the right hand
side of equation (5.10) have a sum greater than unity. In that case, the process
u∆x v∆x
will diverge. The quantities, such as and play a vital role in such
ν ν
cases. As a matter of fact, these quantities are mesh Reynolds number. By
taking very small grid size, they may be kept below a desired small value of 2.
Another remedy which could be applied along with, is to introduce a upwind
bias. So the difference quotients with respect to convective components become
 
∂ω ωi, j − ωi−1, j
u = ui, j for ui, j > 0 (5.26)
∂x ∆x
 
ωi+1, j − ωi, j
= ui, j for ui, j < 0 (5.27)
∆x
 Vorticity-Stream Function Approach 5.7

∂ω
In a similar way, v can be evaluated. These are then used in equation (5.9)
∂y
for the subsequent development of equation (5.10).

Regarding the vorticity boundary condition given by equation (5.14), it can

be said that the condition is first order accurate. The accuracy can be increased
and we can try for a second order accurate boundary condition.
∂ψ
Let us consider equation (5.13) once again. is zero along the solid boundary
∂y
∂2ψ
and ωi, j =− . If we want to increase the accuracy, we have to retain the
∂y 2
fourth term on the right hand side of equation (5.13). Now we shall try to
evaluate the fourth term. We know from the definition of vorticity:
∂ 2v ∂2u ∂3ψ
   
∂ω ∂ ∂v ∂ψ
= − 2 = − because u =
∂y ∂x∂y ∂y ∂x ∂y ∂y 3 ∂y
or,
∂ 2u ∂ 3 ψ
 
∂ω ∂v ∂u
=− 2 − from continuity, =−
∂y ∂x ∂y 3 ∂y ∂x
or,
∂3ψ
 2 
∂ω ∂ u
=− = 0, at the wall (5.28)
∂y at w ∂y 3 at w ∂x2
∂ψ ∂2ψ
Invoking (5.28) and substituting the values for and in equation (5.13),
∂y ∂y 2
we get
∆y 2 ∆y 3 ∂ω
ψw+1 = ψw − ω|at w − (5.29)
2 6 ∂y at w
Instead of ∆y we shall write ∆n to mean normal direction from the wall. We
∂ω
shall also substitute by a forward difference quotient as
∂y
(∆n)2 (∆n)3 ω|at w+1 − ω|at w
 
ψw+1 = ψw − ω|at w − x
2 6 (∆n)
or,
(∆n)2 (∆n)2
ψw+1 − ψw + ω|at w+1 = − ω|at w
6 3
or,
−3(ψw+1 − ψw ) 1
ω|at w = − ω|at w+1 (5.30)
(∆n)2 2
Equation (5.30) is second order accurate boundary condition but many a times
it does not lead to stable computation for high Reynolds number flows.
5.8 Computational Fluid Dynamics

5.5 Calculation of Pressure

Upon completing the steps enumerated earlier, the velocity components are
determined at each grid point. In order to determine the pressure at each point,
it is necessary to derive Poission equation for pressure. This equation is derived
by differenciating equation (5.1) by x and equation (5.2) by y and then adding
them
 2  2  2
∂2v
   
∂u ∂v ∂v ∂u ∂u
+ +2 +u +
∂x ∂y ∂x ∂y ∂x2 ∂x∂y
 2
∂2v
  
∂ u 1 2 ∂ ∂
∇2 u + ∇2 v



+v + 2 =− ∇ p+ν (5.31)
∂x∂y ∂y ρ ∂x ∂y

By making use of continuity equation (5.31) can be reduced to

 
2 ∂u ∂v ∂u ∂v
∇ p = 2ρ −
∂x ∂y ∂y ∂x
or,

" 2 #
∂2ψ ∂2ψ ∂2ψ
  
2
∇ p = 2ρ − (5.32)
∂x2 ∂y 2 ∂x∂y

Thus we have obtained a Poisson equation for pressure which is analogous to

equation(5.6). The Poisson equation for pressure can be solved over the entire
domain using appropriate conditions at the confining surfaces.

5.6 Application to Curvilinear Geometries

Let us consider steady, incompressible flow over a long cylinder of radius R0 .
The stream function and vorticity transport equations for this flow geometry
are written using cylindrical polar coordinates (r, θ) system as:
 2
1 ∂2ω

∂ω vθ ∂ω ∂ ω 1 ∂ω
vr + =ν + + 2 (5.33)
∂r r ∂θ ∂r2 r ∂r r ∂θ2

The vorticity is given by

1 ∂ 1 ∂vr
ω= (rvθ ) −
r ∂r r ∂θ
The radial and tangential velocity components are defined as

1 ∂ψ ∂ψ
vr = and vθ = − (5.34)
r ∂θ ∂r
5.8 Computational Fluid Dynamics

5.5 Calculation of Pressure

Upon completing the steps enumerated earlier, the velocity components are
determined at each grid point. In order to determine the pressure at each point,
it is necessary to derive Poission equation for pressure. This equation is derived
by differenciating equation (5.1) by x and equation (5.2) by y and then adding
them
 2  2  2
∂2v
   
∂u ∂v ∂v ∂u ∂u
+ +2 +u +
∂x ∂y ∂x ∂y ∂x2 ∂x∂y
 2
∂2v
  
∂ u 1 2 ∂ ∂
∇2 u + ∇2 v



+v + 2 =− ∇ p+ν (5.31)
∂x∂y ∂y ρ ∂x ∂y

By making use of continuity equation (5.31) can be reduced to

 
2 ∂u ∂v ∂u ∂v
∇ p = 2ρ −
∂x ∂y ∂y ∂x
or,

" 2 #
∂2ψ ∂2ψ ∂2ψ
  
2
∇ p = 2ρ − (5.32)
∂x2 ∂y 2 ∂x∂y

Thus we have obtained a Poisson equation for pressure which is analogous to

equation(5.6). The Poisson equation for pressure can be solved over the entire
domain using appropriate conditions at the confining surfaces.

5.6 Application to Curvilinear Geometries

Let us consider steady, incompressible flow over a long cylinder of radius R0 .
The stream function and vorticity transport equations for this flow geometry
are written using cylindrical polar coordinates (r, θ) system as:
 2
1 ∂2ω

∂ω vθ ∂ω ∂ ω 1 ∂ω
vr + =ν + + 2 (5.33)
∂r r ∂θ ∂r2 r ∂r r ∂θ2

The vorticity is given by

1 ∂ 1 ∂vr
ω= (rvθ ) −
r ∂r r ∂θ
The radial and tangential velocity components are defined as

1 ∂ψ ∂ψ
vr = and vθ = − (5.34)
r ∂θ ∂r
 Vorticity-Stream Function Approach 5.9

Figure 5.2: Finite Difference Grid for Flow Over a Cylinder

Substituting velocity components in vorticity, we get

∂ 2 ψ 1 ∂ψ 1 ∂2ψ
 
+ + 2 = −ω (5.35)
∂r2 r ∂r r ∂θ2

The boundary conditions for stream function and vorticity have to be derived.
The stream function is a constant on the axis (θ = 0 and θ = π). This follows
vθ = 0 by the virtue of axisymmetry. On the solid surface, at r = R0 , a
constant stream line is obtainable following the normal velocity vr = 0 on the
cylinder.This flow configuration demands that the axis and the cylinder surface
are the parts of the same streamline. The value of the stream function can be
set equal to zero on these surfaces (Figure 5.2). Therefore, we can write

ψ = 0 on θ = 0, r ≥ R0
ψ = 0 on θ = π, r ≥ R0
ψ = 0 on r = R0 , 0 ≤ θ ≤ π (5.36)

In the far field, velocity becomes uniform entailing U0 as r → ∞. In terms of

stream function, this condition is expressed as:

r → ∞ ψ → U0 r sin θ (5.37)

Usually, a large cylindrical surface (r = R∞ ) is considered for implementing

the far field condition (Figure 5.2). Therefore, the far field boundary condition
(Equation (5.37)) is approximated as:

ψ = U0 R∞ sin θ at r = R∞ (5.38)
5.10 Computational Fluid Dynamics

The far field boundary condition for the vorticity is given by

ω = 0 at r = R∞ (5.39)

The velocity is uniform and velocity gradients are zero, the above-mentioned
condition if justified. On the axis (θ = 0 and π), ∂vr /∂θ = 0 due to axisymmetry
and ψ = 0 for all r.
Thus, from Equations (5.34) and (5.35) it can be seen that

ω = 0 on θ = 0, r ≥ R0
ω = 0 on θ = π, r ≥ R0 (5.40)

At r = R0 , the vorticity is not zero. Here, the wall vorticity has to be derived in
the same way as in the backward facing step problem. Equation (5.35) together
with the impervious and no-slip conditions of vr = 0 and vθ = 0 (or ψ = 0 and
∂ψ/∂r = 0 for all θ), yields
∂2ψ
= −ω (5.41)
∂r2
Now, the wall vorticity can be calculated using the same strategy as (5.14).
Considering a cylindrical-polar grid as shown in Figure 5.2, the derivatives of
stream function as well as vorticity can be replaced by their central difference
forms as:

∂ψ ψi+1, j − ψi−1, j ∂ψ ψi, j+1 − ψi, j−1

= ; =
∂θ 2∆θ ∂r 2∆r
(5.42)
∂ 2ψ ψi+1, j − 2ψi, j + ψi−1, j ∂ 2 ψ ψi, j+1 − 2ψi, j + ψi, j−1
= ; =
∂θ2 ∆θ2 ∂r2 ∆r2

Similarly, expressions for the derivatives of vorticity can also be obtained. Sub-
stituting these expressions into the governing equations (5.33) and (5.34) the
nodal equations of all interior nodes (i = 2, ... imax − 1, j = 2, ...jmax − 1) can
be derived. Finally, the set of discretized equations and boundary conditions are
solved by iterative methods. An upwinding strategy for the convective terms
can be implemented based on the magnitudes of vr and vθ , in a similar manner
as described earlier.

The important fact to be kept in mind while simulating flows in cylindrical-

polar geometry is that due to the involvement of (1/r) term, certain terms in
0
the governing equations or boundary conditions may take a form at r = 0.
0
0
In such cases, one can resolve the form using the L’ Hospital rule. More de-
0
tailed treatment on stream function vorticity formulation is available elsewhere
(Sundararajan, 2003).
 Vorticity-Stream Function Approach 5.11

References
1. Fromm, J.E. and Harlow, F.H., Numerical Solution of the Problem of
Vortex Street Development, The Physics of Fluids, Vol. 6, pp. 975-982,
1963.
2. Greenspan, D., Numerical Studies of Prototype Cavity Flow Problems,
The Computer J., Vol. 12, pp. 88-93, 1969.
3. Pao, Y.H. and Daugherty, R., Time Dependent Viscous Incompressible
Flow Past a Finite Flat Plate, Boeing Scientific Research Laboratories,
D1 - 82 - 0822, 1969.
4. Sundararajan, T., Solution of Viscous Incompressible Flows by the Stream
Function- Vorticity Formulation, in Computational Fluid Flow and Heat
Transfer, 2nd Edition, pp. 209- 227 ed. K Muralidhar and T. Sundarara-
jan,Narosa Publishing House, India, 2003.
5. Thoman, D.C. and Szewczyk, A. A., Numerical Solutions of Time De-
pendent Two-dimensional Flow of a Viscous Incompressible Fluid Over
Stationary and Rotating Cylinders, Tech. Rept. 66-14, Heat Transfer and
Fluid Mechanics Lab., University of Notre Dame, Indiana, USA, 1966.

Problems
1. Suppose you are modeling laminar flow by Stream-function-vorticity method
in a method in a 2D plane channel as shown below. Write the boundary
conditions for u, v, ψ and ω on AB, BC, CD and AD sides (please follow
the order in which they appear)

Figure 5.3: Symmetric flow in a channel

2. Drive first order and second order accurate vorticity boundary conditions
on a solid surface.