Numerical Simulation of the Kelvin Wake Patterns

: The ship wave is of great interest for wave drag and coastal erosion. This paper proposes a mechanism of ship wave transformation to explore the effects of ship speed and ship size on the waveform. Firstly, based on the theory of potential ﬂow, the boundary integral equations for the Kelvin ship waves are obtained by deploying the different Kelvin sources or Rankine sources. Then, these integral equations are numerically discretized to a set of nonlinear equations. Finally, the Jacobian − free Newton–Krylov method with a preconditioner is adopted to solve the nonlinear equations. Though imitating plenty of different Kelvin wave patterns, the mechanism of ship wave transformation is proposed to conveniently generate the polymorphic Kelvin wake patterns. The above numerical simulation scheme is veriﬁed by comparing simulation results with real ship waves. After that, the wake angle is discussed with the effects of Froude number, source strength and source type by following the mechanism of ship wave transformation. The results show that the wake angle tends to decrease with ship speed but increase with ship size. In addition, for high ship speeds, the effect on the wake angle can be more dramatic.


Introduction
When a vessel or a body steadily translates into calm water, free surface waves are generated and these profiles are referred to as "Kelvin ship waves" [1]. The larger the Kelvin wake generated by a ship, the larger the drag force will be on the ship [2]. Meanwhile, the ship waves contain energy, which impacts the coastline and riverbank [3]. Calculating how much the wave drag and knowing how the energy distribution of a ship wave changes with respect to certain properties (e.g., ship speed, ship size and ship shape) can help design the ship hull and inform environmentally friendly shipping policies. Inferring characteristics of wake waveform can help identify ships that have been detected by monitoring instruments [4][5][6].
Initially, Kelvin [7] was fascinated by the V−shaped wake created by a ship advancing through water; one of his most remarkable results concerned the well−known constant wake angle 19.47 • , which is known in the community as the Kelvin angle. Lighthill [8] suggested a simplified stationary phase analysis of the dispersion relation, showing that the open dispersion curve is divided by an inflexion into two branches corresponding to transverse waves and divergent waves, respectively. Rabaud [9] noted that for sufficiently fast−moving ships, contrary to commonly held views, the wake angle that is observed behind a steadily moving ship is less than the well−known Kelvin angle, according to simulation results. This finding aroused the interest of many academics. Verberck [10] utilized a steadily moving pressure disturbance to analyze the variation in wake angle for large Froude numbers. Darmon [11] provided an explanation for this phenomenon, namely that the apparent wake angle is naturally provided by lining up the points of maximum wave amplitude along the divergent waves. Furthermore, this is similar for

Problem Formulation
Based on potential flow theory, this paper considers the irrotational flow of an inviscid, incompressible fluid of infinite depth and ignores the effects of surface tension. The fluid possesses an upper free boundary and is subject to the downward acceleration g of gravity. A Cartesian coordinate is defined such that the x− and y−axes in the plane of the free surface and z−axis point vertically with the origin at (0,0,0). The flow is directed along the positive x−axis with uniform speed U. Meanwhile, a source singularity whose strength is m or a Rankine source whose strength is κ is introduced at a distance L below the surface, as illustrated in Figure 1. The source strength represents ship size, a higher source strength means a larger ship size. The transient waves will be generated on the free surface, on account of the disturbance caused by a source. This is the steady−state problem, considering the wave pattern does not change over time. The free surface location is described by the equation z = ζ (x , y ), with labelling the velocity potential Φ (x , y , z ); here, the apostrophe means the dimensionless variable. Meanwhile, dimensionless variables can be defined by scaling all lengths with respect to L and all speeds with respect to U. Appl. Sci. 2022, 12, x FOR PEER REVIEW 3 of 18 free surface and −axis point vertically with the origin at (0,0,0). The flow is directed along the positive −axis with uniform speed . Meanwhile, a source singularity whose strength is or a Rankine source whose strength is is introduced at a distance below the surface, as illustrated in Figure 1. The source strength represents ship size, a higher source strength means a larger ship size. The transient waves will be generated on the free surface, on account of the disturbance caused by a source. This is the steady−state problem, considering the wave pattern does not change over time. The free surface location is described by the equation = ζ′( ′, ′) , with labelling the velocity potential Φ′( ′, ′, ′); here, the apostrophe means the dimensionless variable. Meanwhile, dimensionless variables can be defined by scaling all lengths with respect to and all speeds with respect to . Assuming the fluid is inviscid and incompressible, Laplace's equation can be satisfied, as follows: By deploying a source point, the velocity potential can be described, as follows: Φ ~ − 4 + + ( + 1) as ( , , ) → (0,0, −1).
As for the neighbourhood of the Rankine source point, the appropriate limiting behaviour can be described, as follows: Φ′~′ + ′ ′ 4 ( ′ + ′ + ( ′ + 1)) as ( ′, ′, ′) → (0,0, −1), where ′ and are the Kelvin source strength and Rankine source strength in the form: Since the kinematic and dynamic boundary conditions must be satisfied on the free surface, the nonlinear conditions are adopted, as follows: where ′ is the Froude number in the form: Assuming the fluid is inviscid and incompressible, Laplace's equation can be satisfied, as follows: By deploying a source point, the velocity potential can be described, as follows: As for the neighbourhood of the Rankine source point, the appropriate limiting behaviour can be described, as follows: 4π(x 2 + y 2 + (z + 1)) 3 2 as (x , y , z ) → (0, 0, −1), where and µ are the Kelvin source strength and Rankine source strength in the form: Since the kinematic and dynamic boundary conditions must be satisfied on the free surface, the nonlinear conditions are adopted, as follows: where F is the Froude number in the form: Finally, the flow will approach the free stream both far upstream and infinitely far below the free surface, the final two conditions are provided, as follows: The problem defined by Equations (1)-(10) can be solved by the boundary−integral method. The field region can be regarded as a half−sphere V whose boundaries are S ∞ and S T , except for a small sphere with surface S 1 about the source point (0, 0, −1) and a small hemisphere with surface S Q about typical point Q(x , y , z ) on the free surface, as shown in Figure 2. Thus, the boundary of the volume V can be written, as follows: Appl. Sci. 2022, 12, x FOR PEER REVIEW 4 of 18 Finally, the flow will approach the free stream both far upstream and infinitely far below the free surface, the final two conditions are provided, as follows: The problem defined by Equations (1)-(10) can be solved by the boundary−integral method. The field region can be regarded as a half−sphere whose boundaries are and , except for a small sphere with surface about the source point (0, 0, −1) and a small hemisphere with surface about typical point ( ′, ′, ′) on the free surface, as shown in Figure 2. Thus, the boundary of the volume can be written, as follows: Given random point ( ′, ′, ′) within volume , the distance between and is ′ , which can be described as According to Equation (9), the function Φ′ − ′ is harmonic within volume . As ( ′, ′, ′) is any point within volume , it can be proved that the function is harmonic. Therefore, The Green's second formula can be satisfied, as follows: where denotes the unit normal to the surface. As shown in Equation (11), the surface integral equation can be divided into four segments. Firstly, when surface expands to infinity, the result of this partial integration can be described as Secondly, the hemispherical surface shrinks to the surface point ( ′, ′, ′) , which means ′ tends to zero. Utilizing the mean value theorem of integrals, the result can be described as Given random point P(ρ , σ , τ ) within volume V, the distance between P and Q is R PQ , which can be described as According to Equation (9), the function Φ − x . is harmonic within volume V. As P(ρ , σ , τ ). is any point within volume V, it can be proved that the function 1 R PQ is harmonic. Therefore, The Green's second formula can be satisfied, as follows: where n denotes the unit normal to the surface. As shown in Equation (11), the surface integral equation can be divided into four segments. Firstly, when surface S ∞ expands to infinity, the result of this partial integration can be described as Secondly, the hemispherical surface S Q . shrinks to the surface point Q(x , y , z ), which means R PQ tends to zero. Utilizing the mean value theorem of integrals, the result can be described as Appl. Sci. 2022, 12, 6265 Thirdly, as the spherical surface S 1 collapses to the source point (0, 0, −1), the result of this partial integration can be described as Finally, Equation (13) yields where point P now lies on the punctures surface S T . The integral in Equation (17) can be rendered less singular by subtraction of the following term, which is equal to zero by utilizing Gauss flux theorem, Thus, Equation (17) can be re−written, as follows: where the symbol n denotes the downward pointing unit normal vector to the free surface and can be written, as follows: Then the surface integrals become the double integrals by projection onto the x − y plane using the following formula: On the free surface, the velocity potential can be regarded as a function simply of the two independent variables x and y and will be denoted by the symbol φ (x ,y ) = Φ (x , y , ζ (x , y )); the final boundary integral equation can be described, as follows: where the kernel functions K (1) (ρ , σ , x , y ) and K (2) (ρ , σ , x , y ) are described, as follows: Moreover, the free surface conditions can be simplified by the symbol φ (x , y ). Then, the kinematic and dynamic boundary conditions of the free surface are combined to be In conclusion, the boundary integral method is used to obtain the boundary integral Equation (22), which identically satisfies Laplace's equation, the limiting condition (2) and the far−field conditions (9) and (10). Meanwhile, the boundary integral equation should satisfy the combined free surface condition (23); thus, we are left to solve Equations (22) and (23), and the functions ζ and φ constitute the solution of the present problem.

Numerical Discretization
To solve the nonlinear problem numerically, the x−coordinates and y−coordinates of nodes in the N × M mesh are x 1 , x 2 , . . . , x N and y 1 , y 2 , . . . , y N with intervals ∆x and ∆y in the x− and y−directions, respectively. The free surface elevation ζ (x , y ) and the velocity potential φ (x , y ) are represented by discrete point values ζ k,l and φ k,l at the nodes x k , y l , k = 1, . . . , N, l = 1, . . . , M. Further, this paper chooses the x−derivatives of the functions φ and ζ as the basis for the solutions, together with the values of φ and ζ at the upstream boundary of the truncated domain, resulting in the vector of 2(N + 1)M unknowns u to be Figure 3 shows the distribution of u at an N × M mesh, where φ is located at the left side of the mesh, represented by these black dots; φ x is located at each node, represented by these congruent triangles. The distribution of ζ and ζ x . works in the same way. Moreover, the free surface conditions can be simplified by the symbol ′( , ). Then, the kinematic and dynamic boundary conditions of the free surface are combined to be In conclusion, the boundary integral method is used to obtain the boundary integral Equation (22), which identically satisfies Laplace's equation, the limiting condition (2) and the far−field conditions (9) and (10). Meanwhile, the boundary integral equation should satisfy the combined free surface condition (23); thus, we are left to solve Equations (22) and (23), and the functions ′ and ′ constitute the solution of the present problem.

Numerical Discretization
To solve the nonlinear problem numerically, the −coordinates and −coordinates  Figure 3 shows the distribution of at an × mesh, where ′ is located at the left side of the mesh, represented by these black dots; is located at each node, represented by these congruent triangles. The distribution of ′ and ′ works in the same way. To solve this problem, the vector u is firstly initialized. Then, the remaining ζ can be obtained by trapezoidal−rule integration using ζ x : Once ζ is determined, it is facile to compute ζ y by fitting a cubic spline through the nodes ζ k,1 , . . . , ζ k,M for k = 1, . . . , N. With the same method, φ and φ y at each node can be calculated by using φ x 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: (27) and where the kernel function S (2) (ρ , σ , x , y ) can be described, as follows: Now, the integral I (Q) contains the singularity and is computed exactly in terms of logarithms: 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 x f x + n f = 0 along the boundary x = x 1 for the horizontal components of velocity and free surface elevation decaying exponentially. There are 4M equations by applying the radiation condition for l = 1, . . . , M, as follows: 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.

Jacobian−Free Newton-Krylov Method
After numerical discretization, a nonlinear system of equations could be obtained, as follows: where u is the vector of unknowns of the length 2(N + 1)M. This nonlinear system can be solved by Newton's method [35]. The iterative process is as follows: set an initial guess u 0 , for t = 0, 1, 2, . . ., until convergence: where J(u t ) = F(u t ) is the system Jacobian [38]. However, Newton's method cannot solve large nonlinear systems efficiently. Therefore, based on Newton's method, this paper adopts the Jacobian−free Newton-Krylov method with a preconditioner, which is carefully constructed with entries taken from the Jacobian of the linearized problem. There is one more inner iteration, since Equation (36) is solved by using the iterative Generalized Minimum Residual algorithm [39]. The GMRES method is one of the Krylov subspace methods that are attractive as linear solvers in the context of nonlinear Newton iteration. Since each iteration does not require the exact value of δu t in solving nonlinear equations, the GMRES method can increase computation speed. Firstly, the approximate solution of δu t is found by projecting obliquely onto the Krylov subspace where m is the value of the subspace dimension and the accuracy of the solution increases with the subspace dimension, J t = J(u t ),F t = F(u t ). The matrix P ≈ J t is the preconditioner matrix, whose purpose is to construct an approximation to the Jacobian J t , which is cheap to form and to factorize. The calculation speed of the GMRES method can be significantly improved with a preconditioner matrix, because the spectrum of the preconditioned Jacobian J t P −1 exhibits a clustering of eigenvalues [40]. An initial linear residual r 0 is defined, given an initial guess δu 0 , for the Newton correction, Subsequently, the GMRES iteration minimizes r t 2 to a suitable value. These Jacobian− vector products can be approximated by applying first−order difference quotients: where v represents an arbitrary vector used in building the Krylov subspace [41], and the h is a small perturbation [40] In order to calculate Equation (36) conveniently, the preconditioner can be divided up into four equal submatrices and factorized using the block decomposition, where I is the unit matrix, A, B, C and D are the four submatrices, which are constructed in the base of Jacobian J t . It is possible to solve the system Pr = v by performing the following operations, According to the progressive order from Equations (43)-(45), the calculation of P −1 v in Equation (40) can be facilitated. The reasons are as follows: the submatrix A is tridiagonal, allowing for easy storage and fast factorization; the submatrices B and C are only used in matrix vector multiplication operations and, thus, can be implemented as functions that perform these operations rather than stored as matrices.

Simulation Results
According to the solution of the system of nonlinear equations, a perspective view of the nonlinear Kelvin wake pattern is given in Figure 4, for the flow past the Ranking source with F = 0.7 and µ = 0.3, computed on a 151 × 51 mesh with ∆x = 0.3, ∆y = 0.3. As shown in this free surface pattern, there are two different types of waves. One type of wave is the transverse wave, whose direction of propagation is consistent with the ship heading and the crests are slightly lower, being the main part of the whole pattern. The other waves are divergent waves, distributed symmetrically along the ship heading, whose crests appear to form ridges pointing diagonally away from the source, with the large steepness. The amplitude of the ship waves decays in proportion to the square root of distance along the x−axis, forming the well−known V−shaped Kelvin ship wave. In the numerical scheme proposed in this paper, there are three parameters to regulate the simulation results, with the polymorphic ship wave patterns appearing. The parameters include Froude number, source strength and source type. By adjusting these parameters, the wake angle, wave amplitude, prominent waves and wavelength can be changed. After plenty of experiments on the wake waveform, this paper proposed the following mechanism of ship wave transformation. The wake angle tends to decrease with Froude number but increase with the strength of the disturbance; the divergent waves will be prominent with the sufficiently large Froude number or the Rankine source; the wavelength dramatically increases with Froude number. It is suitable to choose the large Froude number, low source strength and Rankine source to simulate the ship wave with prominent divergent waves, small wake angle and long wavelength. In the opposite case the ship wave with prominent transverse waves, large wake angle and short wavelength will be simulated. Here, three kinds of real ship wave patterns are used as examples shown in Figure 5.  In the numerical scheme proposed in this paper, there are three parameters to regulate the simulation results, with the polymorphic ship wave patterns appearing. The parameters include Froude number, source strength and source type. By adjusting these parameters, the wake angle, wave amplitude, prominent waves and wavelength can be changed. After plenty of experiments on the wake waveform, this paper proposed the following mechanism of ship wave transformation. The wake angle tends to decrease with Froude number but increase with the strength of the disturbance; the divergent waves will be prominent with the sufficiently large Froude number or the Rankine source; the wavelength dramatically increases with Froude number. It is suitable to choose the large Froude number, low source strength and Rankine source to simulate the ship wave with prominent divergent waves, small wake angle and long wavelength. In the opposite case, the ship wave with prominent transverse waves, large wake angle and short wavelength will be simulated. Here, three kinds of real ship wave patterns are used as examples, shown in Figure 5.

1.
The ship waves have small wake angle, prominent divergent waves and long wavelength. In this case, the vessel's overall length is short but its speed is high. In Figure 5a, there is a wake wave from the speedboat. Therefore, the Rankine source is utilized to make divergent waves conspicuous, then continuously increasing the Froude number until the wake angle is suitable. The simulation result is consistent with the real ship wave pattern.

2.
The wake angles are the same as Kelvin angle, with prominent transverse waves and either the value of overall length or speed is not large, like the pilot boat in Figure 5c (note that: due to the far perspective, the entire wake pattern looks slightly smaller. However, the primary waveform is also consistent with Figure 5d). It is better to choose the Kelvin source, low source strength and the moderately small value of the Froude number. The perspective view of the pattern for F = 0.7 and = 0.5 is presented in Figure 5d. 3.
The ship waves have large wage angle, prominent divergent waves and short wavelength. Generally, only the large vessels can profile these wakes (Figure 5e). It is right to enhance the source strength to increase the wake angle, and then utilize the Rankine source to make the divergent waves prominent. The perspective view of the pattern for F = 0.7 and µ = 0.5 is presented in Figure 5f.
In conclusion, the simulation method adopted in this paper can simulate the polymorphic ship waves, and the simulation results are consistent with real ship waves. Therefore, it is feasible to discuss the wake angle on that basis. will be prominent with the sufficiently large Froude number or the Rankine source; the wavelength dramatically increases with Froude number. It is suitable to choose the large Froude number, low source strength and Rankine source to simulate the ship wave with prominent divergent waves, small wake angle and long wavelength. In the opposite case, the ship wave with prominent transverse waves, large wake angle and short wavelength will be simulated. Here, three kinds of real ship wave patterns are used as examples, shown in Figure 5.

Discussions of Kelvin Wake Angle
In Section 4, this paper proposes the approximate variation tendency of the waveform by adjusting the parameters. Then, in this section, the wake angle is discussed for Froude number, source strength and source type in more detail.

The Effect of Froude Number on Wake Angle
To demonstrate the effect of Froude number on apparent wake angle, Figure 6 shows three free surface patterns for flow past the Rankine source, computed for F = 0.7, 2.5, 5.0 and µ = 0.4 on a 181 × 61 mesh with ∆x = 0.3, ∆y = 0.3. The wake angle decreases as Froude number increases. The reason is that the highest peaks do not locate in the outermost divergent under enough large Froude numbers, and then the wake angle gets smaller apparently. In Figure 6a, the Froude number is a moderately small value and the wake angle is almost the same as Kelvin angle. However, in Figure 6b,c, the wake angle is much smaller than the Kelvin angle, ending up with the wake angle and Froude approximately inversely proportional. Appl Therefore, we should use specific data to obtain an accurate conclusion. The wake angle is measured in the contour plot, whose highest peaks are easily found. As shown in Figure 7, the angle between the line connecting these peaks and the ship heading is the wake angle . To clearly illustrate the effect of Froude number ′, this paper keeps the source strength and source type invariant, obtains a series of data of wake angles with changing Froude number and sets two groups of data for comparing. The Cartesian coordinates of ′ and are established, as shown in Figure 8. The Froude number ′ ranges from 0.8 to 5.0, the trend of the wake angle obviously tends to decrease with Froude number. However, gradients are obviously different between the range of Froude number from 0.8 to 2.0 and the range of Froude number from 3.0 to 5.0. In the first interval, the wake angle decreases slowly with Froude number and the value of the wake angle fluctuates up and down; in the second interval, wake angle decreases dramatically with Froude number. In sum, the downtrend of the wake angle with Froude number has two stages; it is only when the Froude number reaches a threshold that the wake angle is inversely proportional to the Froude number; the threshold is independent of the source strength, around 2.0 when the Rankine source is used as a disturbance. This finding can be reflected in real situations; only when the ship speed is high enough will the transformation of wake angle be noticeable.  Therefore, we should use specific data to obtain an accurate conclusion. The wake angle θ is measured in the contour plot, whose highest peaks are easily found. As shown in Figure 7, the angle between the line connecting these peaks and the ship heading is the wake angle θ. To clearly illustrate the effect of Froude number F , this paper keeps the source strength and source type invariant, obtains a series of data of wake angles with changing Froude number and sets two groups of data for comparing. The Cartesian coordinates of F and θ. are established, as shown in Figure 8. The Froude number F ranges from 0.8 to 5.0, the trend of the wake angle θ obviously tends to decrease with Froude number. However, gradients are obviously different between the range of Froude number from 0.8 to 2.0 and the range of Froude number from 3.0 to 5.0. In the first interval, the wake angle decreases slowly with Froude number and the value of the wake angle fluctuates up and down; in the second interval, wake angle decreases dramatically with Froude number. In sum, the downtrend of the wake angle with Froude number has two stages; it is only when the Froude number reaches a threshold that the wake angle is inversely proportional to the Froude number; the threshold is independent of the source strength, around 2.0 when the Rankine source is used as a disturbance. This finding can be reflected in real situations; only when the ship speed is high enough will the transformation of wake angle be noticeable. Therefore, we should use specific data to obtain an accurate conclusion. The wake angle is measured in the contour plot, whose highest peaks are easily found. As shown in Figure 7, the angle between the line connecting these peaks and the ship heading is the wake angle . To clearly illustrate the effect of Froude number ′, this paper keeps the source strength and source type invariant, obtains a series of data of wake angles with changing Froude number and sets two groups of data for comparing. The Cartesian coordinates of ′ and are established, as shown in Figure 8. The Froude number ′ ranges from 0.8 to 5.0, the trend of the wake angle obviously tends to decrease with Froude number. However, gradients are obviously different between the range of Froude number from 0.8 to 2.0 and the range of Froude number from 3.0 to 5.0. In the first interval, the wake angle decreases slowly with Froude number and the value of the wake angle fluctuates up and down; in the second interval, wake angle decreases dramatically with Froude number. In sum, the downtrend of the wake angle with Froude number has two stages; it is only when the Froude number reaches a threshold that the wake angle is inversely proportional to the Froude number; the threshold is independent of the source strength, around 2.0 when the Rankine source is used as a disturbance. This finding can be reflected in real situations; only when the ship speed is high enough will the transformation of wake angle be noticeable.

The Effect of Source Strength on Wake Angle
In this paper, different source strengths = 1.0, 1.5, 2.5 are utilized to generate different wave patterns for the flow past Kelvin source with ′ = 0.8, as shown in Figure 9.
It is well known that the Kelvin angle is 19.47°. Comparing the wake angles with the Kelvin angle, the wake angle is slightly smaller than the Kelvin angle in Figure 9a. However, the Kelvin angle in Figure 9c is larger than the Kelvin angle. Although the wake angle change is not significant, the wake angle is tending to slowly increase with the source strength. In other words, the wake angle is tending to slowly increase with the ship size. Meanwhile, Cartesian coordinates of and are established by keeping Froude number and source type invariant and adjusting source strength, obtaining two sets of data with ′ = 0.8 and ′ = 1.2, shown in Figure 10. The source strength ranges from 1.0 to 3.0 and the wake angle increases with the source strength. However, the two sets of data show slightly different trends. With ′ = 0.8, the wake angle increases at an almost steady rate with source strength; with ′ = 1.2, the wake angle barely changes as source strength increases at first, then increases at a steady rate after source strength reaching 1.8. In a word, compared with Froude number, although the source strength cannot greatly affect the wake angle, it is clear that the wake angle increases with source strength.

The Effect of Source Strength on Wake Angle
In this paper, different source strengths = 1.0, 1.5, 2.5 are utilized to generate different wave patterns for the flow past Kelvin source with F = 0.8, as shown in Figure 9. It is well known that the Kelvin angle is 19.47 • . Comparing the wake angles with the Kelvin angle, the wake angle is slightly smaller than the Kelvin angle in Figure 9a. However, the Kelvin angle in Figure 9c is larger than the Kelvin angle. Although the wake angle change is not significant, the wake angle is tending to slowly increase with the source strength. In other words, the wake angle is tending to slowly increase with the ship size.

The Effect of Source Strength on Wake Angle
In this paper, different source strengths = 1.0, 1.5, 2.5 are utilized to generate different wave patterns for the flow past Kelvin source with ′ = 0.8, as shown in Figure 9.
It is well known that the Kelvin angle is 19.47°. Comparing the wake angles with the Kelvin angle, the wake angle is slightly smaller than the Kelvin angle in Figure 9a. However, the Kelvin angle in Figure 9c is larger than the Kelvin angle. Although the wake angle change is not significant, the wake angle is tending to slowly increase with the source strength. In other words, the wake angle is tending to slowly increase with the ship size. Meanwhile, Cartesian coordinates of and are established by keeping Froude number and source type invariant and adjusting source strength, obtaining two sets of data with ′ = 0.8 and ′ = 1.2, shown in Figure 10. The source strength ranges from 1.0 to 3.0 and the wake angle increases with the source strength. However, the two sets of data show slightly different trends. With ′ = 0.8, the wake angle increases at an almost steady rate with source strength; with ′ = 1.2, the wake angle barely changes as source strength increases at first, then increases at a steady rate after source strength reaching 1.8. In a word, compared with Froude number, although the source strength cannot greatly affect the wake angle, it is clear that the wake angle increases with source strength.  Meanwhile, Cartesian coordinates of θ and are established by keeping Froude number and source type invariant and adjusting source strength, obtaining two sets of data with F = 0.8. and F = 1.2, shown in Figure 10. The source strength . ranges from 1.0 to 3.0 and the wake angle increases with the source strength. However, the two sets of data show slightly different trends. With F = 0.8, the wake angle increases at an almost steady rate with source strength; with F = 1.2, the wake angle barely changes as source strength increases at first, then increases at a steady rate after source strength reaching 1.8. In a word, compared with Froude number, although the source strength cannot greatly affect the wake angle, it is clear that the wake angle increases with source strength.

The Effect of Source Type on Wake Angle
The two free surface patterns for flow past the Kelvin source and Rankine source, respectively, are computed for the Froude number ′ = 0.8 . The amplitudes of ship waves, generated by two source types, should be adjusted to equalize. A wave profile computed for Kelvin source with ′ = 0.9 is present in Figure 11a and the wave profile computed for the Rankine source with ′ = 0.5 is given in Figure 11b. After measurement, the wake angles in the two free surface patterns are approximately equal, which means that source type does not affect the variation in wake angle. Furthermore, the wake angle has the same trend, no matter whether there is a Kelvin source or Rankine source as a disturbance. In Figure 12, there are trends of wake angle with the effect of Froude number for Kelvin source and Rankine source; both downtrends have two stages, but the first interval has a longer range in Figure 12a. The wake angle sharply decreases at Froude number 2.8 when the flow passes the Kelvin source. However, the wake angle sharply decreases at Froude number 2.0 when the flow passes the Rankine source. Therefore, the source type has an effect on the threshold about Froude number.
On the other hand, the Rankine source possesses a better capacity to avoid introducing very small spurious waves throughout the domain. Focusing on the upstream domain in Figure 11a, the spurious waves can be found, so that the whole free surface elevation is touched by spurious waves. Nevertheless, if we want to simulate the ship wave accompanied with sea clutters, introducing small spurious waves is a reasonable method.

The Effect of Source Type on Wake Angle
The two free surface patterns for flow past the Kelvin source and Rankine source, respectively, are computed for the Froude number F = 0.8. The amplitudes of ship waves, generated by two source types, should be adjusted to equalize. A wave profile computed for Kelvin source with = 0.9 is present in Figure 11a and the wave profile computed for the Rankine source with µ = 0.5 is given in Figure 11b. After measurement, the wake angles in the two free surface patterns are approximately equal, which means that source type does not affect the variation in wake angle. Furthermore, the wake angle has the same trend, no matter whether there is a Kelvin source or Rankine source as a disturbance. In Figure 12, there are trends of wake angle with the effect of Froude number for Kelvin source and Rankine source; both downtrends have two stages, but the first interval has a longer range in Figure 12a. The wake angle sharply decreases at Froude number 2.8 when the flow passes the Kelvin source. However, the wake angle sharply decreases at Froude number 2.0 when the flow passes the Rankine source. Therefore, the source type has an effect on the threshold about Froude number. Appl. Sci. 2022, 12, x FOR PEER REVIEW 14 of 18 Figure 10. A plot of apparent wake angle against Froude number for a flow past the Kelvin source, with the source strength ′ = 0.8 (the red circles) and ′ = 1.2 (the blue circles).

The Effect of Source Type on Wake Angle
The two free surface patterns for flow past the Kelvin source and Rankine source, respectively, are computed for the Froude number ′ = 0.8 . The amplitudes of ship waves, generated by two source types, should be adjusted to equalize. A wave profile computed for Kelvin source with ′ = 0.9 is present in Figure 11a and the wave profile computed for the Rankine source with ′ = 0.5 is given in Figure 11b. After measurement, the wake angles in the two free surface patterns are approximately equal, which means that source type does not affect the variation in wake angle. Furthermore, the wake angle has the same trend, no matter whether there is a Kelvin source or Rankine source as a disturbance. In Figure 12, there are trends of wake angle with the effect of Froude number for Kelvin source and Rankine source; both downtrends have two stages, but the first interval has a longer range in Figure 12a. The wake angle sharply decreases at Froude number 2.8 when the flow passes the Kelvin source. However, the wake angle sharply decreases at Froude number 2.0 when the flow passes the Rankine source. Therefore, the source type has an effect on the threshold about Froude number.
On the other hand, the Rankine source possesses a better capacity to avoid introducing very small spurious waves throughout the domain. Focusing on the upstream domain in Figure 11a, the spurious waves can be found, so that the whole free surface elevation is touched by spurious waves. Nevertheless, if we want to simulate the ship wave accompanied with sea clutters, introducing small spurious waves is a reasonable method.

Conclusions
In order to elucidate the effect of ship size on the Kelvin ship waves, this paper implements the numerical scheme, which combines the boundary integral method with the Jacobian-free Newton-Krylov method, for the simulation of nonlinear ship waves on an infinite fluid. Though imitating plenty of different Kelvin wave patterns, the mechanism of ship wave transformation is proposed, and then it is verified by comparing simulation results with real ship waves. Finally, the wake angle is discussed with the effects of Froude number, source strength and source type and the following conclusions can be drawn.
1. The wake angle tends to decrease with Froude number and the downtrend has two stages. The wake angle is inversely proportional to Froude number and decreases dramatically after Froude number reaches a threshold, which is around 2.0 for Rankine source and is around 2.8 for Kelvin source. 2. The wake angle tends to increase with source strength, meaning that ship size can affect the ship waves; the larger the ship size, the larger the wake angle generated will be. 3. Because the gradient of source strength is lower than the gradient of Froude number, the wake angle change caused by ship size is not as visually obvious as the ship speed. Meanwhile, it is hard for source strength to increase the wake angle when the Froude number is large. 4. With either the Kelvin source or Rankine source as a disturbance, the variation trends of wake angle are approximately identical with the effects of Froude number and source strength. However, a more accurate solution for nonlinear ship waves can be solved when the free surface flow passes Rankine source.
This paper suggests that there will be a sudden change in the downtrend of wake angle when the ship speed reaches a threshold. The next stage in research will pay attention to exploring the cause of the sudden change and investigating the influencing factors of threshold, e.g., water depth, ship shape, etc. The theoretical conclusion about threshold value will be verified by model test. This study will provide reference for informing environmentally friendly shipping policies, e.g., speed limits, etc. On the other hand, the Rankine source possesses a better capacity to avoid introducing very small spurious waves throughout the domain. Focusing on the upstream domain in Figure 11a, the spurious waves can be found, so that the whole free surface elevation is touched by spurious waves. Nevertheless, if we want to simulate the ship wave accompanied with sea clutters, introducing small spurious waves is a reasonable method.

Conclusions
In order to elucidate the effect of ship size on the Kelvin ship waves, this paper implements the numerical scheme, which combines the boundary integral method with the Jacobian-free Newton-Krylov method, for the simulation of nonlinear ship waves on an infinite fluid. Though imitating plenty of different Kelvin wave patterns, the mechanism of ship wave transformation is proposed, and then it is verified by comparing simulation results with real ship waves. Finally, the wake angle is discussed with the effects of Froude number, source strength and source type and the following conclusions can be drawn.

1.
The wake angle tends to decrease with Froude number and the downtrend has two stages. The wake angle is inversely proportional to Froude number and decreases dramatically after Froude number reaches a threshold, which is around 2.0 for Rankine source and is around 2.8 for Kelvin source.

2.
The wake angle tends to increase with source strength, meaning that ship size can affect the ship waves; the larger the ship size, the larger the wake angle generated will be. 3.
Because the gradient of source strength is lower than the gradient of Froude number, the wake angle change caused by ship size is not as visually obvious as the ship speed. Meanwhile, it is hard for source strength to increase the wake angle when the Froude number is large.

4.
With either the Kelvin source or Rankine source as a disturbance, the variation trends of wake angle are approximately identical with the effects of Froude number and source strength. However, a more accurate solution for nonlinear ship waves can be solved when the free surface flow passes Rankine source.
This paper suggests that there will be a sudden change in the downtrend of wake angle when the ship speed reaches a threshold. The next stage in research will pay attention to exploring the cause of the sudden change and investigating the influencing factors of threshold, e.g., water depth, ship shape, etc. The theoretical conclusion about threshold value will be verified by model test. This study will provide reference for informing environmentally friendly shipping policies, e.g., speed limits, etc.