Advanced Mathematics and Mechanics Applications Using MATLAB phần 5 ppsx

221: axis([uwmin,uwmax,ywmin,ywmax]);

222: axis off; hold on;

223: title(’Trace of Linearized Cable Motion’); 224:

225: % Plot successive positions 226: for j=1:ntime

227: ut=u(j,:); plot(ut,y,’-’);

228: figure(gcf); pause(.5); 229:

230: % Erase image before next one appears

231: if rubout & j < ntime, cla, end 232: end

7.2 Direct Integration Methods

Using stepwise integration methods to solve the structural dynamics equation pro￾vides an alternative to frequency analysis methods. If we invert the mass matrix and

save the result for later use, the n degree-of-freedom system can be expressed con￾cisely as a first order system in 2n unknowns for a vector z = [x; v], where v is the

time derivative of x. The system can be solved by applying the variable step-size

differential equation integrator ode45 as indicated in the following function:

function [t,x]=strdynrk(t,x0,v0,m,c,k,functim)

% [t,x]=strdynrk(t,x0,v0,m,c,k,functim)

global Mi C K F n n1 n2

Mi=inv(m); C=c; K=k; F=functim;

n=size(m,1); n1=1:n; n2=n+1:2*n;

[t,z]=ode45(@sde,t,[x0(:);v0(:)]); x=z(:,n1);

%================================

function zp=sde(t,z)

global Mi C K F n n1 n2

zp=[z(n2); Mi*(feval(F,t)-C*z(n2)-K*z(n1))];

%================================

function f=func(t)

% m=eye(3,3); k=[2,-1,0;-1,2,-1;0,-1,2];

% c=.05*k;

f=[-1;0;1]*sin(1.413*t);

In this function, the inverted mass matrix is stored in a global variable M i, the

damping and stiffness matrices are in C and K, and the forcing function name is

stored in a character string called functim. Although this approach is easy to im-

© 2003 by CRC Press LLC

plement, the resulting analysis can be very time consuming for systems involving

several hundred degrees of freedom. Variable step integrators make adjustments to

control stability and accuracy which can require very small integration steps. Con￾sequently, less sophisticated formulations employing fixed step-size are often em￾ployed in finite element programs. We will investigate two such algorithms derived

from trapezoidal integration rules [7, 113]. The two fundamental integration formu￾las [26] needed are:

b

a

f(t)dt = h

2

[f(a) + f(b)] − h3

12 f(1)

and

b

a

f(t)dt = h

2

[f(a) + f(b)] + h2

12 [f

(a) − f

(b)] + h5

720f(4)(2)

where a<i < b and h = b−a. The first formula, called the trapezoidal rule , gives

a zero truncation error term when applied to a linear function. Similarly, the second

formula, called the trapezoidal rule with end correction , has a zero final term for a

cubic integrand.

The idea is to multiply the differential equation by dt, integrate from t to (t + h),

and employ numerical integration formulas while observing that M, C, and K are

constant matrices, or

M

t+h

t

V dt ˙ + C

t+h

t

X dt ˙ + K

t+h

t

X dt =

t+h

t

P(t) dt

and

t+h

t

X dt ˙ =

t+h

t

V dt.

For brevity we utilize a notation characterized by X(t) = X0, X(t + h) = X1,

X˜ = X1 − X0. The trapezoidal rule immediately leads to

M +

h

2

C +

h2

4

K

V˜ =

t+h

t

P(t)dt − h

CV0 + K(X0 +

h

2

V0)

+ O(h3).

The last equation is a balance of impulse and momentum change involving the effec￾tive mass matrix

Me =

M +

h

2

C +

h2

4

K

which can be inverted once and used repeatedly if the step-size is not changed.

To integrate the forcing function we can use the midpoint rule [26] which states

that

b

a

P(t) dt = hP a + b

2



+ O(h3).

© 2003 by CRC Press LLC

Solving for V˜ yields

V˜ =

M +

h

2

C +

h2

4

K

−1 

P



t +

h

2



− CV0 − K



X0 +

h

2

V0



h

+ O(h3).

The velocity and position at (t + h) are then computed as

V1 = V0 + V,X ˜ 1 = X0 +

h

2

[V0 + V1] + O(h3).

A more accurate formula with truncation error of order h 5 can be developed from

the extended trapezoidal rule. This leads to

MV˜ + CX˜ + K

h

2

(X˜ + 2X0) − h2

12V˜

=

t+h

t

P(t)dt + O(h5)

and

X˜ = h

2

[V˜ + 2V0] + h2

12 [V˙

0 − V˙

1] + O(h5).

Multiplying the last equation by M and employing the differential equation to reduce

the V˙

0 − V˙

1 terms gives

MX˜ = h

2

M[V˜ + 2V0] + h2

12 [−P˜ + CV˜ + KX˜ ] + O(h5).

These results can be arranged into a single matrix equation to be solved for X˜ and

V˜ :



−(h

2M + h2

12 C) (M − h2

12 K)

(M − h2

12 K) (C + h

2 K)

 

V˜

X˜

=

hMV0 + h2

12 (P0 − P1)



P dt − hKX0

+ O(h5).

A Gauss two-point formula [26] evaluates the force integral consistent with the de￾sired error order so that

t+h

t

P(t)dt = h

2 [P(t + αh) + P(t + βh)] + O(h5)

where α = 3−√3

6 and β = 3+√3

6 .

7.2.1 Example on Cable Response by Direct Integration

Functions implementing the last two algorithms appear in the following program

which solves the previously considered cable dynamics example by direct integra￾tion. Questions of computational efficiency and numerical accuracy are examined for

two different step-sizes. Figures 7.7 and 7.8 present solution times as multiples of

the times needed for a modal response solution. The accuracy measures employed

© 2003 by CRC Press LLC

0 5 10 15 20 25 30 35 40 45 50 0

0.05

0.1

0.15

0.2

0.25

Solution Error For Implicit 2nd Order Integrator

time

solution error measure

h= 0.04, relative cputime= 34.6721

h= 0.08, relative cputime= 17.5615

Figure 7.7: Solution Error for Implicit 2nd Order Integrator

are described next. Note that the displacement response matrix has rows describ￾ing system positions at successive times. Consequently, a measure of the difference

between approximate and exact solutions is given by the vector

error_vector = \bsqrt(\bsum(((x_aprox-x_exact).à2)í));

Typically this vector has small initial components (near t = 0) and larger compo￾nents (near the final time). The error measure is compared for different integrators

and time steps in the figures. Note that the fourth order integrator is more efficient

than the second order integrator because a larger integration step can be taken with￾out excessive loss in accuracy. Using h = 0.4 for mckde4i achieved nearly the same

accuracy as that given by mckde2i with h = 0.067. However, the computation time

for mckde2i was several times as large as that for mckde4i.

In the past it has been traditional to use only second order methods for solving

the structural dynamics equation. This may have been dictated by considerations on

computer memory. Since workstations widely available today have relatively large

memories and can invert a matrix of order two hundred in about half a second, it

appears that use of high order integrators may gain in popularity.

The following computer program concludes our chapter on the solution of linear,

© 2003 by CRC Press LLC

0 5 10 15 20 25 30 35 40 45 50 0

0.05

0.1

0.15

0.2

0.25

Solution Error For Implicit 4th Order Integrator

time

solution error measure

h= 0.2, relative cputime= 13.9508

h= 0.4, relative cputime= 7.2951

Figure 7.8: Solution Error for Implicit 4th Order Integrator

constant-coefficient matrix differential equations. Then we will study, in the next

chapter, the Runge-Kutta method for integrating nonlinear problems.

© 2003 by CRC Press LLC

MATLAB Example

Program deislner

1: sfunction deislner

2: %

3: % Example: deislner

4: % ~~~~~~~~~~~~~~~~~~

5: % Solution error for simulation of cable

6: % motion using a second or a fourth order

7: % implicit integrator.

8: %

9: % This program uses implicit second or fourth

10: % order integrators to compute the dynamical

11: % response of a cable which is suspended at 12: % one end and is free at the other end. The

13: % cable is given a uniform initial velocity. 14: % A plot of the solution error is given for

15: % two cases where approximate solutions are

16: % generated using numerical integration rather 17: % than modal response which is exact.

18: %

19: % User m functions required: 20: % mckde2i, mckde4i, cablemk, udfrevib,

21: % plterror 22:

23: % Choose a model having twenty links of

24: % equal length 25:

26: fprintf(...

27: ’\nPlease wait: solution takes a while\n’) 28: clear all

29: n=20; gravty=1.; n2=1+fix(n/2); 30: masses=ones(n,1)/n; lengths=ones(n,1)/n;

31:

32: % First generate the exact solution by 33: % modal superposition

34: [m,k]=cablemk(masses,lengths,gravty);

35: c=zeros(size(m)); 36: dsp=zeros(n,1); vel=ones(n,1);

37: t0=0; tfin=50; ntim=126; h=(tfin-t0)/(ntim-1); 38:

39: % Numbers of repetitions each solution is

40: % performed to get accurate cpu times for

© 2003 by CRC Press LLC