Validation of OpenSees for Tsunami Loading on Bridge Superstructures

The performance of bridges prone to tsunami hazards is critical to the resilience of coastal communities. An important component of resilience is the development of design guidelines, which requires simulating the hydrodynamic loads imparted on bridges by tsunami waves. To support such simulation efforts, data collected from wave flume experiments were used to validate the particle finite-element method (PFEM) implementation in OpenSees, which has been widely used for earthquake engineering simulation of bridges. The numerical results agree well with the experimental data, showing the potential use ofOpenSees for determining hydrodynamic loads on engineered structures. Given the wide-ranging earthquake engineering simulation capabilities ofOpenSees, these validations form a basis for simulating bridge responses to multihazard earthquake and tsunami scenarios. DOI: 10.1061/(ASCE)BE.1943-5592.0001221. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, http://creativecommons.org/licenses/by/4.0/. Author keywords: Simulation; Tsunamis; Wave loads; Finite elements; Particle methods.


Introduction
In recent years, major tsunamis have taken many lives and caused severe damage to coastal infrastructure. The 2004 tsunami in the Indian Ocean killed 230,000 people, and the 2011 Tohoku tsunami in the Pacific Ocean killed more than 30,000 people and destroyed more than 125,000 buildings. In both of these events, earthquake and tsunami damage was observed for structures of various construction types, including steel, reinforced concrete, and timber Saatcioglu et al. 2006;Fraser et al. 2013). In addition to building damage, the Tohoku tsunami washed away more than 300 bridges (Suppasri et al. 2014). The loss of functionality of these bridges severely hindered emergency response and post-event recovery. Aside from the loss of human lives, the Tohoku event was a stark reminder that areas with well-engineered infrastructure designed to resist high seismic loads are not immune to the hydrodynamic loads that can be generated by tsunamis.
The last major tsunami to affect the United States resulted from the magnitude 9.2 Alaskan earthquake in 1964. In addition to damage in Anchorage, tsunami waves reached the Pacific Northwest and northern California, leading to damage and casualties. Hawaii also faces tsunami hazards from events around the Pacific rim; for example, a magnitude 8.3 earthquake in Chile in 1960 caused a tsunami that led to severe damage in the town of Hilo. Coastal communities in the Pacific Northwest face tsunami hazards from the Cascadia Subduction Zone (CSZ), which has a 40% chance of a major earthquake and subsequent tsunami occurring in the next 50 years (Goldfinger et al. 2012). Many bridges constructed during the 1950s and 1960s, with little to no seismic provisions and thus very little resistance to lateral loads, are still in use today along the Oregon and Washington coasts. Some of these bridges serve as the only regional transportation lifeline for coastal communities. Determining tsunami loads is critical to the design and assessment of these, and other, coastal bridges.
Various physical experiments have been conducted to assess tsunami loading on bridges. In response to the 2011 Tohoku event, the Public Works Research Institute (PWRI) in Japan conducted wave flume experiments to determine hydrodynamic loads imparted on bridge superstructures (Hoshikuma et al. 2013). In addition, Hayatdavoodi et al. (2014) and Seiffert et al. (2015) performed experiments and computational simulations for solitary wave loads on deck-girder bridges. The results indicated that trapped air and the number of girders have a significant effect on the forces imparted on elevated deck-girder superstructures. Istrati et al. (2016) conducted experiments to investigate forces on bridge-deck-to-substructure connections during tsunami loading, where it was found that flexible connections can lead to significant reductions in superstructure forces.
With experimental data, numerical simulation models of bridge response to tsunami loading can be validated to explore mitigation strategies and guide the design of additional experiments, as well as to develop guidelines for tsunami-resistant bridge design. Analysis of bridge deck failure resulting from tsunami loading was conducted by Bricker et al. (2012) using Reynolds-averaged Navier-Stokes (RANS) models in the OpenFOAM computational fluid dynamics program. Motley et al. (2016) also used OpenFOAM to investigate the effects of bridge skew on the hydrodynamic forces imparted by tsunami bores. A systematic integration of numerical models for tsunami generation, propagation, and inundation was developed by Azadbakht and Yim (2016), who used site-specific flow height and velocity to compute bridge forces using LS-DYNA. Numerical studies of tsunami loading on bridge piers were carried out by Wei et al. (2015) using smoothed particle hydrodynamics (SPH). In addition, SPH was used to investigate tsunami loads on bridge superstructures with and without protective breakwaters (Wei and Dalrymple 2016).
Simulation of the sequential effects of seismic loading followed by hydrodynamic loading motivates the continued development and validation of OpenSees for fluid-structure interaction (FSI) simulations. OpenSees was developed for the simulation of structural and geotechnical system response to earthquake loading, and it has been used extensively for earthquake engineering simulations of bridge structures (Aygün et al. 2011;Seo and Linzell 2013;Liang et al. 2016). Recently, OpenSees was extended to the particle finite-element method (PFEM) (Zhu and Scott 2014b), a Lagrangian formulation that is able to simulate FSI with the wide range of nonlinear finite-element formulations and constitutive models in OpenSees.
In this article, FSI simulations using the PFEM in OpenSees are validated against wave flume experiments conducted at PWRI on scaled models of various bridge deck configurations (Hoshikuma et al. 2013). After a brief overview of the PFEM approach to FSI, the PWRI wave flume experiments are described. Wave height and velocity in the flume, along with pressures and forces exerted on the bridge decks, are compared with the experimental results. The validated FSI simulations make OpenSees available for a number of wide-ranging applications involving sequential earthquake and tsunami hazards.

FSI Simulation Using the PFEM
Various simulation approaches have been developed for FSI with incompressible Newtonian flows using either Lagrangian, Eulerian, or arbitrary Lagrangian-Eulerian (ALE) formulations (Girault and Raviart 1986;Gunzburger 1989;Baiges and Codina 2010;Radovitzky and Ortiz 1998;Tezduyar et al. 1992). All three formulations are suited to particular classes of problems for the simulation of solids, fluids, FSI, and contact. Important advantages of the Lagrangian formulation with respect to FSI include easy tracking of the fluid free surface and conformity with the Lagrangian formulation of structural mechanics. The PFEM (Oñate et al. 2004;Carbonell et al. 2010;Cremonesi et al. 2011a, b;Oñate et al. 2010;Sabel et al. 2014) is one approach that uses a Lagrangian formulation of fluid response. The motion of fluid and structural particles is tracked at all time steps, and the fluid mesh is updated at every time step to maintain numerical accuracy. With Lagrangian formulations of the fluid and solid response, the PFEM equations of FSI are monolithic, leading to a more efficient implementation than staggered approaches with disparate modules for fluid and structure response.
Initial implementations of the PFEM used a moving-mesh formulation where all meshes of the fluid and structure domains are updated at each time step (Oñate et al. 2004). This movingmesh approach makes it easy to track the fluid free surface and to have a high-quality fluid mesh at the start of each time step. However, the mesh is not updated during iterations within the time step, which can require small time steps to prevent illconditioning of the global system of equations. A further result of the ill-conditioning is a very significant oscillation in the computed pressures and forces on the structure, making it difficult to quantify engineering demand parameters (Cressman 2016). To overcome these numerical difficulties, subsequent versions of the PFEM used a fully fixed mesh with Lagrangian particles (Idelsohn et al. 2014), similar to particle-in-cell methods (Brackbill and Ruppel 1986;Brackbill et al. 1988) and the material point method (Więckowski 2004). The fixed-mesh approach is more computationally efficient than the moving-mesh approach because remeshing is not required. This approach has been used for multifluid flows and FSI (Becker et al. 2015); however, it is limited to specific types of solid elements not typically employed in structural engineering applications.
The approach taken herein to validate OpenSees PFEM simulations of tsunami loading on bridges is the background mesh method, which strikes a balance between the moving-mesh and fixed-mesh formulations (Becker et al. 2015). In this approach, the fluid mesh is fixed, except at the interface of the fluid and structure domains, where a local moving mesh is created at each time step. The structural mesh is not constrained to that of the fixed fluid mesh, and as a result, any line or solid finite element based on a Lagrangian formulation can be used for FSI simulations. Complete details of the PFEM fixed-mesh method are described by Becker et al. (2015); however, a brief overview of the computational steps taken during one time step are described as follows: 1. The initial states of a set of Lagrangian fluid particles (dots) are shown in Fig. 1(a), along with a structural mesh (shaded area) that does not change during the FSI simulation. The location and number of particles are specified by the analyst. The FSI area is identified as the white area surrounding the structure. 2. The fixed fluid mesh is created except in the FSI area. Only elements containing particles are created. Both nodes and elements in the fixed mesh are uniformly distributed and not distorted, as shown in Fig. 1(b), which leads to good computational conditions. 3. Delaunay triangulation (Shewchuk 1997) is performed for nodes on the structural exterior and the surrounding nodes of the fixed mesh. Interface fluid elements are created if they connect at least one grid node [ Fig. 1(c)]. 4. Based on the newly created mesh, the governing PFEM equations for velocity and pressure are formed and solved (Zhu and Scott 2014a). The state of each particle is interpolated from the displacements, velocities, and pressures of its containing element. Then particles are moved to their new positions [ Fig. 1(d)], and the analysis proceeds to the next time step.
With this computational approach, the entire mesh is neither fully fixed nor fully moving because it is only updated in the local area around the structural domain. As a result, the computational cost is significantly reduced compared with the full moving-mesh method. In addition, this approach gives flexibility of the mesh updating in the FSI area because structural nodes are not always on the grid points.

Description of PWRI Experiments
To address bridge failures observed in the 2011 east Japan earthquake and tsunami, PWRI performed a series of experiments in a 1-m-wide, 30-m-long wave flume, as shown in Fig. 2. A gate served as a dam to a 12-m-long reservoir filled to height h res , and bridge deck models were placed 7.5 m downstream from the gate on a pier 20 cm above the base of the flume. Two wave gauges located 1 and 2.5 m in front of the bridge deck recorded the flow height and velocity of the approaching wave.
The test cases consisted of various bridge deck configurations, standing water levels, and wave heights for which hydrodynamic forces imparted on the bridge decks were recorded (Hoshikuma et al. 2013). The three cases shown in Fig. 3 were selected for validation of the PFEM background mesh implementation in OpenSees. The model bridge deck was developed at 1/20 length scale, with cross sectional dimensions of 50-cm width, 10-cm height, and 5-cm overhang. These dimensions correspond to a prototype bridge deck of 10-m width, 2-m height, and 1-m overhang. Using Froude similitude laws with the 1/20 length scale, the experimental results for time and force were scaled by ffiffiffiffiffi 20 p and 20 3 , respectively, to the prototype bridge. The transverse dimension of the model bridge deck was 0.985 m.
As shown in Fig. 3, the standing water level, h, was 10 cm (2 m at prototype scale), and the clearance was 10 cm (2 m prototype) for Case 1. The reservoir depth, h res , was 61.7 cm, corresponding to a wave that would be expected to reach the top of the deck. Case 2 differed from Case 1 in its standing water level of 15 cm (3 m prototype), clearance of 5 cm (1 m prototype), and reservoir depth of 51.8 cm, which was also expected lead to the wave reaching the top of the deck. All conditions were the same for Case 3 as in Case 1, except for the addition of a protective fairing on the seaward side of the bridge. It was anticipated that the fairing would divert flow from the deck, thereby reducing hydrodynamic loads on the bridge deck.

Numerical Simulation Results
Numerical simulations were carried out in OpenSees for the three cases shown in Fig. 3. In the experiments, the dam gate was   Fig. 4, were compared with recorded values. Finally, the simulated total horizontal and vertical force were compared to the forces recorded by biaxial load cells at the four girder-to-pier connections. The simulated forces were obtained by summing the horizontal and vertical reactions, including inertial effects, at the connection locations.  Additional details on the experimental setup and recordings can be found in Hoshikuma et al. (2013).
All three cases used a 0.02-m mesh size for the background grid, in which three or four particles were placed in each cell. Physical constants for the fluid were density r = 1,000.0 kg/m 3 and viscosity m = 10 -3 Pa-s. To represent the deck geometry, the bridge was discretized using 158 structural nodes with 204 triangular plane-stress finite elements of thickness 0.985 m. The bridge deck material had density r = 600.0 kg/m 3 , elastic modulus E = 10 GPa, and Poisson's ratio = 0.3. Because the exact locations of the pressure sensors in Fig. 4 were unknown, nodal pressures were averaged in the area of the approximate sensor location. Horizontal and vertical forces were obtained from the simulation as the dynamic reaction forces at the girder-to-pier connections.
The PFEM simulations in OpenSees use backward Euler time integration with a time step of 0.006 s. Convergence in a time step is achieved when the norms of both the incremental velocity and incremental pressure vectors are less than 10 -6 . A single-core Intel Xeon central processing unit (E3-1225 v3) at 3.20 GHz was used for the simulation, and the real time for a simulation time step averaged 4.7 sec. All simulations were performed in two dimensions with single-phase flow, modeling assumptions that will be addressed in the ensuing discussions of simulation results.
Case 1: 10-cm Standing Water, 20-cm Wave Height To simulate the wave flow and bridge response for Case 1 with a 0.02-m mesh size, the fluid domain was comprised of 24,352 grid nodes, 24,352 pressure nodes, and 45,668 fluid elements. The tank contained 1,140 Â 58 particles, and the standing water contained 1,700 Â 9 particles.
Comparisons of the simulated and experimental wave height and velocity at the 1-m wave gauge are shown in Fig. 5. There is generally good agreement of the wave height and velocity even though the mechanisms of reservoir release differ between the simulation and experiment and the simulated velocity was obtained by interpolation of surrounding particle velocities. Snapshots of wave impact and surge on the bridge deck for Case 1 are shown in Fig. 6 for both the numerical model and the experiment. The fluid pressures recorded on the seaward side of the bridge deck generally match the experimental results, as shown in Fig. 7; however, the single-phase simulation was not able to capture the initial pressure shocks recorded in the experiment. In addition, compared with the experiments shown in Fig. 6, there appeared to be a larger volume of overtopping water in the simulation, which also contributed to a larger simulated pressure at location p 2 .
The simulated total horizontal and vertical forces are compared with the experiment in Fig. 8. Although there is good agreement of the peak horizontal force and peak downward vertical force, the simulation was not able to capture the peak vertical uplift force. Significant differences between simulated and experimental peak uplift force were also noted by Motley et al. (2016) and Wei and Dalrymple (2016) in simulations of this same PWRI experiment. The three-dimensional (3D), single-phase smooth particle hydrodynamics (SPH) simulation by Wei and Dalrymple (2016) underestimated the maximum uplift force despite good agreement with the first peak when using a smaller particle size. Motley et al. (2016) overestimated the peak uplift force in two-dimensional (2D), two-phase simulations but was able to obtain a better prediction using 3D, two-phase flow simulations with a very fine mesh in OpenFOAM. These results reveal the importance of 3D simulations and fine resolution to capture localized flows and two-phase flow to capture the effects of trapped air.  Although there were discrepancies in the steady-state response, the OpenSees PFEM simulations captured the basic trends of the horizontal and vertical force steady-state response. The magnitude of the vertical force was significantly affected by the exclusion of trapped air with single-phase flow (Bricker and Nakayama 2014;Seiffert et al. 2015). Compared with the experiment shown in Fig. 8, there was a large simulated downward vertical force because trapped air was not present to counter the weight of the overtopping water, the volume of which was also overpredicted in the simulation. In addition, the lack of fluid flow behind the bridge contributed to the simulation's overprediction of the steady-state horizontal force.
A quantitative comparison using normalized root-mean-square deviation (RMSD) of the simulated and experimental time histories of total horizontal and vertical forces is shown in Table 1. The RMSD was less than 20% for both force-time histories in Case 1.

Case 2: 15-cm Standing Water, 15-cm Wave Height
The experimental model for Case 2 differed from Case 1 in the standing water level and the wave height, as indicated in Fig. 3. As in Case 1, the wave was expected to reach the top of the bridge deck (30-cm model scale); however, the wave would be slower in this case. For the numerical simulation, the reservoir contained 1,068 Â 46 particles, and the standing water contained 1,593 Â 14 particles from a 0.02-m mesh size. The total number of particles was slightly less than in Case 1.
The simulated wave height and velocity are in good agreement with the experimental results, as shown in Fig. 9. Snapshots of the experiment are not available for this case. Comparisons of the simulated and experimentally recorded pressures on the seaward side of the bridge deck are shown in Fig. 10. The steadystate pressure generally agrees well except for the location on top of the bridge deck (p 2 ), which could be attributed to differences  in the amount of overtopping water between the experiment and simulation.
The total horizontal and vertical forces for Case 2 are compared for the numerical and experimental results in Fig. 11. The simulation underestimated the peak horizontal force but showed less discrepancy with the steady-state horizontal force than that observed in Case 1. Similar to Case 1, the simulation underpredicted the peak uplift force; however, unlike Case 1, the steadystate vertical downward force shows generally good agreement. These observations are corroborated by lower RMSD values (compared with Case 1) for the force-time histories, as shown in Table 1. The experimental setup for Case 3 was the same as in Case 1 except for the addition of a triangular fairing to divert wave impact from the bridge deck. With the same wave height and velocity as Case 1, the experimental and simulated results can be compared to evaluate the effectiveness of the fairing in reducing the forces on the bridge deck. The same mesh size of 0.02 m and fluid and bridge materials were also used. The tank contained 1,068 Â 55 particles, and the standing water contained 1,593 Â 8 particles.
The time histories of wave height and velocity were similar to those of Case 1, as expected (Fig. 12), but with some slight discrepancies after wave impact. Whereas the simulated flow conditions were similar between Cases 1 and 3, there were differences in the experimentally recorded flow conditions, emphasizing the difficulty of repeating tsunami bore experiments. Snapshots of the wave impact and surge on the bridge deck for Case 3 are shown in Fig. 13 along with the experimental photos. Pressures at the four locations shown in Fig. 4 were not recorded during the experiment because of the addition of the fairing, which covered the sensor locations shown in Fig. 4.
The effectiveness of the fairing was quantified by comparing the total forces for Case 3 (Fig. 14) with those from Case 1 (Fig. 8), where the peak horizontal force was reduced by 25% and the peak uplift force by 50%. The peak horizontal force was reduced because the fairing diverted fluid away from oblique contact with the bridge deck, and the peak uplift force was also reduced. Both the numerical and experimental results reflect these reductions in horizontal and vertical force. The total simulated downward force and the experimental downward force at the steady state were roughly the same for the two cases because the fairing did not significantly affect the amount of water flowing over the deck or the amount of trapped air.
The numerical results generally match the experiments; however, the downward forces were overestimated by the simulation for Case 3. As shown in Table 1, the deviation of vertical forces was more than 30% from the experimental results. A reasonable explanation is found in Seiffert et al. (2015) and Bricker and Nakayama (2014) that the uplift forces with fully trapped air can be as much as six times larger than the case without trapped air. Although the simulated wave height and velocity were similar between Case 1 and 3, there was natural variation in the experimentally recorded flow conditions for these cases that the simulations could not capture.

Conclusion
Validations of the PFEM in OpenSees for FSI simulations of tsunami loading on bridge decks were presented using data from wave  flume experiments conducted at PWRI in Japan. Comparisons of wave height and velocity in the flume show that the OpenSees PFEM implementation is able to simulate tsunami bore flow as it approaches the bridge. In addition, the OpenSees simulations are able to capture pressure at various locations on the bridge deck. The simulations reproduce the peak horizontal forces acting on the bridge decks, including the case where fairing protects the seaward side of the bridge. The simulations, however, tended to underestimate the peak uplift force and overestimate the steady state downward force. These results support the importance of 3D simulation and two-phase flow reported in previous research on simulations of tsunami-bridge interaction.
In future developments of the PFEM in OpenSees, 3D and two-phase simulations will be carried out to capture pressure shocks and the localized air pockets and flows observed between the girders of the bridge superstructure. In addition, for computational efficiency, an approach using a variable mesh size, where a coarse mesh is used for the fluid domain in conjunction with a refined mesh in and around the structural domain, will be addressed in OpenSees. These outcomes have important implications for the development of guidelines and loading equations for tsunami-resistant bridge design and the design of additional experiments to investigate protective strategies. Finally, validated FSI models in OpenSees open this software framework to additional applications for sequential earthquake and tsunami hazards (Carey et al. 2014;Kim et al. 2017;Scott and Mason 2017).