Просмотреть файл в отдельном окне: 0adb15a74268b47266628354ff227630.pdf

8

Partial Differential

Equations

8.1 INTRODUCTION

Different engineering disciplines solve different types of problems in their respective

fields. For mechanical engineers, they may need to solve the temperature change

within a solid when it is heated by the interior heat sources or due to a rise or

decrease of its boundary temperatures. For electrical engineers, they may need to

find the voltages at all circuit joints of a computer chip board. Temperature and

voltage are the variables in their respective fields. Hence, they are called field

variables. It is easy to understand that the value of the field variable is spacedependent and time-dependent. That is to say, that we are interested to know the

spatial and temporal changes of the field variable. Let us denote the field variable

as . and let the independent variables be xi which could be the time t, or, the space

coordinates as x, y, and z. In order not to overly complicate the discussion, we

introduce the general two-dimensional partial differential equation which governs

the field variable in the form of:

A(x1, x 2 )

?2 ?

?2 ?

?2 ?

+ B(x1, x 2 )

+ C(x1, x 2 ) 2

2

?x1

?x1?x 2

?x 2

?

?? ?? ?

= F? x1, x 2 , ?,

,

?x1 ?x 2 ??

?

(1)

where the coefficient functions A, B, and C in the general cases are dependent on

the variables x1 and x2, and the right-hand-side function F, called forcing function

may depend not only on the independent variables x1 and x2 but may also depend

on the first derivatives of . There are innumerable of feasible solutions for Equation

1. However, when the initial and/or boundary conditions are specified, only particular

solution(s) would then be found appropriate.

In this chapter, we will discuss three simple cases when A, B, and C are all

constants. The first case is a two-dimensional, steady-state heat conduction problem

involving temperature as the field variable and only the spatial distribution of the

temperature needs to be determined, Equation 1 is reduced to a parabolic partial

differential equation named after Poisson and Laplace when the forcing function F

is not equal to, or, equal to zero, respectively. This is a case when the coefficient

functions in Equation 1 are related by the condition B2–4AC0.

The reason that these problems are called parabolic, elliptical, and hyperbolic

is because their characteristic curves have such geometric features. Readers interested in exploring these features should refer to a textbook on partial differential

equations.

Details will be presented regarding how the forward, backward, and central

differences discussed in Chapter 4 are to be applied for approximating the first and

second derivative terms appearing in Equation 1. Repetitive algorithms can be

devised to facilitate programming for straight-forward computation of the spatial

and temporal changes of the field variable. Numerical examples are provided to

illustrate how these changes can be determined by use of either QuickBASIC,

FORTRAN, MATLAB, or, Mathematica programs.

Although explanation of the procedure for numerical solution of these three

types of problems is given only for the simple one- and two-dimensional cases, but

its extension to the higher dimension case is straight forward. For example, one may

attempt to solve the transient heat conduction problem of a thin plate by having two

space variables instead of one space variable for a long rod. The steady-state heat

conduction problem of a thin plate can be extended for the case of a three-dimensional solid, and the string vibration problem can be extended to a two-dimensional

membrane problem.

8.2 PROGRAM PARABPDE — NUMERICAL SOLUTION

OF PARABOLIC PARTIAL DIFFERENTIAL EQUATIONS

The program ParabPDE is designed for numerically solving engineering problems

governed by parabolic partial differential equation in the form of:

??

?2?

=a 2

?t

?x

(1)

and is a function of t and x and satisfies a certain set of supplementary conditions.

Equation 1 is called a parabolic partial differential equation. For example, could

be the temperature, T, of a longitudinal rod shown in Figure 1 and the parameter a

in Equation 1 could be equal to k/c where k, c, and are the thermal conductivity,

specific heat, and specific weight of the rod, respectively. To make the problem more

specific, the rod may have an initial temperature of 0°F throughout and it is completely insulated around its lateral surface and also at its right end. If its left end is

to be maintained at 100°F beginning at the time t = 0, then it is of interest to know

© 2001 by CRC Press LLC

FIGURE 1. could be the temperature, T, of a longitudinal rod.

how the temperatures along the entire length of the rod will be changing as the time

progresses. This is therefore a transient heat conduction problem. One would like

to know how long would it take to have the entire rod reaching a uniform temperature

of 100°F.

If the rod is made of a single material, k/c would then be equal to a constant.

Analytical solution can be found for this simple case.1 For the general case that the

rod may be composed of a number of different materials and the physical properties

k, c, and would not only depend on the spatial variable x but may also depend on

the temporal variable t. The more complicated the variation of these properties in x

and t, the more likely no analytical solution is possible and the problem can only

be solved numerically. The finite-difference approximation of Equation 1 can be

achieved by applying the forward difference for the first derivative with respective

to t and central difference for the second derivative with respect to x as follows (for

t at ti and x at xj):

?T Ti +1, j ? Ti, j

=?

?t

?t

and

?2 T Ti, j?1 ? 2 Ti, j + Ti, j+1

=?

?x 2

( ?x)2

If k/c is changing in time and also changing from one location to another, we

could designate it as ai,j. As a consequence, Equation 1 can then be written as:

Ti +1, j ? Ti, j

?t

= a i, j

Ti, j?1 ? 2 Ti, j + Ti, j+1

( ?x)2

(2)

Since the initial temperature distribution T is known, the above expression

suggests that for a numerical solution we may select an appropriate increment in t,

t, and the temperature be determined at a finite number of stations, N. It is advisable

to have these stations be equally spaced so that the increment x is equal to L/(N–1)

where L is the length of the rod, and the instants are to be designated as t1 = 0, t2 =

t,…, ti = (i–1)t, and the stations as x1 = 0, x2 = x,…, xj = (j–1)?x,…, and xN =

(N–1)x = L. The task at hand is then to find T(ti,xj) for i = 1,2,… and j = 1,2,…,N.

It can be noticed from Equation 2 that the there is only one temperature at ti + 1 and

can be expressed in terms of those at the preceding instant ti as:

© 2001 by CRC Press LLC

Ti +1, j = Ti, j +

a i, j?t

( ?x)2

(T

i , j?1

? 2 Ti, j + Ti, j+1

)

(3)

Equation 3 is to be used for j = 2 through j = N–1. For the last station, j = N,

which is insulated,the temperatures on both side of this station can be assumed to

be equal (the station N + 1 is a fictitious one!). The modified equation for this

particular station is:

Ti +1,N = Ti,N +

2a i,N ?t

( ?x)2

(T

i ,N ?1

? Ti,N

)

(4)

For generating the temperatures of the rod at N stations for any specified time

increment t until the temperatures are almost all equal to 100°F throughout, the

program ParabPDE has been applied. It is listed below along with a typical printout

of the results.

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Output

© 2001 by CRC Press LLC

QUICKBASIC VERSION

MATLAB APPLICATIONS

A MATLAB version of ParabPDE can be created easily by converting the

QuickBASIC program. The m file may be arranged as follows:

© 2001 by CRC Press LLC

For solving the sample transient temperature problem, this m file can be called

and interactive MATLAB instructions can be entered through keyboard to obtain

the temperature distribution of the rod at various times:

© 2001 by CRC Press LLC

FIGURE 2. A composite graph with axes labels and markings of the curves by making use

of the MATLAB commands xlabel, ylabel, and text.

Notice that results of temperature distributions which have been terminated using

five differentials of 1, 0.1, 0.01, 0.001, and 0.0001 Fahrenheit are saved in Tsave

and then later plotted. Since T is a row matrix, the transpose of T, T, is stored in

an appropriate column of Tsave. If the temperatures were rounded, the rod reaches

a uniform temperature distribution of 100°F when the required temperature differential is selected to be 0.001°F. If the fifth curve for the temperature differential

equal to 0.0001 is plotted, it will be too close to the forth curve and also the plot

function provides only four line types (solid, broken, dot, and center lines), the fifth

set of results is therefore not saved in Tsave. Figure 2 shows a composite graph with

axes labels and markings of the curves by making use of the MATLAB commands

xlabel, ylabel, and text.

MATHEMATICA APPLICATIONS

The heat conduction problem of an insulated rod previously discussed in the

versions for FORTRAN, QuickBASIC, and MATLAB can be solved by application

of Mathematica as follows:

© 2001 by CRC Press LLC

In[1]: = (n = 11; rk = 0.037; c = 0.212; rho = 168; dt = 1; dx = 0.1;

c1 = rk/c/rho*dt/dx^2; nm1 = n–1;)

In[2]: = (t = Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100; tdf = 1;

tm = 0; flag1 = 0;

In[3]: = (While[flag1 = = 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),

{j,2,nm1}];

tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);tm = tm + dt;flag1 = 1;

Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,

Continue],{i,n}];

Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm]; Print[N[tn,3]])

Out[3] = t = 24

{100., 65.7, 37.5, 18.5, 7.83, 2.85, 0.892, 0.241, 0.0561, 0.0116, 0.00394}

The temperature distribution of the heated rod after 24 seconds is same as

obtained by the FORTRAN, QuickBASIC, and MATLAB versions when every

component of two consecutive temperature distribution, kept as {t} and {tn}, differ

no more than the allowed set value of tdf = 1 degree in In[2]. Any component has

a difference exceeding the tdf value will cause flag1 to change from a value of 1 to

0 and the Break command in the second Do loop in In[2] to exit and to continue

the iteration. flag1 is created to control the While command which determines when

the iteration should be terminated. The N[tn,3] instructs the components of {tn} be

printed with 3 significant figures.

When tdf is changed to a value of 0.5, Mathematica can again be applied to yield

In[4]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;

tdf = 0.5; tm = 0; flag1 = 0;)

In[5]: = (While[flag1 = = 0,Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),

{j,2,nm1}];

tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);tm = tm + dt;flag1 = 1;

Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,

Continue],{i,n}];

Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm]; Print[N[tn,3]])

Out[5]: = t = 49

{100., 75.5, 53.2, 34.9, 21.3, 12., 6.24, 3., 1.35, 0.617, 0.413}

Notice that Mathematica takes only 49 seconds, one second less than that

required for the FORTRAN, QuickBASIC, and MATLAB versions. The reason is

that Mathematica keeps more significant digits in carrying out all computations.

To show more on the effect of changing the tdf value, the following Mathematica

runs are provided:

© 2001 by CRC Press LLC

In[6]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;

tdf = 0.1; tm = 0; flag1 = 0;)

In[7]: = (While[flag1 == 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),

{j,2,nm1}];

tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);

tm = tm + dt; flag1 = 1;

Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,

Continue],{i,n}];

Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm];

Print[N[tn,3]])

Out[7]: = t = 462

{100., 93.9, 88., 82.3, 77.1, 72.5, 68.5, 65.3, 63., 61.6, 61.1}

In[8]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;

tdf = 0.05; tm = 0; flag1 = 0;)

In[9]: = (While[flag1 == 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),

{j,2,nm1}];

tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);

tm = tm + dt; flag1 = 1;

Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,

Continue],{i,n}];

Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm];

Print[N[tn,3]])

Out[9]: = t = 732

{100., 97., 94., 91.2, 88.5, 86.2, 84.2, 82.6, 81.5, 80.8, 80.5}

In[10]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;

tdf = 0.0001; tm = 0; flag1 = 0;

In[11]: = (While[flag1 == 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),

{j,2,nm1}];

tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);

tm = tm + dt; flag1 = 1;

Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,

Continue],{i,n}];

Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm]; Print[N[tn,3]])

Out[11]: = t = 3159

{100., 100., 100., 100., 100., 100., 100., 100., 100., 100., 100.}

For tdf = 0.1 and tdf = 0.05, Mathematica continues to take lesser time than the

other version; but when tdf = 0.0001, Mathematica needs two additional seconds. The

© 2001 by CRC Press LLC

reason is that seven significant figures are required in the last case, rounding may

have resulted in earlier termination of the iteration when the FORTRAN, QuickBASIC, and MATLAB versions are employed.

8.3 PROGRAM RELAXATN — SOLVING ELLIPTICAL PARTIAL

DIFFERENTIAL EQUATIONS BY RELAXATION METHOD

The program Relaxatn is designed for solving engineering problems which are

governed by elliptical partial differential equation of the form:

?2 ? ?2 ?

+

= F(x, y)

?x 2 ?y 2

(1)

where is called the field function and F(x,y) is called forcing function. When the

steady-state heat conduction of a two-dimensional domain is considered,2 then the

field function becomes the temperature distribution, T(x,y), and the forcing function

becomes the heat-source function, Q(x,y). If the distribution of T is influenced only

by the temperatures at the boundary of the domain, then Q(x,y) = 0 and Equation

1 which is often called a Poisson equation is reduced to a Laplace equation:

?2 T ?2 T

+

=0

?x 2 ?y 2

(2)

The second-order, central-difference formulas (read the program DiffTabl) can

be applied to approximate the second derivatives in the above equation at an arbitrary

point x = xi and y = yj in the domain as:

T ? 2 Ti, j + Ti +1, j

?2 T

at x i , y j =? i ?1, j

2

?x

( ?x)2

and

T ? 2 Ti, j + Ti, j+1

?2 T

at x i , y j =? i, j?1

2

?y

(?y)2

By substituting both of the above equations into Equation 2 and taking equal

increments in both x- and y-directions, the reduced equation is, for x = y,

Ti, j =

(

1

T + Ti +1, j + Ti, j?1 + Ti, j+1

4 i ?1, j

)

(3)

The result is expected when the temperature distribution reaches a steady state

because it states that the temperature at any point should be equal to the average of

its surrounding temperatures.

© 2001 by CRC Press LLC

FIGURE 3. A plate, which initially (at time t = 0) had a temperature equal to 0°F throughout,

insulated along a portion of its boundary and suddenly heated at its upper left boundary to

maintain a linearly varying temperature Tb.

Before we proceed further, it is appropriate at this time to introduce a numerical

case. Suppose that a plate which initially (at time t = 0) has a temperature equal to

0°F throughout and is insulated along a portion of its boundary, is suddenly heated

at its upper left boundary to maintain a linearly varying temperature Tb as shown in

Figure 3. If this heating process is to be maintained, we are then interested in

knowing the temperature distribution, changed from its initial state of uniformly

equal to 0°F if given sufficient time to allow it to reach an equilibrium (steady) state.

Numerically, we intend to calculate the temperatures, denoted as a matrix [T], at a

selected number of locations. Therefore, the plate is first divided into a gridwork of

M rows and N columns along the x- and y-directions, respectively, as indicated in

Figure 1. The directions of x- and y-axes are so selected for the convenience of

associating them with the row and column of the temperature matrix [T] which is

of order M by N. The values of M and N should be so decided such that the

increments x and y are equal in order to apply Equation 3. To be more specific,

let M = 10 and N = 20 and the linear temperature variation along the upper left

boundary Tb be:

Ti,1 = 10(i ? 1)

for

i = 1, 2,…, 6

(4)

Here, Ti,j is to be understood as the temperature at the location (xi,yj). Equation

4 describes the temperature along the left boundary y = y1 but only for xi in the

range of i = 1 to i = 6.

Since the temperatures at the stations which are on the insulated boundaries of

the plate are also involved, these unknown temperatures need to be treated differently.

By an insulated boundary, it means that there is no heat transfer normal to the

boundary. Since the heat flow is depended on the temperature difference across that

© 2001 by CRC Press LLC

insulated boundary, mathematically it requires that T/n = 0 there, n being the

normal direction. At the vertical boundaries x = xM, we have T/n = T/x = 0 since

x is the normal direction. Based on the central difference and considering two

increments in the x direction, at a yjth station we can have:

T

? TM ?1, j

?T

at x M , y j =? M +1, j

=0

?x

2( ?x)

(

)

or

TM +1, j = TM ?1, j

(5)

Since xM + 1 is below the bottom boundary of the plate shown in Figure 3, there

is no need to calculate the temperatures there, however Equation 5 enables the

temperatures at the stations along the bottom boundary of the plate TM,j for j = 1 to

N to be averaged. Returning to Equation 4, we notice that if Equation 5 is substituted

into it, the resulting equation which relates only to three neighboring temperatures is:

TM, j =

(

1

2 TM ?1, j + TM, j?1 + TM, j+1

4

)

for j = 2, 3,…, N ? 1

(6)

Notice that j = 1 and j = N are not covered in Equation 6. These two cases

concerning the insulated stations at the left and right bottom corners of the plate

will be discussed after we address the two vertical, insulated boundaries y = y1 and

y = yN.

For the boundaries y = y1, ?T/?n becomes ?T/?y. Again, we can apply the central

difference for double y increments to obtain:

T ? Ti,2

?T

at (x1, y1 ) =? i,0

=0

?y

2( ?y)

or

Ti,0 = Ti,2

(7)

Thus, the modified Equation 3 for the left insulated boundary is:

Ti,1 =

(

1

T + Ti +1,1 + 2 Ti,2

4 i ?1,1

)

for

i = 7, 8, 9

In a similar manner, we can derive for the right insulated boundary y = yN

T

? Ti,N ?1

?T

at (x i , y N ) =? i,N +1

=0

?y

2( ?y)

or

© 2001 by CRC Press LLC

(8)

Ti,N +1 = Ti,N ?1

(9)

and

Ti,N =

(

1

T

+ Ti +1,N + 2 Ti,N ?1

4 i ?1,N

)

for i = 8, 9

(10)

Having derived Equations 6, 8, and 10, it is easy to deduce the two special

equation for the corner insulated stations to be:

TM,1 =

(

1

T + TM ?1,1

2 M ,2

)

(11)

and

TM,N =

(

1

T

+ TM ?1,N

2 M,N ?1

)

(12)

We have derived all equations needed for averaging the temperature at any station

of interest including those at the insulated boundaries by utilizing those at its

neighboring stations. It suggests that a continuous upgrading process can be developed which assumes that the neighboring temperatures are known. This so-called

relaxation method starts with an initial assumed distribution of temperature [T(0)]

and continues to use Equations 3 and 6 to 12 until the differences at all locations

are small enough. Mathematically, the process terminates when:

M

N

i =1

j=1

?? T

( k +1)

i, j

? Ti(,kj ) < ?

(13)

where is a prescribed tolerance of accuracy and k = 0,1,2,… is the number of

sweeps in upgrading the temperature distribution. Superscripts (k + 1) and (k) refer

to the improved and previous distributions, respectively. The order of sweep will

affect how the temperatures should be upgraded. For example, if the temperatures

are to be re-averaged from top to bottom and left to right, referring to Figure 1, then

Equation 3 is to be modified as:

Ti(,kj +1) =

(

1 ( k +1)

T

+ Ti(+k1), j + Ti(,kj?+11) + Ti(,kj+)1

4 i ?1, j

)

(14)

Notice that the neighboring temperatures in the row above, i–1, and in the column

to the left, j–1, have already been upgraded while those in the row below, i + 1, and

in the column to the right, j + 1, are yet to be upgraded. Similar modifications are

to be made to Equations 6 to 12 during relaxation.

© 2001 by CRC Press LLC

The program Relaxatn is developed according to the relaxation method

described above. For solving the problem shown in Figure 3, both FORTRAN and

QuickBASIC versions are made available for interactively specifying the tolerance.

Sample results are presented below.

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Results

The program Relaxatn is first applied for an interactively entered value of equal

to 100. Only one relaxation needs to be implemented as shown below. The temperature distribution for Sweep #1 is actually the initial assumed distribution. One

cannot assess how accurate this distribution is. The second run specifies that be

© 2001 by CRC Press LLC

equal to 1. The results show that 136 relaxation steps are required. For giving more

insight on how the relaxation has proceeded, Sweeps #10, #30, #50, #100, and #137

are presented for interested readers. It clearly indicates that a tolerance of equal to

100 is definitely inadequate.

© 2001 by CRC Press LLC

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Results

Irregular Boundaries

Practically, there are cases where the domain of heat conduction have boundaries

which are quite irregular geometrically as illustrated in Figure 4. For such cases, the

equation derived based on the relaxation method, Equation 3, which states that the

temperature at any point has the average value of those at its four neighboring points

if they are equally apart, has to be modified. The modified equation can be derived

using a simple argument applied in both x and y directions. For example, consider

the temperature at the point G, TG, in Figure 4. First, let us investigate the horizontal,

y direction (for convenience of associating x and y with the row and column indices

i and j, respectively as in Figure 2). We observe that TG is affected more by the

temperature at the point C, TC, than by that at the point I, TI because the point C is

closer to the point G than the point I. Since the closer the point, the greater the

influence, based on linear variation of the temperature we can then write:

TG =

?y

???y

1

??

TC +

TI =

TC +

T

?y + ???y

?y + ??y

1 + ??

1 + ?? I

(15)

where the increment from point I to point G is the regular increment y while that

between G and C is less and equal to y with having a value between 0 and

1. Similarly, along the vertical, x direction and considering the points B, G, and H

and a regular increment x, we can have:

TG =

© 2001 by CRC Press LLC

?x

? ??x

1

??

TB +

TH =

TB +

T

?x + ? ??x

?x + ? ??x

1 + ??

1 + ?? H

(16)

FIGURE 4. There are cases where the domain of heat conduction have boundaries which

are quite irregular geometrically.

where like , has a value between 0 and 1. As often is the case, the regular

increments x and y are taken to be equal to each other for the simplicity of

computation. Equations 15 and 16 can then be combined and by taking both x and

y directions into consideration, an averaging approach leads to:

TG =

?

1? 1

??

1

??

TB +

TH +

TC +

TI ?

?

2 ? 1 + ??

1 + ??

1 + ??

1 + ?? ?

(17)

For every group of five points such as B, C, G, H, and I in Figure 4 situated at

any irregular boundary, the values of and have to be measured and Equation

17 is to be used during the relaxation process if the boundary temperatures are

known.

If some points along an irregular boundary are insulated such as the points B

and C in Figure 4, we need to derive new formula to replace Equation 6 or Equation

10. The insulated condition along BC requires Tn = 0 where n is the direction

normal to the cord BC when the arc BC is approximated linearly. If the values of ‘

and ‘ are known, we can replace the condition T/n = 0 with

?T ?TdX ?TdY ?T

?T

?T

?T

=

+

=

+ ??

=0

sin ? +

cos ? = ? ?

?n ?Xdn ?Ydn ?X

?Y

?X

?Y

(18)

The remainder of derivation is left as a homework problem.

MATLAB APPLICATION

A Relaxatn.m file can be created to perform interactive MATLAB operations

and generate plots of the temperature distributions during the course of relaxation.

This file may be prepared as follows:

© 2001 by CRC Press LLC

This file can be applied to solve the sample problem run by first specifying the

boundary temperatures described in Equation 4 to obtain an initial distribution by

entering the MATLAB instructions:

© 2001 by CRC Press LLC

The fprintf command enables a label be added, in which the format %3.0f

requests 3 columns be provided without the decimal point for printing the value of

NR, and \n requests that next printout should be started on a new line. The Relaxatn.m

can now be utilized to perform the relaxations. Let first perform one relaxation by

entering

>> NR = NR + 1; fprintf(‘Sweep # %3.0f \n’,NR), [D,T] = feval(A:Relaxatn’,T]

The resulting display of the error defined in Equation 13 and the second temperature distribution is:

© 2001 by CRC Press LLC

In case that we need to have the 30th temperature distribution by performing

29 consecutive relaxations, we enter:

>> for NR = 3:30; [D,T] = feval(A:Relaxatn',T]; end

>> fprintf('Sweep # %3.0f \n',NR),D,T

The resulting display of the error defined in Equation 13 and the 30th temperature

distribution is:

To obtain a plot of this temperature distribution after the initial temperature

distribution has been relaxed 29 times, with gridwork and title as shown in Figure 5a,

the interactive MATLAB instructions entered are:

>> V = 0:1:50; contour(T,V’), grid, title(‘* After 30 relaxations *’)

Notice that 51 contours having values 0 through 50 with an increment of 1

defined in the row matrix V. In Figure 5a, the contour having a value equal to 0 is

© 2001 by CRC Press LLC

FIGURE 5. After 38 relaxations (a) and after 137 relaxations (b).

© 2001 by CRC Press LLC

along the upper edge (Y = 10) and right edge (X = 20), the first curved contour of

the right has a value equal to 1, and the values of the contours are increased from

right to left until the point marked “5” which has a temperature equal to 50 is

reached. It should be noted that along the left edge, the uppermost point marked

“10” has a temperature equal to zero and the temperatures are increased linearly (as

for the initial conditions) to 50 at the point marked with “5”, and from that point

down to the point marked “1” the entire lower portion of the left edge is insulated.

For obtaining the 137th sweep, we can continue to call the service Relaxatn.m

by similarly applying the MATLAB instructions as follows:

>> for NR = 31:137; [D,T] = feval(A:Relaxatn',T]; end

>> fprintf('Sweep # %3.0f \n',NR),D,T

The resulting display of the error defined in Equation 13 and the 137th temperature distribution is:

© 2001 by CRC Press LLC

Figure 5b shows the 137th temperature distribution when the interactive MATLAB instructions entered are:

>> V = 0:1:50; contour(T,V’), grid, title(‘* After 137 relaxations *’)

Notice that area near the insulated boundaries at the right-lower corner has finally

reached a steady-state temperature distribution, i.e., changes of the entire temperature

distribution will be insignificant if more relaxations were pursued.

MATHEMATICA APPLICATIONS

To apply the relaxation method for finding the steady-state temperature distribution of the heated plate already solved by the FORTRAN, QuickBASIC, and

MATLAB versions, here we make use of the Do, If, and While commands of

Mathematica to generate similar results through the following interactive operations:

In[1]: = t = Table[0,{10},{20}]; eps = 100; nr = 0; d = eps + 1;

In[2]: = Do[t[[i,1]] = (I–1)*10,{i,2,6}];

In[3]: = (While[d>eps, d = 0;nr = nr + 1;

Do[Do[ts = t[[i,j]]; t[[i,j]] = .25*(t[[I–1,j]] + t[[I + 1,j]]

+ t[[i,j–1]] + t[[i,j + 1]]);

d = Abs[ts-t[[i,j]]] + d;,{j,2,19}],{i,2,6}];

Do[ts = t[[i,1]]; t[[i,1]] = .25*(t[[I–1,1]] + t[[I + 1,1]]

+ 2*t[[i,2]]); d = Abs[ts-t[[i,1]]] + d;

Do[ts = t[[i,j]]; t[[i,j]] = .25*(t[[I–1,j]] + t[[I + 1,j]]

+ t[[i,j–1]] + t[[i,j + 1]]);

d = Abs[ts-t[[i,j]]] + d;,{j,2,19}];

If[i = = 7,Continue, ts = t[[i,20]];

t[[i,20]] = .25*(t[[I–1,20]] + t[[I + 1,20]] + 2*t[[i,19]]);

d = Abs[ts-t[[i,20]]] + d;],{i,7,9}];

ts = t[[10,1]]; t[[10,1]] = .5*(t[[9,1]] + t[[10,2]]);

d = Abs[ts-t[[10,1]]] + d;

Do[ts = t[[10,j]]; t[[10,j]] = .25*(2*t[[9,j]] + t[[10,j–1]]

+ t[[10,j + 1]]); d = Abs[ts-t[[10,j]]] + d;,{j,2,19}];

ts = t[[10,20]]; t[[10,20]] = .5*(t[[9,20]] + t[[10,19]]);

d = Abs[ts-t[[10,20]]] + d;])

In[4]: = Print[“Sweep #”,nr]; Round[N[t,2]]

Out[4] = Sweep #2

{{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{10, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{20, 9, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

© 2001 by CRC Press LLC

{30, 13, 5, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{40, 18, 7, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{50, 20, 8, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{17, 11, 5, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{ 6, 5, 3, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{ 2, 2, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}

Notice that In[2] initializes the boundary temperatures, nr keeps the count of

how many sweeps have been performed, and the function Round is employed in

In[4] to round the temperature value to a two-digit integer. When the total temperature differences, d, is limited to eps = 100 degrees, the t values obtained after two

sweeps are slightly different from those obtained by the other versions, this is again

because Mathematica keeps more significant digits in all computation steps than

those in the FORTRAN, QuickBASIC, and MATLAB programs. By changing the

eps value from 100 degrees to 1 degree, Mathematica also takes 137 sweeps to

converge as in the FORTRAN, QuickBASIC, and MATLAB versions:

In[5]: = t = Table[0,{10},{20}]; eps = 1; nr = 0; d = eps + 1;

In[6]: = Do[t[[i,1]] = (I–1)*10,{i,2,6}];

In[7]: = Print[“Sweep #”,nr]; Round[N[t,2]]

Out[7] = Sweep #137

{{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{10, 8, 7, 5, 4, 3, 3, 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0},

{20, 16, 13, 10, 8, 7, 5, 4, 3, 3, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0},

{30, 24, 19, 15, 12, 9, 8, 6, 5, 4, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0},

{40, 30, 23, 18, 15, 12, 10, 8, 6, 5, 4, 4, 3, 2, 2, 1, 1, 1, 0, 0},

{50, 34, 26, 20, 16, 13, 11, 9, 7, 6, 5, 4, 3, 3, 2, 2, 1, 1, 0, 0},

{35, 30, 25, 21, 17, 15, 12, 10, 8, 7, 6, 5, 4, 3, 3, 2, 2, 1, 1, 0},

{29, 27, 24, 21, 18, 15, 13, 11, 9, 7, 6, 5, 4, 3, 3, 2, 2, 1, 1, 1},

{26, 26, 23, 21, 18, 16, 13, 11, 9, 8, 6, 5, 4, 4, 3, 2, 2, 1, 1, 1},

{26, 25, 23, 21, 18, 16, 13, 11, 9, 8, 7, 5, 4, 4, 3, 2, 2, 2, 1, 1}}

The above results all agree with those obtained earlier.

WARPING ANALYSIS

OF A

TWISTED BAR

WITH

NONCIRCULAR CROSS SECTION

As another example of applying the relaxation method for engineering analysis,

consider the case of a long bar of uniform rectangular cross section twisted by the

two equal torques (T) at its ends. The cross section of the twisted bar becomes

warped as shown in Figure 6. If z-axis is assigned to the longitudinal direction of

the bar, to find the amount of warping at any point (x,y) of the cross sectioned

surface, denoted as W(x,y), the relaxation method can again be employed because

© 2001 by CRC Press LLC

FIGURE 6. In the case of a bar of uniform rectangular cross section twisted by the two

equal torques (T) at its ends, the cross section of the twisted bar becomes warped as shown.

W(x,y) is governed by the Laplace equation.3 Due to anti-symmetry of W(x,y), that

is W(–x,y) = W(x,–y) = –W(x,y), only one-fourth of the cross section needs to be

analyzed. Let us consider a rectangular region 0?x?a and 0?y?b. It is obvious that

the anti-symmetry leads to the conditions W(x,0) = 0 and W(0,y) = 0 along two of

the four linear boundaries. To derive the boundary conditions along the right side

x = a and the upper side y = b, we have to utilize the relationships between the

warping function W(x,y) and shear stresses

xz and

yz which are:

? ?W

?

? xz = G??

? y?

? ?x

?

and

?

? ?W

? yz = G??

+ x?

?

? ?y

(19,20)

where G and are the shear rigidity and twisting angle of the bar, respectively.

Along the boundaries x = a and y = b, the shear stress should be tangential to the

lateral surface of the twisted bar requires:

? dx

? ?W

? dy ? ?W

? y?

??

=0

+ x?

?

? ?x

? ds ? ?y

? ds

(21)

where s is the variable changed along the boundary. Along the upper boundary y =

b, dy/ds = 0 and dx/ds = 1 and along the right boundary x = a, dx/ds = 0 and dy/ds =

1. Consequently, Equation 19 reduces to:

?W

= ?x

?y

for y = b

and

?W

=y

x

for x = a

(22,23)

To apply the relaxation method for solving W(x,y) which satisfies the Laplace

equation for 0

Partial Differential

Equations

8.1 INTRODUCTION

Different engineering disciplines solve different types of problems in their respective

fields. For mechanical engineers, they may need to solve the temperature change

within a solid when it is heated by the interior heat sources or due to a rise or

decrease of its boundary temperatures. For electrical engineers, they may need to

find the voltages at all circuit joints of a computer chip board. Temperature and

voltage are the variables in their respective fields. Hence, they are called field

variables. It is easy to understand that the value of the field variable is spacedependent and time-dependent. That is to say, that we are interested to know the

spatial and temporal changes of the field variable. Let us denote the field variable

as . and let the independent variables be xi which could be the time t, or, the space

coordinates as x, y, and z. In order not to overly complicate the discussion, we

introduce the general two-dimensional partial differential equation which governs

the field variable in the form of:

A(x1, x 2 )

?2 ?

?2 ?

?2 ?

+ B(x1, x 2 )

+ C(x1, x 2 ) 2

2

?x1

?x1?x 2

?x 2

?

?? ?? ?

= F? x1, x 2 , ?,

,

?x1 ?x 2 ??

?

(1)

where the coefficient functions A, B, and C in the general cases are dependent on

the variables x1 and x2, and the right-hand-side function F, called forcing function

may depend not only on the independent variables x1 and x2 but may also depend

on the first derivatives of . There are innumerable of feasible solutions for Equation

1. However, when the initial and/or boundary conditions are specified, only particular

solution(s) would then be found appropriate.

In this chapter, we will discuss three simple cases when A, B, and C are all

constants. The first case is a two-dimensional, steady-state heat conduction problem

involving temperature as the field variable and only the spatial distribution of the

temperature needs to be determined, Equation 1 is reduced to a parabolic partial

differential equation named after Poisson and Laplace when the forcing function F

is not equal to, or, equal to zero, respectively. This is a case when the coefficient

functions in Equation 1 are related by the condition B2–4AC0.

The reason that these problems are called parabolic, elliptical, and hyperbolic

is because their characteristic curves have such geometric features. Readers interested in exploring these features should refer to a textbook on partial differential

equations.

Details will be presented regarding how the forward, backward, and central

differences discussed in Chapter 4 are to be applied for approximating the first and

second derivative terms appearing in Equation 1. Repetitive algorithms can be

devised to facilitate programming for straight-forward computation of the spatial

and temporal changes of the field variable. Numerical examples are provided to

illustrate how these changes can be determined by use of either QuickBASIC,

FORTRAN, MATLAB, or, Mathematica programs.

Although explanation of the procedure for numerical solution of these three

types of problems is given only for the simple one- and two-dimensional cases, but

its extension to the higher dimension case is straight forward. For example, one may

attempt to solve the transient heat conduction problem of a thin plate by having two

space variables instead of one space variable for a long rod. The steady-state heat

conduction problem of a thin plate can be extended for the case of a three-dimensional solid, and the string vibration problem can be extended to a two-dimensional

membrane problem.

8.2 PROGRAM PARABPDE — NUMERICAL SOLUTION

OF PARABOLIC PARTIAL DIFFERENTIAL EQUATIONS

The program ParabPDE is designed for numerically solving engineering problems

governed by parabolic partial differential equation in the form of:

??

?2?

=a 2

?t

?x

(1)

and is a function of t and x and satisfies a certain set of supplementary conditions.

Equation 1 is called a parabolic partial differential equation. For example, could

be the temperature, T, of a longitudinal rod shown in Figure 1 and the parameter a

in Equation 1 could be equal to k/c where k, c, and are the thermal conductivity,

specific heat, and specific weight of the rod, respectively. To make the problem more

specific, the rod may have an initial temperature of 0°F throughout and it is completely insulated around its lateral surface and also at its right end. If its left end is

to be maintained at 100°F beginning at the time t = 0, then it is of interest to know

© 2001 by CRC Press LLC

FIGURE 1. could be the temperature, T, of a longitudinal rod.

how the temperatures along the entire length of the rod will be changing as the time

progresses. This is therefore a transient heat conduction problem. One would like

to know how long would it take to have the entire rod reaching a uniform temperature

of 100°F.

If the rod is made of a single material, k/c would then be equal to a constant.

Analytical solution can be found for this simple case.1 For the general case that the

rod may be composed of a number of different materials and the physical properties

k, c, and would not only depend on the spatial variable x but may also depend on

the temporal variable t. The more complicated the variation of these properties in x

and t, the more likely no analytical solution is possible and the problem can only

be solved numerically. The finite-difference approximation of Equation 1 can be

achieved by applying the forward difference for the first derivative with respective

to t and central difference for the second derivative with respect to x as follows (for

t at ti and x at xj):

?T Ti +1, j ? Ti, j

=?

?t

?t

and

?2 T Ti, j?1 ? 2 Ti, j + Ti, j+1

=?

?x 2

( ?x)2

If k/c is changing in time and also changing from one location to another, we

could designate it as ai,j. As a consequence, Equation 1 can then be written as:

Ti +1, j ? Ti, j

?t

= a i, j

Ti, j?1 ? 2 Ti, j + Ti, j+1

( ?x)2

(2)

Since the initial temperature distribution T is known, the above expression

suggests that for a numerical solution we may select an appropriate increment in t,

t, and the temperature be determined at a finite number of stations, N. It is advisable

to have these stations be equally spaced so that the increment x is equal to L/(N–1)

where L is the length of the rod, and the instants are to be designated as t1 = 0, t2 =

t,…, ti = (i–1)t, and the stations as x1 = 0, x2 = x,…, xj = (j–1)?x,…, and xN =

(N–1)x = L. The task at hand is then to find T(ti,xj) for i = 1,2,… and j = 1,2,…,N.

It can be noticed from Equation 2 that the there is only one temperature at ti + 1 and

can be expressed in terms of those at the preceding instant ti as:

© 2001 by CRC Press LLC

Ti +1, j = Ti, j +

a i, j?t

( ?x)2

(T

i , j?1

? 2 Ti, j + Ti, j+1

)

(3)

Equation 3 is to be used for j = 2 through j = N–1. For the last station, j = N,

which is insulated,the temperatures on both side of this station can be assumed to

be equal (the station N + 1 is a fictitious one!). The modified equation for this

particular station is:

Ti +1,N = Ti,N +

2a i,N ?t

( ?x)2

(T

i ,N ?1

? Ti,N

)

(4)

For generating the temperatures of the rod at N stations for any specified time

increment t until the temperatures are almost all equal to 100°F throughout, the

program ParabPDE has been applied. It is listed below along with a typical printout

of the results.

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Output

© 2001 by CRC Press LLC

QUICKBASIC VERSION

MATLAB APPLICATIONS

A MATLAB version of ParabPDE can be created easily by converting the

QuickBASIC program. The m file may be arranged as follows:

© 2001 by CRC Press LLC

For solving the sample transient temperature problem, this m file can be called

and interactive MATLAB instructions can be entered through keyboard to obtain

the temperature distribution of the rod at various times:

© 2001 by CRC Press LLC

FIGURE 2. A composite graph with axes labels and markings of the curves by making use

of the MATLAB commands xlabel, ylabel, and text.

Notice that results of temperature distributions which have been terminated using

five differentials of 1, 0.1, 0.01, 0.001, and 0.0001 Fahrenheit are saved in Tsave

and then later plotted. Since T is a row matrix, the transpose of T, T, is stored in

an appropriate column of Tsave. If the temperatures were rounded, the rod reaches

a uniform temperature distribution of 100°F when the required temperature differential is selected to be 0.001°F. If the fifth curve for the temperature differential

equal to 0.0001 is plotted, it will be too close to the forth curve and also the plot

function provides only four line types (solid, broken, dot, and center lines), the fifth

set of results is therefore not saved in Tsave. Figure 2 shows a composite graph with

axes labels and markings of the curves by making use of the MATLAB commands

xlabel, ylabel, and text.

MATHEMATICA APPLICATIONS

The heat conduction problem of an insulated rod previously discussed in the

versions for FORTRAN, QuickBASIC, and MATLAB can be solved by application

of Mathematica as follows:

© 2001 by CRC Press LLC

In[1]: = (n = 11; rk = 0.037; c = 0.212; rho = 168; dt = 1; dx = 0.1;

c1 = rk/c/rho*dt/dx^2; nm1 = n–1;)

In[2]: = (t = Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100; tdf = 1;

tm = 0; flag1 = 0;

In[3]: = (While[flag1 = = 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),

{j,2,nm1}];

tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);tm = tm + dt;flag1 = 1;

Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,

Continue],{i,n}];

Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm]; Print[N[tn,3]])

Out[3] = t = 24

{100., 65.7, 37.5, 18.5, 7.83, 2.85, 0.892, 0.241, 0.0561, 0.0116, 0.00394}

The temperature distribution of the heated rod after 24 seconds is same as

obtained by the FORTRAN, QuickBASIC, and MATLAB versions when every

component of two consecutive temperature distribution, kept as {t} and {tn}, differ

no more than the allowed set value of tdf = 1 degree in In[2]. Any component has

a difference exceeding the tdf value will cause flag1 to change from a value of 1 to

0 and the Break command in the second Do loop in In[2] to exit and to continue

the iteration. flag1 is created to control the While command which determines when

the iteration should be terminated. The N[tn,3] instructs the components of {tn} be

printed with 3 significant figures.

When tdf is changed to a value of 0.5, Mathematica can again be applied to yield

In[4]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;

tdf = 0.5; tm = 0; flag1 = 0;)

In[5]: = (While[flag1 = = 0,Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),

{j,2,nm1}];

tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);tm = tm + dt;flag1 = 1;

Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,

Continue],{i,n}];

Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm]; Print[N[tn,3]])

Out[5]: = t = 49

{100., 75.5, 53.2, 34.9, 21.3, 12., 6.24, 3., 1.35, 0.617, 0.413}

Notice that Mathematica takes only 49 seconds, one second less than that

required for the FORTRAN, QuickBASIC, and MATLAB versions. The reason is

that Mathematica keeps more significant digits in carrying out all computations.

To show more on the effect of changing the tdf value, the following Mathematica

runs are provided:

© 2001 by CRC Press LLC

In[6]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;

tdf = 0.1; tm = 0; flag1 = 0;)

In[7]: = (While[flag1 == 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),

{j,2,nm1}];

tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);

tm = tm + dt; flag1 = 1;

Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,

Continue],{i,n}];

Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm];

Print[N[tn,3]])

Out[7]: = t = 462

{100., 93.9, 88., 82.3, 77.1, 72.5, 68.5, 65.3, 63., 61.6, 61.1}

In[8]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;

tdf = 0.05; tm = 0; flag1 = 0;)

In[9]: = (While[flag1 == 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),

{j,2,nm1}];

tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);

tm = tm + dt; flag1 = 1;

Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,

Continue],{i,n}];

Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm];

Print[N[tn,3]])

Out[9]: = t = 732

{100., 97., 94., 91.2, 88.5, 86.2, 84.2, 82.6, 81.5, 80.8, 80.5}

In[10]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;

tdf = 0.0001; tm = 0; flag1 = 0;

In[11]: = (While[flag1 == 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),

{j,2,nm1}];

tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);

tm = tm + dt; flag1 = 1;

Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,

Continue],{i,n}];

Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm]; Print[N[tn,3]])

Out[11]: = t = 3159

{100., 100., 100., 100., 100., 100., 100., 100., 100., 100., 100.}

For tdf = 0.1 and tdf = 0.05, Mathematica continues to take lesser time than the

other version; but when tdf = 0.0001, Mathematica needs two additional seconds. The

© 2001 by CRC Press LLC

reason is that seven significant figures are required in the last case, rounding may

have resulted in earlier termination of the iteration when the FORTRAN, QuickBASIC, and MATLAB versions are employed.

8.3 PROGRAM RELAXATN — SOLVING ELLIPTICAL PARTIAL

DIFFERENTIAL EQUATIONS BY RELAXATION METHOD

The program Relaxatn is designed for solving engineering problems which are

governed by elliptical partial differential equation of the form:

?2 ? ?2 ?

+

= F(x, y)

?x 2 ?y 2

(1)

where is called the field function and F(x,y) is called forcing function. When the

steady-state heat conduction of a two-dimensional domain is considered,2 then the

field function becomes the temperature distribution, T(x,y), and the forcing function

becomes the heat-source function, Q(x,y). If the distribution of T is influenced only

by the temperatures at the boundary of the domain, then Q(x,y) = 0 and Equation

1 which is often called a Poisson equation is reduced to a Laplace equation:

?2 T ?2 T

+

=0

?x 2 ?y 2

(2)

The second-order, central-difference formulas (read the program DiffTabl) can

be applied to approximate the second derivatives in the above equation at an arbitrary

point x = xi and y = yj in the domain as:

T ? 2 Ti, j + Ti +1, j

?2 T

at x i , y j =? i ?1, j

2

?x

( ?x)2

and

T ? 2 Ti, j + Ti, j+1

?2 T

at x i , y j =? i, j?1

2

?y

(?y)2

By substituting both of the above equations into Equation 2 and taking equal

increments in both x- and y-directions, the reduced equation is, for x = y,

Ti, j =

(

1

T + Ti +1, j + Ti, j?1 + Ti, j+1

4 i ?1, j

)

(3)

The result is expected when the temperature distribution reaches a steady state

because it states that the temperature at any point should be equal to the average of

its surrounding temperatures.

© 2001 by CRC Press LLC

FIGURE 3. A plate, which initially (at time t = 0) had a temperature equal to 0°F throughout,

insulated along a portion of its boundary and suddenly heated at its upper left boundary to

maintain a linearly varying temperature Tb.

Before we proceed further, it is appropriate at this time to introduce a numerical

case. Suppose that a plate which initially (at time t = 0) has a temperature equal to

0°F throughout and is insulated along a portion of its boundary, is suddenly heated

at its upper left boundary to maintain a linearly varying temperature Tb as shown in

Figure 3. If this heating process is to be maintained, we are then interested in

knowing the temperature distribution, changed from its initial state of uniformly

equal to 0°F if given sufficient time to allow it to reach an equilibrium (steady) state.

Numerically, we intend to calculate the temperatures, denoted as a matrix [T], at a

selected number of locations. Therefore, the plate is first divided into a gridwork of

M rows and N columns along the x- and y-directions, respectively, as indicated in

Figure 1. The directions of x- and y-axes are so selected for the convenience of

associating them with the row and column of the temperature matrix [T] which is

of order M by N. The values of M and N should be so decided such that the

increments x and y are equal in order to apply Equation 3. To be more specific,

let M = 10 and N = 20 and the linear temperature variation along the upper left

boundary Tb be:

Ti,1 = 10(i ? 1)

for

i = 1, 2,…, 6

(4)

Here, Ti,j is to be understood as the temperature at the location (xi,yj). Equation

4 describes the temperature along the left boundary y = y1 but only for xi in the

range of i = 1 to i = 6.

Since the temperatures at the stations which are on the insulated boundaries of

the plate are also involved, these unknown temperatures need to be treated differently.

By an insulated boundary, it means that there is no heat transfer normal to the

boundary. Since the heat flow is depended on the temperature difference across that

© 2001 by CRC Press LLC

insulated boundary, mathematically it requires that T/n = 0 there, n being the

normal direction. At the vertical boundaries x = xM, we have T/n = T/x = 0 since

x is the normal direction. Based on the central difference and considering two

increments in the x direction, at a yjth station we can have:

T

? TM ?1, j

?T

at x M , y j =? M +1, j

=0

?x

2( ?x)

(

)

or

TM +1, j = TM ?1, j

(5)

Since xM + 1 is below the bottom boundary of the plate shown in Figure 3, there

is no need to calculate the temperatures there, however Equation 5 enables the

temperatures at the stations along the bottom boundary of the plate TM,j for j = 1 to

N to be averaged. Returning to Equation 4, we notice that if Equation 5 is substituted

into it, the resulting equation which relates only to three neighboring temperatures is:

TM, j =

(

1

2 TM ?1, j + TM, j?1 + TM, j+1

4

)

for j = 2, 3,…, N ? 1

(6)

Notice that j = 1 and j = N are not covered in Equation 6. These two cases

concerning the insulated stations at the left and right bottom corners of the plate

will be discussed after we address the two vertical, insulated boundaries y = y1 and

y = yN.

For the boundaries y = y1, ?T/?n becomes ?T/?y. Again, we can apply the central

difference for double y increments to obtain:

T ? Ti,2

?T

at (x1, y1 ) =? i,0

=0

?y

2( ?y)

or

Ti,0 = Ti,2

(7)

Thus, the modified Equation 3 for the left insulated boundary is:

Ti,1 =

(

1

T + Ti +1,1 + 2 Ti,2

4 i ?1,1

)

for

i = 7, 8, 9

In a similar manner, we can derive for the right insulated boundary y = yN

T

? Ti,N ?1

?T

at (x i , y N ) =? i,N +1

=0

?y

2( ?y)

or

© 2001 by CRC Press LLC

(8)

Ti,N +1 = Ti,N ?1

(9)

and

Ti,N =

(

1

T

+ Ti +1,N + 2 Ti,N ?1

4 i ?1,N

)

for i = 8, 9

(10)

Having derived Equations 6, 8, and 10, it is easy to deduce the two special

equation for the corner insulated stations to be:

TM,1 =

(

1

T + TM ?1,1

2 M ,2

)

(11)

and

TM,N =

(

1

T

+ TM ?1,N

2 M,N ?1

)

(12)

We have derived all equations needed for averaging the temperature at any station

of interest including those at the insulated boundaries by utilizing those at its

neighboring stations. It suggests that a continuous upgrading process can be developed which assumes that the neighboring temperatures are known. This so-called

relaxation method starts with an initial assumed distribution of temperature [T(0)]

and continues to use Equations 3 and 6 to 12 until the differences at all locations

are small enough. Mathematically, the process terminates when:

M

N

i =1

j=1

?? T

( k +1)

i, j

? Ti(,kj ) < ?

(13)

where is a prescribed tolerance of accuracy and k = 0,1,2,… is the number of

sweeps in upgrading the temperature distribution. Superscripts (k + 1) and (k) refer

to the improved and previous distributions, respectively. The order of sweep will

affect how the temperatures should be upgraded. For example, if the temperatures

are to be re-averaged from top to bottom and left to right, referring to Figure 1, then

Equation 3 is to be modified as:

Ti(,kj +1) =

(

1 ( k +1)

T

+ Ti(+k1), j + Ti(,kj?+11) + Ti(,kj+)1

4 i ?1, j

)

(14)

Notice that the neighboring temperatures in the row above, i–1, and in the column

to the left, j–1, have already been upgraded while those in the row below, i + 1, and

in the column to the right, j + 1, are yet to be upgraded. Similar modifications are

to be made to Equations 6 to 12 during relaxation.

© 2001 by CRC Press LLC

The program Relaxatn is developed according to the relaxation method

described above. For solving the problem shown in Figure 3, both FORTRAN and

QuickBASIC versions are made available for interactively specifying the tolerance.

Sample results are presented below.

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Results

The program Relaxatn is first applied for an interactively entered value of equal

to 100. Only one relaxation needs to be implemented as shown below. The temperature distribution for Sweep #1 is actually the initial assumed distribution. One

cannot assess how accurate this distribution is. The second run specifies that be

© 2001 by CRC Press LLC

equal to 1. The results show that 136 relaxation steps are required. For giving more

insight on how the relaxation has proceeded, Sweeps #10, #30, #50, #100, and #137

are presented for interested readers. It clearly indicates that a tolerance of equal to

100 is definitely inadequate.

© 2001 by CRC Press LLC

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Results

Irregular Boundaries

Practically, there are cases where the domain of heat conduction have boundaries

which are quite irregular geometrically as illustrated in Figure 4. For such cases, the

equation derived based on the relaxation method, Equation 3, which states that the

temperature at any point has the average value of those at its four neighboring points

if they are equally apart, has to be modified. The modified equation can be derived

using a simple argument applied in both x and y directions. For example, consider

the temperature at the point G, TG, in Figure 4. First, let us investigate the horizontal,

y direction (for convenience of associating x and y with the row and column indices

i and j, respectively as in Figure 2). We observe that TG is affected more by the

temperature at the point C, TC, than by that at the point I, TI because the point C is

closer to the point G than the point I. Since the closer the point, the greater the

influence, based on linear variation of the temperature we can then write:

TG =

?y

???y

1

??

TC +

TI =

TC +

T

?y + ???y

?y + ??y

1 + ??

1 + ?? I

(15)

where the increment from point I to point G is the regular increment y while that

between G and C is less and equal to y with having a value between 0 and

1. Similarly, along the vertical, x direction and considering the points B, G, and H

and a regular increment x, we can have:

TG =

© 2001 by CRC Press LLC

?x

? ??x

1

??

TB +

TH =

TB +

T

?x + ? ??x

?x + ? ??x

1 + ??

1 + ?? H

(16)

FIGURE 4. There are cases where the domain of heat conduction have boundaries which

are quite irregular geometrically.

where like , has a value between 0 and 1. As often is the case, the regular

increments x and y are taken to be equal to each other for the simplicity of

computation. Equations 15 and 16 can then be combined and by taking both x and

y directions into consideration, an averaging approach leads to:

TG =

?

1? 1

??

1

??

TB +

TH +

TC +

TI ?

?

2 ? 1 + ??

1 + ??

1 + ??

1 + ?? ?

(17)

For every group of five points such as B, C, G, H, and I in Figure 4 situated at

any irregular boundary, the values of and have to be measured and Equation

17 is to be used during the relaxation process if the boundary temperatures are

known.

If some points along an irregular boundary are insulated such as the points B

and C in Figure 4, we need to derive new formula to replace Equation 6 or Equation

10. The insulated condition along BC requires Tn = 0 where n is the direction

normal to the cord BC when the arc BC is approximated linearly. If the values of ‘

and ‘ are known, we can replace the condition T/n = 0 with

?T ?TdX ?TdY ?T

?T

?T

?T

=

+

=

+ ??

=0

sin ? +

cos ? = ? ?

?n ?Xdn ?Ydn ?X

?Y

?X

?Y

(18)

The remainder of derivation is left as a homework problem.

MATLAB APPLICATION

A Relaxatn.m file can be created to perform interactive MATLAB operations

and generate plots of the temperature distributions during the course of relaxation.

This file may be prepared as follows:

© 2001 by CRC Press LLC

This file can be applied to solve the sample problem run by first specifying the

boundary temperatures described in Equation 4 to obtain an initial distribution by

entering the MATLAB instructions:

© 2001 by CRC Press LLC

The fprintf command enables a label be added, in which the format %3.0f

requests 3 columns be provided without the decimal point for printing the value of

NR, and \n requests that next printout should be started on a new line. The Relaxatn.m

can now be utilized to perform the relaxations. Let first perform one relaxation by

entering

>> NR = NR + 1; fprintf(‘Sweep # %3.0f \n’,NR), [D,T] = feval(A:Relaxatn’,T]

The resulting display of the error defined in Equation 13 and the second temperature distribution is:

© 2001 by CRC Press LLC

In case that we need to have the 30th temperature distribution by performing

29 consecutive relaxations, we enter:

>> for NR = 3:30; [D,T] = feval(A:Relaxatn',T]; end

>> fprintf('Sweep # %3.0f \n',NR),D,T

The resulting display of the error defined in Equation 13 and the 30th temperature

distribution is:

To obtain a plot of this temperature distribution after the initial temperature

distribution has been relaxed 29 times, with gridwork and title as shown in Figure 5a,

the interactive MATLAB instructions entered are:

>> V = 0:1:50; contour(T,V’), grid, title(‘* After 30 relaxations *’)

Notice that 51 contours having values 0 through 50 with an increment of 1

defined in the row matrix V. In Figure 5a, the contour having a value equal to 0 is

© 2001 by CRC Press LLC

FIGURE 5. After 38 relaxations (a) and after 137 relaxations (b).

© 2001 by CRC Press LLC

along the upper edge (Y = 10) and right edge (X = 20), the first curved contour of

the right has a value equal to 1, and the values of the contours are increased from

right to left until the point marked “5” which has a temperature equal to 50 is

reached. It should be noted that along the left edge, the uppermost point marked

“10” has a temperature equal to zero and the temperatures are increased linearly (as

for the initial conditions) to 50 at the point marked with “5”, and from that point

down to the point marked “1” the entire lower portion of the left edge is insulated.

For obtaining the 137th sweep, we can continue to call the service Relaxatn.m

by similarly applying the MATLAB instructions as follows:

>> for NR = 31:137; [D,T] = feval(A:Relaxatn',T]; end

>> fprintf('Sweep # %3.0f \n',NR),D,T

The resulting display of the error defined in Equation 13 and the 137th temperature distribution is:

© 2001 by CRC Press LLC

Figure 5b shows the 137th temperature distribution when the interactive MATLAB instructions entered are:

>> V = 0:1:50; contour(T,V’), grid, title(‘* After 137 relaxations *’)

Notice that area near the insulated boundaries at the right-lower corner has finally

reached a steady-state temperature distribution, i.e., changes of the entire temperature

distribution will be insignificant if more relaxations were pursued.

MATHEMATICA APPLICATIONS

To apply the relaxation method for finding the steady-state temperature distribution of the heated plate already solved by the FORTRAN, QuickBASIC, and

MATLAB versions, here we make use of the Do, If, and While commands of

Mathematica to generate similar results through the following interactive operations:

In[1]: = t = Table[0,{10},{20}]; eps = 100; nr = 0; d = eps + 1;

In[2]: = Do[t[[i,1]] = (I–1)*10,{i,2,6}];

In[3]: = (While[d>eps, d = 0;nr = nr + 1;

Do[Do[ts = t[[i,j]]; t[[i,j]] = .25*(t[[I–1,j]] + t[[I + 1,j]]

+ t[[i,j–1]] + t[[i,j + 1]]);

d = Abs[ts-t[[i,j]]] + d;,{j,2,19}],{i,2,6}];

Do[ts = t[[i,1]]; t[[i,1]] = .25*(t[[I–1,1]] + t[[I + 1,1]]

+ 2*t[[i,2]]); d = Abs[ts-t[[i,1]]] + d;

Do[ts = t[[i,j]]; t[[i,j]] = .25*(t[[I–1,j]] + t[[I + 1,j]]

+ t[[i,j–1]] + t[[i,j + 1]]);

d = Abs[ts-t[[i,j]]] + d;,{j,2,19}];

If[i = = 7,Continue, ts = t[[i,20]];

t[[i,20]] = .25*(t[[I–1,20]] + t[[I + 1,20]] + 2*t[[i,19]]);

d = Abs[ts-t[[i,20]]] + d;],{i,7,9}];

ts = t[[10,1]]; t[[10,1]] = .5*(t[[9,1]] + t[[10,2]]);

d = Abs[ts-t[[10,1]]] + d;

Do[ts = t[[10,j]]; t[[10,j]] = .25*(2*t[[9,j]] + t[[10,j–1]]

+ t[[10,j + 1]]); d = Abs[ts-t[[10,j]]] + d;,{j,2,19}];

ts = t[[10,20]]; t[[10,20]] = .5*(t[[9,20]] + t[[10,19]]);

d = Abs[ts-t[[10,20]]] + d;])

In[4]: = Print[“Sweep #”,nr]; Round[N[t,2]]

Out[4] = Sweep #2

{{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{10, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{20, 9, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

© 2001 by CRC Press LLC

{30, 13, 5, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{40, 18, 7, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{50, 20, 8, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{17, 11, 5, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{ 6, 5, 3, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{ 2, 2, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}

Notice that In[2] initializes the boundary temperatures, nr keeps the count of

how many sweeps have been performed, and the function Round is employed in

In[4] to round the temperature value to a two-digit integer. When the total temperature differences, d, is limited to eps = 100 degrees, the t values obtained after two

sweeps are slightly different from those obtained by the other versions, this is again

because Mathematica keeps more significant digits in all computation steps than

those in the FORTRAN, QuickBASIC, and MATLAB programs. By changing the

eps value from 100 degrees to 1 degree, Mathematica also takes 137 sweeps to

converge as in the FORTRAN, QuickBASIC, and MATLAB versions:

In[5]: = t = Table[0,{10},{20}]; eps = 1; nr = 0; d = eps + 1;

In[6]: = Do[t[[i,1]] = (I–1)*10,{i,2,6}];

In[7]: = Print[“Sweep #”,nr]; Round[N[t,2]]

Out[7] = Sweep #137

{{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

{10, 8, 7, 5, 4, 3, 3, 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0},

{20, 16, 13, 10, 8, 7, 5, 4, 3, 3, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0},

{30, 24, 19, 15, 12, 9, 8, 6, 5, 4, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0},

{40, 30, 23, 18, 15, 12, 10, 8, 6, 5, 4, 4, 3, 2, 2, 1, 1, 1, 0, 0},

{50, 34, 26, 20, 16, 13, 11, 9, 7, 6, 5, 4, 3, 3, 2, 2, 1, 1, 0, 0},

{35, 30, 25, 21, 17, 15, 12, 10, 8, 7, 6, 5, 4, 3, 3, 2, 2, 1, 1, 0},

{29, 27, 24, 21, 18, 15, 13, 11, 9, 7, 6, 5, 4, 3, 3, 2, 2, 1, 1, 1},

{26, 26, 23, 21, 18, 16, 13, 11, 9, 8, 6, 5, 4, 4, 3, 2, 2, 1, 1, 1},

{26, 25, 23, 21, 18, 16, 13, 11, 9, 8, 7, 5, 4, 4, 3, 2, 2, 2, 1, 1}}

The above results all agree with those obtained earlier.

WARPING ANALYSIS

OF A

TWISTED BAR

WITH

NONCIRCULAR CROSS SECTION

As another example of applying the relaxation method for engineering analysis,

consider the case of a long bar of uniform rectangular cross section twisted by the

two equal torques (T) at its ends. The cross section of the twisted bar becomes

warped as shown in Figure 6. If z-axis is assigned to the longitudinal direction of

the bar, to find the amount of warping at any point (x,y) of the cross sectioned

surface, denoted as W(x,y), the relaxation method can again be employed because

© 2001 by CRC Press LLC

FIGURE 6. In the case of a bar of uniform rectangular cross section twisted by the two

equal torques (T) at its ends, the cross section of the twisted bar becomes warped as shown.

W(x,y) is governed by the Laplace equation.3 Due to anti-symmetry of W(x,y), that

is W(–x,y) = W(x,–y) = –W(x,y), only one-fourth of the cross section needs to be

analyzed. Let us consider a rectangular region 0?x?a and 0?y?b. It is obvious that

the anti-symmetry leads to the conditions W(x,0) = 0 and W(0,y) = 0 along two of

the four linear boundaries. To derive the boundary conditions along the right side

x = a and the upper side y = b, we have to utilize the relationships between the

warping function W(x,y) and shear stresses

xz and

yz which are:

? ?W

?

? xz = G??

? y?

? ?x

?

and

?

? ?W

? yz = G??

+ x?

?

? ?y

(19,20)

where G and are the shear rigidity and twisting angle of the bar, respectively.

Along the boundaries x = a and y = b, the shear stress should be tangential to the

lateral surface of the twisted bar requires:

? dx

? ?W

? dy ? ?W

? y?

??

=0

+ x?

?

? ?x

? ds ? ?y

? ds

(21)

where s is the variable changed along the boundary. Along the upper boundary y =

b, dy/ds = 0 and dx/ds = 1 and along the right boundary x = a, dx/ds = 0 and dy/ds =

1. Consequently, Equation 19 reduces to:

?W

= ?x

?y

for y = b

and

?W

=y

x

for x = a

(22,23)

To apply the relaxation method for solving W(x,y) which satisfies the Laplace

equation for 0