Hydro-elastic analysis of standing submerged structures under seismic excitations with SPH-FEM approach

https://doi.org/10.1590/1679-78256266 Abstract In this paper, a fully coupled fluid structure interaction (FSI) method is used to investigate the hydro-elastic response of a fully submerged standing structure under seismic excitations. Two different domains (solid and fluid) are modelled by mesh based and meshless methods, respectively. Solid domain is modeled by finite element method (FEM) and fluid domain is modeled by smoothed particle hydrodynamics (SPH). Coupling between FEM and SPH is implemented via contact mechanics, and this method differs from others in the way of coupling mechanism. Invading SPH particles are solved together with finite elements by using contact mechanics, then, a fully coupled method is achieved. In the scope of this research, different seismic excitations are applied to a rectangular tank. Half of the tank is filled with water and a submerged rubber plate is attached to its mid bottom. Thus, two-dimensional motion of rubber plate and water is investigated experimentally and simulated fault Graphical


INTRODUCTION
Many natural effects are considered for design of structures: Base excitations, dead and live loads, wind loads and even blast loads are critical for structures. Most of them are basically formulated for design of structures by engineers. In design phase, it is very important to determine the loads applied to structure and it depends to the nature. Sometimes it is simple to predict and sometimes it is not as in earthquake. The most ambiguous phenomenon is the earthquake in structural design and it will be more complex phenomenon for submerged structures with additional fluid effect.
The additional effect coming from fluid domain is simulated as additional masses applied to the structures submerged in the fluid. According to this theory, the effect from fluid is a force assumed to be equal to the hypothetical mass of fluid in motion times acceleration of structure (Han and Xu 1996). This theory is used by many researchers as coupling method of fluid and structure (Liang et al. 2001). Besides, theories (Housner 1963) for simulation of sloshing effect was proposed in early researches. Basically, these problems can be classified into two as structure in the fluid and fluid in the structure (Jing, Feng, and Cheng 2019); both theories use the same approach and these formulations are very useful for engineers.
Nowadays, researchers propose new methods or combined methods Dinçer et al. 2019;Fernández-Méndez, Bonet, and Huerta 2005) already in use for better simulations. Some researchers use mesh-based methods (Cho and Lee 2004;Ghalandari et al. 2019) to take the advantage of robustness. In contrast, some use particle-based methods to take the advantage of suppleness. Popular ones are Smoothed Particle Hydrodynamics (SPH) (Dinçer, Bozkuş, and Tijsseling 2018;Shao et al. 2012), Moving Particle Semi-Implicit (MPS) Zhang, Wan, and Hino 2014) and Discrete Finite Element Method (DFEM) (Chen and Wang 2019). There are many method combinations (Fourey et al. 2014;Hu et al. 2014;Jiang et al. 2018;Liu and Liu 2019;Salih et al. 2019;Yildizdag et al. 2019;Zhang, Qiang, and Gao 2011) in the literature. All has advantages and disadvantages when compared. Most popular one is SPH-FEM. Although it was firstly applied to structure to structure impact problem (Attaway, Heinstein, and Swegle 1994), generally this method is used to solve FSI problem in which SPH is used to simulate the fluid domain and FEM is used to simulate the solid domain. Main aim is to prevent the overlap of domains in FSI problems, many methods were proposed to prevent the overlap. In master slave algorithms, SPH particles are prevented to penetrate into finite element meshes by applying contact forces (Groenenboom and Cartwright 2010;Zhang et al. 2011). In some of them, finite element nodes are treated as SPH particles as boundary (Fernández-Méndez et al. 2005). However, in the proposed method, SPH particles in contact with the structure is solved together with the structure. This is satisfied by contact mechanics. The main difference of the proposed method is solving invading SPH particles together with finite elements of structure.
In the literature, some researchers deal with determination of sloshing effect (Rawat et al. 2019) and some determined the baffle effect on sloshing (Cheng, Jing, and Gong 2019;Sanapala et al. 2018;Zheng et al. 2018). In some researches, dam break problems were investigated Dinçer et al. 2019; and in some others submerged elastic structures were investigated (Ghalandari et al. 2019). While some researchers were interested in a part of a structure submerged in water (Liang et al. 2001) or under effect of fluid (Yu et al. 2019), some others simulated structure completely (Ding et al. 2018;Liaw and Chopra 1974;Williams and Moubayed 1990). Problem classification can be made for FSI problems in many ways: type of structure, fluid, loading or even submergence of the structure. Except the closed form solutions and some methods having theoretical simplifications, a fully coupled FSI method can solve all type of problems, properly.
There are many submerged structures like relics or newly developed technologic structures like Submerged Floating Tunnels -SFT (Jin and Kim 2018) or offshore towers. Although, these structures are supported from ground, they are mostly affected from the fluid under a base excitation. This is the first time a recently proposed fully coupled SPH-FEM based FSI method is used for hydro-elastic analysis of a submerged structure under different base excitations. Smoothed particle hydrodynamics (SPH) is used for simulation of the fluid part and finite element method (FEM) is used for solid part of the model. Coupling between SPH and FEM is implemented via contact mechanics. In the scope of this research, series of experiments are carried out for defined problem and detailed numerical results are shown for structural part of the problem with experimental results.

NUMERICAL METHODS
Increasing number of type of methods and combinations of them make it possible to simulate an FSI problem in various ways. Indeed, finite element method is indispensable of them. However, nowadays it loses its popularity against smoothed particle hydrodynamics for analysis of fluid domains. Researchers uses particle-based methods like SPH or MPS even for crack propagation problems for solid structures. Nevertheless, FEM preserves its significance among all other methods. In this paper, most preferential methods among researchers are used to simulate the FSI problem which was proposed by author in his recent studies Dinçer et al. 2019). Method was verified in these researches with many benchmark problems. In this method, solid part of the problem is modelled by FEM, fluid part of the problem is modelled by SPH method and coupling between solid and fluid part is satisfied via contact mechanics.

Condensed form of governing equation of motion for a typical finite element is defined in
in which, mass M , damping C , and stiffness matrix K and force vector are defined below.
where  is the density, N is the volume-displacement interpolation matrix,  is the mass proportional definition constant, L B is the linear strain-displacement transformation matrix, E is the constitutive matrix, S is the matrix of second Piola-Kirchhoff stress, NL B is the non-linear strain-displacement transformation matrix and f is the body force vector. More detailed derivation can be followed from meritorious works (Bathe 2005). The time domain is discretized using Wilson-θ method (Wilson, Farhoomand, and Bathe 1972) in Eq. (2).
( 2) where effective stiffness matrix and effective force vector are defined below.
where  is a constant equal to , is the time increment, is the constant for Wilson-θ method which is ideally taken as 1.42. Accordingly, velocity and acceleration are achieved by Eq.
(3) and Eq. (4), respectively and iterative solutions are carried out for next time step. (3) Rayleigh damping theory is used to determine the damping matrix C ).

Fluid
SPH is a particle based Lagrangian approach ). Kernel function is used to transform the partial differential equations into ordinary differential equations. Time variation of density , velocity i V and position of i th particle is given through Eq. (5) to Eq. (7).
Latin American Journal of Solids and Structures, 2020, 17(9), e333 4/14 where is the number of neighbour particles, pressure is defined depending on initial pressure and initial density 0  as in Eq. (8), i and j denotes the neighbor particles, m is the mass, ij V is the velocity difference between neighbor particles, is the body force including gravitational force, is the contact force applied by the structure, and W is the cubic spline kernel defined in Eq. (9) in two dimensional. (8) where is the smoothing length taken as 1.33 times the initial particle spacing, ij r is the position difference between neighbor particles. Actually, time variation of position of a particle should be the velocity of the particle, however a velocity stabilizer is added as shown in Eq. (7) in addition to velocity of the particle. An artificial viscosity , which is defined in Eq. (10), is also added to the pressure terms as shown in Eq. (6) for stabilizing the numerical algorithm (Monaghan 1992).
where and are empirical coefficients taken as 0.2 and 1, respectively. is the artificial viscosity which is defined in Eq. (11) and c is the artificial sound speed.
The leap-frog algorithm is used for time stepping and maximum step size is calculated from the Courant-Friedrichs-Lewy (CFL) condition (Anderson 1995;Hirsch 2007). Free slip boundary condition on fixed walls is satisfied by using mirror particle method. The boundary between solid and fluid domain is satisfied by contact mechanics defined as follows.

Coupling
Potential of a contact (Bathe and Chaudhary 1985;Dinçer et al. 2019) consisting of a node and a line is defined in Eq. (12). This potential is due to the invasion of boundary defined as a line in 2D geometry by the node.
where denotes the invader node, is the contact force at time , is the incremental contact force, where is the unit normal vector of the line element, l is the length of the line element, and is position vector of particle and node F at time , respectively.
A new potential is achieved by subtracting the summation of contact potential due to each invader particle from its existing potential . Substituting Eq. (12) into the new potential , incremental finite element equations including contact conditions is achieved as in Eq. (14).
where is the stiffness matrix of solid structure, is the contact stiffness matrix, is the incremental displacement vector, total applied external load vector, is the equivalent nodal forces vector, is the contact force, the elements of which are calculated using position of contact point which is defined by  .
At this stage, mass participation of invading particles should be added to the set of equations (Eq. (14)) to solve the structure with invaders together. Newmark method (Newmark 1959) is used for this purpose by taking constants α and β as 0.5 and 0.25, respectively. Accordingly, mass matrix of invaders is multiplied by before substitution.
More detailed derivation of a node to line  or line to line (Bathe and Chaudhary 1985) contact potential can be followed from previous researches

EXPERIMENTS
Experiments are designed to investigate the motion of standing submerged highly elastic structures. It is preferred to use structure with low Young's modulus. Because, more displacement will make interaction more difficult, so the method has been tested for harsh cases. Accordingly, a rubber plate in a tank filled with water is excited by different motions.
The setup is given in Figure 1(a). Dimensions of tank are 30 cm by 25 cm. Tank is made of plexiglas having 1 cm thickness which is enough for rigid solid assumption. Water height is fixed to 10 cm for all experiment cases. A thin rubber plate having 3 mm thickness and 8 cm height is fixed at the mid bottom of the tank. A high-speed camera is deployed at the center axis of the tank, in front of the rubber plate. Both camera and tank are placed on the uniaxial shake table. Photo of the setup is also given in Figure 1  Two types of excitations are applied to the setup. Both of them are transient motions. First one is pulse type motion or near fault type ground motion. The second one is a real earthquake motion. Numerical model of structure is discretized by 160 four node quadrilateral bilinear finite elements. Young's modulus of rubber plate is determined as 10 MPa. Poisson's ratio of it is taken as 0.3 and damping coefficient is determined as %8. Damping ratio is achieved by vibration tests. In this test, a tip displacement is applied to the same rubber plate and damped motion of it is investigated. Results of the vibration test is given in Figure 2. Period of the motion is determined as 0.114 s for both experimental and numerical analysis. An average damping ratio is calculated from the vibration experiment as %8 and this is used for simulation of vibration test and for other proceeding simulations. Density of rubber plate is 2000 kg/m3 whereas density of water is 1000 kg/m3. 13146 particles are used for fluid (water) part. Initial particle spacing for SPH simulations is taken as 1.5 mm. Time interval used for analysis is calculated from CFL condition as 0.016 ms. It should be noted that a convergence analysis has been conducted for all the simulations and an optimum number of finite elements and particles are used for solid and fluid domains, respectively.

Pulse type excitations
The characteristics of near fault excitations are different from far ones. The energy applied to the structure is much higher for near fault excitations. Therefore, near fault earthquakes are more destructive and hence interest of many researchers. These motions are characterized in many forms. In general, these motions are long acceleration pulses and classified into three as Type A, B and C1 (Makris and Black 2004a). Type A is one sine pulse, Type B and C1 are cosine pulses. Near fault earthquakes consist these types of pulses at the very initial stage of the motion. For example, Chi-Chi earthquake in 1999 consists an A type pulse and Erzincan earthquake in 1992 consists a B type pulse (Makris and Black 2004b).
In this case, three different pulse type excitations are applied to the setup. These pulses are Type A, B and C1. Acceleration, velocity and displacement of these pulses are given in Figure 3. Experimentally and numerically achieved tip displacement of rubber plate under excitation of three different pulses are given in Figure 5. Corresponding kinetic energy histories per total masses are given for each type of excitation in Figure 4. Although total applied energies (Fajfar and Krawinkler, 1992) of type B and C1 are equal and greatest kinetic energy is seen in type A excitation, greatest tip displacement is seen in type C1 excitation. This is mainly due to the time of energy transfer. If the period of pulse is equal to the period of the motion, amplitude of the motion will increase, which is a well-known phenomenon named as resonance.  Although, tip displacements occur at the same time for Type B and C1 as expected, Type A has delayed peaks. While, first five peak tip displacements of Type B and C1 occur at 0.16 s, 0.42 s, 0.69 s, 0.98 s and 1.30 s, type A has peak tip displacements at 0.26 s, 0.53 s, 0.85 s, 1.17 s, 1.48 s (see Figure 5). Accordingly, period of motion of rubber plate is 0.568 s for Type B and C1 and it is 0.608 s for Type A. Both of them are close to the period of the applied pulse type motions which is 0.6 s (see Figure 3). Free surface profiles are given in Figure 6 for each pulse at the time of first five peak tip displacements of rubber plate. Additionally, stresses at peak tip displacements of rubber plate are illustrated with the free surface profile of water in Figure 7.  In order to provide more comparable results for free surface profiles, water level histories at the walls of the tank are given for each near fault excitations in Figure 8. System is run once again for the rubber plate itself, which is a dry case for near fault excitations. History of tip displacement of dry rubber plate is given in Figure 9. All types have almost same peak tip displacement. Type B and C1 has same tip displacement history till to the 0.6 s however Type A is different from them. It is seen that displacements of dry case are approximately 1/30 of wet case. It shows that motion of standing elastic structure is mostly driven by fluid. It is a known fact that a baffle type structure placed at the mid-centre of the tank will affect the sloshing behaviour. However, how much is the effect for a structure having low elasticity? In order to provide the effect on sloshing, comparison of water levels at the walls of tank for baffled and un-baffled cases are given in Figure 10. It seen that initial peak water levels of both cases are almost same; however, water levels of baffled cases are lower than the un-baffled cases for progressing peaks. It should be noted that, first peak water levels occur during the applied excitations, however second peaks of Type A and B occur after the excitations while it happens during the excitation of Type C1. It can be generalized as oscillations are damped more for baffled case after the end of pulse type excitations. Nevertheless, it can't be said that elastic baffle has much effect on sloshing during pulse type excitations. Beside the comparison of peaks of water levels, it is seen that baffle shortens the period of the oscillations of water.

Earthquake excitations
Earthquakes are very transient motions. Although, some researchers try to classify them e.g. near fault ones (Makris and Black 2004b), a general characterization is not possible for them. In this case, one of the most destructive earthquake affecting the El Centro city in southern California in 1940 is selected to apply to the system. Acceleration velocity and displacement time histories of this ground motion are given in Figure 11. Displacement history of El Centro earthquake (Fan and Ahmadi 1990) as given in Figure 11c is applied to the shake table and tip displacement histories of rubber plate is achieved both experimentally and numerically ( Figure 12). Considering the assumptions made for numeric analysis and inevitable irregularities of experimental setup, good match is obtained even for a true-life earthquake. A more detailed comparison is illustrated in Figure 13.  System is run once again for the rubber plate itself, which is a dry case. It is seen that there is almost no motion except vibrations. Maximum vibration is recorded at 2.75 s as 2.4 mm. History of tip displacement of dry rubber plate is given in Figure 14.

CONCLUSIONS
In the scope of this research, a submerged standing structure is investigated both experimentally and numerically under seismic excitations to show effects of seismic excitations to these types of structures. Although, some general classifications as far or near field excitations can be made for earthquakes, true-life earthquakes are unclassifiable seismic excitations. Therefore, two types of seismic excitations are selected to apply in this research. First one is the near field excitation which is applied in three different variations as type A, B and C1. In all types, results are satisfying. It is seen that all types lead to approximately same tip displacement pattern. However, applied energies are equal for type A and B and greater for type C1. This shows that applied energy amount is not the state of destructivity for a submerged structure. When the near fault excitations are applied to water itself, it is seen that water levels at the walls of the tank are same with the baffled (rubber plate) case during excitation. Baffle damps the oscillation of water after the end of excitation. Second excitation is a real-life earthquake which razed El Centro city in 1940. Result is also satisfying for a real case. Although same PGA is applied to both El Centro earthquake and near fault excitations, maximum tip displacement of rubber plate under near fault excitations is greater than El Centro earthquake.
Although, a submerged structure supported at the ground receives the excitation directly from the ground, not from the water (indirectly), motion of it is not independent of water. In fact, there will be even no tip displacement other than small vibrations for same setup with no water, even for a rubber plate having low Young's modulus as in this paper. It shows that interaction between water and structure forms the motion of submerged structure and a fully coupled FSI method as proposed in this paper should be used to analyse the system.
A fictious 2D FSI problem was investigated in this research to lead new studies in this area. The geometry of the tank and the structure can be reformed for further studies. Besides, data achieved from the submerged structure will be hard but precious.