OpenSees Three-Dimensional Computational Modeling of Ground-Structure Systems and Liquefaction Scenarios

The OpenSees computational platform has allowed unprecedented opportunities for conducting seismic nonlinear soil-structure interaction simulations. On the geotechnical side, capabilities such as coupled solid-fluid formulations and nonlinear incrementalplasticity approaches allow for representation of the involved dynamic/seismic responses. This paper presents recent research that facilitated such endeavors in terms of response of ground-foundation-structure systems using advanced material modeling techniques and high-performance computing resources. Representative numerical results are shown for large-scale soil-structure systems, and ground modification liquefaction countermeasures. In addition, graphical user interface enabling tools for routine usage of such 3D simulation environments are presented, as an important element in support of wider adoption and practical applications. In this context, Performance-Based Earthquake Engineering (PBEE) analysis of bridge-ground systems is highlighted as an important topical application.


PressureIndependMultiYieldSoftening material
Based on the existing multi-yield surface J2 associative plasticity formulation [Prevost (1977); Elgamal, Yan, Yang et al. (2008);Lu, Elgamal, Yan et al. (2011)], the PressureIndependMultiYield material ( Fig. 1) has been implemented in OpenSees for simulating elasto-plastic undrained clay-type shear response. In this model, a hyperbolic backbone curve 1 / r Gγ τ γ γ = + can be employed to describe the backbone shear stressstrain curve ( Fig. 1(b)), where G is the low-strain shear modulus, τ and γ denote the octahedral shear stress and shear strain, respectively, and r γ is the reference shear strain (computed as , in which max τ is the peak shear strength that corresponds to the shear strain max γ ). Alternatively, conventional shear modulus reduction (G/Gmax) curves may be conveniently employed for the specification of yield surfaces and their characteristics. Further details and model implementation specifics can be found in Prevost [Prevost (1977)  Extending the PressureIndependMultiYield material with a newly developed strain softening logic, a PressureIndependMultiYieldSoftening material was recently developed and calibrated based on a series of laboratory sample test data [Qiu (2020)]. For this model, Fig. 2(a) presents a range of response scenarios that can be generated by the incorporated softening logic [Qiu (2020)]. As such, representation of strain softening is an important consideration for a wide range of ground formations involving sensitive clays, cemented soils, over-consolidated soils, very dense soils, or frozen soils among others. Fig. 2(b) depicts the model response under cyclic loading with an imposed static shear stress bias. As can be seen in Fig. 2(b), the model reproduces the downslope cycle-by-cycle accumulation of shear deformations and the degradation of shear stiffness and strength.

PressureDependMultiYield03 material
The PressureDependMultiYield02 material (Fig. 3) has mechanisms to simulate the dilatancy and cyclic mobility of pressure sensitive soils ; Yang, Elgamal and Parra et al. (2003)]. Calibration was performed based on a series of laboratory and centrifuge tests, and the model parameters were provided for sands with different relative densities. Representative simulation results using the PressureDependMultiYield02 material for a symmetric stress-controlled cyclic shear loading test are shown in Fig. 3

LadeDuncanMultiYield material
To allow for further accuracy in capturing 3D shear response, the PressureDependMultiYield03 material has been extended to incorporate the Lode angle effect ( Fig. 5(a)) by employing the Lade-Duncan failure criterion as the yield function PDMY03 Dr = 75% [Lade and Duncan (1975); ; Qiu (2020)]. This failure criterion is represented by Chen et al. [Chen and Mizuno (1990)]: where, I1 is volumetric stress, J2 and J3 are second and third deviatoric stress invariants, respectively, parameter k1 (>27) is related to soil shear strength (or friction angle φ). A typical yield surface fm ( Fig. 5(a)) is defined by Yang et al. ]: where ηm is a normalized yield surface size (0<ηm<1) and material parameter 1 =  ]. In Eq.
(2), 2 1 : I is 2 nd unit tensor, deviatoric tensor a is back stress (yield surface center), and the operators "·"and ":" denote single and double contraction of two tensors, respectively. For illustration, the LadeDuncanMultiYield material was calibrated using laboratory cyclic undrained stress-controlled test data of a reclaimed gravely fill Masado soil [Hatanaka, Uchida and Ohara (1997)]. It can be seen that the soil gradually loses its effective confinement during cyclic triaxial loading, and the axial strain accumulates on a cycle-bycycle basis (Fig. 5(b)) after liquefaction, with a strong dilative tendency. It is noted that this LadeDuncanMultiYield material was recently updated [Qiu (2020)] to include the liquefaction triggering logic proposed by Idriss et al. [Idriss and Boulanger (2008)].
3 Large-scale numerical modeling of 3D soil-structure interaction Seismic response of structures is dictated by soil-structure interaction (SSI) during earthquake loadings. To reproduce such SSI effects computationally, a large domain of the soil surrounding the structure is often necessary to be modeled [Chae, Ugai and Wakai (2004) (2013)]. The observed response was often noted to be highly influenced by the global bridge-ground overall characteristics as an integral system.  Figure 6: A bridge-ground system [Qiu (2020)]: (a) Ground configuration; (b) Deformed mesh at the end of shaking (factor=15) On this basis, a representative 3D FE model underlain by liquefiable soil was developed to numerically investigate the consequences of liquefaction-induced ground deformation [Qiu (2020)]. The overall ground configuration is shown in Fig. 6(a). In general, the site soil profile consists of shallow, fine-grained soils underlain by medium dense to dense sands, with loose sand layers (red zone in Fig. 6(a)). Water table was prescribed at the elevation of 30.5 m. The 18-span reinforced concrete bridge is approximately 197 m long and 9.9 m wide, composed of a reinforced concrete deck on vertical pier walls with pile groups and/or single piles. Two expansion joints (with a 0.05 m gap each, Fig. 6(a)) are located adjacent to the piers. In the transverse direction, a 1.83 m wide 3D slice was taken as a representative of the bridge's geometric configuration. As such, the FE model of actual 3D geometric pile layout was presented, with the mesh composed of 17,415 brick elements. Along both side mesh boundaries, 2D plane strain soil columns of large size and depth (not shown) are included to minimize boundary effects, efficiently reproduce the desired freefield site response at these locations [Su, Lu, Elgamal et al. (2017)]. Along the FE model base, the Lysmer-Kuhlemeyer absorbing boundary [Lysmer and Kuhlemeyer (1969)] was applied. Input motion was selected as that of the 1940 Imperial Valley earthquake at El Centro Station (Component S00E), scaled to a peak acceleration of about 0.4 g. Using this excitation, via deconvolution by Shake91 [Idriss and Sun (1993)], an incident earthquake motion was derived and imparted [Elgamal, Yan, Yang et al. (2008)] in the x-direction (i.e., input shaking only in the longitudinal direction), along the base of the FE model (elevation 0.0 in Fig. 6(a)). Fig. 6(b) displays the deformed FE mesh at the end of shaking (colors show horizontal displacement). Away from the bridge, lateral deformation above the liquefied strata was to the right, about 0.2 m on the left-side and only 0.008 m on the right-side. Superposed on this global displacement pattern (to the right), local downslope deformations are seen at both ends of the bridge (Fig. 6(b)). Near the abutments, additional downslope deformation relative to the surrounding ground was about 0.35 m (left-side), and 0.4 m (right-side). Due to the downslope deformations at both ends, expansion joint gaps were closed, with the bridge deck through its piers and foundations acting as a strut that provides added restraint to the slope deformations.  (2013)]. The observed seismic performance of pile-supported wharves is significantly dictated by the behavior of the surrounding soil. To address this consideration, a ground-structure 3D FE model was developed [Su, Lu, Elgamal et al. (2017)] to capture and further elucidate many salient characteristics, such as pile-to-pile and pile-ground interaction mechanisms. On this basis, a representative 6.3 m wide 3D slice of a wharf-ground configuration ( Fig. 7(a)) was presented [EMI (2001); Su, Lu, Elgamal et al. (2017)]. As such, the FE mesh of this pile-supported wharf included 9,649 3D brick elements ( Fig.  7(a)). The input motion was selected as 1994 Northridge earthquake (Component S48W) at the Rinaldi Receiving Station, scaled down to half amplitude [Su, Lu, Elgamal et al. (2017)]. Along both side mesh boundaries, 2D plane strain soil columns of large size and depth (not shown) are included to minimize boundary effects ( Fig. 7(a)), efficiently reproduce the desired free-field site response at these locations [Su, Lu, Elgamal et al. (2017)]. A parallel version of OpenSees (OpenSeesMP) was employed to conduct the 3D simulations of this wharf-ground configuration ( Fig. 7(a)). OpenSeesMP was created for running input scripts with user-specified subdomains, parallel equation solver and parallel numberer [McKenna (2011)]. As seen in Fig. 7(a), the FE mesh of the pile-supported wharf model was divided into three zones (handled simultaneously by three different processors). With the computations conducted in parallel (implicit time integration), execution time was reduced by more than 50% [Lu, Peng, Elgamal et al. (2004)]. Fig. 7(b) displays the deformed mesh at maximum deck displacement. In this figure, the largest lateral displacements showed in the upper layers, mainly due to the shear deformation in the relatively weak clay stratum ( Fig. 7(a)). Away from the slope, lateral deformation on the landside above the clay was about 0.25 m towards the waterside. Superposed on this landside displacement pattern, local downslope deformation is seen in the slope zone ( Fig.  7(b)), reaching a peak value of about 0.45 m relative to the surrounding ground.

OpenSeesPL, ground modification, and pile/pile group analyses
A PC-based graphical pre-and post-processor (user interface) OpenSeesPL has been developed to facilitate efficient execution of 3D ground-foundation FE simulations [Lu (2006); Lu, Yang and Elgamal (2006); ]; http://www.soilquake.net/openseespl/) [Lu (2006); Lu, Yang and Elgamal (2006); ]. Currently, OpenSeesPL permits analyses of footings, piles and pile groups (Fig. 8) under static and earthquake loading conditions [Wang, Lu and Elgamal et al. (2013)]. Various scenarios of ground modification can be also studied by specifying the material within the pile zone [Elgamal, Lu and Forcellini (2009) OpenSeesPL includes options for (1) defining the pile geometry (circular or square pile) and soil material properties (linearly elastic or nonlinear), (2) defining the 3D spatial soil domain, (3) defining the boundary conditions and ground input excitation or push-over analysis parameters, and (4) selecting soil materials from an available menu of cohesive and cohesionless soil materials. The menu of soil materials includes a complementary set of soil modeling parameters representing soft, medium, and stiff cohesive soils, and loose, medium, and dense cohesionless soils (with gravel, sand, or silt permeability). OpenSeesPL also allows convenient postprocessing and graphical visualization of the FE analysis results including ground response time histories, pile responses (profiles and time histories), and the deformed mesh ( Fig. 8(b)). The user interface makes it possible for users to quickly build a FE model, run the FE analysis, and evaluate the performance of the groundfoundation system [Lu, Yang and Elgamal (2006)]. Representative studies employing OpenSeesPL are presented below.

Stone column mitigation of liquefaction-induced lateral ground deformations
Liquefaction-induced lateral spreading may result in disruption of function, significant damage, and considerable replacement cost for existing foundations, buildings, and structures. This type of adverse response has been observed in many earthquakes, such as Niigata, Japan 1964 [Seed and Idriss (1967)], Chi-Chi, Taiwan 1999[EERI (2000], Kocaeli, Turkey 1999[EERI (2001] and Christchurch, New Zealand 2011 [Cubrinovski, Bradley, Wotherspoon et al. (2011);Cubrinovski (2013)]. The risk of liquefaction and associated ground deformation may be reduced by ground improvement approaches including solidification (cementation), densification, and gravel drains or stone columns (SCs). Use of gravel drains/SCs is a relatively recent development compared with more traditional soil remediation approaches [Elgamal, Lu and Forcellini (2009)].

Mitigation by stone columns
Using OpenSeesPL, a mildly inclined 4° saturated sand layer of 10 m thickness was studied (permeability k=6.6×10 −5 m/s), based on Nevada sand properties at a relative density Dr of about 40% ; Elgamal, Yang, Parra et al. (2003); Elgamal, Lu and Forcellini (2009)]. Dense sand properties with a representative gravel permeability k=1.0×10 −2 m/s were employed to represent the SC properties. In the remediated cases, the area replacement ratio = = 2 4 2 is conventionally defined as the ratio of the SC area Ar over the tributary area A ( Fig. 9(a)), where d=diameter of the SC, and s=distance (spacing) between SC centers. Diameter of the SC in all remediated cases was kept at 0.6 m. A total of five simulations were conducted to explore Arr values from 10% up to 50% with base input motion being that recorded at the Wildlife site during the 1987 Superstition Hills earthquake [Elgamal, Lu and Forcellini (2009)]. The medium sand case (MS) represents the original benchmark unremediated scenario. As such, the MS results provide a reference of free-field site response for comparing with the remaining SC remediated cases ( Fig. 9(b)). When the SC permeability k=1.0×10 -2 m/s, the permanent displacement at the ground surface fell below the large displacement range (0.3 m or 1 ft according to Youd 1989 for Arr=20%-30%) [Youd (1989)]. To achieve a lower permanent displacement [CGS (2008)], an Arr=40% is required.

Influence of stone column spatial distribution
In order to explore influence of SC spatial distribution on liquefaction-induced lateral deformation, a number of SC configurations (1×1, 2×2 and 4×4) for a given Arr=10% were studied ( Fig. 10(a)). In these configurations [Lu, Kamatchi and Elgamal (2019)], the same Arr is maintained by deploying smaller diameter columns, evenly distributed spatially, so as to increasingly reduce the drainage path ( Fig. 10(b)). On this basis, a total of seven simulations were performed exploring different combinations of SC geometric configurations and permeabilities. As mentioned above, the MS (unremediated) case represents the original benchmark situation to provide a reference of free-field site response for comparing with the remaining remediated cases. Figs. 10(c) and 10(d) show the lateral displacement time histories at the ground surface. Permanent ground surface displacement for the unremediated case is 1.15 m. A major reduction in ground deformation resulted because of the increased SC permeability. The more evenly distributed (or smaller-diameter) SC configuration (e.g., the 2×2 pattern compared to the 1×1, or the 4×4 compared to the 2×2) is seen to be more effective in reducing the lateral deformation of the sand stratum (Figs. 10(c) and 10(d)). As mentioned above, this is mainly attributed to the shorter drainage path in the smaller-diameter configurations. With a higher SC permeability k=0.1 m/s, permanent ground surface deformation falls below the large displacement range (0.3 or 1 ft m according to Youd 1989) [Youd (1989)], with substantial reductions in the 2×2 and 4×4 configurations compared to the 1×1 scenario [Lu, Kamatchi and Elgamal (2019)].

Lateral load on a large pile group
Using OpenSeesPL, a systematic 3D FE study of a large pile group subjected to lateral loading was performed (Figs. 11(a) and 11(b)). The FE model was representative of the Dumbarton Bridge (California) Pier 23 pile-group foundation geometry [Wang, Lu and Elgamal et al. (2013); Wang (2015)]. The pile group is configured in an 8×4 arrangement with a longitudinal spacing of 2 pile diameters and a transversal spacing of 2.15 pile diameters on center. The diameter and length for each pile were 1.37 m and 30.8 m, respectively. The group is rigidly connected by a pile cap, at a substantial distance of 14.3 m above the mudline. (c) Figure 11: A large pile group subjected to lateral loading [Wang, Lu and Elgamal et al. (2013)]: (a) Close-up of FE mesh (created in OpenSeesPL); (b) Layout of pile group; (c) Axial forces, initial, and at 0.3 m pile cap deflection The pile response was assumed to be linear elastic and the bending stiffness for each pile was EI=2×10 6 kN-m 2 [Wang, Lu and Elgamal (2013)]. A summary of the site soil profile is reported in Wang (2015) [Wang (2015)]. In the employed OpenSeesPL half mesh of Fig.  11 (due to symmetry), the vertical dead load (28,900 kN representing the tributary own weight of the bridge deck.) was imposed initially (after application of soil own weight). A lateral pile cap longitudinal displacement was then applied (at the center of the pile cap) up to a maximum of 0.30 m. Fig. 11(c) shows the axial force distribution between piles (1-16). In this figure, a wide range in the share of load is seen in each pile. At the initial static state, the corner piles (1, 8, 9 and 16) carry the highest force. At the 0.3 m pile cap displacement, the compressive axial forces of front piles (6-8 and 14-16) increase significantly, conversely the back piles (1, 2 and 9, 10) experience tensile forces. Such variation in axial forces is influenced by the relatively large distance between the pile cap and the mudline, exerting a significant moment on the overall pile foundation along with the pile lateral forces.

Performance-based earthquake engineering (PBEE) for bridge systems and
MSBridge PBEE aims to quantify the seismic performance of engineered facilities, such as bridges, using metrics that are compatible with network risk assessment and of interest to engineers and stakeholders alike [Mackie, Lu and Elgamal (2012)]. Using a fully developed bridge PBEE assessment framework [Mackie, Lu and Elgamal (2012)], a pre-and post-processor user interface (Fig. 12) was developed for application within the simple environment of a 2 span, single column bridge founded on a 3D ground soil mesh. This interface (https://peer.berkeley.edu/bridgepbee/) integrates all elements of the underlying FE simulations along with the subsequent PBEE calculations [Mackie, Lu and Elgamal (2012)]. In the time history analysis phase, the bridge column is modeled by nonlinear fiber section beam-column elements. Abutment models considering passive soil resistance, bearing pads, and shear keys are included. If a computation cluster is available, a parallel version of OpenSees (OpenSeesSP) can be used to efficiently complete all individual shaking event simulations simultaneously (Fig. 13), in a couple of hours versus a couple of days for instance.   The bridge-foundation-ground mesh is shown in Fig. 14(a) for illustration. Response of the same bridge is explored under stiff and soft supporting soil conditions [Mackie, Lu and Elgamal (2012)]. Fig. 14(b) shows the corresponding mean repair cost ratio (RCR) loss models. The consequences of shaking and repair do not begin to accumulate until a PGV of approximately 20 cm/s is reached, and repair costs are clearly higher for the soft soil scenario. Building on the above developments, systematic PBEE analyses were recently implemented for multi-span reinforced concrete bridges [Almutairi, Lu, Elgamal et al. (2018)]. In this regard, the bridge seismic response is handled by the new highly versatile pre-and post-processor graphical user interface MSBridge [Elgamal, Lu, Almutairi et al. (2017); Almutairi, Lu, Elgamal et al. (2018)], http://www.soilquake.net/msbridge/) [Elgamal, Lu, Almutairi et al. (2017); Almutairi, Lu, Elgamal et al. (2018)]. Using OpenSees, this interface is focused on efficiently conducting nonlinear time history analyses for a wide range of multi-span bridge configurations. In this tool ( Fig. 15(a)), bridge structures, abutments, and foundation response mechanisms are integrated within a unified framework. Furthermore, MSBridge allows for addressing possible variability in the bridge deck, bent cap, column, foundation, and/or ground configuration properties (on a bent-by-bent basis). As such, MSBridge permits the simulation of key scenarios of significance for bridge upgrades, widening, extensions, and retrofits. Convenient postprocessing capabilities allow graphical visualization of the results including the deformed mesh ( Fig. 15(b)), as well as column/pile response. Using MSBridge, it is possible for structural engineers/researchers to rapidly build models of complex multi-span bridges, run the FE analysis, evaluate performance of the bridge-ground system, and assess the PBEE outcomes [Elgamal, Lu and Mackie (2014); Almutairi, Lu, Elgamal et al. (2018)].
(a) (b) Figure 14: PBEE analysis of a bridge-ground system [Mackie, Lu and Elgamal (2012)]: (a) FE mesh (created in BridgePBEE); (b) Mean repair cost ratios (case 1 represents stiff ground, and Case 2 represents soft ground) In summary, the main capabilities of MSBridge include: i) horizontal and vertical alignments, with different skew angles for bents/ abutments; ii) beam-column element with fiber section for bridge columns and piles; iii) deck hinges, isolation bearings, steel jackets, and abutment models; iv) foundation represented by Foundation Matrix or soil springs (py, t-z and q-z); v) analysis for suites of ground motions (built-in and/or user-defined); vi) PBEE outcomes in terms of repair cost, time and carbon footprint; and vii) pushover and mode shape analysis options.

Summary and conclusions
Recent research efforts that address the static and dynamic/seismic response of groundfoundation-structure systems have been presented. Specifically, developments within the computational platform OpenSees were outlined. Representative numerical results were included for the response of large foundation-ground-structure systems, ground modification liquefaction countermeasures, and PBEE analyses. To facilitate 3D computations, graphical user interfaces (OpenSeesPL, BridgePBEE and MSBridge) were shown to be useful tools for conducting simulations of ground modification, ground-pile SSI, and multi-span bridge-ground systems with PBEE assessments. Overall, the presented studies aim to illustrate the potential for further reliance on computer simulation in the assessment of nonlinear SSI response. Challenges in calibration and in high fidelity modeling are being gradually overcome. With careful attention to the involved modeling details, effective insights may be gleaned for a wide range of practical applications.