A parametric calculation method for wave dissipation by flexible vegetation

Abstract

The invention discloses a parameterized calculation method for flexible vegetation wave elimination, which relates to the technical field of numerical calculation methods for flexible vegetation wave elimination, and comprises the following steps: S2: determining a set of grids in a vegetation area according to the total height of vegetation; S3: generating a set of grid surfaces; S4: obtaining a volume scalar matrix; S5: reading parameters in a dictionary file; S6: constructing interpolation sampling lines according to each element in the volume scalar matrix of S4; S7: calculating the instantaneous vertical average water flow velocity and instantaneous wave surface displacement on each interpolation sampling line, and calculating the equivalent height of vegetation at each interpolation sampling line; S8: obtaining a volume scalar matrix according to the equivalent height of vegetation at each interpolation sampling line; S9: calculating the drag coefficient through a grid whose volume scalar median value in S8 is 1. The method of the invention can calculate the instantaneous equivalent height and the drag coefficient of vegetation during the calculation process.

Description

Technical Field

[0001] The invention relates to the technical field of numerical calculation methods for flexible vegetation wave dissipation, and in particular to a parameterized calculation method for flexible vegetation wave dissipation. Background Art

[0002] The numerical calculation methods of flexible vegetation wave elimination based on computational fluid dynamics (CFD) model can be mainly divided into direct calculation method and parameterized calculation method. In recent years, some scholars have realized the fluid-solid coupling calculation of flexible vegetation and fluid interaction by combining computational fluid dynamics method and finite element method. The above direct calculation method needs to describe the detailed characteristics of vegetation or calculate the fluid-solid coupling process, resulting in excessive calculation amount, making it difficult to realize the calculation of high-density vegetation wave elimination and unable to meet engineering needs. The parameterized calculation method generalizes the vegetation area into a resistance area, and realizes the calculation of vegetation wave elimination by incorporating the resistance source term of vegetation into the control equation. In this method, there is no need to describe each plant, so it has high calculation efficiency. The control equation of the parameterized calculation method is mainly the NS equation (Navier-Stokes equations), and the VOF equation (Volume of Fluid equation) is combined to calculate the fluctuation of the free liquid surface.

[0003] The current parameterized calculation method in CFD calculation generally first delineates the vegetation area (specify which grids belong to the vegetation area before calculation), and applies resistance and inertia terms to the grids in the vegetation area during the calculation process. The range of the vegetation area is immutable during the calculation process (as shown in Figure 1, parameterized vegetation area (red area) and wave surface line (black line)), so the vegetation wave dissipation parameterized model in CFD is mainly used to calculate the wave dissipation of rigid vegetation. In addition, the drag coefficient in the entire vegetation area Also remains unchanged, this value is mainly determined by some empirical formulas or by calibration with experimental data.

[0004] Some scholars have also applied parametric calculation methods to the calculation of flexible vegetation by introducing the concept of equivalent height. Specifically, the phenomenon of lowering the height of vegetation caused by the bending of flexible vegetation under the action of waves is considered, and the bent flexible vegetation is regarded as rigid vegetation at a certain height (this height is usually less than the height of the flexible vegetation itself under still water conditions), so that the flexible vegetation and the rigid vegetation at this height have the same wave-breaking ability, and this height becomes the equivalent height.

[0005] However, the above-mentioned parameterized calculation method still has the following problems when calculating flexible vegetation: 1. In the current equivalent height method, the equivalent height calculated in advance is used to determine the vegetation area in the wave-breaking calculation process, and the determined flexible vegetation area remains unchanged during the calculation process. This point is inconsistent with the actual situation. In practice, flexible vegetation will continue to fluctuate with the movement of waves. In addition, the formula for calculating the equivalent height requires the maximum horizontal velocity of the wave and the maximum water particle displacement of the wave. The existing methods mostly select the wave height and period of the leading edge of the vegetation area to determine the wave velocity and water particle displacement. The equivalent height calculated by this data is applied to the entire vegetation area during calculation. However, in actual situations, if the vegetation area is long, the waves are constantly attenuating in the vegetation area, so the equivalent height of the vegetation should also be constantly changing. The existing model cannot accurately reflect the changes in the instantaneous equivalent height of the flexible vegetation. 2. The determination of the drag coefficient of the vegetation area also has similar problems as the equivalent height, that is, The value of remains unchanged in the vegetation area, but according to the calculation formula of the resistance source term The value of is related to the instantaneous force on the vegetation and the instantaneous local flow velocity (e.g. Figure 2 As shown in the literature, from Hu Z, Suzuki T, Zitman T, et al. Laboratory study onwave dissipation by vegetation in combined current–wave flow[J]. Coastal Engineering, 2014, 88: 131-142.), constant It is difficult to accurately reflect the change characteristics of the resistance encountered by waves in the vegetation area, so it is difficult to accurately calculate the change characteristics of the flow field in the vegetation area (especially when there is sediment transport in the vegetation area); in addition, in order to obtain better calculation results of wave height change in most models, The values ​​of are calibrated by comparing with experimental data, but for calculations lacking experimental data (such as calculations on prototype scales), It is difficult to determine.

[0006] In view of this, this invention is proposed. Summary of the invention

[0007] The purpose of the present invention is to solve the shortcomings existing in the prior art, and a parameterized calculation method for flexible vegetation wave breaking is proposed. During the calculation process, the instantaneous equivalent height and drag force coefficient of the flexible vegetation can be calculated according to the instantaneous vertical average flow velocity and wave surface displacement. It has been verified that the method of the present invention can better predict the wave breaking process in the flexible vegetation area.

[0008] In order to achieve the above object, the present invention adopts the following technical solutions:

[0009] A parameterized calculation method for flexible vegetation wave elimination, characterized in that it comprises the following steps:

[0010] S1: Use OpenFOAM to establish the computational area grid;

[0011] S2: Total height of vegetation under still water conditions , determine the set cellZone of the grids in the vegetation area through the setFields and topoSet tools;

[0012] S3: Generate mesh face set faceZone through topoSet and boxToFace tools;

[0013] S4: traverse the grid center points after cellZone projection through the pre-processing program cellFaceProject, obtain the volume scalar matrix cellFaceAddressing, and write cellFaceAddressing into the initial moment folder;

[0014] S5: read the parameters in the dictionary file, including: vegetation stem width, vegetation stem height, density of vegetation stems per unit area, vegetation leaf width, vegetation leaf thickness, vegetation leaf length, number of leaves on each stem, Young's modulus of vegetation stems and leaves, still water depth, sampling line length, sampling interpolation method, coefficients A, B and C in the drag force calculation formula;

[0015] S6: construct an interpolation sampling line according to the value of each element in cellFaceAddressing;

[0016] S7: Calculate the instantaneous vertical average water flow velocity and instantaneous wave surface displacement on each interpolation sampling line, and calculate the equivalent height of vegetation at each interpolation sampling line ;

[0017] S8: Based on the equivalent height of vegetation at each interpolation sampling line Get the volume scalar matrix vegetationCell;

[0018] S9: traverse cellZone, calculate the drag force coefficient through the grid with a value of 1 in vegetationCell, and solve the NS equation and VOF equation based on the calculated drag force coefficient.

[0019] Further, the S4 comprises the following steps:

[0020] S41: Project the vertex coordinates of each face in faceZone to the horizontal plane of z=0 to obtain a number of projected faces;

[0021] S42: Project the coordinates of each grid center point in cellZone to the horizontal plane of z=0 to obtain a number of projected grid center points;

[0022] S43: construct a volume scalar matrix cellFaceAddressing, wherein each element corresponds one-to-one to the center point of the grid after the projection of cellZone in S42, and the initial value of each element is -1;

[0023] S44: traverse the center point of the grid after cellZone projection in S42, and determine whether it is located on any surface after faceZone projection in S41 by using the ray method. If so, assign the element corresponding to the center point of the grid after projection in cellFaceAddressing to the number value of the surface after faceZone projection corresponding to it;

[0024] S45: After traversing the center points of each grid after the projection of cellZone in S42, the value on the cellFaceAddressing boundary is updated through the zero gradient boundary condition to obtain the final cellFaceAddressing, and write it into the initial moment folder.

[0025] Further, the S6 comprises the following steps:

[0026] S61: Use the std::map function in the C++ standard library to establish a mapping between the elements in cellFaceAddressing and the coordinates of the face center of the bottom boundary grid of the computational domain;

[0027] S62: Calculate the coordinate average of the center coordinates of the bottom boundary grids with the same mapping value, and use the coordinate average as the bottom coordinate of the interpolation sampling line;

[0028] S63: translating vertically upward by a distance within a preset range as the top coordinate of the interpolation sampling line;

[0029] S64: Construct the interpolation sampling lines of each mesh cluster through the midPointAndFaceSet class.

[0030] Furthermore, before S61, the following steps are also included:

[0031] S60: Read the value of each element in cellFaceAddressing and calculate the coordinates of all mesh face centers on the bottom boundary of the regional mesh.

[0032] Furthermore, before S61, the following steps are also included:

[0033] S601: Obtain the value of each element in cellFaceAddressing on each processor and the coordinates of all mesh face centers on the bottom boundary of the calculation area mesh through the Pstream::gatherList function;

[0034] In the S61, a mapping between the elements in cellFaceAddressing and the mesh face center coordinates is established on the main processor through the std::map function in the C++ standard library;

[0035] In the step S62, the set of bottom coordinates of the interpolation sampling line in the main processor is sent to each of the other processors;

[0036] After S64, the following steps are also included:

[0037] S65: The gatherOnly function in the globalIndex class is used to obtain the sampling point coordinate set of the interpolation sampling line stored in the other processors on the main processor, and the gatherSort function in the coordSet is used to establish a sorting map of the top coordinates and the bottom coordinates of each sampling point on the interpolation sampling line on the main processor.

[0038] Further, the S7 comprises the following steps:

[0039] S71: at each time step, interpolating the flow velocity u and the liquid volume fraction α on each interpolation sampling line by using a set sampling interpolation method;

[0040] S72: Integrate u and α by trapezoidal integration method to obtain the instantaneous vertical average flow velocity u(t) and instantaneous water depth on each interpolation sampling line, and take the absolute value of the difference between the instantaneous water depth and the static water depth as the instantaneous wave surface displacement A(t);

[0041] S72: Calculate the equivalent height of vegetation at each interpolation sampling line based on the instantaneous vertical average flow velocity and instantaneous wave surface displacement on each interpolation sampling line .

[0042] Further, the S72 calculates the equivalent height of the vegetation through the following steps: :

[0043] S721: Calculate the Cauchy number of vegetation stems and the Cauchy number of vegetation leaves , the calculation formula is as follows:

[0044] ,

[0045] ,

[0046] ,

[0047] ,

[0048] In the formula, is the density of water, is the width of the vegetation stem, is the width of vegetation leaves, is the height of the vegetation stem, is the height of vegetation leaves, is the Young's modulus of the vegetation stem, is the Young's modulus of vegetation leaves, is the thickness of vegetation leaves;

[0049] S722: Calculate the equivalent height of vegetation stems and the equivalent height of vegetation leaves , the calculation formula is as follows:

[0050] ,

[0051] ,

[0052] ,

[0053] ,

[0054] In the formula, K is a constant, which can generally be taken as 0.94;

[0055] S373: Calculate the equivalent height of vegetation , the calculation formula is as follows:

[0056] ,

[0057] In the formula, is the number of leaves on a single stem.

[0058] Further, the S71 includes the following steps:

[0059] S711: at each time step, interpolating the flow velocity u and the liquid volume fraction α on each interpolation sampling line by using a set sampling interpolation method;

[0060] S712: on the main processor, obtain the interpolation results on the other processors and the coordinates of the interpolation points on the corresponding sampling lines through the gatherInplace function in the globalIndex class;

[0061] S713: on the main processor, the interpolation results and the corresponding interpolation point coordinates obtained on the remaining processors are sorted by using the sorting mapping in the gatherSort function in the coordSet class;

[0062] S72 includes the following steps:

[0063] S721: Integrate u and α on the main processor by using the trapezoidal integration method to obtain the instantaneous vertical average flow velocity u(t) and the instantaneous water depth on each interpolation sampling line, and take the absolute value of the difference between the instantaneous water depth and the static water depth as the instantaneous wave surface displacement A(t);

[0064] S722: Send u(t) and A(t) calculated on the main processor to the remaining processors.

[0065] Further, the S8 comprises the following steps:

[0066] S81: construct a body scalar matrix vegetationCell, in which each element corresponds one-to-one to each element in cellFaceAddressing, and the initial value of each element is 0;

[0067] S82: Determine whether the height of the center of the grid in the grid cluster corresponding to each element in cellFaceAddressing is less than the equivalent height of the vegetation at the corresponding interpolation sampling line If so, set the corresponding element of vegetationCell to 1.

[0068] Further, in S9, the drag coefficient is calculated by the following formula: :

[0069]

[0070] ,

[0071] In the formula, A, B and C are constants, and μ is the dynamic viscosity of water.

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

[0073] 1. During the calculation process, the instantaneous equivalent height of the flexible vegetation is calculated in cellZone according to the instantaneous vertical average flow velocity and wave surface displacement. The flexible vegetation area that plays a wave-breaking role in the calculation area grid can be determined according to the instantaneous equivalent height. In addition, the drag coefficient is calculated in cellZone according to the instantaneous vertical average flow velocity, which can realize high-efficiency and high-precision calculation of flexible vegetation.

[0074] 2. It has been verified that the calculated effective wave height and drag coefficient are consistent with the actual vegetation movement situation, can more accurately calculate the trend of wave height attenuation along the vegetation area, and can better predict the wave dissipation process in the flexible vegetation area. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 It is an example diagram of vegetation area of ​​parameterized calculation method in the prior art;

[0076] Figure 2 Variation diagram of (a) instantaneous local flow velocity, (b) instantaneous resistance and (c) calculated instantaneous drag coefficient for measured vegetation;

[0077] Figure 3 It is a flow chart of the parameterized calculation method of flexible vegetation wave elimination in Example 1;

[0078] Figure 4 It is a structured grid distribution map of vegetation area;

[0079] Figure 5 It is the unstructured grid distribution map of vegetation area;

[0080] Figure 6 It is a faceZone diagram in the vegetation area;

[0081] Figure 7 It is the distribution map of cellFaceAddressing values ​​at different locations in the vegetation area;

[0082] Figure 8 Schematic diagram of the flume experiment arrangement of Liu et al. (2023);

[0083] Fig. 9 It is the change process of the flexible vegetation area determined by the calculated effective height within a wave cycle;

[0084] Fig.10 for Fig. 9 The change process of the drag force coefficient at the corresponding moment;

[0085] Fig.11 This is a comparison chart of the changes in measured wave height and calculated wave height near the vegetation area under different working conditions. DETAILED DESCRIPTION

[0086] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

[0087] Embodiment 1:

[0088] A parameterized calculation method for flexible vegetation wave breaking, such as Figure 3 As shown, the following steps are included:

[0089] S1: Use OpenFOAM to establish the computational area grid.

[0090] In this embodiment, the calculation area grid is the calculated water tank area. In the following, x is the length direction of the area, y is the width direction of the area, and z is the vertical direction of the area.

[0091] S2: Total height of vegetation under static water conditions , determine the collection cellZone of grids in the vegetation area through the setFields and topoSet tools.

[0092] In this embodiment, cellZone is a grid set of all grids in the vegetation area.

[0093] S3: Generate mesh face set faceZone through topoSet and boxToFace tools.

[0094] In this embodiment, when generating faceZone, it is ensured that the parallelism between faceZone and the bottom boundary of the calculation area grid is less than a preset threshold, and each face in faceZone is larger in vertical area than any face parallel to it in cellZone, that is, each face in faceZone is the face with the largest area in the vertical direction, and the scale of faceZone in the x and y directions is the same as that of cellZone.

[0095] S4: Use the self-developed pre-processing program cellFaceProject to traverse the grid center points after cellZone projection, obtain the volume scalar matrix cellFaceAddressing, and write cellFaceAddressing into the initial moment folder.

[0096] S5: Read the parameters in the dictionary file, including: vegetation stem width, vegetation stem height, density of vegetation stems per unit area, vegetation leaf width, vegetation leaf thickness, vegetation leaf length, number of leaves on each stem, Young's modulus of vegetation stems and leaves, still water depth, sampling line length, sampling interpolation method, and coefficients A, B, and C in the drag force calculation formula.

[0097] In this embodiment, the parameters in the dictionary file are read, and the parameters that must be read include: vegetation stem width, vegetation stem height, density of vegetation stems per unit area, Young's modulus of vegetation stems, still water depth, sampling line length, sampling interpolation method, and coefficients A, B, and C in the drag force calculation formula; wherein the sampling interpolation method can adopt the interpolation method provided in OpenFOAM, such as the "cellPoint" method; in addition, the parameters that can be read optionally or conditionally include: vegetation leaf width, vegetation leaf thickness, vegetation leaf length, Young's modulus of vegetation leaves, and the number of leaves on each stem.

[0098] S6: Construct an interpolation sampling line according to the value of each element in cellFaceAddressing.

[0099] S7: Calculate the instantaneous vertical average water flow velocity and instantaneous wave surface displacement on each interpolation sampling line, and calculate the equivalent height of vegetation at each interpolation sampling line .

[0100] S8: Based on the equivalent height of vegetation at each interpolation sampling line Get the volume scalar matrix vegetationCell.

[0101] S9: traverse cellZone, calculate the drag force coefficient through the grid with a value of 1 in vegetationCell, and solve the NS equation and VOF equation based on the calculated drag force coefficient.

[0102] In this embodiment, the VOF equation and the NS equation may be solved by using the existing method, which will not be elaborated here.

[0103] In this embodiment, S7-S9 are executed cyclically at each time step until the calculation is completed.

[0104] The parameterized calculation method for flexible vegetation wave breaking in this embodiment, during the calculation process, calculates the instantaneous equivalent height of the flexible vegetation in cellZone according to the instantaneous vertical average flow velocity and the wave surface displacement, and the flexible vegetation area that plays a wave breaking role in the calculation area grid can be determined according to the instantaneous equivalent height; in addition, the drag force coefficient is calculated in cellZone according to the instantaneous vertical average flow velocity, which can achieve high-efficiency and high-precision flexible vegetation calculation.

[0105] In an optional embodiment, S4 includes the following steps: S41: projecting the vertex coordinates of each face in faceZone to the horizontal plane of z=0 to obtain a number of projected faces; S42: projecting the coordinates of each grid center point in cellZone to the horizontal plane of z=0 to obtain a number of projected grid center points; S43: constructing a volume scalar matrix cellFaceAddressing, wherein each element corresponds to a grid center point after cellZone projection in S42, wherein the initial value of each element is -1; S44: traversing the grid center points after cellZone projection in S42, and judging whether it is located on any face after faceZone projection in S41 by ray method, if so, assigning the element corresponding to the grid center point after projection in cellFaceAddressing to the number value of the faceZone face corresponding to it. S45: after traversing each grid center point after cellZone projection in S42, updating the value on the cellFaceAddressing boundary by zero gradient boundary condition, obtaining the final cellFaceAddressing, and writing it to the initial moment folder.

[0106] In this optional embodiment, in S44, if a certain projected grid center point in cellZone is not located on any projected surface of faceZone, the element corresponding to the projected grid center point in cellFaceAddressing is not assigned a value, that is, the element corresponding to the projected grid center point in cellFaceAddressing is an initial value of -1; the number value of the faceZone surface is the faceId value in OpenFOAM.

[0107] In the flexible vegetation calculation process, the grids at different x, y positions and on the same vertical line are clustered. If the grids in the vegetation area are structured grids and are evenly distributed in the vertical direction (as shown in Figure 4), the grids in the cellZone can be traversed according to the coordinates of the center of the bottom boundary grid of the vegetation area, and the grids with the same x, y coordinates can be clustered. In actual calculations, the amount of calculation using structured grids is usually large and is generally affected by factors such as complex terrain. Therefore, unstructured grids are often used in the calculation area (as shown in Figure 5, Figure 4 and 5 In the figure, the parameterized vegetation area is the red area, the static water level is the black line area, and the computational grid is the yellow line area). However, the center of the mesh surface in the bottom boundary encryption area cannot correspond one-to-one with the center of the mesh body inside the computational area grid. Therefore, the method of traversing and comparing coordinates is not applicable to the calculation of unstructured grids.

[0108] Based on this, S3 and S41-S45 are designed in this embodiment and the optional embodiment, which are applicable to the vegetation area clustering method of structured and unstructured grids. The faceZone schematic diagram in the vegetation area is as follows: Figure 6 As shown in the figure, faceZone is the white line in the figure, and cellFaceAddressing at different locations in the vegetation area is as follows Figure 7 As shown, the grid lines in the figure are black lines.

[0109] In an optional embodiment, S6 includes the following steps: S61: establishing a mapping between the elements in cellFaceAddressing and the coordinates of the center of the bottom boundary grid of the calculation domain through the std::map function in the C++ standard library; S62: calculating the coordinate average of the center of the bottom boundary grid with the same mapping value, and using the coordinate average as the bottom coordinate of the interpolation sampling line; S63: translating vertically upward by a distance within a preset range as the top coordinate of the interpolation sampling line; S64: constructing the interpolation sampling lines of each grid cluster through the midPointAndFaceSet class.

[0110] In this optional embodiment, in S62, the mapping value of the grid face center coordinate is the value of the element in cellFaceAddressing mapped thereto; in S63, the translation distance is represented by sampLength, and its minimum length needs to be greater than the water depth + 0.5 times the wave height.

[0111] In an optional embodiment, the following steps are further included before S61:

[0112] S60: Read the value of each element in cellFaceAddressing and calculate the coordinates of all mesh face centers on the bottom boundary of the regional mesh.

[0113] In this optional embodiment, it is used in a single-core computing scenario.

[0114] In an optional embodiment, the following steps are also included before S61: S601: Obtain the value of each element in cellFaceAddressing on each processor and the coordinates of all mesh face centers on the bottom boundary of the calculation area mesh through the Pstream::gatherList function; in S61, establish a mapping between the elements in cellFaceAddressing and the mesh face center coordinates on the main processor through the std::map function in the c++ standard library; in S62, send the set of bottom coordinates of the interpolation sampling lines in the main processor to the remaining processors, and the remaining processors refer to processors other than the main processor. The following steps are also included after S64: S65: Through the gatherOnly function in the globalIndex class, the main processor obtains the set of sampling point coordinates of the interpolation sampling lines stored in the remaining processors, and uses the gatherSort function in coordSet to establish a sorted mapping of the coordinates of each sampling point on the interpolation sampling line on the main processor.

[0115] In this optional embodiment, a parallel computing scenario is used.

[0116] In an optional embodiment, the S7 includes the following steps: S71: at each time step, the flow velocity u and the liquid volume fraction α are interpolated on each interpolation sampling line by a set sampling interpolation method; S72: u and α are integrated by a trapezoidal integration method to obtain the instantaneous vertical average flow velocity u(t) and the instantaneous water depth on each interpolation sampling line, respectively, and the absolute value of the difference between the instantaneous water depth and the static water depth is taken as the instantaneous wave surface displacement A(t); S72: the equivalent height of the vegetation at each interpolation sampling line is calculated according to the instantaneous vertical average flow velocity and the instantaneous wave surface displacement on each interpolation sampling line. .

[0117] In this optional embodiment, it is used in a single-core computing scenario.

[0118] In an optional embodiment, the S72 calculates the equivalent height of the vegetation by the following steps: :

[0119] S721: Calculate the Cauchy number of vegetation stems and the Cauchy number of vegetation leaves , the calculation formula is as follows:

[0120] ,

[0121] ,

[0122] ,

[0123] ,

[0124] In the formula, is the density of water, is the width of the vegetation stem, is the width of vegetation leaves, is the height of the vegetation stem, is the height of vegetation leaves, is the Young's modulus of the vegetation stem, is the Young's modulus of vegetation leaves, is the thickness of vegetation leaves, where , , , , Read directly during calculation.

[0125] S722: Calculate the equivalent height of vegetation stems and the equivalent height of vegetation leaves , the calculation formula is as follows:

[0126] ,

[0127] ,

[0128] ,

[0129] ,

[0130] In the formula, K is a constant, which can generally be taken as 0.94.

[0131] S373: Calculate the equivalent height of vegetation , the calculation formula is as follows:

[0132] ,

[0133] In the formula, is the number of leaves on a single stem, read directly before calculation.

[0134] In an optional embodiment, the S71 includes the following steps: S711: at each time step, the flow velocity u and the liquid volume fraction α are interpolated on each interpolation sampling line by the set sampling interpolation method; S712: on the main processor, the interpolation results on each other processor and the coordinates of the interpolation points on the corresponding sampling lines are obtained by the gatherInplace function in the globalIndex class; S713: on the main processor, the interpolation results on each other processor and the corresponding interpolation point coordinates are re-sorted by the sorting mapping in the gatherSort function in the coordSet class; S72 includes the following steps: S721: on the main processor, u and α are integrated by the trapezoidal integration method to obtain the instantaneous vertical average flow velocity u(t) and instantaneous water depth on each interpolation sampling line respectively, and the absolute value of the difference between the instantaneous water depth and the static water depth is used as the instantaneous wave surface displacement A(t); S722: u(t) and A(t) calculated on the main processor are sent to each other processor.

[0135] In this optional embodiment, a parallel computing scenario is used.

[0136] In an optional embodiment, the S8 includes the following steps: S81: constructing a scalar matrix vegetationCell, wherein each element corresponds to each element in cellFaceAddressing one by one, wherein the initial value of each element is 0; S82: determining whether the height of the center of the grid in the grid cluster corresponding to each element in cellFaceAddressing is less than the equivalent height of the vegetation at the corresponding interpolation sampling line If so, set the corresponding element of vegetationCell to 1.

[0137] In an optional embodiment, the drag force coefficient is calculated in S9 by the following formula: :

[0138] ,

[0139] ,

[0140] In the formula, A, B and C are constants, μ is the dynamic viscosity of water, where A can be 0.57, B can be 1546, and C can be 1.11.

[0141] The parameterized calculation method of flexible vegetation wave breaking in this embodiment is verified with reference to the test results of the flexible vegetation wave breaking flume conducted by Liu et al. (2023).

[0142] The actual test was conducted in a wave flume 80m long, 0.5m wide and 1.2m deep. The total length of the vegetation area was 3.16m, located in the middle of the flume. Three types of cylindrical vegetation were installed in the vegetation area in a staggered arrangement across the width of the flume. The test used a flexible material, polytetrafluoroethylene (PTFE) foam tube, whose Young's modulus was measured to be 1.34Gpa. Under static water conditions, the height and diameter of the vegetation were and The stem density in the vegetation area is Plant / m 2 In the experiment, the wave height at 7 positions was measured using a resistive wave height meter. The measurement positions were x = 0m, 0.6m, 1.2m, 1.8m, 2.4m, 3m and 3.6m, respectively. x = 0m is at the front of the vegetation area. The layout and experimental settings of the wave flume are as follows: Figure 8 shown.

[0143] Four working conditions under different water depths and period conditions in Liu et al. (2023) were selected to verify the parameterized calculation method of flexible vegetation wave elimination in this embodiment. The specific parameter settings are shown in Table 1.

[0144] Table 1

[0145] Working conditions Water depth d / m <![CDATA[Incident wave height H 0 / m]]> Wave period T / s one 0.3 0.08 1.4 two 0.35 0.08 1.4 three 0.4 0.08 1.6 Four 0.5 0.08 1.6

[0146] Based on the above experimental settings, a numerical wave flume was established in OpenFOAM. The length (x) of the flume was 19m, the height (z) was 0.8m, the width (y) was 0.1m, and the vegetation area was located 5m away from the wave-making boundary. The grid length in the x and z directions of the wave flume was 0.02m, and the y direction was set to 1 grid. The grid was encrypted within the wave height and vegetation range, and the grid scale in the x and z directions was 0.005m. The calculated Courant number was controlled below 0.3.

[0147] According to the parameterized calculation method of flexible vegetation wave elimination in this embodiment, Fig. 9 is the change process of the flexible vegetation area determined by the calculated effective height within one wave cycle under condition 4. Fig.10 is the change process of the drag coefficient at the corresponding moment. It can be found that at the crest and trough positions, the flexible vegetation area is lower due to the larger average flow velocity of the cross section and the larger wave surface displacement, while the flexible vegetation area is higher at the wave equilibrium position. This law is consistent with the formula calculation results and the actual vegetation movement. In addition, the change law of the drag coefficient also conforms to Figure 2The change trend of the drag coefficient in the wave surface is shown in Figure 11, and the drag coefficient is larger at the wave surface equilibrium position, while the drag coefficient is relatively smaller at the wave crest and trough positions. Figure 11 is a comparison diagram of the measured wave height and the calculated wave height changes near the vegetation area under different working conditions. It can be found that the parameterized calculation method of the flexible vegetation wave elimination in this embodiment can more accurately calculate the change trend of the wave height attenuation along the vegetation area, and can better predict the wave elimination process in the flexible vegetation area.

[0148] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical scheme and inventive concept of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. A parameterized calculation method for flexible vegetation wave elimination, characterized in that: The steps include: S1: Use OpenFOAM to establish the computational area grid; S2: Total height of vegetation under still water conditions , determine the set cellZone of the grids in the vegetation area through the setFields and topoSet tools; S3: Generate mesh face set faceZone through topoSet and boxToFace tools; S4: traverse the grid center points after cellZone projection through the pre-processing program cellFaceProject, obtain the volume scalar matrix cellFaceAddressing, and write cellFaceAddressing into the initial moment folder; S5: read the parameters in the dictionary file, including: vegetation stem width, vegetation stem height, density of vegetation stems per unit area, vegetation leaf width, vegetation leaf thickness, vegetation leaf length, number of leaves on each stem, Young's modulus of vegetation stems and leaves, still water depth, sampling line length, sampling interpolation method, coefficients A, B and C in the drag force calculation formula; S6: construct an interpolation sampling line according to the value of each element in cellFaceAddressing; S7: Calculate the instantaneous vertical average water flow velocity and instantaneous wave surface displacement on each interpolation sampling line, and calculate the equivalent height of vegetation at each interpolation sampling line ; S8: Based on the equivalent height of vegetation at each interpolation sampling line Get the volume scalar matrix vegetationCell; S9: traverse cellZone, calculate the drag force coefficient through the grid with a value of 1 in vegetationCell, and solve the NS equation and VOF equation based on the calculated drag force coefficient.

2. A parameterized calculation method for flexible vegetation wave elimination according to claim 1, characterized in that: The S4 comprises the following steps: S41: Project the vertex coordinates of each face in faceZone to the horizontal plane of z=0 to obtain a number of projected faces; S42: Project the coordinates of each grid center point in cellZone to the horizontal plane of z=0 to obtain a number of projected grid center points; S43: construct a volume scalar matrix cellFaceAddressing, where each element corresponds one-to-one to the center point of the grid after the projection of cellZone in S42, and the initial value of each element is -1; S44: traverse the center point of the grid after the projection of cellZone in S42, and determine whether it is located on any surface after the projection of faceZone in S41 by using the ray method. If so, assign the element corresponding to the center point of the grid after the projection in cellFaceAddressing to the number value of the corresponding faceZone surface; S45: After traversing the center points of each grid after the projection of cellZone in S42, the value on the cellFaceAddressing boundary is updated through the zero gradient boundary condition to obtain the final cellFaceAddressing, and write it into the initial moment folder.

3. A parameterized calculation method for flexible vegetation wave elimination according to claim 2, characterized in that: The S6 comprises the following steps: S61: Use the std::map function in the C++ standard library to establish a mapping between the elements in cellFaceAddressing and the coordinates of the face center of the bottom boundary grid of the computational domain; S62: Calculate the coordinate average of the center coordinates of the bottom boundary grids with the same mapping value, and use the coordinate average as the bottom coordinate of the interpolation sampling line; S63: translating vertically upward by a distance within a preset range as the top coordinate of the interpolation sampling line; S64: Construct the interpolation sampling lines of each mesh cluster through the midPointAndFaceSet class.

4. A parameterized calculation method for flexible vegetation wave elimination according to claim 3, characterized in that: Before S61, the following steps are also included: S60: Read the value of each element in cellFaceAddressing and calculate the coordinates of all mesh face centers on the bottom boundary of the regional mesh.

5. A parameterized calculation method for flexible vegetation wave elimination according to claim 3, characterized in that: Before S61, the following steps are also included: S601: Obtain the value of each element in cellFaceAddressing on each processor and the coordinates of all mesh face centers on the bottom boundary of the calculation area mesh through the Pstream::gatherList function; In the S61, a mapping between the elements in cellFaceAddressing and the mesh face center coordinates is established on the main processor through the std::map function in the C++ standard library; In the step S62, the set of bottom coordinates of the interpolation sampling line in the main processor is sent to each of the other processors; After S64, the following steps are also included: S65: The gatherOnly function in the globalIndex class is used to enable the main processor to obtain the sampling point coordinate set of the interpolation sampling line stored in each of the other processors, and the gatherSort function in the coordSet is used to establish a sorting map of the coordinates of each sampling point on the interpolation sampling line on the main processor.

6. A parameterized calculation method for flexible vegetation wave elimination according to claim 3, characterized in that: The S7 comprises the following steps: S71: at each time step, interpolating the flow velocity u and the liquid volume fraction α on each interpolation sampling line by using a set sampling interpolation method; S72: Integrate u and α by trapezoidal integration method to obtain the instantaneous vertical average flow velocity u(t) and instantaneous water depth on each interpolation sampling line, and take the absolute value of the difference between the instantaneous water depth and the static water depth as the instantaneous wave surface displacement A(t); S72: Calculate the equivalent height of vegetation at each interpolation sampling line based on the instantaneous vertical average flow velocity and instantaneous wave surface displacement on each interpolation sampling line .

7. A parameterized calculation method for flexible vegetation wave elimination according to claim 6, characterized in that: The S72 calculates the equivalent height of vegetation by the following steps: : S721: Calculate the Cauchy number of vegetation stems and the Cauchy number of vegetation leaves , the calculation formula is as follows: ， ， ， ， In the formula, is the density of water, is the width of the vegetation stem, is the width of vegetation leaves, is the height of the vegetation stem, is the height of vegetation leaves, is the Young's modulus of the vegetation stem, is the Young's modulus of vegetation leaves, is the thickness of vegetation leaves; S722: Calculate the equivalent height of vegetation stems and the equivalent height of vegetation leaves , the calculation formula is as follows: ； ； ； ； In the formula, K is a constant, which can generally be taken as 0.94; S373: Calculate the equivalent height of vegetation , the calculation formula is as follows: ； In the formula, is the number of leaves on a single stem.

8. A parameterized calculation method for flexible vegetation wave elimination according to claim 6, characterized in that: The S71 includes the following steps: S711: at each time step, interpolating the flow velocity u and the liquid volume fraction α on each interpolation sampling line by using a set sampling interpolation method; S712: on the main processor, obtain the interpolation results on the other processors and the coordinates of the interpolation points on the corresponding sampling lines through the gatherInplace function in the globalIndex class; S713: on the main processor, the interpolation results and the corresponding interpolation point coordinates obtained on the remaining processors are sorted by using the sorting mapping in the gatherSort function in the coordSet class; S72 includes the following steps: S721: Integrate u and α on the main processor by using the trapezoidal integration method to obtain the instantaneous vertical average flow velocity u(t) and the instantaneous water depth on each interpolation sampling line, and take the absolute value of the difference between the instantaneous water depth and the static water depth as the instantaneous wave surface displacement A(t); S722: Send u(t) and A(t) calculated on the main processor to the remaining processors.

9. A parameterized calculation method for flexible vegetation wave elimination according to claim 6, characterized in that: The S8 comprises the following steps: S81: construct a body scalar matrix vegetationCell, in which each element corresponds one-to-one to each element in cellFaceAddressing, and the initial value of each element is 0; S82: Determine whether the height of the center of the grid in the grid cluster corresponding to each element in cellFaceAddressing is less than the equivalent height of the vegetation at the corresponding interpolation sampling line If so, set the corresponding element of vegetationCell to 1.

10. A parameterized calculation method for flexible vegetation wave elimination according to claim 9, characterized in that: In S9, the drag coefficient is calculated by the following formula: : ， ， In the formula, A, B and C are constants, and μ is the dynamic viscosity of water.

Images