CN113204905B - Numerical simulation method for finite element by contact induced polarization method - Google Patents
Numerical simulation method for finite element by contact induced polarization method Download PDFInfo
- Publication number
- CN113204905B CN113204905B CN202110497481.1A CN202110497481A CN113204905B CN 113204905 B CN113204905 B CN 113204905B CN 202110497481 A CN202110497481 A CN 202110497481A CN 113204905 B CN113204905 B CN 113204905B
- Authority
- CN
- China
- Prior art keywords
- grid
- finite element
- polarization
- induced polarization
- subdivision
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Geometry (AREA)
- Theoretical Computer Science (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Computer Graphics (AREA)
- Remote Sensing (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention belongs to the technical field of engineering, hydrology and environmental geophysical exploration, and discloses a finite element numerical simulation method of a contact induced polarization method. The invention solves the difficult problem of the traditional finite element numerical simulation calculation method of the lossy dielectric contact induced polarization method with any shape, has simple algorithm structure, simple and convenient realization process, reasonable design and stable and reliable calculation result, can be widely applied to the numerical simulation calculation of the exploration polarizability of the lossy dielectric with any shape related to the engineering, hydrology and environmental geophysical exploration field, and provides method support for the engineering, hydrology and environmental geophysical exploration numerical simulation, inversion and explanation based on complex lossy dielectric power supply.
Description
Technical Field
The invention belongs to the technical field of engineering, hydrology and environmental geophysical exploration, and particularly relates to a finite element numerical simulation method of a contact induced polarization method.
Background
At present, the charging method utilizes natural or artificially exposed dew point of good conductor and dew point of underground water to directly connect a power supply electrode (generally positive electrode), and simultaneously places another power supply electrode at a position meeting the requirement of infinity, and observes the change and distribution rule of a charging electric field through two measuring electrodes to estimate the physical property and the spatial distribution of the good conductor. Because the shape of the charging medium range is similar to that of the observed potential abnormal range, the actual charging medium range can be inferred according to the observed potential abnormality, and the method is often applied to the exploration fields of metal exploration detail and exploration stages, underground fluid chase in hydrogeological engineering geological investigation and the like. However, the conventional charging principle is based on the precondition of stabilizing the equipotential body of an ideal conductor (resistivity of zero) in an electric current field, such as a metal ore body or hypersalinity groundwater, which has a low resistivity relative to surrounding rocks, and can be regarded as an ideal conductor approximately, and also the non-ideal conductor (non-equipotential body or lossy medium) is generally regarded as an ideal conductor. However, the resistivity of the charging medium in practical application is usually not zero, that is, the practical application relates to a lossy medium target body, and a method for effectively simulating the numerical value of a finite element by a lossy medium contact induced polarization method is not available at present.
In practical exploration application, if the scale of the charging body is smaller or the central burial depth of the charging body is larger, the charging electric field of the charging body is similar to the point power supply electric field at the center, so that the central burial depth can be deduced by using the point power supply electric field potential curve distribution. When the charging body is large in scale or the center of the charging body is shallow in burial depth, the distortion of the charging electric field observed on the ground surface is obviously different from the electric field distribution of the power supply of the underground point. In general, in the conventional method, a non-ideal conductor (non-equipotential body) is determined based on the non-overlapping relationship between the maximum value of the potential curve and the projected position of the charging point on the ground. And the range inference of the allelotype is to propose the power supply of different charging points, and comprehensively judge the range of the allelotype by matching with other geophysical prospecting methods. Therefore, the development of the finite element numerical simulation method of the contact induced polarization method for the lossy medium with any shape is very important in practical significance, and the theoretical support of the charge induced polarization method of the lossy medium with a complex shape is provided for practical exploration application.
The technical field of engineering, hydrology and environmental geophysical exploration does not see research reports on finite element numerical simulation of a lossy dielectric contact induced polarization method in any shape, so that development of an innovative finite element numerical simulation method of the lossy dielectric contact induced polarization method in any shape is urgently needed for detecting and researching target body distribution under the complex lossy dielectric power supply condition related to engineering, hydrology and environmental geophysical exploration.
Through the above analysis, the problems and defects existing in the prior art are as follows:
the existing contact induced polarization method is poor in theoretical applicability due to the fact that a complex lossy medium related to actual exploration lacks an effective numerical simulation method.
Disclosure of Invention
Aiming at the problems existing in the prior art, the invention provides a finite element numerical simulation method of a contact induced polarization method.
The invention is realized in such a way that a method for simulating the numerical value of a finite element by a contact induced polarization method comprises the following steps:
step one, partitioning a lossy dielectric region with any shape: splitting a tetrahedral staggered grid aiming at a three-dimensional lossy dielectric region with any shape, and dispersing a contact type charging three-dimensional region into a plurality of tetrahedral grids;
step two, solving a finite element method numerical equation: establishing an equation set which meets the potential of the contact induced polarization method, and solving the equation set to obtain the total potential distribution of the subdivision grid;
step three, calculating contact type polarizability: converting the corresponding polarization value by using the equivalent resistivity of each subdivision grid, and solving and calculating by using a contact induced polarization finite element method of all subdivision grid units again to obtain the total potential distribution under the condition of considering contact induced polarization;
step four, calculating the visual polarization rate by a contact induced polarization method: and calculating the visual polarization value of each ground measuring point according to the visual polarization value calculation formula.
Further, in the first step, the arbitrarily shaped lossy dielectric region division specifically includes:
a Cartesian rectangular coordinate system is adopted, the Z axis of the coordinate system is taken to vertically upwards, and the X axis and Y axis directions are determined according to a right-hand spiral rule; determining a regular space distribution range of a three-dimensional lossy medium and surrounding rock medium calculation area containing any shape in the direction of a coordinate axis X, Y, Z, dividing the regular calculation area into M grid nodes according to a certain interval in the X direction by adopting regular cuboid grids, dividing the regular calculation area into N grid nodes according to a certain interval in the Y direction, dividing the regular calculation area into Q grid nodes according to a certain interval in the Z direction, and dividing each regular cuboid grid into staggered tetrahedral grids.
Further, the node lengths of the three-axis direction subdivision grid of the coordinate system are allowed to be different from each other, and the interval subdivision length d is required to satisfy d less than or equal to L/10, wherein L is the spreading length of the three-dimensional lossy medium with any shape in the directions of all coordinate axes, and thus, the calculation area containing the three-dimensional lossy medium with any shape and the surrounding rock medium is subdivided into (M-1) x (N-1) x (Q-1) x 5 grid units;
according to the spatial distribution of the three-dimensional lossy dielectric with arbitrary shape in the coordinate system, nm grid cells which correspond to the spatial distribution of the three-dimensional lossy dielectric with arbitrary shape in (M-1) x (N-1) x (Q-1) x 5 grid cells are determined, and Nm three-dimensional lossy dielectric subdivision grid cells with arbitrary shape correspond to Nm cell point power supplies.
Further, the node coordinate of the three-dimensional lossy medium with any shape distributed in the X direction of the coordinate axis is X i The coordinates of the distributed nodes in the Y direction are Y j The coordinate of the distributed nodes in the Z direction is Z k Wherein i=1, …, mx; j=1, …, ny; z=1, …, qz, mx, ny and Qz are X, Y and the number of Z-direction distributed nodes, respectively, and the arbitrary-shape three-dimensional lossy dielectric subdivision cell grid center coordinates are (x i +d x,i /2,y j +d y,j /2,z k +d z,k /2)，d x,i 、d y,j D z,k The coordinate axis X, Y and the ith, j and k mesh intervals in the Z direction are divided;
the charging point position S coordinates are (x s ,y s ,z s ) The power supply current intensity of the power supply point is I total The subdivision resistivity of the three-dimensional lossy dielectric region with any shape can be allowed to change along with the change of the position of the subdivision cell grid, which is marked as rho q ，q＝1,…,Nm；ρ b Is the resistivity of surrounding rock medium; pi= 3.1415926; the subdivision unit grid is regarded as a point power supply, and the current intensity is recorded as I q Q=1, …, nm, the ground observation point is a coordinate (x a ,y a ,z a ) A=1, …, a; a is the number of observation points; u is the total potential value caused by the condition of the contact type power supply of the mesh nodes of the subdivision unit, namely the total potential value of all the mesh nodes of the calculation area of the three-dimensional lossy medium with any shape to be solved.
Further, in the second step, the finite element method numerical equation solving specifically includes:
calculating the power supply current intensity I according to the formula (1) total During charging, the grid point power supply current intensity coefficient I of three-dimensional lossy dielectric area splitting unit with arbitrary shape q ，q＝1,…,Nm：
Wherein r is s,q The distance from the charging point to the power center of the grid point of the q-th split unit is set;
wherein omega m For region Ω versus current source I m Solid angle of point opening Γ s And Γ ∞ The ground boundary and the underground infinity boundary of the region omega are respectively, n is the normal vector direction of the boundary, delta is the dirac function, and r 1 And r 2 To calculate the spatial distance of the point to the supply point a/B,
is->
Vector sum->
A directional cosine function of the vector;
the linear interpolation, the unit integration and the boundary integration by adopting the finite element method are adopted to finally synthesize the total and then spread the total into a matrix formed by all nodes, and then all the units are added to obtain
Where u is a column vector of u of all nodes; k (K) e ＝1/ρ(K 1e +K 2e )，
Is K e Is used for the expansion matrix of the (c),
p＝(0…I 1 …I Nm …0) T only the potential corresponding to the node where the charging point source is located in p, and the rest are zero.
Let the variation of equation (3) be 0, to obtain a system of linear equations:
Ku＝p (4)
and solving the equation set to obtain the potential of each node.
Further, in the third step, the calculating of the contact polarizability specifically includes:
converting the excitation effect of the considered polarization rate model into resistivity values of all the subdivision units according to an equivalent resistivity formula (5) according to the total potential of the grid points of the three-dimensional lossy dielectric area subdivision units with the arbitrary shape calculated in the second step:
wherein η is the model polarizability value;
further, in the fourth step, the calculating of the visual polarizability by the contact induced polarization method specifically includes:
obtaining the total potential u of grid points of all subdivision units considering the polarizability model according to the third step 2 Obtaining a visual polarization value of the ground observation by a conversion contact induced polarization method according to a limit polarization ratio by a formula (6):
wherein Deltau 2 And delta u is the total potential, eta, of two adjacent grid nodes on the ground, which are solved by a finite element method aiming at a model with the polarizability being considered and a model without the polarizability being considered s And the visual polarization value is corresponding to the intermediate point of the two adjacent grid nodes.
Another object of the present invention is to provide a contact induced polarization finite element numerical simulation system, including:
the lossy dielectric region subdivision module is used for subdividing a three-dimensional lossy dielectric region with any shape by adopting tetrahedral staggered grids, and dispersing the contact type charging three-dimensional region into a plurality of tetrahedral grids;
the finite element method numerical equation solving module is used for establishing an equation set for meeting the potential of the contact induced polarization method, and solving the equation set to obtain the total potential distribution of the subdivision grid;
the contact type polarizability calculation module is used for converting the corresponding polarizability value with the equivalent resistivity of each split grid, solving and calculating the contact type induced polarization finite element method of all the split grid units again, and obtaining the total potential distribution under the condition of considering the contact type induced polarization;
the contact induced polarization method visual polarization rate calculation module is used for calculating the visual polarization rate value of each ground measuring point according to a visual polarization rate calculation formula.
Another object of the present invention is to provide a computer readable storage medium storing a computer program which, when executed by a processor, causes the processor to perform the steps of:
step one, partitioning a lossy dielectric region with any shape: splitting a tetrahedral staggered grid aiming at a three-dimensional lossy dielectric region with any shape, and dispersing a contact type charging three-dimensional region into a plurality of tetrahedral grids;
step two, solving a finite element method numerical equation: establishing an equation set which meets the potential of the contact induced polarization method, and solving the equation set to obtain the total potential distribution of the subdivision grid;
step three, calculating contact type polarizability: converting the corresponding polarization value by using the equivalent resistivity of each subdivision grid, and solving and calculating by using a contact induced polarization finite element method of all subdivision grid units again to obtain the total potential distribution under the condition of considering contact induced polarization;
step four, calculating the visual polarization rate by a contact induced polarization method: and calculating the visual polarization value of each ground measuring point according to the visual polarization value calculation formula.
Another object of the present invention is to provide an information data processing terminal including a memory and a processor, the memory storing a computer program which, when executed by the processor, causes the processor to execute the touch induced polarization finite element numerical simulation method.
By combining all the technical schemes, the invention has the advantages and positive effects that:
compared with the prior art, the numerical simulation method for developing the contact induced polarization method aiming at the complex lossy medium in the geophysical exploration field is one of the invention points of the patent, and provides an effective numerical simulation method for the geophysical exploration field of the complex lossy medium power supply.
The invention is a novel numerical simulation algorithm, and according to the limit value problem of the finite element method of the lossy dielectric subdivision grid, which is designed in the step two and meets the special charging current intensity, the feasibility and the precision of the numerical simulation of the complex lossy dielectric development contact induced polarization method in the geophysical exploration field are unified.
The invention solves the difficult problem of the traditional finite element numerical simulation calculation method of the lossy dielectric contact induced polarization method with any shape, has simple algorithm structure, simple and convenient realization process, reasonable design and stable and reliable calculation result, can be widely applied to the numerical simulation calculation of the exploration polarizability of the lossy dielectric with any shape related to the engineering, hydrology and environmental geophysical exploration field, and provides method support for the engineering, hydrology and environmental geophysical exploration numerical simulation, inversion and explanation based on complex lossy dielectric power supply.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present application, the following description will briefly explain the drawings needed in the embodiments of the present application, and it is obvious that the drawings described below are only some embodiments of the present application, and that other drawings can be obtained according to these drawings without inventive effort for a person skilled in the art.
Fig. 1 is a schematic diagram of a calculation flow of a finite element numerical simulation method of a contact induced polarization method according to an embodiment of the present invention.
FIG. 2 is a schematic three-dimensional view of a design model provided by an embodiment of the present invention.
FIG. 3 is a three-dimensional grid cut-away view of a design model provided by an embodiment of the present invention.
FIG. 4 is a schematic diagram of a design model and an observation system according to an embodiment of the present invention.
Fig. 5 is a graph showing the comparison between the Numerical simulation result (numerical_new_group) and the Physical simulation result (physical_new_group) of the design model contact induced polarization method and the conventional induced polarization method (physical_old_group) according to the embodiment of the present invention.
Detailed Description
The present invention will be described in further detail with reference to the following examples in order to make the objects, technical solutions and advantages of the present invention more apparent. It should be understood that the specific embodiments described herein are for purposes of illustration only and are not intended to limit the scope of the invention.
Aiming at the problems existing in the prior art, the invention provides a finite element numerical simulation method of a contact induced polarization method, and the invention is described in detail below with reference to the accompanying drawings.
As shown in fig. 1, the method for simulating the numerical value of the finite element by the contact induced polarization method provided by the embodiment of the invention comprises the following steps:
step one: arbitrarily shaped lossy dielectric region subdivision:
referring to fig. 2, the cartesian rectangular coordinate system is adopted, and the directions of the X axis and the Y axis are determined according to the right-hand screw rule by taking the coordinate system with the Z axis vertically upward. Referring to fig. 3, it is determined that a regular spatial distribution range of a three-dimensional lossy medium and a surrounding rock medium calculation region including any shape in the direction of a coordinate axis X, Y, Z is 0 to 100cm, 0 to 120cm, 0 to 30cm, the regular calculation region is divided into m=58 mesh nodes at a certain pitch in the X direction by using regular rectangular parallelepiped meshes, n=58 mesh nodes are divided into n=58 mesh nodes at a certain pitch in the Y direction, q=17 mesh nodes are divided into a certain pitch in the Z direction, and then staggered tetrahedral mesh division is performed for each regular rectangular parallelepiped mesh. The node lengths of the three-axis direction subdivision grids of the coordinate system are allowed to be different from each other, and the interval subdivision length d is required to be less than or equal to d and less than or equal to L/10, wherein L is the spreading length of the three-dimensional lossy medium with any shape in the directions of all coordinate axes, and thus, the calculation area containing the three-dimensional lossy medium with any shape and the surrounding rock medium is subdivided into (M-1) x (N-1) x (Q-1) x 5 grid units.
According to the spatial distribution of the three-dimensional lossy medium with any shape in the coordinate system, nm= 57188 grid cells corresponding to the spatial distribution of the three-dimensional lossy medium in (M-1) x (N-1) x (Q-1) x 5 grid cells are determined, and Nm three-dimensional lossy medium with any shape is split into grid cells corresponding to Nm cell point power supplies.
The node coordinate of the distribution node of the three-dimensional lossy medium with any shape in the X direction of the coordinate axis is X i The coordinates of the distributed nodes in the Y direction are Y j The coordinate of the distributed nodes in the Z direction is Z k Wherein i=1, …, mx; j=1, …, ny; z=1, …, qz, mx, ny and Qz are X, Y and the number of Z-direction distribution nodes, respectively. The center coordinate of the grid of the three-dimensional lossy dielectric subdivision unit with any shape is (x) i +d x,i /2,y j +d y,j /2,z k +d z,k /2)，d x,i 、d y,j D z,k The mesh spacing is divided for coordinate axis X, Y and Z-direction i, j and k.
The charging point position S coordinates are (x s ,y s ,z s ) The power supply current intensity of the power supply point is I total With reference to fig. 4, arbitrarily shaped three-dimensional lossy dielectric region subdivision resistivity may be allowed to vary with subdivision cell grid position, denoted as ρ =1 amp q ＝15Ω.m，q＝1,…,Nm；ρ b =50Ω.m is the surrounding rock medium resistivity; pi= 3.1415926; the subdivision unit grid is regarded as a point power supply, and the current intensity is recorded as I q Q=1, …, nm, the ground observation point is a coordinate (x a ,y a ,z a ) A=1, …, a; a=30 is the number of observation points; u is the total potential value caused by the condition of the contact type power supply of the mesh nodes of the subdivision unit, namely, step three is to solve all three-dimensional lossy media with arbitrary shapes in the calculation areaTotal potential value of the grid node.
Step two: finite element method numerical equation solving:
calculating the power supply current intensity I according to the formula (1) total During charging, the grid point power supply current intensity coefficient I of three-dimensional lossy dielectric area splitting unit with arbitrary shape q ，q＝1,…,Nm：
Wherein r is s,q The distance from the charging point to the power center of the grid point of the q-th split unit is set.
Wherein omega m For region Ω versus current source I m Solid angle of point opening Γ s And Γ ∞ The ground boundary and the underground infinity boundary of the region omega are respectively, n is the normal vector direction of the boundary, delta is the dirac function, and r 1 And r 2 To calculate the spatial distance of the point to the supply point a/B,
is->
Vector sum->
The directional cosine function of the vector.
The linear interpolation, the unit integration and the boundary integration by adopting the finite element method are adopted to finally synthesize the total and then spread the total into a matrix formed by all nodes, and then all the units are added to obtain
Where u is a column vector of u of all nodes; k (K) e ＝1/ρ(K 1e +K 2e )，
Is K e Is used for the expansion matrix of the (c),
p＝(0…I 1 …I Nm …0) T only the potential corresponding to the node where the charging point source is located in p, and the rest are zero.
Let the variation of equation (3) be 0, to obtain a system of linear equations:
Ku＝p (4)
and solving the equation set to obtain the potential of each node.
Step three: contact type polarizability calculation
Converting the excitation effect of the considered polarization rate model into resistivity values of all the subdivision units according to an equivalent resistivity formula (5) according to the total potential of the grid points of the three-dimensional lossy dielectric area subdivision units with the arbitrary shape calculated in the second step:
wherein, eta is the model polarization value of 6.8 percent and the background medium polarization value of 1 percent. Solving according to formula (5) to obtain total potential u of grid points of all split units considering the polarizability model 2 Obtaining a visual polarization value of the ground observation by a conversion contact induced polarization method according to a limit polarization ratio by a formula (6):
wherein Deltau 2 And delta u is the total potential, eta, of two adjacent grid nodes on the ground, which are solved by a finite element method aiming at a model with the polarizability being considered and a model without the polarizability being considered s Corresponding to the intermediate point of the two adjacent grid nodesVisual polarization value.
The invention also provides a contact induced polarization method finite element numerical simulation system, which comprises:
the lossy dielectric region subdivision module is used for subdividing a three-dimensional lossy dielectric region with any shape by adopting tetrahedral staggered grids, and dispersing the contact type charging three-dimensional region into a plurality of tetrahedral grids;
the finite element method numerical equation solving module is used for establishing an equation set for meeting the potential of the contact induced polarization method, and solving the equation set to obtain the total potential distribution of the subdivision grid;
the contact type polarizability calculation module is used for converting the corresponding polarizability value with the equivalent resistivity of each split grid, solving and calculating the contact type induced polarization finite element method of all the split grid units again, and obtaining the total potential distribution under the condition of considering the contact type induced polarization;
the contact induced polarization method visual polarization rate calculation module is used for calculating the visual polarization rate value of each ground measuring point according to a visual polarization rate calculation formula.
The visual polarization rate curve calculated by numerical simulation aiming at the soil tank lossy dielectric contact induced polarization method and the physical simulation observation result of the same model (fig. 2, 3 and 4) are adopted. Referring to fig. 5, a graph of the Numerical simulation result (numerical_new_group) and the physical_new_group (physical_group) of the conventional induced polarization method (physical_old_group) of the present invention is shown. In the figure, cm represents cm, and Chargeability represents polarizability. Compared with the prior art, the method has the advantages that the calculation result is better matched with the physical simulation observation result from the aspect of morphology, and only weak difference exists near the boundary of the model, which is caused by the fact that the physical simulation environment and the numerical simulation environment cannot be completely consistent. The fitting relative error of the two data is used for measuring the degree of coincidence of the calculated results, the fitting relative error of the potential curve shown in fig. 5 is 2%, and meanwhile, compared with the traditional induced polarization method, the contact induced polarization method observation result is much higher in abnormal amplitude, so that the accuracy and feasibility of the embodiment of the invention are verified.
The foregoing is merely illustrative of specific embodiments of the present invention, and the scope of the invention is not limited thereto, but any modifications, equivalents, improvements and alternatives falling within the spirit and principles of the present invention will be apparent to those skilled in the art within the scope of the present invention.
Claims (8)
1. The method for simulating the numerical value of the finite element by the contact induced polarization method is characterized by comprising the following steps of:
step one, partitioning a lossy dielectric region with any shape: splitting a tetrahedral staggered grid aiming at a three-dimensional lossy dielectric region with any shape, and dispersing a contact type charging three-dimensional region into a plurality of tetrahedral grids;
step two, solving a finite element method numerical equation: establishing an equation set which meets the potential of the contact induced polarization method, and solving the equation set to obtain the total potential distribution of the subdivision grid;
step three, calculating contact type polarizability: converting the corresponding polarization value by using the equivalent resistivity of each subdivision grid, and solving and calculating by using a contact induced polarization finite element method of all subdivision grid units again to obtain the total potential distribution under the condition of considering contact induced polarization;
step four, calculating the visual polarization rate by a contact induced polarization method: calculating the visual polarization value of each ground measuring point according to the visual polarization value calculation formula;
in the second step, the finite element method numerical equation solving specifically includes:
calculating power supply current intensity coefficients Iq, q=1, …, nm of grid points of the three-dimensional lossy-dielectric-area splitting unit of arbitrary shape at the time of charging of the power supply current intensity Itotal according to the formula (1):
wherein r is s,q The distance from the charging point to the power center of the grid point of the q-th split unit is set;
wherein ωm is the region Ω versus the current source I m Solid angle of point opening Γ s And Γ ∞ The ground boundary and the underground infinity boundary of the region omega are respectively, n is the normal vector direction of the boundary, delta is the dirac function, and r 1 And r 2 To calculate the spatial distance of the point to the power supply point A, B,
is->
Vector sum->
A directional cosine function of the vector;
the linear interpolation, the unit integration and the boundary integration by adopting the finite element method are adopted to finally synthesize the total and then spread the total into a matrix formed by all nodes, and then all the units are added to obtain
Where u is a column vector of u of all nodes; k (K) e ＝1/ρ(K 1e +K 2e )，
Is K e Is used for the expansion matrix of the (c),
p＝(0…I 1 …I Nm …0) T only the potential corresponding to the node where the charging point source is located in p, and the rest is zero;
let the variation of equation (3) be 0, to obtain a system of linear equations:
Ku＝p (4)
and solving the equation set to obtain the potential of each node.
2. The method for simulating the numerical value of a finite element by a contact induced polarization method according to claim 1, wherein in the first step, the arbitrarily shaped lossy dielectric region division specifically comprises:
a Cartesian rectangular coordinate system is adopted, the Z axis of the coordinate system is taken to vertically upwards, and the X axis and Y axis directions are determined according to a right-hand spiral rule; determining a regular space distribution range of a three-dimensional lossy medium and surrounding rock medium calculation area containing any shape in the direction of a coordinate axis X, Y, Z, dividing the regular calculation area into M grid nodes according to a certain interval in the X direction by adopting regular cuboid grids, dividing the regular calculation area into N grid nodes according to a certain interval in the Y direction, dividing the regular calculation area into Q grid nodes according to a certain interval in the Z direction, and dividing each regular cuboid grid into staggered tetrahedral grids.
3. The method for simulating numerical values of finite elements by using a contact induced polarization method according to claim 2, wherein the lengths of nodes of the split grids in the three-axis direction of the coordinate system are allowed to be different from each other, and the split length d between the nodes is required to satisfy d.ltoreq.L/10, wherein L is the spread length of the three-dimensional lossy medium in each coordinate axis direction in any shape, and thus, the calculation area containing the three-dimensional lossy medium in any shape and the surrounding rock medium is split into (M-1) x (N-1) x (Q-1) x 5 grid elements;
according to the spatial distribution of the three-dimensional lossy dielectric with arbitrary shape in the coordinate system, nm grid cells which correspond to the spatial distribution of the three-dimensional lossy dielectric with arbitrary shape in (M-1) x (N-1) x (Q-1) x 5 grid cells are determined, and Nm three-dimensional lossy dielectric subdivision grid cells with arbitrary shape correspond to Nm cell point power supplies.
4. A joint as claimed in claim 3The finite element numerical simulation method of the touch induced polarization method is characterized in that the distribution node coordinate of the three-dimensional lossy medium with any shape in the X direction of a coordinate axis is xi, the distribution node coordinate in the Y direction is yj, and the distribution node coordinate in the Z direction is zk, wherein i=1, … and Mx; j=1, …, ny; z=1, …, qz, mx, ny and Qz are X, Y and the number of Z-direction distributed nodes, respectively, and the arbitrary-shape three-dimensional lossy dielectric subdivision cell grid center coordinates are (x i +d x ,i/2,y j +d y ,j/2,z k +d z ,k/2)，d x,i 、d y,j D z,k The coordinate axis X, Y and the ith, j and k mesh intervals in the Z direction are divided;
the charging point position S coordinates are (x s ,y s ,z s ) The power supply current intensity of the power supply point is Itotal, the subdivision resistivity of the three-dimensional lossy dielectric area with any shape can be allowed to change along with the change of the mesh position of the subdivision unit, and the subdivision resistivity is marked as ρq, q=1, … and Nm; ρb is the resistivity of the surrounding rock medium; pi= 3.1415926; the split unit grid is regarded as a point power supply, the current intensity is recorded as Iq, q=1, …, nm, and the ground observation point is a coordinate (x a ,y a ,z a ) A=1, …, a; a is the number of observation points; u is the total potential value caused by the condition of the contact type power supply of the mesh nodes of the subdivision unit, namely the total potential value of all the mesh nodes of the calculation area of the three-dimensional lossy medium with any shape to be solved.
5. The method for simulating the numerical value of a finite element by a contact induced polarization method according to claim 1, wherein in the third step, the calculation of the contact polarization ratio specifically includes:
converting the excitation effect of the considered polarization rate model into resistivity values of all the subdivision units according to an equivalent resistivity formula (5) according to the total potential of the grid points of the three-dimensional lossy dielectric area subdivision units with the arbitrary shape calculated in the second step:
wherein η is the model polarizability value.
6. The method for simulating the numerical value of a finite element by using a contact induced polarization method according to claim 1, wherein the calculating of the apparent polarization ratio by using the contact induced polarization method specifically comprises:
obtaining the total potential u2 of grid points of all subdivision units considering the polarizability model according to the third step, and obtaining according to the limit polarizability
Taking a formula (6) to convert the visual polarization value of the ground observation by the contact induced polarization method:
wherein Deltau 2 And delta u is the total potential solved by the finite element method of two adjacent grid nodes on the ground aiming at the mode of considering the polarization rate and the mode of not considering the polarization rate, and eta is the visual polarization rate value of the middle point of the two adjacent grid nodes.
7. An information data processing terminal, characterized in that the information data processing terminal comprises a memory and a processor, the memory storing a computer program which, when executed by the processor, causes the processor to execute the touch induced polarization finite element numerical simulation method according to any one of claims 1 to 6.
8. A contact induced polarization finite element numerical simulation system for implementing the simulation method of claim 1, wherein the contact induced polarization finite element numerical simulation system comprises:
the lossy dielectric region subdivision module is used for subdividing a three-dimensional lossy dielectric region with any shape by adopting tetrahedral staggered grids, and dispersing the contact type charging three-dimensional region into a plurality of tetrahedral grids;
the finite element method numerical equation solving module is used for establishing an equation set for meeting the potential of the contact induced polarization method, and solving the equation set to obtain the total potential distribution of the subdivision grid;
the contact type polarizability calculation module is used for converting the corresponding polarizability value with the equivalent resistivity of each split grid, solving and calculating the contact type induced polarization finite element method of all the split grid units again, and obtaining the total potential distribution under the condition of considering the contact type induced polarization;
the contact induced polarization method visual polarization rate calculation module is used for calculating the visual polarization rate value of each ground measuring point according to a visual polarization rate calculation formula.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110497481.1A CN113204905B (en) | 2021-05-08 | 2021-05-08 | Numerical simulation method for finite element by contact induced polarization method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110497481.1A CN113204905B (en) | 2021-05-08 | 2021-05-08 | Numerical simulation method for finite element by contact induced polarization method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113204905A CN113204905A (en) | 2021-08-03 |
| CN113204905B true CN113204905B (en) | 2023-06-20 |
Family
ID=77029736
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110497481.1A Active CN113204905B (en) | 2021-05-08 | 2021-05-08 | Numerical simulation method for finite element by contact induced polarization method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113204905B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113624634B (en) * | 2021-08-11 | 2022-04-22 | 北京师范大学 | Method for estimating content of metal elements in buried environment |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105510982A (en) * | 2015-12-16 | 2016-04-20 | 武汉工程大学 | TBM construction tunnel focusing type forward three-dimensional multi-electrode online detection system based on induced polarization method |
| CN108287371A (en) * | 2018-01-31 | 2018-07-17 | 中南大学 | Background grid Adaptive meshing method in dc resistivity element-free menthod |
Family Cites Families (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4781797B2 (en) * | 2005-11-29 | 2011-09-28 | 株式会社 日立ディスプレイズ | Liquid crystal display |
| CN107742015B (en) * | 2017-09-30 | 2021-04-23 | 中南大学 | DC induced polarization method three-dimensional numerical simulation method based on arbitrary dipole-dipole device |
| CN109470145A (en) * | 2018-12-07 | 2019-03-15 | 哈尔滨工业大学 | Polarization Modulation high resolution Stereo Vision Measurement System and method |
| CN209460416U (en) * | 2019-04-22 | 2019-10-01 | 陕西地矿第二综合物探大队有限公司 | Well-ground induced polarization method measuring device |
| CN110068867B (en) * | 2019-05-08 | 2020-12-25 | 桂林理工大学 | Method for monitoring heavy metal sewage leakage by induced polarization method with embedded measuring electrode |
| CN110007351B (en) * | 2019-05-08 | 2024-06-21 | 桂林理工大学 | Induced polarization method for detecting heavy metal sewage |
| CN110068873B (en) * | 2019-05-10 | 2020-09-25 | 成都理工大学 | Magnetotelluric three-dimensional forward modeling method based on spherical coordinate system |
| CN111755945B (en) * | 2020-06-03 | 2021-06-01 | 郑州大学 | Tunable dual-wavelength plasma nano laser and optical material thereof |
| CN111767669A (en) * | 2020-07-08 | 2020-10-13 | 湖南省有色地质勘查研究院 | Novel pseudo-random induced polarization finite element numerical simulation method and system |
-
2021
- 2021-05-08 CN CN202110497481.1A patent/CN113204905B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105510982A (en) * | 2015-12-16 | 2016-04-20 | 武汉工程大学 | TBM construction tunnel focusing type forward three-dimensional multi-electrode online detection system based on induced polarization method |
| CN108287371A (en) * | 2018-01-31 | 2018-07-17 | 中南大学 | Background grid Adaptive meshing method in dc resistivity element-free menthod |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113204905A (en) | 2021-08-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Rücker et al. | Three-dimensional modelling and inversion of dc resistivity data incorporating topography—I. Modelling | |
| CN107742015B (en) | DC induced polarization method three-dimensional numerical simulation method based on arbitrary dipole-dipole device | |
| CN108388707B (en) | Direct-current magnetic bias calculation method based on field-circuit coupling under three-dimensional asymmetric structure soil model | |
| CN112949134B (en) | Earth-well transient electromagnetic inversion method based on non-structural finite element method | |
| CN110346835B (en) | Magnetotelluric forward modeling method, forward modeling system, storage medium, and electronic device | |
| CN112287534A (en) | NUFFT-based two-dimensional magnetic anomaly fast forward modeling method and device | |
| CN110346834B (en) | Forward modeling method and system for three-dimensional frequency domain controllable source electromagnetism | |
| CN113255230B (en) | Gravity model forward modeling method and system based on MQ radial basis function | |
| CN105427190B (en) | Calculation method for ground three-dimensional power frequency electric field below UHVAC power transmission line in complex terrain | |
| CN113051779B (en) | Numerical simulation method of three-dimensional direct-current resistivity method | |
| CN114970290B (en) | Ground transient electromagnetic method parallel inversion method and system | |
| CN113204905B (en) | Numerical simulation method for finite element by contact induced polarization method | |
| CN115238550A (en) | Self-adaptive unstructured grid landslide rainfall geoelectric field numerical simulation calculation method | |
| CN111638556A (en) | Magnetotelluric forward modeling method and device based on geospatial solution strategy and storage medium | |
| CN106019394A (en) | Three-dimensional parallel inversion method for nonlinear conjugate gradient of ocean magnetotelluric field | |
| CN117538945A (en) | Three-dimensional magnetotelluric multi-resolution inversion method, device, equipment and medium | |
| CN111103628B (en) | Two-dimensional inversion method and device for electric field data in magnetotelluric (TE) polarization mode | |
| CN112163332A (en) | 2.5-dimensional natural electric field numerical simulation method based on infinite element coupling of natural units | |
| CN108873084B (en) | It is a kind of based on Partition of Unity integral dc resistivity without unit forward modeling method | |
| CN114970289B (en) | Three-dimensional magnetotelluric anisotropy forward modeling numerical simulation method, equipment and medium | |
| CN114047554B (en) | Earth resistivity model modeling method, apparatus, computer device and storage medium | |
| Cheruvu et al. | A spectral finite volume transport scheme on the cubed-sphere | |
| CN113238285B (en) | Resistivity calculation method, system and terminal for geophysical charging method exploration | |
| CN113779818B (en) | Three-dimensional geologic body electromagnetic field numerical simulation method, device, equipment and medium thereof | |
| Gajda-Zagórska et al. | Multi-objective hierarchic memetic solver for inverse parametric problems |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |