Tsunami simulation using particle method

Tsunami is a natural disaster that have resulted in dreadful damages over time. Extensive researches have been conducted to scrutinize and counteract the natural hazard using three major research components which are: field monitoring, laboratory tests, and numerical methods. However, laboratory tests are high-priced and arduous. Numerical simulation overcomes these drawbacks and can be utilized in collaboration with laboratory tests. Recently, newly introduced meshless Lagrangian particle method called Smoothed Particle Hydrodynamics (SPH) has gained attention. In this paper, SPH method has been employed to simulate tsunami. A SPH code is developed from scratch. To validate the code, a traditional dam break simulation is conducted. Lastly, a tsunami model is simulated using the developed SPH code and compared with past experimental data. The results indicate that the code is in accordance with previous experimental data and numerical simulation. Whereby, there’s been a slight deviation arises in tsunami simulation. The velocity of the code is relatively less to that of the experimental data. Such inconsistencies could emerge due to a number of reasons, i.e. the choice of the SPH parameters and model simplification. Generally, the developed SPH code had a satisfactory performance to model tsunami and dam-break problem.


Introduction
On Sunday, 26 December 2004, almost 230,000 people lost their lives within several hours due to the Indian Ocean tsunami [1][2][3].The tsunami disaster hit 15 Asian and African countries [2].It had the most intense impacts in Indonesia, Aceh and Nias.About hundreds of kilometres of the coastline and thousands of hectares of land were wrecked [4].The aggregate financial damage was estimated to be about USD 4.5 billion [4,5].On 17 July 2006, tsunami killed approximately 802 people were either killed or declared missing as a result of a tsunami that occurred in Tasikmalaya, Indonesia.The disaster was triggered by a 7.7 Mw earthquake [6].
With respect to historical records, Indonesia is listed as one of the top two countries in Asia experiencing the most tsunami disaster after Japan [2,6].The two regions in Asia facing the highest tsunami risk are South-eastern and Southern Asia [2].In the last four decades, tsunami disasters in Indonesia have escalated by a massive 80 % [6].
Tsunami is a phenomenon related to a rapid displaced water column [7].A number of mechanisms have the capability of triggering a tsunami, and they are: earthquake, submarine landslide, and meteorite impact [7].In Indonesia, 88% of tsunami disasters have originated from earthquakes [6].Furthermore, Indonesia is named as one of the countries with the highest seismic activity [3].During 2004 -2007, six devastating earthquakes befell Indonesia [3].
Attributable to the recent tsunami phenomenon, researches about tsunami have dramatically escalated in the last decade [7].The foremost issues of the tsunami research are numerical and experimental modelling, field investigation of tsunami impacts, and tsunamigenic research [7].In general terms, the laboratory tests are usually expensive and difficult [8].Numerical simulation can work together with this laboratory test to overcome the drawbacks.In hydraulic engineering, Finite Different Method (FDM) is favoured, but several problems are encountered when dealing with rapid change water level in time and space, breaking wave, etc. [8].The recently introduced Smoothed Particle Hydrodynamics (SPH) has done well to overcome these challenges.
SPH is a mesh-less particle method that has proven to solve numerous engineering complications.In this paper, the tsunami phenomenon is simulated using SPH.By means of C++ programming language, a code to resolve hydrodynamic problems is developed.The dam-break simulation is utilized to validate the code.Lastly, the developed code is employed to simulate a tsunami phenomenon.The result is set side by side with previously collated experimental data of tsunami from a different researcher to put the accuracy of the SPH code to a test.

Smoothed particle hydrodynamics (SPH)
Smoothed Particle Hydrodynamics (SPH) is a mesh-free lagrangian method to approximate solutions of equations of fluid dynamics by discretizing continuum material into a set of particles [9][10][11][12].It was initially proposed by [13] and [14] to solve the astrophysics problem.In SPH method, each particle keeps properties of the material, i.e., mass, density, pressure, etc.These properties are brought when the particles move with their velocity [10,11,15].To overcome the inaccuracy of large displacement material, the meshing technique is avoided [10].Material properties of each particle are obtained using the interpolation technique over its neighboring particles [11,15].The interpolation technique utilizes an integral function which is described in several texts [15].The discretized equation of integral function is: where the subscripts a and b denote integration point and neighbouring particles respectively.
The density and mass of the particle are denoted by m, respectively.r ab is the distance of particle a (integration point) and neighbouring particle b and W is the kernel (smoothing function).
In this study, the cubic-spline kernel function proposed by [16] is used: where  d is 10/(7h 2 ) and 1/(h 3 ) in 2D and 3D respectively.
Conservation of mass using artificial viscosity [17] and equation of mass can be written in SPH equation as: where P is pressure, g is the gravitational acceleration, and v ab = v a -v b . ab is viscosity term obtained in [17].As fluid is considered weakly compressible, an equation of state can be used to calculate the fluid pressure: where = 7, � � � � �  � /,  0 = 1000 kg/m 3 .
Particles move using averaged neighbouring particles called XSPH: where  = 0.5 and Eq. ( 3) and ( 4) are integrated using the Leap-Frog scheme:

Programming
A SPH code called UNSph (Universitas Sebelas Maret Smoothed Particle Hydrodynamics) code is developed from scratch making use of C++ programming language.To accelerate calculation processes, it employs multithreading programming.Linked list algorithm is implemented for neighbour searching.Therefore, the computation domain is divided into square cells with the width of 2h.Each particle in a cell only forms an interaction with other particles within the confines of neighbouring cells.Time-integration is implemented using the Leap-Frog scheme, whereas the time stepping complies with the CFL (Courant-Friedrichs-Lewy) conditions.

Breaking dam
UNSph code is validated using a dam-break simulation.A dam-break phenomenon is described as the surge of fluid after a dam is broken [8].Similar to the past SPH researches [8,17,18], the research also simulates the collapsed column of water after the dam is suddenly removed.The result is compared with [18,19] to investigate the performance of the UNSph code.https://doi.org/10.1051/matecconf/201819505013ICRMCE 2018 The set-up of the breaking dam simulation in the top-opened tank is shown in Fig. 1.The initial width and height of the water column are a and 2 respectively with a = 0.146 m.The water column with the length of 4 is placed in the tank.The simulation condition is adopted from [18], except that the height of vertical wall of the tank is 4, while in [18] is 2.5.At � � �, the water is allowed to flow.When the particles of water exceed the top boundary, the particles are prevented to re-entering the tank.The particle spacing of 0.002 m, � � ���� kg/m 3 , � � ��.5� m/s, � � �.2� are employed.

Tsunami simulation
There are many experimental modelling methods to generate long waves, e.g.piston-type wave, dam-break, and a constant flow pump.A piston-type wave-maker is a horizontally moving mechanical paddle, which is driven in piston-like motions [20][21][22].A constant flow pump wave generator is a wave maker employing a pump to discharge a constant flow, which produces a gradual increase of water level [21].The dam break-induced tsunami is generated when a tight-door of the column of the water is suddenly removed.
In this research, a dam-break analogy for wave generation is picked to simulate tsunami.The schematic figure of tsunami basin utilized in the numerical analysis is shown in the Fig. 2. The 24.7 m in length and 0.8 m in height of tsunami tank is used.The wall is located in the left boundary, while the right boundary is opened to let the water leave the computational domain.The depth of the water is 0.4 m.The 0.25 m high and 7.9 m wide water column is set up above the 0.4 m depth of the water and attached next to the left boundary.The slope of the basin is 1:40.This set up is adopted from [23].In the original experiment of [23], the length of tsunami basin is 44 m.It has two slopes, i.e. 1:40 and 1:100.Due to the fact that the experimental set up is located in the 1:40 slope section and the computational coast is needed to be reduced, the numerical simulation is only focused on the 1:40 slope region.Thus, the total length of the tsunami basin in the numerical model is 24.7 m.In the simulation, all of the dam-break's particle properties are adopted, except the particle spacing, which is 0.01 m.

Breaking dam
To investigate the performance of the UNSph code, the result obtained from numerical analysis is compared with that of the experimental test [19] and the past numerical analysis carried out in previous times [18].In [19], the experimental tests are conducted using various n 2 and a, whereas in this paper, n 2 = 2.In the initial condition of the experiment, the paper diaphragm is used to prevent the water column from flowing out.Then, the electric current is streamed through the paper to break it and collapse the water column collapsed.The motions of the water at several intervals of times are captured.The water surge profile after the shatter of the dam breaks is illustrated in Fig. 3.At the inceptive stages of the simulation, the water column stands up at the initial condition (Fig. 3.a).Suddenly, the water is allowed to flow to the right boundary, while the left-hand side of the water column is confined by the left wall boundary (Fig. 3.a -d).At about t = 3 s, the surge front smashes into the right wall boundary (Fig. 3.e), as a result, deflecting the surge is deflected (Fig. 3.f).To validate the SPH code, the dimensionless of remaining water column height H and front position of the flow Z are introduced [19]: where a, h, z, and n 2 are the width of the water column, the top of the column, the surge front, and the ration between the initial height and width of the water column, respectively.The dimensionless H and Z with its corresponding normalized time  and T are plotted in Fig. 4. In Fig. 4.a, the result of the UNSph code is compared with the experimental data from [19], as well as the numerical analysis form [18].According to the figure, they are no significant differences between the result of the UNSph code, [19], and [18].The small deviation arises in the result of the surge front (Fig. 4.b).The slope of the graph represents the velocity of the surge front.In the initial normalized time, the slope is zero.The velocity is increased with the increase of the time until it approaches a constant value in T = 1.25 -1.5.In the experimental test, there are two different experimental data plotted, i.e., experimental test using the width of column a = 1.125 inches and 2.25 inches.
At the beginning of the test, the two data are almost similar.The deviation becomes apparent after normalized time T is about 1.5.The velocity of surge front of experimental data with a = 1.125 tends to be higher than another experimental data.The surge front of the UNSph is slightly higher than the other data [18,19] at the beginning of the test.The surge front of the UNSph approaches the value of the past experiment data where a = 1.125 inches data, while it remains to have a small deviation in comparison with a = 2.25 inches experimental data and past numerical analysis.Since there are no significant differences between the results of UNSph, experimental test, and past numerical analysis, the fluid motion is correctly simulated using UNSph code.

Tsunami simulation
The numerical simulation of the tsunami is conducted using the dam-break analogy to generate the wave.A 0.25 m height and 7.9 m width of the water column is permitted to flow.The flow of the breaking column of the water creates a solitary wave with a long period.The wave travels along the tsunami basin with a slope of 1:40 (Fig. 5).The wave then starts to break (Fig. 5.c) and approaches the observation point located 20.5 m from the left boundary (Fig. 5.d).According to the Fig. 5, the UNSph code can simulate the breaking wave phenomenon.It is also found that the UNSph code can simulate the wave run-up correctly.
To validate the simulation, the result is compared with the experimental data [23].Fig. 6 shows the result of tsunami simulation of the UNSph code and the experimental test from [23].Fixed particles having zero velocity represent the wall and the floor.In the UNSph simulation, the generated tsunami wave is over 10 cm high (Fig. 5.d and Fig. 6.a), it means a 15 m high tsunami wave in the full-scale.It has good agreement with the experimental test [23].
The velocity of the experimental test in [23] is measured at 0.02 m below the free surface water in calm condition.This point is separated about 20.5 m far from the left boundary.The maximum velocity recorded in [23] is around 1 m/s.In the UNSph simulation, the velocity of 1 m/s is located near the free surface.This condition affects the arrival time of the surge in the observed point.In [23], the wave reaches the observed point in approximately around 5 -5.5 seconds after the watertight door is opened, while in the current research it is about 6.75 seconds.
The deviations happen due to several reasons, such as errors due to the simplification of the numerical analysis and the selection of the SPH parameters.For example, the numerical analysis ignores the motion of the watertight door when it is opened.In the experiment of [23], when the watertight door is opened, small waves are generated.The accuracy of the SPH method is influenced by the parameters selection of kernel, smoothing length, viscosity equation, background pressure, and speed of sound, etc. [24].Overall, the UNSph has a good agreement with the experimental test of [23], although the velocity of the wave simulated by the UNSph code is smaller than the experimental code.Generally speaking, the UNSph code has an excellent performance to model a tsunami and dam-break problem.

Conclusions
The recently emerged numerical method overcoming the large displacement problems is implemented from the scratch and named UNSph.The code is validated with the use of breaking dam and long period tsunami wave problem.During the dam-break simulation, the UNSph code has correctly aligned with the compared previous experiment and numerical analysis.The surge motion of water and the top of remaining water column are simulated without a significant deviation from the preceding studies.The tsunami simulation is also in accordance with the preceding experimental test, although the velocity of the wave simulated by the UNSph code has a small deviation in comparison with the earlier study.This small deviation emerges due to several factors, i.e. the simplification of the numerical analysis and the selection of the SPH parameters.
This research was supported by PNBP UNS 2018 Grant led by Agus P. Saido.

Fig. 1 .
Fig. 1.Set up of the breaking dam simulation.

Fig. 2 .5
Fig. 2. The sketch of the tsunami basin (not to scaled).

Fig. 3 .
Fig. 3. Water profile after the dam is broken.

Fig. 4 .
Fig. 4. Comparison of dam-break simulation using UNSph code and past researches: a) Normalized remaining column height Vs.Normalized Time , b) Normalized surge front Vs.Normalized Time T.