# On the Efficacy of Turbulence Modelling for Sloshing
## 3 Methodology  
The multiphase VoF fluid solver interFoam, from the OpenFOAM v8 [60] code base is employed to simulate a vertically oscillated sloshing tank with harmonic motion. This solves the incompressible Navier-Stokes equations expressed as an arbitrary moving volume to account for moving boundaries [61; 62]. These are in the form  
\[\
abla\cdot\mathbf{U}=0, \tag{8}\]  
and  
\[\frac{\partial\rho_{f}\mathbf{U}}{\partial t}+\
abla\cdot[\mathbf{U}\cdot\rho_{f}(\bm {U}-\mathbf{U}_{s})]=-\
abla p+\
abla\cdot[(\mu+\mu_{t})\
abla\mathbf{U})]+\mathbf{s_{ \phi}}, \tag{9}\]  
where \(\rho_{f}\) is the fluid density, \(t\) is the time, \(\mathbf{U}_{s}\) is the velocity of the surface \(S\) (i.e., the grid velocity), \(\mu\) is the dynamic viscosity of the fluid, \(p\) is the fluid pressure and \(\mathbf{s_{\phi}}\) is a source term that includes both surface tension and gravity. The VoF method is used to capture the liquid-gas interface. In this approach, the density and viscosity are determined locally as a ratio of the occupying fluid:  
\[ \begin{split}\rho_{f}&=(1-\alpha)\rho_{g}+\alpha\rho _{l},\\ \mu&=(1-\alpha)\mu_{g}+\alpha\mu_{l}\end{split}  \tag{10}\]  
where the subscripts \(g\) and \(l\) refer to the properties of the gas and liquid (here, properties are for air and water). \(\alpha\) is the volume fraction of the liquid. The interface is constructed with the Multidimensional Universal Limiter for Explicit Solution (MULES) scheme, and artificial interface compressibility is employed to prevent excessive numerical diffusion of the phase fraction field [63]. In our work herein, the extent of the interface compression is set conservatively, such that spurious deformation is prevented, while still being effective in reducing interface numerical diffusion. The OpenFOAM keyword CAlpha controls the extent of this compression and is set to the widely adopted value of 1 in all simulations presented [37]. The turbulent viscosity \(\mu_{t}\), is computed as a function of the turbulence kinetic energy, \(k\), and the dissipation rate, \(\epsilon\), or the specific dissipation rate, \(\omega\). Additional transport equations are solved to approximate these values, with each having different forms in each turbulence model. Three models are tested here, the \(k\)-\(\epsilon\) model [64; 65], the Re-Normalization Group (RNG) \(k\)-\(\epsilon\) model [66] and the \(k\)-\(\omega\) shear-stress transport (SST) model [67]. These models are widely used, and are available in many commercial and open-source CFD codes. The models are well-described in the literature, and hence not presented in detail here. The damping functions at the modelled walls are used according to the original formulation of each model. To test the effectiveness of these models in providing an accurate representation of sloshing dynamics, no-model simulations that do not use a turbulence model are also conducted and compared with the experimental data. Within the chosen fluid solver, the oscillatory movement of the tank is imposed by displacing the domain using its dynamic mesh functionality and applying no-slip boundary conditions on the walls. In OpenFOAM, this is achieved via its dynamic mesh library dynamicMotionSolverFvMesh. There are two time-integration schemes available in interFoam, the first-order Euler scheme and the second-order Crank-Nicolson scheme. It is found that the Crank-Nicolson scheme causes instabilities and spurious oscillation in the volume fraction \(\alpha\), and hence introduces an unrealistic oscillation in the computed wall stresses. The Euler scheme produces smoother results compared to the higher order method. Similar observations are reported in [32]. Therefore, the first-order Euler scheme is used in the simulations presented, and is combined with a variable time step size and a maximum Courant number of 0.25. Regarding spatial discretisation, a uniform hexahedral mesh is used and three mesh resolutions are considered.