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)

1/--страниц