# Numerical Study on Airfoil Aerodynamics in Proximity to
## 2 Numerical Method and Validation  
In this research, the NACA4412 airfoil of the chord length \(c=0.5\) m is investigated. The freestream velocity \(V_{\infty}\) is 30 m/s, and the Reynolds number \(Re\) based on the chord length is \(3\times 10^{5}\). The lift coefficient of the airfoil in proximity to the wavy water surface is the main concern, and it is calculated at four angles of attack: \(\alpha=0^{\circ}\), \(1^{\circ}\), \(2^{\circ}\) and \(4^{\circ}\). The calculation domain is presented in Figure 1. The inlet, upper, and lower boundaries are positioned \(20c\) away from the airfoil, and the outlet boundary is located at \(25c\) from the trailing edge of the airfoil. A pressure outlet boundary condition is set for the upper and outlet boundaries. The inlet boundary and lower boundary are set as velocity inlet condition and nonslip wall condition respectively. In addition, the airfoil surface is set as nonslip wall condition. Figure 1 also gives the definition of the ride height \(h\) and wave parameters. The ride height is defined as the distance from the trailing edge of the airfoil to the still water level and set as \(0.2c\). The wavelength \(\lambda\) is \(5c\) which refers to existing studies and for each angle, the wave amplitude \(a\) varies from \(0.05c\) to \(0.1c\). The numerical wave is generated by using the volume of fluid (VOF) technique [28]. Based on the first-order linear wave theory [29], the wave propagating along the \(x\)-axis is generated by applying the inflow wave making method based on the stream function [30]. The wave surface elevation and the velocity components are specified at the velocity inlet, and they can be calculated as  
\[\eta=a\cos(kx-\omega t) \tag{1}\]  
\[V_{x}=a\omega\frac{\cosh k(\eta+d)}{\sinh kd}\cos(kx-\omega t) \tag{2}\]  
\[V_{z}=a\omega\frac{\sinh k(\eta+d)}{\sinh kd}\sin(kx-\omega t) \tag{3}\]  
where \(k=\) wave number, with \(k=2\pi/\lambda\); \(\omega=\) wave frequency; \(t=\) time after the wave crest passes through the original point; and \(d=\) water depth. The structural grid method is employed for numerical simulations, as shown in Figure 2. When simulating the high-amplitude wave, refinement of the grids around the water surface is required to achieve better computational accuracy. The grid sizes around the water surface in the wavelength and wave height directions are \(1/10a\) and \(1/25a\), respectively. There are 405 nodes on the airfoil, and the y+ value for the first mesh point away from the airfoil is approximately 1. The total number of mesh cells is 0.39 million. The numerical simulations are conducted by using the commercial software STAR-CCM + C. The unsteady flow is numerically simulated by solving the unsteady incompressible Reynolds averaged Navier-Stokes equations (RANS) with the \(k\)-\(\omega\) shear stress transport (SST) turbulence model. The transport equations of the model are  
\[\frac{\partial}{\partial t}(\rho k)+\
abla\cdot(\rho k\overline{\psi})=\
abla \cdot((\mu+\sigma_{k}\mu_{t})\
abla k)-\rho\beta*f_{\beta*}(\omega k-\omega_{0} k_{0})+G_{k}+G_{nl}+G_{b}+S_{k} \tag{4}\]  
\[\frac{\partial}{\partial t}(\rho\omega)+\
abla\cdot(\rho\omega\overline{\psi} )=\
abla\cdot((\mu+\sigma_{\omega}\mu_{t})\
abla k)-\rho\beta f_{\beta}\Big{(} \omega^{2}-\omega_{0}^{2}\Big{)}+G_{\omega}+D_{\omega}+S_{\omega} \tag{5}\]  
where \(f_{\beta^{*}}\) represents the free-shear modification factor, \(f_{\beta}\) represents the vortex-stretching modification factor, \(k_{0}\) and \(\omega_{0}\) are the ambient turbulence values that counteract turbulence decay, \(G_{k}\) represents the turbulent production, \(G_{nl}\) represents the \"non-linear\" production, \(G_{b}\) represents the buoyancy production, \(G_{\omega}\) represents the specific dissipation production, \(D_{\omega}\) represents the cross-diffusion term, \(S_{k}\) and \(S_{\omega}\) are source terms.