An Improved Linear Complementarity Solver for the Dynamic Analysis of Blast Loaded Structures

The linear complementarity problem (LCP) approach, expedited by using the simple rigid–plastic theory, has been utilized successfully in predicting the numerical response of the ductile steel or concrete structures subjected to short-duration, high-intensity dynamic loads. The current study attempts to improve the computational stability of this powerful technique while determining the response of skeletal structures under blast loading. The performance of the Lemke LCP solver is amplified by introducing an automatic time-stepping scheme to efficiently trace the complex dynamic response using either lumped mass or continuous mass discretization. The computational efficiency of this solver is tested against carefully chosen three numerical examples, and the acquired results are in good agreement with the derived closed-form solution and results from other sources.

Most of the numerical models used in research require a high level of expertise and are computationally expensive (Jones, 1976). This demands a simple numerical method that can predict the blast loading response with reasonable accuracy.
The simplicity of rigid-plastic approximation has been known for several decades as an analytical tool defining and crystallizing concepts that can be extended to a wide class of problems (Cennamo et al., 2017;Jones, 1986;Lee & Symonds, 1952a;Taylor, 1948). The rigid-plastic theory was first applied to dynamic problems in 1949 (Taylor, 1948) and systematically studied in 1952 (Lee & Symonds, 1952b). These studies subsequently yielded vast literature on the investigations of structures submitted to extreme dynamic loading (Bleich & Shaw, 1960;Jones & Shen, 1993;Ling et al., 2017;Lowe et al., 1972;Mehreganian et al., 2019;Menkes & Opat, 1973;Parkes, 1955Parkes, , 1958Symonds, 1967;Symonds & Frye, 1988). Yet, it is noteworthy that each closed-form theoretical solution requires postulating a kinematically admissible velocity profile for the evolution of displaced configuration.
Therefore, most procedures incorporating the simple rigid-plastic theory are oriented towards obtaining specific results for specific problems and demand personal judgment for getting the best results.
The numerical investigation of the dynamic response of rigid-plastic structures appears not to have been examined extensively. The commercially available finite element softwares, such as ABAQUS, ANSYS, etc., experience severe mathematical complications when a purely rigid-plastic model is adopted. From this precedent, mathematical programming provides an approach to incorporate the rigid-plastic theory into a unified formalism that is not specific and problem-oriented. Furthermore, mathematical programming (Dantzig, 1998;Dantzig et al., 1955;Guzas & Earls, 2011;Lemke, 1978;Murty, 1983) has a broader application in various specialized fields of engineering, such as robotics (Nuseirat & Stavroulakis, 2000), fluid simulation (Andersen et al., 2017), and agriculture (Garrido et al., 2020). It has the potential to proffer a finite element-based numerical formulation (Maier, 1984;Martin, 1964) that facilitates any distribution of mass, spatial placement of applied loading, and temporal variation of the associated load pulses. Moreover, once the physical modeling decisions have been taken, a complete algorithmic or completely automatic solution procedure can be obtained. The benefits of mathematical programming have been recognized for more than 60 years. Extensive surveys of the use of the mathematical programming application in engineering plasticity have been reported by Maier (Maier, 1984), Maier and Munro (Maier & Munro, 1982) and Maier and Lloyd Smith (Maier & Munro, 1982;Smith, 1974).
Tamuzh (Tamuzh, 1962) first offered a procedure to determine the response of rigid, perfectly plastic continua using kinematic minimum principle. Later Capurso (Capurso, 1972) used this principle to formulate a mathematical programming problem capable of tracing the response of rigid-plastic framed structures submitted to short-duration pulse loads. The capability of such mathematical programming problem can be enhanced by incorporating the effects of strain rate and large displacements in rigid-plastic material model (Cennamo et al., 2017;Jones, 1986;Lee & Symonds, 1952a;Taylor, 1948). In recent times, Patsios and Spiliopoulos (Patsios & Spiliopoulos, 2018) have analyzed structural frames using the mathematical programming method. Moreover, Milani and Tralli (Milani et al., 2009) have employed mathematical programming to model the behavior of masonry walls. Later, Portioli has used a similar model for the dynamic and pushover analysis of masonry structures (Portioli, 2020). Although mathematical programming has the potential to solve problems in dynamic plasticity (Khan et al., 2013;Rodigari et al., 2019), this tool has not been exploited to any great extent (Milani et al., 2009;Wu et al., 2020).
In a companion paper by Khan et al. (Khan et al., 2013), a mathematical programming formulation, called the Linear Complementarity Problem LCP (Lloyd Smith & Sahlit, 1991;Smith, 1990), was developed. In that study, the LCP solution was employed as a simplified design procedure for the assessment of structures under impact loading. This type of loading required an initial velocity profile for starting the Lemke Algorithm (Khan et al., 2021). However, it was shown (Khan et al., 2013) that such problems posed considerable numerical difficulties in finding solutions for simple bending-only models. These numerical problems were avoided by considering the bending-shear interaction in the model. In a recent research (Khan et al., 2021), the Lemke algorithmic solution was tested by investigating several engineering problems; these included a simply supported beam subjected to rectangular pulse load and a portal frame subjected to triangular pulse load. A comparison of the LCP approach with the theoretical and the numerical solutions using lumped mass discretization showed excellent agreement. Nevertheless, there were occasions, particularly in the analysis of the portal frame problem, when the algorithm would become unstable producing illogical results. One remedy to reduce this instability was to use a smaller time step that is specified at the outset of the evolutionary process. In the current study, this numerical instability of the Lemke solver is investigated, and remedial measures are proposed. Specifically, the instability becomes more pronounced with the continuous mass discretization and mesh refinement. It turns out that this instability can disappear if the time step is allowed to adjust automatically, especially in highly nonlinear dynamic problems. A newly developed time-step controller subroutine is incorporated within a MATLAB program that reduces the increment size repeatedly until stability is achieved. It also overcomes the instability problems encountered when refining the finite element mesh (Khan et al., 2013). Three illustrative examples are carefully selected to test the accuracy and efficiency of the improved LCP solver. The first two examples are comparable to that of the recent investigation (Khan et al., 2021), but with the increased non-linearity to the dynamic response. The third verification example investigates the experimental study of Zhang et al. (Zhang et al., 2013) involving the dynamic response of a semi-clamped RCC beam.

Research Significance and Assumptions
Computation of structural response under blast loading is important for reducing the risk and ensuring the safety of people. This demands a fast algorithm that can efficiently handle material nonlinearities and give accurate results. There exist limited numerical methods that can be applied to capture the complex dynamic behavior of structures subjected to blast loading. Commonly used methods of determining the response of structures to blast loading are simulations in commercial finite element softwares, such as ABAQUS, ANSYS, etc. However, due to their higher computational cost and requirement of skills to simulate complex models, an efficient formulation is needed that has the important practical benefit of speed and simplicity. Inspired by the promise shown by the previous investigations, rigid-plastic structures under impact (Khan et al., 2013) and Lemke Algorithm for the rigid-plastic response (Khan et al., 2021), the present study is a more systematic exploration of this powerful approach to efficiently predict the response of skeletal structures under blast loading. A predicament encountered previously (Khan et al., 2021) was a tendency of the LCP algorithm to crash when a large number of continuous mass finite elements were employed in the modeling. Such problems are resolved in the current study when a newly developed automatic time-step controller is incorporated. Three computational examples are presented to prove the accuracy and efficiency of the updated LCP solver with low-order finite elements.
The following assumptions are made in the proposed LCP formulation: 1. The initial energy is so large that the elastic response can be neglected 2. A skeletal structure is envisaged as an assemblage of elements connected with other elements at stations called nodes. 3. The development of plasticity is confined to discrete points or nodes on the structural element. No account is taken of the spreading of plasticity, either in the direction of the locus of cross section or within the structural element between the discrete points.
4. The mass properties of a structure are defined at discrete points. Thus, it is idealized that massless elements connect these lumped masses. 5. The structural system is considered as rigid-plastic, incompressible, and inextensible. Therefore, the effect of strain rate and strain hardening are ignored. 6. The deflection of the structural system remains small during the dynamic response.

Proposed Dynamic Rigid-Plastic Model
In the current section, the fundamental vectorial conditions, namely the kinetics, the kinematics, and the material constitution, characterizing the behavior of rigid, perfectly plastic structural systems undergoing dynamic disturbances due to blast loading, are combined consistently. The structure is envisaged as an assembly of discrete finite beam elements, such as a discretized beam shown in Fig. 1, whose structural mass is either concentrated at the boundaries of the elements or continuously distributed along the elements. These continuous mass elements can cause instability issues in the recently reported Lemke solver (Khan et al., 2021). Thus the efficiency and stability of the solver are improved by implementing an automatic time-stepping algorithm with the Lemke solver.

Representation of Kinetics and Kinematics as Nodal Description
The structure in Fig. 1 is formed from nodal substructure having the generalized stresses S (i.e., bending moment M) and the generalized strain rates ṡ (i.e., rotation rates θ ) imposed at the critical sections. The structure is subdivided into N finite elements, in which the independent movements of the interconnecting nodes are governed by β degrees of freedom. Any kinematically Page 4 of 25 Khan et al. Int J Concr Struct Mater (2022) 16:47 consistent velocity distribution or profile may be specified completely in terms of β independent nodal velocitiesq j j = 1, 2, . . . . . . , β . For an assembly of the inextensible planar elements, with α static indeterminacy and S plastic rotational deformations occurring at the element extremities, the kinematic indeterminacy number can be established asβ = S − α.
When each of the β independent nodal velocities q is released, a velocity profile is generated; Figs. 1, 2, 3. Geometric consideration of the velocity profile can help in determining the independent member deformation rates ẋ h (h = 1, 2, . . . .., 2N ) , the velocities related to the center of gravity of mass u k (k = 1, 2, . . . .., γ ) , and the load-point velocities δ l (l = 1, 2, . . . .., n) . Hence, the nodal representation of the kinematic equations has the form: where the coefficient matrix is constant, provided that the motion falls within small displacements.
Let the structure be subjected to n discrete timedependent loads l (l = 1, 2, . . . .., n) applied at the nodes. By employing the D' Alembert principle, during every instant of the accelerated motion of a structure, the applied loads and the inertia forces µ k (k = 1, 2, . . . .., γ ) are in equilibrium with the independent member forces X h (h = 1, 2, . . . .., 2N ) . Corresponding to the independent nodal displacements, the nodal forces of constraint Q j j = 1, 2, . . . .., β are applied. For the satisfaction of the dynamic equilibrium, it is necessary that the constraints Q j must vanish, giving the nodal kinetics description for the assembly of all elements: where T denotes transposed coefficient matrix. It may be observed that (1) and (2) satisfy the adjoint relationship of kinetic-kinematic duality. (1) The independent relations (1) and (2) have no causeeffect relationship between the kinetic and kinematic variables. Nevertheless, the relation implicitly links the inertia forces µ k (k = 1, 2, . . . .., γ ) , located at the mass centroid, with the corresponding centroidal accelerations ü k (k = 1, 2, . . . .., γ ) of the actual motion of the system. In this inertial law, the diagonal mass matrix m k (k = 1, 2, . . . .., γ ) constitutes the mass or moment of inertia related to the corresponding centroidal accelerations.

Material Model
The cause-effect relation between the stress-resultant S i 1 (bending moment M i ) and its dual strain-resultant rate ṡ i 1 (rotation rate θ i ) at critical section i, (i = 1, 2, . . . , χ) , is illustrated in Fig. 4. The yielding at the critical section i is defined by two variables, i.e., the plastic potential y * and the plastic multiplier rate ẋ * . Notably, X * collects the moment capacities for positive or negative bending. Suppose if X +i * ≥ 0 , so that the critical section i has deformed plastically ( y +i * = 0,ẋ +i * ≥ 0 ) and therefore the stress-resultant S i 1 is positive. A similar argument applies to X −i * ≥ 0, y −i * ≥ 0,ẋ −i * ≥ 0 when S i 1 is negative. ( Centroidal velocities in a lumped mass system.
Centroidal velocities in a continuous mass system.

Fig. 4
Variables describing a simple flexural plastic hinge.
Page 5 of 25 Khan et al. Int J Concr Struct Mater (2022) 16:47 The plastic potential y * and the plastic multiplier rate ẋ * are coupled together via complementarity condition ( y +i * ẋ +i * = 0 ). This condition ensures the irreversible nature of plasticity. It means that the plastic deformation will only occur (say ẋ +i * > 0 ) if the corresponding yield limit is attained y +i * = 0 . Whereas, if the yield limit is not attained (say y +i * > 0 ), then the plastic deformation cannot be active ẋ +i * = 0 . Now the rigid-perfectly plastic constitutive relation for each critical section i is written in matrix notation: where N is the matrix defining the exterior unit normal to the yield function. More specifically, N = I −I and I is the identity matrix.
It is of interest to develop a governing mathematical system that couples the nonholonomic (Maier & Munro, 1982) constitutive relations, (4) to (7), with the kinetic and the kinematic relations, (1) to (3). Fig. 1 clearly illustrates that the independent member forces X and the independent member deformation rates ẋ can be defined by the respective stress-resultants S and the strain-resultant rates ṡ . These can be collected for all the constituent finite elements:

The Mathematical Formulation
The vectorial relations (1) , (2), (3), together with the triad of complementarity conditions (5), (6), (7), can be combined into a set of second-order differential equation with respect to time. Nevertheless, this set is made more complex by the complementarity conditions. As no closed-form solution is known to this kind of mathematical problem, adopting a numerical solution appears reasonable. Therefore, a time marching scheme is introduced in order to allow the solution to be advanced from a time station t n to t n+1 = t n + t , where subscript n is an integer defining successive discrete time stations and t is the intervening time increment. Then the centroidal velocities and accelerations can be expressed in the Newmark's time-integration scheme: and in which integration constants are It is found after thorough investigations (Khan et al., 2013) that suitable results are obtained for rigid-plastic dynamics if α = 0.25 and γ = 0.5.
The right-hand side sub-vector Y n+1 of (17) is given by: and the mass matrix M q is given by relation: The approximating governing system (17), (18), (19), (20) has a mathematical structure of a linear complementarity problem (LCP). It may be noticed that the variables [ẋ * , y * ] are restrained into the complementary pairs, whereas the leading sub-matrix related to variables [q, X] is negative semi-definite. In this work, the governing system is solved efficiently by Lemke algorithm due to its simplicity and robustness.

Initiation
The incremental numerical process shown in (17), (18), , (20), representing the evolutive sequence of the dynamic response, is not self-starting. Therefore, it is necessary to establish a subroutine to calculate the relevant accelerations at a certain instant. These accelerations are of particular relevance at the commencement of the motion and the deactivation of the previously active section. At the initiation of the response, the vector of plastic multiplier rates ẋ * is separated into Y yielded sections and R rigid sections. So, these subsets of the multiplier rates ẋ * can be expressed as: Using Tamuzh principle (Tamuzh, 1962), the discrete law (4) to (7) can be written in the form Notably, the relation (25) to (29) represents the admissible acceleration field that is derived from differentiating the rigid-perfectly plastic constitutive relation (4) to (7).
Differentiating with respect to time the kinematic relation (1), and together with the admissible acceleration field (25) to (29), Sahlit (Lloyd Smith & Sahlit, 1991) re-established the governing system in terms of accelerations: with variables q,ẍ * y , X unrestricted. The previous study (Maier, 1984) of impact required an initial velocity field to coerce the dynamically loaded structure into motion, so the vector of initial loading was null 0 = 0 at the start of motion. The current investigation has a load vector = 0 , while the prescribed initial velocity is null q 0 = 0.

Automatic Time-Stepping Algorithm
The LCP formulation (17), (18), (19), (20) has stability restrictions on the size of the time step that is specified at the outset of the evolutive response. Lemke algorithm (Khan et al., 2021) solves the LCP formulation very efficiently and robustly, provided the fixed time-step size is small enough to ensure stability. A smaller time step defined at the outset can provide a detailed description of the whole response, but is computationally expensive. Selecting a fixed step size is very challenging (Khan et al., 2013) in the case of the continuous mass model of Fig. 3, especially as the mesh is refined. The continuous model with a bigger time step can cause the LCP solver either to crash or converge to a physically meaningless solution. A possible solution to overcome this problem is an automatic selection of time-step size, which can maximize the numerical accuracy while minimizing the computational cost.
The proposed automatic time-stepping algorithm is implemented via a routine in the LCP iterative procedure as shown in Fig. 5. This scheme adjusts the time step when it detects the anomalous value of plastic multiplier rates ẋ * occurring over several consecutive time intervals. These anomalous values of plastic multiplier rates are given in Fig. 5. Thus, the evolutive sequence of the dynamic response is terminated temporarily at the instant when this anomaly is detected and restarted with the increment size shortened. The stepping algorithm can reduce the increment size repeatedly if the anomaly remains. Once the anomaly is resolved, the time increment is lengthened to its original value. This automatic control allows substantial gains in the numerical stability of the LCP algorithm. If unstressing is detected within the increment ∆t, another subroutine calculates the unstressing time instant t n+ε within the interval ∆t (Khan et al., 2013). Subsequently, the structural variables at t n+ε are determined, and the evolutive process of the LCP system (17), (18), (19), (20) is reinitiated with t n+ε as the starting-time.

Case Study: Triangular Pressure Pulse on a Propped Cantilever Beam
The underlying behavior of a pulse-loaded propped cantilever beam is illustrated in this section; Fig. 6. A recent study (Khan et al., 2021), in which a rectangular pulseloaded simply supported beam was examined using the Lemke Algorithm, is extended to a linearly decaying pulse and asymmetric support condition. The primary aim is to validate the automatic time-stepping algorithm by reference to a suitably complex problem. More specifically, the support and the loading condition of this case study problem are chosen to enhance the complexity of  Page 8 of 25 Khan et al. Int J Concr Struct Mater (2022) 16:47 the traveling plastic hinge phases. The exact solution is derived here to give insight into the associated physical phenomena and to demonstrate the efficiency of the improved Lemke algorithm.

Problem Statement
Let a rigid-perfectly plastic propped cantilever beam of length L and uniform mass per unit length m be subjected to a uniformly distributed linearly decaying pulse, shown in Fig. 6. The pulse can be written as: Let the fully plastic bending moment of the beam be M p , while the effect of shear force on yielding be neglected. A sequence of the response phases can be performed based on the magnitude of the force P 0 ; Fig. 7. Whenever the applied load is less than the static collapse load P c , the rigid-perfectly plastic beam remains static and undeformed (Jones, 1990).
The proposed theoretical solution of the problem shown in Fig. 6 is discussed below. The description of the phases is also given in Fig. 7.

Theoretical Response of the Beam, P c ≤ P 0 ≤ 3P c
Phase 1 of motion An analytical solution is desired using the transverse velocity distribution shown in Fig. 8. The beam is idealized as two rigid arms connected by a stationary plastic hinge located at x = ξ * from the clamped end. It is found, if P c < P 0 < 2P c , the motion commences at t = 0 , and then ceases before the pulse terminates; that is, t < τ.
The bending moment is given by and (35) A stationary hinge is formed at ξ * = 0.585786438L, Fig. 8, which is similar to the static collapse case. Within the phase the velocity at ξ * is given by As mentioned earlier, if P c < P 0 < 2P c , the motion ceases before the pulse terminates. The beam will reach its final position when the central velocity Ẇ 1 = 0 , which occurs when: However, if 2P c < P 0 < 3P c , the motion continues until t > τ . Therefore, the first phase of motion concludes at t = τ , when the beam has an associated maximum transverse velocity: for t = τ and 2P c < P 0 ≤ 3P c , and the peak transverse displacement: for t = τ and 2P c < P 0 ≤ 3P c .

Phase 2 of Motion
The beam will be unloaded when the external pressure terminates at t = τ ; Fig. 6. If the transverse velocity of the beam at time t = τ conforms to (42), then the beam will continue to deform after time t ≥ τ until the remaining kinetic energy is dissipated at the plastic hinges. After unloading, the beam continues to deform and the velocity at the center of the beam is represented as with the corresponding maximum displacement is: The motion comes to stand still when Ẇ 2 = 0 , which occurs when:

LCP Prediction of Beam Response
Consider the propped cantilever beam shown in Fig. 6 that is impelled by a uniformly distributed linearly decaying force pulse. Two amplitudes of pressure pulses were considered, one with force magnitude P 0 = 1.5P c and non-dimensional duration τ = M p τ/IL = 0.458 , and the other with force magnitude P 0 = 2.5P c and non-dimensional duration τ = M p τ/IL = 0.275 , where I = 1/2P 0 L is the total impulse of the load. The rigid beam was discretized into 100 finite elements by idealizing the structure into either lumped mass elements (Fig. 2) or continuous mass elements (Fig. 3).
Investigation showed that the LCP offered promising results having small errors in most of the examined quantities. Table 1   Page 10 of 25 Khan et al. Int J Concr Struct Mater (2022) 16:47 both calculated for the force magnitude P 0 = 1.5P c and P 0 = 2.5P c . It is evident from the table that for P c < P 0 < 2.0P c the motion ceases before the pulse terminates at τ = M p τ/IL = 0.275 ; whereas, for 2P c < P 0 < 3P c the motion continues until after the pulse terminates. Hence, the LCP solution confirms that the dynamic response is characterized by the pulse magnitude. Fig. 9 shows the evolution of non-dimensional central displacement W 1 = (W /L).(mL)M p /I 2 at ξ * (Fig. 8) when P 0 = 2.5P c . It is apparent from Fig. 9 that the LCP results concur with the theoretical solution.
Figs. 10 and 11 show the bending moment distribution across the non-dimensional distance ζ /L from the clamped support of the beam. Once again, the accuracy and efficacy of the LCP solution are demonstrated.
All these results show that both the lumped mass and the continuous mass discretization substantially agree with the closed-form solution. Previously (Khan et al., 2013), the continuous mass elements incurred extreme difficulty, especially after mesh refinement. Therefore, the recent study (Khan et al., 2021) used only lumped mass discretization to solve various numerical examples of blast-loaded skeletal systems. The current investigation clearly shows that regardless of lumped or continuous mass discretization, mesh refinement is not an issue when an automatic time-stepping scheme is adopted. However, a suitable complex test involving the phase of  Fig. 9 Evolution of central displacement at ξ * (P 0 = 2.5P c ). Page 11 of 25 Khan et al. Int J Concr Struct Mater (2022) 16:47 traveling plastic hinges can provide a better opportunity for the robustness of the LCP method.

Theoretical Response of the Beam, P 0 ≥ 3P c
In order to furnish an elaborate numerical test of the improved LCP method, the phase of the traveling plastic hinges is investigated. This phase is fostered when the applied pressure pulse is of higher magnitude, i.e., P 0 ≥ 3P c . It transpires that for pressure pulses 3P c < P 0 < 6P c two traveling plastic hinges appear within the beam span and move inward; Fig. 12. These hinges coalesce before the pulse terminates. However, if the pressure pulse P 0 > 6P c , then the traveling hinges do not coalesce until after the pulse terminates.

Phase 1 of Motion
Initially, at t = 0 , it is postulated that three plastic hinges form, two within the beam span and the third at the clamped support. The transverse velocity distribution is in Fig. 12. Thus, for 0 ≤ x ≤ ξ * 1 , and With the imposition of the triangular pulse load on the beam two plastic hinges are formed, Fig. 12, at a distance of ξ 0 1 .
from the root of the cantilever, and ξ 0 2 (Fig. 12) from the other side.    Equations (51) and (52) show that the extent of the central zone ξ 0 1 ≤ x ≤ ξ 0 2 depends on the pressure pulse magnitude and its variation with time. Within this phase, the velocity is and the displacement is given by The bending moment is given by for 0 ≤ x ≤ ξ * 1 , and The beam continues to deform and the plastic hinges now move inwards. The velocity of these traveling plastic hinges at any time t can be determined from (55) Page 14 of 25 Khan et al. Int J Concr Struct Mater (2022) 16:47 for ξ * 1 ≤ x ≤ ξ * and for ξ * ≤ x ≤ ξ * 2 . (58) It is found that for pressure pulses 3P c < P 0 < 6P c the two traveling plastic hinges coalesce before the pulse terminates; whereas, if pressure pulse P 0 > 6P c , then the traveling hinges coalesce after the pulse terminates at t = τ . In the former case, Eqs. (50), (53), and (54) can be used; whereas, in the latter case, again (50), (53), and (54) can be used to get the transverse acceleration, velocity and displacement between 0 < t ≤ τ . For the subsequent interval τ ≤ t ≤ T 1 , when the pressure pulse is removed, the acceleration d 2 W 2 /dt 2 = 0 . Thus, Page 15 of 25 Khan et al. Int J Concr Struct Mater (2022) 16:47 at τ ≤ t ≤ T 1 and P 0 > 6P c , and the transverse displacement: at τ ≤ t ≤ T 1 and P 0 > 6P c . For 3P c < P 0 < 6P c , the juncture of the traveling plastic hinges may be found as ξ * = 0.585786438L with the velocity of the traveling plastic hinges made zero. These hinges coalesce at the interval T 1 .
When P 0 > 6P c , Eqs. (58) and (59) may be used to give the location ξ # 1 and ξ # 2 of the plastic hinges at time t = τ when the pressure pulse is removed. The subsequent position of the traveling plastic hinges can be found from the following equations: for ξ * 1 ≤ x ≤ ξ * and P 0 > 6P c at time τ ≤ t ≤ T 1 , and for ξ * ≤ x ≤ L and P 0 > 6P c at time τ ≤ t ≤ T 1 . See Figs. 8 and 12 for further clarity. The differential Eqs. (58) and (59) are solved using the Runge-Kutta fourth-order method. A little finesse is required to initiate the numerical solution, since the time t = 0 is a singular point at which the travel velocity of both the span hinges is infinite.

Phase 2 of Motion
Phase 2 initiates when the traveling plastic hinges meet at ξ * = 0.585786438L as presented in Fig. 8. Therefore, the remaining kinetic energy is consumed at hinges located at ξ * and the fixed end. This phase ends when the beam come to standstill, so the beam kinetic energy is fully dissipated.
Equations (64) and (65) can be rewritten for the interval when the pulse is removed: However, if P 0 > 6P c , the pulse terminates before the traveling hinges coalesce. So, Consequently, the transverse acceleration is: Page 17 of 25 Khan et al. Int J Concr Struct Mater (2022) 16:47 for P 0 ≥ 6P c when T 1 ≤ t ≤ T f , and The motion comes to a standstill when Ẇ 3 = 0 , which occurs when: for 3P c < P 0 < 6P c , and for P 0 ≥ 6P c .

LCP Prediction of the Beam Response, P 0 ≥ 3P c
For the LCP solution, the beam shown in Fig. 6 was discretized into n = 100 finite elements, for both lumped mass and continuous mass element types. Again, two amplitudes of pressure pulses were examined, one with force magnitude P 0 = 3.5P c and non-dimensional duration τ = M p τ/IL = 0.196 , and the other with force magnitude P 0 = 12.5P c and non-dimensional duration τ = M p τ/IL = 0.0549.
The LCP solution showed encouraging results with a small error in most of the examined quantities. Table 2 gives results for the non-dimensional central displacement W = (W /L).(mL)M p /I 2 and the non-dimensional time t = M p T /IL both calculated at various phase transitions. Fig. 13 indicates the evolution of the non-dimensional maximum central displacement at ξ * (Fig. 8) throughout Phases 1 and 2 of motion for the pulse magnitude P 0 = 12.5P c . It is evident that the LCP results agree substantially with the theoretical solution.
When the pulse load is applied suddenly att = 0 , two plastic hinges are formed simultaneously at x = ξ 0 1 and x = ξ 0 2 , with ξ 0 1 /L = 0.28706 andξ 0 2 /L = 0.79702 , Fig. 12, while a stationary plastic hinge is formed at the clamped end. Both sets of hinges at ξ 0 1 and ξ 0 2 start moving inwards towards a stationary pointξ * /L = 0.58579 . Figs. 14 and 15 show the path traced by the traveling plastic hinges for the applied pulseP 0 = 12.5P c . This phase terminates when the traveling hinges reach ξ * /L = 0.58579 at nondimensional timet 1 = M p T 1 Fig. 16 indicates the bending moment evolution of the beam at phase transition for both the theoretical and the LCP models. The bending moment distribution is for the decaying pulse with initial amplitude P 0 = 12.5P c . Again, both the LCP solution and the theoretical solution show substantially close agreement from the start until the end.

Performance of Improved Lemke Solver
The earlier investigation (Khan et al., 2013) exhibited a substantial error in various response parameters for the continuous mass element modeling. This error was attributed to the approximate velocity profile supplied to the LCP algorithm. The recent pulse-loaded beam investigation (Khan et al., 2021) was carried out for the lumped mass discretization only. Both of these investigations suggested that the spatial discretization requires sufficiently small time step and tolerance so as to allow the dynamic response adequately traced by the Lemke solver. The current implementation of the automatic time-stepping algorithm makes the solver efficient and stable by adjusting the time step whenever the instability occurs. The performance of the algorithm can be appreciated that both the lumped mass and continuous mass elements provide workable results with less than 1% error in most examined quantities. This reduced error in the case of continuous mass elements undoubtedly occurs because of the known initial load required to initiate the LCP algorithm. The LCP solution of the beams, investigated above in Subsections 4.3 and 4.5, is obtained considering the minimum tolerance of 1 × 10 -6 and the maximum time step of 4 × 10 -5 for capturing all phases of motion. The calculation time of the Lemke Solver is on average about 70 s. The problems without plastic hinge travel phase (Subsection 4.3) results in stable execution of the algorithm employing the initially specified time increment. However, the Lemke algorithm instability is prevalent if the motion involves a plastic hinge travel phase, thus requiring variation in time-step size to overcome this numerical difficulty. Fig. 17 shows the variation in the size of the time step during the dynamic response of the continuous mass discretized beam under the applied pulse load magnitude P 0 = 12.5P c . It transpires that three different time increment cases are possible, namely the initially specified time increment, unstressing time increment, and the automatic time increment; Fig. 17. The figure shows that the initially supplied maximum time step ∆t max is set equal to 4 × 10 -5 s, which sharply reduces to 4 × 10 -6 s at various intervals when numerical difficulties are detected by the automatic time controlling algorithm. All the other reductions in the time-step size indicate the instants of unstressing. The solver tends to reach the initial time step once the numerical instability is resolved.

Case Study: Portal Frame Subjected to Triangular
Pulse Load

Problem Statement
The second example of a single-story fixed-based portal frame is borrowed from the previous investigation (Khan et al., 2021), which compared the LCP solution with a finite element model through ABAQUS. The LCP solution of this problem is reexamined with continuous mass discretization and a refined mesh. Similarly, the ABAQUS model is revamped to reflect the actual behavior of a rigid frame. A schematic of the portal frame is shown in Fig. 18, which is fabricated from a UKB 610 × 305 × 179 in S355 steel, having a unit mass of 179 kg/m and a plastic moment capacity of 1910 kNm. A triangular pulse load of 6.55 N/mm 2 is uniformly applied on the left column of the frame for a duration of 0.0082 s.

LCP Response of Frame
The earlier (Khan et al., 2021) LCP investigation considered this portal frame discretized into 30 lumped mass elements. It was found that mesh refinement leads to the Lemke solver instability if the fixed time increment is not chosen appropriately. The current investigation has seen a marked increase in the solver instability if continuous mass discretization is adopted. This numerical difficulty is remedied by the automatic time-stepping algorithm that suitably changes the increment size based on the changes in dynamic response. Thus, the portal frame problem is revisited to test the efficacy of the time-stepping algorithm. This time the mesh is refined into 20 elements per member (60 elements altogether) by employing either the lumped mass discretization or continuous mass discretization. The LCP response obtained is shown in Fig. 19. The plastic hinges in this figure are expressed symbolically in Fig. 19. The motion evolves from the initial profile of Fig. 20a to the final profile of Fig. 20(d) and ceases before the termination of the pulse. Essentially, both the lumped mass and the continuous mass analyses run smoothly, returning the sway displacement of 39.3 mm and 41.5 mm, respectively, thus showing a coherent and consistent solution. Four critical sections are always active during each phase of motion, leading to a constant bending moment at every critical section; Fig. 20. Finally, it is essential to remark that the initially supplied maximum time step ∆t max to capture the dynamic response is set at 8 × 10 -5 s, Fig. 21. During the recurrent sequence of the Lemke Solver, this increment size is adjusted twice to 1.6 × 10 -6 s for ensuring algorithm stability. The overall MATLAB execution time for the frame problem is 15 s.

Response of the Frame Predicted by ABAQUS
Previously reported (Khan et al., 2021) ABAQUS model of the portal frame was developed by prescribing the maximum possible Young modulus to achieve rigidity of the connections and members. This approach of employing a high-magnitude Young modulus resulted in an anomalous shear failure at one of the hinges at the portal base. The 3D numerical model is redeveloped in ABAQUS Explicit by providing rigidity of the connections through suitable stiffeners; Fig. 22. The model is constructed using tetrahedral brick elements while considering four bricks along the thickness of the member to capture the bending behavior. The material properties assumed for the steel frame are density = density of 7850 kg/m 3 , the Poisson's ratio = 0.3 and Young modulus = 200 Gpa. The validity of the results is examined by making a comparison between the displacement time histories of the LCP frame with the ABAQUS approach, Fig. 22. In general, the LCP sway displacement of 39.3 mm in the case of lumped mass elements and 41.5 mm in the case of continuous mass elements concur with the results of ABAQUS 40 mm (Fig. 23). It is also observed that the location of plastic hinges in the final phase of motion (Fig. 20d) coincides with the Misses stress field of Fig. 24. Interestingly, the type of failure obtained from the ABAQUS model involves bending and shear effects at the portal base. This clearly shows that the transverse shear force has an important influence on the dynamic plastic behavior when the amplitude of the applied load is high. Thus the LCP solution involving bending and shear deformation (Khan et al., 2013) at the support base can be considered closer to the actual behavior shown in the ABAQUS model. However, the comparison of the LCP and ABAQUS sway displacements shows that the assumption of infinite resistance to shear deformation in the LCP model leads to satisfactory results. Moreover, the computational time required to solve the frame problem in ABAQUS is about 20 min, which is more computationally costly than the 15-s run time of LCP solver.

Case Study: LCP Prediction of Clamped Reinforced Concrete Beam
The dynamic response of a clamped reinforced concrete beam subjected to blast loading has been examined experimentally by many researchers. It is aimed in this section to compare the existing experimental results with the predicted results of the LCP. The LCP results are further validated with ABAQUS software. To be specific, particular attention has been focused on the LCP mid-span deflection at standstill and the test results are borrowed from the experimental study of Zhang et al. (Zhang et al., 2013).

Experimental Program of Zhang (Zhang et al., 2013)
As part of the Zhang et al. (Zhang et al., 2013) experimental program, a series of 1100 mm long beams were tested under a variable TNT mass and identical standoff distance of 400 mm, as shown in Fig. 25a

Response of the RC Beam Predicted by ABAQUS
The primary aim of this subsection is to validate the RC beam against the full 3D numerical models in finite element commercial software ABAQUS Explicit. Mechanical properties of concrete are modeled using Concrete Damage Plasticity; Table 3. C3D8R brick elements are used to model the concrete part of the beam, and T3D2 wire elements are used to model the steel in the beam shown in Fig. 25(a). The concrete and steel are linked through embedded region interaction. The concrete has a density of 2400kg/m 3 , the Poisson's ratio of 0.2 and the Young modulus of 29.725 GPa. Similarly, the reinforcement has a density of 7850kg/m 3 , the Poisson's ratio of  Page 20 of 25 Khan et al. Int J Concr Struct Mater (2022)   Page 21 of 25 Khan et al. Int J Concr Struct Mater (2022) 16:47 the Young modulus of 10 MPa is provided between the top steel plate and the beam to ensure partial support fixity. The blast load is applied on the beam using built-in module of 'CONWEP' by assigning an equivalent TNT explosive. A convergence study is carried out to optimize the mesh size. The ABAQUS results show a maximum deflection of 8.44 mm for the Beam B2-1, which is very close to the residual experimental deflection of 9.0 mm. The deformed configuration of the beam is shown in Fig. 27. The maximum deflection occurs at the mid-span, conforming to that of the maximum plastic strains in Fig. 28. It is observed from Fig. 28 that the plasticity spread over the central half of the beam. The overall computational time required to solve the beam problem in ABAQUS is about 15 min.

Response of the RC Beam Predicted by LCP Solution
The dynamic response of the beam is predicted using a rigid-plasticity based LCP model; Fig. 23b. The LCP beam is investigated for both 30 contiguous lumped mass elements and a similar number of continuous mass elements, obtained from the beam's mesh convergence tests. The RC beam has a unit mass of 24 kg/m and the plastic moment of resistance is 1710 Nm. Partial fixity of the beam is ensured by assigning one-third of this moment resistance to the supports. The investigated blast loading scenario is shown in Fig. 25b. Using the scaled distance and charge mass, Kinney formulation is employed for computing peak overpressure value (Karlos & Solomon, 2013). Table 4 and Fig. 29 show the mid-span displacements of the RC beam determined from the LCP, ABAQUS, and the experimental study of Zhang (Zhang et al., 2013). Table 4 gives the result of the final central deflection calculated at the respective time. It can be seen from the table that the LCP results are in good agreement with the ABAQUS and experimental results. Fig. 29 shows the evolution of the central displacement, where the LCP results present a small error compared to the ABAQUS and experimental results. The LCP solution showed two traveling plastic hinges within the beam span, which coalesced at the mid-span before the termination of the applied pulse load. At the motion    Khan et al. Int J Concr Struct Mater (2022) 16:47 termination, the central half of the beam becomes plastic, similar to the ABAQUS result of Fig. 28. In the LCP analysis, the initially supplied maximum time step ∆t max to capture the dynamic response is set at 8 × 10 -5 s; Fig. 30. During the recurrent sequence of the Lemke Solver, this increment size is adjusted once to 4.0 × 10 -6 s for ensuring algorithm stability. The overall MATLAB execution time for the frame problem is 10 s, which is substantially less than the ABAQUS run time of 15 min.

Results Discussion and Future Work
The important feature of the numerical implementation developed in this work is its provision for the automatic change of the magnitude of the time increment. This feature is highly advantageous for tracing the true behavior of the concrete and steel structures. Three examples are reported to demonstrate the accuracy and computational efficiency of the LCP formulation, considering the time increment algorithm. Interestingly, the numerical complications identified in the previous study (Khan et al., 2013) of impact problems are nonexistent in the current case of pulse loading. This is because the impact problems required adjustment of the initial velocity profile for starting the LCP formulation; however, the current loading case requires an initial load to initiate the LCP. Moreover, it is found that the Lemke Algorithm (Khan et al., 2021) can break down in the case of blast loading problems if finely discretized continuous mass elements are adopted. Again, this numerical difficulty can be resolved by allowing automatic control of the time increment. However, the error in the results for the continuous element procedure is marginally greater than in the lumped mass discretization, which is contrary to expectation. There is a possibility that the large size of the LCP matrix for the continuous mass model may have contributed to this in the form of round-off errors. The following three examples illustrate the effectiveness of the automatic time control algorithm. The first investigation involves a theoretical study to explore the underlying mechanics of a rigid-plastic propped cantilever beam acted on by a linearly decaying triangular pulse. The nonlinear effects resulting from the asymmetry in the support condition and the decaying pulse load provide a better opportunity to test the improved LCP solver, compared to the relatively simple problem examined recently (Khan et al., 2021). It is seen that the time-step controller subroutine automatically identifies the instant where instability occurs, reduces the time-step size until the instability is resolved, and thereafter continues the evolutionary process. With either the masses being lumped at the element ends or uniformly distributed along the element, the LCP solution traces the dynamic response that matches the derived continuum solution remarkably well.
A second investigation of a pulse-loaded steel portal frame is borrowed from the recent study (Khan et al., 2021) showing that the LCP formulation can effectively capture the complex dynamic response through lumped mass discretization of the frame. In the current study, the portal frame is discretized into continuous mass elements and then the Lemke solver is examined for its computational merit. It is found that mesh refinement using continuous mass elements is a source of algorithmic instability. However, the LCP solver supplemented by the automatic time-stepping algorithm is computationally advantageous. The previous ABAQUS model (Khan et al., 2021) is further improved by the addition of stiffeners at various locations of beam and columns to make the assembly rigid, unlike previously reported analysis (Khan et al., 2021) where the modulus of elasticity was significantly enhanced to ensure the rigidity of the structure. The comparison of LCP results with the ABAQUS solution indicates that the sequence of mechanism movements is difficult to identify in the ABAQUS output because of the interspersed elastic and plastic deformations. In addition, the dynamic response of the ABAQUS frame shows the influence of bending and shear deformations, which is absent in the LCP model since that is based on the idealization of infinite resistance to shear deformation. However, the final deformation mode of the ABAQUS frame model is broadly similar to that of the LCP model. It also turned out that the problems encountered by continuous mass discretization (Khan et al., 2013), such as spurious oscillations in stress-resultants, are eliminated in solving the blast loading problems.
The third benchmark test is the experimental results (Zhang et al., 2013) against which is compared the LCP mid-span displacement for a blast-loaded RC beam with both ends semi-clamped. This experimental result Page 23 of 25 Khan et al. Int J Concr Struct Mater (2022) 16:47 is borrowed from the literature (Zhang et al., 2013). ABAQUS model of the beam is also developed and calibrated on the experimental results. The concrete damage plasticity model is utilized for simulating concrete behavior in ABAQUS. The calibration is carried out until an acceptable level of error is reached in the mid-span displacement. After this calibration, the plastic hinge distribution, evolution of displacement, and strain time history are compared with the LCP solution. It is shown that the LCP solution for both the lumped mass and continuous mass discretization of 50 elements can capture the mid-span deflection satisfactorily. In the case of continuous mass discretization, the time-stepping algorithm again shows exceptional promise in ensuring numerical stability. Finally, it is worth asserting that the proposed automatic time-stepping scheme has made the LCP formulation an efficient computer-based tool, which can be applied to size the skeletal structures subjected to blast loads. This simple finite element tool can predict the main features of behavior for a particular structural design rather quickly for both lumped mass and continuous mass discretization, thereby avoiding the use of advanced finite element analysis methods that impose considerable demand on computational resources. For problems where increased accuracy or further details of the structural behavior is required, then the present formulation may be straightforwardly extended to incorporate additional aspects, such as bending-shear interaction, strain rate, strain hardening, and large displacements. The LCP model involving large displacement of skeletal structures is in progress and the analysis will be reported in due course.

Conclusions
• An automatic time-stepping algorithm has been successfully implemented in the previously developed Lemke solver (Khan et al., 2021) to remedy the wild spurious oscillations in the response variables of pulse-loaded skeletal structures. This improved timestepping strategy allows finer mesh configuration for: (i) lumped mass discretization, and (ii) continuous mass discretization. A series of three examples are presented to demonstrate the accuracy, effectiveness, and versatility of the developed methodology • The first investigation on a propped cantilever beam subjected to a decaying triangular pulse provides evidence that the time-stepping algorithm can stabilize the Lemke solver during various phases of motion. The theoretical response phases are derived herein since, to the authors' knowledge, no closedform solutions are available for this problem. Com-parison of the numerical solution with the theoretical response parameters indicates accurate predictions. For instance, the errors in the central displacement, the bending moment distribution, and the plastic hinge travel are less than 1.5%. • A second investigation of a pulse-loaded steel portal frame is discussed. The numerical difficulties related to the continuous mass discretization are accommodated with the introduction of the automatic time-stepping algorithm. The investigation of 60 elements (lumped mass and continuous mass) LCP frame showed less than 1% error in sway displacement compared to ABAQUS tetrahedron brick element frame. It is also noted that the execution time of the ABAQUS simulation is 20 min, whereas the time for LCP is just 15 s indicating the computational efficiency of the improved LCP solver. • To demonstrate the reliability of the improved LCP solver, the third numerical test is performed on blastloaded reinforced concrete beams. ABAQUS model is developed and validated against the experimental deflection of RC beam subjected to triangular pulse load. Having obtained a good correlation between the ABAQUS model and experimental data, a further comparative study is carried out to test the validity of the LCP model. It is shown by the time history of deflection and strain rates at a control point that the LCP solution for both the lumped mass and continuous mass discretization of 50 elements can capture the response satisfactorily. In the case of LCP continuous mass discretization, the time-stepping algorithm is again able to enhance the numerical stability. • It is asserted that the time-stepping algorithm can regularize the highly oscillatory response of the rigid-plastic skeletal structures when discretized into continuous mass elements.