Dynamics of a circular membrane with an eccentric circular areal constraint: Analysis and accurate simulations

Simulation Modelling Practice and Theory 31 (2013) 149–168

Contents lists available at SciVerse ScienceDirect

Simulation Modelling Practice and Theory journal homepage: www.elsevier.com/locate/simpat

Dynamics of a circular membrane with an eccentric circular areal constraint: Analysis and accurate simulations Assaad Alsahlani, Ranjan Mukherjee ⇑ Department of Mechanical Engineering, 2555 Engineering Building, Michigan State University, East Lansing, MI 48824-1226, United States

a r t i c l e

i n f o

Article history: Received 4 June 2012 Received in revised form 10 October 2012 Accepted 12 October 2012 Available online 23 December 2012 Keywords: Circular membrane Constrained membrane Orthogonality Vibration Accuracy Simulations Animations

a b s t r a c t The dynamics of a circular membrane with an eccentric circular areal constraint is studied under arbitrary initial conditions. The symmetric and antisymmetric modes of vibrations are investigated and the results are compared with those in the literature. The accuracy of the eigenfrequencies and mode shapes are studied for different sizes and locations of the constraint and it is shown that proper choice of the number of angular modes is critical to accurate computation of the mode shapes; fewer or more number of angular modes can provide relatively accurate eigenfrequencies but completely inaccurate mode shapes. The orthogonality of distinct modes of constrained membranes is mathematically established in this work. The orthogonality property of the modes is used to compute the modal coefficients for simulation of the dynamics. In comparison to prior work, where the objective has been limited to computing the first few eigenfrequencies, this work presents a method for accurately computing the eigenfrequencies, mode shapes, and modal coefficients needed for dynamics simulation. Numerical simulations and video animations of vibration of constrained membranes are presented for arbitrary initial conditions. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction The dynamics of a membrane is a classical problem in mechanical vibrations. Among the different types of membranes, circular membranes are the most commonly studied due to their large number of applications. From the study of musical notes of percussion instruments, circular membranes have been used to design diaphragms for condenser microphones, model the dynamics of the human ear [1], understand the vibration characteristics of membrane mirrors and gossamer structures [2], measure surface tension [3], and design ink-jet printers [4]. For a majority of these applications, unconstrained membrane models have been used. A constrained membrane is one in which a portion of the membrane does not undergo vibration, and an application where constrained membrane models have been used is waveguides. Early work on this subject dates back to the 1960s and 1970s. Yee and Audeh [5,6] used the point-matching method to determine eigenfrequencies of waveguides with eccentric cross-section. Davies and Muilwyk [7] and Steele [8] used the finite difference method to compute eigenvalues of waveguides with arbitrary cross-section. The Helmholtz equation, which describes the dynamics of waveguides, was solved by Arlett et al. [9] and Gass [10] using finite element methods, by Laura [11] and Hine [12] using the Galerkin method, and by Bulley and Davies [13] using the Rayleigh–Ritz method. The similarity between the differential equations of membranes and waveguides motivated the study of circular membranes with constraints in the 1970s and 1980s. The eigenvalues of a circular membrane with an internal eccentric circular areal constraint were first computed by Nagaya [14] using Fourier series. Another solution to this problem was presented by Lin [15], who used Graf’s addition theorem to transform the Bessel functions from one set of polar coordinates to another. Lin ⇑ Corresponding author. Tel.: +1 517 355 1834. E-mail address: [email protected] (R. Mukherjee). 1569-190X/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.simpat.2012.10.008

150

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

[15] missed a term in the final expression, which affected the results of calculation. The missing term was identified by Singh and Kothari [16] and Singh et al. [17], who used Nagaya’s approach together with Graf’s addition theorem. Later, Nagaya and Hai [18] studied composite membranes with arbitrarily shaped inner and outer boundaries. The vibration of membranes remains an active subject of research with recent work focusing on complex boundary conditions and shapes, inhomogeneities, accuracy of solutions, and applications. Investigations of membranes with complex boundary conditions and shapes include the work by Laura and Gutierrez [19], Wu and Wang [20], Kang and Lee [21], and Chen et al. [22]. Investigations of non-homogeneous membranes include the work by Cortinez and Laura [23] and Cap [24], and that of accuracy of solutions include the work by Zhao and Stevens [25]. Applications of membranes are diverse and membrane vibrations have been investigated in the context of micro-air vehicle wings [26–28], energy harvesting [29,30], and vibration suppression [31]. A thorough review of the literature indicates that investigations of constrained membranes have been limited to computing the eigenfrequencies of vibration. Although accurate computation of the eigenfrequencies may be sufficient for many applications, it is not sufficient for accurate computation of the mode shapes and/or modal coefficients of a constrained membrane. Simulation of the dynamics requires computation of the modal coefficients using orthogonality property of the modes. The orthogonality property of modes is well-known but it has not been established mathematically for constrained membranes. We establish the orthogonality property of all modes of a circular membrane with an internal circular areal constraint for the first time in this paper. The computation of the modal coefficients require accurate computation of the mode shapes and this requires proper choice of the number of angular modes. A small number of modes result in mismatch of boundary conditions whereas a large number results in the eigenmatrix having columns of very small and very large numbers due to high sensitivity of the Bessel functions. We present an algorithm that determines the appropriate number of angular modes for arbitrary size and location of the constraint. In addition to this algorithm, which computes the mode shapes, we provide an algorithm for accurate computation of the modal coefficients for arbitrary initial conditions. Although simulations of unconstrained membranes have been presented in the literature [32–35], our paper is the first effort that provides a method for accurate simulation of constrained membrane dynamics. Our paper is organized as follows. A formal problem statement and a list of assumptions are provided in Section 2. The expression for the symmetric and antisymmetric modes are derived in Section 3 using Lin’s analysis [15]. In Section 4 we establish orthogonality between the symmetric modes, the antisymmetric modes, and between the symmetric and antisymmetric modes of a constrained membrane. In Section 5 we provide the algorithm for accurate computation of the mode shapes, and in Section 6 we provide the algorithm for computation of the modal coefficients. Simulations and video animations of the dynamics of a membrane with different initial conditions and constraints are also presented in Section 6. Section 7 provides concluding remarks. 2. Problem statement and assumptions Consider the circular membrane of radius R and outer boundary Co, shown in Fig. 1. A circular area of radius a and boundary Ci , located at a distance d from the center of the membrane, has zero displacement at all times; this constrains the dynamics of the membrane. To study the vibration of the active region of the membrane, we make the following assumptions: 1. The membrane is homogeneous and has a constant mass per unit area equal to l. The tension in the membrane is equal to T and remains constant at all times. The membrane undergoes transverse vibration and it is not affected by gravity. 2. The amplitude of oscillation of the membrane is small and its motion can be expressed by the standard relation in polar coordinates [36]:

Fig. 1. A circular membrane with an inner circular areal constraint.

151

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

@ 2 Z 1 @Z 1 @ 2 Z 1 @2Z þ þ 2 2¼ 2 2; 2 @r r @r r @h c @t

c¼

pffiffiffiffiffiffiffiffiffi T=l

ð1Þ

where Z = Z(r, h, t) is the transverse displacement. The distance of any point P is measured from the center of the inner circle O and is denoted by r; h denotes the angle that vector r makes with the x-axis. 3. The transverse displacement of all points in the area enclosed by the inner circle is zero, i.e., Z(r, h, t) = 0 for r 2 [0, a], h 2 [0, 2p]. The displacement of the membrane is also zero along the outer boundary Co. 3. Background – modes of vibration 3.1. General solution A general solution to Eq. (1) is obtained by separation of variables. Assuming Z(r, h, t) = u(r, h) F(t) and rearranging terms, we get

( ) @2u 1 @u 1 @2u 1 1 @2F ¼ 2 þ þ 2 2 2 r @r r @h c F @t 2 u @r 1

ð2Þ

Let k2 be the separation constant. Then, the right-hand side of Eq. (2) gives

1 1 @2F 2 ¼ k ; c2 F @t 2

or

@2F 2 þ c2 k F ¼ 0 @t 2

If we define the natural frequency x , c k,1 we get

FðtÞ ¼ A cos xt þ B sin xt

ð3Þ 2

where A and B are constants whose values depend on initial conditions. Equating the left-hand side of Eq. (2) to k and rearranging terms, we get

r2

@2u @u @2u 2 þr þ k r2 u ¼  2 2 @r @r @h

ð4Þ

To separate the radial and angular terms in Eq. (4), we substitute u(r, h) = U(r) V(h); this yields

( ) 1 2 @2U @U 1 @2V 2 2 ¼ r þ k þ r r 2 U @r @r V @h2

ð5Þ

By choosing m2 to be the separation constant, we get

1 @2V ¼ m2 ; V @h2

)

VðhÞ ¼ C m cos mh þ Dm sin mh

ð6Þ

where Cm and Dm are constants whose values depend on initial conditions. Since we have

VðpÞ ¼ VðpÞ;

@V @V ðpÞ ¼ ðpÞ @h @h

it can be shown2 that m takes integer values 0, 1, 2, . . .. From the left-hand side of Eq. (5) we get

r2

@2U @U 2 þ ðk r 2  m2 ÞU ¼ 0 þr @r 2 @r

ð7Þ

For any integer value of m, the general solution to Eq. (7) is as follows:

UðrÞ ¼ bm J m ðkrÞ þ cm Y m ðkrÞ where bm and cm are constants, and Jm(kr) and Ym(kr) are Bessel functions [36] of the first and second kind, respectively. The general solution to the equation of motion of the circular membrane can now be written as

Zðr; h; tÞ ¼

1 X fbm J m ðkrÞ þ cm Y m ðkrÞgfC m cos mh þ Dm sin mhgfA cos xt þ B sin xtg

ð8Þ

m¼0

The solution in Eq. (8) can be rewritten as

1 The separation constant k is related to the natural frequency by the relation k = x/c and has the unit of rad/m. Since c is a constant, the value of k is proportional to the frequency. In the sequel we will refer to k as the eigenfrequency. 2 The integer m is the number of angular nodes. For an unconstrained membrane, it is equal to the number of diametrical lines along which the membrane has zero displacement.

152

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

Zðr; h; tÞ ¼ fus ðr; hÞ þ uc ðr; hÞgfA cos xt þ B sin xtg

ð9Þ

where

us ðr; hÞ ¼ uc ðr; hÞ ¼

1 X fb1m J m ðkrÞ þ c1m Y m ðkrÞg sin mh m¼1 1 X

ð10Þ

fb2m Jm ðkrÞ þ c2m Y m ðkrÞg cos mh

m¼0

and b1m = bmDm, c1m = cmDm, b2m = bmCm and c2m = cmCm. The constants in Eq. (9), namely, k, A, B, m, b1m, c1m, b2m, c2m, can be determined by applying boundary conditions and initial conditions. 3.2. Boundary conditions The boundary conditions are: 1. The membrane has zero displacement for r = a. From Eq. (9) we get

Zða; h; tÞ ¼ fus ða; hÞ þ uc ða; hÞgfA cos xt þ B sin xtg ¼ 0 This implies {us(a, h) + uc(a, h)} = 0, or 1 X

fb1m J m ðkaÞ þ c1m Y m ðkaÞg sin mh þ

m¼1

1 X fb2m J m ðkaÞ þ c2m Y m ðkaÞg cos mh ¼ 0

ð11Þ

m¼0

2. The membrane has zero displacement for all r values for which q = R. This implies 1 X

fb1m J m ðkrjq¼R Þ þ c1m Y m ðkrjq¼R Þg sin mh þ

m¼1

1 X fb2m J m ðkrjq¼R Þ þ c2m Y m ðkrjq¼R Þg cos mh ¼ 0

ð12Þ

m¼0

The solutions to Eqs. (11) and (12) were investigated by several researchers [14–17] and the eigenfrequencies were determined for the even (symmetric) and the odd (antisymmetric) modes.

3.3. Even or symmetric modes The symmetric modes are associated with the cosine term cos mh in Eq. (10) and are denoted by uc(r, h). The superscript ‘‘c’’ is used to represent the cosine terms. Since the symmetric modes satisfy the differential equation independently, the boundary conditions in Eqs. (11) and (12) give 1 X fb2m Jm ðkaÞ þ c2m Y m ðkaÞg cos mh ¼ 0 m¼0 1 X

ð13Þ

fb2m J m ðkrjq¼R Þ þ c2m Y m ðkrjq¼R Þg cos mh ¼ 0

m¼0

Eq. (13) cannot be solved directly for k since the second equation contains rjq=R, which is a function of h. From Fig. 1, we can show

rjq¼R ¼ rðhÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 R2 þ d  2Rd cos w ¼ R2  d sin h  d cos h

ð14Þ

To solve Eq. (13) analytically, Lin [15] used Graf’s addition formula [37] to transform the Bessel functions from (r, h) coordinates to (q, w) coordinates. This transformation is given below

J m ðkrÞeimh ¼

1 X

Jmþq ðkqÞJ q ðkdÞeiðmþqÞw

q¼1

Y m ðkrÞeimh ¼

1 X

ð15Þ

Y mþq ðkqÞJ q ðkdÞeiðmþqÞw

q¼1

Using the expressions in Eq. (15), Eq. (13) is rewritten as follows 1 X fb2m Jm ðkaÞ þ c2m Y m ðkaÞg cos mh ¼ 0 m¼0 1 X

1 X

m¼0 q¼1

ð16Þ fb2m J mþq ðkRÞ þ c2m Y mþq ðkRÞgJ q ðkdÞ cosðm þ qÞw ¼ 0

153

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

Eq. (16) derived by Lin [15] was revised by Singh et al. [17]; they expanded the inner summation and used the property of Bessel functions

J m ðxÞ ¼ ð1Þm Jm ðxÞ;

Y m ðxÞ ¼ ð1Þm Y m ðxÞ

to obtain 1 X fb2n J n ðkaÞ þ c2n Y n ðkaÞg cos nh ¼ 0 n¼0 1 X 1 X

ð17Þ J n ðkRÞ þ c2m Y n ðkRÞgfJ nm ðkdÞ þ ð1Þm J nþm ðkdÞg cos nw ¼ 0

n fb2m

m¼0 n¼0

where

n is defined as

n ¼



1=2 if n ¼ 0 1

if n > 0

Expanding Eq. (17) from m = 0 to N horizontally and from n = 0 to N vertically, gives the matrix equation

Pc X c ¼ 0 c

ð18Þ

(2N+2)(2N+2)

where P 2 R

2

J 0 ðf1 Þ

6 1 J ðf ÞHc 6 2 0 2 00

6 6 0 6 6 J ðf ÞHc c 6 2 1 10 P ¼6 6 .. 6 . 6 6 4 0

J N ðf2 ÞHcN0

X c ¼ ½ b20

(2N+2)

and A 2 R

are given by the relations

Y 0 ðf1 Þ

0

0



0

0

1 Y ðf ÞHc00 2 0 2

1 J ðf ÞHc01 2 0 2

1 Y ðf ÞHc01 2 0 2



1 J ðf ÞHc0N 2 0 2

1 Y ðf ÞHc0N 2 0 2

0

J 1 ðf1 Þ

Y 1 ðf1 Þ



0

0

Y 1 ðf2 ÞHc10 .. .

J 1 ðf2 ÞHc11 .. .

Y 1 ðf2 ÞHc11 .. .

 .. .

J N ðf2 ÞHc1N .. .

Y N ðf2 ÞHc1N .. .

0 Y N ðf2 ÞHcN0

0 J N ðf2 ÞHcN1

0 Y N ðf2 ÞHcN1

 

J N ðf1 Þ J N ðf2 ÞHcNN

Y N ðf2 ÞHcNN

Y N ðf1 Þ

3 7 7 7 7 7 7 7 7 7 7 7 7 5

c20 b21 c21    b2N c2N T c

and where f1 = ka, f2 = kR, and Hcnm ¼ fJ nm ðkdÞ þ ð1Þm J nþm ðkdÞg. Eq. (18) is used to solve for the eigenfrequencies ki and eigenvectors X ci for i = 1, 2, . . .. Using Eq. (10), the symmetric modes can be subsequently obtained as follows:

uci ðr; hÞ ¼ b20;i

N  X b2m;i m¼0

b20;i

c

J m ðki rÞ þ

c2m;i b20;i

 c Y m ðki rÞ cos mh

ð19Þ

where the terms (b2m, i/b20,i), (c2m, i/b20,i), m = 0, 1, . . . , N are entries of the ith eignevector X ci , and b20,i is a scaling factor for the ith symmetric mode. 3.4. Odd or antisymmetric modes The antisymmetric modes are associated with the sine term sin mh in Eq. (10) and are denoted by us(r, h); the superscript ‘‘s’’ is used to represent the sine terms. For the symmetric modes, the boundary conditions in Eqs. (11) and (12) give 1 X fb1m J m ðkaÞ þ c1m Y m ðkaÞg sin mh ¼ 0 m¼1 1 X

ð20Þ

fb1m J m ðkrjq¼R Þ þ c1m Y m ðkrjq¼R Þg sin mh ¼ 0

m¼1

Using Graf’s addition formula [37], the boundary conditions in Eq. (20) can be rewritten as 1 X fb1n J n ðkaÞ þ c1n Y n ðkaÞg sin nh ¼ 0 n¼1 1 X 1 X fb1m J n ðkRÞ þ c1m Y n ðkRÞgfJ nm ðkdÞ  ð1Þm J nþm ðkdÞg sin nw ¼ 0

ð21Þ

m¼1 n¼1

Expanding Eq. (21) from m = 1 to M and n = 1 to M, M = (N + 1), gives the matrix equation

Ps X s ¼ 0

ð22Þ

154

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

where Ps 2 R2M2M and B 2 R2M are given by the relations

2

J 1 ðf1 Þ

6 J ðf ÞHs 6 1 2 11 6 6 0 6 6 J ðf ÞHs s 6 2 21 P ¼6 2 6 .. 6 . 6 6 4 0

J M ðf2 ÞHsM1

Y 1 ðf1 Þ

0

0



0

Y 1 ðf2 ÞHs11

J1 ðf2 ÞHs12

Y 1 ðf2 ÞHs12



J1 ðf2 ÞHs1M

0

J 2 ðf1 Þ

Y 2 ðf1 Þ



0

Y 2 ðf2 ÞHs21 .. .

J2 ðf2 ÞHs22 .. .

Y 2 ðf2 ÞHs22 .. .

 .. .

J2 ðf2 ÞHs2M .. .



J M ðf1 Þ

0

0

0

Y M ðf2 ÞHsM1

J M ðf2 ÞHsM2

Y M ðf2 ÞHsM2

   J M ðf2 ÞHsMM

0

3

Y 1 ðf2 ÞHs1M 7 7 7 7 0 7 Y 2 ðf2 ÞHs2M 7 7 7 7 .. 7 . 7 7 Y M ðf1 Þ 5

Y M ðf2 ÞHsMM

c11 b12 c12    b1M c1M T

X s ¼ ½ b11

s

and where f1 = ka, f2 = kR, and Hsnm ¼ fJ nm ðkdÞ  ð1Þm J nþm ðkdÞg. Eq. (22) is used to solve for the eigenfrequencies ki and eigenvectors X si for i = 1, 2, . . .. Using Eq. (10), the symmetric modes can be subsequently obtained as follows:

usi ðr; hÞ ¼ b11;i

M  X b1m;i m¼1

b11;i

  s  c1m;i  s  J m ki r þ Y m ki r sin mh b11;i

ð23Þ

   where the terms b1m;i =b11;i ; c1m;i =b11;i ; m ¼ 1; 2; . . . ; M are entries of the ith eignevector X si , and b11;i is a scaling factor for the ith antisymmetric mode. From the expressions in Eqs. (19) and (23) it can be seen that each mode shape of a constrained membrane ui is a summation of terms that are a function of m, i.e., each mode shape depends on all the angular nodes. This is different from unconstrained membranes where each mode shape is associated with a specific angular node only and is denoted by ui,m, i.e., a pair of subscripts are used to identify the mode shape. Since each mode shape is associated with an eigenfrequency, the eigenfrequencies of unconstrained membranes are denoted by ki,m whereas the eigenfrequencies of constrained membranes are denoted by ki. In the literature, the subscripts for frequencies and mode shapes of unconstrained membranes are numbered starting with zero, i.e., i, m = 0, 1, 2, . . .. In a slight deviation from this tradition, the frequencies and mode shapes of constrained membranes have been numbered starting with one in this paper, i.e., i = 1,2,3, . . .. 4. Orthogonality of modes The differential equation in Eq. (2) is related to the eigenfrequency k as follows

( ) @2u 1 @u 1 @2u 1 1 @2F 2 ¼ 2 þ ¼ k þ 2 2 2 r @r r @h c F @t2 u @r 1

and can be rewritten in the form

  1 @ @u 1 @2u 2 r þk u¼0 þ 2 r @r @r r @h2

ð24Þ

Let ki and kj be the ith and jth eigenfrequency associated with the modes ui and uj, respectively. Each of the modes ui and uj may be symmetric or anti-symmetric. For these two modes, Eq. (24) can be written as follows

  1 @ @ ui 1 r þ 2 r @r r @r   @ uj 1 @ 1 r þ 2 r @r r @r

@ 2 ui @h2 @ 2 uj @h2

2

ð25Þ

2

ð26Þ

þ ki ui ¼ 0 þ kj uj ¼ 0

Multiplying Eq. (25) by uj and Eq. (26) by ui and subtracting one from the other, we get

( )     

n o @ uj @ 2 uj 1 @ @ ui @ 1 @ 2 ui r r uj  ui þ 2 uj  2 ui þ k2i  k2j ui uj ¼ 0 2 r @r @r r @r @r @h @h

ð27Þ

Integrating Eq. (27) over the active region of membrane now gives

Z 2p Z 0

)    

Z 2p Z rjq¼R ( 2 n o @ uj @ 2 uj @ @ ui @ 1 @ ui uj  ui dr dh þ uj  2 ui dr dh þ k2i  k2j r r 2 @r @r r @h @r @r @h 0 a



rjq¼R

a



Z 2p Z 0

a

rjq¼R

ui uj r dr dh ¼ 0

ð28Þ

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

155

Integration of the first term by parts gives

Z 2p Z 0

rjq¼R

a



rjq¼R    

Z 2p  @ uj @u @ @ ui @ @ ui r r uj  ui dr dh ¼ r uj  r j ui dh @r @r @r @r @r @r 0 r¼a  Z 2p Z rjq¼R  @ ui @ uj @ uj @ ui r dr dh   @r @r @r @r 0 a

Since both ui and uj are zero for r = a and r values where q = R, the first term on the right-hand side of the above equation is zero. The second term on the right-hand side of the above equation is zero trivially and this implies

Z 2p Z 0

rjq¼R

a

( ) n o Z 2p Z rjq¼R @ 2 uj 1 @ 2 ui 2 2 u  u ui uj r dr dh ¼ 0 j i dr dh þ ki  kj 2 2 r @h @h 0 a

ð29Þ

The first term in Eq. (29) is given by the expression

Z 2p Z 0

rjq¼R

a

( ) @ 2 uj 1 @ 2 ui u  u i dr dh r @h2 j @h2

ð30Þ

where the limits of integration correspond to the unconstrained area of the membrane in Fig. 1. In Eq. (30), the integration is first carried out with respect to r and then with respect to h. For the first integration, the upper limit of r is expressed as a function of h; this functional dependence of r on h is given by Eq. (14) and shown in Fig. 2a. We evaluate the expression in Eq. (30) by changing the order of integration. This requires h to be expressed as a function of r, which is obtained from the geometry in Fig. 1:

" 1

h ¼ hðrÞ ¼ cos

# 2 R2  r2  d ; 2rd

h 2 ½0; p

ð31Þ

The area of integration is redrawn in Fig. 2b and the integral in Eq. (30) is expressed as the sum of the following three terms:

) ) Z Rþd Z hðrÞ ( 2 Z þp ( 2 @ 2 uj @ 2 uj 1 @ ui 1 @ ui uj  2 ui dh dr þ u  2 ui dh dr r @h2 j @h2 @h @h a Rd p p r ( ) Z Rþd Z þp 2 2 @ u 1 @ ui þ uj  2 j ui dh dr 2 r @h @h Rd þhðrÞ

Z

Rd

ð32Þ

where h(r) is defined by Eq. (31). For each of the three terms in Eq. (32), the inner integration with respect to h by parts, gives

Z " 2 @ ui h

2

@h

uj 

@ 2 uj 2

@h

#

ui dh ¼

 hu Z   hu  @u @u @ ui @ ui @ uj @ uj @ ui @ ui uj  j ui  uj  j ui  dh ¼ @h @h @h @h @h @h @h @h h hl hl

ð33Þ

where hl and hu denote the lower and upper limits of h. Using Eq. (33), the three terms in Eq. (32) can be simplified to the form

Z a

Rd

 þp hðrÞ þp Z Rþd  Z Rþd  @u @u @u 1 @ ui 1 @ ui 1 @ ui uj  j ui dr þ uj  j ui dr þ uj  j ui dr r @h @h @h @h @h @h Rd r Rd r p p þhðrÞ

(a)

ð34Þ

(b)

Fig. 2. Shaded region showing the area of integration for the two cases where the integration is done (a) first with respect to r and then with respect to h and (b) first with respect to h and then with respect to r.

156

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

Since h = +h(r) and h = h(r) describe the outer boundary of the membrane where ui = uj = 0 for all i and j, the above integrals simplify to the form

Z

Rd

a

 þp þp þp Z Rþd  Z Rþd  @u @u @u 1 @ ui 1 @ ui 1 @ ui uj  j ui dr þ uj  j ui dr ¼ uj  j ui dr r @h r @h @h @h @h @h Rd r a p p p

ð35Þ

Now consider the three separate cases: 1. Modes ui and uj are both symmetric: Using Eq. (19), both (@ ui/@ h) and (@ uj/@h) can be expressed in the form 1 @ ui X ðÞ sin mh; ¼ @h m¼0

1 @ uj X ðÞ sin mh ¼ @h m¼0

Since their values are zero at h = ±p, the integral in Eq. (35) is identically zero. 2. Modes ui and uj are both anti-symmetric: It can be see from Eq. (23) that ui and uj have the form

ui ¼

1 X ðÞ sin mh; m¼0

uj ¼

1 X ðÞ sin mh m¼0

Since their values are zero at h = ±p, the integral in Eq. (35) is identically zero. 3. Mode ui is symmetric and mode uj is anti-symmetric: Using Eqs. (19) and (23) it can be shown that both the terms (@ ui/ @h)uj and (@ uj/@h)ui are zero at h = ±p, and therefore the integral in Eq. (35) is identically zero. This proves that the integral in Eq. (30) is zero for any two distinct modes ui and uj, irrespective of whether they are both symmetric, both antisymmetric, or one of them is symmetric and the other is antisymmetric. From Eq. (29) we can now conclude

Z 2p Z 0

a

rjq¼R

ui uj r dr dh ¼ 0; i – j

ð36Þ

which is the orthogonality condition for a pair of distinct modes. 5. Accurate computation of eigenfrequencies and mode shapes 5.1. Computational algorithm We present an algorithm for numerical computation of the eigenfrequencies and their corresponding mode shapes. The algorithm is identical for symmetric and antisymmetric modes and therefore we discuss it here for the symmetric modes c only. The algorithm first computes all eigenfrequencies ki ; i ¼ 1; 2; . . ., that are less than a user-specified maximum eigen c c frequency kmax . The interval 0; kmax is divided into small segments of length Dk and the eigenfrequencies within each segment are computed separately. The number of eigenfrequencies that show up within each segment depend on the value of N, which is related to the size of the matrix Pc in Eq. (18). Starting from zero, increasing the value of N typically results in larger number of zero crossings of the determinant Pc (see Fig. 3), and hence a larger number of eigenfrequencies within a specific segment. Beyond a certain value of N, the number of eigenfrequencies in a given segment will not change. At this point, increasing the value of N produces very small changes in the values of the eigenfrequencies but such changes can result in significant improvement or loss of accuracy of the associated mode shapes. The flowchart shown in Fig. 4 provides the

Fig. 3. The eigenfrequencies are the zero crossings of the determinant of the matrix Pc for symmetric modes, and Ps for antisymmetric modes.

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

157

 c algorithm for computing all eigenfrequencies in the interval 0; kmax , and accurately computing the mode shape associated with each eigenfrequency  through proper choice of N. A brief summary of the algorithm is provided next. c The interval 0; kmax is divided into small segments of length Dk. For any given segment, the boundaries of the segment are denoted by d0 and df. The value of k is changed in increments of d (small number) from d0 to df, and zero crossings of j Pc j are stored as potential eigenfrequencies. For a given segment, this procedure is repeated by increasing N, N 2 {0, 1, 2, . . . , Nmax}, until no additional zero crossings occur and the mode shapes associated with the eigenfrequencies have been computed accurately. The mode shapes are computed using Eq. (19) and are checked for accuracy by verifying uci ¼ 0 along the inner and outer boundary of the membrane. The procedure is repeated for each segment until the maximum userspecified value of k = kmax is reached. In the literature, several methods have been proposed for computing the eigenfrequency and mode shapes of constrained membranes. Most of these methods are well suited to computing the first few eigenfrequencies and mode shapes only. For simulation of the dynamics of a constrained membrane, it is necessary to accurately compute many frequencies and mode shapes; this is difficult to compute using existing methods due to singularity. For example, Nagaya [14], Singh and Kothari [16], and Singh et al. [17] proposed to combine the two equations for symmetric modes in Eq. (13) into the following single equation



1 X J ðkaÞ b2m J m ðkrjq¼R Þ  m Y m ðkrjq¼R Þ cos mh ¼ 0 Y m ðkaÞ m¼0

ð37Þ

This procedure reduces the dimension of the coefficient matrix by a factor of two and reduces the computational burden but the presence of the Bessel function of the second kind, Ym(ka), in the denominator makes it prone to singularity. This problem becomes critical during the computation of higher frequencies.

Fig. 4. Algorithm for computing eigenfrequencies and mode shapes for symmetric modes.

158

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

The eigenfrequencies and mode shapes of the antisymmetric modes can be computed using the algorithm in Fig. 4 by c s replacing Pc ; ki , and uci with P s ; ki , and usi , respectively, and N with M. Ps is defined in Eq. (22) and usi is defined in Eq. (23). 5.2. Numerical results We use the algorithm in Section 5.1 to accurately compute the eigenfrequencies and mode shapes of a constrained membrane with R/a = 0.25 and d/a = 1.0. This problem was considered in the literature [6,14,16] and the results were summarized by Singh et al. [17]. We compare these results with those obtained using our algorithm, implemented with kmax = 10a and Dk = kmax. Table 1 shows the non-dimensional fundamental antisymmetric frequencies for the constrained membrane obtained using our algorithm; the value of M is gradually increased until all modes are accurately computed. As mentioned in Section 5.1, the accuracy of each mode shape is verified by computing the displacement of 100 points along the inner and outer membrane boundaries. For accuracy, the displacement of the membrane at all these points must satisfy  6 104. In Table 1, eigenfrequencies associated with accurate mode shapes (this depends on correct choice of M) have a star superscript. The eigenfrequencies obtained from the literature are presented in separate rows at the bottom of the table and are placed within parenthesis. The following observations can be made from the results presented in Table 1: 1. It is not possible to identify all the eigenfrequencies in an interval using a small value of M. Increasing the value of M helps in identifying more eigenfrequencies and M P 6 identifies all eight eigenfrequencies in this particular example. 2. For M = 6, all eigenfrequencies are identified but only the first eigenfrequency has an accurate mode shape. Increasing the value of M results in a larger number of accurate mode shapes and M = 12 results in accurate computation of all mode shapes. 3. A large value of M is not guaranteed to result in accurate mode shapes. For example, M = 13 gives all accurate mode shapes but none of the mode shapes are accurate for M = 40. This is because of the sensitivity of the Bessel functions. This justifies the need for increasing M gradually rather than choosing a pre-specified large integer. 4. Although all the mode shapes can be accurately computed using M = 12, the computational burden can be significantly reduced by using the lowest value of M for each eigenfrequency that gives an accurate mode shape. For example, the c c mode shape for k1 a should be computed using M = 5, the mode shape for k2 a should be computed using M = 7, and so on. 5. The eigenfrequencies in the literature have values that are very close to those obtained by our algorithm. However, these small errors result in large errors in the shape of the modes, which will violate boundary and orthogonality conditions. This will be discussed further later. We next present results for different sizes and locations of the constraint in the membrane. The specific values of a/R and d/ a used are taken from Singh et al. [17]. Table 2 presents the first five non-dimensional frequencies of symmetric modes of the constrained membranes and values of N used to compute the corresponding mode shapes. For a given set of a/R and d/a values, the first row contains eigenfrequencies and N values that correspond to accurate mode shapes. The procedure for obtaining accurate eigenfrequencies and mode shapes was discussed in Section 5.1. The second row contains eigenfrequencies that are closest to values presented in the literature; these eigenfrequencies are obtained by gradually increasing N until the error between the computed value and the value in the literature is minimum. In all cases, the error is minimum for a value of N for which the mode shape is inaccurate – this is indicated by the underlined eigenfrequencies. The eigenfrequencies in parenthesis are taken from the literature [17] and are presented in the third and fourth rows. The following observations can be made from the results presented in Table 2 pertaining to accurate eigenfrequencies and mode shapes: Table 1 First eight non-dimensional frequencies of antisymmetric modes of constrained membrane with a/R = 0.25 and d/a = 1.0. The results in parenthesis are from [6,14,17], and the star superscript indicates frequencies with accurate mode shapes computed using our algorithm. c

c

c

c

c

c

c

c

M

k1 a

k2 a

k3 a

k4 a

k5 a

k6 a

k7 a

k8 a

1 2 3 4 5 6 7 8 9 10 11 12 13 40 – –

1.1119 1.0636 1.0661 1.0660 1.0660⁄ 1.0660⁄ 1.0660⁄ 1.0660⁄ 1.0660⁄ 1.0660⁄ 1.0660⁄ 1.0660⁄ 1.0660⁄ 1.0647 (1.0659) (1.0500)

2.1342 1.4307 1.3811 1.3850 1.3848 1.3848 1.3848⁄ 1.3848⁄ 1.3848⁄ 1.3848⁄ 1.3848⁄ 1.3848⁄ 1.3848⁄ 1.3822 (1.3848) (1.3800)

– 1.9199 1.6875 1.6557 1.6591 1.6589 1.6589 1.6589⁄ 1.6589⁄ 1.6589⁄ 1.6589⁄ 1.6589⁄ 1.6589⁄ 1.6572 (1.6589) –

– – 1.9290 1.9288 1.9275 1.9284 1.9284 1.9284 1.9284⁄ 1.9284⁄ 1.9284⁄ 1.9284⁄ 1.9284⁄ 1.9272 (1.9283) –

– – 2.3177 1.9473 1.9296 1.9312 1.9311 1.9311 1.9311⁄ 1.9311⁄ 1.9311⁄ 1.9311⁄ 1.9311⁄ 1.9297 – –

– – – 2.3280 2.2195 2.2083 2.2101 2.2099 2.2099 2.2099 2.2099⁄ 2.2099⁄ 2.2099⁄ 2.2072 – –

– – – – 2.3292 2.3288 2.3288 2.3288 2.3288 2.3288⁄ 2.3288⁄ 2.3288⁄ 2.3288⁄ 2.3272 – –

– – – – – 2.4982 2.4916 2.4928 2.4927 2.4927 2.4927 2.4927⁄ 2.4927⁄ 2.4922 – –

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

159

Table 2 First five non-dimensional frequencies of symmetric modes of constrained membranes and values of N used to compute the mode shapes. For a given set of a/R and d/a values, the first row contains values that correspond to accurate mode shapes. The second row contains values that are closest to those presented in the literature – the eigenfrequencies are underlined to indicate that the corresponding mode shapes are inaccurate. The results from the literature [6,14,17] are presented in the following rows with eigenfrequencies in parenthesis.

1. Lower values of N are required to compute lower eigenfrequencies, and vice versa. This can also be observed from the data presented in Table 1. 2. Lower values of N are required for lower values of d, i.e., constraints located farther from the center of the membrane, and vice versa. 3. Higher values of N are required for smaller values of a, i.e., constraints of smaller size, and vice versa. The effect of N on the accuracy of eigenfrequencies was first observed by Nagaya [14], who used 8 6 N 6 12 to confine the error to 0.5%. Such choice of N may indeed result in accurate computation of the eigenfrequency but there is no guarantee that the mode shape will be accurate. The accuracy of the mode shape depends on the size and location of the constraint, as discussed above, and incorrect choice of N can result in significant error in the mode shape even when the eigenfrequency is c computed accurately. This is illustrated with the help of Fig. 5, which shows the accurate and inaccurate mode shapes for k3 a for a/R = 0.5 and d/a = 0.6, 0.8, and 1.0. It should be noted that the eigenfrequencies corresponding to these accurate and inaccurate mode shapes, shown in Table 2, indicate negligible error. For example, for a/R = 0.5 and d/a = 0.8, the error in the eigenfrequency is only

3:5975  3:5620  100  1% 3:5975 but the mode shapes in Fig. 5b are significantly different; and of course, the inaccurate mode shape does not satisfy the boundary conditions. Similar observations can be made for the cases a/R = 0.5, d/a = 0.6 and a/R = 0.5, d/a = 1.0. For these cases, the error in the eigenfrequency is   0.06% and 0.02%, respectively but the corresponding mode shapes are significantly different. We complete this section with numerical results for an example where the size of the constraint is small and the location of the constraint is far away from the center of the membrane, i.e., a/R is small and d/a is large. For such cases, it is difficult to obtain the eigenfrequencies and mode shapes accurately because the Bessel functions are sensitive to large values of d and small values of a. Despite this sensitivity, our algorithm was able to compute the first twenty modes of a membrane with a/R = 0.1 and d/a = 9.0. The first symmetric mode and the sixteenth symmetric mode of the membrane are shown in Fig. 6

160

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

(a)

(b)

(c)

c

Fig. 5. Modes shapes k3 a for constrained membranes in Table 2 with a/R = 0.5 and (a) d/a = 0.6, (b) d/a = 0.8, and (c) d/a = 1.0. For each of these three cases, the inaccurate mode shapes are shown above the accurate mode shapes.

together with the non-dimensional eigenfrequencies and the required value of N. Note that N = 18 was required for computing the sixteenth mode accurately. 6. Simulation of constrained membrane dynamics 6.1. Computation of modal coefficients From the general form of the dynamics of the membrane given in Eq. (9), the complete solution can be written as follows

Zðr; h; tÞ ¼

1 X







uci kci r; h Aci cos xci t þ Bci sin xci t þ

i¼1

1 X



h

usj ksj r; h Asj cos xsi t þ Bsj sin xsj t

i

ð38Þ

j¼1

(a)

(b)

Fig. 6. First symmetric mode and sixteenth symmetric mode of a constrained membrane with a/R = 0.1 and d/a = 9.0. For both these cases, the isometric view is placed above the top view and the value of N required for accuracy is shown.

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168 c

161

s

where uci , usi are defined in Eqs. (19) and (23), respectively, and xci ¼ cki , xsj ¼ ckj ; i; j ¼ 1; 2; . . .. The initial conditions are assumed to be

Zðr; h; 0Þ ¼ Gðr; hÞ;

_ h; 0Þ ¼ Hðr; hÞ Zðr;

Applying the displacement initial conditions to Eq. (38), we get 1 1 X X Aci uci þ Asj usj ¼ Gðr; hÞ i¼1

ð39Þ

j¼1

Multiplying Eq. (39) by uci , integrating both sides over the unconstrained area of the membrane, and using the orthogonality condition in Eq. (36), we get

Aci ¼

 Z 2p Z 0

rjq¼R

a

Gðr; hÞuci r dr dh

 Z 2p Z 0

rjq¼R



a

uci

 r dr dh ;

2

i ¼ 1; 2; . . .

ð40Þ

Multiplying Eq. (39) by usi , integrating both sides over the unconstrained area of the membrane, and using the orthogonality condition in Eq. (36), we get

Asj ¼

 Z 2p Z 0

rjq¼R

a

Gðr; hÞ usi r dr dh

 Z 2p Z 0

rjq¼R



usi

a

2

 r dr dh ;

j ¼ 1; 2; . . .

ð41Þ

Applying the velocity initial conditions to Eq. (38), we get 1 1 X X Bci xci uci þ Bsj xsj usj ¼ Hðr; hÞ i¼1

ð42Þ

j¼1

By repeating the same procedure as above, the remaining modal coefficients can be obtained as

Bci ¼ Bsj ¼

 Z 2p Z

1

x

c i

0

 Z 2p Z 1

xsi

0

rjq¼R

a rjq¼R a

Hðr; hÞ uci r dr dh Hðr; hÞusi r dr dh

 Z 2p Z 0



a

 Z 2p Z 0

rjq¼R

a

uci

2

 r dr dh ;

 rjq¼R   usi 2 r dr dh ;

i ¼ 1; 2; . . . j ¼ 1; 2; . . .

ð43Þ ð44Þ

The integrals in Eqs. (40), (41), (43) and (44) are of the form

Z 2p Z 0

gðhÞ

f ðr; hÞdr dh

a

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 where gðhÞ ¼ rjq¼R ¼ R2  d sin h  d cos h. Through the change of variables x = h and r = [(1  y)a + y g(x)] the integrals are changed to the form

Z 2p Z 0

1

f ½ð1  yÞa þ ygðxÞ; x½gðxÞ  ady dx ¼

Z 2p Z 0

0

1

hðy; xÞ dy dx

0

The integrations are performed numerically using a n-point Gaussian quadrature rule as follows

Z 2p Z 0

0

1

hðy; xÞ dy dx ¼

p 2

Z

1

1

Z

1

hðf; gÞ dg df ¼

1

n X n pX

2

wi wj hðfi ; gj Þ

ð45Þ

i¼1 j¼1

where x = p(g + 1), y = (f + 1)/2, wi and wj are the weights, and fi and gj are the integration points respectively. Since Bessel functions are very sensitive to numerical integration, accurate computation of the modal coefficients requires correct choice of the number of integration points in Eq. (45). An algorithm for properly choosing n and accurately computing the modal coefficients is described by the flow chart in Fig. 7, which is self-explanatory. 6.2. Numerical results verifying orthogonality In Section 5.2 we showed that results in the literature provide reasonably accurate eigenfrequencies but the associated mode shapes are inaccurate and as such do not satisfy the boundary conditions. In this section we show that these inaccurate modes do not satisfy the orthogonality condition; the accurate modes satisfy the orthogonality conditions but require proper choice of the number of integration points n. The number of integration points can be chosen properly by following the algorithm described by the flowchart in Fig. 7. We consider constrained membranes with a/R = 0.5 and d/a = 0.6, 0.8, 1.0; these cases were presented in Table 2. For this membrane, the third symmetric mode uc3 and fourth symmetric mode uc4 are used to check for orthogonality. In particular, we check the orthogonality of uc4 (accurate mode shape) with both uc3 (accurate mode shape) and uc3 (inaccurate mode shape). The results are presented in Table 3. To show the orthogonality between symmetric–symmetric and symmetric–antisymmetric modes, we consider a membrane with a/R = 0.1 and d/a = 0, 4, 8. The results are shown in Table 4. It is clear that accurate computation of the mode

162

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

Fig. 7. Algorithm for computing modal coefficients.

Table 3 Checking orthogonality between the accurate and inaccurate third symmetric mode and the accurate fourth symmetric mode of constrained membranes with a/R = 0.5 and d/a = 0.6, 0.8, 1.0. The inaccurate third symmetric mode is underlined. d/a n

0.6 D

4 8 12 16 20 24 28 32

0.0091 0.0366 0.0142 0.0056 0.0013 0.0015 0.0016 0.0016

c 3;

c 4

u u

E



uc3 ; uc4



0.0149 0.0821 0.0080 0.0081 0.0011 0.0000 0.0000 0.0000

0.8 D

c 3;

c 4

u u

E

0.0221 0.0621 0.0069 0.0187 0.0080 0.0061 0.0068 0.0067



uc3 ; uc4



0.0465 0.1769 0.0514 0.0427 0.0184 0.0025 0.0002 0.0000

1.0 D

uc3 ; uc4

E

0.0560 0.0171 0.0107 0.0026 0.0030 0.0019 0.0021 0.0022



uc3 ; uc4



0.3005 0.1715 0.0327 0.0222 0.0244 0.0058 0.0006 0.0000

Table 4 Number of integration points n required to achieve orthogonality between symmetric–symmetric and symmetric–antisymmetric modes of constrained membranes with a/R = 0.1 and d/a = 0,4,8. d/a n

0 

4 8 12 16 20 24 28 32 36

0.0013 0.0013 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

c 4;

c 5

u u





c 4;

s 5

u u

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000



4 

c 4;

c 5

u u



0.3009 0.2686 0.1659 0.0741 0.0086 0.0004 0.0000 0.0000 0.0000



c 4;

s 5

u u

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000



8 

uc4 ; uc5



0.0746 0.0215 0.0293 0.0182 0.0224 0.0100 0.0020 0.0002 0.0000



uc4 ; us5

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000



A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

163

shapes requires proper choice of the number of integration points n. The number of integration points depends on the size and location of the constraint and can be chosen using the flowchart in Fig. 7. It may be mentioned that symmetric–antisymmetric modes are orthogonal for small values of n and higher values of n are required to satisfy orthogonality between symmetric–symmetric and antisymmetric–antisymmetric modes. 6.3. Simulation of free vibration We simulate the motion of a membrane subjected to different constraints and initial conditions. The membrane is assumed to have the following specifications:

l ¼ 0:25 kg m2 ; T ¼ 1 N m1

R ¼ 1 m;

Three cases are simulated using 16 symmetric and 16 antisymmetric modes. The constraints and initial conditions for these three cases are described below:  Case 1: The membrane has nonzero initial displacement, zero initial velocity, and the constraint is located at the boundary – see Fig. 8. The constraint location and the initial displacement are described by the expressions:

a=R ¼ 0:12;

d=R ¼ 0:88

8

Gðr; hÞ ¼ 0:1ðr  aÞ ½r  rðhÞ;

Hðr; hÞ ¼ 0;

rPa

 Case 2: The membrane has nonzero initial displacement, zero initial velocity, and the constraint is located at some distance from the boundary – see Fig. 9. The constraint location and the initial displacement are described by the expressions:

a=R ¼ 0:12; d=R ¼ 0:66 3 Gðr; hÞ ¼ ðr  aÞ½r  rðhÞ8 ; Hðr; hÞ ¼ 0; 5

rPa

 Case 3: In this case, the membrane is initially assumed to be unconstrained and vibrating in its first mode. As it passes through the mean position, a constraint with the following specifications:

a=R ¼ 0:1;

d=R ¼ 0:4

is applied and held fixed for all future time. The position and velocity of the unconstrained membrane are given by the expressions

(a)

(b)

Fig. 8. Initial displacement of the membrane in Case 1: (a) Isometric view. (b) Side view.

(a)

(b)

Fig. 9. Initial displacement of the membrane in Case 2: (a) Isometric view. (b) Side view.

164

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

(a)

(b)

Fig. 10. Snapshots of the vibrating membrane in Case 1 over 1 s of its motion. Time progresses from the top snapshot to the bottom snapshot in column (a) and then in column (b).

165

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

(a)

(b)

(c)

Fig. 11. Snapshots of the vibrating membrane in Case 2 over 1 s of its motion. Time progresses from the top snapshot to the bottom snapshot in column (a), then column (b) and column (c).

 c  Zðr; h; tÞ ¼ J 0 k1 q cos xc1 t;

  _ h; tÞ ¼ xc J kc q sin xc t Zðr; 1 1 0 1

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 c wherepqffiffiffiffiffiffiffiffi , ffi r2 þ d þ 2rd cos h is defined in Fig. 1, k1 R ¼ 2:405 and is obtained from the first zero of the Bessel function J 0 , c c x1 , T=qk1 , and t = 0 denotes the initial time when the membrane is released from rest. At t ¼ p=2xc1 , the membrane passes through the mean position and the constraint is applied. The initial conditions for the constrained membrane are therefore

166

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

Fig. 12. Snapshots of the vibrating membrane in Case 3 during one cycle of its vibration.

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

( Gðr; hÞ ¼ 0;

Hðr; hÞ ¼

167

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  c  2 2 xc1 J 0 k1 q : a < r 6 R2  d sin h  d cos h 0

:06r6a

Fig. 10 shows snapshots of the membrane over 1 s of its motion for Case 1; both isometric and side views are included. Fig. 11 shows the motion of the membrane for Case 2. In Figs. 10 and 11, time progresses from the the top snapshot to the bottom snapshot in each column, and from the left column to the right column. Fig. 12 shows the motion of the membrane for Case 3; the isometric view is presented together with two side views. The simulations show plausible membrane motion including reflection of waves at the boundaries. Both internal and external boundary conditions are satisfied and this is indicative of the accuracy of our simulations. Prolonged duration video animations of all three cases are included as supplemental material. To the best of our knowledge, such simulations of constrained membranes have not been presented in the literature. 7. Conclusion We investigated the dynamics of a circular membrane with an internal eccentric circular areal constraint. The membrane is assumed to be fixed at its outer boundary and the areal constraint is assumed to impose zero displacement over the entire internal circular area. This problem has been investigated by several researchers but prior efforts have been limited to solving for the first few eigenfrequencies; the orthogonality property of the modes have not been established and the procedure for computing the mode shapes or simulating the dynamics has not been presented. In this paper we establish the orthogonality property of the symmetric and antisymmetric modes of the constrained membrane and provide a systematic method for computing the eigenfrequencies and mode shapes accurately. It is shown that the number of terms in the series expansion plays a critical role in accurate computation of the mode shapes. While fewer terms result in truncation errors, too many terms lead to numerical errors due to high sensitivity of the Bessel functions. An algorithm is presented to choose the appropriate number of terms and the importance of the algorithm is demonstrated through examples where fewer or more terms result in accurate eigenfrequencies but inaccurate mode shapes. The orthogonality conditions established in this paper are used to determine modal coefficients for dynamics simulations. This however requires proper choice of the number of numerical integration points. An algorithm is presented to determine the appropriate number of numerical integration points which is shown to depend on the size and location of the areal constraint. Using the two algorithms for determination of the mode shapes and computation of the modal coefficients, the dynamics of a constrained membrane is accurately simulated; the accuracy of the simulations can be verified from the boundary conditions. The tools developed for simulation of constrained membrane dynamics, presented in this paper, can be used to study the vibration of membrane structures and the efficacy of control methodologies for vibration mitigation. Acknowledgement The authors gratefully acknowledge the support provided by Air Force Office of Scientific Research, AFOSR Grant FA955010-1-0500, for this work. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.simpat.2012.10.008. References [1] M. Sosa, A.A.O. Carneiro, J.F. Colafemina, O. Baffa, A new magnetic probe to study the vibration of the tympanic membrane, Journal of Magnetism and Magnetic Materials 226–230 (2) (2001) 2067–2069. [2] E.J. Ruggiero, D.J. Inman, Modeling and vibration control of an active membrane mirror, Smart Materials and Structures 18 (2009) 095027. [3] K. Miyano, Surface tension measured by vibrating membranes: an application to smectic-A and -B phases, Physical Review A 26 (3) (1982) 1820–1823. [4] G. Percin, T.S. Lundgren, B.T. Khuri-Yakub, Controlled ink-jet printing and deposition of organic polymers and solid particles, Applied Physics Letters 73 (16) (1998) 2375–2378. [5] H.Y. Yee, N.F. Audeh, Uniform waveguides with arbitrary cross sections considered by the point-matching method, Institute of Electrical and Electronic Engineers Transactions MTT13 (1965) 847–851. [6] H.Y. Yee, N.F. Audeh, Cut-off frequencies of eccentric waveguides, Institute of Electrical and Electronic Engineers Transactions MTT14 (1966) 487–493. [7] J.B. Davies, C.A. Muilwyk, Numerical solution of uniform hollow waveguides with boundaries of arbitrary shape, Proceedings of the Institute of Electrical and Electronic Engineers 113 (1966) 277–284. [8] C.W. Steele, Numerical computation of electric and magnetic field in a uniform waveguide of arbitrary cross section, Journal of Computational Physics 3 (1968) 148–153. [9] P.L. Arlett, A.K. Bahrani, O.C. Zienkiewicz, Application of finite elements to the solution of Helmholtz’s equation, Proceedings of the Institute of Electrical and Electronic Engineers 155 (1968) 1762–1766. [10] N. Gass, Solution of the three dimensional scalar Helmholtz equation by a finite element formulation, ZAMM – Journal of Applied Mathematics and Mechanics 57 (2) (1977) 95–99. [11] P.A.A. Laura, A simple method for the calculation of cut-off frequencies of waveguides with arbitrary cross section, Proceedings of the Institute of Electrical and Electronic Engineers 54 (1966) 1496–1497.

168

A. Alsahlani, R. Mukherjee / Simulation Modelling Practice and Theory 31 (2013) 149–168

[12] M.J. Hine, Eigenvalues for a uniform fluid waveguide with an eccentric-annulus cross-section, Journal of Sound and Vibration 15 (3) (1971) 295–305. [13] R.M. Bulley, J.B. Davies, Computation of approximate polynomial solutions to TE modes in an arbitrarily shaped waveguide, Institute of Electrical and Electronic Engineers Transactions MTT17 (1969) 440–446. [14] K. Nagaya, Vibration of a membrane having a circular outer boundary and an eccentric circular inner boundary, Journal of Sound and Vibration 50 (1979) 545–551. [15] W.H. Lin, Free transverse vibrations of uniform circular plates and membranes with eccentric circular holes, Journal of Sound and Vibration 81 (3) (1982) 425–435. [16] G.S. Singh, L.S. Kothari, On the solution of the two-dimensional Helmholtz equation, Journal of Mathematical Physics 25 (1984) 810–811. [17] G.S. Singh, A. Vishen, G.P. Srivastava, L.S. Kothari, Comments on free transverse vibrations of uniform circular plates and membranes with eccentric holes, Journal of Sound and Vibration 100 (1) (1985) 141–145. [18] K. Nagaya, Y. Hai, Free vibration of composite membranes with arbitrary shape, Journal of Sound and Vibration 100 (1) (1985) 123–134. [19] P.A.A. Laura, R.H. Gutierrez, Determination of the fundamental eigenvalue of a rhombic membrane with a concentric circular hole, Journal of Sound and Vibration 153 (1) (1992) 176–179. [20] L.H. Wu, C.Y. Wang, Fundamental frequencies of a circular membrane with a centered strip, Journal of Sound and Vibration 239 (2) (2011) 363–368. [21] S.W. Kang, J.M. Lee, Free vibration analysis of composite rectangular membranes with an oblique interface, Journal of Sound and Vibration 251 (3) (2002) 505–517. [22] J.T. Chen, S.K. Kao, W.M. Lee, Y.T. Lee, Eigensolutions of the Helmholtz equation for a multiply connected domain with circular boundaries using the multipole Trefftz method, Engineering Analysis with Boundary Elements 34 (5) (2010) 463–470. [23] V.H. Cortinez, P.A.A. Laura, Vibrations of non-homogeneous rectangular membranes, Journal of Sound and Vibration 156 (2) (1992) 217–225. [24] F.F. Cap, Eigenfrequencies of membranes of arbitrary boundary and with varying surface mass density, Applied Mathematics and Computation 124 (2001) 319–329. [25] C. Zhao, G.P. Stevens, A practical error estimator for finite element predicted natural frequencies of membrane vibration problems, Journal of Sound and Vibration 195 (5) (1996) 739–756. [26] B. Stanford, P. Ifju, R. Albertani, W. Shyy, Fixed membrane wings for micro air vehicles: experimental characterization, numerical modeling, and tailoring, Progress in Aerospace Sciences 44 (2008) 258–294. [27] A. Song, X. Tian, E. Israeli, R. Galvao, K. Bishop, S. Swartz, K. Breuer, Aeromechanics of membrane wings with implications for animal flight, AIAA Journal 46 (8) (2008). [28] P. Rojratsirikul, M.S. Genc, Z. Wang, I. Gursul, Flow-induced vibrations of low aspect ratio rectangular membrane wings, Journal of Fluids and Structures 27 (8) (2011) 1296–1309. [29] D.-A. Wang, H.-H. Ko, Piezoelectric energy harvesting from flow-induced vibration, Journal of Micromechanics and Microengineering 20 (2) (2010). [30] J.J. Allen, A.J. Smits, Energy Harvesting Eel, Journal of Fluids and Structures 15 (3) (2011) 629–640. [31] P.A. Tarazaga, M.E. Johnson, D.J. Inman, Vibro-acoustics of a pressurized optical membrane, Mechanical Systems and Signal Processing 30 (2012) 373– 392. [32] A. Arokiasamy, Simulation of vibrations of rectangular and square membranes using computer graphics, Simulation Modeling Practice and Theory 10 (1–2) (2002) 3–12. [33] E. Alvarado-Anell, M. Sosa, M.A. Moreles, Computational study of forced oscillations in a membrane, Revista Mexicana de Fisica 51 (2) (2005) 102–107. [34] E. Alvarado-Anell, M. Sosa, M.A. Moreles, Numerical simulation of the dynamical properties of the human tympanum, Revista Mexicana de Fisica 54 (2) (2008) 135–140. [35] A. Bhide, Computer simulation of human tympanic membrane, The Journal of Laryngology and Otology 106 (1992) 878–881. [36] E. Kreyszig, Advanced Engineering Mathematics, ninth ed., John Wiley and Sons, Hoboken, NJ, 2006. [37] G.N. Watson, A Treatise on the Theory of Bessel Functions, Cambridge University Press, 1958.