Transcript Partial Differential Equations

PARTIAL DIFFERENTIAL EQUATIONS
ENGR 351
Numerical Methods for Engineers
Southern Illinois University Carbondale
College of Engineering
Dr. L.R. Chevalier
Dr. B.A. DeVantier
Copyright© 2003 by Lizette R. Chevalier and Bruce A. DeVantier
Permission is granted to students at Southern Illinois University at Carbondale
to make one copy of this material for use in the class ENGR 351, Numerical
Methods for Engineers. No other permission is granted.
All other rights are reserved. No part of this publication may be reproduced,
stored in a retrieval system, or transmitted, in any form or by any means,
electronic, mechanical, photocopying, recording, or otherwise, without
the prior written permission of the copyright owner.
Partial Differential Equations
 An equation involving partial derivatives of
an unknown function of two or more
independent variables
 The following are examples. Note: u depends
on both x and y
3
 u 
 u
 u
 3u
 2 xy 2  u  1  2   6
x
2
2
x
y
xy
 x 
 2u
 2u
 2u
u
 x 2  8u  5 y
 xu  x
2
xy
y
x
y
2
2
2
Partial Differential Equations
 Because of their widespread application in
engineering, our study of PDE will focus on
linear, second-order equations
 The following general form will be evaluated
for B2 - 4AC
 2u
 2u
 2u
A 2 B
C 2 D 0
x
xy
y
B2-4AC
<0
Category
Elliptic
Example
Laplace equation (steady state with
2 spatial dimensions
 2T  2T
 2 0
2
x
y
=0
Parabolic
Heat conduction equation (time variable
with one spatial dimension
 2 T T
k 2 
x
t
>0
Hyperbolic
Wave equation (time-variable with one
spatial dimension
2y 1 2y

x 2 c2 t 2
Scope of Lectures on PDE
 Finite Difference: Elliptic



The Laplace Equation
Finite difference solution
Boundary conditions
 Finite Difference: Parabolic



Heat conduction
Explicit method
Simple implicit method
Specific Study Objectives
 Recognize the difference between elliptic,
parabolic, and hyperbolic PDE
 Recognize that the Liebmann method is equivalent
to the Gauss-Seidel approach for solving
simultaneous linear algebraic equations
 Recognize the distinction between Dirichlet and
derivative boundary conditions
 Know the difference between convergence and
stability of parabolic PDE
Finite Difference: Elliptic Equations
B2- 4AC < 0
 Typically used to characterize steady-state
boundary value problems
 Before solving, the Laplace equation will be
solved from a physical problem
 2u
 2u
 2u
A 2 B
C 2 D0
x
xy
y
 2u  2u
 2 0
2
x
y
The Laplace Equation
 u  u
 2 0
2
x y
2
2
 Models a variety of problems involving the potential
of an unknown variable
 We will consider cases involving thermodynamics,
fluid flow, and flow through porous media
The Laplace Equation
 2T  2T
 2 0
2
 x y
 Let’s consider the case of a plate heated
from the boundaries
 How is this equation derived from basic
concepts of continuity?
 How does it relate to flow fields?
Consider the plate below, with thickness Dz.
The temperatures are known at the boundaries.
What is the temperature throughout the plate?
T = 400
T = 200
T= 200
T = 200
Divide into a grid, with increments by Dx and Dy
y
T = 400
T = 200
T= 200
T = 200
x
Divide into a grid, with increments by Dx and Dy
y
T = 400
T = 200
T= 200
T = 200
x
What is the temperature here, if
using a block centered scheme?
y
T = 400
T = 200
T= 200
T = 200
x
What is the temperature here, if
using a grid centered scheme?
y
T = 400
T = 200
T= 200
T = 200
x
Consider the element shown below on the face of
a plate D z in thickness.
The plate is illustrated everywhere but at its edges or
boundaries, where the temperature can be set.
y
Dy
x
Dx
q(y + D y)
q(x + D x)
q(x)
Consider the steady state heat flux q
in and out of the elemental volume.
q(y)
By continuity, the flow of heat in must equal the flow of heat
out.
q x DyDzDt  q y DxDzDt 
q x  Dx DyDzDt  q y  Dy DxDzDt
q x DyDzDt  q y DxDzDt 
q x  Dx DyDzDt  q y  Dy DxDzDt
Rearranging terms
[ q x  Dx DyDzDt  q x DyDzDt ]
 [ q y  Dy DxDzDt  q y DxDzDt ]  0
Dividing by Dx, Dy, Dz and D t :
q x  Dx   q x  q y  Dy   q y 

0
Dx
Dy
As Dx & Dy approach zero, the equation
reduces to:
q x q y

0
x y
Again, this is our
continuity equation
 q q

0
x y
Equation A
The link between flux and temperature is
provided by Fourier’s Law of heat conduction
T
qi   kC
i
Equation B
where qi is the heat flux in the direction i.
Substitute B into A to get the Laplace
equation
 q q

0
x y
qi   kC
T
i
Equation A
Equation B
 q q  
T   
T 


   kC
    kC
x y x 
x  y 
y 
 2T  2T

0
2
2
x
y
Consider Fluid Flow
In fluid flow, where the fluid is a liquid or a gas, the
continuity equation is:
Vx Vy

0
x
y
The link here can by either of the following sets of equations:
The potential function:
Stream function:

Vx 
x

Vy 
y

Vx 
y

Vy  
x
Vx Vy

0
x
y

x
Vy 

y

y
Vy  

x
Vx 
Vx 
The Laplace equation is then
 2  2
 2 0
2
x y
 2  2
 2 0
2
x
y
Flow in Porous Media
 q q

0
x y
H
qi   K
i
Darcy’s Law
The link between flux and the pressure head is provided by
Darcy’s Law
 2h  2h
 2 0
2
x y
Poisson Equation
For a case with sources and sinks within the 2-D
domain, as represented by f(x,y), we have the
Poisson equation.
 2u  2u
 2  f (x, y)
2
x
y
Now let’s consider solution techniques.
Evaluate these equations based on the grid and
central difference equations
 2u ui 1, j  2ui , j  ui 1, j

2
Dx 2
x
 2u ui , j 1  2ui , j  ui , j 1

2
Dy 2
y
(i,j+1)
(i,j)
(i+1,j)
(i-1,j)
(i,j-1)
ui 1, j  2ui , j  ui 1, j
Dx
2

ui , j 1  2ui , j  ui , j 1
Dy
2
0
If D x = D y
we can collect the terms
to get:
ui 1, j  ui 1, j  ui , j 1  ui , j 1  4ui , j  0
(i,j+1)
(i,j)
(i+1,j)
(i-1,j)
(i,j-1)
ui 1, j  ui 1, j  ui , j 1  ui , j 1  4ui , j  0
This equation is referred
to as the Laplacian difference
equation.
(i,j+1)
(i,j)
(i+1,j)
(i-1,j)
It can be applied to all
interior points.
(i,j-1)
We must now consider
what to do with the
boundary nodes.
Boundary Conditions
 Dirichlet boundary conditions: u is specified at the
boundary
 Temperature
 Head
 Neumann boundary condition: the derivative is
specified
 qi
h
T
or
xi
xi
 Combination of both u and its derivative (mixed
boundary condition)
The simplest case is where the boundaries are
specified as fixed values.
This case is known as the Dirichlet boundary
conditions.
u2
u1
u3
u4
Consider how we can deal with the lower node shown, u1,1
u2
u1
1,2
u3
1,1
Note:
This grid would result
in nine simultaneous
equations.
2,1
u4
-4u1,1 +u1,2+u2,1+u1 +u4 = 0
Let’s consider how to model the Neumann boundary condition
 u ui 1, j  ui 1, j

x
2 Dx
(0,100)
h
0
x
y
centered finite divided difference
approximation
h = 0.05x + 100
2h 2h

0
2
2
x
y
(200,100)
h
0
x
x
(0,0)
h
0
y
(200,0)
suppose we wanted to consider this end
grid point
h
0
x
The two boundaries are
consider to be
symmetry lines due to
the fact that the BC
translates in the finite
difference form to:
1,2
1,1
h
0
y
2,1
h
i+1,j
= h i-1,j
and
h
i,j+1
= h i,j-1
h1,1 = (2h1,2 + 2 h2,1)/4
h
0
x
h1,2 = (h1,1 + h1,3+2h22)/4
1,2
1,1
h
0
y
2,2
2,1
Example
The grid on the next slide is designed to solve the
LaPlace equation
 2  2
 2 0
2
x
y
Write the finite difference equations for the nodes
(1,1), (1,2), and (2,1). Note that the lower
boundary is a Dirichlet boundary condition, the left
boundary is a Neumann boundary condition, and Dx
= Dy.
The Liebmann Method
 Most numerical solutions of the
Laplace equation involves systems that
are much larger that the general
system we just evaluated
 Note that there are a maximum of five
unknown terms per line
 This results in a significant number of
terms with zero’s
The Liebmann Method
 In addition to the fact that they are prone to
round-off errors, using elimination methods
on such sparse system waste a great amount
of computer memory storing zeros
 Therefore, we commonly employ approaches
such as Gauss-Seidel, which when applied to
PDEs is also referred to as Liebmann’s
method.
The Liebmann Method
 In addition the equations will lead to a
matrix that is diagonally dominant.
 Therefore the procedure will converge to a
stable solution.
 Over relaxation is often employed to
accelerate the rate of convergence
ui 1, j  ui 1, j  ui , j 1  ui , j 1  4ui , j  0
new
old
uinew


u

1


u

 i, j
,j
i, j
ui 1, j  ui 1, j  ui , j 1  ui , j 1  4ui , j  0
new
old
uinew


u

1


u

 i, j
,j
i, j
As with the conventional Gauss Seidel method, the
iterations are repeated until each point falls below a prespecified tolerance:
s 
old
uinew

u
,j
i, j
new
i, j
u
 100
Groundwater Flow Example
Modeling 1/2 of the system shown, we can develop the
following schematic where Dx = Dy = 20 m
(0,100)
h
0
x
y
h = 0.05x + 100
2h 2h

0
2
2
x
y
(200,100)
h
0
x
x
(0,0)
h
0
y
(200,0)
The finite difference equations can be solved using a
a spreadsheet. This next example is part of the PDE
example you can download from my homepage.
CAN USE EXCEL DEMONSTRATION
100
=A1+0.05*20
=B1+0.05*20
=C1+0.05*20
=D1+0.05*20
=E1+0.05*20
=F1+0.05*20
=G1+0.05*20
=H1+0.05*20
=I1+0.05*20
=J1+0.05*20
=(A1+2*B2+A3)/4
=(B1+C2+B3+A2)/4
=(C1+D2+C3+B2)/4
=(D1+E2+D3+C2)/4
=(E1+F2+E3+D2)/4
=(F1+G2+F3+E2)/4
=(G1+H2+G3+F2)/4
=(H1+I2+H3+G2)/4
=(I1+J2+I3+H2)/4
=(J1+K2+J3+I2)/4
=(K1+K3+2*J2)/4
=(A2+2*B3+A4)/4
=(B2+C3+B4+A3)/4
=(C2+D3+C4+B3)/4
=(D2+E3+D4+C3)/4
=(E2+F3+E4+D3)/4
=(F2+G3+F4+E3)/4
=(G2+H3+G4+F3)/4
=(H2+I3+H4+G3)/4
=(I2+J3+I4+H3)/4
=(J2+K3+J4+I3)/4
=(K2+K4+2*J3)/4
=(A3+2*B4+A5)/4
=(B3+C4+B5+A4)/4
=(C3+D4+C5+B4)/4
=(D3+E4+D5+C4)/4
=(E3+F4+E5+D4)/4
=(F3+G4+F5+E4)/4
=(G3+H4+G5+F4)/4
=(H3+I4+H5+G4)/4
=(I3+J4+I5+H4)/4
=(J3+K4+J5+I4)/4
=(K3+K5+2*J4)/4
=(A4+2*B5+A6)/4
=(B4+C5+B6+A5)/4
=(C4+D5+C6+B5)/4
=(D4+E5+D6+C5)/4
=(E4+F5+E6+D5)/4
=(F4+G5+F6+E5)/4
=(G4+H5+G6+F5)/4
=(H4+I5+H6+G5)/4
=(I4+J5+I6+H5)/4
=(J4+K5+J6+I5)/4
=(K4+K6+2*J5)/4
=(2*A5+2*B6)/4
=(2*B5+C6+A6)/4
=(2*C5+D6+B6)/4
=(2*D5+E6+C6)/4
=(2*E5+F6+D6)/4
=(2*F5+G6+E6)/4
=(2*G5+H6+F6)/4
=(2*H5+I6+G6)/4
=(2*I5+J6+H6)/4
=(2*J5+K6+I6)/4
=(2*K5+2*J6)/4
100
=A1+0.05*20
=B1+0.05*20
=C1+0.05*20
=D1+0.05*20
=E1+0.05*20
=F1+0.05*20
=G1+0.05*20
=H1+0.05*20
=I1+0.05*20
=J1+0.05*20
=(A1+2*B2+A3)/4
=(B1+C2+B3+A2)/4
=(C1+D2+C3+B2)/4
=(D1+E2+D3+C2)/4
=(E1+F2+E3+D2)/4
=(F1+G2+F3+E2)/4
=(G1+H2+G3+F2)/4
=(H1+I2+H3+G2)/4
=(I1+J2+I3+H2)/4
=(J1+K2+J3+I2)/4
=(K1+K3+2*J2)/4
=(A2+2*B3+A4)/4
=(B2+C3+B4+A3)/4
=(C2+D3+C4+B3)/4
=(D2+E3+D4+C3)/4
=(E2+F3+E4+D3)/4
=(F2+G3+F4+E3)/4
=(G2+H3+G4+F3)/4
=(H2+I3+H4+G3)/4
=(I2+J3+I4+H3)/4
=(J2+K3+J4+I3)/4
=(K2+K4+2*J3)/4
=(A3+2*B4+A5)/4
=(B3+C4+B5+A4)/4
=(C3+D4+C5+B4)/4
=(D3+E4+D5+C4)/4
=(E3+F4+E5+D4)/4
=(F3+G4+F5+E4)/4
=(G3+H4+G5+F4)/4
=(H3+I4+H5+G4)/4
=(I3+J4+I5+H4)/4
=(J3+K4+J5+I4)/4
=(K3+K5+2*J4)/4
=(A4+2*B5+A6)/4
=(B4+C5+B6+A5)/4
=(C4+D5+C6+B5)/4
=(D4+E5+D6+C5)/4
=(E4+F5+E6+D5)/4
=(F4+G5+F6+E5)/4
=(G4+H5+G6+F5)/4
=(H4+I5+H6+G5)/4
=(I4+J5+I6+H5)/4
=(J4+K5+J6+I5)/4
=(K4+K6+2*J5)/4
=(2*A5+2*B6)/4
=(2*B5+C6+A6)/4
=(2*C5+D6+B6)/4
=(2*D5+E6+C6)/4
=(2*E5+F6+D6)/4
=(2*F5+G6+E6)/4
=(2*G5+H6+F6)/4
=(2*H5+I6+G6)/4
=(2*I5+J6+H6)/4
=(2*J5+K6+I6)/4
=(2*K5+2*J6)/4
100
=(A1+2*B2+A3)/4
=(A2+2*B3+A4)/4
=(A3+2*B4+A5)/4
=(A4+2*B5+A6)/4
=(2*A5+2*B6)/4
You will get an
error message in
Excel that state that
it will not resolve
a circular reference.
100
=A1+0.05*20
=B1+0.05*20
=C1+0.05*20
=D1+0.05*20
=E1+0.05*20
=F1+0.05*20
=G1+0.05*20
=H1+0.05*20
=I1+0.05*20
=J1+0.05*20
=(A1+2*B2+A3)/4
=(B1+C2+B3+A2)/4
=(C1+D2+C3+B2)/4
=(D1+E2+D3+C2)/4
=(E1+F2+E3+D2)/4
=(F1+G2+F3+E2)/4
=(G1+H2+G3+F2)/4
=(H1+I2+H3+G2)/4
=(I1+J2+I3+H2)/4
=(J1+K2+J3+I2)/4
=(K1+K3+2*J2)/4
=(A2+2*B3+A4)/4
=(B2+C3+B4+A3)/4
=(C2+D3+C4+B3)/4
=(D2+E3+D4+C3)/4
=(E2+F3+E4+D3)/4
=(F2+G3+F4+E3)/4
=(G2+H3+G4+F3)/4
=(H2+I3+H4+G3)/4
=(I2+J3+I4+H3)/4
=(J2+K3+J4+I3)/4
=(K2+K4+2*J3)/4
=(A3+2*B4+A5)/4
=(B3+C4+B5+A4)/4
=(C3+D4+C5+B4)/4
=(D3+E4+D5+C4)/4
=(E3+F4+E5+D4)/4
=(F3+G4+F5+E4)/4
=(G3+H4+G5+F4)/4
=(H3+I4+H5+G4)/4
=(I3+J4+I5+H4)/4
=(J3+K4+J5+I4)/4
=(K3+K5+2*J4)/4
=(A4+2*B5+A6)/4
=(B4+C5+B6+A5)/4
=(C4+D5+C6+B5)/4
=(D4+E5+D6+C5)/4
=(E4+F5+E6+D5)/4
=(F4+G5+F6+E5)/4
=(G4+H5+G6+F5)/4
=(H4+I5+H6+G5)/4
=(I4+J5+I6+H5)/4
=(J4+K5+J6+I5)/4
=(K4+K6+2*J5)/4
=(2*A5+2*B6)/4
=(2*B5+C6+A6)/4
=(2*C5+D6+B6)/4
=(2*D5+E6+C6)/4
=(2*E5+F6+D6)/4
=(2*F5+G6+E6)/4
=(2*G5+H6+F6)/4
=(2*H5+I6+G6)/4
=(2*I5+J6+H6)/4
=(2*J5+K6+I6)/4
=(2*K5+2*J6)/4
=B1+0.05*20
=(C1+D2+C3+B2)/4
=(C2+D3+C4+B3)/4
=(C3+D4+C5+B4)/4
=(C4+D5+C6+B5)/4
=(2*C5+D6+B6)/4
After selecting the
appropriate command,
EXCEL with perform
the Liebmann method
for you.
In fact, you will be able
to watch the iterations.
Table 2: Results of finite difference model.
A
1
2
3
4
5
6
100
101.6
102.5
103
103.3
103.3
B
101
102
102.7
103.1
103.3
103.4
C
D
E
102
103
104
102.6 103.4 104.2
103.1 103.7 104.3
103.4 103.9 104.4
103.6
104 104.5
103.7
104 104.5
F
G
105
105
105
105
105
105
H
I
J
K
106
107
108
109
110
105.8 106.6 107.4
108 108.4
105.7 106.3 106.9 107.3 107.5
105.6 106.1 106.6 106.9
107
105.5
106 106.4 106.7 106.7
105.5
106 106.3 106.6 106.7
Table 2: Results of finite difference model.
A
1
2
3
4
5
6
100
101.6
102.5
103
103.3
103.3
B
101
102
102.7
103.1
103.3
103.4
C
D
E
102 103 104
102.6 103.4 104.2
103.1 103.7 104.3
103.4 103.9 104.4
103.6 104 104.5
103.7 104 104.5
F
G
105
105
105
105
105
105
H
I
J
K
106 107 108 109 110
105.8 106.6 107.4 108 108.4
105.7 106.3 106.9 107.3 107.5
105.6 106.1 106.6 106.9 107
105.5 106 106.4 106.7 106.7
105.5 106 106.3 106.6 106.7
• Assuming K = 5 m/day, we can calculate the q’s
• Use centered differences, for example @ cell B2
• qx = -K h/x = -(5m/d) (102.6m-101.6m) /(2*20m) = -0.125 m/d
• qy = -K h/y = -(5m/d) (101 m-102.7 m) /(2*20m) = 0.2125 m/d
...end of problem.
Secondary Variables
 Because its distribution is described by the
Laplace equation, temperature is considered
to be the primary variable in the heated
plate problem
 A secondary variable may also be of interest
 In this case, the second variable is the rate
of heat flux across the place surface
T
qi   kC
i
Secondary Variables
Secondary Variables
T
qi   kC
i
q x  k '
q y  k '
Ti 1. j  Ti 1, j
2Dx
Ti. j 1  Ti , j 1
2Dy
FINITE DIFFERENCE
APPROXIMATION
BASED ON RESULTS
OF TEMPERATURE
DISTRIBUTION
The Resulting Flux Is A Vector With
Magnitude And Direction
q n  q x2  q y2
 qy
  tan 
 qx
for q x  0



 qy
  tan 
 qx
for q x  0

  180

1
1
Note:
 is in degrees
If qx=0,  is 90 or 270 
depending on whether
qy is positive or
negative, respectively
Finite Difference: Parabolic
Equations B2- 4AC = 0
 2u
 2u
 2u
A 2 B
C 2 D0
x
xy
y
These equations are used to characterize
transient problems.
We will first study this in one spatial direction
then we will discuss the results in 2-D.
Finite Difference: Parabolic
Equations B2- 4AC = 0
Consider the heat-conduction equation
 2 T T
k 2 
x
t
As with the elliptic PDEs, parabolic equations
can be solved by substituting finite difference
equations for the partial derivatives.
However we must now consider changes in
time as well as space.
y
t
x
temporal
{
l
i
u
{
x
spatial
 T T
k 2 
x
t
2
T  2Ti  T
k
2
Dx 
l
i 1
l
l
i 1

Ti
l 1
 Ti
Dt
l
We can further reduce the equation:
Ti l 1  2Ti l  Ti l 1 Ti l 1  Ti l
k

2
Dt
Dx 
Ti l 1  Ti l   Ti l 1  2Ti l  Ti l 1 
where
  kDt
Dx 
2
NOTE:
Now the temperature
at a node is estimated
as a function of the
temperature at the
node, and surrounding
nodes, but at a previous
time
Example
Consider a thin insulated rod 10 cm long with
k = 0.835 cm2/s
Let D x = 2 cm and D t = 0.1 sec.
cold
hot
At t=0 the temperature of the rod is zero.
Convergence and Stability
 Convergence means that as D x and D t
approach zero, the results of the numerical
technique approach the true solution
 Stability means that the errors at any stage
of the computation are attenuated, not
amplified, as the computation progresses
 The explicit method is stable and convergent
if
  12
TL
To
Derivative Boundary Conditions
In our previous example To and TL were constant values.
However, we may also have derivative boundary conditions
T0i  1  T0i    T1i  2T0i  Ti1
Thus we introduce an imaginary point at i = -1
This point provides the vehicle for providing the derivative BC
For the case of qo = 0, so T-1 = T1 .
In this case the balance at node 0 is:
T0l 1  T0l  2T1l  2T0l 
TL
q0 = 0
Derivative Boundary Conditions
TL
q0 = 10
Derivative Boundary Conditions
For the case of qo = 10, we need to know k’ [= k/(C)].
Assuming k’ =1, then 10 = - (1) dT/dx,
or dT/dx = -10
Derivative Boundary Conditions
T1l  Tl1
10   k
2 Dx
Dx
l
l
T1  T1  20
k
 l
Dx
l 1
l
l
T0  T0   2T1  20
 2T0 
1


Implicit Method
 Explicit methods have problems
relating to stability
 Implicit methods overcome this but at
the expense of introducing a more
complicated algorithm
 In this algorithm, we develop
simultaneous equations
Explicit
Ti l 1  2Ti l  Ti l 1
Dx 2
Explicit
Ti l 1  2Ti l  Ti l 1
Dx 2
Implicit
Ti l 11  2Ti l 1  Ti l 11
Dx 2
Explicit
Ti l 1  2Ti l  Ti l 1
Dx 2
Implicit
Ti l 11  2Ti l 1  Ti l 11
2
Dx 
grid point involved with space difference
Explicit
Ti l 1  2Ti l  Ti l 1
Dx 2
Implicit
Ti l 11  2Ti l 1  Ti l 11
Dx 2
grid point involved with space difference
grid point involved with time difference
Explicit
Implicit
T  2Ti  T
Dx 2
l
i 1
l
l
i 1
Ti l 11  2Ti l 1  Ti l 11
Dx 2
With the implicit method, we develop a set of simultaneous
equations at step in time
Implicit Method
Ti l 11  2Ti l 1  Ti l 11 Ti l 1  Ti l
k

2
Dt
 Dx 
which can be expressed as:
 Ti l 11  1  2 Ti l 1  Ti l 11  Ti l
For the case where the temperature level is given at
the end by a function f0 i.e. x = 0
T0l 1  f 0 t l 1 
Implicit Method
Substituting
 Ti l 11  1  2 Ti l 1  Ti l 11  Ti l
T0l 1  f 0 t l 1 
1  2 Ti l 1  Ti l 11  Ti l  f 0 t l 1 
In the previous example problem, we get a 4 x 4 matrix
to solve for the four interior nodes for each time step