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© 2000 by Lizette R. Chevalier
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
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
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
 2 2
2
x
c t
y or t
set up a grid
estimate the dependent
variable at the center
x
y or t
set up a grid
estimate the dependent
variable at the center
or intersections
of the grid
x
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
Specific Study Objectives
• Know the difference between convergence
and stability of parabolic PDE
• Understand the difference between explicit
and implicit schemes for solving 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
Elliptic Equations
2
B - 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
First, recognize how the shape can be set
in an x-y coordinate system
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 by at its edges or
boundaries, where the temperature can be set.
y
Dy
Dx
x
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
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
x
Dx 2
 2u ui , j 1  2ui , j  ui , j 1

2
y
Dy 2
(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
If D x = D y
we can collect the terms
to get:
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
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
y
(0,100)
h
0
x
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
2,1
h i+1,j = h i-1,j
and
h
0
y
h i,j+1 = h i,j-1
h1,1 = (2h1,2 + 2 h2,1)/4
h
0
x
1,2
2,2
1,1
2,1
h
0
y
h1,2 = (h1,1 + h1,3+2h22)/4
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
i, j
u
 u
new
i, j
 1   u
old
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 pre-specified
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 D x = D y = 20 m
y
(0,100)
h
0
x
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.
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
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
Note:
q
 is in degrees
1  y 
  tan  
If qx=0,  is 90 or 270 
 qx 
depending on whether qy
for q x  0
is positive or negative,
q
1  y 
  tan    180 respectively
 qx 
for q x  0
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.
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
spatial
x
 2T T
k 2 
x
t
Ti l 1  2Ti l  Ti l 1 Ti l 1  Ti l
k

2
DT
Dx 
Centered finite divided difference
Forward finite divided
difference
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  2Ti1l  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
Implicit
T  2Ti  T
Dx 2
l
i 1
l
l
i 1
Ti l 11  2Ti l 1  Ti l 11
Dx 2
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
grid point involved with space 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
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
grid point involved with space difference
grid point involved with time difference
With the implicit method, we develop a set of simultaneous
equations at step in time
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 
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