close

Вход

Забыли?

вход по аккаунту

?

529

код для вставкиСкачать
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 41, 1277–1296 (1998)
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
AN M(M−1 K)m PROPORTIONAL DAMPING IN EXPLICIT
INTEGRATION OF DYNAMIC STRUCTURAL SYSTEMS
A. MUNJIZA1;∗ , D. R. J. OWEN 2 AND A. J. L. CROOK 2
1Department
of Engineering, QMW, University of London, London, U.K.
of Civil Engineering, University of Wales, Swansea, U.K.
2Department
ABSTRACT
Algorithmic aspects of energy dissipation mechanisms of dynamic structural systems in conjunction with
central dierence time integration method are investigated and damping proportional to M(M−1 K)m (where
K is the stiness matrix, M is the mass matrix and m is a damping parameter) is proposed. Detailed algorithms
for M(M−1 K)m proportional damping for the central dierence time integration method are presented together
with stability criteria and numerical test problems. ? 1998 John Wiley & Sons, Ltd.
KEY WORDS: structural dynamics; damping; discrete elements
1. INTRODUCTION
With increasing computing power becoming available, ever larger and more complex problems are
being routinely considered in computational mechanics. In particular, practical shell (2 · 5 − D) and
3-D structures are being increasingly considered and such large scale nite element simulations
give rise to large, sparse matrices. The cost of solution and storage requirements of direct methods
increase dramatically for 3-D problems, as compared to 2-D, which makes direct solvers (in
static analysis) and implicit time integration methods (in dynamic analysis) less competitive in
comparison to iterative solvers and explicit and mixed time integration procedures.1
While in 2-D problems the ratio between the highest and the lowest eigenvalue of the stiness
matrix for a homogeneous mesh is proportional to N , in 3-D problems it is only proportional to
N 2=3 , where N is the number of elements. In this context, iterative solvers, explicit and mixed
explicit time integration methods and explicit codes in general are more competitive than implicit versions. This fact is recognized in many elds such as computational uid dynamics and
semiconductor device simulation, and more recently in structural problems.
In structural problems direct methods are still commonly used, but robust explicit formulations
for both linear and nonlinear problems are emerging (at present explicit formulations are usually
included as an option in large nite element systems1 ).
∗
Correspondence to: A. Munjiza, Department of Engineering, Queen Mary and Westeld College, Mile End Road, London
E1 4NS, U.K. E-mail: a.munjiza@qmw.ac.uk.
CCC 0029–5981/98/071277–20$17.50
? 1998 John Wiley & Sons, Ltd.
Received 7 March 1995
Revised 22 May 1997
1278
A. MUNJIZA, D. R. J. OWEN AND A. J. L. CROOK
While advanced implicit time integration methods have parametric control of numerical damping
of high frequencies, it is not the case with most explicit methods; especially, not with the central
dierence method, which does not have any numerical damping.21
In this work, energy-dissipation mechanisms in structural systems and their relation to the central
dierence method are investigated. In this context, a so-called M(M−1 K)m proportional damping
is introduced into the central dierence time integration method.
Theoretical background for M(M−1 K)m damping together with detailed implementation algorithm and linear test problems is presented in this paper, while application to non-linear problems
is limited to examples only (without any investigation into the accuracy or stability)—presented
in the appendix.
2. DISCRETIZATION IN SPACE
The dynamic response of a structural system is in mathematical terms an initial boundary value
problem, which comprises a set of partial dierential equations, boundary conditions and initial
values.
One approach towards formulating algorithms for solving partial dierential equations associated
with the problem is to rst discretize the spatial domain, which results in a system of ordinary
dierential equations with time as the independent variable
Kx + C ẋ + M x = p
(1)
where K is the stiness matrix, C is the damping matrix, M is the mass matrix and p is the vector
of external loads. For discretization in space, nite dierences, nite elements, nite volumes,
boundary solutions and Galerkin methods,13 in general, have all been applied in dierent situations.
However, nite elements are by far the most commonly used. Even in the context of discrete
elements,4 – 8 nite elements have been introduced to capture the deformability of each individual
body and as a result combined nite=discrete element approaches have been recently developed
for a range of industrial problems.8
3. DISCRETIZATION IN TIME
A wide range of techniques are available to solve the problem dened by equation (1). Modal
analysis is traditionally used for systems dominated by low frequencies, whilst direct integration
techniques are better suited for the transient response of systems with a wide spectrum of frequencies.
In general, direct integration methods are formulated using nite dierence or nite-elementbased discretizations in time9; 10 and are classied as linear (linear single-step or linear multi-step
methods2 ) and non-linear methods (such as Runge–Kutta methods11 and rational Runge–Kutta
methods). Alternatively, direct integration methods can be classied into implicit, explicit and
mixed.
Among the well-known implicit methods are Park’s method,12 collocation methods13 and
-methods.14 These techniques require solution of a system of algebraic equations at each time
step, and are all unconditionally stable, second-order-accurate and dissipative (collocation and
-methods permit parametric control of dissipation). The frequently employed Newmark method
can also be classied as a special case of the -method.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
AN M(M−1 K)m PROPORTIONAL DAMPING
1279
Mixed methods are based on either domain or operator-splitting techniques (although it has been
proved that both approaches are equivalent15 ). Domain-splitting schemes use dierent methods for
dierent regions of the domain,15; 16 while operator-splitting employs decomposition of the stiness
and damping matrices.17 A special class of mixed explicit methods are those methods based on
subcycling.18
Explicit methods do not require the solution of a system of algebraic equations; however, in
general they are only conditionally stable. The most popular explicit method is the central dierence
method, which is conditionally stable for linear problems, while slight energy imbalance has been
experienced with some non-linear applications.19
There have been considerable eorts to improve the performance of explicit schemes,20 – 24 and
the developed versions have in general accuracy of order greater or equal to two. However, this
has resulted in either a reduced critical time step or a need for residual forces to be calculated
more than once during each time step.
In the context of discrete element methods, explicit methods are commonly used in solution. For
combined nite=discrete element models, the authors have also employed an explicit formulation
based on the central dierence scheme, but have found that a well-behaved damping model is
essential for some problems in this class. For instance, in order to obtain the fracture pattern of
brittle rock, one has to carefully control the damping of higher frequencies, while to predict the
muck-pile in blasting simulations8 the damping of low frequencies must be controlled.
4. ENERGY DISSIPATION
Practical dynamic systems dissipate energy through damping mechanisms such as elastic hysteresis,
internal friction, plastic deformation and interaction with surrounding solids or uids. In nite
element problems damping, for instance, controls resonance, while in discrete element problems the
state of rest and the time to steady state are greatly inuenced by the damping model employed.8
Spatial discretization of a damped structural system can be written as
Kx + C ẋ + M x = p
(2)
Several cases of practical importance are listed below:
(1) Undamped system
C=O
(3)
(2) Viscous damping at small velocities
C=D
(4)
where D is a diagonal matrix
(3) Viscous damping at large velocities25
C = D ẋ
(5)
(4) Rayleigh damping
C=
? 1998 John Wiley & Sons, Ltd.
2
2
M+
K
!1
!n
(6)
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
1280
A. MUNJIZA, D. R. J. OWEN AND A. J. L. CROOK
Figure 1. Rayleigh damping
where and are damping parameters. Rayleigh damping is very often employed in engineering
applications for its exibility in the sense that it can be tuned to predominantly damp high or low
frequencies. In Figure 1 the damping ratio as a function of scaled frequency (actual frequency
divided by the maximum frequency of the system) is shown for dierent damping parameters and .
5. HIGH FREQUENCIES IN THE CENTRAL DIFFERENCE METHOD
The central dierence time integration method does not introduce any numerical damping. Additionally, the stress and strain elds due to high frequencies are overestimated. This is illustrated
in Figures 2 and 3, where free oscillations of the initial value problem
x + x = 0;
x(0) = 1;
ẋ(0) = 0
(7)
are integrated by the central dierence method by employing a time step 0·7 s (Figure 2) and 1·9 s
(Figure 3).
Figure 4 shows the minimum value of x(t) in the interval 0¡t¡∞ as a function of the time
step employed. Integration errors of this type can cause problems in some applications, such as
the simulation of progressive fracturing of brittle rock.
If K proportional Rayleigh damping is introduced, the high frequencies are eectively damped.
However, while M proportional Rayleigh damping can be implicitly formulated in the central
dierence method, it is not the case with K proportional damping,21 which means that the critical
time step of the central dierence method is decreasing with increasing damping. In Figure 5, a
so-called step ratio (the critical time step for a damped system divided by the critical time step
for an undamped system) is shown. It can be seen that the critical time step is reduced by more
than a factor of two if the highest frequency of the system is critically damped.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
AN M(M−1 K)m PROPORTIONAL DAMPING
1281
Figure 2. SDOF system integrated by the central dierence scheme and with a time step 35 per cent of critical
Figure 3. SDOF system integrated by the central dierence scheme and with a time step 95 per cent of critical
6. DEFINITION OF M(M−1 K)m PROPORTIONAL DAMPING
As mentioned earlier, the central dierence explicit time integration method in its standard form is
employed in this work (in this context in further text it is assumed that K is positive-denite symmetric real matrix, while M is positive-denite diagonal real matrix). The only novelty introduced
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
1282
A. MUNJIZA, D. R. J. OWEN AND A. J. L. CROOK
Figure 4. Eect of time step on the accuracy of the central dierence scheme
Figure 5. Eect of K proportional damping on the critical time step
in this work is damping, i.e. the matrix C dened by
C = 2M(M−1 K)m
(8)
where is the damping parameter and m (restricted to m¿0) is a so-called damping exponent.
Thus, system (2) can be written as
Kx + 2M(M−1 K)m ẋ + M x = p
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
AN M(M−1 K)m PROPORTIONAL DAMPING
1283
M−1 Kx + 2(M−1 K)m ẋ + x = M−1 p
(9)
C = 2M
(10)
or
For m = 0
and a mass proportional damping case is recovered. In a similar way, if m = 1 a stiness proportional Rayleigh damping is achieved, while all frequencies are damped critically if m = 0·5.
7. REDUCTION TO A SDOF PROBLEM
The modal decomposition of the system (9) can be performed by solving the undamped free
vibration equation
M−1 Kx + x = 0
or
Kx + M x = 0
(11)
which leads to the eigenvalue problem
K = M
(12)
solution of which are n eigenvalues i , and n eigenvectors
T
i M
T
i M
T
i K
T
i K
j
= 0;
i 6= j
j
= 1;
i=j
= 0;
i 6= j
j
j = i ;
i
where
i=j
(13)
A modal decomposition procedure applied to system (9) yields the corresponding single degree
of freedom (SDOF) equations which appear as
!i2 qi + 2!i2m q̇i + qi = i
where !i is the frequency of the ith mode shape
i,
(14)
while qi is the mode amplitude.
Proof. Assume the solution in the form
x=
Substitute into (9)
K
n
P
j=1
n
P
j=1
!
qj
j
? 1998 John Wiley & Sons, Ltd.
+ 2[ M(M
−1
m
K) ]
qj
j
n
P
j=1
!
q̇j
j
+M
n
P
j=1
!
qj
j
=p
(15)
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
1284
A. MUNJIZA, D. R. J. OWEN AND A. J. L. CROOK
Multiplication by
T
i K
n
P
j=1
T
i
qj
leads to
!
T
−1
m
i [M(M K) ]
+ 2
j
Now from (12) and (13)
T
i K
T
i M
n
P
j=1
n
P
j=1
(M
2
T
−1
m
i [M(M K) ]
n
P
j=1
j=1
!
q̇j
+
j
T
i M
n
P
j=1
!
qj
=
j
T
i p
!
qj
= i qi = !i2 qi
j
!
qj
m
K)
q̇j
= qi
j
(M−1 K)
−1
n
P
j
j
!
j
= M−1 (K j ) = M−1 !j2 M
=
!j2m j
= 2
T
i M
n
P
j=1
j
= !j2
j
!
q̇j !j2m j
= 2!i2m q̇i
Substitution into (15) results in equation (14).
If index i is suppressed, we can write
!2 q + 2!2m q̇ + q = (16)
Thus, the behaviour of M(M−1 K)m damped systems can be investigated on a SDOF problem
(16), for which purpose equation (16) can be written as the homogeneous form
q + 2!2m q̇ + !2 q = 0
(17)
Critical damping for frequency ! is 2!, thus the damping ratio z, i.e. actual damping divided by
critical damping is
z=
2!2m
= !2m−1
2!
(18)
In Figure 6, damping ratio z as a function of scaled frequency !=!n (where !n is the highest
frequency in the system) is shown. In each case the parameter is chosen so that the frequency
!n is critically damped. If m = 0·5, damping ratio is constant over whole range of frequencies. If
m¿0·5, damping ratio for higher frequencies is larger than damping ratio for smaller frequencies.
If m¡0·5, damping ratio for higher frequencies is smaller than damping ratio for smaller frequencies. Thus damping over whole range of frequencies can be controlled by changing m.
8. CALCULATION OF M(M−1 K)m DAMPING
It is noted that no new time integration method is introduced in this work. The central dierence
time integration method is employed instead. However, at each time step damping vector d is
calculated from the known velocity of the system ẋ
d = 2M(M−1 K)m ẋ
? 1998 John Wiley & Sons, Ltd.
(19)
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
AN M(M−1 K)m PROPORTIONAL DAMPING
1285
Figure 6. Damping ratio for M(M−1 K)m proportional damping
To compute damping vector exactly would be prohibitive in terms of CPU time, thus an approximate solution is introduced instead. Each component di of the damping vector d is approximated by
di = mi 2(1m ai + nm bi )
(20)
where mi is lumped mass, 1 and n are, respectively, the lowest and the highest eigenvalue of
the system (11), while ai and bi are unknown parameters that can be calculated if collocation at
m = 0 and m = 1 is assumed
2mi dˆi = mi 2(10 ai + n0 bi )
2mi d˜i = mi 2(11 ai + n1 bi )
(21)
where dˆi is a component of vector d̂ = (M−1 K)0 ẋ = ẋ and d˜i is a component of vector d̃ =
(M−1 K)1 ẋ = M−1 Kẋ From equations (21), parameters ai and bi are calculated for each degree of
freedom
ai =
dˆi n − d˜i
n − 1
(22)
bi =
d˜i − dˆi 1
n − 1
(23)
The above approximation is equivalent to the approximation of the function m , by a straight
line as illustrated for m = 0·5 in Figure 7. For frequencies !1 and !2 the above introduced approximation yields exact result, while for all other frequencies the damping is underestimated, thus
for = 1 and m = 0·5 only the highest and the lowest frequencies are critically damped while all
other frequencies are underdamped.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
1286
A. MUNJIZA, D. R. J. OWEN AND A. J. L. CROOK
Figure 7. Approximation of M(M−1 K)m damping
9. CENTRAL DIFFERENCE METHOD AND M(M−1 K)m DAMPING
As explained earlier, the central dierence time integration method is employed in this work. At
each time step damping vector is calculated component by component as given by equation (20).
During implementation it was found convenient to write equation (20) in such a way that damping parameters supplied as input bear physical meaning of damping ratio for lowest frequency
z1 = !12m−1 and damping ratio for the highest frequency zn = !n2m−1 . With this modications
equation (20) becomes
di = (2z1 !1 ai + 2zn !n bi )mi
(24)
If at time step i − 1 time is t and at time step i time is t + t, the velocity and displacements
at time step i can be expressed in terms of velocity and displacements at time step i − 1 in the
standard central dierence form with implicit damping term21
ẋi =
bi−1
ri−1
ai−1
+
+
t
1 + 2z1 !1 t
1 + 2zn !n t
mi
xi = xi−1 + ẋi t
or
1
ẋi =
(n − 1 )
d˜i−1 − dˆi−1 1
dˆi−1 n − d˜i−1
+
1 + 2z1 !1 t
1 + 2zn !n t
xi = xi−1 + ẋi t
!
+
ri−1
t
mi
(25)
where ri is a component of vector of residual forces r = p − Kx.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
AN M(M−1 K)m PROPORTIONAL DAMPING
1287
In the limit
xi − xi−1
xi−1 + ẋi t − xi−1
= lim
t→0
t→0
t
t
= ẋi
bi−1
ri−1
ai−1
+
+
t − (ai−1 + bi−1 )
ẋi − ẋi−1
1 + 2z1 !1 t
1 + 2zn !n t
mi
= lim
lim
t→0
t→0
t
t
−2zn !n bi−1
ri−1
−2z1 !1 ai−1
+
+
= lim
t→0 1 + 2z1 !1 t
1 + 2zn !n t
mi
ri−1
= −2z1 !1 ai−1 − 2zn !n bi−1 +
mi
di−1
ri−1
=−
+
mi
mi
lim
the original dierential equation is recovered.
10. NUMERICAL STABILITY
Numerical stability of the central dierence method in conjunction with M(M−1 K)m proportional
damping can be investigated on SDOF system for frequency !i with index i suppressed (16),
!2 q + 2!2m q̇ + q = 0
(26)
The algorithmic procedure (25) when applied to this system results in
q̇i−1 n − q̇i−1 !2
n − 1
q̇i−1 !2 − q̇i−1 1
=
n − 1
ai−1 =
bi−1
and thus yields equivalent integration procedure for each modal form
q̇i−1
!2 − 1
n − ! 2
+
− !2 qi−1 t
q̇i =
(n − 1 ) 1 + 2z1 !1 t
1 + 2zn !n t
qi = qi−1 + q̇i t
(27)
or
q̇i = (1 − )q̇i−1 − !2 tqi−1
qi = qi−1 + q̇i t = qi−1 + [(1 − )q̇i−1 − !2 tqi−1 ]t = (1 − )q̇i−1 t + (1 − !2 t 2 )qi−1
(28)
where
=1 −
? 1998 John Wiley & Sons, Ltd.
!2 − 1
n − !2
+
1 + 2z1 !1 t
1 + 2zn !n t
1
n − 1
(29)
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
1288
A. MUNJIZA, D. R. J. OWEN AND A. J. L. CROOK
Figure 8. Eect of M(M−1 K)m proportional damping on critical time step of the central dierence scheme
Since 1 6!2 6n , parameter is bounded by 0661 and (28) can be written as
1 − −!2 t 2
q̇i−1 t
q̇i t
=
1 − 1 − !2 t 2
qi
qi−1
(30)
Parameter is in some sense scaling the velocity at each time step, so in further discussion it
is referred to as the ‘velocity scaling parameter’, while in a similar way the parameter s = 2t=!
is referred to as the ‘time step scaling parameter’. By checking the spectral radius of the matrix
1 − −!2 t 2
A=
1 − 1 − !2 t 2
the time step scaling parameter s can be expressed as a function of the velocity scaling parameter (Figure 8). It can be seen from the gure that the critical time step for the central
dierence method in conjunction with the M(M−1 K)m damping is in fact never reduced by more
then 30 per cent compared to the critical time step for the central dierence method applied to
the undamped system. In other words, a wide range of damping parameters is available with very
little eect on the stability of the central dierence time integration scheme.
11. LINEAR TEST PROBLEMS
In this section, we consider direct integration of the system (9) by using the algorithm for
M(M−1 K)m damping in an explicit code. Since the central dierence time integration employed
(25) is from the class of linear time integration methods, the same result would be obtained if
one rstly performed modal decomposition and then integrated each mode or if direct integration
was performed without modal decomposition.21
In this context, we assemble a numerical example of the system (9) such that the highest
frequency of the system is !n = 1 and the lowest frequency of the system is !1 = 0·01. The
damping of modes, ! = 1, ! = 0·1 and ! = 0·01, with the assumption that the initial displacement
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
AN M(M−1 K)m PROPORTIONAL DAMPING
1289
Figure 9. Displacements for the frequency ! = 0·01 and damping parameters z1 = 0 and zn → ∞
and initial velocity for all modes is q(0) = 1·0 and q̇(0) = 0·0, is investigated. Damping parameters
z1 and zn are variable.
The damping of mode ! = 0·01 (the lowest frequency of the system) is rstly investigated.
In Figure 9 displacements for z1 = 0 and zn → ∞ are shown. Actually, this is the case of free
oscillations. In other words, damping parameter zn does not have any eect on the damping of the
lowest frequency. The opposite is also true, i.e. damping parameter z1 does not have any eect on
the damping of the highest frequency !n . If z1 is increased to z1 = 0·25, underdamped oscillations
of the mode ! = 0·01 are obtained, while parameter z1 = 1 results in critical damping (Figure 10),
and z1 ¿1 results in overdamping.
Secondly, damping of mode ! = 1 (the highest frequency of the system) is studied. As mentioned
above, damping of this mode is not aected by the value of the parameter z1 . However, parameter
zn has a great inuence on the damping of this mode. Figure 11 shows the underdamped case
with zn = 0·5 and also the overdamped cases with zn = 2·5 and zn → ∞.
Finally, damping of the mode ! = 0·1 is presented. This mode has a frequency which is larger
than the lowest and smaller than the highest frequency of the system, thus both parameters z1 and zn
have an eect on its damping. Displacements for z1 = 0 combined with zn = 0·5, zn = 2·5 and
zn → ∞ are presented in Figure 12. Displacements for z1 = 1 combined with zn = 0·5, zn = 2·5 and
zn → ∞ are presented in Figure 13. In all cases the oscillations are underdamped.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
1290
A. MUNJIZA, D. R. J. OWEN AND A. J. L. CROOK
Figure 10. Displacements for the frequency ! = 0·01 and any value of zn
Figure 11. Displacements for the frequency ! = 1·0 and any value of z1
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
AN M(M−1 K)m PROPORTIONAL DAMPING
1291
Figure 12. Displacements for the frequency ! = 0·1 and z1 = 0
Figure 13. Displacements for the frequency ! = 0·1 and z1 = 1
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
1292
A. MUNJIZA, D. R. J. OWEN AND A. J. L. CROOK
All the results presented are obtained by employing a time step t = 0·5 s. The results show
that by varying damping parameters z1 and zn a reasonable control over damping for the entire
range of frequencies can be achieved without any signicant cost in terms of CPU or algorithmic
complexity, which in short is a major contribution presented in this work.
12. CONCLUDING REMARKS
To summarize the results, the approach to the modelling of structural damping in an explicit code
advocated in this paper provides a remarkably simple methodology for
1. damping high frequencies, or
2. damping low frequencies, or
3. damping both high and low frequencies.
The parameters that control the proposed M(M−1 K)m damping have a physical interpretation
as the damping ratio of the highest and damping ratio of the lowest frequency and thus can be
tuned to achieve the required damping over a range of frequencies.
The algorithm proposed does not signicantly increase the processing time at each time step.
Our estimation is that for linear structural dynamics problems, the algorithm with M(M−1 K)m
damping may consume up to 5 per cent more time in comparison to the same algorithm without
any damping.
The reduction of the critical time step due to the presence of M(M−1 K)m damping does not
exceed 30 per cent for the entire range of damping parameters that can be supplied (0 6 z1 6∞
and 0 6 zn 6∞), which is a small decrease in comparison to K proportional damping.
APPENDIX
Non-linear applications
M(M−1 K)m proportional damping has been developed for linear problems. However, it has
been also incorporated into the combined nite-discrete element method, which had central dierence time integration method already implemented. The only change needed is calculation of the
damping vector at each time step as explained earlier.
The stability and accuracy of the central dierence time integration method in such highly nonlinear problems is an open question. Despite this the central dierence time integration method remains widely used in both linear and non-linear problems. Incorporation of M(M−1 K)m
proportional damping into the method changes only the way in which components of the damping
vector are calculated at each time step. In this sense introduction of M(M−1 K)m proportional
damping is equivalent to changing from M to K proportional damping. Thus, it is felt important
to add numerical examples concerning non-linear analysis. In this case typical combined nite
discrete element problems are shown. It is beyond scope of this paper to go into details of the
combined nite discrete element method (see Reference 8).
Example 1. In Figure 14 an initial mesh for the impact of a square bar is shown. The bar
is 500 mm long and 100 mm deep. The support is provided through contact with steel bolts. The
bottom supports are at distance 97·6 mm from the bar ends, while the upper supports are at distance
46·8 mm from the ends. The bar is made of concrete with density 2340 kg=m3 , compression strength
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
AN M(M−1 K)m PROPORTIONAL DAMPING
1293
Figure 14. Bar impact problem
31 MPa, tensile strength 3·15 MPa, elastic modulus 25·8 GPa, Poisson ratio 0·17 and fracture energy
0·3 kJ=m2 .
The apparent contact surface at the moment of impact is 100 × 100 mm, with an impact velocity
of 5·85 m=s and impact mass of 30 kg. The problem is solved as a plane stress problem, with impact
conditions achieved through supplying a uniform density of the impactor, which is assumed to be
elastic with elastic modulus of 210 GPa and Poisson ratio of 0·33.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
1294
A. MUNJIZA, D. R. J. OWEN AND A. J. L. CROOK
Figure 15. Cleavage problem
The motion and fracturing sequence is shown in Figure 14. The rst cracks appear approximately
on lines connecting the bottom supports and the bottom edges of the impactor. This cracks widen
and new cracks open until the bar is transformed into a kinematic mechanism with individual
fragments moving apart.
It should be emphasized that proper momentum balance (in terms of damping being a result of
internal forces as in the case of K proportional damping or external forces as in the case of M
proportional damping) is of paramount importance for at least two reasons: rst, an imbalance in
momentum can result in an incorrect crack pattern and second, postfracturing motion is a result
of kinetic energy and momentum transferred to individual fragments through prefracturing plastic
mechanisms and through the fracturing process itself.
In this problem the momentum balance is preserved through M(M−1 K)m proportional damping by specifying damping parameters z1 = 0 and zn = 1. The parameter z1 = 0 ensures the
preservation of momentum, while the parameter zn = 1 silences high frequencies and thus contributes to stabilization of the fracturing process in a sense that the main cracks propagate, while
the cracks in close proximity to the main cracks close through energy dissipation and stress
release.
Example 2. In Figure 15 an initial nite element mesh for the cleavage of a concrete specimen
100 × 100 mm in size is shown. A plane strain situation is assumed with elasticity modulus E =
20 GPa, Poisson ratio = 0·17, tensile strength ft = 2·0 MPa and fracture energy Gf = 0·2 kJ=m2 .
The impactor object is assumed to move with a constant velocity of 1 m=s.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
AN M(M−1 K)m PROPORTIONAL DAMPING
1295
It can be seen that although, the main crack starts at the top of the specimen and propagates
towards the bottom, some parallel cracks in close proximity to the main crack also appear. This is
not only due to the kinematics and nature of the fracture processing algorithm itself,8 but also due
to the dynamic nature of the problem, in which a M(M−1 K)m proportional damping has benecial
eect through silencing high frequencies. For this purpose damping parameters z1 = 0 and zn = 1
are assumed, i.e. the rigid motion is not damped at all, while the highest frequency is damped
critically.
REFERENCES
1. K. J. Bathe, J. Walczak and H. Zhang, ‘Some recent advances for practical nite element analysis’, Comput. Struct.,
47, 511–521 (1993).
2. T. J. R. Hughes, ‘Analysis of transient algorithms with particular reference to stability behaviour’, in T. Belytschko
and T. J. R. Hughes (eds.), Computational Methods in Transient Analysis, North-Holland, Amsterdam, 1983.
3. C. A. J. Fletcher, Computational Galerkin Methods, Springer, Berlin, 1984.
4. J. Lemos, R. D. Hart and P. A. Cundall, ‘A generalized distinct element program for modelling jointed rock mass’,
in O. Stephansson (ed.), Fundamentals of Rock Joints, Centek Publishers, Lulea, 1985.
5. J. R. Williams, G. Hocking and G. G. W. Mustoe, ‘The theoretical basis of discrete element methods’, NUMETA ’85
Numerical Methods in Engineering; Theory and Applications, Balkema, Rotterdam, 1985.
6. J. R. Williams and G. Mustoe, Proc. 2nd U.S. Conf. on Discrete Element Methods, MIT, MA, 1993.
7. G. Mustoe and J. R. Williams, Proc. 1st U.S. Conf. on Discrete Element Methods, Golden, CO, 1989.
8. A. Munjiza, D. R. J. Owen and N. Bicanic, ‘A combined nite-discrete element method in transient dynamics of
fracturing solids’, Int. J. Engng. Comput., 12, 145 –174 (1995).
9. G. M. Hulbert, ‘Time nite element methods for structural dynamics’, Int. J. Numer. Meth. Engng., 33, 307–331
(1992).
10. T. J. R. Hughes and A. Brooks, ‘A theoretical framework for Petrov-Galerkin methods with discontinuous weighting
functions: Application of the streamline upwind procedure’, in R. H. Gallagher (ed.), Finite Elements in Fluids, Vol. 4,
Wiley, London, 1982.
11. R. Cortell, ‘Application of the fourth-order Runge–Kutta method for the solution of high-order general initial value
problems’, Comput. Struct., 49, 897–900 (1993).
12. K. C. Park, ‘Evaluating time integration methods for nonlinear dynamic analysis’, in T. Belytschko, J. R. Osias and
P. V. Marcal (eds.), Finite Element Analysis of Transient Non-Linear Behaviour, Applied Mechanics Simposia Series,
ASME, New York, 1975.
13. H. M. Hilber and T. J. R. Hughes, ‘Collocation, dissipation and overshoot for time integration schemes in structural
dynamics’, Earthquake Engng. Struct. Dyn., 6, 99 –118 (1978).
14. H. M. Hilber, T. J. R. Hughes and R. L. Taylor, ‘Improved numerical dissipation for time integration algorithms in
structural dynamics’, Earthquake Engng. Struct. Dyn., 5, 283 –292 (1977).
15. F. Armero and J. C. Simo, ‘A new unconditionally stable fractional step method for non-linear coupled
thermomechanical problems’, Int. J. Numer. Meth. Engng., 35, 737–766 (1992).
16. K. C. Park and C. A. Felippa, ‘Partitioned analysis of coupled problems’, in T. Belytschko and T. J. R. Hughes (eds.),
Computational Methods in Transient Analysis, North-Holland, Amsterdam, 1983.
17. D. M. Trujillo, ‘An unconditionally stable explicit algorithm for structural dynamics’, Int. J. Numer. Meth. Engng.,
11, 1579 –1592 (1977).
18. T. Belytschko and N. D. Gilbertsen, ‘Implementation of mixed time integration techniques on a vectorised computer
with shared memory’, Int. J. Numer. Meth. Engng., 35, 1803 –1828 (1992).
19. A. Munjiza, N. Bicanic, D. R. J. Owen and Z. Ren, ‘The central dierence time integration scheme in contact impact
problems’, in Bicanic et al. (eds.), Proc. NEC-91; Int. Conf on Nonlinear Engineering Computations, Pineridge Press,
Swansea, 1991.
20. J. Kujawski, C. Miedzialowski, L. Krasuski and R. H. Gallagher, ‘Analysis of higher order explicit time integration
operators for non-linear problems’, Proc. NUMETA 87, Swansea, 1987.
21. P. I. Kimmark and W. G. Gray, ‘Fourth-order accurate one-step integration methods with large imaginary stability
limits’, Numer. Methods Part. Di. Equat., 2, 63 –70 (1986).
22. O. C. Zienkiewicz, W. L. Wood, N. W. Hine and R. L. Taylor, ‘A unied set of single step algorithms, Part 1: general
formulation and applications’, Int. J. Numer. Meth. Engng., 20, 1529 –1552 (1984).
23. C. Ho and R. L. Taylor, ‘Higher derivative explicit one step methods for non-linear dynamic problems, Part 1: design
and theory’, Int. J. Numer. Meth. Engng., 29, 275 –290 (1990).
24. R. M. Thomas and I. Gladwell, ‘Variable-order variable-step algorithms for second-order systems, Part 1: the methods’,
Int. J. Numer. Meth. Engng., 26, 55 – 80 (1988).
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
1296
A. MUNJIZA, D. R. J. OWEN AND A. J. L. CROOK
25. G. D. Stefanov, ‘Investigation into the calculation of aerodynamic damping of structural assemblies’, Comput. Struct.,
46, 397– 401 (1993).
26. D. Gottlieb and S. A. Orszag, Numerical Analysis of Spectral Methods; Theory and Applications, SIAM, Philadelphia,
1977.
27. J. C. Simo and K. K. Wong, ‘Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy
and momentum’, Int. J. Numer. Meth. Engng., 31, 19 –52 (1992).
28. D. Greenspan, ‘Conservative numerical methods for x = f(x)’, J. Comput. Phys., 56, 28 – 41 (1984).
29. G. Dahlquist, ‘A special stability problem for linear multistep methods’, BIT, 3, 27– 43 (1963).
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng., 41, 1277–1296 (1998)
Документ
Категория
Без категории
Просмотров
3
Размер файла
341 Кб
Теги
529
1/--страниц
Пожаловаться на содержимое документа