# Numerical Bichromatic Wave Generation Using Designed Mass Source Function

Authors: Min-yi Chen; Hong-sheng Zhang; En-xian Zhou; Da-li Xu
Journal: Journal of Mathematics
 (2021)

---

## Abstract

A mass source wave-maker method is generalized as the two-wave-source wave-maker method to generate bichromatic waves in the numerical model, whose governing equations are Navier–Stokes equations with the continuity equation. The Fluent software is taken as the calculation platform. In the numerical model, the waves at both the left and right ends of the numerical wave flume are absorbed with the momentum sources added in Navier–Stokes equations. The numerical simulation of bichromatic waves propagation with different frequencies in uniform deep, intermediate, and shallow water has been conducted. The numerical solutions are compared with the theoretical solutions obtained on the basis of Stokes waves theory. The frequency spectrum analyses of the results are conducted and discussed, and the differences between the weakly nonlinear theoretical solutions and the fully nonlinear numerical results are investigated in detail. It is found that the numerical model can effectively simulate the nonlinear effect of bichromatic waves in water with different depths, and the theoretical solutions only adapt the deep and intermediate water. The results indicate that the present numerical model is valuable in the aspect of practical application.


---

## Body

## 1. Introduction

Nonlinear problem is one of the important subjects of gravity surface waves. Phillips [1] made a pioneering work for the study of wave interaction dynamics in the random gravity wave fields with finite amplitude by analyzing the third-order interaction of waves. After that, theoretical studies about the nonlinear wave-wave interaction have been conducted [2–4]. The studies of the nonlinear effect and evolution of bichromatic waves are included in the work of Madsen and Fuhrman [5], Lin at al. [6], Halfiani and Ramli [7], and others. Hong [8] derived a fourth-order approximate theoretical solution of the nonlinear interaction of surface gravity waves in water with uniform depth based on the Stokes finite amplitude wave theory by the application of the perturbation method.With the development of computer technology, the numerical simulation has gradually become an important approach to study the nonlinear interaction problem, and many researchers have provided different numerical models to study this problem. With the nonlinear Boussinesq-type equations used as the governing equations, different numerical models are provided by Madsen and Sørensen [9] and Zhang et al.[10] to solve the nonlinear interaction problems. Ohyama and Nadaoka [11] and Ohyama at al. [12] developed a nonreflective numerical wave flume to analyze the nonlinear wave fields, and the numerical model was then applied to study the decomposition phenomenon of waves when they passed through a submerged shelf without breaking, and subsequently compared the effectiveness of three different numerical models, which are based on different nonlinear water wave theories, respectively. Different nonhydrostatic models were proposed and developed by Young and Wu [13], Dong et al.[14], and Ai et al.[15] to simulate the propagation of nonlinear waves trains.The source wave-maker method can be divided into two types according to the different governing equations. Adding the source function to the continuity equation is called the mass source wave-maker, and adding the source function to the momentum equation is called the momentum source wave-maker. There are different approaches to generating bichromatic waves in numerical wave flume (Lin and Liu [16] and others). However, some specific assumptions are needed if the bichromatic waves are generated directly at a fixed location. In this study, to simulate a more realistic wave field, the traditional mass source wave-maker method is revised as the two wave sources wave-maker method to generate the bichromatic waves in the numerical wave flume whose governing equations are the incompressible N–S. equations with the continuity equation. The VOF method is used in the numerical model, where the Fluent software is taken as the calculation platform. The present numerical model is completely nonlinear since its governing equations are incompressible N–S. equations with the continuity equation. The present numerical model is used to study the interactions of bichromatic waves, the numerical results are compared with theoretical solutions, and the differences between the numerical results and the theoretical results are analyzed in detail.


## 2. Numerical Model

### 2.1. The Governing Equations

The numerical wave flume was built with the Fluent software in this study. The flow motion of incompressible fluid is described by the continuity equation and Navier–Stokes equations:(1)∂u∂x+∂w∂z=Sm,ρ∂u∂t+u∂u∂x+w∂u∂z=−∂p∂x+μ∂2u∂x2+∂2u∂z2+Sx,ρ∂w∂t+u∂w∂x+w∂w∂z=ρg−∂p∂z+μ∂2w∂x2+∂2w∂z2+Sz,x and z form the rectangular coordinate-system with z measured vertically upwards from still water level; t is the time; ρ is the fluid density; u and w are the velocity components in the x- and z-directions, respectively; p is the pressure; μ is the hydrodynamic viscosity coefficient (μ = 0 in this study); g is the acceleration of gravity; Sm is the additional mass source term; Sx and Sz are the additional momentum source terms in the x- and z-directions, respectively.


### 2.2. The Formula of Mass Source Function

Lin and Liu [16] proposed an approach to generating waves by using designed mass source functions for the mass conservation equation in the internal flow region. The mass source term function Sm is expressed as follows:(2)Smt=cHAΩcoskx−σt,where H is the wave height, k is the wave number, σ is wave circular frequency, c is the wave velocity, and AΩ is the total area of the grid covered by the wave-maker region.Two mass sources will be set in the wave flume, and the approach to generating waves by using designed mass source function will thus be generalized. The improved approach will be applied to study the wave-wave interaction in the shallow or deep water. The set of the numerical wave flume is shown in Figure1, where xi (i = 1, 2) is the horizontal coordinate of the center point of wave-maker region Ω. The horizontal coordinate of the center point for the first source, x1, equates 0, and that for the second source is taken as x2. The functions for the two mass sources are expressed as follows:(3)Sm1=cHAΩcosσt,Sm2t=cHAΩcoskx2−σt,where Sm1 and Sm2 are the mass source functions of S1 and S2, respectively.Figure 1
Sketch of the numerical wave flume.The source region is set inside the computational domain and the source with finite volume is used. Since the two-dimensional case is discussed only in this study, the shape of source region is set as a rectangle, which is equivalent to a source or a sink. The VOF method is used to track the position of the free surface.


### 2.3. Momentum Source Functions for Wave Absorption

Wave absorption is required at the left and right boundaries of the model when numerical simulation of wave propagation is carried out. Artificial wave absorption domains are set at the left and right ends of numerical wave flume, respectively. The wave absorption effect is related to the width of the wave absorption,Ls. In this study, Ls is chosen as the maximum wavelength of the two incident waves. Wave absorption is realized by adding source functions to the momentum equations. The functions for the source terms are expressed as follows:(4)Sx=−ρξxu,Sz=−ρξxw,where(5)ξx=αx−xslxsr−xsl.α is the damping coefficient. And after modelling test, it is taken as 6.0 in this study; xsl is the horizontal coordinate of the head of the wave absorption region, while xsr is the far end of the wave absorption region. A similar expression to Equation (5) is adopted at the front end of the numerical wave flume.


### 2.4. Boundary Conditions

Because the Fluent software is used to establish the numerical wave flume, the secondary development of Fluent software is required in the wave-maker and wave absorption regions, and the source terms (Sm, Sx, and Sz) in the governing equations are defined, respectively. The left, right, and bottom boundaries of the wave flume are all set as the wall boundaries, and the top boundary is set as the pressure inlet.


## 2.1. The Governing Equations

The numerical wave flume was built with the Fluent software in this study. The flow motion of incompressible fluid is described by the continuity equation and Navier–Stokes equations:(1)∂u∂x+∂w∂z=Sm,ρ∂u∂t+u∂u∂x+w∂u∂z=−∂p∂x+μ∂2u∂x2+∂2u∂z2+Sx,ρ∂w∂t+u∂w∂x+w∂w∂z=ρg−∂p∂z+μ∂2w∂x2+∂2w∂z2+Sz,x and z form the rectangular coordinate-system with z measured vertically upwards from still water level; t is the time; ρ is the fluid density; u and w are the velocity components in the x- and z-directions, respectively; p is the pressure; μ is the hydrodynamic viscosity coefficient (μ = 0 in this study); g is the acceleration of gravity; Sm is the additional mass source term; Sx and Sz are the additional momentum source terms in the x- and z-directions, respectively.


## 2.2. The Formula of Mass Source Function

Lin and Liu [16] proposed an approach to generating waves by using designed mass source functions for the mass conservation equation in the internal flow region. The mass source term function Sm is expressed as follows:(2)Smt=cHAΩcoskx−σt,where H is the wave height, k is the wave number, σ is wave circular frequency, c is the wave velocity, and AΩ is the total area of the grid covered by the wave-maker region.Two mass sources will be set in the wave flume, and the approach to generating waves by using designed mass source function will thus be generalized. The improved approach will be applied to study the wave-wave interaction in the shallow or deep water. The set of the numerical wave flume is shown in Figure1, where xi (i = 1, 2) is the horizontal coordinate of the center point of wave-maker region Ω. The horizontal coordinate of the center point for the first source, x1, equates 0, and that for the second source is taken as x2. The functions for the two mass sources are expressed as follows:(3)Sm1=cHAΩcosσt,Sm2t=cHAΩcoskx2−σt,where Sm1 and Sm2 are the mass source functions of S1 and S2, respectively.Figure 1
Sketch of the numerical wave flume.The source region is set inside the computational domain and the source with finite volume is used. Since the two-dimensional case is discussed only in this study, the shape of source region is set as a rectangle, which is equivalent to a source or a sink. The VOF method is used to track the position of the free surface.


## 2.3. Momentum Source Functions for Wave Absorption

Wave absorption is required at the left and right boundaries of the model when numerical simulation of wave propagation is carried out. Artificial wave absorption domains are set at the left and right ends of numerical wave flume, respectively. The wave absorption effect is related to the width of the wave absorption,Ls. In this study, Ls is chosen as the maximum wavelength of the two incident waves. Wave absorption is realized by adding source functions to the momentum equations. The functions for the source terms are expressed as follows:(4)Sx=−ρξxu,Sz=−ρξxw,where(5)ξx=αx−xslxsr−xsl.α is the damping coefficient. And after modelling test, it is taken as 6.0 in this study; xsl is the horizontal coordinate of the head of the wave absorption region, while xsr is the far end of the wave absorption region. A similar expression to Equation (5) is adopted at the front end of the numerical wave flume.


## 2.4. Boundary Conditions

Because the Fluent software is used to establish the numerical wave flume, the secondary development of Fluent software is required in the wave-maker and wave absorption regions, and the source terms (Sm, Sx, and Sz) in the governing equations are defined, respectively. The left, right, and bottom boundaries of the wave flume are all set as the wall boundaries, and the top boundary is set as the pressure inlet.


## 3. Numerical Simulation of Bichromatic Waves on Constant Water Depth

### 3.1. Setup of Numerical Model

Using the perturbation method, Hong [8] derived the theoretical solution for the nonlinear interaction of surface gravity waves. In this study, we compared the numerical results with the second-order theoretical solutions, since the amplitudes of higher-order harmonics are very small.Four combinations of incident bichromatic waves with different periods have been calculated in this study. The amplitudes of incident waves were all set as 0.01 m. The parameters of incident waves are shown in Table1, where μ (μ = h/L) denotes the ratio of water depth to the wave length. The effective length of the computation domain of the numerical water flume was set as 20.0 times the maximum wavelength of the two incident waves. The water depths of the case C1 and C2 were both set as 0.45 m. The sizes of the two wave sources were the same, and the width and height were set as 0.2 m and 0.44 m, respectively. The center of the wave source region was located at almost 0.5 times the water depth below the still water surface. The left end of S1 was 2.0 times the maximum wavelength away from the front wave absorption region, S2 was located at the right side of S1, and the horizontal distance between the centers of the two wave source regions was 2.0 m. For cases C1 and C2, the entire computation domain was uniformly discretized by a grid system with ∆x = 0.05 m in the horizontal direction and ∆z = 0.01 m in the vertical direction. For cases D1 and D2, the water depths were both set as 4.0 m and the width and height of the two wave source regions were set as 0.1 m and 0.74 m, respectively. The distance between the top of the wave source and the still water surface was the same as the wave height. The grid size in the horizontal direction, ∆x, was taken as 0.05 m near the wave source region, and taken as 0.08 m in the rest of the computation domain; ∆x was nonuniform in the sponge-layer, which was fine near the front end of the absorption and coarse near the trailing end. The grid size in the vertical direction, ∆z, was taken as 0.025 m near the free surface, and taken as 0.18 m near the bottom. The other relevant parameters were the same as those of the cases C1 and C2, and the time steps, ∆t, for all cases were set as 0.01 s. The frequencies of the two incident waves were f1 and f2, respectively. The nonlinear interactions between waves resulted in the sum frequency components, difference frequency components, and higher harmonic components. The different frequency components are shown in Table 2.Table 1
Table of the parameters of different combinations of incident waves.
CasesIncident wave1Incident wave2T1/sL1/mμ1T2/sL2/mμ2C123.880.122.55.000.09C223.880.124.28.670.05D126.240.642.59.650.41D226.240.644.222.310.18Table 2
Table of main components of frequency spectrum analysis for different cases (HZ).
f1f22f12f2f1 + f2f1 + 2f22f1 + f22f1 + 2f2f1 − f22f2 − f12f1 − f22f1 − 2f2C10.500.401.000.800.901.301.401.800.100.300.600.20C20.500.241.000.480.740.981.241.480.260.020.760.52D10.500.401.000.800.901.301.401.800.100.300.600.20D20.500.241.000.480.740.981.241.480.260.020.760.52


### 3.2. Calculation Results

Figures2 and 3 show the comparisons of wave profiles between the numerical solutions and the theoretical solutions for cases C1 and C2, respectively. Case C1 considers the bichromatic wave combination of two transition water waves, and case C2 considers the bichromatic wave combination of a transition water wave and a shallow water wave. It can be seen from Figures 2 and 3 that after the two wave trains are superimposed, they still show obvious periodicity; the amplitude, however, is changed. Comparing Figure 2(b) with 2(b), it can be observed that the superposed waves propagate forward steadily. It can be seen from Figure 3 that the superimposed wave crest of case C2 becomes steeper, and a second wave crest appears. Comparing Figure 2 with 3, the calculated wave height in Figure 3 is smaller than that in Figure 2, indicating that the nonlinearity for case C2 is stronger than that for case C1. Due to the difference of wavelengths between the two incident waves of case C2 being large, the effect of the relatively short wave (the wave with shorter wavelength of the two incident waves) is equivalent to adding a disturbance to the relatively long wave (the wave with longer wavelength of the incident waves). The wave profile shown in Figure 3 also confirms this. The wave heights of the numerical solutions are slightly smaller than those of the theoretical numerical solutions. This is because the wave energy can be more effectively redistributed in the numerical solutions between different orders of harmonics due to the nonlinearity.Figure 2
Comparisons of the wave profiles between the calculation results (solid line) and theoretical results (dashed line) at different gauge locations for case C1.
(a)(b)Figure 3
Comparisons of the wave profiles between the calculation results (solid line) and theoretical results (dashed line) at different gauge locations for case C2.
(a)(b)Spectrum analyses of the calculation results are conducted to analyze the nonlinear interactions between different wave components. The frequency components are obtained by Fourier transformation, the time domain data of the wave profile at different locations. Figures4 and 5 show the spatial variations of the amplitudes of each frequency component for cases C1 and C2, respectively.Figure 4
Comparisons of the spatial variation of amplitudes of each frequency component between calculation and theoretical results for case C1., ,  calculation solutions; , ,  theoretical solutions.
(a)(b)(c)(d)(e)Figure 5
Comparisons of the spatial variation of amplitudes of each frequency component between calculation and theoretical results for case C2., ,  calculation solutions; , ,  theoretical solutions.
(a)(b)(c)(d)(e)It can be seen from Figure4(a) that the slight fluctuations appear in the numerical solutions of the two primary frequencies along the wave flume. It should be noted that the calculated amplitude of the primary frequency component f1 is slightly smaller than that of f2 along the wave flume, which indicates that the relatively long wave is slightly affected by the nonlinear effect, and thus it further verifies that the relatively long wave in Figure 2 plays a dominant role in the wave profile in the process of nonlinear interactions. It is shown in Figure 4(b) that the numerical solutions of the difference frequency components f1 − f2 oscillate along the wave flume, and the amplitudes of numerical solutions are greater than those of the theoretical solutions. In Figure 4(c), the calculated second harmonic component 2f2 exhibits periodic spatial variation as follows:(6)λ2=2πk2−2k1.k1 and k2 satisfy the following linear dispersion relations, respectively:(7)k2tanhk2h=4σ2g,k1tanhk1h=σ2g.Equation (6) gives the expression of the beat length λ2 of the second harmonic component 2fi (i = 1, 2), where k1 and k2 are the wave numbers of the primary wave component and second harmonic component, respectively [17]. The beat length of the second harmonic component 2f2, λ2 14.4 m, obtained from Equation (6) is consistent with the result shown in Figure 4(c). The second harmonic component 2f1 of another primary wave component f1 oscillates around the amplitude of the theoretical solution in the process of propagation. The amplitude of the sum frequency component f1 + f2 of the numerical solutions is slightly larger than that of 2f1, and smaller than that of 2f2, which is consistent with the order of magnitude between theoretical solutions. Comparing Figures 4(c) and 4(d), it is found that the calculated difference of frequency components (2f1 − f2, f1 − 2f2, and 2f1 − 2f2) between the primary frequencies and the second harmonics and the second harmonics of each primary wave become obvious in the second half of the wave flume, which satisfies the relationship of wave-energy conservation. During the process of wave propagation, the nonlinear effect between different components continues to deepen, and the wave energy transfers from low-frequency components to higher-frequency components. It can be seen from Figure 4(e) that the calculated sum frequency components (2f1 + f2, f1 + 2f2 and 2f1 + 2f2) between the primary frequencies and the second harmonics and the second harmonics of each primary wave are all extremely small.It is shown in Figure5(a) that the calculated magnitude of the primary frequency component f1 is reduced while that of f2 remains unchanged generally. Comparing Figure 4(a) with Figure 5(a), it is found that, when the wavelength of incident wave 1 remains unchanged, a greater wavelength of incident wave 2 results in a greater reduction of the amplitude of the primary frequency component f1, which indicates that a greater difference of wavelength between two incident waves results in a stronger nonlinear effect. As shown in Figure 5(b), the numerical solution of the difference frequency component f1 − f2 appears spatial modulations and is obviously greater than the theoretical solution. This is because the incident wave 2 of case C2 is a shallow water wave, and the component f1 − f2 is obviously out of the application scope of Stokes wave theory. Figure 5(c) shows that the sum frequency component f1 + f2 and the second harmonic component 2f2both exhibit periodic variation along the wave flume. The beat length of the second harmonic component calculated with Equation (6), λ2, equates 69.813 m, and it is almost consistent with the result shown in Figure 5(c). Comparing Figures 4(d) and 4(e) with Figures 5(d) and 5(e), it can be found that, for case C2, the amplitudes of different frequency components and sum frequency components between the primary frequencies and the second harmonics, and the second harmonics of each primary frequency are obviously greater than the corresponding amplitudes of case C1, which further illustrates that the nonlinear effect between waves increases as the wavelength difference increases.Cases D1 and D2 consider the combinations of deep-water waves and intermediate water waves with different dispersion characteristics and are used to study the nonlinear effects between different types of bichromatic waves. The time series of wave profiles of calculated and theoretical solutions at different locations are selected for comparative analyses, and the results are shown in Figures6 and 7, respectively.Figure 6
Comparisons of the wave profiles between the calculation results (solid line) and theoretical results (dashed line) at different gauge locations for case D1.
(a)(b)Figure 7
Comparisons of the wave profiles between the calculation results (solid line) and theoretical results (dashed line) at different gauge locations for case D2.
(a)(b)Figure6 shows that the calculation solutions of case D1 are in good agreement with the theoretical solutions, and the time series of wave surfaces also coincides basically. This is because μ2 of this case is 0.41, which is close to the deep-water wave, and the two incident waves can thus be both considered as the deep-water waves whose nonlinear interactions are weak and correspond to the application scope of Stokes wave theory. The wave field of case D2 is similar to that of case C2, μ2 of case D2 is 0.18, and the wavelength of the incident wave 2 is long. The nonlinear interactions between two wave trains of case D2 are strong. The theoretical solutions cannot correctly describe the nonlinear wave field of case D2. This is why the wave heights of the theoretical solutions are different from those of the numerical solutions in Figure 7. Figures 6 and 7 show that the numerical model established in this study can also effectively simulate the wave propagation of bichromatic waves in deep-water or in medium-water depth.Similarly, to analyze the wave field due to the nonlinear interaction of two waves, the spatial variations of the amplitudes of each frequency component at different locations are shown in Figures8 and 9. The calculated sum of frequency components (2f1 + f2, f1 + 2f2 and 2f1 + 2f2) almost equates 0, and they would not be shown in Figures 8 and 9. In view of the fact that case D1 is a combination of a deep-water wave and a transition water wave, Figure 8(a) shows that the numerical solutions of the two primary frequency components are both reduced and that the amplitude of the primary frequency f1 is smaller than that of f2. It is shown that the relatively short-wave amplitude is reduced more than the relatively long-wave amplitude, which indicates that the relatively short wave is more sensitive by the nonlinear interactions.Figure 8
Comparisons of the spatial variation of amplitudes of each frequency component between calculation and theoretical results for case D1., ,  calculation solutions; , ,  theoretical solutions.
(a)(b)(c)(d)Figure 9
Comparisons of the spatial variation of amplitudes of each frequency component between calculation and theoretical results for case D2., ,  calculation solutions; , ,  theoretical solutions.
(a)(b)(c)(d)Comparing Figures8(b)–8(d), it can be found that the magnitudes of the numerical solution and theoretical solution of each harmonic component are both small, which indicates that the nonlinear effect between the two wave trains of case D1 is weak. This is consistent with the results shown in Figure 6. It can be seen from Figure 8(c) that the sum frequency component (f1 + f2) and the second harmonic components (2f2 and 2f1) can be distinguished and the slight fluctuations appear near the wave source region; however, the fluctuations begin to disappear as the distance of wave propagation increases.Case D2 considers the combination of a deep-water wave and a transition water wave. It is shown in Figure9(a) that the amplitude of the relatively short wave is reduced while the amplitude of the relatively long wave is consistent with the theoretical solution, which indicates that the nonlinearity weakens the energy of the relatively short wave. It can be seen from Figures 9(b) and 9(c) that the difference frequency component (f1 − f2), the sum frequency component (f1 + f2), and the second harmonic components (2f2 and 2f1) present obvious magnitudes. The difference frequency component (f1 − f2) seems to propagate with a certain period. The second harmonic component 2f2 is spatially modulated, and it changes obviously periodically. Under the disturbance of short wave, the amplitude of the second harmonic component of the long-wave changes continuously, showing that the relatively long wave, wave 2, is more affected by wave 1. In view of the difference frequency components between the primary frequencies and the second harmonics, the second harmonics of each primary frequency (2f1 − f2, f1 − 2f2 and 2f1 − 2f2) are distinguishable during the propagation process.In general, the numerical results of the sum frequency components and difference frequency components between the higher harmonics and the higher harmonic and the primary frequency are greater than those of the theoretical solutions. This is because the theoretical solution is derived from the weakly nonlinear Stokes wave theory, while the numerical model established in this study is fully nonlinear. The fact that the obtained numerical solutions can better describe the nonlinear effect between wave components than the theoretical solutions further illustrates that the numerical model of wave-maker by setting two wave source regions can effectively simulate the nonlinear effect between two incident waves.


## 3.1. Setup of Numerical Model

Using the perturbation method, Hong [8] derived the theoretical solution for the nonlinear interaction of surface gravity waves. In this study, we compared the numerical results with the second-order theoretical solutions, since the amplitudes of higher-order harmonics are very small.Four combinations of incident bichromatic waves with different periods have been calculated in this study. The amplitudes of incident waves were all set as 0.01 m. The parameters of incident waves are shown in Table1, where μ (μ = h/L) denotes the ratio of water depth to the wave length. The effective length of the computation domain of the numerical water flume was set as 20.0 times the maximum wavelength of the two incident waves. The water depths of the case C1 and C2 were both set as 0.45 m. The sizes of the two wave sources were the same, and the width and height were set as 0.2 m and 0.44 m, respectively. The center of the wave source region was located at almost 0.5 times the water depth below the still water surface. The left end of S1 was 2.0 times the maximum wavelength away from the front wave absorption region, S2 was located at the right side of S1, and the horizontal distance between the centers of the two wave source regions was 2.0 m. For cases C1 and C2, the entire computation domain was uniformly discretized by a grid system with ∆x = 0.05 m in the horizontal direction and ∆z = 0.01 m in the vertical direction. For cases D1 and D2, the water depths were both set as 4.0 m and the width and height of the two wave source regions were set as 0.1 m and 0.74 m, respectively. The distance between the top of the wave source and the still water surface was the same as the wave height. The grid size in the horizontal direction, ∆x, was taken as 0.05 m near the wave source region, and taken as 0.08 m in the rest of the computation domain; ∆x was nonuniform in the sponge-layer, which was fine near the front end of the absorption and coarse near the trailing end. The grid size in the vertical direction, ∆z, was taken as 0.025 m near the free surface, and taken as 0.18 m near the bottom. The other relevant parameters were the same as those of the cases C1 and C2, and the time steps, ∆t, for all cases were set as 0.01 s. The frequencies of the two incident waves were f1 and f2, respectively. The nonlinear interactions between waves resulted in the sum frequency components, difference frequency components, and higher harmonic components. The different frequency components are shown in Table 2.Table 1
Table of the parameters of different combinations of incident waves.
CasesIncident wave1Incident wave2T1/sL1/mμ1T2/sL2/mμ2C123.880.122.55.000.09C223.880.124.28.670.05D126.240.642.59.650.41D226.240.644.222.310.18Table 2
Table of main components of frequency spectrum analysis for different cases (HZ).
f1f22f12f2f1 + f2f1 + 2f22f1 + f22f1 + 2f2f1 − f22f2 − f12f1 − f22f1 − 2f2C10.500.401.000.800.901.301.401.800.100.300.600.20C20.500.241.000.480.740.981.241.480.260.020.760.52D10.500.401.000.800.901.301.401.800.100.300.600.20D20.500.241.000.480.740.981.241.480.260.020.760.52


## 3.2. Calculation Results

Figures2 and 3 show the comparisons of wave profiles between the numerical solutions and the theoretical solutions for cases C1 and C2, respectively. Case C1 considers the bichromatic wave combination of two transition water waves, and case C2 considers the bichromatic wave combination of a transition water wave and a shallow water wave. It can be seen from Figures 2 and 3 that after the two wave trains are superimposed, they still show obvious periodicity; the amplitude, however, is changed. Comparing Figure 2(b) with 2(b), it can be observed that the superposed waves propagate forward steadily. It can be seen from Figure 3 that the superimposed wave crest of case C2 becomes steeper, and a second wave crest appears. Comparing Figure 2 with 3, the calculated wave height in Figure 3 is smaller than that in Figure 2, indicating that the nonlinearity for case C2 is stronger than that for case C1. Due to the difference of wavelengths between the two incident waves of case C2 being large, the effect of the relatively short wave (the wave with shorter wavelength of the two incident waves) is equivalent to adding a disturbance to the relatively long wave (the wave with longer wavelength of the incident waves). The wave profile shown in Figure 3 also confirms this. The wave heights of the numerical solutions are slightly smaller than those of the theoretical numerical solutions. This is because the wave energy can be more effectively redistributed in the numerical solutions between different orders of harmonics due to the nonlinearity.Figure 2
Comparisons of the wave profiles between the calculation results (solid line) and theoretical results (dashed line) at different gauge locations for case C1.
(a)(b)Figure 3
Comparisons of the wave profiles between the calculation results (solid line) and theoretical results (dashed line) at different gauge locations for case C2.
(a)(b)Spectrum analyses of the calculation results are conducted to analyze the nonlinear interactions between different wave components. The frequency components are obtained by Fourier transformation, the time domain data of the wave profile at different locations. Figures4 and 5 show the spatial variations of the amplitudes of each frequency component for cases C1 and C2, respectively.Figure 4
Comparisons of the spatial variation of amplitudes of each frequency component between calculation and theoretical results for case C1., ,  calculation solutions; , ,  theoretical solutions.
(a)(b)(c)(d)(e)Figure 5
Comparisons of the spatial variation of amplitudes of each frequency component between calculation and theoretical results for case C2., ,  calculation solutions; , ,  theoretical solutions.
(a)(b)(c)(d)(e)It can be seen from Figure4(a) that the slight fluctuations appear in the numerical solutions of the two primary frequencies along the wave flume. It should be noted that the calculated amplitude of the primary frequency component f1 is slightly smaller than that of f2 along the wave flume, which indicates that the relatively long wave is slightly affected by the nonlinear effect, and thus it further verifies that the relatively long wave in Figure 2 plays a dominant role in the wave profile in the process of nonlinear interactions. It is shown in Figure 4(b) that the numerical solutions of the difference frequency components f1 − f2 oscillate along the wave flume, and the amplitudes of numerical solutions are greater than those of the theoretical solutions. In Figure 4(c), the calculated second harmonic component 2f2 exhibits periodic spatial variation as follows:(6)λ2=2πk2−2k1.k1 and k2 satisfy the following linear dispersion relations, respectively:(7)k2tanhk2h=4σ2g,k1tanhk1h=σ2g.Equation (6) gives the expression of the beat length λ2 of the second harmonic component 2fi (i = 1, 2), where k1 and k2 are the wave numbers of the primary wave component and second harmonic component, respectively [17]. The beat length of the second harmonic component 2f2, λ2 14.4 m, obtained from Equation (6) is consistent with the result shown in Figure 4(c). The second harmonic component 2f1 of another primary wave component f1 oscillates around the amplitude of the theoretical solution in the process of propagation. The amplitude of the sum frequency component f1 + f2 of the numerical solutions is slightly larger than that of 2f1, and smaller than that of 2f2, which is consistent with the order of magnitude between theoretical solutions. Comparing Figures 4(c) and 4(d), it is found that the calculated difference of frequency components (2f1 − f2, f1 − 2f2, and 2f1 − 2f2) between the primary frequencies and the second harmonics and the second harmonics of each primary wave become obvious in the second half of the wave flume, which satisfies the relationship of wave-energy conservation. During the process of wave propagation, the nonlinear effect between different components continues to deepen, and the wave energy transfers from low-frequency components to higher-frequency components. It can be seen from Figure 4(e) that the calculated sum frequency components (2f1 + f2, f1 + 2f2 and 2f1 + 2f2) between the primary frequencies and the second harmonics and the second harmonics of each primary wave are all extremely small.It is shown in Figure5(a) that the calculated magnitude of the primary frequency component f1 is reduced while that of f2 remains unchanged generally. Comparing Figure 4(a) with Figure 5(a), it is found that, when the wavelength of incident wave 1 remains unchanged, a greater wavelength of incident wave 2 results in a greater reduction of the amplitude of the primary frequency component f1, which indicates that a greater difference of wavelength between two incident waves results in a stronger nonlinear effect. As shown in Figure 5(b), the numerical solution of the difference frequency component f1 − f2 appears spatial modulations and is obviously greater than the theoretical solution. This is because the incident wave 2 of case C2 is a shallow water wave, and the component f1 − f2 is obviously out of the application scope of Stokes wave theory. Figure 5(c) shows that the sum frequency component f1 + f2 and the second harmonic component 2f2both exhibit periodic variation along the wave flume. The beat length of the second harmonic component calculated with Equation (6), λ2, equates 69.813 m, and it is almost consistent with the result shown in Figure 5(c). Comparing Figures 4(d) and 4(e) with Figures 5(d) and 5(e), it can be found that, for case C2, the amplitudes of different frequency components and sum frequency components between the primary frequencies and the second harmonics, and the second harmonics of each primary frequency are obviously greater than the corresponding amplitudes of case C1, which further illustrates that the nonlinear effect between waves increases as the wavelength difference increases.Cases D1 and D2 consider the combinations of deep-water waves and intermediate water waves with different dispersion characteristics and are used to study the nonlinear effects between different types of bichromatic waves. The time series of wave profiles of calculated and theoretical solutions at different locations are selected for comparative analyses, and the results are shown in Figures6 and 7, respectively.Figure 6
Comparisons of the wave profiles between the calculation results (solid line) and theoretical results (dashed line) at different gauge locations for case D1.
(a)(b)Figure 7
Comparisons of the wave profiles between the calculation results (solid line) and theoretical results (dashed line) at different gauge locations for case D2.
(a)(b)Figure6 shows that the calculation solutions of case D1 are in good agreement with the theoretical solutions, and the time series of wave surfaces also coincides basically. This is because μ2 of this case is 0.41, which is close to the deep-water wave, and the two incident waves can thus be both considered as the deep-water waves whose nonlinear interactions are weak and correspond to the application scope of Stokes wave theory. The wave field of case D2 is similar to that of case C2, μ2 of case D2 is 0.18, and the wavelength of the incident wave 2 is long. The nonlinear interactions between two wave trains of case D2 are strong. The theoretical solutions cannot correctly describe the nonlinear wave field of case D2. This is why the wave heights of the theoretical solutions are different from those of the numerical solutions in Figure 7. Figures 6 and 7 show that the numerical model established in this study can also effectively simulate the wave propagation of bichromatic waves in deep-water or in medium-water depth.Similarly, to analyze the wave field due to the nonlinear interaction of two waves, the spatial variations of the amplitudes of each frequency component at different locations are shown in Figures8 and 9. The calculated sum of frequency components (2f1 + f2, f1 + 2f2 and 2f1 + 2f2) almost equates 0, and they would not be shown in Figures 8 and 9. In view of the fact that case D1 is a combination of a deep-water wave and a transition water wave, Figure 8(a) shows that the numerical solutions of the two primary frequency components are both reduced and that the amplitude of the primary frequency f1 is smaller than that of f2. It is shown that the relatively short-wave amplitude is reduced more than the relatively long-wave amplitude, which indicates that the relatively short wave is more sensitive by the nonlinear interactions.Figure 8
Comparisons of the spatial variation of amplitudes of each frequency component between calculation and theoretical results for case D1., ,  calculation solutions; , ,  theoretical solutions.
(a)(b)(c)(d)Figure 9
Comparisons of the spatial variation of amplitudes of each frequency component between calculation and theoretical results for case D2., ,  calculation solutions; , ,  theoretical solutions.
(a)(b)(c)(d)Comparing Figures8(b)–8(d), it can be found that the magnitudes of the numerical solution and theoretical solution of each harmonic component are both small, which indicates that the nonlinear effect between the two wave trains of case D1 is weak. This is consistent with the results shown in Figure 6. It can be seen from Figure 8(c) that the sum frequency component (f1 + f2) and the second harmonic components (2f2 and 2f1) can be distinguished and the slight fluctuations appear near the wave source region; however, the fluctuations begin to disappear as the distance of wave propagation increases.Case D2 considers the combination of a deep-water wave and a transition water wave. It is shown in Figure9(a) that the amplitude of the relatively short wave is reduced while the amplitude of the relatively long wave is consistent with the theoretical solution, which indicates that the nonlinearity weakens the energy of the relatively short wave. It can be seen from Figures 9(b) and 9(c) that the difference frequency component (f1 − f2), the sum frequency component (f1 + f2), and the second harmonic components (2f2 and 2f1) present obvious magnitudes. The difference frequency component (f1 − f2) seems to propagate with a certain period. The second harmonic component 2f2 is spatially modulated, and it changes obviously periodically. Under the disturbance of short wave, the amplitude of the second harmonic component of the long-wave changes continuously, showing that the relatively long wave, wave 2, is more affected by wave 1. In view of the difference frequency components between the primary frequencies and the second harmonics, the second harmonics of each primary frequency (2f1 − f2, f1 − 2f2 and 2f1 − 2f2) are distinguishable during the propagation process.In general, the numerical results of the sum frequency components and difference frequency components between the higher harmonics and the higher harmonic and the primary frequency are greater than those of the theoretical solutions. This is because the theoretical solution is derived from the weakly nonlinear Stokes wave theory, while the numerical model established in this study is fully nonlinear. The fact that the obtained numerical solutions can better describe the nonlinear effect between wave components than the theoretical solutions further illustrates that the numerical model of wave-maker by setting two wave source regions can effectively simulate the nonlinear effect between two incident waves.


## 4. Conclusions

The Navier–Stokes equations and the continuity equation were employed as the governing equations, and the mass source wave-maker method was used to establish a numerical wave model based on the Fluent software in this study. Two wave sources were set up inside the numerical wave flume to simulate the propagation of bichromatic waves, and the mass source wave-maker method was thus generalized. The comparisons between the calculation solutions and the theoretical solutions show that the generalized mass source wave-maker method can effectively describe the propagation of the bichromatic waves.The results of frequency spectrum analysis show that the nonlinear interactions between two incident wave trains can result in the generation of higher-order harmonics, and then the nonlinear effect will occur between higher-order harmonics and primary wave components and between higher harmonics themselves. The strength of the nonlinear effect is related to the incident wave parameters. The comparisons between the theoretical solutions and the numerical solutions of the frequency spectrum analysis show that the numerical solution can better describe the nonlinear interactions between wave components than the theoretical solutions of Stokes theory. This is because the numerical model established in this study is fully nonlinear while the theoretical solutions of Stokes theory are weakly nonlinear. It is indicated that the numerical model can realistically perform the nonlinear interaction between two wave trains in real situations.In view of the fact that the actual wave fields are complicated, the nonlinear interactions between different types of waves will be encountered. It can be deduced from the above calculation that the generalized source wave-maker method can simulate the wave field more realistically. We prepare to apply the present generalized mass source method to study the relevant nonlinear effects between different types of incident waves or wave and current in the near future. For example, with the type of the input mass source function changed, the numerical model established in this study may also be applied to the numerical simulation of different types of wave interactions, such as the interactions between Stokes waves and cnoidal waves. On the other hand, with the directions of two incident waves being changed, the model may also be applied to study the interactions between two wave trains with different directions. As one of the wave sources is set as the current source, it may be used to study the wave-current interaction problems.



---

*Source: 8514751-2021-10-31.xml*