## 2 Numerical Model and Verification  
This study uses the Lagrangian-Lagrangian approach to model the dynamic behavior of liquid storage tanks, which is based on the total energy and stain energy of liquids. It has been demonstrated that this method is capable of modeling fluid behavior in rigid tank systems [52; 54; 55]. This method is applied in this study by taking into account the flexibility of the tank's wall. The fluid potential energy (\(U\)), in this approach, is equal to the sum of strain energy and increasing potential energy (\(\Pi_{\rm s}\)), as shown in Equation (1).  
\[U=\Pi_{\epsilon}+\Pi_{\rm s} \tag{1}\]  
where \(\Pi_{\epsilon}\) and \(\Pi_{s}\) are the strain and potential energies, respectively.  
The strain energy equation is defined as follows:  
\[T=\frac{1}{2}\!\int\rho v^{2}vdv \tag{2}\]  
In Equation (2), \(v\) is the velocity vector, \(\rho\) is the mass density of the fluid, and \(T\) is the strain energy. The Lagrange equation can be therefore written as follows:  
\[\frac{d}{dt}\left(\frac{dT}{\dot{u}_{j}}\right)-\frac{\partial T}{\partial u_ {j}}+\frac{\partial U}{\partial u_{j}}=F_{j} \tag{3}\]  
where \(u_{j}\) and \(F_{j}\) are the displacement and force vectors of the \(j\)th component.  
If \(S\) is the rigidity matrix for the surface elements, and \(K\) and \(M\) are the stiffness and mass matrixes of the fluid elements, Equations (1) and (2) can be rewritten as Equations (4) and (5):  
\[U=\Pi_{\epsilon}+\Pi_{s}=\frac{1}{2}u^{T}Ku+\frac{1}{2}u_{s}^{T}Su_{s} \tag{4}\]  
\[T=\frac{1}{2}v^{T}Mv \tag{5}\]  
The rigidity matrix (\(S\)), stiffness matrix (\(K\)), and mass matrix (\(M\)) are given below:  
\[M=\rho\int Q^{T}QdV \tag{6}\]  
\[K=\int B^{T}EBdV \tag{7}\]  
\[S=\rho g\!\int{Q_{s}}^{T}Q_{s}dA \tag{8}\]  
where \(Q\) and \(Q_{s}\) are the interpolation functions for three-dimensional elements and two-dimensional surface elements, respectively. In fluid meshes, the interpolation matrix maps degrees of freedom to degrees of freedom in structural meshes. In FSI simulations, information is transferred between fluid and structural domains by the interpolation matrix.  
If Equations (4) and (5) are substituted in the Lagrangian equation (Equation (3)), the governing equation of motion can be written as follows (Equation (9)):  
\[Ma+Ku+Su_{s}=R \tag{9}\]  
where \(a\) is the acceleration vector of the nodes at the structural domain's boundary elements, and \(R\) is the load vector, which generally varies in each time interval.  
To calculate the pressure and displacement of the system and their derivatives, the time-stepping Newmark method is employed, as shown in Equations (10) and (11).  
\[\begin{bmatrix}\dot{U}\\ \dot{\rm P}\end{bmatrix}_{i+1}=\begin{bmatrix}\dot{U}\\ \dot{\rm P}\end{bmatrix}_{i}+[(1-\gamma)\Delta t]\begin{bmatrix}\ddot{U}\\ \ddot{\rm P}\end{bmatrix}_{i}+(\gamma\Delta t)\begin{bmatrix}\ddot{U}\\ \ddot{\rm P}\end{bmatrix}_{i+1} \tag{10}\]\[\begin{bmatrix}U\\ P\end{bmatrix}_{i+1}=\begin{bmatrix}U\\ P\end{bmatrix}_{i}+(\Delta t)\begin{bmatrix}\dot{U}\\ \dot{\mathrm{p}}\end{bmatrix}_{i}+[(0.5-\beta)(\Delta t)^{2}]\begin{bmatrix} \ddot{U}\\ \ddot{\mathrm{p}}\end{bmatrix}_{i}+[\beta(\Delta t)^{2}]\begin{bmatrix}\ddot{U} \\ \ddot{\mathrm{p}}\end{bmatrix}_{i+1} \tag{11}\]  
In Equations (10) and (11), \(\beta\) and \(\gamma\) are the Newmark parameters allowing engineers to control the trade-off between numerical stability and accuracy. The Newmark formulation is implemented in the ANSYS Mechanical Release 13.0. The typical value for \(\gamma\) is 0.5, and for \(\beta\) it is between 0.17 and 0.25 [56].