Пресс-релиз популярных книг

Авторы: 111 А Б В Г Д Е Ж З И Й К Л М Н О П Р С Т У Ф Х Ц Ч Ш Щ Э Ю Я

Книги: 164 А Б В Г Д Е Ж З И Й К Л М Н О П Р С Т У Ф Х Ц Ч Ш Щ Э Ю Я

На сайте 111 авторов, 92 книг, 72 статей, 5913 глав.

2.11 Dynamics: The α-Method

Back

Instead of using modal dynamics to solve Equation (2.27), a direct solution in the space–time

domain is also feasible. Usually, the space domain is meshed with finite elements, whereas

the time domain is discretized with finite differences. The α-method, which will be presented

here because of its excellent performance, is a further development of the Newmark algorithm

(Zienkiewicz and Taylor 1989). Major contributions to the α-method were made in Hilber’s

Ph.D. Thesis (Hilber 1976) and in (Hilber et al. 1977), (Hilber and Hughes 1978), (Hulbert

and Hughes 1987) and (Miranda et al. 1989). Here, the α-method is introduced in its classical

form. For extensions of the α-method and other time-integration schemes see (Hughes 2000),

(Mu˘gan and Hulbert 2001a), (Mu˘gan and Hulbert 2001b) and (Chung et al. 2003).

Important criteria to evaluate a numerical procedure are accuracy, consistency, stability

and high-frequency dissipation. In general, second-order accuracy is strived at. This

means that the error in each iteration is O(_t2) where _t is the size of a time increment.

Accordingly, the error decreases as _t → 0, which also implies consistency. The stability

issue is generally linked to the size of _t. For some algorithms, the solution grows out of

bounds for large values of _t. This is, especially in explicit codes, a matter of concern.

If the size of _t does not matter, the algorithm is called unconditionally stable. Stability

and consistency together imply convergence. High-frequency dissipation is related to the

wish to attenuate high frequencies, which are generally less accurate due to the limited

resolution of the finite element mesh.

2.11.1 Implicit formulation

The equation to be solved is Equation (2.27):

_K_ _U_ + _M_ _ ¨U _ = _F_ . (2.323)

Finite difference discretization in time means that the finite element variables are calculated

at discrete times, for example, t = t0, t1, . . . , tn, tn+1, . . . . For simplicity, let us focus on

a material point with displacement u, velocity v and acceleration a. Integration of

a = ˙v (2.324)

yields

v(t) = vn + _ t

a(ξ ) dξ (2.325)

or, for t = tn+1,

vn+1 = vn + _ tn+1

a(ξ) dξ. (2.326)

The integral on the right-hand side of Equation (2.326) cannot be solved since a is unknown

except at t = tn and t = tn+1. However, if we approximate a(ξ ) by a linear combination

of an and an+1,

a(t) ∼ (1 − γ )an + γan+1, t∈ [tn, tn+1] (2.327)

LINEAR MECHANICAL APPLICATIONS 121

we can perform the integration leading to

vn+1 = vn + _t[(1 − γ )an + γ an+1]. (2.328)

Similarly for u

˙u

= v. (2.329)

Hence,

u(t) = un + _ t

v(η) dη. (2.330)

Substituting Equation (2.325) into Equation (2.330) yields

u(t) = un + vn(t − tn) + _ t

_ η

a(ξ ) dξ dη. (2.331)

Assuming a in the interval [tn, tn+1] to be a linear combination of an and an+1 (not

necessarily the same as in Equation (2.328), i.e. 2β _= γ in general),

a(t) = (1 − 2β)an + 2βan+1, t∈ [tn, tn+1] (2.332)

one obtains

u(t) = un + vn(t − tn) + 1

(t − tn)2[(1 − 2β) an + 2βan+1] (2.333)

or, for t = tn+1,

un+1 = un + _tvn + 1

_t2[(1 − 2β) an + 2βan+1]. (2.334)

Notice that un + _tvn is the displacement which applies if the acceleration is zero. Denoting

_V _ := _ ˙U _ and _A_ := _ ¨U _, and letting

_A_n+1

= _A_n

+ __A_ (2.335)

Equations (2.328) and (2.334) can be written for the complete mesh in the form

_V _n+1

= _V _n

+ _t _(1 − γ ) _A_n

+ γ _A_n+1

_ (2.336)

_U_n+1

= _U_n

+ _t _V _n

+ 1

(_t)2 _(1 − 2β) _A_n

+ 2β _A_n+1

_ . (2.337)

Equation (2.323) has to be satisfied at tn+1. Hence,

_M_ _A_n+1

+ _C_ _V _n+1

+ _K_ _U_n+1

= _F_ext

n+1 . (2.338)

Here, a friction term was inserted (cf Section 2.9.2), and the index “ext” stands for external

force. Substitution of Equations (2.336) and (2.337) into Equation (2.338) leads to the

Newmark algorithm. However, experience has shown that the high-frequency dissipation

122 LINEAR MECHANICAL APPLICATIONS

can be much improved if all terms in Equation (2.338) except the acceleration term are

evaluated at an intermediate position between tn and tn+1:

_M_ _A_n+1

+ (1 + α) _C_ _V _n+1

− α _C_ _V _n

+ (1 + α) _K_ _U_n+1

− α _K_ _U_n

= (1 + α) _F_ext

n+1

− α _F_ext

n , −1 ≤ α ≤ 0. (2.339)

Since the stiffness term and friction term can also be considered as an internal force:

_F_int

n+1 := _C_ _V _n+1

+ _K_ _U_n+1 (2.340)

Equation (2.339) amounts to

_M_ _A_n+1

+ (1 + α) _F_int

n+1

− α _F_int

= (1 + α) _F_ext

n+1

− α _F_ext

n . (2.341)

Substitution of Equations (2.336) and (2.337) into Equation (2.340) yields the α-method

(called after the parameter α). In the next sections, it will be proven that if β and γ satisfy

β = 1

(1 − α)2 (2.342)

γ = 1

− α (2.343)

the algorithm is second-order accurate and unconditionally stable for α ∈ [−1/3, 0]. Maximum

high-frequency dissipation is obtained for α = −1/3. For α = 0, there is no highfrequency

dissipation and the α-method reduces to a special case of the Newmark method

(also called the average acceleration method since γ = 1/2 and β = 1/4 corresponds to

taking the mean of _A_n and _A_n+1). Defining

_ ˜V _n+1

= _V _n

+ (1 − γ )_t _A_n (2.344)

_ ˜U _n+1

= _U_n

+ _t _V _n

+ 1

(_t)2(1 − 2β) _A_n (2.345)

Equations (2.336) and (2.337) yield

_V _n+1

= _ ˜V _n+1

+ γ_t _A_n+1 (2.346)

_U_n+1

= _ ˜U _n+1

+ β(_t)2 _A_n+1 . (2.347)

The quantities _ ˜U _n+1 and _ ˜V _n+1 can be considered as predictor values and depend on

values at time tn only. Substitution of Equations (2.346) and (2.347) into Equation (2.339)

yields

__M_ + (1 + α) _C_ _tγ + (1 + α) _K_ (_t)2β_ _A_n+1

+ (1 + α) _C_ _ ˜V _n+1

− α _C_ _V _n

+ (1 + α) _K_ _ ˜U _n+1

− α _K_ _U_n

= (1 + α) _F_ext

n+1

− α _F_ext

n (2.348)

LINEAR MECHANICAL APPLICATIONS 123

which is equivalent to

∗_ _A_n+1

= _F_ (2.349)

where

∗_ = _M_ + (1 + α) _C_γ_t + (1 + α) _K_ β(_t)2 (2.350)

_F_ = (1 + α) _F_ext

n+1

− α _F_ext

− (1 + α) _ ˜ F_

int

n+1

+ α _F_int

n (2.351)

_F_int

= _C_ _V _n

+ _K_ _U_n (2.352)

_ ˜ F_

int

n+1

= _C_ _ ˜ V _n+1

+ _K_ _ ˜U _n+1 . (2.353)

After solving for _A_n+1, _V _n+1 and _U_n+1 can be determined using Equations (2.346)

and (2.347).

2.11.2 Extension to nonlinear applications

The procedure can also be extended to nonlinear problems, in which _C_ and _K_ are

nonlinear functions of _U_ and _V _. Now, Equation (2.339) is replaced by

_M_ _A_n+1

+ (1 + α) _C(_U_n+1 , _V _n+1)_ _V _n+1

− α _C(_U_n , _V _n)_ _V _n

+ (1 + α) _K(_U_n+1 , _V _n+1)_ _U_n+1

− α _K(_U_n , _V _n)_ _U_n

= (1 + α) _F_ext

n+1

− α _F_ext

. (2.354)

Defining

_K_n := _K(_U_n , _V _n)_ (2.355)

_C_n := _C(_U_n

, _V _n

)_ (2.356)

Equations (2.349) and (2.351) still apply but Equations (2.350) and (2.352) to Equation

(2.353) now yield

∗_ = _M_ + (1 + α) _C_n+1 γ_t + (1 + α) _K_n+1 β(_t)2 (2.357)

_F_int

= _C_n _V _n

+ _K_n _U_n (2.358)

_ ˜ F_

int

n+1

= _C_n+1 _ ˜ V _n+1

+ _K_n+1 _ ˜U _n+1 . (2.359)

_K_n+1 and _C_n+1 are, however, not a priori known. Therefore, they are calculated on the

basis of _ ˜U _n+1 and _ ˜V _n+1 (the predictor values) and iterations are run till convergence.

Denoting the predictor values by

_V _(1)

n+1

= _ ˜ V _n+1

= _V _n

+ (1 − γ )_t _A_n (2.360)

_U_(1)

n+1

= _ ˜U _n+1

= _U_n

+ _t _V _n

+ 1

(_t)2(1 − 2β) _A_n (2.361)

124 LINEAR MECHANICAL APPLICATIONS

the force in iteration 1 amounts to

_F_(1) = (1 + α) _F_ext

n+1

− α _F_ext

− (1 + α) _F_(1)int

n+1

+ α _F_int

n (2.362)

where

_F_int

= _C_n _V _n

+ _K_n _U_n (2.363)

_F_(1)int

n+1

= _C_(1)

n+1 _V _(1)

n+1

+ _K_(1)

n+1 _U_(1)

n+1 . (2.364)

In Equation (2.364) use was made of the following abbreviations:

_C_(i)

n+1 := _C(_U_(i)

n+1 , _V _(i)

n+1)_ (2.365)

and

_K_(i)

n+1 := _K(_U_(i)

n+1 , _V _(i)

n+1)_ (2.366)

Equation (2.349) now amounts to

∗_(1) _A_(2)

n+1

= _F_(1) (2.367)

where

∗_(1) = _M_ + (1 + α) _C_(1)

n+1 γ_t + (1 + α) _K_(1)

n+1 β(_t)2. (2.368)

The corrected velocity and displacement satisfy

_V _(2)

n+1

= _V _(1)

n+1

+ γ_t _A_(2)

n+1 (2.369)

_U_(2)

n+1

= _U_(1)

n+1

+ β(_t)2 _A_(2)

n+1 . (2.370)

Now one can use the corrected displacements and velocities to calculate new values for

∗_ and _F_ and iterate until _A_n+1 has converged. For iteration i, one obtains

_F_(i) = (1 + α) _F_ext

n+1

− α _F_ext

− (1 + α) _F_(i)int

n+1

+ α _F_int

n (2.371)

_F_(i)int

n+1

= _C_(i)

n+1 _V _(1)

n+1

+ _K_(i)

n+1 _U_(1)

n+1 (2.372)

∗_(i) = _M_ + (1 + α) _C_(i)

n+1 γ_t + (1 + α) _K_(i)

n+1 β(_t)2 (2.373)

∗_(i) _A_(i+1)

n+1

= _F_(i) (2.374)

_V _(i+1)

n+1

= _V _(1)

n+1

+ γ_t _A_(i+1)

n+1 (2.375)

_U_(i+1)

n+1

= _U_(1)

n+1

+ β(_t)2 _A_(i+1)

n+1 . (2.376)

This scheme has the disadvantage that _V _(1)

n+1 and _U_(1)

n+1 have to be stored. This can

be avoided by writing Equations (2.375) and (2.376) for i and solving for _U_(1)

n+1 and

_V _(1)

n+1 :

_V _(1)

n+1

= _V _(i)

n+1

− _tγ _A_(i)

n+1 (2.377)

_U_(1)

n+1

= _U_(i)

n+1

− (_t)2β _A_(i)

n+1 . (2.378)

LINEAR MECHANICAL APPLICATIONS 125

Substituting these equations into Equations (2.371) to (2.374) one obtains

,_M_ + (1 + α) _C_(i)

n+1 γ_t + (1 + α) _K_(i)

n+1 β(_t)2- _A_(i+1)

n+1

− ,(1 + α) _C_(i)

n+1 γ_t + (1 + α) _K_(i)

n+1 β(_t)2- _A_(i)

n+1

= (1 + α) _F_ext

n+1

− α _F_ext

+ α _F_int

− (1 + α) ,_C_(i)

n+1 _V _(i)

n+1

+ _K_(i)

n+1 _U_(i)

n+1

- . (2.379)

By defining

__A_(i) = _A_(i+1)

n+1

− _A_(i)

n+1 (2.380)

_A_(1)

n+1

= _0_ (2.381)

and

_R_(i) = _F_(i)∗ − _M_ _A_(i)

n+1 (2.382)

where

_F_(i)∗ = (1 + α) _F_ext

n+1

− α _F_ext

− (1 + α) _F_(i)int∗

n+1

+ α _F_int

n (2.383)

and

_F_(i)int∗

n+1

= _C_(i)

n+1 _V _(i)

n+1

+ _K_(i)

n+1 _U_(i)

n+1 (2.384)

Equation (2.379) can be reduced to

∗_(i) __A_(i) = _R_(i)

. (2.385)

Furthermore, Equations (2.375) to (2.378) yield

_V _(i+1)

n+1

= _V _(i)

n+1

+ γ_t __A_(i) (2.386)

_U_(i+1)

n+1

= _U_(i)

n+1

+ β(_t)2 __A_(i)

. (2.387)

_R_(i) is the residual in iteration i and its size can be used as a criterion to exit the loop.

Summarizing,

1. Calculate the predictor values: Equations (2.360), (2.361)

2. Loop i, i = 1, 2, . . . : multicorrector step

(a) Calculate _F_(i)∗

(Equation (2.383)), _M

∗_(i) (Equation (2.373)), _R_(i)

(Equation (2.382)). If _R_(i) < _: exit.

(b) Calculate __A_(i) (Equation (2.385).

126 LINEAR MECHANICAL APPLICATIONS

Those degrees of freedom for which boundary values are given are not a part of the

solution system. For these degrees of freedom, either _U_n+1, _V _n+1 or _A_n+1 is given.

The predictor step yields _U_(1)

n+1, _V _(1)

n+1 and _A_(1)

n+1 (Equations (2.360), (2.361) and

(2.381)) as usual. However, now only one corrector step is necessary. We set

_U_(2)

n+1

= _U_n+1 (2.388)

_V _(2)

n+1

= _V _n+1 (2.389)

_A_(2)

n+1

= _A_n+1 . (2.390)

Then, Equations (2.387), (2.386) and (2.380) yield for i = 1

__A_(1) = 1

β(_t)2

__U_n+1

− _U_(1)

n+1

_ (2.391)

__A_(1) = 1

γ_t

__V _n+1

− _V _(1)

n+1

_ (2.392)

__A_(1) = _A_n+1 . (2.393)

Substituting __A_(1) into Equation (2.386) or Equation (2.387) yields the remaining

unknowns (displacement if the velocity is given, velocity if the displacement is given and

displacement and velocity if the acceleration is given).

The assumption that a in [tn, tn+1] is a linear combination of an and an+1 works

well unless a changes discontinuously. For instance, if the external force jumps at t = t

n ,

the acceleration a

n has to be adjusted accordingly to get accurate results. The correct

acceleration is obtained by using Equation (2.338):

_M_ _A_n+

= _F_ext

− _F_int

= _F_ext

+ __F_ − _F_int

n (2.394)

where __F_ is the force jump. Consequently, the acceleration jump amounts to

_M_ ._A_n+

− _A_n

/ = __F_ . (2.395)

2.11.3 Consistency and accuracy of the implicit formulation

In order to examine the consistency, accuracy and stability of the implicit scheme,

Equation (2.349) together with Equations (2.344) to (2.345) and Equations (2.350) to

(2.353) are written for a homogeneous (no external force) single degree of freedom

LINEAR MECHANICAL APPLICATIONS 127

system:

[m + (1 + α)cγ_t + (1 + α)kβ(_t)2]an+1

= αcvn + αkun − (1 + α)c[vn + (1 − γ )_tan]

− (1 + α)k[un + _tvn + 1

(_t)2(1 − 2β)an] (2.396)

[m + (1 + α)cγ_t + (1 + α)kβ(_t)2]an+1

= −(1 + α)[c(1 − γ )_t + 1

k(_t)2(1 − 2β)]an

− [c + (1 + α)k_t ]vn − kun. (2.397)

Defining the frequency _ and the friction coefficient ξ by

_2 := k

(_t)2 (2.398)

ξ := c

2m_

_t (2.399)

which is equivalent to

k = m_2

(_t)2 (2.400)

c = 2m_ξ

(2.401)

one can replace k and c in Equation (2.397), yielding

[1 + (1 + α)2_ξγ + (1 + α)_2β]an+1(_t)2

= −(1 + α)[(1 − γ )2_ξ + 12

(1 − 2β)_2]an(_t)2

− [2_ξ + (1 + α)_2]vn_t − _2un. (2.402)

The one-dimensional equivalent of Equations (2.344) to (2.347) yields

_tvn+1 = _tvn + (1 − γ )(_t)2an + γ (_t)2an+1 (2.403)

un+1 = un + _tvn + 1

(1 − 2β)(_t)2an + β(_t)2an+1. (2.404)

After substitution of Equation (2.402) into the right-hand side of Equations (2.403) and

(2.404), these three equations form a system expressing un+1, _tvn+1 and (_t)2an+1 in

terms of un, _tvn and (_t)2an:





un+1

(_t)vn+1

(_t)2an+1





= _A_





(_t)vn

(_t)2an





(2.405)

128 LINEAR MECHANICAL APPLICATIONS

where

_A_ =





1 + βP 1 + βQ 12

+ β(R − 1)

γP 1 + γQ 1 + γ (R − 1)

P Q R





(2.406)

and

P = −_2/d (2.407)

Q = −[2ξ_ + (1 + α)_2]/d (2.408)

R = −(1 + α)[(1 − γ )2_ξ + 12

(1 − 2β)_2]/d (2.409)

d = 1 + (1 + α)2_ξγ + (1 + α)_2β. (2.410)

Equation (2.405) is a set of three homogeneous finite difference equations of the first order

with constant coefficients. Solutions are obtained by setting





(_t)vn

(_t)2an





= λn _X_ (2.411)

leading to

_A_ _X_ = λ _X_ (2.412)

which is a classical eigenvalue problem. Solutions exist if the characteristic equation is

satisfied

λ3 − I1λ2 + I2λ − I3 = 0 (2.413)

where I1, I2 and I3 are the invariants of A (the elements of A are denoted by A11, . . .) :

I1 = tr A = A11 + A22 + A33 (2.414)

I2 =_

A22 A23

A32 A33

____

A11 A13

A31 A33

____

A11 A12

A21 A22

____

(2.415)

I3 = det A. (2.416)

Equation (2.413) is also the characteristic equation of the following third-order equation:

D = un+1 − I1un + I2un−1 − I3un−2 = 0 (2.417)

which is equivalent to Equation (2.405). After some algebra, one obtains the invariants of

Equation (2.406):

I1 = 2 − ,2ξ__1 + α(1 − γ )_ + _2 _(1 + α)(γ + 12

) − βα_- /d (2.418)

I2 = 1 − ,2ξ__1 + 2α(1 − γ )_ + _2 _γ − 1/2 + 2α(γ − β)_- /d (2.419)

I3 = −_2ξ_α(1 − γ ) + _2α(γ − β − 12

_ /d. (2.420)

LINEAR MECHANICAL APPLICATIONS 129

To examine the consistency and accuracy, Equation (2.417) is expanded by Taylor series

about un = u:

un+1 = un + ˙u_t + ¨u

(_t)2

+ ...

(_t)3

+· · · (2.421)

and similarly for un−1 and un−2. Collecting terms, one obtains

D = (1 − I1 + I2 − I3)u + (1 − I2 + 2I3) ˙ u_t

+ (1 + I2 − 4I3)¨u

(_t)2

+ (1 − I2 + 8I3)

...

(_t)3

+· · · (2.422)

Defining

I i := Ii(d = 1) (2.423)

one obtains after some algebra

1 − I 1 + I 2 − I 3 = _2 (2.424)

1 − I 2 + 2I 3 = _

2 [4ξ + (−1 + 2α + 2γ )_] (2.425)

1 + I 2 − 4I 3 = 12

,4 − 4[1 + 2α(γ − 1)]ξ_ + (1 − 4α − 4αβ − 2γ + 4αγ )_2-

(2.426)

1 − I 2 + 8I 3 = −12

_{−4[1 + 6α(γ − 1)]ξ + [1 − 4α(2 + 3β − 3γ ) − 2γ ]_} (2.427)

and

1 − I1 + I2 − I3 = (1 − I 1 + I 2 − I 3)/d (2.428)

1 − I2 + 2I3 = (1 − I 2 + 2I 3)/d (2.429)

1 + I2 − 4I3 = [2(d − 1) + 1 + I 2 − 4I 3]/d (2.430)

1 − I2 + 8I3 = (1 − I 2 + 8I 3)/d. (2.431)

Now, we will assume there is no friction: ξ = 0. Defining ω by

_ := ω_t (2.432)

collecting terms in Equation (2.422) yields

D = (_t)2(ω2u+ ¨u) + (_t)3 _12

(−1 + 2α + 2γ )ω2˙u_ + O(_t)4 (2.433)

¨u

+ ω2u = D

(_t)2

− _t _(α + γ − 12

)ω2˙u_ + O(_t)2. (2.434)

The left-hand side is the governing equation without friction, the first term on the righthand

side is its approximation. For _t →0 both are identical. Consequently, the numerical

scheme is consistent. The accuracy is given by the power of the subsequent _t terms on

the right-hand side. In general, the scheme is first order. However, if

γ = 12

− α (2.435)

the scheme is second order. This will be assumed to be the case in the following section.

130 LINEAR MECHANICAL APPLICATIONS

2.11.4 Stability of the implicit scheme

A numerical scheme is stable if its amplification matrix _A_, Equation (2.406), has no

eigenvalues with size greater than one. Indeed, if |λ| > 1, the homogeneous solution,

Equation (2.411) diverges. Since a particular solution (i.e. a solution of the inhomogeneous

problem) augmented by a linear combination of the homogeneous solutions is a

particular solution as well, any instable homogeneous solution will lead to an instability

of the inhomogeneous scheme. Accordingly, we have to solve Equation (2.413) for λ and

check that |λ| ≤ 1.

Since the coefficients I1, I2 and I3 are a function of _, Equations (2.418) to (2.420),

so is its solution λ. Recall that the parameter γ is defined by Equation (2.435), α and β

are still adjustable. The strategy that will be adopted here follows the treatise in (Hilber

1976): first, we check the value of |λ| for _ = 0, then we examine _→∞ and finally

look at the values in-between.

1. For _ = 0 we have I1 = 2, I2 = 1 and I3 = 0. The roots are λ = 0 and λ = 1, the

latter is a double root. The condition |λ| ≤ 1 is satisfied.

2. For _→∞I1, I2 and I3 reduce to

I1 ∼ 2 − (1 + α)(γ + 12

) − βα

(1 + α)β

, _→∞ (2.436)

I2 ∼ 1 − γ − 12

+ 2α(γ − β)

(1 + α)β

, _→∞ (2.437)

I3 ∼ α(γ − β − 12

)

(1 + α)β

, _→∞. (2.438)

Notice that ξ drops out for _→∞. Substituting γ = 12

− α and rearranging,

Equation (2.413) yields

α(α + β) − (2α2 + β(1 + α) + 2αβ)λ

+ [2β(1 + α) + αβ − 1 + α2]λ2 − λ3(1 + α)β = 0. (2.439)

Since the term (1 + α)/β that is multiplying λ3 originates from the denominator

of Equations (2.436) to (2.438), α = −1 and β = 0 must be excluded. The solution

of Equation (2.439) is (by inspection or by using a symbolic mathematical program)

λ3 = α

1 + α

(2.440)

λ1,2 = 1 − 1 − α

2β

± 1

2β

_(1 − α)2 − 4β. (2.441)

Let us first have a closer look at λ1 and λ2, Figures 2.25 and 2.26.

LINEAR MECHANICAL APPLICATIONS 131

−2

−1

|λ1| < 1 |λ1| = 1 |λ1| > 1

λ1 > 1

λ1→∞

λ1 →−∞

λ1 < −1

λ1 = −1

−1 < λ1 < 1

β = (1 − α)2/4

β = (1 − 2α)/4

Figure 2.25 Evaluation of λ1

(a) For (1 − α)2 < 4β, λ1 and λ2 are complex and

|λ1| = |λ2| =

0α + β

. (2.442)

|λ1| = |λ2| = 1 for α = 0, to the left of the β-axis we have |λ1| < 1, |λ2| < 1,

to the right we obtain |λ1| > 1, |λ2| > 1.

(b) On the parabola (1 − α)2 = 4β and one finds

λ1 = λ2 = −1 + α

1 − α

(2.443)

leading to

α < 0⇒−1 < λ1,2 < 1

α = 0 ⇒ λ1,2 = −1

132 LINEAR MECHANICAL APPLICATIONS

−2

−1

|λ2| < 1 |λ2| = 1 |λ2| > 1

λ2 > 1

λ2→∞

λ2 →−∞

λ2 < −1

λ2 = −1

−1 < λ2 < 1

β = (1 − α)2/4

β = (1 − 2α)/4

Figure 2.26 Evaluation of λ2

0 < α < 1 ⇒ λ1,2 < −1

α > 1 ⇒ λ1,2 > 1. (2.444)

λ1 _= 1, λ2 _= 1 (2.445)

λ1 = −1 ⇔ β = 1 − 2α

, α≥ 0 (2.446)

λ2 = −1 ⇔ β = 1 − 2α

, α≤ 0. (2.447)

The straight line β = (1 − 2α)/4 is tangent to the parabola at α = 0. For β →0,

λ1,2 is not properly defined by Equation (2.441) and an asymptotic expansion

must be developed

λ1,2 = 1 − 1 − α

2β

|1 − α|

2β

1 − 4β

(1 − α)2

LINEAR MECHANICAL APPLICATIONS 133

= 1 − 1 − α

2β

|1 − α|

2β

_1 − 2β

(1 − α)2

+ O(β2)

, β→0

= 1 − 1 − α

2β

|1 − α|

2β

∓

|1 − α|

(1 − α)2

+ O(β), β → 0. (2.448)

This yields for λ1

λ1 = 1 − 1

1 − α

+ O(β), β → 0, α < 1

λ1 = α − 1

+ O(1), β → 0, α > 1

(2.449)

and for λ2

λ2 = −1 − α

+ O(1), β → 0, α < 1

λ2 = 1 + 1

α − 1

+ O(β), β → 0, α > 1

. (2.450)

Summarizing, the straight line β = (1 − 2α)/4 divides the region under the parabola

into three zones. Only for α ≤ 0 and β ≥ (1 − 2α)/4 we have |λ1| ≤ 1 and |λ2| ≤ 1.

λ1 and λ2 are sometimes called the principal roots.

For the third eigenvalue (sometimes called the spurious root), one obtains

|λ3| = 1 ⇔ λ3 = −1 ⇔ α = −1

. (2.451)

However (see Figure 2.27),

|λ3| ≤ 1 ⇔ α ≥ −1

|λ3| > 1 ⇔ α < −1

(2.452)

Accordingly, the implicit scheme is unconditionally stable (i.e. |λ1| ≤ 1, |λ2| ≤ 1

and |λ3| ≤ 1) at high frequencies only if

−1

≤ α ≤ 0

β ≥ 1 − 2α

(2.453)

Now, for _→∞, we would like to maximize the dissipation to get rid of spurious

high-frequency effects, that is, we seek to minimize max(|λ1|, |λ2|, |λ3|). How does

β affect max(|λ1|, |λ2|, |λ3|)? Since λ3 is not a function of β, we focus on λ1 and

λ2. One can prove that for fixed α ≥ −12

, max(|λ1|, |λ2|) attains a minimum for

134 LINEAR MECHANICAL APPLICATIONS

−2

−1.5

−1

−0.5

0.5

1.5

2.5

λ3(−)

α(−)

Figure 2.27 Evaluation of λ3

β = (1 − α)2/4, that is, on the parabola, satisfying

|λ1| = |λ2| = 1 + α

1 − α

. (2.454)

Indeed, for β > (1 − α)2/4 we deduced Equation (2.442) which is a monotonic

increasing function of β toward 1 for β →∞. The derivative of Equation (2.441)

with respect to β satisfies

∂λ1,2

∂β

= 1 − α

2β2

∓ 1

2β2

_(1 − α)2 + 4β ∓ 1

β_(1 − α)2 + 4β

. (2.455)

For β = (1 − α)2/4 one gets

λ1 = λ2 = −1 + α

1 − α

<0 forα > −1. (2.456)

Now, ∂λ2

∂β > 0 (Equation (2.455)), hence, λ2 decreases with decreasing β. Accordingly,

|λ2| and max(|λ1|, |λ2|) increase with decreasing β in a neighborhood of

β = (1 − α)2/4. This completes the proof.

Consequently, we maximize the high-frequency dissipation if we take

β = (1 − α)2

(2.457)

LINEAR MECHANICAL APPLICATIONS 135

and the only parameter left is α. Summarizing, for −12

≤ α ≤ 0:

|λ1| = |λ2| = 1 + α

1 − α

≤ 1 (2.458)

|λ3| =

−α

1 + α

≤ 1 (2.459)

|λ1| = |λ2| = |λ3| = 12

for α = −13

. |λ________1| = |λ2| is a monotonic decreasing function

of α and |λ3| is monotonic increasing. Consequently, max(|λ1|, |λ2|, |λ3|) has a minimum

at λ = −1/3. The complete range [−12

, 0] for max(|λ1|, |λ2|, |λ3|) is covered

if one takes α ∈ [−12

, 0]. This concludes the treatment for _→∞.

3. To assure that max(|λ1|, |λ2|, |λ3|) ≤ 1 for 0 < _ < ∞, the feasibility of |λ| = 1 as

solution of Equation (2.413) will be checked. In general, the solution can be complex.

Substituting λ = eiϕ in Equation (2.413) and separating the real and imaginary part

of the equation yields

cos 3ϕ − I1 cos 2ϕ + I2 cos ϕ − I3 = 0 (2.460)

sin 3ϕ − I1 sin 2ϕ + I2 sin ϕ = 0. (2.461)

Both equations must be satisfied. Expanding sin 2ϕ = 2 sinϕ cos ϕ and sin 3ϕ =

−3 cos ϕ + 4 cos3 ϕ in Equation (2.461) yields

sin ϕ(−1 + 4 cos2 ϕ − 2I1 cos ϕ + I2) = 0 (2.462)

which implies

sin ϕ = 0 ⇔ ϕ = 0, π (2.463)

4 cos2 ϕ − 2I1 cos ϕ + I2 = 1. (2.464)

(a) For ϕ = 0, Equation (2.460) yields

1 − I1 + I2 − I3 = 0 (2.465)

(b) for ϕ = π one obtains

1 + I1 + I2 + I3 = 0. (2.466)

4 cos3 ϕ in Equation (2.460), we get

cos ϕ(−3 + 4 cos2 ϕ − 2I1 cos ϕ + I2) + I1 − I3 = 0 (2.467)

and Equation (2.467) can finally be transformed into

cos ϕ = I1 − I3

. (2.468)

Substitution of Equation (2.468) into Equation (2.464) finally yields

I3(I3 − I1) + I2 − 1 = 0. (2.469)

136 LINEAR MECHANICAL APPLICATIONS

Equations (2.465), (2.466) and (2.469) cover all cases for which |λ| = 1.

(a) Equation (2.465) leads to (Equation (2.428) and Equation (2.424))

_2 = 0 (2.470)

which yields a double root for _ = 0. This was already covered previously.

(b) Substituting γ = 12

− α and β = (1 − α)2/4 in Equations (2.418) to (2.420)

and using these in Equation (2.466) yields

α2(1 + 2α)_2 − 8α(1 + 2α)ξ_ + 4 = 0 (2.471)

with roots

_1,2 = 2ξ

1 ±

1 − 1

4ξ 2(1 + 2α)

. (2.472)

For 4ξ 2(1 + 2α) < 1, there are no real solutions, for 4ξ 2(1 + 2α) ≥ 1 the solutions

are both negative since α < 0 (for α = 0, Equation (2.471) has no solution

either). Accordingly, there are no positive real solutions of Equation (2.471).

_[8ξ + 8ξ 2_ + 2ξ(1 + α2)_2 − α(1 + α)2_3] = 0. (2.473)

_ = 0 was already covered. For ξ _= 0 and α _= 0 all the coefficients in

Equation (2.473) are strictly positive (ξ > 0, −12

≤ α < 0) and _ = 0 is the

only solution. If one of them is zero but not both, the same reasoning applies

to the nonzero terms. If ξ = α = 0 the equation is satisfied for all _. This

corresponds to the classical Newmark algorithm.

Summarizing, for the parameter combinations γ = 12

− α, β = (1 − α)2/4, α ∈ [−13

, 0],

the implicit scheme is unconditionally stable and second-order accurate. The spectral radius

for different values of α and ξ = 0 is shown in Figure 2.28, whereas the effect of ξ for

α = −13

is plotted in Figure 2.29. The figure shows that there is nearly no numerical dissipation

for small _-values. For increasing values of _ the dissipation gradually increases.

In the previous derivation, it was shown that the solutions of the characteristic equation can

be complex. This leads to oscillatory damping and a corresponding (small) period error of

the solution. For more information the reader is referred to (Miranda et al. 1989), (Hilber

and Hughes 1978) and (Bathe 1995).

2.11.5 Explicit formulation

The method in Section 2.11.1 is essentially implicit owing to the formulation of _M

∗_ in

Equation (2.350) and Equation (2.357). If _M

∗_ is diagonal, the method is explicit. The

mass matrix can be made diagonal by lumping. The problem is the damping matrix _C_ and

the stiffness matrix _K_, which are usually not diagonal. In the explicit predictor–corrector

LINEAR MECHANICAL APPLICATIONS 137

0 0.5 1

1.5 2 2.5 3 3.5 4

_(−)

max(|λ1|, |λ2|, |λ3|)(−)

α = 0.0

α = −0.1

α = −0.2

α = −0.3

0.8

0.9

1.05

0.95

0.85

Figure 2.28 Spectral radius for ξ = 0

0 0.5 1

1.5 2 2.5 3 3.5 4

_(−)

max(|λ1|, |λ2|, |λ3|)(−)

ξ = 0.0

ξ = 0.1

ξ = 0.2

ξ = 0.3

0.7

0.8

0.9

1.05

0.95

0.85

0.75

0.65

Figure 2.29 Spectral radius for α = −0.3

138 LINEAR MECHANICAL APPLICATIONS

procedure introduced here (Miranda et al. 1989), the diagonalization of _M

∗_ is achieved

by simply dropping the last two terms in Equation (2.350), that is, by setting

∗_ = _M_ (2.474)

and lumping _M_. This really amounts to replacing Equation (2.339) by

_M_ _A_n+1

+ (1 + α) _C_ _ ˜V _n+1

− α _C_ _V _n

+ (1 + α) _K_ _ ˜U _n+1

− α _K_ _U_n

= (1 + α) _F_ext

n+1

− α _F_ext

, −1 ≤ α ≤ 0. (2.475)

It can be shown that this iterative scheme is second-order accurate if Equation (2.342)

and Equation (2.343) are satisfied. Furthermore, high-frequency dissipation is achieved

for α < 0. However, the explicit scheme is not unconditionally stable. Indeed, the onedimensional

equivalent of the explicit scheme corresponds to the equivalent model of

the implicit scheme in which the parameter d defined in Equation (2.410) is replaced

by 1. Equation (2.433) still applies and the explicit scheme is consistent and second-order

accurate for ξ = 0 if γ = 12

− α. However, the explicit scheme is not stable for _→∞.

To check stability, Equation (2.465), (2.466) and (2.469), which still apply, are analyzed.

β = (1 − α)2/4 is assumed throughout. Equation (2.465) reduces to _2 = 0 and deserves

no further attention. Equation (2.466) now yields

(1 − α − 2α2 − α3)_2 + 4(1 + α + 2α2)ξ_ − 4 = 0 (2.476)

leading to

_1,2 = 2(1 + α + 2α2)

1 − α − 2α2 − α3





ξ 2 + 1 − α − 2α2 − α3

(1 + α + 2α2)2

− ξ



. (2.477)

Since 1 − α − 2α2 − α3 > 0 and 1+ α + 2α2 > 0 for −12

≤ α ≤ 0, the positive root in

Equation (2.477) marks a relevant crossing of the |λ| = 1 line. Figure 2.30 shows _1 as a

function of α for different ξ values.

Equation (2.469) reduces to

_{−8ξ − 8α(1 + 2α)ξ2_ + 2ξ [−α(1 + 2α) + α3]_2 + α(1 + α)2_3} = 0. (2.478)

For ξ = 0 the only solution is _ = 0. A numerical analysis shows that also for ξ > 0

there are no positive real roots of Equation (2.478) smaller than _1 in Equation (2.477).

Summarizing, the stable regime is limited by a critical “frequency” value given by _1 in

Equation (2.477) and plotted in Figure 2.30.

2.11.6 The consistent mass matrix

The consistent mass matrix is obtained by evaluating Equation (2.24). This is usually performed

by one of the integration schemes from Section 2.3. In the previous sections, it was

explained that a force jump leads to a jump in the acceleration. To this end Equation (2.395)

has to be evaluated. This is a system of equations with the mass matrix on the left-hand side.

LINEAR MECHANICAL APPLICATIONS 139

1.5

_cr(−)

ξ = 0.0

ξ = 0.1

ξ = 0.2

ξ = 0.3

α(−)

1.9

1.8

1.7

1.6

1.4

−0.3 −0.25 −0.2 −0.15 −0.1 −0.05

Figure 2.30 Critical frequency for the explicit scheme

Equation (2.395) cannot be solved if the mass matrix is singular. Physically, this situation

cannot arise since the mass matrix is positive-definite. Indeed, the kinetic energy satisfies

K = 12

_V _T _M_ _V _ = 0 ⇔ _V _ = 0. (2.479)

Consequently,

_M_ _V _ = 0 ⇒ _V _ = 0 (2.480)

and _M_ is regular. However, _M_ can become singular because of the numerical integration.

To realize this, consider just one finite element with constant initial density and recall

that the numerical integration of Equation (2.24) amounts to

Mij ≈ ρ0

k=1

wkϕikϕjkJ

∗

k (2.481)

where Mij denotes the entry in row i and column j of the matrix _M_, wk are the weighting

functions, J

∗

k is the Jacobian determinant of the global–local coordinate transformation

and the indices K and M in Equation (2.24) were dropped since the mass matrix does

not depend on them (the mass matrix really consists of three identical submatrices, one

for each coordinate direction). The size of the matrix in Equation (2.481) is equal to the

number of nodes in the element. A matrix _a_ is singular if its determinant vanishes.

Suppose there is only one integration point. This is the case for the 8-node brick element

with reduced integration and the four-node tetrahedral element with standard integration.

Accordingly, Equation (2.481) reduces to

Mij ≈ ρ0w1ϕi1ϕj1J

∗

1 . (2.482)

140 LINEAR MECHANICAL APPLICATIONS

Furthermore, the shape functions and the location of the integration point for these elements

are such that

ϕi1 = ϕj1 =: ϕ1 ∀i, j (2.483)

and we get

Mij ≈ ρ0w1ϕ1ϕ1J

∗

1 . (2.484)

All entries in the matrix are identical: the matrix is singular. Also, 20-node brick elements

with reduced integration frequently lead to badly conditioned mass matrices. Therefore, it

is advisable to use the higher-order schemes to integrate the mass matrix.

2.11.7 Lumped mass matrix

In the explicit formulation the mass matrix is reduced to a diagonal form. This can be

performed in several ways (Zienkiewicz and Taylor 1989). Here, only one method will be

discussed, which is used in the CalculiXcode (CalculiX 2003). In this method, the lumped

mass matrix is obtained by scaling the diagonal terms of the consistent mass matrix such

that the total mass is recovered. Denoting the consistent element mass matrix by _MCij _

and the lumped element mass matrix by _MLij _ one finds

MLii = MCii

j=1MCjj

(2.485)

where Me is the total mass of the matrix, that is,

Me =

i=1

j=1

MCij (2.486)

and n is the number of nodes in the element. This rule is applied to linear elements. For

quadratic elements, a distinction is made between vertex node contributions and midside

node contributions. Denoting the set of vertex nodes by VN and the set of midside nodes

by MN, we define

α := 1i∈VN 2Ve

ϕ2

i dV

1j∈MN 2Ve

ϕ2

i dV

. (2.487)

α is a measure for the mass concentrated in the vertex nodes relative to the mass in the

midside nodes. The integration in Equation (2.487) is performed in local coordinates. The

lumped mass entries are now obtained by

MLii = MCii

_ Me

1j∈VNMCjj

__ α

1 + α

, i∈ VN (2.488)

MLii = MCii

_ Me

1j∈MNMCjj

__ 1

1 + α

, i∈ MN. (2.489)

Summing the masses in Equation (2.488) and Equation (2.489) readily shows that the total

element mass is correctly reproduced. The factor α for some widely used quadratic elements

is listed in Table 2.4.

LINEAR MECHANICAL APPLICATIONS 141

Table 2.4 The lumping factor α for

several element types.

Element type α

20-node brick element 0.2917

10-node tetrahedral element 0.1203

15-node wedge element 0.2141

p A B

p = 10 MPa

E = 210, 000 MPa

ν = 0.3

ρ = 7800 kg/m3

10 mm 10 mm

10 mm 30 mm

Figure 2.31 Geometry of the spherical shell and material data

2.11.8 Spherical shell subject to a suddenly applied uniform pressure

Consider the thick spherical shell in Figure 2.31 (only one-eighth is shown). At t = 0, a

pressure p is applied and we are interested in the radial stresses as a function of time. It is

known that the ensuing pressure waves travel at a speed c1 satisfying (Graff 1975)

c1 =

0λ + 2μ

= 6.0202 × 106 mm/s (2.490)

where λ and μ are Lam´e’s constants and ρ is the density of the material (λ and μ can

be calculated from Young’s modulus E and the Poisson coefficient ν in Figure (2.31) by

use of Equations (1.450) and (1.451)). This means that they reach the outer surface of the

shell after 8.3 × 10−6 s. One-eighth of the shell is meshed by 10 rows of 20-node brick

elements with reduced integration across the thickness and 75 elements in circumferential

direction, resulting in 750 elements.

Figure 2.32 shows the radial stress in points A and B by using the implicit α-method

with α = −0.05 and compares these results with the analytical solution (dashed lines,

142 LINEAR MECHANICAL APPLICATIONS

−5

−4

−3

−2

−1

5 10 15 20

t (10−6 s)

σrr (MPa)

r = 20 mm

r = 30 mm

Figure 2.32 Radial stress after pressure surge

(Eringen 1980)). The analytical solution applies to a hole in infinite space, and therefore

no reflection takes place. Since the radial wave equation does not exhibit dispersion (Graff

1975) the signal form is kept during propagation, although its amplitude changes. The

numerical solution hits the extremal values of the analytical solution and the time at which

they occur well. The finite element results are nodal values and therefore they are smeared

out due to the extrapolation within the element (cf Section 2.4). For times exceeding the

transversal time of the wall thickness, the wave is reflected leading to maxima at t =

1.4 × 10−5 s in B and t = 1.5 × 10−5 s in A. The subsequent maximum in A (after two

reflections) takes place at t = 1.82 × 10−5 s. The results also show that quadratic elements

tend to lead to oscillatory solutions for short-time calculations.

Пресс-релиз популярных книг

2.11 Dynamics: The α-Method

Популярные книги

Популярные статьи

Навигация