# Numerical Simulation of the Kelvin Wake Patterns
## 3 Numerical Methods
### Numerical Discretization
\begin{array}{ll}S^{(2)}(\rho^{\prime},\sigma^{\prime},x^{\prime},y^{\prime })&=\frac{1}{\sqrt{A^{\prime}(\rho^{\prime}-x^{\prime})^{2}+B^{\prime}(\rho^ {\prime}-x^{\prime})(\sigma^{\prime}-y^{\prime})+C^{\prime}(\sigma^{\prime}- y^{\prime})^{2}}}\\ &+\frac{1}{\sqrt{A^{\prime}(\rho^{\prime}-x^{\prime})^{2}-B^{\prime}(\rho^{ \prime}-x^{\prime})(\sigma^{\prime}+y^{\prime})+C^{\prime}(\sigma^{\prime}+y^ {\prime})^{2}}},\\ &A^{\prime}=1+{\zeta^{\prime}}^{2}_{x}(Q),\\ &B^{\prime}=2{\zeta^{\prime}}^{2}_{x}(Q)\zeta^{\prime}_{y}(Q),\\ &C^{\prime}=1+{\zeta^{\prime}}^{2}_{y}(Q).\end{array}  \tag{29}\]  
Now, the integral \(I^{\prime}(Q)\) contains the singularity and is computed exactly in terms of logarithms:  
\[ \begin{array}{ll}\int\frac{dsdt}{[As^{2}+Bst+C^{2}]}=&\frac{t}{A^{\frac{1}{2 }}}\ln\biggl{(}2 As+Bt+2\bigl{[}A\bigl{(}As^{2}+Bst+Ct^{2}\bigr{)}\bigr{]}^{ \frac{1}{2}}\biggr{)}\\ &+\frac{S}{C^{2}}\ln\biggl{(}2 Ct+Bs+2\bigl{[}C\bigl{(}As^{2}+Bst+Ct^{2}\bigr{)} \bigr{]}^{\frac{1}{2}}\biggr{)}+h_{1}(s)+h_{1}(t).\end{array}  \tag{30}\]  
Equation (9) is called the radiation condition; it can select the particular solution that does not contain the small spurious waves. By using the different radiation condition formulations, the results of different accuracies can be obtained [37]. Generally, plausible candidates for the radiation condition include requiring that far upstream:  
1. The free surface elevation is zero;
2. The vertical component of velocity is zero;
3. The free surface elevation decays exponentially;
4. The vertical component of velocity decays exponentially. By investigating these radiation condition formulations, the radiation conditions of class 4 are more successful in selecting the particular solution, which does not contain the small spurious waves. The radiation condition can be enforced near the upstream end of the computational domain. However, in the present problem, only horizontal components need to be solved. The vertical component of velocity is the function that consists of the horizontal components of velocity and elevation, as shown in Equation (6). Thus, the radiation condition formulations of the vertical component of velocity decaying exponentially should be reasonably altered. This paper enforces an equation in the form \(xf_{x}+nf=0\) along the boundary \(x^{\prime}=x_{1}^{\prime}\) for the horizontal components of velocity and free surface elevation decaying exponentially. There are \(4M\) equations by applying the radiation condition for \(l=1,\ \ldots,\ M\), as follows:  
\[x_{1}^{\prime}\Big{(}\big{(}\phi_{x}^{\prime}\big{)}_{1,l}-1\Big{)}+n\Big{(} \phi_{1,l}^{\prime}-x_{1}^{\prime}\Big{)}=0, \tag{31}\]  
\[x_{1}^{\prime}\big{(}\phi_{xx}^{\prime}\big{)}_{1,l}+n\Big{(}\big{(}\phi_{x}^{ \prime}\big{)}_{1,l}-1\Big{)}=0, \tag{32}\]  
\[x_{1}^{\prime}\big{(}\zeta_{x}^{\prime}\big{)}_{1,l}+n\zeta_{1,l}^{\prime}=0, \tag{33}\]  
\[x_{1}^{\prime}\big{(}\zeta_{xx}^{\prime}\big{)}_{1,l}+n\big{(}\zeta_{x}^{ \prime}\big{)}_{1,l}=0, \tag{34}\]  
where \(n\) is the decay coefficient. In conclusion, there are totally \(2(N+1)M\) equations to be solved for the \(2(N+1)M\) unknowns in the nonlinear ship wave problem. The boundary integral equation has been discretized to the \((N-1)M\) nonlinear equations, and additional \((N-1)M\) nonlinear equations are given by evaluating the free surface condition at the half mesh points, and \(4M\) equations are provided through applying the radiation condition, total \(2(N+1)M\) equations.