SPH Numerical Modeling for the Wave–Thin Structure Interaction

In this paper, a numerical model of 2D weakly compressible smoothed particle hydrodynamics (WCSPH) is developed to simulate the interaction between waves and thin structures. A new color domain particle (CDP) technique is proposed to overcome difficulties of applying the ghost particle method to thin structures in dealing with solid boundaries. The new technique can deal with zero-thickness structures. To apply this enforcing technique, the computational fluid domain is divided into sub domains, i.e., boundary domains and internal domains. A color value is assigned to each particle, and contains the information of the domains in which the particle belongs to and the particles can interact with. A particle, nearby a thin boundary, is prevented from interacting with particles, which should not interact with on the other side of the structure. It is possible to model thin structures, or the structures with the thickness negligible with this technique. The proposed WCSPH module is validated for a still water tank, divided by a thin plate at the middle section, with different water levels in the subdomains, and is applied to simulate the interaction between regular waves and a perforated vertical plate. Finally, the computation is carried out for waves and submerged twin-horizontal plate interaction. It is shown that the numerical results agree well with experimental data in terms of the pressure distribution, pressure time series and wave transmission.


Introduction
For wave-structure interaction problems, the interface between fluid and the structure always causes more attention in numerical simulations. Day and Potts (1994) reported that an unsuitable size of elements adjacent to the fluid-structure interface may cause ill-conditioning of the stiffness matrix and high stress gradients in the FEM (Finite Element Method). The IB (Immersed-Boundary) method is imposed to model interfaces in the FDM (Finite Difference Method) (Fadlun et al., 2000), Karimanal and Nair (2000) evaluated the accuracy when a zero-thickness conducting plate was used to model a thick conducting plate in the FVM (Finite Volume Method). Nowadays, the Smoothed Particle Hydrodynamics (SPH), a meshless Lagrangian method, is widely adopted for the simulations of the fluid flow and fluid-structure interaction. The solid boundary treatment is also the key problem in applying the SPH. A few techniques in dealing with the boundaries are available, including the repulsive force boundary condition, dynamic boundary condition, standard ghost particle boundary treat-ment method, and fixed ghost particle boundary treatment method. But when dealing with very thin plate whose thickness is negligible, all of the techniques mentioned above may fail due to the fact that the particle supporting domain on one side of the plate may cover a domain on the other side, leading to the interaction of fluid particles on each side which does not exist physically. Monaghan (1994) introduced a repulsive force boundary condition (RBC), in which a group of virtual boundary particles with the repulsive force to fluid particles near the boundary are applied to prevent fluid particles passing through the boundary. Crespo et al. (2007) developed a dynamic boundary condition (DBC), in which the boundary particles share the same equations of the continuity and state as particles placed inside the domain, whereas their positions and velocities remain unaltered or prescribed externally. Another method to model the solid boundary is the ghost particles approach which exactly enforces the boundary conditions in inviscid free-surface incompressible flows (Colagrossi and Landrini, 2003). In the ghost particle tech-nique, the fluid particles approaching the solid boundary are mirrored in respect to the body profile, in a layer with the size being equal to the adopted kernel radius. This method can be easily applied to straight profiles and right angles, but it becomes more complicated when a generic solid shape is considered. Shao and Lo (2003) adopted ghost particles to overcome the truncation by the solid boundary. Lee et al. (2008) introduced uniform distributed ghost particles to handle the solid boundary. Marrone et al. (2011) proposed an enhanced version of this method, introducing fixed ghost particles. In this case, the ghost particles are fixed in a reference frame of the body and the values attributed to these particles are calculated at their interpolation nodes located inside the fluid domain. The main advantage of fixed ghost particle technique over the standard ghost particle method is that their distribution is always uniform because it does not depend on the fluid particle positions.
For wave-structure interaction problems, the SPH modeling have been addressed to study horizontal decks (Gómez-Gesteira et al., 2005), vertical wall breakwater (Didier et al., 2014), porous submerged breakwaters (Ren et al., 2014), freely floating body (Ren et al., 2015), and oscillating water chamber (Didier et al., 2016). Nearly all structure boundaries in their model are surrounded by fluid in one side, or the structure has a significant thickness. Thus, there is no need to pay much attention to the structure thickness. But for the structures having a very small value of the thickness over the width ratio (thin structure or zero-thickness structure), great care should be taken in dealing with the boundary conditions. Meringolo et al. (2015) mentioned that the numerical simulation of a thin structure leads to a small initial spatial resolution, resulting in an increase of the computational cost. They introduced the multi-node fixed ghost particle method to model solid boundaries in which the thickness of the body is small, and the fluid surrounds the wall from more than one side. Thus, the particle number in the thickness direction is reduced to half that of the traditional method. With this method the computational cost decreased significantly.
The paper is organized as follows. The SPH method is briefly described in Section 2. A particular focus is on a color domain particle (CDP) technique, a new thin structure boundary treatment technique. The results of the model calibrations and case studies are then presented in Section 3. We finish with the conclusions in Section 4.

Governing equations
The governing equations for simulating the viscous flows are the mass and momentum conservation equations. With regard to a fluid particle in the SPH framework, the Lagrangian form of the Navier-Stokes equations are written as follows: where ρ is the fluid particle density; t is time; u is the particle velocity; p is pressure; ν 0 is the kinetic viscosity of laminar flow (10 -6 m 2 /s), g is the gravitational acceleration; f ext is the external body force, τ is the sub-particle scale (SPS) stress tensor, and can be expressed as follows: where is the turbulence eddy viscosity, k is the turbulence kinetic energy, C s is the Smagorinsky constant (0.12), C I =0.0066, Δl is the particle-particle spacing, , and S ij the element of the SPS strain tensor, where i and j denote the coordinate indexes. For the WCSPH model, the state equation is used to determine the fluid pressure: where c 0 is the artificial sound speed, ρ 0 is the reference density, and γ is the polytrophic constant (γ=7 for water). For computational reasons, it is a common practice to assume that c 0 takes the value obtained from the following constraint (e.g., Monaghan, 1994;Colagrossi and Landrini, 2003): where d is the water depth. For water wave problems, the characteristic velocity u in Eq. (6) is set to be c w , which is the wave celerity . Although the physical sound velocity is much larger than the value chosen for c 0 , such a condition does not affect the actual flow dynamics and it ensures that the density oscillations would not exceed 1% (Monaghan, 2005).

SPH discrete formula
The SPH method is based on the interpolation theory, in which the fundamental principle is to estimate any function A(r) by collection of N particles interacting with each other.
where i and j are the particle indexes, r is the position vector, r ij is the relative vector from position i to j, h is the smooth length, and is the kernel function. Following Colagrossi and Landrini (2003) and Cherfils et al. (2012), the discretized SPH equations are described as: where , F ij is the relationship function, which will be discussed later in Section 2.3, while denotes the gradient taken with respect to the coordinates of Particle i, and is the kernel function in an abbreviated form. In current SPH model, a Gaussian kernel, with a cutoff radius set at 3h is used, h/Δx =4/3 and Δx is the initial inter particle distance. The adopted kernel is expressed as follows: where α D is 1/(πh 2 ) in our 2D simulations.
In Eq. (8), the continuity equation, a Rusanov flux has been added, following the work by Ferrari et al. (2009). It helps to increase the scheme stability by reducing the density fluctuations, which are often observed with weakly compressible models. The Rusanov flux involves the effective numerical sound speed, , which can be given by: n ij is the normal vector between particles i and j, and can be given by: 2.3 Boundary conditions and colored domain particle technique Nearly all SPH methods need to build up the particle link list. The link is built up on the particles, with a distance between them smaller than R=σh, and σ is the coefficient determined by the kernel type and simulation problem characteristics. A particle's support domain is Φ(r ij <σh) (see Fig. 1), following Monaghan and Lattanzio (1985).
There is no problem using traditional boundary treatment technique when the particles are far from the boundaries. However, when it comes to a boundary, the support domain is truncated. A boundary treatment is needed. Here, for all default boundaries, following Marrone et al. (2011), the fixed ghost particle method is applied. According to Cher-fils et al. (2012), to prevent the particles from the penetrating boundaries due to the impact or the time step not sufficiently small to keep all the fluid particles inside the domain, an extra force, as described by Eq. (14), is added. In the equation, D is the fluid height in the simulated case, r ib and x ib are the distance and vector from the corresponding boundary point to particle i respectively.
Another difficulty we are facing now is a thin structure inside the computation domain (Fig. 2). Applying traditional fixed ghost particle method could lead to an error link between different sides of the structure. For example, the particles in Ω 0 should not interact with particles in Ω 1 directly. But the default nearest neighbor particle searching procedure will setup a link between Particles i and j, as illustrated in Fig. 2, which is physically impossible. Therefore, we need to distinguish the particles from different regions, and we also need to know how the particles from different regions influencing each other. The solution is a relationship function for the particles from different regions. Firstly, a particle needs to contain a region list which informs other particles whether it can interact with itself (whitelist) or not (blacklist). Here we setup the region list that a particle has relationship with (whitelist). To implement this idea, the computational domain is divided into sub-domains -the internal domains and near boundary domains. The thickness of near boundary domains is usually  REN Xi-feng et al. China Ocean Eng., 2018, Vol. 32, No. 2, P. 157-168 equal to the support domain radius, R. In this way, we can control the relationship of the particles in different domains sufficiently and efficiently. A color value is assigned to each particle. The value contains the information of the domain that the particle belongs to and the domain whitelist that the particle can interact with. Moreover, with a particle color value, we setup a function to check the relationship between the particles, and drop the link if the two particles should not interact with, although they may be very close to. The same rule is also applied to ghost particles.
To minimize the additional computational consumption, the following relationship function is used: F(i, j) is the relationship function; C i and C j are colors of Particles i and j, respectively; & is the bitwise AND operation. The particle's color can be assigned as: C ΩN is the color of domain Ω N , taking the value of I ΩN , I ΩN is the domain ID of domain Ω N , which can be given as required. Here, we use a 32bit integer to save a particle's color. Therefore, a particle, not near the boundaries, which are called the background particles, should interact with all particles. Its color is assigned a hex value #FFFFFFFF. Thus, all the bits of the color are 1. When we do the bitwise operation AND on the colors of the background particle and other particles, none zero value is induced. It should be noticed that the domain ID has a range from 0 to 31, because we have used a 32bit integer representing a particle's color.
is the domain ID of the particle i belongs to. is the domain ID that the particles, in domain i, may interact with. | in Eq. (16) is the bitwise OR operator. If the particles within domain i interact with multi domains, more bitwise OR operation should be taken. Thus, the particle color contains a relationship white list. Fig. 3 shows a more complex thin structure boundary domain setup. With the thin structure considered, the particles in Ω 0 may interact with the particles in Ω 2 , besides the background particles, and they should be assigned with a color value, 5. Particles in Ω 2 may interact with the particles in Ω 0 . They should be assigned with , also. It should be noted that if the domains are far from each other, like Ω 0 and Ω 5 , the domain ID could be reused. Therefore, the domain number is not limited to 32.
Another instance of applying the color particle domain technique is the geometrical singularities of the plate ends. The ends of a thin structure need to be paid more attention because of the connection of two sides boundaries of the plate. Its fluid domain was divided into three sub-domains as illustrated in Fig. 4. The interaction lists of each domain near singular end are given in Table 1. In Table 1, the domain is the background domain, in which the particles can interact with the particles from any domain. N' means the corresponding ghost boundary domain of the near boundary domain Ω N with respect to the structure body profile, and it shares the same domain ID with Domain Ω N . Another thing we should consider is that when a particle from Domain Ω 2 goes into Ω 1 or Ω 3 at a very tiny distance from the singular end, the tiny distance may result in a very large repulsive force deduced by Eq. (14). To prevent a particle from too close to the singular end, the repulsive force by Eq. (14) is also applied between a particle in Ω 2 and the singular end.

Free surface detection
For an absorbing wave maker, the free surface in front of the wave maker needs to be determined dynamically (Gao et al., 2012). There are two main ways to obtain the surface elevation. One way is to integral the particles volumes within a snap range along x axis (Eq. (17)).   where x lb and x ub are the selected integral gap lower and upper bound, respectively. Δx is the initial particle space. The other way is interpolating the surface particle positions. To estimate the water surface elevation at Position x 0 in a water flume, the nearby surface particles are detected first by the fast free-surface detection method (Marrone et al., 2010), within a range of the distance to x 0 not larger than σh. Then the spline interpolation is used to obtain the value at x 0 (Fig. 5). A comparison of the surface elevation time series by the interpolation of the surface particles and by Eq. (17) is shown in Fig. 6. It is obvious that the interpolated time series is smoother and does not shift from the still water level, compared with the time series deduced by Eq. (17).

Hydrostatic test case
As described by Meringolo et al. (2015), without any additional treatment, the fixed ghost particle method (Marrone et al., 2011) and even though enforced by the multinode approach technique, cannot model a zero-thickness plate. Because when the ghost particle boundary method is used, if a thin structure is immersed in a fluid mass, the initial particle space is driven by the thickness of the structure. For the original fixed ghost particle method, the initial particle space Δx≤h/8, and for the multi-node enhanced fixed ghost particle method, Δx≤h/4. The number of total particles presents a power function type growth, the expo-nent is 2 for a 2D case, and 3 for a 3D case. When h is close to zero, it approaches infinity. The ghost particle method also in need of the initial particle space is set to Δx≤h/8, without the color domain particle enhancement.
The hydrostatic test involves a tank, divided into two subdomains by a zero-thickness plate at the middle of the tank, containing still water at different levels in the subdomains. In particular, the length of the tank is L=1.0 m, the plate is placed at x=0.5 m, and the water level in the left subdomain is d 1 =0.4 m, while in the right subdomain it is 0.1 m, see Fig. 7. The adopted resolution is Δx=0.01 m.
Firstly, the thin plate is modeled with the standard fixed ghost particle method. Following the basic process of the fixed ghost particle method, the ghost particles are generated for all boundaries and the thin plate. As foreseeable, the generated ghost particles may lie in the other subdomain or override other fixed ghost particles (Fig. 8). Fig. 9 shows the unphysical pressure shock with a consequent appearance of a relevant velocity field, which is due to the wrong relationships between the fluid particles in the left tank and right tank, the wrong relationship between the structure bound ghost particles and fluid particles. When the color domain particle enforcement is employed, the correct solution is conducted. In this case, the colored fluid domains and boundary domains are illustrated in Fig. 10 respectively. To be succinct, the default fluid domain with ID=0 is omitted. Although there is a simpler way    to deal with this problem (We just need to color the left tank and boundaries generated for the left tank with one value, color the right tank and boundaries generated for the right tank with another value, and we set them not to interact with each other, the simulation will go as expected), in order to present the general process of applying the color domain method, we still use the generalized color domain process that may be somewhat more complicated. Firstly, the default fluid domain, ID=1, is assigned along the water tank, as well as the default boundary domain, also ID=1, is assigned synchronously. As ordinary domains, the default fluid domain and boundary domain can interact with any other domain. When the middle thin plate is placed in, the left and right side fluid domains are colored with different values to distinguish them from each other, also the boundary domain is colored correspondingly. When a structure is connected to the bottom edge, both the fluid domain and boundary domain require some extra trick. First, the boundary domains will overlap due to the presence of the corners (Fig.  10), thus, we color them with different values and setup the relationship list. Then, come to the fluid domains, since we assume that the default fluid domains and boundary domain can interact with any domain basically, to prevent the fluid particle in one side from interacting with the boundary particles generated for another domain, the fluid domains, close to the bottom, are increased to two times the width of the support domain radius, see Fig. 11. The coordinating domain interaction whitelist is given in Table 2.
It also should be noted that the division of the domain is not unique, and it just needs to be coordinated with the domain interaction whitelists.

Hydrodynamic test case
To estimate the effect of the current boundary treatment method on hydrodynamic problems, a more complex structure, a fully perforated breakwater with zero-thickness was studied here.
The schematic representation of the full perforated breakwater is illustrated by Fig. 13, where L=4 m, d=0.4 m and the chamber width, L b =0.54 m. For a 2D problem, the

Region ID
Interaction region ID whitelist 0 0, 1, 2, 3 1 0, 1, 2, 3 2 0, 1, 2 3 0, 1, 3 ε = d 2 / (d 1 + d 2 ) = 0.25 porosity of the breakwater is , where d 2 is the height of the slot and d 1 is the height of the solid part. To compare with literature results (Meringolo et al., 2015), the height and period of the incident regular wave are respectively H=0.08 m and T=1.03 s, resulting in the second-order Stokes wave trains under the intermediate water depth conditions.
To model the zero-thickness perforated wall of the caisson, the colored fluid particle method and single point treatment are applied. The singular point treatment described in Section 2.3 is also used here. The colored fluid domains and boundary domains are sketched in Fig. 14 respectively.
Although we could use any particle size by involving our colored-particle method, the particle size is also important for detailed velocity field and pressure field. By considering the computational domain size and the size of the hole on a perforated wall, the particle size, we used here, is Δx=0.005 m. Thus, the total number of the fluid particles is 64000. On an Intel E5-2690@3.0 GHz with 64 GB RAM work station, it takes about 8 hours for a simulation of 10 s.
In Fig. 15-Fig. 17, three significant time instants of a wave interacted with the fully perforated breakwater from Meringolo et al. (2015) (powered by the multi-node boundary enhancement), current SPH model and VOF mode are displayed. The first two instants correspond to the maximum dynamic pressure occurring at the SWL for the front and rear faces, respectively. The last instant refers to the occurrence of the wave trough at the perforated front wall. These three instants are chosen to verify the stability of different approaches to the boundaries. The results show that no spurious flow process appears near the slotted front wall at any instants. From the phenomenon of different models, the SPH models show superiority to the VOF model when dealing with the free surface motion. The pressure field from different models show good agreement with each other. In the horizontal velocity field, the current SPH model shows better agreement with the VOF model than Meringolo et al. (2015) SPH model in the distribution trend. But the SPH model used by Meringolo et al. (2015) has a more continuous distribution than the current model. The good agreement of the current model with the VOF results and literature results indicates that the treatment of singular point is efficient. Fig. 18 gives the results in terms of the spatial distribution of the dynamic pressure along the reference walls. It is very clear that the current SPH model has a better agree-   REN Xi-feng et al. China Ocean Eng., 2018, Vol. 32, No. 2, P. 157-168 ment with the empirical formula of Tabet-Aoul and Lambert (2003) than the Molteni and Colagrossi (2009) model used by Meringolo et al. (2015). The disagreement of the pressures near the slots between numerical result and the empirical formula is that the empirical formula assumes a linear variation between the control points and does not consider the effects induced by the presence of the slots on the wall. The values of the Mean Square Error Percent (MSEP) of different models at the front and rear walls of the perforated breakwater are shown in Fig. 19. The results show that the current SPH model has a lower MSEP value on the both walls compared with the Molteni and Colagrossi model (δ=0.2) used by Meringolo et al. (2015). On the rear wall, the value is between the values of the models of Molteni and Colagrossi (δ=0.1) and Molteni and Colagrossi (δ=0.2).

Submerged twin-horizontal plate breakwater
The focus of the current SPH model is to handle the structures surrounded by fluid. Thus, here the submerged twin-horizontal plate breakwater model is applied. The experimental results from Li et al. (2010) are used for the  comparison. The numerical computation domain is mainly the same as the experiment setup, except the distance from the wavemaker to the upwind side of the plates and the distance from the lee side of the plates to the end of the wave flume, see Fig. 20. The plate width, B, is 1 m, the space between the two plates, S, is 0.1 m, the water depth, d, is 0.48 m and the submerged depth, d', is 0.08 m. The pressure transducers are located on the upper and lower sides of each plate, detail information is given in Fig. 21.
For a long term simulation, the flume must absorb the reflected wave and transmitted wave. Thus, an active absorption wave maker is needed. According to Gao et al. (2012), for a piston-type wave maker, the horizontal velocity of the wave maker is written as: where x 0 is the wave maker position, U is the horizontal velocity, η 0 and η are the target wave surface elevation and local surface elevation, respectively, and ω is the wave circular frequency. Φ is the transfer function, according to Hughes (1993), using the first order theory, it is given as: where κ is the wave number. The wave damping area is enforced by applying a damping body force to the particle in the damping area. The external damping force f ext , in Eq. (9), can be given as: where α is a shape coefficient, equal to 3.5; β is the damping speed coefficient, which can be determined by numerical tests, here it is 50; x 0 is the damping zone starting point; x is the particle position; L d is the damping zone length; u i is the particle velocity.
In the numerical computations, the initial spatial resolu-   tion is derived by the space between the two plates. To obtain reasonable results, sufficient particles are needed, Thus, the initial particle space Δx=0.01 m is used. The distance from the wavemaker to the upwind side of the plates is 2L, and the distance from the lee side of the plates to the start of the damping zone is 2L. The damping area length equals to L, and L is the wave length. Therefore, this results in a total of 90000 fluid particles. For a typical computation situation, with L=2.5 m, a 40 s simulation takes about 7 hours on a work station, power by Intel ® Xeon E5 2690@3.0GHz and 64GB RAM. Fig. 22 shows the surface elevation time series of the numerical and experimental results for the case characterized by B/ L=0.40, H/d=0.17, S/d=0.21 and d'/d=0.17, the gauge is placed at the distance of 1.58L from the lee side of the plates. The current SPH mode results show good agreement with the experimental results, especially the frequency of the wave peaks and the amplitude of the wave trough.
The reflected and transmitted waves are very important for a breakwater. The incident and reflected wave spectra, deduced from the application of the split method by Goda and Suzuki (1977), are described in Fig. 23. As noticeable, the incident wave spectrum has only one wave harmonic, the fundamental peak frequency, as the destination wave is a linear wave, according to Le Méhauté (1976). Whereas, the spectral density of the reflected wave spectrum shows an additional secondary wave harmonic at two times the fundamental frequency besides the base one. Fig. 24 shows the amplitude spectra of the transmitted wave by the current SPH model besides the experimental results of Li et al. (2010). The spectral energy is concentrated in the fundamental incident wave peak frequency, and two and three times the frequency. And the energy at the secondary wave frequency is the largest. Fig. 25 shows the comparisons on the pressure series of the numerical and experimental results (Li et al., 2010) of the pressure transducers, which are placed on each side of the two plates, and the marker and position are detailed in Fig. 21. In general, the numerical results show good agreement with the experimental results, either the pressure oscillation amplitude or the frequency. But the pressure series of #1 and #5 pressure transducers, which are located on the top of the top plate, and are the nearest to the upwind side and nearest to the lee side, respectively, is not as much conform as others. This is likely due to the complexity of the wave break on the top of the plate.
The transmission coefficients, K t , of the numerical simu-    lations and experiments by Li (2014) are given in Fig. 26. An overestimation of K t obtained by the SPH method can be observed for most B/L values. The minimum value of the transmission coefficient occurs near B/L=0.4, and the value is 0.327, which is very close to the experimental value, 0.3.

Conclusions
A new boundary treatment technique, Color Domain Particle (CDP), is proposed and implemented to deal with zero-thickness structures surrounded by fluid. With the help of CDP, the particles nearby a zero-thickness structure are prevented from interacting with the particles on the other side. To validate the improved SPH model, a static water tank with a thin plate at the middle resulting in two subdomains with different still water levels was examined. The results show that the pressure and the velocity field are exactly the same as the analytical solution, and no mixing occurs. For further evaluations, a more complex perforated breakwater problem is studied. The pressure field and horizontal velocity field show good agreement with literature results. The current SPH model delivers more accurate pressures on the front wall of a fully perforated breakwater. The results indicate that the CDP boundary treatment technique works well with thin structures immersed in water. At last, a study on a submerged twin-horizontal plate breakwater is  carried out. The transmitted wave time series agree with the experimental results very well. The pressure values at either side of the plates also show good agreement with the experiment values. The transmission coefficients of different wave lengths are also compared with the experimental results. The minimum transmission coefficient value occurs near B/L=0.4 and is very close to the experimental value. It should also be noted that the new boundary treatment method can decrease the computation cost remarkably, because the thickness of the thin structure is not to be considered when using the new boundary treatment technique.