A method for generating waves with accurate dispersion

Abstract

This invention discloses a method for generating waves with accurate dispersion, comprising: determining model calculation parameters based on the basic parameters of the ocean region to be simulated; constructing a numerical model of ocean waves with accurate dispersion based on the model calculation parameters; and generating ocean waves with accurate dispersion using the numerical model of ocean waves. The numerical model of ocean waves includes: obtaining an integral-differential governing equation set based on the Euler equations and free surface boundary conditions through continuous Taylor expansion and Fourier transform; the governing equation set includes a continuity equation and a momentum equation; embedding physical process correction terms into the momentum equation based on the governing equation set with embedded physical process correction terms; and performing numerical discretization using the finite difference method based on the governing equation set after embedding the physical process correction terms, wherein time discretization uses a prediction-correction scheme and spatial discretization uses a higher-order difference scheme to achieve the generation of waves with accurate dispersion.

Description

Technical Field

[0001] This invention belongs to the field of fluid dynamics technology in marine engineering, and particularly relates to a method for generating waves with precise dispersion. Background Technology

[0002] Waves are a widespread natural physical phenomenon in the marine environment, formed by multiple factors including wind drive, tidal fluctuations, and air pressure changes. The term "sea wave" usually refers specifically to wind-induced waves, whose development scale is jointly controlled by the duration, speed, and distance of wind action. Wind waves belong to the category of irregular waves, characterized by chaotic waveforms and random propagation directions. When these waves continue to spread outward beyond the wind's influence area, they transform into swells. Significantly different from wind waves, swells are free-flowing waves with regular, smooth crests and troughs, and their waveform regularity increases with distance from the source wind region.

[0003] Wave motion has a crucial impact on sediment transport processes and shoreline geomorphological evolution in nearshore areas. As waves propagate from the open sea to the nearshore, they undergo a series of dynamic responses, including morphological deformation, energy reflection, propagation direction refraction, and wave breaking. After wave breaking, nearshore circulation systems such as coastal currents and rift currents are generated. These systems, combined with wave motion and tidal currents, constitute a complex wave-current coupled dynamic environment. This coupled environment not only profoundly restricts the diffusion paths of nearshore pollutants but also has significant impacts on practical applications such as port construction and nearshore aquaculture.

[0004] In the field of marine engineering, accurately characterizing the propagation and evolution of waves is a core research topic. Currently, the academic community has developed various mathematical models for wave simulation, such as the Boussinesq model and the gentle slope equation model. However, these models still have significant shortcomings in terms of dispersion characteristic representation, nonlinear effect description, and the applicability of wave spectra. For example, the Boussinesq model uses an approximate calculation method to describe dispersion characteristics, resulting in significant limitations in high-frequency wave simulation scenarios; the gentle slope equation model is only suitable for simulating regular waves or narrow-spectrum irregular waves. The simulation requirements for high-frequency waves and wide-spectrum irregular waves in actual marine environments place higher demands on models, including high-precision nonlinear characterization, accurate reproduction of dispersion characteristics, and wide-spectrum applicability. Therefore, developing a new wave mathematical model with high-precision simulation capabilities has become an urgent technical challenge in the field of marine engineering. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention proposes a method for generating waves with accurate dispersion, which overcomes the issues of low computational accuracy and limited applicability of existing ocean wave numerical models.

[0006] To achieve the above objectives, the present invention provides a method for generating waves with precise dispersion, comprising: Based on the basic parameters of the ocean area to be simulated, determine the model calculation parameters; Based on the calculated parameters of the model, a numerical model of ocean waves with accurate dispersion is constructed, and ocean waves with accurate dispersion are generated using the numerical model of ocean waves. The construction of the ocean wave numerical model includes: based on the Euler equations and free surface boundary conditions, obtaining an integral-differential control equation set through continuous Taylor expansion and Fourier transform. The control equation set includes a continuity equation and a momentum equation, wherein the continuity equation introduces an integral-differential operator E to achieve a broadband wave description. According to the governing equations, physical process correction terms are embedded in the momentum equation, including underwater friction terms, turbulent mixing terms, and wave breaking energy dissipation terms. Based on the control equations with embedded physical process correction terms, the finite difference method is used for numerical discretization. The time discretization adopts the prediction-correction scheme, and the spatial discretization adopts the high-order difference scheme. The integral operation is accelerated by the CUDA-GPU heterogeneous parallel computing architecture to obtain the spatiotemporal distribution of free surface rise and realize the generation of waves with accurate dispersion.

[0007] Optionally, the basic parameters include water depth conditions, spectral peak frequency, spectral peak period, statistical characteristic wave height of incident waves, and incident azimuth angle of long-period waves. The model calculation parameters include computational domain range parameters, spatiotemporal resolution parameters, bottom friction coefficient, lateral mixing coefficient, and wave breaking parameters.

[0008] Optionally, the ocean wave numerical model with accurate dispersion is: ; ; in, For the free surface to rise, For time, For the horizontal gradient operator, Because the local water is still and deep, It is the acceleration due to gravity. Let be the horizontal velocity vector at the free surface. For the integral-differential operator.

[0009] Optionally, the integral-differential operator for: ; in, k For wave number, h Because the local water is still and deep, For the horizontal gradient operator, x For the coordinates of the control points, Let i be the coordinates of the integration point, and i be the imaginary unit.

[0010] Optionally, embedding a seabed friction term into the momentum equation according to the governing equations includes: , ; in, The coefficient of friction at the bottom of the water is a constant. for x underwater friction in the direction of the water. for y underwater friction in the direction of the water. For the total water depth, for x Average speed in direction over time for y Average speed over time in the direction.

[0011] Optionally, embedding turbulent mixing terms into the momentum equation according to the governing equations includes: ; ; in, Where is the eddy viscosity coefficient. for x Turbulent mixing forces in the direction of flow for y Turbulent mixing forces in the direction of flow.

[0012] Optionally, embedding wave breaking energy dissipation terms into the momentum equation according to the governing equations includes: ; ; in, for x The direction of wave breaking force, for y The direction of wave breaking force, The value is the eddy viscosity coefficient caused by wave breaking.

[0013] Optionally, the numerical discretization solution using the finite difference method is performed based on the set of governing equations after embedding the physical process correction terms, including: The prediction step is calculated based on the fifth-order Adams-Bashforth scheme to obtain the prediction value for the next time level; The correction step is calculated using the sixth-order Adams-Moulton scheme to obtain the correction value for the next time layer. The convergence of the correction step is determined based on the preset error limit. When the number of iterations exceeds the threshold, the weighted average method is used to correct the two correction values ​​before and after to obtain the optimized correction value of the current time step. The spatial derivative term is discretized using the five-point central difference scheme, and the infinite integration interval is transformed into a finite interval for numerical integration based on the rapid decay characteristic of the kernel function.

[0014] Compared with the prior art, the present invention has the following advantages and technical effects: 1. The new ocean wave model proposed in this invention can accurately describe the strong nonlinear characteristics of wave propagation and breaking on submerged breakwaters and slopes, and compared with other models, it can give the spatial distribution characteristics of each harmonic more closely with experimental results.

[0015] 2. The momentum equation of this invention is based on the Euler equation. The derivation process is completely accurate without any approximations, which is very different from other models. This also lays the foundation for the high-precision numerical algorithm of this model.

[0016] 3. The continuity equation of this invention is based on free surface boundary conditions, which can account for vortex-induced rotational flow and has completely accurate dispersiveness and shallowing properties. The equation does not contain any wavenumber-related coefficients, and there are no restrictions on water depth during the derivation process; therefore, irregular waves at arbitrary water depths can be studied.

[0017] 4. This invention uses a high-order numerical scheme to solve the equation numerically, so that the accurate dispersion of the equation can be realized in numerical calculation. Attached Figure Description

[0018] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments and descriptions of this application are used to explain this application and do not constitute an undue limitation of this application. In the drawings: Figure 1 This is a diagram showing the configuration of the computational domain and virtual mesh according to an embodiment of the present invention; Figure 2 This is a distribution diagram of the kernel functions in an embodiment of the present invention; Figure 3 This is a parallel flowchart of an embodiment of the present invention; Figure 4 This is a comparison diagram of the wavefront elevation of the computational model in this embodiment of the invention and the theoretical solution; Figure 5 This is a comparison graph of wavefront elevation and experimental results for different models in embodiments of the present invention; Figure 6 This is a topographic map of a double-groove sandbar according to an embodiment of the present invention; Figure 7 These are comparison diagrams of the time-averaged flow field in embodiments of the present invention, where (e) and (f) are numerical results obtained by applying the model under different side mixing coefficients, respectively. Figure 8 These are wave height distribution diagrams for different models in embodiments of the present invention, wherein (a) is a comparison result diagram for x=16m, (b) is a comparison result diagram for x=15.2m, and (c) is a comparison result diagram for x=14m; Figure 9 These are distribution diagrams of the time-averaged velocity of different models in the embodiments of the present invention. (a) is a comparison result diagram of x=16m, (b) is a comparison result diagram of x=15.2m, and (c) is a comparison result diagram of x=14m. Figure 10 This is a distribution diagram of wave height and water level changes in the direction perpendicular to the shoreline according to an embodiment of the present invention. Among them, (a) is a comparison diagram of wave height in the direction of y=9.2m, (b) is a comparison diagram of wave height in the direction of y=4.6m, (c) is a comparison diagram of water level changes in the direction of y=9.2m in the direction perpendicular to the shoreline, and (d) is a comparison diagram of water level changes in the direction of y=4.6m in the direction perpendicular to the shoreline. Figure 11 This is a flowchart of a method for generating waves with precise dispersion according to an embodiment of the present invention. Detailed Implementation

[0019] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0020] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0021] This embodiment proposes a method for generating waves with accurate dispersion, such as... Figure 11 As shown, the specific steps include: Based on the basic parameters of the ocean area to be simulated, determine the model calculation parameters; Based on the calculated parameters of the model, a numerical model of ocean waves with accurate dispersion is constructed, and ocean waves with accurate dispersion are generated using the numerical model of ocean waves. The construction of the ocean wave numerical model includes: based on the Euler equations and free surface boundary conditions, obtaining an integral-differential control equation set through continuous Taylor expansion and Fourier transform. The control equation set includes a continuity equation and a momentum equation, wherein the continuity equation introduces an integral-differential operator E to achieve a broadband wave description. According to the governing equations, physical process correction terms are embedded in the momentum equation, including underwater friction terms, turbulent mixing terms, and wave breaking energy dissipation terms. Based on the control equations with embedded physical process correction terms, the finite difference method is used for numerical discretization. The time discretization adopts the prediction-correction scheme, and the spatial discretization adopts the high-order difference scheme. The integral operation is accelerated by the CUDA-GPU heterogeneous parallel computing architecture to obtain the spatiotemporal distribution of free surface rise and realize the generation of waves with accurate dispersion.

[0022] Specifically, to fundamentally address the problem of low dispersion in numerical calculations, this invention first establishes a new ocean wave equation with accurate dispersion based on the Euler equation and free surface boundary conditions through continuous Taylor expansion. Then, using numerical calculation methods, this equation is modeled, incorporating physical phenomena such as seabed friction, wave breaking, turbulent mixing, wave generation, and wave ramping. This model can accurately describe the strongly nonlinear characteristics of wave propagation and breaking on submerged breakwaters and slopes, and can provide spatial distribution characteristics of harmonics on slopes that closely approximate experimental results. The operational steps are given below: Step 1, Calculation model settings: Define the basic information of the ocean area to be calculated, specifically covering core elements such as water depth conditions, spectral peak frequency, spectral peak period, statistical characteristic wave height of incident waves, and incident azimuth angle of long-period waves. Based on this, specifically adjust the model calculation parameters: The high-precision fully dispersive ocean wave numerical model used in this embodiment is an integral-differential wave control equation system derived from Euler equations and free surface kinematics and dynamic boundary conditions through continuous Taylor expansion and Fourier transform / inverse transform. It has the characteristics of no approximate accurate dispersion and strong nonlinear characterization capability.

[0023] The core features of the model: Theoretical basis: The Euler equations for a uniform, ideal, incompressible fluid are used as the starting equations. The derivation process does not make any approximations, and the momentum equation is completely accurate. The continuity equation is constructed based on the boundary conditions of a free surface. Broad-spectrum waves (including high-frequency waves and irregular waves) are described by kernel function integration. Rotational flow caused by vertical vortices can be considered.

[0024] Physical module integration: By embedding mathematical expressions of key physical processes such as seabed friction, turbulent mixing, wave breaking, wave generation, and wave ramping into the basic control equations, it can simulate the dynamic response of waves from the open sea to the nearshore (deformation, refraction, breaking, nearshore circulation generation, etc.).

[0025] Numerical solution method: Spatial discretization is performed using the finite difference method, and temporal discretization uses the prediction-correction scheme; To meet the need for efficient integral solving, a CUDA-GPU heterogeneous parallel computing architecture is introduced to accelerate integral operations and ensure the efficiency of large computational domain and high-precision simulation.

[0026] Scope of application: No water depth limitation, applicable to any working condition such as deep water, shallow water, and transitional water depth; supports simulation of various wave types such as regular waves, broadband irregular waves, and wind-wave-surge coupling, especially suitable for wave propagation and breaking simulation under complex terrains such as submerged dikes, slopes, and sandbars-ditches.

[0027] Specific debugging process for model calculation parameters: The purpose of adjusting the model calculation parameters is to ensure the accuracy and stability of the simulation results. It needs to be carried out in accordance with the principle of "initial value setting → step-by-step debugging → verification and optimization → comprehensive confirmation" in combination with the specific research conditions (wave type, terrain conditions, simulation target). The specific process is as follows: 1: Pre-debugging phase—basic parameter review and initial value setting: Define the core operating parameters: Based on measured data or experimental design of the ocean area to be simulated, determine the baseline values ​​of the basic input parameters, including: Topographic parameters: spatial distribution of still water depth h, topographic slope (e.g., submerged breakwater slope 1:20, sandbar slope 1:30); Wave parameters: peak period Tp, peak frequency fp, characteristic wave height Hs, incident azimuth; Spatiotemporal resolution parameters: The spatial step size is set according to the wavelength ratio, with an initial value of λ / 40 (λ is the incident wave wavelength); for complex terrain, λ / 5 can be used; the time step size satisfies the CFL stability condition, with an initial value of T / 80 (T is the incident wave period). Bottom friction parameters: 0.01~0.03 for regular wave conditions; 0.03~0.08 for irregular / fragmented wave conditions; can be appropriately increased for rough terrain. Turbulent mixing parameters: Based on the experience of the implementation examples, the initial value is 1.0; when simulating strong turbulent flow fields such as cracked flow, it can be 0.25~0.5, and when simulating weak turbulence, it can be 2.0~3.0; 2: Step-by-step debugging – targeted optimization of core parameters: (1) Debugging the computational domain range: Incident boundary: A wave energy development zone of at least 5 times the wavelength must be reserved to ensure that the incident wave reaches a stable state before entering the target simulation area; Outgoing boundary: A wave-absorbing zone of at least 3 times the wavelength must be reserved, or a sponge layer wave-absorbing technology should be used to prevent waves from reflecting back into the simulated area; Lateral boundary: Adjusted according to the azimuth of wave incidence. If it is vertical incidence, the distance between the lateral boundary and the target area should not be less than twice the wavelength. If it is oblique incidence, the lateral range should be expanded according to the refraction angle.

[0028] Verification criteria: Monitor the stability of wave parameters (wave height, period) within the target area. If the parameter fluctuation is less than 5%, the calculation domain range is reasonable.

[0029] (2) Spatiotemporal resolution adjustment: Spatial step convergence test: With the time step and other parameters fixed, simulations are performed with different spatial steps in sequence, and the deviation between the wavefront rise time series of the target region and the theoretical solution (such as the fifth-order Stokes wave solution) is compared. Time step convergence test: Fix the optimal spatial step size and monitor computational stability (no numerical oscillation) and result accuracy; Verification criteria: When the spatial step size is reduced to λ / 40, the deviation of wavefront rise is less than 3%, and the deviation does not decrease significantly when the step size is further reduced. This step size is the optimal spatial step size. The time step size must meet the CFL condition, and the simulation results must be free of oscillations and have smooth and continuous waveforms.

[0030] (3) Adjustment of bottom friction coefficient: Select a flat-bottomed terrain without wave breakage, set different bottom friction coefficients for simulation, and compare the simulated wave height attenuation rate along the path with the measured / experimental attenuation rate. If the simulated attenuation rate is less than the measured value, increase fw; if it is greater than the measured value, decrease fw.

[0031] (4) Adjustment of lateral mixing coefficient: For nearshore circulation simulation conditions such as rift flow and coastal flow, multiple sets of cm (e.g., 0.25, 0.5, 1.0, 2.0) were set; the consistency between the simulated time-averaged flow field (velocity magnitude and direction) and experimental data (e.g., Haller double-groove sandbar experiment) was compared. (5) Wave breaking parameter adjustment: Adjusting the crushing judgment parameter B: Increasing B will delay the occurrence of crushing, while decreasing B will cause crushing to occur earlier. The optimal value is determined by comparing the simulated crushing position with the experimental observation position. Adjusting the crushing strength parameter α: The larger α is, the stronger the crushing energy dissipation. It needs to be matched with the attenuation amplitude of the crushing wave height in the experiment.

[0032] 3: Comprehensive verification and iterative optimization: Multi-parameter combination verification: Combine the optimal parameters obtained from each step of debugging, conduct full-condition simulation, and compare the overall consistency of the core output indicators (wave rise, wave height distribution, and time-averaged flow field) with the experimental / measured data.

[0033] Sensitivity analysis: For key parameters (such as lateral mixing coefficient and crushing strength parameter), a disturbance of ±20% is applied to analyze the variation range of the simulation results and clarify the sensitivity level of the parameter; for highly sensitive parameters, further fine-tuning is required.

[0034] Iterative optimization: If the deviation between the simulation results and the experimental data exceeds the allowable range, return to the step-by-step debugging stage of the corresponding parameters, adjust the initial values ​​and recalculate until the accuracy requirements are met.

[0035] 4: Parameter Fixation and Output: The finalized optimal parameter set will be compiled and archived, including the computational domain range, spatiotemporal resolution, bottom friction coefficient, lateral mixing coefficient, wave breaking parameters, etc., as the standard parameter configuration for this type of working condition; at the same time, key data during the parameter debugging process (such as convergence test curves and deviation comparison tables) will be recorded to provide a reference for the simulation of similar working conditions.

[0036] The computational domain is defined based on actual working conditions, and the spatial and temporal resolutions are optimized. At the same time, key indicators such as bottom friction coefficient, lateral mixing coefficient, and wave breaking parameters are flexibly adjusted according to simulation requirements.

[0037] Step 2: Perform wave numerical simulation using the new computational model: Step 1, setting up the computational model, and step 2, using the new computational model to perform wave numerical simulation, are an organic whole that is interconnected, mutually supportive, and has a unified goal. Together, they constitute a complete technical chain of "preparation-execution".

[0038] One of the core tasks of Step 1 is to identify the core basic parameters of the ocean area to be simulated. These parameters are essential input conditions for numerical simulation in Step 2. In Step 2, the basic parameters determined in Step 1 need to be substituted into the governing equations (1)-(3) to assign specific operating conditions to the equations. Without the parameter input from Step 1, equations (1)-(3) are merely general mathematical expressions and cannot be used for accurate simulation of specific sea areas and specific wave conditions.

[0039] The second core task of Step 1 is to perform targeted debugging of the model's computational parameters. The debugging results directly determine the accuracy, stability, and reliability of the numerical simulation in Step 2. A reasonable setting of the computational domain can avoid boundary reflection interference during the simulation in Step 2, ensuring that the wave propagation process conforms to actual physical laws. Optimized configuration of the spatiotemporal resolution can eliminate numerical dispersion errors in solving the equations in Step 2, ensuring the good agreement between the simulated waveform and the theoretical solution / experimental results. Debugging parameters such as bottom friction, turbulent mixing, and wave breaking can ensure that Step 2 closely matches actual physical phenomena when simulating bottom energy dissipation, nearshore turbulence effects, and wave breaking processes.

[0040] The logical sequence of the model definition: Step 1 clarifies the core characteristics (theoretical basis, physical module integration, numerical solution method, and applicable scope) of the high-precision fully dispersive ocean wave numerical model used in this method. This is the logical basis for selecting and using the model in Step 2 to conduct simulations. The governing equations (1)-(3) used in Step 2 are the mathematical carriers of the model defined in Step 1; the physical phenomena mentioned in Step 2, such as "bottom friction, wave breaking, and turbulent mixing," also completely correspond to the content of "model integration of key physical modules" in Step 1. The two are a direct correspondence between "theoretical definition of the model" and "mathematical implementation of the model."

[0041] Substitute the parameters determined in step 1 into the following calculation model: The parameters specifically include: Topographic parameters: spatial distribution of local still water depth, topographic slope (e.g., submerged breakwater slope 1:20, sandbar slope 1:30); Wave parameters: peak frequency, peak period, statistical characteristic wave height of incident waves, and incident azimuth of long-period waves.

[0042] Computational domain configuration parameters: range of incident boundary wave energy development zone, range of outgoing boundary wave absorption zone, and range of lateral boundary extension; Spatiotemporal resolution parameters: spatial step size optimized after convergence testing, and time step size that satisfies the CFL stability condition; Physical process parameters: bottom friction coefficient, lateral mixing coefficient, wave breakage judgment parameter B, breakage strength parameter α, etc.

[0043] (1); (2); In the formula, For the free surface to rise, For time, For the horizontal gradient operator, Because the local water is still and deep, It is the acceleration due to gravity. Let be the horizontal velocity vector at the free surface. Integral-differential operators are given below. The expression: (3); In the formula, k Let be the wave number. When applying the above equations (1)-(2) to the simulation of actual ocean waves, it is also necessary to consider physical phenomena including underwater friction, wave breaking, turbulent mixing, wave generation, and wave rampage.

[0044] More specifically, equations (1) and (2) are the core coupled control equations of the high-precision fully dispersive ocean wave numerical model constructed in this invention. They are derived based on the Euler equation and the free surface boundary conditions, and are mutually premised and mutually constrained, jointly describing the mass conservation and momentum conservation laws of wave motion.

[0045] Equation (1) is the continuity equation: it is derived based on the kinematic boundary conditions of the free surface, combined with Taylor expansion and Fourier transform / inverse transform, and introduces the integral-differential operator E. Its core function is to describe the mass transfer relationship between the free surface rise and the horizontal velocity vector at the free surface, reflecting the mass conservation characteristics of wave motion. This equation describes a wide spectrum of waves through kernel function integration, and has the advantages of no water depth limitation and the ability to consider vortex flow.

[0046] Equation (2) is the momentum equation: it is directly derived from the Euler equation. The derivation process does not make any approximations. Its core function is to describe the spatiotemporal variation of the horizontal velocity vector at the free surface and to reflect the momentum conservation characteristics of wave motion. This equation can be embedded in the mathematical expression of physical processes such as underwater friction, turbulent mixing, and wave breaking, and accurately characterizes the dynamic response in wave propagation.

[0047] The coupling relationship between equations (1) and (2) in the solution process needs to be solved in a coupled manner. From equation (1), the expression for the horizontal velocity vector at the free surface can be derived from the free surface rise through the operation of the integral-differential operator E. Substituting this velocity expression into equation (2), the spatiotemporal distribution of the free surface rise can be further solved. The two derivations form a closed solution system. Combined with the parameters determined in step 1 and the physical process correction terms, a high-precision simulation of the entire process of wave propagation, deformation, and breaking is finally achieved. The integral-differential form of equation (1) ensures the complete and accurate dispersion characteristics of the model, which can accurately reproduce the propagation law of broadband waves. The approximation-free derivation of equation (2) ensures the strong nonlinear characterization ability of the model, which can accurately characterize the breaking of waves and the nearshore circulation generation process under complex terrain.

[0048] Furthermore, the basic parameters include water depth conditions, spectral peak frequency, spectral peak period, statistical characteristic wave height of incident waves, and incident azimuth angle of long-period waves; The model calculation parameters include computational domain range parameters, spatiotemporal resolution parameters, bottom friction coefficient, lateral mixing coefficient, and wave breaking parameters.

[0049] Furthermore, the numerical model of ocean waves with accurate dispersion is as follows: (4); (5); in, For the free surface to rise, For time, For the horizontal gradient operator, Because the local water is still and deep, It is the acceleration due to gravity. Let be the horizontal velocity vector at the free surface. For the integral-differential operator.

[0050] Furthermore, the integral-differential operator E is: (6); in, k For wave number, h Because the local water is still and deep, For the horizontal gradient operator, x For the coordinates of the control points, Let i be the coordinates of the integration point, and i be the imaginary unit.

[0051] Furthermore, when using the new ocean wave model for simulation calculations, the physical phenomenon of seabed friction is considered. To address this phenomenon, a seabed friction term is added to the right-hand side of the model's momentum equation. The components of the expression are as follows: , (7); in, The coefficient of friction at the bottom of the water is a constant. for x underwater friction in the direction of the water. for y underwater friction in the direction of the water. For the total water depth, for x Average speed in direction over time for y Average speed over time in the direction.

[0052] Furthermore, when using the new ocean wave model for simulation calculations, the turbulent mixing and dissipation phenomenon of waves must be considered. To address this phenomenon, a turbulent mixing term is added to the right-hand side of the model's momentum equation. The specific solution process is as follows: This embodiment uses a Smagorinsky-type subgrid model to consider the influence of eddy viscosity on turbulent mixing in the flow field, as expressed below: (8); (9); in, Where is the eddy viscosity coefficient. for x Turbulent mixing forces in the direction of flow for y Turbulent mixing forces in the direction of flow.

[0053] (10); In the formula, and They are respectively x and y The average velocity at the still water level in the direction of the current. This is the lateral mixing coefficient, which is generally taken as 0.05-3.0.

[0054] Furthermore, when using the new ocean wave model for simulation calculations, the physical phenomenon of wave breaking must be considered. To address this phenomenon, a wave breaking energy dissipation term is added to the right-hand side of the model's momentum equation. By introducing eddy viscosity terms To address the energy dissipation caused by wave breaking, an extension term is added to the right-hand side of the momentum equation. Its expression is: (11); (12); in, for x The direction of wave breaking force, for y The direction of wave breaking force, The eddy viscosity coefficient is induced by wave breaking. ; in Wave breaking intensity, To determine the parameters of wave breaking, when At that time, the waves did not break. The possible values ​​are as follows: (13); in, Used to define the start and end times of wave breaking: (14); In the formula, The start time of the breakage. For the duration of wave breaking, This is an adjustable parameter, where t is the current time. and These are the crushing start parameter and the crushing termination parameter, respectively, and their expressions are: (15); in, and It is a constant.

[0055] Furthermore, based on the governing equations after embedding the physical process correction terms, the numerical discretization solution is performed using the finite difference method, including: The prediction step is calculated based on the fifth-order Adams-Bashforth scheme to obtain the prediction value for the next time level; The correction step is calculated using the sixth-order Adams-Moulton scheme to obtain the correction value for the next time layer. The convergence of the correction step is determined based on the preset error limit. When the number of iterations exceeds the threshold, the weighted average method is used to correct the two correction values ​​before and after to obtain the optimized correction value of the current time step. The spatial derivative term is discretized using the five-point central difference scheme, and the infinite integration interval is transformed into a finite interval for numerical integration based on the rapid decay characteristic of the kernel function.

[0056] Example 1: This example provides the derivation process of a new ocean wave equation. First, based on the Euler equation and free surface boundary conditions, the equation is established using velocity and wave rise at the free surface as variables. This equation is a continuous broadband wave equation; the continuous spectrum is achieved through integration of the kernel function and possesses accurate dispersion. Compared to narrow-spectrum wave equations, this equation has no frequency limitation and can describe free waves and broadband irregular waves of arbitrary frequencies. This equation overcomes the limitation of narrow-spectrum equations that cannot consider swirling wave motions, and can be used to study wave phenomena such as coastal currents and nearshore circulation. Specifically, it includes the following steps: Step 1, Derivation of the momentum equation: The origin of the coordinate system is located at the still water surface, and the z-axis is set vertically upward. The fluid under study is a homogeneous ideal fluid with incompressible and inviscid properties, and the influence of surface tension is neglected. To ensure that the derived equations are applicable to swirling flow scenarios, the derivation process will be based on the Euler equations and the kinematic boundary conditions of free surfaces.

[0057] (16); (17); (18); Where w is the vertical velocity. Let be the density and p be the pressure. The equation is expressed at a free surface using the relationship: (19); In the formula, The variable is located at the free surface. Equations (16) and (17) are expressed at the free surface, and the pressure term is eliminated. The result is then simplified to: (20); Equation (20) is completely accurate because it makes no assumptions in the derivation process.

[0058] Step 2, give the governing equations: The continuity equation is established based on the kinematic boundary conditions (18), which requires... Using the vertical velocity at the still water surface and horizontal speed To express, and then use To express. and Perform a Taylor expansion at a still water surface, that is: , (twenty one); In the formula Let be the horizontal coordinate vector. Using the continuity equation and the irrotational equation: , (twenty two); This can be expressed as: (twenty three); (twenty four); The following is based on the underwater boundary conditions. and Relationship: (25); Substituting equations (23) and (24) into the equations, we get: (26); In the formula, (27); After sorting, we get: (28); Using the principle of Fourier transform, let: (29); Substituting equation (29) into equation (28), we get: (30); The following inverse Fourier transform yields the following: (31); In the formula, . and These are the zeroth-order and first-order kernel functions, respectively, and their expressions are: , (32); make: (33); In the formula, express The water depth gradient term. Combining all terms in the above equation, we can simplify it to: (34); In the formula, the expression for the integral-differential operator E is: (35); The above relationship can be used to establish and and The relationship can be obtained through continuous iteration: (36); Substituting equation (36) into equation (18) yields the expression for the continuity equation. Therefore, the wave equation is: (37); (38); The above (37) and (38) are the newly established ocean wave equations.

[0059] Example 2 – Solving the Equation: When applying the ocean wave equations derived in Example 1 to actual ocean wave simulation scenarios, the influence of physical phenomena such as seabed friction, wave breaking, turbulent mixing, wave generation, and wave ramping needs to be incorporated into the equation system. The following is a brief explanation of how to handle these physical phenomena, followed by a detailed description of the specific numerical solution process.

[0060] (1) Underwater friction treatment: With the addition of the underwater friction term to the equation, the components of the expression are as follows: , (39); In the formula, For the total water depth, The coefficient of friction is the bottom friction factor, and its specific value depends on conditions such as wave type.

[0061] (2) Turbulent mixing: The turbulent mixing dissipation term is the eddy viscosity dissipation induced by wave breaking. This embodiment uses a Smagorinsky-type subgrid model to consider the influence of eddy viscosity on turbulent mixing in the flow field, expressed as: (40); (41); in, The eddy viscosity coefficient is given by the following formula: (42); In the formula, and Let be the average velocities at still water level in the x and y directions, respectively. This is the lateral mixing coefficient, which is generally taken as 0.05-3.0.

[0062] (3) Wave breaking term: The expression for energy dissipation caused by wave breaking is: (43); (44); In the formula, The eddy viscosity coefficient caused by wave breaking is expressed as follows: (45); in Let B be the wave breaking intensity, and let B be the parameter for judging wave breaking. When B=0, the wave does not break. The values ​​of B are as follows: (46); in, Used to define the start and end times of wave breaking: (47); In the formula, The start time of the breakage. For the duration of wave breaking, This is an adjustable parameter, where t is the current time. and These are the crushing start parameter and the crushing termination parameter, respectively, and their expressions are: (48); (4) Numerical solution of the model: This embodiment uses the finite difference method as its core numerical solution system, specifically for the numerical discretization and solution of wave equations. From the perspective of the characteristics of the governing equations, the continuity equation has explicit solution conditions, while the momentum equation has a mixed derivative term of velocity with respect to time on the right-hand side. The presence of this term prevents the momentum equation from being directly solved explicitly. Therefore, it is necessary to specifically extract and separate the time-mixed derivative component in the momentum equation to decouple the equations before proceeding with subsequent calculations. The numerical scheme design considers two core requirements: first, ensuring high accuracy in equation discretization; and second, guaranteeing stable convergence in the calculation process. The time layer uses the fifth-order forecast and sixth-order correction Adams-Bashforth-Moulton (ABM) scheme, and the spatial layer uses a five-point fourth-order difference scheme. Detailed numerical solution steps are as follows.

[0063] The momentum equation (23) needs to be decoupled before numerical solution. This is achieved by simulating the two component equations. and The momentum equation (23) is decomposed into component equations as the variables to be solved. The momentum equation (23) is then expressed in component form: (49); (50); Moving the velocity-time first derivative terms in equations (49) and (50) to the left side of the equations and the remaining terms to the right side, we get: , (51); In the formula, , , (52); (53); (54); Next, separate the first time derivative term containing velocity in the formula, that is: , (55); Equation (55) is the momentum equation after decoupling, which can be solved explicitly.

[0064] Numerical discretization of the governing equations: After the above decoupling, both the continuity equation and the momentum equation can be solved explicitly. To extend the applicability of the equations, energy dissipation terms caused by bottom friction, wave breaking, and the mixed subgrid effect are added to the momentum equation. Therefore, the governing equations are written in the following form: (56); (57); (58); In the formula, , and These are the bottom friction term, the turbulent mixing term, and the energy dissipation term caused by wave breaking, respectively. The governing equations employ a fifth-order forecast, sixth-order corrected ABM scheme at the time level. According to the numerical scheme requirements, the time level is divided into six layers, each layer being... and .in Each physical quantity on the layer is responsible for calculating the right-hand side of the equation, and the result can be obtained. The predicted value for the layer, i.e.: (59); In the formula, the subscript p represents the prediction step. Represents the left-hand side of equations (56) to (58) , For all terms on the right-hand side of the equation. Using... and The physical quantities on the layer can be calculated. The values ​​of each physical quantity at each time step of the correction are: (60); Calculation results The following error requirements must be met: (61); In the formula, the subscript k represents the number of iterations, and the error limit is... When waves exhibit strong nonlinear characteristics, the number of iterations in the aforementioned numerical iteration process increases significantly, affecting computational efficiency and potentially increasing the risk of convergence fluctuations. Therefore, this paper specifically employs an accelerated iterative convergence method. The specific rule is as follows: when the number of correction step iterations exceeds a preset threshold, the current correction step calculation result is determined to have deviated from the true value. Blind iteration is discontinued, and instead, a weighted average of the calculation values ​​from the two preceding and following correction steps is used to reassign the correction step value. The expression for this is: (62); In the formula, This is the correction value for the current time step. Previous time step correction value, This is the corrected value for the current time step.

[0065] The right-hand side of the momentum equations (57) and (58) contains the first and second time derivative terms, which are calculated using the following formula.

[0066] Prediction Step: (63); (64); Calibration step: (65); (66); In the formula, the subscript represents the grid space marker, and the superscript represents the time layer.

[0067] In this embodiment, the spatial discretization process of the equations is implemented based on a non-interlaced rectangular grid system. Considering the technical requirements of applying the central difference scheme to all nodes in the computational domain, a method of adding virtual grids is adopted. The control equations constructed in this embodiment are integral-differential equations. This characteristic determines that the virtual grid outside the boundary must not only meet the accuracy requirements of differential discretization, but also conform to the numerical conditions of integral operations. Based on this, kb layers of virtual grids are extended outside the boundary of the computational domain. Comparing the grid requirements for differential solutions and integral operations, it can be seen that the number of grids required for integral operations is greater. Therefore, the number of kb layers needs to be determined by the number of virtual grids corresponding to integral operations. The specific mathematical expression of kb will be derived in detail in the following content. Among them, the values ​​of the nodes inside the computational domain can be directly obtained through iterative solution of the equations, while the values ​​of the virtual grid nodes are assigned according to the preset boundary conditions. The spatial layout of the computational domain and the virtual grid is shown in Figure 1. For the internal nodes in the computational domain, a five-point central difference scheme is used for spatial discretization. For example, its corresponding difference scheme is: (67); (68); (69); The equation also includes mixed derivatives, and its difference scheme is as follows: (70); (71); Calculation of the integral kernel function: Continuity equations are differential-integral equations. Solving for the integral term is an important step in solving the equation, and it can be written in the following form: (72); Analysis of the above equation reveals that the integral operator T belongs to the improper integral, which lacks the conditions for direct solution. Therefore, it is necessary to first solve the numerical computation problem of improper integrals. The core idea is to systematically analyze the kernel function in the equation. Figure 2, combining the specific expressions of the zeroth-order and first-order kernel functions, details the spatial distribution characteristics of the kernel function. From the distribution results, it can be observed that although... The integral is a improper integral, but the kernel function has... The kernel function exhibits a rapid decay over time, with its value approaching zero after decay. Utilizing this key property, the original infinite integration interval can be transformed into a finite interval for numerical integration, significantly reducing computational complexity. Based on these properties of the kernel function, control points... The integral of position can be transformed into the following form: (73); (74); In the formula, The parameters are defined as follows: Control points The actual water depth at the location is the integration limit parameter, and in this embodiment, the value of this parameter is set to [value missing]. It is important to note that the second-order kernel function differs significantly from the zeroth-order and first-order kernel functions; its... Zero point position ( A singularity exists, which manifests as a point at zero. The kernel function value approaches infinity. Considering the differences and common characteristics of various kernel functions, this embodiment uniformly uses a logarithmic (ln function) polynomial for fitting the zeroth, first, and second-order kernel functions. The fitting relationship is as follows: , (75); From this fitting formula, it can be further deduced that the range of the integration interval is determined by the water depth, exhibiting a significant positive correlation: the greater the water depth, the larger the required integration range; the smaller the water depth, the smaller the integration range. In subsequent numerical calculations, the number of virtual mesh layers outside the boundary is determined precisely by... Parameters and Parameter coordinated control ensures integration accuracy and computational stability.

[0068] Parallel algorithms: The wave control equations constructed above are integral-differential equations. For continuous equations, the efficiency of integration directly determines the overall solution efficiency. To clearly explain the integration optimization approach, this section selects a simplified scenario of first-order equations and only considers the flat-bottom term. This numerical calculation uses a horizontal two-dimensional spatial discretization mode. We define the number of grid nodes in the x-direction as M and the number of grid nodes in the y-direction as N. The total number of grid nodes in the entire computational domain is then... The corresponding first-order governing equations for continuity are expressed as follows: (76); Analysis of the above formula shows that each grid coordinate point within the computational domain needs to independently complete a double integration operation, corresponding to a total number of calculations within a single time step. ,in and Represent direction and The number of integral iterations in the direction, this and The values ​​of the two parameters are closely related to the water depth of the computational region, showing a clear positive correlation. That is, the greater the water depth, the wider the integration coverage, and the more integration iterations are required.

[0069] If the integration calculation relies solely on CPU serial computation, it will consume a significant amount of computation time, making it difficult to meet the requirements for efficient solution. Crucially, the continuity equation constructed in this paper approximates the nonlinear terms to fifth-order precision, and four new second-order water depth gradient-related terms have been added to the integration module, resulting in a multiplicative increase in overall computational load. Therefore, a heterogeneous parallel computing model combining GPU and CPU is necessary to handle the integration task. This model uses the CUDA-GPU parallel computing algorithm to achieve efficient integration. The core implementation process of this algorithm includes: allocating dedicated storage space in the GPU device's memory, initializing the computational data and transferring it between the host and device, declaring data attributes and defining the valid data region and lifecycle, and after the integration calculation is completed, sending the result data back to the CPU host and releasing the memory space occupied by the GPU device. From the perspective of GPU's physical hardware architecture design logic, its core positioning is to support large-scale data parallel computing. The larger the data scale, the more significant the acceleration effect of parallel computing. Conversely, if the computing task only needs to call a small number of computing cores, resulting in most cores being idle, then such tasks are more suitable for CPU serial execution, which can effectively improve overall computing efficiency. In addition, auxiliary operations such as parameter initialization, output of calculation results, and persistent storage are also not suitable for GPU execution due to their small computational load and strong interactivity. Based on the above characteristics, the porting of GPU parallel computing needs to follow three core steps: first, systematically analyze the algorithm architecture of the original program and identify the parallelization adaptation nodes; second, implement parallelization reconstruction for core computing modules; and third, optimize the data transmission process between the host and the device and the GPU computing core scheduling strategy.

[0070] The core computational load of this model is concentrated in the integration module, and its parallel computing process is described in detail below: First, the CPU performs full parameter initialization, including reading the initial input file, wavenumber calculation, computational domain terrain generation, and solving for correlation coefficients in the integration operation. Second, the various parameters required for GPU parallel computing are transferred from the CPU host to the GPU device, where the integration term is decomposed into... , , and The four-layer architecture allocates threads based on the number of CUDA cores in the GPU, with layers 1 and 2 assigned to the gank computing module. and The computation is distributed to the worker computing modules in two layers. The number of threads in each module must match the model of the graphics card used. For example, the RTX 4090 graphics card is equipped with 10752 CUDA cores, and its theoretical parallel computing capability is equivalent to 336 times that of a 16-core 23-thread CPU, which can significantly improve the efficiency of integral calculation. In the third step, after the GPU completes the parallel calculation of the integral term and the spatial derivative term, the calculation result is sent back to the CPU host. In the fourth step, the continuity equation and momentum equation are coupled and solved on the CPU, and the calculation result is output and stored. The overall calculation process is shown in Figure 3.

[0071] Example 3 – Comparison with Stokes wave theory solution: This embodiment compares the wave equation with perfect dispersion obtained in Example 1 with the theoretical solution of a fifth-order Stokes wave to prove the model's perfect dispersion characteristics. Specifically, it includes the following steps: Step 1, Calculation model settings: The numerical terrain was set with a water depth of 0.4 m, an incident wave height of 0.04 m, and periods of 0.5 s and 1.0 s, respectively. The wave rise of the two models was compared with the theoretical solution of a fifth-order Stokes wave. A fifth-order Stokes wave was used as the incident wave at the incident boundary. In the numerical simulation, the spatial step size was taken as... =L / 40, time step is taken as... =T / 80.

[0072] Step 2: Simulate using a wave numerical model: Figure 4 A comparison between the wavefront rise of the calculated model and the theoretical solution is given in the figure. It can be seen from the figure that the calculated model is in good agreement with the theoretical solution. The dispersion of this model is completely accurate, so there is no error in the calculated wavelength even at high frequencies such as T=0.5s.

[0073] Example 4 – Simulation of waves crossing a submerged breakwater: The following simulations demonstrate the crossing of submerged breakwaters by both regular and irregular waves. To prove that this model is more accurate than the higher-order Boussinesq equations in simulating wave propagation in complex terrain, simulations of submerged breakwaters with different slopes were conducted. The specific steps include: Step 1, Calculation model settings: The experimental results are derived from Luth's (Luth HR, Klopman G, Kitou N. Kinematics of wavesbreaking partially on an offshore bar, LDV measurements of waves with and without a net onshore current. Report H-1573. Delft Hydraulics, Delft, The Netherlands. 1994.) physical model experiment, with the experimental setup as follows: Figure 5 As shown, the water depth at the bottom of the dike is 0.4m, the water depth at the top of the dike is 0.1m, the experimental wave conditions are a wave height of 0.02, a period of 2.02s, a slope of 1:20 in front of the dike, and a slope of 1:10 behind the dike. In the numerical simulation, the spatial step size... =0.02m, time step is =0.01s.

[0074] Step 2: Simulate using a wave numerical model: Figure 5 A comparison between the wavefront rise of the calculated model and the theoretical solution is given in the figure. It can be seen from the figure that the calculated model is in good agreement with the theoretical solution. The dispersion of this model is completely accurate, so there is no error in the calculated wavelength even at high frequencies such as T=0.5s.

[0075] Example 5 – Numerical Simulation of Rift Flow in a Dual-Channel Sandbar Topography: This embodiment conducts a systematic numerical simulation study on Haller's (Haller, Merrick C. Experimental study of nearshore dynamics on a barred beach with rip channels[J]. Journal of Geophysical Research Oceans, 2002, 107(C6):1-21.) classic double-groove sandbar topographic rift flow experiment. The study first verifies the applicability of the established model to nearshore strongly nonlinear dynamic processes such as rift flow by quantitatively comparing the numerical calculation results with the measured experimental results; at the same time, it supplements the comparison and analysis of the calculation results of this model with the results obtained by Haas (Haas KA, Warner J C. Comparing a quasi-3D to a full 3D nearshorecirculation model: SHORECIRC and ROMS[J]. Ocean Modelling, 2009, 26(1):91-103.) using the SHORECIRC model, further clarifying the improvement in calculation accuracy and the difference in results of the model in this paper compared with the traditional phase-averaged model.

[0076] Step 1, Calculation model settings: Haller's experimental terrain design is typical, consisting of a 1:5 steep slope section connected to a 1:30 gentle slope section. The 1:30 gentle slope area features long, narrow sandbars. On these sandbars, 3.6 m from each side of the shoreline, two 1.8 m long trenches are specifically designed. The complete experimental terrain is as follows: Figure 6 As shown. This numerical simulation uses the Test B standard example published by Haller, with clearly defined core parameters: water depth h = 0.698 m, incident regular wave height H = 5.12 cm, and wave period T = 1 s. During the numerical calculation, the spatial and time steps were reasonably set to balance accuracy and efficiency. The model was vertically discretized into a 5-layer structure, and the boundary conditions on both sides of the coast were non-slip solid boundary conditions. The key calculation results, such as the time-averaged flow field distribution, wave height variation, and water level increase / decrease presented below, are all taken from the surface calculation output data.

[0077] Step 2: Simulate using a wave numerical model: Figure 7 A comparison of the time-averaged flow fields for different models is presented. Figure 7 (e) and Figure 7 (f) are the numerical results of applying this model under different side mixing coefficients. Figure 7 (e) Using a side-mixing coefficient cm=2.0, Figure 7 (f) Using a side mixing coefficient cm=0.25, therefore Figure 7 The results in (e) show that the flow velocity at both fracture channels propagates positively toward the shore. Figure 7 The results in (f) show that the flow velocity at both split channels exhibits biased propagation, and this simulation result is closer to the experimental result.

[0078] Figure 8 visually presents the comparison results of wave height and the spatial distribution of wave increase and decrease in water volume between this model and the comparative model. (For...) Figure 8 (a) In the key area of ​​the split channel, the wave height value calculated by this model is higher than that of the SHORECIRC model, and the simulation results of the two models are in high agreement with the measured experimental data, with good overall fit. Figure 8 (b) In the multi-section comparison, the results calculated by this model are in better agreement with the experimental results. This advantage is particularly evident at the three typical sections of y=9 m, 11 m and 15 m. Figure 8 The comparison results in (c) are the opposite; the simulation results of the SHORECIRC model fit the experimental values ​​more closely, while the calculation results of this model show a slight deviation, indicating slightly insufficient simulation accuracy. Figure 9 shows a comparison of the distribution characteristics of time-averaged flow velocities U and V for different models. Figure 9 (a) shows that there are significant differences in the calculation results of the two models in the single-sided split channel region, specifically that the time-averaged velocity U simulated by the SHORECIRC model is greater than that of this model. Figure 9 In (b), the time-averaged velocity U spatial distribution simulated by this model is more uniform. In summary, the calculated results of this index by the SHORECIRC model are in better agreement with the experimental results. Figure 9 In (c), the time-averaged flow velocity calculation results of the two models are quite similar and both fit well with the experimental measured data, thus verifying the reliability of the simulation.

[0079] Since Haller's original experiment did not provide data on wave height, wave rise / fall, or flow rate in the direction perpendicular to the shoreline, nor did it provide experimental results on the vertical distribution of the time-averaged flow velocity U, it was impossible to directly verify the accuracy of the model's simulation of these elements through experimental data. Therefore, this embodiment further selects the calculation results obtained by Michaud et al. using the SYMPHONIE model as a comparison basis to corroborate the accuracy and rationality of the model's simulation of wave height, wave rise / fall, and flow velocity U. Figure 10 shows the distribution patterns of wave height and wave rise / fall in the direction perpendicular to the shoreline, clearly indicating that y=9.2 m is the center of the sandbar and y=4.6 m is the center of the fracture channel, providing a clear reference for interpreting the results. From Figure 10(a) and Figure 10 As can be seen from the comparison in (b), the calculation results of this model are basically consistent with those of the SYMPHONIE model, and the core laws are highly consistent. Figure 10 (c) and Figure 10 (d) reveals the differences between the two models: This model exhibits significant wave reduction in the x=12 m-14 m range, while the SYMPHONIE model fails to simulate this characteristic; in the x=15 m-19 m range, the wave increase calculated by this model is slightly greater than that of the SYMPHONIE model. These differences align with the core theory of wave motion, specifically: outside the break zone, wave height gradually increases before breaking, and the corresponding wave reduction increases synchronously with the wave height towards the coast, reaching its maximum at the break point; after entering the break zone, wave increase exhibits a linear change, rising towards the shore from the break point; when the wave propagates to the sandbar area, a wave recovery effect occurs, the wave height no longer decreases, and the corresponding wave increase / decrease remains stable; after secondary breaking, the wave height continues to decrease, and wave increase still maintains a linear growth trend. In summary, the numerical simulation results of this model are completely consistent with wave motion theory, supporting the rationality of the simulation.

[0080] The above are merely preferred embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for generating waves with precise dispersion, characterized in that, include: Based on the basic parameters of the ocean area to be simulated, determine the model calculation parameters; Based on the calculated parameters of the model, a numerical model of ocean waves with accurate dispersion is constructed, and ocean waves with accurate dispersion are generated using the numerical model of ocean waves. The construction of the ocean wave numerical model includes: based on the Euler equations and free surface boundary conditions, obtaining an integral-differential control equation set through continuous Taylor expansion and Fourier transform. The control equation set includes a continuity equation and a momentum equation, wherein the continuity equation introduces an integral-differential operator E to achieve a broadband wave description. According to the governing equations, physical process correction terms are embedded in the momentum equation, including underwater friction terms, turbulent mixing terms, and wave breaking energy dissipation terms. Based on the control equations with embedded physical process correction terms, the finite difference method is used for numerical discretization. The time discretization adopts the prediction-correction scheme, and the spatial discretization adopts the high-order difference scheme. The integral operation is accelerated by the CUDA-GPU heterogeneous parallel computing architecture to obtain the spatiotemporal distribution of free surface rise and realize the generation of waves with accurate dispersion.

2. The method for generating waves with precise dispersion according to claim 1, characterized in that, The basic parameters include water depth conditions, spectral peak frequency, spectral peak period, statistical characteristic wave height of incident waves, and incident azimuth angle of long-period waves. The model calculation parameters include computational domain range parameters, spatiotemporal resolution parameters, bottom friction coefficient, lateral mixing coefficient, and wave breaking parameters.

3. The method for generating waves with precise dispersion according to claim 1, characterized in that, The ocean wave numerical model with accurate dispersion is as follows: ； ； in, For the free surface to rise, For time, For the horizontal gradient operator, Because the local water is still and deep, It is the acceleration due to gravity. Let be the horizontal velocity vector at the free surface. For the integral-differential operator.

4. The method for generating waves with precise dispersion according to claim 1, characterized in that, The integral-differential operator for: ； in, k For wave number, h Because the local water is still and deep, For the horizontal gradient operator, x For the coordinates of the control points, Let i be the coordinates of the integration point, and i be the imaginary unit.

5. The method for generating waves with precise dispersion according to claim 4, characterized in that, According to the governing equations, embedding the seabed friction term into the momentum equation includes: ， ； in, The coefficient of friction at the bottom of the water is a constant. for x underwater friction in the direction of the water. for y underwater friction in the direction of the water. For the total water depth, for x Average speed in direction over time for y Average speed over time in the direction.

6. The method for generating waves with precise dispersion according to claim 5, characterized in that, According to the governing equations, embedding turbulent mixing terms into the momentum equation includes: ； ； in, Where is the eddy viscosity coefficient. for x Turbulent mixing forces in the direction of flow for y Turbulent mixing forces in the direction of flow.

7. The method for generating waves with precise dispersion according to claim 6, characterized in that, According to the governing equations, embedding wave breaking energy dissipation terms into the momentum equation includes: ； ； in, for x The direction of wave breaking force, for y The direction of wave breaking force, The value is the eddy viscosity coefficient caused by wave breaking.

8. The method for generating waves with precise dispersion according to claim 7, characterized in that, Based on the set of governing equations with embedded physical process correction terms, the numerical discretization solution using the finite difference method includes: The prediction step is calculated based on the fifth-order Adams-Bashforth scheme to obtain the prediction value for the next time level; The correction step is calculated using the sixth-order Adams-Moulton scheme to obtain the correction value for the next time layer. The convergence of the correction step is determined based on the preset error limit. When the number of iterations exceeds the threshold, the weighted average method is used to correct the two correction values ​​before and after to obtain the optimized correction value of the current time step. The spatial derivative term is discretized using the five-point central difference scheme, and the infinite integration interval is transformed into a finite interval for numerical integration based on the rapid decay characteristic of the kernel function.

Images