CN115238550B - Land electric field numerical simulation calculation method for landslide rainfall of self-adaptive unstructured grid - Google Patents

Land electric field numerical simulation calculation method for landslide rainfall of self-adaptive unstructured grid Download PDF

Info

Publication number
CN115238550B
CN115238550B CN202210880653.8A CN202210880653A CN115238550B CN 115238550 B CN115238550 B CN 115238550B CN 202210880653 A CN202210880653 A CN 202210880653A CN 115238550 B CN115238550 B CN 115238550B
Authority
CN
China
Prior art keywords
electric field
grid
calculating
landslide
saturation
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
Application number
CN202210880653.8A
Other languages
Chinese (zh)
Other versions
CN115238550A (en
Inventor
鲁光银
白冬鑫
朱自强
张升彪
蒋召杰
彭德明
刘智
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hunan Zhili Engineering Science And Technology Co ltd
Central South University
Original Assignee
Hunan Zhili Engineering Science And Technology Co ltd
Central South University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hunan Zhili Engineering Science And Technology Co ltd, Central South University filed Critical Hunan Zhili Engineering Science And Technology Co ltd
Priority to CN202210880653.8A priority Critical patent/CN115238550B/en
Publication of CN115238550A publication Critical patent/CN115238550A/en
Application granted granted Critical
Publication of CN115238550B publication Critical patent/CN115238550B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13Differential equations
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/10Geometric CAD
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20Finite element generation, e.g. wire-frame surface description, tesselation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling

Abstract

The invention discloses a land electric field numerical simulation calculation method of landslide rainfall of a self-adaptive unstructured grid, which comprises the following steps: an ERT monitoring device is deployed on the landslide, and a geometric model of the landslide is built; performing initial subdivision of triangular meshes of a ground electric field and a seepage field in a geometric model; establishing a landslide rainfall infiltration control equation and carrying out iterative solution to enable the triangular mesh of the initial seepage field to be self-adaptively and dynamically adjusted to meet the posterior error requirement, and calculating to obtain the saturation at all moments; calculating the saturation of each node of the ground electric field grid by using the seepage field grid nodes, and then calculating to obtain the conductivity of each unit of the ground electric field grid; establishing a control equation of a three-dimensional power supply electric field, solving in an iterative manner, traversing electrode arrangement to obtain apparent resistivity of all point power supplies at all moments, and enabling an initial current field triangular grid to be self-adaptively and dynamically adjusted to meet the posterior error requirement; and (5) performing visual display. Therefore, the space-time evolution process of the water content of the landslide in the rainfall process is carried out through numerical simulation.

Description

Land electric field numerical simulation calculation method for landslide rainfall of self-adaptive unstructured grid
Technical Field
The invention relates to the technical field of geophysical monitoring, in particular to a land electric field numerical simulation calculation method of landslide rainfall of a self-adaptive unstructured grid.
Background
Rainfall infiltration is an important inducing factor causing landslide damage, and effective mastering of the migration rule of groundwater in the rainfall process is important for understanding the inherent mechanism of landslide damage and performing early warning treatment. In the rainfall process, the hydrologic characteristics of the saturated region and the unsaturated region of the landslide soil body are continuously changed, and the process is saturated-unsaturated seepage. In the past, the space-time evolution monitoring of water in the saturated-unsaturated seepage process often adopts a mode of deploying a hydrologic monitoring sensor, and the point-based mode only can reflect hydrologic conditions of partial positions, is not necessarily representative, and often needs to drill holes, so that the cost is high, and the original structure and seepage path of soil are damaged.
Because of the correlation between soil resistivity and water content, geophysical techniques, represented by electric methods, are also widely used for detection of groundwater and hydrologic environments, and are mainly used in comparison with direct current methods, natural electric methods, transient electromagnetic methods, ground penetrating radars, nuclear magnetic resonance and the like. The direct current method is typically represented by a resistivity tomography (ERT) technology, the technology reflects the internal environment by acquiring a apparent resistivity image on the ground surface, and the earth electric field response acquired by the ERT technology can effectively represent the hydrologic process in the landslide soil body, so that the monitoring purpose is achieved. The core is the water-electricity relationship with the middle as a tie, which is important for understanding the whole process. The current research on hydropower relations is mostly based on experience, wherein an Archie formula is representative, but a definite constitutive model is not yet available, because the soil resistivity parameter is related to soil structure, water content, soil composition and the like, and the uncertainty exists. This uncertainty directly affects the link between the seepage field and the ground field, for which many studies have conducted many controlled indoor experiments. Numerical modeling is also an available research method in addition to experiments, but currently less research is being done in this regard.
There have been studies to simulate the water permeation process by decreasing the resistivity of a portion of the model. The mode can better reflect the earth electric field response of the seepage process, but simplifies the space-time evolution process of the moisture. The numerical simulation basically adopts fixed structured grids, which are not applicable to landslide research, firstly, because the landslide model is a very complex structure, and the structured grids do not have good adaptability to the complex structure. Secondly, the influence of the fixed grid on the solving precision is quite large, and if the initial given grid is unsuitable, the solving precision is difficult to improve.
Disclosure of Invention
First, the technical problem to be solved
Based on the problems, the invention provides a land electric field numerical simulation calculation method of self-adaptive unstructured grid landslide rainfall, which solves the problem of the space-time evolution process of the water content of the landslide in the rainfall process through numerical simulation, and in the simulation process, the self-adaptive unstructured grid is adopted, so that the space-time evolution process of the water content of the landslide in the rainfall process is more real and more suitable for reflecting.
(II) technical scheme
Based on the technical problems, the invention provides a land electric field numerical simulation calculation method of landslide rainfall of a self-adaptive non-structural grid, which comprises the following steps:
s1, deploying ERT monitoring equipment on a landslide, collecting monitoring data, collecting rock parameters and geological data on the landslide, and establishing a geometric model of the landslide through software; performing initial subdivision of the triangular mesh of the ground electric field in the geometric model through software;
s2, carrying out subdivision of an initial seepage field triangular mesh in the geometric model through software, establishing a landslide rainfall infiltration control equation and carrying out iterative solution, enabling the initial seepage field triangular mesh to be self-adaptively and dynamically adjusted to meet the posterior error requirement, and calculating to obtain the total water head, pore water pressure, saturation and volume water content at all moments;
s3, calculating the saturation of each node of the ground electric field grid by using the saturation data of the grid nodes of the seepage field grid as a reference through an inverse distance interpolation algorithm, and then sequentially calculating the resistivity of each node of the ground electric field grid and the conductivity of each unit of the ground electric field grid by using the saturation of each node of the ground electric field grid;
s4, establishing a control equation of a three-dimensional power supply electric field, carrying out iterative solution, calculating electric potential according to the conductivity of each unit of the ground electric field grid, traversing electrode arrangement to obtain apparent resistivity of all point power supplies at all moments, and enabling the initial current field triangular grid to be self-adaptively and dynamically adjusted to meet the requirement of posterior error;
and S5, visually displaying the saturation at all the moments and the apparent resistivity of the power supply at all the points at all the moments.
Further, the step S2 includes:
s21, performing initial triangle mesh subdivision of a seepage field in a geometric model through software;
s22, establishing a landslide rainfall infiltration control equation, converting into a synthesized matrix equation about the total water head according to a linear interpolation shape function and finite element analysis of triangle unit discrete subdivision, calculating to obtain a matrix coefficient K, M, F of the matrix equation about the total water head at the current moment by using an initial permeability coefficient, and calculating the total water head H at the next moment according to an implicit differential equation t+Δt
The landslide rainfall infiltration control equation:
boundary conditions:
wherein H is the total water head, ρ w G is gravity acceleration, k is liquid density x ,k y ,k z The permeability coefficients in the x, y and z directions,slope of characteristic curve of soil and water, theta is volume water content, mu a Is the air pressure of the gap, mu w Mu, pore water pressure aw For matrix suction, in saturation state +.>Equal to 0; Γ -shaped structure 1 Is the boundary section of the known head, H 1 Is a known head; Γ -shaped structure 2 For boundary segments of known outgoing or incoming flows, q is the boundary flow, +.>Normal vector of boundary;
the linear interpolation shape function:
wherein delta is the area of the triangle unit, a i =x j y k -x k y j ,b i =y j -y k ,c i =x k -x j ,a j =x k y i -x i y k ,b j =y k -y i ,c j =x i -x k ,a k =x i y j -x j y i ,b k =y i -y j ,c k =x j -x i
And the landslide rainfall infiltration control equation after the finite element analysis:
wherein Ω is a calculation region, Γ is a traffic boundary, 1 if a node is a boundary node, or 0 if not;
the matrix equation for total head:
the implicit differential equation:
s23, judging whether the posterior error meets the requirement, if so, entering a step S24, and if not, updating the permeability coefficient k x ,k y ,k z Returning to the previous step for recalculation;
s24, calculating the pore water pressure, the volume water content and the saturation at the current moment according to the total water head at the current moment, and continuously calculating the pore water pressure, the volume water content and the saturation at the next moment until all the time steps.
Further, the calculation formula of the air gap water pressure is as follows: mu (mu) w =ρ w g(H-z),
The calculation formula of the volume water content is as follows:
the saturation is equal to the volume moisture content divided by the known porosity,
in theta r For residual water content, θ s And (3) the saturated water content, and a, n and m are model fitting parameters.
Further, the step S3 includes the following steps:
s31, calculating the saturation of each node of the ground electric field grid by using the saturation data of the grid nodes of the seepage field grid as a reference through an inverse distance interpolation algorithm;
s32, calculating the resistivity of each node of the ground electric field grid according to a landslide water-electricity conversion formula by using the saturation of each node of the ground electric field grid; the landslide water-electricity conversion formula is as follows:
wherein ρ is 0 For initial resistivity S 0 The saturation is the initial state saturation, n is the saturation index, and parameters K and n are obtained through fitting measured data;
s33, calculating the resistivity and the conductivity of each unit of the ground electric field grid according to the resistivity of each node of the ground electric field grid, wherein the average value of the resistivity of three nodes of one unit is the resistivity of the unit, and the conductivity of the unit is the inverse of the resistivity of the unit.
Further, the step S4 includes the following steps:
s41, calculating electric field distribution of all point power supplies according to the conductivity of each unit of the ground electric field grid: establishing a control equation of a three-dimensional power supply electric field, converting the control equation into a synthesized matrix equation related to the auxiliary potential according to a linear interpolation shape function and finite element analysis of discrete subdivision of a triangle unit, solving the auxiliary potential of the optimal wave number by adopting a differential evolution algorithm, and then obtaining a real potential by Fourier inverse transformation;
the control equation of the three-dimensional power supply electric field is as follows:
wherein sigma is medium conductivity, namely conductivity of each unit of the ground electric field grid, mu is potential, I is current intensity, (x) 0 ,y 0 ,z 0 ) Delta is a dirac function related to a power supply point, r is a boundary-to-point power supply distance, c is a constant, n is an out-of-boundary normal vector omega is a research area, and Γ is s Γ is the surface boundary for contact with air Is an infinity boundary;
assuming that the medium conductivity does not change in the y-direction, it is transformed by fourier:
wherein K is wave number, U is a pair-type potential, K 1 Correcting Bessel function, K for the first order of the second class 0 Correcting a Bessel function for the second class zero order;
the boundary value problem of the formula is:
the method comprises the following steps of performing interpolation calculation and conversion to obtain:
the matrix equation for the paraelectric potential is: ku=p
Wherein K is an overall coefficient matrix, U is an unknown auxiliary potential vector to be calculated, and P is a column vector related to a point power supply;
the objective function is minimum:
wherein r is j (j=1, 2, …, M) is the electrode pitch, M is the number of electrode pitches, N is the number of wavenumbers, k i ≥0,g i ≥0,1≤i≤N;
Solving the equation to calculate the paid potential of a certain wave number, and then carrying out Fourier inverse transformation by adopting the following equation to obtain the real potential:
k in i G is the Fourier transform coefficient i Is an inverse Fourier transform coefficient, μ (x, y 0 Z) is a section y=y 0 Is a true potential of (2);
s42, traversing the electrode arrangement to obtain total potential distribution of the electrode arrangement, measuring the potential difference of the electrodes, and calculating apparent resistivity to obtain apparent resistivity distribution at the current t moment;
wherein ρ is s In order to achieve a visual resistivity of the glass,k is the device coefficient, Δμ is the potential difference;
s43, judging whether the posterior error meets the requirement, if so, entering a step S5, otherwise, updating the wave number k, and returning to the step S41 for recalculation.
Further, the determination of the posterior error in steps S23 and S43 adopts the following method:
in the above description, e is k As a posterior error of the kth cell,for a true potential gradient, +.>For the potential gradient after recovery, +.>N i For the form function, P is defined within a small block of the same conductivity around the current cell, V is the number of cells around the current cell, and the small block is defined as follows:
Λ=Ω k ∪∑Ω jk =σ j
Ω k omega is the current cell j Is omega k Surrounding cells of equal resistivity to them,
the least square method is adopted to solve:
further, the visual presentation includes drawing a contour map, a change curve, or a dynamic video.
The invention also discloses a land electric field numerical simulation calculation system of the landslide rainfall of the self-adaptive unstructured grid, which comprises:
at least one processor; and at least one memory communicatively coupled to the processor, wherein:
the memory stores program instructions executable by the processor, the processor invoking the program instructions to perform the method.
A non-transitory computer-readable storage medium stores computer instructions that cause the computer to perform the method.
(III) beneficial effects
The technical scheme of the invention has the following advantages:
(1) According to the invention, saturation of grid nodes of the seepage field is obtained through the self-adaptively adjusted non-structural grid, conductivity of each unit of the earth electric field is obtained according to the relation between the seepage field and the earth electric field, apparent resistivity of each point power supply in the earth electric field is obtained through the self-adaptively adjusted non-structural grid, a space-time evolution process of moisture in a landslide rainfall process is obtained through a forward thought, and the self-adaptively dynamically adjusted non-structural grid is used for making the device suitable for a landslide model and high in precision;
(2) According to the invention, the seepage field and the ground electric field are divided into initial triangular grids, and then iteration solution is carried out through control equation establishment, interpolation calculation and finite element analysis until the requirements of posterior error are met, and the initial triangular grids which are not fixedly structured are subjected to self-adaptive dynamic adjustment, so that the initial triangular grids are continuously encrypted by grids close to the real seepage field or the ground electric field, thereby better approaching the real landslide topography and structure, and improving the solution precision;
(3) When the invention calculates the earth electric field, firstly, the earth electric field generated by the electrode is calculated by taking all electrodes as point sources, then all electrode arrangements are traversed, the earth electric field distribution of all power sources in each arrangement is obtained and added to be the earth electric field distribution of the arrangement, and finally, the apparent resistivity is calculated, so that each electrode taken as a point source only carries out one finite element, and the calculation efficiency is greatly improved; the potential distribution generated by the power supply at each point is also mutually independent, and the calculation process at each moment can be independently operated when the electric field calculation is carried out, so that the distributed calculation can be carried out through a plurality of computers, and the purpose of accelerating the calculation is achieved.
Drawings
The features and advantages of the present invention will be more clearly understood by reference to the accompanying drawings, which are illustrative and should not be construed as limiting the invention in any way, in which:
FIG. 1 is a flow chart of a method for simulating and calculating the value of a ground electric field of landslide rainfall of a self-adaptive unstructured grid according to the embodiment of the invention;
FIG. 2 is a schematic diagram of a geometric model of a landslide constructed in accordance with an embodiment of the present invention;
FIG. 3 is a schematic diagram of an unstructured grid, i.e., a triangular grid, dissected in accordance with an embodiment of the present invention;
FIG. 4 is a schematic diagram showing saturation changes at different moments during a rainfall process according to an embodiment of the present invention;
FIG. 5 is a graph showing changes in apparent resistivity at different times during a rainfall process according to an embodiment of the present invention.
Detailed Description
The following describes in further detail the embodiments of the present invention with reference to the drawings and examples. The following examples are illustrative of the invention and are not intended to limit the scope of the invention.
The embodiment of the invention discloses a land electric field numerical simulation calculation method of landslide rainfall of a self-adaptive unstructured grid, which is shown in figure 1 and comprises the following steps:
s1, deploying ERT monitoring equipment on a landslide, collecting monitoring data, collecting rock parameters and geological data on the landslide at the same time, establishing a geometric model of the landslide through software, and giving material properties, wherein the parameters are shown in the following table; performing initial ground electric field triangulation mesh generation in a geometric model through software, performing mesh encryption on the ground surface, particularly near electrodes, determining an electrode arrangement mode as shown in fig. 3, and giving an initial total water head distribution, a current I, a time interval delta t and the like;
s2, carrying out subdivision of an initial seepage field triangular mesh in the geometric model through software, establishing a landslide rainfall infiltration control equation and carrying out iterative solution, enabling the initial seepage field triangular mesh to be self-adaptively and dynamically adjusted to meet the posterior error requirement, and calculating to obtain the total water head, pore water pressure, saturation and volume water content at all moments;
s21, performing initial triangle mesh subdivision of a seepage field in a geometric model through software;
s22, establishing a landslide rainfall infiltration control equation, converting into a synthesized matrix equation about the total water head according to a linear interpolation shape function and finite element analysis of triangle unit discrete subdivision, calculating to obtain a matrix coefficient K, M, F of the matrix equation about the total water head at the current moment by using an initial permeability coefficient, and calculating the total water head H at the next moment according to an implicit differential equation t+Δt
The essence of the landslide rainfall infiltration process is saturated-unsaturated seepage, and the control equation is that
Wherein H is the total water head, ρ w Is the density of liquid, k x ,k y ,k z Is the permeability coefficient in three directions.Is the slope of the characteristic curve of soil and water, theta is the volume water content,μ a is the air pressure of the gap, mu w Mu, pore water pressure aw For matrix suction, in saturation state +.>Equal to 0.
The initial conditions for this equation are:
H(x,y,z,t)=H 0 (x,y,z,t 0 ) (2)
h in 0 The water head distribution at the initial moment;
boundary conditions mainly include two classes, head boundary conditions and flow boundary conditions, head boundary conditions are also first class boundary conditions, meaning that the head distribution over the boundary is known, expressed as:
wherein F 1 Is the boundary section of the known head, H 1 Is a known head.
Traffic boundaries are a second type of boundary condition, i.e. the flux at the boundary is known, expressed as:
in the middle ofIs the normal vector of the boundary Γ 2 For a boundary segment of known outgoing or incoming traffic, q is the boundary traffic.
Discrete subdivision of the investigation region is performed using triangle units, three vertices i (x i ,y i ),j(x j ,y j ),k(x k ,y k ) In reverse order, linear interpolation is adopted in the units, and the shape function is
Wherein delta is the area of the triangle unit, a i =x j y k -x k y j ,b i =y j -y k ,c i =x k -x j ,a j =x k y i -x i y k ,b j =y k -y i ,c j =x i -x k ,a k =x i y j -x j y i ,b k =y i -y j ,c k =x j -x i
The Galerkin method is adopted to carry out finite element analysis, and the finite element format of the saturated-unsaturated seepage equation can be obtained as follows:
where Ω is the calculation region, Γ is the flow boundary,
in the formula, if the node is a boundary node, the node is 1, otherwise, the node is 0.
The triangle units are synthesized in a total way, and then
For the iterative solution of the time differential term in the above formula by adopting a differential mode, the invention adopts an implicit differential mode as follows
S23, judging whether the posterior error meets the requirement, if so, entering a step S24, and if not, updating the permeability coefficient k x ,k y ,k z Returning to the previous step for recalculation;
and (3) performing posterior error judgment once through each calculation of the total water heads in the steps S22 and S23, and performing next calculation of the total water head after the grids are encrypted, so that each node of the triangular grid of the initial seepage field is sequentially adjusted, and the grid encryption enables the triangular grid of the seepage field to be more in line with the actual seepage field.
S24, calculating the pore water pressure, the volume water content and the saturation at the current moment according to the total water head at the current moment, and continuously calculating the pore water pressure, the volume water content and the saturation at the next moment until all the time steps;
the total water head of each time step can be calculated by solving the equation, and then the pore water pressure is calculated by the following formula:
μ w =ρ w g(H-z) (12)
after the pore water pressure is obtained, the volume water content can be solved according to the soil water characteristic curve, then the saturation can be calculated according to the porosity, and the saturation is equal to the volume water content divided by the known porosity. The soil-water characteristic curve is fitted by adopting a VG model, and the calculation formula is as follows:
in the above formula, theta is the volume water content to be calculated, and theta r For residual water content, θ r Is saturated water content, a, n, m are model fitting parameters, ρ w Density of water, g is gravity acceleration, mu w Is pore water pressure.
S3, calculating the saturation of each node of the ground electric field grid by using the saturation data of the grid nodes of the seepage field grid as a reference through an inverse distance interpolation algorithm, and then sequentially calculating the resistivity of each node of the ground electric field grid and the conductivity of each unit of the ground electric field grid by using the saturation of each node of the ground electric field grid;
s31, calculating the saturation of each node of the ground electric field grid by using the saturation data of the grid nodes of the seepage field grid as a reference through the existing inverse distance interpolation algorithm;
s32, calculating the resistivity of each node of the ground electric field grid according to a landslide water-electricity conversion formula by using the saturation of each node of the ground electric field grid;
s33, calculating the resistivity and the conductivity of each unit of the ground electric field grid according to the resistivity of each node of the ground electric field grid, wherein the average value of the resistivity of three nodes of one unit is the resistivity of the unit, and the conductivity of the unit is the inverse of the resistivity of the unit;
as a bridge between the percolation field and the ground electric field, the relationship between saturation and resistivity can be expressed as:
ρ in the formula (13) is the resistivity of the rock-soil body;is porosity; m becomes the cementation index; ρ w For the pore water resistivity, S is the saturation, n is the saturation index, and a is the fitting parameter. For the same rock-soil body, the pore structure and the physical property of pore water are not changed in the seepage process, and the Archie formula can be simplified as follows under the condition of only considering the saturation change:
/>
ρ in formula (14) 0 For initial resistivity S 0 For the saturation of the initial state, n is a saturation index, and for the measured data, the relation between the saturation and the resistivity can be obtained by fitting calculation parameters K and n.
S4, establishing a control equation of a three-dimensional power supply electric field, carrying out iterative solution, calculating electric potential according to the conductivity of each unit of the ground electric field grid, traversing electrode arrangement to obtain apparent resistivity of all point power supplies at all moments, and enabling the initial current field triangular grid to be self-adaptively and dynamically adjusted to meet the requirement of posterior error;
s41, calculating electric field distribution of all point power supplies according to the conductivity of each unit of the ground electric field grid: establishing a control equation of a three-dimensional power supply electric field, converting the control equation into a synthesized matrix equation related to the auxiliary electric potential according to a linear interpolation shape function and finite element analysis of discrete subdivision of a triangle unit, directly calculating a group of optimal wave numbers by adopting a differential evolution algorithm, calculating auxiliary electric potentials of the optimal wave numbers, and then obtaining real electric potentials by Fourier inverse transformation;
the boundary problem of the three-dimensional point source electric field is as follows
Wherein sigma is medium conductivity, namely conductivity of each unit of the ground electric field grid, mu is potential, I is current intensity, (x) 0 ,y 0 ,z 0 ) Delta is a dirac function related to a power supply point, r is a boundary-to-point power supply distance, c is a constant, n is an out-of-boundary normal vector omega is a research area, and Γ is s Γ is the surface boundary for contact with air Is an infinity boundary.
The three-dimensional problem described above can be converted to a two-dimensional problem, commonly referred to as 2.5-dimensional, by fourier transformation, assuming that the medium conductivity does not change in the y-direction. The 2.5-dimensional control equation is as follows
Wherein K is wave number, U is a pair-type potential, K 1 Correcting Bessel function, K for the first order of the second class 0 The Bessel function is modified for the second zero order, and A is.
The boundary value problem of the formula is:
still adopt triangle unit to split area, the interpolation process is consistent with step S22, the first term integral in the formula (15) is calculated
In the formula (18)
Second term integral in calculation (17)
In (19)
Calculation of the third term integral of (17)
Calculation of the fourth integral of (17)
In the formula (21)
So equation (16) can be expressed as:
K 1e is related to the unit structure only, and is irrelevant to the wave number, boundary and point power supply position, so that K can be firstly calculated 1e Is then calculated for each permutation only K 2e And K is equal to 3e The K2 matrix and the K3 matrix are formed, and the purpose of accelerating calculation time is achieved by calculating the K1 matrix only once.
The solution to the variation can be found as a system of equations:
KU=P (24)
in the formula (23), K is an overall coefficient matrix, U is an unknown pair-type potential vector to be calculated, and P is a column vector related to a point power supply. Solving the equation to calculate the pair-type potential of a certain wave number, and then carrying out Fourier inverse transformation by adopting the following equation to obtain the real potential.
K in i G is the Fourier transform coefficient i Is an inverse Fourier transform coefficient, μ (x, y 0 Z) is a section y=y 0 Is a true potential of (c).
In the 2.5-dimensional problem, the choice of wavenumbers determines the accuracy of the calculation results, and generally the smaller the model, the larger the wavenumber. The calculation method of the invention is to select a group of proper electrode distances r j (j=1, 2, …, M) a set of optimal wavenumbers is directly calculated using a differential evolution algorithm to minimize the objective function in the equation.
In the formula (25), M is the electrode distance number, N is the wave number, and k i ≥0,g i ≥0,1≤i≤N,Is the objective function value.
S42, traversing the electrode arrangement to obtain total potential distribution of the electrode arrangement, measuring the potential difference of the electrodes, and calculating apparent resistivity to obtain apparent resistivity distribution at the current time t.
After the potential distribution is obtained, the apparent resistivity response can be calculated by
Rho in the above s For apparent resistivity, K is the device coefficient, Δμ is the potential difference, and I is the supply current, typically 1A.
Firstly, calculating the generated earth electric field of all electrodes serving as point sources, traversing all electrode arrangements, obtaining the earth electric field distribution of all power supplies in each arrangement, adding the earth electric field distribution to obtain the earth electric field distribution of the arrangement, and finally, calculating apparent resistivity, wherein each electrode serving as a point source only carries out one finite element, thereby greatly improving the calculation efficiency.
The calculation process at each moment can be independently operated when the electric field calculation is carried out, so that the distributed calculation can be carried out through a plurality of computers, and the purpose of accelerating the calculation is achieved. In addition, the potential distribution generated by the power supplies at all points is calculated independently, so that the program uses OpenMP to perform parallelization processing, and the calculation speed is greatly increased.
S43, judging whether the posterior error meets the requirement, if so, entering a step S5, if not, updating the wave number k, and returning to the step S41 for recalculation;
unlike step S2, after the apparent resistivity of all the point power supplies is calculated in step S41 and step S42, the posterior error is determined in step S43, because the calculation of the electric field can be performed on all the point power supplies in a parallel manner according to the superposition principle of the electric field, which is beneficial to accelerating the calculation speed.
The posterior error determination in steps S23 and S43 adopts the following method:
the core of the self-adaptive subdivision is to evaluate the posterior error of the unit, the SPR technology is adopted in the research, the posterior error of the unit is considered to be the L2 norm of the gradient and the recovery gradient, and the calculation formula of the posterior error of the unit is as follows:
in the above description, e is k As a posterior error of the kth cell,for a true potential gradient, +.>Is the potential gradient after recovery. Wherein the restored gradient calculation formula is
/>
N i For the form function, P is defined within a small block of the same conductivity around the current cell, V is the number of cells around the current cell, and the small block is defined as follows:
Λ=Ω k ∪∑Ω jk =σ j (30)
omega in the above k Omega is the current cell j Is omega k Surrounding cells of equal resistivity, P can be written as
The method can adopt a least square method to solve:
s5, carrying out visual display on the saturation at all the moments and the apparent resistivity of the power supply at all the points at all the moments;
drawing a specific key point groundwater hydrological parameter for the saturation contour map of each moment obtained in S2, and a ground electrical parameter change curve along with time for the apparent resistivity contour map of each moment obtained in S4, and making a saturation and apparent resistivity response dynamic video, wherein (a) in FIG. 4 is saturation of an initial state, (b) in FIG. 4 is saturation of a first day, (c) in FIG. 4 is saturation of a second day, (d) in FIG. 4 is saturation of a third day, and (e) in FIG. 4 is saturation of a fourth day, and (f) in FIG. 4 is a change curve along with time of saturation of four monitoring points ABCD in FIG. 2; as shown in fig. 5, the apparent resistivity changes at different times during the rainfall process are shown, wherein (a) in fig. 5 is the apparent resistivity of the initial state, (b) in fig. 5 is the apparent resistivity of the first day, (c) in fig. 5 is the apparent resistivity of the second day, (d) in fig. 5 is the apparent resistivity of the third day, and (e) in fig. 5 is the apparent resistivity of the fourth day, and (f) in fig. 5 is the variation curve of the apparent resistivity of the ABCD four monitoring points in fig. 2 with time.
Wherein, for the calculation of the ground electric field, if the generated power ground electric field of each arrangement is calculated every time the arrangement is calculated, most point power sources are repeatedly calculated many times, which causes great waste of calculation resources, and this process is time-consuming because the electric field calculation of each point power source is a finite element process. Therefore, firstly, the earth electric field generated by all electrodes is calculated as point sources, then all electrode arrangements are traversed, the earth electric field distribution of all power sources in each arrangement is obtained and added to be the earth electric field distribution of the arrangement, and finally, the apparent resistivity is calculated, thus each operation is carried outThe electrode of the point source is subjected to finite element only once, so that the calculation efficiency is greatly improved; second K 1e Is related to the unit structure only, and is irrelevant to the wave number, boundary and point power supply position, so that K can be firstly calculated 1e Is then calculated for each permutation only K 2e And K is equal to 3e The K2 matrix and the K3 matrix are formed, and the purpose of accelerating calculation time is achieved by calculating the K1 matrix only once. The calculation process at each moment can be independently operated when the electric field calculation is carried out, so that the distributed calculation can be carried out through a plurality of computers, and the purpose of accelerating the calculation is achieved. In addition, the potential distribution generated by the power supplies at all points is calculated independently, so that the OpenMP is used for parallelization, and the calculation speed is greatly increased.
Finally, it should be noted that the above-mentioned method for simulating and calculating the value of the earth electric field may be converted into software program instructions, which may be implemented by using a computing system including a processor and a memory, or may be implemented by computer instructions stored in a non-transitory computer readable storage medium. The integrated units implemented in the form of software functional units described above may be stored in a computer readable storage medium. The software functional unit is stored in a storage medium, and includes several instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) or a processor (processor) to perform part of the steps of the methods according to the embodiments of the present invention. And the aforementioned storage medium includes: a U-disk, a removable hard disk, a Read-Only Memory (ROM), a random access Memory (Random Access Memory, RAM), a magnetic disk, or an optical disk, or other various media capable of storing program codes.
In summary, the method for simulating and calculating the land electric field value of the landslide rainfall of the self-adaptive unstructured grid has the following beneficial effects:
(1) According to the invention, saturation of grid nodes of the seepage field is obtained through the self-adaptively adjusted non-structural grid, conductivity of each unit of the earth electric field is obtained according to the relation between the seepage field and the earth electric field, apparent resistivity of each point power supply in the earth electric field is obtained through the self-adaptively adjusted non-structural grid, a space-time evolution process of moisture in a landslide rainfall process is obtained through a forward thought, and the self-adaptively dynamically adjusted non-structural grid is used for making the device suitable for a landslide model and high in precision;
(2) According to the invention, the seepage field and the ground electric field are divided into initial triangular grids, and then iteration solution is carried out through control equation establishment, interpolation calculation and finite element analysis until the requirements of posterior error are met, and the initial triangular grids which are not fixedly structured are subjected to self-adaptive dynamic adjustment, so that the initial triangular grids are continuously encrypted by grids close to the real seepage field or the ground electric field, thereby better approaching the real landslide topography and structure, and improving the solution precision;
(3) When the invention calculates the earth electric field, firstly, the earth electric field generated by the electrode is calculated by taking all electrodes as point sources, then all electrode arrangements are traversed, the earth electric field distribution of all power sources in each arrangement is obtained and added to be the earth electric field distribution of the arrangement, and finally, the apparent resistivity is calculated, so that each electrode taken as a point source only carries out one finite element, and the calculation efficiency is greatly improved; the potential distribution generated by the power supply at each point is also mutually independent, and the calculation process at each moment can be independently operated when the electric field calculation is carried out, so that the distributed calculation can be carried out through a plurality of computers, and the purpose of accelerating the calculation is achieved.
Finally, it should be noted that: the above embodiments are only for illustrating the technical solution of the present invention, and are not limiting; although embodiments of the present invention have been described in connection with the accompanying drawings, various modifications and variations may be made by those skilled in the art without departing from the spirit and scope of the invention, and such modifications and variations fall within the scope of the invention as defined by the appended claims.

Claims (7)

1. A land electric field numerical simulation calculation method of self-adaptive unstructured grid landslide rainfall is characterized by comprising the following steps:
s1, deploying ERT monitoring equipment on a landslide, collecting monitoring data, collecting rock parameters and geological data on the landslide, and establishing a geometric model of the landslide through software; performing initial subdivision of the triangular mesh of the ground electric field in the geometric model through software;
s2, carrying out subdivision of an initial seepage field triangular mesh in the geometric model through software, establishing a landslide rainfall infiltration control equation and carrying out iterative solution, enabling the initial seepage field triangular mesh to be self-adaptively and dynamically adjusted to meet the posterior error requirement, and calculating to obtain the total water head, pore water pressure, saturation and volume water content at all moments;
s21, performing initial triangle mesh subdivision of a seepage field in a geometric model through software;
s22, establishing a landslide rainfall infiltration control equation, converting into a synthesized matrix equation about the total water head according to a linear interpolation shape function and finite element analysis of triangle unit discrete subdivision, calculating to obtain a matrix coefficient K, M, F of the matrix equation about the total water head at the current moment by using an initial permeability coefficient, and calculating the total water head H at the next moment according to an implicit differential equation t+Δt
The landslide rainfall infiltration control equation:
boundary conditions:(x,y,z)∈Γ 1
(x,y,z)∈Γ 2
wherein t represents the current time, t+Δt represents the next time, H is the total head, ρ w G is gravity acceleration, k is liquid density x ,k y ,k z The permeability coefficients in the x, y and z directions,slope of characteristic curve of soil and water, theta is volume water content, mu a Mu, pore air pressure w Mu, pore water pressure aw For matrix suction, in saturation state +.>Equal to 0; Γ -shaped structure 1 Is the boundary section of the known head, H 1 Is a known head; Γ -shaped structure 2 For boundary segments of known outgoing or incoming flows, q is the boundary flow, +.>Normal vector of boundary;
the linear interpolation shape function:r=i,j,k
wherein delta is the area of the triangle unit, a i =x j y k -x k y j ,b i =y j -y k ,c i =x k -x j ,a j =x k y i -x i y k ,b j =y k -y i ,c j =x i -x k ,a k =x i y j -x j y i ,b k =y i -y j ,c k =x j -x i
And the landslide rainfall infiltration control equation after the finite element analysis:
wherein Ω is a calculation region, Γ is a traffic boundary, 1 if a node is a boundary node, or 0 if not;
the matrix equation for total head:
the implicit differential equation:
s23, judging whether the posterior error meets the requirement, if so, entering a step S24, and if not, updating the permeability coefficient k x ,k y ,k z Returning to the previous step for recalculation;
s24, calculating the pore water pressure, the volume water content and the saturation at the current moment according to the total water head at the current moment, and continuously calculating the pore water pressure, the volume water content and the saturation at the next moment until all the time steps;
s3, calculating the saturation of each node of the ground electric field grid by using the saturation data of the grid nodes of the seepage field grid as a reference through an inverse distance interpolation algorithm, and then sequentially calculating the resistivity of each node of the ground electric field grid and the conductivity of each unit of the ground electric field grid by using the saturation of each node of the ground electric field grid;
s4, establishing a control equation of a three-dimensional power supply electric field, carrying out iterative solution, calculating electric potential according to the conductivity of each unit of the ground electric field grid, traversing electrode arrangement to obtain apparent resistivity of all point power supplies at all moments, and enabling an initial current field triangular grid to be self-adaptively and dynamically adjusted to meet the requirement of posterior error;
s41, calculating electric field distribution of all point power supplies according to the conductivity of each unit of the ground electric field grid: establishing a control equation of a three-dimensional power supply electric field, converting the control equation into a synthesized matrix equation related to the auxiliary potential according to a linear interpolation shape function and finite element analysis of discrete subdivision of a triangle unit, solving the auxiliary potential of the optimal wave number by adopting a differential evolution algorithm, and then obtaining a real potential by Fourier inverse transformation;
the control equation of the three-dimensional power supply electric field is as follows:
wherein sigma is medium conductivity, namely conductivity of each unit of the ground electric field grid, mu is potential, I is current intensity, (x) 0 ,y 0 ,z 0 ) Delta is a dirac function related to a power supply point, r is a boundary-to-point power supply distance, c is a constant, n is an out-of-boundary normal vector, Ω is a research area, and Γ is s Γ is the surface boundary for contact with air Is an infinity boundary;
assuming that the medium conductivity does not change in the y-direction, it is transformed by fourier:
wherein K is wave number, U is a pair-type potential, K 1 Correcting Bessel function, K for the first order of the second class 0 Correcting a Bessel function for the second class zero order;
the boundary value problem of the formula is:
the method comprises the following steps of performing interpolation calculation and conversion to obtain:
the matrix equation for the paraelectric potential is: ku=p
Wherein K is an overall coefficient matrix, U is an unknown auxiliary potential vector to be calculated, and P is a column vector related to a point power supply;
the objective function is minimum:
wherein r is j For the electrode spacing, j=1, 2, the contents of M, M is the number of electrode distances, N is the number of wave numbers, k i ≥0,g i ≥0,1≤i≤N;
Solving the equation to calculate the paid potential of a certain wave number, and then carrying out Fourier inverse transformation by adopting the following equation to obtain the real potential:
k in i G is the Fourier transform coefficient i Is an inverse Fourier transform coefficient, μ (x, y 0 Z) is a section y=y 0 Is a true potential of (2);
s42, traversing the electrode arrangement to obtain total potential distribution of the electrode arrangement, measuring the potential difference of the electrodes, and calculating apparent resistivity to obtain apparent resistivity distribution at the current t moment;
wherein ρ is s The apparent resistivity, K is the device coefficient, and Δμ is the potential difference;
s43, judging whether the posterior error meets the requirement, if so, entering a step S5, if not, updating the wave number k, and returning to the step S41 for recalculation;
and S5, visually displaying the saturation at all the moments and the apparent resistivity of the power supply at all the points at all the moments.
2. The method for simulating and calculating the ground electric field value of landslide rainfall of a self-adaptive unstructured grid according to claim 1, wherein the calculation formula of the pore water pressure is as follows: mu (mu) w =ρ w g(H-z),
The calculation formula of the volume water content is as follows:
the saturation is equal to the volume moisture content divided by the known porosity,
in theta r For residual water content, θ s And (3) the saturated water content, and a, n and m are model fitting parameters.
3. The method for simulating and calculating the ground electric field value of the landslide rainfall of the self-adaptive unstructured grid according to claim 2, wherein the step S3 comprises the following steps:
s31, calculating the saturation of each node of the ground electric field grid by using the saturation data of the grid nodes of the seepage field grid as a reference through an inverse distance interpolation algorithm;
s32, calculating the resistivity of each node of the ground electric field grid according to a landslide water-electricity conversion formula by using the saturation of each node of the ground electric field grid; the landslide water-electricity conversion formula is as follows:
wherein ρ is the resistivity of the rock-soil body, S is the saturation, ρ 0 For initial resistivity S 0 The saturation is the initial state saturation, n is the saturation index, and parameters K and n are obtained through fitting measured data;
s33, calculating the resistivity and the conductivity of each unit of the ground electric field grid according to the resistivity of each node of the ground electric field grid, wherein the average value of the resistivity of three nodes of one unit is the resistivity of the unit, and the conductivity of the unit is the inverse of the resistivity of the unit.
4. The method for simulating and calculating the ground electric field value of the landslide rainfall of the self-adaptive unstructured grid according to claim 3, wherein the posterior error is judged in the steps S23 and S43 by adopting the following method:
in the above description, e is k As a posterior error of the kth cell,for a true potential gradient, +.>For the potential gradient after recovery, +.>N i For the form function, P is defined within a small block of the same conductivity around the current cell, V is the number of cells around the current cell, and the small block is defined as follows:
Λ=Ω k ∪∑Ω jk =σ j
Ω k omega is the current cell j Is omega k Surrounding cells of equal resistivity to them,
the least square method is adopted to solve:
5. the method for simulating and calculating the ground electric field value of the landslide rainfall of the self-adaptive unstructured grid according to claim 4, wherein the visual display comprises drawing a contour map, a change curve or a dynamic video.
6. An electric field numerical simulation calculation system for landslide rainfall of a self-adaptive unstructured grid is characterized by comprising:
at least one processor; and at least one memory communicatively coupled to the processor, wherein:
the memory stores program instructions executable by the processor, the processor invoking the program instructions to perform the method of any of claims 1-5.
7. A non-transitory computer readable storage medium storing computer instructions that cause the computer to perform the method of any one of claims 1 to 5.
CN202210880653.8A 2022-07-25 2022-07-25 Land electric field numerical simulation calculation method for landslide rainfall of self-adaptive unstructured grid Active CN115238550B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210880653.8A CN115238550B (en) 2022-07-25 2022-07-25 Land electric field numerical simulation calculation method for landslide rainfall of self-adaptive unstructured grid

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210880653.8A CN115238550B (en) 2022-07-25 2022-07-25 Land electric field numerical simulation calculation method for landslide rainfall of self-adaptive unstructured grid

Publications (2)

Publication Number Publication Date
CN115238550A CN115238550A (en) 2022-10-25
CN115238550B true CN115238550B (en) 2024-02-13

Family

ID=83675777

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210880653.8A Active CN115238550B (en) 2022-07-25 2022-07-25 Land electric field numerical simulation calculation method for landslide rainfall of self-adaptive unstructured grid

Country Status (1)

Country Link
CN (1) CN115238550B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116227308B (en) * 2023-05-09 2023-07-18 广东石油化工学院 Numerical simulation method and system for shallow logging natural electric field
CN116502567B (en) * 2023-06-28 2023-09-12 中国空气动力研究与发展中心计算空气动力研究所 Interpolation solving method, device, equipment and medium of unstructured grid flow field
CN117113785B (en) * 2023-10-24 2024-01-09 长江大学武汉校区 Pore water pressure measuring method and device, electronic equipment and storage medium

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4792761A (en) * 1987-04-06 1988-12-20 King Thomas C Geophysical prospecting with collimated magnetotelluric fields
CN113406150A (en) * 2021-06-18 2021-09-17 湖南致力工程科技有限公司 Rainfall process active ground electric field time-space response simulation experiment device

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4792761A (en) * 1987-04-06 1988-12-20 King Thomas C Geophysical prospecting with collimated magnetotelluric fields
CN113406150A (en) * 2021-06-18 2021-09-17 湖南致力工程科技有限公司 Rainfall process active ground electric field time-space response simulation experiment device

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
LANCZOS迭代算法及其在三维地电场数值模拟计算中的应用研究;宛新林, 席道瑛, 高尔根;《物探化探计算技术》;第26卷(第1期);全文 *
基于对偶加权后验误差估计的2.5维直流电阻率自适应有限元正演;严波 等;《物探与化探》;第38卷(第1期);全文 *
稳定地电场三维有限差分正演模拟;邓正栋 等;《石油物探》;第40卷(第1期);全文 *

Also Published As

Publication number Publication date
CN115238550A (en) 2022-10-25

Similar Documents

Publication Publication Date Title
CN115238550B (en) Land electric field numerical simulation calculation method for landslide rainfall of self-adaptive unstructured grid
US10048226B2 (en) Imaging method and apparatus based on magnetic flux leakage testing
CN113221393B (en) Three-dimensional magnetotelluric anisotropy inversion method based on non-structural finite element method
CN112949134B (en) Earth-well transient electromagnetic inversion method based on non-structural finite element method
CN104835202A (en) Quick three-dimensional virtual scene constructing method
Nan et al. Groundwater parameter estimation using the ensemble Kalman filter with localization
Peters Modified conceptual model for compensated root water uptake–A simulation study
Pryet et al. 3D resistivity gridding of large AEM datasets: A step toward enhanced geological interpretation
Jeannot et al. A low-dimensional integrated subsurface hydrological model coupled with 2-D overland flow: Application to a restored fluvial hydrosystem (Upper Rhine River–France)
Yang et al. A novel algorithm with heuristic information for extracting drainage networks from raster DEMs
Xing et al. Bathymetry inversion using the modified gravity-geologic method: application of the rectangular prism model and Tikhonov regularization
CN114970289B (en) Three-dimensional magnetotelluric anisotropy forward modeling numerical simulation method, equipment and medium
Camporese et al. Comparison of data assimilation techniques for a coupled model of surface and subsurface flow
CN114155354B (en) Method and device for reconstructing capacitance tomography based on graph convolution network
Heidari et al. History matching of reservoir models by ensemble Kalman filtering: The state of the art and a sensitivity study
CN115508900A (en) Ground dragging type transient electromagnetic imaging method and system
CN115238565A (en) Resistivity model reconstruction network training method, electromagnetic inversion method and device
Erdal et al. The value of simplified models for spin up of complex models with an application to subsurface hydrology
Neumann et al. Comparing the" bathtub method" with Mike 21 HD flow model for modelling storm surge inundation
Young et al. Initialization and setup of the Coastal Model Test Bed: integrated bathymetry
Gáspár et al. Lagrangian modelling of the convective diffusion problem using unstructured grids and multigrid technique
CN111597752A (en) Cross-hole resistivity CT deep learning inversion method and system for balancing inter-hole sensitivity
Kumar et al. A GIS-based methodology for assigning a flux boundary to a Numerical Groundwater Flow Model and its effect on model calibration
CN111126439A (en) Method for generating coal mine underground two-dimensional mine earthquake imaging training data set
CN104764939B (en) The big plane iterative method of the upward depth conversion of ship underwater static electric field in deep-sea

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