# Numerical Method for Calculating Fourier Coefficients and Properties of Water Waves with Shear Current and Vorticity in Finite Depth

## Article information

J. Ocean Eng. Technol. 2023;37(6):256-265
Publication date (electronic) : 2023 December 20
doi : https://doi.org/10.26748/KSOE.2023.034
1Engineer, Offshore structure design department, Hanwha Ocean Co., LTD, Geoje, Korea
Corresponding author JangRyong Shin: +82-55-735-5117, jrshin9797@hanwha.com
Received 2023 October 13; Revised 2023 December 01; Accepted 2023 December 08.

## Abstract

Many numerical methods have been developed since 1961, but unresolved issues remain. This study developed a numerical method to address these issues and determine the coefficients and properties of rotational waves with a shear current in a finite water depth. The number of unknown constants was reduced significantly by introducing a wavelength-independent coordinate system. The reference depth was calculated independently using the shooting method. Therefore, there was no need for partial derivatives with respect to the wavelength and the reference depth, which simplified the numerical formulation. This method had less than half of the unknown constants of the other method because Newton’s method only determines the coefficients. The breaking limit was calculated for verification, and the result agreed with the Miche formula. The water particle velocities were calculated, and the results were consistent with the experimental data. Dispersion relations were calculated, and the results are consistent with other numerical findings. The convergence of this method was examined. Although the required series order was reduced significantly, the total error was smaller, with a faster convergence speed.

## 1. Introduction

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−3 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, −LxL. 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.

Typical wave profile in the coordinate system

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).

## 2. The Complete Solution to the Euler Equations

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 except for those defined separately, such as S and θ. Using the notation, the stream function is defined by ψ = ψ̄ω/k2.The horizontal and vertical velocities are defined by u = Cū and v = Cv̄, respectively, where C = ω/k is the celerity. The vorticity is defined by Ω = ωΩ̄. The pressure is defined by p = ρC2, where ρ is the water density. The dimensionless wave height (steepness) is defined by S = kH. The linear steepness θ = ω2H/g (where g is the gravity) is a constant for a particular wave. When h → ∞ and H → 0, Sθ, which gives the dispersion relation (ω2 = gk) of deep water linear waves. The stream function to satisfy the incompressible Euler equations is as follows:

(1) ψ¯=αc0sinh0α0+n=1NcnsinhnαcoshnDcosnβ
where ndef̳n2+ε; N is the series order; and cn 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:
(2) Ω¯=ε{c0sinh0α0n=1NcnsinhnαcoshnDcosnβ}

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 Ω̄(β,γ) = C1ε, where C1 is a constant. Substituting Ω̄(β,γ) = C1ε into Eq. (2), the wave profile is determined to satisfy the KFSBC (Kinematic free surface boundary condition) as follows:

(3) c0(sinh0γsinh0D)0=n=1Ncn{sinhnγcoshnDcosnβtanhnDcosnπ2}

The constant C1 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.

(4) u¯=1c0cosh0α+n=1NncncoshnαcoshnDcosnβ

The vertical velocity is as follows:

(5) v¯=n=1NncnsinhnαcoshnDsinnβ

Substituting Eqs. (1)(5) to the Euler momentum equations, Bernoulli’s principle for rotational flow can be calculated as follows.

(6) U¯22+Sαθ+p¯=Q
where Ū2 = V̄2ε (αψ̄)2, 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 Ū2 = V̄2ε (αψ̄)2, 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, Ū2 = V̄2 = (1 − ū)2 + v̄2, 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, Ū2 = V̄2ε(αψ̄)2, the dynamic pressure (kinetic energy) Ū2 is as follows:
(7) U¯2=c022n=1NΛ0ncn+n=1Ncn{m=1NΛnmcm}
where
(8) Λ0n(β,α)def̳c0cosnβcoshnD{ncosh0αcoshnα0sinh0αsinhnα}
and
(9) Λnm(β,α)def̳Z1cosh(n+m)αcos(nm)β+Z2cosh(nm)αcos(nm)β+Z3cosh(n+m)αcos(n+m)β+Z4cosh(nm)αcos(n+m)β
where Z1 = (〈n〉〈m〉 + nm − ε)/4W; Z2 = (〈n〉〈m〉 − nm + ε)/4W; Z3 = (〈n〉〈m〉 − nm − ε)/4W; Z4=nm+nm+ε4W; and W = cosh〈nD cosh〈mD. If Λ0n is the component of column vector {Λ0n} and cn is the component of the column vector {cn}, 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 {cn}t[Λnm]{cm}, where {cn}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 c0 is determined as follows:
(10) c0=1U¯b

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

(11) σ=1khcos1(1U¯s1U¯b)for|1U¯s1U¯b|1

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., (π/2, D) = 0 because the point (π/2, D) is on the free surface. Therefore, the pressure field is

(12) p¯=12{U¯2(β,α)U¯2(π2,D)}S(αD)θ

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

(13) Sζθ=12{U¯2(β,γ)U¯2(π2,D)}

The wave height condition is defined as

(14) S=ζ(0)ζ(π).

The water depth condition is defined as

(15) kh=12πππγdβ

## 3. Newton’s Method for Coefficients and Steepness

The solution contains N+2 unknown constants: N Fourier coefficients, cn, 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=0cncosnx, 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, cn, because cosnxm and f(xm) 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, Xm, are introduced in this chapter to simplify the numerical formulation. Unlike Fenton’s method, Xm are merely parameters for calculating the coefficients and the steepness in this chapter. After determining the coefficients, Xm are no longer used. Hence, the wave elevation is denoted with Xm 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, β1 = 0 and βN = π. Denoting ζ(βm) as Xm, i.e., Xm = ζ(βm), Eq. (3) is presented.

(16) Kmncn=Ym

The summation convention is considered in Eq. (16). The repeated subscript “n” is a dummy subscript, whereas Kmn is a second-order tensor, and the terms cn and Xm are vectors in N-dimensional space. From Eq. (3), the component of the tensor Kmn is presented as

(17) Kmn=sinhnX¯mcoshnDcosnβmtanhnDcosnπ2
where m = Xm + D. Because 〈0〉= σi, the component of the first-order tensor, Ym is presented as
(18) Ym=c0{sinσX¯msinσD}σ

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

(19) cn=GnmYm

Eq. (13) is also evaluated at the same phase. A set of N nonlinear equations for calculating the N parameters Xn 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 = X1XN in this approach. For Newton’s method, the error vector Em is defined from Eq. (13) as follows:

(20) Epdef̳SXpθ+Apncn12Bpnmcncm

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

(21) Apndef̳Λ0n(βp,X¯p)Λ0n(π2,D)
(22) Bpnmdef̳Λnm(βp,X¯p)Λnm(π2,D)

Substituting Eq. (19) into Eq. (20)Ep = 0 gives a set of nonlinear equations for calculating the parameters, Xq. 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 Xq points in terms of ΔXq and ignoring higher-order terms, the following set of linear equations is obtained:

(23) Ep,qΔXq=Ep

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

(24) Xq(r+1)=Xq(r)+ΔXq(r)

The superscript (r) means the step of Newton’s method. It is clear that all the steps are rth 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 Xq, the partial derivative Ep,q of the error vector is

(25) Ep,q=S,qXpθSδpqθ+Apn,qcn+Apncn,q12Bpnm,qcncmBpnmcn,qcm.

Note that cn,qcm = cncm,q Because S = X1XN,

(26) S,q={1forq=10forq1orN1forq=N}.

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

(27) Kmn,pcn+Kmncn,p=Ym,p.

Multiplying Eq. (27) by the tensor Gim, the partial derivative of the coefficient is determined as follows.

(28) ci,p=Gim(Ym,pKmn,pcn),
where GimKmn = δin and δincn,p = ci,p and δin is the second-order isotropic tensor. Differentiating Eq. (17) with respect to Xp, the components of the third-order tensor Kmn,p are presented as
(29) Kmn,p=δmpncoshnX¯mcoshnDcosnβm

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

(30) Ym,p=δmpc0cosσX¯m.

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

(31) Apn,q=δpqn2sinhnX¯pcosσX¯pcosnβpcoshnD

Differentiating Eq. (22) with respect to Xq, the component of the forth-order tensor Bpnm,q is

(32) Bpnm,q=δpq[Z1*sinh(n+m)X¯pcos(nm)βp+Z2*sinh(nm)X¯pcos(nm)βp+Z3*sinh(n+m)X¯pcos(n+m)βp+Z4*sinh(nm)X¯pcos(n+m)βp].
where Z1*=Z1(n+m); Z2*=Z2(nm); Z3*=Z3(n+m); and Z4*=Z4(nm).

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, Em by N, the mean squared error (MSE) in the DFSBC is defined as

(33) MSE=EmEmN

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 (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:

(34) Xm(0)=S(0)21λ2λ(1λcosβm1+λ22λcosβm11+λ2)
where the steepness S(0) 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.

## 4. Shooting Method for Reference Depth

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:

(35) (kη0)(q)=12Mi=1M(ζi(q)+ζi+1(q))
where (q) stands for the q’th step of the secant method, and ζi = ζ(βi). ζ is an even function. Hence, β1 = 0; βi =(i − 1)π/M; and βM+1 = π. In Fenton (1988), Eq. (35) must be calculated using only N data of Xi 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:
(36) h(q)=1k{D(q)(kη0)(q)}

The reference depth is calculated using the secant method as follows.

(37) D=limqD(q+1)
where
(38) D(q+1)=D(q1)+D(q)D(q1)h(q)h(q1){hh(q1)}
where D(0) 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(0). Hence, D(1) is selected with a value close to D(0), e.g., D(1) = 1.001D(0). For each step, the error in the water depth is as follows:
(39) Error(q)=|hh(q)|h

## 5. Calculation Procedure and Results

### 5.1 Calculation Procedure

Fig. 2 presents a flowchart. While there is only one do-loop in Fenton’s method, the chart comprises three do-loops. The outer do-loop implements the secant method for calculating the reference depth, as described in Ch. 4. The first inner do-loop implements the Newton method for calculating the coefficients and steepness, as described in Ch. 3. The other inner do-loop implements Newton’s method for calculating the wave profile, as described in the Appendix.

Calculation procedure

Initial approximations of the shape factor, steepness, and reference depth are first calculated by substituting the period, height, and depth into the regression analysis results presented by Shin (2019). The initial approximations of the wave profile are calculated using Eq. (34).

The first inner do-loop is performed. The coefficients, cn(r,q), are calculated using Eq. (19), where “r” represents the step of the Newton method in Ch. 3 and “q” represents the step of the secant method in Ch. 4. Substituting cn(r,q) and D(q) into Eqs. (20) and (25), ΔXi(r,q) are calculated using Eq. (23) and then Xi(r+1,q) are calculated using Eq. (24). The MSE is calculated with Eq. (33). If the MSE is greater than the tolerance, Xi(r,q) is replaced with Xi(r+1,q). The do-loop continues until Newton’s method converges to the complete solution within the specified tolerance.

The second inner do-loop is then performed to calculate the free surface elevations at more phases than those considered in the first inner do-loop. The elevations, ζi(r+1,q), are calculated using Newton’s method in the Appendix, where “r” represents the step of Newton’s method.

By substituting ζi(q)(:limrζi(r+1,q)) into Eq. (35), the water depth h(q) is calculated using Eq. (36). The reference depth D(q+1) is calculated using Eq. (38). Replacing D(q) with D(q+1), the outer do-loop should be repeated until the calculated depth converges to the actual depth. The first step solution of Newton’s method in the q+1 step of the secant method, Xm(0,q+1) is determined using the converged value, Xm(,q)def̳limrXm(r,q) in the q step of the secant method, i.e., Xm(0,q+1)=Xm(,q). Because D(0) is accurate, the secant method converges within three steps to the acceptable value whose error is less than 10−4 percent. The outputs are the coefficients, the reference depth, and the steepness.

### 5.2 Application Range of this Study

The coefficients and wave properties are functions of two variables: the reference depth and the linear steepness. The functions are presented with closed-form solutions in Shin (2019), which are the regression results determined using numerically calculated data for 2091 waves, as shown below. The domain of the functions is {(θ, D) | 0 ≤ θ < 1 and 0 < D}. A subset of the domain, {(θ, D) | 0 ≤ θ < 1 and 0 < D < 5}, was considered instead of the entire domain because waves are similar to deep water waves when D > 5. Seventeen cases of the reference depth were considered by dividing the domain by D = 0.1, 0.15, 0.2, ..., 1.0, 1.2, 1.5, 2, 3, 4, and 5. For each reference depth, the linear steepness was considered from 0 to 0.98 at 0.01 intervals, from 0.98 to 0.989 at 0.001 intervals, from 0.989 to 0.9899 at 0.0001 intervals, and from 0.9899 to 0.98997 at 0.00001 intervals, resulting in 123 cases of linear steepness. Therefore, 17 × 123 = 2091 irrotational waves were calculated with N = 35 and M = 180. The limits can be determined easily because the Newton method presented in the Appendix diverges at the crest when the steepness reaches its limit. These limits were represented by 17 red solid circle symbols in Fig. 3, which are calculated with D = 0.1, 0.15, 0.2, ..., 1.0, 1.2, 1.5, 2, 3, 4, and 5 in sequence from left to right. The solid curve is the regression analysis result determined using the 17 data. As shown in Fig. 3, the breaking limit is in good accordance with the Miche formula. Because the 2091 waves are dimensionless waves, they represent all waves below the breaking limit. Therefore, this study confirms stable convergence characteristics for all water waves.

Application range of this study cited from Shin (2019)

### 5.3 Water Particle Velocities Under the Crest

The eight waves considered in Le Méhauté et al. (1968) were calculated. Fig. 3 and Table 1 present the waves with solid diamond symbols. The results were compared with experimental data in Fig. 2 by Shin (2022). Fig. 4 shows the wave (h) during the eight waves. Referring to Fig. 2 of Shin (2022), the other waves give the same pattern, as shown in Fig. 4. The irrotational results were calculated with ε = 0. The regression analysis results were calculated by Shin (2019).

Eight waves considered in Le Méhautéet al. (1968)

Water particle velocities under the crest of the wave (h)

Rienecker and Fenton (1981) also calculated the waves, as shown in Fig. 2 of that work. The waves were also calculated by Fenton (1990), and the results are shown in Fig. 3(c) of that work. The cnoidal theory is the new cnoidal theory by Fenton (1990), which closely agrees with the Fourier approximation by Fenton (1988). Tao et al. (2001) also calculated the waves; the results are shown in Fig. 2 of that work. HAM (Homotopy analysis method) by Tao et al. (2001) results closely agree with the Fourier approximation reported by Rienceker et al. (1981). Therefore, the Fourier approximation and HAM results are not presented in Fig. 4. This study agrees with the experimental data, unlike those obtained using HAM or Fenton’s method.

### 5.4 Dispersion Relation

Table 3 of Tao et al. (2007) provides a detailed comparison of the HAM solutions of waves in finite water depth and the results of Cokelet (1977) and Rienecker and Fenton (1981). The results are given in terms of the non-dimensional phase speed squared kC2/g at various values of H/h, for a constant value of exp(−kR/C) = 0.5 (R is the unit span denoting the total volume rate of flow underneath the steady wave per unit length in a direction normal to the x, z plane), corresponding to a wavelength to water depth ratio of L/h ≈ 9 (Tao et al., 2007). Table 1 of Rienecker and Fenton (1981) presents the same comparison between the solutions of Rienecker and Fenton (1981) for N = 16, 32, and 64, the results of Cokelet (1977), and the results of Vanden-Broeck & Schwarts.

Errors of the first step

The non-dimensional phase kC2/g is presented with physical quantities defined in this study as follows:

(40) kC2g=kg(ωk)2=ω2gk=θkH=θS

H/h is also presented with physical quantities defined in this study.

(41) Hh=kHkh=Skh
kh = 2π/9 in Eq. (41) because L/h ≈ 9. Therefore, the comparison shows the relationship between steepness, S, and linear steepness, θ for the relative water depth, kh = 2π/9. Hence, the comparison shows the dispersion relation of waves in a finite water depth, kh = 2π/9. Table 2 lists the relation. The data in the first and second columns were cited from Rienecker and Fenton (1981), which are the results calculated using N = 64. The data in the fourth column was obtained by multiplying the data in the second column by 2π/9, as stated in Eq. (41). The data in the third column was obtained by multiplying the data in the first column by the data in the fourth column, as stated in Eq. (40). Shin (2018, 2019) also calculated the relation for various water depths, D = 0.1, 0.15, 0.2, ... , 1.0, 1.2, 1.5, 2, and ∞. For a given reference depth, the linear steepness θ was considered with a highly dense interval from 0 to the application limit of this study. Fig. 5 presents the calculation results. The 15 curves are the results calculated with D = 0.1, 0.15, 0.2, ... , 1.0, 1.2, 1.5, 2, and ∞ in sequence from left to right. The peaks of each curve represent the limit of wave breaking. The nine data in Table 2 are also presented with solid circles in Fig. 5, which agree with the curve for D = 0.7. The small difference results from the difference between kh = 2π/9 and D = 0.7. The relative depth for D = 0.7, is presented as kh = D0 = 0.7 − 0 > 2π/9. because D = kh + 0 and kη0 is negative. Therefore, the D = 0.7 curve lies slightly to the right of the eight data because 0 is very small.

Dispersion relation by Rienecker and Fenton (1981)

Dispersion relation for irrotational waves in the reference depths, D = 0.1, 0.15, 0.2, ... , 1.0, 1.2, 1.5, 2, and ∞ in sequence from left to right

### 5.5 Error Check of this Study

The numerical method focuses on minimizing the errors in the DFSBC defined in Eq. (20) and the water depth conditions defined in Eq. (15). The convergence of this study according to the series order was analyzed by calculating the eight waves considered in Le Méhauté et al. (1968) using N = 8, 10, 16, and 35. Fig. 6 presents the calculated wave profiles for the waves. The abscissa represents the phase β = kxωt in −180° ≤ β ≤ 180°, while the ordinate represents the dimensionless elevation, . The errors are summarized in Tables 3 and 4. The first step is the result of Shin (2019). Table 3 lists the error of Shin (2019). There are three steps for calculating the water depth condition using the secant method and three steps for the Newton method at each secant step. Therefore, there are 10 steps; Table 4 presents the errors for step 10. The blanks indicate that they are unsuitable. The series order, N = 8, is unsuitable for waves (c) and (d), and the series order, N = 10, is unsuitable for wave (d) because they do not satisfy the condition for regular waves. Hence, there is no local extrema between the crest and the trough in a regular wave. The deeper the water depth and the lower the wave height, the smaller the number of series orders required. If the series order is increased beyond the optimal number, the convergence of this study decreases, as presented by Shin (2023). Even for very long waves, the convergence is very stable, unlike other methods.

Wave profiles (The legend was arranged in descending order of the crest height. The wave with the highest crest is wave (e), and the wave with the lowest crest is wave (d).)

Errors of the step 10

## 6. Conclusions

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.

## Notes

The authors declare that they have no conflict of interest.

## References

Baddour R. E, Song S. W. 1988;The rotational flow of finite amplitude periodic water waves on shear currents. Applied Ocean Research 20(3):163–171. https://doi.org/10.1016/S0141-1187(98)00015-7 .
Chaplin J. R. 1979;Developments of stream-function wave theory. Coastal Engineering 3:179–205. https://doi.org/10.1016/0378-3839(79)90020-6 .
Chappelear J. E. 1961;Direct numerical calculation of wave properties. Journal of Geophysical Research 66(2):501–508. https://doi.org/10.1029/JZ066i002p00501 .
Cokelet E. D. 1977;Steep gravity waves in water of arbitrary uniform depth. Philosophical Transaction of the Royal Society of London. Series A, Mathematical and Physical Sciences 286(1335):183–230. https://doi.org/10.1098/rsta.1977.0113 .
Dean R. G. 1965;Stream function representation of nonlinear ocean waves. Journal of Geophysical Research 70(18):4561–4572. https://doi.org/10.1029/JZ070i018p04561 .
Dean R. G, Dalrymple R. A. 1984. Water wave mechanics for engineers and scientists Prentice-Hall, Inc. p. 368.
Fenton J. D. 1988;The numerical solution of steady water wave problems. Computers & Geosciences 14(3):357–368.
Fenton J. D. 1990. Nonlinear wave theories. In : Le Méhauté B, Hanes D. M, eds. The Sea, 9: Ocean Engineering Science, Part A Wiley.
Le Méhauté B, Divoky D, Lin A. 1968. Shallow water waves: A comparison of theories and experiment. Proceedings of 11th Conference Coastal Engineering London, UK: ASCE. 1p. 86–107.
Rienecker M. M, Fenton J. D. 1981;A fourier approximation method for steady water waves. Journal of Fluid Mechanics 104:119–137. https://doi.org/10.1017/S0022112081002851 .
Shin J. 2018. An analytical approximation and a numerical method to progressing periodic water waves on irrotational flow over a horizontal bed. The 28th ISOPE Conference, Sapporo, Japan (ISOPE-I-18-567)
Shin J. 2019;A regression analysis result for water waves on irrotational flow over a horizontal bed. International Journal of Offshore and Polar Engineering 29(4):461–466. https://doi.org/10.17736/ijope.2019.hc17 .
Shin J. 2022;Solution for water waves with a shear current and vorticity over a flat bed. International Journal of Offshore and Polar Engineering 32(4):418–423. https://doi.org/10.17736/ijope.2022.hc29 .
Shin J. 2023;Newton’s method to determine Fourier coefficients and wave properties for deep water waves. Journal of Ocean Engineering and Technology 37(2):49–57. https://doi.org/10.26748/KSOE.2022.037 .
Tao L, Song H, Chakrabarti S. 2007;Nonlinear progressive waves in water of finite depth -an analytic approximation. Coastal Engineering 54(11):825–834. https://doi.org/10.1016/j.coastaleng.2007.05.008 .
Vanden-Broeck J. M, Schwarts L. W. 1979;Numerical computation of steep gravity waves in shallow water. Physic and Fluids 22(10):1868–1871. https://doi.org/10.1063/1.862492 .

## Appendices

Appendix

#### Newton’s Method for the Wave Profile

Note that the coefficient cn and the reference depth D are known in this Appendix. An implicit function f(β, ζ) = 0 can be considered as an equation with respect to the dependent variable ζ because the independent variable β is known. Thus, using Newton’s method, the explicit representation of the function is

(A1) ζ(β)=limrζ(r+1)(β)
where (r) stands for the step of Newton’s method in this Appendix and
(A2) ζ(r+1)(β)=ζ(r)(β)f(β,ζ(r))f(β,ζ(r))
where f(β, ζ) is defined as follows using Eq. (3):
(A3) f(β,ζ)=c0(sinσγsinσD)σ+n=1Ncn{sinhnγcoshnDcosnβtanhnDcosnπ2}

Because 0 ≤ |ζ| < 1 for all waves, the first step solution ζ(1) is calculated with the power series expansion of f in ζ and ignoring the higher order terms than the second order. As f(β, ζ)= 0, the power series expansion is a quadratic equation with respect to ζ. Then we have

(A4) ζ(1)=f(β,0){f(β,0)}22{f(β,0)}{f(β,0)}f(β,0)

The other root of the quadratic equation is trivial. Eq. (A4) has the limit that is determined with the discriminant of the quadratic equation. On the other hand, the limit is greater than the wave-breaking limit. Therefore, Eq. (A4) is valid for all waves. Iso-streamlines and iso-vorticity lines are calculated by applying the method presented in this Appendix to Eq. (1) and Eq. (2), respectively. An iso-streamline for the specified value of ψ̄ = constant can be obtained by applying Eq. (A5) to the Newton method instead of Eq. (A3) and replacing ζ with α.

(A5) f(β,α)=ψ¯+αc0sinh0α0+n=1NcnsinhnαcoshnDcosnβ

The iso-vorticity line for the specified value of Ω̄ = constant can also be obtained by applying Eq. (A6) to the Newton method instead of Eq. (A3) and replacing ζ with α.

(A6) f(β,α)=Ω¯+ε{c0sinh0α0n=1NcnsinhnαcoshnDcosnβ}

## Article information Continued

### Fig. 1.

Typical wave profile in the coordinate system

### Fig. 2.

Calculation procedure

### Fig. 3.

Application range of this study cited from Shin (2019)

### Fig. 4.

Water particle velocities under the crest of the wave (h)

### Fig. 5.

Dispersion relation for irrotational waves in the reference depths, D = 0.1, 0.15, 0.2, ... , 1.0, 1.2, 1.5, 2, and ∞ in sequence from left to right

### Fig. 6.

Wave profiles (The legend was arranged in descending order of the crest height. The wave with the highest crest is wave (e), and the wave with the lowest crest is wave (d).)

### Table 1.

Eight waves considered in Le Méhautéet al. (1968)

 Wave (a) (b) (c) (d) H/h 0.434 0.420 0.389 0.433 Tg/h 8.590 15.870 22.490 27.240 Wave (e) (f) (g) (h) H/h 0.499 0.522 0.493 0.548 Tg/h 8.590 15.870 22.502 27.270

### Table 2.

Dispersion relation by Rienecker and Fenton (1981)

kC2/g H/h θ S
0.1729974 0.615059 0.074284 0.120775
0.2526308 0.631112 0.111309 0.176370
0.3802643 0.666501 0.176939 0.265475
0.4944549 0.706443 0.243860 0.345195
0.6024470 0.748230 0.314696 0.420587
0.6512510 0.764403 0.347543 0.454659
0.6721430 0.767760 0.360267 0.469244
0.6832000 0.767030 0.365845 0.476964
0.6908000 0.763000 0.367972 0.482269

### Table 3.

Errors of the first step

Wave Root mean squire error in the DFSBC (%)

N = 8 N = 10 N = 16 N = 35
a 2.30E-04 2.31E-04 2.31E-04 2.31E-04
b 4.03E-03 4.00E-03 4.00E-03 4.00E-03
c 1.19E-03 1.18E-03 1.18E-03
d 1.18E-03 1.19E-03
e 1.14E-03 1.14E-03 1.14E-03 1.14E-03
f 3.66E-03 3.57E-03 3.55E-03 3.55E-03
g 1.91E-03 1.57E-03 1.58E-03
h 1.68E-03 1.77E-03
Wave Relative error in the water depth condition (%)

N = 8 N = 10 N = 16 N = 35
a 1.35E-04 5.50E-07 9.45E-11 2.73E-12
b −3.32E-01 −3.30E-01 −3.29E-01 −3.29E-01
c 5.54E-03 −1.57E-05 −1.76E-08
d −8.19E-04 3.32E-06
e 2.79E-04 −2.74E-06 7.34E-10 2.62E-11
f −1.05E-02 −5.16E-03 1.36E-05 −9.90E-12
g 2.40E-02 −6.74E-04 −2.79E-07
h −8.60E-03 4.36E-05

### Table 4.

Errors of the step 10

Wave Root mean squire error in the DFSBC (%)

N = 8 N = 10 N = 16 N = 35
a 1.78E-09 2.88E-14 8.27E-20 1.02E-21
b 5.19E-07 2.41E-07 6.85E-12 3.62E-25
c 5.02E-06 8.41E-10 1.78E-19
d 6.11E-08 3.11E-15
e 2.37E-08 6.18E-11 2.69E-15 8.91E-21
f 1.44E-05 4.17E-07 7.85E-10 3.28E-20
g 1.01E-04 6.97E-09 1.83E-16
h 3.63E-08 1.19E-12
Wave Relative error in the water depth condition (%)

N = 8 N = 10 N = 16 N = 35
a 9.60E-05 9.59E-05 9.60E-05 9.60E-05
b 2.12E-04 2.25E-04 2.09E-04 2.09E-04
c 2.11E-04 4.20E-04 4.20E-04
d 2.75E-05 2.68E-05
e 1.54E-05 1.51E-05 1.51E-05 1.51E-05
f 1.57E-04 7.80E-05 9.23E-05 9.26E-05
g 2.18E-05 4.08E-04 4.05E-04
h 9.40E-05 8.56E-05