Two degree of freedom numerical simulation of faraday waves

TCcfi BVET

abbraum voor

Stheopshydrmechanlca

Aichief

Mekelweg 2,2828 CD Delit

TWO DEGREE-OF-FREEDOM NUMERICAL SIMULATION OF FARADAY WAVES

by

Cbñstopher Wayne Carter

A Thesis submitted in partial fulfillment of the requirements for the degree of

Master of Science

(Naval Architecture and Marine Engineering) at the University of Michigan

1998

Thesis Advisor

Iwould like to dedicate this thesis toy mother, for without her I would have never accomplished this much,

and to my soon-to-be-Wife, Karlie, who always supported me while I was at graduate school,

even through the toughest of times, andencouraged

ACKNOWLEDGMENTS

I would like to thank my thesis advisor, Associate Professor Dale 0._{IaiT His guidance} and patience in completing this thesis was invaluable _{I can say, with out a doubt, that without his} support, this thesis would have never been completed In addition, I would like to thank Professor Marc Perhn for his assistance with this thesis He was my main source of information on the Faraday wave expenments Without his help, much of the insight about the prior research would have gone unnoficed.

I would like to thank the faculty and staff of the Department of Na'al Architecture and Marine Engineering who all had some impact on my time at the University of Michigan I would also like to give a special thanks to Ms Colleen Vogler Without hersuperior administrative

skills, I would have never gotten anything done while in graduateschool.

Finally, I want to thank my family and friends whO supported me during my graduate

TABLE OF CONTENTS DEDICATION 2

ACNOtIEDGME1S

.1

3 LIST OF FIGURES 5 ABSTRACT 8

CHAR

L INTRODUCTION 9

U. SYSTEM MODEL AND DERIVATION OF EQUATIONS 10

rn. NUMERICAL ANALYSIS OF EQUATIONS OF MOTION 22

3.1 State space representation 22

32 Fortran program 25

IV CALCULATIONS OF PARAMETERS 29

4.1 Mass and rod length - 30

4.2 DerivatiOn of equations of motion Using linear theory ... 31 4.3 Linearization of non-linear equations of motion ... 34

4.4 Calculations of cOnstants 38

V. NONLINEAR SIMULATION RESULTS 49

V CONcLUSION - 57

APPENDDCA 58

LIST OF FIGURES

Figure

1 Sketch of Faraday wave with proposed model

2 Sketch of three degree-of-freedom model of system

3 _{Sketch of simplified two degreeof-freedom model of system} 4 _{Sketch of forces acting on the model used to find the moment about}

the pin connection R

Sketch of forces acting on the model used to find the moment about the first point mass M3

6 Sketch of the di.tnensions of tank

7 _{Sketch of model showing coordinate transfOrmation} 8 Graph of the linear response of system of first mode shape

0= l.5717963268rad,$=1.5707963268rad, O =4)

=0

9 _{Graph of the linear response of system of second. mode shape}

0= 1.5707963268 r _{$= 1.5717963268} O _{= 4}

=0

10 _{Graph of the linear response of system with both M1 and M2 displaced}

0=1.5687963268rad,4)=1.57l7963268rad, O ₌

=0

11 _{Graph of damped sine wave With envelope curve} (D=lOandI3=O.3/sec

12 _{Graph of the Wave-peak decay with least-square approximation}

0.3Isec

13 Graph of linear response of M1 with envelope curve 0= l.5717963268rad,$= 1.5707963268rad,

₌₄

=0

44 45 46 46 47 pg. # 10 11 12 12 12 30 31

14 Graph of the wave-peak decay with least-square approximation

= 0.05/sec 47

15 Graph of System response wIth period one standing wave 0= 1.9394222 tad, 4, = 0.75032909 tad, O = 7.6655593 tad/sec.

4) = 7.3153634rad/sec, amp = 0.006 meters, CO = 10.2324844 rad/sec 50

16 Graph of phase plane of system period one

o = 1.9394222 tad, 4, = 0.75032909 rad, 0 = 7.6655593 rad/sec,

4) = 73153634 ra /sec, amp =0.006 meters, w= 10.2324844 tad/sec 50

17 Graph of system response, from rest, of period one standing waves o 1.5707963268 rad, 4,= 1.5707963268rad, O =4)

=0,

amp = 0.00265 meters, w = 10.053 rad/sec 51 18 Graph of system rsponse, from.rest, of period one standing wave

from 190 sec to 200 sec: 0= 1.5707963268 rad, $ = 1.5707963268 rad,

O=4)=0,amp=0.Oo265meters,co=1O.O53radIseC 51

19 Graph ofphase plane of system, from rest, from 190 sec to 200 sec:penod one 0= 1.5701963268 rad, 4, = 1.5707963268 rad, O _{= 4)} = 0,

amp = 0.00265 meters, co = 10.053 rad/sec 51 20 Sketch of dimpled standing wave with idealized model 52 21 Graph of system response of period one dimple

o = 1.5707963268 rad, 4, = 1.5707963268 tad, 0 = 4) = 0,

amp 0.003705 meters, o 10.053 tad/sec

::

52

22 Graph of system response with period three standing wave o = 1.9817702 rad, 4, = 0.66793155 rad, é = 7.6345674 rad/sec,

4) = 7.2529606 rail/sec. amp 0.006 meters, w = 10.18336536 rail/sec 53 23 Graph of phase plane of system: period three

o = 1.9811702 rad, $ = 0.66793155 ra, 0 = 7.6345674 rail/see,

4) = 7.2529606 rad/sec, amp 0.006 meters,0 = 10.18336536 rail/sec

53,

24 Graph of system response, from rest, of penod three standing wave,

0= 1.5707963268rad,4,= 1.5707963268rad, 0 = 4)

=0,

25 Graph of system response, from rest, of period three standing wave from 240 sec to 250 ec: 0= 1.5707963268 rad, = 13707963268 rad,

O==0,arnp=O.O038meters,co1O.O53rad/sec

26 Graph of phase, plane of system, from rest, from 240 sec to 250 sec: period

three 0= 1.5701963268 rad, $ = 1.5707963268 rad, 6 =

=0,

amp = 0.0038 meters, o = 10.053 tad/sec =

27 _{Graph of system demonstrating breaking Waves for height of M1} 8= 1.5707963268 tad, = 1.5707963268 rad, O = 4, =

amp=0.0045 metes, (0= I0.OS3rad/sec

28 Gtaph of system demonstrating breaking Waves for angle of O 0= 1.570796328 rad,4 = 1.5707963268 rad, = 4, =0, amp = 0.004 neters, o = 10.053 rad/sec

29 Graph of phase plane of system demonstrating braking waves (0 vs. 0)

0= 1.5707963268rad,4=1.5707963268rad, 0 =4,

=0,

amp = 0.0045 neters, (0= 10.053 rad/sec

54 54 55 55 55 3

ABSTRACT

Previous research at the University of Michigan included experiments with Faraday waves. Researchers were able to excite the system and find both period one and period three breaking standing waves. Further excitation led to increased breaking waves. During that research, questions arose as to the behavior that the Faraday wave was exhibiting. One method to help explain that behavior is to create a numerical model of the Faraday wave. The purpose of this thesis is to create a discrete model of the system and test it to ensure that the model has the same characteristics as the Faraday wave in the experiment. A simplified two degree-of-freedom

model is designed and its differential equations of motion derivecL The equations Of motion are

programmed in Fortran '71 and solved with the aid of the "isoda" software package. Values for constant parameters are calculated using a linear model of the system and linear response data from the experiment. The model is successfully excited to both a period one and period three

standing wave.

CHAPTER I

INTRODUCI1ON

Previous research at the University of Michigan included laboratory expenments on

Fara-day waves. In a study by Jiang, Ting, Perlin and Schultz [1.], a low aspect ratiotanic, flllçd with water, was oscillated vertically. At a certain frequency and amplitude, the water began to form a standing wave of the first subharmonic. The first subhaflnomc wave length wasthe length of the tank This wave was found to be a.period one Standing wave (11. In a second paper by Jiang, Per-li and Schultz [2], the researchers found that by increasing the ampPer-litude and changingthe fre-quency slightly, they were able to find a period three standing wave. After increasing the amplitude further, the wave began to break and exhibit chaotic behavior (2]. This is not consid-ered a northal route to chaos. The normal route to chaos is through period doubling [3]. So the question arises, is the route through period three evidence that this System may becomechaotic. In the se ond paper, the authors stated that they wanted to çxplore numerical methods to qualita-tively describe period tripling [2].

The purpose of this thesis is to create a non-linear discrete spring-mass-dashpot model that

will simulate the behavior of the Faraday wave. The idea is to cteate a reduced finite degree-of-freedom tractable model that will make this complex system easier to analyze. Thismodel would

then help explain some of the physics of the Faraday wave and maybe explain how period tripling

could be. a possible route to chaos. Though the final goal of all the research is to describe the physics of the system, this thesis deals with creating the model and fliIcing it match the behavior of the Faraday waves. This includes the period one and period three response of the system. Fur-ther research is necessary to better understand the Faraday vrave.

0.5

-0

0.5

1

-CHAPTER II

SYSTEM MODEL AND DERiVATION OF EQUATIONS

For the purposes of thi thesis, the simplest discrete model was sought which could cap-ture the salient feacap-tures of the actual system. It is known that the type of wave of interest is a

standing wave at, the fi±st sub.harmothc frequency. For a standing wave, there are two nodes that

will result from the resonant frequency. Becatise of this, the system will be modeled by using only the section of the standing wave between the two nodes. Figure 1 is a diagram of an ideal standing wave used as a visual example. Also in Figure 1 is a rough sketch of the systembeing considered for simulation to capture important results from the original Farathy experiments of references [1] and [2]. Superimposed on the standing wave is a three degree of freedom

oscillat-ing system.

0 1 2 3 4

Figure 1

Sketch of Faraday wa'e with proposed model

The Faraday wave can be modeled approximately using this three degree of freedom spring-dashpot system. The system consists Of three point masses in conjunction with a spring

and damper for each mass. The outside edges of the model are connected using pin connections, which are the nodes as indicated. This will ensure that the pin connections do_{not move relative to} the base. The system is forced using abase excitation. The forced motion of the base,flz),is

rep-resented using a sinusoidal wave with a constant amplitude and frequency. Figure 2 is a more detailed sketch of the model being considered.

-

_____

_{- _/t}

_{- - _____}

ftt)J

Figure 2

Three degree-of-freedom model of system

The three point masses are connected by massless rods Each point mass is of equal weight and has a respective spring coefficient,A; and damping coefficient, c. The length of the total system is equal to half the length of the tank Advantage is taken of the _{symthetry of the} sys-tem to create a two degree-of-freedom simplified syssys-tem. Figure 3 is a sketch of the_{new and final} system that will be smulated.

The length of the system when horizontal is 1, which is equivalent to a quarter of the length

x2

Figure4

Sketch of forces ting on the model Used to

flid the moment about the pin connection R

Figure 3

Simplified two degree-of-freedom model of system

same as they were with the full system. The valtiesof m2, k2 c2ar considered to be half of the original values m order to rnalntain the use of symrnetly. In oilier worth, when the model was cut in half, thó Values for the middle components werealso cut in half. TO compensate for the

effects of the right part of the system, k3 and C3 are added.

Figure 5

Sketch of forces acting on tl model used to find the mnt about the first point mass M1

Two equations of motion are needed to properly model a two degree of freedom_system. The moment about the pin connection R and the moment about the first pointmass, m1, are used to find the equations of motion. A simplified model of the system is shown in Figure 4 which will

help in finding the moment about R; Figure 4 shows the forces acting on the model that are

rele-vant to determining the moment about the pinned connection, R.

Figure 5 shows a sub-section model used in determining the moment about m1. Figure 5

is a sketch of the system that shows the forces acting on the model that are relevant to determining

the moment about the flit pointmass, m1.

Before the equations of motion are determined, the following kinematic relationsare

defined. Referring to Figure 3 again, x1 is the displacement ofm1 in the x direction relative to the

pinned connection. The lisplacement is related to the angle e by

=

The velocity and acceleration are calculated by finding the first and second derivative ofx1 with respect to time (as indicated by x and 1 respectively):

1. = OcosO

=

_{[ésin(0)+cos(0)]}

Also, y is defined as the displacement of the first point_{mass, m1, with respect to the base of the} system in the y direction. The displacementy , the velocity and the acceleratiOn $7 are defined respectively as follows:

Yi = cosO

yI 1[O2cos(0)+Osin(0)]

Finally, z1 is defined as the vertical displacement of in1 relative to afixed datum. The position of

the base is given by the displacement forcing function,itt). Hence, z1 and its derivatives are

defined as follows:

z1 =

f(t)+y1

=

The relative acceleration for mass One in the y direction1 j' , wasdefined above. Substituting this

into the acceleration gives

f(t)_(è2cosO+OsmO]

Next the displacement variables for mass two, in2, will be defined. x2 is the displacement

of the point mass,_rn2,in the x direction relative to the pin contion. x2 is equal to the displace-ment of the fit mass, x1, plus the horizontal distance betweenthe two masses.

1.

1

+sm$= (sm8+sin)

The velocity and acceleration are found by calculating the first and second derivative of the

dis-placement.x2.

Also,_{Y2 is the vertical displacement of m2 with respectto the base of the system. Y2 is equal to the}

displacement of the flit mass, Yi plus the vertical distance between m1 and m2. The relative dis. placement velocity and accelerationare given as foI1Ows

Y2

y1+cosO=(cos6+cos)

92 =

[Osin9-4sin]

=

The last definition needed is_{z2 which is the veitical displacement between m2 and the fixed}

datum. z2 is equal the displacement of the secOnd_{mass with respect to the base,Y2 plus the}

dis-tance from the datum to the base. z2 and its derivatives are thus as follows.

z2 =f(t)+y2

Z2 = f(r)+92

= /(t)+y2

Substitution of the relathre acceleration of_{mass two, y2, into the acceleration gives the equation:}

12

f(t)=[2cosO+OsinO+$2cos$+sm]

With all the necessaiy variables defined, the two differential equations_{of motion can be}

determined. The first equation is found using Figure 4 as a guide. The first equation is

deter-mined by Summing the moments about the pinned connection, R.

= 0 = x1[m1z1 +m1g +k1(y1 y)+c191]y1m1X1 +x2[rn212+m2g '

₍₁₎

+ k2(y2

-

+ c,j2] y2[m212 + k3(x2 x) + c3x]

moments about the point mass m1.

Mm = 0 =

c4(_O)+k4(4'_O°)+(y2yi)[m2I2+k3(Xz4)+C32]

(2)

(x2x1)(m212 +

m2g

+ k2(y2 y) +

c22]

These are the two differential equations of motion that will be used to model the system. The next step is to substitute the values of the displacements, velocities and accelerations that

Were determined from above into the equations of motion so as to representthe model in terms of

9 and 4'and their derivatives1 Substituting these values into equation (1), gives the following:

sin8[mi(j(t)_(O2cos9

₊ sin9))+ mg +

ki(cos9 p;)

+

ci_ôsiñO)]

-

(3)

cose(m1)[O sinO +

OcosO]

+ (sm9 + sin$)

[m2(f(:) - (O2cosO + sin9 + 42cos$ + sin,)) +m2g +

k2((cosO

+ cos4') y;)

+ c2(Osm9

-

4sin4')]_ (cosO + cos$)[m:z(_O2sin9 + cosO

+cos4')

+

k3((sinO

+ sin$)

+ c3(cosO +

4'cos$)] o

After algebraic mmipu1atiOn, the first equation of motion can be rewritten:

01(mi +m2)+

m2(sin4sinO +

cos$cosO)+2m2(sinOcos4'_

cosOsin4')J(t)

(4)

sinO(m1 + rn2) - gsin9(m1 + m)k1

sinO(cosOY) k4sinO(cosO

+ cos4') + Ic2

.v sinO +

k3çcoso(sine

+ sin4') k3 x cosO +

cifO(sinO)2

+

c2(O(sinO)2

_{4'sinesin4')+ c3f(O(cose)2}

+$cosOcos$)

=

_[m4sinOsin

+ cosOcos$) + 2m2(cosOsin4'

- cos$sinO) + m-

J(t)m2

Sin4' - gm2sin$ - k2çsin,(coso

+ cost) + k2

v sin$ +

k3fcos$(sinO + sin$)

Similarly from equation (2) we find:

C4(4 O) ₊k9(4 =

8

0)+

((cos8

+

cos$).

cos9) ₍₅₎

[m2(o2sino

_+OcosO

_-

_{2sin+ cos)k3((sinO}

₊

_sm$)

_x;)

+ c(OcosO

+

4cos)]

-

((smO

+

sink) -

sinO)

[m2(f(t)

(O2cos8+ sm8+ cos +

si$])

+ m2g

+

k2((cos9

+

cos) y

)+

c4(=OsinO

sin)]

= o

Again, common algebraic expressions are collected to yield:

c4( - O)k9(Ø

°)

=O(m cos$cos8 ₍₆₎

+O2(=m22cossin8+

nzsincosO)+ (m2(cos$)2

+

m2(sm$)2)

+

2(_m2cosqsin

_{+rn2çsinscos.)+}

_{f(t)(_msin$)}

₊

_{g(_m2sin)}

k2(_sin(cosO

+

cos.))

+ k2

y sine

+

k3cos(sinO

+

sin)

1c3 xcos

+

c4OsinOsinl)

+

4(Sin)2)

+c3(OcosOcos4 +

(cos)2)

Notice that the nght hand side of equation (4) s equal to the right hand side of equation (6). Equation (6) an therefore be substituted into equation_{(4) Also the trigonometric identities}

of [sin( -0)

= cosOsin$ cos4säiOJ and [cos( -0) =cos4,cos0+sin$sinO] can be Used to further

simplify the equations. Making these substitutions into equations (4) and (6) gives the following

+ m2) + _{m2cos( 0)}

02r.

-0)-f(t)sm9(,n1

+

m2)-gsm0

(7)

(m1 +m2)-/c1

sine(cose_YT)k2...Sifle(_...(COSe+cosØ)+

+k3cose(L(sifle+

sin$)_ r;)+k4(,_e +

O)

cO(sjne)2c2.L

(O(sinQ)2

_{+ )sin0sjn$) + c3(OcosOcos}

(cos$)2)

+ C4( -

= 0

(cos(,

-0)

m2+O2m2!sjn,

_-0)

_{- f(t)rn2sjn4)}

_-g?nsin$

+k2i (8) sinO(_

_{(cos0 + cos)}

₊

2)+

_{k3cos$(_(sm0}

_{+ sin$) -}

x)

_-0-

0)

+ c2sin*(0Sm9 + $sin)

_{+ c3cos$(9cos0 + $cos$) + c4($ -0)}

= 0

Next, the equations _{of mot on are simplified further for ease of}

use. A new variable M is

defined as the total_{mass of the system and the}

two point masses_{are flOimSI ze with respect to}

the total mass to_{create M1 and M2.}

M = nz1+m2

m1 andin2_{can therefore be defined as follows:}

m1 =

M1 M

?n2 =

hi2 M.

Substituting the definitions of m1 and in2 into_{equation (7) and(8)yields:} -.

and

efM

+ _fMM2cos(

9)

MM2sin(

_{- 0)f(MsinOgMsin8}

(9) - k1

_{siiio(cose_y;)+}

k2sin0(-(cos0

+

_cos)+

+

_k3cos8((sin9

_{+sin4,) -}

₎₊

k4(

- 9+ °) +

_{c1çO(sili0)2+ c2}

(O(sin0)2

+ sinOsin) +

_{c3f(é(cose)2 + cOs9cos$)}

+ c(+O)

=0

for the first equation_{of motion and,}

12 _.

_12.2

12 _..

_i

i

OMM2cos(d)

- 9)+ MM2

+0

_{MM2rifl( - 0)}

f(t)MM2sin*gMM2-2sin

(10)

+ k2.sin4(_ (cos0

_{+.cos) +}

₂₎

_{+ k3cOs$((sm0}

+ sink)

-

_x;)

_{+ k4($}

°)'

+ c2sin$(0 sinG

_{+ sin) + c cos(0cos9}

cos) + c ( 0)

= 0

4.

_.34

4 for the second_{equation of motion.}

Next multiplying_{equations (9) and(10)by}

j2 we get

+ M2cos($

_{- 0) 42M2sin($ - 0)f(t)sin9gsin0}

-

_{sin8(cos0_.)}

_{sin0(._}

(cosO + cos4i) +

+

_cos9(si9

_{+ sink}

_22.]+

8

_O)

+1è(sjfl0)2

+

_{sin0(sinO + 4sin)}

_{+ f-cos0(Ocos0}

+ 4cos) +

_{0) =0}

M

. M

-

Cl

c1=Ii

-

C2

c2=

f =

Asin(o)t)

-6M2cos( 0)

+ +O2M2sm(4

_{- e) - JO)M2sin$ gM2sin +}

(12) sin

_{(cosO + cos) +}

2L)

_{cos4{sin0+ sin_22.)+!2(,_e_o)}

+

_{sin$(ésine +4sin$)+}

_{cosq(écos0 + $cos$)+}

_{t(Q_ O)} 0

M _M

Mi2 To further simplify the_{equations of motion, the}

constant variables in the equation can be

norniliwj First, the spring

_{coefficients can be divided by thetotal mass:}

-

4k4.

K4---iF'

Mi2

The damping rates can be normalized with tespecuo the _{total mass by dividing the damping} val-ues by M:

-

4C4

c4=___

'.4 _Mi2

The forcing frequency is a sinusoidal wave in Wbih the amplitude_{and frequency can be}

speci-fied, f

_{Asin(wt). The forcing frequencyand the gravity term}

are nOrn11i7ed with respect to

the length of the massless bar that connects the two_{masses together. Dividing the forcing} fre-quency and gravity by this length_gives:

__2g

g The last normalizetj_{values to be defihed}

are those for the initial positions of the masses.

hpartic-ular

_{y, ,;}

, and _{These values are normalized with}

respect to the length of the massless bars.

=

2.

2

Substituting these definitjons into_{equations (11) and (12) yields the final}

differential equations_of motion that are used to model the _system:

o_{+ ØM2co(}

-

-2M2sm(

_{- O)f(t)sjn9sm}

1sinO(cose_y1)2sjflo(.., cosO Cos+2)

+K3cosO(sifle+sifl$)($0e_.2

C2SinO(sjn9 +4sin4)

_C3cos9(9cosO

_{+$cosØ)+ C4( 0)}

=

0

iJ2cos(4,-. O)M2

_{+ O2M2sin(?p_}

O)-1(t)M2sin_M2sin,+R2

Sifl(_COSO_COS+y2)+COS$(SjflO+S.._)A;

=

0

(13)

CHAFrER ifi

NUMERICAL ANALYSIS OF EQUATIONSOF MOTION

It would beveiy diffictilt, if not impossible, to find a closed form solution for the equations

of motion of the system being modeled. Becauseof this, a numerical simulation program to com-pute the response of the system is used. A powerful software package named "isoda"[4] is used

to simulate the system response. "Lsoth" is a computer simulation that solves diffçrential equa-tions numerically. A subroutine must be Written to properly call up the Isoda package. The pro-gram and subroutine are given at the end of this sçction. To use "isoda", the deferential equations of motion must be in state space form because Isçda can n1y solve first order differential equa-tions and the model has second order differential equaequa-tions. Therefore, the equaequa-tions of motion derived chapter 2 must be rewritten in state space form.

3.1 StatE space representalion

The state space representation is the mathematical representation of the system such that the state of the system is described by the Set offirst-order differential equations written in tenDs of a common state variables [5]. The equationsof motion derived in the previous sections is one mathematical representation of the dynamical system. However, "isoda" needs the equations of motion in state space representation. First recall the original equations of motion equations (13) and (14):

O + M2cos($

-' 0)

$2M2sin(4,

- 0)f t)sinejsinO

₍₁₅₎

X sin0(cos0

_{)+ siñ0( cosO. cOs +)}

A3cós9(sin9 + sn;) +k(o

_{+°)+ C1(sin0)2}

+ C2s1fl0(OsinO + sin) + 3cosO(Ocos4 +4cosØ)

_{+ C4( 0)}

₌₀

for the first equation and the second_{equation is:}

M2cos( 0) + $M2 + 2M2sin( 0) 7(t)M2sm

_{M2sin + K2}

₍₁₆

sin$(cOsOcos$+2) +K3cOs$(sin0+

_{sin$-;)-.4($-e-°)}

= 0 The newstate Space variables are defined as the following:

Yi 0.

Y2 =

y3 =$

y4 Taking the derivative of each gives:

91 = ₌ Y2

Y2

y3 = =

y4 =

The two equations, (15) and (16),_{can, be used to solve for 0 and . Using (15) and (16) pluS the}

y1 = Y2

=

[1

2][M2sin(y3_yi)(ycos(y3_y1)+y)

I M2(cos(y3y1))

+f(t)(M2cos(y3 y1)sin(y3)+ sin(y1))

+(M2COS(Y3_Yl)sin(y3)+.sm(y))_isin(y)(cos(y).)

K2(cos(y1)

cOS(y3)+92)(_Sin(y3)cos(y3_y1)+ sin(y1)) A3(sin(y1)+ sin(Y3);)(CO$(Y3)CO(3_1) cos(y1)) K4(y1y3+)(cos(y3_y1)+ I)-1(y2)(sin(y1))2

C2((y2)sin(y1) _{+ (y4)sin(y3))(sin(y3)cos(y3 y1)}

+ sn(y1)) C3((Y2)cOs(y1)_{+ (Y4)cos(y3))(cos(y3)cos(y3} y1) +

cos(y))

C4(y4+y2)(COs(y3_y1)+1)] y3 = y4 2 M2-M2(cos(y3 y1))

J(t)M2(cos(y3 y1)sm(y1)- sin(y3))

s _{(y3))) + R1M2 sin(y1)(=cos(y1)}

cos(y3 y1) +_K2(COS(y1)

-

_cos(y3)

+92)(M2sin(y1)cos(y3 y1) -

sn(y3))

+K3(sm(YI)+S1n(y3);)(M(Y)COS(YY))

+C2((y2)i(y1)₊

(Y4)sin(y3))(M2sin(y1)cos(y3y1) - sin(y3)) + C3((y2)cos(y1) +

(Y4)Cos(y3))(M2cos(y1)cos(y3._y1)_ cos(y3))

C C

C

These are the final _{equations, which are in state space form,}

that are Used in the "Isoda" package.

3.2

_{For an.}

_program

To use the isocla_{software, an intetface fortran progr}

m must be Written to call _{Isoda. This} prograni is given below.

The state space equations are entered into _{the Isoda software package} using a subroutine.

c2345678901234567890123456789012345

PROGRAM MAIN

C

c _{This program calls up the isoda}

c _equatjo

of motion for _{the two}

c model. _{The program uses a subr}

c _{of motion in}

state, space form.

C C

C to COmpile:

c f77 -cg89

-o faraday.exe faraday.f _{-L/usr/fth,lib}

-lcmix -lslatec _-lb].as

C

C

real*8 kone, ktwo, kthree,_{kfour, cone}

CtWo,ctee cfour

real*8 mone, mtwo _{length, gray,}

amp, w, xinot

implicit real*8 (a-h, _o-z)

dimension y(4), _{ato].(4), rwor]c(182)}

iwork(30) v(4)

C

coxon /great/ _{mone,mtwo, length,}

gray, anp,w,xot common /damp/ kone,ktwo

kthree, kfou.r,_{cone Ctwo, cthree}

four

external fcñ

C

neq = 4

software package to solve the

degree of freedom _{Faraday wave}

outjne that has _{the equations}

OPENI1, PILE= 'fa.raday.dat') OPEN(2, FILE=_{'Statespace out )}

OPEN(3,FILE=thetaajg.)

OPDi(4,ILEphi.glea)

OPEN(5SFILE= 'theta Phase'

°P6IPILE'phi.phasel)

OPEN(7,pi=_{'massl.heightu)}

OPEN(8.FILE1ss2'.height.)

OPEN(9DFILE=Imesstàtal.height)

OPEN(lO, PILE _{'massi .phase')}

ON(1l,FILE=Imess2phasel)

C C C C C C C C C C C

arid the initial conditions for the (v) vector

are given in v(i), 1=1,2,3,4 _{Here Cv) is}

determined from the equations _{of motion:}

(vdot) = [H]Cv} + (g) {f)

Use English units consistent _{with the input of}

ubroutine indat read (1, *) read(1, *) read(1, *) read(1, *) t = 0.OdO tout = 0..OdO itol = 2 rtol = 1.Od-12 atol(l) = 1.Qd-12 ato].(2) _1.Od-12 atol(3) = 1.Od-12 ató].(4) = l.Od-j2 itask = 1 istate = 1 iopt = 0 lrw = 182 liw = 30 jt = 2 do 40 iout = l,ntjines do -ideq = 1,neg y(ideq) _v(ideg) end do C C

1 _{iopt,rwork, lrw, iwork, liW,}

jdurn, jt)

call lsoda(fcri,neq,y t, tout, itol,rtol,atol,jtaskjstate

20 _{for2nat(5ej8.j)}

write (2,20)t,y(1) ,y(2),y(3) ,y(4)

write (3, 20) t,y(1) write(4,20)t,y(3) write(5,20)y(1) ,y(2) Write(6,20)y(3) ,y(4) ethl=.075*cos(y(1)) eth3=.o75*COS((3)) eth2=ethl+eth3 write(7,20) t, ethl write(8,20) t, eth2 write(9,20)t, eth3 write(1O,20)ethjy(2) write(lj., 20)eth3,y(4) ntiines, delt (V(i), i1,4)

morxé,mtwo, length, gray, amp,w,xino

40 _{tout = tout+delt}

stop

end

C _*

subroutine _{fcn(neq, t,y, f)}

C

c _{this subroutine}

defines the _{state space}

representation c _{for the}

equations of

motion of the modeled _system

C implicit rea].*8 (a-h, o-z) C real*8 mone,mtwo,length,grav,amp,w,xinot

reàl*8 kone,_{ktwó, kthree,}

kfour, cone,

ctwo,cthree, cfour

real*8

m2,x2bar,ylbar,y2bar,_{gbar, fbar}

real*8

k1,k2,k3,k4,ci,c2,c3c4 dimension

y(4),f(4),g(4)

common /great/

mone, mtwo, length,

gray, amp, w, xinot common /danp/

kone, ktwo, kthree, kfour,

cone, ctwo, cthree,_cfour

C m2 = mtwo/ (mone+flitwo) c2bar=2.0 ylbar= (2*mone*grav) / (kone*length)

y2bar= (2*mtwo*grav)_I(ktwot

length) gbar= 2*grav/length fbar= 2*amp/length k1 konef(mone+mtwo) k2= ktwo/ (mone+rntwo) k3= kthree/ (mone+mtwo) k4 _{4*kfou.r/ (length*} (mone+mtwo)) cl= cone/_(mone+mtwo,) c2= ctwo/ (mone+mtwo) c3= cthree/_(ne+mtwo) c4= 4*cfôur/_(].ength* (mone+mtwo)) C C C sy3yl=dsin(y(3)-y(l)) cy3yl=dcos(y(3)-y(].)) syl=dsirt(y(1)) sy3=dsin(y(3)) cyl=dcos(y(1)) cy3=dcos(y(3)) force=_fbar*w*w*sin(w*t) denoml=1-in2 *cy3y1*cy3yl denoni2=m2 -m2* _{tcy3yl *cy3yl}

*sy3yl*(y(2)*y(2)*cy3yl+y(4)*y(4))

h2= force*syl__{force*m2 *cy3y1} tsy3

h3=gbar* syl _gbar *m2 *Cy3yj *y3

h4=kl*syl*

(-cy1+ylbar)

hS=k2 * (-cy].

-cy3+y2bar) * (y3_*cy3y1+sy1)

h6=k3 *

(syl+sy3-x2 bar) *_{(y3 cy3yl+cyl)}

h7k4*

(y(1)-y(3)4xinot) *

(cy3yl+1)

h8=cl*y(2)_*syl*syl

C

hlO=c3(y(2) *cyl+y(4) *cy3) *_(j3*y3y1+,j)

h11c4*(_y(4)+y(2) )(cy3yi+i)

pl=m2*sy3yl,* (y(4) *y(4) *1fl2*cy3y1+y(2) *y(2))

p2force*xn2* (syltcy3yl_sy3)

p3g**

_{(syl*cy3y1_y3)}

p4=kl*m2*syl* (-cyl+ylbar) *cy3yl

p5=k2* (-cy1-cy3+y2) * (*Cy3y1*Sy1....Sy3) p6=k3* (syl+sy3-x21,ar) * p8=cl*y(2) *syl*sy**y3 p9=c2*(y(2) PlOC3*(y(2)*cy1+y(4)*cy3)*(m2*cy3y1*l_Cy3) P1lc4*(y(2)_y(4))*(m2*cy3yl+1) ht=hl+h2+h3 -h4-hh6_h7_h8...h9_h1O...h11 Pt=-pi-p2-p3 +P4+PS+p6+p7+pB+p9+plO+pll g(2)=ht/deno g (4) Pt/denom2 g(1) = O.dO g(3) O.dO = y(2) = g(2) = y(4) = g(4) return end

CHAPTIR lv

CALCULATIONS OF PARAMFTERS

Now that the equations of_{motion have been derived and the computer simulation code has} been written, numerical values_{for the constant Variables of the system need to be determined}

so as to match the characteristics of the Faraday wave. These constants are the_{masses, dashpots,} springs and length of the rods. Data from the Faraday _{eperiment [1] is used tocompute the} cOn-stants. information about the tank size and water volume is used_{to determine the values for the} mpsses and the length of the rods. Linear Faraday wave data_{axe used to determine the values for} the spring and damping coefficients.

The determination of these values_{can be very difficult in a non-linear problem such as this} one. Simpliflcations must be used to create a physical model whose_{constants can be easily} deter-mined. In the Faraday experiment [I], to measure damping, the_{tank was forced at a resonate} fre-quency to generate a standing wave. The excitation was stopped, and the_{standing wave was}

allowed top out completely.

_{Wave amplitude data}

were collected to compute the damping

effects of the Faraday_{wave. Reference [1] has a graph of the}

amplitudes of the wave on page 286,

figure 4. The peak amplitudes show a decay with a calculated decay rate of 0.050 1/s. This_value is used to determine the damping_{and spring constants.}

It is noted that the decay rate was determined when the amplitude was small Thiscan be

considered as part of the linearregion of the system. Therefore, _{to simplify the calculations of the} constants, it is beneficial to linearize the current model and then use the_{experimental data to} determine the values of the constants. The equations 6f motion for_{the linear model must now be} derived. The linear model will be derived in two ways. First, a linear model will be used_{and the} equations of motion will be determined from this model. Second,_{to serve as a check the} n-lin-ear equations of motion will be linn-lin-earized to find the inn-lin-ear equations of_motion.

4.1 _{Mass and rod length}

Before deriving_{the linear and}

linearized equations_{of motion, the}

values for each point mass and the values of the niassless_{rods that}

connects each point mass are determined. _{Figure 6} shows a sketch of the tank, its_{dimensions, and the}

model.

Figure 6

Sketch of the_{difliensions of tank} Geometry is used_{to help find the values of thC point}

masses and the lengths_{of the rods. The} length of the tank is 600 mm. The_{three degree of}

freedom system goes from node to node, hence it is equal to half the length of the tank fora linear

wave. The length of the two degree of freedom model is equal to half the length of_{the three degree}

of freedom model. Therefore,_{the System} being modeled is 150 mm long and each_{massless rod is 75}

mm long.

For the masses, let each one_{represents a certain section}

of the tank. If_{the volume ofeach}

section is the same, then the weight for_{each mass is the}

same. When the model _{Was simplifiedto}

a tWOd

_{-of.frjom system,}

the system_{mass was bayed. This}

means that the mass, M2,_is

also half its original mass. For_{simplicity, the total}

unity. This is acceptable because the system was normalized with respect to the total mass. This meansthatM1 will be 2/3 of the total _{mass and M2 will be 1/3 of the totalmass Actual values for}

the masses can be found by_{using the density of water, but this is unnecessary because the}

system

is normalized with respect to the total mass Only proper ratios_{are necessary.}

4.2 _{Derivation of equations of motion using linear theory}

The first derivation of the linear_{equations of mOtion is from a lineat model of the system.} There are a few assumptions that are used in liflear theory that will be discussed_{1ater To start, it} is assumed that the model is lying horizontally Two new angles are used to represent the_two

degrees of freedom.

Figure 7

Sketch of model showing coordinate transfonnation

Figure 7 shows a simplified model of the system with the two new angles. For the linear derivation, the equations of motion are derived from these angles. Using linear theory, thesmall

angle approximation is used. This approxixnadon assumes that the angles O and 02.are_small. This allows the following simplifications to the trigonometric values of the angles

Again, two equations of motion

are needed to describe the two degree of_{freedom system.} Two_{moments are used}

as before to find the equations of motion:_{the moment about}

the pin con-nection, R, and the mOmer about M1. In_{the linear model, c3 and k3 do not have}

an effect on the

system fOr small_angles

ofpeturbatio _{therefore, they}

can be ignoredj It will be seen, later in

this paper, that k4 and c4 are not needetj in this_{model and are not used in the linear}

model. c4 and are retained to allow for their_{use in future research}

projects. First, the moment about the pin

Connection, R, is found to get the following:

MR = 0 = kI6I+k2(oe)cll2

Rearranging this _{equation gives:}

12i..

2,

1121r?

{mi+rn2_]O1

_{+{rn2Je2+[c_.cjo}

+{c2jO2+[kl_+k...jo

+ Next the moment about the first_{mass, M1, is found}

to get the f011owing.

M1=-

rn1

(rn1 +rn2) k1

+rn2)

m(61

+e2)

= 0 (17)

1Mm

_{+ 62) +}

_c(e

_{+ O2) + m2(O1} Again, rearranging this equation_produces:

[rn2çjo1

_{+[rn2çJ2+{c2çje1}

+IIc2c]92+fk2cjel+{k2]e2

= 0 (18)

Equations (17) and_{(18) are the equations of motion for the}

linear system. To simplify the system and to match the equations of motion_{to what were derived}

from before, the_{following termsare} defined;

k2

+ m2)

C2_(_+)

-Substitute these definitions_{into equations (17)}

and (18), and multiply each by _{to get:}

12(m1 +m2)

[M1 + 2M2Je1 + [2M2]O2 + ₊

+ [22]O2 +[R +2A2jo1 + [2k2]02 = 0

(19)

.+ [M2J2 + [C2J01

_{+ [2]}

_{+ [R2]81 + [A2]02 = 0}

(20) Next, rewrite_{equation (19) to get the following form.}

[M1 + M]O1 + [M2]e2 + [C1 + c2191 + [c2je2+ [K1 +_{K2J0, + [K2]O}

= [[M2101 +[M2]02 + [c2]e1 _{+ [C2 102 + [K2]01 +}

[K2J02]

Observe that the right_{hand side ofequation}

(21) is equal to the left_{hand side of equation (20)} which is equal to_{zero. This means that the}

right hand side of equation (21) is equalto zero. Also

note that M1 is twice the_{mass of M2. Recall the}

mass ratios discussed from before where M1 is

equal to 2/3 and M2 is equal to 1/3. Therefore,

the total mass is equal to one or M1 +M2 ₌

Equations (20) and (21) can be written as:

+ [M]02 + E:c1 + c2je1 _{+[c2102 + (K1 + K2J01}

+ [2]02 = 0

(22)

[M2j1 + [M2]02 + [2JO1 + [2]O2 + [K2]O,

_{+ [Rj02 = 0}

(23)

Equations (22) and (23) are the final equations of_{motion for the linear modeL}

4.3 _{Linearization of}

flon-Jine _{equatio of motion}

Now the _{non-linear equations}

of motion, that_{were derived in chapter 2,}

will be linearjed

to derive the linear_{eqUations of motion.}

This will verify_{that the linear equations that were} derived using_{a linear model are correct.}

Refer back_{to Figure 7. It}

can be Seen that if the model is to be horizontal,_{the angles 0}

and must have the value

_of

. Therefore, the angles 0 and _{canbe converted}

to the angles

_{0 and}

e2 by using the

foUowingdefinjtj5

Next, these_{definitions are used}

to make a

_{transfomajo of the}

tiigonoznetrjc functions from the original non-linear equalion of motion_{to angles in the}

linearized equations Of motion in the_new

coordinate system. The _{definition for the sine}

and cosine of 0 and_{$ are converted to angles}

of 0

and 02, and are linearized using

_{a small angle approxij0}

The linearized definitions of the

sin(0) and the cos(0) are as follows:

sin(_O1)

cos(01)

= 1

_{cos(_o1)}

= sin(01) = 8,

On a side note, the following_{trigonometric identities}

ate necessazy to complete the_conversions.

sin(xy)

_{= sinxc,sy... cosxsiny}

Cos(xy)

_{cosxcosy+ sinxsiny}

The final

conversion definitions for sin(0) and cos(0)_{are found using the above}

trigonometric identities:

.II,

,$i IC.

m_01)

₌

= cos91

cos(_91)

_{= coscos01}

₊

_si4sine,

Therefore:

sin01

cOs001

The same logic is used _{tofind the final conversion definitions}

for the sin($) and the cos(0).

sm4)l

cos$=02

Looking at the non-linear equations of motion, note that_{there are terms of sm(-0) and cos(40)} in the equations. To convert these to the linearized form _{in the flew coordinate system, the above} trigonometric identities are used. Then, using the srtiall_{angle .approxim2tjon the sine and cosine}

terms are converted.

s'n(-0)

=

_{sin$cOsOcos$sjn0}

_{= 01+02}

cos($ 0)

=

_{coscosO + sinsin0}

_{= °I°2}

_{+ 1 =}

1aOTo

The cos($.0)_{conversion creates a term that has}

the product of two angles. This value is a higher order term (H.O.T) _{and is considered srnall}

enough to ignore. Next the conversions for the deny-athres of the angle 0_{need to be determined and}

expressed in terms of the fltt _{and secOnd} deriva-tive of 01.

Notice that the first_{and second derivative of 0 is}

equal to the negative ofthe derivatives of 0.

The. same derivative

conversions are Computed for theangle$.

$ =

Now that all of the_{conversion de nitions have been derived, all of these}

conversions.

be used to linearize_{the non-linear equations}

of motion into the new oordinate system. _Referring to equation (13) again from the original non4inear equations:

o + 4M2cos(q. 8)-2M2sin($_

O)f(t)sjnO_sm9

_K!sine(cose_yI)+;2sme(cosO

_{= COS 'Y2)}

+K3coso(sjo+

* c2slno(esine

_{+ $sin) + 3coS9(co}

+ cos) + C4(= 8)

_{= 0}

Substitute the

Conversion lefinitions derived above into_{thi equation:}

(24)

Recall, in chapter 2, some of_{the defiñitio that were made}

during the original derivation of the quaxions of motion.

Simp1iflcation.o the linearized equation of_{motion can be made using these} definitions. Note the_following:

x2

Notice that; is in the 1inearjj

_{equation above as part of the A'}

term. With x2 _{equalto2, the} K3 term is then equal to zero.

= 0

Note also that (O2)2, (O1)2

_{(OiO) and}

(2O2) are áIlhigherorder _{When looking}

_atjust

the linear case, these higher_{order terms do}

not contribute Significantly_{to the equations and can}

therefore be ignorel_{as was stated above.}

Remember that _{5 was defined as follows:}

For this system,

_y

Now for the _{and A2 terms, recall:}

Usingthese_definitions,

and .v2 can be multiplied by and

A.

2m1g

1

l(rn1+m2) K2y2=

l(m1 +m2)

Adding these two prOducts together gives:

-

2(m1+rn2)g K1y1+K2y2=

l(m1+m2) 1

Also recall = _{Using all the above definitions}

and substituting them into equatiOn (24)

gives:

01

_{+M202+(clc2)9Jc?82(KK)e+K$}

= 0

₍₂₅₎

Equation (25) is flu_{al linearized equation of}

motion for the first _{non-linear equation. Now reca1i}

equation (14) which_{was the second original non-linear equation of motion.}

yI Likewise

-

2m2g Y2 - IL

"2

2m1g 1k1

OM2cOs(s_e)+sM2+e2Msn(.e)y(t)M,

sin(_ cosO

cos4

Y2)+ K3cosO(srne +

+C2s1n$(6sine+qsin4)+cos4,(Oe._

_o

First substitute _{the conversion}

for the angle _factors.

81M2--02M2+(_eI)2M2(91

_{+O2)M2+2(...O.}

02Y2)+K392(1 + I

X)

= o

Note that 2gM2

_-2m2g

jM2 =

_{= l(m1+m2)}

= K2y2

Using this definition and the_{definitions from above,}

the final linearized_{equation of motion}

can be calculated for the second non-iinear equation:

M2OI+M2e2+C28I+.9+KOKe

= o

₍₂₆₎

Now it can be seen that the linearized

equations of motion and the linear_{equations of motionare} exactly the same. With these

equations, values for the spring and damping_{coefficients can be} derived to match_{the ch ctensties of the Faraday}

waves. 4.4 _{Calculations of constants}

Now that_{a linear model of the}

system has been derived,

values for the constaflt variables

can be found that_{match the characterjstic}

of the Faraday wave. A modal_{atialysis of the linear} equations of motion_{will allow for the}

determination of the_{mode shapes that will give}

a linear response. This linear response will allow for_{a comparison of the}

response of the computer

sixnu-lation to the response of the experiment _{data at small}

angles. This_{companson will verify if the}

modeled system has the same characteristi

as the actual Faraday wave. Recall the linear equa-lions of motion for the thodeled_{system equations (25) and (26).}

01 +M20 +(C1+

C2)O1 + C202+(K1 +K2)91 + K202 = 0

M201+M202+c201+c282+K291+K202

0

These equations can be rewri ten to take thefollowing form:

M22 + C202 + K282 -(11201 + C20i + 28i1, (27) and

M202 + C202 + K262 = [Gi + (C1 + C2)0i +(K1 + K2)01] (28) Notice that the right hand side of equations (21) and (28) are equal toeach Other Setting the left

hand side Of eachequation equal to each othergives:

O1+(C1+C2)0i+(Kl+t2)O1 =

M201+C201+K201

After some algebraicmanipulations the following tWOequations of motion are ascertained:

.u_M2)ei+clel+KIel

0 . (29)

and

M20I+M202+C20l+C202tK2O12O2 =

0

(30)

In matrix form, these arewritten;

[1 M

011011 +

1i

o

[o1 + fR

011011 = rol L M2 M2J 1.oI L2 C2J Le2] [2 K2] L02i L0J

Next, rewriting eqIatiOn (29) to get

-1

= (1_M2)61(1_M2) I

Equation (31) is then substituted into equation (30):

1-I

2L(l_M2)01_(l_M)01]+M202+C201+C202+K201+K202 = 0

(32)

Rearranging equation (32) gives:

M2O2₊

(2

(1 M))01 + C202+

(2

-

(iM2))°' + K202 = 0

(33) Eqüations(29) and (33) are used to find the equations of motion in matrix form.

11M2 o111

0 1e11

0

Fe11 lol

[

0

Mj

Lo2i

c2(-(

C Le2i

+ K2(l)

K2 Loj = Loi

From the previous sections, it is known that the masses, spring coefficients, and

danipingcoeffi-cientsmaintaincertainratiostoeachother. Thisprovidesthat M1

,M2 =

andthat

2R2 =

C1 = 22. Substituting these values into the above matrix gives:

o11oI1+12X2 01181 _rol

o

LJ

0 c2j 1e21

[ 0 K2j

L02]

iq

Multiplying the first row by and multiplying the second row by3 gives:

F101111+32 0

Fei LO 1J[o2j

0 32 [e2

i32 ol[ei

L 0 3R2j [02 -.101 LqI -4' This produces two decoupled, single degree-offreedoni equations of motion that have the form:

01+3C201+3K201 = 0 - (34)

02+3C202+3K292 0 (35)

This shows that the masses are uncoupled. In other words, in the linear case, the position and velocity of one mass does not, have an effect on the position of the second mass. One of the m2in reasons is that in the linear case, K3 and C3 have no effe t on the system. Without these, the total

system is juM tWo uncoupled systems.

Data that was derived, in the original Faraday experiment [1], is used to calculate the

-spring and damping values. The forcing frequency of oscillation for the period one Wave was 3.23

H2 [1] The frequency of the period one wave was -at half this 'alue. Therefore, the frequency was 10.147 rad/sec. This is the natural damped frequency of the model, Wd. It is also known that

the. ihmping tate of the faraday wave was 0.O5 [1]. This informati- on is used to compute a value

for the spring and damping coefficients

Under standard vibration theory s the damping tate of a one degree of freedom system usually defined as [7]:

and the damping ratio is:

r

2Ma) Therefore

, C 3 2

Because the damping rate beta is equal_{to 0.05 the damping values}

0 the _{system is calculated to} betobe, C2

_{= 0.o. Because}

istwicethevalueof2

_{= 0.0g.}

The damped_{frequency for a single degree of freedom}

spring-dashpot system is [6]:

=

Jj2

Using the definitions

_{of and o, the}

equation for the damped_{frequency can be modified}

to be:

£

Knowing that _{0d= 10.147 and = 0.05, simple}

calculation provide that A;2

_{= 34.323. This}

means that A;1 = 68.64& _{Tis information}

can be used to simulate_{the free vibration of the sys.}

tern.

When simulating_{the free vibration,it is important}

to note the value of the dampingrate. To see if the _{system matches the llneai}

response produced by the_{faraday wave, the computer}

sini-ulation has to be tested using small angles as initial_conditions.

One other_{problem that arises is that}

the system has _{two degrees of freedom}

and not one.

This means that the system has two mode shapes that_{need to be found}

to simulate the linear damping. The mode shapes are found by finding_{the eIgenvecto}

of the linear system.

The mode shapes_{of the system vill be}

calculated by finding the mode Shapes_{of the linear} form of the model.

Recall equations (25) and (26) which are the linear equations_{of motion.} Ignoring the damping coefficients, these equations_{in matrix form are:}

11 MJ Foil

[M2

M2JLJ

For simplification,_{use the following}

nOtation [61:

112211°1l

= fol

I

2 A;2][02] _[.0I

Me+X0

= 0

This can be rewritten to be [6]:

Mke

The linear equations are. now in state space form. From here_{define the A matrix to be [6]:}

=

To find the mode shapes ofthe linear system, flhd the_{eigenvectors of the A matrix [6].}

The mode shapes for this system_are:

['liol

The first mode shape will be sun lated as a checL The initial _{angular velOcities for each}

angleiszero. Thismeansthem

_{esstartfromrest Letthein}

angu1arpositionofo...

+ 0.001 _{1.5717963268. Because the initial}

angular position of the second_{mass is uncoupled, it}

does not matter where the mass is placed. For this test, let _{have an initial angular position}

_of

= _{= 1.5707963268. This is a value that has the}

second mass flat with respect to the horizontal.

Figure 8 shows the results of this simulation. Notice that_{even though mass one is oscillating,}

mass to remains horizontal. This_{shows that the}

movement of mass One does not affect the

movement of mass two. Next, the_{same Simulation will be tried but in}

reverse order. Begin, again, with zero initial angular velocity.This time 0_{will have a hOrizontal initial condition, or 0=}

3

= 1.57079632, and

_{will have an initial position of =}

+ 0.001 = 15717963268.

Figure 9 shows the result of this simulation. Again, the _{Same results as in the first simula} don are shown in the_{second simulation. For}

3

II

--height of meBa I'andmaee

2(metere) ---'. .4... -S ...p _{S...IP St...} I-. Lii -.1 II

jTf'':

1:1

0

o

height of mesa land

rnass2. (meters)

21

I

which both masses are going to oscillate. Again stall frothrest by having both_{angular velocities} zero. Consider initial angular positions

_{of 0=}

- 0.002= 1.5687963268, and $ = + 0.001 =

1.5717963268. Figure 10 shows this _{results of' this simulation.}

0.00015 0.0001 os -0.00015 0

Ii

2 4 $ (sec) A 9 10 Figure 10

Linear response of system with both_{M1 and 142 displaced} 0 = 1.5687963268 rad, _{= 1.5717963268 rad, 0 =$= 0}

With this simulation, it is again verified that the two degrees of freedom _{are indeed uncoupled.} The position of mass one was initially started at twice the_{height of mass two. Throughout the} oscillation, mass one was always double the height of mass two. The masses also_{maintained the} same damped frequency. The movement of one mass did not affect the_{movement of the other} mass.

It needs to be determined whether the response of the _{system matches the time histoty}

of

the Faraday wave. A fortran program_{was written that is able to determine the time and heights of} the peak elevation for each oscillation, With this, the_{damped frequency and the damping rate}

can be calculated. To find the damping frequency, the period of each_{rnaxnnum needs to be} deter-mined and then converted to rad/sec. To filid the damping_{rate, the natural log oEeach maximum} is calculated and plotted, then the least_{square line that best approxirn2u's the}

points is found. The

lIlItIf H III! 1111i1ii

slope of the line wiljb

_{the am,ing rate}

To assist in the visualization Of this, the following_{example was created. The following} formula is used to demonstrate a damped oscillating _{Wave: y(t) = esin(cot). Figure}

11 shows the

time history of the_{damped sine wave using}

a natural frequency of o = 10 rad/sec and a damping

rate of 0.3/sec.

-i

Figure 11 _{Figure 12}

Plot of damped sine wave with envelope curve _{Plot of the wave-peak decay}

with least-square

w lOandP=O.3/sec

approximation 5=O3Isec The wave is damped_{by the e terni and forms the envelope}

curve. is the damping_{rate that}

was discussed previously in this section. From this, the peak_{values that form the envelope curve} are found. Two of these peak values can be used to find the_damp

frequency. The damped fre-quency is found by:

COd(.)(

Remember that the damped_{frequency is not the Same}

as the natural frequency. From before, the damped frequency is_{equal to}

0d=

With a damping_{rate of 03 arni a natural frequency of}

10 rad/sec, the damped_{frequency should be} about9.995 _{rad/sec. Using the data from the computer simulation,}

calculated to be 9973 rad/sec. This is very close to the atual _value.

To calculate the damping_{rate of the Simulation, begin by looking at the part of the}

equa-tion that has. the damping_{rate by iseif. The equation isas follows:}

y(t) =: Ae

Taking the naturaj logarithm of each side and doing_{some algebraic manipulation, gives the} fol-lowing:

ln(y) = ln(A)z

Notice that this fonxaila has._{the same format asa for ula for a line in Which}

is the slope.

Divide each peak value by_{the value of the very first peak.and then take the iaturai logarithm} of it.

A.=

(H

Figure 12 show a plot of all these points. Using _{a least squares approxinition, the slope} Of the line can. be found and this_{will give the value of the damping}

rate. For the example problem, the damping rate is 0.300844, which is very close to the original damping rate of 03.

figure 13

Linear response of M1 with envelope curve

0=1.5717963268 md, _{=1.5707963268 rad, 0=41 = 0}

figure 14

Plot of the wave-peak decay with least-square

approximation = 0.05/sec

Now look at the_{damping rate for the system that}

was originally modeled. The flstmode

shape is used and only_{the response of the first}

mass isplottecL Figures 13 and 14 show_the response of the system and the least_{squares approximation. For this Simulation,}

the damped fre-quency was calculated to be 10._{150 rad/sec which is very close to the Faraday waves 10. 147 rad/} sec. The damping rate_{was found to be O.OS00008/sec which}

is the same as the Faraday_waves

O.OSIsec. This shows that_{the computer modelhas the same linear characteristics as the linear} rsponse of the Faraday wave.

Next, the system will be _{forced at the proper frequency}

and amplitude to get a period one

CHAPTERS V

NONLINEAR SIMULATION RESULTS

Now that a computer simulation of the model that matches the characteriStics of the Fara-day wave exists for the linear case, tests need to be performed to ensure that the model matches the character tics for forced amplification. Simulations need to be completed to find period one and period three waves. The period one and penod three waves are found through trial and error by running numerous Simulations and changing initial, conditions., Te initial conditions that are

changed are the angular pOsition and angular Velocity of the two degrees of freedom 0 and _. The forcing frequency and forcing amplitUde is also changed. After this, 'it is important to see if the

period one and period three wave can be attained by starting the system from rest in the horizontal

position and then using oüly the forcing frequency to excite the system to a period one and period

three wave.

There are numerous initial conditions that can be used to produce a period one wave. Some of the waves initiate as period one and then become unstable overtime. Further research would be to find all the stable and unstable period one waves and determine areas of stability.

To demonstrate the period one wave, the following initial conditions are chosen. For this example, the initial condition ofOis 0 = 1.9394222 rad.. This is equivalent to starting at 0.027025 meters below the horizontal. The initial condition of is Q 0.75032909 rad which is 0.05486 meters above the horizontal. The initial angular velocities are 7.665593 rad/sec for 0 and 7.3153634 rad/sec for 4,. The forcing function has an amplitude of 0.006 meters and a frequency of 10.2324844 tad/sec. Figure 15 exhibits the time history of the response Of the System. The first plot is the response of just mass one. This is the height of the point mass M1 relative to the horizontaL The second line represents the height of the.second point mass M2 relative toM1. The

third response is a time history of the _{total height of the system relativeto the horizontal. When}

a

standing wave is being simulated,_{as in this case, the total height is the maximum amplitude} of the system.

-Figure 15

System response with period one standingwave

Figure 16

Phase plane of system: period_one

Figure 16 is the phase plane plot of the first_{point mass M1. This shows that the}

system is a stable period one standing wave.

Obviously, trying to start a system like this in a laboratory environment using these initial conditions is very unrealistic. In the original Faraday_{experiments, the tank was started from}

rest. Over time, the water in the tank evolved to a standing wave of period one. For the computer model, it is desired to match the results of the experiment In the experiment, an amplitude of 0.00265 meters and a frequency of 20.106 rad/sec [2]. This gave a damped frequency of 10.053 rad/sec (the Faraday wave is_{a first subharmonic, which means that the forcing frequency}

is twice the frequency of the Faraday wave). The model in this thesis is not garnered in the same manner as the Faraday wave. The forcing frequency for all of the computer simulations will be the same as the damped frequency from the experiments. For this example the frequency used is 10.053 rad/sec. The initial angular velocity for_{both B and is set to zero and the system isstarted in the} horizontal position, which is _{for both angles. Figure 17 shows the time history of the} simula-tion for these condisimula-tions.

53 t 1.1

E 0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 0 20 40 60 80 100 me (see) Figure 17

System response, from rest, of periodone standing waves amp =0.00265metezs, w = 10.053 rad/sec

Note that after about 150 seconds the system is in period one motion. To better analyze the sys-tem, figure 18 shows the response of the M1 component of the system from 190 seconds to 200

seconds.

Figure 18

System response, from rest, of periodone

standing wave from 190 sec to 200 sec amp = 0.00265 meters, w = 10.053 red/sec

Figure 19

Phase plane of system, from rest from l9Osect000secpeiiodcee amp = 0.00265 meters, 0=10.053 md/sec

I'S IN * , 120 140 160 180

-I

From figure 18, it_{appears that the system is exhibiting periodic motion and specifically penod}

one motion. To check this, the phase of the_{system is plotted using only data points from 190} sec-onds to 200 secsec-onds. Figure 19 shows_{the phase plane and it is seen that thesystem is period one.} These simulations demonstrate that the computermodel can be started from rest and resonated to a period one standing wave using the same amplitude and frequency as in the physical Faraday

wave experiments.

One of the period_{one responses of the Faraday wave in the original experiment had a} dim-ple in the free surface _{[1] If this dimpled crest can be simulated, it will provide cv dence that the}

model matches some features of the Fataday wave. The dimple effect.is Illustrated using Figure 20. The figure shows that as the wave is at its max mum elevation, a dimpled crest forms at the peak of the wave [1]. In the experiment,_{the dimpled wave was found usinga forcing amplitude} of 0.00385 meters. Figure 21 shows the response of the system using a forcing amplitude of

0.003705 meters and the same forcing frequency, 10.053 rad/sec.

Figure 20

Sketch of dimpled standing_{wave with idealized model}

a

jut:

V I I I V

ti

t q

___=

a I Figure 21

System response of periOd one dimple amp = 0.003705 meters,

0

_{= 10.053 rad/sec}

Figure 21 is the corresponding system response andit shOws that the system exhibits the dimpled response. The model shows that the dimple _{occurs at both the maximupi and minimum elevations}

unlike the Faraday Wave experiment.

To further test the model, it is desirable to match the period three response of the experi-mental system. Again, trial and error was used to find initial conditions that have a period three

$ _I ¶

response. The initial position of 0 was 0 = 1.9817702 rad and fore was$ =0.66793155 rad. These are equivalent to placing M1 at 0.02996267_{meters below the horizontal and placing M2 at} 0.0588828 meters above the horizontal. The initial angular velocities are 0 = 7.6345674 rad/sec

and _{4) = 7.2529606 rad/sec. The forcing function has an amplitude of 0.006meters at a frequency} of 10.18336536 rad/sec. Figure 22 shows the_{response of the system using these initial} condi-tions. Figure 23 is the phase plane of the system forM1.

2 $ Figure 22

-at -to S Figure 23

System response with peiiod three standing wave _{Phase plane of system: penod three}

From the phase plane, it can be seen that the _{system is a period three standing wave. Again, these} initial conditions are not easily implemented in a laboratory setting. The system should be tested from rest, using a forcing frequency and amplitude,_{and excited to a period three wave. This}

proved to be a difficult task.

As with the period one standing wave, the initial_{angular velocity is zero for both 0 and} . The initial angular positions for e and are _{, or horizontal. In the original Faraday wave} exper-iments, the period three standing wave was obtained using a forcing frequency of 10.053 Hz_and an amplitude of 0.0046 meters [21.

Using these values in the simulation did_{not generate a period three standingwave. The} model showed this wave to be breaking. The_{explanation of the simulated breaking wave will be} explained shortly. To find the period three standing_{wave, the value of the amplitude of the}

fore

ing was varied,_{while the forcing frequency was held constant. The amplitude was changed until} period three was found. At an amplitude of 0.0038meters, a period three standing wave was

found using the computer simulation. Figure 24 shows the time histoiy of the response.

E

Figure 24

System response, from rest, of period three standing wave, amp= 0.0038 meters,₀= 10.053 rad/sec

.0 S

-Figure 25 _{FIgure 26}

System response, from rest, of period three Phase plane of system, from rest,

standing wave from 240 sec to 250 sec _{from 240 sec to 250 sec: period three} amp = 0.0038 meters, o = 10.053 rad/sec _{amp =0.0038 meters, 0=10.053 radFsec}

¶

I

I

3 3 3 S

To facilitate viewing the period three motion, Figure 25 is a graph of the system from 240 seconds

to 250 seconds. Figure 26 shows the phase plane of the system with these sair data. It is clear that the system is exhibiting period three behavior. Note that this is a very stable period three

standing wave.

Figure 27 Figure 28

Plot of system demonstrating breaking waves Plc of system demonstrating breaking waves for height of M1: amp = 0.0045 meters, for angle of 0: amp = 0.0045 meters,

W= 10.OS3rad/sec (U= 10.053 rad/sec

IO 0 10 20 30 40 70

Figure 29

Increasing the forcing amplitude to 0.0045_{meters was done next to demonstrate chaotic} behavior. Figure 27 shows the results of the smulation. Itis difficult to see whether the system is

exhibiting breaking wave behavior in Figure 27. To facilitate_{this, Figure 28 is a plot of the} angu-lar position Of 0 for the same simulation. It is clear that the mo el is rotating about the pin con-nect.ion. For this model, this behavior is considered wave breaking. A plot of the phase is shown in Figure 29. This shows that the system does not become periodic. The system shows contiñu-ous breaking wave behavior. It is easily seen that increasing the amplitude of the forcing gener-ates more breaking waves and chaotic-like behavior of the system. This then shows considerable similarity with the results of the Faraday experiment in that the_{system goes from period one, to} periOd three, to bteaking.