Comparison between recent implicit time integration methods with frequency dissipation for nonlinear structural applications

Abstract

The present paper aims to test recent (Truly self-starting two sub-step method and three-parameter single-step implicit method) and classical (Generalized-α, HHT-α, and WBZ-α methods) time integration methods using the geometrically nonlinear Positional Finite Element Method (PFEM). The numerical formulation is based on the total Lagrangian approach and uses the Hessian matrix to obtain the response. The mixed hardening inelastic model applied to PFEM is also presented. Two examples validate the time integration algorithms and the inelastic model. In the first example, the mixed hardening inelastic model is compared with the the bilinear stress-strain model and the elastic-perfectly plastic hinge model, and aspects such as amplitude decay and period elongation are discussed. In the second example, the implemented algorithms are verified in a severe geometrically nonlinear example, considering the influence of numerical dissipation, time interval, and the number of elements in the response. Results show the relevance of numerical damping for numerical stabilization and the good performance of the Generalized-α algorithm.

Keywords
Geometrical Nonlinearity; Mixed Hardening Model; Time Integration Methods; Positional Finite Element Formulation; Numerical Dissipation Control; Structural Behavior

Graphical Abstract

1 INTRODUCTION

The response of the equation of motion in structural systems under arbitrary forces can be obtained by modal superposition techniques and methods of direct time integration. The superposition techniques use linear modal analyses and express the response in a series expansion of the vibration modes and modal participation contribution. Its effectiveness is remarkable since the fundamental modes are predominant; however, they are limited to linear cases (Géradin and Rixen, 2015).

In this context, direct integration methods are an alternative, as they allow the consideration of the contribution of high frequencies in the response directly. In a complementary way, they can be applied in linear and nonlinear analyses (Géradin and Rixen, 2015). According to Clough (2003), direct integration methods provide the only complete general approach for the analysis of the nonlinear response and are equally valid in the linear cases, regardless of the structural behavior. These methods are so effective and convenient that time-domain analyses almost always are done by some form of step-by-step analysis, regardless of whether or not the response behavior is linear.

The term “direct” means that no transformation of the equations of motion is performed before the process of numerical integration. According to Bathe (2014), direct numerical integration is based on two ideas. In the first one, the equation of motion need not be satisfied at any time t but at discrete time intervals. The equilibrium equation (including inertial forces and damping) must be satisfied at discrete points in the solution interval. In the second idea, a variation must be assumed for the displacements (or positions), velocities, and accelerations within each time interval. This assumption is determinant in the precision, stability, and computational resource used in the solution procedure.

In the last two decades, researchers proposed improvements in time integration methods or even presented new ones. Chung and Hulbert (1993) presented the family of algorithms called Generalized-α with the numerical dissipation being controlled by the user. The parameter allows the algorithm to dissipate the high-frequencies and minimize the low-frequency dissipation. Besides, the proposed algorithm with the appropriate choice of parameters results in use of the following methods: Newmark, Hilber-Hughes-Taylor-α (HHT-α), and Wood-Bossak-Zienkiewicz-α (WBZ-α).

Kontoe et al. (2004) compared the Generalized-α, Newmark, HHT-α, and WBZ-α time integration schemes for a boundary value problem with a deep foundation subjected to earthquakes. The comparison considered precision, CPU time, and the numerical dissipation control of high frequencies. Their preliminary results demonstrate the superiority of the Generalized-α algorithm in earthquake transient analyses. A simple implicit scheme of time integration for systems under large deformations and long duration of time was proposed by Bathe (2007) for cases in which the Newmark method does not conserve the system energy and is unstable. The method is about twice as expensive as the trapezoidal rule per time step, but less steps can frequently be used, and the method remains stable when the trapezoidal rule fails. Kuo et al. (2012) developed a direct integration algorithm for a nonlinear dynamic analysis of structures whose accuracy is of the fourth-order time and self-starting; this is so that errors caused by the initial estimate of acceleration are eliminated. Still, according to the authors, it is possible to use large time steps. Rostami et al. (2013) presented an algorithm where the cubic B-spline method is developed for systems with multiple degrees of freedom. The algorithm uses periodic cubic base functions over time and is conditionally stable. These authors performed analyses of stability and precision, and concluded that the results of the proposed method are coincident with those from the linear acceleration method; however, they are obtained in a faster and more efficient form since the displacement, velocity, and acceleration vectors are calculated independently.

A new family of self-initializing algorithms for dynamic analysis was proposed by Soares Jr. (2015) with the implicit formulation being unconditionally stable and the explicit formulation having conditioned stability. Moreover, the method has a second-order precision and a parameter for numerical dissipation control. Noh and Bathe (2019a) sought to improve the ρ_∞-Bathe method based on the rate of division of the time step and the spectral radius over long time intervals using the least possible number of parameters. The scheme is effective to prescribe in a smooth manner from no amplitude decay to very large amplitude decays, with a correspondingly small period elongation to very large period elongations while maintaining second-order accuracy.

Malakiyeh et al. (2019) considered the Bathe method in the wave propagation problem. They also pointed out that this method is unconditionally stable and often used without parameter adjustment. However, one can make use of these parameters although detailed numerical experiments are necessary. For Zhang and Xing (2019), time integration methods with three parameters and with a solution step (like the Generalized-α) lose precision when interpolating the external force vectors; this is despite improving numerical dissipation when modifying the equilibrium equation at discrete time points. To solve this problem and still allow the application of these methods in the solution of second-order nonlinear differential equations, the authors propose a new method (also with three parameters and a solution step) called the Three-Parameter Single-Step Implicit Method with an additional variable inserted in the update of equations.

Zhang (2020) tried to carry out a comprehensive study about the Optimal A-Stable Linear Two-Step (OALTS) time integration method. This scheme is implicit with second-order precision for displacement, velocity, and acceleration simultaneously with controllable numerical dissipation. Li and Yu (2020) developed an implicit, self-initializing, unconditionally stable, and second-order precision algorithm for solving nonlinear dynamic problems. The response is obtained through two substeps and uses two parameters to control the numerical dissipation at high and low frequencies.

Borrvall and Lawson (2020) pointed out that the Newmark method is prone to numerical instability when the models are submitted to rotational movements despite its unconditional stability and energy conservation. They used algorithms implemented in LS-DYNA®, such as Bathe, HHT-α, and Finite Rotational Dynamics (FRD), to avoid these instabilities and indicate the most suitable algorithms. Behnoudfar et al. (2020) extended the Generalized-α method whose precision is in the range of second-order to third-order but controlled by a single numerical dissipation parameter. The authors state that the method is unconditionally stable and provide a direction for generalizing the method to higher-order schemes. Kim (2020) reviewed and critically analyzed explicit and implicit numerical integration algorithms, including an explicit algorithm implemented by the author himself. Computational aspects, similarities, and orientations on their applications in the transient analysis were detailed from examples.

Researches have been developed using mixed (or combined) hardening models to consider the behavior of materials based on the linear or nonlinear back-stress approach. Seisenbacher et al. (2018) performed a study of parameters for the combined model of isothermal components mechanically loaded (based on Chaboche's model), which are validated by low cycle fatigue (LCF) tests. The authors state that the model allows simulations independent of the load cycles number and the chosen strain amplitude. In the study by Geng et al. (2018), the mixed hardening model is adopted to analyze the residual stresses of butted weld joint of a huge cylinder with an ultra-thick wall via Lemaitre-Chaboche's formulation. Jiang et al. (2021) carried out a similar investigation with Chaboche's formulation but to determine the fatigue life prediction of the 316L stainless steel weld joint.

Zhao et al. (2019) investigated the low cycle fatigue (LCF) and LCF-creep behaviors of P92 steel with strain-control mode and fully inverse loading waveforms at a specific temperature using the mixed hardening model. Shi et al. (2020) adopted three hardening models to investigate the size effect on metal sheets on the subsequent yield. They identified that the mixed hardening model presented better results compared to the pure isotropic and pure kinematic hardening. The mixed hardening model was also used by Liu et al. (2021) to verify the influence of lubrication conditions in the metal forming process of tubes by rotary swaging. Yang et al. (2020) introduced a methodology to calibrate Chaboche's mixed hardening formulation parameters. This model combines isotropic hardening with kinematic hardening, which is widely used to describe the change in the yield surface. Overall, mixed models are often considered for describing both the expansion of the yield surface and its translation, which, for example, allows the capturing of the Bauschinger effect. Thus, a more accurate characterization of the simulations is guaranteed.

The present work proposes a comparison of numerical time integration algorithms implemented using Positional Finite Element Method (PFEM) formulation, which is based on the total Lagrangian approach and uses the Hessian matrix (instead of stiffness matrix) to obtain the positional response. The formulation’s main advantage in this case is that the chain rule is not explicitly applied to establish the strain gradient (Coda and Paccola, 2008; Coda, 2009; Coda and Paccola, 2010) occurring as a simple numerical matrix inversion. As a result, PFEM can be used to solve any geometrically nonlinear problem (Coda and Paccola, 2010; Coda and Paccola, 2011; Nogueira et al., 2013; Silva and Coda, 2012; Pascon and Coda, 2013). The mixed hardening model is presented in detail and adapted for the proposed formulation to consider the physical nonlinearity. Two highly nonlinear transient examples are used to demonstrate that the PFEM formulation, with the time integration algorithms and the mixed hardening model, can be a significant and viable alternative to geometrical and physical nonlinear analysis of structural plane frames.

2 THE POSITIONAL FINITE ELEMENT METHOD FORMULATION

The PFEM formulation applied to dynamic problems have been developed by Greco (2004), Coda and Greco (2004), and Greco and Coda (2006). Figure 1 presents the generic parameterized curve representing the element, and its geometry is described as a function of the dimensionless variable ξ (which varies from 0 to 1) as well as the global coordinate system.

Fig. 1
Generic parameterized curve representing the plane frame finite element.

The conservation of energy in a mechanical system occurs if the input and output of energy are at balance. If there is some kind of dissipation, the total energy of the system changes along time. The total potential energy of a system can be described as follows:

${\Pi }_0=\Pi +Q(X,t)$ (1)

where Q(X,t) is the quantity of energy withdrawn from the simple conservative idealized energy Π₀, Π is the remaining (actual) mechanical energy of the system, X is the coordinate parameter, and t is the time. For a structural problem associated with a fixed reference system (Figure 1), the ideal potential energy function can be written as the composition of the strain energy U_t, the potential energy of applied forces P, the kinetic energy K, and dissipation Q as follows:

${\Pi }_0={\text{U}}_{\text{t}}+P+\text{K}+Q$ (2)

The strain energy function of the frame U_t is considered stored in the initial volume of the body V₀ and is written as an integral of a specific strain energy value u_e:

$\begin{array}{l}{\text{U}}_{\text{t}}={\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{u}_e}{\text{dV}}_{\text{0}}=\frac{\text{E}}{\text{2}}{\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{\epsilon }^{\text{2}}}{\text{dV}}_{\text{0}}-\text{E}{\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }\epsilon }{\epsilon }^{\text{p}}{\text{dV}}_{\text{0}}\\ \epsilon ={\epsilon }_{\text{med}}+\frac{\text{1}}{\text{r}}\text{z}=\left[\frac{\text{l}}{{\text{l}}_{\text{0}}}\sqrt{{\left({\text{l}}_{\text{x}}\right)}^{\text{2}}+{\left({\text{3c}}_{\text{1}}{\xi }^{\text{2}}+{\text{2c}}_{\text{2}}\xi +{\text{c}}_{\text{3}}\right)}^{\text{2}}}-1\right]+\frac{\text{1}}{\text{r}}\text{z}\\ \frac{\text{1}}{\text{r}}=\frac{{\text{l}}_{\text{x}}\left({\text{6c}}_{\text{2}}\xi +{\text{2c}}_{\text{3}}\right)}{{\left(\sqrt{{\left({\text{l}}_{\text{x}}\right)}^{\text{2}}+{\left({\text{3c}}_{\text{1}}{\xi }^{\text{2}}+{\text{2c}}_{\text{2}}\xi +{\text{c}}_{\text{3}}\right)}^{\text{2}}}\right)}^{\text{3}}}\end{array}$ (3)

where c₁, c₂ and c₃ are functions of θ₁, θ₂, l_x(=X₂-X₁) and l_y(=Y₂-Y₁); 1/r is the curvature, E is the Young modulus, r is the curvature radius, z is the height of equivalent cross-section, ε_med is the engineering strain, ε is the strain measure for geometric nonlinearity, ε^p is the plastic strain, and ξ is a dimensionless variable (from 0 to 1).

In the initial position (non-deformed), the strain energy is assumed to be zero. The work of the applied conservative concentrated forces P is written as follows:

$P=-F_i{X}_i$ (4)

where F_i represents forces (or moments) applied in i direction and X_i is the ith coordinate parameter of the point where the load is applied. The kinetic energy is given by the following:

$$\text{K}=\frac{1}{2}{\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{\rho }_0}{\stackrel{.}{\text{X}}}_{\text{i}}{\stackrel{.}{\text{X}}}_{\text{i}}{\text{dV}}_{\text{0}}$$ (5)

where ${\stackrel{.}{\text{X}}}_{\text{i}}$ is the ith velocity and ρ₀ is the specific mass. The dissipative term is written in its differential form as follows:

$$\frac{\partial }{\partial {\text{X}}_{\text{i}}}\text{Q}(X,t)={\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }\frac{\partial }{\partial {\text{X}}_{\text{i}}}q(\text{X},t)}{\text{dV}}_{\text{0}}={\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{\lambda }_m{\stackrel{.}{\text{X}}}_i}{\text{dV}}_{\text{0}}$$ (6)

where q(X,t) is the specific dissipative functional, λ_m is a proportionality constant, and X_i is the position of any specific point (for PFEM, it is a nodal position). The potential energy of the applied forces may not be zero in the reference configuration (Coda and Greco, 2004). Eq. (2) can be rewritten as follows:

$${\Pi }_0={\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{u}_e}{\text{dV}}_{\text{0}}-F_i{X}_i+\frac{1}{2}{\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{\rho }_0}{\stackrel{.}{\text{X}}}_{\text{i}}{\stackrel{.}{\text{X}}}_{\text{i}}{\text{dV}}_{\text{0}}+Q$$ (7)

The energy function can be evaluated by the following approximation:

${\Pi }_0={\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{u}_e}\text{(}\xi {\text{,X}}_{\text{i}}{\text{)dV}}_{\text{0}}-F_i{X}_i+\frac{1}{2}{\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{\rho }_0}{\stackrel{.}{\text{X}}}_i^2\text{(}\xi {\text{,X}}_{\text{i}}{\text{)dV}}_{\text{0}}+Q\text{(}\xi {\text{,X}}_{\text{i}}\text{)}$ (8)

The principle of minimum potential energy consists of the real values ​​that satisfy the equilibrium equations and make the potential energy stationary for all admissible displacements in a certain system. This principle is used in _Π0 by differentiating Eq. (8) regarding a generic nodal position X_s, which results in the following:

$$\begin{array}{l}\frac{\partial {\Pi }_0}{\partial X_s}={\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }\frac{u_e\text{(}\xi {\text{,X}}_{\text{i}}\text{)}}{\partial X_s}}{\text{dV}}_{\text{0}}-F_s+{\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{\rho }_0}{\stackrel{.}{X}}_i\text{(}\xi {\text{,X}}_{\text{i}}\text{)}\frac{\partial {\stackrel{.}{X}}_i\text{(}\xi {\text{,X}}_{\text{i}}\text{)}}{\partial X_s}{\text{dV}}_{\text{0}}\\ +{\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{\lambda }_m{\rho }_0{\stackrel{.}{\text{X}}}_s}\text{(}\xi {\text{,X}}_{\text{i}}{\text{)dV}}_{\text{0}}=0\end{array}$$ (9)

In Eq. (9), ${\stackrel{.}{\text{X}}}_{\text{s}}$ is a generic nodal velocity. In a vector form, Eq. (9) can be rewritten as follows:

$\frac{\partial {\Pi }_0}{\partial X}=\frac{\partial {\text{U}}_{\text{t}}}{\partial X}-F+F_{inert.}+F_{damp.}=0$ (10)

where the internal forces vector ∂U_t/∂X, the inertia forces F_inert., and the damping forces F_damp. are given, respectively, by the following:

$\frac{\partial {\text{U}}_{\text{t}}}{\partial X}={\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }\frac{u_e\text{(}\xi {\text{,X}}_{\text{i}}\text{)}}{\partial X_s}}{\text{dV}}_{\text{0}}$ (11)

$F_{inert.}={\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{\rho }_0}{\dot{X}}_i\text{(}\xi {\text{,X}}_{\text{i}}\text{)}\frac{\partial {\stackrel{.}{X}}_i\text{(}\xi {\text{,X}}_{\text{i}}\text{)}}{\partial X_s}{\text{dV}}_{\text{0}}=M\ddot{X}$ (12)

$F_{damp.}={\displaystyle \underset{{\text{V}}_{\text{0}}}{\int }{\lambda }_m{\rho }_0{\stackrel{.}{\text{X}}}_s}\text{(}\xi {\text{,X}}_{\text{i}}{\text{)dV}}_{\text{0}}=C\dot{X}$ (13)

In Eq. (12)-(13), M, C, $\dot{X}$_, and $\ddot{X}$ are mass matrix, damping matrix, velocity vector, and acceleration vector, respectively. Then, Eq. (10) becomes:

$\frac{\partial {\Pi }_0}{\partial X}=\frac{\partial {\text{U}}_{\text{t}}}{\partial X}-F+M\ddot{X}+C\dot{X}=0$ (14)

In Eq. (14), time t and position vector X are the variables, and F is the external loads vector. It is necessary to discretize the system to perform the numerical integration in the time domain. Rewriting for time t+Δt:

${\left.\frac{\partial {\Pi }_0}{\partial X}\right|}_{\text{t}+\text{Δt}}={\left.\frac{\partial {\text{U}}_{\text{t}}}{\partial X}\right|}_{\text{t}+\text{Δt}}-F_{\text{t}+\text{Δt}}+M{\ddot{X}}_{\text{t}+\text{Δt}}+C{\stackrel{.}{X}}_{\text{t}+\text{Δt}}=0$ (15)

which represents the condition of geometric nonlinear dynamic equilibrium equation of motion, whose solution can be obtained using an iterative process.

3 NONLINEAR POSITIONAL FINITE ELEMENT FORMULATION FOR SOME TIME INTEGRATION METHODS

The following algorithm schemes use the spectral radius ρ_∞ [0; 1] as a numerical dissipation control parameter for the contribution of high frequencies to the response. In the lower limit, there is a complete dissipation of the high-frequency contribution response. On the other hand, there is no dissipation in the upper limit.

3.1 The Generalized-α, HHT-α, and WBZ-α methods

These methods have the following characteristics: implicit, second-order time precision, unconditional stability for the linear case, and numerical dissipation control of high frequencies (Chung and Hulbert, 1993). Eq. (15) can be rewritten as follows:

${\left.\frac{\partial {\Pi }_0}{\partial X}\right|}_{\text{t}+\text{Δt}}={\left.\frac{\partial {\text{U}}_{\text{t}}}{\partial X}\right|}_{\text{t}+\text{Δt}-{\alpha }_{\text{f}}}-F_{\text{t}+\text{Δt}-{\alpha }_{\text{f}}}+M{\ddot{X}}_{\text{t}+\text{Δt}-{\alpha }_{\text{m}}}+C{\stackrel{.}{X}}_{\text{t}+\text{Δt}-{\alpha }_{\text{f}}}=0$ (16)

The position, velocity, acceleration, and external force vectors are given, respectively, by the following:

$X_{\text{t}+\text{Δt}-{\alpha }_{\text{f}}}=\left(\text{1}-{\alpha }_{\text{f}}\right)X_{\text{t}+\Delta \text{t}}+{\alpha }_{\text{f}}{X}_{\text{t}}$ (17)

${\dot{X}}_{\text{t}+\text{Δt}-{\alpha }_{\text{f}}}=\left(\text{1}-{\alpha }_{\text{f}}\right){\stackrel{.}{X}}_{\text{t}+\Delta \text{t}}+{\alpha }_{\text{f}}{\stackrel{.}{X}}_{\text{t}}$ (18)

${\ddot{X}}_{\text{t}+\text{Δt}-{\alpha }_{\text{m}}}=\left(\text{1}-{\alpha }_{\text{m}}\right){\ddot{X}}_{\text{t}+\Delta \text{t}}+{\alpha }_{\text{m}}{\ddot{X}}_{\text{t}}$ (19)

$F_{\text{t}+\Delta \text{t}-{\alpha }_{\text{f}}}=\left(\text{1}-{\alpha }_{\text{f}}\right)F_{\text{t}+\Delta \text{t}}+{\alpha }_{\text{f}}{F}_{\text{t}}$ (20)

with acceleration and velocity at t+Δt, respectively, equal to:

${\ddot{X}}_{\text{t}+\text{Δt}}=-\frac{\text{1}}{\beta \Delta \text{t}}{\dot{X}}_{\text{t}}-\frac{\text{1}}{\beta \Delta {\text{t}}^{\text{2}}}{X}_{\text{t}}+\frac{\text{1}}{\beta \Delta {\text{t}}^{\text{2}}}{X}_{\text{t}+\text{Δt}}-\frac{\text{1}}{\text{2}\beta }{\ddot{X}}_{\text{t}}+{\ddot{X}}_{\text{t}}$ (21)

${\dot{X}}_{\text{t}+\text{Δt}}=-\frac{\gamma }{\beta }{\stackrel{.}{X}}_{\text{t}}+{\stackrel{.}{X}}_{\text{t}}-\frac{\gamma }{\beta \Delta \text{t}}{X}_{\text{t}}-\frac{\gamma \Delta \text{t}}{\text{2}\beta }{\ddot{X}}_{\text{t}}+\frac{\gamma }{\beta \Delta \text{t}}{X}_{\text{t}+\Delta \text{t}}+\Delta \text{t}{\ddot{X}}_{\text{t}}$ (22)

The equation of motion can be described as follows:

$\begin{array}{c}{\left.\frac{\partial {\Pi }_0}{\partial X}\right|}_{\text{t}+\text{Δt}}=\left({\alpha }_{\text{f}}-\text{1}\right){\left.\frac{\partial {\text{U}}_{\text{t}}}{\partial X}\right|}_{\text{t}+\text{Δt}}-\left(\left(\text{1}-{\alpha }_{\text{f}}\right)F_{\text{t}+\Delta \text{t}}+{\alpha }_{\text{f}}{F}_{\text{t}}\right)+\frac{\left({\alpha }_{\text{m}}-\text{1}\right)}{\beta \Delta {\text{t}}^{\text{2}}}M{X}_{\text{t}+\text{Δt}}\\ +M{D}_{\text{M}}+C{D}_{\text{C}}+\frac{\gamma \left({\alpha }_{\text{f}}-\text{1}\right)}{\beta \Delta \text{t}}C{X}_{\text{t}+\text{Δt}}=0\end{array}$ (23)

$D_{\text{M}}={\ddot{X}}_{\text{t}}+\frac{\left({\alpha }_{\text{m}}-\text{1}\right)}{\beta \Delta {\text{t}}^{\text{2}}}{X}_{\text{t}}+\frac{\left({\alpha }_{\text{m}}-\text{1}\right)}{\beta \Delta \text{t}}{\stackrel{.}{X}}_{\text{t}}+\frac{\left({\alpha }_{\text{m}}-\text{1}\right)}{\text{2}\beta }{\ddot{X}}_{\text{t}}$ (24)

$$D_{\text{C}}={\stackrel{.}{X}}_{\text{t}}+\frac{\left({\alpha }_{\text{f}}-\text{1}\right)\gamma }{\beta }{\stackrel{.}{X}}_{\text{t}}+\frac{\left({\alpha }_{\text{f}}-\text{1}\right)\Delta \text{t}\gamma }{\text{2}\beta }{\ddot{X}}_{\text{t}}+\left(\text{1}-{\alpha }_{\text{f}}\right)\Delta \text{t}{\ddot{X}}_{\text{t}}$$ (25)

where D_M and D_C are the vectors related to the mass and damping matrices, respectively, with the known variables at time t. The Hessian matrix is given by the following:

${\left.\frac{{\partial }^{\text{2}}{\Pi }_0}{\partial X^{\text{2}}}\right|}_{\text{t}+\text{Δt}}=\left({\alpha }_{\text{f}}-\text{1}\right){\left.\frac{{\partial }^{\text{2}}{\text{U}}_{\text{t}}}{\partial X^{\text{2}}}\right|}_{\text{t}+\text{Δt}}+\frac{\left({\alpha }_{\text{m}}-\text{1}\right)}{\beta \Delta {\text{t}}^{\text{2}}}M+\frac{\left({\alpha }_{\text{f}}-\text{1}\right)\gamma }{\beta \Delta \text{t}}C$ (26)

with the following parameters:

$\beta =\frac{{\left(\text{1}-{\alpha }_{\text{m}}+{\alpha }_{\text{f}}\right)}^{\text{2}}}{\text{4}};\gamma =\frac{\text{1}}{\text{2}}-{\alpha }_{\text{m}}+{\alpha }_{\text{f}};{\alpha }_{\text{m}}=\frac{\text{2}{\rho }_{\infty }-\text{1}}{{\rho }_{\infty }+\text{1}};{\alpha }_{\text{f}}=\frac{{\rho }_{\infty }}{{\rho }_{\infty }+\text{1}}$ (27)

Eq. (27) presents the algorithm parameters: α_m corresponds to the mass terms, and α_f is related to the other terms. In addition, β and γ are complementary parameters of the method. The HHT-α and WBZ-α methods, which can be found in the research of Hilber et al. (1977) and Wood et al. (1980), are obtained considering α_m=0 and α_f=0, respectively.

3.2 Truly Self-Starting Two Sub-Step method

The method proposed by Li and Yu (2020) has significant features. Acceleration is not necessary to initialize the algorithm (truly self-start), but it can be achieved at each time step. Besides, the displacement (or position), velocity, and acceleration are obtained with second-order precision. The method has unconditional stability for linear elastic problems.

3.2.1 First sub-step

The first variation of the total potential energy, Eq. (14), can be approximated by the following:

${\left.\frac{\partial {\Pi }_0}{\partial X}\right|}_{\text{t}+\gamma \text{Δt}}={\left.\frac{\partial {\text{U}}_{\text{t}}}{\partial X}\right|}_{\text{t}+\gamma \text{Δt}}-F_{\text{t}+\gamma \text{Δt}}+M{\ddot{X}}_{\text{t}+\gamma \text{Δt}}+C{\stackrel{.}{X}}_{\text{t}+\gamma \text{Δt}}=0$ (28)

and position, velocity and external force vectors can be written, respectively, as follows:

$X_{\text{t}+\gamma \text{Δt}}=X_{\text{t}}+\Delta \text{t}\gamma {\stackrel{.}{X}}_{\text{t}+\gamma \text{Δt}}$ (29)

${\dot{X}}_{\text{t}+\gamma \text{Δt}}={\stackrel{.}{X}}_{\text{t}}+\Delta \text{t}\gamma {\ddot{X}}_{\text{t}+\gamma \text{Δt}}$ (30)

$F_{\text{t}+\gamma \text{Δt}}=\left(\text{1}-\gamma \right)F_{\text{t}+\gamma \text{Δt}}+\gamma F_{\text{t}}$ (31)

Thus, the equation of motion is obtained by the following:

$\begin{array}{c}{\left.\frac{\partial {\Pi }_0}{\partial X}\right|}_{\text{t}+\gamma \text{Δt}}={\left.\frac{\partial {\text{U}}_{\text{t}}}{\partial X}\right|}_{\text{t}+\gamma \text{Δt}}-\left[\left(\text{1}-\gamma \right)F_{\text{t}+\gamma \text{Δt}}+\gamma F_{\text{t}}\right]+\frac{\text{1}}{\gamma \Delta {\text{t}}^{\text{2}}}M{X}_{\text{t}+\gamma \text{Δt}}\\ +M{D}_{\text{M}}+C{D}_{\text{C}}+\frac{\text{1}}{\gamma \Delta \text{t}}C{X}_{\text{t}+\gamma \text{Δt}}=0\end{array}$ (32)

$D_{\text{M}}=-\frac{\text{1}}{\gamma \Delta \text{t}}{\stackrel{.}{X}}_{\text{t}}-\frac{\text{1}}{{\left(\gamma \Delta \text{t}\right)}^{\text{2}}}{X}_{\text{t}}$ (33)

$D_{\text{C}}=\frac{\text{1}}{\gamma \Delta \text{t}}{X}_{\text{t}}$ (34)

and the second variation of the total potential energy leads to the Hessian Matrix:

${\left.\frac{{\partial }^{\text{2}}{\Pi }_0}{\partial X^{\text{2}}}\right|}_{\text{t}+\gamma \Delta \text{t}}={\left.\frac{{\partial }^{\text{2}}{\text{U}}_{\text{t}}}{\partial X^{\text{2}}}\right|}_{\text{t}+\gamma \Delta \text{t}}+\frac{\text{1}}{{\left(\beta \Delta \text{t}\right)}^{\text{2}}}M+\frac{\text{1}}{\gamma \Delta \text{t}}C$ (35)

3.2.2 Second sub-step

Similarly, the Eqs. (28)-(35) for the second sub-step are given by the following:

${\left.\frac{\partial {\Pi }_0}{\partial X}\right|}_{\text{t}+{\gamma }_1\Delta \text{t}}={\left.\frac{\partial {\text{U}}_{\text{t}}}{\partial X}\right|}_{\text{t}+{\gamma }_1\Delta \text{t}}-F_{\text{t}+{\gamma }_{\text{1}}\text{Δt}}+M{\ddot{X}}_{\text{t}+{\gamma }_1\Delta \text{t}}+C{\stackrel{.}{X}}_{\text{t}+{\gamma }_1\Delta \text{t}}=0$ (36)

$X_{\text{t}+{\gamma }_1\Delta \text{t}}=X_{\text{t}}+\Delta \text{t}\left[\left({\gamma }_{\text{1}}-\alpha \right){\stackrel{.}{X}}_{\text{t}+\gamma \text{Δt}}+\alpha {\stackrel{.}{X}}_{\text{t}+{\gamma }_{\text{1}}\text{Δt}}\right]$ (37)

${\dot{X}}_{\text{t}+{\gamma }_1\Delta \text{t}}={\dot{X}}_{\text{t}}+\Delta \text{t}\left[\left({\gamma }_{\text{1}}-\alpha \right){\ddot{X}}_{\text{t}+\gamma \text{Δt}}+\alpha {\ddot{X}}_{\text{t}+{\gamma }_1\Delta \text{t}}\right]$ (38)

$F_{\text{t}+{\gamma }_1\Delta \text{t}}=\left(\text{1}-{\gamma }_{\text{1}}\right)F_{\text{t}+{\gamma }_1\Delta \text{t}}+{\gamma }_{\text{1}}{F}_{\text{t}}$ (39)

$\begin{array}{c}{\left.\frac{\partial {\Pi }_0}{\partial X}\right|}_{\text{t}+{\gamma }_1\Delta \text{t}}={\left.\frac{\partial {\text{U}}_{\text{t}}}{\partial X}\right|}_{\text{t}+{\gamma }_1\Delta \text{t}}-\left[\left(\text{1}-{\gamma }_{\text{1}}\right)F_{\text{t}+{\gamma }_1\Delta \text{t}}+{\gamma }_{\text{1}}{F}_{\text{t}}\right]+\frac{\text{1}}{{\left(\alpha \Delta \text{t}\right)}^{\text{2}}}M{X}_{\text{t}+{\gamma }_1\Delta \text{t}}\\ +M{D}_{\text{M}}+C{D}_{\text{C}}+\frac{\text{1}}{\alpha \Delta \text{t}}C{X}_{\text{t}+{\gamma }_1\Delta \text{t}}=0\end{array}$ (40)

$\begin{array}{c}{D}_{\text{M}}=-\frac{\text{1}}{\alpha \Delta \text{t}}{\dot{X}}_{\text{t}}-\frac{{\gamma }_{\text{1}}}{\alpha }{\ddot{X}}_{\text{t}+\gamma \Delta \text{t}}+{\ddot{X}}_{\text{t}+\gamma \Delta \text{t}}-\frac{{\gamma }_{\text{1}}}{{\alpha }^{\text{2}}\Delta \text{t}}{\stackrel{.}{X}}_{\text{t}+\gamma \Delta \text{t}}\\ +\frac{\text{1}}{\alpha \Delta \text{t}}{\stackrel{.}{X}}_{\text{t}+\gamma \Delta \text{t}}-\frac{\text{1}}{{\left(\alpha \Delta \text{t}\right)}^{\text{2}}}{X}_{\text{t}}\end{array}$ (41)

$D_{\text{C}}=-\frac{{\gamma }_{\text{1}}}{\alpha }{\stackrel{.}{X}}_{\text{t}+\gamma \Delta \text{t}}+{\stackrel{.}{X}}_{\text{t}+\gamma \Delta \text{t}}-\frac{\text{1}}{\alpha \Delta \text{t}}{X}_{\text{t}}$ (42)

${\left.\frac{{\partial }^{\text{2}}{\Pi }_0}{\partial X^{\text{2}}}\right|}_{\text{t}+\gamma \Delta \text{t}}={\left.\frac{{\partial }^{\text{2}}{\text{U}}_{\text{t}}}{\partial X^{\text{2}}}\right|}_{\text{t}+\gamma \Delta \text{t}}+\frac{\text{1}}{{\left(\beta \Delta \text{t}\right)}^{\text{2}}}M+\frac{\text{1}}{\gamma \Delta \text{t}}C$ (43)

3.2.3 Updating the variables

After the sub-steps, the updated vectors for position and velocity are given, respectively, by the following:

$X_{\text{t}+\text{Δt}}=X_{\text{t}}+\Delta \text{t}\left[\left(\text{1}-\beta \right){\stackrel{.}{X}}_{\text{t}+\gamma \text{Δt}}+\beta {\stackrel{.}{X}}_{\text{t}+{\gamma }_1\Delta \text{t}}\right]$ (44)

${\stackrel{.}{X}}_{\text{t}+\text{Δt}}={\stackrel{.}{X}}_{\text{t}}+\Delta \text{t}\left[\left(\text{1}-\beta \right){\ddot{X}}_{\text{t}+\gamma \text{Δt}}+\beta {\ddot{X}}_{\text{t}+{\gamma }_1\Delta \text{t}}\right]$ (45)

and the acceleration equation assumes the form:

${\ddot{X}}_{\text{t}+\text{Δt}}=\left(\text{1}-\chi \right){\ddot{X}}_{\text{t}+\gamma \text{Δt}}+\chi {\ddot{X}}_{\text{t}+{\gamma }_1\Delta \text{t}}$ (46)

The parameters of the method are presented in Eq. (47):

$\alpha =\frac{\text{1}-\text{2}\gamma }{\text{2}\left({\rho }_{\infty }\gamma -\gamma +\text{1}\right)};{\gamma }_{\text{1}}=\alpha +\gamma +\sqrt{{\alpha }^{\text{2}}+{\gamma }^{\text{2}}};\chi =\frac{\gamma -\text{1}}{\gamma -{\gamma }_{\text{1}}};\beta =\frac{\text{2}\gamma -\text{1}}{\text{2}\left(\gamma -{\gamma }_{\text{1}}\right)}$ (47)

where γ and ρ_∞ are the coefficients to control numerical dissipation at low and high frequencies, respectively. The dependent parameters β, χ, α and γ₁ are used to ensure accuracy requirements of the present scheme.

3.3 Three-Parameter Single-Step Implicit Method

The method proposed by Zhang and Xing (2019), called the Three-Parameter Single-Step method, with an implicit form strictly satisfies the equilibrium equation in each time step and has a desirable performance according to the authors. The position, velocity, acceleration, and the arbitrary vector θ_t+Δt are defined as follows:

$X_{\text{t}+\text{Δt}}=X_{\text{t}}+\text{Δt}{\stackrel{.}{X}}_{\text{t}}+\frac{{\text{Δt}}^{\text{2}}}{\text{2}}{\ddot{X}}_{\text{t}}+\frac{{\text{Δt}}^{\text{3}}}{\text{6}}\left[\left(\text{1}-\alpha \right){\theta }_{\text{t}}+\alpha {\theta }_{\text{t}+\text{Δt}}\right]$ (48)

${\stackrel{.}{X}}_{\text{t}+\text{Δt}}={\stackrel{.}{X}}_{\text{t}}+\text{Δt}{\ddot{X}}_{\text{t}}+\frac{{\text{Δt}}^{\text{2}}}{\text{2}}\left[\left(\text{1}-\delta \right){\theta }_{\text{t}}+\delta {\theta }_{\text{t}+\text{Δt}}\right]$ (49)

${\ddot{X}}_{\text{t}+\text{Δt}}={\ddot{X}}_{\text{t}}+\text{Δt}\left[\left(\text{1}-\gamma \right){\theta }_{\text{t}}+\gamma {\theta }_{\text{t}+\text{Δt}}\right]$ (50)

${\theta }_{\text{t}+\text{Δt}}=-\frac{\text{6}}{\alpha {\text{Δt}}^{\text{2}}}{\stackrel{.}{X}}_{\text{t}}-\frac{{\theta }_{\text{t}}}{\alpha }+{\theta }_{\text{t}}-\frac{\text{6}}{\alpha {\text{Δt}}^{\text{3}}}{X}_{\text{t}}-\frac{\text{3}}{\alpha \text{Δt}}{\ddot{X}}_{\text{t}}+\frac{\text{6}}{\alpha {\text{Δt}}^{\text{3}}}{X}_{\text{t}+\text{Δt}}$ (51)

The arbitrary vector can be defined as θ_t=0 for t=0 and must be updated at each step. The equation of motion is defined to be:

$\begin{array}{l}{\left.\frac{\partial {\Pi }_0}{\partial X}\right|}_{\text{t}+\text{Δt}}={\left.\frac{\partial {\text{U}}_{\text{t}}}{\partial X}\right|}_{\text{t}+\text{Δt}}-F_{\text{t}+\Delta \text{t}}+\frac{\text{6}\gamma }{\alpha \Delta {\text{t}}^{\text{2}}}M{X}_{\text{t}+\text{Δt}}\\ +M{D}_{\text{M}}+C{D}_{\text{C}}+\frac{\text{3}\delta }{\alpha \Delta \text{t}}C{X}_{\text{t}+\text{Δt}}=0\end{array}$ (52)

$D_{\text{M}}=\frac{\text{6}\gamma }{\alpha \text{Δt}}{\dot{X}}_{\text{t}}+\frac{\text{Δt}\gamma }{\alpha }{\theta }_{\text{t}}-\text{Δt}{\theta }_{\text{t}}+\frac{\text{6}\gamma }{\alpha {\text{Δt}}^{\text{2}}}{X}_{\text{t}}+\frac{\text{3}\gamma }{\alpha }{\ddot{X}}_{\text{t}}-{\ddot{X}}_{\text{t}}$ (53)

$D_{\text{C}}=\frac{\text{3}\delta }{\alpha }{\dot{X}}_{\text{t}}-{\dot{X}}_{\text{t}}+\frac{\delta {\text{Δt}}^{\text{2}}}{\text{2}\alpha }{\theta }_{\text{t}}-\frac{{\text{Δt}}^{\text{2}}}{\text{2}}{\theta }_{\text{t}}+\frac{\text{3}\delta }{\alpha \text{Δt}}{X}_{\text{t}}+\frac{\text{3}\delta \text{Δt}}{\text{2}\alpha }{\ddot{X}}_{\text{t}}-\text{Δt}{\ddot{X}}_{\text{t}}$ (54)

The Hessian Matrix is such that:

${\left.\frac{{\partial }^{\text{2}}{\Pi }_0}{\partial X^{\text{2}}}\right|}_{\text{t}+\text{Δt}}={\left.\frac{{\partial }^{\text{2}}{\text{U}}_{\text{t}}}{\partial X^{\text{2}}}\right|}_{\text{t}+\text{Δt}}+\frac{\text{6}\gamma }{\alpha \Delta {\text{t}}^{\text{2}}}M+\frac{\text{3}\gamma }{\alpha \Delta \text{t}}C$ (55)

and the parameters are obtained by:

$\alpha =\frac{\text{6}}{{\left({\rho }_{\infty }+\text{1}\right)}^{\text{3}}};\delta =\frac{\text{2}{\rho }_{\infty }^{\text{2}}-\text{5}{\rho }_{\infty }+\text{11}}{\text{3}{\left({\rho }_{\infty }+\text{1}\right)}^{\text{2}}};\gamma =\frac{\text{2}-{\rho }_{\infty }}{{\rho }_{\infty }+\text{1}}$ (56)

The spectral radius ρ_∞ is the only independent parameter, with δ, α, and γ as a dependent parameters to optimize the implicit method.

3.4 Research Methodology Flowchart

The following flowchart algorithm (Figure 2) shows the research methodology of the present work. The standard Newton-Raphson is used as an iterative process and the Euclidean norm as a convergence criterion.

Fig. 2
Flowchart of the algorithm proposed in this work.

4 MIXED HARDENING MODEL APPLIED TO PFEM

A mixed hardening model combines isotropic hardening (expansion of the elastic surface) and kinematic hardening (translation of the elastic domain center). It is simple and has proximity to the experimental occurrence, particularly for cases of plasticity in metals. Greco (2004) presented the implementation of this inelastic model in the PFEM formulation. Figure 3 shows a graphical representation for a one-dimensional element.

Fig. 3
Evolution of the yield surface and yield stress.

The yield function f and yield criterion used are presented, respectively, by the following:

$\text{f}=\left(\xi \right)-\left({\sigma }_{\text{y}}+\text{K}\alpha \right)\text{,}\xi =\sigma -\text{q}$ (57)

$\text{f}\le \text{0}$ (58)

In Eq. (57), K is the plastic module. The measures α and q (back-stress) are internal variables that measure the expansion and change of the center of the elastic domain, respectively. The yield stress of the material is σ_Y. A state is said to be admissible when Eq. (58) is verified.

By assumption, considering a strain decomposition, the stress is given by the following:

$\sigma =\text{E}\left(\epsilon -{\epsilon }^{\text{p}}\right)$ (59)

The evolution of the plastic strain (or plastic flow), which is described by Eq. (60), and of the internal hardening variable, Eq. (61), are given by the following:

${\dot{\epsilon }}^{\text{p}}=\gamma \text{sign}\left(\xi \right)\text{,}\gamma \ge \text{0}$ (60)

$\dot{\alpha }=\left|{\dot{\epsilon }}^p\right|=\gamma$ (61)

where γ is called plastic multiplier. The evolution of the internal variable $\dot{q}$ is given by Ziegler's Law:

$\dot{q}=H{\dot{\epsilon }}^p=\gamma H\text{sign}\left(\xi \right)$ (62)

and H is the kinematic hardening module (constant).

Two conditions govern the admissibility of stress states and the evolution of plastic deformation. The first one, called the Kuhn-Tucker condition, is given by the relationships in Eq. (63):

$\gamma \ge \text{0}\text{,}\text{f}\left(\sigma \text{,q}\text{,}\alpha \right)\le \text{0}\text{,}\gamma \text{f}\left(\sigma \text{,q}\text{,}\alpha \right)=\text{0}$ (63)

which guarantees an admissible state. In this case, a plastic deformation occurs only when $\text{f}\left(\sigma \text{,q}\text{,}\alpha \right)=\text{0}$, i.e.,$\gamma >\text{0}$.

The second one, called consistency condition, is described as follows:

$\gamma \stackrel{.}{\text{f}}\left(\sigma \text{,q}\text{,}\alpha \right)=\text{0}\text{,}\text{if}\stackrel{.}{\text{f}}\left(\sigma \text{,q}\text{,}\alpha \right)=\text{0}$ and $\gamma >\text{0}$ (64)

and requires for ${\dot{\epsilon }}^p$ that $\gamma >\text{0}$. It must persist on the yield surface so that $\stackrel{.}{\text{f}}\left(\sigma \text{,q}\text{,}\alpha \right)=\text{0}$.

Based on the two conditions, Eq. (63) and (64), and the yield criterion function, Eq. (57), it is possible to write:

$\gamma =\frac{\text{sign}\left(\xi \right)\text{E}\dot{\epsilon }}{\text{E}+\text{H}+\text{K}}$ (65)

When K=H=0, the model becomes perfectly plastic. For K>0 and H=0, pure isotropic hardening is obtained.

4.1 Solution based on incremental terms of the variables

The problem solution involves the consistency condition and the Kuhn-Tucker condition. The insertion of a trial stress state (supposedly elastic) allows writing the solution in incremental terms:

${\text{f}}_{\text{n}+\text{1}}^{\text{trial}}=\left({\xi }_{\text{n}+\text{1}}^{\text{trial}}\right)-\left[{\sigma }_{\text{Y}}+\text{K}{\alpha }_n\right]$, ${\xi }_{\text{n}+\text{1}}^{\text{trial}}={\sigma }_{\text{n}+\text{1}}^{\text{trial}}-{\text{q}}_{\text{n}}$(66)

$\Delta \gamma =\frac{{\text{f}}_{\text{n}+\text{1}}^{\text{trial}}}{\left(\text{E}+\text{K}+\text{H}\right)}>\text{0}$ (67)

${\sigma }_{\text{n}+\text{1}}={\sigma }_{\text{n}+\text{1}}^{\text{trial}}-\Delta \gamma \text{Esign}\left({\xi }_{\text{n}+\text{1}}^{\text{trial}}\right)$ (68)

${\epsilon }_{\text{n}+\text{1}}^{\text{p}}={\epsilon }_{\text{n}}^{\text{p}}+\Delta \gamma \text{sign}\left({\xi }_{\text{n}+\text{1}}^{\text{trial}}\right)$ (69)

${\text{q}}_{\text{n}+\text{1}}={\text{q}}_{\text{n}}+\Delta \gamma \text{Hsign}\left({\xi }_{\text{n}+\text{1}}^{\text{trial}}\right)$ (70)

${\alpha }_{\text{n}+\text{1}}={\alpha }_{\text{n}}+\Delta \gamma$ (71)

with the initial conditions ${\epsilon }_{\text{n}}^{\text{p}}$, ${\alpha }_{\text{n}}$, and ${\text{q}}_{\text{n}}$ to solve the differential equation problem.

The following relations are based on an incremental trial state and define whether plastic correction is necessary as is shown below in Eq. (72):

${\text{f}}_{\text{n}+\text{1}}^{\text{trial}}=\left\{\begin{array}{c}\le \text{0}\text{,}\Delta \gamma =\text{0}\text{then}\text{elastic}\text{step}\\ >\text{0}\text{,}\Delta \gamma >\text{0}\text{then}\text{plastic}\text{step}\end{array}\right\}$ (72)

Chart 1, which is based on Eqs. (66)-(72), shows the step-by-step solution for the incremental plastic correction. For details regarding the algorithm, see Greco (2004) and Simo and Hughes (2006). This incremental plastic correction is updated at each time step.

Chart 1
Step-by-step algorithm for mixed hardening model.

1. Given the initial variables ${\epsilon }_{\text{n}}^{\text{p}}=\text{0}$, ${\alpha }_{\text{n}}=\text{0}$, ${\text{q}}_{\text{n}}=\text{0}$;

2. Given the elastic trial stress (Eq. 58) in a point, calculate by PFEM and Gauss Quadrature: ${\sigma }_{\text{n}+\text{1}}^{\text{trial}}=\text{E}\left({\epsilon }_{\text{n}+\text{1}}-{\epsilon }_{\text{n}}^{\text{p}}\right)$;

3. Compute the trial yield function:

3.1. Eq. (66);

4. Verify the yield criterion:

4.1. If ${\text{f}}_{\text{n}+\text{1}}^{\text{trial}}\le \text{0}$, then elastic step;

4.2. Otherwise, plastic correction:

4.2.1. Calculate returning-mapping:

4.2.1.1. Update incremental plastic multiplier, Eq. (67);

4.2.1.2. Update the stress state, Eq. (68);

4.2.1.3. Update the strain state, Eq. (69);

4.2.1.4. Update the parameter $q$, Eq. (70);

4.2.1.5. Update the parameter α, Eq. (71).

5. End plastic correction.

5 EXAMPLES

The purpose of the two numerical examples is to show the limitations of recent time integration schemes to solve dynamic nonlinear problems. Thus, the first numerical example is essentially physically nonlinear, and the second one is geometrically nonlinear only. In both examples, the spectral radii were considered as ρ_∞=0 (complete numerical dissipation) and ρ_∞=0.9 (to avoid convergence problems for most algorithms). Specifically, γ=0.1 and γ=0.4 were applied to the Li and Yu (2020) method since, for γ=0.0, the D_M and D_C vectors would be undefined, because γ appears in the divisor of these expressions, Eqs (32)-(35), in the first sub-step. Similarly, for γ=0.5, α assumes the value 0.0, and the vectors D_M and D_C become undefined again, because α is in the divisor of these expressions, Eqs (40)-(42), in the second sub-step. Both situations are independent of the ρ_∞, even for cases without physical damping (C=0). In all cases, the γ and β parameters for the Newmark method are, respectively, 0.25 and 0.5.

5.1 Toridis-Khozeimeh plane frame

The first example consists of a portal steel plane frame with fixed ends (Figure 4) submitted to a lateral step load (P₀=444.82kN) at the top left corner. Toridis and Khozeimeh (1971), Marur and Kant (1994), Toi and Isobe (1996), Chan and Chui (2000), and Silva et al. (2011), among others, studied this case. Beam and columns are W16x36 and W12x120 steel I profiles, respectively, discretized in three finite elements per bar. Young modulus and yield stress are equal to 207GPa and 470MPa, respectively. The mass of steel profile members is lumped at element nodes and are multiplied by 625. In this work, tolerance for the iterative process and time step is considered 10^-5 and Δt=0.001s, respectively.

Fig. 4
Toridis-Khozeimeh steel plane frame: geometry (left) and load (right).

Results presented in Figure 5, position variation X(t) on applied force degree of freedom, indicate the need of use of at least three finite elements per bar. Amplitude range of position variation indicates inelastic effects. Both classical and recent direct time integration schemes presented practically same results. No significant effect of numerical dissipation (based on parameter ρ_∞) was observed in these responses. The amplitude of vibration with the bilinear stress-strain model with an inelastic branch slope equal to 1/10 of the elastic branch (Toi and Isobe, 1996, 4 elem.) is less than the mixed hardening model (present work, 3 elem.) combined with other time integration methods (Figure 5). It should be emphasized that Toi and Isobe (1996) employed the Adaptively Shifted Integration (ASI) technique to the cubic element based on the Bernoulli-Euler hypothesis. Thus, there is the possibility of the formation of a plastic hinge at Gauss point or even at the ends of the element, which occurred in this example. Therefore, the formulation used by Toi and Isobe (1996) allows more energy dissipation than the one used in the present work. This fact justifies the difference in the results for the stationary vibration period. The other responses coincide with Newmark (Figures 5 and d). For ρ_∞=0.9, the Li and Yu (2020) scheme has a period elongation about Newmark (Figure 5). The WBZ-α and Zhang and Xing (2019) presented convergence problems at t~0.8s (Figure 5) and t~0.5s (Figure 5), respectively, for the same ρ_∞.

Fig. 5
Position versus time (top left corner) for steel plane frame with the mixed hardening model (present work) and bilinear stress-strain model (Toi and Isobe, 1996).

Similar behavior is observed in Figure 6. Results show a close agreement with Chan and Chui (2000) and Silva et al. (2011); see Figures 6a, c, and d. However, the contribution of the high modes in response (ρ_∞=0.9) increased the period for the Li and Yu (2020) algorithm (Figure 6b). It is worth mentioning that inelastic effects predominate over geometrical nonlinearity in this example. Therefore, introducing the contribution of high frequencies to the response (ρ_∞=0.9) did not significantly change the amplitude of positions or even the number of iterations per time step, compared to ρ_∞=0.0.

Fig. 6
Position versus time (top left corner) for steel plane frame with the mixed hardening model (present work) and elastic-perfectly plastic hinge model (Chan and Chui, 2000; Silva et al., 2011).

Two aspects must be highlighted in Figure 7 and 8. The amplitude decay of the Generalized-α and the HHT-α (Figure 7e and f) is caused by the contribution of high frequencies in the response. In addition, there is a period elongation for Li and Yu (2020) response, with a half wavelength lag in t~17s (Figure 8a). The other time integration methods do not present these behaviors. Figures 8b, c, and d indicate an oscillatory behavior of the response obtained from the Li and Yu (2020) and Newmark algorithms, apparently indicating some instability for the physical nonlinearity. However, in the Generalized-α and HHT-α schemes, the amplitude decay of the response over time results in an elliptical spiral (Figures 8c and d) similar to those found in systems with physical damping.

Fig. 7
Amplitude decay (top left corner) for steel plane frame with the mixed hardening model.

Fig. 8
Period elongation and phase space (top left corner) for steel plane frame with the mixed hardening model.

5.2 Three-spring plane frame

The structural system consists of a steel plane frame with four steel bars (Figure 9). The bars have a circular hollow section with 100mm external diameter and 5mm thickness. The system is pinned at the ends and has three rotational springs. Young modulus, specific mass, and spring stiffness are E=200GPa, ρ=7850kg/m³, and k=16881.15177N/m, respectively. One finite element per bar was used, and the magnitude of the periodic triangular load is P₀=30kN. The convergence tolerance for the iterative process is 10^-7. It is intended to evaluate the influence of Δt in the time evolution of the algorithms without viscous damping.

Fig. 9
Three-springs steel plane frame: geometry (left), cross-section (center) and periodic triangular load (right).

Results presented in Figure 10 are related to the vertical position variation on applied force degree of freedom. Initially, only one finite element per bar was used in this numerical example. Amplitude range of position variation indicates strong nonlinear geometrical effects. Figure 10b indicates the need of use of numerical dissipation (based on parameter ρ_∞). It is worth to remember that the smaller is the value of parameter ρ_∞, the higher is the numerical dissipation. It can be seen in Figure 10a that the Generalized-α, Li and Yu (2020), and WBZ-α methods provided a complete time-history response for ρ_∞=0, considering Δt=5x10^-3s. The amplitude of vibration for WBZ-α is the smallest in the reloading interval. For the Li and Yu (2020) method, different range values provided different amplitudes in the reloading stage. For ρ_∞=0.9, no method provided a complete time-history (Figure 10b).

Fig. 10
Vertical position versus time for central spring (Δt=5x10^-3s; 1 element per bar).

On the other hand, the Generalized-α, Li and Yu (2020), γ=0.1, and WBZ-α methods exhibited time-history up to 10s (ρ_∞=0.0; Δt=1x10^-3s). The response curves for Generalized-α and WBZ-α are superimposed over almost the entire analysis time (Figure 11a). For ρ_∞=0.9, only Generalized-α and HHT-α provided a complete time-history response (Figure 11b). The numerical damping in the reloading stage concerning the previous cases is remarkable along with that in the free vibration interval (9s<t<10s).

Fig. 11
Vertical position versus time for central spring (Δt=1x10^-3s; 1 element per bar).

Figure 12a shows that only Generalized-α and WBZ-α provided the complete time-history response and are superimposed before the snap-through (t<2.2s) and at the start of reloading (t>5s). However, for ρ_∞=0.9, only Zhang and Xing (2019) and WBZ-α did not show complete time-history (Figure 12b). It is noteworthy that in free vibration (9s<t<10s), the Generalized-α does not present large amplitudes like the other methods (Figure 12b).

Fig. 12
Vertical position versus time for central spring (Δt=5x10^-4s; 1 element per bar).

For the shortest time interval of this example (Δt=1x10^-4s), two algorithms provided the complete time history for ρ_∞=0: Generalized-α and WBZ-α (Figure 13a). For ρ_∞=0.9, the Generalized-α and HHT-α responses showed a lag between them in reloading, despite the amplitude of the position in both responses being equivalent (Figure 13b). Comparing Figure 10, 12, 13, higher modes contribute significantly to the response by decreasing Δt (ρ_∞=0.0). However, the amplitude increases, diverging from the response for the refined mesh model (5<t<9s).

Fig. 13
Vertical position versus time for central spring (Δt=1x10^-4s; 1 element per bar).

Now, considering a mesh refinement (5 and 10 elements per bar) for the Newmark algorithm, the time-history responses are presented in Figure 14. The spectral radius reduces the response amplitude (5<t<9s) by bringing it closer to the refined mesh model despite the need for a greater number of iterations per time step.

Fig. 14
Newmark and Generalized-α: position versus time (Δt=1x10^-4s).

Table 1 presents the time march progress for the implemented algorithms. Newmark, widely used in practice, did not present a complete time history in any of the evaluated Δt (considering one finite element for each bar), reaching a maximum time march progress of 57.06%. The least efficient method in this case was that of Zhang and Xing (2019),with no complete time history in any of the proposed Δt and with a maximum time march progress of 5.42%. The most efficient method was the Generalized-α, presenting a complete time history in all cases evaluated, except for Δt equal to 5x10^-3s, with the contribution of high frequencies in the response.

Table 1
Time march progress for time integration algorithms (t_f=10s) – three spring plane frame.

Δt (s)	ρ∞		Newmark	Generalized-α	HHT-α	Li and Yu (2020) (γ=0.1)	Li and Yu (2020) (γ=0.4)	Zhang and Xing (2019)	WBZ-α
5x10-3	0.0	t reached (s)	3.395	10	3.375	10	10	0.085	10
		TMP (%)	33.95	100.00	33.75	100.00	100.00	0.85	100.00
	0.9	t reached (s)	-	2.37	2.475	3.79	3.65	0	0.12
		TMP (%)	-	23.70	24.75	37.90	36.50	0.00	1.20
1x10-3	0.0	t reached (s)	5.022	10	4.367	10	5.184	0.055	10
		TMP (%)	50.22	100.00	43.67	100.00	51.84	0.55	100.00
	0.9	t reached (s)	-	10	10	4.144	5.313	0.05	0.049
		TMP (%)	-	100.00	100.00	41.44	53.13	0.50	0.49
5x10-4	0.0	t reached (s)	4.2555	10	5.581	5	5.349	0.304	10
		TMP (%)	42.56	100.00	55.81	50.00	53.49	3.04	100.00
	0.9	t reached (s)	-	10	10	10	10	0.268	0.053
		TMP (%)	-	100.00	100.00	100.00	100.00	2.68	0.53
1x10-4	0.0	t reached (s)	5.7063	10	4.8495	4.7599	4.9109	0.0685	10
		TMP (%)	57.06	100.00	48.50	47.60	49.11	0.69	100.00
	0.9	t reached (s)	-	10	10	5.3361	5.3423	0.542	0.0792
		TMP (%)	-	100.00	100.00	53.36	53.42	5.42	0.79

As the system is geometrically nonlinear, subjected to forced vibration, and numerically damped, the phase space presents characteristics of systems with chaotic behavior (Figure 15). For ρ_∞=0.0, the range of velocities and positions (concerning the undeformed configuration) indicate accumulation of kinetic and strain energy for the Generalized-α and HHT-α algorithms, without amplitude decay effect. The same occurs with Newmark. There is a significant reduction in velocities and positions, indicating energy loss caused by numerical damping (ρ_∞=0.9). The phase space diagrams demonstrate that the responses are equal for Generalized-α and WBZ-α (ρ_∞=0.0) and for Generalized-α and HHT-α (ρ_∞=0.9). Therefore, the Generalized-α algorithm seems to incorporate the best numerical dissipation characteristics of the HHT-α and WBZ-α methods for the evaluated conditions.

Fig. 15
Phase space for Newmark, Generalized-α, HHT-α, and WBZ-α methods at central spring (Δt=1x10^-4s, 1 element per bar).

6 CONCLUSION

This work presented a comparison of the time integration algorithms implemented using the Positional Finite Element Method (PFEM) formulation and the mixed hardening model scheme applied to the same formulation. The mixed hardening model showedgood agreement with the inelastic models of the indicated literature (Toridis and Khozeimeh, 1971; Marur and Kant, 1994; Toi and Isobe, 1996; Chan and Chui, 2000; Silva et al., 2011) besides its simplicity of implementation and good agreement with experimental results (Greco, 2004). Regarding the time integration algorithms, it was observed that the discretization of a structure with few finite elements can lead to a poor model for the higher modes of vibration, and its influence on the response is highly significant. Numerical damping becomes very useful to reduce this imprecision in the response despite the need for more iterations per time step.

From the two numerical results evaluated, the first mainly ruled by inelastic behavior and the second ruled by nonlinear geometrical behavior, it is possible to conclude that the Generalized-α is the most efficient algorithm to deal with these problems (i.e., it is stable and simpler to implement than newer algorithms), and WBZ-α and three-parameter single-step implicit method are the least efficient algorithms to deal with these problems. In addition, the mixed hardening model presented a good agreement with the proposed example, but it is recommended to carry out further analyses considering the two forms of nonlinearity.

Therefore, the PFEM formulation with the time integration algorithms and the mixed hardening model is a significant and viable alternative to the geometrical and physical nonlinear analysis of plane frames. Finally, it is worth mentioning that considerable contributions have been made recently in terms of new time integration algorithms. However, as is the case with any new techniques, these algorithms must be tested for more complex phenomena. The present paper tested some of recent time integration algorithms for two severe nonlinear conditions, and the results obtained showed that classic time integration algorithms still maintain their value.

Acknowledgments

The authors are grateful to CAPES (Coordination of Improvement of Higher Education Personnel), CNPq (National Council for Scientific and Technological Development), FAPEMIG (Foundation for Research Support of Minas Gerais State), UFMG (Federal University of Minas Gerais), PUC Minas (Pontifical Catholic University of Minas Gerais) and UFOP (Federal University of Ouro Preto) for the opportunity to develop this work.

• https://doi.org/10.1590/1679-78256973

References

• Bathe, K. J. (2007). Conserving energy and momentum in nonlinear dynamics: a simple implicit time integration scheme. Comp. and Str. 85(7-8):437-445.
• Bathe, K. J. (2014). Finite Element Procedures. Prentice Hall, Watertown.
• Behnoudfar, P., Deng, Q., Calo, V. M. (2020). High-order generalized-alpha method. Applic. Eng. Sci. 4:100021.
• Borrvall, T., Lawson, M. R. (2020). On Accuracy and Stability of Implicit Time Integration Schemes for Rotating Structures. 16th Int. LS-DYNA® Users Conference, Virtual Event.
• Chan, S. L., Chui, P. T. (2000). Non-linear static and cyclic analysis of steel frames with semi-rigid connections. Elsevier, Oxford.
• Chung, J., Hulbert, G. M. (1993). A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalized-α Method. J. Appl. Mech. 60:371-375.
• Clough, R. W. (2003). Dynamics of Structures. Computers and Structures, Berkeley.
• Coda, H. B., Greco, M. (2004). A simple FEM formulation for large deflection 2D frame analysis based on position description. Comp. Meth. Appl. Mech. and Eng. 193:3541-3557.
• Coda, H. B., Paccola, R. R. (2008). A positional FEM formulation for geometrical non-linear analysis of shells. Latin American Journal of Solids and Structures, 5(3):205-223.
• Coda, H. B. (2009). Two dimensional analysis of inflatable structures by the positional FEM. Latin American Journal of Solids and Structures, 6(3):187-212.
• Coda, H. B., Paccola, R. R. (2010). Improved finite element for 3D laminate frame analysis including warping for any cross-section. Applied Mathematical Modelling, 34(4):1107-1137.
• Coda, H. B., Paccola, R. R. (2011). A FEM procedure based on positions and unconstrained vectors applied to non-linear dynamic of 3D frames. Finite Elements in Analysis and Design, 42(4):319-333.
• Geng, L., Tu, S. T., Gong, J., Jiang, W., Zhang, W. (2018). On residual stress and relief for an ultra-thick cylinder weld joint based on mixed hardening model: numerical and experimental studies. J. Press. Ves. Tech. 140(4):041405.
• Géradin, M., Rixen, D. J. (2015). Mechanical Vibrations: Theory And Application To Structural Dynamics, 3th ed. John Wiley and Sons, United Kingdom.
• Greco, M. (2004). Análise de problemas de contato/impacto em estruturas de comportamento não linear pelo método dos elementos finitos, PhD Thesis (in Portuguese) University of São Paulo, Brazil.
• Greco, M., Coda, H. B. (2006). Positional FEM formulation for flexible multi-body dynamic analysis. Journal of Sound and Vibration, 290:1141-1174.
• Hilber, H. M., Hughes T. J. R., Taylor, R. L. (1977). Improved numerical dissipation for time integration algorithms in structural dynamics. Earth. Eng. and Str. Dyn. 5(3):283-292.
• Jiang, W., Xie, X., Wang, T., Zhang, X., Tu, S. T., Wang, J., Zhao, X. (2021). Fatigue life prediction of 316L stainless steel weld joint including the role of residual stress and its evolution: Experimental and modelling. Int. J. Fatig. 143:105997.
• Kim, W. (2020). A Comparative Study of Implicit and Explicit Composite Time Integration Schemes. Int. J. Str. Stab. and Dyn. 20(13):25.
• Kontoe, S., Zdravković, L., Potts, D. M. (2004). Performance of the generalized-α integration method in dynamic geotechnical problems. Ninth International Symposium on Numerical Models in Geomechanics - NUMOG IX:211-216.
• Kuo, S., Yau, J. D., Yang, Y. B. (2012). A Robust Time-integration Algorithm For Solving Nonlinear Dynamic Problems With Large Rotations And Displacements. Int. J. Str. Stab. Dyn. 12(6):24.
• Li, J., Yu, K. (2020). A truly self-starting implicit family of integration algorithms with dissipation control for nonlinear dynamics. Nonlin. Dynam. 102:2503-2530.
• Liu, Y., Liu, J., Herrmann, M., Schenck, C., Kuhfuss, B. (2021). Material Flow in Infeed Rotary Swaging of Tubes. Materials. 14(1):58.
• Malakiyeh, M. M., Shojaee, S., Bathe, K. J. (2019). The Bathe time integration method revisited for prescribing desired numerical dissipation. Comp. and Str. 212:289-298.
• Marur, S. R., Kant, T. (1994). A stress correction procedure for the analysis of inelastic frames under transient dynamic loads, Comp. and Str. 50(5):603-613.
• Noh, G., Bathe, K. J. (2019a). The Bathe time integration method with controllable spectral radius: The ρ_∞-Bathe method. Comp. and Str. 212:299-310.
• Nogueira, C. G., Venturini, W. S., Coda, H. B. (2013). Material and geometric nonlinear analysis of reinforced concrete frame structures considering the influence of shear strength complementary mechanisms. Latin American Journal of Solids and Structures, 10(5):953-980.
• Pascon, J. P., and Coda, H. B. (2013). A shell finite element formulation to analyze highly deformable rubber-like materials. Latin American Journal of Solids and Structures, 10(6):1177-1209.
• Rostami, S., Shojaee, S., Saffari, H. (2013). An explicit time integration method for structural dynamics using cubic B-spline polynomial functions. Scient. Iran. A, 20(1):23-33.
• Seisenbacher, B., Winter, G., Grün, F. (2018). Improved Approach to Determine the Material Parameters for a Combined Hardening Model. Mat. Sci. and Applic. 9:357-367.
• Shi, J. J., Meng, B., Cheng C., Wan, M. (2020). Size effect on the subsequent yield and hardening behavior of metal foil. Int. J. Mech. Sci. 105686.
• Silva, A. R. D., Fernandes, W. L., Silveira, R. A. M., Gonçalves, P. B. (2011). Elastoplastic analysis of plane steel frames under dynamic loading. XI International Conference on Computational Plasticity - Fundamentals and Applications (COMPLAS XI), Barcelona, Spain.
• Silva, W. Q., Coda, H. B. (2012). Numerical combination for nonlinear analysis of structures coupled to layered soils. Latin American Journal of Solids and Structures, 9(2):235-257.
• Simo, J. C., Hughes, T. J. R. (2006). Computational inelasticity, Springer Science, New York.
• Soares Jr., D. (2015). A simple and effective new family of time marching procedures for dynamics. J. Comp. Meth. Appl. Mech. and Eng. 283:1138-1166.
• Toi, Y., Isobe, D. (1996). Finite element analysis of quasi-static and dynamic collapse behaviors of framed structures by the adaptively shifted integration technique, Comp. and Str. 58(5):947-955.
• Toridis, T. G., Khozeimeh, K. (1971). Inelastic response of frames to dynamic loads. J. Eng. Mech. Div. 97(3):847-863.
• Wood, W. L., Bossak, M., Zienkiewicz, O. C. (1980). An alpha modification of Newmark's method. Int. J. for Num. Meth. Eng. 15(10):1562-1566.
• Yang, H., Zhang, W., Zhuang, X., Zhao, Z. (2020). Calibration of Chaboche Combined Hardening Model for Large Strain Range. Proc. Manufact. 47:867-872.
• Zhang H., Xing, Y., (2019). A three-parameter single-step time integration method for structural dynamic analysis. Act. Mech. Sin. 35(1):112-128.
• Zhang, J. (2020). A-stable two-step time integration methods with controllable numerical dissipation for structural dynamics. Int. J. Num. Meth. Eng. 121:54-92.
• Zhao, L., Xu, L., Han, Y., Jing, H., Gao, Z. (2019). Modelling creep-fatigue behaviors using a modified combined kinematic and isotropic hardening model considering the damage accumulation. Int. J. Mech. Sci. 161:105016.

Edited by

Editor: Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  04 May 2022
• Date of issue
  2022

History

• Received
  21 Jan 2022
• Reviewed
  31 Mar 2022
• Accepted
  06 Apr 2022
• Published
  08 Apr 2022