A hybrid numerical model based on FNPT-NS for the estimation of long wave run-up
Introduction
The extreme wave events such as cyclonic storms, coastal storms, and tsunamis pose a severe risk to coastal regions. These events, along with the sea-level rise and climate change, threaten the very existence of communities that live along the coast. Predicting and mitigating these threats is essential for disaster mitigation and management in these regions. It is necessary to understand these phenomena before taking up any developmental activity in the coastal zone.
Tsunamis are extreme waves generated by an earthquake or any other disturbances on the ocean floor or landmass close to the sea. These waves trace its origin to deeper waters and usually travel for very long distances. Some of the characteristics of a tsunami are long wave periods and tremendous wave heights (as they arrive in the nearshore region). The behaviour of tsunami waves changes as they propagate from deep water to shallow waters. Over the years, there has been a significant number of tsunami events recorded in many parts of the world. Our previous encounters with these waves have helped us to understand how they affect the coastal region. All of these make the study of these extreme waves an exciting topic for scientists and engineers.
Although tsunamis are often considered in the framework of the nonlinear shallow-water theory, they may also be affected by dispersive effects, when propagating over a long distance. The latter may result in disintegration of the initial wave into a series of solitons of smaller period and increase in wave steepness due to exchange of energy (see, for instance, Akrish et al., 2016; Brühl and Oumeraci, 2016; Sriram et al., 2016; Ezersky et al., 2013; Herterich and Dias, 2019; Li et al., 2019).
One of the first models proposed for the laboratory study of the tsunami was a solitary wave. A solitary wave has just a crest, which makes them look like an elongated pulse signal. They are defined with a fixed wave period and wave height. These waves were used in Synolakis (1986, 1987). The experimental results were compared to an in-house numerical model. The results were used to develop the widely popular analytical expressions for calculating wave run-up height on a plane slope.
It was often observed that during the onset of a tsunami event, the shoreline recedes. This behaviour is not characteristic of the elongated pulse/solitary wave model. The receding shoreline phenomena were found to cause an increase in the number of causalities during the recent tsunami events. Madsen et al. (2008) explained the solitary wave paradigm for the tsunami. They observed that most of the globally accepted models for tsunamis consist of either solitary wave or cnoidal wave models. They noted that, even though these formulations are useful for validation and verification of the numerical models, the relation between these models and geophysical tsunami was not appropriately established.
The N-wave is another model of tsunami waves (Tadepalli and Synolakis (1994)). In comparison to the solitary wave model, N-wave has a distinctive trough and crest. Based on the relative shape of the crest and trough, the N-wave can be symmetric or asymmetric. In our studies, Manoj Kumar et al. (2017), we observed that a leading trough N-wave model could accurately model the receding shoreline phenomenon seen during some of the field observations. We also noted that the breaking characteristics and effective run-up of both solitary and N-wave models were quite different even when they had the same wave height.
Experiments in the model scale, analytical formulations, and numerical models are some of the most widely used tools when it comes to the study of wave evolution and run-up. Analytical models and empirical models are used for initial assessment of maximum inundation and run-up (Van der Meer et al. (2016); Synolakis (1987)). Schimmels et al. (2016) and Sriram et al. (2016) used a piston-type wavemaker to generate and study the evolution of tsunami in a large-scale experimental facility. The experiments were carried out at the Large Wave Flume (Groβer Wellenkanal, GWK) in Hannover, Germany. The above papers brought to light the importance of modelling long wave propagation in large-scale testing facilities. Further, the run-up and inundation of the generated waves were studied using the analytical model developed in Didenkulova (2009). It was noted that the run-up on a mild slope (<1:18 concerning the study) was challenging to estimate due to the phenomenon of wave breaking during the run-up.
Empirical models have been proposed by the researchers to estimate the run-up of tsunami-like waves. The expressions were developed based on rigorous model studies performed in the laboratory using solitary, cnoidal, and N-wave models. In many of the cases, the results were supported by numerical studies based on the shallow water equations, Boussinesq equations, potential flow theory, and Navier Stokes equations. These estimates help in providing a starting point for researchers when it comes to the study of extreme wave events. The analytical expressions are usually developed as closed-form solutions. During the wave evolution, the chances of wave deformation are higher. The closed-form solutions tend to fail here. The alternate approach would be to use the Navier-Stokes or potential flow models. A good number of models are found to give excellent results.
Over the years, numerical models have gained popularity in studies related to the behaviour of long waves and their propagation characteristics. They provide flexibility to the researchers by simulating many different test cases, which would be too expensive to perform in physical model tests. Mesh-based methods and particle methods are now being extensively used in this area. Smooth Particle Hydrodynamics (SPH) is a popular particle method for simulating the free-surface flows. Gingold and Monaghan (1977) and (Lucy, 1977) developed this method for problems in astrophysics. Since then, the meshless nature of the technique has made it popular among researchers working on a broad spectrum of topics. The original formulation has been improved considerably over the years to accommodate the needs of specific problem requirements. Gotoh et al. (2004) used this method to investigate the interaction between regular waves and partially immersed breakwater. The study considered both overtopping and non-overtopping cases. Large Eddy Simulation (LES) was used in the formulation to account for the viscous effects. A similar approach was used in Rogers and Dalrymple (2005) to model breaking waves. In both these studies, the test cases include 2D and 3D models. Dalrymple et al. (2006) presented the challenges posed by different mathematical approaches when used for modelling the tsunamis. Emphasis was given to the offshore and onshore phenomenon caused by tsunamis. The authors commented that no models were able to provide a real-time model of the 2004 Tsunami event.
Crespo et al. (2007) used the SPH-LES formulation by Rogers and Dalrymple (2005) to study the effect of large waves on coastal structures in the presence of a dike. Crespo et al. (2008) used the SPH technique to simulate the dam break problem. Instead of using the modified formulation for viscosity, the artificial viscosity formulation proposed by Monaghan (1992) was used. Dalrymple et al. (2010) compiled the developments in SPH over the years in simulating water waves. Special consideration was given to different formulation and solution techniques to obtain a stable solution. Rogers and Dalrymple (2010) used the SPH method to model Tsunami waves. The study was based on the model proposed in Rogers and Dalrymple (2005). Capone et al. (2010) used SPH to simulate the waves caused by submarine landslides. The governing equation used was based on the rheological non-Newtonian Bingham model. The model was able to represent the landslide deformation and its interaction with water. The SPH formulation based on the artificial viscosity term was used.
In all these studies, the researchers have noted that the 3D simulations need a large number of particles. The compressible fluid model affects the spatial and temporal resolution requirements. Another observation was that the models were good at simulating relatively smaller regions with a reasonable degree of accuracy and not suitable for large areas. The high resolution is directly proportional to computational costs. Parallelization and hybridization were proposed as solutions. GPU-SPH (Hérault, 2010) was introduced as a faster alternative that can accommodate a large number of particles and still be computationally effective. Farahani et al. (2013); Wei et al. (2016, 2015); Wei and Dalrymple (2016) have effectively used this model. The studies were able to reproduce phenomena such as mean wave-induced nearshore circulation created by a breaking wave, wave plunging, collapsing, and run-up and tsunami forces on bridges, respectively. Gotoh and Khayyer (2018) have presented a detailed review of particle methods. Even though particle methods are good alternatives to mesh-based methods, most of the popular techniques are computationally expensive even for small wave steepness simulations. A feasible option is the use of hybrid models.
A hybrid model combines two or more models and solution techniques to build a robust method to solve a single physical problem. They combine the advantage of different models and help to use the best possible models for a physical problem. The solution and coupling strategy adopted has to be robust to make sure that the domains communicate effectively without sabotaging the integrity of the flow problem. Depending on how the domains are coupled at the interface, we have two specific classes of solution approaches: weak coupled approach and strongly coupled approach. In the weakly coupled method, the information exchange at the interface occurs only in one direction. The domain and mathematical models are independent of each other. There is no feedback. Whereas in the strong coupled approach, there is information exchange at the interface between the two models in both directions at every time step. The domains are coupled in both space and time.
The flexibility and higher computational costs associated with the particle methods make them an ideal candidate for the hybrid models. There have been many instances in recent years where hybrid models have been used to study hydrodynamic and coastal engineering problems. SPH (Smooth Particle Hydrodynamics) is being widely used in hybrid models. Altomare et al. (2018); Mogan et al. (2018); Napoli et al. (2016); Narayanaswamy et al. (2010) are some of the recent examples. The mesh-based method was combined with SPH (along with parallel computing and CUDA capabilities) to solve wave problems. The authors have claimed improvement in total simulation duration. Physical problems such as wave overtopping and standing waves were considered in these studies.
In the present paper, the experimental data from the Large Wave Flume (Groβer Wellenkanal, GWK), Hannover, Germany, is used as an input to the numerical model to simulate wave propagation and its run-up characteristics. The numerical method employed is based on Sriram et al. (2014). This method combines a potential flow theory based on the Finite Element Method and the Navier-Stokes equations based on Meshless Local Petrov-Galerkin Method.
The paper is organized in the following manner. Firstly, the mathematical formulation and a brief explanation of the solution technique used are provided. Secondly, validation of the model against experimental data of long wave propagation and run-up from large scale experiments is reported. Thirdly, a comparison of the small-amplitude wave run-up (representing the rapidly raising tide) with the analytical model is reported. And finally, the simulations for the undular bore formation and its plunging breaking over the slope are presented.
The strongly coupled 2D numerical model is created by combining two subdomains. One based on the fully nonlinear potential flow theory (FNPT) solved using the Finite Element Method (FEM) and the second based on the Navier-Stokes (NS) flow model solved by using Meshless Local Petrov-Galerkin Method based on Rankine source solution (MLPG_R) method. The detailed description and discussion on this method and formulation can be found in Sriram et al. (2014). A brief description of each formulation is given below, along with the coupling strategy adopted.
The fluid is assumed to be incompressible, inviscid, and flow irrotational. The flow is defined in terms of the velocity potential using the Laplace equation:
The boundary condition at the bottom impermeable boundary (flat - bottom) isThe wave paddle motion is prescribed byHere is the time history of paddle motion at . denotes the wave paddle boundary. The nonlinear dynamic free surface boundary condition in the Lagrangian system is given byAt the free surface, pressure, is set to zero and , where is the water surface oscillations. The potential inside an element can be expressed in terms of nodal potential, as is the shape function, and is the number of nodes in an element. The finite element formulation for the system is given in the following form, is the total number of nodes in the domain, is the free surface boundary. A detailed account of solution procedures can be found in Sriram (2008).
The governing equations consist of the continuity equation and momentum equations, as given below: is the gravitational acceleration, is the fluid velocity vector, is the pressure, is the density of the fluid, and is the kinematic viscosity. The solution procedure is supported by the boundary conditions specified in the problem domain. The Lagrangian form of kinematic and dynamic boundary conditions on the free surface is given by
In Eq. (10), denotes the position vector for the particle position. The boundary condition on the rigid boundary is given by denotes the unit normal vector from the rigid boundary, and are velocity and acceleration of the solid boundary. The bottom boundary was specified using a slip boundary condition, neglecting the boundary layer effects.
To evaluate the system of equation, we use the Chorin's time split algorithm given in Chorin (1968). The solution procedure starts from the stage where the state values are all known. The values at state will be estimated based on the algorithm. First, we evaluate the intermediate velocity () and position () based on the values at :Here, in Eq. (14), the term is determined based on Koshizuka et al. (1998). The pressure, is estimated using the semi-implicit Poisson's equation given by Ma and Zhou (2009). here is the artificial coefficient with values between 0 and 1. We use the value of = 0 (Sriram and Ma (2012)). This gives us the equation below.The boundary conditions on the surface and solid boundaries were set as and respectively. The estimation makes use of the MLPG_R. The local weak formulation for a circular support domain is given byIn Eq. (18), is unit normal pointing outside of the subdomain, is the solution for Rankine source in an unbounded 2D domain with , the distance between point under consideration and centre of local subdomain .
The pressure is estimated from the local weak formulation (Eq. (18)). The velocity is then calculated based on the pressure from the pressure Poisson equation. Based on the velocity and the initial position of the particle, the particle is moved to a new position. The calculation is carried out for all the particles in the support domain. The coupling technique adopted, and the algorithm employed is discussed in the next section.
Two methods to achieve coupling in the hybrid models are weak coupling and strong coupling. In weakly coupled models, the domains are independent of each other, except for the boundaries. The information is exchanged from one model to the other without any feedback. This approach allows the simulation to be run stand-alone. In a strongly coupled model, the domains communicate with each other at the interface during all time instances. The boundary information is shared at the interface. The strong coupling has to be implemented in two ways (i) in space and (ii) in time. For coupling in space, we use moving overlapping zone method. Fig. 1 below shows the concept.
In this method, boundaries can move as the simulation advances. The line denotes the boundary for the FNPT domain and always moves along with the free surface and is kept straight (adopting quasi-Eulerian and Lagrangian approach). The line represents the boundary for the NS domain, and it can change shape. Here and are lengths of the fluid domain and the overlapping zone respectively at time . As the simulation advances to the time , the value changes to and . The next important step is to understand the treatment of different boundaries during the simulation.
In effect, there are three boundaries, as described in Fig. 2. The boundary marks the start of the NS domain, boundary , marks the end of FNPT domain, and boundary corresponds to the region where the feed particle starts. Special treatments were used to incorporate and extract these values into the simulation. There are two distinct formulation and solution methods for the two different domains here. The approach adopted is described in the algorithm given in Fig. 3. The strong coupling in space for the fluid-fluid subsystem is based on Sriram et al. (2014). The algorithm operates from the time step , where all the state variables are known. The intermediate values of pressure, is estimated by solving the boundary value problem for . The estimation is done based on the values and . Where, and denotes velocity potential and position, respectively. From the intermediate value of potential, the intermediate value of pressure, on the boundary is estimated. (Section B, Fig. 3). is the boundary (on which the values are calculated) of the MLPG domain and denotes that the operation has been carried out in the FEM domain. Based on the boundary values of pressures obtained from the simulations carried out in the potential flow domain, the boundary value problem for pressure was solved in the Navier Stokes domain using the MLPG_R method. The boundary value on the domain is updated based on the estimates from the MLPG_R (Section A, Fig. 3).
To ensure that the pressure gradient is evaluated correctly on the boundary , feeding particles are used. For calculating the pressure gradient, the Simplified Finite Difference scheme (SFDI) given in Ma (2008) is used. The estimation based on SFDI depends on the local support domain for each particle. If a particle is on the boundary , it doesn't have enough number of particles for calculation. Additional particles were introduced between the boundary and to address this issue. In the simulations presented here, ten layers of feeding particles were used.
The FNPT and NS domains share the region between boundaries and . This region is called the two-layered particle (TLP) zone. Since two different solvers are used, the field variables estimated in this region will be similar but not the same. Using two solvers on the same domain create continuity issues along the boundaries and affects the stability of the solution process. A smoothing technique was employed to address this problem. The smoothing function uses a location-based approach to select the values from both the domains and applies weights to smooth the velocity in the region shared between the two solvers. The technique has been tested and found to work well.
Another critical factor to be kept in mind here is that there are two different solvers using two different time integration approaches; hence, coupling in time needs to be carried out. The FEM uses explicit time integration while the MLPG_R uses fully implicit time integration. To connect both the models, RK2 predictor-corrector approach was employed. The process is illustrated in the flow chart given in Fig. 3.
The operation starts from a physical state at a time, , where the free surface elevation and velocity potential are known. The Runge-Kutta 2nd order time integration technique was employed to find out the velocity potential and new position for the next time step using Eqs. (4), (5), (6) and Eq. (19a) –(19d):Here the values prefixed with in the superscript denote intermediate values of position and velocity potential. A detailed description of the procedures and validation cases of the hybrid model can be found in Sriram et al. (2014).
Section snippets
Results and discussion
In Sriram et al. (2016), three different scenarios for the tsunami approach based on Shuto (1985) classification was considered.
- 1.
Non-breaking wave surging that act as a rapidly rising tide, observed during small and moderate tsunami events after short distance propagation;
- 2.
Breaking bore or hydraulic jump (“wall of water”), observed as a result of wave breaking during a large tsunami event after short distance propagation;
- 3.
Undular bore, observed after the long-distance propagation of tsunami (in
Conclusion
The estimation of wave run-up has a critical role in any infrastructure development projects in proximity to the coastal zones. In the present study, a 2D hybrid numerical model was used to estimate the run-up of tsunami-like waves and to understand their propagation and evolution over a long distance. The method was able to simulate the propagation of long elongated pulse waves and split-up for small and steep waves. Good agreement with the laboratory measurements is obtained. Further, the
CRediT authorship contribution statement
G. Manoj Kumar: Conceptualization, Formal analysis, Writing - original draft. V. Sriram: Conceptualization, Formal analysis. I. Didenkulova: Conceptualization, Formal analysis.
Acknowledgements
This paper presents the results obtained through the collaborative research project work between the Indian Institute of Technology, Madras, and Nizhny Novgorod State Technical University n.a. R.E. Alekseev, Russia on ‘Impact of waterborne debris on the nearshore structures during extreme coastal floods’ funded by Department of Science and Technology, Government of India (Grant agreement no. INT/RUS/RFBR/P-203) and RFBR (15-55-45053). VS also acknowledges the support from DST-INSPIRE. ID
References (46)
- et al.
Improved relaxation zone method in SPH-based model for coastal engineering applications
Appl. Ocean Res.
(2018) - et al.
Analysis of long-period cosine-wave dispersion in very shallow water using nonlinear Fourier transform based on KdV equation
Appl. Ocean Res.
(2016) - et al.
Hydrodynamic analysis and optimization of the Titan submarine via the SPH and Finite–Volume methods
Comput. Fluid
(2018) - et al.
A coupled finite volume–smoothed particle hydrodynamics method for incompressible flows
Comput. Methods Appl. Mech. Eng.
(2016) - et al.
Tsunami generation in a large scale experimental facility
Coast Eng.
(2016) - et al.
Tsunami evolution and run-up in a large scale experimental facility
Coast Eng.
(2016) - et al.
Improved MLPG_R method for simulating 2D interaction between violent waves and elastic structures
J. Comput. Phys.
(2012) - et al.
A hybrid method for modelling two dimensional non-breaking and breaking waves
J. Comput. Phys.
(2014) - et al.
SPH modeling of dynamic impact of tsunami bore on bridge piers
Coast Eng.
(2015) - et al.
Steepness and spectrum of nonlinear deformed shallow water wave
Ocean Eng
(2008)
Extreme run-up events on a vertical wall due to nonlinear evolution of incident wave groups
J. Fluid Mech.
SPH modelling of water waves generated by submarine landslides
J. Hydraul. Res.
On the runup of long waves on a plane beach
J. Geophys. Res.: Oceans
Numerical solution of the Navier-Stokes equations
Math. Comput.
3D SPH simulation of large waves mitigation with a dike
J. Hydraul. Res.
Modeling dam break behavior over a wet bed by a SPH technique
J. Waterw. Port, Coast. Ocean Eng.
Smoothed Particle Hydrodynamics for Water Waves
Tsunamis and challenges for accurate modeling
Oceanography
“New Trends in the Analytical Theory of Long Sea Wave runup.” Applied Wave Mathematics: Selected Topics in Solids, Fluids, and Mathematical Methods
Resonance phenomena at the long wave run-up on the coast
Hazards Earth Syst. Sci.
Three-Dimensional SPH modeling of a bar/rip channel system
J. Waterw. Port, Coast. Ocean Eng.
Smoothed particle hydrodynamics - theory and application to non-spherical stars
Mon. Not. Roy. Astron. Soc.
On the state-of-the-art of particle methods for coastal and ocean engineering
Coast Eng. J.
Cited by (11)
A three dimensional hybrid fully nonlinear potential flow and Navier Stokes model for wave structure interactions
2022, Ocean EngineeringCitation Excerpt :viscous solvers based on the Navier-Stokes (NS) equations. Algorithms based on inviscid solvers are fast and energy preserving making them readily applicable to long duration wave simulations in large domains (Manoj Kumar et al., 2020). However, owing to an inviscid and irrotational framework of the governing equations, viscous or aeration effects in the near-field of a structure (particularly at low Keulegan-Carpenter numbers) or in the surf zone cannot be accounted for.
Three-dimensional coupling between Boussinesq (FEM) and Navier–Stokes (particle based) models for wave structure interaction
2022, Ocean EngineeringCitation Excerpt :The 2D MLPG_R model was coupled with a fully nonlinear potential theory (FNPT) model for simulation of breaking waves (Sriram et al., 2014). This hybrid model was used for simulating run-up of long waves in Manoj Kumar et al. (2020) and impact of breaking waves on elastic wall in Kumar and Sriram (2020). Finally, the MLPG_R model in 3D was validated against experimental results for interaction of cylinder with focusing waves in Agarwal et al. (2021, 2019a).
Three-dimensional numerical simulations for mitigation of tsunami wave impact using intermittent sea dikes
2022, Ocean EngineeringCitation Excerpt :Jiang et al. (2017) assumed that a solitary wave generated through a single forward stroke of piston/paddle wave maker may represent a tsunami wave. A typical characteristic observed during the impact of several tsunami waves on the on-shore buildings is an aerated bore front impacting on the structure, followed by the continuous incoming bore for some time (Manoj Kumar et al., 2020; Koyano et al., 2021). In the initial stages of impact, the coastal structures experience an impulsive force due to impact of aerated bore front.
Assessment of one-way coupling methods from a potential to a viscous flow solver based on domain- and functional-decomposition for fixed submerged bodies in nonlinear waves
2022, European Journal of Mechanics, B/FluidsCitation Excerpt :In [30], the FNPF model, solved with a BEM, is strongly coupled with a NS solver based on a Finite Element Method (FEM), and applied to a wave breaking problem. This numerical model is later applied in [31] to the estimation of long wave run-up. A solver denoted “qaleFOAM” was also developed, coupling a FNPF solver using a Quasi Lagrangian Eulerian FEM with OpenFOAM®.
Tsunami-like flow induced forces on the structure: Dependence of the hydrodynamic force coefficients on Froude number and flow channel width in quasi-steady flow phase
2022, Coastal EngineeringCitation Excerpt :One of the typical characteristics of a tsunami inundating the coast is the initial aerated bore front followed by the continuous incoming bore for few minutes to hours (Sundar et al., 2007; Mori et al., 2011; Rossetto et al., 2011; Goseberg et al., 2013; Sriram et al., 2016; Manoj Kumar et al., 2020; Koyano et al., 2021).
Review on the local weak form-based meshless method (MLPG): Developments and Applications in Ocean Engineering
2021, Applied Ocean Research