one-up PDF - Portland State University

The Control-Volume Finite-Difference
Approximation to the Diffusion Equation
ME 448/548 Notes
Gerald Recktenwald
Portland State University
Department of Mechanical & Materials Engineering
[email protected]
21 January 2014
ME 448/548: 2D Diffusion Equation
Motivation
Numerical solution of the 2D Poisson equation is the next step in developing our
knowledge of CFD technique.
• Introduce the Finite Volume Method
. Naturally deal with material discontinuity
. Can naturally enforce conservation of mass and energy
. Core idea in many (not all) commercial CFD codes
• Extend analysis to two spatial dimensions
. More interesting practical applications
. More complex data structures
. More complex procedures to solve the Ax = b problem
ME 448/548: 2D Diffusion Equation
1
Overview
Goals for this unit
• Introduce the Poisson equation
. A model of steady heat conduction with a source term
. Form of the pressure equation for incompressible flow
. Precursor to the generalized advection diffusion equation
•
•
•
•
•
•
Use the finite volume method to obtain discrete equations
Allow for non-uniform mesh, diffusion coefficient, and source term.
Introduce a set of Matlab codes for 2D Control Volume Finite Difference (CVFD)
Demonstrate solutions for three model problems.
Measure truncation error for a model problem with a simple solution
Apply to fully-developed flow in rectangular ducts
ME 448/548: 2D Diffusion Equation
2
Model Problem
The two dimensional diffusion equation in Cartesian coordinates is
∂
∂φ
∂
∂φ
Γ
+
Γ
+S =0
∂x
∂x
∂y
∂y
(1)
where φ is the scalar field, Γ is the diffusion coefficient, and S is the source term.
Example: Heat conduction in a rectangular domain
Γ = k, thermal conductivity or Γ =
k
=α
ρcp
S = volumetric heat source
ME 448/548: 2D Diffusion Equation
3
2D Cartesian Finite Volume Mesh
xu(1) = 0
x(i)
xu(nx+1)
xu(i)
yv(ny+1)
ny
nx ny
y(j)
3
yv(j)
Interior node
Boundary node
2
1
y
j=0
i=0
nx+1
nx+2
1
2
1
2
Ambiguous corner node
nx
3
yv(1) = 0
nx
x
ME 448/548: 2D Diffusion Equation
4
Finite Volume Mesh: Compass Point Notation
Nodes are labelled with their position
relative to the typical point “P”.
N, S, E, W are the names of north,
south, east and west neighbors.
N
W
E
y
x
S
Capital letters designate node locations.
φN , φS , φE , φW are the values of φ at
the north, south, east and west neighbor
nodes.
Use of compass point names simplifies
the algebra. For example the value of φ
at point N is φN or φi,j+1, but φN is
simpler to write.
Compass point notation is an historic
convention that does not extend to
unstructured meshes.
ME 448/548: 2D Diffusion Equation
5
Finite Volume Mesh
Detailed notation for a typical 2D control
volume
xw
∆x
xe and xw are the x coordinates of the
east and west faces of the control
volume around P
xe
N
δyn
yn
P
W
Lowercase letters designate the interface
between the nodes.
E
yn and ys are the y coordinates of the
north and south faces of the control
volume around P
∆y
δys
ys
δxw
δxe
Similarly, φn, φs, φe, φw are the values
of φ at the north, south, east and west
control volume interfaces.
S
ME 448/548: 2D Diffusion Equation
6
Finite Volume Mesh
Detailed notation for a typical 2D control
volume
xw
∆x
Regardless of whether the control
volumes sizes are uniform, node P is
always located in the geometric center
of each control volume.
xe
Thus,
N
δyn
xP − xw
=
yP − ys
=
yn
P
W
E
∆y
∆x
2
∆y
yn − yP =
2
xe − xP =
for uniform or non-uniform meshes
δys
ys
δxw
δxe
S
ME 448/548: 2D Diffusion Equation
7
Finite Volume Mesh
Detailed notation for a typical 2D control
volume
In contrast, the distances between the
nodes is not assumed to be uniform
δxe = xE − xP 6= δxw = xP − xW
xw
∆x
xe
δyn = yN − yP 6= δys = yP − yS
N
δyn
yn
P
W
E
δxe = xE − xP
6=
δxw = xP − xW
δyn = yN − yP
6=
δys = yP − yS
∆y
δys
ys
δxw
δxe
Note the use of upper and lower case
subscripts: Upper case refers to nodes;
lower case refers to interfaces.
S
ME 448/548: 2D Diffusion Equation
8
Convert the Differential Equation to a Discrete Equation
Integrate over the control volume
xw
Z
yn
ys
Z
xe
xw
∂
∂x
Z
∂φ
Γ
∂x
yn
=
ys
∆x
xe
dx dy
∂φ
Γ
∂x
N
δyn
−
e
∂φ
Γ
∂x
yn
P
W
E
∆y
dy
w
∂φ
∂φ
Γ
≈
− Γ
∆y
∂x e
∂x w
δys
ys
δxw
δxe
S
φE − φP
φP − φW
∆y
≈ Γe
− Γw
δxe
δxw
The final step is obtained by using central difference approximations for the derivatives at
the interfaces.
ME 448/548: 2D Diffusion Equation
9
Convert the Differential Equation to a Discrete Equation
Integrate the source term
Z
xe
Z
yn
S dy dx ≈ SP ∆x ∆y
xw
(2)
ys
Note: The Control Volume Finite Difference (CVFD) method treats the source term and
diffusion coefficients as piecewise constants. This is a rather crude approximation, say,
compared to allowing the source term and diffusion coefficient to vary linearly within the
control volume.
However, piecewise constant profiles of S and Γ allow the method to be conservative,
i.e., conserving mass or energy, automatically. The conservative nature of the CVFD
method is one of its primary strengths.
ME 448/548: 2D Diffusion Equation
10
Convert the Differential Equation to a Discrete Equation
Putting pieces back into the model equation gives the discrete system of equations
−aS φS − aW φW + aP φP − aE φE − aN φN = b
(3)
where
aE =
Γe
,
∆x δxe
aW =
Γw
,
∆x δxw
aN =
Γn
,
∆y δy n
aP = aE + aW + aN + aS
b = SP
aS =
Γs
∆y δy s
(4)
(5)
This is not a tridiagonal system of equations.
ME 448/548: 2D Diffusion Equation
11
Non-uniform Γ
Continuity of fluxes at the interface requires
∂φ ∂φ ∂φ ΓP
= ΓE
= Γe
∂x xe−
∂x xe+
∂x xe
material 1
material 2
E
P
Use central difference approximations
Γe
Γe
φE − φP
δxe
=
φE − φP
δxe
=
ME 448/548: 2D Diffusion Equation
ΓP
ΓE
φe − φP
δxe−
(6)
φE − φe
δxe+
(7)
δxe-
δxe+
δxe
12
Non-uniform Γ
Equations 6 and 7 can be rearranged as
material 1
φe − φP
=
φE − φe
=
δxe− Γe
(φE − φP )
ΓP δxe
δxe+ Γe
(φE − φP )
ΓE δxe
material 2
(8)
E
P
(9)
Add Equation 8 and Equation 9
δxe-
Γe
δxe−
δxe+
φE − φP =
(φE − φP )
+
.
δxe
ΓP
ΓE
δxe+
δxe
Cancel the factor of (φE − φP ) and solve for Γe/δxe to get
Γe
=
δxe
ME 448/548: 2D Diffusion Equation
δxe+
δxe−
+
ΓP
ΓE
−1
=
ΓE ΓP
.
δxe−ΓE + δxe+ΓP
13
Non-uniform Γ
Thus, the diffusion coefficient at the interface that results in flux continuity is
Γe =
ΓE ΓP
βΓE + (1 − β)ΓP
(10)
β≡
xe − xP
δxe−
=
δxe
xE − xP
(11)
where
An analogous derivation gives formulas for Γw , Γn, and Γs.
ME 448/548: 2D Diffusion Equation
14
Solving the System of Equations
Regardless of uniform or variable Γ, the discrete equation has a five-point stencil, and the
discrete equation for any interior node can be written.
−aS φS − aW φW + aP φP − aE φE − aN φN = b
(12)
To set up the matrix for this system of equations, we need to re-number the unknowns.
Important: i and j subscripts for the mesh are not the same as the row and
column indices in the system Ax = b.
ME 448/548: 2D Diffusion Equation
15
Solving the System of Equations
Order the nodes:
xu(1) = 0
x(i)
xu(nx+1)
xu(i)
n = i + (j − 1)nx
yv(ny+1)
ny
n is the node number (row number)
in Ax = b. i and j are the mesh
indices corresponding to xi and yj .
nx ny
y(j)
3
yv(j)
Interior node
Boundary node
With natural ordering the neighbors
in the compass point notation have
these indices
np = i + (j-1)*nx
ne = np + 1
nw = np - 1
nn = np + nx
ns = np - nx
ME 448/548: 2D Diffusion Equation
2
1
y
j=0
i=0
nx+1
nx+2
1
2
1
2
Ambiguous corner node
nx
3
yv(1) = 0
nx
x
16
Solving the System of Equations
This leads to a vector of unknowns

φ1,1
 φ2,1

...


 φ
nx,1


 φ1,2

 φ2,2

...


 φi,j

...

φnx,ny
ME 448/548: 2D Diffusion Equation


































⇐⇒
φ1
φ2
...
φnx
φnx+1
φnx+2
...
φn
...
φN

















(13)
17
Solving the System of Equations
The coefficient matrix has 5 non-zero diagonals
A=
ME 448/548: 2D Diffusion Equation
aS
a W ap
aE
aN
18
Algorithm for obtaining the numerical solution
1. Define physical parameters: Lx, Ly , Γ(x, y) and boundary conditions
2. Define the mesh: nx and ny if uniform
3. Compute the coefficient matrix
4. Solve the system of equations Aφ = b, where φ is the vector of unknowns.
5. Post-process to visualize the solution
The Poisson equation is steady. Each step is performed only once.
ME 448/548: 2D Diffusion Equation
19
Structured Mesh
The CVFD Matlab codes use structured meshes. in size.
•
•
•
•
•
The cells in the domain are topologically equivalent to a rectangular array
The cells need not be uniform in size.
Each cell not on a boundary touches four other cells
Each row has the same number of cells
Each column has the same number of cells
Uniform Mesh
Block-Uniform Mesh
Ly3, ny3
∆y
Ly2, ny2
Ly1, ny1
y
∆x
Lx1, nx1
Lx2, nx2
x
ME 448/548: 2D Diffusion Equation
20
Boundary Conditions (part 1)
Boundary
type
Boundary
Condition
Post-processing in fvpost
1
Specified T
Compute q 00 from discrete approximation to Fourier’s law.
T − Ti
00
q =k b
xb − xi
where Ti and Tb are interior and boundary temperatures, respectively.
2
Specified q 00
Compute Tb from discrete approximation to Fourier’s law.
00 xb − xi
Tb = Ti + q
k
where Ti and Tb are interior and boundary temperatures, respectively.
ME 448/548: 2D Diffusion Equation
21
Boundary Conditions (part 2)
Boundary Boundary
type
Condition
3
Convection
Post-processing in fvpost
From specified h and T∞ , compute boundary temperature and heat flux
through the cell face on the boundary. Continuity of heat flux requires
T − Ti
−k b
= h(Tb − Tamb )
x b − xi
which can be solved for Tb to give
Tb =
hTamb + (k/δxe )Ti
h + (k/δxe )
where δxe = xb − xi
4
Symmetry
q 00 = 0. Set boundary Tb equal to adjacent interior Ti .
ME 448/548: 2D Diffusion Equation
22
Matlab codes for obtaining the numerical solution
A set of general purpose codes has been written to facilitate experimentation with the
CVFD method.
Algorithm Tasks
Core Routines
Define the mesh
fvUniformMesh or
fvUniBlockMesh
Define boundary conditions
Compute finite-volume coefficients for interior cells
fvcoef
Adjust coefficients for boundary conditions
fvbc
Solve system of equations
Assemble coefficient matrix
fvAmatrix
Solve
Compute boundary values and/or fluxes
fvpost
Plot results
ME 448/548: 2D Diffusion Equation
23
Model Problem 1
Choose a source term that may be physically unrealistic, but one that gives an exact
solution that is easy to evaluate
" 2 2#
2π
πx
2πy
π
+
sin
S=
sin
Lx
Ly
Lx
Ly
The exact solution is
φ = sin
πx
Lx
sin
2πy
Ly
Main code to solve this problem is in demoModel1.m
ME 448/548: 2D Diffusion Equation
24
Model Problem 1
The exact solution is
2
φ = sin
πx
Lx
sin
2πy
Ly
1.5
1
0.5
0
1
0.5
0
x
ME 448/548: 2D Diffusion Equation
1
0.8
0.6
0.4
0.2
0
y
25
Solutions to Model Problem 1 Show Correct Truncation Error
kek2
=
N
qP
N
e2i
√
e¯
N e¯2
∼
=√ .
N
N
−1
10
Measured
Theoretical error ~ (∆ x)3
−2
10
−3
10
Measured error
The local truncation error at each node is
ei ∼ O(∆x2).
Since the exact solution is known we can
compute
−4
10
−5
10
−6
10
−7
where N = nxny is the total number of
interior nodes in the domain, and e¯ is the
average truncation error per node.
10
−8
10
−3
10
−2
−1
10
10
0
10
∆x
Since ei ∼ O(∆x2), N ∼ n2x, and ∆x = Lx/(nx + 1), we can estimate
3
O L2x/(nx + 1)2
kek2
e¯
O(∆x2)
1
3
∼√ =
=
∼O
= O(∆x ).
N
nx
nx
nx
N
ME 448/548: 2D Diffusion Equation
26
Model Problem 2
Uniform source term: S = 1.
0.08
Analytical solution is an infinite
series
Code in demoModel2.m
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
1
0.5
0
x
ME 448/548: 2D Diffusion Equation
1
0.8
0.6
0.4
0.2
0
y
27
Model Problem 3
Heat conduction in a rectangle consisting of
two material regions.
• Inner rectangle with high conductivity
• Outer rectangle with low conductivity
• Inner region has uniform heat source
Ly
0.5Ly
Γ2 = αΓ1
Γ1 = 1
S2 = 1000
S1 = 0
0.25Ly
0.25Lx
0.5Lx
Lx
The discontinuity and difference in material properties can be used to stress the solution
algorithm.
Γ2
α=
Γ1
The analytical solution does not exist. Code in demoModel3.m
ME 448/548: 2D Diffusion Equation
28
Model Problem 3
Solution with α = 100
35
Ly
0.5Ly
Γ2 = αΓ1
Γ1 = 1
S2 = 1000
S1 = 0
30
25
0.25Ly
20
0.25Lx
0.5Lx
15
Lx
10
5
0
1
0.5
0
x
ME 448/548: 2D Diffusion Equation
1
0.8
0.6
0.4
0.2
0
y
29
Model Problem 4: Fully Developed Flow in a Rectangular Duct
For simple fully-developed flow the governing equation for the axial velocity w is
"
#
2
2
∂ w
∂ w
dp
µ
+
=0
−
∂x2
∂y 2
dz
This corresponds to the generic model equation with
φ = w,
ME 448/548: 2D Diffusion Equation
Γ = µ (= constant),
S=−
dp
.
dz
30
Model Problem 4: Fully Developed Flow in a Rectangular Duct
The symmetry in the problem allows alternative ways of defining the numerical model
Full Duct
Quarter Duct
Ly
Ly
y
x
y
Lx
x
Lx
For the full duct simulation depicted on the left hand side, the boundary conditions are no
slip conditions on all four walls.
w(x, 0) = w(x, Ly ) = w(0, y) = w(Lx, y) = 0.
ME 448/548: 2D Diffusion Equation
(full duct)
31
Model Problem 4: Fully Developed Flow in a Rectangular Duct
The symmetry in the problem allows alternative ways of defining the numerical model
Full Duct
Quarter Duct
Ly
Ly
y
x
y
Lx
x
Lx
For the quarter duct simulation depicted on the right hand side, the boundary conditions
are no slip conditions on the solid walls (x = Lx and y = Ly )
w(Ly , y) = w(x, Lx) = 0
(quarter duct)
and symmetry conditions on the other two planes
∂u ∂u =
= 0.
∂x x=0
∂y y=0
ME 448/548: 2D Diffusion Equation
32

Download Report