General elliptic partial differential equations all have appropriate energy extremum problems, that is, variational problems are equivalent to them.
Take the right-hand rectangular coordinate system, with the X and Y axes on the ground and the Z axis vertically downward, so that the X axis direction is consistent with the structural strike, as shown in Figure 9-5. The conductivity of the medium is only a function of y and z, and it is assumed that the field source is a plane electromagnetic wave propagating vertically from top to bottom, thus ensuring that the boundary condition is also two-dimensional, that is, independent of the X axis. The actual unit system is adopted, the electromagnetic wave is taken as the harmonic field EI Ω T, and the three-dimensional wave equation is used.
Figure 9-5 Right Hand Cartesian Coordinate System
Solid geophysics: seismology, geoelectricity and geothermal energy
become
Solid geophysics: seismology, geoelectricity and geothermal energy
Where k is the wave number, u is y, the electric field e or magnetic field h in the region s on the z plane, and the boundary conditions are
Solid geophysics: seismology, geoelectricity and geothermal energy
г is the boundary of the domain S, that is, the field value on the boundary is called u0. In other words, its definite solution problem is
Solid geophysics: seismology, geoelectricity and geothermal energy
Mathematically, it can be proved that the variational problem equivalent to the above definite solution problem is: for functional
Solid geophysics: seismology, geoelectricity and geothermal energy
Minimize, that is, δ j (u) = 0.
(2) Subdivision unit
Here, dividing the region into small units is more complicated than one-dimensional problem. Usually, the S region is divided into many triangular units by some dividing lines. Two adjacent triangles have a pair of common vertices and a common edge. When a curve boundary is encountered, it can be approximated by the straight edge of a triangle. If there are different media in the domain, you can choose the dividing line as the dividing line, which is called the natural boundary.
The requirements of subdivision are: (1) Each cell can only intersect with one vertex, as shown in Figure 9-6, and point A is not allowed, because point A in △bcf is not a vertex; ② The field values of the same vertex in different units should be equal; (3) The vertex of each unit is called a node, and the node numbering sequence of each unit should be consistent, generally counterclockwise, and each unit also has a number.
Figure 9-6 Rules for Dividing Units
The value of the function u(y, z) at the node is called the node parameter.
So far, we have discretized the whole field S, so that we can get equation (9-60).
Solid geophysics: seismology, geoelectricity and geothermal energy
here
Solid geophysics: seismology, geoelectricity and geothermal energy
Where m is the number of cells.
(3) unit analysis
1. area coordinates
In order to relate the field value of any point with the coordinates of that point on the plane by simple relations (such as linear relations), we introduce the concept of area coordinates.
Investigate how the position of any point P in a triangle is expressed by area.
Let the three vertices of a triangle on the YZ plane be
I(yi, zi), j(yi, zj) and k(yk, zk), as shown in figure 9-7.
Figure 9-7 Concept of Area Coordinates
Now define
Solid geophysics: seismology, geoelectricity and geothermal energy
Among them:
Solid geophysics: seismology, geoelectricity and geothermal energy
Is the area of a triangle; △i is the area of the small triangle pjk corresponding to point I; Δ j and Δ k are similar to Δ i, which are the areas of triangle ipk and triangle ijp corresponding to points J and K respectively.
Therefore, the position of the point P in the triangular element can be represented by φi, φj and φk, and the sizes of φi, φj and φk determine the position of the point in the element. Therefore, φi, φj and φk are called area coordinates, and they have the following relationship:
Solid geophysics: seismology, geoelectricity and geothermal energy
In addition, after the triangle area is expressed by determinant, the area coordinates can be expressed as
Solid geophysics: seismology, geoelectricity and geothermal energy
These include:
Solid geophysics: seismology, geoelectricity and geothermal energy
When point P falls on point I, that is, when the coordinate of point P is (yi, zi), there are
Solid geophysics: seismology, geoelectricity and geothermal energy
If point P falls on point J or point K, there are respectively
Solid geophysics: seismology, geoelectricity and geothermal energy
Their images are shown in Figure 9-8.
Equations (9-67), (9-68) and (9-69) can be summarized as follows.
Solid geophysics: seismology, geoelectricity and geothermal energy
Figure 9-8 Schematic diagram of unit area coordinates
2. Representation of field values in cells
Let the field value on each node be ui (I = 1, 2, ..., n), then in any cell E, the node number is I, J, K, and the node parameter is UI, uj, uk, so the field value of any point in cell E can be expressed by the node field value and regional coordinates:
Solid geophysics: seismology, geoelectricity and geothermal energy
here
Solid geophysics: seismology, geoelectricity and geothermal energy
The partial derivative of the field can be expressed as
Solid geophysics: seismology, geoelectricity and geothermal energy
3. Calculation of JE (u)
It can be obtained from formula (9-70) and formula (9-7 1).
Solid geophysics: seismology, geoelectricity and geothermal energy
These include:
Solid geophysics: seismology, geoelectricity and geothermal energy
therefore
Solid geophysics: seismology, geoelectricity and geothermal energy
These include:
Solid geophysics: seismology, geoelectricity and geothermal energy
K is the wave number.
At this point, the element stiffness matrix Je(u) can be obtained, that is
Solid geophysics: seismology, geoelectricity and geothermal energy
4. Comprehensive integration
In order to "integrate", it is necessary to expand the element stiffness matrix into a large matrix of N×N order (n is the number of nodes), and then put each element stiffness matrix in the corresponding position.
Solid geophysics: seismology, geoelectricity and geothermal energy
5. Solve the variational equation
The variational problem of δ j (u) = 0 can be approximated as
Solid geophysics: seismology, geoelectricity and geothermal energy
but
Solid geophysics: seismology, geoelectricity and geothermal energy
Considering that u(y, z) in the original problem should satisfy the boundary condition u(y, z) | г = u0, the matrix u can be written in block form:
Solid geophysics: seismology, geoelectricity and geothermal energy
Solid geophysics: seismology, geoelectricity and geothermal energy
Write q and t in the form of blocks accordingly, and you have it.
Solid geophysics: seismology, geoelectricity and geothermal energy
in consideration of
Solid geophysics: seismology, geoelectricity and geothermal energy
If the matrix p is symmetric, then
Solid geophysics: seismology, geoelectricity and geothermal energy
So,
Solid geophysics: seismology, geoelectricity and geothermal energy
Because q and t are symmetrical, there is
Solid geophysics: seismology, geoelectricity and geothermal energy
therefore
Solid geophysics: seismology, geoelectricity and geothermal energy
manufacture
Solid geophysics: seismology, geoelectricity and geothermal energy
Then the above formula can be written as
Solid geophysics: seismology, geoelectricity and geothermal energy
Represented by a unified system of linear equations
Solid geophysics: seismology, geoelectricity and geothermal energy
The coefficient matrix of the above linear equations is symmetric positive definite matrix and sparse matrix, and the specific solution can be solved by direct method, so I won't repeat it here.
This is the main step of the finite element method, and here is the simplest method. The unit division method can also adopt the rectangular quadrilateral division method; The simple relationship between field quantity and coordinates can be approximated by quadratic curve besides linear relationship.
6. Treatment of boundary
Because there are some conceptual differences between the boundary value problems proposed in the finite element method and the boundary conditions used in the electrical method, it is necessary to mention them first in order to avoid confusion in expression. Judging from the customary formulation of electrical prospecting, the so-called boundary generally refers to the interface with different electrical properties. As for the boundary of the studied domain, it is generally called limit condition, such as infinite condition, field source condition, surface condition and so on. However, in the finite element method, boundary problems are usually divided into two categories: one is the physical interface in the research domain, and all those belonging to this boundary are called free boundary conditions. Because the physical interface must be a boundary of the unit when dividing the unit, this processing method ensures the continuity of the field values on both sides of the electric boundary, that is to say, the division and interpolation have ensured the conditions of the field on the boundary, so it is not necessary to take special consideration in the specific operation. According to the finite element method, the physical requirements on the boundary will naturally be met. The other is the so-called mandatory boundary condition, which actually refers to the boundary value of the research domain. At this time, when using boundary conditions, it can be divided into three situations: ① Take the first kind of boundary value, that is, Dirichlet problem. The field value on the boundary is known in a specific solution; ② Take the second boundary value, that is, the Neumann problem, and know the normal derivative of the field value on the boundary; (3) The mixture of the above two is the so-called mixed boundary value problem.
Figure 9-9 Schematic Diagram of Boundary Selection Principles
When we consider the numerical calculation of magnetotelluric field, we mainly want to solve the problem of uneven level. At this time, we regard the horizontal layer as the background field, and the field caused by the inhomogeneity in the horizontal layer stratum is called the abnormal field, or the secondary field. In this way, we take the horizontally layered space field value that has been solved analytically as the boundary value for calculating the abnormal secondary field. As shown in Figure 9-9, the selection principle of boundaries at this time should make the boundaries such as AC, CD and DB far enough away from the abnormal body, and consider that the secondary field is zero on these boundaries. There is no doubt that this assumption is easy to do for AC, CD and DB boundaries. For AB surface (that is, the surface), this selection requirement is meaningless, because for the purpose of detection, we can only understand the underground anomalies by observing the surface anomalies, so it is unrealistic to assume that the surface anomalies are zero. At this time, the method we deal with is that the fixed solution domain of AB edge extends upward to a height of several kilometers from the surface. We think that the abnormal secondary field is zero above this height, which is equivalent to an extra layer of air in the studied profile.
When considering the boundary treatment of magnetotelluric field, it is assumed that the field source is a plane electromagnetic wave with vertical incidence, and the section studied is a two-dimensional horizontal inhomogeneous body, which is included in a certain layer of the horizontal layered section.