# Newton’s Method to Determine Fourier Coefficients and Wave Properties for Deep Water Waves

## Article information

J. Ocean Eng. Technol. 2023;37(2):49-57
Publication date (electronic) : 2023 March 6
doi : https://doi.org/10.26748/KSOE.2022.037
1Engineer, Offshore structure design department, Daewoo Shipbuilding & Marine Engineering co., LTD, Geoje, Korea
Corresponding author JangRyong Shin: +82-55-735-5117, jrshin@dsme.co.kr
Received 2022 November 17; Revised 2023 February 2; Accepted 2023 February 13.

## Abstract

Since Chappelear developed a Fourier approximation method, considerable research efforts have been made. On the other hand, Fourier approximations are unsuitable for deep water waves. The purpose of this study is to provide a Fourier approximation suitable even for deep water waves and a numerical method to determine the Fourier coefficients and the wave properties. In addition, the convergence of the solution was tested in terms of its order. This paper presents a velocity potential satisfying the Laplace equation and the bottom boundary condition (BBC) with a truncated Fourier series. Two wave profiles were derived by applying the potential to the kinematic free surface boundary condition (KFSBC) and the dynamic free surface boundary condition (DFSBC). A set of nonlinear equations was represented to determine the Fourier coefficients, which were derived so that the two profiles are identical at specified phases. The set of equations was solved using Newton’s method. This study proved that there is a limit to the series order, i.e., the maximum series order is N=12, and that there is a height limitation of this method which is slightly lower than the Michell theory. The reason why the other Fourier approximations are not suitable for deep water waves is discussed.

## 1. Introduction

The undulations on the water surface adopt inherently beautiful forms, such as the smooth regular features of progressive waves. On the other hand, despite the innate curiosity concerning water wave motion, it was only in the 19th century by Airy that efforts to provide a mathematical description of wave motion began in earnest (Henry, 2008). Since the early development of the theory describing the waves using the perturbation method by Stokes (1847), many nonlinear solutions were obtained, both analytically and numerically.

Based on the use of truncated Fourier expansions for flow field quantities, Chappelear (1961) and Dean (1965) developed a ‘Fourier approximation’ (Tao et al., 2007). By choosing the expansions to satisfy the Laplace equation and the BBC, the problem was reduced to solving a set of nonlinear equations for each of the Fourier coefficients, and the wave properties (Tao et al., 2007). As the set of nonlinear equations is derived from the finite series, the expansions were represented by truncated Fourier expansions, which is why they are coined as Fourier approximations. Considerable research (Chappelear, 1961; Dean, 1965; Chaplin, 1979; Rienecker and Fenton, 1981; Fenton, 1988) has been conducted on this approach.

The steam function theory was developed by Dean (1965). Dean (1965) proposed the use of the root mean square errors (RMSE) in the dynamic free surface boundary condition (DFSBC) and the kinematic free surface boundary condition (KFSBC) as an error evaluation criterion for wave theories. Dean (1965) also defined the Lagrangian function by introducing the linear summation of the two errors on the free surface with a Lagrange multiplier. Fourier coefficients were determined to minimize the Lagrangian function. The required order (N) of the stream function theory is determined by the wave steepness and shallow water parameter. For N = 1, the stream function theory reduces to the Airy wave theory. As the breaking wave height is approached, more terms are required to accurately represent the wave. Det Norske Veritas (DNV, 2010) described that stream function wave theory has an error of more than 1% in deep water waves whose height is over 90% of the breaking limit even though the required order is increased (See Figs. 34 in DNV (2010)). Chaplin (1979) applied the Schmidt orthogonalization process as an alternative to Dean’s method and obtained improved results (Tao et al., 2007).

Flow chart for numerical analysis.

Dimensionless wave profiles calculated with N = 12 (ordinate is the dimensionless elevation from the still water line (SWL), )

Rienecker and Fenton (1981) also adopted the stream function introduced by Dean (1965). The unknown constants were calculated directly with Newton’s method in Rienecker and Fenton (1981). Therefore, the method required partial derivatives with regard to unknown constants, which complicated the numerical formulation of Newton’s method. Fenton (1988) obtained all the partial derivatives for simplification of this method numerically. Because the reference depth parameter is simultaneously determined with the Fourier coefficients, the required order should be increased to reduce the numerical error in the water depth condition. Newton’s method is unstable because the number of unknown coefficients increases when the required order increases. Therefore, the Newton method requires an initial solution closer to the exact solution. To solve it, in common with other versions of the Fourier approximation, it is sometimes necessary to solve a sequence of lower waves, extrapolating forward in height steps until the desired height is reached (Fenton, 1988). This problem can be avoided using the sequence of height steps (Fenton, 1990). As Rienecker and Fenton (1981) pointed out, several aspects of Fourier approximations beyond just the truncation of the series inhibited its widespread use (Tao et al., 2007). Neither of the above stream function approaches can be applied to waves in deep water because the stream function expansions contain hyperbolic functions (Fenton, 1988).

Shin (2016) reported that most numerical errors resulted from denominators in the stream function, the numerical integration of the water depth condition, and problems involved in Newton’s method. Shin (2016) solved the problems by separately calculating the water depth condition and introducing deep water velocity potential without denominators and a dimensionless coordinate system. The flow field was represented with a velocity potential to satisfy the Laplace equation and the bottom boundary condition (BBC). Two wave profiles were calculated by applying the velocity potential to the KFSBC and the DFSBC. The potential contains N+2 unknown constants, which are N Fourier coefficients, the wave steepness, and the reference depth parameter, while the stream function contains 2N+6 unknown constants (Rienecker and Fenton, 1981; Fenton, 1988). The wavelength and the reference depth parameter are not coupled to the Fourier coefficients because the dimensionless coordinate system was defined with the phase and the dimensionless elevation. Therefore, the wavelength and the reference depth parameter were determined independently. As a result, the required order was reduced in Shin (2016) and Shin (2021)

All the above-mentioned advantages resulted from the dimensionless coordinate system (Refer to Figs. 1 and 2). In addition to the advantages, there are some advantages to the coordinate system in Fig. 2. The independent variable is the relative horizontal position with regard to the wavelength in the moving coordinate system (Rienecker and Fenton, 1981; Fenton, 1988). The horizontal positions are also unknown values in Fenton’s method because the wavelength is an unknown constant, which makes the numerical procedure much more complicated. On the other hand, the independent variable is the phase in the range [0, π] in Shin (2016), which is a known value. In addition, partial derivatives with regard to the wavelength are unnecessary, unlike Fenton’s method (Rienecker and Fenton, 1981; Fenton, 1988). All the partial derivatives necessary can be calculated with tensor analysis, while they are calculated numerically in Fenton’s method. Therefore, there are no errors in the partial derivatives, unlike Fenton’s method. The numerical formulation of Newton’s method was also represented by tensor analysis.

Conventional coordinate system for a progressive wave (H: wave height, T: wave period, L : wavelength, η(t, x): free surface elevation from the still water line (SWL)).

Dimensionless coordinate system for a progressive wave.

Some waves near the breaking limit were calculated for verification, and their errors were checked with Dean’s criteria. The deep water breaking limit was checked. According to the required order, the trend of numerical error was checked. There is a limit to the series order, i.e., N ≤ 12. This is why the other approximations are unsuitable for deep water waves.

## 2. Complete Solution

The complete solution was presented by Shin (2016), and the results are summarized in this chapter. Two coordinate systems were adopted to describe a progressive water wave. One is the conventional coordinate system (t, x, y) shown in Fig. 1. The origin is located on the still water line. The x-axis is in the direction of the wave propagation, the y-axis points upwards, and t is time. The fluid domain is bounded by the free surface y = η(t,x).

The other coordinate system (β, α) is the dimensionless stationary frame shown in Fig. 2. The origin is located at the point under the crest on the reference line, which is the horizontal line passing through two points at ηo = η (±90°) on the free surface. Therefore, the wave profile is a fixed, periodic, and even function in the system. The horizontal axis is the phase, β = kxωt, whose domain is −πβπ, where k is the wave number defined by kdef̳2π/L and ω is the angular frequency defined by ωdef̳2π/T . The vertical axis is the dimensionless elevation from the reference line, α = k(yη0) in αζ, where ζ = k(ηη0) is the dimensionless free surface elevation from the reference line.

The velocity potential in finite depth is represented in Shin (2019). The hyperbolic functions in the denominator are divergent as the water depth increases. Therefore, the potential is represented by the original deep water form without the hyperbolic functions.

(1) ϕ=ωk2n=1Nanenαsinnβ
where an is a Fourier coefficient, which is dimensionless. The horizontal water-particle velocity is
(2) u(β,α)=cn=1Nnanenαcosnβ
where c = ω/k is the wave celerity. The vertical water-particle velocity is
(3) v(β,α)=cn=1Nnanenαsinnβ

limαv=0, which is the BBC. The water-particle accelerations in the horizontal direction and the vertical direction are represented as follows:

(4) ut=ω2kn=1Nn2anenαsinnβ
(5) vt=ω2kn=1Nn2anenαcosnβ

The KFSBC is presented by an ordinary differential equation in the coordinate system. Solving the differential equation, the wave profile was calculated as follows (Refer to Appendix A):

(6) ζ=n=1Nan{enζcosnβcosnπ2}

The integral constant was determined by the definition of the reference line, i.e., ζ(±π/2) = 0. Eq. (6) contains only N coefficients, while Rienecker et al. (1981) contains N+5 unknown constants, which are N+1 coefficients, the wave number, celerity, volume rate, and Bernoulli’s constant. Compared to Rienecker et al. (1981), Eq. (6) is simple. Unknown constants are reduced dramatically by introducing the coordinate system in Fig. 2. Bernoulli’s equation is represented in the dimensionless coordinate system as follows:

(7) p(β,α)ρc2=U(β,α)U012{U2(β,α)+V2(β,α)}+12(U02+V02)Sαθ
where ρ is the water density, and U and V are dimensionless horizontal and vertical velocities defined by Udef̳u/c and Vdef̳v/c . Bernoulli’s constant is eliminated by the definition of the reference line, i.e., p(±π/2,0) = 0. U0 = U(±π/2,0) and V0 = V(±π/2,0) are velocities at the phase of β = ±π/2 on the free surface. The linear steepness is defined with θdef̳(ω2H)/g where g is gravity and H is the wave height. The linear steepness is a constant for a particular wave. The other wave profile is calculated by applying the DFSBC (p(β, ζ = 0) to Eq. (7), as follows:
(8) ζSθ=U(β,ζ)U012{U2(β,ζ)+V2(β,ζ)}+12(U02+V02)

Eq. (8) contains only N coefficients and wave steepness. In addition, Eq. (8) also satisfies that ζ(±π/2) = 0. Dimensionless wave height (steepness) is defined by Sdef̳kH, which provides the dispersion relation because the wave number is calculated as k = S/H. The reference line is determined by the water depth condition presented as follows:

(9) kη0=-12π-ππζdβ

The steepness is determined by the wave height condition presented as follows:

(10) S=ζ(0)ζ(π)

For N → ∞, Eqs. (1)(10) provide the complete solution to irrotational deep water waves. The solution contains N+2 unknown constants, i.e., the Fourier coefficient, an, wave steepness, S, and reference depth parameter, 0.

This study aims to present a method to numerically determine the constants. N+2 equations are necessary to determine the unknown constants. The two profiles in Eqs. (6) and (8) should be identical for all phases. Hence, a set of N equations are obtained so that Eqs. (6) and (8) are equal at N phases. In addition, there are two equations, i.e., Eqs. (9) and (10). As a result, there are N+2 equations to determine the unknown constants.

An example for N = 1 and for N = 2 is presented in Appendix B and C, and Shin (2016) and Shin (2021) presented an example for N = 3. This method is further simplified and generalized by tensor analysis in the next chapter.

## 3. The Fourier Coefficients and the Steepness

A method to determine the Fourier coefficients and the steepness is presented for an arbitrary number of N.

In the other Fourier approximations (Rienecker and Fenton (1981), Fenton (1988)), the profile, η is directly calculated instead of ζ. Therefore, Eq. (9) must be simultaneously calculated with the coefficients. In this study, however, the profile η was calculated with ζ. Therefore, Eq. (9) was not coupled with the coefficients in this study and was calculated after they were determined. Hence, Eq. (9) is not considered in this chapter.

When the wave profile, ζ is prescribed in advance, it is possible to convert Eq. (6) to a set of linear equations to determine the coefficients. Because the wave profile is an even function, the phases βm for m = 1,2, ...,N are considered in the range, 0 ≤ βmπ and βmπ/2 because Eqs. (6) and (8) already satisfy that ζ(±π/2) = 0. β1 = 0 and βN = π. Letting Xm = ζ(βm), Eq. (6) is represented at the phase βm as follows:

(11) Kmnan=Xm

The summation convention is considered in Eq. (11). The repeated subscript “n” is a dummy subscript. Kmn is a second-order tensor, an and Xm are vectors in N-dimensional space. The component of the tensor Kmn is presented as follows:

(12) Kmn=enXmcos(nβm)cos(nπ2)

Using the inverse tensor Rmn of the tensor Kmn, the solution to Eq. (11) is determined easily as follows:

(13) an=RnmXm

Using Eq. (8), the error vector Em in the DFSBC is defined as

(14) Em=XmSθ+Umnan12εmpqUpnanU¯qjaj12εmnqVpnanV¯qjaj
where S = X1 — XN from the wave height condition in Eq. (10). Therefore, the wave height condition is automatically integrated in Eq. (14). The components of the second-order tensors in Eq. (14) are represented as follows:
(15) Umn=n{enXmcos(nβm)cos(nπ2)}
(16) U¯mn=n{enXmcos(nβm)+cos(nπ2)}
(17) Vmn=n{enXmsin(nβm)sin(nπ2)}
(18) V¯mn=n{enXmsin(nβm)+sin(nπ2)}

The third-order tensor εijp is defined as follows:

(19) εpqr={1for p=q=r0fortheothercase

A set of N nonlinear equations represented by Em = 0 for Xm are obtained by substituting Eq. (13) into Eq. (14). The set of equations is solved with Newton’ s method. A set of linear equations can be derived by denoting partial derivatives of a tensor using commas and indices as ()Xi=(),i and applying Newton’s method to Eq. (14).

(20) Em,iΔXi=Em

Because ΔXi[r]=Xi[r+1]Xi[r], the solution in the next step is

(21) Xi[r+1]=Xi[r]+ΔXi[r]

The superscript [r] with [ ] means the step of Newton’s method. All the steps are rth in all equations except Eq. (21) in this chapter. Therefore, for the simplification of equations, the superscript [r] was omitted in all the equations except Eq. (21). Differentiating Eq. (14) with respect to Xi, the partial derivative Em,i of the error vector is

(22) Em,i=δmiSθXmS,iθ+Umn,ian+Umnan,i12εmpqUpn,ianU¯qjaj12εmpqUpnan,iU¯qjaj12εmpqUpnanU¯qj,iaj12εmpqUpnanU¯qjaj,i12εmpqVpn,ianV¯qjaj12εmpqVpnan,iV¯qjaj12εmpqVpnanV¯qj,iaj12εmpqVpnanV¯qjaj,i

Differentiating Eqs. (15)(18) with respect to Xi, the components of the third-order tensors are

(23) Umn,q=U¯mn,q=n2δmqenXmcos(nβm)
(24) Vmn,q=V¯mn,q=n2δmqenXmsin(nβm)

Because S = X1XN,

(25) S,i={1fori=10fori1orN1fori=N}

Differentiating Eq. (11) with respect to Xp, we have

(26) Kmn,pan+Kmnan,p=δmp
where Xm p = δmp and δmp is the second-order isotropic tensor. Multiplying Eq. (26) by the tensor Rim, the partial derivative of the coefficient is easily determined as follows:
(27) ai,p=RipRimKmn,pan
where RimKmn = δin and δinan,p = ai,p. Differentiating Eq. (12) with respect to Xp,
(28) Kmn,p=δmpnenXmcosnβm

Therefore, all the partial derivatives necessary were easily obtained without errors, unlike Fenton (1988).

Representing an individual component of a tensor as shown in Eqs. (12), (15)(18), (23)24), and (28), the summation convention was not considered. As a criterion of convergence, the Root Mean Square Error in the DFSBC proposed by Dean (1965) was adopted. Dividing the inner product of the error vector Em by N, the RMSE in the DFSBC (Dean, 1965) is defined as follows:

(29) RMSE=EmEmN×100%

Newton’s method was rapidly convergent to the complete solutions (whose RMSE ≤ 10−13% within three steps) when the result of Shin (2021) was used as the first step solution. There are some differences in the above approach compared to Fenton’s method (Rienecker et al., 1981; Fenton, 1988).

• Eq. (11) is linear and decoupled from Eq. (20).

• All the partial derivatives with respect to wavelength are not required.

• All equations can be formulated with tensor analysis.

• Xi are merely parameters for calculating the coefficients and the steepness in this study. After determining the coefficients, Xt are no longer used. Hence, the wave elevation is denoted with Xi instead of ζi in Eqs. (11)(28).

• After determining the coefficients and the steepness, the wave profile and water depth condition were calculated separately, as presented in the next chapter.

## 4. Wave Profile and Reference Depth Parameter

The Fourier coefficients and the steepness were determined in advance. Therefore, the wave profile can be calculated with Eq. (6) as follows: Eq. (6) is represented as

(30) F(β,ζ)=ζ+n=1Nan{enζcosnβcosnπ2}
where F(β, ζ) = 0. The power series of F(β, ζ) is represented as follows:
(31) F(β,ζ)=q=0F(q)(β,0)q!ζq
where F(q)=qFζq.

Because F(β ζ) = 0, Eq. (31) provides a polynomial equation with regard to ζ. For q = 1, Eq. (31) provides a linear equation. Eq. (31) provides a quadratic equation for q = 2, and Eq. (31) provides a cubic equation for q = 3. The equations provide good approximations to the wave profile because |ζ| < 1 for all cases. For q=2, the following approximation is obtained:

(32) ζ[1]=F(1)(β,0){F(1)(β,0)}22F(β,0)F(2)(β,0)F(2)(β,0)

The other root is a trivial solution because it does not satisfy that |ζ| < 1. The approximation is used as the first step in the Newton’s method defined as follows. The superscript [n] means the nth step of Newton’s method (where n is a natural number), while the superscript (n) means the nth order partial derivative with respect to ζ. The wave profile is calculated as follows by applying Newton’s method:

(33) ζ=limmζ[m+1]
where
(34) ζ[m+1]=ζ[m]F(β,ζ[m])F(β,ζ[m])

|F‘(β, ζ)| ≥ |F’(0, ζc)| for all waves and |F‘(0, ζc)| = 0 is the breaking condition proposed by Stokes (Chakrabarti, 1987). Therefore, Eq. (34) is valid for all waves except breaking waves. The Newton method rapidly converges to the complete solution. The wave elevations ζi = ζ (βi) were calculated using Eq. (33). Note that ζi stands for the free surface elevation at phase βi. The integral is numerically calculated by substituting the results in water depth conditions as follows:

(35) kη0=12Mi=1M(ζi+ζi+1)

Because ζ is an even function, β1 = 0, βi = (i − 1)π/M and βM+1 = π. Note that M is independent of N. When M is increased, Eq. (35) can be calculated more accurately rather than Rienecker et al. (1981) and Fenton (1988).

## 5. Numerical Analysis Procedure

Fig. 3 presents a flow chart of the procedure. The chart comprises two DO-loops because the Fourier coefficients and wave profile are independently calculated. The first step, Xm[1], was calculated using the result in Shin (2021). The first DO-loop calculates the coefficients and the steepness. The coefficients, an[r] are calculated by Eq. (13). Substituting, an[r] into Eq. (20), ΔXi[r] are calculated by Eq. (20) and then Xi[r+1] are calculated by Eq. (21). The RMSE is calculated with Eq. (29). If the RMSE is greater than the tolerance (1 × 10−13%), Xi[r] is replaced with Xi[r+1]. The DO-loop continues until Newton’s method converges to the complete solution within the tolerance.

In this DO-loop, the Fourier coefficients and the steepness are determined. The wave number is calculated as follows:

(36) k=SH

The second DO-loop is used to calculate the free surface elevations at more phases than the phases considered in the first DO-loop. The elevations, ζi[m+1] are calculated using Newton’s method in Eq. (34). The DO-loop continues until Newton’s method converges to the complete solution within the tolerance (|F(βi, ζi)| ≤ 1 × 10−17) . Substituting ζi into Eq. (35), the reference depth parameter o was calculated. With the results, the wave profile was calculated as follows:

(37) ηi=ζi+kη0k

The water particle velocities were calculated by substituting the Fourier coefficients into Eqs. (2)(3), and the accelerations were calculated by substituting the Fourier coefficients into Eqs. (4)(5). The pressure field was calculated by substituting the Fourier coefficients into Eq. (7). The other wave profile was calculated by substituting the results into the right side of Eq. (8). The error in DFSBC can be checked by comparing the two profiles.

## 6. Results

The Fourier coefficients, an, the steepness, S, and the reference depth parameter, η0 are functions of one variable whose independent variable is the linear steepness, θ. Consequently, Shin (2019) numerically calculated some data for the coefficients, steepness, and reference depth parameters in the range, 0 < θ < 1. By curve fitting the data, the Fourier coefficients, steepness, and reference depth parameter are represented by Newton’s polynomials in Shin (2021), which give closed-form solutions. In this calculation, N = 3 was considered. The results satisfied the Laplace equation, the BBC, and the KFSBC. The RMSE in the DFSBC was less than 1% in the range, H/Lo ≤ 0.142 where Lo is the linear wavelength defined as follows:

(38) Lo=2πgω2

Because θ = 2πH/Lo, H/Lo = 0.142 corresponds to θ = 0.892. As a result, Shin (2016, 2021) reported a good approximation whose error was less than 1% in the range, 0 ≤ θ ≤ 0.892. Although the error increases, Shin (2021) still applies to the waves in the range, 0.892 ≤ θ < 0.999. When θ > 0.999, wave profile by Shin (2021) does not converge. Each researcher has a different criterion regarding deep water breaking limitations. Chakarabarti (1987) proposes H/Lo = 0.142. According to Dean and Dalrymple, (1984), Michell theory is H/Lo = 0.17, which corresponds to H/L = 0.142. According to Stokes, the breaking criterion is u = c at the crest. In this study, the limitation was checked and the errors according to the required order were also checked. It is also discussed why the other Fourier approximations are unsuitable for deep water waves. For the verification, three waves with a period of 6 s are considered, which were tabulated in Table 1. The following series order was considered: N = 1; N = 3, 6, 10, 12, 13, and 35; M = 180. The RMSE in the DFSBC was calculated and tabulated in Tables 2 and 3. As the required order is increased, the error is decreased. When N ≥13, Newton’s method in Ch. 3 did not converge because e is very large and aN is very small. As a result, this study is available for N ≤12.

Test waves with period T = 6 s

RMSE (%)

RMSE (%) as per Newton’s method step for N = 12

Eq. (35) was coupled with the other equations to calculate the Fourier coefficients, N = M and Eq. (35) is calculated with kη0=12Ni=1N(Xi+Xi+1) in Fentons’ method. Therefore, the required order should be increased in Fenton’s method. Because most errors of Fenton’s method resulted from the numerical integration of the water depth condition, M should be increased to reduce it. Fenton’s method, M = N = 64, is greater than N=13. Therefore, Fenton’s method is unsuitable for deep water waves.

The wave height to wavelength ratio was calculated, as listed in Table 4. The wavelengths were similar to the 5th-order Stokes wave. When θ > 1.028, this study does not converge. More precisely, Newton’s method does not converge because the wave profile at the crest is sharp when θ > 1.028. Therefore, θ = 1.028 is the limitation of this study, which corresponds H/L = 0.137 and is slightly less than 0.142 according to Dean et al., (1984).

Wave height to wavelength ratio, H/L.

The convergent speed decreased when the required order increased. For θ ≤ 0.892, Shin (2021) is acceptable as the first step solution in Newton’s method. For θ > 0.892, the method reported by Shin (2021) is unsuitable for the first step solution because the Newton method does not converge. This problem can be avoided using the sequence of height steps. For waves in the range 0.892 < θ ≤ 0.999, the profile for θ = 0.892 is used as the first step solution, and for waves in the range 0.999 < θ ≤ 1.028, the profile for θ = 0.999 is used as the first step solution. The profiles are shown in Fig. 4; the horizontal velocities are shown in Fig. 5; the vertical velocities are shown in Fig. 6. The relative horizontal velocity at the crest is 0.782 for wave (c), which is less than Stokes’ criteria, 1. The Bernoulli’s constant in Fig. 4 is calculated for wave (c).

Relative horizontal velocities calculated with N = 12 (the ordinate is the dimensionless velocity, u/c)

Relative vertical velocities calculated with N = 12 (the ordinate is the dimensionless velocity, v/c)

## 7. Conclusions

This study aimed to provide a numerical method to calculate the Fourier coefficient, steepness, and reference depth parameter. The numerical procedure is much more simplified than Fenton’s method (Rienecker et al., 1981). The maj or simplifications are as follows:

1. Although the moving coordinate system by Dean was adopted in Fenton’s method, a dimensionless coordinated system was adopted in this study. As a result, all the partial derivatives with respect to wavelength are not required. Some parameters were eliminated.

2. All equations were formulated by tensor analysis. Therefore, numerical equations were much more simplified and had no errors.

3. In the other Fourier approximation, the required order was determined to reduce the error of the numerical integration in the water depth condition because it was solved simultaneously with the Fourier coefficients. On the other hand, the water depth condition was calculated independently. Therefore, the required order was reduced dramatically.

The error was reduced when the required order was increased. Nevertheless, there is a limit to the required order. This study is not applicable when N ≥ 13. Deep water breaking limitation was checked. This study is valid for waves in the range, 0 < θ ≤ 1.028. The limitation 1.028 correspond to H/L = 0.137, which is slightly lower than the Michell theory.

## Notes

The authors declare that they have no conflict of interest.

## References

Chakrabarti S. K. 1987. Hydrodynamics of offshore structures London, UK: Computational Mechanics Publications.
Chaplin J. R. 1979;Developments of stream-function wave theory. Journal of 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.
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.
Det Norske Veritas (DNV). 2010;Environmental conditions and environmental loads. Recommended Practice DNV-RP- C205
Fenton J. D. 1988;The numerical solution of steady water wave problems. Computers & Geosciences 14(3):357–368. https://doi.org/10.1016/0098-3004(88)90066-0.
Fenton J. D. 1990. Nonlinear wave theories. In : Le Mehaute B, Hanes D.M, eds. Ocean Engineering Science, The Sea 9Wiley.
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. 2016;Analytical Approximation in Deep Water Waves. Journal of Advanced Research in Ocean Engineering 2(1):1–11. https://doi.org/10.5574/JAROE.2016.2.1.001.
Shin J. 2019. A Regression Analysis Result for Water Waves on Irrotational Flow over a Horizontal Bed. International Journal of Offshore Polar Engineering. IJOPE 29(4)461–466. https://doi.org/10.17736/ijope.2019.hc17.
Shin J. 2021. A Fourier series approximation for deep-water waves. Journal of Ocean Engineering and Technology 36(2)101–107. https://doi.org/10.26748/KSOE.2021.092.
Stokes G. G. 1847;On the Theory of Oscillatory Waves. Transactions of the Cambridge Philosophical Society 8:441–473.
Tao L, Song H, Chakrabarti S. 2007;Nonlinear progressive waves in water of finite depth -an analytic approximation. Coastal Eng 54(11):825–834.

## Appendices

Appendix A.

#### Solution to KFSBC

The KFSBC is presented in the conventional coordinate system as follows:

(A1) v=ηt+uηx

Because ζ = k(ηη0), we have

(A2) ηt=ωβ{ζk+η0}=ωkζβ=cζβ=cζdβ

As ζ = ζ (β), then we have

(A3) ζβ=dζdβ

And

(A4) ηx=kβ{ζk+η0}=ζdβ=dζdβ

Substituting Eqs. (A2) and (A4) into Eq. (A1), the KFSBC is presented in the dimensionless coordinate system as follows:

(A5) v=cdζdβ+udζdβ

Substituting Eqs. (2) and (3) into Eq. (A5),

(A6) cn=1Nnanenζsinnβ=cdζdβ+{cn=1Nnanenζcosnβ}dζdβ

Because α = ζ on the free surface, then the water particle velocities are presented as follows:

(A7) u(β,ζ)=cn=1Nnanenζcosnβ
and
(A8) v(β,ζ)=cn=1Nnanenζsinnβ

Dividing Eq. (A6) by the celerity and multiplying the result by

(A9) {n=1Nnanenζsinnβ}dβ=dζ+{n=1Nnanenζcosnβ}dζ

The above is presented as follows

(A10) dζ={n=1Nnanenζcosnβ}dζ{n=1Nnanenζsinnβ}dβ

Where the first term on the right-hand side of Eq. (A10) is presented as follows:

(A11) {n=1Nnanenζcosnβ}dζ=dζ{n=1Nanenζcosnβ}dζ

The second term on the right-hand side of Eq. (A10) is presented as follows:

(A12) {n=1Nnanenζsinnβ}dβ=dβ{n=1Nanenζcosnβ}dβ

The right-hand-side of Eq. (A10) is presented as follows:

(A13) dζ=d{n=1Nanenζcosnβ}

Integrating the above equation,

(A14) ζ=n=1Nanenζcosnβ+C1

Using the definition of the reference line (Refer to Fig. 2), we have ζ(±π2)=0, then the integral constant is determined as follows:

(A15) C1=n=1Nancosnπ2

Substituting Eq. (A15) into Eq. (A14) results in Eq. (6).

Appendix B.

#### Solution lor N=1

For N = 1, the wave profile in Eq. (6) is represented as follows:

(A16) ζ=a1eζcosβ

From Eq. (8), the other wave profile is represented as follows:

(A17) ζSθ=a1eζcosβa12e2ζ2+a122

From Eq. (A16), crest height, i.e., the elevation at β = 0 is determined as follows:

(A18) ζc=a1eζc

Therefore, the Fourier coefficient is determined as follows:

(A19) a1=ζceζc

From Eq. (A16), trough depth, i.e., the elevation at β = π is determined as follows:

(A20) ζt=a1eζt

By substituting Eq. (A19) into Eq. (A20) and substituting Eq. (10) into the result, ζt = −ζceS. Substituting it into Eq. (10) allows the crest height to be determined as follows:

(A21) ζc=S1+eS

The trough depth is determined by substituting Eq. (A21) into ζt = −ζceS:

(A22) ζt=SeS1+eS

The wave steepness is determined so that Eq. (A16) and Eq. (A17) meet at phase β = 0. Therefore the following equation can be used:

(A23) S1S2(1+eS){1exp(2S1+eS)}=θ

Eq. (A23) is the dispersion relation. For small amplitude, an approximation to Eq. (A23) is S = θ, which gives the linear dispersion relation ω2 = gk . Applying Newton’s method to Eq. (A23), the steepness was calculated as follows

(A24) S=limmS[m+1]

Where S[1] = θ and

(A25) S[m+1]=S[m]f(S[m])f(S[m])

From Eq. (A23),

(A26) f(S)=2S(1+eS)θ[2(1+eS)S{1exp(2S1+eS)}]

And differentiating Eq. (A26) with regard to S,

(A27) f(S)=2(1+eS)2SeSθ[2eS{1exp(2S1+eS)}2{1+eS+SeS(1+eS)2}exp(2S1+eS)]

The wave profile is calculated with the method in Ch. 3, in which F(β, ζ) = −ζ + a1eζ cos β. This result has an error of less than 1.165% for waves in the range, H/Lo ≤ 0.142, which is less than the error of the 5th order Stokes’ wave theory. Wave (c) was calculated with N = 1, N= 12, and the 5th Stokes theory, and the results were compared in Figs. A1A3. N = 1 gave similar results as the 5th Stokes theory. Comparing Fig. 4 and Fig. A1, the wave profiles calculated by N = 1, and the 5th Stokes theory was less sharp than those calculated by N = 12.

Dimensionless wave profiles calculated with N = 1, N = 12, and 5th Stokes theory (the ordinate is the dimensionless elevation from SWL, )

Relative horizontal velocities calculated with N = 1, N = 12, and 5th Stokes theory (the ordinate is the dimensionless velocity, u/c)

Relative vertical velocities calculated with N = 1, N = 12, and 5th Stokes theory (the ordinate is the dimensionless velocity, v/c)

Appendix C.

#### Solution for N=2

For N=2, Eq. (6) is represented as follows:

(A28) ζ=a1eζcosβ+a2e2ζcos2β+a2

Eq. (8) is represented as follows:

(A29) ζSθ=a1eζcosβ+2a2e2ζcos2βa12e2ζ22a22e2ζ2a1a2e3ζcosβ+2a2+a122+2a22

The coefficients are determined so that Eqs. (A28) and (A29) are equal to each other at phase β = 0 and β = ±π and to satisfy Eq. (10). Eq. (A28) is represented for X1 = ζ(0) as follows:

(A30) X1=a1eX1+a2(e2X1+1)

Eq. (A28) is represented for X2 = ζ(π) as follows:

(A31) X2=a1eX2+a2(e2X2+1)

Solving Eqs. (A30)A31) for a1 and a2, the two coefficients are determined as follows:

(A32) a1=X1(e2X2+1)X2(e2X1+1)eX1(e2X2+1)+eX2(e2X1+1)
and
(A33) a2=X2eX1+X1eX2eX1(e2X2+1)+eX2(e2X1+1)

The error of Eq. (A29) for X1 = ζ (π) is represented as follows:

(A34) E1=X1(X1X2)θ+a1eX1+2a2e2X1a12e2X122a22e2X12a1a2e3X1+2a2+a122+2a22

The error of Eq. (A29) is represented for X2 = ζ (π) as follows:

(A35) E2=X2(X1X2)θ+a1eX2+2a2e2X2a12e2X222a22e2X22a1a2e3X2+2a2+a122+2a22

Therefore, there are five equations, i.e., Eqs. (A32)A35) and Eq. (10) to determine the unknown constants, a1, a2, X1, X2, and S. The set of equations was solved with Newton’s method presented in Ch. 3. The reference depth parameter and wave profile were determined using the method presented in Ch. 4.

## Article information Continued

### Fig. A1

Dimensionless wave profiles calculated with N = 1, N = 12, and 5th Stokes theory (the ordinate is the dimensionless elevation from SWL, )

### Fig. A2

Relative horizontal velocities calculated with N = 1, N = 12, and 5th Stokes theory (the ordinate is the dimensionless velocity, u/c)

### Fig. A3

Relative vertical velocities calculated with N = 1, N = 12, and 5th Stokes theory (the ordinate is the dimensionless velocity, v/c)

### Fig. 1.

Conventional coordinate system for a progressive wave (H: wave height, T: wave period, L : wavelength, η(t, x): free surface elevation from the still water line (SWL)).

### Fig. 2.

Dimensionless coordinate system for a progressive wave.

### Fig. 3.

Flow chart for numerical analysis.

### Fig. 4.

Dimensionless wave profiles calculated with N = 12 (ordinate is the dimensionless elevation from the still water line (SWL), )

### Fig. 5.

Relative horizontal velocities calculated with N = 12 (the ordinate is the dimensionless velocity, u/c)

### Fig. 6.

Relative vertical velocities calculated with N = 12 (the ordinate is the dimensionless velocity, v/c)

### Table 1.

Test waves with period T = 6 s

Wave (a) Wave (b) Wave (c)
θ 0.892 0.999 1.028
H(m) 7.980 8.937 9.196
H/Lo 0.142 0.159 0.164

### Table 2.

RMSE (%)

Theories θ = 0.892 θ = 0.999 θ = 1.028
5th Stokes 2.492 3.950 4.459
N = 1 1.758 2.226 2.360
Shin (2021) 8.26E-2 1.13 Not applicable
N = 6 1.08E-2 9.38E-2 2.15E-1
N = 10 2.78E-4 8.21E-3 2.94E-2
N = 12 6.05E-5 4.85E-3 2.04E-2

### Table 3.

RMSE (%) as per Newton’s method step for N = 12

Wave (a) (b) (c)
Step 1 0.207 2.341) 0.6281)
Step 2 3.27E-3 0.285 3.15E-2
Step 3 6.06E-5 2.01E-2 2.04E-22)
Step 4 6.05E-5 4.85E-3 2.55E-22)

Note that 1En = 1 × 10n

1)

Note that 0.628% is less than 0.207% or 2.34%. The reason is that the first step for wave (c) is the wave profile of wave (b) in step 4, and the first step for wave (b) is the wave profile of wave (a) in step 4, while the first step for wave (a) is the result by Shin (2021) with N = 3.

2)

Note that 2.55E-2 is greater than 2.04E-2. It means that the minimum error is 2.04E-2. As the step is increased, the error oscillates between the two values.

### Table 4.

Wave height to wavelength ratio, H/L.

Theories θ = 0.892 θ = 0.999 θ = 1.028
Airy 0.142 0.159 0.164
5th Stokes 0.123 0.135 0.138
N = 1 0.119 0.129 0.131
Shin (2021) 0.123 0.135 Not applicable
N = 6 0.123 0.134 0.137
N = 10 0.123 0.134 0.137
N = 12 0.123 0.134 0.137