CN113204905A - Contact induced polarization finite element numerical simulation method - Google Patents
Contact induced polarization finite element numerical simulation method Download PDFInfo
- Publication number
- CN113204905A CN113204905A CN202110497481.1A CN202110497481A CN113204905A CN 113204905 A CN113204905 A CN 113204905A CN 202110497481 A CN202110497481 A CN 202110497481A CN 113204905 A CN113204905 A CN 113204905A
- Authority
- CN
- China
- Prior art keywords
- grid
- subdivision
- induced polarization
- finite element
- contact
- 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.)
- Granted
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)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Remote Sensing (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 contact induced polarization method finite element numerical simulation method. The method solves the problem of the traditional finite element numerical simulation calculation method of the lossy dielectric contact induced polarization method in 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 lossy dielectric charging method exploration polarizability in any shape in the field of engineering, hydrology and environmental geophysical exploration, and provides method support for the numerical simulation, inversion and explanation of the engineering, hydrology and environmental geophysical exploration based on a 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 contact induced polarization finite element numerical simulation method.
Background
At present, the charging method is to directly connect a power supply electrode (generally a positive electrode) to a natural or artificially exposed good conductor outcrop and an underground water outcrop, and to place another power supply electrode at a position satisfying the requirement of 'infinity', and to estimate the physical properties and spatial distribution of the good conductor by observing the change and distribution rule of a charging electric field through two measuring electrodes. Because the charging medium range and the observation potential abnormal range are similar in shape, the actual charging medium range can be inferred according to the observation potential abnormal, and the method is often applied to the exploration fields of metal exploration in detail investigation and exploration stages, underground fluid exploration in hydrogeological engineering geological survey and the like. However, the conventional charging method principle is based on the premise that the equivalent volume of an ideal conductor (with the resistivity being zero) in a stable current field, such as a metal ore body or high-salinity underground water, has a very low resistivity relative to surrounding rocks, and can be approximately regarded as an ideal conductor, and a non-ideal conductor (non-equivalent volume or lossy medium) is generally also approximately 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 finite element numerical value by a lossy medium contact induced polarization method is not available at present.
In practical exploration application, if the charging body is small in scale or the center of the charging body is large in buried depth, the charging electric field is close to the point power supply electric field in the center, and therefore the center buried depth can be inferred by using the potential curve distribution of the point power supply electric field. 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 earth surface is obviously different from the distribution of the power supply electric field of the underground point. In general, in the conventional method for solving the above problem, a non-ideal conductor (non-equivalent body) is determined based on a relationship that the maximum value of the potential curve and the projected position of the charging point on the ground do not coincide with each other. And for the range inference of unequal bits, the power supply at different charging point positions is suggested, and the range of the unequal bits is comprehensively judged by matching with other geophysical prospecting methods. Therefore, the method for effectively simulating the finite element numerical value of the contact induced polarization method aiming at the lossy medium with any shape is developed, the theoretical support of the charging induced polarization method of the lossy medium with a complex shape is provided for practical exploration and application, and the method has very important practical significance.
The technical field of current engineering, hydrological and environmental geophysical exploration does not have any research report for numerical simulation of finite elements of lossy medium contact induced polarization methods in any shapes, so that an innovative numerical simulation method for finite elements of lossy medium contact induced polarization methods in any shapes is urgently needed to be developed and used for detecting and researching target body distribution under the condition of complex lossy medium power supplies related to engineering, hydrological and environmental geophysical exploration.
Through the above analysis, the problems and defects of the prior art are as follows:
the complex lossy medium related to the actual exploration by the existing contact induced polarization method lacks an effective numerical simulation method, so that the theoretical applicability of the method is poor.
Disclosure of Invention
Aiming at the problems in the prior art, the invention provides a finite element numerical simulation method by a contact induced polarization method.
The invention is realized in this way, and a contact induced polarization finite element numerical simulation method comprises the following steps:
step one, dividing a lossy medium region in any shape: aiming at a three-dimensional lossy medium region in any shape, dividing by using a tetrahedral staggered grid, and dispersing a contact type charging three-dimensional region into a plurality of tetrahedral grids;
step two, solving a numerical equation of a finite element method: establishing an equation set which is satisfied by 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 the contact polarizability: converting the corresponding polarizability value by the equivalent resistivity of each subdivision grid, solving and calculating by a unit contact induced polarization finite element method of all subdivision grids again, and acquiring 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 apparent polarization rate value of each ground measuring point according to an apparent polarization rate solving formula.
Further, in the first step, the dividing of the lossy medium region of any shape specifically includes:
a Cartesian rectangular coordinate system is adopted, the Z axis of the coordinate system is taken to be vertical upwards, and the directions of the X axis and the Y axis are determined according to a right-hand screw rule; determining the regular spatial distribution range of a calculation region containing three-dimensional lossy media and surrounding rock media in any shape in the direction of a coordinate axis X, Y, Z, dividing the regular calculation region into M grid nodes in the X direction according to a certain distance by adopting regular cuboid grids, dividing the regular calculation region into N grid nodes in the Y direction according to a certain distance, dividing the regular calculation region into Q grid nodes in the Z direction according to a certain distance, and then carrying out staggered tetrahedral grid division on each regular cuboid grid.
Furthermore, the lengths of the nodes of the subdivision grids in the three-axis direction of the coordinate system are allowed to be different, and the distance subdivision length d of the nodes is equal to or less than L/10, wherein L is the spreading length of the three-dimensional lossy medium in any shape in each coordinate axis direction, so that a calculation area containing the three-dimensional lossy medium in any shape and the surrounding rock medium is divided 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 a coordinate system, Nm grid cells corresponding to the spatial distribution of the three-dimensional lossy medium with any shape in (M-1) x (N-1) x (Q-1) x 5 grid cells are determined, and Nm cell point power supplies corresponding to the Nm grid cells are divided by the Nm three-dimensional lossy medium with any shape.
Furthermore, the coordinate of the distributed node of the three-dimensional lossy medium with any shape in the X direction of the coordinate axis is XiAnd the coordinate of the distributed nodes in the Y direction is YjThe coordinate of the distributed nodes in the Z direction is ZkWherein i is 1, …, Mx; j is 1, …, Ny; z is 1 and …, Qz, Mx, Ny and Qz are respectively X, Y and the number of distributed nodes in the Z direction, and the grid center coordinate of the arbitrarily-shaped three-dimensional lossy medium subdivision unit is (x)i+dx,i/2,yj+dy,j/2,zk+dz,k/2)，dx,i、dy,jAnd dz,kThe coordinate axis X, Y and the i, j and k subdivision grid intervals in the Z direction;
the position S coordinate of the charging point is (x)s,ys,zs) Intensity of supply current to supply pointIs ItotalAnd the subdivision resistivity of the three-dimensional lossy medium region with any shape can be allowed to change along with the change of the grid position of the subdivision unit and is recorded as rhoq，q＝1,…,Nm；ρbIs the surrounding rock medium resistivity; PI is 3.1415926; the subdivision unit grid is taken as a point power supply, and the current intensity is recorded as IqQ is 1, …, Nm, and the ground observation point is the coordinate (x)a,ya,za) A is 1, …, a; a is the number of observation points; and u is a total potential value caused under the condition of the subdivision unit grid node contact type power supply, namely the total potential value of all grid nodes of the three-dimensional lossy medium with any shape to be solved in the calculation area.
Further, in the second step, the solving of the finite element method numerical equation specifically includes:
calculating the current intensity I of the power supply according to the formula (1)totalWhen charging, the grid point power supply current intensity coefficient I of the arbitrary-shape three-dimensional lossy medium area subdivision unitq，q＝1,…,Nm：
Wherein r iss,qThe distance from the charging point to the power center of the mesh point of the q-th subdivision unit is obtained;
wherein ω ismFor region omega to current source ImSolid angle spanned by dots, ΓsAnd Γ∞The ground boundary and the underground infinite boundary of the region omega are respectively, n is the normal vector direction of the boundary, delta is the Dirac function, r1And r2To calculate the spatial distance of the point to the supply point a/B,
is composed of
Vector sum
The directional cosine function of the vector;
linear interpolation, unit integration and boundary integration by using finite element method, finally totally synthesizing and re-expanding into a matrix composed of all nodes, and further adding all elements to obtain
Wherein u is a column vector consisting of u of all nodes; ke＝1/ρ(K1e+K2e)，
Is KeThe spreading matrix of (a) is set,
p＝(0…I1…INm…0)Tand only the potential corresponding to the node where the charging point source is located in p is zero, and the rest are zero.
Let the variation of equation (3) be 0, resulting in a linear equation set:
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:
and according to the total potential of the grid points of the subdivision units of the three-dimensional lossy medium region with any shape obtained by calculation in the step two, converting the excitation effect of the model considering the polarizability into the resistivity value of each subdivision unit according to an equivalent resistivity formula (5):
wherein η is the model polarizability value;
further, in the fourth step, the calculation of the apparent polarizability by the contact induced polarization method specifically includes:
obtaining the total grid point potential u of all subdivision units considering the polarizability model according to the third step2And calculating a formula (6) according to the limit polarizability to convert the apparent polarizability value observed on the ground by the contact induced polarization method:
wherein, Δ u2And delta u is the total potential, eta, of two adjacent grid nodes on the ground, which is solved by a finite element method when the polarizability model is considered and the polarizability model is not consideredsCorresponding to the apparent polarization rate value of 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, which includes:
the lossy medium region subdivision module is used for subdividing the three-dimensional lossy medium region in any shape by adopting tetrahedral staggered grids and dispersing the contact 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 which is satisfied by 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 calculating module is used for converting the corresponding polarizability value by using the equivalent resistivity of each subdivision grid, solving and calculating the unit contact type induced polarization finite element method of all subdivision grids again, and acquiring the total potential distribution under the condition of considering the contact type induced polarization;
and the contact induced polarization method apparent polarization rate calculation module is used for calculating the apparent polarization rate value of each ground measuring point according to an apparent polarization rate solving formula.
It is another object of the present invention 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, dividing a lossy medium region in any shape: aiming at a three-dimensional lossy medium region in any shape, dividing by using a tetrahedral staggered grid, and dispersing a contact type charging three-dimensional region into a plurality of tetrahedral grids;
step two, solving a numerical equation of a finite element method: establishing an equation set which is satisfied by 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 the contact polarizability: converting the corresponding polarizability value by the equivalent resistivity of each subdivision grid, solving and calculating by a unit contact induced polarization finite element method of all subdivision grids again, and acquiring 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 apparent polarization rate value of each ground measuring point according to an apparent polarization rate solving 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, the computer program, when executed by the processor, causing the processor to execute the contact 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 invention discloses a contact induced polarization method numerical simulation method for developing complex lossy media in the field of geophysical exploration, and provides an effective numerical simulation method for the field of complex lossy media power supply geophysical exploration.
The invention is a novel numerical simulation algorithm, and the unification of feasibility and precision of developing the contact induced polarization method numerical simulation aiming at the complex lossy medium in the geophysical exploration field is realized according to the boundary value problem of the lossy medium subdivision grid finite element method which is designed according to the step two and meets the special charging current intensity.
The method solves the problem of the traditional finite element numerical simulation calculation method of the lossy dielectric contact induced polarization method in 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 lossy dielectric charging method exploration polarizability in any shape in the field of engineering, hydrology and environmental geophysical exploration, and provides method support for the numerical simulation, inversion and explanation of the engineering, hydrology and environmental geophysical exploration based on a complex lossy dielectric power supply.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings needed to be used in the embodiments of the present application will be briefly described below, and it is obvious that the drawings described below are only some embodiments of the present application, and it is obvious for those skilled in the art that other drawings can be obtained from the drawings without creative efforts.
Fig. 1 is a schematic calculation flow diagram of a finite element numerical simulation method by a contact induced polarization method according to an embodiment of the present invention.
FIG. 2 is a three-dimensional schematic diagram of a design model provided by an embodiment of the invention.
Fig. 3 is a three-dimensional mesh split view of a design model provided by an embodiment of the 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 comparison graph of the Numerical simulation calculation result (Numerical _ new _ group) of the design model contact induced polarization method provided by the embodiment of the present invention, the Physical simulation result (Physical _ new _ group) and the conventional induced polarization method (Physical _ old _ group).
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described in detail with reference to the following embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
In view of the problems in the prior art, the present invention provides a finite element numerical simulation method by a contact induced polarization method, and the present invention is described in detail below with reference to the accompanying drawings.
As shown in fig. 1, a method for finite element numerical simulation by a contact induced polarization method according to an embodiment of the present invention includes:
the method comprises the following steps: dividing lossy medium regions in any shapes:
referring to fig. 2, a cartesian rectangular coordinate system is adopted, and the Z-axis of the coordinate system is taken to be vertically upward, and then the X-axis and Y-axis directions are determined according to the right-hand screw rule. Referring to fig. 3, determining that the regular spatial distribution range of a calculation region containing a three-dimensional lossy medium and a surrounding rock medium in any shape in the direction of a coordinate axis X, Y, Z is 0-100 cm, 0-120 cm and 0-30 cm, dividing the regular calculation region into M-58 grid nodes in the X direction at a certain interval by adopting regular cuboid grids, dividing the regular calculation region into N-58 grid nodes in the Y direction at a certain interval, dividing the regular calculation region into Q-17 grid nodes in the Z direction at a certain interval, and then performing staggered tetrahedral grid division on each regular cuboid grid. The lengths of the nodes of the division grids in the three-axis direction of the coordinate system are allowed to be different, and the spacing division length d of the division grids is required to meet the condition that d is not more than L/10, wherein L is the distribution length of the three-dimensional lossy medium in any shape in each coordinate axis direction, so that a calculation area containing the three-dimensional lossy medium in any shape and the surrounding rock medium is divided into (M-1) × (N-1) × (Q-1) × 5 grid units.
According to the spatial distribution of the three-dimensional lossy medium with the arbitrary shape in the coordinate system, the Nm (corresponding to the self spatial distribution) of the three-dimensional lossy medium in (M-1) x (N-1) x (Q-1) x 5 grid cells is determined to be 57188 grid cells, and the Nm three-dimensional lossy medium in the arbitrary shape divides the grid cells into Nm cell point power supplies.
The coordinate of a distributed node of the three-dimensional lossy medium with any shape in the X direction of the coordinate axis is XiAnd the coordinate of the distributed nodes in the Y direction is YjThe coordinate of the distributed nodes in the Z direction is ZkWherein i is 1, …, Mx; j is 1, …, Ny; z is 1, …, Qz, Mx, Ny and Qz are X, Y and the number of distributed nodes in the Z direction, respectively. Random shapeThe grid center coordinate of the subdivision unit of the three-dimensional lossy medium is (x)i+dx,i/2,yj+dy,j/2,zk+dz,k/2)，dx,i、dy,jAnd dz,kThe i, j and k grid spacings are plotted along the axis X, Y and the Z direction.
The position S coordinate of the charging point is (x)s,ys,zs) The power supply point has a power supply current intensity of ItotalReferring to fig. 4, the allowable arbitrary-shaped three-dimensional lossy medium region subdivision resistivity varies with the position of the subdivision cell grid, denoted as ρ, at 1 ampereq＝15Ω.m，q＝1,…,Nm；ρbM is the surrounding rock dielectric resistivity; PI is 3.1415926; the subdivision unit grid is taken as a point power supply, and the current intensity is recorded as IqQ is 1, …, Nm, and the ground observation point is the coordinate (x)a,ya,za) A is 1, …, a; a is 30, namely the number of observation points; and u is a total potential value caused under the condition of the subdivision unit grid node contact type power supply, namely the total potential value of all grid nodes of the three-dimensional lossy medium with any shape in the calculation area to be solved in the step three.
Step two: solving numerical equations by a finite element method:
calculating the current intensity I of the power supply according to the formula (1)totalWhen charging, the grid point power supply current intensity coefficient I of the arbitrary-shape three-dimensional lossy medium area subdivision unitq，q＝1,…,Nm：
Wherein r iss,qAnd the distance from the charging point to the power center of the mesh point of the q-th subdivision unit is shown.
Wherein ω ismFor region omega to current source ImSolid angle spanned by dots, ΓsAnd Γ∞Ground of region omega respectivelyBoundary and underground infinite boundary, n is the normal vector direction of the boundary, delta is the Dirac function, r1And r2To calculate the spatial distance of the point to the supply point a/B,
is composed of
Vector sum
The direction cosine function of the vector.
Linear interpolation, unit integration and boundary integration by using finite element method, finally totally synthesizing and re-expanding into a matrix composed of all nodes, and further adding all elements to obtain
Wherein u is a column vector consisting of u of all nodes; ke＝1/ρ(K1e+K2e)，
Is KeThe spreading matrix of (a) is set,
p＝(0…I1…INm…0)Tand only the potential corresponding to the node where the charging point source is located in p is zero, and the rest are zero.
Let the variation of equation (3) be 0, resulting in a linear equation set:
Ku＝p (4)
and solving the equation set to obtain the potential of each node.
Step three: contact polarizability calculation
And according to the total potential of the grid points of the subdivision units of the three-dimensional lossy medium region with any shape obtained by calculation in the step two, converting the excitation effect of the model considering the polarizability into the resistivity value of each subdivision unit according to an equivalent resistivity formula (5):
wherein eta is the model polarizability value of 6.8%, and the background medium polarizability value is 1%. Solving according to the formula (5) to obtain the total grid point potential u of all subdivision units considering the polarizability model2And calculating a formula (6) according to the limit polarizability to convert the apparent polarizability value observed on the ground by the contact induced polarization method:
wherein, Δ u2And delta u is the total potential, eta, of two adjacent grid nodes on the ground, which is solved by a finite element method when the polarizability model is considered and the polarizability model is not consideredsCorresponding to the apparent polarization rate value of the intermediate point of the two adjacent grid nodes.
The invention also provides a contact induced polarization method finite element numerical simulation system, which comprises:
the lossy medium region subdivision module is used for subdividing the three-dimensional lossy medium region in any shape by adopting tetrahedral staggered grids and dispersing the contact 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 which is satisfied by 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 calculating module is used for converting the corresponding polarizability value by using the equivalent resistivity of each subdivision grid, solving and calculating the unit contact type induced polarization finite element method of all subdivision grids again, and acquiring the total potential distribution under the condition of considering the contact type induced polarization;
and the contact induced polarization method apparent polarization rate calculation module is used for calculating the apparent polarization rate value of each ground measuring point according to an apparent polarization rate solving formula.
The apparent polarizability curve calculated by numerical simulation aiming at the lossy medium contact induced polarization method of the soil tank and the physical simulation observation result of the soil tank of the same model (figure 2, figure 3 and figure 4) are used. Referring to fig. 5, a comparison graph of the Numerical simulation calculation result (Numerical _ new _ ground) of the contact-type induced polarization method of the design model of the invention, the Physical simulation result (Physical _ new _ ground) and the conventional induced polarization method (Physical _ old _ ground) is shown. In the figure, cm represents cm, and Chargeability represents polarizability. As can be seen from the comparison, the calculation result of the invention is well matched with the physical simulation observation result in form, and only weak difference exists near the model boundary, 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 coincidence degree of the calculated result, the fitting relative error of the potential curve shown in FIG. 5 is 2%, and meanwhile, the contact induced polarization method observation result is much higher than the abnormal amplitude of the traditional induced polarization method, so that the accuracy and the feasibility of the embodiment of the invention are verified.
The above description is only for the purpose of illustrating the present invention and the appended claims are not to be construed as limiting the scope of the invention, which is intended to cover all modifications, equivalents and improvements that are within the spirit and scope of the invention as defined by the appended claims.
Claims (10)
1. A contact induced polarization finite element numerical simulation method is characterized by comprising the following steps:
step one, dividing a lossy medium region in any shape: aiming at a three-dimensional lossy medium region in any shape, dividing by using a tetrahedral staggered grid, and dispersing a contact type charging three-dimensional region into a plurality of tetrahedral grids;
step two, solving a numerical equation of a finite element method: establishing an equation set which is satisfied by 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 the contact polarizability: converting the corresponding polarizability value by the equivalent resistivity of each subdivision grid, solving and calculating by a unit contact induced polarization finite element method of all subdivision grids again, and acquiring 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 apparent polarization rate value of each ground measuring point according to an apparent polarization rate solving formula.
2. The contact induced polarization finite element numerical simulation method of claim 1, wherein in the first step, the subdivision of the lossy medium region of any shape specifically comprises:
a Cartesian rectangular coordinate system is adopted, the Z axis of the coordinate system is taken to be vertical upwards, and the directions of the X axis and the Y axis are determined according to a right-hand screw rule; determining the regular spatial distribution range of a calculation region containing three-dimensional lossy media and surrounding rock media in any shape in the direction of a coordinate axis X, Y, Z, dividing the regular calculation region into M grid nodes in the X direction according to a certain distance by adopting regular cuboid grids, dividing the regular calculation region into N grid nodes in the Y direction according to a certain distance, dividing the regular calculation region into Q grid nodes in the Z direction according to a certain distance, and then carrying out staggered tetrahedral grid division on each regular cuboid grid.
3. The contact-type induced polarization finite element numerical simulation method of claim 2, wherein the node lengths of the division grids in the three-axis direction of the coordinate system are allowed to be different, and the division length d is equal to or less than L/10, wherein L is the distribution length of the three-dimensional lossy medium with any shape in each coordinate axis direction, so that the calculation area containing the three-dimensional lossy medium with any shape and the surrounding rock medium is divided into (M-1) × (N-1) × (Q-1) × 5 grid elements;
according to the spatial distribution of the three-dimensional lossy medium with any shape in a coordinate system, Nm grid cells corresponding to the spatial distribution of the three-dimensional lossy medium with any shape in (M-1) x (N-1) x (Q-1) x 5 grid cells are determined, and Nm cell point power supplies corresponding to the Nm grid cells are divided by the Nm three-dimensional lossy medium with any shape.
4. The contact induced polarization finite element numerical simulation method of claim 3, wherein the coordinate of the distributed node of the arbitrary-shaped three-dimensional lossy medium in the X direction of the coordinate axis is XiAnd the coordinate of the distributed nodes in the Y direction is YjThe coordinate of the distributed nodes in the Z direction is ZkWherein i is 1, …, Mx; j is 1, …, Ny; z is 1 and …, Qz, Mx, Ny and Qz are respectively X, Y and the number of distributed nodes in the Z direction, and the grid center coordinate of the arbitrarily-shaped three-dimensional lossy medium subdivision unit is (x)i+dx,i/2,yj+dy,j/2,zk+dz,k/2)，dx,i、dy,jAnd dz,kThe coordinate axis X, Y and the i, j and k subdivision grid intervals in the Z direction;
the position S coordinate of the charging point is (x)s,ys,zs) The power supply point has a power supply current intensity of ItotalAnd the subdivision resistivity of the three-dimensional lossy medium region with any shape can be allowed to change along with the change of the grid position of the subdivision unit and is recorded as rhoq，q＝1,…,Nm；ρbIs the surrounding rock medium resistivity; PI is 3.1415926; the subdivision unit grid is taken as a point power supply, and the current intensity is recorded as IqQ is 1, …, Nm, and the ground observation point is the coordinate (x)a,ya,za) A is 1, …, a; a is the number of observation points; and u is a total potential value caused under the condition of the subdivision unit grid node contact type power supply, namely the total potential value of all grid nodes of the three-dimensional lossy medium with any shape to be solved in the calculation area.
5. The contact induced polarization finite element numerical simulation method of claim 1, wherein in the second step, the solution of the finite element numerical equation specifically comprises:
calculating the current intensity I of the power supply according to the formula (1)totalWhen charging, the grid point power supply current of the subdivision unit of the three-dimensional lossy medium region with any shape is strongCoefficient of degree Iq，q＝1,…,Nm：
Wherein r iss,qThe distance from the charging point to the power center of the mesh point of the q-th subdivision unit is obtained;
wherein ω ismFor region omega to current source ImSolid angle spanned by dots, ΓsAnd Γ∞The ground boundary and the underground infinite boundary of the region omega are respectively, n is the normal vector direction of the boundary, delta is the Dirac function, r1And r2To calculate the spatial distance of the point to the power supply point A, B,
is composed of
Vector sum
The directional cosine function of the vector;
linear interpolation, unit integration and boundary integration by using finite element method, finally totally synthesizing and re-expanding into a matrix composed of all nodes, and further adding all elements to obtain
Wherein u is a column vector consisting of u of all nodes; ke＝1/ρ(K1e+K2e)，
Is KeThe spreading matrix of (a) is set,
and only the potential corresponding to the node where the charging point source is located in the p is zero, and the rest of the potentials are zero.
Let the variation of equation (3) be 0, resulting in a linear equation set:
Ku＝p (4)
and solving the equation set to obtain the potential of each node.
6. The contact induced polarization finite element numerical simulation method of claim 1, wherein in step three, the contact polarizability calculation specifically comprises:
and according to the total potential of the grid points of the subdivision units of the three-dimensional lossy medium region with any shape obtained by calculation in the step two, converting the excitation effect of the model considering the polarizability into the resistivity value of each subdivision unit according to an equivalent resistivity formula (5):
wherein η is the model polarizability value.
7. The method for finite element numerical simulation by contact induced polarization method according to claim 1, wherein the step four, the calculation of the apparent polarizability by contact induced polarization method specifically comprises:
obtaining the total grid point potential u of all subdivision units considering the polarizability model according to the third step2And calculating a formula (6) according to the limit polarizability to convert the apparent polarizability value observed on the ground by the contact induced polarization method:
wherein, Δ u2And delta u is the total potential, eta, of two adjacent grid nodes on the ground, which is solved by a finite element method when the polarizability model is considered and the polarizability model is not consideredsCorresponding to the apparent polarization rate value of the intermediate point of the two adjacent grid nodes.
8. A contact induced polarization finite element numerical simulation system, comprising:
the lossy medium region subdivision module is used for subdividing the three-dimensional lossy medium region in any shape by adopting tetrahedral staggered grids and dispersing the contact 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 which is satisfied by 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 calculating module is used for converting the corresponding polarizability value by using the equivalent resistivity of each subdivision grid, solving and calculating the unit contact type induced polarization finite element method of all subdivision grids again, and acquiring the total potential distribution under the condition of considering the contact type induced polarization;
and the contact induced polarization method apparent polarization rate calculation module is used for calculating the apparent polarization rate value of each ground measuring point according to an apparent polarization rate solving formula.
9. 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, dividing a lossy medium region in any shape: aiming at a three-dimensional lossy medium region in any shape, dividing by using a tetrahedral staggered grid, and dispersing a contact type charging three-dimensional region into a plurality of tetrahedral grids;
step two, solving a numerical equation of a finite element method: establishing an equation set which is satisfied by 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 the contact polarizability: converting the corresponding polarizability value by the equivalent resistivity of each subdivision grid, solving and calculating by a unit contact induced polarization finite element method of all subdivision grids again, and acquiring 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 apparent polarization rate value of each ground measuring point according to an apparent polarization rate solving formula.
10. An information data processing terminal, characterized in that the information data processing terminal comprises a memory and a processor, the memory stores a computer program, and the computer program, when executed by the processor, causes the processor to execute the contact induced polarization finite element numerical simulation method according to any one of claims 1 to 7.
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 true CN113204905A (en) | 2021-08-03 |
| CN113204905B 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) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113624634A (en) * | 2021-08-11 | 2021-11-09 | 北京师范大学 | Method for estimating content of metal elements in buried environment |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070121029A1 (en) * | 2005-11-29 | 2007-05-31 | Koichi Fukuda | Liquid crystal display device |
| 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 |
| CN107742015A (en) * | 2017-09-30 | 2018-02-27 | 中南大学 | DC lasering electric method Three-dimensional Numerical Simulation Method based on any dipole dipole device |
| CN108287371A (en) * | 2018-01-31 | 2018-07-17 | 中南大学 | Background grid Adaptive meshing method in dc resistivity element-free menthod |
| CN109470145A (en) * | 2018-12-07 | 2019-03-15 | 哈尔滨工业大学 | Polarization Modulation high resolution Stereo Vision Measurement System and method |
| CN110007351A (en) * | 2019-05-08 | 2019-07-12 | 桂林理工大学 | A kind of induced polarization method detecting heavy metal containing sewage |
| CN110068873A (en) * | 2019-05-10 | 2019-07-30 | 成都理工大学 | A kind of magnetotelluric D integral pin-fin tube method based on spherical coordinate system |
| CN110068867A (en) * | 2019-05-08 | 2019-07-30 | 桂林理工大学 | A kind of induced polarization method monitoring heavy metal containing sewage leakage method of pre-buried measuring electrode |
| CN209460416U (en) * | 2019-04-22 | 2019-10-01 | 陕西地矿第二综合物探大队有限公司 | Well-ground induced polarization method measuring device |
| CN111755945A (en) * | 2020-06-03 | 2020-10-09 | 郑州大学 | 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 (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070121029A1 (en) * | 2005-11-29 | 2007-05-31 | Koichi Fukuda | Liquid crystal display device |
| 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 |
| CN107742015A (en) * | 2017-09-30 | 2018-02-27 | 中南大学 | DC lasering electric method Three-dimensional Numerical Simulation Method based on any dipole dipole device |
| CN108287371A (en) * | 2018-01-31 | 2018-07-17 | 中南大学 | Background grid Adaptive meshing method in dc resistivity element-free menthod |
| 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 |
| CN110007351A (en) * | 2019-05-08 | 2019-07-12 | 桂林理工大学 | A kind of induced polarization method detecting heavy metal containing sewage |
| CN110068867A (en) * | 2019-05-08 | 2019-07-30 | 桂林理工大学 | A kind of induced polarization method monitoring heavy metal containing sewage leakage method of pre-buried measuring electrode |
| CN110068873A (en) * | 2019-05-10 | 2019-07-30 | 成都理工大学 | A kind of magnetotelluric D integral pin-fin tube method based on spherical coordinate system |
| CN111755945A (en) * | 2020-06-03 | 2020-10-09 | 郑州大学 | 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 |
Non-Patent Citations (2)
| Title |
|---|
| 王广涛: "矿井时间域激电场响应特征研究", 《中国优秀硕士学位论文全文数据库 (工程科技Ⅰ辑)》 * |
| 王广涛: "矿井时间域激电场响应特征研究", 《中国优秀硕士学位论文全文数据库 (工程科技Ⅰ辑)》, 15 March 2018 (2018-03-15), pages 021 - 22 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113624634A (en) * | 2021-08-11 | 2021-11-09 | 北京师范大学 | Method for estimating content of metal elements in buried environment |
| CN113624634B (en) * | 2021-08-11 | 2022-04-22 | 北京师范大学 | Method for estimating content of metal elements in buried environment |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113204905B (en) | 2023-06-20 |
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 | |
| Bing et al. | Finite element three dimensional direct current resistivity modelling: accuracy and efficiency considerations | |
| Blome et al. | Advances in three-dimensional geoelectric forward solver techniques | |
| Zhou | Electrical resistivity tomography: a subsurface-imaging technique | |
| Papadopoulos et al. | An algorithm for fast 3D inversion of surface electrical resistivity tomography data: application on imaging buried antiquities | |
| CN112949134B (en) | Earth-well transient electromagnetic inversion method based on non-structural finite element method | |
| CN115238550A (en) | Self-adaptive unstructured grid landslide rainfall geoelectric field numerical simulation calculation method | |
| CN113255230B (en) | Gravity model forward modeling method and system based on MQ radial basis function | |
| CN113051779B (en) | Numerical simulation method of three-dimensional direct-current resistivity method | |
| CN111638556A (en) | Magnetotelluric forward modeling method and device based on geospatial solution strategy and storage medium | |
| CN106019386A (en) | Method for measuring apparent resistivity through ratios in loop source | |
| CN116842813A (en) | Three-dimensional geoelectromagnetic forward modeling method | |
| CN113204905A (en) | Contact induced polarization finite element numerical simulation method | |
| Drahor et al. | 3D resistivity imaging from an archaeological site in south‐western Anatolia, Turkey: a case study | |
| CN114047554B (en) | Earth resistivity model modeling method, apparatus, computer device and storage medium | |
| De Giorgi et al. | Passive and active electric methods: new frontiers of application | |
| CN113238285B (en) | Resistivity calculation method, system and terminal for geophysical charging method exploration | |
| Ye et al. | 2.5 D induced polarization forward modeling using the adaptive finite-element method | |
| CN113221412B (en) | Method, system, terminal and medium for calculating charging potential data of lossy medium | |
| CN113221411B (en) | Charging potential numerical simulation method, system and terminal for lossy medium with any shape | |
| Jia et al. | Modeling of complex geological body and computation of geomagnetic anomaly | |
| Hernández Contreras et al. | A computational package to compute the electrical resistivity tomography response for regular bodies immersed in a homogeneous half-space | |
| Turarova et al. | Evaluation of the 3D Topographic Effect of Homogeneous and Inhomogeneous Media on the Results of 2D Inversion of Electrical Resistivity Tomography Data | |
| CN111597752A (en) | Cross-hole resistivity CT deep learning inversion method and system for balancing inter-hole sensitivity |
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 |