THREE-DIMENSIONAL NUMERICAL ANALYSIS ON DEFORMATION OF RUN-UP TSUNAMI AND TSUNAMI FORCE ACTING ON SQUARE STRUCTURES

A three-dimensional two-way coupled fluid-sediment interaction model (FSM) is applied to investigate run-up tsunami deformation and tsunami force acting on square structures on land. The FSM consists of a generalized Navier-Stokes solver (GNS) for multi-phase flow including porous flow, a volume of fluid module (VFM) for airwater interface tracking, and a sediment transport module (STM) for fluid-sediment interface tracking. In the FSM, a two-way coupling procedure is implemented at each time step to connect the GNS with the VFM and the STM. The predictive capability of the FSM is demonstrated through comparison between numerical results and experimental data in terms of water surface elevation, inundation depth, and tsunami force. The process of tsunami run-up in the presence of square structures is investigated in terms of vortex structures. The result shows that the FSM is a useful tool providing detailed information in discussing run-up tsunami deformation and tsunami force.


INTRODUCTION
Tsunami-resistant structures have been constructed to serve as emergency evacuation sites in case of tsunami attack.Such structures are also expected to reduce inundation depth in coastal areas and tsunami force acting on inland buildings.Tsunami force has been investigated in many experimental and numerical studies (e.g., Asakura et al. 2000, Ikeno and Tanaka 2003, Ikeya et al. 2005, Nakamura et al. 2007).However, most of the studies have been limited to investigations into a single structure.Simamora et al. (2007) carried out a series of hydraulic experiments on tsunami force in the presence of multiple cubic structures, and demonstrated the efficiency of tsunami-resistant structures to reduce tsunami force on inland structures.Although the tsunami force was analyzed in terms of inundation depth around the structures, the mechanism of reducing the tsunami force was inadequately discussed because of the lack of enough experimental data.Furthermore, little research has been undertaken to predict and investigate tsunami force in the presence of multiple structures with a three-dimensional numerical model.
In this study, run-up tsunami deformation and tsunami force in the presence of multiple structures (Simamora et al. 2007) are investigated with a three-dimensional two-way coupled fluid-sediment interaction model (hereinafter referred to as FSM; Nakamura and Yim 2010).In the following section, brief explanation of the FSM is presented for completeness.Next, the predictive capability of the FSM is verified through comparison between numerical results and experimental data.Finally, run-up tsunami deformation in the presence of structures is discussed in terms of vortex structures to demonstrate the usefulness of the FSM in discussing run-up tsunami deformation and tsunami force.GNS and the VFM are briefly explained in the following sections.Note that the STM is not explained below because sediment transport and resulting seabed change are not dealt with in this paper.(Brackbill et al. 1992); i R  laminar and turbulent drag force vector due to porous media (Mizutani et al. 1996

Assuming
strain rate tensor; * q  intensity of wave generation source/sink per unit time (Kawasaki 1999); i Q  wave generation source/sink vector; and ij   artificial damping factor matrix.In the FSM, the x axis is the cross-shore direction, the y axis is the long-shore direction, and the z axis is the upward direction.In deriving Eq. ( 2), the spatial variation in the porosity is taken into account because of possible sharp change in the porosity around the surface of porous media.However, based on the formulation of CADMAS-SURF (CDIT 2001), the spatial variation in the porosity is assumed to be negligible only in deriving the pressure gradient term of Eq. (2) (the first term of the right hand side of Eq. ( 2)).This is to ensure the equilibrium between the pressure gradient term and the gravitational acceleration term (the second term of the right hand side of Eq. ( 2)) in still water regardless of the spatial change in the porosity.In Eq. ( 2), the formulation of s i f , i R , i Q , and ij  are respectively expressed as ( , , ) in artificial damping zone 0 o t h e r w i s e x y z   artificial damping factor.As expressed in Eq. ( 7), the artificial damping factor matrix ij  is introduced to dissipate outgoing waves in an artificial damping zone.Based on the formulation of ( , , ) x y z  used in Hinatsu (1992) and Cruz et al. (1993), ( , , ) x y z  in Eq. ( 7) is given as ( , , ) sin min ,1 2 in which C   artificial damping factor coefficient; h  still water depth;   length of the damping zone; x   horizontal distance from the inside boundary of the damping zone; z   vertical distance from the bottom of the computational domain; and   water surface elevation at x  (see Fig. 2).
Based on the dynamic two-parameter mixed model (DTM; Salvetti and Banerjee 1995), the anisotropic part of the turbulent stress tensor a ij  in Eq. ( 2) is given as (Morinishi and Vasilyev 2001) ; and   filter size).The superscripts  and ~ represent the values in the grid scale and the test scale, respectively.In the FSM, a box filter in physical space is adopted as the grid-scale filter and the testscale filter.The grid-scale filter size in each direction  is set at the size of numerical cells in the corresponding direction.Once the value of  is determined, the value of a ij  is computed from Eqs.
(9) to (11) with only i v .In this study, the value of  is set at 2.0 (Germano et al. 1991).
The values of all parameters used in the FSM are listed in Table 1.As explained later, neither porous media nor artificial damping zones are set in a computational domain in this study.

Accordingly, the values of
C , and C  are not required to be set.However, 0.04  (1996), and the value of C  is determined at 3.0 through trial and error in preliminary runs.

Computational Schemes of GNS and VFM
In the FSM, the simplified marker and cell (MAC) method (Amsden and Harlow 1970) is used to discretize the time derivative of the continuity equation (Eq.( 1)) and the generalized NS equation (Eq.( 2)).Based on Kajishima (1999), the pressure gradient term and the gravitational acceleration term (the first and second terms of the right hand side of Eq. ( 2)) are discretized with the first-order Euler forward-difference scheme.To ensure the stable calculation, the laminar drag force term (the first term of the right hand side of Eq. ( 5)) is discretized with the second-order Crank-Nicolson scheme.The other terms are discretized with the third-order Adams-Bashforth scheme for the accurate calculation.As for the space derivative, the central-difference scheme is used for Eqs. ( 1) and ( 2) except that the convective term of Eq. ( 2) (the last term of the left hand side of Eq. ( 2)) is discretized with the thirdorder total variation diminishing (TVD) scheme (Osher and Chakravarthy 1984).Accordingly, the difference equations at the prediction step and the correction step are respectively expressed as in which p i v  predicted flow/seepage velocity vector; the superscript n  time step number; 1/ 2 n t    time increment between the n th time step and the   from the following Poisson equation, which is derived by differentiating Eq. ( 13) with respect to i x .
In Eq. ( 12) to ( 14), the values of 0 A , and n B are respectively given as To track air-water interface motion, the advection equation of the VOF function (Eq.( 3)) is solved with the multi-interface advection and reconstruction solver (MARS; Kunugi 2000).The MARS is one of the piecewise linear interface calculations (PLICs), in which the air-water interface in each cell is expressed as an inclined plane.Detailed explanation of the MARS can be found in Kunugi (2000).

NUMERICAL CONDITIONS
Run-up tsunami deformation and tsunami force in the presence of multiple cubic structures (Simamora et al. 2007) are analyzed with the FSM. Figure 3 shows a schematic of a computational domain.The domain has the same dimensions as the experimental setup except that only a half side of the domain ( 0 y  ) is computed to reduce the computational effort because the experimental setup is practically symmetric with respect to 0 y  .As shown in Fig. 3, a 3700 mm long, 1925 mm wide, and 620 mm high impermeable horizontal land is set at the landward end of the domain, and an 1810 mm long, 2925 mm wide, and 570 mm high impermeable horizontal bed and an inclined impermeable flat bed with a slope of approximately 1/3 are set in front of the land.
The size of numerical cells is determined to ensure an appropriate balance between the predictive accuracy and the computational effort.Figure 4 shows the numerical cells adopted in this study.The domain of 0 700mm x   , 0 150mm y   , and 50 100 mm z    is discretized with uniform cells of 20 12.5 10 mm   .The remainder of the entire domain is discretized with non-uniform cells whose width increases by approximately 5 to 20% with the distance from the domain discretized with the uniform cells.The non-slip condition is used for the impermeable beds and the side walls.The gradient-free condition is used for the landward boundary and 0 y  .To generate a tsunami, the flow velocity at the seaward boundary is specified to be the speed of a wave paddle recorded in the hydraulic experiments.As for the VOF function, all boundaries are exposed to the gradient-free condition.
A total of seven numerical runs are carried out.In Case 1, no structure is set on the land.In Cases 2 to 7, one or two cubic structures with the dimensions of 100 100 100 mm   are set on the land.The arrangement of the structures is shown in Fig. 5.Among Cases 2 to 7, the landward structure is shielded from a run-up tsunami by the seaward structure in Cases 4 to 7, while no shielding structure exists in Cases 2 and 3.In all cases, the still water depth is set at the constant value of 600 mm.
To verify the validity of the FSM, numerical results are analyzed and compared with experimental data measured in Simamora et al. (2007).In the hydraulic experiments, water surface elevation ( x  -4300, -3800, -3300, -2800, -2300, -1750, -1250, -750, and -250 mm on  0 mm) and inundation depth ( x  200, 300, 400, 500, 600, 700, and 800 mm on y  0 mm) were measured in Case 1.The measurement of inundation depth in front of the landward structure ( x  200 mm, y  0 mm in Cases 2 and 4, and x  600 mm, y  0 mm in Cases 3, 5, 6, and 7) and tsunami force acting on it were carried out in the other cases.The comparison of the water surface elevation, the inundation depth, and the tsunami force is discussed in the following sections.

Predictive Capability of FSM
To demonstrate the predictive capability of the FSM, numerical results computed with FSM are compared with experimental data measured in the hydraulic experiments (Simamora et al. 2007).Figures 6 and 7 show the comparison of the water surface elevation and the inundation depth in the absence of structures (Case 1).In the figures, the solid lines represent numerical results, and the circles represent experimental data (Simamora et al. 2007).As shown in Fig. 6, the propagation of an incident tsunami and subsequent reflected waves is predicted well.Note that the FSM overestimates the backwash on the uniform shallow water area ( x  -1750 to -250 mm).The experimental data in Fig. 6 show that the backwash and its landward propagation are observed in the offshore side of x  -2300 mm, while the backwash less than -4 mm is not observed in the onshore side of x  -1750 mm.This is probably because a decrease in the water level was not measured in the hydraulic experiments because of the very shallow water depth of 30 mm and the limitation of measurement devices.As for the inundation depth, Fig. 7 shows that the numerical results also agree with the experimental data except that the maximum inundation depth in the offshore side of x  400 mm is slightly underestimated.As explained later, no vortex is formed on the onshore side of the land, while strong vortices with air bubble entrainment are formed around the seaward edge of the land.Because of the disturbed flow, the increase in the inundation depth is not tracked near the sea.
Figures 8 to 12 show the comparison of the inundation depth and the tsunami force in the presence of structures (Cases 2 to 7).In Cases 2 and 3 for the absence of shielding structures, Figs. 8 and 9 show that the numerical results reasonably agree with the experimental data.Specifically, the rise time of the inundation depth and the tsunami force is predicted well, while the maximum tsunami force is  slightly overestimated in both cases.As indicated in Figs. 10 to 13, both the inundation depth and the tsunami force are in agreement between the numerical results and the experimental data in Cases 4 to 7 for the presence of the shielding structures.Figures 10 and 11 also show that the tsunami force is slightly underestimated for the presence of the shielding structure located 100 mm offshore from the onshore structure.In addition, there is a slight difference in the rise time of the inundation depth in Fig. 10(a).However, Fig. 13 shows that the FMS gives the good predicted values of both the inundation depth and the tsunami force in spite of the complicated arrangement of the structures.
As a result, the predictive capability of the FSM is verified in terms of the water surface elevation, the inundation depth, and the tsunami force regardless of the presence of the structures.

Run-Up Tsunami Deformation and Vortex Structure
Figure 14 shows the process of tsunami run-up in the presence of the structure located 600 mm onshore from the sea (Cases 3, 5, and 7).Among three cases, no shielding structure exists in Case 3, while the shielding structures are set in front of the landward structure in Cases 5 and 7.In Fig. 14, vortex structures are visualized with the 2  definition (Jeong and Hussain 1995), in which vortices are defined as 2   0. As pointed out in Miura and Kida (1998), it is difficult to identify vortex structures with 2   0 because a wide area of the computational domain is possibly covered by isosurfaces of 2   0. To identify vortex structures more clearly, isosurfaces of 2  -30 are presented in Fig. 14.
Note that the tsunami begins to be generated at t  0.0 s.As shown in Fig. 14(a), the bore-like tsunami propagates close to the land at t  8.90 s.At the same time, vortex cores are formed under the crown of the tsunami.Note that the wave field including the vortex cores is not uniform in the cross-shore direction because the water area slightly narrows near the toe of the slope (see Figs. 3 and 4).In addition, it is also observed from Fig. 14(a) that there is a slight difference in the wave field among three cases because the air flow over the wave is affected by the arrangement of the structures.Figures 14(b) and (c) show that large vortex cores are formed around the seaward edge of the land early in the tsunami run-up.In Case 7, the tip of the run-up tsunami reaches the seaward structure at t  9.35 s.After that, the run-up tsunami subsequently propagates in between and around the structures in Case 7 (Fig. 14(d)).At this time, vortex cores appear around the seaward edges of the structures, probably resulting in the dissipation of the tsunami energy.As shown in Fig. 14(e), the propagating tsunami also reaches the landward structure at t  9.75 s.Note that the tsunami propagating between the seaward structures hits the seaward surface of the landward structure in Case 7, while the tsunami propagating around the seaward structure hits the seaward edges of the landward structure in Case 5. Consequently, the maximum tsunami force acting on the landward structure in Case 5 is smaller than that in Case 7 (see Figs. 11(b) and 13(b)).After that, the tsunami propagates farther around the structure (Figs.14(f) and (g)).At this time, the inundation depth in front of the landward structure remains high in Case 3, while the inundation depth gradually decreases in Cases 5 and 7 because the shielding structures prevent the subsequent run-up tsunami from hitting the landward structure directly.As shown in Figs. 9, 11, and 13, the inundation depth and the tsunami force are therefore reduced in Cases 5 and 7 because of the presence of the shielding structures.Finally, the run-up tsunami begins to be subsided, as shown in Fig. 14(h).
From the results explained above, it is found that the FSM is a useful tool providing detailed information including vortex structures in discussing run-up tsunami deformation and tsunami force.

CONCLUDING REMARKS
A three-dimensional two-way coupled fluid-sediment interaction model (FSM; Nakamura and Yim 2010) is applied to predict hydraulic experiments on run-up tsunami deformation and tsunami force in the presence of multiple cubic structures on the land (Simamora et al. 2007).The comparison between numerical results computed with the FSM and experimental data measured in the hydraulic experiments demonstrates the predictive capability of the FSM in terms of water surface elevation, inundation depth, and tsunami force.The process of tsunami run-up in the presence of the structures is investigated in terms of vortex structures visualized with the 2  definition, and it is found that the FSM is a useful tool in discussing run-up tsunami deformation and tsunami force.However, there is a slight difference in inundation depth and tsunami force depending on the arrangement of structures.To address the issues, further studies are required to be conducted to improve the predictive accuracy of the FSM.

Figure 1 .
Figure 1.Typical computational domain of the FSM.
sufficiently small temporal variation in the porosity m representing the volume fraction of void space in each cell ( pure fluids), Nakamura and Yim (2010) derived a continuity equation, a generalized Navier-Stokes (NS) equation, and an advection equation of the volume of fluid (VOF) function as follows: 3) in which i v  fluid/seepage flow velocity vector; p  pressure; F  VOF function representing the volume fraction of water in each cell ( water surface; and 1 F  for water); i x  position vector; t  time; viscosity of fluid ( w  and a   those of water and air, respectively); A C  added mass coefficient; s i f  surface tension force vector based on the continuum surface force (CSF) model

Figure 2 .
Figure 2. Artificial damping zone in a typical computational domain of the FSM.Table 1. Parameters in the FSM.Gravitational acceleration g

Figure
Figure 5. Structure arrangement in Cases 2 to 7.

Figure 3 .
Figure 3. Schematic of the computational domain.

Figure 4 .
Figure 4. Numerical cells in the computational domain.

Figure 6 .
Figure 6.Comparison of the water surface elevation in the water area (Case 1; no structure).

Figure 7 .
Figure 7.Comparison of the inundation depth on the land (Case 1; no structure).
of run-up tsunami deformation and vortex cores.