E L B 

 V

D H  C

 A  K

R H B F

PERIODIC SOLUTIONS OF A MULTI-DOF BEAM SYSTEM WITH IMPACT

The steady state behaviour is analyzed of a periodically driven multi-DOF beam system which has an elastic stop at its middle. The elastic stop is modelled in a continuous way by using the contact law of Hertz. The beam is modelled by using finite elements and subsequently reduced by using a component mode synthesis method. The steady state behaviour of the system reduced to one, two and four DOF's is investigated by calculating periodic solutions at varying excitation frequency. Periodic solutions are calculated by solving two-point boundary value problems by the multiple shooting method in combination with a path-following technique. It is shown that models with more than one DOF are required for a good assessment of the long-term behaviour of the system.

INTRODUCTION

In many practical engineering applications of mechanical systems, impacts at stops occur. It is important to carry out a dynamical analysis of such systems, to identify and subsequently reduce the noise and wear caused by repeated unacceptably large impacts. Examples of systems with impacts are gear rattle, heat exchanger tube wear in nuclear power stations and ships colliding against fenders. Systems with stops are typical examples of systems with strong local non-linearities. Although the non-linearity is local, the overall dynamic response of the system may change drastically. A system with stops cannot be linearized, so it is very difficult to predict the system response without a non-linear analysis.

In general, two different approaches are possible for analysis of the impact in the systems outlined above. One approach is to assume that the impact occurs instantaneously. The analysis is thus divided into two intervals: before and after the impact. The interconnection between the two intervals is made by a momentum balance and a parameter representing the amount of energy dissipation in the impact, called the coefficient of restitution. This approach is used in many analyses of impact oscillators and has led to the identification of bifurcational behaviour which is not found in smooth dynamical systems. Impact rules for such assumed instantaneous reversal of velocity induce discontinuities in the derivatives of maps defined by taking a section through the full phase space. In fact, these discontinuities yield the ''new'' bifurcational behaviour, called grazing or border-collision bifurcations [1][START_REF] N | Border-collision bifurcations: An explanation for observed bifurcation phenomena[END_REF][START_REF] C | Grazing bifurcations in impact oscillators[END_REF].

The other approach is to let the collision forces act in a continuous matter. The impact analysis of a system of colliding bodies is hence performed simply by including the collision of the contact forces in the system equations of motion during the contact period. A suitable model for the contact forces during the contact period is the non-linear Hertzian force-displacement law [START_REF] H | Schriften vermischten Inhalts[END_REF]. This model itself does not represent the energy dissipation process. Although the theory of Hertz is based on elasticity, some studies have been performed to extend the theory to include energy dissipation [START_REF] G | Impact: the Theory and Physical Behaviour of Colliding Solids[END_REF][START_REF] L | Continuous contact force models for impact analysis in multibody systems[END_REF]. In the paper by Foale and Bishop [7] bifurcations occurring in single-DOF impact oscillators in which the impact is modelled by using the Hertzian contact law are compared with ''grazing'' bifurcations which occur in impact oscillators where the impact is modelled in a discontinuous way.

The main focus of research into impact oscillators in recent years has been into single-degree-of-freedom systems. In this paper we will focus on the influence on the system response of adding more degrees of freedom (DOF's). A beam system is investigated which has an elastic stop at its middle and is periodically excited in the middle. A less complicated version of this beam system, supported in the middle by a massless one-sided linear spring, was investigated earlier by Fey [START_REF] F | Steady-state behaviour of reduced dynamic systems with local nonlinearities[END_REF], Fey et al. [START_REF] F | Nonlinearity and Chaos in Engineering Dynamics[END_REF] and Van de Vorst et al. [START_REF]  | [END_REF]. This beam system can be divided into a linear and a non-linear component. The linear component (the beam) is modelled by means of the finite element method and consequently has many more DOF's than the non-linear one. Because the numerical analysis of the resulting non-linear system is very expensive from a computational point of view, in particular for increasing number of DOF's, it is worthwhile to keep the number of DOF's as low as possible. This can be achieved by applying a reduction method to the finite element model of the linear component (the beam). The particular reduction method applied is the component mode synthesis method [START_REF] J | Combined Experimental/Analytical Modeling of Dynamic Structural Systems Using Substructure Synthesis[END_REF][START_REF] F | Steady-state behaviour of reduced dynamic systems with local nonlinearities[END_REF] and offers the possibility for a considerable reduction of the DOF's and for a systematic investigation of the influence of adding or deleting DOF's. Moreover, the component mode synthesis method can easily be used for more complex systems. The component mode synthesis method utilized in this paper uses free interface eigenmodes up to a cut-off frequency and residual flexibility modes to approximate the dynamic and static behaviour of the linear component. After reduction of the linear component, the non-linearity is added, resulting in a reduced non-linear system which will be valid for frequencies up to the cut-off frequency used in the reduction procedure. Because generically non-linear systems generate higher frequencies than their excitation frequency, the cut-off frequency has to be chosen much higher than the maximum excitation frequency.

The elastic stop is modelled by using the contact law of Hertz without energy dissipation, as mentioned earlier. Hence, the collision forces act in a continuous way. Periodic solutions of the reduced system are calculated for varying excitation frequency by solving two-point boundary value problems by the multiple shooting mehtod [START_REF] A | Numerical Solution of Boundary Value Problems for Ordinary Differential Equations[END_REF] in combination with a path-following method [START_REF] F | Steady-state behaviour of reduced dynamic systems with local nonlinearities[END_REF]. The response of the system is investigated for models with one, two and four DOF's. Also, different modal damping levels are applied. All calculations presented in this paper were carried out by using a development release of the finite element package DIANA [START_REF]DIANA User's Manual[END_REF].

BEAM SYSTEM WITH AN ELASTIC STOP AT ITS MIDDLE

The beam system which is analyzed numerically is shown in Figure 1. The beam is supported at both ends by leaf springs. If the displacement y of the middle of the beam positive, the beam hits a spherical elastic contact with radius r = 0•005 m, which is connected rigidly to the outside world. In the middle of the beam also a rotating mass is present generating a periodic force. Energy dissipation is established by a linear damper (not shown in Figure 1) at the middle of the beam (see also section 3).

The beam and leaf springs are modelled by using the finite element method and reduced by using the component mode synthesis method mentioned in section 1. Because of symmetry, only half the beam is modelled. (The elastic stop and the excitation force are positioned exactly at the middle of the beam. Hence the asymmetric model is linear, which means that symmetry breaking bifurcations cannot occur.) As mentioned in section 1, the contact force F s between the beam and the spherical elastic contact is modelled by using Hertz's law [START_REF] H | Schriften vermischten Inhalts[END_REF][START_REF] G | Impact: the Theory and Physical Behaviour of Colliding Solids[END_REF]: that is,

F s (y) = 6 k s y 3/2 , 0, y e 0 y Q 07 . ( 2.1)
In equation (2.1) the parameter k s is taken as 1•034 × 10 10 Nm 2/3 for a contact radius of r = 5 mm and E = 2•1 × 10 11 N/m 2 for the spherical contact.

SINGLE-DOF MODEL

After reduction of the finite element model, the equation of motion for the single-DOF model may be obtained as my¨+ by˙+ ky + k p y 3/2 = m e r e v 2 cos (vt).

(3.1)

Here

k p = 6 k s , 0, y e 0 y Q 07 and k s = 5•17 × 10 9 Nm 2/3 , v = 2pf e , b = 2jzmk, j = 0•0142, m = (8 T 1 M8 1 )/8 2 1y = 2•48 kg, k = (8 T
1 K8 1 )/8 2 1y = 16 877 N/m, m e = 1•0 g and r e = 4•92 mm. M and K are the mass and the stiffness matrices of the linear part of the system, i.e., the beam without the one-sided spring. 8 1 is the first eigenmode of the linear system with an eigenfrequency of 13•1 Hz. 8 1y is the element of 8 1 corresponding to the displacement of the middle of the beam y.

Branches of periodic solutions of this system have been calculated for varying excitation frequency f e . The maximum absolute displacements =y = max occurring in these periodic solutions are shown in Figure 2. Besides the harmonic resonance peak near 25•3 Hz, the figure also shows subharmonic resonance peaks. In Figure 3 the periodic solution at f e = 25•3 Hz is shown. Notice that, compared to the linear system, the first harmonic resonance peak has increased from 13•1 Hz to 25•3 Hz. Shaw and Holmes [START_REF] S | [END_REF] have given the following expression for the bilinear eigenfrequency f b of an undamped bilinear single-DOF system:

f b = {2z1 + a > (1 + z1 + a)}f 1 .
(3.2) Here a = k b /k, with k b the stiffness of the bilinear spring. In our system the contact forces are non-linear. If, however, a is calculated from equation (3.2), with f b = 25•3 Hz and f 1 = 13•1 Hz, this results in a 1 790 and by using this a the contact forces resulting from Hertz's law and from the bilinear system can be compared. In the bilinear case the contact force is F s b (y) = aky = 790 × 16 877y = 1•33 × 10 7 y N. In Hertz's law the contact force is F s H (y) = k s y 3/2 = 5•17 × 10 9 y 3/2 N. This means that F s b q F s H for 0 Q y Q 6•6 × 10 -6 and F s b Q F s H for y q 6•6 × 10 -6 . In our case (Figure 2) the penetration of the contact in the beam is y Q 4•5 × 10 -6 .

In Figure 2 it is shown that the influence of the stiffening character of Hertz's contact law is small because the high resonance peaks do not bend off much. Only in the top of the 1/3 subharmonic resonance peak near 76 Hz a small part of the 1/3 subharmonic branch is unstable, caused by the stiffening effect. The influence of the stiffening character of the contact law is low because the maximum penetration of the elastic contact in the beam does not vary very much compared to the maximum amplitude of the system, and the contact time of the beam and the elastic contact is low compared to the period of the periodic solutions.

The highest 1/2 and 1/3 subharmonic resonance peaks in Figure 2 are related to the harmonic resonance peak. The lower 1/3 (near 37•1 Hz) and 1/5 (near 62•2 Hz) subharmonic solutions are related to the second superharmonic resonance peak near 12•4 Hz. The 1/2 (near 16•1 Hz), 1/4 (near 33 Hz), 1/7 (near 58 Hz) and 1/10 (near 83 Hz) subharmonic solutions are related to the third superharmonic resonance peak near 8•1 Hz. The 1/3 (near 18•4 Hz), 1/5 (near 31 Hz), 1/9 (near 56 Hz) and 1/13 (near 81 Hz) subharmonic solutions are related to the fourth superharmonic resonance peak near 6•0 Hz (not visible in Figure 2).

On the harmonic branch and every subharmonic branch one or more periodic doubling routes leading to chaos exist. (These are not shown in Figure 2; the periodic doubling routes exist at the unstable parts of the branches marked by two flip bifurcations.) On the 1/2 subharmonic branch the qualitative dynamics change very much for small frequency variations near 40 Hz. This is caused by the fact that near 40 Hz grazing impact occurs: the beam just touches the elastic contact. For frequencies lower than 40 Hz in the 1/2 subharmonic solutions the beam hits the contact twice per period, whereas for frequencies larger than 40 Hz in the 1/2 subharmonic solutions the beam hits the contact once per period. Because the contact law is continuous, the system response does not change dramatically at grazing impact; no grazing bifurcations can be found in this system. However, as for discontinuous impacts, one would expect that the system response changes very much in the neighbourhood of grazing impact: i.e., if the penetration increases. Because of the grazing impact near 40 Hz, the maximum displacements of the system increase and two periodic doubling routes leading to chaos were found in that frequency area. Only the 1/4 subharmonic solutions of these period doubling routes are shown in Figure 2. Grazing impact also occurs at the 1/3 subharmonic resonance peak near 66 Hz.

TWO-DOF MODEL

By using the component mode synthesis method, the system can be reduced to two DOF's. Then, in the model one residual flexibility mode and one free-interface eigenmode ( f 1 = 13•1 Hz) are included. The two-DOF model has two eigenfrequencies: f 1 = 13•1 Hz and f 2 = 125•1 Hz. Notice that because only one free-interface eigenmode is included, the second eigenfrequency is inaccurate. Because the second eigenfrequency of the beam, f 2 = 117•2 Hz, lies close to the second eigenfrequency of the two-DOF model, still this model is suitable to show the effect of adding one DOF to the system on the system response. In the two-DOF model, besides the linear damper in the middle, extra damping is added by means of modal damping of the two eigenmodes with modal damping coefficient j m = 0•01. Extra modal damping is added because it is impossible to calculate periodic solutions without modal damping due to the low damping level of the system (without modal damping all periodic solutions are numerically unstable) and in practice also modal damping is present.

In Figure 4 are shown the maximum absolute displacements =y = max of the middle of the beam, occurring in the periodic solutions for varying excitation frequency of the two-DOF model. The figure shows again the harmonic, 1/2 subharmonic and 1/3 subharmonic solutions, now with extra resonance peaks. The extra resonance peaks on the harmonic branch are superharmonic resonance peaks of the second eigenfrequency of the model: a fifth superharmonic resonance peak (near 26 Hz), a fourth superharmonic resonance peak (near 32 Hz), a third superharmonic resonance peak (near 44 Hz) and a second superharmonic resonance peak (near 64 Hz). These superharmonic resonance peaks cause an additional resonance peak on the 1/2 subharmonic and 1/3 subharmonic branches. These subharmonic branches have three peaks, where the third peak is related to the fifth superharmonic resonance peak (near 26 Hz) of the second eigenfrequency of the model.

The vertical dashed lines in Figure 4 show the frequency at which the superharmonic resonance peaks of the second eigenmode are expected (for f 2 = 125•1 Hz). These dashed lines show that the fifth, fourth, third and second superharmonic resonances occur at a higher frequency than expected. The second harmonic resonance peak (not visible in Figure 4) occurs near 126•5 Hz, so the second ''eigenfrequency'' † increases from 125•1 Hz to 126.5 Hz. This does not explain the higher frequencies for which superharmonic resonances occur; the fifth superharmonic resonance peak occurs near 26•5 Hz and 26•5 × 5 = 132•5 Hz q 126•5 Hz. In a bilinear system with a = 790, which shows a similar shift of the harmonic resonance peaks as in our system (see section 3), the superharmonic resonances occur at the same higher frequencies. This means that the higher superharmonic resonance frequencies are not caused by the stiffening character of the contact law of Hertz. Investigation of the periodic solutions in the superharmonic resonance peaks shows that the closer the superharmonic resonance peak lies to the first harmonic resonance peak, the greater the penetration of the elastic contact in the beam and the lower the relative time per period for which the elastic contact is ''connected'' to the beam. Furthermore, the superharmonic resonance peaks which lie close to the first harmonic resonance peak are not dominated only by the second ''eigenfrequency'' but also by the first ''eigenfrequency''. This means that the closer a nth superharmonic resonance peak lies to the first harmonic resonance peak, the more the nth superharmonic resonance frequency will be increased from f 2 /n to (

f 1 + f 2 )/n.
Near the first harmonic resonance peak (23 Hz) two separated harmonic branches coexist ((a) and (b)). Both branches have a harmonic resonance peak. One branch (b) is unstable at the top. The unstable part on this branch is marked by two cyclic fold bifurcations and underneath the unstable part a stable branch (c) exists with a period doubling route to chaos. Investigation of the periodic solutions shows that in the upper (unstable) part of the branch (b) the beam hits the contact twice per period and in the lower (stable) part of the branch (c) the beam hits the contact once per period (see also the later Figure 6 in this section). Near the two cyclic fold bifurcations grazing impact occurs. At the top of the other (stable) harmonic branch (a) the beam hits the contact twice per period. Furthermore, the periodic solutions of the upper unstable part (b) are dominated by only the first ''eigenfrequency'' of the system, where the periodic solutions of the other branches (a) and (c) are dominated by both the first and the second ''eigenfrequency''. Apparently, because of the greater velocity of the beam at impact, the higher modes of the system are excited. As mentioned before, the superharmonic resonances near the first harmonic resonance peak occur at a higher frequency than expected. Because of this, the sixth superharmonic resonance peak also occurs at a higher frequency than expected. Since branch (b) is dominated by the first ''eigenfrequency'' only and branch (a) is dominated by both the first and the second ''eigenfrequency'', it can be concluded that the resonance peak on the stable branch (a) near 23 Hz is the sixth superharmonic resonance peak (( f 1 + f 2 )/n = (13•1 + 125•1)/6 = 23 Hz) and the resonance peak on the unstable branch (b) is the first harmonic resonance peak.

In the 1/2 and 1/3 subharmonic resonance peaks the beam also hits the contact twice per period, whereas in the single-DOF model the beam hits the contact once per period near the harmonic, 1/2 subharmonic and 1/3 subharmonic resonance peaks. In Figure 4 it is shown that near the 1/2 subharmonic and 1/3 subharmonic resonance peaks no separated branches coexist as they do near the harmonic resonance peak. However, near these subharmonic peaks similar behaviour is found as that near the harmonic resonance peak. With an increase of frequency from 45 Hz to 47 Hz on the 1/2 subharmonic branch (a), the maximum displacement increases until the top whereafter it decreases until the branch becomes unstable via a cyclic fold bifurcation. The branch (b) turns back and simultaneously the maximum displacement increases again. Near the cyclic fold bifurcation † Because the model is non-linear, one cannot speak of the eigenfrequencies of the model. Here, ''eigenfrequency'' is used for the frequency at which the model resonates.

grazing impact occurs; however, in both the stable (a) and the unstable (b) branch the beam again hits the contact twice per period. Further decreasing of the frequency results in a decrease of amplitude of the unstable branch which then becomes stable via a cyclic fold bifurcation at 44•5 Hz and turns back again with increasing frequency. In this stable branch (c) the beam hits the contact once per period (near 44•5 Hz grazing impact occurs). Further increasing the frequency leads to a resonance peak where the branch is unstable marked by two flip bifurcations. Here a period doubling route leading to chaos exists. Only the 1/4 subharmonic solutions of this period doubling route are shown in Figure 4. At the top of the 1/2 harmonic branch (c) at 47 Hz again grazing impact occurs, although for decreasing and increasing frequency in the 1/2 subharmonic solutions the beam hits the contact once per period.

The 1/2 subharmonic resonance peak near 47 Hz on branch (c) cannot be related directly to a superharmonic resonance peak. Investigation of the period solutions shows that again the stable 1/2 subharmonic branches (a) and (c) are dominated by both the first and the second ''eigenfrequencies'' while the unstable branch (b) (marked by two cyclic fold bifurcations) is dominated by the first ''eigenfrequency'' only. Near the 1/3 subharmonic resonance peak similar behaviour was found.

The subharmonic solutions related to the superharmonic resonance peaks of the first resonance peak, which dominated the response of the single-DOF model between the highest harmonic, 1/2 and 1/3 subharmonic peaks, could not be found in the two-DOF model. Only the 1/3 subharmonic solution related to the second superharmonic resonance peak of the first eigenmode could be found near 33 Hz and this branch is fully unstable now. Furthermore, now more period doubling routes leading to chaos are found on the harmonic and subharmonic branches, and in the frequency ranges 30-37 Hz and 59-66 Hz chaotic behaviour was found.

Notice that, in contrast to the frequency response of the single-DOF model, the 1/2 subharmonic and 1/3 subharmonic resonance peaks are now lower than the harmonic resonance peak. This is caused by the anti-resonance near 85 Hz. Furthermore, the frequencies at which the harmonic and subharmonic resonance peaks occur are decreased by 2 Hz. The reason for this is that the beam now hits the contact twice per period in the peaks, whereas the response at the peaks is dominated by the first two ''eigenfrequencies''. In comparison with the linear system, the first eigenfrequency is moved from 13•1 Hz to 23•1 Hz. The second eigenfrequency has not changed very much (125•1 Hz:126•5 Hz).

In Figure 5 are shown the maximum displacements occurring in the periodic solutions if the model damping coefficient is decreased to j m = 0•001. The response of the system is similar to the response for j m = 0•01. However, now the resonance peaks are higher. Nevertheless, in the highest resonance peaks the beam still hits the contact twice per period. Furthermore, the unstable parts on the branches (most of them with period doubling routes) have become larger. Again, the subharmonic solutions caused by superharmonic resonance peaks of the first harmonic resonance peak could not be found. One additional 1/14 unstable subharmonic branch was found near 57 Hz.

The periodic solutions of the two-DOF model at f e = 23•0 Hz are shown in Figure 6. The figure shows that the periodic solutions of branches (a) and (b) hit the contact twice and that, in comparison to the periodic solution of the single-DOF model (Figure 3), in the periodic solutions of the two-DOF model higher frequencies are present.

FOUR-DOF MODEL

In this section the system is reduced to four DOF's. Simulations with the system reduced to five DOF's showed no major differences, so it may be concluded that the four-DOF In Figure 7 are shown the maximum absolute displacements =y = max of the middle of the beam occurring in the periodic solutions for varying excitation frequency of the four-DOF model with j m = 0•01. Again, the response is dominated by the harmonic resonance peak and 1/2 and 1/3 subharmonic resonance peaks which are related to the first harmonic resonance peak. Also again, superharmonic resonances exist at the harmonic branch which are related to the second and third ''eigenfrequencies'' of the system and again these superharmonic resonance peaks occur at a higher frequency than expected just as in the two-DOF model. These superharmonic resonance peaks result in additional resonance peaks on the 1/2 and 1/3 subharmonic solutions. The first (highest) peak on the branches with subharmonics is related to the first harmonic resonance peak. The other two peaks on the 1/2 and 1/3 subharmonic branches are dominated by both the second and third ''eigenfrequencies''. In the previous section we saw that every resonance peak near the first harmonic resonance peak is projected on the 1/2 and 1/3 subharmonic branches. If the extra two peaks on the 1/2 and 1/3 subharmonic branches are projected back on the first harmonic resonance peak, these two peaks have to be caused by the fifth superharmonic resonance peak near 27 Hz of the second ''eigenfrequency'' (see the inset of Figure 7). Because of superharmonic resonance of the third ''eigenfrequency'', the fifth superharmonic resonance peak is divided into two peaks: one peak is high and one peak is much lower (with two peaks on it). Because the fifth superharmonic resonance peak is divided into two peaks, also two extra peaks occur on the 1/2 and 1/3 subharmonic branches.

In contrast to the two-DOF model, in the four-DOF model no second harmonic branch was found near the first harmonic resonance peak. Still, there might be a second harmonic branch but the 1/2 and 1/3 subharmonic resonance peaks do not indicate that one exists. Also, because in the four-DOF model the second eigenfrequency is moved from 125•1 Hz to 117•5 Hz, in comparison with the two-DOF model, now the frequency for which the sixth superharmonic resonance peak occurs has also decreased. In the two-DOF model the sixth superharmonic resonance peak was responsible for the existence of the second harmonic branch. The harmonic, 1/2 and 1/3 subharmonic resonance peaks are stable and the beam hits the contact twice per period. In the two extra peaks on the 1/2 and 1/3 subharmonic solutions the beam hits the contact once per period.

In Figure 7 it is shown that fewer period doubling routes on the harmonic, 1/2 and 1/3 subharmonic branches exist than in the two-DOF model with j m = 0•01. This is caused by the fact that because now four eigenmodes are damped with modal damping coefficient j m = 0•01, and hence the amount of damping in the system is increased which results in less non-linear phenomena. Nevertheless, in the frequency areas 30-38 Hz and 58-63 Hz chaotic behaviour was found.

The response of the system if the modal damping coefficient is decreased to j m = 0•001 is shown in Figure 8. Compared to Figure 7, now all sub-and superharmonic resonance peaks increase and in the frequency domain in which the harmonic resonance peak exists many superharmonic resonance peaks can be seen. These superharmonic resonance peaks are not related to one ''eigenfrequency'' of the model. The superharmonic resonance peaks in the frequency domain 24-26 Hz are again dominated by the second and third ''eigenfrequency'' of the model, and they are responsible for the two extra resonance peaks on the 1/2 and 1/3 subharmonic branches. Again in these peaks the beam hits the contact once per period. The higher peaks in the frequency domain 22-23 Hz are dominated by the first and second ''eigenfrequencies'' of the model. Some peaks are unstable and here period doubling routes leading to chaos exist. In the stable peaks the beam hits the contact once per period, in the unstable peaks the beam hits the contact three times per period (see Figure 9). Also on other parts of the branches additional period doubling routes leading to chaos exist than for j m = 0•01.

CONCLUSIONS

Periodic solutions have been calculated for a multi-DOF beam system with an elastic stop at its middle for varying excitation frequency. The linear beam was modelled by using finite elements and reduced by using component mode synthesis. The elastic contact as modelled by using the contact law of Hertz and the periodic solutions were calculated by solving two-point boundary value problems by the multiple shooting method. By using this method it is possible to calculate stable and (very) unstable periodic solutions and in combination with a path-following method very complex branches of periodic solutions can be calculated.

The results show that the behaviour of the single-DOF system can be globally predicted in advance if the frequency for which the harmonic resonance peak occurs is known. All the superharmonic and subharmonic resonance peaks can be predicted by using the harmonic resonance frequency. However, because on every harmonic and subharmonic branch period doubling routes leading to chaos exist, the behaviour of the single-DOF system is dominated by chaotic behaviour between the highest harmonic, 1/2 and 1/3 subharmonic resonance peaks.

The two-DOF model shows globally the same maximum amplitudes as the single-DOF system. However, because of the second eigenmode an additional peak is found on the 1/2 and 1/3 subharmonic branches. In contrast to the single-DOF model, where the beam hits the elastic contact once per period in the highest resonance peaks, now the beam hits the elastic contact twice per period in these peaks. This means that the periodic solutions near the resonance peaks are very different from those of the single-DOF model. Between the highest harmonic, 1/2 and 1/3 subharmonic resonance peaks the behaviour is again dominated by chaotic behaviour, which is now caused by additional period doubling routes on the harmonic, 1/2 and 1/3 subharmonic resonance branches. The subharmonics caused by superharmonic resonances of the first ''eigenfrequency'', which dominated the response of the single-DOF system, have disappeared.

The response of the four-DOF model shows that although the third linear eigenfrequency of the beam is four times the highest excitation frequency under investigation, this ''eigenfrequency'' still has a large influence on the system behaviour. The results show that the amount of damping in the model has a large influence on the frequency responses. Increasing the modal damping leads to fewer period doubling routes and smaller frequency intervals with chaotic behaviour. However, also at higher damping levels 1/2 and 1/3 subharmonics dominate the frequency response and the influence of the third linear eigenfrequency is still large.

Because the impact contact has been modelled using a continuous contact law, no grazing bifurcations were found. In the neighbourhood of grazing impact, i.e., the beam just touches the elastic contact, in many cases period doubling routes occur. However, this is not always the case, and not every period doubling route is in the neighbourhood of a frequency for which grazing impact occurs.

Figure 1 .

 1 Figure 1. The beam system. E = 2•1 × 10 11 N/m 2 , r = 0•005 m, h = 0•01 m, hb = 0•001 m, b = 0•09 m, bb= 0•075 m, r = 7746 kg/m 3 , rb = 7713 kg/m 3 .

Figure 2 .

 2 Figure 2. The maximum displacements of the periodic solutions of the 1-DOF model. • • • • , Unstable; --, stable.

Figure 3 .

 3 Figure 3. The periodic solution of the 1-DOF model at fe = 25•3 Hz.

Figure 4 .

 4 Figure 4. The maximum displacements of the periodic solutions of the 2-DOF model; jm = 0•01. • • • •, Unstable; --, stable.

Figure 5 .Figure 6 .

 56 Figure 5. As Figure 4, but jm = 0•001.

Figure 7 .

 7 Figure 7. The maximum displacements of the 4-DOF model; jm = 0•01. • • • •, unstable; --, stable.

Figure 8 .

 8 Figure 8. As Figure 7, but jm = 0•001.

Figure 9 .

 9 Figure 9. The periodic solution of the 4-DOF model at fe = 22•5 Hz; jm = 0•001.

ACKNOWLEDGMENT

This work was supported by the Centre for Mechanical Engineering of the TNO Building and Construction Research, Delft, The Netherlands.


J M T Thompson

R Ghaffari

Chaotic dynamics of an impact oscillator

An impact oscillator is shown to exhibit complex dynamics. Its resonance response contains regions where, after an infinite cascade of period-doubling bifurcations, chaotic motion typical of a strange attractor is observed. The regions are bounded on both sides by subharmonic resonances. Quantitative agreements are obtained with the Feigenbaum scenario of chaos. This novel feature of a substantial marine technology program may be of general cross-disciplinary interest to mathematicians and physicists.

There is much current interest in physics in chaotic motions governed by strange attractors. We present a new example of a system exhibiting such phenom- ena, and make quantitative comparisons with a scenario of period doubling as a route to chaos. 'T here is good agreement, supporting the hope that the behavior of many nonlinear dynamical systems can be understood by means of a universal model.

The present work arose in marine engineering, where the slackening of a mooring can introduce a discontinuity in stiffness of, for example, an articulat- ed oil-loading tower. ' Such systems can sometimes be modeled as an impact oscillator, that rebounds elastically whenever the displacement Xdrops to zero, subjected to sinusoidal forcing. This simple system is of wide interest beyond marine technology.

The nondimensionalized equation of motion of our impact oscillator can be written as' X+ -X+ 2 X= -2sinr, X &0 4q where a dot denotes differentiation with respect to time r, i; is the damping ratio defined with respect to the effective natural frequency of the unforced and undamped rebounding system, and q is the ratio of the forcing frequency to this effective frequency.

Defining y as half the maximum value of X, the reso- nance response curve has been determined as y(7i)

in Fig. 1 for the fixed damping ratio shown. Well- defined resonant peaks are seen, corresponding to a fundamental response ( n = 1) at q = l and subhar- monic resonances of order n =2, 3, 4, . . . at q = n.

The areas shown in Fig. 1 correspond to solutions with one impact per response cycle, but between the peaks more than one impact is observed. A precise digital computer program' has been used to explore carefully the region between two adjacent peaks, us- ing the Poincare mapping points (X,X) at r equal to multiples of 2n. A cascade of period-doubling bifur- cations is observed (Fig. 2) leading to a chaotic solu- tion at q =4.5. Notice that the chaos is located between an n =4 cascade deriving from the n = 4 resonance, and an n = 5 cascade originating from the n =5 resonance.

This approach to chaos via period doubling is in line with the Feigenbaum scenario which draws on the behavior of a universal quadratic map. The bi- furcations of this map have the property that the con- secutive control parameter intervals tend to a fixed ratio of 5 =4.6692. . . as an accummulation point is approached.

In our problem, we define q"asthe range of q over which the subharmonic of order n is observed. We then find from a refined version of Fig. 2 that qs/gi6=4 56, g&6/q32 469 and 7i32/q64=4. 64.

These agree very well with Feigenbaum's number, which relates strictly to the limit as n tends to infinity, A standard test for the chaotic motions of a strange attractor is that solutions from adjacent starts should diverge exponentially until they become completely uncorrelated. This behavior for the impact oscillator riou ing i urcations leading to chaos at g=4.5.

is shown in Fig. 3 for three starting separations on a plot of -log1o 8 against steps of Ar =2m. . The noisy straight lines confirm that 8 varies as AoN', where i is the number of steps and N is the Liapunov number.

From a set of graphs similar to Fig. 3, we estimate the value of N =1.17 for our impact oscillator. This compares with N =1.26 for regions of the quadratic map, and N =1.52 for the Henon strange attrac- tor 7 11

We note, finally, that we have observed the coexistence of an n =10 subharmonic with a presumed strange attractor at a value of q =4.55 as shown in Fig. 4. The solution obtained here depends only on the starting conditions of the time integration. This discovery of period-doubling bifurcations, chaos, and strange attractors in the resonance of a simple impact oscillator may be of interest to en- gineers, physicists, and mathematicians alike. Full details of the marine study are planned to be published elsewhere. ' -tog R-

  FIG. 1. Resonance response curve for the impact oscilla- tor.

  FIG. 2. Sequence of steady-state Poincare maps showing period-d bl' b'f

  . Divergence study showing a noisy exponential growth of the separation between close starts.

FIG. 4 .

 4 FIG. 4. Attractors of two coexisting multiple solutions for a fixed system. The attractor observed depends on the starting conditions.


J C Simo

N Tarnow

The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics 1

Dedicated to K. Kirchgiissner on the occasion of his 60th birthday

Lagrangian description

§ 1. Introduction and motivation

Classical nonlinear continuum mechanics is built on the balance laws of linear and angular momentum, supplemented by balance of energy and a statement of the second law. These two additional laws are satisfied if the internal dissipation in the material ~~nt, understood as the difference between the stress power and the rate of change of the internal energy, is positive. 3 In particular, for nonlinear elastodynamics with conservative loading, the internal dissipation vanishes for smooth solutions and the total energy in the system E~ot is exactly conserved. For non-smooth solutions; i.e. in the presence of shocks, dissipation arises and the energy in the system is no longer conserved; see [START_REF] Knowles | On the dissipation associated with equilibrium shocks in finite elasticity[END_REF]. In what follows, we shall only be concerned with smooth solutions.

In the absence of external loading or in the presence of symmetries, the dynamics of a continuum exhibits two fundamental conservation laws in the pure traction initial boundary value problem (IBVP): The conservation laws of total linear momentum L 1 and total angular momentum 1 1 • In addition, the condition of positive dissipation yields an a-priori estimate on the rate of change of the total energy, which decays exactly by minus the dissipation. In summary: 3

Implicit here is the assumption that the entropy production due to heat conduction is positive; equivalently, the heat flux does not oppose the temperature gradient; see [START_REF] Truesdell | The Nonlinear Field Theories of Mechanics[END_REF].

These fundamental properties of the continuum dynamics play a key role both in purely mathematical aspects of the IBVP as well as in many engineering applications. For instance, the realization of attitude controls for orbiting flexible satellites exploits these properties in a crucial manner.

Temporal and spatial finite difference/finite element discretizations of the continuum dynamics need not, and in general will not, inherit the conservation of momentum properties ( 1.1 ) 1 , 2 and the a-priori estimate ( 1.1 ) 3 • The central goal of this work is the systematic construction of algorithmic approximations to the continuum dynamics that will inherit, by design, these conservation properties and the a-priori estimate. There is a sizable body of literature concerned with the design of algorithms possessing these conservation properties; see e.g. [START_REF] Bayliss | How to make your algorithm conservative[END_REF] and Labudde and Greenspan [1976a,b]. It is often the case, however, that algorithms designed to preserve the a-priori energy estimate result in violation of momentum conservation and vice versa. For instance, the conservation form of the mid-point rule is an exact momentum conserving algorithm which does not conserve energy for autonomous Hamiltonian systems, except in the linear regime; see Simo, Tarnow and Wong [ 1992]. On the opposite side, the method of [START_REF] Hughes | Transient finite element formulations that preserve energy[END_REF] yields an energy conserving algorithm which does not preserve angular momentum. For rigid body mechanics, on the other hand, [START_REF] Simo | Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and angular momentum[END_REF] show that the conservation form of the mid-point rule in body coordinates, together with the exponential mapping in the rotation group, define an exact energy and momentum conserving algorithm.

For smooth solutions in nonlinear elastodynamics, the model problem of interest here, the algorithmic treatment proposed below achieves exact conservation of both energy and total linear and angular momentum. Remarkably, the scheme amounts to a simple modification of the conservation form of the mid-point rule, according to the following prescription:

Step 1. Write a mid-point rule approximation to the IBVP for nonlinear elastodynamics, written in conservation form, with the symmetric stress tensor left open. We show that this step defines a 6-parameter family of exact momentum conserving algorithms, regardless of the specific form of the constitutive equation ..

Step 2. For elastodynamics define the convected symmetric stress tensor by evaluating the gradient of the stored energy function at a convex combination of the right Cauchy -Green tensors at the initial and final times of a typical time increment.

From a mechanical standpoint, Step 1 is equivalent to establishing balance of linear momentum (in conservation form) on the mid-point configuration defined as the average of the initial and final configurations in a time step.

It is somewhat surprising that global balance of momentum is satisfied for any symmetric stress tensor. For the Saint Venant-Kirchhoff model Step 2 has this remarkably simple interpretation: Evaluate the elastic constitutive equation with the average Lagragian strain tensor and not with the Lagrangian strain tensor of the average configuration. We remark that convex combinations of the right Cauchy-Green tensor make geometric since positive definiteness is preserved and, more importantly, all the tensors involved live in the same (reference) configuration.

Although nonlinear elastodynamics is an infinite dimensional Hamiltonian system, in the sense described in Simo, Marsden and Krishnaprasad [ 198~], the algorithmic treatment proposed herein does not rely on this Hamiltonian structure and is also applicable to dissipative systems. In fact, in contrast with the conservation form of the mid-point rule, it can be shown that the energy-momentum algorithms described herein do not define a symplectic transformation in phase space. See Arnold [ 1988] for an explanation of this terminology and De Vogelaere [1956], Feng [START_REF] Kan | Difference schemes for Hamiltonian formalism and symplectic geometry[END_REF], [START_REF] Zhong | Lie -Poisson Hami/ton -Jacobi theory and Lie -Poisson integrators[END_REF], [START_REF] Lasagni | Canonical Runge-Kutta methods[END_REF], [START_REF] Sanz-Serna | Runge-Kutta schemes for Hamiltonian systems[END_REF] and Simo, Tarnow and Wong [ 1992] for the algorithmic implications. Nevertheless, the property of exact energy conservation induces a natural notion of algorithmic stability not present in symplectic schemes such as the mid-point rule. In fact, as demonstrated by our numerical simulations, both the mid-point and trapezoidal rules may produce numerical solutions exhibiting blow-up for sufficiently large step-sizes, while the proposed scheme renders solutions that remain stable. §2. The continuum problem. Conservation laws

We summarize below basic results on continuum mechanics needed for our subsequent algorithmic developments presented in Section 3. This methodology is by no means restricted to the ( dissipationless) Hamiltonian case, as exemplified in the continuum setting by nonlinear elastodynamics. In fact, our construction of exact momentum preserving algorithms is independent of the specific form adopted by the constitutive equations.

In what follows 91 c ~ndim will denote the reference placement occupied by a continuum body where 1 < ndim < 3 is the spatial dimension. The set Pi is assumed open, bounded, with smooth boundary oPA. In agreement with a Lagrangian description of the motion material points are labeled throughout by their position X E 91 in the reference placement.

Some basic notation

Let the reference boundary of!A be partitioned into disjoint subsets such that ofJl = r Q) u r u with r Q) n r u = 0. The (smooth) manifold Q of admissible configurations is then defined by (2. [START_REF] Hughes | Transient finite element formulations that preserve energy[END_REF] where iP is specified data on r cp• For elastodynamics one typically assumes that Q c [Hs(BI)]ndim with s > ndim/2

+ 1. A motion in some time interval ~ is a C 1 map t E ~ ~ tp 1 E Q, with associated material velocity field V 1 = q,t and canonical momenta n 1 =eo V 1 ,
DqJtStN = t; on r 0" X ~' (2.3)
where N: o!!J ~ S 2 is the unit normal field to the reference boundary of!A and t;: r 11 x ~ ~ !R.ndim is the prescribed nominal traction vector. In addition to the boundary tractions it will be assumed that the body is acted upon by a body force ft : PJ x ~ -t !R.ndim.

The weak form of balance of momentum

Let "// c [Hs(PJ)]ndim denote the space of test functions associated with admissible configurations in Q, defined in the usual fashion as (2.5c) Equations (2.5a,b,c) define• a nonlinear initial boundary value problem for the motion t E 0 H qJ 1 E Q of the continuum system once a constitutive equation relating the convected stress tensor St to the motion (/) 1 is specified.

Constitutive equations. Positive internal dissipation

As pointed out above, the algorithmic approach described below encompasses a general class of constitutive relations, not necessarily restricted to elastodynamics, subject to the only requirement of positive dissipation. This restriction leads naturally to meaningful definitions of nonlinear algorithmic stability.

Energy balance: Positive dissipation

Let K and gp~nt denote the kinetic energy and the stress power of the continuum system, respectively, as defined by the relations Kt= ~(no Vt) = ~!Ieo Vt 11 2 and gp~nt = (DqJrSo grad[Vt]).

(2.6)

We further denote by _gp~xt the expended power of the external loading, defined by the expression _qp~xt := (f;, V,)+(!;, vt )r t E 0.

(2.7) For fixed time t E ~ the material velocity field defines an admissible test function. Therefore, setting 1J = V 1 in (2.5) and making use of the notations above yields the result dKt .

-+ ,q;mt = _gpext for any t E 0,

dt t t '
(2.8) known as the theorem of expended power.

Now let

E~nt := JBl e df!J be the internal energy function which specifies the stored energy in the continuum system and define the internal dissipation function !0 ~nt of the system as the difference between the total stress power ,q;~nt and the rate of change of the internal energy. We restrict our attention to constitutive equations obeying the following reduced dissipation inequality . .

dE~nt gmt ==.?Jmt- > 0 for any t E 0 t t dt '
(2.9) which (in the absence of heat conduction) can be viewed as an equivalent statement of the Clausius-Duhem form of the second law. Combining (2.8) and (2.9) yields the following result for the rate of change of the kinetic and internal energies:

(2.10)

This relation plays a crucial role in our subsequent algorithmic developments.

Constitutive equations: Examples

Two representative and fairly classical examples which fall within the preceding framework are: i. Nonlinear elastodynamics. Here the internal energy is defined locally in terms of the deformation gradient D(/) 1 of the motion by means of the frame-invariant stored energy function e = e(Dlp i Dlf't). The convected stress tensor st and the internal dissipation then become st =2Ve(DqJ'{Dlpt) and .@~nt _o Yt E ~-

(2.11)

The Cauchy stress tensor a, is obtained from (2.11) 1 via the Piola-transformation (2.2). The fact that .@~nt vanishes for an elastic material justifies the denomination 'internal dissipation' assigned to .@~nt. For the incompressible problem, obtained by appending the constraint J(qJ,) = 1 in n x ~' the Cauchy stress tensor is defined up to a spherical part pl involving the hydrostatic pressure p. Nonlinear elastodynamics possesses a well-understood Hamiltonian structure; see Simo, Marsden and Krishnaprasad [ 1988] for detailed exposition.

ii. Incompressible Navier-Stokes system. Let Vt == vt 0 ( / ) ; 1 be the spatial velocity field with spatial gradient

Vv 1 = (grad[ V 1 ]Dlf'; 1 ) o qJ; 1
• In the linear incompressible viscous model the Cauchy stress is defined by the constitutive relation

a 1 = -pl + 2}1 sym[Vv 1 ] with J(qJ 1 ) = 1 <=> div[v 1 ] = 0 Vt EL (2.12)
The convected stress tensor S 1 is obtained from (2.12) via the Piola-transformation (2.2). The internal dissipation function is

~~nt = f1 J. Vv 1 • Vvt dx > 0 for Jl > 0, tpl(fJI)
(2.13) while th'e internal energy vanishes identically; i.e., E~nt = 0.

Several classical models in nonlinear solid mechanics exhibit positive dissipation while retaining a non-vanishing stored energy function; 1.e., finite plasticity.

The unforced problem: Balance laws and a-priori estimate

A number of important applications involve the solution of the initial boundary value problem with prescribed zero loading and driven, therefore, solely by initial conditions; i.e., the unforced problem. This situation is also relevant to the nonlinear stability analysis of algorithmic approximations (see [START_REF] Simo | Nonlinear stability of the time discrete variational problem in nonlinear heat conduction plasticity and viscoplasticity[END_REF]) and will be considered in what follows.

Conservation laws and a-priori estimate

Recall that the resultant force F~xt and the resultant moment T~xt of the prescribed external loading (relative to some inertial frame) are defined by the expressions F~xt := f j; df!J + I t; dr and r~xt == r CfJt X j; df!J + I f/Jt X ~ dr.

Je~

Jr a

Je~

Jr a

(2.14)

The following assumptions specify the framework relevant to the analysis of algorithms designed to inherit conservation laws present in the continuum problem:

(A1) Assume that p~xt = 0, r~xt = 0 and q;~xt = 0, Vt E 0.

(A2) Consider pure Neumann data; equivalently, assume r q> = 0.

From a physical standpoint these two conditions can be realized by considering the motion of a continuum body subject to no boundary restrictions, with prescribed zero loading, driven by specified initial conditions. If Lt and 1 1 denote the total linear and total angular momentum defined by the standard expressions L, •= f .. n, dt4 and J, •= L qJ, x n, dt4, (2.15) under assumptions (AI) and (A2) above the ensuing motion possesses two fundamental conservation laws and one a-priori estimate:

dLt = O dJ 1 = O dt ' dt (2.16)
The a-priori estimate (2.16) 3 follows at once from relation (2.10) since q;~xt = 0. Relations (2.16) 1 , 2 are the familiar conservation laws of total linear and total angular momentum, respectively. These latter two results can be derived directly from the weak form (2.5) of the momentum equation by the argument given below. The same approach is exploited in the algorithmic analysis of the momentum conserving properties presented in Section 3 .2.

Conservation of momentum and the weak form

Since r fP = 0 by assumption, it follows that test functions in "Y are unrestricted on the boundary. In particular, two specific class of test functions are now contained in "Y:

i. Infinitesimal translations, defined by q(X) = ~' 'if X E ~'
where e E rR_ndim is a fixed vector. For this class of test functions we have grad[q] -0. Expression (2.5) along with definitions (2.14) 1 and (2.15) 1 give (2.17) which implies F;xt -L 1 = 0 since ( 2.17) holds for any e E rR 3 and collapses to (2.16) 1 since F~xt = 0.

ii. Infinitesimal rotations (ndim = 3), defined by '1 = e x (/Jc where e E rR 3 is a fixed vector defining the rotation axis. Following a standard convention let ( denote the associated skew-symmetric matrix defined by the relation (a:= e x a for any vector a E IR 3 • For this class of test functions a direct calculation gives grad[q] = (DqJ 1 • Inserting this result in (2.5) and using (2.14) 2 along with (2.15h yields, after a straightforward manipulation, the result

(DqJ 1 Sn (D(/Jc) = (cpt xft, ~) + (cpt X 1;, ~)r-(({Jt X Ttn ~) = [T;xt -jt] , ~.
(2.18)

Since the Cauchy stress tensor ac is symmetric, from the Piola transformation (2.2) we conclude that Dcp 1 S 1 DqJ T is also symmetric and, therefore,

"" r "" _ • "" ""r (DtptSn enlf't > = (DtpcScD(/J t ' ~ > = 0 Since ~ = -' . (2.19)
Since e E ~3 is arbitrary combining (2.19) and (2.18) yields r;xtj = 0 which implies (2.16) 1 since T~xt = 0. §3. Time discretization. Conserving schemes

The dynamics generated by a number of important continuum system, including nonlinear elastodynamics, finite plasticity and the Navier-Stokes system among many others, obeys the conservation laws (2.16) 1 , 2 anda-priori estimate (2.16) 3 provided assumptions (Al) and(A2) hold. Our main goal in this Section is the design of algorithmic approximations which will inherit these crucial properties of the continuum dynamics.

The algorithmic framework. Conserving schemes

Consider a partition fl = U~=o [tn, tn+ d of the time interval fl of interest, focus attention on a typical time interval [tn, tn + d and let ilt == tn + 1 -tn > 0 denote the corresponding time-step size. Following a standard convention, we shall denote by either ( • )n or ( • )n + 1 the algorithmic approximations at times tn and tn+ 1 to the continuum (time-dependent) variable (•)t. In our algorithmic analysis we use the following restatement of condition (Al) of Section 2 (Al *) For the unforced problem, assume that F~xt = 0, T~xt = 0 and

q>~xt = 0, for all t E [tn, tn + tl•
We assume all the variables of interest specified at time tn; in particular, the initial data (cpn, nn) and the symmetric Piola-Kirchhoffstress .. c • (X), for all X E ~. Our goal is the design of a consistent approximation to lile IBVP defined by (2.5a,b,c), along with a suitable constitutive equation of the type described in Section 2.2, possessing the following properties whenever assumptions (A 1 *) and ( A2) hold: i. Discrete momentum conservation in the sense that (3.la)

ii. Discrete energy conservation, in the sense that

E tot < Etot d Etot = Etot I'M int = 0 _ 0 1 N n + I - n an n + 1 - n ~ ;;z; t - , n -, , • • • , • (3.1b)
The requirement E~o~ 1 = E~ot whenever .@~nt = 0 in D corresponds to exact energy conservation in the Hamiltonian case and will be the only situation addressed below within the specific context of elastodynamics. Algorithmic. approximations possessing the properties in (3.1) will be referred to as energy-momentum conserving schemes in what follows.

For nonlinear elastodynamics, the design of a conserving scheme is . accomplished in two steps. First, we identify a class of one-step implicit algorithms which, remarkably, achieve exact momentum conservation with total independence of the approximation used for the convected stress tensor. In fact, the property of exact momentum conservation holds for any constitutive model. Second, we specialize this class of algorithms to nonlin-ear elasticity and show that an approximation to the constitutive equation (2.11) 1 can always be constructed such that exact conservation of energy is achieved. From a geometric standpoint this two-step construction is analogous to the approach recently advocated in Simo, Tarnow and Wong [ 1992]. The actual definition of the two steps, however, is entirely different.

A class of exact momentum conserving algorithms

Consider a one-parameter family of configurations (/)n + rx E Q, with associated velocity field Vn +a., defined by the convex combination

+ <Dcpn + rxS, grad['l]) = lfn+rx' tJ) +<~+a' tJ)r VtJ E "//, (3.3a)
where S = sr is a completely arbitrary symmetric tensor field defined on 86.

The algorithmic weak form (3.3a) is supplemented by the local relations

(3.3b)
which cotnplete the algorithmic approximation to equations (2.5a,b). We emphasize that the symmetric tensor S is regarded here as defining a 6-parameter family of algorithms to be further specialized below. Our objective is to examine the conservation properties inherited by the time discretization (2.5a,b ).

Discrete momentum conservation

Let p~x~ iX and T~x~ a denote the total resultant force and resultant moment defined by (2.14) and evaluated at time tn+ 1 • We show below that ifS is symmetric but otherwise completely arbitrary the algorithmic approximation defined by (3.3a,b) yields the following relations

Ln + 1 -Ln = ~tF~x! a for any a E (0, 1],} J n + 1 --Jn = ~t T~~ rx only if a = ~ . ( 3
'•(In+ I -In)=<,, (/)n+ I X 1tn+ I -(/)n X 1tn) = ( ~' lfJn + I X (1tn + I -1tn) + ( lfJn + 1 -lfJn) X 1tn)) = ( c;, CfJn + a X (1tn + I -1tn)) + ( ,, (lfJn + 1 -lfJn) X (ct1tn + ( 1 -a)1tn + 1) ).
(3.5)

The first term of the last right hand side of (3.5) can be evaluated as follows.

We observe that '1 =; x lfJn + a is an admissible test function in "// for arbitrary (constant)~' with gradient grad[q] = (DlfJn+a• In view of (3.3) and the symmetry of S we conclude that ( ~' lfJn +a X (nn + 1 -1tn)) = ( ~ X (/Jn +a' 1tn + 1 -1tn)

= L\t(~, lfJn+a Xfn+a) + (~, lfJn+a X ln+a)r -L\t (grad[ ~ X CfJn +a], D(/Jn +aS) = L\t~ ' T~x! a -L\t ( !, Dcpn + aSDcp ~+a.) (3.6)
The second term in (3.5) can be further simplified with the aid of (3.3b) to obtain

(~,_ (q>n+l-q>n) X (ann +(1-a)nn+l)) = L\t(~, Jln+a X (21tn+(1/2) -1tn+a)) = 2L\t( ~' eo Vn +IX X Jln + (1 /2) ).
Inserting the last two expressions into ( 3.5) gives the result

(3.7) ~. (Jn+ I -Jn) = L\t~. T~~a + 2L\t(~, eo Jln+a. X Vn+(l /2) ), (3.8)
which collapses to (3.4) 2 if a =! and completes the proof.

Observe that exact conservation of the total linear momentum and the total angular momentum alone does not even guarantee consistency of the algorithm. As the preceding proof clearly shows, one can assign arbitrary values to the convected stress tensor S without altering the linear and the angular momentum within the step.

Algorithmic energy balance within a time step

The following results are independent of the specific form of the constitutive equation for the convected stress tensor S 1 • Within a typical time interval [tn, tn + d, the class of exact momentum-conserving algorithms obtained by setting a=! in (3.3a,b) satisfies the (finite) incremental relation

[ Kn + 1 -Kn] + < S, E( lf'n + 1) -E( lf'n)) A Kinetic Energy A Internal Stress Work L\ Potential Energy of Loading (3.9)
which can be viewed as the algorithmic counterpart of the expended power theorem (2.8) in the interval [tn, tn + 1 ]. By comparing the continuum relation

(2.8) with ( 3.9) we concluded that the change in the potential energy of the loading predicted by the algorithm is exact for dead loads. Likewise, the predicted change in kinetic energy is also exact. On the other hand, the change internal energy predicted by the algorithm depends on the form of the algorithmic approximation adopted for the constitutive equation which, so far, remains unspecified. To prove (3.9) we make use of the identity (3.10) which is easily verified by a direct computation, and use (2.5b) to derive the relation

sym[DqJ~+ ( 1 ; 2 )grad[~tVn + (1 /2)]] = sym[Dq>~+ (1/2)(Dq>n+ 1 -Dq>n)] =~[C(q>n+d -C(q>n)] =E(q>n+ 1) -E(q>n).
( 3.11)

Setting fJ = Vn + (I / 2 ) in the algorithmic weak form (3.3) and making use of these two results yields (3.9). Finally, let ~!0 ~~:. 1 n + 11 > 0 denote the algorithmic approximation to the time integral J;: + 1 !0~nt dt of the internal dissipation function. The algorith- mic approximation (3.3) is then said to be consistent with the reduced dissipation inequality if the following discrete counterpart of (2.9) holds: (3.12)

By analogy with the continuum problem, we view condition (3.12) as a restriction placed on the algorithmic approximation to the constitutive equation for S. In particular, in the absence of external loading [i.e., assuming condition (Al *)], if (3.12) holds then by combining (3.9) and (3.12) we obtain

£~0~ 1 -E~ot := [Kn + 1 + E~~.] -[Kn + £~nt] = -~!0 ~~nt• tn + 11 < 0, (3.13)
which implies the discrete conserving condition (3.1). For elastodynamics, the case considered below, ~!0 i~:. tn + d = 0 and (3.13) yields the condition of exact energy conservation. Consequently, (3.12) is the crucial condition to be verified in the algorithmic approximation to the constitutive equations.

An exact energy -momentum scheme for (smooth) elastodynamics

We show below that for elastodynamics the constraint to designate the right Cauchy-Green tensor and the Lagrangrian strain tensors associated with given configuration, respectively. To motivate the treatment of general constitutive models we consider first the simpler situation afforded by a quadratic potential.

Elastodynamics: The Saint Venant-Kirchhoff model

A key insight into the structure of the algorithmic constitutive relations is gained by considering first the case of a quadratic stored energy function of the form

e(C(tp)) = ~E(q>) • CE(tp), where C:=A1 ® 1 + 2!ll, (3.15)
which defines the so-called Saint Venant-Kirchhoff material. The wellknown limitations of this model include failure of the poly-convexity condition and incorrect behavior for extreme strains; see e.g., [START_REF] Ciarlet | Mathematical Elasticity Volume I. The Three Dimensional Theory[END_REF].

For the purpose of the present discussion, however, the model has the key property of leading to the stress-strain relation S 1 = CE(lf't) linear in the Lagrangian strain. As a result, using the identity

E(lf'n+ .) • CE(lf'n+ .) -E(q>n) • CE(lf'n) = [E(lf'n+.)-E(lf'n)] • C[E(lf'n+ .) + E(q>n)],
(3.16) the change in internal energy within a time step [tn, tn + d is expressed as

E~•~ 1-E~•••= t [e(C(tp.+ I))-e(C(tp.)) df!4 = ~ <[E(lf'n+ .) -E(lf',J], C[E(lf'n+.) + E(q>n)]). . (3.17)
In view of (3.12) from (3.17) we conclude that for a stored energy function quadratic in the (nonlinear) Lagrangian strains exact energy conservation for dead loading is achieved if

E~n~ 1 -E~nt = (S, E((/)n + 1) -E(q>n)) <=>IS= C[E(q>n + I) + E((/Jn)]/2.1 (3.18)
This result is unexpected. Since the internal energy term is evaluated at the mid-point configuration tl'n + (t/ 2 ) E Q it would be natural to adopt the constitutive equation S = CE(q>n + < 1 ; 2 ))• However, in sharp constrast with the result obtained in (3.18), this choice would result in failure of the property exact energy conservation in an algorithmic setting. In other words, for the Saint Venant-Kirchhoff model, use of the average of the momentum, by adopting the algorithmic constitutive equation S = 2Ve(Cn + Po) with /3 0 E (0, 1) such that: e(Cn+ de(Cn) = S •1(Cn+ 1 -Cn).

(3.23)

For general constitutive models it is shown below that the actual determination of {3 0 can be performed essentially at no additional cost as part of the standard iterative solution process for the (finite element form) of the momentum equation.

Remarks 3.1.

1. The existence of the required /3 0 E (0, 1) in (3.23) is guaranteed by the mean value theorem. In particular, if e( •) is quadratic we have /3 0 = ! is agreement with our analysis for the Saint Venant-Kirchhoff modcl.

• 2. Second order accurate algorithms. An accuracy analysis identical to that presented in Appendix I shows that the algorithm defined by (3.3a,b) along with (3.23) is second order accurate if and only if a = {3 0 = !. Thus, if the algorithmic constitutive equation (3.23) is adopted, exact energy conservation and second order accuracy can only be achieved for the Saint Venant-Kirchhoff model. This limitation is removed as follows. Define the function g: [0, 1]---+ lR by setting (see Figure 3.2) g(/3) ==~[e(Cn+P)-e(Cn + <t -m)J.

(3.24) Po E (0, ~] always exists such that g(l) -g(O) = g '(f3o) . . ',

• .

g(l)

Clearly, g( •) is odd about fJ = ~; i.e., it satisfies the property By the mean value theorem, there exists {3 0 E (0, 1) such that g( 1) -g(O) =

g(f3) = -g(l-{3) 'i/{3 E [0, 1
g'([J 0 ). Property (3.25) 1 then implies that the mean value theorem is also satisfied for 1 -{3 0 E (0, 1). Setting

(3.26) from (3.25) 1 and relation g( 1)-g(O) = g'({3 0 ) = g'( 1 -/3 0 ) we conclude that e(Cn+ J)e(Cn) = s. Hen+ I -Cn]. Therefore, the algorithmic constitutive equation also yields exact energy conservation [for a suitable {3 0 E (0, 1)] but, in contrast with relation (3.23), the resulting exact energy-momentum conserving algorithm retains the second order accuracy property even if {3 0 # ~. See Appendix I for a detailed accuracy analysis.

3. The energy-momentum conserving algorithm described above is clearly consistent. Unconditional (formal) nonlinear stability follows from exact energy conservation since, formally, the total energy defines a Lyapunov function for the dynamics.

D §4. Spatial discretization: Galerkin FEM projection

We examine below the Galerkin projection of the algorithmic weak form (3.3a,b) of the temporal discretization of the IBVP onto a finite dimensional phase and show that this projection preserves the momentum maps. Moreover, we show that the time discretization of the weak form commutes with the spatial Galerkin projection. Consequently, all the conservation properties derived in the preceding Section are inherited by the discrete (algebraic) finite element problem. For completeness, key aspects involved in the actual finite element implementation of the exact-momentum conserving algorithm for elastodynamics are briefly summarized.

The Galerkin FEM spatial projection

Let ? denote the (infinite dimensional) state space for the problem at hand, consisting of admissible configurations and material velocity fields on !16, defined as where (lfJh, Vh) E gh and (lfJ, V) E !7 is an arbitrary point in the state space.

The dimension of gh is nndof X nndof where nndor==ndim X nnode• Similarly, the finite element subs pace "// h c "// of admissible test functions is defined by

{ nnode } fh:= 'lh E fl'lh = L NA(X)cA, CA E (Rndim for A= 1, ... 'ndim • A= l ( 4.3)
In what follows we shall use the notation MAB== (NA, eoNB) = MBA to designate the entries of the "mass matrix" of the finite element discretization.

Projected algorithmic dynamics and momentum maps

By inserting the spatial interpolations defined by ( 4.2) and ( 4.3) into the time-discretization (3.3a) of the weak form of balance of momentum we obtain a nonlinear algebraic equation defined by the expression Furthermore, this diagram preserves the momentum maps in the sense that the spatial and temporal discretizations of the conservation laws of momentum, applied in any order, define the conservation laws of the discrete algebraic problem (for rx = 1).

1 fit [dn + lA -dnA] = Vn+ocA' ( 4.4) ...., nnode } ~ MAB (V V ] + Fint A Fext A i..J A n+ 1B-nB n+oc = n+oo B = l LJ.l ...
L ~ + 1 -L ~ = ~t L F~x! ~ for any a E [0, 1], A = 1

Algorithmic dissipation and energy balance within a time step

It can be easily shown that the incremental form (3.9) of the expended power theorem also holds in [/h. Again the result is independent of the order in which the Galerkin interpolation and the (conserving) time discretization are applied. Similarly, the algorithmic approximation restricted to [/ h is said to ~e consistent with the reduced dissipation inequality if the counterpart of (3.12) holds; i.e., c ~.@1~:.~n

+ 11 := <S\ E(q>~ + d -E(q>~))-[E~n~h1 -E~nt h] > 0. ( 4.10)
As in the time semi-discrete case, in the absence of external loading [i.e., assumption (A 1 *)] the total energy E 101 h == Kh + Eint h computed with an algorithm for which condition ( 4.10) holds satisfies the a-priori estimate

Etoth_Etoth:=-~.@inth <0 in!/hx[t t ] n + 1 n lm tn + 1] - n' n + 1 • ( 4.11)
For elastodynamics the total energy is conserved if ( 4.10) holds with L\.@1~:.~n + Il = 0. This is the case for the algorithmic constitutive equations (3.23) now evaluated on Qh. In summary, for elastodynamics the Galerkin finite element projection of the dynamics onto a finite dimensional phase space f/h preserves all the algorithmic conservation properties derived in the preceding section and commutes with the time discretization. The situation is not nearly as transparent in other models of continuum physics, such as nonlinear rods and shells, where straightforward Galerkin projection may not preserve the momentum maps and spoil, therefore, algorithmic conservation properties. We refer to [START_REF] Simo | On a stress resultant geometrically shell model . .e_art VI: Conserving algorithms for nonlinear dynamics[END_REF] for a detailed analysis of these and related issues in the context of nonlinear shells.

Implementation. Remarks on computational aspects

We outline below the steps within a typical time subinterval [tn, tn + d involved in the actual implementation of the exact energy-momentum conserving algorithm for nonlinear elastodynamics. Remarkably, this implementation is essentially identical to that of standard algorithms; The only additional effort required being the local solution (at quadrature points) of a scalar nonlinear equation. Clearly this addition involves negligible cost.

Assume that (dnA, VnA) are prescribed initial data at the nodal points and time tn, for A = 1, 2, ... , nnode. Since the discrete problem ( 4.4) along with constitutive equations (3.23) is nonlinear, an iterative solution strategy must be adopted. In what follows the approximate solution of ( 4.4) at the kth iteration within the time will be denoted simply by (dn+lA' vn+JA).

Step 1. Given (dn+ lA' vn+ lA), (A = 1, 2, ... 'nnode), compute the pro- jected deformation gradients D<p~, D<p~ + (t/ 2 ) and D<p~ + 1 by evaluation of the general formula

nnode Dtp~ + a = L [cxdn +!A + ( 1 -cx)dnA] ® grad[NA]
( 4.12)

A= 1

for a = 0, a = ~ and rx = 1, respectively. Then compute the right Cauchy-Green tensors: C~ = D<p~rD<p~ and C~+ 1 =

D hT D h lf'n+ 1 lf'n+ I•
Step 2. For fixed C~ + 1 (defined above) and given C~ define the (generally) nonlinear function g': (0, 1) ~ IR as g'(fJ) := e(C~ + 1) -e(C~) ( 4.13)

where C~ + a == rxC~+ 1 + ( 1 -rx)C~ is evaluated at a= f3 and cx = 1 -P in (0, 1), respectively. Then find the value {3 0 E (0, 1)

closest to ~ such that g'(/3 0 ) = 0. For the Saint Venant-Kirchhoff model /3 0 = t.

Step 3. Evaluate the convected stress tensor and the convected moduli by the algorithmic formulae Remarks 4.1

Sh== Ve(C~+Po) + Ve(C~+(J -Po>), } Ch==2[fJo V 2 [e(C~+Po) + (1-{J 0 )V 2 e(C~+O-Po>)].
1. The scalar equation g'(P) = 0 in Step 2 is easily solved since the solution Po E (0, 1) is precisely bounded. Thus any variant of Newton's method could be used, in conjunction with line search due to the lack of convexity of the stored energy e( •). Observe that the derivative of g"(P) is trivially computed by the closed-form formula g"(fl) = -1(C~+ 1-C~) • [V 2 e(C~+P)-V 2 e(C~+(1-fJ))](C~+ 1-C~).

( 4.17) 2. Strictly speaking, the matrix KAB obtained by exact linearization of the internal force vector F~n~1 112 ) is non-symmetric. In fact a standard calculation readily gives the expression x<kJAB =f.., [(grad[NA] • Sh grad[N"])l Geometric Term ( 4.18) Material Term which shows that the lack of symmetry is associated with the material term arising for the linearization of the algorithmic constitutive equation ( 4.14) 1 • We observe that this lack of symmetry is 'mild' and disappears as At~o.

• §5. Representative numerical simulations

The goal of this Section is to assess the performance exhibited by energy-momentum algorithm described above relative to two representative algorithms: The widely used trapezoidal rule and the symplectic midpointrule, shown in Simo, Tarnow and Wong [ 1992] to define an exact momentum preserving scheme. The trapezoidal rule (note the mid-point rule) can be obtained by setting {3 0 = ~ and y = ~ the classical Newmark family. We consider two representative simulations that exhibit finite deformations coupled with large overall rigid body motion.

Dynamics of a tumbling L-shaped block

In this first example we consider the dynamics of an L-shaped block subject to no boundary restrictions with material response governed by a Saint Venant-Kirchhoff elastic model. After an initial loading phase the block tumbles free of external forces in the ambient space undergoing finite deformation accompanied by large overall rotations and translations. All the data required to define the problem and perform the simulation are summarized in Figure 5.1.

According to our analysis in Section 2 the algorithmic constitutive

equation Sh = C~[E(cp~+ 1 ) + E(cp~)]
for the second Piola-Kirchhoff tensor along with the algorithmic weak form (3.3a,b) particularized at a=~ defines an exact energy-momentum conserving scheme. A main goal of this simulation is to verify numerically these properties. The sequence of deformed shapes obtained in the course of the simulation, shown in Figure 5.2 without any magnification of the actual deformations, clearly demonstrates the presence of finite deformations accompanied by a large overall rigid motion.

The sequence of plots summarized in Figure 5.3 contain the time histories total energy and the three components of the total angular momentum computed with the energy-momentum algorithm, the symplectic midpoint-rule and the trapezoidal rule of the Newmark family, respectively. The time history of the total linear momentum is not reported since these

Fx = 1600. • p(t) Fy = 800. • p(t) Fz = 1600. • p(t) Mx =0. My =0. Mz = -1200. • p(t) Fx = 800. • p(t) Fy = -1600. • p(t) Fz = -800. • p(t) Mx = -1200. • p( i) My =0.
Mz =0.

for t < 0.5 for 0.5 < t < 1.0 for t > 1.0

Geometry and finite element mesh for the simulation of the dynamics of a 'tumbling' L-shape elastic block with material response governed by the Saint Venant-Kirchhoff elastic model. The body is initially at rest and then subjected to the loading given in the figure. three algorithms, as well as most of the commonly use~ time stepping schemes, trivially preserve the conservation law of total linear momentum.

The results shown in Figure 5.3 include computations with several stepsizes and demonstrate the excellent performance exhibited by the energymomentum algorithm. In particular, the exact conservation of total energy yields an unconditionally stable scheme. In sharp contrast with this performance, both the symplectic mid-point rule and the trapezoidal rule exhibit a significant energy growth for the larger time-step sizes in long-term calculations, which leads to an eventual blow-up of the algorithmic solution and strongly suggests a loss of the unconditional stability property of these methods in the nonlinear regime. We remark that the exact conservation of total angular momentum exhibited by the mid-point rule is in agreement with the analysis in . Simo, Tarnow and Wong [ 1992]. Although the analysis in this latter reference shows that the trapezoidal rule fails to conserve total angular momentum, the results in Figure 3.3 indicate an excellent overall behavior. We attribute this good performance to the small rotational increments present in this problem as a result of the step-sizes employed. Nevertheless, failure of trapezoidal rule to conserve angular momentum may become significant in other situations.

Dynamics of a flexible satellite-like structure

In the second example we present the simulation of dynamics of a satellite~like structure with constitutive response again modeled by a Saint Venant-Kirchhoff material. pertaining to the geometry of the problem, initial conditions and material properties. As in the preceding example, after an initial loading phase the loads are removed and the subsequent motion of the structure in the ambient space exactly conserves momentum and energy. A representative deformed shape

I Fx = 9600. • p(t) load on = F _ 0 conic section Y - • Fz = 4800. • p(t)
{ 0.1 • t for t < 0.15 p(t) = 0.015-0.1 • t for 0.15 < t < 0.3 0. for t > 1.0 Geometry and finite element mesh for the simulation of the dynamics of a flexible satelite-like structure, with material response governed by the Saint Venant-Kirchhoff elastic model. The body is initially at rest and then subjected to the loading given in the figure. We have presented a new methodology for the construction of time stepping algorithms that inherit, by design, the conservation laws of momentum along with an a-priori estimate on the rate of decay of the total energy. For elastodynamics, the case examined in detail in this paper, we have presented a new class of second order accurate, unconditionally stable, exact energy and momentum conserving algorithms. The excellent performance exhibited by this conserving scheme has been demonstrated in two numerical simulations involving true finite strains accompanied by large overall rotations. It has been shown that the proposed method remains stable while two widely used algorithms, the Newmark trapezoidal rule and the symplectic mid-point rule, exhibit blow-up for finite step sizes. Remarkably, the actual implementation of these energy-momentum schemes is essentially identical to the one of standard schemes, and involves negligible additional cost.

The simulations in Section 5 demonstrate how algorithms possessing the same asymptotic order of accuracy and identical properties in the linear regime can exhibit drastically different responses in the nonlinear regime. In this connection, it should be emphasized that the three algorithms compared in Section 5 do become identical when restricted to the linear regime. Altho~gh these schemes are unconditionally stable in the linear regimes, their performance in the nonlinear regime is vastly different. The numerical result in Section 5 appear to suggest, therefore, that algorithmic conservation properties play an important role in the practical assessment of the actual performance to be expected from a given algorithm in the nonlinear regime. In particular, conservation of the total energy for a Hamiltonian system such as elastodynamics is not only desirable from a physical standpoint, but also a manifestation of unconditionally ( nonlinear) stability of the algorithm.

The ideas set forth herein are by no means restricted to elastodynamics and can be applied, in particular, to dynamical systems exhibiting (positive) internal dissipation. In future publications we plan to exploit this framework in the design of algorithms for dissipative systems such as the Navier-Stokes models. In addition, these ideas can be easily extended to Hamiltonian systems possessing a significantly more complex geometric structure, such as rod and shells. Finally, we plan to address the design of conserving schemes possessing higher order accuracy. unsubstantiated, that preserving the symplectic character of the Hamiltonian flow will enhance the long-term characteristics of an algorithm. In contrast with the notions of energy and momentum, which remain applicable to dissipative systems, the symplectic property is clearly restricted to Hamiltonian system. If the Hamiltonian system is non-integrable, a result of [START_REF] Zhong | Lie -Poisson Hami/ton -Jacobi theory and Lie -Poisson integrators[END_REF] implies that an energy preserving algorithm which is also symplectic can differ from the exact Hamiltonian flow only up to a time reparametrization and, therefore, it is unlikely to exist. Since nonlinear elastodynamics is a non-integrable system, the conclusions of the analysis given below are consistent with this observation. We show that the energymomentum conserving scheme defines a symplectic transformation only up to terms of order 11.t 3

• Additionally, we show that the mid-point rule is the only member of the one parameter family of generalized mid-point rule schemes which is symplectic. We remark that our numerical simulations in Section 5 strongly suggest that exact momentum and energy conser_vation result in significantly better performance in long-term dynamic calculations than exact conservation of the symplectic structure.

Discrete canonical Hamiltonian structure

We outline below the Hamiltonian structure of the finite dimensional dynamical system obtained via a Galerkin discretization on nonlinear elastodynamics, as described in Section 4. This Hamiltonian structure is considered in its own right, with independence of the Hamiltonian structure of the infinite dimensional dynamical system, as described in Simo, Marsden and Krishnaprassad [ 1988], which plays no role in these developments. Accordingly, given the discrete state space f/h described in Section 4, define pA:= nre MA 8 Vn with MA 8 := I eoNA(X)N 8 (X)d~ (II.4)

A= 1
Using if necessary a numerical quadrature formula with total number of quadrature points nquad weights g 1 , I= 1, ... , nquad the potential energy function becomes (11.5) where cp7 denotes the configuration cph evaluated at the quadrature point X 1 • Clearly, (11.5) defines a function of q E IR 3 nnode. Finally, let ~ denote the symplectic ndof x ndof defined by (11.6) where 0 and 1 are the zero and unit ndo f x ndof matrices, respectively. Recall that Jl possess the following well-known properties Jl = -~r = -Jl-1 and det[Jl] = 1.

(11.7)

With the foregoing notation in hand, the Galerkin discretization of nonlinear elastodynamics possess a canonical Hamiltonian structure with i. (Nodal) canonical phase space P = IR 2 ndof ii. (Separable) Hamiltonian H(z) = K(p) + U (q) iii. (Canonical) symplectic two-form w(bzh <5z 2 ) ==bz[ Jlbz 2 , for bzb bZ2 E TZP.

Here TzP denotes the tangent space at z E P of admissible variations. For such a Hamiltonian system, the canonical equations of motion take the standard form ! (z(t)) = JV H(z(t))} with V H(z) •= {~(~} .

(II.8)

z( t) It = o = Zo
We shall denote by !F 1 : P ~ P the flow generated by equations (11.8), defined by the standard relation z(t) = IF 1 (Z 0 ) .

l/.2. Symplectic algorithms

Recall that the symplectic two-form defined above is conserved by the Hamiltonian flow !F 1 ( •) In general, an arbitrary mapping G : P ~ P which satisfies condition (11.9) [or (11.10) with IF 1 replaced by G] is called a symplectic transformation. Symplectic maps are canonical transformations which preserve the form of Hamilton's equations, the volume in phase space and the integral invariants. In particular, an algorithm defines a time-one-map in phase space which advances the initial conditions Zn E P at time tn to the algorithmic solution Zn + 1 E P according to a generally implicit rule G(zn + 1 , Zn) = 0 consistent with Hamilton's equations. The algorithm is thus said to be symplectic if this time-one-map is a symplectic transformation. A direct verification of this condition involves considering the linearized algorithm bzn+ 1 = Abzn• The matrix A is called the linearized amplification matrix of the scheme and is given by

A =A! 1 A2
where A 1 == [oG(zn+ 1, Zn)J and A 2 == -[oG(zn + h Zn)J. It follows that the time-stepping scheme is symplectic if AT JJA = JJ. From a computational stand point, an alternative form of this condition that avoids inversion of the matrix A 1 , and hence much easier to check in applications, is derived as follows. From (11.11) we conclude that (11.12) and using properties (11.6) we arrive at the equivalent characterization I ATJJA = JJ.., A, JJA r-A2JJAI = o.l

We remark that the spectral stability properties of the algorithm can be inferred from the roots of the characteristic polynomial p(A.) := det[A -21 2 n].

Since det[ A] = 1 one concludes that the spectral radius (} of symplectic integrators must be (} = 1. As a result, a spectrally stable algorithm which exhibits numerical dissipation cannot be symplectic.

11.3. Approximate symplectic character of the energy-momentum scheme

First, we compute the amplification matrix of the one-parameter family of generalized mid-point rule algorithms and use condition (11.12) to show that the mid-point rule is the only symplectic member of this family. This result appears to be not well-known although the symplectic character of the mid-point rule is well-established, see Feng [START_REF] Kan | Difference schemes for Hamiltonian formalism and symplectic geometry[END_REF]. Subsequently, we exploit this result in an expansion of the discrete energy-momentum algorithm in powers of the time-step At to conclude that this latter scheme fails to be symplectic, at worst, to order (L\t) 3 . With these expressions in hand, condition (II.12) then becomes (11.15) which shows that, in general, the class of algorithm (11.13) is symplectic only if a = 1.

For the Galerkin finite element discretization of nonlinear elastodynamics in Section 4, V H(zn +IX) is computed from (11.8h evaluated at Zn + ex with aqAU(q) =f.., Dfl'!+(l/2)s::,pgrad[NA(X)] df-4, s~p := 2v e( cc cp~ + ( l /2))), (11.16)

1

  Invited Lecture, Presented at Oberwolfach January 1992. 2 Supported by AFOSR under Grant No. 2-DJA-826 with Stanford University.

  where the superposed dot denotes time differentiation and eo: 11 ~ !R.+ is the reference density. Throughout our developments we shall denote by S 1 the symmetric Piola-Kirchho.ff or convected stress tensor associated with the motion (/)r and related to the Cauchy stress tensor a 1 via the well-known relation (2.2) known as the Piola-transformation. In terms of the convected stress tensor the prescribed traction boundary condition takes the form

  we denote by < •, •) the L 2 -inner product on f1l of scalars, vectors or tensors, depending on the context. Setting (', )r:= Jra-(•)(•) dr, the weak form of the balance of momentum at time t E ~ takes the form (2.5a) Here grad[ •] is the gradient operator of vector fields on fJ.6 relative to . a (material) coordinate system {XA} (A= 1, ... , ndim) on BB; i.e., in /coordinates we have (grad[q])A = or,ajoXA where 'la (a= 1, ... , nctim) are the components of 'I E "//. DqJt, on the other hand, denotes the derivative of point maps relative to material coordinates; i.e. (DqJ 1 )~ = acpa;axA. In IR\ D[ •] and grad[ •] can be identified. The variational equation (2.5a) is supplemented by the local relations vt = lPt and nt = Qo vt in 1J X L along with the initial conditions (2.5b) (/)t lt=O = (/)o and 1tr lt=O = Qo Vo in~.

  Figure 3.1 Illustration of the three configurations involved in the formulation of the algorithm within the step [tn, tn + tJ. Note that all the convected stress and right Cauchy-Green tensors live in the reference configuration &1.

  balance of momentum 1 ~t <nn + 1 -1tn, '1)

  .4) In particular, since F~x! a = 0 and T~x! a = 0 if conditions (AI*) and (A2) hold, relations (3.4) yield the discrete momentum conservation law (3.1) within the interval [tn, ln+ d provided that a=-!. Relation (3.4) 1 follows immediately from (3.3a) merely by setting 11 =' (constant) in !!4 and noting that grad['] = 0. To prove relation (3.4h we take ' E IR 3 and use definition ( 3.2) to compute

E~n~ 1 -

 1 E~nt = (S, £( (/}n + 1) -£( (/}n)) -~!01~, tn + tJ ' with ~!0 ri~t r 11 > 0.n• n +

  ( 3.1) of exact energy conservation provides an explicit algorithmic constitutive equation for the convected stress S within a typical time interval [tn , tn + 1 ] that com- pletely defines the algorithmic approximation (3.3a,b). In what follows we use the notation C(q>) :=Dq>TDq> and E(q>) :=HC(q>) -1] for q> E Q, (3.14)

Figure 3. 2

 2 Figure 3.2 The mean value theorem applied to the odd function g(-) defined in Remark 3.1.2. Here

  ] :::> g(l)g(O) = e(Cn+ I)e(Cn).(3.25)

( 4. 1 )

 1 We denote by NA: 11 ~ ~ the prescribed nnode global finite element shape functions of a finite element discretization of the reference placement !16, with nodal points XA E dJ, A = 1, ... , nnode• As usual we assume the (completeness) condition NA(XB) = b~. The Galerkin projection of the state space !/ onto a finite dimensional state subspace space gh c !7 is defined via the formulae nnode cph(X) = L NA(X) dA and nnode Vh(X) = L NA(X)VA, vx E PJ, (4.2) A= l A= 1

for A = 1

 1 , 2, ... , nnode, where F~n~-! and F~x! ~ denote the internal and external forces associated with node X A at time tn +eo as defined by the relations. F:;'~~ •= L. Dli'~+•Sh grad[NA] ~, F~x!~ := r NAfn +oc dPA + r NAln+rx dr. J8l Jr u (4.5) Consider next the momentum conservation properties inherited by the algebraic problem ( 4.4). By inserting the interpolations defined by ( 4.2) into definitions ( 2.15) gives the following expressions for the projected total linear momentum and the projected total angular momentum: nnode nnode nnode L7 := L <eo, NA> vtA and J7 := L L

  ~~=~~ L dn+aA X F~x!~ only if (X=~• A= I In particular, in the absence of external loads [i.e., under assumption (Al *)] for the value r:t = ~ one recovers the discrete conservation laws of linear and angular momentum; namely, L~+ 1 = L~ and J~+ I = J~. The proof of the discrete relations ( 4.7) relies crucially on the following invarithe symmetry of Sh since ~ • J: dn +>A X F;ntA = (~:[(~X dn+•A) ®grad[NA]J, Dlp~+.sh)= ( e, Dcp~ +a Sh Dcp~: a) = 0, ( 4.9) for any vector ~ E IR 3 . With this result in hand, the proof of ( 4. 7) follows along the same lines as the argument presented in Section 3.2 for the time-discrete continuous problem and further details are omitted. Recall that Sh = ShT remains so far arbitrary.It can be easily verified that a result identical to that in ( 4.4) is obtained if the Galerkin projection is applied first to the continuous weak form (2.5a,b) of the momentum equation, followed by the time discretization. Therefore, the following diagram commutes:

(

  F~n~1 1 12 > and F~x!t-1 112 ) via formulae ( 4.5) and evalu- ate the dynamic residual as If IlRA/I < TOL then (dn+IA' vn+lA) is the sought solution.Step 4. If )jRAJJ > TOL then solve for 11dn+JA the linearized system n i e [A 2 2 MAB1 + KABJ /1dn + lB +RA= 0.

  ( 4.16)B = I LllUpdate the solution as dn +lA +--dn + lA + fld 11 + lA and return toStep 1.

  Figure 5.1

Figure 5. 2

 2 Figure 5.2 Sequence of deformed shapes of the L-shaped block obtained with the exact energy-momentum conserving algorithm for nonlinear elasto-dynamics.

Figure

  Figure 5.3 Time histories of the total angular momentum and the total energy computed with the exact energymomentum conserving algorithm, the symplectic mid-point rule and the trapezoidal rule.

Fx

  Figure 5.4

Figure 5. 5

 5 Figure 5.5 Finite element mesh and representative deformed shape obtained in the course of the full dynamic response of a flexible satellite-like structure, with material response governed by the Saint Venant-Kirchhotf elastic model.

  Figure5.6 Time histories of the total angular momentum and the total energy computed with the exact energymomentum conserving algorithm, the symplectic mid-point rule and the trapezoidal rule.

(

  II.3) Next, one defines the potential energy function U: ~ndor --t !R as follows. Recall that the deformation gradient Dtph and the right Cauchy-Green nnode Dcph = L dA ® grad[NA(X)] and C(cph) == [Dcph]TDcph.

  in the sense that the relation (II.9) holds for arbitrary variations Jzb Jz 2 E P. Here DIF 1 (Z) denotes the derivative of IF 1 at z E P defined in canonical coordinates by the 2ndof x 2ndor matrix of partial derivatives. Condition (11.9) is equivalent to the requirement I [DIF,(z)VJIDIF,(z) = Jl; I (II.IO) in other words, the matrix D [F 1 (z) is symplectic. Condition (11.9) has a number of well-known consequences. In particular, from (II.IO) one concludes that the Hamiltonian flow preserves the volume in phase space (Liouville's theorem) since det[D IF 1 (z)] = 1.

  1/.3.1. The symplectic midpoint ruleFor a general canonical Hamiltonian system, the one-parameter family of generalized mid-point rule algorithms takes the formG(zn + b Zn) := Zn + 1 -Zn -AtJJV H(zn +IX) = 0,} Zn+IX = azn+ 1 + (1-!X)Zn for IX E [0, 1]. (11.13)Differentiating expression (II.13) 1 and using the notation in(11.11) gives the resultA 1 =1-aL\tJJV 2 H(zn+ct) } A 2 = 1 + (1-a)L\tJJV 2 H(zn+ct)(II.14)

  Jn+(l /2) ==2 L..A = IN (X nA + n+ lA lS t e ffil .. point con guratlon and C(qJ~ + (t / 2 ) ) is the associated right Cauchy-Green tensor defined by (II.4h.

  M AB dtA X vtB' V t E L Now let L~, L~+ 1 and J~, J~+ 1 denote the total linear and total angular momentum at times tn and tn + 1 , respectively, obtained by evaluating expressions ( 4.6) for the nodal values (dn, Vn) and (dn+ b Vn+ I). If (dn +I, Vn + I) E !Rnnode X !Rnnode satisfies the algebraic problem ( 4.4) then for the Neumann problem [i.e., under assumption (A2)] the following relations holds nnode

	A= 1	A=lB=1
		(4.6)

  X nnode• It follows that z E p == !R 271 dof. We shall refer to p as the nodal phase space associated with Galerkin finite element discretizations of nonlinear elastodynamics. On P one defines the kinetic energy function via the expression

	B = l	J~		(ILl)
	and introduce the following vector notation no e an z == q == (d'l' d2, ... 'dn d) E [Rndof, d p := (p',p2, ... ,pnnode) E ~ndof,	{q} p •	(II.2)
	where ndor=== 3

K(p) == ~p . M -1 p > 0, with M= [MA 8 13 X 3] E ~ndof X ~ndof.

Acknowledgements

We are indebted to D. Lewis and K. Wong for helpful comments and many discussions on the subject of this paper. Support for this research was provided by AFOSR under contract No. 2-DJA-826 with Stanford University. This support and the interest and comments of Dr. Marc Jacobs are gratefully acknowledged.

Appendix I. Accuracy analysis of energy-momentum schemes

In this Appendix we outline the accuracy analysis of the discrete energy-momentum schemes presented in Section 4. Recall that the time discretization of the algorithmic weak form is ( Qo(Vn + 1 -Vn), tl ) + ~t (DlfJn+ ( 1 /2)S, grad [t~] ) :=rAt -0, tit lfJn+ I -lfJn-l [Vn+ 1 + Vn] ==!M-0, (1.1) where we have assumed for simplicity zero external loading. We restrict our attention to the case in which the algorithmic constitutive equation for S is defined as (1.2) with f3 E (0, 1) defined so that energy is exactly conserved; as explained in Section 4. The analysis of the first order accurate scheme, corresponding to S = 2Ve(Cn + 13 ) follows along identical lines and further details are omitted.

Let ( lfJ( •, tn + 1 ), V( •, tn + 1 )) denote the exact solution at time tn + 1 defined on f?J. A Taylor series expansion about tn then yields

(1.3) Setting Fn + 1 := Dl{J( •, tn + 1 ) and Fn : = DqJ( •, tn) a Taylor series expansion of en + 13 yields

(1.4)

Inserting this last result into the Taylor series expansion of the algorithmic constitutive equation ( 1.2) about tn gives the following result

= 2Ve(Cn) + AtV 2 e(Cn): Cn + lP(At 2 ).

(1.5)

Now let L(At) denote the result obtained by inserting the exact solution into the left-hand-side of the algorithmic equation (1.1 ) 1 • Making use of the expansions (1.5) and ( 1.3) associated with the exact solution we obtain the following result for (1.1)h ordered in increasing powers of At:

Similarly, denoting by -r(At) the result obtained by inserting the exact solution into the left-hand side of the algorithmic equation ( 1.1 ) 2 we find

By assumption, the weak form of the equations and its time differentiation hold at time tn. Consequently, in view of the form of (1.6a,b) it follows that (1.7)

Recalling the classical definition of the local truncation error (see e.g., Richtmyer and Morton [ 1967]) we conclude that the proposed energymomentum algorithm is second order accurate for any value {1 E [0, 1].

Appendix II. Symplectic integration of Hamiltonian systems

This Appendix summarizes the Hamiltonian structure underlying Galerkin finite element discretizations on nonlinear elastodynamics. The main goal is to investigate to what extent the proposed discrete energymomentum conserving algorithm defines a symplectic scheme. Recall that an algorithm is called symplectic if the (generally implicit) map that advances the solution in phase space from time tn to time ln + 1 defines a symplectic transformation. This type of algorithms has received a considerable amount of attention in recent years; see e.g., [START_REF] Ruth | A canonical integration technique[END_REF], [START_REF] Channell | Symplectic Integration Algorithms[END_REF]~ [START_REF] Lasagni | Canonical Runge-Kutta methods[END_REF], Sanz-Serna [ 1988] and the recent review paper of [START_REF] Scovel | Symplectic numerical integration of Hamiltonian Systems[END_REF]. The underlying motivation is the belief, to some extent

The discrete energy-momentum method

The algorithmic equations that define the exact energy-momentum conserving algorithm now take the following form ~~'(

Thus, in contrast with the symplectic mid-point rule which is defined by (11.13) with rx =!along with (11.15), the discrete energy-momentum method uses F~n~(l/ 2 ) in place of VU(q). We recall the expression (11.20)

As a result, we obtained the following estimate for condition (11.12):

AimJI[Aim]r-A2mJl[A2m]r = [(9(L\t3)], (11.21) which shows that the discrete energy-momentum method defines a symplectic transformation up to terms of order (9(L\t 3 ).