# A Fourier Series Approximation for Deep-water Waves

## Article information

## Abstract

Dean (1965) proposed the use of the root mean square error (RMSE) in the dynamic free surface boundary condition (DFSBC) and kinematic free-surface boundary condition (KFSBC) as an error evaluation criterion for wave theories. There are well known wave theories with RMSE more than 1%, such as Airy theory, Stokes theory, Dean’s stream function theory, Fenton’s theory, and trochodial theory for deep-water waves. However, none of them can be applied for deep-water breaking waves. The purpose of this study is to provide a closed-form solution for deep-water waves with RMSE less than 1% even for breaking waves. This study is based on a previous study (Shin, 2016), and all flow fields were simplified for deep-water waves. For a closed-form solution, all Fourier series coefficients and all related parameters are presented with Newton’s polynomials, which were determined by curve fitting data (Shin, 2016). For verification, a wave in Miche’s limit was calculated, and, the profiles, velocities, and the accelerations were compared with those of 5^{th}-order Stokes theory. The results give greater velocities and acceleration than 5^{th}-order Stokes theory, and the wavelength depends on the wave height. The results satisfy the Laplace equation, bottom boundary condition (BBC), and KFSBC, while Stokes theory satisfies only the Laplace equation and BBC. RMSE in DFSBC less than 7.25×10^{−2}% was obtained. The series order of the proposed method is three, but the series order of 5^{th}-order Stokes theory is five. Nevertheless, this study provides less RMSE than 5^{th}-order Stokes theory. As a result, the method is suitable for offshore structural design.

**Keywords:**Stokes wave theory; Deep-water waves; Fourier series approximation; Newton’s polynomial; Divided differences method

## 1. Introduction

The Brazilian multinational petroleum corporation Petrobras and partners are developing the Buzios oil field, which is approximately 210 km offshore of Brazil. The Floating production storage and offloading (FPSO) is spread-moored in a maximum water depth of 2,030 m and has facilities to receive oil from sub-sea wells. It also has production plant facilities to process fluids, stabilize them, and separate produced water and natural gas, which is re-injected into a dedicated reservoir. Processed liquids are metered, stored in the FPSO cargo storage tanks, and offloaded to export tankers. The design life of the FPSO is 30 years.

The relative water depth is defined as *kh = ω*^{2}*h*/*g*, where *k =* 2*π*/*L* is the wave number, *ω =* 2*π*/*T* is the angular frequency, *L* is the wavelength, *T* is the wave period, *h* is water depth, and *g* is gravity. According to DNV (2010a), it is normally not necessary to investigate wave periods longer than 18 s. Therefore, the relative water depth is greater than 25 for the project. When the relative water depth is greater than 2, deep-water wave theories are applicable (Chakrabarti, 1987; DNV, 2010b; Shin, 2019).

Well-known wave theories include Airy theory, Stokes theory, Dean’s stream function theory, Fenton’s theory, and trochodial (Gerstner) theory for deep-water waves in offshore structure design. Trochodial theory is an exact solution of the Euler equation with vorticity. The first rotational solution was described by Gerstner in 1802 and was independently rediscovered later by Rankine (1863). A mathematical analysis of trochodial theory was performed by Constantin (Henry, 2008). The wavelength is independent of the trochoidal wave’s height, unlike in Stokes’ wave theory and observations. The trajectories of a water particle are closed circles, in contrast with the usual experimental observation of Stokes drift associated with wave motion. Therefore, trochodial theory is of limited use for offshore structure design.

Airy theory, Stokes theory, Dean’s stream function theory, and Fenton’s theory are irrotational wave theories, unlike Trochodial theory. The wavelength is also independent of Airy wave’s height, and Airy theory is applicable for *kH* ≤1/25 (Chakrabarti, 1987), where *H* is the wave height. Therefore, Airy theory is unsuitable for describing waves near the Miche limit (DNV, 2010b); i.e., *H*/*L _{o} =* 0.14, where

*L*is the wavelength calculated by Airy theory.

_{o}The first nonlinear solution was reported by Stokes (1847), who calculated waves in infinite/finite depth. De (1955) and Skjelbreia et al. (1960) obtained the 5^{th}-order solutions for waves in finite depth. Schwartz (1974) and Cokelet (1977) obtained solutions of higher-order expansions. Stokes theory is based on a perturbation method and provides a closed-form solution, but it does not satisfy two boundary conditions on the free surface. The 5^{th}-order theory (De, 1955; Skjelbreia and Hendrickson, 1960) is applicable for all deep waves (Chakrabarti, 1987), but the error of the theory is more than 1% for waves near the Miche limit (DNV, 2010b). Stokes theory gives non-conservative results for offshore structure design.

Based on the use of truncated Fourier expansions for field quantities, Chappelear (1961) and Dean (1965) developed a Fourier approximation method. By choosing expansions to satisfy the Laplace equation and the BBC, the problem is reduced to solving a set of nonlinear equations for each of the Fourier coefficients. Dean’s stream function theory (Dean, 1965) is a purely numerical procedure (DNV, 2010b). It requires order higher than 11 for describing waves near the Miche limit (DNV, 2010b). Nevertheless, its error is more than 1% (DNV, 2010b).

Rienecker and Fenton (1981) adopted the stream function expansion (Dean, 1965) and directly solved for the coefficients by Newton’s method. This method was further simplified by Fenton (1988). The major simplification is that all the partial derivatives are obtained numerically. The most error results from the reason that the wave profile and the Fourier coefficients are simultaneously calculated in Fenton’s theory (Fenton, 1988). To reduce the error, Fenton’s theory requires more series order than 64 (Rienecker and Fenton, 1981). However, since the stream function expansions contain hyperbolic functions, neither of these stream function approaches can be applied for deep-water waves.

Shin (2016) used a Fourier approximation and directly solved the coefficients by Newton’s method, but there are some differences between Shin’s (2016) and Fenton’s theories (Rienecker et al., 1981; Fenton 1988). While Fenton’s theory (Rienecker and Fenton, 1981; Fenton 1988) adopted the stream function (Dean, 1965), Shin (2016) adopted the velocity potential. Fenton’s theory adopted the moving coordinate system introduced by Dean (1965), but Shin (2016) adopted a dimensionless coordinate system. The moving coordinate system requires partial derivatives with regard to the wavelength since the wavelength is coupled to the Fourier coefficients, but the dimensionless coordinate system does not. Therefore, Newton’s method was further simplified (Shin, 2016). While the wave profile and the Fourier coefficients are simultaneously calculated in Fenton’s theory, they are independently calculated by Shin (2016). As a result, the required series order and total error are dramatically reduced (Shin, 2016).

These problems were discussed with Petrobas for the project. It was agreed that Shin’s (2016) wave theory would be applied for calculating wave loads on hard piping riser structures because the theory has less error and provides greater wave loads than 5^{th}-order Stokes theory (Chakrabarti, 1987; De, 1955; Stokes, 1847; Stokes, 1880; DNV, 2010b; Skjelbreia and Hendrickson, 1960). For the application, the results from Shin (2016) were further simplified in this study. The Fourier series is represented for deep-water waves, and all Fourier coefficients and all related parameters are represented with Newton’s polynomials. As a result, this study provides a closed-form solution for deep-water waves. A wave in Miche’s limit was calculated and compared to that of 5^{th}-order Stokes theory. The profiles, velocities, and accelerations were also compared to each other for verification.

## 2. Coordinate Systems

For the description of a progressive water wave, two coordinate systems are adopted. 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 wave propagation, the *y*-axis points upwards, and *t* is time. The fluid domain is bounded by a 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 and two phases of

*β*= ±90°. Therefore, the wave profile is a fixed, periodic, even function in the system. The horizontal axis is the phase,

*β = kx*−

*ω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 dimensionless coordinate system provides several advantages besides the one mentioned in previous section. One independent variable is removed in the coordinate system. Therefore, the Korteweg-de Vries equation and the Benjamin-Bona equation are presented with ordinary differential equations in the coordinate system. The KFSBC is presented with an exact differential equation in the coordinate system.

One parameter is removed in this coordinate system. While deep-water waves are described with two parameters (i.e., wave height and wave period) in a conventional coordinate system or in a moving coordinate system (Dean, 1965), they are described with only one parameter, the linear steepness *θ =* (*ω*^{2}*H*)/*g*, in the dimensionless coordinate system. Furthermore, the linear steepness is in the range of 0≤*θ*<1 since there is a wave-height limitation. Therefore, it is possible to convert the result from Shin (2016) to a closed-form solution.

## 3. Fourier Series Approximation

For deep-water waves, the solution is represented as follows. By using Fourier series approximation, the velocity potential is presented as follows:

*a*is a Fourier coefficient, and

_{n}*N*is the required Fourier series order. The velocity potential satisfies the Laplace equation and the BBC. After differentiating Eq. (1) with respect to

*x*and

*y*, the horizontal water-particle velocity is:

The water-particle accelerations in the horizontal direction and the vertical direction are represented by:

From the KFSBC, the wave profile is calculated as follows:

Using Bernoulli’s equation, the pressure is calculated as follows:

*ρ*is the water density, and

*U*and

*V*are dimensionless horizontal and vertical velocities defined by

*U*

_{0}

*=U*(

*π*/2,0) and

*V*

_{0}

*=V*(

*π*/2,0) are velocities at the phase of

*β*=±

*90°*on the free surface. By applying the DFSBC in Eq. (7), the other wave profile is calculated as follows:

The linear steepness *θ = ω*^{2}*H*/*g* is a known constant for a particular wave. Dimensionless wave height (steepness) is defined by *S = kH* (*H*: wave height), which provides the dispersion relation because the wave number is calculated as *k =S*/*H*.

## 4. Wave Profile

An implicit function *f*(*β*,*ζ*)*=* 0 can be considered as an equation with respect to the dependent variable *ζ* because the independent variable *β* is known. Therefore, using Newton’s method, the explicit representation of the function is presented:

*n*stands for a step of Newton’s method, and

*ζ*

_{n}_{+ 1}is:

*f*(

*β*,

*ζ*) is the error in the KFSBC defined by Eq. (6) as follows:

*f*′(

*β*,

*ζ*) is the first-order partial derivative with respect to the variable

*ζ*. Therefore, we have:

Because |*ζ*|<1 for all waves, the first-step solution *ζ*_{1} is calculated with the power series expansion of *f* in *ζ*, i.e.,
*f*(*β*,*ζ*)*=* 0 for all phases, the power series expansion is a linear equation with respect to *ζ*, i.e., *f*(*β*,0)+*ζf’*(*β*,0)*=* 0. Then we have:

When *θ*≤0.98992, Newton’s method rapidly and absolutely converges to a complete solution. Therefore, *θ*_{max}=0.98992 is the limit in this study. Since the limit is greater than the Miche limit (DNV, 2010b; Shin, 2019; Chakrabarti, 1987), which is 0.28*π*, the proposed method can be used for all deep-water waves.

## 5. Verification

Fourier coefficients were determined so that the two profiles meet at four phases: *β* = 0°, 30°, 90°, and 180°. Therefore, Eq. (6) satisfies the KFSBC for all phases and satisfies the DFSBC at only the four phases. Eq. (8) satisfies the DFSBC for all phases and satisfies the KFSBC at only the four phases, but the difference of the two profiles is very small.

The other wave profile can be calculated by applying the same method presented in section 4 to Eq. (8) or by substituting the result of Eq. (9) into the right side of Eq. (8). A wave with a period of 6 s and height of 7.87 m was calculated by the proposed method and by 5^{th}-order Stokes theory. The wave is in the Miche limit. Wave profiles are compared in Fig. 3. The red solid line is the profile from KFSBC and was calculated with Eq. (6). The green dotted line is the profile from DFSBC and was calculated with Eq. (8). It is extremely difficult to distinguish the two curves. The purple dashed line is the profile calculated with 5^{th}-order Stokes theory. Bernoulli’s constants were also calculated on the free surface. The blue solid line is Bernoulli’s constant from the proposed method, and the blue dashed line is Bernoulli’s constant calculated with 5^{th}-order Stokes theory.

The constants represent error in the DFSBC. If there is no error in the DFSBC, the constant is presented with a horizontal straight line like the blue solid line from Bernoulli’s principle. Therefore, the blue dashed line shows that Stokes theory has some error in the DFSBC. The blue solid line has a total RMSE (Chakrabarti, 1987) of 7.25×10^{−2}%, but the blue dashed line has a total RMSE of 0.37%. In addition to the error, Stokes theory has 1.42% error in KFSBC, unlike the proposed method. Therefore, the total error of Stokes theory is 1.79%.

The wavelength is 64.8 m in the proposed method and 64.52 m in Stokes theory. *L*_{0}*=* 56.21 m is obtained by Airy theory. The wavelength calculated by the proposed method depends on the wave height, unlike Airy theory or trochodial theory.

Water particle velocities on the free surface are presented in Fig. 4. The two solid lines were calculated by the proposed method. The two dashed lines were calculated by 5^{th}-order Stokes theory. The two symmetric curves represent horizontal velocities. There is a big difference in the horizontal velocity under the crest, unlike the wave profiles shown in Fig. 3. As a result, there is more difference in the drag force since the force is proportional to the square of the velocities.

Water particle accelerations on the free surface are presented in Fig. 5. The two solid lines were calculated by the proposed method. The two dashed lines were calculated by 5^{th}-order Stokes theory. The two symmetric curves represent vertical accelerations. There are big differences in the acceleration. As a result, there are big differences in the inertia force since the force is proportional to the acceleration.

Figs. 4–5 show that Stokes theory gives non-conservative results for offshore structure design. The analysis results are summarized in Table 1.

Even though the series order of the proposed method is less than in Stokes theory’s, the error of the proposed method is lower.

## 6. Solution Method and Curve Fitting Results

This study is based on a Fourier approximation. There are differences in numerical methods to determine the Fourier coefficients and the related parameters between Shin’s (2016) method and the other Fourier approximations (Dean, 1965; Chappeler, 1961; Chaplin, 1979; Rienecker and Fenton, 1981; Fenton, 1988). The result from Shin (2016) is summarized for a closed-form solution in the following.

The solution contains 5 unknown constants: 3 Fourier coefficients *a _{n}*, the steepness

*S*, and the reference line −

*kη*

_{0}. Two wave profiles from the KFSBC and the DFSBC were presented with implicit functions in Eqs. (6) and (8). Considering the dimensionless coordinate system, the two wave profiles should be even functions. The Fourier series of a periodic even function is presented as

When *f*(*x*) is an implicit function like Eqs. (6) and (8), the method is impossible. A method to solve the problem is to convert the series to a set of algebraic equations that are obtained when we calculate the series at some phases instead of all phases by replacing the infinite series with a truncated series:
*m=*1, ... , *N*.

Because cos*nx _{m}* and

*f*(

*x*) are known, we have N algebraic equations for calculating the coefficients,

_{m}*a*. Therefore, we have a set of

_{n}*N*equations such that Eqs. (6) and (8) are equal to each other at

*N*phases. In this study,

*N=*3 was considered. Additionally, we have two equations: the wave height condition and water depth condition. As a result, we have 5 equations to determine the unknown constants.

Assuming that *ζ _{c} = ζ*(0

*),*

^{o}*ζ*

_{30}

*= ζ*(30

*), and*

^{o}*ζ*(180

_{t}= ζ*) are known values and applying the idea above to Eq. (6), we have a set of linear algebraic equations with regard to the three coefficients. From the set of equations, the Fourier coefficient*

^{o}*a*

_{2}is determined as follows:

The Fourier coefficient *a*_{1} is determined as follows:

The Fourier coefficient *a*_{3} is determined as follows:

Applying the same idea to Eq. (8) and substituting Eqs. (14), (17), and (18) in Eq. (8), we have a set of three nonlinear equations with regard to *ζ _{c}*,

*ζ*

_{30}, and

*ζ*. From the wave height condition, we have the following:

_{t}We also have an equation from the water depth condition, in which the still water line is the average of the wave profile during a wave period or over a wavelength. Because the horizontal axis of the conventional coordinate system is on the still water line, the reference line is represented as:

Therefore, we have five nonlinear equations, which were solved by Newton’s method by Shin (2016). All the unknown constants are functions of the linear steepness *θ*, which is a known constant for a particular wave in the range of 0≤*θ*<1. They were calculated with a highly dense interval of the linear steepness, and the results were shown by Shin (2016). The constants are represented with Newton’s polynomials in this study. For the curve fitting, the divided differences method (Gerald and Wheatley, 2006) is considered. For *θ*≤0.9, the steepness is presented as follows:

*F*is a Newton basis polynomial defined as

_{i}*θ*<1, the steepness is:

*k =S*/

*H*. For

*θ*≤0.9, the ratio of the trough depth to crest height is presented as follows:

For 0.9<*θ*<1, the ratio of the trough depth to crest height is presented as follows:

For *θ*≤0.9, the ratio of the elevation *ζ*_{30} at a phase of 30° to the crest height is presented as follows:

For 0.9<*θ*<1, the ratio of the elevation *ζ*_{30} at phase of 30° to the crest height is presented as follows:

By using the definition in Eq. (23) and Eq. (19), we have:

Therefore, by substituting Eqs. (21) or (22) and (23) or (24) in Eq. (27), the crest elevation is determined. By using the definition in Eq. (23), we have:

By substituting Eqs. (23) or (24) and (27) in Eq. (28), the trough elevation is determined. By using the definition in Eq. (25), we have:

By substituting Eqs. (25) or (26) and (27) in Eq. (29), the wave elevation at a phase of *β* = 30° is determined. The reference line is numerically calculated by applying Eq. (9) in Eq. (20). The reference line was also calculated with a highly dense interval of the linear steepness. The result is represented with Newton’s polynomials as follows:

For *θ*≤0.9, the parameter *I _{k}* is presented as:

For 0.9<*θ*<1*,* the parameter *I _{k}* is presented as:

The dimensionless elevation is determined by Eq. (9), the reference line *kη*_{0} is determined by Eq. (30), wave number is determined by Eq. (21) or Eq. (22) and then, substituting the three results in the definition of the dimensionless elevation, *ζ = k*(*η*− *η*_{0}), the wave profile is calculated as:

*ζ _{c}*,

*ζ*

_{30}, and

*ζ*are merely parameters to calculate the Fourier coefficients. The wave profile is calculated with Eqs. (9) and (33). This is one of the major differences between Shin’s (2016) theory and Fenton’s theory (Rienecker et al., 1981; Fenton 1988).

_{t}## 7. Conclusions

The purpose of this study was to provide a closed-form solution for deep-water waves with error less than 1%. For this purpose, the result from Shin (2016) was simplified. The major simplification is that all Fourier coefficients and all parameters are presented with Newton’s polynomials for a closed-form solution and engineering application. Therefore, a numerical procedure is not necessary for calculating the coefficients and the parameters any more.

The required Fourier series order is *N =* 3 in this study, while more than 11 is required by Dean’s stream function theory (DNV, 2010b), and more than *N =* 64 is required by Fenton’s theory (Rienecker and Fenton, 1981). The series order of 5^{th}-order Stokes theory is 5. Therefore, the required series order is dramatically reduced. The result has less error, is simpler, and gives greater water particle velocities and water particle accelerations than the 5th-order Stokes wave theory. Therefore, the Stokes wave theory is not conservative, and the proposed method is more suitable for offshore structure design.

## Acknowledgements

We would like to thank Petrobras for applying the proposed method in the FPSO project.

## Notes

**Conflict of Interest**

No potential conflict of interest relevant to this article was reported.