CN111210522B - Method for tracking streamline distribution in three-dimensional unstructured grid flow field by using FEM (finite element modeling) - Google Patents
Method for tracking streamline distribution in three-dimensional unstructured grid flow field by using FEM (finite element modeling) Download PDFInfo
- Publication number
- CN111210522B CN111210522B CN202010035803.6A CN202010035803A CN111210522B CN 111210522 B CN111210522 B CN 111210522B CN 202010035803 A CN202010035803 A CN 202010035803A CN 111210522 B CN111210522 B CN 111210522B
- Authority
- CN
- China
- Prior art keywords
- coordinate
- master element
- coordinates
- real
- space
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- 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
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T19/00—Manipulating 3D models or images for computer graphics
- G06T19/20—Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- Architecture (AREA)
- Computer Hardware Design (AREA)
- General Engineering & Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a method for utilizing FEM method for tracing streamline distribution in three-dimensional unstructured grid flow field, establishing compact sandstone reservoir numerical simulation model of three-dimensional unstructured grid, locating the position of the initial point of any streamline in tetrahedral grid system, and calculating the initial point Pi0Coordinates and velocities in Master Element space; obtaining the time and the exit coordinates of the particles leaving the tetrahedral mesh in the Master Element space; converting the outlet coordinates into real space coordinates, and connecting a starting point and an end point in the real space to obtain a streamline line segment in the grid; judgment of PeIf the flow reaches the sink in the flow field, the next flow line can be tracked if the flow reaches the sink, which indicates that the tracking of the flow line is finished, otherwise, P is usedeReplacement of Pi0Repeating the above steps to continue tracing the flow line until the flow line reaches the sink, and then tracing the next flow line. The method can accurately track the streamline distribution condition in the tetrahedral unstructured grid system.
Description
Technical Field
The invention relates to a method for tracking streamline distribution in a three-dimensional unstructured grid flow field by using FEM (finite element modeling), belonging to the technical field of unconventional oil and gas exploration and development.
Background
The numerical simulation based on the streamline is an engineering method for describing a dynamic physical process of fluid flow in a reservoir, and the method is widely applied to the aspects of petroleum engineering, underground water problems and the like. From the last 50 s, streamline simulations were introduced into the oil field, primarily to track fluid flow paths and visualize velocity distributions within the flow field. The accuracy of streamline tracing is important to solving the mass transfer problem along the streamline. In the past decades, streamline simulation methods have found wide application in numerical reservoir simulation, production optimization, history fitting, prediction of EOR production, etc. in the petroleum industry.
There are two general types of methods used to trace streamlines: Runge-Kutta streamline numerical tracking method and Pollock semi-analytic streamline tracking method. However, the Runge-Kutta streamline numerical tracking method needs to iteratively calculate the streamline after a given time step length is given, so that the Runge-Kutta streamline tracking efficiency is very low, a TOF grid is not easy to generate, and the running of subsequent streamline-based numerical simulation work is not facilitated. Therefore, Pollock is a streamline tracking method which is widely applied, wherein the Pollock tracking streamline is generally divided into two steps, and the pressure in each grid is solved through a numerical simulation method to obtain a pressure field and a speed field; then setting the starting point of the streamline, obtaining the speed of the starting point according to the speed field, calculating the time for the starting point to leave the grid, then tracking the next grid until the tracing of the streamline is finished, and connecting the starting point and the end point in all the grids, thus obtaining the streamline by tracing. Pollock proposes a piecewise analytical solution, derived by TOF, when assuming that the velocity of each direction in the velocity field varies linearly, and is independent of the velocities of other directions. Given the start point coordinates and velocity for any hypothetical fluid particle, the algorithm can determine the exit coordinates from the grid and the time it takes to pass through the grid (this is TOF, is the residence time of the particle within a grid). This class of methods relies on the structure, computational accuracy and efficiency of the grid. However, Pollock is suitable for structured grids, and cannot track streamlines when there are sources or sinks in the grid. Although streamline simulation in unstructured grids has been studied by many scholars, these algorithms suffer from at least one of the following disadvantages:
(1) these formulas rely on the concept of a control volume, requiring flow to be applied to each face of the mesh, and a certain flow in each mesh;
(2) in the process of tracing the streamline, due to the numerical solution of a differential equation, TOF grids cannot be effectively calculated and obtained;
(3) most of the unstructured grids are two-dimensional, and streamline tracing algorithms in the three-dimensional unstructured grids are few.
Therefore, the streamline obtained by the current method can only be used for visualizing the velocity field, and the one-dimensional mass transfer equation cannot be mapped on the streamline, so that the subsequent expansion of numerical simulation work based on the streamline is not facilitated. In order to carry out numerical simulation work based on flow lines, a new method for tracking the flow lines in a three-dimensional unstructured grid (tetrahedral) flow field needs to be developed.
Disclosure of Invention
Aiming at the problems, the invention mainly overcomes the defects in the prior art and provides a method for tracking streamline distribution in a three-dimensional unstructured grid flow field by using FEM.
The technical scheme provided by the invention for solving the technical problems is as follows: the method for tracking streamline distribution in the three-dimensional unstructured grid flow field by using FEM comprises the following steps:
s1, establishing a numerical simulation model of the compact sandstone reservoir of the three-dimensional unstructured grid by using a finite element method, and acquiring a pressure field and a velocity field of the compact sandstone reservoir according to parameters of the compact sandstone reservoir;
step S2, setting the number of flow lines, and determining the distribution of initial points according to the positions of the sources in the flow field;
step S3, for the starting point P of any streamline ii0(x0,y0,z0) Locating its position within the tetrahedral mesh system to obtain Pi0(x0,y0,z0) The Cell ID of the tetrahedral mesh is obtained, and the coordinates (x) of the tetrahedral mesh vertex of the mesh are obtained1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Velocity [ v ] corresponding to tetrahedral mesh vertexx1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]；
Step S4, converting the starting point P by the coordinate conversion formulai0(x0,y0,z0) And tetrahedral mesh vertex coordinates (x)1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Conversion to coordinates in Master Element ([ xi ])1,η1,ζ1)，(ξ2,η2,ζ2)，(ξ3,η3,ζ3)，(ξ4,η4,ζ4) Velocity [ v ] corresponding to tetrahedral mesh vertex through Jacobian matrixx1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]Conversion to speed in Master Element ([ v ]ξ1,vη1,vζ1],[vξ2,vη2,vζ2],[vξ3,vη3,vζ3],[vξ4,vη4,vζ4])；
Step S5, calculating the starting point P by using Shape Function of tetrahedral mesh and velocity of tetrahedral mesh vertex in Master Elementi0(x0,y0,z0) The speed of (d);
step S6, in the Master Element, solving an equation by using Shape Function and establishing an arbitrary particle velocity and position relation equation to obtain a track equation of the particles in the Master Element;
step S7, solving a track equation in the Master Element through a Newton-Raphson iterative algorithm in the Master Element to obtain four tetrahedral unit arrival particlesThe time of the surface is taken as the minimum positive value as the flight time delta tau in the unit, and the outlet coordinate P is calculatede *(ξe,ηe,ζe)；
Step S8, using Shape Function to coordinate the exit Pe *(ξe,ηe,ζe) Conversion to real space coordinates Pe(xe,ye,ze) Accumulating delta tau to establish TOF coordinates, and connecting a starting point and an end point in a real space to obtain a streamline line segment in the grid;
step S9, judgment Pe(xe,ye,ze) Whether to a sink within the flowfield;
step S10, if it is reached, it indicates that the tracing of the flow line is finished, then the next flow line can be traced, otherwise, P is usede(xe,ye,ze) Replacement of Pi0(x0,y0,z0) Repeating the steps S3-S9 to continue tracing the flow line until the flow line reaches the sink, and then tracing the next flow line.
The further technical proposal is that the vertex coordinates (x) of the tetrahedral mesh in the step S41,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Converting between a real space and a MasterElement space;
the coordinate conversion formula from the MasterElement space to the real space is as follows:
x＝x1(1-ξ-η-ζ)+x2ξ+x3η+x4ζ
y＝y1(1-ξ-η-ζ)+y2ξ+y3η+y4ζ
z＝z1(1-ξ-η-ζ)+z2ξ+z3η+z4ζ
the coordinate conversion formula from the real space to the MasterElement space is as follows:
wherein
A1＝(x1-x2),A2＝(x1-x3),A3＝(x1-x4)
B1＝(y1-y2),B2＝(y1-y2),B3＝(y1-y4)
C1＝(z1-z2),C2＝(z1-z3),C3＝(z1-z4)
In the formula: x is an x coordinate of a real space, y is a y coordinate of the real space, z is a z coordinate of the real space, and xi is a coordinate corresponding to the x coordinate in the MasterElement; eta is a coordinate corresponding to the y coordinate in the MasterElement; zeta is the coordinate corresponding to the z coordinate in MasterElement; (x)1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) The coordinates of four vertices in the tetrahedron in the real space are respectively converted into (0,0,0), (1,0,0), (0,1,0), (0,0,1) in the MasterElement.
The further technical scheme is that the velocity [ v ] corresponding to the vertex of the tetrahedral mesh in the step S4x1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]And (3) converting the real space into a MasterElement space, wherein the conversion formula is as follows:
wherein the Jacobian matrix is:
in the formula: upsilon isx、υyAnd upsilonyVelocity in the x, y and z directions, respectively, in real space, uξ、υηAnd upsilonζThe velocities in the ξ, η, and ζ directions in the MasterElement.
The further technical scheme is that the specific process of the step S7 is as follows:
step S71, starting from the starting point, the mass point is set to be able to leave the tetrahedron from any surface of the tetrahedron, and P is calculated separatelyi0(x0,y0,z0) From side 1 (P)2-P1-P4) Noodle 2 (P)2-P1-P3) Side 3 (P)3-P1-P4) Side 4 (P)2-P3-P4) Is a departure time Δ τ1、Δτ2、Δτ3、Δτ4；
Step S72, taking the leaving time Delta tau1、Δτ2、Δτ3、Δτ4Is the time of flight Δ τ within the cell;
step S73, calculating the exit coordinates in the MasterElement according to the flight time delta tau in the unit, and substituting the flight time delta tau in the unit into a case1 formula to obtain P through calculatione *(ξe,ηe,ζe)；
Where case1 is formulated as follows:
the invention has the following beneficial effects:
(1) in the streamline tracing process, only the attributes (coordinates and flow) of the vertices of the tetrahedral mesh need to be obtained, and the flow on each surface of the tetrahedron is not needed;
(2) the coordinates of the outlet of the flow line in the tetrahedron can be determined according to the coordinates of the inlet of the flow line, and the characteristic is similar to the Pollock algorithm;
(3) the invention provides an analytic solution of an equation describing a streamline track in a tetrahedral mesh, can accurately and efficiently determine the outlet coordinates of a streamline in a surface body mesh and the time of leaving the mesh according to the inlet coordinates, and accumulate delta tau so as to establish a TOF coordinate system;
(4) streamline tracing within a three-dimensional unstructured mesh (tetrahedral mesh) flow field can be achieved.
Drawings
FIG. 1 is a conversion of real space coordinates to MasterElement space coordinate codes in a FENICS system;
FIG. 2 is a block diagram of the conversion of Master Element space coordinates to real space coordinate codes in the FENICS system;
FIG. 3 is a conversion of real space velocity to MasterElement space velocity code in the FENICS system;
FIG. 4 is a tetrahedral entity map and coordinate distribution in real space and Master Element space;
FIG. 5 is a technical route diagram of a tight sandstone reservoir utilizing a finite element tracing streamline method in the embodiment of the invention;
FIG. 6 is a three-dimensional flow field discrete grid diagram in an embodiment of the present invention;
FIG. 7 is a pressure field diagram of a three-dimensional flow field in an embodiment of the present invention;
FIG. 8 is a diagram of a distribution of flow lines within a three-dimensional unstructured mesh (tetrahedral) flow field (three-dimensional view) in an embodiment of the present invention.
Detailed Description
The present invention will be further described with reference to the following examples and the accompanying drawings.
Example 1
The invention discloses a method for tracking streamline distribution in a three-dimensional unstructured grid flow field by using FEM, which comprises the following steps:
s1, establishing a Darcy flow numerical simulation model of the three-dimensional unstructured grid by using a finite element method, and acquiring a pressure field and a speed field of a three-dimensional flow field;
step S2, setting the number of the streamlines as 100, and uniformly distributing the starting points of the streamlines around the source (the radial quantity is 10, and the normal quantity is 10);
step S3, for the starting point P of any streamline ii0(x0,y0,z0) Locating its position within the tetrahedral mesh system to obtain Pi0(x0,y0,z0) The mesh number CellID of the tetrahedral mesh is obtained, and the tetrahedral mesh vertex coordinates (x) of the mesh are obtained1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) And corresponding velocity [ v ]x1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]；
Step S4, converting the starting point P by the coordinate conversion formulai0(x0,y0,z0) And tetrahedral mesh vertex coordinates (x)1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Conversion to coordinates in MasterElement (ξ)1,η1,ζ1)，(ξ2,η2,ζ2)，(ξ3,η3,ζ3)，(ξ4,η4,ζ4) And converting the corresponding velocity of the tetrahedral mesh vertex into the velocity ([ v ] in the MasterElement) through a Jacobian matrixξ1,vη1,vζ1],[vξ2,vη2,vζ2],[vξ3,vη3,vζ3],[vξ4,vη4,vζ4]) The translation code is as shown in FIG. 3;
the tetrahedral mesh vertex coordinates need to be transformed into each other in real space and MasterElement space,
the formula of coordinate conversion from MasterElement space to real space (the conversion code is shown in FIG. 2):
x＝x1(1-ξ-η-ζ)+x2ξ+x3η+x4ζ
y＝y1(1-ξ-η-ζ)+y2ξ+y3η+y4ζ
z＝z1(1-ξ-η-ζ)+z2ξ+z3η+z4ζ (1)
and the coordinate conversion formula from the real space to the MasterElement space is (the conversion code is shown in fig. 1):
wherein
A1＝(x1-x2),A2＝(x1-x3),A3＝(x1-x4)
B1＝(y1-y2),B2＝(y1-y2),B3＝(y1-y4)
C1＝(z1-z2),C2＝(z1-z3),C3＝(z1-z4)
In the formula: x is an x coordinate of a real space, y is a y coordinate of the real space, z is a z coordinate of the real space, and xi is a coordinate corresponding to the x coordinate in the MasterElement; eta is a coordinate corresponding to the y coordinate in the MasterElement; zeta is the coordinate corresponding to the z coordinate in MasterElement; (x)1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Coordinates of four vertices in a tetrahedron in real space are respectively converted into MasterElement coordinates of (0,0,0), (1,0,0), (0,1,0), (0,0,1), and (0,0,1), as shown in fig. 4;
in the streamline tracing process, the velocity of the tetrahedral mesh vertex only needs to be converted into Master Element space from real space, and the conversion formula is
Wherein the Jacobian matrix is:
in the formula: (x)1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Coordinates of four vertexes of a tetrahedron in a real space are respectively; upsilon isx、υyAnd upsilonyVelocity in the x, y and z directions, respectively, in real space, uξ、υηAnd upsilonζThe velocities in the xi, eta and zeta directions in the MasterElement;
step S5, calculating the starting point P in the MasterElement by using Shape Function of the tetrahedral mesh and the velocity of the vertex of the tetrahedral meshi0(x0,y0,z0) The speed of (d);
step S6, in the Master Element, establishing an arbitrary particle velocity and position relation equation by using Shape Function, and solving the equation to obtain a track equation of particles in the Master Element;
the trajectory equation solving process is as follows:
the Master element space is a regular tetrahedron, and P is given after any starting pointi0(x0,y0,z0) Then, Pi0(x0,y0,z0) Will leave from one face of a regular tetrahedron, the time taken is the time of flight TOF leaving the mesh, denoted Δ τ;
so the first-order lagrange shape function:
S1(ξ,η,ζ)＝1-ξ-η-ζ
S2(ξ,η,ζ)＝ξ
S3(ξ,η,ζ)＝η
S4(ξ,η,ζ)＝ζ (4)
velocity upsilon of arbitrary point P (xi, eta)ξIs composed of
υξ＝S1(ξ,η,ζ)υξ1+S2(ξ,η,ζ)υξ2+S3(ξ,η,ζ)υξ3+S4(ξ,η,ζ)υξ4
υη＝S1(ξ,η,ζ)υη1+S2(ξ,η,ζ)υη2+S3(ξ,η,ζ)υη3+S4(ξ,η,ζ)υη4
υζ＝S1(ξ,η,ζ)υζ1+S2(ξ,η,ζ)υζ2+S3(ξ,η,ζ)υζ3+S4(ξ,η,ζ)υζ4 (5)
Order to
So the equation can be written as
Conversion to matrix form
The general solution form of equation 8 is
X(t)＝E1Φ1(t)+E2Φ2(t)+E3Φ3(t)+C (9)
Wherein the characteristic function phi1(t),Φ2(t),Φ3(t) and coefficient E1,E2,E3and C is determined according to the eigenvalue of the matrix A, and X (t) has 9 different conditions according to the eigenvalue distribution condition of A;
case1 three unequal non-zero real number feature roots: (0 ≠ r)1≠r2≠r3≠0)
Case2 three non-zero real feature roots, where two roots are equal: (0 ≠ r)1≠r2＝r3≠0)
APc+B＝0 (11)
Case3 three equal non-zero real feature roots: (0 ≠ r)1≠r2＝r3≠0)
E1＝P0-Pc
E2＝(A-r1I)(P0-Pc)
APc+B＝0 (12)
Case4 one non-zero real feature root and two complex roots: (r)1≠0,μ+λi,λ≠0)
APc+B＝0 (13)
Case5 has a zero real feature root and two complex roots: (r)1＝0,μ+λi,λ≠0)
G2＝(μA-(μ2-λ2)I)
Case6 has one zero real characteristic root and two non-zero real roots: (r)1＝0,0≠r2≠r3≠0)
Case7 has two zero real characteristic roots and one non-zero real root: (r)1＝0,r2＝0,r3≠0)
E1＝AP0+B
Case8 has only a zero solution: (r)1＝0,r2＝0,r3＝0)
X(t)＝E1t+E2t2+E3t3+P0
E1＝AP0+B
Case9 has one real feature root of zero and two equal and non-zero real feature roots: (r)1＝0,r2＝0,r3＝0)
Step S7, solving a trajectory equation in the Master Element through a Newton-Raphson iterative algorithm in the Master Element to obtain the time of the particles to reach four surfaces of the tetrahedral unit, taking the minimum positive value as the flight time delta tau in the unit, and calculating an exit coordinate Pe*(ξe,ηe,ζe)；
The Master element is a regular tetrahedron (as shown in FIG. 4), and four vertices are P1(0,0,0),P2(1,0,0),P3(0,1,0),P4(0,0,1), four faces are: noodle 1 (P)2-P1-P4) Flour 2 (P)2-P1-P3) Flour 3 (P)3-P1-P4) Flour 4 (P)2-P3-P4) Particle from origin Pi0(x0,y0,z0) After a certain time has elapsed since the start of the operation, a point P on one face of the tetrahedron is reachedeThe time taken for the tetrahedron to leave is TOF, denoted by Δ τ;
according to the locus equation of the particle at Master Element, the outlet coordinate Pe*(ξe,ηe,ζe) And Δ τ can be calculated by the following four steps:
the first step is as follows: starting from the starting point, a particle is calculated P separately, assuming that the particle can depart from any one of the faces of a tetrahedroni0(x0,y0,z0) From side 1 (P)2-P1-P4) Flour 2 (P)2-P1-P3) Flour 3 (P)3-P1-P4) Flour 4 (P)2-P3-P4) Is a departure time Δ τ1、Δτ2、Δτ3、Δτ4；
Taking Case1 as an example, given any position, the starting point P can be obtained by a Newton-Raphson algorithmi0(x0,y0,z0) When this position is reached, given ξ ═ 0, Δ τ can be determined1(ii) a η ═ 0, Δ τ can be determined2(ii) a ζ is 0, Δ τ can be determined3(ii) a ξ + η + ζ ═ 1, Δ τ can be determined4；
The second step is that: based on the calculation result (departure time Δ τ)1、Δτ2、Δτ3、Δτ4) Finding out real delta tau; under the condition of giving the vertex coordinates and the speed of the tetrahedral mesh, TOF calculated according to a trajectory equation can only be a minimum positive value, and the rest delta tau is invalid;
the time of flight Δ τ in the cell is therefore min (Δ τ)1，Δτ2，Δτ3，Δτ4)；
The third step: calculating the exit coordinates in Master Element according to the flight time delta tau in the unit, and substituting the flight time delta tau in the unit into case1 formula to obtain Pe*(ξe,ηe,ζe)；
S8, mixing Pe*(ξe,ηe,ζe) Converting the actual space into the Shape Function to obtain an exit coordinate P in the actual spacee(xe,ye,ze) And accumulating the time of flight in the unit to obtain TOF (time of flight)；
The exit coordinate P in the Master Element spacee*(ξe,ηe,ζe) Conversion to real space exit coordinates Pe(xe,ye,ze) The formula of (a):
xe＝x1(1-ξe-ηe-ζe)+x2ξe+x3ηe+x3ζe
ye＝y1(1-ξe-ηe-ζe)+y2ξe+y3ηe+y3ζe
ze＝z1(1-ξe-ηe-ζe)+z2ξe+z3ηe+z3ζe；
s9, judgment Pe(xe,ye,ze) Whether to a sink within the flowfield;
s10, if the flow line is traced, the next flow line can be traced, otherwise, P is usede(xe,ye,ze) Replacement of Pi0(x0,y0,z0) Repeating steps S3 through S9 continues to trace the flow line until the flow line reaches a sink, and then subsequently tracing the next flow line.
The above steps can be summarized as a technical route map as shown in fig. 5, after a three-dimensional flow field is obtained (fig. 6 is a discretized grid map of the three-dimensional flow field, and fig. 7 is a pressure field map of the three-dimensional flow field), 100 streamlines are traced by taking an injection well as a starting point, and the streamline result is shown in fig. 8.
Although the present invention has been described with reference to the above embodiments, it should be understood that the present invention is not limited to the above embodiments, and those skilled in the art can make various changes and modifications without departing from the scope of the present invention.
Claims (3)
1. The method for tracking streamline distribution in the three-dimensional unstructured grid flow field by using FEM is characterized by comprising the following steps of:
s1, establishing a numerical simulation model of the compact sandstone reservoir of the three-dimensional unstructured grid by using a finite element method, and acquiring a pressure field and a velocity field of the compact sandstone reservoir according to parameters of the compact sandstone reservoir;
step S2, setting the number of flow lines, and determining the distribution of initial points according to the positions of the sources in the flow field;
step S3, for the starting point P of any streamline ii0(x0,y0,z0) Locating its position within the tetrahedral mesh system to obtain Pi0(x0,y0,z0) The Cell ID of the tetrahedral mesh is obtained, and the coordinates (x) of the tetrahedral mesh vertex of the mesh are obtained1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Velocity [ v ] corresponding to tetrahedral mesh vertexx1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]；
Step S4, converting the starting point P by the coordinate conversion formulai0(x0,y0,z0) And tetrahedral mesh vertex coordinates (x)1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Conversion to coordinates in Master Element ([ xi ])1,η1,ζ1)，(ξ2,η2,ζ2)，(ξ3,η3,ζ3)，(ξ4,η4,ζ4) The corresponding speed of the tetrahedral mesh vertex is determined by the Jacobian matrixDegree [ v ]x1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]Conversion to speed in Master Element ([ v ]ξ1,vη1,vζ1],[vξ2,vη2,vζ2],[vξ3,vη3,vζ3],[vξ4,vη4,vζ4])；
Step S5, calculating the starting point P by using Shape Function of tetrahedral mesh and velocity of tetrahedral mesh vertex in Master Elementi0(x0,y0,z0) The speed of (d);
step S6, in the Master Element, solving an equation by using Shape Function and establishing an arbitrary particle velocity and position relation equation to obtain a track equation of the particles in the Master Element;
the trajectory equation solving process is as follows:
the Master element space is a regular tetrahedron, and any starting point P is giveni0(x0,y0,z0) Then, Pi0(x0,y0,z0) Will leave from one face of a regular tetrahedron, the time taken is the time of flight TOF leaving the mesh, denoted Δ τ;
so the first-order lagrange shape function:
S1(ξ,η,ζ)＝1-ξ-η-ζ
S2(ξ,η,ζ)＝ξ
S3(ξ,η,ζ)＝η
S4(ξ,η,ζ)＝ζ (4)
in the formula: xi is a coordinate corresponding to the x coordinate in the Master Element; eta is a coordinate corresponding to the y coordinate in the Master Element; zeta is the coordinate corresponding to the z coordinate in the Master Element;
speed of arbitrary point P (xi, eta, zeta) is
υξ＝S1(ξ,η,ζ)υξ1+S2(ξ,η,ζ)υξ2+S3(ξ,η,ζ)υξ3+S4(ξ,η,ζ)υξ4
υη＝S1(ξ,η,ζ)υη1+S2(ξ,η,ζ)υη2+S3(ξ,η,ζ)υη3+S4(ξ,η,ζ)υη4
υζ＝S1(ξ,η,ζ)υζ1+S2(ξ,η,ζ)υζ2+S3(ξ,η,ζ)υζ3+S4(ξ,η,ζ)υζ4 (5)
Order to
So the equation can be written as
Conversion to matrix form
The general solution form of equation (8) is
X(t)＝E1Φ1(t)+E2Φ2(t)+E3Φ3(t)+C (9)
Wherein the characteristic function phi1(t)，Φ2(t)，Φ3(t) and coefficient E1，E2，E3C is determined according to the eigenvalue of the matrix A, and X (t) has 9 different conditions according to the eigenvalue distribution condition of the eigenvalue of A;
case1 three unequal non-zero real number feature roots:
case2 three non-zero real feature roots, where two roots are equal:
APc+B＝0 (11)
case3 three equal non-zero real feature roots:
E1＝P0-Pc
E2＝(A-r1I)(P0-Pc)
APc+B＝0 (12)
case4 one non-zero real feature root and two complex roots:
APc+B＝0 (13)
case5 has a zero real feature root and two complex roots:
G2＝(μA-(μ2-λ2)I)
case6 has one zero real characteristic root and two non-zero real roots:
case7 has two zero real characteristic roots and one non-zero real root:
E1＝AP0+B
case8 has only a zero solution:
X(t)＝E1t+E2t2+E3t3+P0
E1＝AP0+B
case9 has one real feature root of zero and two equal and non-zero real feature roots:
step S7, solving a trajectory equation in the Master Element through a Newton-Raphson iterative algorithm in the Master Element to obtain the time of the particles to reach four surfaces of the tetrahedral unit, taking the minimum positive value as the flight time delta tau in the unit, and calculating an exit coordinate Pe *(ξe,ηe,ζe)；
In step S71, the mass point is set to be able to move from any surface of the tetrahedron to four surfaces starting from the originBody, and separately calculate Pi0(x0,y0,z0) Time Δ τ of departure from surface 1, surface 2, surface 3, and surface 41、Δτ2、Δτ3、Δτ4；
Step S72, taking the leaving time Delta tau1、Δτ2、Δτ3、Δτ4Is the time of flight Δ τ within the cell;
step S73, calculating an exit coordinate in Master Element according to the flight time delta tau in the unit, and substituting the flight time delta tau in the unit into a case1 formula to obtain P through calculatione *(ξe,ηe,ζe)；
Step S8, using Shape Function to coordinate the exit Pe *(ξe,ηe,ζe) Conversion to real space coordinates Pe(xe,ye,ze) Accumulating delta tau to establish TOF coordinates, and connecting a starting point and an end point in a real space to obtain a streamline line segment in the grid;
step S9, judgment Pe(xe,ye,ze) Whether to a sink within the flowfield;
step S10, if it is reached, it indicates that the tracing of the flow line is finished, then the next flow line can be traced, otherwise, P is usede(xe,ye,ze) Replacement of Pi0(x0,y0,z0) Repeating the steps S3-S9 to continue tracing the flow line until the flow line reaches the sink, and then tracing the next flow line.
2. The method for tracking streamline distribution in three-dimensional unstructured mesh flow field by FEM according to claim 1, wherein the tetrahedral mesh vertex coordinates (x) in the step S41,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Performing mutual conversion between a real space and a Master Element space;
the coordinate conversion formula from the Master Element space to the real space is as follows:
x＝x1(1-ξ-η-ζ)+x2ξ+x3η+x4ζ
y＝y1(1-ξ-η-ζ)+y2ξ+y3η+y4ζ
z＝z1(1-ξ-η-ζ)+z2ξ+z3η+z4ζ
the coordinate conversion formula from the real space to the Master Element space is as follows:
wherein
A1＝(x1-x2),A2＝(x1-x3),A3＝(x1-x4)
B1＝(y1-y2),B2＝(y1-y2),B3＝(y1-y4)
C1＝(z1-z2),C2＝(z1-z3),C3＝(z1-z4)
In the formula: x is the x coordinate of real space; y is the y coordinate of real space; z is the z coordinate of real space; xi is a coordinate corresponding to the x coordinate in the Master Element; eta is a coordinate corresponding to the y coordinate in the Master Element; zeta is the coordinate corresponding to the z coordinate in the Master Element; (x)1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Respectively the coordinates of four vertexes in a tetrahedron in the real space.
3. The method for tracking streamline distribution in three-dimensional unstructured mesh flow field by FEM according to claim 2, wherein the velocities [ v ] corresponding to the vertices of the tetrahedral mesh in the step S4x1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]Converting the real space into the Master Element space, wherein the conversion formula is as follows:
wherein the Jacobian matrix is:
in the formula: upsilon isx、υyAnd upsilonzVelocity in the x, y and z directions, respectively, in real space, uξ、υηAnd upsilonζThe velocities in the ξ, η and ζ directions in the Master Element.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010035803.6A CN111210522B (en) | 2020-01-14 | 2020-01-14 | Method for tracking streamline distribution in three-dimensional unstructured grid flow field by using FEM (finite element modeling) |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010035803.6A CN111210522B (en) | 2020-01-14 | 2020-01-14 | Method for tracking streamline distribution in three-dimensional unstructured grid flow field by using FEM (finite element modeling) |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111210522A CN111210522A (en) | 2020-05-29 |
| CN111210522B true CN111210522B (en) | 2021-03-19 |
Family
ID=70786736
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010035803.6A Active CN111210522B (en) | 2020-01-14 | 2020-01-14 | Method for tracking streamline distribution in three-dimensional unstructured grid flow field by using FEM (finite element modeling) |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111210522B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116861753B (en) * | 2023-07-28 | 2024-05-03 | 长江大学 | Novel oil-water two-phase streamline simulation method based on finite difference simulation method |
Family Cites Families (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102156779B (en) * | 2011-04-13 | 2013-03-20 | 北京石油化工学院 | Subsurface flow simulating and predictive analysis method |
| CN102323965A (en) * | 2011-08-25 | 2012-01-18 | 浙江大学 | Intelligent analysis method oriented to unsteady three-dimensional flow field vortex structure |
| CN102495427B (en) * | 2011-12-10 | 2013-06-19 | 北京航空航天大学 | Interface perception ray tracing method based on implicit model expression |
| US9071150B2 (en) * | 2013-05-07 | 2015-06-30 | University Of Central Florida Research Foundation, Inc. | Variable frequency iteration MPPT for resonant power converters |
| CN105069826B (en) * | 2015-08-26 | 2018-06-29 | 中国科学院深圳先进技术研究院 | The modeling method of elastomeric objects amoeboid movement |
| EP3333738A1 (en) * | 2016-12-06 | 2018-06-13 | Fujitsu Limited | Streakline visualization apparatus, method, and program |
| KR20180128804A (en) * | 2017-05-24 | 2018-12-04 | 필드지 주식회사 | Ship including vortex generating fin |
| CN108241777B (en) * | 2017-12-27 | 2020-03-27 | 青岛海洋地质研究所 | Method for calculating seepage velocity field in hydrate sediment based on unstructured grid finite element method |
| CN108875140B (en) * | 2018-05-24 | 2022-06-03 | 西安石油大学 | Method for simulating asphaltene deposition adsorption damage of heavy oil reservoir based on digital core model |
| CN108846224B (en) * | 2018-06-27 | 2019-07-12 | 中国人民解放军国防科技大学 | Supersonic flow channel design method and device |
| CN109635353B (en) * | 2018-11-17 | 2023-04-28 | 天津大学 | Evaluation method for influence of chef operation on lampblack trapping efficiency of range hood |
-
2020
- 2020-01-14 CN CN202010035803.6A patent/CN111210522B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN111210522A (en) | 2020-05-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111709171B (en) | Isogeometric solving and heat dissipation topology generation method for heat flow strong coupling problem | |
| US10296672B2 (en) | Generating inviscid and viscous fluid-flow simulations over an aircraft surface using a fluid-flow mesh | |
| Rane et al. | Grid deformation strategies for CFD analysis of screw compressors | |
| CN114077802B (en) | Particle modeling method using shape function interpolation to replace kernel function approximation | |
| CN105718634A (en) | Airfoil robust optimization design method based on non-probability interval analysis model | |
| CN111210522B (en) | Method for tracking streamline distribution in three-dimensional unstructured grid flow field by using FEM (finite element modeling) | |
| CA2964250A1 (en) | Junction models for simulating proppant transport in dynamic fracture networks | |
| Sitaraman et al. | Progress in strand mesh generation and domain connectivity for dual-mesh cfd simulations | |
| CN107341322A (en) | A kind of method for monitoring core level equipment and pipeline fatigue damage on-line | |
| Pandya et al. | Accuracy, Scalability, and Efficiency of Mixed-Element USM3D for Benchmark Three-Dimensional Flows | |
| Sergio et al. | Optimization methodology assessment for the inlet velocity profile of a hydraulic turbine draft tube: part II—performance evaluation of draft tube model | |
| Cummings et al. | Supersonic, turbulent flow computation and drag optimization for axisymmetric afterbodies | |
| US9087165B2 (en) | Automatic extremum detection on a surface mesh of a component | |
| Mei et al. | Implicit numerical simulation of transonic flow through turbine cascades on unstructured grids | |
| CN110717295B (en) | Method for tracking streamline distribution of tight sandstone reservoir by using finite element method | |
| CN104992046A (en) | Computing system and method of fluid mechanics | |
| CN106294913A (en) | The method improving parts CALCULATION OF THERMAL result reliability | |
| Saleel et al. | On simulation of backward facing step flow using immersed boundary method | |
| Stern et al. | Effects of waves on the wake of a surface-piercing flat plate: experiment and theory | |
| CN106066912A (en) | A kind of generation method of multi partition structured grid | |
| Lakshminarayan et al. | Fully Automated Surface Mesh Adaptation in Strand Grid Framework | |
| CN111008492B (en) | High-order unit Euler equation numerical simulation method based on Jacobian-free matrix | |
| Sharbatdar | Error estimation and mesh adaptation paradigm for unstructured mesh finite volume methods | |
| CN114707384B (en) | Fusion reactor cladding unified three-dimensional model unstructured grid nuclear thermal coupling calculation method | |
| CN114036815B (en) | True and false dual particle model modeling method for coupling physical field quick solving |
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 |