Multi-Scale Bridge Wash Out Simulation During Tsunami by Using a Particle Method

In 2011, the huge tsunami caused by the great east Japan earthquake devastated many infrastructures in pacific coast of north eastern Japan. Particularly, collapse of bridges caused a traffic disorder and these collapse behaviors led to delay of recovery after the disaster. In this study, the bridge wash away accident is selected as a target issue, and it is represented by a numerical simulation. For this purpose, Smoothed Particle Hydrodynamics (SPH) Method, which is one of the pure mesh free methods, is utilized for the rigid body motion simulation. In this study, rigid body motion is introduced for the fluid-rigid interaction behavior during bridge wash away simulation. In the numerical analysis, the upper bridge structure is washed away by receiving an impact fluid force. The wash away simulation of two types of the bridge girder showed good agreement with the real accident on the great east Japan earthquake tsunami.


Introduction
On March 11, 2011, the huge tsunami caused by the great east Japan earthquake devastated many infrastructures in pacific coast of north eastern Japan.Particularly, the damage of outflow of bridge girders caused a traffic disorder and these collapse behaviors led to delay of recovery after the disaster.After 2011 tsunami, disaster prevention and mitigation techniques are actively developing in coastal infrastructures and establishing prediction method for tsunami disaster is one of the severe issues toward the next millennium tsunami.
In this study, the bridge wash out accident is selected as a target issue, and we try to represent these accidents by using a numerical analysis.For this purpose, one of the mesh free methods; Smoothed Particle Hydrodynamics (SPH) Method [1,2] is utilized for Tsunami flow.The main advantage of SPH and the other particle simulation is the absence of a computational grid or mesh since it is spatially discretized into Lagrangian moving particles.This allows the possibility of easily modeling flows with a complex geometry or flows where large deformations or the appearance of a free surface occurs.Recently, this method is widely used in field of fluid and solid dynamics.A stabilized ISPH [3], which is one of the modified versions of the SPH and can evaluate much smoother pressure distribution, has been developed by our research group.Then a fluid-solid interaction algorithm including rigid body motion is developed in this study.

SPH Formulation
In this paper, SPH method was adopted.A basic concept in SPH method is that any function φ attached to particle "i" at a position r i is written as a summation of contributions from neighbor particles Note that, the triangle bracket <ϕ i > means SPH approximation of a function φ.The divergence of a vector function can be assumed by using the above defined SPH approximation as follows: and the expression for the gradient can be represented by: In this study, a stabilized ISPH method developed in the incompressible fluid analysis was adopted.In this method, the pressure is calculated implicitly and the velocity is calculated explicitly.

Governing equations
The governing equations are the continuum equation and the Navier-Stokes equation.These equations for the flow are represented as: wthere, ρ and ν are density and kinematic viscosity of fluid, u and p are the velocity and pressure vectors of fluid respectively.F is external force, and t indicates time.The turbulence stress τ is necessary to represent the effects of turbulence with coarse spatial grids.In the most general incompressible flow approach, the density is assumed by a constant value with its initial value.

Modification in the source term of pressure Poisson equation
The main concept in an incompressible SPH method is to solve a discretized pressure Poisson equation at every time step to get the pressure value.In a sense of physical observation, physical density should keep its initial value for incompressible flow.However, during numerical simulation, the 'particle' density may change slightly from the initial value because the particle density is strongly dependent on particle locations in the SPH method.If the particle distribution can keep almost uniformity, the difference between 'physical' and 'particle' density may be vanishingly small.In other words, accurate SPH results in incompressible flow need to keep the uniform particle distribution.For this purpose, the different source term in pressure Poisson equation can be derived using the 'particle' density.The SPH interpolations are introduced into the original mass conservation law before the perfect compressibility condition is applied.
Then, the pressure Poisson equation reformulated as: where, α is relaxation coefficient, u i * is temporal velocity and triangle bracket < > means SPH approximation.Note that this relaxation coefficient is strongly dependent on the time increment and the particle resolution.Then, the reasonable value can be estimated by the simple hydrostatic pressure test using the same settings on its time increment and the resolution.

Treatment of moving rigid body
In this study, the general momentum conservation law of the rigid body is solved numerically with the external forces including hydrodynamic as a fluid-rigid interaction formulation.Figure 1 shows the rigid motion algorithm.The rigid body is discretized to particles at the beginning.Next, it is postulated that the external force is calculated from the gravity g, hydrodynamic force at the rigid surface F f and contact force between rigid body and fixed boundary F e , and these value are calculated by following equations.
where, P i is the pressure of surface rigid particle, and ΔS i is the area of surface rigid particle.And the contact force is modeled by the elastic contact theory by Hertz.In the above equation, l is the distance between rigid particle and the fixed boundary.In above equation, ν i , E i , E j is Poisson ratio, elastic modulus of the particle coming into contact, elastic modulus of the contacted particle.And then the hydrodynamic moment M f and contact moment M e are calculated by following equation.And then, translational velocity T and angular velocity vector ω are evaluated by following equations.

( )
Finally, the velocity of each particle in the solid body can be expressed as :

Bridge Wash Out Simulation
Figure 2 shows Utatsu Bridge after the Tsunami impact, the bridge is in Minamisanriku-cho, Miyagi prefecture.From this picture, high height girders turned over and low height girders are pushed out without rotation.In this study, we have selected two types of bridge girder model; high girder and low girder model as a target model and compare the simulation results and disaster report.The detail of analysis model and girder model are shown in Figure 3 to Figure 5.The density of particles of the girder model is 24.5 kg/m 3 referring to the density of reinforced concrete.And the wave is modeled for a gentle stream.Initial water level is set to be 8m and the initial velocity of the wave is set 10m/s.In addition, 10m/s is continuously given at the position of 20m from the left corner of the water storage.The depth in the models is 11m.The particle distance d 0 = 0.2m, time increment Δt = 0.001s and the total number of particles is about one million.
The results of bridge wash out behavior are shown in Figure 6.From these results, the high girder model case shows a higher rotation than the low girder model case.The disaster report say that the high height girders turned over and low girders just wash out without rotation, and our simulation results show a similar tendency to this disaster report.Our fluid-rigid coupling algorithm can expressed the change of wash out behavior depending on the shape of the bridge girder.However, turn-over of the bridge girder cannot be represented because of lack of several numerical settings.One of the causes is considered that the inflow condition is not correct as a real tsunami accident.Then, we have conducted a tsunami innundation simulation as a next topics.

Tsunami Run-up Simulation
According to the conclusion at previsou section, the tsunami run-up simulation is conducted to reconsider the inflow condition.The Utatsu area where Utatsu Bridge is in is selected as a simulation area.Figure 7 shows the analysis area and Utatsu Bridge.In addition, the Tsunami area is modeled and initial water level is set to be 10m.And then, 10m/s is continuously given for water particles which pass through the inflow area.The particle distance d 0 = 3.0m, time increment Δt = 0.001s and the total number of particles is about seven million.
Figure 8 shows the comparison of tsunami inundation area between disaster report and the simulation result.Although the simulation result shows a similar tendency in the comparison of inundation area, the details of inundation area shows several mismatch between them.The simulated result is overestimated into the northern direction, the other hand the inundation area into the western direction is smaller than disaster report [4].It's because the inflow velocity into northern direction is faster than actual velocity in this simulation.Figure 9 shows velocity distribution around the Utatsu Bridge by our simulation.In comparison in these two areas, the Tsunami velocity around the high height girder (G8-G10) is larger than the velocity around the low height girder (G3-G7).The above simple inflow conditon may be replaced by a consistent inflow condition based on a large area tsunami propagation simulation based on the shallow water equation and/or observation data.

Conclusion
In this study, fluid-rigid interaction formulation is conducted by introducing rigid motion algorithm into ISPH.From bridge wash out simulation, we could express the change of wash out behavior depending on the shape of the bridge girder and it shows good agreements with the disaster report on

Figure 7 .
Figure 7. Analysis area and Tsunami are at Utatsu.

Figure 8 .
Figure 8.Comparison of tsunami inundation between disaster report and simulation result.

Figure 9 .
Figure 9.The velocity distribution around the Utatsu Bridge