A finite-difference solution to a two-dimensional vibrating membrane problem using a spreadsheet program

Advances in Engineering Software 23 (1995) 115- 120

Copyright 0 1995 Elsevier ScienceLimited Printed in Great Britain. All rights reserved 0965-9978/95/$09,50

0965-9978(95)00024-O

ELSEVIER

Technical Note A finite-difference solution to a two-dimensional vibrating membrane problem using a spreadsheet program Abdelwahab Kharab Mathematics Department, King Fahd University of Petroleum & Minerals, Dhahran 31261, Saudi Arabia

(Received7 December1994;accepted1 May 1995) A complete 123macro that finds the numerical solution of a two-dimensional vibrating membraneproblemis described.The numericalmethodis basedon an explicit finite differenceapproximation. The 123 program allows initial and boundaryconditionsto be entereddirectly into the spreadsheet and the numerical solution is obtainedby pressinga singlekey. This important featurewill allow the user to easily experimentthe model problem with various initial and boundary conditions. solutions consisted of a number of steps which the user had to execute, our approach proposes a solution which is executed in a single 123 macro program. This would help to give students a better understanding of the numerical method and leave them free to focus on the unfolding of the solution.

1 INTRODUCTION

The wave equation has many applications in scienceand engineering. It applies to problems in such areas as vibrations, electrostatistics, gas dynamics and acoustics. The present article describes a method for using spreadsheets to find the numerical solution of a twodimensional wave equation that applies to a vibration of a square membrane. Spreadsheets have been a power tool for business calculations, but they can also be a powerful tool for tackling many engineering problems. More recently, spreadsheetprograms run on microcomputers have been widely used to find the numerical solution of many problems. Spreadsheetsprograms can provide rapid and simple numerical solutions for boundary value problems in one or two dimensions.“* A method for using spreadsheets to perform a numerical analysis of a steady-state conduction is presented in Eid.3 Kerkes4 used it to obtain the numerical solution for transient fluid flow in a porous medium with a moving drainage boundary. In engineering they have been used for analysing logical networks’ and solving differential equations.’ Crow6 describes the spreadsheet solution to Laplace’s equation using a finite difference method. Unlike the approach adopted in Refs. 1, 3-6 where the

2 THE MODEL

PROBLEM

Consider a thin square membrane of mass per unit area p stretched by a tensile force per unit length of T as shown in Fig. 1.

Fig. 1. Physicalmodel. 115

116

A. Kharab

Assuming that all motion is at right angles to the xy plane, the membrane is tightly stretched so that the force T can be regarded as constant and there is no external force applied to the membrane, the equation of the displacement u(x, y, t) is given by the two-dimensional wave equation a% a% @+@=&g.

1 a2u

OO.

(1)

(2) At the boundary of the square 0 5 x, y 5 a, the displacement is taken to be zero, that is @,Y, 4 = 0,

4a, Y, t) = 0,

t > 0,

u(x, 0, t) = 0,

u(x, a, t) = 0,

t > 0.

(3)

3 NUMERICAL METHOD AND SPREADSHEET DESCRIPTION The explicit form of the finite difference equations for the two-dimensional wave equation is described in many numerical analysis books.7s Development of a finite difference method for eqn (1) necessitatesthe introduction of a square mesh of spacing h in the square region 0 5 x, y 5 a. The nodal points are defined by i, j=O,l,...,

M

h = a/M.

We define by k the time-step size with t, = nk, n=0,1,2 ).... Approximating the space derivative by the central second difference and the time derivative by a forward difference at (xi, yj, tn), results in the difference equation U~j+l-2U~j+U~j-l h2 =

24~~’ -2qj

+4+l,j-2u2j+4-l,j

h2 + d/J1

k2

’

Solving eqn (4) for utf t , we obtain the explicit equation for ‘marching’ ahead of time !A:;’ = 2(1 - 2X2)ZJj +x2(4+l,j+utj+l

+4,j-1

+&l,j)

-l&l, i, j=

n+l ui, 0

fl+1

=UiM

_ -O,

j=O,l,...,

M,

n=O,l,...

=O,

i=O,l,...,

M,

n=O,l,... (6)

1,2 ,..., M-

1, n= 1,2,...

where r&j = U(xi,yj, tn) and X = cYk/h.

i, j=O,l,...,

M.

(7) The initial velocity condition (4) is replaced by the a central-difference, giving utj

=A,j,

UiTj’ = U!,j

Utj

=

-

2kgi,j,

i, j= 0,l ,.“, M.

(8)

(1 - 2X2)h,j +$f

Here a2 = T/p, a are constants andf(x) and g(x) are given functions.

with

n+l = uA4,j

By setting IZ= 0 in (5) and using (8) we get

0 < x,y 5 a.

U(X,Y,O) =f(x,y),

y=jh

n+l uO,j

and the initial condition (2) implies

The initial velocity and displacement are given by

Xi = ih,

The boundary conditions (3) give

(5)

i, j=

z+l,j +h,j+l 1,2 ,..., M-

+f;:,j-1 +h-l,j) 1

+kgiJ, (9)

where&j =f(xi,Yj) and gi,j = g(xi,Yj)We used an explicit method of the approximation to the governing eqn (1) because it involves a simple algorithm in the spreadsheet.The only restriction is that the mesh ratio X I l/a to ensure the stability of the numerical process. A spreadsheet is a two-dimensional grid, made up of individual cells. Each intersection of a row and a column makes up a cell. The spreadsheet structure is quite well suited for finite difference methods since each nodal point can be represented by a cell. The 123 macro that gives the numerical solution using the above method is shown in Fig. 2 with somecomments given on the left-hand side of the range where the macro is written. The program follows the same format as programs written in FORTRAN or BASIC languages. We will describe the main parts of the 123 macro. At the beginning, the program defines all variables to be used and sets the values of X (LAMDA), h (DX) and K. These commands are contained in cells ID4..IDll. There are five FOR loops in the program. Loop 1 and 2 are used to set the values of the initial displacement u(x, y, 0) =f(x, y). The functionf(x, y) is defined by the variable F(XY) and g(x, y) by the variable G(XY). At the end of the FOR loop 1, the values of u(x, y, 0) are obtained as shown in Fig. 3. The FOR loops 3, 4 and 5 are the ones used to compute u from eqns (5) and (9) (seeFig. 4). Loop 3 sets the counter for the time steps and loop 4 and 5 for the spacesvariables x and y. Equations (5) and (9) reside in cells IG58 and IG57 and are written using the 123 function “@INDEX(UW,J,I)” which plays the same role as arrays do in FORTRAN. Note that at each time step t = (n + l)k, u(x, y, (n + 1)k) is computed using its values at the two previous time steps, that is at t = nk and t = (n - 1)k. The values of u at t = (n - 1)k and t = nk are saved respectively in the ranges ‘US’ and ‘U’.

Finite-diflerence solution to a 2-D vibrating membraneproblem IH

A : 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 10 19 20 21 22 23 24 25 26 27 29 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 40 49 50 51 52 53 54 55 56 57 58 59 60

“&vi cxnioi es& ,&&

,

‘~‘OR’&op

:Kw!T) to ‘US’ i -4

. ,[G~T~JL+(~IGfjT [LtgT Y,SOX)JFfJF. I.O..M.l,!F32)-

II

.__

.M

IK

_ ,. fhi

J)

.: f=ii ‘W with 7.8R FOR !m? 4 ---’ ....Movear~to’u .. ~Ft~e~e,rang~ u\I :

If n=1 set yzlwx Mow culsor to ‘U FOR LOOP 5 --a

(IF N=l}(LET Y.J’DX)(GOTO)UW-(RIGHT J) (FOR I.l.M-1 ,l,lG53)(RETURN)

Move C”c3of dowll If l-i=1 set y&IX

Compute Compute

u from (9) u from (5)

WW (IF N=l)(LET (IF N-l)(w (IF N>l)&Et’ (IF N=l)(LET (IF N>l)(LET +IMl-IN-(RETURN)

X.I’DX)I~OINQEX(USJ.I-t)+0lNDU((US,J~+l)+OINDEX(US,J+1,I)+0tNDEX(US,J-1,I))IM2.O~DM(U,J.I-1)+OINDU((U,J,l+l)+OI~(U.J+1.I)+OINDU((U.J-l,I))IMl.(l-~D~6lNbW(USJ,l)~A~~+~Q(~)IM1,2’(1-2~~MD~2)~OINDEX(U.J,I)+VIMDA”TIM2-OlNDEX(US,J.I))-

Fig. 2. 123 macro.

1 2 3 4 5 6 7 a 9 10 11 12 13 14 15 16 17

IL

117

A

Fig. 3. Initial temperature for the numerical example.

IN

,..

10

118

A. Kharab

lc

A 25 26 27 28 29 30

31 32 33 34 35 36 37 36 39 40 41 42 43 44 45 46 47 40 49 50 51 52 53 54 55 56 57 56 59 60

I

t

,I[FOR I P3.m 8

m

I/ I I

I

IF

F

N,l,~l,IE37) I I

I

I I

I

I

I

I

I I(IF N=lllLB

I

I

I I I

I

4

I

Y.J’DX)-

r--llFORl.l.M-,‘l,1,IG53l IOI ,

I

Fig. 4. The thi:eeFOR loopsusedto computeu for n =

I

IU(x,$) at tc !(n-1)k I I j .:,: 1. ;, I: f

I f :!

I

The values of u at the actual time stepsare written in the range ‘UW’. In order to understand how the values of u are computed at each time step, we now follow the 123 macro from the beginning of loop 3 residing in cell ID25 At this step the initial value of u had already been calculated and written in range B7..L17 as shown in Fig. 3. The FOR loop 3 setsthe counter for the time step and branches the macro to cell IE37. The command in cell IE 37 writes the actual time two rows above the range UW. The next command in cell IE38 fills all cells in range UW with zeros and at the same time writes the boundary conditions. The inner FOR loop 4 sets the counter for the y-space variable and branches the macro to cell IF48. The command in cell IF48 computes the values of y to be used for computing the value of g in eqn (9). The command in cell IF49 just moves the cursor to the top-cell in the range UW and J columns to the right. The last inner FOR loop 5 sets the counter for the x-space variable and branches the macro to cell IG53 where the macro computes the value of u(x,y, t) using eqns (9) for n = 1 and (5) for n > 1. Figure 5 shows the computational molecule of eqn (5) in the spy.eadsheet.

n

I I

I I I

1 I I

.o.i:‘:-:;l:.:“c~.:.. !:

1,2, . . . . I

:,

.:...~.:.:;:.,-~:.;‘Giii I

-iU(x,$) at tn /nk i ! j &:.::op j @/.‘i:;.l:. ;. j., i o,.::‘i ,. i. i-n.ic

I

I

I

I:‘.::, .,,,aF ..,.bj.;., d.j,

/I

I:

!

@ji : o,’ :, o i

I

I

I

i j

i /

I

! ,

j

;

: /

1

1

i j j

:j /

Ran e ‘U’

E

U(x,i) at t= j(n+l)k : j 1 O.‘...O! .---i 0; ./, i I)/:. ,‘,I

I

/

i i

I

.I

; .; ::{..: i;.

/U(lh,Jk(n+l)k) = 2(1-2LAM~A~2~0INDU((U,J,I]*LAMDA~2~(OIN~(U,J+1,I)~ I ; +OIN~EX(U,J,I+1)*OIND~(U,+1,I)cOlNDU((~,J,Cl)~OlNDU((US~J.I)

Fig. 5.

Computationalmoleculeof eqn (5).

; :

j / I

]

Finite-difference solution to a 2-D vibrating membrane problem 1

A

3 4 5 6 7

A M

I

TMAX

0

c

B I 101

I 1.01

I

07

I 0.05 I

I

MA

F I 11

I

a

n

l= I 11

P(w)

I

f(x.vl

w 1

I I

o~oslN(oPl~oslN(oplr, I 1

119 .I

K

L

I

I

I

I

I

I

8 IQ 11 12 13 14 15 16 17 10 19 20 21 22 23 24 25 26 27 26 29 30 31 32 33 34 35 36 37 36 39 40 41 42 43 44 45 46 47 46 49 50 51 52 53 54 55 56 57 56 59 80 61 62 83 64 66 66 67 : 70 71 72 73

Fig. 6. Output of the 123macro for our model problem at t = 0, 0.2, 0.4, 0.8 and 1.0. Before moving to the next time level (n + 2)k, the range name US is deleted, the range U is renamed US, the range UW is renamed U and a new range below the range U, where the new output of u is to be written for the next time step, is named UW. These commands are contained in cells IE42..IE45.

4 NUMERICAL RESULTS Example As an example we consider a unit square membrane (a=l)witha=l.Th e initial displacement is given by

A. Kharab

Fig. 7. 3D view of the initial displacement f(x, y).

Fig. 8. 3D view of u(x, y, t) at t = 1.0.

the function f(x,r) = sin rxsinr~~ and the initial velocity g(x, y) is taken to be zero. One of the main advantages of using spreadsheetsis that data can be easily entered. For our example enter the data as follows: Move the cursor to cell A2 and type 10 to enter M, then move the cursor to cell B2 and type 3 to enter TMAX and so on as shown in Fig. 5. To enter the function f(x, y), move the cursor to cell G2 and type @sin(@PI*x)*@sin(@PI*y). Figure 3 shows the output of the 123 macro for our model problem with M = 10, k = 0.05, TMAX = 3. Figure 6 shows a 3-D mesh of the displacement at t = 1.0. To run the 123 macro one needs to name it first by moving the cursor to the first cell where the macro starts and use the command /rnc. In our case since the 123 macro starts in cell IDI, we use the command /rnc\A”IDl” to name it \A. Invoke now the 123 macro by pressing Alt-A. The 123 macro was run using Lotus 123 version 3.4 for DOS and the total execution time was 1min and 10s using a GATEWAY 2000 486/33 Personal computer.

built-in graphics capabilities that permit graphs to be viewed on the screen or printed with little effort. Spreadsheets methods can be extended to solve more advanced multidimensional problems for student research projects in engineering and sciences.

5 CONCLUSION

A spreadsheet program has been used to find the numerical solution of a two-dimensional wave equation arising in the vibration of a square membrane. The sample problem used in the paper shows that spreadsheet are not only limited to do basic computations but through user ‘macros’, they can function in the same way as programs written in high-level language without the need of compiling the program. Data and functions are easily entered by simply moving the cursor to the appropriate cells. Moreover, spreadsheetsinclude

ACKNOWLEDGEMENT

The author wishes to acknowledge the support of King Fahd University of Petroleum & Minerals, Dhahran, Saudi Arabia.

REFERENCES 1. Hagler, M. Spreadsheet solution of partial differential equations. IEEE Trans. Educ., 1987, E30, 130-34. 2. Hart, F. X. Validating spreadsheetsolutions to Laplace’s equation. Am. J. Phys., 1987, 57, 1027-34. 3. Eid, J. C. A methodology and tutorial for thermal modeling with PC spreadsheets. Heat transfer Engng, 1987,8,95-107. 4. Kerkes, D. J. Spreadsheetsolution for transient fluid flow in a porous media with a moving drainage boundary. Computational Mechanics Publications. Proc. First Znt. Co&, Southampton U.K., Vol. 1, pp. 49-62, 1991. 5. Crow, T. T. Solutions to Laplace’s equation using spreadsheetson a personal computer. Am. J. Phys., 1987, 55, 817-23. 6. Hoffman, J. D. Numerical Methods for Engineers and Scientists. McGraw-Hill, Singapore, 1992. 7. Ames, W. F. Numerical Methods for Partial D@erential Equations (second edition), Academic Press, New York, 1977. 8. Peabody, F., Nyberg, D. W. & Dtmford W. G. The use of spreadsheetprogram to design motors on 1987, a personal computer. IEEE Trans. Znd. Appl., 1987, I&23(3), 520-5. 9. Rao, N. D. Typical applications of microcomputer spreadsheets to electrical engineering problems. IEEE Trans. Educ., 1984, E-27, 237-42.