Flexure and shear response of an impulsively loaded rigid-plastic beam by enhanced linear complementarity approach

The linear complementarity approach has been utilized as a systematic and unified numerical process for determining the response of a rigid-plastic structure subjected to impulsive loading. However, the popular Lemke Algorithm for solving linear complementarity problems (LCP) encounters numerical instability issues whilst tracing the response of structures under extreme dynamic loading. This paper presents an efficient LCP approach with an enhanced initiation subroutine for resolving the numerical difficulties of the solver. The numerical response of the impulsively loaded structures is affected by the initial velocity profile, which if not found correctly can undermine the overall response. In the current study, the initial velocity profile is determined by a Linear Programming (LP) subroutine minimizing the energy function. An example of a uniform impulsively loaded simply supported beam is adduced to show the validity and accuracy of the proposed approach. The beam is approximated with bending hinges having infinite resistance to shear. Comparison of the numerical results to the available closed-form solution confirms the excellent performance of the approach. However, a subsequent investigation into a beam having the same support conditions and the applied loading, but with bending and shear deformation, results in numerical instability despite optimizing the initial velocity profile. Thus a more generic description of kinetics and kinematics is proposed that can further enhance the numerical efficiency of the LCP formulation. The ensuing numerical results are compared with the available close form solution to assess the accuracy and efficiency of the developed approach.

www.nature.com/scientificreports/ loaded rigid plastic structures is presented by Kaliszky and Logo 26 . A further extension of the quadratic program has been rigid plastic modeling of walls subjected to extreme dynamic loading 27 . Essentially, the concept of mathematical programming has been successfully employed by Smith and Sahlit 28,29 and Khan et al. 30 . There seems to be considerable interest in the study of rigid-plastic structures under extreme dynamics loading as is reflected by the large body of work in the area; Lee and Symonds 31 , Symonds 32 , and Jones 33,34 . The idea of utilizing rigid-plastic idealization was suggested by Lee and Symonds 31 while investigating the response of free-free beam subject to blast loading. Later, Parkes 27 devised the concept of the traveling plastic hinge, thus motivating series of inspired works 35,36 . More recently, this simple theory has been extended to incorporate the effects of the bending-shear interactions, strain rate and large displacements [37][38][39][40][41][42][43] .
The effect of shear forces appears to be important in evaluating the quasit-static behavior of structures [43][44][45][46] , but more significant in the structures submitted to extreme dynamic loading, Bleich and Shaw 47 . The shear force effect can markedly influence the behavior of those beams in which the dynamic response involves higher modes, Symonds 48,49 and Jones 50 . It transpires that the shear effects are more important for the I beams than for the rectangular beams of the same length; Menkes and Opat 51 , Jones 52 and Karunes and Onat 53 . Regardless of the properties of any cross-section, the interaction curve relating the shear force and the bending moment cannot be uniquely defined. Mathematical convenience has motivated researchers to approximate the interaction relation to rectangular yield criterion, Symonds 48,49,54 . This approach has lead to the notion of shear sliding within the sheared zone is commonly adopted 50,52,[55][56][57][58][59] , which can be conveniently modeled by 1D rigid-plastic solid elements 30 .
A detailed 1D rigid-plastic modeling strategy has been developed previously 25,30 which was validated against the continuum solution of the impacted beam. This numerical approach, formulated as a Linear Complementarity Problem (LCP), led to fast and accurate dynamic response prediction, unlike the commercial finite element softwares that require extensive computer memory and a longer execution time 60 . However, there were instances when the accuracy of the LCP solver 25 was undermined by issues concerning numerical instability. These problems were attributed to the approximate initial velocity profile, which is required to initiate the numerical process.
To achieve the numerical stability of the LCP approach, the current study proposes a new initiation technique that utilizes a linear program (LP) to determine the actual initial velocity profile. Lumped mass and continuous mass elements are employed to model simply supported beam having infinite resistance to shear deformation. The beam is subjected to impulsive loading and the initial velocity profile is determined from the newly developed LP, which is implemented in the LCP approach 30 . The accuracy and numerical efficiency of the results are confirmed by comparing the results with the continuum solution of the beam problem. However, the numerical difficulties reappear in the subsequent investigation of the impulsively loaded beam in which shear deformation is allowed at the support location. An improved kinetic and kinematic description is developed to circumvent these difficulties. Results from the shear critical simply supported beam ensure the computational robustness of the LCP using this general description.

Improved dynamic rigid-plastic model
In the previous research 30 , the formulation in the form of a Linear Complementarity Problem (LCP) has been developed to determine the response of rigid-plastic beams under impact loading. The initiation sub-routine of that formulation is modified herein to resolve the numerical difficulties that appeared previously. Besides, the formulation is extended to beams under impulsive loading, which shows pronounced sensitivity with the current semi-definite LCP 28,29 . Representation of kinetics and kinematics as nodal description. The kinetic and kinematic setting of the actual continuous structure of Fig. 1 is revisited herein. This structure is replaced with a system of N discretized regions, whose nodal movements are controlled by the β degrees of freedom. Considering α static indeterminacy of the structure and S plastic rotational deformation at the nodes, the kinematic indeterminacy number takes the form β = S − α.  www.nature.com/scientificreports/ It can be shown that, when β independent nodal velocities q j j = 1, 2, . . . . . . , β are released, there exists a kinematically consistent independent member deformation rates ẋ h (h = 1, 2, . . . .., 2N) ; indicated in Fig. 2, velocities related to center of gravity of mass u k (k = 1, 2, . . . .., γ ) ; indicated in Figs. 3 and 4, and load point velocities δ = 1, 2, . . . .., n . Hence, the nodal setting of kinematics must be written: where the coefficient matrix is constant, provided that the motion falls within small displacements.
Material model. The constitutive relation between the stress-resultant S i 1 (bending moment M i ) and the strain-resultant rate ṡ i 1 (rotation rate θ i ) at the control node section i(i = 1, 2, . . . , χ ) of the discrete model is shown in Fig. 4. The vectors y * and ẋ * are the respective plastic potentials and the plastic multiplier rates, both controlling the plastic yielding at the critical section. X * defines the vector containing the plastic capacities of all critical sections. Finally, it must be noted that the plastic potentials, the plastic multiplier rates, and the plastic capacities are non-negative at +i or −i direction.
In view of Fig. 4 and the above description, the nonholonomic plasticity relations can be straightforwardly approximated 61 as where N is the vector of identity matrix I in the form I −I .
To develop a governing mathematical system, it is necessary to provide the following transformations that link the independent member forces X and the stress-resultants S, and their dual deformation rates ẋ and strain resultants ṡ: The mathematical formulation. Before developing a governing system based upon the foregoing relations, it is essential to select an approximate time-integration scheme that can determine a numerical solution sequentially. Hence the independent nodal velocities and the nodal displacements are advanced from time t n . to t n+1 using Newmark time integration scheme, that is, in which α , γ are the integration constant affecting the stability and the accuracy of the numerical process. By rewriting (10a) where the integration constants are: Stability and robustness of rigid plastic dynamics seem to be ensured when the average acceleration scheme is constant, in which α = 0.25 and γ = 0.5.
Considering Eqs. (1) to (9) at time t = t n+1 , and introducing Newmark scheme (11) to (13) the governing system can be written in the form The right-hand side sub-vector Y n+1 of (14) and the mass matrix M q are given by: The formalism (14) to (17) has a mathematical structure of a linear complementarity problem (LCP). It is important to appreciate that the velocities q n+1 for the next iteration are given directly, the accelerations q n+1 should not be found from (11) since it does not satisfy the kinetic law (2). Instead, by continuing to satisfy that law, together with (1) and (3), www.nature.com/scientificreports/ This can be solved uniquely for q n+1 for appropriate modeling having positive definite M q . Of course, for motion caused by the application of all impulses at the same time, t = t 0 , the load vectors in (18) and (20) are Initiation. For the impulsive loading, the LCP formulation (14) to (17) requires a velocity profile to initiate the evolutionary process. Although the initial displacements q 0 at time t = t 0 are directly prescribed, and usually q 0 = 0 , the initial velocities q 0 must be inferred from the prescribed magnitude and distribution of the initial loading impulses . Integrating (2) with respect to time, over the small interval [t 0 − ε, t 0 + ε] containing the impulsive loading, substituting (1) and (3), and taking the limit The first integral vanishes because the element end moments X for all elements are assumed to be bounded by the finite plastic moment values X * . The second integral gives the velocities q 0 at t = t 0 acquired from rest due to the initial impulses , which are given by the third integral. Thus which allows q 0 to be solved uniquely for modeling having positive definite M q . Now consider the evolutive process (14) to (17), which is not self-starting. So a self-starting subroutine is developed which calculates the relevant accelerations in (18). For determining these accelerations, the vector of plastic multiplier rate ẋ * is partitioned into the subset Y yielded (active) section and the subset R rigid (nonactive) sections as.
The yielded y and rigid r subscripts, shown in (23) and (24), can be employed to derive an acceleration relation of Tamuzh 62 Thus, in view of (1) and (25) to (29), the governing system 28,29 of acceleration can be established as If the structure is coerced into motion by an initial impulse or by an initial velocity field, the vector of initial loading is null 0 = 0 at the start of motion. Accordingly, the set of Y active plastic hinges ẋ * y > 0 can be easily deduced.
To establish the initial governing system (30) to (33) of the impulsively loaded structure, the yield surfaces activated by the initial velocities q 0 must be determined. More formally, it follows from (4) to (7), if q 0 is found from (22), ẋ * can be defined using the linear program Notably, by incorporating the above linear program in a newly developed subroutine, it turns out that most of the numerical difficulties that followed previously 30 due to the refinement of mesh disappear.
Plastic unstressing. The evolutive formulation (14) to (17) is unable to determine plastic unstressing within an increment ∆t. In order to sustain the accuracy and stability of the numerical process, another subroutine is included that determines the unstressing time within the interval ∆t [28][29][30] . Once the unstressing is detected, the www.nature.com/scientificreports/ evolutionary process is temporarily terminated, the multiplier rates are partitioned into active Y and in-active R sections, and the system (14) to (17)  Problem statement. Figure 5 shows the beam subjected to the action of a uniformly distributed vertical impulse of total value I which produces a uniform initial transverse instantaneous velocity V 0 = I/2mL . Let the length of the beam be 2L with uniform mass per unit length m and let the plastic moment capacity be M p for bending of either sense in the vertical plane.
Theoretical beam response. A continuum solution for the rigid-plastic beam of Fig. 5 has been presented by Symonds 63 , and in some detail by Jones 34 . Just after the instant of loading, the cross-sections directly over the supports have no initial velocity and no initial displacement; the rest of the beam acquires an initial velocity V 0 at no initial displacement, Fig. 5. Phase 1 of the motion begins with plastic hinges forming within the beam span and then traveling towards the midspan section. This velocity profile is illustrated in Fig. 6. Phase 2 of the motion initiates when the traveling plastic hinges meet at the midspan and this phase terminates when all the kinetic energy is consumed by the stationary hinge at the midspan. The velocity distribution of Phase 2 is shown in Fig. 7.
A typical requirement of the theoretical solution is the knowledge of the velocity profile during any phase of motion. The proposed linear program (34) identifies and calculates the velocity profile automatically.   To initiate the numerical sequence of LCP formulation, the initial velocity profile is required. In the case of lumped mass elements, this is rectangular shaped having a value of V 0 = I/2mL at all the masses, excepting those directly over the supports. In the case of uniform mass element, the vector q 0 of initial velocities is calculated from (22) with a discrete vertical impulse I/N applied at the centroid of each finite element. Given these initial velocities, the linear program (34) generates an optimum velocity profile shown in Fig. 8, which the LCP solver can handle very effectively.
By examining the results offered by LCP obtained for both lumped-mass elements (Table 1) and uniform-mass elements (Table 2), it is clear that the lumped -mass element presented a very small error. These tables present the results for the non-dimensional central displacement W f = W f /L (mL)M p /I 2 and the non-dimensional cessation time t = M p t 1 /IL for a range of finite elements. Likewise, the range of finite elements is employed to achieve convergence, Fig. 9.     Figure 10 shows the evolution of ξ , the distance of traveling plastic hinges from midspan, throughout Phase 1. It is observed that the LCP results agree closely with the theoretical solution (Jones 34 ) for the first phase of motion.
The ensuing dynamic bending moment, initiated from the stationary conditions, is shown in Fig. 11. Figure 11a shows a very small error in the bending moment profile between the lumped mass element and continuous mass elements. For interest, it may be noted that in the previous study 30 the non-dimensional bending moment showed substantial error at the start, which occurred due to considering an incorrect velocity profile. This error is significantly reduced by using a compatible velocity profile that is generated from proposed the linear program (34). Phase 2 begins at t = M p t/IL = 0.083 when the traveling plastic hinges meet at midspan. The corresponding bending moment distribution, Fig. 11d, then remains unchanged throughout Phase 2, in which the motion results from the active hinge only at midspan. The motion terminates at t = M p t/IL = 0.25 when the bending moment distribution jumps to the null state. In contrast to the result obtained from the uniform-mass model, it is seen that the lumped-mass model tracks the changes of the dynamic bending moment accurately throughout the whole motion.

Effect of transverse shear on simply supported beam
The response of the simply supported beam, Fig. 5, is analyzed by employing finite plastic shear deformation at the support locations. This problem is used as a reference solution to highlight some limitations of the kinetic and kinematic descriptions of (1) and (2) respectively. These descriptions are improved and the numerical solution is compared with the exact solution of Symonds 50 Bending-shear material model. A special examination is made of the effect of shear and bending response controlled by the rectangular yield curve. Figure 12 shows the rotational and the shear deformationrates and their representation by the independent deformation-rates.
The plasticity relation (4) to (7) is still applicable when the shear force interacts with the bending moment, Maier 61 . It is convenient to develop the interactive rectangular yield criterion at a critical section i in a threefold sense: The yield criterion;      www.nature.com/scientificreports/ The Association Rule: Figure 13 relates the bending moment S 1 and the shear force S 2 , at a single critical section i , using an approximate rectangular yield curve. The simplicity of the rectangular yield criterion is that, if only y +i * 1 = 0 , a simple plastic hinge having a positive value of θ i is activated, and if only y −i * 2 = 0 , a shear slide having a negative value of γ i is activated. Interaction occurs only at a corner of the yield surface, such as y +i * 1 = 0, y −i 2 = 0, for which plastic deformation may occur simultaneously at the section in a plastic hinge with positive θ i and a shear slide with negative γ i .

Representation of kinetics and kinematics as the nodal description.
If the total number of plastic rotations and shear deformations that can be involved in the motion of the beam is S , then the kinematic indeterminacy of the general velocity profile is β = S − α . Unlike the previous investigation 30 , it may happen that β is greater than the number of independent velocities associated with the motion of the nodes of the structural model, as may now be demonstrated. Surely, this investigation extends the previous study 30 toward a more general description of kinetics and kinematics.
Equations (8), (9) need some modification to link the bending and the shear forces and corresponding strain resultant rates S and ṡ to the independent element forces and deformation rates X and ẋ . For the element, with critical Sects. 1 and 2, as in Figs. 1 and 12, or, more generally, for element m containing critical sections i and j , the contributions of this element to the assembly of all elements in Eqs. (8), (9) are   (36) and (38), the matrix product TN for the whole element assembly is obtained.
If all system nodes are constrained not to move, then ẋ m = 0 for all elements. Despite this, the n = 4 plastic strain-rates for element m , ṡ i and ṡ j in (38), can still be chosen with some freedom to permit rigid-body motion of that element. The rank of the matrix in (38) is r = 2 , and the number of degrees of freedom in the rigid-body motion is given by β m = n − r = 2 . If an element has fewer than n = 4 possible strain-rates, and the relevant matrix column(s) are deleted, then the rank and β m can be similarly determined. It follows that β = β N + β M , where β N is the number of independent velocities of the nodes, those velocities being the components of q N , and β M is the number of independent velocities in the elements permitting internal rigid-body motions, those velocities being the components of q M . Now, the vector of strain-rates ṡ at all critical sections can be expressed in terms of the complete set of β independent velocities q . Thus, ẋ in (1) can be replaced by ṡ, and so, by kinetic-kinematic duality, X in the β = β N + β M Eq. (2) can be replaced by S, Equations (40) and (41) could be used to construct LCP (14) in terms of S n+1 . However, reverting from ṡ to ẋ , with T obtained from collecting together (38), for all elements m, which follows from the fact that ẋ = 0 whenever q N = 0 . With (42), the second equation set in LCP (14) is uncoupled from the velocities q M . Also, replacing S by X in (41), using (42), In general, when β M > 0 , the first set of equations in LCP (14) consists of β equations involving the components of both q N and q M . However, when the masses are lumped at the nodes, submatrix A T dM in (43) is null. When the loads are also applied at the nodes, submatrix A T 0M in (43) is also null. Then the first set of equations in LCP(41) is automatically condensed to a set of β N equations, uncoupled from the velocities q M . Automatic condensation does not apply when the masses are distributed in any way within the elements.
Problem statement. The problem set out in Fig. 5 is revisited in Fig. 14, but now the beam is assumed to have a finite transverse shear strength V p and plasticity to be governed by the rectangular yield criterion and associated flow rule of Fig. 13. Once more, motion commences with the uniform initial velocity distribution, V 0 = I/2mL.

Theoretical solution.
A theoretical treatment of a shear flexible simply supported beam of Fig. 14 has been given by Nonaka 64 and Jones 34 , following an earlier discussion of the problem by Symonds 49 . The character of the solution is much influenced by shear bending ratiȯ Finite shear resistance only at supports Figure 14. Impulsively loaded simply supported beam including finite shear resistance at the support location. www.nature.com/scientificreports/ The traveling plastic hinge phase is a suitable problem to illustrate the performance of the LCP solver. According to the theoretical solution, the overall response becomes more bending dominant as ν increases. For values of ν > 1.5 , the theory predicts a motion involving traveling plastic hinges in one of three distinct velocity phases.
The first phase of motion commences with plastic shear deformation at the supports of velocity Ẇ s , together with a pair of stationary plastic bending hinges at location x = ±ξ 0 , Fig. 15, where ξ 0 /L = 1 − 3/2ν . It is noted that the choice of ν > 1.5 would make 0 < ξ 0 < L . The phase terminates at time t 1 = 3mL 2 V 0 / 8ν 2 M P = 3IL/ 16ν 2 M P when Ẇ s = 0 and shear-sliding ceases.
The third phase of motion then commences when the ensuing deformation is concentrated in a stationary plastic bending hinge at the center of the beam, Fig. 17. This phase terminates at time LCP Prediction of response. The numerical LCP simulation is conducted using lumped mass elements and uniform mass elements. The beam has been discretized into 100 finite element mesh chosen to capture more faithfully the second phase of motion, the phase of traveling plastic hinges. The strength ratio ν = V p L/2M p = 1.72209 has been back-calculated from the value of ξ 0 . Table 3 compares the non-dimensional final central displacement W f = (W f /L)(mL)(M p /I 2 ) and the nondimensional time t = M p t/IL , for the both lumped mass and the uniform mass finite elements. The agreement between the theoretical and LCP results is quite reasonable, particularly when considering the lumped mass discretization. Figures 18,19,20,21 show graphically the evolutionary development in time t of the shear deformation W s , the central transverse velocity W , the central transverse displacement W , and the position ξ from midspan of the traveling hinges, all in non-dimensional form. Each figure shows the general ability of the LCP method to seek out automatically and follow the changing velocity phases over the whole extent of the motion. In Fig. 19, it will be observed that the lumped elements LCP solution for the transverse central velocity wanders slightly off track when the process is changing from the Phase 2 velocity profile to that of Phase 3. For the uniform mass LCP solution, this difference is hardly visible. In Fig. 21, the LCP solutions identify the position of the stationary plastic hinges in Phase 1 as ξ 0 = ξ 0 /L = 0.12897 and give the end time t 1 of the phase with good accuracy   www.nature.com/scientificreports/ (Table 3). This figure then shows the effective tracking by the LCP solutions of the position ξ of these hinges as they travel towards midspan in Phase 2, but it also shows the build-up to a small but slightly more significant error in the predicted end time t 2 of the phase ( Table 3). The bending moment distribution along the non-dimensional distance ζ /L of the beam for N = 100 elements is compared with the theoretical result, Fig. 22. Excellent agreement is observed for both choices of discretization, unlike the previous distribution in Fig. 11a. As in the previous case, the motion ceases at the time t = M p t/IL = 0.25.

Remarks on rigid plastic analysis and Lemke's Algorithm
It may be needless to say that the rigid-plastic dynamic analysis is adequately expressive of the true behavior of the concrete and the steel structures, provided that the total input energy transmitted to the structure is significantly larger than the maximum stored elastic strain energy. By excluding all the elastic deformation, rigid-plastic dynamics allow gaining insights into the precise mechanisms responsible for dissipating the plastic energy. The rigid-plastic theory constitutes an intuitive appreciation of plastic deformation, whereas the elasto-plastic theory, www.nature.com/scientificreports/ although embracing all relevant parameters, tends to define less sharply the evolution of the plastic deformation mechanism. Due to numerical instability, not many computational procedures are available that can determine an intuitive appreciation of the evolutive dynamic rigid-plastic deformation. Most of the commercially available finite element analysis software, although embracing all relevant parameters, tends to offer the evolution of the plastic deformation that is smeared by elasticity. In this context, the current study is concerned with adjusting and improving the previously developed systematic numerical procedure in the form of a linear complementarity problem (LCP) 30 for the impulsively loaded rigid-plastic beam problems.
In the case of an impulsively loaded beam subdivided into continuous mass elements having infinite resistance to shear deformation, the shear force causes the initial velocity profile to zig-zag like a tooth profile of a saw. Such a profile happens to be taxing on the LCP solver presenting considerable numerical difficulties, particularly with larger discretization. A workaround proposed previously 30 to overcome this difficulty was to use the lumped mass velocity profile instead. The proposed linear program (34) manages to overcome this intrinsic limitation of the LCP approach by providing a geometrically compatible velocity for the stable implementation of the algorithm. This implementation of the linear program in the Lemke solver 65-68 has proven to be a much more robust technique for simulating refined beam discretization.
Nevertheless, numerical instability issues re-emerged in the second investigation of the impulsively loaded simply supported beam in which the dynamic response involves shear deformation at the simple supports and bending deformation elsewhere. The previous description of kinetics (1) and kinematics (2) suffered irregular convergence when applied to the bending-shear interaction problem due to unconstrained rigid body motion. Therefore, modification is made to the kinetic and kinematic description, shown in (40) and (41) respectively. This strategy helped to significantly improve the convergence of the LCP solver.

Conclusions
This study is concerned with adjusting and improving the previously developed numerical procedure in the form of a linear complementarity problem (LCP) 30 . Numerical difficulties concerning the previously investigated impact problem 30 can be summarized as follows: • The use of the actual initial velocity profile in the continuously discretized beam can lead to severe numerical difficulties and inaccurate prediction of the impacted beam having infinite resistance to shear deformation  www.nature.com/scientificreports/ • Rigorous models involving bending-shear interaction can lead to numerical instability issues if the kinematic descriptions are not defined in terms of independent velocities of the nodes and independent velocities permitting the internal rigid body motions.
In the context of the above limitations, extra support is provided to the LCP formulation to make it computationally robust. Thus, a linear program (LP) is developed that determines an optimized initial velocity profile, which is necessary for initiating the dynamic response. For illustrative purposes, the improved LCP formulation is applied to an example of an impulsively loaded simply supported beam having finite bending strength and infinite shear strength. In comparison with the available continuum solution, the LCP analysis generates the complete evolution of all response parameters with reasonably good accuracy. However, the LCP solver becomes susceptible to breakdowns when a second illustrative numerical test is performed on the same beam, but with finite shear resistance at supports. Enhancing the kinetic and kinematic description to permit the internal rigid body motion renders the LCP solver much more robust and overcomes all the numerical difficulties. The numerical results are again compared to the continuum solution of the interaction problem borrowed from literature. The accuracy and computational efficiency of the LCP bending shear interaction model are demonstrated by accurate prediction of the bending moments, deflections, and failure mechanisms. In this regard, it is worth stressing, that the scope of this enhanced formulation can be widened to include the interaction effects of torsion and axial force.