Method and program product for automatic generation of a three-dimensional geological finite element model

Abstract

The application discloses a method and a program product for automatically generating a three-dimensional geological finite element model, which determines a horizontal reference surface based on the maximum height difference of the ground surface, determines a first grid size and a second grid size based on a shear wave velocity, and then comprehensively determines the number of ground surface part grids and divides hexahedral units; according to the first grid size, equally divides a base part into horizontal slices with the same height, and determines the grid size of the horizontal slices in combination with a corrected minimum shear wave velocity; according to the grid size and thickness requirement, combines adjacent horizontal slices to obtain a final rock layer division scheme, and then comprehensively determines the number of rock layer grids and divides the rock layer into hexahedral units; and according to the three-dimensional center coordinates of each hexahedral unit, compares the three-dimensional center coordinates with an original geological model, and sets the stratum category of the hexahedral unit. The application greatly reduces the modeling time and cost of a three-dimensional non-uniform field finite element model of a complex terrain, has high grid quality, and is beneficial to the application of the three-dimensional geological model in the simulation and analysis of earthquake disasters.

Description

Technical Field

[0001] This invention relates to the field of finite element preprocessing technology in geotechnical engineering, and in particular to a method and program product for automatically generating three-dimensional geological finite element models. Background Technology

[0002] Existing earthquake damage surveys and studies have shown that local topographic relief and media inhomogeneity significantly impact seismic wave propagation and building damage. Meanwhile, many cities have established 3D geological models for urban planning and disaster prevention. However, these 3D geological models cannot be directly used for finite element dynamic analysis. Current methods involve importing the 3D geological model into finite element preprocessing software (such as Hypermesh, SolidWorks, and Midasgts) for mesh generation. However, due to the difficulty in controlling the quality of 3D meshes under complex terrain conditions, this preprocessing process requires a high level of technical skill from operators, is cumbersome and complex, difficult to automate, and time-consuming. Summary of the Invention

[0003] This invention provides a method and program product for automatically generating three-dimensional geological finite element models, which can solve the problem of automatic mesh generation in three-dimensional non-uniform site finite element modeling under complex terrain conditions.

[0004] To achieve the above objectives, the present invention adopts the following technical solution:

[0005] The specific steps for automatically generating three-dimensional geological finite element models are as follows:

[0006] Step 1: Based on the maximum elevation difference of the surface in the 3D geological model, determine the horizontal datum and the 3D coordinate system; where the area above the horizontal datum is the surface part, and the area below the horizontal datum is the base part.

[0007] Step 2: Based on the shear wave velocities of different strata in the 3D geological model, determine the first grid size by the minimum shear wave velocity and the second grid size by the maximum shear wave velocity; based on the first and second grid sizes and combined with the maximum surface elevation difference, comprehensively determine the number of grids in the x, y, and z directions for the surface portion.

[0008] Step 3: Divide the surface portion of the 3D geological model into rock columns based on the number of grids in the x and y directions. Then, combining the number of grids in the z direction, divide each rock column into eight-node hexahedral elements from top to bottom; this completes the meshing of the surface portion.

[0009] Step 4: Based on the first grid size of the surface part, divide the base part into horizontal slices of the same height; count the minimum shear wave velocity of different strata in each horizontal slice from top to bottom and sort them; correct the sorted minimum shear wave velocity so that the minimum shear wave velocity of each layer is not lower than the minimum shear wave velocity of the layer above.

[0010] Step 5: Determine the horizontal slice grid size based on the corrected minimum shear wave velocity; merge adjacent horizontal slices according to the horizontal slice grid size to obtain the final rock layer division scheme for the base section;

[0011] Step 6: Based on the thickness and grid size of each rock layer in the base section and combined with the first grid size of the surface section, comprehensively determine the number of grids in the x, y, and z directions for each rock layer;

[0012] Step 7: Divide each rock layer in the base into equal parts in the x, y, and z directions according to the determined number of meshes, and divide them into eight-node hexahedral elements; then, use a two-in-one method to divide transition elements at the top of each rock layer to complete the transition of meshes of different sizes in the x and y directions.

[0013] Step 8: Based on the three-dimensional centroid coordinates of each hexahedral element, compare them with the original geological model and set its stratigraphic category to be the same as the stratigraphic category corresponding to the centroid position of each element in the original geological model; thus, the automatic generation of the three-dimensional geological finite element model is completed.

[0014] Furthermore, in step one, the maximum surface elevation difference H1 = H is calculated based on the three-dimensional geological model. max -H min ;

[0015] Among them, H max and H min These are the elevations of the highest and lowest points on the Earth's surface, respectively.

[0016] Select an altitude of H min The horizontal plane of H1 is used as the reference plane, and a three-dimensional Cartesian coordinate system is established with this plane as the xy plane. The z-axis is positive upwards and negative downwards.

[0017] Furthermore, in step two, based on the minimum shear wave velocity V of different strata on the surface... s,min Calculate the first mesh size S1:

[0018] S1=k1×V s,min / f max

[0019] k1 is an empirical constant, taken as 1 / 8; f max To effectively simulate the upper limit of the frequency band;

[0020] Then, based on the maximum shear wave velocity V of different strata on the surface... s,max Calculate the second mesh size S2:

[0021]

[0022] Indicates rounding down;

[0023] The number of grids N in the x, y, z directions of the surface region is calculated using the following formula. x1 N y1 and N z1 :

[0024]

[0025]

[0026]

[0027] W and T represent the length and width of the three-dimensional geological model in the x and y directions, respectively. This indicates rounding up to the nearest integer.

[0028] Furthermore, for step three, the three-dimensional geological model is divided into N equal parts in the xy plane. x1 ×N y1 Each rock column was analyzed, and the coordinates of the four vertices on the top surface of each column were obtained by interpolating the elevation data from the ground surface. Then, the four sides of each rock column were divided into N equal parts. z1 Each rock column is divided into N sections from top to bottom. z1 Each element consists of an eight-node hexahedral cell; once completed, the surface portion is divided into N... x1 ×N y1 ×N z1 An eight-node hexahedral element.

[0029] Furthermore, in step four, the base portion is divided into m horizontal slices of equal height in the z-direction. m can be calculated using the following formula.

[0030]

[0031] Z b Let V be the z-coordinate of the bottom of the 3D geological model; then, from top to bottom, the minimum shear wave velocity of different strata in each horizontal slice is counted and denoted as V. s,1 V s,2 V s,3 ...V s,m .

[0032] Furthermore, regarding V s,1 V s,2V s,3 ...V s,m The values ​​of V are corrected: starting from the bottom layer, compare V... s,m and V s,m-1 The size of V s,m-1 >V s,m Then V s,m-1 The value is corrected to V s,m The value of V; if V s,m-1 ≤V s,m If so, no adjustment is needed; then compare V. s,m-1 and V s,m-2 The size of V s,m-2 >V s,m-1 Then V s,m-2 The value is corrected to V s,m-1 The value of V; if V s,m-2 ≤V s,m-1 If the value is not specified, no adjustment is needed; and so on, to ensure that the corrected shear wave velocity satisfies V. s,1 ≤V s,2 ≤V s,3 ≤...≤.V s,m .

[0033] Furthermore, in step five, the grid size L of each horizontal slice is determined according to the following formula. i :

[0034]

[0035] Next, adjacent horizontal slices with the same grid size are merged into rock layers; then, from top to bottom, rock layers with a thickness less than three times the grid size are successively incorporated into the rock layers above them, resulting in the final rock layer division scheme for the base part.

[0036] Furthermore, in step six, the number of grids for each rock layer in the x, y, and z directions is calculated using the following formula:

[0037] N x2 =N x1 ×S1 / L

[0038] N y2 =N y1 ×S1 / L

[0039]

[0040] Where D and L represent the rock layer thickness and the rock layer grid size, respectively.

[0041] A computer program product includes a computer program / instructions that, when executed by a processor, implement the above-described method steps for automatically generating a three-dimensional geological finite element model.

[0042] The beneficial effects of this invention are reflected in:

[0043] 1) By dividing the horizontal reference plane, this invention facilitates the division of the three-dimensional geological model into two parts: the surface and the base. The two parts can be meshed separately, which takes into account both the meshing quality and the simplification of meshing.

[0044] 2) This invention, through the division of the surface part into grids, based on shear wave velocity sorting, horizontal rock column division, and vertical division, ensures the accuracy and operability of the surface part grid division into hexahedral units.

[0045] 3) This invention, through the division of the base part of the mesh, based on shear wave velocity correction, rock layer merging and other technologies, ensures the simplification and effective transition of the mesh, and ensures the rapid division and convenient operation of hexahedral elements.

[0046] This invention significantly reduces the time and cost of finite element modeling for 3D non-uniform sites under complex terrain conditions, facilitating the application of 3D geological models in earthquake disaster simulation and analysis. Furthermore, the 3D finite element model obtained by this method is entirely composed of eight-node hexahedral elements, resulting in high mesh quality and reducing the likelihood of distortion and instability. Other features and advantages of this invention will be set forth in the following description and will be apparent in part from the description, or may be learned by practicing the invention; the main objectives and other advantages of the invention may be realized and obtained by means of the methods particularly pointed out in the description. Attached Figure Description

[0047] Figure 1 This is a three-dimensional geological map of a river valley;

[0048] Figure 2 It is an automatically generated 3D finite element model corresponding to a 3D geological map of a river valley. Detailed Implementation

[0049] Taking a river valley site with a horizontal projected area of ​​500m × 500m as an example, such as Figure 1 As shown, the three-dimensional geological model of the site contains three strata: Quaternary sediments in the first layer, silty mudstone in the second layer, and slate in the third layer, with shear wave velocities of 400 m / s, 800 m / s, and 1600 m / s, respectively. The maximum and minimum elevations at the surface are 11.65 m and 73.15 m, respectively, and the elevation at the bottom of the model is -149.85 m.

[0050] Combination Figure 1 and Figure 2 As shown, the method for automatically generating a three-dimensional geological finite element model is further explained. The specific steps are as follows:

[0051] Step 1: Based on the maximum elevation difference of the surface in the 3D geological model, determine the horizontal datum and the 3D coordinate system; where the part above the horizontal datum is the surface part, and the part below the horizontal datum is the base part.

[0052] In step one, the maximum surface elevation difference H1 = H is calculated based on the three-dimensional geological model. max -H min ;

[0053] Among them, H max and H min These are the elevations of the highest and lowest points on the Earth's surface, respectively.

[0054] Select an altitude of H min The horizontal plane of H1 is used as the reference plane, and a three-dimensional Cartesian coordinate system is established with this plane as the xy plane. The z-axis is positive upwards and negative downwards.

[0055] Based on step one, H1 = 61.5m can be determined, and a horizontal plane with an altitude of -49.85m is selected as the reference plane.

[0056] Step 2: Based on the shear wave velocities of different strata in the 3D geological model, determine the first grid size by the minimum shear wave velocity and the second grid size by the maximum shear wave velocity; based on the first and second grid sizes and combined with the maximum surface elevation difference, comprehensively determine the number of grids in the x, y, and z directions for the surface portion.

[0057] In step two, based on the minimum shear wave velocity V of different strata on the surface... s,min Calculate the first mesh size S1:

[0058] S1=k1×V s,min / f max

[0059] k1 is an empirical constant, taken as 1 / 8; f max To effectively simulate the upper limit of the frequency band;

[0060] Then, based on the maximum shear wave velocity V of different strata on the surface... s,max Calculate the second mesh size S2:

[0061]

[0062] Indicates rounding down;

[0063] The number of grids N in the x, y, z directions of the surface region is calculated using the following formula. x1 N y1 and N z1 :

[0064]

[0065]

[0066]

[0067] W and T represent the length and width of the three-dimensional geological model in the x and y directions, respectively. This indicates rounding up to the nearest integer.

[0068] According to step two, f max With an effective frequency band upper limit of 10Hz, the first and second grid sizes can be determined as S1 = 5m and S2 = 20m, respectively. The number of grids N in the x, y, and z directions for the surface portion is... x1 =100, N y1 =100 and N z1 =25.

[0069] Step 3: Divide the surface portion of the 3D geological model into rock columns based on the number of grids in the x and y directions. Then, based on the number of grids in the z direction, divide each rock column into eight-node hexahedral units from top to bottom; this completes the grid division of the surface portion.

[0070] For step three, the three-dimensional geological model is divided into N equal parts in the xy plane. x1 ×N y1 Each rock column was analyzed, and the coordinates of the four vertices on the top surface of each column were obtained by interpolating the elevation data from the ground surface. Then, the four sides of each rock column were divided into N equal parts. z1 Each rock column is divided into N sections from top to bottom. z1 Each element consists of an eight-node hexahedral cell; once completed, the surface portion is divided into N... x1 ×N y1 ×N z1 Each hexahedral element is divided into 100×100×25 eight-node hexahedral elements. Each hexahedral element has a size of approximately S1 in the x and y directions, and a size in the z direction between 0.5S1 and S1.

[0071] Step 4: Based on the first grid size of the surface portion, divide the base portion into horizontal slices of equal height; count and sort the minimum shear wave velocities of different strata in each horizontal slice from top to bottom, and correct the sorted minimum shear wave velocities so that the minimum shear wave velocity of each layer is not lower than the minimum shear wave velocity of the layer above.

[0072] In step four, the base portion is divided into m horizontal slices of equal height along the z-direction. m can be calculated using the following formula.

[0073]

[0074] Z b Let V be the z-coordinate of the bottom of the 3D geological model; then, from top to bottom, the minimum shear wave velocity of different strata in each horizontal slice is counted and denoted as V. s,1 V s,2 V s,3 ...V s,m .

[0075] For V s,1 V s,2 V s,3 ...V s,m The values ​​of V are corrected: starting from the bottom layer, compare V... s,m and V s,m-1 The size of V s,m-1 >V s,m Then V s,m-1 The value is corrected to V s,m The value of V; if V s,m-1 ≤V s,m If so, no adjustment is needed; then compare V. s,m-1 and V s,m-2 The size of V s,m-2 >V s,m-1 Then V s,m-2 The value is corrected to V s,m-1 The value of V; if V s,m-2 ≤V s,m-1 If the value is not specified, no adjustment is needed; and so on, to ensure that the corrected shear wave velocity satisfies V. s,1 ≤V s,2 ≤V s,3 ≤...≤.V s,m .

[0076] Among them, Z b =-100m. According to the above formula, we can get m=20. After cutting the base part into 20 horizontal slices, we can find that the minimum shear wave velocities of different strata in the upper 8 horizontal slices and the lower 12 horizontal slices are 800m / s and 1600m / s, respectively.

[0077] Step 5: Determine the horizontal slice grid size based on the corrected minimum shear wave velocity; merge adjacent horizontal slices according to the horizontal slice grid size to obtain the final rock layer division scheme for the base section.

[0078] In step five, the grid size L of each horizontal slice is determined according to the following formula. i :

[0079]

[0080] Next, adjacent horizontal slices with the same grid size are merged into rock layers; then, from top to bottom, rock layers with a thickness less than three times the grid size are successively incorporated into the rock layers above them, resulting in the final rock layer division scheme for the base part.

[0081] Step Six: Based on the thickness and grid size of each rock layer in the base section, and combined with the first grid size of the surface section, determine the number of grids in the x, y, and z directions for each rock layer. Calculate the number of grids in the x, y, and z directions for each rock layer using the following formula:

[0082] N x2 =N x1 ×S1 / L

[0083] N y2 =N y1 ×S1 / L

[0084]

[0085] Where D and L represent the rock layer thickness and the rock layer grid size, respectively.

[0086] In this embodiment, according to step five, the base portion comprises two rock layers. The rock layer thicknesses from top to bottom are 40m and 60m, respectively, corresponding to rock layer grid sizes of 10m and 20m. According to the formula in step six, the number of grids in the x, y, and z directions for the rock layer with a grid size of 40m are 50, 50, and 4, respectively; the number of grids in the x, y, and z directions for the rock layer with a grid size of 60m are 25, 25, and 3, respectively.

[0087] Step 7: Divide each rock layer in the base into equal parts in the x, y, and z directions according to the determined number of meshes, and divide them into eight-node hexahedral elements; then, use a two-in-one method to divide transition elements at the top of each rock layer to complete the transition of meshes of different sizes in the x and y directions.

[0088] Step 8: Based on the three-dimensional centroid coordinates of each hexahedral element, compare them with the original geological model and set its stratigraphic category to be the same as the stratigraphic category corresponding to the centroid position of each element in the original geological model; thus completing the automatic generation of the three-dimensional geological finite element model. The final three-dimensional finite element model is obtained based on steps 7 and 8, as follows: Figure 2 As shown. The three-dimensional finite element model obtained by this method is composed entirely of eight-node hexahedral elements, and the ratio of the longest side to the shortest side does not exceed 2. The mesh quality is high and it is not prone to distortion and instability.

[0089] In this embodiment, computer program code for performing the operations of the present invention can be written in one or more programming languages ​​or a combination thereof. Programming languages ​​include object-oriented programming languages—such as Java, Smalltalk, and C++—as well as conventional procedural programming languages—such as the "C" language or similar programming languages. The program code can be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computer (e.g., via the Internet using an Internet service provider).

[0090] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for automatically generating three-dimensional geological finite element models, characterized in that, The specific steps are as follows: Step 1: Based on the maximum elevation difference of the surface in the 3D geological model, determine the horizontal datum and the 3D coordinate system; where the area above the horizontal datum is the surface part, and the area below the horizontal datum is the base part. Step 2: Based on the shear wave velocities of different strata in the 3D geological model, determine the first grid size using the minimum shear wave velocity and the second grid size using the maximum shear wave velocity; based on the first and second grid sizes and combined with the maximum surface elevation difference, comprehensively determine the surface portion... x , y , z The number of grid cells in the three directions; In step two, based on the minimum shear wave velocity of different strata on the surface... V s,min Calculate the first grid size S 1: S 1= k 1× V s,min / f max ； k 1 is an empirical constant, taken as 1 / 8; f max To effectively simulate the upper limit of the frequency band; Then, based on the maximum shear wave velocity of different strata on the surface... V s,max Calculate the second grid size S 2: ； " " indicates rounding down; Step 3, according to x and y The number of grids in each direction divides the surface portion of the 3D geological model into equal rock columns, and then combines... z The number of grids in each direction divides each rock column into eight-node hexahedral units from top to bottom; this completes the grid division of the surface portion. Step 4: Based on the first grid size of the surface part, divide the base part into horizontal slices of the same height; count the minimum shear wave velocity of different strata in each horizontal slice from top to bottom and sort them; correct the sorted minimum shear wave velocity so that the minimum shear wave velocity of each layer is not lower than the minimum shear wave velocity of the layer above. Step 5: Determine the horizontal slice grid size based on the corrected minimum shear wave velocity; merge adjacent horizontal slices according to the horizontal slice grid size to obtain the final rock layer division scheme for the base section; In step five, the grid size for each horizontal slice is determined using the following formula. L i : ； Next, adjacent horizontal slices with the same grid size are merged into rock layers; then, from top to bottom, rock layers with a thickness less than three times the grid size are successively merged into the rock layers above them, resulting in the final rock layer division scheme for the base part; Step Six: Based on the thickness and grid size of each rock layer in the base section, and combined with the first grid size of the surface section, comprehensively determine the location of each rock layer. x , y and z Number of grid cells in a direction; Step 7: For the base section, each rock layer is arranged according to the determined grid number. x , y , z The rock strata are divided into eight-node hexahedral elements by equal division along the direction; then, transition elements are divided at the top of each stratum using a two-in-one method, completing the division of elements of different sizes in... x and y Transition of grid direction; Step 8: Based on the three-dimensional centroid coordinates of each hexahedral element, compare them with the original geological model and set its stratigraphic category to be the same as the stratigraphic category corresponding to the centroid position of each element in the original geological model; thus, the automatic generation of the three-dimensional geological finite element model is completed.

2. The method for automatically generating a three-dimensional geological finite element model as described in claim 1, characterized in that, In step one, the maximum surface elevation difference is calculated based on the three-dimensional geological model. H 1= H max - H min ; in, H max and H min These are the elevations of the highest and lowest points on the Earth's surface, respectively. Select altitude as H min - H The horizontal plane is taken as the reference plane, and this plane is used as... xy Establish a three-dimensional Cartesian coordinate system on a plane. z Upward is positive, downward is negative.

3. The method for automatically generating a three-dimensional geological finite element model as described in claim 2, characterized in that, In step two, the surface portion is calculated according to the following formula. x , y , z Number of grids in the direction N x1 , N y1 and N z1 : ； W and T These respectively represent the three-dimensional geological model in x and y The length and width in the direction, with "" indicating rounding up.

4. The method for automatically generating a three-dimensional geological finite element model as described in claim 3, characterized in that, For step three, in xy Divide the three-dimensional geological model into equal parts in the plane. N x1 × N y1 Each rock column was analyzed, and the coordinates of the four vertices on the top surface of each column were obtained by interpolating the elevation data from the ground surface. Then, the four sides of each rock column were divided into equal parts. N z1 Divide each rock column into sections, and then divide each rock column from top to bottom into sections. N z1 An eight-node hexahedral element; Once everything is complete, the surface portion will be divided into... N x1 × N y1 × N z1 An eight-node hexahedral element.

5. The method for automatically generating a three-dimensional geological finite element model as described in claim 3, characterized in that, In step four, z The base portion is divided into sections in the direction of the direction. m A horizontal slice of the same height, m It can be calculated using the following formula ； Z b The bottom of the three-dimensional geological model z Coordinates; then, from top to bottom, the minimum shear wave velocities of different strata in each horizontal slice are counted and denoted as follows: V s,1 , V s,2 , V s,3 ... V s,m .

6. The method for automatically generating a three-dimensional geological finite element model as described in claim 5, characterized in that, right V s,1 , V s,2 , V s,3 ... V s,m The values ​​are corrected: starting from the bottom layer, comparisons are performed. V s,m and V s,m-1 The size, if V s,m-1 > V s,m Then V s,m-1 The value was corrected to V s,m The value; if V s,m-1 ≤ V s,m If so, no adjustment is needed; then compare again. V s,m-1 and V s,m-2 The size, if V s,m-2 > V s,m-1 Then V s,m-2 The value was corrected to V s,m-1 The value; if V s,m-2 ≤ V s,m-1 If the value is 0, no adjustment is needed; and so on, to ensure that the corrected shear wave velocity meets the requirements. V s,1 ≤ V s,2 ≤ V s,3 ≤...≤. V s,m .

7. The method for automatically generating a three-dimensional geological finite element model as described in claim 6, characterized in that, For step six, calculate the values ​​of each rock layer according to the following formula. x , y and z Number of grids in each direction: N x2 = N x1 × S 1 / L N y2 = N y1 × S 1 / L ； in, D and L These represent the rock layer thickness and the rock layer grid size, respectively.

8. A computer program product comprising a computer program / instructions, characterized in that: When the computer program / instructions are executed by the processor, they implement the steps of the method described in any one of claims 1 to 7.

Images