Transient diffusion simulation method and system based on tetrahedron unit

Abstract

The invention provides a transient diffusion simulation method and system based on tetrahedron units, and belongs to the technical field of transient diffusion simulation, and the method comprises the steps: selecting a simulation region, and dividing the simulation region into tetrahedron grids; setting initial conditions; generating a unit matrix; generating an overall matrix; advancing time; substituting the concentration at the previous moment into the overall matrix, substituting the generation item into the overall matrix, generating an equivalent node load, and generating a simultaneous equation set; importing boundary conditions, and solving the simultaneous equation set to obtain a concentration Cn; and judging whether to end or not according to time, and if not, continuing to advance the time until the end. According to the method, aiming at the tetrahedral element, the problem that when the Galerkin finite element method is investigated by diffusion flux with second-order precision, inconsistency exists between a system considering in and out and a system considering transient terms and term generation is solved, so that the dependency of a grid mode is improved, and the precision of a simulation solution obtained by a low-quality element is improved.

Description

Technical Field

[0001] The invention relates to the technical field of transient diffusion simulation, and in particular to a transient diffusion simulation method and system based on tetrahedron units. Background Art

[0002] The Galerkin finite element method is a widely used numerical solution method for differential equations, with its clear error evaluation and wide application range. However, its problems such as mesh pattern dependence and increased simulation error due to low element quality have not yet been solved.

[0003] When considering the inflow and outflow of diffusion flow, the Galerkin finite element method uses a second-order accuracy system that changes due to the unit shape. However, when considering transient terms and source terms, it uses another system that is nearly 1 / 4 of the unit volume, which cannot satisfy the conservation law with second-order accuracy. This system difference problem is called inconsistency. For example, for Figure 2 As for node 1 of the tetrahedral element, as shown in Table 1, the flow area of ​​the Galerkin finite element method when considering the inflow and outflow caused by diffusion is the sum of the diffusion inflow and outflow areas of the other three nodes, but when considering the transient term and the generated term, only the contribution from nearly 1 / 4 of the unit volume is obtained. Similarly, there are inconsistencies for other nodes. By the way, this inconsistency problem does not exist for regular tetrahedral elements. When the diffusion flux with a lower first-order accuracy is used to consider the diffusion inflow and outflow area, this inconsistency problem does not exist for elements of other shapes.

[0004] However, there are no effective improvement measures for the inconsistency problem in the Galerkin finite element method from the perspective of second-order accuracy. Summary of the invention

[0005] The present invention provides a transient diffusion simulation method and system based on tetrahedral units. For tetrahedral units, the inconsistency problem between the system considering entry and exit and the system considering transient terms and generation terms generated when the Galerkin finite element method is examined with diffusion flux of second-order accuracy is solved, thereby improving the dependence of the grid mode and improving the accuracy of the simulation solution obtained by low-quality units.

[0006] This specification discloses a transient diffusion simulation method based on tetrahedral units, including:

[0007] Select a simulation area and divide the simulation area into tetrahedral grids;

[0008] Set the initial conditions, t 0 =0, n=1;

[0009] Generate a cell matrix:

[0010]

[0011] Where C represents concentration, Δt is the time step, t represents time, and C n represents the concentration at the nth moment, C n-1 represents the concentration at the n-1th moment, k represents the diffusion coefficient, ∫ e de represents the integral over element e, where e represents the tetrahedral element around node P, and W P represents the Galerkin weight function, the subscript O represents the value at the geometric center of the unit, Vn P is the volume of the virtual node domain of the object node P in the unit, V represents the volume of the unit, is the time parameter, represents the partial derivative of x, represents the partial derivative of y, It means partial derivative of z, x, y, z represent three-dimensional rectangular coordinates;

[0012] Generate the overall matrix:

[0013]

[0014] Among them, Σ represents the sum of the integrals of each unit;

[0015] Advance time: t n =t n-1 +△t;t n represents the nth moment, t n-1 represents the n-1th moment;

[0016] Substituting the concentration at the previous moment into Substitute the generated term S into Generate equivalent nodal loads and generate a system of simultaneous equations:

[0017]

[0018]

[0019] Among them, S n is the generated term at the nth moment, S n-1 is the generated item at the n-1th moment;

[0020] Import the boundary conditions and solve the simultaneous equations to obtain the concentration C n ;

[0021] Use time to determine whether it is finished. If it is not finished, then n=n+1 and continue to advance the time until it is finished.

[0022] In this specification, the process of generating the simultaneous equations is as follows:

[0023] Applying the Galerkin finite element method to equation (1), we obtain equation (2) for the object node P considered;

[0024]

[0025] is a transient term, Ω represents the area around the node P consisting of tetrahedral units, ∫ Ω dΩ represents the integral over the region Ω, W P is the Galerkin weight function; the integral on Ω is the sum of the integrals of the surrounding units, and equation (2) can be rewritten as equation (3);

[0026]

[0027] Applying the integration by parts theorem and Gauss-Green theorem to the unit integration on the left side of equation (3), we obtain equation (4) which is the discretization of the Galerkin finite element method;

[0028]

[0029] Use equation (12) to replace equation (4) of the Galerkin method;

[0030]

[0031] in, It is a correction term introduced into the Galerkin finite element method;

[0032] For a unit of any shape, the sum of the volumes of the virtual node domains of all nodes of the same unit satisfies equation (13);

[0033]

[0034] It shows that conservation laws are satisfied at the unit level;

[0035] Formula (12) can be rewritten as formula (19) after moving the terms;

[0036]

[0037] Formula (19) is valid at any time, so we get formula (20);

[0038]

[0039] Substituting equation (19) into equation (20), we get equation (21);

[0040]

[0041] Substituting the mean value theorem of equation (22) into equation (21), the partial derivative with respect to time vanishes, and the finite element discretization result of the transient diffusion equation is obtained, that is, a system of simultaneous equations is obtained.

[0042] In this manual, Figure 2 The three faces shown intersect vertically at node 1, the coordinate origin is placed at node 1, and nodes 2, 3, and 4 are placed on the positive x, y, and z axes, respectively. The unit integral on the left side of equation (4) is calculated as equation (5);

[0043]

[0044] The numerical subscript indicates the node number in the unit; let the node number of the object node P in the unit be 1; ij represents the midpoint of the line segment connecting nodes i and j; k(C1-C2) / x2 in the first term on the right side of equation (5) is the midpoint I 12 The diffusion flux from node 1 to node 2 is approximated by the second-order accuracy, and its coefficient y3z4 /

[0045] 6 is the flow area of ​​the diffusion flux, which passes through point I 12 and perpendicular to the diffusion flux; the physical meaning of the first term on the right side of equation (5) is interpreted as the use point I 12 The second-order accuracy diffusion flux at the point is used to evaluate the diffusion flow rate of the flow area y3z4 / 6; the second term on the right side of equation (5) uses the midpoint I 13 The second-order diffusion flux k(C1-C3) / y3 from node 1 to node 3 is used to evaluate the diffusion flow of the flow area x2z4 / 6; the area x2z4 / 6 passes through point I 13 and is perpendicular to the diffusion flux k(C1-C3) / y3; the third term on the right side of equation (5) uses the midpoint I 14 The second-order accuracy diffusion flux k(C1-C4) from node 1 to node 4 /

[0046] z4 is used to evaluate the diffusion flow of the flow area x2y3 / 6; the area x2y3 / 6 passes through point I 14 And it is perpendicular to the diffusion flux k(C1-C4) / z4; Assume that the node numbers of the object node P in the unit are 2, 3, and 4 respectively, and the unit integrals on the left side of equation (4) are calculated as equations (34), (35), and (36) respectively;

[0047]

[0048] The flow areas of the second-order diffusion flux from nodes 2, 3, and 4 to node 1 are y3z4 / 6, x2z4 / 6, and x2y3 / 6, respectively. The flow area of ​​the second-order diffusion flux pointing from nodes 2, 3, and 4 to each other is zero.

[0049] In this specification, the element integral on the right side of the tetrahedral element (4) of arbitrary shape is calculated, and the contributions from the transient term and the generated term are almost equally distributed to the four nodes of the element in the form of equivalent node loads; when When and S are constants inside the unit, the unit integral is written as (6);

[0050]

[0051] For node 1, the flow area of ​​the diffusion flux is the sum of the flow areas of the diffusion flux of the other three nodes, and the transient term and generated term obtained only come from nearly 1 / 4 of the unit volume.

[0052] In this specification, for a tetrahedral unit of any shape, the node number of the object node P in the unit is set to 1, and the unit integration result on the left side of equation (4) is written as equation (8);

[0053]

[0054]

[0055] A 1i is a coefficient that depends only on the coordinates of the element nodes. Its first subscript 1 is due to the fact that the node number of the object node P in the unit is 1. The representation of the second line of formula (8) uses the theoretically proven formula (9).

[0056]

[0057] In formula (8), L is introduced 1i (i=2,3,4), let L 1i =[(x i -x1) 2 +(y i -y1) 2 +(z i -z1) 2 ] 1 / 2 , that is, it is defined as the distance between node 1 and node i; k(C1-C i ) / L 1i is the midpoint I of the line segment connecting nodes 1 and i 1i The diffusion flux from node 1 to i at the second-order accuracy, its coefficient (-A 1i L 1i ) is the area; the area passes through the midpoint I 1i and the diffusion flux k(C1-C i ) / L 1i The two are multiplied, and the right side of equation (8) becomes the second-order accuracy evaluation of the flow area (-A 1i L1i ) diffusion flow; Figure 3 Shown as area (-A 1i L 1i ) as the bottom surface, node 1 as the vertex, and L 1i / 2 is the height of the three cones, and their volumes are Vn 1i To express, it is defined by formula (10);

[0058]

[0059] Vn P =Vn1=Vn 12 +Vn 13 +Vn 14 (11);

[0060] The volume of the virtual node domain Vn1 of node 1 is the sum of the volumes of the three cones. Since node P and node 1 are identical, Vn P It is also calculated by formula (11).

[0061] In this specification, for a regular tetrahedral unit, Vn P =V / 4, correction term is zero; for units of other shapes, the unit integral on the right side of equation (4), that is, the first term in the brackets on the right side of equation (12) The contribution to the object node is considered to come from the unit volume V being P The weighted transient term and generated term; W P The unit integral itself is equal to V / 4. The correction term introduced by the Galerkin finite element method is the volume of the virtual node domain Vn P After V / 4 is subtracted from the equation, it is multiplied by the transient term and generated term at the geometric center. When the transient term and generated term are constants in the unit, the square brackets on the right side of equation (12) become That is, the volume of the system after correction considering transient terms and generated terms is equal to the volume of the system when the inlet and outlet are considered with the diffusion flow of second-order accuracy.

[0062] In this specification, based on eight shapes, the three faces passing through the node 1 are perpendicular to each other, and the vertical intersecting faces have a common side length of h and a volume of V. r The units represented by Figure 5a The grid is integrated by Galerkin's method (4) to obtain (14);

[0063]

[0064] Assume that the transient term and the generated term are constants in the grid; the subscript A represents the value of C at node A, where node A is the center point of the grid and becomes the object node P, corresponding to node 1; the subscripts -x, +x, -y, +y, -z, +z represent the value of C at the position where point A moves -h, h in the x direction, -h, h in the y direction, and -h, h in the z direction, respectively; the left side of equation (14) is the second-order accuracy discretization of the second-order differential term in equation (1) multiplied by 4V r , while the transient term and generated term on the right side are the corresponding terms in equation (1) multiplied by 2V r , does not converge to the original differential equation; from equations (10) and (11), we can get the virtual node domain Vn of node A in each unit: P =V r / 2, as a result, the right side of equation (12) becomes equation (15);

[0065]

[0066] It is confirmed that equation (16) obtained from equation (12) obviously converges to the original differential equation;

[0067] Based on 24 shapes whose three faces passing through node 1 are perpendicular to each other, the perpendicular intersecting faces have a total side length of h and a volume of V r The unit represented by and the area enclosed by 8 regular tetrahedral units Figure 5b The grid is integrated by Galerkin's method (4) to obtain (17);

[0068]

[0069] Subscript B represents the value of C at node B, where node B is the center point of the grid and becomes the object node P; subscripts -x, +x, -y, +y, -z, +z represent the value of C at the position where point B moves -h, h in the x direction, -h, h in the y direction, and -h, h in the z direction, respectively; the first term on the left side of equation (17) comes from the calculation of 24 volumes V r The second term on the left side of formula (17) comes from the sum of the integrals of the eight regular tetrahedral units, and uses the Taylor expansion formula; the volume of the regular tetrahedral unit is equal to 2V r The sum of the two terms on the left is equivalent to the second-order differential term in equation (1) multiplied by 8V r The transient term and generated term on the right side of equation (17) are the corresponding terms in equation (1) multiplied by 10V. r , does not converge to the original differential equation; from equations (10) and (11), the right side of equation (12) about node B is obtained as Among them, 24×V r / 6 is the node B in 24 with volume Vr The sum of the virtual node domains of the corresponding nodes in the unit, 8×(2V r ) / 4 is the sum of the virtual node domains of the nodes corresponding to node B in the 8 regular tetrahedral units; the sum of the two is As a result, equation (18) converges to the original differential equation from equation (12);

[0070]

[0071] In this specification, when the simulation domain is a cubic region, 0≤x≤1, 0≤y≤1, 0≤z≤1, k=1; for the numerical error in the forced boundary condition problem, the generated term is taken as equation (24), the initial condition is taken as equation (25), and the first type of boundary condition is taken as equation (26), and the analytical solution is obtained as equation (27);

[0072]

[0073] C 0 =sin(πx)sin(πy)sin(πz)(t=0) (25);

[0074] C=0(x=0,1, or y=0,1, or z=0,1) (26);

[0075]

[0076] The numerical error is evaluated using formula (33);

[0077]

[0078] Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the simulation solution of node i, and the maximum value of the analytical solution in the simulation area, respectively; the superscript end represents the end time of the calculation, the numerical simulation is performed with double precision, and the unit integration is performed through four-point numerical integration; End time t end =1.

[0079] In this specification, when the simulation domain is a cubic region, 0≤x≤1, 0≤y≤1, 0≤z≤1, k=1; for the numerical error in the mixed boundary condition problem, the generated term is taken as equation (28), the initial condition is taken as equation (29), the first type of boundary condition is taken as equation (30), and the second type of boundary condition is taken as equation (31), and the analytical solution is obtained as equation (32);

[0080]

[0081]

[0082] In formula (31), n ​​is the external normal of unit length;

[0083]

[0084] The numerical error is evaluated using formula (33);

[0085]

[0086] Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the simulation solution of node i, and the maximum value of the analytical solution in the simulation area, respectively; the superscript end represents the end time of the calculation, the numerical simulation is performed with double precision, and the unit integration is performed through four-point numerical integration; End time t end =1.

[0087] This specification also discloses a transient diffusion simulation system based on tetrahedron units, which is used to implement any of the above-mentioned transient diffusion simulation methods based on tetrahedron units. The transient diffusion simulation system based on tetrahedron units includes:

[0088] A partitioning module is used to select a simulation area and divide the simulation area into tetrahedral grids;

[0089] Setting module, used to set initial conditions, t 0 =0, n=1;

[0090] Unit matrix generation module, used to generate unit matrix:

[0091]

[0092] Where C represents concentration, Δt is the time step, t represents time, and C n represents the concentration at the nth moment, C n-1 represents the concentration at the n-1th moment, k represents the diffusion coefficient, ∫ e de represents the integral over element e, where e represents the tetrahedral element around node P, and W P represents the Galerkin weight function, the subscript O represents the value at the geometric center of the unit, Vn P is the volume of the virtual node domain of the object node P in the unit, V represents the volume of the unit, is the time parameter, represents the partial derivative of x, represents the partial derivative of y, It means partial derivative of z, x, y, z represent three-dimensional rectangular coordinates;

[0093] The overall matrix generation module is used to generate the overall matrix:

[0094]

[0095] Among them, Σ represents the sum of the integrals of each unit;

[0096] Time advancement module, used to advance time: t n =t n-1 +△t;t n represents the nth moment, t n-1 represents the n-1th moment;

[0097] The simultaneous equations generation module is used to substitute the concentration at the previous moment into Substitute the generated term S into Generate equivalent nodal loads and generate a system of simultaneous equations:

[0098]

[0099] Among them, S n is the generated term at the nth moment, S n-1 is the generated item at the n-1th moment;

[0100] The simultaneous equations solving module is used to import boundary conditions, solve the simultaneous equations, and obtain the concentration C n ;

[0101] The end module is used to determine whether it is ended based on time. If it is not ended, n=n+1 and the time continues to advance until it is ended.

[0102] The embodiments of this specification can at least achieve the following beneficial effects:

[0103] The present invention uses tetrahedral units to perform finite element simulation of transient diffusion problems. By introducing correction terms based on virtual node domains, transient terms and generated terms are corrected to ensure that conservation laws are satisfied at the node level with second-order accuracy. The introduction of the correction terms not only ensures that conservation laws are satisfied at the node level, but also satisfies conservation laws at the unit level, and effectively improves the convergence of finite element discretization to the original differential equation.

[0104] The numerical simulation results using a variety of grid modes show that compared with the traditional Galerkin finite element method, the present invention significantly improves the simulation accuracy and improves the grid mode dependency problem. BRIEF DESCRIPTION OF THE DRAWINGS

[0105] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for describing the embodiments are briefly introduced below.

[0106] Figure 1 It is a schematic diagram of the transient diffusion simulation method based on tetrahedral units involved in the present invention.

[0107] Figure 2It is a schematic diagram of the diffusion flow area of ​​the node 1 of the tetrahedron with three perpendicular faces involved in the present invention.

[0108] Figure 3 It is a schematic diagram of the virtual node domain of node 1 involved in the present invention.

[0109] Figure 4 It is a schematic diagram of the diffusion flow area of ​​the node 1 of the regular tetrahedron involved in the present invention.

[0110] Figure 5a Schematic diagram of node A and its surrounding units involved in the present invention.

[0111] Figure 5b Schematic diagram of node B and its surrounding units involved in the present invention.

[0112] Figure 6a It is a schematic diagram of a cube composed of five units involved in the present invention.

[0113] Figure 6b It is a schematic diagram of a cube composed of 6 units involved in the present invention.

[0114] Figure 6c It is a schematic diagram of a rectangular parallelepiped composed of 14 units involved in the present invention.

[0115] Figure 6d It is a schematic diagram of a tetrahedron involved in the present invention divided into four units.

[0116] Figure 6e It is a schematic diagram of the irregular unit division involved in the present invention.

[0117] Figure 7 Schematic diagram of the numerical error of the forced boundary condition problem involved in the present invention.

[0118] Figure 8 Schematic diagram of the numerical error of the mixed boundary condition problem involved in the present invention. DETAILED DESCRIPTION

[0119] The embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0120] like Figure 1 As shown, this embodiment provides a transient diffusion simulation method based on tetrahedral units, including:

[0121] Select a simulation area and divide the simulation area into tetrahedral grids;

[0122] Set the initial conditions, t 0 =0, n=1;

[0123] Generate a cell matrix:

[0124]

[0125]

[0126] Where C represents concentration, Δt is the time step, t represents time, and C n represents the concentration at the nth moment, C n-1 represents the concentration at the n-1th moment, k represents the diffusion coefficient, ∫ e de represents the integral over element e, where e represents the tetrahedral element around node P, and W P represents the Galerkin weight function, the subscript O represents the value at the geometric center of the unit, Vn P is the volume of the virtual node domain of the object node P in the unit, V represents the volume of the unit, is the time parameter, represents the partial derivative of x, represents the partial derivative of y, It means partial derivative of z, x, y, z represent three-dimensional rectangular coordinates;

[0127] Generate the overall matrix:

[0128]

[0129] Among them, Σ represents the sum of the integrals of each unit;

[0130] Advance time: t n =t n-1 +△t;t n represents the nth moment, t n-1 represents the n-1th moment;

[0131] Substituting the concentration at the previous moment into Substitute the generated term S into Generate equivalent nodal loads and generate a system of simultaneous equations:

[0132]

[0133] Among them, S n is the generated term at the nth moment, S n-1 is the generated item at the n-1th moment;

[0134] Import the boundary conditions and solve the simultaneous equations to obtain the concentration C n ;

[0135] Use time to determine whether it is finished. If it is not finished, then n=n+1 and continue to advance the time until it is finished.

[0136] In some embodiments, the process of generating the set of simultaneous equations is as follows:

[0137] Applying the Galerkin finite element method to equation (1), we obtain equation (2) for the object node P considered;

[0138]

[0139] is a transient term, Ω represents the area around the node P consisting of tetrahedral units, ∫ Ω dΩ represents the integral over the region Ω, W P is the Galerkin weight function; the integral on Ω is the sum of the integrals of the surrounding units, and equation (2) can be rewritten as equation (3);

[0140]

[0141] Applying the integration by parts theorem and Gauss-Green theorem to the unit integration on the left side of equation (3), we obtain equation (4) which is the discretization of the Galerkin finite element method;

[0142]

[0143] Use equation (12) to replace equation (4) of the Galerkin method;

[0144]

[0145] in, It is a correction term introduced into the Galerkin finite element method;

[0146] For a unit of any shape, the sum of the volumes of the virtual node domains of all nodes of the same unit satisfies equation (13);

[0147]

[0148] It shows that conservation laws are satisfied at the unit level;

[0149] Formula (12) can be rewritten as formula (19) after moving the terms;

[0150]

[0151] Formula (19) is valid at any time, so we get formula (20);

[0152]

[0153] Substituting equation (19) into equation (20), we get equation (21);

[0154]

[0155] Substituting the mean value theorem of equation (22) into equation (21), the partial derivative with respect to time vanishes, and the finite element discretization result of the transient diffusion equation is obtained, that is, a system of simultaneous equations is obtained.

[0156] In some embodiments, for Figure 2 The three faces shown intersect vertically at node 1, the coordinate origin is placed at node 1, and nodes 2, 3, and 4 are placed on the positive x, y, and z axes, respectively. The unit integral on the left side of equation (4) is calculated as equation (5);

[0157]

[0158] The numerical subscript indicates the node number in the unit; let the node number of the object node P in the unit be 1; ij represents the midpoint of the line segment connecting nodes i and j; k(C1-C2) / x2 in the first term on the right side of equation (5) is the midpoint I 12 The second-order accuracy approximation of the diffusion flux from node 1 to node 2 is given by the coefficient y3z4 / 6, which is the flow area of ​​the diffusion flux. This area passes through point I 12 and perpendicular to the diffusion flux; the physical meaning of the first term on the right side of equation (5) is interpreted as the use point I 12 The second-order accuracy diffusion flux at the point is used to evaluate the diffusion flow rate of the flow area y3z4 / 6; the second term on the right side of equation (5) uses the midpoint I 13 The second-order diffusion flux k(C1-C3) / y3 from node 1 to node 3 is used to evaluate the diffusion flow of the flow area x2z4 / 6; the area x2z4 / 6 passes through point I 13 and is perpendicular to the diffusion flux k(C1-C3) / y3; the third term on the right side of equation (5) uses the midpoint I 14 The second-order diffusion flux k(C1-

[0159] C4) / z4 to evaluate the diffusion flow of the flow area x2y3 / 6; the area x2y3 / 6 passes through point I 14 And it is perpendicular to the diffusion flux k(C1-C4) / z4; Assume that the node numbers of the object node P in the unit are 2, 3, and 4 respectively, and the unit integrals on the left side of equation (4) are calculated as equations (34), (35), and (36) respectively;

[0160]

[0161] The flow areas of the second-order diffusion flux from nodes 2, 3, and 4 to node 1 are y3z4 / 6, x2z4 / 6, and x2y3 / 6, respectively. The flow area of ​​the second-order diffusion flux pointing from nodes 2, 3, and 4 to each other is zero.

[0162] In some embodiments, for a tetrahedral element of arbitrary shape, the element integral on the right side of equation (4) is calculated, and the contributions from the transient term and the generated term are distributed almost equally to the four nodes of the element in the form of equivalent node loads; when When and S are constants inside the unit, the unit integral is written as (6);

[0163]

[0164] For node 1, the flow area of ​​the diffusion flux is the sum of the flow areas of the diffusion flux of the other three nodes, and the transient term and generated term obtained only come from nearly 1 / 4 of the unit volume.

[0165] In some embodiments, for a tetrahedral unit of any shape, the node number of the object node P in the unit is set to 1, and the unit integration result on the left side of equation (4) is written as equation (8);

[0166]

[0167] A 1i is a coefficient that depends only on the coordinates of the element nodes. Its first subscript 1 is due to the fact that the node number of the object node P in the unit is 1. The representation of the second line of formula (8) uses the theoretically proven formula (9).

[0168]

[0169] In formula (8), L is introduced 1i (i=2,3,4), let L 1i =[(x i -x1) 2 +(y i -y1) 2 +(z i -z1) 2 ] 1 / 2 , that is, it is defined as the distance between node 1 and node i; k(C1-C i ) / L 1i is the midpoint I of the line segment connecting nodes 1 and i 1i The diffusion flux from node 1 to i at the second-order accuracy, its coefficient (-A 1i L 1i ) is the area; the area passes through the midpoint I 1i and the diffusion flux k(C1-C i ) / L 1i The two are multiplied, and the right side of equation (8) becomes the second-order accuracy evaluation of the flow area (-A 1i L 1i ) diffusion flow; Figure 3 Shown as area (-A1i L 1i ) as the bottom surface, node 1 as the vertex, and L 1i / 2 is the height of the three cones, and their volumes are Vn 1i To express, it is defined by formula (10);

[0170]

[0171] Vn P =Vn1=Vn 12 +Vn 13 +Vn 14 (11);

[0172] The volume of the virtual node domain Vn1 of node 1 is the sum of the volumes of the three cones. Since node P and node 1 are identical, Vn P It is also calculated by formula (11).

[0173] In some embodiments, for a regular tetrahedral unit, Vn P =V / 4, correction term is zero; for units of other shapes, the unit integral on the right side of equation (4), that is, the first term in the brackets on the right side of equation (12) The contribution to the object node is considered to come from the unit volume V being P The weighted transient term and generated term; W P The unit integral itself is equal to V / 4. The correction term introduced by the Galerkin finite element method is the volume of the virtual node domain Vn P After V / 4 is subtracted from the equation, it is multiplied by the transient term and generated term at the geometric center. When the transient term and generated term are constants in the unit, the square brackets on the right side of equation (12) become That is, the volume of the system after correction considering transient terms and generated terms is equal to the volume of the system when the inlet and outlet are considered with the diffusion flow of second-order accuracy.

[0174] In some embodiments, based on eight shapes in which three faces passing through node 1 are perpendicular to each other, the perpendicular intersecting faces have a common side length of h and a volume of V. r The units represented by Figure 5a The grid is integrated by Galerkin's method (4) to obtain (14);

[0175]

[0176] Assume that the transient term and the generated term are constants in the grid; the subscript A represents the value of C at node A, where node A is the center point of the grid and becomes the object node P, corresponding to node 1; the subscripts -x, +x, -y, +y, -z, +z represent the value of C at the position where point A moves -h, h in the x direction, -h, h in the y direction, and -h, h in the z direction, respectively; the left side of equation (14) is the second-order accuracy discretization of the second-order differential term in equation (1) multiplied by 4V r , while the transient term and generated term on the right side are the corresponding terms in equation (1) multiplied by 2V r , does not converge to the original differential equation; from equations (10) and (11), we can get the virtual node domain Vn of node A in each unit: P =V r / 2, as a result, the right side of equation (12) becomes equation (15);

[0177]

[0178] It is confirmed that equation (16) obtained from equation (12) obviously converges to the original differential equation;

[0179] Based on 24 shapes whose three faces passing through node 1 are perpendicular to each other, the perpendicular intersecting faces have a total side length of h and a volume of V r The unit represented by and the area enclosed by 8 regular tetrahedral units Figure 5b The grid is integrated by Galerkin's method (4) to obtain (17);

[0180]

[0181] Subscript B represents the value of C at node B, where node B is the center point of the grid and becomes the object node P; subscripts -x, +x, -y, +y, -z, +z represent the value of C at the position where point B moves -h, h in the x direction, -h, h in the y direction, and -h, h in the z direction, respectively; the first term on the left side of equation (17) comes from the calculation of 24 volumes V r The second term on the left side of formula (17) comes from the sum of the integrals of the eight regular tetrahedral units, and uses the Taylor expansion formula; the volume of the regular tetrahedral unit is equal to 2V r The sum of the two terms on the left is equivalent to the second-order differential term in equation (1) multiplied by 8V r The transient term and generated term on the right side of equation (17) are the corresponding terms in equation (1) multiplied by 10V. r , does not converge to the original differential equation; from equations (10) and (11), the right side of equation (12) about node B is obtained as Among them, 24×V r / 6 is the node B in 24 with volume Vr The sum of the virtual node domains of the corresponding nodes in the unit, 8×(2V r ) / 4 is the sum of the virtual node domains of the nodes corresponding to node B in the 8 regular tetrahedral units; the sum of the two is As a result, equation (18) converges to the original differential equation from equation (12);

[0182]

[0183] In some embodiments, when the simulation domain is a cubic region, 0≤x≤1, 0≤y≤1, 0≤z≤1, k=1; for the numerical error in the forced boundary condition problem, the generated term is taken as equation (24), the initial condition is taken as equation (25), and the first type of boundary condition is taken as equation (26), and the analytical solution is obtained as equation (27);

[0184]

[0185] C 0 =sin(πx)sin(πy)sin(πz)(t=0) (25);

[0186] C=0(x=0,1, or y=0,1, or z=0,1) (26);

[0187]

[0188] The numerical error is evaluated using formula (33);

[0189]

[0190] Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the simulation solution of node i, and the maximum value of the analytical solution in the simulation area, respectively; the superscript end represents the end time of the calculation, the numerical simulation is performed with double precision, and the unit integration is performed through four-point numerical integration; End time t end =1.

[0191] In some embodiments, when the simulation domain is a cubic region, 0≤x≤1, 0≤y≤1, 0≤z≤1, k=1; for the numerical error in the mixed boundary condition problem, the generated term is taken as equation (28), the initial condition is taken as equation (29), the first type of boundary condition is taken as equation (30), and the second type of boundary condition is taken as equation (31), and the analytical solution is obtained as equation (32);

[0192]

[0193]

[0194] In formula (31), n ​​is the external normal of unit length;

[0195]

[0196] The numerical error is evaluated using formula (33);

[0197]

[0198] Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the simulation solution of node i, and the maximum value of the analytical solution in the simulation area, respectively; the superscript end represents the end time of the calculation, the numerical simulation is performed with double precision, and the unit integration is performed through four-point numerical integration; End time t end =1.

[0199] The embodiment of this specification also discloses a transient diffusion simulation system based on tetrahedron units, which is used to implement any of the above-mentioned transient diffusion simulation methods based on tetrahedron units. The transient diffusion simulation system based on tetrahedron units includes:

[0200] A partitioning module is used to select a simulation area and divide the simulation area into tetrahedral grids;

[0201] Setting module, used to set initial conditions, t 0 =0, n=1;

[0202] A unit matrix generation module, used to generate a unit matrix;

[0203] An overall matrix generation module, used to generate an overall matrix;

[0204] Time advancement module, used to advance time: t n =t n-1 +△t;

[0205] A simultaneous equations generation module is used to substitute the concentration at the previous moment into the overall matrix, substitute the generated term S into the overall matrix, generate equivalent node loads, and generate simultaneous equations;

[0206] The simultaneous equations solving module is used to import boundary conditions, solve the simultaneous equations, and obtain the concentration C n ;

[0207] The end module is used to determine whether it is ended based on time. If it is not ended, n=n+1 and the time continues to advance until it is ended.

[0208] In order to further illustrate the technical solution of the present invention, the technical concept of the present invention is as follows:

[0209] The governing equation for the transient diffusion problem is shown in equation (1).

[0210] Applying the Galerkin finite element method to equation (1), the integral equation of equation (2) is obtained for the object node P considered.

[0211] The integral over Ω is the sum of the integrals of the surrounding units, so equation (2) can be rewritten as equation (3).

[0212] Applying the integration by parts theorem and Gauss-Green theorem to the unit integration on the left side of equation (3), we obtain equation (4) which is the discretization of the Galerkin finite element method.

[0213] Since this manual does not make any improvements to the second type of boundary conditions, it will not be described again.

[0214] for Figure 2 The three faces shown intersect vertically at node 1, the coordinate origin is placed at node 1, and nodes 2, 3, and 4 are placed on the positive axes of x, y, and z respectively. The unit integral on the left side of equation (4) can be calculated as equation (5).

[0215] Here, the numerical subscript represents the node number in the unit. And let the node number of the object node P in the unit be 1. ij represents the midpoint of the line segment connecting nodes i and j. Therefore, k(C1-C2) / x2 in the first term on the right side of equation (5) is the midpoint I 12 The second-order accuracy approximation of the diffusion flux from node 1 to node 2 is given by the coefficient y3z4 / 6, which is the flow area of ​​the diffusion flux. This area passes through point I 12 and perpendicular to the diffusion flux. The physical meaning of the first term on the right side of equation (5) can be interpreted as the use point I 12 The second-order accuracy diffusion flux at the midpoint I is used to evaluate the diffusion flow rate of the flow area y3z4 / 6. Similarly, the second term on the right side of equation (5) uses the midpoint I 13 The diffusion flux k(C1-C3) / y3 from node 1 to node 3 is used to evaluate the diffusion flow of the flow area x2z4 / 6. The area x2z4 / 6 passes through point I 13 and is perpendicular to the diffusion flux k(C1-C3) / y3. The third term on the right side of equation (5) uses the midpoint I 14 The diffusion flux k(C1-C4) / z4 from node 1 to node 4 is used to evaluate the diffusion flow of the flow area x2y3 / 6. The area x2y3 / 6 passes through point I 14 And it is perpendicular to the diffusion flux k(C1-C4) / z4. The flow surfaces of these three diffusion fluxes do not form a butt joint state, and their areas are as mentioned above, but their shapes cannot be defined. Figure 2 The flow surfaces are shown as schematic representations and do not represent their shapes.

[0216] for Figure 2The three faces shown intersect vertically at node 1, the coordinate origin is placed at node 1, nodes 2, 3, and 4 are placed in the unit on the positive x, y, and z axes respectively, and the node numbers of the object node P in the unit are 2, 3, and 4 respectively. The unit integrals on the left side of equation (4) are calculated as equations (34), (35), and (36) respectively. The flow areas of the second-order diffusion flux from nodes 2, 3, and 4 to node 1 are y3z4 / 6, x2z4 / 6, and x2y3 / 6 respectively, and the flow areas of the second-order diffusion flux pointing to each other between nodes 2, 3, and 4 are zero. The flow areas between nodes are listed in Table 1.

[0217] Table 1 Diffusion flow area between two nodes

[0218]

[0219] On the other hand, for the tetrahedral element of arbitrary shape, the element integral on the right side of equation (4) is calculated, and the contributions from the transient term and the generated term are almost equally distributed to the four nodes of the element in the form of equivalent node loads. In particular, when When and S are constants inside the unit, the unit integral can be simply written as formula (6).

[0220] When observing the same system between two moments in time, the conservation law of physical quantities becomes relation (7).

[0221]

[0222] And for Figure 2 For node 1 of the tetrahedral element, as shown in Table 1, the flow area of ​​the Galerkin finite element method when considering the inflow and outflow caused by diffusion is the sum of the diffusion inflow and outflow areas of the other three nodes. However, when considering transient terms and generated terms, only contributions from nearly 1 / 4 of the unit volume are obtained, and the conservation law cannot be satisfied with second-order accuracy. This system difference problem is called inconsistency. Similarly, there are inconsistency problems for other nodes. By the way, this inconsistency problem does not exist for regular tetrahedral elements. When the diffusion flux with lower first-order accuracy is used to consider the diffusion inflow and outflow area, this inconsistency problem does not exist for elements of other shapes.

[0223] However, there is no effective improvement measure for the inconsistency problem in the Galerkin finite element method from the perspective of second-order accuracy. In simulation practice, in order to reduce the error caused by inconsistency, human and financial resources are invested in meshing to make the unit shape as close to a regular tetrahedron as possible.

[0224] For tetrahedral units, the present invention aims to solve the inconsistency problem between the system considering entry and exit and the system considering transient terms and generation terms when the Galerkin finite element method is examined with diffusion flux of second-order accuracy, thereby improving the dependence of the grid pattern and improving the accuracy of the simulation solution obtained by low-quality units.

[0225] One of the core parts of the present invention: a new concept of virtual node domain. For a tetrahedral unit of any shape, the node number of the object node P in the unit is set to 1, and the unit integration result on the left side of equation (4) can be written as equation (8).

[0226] The representation of the second line of equation (8) utilizes equation (9) which can be theoretically proven.

[0227] In formula (8), L is introduced 1i (i=2,3,4), let L 1i =[(x i -x1) 2 +(y i -y1) 2 +(z i -z1) 2 ] 1 / 2 , that is, it is defined as the distance between node 1 and node i. The k(C1-C i ) / L 1i is the midpoint I of the line segment connecting nodes 1 and i 1i The diffusion flux from node 1 to i at the second-order accuracy, its coefficient (-A 1i L 1i ) is the area. This area passes through the midpoint I 1i and the diffusion flux k(C1-C i ) / L 1i The two are multiplied, and the right side of equation (8) becomes the second-order accuracy evaluation of the flow area (-A 1i L 1i ) diffusion flow. Figure 3 Respectively represent the area (-A 1i L 1i ) as the bottom surface, node 1 as the vertex, and L 1i / 2 is the height of the three cones, and their volumes are Vn 1i It is represented by , which is defined by equation (10). Due to the shape of the unit, the volume of the cone may be negative, but this does not require any special computational treatment and does not affect the numerical simulation to be implemented.

[0228] The volume of the virtual node domain Vn1 of node 1 is the sum of the volumes of the three cones. Since node P and node 1 are identical, Vn P It is also calculated by formula (11).

[0229] For tetrahedrons other than regular tetrahedrons, the flow surfaces of the diffusion flux cannot be connected to each other to form a closed state. Therefore, for the virtual node domain, the size of its volume is defined by equations (10) and (11), but the shape cannot be defined. In this sense, Figure 3 The illustration shows the volume of the virtual node domains rather than their shapes. In the case of regular tetrahedral elements, Figure 4 As shown, the flow surface of the diffusion flux can be butted and closed. In the figure, O represents the geometric center of the tetrahedron. i Represents the geometric center of the surface triangle opposite to node i. For example, O2 is the geometric center of the triangle opposite to node 2 with nodes 1, 4, and 3 as vertices.

[0230] The second core part of the present invention: Modification of transient terms and generated terms. Here, the processing method of transient terms and generated terms is modified, and the Galerkin method (4) is replaced by (12) to match the in-out evaluation system with second-order accuracy.

[0231] Here, the second term on the right side is the correction term added to the Galerkin method (4) proposed by the present invention, where the subscript O represents the value at the geometric center of the unit, and Vn P It is the virtual node domain of the object node P in the unit obtained by formula (11).

[0232] The meaning of the correction term in formula (12)

[0233] Consistency: For regular tetrahedral units Vn P =v / 4, the correction term is zero, because this type of unit itself does not have the inconsistency problem. For units of other shapes, the unit integral on the right side of equation (4), that is, the first term in the right side of equation (12), can be considered to contribute to the object node from the unit volume V by W P The weighted transient term and generated term. Note that W P The unit integral itself is equal to V / 4. The correction term introduced by the Galerkin finite element method, that is, the second term in the right side of equation (12), is the volume of the virtual node domain Vn P After subtracting V / 4 from the equation, it is multiplied by the transient term and generated term at the geometric center. In particular, when the transient phase and generated term are constants in the unit, the [] on the right side of equation (12) becomes a simple

[0234] That is, after correction, the volume of the system considering transient terms and generated terms is equal to the volume of the system considering entry and exit, so that compared with the Galerkin method, equation (12) reaches consistency at a higher second-order accuracy level and satisfies the conservation law. It can be expected that the numerical analysis accuracy of non-regular tetrahedral units can be improved and the mesh dependence can be improved.

[0235] Conservation at the unit level: For units of arbitrary shapes, it can be theoretically proved that the sum of the volumes of the virtual node domains of all nodes of the same unit satisfies equation (13).

[0236] This shows that although the introduction of the correction term changes the allocation of the transient term and the generated term to the four nodes in the unit, the sum of their allocations does not change, which means that the present invention satisfies the conservation law even at the unit level.

[0237] Convergence: Pass Figure 5a and Figure 5b The grid is given an example to illustrate the improvement effect on convergence. Figure 5a Shown are 8 Figure 2 The unit shown (the three faces passing through node 1 are perpendicular to each other, and the vertical intersecting faces have a common side length of h and a volume of V r The grid composed of (represented) is integrated by Galerkin's formula (4) to obtain formula (14).

[0238] For simplicity, it is assumed that the transient term and the generated term are constant in the grid. The subscript A represents the value of C at node A, which is the center point of the grid and becomes the object node P. Figure 2 The subscripts -x, +x, -y, +y, -z, +z represent the C values ​​at the positions where point A moves -h, h in the x direction, -h, h in the y direction, and -h, h in the z direction, respectively. The left side of equation (14) is the second-order accuracy discretization of the second-order differential term in equation (1) multiplied by 4V r , while the transient term and generated term on the right side are the corresponding terms in equation (1) multiplied by 2V r , obviously does not converge to the original differential equation. From equation (10), equation (11) and Table 1, we can get the virtual node domain Vn of node A (i.e. node P) in each unit: P =V r / 2, as a result, the right side of equation (12) becomes equation (15).

[0239] It can be confirmed that equation (16) obtained from equation (12) obviously converges to the original differential equation.

[0240] on the other hand, Figure 5b The node B shown is used together with the node A in the unit partitioning, and it consists of 24 Figure 2The unit shown (the three faces passing through node 1 are perpendicular to each other, and the vertical intersecting faces have a common side length of h and a volume of V r (denoted by ) and surrounded by 8 regular tetrahedral units. Applying the Galerkin method to node B, equation (4) is integrated to obtain equation (17).

[0241] The subscript B represents the value of C at node B, where node B is the center point of the grid and becomes the object node P. The subscripts -x, +x, -y, +y, -z, and +z represent the value of C at a position where point B moves -h, h in the x direction, -h, h in the y direction, and -h, h in the z direction, respectively. The first term on the left side of equation (17) comes from the calculation of 24 volumes V r The second term on the left side of equation (17) comes from the sum of the integrals of the eight regular tetrahedral units and uses the Taylor expansion formula. Here, the volume of the regular tetrahedral unit is equal to 2V r The addition of the two terms on the left is equivalent to multiplying the second-order differential term in equation (1) by 8V. r The transient term and generated term on the right side of equation (17) are the corresponding terms in equation (1) multiplied by 10V r , obviously does not converge to the original differential equation. From equation (10), equation (11) and Table 1, the right side of equation (12) about node B is Among them, 24×V r / 6 is the node B in 24 with volume V r The sum of the virtual node domains of the corresponding nodes in the unit, 8×(2V r ) / 4 is the sum of the virtual node domains of the nodes corresponding to node B in the 8 regular tetrahedral units. The sum of the two is As a result, equation (18) is obtained from equation (12), which converges to the original differential equation.

[0242] These examples demonstrate the effect of introducing correction terms on improving convergence.

[0243] Regarding the discretization of the time derivative term, the specific application of formula (12) in simulation requires the discretization of the time derivative term. Formula (12) can be rewritten as formula (19) after moving the terms.

[0244] Formula (19) holds true at any time, thus formula (20) is obtained.

[0245] Substituting equation (19) into equation (20) yields equation (21).

[0246] Substituting the mean value theorem of equation (22) into equation (21), the partial derivative with respect to time disappears, and the finite element discretization result of the transient diffusion equation after applying the present invention is obtained, and a set of simultaneous equations is obtained.

[0247] All the terms on the right side of the simultaneous equations are known, that is, the concentration value C of each node at time n n is a system of algebraic equations with unknown quantities.

[0248] The precision improvement effect of the present invention is demonstrated by two simulation examples. The simulation domain is defined as a cube region (0≤x, y, z≤1), and k=1. As a mandatory boundary condition problem, the generated term is equation (24), the initial condition is equation (25), and the first-class boundary condition is equation (26). The analytical solution is equation (27).

[0249] As a mixed boundary condition problem, the generated term is equation (28), the initial condition is equation (29), the first type boundary condition part is equation (30), and the second type boundary condition part is equation (31). Its analytical solution is equation (32).

[0250] The numerical error is evaluated using formula (33).

[0251] Figures 6a to 6e Shown are the five tetrahedral mesh modes used. Figure 6a The grid is first divided into cubes with a side length of L = 1 / 8, and then each cube is divided into 5 tetrahedrons. Figure 6b The grid is composed of each cube divided into 6 tetrahedrons. Figure 6c The grid is formed by dividing each cube into two cuboids, and then dividing each cuboid into 14 tetrahedrons. Figure 6d The grid is in Figure 6c Add a node to the geometric center of each tetrahedron and further divide the tetrahedron into 4 flat tetrahedrons. Figure 6e The grid is an irregular grid with density. Considering the corresponding unit side length, Figure 6a , Figure 6b , Figure 6e The grid is set, the time step is Δt=0.005, Figure 6c The grid is set to Δt = 0.002, Figure 6d The grid is taken as Δt=0.001.

[0252] Figure 7 shows the numerical errors in the problem with forced boundary conditions, Figure 8 The numerical error in the mixed boundary condition problem is shown. "GFE" represents the traditional Galerkin finite element method, and "CSFE" represents the new method of the present invention. The horizontal axes (a), (b), ..., (e) represent the numerical errors in the mixed boundary condition problem. Figures 6a to 6e The corresponding grid in . The error of the vertical axis is defined by equation (33). It can be confirmed that the analytical accuracy of CSFE is greatly improved regardless of the type of grid and boundary conditions.

Claims

1. A transient diffusion simulation method based on tetrahedral units, characterized in that: include: Select a simulation area and divide the simulation area into tetrahedral grids; Set the initial conditions, t 0 =0, n=1; Generate a cell matrix: Where C represents concentration, Δt is the time step, t represents time, and C n represents the concentration at the nth moment, C n-1 represents the concentration at the n-1th moment, k represents the diffusion coefficient, ∫ e de represents the integral over element e, where e represents the tetrahedral element around node P, and W P represents the Galerkin weight function, the subscript O represents the value at the geometric center of the unit, Vn P is the volume of the virtual node domain of the object node P in the unit, V represents the volume of the unit, is the time parameter, represents the partial derivative of x, represents the partial derivative of y, It means partial derivative of z, x, y, z represent three-dimensional rectangular coordinates; Generate the overall matrix: Among them, Σ represents the sum of the integrals of each unit; Advance time: t n =t n-1 +△t;t n represents the nth moment, t n-1 represents the n-1th moment; Substituting the concentration at the previous moment into Substitute the generated term S into Generate equivalent nodal loads and generate a system of simultaneous equations: Among them, S n is the generated term at the nth moment, S n-1 is the generated item at the n-1th moment; Import the boundary conditions and solve the simultaneous equations to obtain the concentration C n ; Use time to determine whether it is finished. If it is not finished, then n=n+1 and continue to advance the time until it is finished.

2. The transient diffusion simulation method based on tetrahedron unit according to claim 1 is characterized in that: The process of generating the simultaneous equations is as follows: Applying the Galerkin finite element method to equation (1), we obtain equation (2) for the object node P considered; is a transient term, Ω represents the area around the node P consisting of tetrahedral units, ∫ Ω dΩ represents the integral over the region Ω, W P is the Galerkin weight function, the integral on Ω is the sum of the integrals of the surrounding units, and equation (2) can be rewritten as equation (3); Applying the integration by parts theorem and Gauss-Green theorem to the unit integration on the left side of equation (3), we obtain equation (4) which is the discretization of the Galerkin finite element method; Use equation (12) to replace equation (4) of the Galerkin method; For a unit of any shape, the sum of the volumes of the virtual node domains of all nodes of the same unit satisfies equation (13); Formula (12) can be rewritten as formula (19) after moving the terms; Formula (19) is valid at any time, so we get formula (20); Substituting equation (19) into equation (20), we get equation (21); Substituting the mean value theorem of equation (22) into equation (21), the partial derivative with respect to time vanishes, and the finite element discretization result of the transient diffusion equation is obtained, that is, a system of simultaneous equations is obtained.

3. The transient diffusion simulation method based on tetrahedron unit according to claim 2 is characterized in that: For the unit where three faces intersect perpendicularly at node 1, the coordinate origin is placed at node 1, and nodes 2, 3, and 4 are placed on the positive x, y, and z axes respectively, the unit integral on the left side of equation (4) is calculated as equation (5); The numerical subscript indicates the node number in the unit; and the node number of the object node P in the unit is set to 1; I ij represents the midpoint of the line segment connecting nodes i and j; k(C1-C2) / x2 in the first term on the right side of equation (5) is the midpoint I 12 The second-order accuracy approximation of the diffusion flux from node 1 to node 2 is given by the coefficient y3z4 / 6, which is the flow area of ​​the diffusion flux. This area passes through point I 12 and perpendicular to the diffusion flux; the physical meaning of the first term on the right side of equation (5) is interpreted as the use point I 12 The second-order accuracy diffusion flux at the point is used to evaluate the diffusion flow rate of the flow area y3z4 / 6; the second term on the right side of equation (5) uses the midpoint I 13 The second-order diffusion flux k(C1-C3) / y3 from node 1 to node 3 is used to evaluate the flow area The diffusion rate of x2z4 / 6; the area x2z4 / 6 passes through point I 13 and is perpendicular to the diffusion flux k(C1-C3) / y3; the third term on the right side of equation (5) uses the midpoint I 14 The second-order precision diffusion flux k(C1-C4) / z4 from node 1 to node 4 is used to evaluate the diffusion flow of the flow area x2y3 / 6; area x2y3 / 6 passes through point I 14 And it is perpendicular to the diffusion flux k(C1-C4) / z4; Assume that the object node P is numbered 2, 3, and 4 in the unit, and the unit integral on the left side of equation (4) is calculated as equations (34), (35), and (36) respectively; The flow areas of the second-order diffusion flux from nodes 2, 3, and 4 to node 1 are y3z4 / 6, x2z4 / 6, and x2y3 / 6, respectively. The flow area of ​​the second-order diffusion flux pointing from nodes 2, 3, and 4 to each other is zero.

4. The transient diffusion simulation method based on tetrahedron unit according to claim 2 is characterized in that: For a tetrahedral element of arbitrary shape, the element integral on the right side of equation (4) is calculated, and the contributions from the transient term and the generated term are almost equally distributed to the four nodes of the element in the form of equivalent node loads; when When and S are constants inside the unit, the unit integral is written as (6); The transient and generated terms assigned to node 1 both come from only 1 / 4 of the element volume.

5. The transient diffusion simulation method based on tetrahedron unit according to claim 2 is characterized in that: For tetrahedral elements of arbitrary shapes, the node number of the object node P in the element is set to 1, and the element integration result on the left side of equation (4) is written as equation (8); A 1i is a coefficient that depends only on the coordinates of the element nodes. Its first subscript 1 is due to the fact that the node number of the object node P in the unit is 1. The representation of the second line of formula (8) uses the theoretically proven formula (9). In formula (8), L is introduced 1i (i=2,3,4), let L 1i =[(x i -x1) 2 +(y i -y1) 2 +(z i -z1) 2 ] 1 / 2 , that is, it is defined as the distance between node 1 and node i; k(C1-C i ) / L 1i is the midpoint I of the line segment connecting nodes 1 and i 1i The diffusion flux from node 1 to i at the second-order accuracy, its coefficient (-A 1i L 1i ) is the area; the area passes through the midpoint I 1i and the diffusion flux k(C1-C i ) / L 1i The two are multiplied, and the right side of equation (8) becomes the second-order accuracy evaluation of the flow area (-A 1i L 1i ) diffusion flow; area (-A 1i L 1i ) as the bottom surface, node 1 as the vertex, and L 1i / 2 is the height of the three cones, and their volumes are Vn 1i To express, it is defined by formula (10); VN P =Vn1=Vn 12 +Vn 13 +Vn 14 (11); The volume of the virtual node domain Vn1 of node 1 is the sum of the volumes of the three cones. Since node P and node 1 are identical, Vn P It is also calculated by formula (11).

6. The transient diffusion simulation method based on tetrahedron unit according to claim 2 is characterized in that: For regular tetrahedral units, Vn P =V / 4, correction term is zero; for units of other shapes, the unit integral on the right side of equation (4), that is, the first term in the brackets on the right side of equation (12) The contribution to the object node is considered to come from the unit volume V being P The weighted transient term and generated term; W P The unit integral itself is equal to V / 4. The correction term introduced by the Galerkin finite element method is the volume of the virtual node domain Vn P After V / 4 is subtracted from the equation, it is multiplied by the transient term and generated term at the geometric center. When the transient term and generated term are constants in the unit, the square brackets on the right side of equation (12) become That is, the system volume considering transient terms and generated terms is equal to the system volume when the inlet and outlet are considered with diffusion flow of second-order accuracy.

7. The transient diffusion simulation method based on tetrahedron units according to claim 2 is characterized in that: Based on eight shapes, the three faces passing through node 1 are perpendicular to each other, and the vertical intersecting faces have a total side length of h and a volume of V r The grid composed of the units represented by is integrated using the Galerkin method (4) to obtain the formula (14); Assume that the transient term and the generated term are constants in the grid; the subscript A represents the value of C at node A, where node A is the center point of the grid and becomes the object node P, corresponding to node 1; the subscripts -x, +x, -y, +y, -z, +z represent the value of C at the position where point A moves -h, h in the x direction, -h, h in the y direction, and -h, h in the z direction, respectively; the left side of equation (14) is the second-order accuracy discretization of the second-order differential term in equation (1) multiplied by 4V r , while the transient term and generated term on the right side are the corresponding terms in equation (1) multiplied by 2V r ; From equations (10) and (11), we can get the virtual node domain Vn of node A in each unit: P =V r / 2, as a result, the right side of equation (12) becomes equation (15); Confirm that equation (16) obtained from equation (12) converges to the original differential equation; Based on 24 shapes whose three faces passing through node 1 are perpendicular to each other, the perpendicular intersecting faces have a total side length of h and a volume of V r The unit represented by and the mesh enclosed by 8 regular tetrahedral units are integrated using the Galerkin method (4) to obtain the equation (17); Subscript B represents the value of C at node B, where node B is the center point of the grid and becomes the object node P; subscripts -x, +x, -y, +y, -z, +z represent the value of C at the position where point B moves -h, h in the x direction, -h, h in the y direction, and -h, h in the z direction, respectively; the first term on the left side of equation (17) comes from the calculation of 24 volumes V r The second term on the left side of formula (17) comes from the sum of the integrals of the eight regular tetrahedral units, and uses the Taylor expansion formula; the volume of the regular tetrahedral unit is equal to 2V r The sum of the two terms on the left is equivalent to the second-order differential term in equation (1) multiplied by 8V r The transient term and generated term on the right side of equation (17) are the corresponding terms in equation (1) multiplied by 10V. r ; From equations (10) and (11), the right side of equation (12) about node B is Among them, 24×V r / 6 is the node B in 24 with volume V r The sum of the virtual node domains of the corresponding nodes in the unit, 8×(2V r ) / 4 is the sum of the virtual node domains of the nodes corresponding to node B in the 8 regular tetrahedral units; the sum of the two is As a result, equation (18) converges to the original differential equation from equation (12); 8. The transient diffusion simulation method based on tetrahedron units according to claim 2 is characterized in that: When the simulation domain is a cubic region, 0≤x≤1, 0≤y≤1, 0≤z≤1, set k=1; for the numerical error in the forced boundary condition problem, the generated term is taken as equation (24), the initial condition is taken as equation (25), and the first type of boundary condition is taken as equation (26), and the analytical solution is obtained as equation (27); C 0 =sin(πx)sin(πy)sin(πz) (t=0) (25); C=0(x=0,1, ory=0,1, orz=0,1)(26); The numerical error is evaluated using formula (33); Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the simulation solution of node i, and the maximum value of the analytical solution in the simulation area, respectively; the superscript end represents the end time of the calculation, the numerical simulation is performed with double precision, and the unit integration is performed through four-point numerical integration; End time t end =1.

9. The transient diffusion simulation method based on tetrahedron units according to claim 2, characterized in that: When the simulation domain is a cubic region, 0≤x≤1, 0≤y≤1, 0≤z≤1, set k=1; for the numerical error in the mixed boundary condition problem, the generated term is taken as equation (28), the initial condition is taken as equation (29), the first type of boundary condition is taken as equation (30), and the second type of boundary condition is taken as equation (31), and the analytical solution is obtained as equation (32); In formula (31), n ​​is the external normal of unit length; The numerical error is evaluated using formula (33); Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the simulation solution of node i, and the maximum value of the analytical solution in the simulation area, respectively; the superscript end represents the end time of the calculation, the numerical simulation is performed with double precision, and the unit integration is performed through four-point numerical integration; End time t end =1.

10. The transient diffusion simulation system based on tetrahedron unit is characterized by: Used to implement the transient diffusion simulation method based on tetrahedron units according to any one of claims 1 to 9, the transient diffusion simulation system based on tetrahedron units comprises: A partitioning module is used to select a simulation area and divide the simulation area into tetrahedral grids; Setting module, used to set initial conditions, t 0 =0, n=1; Unit matrix generation module, used to generate unit matrix: Where C represents concentration, Δt is the time step, t represents time, and C n represents the concentration at the nth moment, C n-1 represents the concentration at the n-1th moment, k represents the diffusion coefficient, ∫ e de represents the integral over element e, where e represents the tetrahedral element around node P, and W P represents the Galerkin weight function, the subscript O represents the value at the geometric center of the unit, Vn P is the volume of the virtual node domain of the object node P in the unit, V represents the volume of the unit, is the time parameter, represents the partial derivative of x, represents the partial derivative of y, It means partial derivative of z, x, y, z represent three-dimensional rectangular coordinates; The overall matrix generation module is used to generate the overall matrix: Among them, Σ represents the sum of the integrals of each unit; Time advancement module, used to advance time: t n =t n-1 +△t;t n represents the nth moment, t n-1 represents the n-1th moment; The simultaneous equations generation module is used to substitute the concentration at the previous moment into Substitute the generated term S into Generate equivalent nodal loads and generate a system of simultaneous equations: Among them, S n is the generated term at the nth moment, S n-1 is the generated item at the n-1th moment; The simultaneous equations solving module is used to import boundary conditions, solve the simultaneous equations, and obtain the concentration C n ; The end module is used to determine whether it is ended based on time. If it is not ended, n=n+1 and the time continues to advance until it is ended.

Images