Vibration

VIBRATION SIMULATION USING MATLAB

By Park, Jeong Gyu

DEPARTMENT OF PRECISION ENGINEERING KYOTO UNIVERSITY KYOTO, JAPAN JAPAN MAY 2003

c Copyright by Park, Jeong Gyu, 2003

Table of Contents Table of Contents

ii

1 Basics of Matlab

1

1.1

Making matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Matrix Manipulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.4

Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.5

Programming in MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.5.1

The m-files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.5.2

Repeating with for loops . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.5.3

If statements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.5.4

Writing function subroutines . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.6

Saving and Loading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.7

Help . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2 Single Degree of Freedom System 2.1 2.1

2.2

2.3

7

Free ree Vibr Vibrat atio ions ns of Sing Single le-D -Deg egre reee-of of-F -Fre reed edom om Syst System emss . . . . . . . . . . . . . . . . . .

7

2.1.1

Viscous Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

Forced Vibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2.1

Direct Force Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2.2

Base Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

Simulation with MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.3.1

12

Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ii

Table of Contents Table of Contents

ii

1 Basics of Matlab

1

1.1

Making matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Matrix Manipulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.4

Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.5

Programming in MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.5.1

The m-files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.5.2

Repeating with for loops . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.5.3

If statements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.5.4

Writing function subroutines . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.6

Saving and Loading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.7

Help . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2 Single Degree of Freedom System 2.1 2.1

2.2

2.3

7

Free ree Vibr Vibrat atio ions ns of Sing Single le-D -Deg egre reee-of of-F -Fre reed edom om Syst System emss . . . . . . . . . . . . . . . . . .

7

2.1.1

Viscous Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

Forced Vibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2.1

Direct Force Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2.2

Base Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

Simulation with MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.3.1

12

Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ii

2.3.2

State Space Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 Multiple Degree of Freedom Systems 3.1 3.1

13 17

Some Some Bas Basics ics Conce oncept ptss fo for Line Linear ar Vibr Vibrat atin ingg Syst Systeem . . . . . . . . . . . . . . . . . . .

17

3.1.1

Eigenvalue Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.1.2

Orthogonality of normal mode odes . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.1.3

Normalization of Mode ode Shapes pes . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.1.4

Modal Coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

3.2

Proportional Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

3.3

Moda odal Analysis of the Force Respon ponse . . . . . . . . . . . . . . . . . . . . . . . . . .

25

3.4

State-Space Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.4.1

Free Vibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.4.2

Forced Vibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

4 Design for Vibration Suppression

31

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

4.2

Vibration Absorber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

4.2.1

SDOF with Undamped DVA . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

4.2.2

SDOF with damped DVA . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

Isolation Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

4.3.1

Passive Isolators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

4.3.2

Skyhook and Active Isolators . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

4.3.3

Semi-active Isolators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

4.3

5 Vibration of strings and ro ds

42

5.1

Strings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

5.2

Rods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

6 Bending of Beam

44

6.1

Equation of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

6.2

Eigenvalue Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

6.2.1

Boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

Some Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

6.3

iii

6.4

Forced Vibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

6.4.1

Point force excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

6.4.2

Moment excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

7 Plate

54

7.1

Plate in Bending . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

7.2

Equation of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

8 Approximate Method

57

8.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

8.2

Rayleigh Ritz Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

9 Finite Element Analysis 9.1

69

Euler-Bernoulli Beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

9.1.1

Basic relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

9.1.2

Finite Element Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

Thin Plate Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

9.2.1

formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

9.3

Finite Element Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

73

9.4

Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

9.2

iv

Chapter 1

Basics of Matlab Matlab R (MathWorks, Inc., )1 is an interactive program for numerical computation and data visualization. It is used extensively by vibration and control engineers for analysis and design. There are many different toolboxes available which extend the basic functions of MATLAB into different areas.



1.1

Making matrix

Matlab uses variables that are defined to be matrices. A matrix is a collection of numerical values that are organized into a specific configuration of rows and columns. Here are examples of matrices that could be defined in Matlab. A = [1 2 3 4;5 6 7 8;9 10 11 12] Transpose of a matrix using the apostrophe B=A’ C=[2,2,3 4,4,6 5,5,8] The colon operation ’:’ is understood by Matlab to perform special and useful operations. If two integer numbers are separated by a colon, Matlab will generate all of the integers between these two integers. a = 1:8 generates the row vector, a=[12345678] 1

see, http://www.mathworks.com/

1

2

If three numbers, integer or non-integer, are separated by two colons, the middle number is interpreted to be a ”range” and the first and third are interpreted to be ”limits”. Thus b = 0.0 : .2 : 1.0 generates the row vector b = [ 0.0 .2 .4 .6 .8 1.0 ] C=linspace(0,10,21) D=logspace(-1,1,10) eye(3) zeros(3,2)

1.2

Matrix Manipulations

Element of matrix A(2,3) Size size(A) length(a) Transpose A’ Column or row components A(:,3) Matrix addition, subtraction and multiplication D=B*C D=C*B If you have a square matrix, like E, you can also multiply it by itself as many times as you like by raising it to a given power. E  3

∧

Element addition, subtraction and multiplication Another option for matrix manipulation is that you can multiply the corresponding elements of two matrices using the .* operator (the matrices must be the same size to do this). E = [1 2;3 4] F = [2 3;4 5] G = E.*F If wanted to cube each element in the matrix, just use the element-by-element cubing. E.

∧3

3

1.3

Functions

Matlab includes many standard functions. In Matlab sin and pi denotes the trigonometric function sine and the constant π. fun=sin(pi/4) To determine the usage of any function, type help function-name [Example] Verify the variables i, j, cos, exp,log, log10 in MATLAB

1.4

Plotting

One of Matlab most powerful features is the ability to create graphic plots. Here we introduce the elementary ideas for simply presenting a graphic plot of two vectors. Example Plot the sin(x)/x in the interval [π/100, 10π] ————————————– >>x=pi/100:pi/100:10*pi >>y=sin(x)./x >>plot(x,y) >>grid ————————————-

1.5 1.5.1

Programming in MATLAB The m-files

It is convenient to write a number of lines of Matlab code b efore executing the commands. Files that contain a Matlab code are called the m-files.

Table 1.1: Basic Matrix Functions Symbol inv det trace

Explanations Inverse of a matrix Determinant of a matrix Summation of diagonal elements of a matrix

4

Figure 1.1: Sin(x)/x

Table 1.2: Basic Plotting Command Command plot plot(x, (x,y) y) subplot loglog semilogx(x,y) semilogy(x,y) titl titlee xlabel ylabel grid grid

Explanations A Car Carte tesi sian an plot plot of the the vec vecto tors rs x aand nd y A plot of log(x) vs log(y) A plot of log(x) vs y A plot of x vs log(y) plac placiing a titl titlee at top top of gra graphic phicss plot plot

Crea reating ting a grid rid on the the grap graphi hics cs plot lot

5

1.5. 1.5.2 2

Repea Repeati ting ng wit with h for for loops loops

• the for loops

Syntax of the for loop is shown below for n=0:10 x(n+1)=sin(pi*n/10) end The for loops can be nested H=zeros(5) for k=1:5 for l=1:5 H(k,l)=1/(k+l-1) end end

1.5. 1.5.3 3

If stat statem emen ents ts

If statements use relational or logical operations to determine what steps to perform in the solution of a problem.

• the general form of the simple if statement

is

if expression commands end In the case of a simple if statement, if the logical expression is true, the commands is executed. Howeve However, r, if the logical expression expression is false, false, the command is bypassed bypassed and the program control control jumps to the statement that follows the end statement. statement • The if-else statement The if-else statement allows one to execute one set of statements if logical expression is true and a different set of statements if the logical statement is false. The general form is if expression commands(evalu commands(evaluated ated if expression expression is true) else commands(evalu commands(evaluated ated if expression expression is false) end

6

1.5.4

Writing riting function function subroutin subroutines es

function [mean,stdev] = stat(x) n = length(x); mean = sum(x) / n;

1.6 1.6

Savi Saving ng and and Loa Loadi ding ng

who. All variables in the workspace can be viewed by command whos or who. To save all variables from the workspace in binary MAT-file save FILENAME FILENAME An ASCII file is a file containing characters in ASCII format, a format that is independent of ’matlab’ or any other executable program. You can save variables from the workspace in ASCII format with option save filename.dat filename.dat variable variable -ascii To load variables you can use load command. load FILENAME To clear variables you can use load command. clear

1.7

Help

help. To learn more about a function you can use help. >> help for If you do not remember the exact name of a function you can use lookfor >>lookfo >>lookforr sv

Chapter 2

Single Degree of Freedom System In this chapter we will study the responses of systems with a single degree of freedom. It is important topic to master, since the complicated multiple-degree-of-freedom systems(MDOF) can often treated treated as if they are simple collections collections of several several single-degre single-degree-of-f e-of-freedo reedom(SDO m(SDOF) F) systems. systems. Once the responses of SDOF are understood, the study of complicated MDOF becomes relatively easy.

2.1 2.1.1 2.1.1

Free Vibrat Vibrations ions of Single-Degr Single-Degree-ofee-of-F Freedom reedom Systems Systems Visco Viscous us Dampi Damping ng

For the free vibration of a single-degree-of-freedom system with mass m, spring constant k, and viscous damping c, the system undergoes a dynamic displacement x(t) measured from the static equilibrium equilibrium position position of the mass. Applying Applying Newton’s law, law, the equation equation of motion of the system system is represented by

x c F m k Figure 2.1: Single degree of freedom system.

7

8

m¨ x + cx˙ + kx = 0

(2.1.1)

subject to the initial conditions x(0) = x0 and x(0) ˙ = v0 . If we divide (2.1.1) by m we can reexpress it in terms as x ¨ + 2ζωn x˙ + ωn2 x = 0 (2.1.2) where ωn =

√ 

 

k/m is natural angular frequency and ζ  = c/(2 km) is the damping ratio.

To solve the damped system of equation (2.1.2), assuming x = Aest

(2.1.3)

Substituting equation (2.1.3) into equation (2.1.2) yields an algebraic equation in the form s2 + 2ζωn s + ωn2 = 0

(2.1.4)

The solutions of equation (2.1.4) are given by s1,2 =

−ζω n ± ωn

There are three possible cases:

 

(ζ 2

− 1)

(2.1.5)

(a) Overdamped Motion In this case, the damping ratio is greater than 1 (ζ > 1). The discriminant of equation (2.1.5) is positive, resulting in a pair of distinct real roots. The solution of equation (2.1.1) then becomes x(t) = e−ζω t (Ae−ω n

n

√ ζ  −1t 2

+ Be ω

n

√ ζ  −1t 2

)

(2.1.6)

which represents that the vibration will not occur since the damping force is so large that the restoration force from the spring is not sufficient to overcome the damping force. (b) Underdamped Motion In this case the damping ratio is less than 1 (0 < ζ < 1) and the discriminant of equation (2.1.5) is negative, resulting in a complex conjugate pair of roots. The solutions for this case can be expressed as x(t) = e−ζω t (Aejω 1−ζ  t + Be −jω 1−ζ  t ) n

n

√ 

= e−ζω t (Aejω n

2

d

t

n

√ 

2

+ Be −jω t ) d

= e−ζω t (C cos ωd t + D sin ωd t)

(2.1.7)

n

= Xe−ζω

where j =

n

t

sin(ωd t + φ)

√ −1, X  and φ are constants. The the damped natural frequency is denoted by ωd = 1 − ζ 2 ωn

 

(c) Critically Damped Motion In this last case, the damping ratio is exactly 1(0zeta = 1). The solution takes the form x(t) = (A + Bt)e−ω

n

t

where the constants A and B are determined by the initial conditions.

(2.1.8)

9

Homework 2.1.1. Evaluate the constants A and B in equation (2.1.7) using the initial conditions x(0) = x0 and v(0) = v0 . Homework 2.1.2. Describe the definition of logarithmic decrement in free vibration.

2.2 2.2.1

Forced Vibration Direct Force Excitation

For a single-degree-of-freedom system with viscous damping and subjected to a forcing function F (t) as shown in figure 2.1, the equation of motion can be written as m¨ x + cx˙ + kx = F (t)

(2.2.1)

The complete solution to equation (2.2.1) consists of two parts, the homogenous solution (the complementary solution) and the particular solution. The homogenous solution is the same as the free vibration which was described in last section. It is often common to ignore the transient part of  the total solution and focus only on the steady-state response. Taking Laplace transformation of a second order differential equation with zero initial conditions, the transfer function is X (s) 1/m = 2 F (s) s + 2ζωn s + ωn2 where ωn =

(2.2.2)

 

k/m, ζ  = c/2mωn .

Substituting jω for s to calculate the frequency response, where j is the imaginary operator: X ( jω) 1/mω2 = [(ωn /ω)2 1] + 2 jζ (ωn /ω) F ( jω)

−

(2.2.3)

Example 2.2.1. Plot the amplitude and phase angle of the single degree of freedom system.  Example MATLAB Code ————————————————————————————————– clf; clear all; m = 1; zeta = 0.1:0.1:1; k = 1; wn = sqrt(k/m);

10

w = logspace(-1,1,400); rad2deg = 180/pi; s = j*w; for cnt = 1:length(zeta) xfer(cnt,:)=(1/m) ./ (s. 2 + 2*zeta(cnt)*wn*s + wn 2);

∧

∧

mag(cnt,:) = abs(xfer(cnt,:)); phs(cnt,:) = angle(xfer(cnt,:))*rad2deg; end for cnt = 1:length(zeta) figure(1) loglog(w,mag(cnt,:),’k-’) title(’SDOF frequency response magnitudes for zeta = 0.2 to 1.0 in steps of 0.2’) xlabel(’Frequency(rad/sec)’) ylabel(’Magnitude’) grid hold on end hold off  for cnt = 1:length(zeta) figure(2) semilogx(w,phs(cnt,:),’k-’) title(’SDOF frequency response phases for zeta = 0.2 to 1.0 in steps of 0.2’) xlabel(’Frequency(rad/sec)’) ylabel(’Phase’) grid hold on end hold off  ——————————————————————————————–

11

Figure 2.2: SDOF

2.2.2

Base Excitation

Often, machines are harmonically excited through elastic mounting, which may be modeled by springs and dashpots. For example, an automobile suspension system is excited by road surface. Consider the single degree of freedom system in Figure 2.3(a). The structure with mass m is connected to the base by stiffness, k, and damping with viscous damping coefficient c. The equation

12

c x

x

m

c

m

k

k y

y (a)

(b)

Figure 2.3: Free diagram of base excited single degree of freedom system. of motion is m¨ x + c(x˙

− y)˙ + k(x − y) = 0

(2.2.4)

(a)Derive the displacement transmissibility , X/Y  and plot the magnitude and phase. (b) The transmitted force by the base excitation to the structure is F T  = k(x y) + c(x˙ y). ˙ The force transmissibility , F T /kY  is defined as the dimensionless relation between maximum base displacement Y  and the transmitted force magnitude F T .Derive the force transmissibility and plot as function of frequency ratio.

−

−

[Homework2]Sky hook damper Consider the single degree of freedom system in figure 2.3(b).The structure with mass m is connected to the base by stiffness, k. Let us suppose that the viscous damping with viscous damping coefficient c is connected to the sky. (a)Derive the displacement transmissibility , X/Y  and plot the magnitude and phase. (b) Derive the force transmissibility, F T  = k(x y), and plot as function of frequency ratio.

−

2.3 2.3.1

Simulation with MATLAB Transfer Function

The linear time invariant(LTI) systems can be specified by transfer functions. The corresponding command is : sys=tf(num,den) The output sys is a model-specific data structure. Example 2.3.1. Sample Matlab code for Bode plot

13

m=1 zeta=0.1 k=1 wn=sqrt(k/m) den=[1 2*zeta*wn wn 2]

∧

num=[1/m] sys=tf(num,den) bode(sys) Example 2.3.2. The function lsim simulates the response to more general classes of inputs. For example, t=0:0.01:50; u=sin(t); lsim(sys,u,t) simulates the response of the linear system sys to a sine wave for a duration of 50 seconds.

Table 2.1: Basic Commends for Time and Frequency Response Command bode(sys) nyquist(sys) step(sys) impulse(sys) initial(sys, x0) lsim(sys,u,t,x0)

2.3.2

Explanations Bode plot Nyquist plot step response impulse response undriven response to initial condition response to input u

State Space Analysis

It is desirable to change the system equation for an n d.o.f system with n second order differential equation to 2n first order differential equations. The first order form of equations for the system is called as state space form. Start by solving equation second order differential equations. m¨ x + cx˙ + kx = F (t)

(2.3.1)

14

we define the state vector as

T 



(2.3.2)

˙ x(t) = Ax(t) + BF 

(2.3.3)

x(t) =



x(t) x(t) ˙

Then, adding the identity x˙ = x, ˙ equation (2.3.1) can be written in the state form as

where the system matrix A and the input matrix B are : A=



0

1

−m−1k −m−1c

and B=



(2.3.4)

  0

(2.3.5)

m−1

Schematically, a Single Input Single Output(SISO) state space system is represented as shown in

Direct Transmission Matrix D dx(t)/dt

f(t) Input

B

+

Integrate

x(t)

C

+

y(t) Output A System Matrix

Figure 2.4: State space system block diagram Figure (2.4). The scalar input u(t) is fed into both the input matrix B and the direct transmission matrix D. The output of the input matrix is a n 1 vector, where n is the number of states. The output is fed into a summing junction to be added to the output of the C matrix.

×

The output of the B matrix is added to the feedback term coming from the system matrix and is fed intro an integrator block. The output matrix has as many rows as outputs, and has as many columns as states, n.

15

To account for the case where the desired output is not just the states but is some linear combination of the states, and output matrix C is defined to relate the outputs to the states. Also, a matrix D, know as the direct transmission matrix, is multiplied by the input F (t) to account for outputs that are related to the inputs but that bypass the states. y(t) = Cx(t) + DF 

(2.3.6)

The output matrix C has as many rows as outputs required and as many columns as states. The direct transmission matrix D has the same number of columns as the input matrix B and as many rows as the output matrix C. Example 2.3.3. Numerically compute the free vibration of mass-spring-damper system using initial function in MATLAB.  Example MATLAB Code m=1;d=0.1;k =1; A=[0 1;-k/m -c/m]; C=[1 0]; sys=ss(A,[],C,[]); x0=[10,0]; initial(sys,x0)

Figure 2.5: Initial condition results The result of free vibration of the one degree of freedom system is shown in figure 2.5.

16

Homework 2.3.4. One of the common excitation in vibration is a constant force that is applied for a short period of time and then removed. Numerically calculate the response of mass-spring-dashpot system to this excitation in MATLAB. m¨ x + cx˙ + kx = F o [1

− H (t − t1)]

where H is Heaviside function. stepfun is useful command to solve this problem.

(2.3.7)

Chapter 3

Multiple Degree of Freedom Systems

3.1 3.1.1

Some Basics Concepts for Linear Vibrating System Eigenvalue Problem

In the previous chapters a single degree of freedom system with a single mass, damper and spring was considered. Real systems have multiple degrees of freedom and their analysis is complicated by the large number of equations involved. To deal with them, matrix are used. The equation of  motion for n-degree of freedom equation can be written as ¨ + [c] x˙ + [k] x = [bf ] f  [m] x

{}

{}

{}

{}

(3.1.1)

where the mass [m], damping [c], and stiffness [k] matrices are symmetric. First consider undamped vibration without excitation force. The system can be solved by assuming a harmonic solution of the form x = uejωt (3.1.2) Here,u is a vector of constants to be determined, ω is a constant to be determined. Substitution of this assumed form of the solution into the equation of motion yields

−

ω 2 M  + K  uejωt = 0



(3.1.3)

Note that the scalar ejωt = for any value of  t and hence equation (3.1.1) yields the fact that ω and u must satisfy the vector equation ( ω 2 M  + K )u = 0



−

17

18

Note that this represents two algebraic equations in the three unknowns; ω, u1 , u2 where u = [ u1 u2 ]T . This equation is satisfied for any u if the determinant of the above equation is zero.

−

ω 2 M  + K  = 0



(3.1.4)

The simultaneous solution of equation (3.1.4) results in the values of parameter ω 2 . The ω is called as eigenvalues  of the problem. Once the value of  ω is established, the value of the constant vector u can be found by solving equation (3.1.3). Example 3.1.1. Consider the system with two masses represented in figure 3.1.

x2

x1 c2

c1

c3

m1

m2

k1

k3

k2 F2

F1

Figure 3.1: 2dof  The equations of motion become



m1

0

0

m1

   x ¨1

c1 + c2

+

x ¨2

x˙ 1

m1

0

0

m1

c2 + c3

x˙ 2

+

k1 + k2

   x ¨1

+

x ¨2

k1 + k2

−k2

−k2

−k2

−k2

   x1

k2 + k3

k2 + k3

x2

=

F 1 (t)



F 2 (t) (3.1.5) To determine the natural frequencies and natural mode shapes of the system, the undamped free vibration of the system is first considered. Thus the equations reduce to



−c2

  

−c2

    x1 x2

=

0 0

(3.1.6)

Consider a numerical example for the system shown in figure (3.1). Let c1 = c2 = c3 = 0, m1 = 5kg, m2 = 10kg, k1 = 2N/m, k2 = 2N/m, k3 = 4N/m. Substituting in equation (3.1.6) yields



5

0

0 10

   x ¨1 x ¨2

+

4

−2

−2 6

    x1 x2

=

0 0

(3.1.7)

19

Assume harmonic responses of the of the form x1 = A1 exp(iωt) and x2 = A2 exp(iωt). Equation (3.1.6) becomes 5 0 A1 4 2 A1 ω2 = (3.1.8) 0 10 A2 2 6 A2



  

−

−

 

Solving Eigenvalue Problem with MATLAB The eigenvalue problem of a matrix is defined as Au = λu

(3.1.9)

Ku = λMu

(3.1.10)

and generalized eigenvalue problem is The eig is subroutine for computing the eigenvalues and the eigenvectors of the matrix A or >>[V,D]=eig(A) >>[V,D]=eig(K,M) The eigenvalues of system are stored as the diagonal entries of the diagonal matrix D and the associated eigenvectors are stored in columns of the matrix V .



Example MATLAB Code

m=[5 0 ;0 10]; k=[4 -2;-2 6]; [v,d]=eig(k,m) The function eig in MATLAB gives unsorted eigenvalues, so it will be help to make sorting the eigenvalues of the system. function [u,wn]=eigsort(k,m); Omega=sqrt(eig(k,m)); [vtem,d]=eig(k,m); [wn,isort]=sort(Omega); il=length(wn); for i=1:il v(:,i)=vtem(:,isort(i)); end disp(’The natural frequencies are (rad/sec)’) disp(’ ’) wn disp(’ ’) disp(’The eigenvectors of the system are’) v

20

—————————————————————— The two natural frequencies are ω1 =0.6325 rad/s, ω2 =1 rad/s The eigenvectors are u1 =

  1 1

T 

, u2 =



T 

1



−0.5

1

m1

1

m2

1

m1

m2 -1/2

Figure 3.2: Mode shapes for the two degree of freedom system

3.1.2

Orthogonality of normal modes

The modes are orthogonal with respect to the mass matrix and stiffness matrix.

{u}T 2 [m]{u}1 = 0 {u}T 2 [k]{u}1 = 0

(3.1.11)

Mass normalizing equation (3.1.11), we can get the general relations as

{u}T i[m]{u}i = mi, i = 1, 2 {u}T i[k]{u}i = mi ωi2 = ki, i = 1, 2

(3.1.12)

where mi and ki is called modal mass and modal stiffness for the i-th modal vector of vibration. The numerical values of the mode shape will be used to determine the modal mass and modal stiffness. The mode shapes were found to be u1 =



T 

1

−0.5



for ω1 = 1 rad/s

Verification with MATLAB u1 = v(:, 1); u2 = v(:, 2); u1 [m] u1 = 15

∗

∗

T 

  1 1

for ω1 =

 

2/5 rad/s, and u2 =

21

u1 [m] u2 [m] k1 = ω12 k2 = ω22

∗ ∗

3.1.3

∗ u2 = 0 ∗ u2 = 7.5 ∗ m1 = 6 ∗ m2 = 15/2 Normalization of Mode Shapes

While above relations are related to the mass and stiffness of the modal space, it is important to remember that the magnitude of these quantities depends upon the normalization of the modal vectors. Therefore, only the combination of a modal vector together with the associated modal mass and stiffness represent a unique absolute characteristic concerning the system being described. When we scaled the eigenvector such that mi = 1, the equation (3.1.12) becomes

{u}T i[m]{u}i = 1, {u}T i[k]{u}i = ωi2,

i = 1, 2

(3.1.13)

i = 1, 2

This meas that mi is not unique. There are several ways to normalize the mode shapes. (1) The mode shapes can be normalized such that the modal mass mi is set to unity. (2) The largest element of the mode shape is set to unity. (3) A particular element of the mode shape is set to unity. (4) The norm of the mode vector is set to unity.

Example 3.1.2. Using the previous two degree of freedom example, normalize the modal vectors

{ }T i[m]{u}i = 1,

such that u

i = 1, 2 The mass normalization of the first and second natural modes

are

{u}1 = {u}2 = √ 1m2

  √    √     −    1

1 m1

=

1

1

1/2

=

1 15

      − 

1 15/2

1 1

1

1/2

The orthogonality of modes permit us to transform the coupled equations of motion defined in physical coordinate to uncoupled system in the modal coordinate.

22

3.1.4

Modal Coordinates

In solving the equations of motion for an undamped system (3.1.6), the major obstacle encountered when trying to solve for the system response x for a particular set of exciting forces and initial conditions, is the coupling between the equations. The coupling is seen in terms of non-zero off  diagonal elements. If the system of equations could be uncoupled, so that we obtained diagonal mass and stiffness matrices, then each equation would be similar to that of a single degree of freedom system, and could be solved independent of each other. The process of deriving the system response by transforming the equations of motion into an independent set of equations is known as modal analysis Thus the coordinated transformation we are seeking, is one that decouples the system. The new coordinate system can be found referring to orthogonal properties of the mode shapes discussed in equation (3.1.12) and (3.1.13). n

{x(t)} =

{ }

u i q i (t)

(3.1.14)

i=1

where the physical coordinate, x(t) are related with the normal modes, decoupled coordinate, q i .

{ }

{u}i

and the normal

Equation (3.1.14) may be written in matrix form as

{x(t)} = [P ]{q (t)}

(3.1.15)

where [P ] is called the modal matrix . Thus, the modal matrix for a 2-DOF system can appear as [P ] = [ u

{ }1 {u}2

]

(3.1.16)

Substituting equation (3.1.16) into the general equation (3.1.1), we obtain as [m][P ] q¨  + [c][P ] q ˙ + [k][P ] q  = f 

{}

{}

(3.1.17)

{} {}

Multiplying on the left by [P ]T , [P ]T [m][P ] q¨  + [P ]T [c][P ] q ˙ + [P ]T [k][P ] q  = [P ]T  f 

{}

{}

{}

{}

(3.1.18)

We know that orthogonality of the modes with respect to mass and stiffness matrices. Assuming that the viscous damping can be decoupled by modal matrix, we obtain q¨i  (t) + 2ζ i ωi q ˙(t) + ωi2 q i (t) = N i (t), where N i (t) is

The ratio in equation (3.1.20)

i = 1, 2,

·· ·

(3.1.19)

{ { u}T  u}T  i {f (t)} i {f (t)} N i (t) = = T  mi {u}i [m]{u}i

(3.1.20)

{u}T i {u}T i[m]{u}i

(3.1.21)

23

is called modal participation factor . The displacement can be expressed as

∞

∞

{u}i{u}T i{f (t)} (3.1.22) [(ωi2 − ω 2 ) + 2iζ i ωi ] m i i=1 i=1 where ωi is the natural frequency in the i-th mode. If eigenvector {u}i is mass normalized, {u}T  i [m]{u}i = x=

{ }  u i q i =

1.

Numerical Simulation with MATLAB

     1 0

[P ]T  [m] [P ] =

∗

∗

[P ]T  [k] [P ] =

∗ ∗

[P ]T  =



0.2582

0 1

0.4 0 0

1

0.2582

0.3651 0.1826 Thus the equations of motion in modal coordinate are

−

    1 0

q¨1 

0 1

q¨2 

0.4 0

+

0

1

   q 1

=

q 2

0.2582f 1 + 0.2582f 2 0.3651f 1

− 0.1826f 2

   

=

f 1 

f 2

(3.1.23)

The matrix equation of (3.1.23) can be written in terms of algebraic differential equations q¨ 1 + 0.4q 1 = f 1 

(3.1.24)

q¨ 2 + q 2 = f 2 

Hence, the system equations have been uncoupled by using the modal matrix as a coordinate transformation. Example 3.1.3. Calculate the response of the system illustrated in figure (3.3) to the initial displacement x(0) =

  1 1

T 

with x(0) ˙ =

  0

T 

using modal analysis.

0

The initial conditions in modal space become

{q (0)} = [P ]T  ∗ {x(0)} = T 



0.5164

˙ }= {q ˙(0)} = [P ] ∗ {x(0) The modal solution of equation (3.1.24) is

0.1826

  0

T 



0

q 1 (t) = q 1 (0) cos(ω1 t) = 0.5164 cos(0.6325t)

T 

24

x1 2

x2 2

4

5 f1

10 f2

q1 2/5

1 f 1’

q2 1

1 f 2’

Figure 3.3: The undamped two degree of freedom system and broken down to two single degree of  freedom systems q 2 (t) = q 2 (0) cos(ω2 t) = 0.1826 cos(t) Using the transformation x(t) = P q (t) yields that the solution in physical coordinates is x(t) =

3.2



0.1333 cos(0.6325t) + 0.0667 cos(t) 0.1333 cos(0.6325t)

− 0.0333 cos(t)



Proportional Damping

Damping is present in all oscillatory systems. As there are several types of damping, viscous, hysteretic, coulomb etc., it is generally difficult to ascertain which type of damping is represented in a particular structure. In fact a structure may have damping characteristics resulting from a combination of all types. In many cases, however, the damping is small and certain simplifying assumptions can be made. The most common model for damping is proportional damping defined as [c] = α[m] + β [k] (3.2.1) where [c] is damping matrix and α, β  are constants. For the purposes of most practical problems, the simpler relationship will be sufficient. Caughey1 showed that there exists a necessary and sufficient condition for system (3.1.1) to be completely uncoupled is that [m]−1 [c] commute with [m]−1 [k]. ([m]−1 [c])([m]−1 [k]) = ([m]−1 [k])([m]−1 [c]) 1

(3.2.2)

T.K. Caughey, ”Classical Normal Modes in Damped Linear Systems”, Journal of Applied Mechanics, Vol 27, Trans. ASME, pp.269-271, 1960

25

Figure 3.4: Time response of mass 1 and mass 2 or [c][m]−1 [k] = [k][m]−1 [c]

3.3

(3.2.3)

Modal Analysis of the Force Response

The forced response of a multiple-degree-of-freedom system can also be calculated by use of modal analysis. Example 3.3.1 (See Example 4.6.1 in Inman, pp.296). For this example, let m1 = 9kg, m2 = 1kg,k1 = 24N/m, and k2 = 3kg. Assume that the damping is proportional with α = 0 and β  = 0.1, so that c1 = 2.4Ns/m, and c1 = 0.3Ns/m. Also assume that F 1 = 0, and F 2 = 3cos2t. Calculate the steady-state response.

         − 9 0

x ¨1

0 1

x ¨2

+

2.7 0.3

−0.3 0.3

       − x˙ 1 x˙ 2

+

27 3

−3 3

            x1 x2

=

0

0

F 1

0

1

F 2 (3.3.1)

26

Numerical Simulation with MATLAB [P ] =

−

0.2357

−0.7071 [P ]T  ∗ [m] ∗ [P ] =

[P ]T  [k] [P ] =

∗

0.7071



    ∗   ∗  − 

[P ]T  [c] [P ] =

∗

−0.2357

[P ]T [B] =

0

1 0 0 1

0.2

0

0

0.4

2 0 0 4

0.7071

0 0.7071 Hence the decoupled modal equations become q¨ 1 + 0.2q + ˙ 2q 1 =

−0.7071 ∗ 3 ∗ cos2t q¨ 2 + 0.4q + ˙ 4q 2 = 0.7071 ∗ 3 ∗ cos2t

(3.3.2)

Comparing the coefficient of q ˙i to 2ζ i ωi yields .2 ζ 1 = 20√  2

ζ 2 = 20∗.22 Thus the damped natural frequencies becoms ωd1 = ω1 ωd2 = ω2

  −   − 1 1

ζ 12

 1.41 ζ 22  1.99

Note that while the force F 2 is applied only to mass m2 , it becomes applied to both coordinate when transformed to modal coordinates. Let the particular solutions of equations (3.3.2) be q 1 p and q 2 p .

x2

x1 c2

c1 m1 k1

m2 k2

F1

F2

Figure 3.5: Damped two-degree-of-freedom system

27

The steady state solution in the physical coordinate system is xss (t) = [P ]q p (t) =

3.4

−

0.2357q 1 p (t)

− 0.2357q 2 p (t)

−0.7071q 1 p (t) + 0.7071q 2 p (t)



State-Space Approach

Simulation by state-space method is a much easier way to obtain the systems response when compared to computing the response by modal analysis. However, the modal approach is needed to perform design and to gain insight into the dynamics of the system. In this section the simulation method for free vibration and forced vibration by state space formulation will be discussed.

3.4.1

Free Vibration

Consider the forced response of a damped linear system. The most general case can be written as ¨ + [c] x˙ + [k] x = 0 [m] x

{}

{}

(3.4.1)

{}

with initial condition x(0) = x0

˙ x(0) = x˙ 0

Again it is useful to write this expression in a state-space form by defining the two n ˙ then the equations (3.4.1) becomes y1 = x and y2 = x, ˙ y(t) = Ay(t) where A=



0

× 1 vectors (3.4.2)

I 

−m−1k −m−1c



(3.4.3)

The eigenvalues λi will appear in complex conjugate pairs in the form λi =

−ζωi − jω i λi = −ζωi + jω i

  −   − 1

ζ i2

1

ζ i2

(3.4.4)

Example 3.4.1. Consider the system shown in figure 3.1. Calculate the response of the system to the initial condition using state-space method. Let c1 = c3 = 0, c2 = 0.2N s/m, m1 = 2kg,

·

m2 = 1kg, k1 = 0.2N/m, k2 = 0.05N/m, k3 = 0.05N/m. The initial condition of  m1 is 0.1m and let the other parameters be all zero. 

Example MATLAB Code

28

—————-dof2ini.m————————m1=2;m2=1; d1=0; d2=0.2; d3=0; k1=0.2; k2=0.05;k3=0.05; m=[m1 0;0 m2]; d=[d1+d2 -d2; -d2 d2+d3]; k=[k1+k2 -k2;-k2 k2+k3]; A=[zeros(2,2),eye(2);-inv(m)*k,-inv(m)*f]; C = [1 0 0 0]; x0=[0.1 0 0 0]; sys=ss(A,[],C,[]) initial(sys,x0) The simulation result is shown in figure (3.6). Response to Initial Conditions 0.1

0.08

0.06

0.04

0.02   e    )    d    1   u    (    t    i    l    Y   p   :   o   m   T    A

0

−0.02

−0.04

−0.06

−0.08

−0.1

0

50

100

150

200

250

300

350

400

450

500

Time (sec)

Figure 3.6: Time response of mass 1

Homework 5 Consider the system shown in figure 3.5. Let m1 = 10kg, m2 = 1kg, k1 = 0.5N/m, k2 = 0.05N/m, and c1 = 0, c2 = 0.2N s/m. The initial condition of  m1 is 0.1m and let the other parameters be all zero. Plot the transient response of mass 1.

·

29

3.4.2

Forced Vibration

¨ + [c] x˙ + [k] x = [bf ] f  [m] x

{}

{}

{}

(3.4.5)

{}

with initial condition ˙ x(0) = x˙ 0

x(0) = x0 The state-space equations are

˙ y(t) = Ay(t) + Bf  where A=



0

I 

−m−1k −m−1c

and B=



0 m−1 bf 

(3.4.6)



(3.4.7)



(3.4.8)

Example 3.4.2. Compare the frequency response function of the two degree of freedom shown in figure (3.1) between c2 = 0 and c2 = 0.2N s/m. The other parameters are as follows. Let

·

m1 = 2kg,m2 = 1kg, k1 = 0.2N/m, k2 = 0.05N/m, k3 = 0.05N/m and c1 = c3 = 0 N s/m and the

·

excitation force F 2 be zero. 

Example MATLAB Code

—————-dof2frf.m————————m1=2;m2=1; c1=0; c2=0.0; c3=0;k1=0.2; k2=0.05;k3=0.05; Bf=[1; 0]; m=[m1 0; 0 m2]; damp=[c1+c2 -c2; -c2 c2+c3]; K=[k1+k2 -k2;-k2 k2+k3]; A=[zeros(2,2),eye(2);-inv(m)*k,-inv(m)*damp]; B=[zeros(2,1); inv(m)*bf];C = [1 0 0 0]; D=zeros(size(C,1), size(B,2)) sys=ss(A,B,C,D) d1=0; d2=0.02; d3=0; damp1=[c1+c2 -c2; -c2 c2+c3]; Adamp=[zeros(2,2),eye(2);-inv(m)*k,-inv(m)*damp1]; sysdamp=ss(Adamp,B,C,D) w=linspace(0.1, 1, 800)

30

bode(sys,sysdamp,w) —————————— the result of simulation is shown in figure (3.7). Bode Diagram 80 60    )    B    d    (   e    d   u    t    i   n   g   a    M

40 20 0 −20 −40 360 270

   )   g   e    d    (   e   s   a    h    P

180 90 0 −90 −180 −1

0

10

10 Frequency (rad/sec)

Figure 3.7: Frequency response fucntion of mass 1

Chapter 4

Design for Vibration Suppression

4.1

Introduction

A Dynamic Vibration Absorber (DVA) is a device consisting of a reaction mass, a spring element with appropriate damping that is attached to a structure in order to reduce the dynamic response of the structure. The frequency of dynamic absorber is tuned to a particular structural frequency so that when that frequency is excited external force. The concept of DVA was first applied by Frahm in 1909 to reduce the rolling motion of ships as well as hull vibrations. A theroy for the DVA was presented later by Ormondroyd and Den Hartog (1928) 1 . The detailed study of optimal tuning and damping parameters was discussed in Den Hartog’s on Mechanical Vibration (1940) book 2 .

4.2

Vibration Absorber

Figure (4.1) shows a SDOF system having mass m and stiffness k, subjected to external forcing. A dynamic absorber with mass m2 , stiffness k2 , and dashpot c2 is attached to the primary mass. What we now have is a 2DOF problem rather than the original SDOF. The m1 , c1 , k1 system is referred to as the primary system, and m2 , c2 , k2 system is known as the secondary system. The displacement of primary mass and absorbing mass are x1 and x2 , respectively. With this notation, the governing equations take the form 1

J.Ormondroyd, and J.P.Den Hartog, ”The theory of the dynamic vibration absorber”, Trans. ASME, 50, 1928, pp. 9-15 2 J.P.Den Hartog, Mechanical vibration, Dover, 4th ed. Reprint, 1984

31

32



m1

0

0

m2

   x ¨1 x ¨2

+

c1 + c2

− c2

− c2 c2

   x˙ 1

+

x˙ 2

k1 + k2

−k2

−k2 k2

    x1

=

x2

F (t) 0

(4.2.1)

We set xi = Re[X i ejωt ] for steady-state response, which leads to the following complex amplitudefrequency equations,

x2

x1 c2

c1 m1 k1

m2 k2

F1 Figure 4.1: SDOF system coupled with a DVA

 4.2.1

[(k1 + k2 ) + jω(c1 + c2 )

−(k2 + jωc2)

− m1ω2]

−(k2 + jωc2) k2 + jωc2 − m2 ω 2

    X 1 X 2

=

F  0

(4.2.2)

SDOF with Undamped DVA

Let us consider the case where damping is negligible, c1 = c2 = 0. We then find from equation (4.2.2) that (k2 m2 ω 2 ) X 1 = F  (4.2.3) [(k1 + k2 ) ω 2 m1 ](k2 m2 ω 2 ) k22

−

X 2 = F 

[(k1 + k2 )

−

−

k2 ω 2 m1 ](k2

−

−

− m2ω2) − k22

(4.2.4)

where the determinant of the system of coefficients in equation is ∆(ω) = [(k1 + k2 )

− ω2m1](k2 − m2ω2) − k22.

(4.2.5)

First, note from equation (4.2.3) that the magnitude of steady-state vibration, x1 becomes zero when the absorber parameters k2 and m2 is chosen to satisfy the tuning condition ω2 = k2 /m2

(4.2.6)

33

In this case the steady-state motion of the absorber mass is calculated from equation (4.2.4) X 2 =

− kF 2

(4.2.7)

With the main mass standing still and the secondary mass having a motion F/k2 exp( jωt) the force in the damper spring varies as Fexp( jωt), which is actually equal and opposite to the external force.

−

−

For simplifications we want to bring equation (4.2.3) and (4.2.4) into a dimensionless form and for that purpose we introduce the following parameters: xst = F/k1 ; static deflection of primary system. ω12 = k1 /m1 ; natural frequency of primary system ω22 = k2 /m2 ; natural frequency of secondary system µ = m2 /m1 ; mass ratio=secondary mass/primary mass With this definitions, also note that

k2 ω2 = µ 22 = µf 2 k1 ω1

(4.2.8)

where frequency ratio f  is f  = ω2 /ω1 . Then equationa (4.2.3)and (4.2.4) becomes X 1 = xst [1 X 2 = xst [1 where frequency ratio β  is β  = ω/ω 2 .

−

1 β 2 β 2 ][1 + µf 2 (f β )2 ]

−

−

− µf 2

1

− β 2][1 + µf 2 − (f β )2] − µf 2

(4.2.9) (4.2.10)

The absolute value of this system is plotted in figure (4.2) for the case µ = 0.25. In fact, if the driving frequency shifts such that X 1 /xs t > 1, the force transmitted to the primary system is amplified and the absorber system is not an improvement over the original design of the primary system.

|

|

Let us consider the case that the frequency ratio f  = 1, (i.e., ω2 = ω1 , or k2 /m2 = k1 /m1 ) For this special case, equations (4.2.9) and (4.2.10) becomes X 1 = xst (1 X 2 = xst (1

− −

1 β 2 β 2 )(1 + µ β 2 )

−

β 2 )(1

1 +µ

−

−µ

− β 2) − µ

(4.2.11) (4.2.12)

The natural frequencies are determined by setting the denominators equal to zero : (1

− β 2)(1 + µ − β 2) − µ = 0 β 4 − f β 2 (2 + µ) + 1 = 0

(4.2.13)

34

β

Figure 4.2: Plot of normalized magnitude of the primary mass versus the normalized driving frequency The solutions are β 2 = (1 + µ/2)

±

 

µ + µ2 /4

(4.2.14)

This relation is plotted graphically in figure (4.3). Note that as µ is increased, the natural frequencies split farther apart. Homework 4.2.1. Inman’s Book Example 5.3.1 Homework 4.2.2. Inman’s Book Example 5.3.2

4.2.2

SDOF with damped DVA

The equations of motion are given in matrix form by equation (4.2.1). Note that these equations cannot necessarily be solved by using the modal analysis technique of Chapter 3 because the equations do not decouple (KM −1 C  = CM −1 K ).



The steady-state response of equation (4.2.1) can be obtained by assuming a solution of the form by assuming a solution of the form

x(t) = X(t)ejωt =

  X 1 X 2

ejωt

(4.2.15)

where X 1 is the amplitude of vibration of the primary mass and X 2 is the amplitude of vibration of  the absorber mass. From the (4.2.2), we obtain that

35

     µ

β

Figure 4.3: Plot of mass ratio versus system natural frequency(normalized to the frequency of the secondary system

X 1 (k2 m2 ω 2 ) + jc 2 ω = F  det([K ] ω 2 [M ] + jω[C ])

(4.2.16)

X 2 k2 + jc 2 ω = F  det([K ] ω 2 [M ] + jω[C ])

(4.2.17)

− − −

which expresses the magnitude of the response of the primary mass and secondary mass, respectively. Note that these values are complex numbers. First, consider the case for which the internal damping of the primary system is neglected ( c1 = 0). Using complex arithmetic, the amplitude of the motion of the primary mass can be written as the real number X 12 = F 2 [(k1

− m1ω2)(k2 −

(k2 m2 ω2 )2 + (c2 ω)2 m2 ω 2 ) m2 k2 ω 2 ]2 + [k1 (m1 + m2 )ω 2 ]2 c22 ω 2

− −

−

(4.2.18)

It is instructive to examine this amplitude in terms of the dimensionless ratios introduced for the undamped vibration absorber. The amplitude x1 is written in terms of the static deflection xst = F/k of the primary system. In addition, consider the mixed damping ratio defined by ζ  = where ω =

 

c2 2m2 ω1

(4.2.19)

k1 /m1 . Using the standard frequency ratio r = ω/ω1 , the ratio of natural frequencies

36

f  = ω2 /ω1 , and the mass ratio µ = m2 /m1 , equation (4.2.18) can be written as X 1 = xst

 

(2ζr)2 (r 2

−

(2ζr)2 + (r2 f 2 )2 1 + µr2 )2 + [µr2 f 2 (r2

−

−

− 1)(r2 − f 2)]2

(4.2.20)

which expresses the dimensionless amplitude of the primary system. Note that the amplitude of the primary system response is determined by four physical parameter values: mass ratio µ, the ratio of the decoupled natural frequencies f  the ratio of the driving frequency to the primary natural frequency r and the damping ratio of absorber ζ . These four parameters can be considered as design variables and are chosen to give the smallest possible value of the primary mass’s response, x1 for a given application. It is instructive to verify this result for several particular cases Homework 4.2.3. Plot the compliance curves when the parameters of absorbers are: 1. ζ 2 =

∞

2. ζ 2 = 0 3. ζ  = 0.10 4. ζ  = 0.10 for mass ratio µ = 1/20 and frequency ratio f  = 1. The dynamic vibration absorber said to be optimally tuned and damped when the two resonance peaks are equal in magnitude. The optimal frequency ratio f  and damping ratio are given as f  = ζ  =

1 1+µ

 

3µ 8(1 + µ)3

(4.2.21) (4.2.22)

Homework 4.2.4. Obtain the optimal frequency ratio f  and optimal damping ratio ζ  when mass ratio is µ = 1/20.

4.3 4.3.1

Isolation Design Passive Isolators

In figure (4.5), a single-degree-of-freedom vehicle model is shown with (a) passive, (b) skyhook suspensions. The passive system using linear elements has the equation of motion. x ¨ + 2ζωn (x˙ 1

− x˙ 0) + ωn2 (x1 − x0) = 0

(4.3.1)

37

16

zeta=inf. 14

=1/20 f=1

12

zeta=0.32

zeta=0 10

zeta=0.10

   t   s

  x 8    /    1   x 6 4 2 0 0.6

0.7

0.8

0.9

1.0

1.1

1.2

1.3

Excitation frequency ratio (r)

Figure 4.4: Amplitudes of the main mass for various values of absorber damping. All curves pass through the fixed points

x1

m k

c x0

Figure 4.5: Schematic of passive isolators

38

√ 

where ωn2 = k/m, and ζ  = c/2 km. For the base excitation problem it is assumed that the base moves harmonically such that x0 = X 0 exp( jω b t) (4.3.2) where X (0) denotes the amplitude of the base motion and ωb represents the frequency of the base oscillation. The displacement of mass divided by the amplitude of base excitation is obtained as X 1 1 + (2ζr)2 = X 0 (1 r2 )2 + (2ζr)2



−

1/2



(4.3.3)

where the frequency ratio r = ωb /ωn . The ratio is called the displacement transmissibility. Another quantity of interest in the base excitation problem is the force transmitted to the mass as the result of a harmonic displacement of the base. Hence the force transmitted to the mass is the sum of the force in the spring and the force in the damper F T (t) = k(x1 x0 ) + c(x˙ 1 x˙ 0 ) (4.3.4)

−

−

The force transmissibility is defined as F T  1 + (2ζr)2 = r2 (1 r2 )2 + (2ζr)2 kX 0



−



1/2

(4.3.5)

Figure (4.6) is a frequency response plot from equation (4.3.3).

4.3.2

Skyhook and Active Isolators

When active suspensions are used the suspension force can be generated based on control strategies. using optimal control theory and a commonly used quadratic performance criterion it was shown 3 that an optimum single-degree-of-freedom isolator must generate suspension force as F a /m =

−2ζω nx˙ 1 − ωn2 (x1 − x0) = 0

(4.3.6)

leading to a sprung mass equation of motion as x ¨ + 2ζωn x˙ 1 + ωn2 (x1

− x0) = 0

(4.3.7)

whici is the same as he governing equation of skyhook isolators shown in figure (4.7). The displacement of mass divided by the amplitude of base excitation is obtained as



X 1 = X 0 (1

−

1 2 2 r ) + (2ζr)2

1/2



(4.3.8)

where the frequency ratio r = ωb /ωn 3

Bender, E.K., Optimum linear preview control with application to vehicle suspension, ASME, Journal of basic engineering, 90(2), June 1968, pp.213-221.

39

ζ ζ ζ ζ ζ

Figure 4.6: Displacement transmissibility for passive isolator

m

c

m

x1

x1

Fa x0

k x0

Figure 4.7: Schematic of passive isolators

40

ζ ζ ζ ζ ζ

Figure 4.8: Displacement transmissibility for skyhook isolator

41

4.3.3

Semi-active Isolators

Semi-active suspensions respresents a compromise between passive and active ones. The concept of  semi-active suspension was first proposed by Crosby and Karnopp 4 in 1973 x1

m k

Fs x0

Figure 4.9: Schematic of semi-active isolators



F s = F s /m =

4



2ζω n , 0,

x˙ 1 (x˙ 1 x˙ 1 (x˙ 1

− x˙ 0) > 0

− x˙ 0) < 0



(4.3.9)

Crosby, M.J., and Karnopp, D.C., The active damper - a new concept for shock and vibration control, The shock and vibration bulletin, 43(4), June 1973, pp.119-133

Chapter 5

Vibration of strings and rods 5.1

Strings

5.2

Rods

Consider the vibation of an elastic rod (or bar) of length L and of varying cross-sectional area shown in figure (5.1). The forces on the infinitesimal element summed in the x direction are

F

x

F+dF

u(x,t) x x+dx

L Figure 5.1: Cantilevered rod in longitudinal vibration along x

F  + dF 

t) − F  = ρA(x)dx ∂u(x, 2 ∂t

(5.2.1)

where u(x, t) is the deflection of the rod in the x direction. From the solid mechanics, F  = EA(x)

42

∂u(x, t) ∂x

(5.2.2)

43

where E  is the Young’s modulus. The differential form of  F  becomes dF  =

∂F  dx ∂x

(5.2.3)

from the chain rule for partial derivatives. Substitution of equation (5.2.2) and (5.2.3) into (5.2.1) and dividing by dx yields ∂ 2 u(x, t) ∂  ρA(x) = 2 ∂t ∂x



∂u(x, t) EA(x) ∂x



(5.2.4)

When A(x) is a constant this equation becomes ∂ 2 u(x, t) = ∂t 2 The quantity v = the rod.

 

E  ∂ 2 u(x, t) ρ ∂x 2



(5.2.5)

E/ρ defines the velocity of propagation of the displacement (or stress wave) in

Chapter 6

Bending of Beam 6.1

Equation of Motion z

y x

       t 

   b

Figure 6.1: Beam The equation of motion of Euler-Bernoulli Beam is m(x)

∂ 2 w ∂w ∂ 4 w + c + EI  = f (x, t) ∂t 2 ∂t ∂x 4

(6.1.1)

where,m is mass per unit length of beam defined as m = ρA. If no damping and no external force is applied so that c = 0, f (x, t) = 0, and EI (x) and m(x) are assumed to be constant, equation (6.1.1) simplifies ∂ 2 w EI  ∂ 4 w + =0 (6.1.2) ∂t 2 m ∂x 4 Note that the free vibration equation (6.1.2) contains four spatial derivatives and hence requires four boundary conditions. The two time derivatives requires that two initial conditions, one for the displacement and one for the velocity.

44

45

6.2

Eigenvalue Problem

For the eigenvalue problem, assume the product solution as w(x, t) = W (x)F (t)

(6.2.1)

where W (x) depends on the spatial position alone and F (t) depends on time alone. Introducing equation (6.2.1) into equation (6.1.2), we can obtain the following equation as d4 W (x) dx4 where β 4 =

ω2 m EI  , 0

− β 4W (x) = 0

(6.2.2)

< x < L.

       t 

x L Figure 6.2: Cramped-Free transverse beam

6.2.1

Boundary condition

Clamped-free The boundary conditions for the clamped-free case are W (0) = 0 dW (x) x=0 = 0 dx 2 d W (x) x=L = 0 dx2 3 d W (x) x=L = 0 dx3

| | |

(6.2.3)

The solution of equation (6.2.2) is W (x) = C 1 sin βx + C 2 cos βx + C 3 sinh βx + C 4 cosh βx

(6.2.4)

Applying the boundary conditions for x = 0, we find C 2 + C 4 = 0 C 1 + C 3 = 0

(6.2.5)

so that the eigenfunction is reduced to x = L, we get C 1 (sin βL + sinh βL) + C 2 (cos βL + cosh βL) = 0 C 1 (cos βL + cosh βL)

− C 2(sin βL − sinh βL) = 0

(6.2.6)

46

Equating the determinant of the coefficients to zero, we obtain the characteristic equation



(sin βL + sinh βL)

(cos βL + cosh βL)

(cos βL + cosh βL)

−(sin βL − sinh βL)

    C 1

C 2

0

=

0

(6.2.7)

The characteristic equation is cos βL cosh βL =

−1

(6.2.8)

From the numerical analysis β 1 L = 1.875, β 2 L = 4.694, β 3 L = 7.855

     

(1.875)2

EI  mL4

rad/sec

ω2 = (4.694)2

EI  mL4

rad/sec

ω3 = (7.855)2

EI  mL4

rad/sec

ω1 =

We obtain the corresponding eigenfunctions

(6.2.9)

sin βL −sinh βL − cosh βx) + C r cos − sinh βx) βL +cosh βL (sin βx = Ar [(sin β r L − sinh β r L)(sin β r x − sinh β r x) + (cos β r L + cosh β r L)(cos β r x − cosh β r x)]

W r (x) = C r (cos βx

(6.2.10)

Example 6.2.1. The geometric and material properties are ρ(Density)

L(Length)

b(Width)

t(Thickness)

E 

2750 kg/m3

340 mm

22 mm

2 mm

7.00 1010 N/m3

×

The natural frequencies are ω1 = 88.6

rad/sec = 14.1 Hz

ω2 = 555.2

rad/sec = 88.4 Hz

(6.2.11)

ω3 = 1554.7 rad/sec = 247.4 Hz Example 6.2.2. Plot the mode shapes of clamped-free beam with the same dimension specified above example. Normalize the eignefunction as L

  0

A1 = 0.56461, A2 = 0.03139, A3 = 0.00133,

W i2 dx = 1

W 1 (x) = A1 [1.72[cos(5.51x)

− cosh(5.51x)] − 1.26[sin(5.51x) − sinh(5.51x)]] W 2 (x) = A2 [1.72[cos(13.81x) − cosh(13.81x)] − 1.75[sin(13.81x) − sinh(13.81x)]] W 3 (x) = A3 [1.71[cos(23.10x) − cosh(23.10x)] − 1.71[sin(23.1x) − sinh(23.1x)]]

[Homework 7] Calculate the natural frequency and plot first four mode shapes of beam with free-free boundary condition .

47

3.5

3

2.5

2

1.5

1

0.5

0.05

0.1

0.15

0.2

0.25

0.3

3

2

1

0.05

0.1

0.15

0.2

0.25

0.3

0.05

0.1

0.15

0.2

0.25

0.3

-1

-2

2

1

-1

-2

-3

Figure 6.3: First three eigenfunctions for clamped-free beam

48

6.3

Some Properties

The normal modes must satisfy the equation of motion and its boundary conditions. The normal L modes W i are also orthogonal functions satisfying the relation 0 m(x)W i (x)W j (x)dx = 0 for j = i and M i for j = i.

 



From the expansion theorem for self-adjoint distributed systems, the solution of equation (6.1.2) can be expressed as w(x, t) =

∞



W i (x)q i (t)

(6.3.1)

i=1

The generalized coordinate q i (t) can be determined from Lagrange’s equation by establishing the kinetic and potential energies.

6.4

Forced Vibration

The forced response of a beam can be calculated using modal analysis just as in the lumped system. The approach again uses the orthogonality condition of the unforced system’s eigenfunctions to reduce the calculation of the response to a system of decoupled modal equations for the time response. m

∂ 2 w ∂w ∂ 4 w + c + EI  = f (x, t) ∂t 2 ∂t ∂x 4

(6.4.1)

First, expand the applied force p(x, t) in terms of the modes f (x, t) =

∞



f i (t)W i (x)

(6.4.2)

i=1

Multiply both sides of this equation by W j and then integrate over the beam span, L

f i (t) =

 

f (x, t)W i (x)dx

(6.4.3)

0

Substituting equation (6.3.1) and (6.4.2) into equation (6.4.1) we obtain

∞





[mW i (x)¨ q i (t) + cW i q ˙(t) + EI W i (x)q i (t)] =

i=1

∞



f i (t)W i (x)

(6.4.4)

f i (t)W i (x)

(6.4.5)

i=1

We know that the modes are satisfying 

EI W i (x) = ωi2 mW i (x) Substituting this relation into equation (6.4.4) leads to

∞

 i=1

[mW i (x)¨ q i (t) + cW i q ˙(t) +

ωi2 mW i (x))q i (t)]

=

∞

 i=1

49

It is convenient to normalize the eigenfunction as L

  0

W i2 (x)dx = 1, i = 1, 2,

(6.4.6)

· ··

Since the eigenfunction W i (t) are not zero, equation (6.4.5) becomes infinite set of independent modal equations: q¨ i (t) + 2ζ i ωi q ˙(t) + ωi2 q i (t) = f i (t)/m, i = 1, 2, (6.4.7)

·· ·

6.4.1

Point force excitation

The excitation force becomes f (x, t) = F (t)δ (x

− x1)

(6.4.8)

Equation (6.4.3) becomes L

f i (t) =

 

F (t)δ (x

0

− x1)W i(x)dx = F (t)W i(x1)

(6.4.9)

Equation (6.4.7) becomes q¨ i (t) + 2ζ i ωi q ˙(t) + ωi2 q i (t) = F (t)W i (x1 )/m,

i = 1, 2,

·· ·

(6.4.10)

0

10

−1

10

−2

10

   ) −3 10    N    /   m    (    F −4    R    F10

−5

10

−6

10

−7

10

0

50

100

150

200

250

300

350

400

450

500

Frequency (Hz)

Figure 6.4: The frequency response function(FRF) obtained from the modal model. Example 6.4.1. Obtain the frequency response function of the cantilever beam. The specification of beam is the same with the Example 6.2.1. The excitation position, xa is 0.34mm and sensing position, xs is 0.2m. The example MATLAB code is printed below and the frf was printed at figure (6.4)

50



Example MATLAB Code

————————–cantifrf.m————————————clear m = 3; n = 1; z = 0.001; rho = 2750; E = 70e9; L = 0.34; b = 0.022; t = 0.002; A = t*b; Is = t ∧ 3*b/12; mass=rho*A; xa = [0.34]; xs = [0.2]; global betaL beta Ar betaL=[1.875104 4.694091 7.854757]; beta=betaL/L; Ar=[0.56461 0.031393 0.00133]; wn=[88.6 555.2 1554.7]; M = eye(m,m); K = diag(wn(1:m).∧ 2,0); Damp=diag(2*z*wn(1:m),0); Bf=zeros(m,n); ys=zeros(m,n); for i = 1:n for r = 1:m Bf(r,i)=cantimode(r,xa)/mass; end end for i = 1:n for r = 1:m ys(r,i)=cantimode(r,xs); end end A = [zeros(m,m),eye(m);-inv(M)*K,-inv(M)*Damp]; B =[zeros(m,1);inv(M)*Bf(:,1)]; C = [ys’,zeros(n,m)]; D =zeros(size(Cc,1),size(B,2)); w = linspace(0,500*2*pi,800); [mag,phs] = bode(A,B,C,D,1,w); semilogy(w/2/pi,mag(:,1)) —————————————————————————-cantimode.m———————————————function y = cantimode(r,x)

51

global betaL beta Ar y=Ar(r)*((sin(betaL(r))-sinh(betaL(r)))*(sin(beta(r)*x)-sinh(beta(r)*x)) +(cos(betaL(r))+cosh(betaL(r)))*(cos(beta(r)*x)-cosh(beta(r)*x))); —————————————————————

Homework 8 Obtain the frequency response function of the beam with free-free boundary condition. The specification of the system is the same with above example.

6.4.2

Moment excitation

The excitation force becomes f (x, t) = M o

∂  [δ (x ∂x

− x2) − δ (x − x1)]

(6.4.11)

Equation (6.4.3) becomes L

 

f i (t) =

0

M o

∂  [δ (x ∂x





− x2) − δ (x − x1)]W i(x)dx = M o [W i (x2) − W i (x1)]

(6.4.12)

Equation (6.4.7) becomes 

q¨ i (t) + 2ζ i ωi q ˙(t) + ωi2 q i (t) = M o [W i (x2 )



− W i (x1)]/m,

i = 1, 2,

· ··

(6.4.13)

Example 6.4.2. Plot the frequency response function of the cantilever beam excited by coupled moment. ρ(Density)

L(Length)

b(Width)

t(Thickness)

E

x1

x12

2750 kg/m3

340 mm

22 mm

2 mm

7.00 1010 N/m3

48 mm

80 mm



W r (x) = Ar [β r (sin β r L

×

− sinh β r L)(cos β r x − cosh β r x) − β r (cos β r L + cosh β r L)(sin β r x + sinh β r x)] (6.4.14)



Example MATLAB Code

clear m = 3; n = 1; z = 0.001; rho = 2750; E = 70e9; L = 0.34; b = 0.022;

⇒ cantipzt.m

52

t = 0.002; A = t*b; Is = t 3*b/12;

∧

mass=rho*A; xa1 = [0.048]; xa2 = [0.080]; xs = [0.340]; global betaL beta Ar betaL=[1.875104 4.694091 7.854757]; beta=betaL/L; Ar=[0.56461 0.031393 0.00133]; wn=[88.6 555.2 1554.7]; M = eye(m,m); K = diag(wn(1:m). 2,0);

∧

Damp=diag(2*z*wn(1:m),0); Bf=zeros(m,n); ys=zeros(m,n); for i = 1:n for r = 1:m Bf(r,i)=(dcantimode(r,xa2)-dcantimode(r,xa1))/mass; end end for i = 1:n for r = 1:m ys(r,i)=cantimode(r,xs); end end A = [zeros(m,m),eye(m);-inv(M)*K,-inv(M)*Damp]; B =[zeros(m,1);inv(M)*Bf(:,1)]; C = [ys’,zeros(n,m)]; D =zeros(size(C,1),size(B,2)); w = linspace(0,500*2*pi,800); [mag,phs]=bode(A,B,C,D,1,w); semilogy(w/2/pi,mag(:,1)) —————function y = dcantimode(r,x) global betaL beta Ar y1=beta(r)*(sin(betaL(r))-sinh(betaL(r)))*(cos(beta(r)*x)-cosh(beta(r)*x)); y2=-beta(r)*(cos(betaL(r))+cosh(betaL(r)))*(sin(beta(r)*x)+sinh(beta(r)*x)); y=Ar(r)*(y1+y2);

53

0

10

−1

10

−2

10

   )    T    /   m    (    F    R    F

−3

10

−4

10

−5

10

−6

10

0

50

100

150

200

250

300

350

400

450

500

Frequency (Hz)

Figure 6.5: The frequency response function(FRF) excited by PZT.

Chapter 7

Plate 7.1

Plate in Bending

The bending behavior of plates can be understood with a direct extension of what we have already learned about the bending of beams. 1 ∂ 2 w = , ρx ∂x 2

1 ∂ 2 w = ρy ∂y 2

(7.1.1)

let u and v be components of displacements at any point in the plate, u(x,y,z) =

, −z ∂w ∂x

v(x,y,z) =

−z ∂w ∂y

(7.1.2)

The strain are given as follows:

{εb } =

 

εx εy γ xy

 

=

 

∂u ∂x ∂v ∂y ∂u ∂y

+

∂v ∂x

 

=

  − z

∂ 2 w ∂x 2 ∂ 2 w ∂y 2 ∂ 2 w 2 ∂x∂y

 

(7.1.3)

Hooke’s law for plane stress relates these strains to the stress resultants, E  σx = (εx + νε y ) = 1 ν 2

−

σy =

E  (εy + νεx ) = 1 ν 2

−

τ xy =

2

Ez 1 ν 2

− −

− 1 Ez − ν 2

 

∂ 2 w ∂ 2 w + ν  ∂x 2 ∂y 2 ∂ 2 w ∂ 2 w + ν  ∂y 2 ∂x 2 2

∂  w Ez ∂  w =− −2Gz ∂x∂y 1 + ν  ∂x∂y

 

(7.1.4) (7.1.5) (7.1.6)

The stress-strain relationships take the matrix form as

{σb } = [Db]{εb} 54

(7.1.7)

55

where [Db ] =

 

E  1 ν 2

−

1

ν 

0

ν  1

0

0

1−ν  2

0

 

(7.1.8)

Each stress resultant is multiplied by its respective moment arm, yielding the following moments h/2

M x =

     

zσ x dz =

−D

zσ y dz =

−D

−h/2 h/2

M y =

−h/2 h/2

M xy =

 

∂ 2 w ∂ 2 w + ν  ∂x 2 ∂y 2 ∂ 2 w ∂ 2 w + ν  ∂y 2 ∂x 2

 

(7.1.9)

(7.1.10)

2

zτ xy dz =

−h/2

∂  w −D(1 − ν ) ∂x∂y

(7.1.11)

where D = Eh 3 /12(1 ν 2 ) is called the flexural rigidity of the plate. Equation (7.1.9)-(7.1.11) relate moments to deflection w.

−

D



∂ 4 w ∂ 4 w ∂ 4 w + 2 + ∂x 4 ∂x 2 ∂y 2 ∂y 4



+ ρh

∂ 2 w = p(x,y,t) ∂t 2

(7.1.12)

or D

7.2

∇4w + ρhw¨ = p(x,y,t)

(7.1.13)

Equation of Motion

The bending energy expression for the thin plate are 1 V  = 2 The bending energy becomes

1 V  = 2

 {

σb

v

 { } b

v

}T {b }dυ

T 

[Db ] b dυ

{ }

(7.2.1)

(7.2.2)

The kinetic energy of the plate is given by 1 T  = 2

 

ρhw˙ 2 dA

(7.2.3)

A

The response of the structure is defined in physical coordinates as a series expansion over the generalized coordinates: N 

w=

 r =1

φr (x, y)q r (t)

(7.2.4)

56

Substituting equation (7.2.4) into equation (7.2.3), one obtains an expression for the entry of the i th row and j th column of the matrix b

M s,ij = ρh

a

    0

φi (x, y)φj (x, y)dxdy

(7.2.5)

0

Substituting equation (7.2.4) into equation (7.2.2), one obtains an expression for the entry of the i th row and j th column of the stiffness matrix b

K s,ij = Ds

a

     0

where Ds =

0

h3 E s 12 1−ν s2

∂ 2 φi ∂ 2 φj ∂ 2 φi ∂ 2 φj + + ν s ∂x 2 ∂x 2 ∂y 2 ∂y 2



∂ 2 φi ∂ 2 φj ∂ 2 φi ∂ 2 φj + ∂x 2 ∂y 2 ∂y 2 ∂x 2



+ 2(1

−

∂ 2 φi ∂ 2 φj ν s ) dxdy ∂x∂y ∂x∂y (7.2.6)



Chapter 8

Approximate Method 8.1

Introduction

It is too difficult to obtain closed form solution for many problems that are more complex than a group of lumped spring-mass systems or a simple continuous system, such as a string. This section presents methods to obtain approximate solutions. With the techniques to be introduced in this chapter, we can analyze quite general systems efficiently and accurately.

8.2

Rayleigh Ritz Method

The Rayleigh Ritz method obtains an approximate solution to a differential equation with given boundary conditions using the functional of the equation. The procedure of this method can be summarized in two steps as given below: 1. Assume an admissible solution which satisfies the geometric boundary condition and contains unknown coefficients. 2. Substitute the assumed solution into the kinetic and potential energy and find the unknown coefficients. Example 8.2.1. A clamped-pinned beam with dynamic vibration absorber. we must select basis functions satisfy the boundary conditions that φ(x) = dφ(x)/dx = 0 at x = 0 and φ(x) = 0 at x = L. Hence the following functions satisfy these conditions. φr (x) =

x rπx sin( ) L L

(8.2.1)

Three-term Ritz series are considered in this example. the transverse component of beam w is

57

58

expressed as summation of Ritz function as 3

w(x, t) =



φr (x)q r (t)

(8.2.2)

r =1

and let the general displacement of  ma be denoted as q 4 . The kinetic energy and potential energy are 1 T  = 2

1 V  = 2

L

1 mb w˙ 2 dx + ma q ˙42 2

      L

0

2

∂ 2 w ∂x 2

E b I b

0

(8.2.3)

1 dx + k[w(L/2, t) 2

− q 4]2

(8.2.4)

where mb is the mass per unit length of beam.

Substituting equation (8.2.2) into equation (8.2.3) and (8.2.4), one obtains 3

3

            T  =

3

3

1 V  = E b I b 2 r=1 s=1

L

0

1 mb 2 r=1 s=1

L

0

1 φr φs dx q ˙r q ˙s + ma q ˙42 2

∂ 2 φr ∂ 2 φs 1 dx q r q s + k 2 2 ∂x ∂x 2

3

−k

  

1 φr (L/2)q r q 4 + kq 42 2 r =1

3

(8.2.5)

3



φr (L/2)φs (L/2) q r q s

r =1 s=1

(8.2.6)

The mass and stiffness matrices are

  − 

0.1413

[M ] = mb L

0.0901

0.019 0

−0.0901 0.1603

−0.0973 0

0.019

0

−0.0973

0

0.1639

0

0

µ

  

L

L/2

k ma

Figure 8.1: A clamped-pinned beam with dynamic vibration absorber

(8.2.7)

59

and [K ] =

E b I b L3

  

43.376 + 0.25κ

−74.570

75.873

368.323 −74.570 75.873 − 0.25κ −459.529 −κ/2 0

− 0.25κ −κ/2

−459.529

0

1559.3 + 0.25κ

κ/2

κ/2

κ

where µ = ma /(mb L) and κ = kL3 /E b I b . The generalized force are

  

(8.2.8)

L

Q(t) =

 

f (x, t)φr (x)dx

(8.2.9)

0

——-pinclamp.nb——Let the non-dimensional resonant angular frequency be w ¯r = when κ = 0,and ma = 0, are

 

mb L4 E b I b wr .

The resulting eigensolution

w ¯1

w ¯2

w ¯3

Ritz Method(N=3)

15.7563

50.6438

109.9433

Vibration Table

15.4182

49.9648

104.2477

The mode shapes are q 1

q 2

q 3

Mode 1

1

0.1561

-0.0021

Mode 2

0.5047

1.0000

0.1723

Mode 3

0.3058

0.6543

1.000

The natural frequency when κ = 15.75,and µa = 0.1,are

file.pc.m

w1

w2

w3

w4

Ritz Method(N=3)

11.1152

17.7539

50.7127

110.0512

The mode shapes are q 1

q 2

q 3

q 4

Mode 1

0.4282

0.0742

-0.0029

1.0000

Mode 2

1.0000

0.1457

0.0003

-0.4992

Mode 3

0.5101

1.0000

0.1731

-0.0110

Mode 4

0.3043

0.6536

1.000

0.0046

60

0.8 0.6 0.4 0.2 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.5

0

−0.5 0 1 0.5 0 −0.5 0

Figure 8.2: Mode shapes of pin-clamped beam/ pcmode1.m

Figure 8.3: Frequency response function of clamped-pinned beam/ pcmode1.m

61

0.5 0 −0.5 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.5 0 −0.5

0.5 0 −0.5

0.5 0 −0.5

Figure 8.4: Mode shapes of pin-clamped beam with dynamic absorber/ pcmode2.m

Figure 8.5: Frequency response function of clamped-pinned beam with dynamic vibration/pcfrf2.m

62

The mode shapes of beam with absorber mass are shown in figure (8.4) Example 8.2.2. Beam with free-free boundary condition. Ten-term Ritz series are considered in this example. the transverse component of tennis racket w is expressed as summation of Ritz function as 10

w(x, t) =



φr (x)q r (t)

(8.2.10)

r =1

A schematic diagram of tennis racket and ball model was shown in figure (8.6).

z

L

L

x Figure 8.6: A schematic diagram of beam with free-free boundary condition. The translate rigid mode can be expressed by φ1 (x) = 1. To represent rotational rigid mode we select the basis function to be φ2 (x) = x/L. Because there are no geometric boundary condition to satisfy, the power series can be selected as basis function.

φr (x) =

x L



r −1

,

r = 1, 2,

· ·· , 10

(8.2.11)

The kinetic energy and potential energy are 1 L T  = mb w˙ 2 dx 2 −L

 

1 L V  = E b I b 2 −L

∂ 2 w ∂x 2

(8.2.12) 2

   

dx

(8.2.13)

Substituting equation (8.2.10) into equation (8.2.12) and (8.2.13), one obtains 10

10

         

 

L 1 T  = mb φr (x)φs (x)dx q ˙r q ˙s 2 r=1 s=1 −L

1 V  = 2

10

10

L

E b I b

r=1 s=1

∂ 2 φr ∂ 2 φs dx q r q s −L ∂x 2 ∂x 2

(8.2.14)

(8.2.15)

63

The elements of the matrices can be computed as follows : M rs =

mb L r +s−1 [1

K rs = 0,

− (−1)r+s−1],

r, s = 1, 2,

· ·· , 10 K rs = 0, s = 1, 2, and r = 1, 2 · ·· , 10 s−1)(s−2) K rs = E  I  (r−L1)((rr−+2)( [1 + (−1)r+s ], r, s = 3, 4, · ·· , 10 s−5) b

r = 1, 2, and s = 1, 2

·· · , 10

(8.2.16)

b

3

When the system parameters are E b I b = 121N m, mb = 0.355kg/0.685m, 2L = 0.685m, the natural frequencies of the free-free beam are ω1

ω2

ω3

ω4

ω5

ω6

Ritz Method(N=6)(Hz)

0

0

116.9

329.3

1157.6

2361.7

Ritz Method(N=8)(Hz)

0

0

116.0

319.7

669.7

116.6

Ritz Method(N=10)(Hz)

0

0

116

320

628

104.2

Vibration Table

0

0

116

319.4

626.6

1035.9

Example 8.2.3. Tennis racket and ball Ten-term Ritz series are considered in this example. A schematic diagram of tennis racket and ball model was shown in figure (8.7). The displacements of  mB 1 and mB 2 are denoted as q 11 and q 12 .

z L

L

Ls

x kS mB kB mB

Figure 8.7: A schematic diagram of tennis racket and ball model

64

The kinetic energy and potential energy are 1 L 1 1 2 2 T  = mb w˙ 2 dx + mB 1 q ˙11 + mB 2 q ˙12 2 −L 2 2

 

1 L V  = E b I b 2 −L

∂ 2 w ∂x 2

2

   

1 dx + kS [w(LS , t) 2

(8.2.17)

− q 11]2 + 12 kB [q N +1 − q N +2]2

(8.2.18)

Substituting equation (8.2.10) into equation (8.2.17) and (8.2.18), one obtains N 

N 

           −



L 1 1 1 2 2 T  = mb φr (x)φs (x)dx q ˙r q ˙s + mB 1 q ˙N  ˙N  +1 + mB 2 q  +2 2 r=1 s=1 2 2 −L

V  =

1 2

N 

N 

L

E R I R

−L

r =1 s=1

∂ 2 φr ∂x 2

∂ 2 φs dx ∂x 2



1 q r q s + ks 2

N 

kS 

1 kB 2 2 φr (Ls )q r q N +1 + ks q N  [q  +1 + 2 2 N +1 r=1

N 

N 

 



(8.2.19)

φr (Ls )φs (Ls ) q r q s

r =1 s=1

(8.2.20)

− 2q N +1q N +2 + q N 2 +2]

The elements of the mass matrix can be computed as follows: M rs =

mb L r +s−1 [1

− (−1)r+s−1],

r, s = 1, 2,

M N +1,N +1 = mB 1

· ·· , N 

(8.2.21)

M N +2,N +2 = mB 2 The elements of the stiffness matrix are K rs,racket =

E b I b (r−1)(r−2)(s−1)(s−2) [1 L3 (r +s−5)

+ ( 1)r+s ] r, s = 3, 4,

K rs,coupled = ks φr (Ls )φs (Ls ), K r,N +1 = K N +1,r =

−

−kS φr (Ls),

r, s = 1, 2

·· · , N  r = 1, 2 · ·· , N 

· ·· , N  (8.2.22)

K N +1,N +1 = kB + kS  K N +1,N +2 = K N +2,N +1 = K N +2,N +2 = kB

−kB

EI  = 121Nm,2L = 0.685m,mB = 0.355kg/0.685m, Ls = 0.1575, mB 1 = mB 2 = 0.028kg, kS  = 4.15e4N/m, kB = 7.98e4N/m. Calculate the natural frequencies of the tennis racket and ball system. ω1

ω2

ω3

ω4

ω5

ω6

Ritz Method(N=8)(rad/s)

0

0

113.5

142.4

325.5

409.9

Ritz Method(N=10)(Hz)

0

0

114

142

325

410

Ritz Method(N=12)(rad/s)

0

0

114

142

325

410

The mode shapes are

65

Mode

q 1

q 2

q 3

q 4

q 5

q 6

q 7

q 8

q 9

q 10

q 11

q 12

Mode 1

-0.2

1

0

0

0

0

0

0

0

0

0.26

0.26

Mode 2

1

-0.2

0

0

0

0

0

0

0

0

0.91

0.91

Mode 3

0.25

-0.01

-1.

-0.07

0.3

0.03

-0.03

0.0

-0.01

0.

0.17

0.2

Mode 4

0.32

-0.04

-0.8

0.56

0.7

-0.28

-0.5

-0.02

0.15

0.05

-0.72

-1

Mode 5

0.

0.37

-0.02

-1.

0.05

0.73

-0.04

-0.27

0.01

0.04

0.02

-0.04

Mode 6

-0.05

-0.37

0.38

0.97

-0.89

-0.85

0.76

0.45

-0.24

-0.11

1

-0.76

Figure 8.8: Mode shapes of tennis racket and ball model

66

Example 8.2.4. A cantilever beam with a surface mounted active material. Calculate the natural frequencies and mode shapes by Ritz method. x rπx sin( ) L 2L

φr =

(8.2.23)

L

z x2 x1

x mb

EbIb

EpIp m p

Figure 8.9: A schematic diagram of beam with configured with a distributed piezoelectric device. The kinetic energy and potential energy are 1 T  = 2 1 V  = 2

L

L

  0

∂ 2 w ∂x 2

1 mb w˙ dx + 2 2

2

    E b I b

0

1 dx + 2

L

 

2m p w˙ 2 (H [x

L

 

− x1] − H [x − x2])dx

(8.2.24)

(H [x

(8.2.25)

0

E  p I eq

0

∂ 2 w ∂x 2

2

 

− x1] − H [x − x2])dx

where mb is the mass per unit length of beam, EI eq is the effective stiffness per unit length of the combined piezoceramic elements. Now, we expand the solution in terms of a finite set of comparison function: N 

w=



φr (x)q r (t)

(8.2.26)

r=1

Substituting equation (8.2.26) into equation (8.2.24) and (8.2.25), one obtains

1 V  = 2

N 

N 

        

1 T  = 2 N 

L

mb

0

r =1 s=1 N 

r =1 s=1

L

0

        x2

φr (x)φs (x)dx

+ 2m p

∂ 2 φr (x) ∂ 2 φs (x) E b I b dx ∂x 2 ∂x 2

φr (x)φs (x)dx

x1

x2

+ 2E  p I eq

x1



q ˙r q ˙s

(8.2.27)



∂ 2 φr (x) ∂ 2 φs (x) dx ∂x 2 ∂x 2

q r q s (8.2.28)

67

L

Msrs = mb

    0

Mprs = 2m p

φr (x)φs (x)dx

x2

φr (x)φs (x)dx

x1

(8.2.29)

L

Ksrs Kprs

∂ 2 φr (x) ∂ 2 φs (x) = E b I b dx ∂x 2 ∂x 2 0 x ∂ 2 φr (x) ∂ 2 φs (x) = 2E  p I eq dx ∂x 2 ∂x 2 x

 

 

2

1

where the equivalent area moment acting on a beam is written as b p t p3 I eq = + b p t p 12



tb t p + 2 2

2

  

mb L4 E b I b wr .

Let the non-dimensional resonant angular frequency be w ¯r = tions of cantilever beam are

(8.2.30) The resulting eigensolu-

w ¯1

w ¯2

w ¯3

w ¯4

Ritz Method(N=4)

3.5222

22.2704

67.41100

309.79199

Ritz Method(N=6)

3.5172

22.1179

62.3783

124.1380

Vibration Table

3.5160

22.0345

61.6972

120.9019

The specification of the cantilever beam are ρ(Density)

L(Length)

b(Width)

t(Thickness)

E b

2750 kg/m3

340 mm

25.4 mm

2 mm

7.00 1010 N/m3

×

The natural frequencies of the cantilever beam are

Ritz Method(N=6)

f 1

f 2

f 3

f 4

14.11 Hz

88.70 Hz

250.16 Hz

497.84 Hz

The specification of the piezo element are given as ρ(Density)

L p (Length)

b p (Width)

t p (Thickness)

8200 kg/m3

32.7 mm

22 mm

0.22 mm

E  p 1.45

×1011N/m3

x1

x2

30 mm

62.7 mm

The natural frequencies of the cantilever beam with piezoelectric materials are

Ritz Method(N=6)

f 1

f 2

f 3

f 4

16.37 Hz

93.24 Hz

250.17 Hz

495.25 Hz

68

Homework 8.2.5. A cantilever beam with lumped mass M  is shown in figure. Calculate the natural frequencies and mode shapes by Ritz method.

z L x M EbIb

mb

Figure 8.10: A schematic diagram of cantilever beam with lumped mass

Chapter 9

Finite Element Analysis

9.1 9.1.1

Euler-Bernoulli Beam Basic relation

The beam with length 2a and constant cross-sectional area A is shown in figure (??).

z

pz

a

a

x

Figure 9.1: Mode shapes of tennis racket and ball model We assume that the stress component σy , σz , τ xy , andτ yz are zero. It also assumes that the plane sections which are normal to the undeformed axis remain plane after bending. With this assumption, the axial displacement u at a distant z from the neutral axis is u(x, z) = The strain component is

−z ∂w ∂x

∂u ∂w 2 = z 2 ∂x ∂x ∂u ∂w γ xz = + =0 ∂z ∂x εx =

−

69

(9.1.1)

(9.1.2) (9.1.3)

70

The strain energy stored in the elements is given by 1 U  = 2

 

σx εx dV 

(9.1.4)

v

The normal stress is given by σx = Eε x

(9.1.5)

Substituting equation (9.1.5)and (9.1.2) into equation (9.1.4) gives , since dV  = dA dx

·

1 U  = 2 where

a

∂w 2 ∂x 2

        EI y

a

I y =

dx

(9.1.6)

z 2 dA

(9.1.7)

ρAw˙ 2 dx

(9.1.8)

A

The kinetic energy is given by

T  =

1 2

a

a

The virtual work done by external force for the element is δW ext

9.1.2

1 = 2

a

 

 pz δwdx

(9.1.9)

a

Finite Element Modeling

The displacement function can be represented by a polynomial having four constants w = α1 + α2 ξ + α3 ξ 2 + α4 ξ 3

(9.1.10)

The expression(9.1.10) can be written in the following matrix form

w = 1, ξ , ξ2  , ξ 3



Differentiating equation (9.1.10) gives aθy = a

w1 aθy 1 w2 aθy 2

α1 α2 α3 α4

  

(9.1.11)

∂w ∂w = = α2 ξ + 2α3 ξ + 3α4 ξ 2 ∂x ∂ξ 

Evaluating (9.1.10) and (9.1.12) at ξ  =

  

   

  

∓1 gives 1 −1

=

  

0

1

1

1

0

1

1

−2

−1 3

1

1

2

3

    

α1 α2 α3 α4

(9.1.12)

  

(9.1.13)

71

Solving for α gives

{}

{α} = [C ]e {w}e

(9.1.14)

{v}T e = w1, θy1, w2, θy2

(9.1.15)

where

ρAa [m]e = 105

  

EI y [k]e = 2a3

78

22a

27

22a

8a2

13a

27

13a

78

−13a −6a2 −22a

   −

3

3a

3a

4a2

3

−3a

3a

2a2

−3 −3a 3

−3a

−13a −6a2 −22a 8a2

3a 2a2

−3a 4a2

  

  

(9.1.16)

(9.1.17)

72

9.2 9.2.1

Thin Plate Theory formulation

The bending energy expression for the thin plate are 1 U  = 2

 {

σb

v

}T {b}dV 

(9.2.1)

where bending stress and strain components are

{σb} = {σx

σy

τ xy

}T ,

{εb } = {εx

εy

γ xy

}T 

(9.2.2)

The stress-strain relationships take the form

{σb } = [Db]{εb} where E  [Db ] = 1 ν 2

−

 

1

(9.2.3)

ν 

0

ν  1

0

0

0

1−ν  2

 

(9.2.4)

In deriving the energy functions for plate bending, the basic assumptions are that the direct stress in the transverse direction, σz , is zero. Also, a straight line normal to the middle surface of the undeformed plate remains normal after deformation. Therefore, the displacements u and v are given by ∂w ∂w u(x,y,z) = z , v(x,y,z) = z (9.2.5) ∂x ∂y

−

−

where w(x, y) denotes the displacement of middle surface in the z-direction. The strain are given as follows: ∂u εx κx ∂x

{εb} =

 

εy γ xy

 

=

 

∂v ∂y ∂u ∂y

+

∂v ∂x

where the curvature vector κ is given as follows:

 

=

2

{κ} = {κx

κy

∂  w κxy }T  = { 2 ∂x

 −  z

∂ 2 w ∂y 2

κy κxy

 

∂ 2 w 2 ∂x∂y

(9.2.6)

}T 

(9.2.7)

Using (9.2.6), the strain matrix can be written in the form

{εb } = −z{κ}

(9.2.8)

Substituting (9.2.8) and (9.2.3) into (9.2.1)and integrating with respect to z gives 1 U  = 2

 

A

h3 κ 12

{ }T [Db ]{κ}dA

(9.2.9)

73

The kinetic energy of the plate is given by 1 T  = 2

 

ρhw˙ 2 dA

(9.2.10)

A

The normal displacement of plate, w, and the two rotations have relation as follows: θx =

9.3

∂w , ∂y

θy =

− ∂w ∂x

(9.2.11)

Finite Element Modeling

In terms of the master coordinates (9.2.11) becomes θx =

1 ∂w , b ∂η

θy =

− a1 ∂w ∂ξ 

(9.3.1)

Since the rectangular element has 12 degrees of freedom, the displacement function can be represented by a polynomial having twelve terms, that is w = α1 + α2 ξ + α3 η + α4 ξ 2 + α5 ξη + α6 η 2 +α7 ξ 3 + α8 ξ 2 η + α9 ξη 2 + α10 η 3 + α11 ξ 3 η + α12 ξη 3

(9.3.2)

The expression(9.3.2) can be written in the following matrix form w = 1, ξ , η , ξ2  , η ξ , η 2 , ξ 3 , η ξ 2 , η 2 ξ, η3 , η ξ 3 , η 3 ξ  α



{ }

(9.3.3)

where

{α} = α1, α2, α3, ·· · , α12 Evaluating w, bθx , and aθy at ξ  = ∓1, η = ∓1gives {w¯}e = [A]e{α}

(9.3.4)

(9.3.5)

where

{w¯}T e = w1, bθx1, aθy1, ·· · , w4, bθx4, aθy4

(9.3.6)

74

and

[A]e =

         

1

−1 −1

1

0

1

0

−1

0

2

1

1

0

0

1

0

−1 −1

0

0

1

−1

0

0 0 1 0 0 1

1 0 0

1

1

−1 −2

−1 −1 −1 −1 0

1

2

3

1

1

−1 −3

−3 −2 −1 0 3 1 1 −1 1 −1 1 1 −1 1 −1 −1 −1 0 1 0 1 −2 0 1 −2 3 1 3 −1 0 −2 1 0 −3 2 −1 0 3 1

1

1

0

1

1

1

1

1

1

1

1

1

0

1

2

0

1

2

3

1

3

0

−3 −1 −1 −3

−1 −1 −3 −1

−2 −1 1 −1 0 −1 2 −1

0 1 2 0

−3 −2 −1 −1 1 −1 0 1 −2 −3 2 −1

1 3 0

Solving (9.3.5) for α gives

{}

         

{α} = [A]−e 1{w¯}e

(9.3.7)

(9.3.8)

where

  −  −        −

2 3 3

A−1 =

1 8

1

−1 −1

−1 1 1

2

1

1

2

3

1

1

3

−1 −1 −1

1

2

1

−3

−1 −1 1

1

1 3 −1 −1 −3 −1 −1 3 0 0 1 0 0 −1 0 0 −1 0 0 1 4 1 −1 −4 −1 −1 4 −1 1 −4 1 1 0 −1 0 0 −1 0 0 1 0 0 1 0 1 0 −1 −1 0 −1 −1 0 −1 1 0 −1 0 0 −1 0 0 1 0 0 −1 0 0 1 0 1 0 0 −1 0 0 1 0 0 −1 0 1 1 0 1 1 0 −1 1 0 −1 1 0 1 0 1 1 0 1 −1 0 −1 1 0 −1 −1 −1 0 1 1 0 −1 1 0 1 −1 0

         

(9.3.9)

75

The shape function in master elements are given as follows:

         

φ1 φ2 φ3 φ4 φ5 φ6 φ7 φ8 φ9 φ10 φ11 φ12

         

1 8

=

         

+ η3

2

η3

ξ + ξ 3

ξ 3

− 3 η − 3 ξ + 4 η ξ  − −η 1 − η − η 2 + η 3 − ξ + η ξ + η 2 ξ  − η 3 ξ  −1 + η + ξ  − η ξ + ξ 2 − η ξ 2 − ξ 3 + η ξ 3 2 − 3 η + η 3 + 3 ξ  − 4 η ξ + η 3 ξ  − ξ 3 + η ξ 3 1 − η − η 2 + η 3 + ξ  − η ξ  − η 2 ξ + η 3 ξ  1 − η + ξ  − η ξ  − ξ 2 + η ξ 2 − ξ 3 + η ξ 3 2 + 3 η − η 3 + 3 ξ + 4 η ξ  − η 3 ξ  − ξ 3 − η ξ 3 −1 − η + η2 + η3 − ξ  − η ξ + η2 ξ + η3 ξ  1 + η + ξ + η ξ  − ξ 2 − η ξ 2 − ξ 3 − η ξ 3 2 + 3 η − η 3 − 3 ξ  − 4 η ξ + η 3 ξ + ξ 3 + η ξ 3 −1 − η + η2 + η3 + ξ + η ξ  − η2 ξ  − η3 ξ  −1 − η + ξ + η ξ + ξ 2 + η ξ 2 − ξ 3 − η ξ 3

         

(9.3.10)

The finite element approximation of the displacement w over a given element with n nodes has the form w = N 1 (ξ, η) N 2 (ξ, η) N 3 (ξ, η) N 4 (ξ, η) w



{ }e = N (ξ, η){w}e

(9.3.11)

(9.3.12)

where

{w}T e = w1

θx1

θy 1

· ··

w4

θx4

bφ2

aφ3

· ··

φ10

bφ11

T e =

1 w˙ 2

θy 4



and N (ξ, η) = φ1



aφ12



(9.3.13)

Substituting (9.3.11)into (9.2.10) gives

{ }T e[m]e {w˙ }e

(9.3.14)

where

       

[m]e =

T N (ξ, η)dξdη

ρh N (ξ, η)

A

1

= ρhab

1

N (ξ, η)

(9.3.15)

T 

 N (ξ, η)dξdη

−1 −1

is the element inertia matrix. Substitution (9.3.11) and (9.2.7) into (9.2.9)gives U e =

1 w 2

{ }T e[k]e{w}e

(9.3.16)

76

where the stiffness element matrix is [k]e =

 

A

h3 [B]T [Db ][B]dA 12

(9.3.17)

and the curvature vector in finite element [B] is

[B] =

 

∂ 2 ∂x 2 ∂ 2 ∂y 2 ∂ 2 2 ∂x∂y

  

N (ξ, η) =



 

1 ∂ 2 a2 ∂ξ 2 1 ∂ 2 b2 ∂η 2 2 ∂ 2 ab ∂ξ∂η

  

N (ξ, η)

(9.3.18)



The expression of curvature vector with respect to the shape function is as follows

[B] =

 

1 a2 1 b2 2 ab

∂ 2 φ1 ∂ξ 2 ∂ 2 η1 ∂η 2 ∂ 2 φ1 ∂ξ∂η

2 b ∂  φ2 2 a ∂ξ 2 2 1 ∂  φ2 b ∂η 2 2 ∂ 2 φ2 a ∂ξ∂η

2 1 ∂  φ3 a ∂ξ 2 2 a ∂  φ3 b2 ∂η 2 2 ∂ 2 φ3 b ∂ξ∂η

1 a2 1 b2 2 ab

· ·· · ·· · ··

∂ 2 φ10 ∂ξ 2 ∂ 2 φ10 ∂η 2 ∂ 2 φ10 ∂ξ∂η

2 b ∂  φ11 2 a ∂ξ 2 2 1 ∂  φ11 b ∂η 2 2 ∂ 2 φ11 a ∂ξ∂η

2 1 ∂  φ12 a ∂ξ 2 2 a ∂  φ12 b2 ∂η 2 2 ∂ 2 φ12 b ∂ξ∂η

where derivatives of shape functions are as follows :

∂ 2 ∂ξ 2

          

φ1 φ2 φ3 φ4 φ5 φ6 φ7 φ8 φ9 φ10 φ11 φ12

          

 −   − −  −  − −  − −   − − − −   − − 3 ξ 

3 η ξ 

0

1

η

3 ξ + 3 η ξ 

3 ξ + 3 η ξ  0

=

1 4

1+η

3 ξ + 3 η ξ 

3 ξ 

3 η ξ 

0

1

η

3 ξ 

3 η ξ 

3 ξ + 3 η ξ  0

1+η

3 ξ 

3 η ξ 

          

 

(9.3.19)

(9.3.20)

77

∂ 2 ∂η 2

          

∂ 2 ∂ξ∂η

9.4

φ1 φ2 φ3 φ4 φ5 φ6 φ7 φ8 φ9 φ10 φ11 φ12

          

          

φ1 φ2 φ3 φ4 φ5 φ6 φ7 φ8 φ9

φ10 φ11 φ12

 −  − −    − −   − −    − − −   − −  −  − −  −  − − −  − −  −  − − −  − −  − 3η

3 η ξ 

1 + 3 η + ξ 

3 η ξ 

0

3 η + 3 η ξ 

1+3 η

=

1 4

ξ + 3 η ξ 

0

3η

3 η ξ 

1 + 3 η + ξ + 3 η ξ  0

3 η + 3 η ξ 

1+3 η

ξ 

3 η ξ 

0

          

3 η2

3 ξ 2

1+2 η

3 η2

4

1

2 ξ + 3 ξ 2

4 + 3 η 2 + 3 ξ 2 1

1 = 8

2 η + 3 η2

1 + 2 ξ + 3 ξ 2 3 η2

4

3 ξ 2

1 + 2 η + 3 η2

1

2 ξ 

3 ξ 2

4 + 3 η 2 + 3 ξ 2 2η

3 η2

1 + 2 ξ 

3 ξ 2

1

                     

(9.3.21)

(9.3.22)

Example

The natural frequency of simple supported thin plate is given by ωmn

m n = π 2 [( )2 + ( )2 ] L1 L2

where DE  = is the flexural rigidity  of the plate.

Eh 3 12(1 υ2 )

−

 

DE  ρh

(9.4.1)

(9.4.2)