## 2 Description and Numerical Model of the Foil under Waves
### Numerical Method and Model of the Foil under Waves  
In these numerical simulations, the flow field is considered viscous and incompressible. Therefore, the continuity and momentum equations are:  
\[\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x_{i}}(\rho u_{i})=0,\,(i=1,2\text{ and }\rho=\text{constant}) \tag{11}\]  
\[\frac{\partial(\rho u_{i})}{\partial t}+\frac{\partial}{\partial x_{j}}\big{(} \rho u_{i}u_{j}\big{)}=-\frac{\partial p_{r}}{\partial x_{i}}+\frac{\partial }{\partial x_{j}}\Bigg{(}\mu\frac{\partial u_{i}}{\partial x_{j}}-\rho\overline {u^{\prime}{}_{i}u^{\prime}{}_{j}}\Bigg{)}+S_{i},\,(i,j=1,2\text{ and }\rho=\text{constant}) \tag{12}\]  
where \(u_{i}\) and \(u_{j}\) are the transient velocity components, while \(u_{i}{}^{\prime}\) and \(u_{i}{}^{\prime}\) are the fluctuating velocity components. \(p_{r}\) is the transient pressure and \(\overline{u^{\prime}_{i}u^{\prime}_{j}}\) represents the mean value of \(u^{\prime}_{i}u^{\prime}_{j}\). \(S_{i}\) is the source item and is the sum of all other force components that do not belong to the transient term, convective term, or diffusive term of the governing equations. \(\mu\) is the kinematic viscosity coefficient. In addition, the subscripts \(i,j=1,2\) in Equations (11) and (12) identify the components along the \(x\) and \(y\) axis, respectively. The standard \(k\)-\(\epsilon\) turbulence model, which has been commonly used, has some deficiencies when applied to the simulation of strongly swirling flow, curved wall flow, and curved streamline flow [34]. The renormalization group (RNG) \(k\)-\(\epsilon\) two-equation model is applied in this paper. This scheme differs from the standard \(k\)-\(\epsilon\) turbulence scheme, as it includes an additional sink term in the turbulence dissipation equation to account for non-equilibrium strain rates and employs different values for the model coefficients [35]. Therefore, the RNG \(k\)-\(\epsilon\) model can deal with the flows with high strain rate and large streamline curvature, such as the flow field of waves and around foils. The length of the NWT is \(L=5\lambda\). The computational domain is divided in several zones, including air zone, stable waveform zone, and artificial damping zone, as shown in Figure 2. The ratio of the stable waveform zone length to the artificial damping zone length is 3:2 and the foil is placed at 2/5 the stable waveform zone length from the left boundary of the NWT. The boundary conditions of the computational domain is also given in Figure 2. The mesh in the whole computational domain is shown in Figure 3. The stable waveform zone is modelled with an unstructured grid in order to validate the dynamic mesh and re-mesh technologies to ensure good-quality mesh during foil motion. Meanwhile, a structured grid is adopted to model the air domain and the artificial damping domain to improve solving efficiency. The grid near the foil and the free wave surface is densified to satisfy computational accuracy. The meshes of the air zone, stable waveform zone, and artificial damping zone have 12,000, 179,000, and 2000 elements, respectively. The mesh scheme was selected based on reference [30] and a sensitivity analysis was also performed. We selected the fine mesh to obtain a more accurate flow field. In Figure 3, a regular wave is produced from the water inlet boundary by defining the inflow velocity component. Wave generation equations are implemented in the Fluent solver by programming a UDF. In addition, the height of the air zone is \(h_{u}=\lambda/4\). The water depth of the NWT is \(h_{w}=3\lambda/4\) and it is deep enough that the tank bottom has no influence on the surface wave profile. Therefore, the velocity potential of 2D regular waves in our simulations can be given as:  
\[\phi=\frac{qg}{\omega_{w}}e^{ky}\sin(kx-\omega_{w}t) \tag{13}\]  
where \(\omega_{w}\) is the circular frequency of the wave's simple harmonic motion, \(\omega_{w}=(gk)^{0.5}\).