# Numerical Simulation of the Kelvin Wake Patterns
## 3 Numerical Methods  
### Numerical Discretization  
To solve the nonlinear problem numerically, the \(x-\)coordinates and \(y-\)coordinates of nodes in the \(N\times M\) mesh are \(x^{\prime}_{1},\ x^{\prime}_{2},\ \ldots,\ x^{\prime}_{N}\) and \(y^{\prime}_{1},\ y^{\prime}_{2},\ \ldots,\ y^{\prime}_{N}\) with intervals \(\Delta x^{\prime}\) and \(\Delta y^{\prime}\) in the \(x-\) and \(y-\)directions, respectively. The free surface elevation \(\zeta^{\prime}(x^{\prime},y^{\prime})\) and the velocity potential \(\phi^{\prime}(x^{\prime},y^{\prime})\) are represented by discrete point values \(\zeta^{\prime}_{k,l}\) and \(\phi^{\prime}_{k,l}\) at the nodes \((x^{\prime}_{k},y^{\prime}_{l})\), \(k=1,\ \ldots,\ N,\ l=1,\ \ldots,\ M\). Further, this paper chooses the \(x-\)derivatives of the functions \(\phi^{\prime}\) and \(\zeta^{\prime}\) as the basis for the solutions, together with the values of \(\phi^{\prime}\) and \(\zeta^{\prime}\) at the upstream boundary of the truncated domain, resulting in the vector of \(2(N+1)M\) unknowns \(\mathbf{u}\) to be  
\[ \begin{array}{ll}\mathbf{u}=&\left[\phi^{\prime\prime\prime}_{1,1},\left(\phi^{ \prime}_{x}\right)_{1,1},\ \ldots,\ \left(\phi^{\prime}_{x}\right)_{N,1},\ \ldots,\ \phi^{\prime}_{1,M},\ \left(\phi^{\prime}_{x}\right)_{1,M},\ \ldots,\ \left(\phi^{\prime}_{x}\right)_{N,M}, \right.\\ &\zeta^{\prime},\ (\zeta^{\prime}_{x})_{1,1},\ \ldots,\ (\zeta^{\prime}_{x})_{N,1},\ \ldots,\ \zeta^{ \prime}_{1,M},\ \left(\zeta^{\prime}_{x}\right)_{1,M},\ \ldots,\ (\zeta^{\prime}_{x})_{N,M} \right]^{\mathrm{T}}.\end{array}  \tag{24}\]  
Figure 3 shows the distribution of \(\mathbf{u}\) at an \(N\times M\) mesh, where \(\phi^{\prime}\) is located at the left side of the mesh, represented by these black dots; \(\phi^{\prime}_{x}\) is located at each node, represented by these congruent triangles. The distribution of \(\zeta^{\prime}\) and \(\zeta^{\prime}_{x}\). works in the same way. Figure 3: Illustration of the discrete unknowns at an \(N\times M\) mesh, the left black dots are the locations of \(\phi^{\prime}\) and \(\zeta^{\prime}\), the congruent triangles are the locations of the \(\phi^{\prime}_{x}\) and \(\zeta^{\prime}_{x}\). To solve this problem, the vector \(\mathbf{u}\) is firstly initialized. Then, the remaining \(\zeta^{\prime}\) can be obtained by trapezoidal\(-\)rule integration using \(\zeta^{\prime}_{x}\):  
\[\zeta^{\prime}_{k+1,l}=\zeta^{\prime}_{k,l}+\frac{1}{2}\Delta x^{\prime((\zeta^ {\prime}_{x})_{k}+(\zeta^{\prime}_{x})_{k+1,l})},\;l=1,\;\ldots,\;M,k=1,\; \ldots,\;N-1. \tag{25}\]  
Once \(\zeta^{\prime}\) is determined, it is facile to compute \(\zeta^{\prime}_{y}\) by fitting a cubic spline through the nodes \(\zeta^{\prime}_{k,1},\;\ldots,\;\zeta^{\prime}_{k,M}\) for \(k=1,\;\ldots,\;N\). With the same method, \(\phi^{\prime}\) and \(\phi^{\prime}_{y}\) at each node can be calculated by using \(\phi^{\prime}_{x}\)  
\[\phi^{\prime}_{k+1,l}=\phi^{\prime}_{k,l}+\frac{1}{2}\Delta x^{\prime((\phi^{ \prime}_{x})_{k}+(\phi^{\prime}_{x})_{k+1,l})},l=1,\;\ldots,\;M,k=1,\;\ldots, \;N-1. \tag{26}\]  
The boundary integral equation should be reformulated again, because there is a singularity in the second integral of Equation (14). It can be rewritten, as follows:  
\[\int_{y^{\prime}_{1}}^{y^{\prime}_{M}}\int_{x^{\prime}_{0}}^{x^{\prime}_{N}}[ \zeta^{\prime}_{\rho}(P)K^{(2)}\big{(}\rho^{\prime},\sigma^{\prime},x^{\prime },y^{\prime}\big{)}-\zeta^{\prime}_{\;x}(Q)\;S^{(2)}\big{(}\rho^{\prime}, \sigma^{\prime},x^{\prime},y^{\prime}\big{)}]d\sigma^{\prime}d\rho^{\prime}+ \zeta^{\prime}_{x}(Q)I^{\prime}(Q), \tag{27}\]  
and  
\[I^{\prime}(Q)=\int_{y^{\prime}_{1}}^{y^{\prime}_{M}}\int_{x^{\prime}_{0}}^{x^{ \prime}_{N}}\;S^{(2)}\big{(}\rho^{\prime},\sigma^{\prime},x^{\prime},y^{\prime }\big{)}d\sigma^{\prime}d\rho^{\prime}. \tag{28}\]  
where the kernel function \(S^{(2)}(\rho^{\prime},\sigma^{\prime},x^{\prime},y^{\prime})\) can be described, as follows:  
\[