Smoothed Particle Hydrodynamics Simulation of the Dynamic Water Pressure inside an Arch Rib of an Arch Bridge Subjected to Seismic Excitation

: The dynamic pressure exerted on an arch rib by inner water under seismic excitation was studied based on the smoothed particle hydrodynamics (SPH) numerical model. Numerical cases of different seismic excitations in both the horizontal and vertical directions and different water depths were carried out. The fast Fourier transform (FFT) of dynamic pressure exerted on the arch rib showed that the spike frequency is related to the water depth. The surface wave frequency of inner water is very close to the ﬁrst spike frequency. The ﬁrst spike frequency comes close to the second and fourth vibration frequencies when the inner water reaches the depth of 18 m and 11.3 m, respectively. This will likely lead to the resonance of the entire arch rib. By considering barely the ﬁrst spike frequency, when designing a new arch bridge, the possibilities of resonance in its second and fourth modes under seismic excitation cannot be ruled out. Based on the ﬁndings of this investigation, the inner water effect should be considered during the arch bridge design stage.


Introduction
Arch bridges are the main bridge type in China. As of 2019, there are 65 arch bridges with a span of 300 m or more in the world, of which 35 are in China. Among all the arch bridges in China, about half of the arch bridges have a deck arch with a boxing ring structural form. This is an economical type of bridge in mountainous areas and can preferably adapt to the topography of mountainous areas [1]. However, the possibility and intensity of earthquakes are high in mountainous areas, especially in southwestern China ( Figure 1, cited from the China Earthquake Administration GB18306-2001, China Standard Press). Earthquake damage caused to a bridge can have severe consequences. Obviously, the collapse of a bridge places people on or below the bridge at risk. A bridge closure, even if temporary, can have tremendous consequences because bridges often provide vital links in a transportation system. It is necessary to consider bridge performance subjected to seismic excitation when we design arch bridges in mountainous areas.
Generally, the arch rib of this boxing ring structural arch bridge is hollow, and water can enter the arch rib in order to balance the water pressure and reduce self-gravity. There is limited literature on the inner water of an arch bridge. Some authors studied the earthquake-induced added hydrodynamic pressure of the inner water in a hollow bridge pier. Lai et al. [2] studied the hydrodynamic pressure due to the outer and inner water effect on a hollow cylindrical pier based on radiation wave theory. Liu et al. [3] found that the inner water in a hollow bridge pier has similar negative effects on the piers as the outer water around the piers under seismic excitation. Li and Yang [4] studied the inner water in hollow piers in deep waters under earthquakes. They believed the free surface wave and its result of added dampening was negligible. However, for an arch bridge, due to the complexity of the hollow arch rib, the inner water effect is more complicated than the inner water effect in a hollow pier. In extreme situations, when an earthquake happens, does the possibility of damage to the arch bridge increase due to the existence of inner water in the arch rib? This should be a potential topic of the seismic damage caused to arch bridges. Intuitively, the vibration of inner water in the arch rib under seismic excitation can be regarded as a sloshing problem. The sloshing dynamic pressure induces the vibration of the arch rib and may cause resonance, which may damage the bridge. In order to study this sloshing problem, a CFD method-the SPH method-was employed for the numerical simulation. Appl. Sci. 2022, 12, 10626 2 of 23 in hollow piers in deep waters under earthquakes. They believed the free surface wave and its result of added dampening was negligible. However, for an arch bridge, due to the complexity of the hollow arch rib, the inner water effect is more complicated than the inner water effect in a hollow pier. In extreme situations, when an earthquake happens, does the possibility of damage to the arch bridge increase due to the existence of inner water in the arch rib? This should be a potential topic of the seismic damage caused to arch bridges. Intuitively, the vibration of inner water in the arch rib under seismic excitation can be regarded as a sloshing problem. The sloshing dynamic pressure induces the vibration of the arch rib and may cause resonance, which may damage the bridge. In order to study this sloshing problem, a CFD method-the SPH method-was employed for the numerical simulation. SPH is a mesh−free method that was originally developed for astrophysical simulations by Lucy [5] and Gingold and Monaghan [6], and was later extended to model fluid flows by Monaghan [7]. SPH is a pure Lagrangian particle method that does not need a grid to discretize the spatial derivatives. The main advantage of this method is that it prevents severe mesh distortions caused by large deformation. The other advantage of SPH is that there is no need to deal with the free surface owing to the Lagrangian nature of the SPH method. This property makes it particularly attractive to model free surface flows, especially when the surface tension is not important, such as in sloshing problems. Colagrossi et al. [8,9] applied SPH to sloshing simulation. Delorme et al. [10] analyzed the characteristics of impact pressure with a rectangular tank in a rolling motion. Colagrossi et al. [11] studied the violent sloshing phenomenon using an improved SPH method with MLS integral interpolators. Marsh et al. [12] studied liquid sloshing in different−shaped containers. Shao et al. [13] proposed a new coupling boundary method to simulate rectangular tanks in the rolling or swaying motion. Gotoh et al. [14] studied violent sloshing by using an enhanced incompressible SPH method. The SPH method is able to compute the dynamic force on structures [15][16][17] and even soil structure interaction [18]. SPH is a mesh-free method that was originally developed for astrophysical simulations by Lucy [5] and Gingold and Monaghan [6], and was later extended to model fluid flows by Monaghan [7]. SPH is a pure Lagrangian particle method that does not need a grid to discretize the spatial derivatives. The main advantage of this method is that it prevents severe mesh distortions caused by large deformation. The other advantage of SPH is that there is no need to deal with the free surface owing to the Lagrangian nature of the SPH method. This property makes it particularly attractive to model free surface flows, especially when the surface tension is not important, such as in sloshing problems. Colagrossi et al. [8,9] applied SPH to sloshing simulation. Delorme et al. [10] analyzed the characteristics of impact pressure with a rectangular tank in a rolling motion. Colagrossi et al. [11] studied the violent sloshing phenomenon using an improved SPH method with MLS integral interpolators. Marsh et al. [12] studied liquid sloshing in different−shaped containers. Shao et al. [13] proposed a new coupling boundary method to simulate rectangular tanks in the rolling or swaying motion. Gotoh et al. [14] studied violent sloshing by using an enhanced incompressible SPH method. The SPH method is able to compute the dynamic force on structures [15][16][17] and even soil structure interaction [18].
Although SPH has already been applied widely in sloshing simulation and the calculation of the dynamic force on structures, it should be noted that there is no application in the simulation of inner water movement and the dynamic pressure exerted on the arch Appl. Sci. 2022, 12, 10626 3 of 22 rib by using SPH. In this paper, we will explore the SPH application for the arch bridge inner water problem. In order to obtain enough accurate pressure on boundary particles, a boundary treatment similar to the treatment approach of [17] is adopted.

SPH Methods
The movement of inner water in an arch rib under seismic vibration is a sloshing problem. The viscosity of water can be ignored in this study. The governing equations of fluid dynamics can be reduced from the Navier-Stokes (N-S) equations to Euler equations. SPH is a purely Lagrangian method, so the following Lagrangian form continuity and momentum equations are often used in SPH representations: where ρ is the fluid particle density; t is the time; u is the particle velocity vector; p is the pressure; and g is the gravitational acceleration; ∇ = i ∂ ∂x + j ∂ ∂y is the Nabla operator, and (i, j) is the orthogonal basis of 2D Cartesian coordinates.
There are various forms of the SPH formulation for the above governing equations, and they demonstrate both advantages and disadvantages depending on the different hydrodynamic applications. In this paper, we adopted a revised SPH representation proposed by Colagrossi and Landrini [19]. The SPH continuity equation and momentum equation are: where a is the number of reference particles and b represents the number of neighboring particles; N is the total number of neighboring particles of the reference particle a; m is the particle mass; u is the velocity vector; u ab u a − u b is the velocity difference between the reference particle a and its neighboring particle b and ∇ a W ab ∇ r W(|r − r b |, h)| r=r a represents the gradient of a kernel taken with respect to a particle a. Here, W is the kernel function for SPH interpolation. In this work, a modified Gaussian kernel proposed by Colagrossi and Landrini [19] is adopted: where s |r − r |; r (x, y) represents the coordinate vector of the reference particle; r is the coordinate vector of the neighboring particle; h is the smoothed length, which is a measure of the support domain of the kernel function and is taken as a multiple of the initial particle spacing; and δ represents the cut-off limit due to the original Gaussian function is not compact support. Here, we choose δ = 3h as the compact support distance. According to our computational experience, this kernel has the advantage of good accuracy, better stability, and reasonable computational cost.
In this work, the standard weakly compressible SPH method (WCSPH) is adopted. In the WCSPH, the fluid pressure is computed using the following equation of state: where γ = 7 and ρ 0 = 1000 kg/m 3 are used for the water. The constant B is related to the sound speed c, which should be ten times or larger than the highest bulk flow velocity. It can be estimated by B = ρc 2 /γ. In the SPH numerical scheme, Π ab in Equation (4) is an artificial viscosity term that increases the stability properties of the algorithm. Here, the artificial viscosity has the following form: where α is a dimensionless coefficient; ρ ab (ρ a + ρ b )/2 is the average value of the density of particles a and b; c ab (c a + c b )/2 is the average value of the sound speed for the two particles; r ab r a − r b is the difference in spatial locations between the two particles; and µ ab is calculated by a modified form of [19] as follows: where is a small number to prevent singularity; k ab (k a + k b )/2, and k is defined as: where S = sym(∇u) is the rate−of−strain tensor and S : S is the tensor operation.
Since SPH can be affected by a lack of stability, introducing artificial viscosity can help reduce this problem, considering an artificial viscosity coefficient term [20] as shown in Equation (10): where c,ρ and h are sound speed, density, smoothing length, and kernel function, respectively, and the subscript 'ij' represents a mean value. A preliminary analysis was carried out to obtain the best set−up of parameters α, β and ε.
As regards artificial viscosity parameters, a series of calculations with different settings was carried out, covering the range of α and β from 0.3 to 0.01; in particular, the values of the two parameters were always kept equal to each other, and it was decided to modify the value of ε by the same percentage, starting from a value equal to 0.1 in correspondence with the maximum value of α and β. In this case of sloshing calculations, satisfactory results were obtained in correspondence to values of α and β lower than 0.1. In contrast, for higher values, calculation results were clearly dampened without capturing any impact phenomenon. In particular, it was found that the best results were obtained with values of α and β set to 0.05 and 0.01, and it was chosen to adopt the latter, which corresponds to low artificial viscosity.
A relatively large number of computations were performed using SPH to analyze the effect of some parameters, particularly artificial viscosity. The results of this analysis are not included in this paper, and only the final set-up results are presented (artificial viscosity parameters set to 0.01). As mentioned earlier, the influence of artificial viscosity is very pronounced, and in particular higher values result in an overdamping of all phenomena with lower pressure peaks (with impact pressures almost vanishing when values of α and β in artificial viscosity are set equal to 0.3) and without overturning waves.
Therefore, the artificial viscosity coefficient is set to a low value of 0.01, which dramatically reduces the artificial viscosity contribution.
In the SPH method, each particle has a constant mass m a . If the number of particles is constant, mass conservation is intrinsically satisfied. However, due to using the continuity Equation (3) for the density calculation, it is difficult to exactly enforce the consistency between particle mass, density, and its occupied area [21]. To alleviate this problem, the density field is periodically re-initialized by the following formula, which was modified by Colagrossi and Landrini [19] from the form proposed by Belytschko et al. [22]: (11) where W MLS ab is defined by: x ab y ab y ab x ab y ab y 2 In this work, the re-initialization procedure was applied every 20-time steps.
The SPH time integration is carried out based on the following frog-leap schemes as follows: where n is the number of time steps. The time step ∆t is constrained by the CFL condition as follows: where CFL is a dimensionless coefficient ranging from 0.0 to 1.0. In order to prevent particle inter-penetration and to regularize the particle motions, Monaghan [23] introduced the XSPH (X is an unknown factor) velocity correction, which takes into account the neighboring velocities through a mean evaluated within the particle support domain.ũ This corrected velocity is only used in the continuity Equation (3) and the position evolution Equation (15).
The accurate treatment of the solid boundary conditions is important for evaluating the pressure and force on the arch during vibration. Here, a boundary treatment proposed by Oger et al. [17] is adopted to calculate the pressure on the boundary particles. In this boundary treatment, the pressure of the boundary particles is obtained by interpolation from the pressure of fluid particles near the boundary particles.

Model Validation
To validate the SPH numerical model, a benchmark test was carried out by using our self−developed SPH simulation code. The computational results were compared with the experimental data of Brizzolara et al. [20].
In the experiments of Brizzolara et al. [20], the tank dimensions were 900 × 508 × 50 mm. Our simulation cases studied a longitudinal section of the tank. The computational water tank is shown in Figure 2. The tank was rolled along the rolling axis (184 mm over the bottom line of the tank) following a sinusoidal type rotating motion whose amplitude is • and with different excitation frequencies. Six sensors are located on the boundary of the tank as shown in Figure 2.
To validate the SPH numerical model, a benchmark test was carried out by using our self−developed SPH simulation code. The computational results were compared with the experimental data of Brizzolara et al. [20].
In the experiments of Brizzolara et al. [20], the tank dimensions were 900 × 508 × 50 mm. Our simulation cases studied a longitudinal section of the tank. The computational water tank is shown in Figure 2. The tank was rolled along the rolling axis (184 mm over the bottom line of the tank) following a sinusoidal type rotating motion whose amplitude is = 4° and with different excitation frequencies. Six sensors are located on the boundary of the tank as shown in Figure 2. Two water filling levels considered in the present analysis are shown in Table 1. For both water levels, the experiments of Brizzolara et al. [20] were carried out with the tank rolling in the resonance period (cases A and B). For case A, the impact pressure was recorded on the sides for sensor 1. For case B, pressure sensors were located at positions 3 and 6.
The SPH computational parameters were used as follows: equation of state coefficient B = 2.8 10 , artificial viscosity coefficient α = 0.01, and particle spacing ∆ = 2 mm. A constant time step of ∆t = 1 10 was used for the frog−leap time integration. For each of cases A and B, as the tank starts to roll, the main water wave generates gradually and travels from one side of the tank to the other. The water wave breaks during its traveling and impacts the wall of the tank. The experiment results show that the water motion in the tank was periodic despite the high dissipation due to wave breaking. Two water filling levels considered in the present analysis are shown in Table 1. For both water levels, the experiments of Brizzolara et al. [20] were carried out with the tank rolling in the resonance period (cases A and B). For case A, the impact pressure was recorded on the sides for sensor 1. For case B, pressure sensors were located at positions 3 and 6. The SPH computational parameters were used as follows: equation of state coefficient B = 2.8 × 10 5 , artificial viscosity coefficient α = 0.01, and particle spacing ∆x = 2 mm. A constant time step of ∆t = 1 × 10 −5 s was used for the frog-leap time integration.
For each of cases A and B, as the tank starts to roll, the main water wave generates gradually and travels from one side of the tank to the other. The water wave breaks during its traveling and impacts the wall of the tank. The experiment results show that the water motion in the tank was periodic despite the high dissipation due to wave breaking. Figure 3 compares some snapshots of water-free surfaces simulated by the SPH model with experimental pictures of case A. Figure 3a-c shows the initial stage of sloshing. SPH results show good agreement with the experimental pictures. Figure 3d-f show the first wave-breaking event. The jet forms and starts to break at 2.9 s. This plunge-type breaker entraps air and air-water mixture impacts on the wall subsequently. The shape and position of the breaker simulated by SPH are quite similar to the experimental snapshot. Figure 3g-i show the following three times the water wave impacted on the side wall of the tank. The shape and height of the water splash of the SPH numerical results accord quite well with the experimental results. ing. SPH results show good agreement with the experimental pictures. Figure 3d-f show the first wave−breaking event. The jet forms and starts to break at 2.9 s. This plunge−type breaker entraps air and air-water mixture impacts on the wall subsequently. The shape and position of the breaker simulated by SPH are quite similar to the experimental snapshot. Figure 3g-i show the following three times the water wave impacted on the side wall of the tank. The shape and height of the water splash of the SPH numerical results accord quite well with the experimental results.  To validate the accuracy of our SPH model, the SPH−computed pressure time series at the sensor points were compared with the experimental data of Brizzolara et al. [20] in Figure 4 for both cases A and B. Figure 4 shows that the simulated pressure time series of our SPH model are in good agreement with experimental data in both amplitude and phase of pressure oscillation, even after a dozen sloshing periods. To validate the accuracy of our SPH model, the SPH-computed pressure time series at the sensor points were compared with the experimental data of Brizzolara et al. [20] in Figure 4 for both cases A and B. Figure 4 shows that the simulated pressure time series of our SPH model are in good agreement with experimental data in both amplitude and phase of pressure oscillation, even after a dozen sloshing periods. Experimental and SPH results showed very good agreement in terms of free surface shape, the impact phenomenon, and impact pressure on the wall, which means that our self−developed SPH model is capable of simulating the sloshing problems and capturing the water impact pressure on the solid wall accurately.

Numerical Simulation and Analysis of Result
The SPH model was employed to study the influence of inner water on the arch bridge under seismic excitation. A sketch of the arch is shown in Figure 5. The length of the arch rib is 140 m, and the vector height of the arch ring is 23.333 m. The arch rib has a catenary curve shape. In order to reduce self−gravity and to balance the water pressure, the arch rib is designed to be hollow, and water can enter the rib. The designed water level is 15.56 m. Most of the time, the inner water exists in the arch rib. Here, we consider the arch rib as a 2D hollow rib as shown in Figure 5 because the width of the bridge is much larger than the depth of the arch rib box. As preliminary research on the inner water in the arch rib, the 2D simulation will reduce the computational cost drastically, and the 2D simulation results can help us obtain a conceptual understanding of the inner water effect on the arch rib. Experimental and SPH results showed very good agreement in terms of free surface shape, the impact phenomenon, and impact pressure on the wall, which means that our self-developed SPH model is capable of simulating the sloshing problems and capturing the water impact pressure on the solid wall accurately.

Numerical Simulation and Analysis of Result
The SPH model was employed to study the influence of inner water on the arch bridge under seismic excitation. A sketch of the arch is shown in Figure 5. The length of the arch rib is 140 m, and the vector height of the arch ring is 23.333 m. The arch rib has a catenary curve shape. In order to reduce self-gravity and to balance the water pressure, the arch rib is designed to be hollow, and water can enter the rib. The designed water level is 15.56 m. Most of the time, the inner water exists in the arch rib. Here, we consider the arch rib as a 2D hollow rib as shown in Figure 5 because the width of the bridge is much larger than the depth of the arch rib box. As preliminary research on the inner water in the arch rib, the 2D simulation will reduce the computational cost drastically, and the 2D simulation results can help us obtain a conceptual understanding of the inner water effect on the arch rib. The initial particle spacing is ∆ . The initial fluid particle pressure is set up following the static pressure distribution. The initial pressure on boundary particles is interpolated from the fluid particles as per the boundary treatment proposed by Oger et al. [17]. The density of the water is taken as 1000 kg/m3. The coefficient in the equation of state is B = 2.8 10 . The artificial viscosity coefficient is α = 0.01. The XSPH coefficient is ε = 0.1. In order to study the pressure evolution on the boundary, we defined some monitoring points on the solid boundary as shown in Figure 6. B1 to B9 were set on the right side of the arch rib, B10 to B18 on the left side of the arch rib, and were symmetrical about the y−axis.
The objective of this investigation was to find out whether there would be any possibility of damage caused to the arch bridge. Several studies [24,25] on the dynamic response analysis of arch bridges have shown that results with high accuracy can be obtained by using a smaller number of beam cells. As the arch bridge's vibration shape and frequency were analyzed using the FEM method, and the finite element model of the arch bridge was established by software (Figure 7), it is composed of 949 beam elements with 1089 nodes. One main arch circle is composed of 6 pieces of the arch box, the cross−section of the arch ring is shown in Figure 8, and the material of the arch circle is C40 concrete; its density is 2420-2440 kg/m³, modulus of elasticity is 2.0 10 Pa, Poisson's ratio 0.2, and damping ratio 0.05. Because the arch bridge adopts an equal cross−sectional suspension line box−shaped hingeless arch, the two ends of the arch circle are cemented to the abutment, so the boundary condition is cemented at both ends. The first−to fourth−order vibration shapes are shown in Figure 9. The frequency and the illustration of each vibration shape are listed in Table 2. In this paper, the dynamic pressure of the inner water is simulated in 2D. Therefore, the inner water's influence could only cause the vibration shape in the vertical plane, which means only second−and fourth−order vibration shapes will be considered in our study. The initial particle spacing is ∆x. The initial fluid particle pressure is set up following the static pressure distribution. The initial pressure on boundary particles is interpolated from the fluid particles as per the boundary treatment proposed by Oger et al. [17]. The density of the water is taken as 1000 kg/m 3 . The coefficient in the equation of state is B = 2.8 × 10 6 . The artificial viscosity coefficient is α = 0.01. The XSPH coefficient is ε = 0.1. In order to study the pressure evolution on the boundary, we defined some monitoring points on the solid boundary as shown in Figure 6. B1 to B9 were set on the right side of the arch rib, B10 to B18 on the left side of the arch rib, and were symmetrical about the y-axis.
The objective of this investigation was to find out whether there would be any possibility of damage caused to the arch bridge. Several studies [24,25] on the dynamic response analysis of arch bridges have shown that results with high accuracy can be obtained by using a smaller number of beam cells. As the arch bridge's vibration shape and frequency were analyzed using the FEM method, and the finite element model of the arch bridge was established by software (Figure 7), it is composed of 949 beam elements with 1089 nodes. One main arch circle is composed of 6 pieces of the arch box, the cross-section of the arch ring is shown in Figure 8, and the material of the arch circle is C40 concrete; its density is 2420-2440 kg/m 3 , modulus of elasticity is 2.0 × 10 10 Pa, Poisson's ratio 0.2, and damping ratio 0.05. Because the arch bridge adopts an equal cross-sectional suspension line box-shaped hingeless arch, the two ends of the arch circle are cemented to the abutment, so the boundary condition is cemented at both ends. The first-to fourth-order vibration shapes are shown in Figure 9. The frequency and the illustration of each vibration shape are listed in Table 2. In this paper, the dynamic pressure of the inner water is simulated in 2D. Therefore, the inner water's influence could only cause the vibration shape in the vertical plane, which means only second-and fourth-order vibration shapes will be considered in our study.

Results of Different Particle Sizes
First, the numerical simulations using different particle initial spacings ∆ were carried out. The simulation results were compared to evaluate the influence of the particle size. The seismic acceleration data of the 2008 Wenchuan earthquake (east-west component) as shown in Figure 10 were used as excitation data in the three test cases. Three particle initial spacings, i.e., ∆ = 0.05 m, . √ m, 0.025 m, were chosen for comparison.
The fluid particle numbers were 36,910, 73,098, and 147,700, respectively. The boundary particle numbers were 18,159, 25,665, and 36,287, respectively. A constant time marching step ∆t = 2 10 was used for all three cases. The simulation was executed at t = 10 s

Results of Different Particle Sizes
First, the numerical simulations using different particle initial spacings ∆ were carried out. The simulation results were compared to evaluate the influence of the particle size. The seismic acceleration data of the 2008 Wenchuan earthquake (east-west component) as shown in Figure 10 were used as excitation data in the three test cases. Three particle initial spacings, i.e., ∆ = 0.05 m, . √ m, 0.025 m, were chosen for comparison.
The fluid particle numbers were 36,910, 73,098, and 147,700, respectively. The boundary particle numbers were 18,159, 25,665, and 36,287, respectively. A constant time marching step ∆t = 2 10 was used for all three cases. The simulation was executed at t = 10 s using 500,000−time steps because here we only focus on the influence of particle spacing.

Results of Different Particle Sizes
First, the numerical simulations using different particle initial spacings ∆x were carried out. The simulation results were compared to evaluate the influence of the particle size. The seismic acceleration data of the 2008 Wenchuan earthquake (east-west component) as shown in Figure 10 were used as excitation data in the three test cases. Three particle initial spacings, i.e., ∆x = 0.05 m, 0.05 √ 2 m, 0.025 m, were chosen for comparison. The fluid particle numbers were 36,910, 73,098, and 147,700, respectively. The boundary particle numbers were 18,159, 25,665, and 36,287, respectively. A constant time marching step ∆t = 2 × 10 −5 s was used for all three cases. The simulation was executed at t = 10 s using 500,000-time steps because here we only focus on the influence of particle spacing.
The pressure fields of the three cases were compared at different time steps. The results show that the pressure fields for different particle sizes are quite similar. The quantitative comparison of pressure evolution at boundary point B1 to B18 shows that both the amplitude and the phase of pressure oscillation are almost the same for ∆x = 0.05/ √ 2 m and ∆x = 0.025 m, while the results of ∆x = 0.05 m are obviously different from those of the two other cases. Figure 11 shows the comparison of pressure evolution only for points B3, B5, and B8, which are at different water depth positions. The pressure fields of the three cases were compared at different time steps. The results show that the pressure fields for different particle sizes are quite similar. The quantitative comparison of pressure evolution at boundary point B1 to B18 shows that both the amplitude and the phase of pressure oscillation are almost the same for ∆ = 0.05/√2 m and ∆ = 0.025 m, while the results of ∆ = 0.05 m are obviously different from those of the two other cases. Figure 11 shows the comparison of pressure evolution only for points B3, B5, and B8, which are at different water depth positions.
It shows that the same results can be obtained using the medium particle size (∆ = 0.0354 m) compared with the finest particle size (∆ = 0.025 m). Therefore, in the following simulations, the medium particle size was adopted to reduce the simulation time.  Figure 11. Comparison of the influence of different particle sizes on the boundary pressure evolution at three boundary points.

Results of Different Seismic Excitations
In this section, different seismic excitations were used for simulation. In addition to the seismic acceleration data of the 2008 Wenchuan earthquake (east-west component) as shown in Figure 10, the 1940 El Centro earthquake data as shown in Figure 12 were also used for simulation. Thus, there are six cases using different seismic excitations in this part: Wenchuan EW, El Centro EW, El Centro NS, El Centro vertical, El Centro EW+vertical, and El Centro NS+vertical. All simulations were executed at t = 40 s using 2,000,000 time steps, which consumed approximately 160 CPU hours in the series run by using a Xeon E5−2697 (2.6 GHz) Workstation (14 cores). Figure 11. Comparison of the influence of different particle sizes on the boundary pressure evolution at three boundary points.
It shows that the same results can be obtained using the medium particle size (∆x = 0.0354 m) compared with the finest particle size (∆x = 0.025 m). Therefore, in the following simulations, the medium particle size was adopted to reduce the simulation time.

Results of Different Seismic Excitations
In this section, different seismic excitations were used for simulation. In addition to the seismic acceleration data of the 2008 Wenchuan earthquake (east-west component) as shown in Figure 10, the 1940 El Centro earthquake data as shown in Figure 12 were also used for simulation. Thus, there are six cases using different seismic excitations in this part: Wenchuan EW, El Centro EW, El Centro NS, El Centro vertical, El Centro EW + vertical, and El Centro NS + vertical. All simulations were executed at t = 40 s using 2,000,000 time steps, which consumed approximately 160 CPU hours in the series run by using a Xeon E5−2697 (2.6 GHz) Workstation (14 cores). In the cases of Wenchuan EW, El Centro EW, and El Centro NS, the excitation consists of only a horizontal component, while in the three other cases (El Centro vertical, El Centro EW + vertical, El Centro NS + vertical), the excitation consists of a vertical component. For the horizontal vibration cases, as shown in Figure 13, the pressure oscillations at the For the horizontal vibration cases, as shown in Figure 13, the pressure oscillations at the symmetrical points on the left and right sides are opposite in phase. We compared all the symmetrical points from B1&B10 to B9&B18, and the pressure oscillations are all opposite in phase for all of these cases. In Figure 13, we only show some of the symmetrical points at different depths and on both the upper and lower surfaces of the arch rib. This type of pressure oscillation will cause the arch rib to undergo vertical antisymmetric bending, which corresponds to the second-order vibration shape shown in Table 2. However, only in the vertical vibration cases (El Centro vertical), as shown in Figure 14, the pressure evolution at the symmetrical points is in the same phase. This type of pressure oscillation will cause the arch rib to undergo vertical symmetric bending, which corresponds to the fourth-order vibration shape shown in Table 2. For the cases including both horizontal and vertical excitations (El Centro EW + vertical, El Centro NS + vertical), both second-and fourth-order vibration shapes should be considered.
Appl. Sci. 2022, 12, 10626 16 of 23 symmetrical points on the left and right sides are opposite in phase. We compared all the symmetrical points from B1&B10 to B9&B18, and the pressure oscillations are all opposite in phase for all of these cases. In Figure 13, we only show some of the symmetrical points at different depths and on both the upper and lower surfaces of the arch rib. This type of pressure oscillation will cause the arch rib to undergo vertical antisymmetric bending, which corresponds to the second−order vibration shape shown in Table 2. However, only in the vertical vibration cases (El Centro vertical), as shown in Figure 14, the pressure evolution at the symmetrical points is in the same phase. This type of pressure oscillation will cause the arch rib to undergo vertical symmetric bending, which corresponds to the fourth−order vibration shape shown in Table 2. For the cases including both horizontal and vertical excitations (El Centro EW + vertical, El Centro NS + vertical), both second−and fourth−order vibration shapes should be considered.  In order to analyze whether the inner water pressure oscillation can cause resonance in the second−or fourth−order vibration shape, a fast Fourier transform (FFT) was carried out for the pressure oscillation data for all boundary monitoring points B1 to B18 for all six cases. Figure 15 shows the fast Fourier transform results of all six cases. The zero−frequency fraction indicates static pressure on the points, which has no influence on the arch rib. The frequencies of the first four spikes are in Figure 15, and cases are listed in Table 3. The frequency of the first spike is 1.125 for all cases. The subsequent spike frequencies are about 3, 5, and 7 times the first spike frequency, respectively, which means that all spikes will contribute to the possible resonance if the first spike frequency causes resonance in the second−or fourth−order vibration shape. Therefore, we only compared the first spike frequency with the natural frequency of the arch bridge in the second−or fourth−order vibration shapes. The first spike frequency (1.125 Hz) is different from the frequency of the second (0.941722 Hz) and fourth (1.661838 Hz) order vibration shapes. It seems there is no risk of resonance at the designed water level for all excitations. In order to analyze whether the inner water pressure oscillation can cause resonance in the second-or fourth-order vibration shape, a fast Fourier transform (FFT) was carried out for the pressure oscillation data for all boundary monitoring points B1 to B18 for all six cases. Figure 15 shows the fast Fourier transform results of all six cases. The zero-frequency fraction indicates static pressure on the points, which has no influence on the arch rib. The frequencies of the first four spikes are in Figure 15, and cases are listed in Table 3. The frequency of the first spike is 1.125 for all cases. The subsequent spike frequencies are about 3, 5, and 7 times the first spike frequency, respectively, which means that all spikes will contribute to the possible resonance if the first spike frequency causes resonance in the second-or fourth-order vibration shape. Therefore, we only compared the first spike frequency with the natural frequency of the arch bridge in the second-or fourth-order vibration shapes. The first spike frequency (1.125 Hz) is different from the frequency of the second (0.941722 Hz) and fourth (1.661838 Hz) order vibration shapes. It seems there is no risk of resonance at the designed water level for all excitations.  The results show that every spike was at the same frequency for all six cases without correlation with either excitation type or excitation frequency, which means the inner water is under free vibration, and the vibration frequency is the natural frequency of the inner water.
The simulation results show that inner free surface waves are generated due to seismic excitation. Figure 16 shows a snapshot of an inner water wave on the right side of the arch rib. The figure shows that there are a series of waves on the water's surface, and a large vortex formed in the inner water. After careful calculation, the wave frequencies for every case are shown in Table 4. It is shown that the inner water wave frequencies for all cases are 1.1-1.15 Hz, which is very close to the first spike frequency (1.125 Hz). This result also implies that the first spike frequency is the natural frequency of the inner water.  The results show that every spike was at the same frequency for all six cases without correlation with either excitation type or excitation frequency, which means the inner water is under free vibration, and the vibration frequency is the natural frequency of the inner water.
The simulation results show that inner free surface waves are generated due to seismic excitation. Figure 16 shows a snapshot of an inner water wave on the right side of the arch rib. The figure shows that there are a series of waves on the water's surface, and a large vortex formed in the inner water. After careful calculation, the wave frequencies for every case are shown in Table 4. It is shown that the inner water wave frequencies for all cases are 1.1-1.15 Hz, which is very close to the first spike frequency (1.125 Hz). This result also implies that the first spike frequency is the natural frequency of the inner water.
The natural frequency of the inner water depends on the depth of the water. Therefore, the simulation of different water depths was carried out and the frequency of pressure oscillation is investigated in the next section.    The natural frequency of the inner water depends on the depth of the water. Therefore, the simulation of different water depths was carried out and the frequency of pressure oscillation is investigated in the next section.  Figure 17 shows the fast Fourier transform of the pressure oscillation at all boundary monitoring points B1 to B18 for water depths of 13 m, 15 m, 17 m, 19 m, and 21 m. It shows that the first spike frequency becomes smaller as water depth increases. The inner water wave frequencies are calculated and listed in Table 5, which shows that the inner water wave frequency is also similar to the first spike frequency of FFT.  Figure 17 shows the fast Fourier transform of the pressure oscillation at all boundary monitoring points B1 to B18 for water depths of 13 m, 15 m, 17 m, 19 m, and 21 m. It shows that the first spike frequency becomes smaller as water depth increases. The inner water wave frequencies are calculated and listed in Table 5, which shows that the inner water wave frequency is also similar to the first spike frequency of FFT.  The expression of the first spike frequencies and the water depth values was obtained by multinomial fit as below and plotted out in Figure 18. And the data before and after the fit and the errors are shown in Table 6.
f (H/g) 1/2 = 4.10239(d/H) 2 − 9.22848(d/H) + 6.06132 (19) in which d is the water depth and H is the vector height of the arch ring. It seems reasonable to estimate by Equation (18) that the resonance in the second−order vibration shape may occur at a water depth of 18 m, and the resonance in the fourth−order vibration shape occurred at 11.3 m as far as the first spike frequency is concerned. However, since the real situation is quite complex, the simplified analysis in this paper shows merely the possibilities of resonance in extreme situations, which indicates that more attention should be paid to the influence of inner water in an arch rib.  The expression of the first spike frequencies and the water depth values was obtained by multinomial fit as below and plotted out in Figure 18. And the data before and after the fit and the errors are shown in Table 6.
/ / = 4.10239 / − 9.22848 / + 6.06132 (19) in which is the water depth and is the vector height of the arch ring. It seems reasonable to estimate by Equation (18) that the resonance in the second−order vibration shape may occur at a water depth of 18 m, and the resonance in the fourth−order vibration shape occurred at 11.3 m as far as the first spike frequency is concerned. However, since the real situation is quite complex, the simplified analysis in this paper shows merely the possibilities of resonance in extreme situations, which indicates that more attention should be paid to the influence of inner water in an arch rib.

Summary and Concluding Remarks
In this paper, we studied the dynamic pressure exerted on an arch rib by the inner water under seismic excitation using the SPH numerical model. First, the particle size independence was studied by simulating using three different initial particle sizes. Then, the numerical simulations of different seismic excitations in both horizontal and vertical directions were carried out. Finally, numerical simulations of different water depths were carried out. The FFT of dynamic pressure exerted on the arch rib shows that the spike frequency is independent of the seismic excitation and is related to the water depth. The spike frequency is the natural frequency of inner water. The inner water wave frequency is similar to the first spike frequency. By considering only the first spike frequency, the possibilities of resonance in second-and fourth-vibration shapes under seismic excitation are not ruled out.
More attention should be paid to the inner water effect on the arch bridge. In the future, the hydrodynamic pressure exerted on the 3D arch rib by inner water will be studied. Further studies should be carried out to find out the effects on superstructures or other parts of a bridge by the hydrodynamic pressure applied on the arch rib caused by inner water. In a conclusion, the influence of the inner water exerted on an arch rib under seismic excitation should be considered when designing an arch bridge. SPH can be a powerful tool for inner water hydrodynamic pressure studies during the arch bridge design stage.