A method to consider rip current risk near multi-segment submerged dikes

Abstract

The present invention provides a method for considering the risk of rip currents near a multi-section submerged dike, including setting a calculation model, simulating the rip current system on the multi-section submerged dike using a simplified model of a nonlinear gentle slope equation accelerated by a GPU, and evaluating the risk of rip currents using the maximum value of the rip current velocity. The model described in the present invention is obtained by discarding high-order nonlinear terms from a fully nonlinear gentle slope equation model and performing GPU parallel accelerated calculation based on CUDA C, so the model has both accurate dispersion and higher computational efficiency. The model was applied to the simulation of rip current experiments on slotted sand dam terrain, proving the applicability and accuracy of the model for simulating rip current systems near submerged dikes.

Description

Technical Field

[0001] The invention belongs to the technical field of marine engineering hydrodynamics, and in particular relates to a method for considering the risk of rip currents near a multi-section submerged dike. Background Art

[0002] Submerged dikes are a commonly used coastal erosion protection structure, and multi-section submerged dikes are often used in coastal protection of tourist coastal beaches due to their unique economic, landscape and ecological advantages. Previous studies have shown that the groove areas between multi-section submerged dikes are high-risk areas for rip current disasters. These rip currents not only carry a large amount of sediment from the beach area to the open sea, but also pose a serious threat to the safety of tourists. At present, the consideration of rip current risks on multi-section submerged dikes is mostly only for fixed submerged dike structural forms, and no clear rip current risk level is given. For the simulation of rip current systems, most of them use wave phase average models or nonlinear shallow water equation models and Boussinesq equation models with low dispersion accuracy, which are insufficient in the simulation accuracy of wave and flow fields in areas with larger water depths. In addition, due to the large-scale characteristics of the marine environment, the computational domain is large, which often leads to a long simulation time, especially when applied to wave phase identification models. This problem is more prominent. The nonlinear mild slope equation is more suitable for rip current system simulation and rip current risk prediction near multi-section submerged dikes in deep offshore areas due to its accurate dispersion and small computational complexity.

[0003] Therefore, it is an issue that needs to be addressed to develop a nonlinear gentle slope equation model that is applicable to greater water depths, has higher calculation accuracy, and faster calculation speed, so as to more systematically and accurately consider the risk of rip currents near multi-segment submerged dikes. Summary of the invention

[0004] In view of the problems existing in the prior art, the purpose of the present invention is to provide a method for considering the risk of rip current near a multi-segment submerged dike. The method applies a simplified model of a time-domain nonlinear gentle slope equation based on GPU parallel acceleration of CUDAC to the simulation study of the rip current system on the multi-segment submerged dike, and evaluates the risk of rip current near the multi-segment submerged dike according to the maximum rip current velocity and the rip current area distribution based on the obtained time-averaged flow field results. The model described in the invention is obtained by discarding high-order nonlinear terms of the fully nonlinear gentle slope equation model and performing GPU parallel acceleration calculation based on CUDAC, so the model has both accurate dispersion and higher computational efficiency. The model is applied to the simulation of the rip current system on the slotted sand bar terrain. The results show that the current model is more consistent with the wave height, water increase and decrease and flow velocity results in the experiment than the commonly used low-order Boussinesq equation model, which proves the applicability and accuracy of the model for the simulation of the rip current system near the submerged dike. The model is used to simulate the rip current near the multi-segment submerged dike under four incident wave conditions: weak waves, normal waves, storm waves and oblique waves, and the influence of wave energy and wave incident angle on the rip current risk level is analyzed. By considering the influence of different submerged dike widths, profile forms, groove widths and groove numbers in the multi-section submerged dike structure, the influence of the lateral and longitudinal spatial configuration of the multi-section submerged dike on the rip current risk level near the submerged dike is explained.

[0005] The objective of the present invention is achieved through the following technical solutions:

[0006] The present invention provides a method for considering the risk of rip currents near a multi-segment submerged dike, comprising the following steps:

[0007] Step 1, calculation model settings:

[0008] First, determine the multi-section submerged dike configuration and wave field data on the beach to be studied;

[0009] Then, the calculation parameters of the model are set, including wave breaking parameters, bottom friction coefficient, lateral mixing coefficient, calculation space and time step, and total calculation time;

[0010] Step 2: Use the GPU-accelerated simplified model of the nonlinear gentle slope equation to simulate the rip current system on a multi-segment submerged dike:

[0011] Under the model settings given in step 1, the GPU-accelerated nonlinear mild slope equation simplified model is used for simulation calculation, including:

[0012] S21, using a Cartesian rectangular coordinate system (x, y, z), where the x-axis is perpendicular to the shoreline direction, the y-axis is parallel to the shoreline direction, and the z-axis is vertically upward. In this coordinate system, the GPU-accelerated nonlinear gentle slope equation simplified model has the following control equation:

[0013]

[0014] In the formula, t is time, and the subscript t represents the derivative with respect to time. is the symbol of the gradient operator in the horizontal direction, U = (U, V) is the horizontal velocity vector at the still water level, U and V are the flow velocities in the x and y directions on the still water surface, η is the wave surface, h is the still water depth, g is the acceleration of gravity, k is the wave number, T is the wave period, c is the wave speed, the coefficient n = 1 / 2 + kh / sinh (2kh), and R is the term that takes into account wave breaking, bottom friction and lateral mixing effects;

[0015] S22, the above equations (1)-(2) are applied to the simulation of rip current system near multi-section submerged dike under the condition of oblique wave incidence. The problem of wave surface and water flow value at the lateral boundary in open water area should be considered. The lateral cyclic calculation boundary condition is used to deal with this problem, that is, three layers of virtual grids are added outside the calculation domain, and the wave surface and flow velocity results on the virtual grid nodes along the coast are cyclically assigned and then solved;

[0016] S23, for the above equations (1)-(2), in order to meet the timeliness of rip current risk prediction, the model calculation efficiency should be improved, and the unified device architecture CUDA and C language are used to perform GPU parallel calculation on the model to improve the efficiency of rip current risk prediction;

[0017] Step 3: Use the maximum rip flow velocity to assess the rip flow risk:

[0018] Using the simulation results of the water velocity distribution on the multi-stage submerged dike given in step 2, the time-averaged and depth-averaged rip current velocity is calculated, and the rip current risk level is evaluated based on the maximum rip current velocity obtained, including:

[0019] S31, when calculating the time-averaged and depth-averaged rip current velocity, it is necessary to consider 1) converting the velocity at the current static water level into the depth-averaged velocity and 2) defining the rip current velocity using the depth-averaged velocity;

[0020] For problem 1), the flow rate is converted using the following formula using the flow conservation principle:

[0021]

[0022] Where <·> represents the average value over time, u=(u,v), u and v are the average water depth velocities in the x and y directions, respectively;

[0023] For question 2), since the danger of rip current includes both the velocity in the offshore direction and the velocity in the coastal direction, the sum of the velocity vectors u and v is used as the basis for judging the magnitude of the rip current velocity. In addition, in order to exclude the influence of the velocity toward the shore, only the positive part of the velocity u perpendicular to the shore, that is, the offshore part, is taken; therefore, the time-averaged and depth-averaged rip current velocity value is:

[0024]

[0025] In the formula, u + It means that the velocity u takes a positive value;

[0026] S32, according to the full-field rip current velocity value results given by the above formula (11), select the maximum rip current velocity value |u| max ,The rip current risk levels are divided to quantitatively assess the risk of rip currents near multi-segment submerged dikes: |u| max <0.3m / s is considered as no risk of rip flow, |u| max Offshore currents ≥0.3m / s are considered rip current hazards. Specifically, 0.3m / s≤|u| max When <0.8m / s, it is considered as low risk of rip current, 0.8m / s≤|u| max When <1.3m / s, it is considered as medium risk of rip current, |u| max When the flow rate is ≥1.3m / s, it is considered as a high risk of rip current.

[0027] Furthermore, in step 1, the multi-section submerged breakwater configuration and wave field data on the beach to be studied include the scope of the beach to be studied, the still water depth of the seabed topography, the profile type, width, length in the coastal direction and the size of the submerged breakwater spacing in the offshore direction, and the wave height, period and incident angle of the incident wave.

[0028] Furthermore, the specific solution process of S22 in step 2 is:

[0029] Consider that the total number of grids in the computational domain is M in the x direction and N in the y direction, where the value range of the grid nodes (i, j) in the computational domain is i = 1 to M and j = 1 to N, and the value of the virtual grid node (i, j) is i = 0, -1, -2 at the shore boundary of the computational domain; i = M + 1, M + 2, M + 3 at the offshore boundary; j = 0, -1, -2 at the left boundary; and j = N + 1, N + 2, N + 3 at the right boundary;

[0030] Therefore, the values ​​of wave surface η, flow velocity U and V on the virtual grid nodes on the left boundary are as follows:

[0031] η(i,1-j)=η(i,N+1-j),j=1,2,3 (3)

[0032] U(i,1-j)=U(i,N+1-j),j=1,2,3 (4)

[0033] V(i,1-j)=V(i,N+1-j),j=1,2,3 (5)

[0034] On the right boundary, the values ​​of the wave surface η, flow velocities U and V in the virtual grid are as follows:

[0035] η(i,N+j)=η(i,j),j=1,2,3 (6)

[0036] U(i,N+j)=U(i,j),j=1,2 ,3 (7)

[0037] V(i,N+j)=V(i,j),j=1,2 ,3 (8)

[0038] When the incident angle α of the offshore wave is known and the number of computational grids N in the y direction along the coast is fixed, in order to meet the lateral circulation boundary conditions, the value of the computational space step length △y in the coastal direction needs to be adjusted so that the coastal direction contains an integer number of oblique incident waves. Therefore, the value of △y should satisfy the following relationship:

[0039]

[0040] Where, L w =2π / k is the deep water wavelength, M w is an adjustable positive integer parameter;

[0041] After obtaining △y according to the above formula, readjust the calculation domain width to (N-1)△y.

[0042] Furthermore, in step (2), S23, the process of performing GPU parallel computing on the model based on CUDA and C language includes: declaration and definition of variables on the host side and the device side, allocation of memory space for variables on the host side and the device side, initialization of variables on the host side, copying the host side data to the device side, calling the kernel functions on the device side in sequence to complete the computing, copying the device side data back to the host side, and finally outputting the results;

[0043] In the process of calling the kernel function, the number of logically enabled threads is organized at two levels. The first level is to specify the grid size, that is, the number of thread blocks contained in the grid in the x and y directions; the second level is to define the size of the thread block, that is, the number of threads contained in a thread block.

[0044] Compared with the prior art, the present invention has the following beneficial effects:

[0045] 1. In the method of considering the risk of rip currents near multi-section submerged dikes described in the present invention, the simplified model of nonlinear gentle slope equation accelerated by GPU parallelization is applied for the first time to the simulation of rip current system on multi-section submerged dike beaches. Since the model described in the present invention is a nonlinear gentle slope equation type model, it has the characteristics of accurate dispersion and can be better applied to the simulation of wave surface and water flow near multi-section submerged dikes far offshore. The model was applied to the simulation of rip current experiments on slotted sand dam coasts, and the results also showed that the current model has better performance than the weakly dispersive Boussinesq equation model.

[0046] 2. Since the model of the present invention removes the high-order nonlinear terms and the weak nonlinear terms containing the bottom slope in the completely nonlinear gentle slope equation model, the model form is simpler; in addition, since CUDA C is used for GPU parallel computing, the current model has high computational efficiency, ensuring the timeliness of large-scale real coastal rip current risk prediction and assessment.

[0047] 3. The model described in the present invention uses a lateral circulation boundary to simulate the lateral wave surface and flow velocity input results on an infinitely long open coast, which can better consider the real marine environment compared with the traditional solid boundary situation.

[0048] 4. The present invention classifies rip current risks according to the maximum rip current velocity, and clarifies the risk level of rip current disasters near multi-section submerged dikes.

[0049] 5. The present invention provides the influence of incident wave energy, incident angle and transverse and longitudinal arrangement of the submerged dike on the rip current risk level, and reveals the main influencing factors of the rip current risk level near the multi-section submerged dike. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] The present invention will be further described below in conjunction with the accompanying drawings and embodiments:

[0051] Figure 1 is the instantaneous wave surface on the multi-section submerged dike when the wave incident angle is α=30° in Example 1;

[0052] Figure 2 This is a schematic diagram of the calculation process of the GPU parallel acceleration model developed based on CUDAC in Example 2;

[0053] Figure 3 is the computation time of the GPU parallel model and the CPU serial model in Example 2 and the speedup ratio of the GPU parallel model;

[0054] Figure 4 Schematic diagram of the plan layout (a) and vertical shoreline section (b) of the Haller experiment;

[0055] Figure 5The distribution of wave height H (left column) and average water surface <η> (right column) at three sections of y = 1.1 m, y = 4.6 m and y = 7.7 m in the Haller experiment; where ○: experimental value; ---: Wang et al. model; —: the model of the present invention;

[0056] Figure 6 is the vertical shore velocity at the four sections x = 10m, x = 11.25m, x = 12.3m and x = 13m in the Haller experiment (left column) and coastal speed <v>(right column); where ○: experimental value; ---: Wang et al. model; —: the model of the present invention;

[0057] Figure 7 The calculation domain plane layout (a), the calculation domain cross-sectional shape (b) and the Gaussian type submerged dike cross-sectional shape (c) of the multi-stage submerged dike beach in Example 3;

[0058] Figure 8 : are the distribution results of wave height (left column) and average water surface (right column) under different wave conditions in Example 3; wherein the dotted line represents the boundary of the multi-segment submerged dike;

[0059] Fig. 9 : is the time-averaged velocity (→) and rip current area (color) distribution results under different wave conditions in Example 3; wherein the dotted line represents the multi-stage submerged dike boundary;

[0060] Fig.10 : is the time-averaged flow velocity (→) and rip current area (color) distribution results at different submerged dike widths and heights in Example 4; wherein the dotted line represents the multi-segment submerged dike boundary;

[0061] Fig.11 is a schematic diagram of different submerged dike profiles used in Example 4;

[0062] Fig.12 : is the time-averaged flow velocity (→) and rip current area (color) distribution results under different submerged dike profile shapes in Example 4; wherein the dotted line represents the multi-segment submerged dike boundary;

[0063] Fig.13 : is the time-averaged flow velocity (→) and rip current area (color) distribution results under different submerged dike groove widths in Example 5; wherein the dotted line represents the multi-segment submerged dike boundary;

[0064] Fig.14 1 is the time-averaged flow velocity (→) and rip current area (color) distribution results for different numbers of submerged dike grooves in Example 5; wherein the dotted line represents the multi-segment submerged dike boundary. DETAILED DESCRIPTION

[0065] The preferred embodiments of the present invention will be described in detail below in conjunction with examples. It should be understood that the following examples are provided only for the purpose of illustration and are not intended to limit the scope of the present invention. Those skilled in the art may make various modifications and substitutions to the present invention without departing from the purpose and spirit of the present invention. If not otherwise specified, the technical means used in the examples are conventional means well known to those skilled in the art, and the equipment and reagents used are all conventionally commercially available.

[0066] Example 1

[0067] This embodiment provides a method for considering the risk of rip currents near a multi-segment submerged dike, including the following steps:

[0068] Step 1, calculation model settings:

[0069] First, determine the multi-section submerged breakwater configuration and wave field data on the beach to be studied, including the scope of the beach to be studied, the seabed topography and still water depth, the offshore profile type, width, along-shore length and spacing of the submerged breakwaters, and the wave height, period and angle of incidence of the incident waves.

[0070] Then, the calculation parameters of the model are set, including wave breaking parameters, bottom friction coefficient, lateral mixing coefficient, calculation space and time steps, and total calculation time.

[0071] Step 2: Use the GPU-accelerated simplified model of the nonlinear gentle slope equation to simulate the rip current system on a multi-segment submerged dike:

[0072] Under the model settings given in step 1, the GPU-accelerated nonlinear mild slope equation simplified model is used for simulation calculation, including:

[0073] S21, using a Cartesian rectangular coordinate system (x, y, z), where the x-axis is perpendicular to the shoreline direction, the y-axis is parallel to the shoreline direction, and the z-axis is vertically upward. In this coordinate system, the GPU-accelerated nonlinear gentle slope equation simplified model has the following control equation:

[0074]

[0075] In the formula, t is time, and the subscript t represents the derivative with respect to time. is the symbol of the gradient operator in the horizontal direction, U = (U, V) is the horizontal velocity vector at the still water level, U and V are the flow velocities in the x and y directions on the still water surface, η is the wave surface, h is the still water depth, g = 9.81 is the gravitational acceleration, k is the wave number, which can be obtained from the accurate dispersion relation ω 2 =gktanhkh for solution, where ω=2π / T is the wave circular frequency, T is the wave period, c=ω / k is the wave speed, the coefficient n=1 / 2+kh / sinh(2kh), and R is the term considering wave breaking, bottom friction and lateral mixing. Its definition is the same as the result given by Xu et al. (Xu J, Zou ZL, Yan SA Set of Fully Nonlinear Mild Slope Equations[J]. Ocean engineering. 2024, 297: 116881.), which will not be repeated here.

[0076] The above equations (1)-(2) are simplified by removing the high-order nonlinear terms and the weak nonlinear terms containing the bottom slope from the fully nonlinear mild slope equation model of Xu et al., so it has the precise dispersion performance unique to mild slope equations, and because its form is simpler, its computational efficiency is higher.

[0077] S22, the above equations (1)-(2) are applied to the simulation of rip current system near multi-section submerged dike under the condition of oblique wave incidence. The problem of wave surface and water flow at the lateral boundary in open water area should be considered. The lateral cyclic calculation boundary condition is used to deal with this problem, that is, three layers of virtual grids are added outside the calculation domain, and the wave surface and flow velocity results on the virtual grid nodes along the coast are cyclically assigned and then solved.

[0078] The specific solution process is:

[0079] Consider that the total number of grids in the computational domain is M in the x direction and N in the y direction, where the value range of the grid nodes (i, j) in the computational domain is i=1~M and j=1~N, and the value of the virtual grid node (i, j) is i=0, -1, -2 at the shore boundary of the computational domain; i=M+1, M+2, M+3 at the offshore boundary; j=0, -1, -2 at the left boundary; and j=N+1, N+2, N+3 at the right boundary.

[0080] Therefore, the values ​​of wave surface η, flow velocity U and V on the virtual grid nodes on the left boundary are as follows:

[0081] η(i,1-j)=η(i,N+1-j), j=1,2,3 (3)

[0082] U(i,1-j)=U(i,N+1-j), j=1,2,3 (4)

[0083] V(i,1-j)=V(i,N+1-j), j=1,2,3 (5)

[0084] On the right boundary, the values ​​of the wave surface η, flow velocities U and V in the virtual grid are as follows:

[0085] η(i,N+j)=η(i,j), j=1,2 ,3 (6)

[0086] U(i,N+j)=U(i,j), j=1,2 ,3 (7)

[0087] V(i,N+j)=V(i,j), j=1,2 ,3 (8)

[0088] When the incident angle α of the offshore wave is known and the number of computational grids N in the y direction along the coast is fixed, in order to meet the lateral circulation boundary conditions, the value of the computational space step length △y in the coastal direction needs to be adjusted so that the coastal direction contains an integer number of oblique incident waves. Therefore, the value of △y should satisfy the following relationship:

[0089]

[0090] Where, L w =2π / k is the deep water wavelength, M w is an adjustable positive integer parameter;

[0091] After obtaining △y according to the above formula, readjust the calculation domain width to (N-1)△y.

[0092] Figure 1 Taking the wave incident angle α=30° as an example, the wave surface results on the multi-segment submerged breakwater under the lateral circulation boundary condition are given. It can be seen that the waves on both sides have good consistency, which verifies the feasibility of using lateral circulation boundaries to simulate open sea areas.

[0093] S23, for the above equations (1)-(2), in order to meet the timeliness of rip current risk prediction, it is necessary to improve the model calculation efficiency. The unified device architecture CUDA and C language are used to perform GPU parallel calculation on the model to improve the efficiency of rip current risk prediction.

[0094] The development platform used is Windows 10+Intel Complier Visual Studio 2017, the language is C++, the CUDA version used is CUDA 10.2, and the CPU and GPU models and parameters used are as follows:

[0095] CPU: Intel i3-9100f, 4 cores and 4 threads, main frequency is 3.6GHz, single-core turbo frequency is 4.2GHz, and multi-core turbo frequency is 4.0GHz.

[0096] GPU: NVIDIA GeForce RTX 2060 graphics card, Turing 104 architecture, 1920 CUDA cores, 1755MHz core frequency, 6GB video memory capacity, GDDR6 video memory type, 192bit video memory bus width, 336GB / s bandwidth, and 1.4GHz video memory frequency.

[0097] like Figure 2 As shown in the figure, the process of GPU parallel computing of the model based on CUDA and C language includes: declaration and definition of host and device variables, allocation of host and device variable memory space, initialization of host variables, copying host data to device, calling device kernel functions in sequence to complete the calculation, copying device data back to host, and finally outputting the results.

[0098] In the process of calling the kernel function, the number of logically enabled threads is organized at two levels. The first level is to specify the grid size, that is, the number of thread blocks contained in the grid in the x and y directions; the second level is to define the size of the thread block, that is, the number of threads contained in a thread block.

[0099] In this model, based on the above development environment, the thread block size is selected as 16×16 after testing, and the size of the grid can be freely adjusted according to the number of grids contained in the computational domain.

[0100] Step 3: Use the maximum rip flow velocity to assess the rip flow risk:

[0101] Using the simulation results of the water velocity distribution on the multi-stage submerged dike given in step 2, the time-averaged and depth-averaged rip current velocity is calculated, and the rip current risk level is evaluated based on the maximum rip current velocity obtained, including:

[0102] S31, when calculating the time-averaged and depth-averaged rip current velocity, it is necessary to consider 1) converting the velocity at the current static water level into the depth-averaged velocity and 2) defining the rip current velocity using the depth-averaged velocity;

[0103] For problem 1), the flow rate is converted using the following formula using the flow conservation principle:

[0104]

[0105] Where <·> represents the average value over time, u=(u,v), u and v are the average water depth velocities in the x and y directions, respectively;

[0106] For question 2), since the danger of rip current includes both the velocity in the offshore direction and the velocity in the coastal direction, the sum of the velocity vectors u and v is used as the basis for judging the magnitude of the rip current velocity. In addition, in order to exclude the influence of the velocity toward the shore, only the positive part of the velocity u perpendicular to the shore, that is, the offshore part, is taken; therefore, the time-averaged and depth-averaged rip current velocity value is:

[0107]

[0108] In the formula, u + It means that the velocity u takes a positive value.

[0109] S32, according to the full-field rip current velocity value results given by the above formula (11), select the maximum rip current velocity value |u| max , and the rip current risk levels are divided according to Table 1 to quantitatively assess the risk of rip currents near multi-segment submerged dikes: |u| max <0.3m / s is considered as no risk of rip flow, |u| max Offshore currents ≥0.3m / s are considered rip current hazards. Specifically, 0.3m / s≤|u| max When <0.8m / s, it is considered as low risk of rip current, 0.8m / s≤|u| max When <1.3m / s, it is considered as medium risk of rip current, |u| max When the flow rate is ≥1.3m / s, it is considered as a high risk of rip current.

[0110] Table 1 Classification of rip current risk levels

[0111]

[0112] Example 2 - GPU parallel acceleration calculation of simplified model of nonlinear gentle slope equation

[0113] This embodiment provides a brief process of GPU parallel accelerated computing of a simplified model of a nonlinear mild slope equation based on CUDAC, and evaluates the GPU parallel accelerated computing performance.

[0114] (1) Introduction to GPU parallel computing process

[0115] The development platform used in this embodiment is Windows 10+Intel ComplierVisual Studio 2017, the language is C++, the CUDA version used is CUDA 10.2, and the CPU and GPU models and parameters used are as follows:

[0116] CPU: Intel i3-9100f, 4 cores and 4 threads, main frequency is 3.6GHz, single-core turbo frequency is 4.2GHz, and multi-core turbo frequency is 4.0GHz.

[0117] GPU: NVIDIA GeForce RTX 2060 graphics card, Turing 104 architecture, 1920 CUDA cores, 1755MHz core frequency, 6GB video memory capacity, GDDR6 video memory type, 192bit video memory bus width, 336GB / s bandwidth, and 1.4GHz video memory frequency.

[0118] The present invention is based on the GPU parallel model developed by CUDAC, which roughly follows Figure 2 The calculation process is shown in Figure 5. The most important step in this process, which is also the main difference from the CPU serial model, is step 5 - calling the kernel function on the device side to complete the calculation. The kernel function is an important concept in CUDA. It is executed by the CUDA core on the device side for matrix operations and is called on the host side. In the process of calling the kernel function, the number of logically enabled threads needs to be organized according to two levels. The first level is to specify the size of the grid, that is, the number of thread blocks contained in the grid in the x and y directions; the second level is to define the size of the thread block, that is, the number of threads contained in a thread block. In this model, based on the above development environment, the thread block size is selected as 16×16 after testing, and the size of the grid can be freely adjusted according to the number of grids contained in the calculation domain.

[0119] (2) GPU parallel computing acceleration performance analysis

[0120] In order to verify the acceleration effect of the current parallel computing method on the simplified model of the nonlinear mild slope equation, for the same example with the same computational amount, the GPU model and CPU model are used to test the terminal time consumed by simulating the same wave propagation time to calculate the acceleration ratio of the GPU model. The specific method is as follows.

[0121] The above equipment is used to simulate the example of wave propagation on a flat coast. The simulation time is 100 seconds. The terminal calculation time T of the GPU and CPU models are calculated respectively. GPU and T CPU , and calculate the speedup ratio A r =T CPU / T CPU The value of .

[0122] In order to study the effect of the number of grids on the speedup ratio, a horizontal two-dimensional wave propagation simulation example is designed to keep the number of grids in the x direction unchanged (M = 801), and to study it by changing the number of grids in the y direction N. The calculation process does not involve the breaking and climbing of waves, and does not consider the effects of bottom friction and lateral mixing. Figure 3 The calculation time of the two models with different numbers of grids and the speedup ratio of the GPU model relative to the CPU model are given for the GPU and CPU models. As can be seen from the figure, the calculation efficiency of the GPU model based on CUDA programming is significantly higher than that of the CPU model, especially when the number of grids is large, the calculation efficiency is improved by about 23 times. When the number of grids is less than 100,000, the speedup ratio of the GPU model is A r It increases with the number of grids until it reaches a maximum value of about 23; when the number of grids is greater than 100,000, the acceleration ratio of the GPU model basically remains unchanged at the maximum value.

[0123] The above results show that the simplified model of nonlinear mild slope equations based on CUDAC GPU parallel acceleration has a significant improvement in computational efficiency compared with the CPU serial model. When the computational domain area is large enough, its computational efficiency can be improved by more than 20 times. The computational efficiency of the model described in the present invention can meet the needs of rip current prediction and risk assessment near large-scale multi-stage submerged dikes.

[0124] Comparative example: simulation of rip current experiment on a slotted sand bar coast

[0125] In order to verify the accuracy of the model described in the current invention in simulating beach rip currents, the model was used to simulate the rip current experiment on the coast of the pool slot sand bar implemented by Haller (Haller MC, Dalrymple RA, Svendsen IA. Experimental Study of Nearshore Dynamics on a Barred Beach with RipChannels[J]. Journal of Geophysical Research. 2002(107):1–21.). This experiment is widely used in the verification of mathematical models, such as Wang et al. (Wang H, Zhu S, Li X, Zhang W, Nie Y. Numerical Simulations of Rip Currents off Arc-Shaped Coastlines[J]. Acta Oceanologica Sinica. 2018(37):21–30.) The experiment was simulated using the FUNWAVE model based on the Boussinesq equation.

[0126] The plan and section layout of the pool in the Haller experiment above is shown in Figure 4 . The pool is 17.2m long and 18.2m wide. There is a wave maker on the left side of the pool, and a flat-slope coast model with a slope of 1:5 and 1:30 made of concrete on the right side. A sand bar model with a width of 1.2m is superimposed at a position 9m away from the foot of the slope. The sand bar is cut into a three-section sand bar shape by two grooves with a width of 1.8m. The lengths of each section of the sand bar are 3.6m, 7.4m and 3.6m respectively. This comparative example only simulates the Group B test of the Haller experiment, in which the flat bottom water depth is 0.7m, and the incident wave height and period are 5.12cm and 1.0s respectively. In the experiment, the free surface rise on three sections perpendicular to the coast at y = 1.1m, y = 4.6m and y = 7.7m was measured to obtain the wave height H and average water surface <η> data on these sections; the flow velocity on four sections parallel to the coast at x = 10m, x = 11.25m, x = 12.3m and x = 13m was measured to obtain the time-averaged flow velocity near the sandbar slot area. and <v>result.

[0127] In the simulation process of this experiment, the wave breaking strength index is taken as σ = 2.0, and the wave breaking start parameter is taken as e I =0.4, the wave breaking termination parameter is taken as e F =0.05, the friction coefficient of the bottom is taken as f w =0.008, the lateral mixing coefficient is C m =0.25. The spatial steps in the x and y directions are △x=0.05m and △y=0.1m respectively, and the time step is △t=0.02s. Since the waves in this comparative example are incident in the normal direction, there is no need to recalibrate the value of △y in the coastal direction. The total simulation time is 300s, of which the simulation results from 30 to 300s are used to calculate the wave height H, the average water surface <η>, and the time-averaged flow velocity and <v>.

[0128] Figure 5 The calculated wave height H and average water surface <η> results of this comparative example at three sections of y=1.1m, y=4.6m and y=7.7m are compared with the experimental results and the results given by Wang et al. Figure 6 The time-averaged flow velocity calculated for this comparative example at four sections x = 10m, x = 11.25m, x = 12.3m and x = 13m is given. and <v>The results of the present invention are compared with the experimental results and the results given by Wang et al. It can be seen from the figure that the results obtained by using the model described in the present invention are generally more consistent with the experimental results than the results given by Wang et al. Figure 5 As can be seen from b, d, and f, the current model obtains a more accurate average water surface result, which is also the main driving force of the flow field. Therefore, in terms of flow field, the current model also gives a more accurate rip current velocity amplitude. For example, at x=12.3m, the maximum rip current measured in the experiment was 15.5cm / s, located at y=4.6m along the coast; the maximum rip current on this section given by Wang et al. was only 8.5cm / s, located at y=4.2m, and the rip current amplitude error reached about 45%; while the maximum rip current given by the current model is 15.9cm / s, located at y=4.8m, and the rip current amplitude error is only 2.6%. In addition, the current model has a more accurate prediction of the coastal flow. <v>The prediction results are also better than those given by Wang et al.

[0129] In summary, the simplified model of the nonlinear gentle slope equation for GPU parallel acceleration based on CUDAC in the present invention benefits from its precise dispersion characteristics and can provide more accurate rip current velocity simulation results than low-order dispersion accuracy models such as the Boussinesq equation model, and is more suitable for the prediction and assessment of rip current risks near multi-stage submerged dikes on beaches.

[0130] Example 3 - Simulation of rip current system near multi-stage submerged dike under different incident wave conditions

[0131] This embodiment uses the method described in the current invention to simulate the rip current system near the multi-section submerged dike under different incident wave conditions. First, the calculation domain setting in the embodiment is given, and then four typical wave conditions are introduced to simulate the rip current system generation process. Then, the wave field characteristics are displayed in combination with the free surface simulation results. Finally, the rip current risk near the multi-section submerged dike is evaluated based on the distribution characteristics of the water velocity field.

[0132] (1) Computational domain settings

[0133] In this embodiment, the calculation domain is considered to have a coastal width of W d = 640m, offshore length is L d =560m, such as Figure 7 As shown. The upper side of the computational domain is the wave-making area, and the lower side is a flat coast model with a slope of 1:40. The still water depth at the flat bottom is set to h0 = 10m, then the still water line is located 400m away from the foot of the slope, and there is a 48m long wave climbing area above the still water line. During the simulation process, the coordinate origin is selected at the intersection of the left boundary and the still water line, the x-axis is located at the left boundary pointing to the offshore direction, and the y-axis is located at the still water line pointing to the coastal direction. On the flat coast, the distance from the still water line X B = 80m superimposed with a Gaussian wide-top and impermeable submerged dike (see Figure 7 In (c), the width of the submerged dike is W B =32m, top height is H B =1.28m, then the water depth at the top of the embankment is h c =1.28m. The height h of the Gaussian submerged dike B Determined by the following formula:

[0134] h B =H B exp[-(xX B ) 2 / δ] (12)

[0135] Where δ is a parameter for determining the width of the submerged dike. In this embodiment, the width of the submerged dike is W. B =32m, δ=96. In order to obtain a multi-section submerged dike, a section with a width of G is cut in the middle of the submerged dike. B = 32m trench, the length of each side of the submerged dike is 304m. Due to the cyclic lateral boundary effect set in the model, the actual length of each section of the multi-section submerged dike is L B =608m, the distance between submerged dikes, i.e. the width of the trench, is G B =32m.

[0136] (2) Incident wave conditions

[0137] Consider the normal wave condition with normal incidence as the typical incident wave condition, that is, the characteristic incident wave height and period at a flat bottom water depth are H i =0.6m and T = 6.0s. Based on this typical incident wave condition, by reducing the incident wave height to H i =0.4m to obtain the weak wave incident wave condition. By increasing the incident wave height to H i =1.6m and wave period T=12.0s to obtain storm wave incident wave conditions, and oblique wave incident wave conditions are obtained by changing the wave incident angle to 31.97°. The above four groups of characteristic test wave conditions are shown in Table 2. This embodiment studies the influence of wave energy and incident angle on rip current risk by considering the four groups of incident wave conditions.

[0138] Table 2 Four groups of characteristic incident wave conditions

[0139]

[0140] In the simulation process of this embodiment, the wave breaking parameters, bottom friction coefficient and lateral mixing coefficient are uniformly taken and consistent with the values ​​in the comparative example. In the calculation, the spatial step lengths in the x and y directions are taken as △x=0.8m and △y=1.6m respectively, and the time step length is taken as △t=0.08s. For the case of oblique incident waves in this embodiment, there is no need to adjust the value of △y under the lateral circulation boundary condition. This is because the circulation boundary condition (Formula (9)) can be automatically satisfied when the wave incident angle is 31.97°. The total duration of the simulation is 20 minutes, of which the simulation results of 2 to 20 minutes are used to calculate the wave height H, the average water surface <η>, and the time-averaged flow velocity and <v>.

[0141] (3) Wave field characteristics

[0142] Figure 8 The wave height and average water surface distribution results under the above four wave conditions are given. Since rip currents mainly occur near the groove area, the figure only shows the results within the range of y=160~480m and x=0~160m. The wave field distribution characteristics under different wave conditions can be seen from the figure. For wave condition A where the positive incident wave is a weak wave, the waves in the submerged dike area mainly break at the top of the dike, while the waves in the groove area break with delay, and due to the effect of wave diffraction, the distribution of breaking waves in this area presents an inverted Y shape; under the influence of the breaking wave, the coastal wave water increase gradient is mainly located in the area close to the coastline. For typical wave condition B where the incident wave is a normal wave, due to the larger incident wave height compared to wave condition A, the wave breaking positions in the groove area and the submerged dike area are both farther from the shore, and the coastal wave water increase gradient is mainly located between the submerged dike and the coastline, but the submerged dike still controls the occurrence of wave breaking. As the incident wave height continues to increase to the storm wave situation, that is, for wave condition C, the control of the submerged dike on wave breaking basically disappears, and the waves mainly break on the seaward side of the submerged dike. The submerged dike only affects the waves in the breaking process, so at this time the coastal water increase gradient is mainly located on both sides of the groove. For oblique incident waves with incident wave heights of normal waves, that is, wave condition D, due to the shielding effect of the submerged dike on the waves in the groove area, the wave breaking delay and coastal water increase gradient caused by the groove mainly occur on the downstream side of the groove. Since the wave height and coastal water increase gradient drive the generation of water flow in the form of radiation stress and water surface pressure respectively, the wave height and average water surface distribution corresponding to the above four different wave conditions will also drive different rip current morphologies and rip current disaster risks.

[0143] (4) Rip current risk assessment

[0144] Fig. 9 The water velocity field and rip current risk area corresponding to the above four wave conditions are given. It can be seen from the figure that the rip current system can be divided into three parts in the offshore direction: the rip current compensation flow between the submerged dike and the shoreline, the rip current neck in the groove area, and the rip current head on the seaward side of the submerged dike. The maximum rip current mainly occurs in the rip current neck area in the groove, which is determined by the above wave field distribution characteristics. Table 3 gives the maximum rip current velocity |u| corresponding to the four wave conditions. max , Offshore location X off As well as the rip current risk level, it can be seen from the table that for weak wave condition A, the maximum rip current velocity is only 0.57m / s, corresponding to the low rip current risk level. At this time, the rip current area is also small and only concentrated in the groove area. Fig. 9 For the normal wave incident wave condition B, due to the larger incident wave height, the maximum rip current velocity increases to 1.07 m / s, corresponding to the medium rip current risk level. At this time, the distribution area of ​​rip current risk also increases to the area including the rip current compensation flow and the rip current head, see Fig. 9 For storm wave condition C, due to the influence of large waves, the maximum rip current velocity in the groove area increases to 1.53 m / s, corresponding to a high rip current risk level. At this time, the rip current occurrence area includes not only the larger rip current compensation flow and rip current head area, but also the area on the top of the submerged dike. The rip current in this part is actually the bottom return flow caused by the large mass transfer flow, see Fig. 9 From the results in the table, it is also noted that oblique wave condition D also caused a high rip current risk level, and its maximum rip current velocity also reached 1.44m / s. The rip current risk area at this time is different from that of the normal incidence condition, mainly including the dam crest area, the groove area, and the seaward side of the submerged dike downstream of the groove, as shown in Fig. 9 This is because when waves are incident obliquely, the strong longshore current driven by wave radiation stress merges into the offshore current, forming a strong rip current risk area.

[0145] The results of this example show that the risk of rip currents in multi-section submerged dikes is mainly controlled by the maximum rip current velocity in the groove area, and the rip current risk level increases with the incident wave energy and wave incident angle. In normal waves and weak waves with normal wave incidence, tourists should try to avoid approaching the interval area of ​​the multi-section submerged dike; and in storm waves and oblique normal waves, they should try to avoid playing near the submerged dike to avoid the risk of drowning in rip currents.

[0146] Table 3 Maximum rip current velocity, occurrence location and rip current risk level under four wave conditions

[0147]

[0148] Example 4 - Simulation of rip current system near submerged dike with different lateral configurations

[0149] In Example 3, the impact of incident wave conditions on rip current risk under the same submerged dike configuration is considered. In this example, under normal wave conditions (i.e., wave condition B), the impact of different lateral (offshore direction) configurations of submerged dikes on the rip current risk level near the submerged dike is considered. The lateral configuration includes factors such as the width of the submerged dike, the height of the dike top, and the profile shape of the submerged dike. The example process is given below in sequence.

[0150] (1) Influence of the width of the submerged dike and the height of the dike crest

[0151] The width and height of the submerged dike directly affect the breaking morphology of the incident wave, thereby affecting the distribution characteristics of the rip current near the submerged dike. The width of the submerged dike in Example 3 is W B =32m, the height of the embankment top is H B =1.28m, corresponding to a shorter wide-top submerged dike. In this embodiment, by changing δ and H in formula (12) B The value of is used to change the submerged dike width and dike crest height, and then study the influence of the submerged dike width and dike crest height on the rip current risk level. Table 4 shows the four new tests with different submerged dike widths and dike crest heights and their corresponding rip current risk levels. In tests 1 and 2, only the value of δ is changed relative to wave condition B, and the submerged dike width of test 1 is W. B = 16m, the dimensionless width relative to the incident wavelength is W B / L w = 0.33, while the width of the submerged bank in test 2 is W B =8m, dimensionless width is W B / L w =0.17; Tests 3 and 4 only changed H relative to condition B B The value of , the height of the dike top in test 3 is H B = 1.5m, then the dimensionless water depth at the top of the embankment relative to the incident wave height is h c / H i =0.83, while the bank top height of test 4 is H B = 1.8m, dimensionless water depth at the top of the levee is h c / H i = 0.33. For the convenience of comparison, the table also gives the dimensionless submerged dike width and dike top water depth and other related parameters of wave condition B in the control group.

[0152] Table 4. Risk level of rip current near the submerged dike under different lateral configurations

[0153]

[0154] Fig.10 The flow field and rip current risk area distribution near the submerged dike of the above four groups of tests are given. From the maximum rip current velocity values ​​and rip current risk levels in this figure and Table 4, it can be seen that the rip current risk decreases with the decrease of the submerged dike width and increases with the increase of the submerged dike height. For the wide-top submerged dike condition B, the maximum rip current velocity is 1.07m / s; as the dimensionless submerged dike width decreases to W B / L w =0.33, the maximum rip current velocity decreases to 0.67m / s, and the rip current risk level also decreases. At this time, the rip current risk area also changes from the rip current neck and rip current head area mainly concentrated near the groove to the rip current compensation flow area including the submerged dike top and the coastline. As the dimensionless submerged dike width continues to decrease to W B / L w =0.17, the rip current risk area near the groove completely disappears, and the rip current risk only exists at the top of the submerged dike and the nearshore compensation flow area. At this time, the maximum rip current value is only 0.4m / s, and it occurs near the shore. The influence of the dike top height on the rip current risk is mainly reflected in the influence on the maximum rip current velocity in the groove, which can be proved by the similar rip current risk areas of Tests 3 and 4 and Wave Condition B. For Test 3, as the dimensionless dam crest water depth decreases to h c / H i =0.83, the maximum velocity of rip current increases to 1.38m / s, and the risk level of rip current reaches high risk; for test 4, as the water depth at the dam crest continues to decrease to h c / H i =0.33, the maximum rip current velocity increases sharply to 1.99m / s, and the rip current risk increases further.

[0155] The above results show that the wider the submerged dike is and the higher the dike top is, the greater the risk of rip currents is. The height of the submerged dike mainly affects the maximum rip current velocity in the trench, while the width of the submerged dike also affects the distribution of rip current risk areas. For a wide-top submerged dike, the rip current risk is mainly located near the trench area, while for a narrow-top submerged dike, it is also necessary to pay attention to the rip current risk in the nearshore area opposite the trench and the top area of ​​the submerged dike;

[0156] (2) Influence of submerged dike morphology

[0157] This embodiment considers the impact of the change in the lateral configuration of the submerged dike caused by different submerged dike shapes on the risk of rip current near the submerged dike. The Gaussian submerged dike shape of wave condition B in Example 3 is still selected as the control group, and four submerged dike shapes commonly used in engineering projects, namely elliptical, triangular, trapezoidal and rectangular, are added as the test group. Fig.11 In order to overcome the influence of the width and height of the submerged dike, the four submerged dike forms are all taken to have the same width and dike top height as the Gaussian submerged dike in wave condition B, that is, W B =32m, H B =1.28m.

[0158] Fig.12 The flow field and rip current risk area under the above four types of submerged dikes are given. It can be seen from the figure that changing the submerged dike shape does not change the rip current risk level, that is, all four types of submerged dikes show medium rip current risk; in addition, the distribution of rip current risk areas has not changed significantly, and still shows an inverted Y shape. This is because the change of the submerged dike shape has little effect on the maximum rip current in the groove. For elliptical, triangular and trapezoidal submerged dikes, the maximum rip current is not much different from the result of wave condition B, which is about 1.1m / s. For rectangular submerged dikes, the maximum rip current increases slightly to 1.27m / s. In general, the submerged dike shape has little effect on the maximum rip current velocity and the distribution of rip current risk areas.

[0159] The results of this embodiment show that the profile of the submerged dike has little effect on the risk of rip current near the submerged dike, and the main factors that determine the risk level of rip current near the submerged dike in the lateral configuration of the submerged dike are the width of the submerged dike and the height of the dike top.

[0160] Example 5 - Simulation of rip current system near submerged dike with different longitudinal configurations

[0161] In Example 4, the influence of the lateral configuration of the submerged dike on the risk level of rip current was considered. This example discusses the influence of the longitudinal (alongshore) configuration of the submerged dike. While keeping the computational domain width of 640m unchanged, this longitudinal configuration is mainly manifested in the width and number of grooves between the submerged dikes. The former directly determines the distribution range of rip currents in the grooves, while the latter affects the number of groove rip currents within a fixed range. This example still uses wave condition B as the control group, and studies the influence of the longitudinal configuration of the submerged dike on the risk of rip currents by adding other test groups. The process of this example is given in sequence below.

[0162] (1) Influence of trench width between submerged dikes

[0163] In order to consider the influence of the groove width between the submerged dikes, four groups of test cases with other groove widths are added for analysis. B = 32m, the dimensionless distance relative to the incident wavelength is G B / L w =0.67, the newly added group only changes the groove width relative to the wave condition B and keeps other parameters unchanged. The groove widths of the four newly added groups of examples in this embodiment are G B =16m, 64m, 128m and 256m, and their corresponding dimensionless trench widths are G B / L w =0.33, 1.34, 2.67 and 5.33.

[0164] Fig.13 The flow field and rip current area distribution under the above four new groups of different groove widths are given. It can be seen from the figure that the change of groove width does not change the rip current risk level. The maximum rip current velocity for different groove widths is always maintained at about 1.1m / s, which is consistent with wave condition B. However, the change of groove width causes the change of rip current risk area. For small grooves (G B =16m), the rip current risk only occurs between the grooves and the short rip current head area on the seaward side of the groove, showing an I-shaped distribution. As the groove width increases to 32m-128m, the rip current risk area shows a typical inverted Y-shaped distribution (see discussion in Example 3). When the groove width increases to 256m, the rip current risk area formed by the compensating flow on both sides is no longer closed, showing an eight-shaped distribution feature, and the maximum rip current occurs on the shore side of the groove boundary. In general, as the groove width increases, the rip current risk area increases.

[0165] The above results show that increasing the width of the grooves between the submerged dikes will not increase the risk level of rip currents, but will increase the risk area of ​​rip currents. For projects with large spacing between multi-section submerged dikes, tourists should pay more attention to the risk of rip currents at the discontinuous edges of the submerged dikes.

[0166] (2) Influence of the number of trenches between submerged dikes

[0167] The change in the number of grooves will cause the length of a single section of the multi-section submerged dike to decrease, resulting in a ratio of the length of a single section to the groove width, L B / G B In wave condition B, there is only one trench within a width of 640 m along the coast. At this time, the ratio of the length of a single submerged dike to the width of the trench is L B / G B =19. This embodiment adds four sets of tests, keeping other parameters unchanged in wave condition B, and setting the number of grooves to 2, 4, 6, and 8 respectively to study the effect of the number of grooves on the risk of rip current near the multi-section submerged dike. Corresponding to these four numbers of grooves, the ratio of the length of a single-section submerged dike to the width of the groove is L B / G B =9, 4, 2.3 and 1.5.

[0168] Fig.14 The flow field and rip current risk area near the multi-section submerged dike with 2, 4, 6, and 8 grooves are given. It can be seen from the figure that when the number of rip current grooves is less than 6, that is, L B / G B ≥2.3, the rip current risk level remains at the medium level. At this time, the maximum rip current speed is consistent with the wave condition B, maintaining at about 1.1m / s. However, as the number of grooves increases, the offshore extension range of the rip current head gradually decreases. When the number of rip currents increases to 8, that is, L B / G B =1.5, the rip head completely disappears, and the rip risk level is increased to a high level. At this time, the maximum rip velocity reaches 1.42 m / s, which is close to the storm wave situation described in Example 3 (1.53 m / s). This is caused by the change in the rip cycle characteristics. B / G B For the case of L ≥ 2.3, since the length of the submerged bank is larger than the width of the groove, the coastal compensation flow is sufficient to supply the rip current, so the rip current circulation characteristics are "jet characteristics", which can be seen from the extension range of the rip current head far from the shore; while for L B / G B =1.5, the length of the submerged dike is only 1.5 times the width of the groove, and the range of the coastal compensation flow is reduced and insufficient to supply the rip current. Therefore, the onshore flow at the top of the dike and the coastal compensation flow are superimposed and flow out in the groove in the form of rip current. The onshore flow and the rip current constitute a circulation across the submerged dike, so the rip current circulation characteristic is the "circulation mechanism", which can be seen from the obvious vortex pairs on both sides of the groove. As the driving force of the rip current, the vortex pair significantly increases the flow velocity of the rip current, thereby leading to a significant increase in the maximum velocity of the rip current. The above results show that the number of grooves in the submerged dike affects the rip current risk level by affecting the rip current circulation characteristics. When the length of the submerged dike is greater than 2.3 times the width of the groove, the rip current behaves as a "jet mechanism", and the change in the number of grooves does not affect the rip current risk level; when the length of the submerged dike is less than 1.5 times the width of the groove, the rip current behaves as a "circulation mechanism", and the rip current risk level is significantly improved.

[0169] The results of this embodiment show that under normal circumstances, the rip current circulation characteristics of multi-section submerged dikes are "jet mechanisms", and the longitudinal configuration has little effect on the rip current risk level, but will change the rip current distribution characteristics and rip current risk distribution area; however, when the rip current circulation characteristics are transformed from "jet mechanisms" to "circulation mechanisms", the rip current risk level will increase significantly, and tourists need to pay more attention to the situation where the ratio of the submerged dike length to the groove width is small.

[0170] Finally, it should be noted that the above is only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred arrangement scheme, a person skilled in the art should understand that the technical solution of the present invention (such as the use of various formulas, the sequence of steps, etc.) can be modified or replaced by equivalents without departing from the spirit and scope of the technical solution of the present invention.< / v> < / v> < / v> < / v> < / v> < / v>

Claims

1. A method for considering the risk of rip currents near a multi-segment submerged dike, characterized in that: The method comprises the following steps: Step 1, calculation model settings: First, determine the multi-section submerged dike configuration and wave field data on the beach to be studied; Then, the calculation parameters of the model are set, including wave breaking parameters, bottom friction coefficient, lateral mixing coefficient, calculation space and time step, and total calculation time; Step 2: Use the GPU-accelerated simplified model of the nonlinear gentle slope equation to simulate the rip current system on a multi-segment submerged dike: Under the model settings given in step 1, the GPU-accelerated nonlinear mild slope equation simplified model is used for simulation calculation, including: S21, using a Cartesian rectangular coordinate system (x, y, z), where the x-axis is perpendicular to the shoreline direction, the y-axis is parallel to the shoreline direction, and the z-axis is vertically upward. In this coordinate system, the GPU-accelerated nonlinear gentle slope equation simplified model has the following control equation: In the formula, t is time, and the subscript t represents the derivative with respect to time. is the symbol of the gradient operator in the horizontal direction, U = (U, V) is the horizontal velocity vector at the still water level, U and V are the flow velocities in the x and y directions on the still water surface, η is the wave surface, h is the still water depth, g is the acceleration of gravity, k is the wave number, T is the wave period, c is the wave speed, the coefficient n = 1 / 2 + kh / sinh (2kh), and R is the term that takes into account wave breaking, bottom friction and lateral mixing effects; S22, the above equations (1)-(2) are applied to the simulation of rip current system near multi-section submerged dike under the condition of oblique wave incidence. The problem of wave surface and water flow value at the lateral boundary in open water area should be considered. The lateral cyclic calculation boundary condition is used to deal with this problem, that is, three layers of virtual grids are added outside the calculation domain, and the wave surface and flow velocity results on the virtual grid nodes along the coast are cyclically assigned and then solved; S23, for the above equations (1)-(2), in order to meet the timeliness of rip current risk prediction, the model calculation efficiency should be improved, and the unified device architecture CUDA and C language are used to perform GPU parallel calculation on the model to improve the efficiency of rip current risk prediction; Step 3: Use the maximum rip flow velocity to assess the rip flow risk: Using the simulation results of the water velocity distribution on the multi-stage submerged dike given in step 2, the time-averaged and depth-averaged rip current velocity is calculated, and the rip current risk level is evaluated based on the maximum rip current velocity obtained, including: S31, when calculating the time-averaged and depth-averaged rip current velocity, it is necessary to consider 1) converting the velocity at the current static water level into the depth-averaged velocity and 2) defining the rip current velocity using the depth-averaged velocity; For problem 1), the flow rate is converted using the following formula using the flow conservation principle: Where <·> represents the average value over time, u=(u,v), u and v are the average water depth velocities in the x and y directions, respectively; For question 2), since the danger of rip current includes both the velocity in the offshore direction and the velocity in the coastal direction, the sum of the velocity vectors u and v is used as the basis for judging the magnitude of the rip current velocity. In addition, in order to exclude the influence of the velocity toward the shore, only the positive part of the velocity u perpendicular to the shore, that is, the offshore part, is taken; therefore, the time-averaged and depth-averaged rip current velocity value is: In the formula, u + It means that the velocity u takes a positive value; S32, according to the full-field rip current velocity value results given by the above formula (11), select the maximum rip current velocity value |u| max ,The rip current risk levels are divided to quantitatively assess the risk of rip currents near multi-segment submerged dikes: |u| max <0.3m / s is considered as no risk of rip flow, |u| max Offshore currents ≥0.3m / s are considered rip current hazards. Specifically, 0.3m / s≤|u| max When <0.8m / s, it is considered as low risk of rip current, 0.8m / s≤|u| max When <1.3m / s, it is considered as medium risk of rip current, |u| max When the flow rate is ≥1.3m / s, it is considered as a high risk of rip current.

2. The method according to claim 1, characterized in that In step 1, the multi-section submerged breakwater configuration and wave field data on the beach to be studied include the scope of the beach to be studied, the still water depth of the seabed topography, the profile type, width, length in the coastal direction and the size of the submerged breakwater spacing in the offshore direction, and the wave height, period and incident angle of the incident wave.

3. The method according to claim 1, characterized in that The specific solution process of S22 in step 2 is: Consider that the total number of grids in the computational domain is M in the x direction and N in the y direction, where the value range of the grid nodes (i, j) in the computational domain is i = 1 to M and j = 1 to N, and the value of the virtual grid node (i, j) is i = 0, -1, -2 at the shore boundary of the computational domain; i = M + 1, M + 2, M + 3 at the offshore boundary; j = 0, -1, -2 at the left boundary; and j = N + 1, N + 2, N + 3 at the right boundary; Therefore, the values ​​of wave surface η, flow velocity U and V on the virtual grid nodes on the left boundary are as follows: η(i,1-j)=η(i,N+1-j), j=1,2,3 (3) U(i,1-j)=U(i,N+1-j), j=1,2,3 (4) V(i,1-j)=V(i,N+1-j), j=1,2,3 (5) On the right boundary, the values ​​of the wave surface η, flow velocities U and V in the virtual grid are as follows: η(i,N+j)=η(i,j), j=1,2,3 (6) U(i,N+j)=U(i,j), j=1,2,3 (7) V(i,N+j)=V(i,j), j=1,2,3 (8) When the incident angle α of the offshore wave is known and the number of computational grids N in the y direction along the coast is fixed, in order to meet the lateral circulation boundary conditions, the value of the computational space step length △y in the coastal direction needs to be adjusted so that the coastal direction contains an integer number of oblique incident waves. Therefore, the value of △y should satisfy the following relationship: Where, L w =2π / k is the deep water wavelength, M w is an adjustable positive integer parameter; After obtaining △y according to the above formula, readjust the calculation domain width to (N-1)△y.

4. The method according to claim 1, characterized in that: In step (2), S23, the process of GPU parallel computing of the model based on CUDA and C language includes: declaration and definition of host and device variables, allocation of host and device variable memory space, host variable initialization, copying host data to device, calling device kernel functions in sequence to complete the calculation, copying device data back to host, and finally outputting the result; In the process of calling the kernel function, the number of logically enabled threads is organized at two levels. The first level is to specify the grid size, that is, the number of thread blocks contained in the grid in the x and y directions; the second level is to define the size of the thread block, that is, the number of threads contained in a thread block.

Images