Complex terrain processing method and system suitable for three-dimensional non-hydrostatic atmospheric model

Abstract

The application discloses a complex terrain processing method and system suitable for a three-dimensional non-hydrostatic atmospheric model, which comprises the following steps: adopting preset regional terrain data and an atmospheric background field, and performing atmospheric model preprocessing on the terrain data; performing key grid identification on the preset region to obtain a grid distribution, marking terrain nearby grid points according to the grid distribution and recording distance and slope information, obtaining a marked grid distribution, dividing the marked grid distribution according to a process number and calculating a sub-region, and performing adaptive grid encryption on the sub-region; updating a forecast variable value at a next moment by updating an atmospheric equation group of the sub-region after high-order numerical discretization, generating a logarithmic velocity profile based on the Monin-Obukhov similarity theory, processing an atmospheric velocity of a boundary layer according to the logarithmic velocity profile to obtain an ideal velocity distribution curve, and outputting atmospheric motion characteristics under a complex terrain by correcting a forecast variable on a grid point near a lower surface boundary.

Claims

1. A method for processing complex terrain in a three-dimensional non-hydrostatic atmospheric model, characterized in that, Includes the following steps: The terrain data is preprocessed using terrain data and atmospheric background field of a preset area; the terrain data includes longitude, latitude, altitude information, and terrain elevation data. The preset area is identified by key grids to obtain the grid distribution. Based on the grid distribution, grid points near the terrain are marked as velocity reconstruction points. Distance and slope information are recorded to obtain the marked grid distribution. The marked grid distribution is divided according to the number of processes and sub-regions are calculated. The high-order numerical discretization of the atmospheric motion equations in the sub-region is performed and updated to obtain the predicted variable values ​​for the next time step. Logarithmic velocity profiles are generated based on the Moning-Obukhov similarity theory. The boundary layer atmospheric velocity is processed according to the logarithmic velocity profiles to obtain an ideal velocity distribution curve. By constructing a local tangential / normal coordinate system, interpolating the interpolation points of the terrain normal azimuth, reconstructing the velocity vector by combining the Moning-Obukhov similarity theory, applying boundary conditions by converting to Cartesian coordinates, and participating in time iteration, the atmospheric motion characteristics under complex terrain are output; the boundary conditions include optimizing near-ground velocity.

2. The method for processing complex terrain in a three-dimensional non-hydrostatic atmospheric model according to claim 1, characterized in that, A method for obtaining the grid distribution by identifying key grids in the preset area includes: Based on the terrain data and latitude and longitude coordinates of the preset area, the terrain height is extracted into an orthogonally distributed Cartesian coordinate system through an interpolation algorithm. The data is then retrieved column by column from low to high along the horizontal direction. Points with a vertical height less than or equal to the terrain height are marked as internal points of the terrain, i.e., solid region nodes. The remaining points are fluid region nodes. Calculate terrain slope and aspect using a window: ； ； in For grid points The slope of the terrain, The difference in elevation between east and west is the topographical elevation difference. The difference in terrain elevation runs from north to south. This represents the ratio of the east-west terrain elevation difference to the lateral grid spacing. This represents the ratio of the north-south terrain elevation difference to the longitudinal grid spacing. For grid points The slope direction; Calculate the normal distance from the center of the non-terrain interior grid point to the terrain surface: ； in To minimize the impact distance, Let be the geometric height of the center of the k-th layer mesh. For grid points The terrain elevation, is the normal distance from the center of the k-th non-terrain interior point grid to the terrain surface; like When, the current grid point is a point inside the terrain; when Then calculate the normal distance. If 0 Then the current grid point is a terrain boundary point, if If the terrain slope is greater than the slope threshold, then the current grid point is a point influenced by the terrain. Output the points inside the terrain, the points on the terrain boundary, and the points affected by the terrain as a grid distribution.

3. The method for processing complex terrain in a three-dimensional non-hydrostatic atmospheric model according to claim 1, characterized in that, The method for obtaining the marked grid distribution includes: Mark points within the terrain based on terrain boundary points and near-terrain influence points. Check the horizontal and vertical neighborhood grids of the terrain boundary points. If the marked grid points are less than [a certain distance] from the terrain location... Then, the upper-level mesh will be marked; where The vertical resolution of the pattern; A marked grid distribution is generated based on grid coordinates, grid type labels, and terrain parameters, including slope, aspect, and normal distance.

4. The method for processing complex terrain in a three-dimensional non-hydrostatic atmospheric model according to claim 1, characterized in that, The method for dividing the marked grid distribution and calculating sub-regions based on the number of processes includes: The entire simulation computation domain is decomposed according to the number of parallel processes allocated, with each process responsible for a sub-region, and information is exchanged between sub-regions through a set grid width; First, obtain the parallel computing parameters, assign load weights to the horizontal grid, and combine terrain complexity and physical process overhead; the expression is: ； in For grid points The slope of the terrain, For key grid density, This is the slope weighting coefficient. Density weighting coefficient, The load weights for the horizontal grid; The total load of the horizontal grid is obtained by summing the load weights assigned to the horizontal grid and then using the ratio of the total load of the horizontal grid to the total number of processes. A one-dimensional scan of the horizontal grid along longitude or latitude is performed to obtain the initial row and column groupings. The load weight is accumulated from the grid starting point. When the accumulated value reaches the total load of the horizontal grid, the partition boundary is marked to generate the initial horizontal sub-region. If the load deviation of the last partition exceeds 10%, a secondary balance is performed by adjusting the boundary grid assignment. Each sub-region is extended by 2 layers of grid in the horizontal direction as an overlap area, where the overlap width is determined according to the horizontal accuracy of the model. The proportion of critical grids in each vertical grid layer is calculated. If the proportion of critical grids in lower layers is high, a higher vertical weight is assigned, and multiple vertical layers are allocated to each process according to their weight ratio. Record the horizontal and vertical coordinate ranges for each process, including overlapping grid areas, and calculate the number of various grid types within each sub-region: ； ； ； in The total number of grids for process p. The number of critical grids for process p. The critical grid percentage for process p. Let P be the starting coordinates of process p. Let p be the termination coordinates. Let p be the vertical layer number of process p. Let p be the starting level of process p. Let p be the termination level of process p. The grid type is marked as 1 for ordinary grids, 1 for terrain boundary points, and 2 for near-terrain influence points; Output the process ID, coordinate range, number of grids, and percentage of critical grids as the configuration information for the sub-region.

5. The method for processing complex terrain in a three-dimensional non-hydrostatic atmospheric model according to claim 1, characterized in that, The method for adaptive mesh densification and load balancing optimization of the sub-region includes: Obtain the encryption level, encryption judgment threshold, and load balancing target; If a grid point within a sub-region is a terrain boundary point, then L1 level densification is triggered; if the horizontal velocity shear of a grid point within a sub-region is greater than or equal to 0.1 / s and the turbulent kinetic energy is greater than or equal to 1m³, then L1 level densification is triggered. 2 / s 2 If the conditions are met, then L1-level encryption is triggered; if the grid points within the sub-region simultaneously satisfy both terrain and physical process conditions, then L2-level encryption is triggered. A multi-layer nested grid algorithm is used to ensure the nesting relationship between the encrypted area and the coarse grid area, and a minimum encryption size is set; the minimum encryption size is that the horizontal range of a single-level encrypted grid is not less than 5×5 grids. For regions that meet the encryption conditions, fine meshes are generated in layers: for L1 level, the horizontal / vertical mesh spacing of the basic mesh is reduced to 1 / 2; for L2 level, it is further reduced to 1 / 2 on the basis of L1, and vertical layering is synchronously encrypted. Two transition zones are set at the junction of coarse and fine grids, and bilinear interpolation is used to transfer fluxes: for momentum and heat fluxes, they are calculated in a mixed manner according to grid level weights within the transition zone; for scalar fields, conservation interpolation is used to ensure that the mass conservation error is less than or equal to the mass conservation error threshold. Calculate the complexity index of the sub-region, divide the initial sub-regions based on the terrain complexity, and allocate more processes to the sub-region when the slope variance is greater than the slope variance threshold; where the terrain complexity is the slope variance. A load check is performed every 10 time steps. If the maximum deviation of the complexity exponent between processes exceeds the deviation threshold, region migration is initiated: The number of grid blocks to be migrated for overloaded processes is calculated; non-critical grids in the overloaded region are migrated to the lightly loaded process, prioritizing the retention of grid blocks containing terrain boundary points; where the overloaded process is... , This represents the average load.

6. The method for processing complex terrain in a three-dimensional non-hydrostatic atmospheric model according to claim 1, characterized in that, A method for updating the forecast variable values ​​for the next time step by discretizing the atmospheric equations for the sub-region using high-order numerical methods includes: Obtain the initial time prediction variables of the sub-region The initial forecast variables include three-dimensional velocity components, potential temperature, and air pressure; the parameters of the complex terrain processing algorithm include horizontal grid spacing, vertical grid spacing, Coriolis parameters, gravitational acceleration, and terrain distance and slope information calculated in the preprocessing stage. Given a set of core governing equations including dynamic equations, continuity equations, and thermodynamic energy equations, where the set of equations is: ； in The three-dimensional velocity components in Cartesian coordinates. and For atmospheric density and atmospheric density disturbances, For potential temperature, For air pressure, It is gravitational acceleration. Let be the basis vectors of the three-dimensional Cartesian coordinate system. Coriolis force, For subgrid turbulent stress, the right-hand side of the equation of motion is The force term represents the force exerted by the terrain boundary on the atmosphere; a multi-moment constrained finite volume scheme is used on the horizontal mesh. The system achieves fourth-order spatial accuracy. To address the complexity of the terrain in the vertical direction, a high-order conserved finite difference scheme is adopted, and the order is automatically reduced near the terrain. For terrain boundary points, the scheme is reduced to a second-order upwind scheme. The time integration method uses an explicit third-order Runge-Kutta method to divide the forecast variable update into three sub-steps: updating arbitrary forecast variables according to the time step size constrained by CFL conditions; The three-dimensional forecast variable field of the sub-region at the next time step , , , , The output is the predicted variable value for the next time step.

7. The method for processing complex terrain in a three-dimensional non-hydrostatic atmospheric model according to claim 1, characterized in that, The method for generating the logarithmic velocity profile includes: Obtain surface observation data, surface parameters, and meteorological constants. Under neutral atmospheric conditions, calculate the friction velocity using wind speed at a reference altitude. Calculate the sensible heat flux using the temperature difference between the reference altitude and the surface. Calculate the Monin-Obukhov length. The expression is: ； ； ； in It is a false fever. For aerodynamic impedance, For friction speed, The specific heat of dry air, is the von Kármán constant. For reference height, The length of the surface roughness. Zero planar displacement height For reference temperature, For sensible heat flux, The length of the Monin-Obukhov pair. It is the acceleration due to gravity. For surface temperature, air density; To calculate the dimensionless height of the target, the basic outline formula is given as follows: ； in The height is dimensionless. For target height; when When, under neutral conditions, the similarity function is ;when When the condition is unstable, the similarity function is: ,when When the stability condition is met, the similarity function is: ; Calculate the wind speed at any altitude: ； in Let z be the wind speed at any height. For the similarity function correction term in integral form, by... Points earned; When the condition is neutral, the similarity function correction term is 0; when the condition is unstable: Under stable conditions: The basic contour line is given based on the similarity function correction term under different conditions.

8. The method for processing complex terrain in a three-dimensional non-hydrostatic atmospheric model according to claim 1, characterized in that, A method for processing boundary layer atmospheric velocities based on the logarithmic velocity profile to obtain an ideal velocity distribution curve includes: First, based on the actual data from the grid near the ground, the velocity vector values ​​at the interpolation points are interpolated. Inverse distance weighted interpolation is used, and the expression is: ； in For the velocity of the interpolation point, It is the velocity at the nodes near the interpolation point. This represents the number of neighboring nodes retrieved by the algorithm; its size varies depending on the terrain distribution. It is the distance from the interpolation point to nearby nodes. The distance to the farthest node. It's weight. ; By modifying the velocity of the first layer nodes outside the boundary using the velocity of the interpolation points, i.e. the velocity vector of the velocity reconstruction points, the velocity near the boundary is made to achieve the expected distribution.

9. A complex terrain processing system suitable for three-dimensional non-hydrostatic atmospheric models, for performing the method according to any one of claims 1-8, characterized in that, It includes a terrain processing module, a parallel processing module, a non-hydrostatic atmospheric module, a boundary layer processing module, and an immersion algorithm module; The terrain processing module automatically marks the location of points near complex terrain and other information based on different grid distributions and algorithms; The parallel processing module divides the computing region according to the number of processes, and each sub-process is responsible for the corresponding sub-region. The non-hydrostatic atmospheric module completes the spatial discretization and time integration steps of the atmospheric equations to obtain the predicted variable values ​​for the next moment. This is done in the spatial discretization module and the time integration module, respectively. Depending on the numerical discretization method, the spatial discretization module can be further divided into a horizontal discretization module and a vertical discretization module. The time integration module is also divided into a horizontal time integration module and a vertical time integration module depending on the time integration scheme. The boundary layer processing module processes the atmospheric velocity in the boundary layer based on the atmospheric model dynamics framework. The immersion boundary algorithm module corrects the predicted variables at grid points near the underlying surface boundary; the sub-modules of the immersion boundary algorithm module include a local coordinate transformation module, a spatial interpolation module, a reconstructed vector module, an applied boundary module, and a time integration module.

Images