Wave Interaction and Overwash with a Flexible Plate by Smoothed Particle Hydrodynamics

The motion of a flexible elastic plate under wave action is simulated, and the well–known phenomena of overwash is investigated. The fluid motion is modelled by smoothed particle hydrodynamics, a mesh-free solution method which, while computationally demanding, is flexible and able to simulate complex fluid flows. The freely floating plate is modelled using linear thin plate elasticity plus the nonlinear rigid body motions. This assumption limits the elastic plate motion to be small but is valid for many cases both in geophysics and in the laboratory. The principal conclusion is that the inclusion of flexural motion causes significantly less overwash than that which occurs for a rigid plate.


Introduction
The interaction of a flexible floating plate with wave forcing has been the subject of extensive research due to the application in offshore engineering, polar engineering, and geophysics [1][2][3][4]. The general assumption is that the wave amplitude is small enough that the problem may be treated as linear. The first solution for a linear plate in two-dimensions was by [5,6] and in three-dimensions for a circle [7,8] and for an arbitrary plate [9]. We refer to the review article [10] for an extensive description of the literature. The primary focus of the research has been the linear response to small amplitude waves, especially to model very large floating structures, such as a floating airport [11][12][13]. There has also been extensive work on the modelling of wave propagation in the frozen ocean [14][15][16].
Recently there has been a strong interest in understanding the nonlinear motion of a floating elastic plate. While the wave amplitude may be assumed small enough that linear wave theory applies [17], the floating elastic plate is generally employed to model bodies which have a minimal freeboard. For example, the floating runway built in Japan has a freeboard of only around 1 m. Similarly, ice floes have a freeboard of approximately 0.1m. This small freeboard means that it is effortless for fluid to overtop the structure, a phenomenon referred to as overwash. This overwash has been the subject of extensive recent research both theoretically and experimentally [18]. In particular, experiments have shown that, even if the linear theory describes well the motion of the plate, there is significant overwash [19,20]. This overwash has been shown to have a substantial effect on the reflection and transmission [21][22][23] Moreover, overwash will considerably enhance the melting of ice. The phenomenon of overwash is not limited to ice floes and appears in other areas of ocean engineering [24].
There is a peculiar paradox when speaking of overwash. It is a feature of experiments conducted in the laboratory [19,25], and often special barriers are required to avoid the phenomenon [25]. However, it has yet to be observed in the field, and there is some speculation that overwash only exists in laboratory experiments. If this were the case, it would appear to contradict the scaling laws of fluid motion. It seems likely that the phenomenon of overwash does exist in nature but that it leads to rapid ice melting, and so is transitory and difficult to observe. We emphasise here that we do not regard overwash as an event restricted to the laboratory setting.
The theoretical work on investigating overwash has focused on the use of either bespoke numerical methods [20] or OpenFOAM [26][27][28]. However, only in the work of [27] was the elastic motion of the floating plate included. This omission of flexure is even though the phenomenon of overwash generally only appears when the body is sufficiently thin that the elasticity is important.
Smoothed particle hydrodynamics (SPH), a Lagrangian meshfree particle method, is now used widely for computer simulations of free-surface and multiphase flows [29][30][31]. The evolution of any interfaces is easily tracked just by the motions of the SPH particles. This characteristic is an advantage of SPH over the conventional mesh methods in which additional tracking algorithms are needed. Additionally, since the material derivative is used, convection terms are implicitly calculated, and hence numerical diffusion is avoided. One of the earliest work on the application of SPH to simulate water wave breaking on a beach, overwashing a flat plate, and from a dam break is due to Dalrymple and Rogers [32]. Although it is just a phenomenological study, this study highlighted the potential of SPH for challenging multiphase problems. Crespo et al. [33] then simulated the interaction between sea waves and an oscillating water column device. Wen et al. [34] simulated wave interaction with a vertically cylindrical breakwater in a water basin with SPH. In those studies, the structures are undeformed and move as rigid bodies. For deformable structures, the problem becomes even more challenging as a separate solver for the structure deformation is needed. Recently, the extension of the SPH method to floating elastic bodies was considered by Zhang et al. [35]. In that study, Zhang et al. investigated the response of an elastic plate in mild wave conditions with wave steepness ranging from ka = 0.1 to 0.15, in which k and a are the wavenumber and wave amplitude, respectively, and the ratio of the plate length to the wavelength is 1.0. Together with the wave steepness, this ratio is one of the most critical factors influencing dynamical behaviour of floating structures [36]. For example, iceberg velocity was found to be close to fluid particle velocity if the ratio is more than 3.0. Likewise, heave amplitude is similar to the wave amplitude if this ratio about 0.1 [37] and heave resonance is also observed if it is larger 12.0 [38].
In this study, water wave interacting with an elastic floating plate in a water basin is investigated. In contrast to [35], wave conditions in this study are more severe with wave steepness up to 0.5, plate stiffness varies in a wide range of values, from 25 MPa to 2.5 GPa, and the ratio of the plate length to the wavelength is less than 11.0. The fluid phase is solved with an incompressible SPH. In contrast, the elastic behaviour of the floating plate is solved using the thin plate theory with the hydrodynamic force input from the fluid solver. In turn, the deformed plate has a forcing effect on the fluid phase. Water waves are generated using a moving paddle at prescribed frequencies and displacement amplitude. Our study has application to the storm damage of very large floating structures and the melting and breakup of sea ice. However, our focus here is on simulating at the laboratory scale where extensive experiments have been conducted. The rest of the paper is organised as follows: wave generation in a water basin is discussed in Section 2, the SPH method and the two-phase SPH model are presented in Section 3, numerical results are discussed in Section 4, and a brief discussion is given in Section 5.

Wave Generation
Water waves are generated in a water basin, as illustrated in Figure 1. The left boundary is moved to generate water waves. Its displacement follows sinusoidal rule, Its velocity is then u mw = Sω sin (ωt) , (2) in which ω = 2π/T mw (where T mw is the period of the moving boundary) and S are frequency and displacement amplitude of the left boundary motion, respectively. From the wave-generation theory [39], the generated water waves have the same period of the moving boundary, i.e., T = T mw , while its celerity is evaluated from in which d is the initial water height. Since λ = cT, wave length of the generated wave can be obtained by numerically solving the non-linear equation Wave amplitude is estimated from The right boundary is a slope, behaving like a beach, to absorb the wave energy and therefore suppress reflecting waves. Beach configuration for wave absorption is commonly used in water-basin experiments due to its simplicity. In simulations, some other boundary conditions, such as opened boundary [40] or sponge-like boundary [41], have also been utilised.

SPH
SPH approximates a function by its convolution with a smoothing function having compact support, The smoothing function, W(r, h) is isotropic and Upon domain partition into equal sub-domains, then referred as SPH particles, the integral in Equation (6) becomes sum of weighted values of SPH particles within the support domain Ω, Derivative of a function is approximated in the same manner, Taking integration by parts and then discretizing, Equation (9) becomes It is noted that ∇ i W ij = ∂W(|X i − X j |, h)/∂X i is gradient of the weight function with respect to X i . The weight function plays a crucial role in precision of the SPH approximations. In the subsequent simulations, the harmonic-like weight function is chosen since it produces the smallest numerical errors for the SPH approximations [42], where η = r ij /h is the dimensionless distance between i and j, α d for two-dimensional support domain is 0.7103, and h = 1.2 x is chosen.

Fluid Phase
The fluid phase is modelled by SPH particles whose motion is governed by Navier-Stokes equations, in which u i , P i , ν, and g are particle velocity, pressure, viscosity, and gravity acceleration, respectively. Equation (12) is coupled with mass conservation equation, Here, an incompressible SPH (ISPH) scheme is utilized. Prediction-correction algorithm is adopted here for solving Equations (12) and (13) [43][44][45]. Accordingly, intermediate value for SPH particle velocity and position, in the prediction step, is estimated first in which u n i and r n i are velocity and position of SPH particles of the previous step. The viscous term in the right hand side of Equation (14) is approximated by in which ν i and ν j are viscosity of particle i and j, and V j is the particle volume. The parameter = 0.001 h is used to prevent numerical singularity if two particles come too close to each other, that would occur sometimes, especially at the free surface.
In the prediction step that the incompressibility has not been gained, the intermediate fluid density is estimated using Equation (13), which is discretized by The pressure field is attained by solving the Poisson's equation using the intermediate values which can be discretized as follows On applying Equation (19) for all SPH particles, a system of linear algebraic equations is established. Here, we adopt an iteration scheme [44] to solve Equation (19), as such with The iteration starts with the pressure field of the previous step. For small time step, the scheme converges quickly after a few iterations. The velocity field and particle location are then updated The pressure gradient is approximated by

Solid Phase
The floating plate moves freely with the fluid phase except in the horizontal (x) direction. Hence, the plate motions, in two dimensions, include heave (in the vertical direction), pitch (rotation about its centre of mass), and elastic deformation. The heavy and pitch motions are determined froṁ in which F(x, t) is hydraulic force acting on the plate; n z and n p are unit vector of the z-axis and the unit normal vector of the plate, respectively; M p and I 1 are the plate mass and moment of inertia for rotation. The plate deformation is determined using the thin plate theory. Accordingly, equation of motion for a thin plate is where E, I b , A, ρ p are Young modulus, bending moment and cross area of the plate, and plate density, respectively. It is noted that the bending moment of the plate is defined by I b = 4W H 3 /3, in which W and H are halves of the plate width and thickness, respectively, as illustrated in Figure A2. For a free-moving plate, the following boundary conditions are applied Variable-separated solution for the plate deformation under time-dependent acting force is then written where ζ n are the modes of vibration of the plate, as given in Appendix A. On substitution into Equation (27), we obtain Then, we multiple by ζ m and integrate from −L to L and obtain a set of second-order ODEs for a n (t), due to the orthogonality of the modes. Equation (31) is then solved with the initial conditions for the plate at rest a n (0) = 0, a n (0) = 0, to determine the coefficient a n (t) of each mode. In simulation, at every time step, the hydraulic force F is known and therefore a n (t) can be obtained numerically. The plate is modelled by a set of SPH particles, whose velocity is in which d i is distance vector from particle i to the centre of mass of the plate. The elastic velocity term is obtained by taking derivative of ζ with respect to time, Then, position of the plate particles in the next time step is updated by SPH particles also model the fixed and moving walls. The velocity of the SPH particle representing the fixed walls is set to be zero, while Equations (1) and (2) prescribe position and velocity for those of the moving wall. The pressure of SPH particles representing the solid phase is obtained by solving Poisson's Equation (18) [46]. Accordingly, wall pressure arises if fluid particles penetrate the solid phase, leading changes in local density.

Free Surface
In the ISPH scheme, Dirichlet boundary condition P = 0 at the free surface, pressure of SPH particles at the free surface is set to be zero [44]. There are a few algorithms used to identify free-surface SPH particles, such as the weighted density less than 0.9 of the fluid density [47], or divergence of the particle coordinates less than 1.5 [48], i.e., Here, we adopt the latter in the simulations. In order to reduce particle clumping as well as stabilize the free surface due to lack of surface tension, an extra viscosity is applied for the free-surface particles [47] in which u max is maximum flow velocity.

Results
The water basin has dimensions of L bottom = 5 m, L beach = 4 m, and H beach = 1.6 m. Accordingly, the inclined angle of the beach is α = 22 o . The water depth is d = 1.0 m. Displacement amplitude of the left moving boundary is S = 0.1 m and its moving period is varied from T 1 = 1.2 s to T 2 = 1.5 s to produce intermediate-water waves. The floating plate has dimension of 2 cm thickness (2H) and 1 m length (2L). The plate density is ρ p = 500 kg/m 3 , which is half of the water density ρ = 1000 kg/m 3 , and hence the draft is half of the plate thickness, i.e., 1 cm. Spatial resolution is set to be x = 5 mm and the corresponding number of SPH particles is about 259, 340. Time step is set to be 10 −5 s. The fluid viscosity is of 1 Pa s. The simulation length is 7 wave periods and it takes approximately 10 hours for each simulation, which is run parallelly using 40 cores of Intel(R) Xeon(R) Gold 6136 CPU @ 3.00GHz. Figure 2 shows the water waves with period of T 1 = 1.2 s at the observing time of t = 3.5T 1 = 4.2 s. We restrict the domain to x < 8 because the fluid run-up does not extend beyond this point. The longitudinal velocity component is positive in the proximity of the crests and negative in the proximity of the troughs. On the other hand, the vertical velocity component is negative at the upside of the waves and positive at the downside. This is expected to be observed in the propagation process of water waves. The wave's amplitude gradually diminishes as the waves enter the beach region. The pressure field continuously increases, from zero at the free surface, towards the tank bottom. For ease of visualisation and comparison, the pressure field is scaled by the hydrostatic pressure of the water column, P hyd = ρgd = 9.8 KPa. The presence of the waves locally alters the water height, that leads to pressure larger than 1 under the lifted regions and vice versa in the trough region. The beach works pretty well in absorbing the energy of the waves and does not create any significant disturbances in the pressure field. The simulated wavelength and amplitude are λ = 2.35 m and A = 0.187 m, respectively. On adopting Equations (4) and (5), theoretical estimation for those are λ theory = 2.23 m and A theory = 0.191 m. In comparison to the theoretical wave magnitude, the simulated one is slightly smaller due to the viscous effect causing energy dissipation, which is ignored in the theoretical estimation, as well as reflections at the boundaries. The simulation is repeated for the wave period T 2 = 1.5 s. The velocities and pressure are plotted in Figure 3. The simulated wave magnitude and wavelength are 0.159 m and 3.38 m, respectively, while they are 0.162 m and 3.35 m from the theoretical estimation. Again, the simulated wave amplitude is smaller than the theoretical one due to the viscous effect. The good agreement also shows that the particle size, which has specific effects on the generated waves [49,50], is fine enough to achieve the numerical convergence. The pressure field is also linearly increasing from the free surface towards the bottom, and local pressure is more considerable for the larger water depth. The waves are absorbed pretty well by the beach without any significant pressure reflections. It is seen that the wavelength is longer, and the wave amplitude is smaller as the wave period is increased to 1.5 s. In the two cases, the ratio of the water depth to the wavelength is 0.43 and 0.3, respectively, and therefore the generated waves are classified as the intermediate water wave. We note that the SPH method is highly computationally demanding and we, therefore, restrict ourselves to the smallest domain possible in which we have relatively uniform waves with little reflection.

Dynamics of the Floating Plate
The centre of mass of the floating plate initially is placed at X 0 = 4.0 m. The plate density is half the water density, so both the draft and freeboard of the plate equal 1 cm. It is recalled that since the longitudinal motion, or surge of the plate, is restricted, the plate's motions include heave, pitch, and elastic deformation. For ease of analysis, the scaled time t * = (t − t 0 )/T, in which T and t 0 = T(X 0 − L)/λ are the wave period and the time for the first generated wave approaching the floating plate, is introduced. Two values of Young modulus, E 1 = 0.025 GPa and E 2 = 2.5 GPa, are considered to study the effect of elastic deformation on the overwash phenomenon. It should be noted that we are simulating at the laboratory scale. To compare with the geophysical scale, careful dimensional analysis is required. It is noted that conditions for the numerical experiments carried out in the current work are different from the ones investigated by Zhang et al. [35], as summarised in Table 1. More specifically, for observing the overwash phenomenon, wave conditions here are more severe, and Young modulus of the plate is also varied in a broader range of values to examine how elastic deformation affects the overwash. Vibration modes are also analysed.
Firstly, considering incident waves with a wavelength of λ = 3.38 m and wave amplitude of A = 0.159 m. Correspondingly, wave steepness is ka = 0.295, and the ratio of the wavelength to the plate length is 3.38. Figure 4 displays the dynamic response of the floating plate with E = 25 MPa under successive attacks of the incident wave. It is seen that, in the early times, when the wave is not well developed yet, the wave mildly passes underneath the plate. The overwash starts occurring subsequently and gets stronger and stronger subsequently. The overwash is observed at both the ends of the plate, but it occurs stronger at the front end. The plate is significantly bent under the hydrodynamic load from the wave. The simulation is then repeated for a stiffer plate with Young's modulus of E = 2.5 GPa, which is two orders of magnitude higher than the previous case. It is observed in Figure 5 that the overwash phenomenon occurs in a reasonably similar way, but it is more severe at both ends. The plate bending is not obviously observed. The results reveal that the elastic bending of the plate has a specific role in reducing the overwash since it partly absorbs the kinetic energy of the incident waves. Table 1. Summary of the parameters used in the current study in comparison with those investigated in Zhang et al. [35].  The plate dynamically moves with the wave. As observed in Figure 6, heave motion of the plate is also periodic, and its period is as same as that of the incident wave. The heave motion of the plate for both the values of Young modulus is just slightly different, but its maximum and minimum values are quite similar. Within the observing period, those values are 0.14 m and −0.1 m, respectively. In other words, the amplitude of the heave is 0.12 m, which is approximately 0.75 of the wave amplitude. This shows that the heave motion of the plate is not affected much by the elastic deformation of the plate. The pitch motion of the plate is also recorded and shown in Figure 6. Similar to the heave, the pitch is also periodic with the same period of the incident wave, but the phase is slightly shifted. In other words, the pitch motion comes after the heave. In the two cases, the pitch angle varies between −15 • and 16 • .   Figure 7. It is seen that the second vibration mode is dominant regardless of the plate stiffness, although its magnitude is reduced for the higher Young modulus. The reduction in the magnitude of the second vibration mode is almost at the same reduction in the order of magnitude of the Young modulus. Interestingly, the second mode is also periodic with the incident wave. The ratio of wavelength to plate length, (in this case 3.38), is an essential factor for dynamics of the plate [36,38,51], and now it also plays a crucial role in the plate vibrations. Now, the wave with a shorter wavelength, λ = 2.35 m, and a larger wave amplitude, A = 0.187 m, is investigated. Accordingly, the wave condition now is more severe than the previous one because the wave steepness increases to ka = 0.5. The plate dimensions are unchanged. Figure 8 displays dynamics of the plate with E = 25 MPa under this wave condition. It is seen that the overwash occurs much stronger at the front end of the plate. Most of the plate surface is covered by water at t * = 5 due to the overwash. The height of the water layer on the plate surface is also significantly larger. The situation becomes even worse for the stiffer plate with E = 2.5 GPa. The plate completely submerged in the water at t * = 5, as observed in Figure 9. In comparison to the overwash phenomenon observed in the numerical experiments carried out by Zhang et al. [35], it occurs much stronger in the current study.

Parameters
The heave motion of the plate is still periodic with the incident waves regardless of the plate stiffness, as seen in Figure 10. Within the observing period, the averaged heave amplitude is about 8 cm, which is slightly smaller than that in the previous wave condition. This could be caused by the water layer on the top of the plate damping its heave motion. Also, due to the strong overwash, the front end of the plate tends to be pushed downward, which is seen in positive pitch angles in most of the observing period. Deviation in the pitch of the plate with different Young moduli is also seen since t * = 4, that could be caused by varying levels of overwash.
The second vibration mode is seen as still periodic and predominant in this wave condition ( Figure 11). Unlike the previous wave condition where the coefficient a 2 is negative at most of the time, i.e., downward bending, the plate now is bent either upward or downward more equally and forcefully, as seen in sign and magnitude changes of the coefficient a 2 . Another critical point is that the first vibration mode is more evident in this case. This could come from the smaller ratio of the wavelength to the plate length, which is 2.38. Physically, high vibration modes start emerging as the wave frequency gets higher.

Conclusions
The dynamics of flexible floating structures interacting with ocean waves has been a subject of extensive research due to its application in offshore engineering, polar engineering, and geophysic. More specifically, the overwash of a floating elastic plate is of interest since it is used for modelling structures with small freeboard, for instance floating runways or ice floes. Here, the problem of overwash is investigated using SPH, a Lagrangian solver ideally for simulating free surface and interfacial flows. Both the fluid and the floating plate are represented by SPH particles. While the incompressible Navier-Stokes equations dictate motions of the SPH particles representing the fluid phase, SPH particles modelling the plate obey the governing equation from linear plate theory. The SPH simulations are carried out for a laboratory-scale of the problem, i.e., in a water basin. The water waves are generated using a paddle, or a wall, moving at the different periods, T = 1.2 s and 1.5 s, and with a fixed amplitude of = 0.1 m, to result in water waves with different amplitudes, wavelengths, and wave steepnesses of ka = 0.5 and 0.295, respectively. This steepness corresponds to moderately severe wave conditions in reality. Two values of Young modulus, E = 25 MPa and 2.5 GPa, are adopted to investigate how the elastic behaviour of the plate affects the overwash. The numerical results show that as the steepness increases, the overwash occurs more strongly, and the more flexible the plate suffers less overtopping. The plate with E = 2.5 GPa was almost submerged under water just after six wave periods for ka = 0.5. This indicates catastrophic destruction for a floating runway or an ice floe. The study revealed that the deformation of the floating structure partly absorbs the kinetic energy of the fluid phase. Consequently, less overwash is observed for a more flexible plate. The finding agrees very well with observations of Huang et al. [52]. Future research is to build a model in which we simulate ice floe breakup and melting.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Vibrating Modes of a Thin Plate
The plate motion is expanded in terms of the free modes of vibration, which satisfy the homogeneous equation  Table A1 lists value of λ n L of the first six vibration modes of the plate. The modes are given by ζ 2n = N 2n cos(λ 2n x) cos(λ 2n L) + cosh(λ 2n x) cosh(λ 2n L) , (A9) ζ 2n+1 = N 2n+1 sin(λ 2n+1 x) sin(λ 2n+1 L) + sinh(λ 2n+1 x) sinh(λ 2n+1 L) , where N 2n and N 2n+1 are determined from the orthonormality of the vibration modes, L −L ζ n (x)ζ m (x)dx = δ mn . (A11) The first six vibration modes are illustrated in Figure A1. Table A1. The first six values of λ n L.