## 2 SPH Methods
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 \(\alpha\) and \(\beta\) 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\({}_{\mathrm{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 [PERSON] and [PERSON] [19] from the form proposed by [PERSON] et al. [22]:  
\[\rho_{\text{a}}=\sum_{\text{b=1}}^{\text{N}}m_{\text{b}}W_{\text{ab}}^{\text{ MLS}} \tag{11}\]  
where \(W_{\text{ab}}^{\text{MLS}}\) is defined by:  
\[W_{\text{ab}}^{\text{MLS}}=(\varphi_{0}+\varphi_{1}x_{\text{ab}}+\varphi_{2}y_ {\text{ab}})W_{\text{ab}} \tag{12}\]  
\[ \begin{bmatrix}\varphi_{0}\\ \varphi_{1}\\ \varphi_{2}\end{bmatrix} =\left(\sum_{\text{b=1}}^{\text{N}}W_{\text{ab}} \begin{bmatrix}1&x_{\text{ab}}&y_{\text{ab}}\\ x_{\text{ab}}&x_{\text{ab}}^{2}&x_{\text{ab}}y_{\text{ab}}\\ y_{\text{ab}}&x_{\text{ab}}y_{\text{ab}}&y_{\text{ab}}^{2}\end{bmatrix}  \right)^{-1} \begin{bmatrix}1\\ 0\end{bmatrix}  \tag{13}\]  
where \(x_{\text{ab}}=x_{\text{a}}-x_{\text{b}}\), \(y_{\text{ab}}=y_{\text{a}}-y_{\text{b}}\). 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:  
\[\rho^{\text{n}+1/2}=\rho^{\text{n}-1/2}+\frac{\Delta t^{\text{n}-1}+\Delta t^ {\text{n}}}{2}\bigg{(}\frac{D\rho}{Dt}\bigg{)}^{\text{n}} \tag{14}\]  
\[\mathbf{u}^{\text{n}+1/2}=\mathbf{u}^{\text{n}-1/2}+\frac{\Delta t^{\text{n}- 1}+\Delta t^{\text{n}}}{2}\bigg{(}\frac{D\mathbf{u}}{Dt}\bigg{)}^{\text{n}} \tag{15}\]  
\[\mathbf{r}^{\text{n}+1}=\mathbf{r}^{\text{n}}+\Delta t^{\text{n}}\mathbf{u}^{ \text{n}+1/2} \tag{16}\]  
where n is the number of time steps. The time step \(\Delta t\) is constrained by the \(CFL\) condition as follows:  
\[\Delta t\leq CFL\cdot\underset{\text{a}}{\text{min}}\bigg{(}\frac{\text{h}_{ \text{a}}}{\text{c}_{\text{a}}+|\mathbf{u}_{\text{a}}|}\bigg{)} \tag{17}\]  
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, [PERSON] [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. \[\mathbf{\tilde{u}_{\text{a}}}=\mathbf{u}_{\text{a}}+\Delta\mathbf{u}_{\text{ a}}\triangleq\mathbf{u}_{\text{a}}+\varepsilon\sum_{\text{b=1}}^{\text{N}}m_{ \text{b}}\frac{\mathbf{u}_{\text{ab}}}{\overline{\rho}_{\text{ab}}}W_{\text{ab}} \tag{18}\]  
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 [PERSON] 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.