Development and Application of a High-Performance Triangular Shell Element and an Explicit Algorithm in OpenSees for Strongly Nonlinear Analysis

: The open-source finite element software, OpenSees, is widely used in the earthquake engineering community. However, the shell elements and explicit algorithm in OpenSees still require further improvements. Therefore, in this work, a triangular shell element, NLDKGT, and an explicit algorithm are proposed and implemented in OpenSees. Specifically, based on the generalized conforming theory and the updated Lagrangian formulation, the proposed NLDKGT element is suitable for problems with complicated boundary conditions and strong nonlinearity. The accuracy and reliability of the NLDKGT element are validated through typical cases. Furthermore, by adopting the leapfrog integration method, an explicit algorithm in OpenSees and a modal damping model are developed. Finally, the stability and efficiency of the proposed shell element and explicit algorithm are validated through the nonlinear time-history analysis of a high-rise building.


Introduction
The performance of structures against extreme hazards has become an important research topic. By discovering the damage evolution process and failure mechanism, the research outcomes will support the identification and optimization of vulnerable structures. In addition to physical experiments, numerical simulations based on the finite element method, as an important and effective approach, have been widely used [Nesnas and Abdul-Latif (2001); Bradford and Pi (2012) ;Lin, Li, Lu et al. (2016)]. Thus far, strongly nonlinear analyses of structures have been performed extensively, and corresponding simulation strategies have been proposed [Lu, Lu, Guan et al. (2013); Lu, Tian, Cen et al. (2018)]. OpenSees, as an open-source finite element software, is now widely used owing to its high transparency and freedom [McKenna, Scott and Fenves (2009)]. For strongly nonlinear problems, on the one hand, the elements must consider the material and geometric nonlinearity simultaneously; on the other hand, the time integration algorithm should be sufficiently stable during the entire computational process. However, further improvement on these two aspects is still required in OpenSees. In terms of element technology, a typical modeling strategy for the nonlinear analysis of buildings is to adopt fiber elements for beams/columns and shell elements for shear walls and coupling beams [Lu and Guan (2017)]. In OpenSees, the collapse simulation of frame structures has been performed successfully by using fiber elements [Lignos, Chung, Nagae et al. (2011);Xie, Lu, Guan et al. (2015)]. However, it would be difficult for the fiber elements to represent the axial-flexural-shear coupled behavior of shear walls. Therefore, based on the generalized conforming theory, Lu et al. [Lu, Tian, Cen et al. (2018)] proposed and successfully implemented a quadrilateral flat shell element, NLDKGQ, into OpenSees. Subsequently, they performed a collapse simulation of a highrise reinforced concrete (RC) frame-core tube building using NLDKGQ. The NLDKGQ element, consisting of the plate element DKQ and the membrane element GQ12, can avoid shear locking. By introducing the updated Lagrangian formulation, NLDKGQ can simulate the geometric nonlinearity and is suitable for large deformation problems. However, NLDKGQ is not adaptable to triangular meshes; therefore, it is not easy to adopt the NLDKGQ element for the cases with complicated boundaries or curved surfaces. In contrary, triangular shell elements are more adaptive to complicated boundaries, and they can effectively solve mesh distortion and warpage problems. Therefore, it is necessary to propose a triangular shell element for OpenSees. In terms of the time integration algorithm, two types of algorithms exist: implicit algorithm and explicit algorithm. The implicit algorithm is typically used, but a convergence test is essential at each time step. It is noteworthy that the implicit algorithm may fail to perform a complete analysis owing to its strong nonlinearity-induced nonconvergence. Therefore, explicit algorithms are preferred for strong nonlinearity [Lu, Lin, Cen et al. (2015); Pham, Tan and Yu (2017)], which can avoid convergence problems. Among the existing explicit algorithms, the central difference method is the most popular one. Theoretically, the central difference method can be highly efficient when the system of equations can be decoupled. However, the decoupling criterion for the central difference method requires a diagonal damping matrix. The mass-proportional damping matrix is diagonal, but it obviously underestimates the damping ratio of high-order modes and consequently does not yield a satisfactory accuracy [Xie (2015)]. In contrast, if the stiffness-proportional damping model is introduced to restrain high-order modes, the system of equations would fail to decouple, leading to an increased computational time.
To solve the problems above, researchers have proposed numerous solutions. For example, Li et al. [Li, Liao and Du (1992)] derived an explicit difference method for viscoelastic dynamic equations; Du et al. [Du and Wang (2000)] derived an explicit integration formula for damped elastic lumped-mass structures. However, although the algorithms proposed by Li et al. [Li, Liao and Du (1992)] and Du et al. [Du and Wang (2000)] can ensure the decoupling of the system of equations, the equations for displacement and velocity are required to be established and solved separately at each time step, which significantly increases the computational time. Consequently, based on the generalized conforming theory and the updated Lagrangian formulation [Long, Cen and Long (2009)], a new triangular shell element NLDKGT is proposed in this work that is suitable for cases with complicated boundary conditions and problems with strong nonlinearities. Furthermore, by adopting the leapfrog integration method, an explicit algorithm in OpenSees and a modal damping model are developed in this work. The explicit integration algorithm can ensure the decoupling of the system of equations. The accuracy and reliability of the triangular shell element are validated through typical cases. Finally, the stability and efficiency of the proposed shell element and the explicit algorithm are validated through the nonlinear time-history analysis of a high-rise building.

Basic formulation under small deformation
To develop a suitable triangular shell element, the triangular planar membrane element GT9 [Xu and Long (1993)] and the triangular thin plate element DKT [Batoz, Bathe and Ho (1980)] were used to construct the new triangular shell element in this work. The planar membrane element GT9 contains three degrees of freedom (DOFs) at each node by introducing a rigid rotational freedom. In addition, a higher accuracy is achieved by defining higher order displacement fields [Xu and Long (1993)]. The plate element DKT is based on the Kirchhoff theory and can effectively avoid shear locking. Fig. 1 illustrates the decomposition of the NLDKGT element. Consisting of GT9 and DKT, the NLDKGT has six DOFs at each node. This will greatly reduce the connection modeling workload among the shell and the beam/column elements. The nodal displacement q is defined as follows: Here, qi m and qi b denote the nodal displacement components of GQ9 and DKT, respectively. They can be expressed as follows: The displacement u of GQ9 can be obtained through superposition of u0 and uθ as shown in Eq. (4). Here, u0 denotes the linear part of the displacement field, while uθ is an additional rotation displacement. Through Eqs. (5) to (8), u0 and uθ can derived [Xu and Long (1993)].
is the coordinate of node i in GT9 element in the local system. Then, Eqs. (9) and (10)  The stiffness matrix of GQ9 is as follows: Here, Dmm represents the material matrix of GQ9 element. Generally, if the element is made of isotropic linear elastic materials, Eq. (12) can be adopted to derive Dmm. In Eq. (12), E represents the elastic modulus, h represents the element thickness, and ν represents Poisson's ratio. Eq.
(3) defines the rotational DOFs of the DKT element [Batoz, Bathe and Ho (1980)]. The relation between nodal displacements q b and rotational strain χb is as follows: Details on H x and H y are available in Batoz et al. [Batoz, Bathe and Ho (1980)].
The stiffness matrix of DKT is as follows: Here, Dbb denotes the material matrix of DKT. Generally, if the element is made of isotropic linear elastic materials, Dbb can be derived as: For small deformations, based on the plate stiffness matrix Kb and the membrane stiffness matrix Km, the local stiffness matrix of the NLDKGT element can be derived according to the DOF sequencing in Eq. (1). Then, the global element stiffness matrix can be obtained through coordinate transformation.

Geometric nonlinearity
At each time step t, through the updated Lagrangian formulation, the current deformation can be adopted to update the stresses and strains in incremental forms. Based on Kirchhoff's and von Karman's assumptions [Podio-Guidugli (1989)], a linear part (∆e) and a nonlinear part (∆η) constitute the shell element strain increment (∆ε). The linear part (∆e) can be derived from plate rotational strain increment ∆χb and membrane strain increment ∆εm, as shown in Eqs. (16) and (17), respectively.
From t to t + dt, Eq. (18) is adopted to update the shell element stress tensor: Here, at time t, Dtan is the tangential constitutive matrix. In the local coordinate system of the shell element, the system of equations using the updated Lagrangian formulation is as follows: Here, t+dtF and tR represents the external and internal force vector, respectively; the subscript represents the time step. In terms of the stiffness matrix, Kl is the linear part and Knl is the nonlinear part, which can be obtained through Eqs. (20) and (21), respectively.

∫∫
Here, by integrating Dtan through the element thickness as shown in Eq. (22), Dmm, Dmb, Dbm, and Dbb can be solved. Using the bending plate element interpolation function, the matrix G can be derived [Batoz, Bathe and Ho (1980)]. Variables corresponding to the membrane element internal force vector constitute the matrix N t (Eq. (23) Eq. (24) can be used to solve the elemental internal force vector tR in Eq. (19):

Implementation in OpenSees
Under the class of shell in OpenSees, a new class named ShellNLDKGT is added. During the implementation of the NLDKGT, nearly no source code change is made beyond the shell element domain. Through the official website of OpenSees (http://opensees.berkeley.edu), users can download corresponding source code of the NLDKGT element. It is worth noting that, the proposed NLDKGT element is compatible with other elements in OpenSees. Thus, for a real finite element model, users can use the NLDKGT elements in complicated boundary areas, and use other four-node shell elements in regular-shaped areas.

Scordelis-Lo roof problem
The Scordelis-Lo roof problem is shown in Fig. 2. The cylindrical panel is loaded vertically by a uniform dead weight of g=90. The panel is supported by end diaphragms but the sides are free. Owing to the symmetry, only one quarter of the panel is established. Three types of meshes were adopted in this analysis, as listed in Tab. 1. The vertical deflection at point A was recorded. For this case, the geometric nonlinearity was not considered. The exact solution of 0.3024 provided by MacNeal and Harder [MacNeal and Harder (1985)] was used as a reference. The results obtained using the DKT-CST-15RB element [Nicholas, Henryk and Ted (1986)] and OLSON element [Olson and Bearden (1979)] were compared with the results obtained using the NLDKGT element. The DKT-CST-15RB element is a superposition of the DKT plate bending element and the CST plane stress element, with 15 DOFs [Nicholas, Henryk and Ted (1986)]. The OLSON element is an 18-DOF flat triangular shell element reformulated by combining a bending triangle with a plane stress triangle incorporating in-plane rotations at each vertex [Olson and Bearden (1979)]. Tab. 1 shows the comparison results. The NLDKGT element is more accurate compared to the other two elements. Figure 2: Scordelis-Lo roof problem  Fig. 3 shows the twisted beam problem [MacNeal and Harder (1985)]. A concentrated load is applied at the tip along the in-plane (P) and out-of-plane (Q) directions, respectively. A mesh of 2×12 was adopted in this problem. Two load cases were analyzed: (1) P=1, Q=0; and (2) P=0, Q=1. The displacement along the loading direction at the tip was recorded. For this case, the geometric nonlinearity was not considered. The exact solutions provided by MacNeal et al. [MacNeal and Harder (1985)] were used as a reference. Tab. 2 shows the results of the comparison and illustrates the accuracy of the NLDKGT element. Figure 3: Twisted beam problem

Large deformation problem of a cantilever beam
To validate the geometric nonlinearity simulation capacity of the NLDKGT element, a cantilever beam subjected to a pure bending load (out-of-plane) is analyzed [Horrigmoe and Bergan (1978); Park, Cho and Lee (1995)], as shown in Fig. 4. The mesh of 1×10 is adopted. Fig. 5(a) shows the relationship between the normalized moment (κ=M/Mmax) and the horizontal and vertical displacements at the loading point. Fig. 5(b) shows the deformed shape of the cantilever beam under different bending moments. The results show that the NLDKGT element can simulate the large deformation and rotation problems with good accuracy, which is similar to the S4 element in ABAQUS. Such a large deformation capacity makes the NLDKGT element highly suitable for geometric nonlinearity problems.

Buckling analysis of an H-shaped beam
An H-shaped beam shown in Fig. 6 is used to demonstrate the buckling analysis. An isotropic elastic material (E=2.06×10 11 Pa, ν=0.3) is used for the beam. Both the NLDKGQ and NLDKGT elements were adopted in this analysis, and the corresponding meshes are shown in Fig. 6. The two ends of the H-shaped beam were simply supported. In finite element simulations, initial defects (e.g., initial bow imperfections leading to additional moment to the middle of components) are theoretically necessary to simulate the buckling phenomenon. According to the recommendations in EN 1993-1-1 [CEN (2005)], a distributed load of p=0.5 N was imposed at each node on the web of the Hshaped beam to simulate the initial defects. Subsequently, a pressure load was applied at the top of the H-shaped beam, and the relation between the vertical load and displacement along the loading direction was recorded, as shown in Fig. 7. As shown in Fig. 7, the model using NLDKGQ fails to converge when the imposed load approaches 200 kN, i.e., when the H-shaped beam just begins to buckle. This phenomenon is due to the warping of the quadrilateral shell element. In contrast, a stable result is obtained using the triangular shell element NLDKGT. Fig. 8 shows the deformation of the H-shaped beam along the Z direction. The blue and red solid lines denote the deformation shape using the NLDKGT and NLDKGQ elements, respectively, at the time step when NLDKGQ fails to converge. The dashed blue line denotes the final deformation of the model using the NLDKGT element. Through the analysis of this case, the NLDKGT element is proven as more stable and reliable than the NLDKGQ element for the buckling analysis.

RC shear wall experiments
To investigate the performance of the NLDKGT element in simulating RC specimens, the hysteretic behavior of two shear wall specimens is analyzed using OpenSees based on the multilayered shell section model ] and the NLDKGT element. The test specimens include one rectangular wall (denoted as SW1-1) [Zhang (2007)], and one coupled wall (denoted as CW-3) [Chen and Lu (2003)]. The meshing schemes and corresponding hysteretic curves are shown in Fig. 9. The comparisons between the test and simulation results indicate that the NLDKGT element can provide satisfactory simulation results in the nonlinear behavior of RC shear walls.  28).
where M and C are the mass and damping matrices of the system, respectively; R and P are the resisting and external force vectors of the system, respectively; u, u  , and u   denote the displacement, velocity, and acceleration, respectively; the subscript denotes the time.
It is difficult to decouple the equations if C is not a diagonal matrix. To avoid this problem, most researchers adopt the mass-proportional damping model for the central difference method. However, mass-proportional damping will underestimate the damping ratio of high-order vibration modes, which sometimes leads to unreasonable results.

Leapfrog integration method
The leapfrog integration method [Hockney (1970)] is an improved format proposed based on the Verlet integration method [Verlet (1967)]. In the leapfrog method, the equations for updating velocity and displacement are as follows: To adopt the leapfrog method, Eq. (27) Substituting Eqs. (29)-(30) and Eq. (32) into Eq. (31) yields Eq. (33) shows that the system of equations can be decoupled when the mass matrix is diagonal. However, in this method, the kinetic and potential energies of the system are not defined at the same time step, leading to a failure in calculating the total energy directly. To solve this problem, certain additional steps are added to revise the algorithm. The entire process of the revised format is as follows [Sandvik (2018)]: (1) First, calculate t t ∆ + u through Eq. (33); (2) Subsequently, calculate t u  again through the central difference method using t t ∆ + u in Step (1); (3) Solve t u   using the newly obtained t u  : (4) Solve t t ∆ + u using Eqs. (29)-(30).
The revised format above was performed through an iterative process. However, the additional computational cost is still relatively small, because the equations are simple. In addition, the revised format will provide the velocity and displacement at the same time step and is convenient to calculate the total energy at each time step directly. Because the backward difference format is adopted for the velocity, the stability criterion of the algorithm is different from the central difference method. Here, the conclusion will be given directly as follows (more details are provided in Appendix A): where, ωn is the highest angular frequency of the system; Tn is the shortest period of the system; ζ is the damping ratio corresponding to ωn. Eq. (35) shows that the numerical stability of the algorithm is not only related to the system frequency but also to the damping ratio.
For a finite element model, the shortest period Tn can be determined by solving generalized eigenvalues of the system. But, to simplify this procedure, an additional method is often adopted: to solve the shortest period of each element (denoted as min(Tn (e) )) [Wang (2003)]. It has been proved that, the substitute period min(Tn (e) ) is always not longer than Tn. The min(Tn (e) ) for each element is usually approximately estimated by using πL/C. Here, L is the characteristic length of the element; C is the wave speed. These parameters may differ for different kinds of elements. For example, for truss and beam elements, L is the length of the element, and C can be taken as ρ E , where E is the Young's modulus, and ρ is the mass density. For shell elements, three kinds of L are provided by Hallquist [Hallquist (2006)], and C can be taken as , where ν is the Poisson's ratio.
Although different estimation methods can be found for min(Tn (e) ), the basic concept is identical: The stable time step will be smaller for models with smaller element sizes and larger stiffness. Therefore, appropriate meshing schemes should be adopted for models using explicit algorithms. The explicit algorithm above was implemented in OpenSees through a new class called Explicitdifference, which falls under the class of Integrator. The new algorithm fits the OpenSees framework. The source code of the method is available at the official website of OpenSees (http://opensees.berkeley.edu).

Damping model adopted in explicit algorithm
To avoid a high computational cost, restrain unreasonably high frequency vibrations, and ensure the equations to be decoupled, it is necessary to use the superposition of the modal damping and mass-proportional damping models. Thus, the modal damping model was implemented in OpenSees. The modal damping can be expressed as follows [Clough and Penzien (2003)]: where Cm is the modal damping matrix; M is the mass matrix; mn, ζn, ωn, and φn are the modal mass, modal damping ratio, natural vibration frequency, and mode shape corresponding to the nth mode, respectively. The damping model that this work adds to OpenSees ensures that the superposition of the damping ratios from the modal damping, and the mass-proportional damping of each vibration mode is equal to the assigned total damping ratio.

Collapse simulation of a high-rise RC frame-core tube building
In this section, a 42-story RC frame-core tube building with a height of 141.8 m (Fig. 10) (denoted as Building 2N by Lu et al. [Lu, Xie, Guan et al. (2015)]) was simulated using OpenSees. More details about this building are provided by Lu et al. [Lu, Xie, Guan et al. (2015)]. The beams and columns were simulated with fiber beams. The shell element combined with the multilayered shell section was adopted to model the shear walls. Hence, both the material and geometric nonlinearity can be considered. The shear walls in this high-rise building are of regular shapes. In this work, coupling beams of the core tube were simulated using NLDKGT, while other shear walls were modeled using NLDKGQ.  ] First, the El-Centro 1940 ground motion was adopted as the input along the X direction. The peak acceleration was adjusted to 5.1 m/s 2 (2% probability of exceedance in 50 years, as defined in the Chinese code [CMC (2010)]). According to the mesh size, material property, and element size, the time step was set to 4×10 -5 s for the explicit algorithm, and 0.01 s for the implicit algorithm. Tab. 3 provides the information of the analyzed cases.   Fig. 11, the superposition of the mass-proportional and modal damping (i.e., the Ex-MS+MD model) will provide similar results to the implicit algorithm using the Rayleigh damping (i.e., the Im-RL model). However, if only the mass-proportional damping model is adopted, the IDR results are much greater.
For large-scale engineering structures, the fundamental periods are relatively long. Consequently, the high-order vibration modes contribute significantly to the structural responses. Thus, the mass-proportional damping model alone, to some extent, is not suitable for large-scale structures. It is more appropriate to adopt the superposition of the modal damping and mass-proportional damping models to avoid unreasonably high frequency vibrations. An incremental dynamic analysis (IDA) was performed using the Ex-MS+MD10 model, and the peak accelerations were adjusted to 5.1 m/s 2 , 20 m/s 2 , 40 m/s 2 , and 50 m/s 2 , respectively. Here, 20% of the initial slope is used to find the collapse intensity [FEMA (2000); Jalayer (2003); Villaverde (2007)]. According to the criterion above, Building 2N will collapse when the peak acceleration of the El-Centro record is larger than 40 m/s 2 .
(a) Inter-story drift envelope (b) Relation between PGA and the maximum IDR Tab. 4 shows the efficiency comparison among different cases in Tab. 3 and an additional case (explicit algorithm+Rayleigh damping). Among all the cases using the explicit algorithm, the Ex-MS, and Ex-MS+MD models both require the least computational time cost. However, when the Rayleigh damping model is adopted for an explicit algorithm, the time cost becomes much larger (approximately three times that of the Ex-MS model).
The primary reason of this phenomenon is that, the Ex-RL model spends significantly more time on the damping matrix contributed by the stiffness matrix at each time step. The total computational cost of the two cases using the implicit algorithm is less than that of the explicit algorithm. It is noteworthy that, when the implicit algorithm is adopted to perform the collapse analysis (i.e., the Im-RL2 model), the computational time will increase significantly. This is because the number of iterations will increase significantly when the structural components enter strong nonlinearity. Even with a larger convergence tolerance, the average time cost at each step is still 2.4 times that of the Im-RL1 model (which has a smaller ground motion intensity). The explicit algorithm does not require any iteration. This advantage means that the time cost of the explicit algorithm is proportional to the number of time steps. Thus, for strongly nonlinear problems, compared with the explicit algorithm, the implicit algorithm requires more time cost, and sometimes demonstrates no satisfactory results because of convergence failure.
Although the Ex-MS model also demands less computational time, its results are not accurate because the mass-proportional damping alone cannot control the unnecessary high-order vibration. Thus, the Ex-MS+MD model is the best option for the collapse analysis of this building.

Conclusions
For strongly nonlinear analysis, the element and time integration algorithm are two important challenges. However, the shell elements and the explicit algorithm in OpenSees still require further improvements. Therefore, a triangular shell element NLDKGT and an explicit algorithm are proposed and implemented in OpenSees in this work. The conclusions are as follows: (1) Through the validation of classical benchmarks, the triangular shell element NLDKGT was proven accurate and reliable. Compared with the quadrilateral element, the NLDKGT element could not only well consider the geometric nonlinearity, but also exhibited great advantages in strong nonlinear and warpage problems, such as buckling analysis. In addition, it is more flexible to use NLDKGT elements in complicated boundary areas to avoid mesh distortion; (2) An explicit algorithm, along with a modal damping model, was implemented into OpenSees based on the leapfrog method. Through the nonlinear time history analysis of a high-rise RC frame-core tube building, the proposed shell element and explicit algorithm demonstrated higher efficiency and more stable results in strong nonlinear problems.