Dynamics of viscoelastic structures: a time-domain finite element formulation

TECHNISCHE U IVERSITEIT FT LUCHTVAART· EN RU fi,TEVAARHECHNIEK

BI LlOTHEEK

Kluyverweg 1 - 2629 HS DELFT

OYNAMICS OF VISCOELASTIC STRUCTURES:

A TIME-DOMAIN FINITE ELEMENT FORMULATION

by

David Frank Golla

DYNAMICS OF VISCOELASTIC STRUCTURES: A TIME-DOMAIN FINITE ELEMENT FORMULATION

by

David Frank Golla

This dissertation is dedicated to my parents

•

ABSTRACT

A new method for incorporating viscoelastic material damping into a time-domain finite element model is developed. With the addition of new dissipation coordinates, viscoelastic element mass, damping and stiffness matrices are defined. The element matrices -- all are symmetrie and have appropriate definiteness properties -- can be assembled to synthesize motion equations for general viscoelastic structures.

Upper Case A B C 0 A 0 0+ DV E Roman

TABLE OF PRINCIPAL SYMBOLS

- plant matrix (also see definition of

I)

- input matrix (also see definition of

I)

- output matrix (also see definition of

I)

- linear viscous damping matrix

E TO E

- augmented damping matrix

- viscoelastic element damping matrix - undamped modal matrix

G - constitutive matrix

G.. - material moduli (elements of

Q)

lJ

[G], G~~ - constitutive operator tensor

~ -lJ

[G], G~~ - equilibrium elasticity tensor lJ

[G(t)], G~f(t) - relaxation function tensor

H - linear hysteretic damping matrix

H

ETH E

-[R(t)], R~;(t)

-

[G(t)]

=

[G]

+ [R(t)] k~ [H(s)], Hij(S) ~ K L _~

dissipation function tensor, [H(s)]

=

[sA(s)] - stiffness operator

- stiffness matrix

- augmented stiffness matrix

- elastic element stiffness matrix viscoelastic element stiffness matrix - linear hereditary stress-strain law

M p Q R T

u

V Z - mass mat ri x

- augmented mass matri x

- elastic element mass matrix - viscoelastic element mass matrix

- pointing matrix used to assemble element matrices _ functional for viscoelastic variational equation

-R=R A-l/2

-

-f-=-,

matrix of eigenvectors of

!,

.!!e

=

_{[R r , B.fJ} - flexible mode eigenvector of

!

- rigid body mode eigenvector of K ki net ic energy

- stress boundary condition components

transformation matrix for stresses and strains owing to rotat i on of 1 ami nate p ri nc i pal axes

- assumed mode shape for Ritz method - potential strain energy

plate constitutive matrix with partitions

z

= [ :

~

J

Lower Case Roman a b

d

_aa f h aa h(s)

_ viscoelastic dissipation pole residue - viscoelastic dissipation pole

'" - diagonal element of 0

- undamped eigenvector for mode a - generalized force vector

- diagonal element of ~ - dissipation function

.,. j r s t u x z Upper V E11' G G 12

.K

M S Lower j U~per l:!,. 1 A Q Case E₂₂ Case Case j 2

=

-1

- generalized coordinate column - position vector

- Laplace variable - time

- displacement vector; input column - state column

- output col umn

- column of dissipation coordinates Script

- global damping mat ri x - Youngls moduli

- global gyri city matri x - shear modulus

- global stiffness mat ri x - global mass mat ri x

- global spin-dissipation matri x Script

- body force distribution Greek

- prescribed displacement boundary conditions - A = b 1 ock di ag

{Q,

!t}

-

~

=

b 1 ock d i ag {Ài } - Wagnerls function - Q = block diag{w;}

Lower Case Greek

a, ~, y, ö - diss;pation pole-pair parameters - generalized modal force vector

[è] , À. 1 ~a: p [0] , w. 1 E •• lJ o .. lJ Speci al Notat i on - strain, tensor

- modal viscous damping ratio - modal coordinate column - ply rotation angle

- nonzero eigenvalue of K - Poisson's ratios

- modal loss factor - mass density - stress tensor

- modal natural frequency

- nonnegative definite

1. INTRODUCTION

Modern trends for spacecraft are towards large, 1ightweight and very f1exib1e structures. This evo1ution'from the sma11 rigid-body sate11ites of the past is creating cha11enging new prob1ems for structura1 dynamicists and control specia1ists. For sate11ites current1y under plan there is

significant interaction between structura1 vibration modes and attitude control systems. Furthermore, for some app1ications there wi11 be automatic control systems dedicated to preserving the shape of a spacecraft

structure.

These new prob1ems have focused recent attention on the important ro1e of energy dissipation in the stabi1ization and control of sate1ites. It is energy dissipation within a structure th at determines the spin axis

stabi1ity of a spinning (or part1y spinning) sate11ite and the benign

behaviour of an automatica11y controlled sate11ite depends critica11y on the level of damping in the structura1 vibration modes.

Another current trend in spacecraft structures is the use of new 1ightweight composite materials such as graphite-reinforced epoxy. These materials, c1assified as viscoe1astic, display a higher level of material damping than do metals. A1so, for harmonie osci11ations, the level of damping depends on the frequency of vibration. For some app1ications, the energy dissipation in these materials is the dominant souree of structura1 damping.

For examp1e, the dynamics of the communications sate11ite shown in Fig.

1 is significant1y affected by the f1exura1 motion of the long slender tower that conneets the antenna reflector to the 'main bus'. This unwanted motion is 1imited in amplitude by the level of damping in the tower, a 1ightweight graphite-reinforced epoxy truss.

Antenna Reflector

Solar Panel

To analyze fully the dynamics of a flexible satellite it is imperative to have an accurate technique for modeling structural damping. Vet this area of modeling has been greatly neglected in the past Whereas

sophisticated numeri cal methods (especially the finite element method) have been developed to accurately calculate mass and stiffness properties of complex structures, no comparable method has been developed for damping properties.

Energy dissipation in structures results from several microscopie processes that occur both within structural materials and at the joints of structural members. Unfortunately, no single damping process dominates for all structures. Damping processes are as varied as structures themselves.

1.1 Material Damping Processes

Material damping processes have been reviewed by Graham (Ref. 1), Bert (Ref. 2), and Lazan (Ref. 3); and damping processes may be classified in several ways. One important distinction is between anelasticity and

viscoelasticity. Anelasticity (sometimes cal led internal friction) refers to the energy dissipation processes in metals and other polycrystalline substances. Viscoelasticity refers to the nature of polymers (of ten

organic, but not necessarily so) which display properties somewhere between a fluid (viscosity) and an ideal solid (elasticity).

Anelastic energy dissipation processes may be separated into two major classes: homogeneous relaxation and inhomogeneous relaxation (Graham, Ref. 1 and Zener, Ref. 4). Homogeneous relaxation processes are:

(i) thermal diffusion (thermoelasticity) (ii) atomie diffusion

(iii) magnetic diffusion (magnetoelasticity) (iv) ordered distributions

(v) preferred distributions

Each of these processes may be further examined. For example, Graham (Ref. 1) describes three different thermoelastic processes.

Inhomogeneous relaxation processes are stress relaxation along

previously formed slip bands and dislocation relaxation effects. The latter energy dissipation processes can be attributed to the movement of

dislocations in two ways: nominal movements of dislocations about pinned positions and the break-away of dislocations from pinned positions (Ref. 1) •

Energy dissipation in viscoelastic materials results from the curling and uncurling of the long molecular chains associated with polymers (Lazan, Ref. 3). Although the primary atomic bonding in polymers is the valence bond, it is the secondary intermolecular forces (also known as van der Waals forces) that determines the viscoelastic nature of polymers (Ref. 5).

Intermolecular forces acting within polymers are dipole forces, induction forces and dispersion forces.

However, most materi al property data for po lymers is phenomeno.l ogi ca 1 (Bert, Ref. 2). Thus viscoelastic material models are not derived from the underlying physical principles of intermolecular forces. Instead, simple hypotheses and basic thermodynamical considerations lead to somewhat sophisticated mathematical models that accurately describe viscoelastic material behaviour. The most important aspect of this behaviour is the dependence of both stiffness and damping on frequency and temperature (Ref.

2).

Material damping processes mayalso be classified according to the effect of strain rate on dissipated energy. Examples of rate-independent processes are plastic strain and magnetoelasticity. Lazan (Ref. 3) calls

Rate-dependent (or dynamic hysteresis) processes are viscoelasticity and anelasticity (excluding magnetoelasticity). For small strains these

processes are linear.

Furthermore, based on this classification, anelasticity may be regarded as a subcase of viscoelasticity. As Graham (Ref. 1) observes, the

phenomenological models for linear viscoelasticity will also accurately model linear anelastic behaviour even though the underlying physical

principles are different. This is verified by Ashley's linear thermoelastic model (Ref. 6). However, a thorough understanding of the physics of

thermoelasticity would be necessary for designing damping into metallic structural members.

1.2 Energy Dissipation in Structural Joints

Energy dissipation processes in structural joints have been reviewed by Ungar (Ref. 7), Beards (Ref. 8), Goodman (Ref. 9), and Hertz and Crawley (Ref. 10). One significant factor in joint damping is lubrication. Damping in lubricated joints depends primarily on the viscosity and quantity of lubricant (Ref. 7). Dry joints, however, are essentially damped by Coulomb friction (also cal led dry friction). We shall consider here only the energy dissipation processes in dry joints.

The basic characteristics of Coulomb friction have been summarized by Beards (Ref. 8). His list is repeated here. The Coulomb friction force is

(i) dependent on the materials and their surface preparation, (ii) proportional to the normal force across the interface,

(iii) substantially independent of the sliding speed and apparent contact area,

Energy dissipation processes may be further classified depending on the magnitude of interfacial pressure. Beards (Ref. 8) defines three processes: macro-slip, micro-slip and plastic deformation.

Macro-slip occurs for low clamping pressures. A relative sliding motion between two surfaces occurs and damping results from Coulomb friction.

At higher clamping pressures micro-slip occurs. Two surfaces become mutually embedded. Only very small relative displacements of surface asperities occur.

Plastic deformation occurs'when the clamping pressure exceeds the yield pressure of the softer material. Relative motion results in further plastic deformation of the asperities. This rate-independent material damping

process becomes the dominant joint damping source.

Different types of structural joints require different analysis techniques. Two classes of joints are lap joints and slip joints. Lap joints are created when two surfaces are overlappedand fastened together with rivets, bolts, welds or adhesives (Ungar, Ref. 7). Slip joints are created when the end or edge of a structural member (such as. a rod) is

inserted into a mating piece (such as a sleeve) and held together by preload pressure, adhesives, or simply by the geometry of the structure.

Ungar (Ref. 7) has studied the energy dissipation of lap joints. for built-up beams and skin-stringer structures. The effects of fasteners

(bolts and rivets), fastener spacfng and clamping pressures are the major considerations in determining the joint damping. It was found to be difficult to apply Coulomb friction modeling, however, as it is extremely difficult to measure and control interface pressures, especially for riveted joints.

skin-stringer vibrations is an air pumping mechanisme This is not a prob1em for spacecraft structures (except perhaps for 1aunch vehic1es).

Slip joints may be c1assified further as rotationa1 slip joints or trans1ationa1 (10ngitudina1) slip joints depending on the type of 10ading and the geometry of the joint. Trans1ationa1 slip joints have been studied by Goodman (Ref. 9) and Mentel (Ref. 11). Goodman considered two cases; one in which a minimum 10ad is required to initiate slip and one in which slip begins as soon as the 10ad is app1ied. In both cases there was found to be an optimum c1amping pressure that maximized energy dissipation. Mentel further considered the effect of a viscous bonding substance between beam and sleeve, and a1so studied the slip joint damping of c1amped plates.

Hertz and Craw1ey (Ref. 10) have restricted their attention to space structure joi nts. In part i cu1 ar, they studi ed a beam-s1 eeve (translat iona1 slip) joint and a truss-pin (rotationa1 slip) joint. In the absence of a significant pre10ad c1amping pressure, the frictiona1 forces arising in the beam-sleeve joint depend on elastic displacement. Specifica11y, the normal force between surfaces varies linear1y with the amount of slip. This leads to a displacement dependent Coulomb friction model with somewhat different properties than conventiona1 Coulomb friction. With this type of damping the 10ss factor is independent of amplitude. Thus this energy dissipation process behaves 1ike anelastic material damping.

1.3 Phenomenologica1 Damping Mode1s

When viewed co11ectively, the microscopic energy dissipation processes present a formidable problem to structura1 dynamicists. Any attempt to model structura1 damping based on all the under1ying physica1 princip1es wou1d undoubted1y produce an intractable system, even when considering the

Structural dynamicists, however, are concerned with the macroscopie response of a structure. Structural damping models are therefore not

derived on the basis of microscopie energy dissipation processes. Instead, phenomenological damping models are used.

Phenomenological damping models have been reviewed hy Hughes (Ref. 12), Crandall (Ref. 13), Ottens (Ref. 14) and Baldacci (Ref. 15). The three basic phenomenological models are linear viscous damping, linear hysteretic damping and linear viscoelastic damping.

1.3.1 Linear Viscous Damping

Perhaps the most frequently used damping model is linear viscous

damping. The discretized motion equations for an undamped structure may be written

(1.1 )

Applying the modal transformation ~

=

~~ yields a set of uncoupled modal equations, ( 1.2) where ET K E

=

Q2

=

d i ag {w 2} _ _ _ - <X and

A damped model for the structure can be found by inserting modal viscous damping terms to yield

0.3)

Even if the modeled structure has no viscous damping mechanisms, this model (Eq. 1.3) will yield excellent results provided the modal damping factors (ç ) _a are found so as to match the predicted resonant amplitudes with the actual resonant amplitudes of the structure (Crandall, Ref. 13).

However, it is clear that the modal damping factors must be either measured or calculated.

Unfortunately, measured data is rarely available for a full-scale structure. Furthermore, in many instances modal damping factors are not even calculated. Hughes (Ref. 12) describes in a satirical manner the all-too-common approach in which modal damping factors are merely assigned quite arbitrarily.

(origin unknown).

The designated damping factor is usually ç

=

0.005

a

As a first step to the proper calculation of modal viscous damping

factors we assume that a global damping matrix can be calculated in the same way that Mand K are calculated. Thus the damped structural motion equation is

(1.4 )

Applying the modal transformation as before yields

( 1.5)

Generally, D is a full matrix and the modal equations are therefore coupled. Necessary and sufficient conditions for D to be diagonalized by the modal transformation (thus leaving uncoupled 'normal' modes) are known

(Refs. 16, 17); however, there are no known physical reasons why these criteria should ever be met.

Since spacecraft are lightly damped structures, however, one can almost always approximate the given system (Eq. 1.5) with a set of uncoupled modal equations (Eq. 1.3) regardless of the damping coupling (Hughes, Ref. 12). Jacobi 's formulas for eigenvalue and eigenvector perturbations show that only the diagonal elements of 0 have an important effect (Ref. 12), andone can calculate the modal viscous damping factors from

" T

2e w _{a a}

=

d _aa

=

_{-a - - a} eDe (1.6)

This conclusion is based on the assumption that the frequencies are

sufficiently separated. If the frequencies are clustered then the diagonal dominance of

Q

can be tested (Refs. 12, 15).

If required, damped modes for the coupled system canbe calculated using the method of Foss (Ref. 18). In this method an equivalent first-order system is written in the form

A

i

=

!

~ + w{t) (1. 7) where

[:

:]

[:

0

]

A

=

B

=

(1.8) K and •

[

;

]

z

=

[~

1

w

=

( 1.9)

[Here we use the notation of Hughes (Ref. 12).] In this formulation both A and ~ are symmetric matrices but neither is positive definite.

Damped modes are found from the solution of the homogeneous eigenproblem

À. Az.

=

B z.

1 - -1 - - 1

and the eigenvectors have the form

r

Àei;~i

J

~i.

=

L

(1.10)

(1.11)

Thus we have the eigenvalues (À

1·) and corresponding eigenvectors (e.) of the -1 second-order eigenproblem

À.2 Me. + À. De. + K e·

=

0

1 - -1 1 - -1 - -1 (1.12)

The eigenvalues (when not real) and the corresponding eigenvectors

appear in complex conjugate pairs. If (Ài , ~i) is a complex eigenvalue-eigenvector pair then so is (Ài *, ~i*). Note that if a mode is not coupled by the damping terms then ~i _{is real and we have (Ài ,} ~i) _{and (Ài}*, ~i) as eigenvalue-eigenvector pairs.

One important question still remains. How do we calculate the global dampi ng mat ri x?

For elastic structures there are many substructural synthesis methods (Ref. 19). Only recently, however, has damping synthesis been included with these procedures (Refs. 20, 21).

Hughes (Ref. 20) compares several substructural damping synthesis

procedures as applied to a flexible communciations satellite. Modal damping factors for the overall spacecraft dep end on both the type of substructure modes used and the assumed substructure damping models.

In every synthesis procedure, however, the burden of measuring, calculating or guessing -- and these are definitely arranged in order of preference -- modal damping factors has been shifted from the overall structure to each of the substructures. Even for the worst case (when substructure modal damping factors must be guessed) one can say (Hughes, Ref. 20) that it is better to make guesses at the substructural level than at the overall spacecraft level.

1.3.2 Linear Hysteretic Damping

Another popular damping model is linear hysteretic damping. Owing to its wide applications, this model has also been called 'structural damping ' and 'material damping ' •

Linear hysteretic damping is a frequency domain model. The uncoupled modal equations (Eq. 1.2) of an elastic structure written in the frequency domain are

[w _a 2 -_. w2J'T) _a: (jw)

=

"a(jw) ( 1.13)

The modal equations with linear hysteretic damping are

[w 2 - w2 _{+ jsgn(w)~ w 2J1) (jw)}

₌

_"a(jw)

a: a: a a (1.14)

sgn(w) = {

~

-1

w

>

0 w = 0

w

<

0

and ~ is the moda1 10ss factor. The chief characteristic of 1inear ex

hysteretic damping is a 10ss factor that is independent of frequency.

(1.15)

Severa1 incorrect time-domain representations of (1.14) have appeared in the 1 iteratu re as noted by Graham (Ref. 1) and Cranda 11 (Ref. 13). Perhaps the most common of these non-equations is

Tl _ex + w _ex 2(1 + j~)Tl _{ex ex}

1

Y _ex (t) ( 1.16)

from which the term 'complex stiffness method' has arisen to describe linear hysteretic damping. Note that it is on1y (1.16) that is objectionable, not the method itself which has been used quite successfu11y for decades by aeroelasticians (Fung, Ref. 22).

The correct time-domain representation of (1.14) is given by

~w2 (I) Tl(-t}

Tl + w 2Tl - ex ex

f

ex d't

=

Y (t)

ex ex ex 1t t-'t ex (1.17)

_ ( I )

(Cranda1l, Ref. 13, and Bronowicki, Ref. 23). Note that this model is non-causal (the upper limit of integration is (I) not t). Thus there wil 1 be

a precognitive response prior to the start of any input forces.

The motion equations for a structure with linear hysteretic damping may be written

App1ying the moda1 transformation

i

= ~

n

yie1ds

(1.19)

In general, the matrix H is fu11 and therefore the moda1 equations are coupled by the damping terms. As for the viscous damping case, however, on1y the diagona1 e1ements of H are significant and moda1 10ss factors can be def i ned by

'"

!; _a

=

h _aa /w _a 2 (1.20)

Furthermore, it can be shown (Hughes, Ref. 12) that for 1ight1y damped structures there is 1itt1e difference between 1inear viscous and 1inear hysteretic damping mode1s provided

(1.21)

Hence equivalent moda1 viscous damping factors can be defined by

r-e

=!; /2 = h /2w 2

a a aa a (1.22)

There are two methods of damping synthesis that have been used to determine moda1 10ss factors and hence equivalent moda1 damping factors. The first method is attributed to Biggs by Ottens (Ref. 14). Assuming that a single 10ss factor is known for each of m substructures then one can ca1cu1ate

m

L

_ j=1 l;a - m

L

j=1 l;. V • J Ja ( 1.23) V. Ja

where V. is the maximum stored strain energy in the jth substructure in

Ja

mode a. Thus a modal loss factor is calculated as a weighted arithmetic average of substructure loss factors. The weights are just the substructure strain energies of the (undamped) elastic mode. This method has been

successfully applied to structures incorporating viscoelastic materials (Johnson et al, Ref. 24) and is popularly called the modal strain energy method.

The calculation of loss factors for laminated composite plates using the modal strain energy method was accomplished by Morison (Ref. 25).

Experimentally determined loss factors for a single orthotropic lamina (graphite-reinforced epoxy) were used to predict loss factors for laminated plates with arbitrary ply angle configurations. Additionally, the effects of frequency dependent material properties, moisture content and long duration exposure to hard vacuum and thermal cycling were included in the damping model ing. In all cases good agreement between predicted and measured loss factors was achieved.

A method of damping synthesis using measured substructure data has been developed by Kana (Refs. 26, 27) and applied successfully to the Space Shuttle. In Kanals method damping energy (energy dissipated per cycle) is experimentally determined as a function of peak kinetic energy and amplitude for each substructure. To facilitate a comparison with the modal strain energy method we assume here a linear hysteretic damping law with damping

energy proportional to peak kinetic energy. In this case Kana's method becomes m

L

/;. T . = j=1 J Ja /;a _m (1.24 )

L

T. j=l Ja

where T. is the peak kinetic energy of the jth substructure in mode a.

Ja Note that m

L

j=1 T. = Ja m

L

j=1 V.

=

w 2/2 Ja a

However, since Tja; Vja the two methods yield different modal loss factors.

1.3.3 Linear Viscoelastic Damping

( 1.25)

Viscoelastic damping models have arisen from the recent use of

viscoelastic materials to enhance passive damping in structures. Motion equations for viscoelastic structures are most of ten written in the Laplace domain as

•

[s2M + ~(s)]g,{s) = .!Js) + M{S~ +~) (1.26)

For a simple structure in which the stiffness results from just a single material modulus this equation can be written

•

and the uncoupled modal equations become

•

(1.28)

Perhaps the most frequently used representation for h{s) is due to Biot (Ref. 28):

u

h{s)

=

L

(1.29)

1=1

With this choice for h{s) there are several solution techniques applicable to dynamic viscoelasticity.

Swanson (Ref. 29) used a numerical Laplace transform inversion technique for studying the wave propagation in a viscoelastic rode The procedure is simple to use and gives good accuracy.

Differentiation (with respect to time) can also be used (u times) to eliminate the dissipation poles [the poles of h{s)] from the motion

equations (Buhariwala, Ref. 30). Using this method the resulting equations can easily be written in the time domain with a state-space representation. Note, however, that the derivatives of the forces

(f

or Ya) must be

known.

An equivalent state-space representation that does not require differentiation is given by Hughes (Ref. 12). For a single mode the augmented state equations are

where

•

x = A x + r y

0 1 0

·

..

0 2 2 2 2 -w -Cw -w

· ..

-w A = ex ex ex ex (1.31) -a: _{0 -bI 0} al

·

..

•

.

·

•

.

•

·

·

0 a 0

·

..

-b 1) 1) and r = [0 1 0

...

O]T (1.32) -0:

Note that a stiffness proportional linear viscous damping tenn (cw _0: 2) is included in this formulation.

At this point it is interesting to point out an analogy with the modeling of unsteaqy aerodynamic forces. The dissipation function (Eq.

1.29) could also represent the unsteady aerodynamic damping forces arising from the plunging motion of an airfoil. [s2~(s) - ~(s) is the Laplace transform of Wagner's function - can be approximated by .

2

a s +

L

o 1=1

(Jonest approximation, Fung, Ref. 22)]. Similar transfer functions can be used to model general transient aerodynamic forces and gust responses

(Etkin, Refs. 31, 32).

Airplane motion equations using aerodynamic transfer functions have been formulated by Roger (Ref. 33) and Etkin (Ref. 31). The Laplace-domain

4 (s2M + sQ + !)S(s) + ~o + s~l +

L

.t=1 4 + (R

0

+

L

S!

~.t .!!..t).9. ( s) = 0 .t=1

The notation used here is essentially that of Etkin (Ref. 31). The matrices fa, fl and Q.t are matrices of constant aerodynamic influence coefficients.

Similarly, ~o and ~.t are constant aerodynamic matrices for gust response. The column ~ contains the gust velocities at panel locations.

Equivalent state-space representations can also be used for these airplane models. Etkin (Ref. 31) augments the time-domain (steady aerodynamic) motion equations with

(.t= 1, ••• ,4)

and

(.t = 1, ••• , 4)

The resulting motion equations are written in the canonical form

where

x

=

co 1 {s,

i,

J!l' ••• , ~'+' ~l' ••• , ~'+}

1.3.4 The Present Work

The aim of this dissertation is to raise the modeling of structural damping up to the same level of sophistication as mass and stiffness

analysis. In particular, we assume that energy dissipation is due to material viscoelasticity. Furthermore, the method developed here is a natural extension of the finite element method for it is fair to say that any methodology for modeling material damping that does not merge naturally with the finite element method will never be incorporated into engineering pract i ce.

A viscoelastic material is characterized by a constitutive relation in which the instantaneous stress depends on not only the instantaneous strain, but the strain history also. It is this latter dependence that provides energy dissipation. Time-domain equations of motion for a viscoelastic structure form a system of integro-differential equations which are difficult to work with numerically and for this reason viscoelastic equations of motion are almost always formulated in the Laplace domain. Unfortunately, they are of ten wanted in the time domain.

An alternate approach is presented in this paper. Time-domain realizations of Laplace-domain motion equations will be shown to yield a linear differential formulation that includes an augmented set of

generalized coordinates. In particular a second-order realization is

developed that retains the important symmetry and definiteness properties of the structural system matrices. The method will also be shown to be

completely amenable to treatment via the finite element methode Not only can the structural equations of motion be written in the time domain in the standard (second-order, linear differential) form, they can also be

synthesized by assembling certain viscoelastic element matrices. The

development of these viscoelastic element matrices is the chief contribution of this dissertation.

Moterlols ScIence

• CONTINUUM MODEL • • LAPLACE DOMAlN •

Flnlte Element Method

• DISCRETIZED MODEL • • LAPLACE DOMAlN •

General Time-Domain Realizations

llneor

System

Theory

• DISCRETIZED MODEL • • TIME DOMAlN •

Assembly Realization Principle

Viscoelostlc Structurol Dynamlcs Model

• VISCOELASTIC ELEMENT MATRICES • • ELEMENT MATRIX ASSEMBL Y PROCEDURE •

Figure 2. Outline of Viscoelastic Element Matrix Development.

and analysis techniques that lead to the development of viscoelastic element matrices. The constitutive relation for a viscoelastic material will be written as a symmetric fourth-order tensor of convolution integrals. Material properties will thus be represented by a set of relaxation functions in addition to a set of elastic constants. The most important property of the relaxation functions is that of dissipativity: Mechanical work must be done on the material for all possible strain histories.

The equations of motion of a viscoelastic structure will be formulated and discretized in the Laplace domain using a variational principle for viscoelastic continua in conjunction with the finite element method (a Ritz technique). The material properties will thus be represented by the Laplace transforms of the relaxation functions, which in turn are written as series of simple rational functions.

The key to the development of viscoelastic element matrices is in returning to the time domain with a linear differential realization. This step will require the introduction of new 'dissipation' coordinates to the equations of motion. The conventional concept of a time-domain realization will then be generalized to include second-order linear differential

systems. Furthermore, two special classes of realizations will be defined: 'symmetrical realizations' and 'linear matrix-second-order (LMSO)

realizations'. This will lead to the definition of 'viscoelastic element matrices' which retain the important symmetry and definiteness properties familiar to, and useful in, structural dynamics.

Finally, the proposed formulation will be illustrated using both longitudinal and transverse vibrations of a simple rod, and also using the transverse vibrations of an orthotropic composite-material plate. In all cases the viscoelastic element matrices are defined and assembled to form

2. EQUATION5 OF MOTION

The equation of motion governing small displacements within a viscoelastic boQy (Figure 3) can be written

where p{~) is the local mass density and ~ is a linear operator. In the case of linear elasticity ~ is cal led the stiffness operator. For the present case, linear viscoelasticity, ~ produces both conservative and dissipative forces. A prescribed set of boundary conditions (also included in the symbol ~) determines a unique solution.

The operator ~ may be factored to reveal the presence of material properties. This factorization is written here in tensor notation.

E ••

=

(u . . + u . . )/2

lJ 1,J J,l (2.2)

(2.3)

(2.4)

The spatial and time dependenee of all quantities (except p which has only spatial dependence) is implied. The boundary conditions are written

Inertial

Reference

Frame

and n is the local outward-normal unit vector defined in Slo Before

proceeding we would like to comment on two assumptions that have been made: small displacements and linear viscoelasticity.

2.1 Smal 1 Displacements

There are three different assumptions in the category of small displacements. We shall discuss these individually.

(i) Smal 1 Displacement Gradients: The Lagrangian strains for arbitrary finite deformations are written

1

El· J.

=

-2 [u. _1,J . + u. . _J,l + _k LUk . uk·] _{,1 ,J} (2.5)

The nonlinear terms in this equation are negligible if we assume that the displacement gradients are smalle Thus we assume

u . . « 1,

1 ,J all and j

(2.6)

and we use Lagrange's infinitesimal strains (Eq. 2.2). This

assumption is the only 'small-displacement' assumption necessary for the validity of the given motion equations (2.1 - 2.4). In practice, however, there are two other small-displacement assumptions of ten made. Although they reduce somewhat the class of problems that can be analysed they do not exclude problems of practical interest and they facilitate solution procedures.

(ii) Smal 1 Relative Displacements: The assumption of small relative displacements is written

for every possible set of two points (denoted by subscripts land 2) in the body. This assumption includes small displacement gradients as a limiting case when the two points are arbitrarily close to each other. Additionally, this condition implies small rigid-body

rotations. To illustrate this we look at a rigid-body rotation about an arbitrary axis (Figure 4) and choose one point on the rotation axis. If we choose the other point so that the difference in

position vectors (~2 - ~l) is perpendicular to the axis of rotation then the magnitude of the angle of rotation

{el

is found from

(2.B)

Since ~l is zero the small-relative-displacement assumption (2.7) yields

2sinl~1

«

1 (2.9)

Therefore this assumption implies small rigid-body rotations.

(iii) Smal 1 Absolute Displacements: There is no mathematically precise way to denote small absolute displacements from kinematical

considerations only. It is tempting to write

lIu 11

iri

«

I

(2.IO)

but this condition is only as good as our choice of reference frame. It is always possible to choose a reference frame so that condition

/

I I \ \

/

/

\

/

/

/

"-"

/

/

/

/

/

/

/

Figure 4. Rigid-Body Rotation.

/

/

/ I

-..

"

_"

/

/

/

/

/

/

/

/

'\.

oxis

\ \

\

I

J /

/

/

2.10 is satisfied for all bounded displacements no matter how large the bound. For example we would not consider displacements of a million kilometers for an earth orbiting satellite to be in any sense small yet they would satisfy condition 2.10 when written in a

heliocentric reference frame. Therefore the question of what is a small absolute displacement is more a matter of engineering judgment. Precise bounds on absolute displacements can, however, be determined when there is a spatially dependent external force field. Note that the right hand side of Eq. 2.1 is more correctly written as

JCe.

+ .!!' t). Thus we have assumed

(2.11)

The validity of this assumption depends on the external force field gradients and the absolute displacements. For a known external force field, small absolute displacements can be defined by

uj(E.

+ .!!' t) - jC~., t) U

«

ujCe,

t) U (2.12)

Therefore it is a spatially dependent external force field that really determines what is a small absolute displacement.

2.2 Linear Viscoelasticity

We shall assume th at the ~~f in the constitutive relation (2.3) are all linear operators. This assumption along with the others discussed above yields a linear motion equation (2.1). We can restrict our analysis in all

Having discussed the assumptions made we would like to focus attention on the constitutive relation (2.3). In the case of linear elasticity the

'operators'

g~f

are simply constants (e.g., Hooke's Law) and the constitutive relation is written

(J • •

=

lJ (2.13)

Furthermore, the elastic moduli G~; satisfy the symmetry relationships

G~ ~ = G~ ~ = G~ = Gij

lJ Jl lJ

u.

(2.14)

and the elasticity tensor is positive definite (written

[G]

>

0).

3. VISCOELASTJC MATERlAL PROPERTJES

The constitutive relation for a linear viscoelastic material is a linear hereditary stress-strain law. For elastic materials instantaneous stress is a function of instantaneous strain. For viscoelastic materials instantaneous stress is a functional of strain history. This type of relation was first formulated by Boltzmann in 1874 for a three-dimensional isotropic material. It was V. Volterra, however, who developed the

generalized equations for anisotropic viscoelastic continua in 1909. Before proceeding with the mathematical modeling of the constitutive relation, we would like to comment on an important assumption. The theory of linear viscoelasticity used throughout this dissertation is more

specifically linear isothermal viscoelasticity. Dissipated energy inevitably assumes the form of heat within a material. This in turn generates a temperature field that has, in general , both spatial and time

dependence. We shall assume that the rate of energy dissipation is

sufficiently small so that there is no significant departure from an initial uniform temperature. From this we shall infer that dissipated energy causes no significant variation in material properties.

A linear hereditary stress-strain relation may be represented most generally by a Riemann-Stieltjes convolution (Ref. 34)

(J • • (t)

lJ

(3.1)

and material properties are represented by the functions

G~J(t).

The

spatial dependence of all tensors shall be implied throughout. For further notational simplicity, time dependence will also be implied where

practicable. The properties of a linear hereditary law and the functions

G~f(t) are described in Appendix A.

To simplify the mathematics somewhat we shall assume a strain history that is continuous and piecewise continuously differentiable (with respect to time) with

- C l )

<

t <: 0 (3.2)

Under these conditions we can replace the Stieltjes convolution with a Riemann convolution and write the constitutive relation as

where G~~ = lim G~~(t) is the equilibrium elasticity and

lJ t+a> lJ

(3.4 )

[Gurtin and Sternberg (Ref. 34) have extended this result to include a step discontinuity at t

=

O. For this case Eq. 3.3 becomes

where ek~(O) is not necessarily zero. Note that this equation and Eq. 3.3

have the same Laplace transforms (remember condition 3.2).]

(i )

(i i )

The functions G~1(t) satisfy the symmetry relationships

G~~(t ) lJ

=

G~~(t) =

Jl ti< G •. (t) lJ (3.5)

The first relationships are consequences of the symmetry of the stress and strain tensors. The second relationship is derived by Biot (Ref. 28) by the application of Onsager's principle for irreversible thermodynamical

processes. LOnsager's principle stat es that the rate of entropy production is a symmetrical quadratic form of the generalized forces (Ref. 35, pp. 387-395).J In the thermodynamical derivation by Christensen (Ref. 36) this relationship is stated only for initial and final values, namely

kl

=

G .. (0»)

lJ (3.6 )

These are necessary (but not sufficient) conditions for a viscoelastic material to dissipate energy (Ref. 37). We shall use the full symmetry

relationships (Eq. 3.5) as did Biot. We shall also assume that the

functions G~~(t) all have Laplace transforms. Therefore we may write the constitutive relation in the Laplace domain

(3.7)

where

H~~(S}

=

sR~~(s}

lJ lJ

and where both LH(s)j and [R(s)] are symmetrical fourth-order tensors. For notational convenience we use here the dependenee on 's' to denote

Laplace-transformed quantities. We now represent viscoelastic material kl

properties by the equilibrium elasticity constants ~ij and the dissipation functions

H~f(S).

Perhaps the most important viscoelastic material property is

dissipativity. Real materials dissipate energy when they are stressed. In order to keep the concepts of conservative systems and dissipative systems distinct, we present the folTowing definition.

Definition 3.1: A constitutive relation shall be cal led dissipative if, for all t

>

0,

for all sufficiently smooth E . . , with

lJ

and Eij (0)

=

0

The function w(t) has the dimensions of an energy density. The total work done on the body by external forces is given by

W(t) =

IJ!

w(t)dV body

For a thermally isolated body this also represents its internal energy. The final conditions that we impose on the functions representing

viscoelastic material properties are those required to produce a dissipative constitutive relation.

Theorem 3.1: Sufficient conditions for the constitutive relation (Eq. 3.3) to be dissipative are

(i )

and (i i) [G] -( 0 (3.8)

and (i i i) [ReR(jw) ] > 0, all w

[Alternatively we could write condition (iii) as [ImH(jw)] > 0, 0 < w<

co.]

Here we use the symbol -( to mean non-negative definite. Proof

Using the constitutive relation (Eq. 3.3) we can write w(t)

=

w1(t) + I w₂(t). The first term, the equilibrium elasticity term, is written

(3.9)

Integrating once by parts and using the symmetry condition 3.8(i) yields

(3.10)

Condition 3.8(ii) clearly gives us w1(t) .. O. The second term may be written (3.11) where 0, _CD

<

't

<

0 a . . ('t) = ae:_ij O ( ' t ( t (3.12) lJ a't 0, t

<

't

<

CD

Note also that the causality property of the linear hereditary stress-strain law (see Appendix A) gives us

R~~(t) = 0,

lJ _CD

<

t

<

0

Now using Parseval 's theorem and the symmetry condition 3.8(i) yields

where superscript

*

denotes the complex conjugate. [The time dependence of the right hand side of Eq. 3.13 is 'buried ' in the definition of Cl • • (Eq.

lJ

3.12).] Condition 3.8(iii) gives w₂(t)

>

0 since the integrand is

everywhere positive. Therefore the three given conditions imply w(t)

>

O. Having discussed the essential properties of the functions used to model a viscoelastic material we turn now to the mathematical representation of these functions. The function H(s) (subscripts are dropped for

convenience) is a unique material property. However the mathematical representation of H(s) may be chosen by the analyst. Table 1 lists a representative sampling of functions that have been used to represent viscoelastic material properties.

Bagley and Torvik . (Refs. 38, 39) have shown that a fractional derivative model based on molecular theory and requiring only a few

parameters has excellent agreement with measured data. The function used by Biot (Ref. 28), on the other hand, is amenable to the classical analysis techniques of more than a century, namely, Laplace domain analysis of

rational transfer functions. Moreover, their equivalence to systems of

Zinear differential equations is invaluable.

Hence we shall use the series representation

H(s) (3.14)

to model viscoelastic materials. With this representation we accept the potential price of generating a high-order system in return for increased analytical capability. This is not without justification. We note firstly that the current trends in the computer industry 'are towards enormous

operation. Thus high-order systems will not necessarily be viewed as either unwieldy or costly in the future. Indeed, similar considerations thirty years ago would have inhibited the development of matrix structural analysis

techniques. In any case, there is already a substantial body of literature on model order reduction techniques for systems of linear differential equations and these techniques can be applied if desired.

4. THE FINITE ELEMENT METHOD

Arguably the most powerful and most popular technique for sol ving equations of motion (e.g., Eq. 2.1) is the finite element method which owes its mathematical justification to the method of Ritz. In the case of linear elasticity we use Hamilton's principle

where t2 ó

f

(T - V + W)dt

=

0 tI T _{= "2"} 1 • _{<~, pü>} V _ 1 _{- "2" ~, ,..,-}K u> w

=

<~,

i>

and where the inner product is defined as

<~,

'P

=

fff

~Tv

dV body

for any two functions u and v. Now consider the Ritz approximation (in space)

The variational principle (4.1) yields the system of equations

m1.J.

=

<U., pU.>, -1 -J f. ₁

=

<U., _-1

_6>,

k. .. = <U., K U.>(1) 1J -1 ~-J ~

=

col {qi (t)}

Here we have assumed zero stress (natural) boundary conditions.

(4.2)

(4.3)

In the special case of the finite element method the mass and stiffness matrices have factorizations that are extremely important, namely,

where

.Ma

= block diag{!:!l, M~, ••• , ~}

~

=

block diag{~l, K~, ••• , ~}

Here M:, K: are element mass and stiffness matrices and P is a pointing

-J -J

(1) This inner product cannot be directly evaluated when using finite element trial functions. The correct evaluation for this case is the subject of Appendix C.

matrix used to assemble the element matrices. It is convenient to introduce a notation to indicate the assembly of element matrices into system matrices

in a general way. Thus we write

m e _{e e Me} M = ~ M. = ~l~ ~2S

•••

cf) - j = l -J -m and (4.4) m K = ~ K~ = ~1(9 ~2 e e (9 ••• ét) Ke - j = l -J -m

The symbols @ and ~ correspond to + and Land indicate matrix assembly rather than matrix addition. The operation of matrix assembly is an important part of the original engineering interpretation of the finite element methode That this formulation is really a Ritz method (and can therefore rest on the foundations of Ritzian theory) is the true beauty of the finite element methode The form (4.4) sterns from the choice of

piecewise polynomial trial functions.

Consider now the case of linear viscoelasticity. The equation of motion in the Laplace domain is written

Alternatively we have

=

(u. . (s) + U. . (s)) /2

1 ,J J ,1 (4.6)

S2 pU .(S) - cr . . . (s) = f.(s) + p(su. +

u. )

(4.8) 1 1 J ,J 1 1 0 1 0 cr· . (s ) n.

=

T. (s) on SI lJ J 1 (4.9) u i (s) = ~i (s) on S 2 _(4.10) where "kJ. [kJ. kJ. ] [kJ. k J . ] G . . (s) = G.. + H .. (s) = G.. + sR .. (s) lJ lJ lJ lJ lJ

Consider the variational equation

60 = 0 (4.11)

where the functional Q is defined by

- Jff

{cr .. (s)e: .. (s) + [cr . . . (s) + f.(s) + p(su. + ü. )]u.(S)}dV

body 1 J 1 J 1 J ,J 1 1 0 1 0 1 (4.12) +

JJ

[cr .. (s)n. - T.]u.(s)dA S lJ J 1 1 1 +

IJ

cr .. (s )n . ~. (s ) dA S 1 J J 1 2

This functional 'parallels' that given by Christensen (Ref. 36, pp.

239-240); however, it is not simply a transformed functional and inertial terms (including initial value terms) have been added.

The first variation of the functional Q is given by

+

Jff

[S2pu.(s) - a·· .(s) - f.(s) - p(su. +

u.

)]öu.(s)dV

body 1 1 J ,J 1 1 0 1 0 1

- Jff

[e ..

(s)öa . . (s) + u.(s)öa . . . (s)]dV body lJ lJ 1 lJ,J (4.13) +

ff

[a .. (s ) n. -

T.]

at . (s ) dA 1 J J 1 1 SI +

ff

[6.(S) - u.(s)]öa .. (s)n .dA S2 1 1 lJ J +

J

J u . (s) öa .. (s )n . dA S +S 1 1 J J 1 2

Note th at the symmetry property

G~f(s)

=

G~i(s)

has been used. Note also that we have added and subtracted the term

IJ

u.(s)öa .. (s)n.dA S2 1 lJ J

Using the divergence theorem we have

JJ

u.(s)öa .. (s)n.dA

S +$ 1 1 J J

1 2

=

JIJ

[u.(s)öa . . . (s) + u . . (s)öa .. (s)]dV

body 1 lJ,J 1,J lJ

1

=

fI

f

[U . (S) öcr. . . (S) + -2 (U. . (S) + U. . (S)) öcr· . (S) ]d V

body 1 1 J ,J 1 ,J J ,1 1 J

The symmetr_{y of the stress tensor [crij (s)}

=

_{crji (s)] yields the last result.}

The first variation of Q now becomes

+

fIf

_{[s2 pui (s) - cr ... (s) - f.(s) - p(su}. +

ü.

)]öu.(s)dV

body 1 J ,J 1 1 0 1 0 1 +

fIf

[21 (u . . (s) + u . . (s)) - E • • {s)]öcr .. {s)dV body 1 ,J J ,1 1 J 1 J (4.14) +

ff

[cr .. {s)n. - T.]öu.{s)dA S lJ J 1 1 1

In this form we can see that the stationary value of Q yields the same solution as Eqs. 4.6 - 4.10. Note th at the independent variables in the variation are the spatial variables only. Note also that the first two terms in Eq. 4.12 are not the transforms of kinetic and potential energy and

Eq. 4.11 should not be regarded as a transformed principle of analytical dynamics. On the contrary Eq. 4.11 is a variational principle of a more elementary nature which quite sufficiently produces the desired solution.

The finite element method can now be used to yield discretized

equations of motion in which the system matrices are assembled from element matrices. Furthermore, the element stiffness matrices may be expanded,

revealing the independent viscoelastic moduli and their dissipation

functions. The result of this procedure is now written using the assembly notation introduced in Eq. 4.4:

[s2M + K(s)],9.(s)

=

!.(s) + M(sSa +

.90)

(i) m M=g M~ (ii) - j=l-J m e fes)

=

g

f· (s) (iii) j =1 ;) v ~(s)

= L

[G k + Hk(S)]K~ -J k=l - t \ u Hk (s)

= L

a l s/(s + _{b 1)} 1=1 (i v) (v) (4.15)

Two comrnents are in order about the notation. First, we. have abandoned the tensor notation of the previous section in Eq. 4.15(iv); all independent

viscoelastic moduli are now distinguished by a single subscript, (

)k·

Second, we have employed an algorithmic notation to further avoid multiple subscripts: all subscripted variables and the upper limits of summation dep end on all previous subscripts. For example, al' bl and u in Eq. 4.15(v) depend also on j and k.

5. TIME-DOMAIN REALIZATIONS

A realization is an important concept in linear system theory. It enables the analyst to represent a given system by an equivalent system.

In their original definition of realization, Zadeh and Desoer made no reference to the state-space representation of a system. Instead, a realization was an interconnection of components of specific classes

(integrators, differentiators, scalors and adders) equivalent to a given system (Ref. 41, p. 95). Since then, however, the standard linear

fi rst-order system of state equati ons has become the usual (a lthough not the necessary) mathematical representation of a realization. Hence the familiar definition:

Definition 5.1: The time-domain system of state equations

x=Ax+Bu

is a realization of the rational transfer matrix l(s) from ~ to ~ if

We wish now to expand the scope of the mathematical representation of a realization. This motivates the following definitions:

Definition 5.2: The time-domain system of linear differential equations

n di

L

A. - . x = B u i =0 -1 dt 1 -n-1 di 't. =

L

C. - . x i =0 -1 dt 1

-is a general realization or order n of the rational transfer matrix

I(s)

n - l . n · I

=

[L

S1C.][

L

slA·r B

. 0 -1 . 0 -1

-1= 1=

Definition 5.3: A general realization of order n shall be called

syrrunetrical i f

A! = A. ,

-1 -1 i

=

0, ••• , n

Definition 5.4: A general second-order realization given by

shall be cal led a linear matrix-second order (LMSO) realization if

VT = V 0( 0, G =-G T

T

S

=-s

We have made the substitutions

~2

=

M, ~l

=

V +

Q,

Aa

=

!

+

s

to enhance the significanee of a LMSO realization. Since any square matrix can be written as the sum of a symmetrie and a skew-symmetric matrix, we have not lost generality with the substitutions AI

=

V

+ G and Aa

=

K

+

s.

sinee we require their symmetrie parts V and ~ to be non-negative definite. Finally, note that a LMSO realization is symmetrieal if G and S are both zero.

6. THE KERNEL OF THE CONCEPT: ONE DISPLACEMENT VARIABLE, ONE ELASTIC MODULUS, AND ONE DISSIPATION-POLE PAIR IN H(s)

The primary purpose of this dissertation is to develop a time domain finite element modeling teehnique for viseoelastie struetures. To

aeeomplish this we develop a symmetrieal LMSO realization of the diseretized Laplaee-domain equation of motion for a viscoelastic structure (Eq. 4.8).

Before attempting the general case, however, it is prudent to analyse the most elementary case -- an equation with one displacement variable (or generalized coordinate in the language of Lagrange):

{S2 M+ [G + H(s)]K}q(S) = f(s) + M[sqo + qo] (6.1 ) where M

>

0, K

>

0, G

>

0 and als a2s H(s) = s + b 1 + S + b₂ with

The dissipation function H(s) mayalso be written

() as2+ys

H s

= _____

1

-S 2 + ~s + 6

where

a = al + a2,

We begin by finding a symmetrical second order realization of H-I(s) in the form

•

[We define f in the Laplace domain by f (s)

=

H(s)q(s).] Taking the

v v

Laplace transform of Eq. 6.3 and equating the transfer functions yields the following set of algebraic equations:

M - m2

=

0

(i)

M + D - 2md

=

0 (i i )

M6 + D~2 + K6 - d2~2 - 2mk6

=

a (i i i ) (6.4)

o

+ K - 2dk

=

y/~6 (i v)

K - k2

=

0 (v)

The solution of this set of algebraic equations is the subject of Appendix B. The major results are summarized here:

(i) A solution does not exist for an arbitrary choice of parameters [(al' a₂, b₁, b2) or (a, y, ~, 6)]. A solution exists on'ly when a 2

<

O.

(ii) A symmetrical LMSO realization exists only when y

=

a~ _{(a 2}

=

-a I b l/b 2)·

(iii) There are no unique solutions. One of the unknowns may always be arbitrari ly chosen.

The realizations found from these solutions will have the same input-output relationship as the given system. These realizations will therefore be zero-state equivalent to the given system. We would also like the realizations to be zero-input equivalent [Ref. 41J. Thus we have the following constraints imposed by initial values:

M - m2

₌

₀

o -

d2 = 0

m = d

and the initial values are chosen so that

Zo

=

-dqo

io

=

-dqo

( i ) (;i) (6.5) ,( i i i ) (i ) (; i ) (6.6)

Adding Eq. 6.5(ii) to our set of Eqs. (6.4) produces a set of six nonlinear algebraic equations in six unknown parameters (M, 0, K, m, d, k) and four known parameters (a, y, ~, ö). As we shall see, Eq. 6.5(iii) is consistent with their solution. Af ter a little algebra these equations reduce to

M _ m2 = _{0 ( i )} 0 _ d2 = _{0 ( i i )} K _ k2 = _{0 (i i i ) (6.7)} M + 0 - 2md = 0 (i v) M + K - 2mk

=

al ö (v)

o

+ K - 2dk = y/~ö (v i)

Eliminating M, D, K from the final three equations yields (m-d)2 = 0 (m-k)2

=

alö (d-k) 2 = y/fJö (i ) (ii) (6.8) (i i i )

and the solutions are now evident. A solution does not exist for an

arbitrary choice of parameters

(a,

y, ~,ö). A solution exists only when y

= afJ. This is precisely the same condition for a symmetrical LMSO realization! Thus all symmetrical LMSO realizations are zero-state and zero-input equivalent. The solution is not unique; one of the unknowns may be arbitrarily chosen. We select d arbitrarily and write the solution

• A

:JGJ

+ a

rl+~)

2

lJ

l+d) where A d = d/öla A Z = z/öla

[Although there are two roots óf Eq. 6.8(ii) we note that the resulting solutions share a simple equivalence relationship and are therefore not fundamentally different.]

This realization of H-l(s) yields a symmetrical LMSO realization of the full system: The given equation

where

M

>

0, K

>

0,

G

>

0

and

H ( s)

=

0:( S 2 + ~s) S 2 + ~s + Ö

has the symmetrical LMSO realization

• M+

DJ

+ D+

DJ+K+DJ=[:J

(6.9) where A M+ ₌

[:

DJ +

oK

o

Ö

[~2

:J

(6.10) A A D+ =

o:ê

K

[~2

:J

(6.11) Ö A A

[

G

:J +

oK [(l+~)

2

_(l:d)

J

K+ ₌ (6.12) 0 (1+d)

and the initial conditions are chosen so that

Zo

=

-dq ₀ (i )

•

A A .

Zo

=

-dq ₀ (i i )

(6.13)

The definiteness properties of the augmented matrices ~+,

Q+,

K+ are easily verified. These matrices are augmented because we have added a new

A

realizations. We shall call these added coordinates dissipation aoordinates.

Note that the dissipation coordinates are associated with damping. The given symmetrical LMSO realization has a simple physical

interpretation when

d

=

o.

For this case we have the simple spring, mass

and dashpot model illustrated in Figure 5. For comparison purposes we have included the model for the first-order dissipation function H(s) = asj(s+b).

The series connection of a spring and dashpot in the first-order model is known as the Maxwell model for linear viscoelasticity. The present example fixes a mass between the spring and dashpot.

Now we turn our attention to the dissipation function

H(s)

=

a(s2 +

ê

S ) s 2 + f3S + Ö

(6.14)

The restriction y

= a~

in terms of the original parameters (al' a2' bI' b2),

is

The criterion for dissipativity

ImH(jw)

>

0, O<W<CD

yields the two additional conditions

a I b I + a 2b 2 ) 0 a 1 b 2 + a 2b 1 ) 0 (i) (i i ) (6.15) (6.16)

<.TI ... aK/b aK First-Order Maxwell Model ClI3K/8 ClK/S ClK M Symmetrical LMSO Realization