Numerical analysis of the counterintuitive dynamic behavior of the elastic-plastic fully-clamped beams under impulsive loading

Mehdi Shams Alizadeh*, Kourosh Heidari Shirazi**, Shapour Moradi***, Hamid Mohammad Sedighi**** *Department of Mechanical Engineering, Shahid Chamran University of Ahvaz, Ahvaz, Iran, P.O.B. 135, E-mail: mshalizadeh@gmail.com **Department of Mechanical Engineering, Shahid Chamran University of Ahvaz, Ahvaz, Iran, P.O.B. 135, E-mail: k.shirazi@scu.ac.ir ***Department of Mechanical Engineering, Shahid Chamran University of Ahvaz, Ahvaz, Iran, P.O.B. 135, E-mail: moradis@scu.ac.ir; ****Department of Mechanical Engineering, Shahid Chamran University of Ahvaz, Ahvaz, Iran, P.O.B. 135, E-mail: h.msedighi@scu.ac.ir


Introduction
In 1985, an interesting fact was discovered in the USA by a well-known specialist on the dynamics of elastoplastic beams, Symonds and his colleague Yu [1].Computations made by the two scholars showed that under certain conditions when a transverse pressure pulse of magnitude sufficient to cause plastic deformation is applied to a elastic-plastic beam whose ends are attached to smooth fixed pins, the permanent deflection of the beam may come to rest in opposite direction of the load rather than in its direction, as intuitively expected and observed in many tests.This phenomenon was named "anomalous behavior" or "counter-intuitive behavior" of the beam.For the particular problem defined by Symonds & Yu [1], solutions were obtained by ten numerical codes of finite element or finite difference type.Some predicted a final elastic vibration between negative displacement limits (implying a final displacement in the direction opposite to the load), while others predicted positive limits.Although later further investigations are also carried out on counterintuitive behavior of pin-ended beams [2][3][4], as it is shown by other studies [5][6][7], the pin-ended condition is not essential and similar behaviors of this sort may be expected to occur in calculations for elastic-plastic beams with fully clamped ends under impulsive loading.Occurrence of this anomalous response requires the beam to have boundary conditions which prevent the axial displacement so that the transverse displacement causes stretching of the middle surface and the plastic strains in effect convert the beam to an arch.The anomalous behavior dealt with here is a consequence of dynamic instabilities akin to snap-buckling.
Using direct differentiation approach and central difference method, Li and Liu [8] applied the method of parameter sensitivity to study this extremely sensitive problem.They [9] also investigated the counter-intuitive response of a beam from the viewpoint of parameter uncertainty.In addition, Li et al. [10] used a 3DOF elasticplastic Shanley-type model to understand the dependence of anomalous region on system parameters.
The counterintuitive behavior of some other structures has also been the focus of attention, by analytical and finite-element simulation, Li et al. [11] investigated the counter-intuitive behavior of elastic-plastic aluminum ring subjected to radial pulses of pressure loading.Flores-Johnson and Li [12] presented a numerical investigation of the Counter-intuitive response in elastic, perfectly plastic square plates subjected to impulsive loading.Dong et al. [13] extended the energy method for plane-stress rings to plane-strain rings and spherical shells for the determination of the counter-intuitive regions.Having employed the law of thermodynamics and the theorem of Lyapunov instability, Zhao et al. [14] investigated the counter-intuitive phenomena of elastic, perfectly circular and square plates.Ma et al. [15] employed the commercial finite element package, LS-DYNA, to simulate the counter-intuitive response of a single-layer reticulated dome which was initially stressed by static preloading to an interior blast.
There is no report on the proof of the Symonds' finding [1] by experiment under the same conditions, whereas strong evidence is provided in ref. [16] for the existence of the counterintuitive phenomenon in an elasticplastic beam under dynamic load.Although the anomalous responses of both share the same characteristics (i.e.their final displacements just opposite in direction to the impact velocity), they have different boundary and loading conditions.The former is a pin-ended beam subjected to a rectangular pressure pulse, while the latter is fully clamped beams subjected to projectile impact at the mid-span.The results of ref. [16] showed that the occurrence probability of the counterintuitive phenomenon is about fifty percent among the total tests and it is, therefore, evident that the observed phenomenon is uncertain and the existing deterministic models cannot predict it.
Kolsky, et al. [17], through some tests on aluminum alloy beam specimens, made another attempt to observe anomalous behavior previously predicted by Symonds and Yu [1].Since satisfactory tests with pinned ends are more difficult to set up, they used beam specimens with fully fixed ends; a simple technique was applied, in which the impact loading was simulated by pulling the specimen to an initial deflection and then releasing it abruptly.The results of this research showed that only two out of 36 such tests led to permanent deflections in the direction opposite to that of the imposed initial deflection.
The anomalous response of the beams found in experiment of ref. [16] is numerically modeled for pinended beams by Li and Qin [18] using ANSYS LS-DYNA.
In the present paper, the behavior of the fully-clamped elastic-plastic beams found in Kolsky, et al. experiment [17] was studied using finite element code ANSYS LS-DYNA and Galerkin's method.We studied in detail the displacement-time history curves of mid-span of the beams and determined the region of the occurrence of the counterintuitive behavior.Furthermore, using finite element code, the energy diagrams of the beams were also investigated.

Description of the problem
The beam found in Kolsky, et al. experiment [17] was applied in the present study; in this experiment, in order to create the initial deflection, a piano wire is attached to the specimen at its midpoint by a simple mechanism.The beam is bent by attaching the other end of the wire to the testing machine and pulling it.When the desired initial deflection is reached, the wire is cut close to the beam and the subsequent beam's response is recorded.Here, too, due to the quasi static initial forces, the midspan deflection increases up to the define peak value, Then, the unloading takes place and the beam carries out nonlinear vibrations and the beam's behavior after unloading is studied.Here, too, similar to what [17] assumes that the beam is fully clamped at two ends and has the span L = 203.2mm and the uniform rectangular cross-section with the width b = 12.7 mm and the height h = 2.54 mm, as shown in Fig. 1, a.
The beam's material is aluminum alloy 6061-T6 which, with reference to stress-strain diagram shown in ref. [17], has approximating linear elastic-perfectly plastic behavior which implies equal yield stresses in tension and compression as well as absence of strain hardening and strain rate effects (as shown in Fig. 1 ANSYS LS-DYNA, an explicit finite element code for analyzing the large deformation and the dynamic response of inelastic solids, was employed in the present study.The beam was modeled by 20-beam element of BEAM161.The boundary conditions of the modeled beam were chosen so that the first and the last nodes have fullycamped conditions; therefore, all degrees of freedom consisting of the displacements UX, UY, UZ together with the rotations ROTX, ROTY, ROTZ were chosen as zero at these support nodes.For the initial conditions, Quasi-static initial forces were applied at the middle node, causing the mid-span deflection to increase up to the definite peak value, post to the occurrence unloading and before the beam's behavior after unloading being investigated.The behavior of the material was assumed as linear elasticperfectly plastic.The selection of either kinematical hardening or isotropic hardening has no influence on the simulation because the tangent modulus of the material equals zero.Details of these inputs as well as loading and boundary conditions are provided in ANSYS LS-DYNA User's Guide [20].

Galerkin's method
The present paper applied Galerkin's Method as an alternative approach to study displacement-time history curves of mid-span of the beams and to determine the region of the occurrence of the counterintuitive behavior.This method is usually regarded as one of the methods of weighted residuals.
The dimensionless equations of motion for the elastic-plastic beam shown in Fig. 1 Where the dimensionless quantities indicated by hats are defined and related to dimensional quantities with the following relations: ; ; ; Where x and z are the axes of coordinates according to Fig. 1; t is time; u is axial displacement; w is deflection; T is axial force, and M is bending moment.
Eq. ( 1) are integrated by the method of Galerkin.Due to symmetry, we shall consider only one half of the beam for which Henceforth, primes and dots denote differentiation with respect to x and ˆt .
Integrating by parts and taking into account the boundary conditions The displacements û and ŵ are assumed in the following form, in which all boundary and symmetry conditions are satisfied: By substituting these expressions into Eqs.( 8) and ( 9), and carrying out the integrations, the following system of equations were obtained: and To calculate the axial forces and the bending moments in the integrals (15), it can proceed with the following manner: assuming that the hypotheses of Kirchhoff hold, the axial deformation can be found from the below formula: In the strain-stress diagram for cyclic loading of the assumed linear elastic-perfectly plastic material shown in Fig. 1, b, Segment (I) is correspond to pure elastic deformation; segments (II) and (IV) are respectively correspond to plastic tensile loading and plastic compressive loading.In the case of segment (III), elastic unloading takes place.
For evaluation of the stresses, we need to know exactly on which segment we find ourselves at the present instant.For this purpose, it is required to record the values of  and ê from which the unloading begins.These are indicated by symbols

ˆˆT x t x z t dz
These integrals as well as integrals (15) shall be evaluated numerically.Next, we calculate from Eqs. (12)(13)(14) In order to start the programming of this method, it is required to indicate the initial values; it is assumed that the initial dimensionless displacement of the beam can be defined using the following relation: where ini ini ŵ w / h  is the initial dimensionless displacement field of the beam and max max ŵ w / h  is the initial dimensionless maximum displacement of mid-span of the beam.
By connecting this relation with the displacement field found from Eq. ( 11), it can be concluded that Later, we shall take: For numerical analysis, by calculating the integrals ( 19), ( 20) and ( 15 and the algorithm of solution can be carried on.

Displacement-time curves for beams under impulsive loading at mid-span
Due to the quasi static initial forces, the mid-span deflection increased up to the defined peak value.Here, special attention was paid to the beam's response after peak deflection.After the peak deflection, unloading occurred and the beam responded with nonlinear vibrations.In order to analyze this phase of the motion, it can comment the mid-span displacement-time diagrams, some of which being presented in Figs.2-8.
In the case of small initial displacement, as shown in Fig. 2, the beam carried out regular vibrations with constant amplitude around the equilibrium position w = 0, showing that the beam's material remains elastic.By increasing the initial displacement, the beam behaves with irregular vibration with decreasing amplitude after one cycle and the constant amplitude thereafter which is, due to the appearance of the plastic deformation and energy dissipation, in the first cycle of the elastic-plastic motion (Fig. 3).By further increasing of the initial displacement, first, as shown in Fig. 4, more irregularities appeared in the vibrations of the motion.Then, these irregularities in compound with the dynamic instabilities occurred in the beam's behavior, surprisingly tending to move the beam to negative region, which is a sign of the occurrence of the counterintuitive behavior in the beam's motion (Fig. 5).Next, as shown in Fig. 6, the beam's behavior became quasi periodic with smaller amplitude in the counterintuitive region.Then, a similar dynamically instable mechanism, which moved the beam's vibration transiently into the counterintuitive region, caused the direction of the vibration to tend backwards to positive region permanently (Fig. 7).Thereafter, the beam vibrates in positive region with quasi periodic small amplitude stable vibrations and the further increase in the initial displacement would not affect the vibrations behavior of the beam and only causes more decrease in stable vibrations amplitude due to further energy dissipation as a results of the increase in the initial plastic strain (Fig. 8).As shown in Figs.2-8, the derived solutions of both the finite element method and the Galerkin's method are capable to predict the counterintuitive behavior of the beam.These solutions coincide with one another in the positive region (i.e. in direction of the loading) and in the counterintuitive region (i.e. in opposite direction of the loading), with high degree of accuracy, while there is some divergence in the transient zone.It was observed; based on author's experience; that the computational cost of the Galerkin's method is less than that of the finite element solution, which is mainly due to the simpler formulation of the Galerkin's method compared to the finite element method.
Also, it should be noted here that Figs. 6 and 7 are respectively related to the initial displacements due to the quasi static concentrated forces 1510N and 1670N at midspan which are corresponding to the initial forces related to Figs. 2 and 3 (Table 1) in Kolsky, et al. experiments [17].As it can be observed here, the numerical simulations similar to the experimental observation predict the counterintuitive behavior for these two cases of the loading; the only difference is the slight difference between the predicted initial deflections by these two forces, which is due to the differences between the numerical and experimental simulations some of which being mentioned in the following section.It can be observed in Fig. 9 that according to the finite element method, the region of counterintuitive behavior exists for the values of initial displacement which are approximately between 8.6 mm and 11.3 mm.These values of initial displacements are due to initial quasi static forces which are between 1225 N and 1800 N.

Determination of the counterintuitive region
Also, for the Galerkin's method, as it is shown in Fig. 10, for the initial displacements which are approximately between 8.7 mm and 11.5 mm, the counterintuitive behavior will occur.These values of initial displacements are due to initial quasi static forces which are between 1230 N and 1850 N, therefore the predicted counterintuitive region by Galerkin's method has a slightly wider band in comparison to the one with finite element method.It is noteworthy that although these numerical simulations show a continuous region within initial loadings of which a counterintuitive behavior results, still, as it mentioned above, the results of Kolsky, et al. experiment [17] shows that among several tests, only a few of them led to permanent deflections in the direction opposite to that of the imposed initial deflection.Differences which exist between the predicted behavior in numerical analysis and the experimental observations may be as a result of the following: 1.The complete end fixity assumed in the finite element code and Galerkin's method cannot be exactly realized in the laboratory.
2. Differences in the shapes and material properties of the test specimens in Kolsky, et al. experiment [9] might be noted; Furthermore, the specimens might have residual stresses that can be responsible for some of the scatter observed in the tests.
3. Reaching the sudden unloading as in the numerical simulations may be difficult in the test process.Results of the complementary numerical analysis show that the increase in unloading time duration from zero will decrease the width of the region of occurrence of the counterintuitive behavior.

Energy diagrams
Energy Diagrams is an approach introduced in [22,23] to display and interpret the response features simpler and more direct than common studies in terms of phase space geometry.In the present class of problems where small plastic deformations play a major role, phase plane diagrams are not very helpful because they indicate these effects only indirectly but the energy plots show directly the occurrence of plastic strain increments, and their effect on the evolving response.This involves a simple idea namely that of plotting together the total energy, internal energy, and kinetic energy history curves.The discontinuous alternation of the final state between positive and negative values, and all other features of the characteristic diagram of this model become immediately understandable from energy diagrams [22].
Using the finite element code, the total energy, internal energy, and kinetic energy history curves for some previous dynamic responses, such as the responses related to Figs. 3, 5 and 8 are respectively shown in Figs.11-13.Using the conservation law of energy and ignoring the heat loss, as it can be seen in the mentioned figures, the total energy consists mainly of the internal energy and kinetic energy; therefore, the increase in the internal energy is equal to the decrease in the kinetic energy.Since in these diagrams, the components of energy all vary with time and their transformations from one form to another are very complicated, diagrams of the maximal kinetic energy after peak deflection and the corresponding internal and total energy with respect to the mid-span peak deflection are presented in Fig. 14.Regarding these diagrams, three distinct regions can be recognized: In the first part of the diagrams which is correspond to small initial deflection, the internal energy value comes close to the kinetic energy and the beam itself mainly vibrates elastically or with small plastic strain around the equilibrium position.In the last region of the diagrams which shows the large initial deflection, the internal energy is much larger than the kinetic energy and in this situation, no sufficient kinetic energy is available to make the beam move to and stay in the opposite direction of the initial deflection.Therefore, in this case, the counter-intuitive behaviors cannot occur.But in the middle part of the diagrams, there are both proper plastic deformation and enough kinetic energy to make the beam move to and finally stay in the opposite direction of the initial deflection, and, in other words, the counterintuitive behaviors will occur.Therefore, the occurrence of the counter-intuitive behaviors is a result of the adequate plastic deformation as well as the proper proportion of the internal and kinetic energies.
It is depicted in Fig. 14 that, for the beams under investigation, the counterintuitive behavior will possibly occur for all the motions with sufficient plastic deformation and with internal energy within the region that is greater than kinetic energy and lower than three times of kinetic energy.

Conclusion
In this study, the counterintuitive dynamic behavior of the elastic-plastic fully-clamped beams formerly found in Kolsky, et al. experiment [17] is numerically simulated.Finite element code ANSYS LS-DYNA and Galerkin's method are applied to numerical modelling of the problem.Special attention is dedicated to the nonlinear vibrations which proceed after the peak value of deflections.Displacement-time history curves of mid-span of the beams were studied in detail, showing that more irregularities in the vibrations of the motion appear around the counterintuitive behavior, and these irregularities in compound with the dynamic instabilities occurred in the beam's behavior, moving the beam's motion to negative region.The results of these investigations also indicated that both numerical methods are capable to predict the counterintuitive behavior of the beam.The predicted responses of the Galerkin's method as well as those of the finite element method coincide with each other both in the positive region (i.e. in direction of the loading) and in the counterintuitive region (i.e. in opposite direction of the loading), with high degree of accuracy; nevertheless, there are some differences in the transient zone which delineate that the counterintuitive responses are dependent on the numerical method of solution.The investigations also showed that the computational cost of the Galerkin's method is less that of the finite element solution.
One of the main goals of this paper was to determine numerically the region for the beams found in Kolsky, et al. experiment [17] within the initial loadings of which a counterintuitive behavior results.Although the experimental observation [17] showed that among several tests only a few of them led to permanent deflections in the direction opposite to that of the imposed initial deflection, still the computations performed by the two abovemotioned numerical methods showed that there is a continuous region for the occurrence of the counterintuitive behavior.The differences which exist between the predicted behavior in numerical analysis and the experimental observations may have different reasons such as the difference between boundary conditions.In addition, the analysis of the energy diagrams showed that for the beams under investigation, the counterintuitive behavior will possibly occur for all the motions with sufficient plastic deformation and with internal energy within the region that is greater than kinetic energy and lower than three times of kinetic energy.
Fig. 1 a -geometry of beam; b -linear elastic-perfectly plastic stress-strain behavior

 1 
the instant t = 0, it takes 0 mm ˆê   .If plastic deformation develops in the beam, the values  we also have elastic deformations.Eq. (17) holds with the only exception that if the value of  calculated from (17) is greater than 1 is anological to case (2).If Eq. (17) gives a value  we have to sub-After calculating the stress field for any point of the beam at the given instant, we shall find the quantities T and M from the below equations: cording to the central difference method: ) we shall find from Eqs. (12we can evaluate from Eqs. (19) and (20) the quantities 12 ,,

Fig. 6 Fig. 7 Fig. 8
Fig. 6 Mid-span displacement time history curve of the beam, 1510 N initial quasi-static force at mid-point (wmax = 10 mm) It is of interest to see how the behavior of the beam depends on the initial displacement.Results of related calculations by the two used methods are presented in Figs. 9 and 10, where curve (a) presents the peak value of mid-span deflection, curves (b) and (c) are respectively the maximal and minimal values of vibrations subsequent to the peak value.To this diagram, the following comments can be made: In the case of small initial displacements the beam material remains elastic and the beam carries out elastic vibrations around the equilibrium position w = 0; the curves (a) and (b) coincide.By increasing initial displacement, the motion is elastic-plastic, but the curve (b) remains in the positive side of the diagram.After that comes a region of initial displacement, where the curve (b) has negative values which show the counter-intuitive behavior of the beam.If it shall increase initial displacement, still more the beam's response changes vary abruptly and the deflections are positive all the time.

Fig. 11
Fig. 11 Energy time history curves of the beam, 800 N initial quasi-static force at mid-point (wmax = 6 mm)