Influence of Tsunami-Driven Shipping Containers’ Layout on Their Motion

Abstract

This study investigated the interaction between containers under extreme hydrodynamic conditions modeled on tsunamis to assess whether the number and layout of containers affect their motion and to guide future studies on numerical simulations describing tsunami-debris motion. The three-dimensional Fluid–Structure–Sediment–Seabed Interaction Model (FS3M), which can compute tsunami–container interaction, was used as a numerical test model, and numerical results for the specific target of the simulation were compared with experimental data to check the validity and computational accuracy of the FS3M. The study showed that the number of rows (N_x), columns (N_y), and stacks (N_z) in the initial arrangement of containers constitute the main factor affecting the area where the containers spread and their drift motion velocity. An increase in N_x and N_z can effectively reduce the container drift velocity. Conversely, as N_y increases, the drift motion velocity of the center of gravity of the entire group increases. The results of this study can facilitate the development of more realistic building structure scenarios in future research that consider the proposed characteristic damage estimation and comprehensive assessment methods laid out herein.

1. Introduction

Tsunamis have caused significant damage to coastal communities in recent years and are generally followed by complicated recovery processes. More recent tsunami disasters have provided considerable information on the extent of damage caused by debris impact on most existing structures, especially critical infrastructure. As shipping containers are located at all ports globally, floating shipping containers constitute a significant component that can cause structural failure within their flow path in many cases. For example, a field survey was carried out on the structural and infrastructural effects of the earthquake that occurred off the coast of central Chile on 27 February 2010 and the subsequent tsunami that struck Chile, Peru, and other countries along the Pacific Ocean [1]. It was discovered that close to 700 shipping containers from the Talcahuano port had drifted shoreward and piled up in front of and inside a building structure. Moreover, drifting shipping containers entrained by the tsunami may lead to secondary destruction of the coastal structure through the impact of tsunami-debris. Similar conclusions were reached after the Samoan tsunami in 2009 [2]. Therefore, destruction of this scale emphasizes the importance of considering tsunami-driven container transport when regulating design guidelines for buildings in coastal communities in order to maximize the effectiveness of coastal disaster mitigation strategies and improve the resilience of coastal communities to future extreme events.

Scientific assessment of disaster debris management is key to improving the effectiveness of disaster mitigation, prevention, and relief strategies and ensuring the sustainable development of society. It is also an essential scientific basis for tsunami relief, recovery, and reconstruction decision-making. Stolle et al. [3] used a probabilistic design approach to analyze shipping container motion in extreme tsunami events and presented a characterization of the stochastic properties of debris transportation. This research aimed to exploit a stochastic framework for debris hazard assessment and to apply the framework to a tsunami-like scenario. The results showed that hydrodynamic conditions, debris geometry, and the initial arrangement configuration of debris were identified as the most critical factors affecting the lateral distribution of debris in tsunami-entrained flows. This study demonstrates a method of expressing this stochasticity in numerical models of a sufficiently large scale. More recently, Park et al. [4] conducted a hydraulic experimental study of multi-debris motion with two materials of different densities, considering the debris components for the different densities. The purpose was to better understand the comprehensive process of multi-debris transport and its distinguishing features and to evaluate how to reduce the damage to the structured environment during debris transport. They found that the density of each debris element is the main factor that determines the probability of collision. However, there is a need to better understand the potential impact of different densities of debris on structures in realistically reflecting structures along urban shorelines and how to provide an initial database for numerical studies of tsunami-induced debris motion.

Over the years, many researchers have focused on debris management assessment of disasters and field reconnaissance of the disaster debris impact to help comprehend tsunami hazards. However, only a few studies have focused on the accumulation of debris and the transport phase in hydraulic experiments and numerical simulation. Several researchers have conducted hydraulic experiments to investigate detailed measurements of the flow field around a single piece of shipping container and experimentally determined how debris impact is affected by the surrounding fluid [5,6]. Yeom et al. [7,8,9] conducted hydraulic experiments and numerical simulations to investigate the drift behavior of a single piece of debris and proposed a calculation scheme capable of investigating the collision forces and collision times resulting from shipping container drift caused by tsunamis. Collision deformation prediction and collision force estimation were also performed to confirm the applicability of collision analysis to reduce such damage. Several researchers have provided insight into numerical simulation studies related to individual debris–fluid-bridge structure interactions and impact loading [10,11]. For practical purposes, the impact of single debris on structures and the development of analytical solutions are valuable. However, there is ample video evidence from recent tsunamis showing that debris tends to propagate in aggregate form. This phenomenon requires more in-depth research and investigation with theoretical support, which may affect current design guidelines. To this end, prospective studies need to consider methods for developing disaster prevention measures through numerical models.

The ASCE/SEI 7-22 standard provides an arrangement procedure for shipping containers before a significant tsunami. If a port is located in the surrounding region of a shipping container yard, this standard recommends determining whether the impact of shipping containers should be considered in the design of the structure and evaluation. Nakamura et al. [12,13,14] considered shipping containers’ drift behavior onshore caused by tsunami rise and fall through numerical simulation models. They consequently determined how to effectively prevent shipping containers located on the apron from falling into the sea along the cross-shore direction. However, the mechanism to determine whether the containers fell into the sea was unsatisfactory owing to a lack of detailed data. Haehnel et al. [15] used a hydraulic experiment to collide wood debris with simplified vertical structures to investigate the maximum impact force of the floating debris on the structures. In addition, a reliable method to estimate the maximum impact force was proposed. Naito et al. [16] proposed a reliable procedure for quantifying the impact potentiality of debris in tsunami-prone sites. They considered the forces and potential damage caused by debris impacts more intuitively in the structural design guidelines. Further, a method for calculating debris transport potential was developed, considering the initial angle and the various drift patterns of the inland wall affecting the container. However, to understand the reasons for the supposition, the fundamental source of the provisions, and the standard method for obtaining background data, it is necessary to improve the effectiveness of debris management assessment of disasters.

The physical phenomena concerning debris transportation in extreme hydrodynamic conditions have not been fully explored, and even less attention has been given to the interaction of multi-debris transport. Although most researchers acknowledge that the interaction of debris transport is crucial, little research has been performed on a probabilistic basis to scrutinize how debris affect other debris as they spreads away from their initial location. Tomita et al. [17] performed numerical simulations and hydraulic model experiments using a tsunami simulator called T-STOC (“Storm Surge and Tsunami simulator in Oceans and Coastal areas”) to analyze the variation in debris (large and small rectangular objects) positions. They found that the center of gravity of the debris flows in approximately the same manner along the inundation flow. However, to increase the accuracy of the numerical computation, it is indispensable to fully contemplate the flow and interaction between the drifting objects when the inundation flow hits the drifting object population. Tomita et al. [18] conducted hydraulic model experiments to investigate the influence of different shipping container arrangement schemes and reported that the debris interaction affected the drift position. Thus, there is a need to consider numerical calculations of tsunami drifters with regard to the interaction between multiple containers and the variation in position under the same conditions.

The physical experimental model poses significant limits on studying new configurations in terms of time and expense. To overcome these limits, Alessandro et al. [19] proposed a new numerical simulation method based on the coupling between the Discrete Element Method (DEM) and the Lattice-Boltzmann Method (LBM). The model can simulate debris flow as a mixed fluid-particle medium with exceptional tracking of the largest grains. In addition, the debris flow and its interaction with the coastal structures were successfully simulated. Istrati et al. [20] used the Finite Element Method (FEM) to numerically investigate debris loads on containers trapped below the decks of a coastal bridge structure. Hasanpour et al. [21] developed a coupled SPH-FEM simulation model. The method can accurately capture the debris–fluid interaction and the effect on the coastal structures. This study aimed to assess the interaction between containers within extreme hydrodynamic conditions using a numerical simulation model to help comprehend how multiple objects (debris) affect their drifting motion through different initial container layouts, whether the initial number and configuration of containers affect their motion, and how and under what conditions the different drifting motion results appear. This is with the objective of improving the motion behavior assessment of tsunami drifters and the initial arrangement of containers in a tsunami-prone site. The numerical test model used was the three-dimensional coupled Fluid–Structure–Sediment–Seabed Interaction Model (FS3M) developed by Nakamura et al. [22], which can compute fluid–structure interaction. The numerical simulation results were compared with hydraulic experimental data [18] from the specific target of this study to check the validity and computational accuracy of the FS3M. Further, we developed a framework to help calibrate and validate the FS3M to guide future numerical simulation experiments describing tsunami-debris motion.

The remainder of this paper is organized as follows. Section 2 briefly describes the hydraulic experiment of Tomita et al. [18]. Section 3 describes the numerical experimental procedure. Section 4 examines and discusses containers’ center-of-gravity drift behavior, container velocity, and spatial container distribution. Section 5 summarizes the general findings of this study and suggests possible future works.

2. Outline of Hydraulic Experiments

As shown in Figure 1, the specific target was the simulation experiments conducted by Tomita et al. [18] in the Department of Civil and Environmental Engineering of Nagoya University, Japan. The concrete wave flume was 25.0 m in length (x-direction), 2.2 m in width (y-direction), and 1.0 m in depth (z-direction). A piston-type wavemaker with a ceiling stroke of 1.50 m was used that can generate 1/4 period sine-wave-type waves in the form of extruded water. The profile of the wave flume consisted of a 4.45 m long horizontal flat section in front of the wavemaker, approximately 1:6 impermeable slope expanding over 4.45–6.44 m from the wavemaker, and a 4.02 m long and 0.32 m high horizontal flat impermeable bed at 6.44–10.46 m from the wavemaker. This experimental setup was selected to assume tsunami inundation with a relatively gentle slope and flat geography to model typical port container facilities of low-lying coastal areas. The impermeable slabs were made of polyvinyl chloride (PVC) boards to create a reduced friction coefficient rather than the concrete wave flume, and this effectively prevented the drift containers from being blocked in the flow motion. Silicone was applied to the joints of the PVC slabs to prevent water intrusion into the joints. A vertical slab waveguide was installed in the part of the flume without PVC slabs in the cross-shore direction, and permeable meshes with a synthetic resin material called Hechimaron (Shinko-Nylon Co., Osaka, Japan) were applied in the gap between the waveguide and the flume. It was confirmed in preliminary experimental runs that the slab waveguide and permeable meshes had almost no effect on inundation waves on the wave flume.

Figure 2 shows that 12.19 m (40 ft) shipping containers were accurately reduced on a length scale of 1:74 according to the Froude similarity law. A standard intermodal unladen container is 12.19 m (40 ft) in length, 2.44 m (8 ft) in width, and 2.93 m (9 ft 6 in.) in height, with a mass of 3800 kg, stiffness was determined to be 29,800 kN/m. The physical experimental model of the shipping containers was 163 mm long, 33 mm wide, and 35 mm tall, with a mass of 88.7 g, and was hollow. The debris was made from acrylic with an approximate draft before the physical experiments. However, a 2 mm foot was attached to the bottom of the shipping container model to reduce the contact area between the drifter container model and the wave flume. A 1/4 period sine-wave-type wave with a period of 20 s was generated using the wavemaker. (The motion of the shipping container model caused by the incident first wave generated by the wavemaker at approximately 5.0 s of the splitting wave of each hydraulic experimental test is explained in detail in Section 3.3). At this moment, the wave caused a slight breaking wave; however, the wave was not affected by the breaking wave but flowed along the horizontal direction of the inundation flow. Therefore, the shipping container model was considered to be in a floating state due to the incident waves.

Two overhead high-speed integrated cameras (HAS-L1, Digital Image Technology, Tokyo, Japan) were mounted above the flume of the physical experiment section. They were used to track drifting container flow behavior and velocity. These high-speed integrated cameras were simultaneously sampled at 300 fps. The cameras were synchronized with the data acquisition system, and tracking was performed via an image analysis system (Dipp-Motion V, Digital Image Technology, Tokyo, Japan) every 1/100 s.

The Tomita et al. [18] simulation experiments had 13 arrangement modes in total. The center-of-gravity position (X_G, Y_G) of the debris group as a population was calculated from each container model’s average center-of-gravity position for each experimental result. The flow behavior of the drifting containers deviated from the results even for experiments conducted under the same conditions; thus, 10 iterations of the experiment were conducted under the same conditions. In addition, the ensemble average (X_GA, Y_GA) of all 10 iterations of experiments was derived. Furthermore, the deviation values of the results of each experiment were evaluated by the standard deviation criteria (X_SA, Y_SA) about the center-of-gravity position of the population. In addition, the deviation of each model in the population was evaluated by the standard deviation (x_S, y_S) of the center-of-gravity position of each container model relative to (X_G, Y_G); subsequently, the ensemble average population value (x_SA, y_SA) of the results of all 10 iterations was derived. The degree of interaction between containers in the first and second rows and containers in the second and third rows was also examined.

3. Numerical Model and Conditions

3.1. Outline of Numerical Model and Its Improvement

The FS3M [22] consists of a main solver and four modules. The main solver is a large-eddy simulation (LES) model based on the continuity equation and Navier–Stokes equation for computing an incompressible viscous air–water two-phase flow considering the motion of movable objects. The volume-of-fluid (VOF) module is based on the multi-interface advection and reconstruction solver (MARS) [23] to track the air–water interface motion. The immersed-boundary (IB) module is based on the body–force type IB method [24] for fluid–structure interaction (FSI) analysis of movable objects. The other two modules include a sediment transport module for bedload and suspended sediment transport, and a finite-element-model (FEM) module for coupled soil–water analysis. In this study, the main solver and the VOF and IB modules in the FS3M were employed to compute the tsunami–container interaction. Full details of the governing equations and computational scheme of the FS3M can be found in Ref. [22]. In addition, open dynamics engine (ODE) v0.16.2 [25], a physics engine, was incorporated into the IB module to deal with the contact between containers.

The FS3M–ODE combination was first applied to experiments on toppling dominos [26] to verify that ODE was properly incorporated into the FS3M. Following the experimental condition, 30 dominos with an x-directional thickness T_d of 7.70 mm, y-directional width of 21.90 mm, z-directional height of 43.20 mm, and a density of 631 kg/m³ were arranged with an x-directional spacing $\lambda$_d of 13.86 mm $\left({\lambda }_d/T_d=1.8\right)$ on the horizontal bottom in a sufficiently large computational domain filled with air. The domain was discretized with uniform rectangular meshes with an x-directional size of 1.54 mm, y-directional size of 3.65 mm, and z-directional size of 4.32 mm. The domino wave was initiated by applying an x-directional horizontal force of 0.05 N to the first domino for the first 0.01 s. The friction pyramid approximation was used to evaluate the Coulomb friction in ODE. The friction coefficient between the dominos, the friction coefficient between the dominos and the bottom, and the restitution coefficient between the dominos were set to the same values as in the experiments, that is, 0.11, 0.90, and 0.54. The restitution coefficient between the dominos and the floor was set to 0.00 owing to the lack of experimental data following Li et al. [27]. The value of the constraint force mixing (CFM) was determined to be $1.0\times {10}^{-5}$ for stable computation. A default value of 0.2 was used as the error reduction parameter (ERP).

Figure 3 compares the toppling dominos. In Figure 3a, which shows the experimental result [26], t is the time from the moment when the eleventh domino from the left contacts the twelfth domino. In Figure 3b, which shows the numerical result, t is the time from the simulation start and the vectors represent the flow velocity on the central cross-section of the domino array. From Figure 3, it is clear that the domino wave is almost the same between the experiment and the simulation. Furthermore, Figure 3b shows that the flow of air induced by the leaning dominos can be computed using the FS3M–ODE combination.

Figure 4 compares the domino wave speed, which is called the intrinsic collision speed in Ref. [26]. As shown in Figure 4, the computed speed of the domino wave before the tenth domino is slightly larger than that of the experiment. This is probably because the speed of the first several dominos is affected by the method of moving the first domino, and it differs between the experiment and the simulation. However, Figure 4 shows that there is good agreement in the almost constant speed of the domino wave after the tenth domino, suggesting that the contact between the dominos can be computed properly. The computational capability of the FS3M–ODE combination was thus verified and was used for the subsequent experiments.

3.2. Computational Domain and Container Specimen Configurations

Figure 5 shows the wave flume reproduced in the FS3M computational domain. The profile of the wave flume contained a 4.45 m flat section in front of the wave generation boundary, with an approximately 1:6 slope protracted from x = −2.49 m to x = −0.50 m. Container configurations started from x = 0.0 m, which was 0.5 m from the end of the slope, and extended a 2.40 m computational domain flat section, with an impermeable bed elevated 0.32 m. To minimize the scope of the computational domain and thereby the computational time, the length of the computational domain L_y was changed depending on the layout quantity of the container in the y-direction. For this numerical simulation study, the container layout in the y-direction was divided into four types with N_y values of 1, 2, 3, and 4. Each layout corresponds to L_y widths of 0.2640 m, 0.3630 m, 0.4950 m, and 0.6270 m to ensure at least a gap of one container between the container and the side boundaries. Because the starting point of the numerical computation domain is positioned centrally in the domain, for the final domain width in the y-direction, the value of L_y needs to be multiplied by two. These details are illustrated in Figure 5a.

The mesh size in the x-direction was applied to 0.00725 m around the container placement frame landward. However, this value was extended to 0.02 m seaward towards the wave generation boundary. The y-direction was kept at 0.00725 m. The z-direction was applied to 0.0025 m from 0.05 m below the initial water level to the highest value of the container stack at 0.148 m, but this value was extended to 0.02 m at another cell size. To capture the free surface motion near the initial water level and container placement frame, a satisfactory horizontal resolution with sufficient mesh size is essential. A non-uniform mesh was used at another cell size, which maintained accuracy without excessive computational cost. The mesh size was specified to maintain computational accuracy while obtaining acceptable computational cost, depending on the proper balance between the two. The mesh size was determined to be of sufficient resolution for preliminary simulation runs by matching the incident waveform verification and comparison of experimental data with numerical results, as discussed in the following sections.

Table 1. Numerical simulation protocol.

Experimental Category	Number of Containers (N)	Container Orientation (θ)	Number of Rows (Nx)	Number of Columns (Ny)	Number of Stacks (Nz)	Interval of Container W (mm)
11190	1	90	1	1	1	33
21190	2	90	2	1	1	33
31190	3	90	3	1	1	33
41190	4	90	4	1	1	33
11290	2	90	1	1	2	33
21290	4	90	2	1	2	33
31290	6	90	3	1	2	33
41290	8	90	4	1	2	33
12190	2	90	1	2	1	33
22190	4	90	2	2	1	33
32190	6	90	3	2	1	33
42190	8	90	4	2	1	33
12290	4	90	1	2	2	33
22290	8	90	2	2	2	33
32290	12	90	3	2	2	33
42290	16	90	4	2	2	33
13190	3	90	1	3	1	33
23190	6	90	2	3	1	33
33190	9	90	3	3	1	33
43190	12	90	4	3	1	33
13290	6	90	1	3	2	33
23290	12	90	2	3	2	33
33290	18	90	3	3	2	33
43290	24	90	4	3	2	33
14190	4	90	1	4	1	33
24190	8	90	2	4	1	33
34190	12	90	3	4	1	33
44190	16	90	4	4	1	33
14290	8	90	1	4	2	33
24290	16	90	2	4	2	33
34290	24	90	3	4	2	33
44290	32	90	4	4	2	33

summarizes all the details for each layout and shows the effects of the number of debris (N), number of drifter containers in rows (N_x), number of drifter containers in columns (N_y), and stacking height (N_z). The entire numerical simulation container orientation (θ) was set at 90 degrees and the container interval W (mm) was set at 33 mm. The five numbers in the numerical simulation name represent the number of containers in the x-direction, y-direction, z-direction, and container orientation (θ). Diverse types of 32 shipping container arrangement groups were chosen and tested in this numerical simulation experiment, highlighting transmutations in shipping container orientation and quantity. These are the guiding parameters most relevant to the process being investigated. Existing standardized shipping container storage operations are capable of stacking heights of a minimum of five layers [28].

Figure 6a,b show an example of the initial configuration container for the experimental category of Configuration 42290 (see Table 1). Nevertheless, the reality is that numerous other shipping containers may be arranged in considerable numbers and on more capacious surfaces in an actual coastal harbor. However, the number of chosen layers for shipping container models was constrained because of experimental shipping constraints. Thus, a maximum of 32 containers could be selected. It can be observed in Figure 7 that N_z = 2 corresponds to shipping containers arranged on top of the same N_z arrangement with no additional connection. As mentioned above, the containers were arranged in an identical position with an equal horizontal interval (0.033 m) for several numerical simulation experimental runs. The containers were arranged in diverse initial locations to analyze the impact of the initial arrangement on container transportation. The initial arrangements were established based on the number of containers (N), number of rows (N_x), number of columns (N_y), and number of stacks (N_z).

3.3. Validation of Incident Waveform without Containers

Figure 8 explains the experimental data incident waveform η of the water surface measured at the position of the flume wave gauge in Figure 1. Before executing the wave-driven container numerical simulations, preliminary research was carried out to identify and calibrate the incident wave condition to match the incident wave, yielding the best-matched results under several trials without containers. This incident wave contains the first wave generated by the wavemaker in the form of extruded water at approximately 5.0 s of the splitting wave. The splitting wave of this incident wave, affected by soliton fission, causes some wave breaking but does not reach whole wave fragmentation. The basic data for the incident waveform from Tomita et al. [18] was derived from images taken from the side of the wave flume with the maximum inundation depth produced by this wave over the model flume.

The incident wave profile is essential because this is the driving force moving shipping containers. For the best matching attributes, different numerical simulation results were derived by changing the wave height and wave period of the numerical simulation model FS3M input. The numerical computation results were compared with the physical experimental data to determine the optimal value of the wave height and wave period. Only in this calibration, the FS3M computational domain was changed to a vertical two-dimensional space to reduce computational time.

Figure 8 compares the numerical simulation and experimental data results for various wave height and wave period values. The peak of each numerical simulation and experimental data occurred at approximately 7.0 s, and it is generally considered that the incident waveform can be distinguished into a first half and second half at the wave peak. It was observed through several numerical simulations that a change in the set value of the wave height affects the change in the first crest (first half), while a change in the value of the wave period affects the change in the second wave crest (the second half). By comparing (b,c) in Figure 8, it can be concluded that when the wave periods are the same and the wave heights are different, the set value of the wave height affects the change in the first crest (first half), as described above. In Figure 8, the wave heights of (a–c) are 0.083 m, 0.082 m, and 0.081 m. It can be observed that the wave height of 0.081 m is the incident wave value that best fits the numerical simulation test. If the wave period is selected in the interval of approximately 13.0–14.5 s, the second crest (second half) is more consistent. However, if the wave period is in the range of approximately 17.0–22.0 s, the first crest (first half) matches better. Nevertheless, the debris began to move due to the initial first crest (first half) in the experiments. The second crest (second half) can be affected by dissipating waves at the onshore boundary between experiments and simulations. As a result, the 17.0–22.0 s wave period range was not selected. Upon comparing Figure 8c–f, it was determined that a 14.4 s wave period best fits the numerical simulation test incident waveform.

4. Results and Discussion

4.1. Comparison between Numerical Results and Experimental Data

As described in Section 3, the FS3M–ODE combination can compute the interactions between containers (contact and friction) and the interactions between containers and the floor (friction) to solve the motion of the drifting containers. However, it is necessary to calibrate the value of the coefficient of friction (FC) between the ground (chloroethylene slabs in the experiment) and the container (acrylic in the experiment). Using the original data of the hydraulic experiment for Configuration 11900 (hydraulic experiment) obtained by Tomita et al. [18], the density of the container was set to 467.9 kg/m³.

Figure 9 shows the same Configuration 11900 (hydraulic experiment) with different FC values utilizing the FS3M to calculate the x-directional center of gravity displacement position. The x-directional center of gravity’s original displacement position of the hydraulic experimental [18] was compared with the FS3M numerical simulation results. If the FS3M results are lesser or greater than the hydraulic experimental original data, FC should, respectively, be decreased or increased and reduplicated for the next numerical simulation test, until an adequate value of the friction coefficient is found. The value of FC is generally less than or equal to 1.0, except in the case when objects stick to each other. The FC was calibrated using values within 1.0, which represents comparison of the FS3M output value FC = 0.30 (purple dashed line), FC = 0.35 (yellow dashed line), FC = 0.40 (green dashed line), and FC = 0.45 (orange dashed line) with Configuration 11900 (hydraulic experiment) original data (blue solid line). Here, t is the time from the beginning of the container motion. From Figure 9 (blue solid line and yellow dashed line), it can be observed that the maximum difference between Configuration 11900’s (hydraulic experiment) center-of-gravity displacement position and the 0.40 FC numerical simulation center-of-gravity displacement position between 0.10 s and 2.00 s is no more than approximately 4 cm. In addition, we can detect the error range in the x-direction center-of-gravity displacement positions of the FS3M numerical simulation via the root mean square error (RMSE) calculated for FC of 0.30, 0.35, 0.40, and 0.45 from t = 0.00–2.00 s, which are, respectively, 7.03 cm, 3.85 cm, 3.01 cm, and 5.30 cm. All the evidence suggests that an FC value of 0.4 is likely the best value for FS3M numerical simulation. This indicates that the theoretical static friction coefficient between the container model and the chloroethylene false slab surface was approximately 0.4, which resembles the theoretic static friction coefficient between wood and plastic [29].

Figure 10 shows the FS3M numerical simulation Configuration 11190 (yellow dashed line) compared to the corresponding hydraulic experimental Configuration 11900 (blue solid line) for the average center-of-gravity drift in the x-direction. It also shows the FS3M numerical simulation Configuration 22190 (orange dashed line) compared to the corresponding hydraulic experimental Configuration 22901 (green solid line) for the average center-of-gravity drift in the x-direction. Figure 10 (yellow dashed line and blue solid line) reveals that there is only a maximum difference of 4 cm in the center-of-gravity displacement position between the numerical simulation and the corresponding hydraulic experiment values for Configuration 11900. The average center-of-gravity flow trajectory, the magnitude of the flow motion, and the time of the drift velocity slowdown occurrence are similar, which demonstrates the validity and computational accuracy of the FS3M numerical simulation model. Although there is only a maximum difference of 4 cm in the center-of-gravity displacement position between Configurations 22190 (numerical simulation) and 22901 (hydraulic experiment) before 1.60 s, a maximum difference of 10 cm occurs after 1.60–2.00 s. This is because the container gradually rotates around the z axis and moves in the y direction during the drift process, and the projected area subjected to the tsunami force decreases; thus, there is a difference in the x-direction motion and a resulting decrease in the drift flow motion of the container. The variation in the drift flow motion of the container is probably due to irregularities of the chloroethylene slabs assembled in the concrete wave flume. Although the joints of the chloroethylene slabs were covered with silicone, it was not possible to verify in the hydraulic experiments whether water intrusion from the joints occurred to cause a variation in the friction coefficient.

In summary, the FS3M numerical simulation model provided a more ideal experimental state as well as increased computational accuracy, which can effectively guide future numerical simulations describing tsunami-debris motion.

4.2. Containers Center-of-Gravity Drift Behavior

Figure 11 compares the location of the center of gravity in the x direction for the entire container group (increasing number of N_y container specimens with the same initial arrangements of N_x and N_z) at 0.1 s after starting the drift: (a) N_x = 1, N_z = 1; (b) N_x = 2, N_z = 1; (c) N_x = 3, N_z = 1; and (d) N_x = 4, N_z = 1. The layout of N_y in each group increases in the order of 1, 2, 3, and 4. Figure 11a–d shows the motion of the mean center of gravity of the containers in the x direction. From Figure 11a, it can be observed that when N_x is 1, the flow behavior of the center of gravity of the whole group in the x direction does not change significantly when N_y increases from 1 to 4. This phenomenon trend corresponds to the hydraulic experimental results of Tomita et al. [18]. Figure 11a–d show that as the number of N_x container pieces increases, the flow behavior of the center of gravity in the x direction will decrease for the entire group. This is because the increase in the number of containers causes them to be in close proximity and alters the local flow field.

Figure 12 compares the location of the center of gravity in the x direction of the entire container group (increasing number of N_x container specimens with the same initial arrangements of N_y and N_z) at 0.1 s after starting the drift. Figure 13 shows the eight container configurations for N_x = 1 and 4 and N_y = 1, 2, 3, and 4, to compare the conclusions more clearly. The purple and blue dashed lines represent the position of the center of the entire container group at 1.60 s and 2.00 s after the start of drift. When N_y = 1 and 2 (Figure 13a,b,e,f), the increase in N_x will hinder the drift of the center of gravity in the x direction of the whole group; this trend corresponds to the hydraulic experimental results of Tomita et al. [18]. However, when N_y = 3 and 4 (Figure 13c,d,g,h), the increase in N_x will still hinder the drift of the center of gravity in the x direction of the whole group, but this is not as apparent as when N_y = 1 and 2. As shown in Figure 12, it can be observed that the increase in the total number of containers arranged does not slow down the drift of the containers in the x direction; further, as N_y increases, the drift motion of the center of gravity of the entire group will accelerate. This is influenced by columns 1 and 4. From comparing Figure 12c,d, it appears that the trajectories of the average centroid are different; however, the container drift velocities in magnitude and flow motion change time are similar. Figure 12c,d outline a specific characteristic container-transport mechanism. The discrepancy between the container-transport trajectories is observable, as the increase in N_y results in the rotation of the shipping container and the container lateral direction is a more significant deviation. Previous research using an identical shipping container category showed that the debris counterpoise orientation was with the long-axis vertical to the direction of wave flow [30]. Regarding the containers originally arranged through the long-axis vertical to the wave-flow direction, the number of N_y column containers is small (N_y = 1 and 2) and the drift of containers is not subject to local flow-field variations. This results in almost no large-scale rotation of the debris, even if the N_x column increases, with comparatively consistent forces acting on the long-axis section of the containers facing the onshore wave flow. Nevertheless, in the case where the long axis of the containers is placed perpendicular to the wave-flow direction, the number of N_y column containers (N_y = 3 and 4) is more significant, with the increased number of N_y column containers issued in the containers near the open boundary being influenced by the intermediate column containers and the containers rotating. This generates heterogeneous distribution forces acting on the long-axis section of the containers, consequently forcing the containers to diverge from the center axis of the numerical simulation area.

Figure 14 compares the location of the center of gravity in the x direction for the entire container group (increasing number of N_y container specimens with the same initial arrangements of N_x and N_z) at 0.1 s after starting the drift. It shows that an increased number of container layers (N_z = 2) corresponds to a decrease in the wave-flow behavior of the center of gravity of the whole group in the x direction, as well as smaller drifting flow motion than when N_z = 1. For trajectories observed in Figure 14b–d, a reduction of approximately half of the longitudinal displacement of the container was measured in comparison with the one-layer cases of Figure 11b–d. Aside from the increase in the total container weight, the heavier container stack increases the force area, which accounts for the decrease in total container arrangement motion. However, as N_y increases, the drift flow motion of the center of gravity of the entire group will accelerate (except for N_x = 1). This is because column containers near the open boundary are influenced by intermediate column containers. In general, an increase in N_x and N_z can effectively slow down the container drift flow motion and can potentially even do so in advance.

4.3. Comparison of Container Position with Velocity

This section presents three systematic analyses and comparisons of mean velocity among the container configuration group.

Figure 15 compares different number groups of container specimens of the mean velocity (x direction) X_av in the leading-edge flow. In general, in seven container specimen configuration cases, the rapid acceleration of the container can be captured in the initial flow phase. This is because of the horizontal force caused by the incidental wave that the static friction force overcomes, which causes the container to accelerate promptly and reach its peak velocity. When the container configuration N_x = 1 is small (Configurations 11190, 14190, and 14290), the centroid of mass can reach an advanced peak velocity, with the peak occurring faster, at approximately 0.30 s (Configuration 14290’s peak velocity occurs at approximately 0.60 s) and maintaining a more stable and uniform velocity and continuously flowing motion in the x direction. In Figure 11a and Figure 14a, the results are compared with the center-of-gravity positions of the entire container group in the x direction after the start of the drift. The figures show that the center-of-gravity motion positions of the entire container group in Configurations 11190, 14190, and 14290 are more consistent for the same drift time, and the containers drift to the 240 cm position in the x direction at approximately 2.20 s, which is consistent with a stable flow velocity (Figure 15). This indicates that the flow velocity of the whole set of containers is fast and relatively uniformly spread for N_x = 1. In contrast, for larger container configurations with N_x = 4 (Configurations 41190, 41290, 44190, and 44290), the increase in container configuration N_x causes the proximity of the containers to change the local flow field, and the velocity slows down significantly compared to the container configuration with N_x = 1, with the peak velocity occurring at approximately 0.60 s, after which the container flow motion decelerates significantly more rapidly than N_x = 1.

However, comparing Configurations 41190 and 44190 in Figure 15, although the velocities are basically the same at 2.40 s, compared over 0.40–2.20 s, Configuration 44190’s velocities are significantly higher than those of Configuration 41190. By interpreting Figure 11a,d, at 2.20 s, the motion positions of Configurations 41190 and 44190 in the x direction do not change much, which causes their lower flow velocity. At the same time, comparing Configurations 41190 and 44190, the motion position of 44190 in the x direction is significantly larger than that of 41190. This is because the increase in the N_y container causes the intermediate row to affect the container near the open boundary, and the change in the local flow field leads to an increase in N_x, which does not adequately result in the container aggregation increasing the force area and the flow velocity. Furthermore, we can conclude by comparing the total mass of different container configurations that when N_y is small, the velocity of Configuration 11190 (one container) is greater than that of Configuration 41190 (four containers), and the velocity of Configuration 41190 (four containers) is greater than that of Configuration 41290 (eight containers). It can be concluded that, while giving priority to N_y container configuration, the larger total mass of the container configuration will slow down the container flow velocity. It can moreover be observed that the velocity of Configuration 14190 (four containers) is greater than that of Configuration 44190 (16 containers), and that the velocity of Configuration 44190 (16 containers) is greater than that of Configuration 44290 (32 containers). This phenomenon shows that the overall velocity is related to the total mass of the container configuration group; the more significant the mass is, the slower the velocity. In addition, the velocity of Configuration 14290 (eight containers) is greater than that of Configuration 14190 (four containers); the result shows that an increase in N_z container configuration would slow down the container velocity in the case of larger N_y container configurations. This is because containers placed in stacked layers can easily topple over during flow and cause pivoting to initiate a torque around the center of gravity. The torque alters the local flow field of the container motion and causes the container to deviate from its previous trajectory, thereby slowing the container velocity.

In particular, Configuration 41290 becomes almost stationary (velocity = almost 0.00 m/s) during the start flow time of 1.60–1.80 s. Figure 14d shows the location of the center of gravity of the entire container group in the x direction after the start of the drift for Configuration 41290 between 1.60–1.80 s, with the x-direction motion position remaining at x = 62.45 cm, z = 0.10 cm, which is because the horizontal force caused by the incident wave does not overcome the static friction and causes the containers to rapidly reach an almost stationary state.

4.4. Influence of the Center of Gravity on Each Row and Column of Entire Container Group

This section concentrates on the influence of the center of gravity on each row and the column of the whole container group. The average center-of-gravity positions of each row and column of the layout are obtained and compared among different arrangement patterns.

Figure 16 shows the average center-of-gravity positions X1, X2, X3, and X4 for columns 1, 2, 3, and 4 (from offshore side to onshore side starting from x = 0.00 m), and the average center-of-gravity positions Y1, Y2, Y3, and Y4 for rows 1, 2, 3, and 4 (from positive y value to negative y value starting from N_y1 = 0.294 m, N_y2 = 0.098 m, N_y3 = −0.098 m, N_y4 = −0.294 m; N_y1, N_y2, N_y3, and N_y4 represent the maximum y value of the initial container area, and the y position is symmetric with respect to y = 0.00 m) for four different container arrangement groups with Configurations 31190, 33190, 13190, and 11190 (numerical simulation). From Figure 16, it can be seen that an increase in N_x causes the entire group’s center of gravity to move in the x direction and the container motion to slow down. The values of the difference between the respective X1, X2, and X3 of Configurations 31190 and 33190 increased with time after t = 0.40 s when the containers started to move, but the values of the differences between the three were almost equal. This means that column 1 is affected by column 2, and column 2 is greatly affected by column 3. Therefore, the motion of the containers in column 1 is indirectly influenced by the containers in column 3. This indicates that there are interfering effects between containers of the same configuration. This trend corresponds to the hydraulic experimental results of Tomita et al. [18].

Figure 17 shows three different container arrangement groups with Configurations 33190 and 13190 (numerical simulation). The values of the difference between the respective Y1, Y2, and Y3 of Configurations 13190 and 33190 increased with time after t = 0.30 s when the containers started to move, but the different values between Y2 and Y1, Y2, and Y3 are larger than those between Y1 and Y3. This means that the flow motion of the containers in rows 1 and 3 is directly influenced by the motion of the containers in row 2. In addition, the flow motion of the containers in rows 1 and 3 of Configuration 13190 is essentially the same after being disturbed by the flow motion of the containers in the second intermediate row. However, the container motion of row 1 and row 3 of Configuration 33190 is slightly larger than that of row 3 after the container motion disturbance of the second intermediate row container. This is because when N_y = 3, row 2 will cause uneven force distribution in the long-axis cross-section of row 1 and row 3, resulting in inconsistent motion and gradual deviation from the central axis of the numerical simulation area while the container is rotating.

Figure 18 shows four different container arrangement groups with Configurations 34290, 34190, 31290, and 31190 (numerical simulation). The values of the difference between the respective X1, X2, and X3 are almost equal. Comparing X1, X2, and X3 for Configurations 34290 and 34190, and comparing X1, X2, and X3 for Configurations 31290 and 31190, the following can be summarized. When N_x and N_y container configurations are the same, increasing N_z to 2 effectively slows down the container group motion. However, this does not mean that an increase in the total mass of the container configuration group will slow down the motion of the configuration group; for example, comparing all 24 containers of Configuration 34290 with the three containers of 31190, their container drift flow to motion is similar, which means that the increase in mass has no effect on the motion of the container configuration.

4.5. Spatial Container Distribution

This section compares container spreading angles and estimates the potential tsunami–container impact and whether the number and layout of the initial container configuration are key variables affecting container motion.

Figure 19 shows the trajectories of flow for multiple container groups, where the 32 container specimen configurations for the numerical simulation tests are shown in 12 groups. They are Configurations (a) 14290; (b) 41290; (c) 22190; (d) 33190; (e) 43190; and (f) 44190. Figure 19 illustrates a hollow quadrilateral box representing the initial areas where the containers were set initially, to provide enough distance reference for the wave deformation around the container and keep the entire footmark of the initial container domain constant. In addition, the overall center-of-gravity point from the initial position of the container configuration is the center (focal point), which is the convergence point of all meta-containers and performs the spreading angle division. Within each panel, the container trajectory for each of the classifications cataloged in Table 1 is contrapositive to the spreading angles (two dash-dotted dark lines for Figure 19 reference) devised by Naito et al. [16], which was suggested as an unsafe region based on the domain survey assessments from two different topographical regions conducted after the 2011 Japan Tohoku earthquake and tsunami:

$$\theta =\pm 22.5°.$$

(1)

Nistor et al. [28] focused on repeated physical experimental research; they considered the spreading angle’s period history and the maximum spreading angle per debris unit used in each physical experiment results. They determined the unsafe region at which debris spreading increases as the quantity of debris components increases. The results revealed that

$$\pm \theta =\pm 3.69°\pm 0.80°\times N$$

(2)

where N denotes the number of container specimen configurations. In general, without discussing N_x and N_y, for the container configuration with N_z = 1, or the container configuration with a smaller N_y but N_z = 2 (Configurations 22190, 33190, 43190, 44190, and 41290), the experimental numerical simulation of container flow motion results in a narrower spreading angle than that proposed by Naito et al. [16]. This is because this experimental numerical simulation essentially presents an idealized simulation case without considering the structures during flow and without topographic features, such as slopes. However, in the case of larger N_y with N_z = 2 (Configuration 14290), the container propagates beyond the exposed disseminating region proposed by Naito et al. [16]. This suggests that the exposed spreading region equation proposed by Naito et al. [16] is conservative in estimating the container spreading angle for the group of container specimen configurations with larger N_y with N_z = 2.

Correspondingly, in Equation (2), the container spreading angle increases equally as the quantity of container specimen configurations increases. However, observing Figure 19a,d–f, the spreading angles of Configurations 14290 are greater than those of Configurations 33190, 43190, and 44190. Thus, the increase in the number of container configurations is not the only condition leading to the increase in the spreading angle, and the exposed spreading region equation proposed by Nistor et al. [28] seems to be sufficiently conservative for Equation (2). In addition to the increase in the number of container configurations, certain basic observations highlighted that it is the difference in the number of initial container configurations N_x, N_y, and N_z that is the main factor affecting the container spreading dynamics. Figure 19 shows the initial number of container arrangements for different N_x, N_y, and N_z. The initial container arrangement N_z of 2 (containers stacked in multiple layers) with N_y ≥ 2 (Configuration 14290) increases the mass for the same contact surface. This will lead to a tremendous initial friction force that must be overcome by hydrodynamic forces to initiate the movement of the container. Nonetheless, once the friction is overcome, the shipping container displacement in the stacked layers topples over, and their pivoting institutes a torque around the center of gravity of the container. This in turn causes the container to deviate from its antecedent trajectory vector to flow over a broader spreading region.

5. Conclusions

In this study, the interaction between containers within extreme hydrodynamic conditions was investigated using a numerical simulation model. The method to improve the assessment technique of the motion behavior of tsunami drifters and the initial arrangement of containers in a tsunami-prone site was examined and discussed in terms of the numerical test model using the three-dimensional coupled Fluid–Structure–Sediment–Seabed Interaction Model (FS3M) developed by Nakamura et al. [22]. This numerical simulation result was compared with the hydraulic experimental data [18] from the specific target of this study to check the validity and computational accuracy of the FS3M. A framework was proposed to help calibrate and validate the FS3M to guide future numerical simulation experiments describing tsunami-debris motion. This research involved considerable data processing and analysis. Particular consideration was paid to the complex movement of multiple container flows. Based on this work, the major contributions and conclusions are summarized as follows:

• In the comparison of numerical simulation results and hydraulic experimental data, it can be seen that the results and data provide evidence of the validity and computational accuracy of the FS3M numerical simulation model. Furthermore, the FS3M numerical simulation model can provide a more ideal experimental state as well as increased computational accuracy; thus, it can effectively guide future studies of numerical simulations describing tsunami-debris motion;
• The increase in N_x (number of drifter containers in rows) and N_z (stacking height) can effectively slow down the container drift flow motion and can potentially do so in advance—even to the point of a noticeable pause. However, the decrease in N_y (number of drifter containers in columns) in containers of the same N_x configuration significantly slows down the container motion of the container configuration arrangement group in the x direction. However, the overall velocity is related to the total mass of the container configuration group, with more significant mass correlated to lower velocity, whereas priority needs to be given to N_y container configuration;
• The increase or decrease in the total mass of the container configuration group does not have any pivotal effect on the container motion. In general, the initial number of N_x, N_y, and N_z container configurations is the crucial factor affecting container group motion;
• The exposed spreading region equation put forward by Naito et al. [16] is conservative in estimating the container spreading angle with larger N_y and N_z = 2. The increase in the trend of the container configurations number is not the only condition leading to the increase in the spreading angle. The exposed spreading region equation proposed by Nistor et al. [28] is also conservative. In summary, the difference in the number of initial container configurations N_x, N_y, and N_z is the main factor affecting the container spreading dynamics.

In general, this paper highlighted the importance of differences in the number of initial container configurations as affecting the dynamics of container spreading. However, this study focused on limited conditions; for example, only one wave condition and only one container weight (empty container). Future work should consider investigating the influence of multiple debris arrangement groups, and debris fields carried by different tsunami waves through hydraulic experiments and numerical simulations, and propose characteristic damage estimation and comprehensive assessment methods for more realistic structural situations. Further improvement of the computational accuracy of the FS3M is also one of the issues to be solved.