Temporal Instabilities and Dissipative Structures in Fatigued Metals
Michael V. Glazov, author
Prof. Campbell Laird, scientific advisor
Part 1. The Phenomenon of the Dislocation Core and Collective Mesoscale Dislocation Effects In Metallic Alloys
Dislocations and grain boundaries to a significant extent affect, and consequently control, the mechanical and physical properties of materials, [2-4]. The importance of the concept of the "dislocation" is difficult to overestimate *, and since its inception independently by Orowan, Polanyi and Taylor in 1934, dislocations have been studied very extensively both experimentally and theoretically. 
During the last 25 years, with the advent of increasingly more powerful computers and reliable potentials of interatomic interactions, atomistic modeling started playing an important role in these efforts, . As often happens, the properties and limitations of an instrument (digital computer) determine the fundamental length and time scales of the atomic systems that could be studied with its help. In molecular dynamics, in order to maintain the required accuracy and stability of integration algorithms, the value of the time step should usually not exceed 10-14 s, [g]. This requirement imposes limitations on the number of particles in such numerical experiments: a typical simulation run on a mini-supercomputer (e.g., "Stardent-3040") takes from 10 to 30-40 hours and usually involves from several thousands to tens of thousands atoms.
Thus it should not be surprising that attention of the computational materials science community has been attracted to studying the atomic structure of single dislocations and to "the phenomenon of the dislocation core", . This approach has led to the impressive achievements in relating mechanical properties of alloys to the structure and properties of individual dislocations. Examples include: 1). non-planar spreading of the cores of 1/2<111> screw dislocations into several planes of the <111> zone and the subsequent non-Schmid behavior and tension-compression asymmetry of the yield stress in bcc metals [9,10]; 2). temperature- and orientation-dependent "sessile-glissile" transformations of the cores of <011> screw dislocations in some L12- intermetallics and the anomalous yield behavior of such compounds, e.g. Ni3AI (or of hexagonal DO19 Ti3AI, ).
Because of their geometric, crystallographic nature and long-range elastic fields the properties of individual dislocations and collective dislocation phenomena are connected, . This statement can be illustrated using our molecular dynamics simulations of dislocation propagation and dynamics in Ni3AI. We can start with an individual dislocation at a low temperature (T=lOK), introduce the appropriate boundary conditions and obtain an equilibrated dislocation core ** ,Fig. 1 In order to study the mobility of this dislocation at different temperatures, a sufficiently large shear stress is applied to the block of some 10,000 atoms under consideration and, using a suitable visualization technique, one can further monitor the behavior of this block. At all temperatures up to 300K the behavior of the dislocation remains "individual:" it propagates (Fig. 2b) and goes out of the crystal on its left side (Fig. 2c), leaving a characteristic surface slip step. However, at temperature 600K copious dislocation multiplication can be observed instead, and the propagation pattern of the solid-state amorphisation very much resembles a percolation process, Fig.3c.
Figure 1. The results of our molecular dynamics simulation of a relaxed stress free configuration of Ni3AI edge dislocation with b = a0 on (001) glide plane. After 300 atomic vibration periods (12,000 MD-steps) of thermal equilibration at T=lOK. Solid circles - Ni atoms, open circles - Al atoms.
Figure 2. Motion of the edge ( b=aO[O10] ) dislocation in Ni3AI at T=300K under a constant shear stress of 10,480 Mpa (8% of C44). Stress was applied counterclockwise on horizontal and clockwise on vertical surfaces. (a) 4.5*10-12 s (9,000 time steps); (b) 5.0*10-12 s (10,000 time steps); (c) 6.0*10-l2 s (12,000 time steps), - dislocation goes out of the crystal.
Figure 3. Dislocation multiplication and amorphisation in Ni3AI at T=600K under a constant shear stress of 10,480 Mpa. (a) 4.95*10-12 s (9,900 steps); (b) 6.0*10-l2 s (12,000 steps); (c) 7.0*10-12 s (14,000 steps).
Our obvious inability to treat instabilities in fatigue from the atomistic point of view resulted in an attempt to find a non-traditional approach to the problem. In this paper the methods of self-organization theory and nonlinear dynamics are applied for the first time to the studying of the PLC-effect in cyclically deformed (fatigued) metals. This approach finds its justification in the fact that many strongly non-equilibrium nonlinear systems of different physico-chemical nature exhibit a tendency to self-organization (patterning) either in the form of emerging temporal rhythms, or spatial scales, or both, [14,15], and we demonstrate that fatigued metallic alloys fall into this behavior as well. Because of the similarity of pattern-formation analysis in chemical systems and in the theory of flames, fluids, ferrofluids, magnetic bubble materials, it is often referred to as "the self organization theory", [ 1, 13-15]
The main goal of this paper was to establish, via nonlinear dynamical modeling, the fundamental connection between microscopic dislocation mechanisms and macroscopic mechanical response of cyclically deformed metallic alloys. In order to achieve this goal we proceed as follows. Firstly, the experimental data on temporal patterning and the PLC effect will be considered and limitations that must be imposed on mathematical modeling will be discussed. Secondly, models of the "classical" PLC-effect will be critically analyzed, and modifications necessary to describe cyclic loading will be introduced. Finally, we will compare the existing experimental data to the results of our computer modeling and make conclusions.
Part 2. Experimental Data: "Classical" PLC-Effect and Its Counterparts In Fatigue: "Yan-Hong-Laird" and "Neumann" Bursts
Several examples of temporal instabilities (or, using Prigogine's terminology, "temporal dissipative structures", ) are given in Fig. 4 a,b,c. The "classical" PLC effect manifests itself in the form of small stress spikes on the stress vs time curves, Fig. 4a. This experimental data demonstrates serrated flow in a Cu-3.3%Sn binary alloy at 255deg. C and different values of applied strain rate, .
Figure 4. Temporal instabilities: (a). The "classical" Portevin - Le-Chatelier effect in a Cu-3.3%Sn binary alloy at 255deg.C and different values of applied strain rate, ; (b). Stress serrations (small spikes on stress vs time curve) in a cyclically deform ed binary alloy Cu-16at.%AI under strain control (i.e. strain is changed harmonically with time), ; (c). Strain ("Neumann") bursts in a ramp-loaded Al crystal, .
Several features critical to successful modeling of this phenomena should be noted upon careful consideration of the experimental data. First, we notice the nonlinear character of the temporal stress serrations (or strain bursts) in all cases. This is especially clear on Fig. 4b: the stress response of the sample is highly nonlinear while strain control imposed on the system "sample+machine" is harmonic; there is also a phase shift between the maxima on the stress vs time and strain vs time curves. Consequently, one should expect that the system of governing equations modeling such effects should be nonlinear, too (in the first approximation the effects of spatial coupling are neglected and yielding is supposed to be homogeneous in the whole volume of the sample). Further examination of experimental data shown in Fig. 4 reveals that there are two very different time scales involved: the periods of bursts of stress or strain are succeeded by time intervals with very gradual changes of these properties in time. In mathematics such multi-scale behavior results in the systems of ordinary differential equations which are usually called "stiff", and special methods should be used to integrate them reliably and within a reasonable computation time.
Finally, in all cases collective mesoscopic dislocation interactions should be considered. Without proper understanding of the underlying dislocation mechanisms of the PLC- effect and its connection to macroscopic stress-strain response of a material, it is impossible to understand the "Yan-Hong-Laird bursts" or the "Neumann bursts" in fatigued metallic alloys.
These considerations stimulated us to carefully examine the existing models for the Portevin - Le-Chatelier effect and to choose the one that would be most appropriate for modeling temporal instabilities in fatigue, with the introduction of the necessary corrections and modifications.
Part 3. Dynamical Modeling of the PLC Effect
Because of the practical importance of the PLC-effect numerous attempts have been made, with varying degrees of success, to model it numerically.
In the works of Kubin and coworkers, [20,21] the onset and disappearance of jerkyflow was connected to the behavior of two dislocation populations - mobile, Pm, and forest dislocations, Pf . The following set of two weakly nonlinear governing equations was proposed, :
where b stands for the modulus of the Burgers vector and the other parameters are related to the multiplication of mobile dislocations (C1), to their annihilation through mutual interactions (C2) to the storage of mobile dislocations through interaction with forest dislocations (C3) and to recovery (C4) d stands for strain increment, [20,21].
This model has been used for making interesting predictions on the critical conditions for the PLC-effect and for determining the strain rate intervals in which jerky flow can be observed. Several experimental systems have been considered: Cu-Mn, Al-Mg, Cu-Zn, Au-Cu, 121], and a good description of the experimental data has been obtained. Physically the model is well justified and employs quite reasonable assumptions about annihilation, trapping and interactions between dislocations. However, from the mathematical point of view it is oversimplified and does not satisfy the important requirement of stability of nonlinear stress serrations with respect to small perturbations in the initial conditions ("structural stability" of temporal dissipative structures, ). Structurally unstable systems under certain conditions may exhibit oscillatory behavior, but the periods and amplitudes of such oscillations will be determined by initial conditions, not by the properties of a nonlinear system itself. Since jerky flow is a manifestation of the fundamental material instability, it is necessary that in a sufficiently broad interval of initial conditions imposed on the dislocation densities, stress serrations of approximately the same frequency and amplitude be obtained.
It has been proved in the theory of nonlinear oscillations that time-periodic solutions satisfying this requirement do exist, and they are called "limit cycles", . According to the fundamental theorem proved by Hanusse, , it is impossible to have a limit cycle surrounding an unstable node or focus in a system of differential equations with only two variables if nonlinearity of the system is quadratic (or lower). Consequently, two choices are available: 1). a system of only two equations with a cubic nonlinearitv, in which case we come to the famous model proposed by Prigogine and called "the Brusselator" (in the theory of nonlinear oscillations it plays the same role as the concept of a harmonic oscillator - in the theory of linear oscillations, or the Heisenberg model - in the theory of ferromagnetism, ); 2). our second choice is to increase the number of equations to three or more, retaining only quadratic nonlinearities.
In order to get cubic nonlinearity in the system of differential equations one must necessarily assume that tri-molecular quasi-chemical reactions among dislocation populations are possible. From chemistry we know that such events are extremely rare, and therefore the cascades of intermediate dislocation reactions must be introduced in order to justify the overall cubic nonlinearity of the governing equations, .
We have thus been motivated to look for a nonlinear model containing more than two equations that can give, under certain conditions, the limit cycle temporal oscillations. A very interesting approach has been proposed by Ananthakrishna et al., [23,24] which satisfies all of the criteria discussed above, and it was this model (with the necessary modifications given below) that we used in order to describe the temporal instabilities in fatigued metals - the Yan-Hong-Laird bursts and the Neumann bursts.
The Ananthakrishna model introduces three dislocation populations - g , mobile (glissile) dislocations, s - immobile (sessile) and i - dislocations with atmospheres of point defects that mimic Cottrell's idea of dynamic strain aging, . The following quasi-chemical reactions among dislocations of these three types were considered:
According to , the reaction (2a) describes the process of dislocation generation by multiple cross-glide mechanism (or any other suitable process such as a Frank-Read, or a Bardeen-Herring source), ~ is the breeding constant and Vg - the average velocity of the glissile dislocations. This process gives a possibility to connect mesoscale dislocation dynamics to macroscopic stress-strain response of the material (see below). The second reaction in the set describes the conversion of two glissile dislocations into two Cottrell type ones. Annihilation of dislocations is taken care of in the equations (2c) and (2d). Finally, the equations (2e), (2f) and (2g) describe the mutual transformations of dislocations ( mobilization due to the applied stress or thermal activation, transformation of mobile dislocations into dislocations with atmospheres and their subsequent trapping). These three equations form the feedback loop
which under certain well established conditions can serve as a physical reason for the occurrence of stress serrations in jerky flow, .
After all the dislocation reactions are written down, one can obtain a set of three nonlinear differential equations for the process rates and to couple these equations to the machine equation representing the load sensed by the load cell. The final set of four differential equations looks as follows, :
In the last equation of the system (4) the coefficient K stands for some generalized effective compliance of the "sample+machine" system, is the imposed strain rate (i.e. the sum of the elastic and plastic strain rates) and - the Burgers vector. The second term in this equation is the plastic strain rate, which is proportional to the Burgers vector , the number of glissile dislocations Ng, and the velocity of gliding . A simple phenomenological power law can be taken in the first approximation to relate it to the effective stress applied to the material:
In our simulations we always assumed for simplicity that m = 2. The effective stress can be related to the applied stress in eqn. (6) in the following simple way, :
where in eqn. (5) is the stress required to induce velocity VO, H is a constant
characteristic of hardening and N is the total dislocation density ( ).
This is the Ananthakrishna model which has been successfully applied to the
description of the classical PLC-effect and even to such an interesting phenomenon as transition to temporal chaos on the stress vs strain curves in a certain interval of applied strain rates, which was confirmed experimentally later, . The last step in the analysis of the system of differential equations (4) is to reduce them to the dimensionless form by introducing the proper scaling of variables (it is much easier to solve numerically a system of equations in which all variables change in the interval close to unity, rather than differ by several orders of magnitude).
The new variables x, y, z, and are introduced as given below,
; ; ; (7)
and the original system (4) can be rewritten in the following scaled form, :
The contribution of the dislocations with atmospheres to the overall dislocation mobility was considered as negligibly small, and we assumed that in eqn. 8d, g = O
We are finally in a position to discuss the modifications that should be introduced into the model in order to describe temporal instabilities in fatigue.
Part 4. Modeling of Temporal Instabilities In Fatigued Metals
In order to determine the area in the parameter space which corresponds to the emergence of the limit cycle oscillations, a standard procedure is usually carried out called "the linear stability analysis", . According to this procedure, the stability of the thermodynamic (i.e. steady-state) branch of solutions is probed, for which the derivatives of x,y,z, and with respect to dimensionless time are equal to zero. In other words, for a linearized equation we look for the region in the parameter space (a, b, c, d, h, k, e, g ) in which eigenvalues 6~), of this linearized form around /~a are complex with Re ~l > 0, - this ensures that the trajectories in the phase space (x,y,z, ) will spiral out, [22,23]. The next step is to look for a surface S surrounding Pa into which all such trajectories enter. These requirements ensure a bounded variation of trajectories in the phase space often leading to a limit cycle. After this area in the parameter space has been established, the set of ordinary differential equations (8) can be integrated numerically for the values of parameters that guarantee the existence of oscillatory solutions, .
This program, clear and straightforward for the original system of equations (8), is unfortunately inapplicable to fatigued metals, because for cyclic loading steady-state solutions simply do not exist. In this situation tedious numerical experimentation remains the only tool available for the studying of such systems.
The modifications that must be introduced into the system (8) are obvious: instead of the constant strain rate, e, in equation (8d) one needs to put:
a). for pure cyclic loading (i.e. applied strain varies harmonically) we get for the strain rate: , where f is the effective loading amplitude and a) is the frequency;
b). for "creep fatigue", when a constant tensile loading is superimposed on the cyclic process (mostly at elevated temperatures), we obviously get:
c). finally, in order to study the "Neumann bursts" in a ramp-loaded metallic sample, the strain rate (the total strain changes as ) with the amplitude growing linearly with time, and .
There is one more property of the system that requires some attention. As it has been mentioned in Part 2 of the paper, because of the multi-scale behavior it is very difficult to integrate the system of equations (8) using conventional algorithms. The extremely small time step necessary in the parts of solution with temporal instabilities make the integration over multiple periods prohibitively slow. To overcome this problem, the 4-th order Rosenbrock algorithm with an automatically adjustable time step was used, .
Part 5. The Results of Computer Modeling
5.1 Harmonic Loading and Creep Fatigue We began by applying pure sinusoidal loading with various forcing frequencies and amplitudes to find a set of forcing parameters that qualitatively reproduced the phenomena observed in experiment. Figure 5a shows the time evolution of stress for a purely harmonic loading. Notice that in the first half period there is some indication on the stress curve of the stress "serrations", but all indications of activity on the stress curve are gone by the second period. Examining the mobile dislocation density as a function of time in Figure 5b substantiates this observation.
Figure 5. (a). Simulation data for stress vs time behavior of the model under purely harmonic loading; (b) Simulation data for mobile dislocation density vs time under the same conditions as (a). The following set of parameters has been used: (a, b, c, d, e, f, h, k, w) = (1.1, 0.012, 0.005, 0.0001, 0.0, 30.0, 0.2, 0.9, 0.0016) The initial densities of mobile, immobile, and Cottrell-type dislocations have been equal to 103 cm-2
After tedious numerical experimentation it has been found that adding a relatively small constant strain rate term to the sinusoidal forcing function caused stable and reproducible stress serrations that last practically indefinitely in time. This type of loading corresponds to the condition of creep fatigue. The first three periods of a "creep fatigue solution" are shown in Fig. 6. A comparison of the stress curve to the mobile dislocation density curve shows a very clear correspondence of peaks in dislocation density and peaks in stress: when stress goes down, the immobilized dislocations break free from the atmospheres of point defects (or any other pinning centers, ) and the mobile dislocation density undergoes a dramatic increase. This physically consistent result indicates that the modified Ananthakrishna model is capable of reproducing even the finest effects of dislocation dynamics in fatigued metals.
Figure 6. Computer simulation results for the case when a small constant strain rate was superimposed on the otherwise harmonic loading: (a). stress vs time, (b) mobile dislocation density vs time. The set of parameters: (a, b, c, e, f) = (0.7, 0.002, 0.008, 5.0, 20.0), the other parameters are as in Fig. 5.
Figure 7. Comparison of simulation data to experimental results. (a). Stress serrations in Cu-16at.%AI single crystal. Number of cycles applied to the sample N= 9750, total strain amplitude - 0.00322, ; (b). our simulation data for stress serrations (see also Fig. 6)
We attribute this constant strain rate requirement to the difference in the frequency spectra of the two different forcing functions, and . For a system of differential equations certain resonant frequencies exist (generally speaking, dependent on time, because the system is nonlinear) The frequency spectrum of pure sinusoidal loading function consists of only one frequency and thus it is very difficult to produce the complex resonance phenomena that, in our opinion might be responsible for the stress serration effect. By adding a small constant strain rate, we effectively expand the frequency spectrum to infinity, enabling the system to "pick out" the frequencies which are close to those at which it resonates. Thus, the stress serration effect can be explained as a non-linear resonance phenomenon in a system of non-linear ordinary differential equations. Analogous effects were observed and studied, mostly via computer experiments, in chemical engineering (in the theory of periodically controlled chemical reactors), .
One important aspect of the creep fatigue solution, however, is that it does not exhibit the serrations in compression that one would expect on the basis of experimental observations, Fig 4b. There is some evidence of serrations in compression if one examines carefully the dislocation activity vs time curve in Fig. 5a (at t=4000 and t=8000), but we could not model a pronounced effect.
The reason why it was not possible to obtain stress serrations in compression numerically might be connected with one fundamental limitation of our model. Indeed, temporal instabilities in fatigued metals can occur only at finite strains. In this situation, depending on the choice of coordinate system, it becomes necessary to introduce different strain measures (e.g., Cauchy-Green, or Almansi strain etc.) which coincide only in the case of infinitely small deformations, . Consequently, the time derivatives of the stress and strain tensors must satisfy the requirement of objectivity where A is the time derivative of Eulerian-objective tensor A, and is the lattice spintensor, . As a result, it is impossible to introduce objective stress or strain rateswithout considering the spatial coupling which was not taken care of in our model. This fundamental limitation might manifest itself in the form of symmetry-breaking behavior of stress serrations in tension and compression.
5.2 Ramp Loading An attempt was also made to reproduce the Neumann bursts in a ramp-loaded metallic system (see Part 4, case "c"). The bursts occur during the gradual increase of the stress amplitude from zero, Fig. 8a. Comparing our solution to the experimental result of Neumann, Fig. 8b, one can see how the dislocation density bursts obtained in numerical simulations correlate to the characteristic stress-strain response observed in experiment, .
Figure 8. Ramp loading: (a). Strain ("Neumann") bursts in a ramp-loaded Al single crystal, ; (b). mobile dislocation density vs time. Only data on the mobile dislocation behavior is given, because the behavior of immobile and Cottrell-type dislocations is highly correlated with (b).
One can attempt to link this result to the frequency spectrum of the forcing function, . The spectrum of the ramp loading function is broader than that of a pure sinusoid. This might be responsible for the initial resonance observed in our computer experiments. However, since the forcing function spectrum is not infinitely broad (as in the case of creep fatigue) the system still might have some difficulty "picking out" frequencies that match the conditions of resonance at any given moment of time. This might serve as a possible explanation of the deterioration of the Neumann bursts - a well-established experimental fact, .
Part 6. Conclusions
* An approach of non-linear dynamics has been used to model temporal instabilities (dissipative structures) in fatigued metallic alloys.
* Well-understood reactions between populations of dislocations and the assumption of homogeneous yielding were used to develop a mathematical model consisting of a system of nonlinear ordinary differential equations.
* Using this model we reproduced and explained the experimentally observed stress serrations under the conditions of 1). pure harmonic loading, 2). creep fatigue and 3). ramp loading, thus establishing for the first time the fundamental connection between the micro- and macro-mechanics of cyclic deformation.
I would like to express my deepest gratitude to my advisor, Prof. Campbell Laird, who helped introduce me to the exciting area of material fatigue, guided me when I was lost and gave me freedom to explore non-traditional approaches to old problems. Discussions with Professors John L.Bassani, Louis A.Girifalco, Charles McMahon, Vaclav Vitek and participation of the undergraduate NSF fellow David R.Williams (Cornell University) in part of this work during the summer of 1994 are also gratefully acknowledged.
1. G.Nicolis, I.Prigogine, Exploring Complexity. W.H.Freeman, NY (1989).
2. V. Vitek, Atomistic Studies of the Structure of Grain Boundaries and Dislocations in Proc. of the NATO Workshop on the Stability of Materials, to be published (1995).
3. P.Hirth and J.Lothe, Theory of Dislocations 2nd edition, Wiley Interscience, NY (1982).
4. A.P.Sutton and R.W.Balluffi, Interfaces in Crystalline Materials, Oxford University Press, Oxford (1995).
5. C.Laird, The Application of Dislocation Concepts in Fatigue: Dislocations in Solids (edited by F.R.N.Nabarro), North Holland, Amsterdam (1983) Chapter 26.
6. D. Hull and D.J.Bacon, Introduction to Dislocations, Pergamon Press, Oxford (1984).
7. J.W.Gibbs, The Collected Works--Vol. 1, Yale University Press, New Haven (1906).
8. Simulation Approach to Solids (edited by M.Ronchetti and G Jacucci) Kluwer Academic Publishers, Dodrecht (1991).
9. V. Vitek, Crystal Lattice Defects, v.5 (1975) 1.
10. M.S.DuesbeITy and Z.S.Basinski, Acta Metall. Mater., v.41 (1993) 643.
11. L.P.Kubin, Dislocation Patterning, in: Materials Science and Technology (a Comprehensive Treatment), ed. by R.W.Cahn, P.Haasen and E.J.Kramer, v.6, VCH, Weinheim (1993) p.137-190.
12. M.V.Glazov, L. Llanes and C.Laird, Phys. Stat. Sol, v.149 (1995) pp.297-321.
13. M.V.Glazov and C.Laird, Acta Metall. Mater. (1995)--in print.
14. G. Nicolis and I.Prigogine, Self-Organization in Non-Equilibrium Systems from Dissipative Structures to Order through Fluctuations). J.Wiley, NY (1977)
15. H. Haken, Advanced Synergetics, Springer-Verlag, Berlin (1983).
16. W.Rauchle, O.Vohringer, E.Macherauch, Mat. Sci.Eng., 12 (1973) 147.
17. S. I. Hong and C. Laird, Mat. Sci. Eng, A~24 (1990) 183.
18. P.Neumann, ActaMet., 17(1969) 1219.
19. W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery, Numerical Recipes in FORTRAN, The Art of Scientific Computing 2nd edn., Cambridge Univ.Press, NY (1992).
20. L.P.Kubin and Yu.Estrin, Acta Metall. Mater., 38 (1990) 697.
21. L.P.Kubin and Yu.Estrin, Phys. Stat. Sol. fB), 172 (1992) 173.
22. P.G.Drazin, Nonlinear Systems, Cambridge University Press, Cambridge (1992).
23. G.Ananthakrishna and M.C.Valsakurnar, J.Phys. D: ,4ppl.Phys. 15 (1982) L171.
24. G.Ananthakrishna, Scripta Metall. Mater., 29 (1993) 1183.
25. H.Neuhauser, Int. J. of Plasticity, 9 (1993) 421.
26. P.Rehmus, J.Ross, Periodically Perturbed Chemical Systems. In: Oscillations and Traveling Waves in Chemical Systems, ed. by R.Field and M.Burger, Wiley Interscience, NY (1985) pp.287-331.
27. C.Truesdell, The Elements of Continuum Mechanics, Springer-Verlag, Berlin (1985).
|Search||Student Member Application||Other Best Papers||Student Members' Page||TMS OnLine|