Since Chappelear (1961) developed an early numerical method, many other numerical methods have been derived (Dean, 1965; Rienecker and Fenton, 1981; Chaplin, 1979; Cokelet, 1977; Fenton, 1988). Rienecker and Fenton (1981) solved the coefficients and some wave properties directly using Newton’s method. This method was further simplified by Fenton (1988). The major modification was that all the required partial differentials were calculated numerically (Tao et al., 2007). Although neither of the above numerical methods could be applied to waves in deep water (Fenton, 1990), Shin (2023) calculated deep water breaking waves with an error of less than 2.04×10⁻³ percent, demonstrating that Fourier approximation is effective even for deep water waves. The objective of this study was to extend Shin’s method (Shin, 2023) to include rotational water waves with a shear current in a finite depth (Shin, 2022). The basic concept of Shin’s method (Shin, 2023) for establishing a set of equations was similar to Fendon’s method (Rienecker and Fenton, 1981; Fenton, 1988), and some issues with Fenton’s method have been addressed. The first issue is inherent in Newton’s method, which necessitates an initial solution in close proximity to the final solution. In common with other versions of the Fourier approximation method (Chappelear, 1961; Dean, 1965; Chaplin, 1979), it is sometimes necessary to solve a sequence of lower waves, extrapolating forward in height steps until the desired height is reached (Fenton, 1988). These methods can occasionally converge to the wrong solution with very long waves (Fenton, 1988). This paper shows that regression analysis can easily avoid this issue (Shin, 2019).

The second issue arises from the coordinate system. Because Fenton’s method adopts the moving coordinate system proposed by Dean (1965), the abscissa is the position from the crest in the range, −L≤x≤L. Therefore, the relative position with respect to the wavelength, L, is unknown because the wavelength is unknown. As a result, the number of unknown variables to be determined in Fenton’s method is more than double the number in this study. Furthermore, a nonlinear equation for calculating the wavelength was coupled to the set of equations for calculating the coefficients. This coupling complicates the problem because Newton’s method requires partial derivatives with respect to the wavelength. The issue can be avoided easily using the dimensionless coordinate system in Fig. 1 because it is independent of the wavelength. The abscissa is the phase in the range −π≤β≤π, and is known. After determining the steepness, S, the wave-number is calculated as k=S/H because wave height, H is known. The steepness is calculated using the wave height condition, which is integrated automatically into the set of equations because it is linear. Therefore, it is not necessary to calculate the partial derivatives with regard to the wavelength in this study.

The third issue arises from the water depth condition used to calculate the reference depth, D, which is unknown. Fenton’s method requires partial derivatives with respect to the reference depth because the water depth condition is coupled to the set of nonlinear equations used to calculate the coefficients. Most errors in Fenton’s method occur when the water depth condition is calculated numerically. The required series order must be increased to reduce the error, and Simpson’s rule was used for the numerical integration. This issue was avoided by calculating the reference depth independently using the shooting method in the study. In numerical analysis, the shooting method is used to solve a boundary value problem by transforming it into an initial value problem. A set of “n” equations with “n” unknowns can be converted to a set of “n−1” equations with “n−1” unknowns because one variable can be calculated independently, assuming it is a known value in advance. The assumed value was determined using the shooting method. Although it has the disadvantage of requiring the calculation of one procedure as two procedures, the numerical formulation simplifies the process and significantly reduces the total number of steps needed to converge to a complete solution. The reference depth is not coupled to the set of equations for calculating the Fourier coefficients because the water depth condition is calculated independently.

The fourth issue is a method for calculating the partial derivatives required in Newton’s method. While all the necessary partial derivatives were obtained numerically in Fenton (1988), they were calculated analytically by tensor analysis in the present study. Therefore, there are no errors in the partial derivatives in this study.

Although the above-mentioned studies are irrational flows, this study is rotational flow based on Shin (2022). The numerical method used to calculate the Fourier coefficients and the related wave properties was not presented by Shin (2022) but is presented in the present study. Shin (2018) presented the numerical method for irrotational waves. The method was extended to rotational waves in this study. Shin (2022) reported the complete solutions to the Euler equations and calculated the waves by Le Mehaute et al. (1968). The fluid velocities calculated by irrotational wave theories (Fenton, 1990; Tao et al., 2007) and by rotational wave theory (Shin, 2022) were compared with the experimental data reported by Le Mehaute et al. (1968) on the same graphs. Unlike the irrotational waves, the results of the rotational wave theory agreed well with the experimental data. The calculated wave-breaking limit was in good accordance with the Miche formula. While Rienecker and Fenton (1981), Fenton (1988), Tao et al. (2007), Cokelet (1977), and Vadden-Broeck and Schwarts (1979) calculated the dispersion relations for some waves in a relative depth of 0.7, this study calculated these relations for many waves with different heights at various depths. The results were then extrapolated to all waves using regression analysis (Shin, 2019). The dispersion relations for a relative depth of 0.7 were compared and showed good accordance with the other results. The convergence of this study was tested according to the series orders. Even for relatively small orders, the error was significantly lower than Fenton’s method and Tao et al. (2007).

This chapter summarizes the solution reported by Shin (2022). The two coordinate systems were adopted. The first was a conventional coordinate system (t, x, y), as reported by Dean et al. (1984), in which the origin is located at a fixed point on the still water line (SWL). The x-axis is in the direction of wave propagation; the y-axis points vertically upward, and t is the time. The fluid domain is bounded by a flat bed, y = −h (water depth), and by a free surface y = η (t, x). The second is a dimensionless coordinate system, (β α), as shown in Fig. 1, which is a stationary frame. The origin is located under the crest of the bed. Therefore, the wave profile is fixed with a periodic, even function in the system. The abscissa is phase β = kx − ωt in −π ≤ β ≤ π, where k = 2π/L is the wave number and ω = 2π/T is the angular frequency. Here, T is the wave period; L is the wavelength, and the ordinate is the dimensionless elevation from the bed, α = k(y + h) in 0 ≤ α ≤ γ, where γ = k(η + h) is the dimensionless free surface elevation from the bed. The term ζ = γ − D is a dimensionless free surface elevation measured from the reference line (the horizontal line passing through two points on the free surface at two phases β = ±π/2. The reference depth D = k(η_o + h) is the dimensionless distance from the bed to the reference line where η_o = η(±π/2). Because the coordinate system is independent of the wavelength, it has several advantages, as presented by Shin (2023). The dimensionless quantity of a flow field f is denoted by f̄ except for those defined separately, such as S and θ. Using the notation, the stream function is defined by ψ = ψ̄ω/k².The horizontal and vertical velocities are defined by u = Cū and v = Cv̄, respectively, where C = ω/k is the celerity. The vorticity is defined by Ω = ωΩ̄. The pressure is defined by p = ρC²p̄, where ρ is the water density. The dimensionless wave height (steepness) is defined by S = kH. The linear steepness θ = ω²H/g (where g is the gravity) is a constant for a particular wave. When h → ∞ and H → 0, S → θ, which gives the dispersion relation (ω² = gk) of deep water linear waves. The stream function to satisfy the incompressible Euler equations is as follows: where 〈n〉def̳n2+ε; N is the series order; and c_n are the Fourier coefficients. The first and second terms describe a current. On the other hand, the two terms can be deleted for an irrotational wave (ε = 0). The vorticity is presented as follows:

This solution also contains irrotational waves, where the constant, ε is referred to as the strength of vorticity because Ω̄ = 0 when ε = 0. When ε is non-positive, the vorticity has the same direction as the motion of the particles. The iso-stream and iso-vorticity lines can be calculated using Newton’s method in the Appendix. The free surface is an iso-vorticity line given by Ω̄(β,γ) = C₁ε, where C₁ is a constant. Substituting Ω̄(β,γ) = C₁ε into Eq. (2), the wave profile is determined to satisfy the KFSBC (Kinematic free surface boundary condition) as follows:

The constant C₁ is determined from the definition of the reference line in Fig. 1; that is, γ(±π/2) = D. The profile is an implicit function that can be calculated using Newton’s method in the Appendix. When ε → 0, Eq. (3) yields the profile defined by Shin (2018). Moreover, the sea bed is an iso-vorticity line given by Ω̄(β,0) = 0. Differentiating Eq. (1) with respect to α and β, the horizontal velocity becomes the following.

The vertical velocity is as follows:

Substituting Eqs. (1)–(5) to the Euler momentum equations, Bernoulli’s principle for rotational flow can be calculated as follows. where Ū² = V̄² − ε (α − ψ̄)², and Q is a constant. A mechanical analog to irrotational and rotational flows was depicted by considering a carnival Ferris wheel in Fig. 2.11 in Dean et al. (1984). While the irrotational motion of chairs involves rotating around the center of the Ferris wheel, the rotational motion of the chairs involves both revolving around the center of the Ferris wheel and rotating around their center. Therefore, the total kinetic energy of the rotational motion is represented by the sum of the revolving kinetic energy around the center of the Ferris wheel and the rotating kinetic energy around their center. Analogous to the Ferris wheel, the total kinetic energy of rotational flow is represented by Ū² = V̄² − ε (α − ψ̄)², where the first term corresponds to the revolving kinetic energy around the center of the Ferris wheel, and the second term corresponds to the rotating kinetic energy around its center. When ε = 0, Ū² = V̄² = (1 − ū)² + v̄², which is the kinetic energy for irrotational flow. In contrast to Baddour et al. (1998), ε is not positive because kinetic energy functions should be positive-definite functions. Therefore, let εdef̳−σ2, where σ is a non-negative real number. Hence, 〈0〉 = σi. Substituting Eqs. (1), (4), and (5) into the relation, Ū² = V̄² − ε(α − ψ̄)², the dynamic pressure (kinetic energy) Ū² is as follows: where and where Z₁ = (〈n〉〈m〉 + nm − ε)/4W; Z₂ = (〈n〉〈m〉 − nm + ε)/4W; Z₃ = (〈n〉〈m〉 − nm − ε)/4W; Z4=〈n〉〈m〉+nm+ε4W; and W = cosh〈n〉D cosh〈m〉D. If Λ_0n is the component of column vector {Λ_0n} and c_n is the component of the column vector {c_n}, the second term on the right side in Eq. (7) becomes an inner product of the two vectors. Accordingly, Λ_nm can be considered a component of the square matrix [Λ_nm] in N-dimensional space. Therefore, the third term on the right side of Eq. (7) is presented as {c_n}^t[Λ_nm]{c_m}, where {c_n}^t is a row vector. This concept helps formulate numerical calculations in this study. Considering the average speed Ū_b over a wave period on the sea bed, the constant c₀ is determined as follows:

For the average speed Ū_s over a wave period on still water, the strength of vorticity is as follows:

From Stokes’ breaking criterion, |Ū_b| ≤ 1 and |Ū_s| ≤ 1. The Bernoulli’s constant Q is determined from the definition of the reference line (i.e., γ(π/2) = D) and the DFSBC (Dynamic free surface boundary condition), i.e., p̄(π/2, D) = 0 because the point (π/2, D) is on the free surface. Therefore, the pressure field is

By applying the DFSBC, i.e., p̄(β, γ) = 0 on α = γ to Eq. (12), the other wave profile can be determined as follows:

The wave height condition is defined as

The water depth condition is defined as

The solution contains N+2 unknown constants: N Fourier coefficients, c_n, steepness, S, and reference depth, D. Eqs. (3) and (13) present two wave profiles with implicit functions. The two wave profiles are even functions considering the coordinate system. The Fourier series of a periodic even function is presented as f(x)=∑n=0∞cncosnx, in −π ≤ x ≤ π. The coefficients are generally determined from the orthogonality of trigonometric functions. The method is inapplicable when f(x) is an implicit function, as expressed in Eqs. (3) and (13). The problem can be solved by converting the series to a set of algebraic equations, which are obtained by calculating the series at some phases instead of all phases. This is done by replacing the infinite series with a truncated series, i.e., ∑n=0Ncncosnxm=f(xm) for m = 1, …, N. There are N algebraic equations for calculating the coefficients, c_n, because cosnx_m and f(x_m) are known. Hence, there are a set of N equations, and Eqs. (3) and (13) are equal to each other at N+1 phases (Note that Eqs. (3) and (13) are already equal to each other at phase ±π/2 because the integral constants were determined to satisfy the definition of the reference line) and two equations, i.e., Eqs. (14) and (15). Therefore, there are N+2 equations to determine the unknown constants.

Referring to Eqs (3)–(5), the denominators are expressed as transcendental functions of the reference depth. When determining the reference depth using Newton’s method in conjunction with the coefficients, it is necessary to calculate the partial derivative with respect to the reference depth. Furthermore, Eq. (15) only permits a numerical approach. This difficulty can be overcome by calculating the reference depth independently using the shooting method, assuming that the reference depth is known. Based on the assumption, Eq. (15) is independent of the other equations and will be calculated in the next chapter. In addition to this assumption, N parameters, X_m, are introduced in this chapter to simplify the numerical formulation. Unlike Fenton’s method, X_m are merely parameters for calculating the coefficients and the steepness in this chapter. After determining the coefficients, X_m are no longer used. Hence, the wave elevation is denoted with X_m instead of ζ_m in this chapter.

When the wave profile ζ and the reference depth D are prescribed, it is possible to convert Eq. (3) into a set of linear equations for the coefficients. Because the wave profile is an even function, the phases β_m(m = 1,2, …, N) are considered in the range 0 ≤ β_m ≤ π and β_m ≠ π/2 because the two wave profiles are already equal at phase β = π/2. In addition, β₁ = 0 and β_N = π. Denoting ζ(β_m) as X_m, i.e., X_m = ζ(β_m), Eq. (3) is presented.

The summation convention is considered in Eq. (16). The repeated subscript “n” is a dummy subscript, whereas K_mn is a second-order tensor, and the terms c_n and X_m are vectors in N-dimensional space. From Eq. (3), the component of the tensor K_mn is presented as where X̄_m = X_m + D. Because 〈0〉= σi, the component of the first-order tensor, Y_m is presented as

The summation convention is not applied when the component of a tensor is presented as in Eqs. (17) and (18). Using the inverse tensor G_mn of the tensor K_mn, the solution to Eq. (16) is determined easily as follows.

Eq. (13) is also evaluated at the same phase. A set of N nonlinear equations for calculating the N parameters X_n is obtained by substituting Eq. (19) into Eq. (13) at the same phase because the wave height condition in Eq. (14) is expressed as S = X₁ − X_N in this approach. For Newton’s method, the error vector E_m is defined from Eq. (13) as follows:

From Eqs. (8) and (9), the second-order tensor A_pn and the third-order tensor B_pnm are defined as

Substituting Eq. (19) into Eq. (20)E_p = 0 gives a set of nonlinear equations for calculating the parameters, X_q. Therefore, unlike Fenton’s method (Rienecker and Fenton, 1981; Fenton, 1988), the set of N nonlinear equations is solved using Newton’s method in this study. Therefore, when using the same series order, the number of unknown constants is less than half of that in Fenton’s method.

Denoting partial derivatives of a tensor by making use of commas and indices as ∂()∂Xq=(),q and expanding Eq. (20) as a Taylor series about X_q points in terms of ΔX_q and ignoring higher-order terms, the following set of linear equations is obtained:

Because ΔXq(r)=Xq(r+1)−Xq(r), the solution in the next step is

The superscript (r) means the step of Newton’s method. It is clear that all the steps are r^th in all equations except Eq. (24). Thus, for the simplification of equations, the superscript (r) was omitted in all the equations except Eq. (24). Differentiating Eq. (20) with respect to X_q, the partial derivative E_p,q of the error vector is

Note that c_n,qc_m = c_nc_m,q Because S = X₁ −X_N,

Differentiating Eq. (16) with respect to X_p results in

Multiplying Eq. (27) by the tensor G_im, the partial derivative of the coefficient is determined as follows. where G_imK_mn = δ_in and δ_inc_n,p = c_i,p and δ_in is the second-order isotropic tensor. Differentiating Eq. (17) with respect to X_p, the components of the third-order tensor K_mn,p are presented as

Differentiating Eq. (18) with respect to X_p, the component of the second-order tensor Y_m,p is presented as

Differentiating Eq. (21) with respect to X_q, the component of the third-order tensor A_pn,q is expressed as follows:

Differentiating Eq. (22) with respect to X_q, the component of the forth-order tensor B_pnm,q is where Z1*=Z1(〈n〉+〈m〉); Z2*=Z2(〈n〉−〈m〉); Z3*=Z3(〈n〉+〈m〉); and Z4*=Z4(〈n〉−〈m〉).

Based on the assumption that the reference depth is known, Eq. (16) is linear and decoupled from Eq. (20), unlike in Fenton’s method (Rienecker and Fenton, 1981; Fenton, 1988). Dividing the inner product of the error vector, E_m by N, the mean squared error (MSE) in the DFSBC is defined as

The MSE is used as a criterion for determining the convergence of the Newton method.

Shin (2019) was used as the first step, Xm(0) in Newton’s method. The author reported that the real part of a geometric series, f(β) with a common ratio, λe^iβ (where “i” is the imaginary unit), closely resembles the shape of a water wave. A reliable approximation of the profile of a water wave can be achieved by shifting the function down by “f(π/2)” to align with the reference line in Fig. 1, normalizing the result by its height and multiplying the result by the actual wave amplitude “S/2” to meet the wave height condition, as reported by Shin (2019).

The profile satisfies the definition of the reference depth, i.e., ζ(±π/2) = 0, the wave height condition in Eq. (14), and the water depth condition in Eq. (15). Replacing ζ with Xm(0), the profile at the phase β_m is presented as follows: where the steepness S⁽⁰⁾ and the shape factor λ are presented in Shin (2019), and (0) means the first step in Newton’s method. The function is much simpler than the other nonlinear waves because it contains only two parameters, i.e., wave steepness, S and the shape factor, λ. When the shape factor approaches 0, it becomes an Airy wave. When the shape factor approaches 0.5, it becomes a fifth-order Stokes wave. When the shape factor approaches 1, it becomes a Solitary wave. The shape factor is a positive real number less than 1. Because the function was derived independently from the Euler equations, the two parameters were determined to minimize the error in two boundary conditions on the free surface using a variational method. Because Eq. (34) is very close to the complete solution, Newton’s method rapidly and converges for all waves.

In contrast to the assumption considered in Ch. 3, the reference depth is closely linked to the coefficients and is unknown. Therefore, while maintaining the validity of the method presented in Ch. 3, the shooting method is suitable for determining the reference depth. The reference depth is adjusted until the calculated water depth from the water depth condition matches the actual value using the shooting method. The secant method is used as a root-finding algorithm. In each step of the secant method, substituting the calculated coefficients and the prescribed reference depth into Eqs. (A1)–(A4), the wave elevations ζi(q) are calculated with highly dense interval using Newton’s method in the Appendix. The integral is numerically calculated by substituting the results into Eq. (15) as follows: where (q) stands for the q’^th step of the secant method, and ζ_i = ζ(β_i). ζ is an even function. Hence, β₁ = 0; β_i =(i − 1)π/M; and β_M+1 = π. In Fenton (1988), Eq. (35) must be calculated using only N data of X_i because Eq. (35) is coupled to the other equations. Therefore, M = N in Fenton (1988). On the other hand, M is independent of N in this study. Eq. (35) can be calculated accurately compared to Fenton (1988) because M can be freely increased in this study. The water depth is calculated as follows:

The reference depth is calculated using the secant method as follows. where where D⁽⁰⁾ is determined using the regression result in Shin (2019), and h is the actual water depth. The actual value is in the vicinity of D⁽⁰⁾. Hence, D⁽¹⁾ is selected with a value close to D⁽⁰⁾, e.g., D⁽¹⁾ = 1.001D⁽⁰⁾. For each step, the error in the water depth is as follows:

This paper reported a numerical method to calculate the Fourier coefficient, steepness, and reference depth parameter. The method is much simpler than Fenton’s method or HAM. The partial derivatives with respect to wavelength and reference depth, as well as certain parameters, were eliminated utilizing a dimensionless coordinate system. The water depth condition was calculated independently using the shooting method in this study.