Current location - Plastic Surgery and Aesthetics Network - Wedding planning company - Two-dimensional finite element method of magnetotelluric field
Two-dimensional finite element method of magnetotelluric field
(A) the establishment of the variational equation

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.