Mechanical Engineering Mechanical Engineering

Mech3750 2D Elliptical Equations

Laplace's Equation and Poisson's Equation

Consider the 2D unsteady heat conduction equation (in non-dimensional form)1

q
t
 = 2 q
x2
 + 2 q
y2
(1)
which includes heat flows in the y-direction as well as the x-direction. For example, this equation would describe the cooling of a hot plate between two insulating layers where heat is lost from the plate only at the edges (in the x and y directions, not in the direction normal to the plate).

Depending on the boundary conditions, a steady state temperature distribution may be reached in the plate.


For example, each edge of the plate may be held at constant temperature. A steady state temperature distribution is shown in the figure: the North boundary is at a temperature 3, the south boundary at temperature 2, the west boundary at temperature 1 and the east boundary at temperature 0. The direction of heat flow is always normal to the lines of constant temperature (normal to the temperature contours) and from high to low temperature.

In the steady state, the derivative q/ t is zero throughout the domain and equation (1) becomes

2 q
x2
 + 2 q
y2
 = 0.
(2)
This is known as Laplace's equation which has many applications in physics. For instance, it can describe the distribution of electrical potential in a two-dimensional conducting plate (the temperature contours in the figure could just as easily be contours of electrical potential when fixed potentials (voltages) are applied to the four sides of the plate). It can be used to describe the flow of a perfect (inviscid) fluid - the lines normal to the contour lines can be interpreted as streamlines of a fluid flow through the region.

Finite difference form of Laplace's equation

Let the x-y plane be divided into a regular grid with node values given by
xj = jDx,  j = 1,2,3 ...    yk = k Dy,  k = 1,2,3 ...
We seek a numerical solution of (2) in the form of discrete node values qj,k.

Write the finite difference form of (2) at node (j,k). Using the central difference approximation

2 q
x2

j,k 
 = qj-1,k - 2qj,k + qj+1,k
Dx2
 + O(Dx2),
and an equivalent expression for 2 q/y2, we have

qj-1,k - 2qj,k + qj+1,k
Dx2
 + qj,k+1 - 2qj,k + qj,k+1
Dy2
 = 0.
(3)
If this equation is satisfied at every node by a set of values qj,k then that set of values is an exact solution of the finite difference form of Laplace's equation (3), and hence an approximate solution of the original Laplace's equation (2).

By multiplying through by Dx2 Dy2 and re-arranging we get

qj,k = Dy2 ( qj-1,k + qj+1,k ) + Dx2 ( qj,k-1 + qj,k+1 )
2 (Dx2 + Dy2)
.
(4)

Iterative solution

The equation (4) can form the basis of an iterative method of solution, provided the values of qj,k are known and constant on the boundary of the region. We can start with a set of interior values qj,k and derive a new set of values by applying the above equation at each internal node to derive a new estimate of the value qj,k. That is, the new value of qj,k is a function, given by (4), of its four neighbouring values qj-1,k, qj+1,k, qj,k-1 and qj,k+1. After all node values have been changed according to (4) the equation will, in general, no longer be valid at all the nodes. However, the procedure can be repeated, successively replacing each node value by the new value calculated by (4), until the change in all the nodes values is very small and equation(4) is satisfied at each node to sufficient accuracy.

For simple grid Dx = Dy

Consider the simple case where Dx = Dy. Then (4) reduces to
qj,k = 1
4
 ( qj-1,k + qj+1,k + qj,k-1 + qj,k+1 )
(5)
A simple code to find the solution for given fixed boundary conditions will look something like this: 
set boundary values of theta
initial guess theta = 0 may be as good as any
repeat 
   for each internal (non-boundary) node j,k
      theta(j,k) = 0.25*(theta(j-1,k) + theta(j+1,k) + ...
                           theta(j,k-1) + theta(j,k+1);
   end
until all changes in theta are small
Note that the above code is written so that the newest node values are always used in the updating formula. This can inhibit vectorisation (and computational efficiency) and sometimes the code is written using two levels or layers of node values. The next estimate of node values is made to depend on the previous level only. A possible coding is as follows: 
last = 1; next = 2;
set boundary values of theta(last,j,k)
set initial values of theta(last,j,k)
repeat 
   for each internal (non-boundary) node j,k
      theta(next,j,k) = 0.25*(theta(last,j-1,k) + theta(last,j+1,k) + ...
                           theta(last,j,k-1) + theta(last,j,k+1);
      temp = last; last = next; next = temp; % swap the 'next' and 'last' levels
   end
until all changes in theta(last,j,k) - theta(next,j,k) are small

Footnotes:

1 That is, the non-dimensional variables are

  1. t = a×time/ Lchar where Lchar is a characteristic length and a is the thermal diffusivity
  2. x and y = length/Lchar
  3. q = (T - Tchar)/ Tchar, where Tchar is a characteristic temperature.

File translated from TEX by TTH, version 2.20.
On 8 May 2000, 09:40.