Quadrilateral unit transient heat conduction simulation method and system based on virtual node domain

Abstract

The invention provides a quadrilateral unit transient heat conduction simulation method and system based on a virtual node domain, and belongs to the technical field of transient heat conduction simulation, and the method comprises the steps: selecting a simulation region, and dividing the simulation region into quadrilateral grids; setting initial conditions; generating a unit matrix; generating an overall matrix; advancing time; the temperature at the previous moment is substituted into the overall matrix, the source item is substituted into the overall matrix, an equivalent node load is generated, and a simultaneous equation set is generated; boundary conditions are imported, the simultaneous equation set is solved, and the temperature Tn is obtained; 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 a quadrilateral unit, the problem that when the Galerkin finite element method is investigated by the heat flux density with second-order precision, inconsistency exists between a system considering entering and exiting and a system considering a transient term and a source term is solved.

Description

Technical Field

[0001] The present invention relates to the technical field of transient heat conduction simulation, and in particular to a method and system for simulating transient heat conduction of quadrilateral units based on virtual node domains. Background Art

[0002] The Galerkin finite element method is one of the most commonly used numerical solutions to differential equations, with its clear error assessment and wide range of applications. However, its mesh pattern dependence and the increase in numerical errors due to poor unit quality have not yet been resolved. For example, when discretizing the Poisson equation on a quadrilateral mesh, the discretization does not converge to the original differential equation due to the unit shape, and the maximum principle cannot be satisfied. Similar problems exist in the finite element method simulation of transient heat conduction in addition to the Poisson equation.

[0003] The Galerkin finite element method is based on a second-order accurate system when considering the inflow and outflow of heat. This is different from the system used when considering transient terms and source terms, and cannot satisfy the conservation law with second-order accuracy. This system difference problem is called inconsistency. For example Figure 4 Node 1. Similarly, there are inconsistency problems for other nodes, and other shapes of units with inconsistency problems can also be listed. Summary of the invention

[0004] The present invention provides a method and system for simulating transient heat conduction of quadrilateral units based on virtual node domains, and solves the inconsistency problem between the system when considering the in and out and the system when considering the transient term and the source term when the Galerkin finite element method is investigated with the heat flux density of the second order accuracy for the quadrilateral units.

[0005] This specification discloses a quadrilateral unit transient heat conduction simulation method based on a virtual node domain, comprising:

[0006] Select a simulation area and divide the simulation area into quadrilateral grids;

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

[0008] Generate a cell matrix:

[0009]

[0010]

[0011] Where Δt is the time step, t is the time, T n represents the temperature at the nth moment, T n-1 represents the temperature at the n-1th moment, λ represents the thermal conductivity, W Prepresents the Galerkin weight function, ξ and η are the horizontal and vertical coordinates of the isoparametric unit respectively; represents the integral over isoparametric units; is the Jacobian determinant; the subscript O indicates the value at the geometric center of the unit, Vn P is the area of ​​the virtual node domain of the target node P in the unit, is the time parameter, represents the partial derivative of x, It means to find the partial derivative of y, x and y represent the plane 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] Substitute the temperature at the previous moment into Substituting the source term Q into

[0017] Generate equivalent nodal loads and generate a system of simultaneous equations:

[0018]

[0019] Among them, Q n is the source term at the nth moment, Q n-1 is the source term at the n-1th moment;

[0020] Import the boundary conditions and solve the simultaneous equations to obtain the temperature T 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 target node P considered;

[0024]

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

[0026]

[0027] e represents the unit around the node P, ∫ e de represents the integral on the element e, Σ represents the sum of the integrals of each element; for the element integral on the left side of equation (3), the integration by parts theorem and Gauss-Green theorem are applied to obtain equation (4) which is the discretization of the Galerkin finite element method;

[0028]

[0029] For quadrilateral elements, the integral of each element in equation (4) becomes the unified equation (5);

[0030]

[0031] Use equation (14) to replace equation (5) of the Galerkin method;

[0032]

[0033] Vn P is the virtual node domain area of ​​the target node P in the unit;

[0034] Theoretically, it can be proved that the sum of the areas of the virtual node domains of all nodes of the same unit satisfies equation (15), and the conservation law is satisfied at the unit level;

[0035]

[0036] Where S represents the area of ​​the unit;

[0037] After moving the terms in equation (14), it can be rewritten as equation (18);

[0038]

[0039] Formula (18) is valid at any time, that is, formula (19) is obtained;

[0040]

[0041] Substituting formula (18) into formula (19) yields formula (20);

[0042]

[0043] Substituting the mean value theorem of equation (21) into equation (20), the time derivative term disappears, and the finite element discretization result of the transient heat conduction equation is obtained, that is, a system of simultaneous equations is obtained.

[0044] In this manual, Figure 4 The parallelogram unit (a parallelogram unit with a height equal to the base and tilted π / 4), the unit integral on the left side of equation (5) is calculated as equation (6);

[0045]

[0046] The node number of the target node P in the unit is 1, and the node number in the unit is represented by a digital subscript; T1, T2, T3 represent the temperature of nodes 1, 2, and 3 respectively, h is the height, I ij represents the midpoint of the line segment connecting nodes i and j; λ(T1-T2) / h in the first term on the right side of equation (6) is the midpoint I 12 The heat flux density from node 1 to node 2 is approximated by the second order accuracy, and its coefficient h / 2 is the heat flux density from node 1 to node 2. 12 To point I 23 The first term on the right side of equation (6) is the distance of point I. 12 The heat flux at the vertical line segment I is evaluated with second-order accuracy. 12 I 23 The second term on the right side of equation (6) uses the second-order accuracy heat flux density from node 1 to node 3 at the geometric center of the unit O to evaluate the heat flux through the vertical line segment I 23 I 14 The heat flow from node 1 to 3 and the line segment I 23 I 14 The right side of equation (6) indicates that the control line I 12 I 23 and I 23 I 14 Outgoing heat flux; the area enclosed by the cell boundary and the control line is a grayed-out quadrilateral 1I 12 I 23 I 14 , whose area is 3 / 8 of the unit area.

[0047] In this manual, Figure 4 The parallelogram element on the right side of equation (5) is calculated by element integration. The contributions from the transient term and the source term are almost equally distributed to the four nodes of the element in the form of equivalent node loads. When and Q are constants inside the unit, the unit integral is as shown in equation (7);

[0048]

[0049] The transient term and source term assigned to node 1 come from only 1 / 4 of the unit area, which is the same as that of quadrilateral 1I. 12 I 23 I14 The areas are quite different, and there is an inconsistency problem from the perspective of conservation law.

[0050] In this specification, for a quadrilateral unit of any shape, the node number of the target node P in the unit is set to 1, and the integral result of the left side of equation (5) is equation (9);

[0051]

[0052] A 1i It is a coefficient that depends only on the coordinates of the element nodes. Its first subscript 1 is due to the node number of the target node P in the unit. Theoretically, it can be proved that formula (10) is valid.

[0053]

[0054] Using formula (10), formula (9) can be rewritten as formula (11);

[0055]

[0056] Introducing L 1i (i=2,3,4), let L 1i =[(x i -x1) 2 +(y i -y1) 2 ] 1 / 2 , that is, it is defined as the distance between node 1 and node i; λ(T1-T i ) / L 1i is the midpoint I of the line segment connecting nodes 1 and i 1i The heat flux density at the second-order accuracy from node 1 to i at the position is calculated based on its coefficient (-A 1i L 1i ) is the length (-A 1i L 1i ) that passes through the midpoint I 1i And it is perpendicular to the heat flux, that is, perpendicular to the line segment connecting nodes 1 and i; multiplying the two, the right side of formula (11) becomes the passing line segment (-A) evaluated with second-order accuracy 1i L 1i ) heat flow; in the unit, the line segment (-A 12 L 12 )、(-A 13 L 13 )、(-A 14 L 14 ) as the base, and the three triangles with node 1 as the vertex, their areas are respectively expressed as Vn 12 、Vn 13 、Vn 14To express, the areas of these three triangles are calculated by formula (12);

[0057]

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

[0059] The area of ​​the virtual node domain Vn1 of node 1 is the sum of the areas of these three triangles. Since node P and node 1 are identical, Vn P It is also calculated by formula (13).

[0060] In this specification, for square cells and rectangular cells, Vn P =S / 4, the correction term is zero; for example Figure 4 The unit shown in the figure has Vn1=3S / 8 obtained by equation (13), which is consistent with the quadrilateral 1I 12 I 23 I 14 coincidence; the right side of formula (5), that is, the first term in the brackets on the right side of formula (14) The contribution of the cell integral to the target node is considered to come from the cell area S divided by W P The weighted transient term and source term; the correction term introduced to the Galerkin finite element method, that is, the second term in the bracket on the right side of equation (14) (Vn P - Its area is from the virtual node domain Vn P Subtracted Then multiply it by the transient term and source term at the geometric center O; when the transient phase and source terms are constants in the unit, the square brackets on the right side of equation (14) become That is, the system area considering transient terms and source terms is equal to the system area considering entry and exit.

[0061] In this manual, there are eight Figure 4 The grid composed of parallelograms is integrated using equation (5) to obtain equation (16);

[0062]

[0063] Here, the transient term and source term are assumed to be constant in the grid, and the subscripts -x, +x, -y, and +y represent Figure 6The corresponding points in (1); the left side of equation (16) is the discretization of the second-order accuracy of the second-order differential term in equation (1), and the transient term and source term on the right side are twice the corresponding terms in equation (1), which obviously does not converge to the original differential equation; the virtual node domain Vn2 = S / 8 of node 2 is obtained from equation (13), and equation (14) after introducing the correction term is used to obtain equation (17) that converges to the original differential equation.

[0064] In this specification, when the simulation domain is a square area, 0≤x≤1, 0≤y≤1, λ=1; for the numerical error in the forced boundary condition problem, the source term is taken as equation (23), the initial condition is taken as equation (24), and the forced boundary condition is taken as equation (25), and the analytical solution is obtained as equation (26);

[0065]

[0066] T 0 = sin(πx)sin(πy) (24);

[0067] T=0(x=0,1, or y=0,1) (25);

[0068]

[0069] The numerical value error is evaluated using formula (32);

[0070]

[0071] Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the numerical 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; the time step is Δt = 0.005, End time t end =1.

[0072] In this specification, when the simulation domain is a square area, 0≤x≤1, 0≤y≤1, λ=1; for the numerical error in the mixed boundary condition problem, equation (27) is used as the source term, equation (28) is used as the initial condition, equation (29) in the mixed boundary condition is used as the forced boundary condition, and equation (30) is used as the natural boundary condition to obtain the analytical solution as equation (31);

[0073]

[0074] In formula (30), n represents the unit external normal direction;

[0075] The numerical value error is evaluated using formula (32);

[0076]

[0077] Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the numerical 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; the time step is Δt = 0.005, End time t end =1.

[0078] This specification also discloses a quadrilateral unit transient heat conduction simulation system based on a virtual node domain, which is used to implement any of the quadrilateral unit transient heat conduction simulation methods based on a virtual node domain described above. The quadrilateral unit transient heat conduction simulation system based on a virtual node domain includes:

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

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

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

[0082]

[0083] Where Δt is the time step, t is the time, T n represents the temperature at the nth moment, T n-1 represents the temperature at the n-1th moment, λ represents the thermal conductivity, W P represents the Galerkin weight function, ξ and η are the horizontal and vertical coordinates of the isoparametric unit respectively; represents the integral over isoparametric units; is the Jacobian determinant; the subscript O indicates the value at the geometric center of the unit, Vn P is the area of ​​the virtual node domain of the target node P in the unit, is the time parameter, represents the partial derivative of x, It means to find the partial derivative of y, x and y represent the plane coordinates;

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

[0085]

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

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

[0088] The module for generating simultaneous equations is used to substitute the temperature at the previous moment into Substituting the source term Q into Generate equivalent nodal loads and generate a system of simultaneous equations:

[0089]

[0090] Among them, Q n is the source term at the nth moment, Q n-1 is the source term at the n-1th moment;

[0091] The simultaneous equations solving module is used to import boundary conditions, solve the simultaneous equations, and obtain the temperature T n ;

[0092] 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.

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

[0094] The present invention proposes a new processing method for transient terms and source terms when using quadrilateral units for finite element simulation, so that the system used when considering them is equivalent to or close to the system when considering entry and exit, and the conservation law is satisfied at the node level with higher accuracy than the traditional Galerkin method, thereby improving the dependency problem of the grid mode and improving the accuracy of the numerical solution obtained by low-quality units. BRIEF DESCRIPTION OF THE DRAWINGS

[0095] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings required for use in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other accompanying drawings can be obtained based on these accompanying drawings without paying creative work.

[0096] Figure 1 It is a schematic diagram of the transient heat conduction simulation method of quadrilateral unit based on virtual node domain involved in the present invention.

[0097] Figure 2 It is a schematic diagram of the integration area of ​​the node P involved in the present invention.

[0098] Figure 3a It is a schematic diagram of the unit coordinates involved in the present invention.

[0099] Figure 3bIt is a schematic diagram of the isoparametric unit coordinates involved in the present invention.

[0100] Figure 4 Schematic diagram of the heat flux investigation of a parallelogram unit with a height equal to the base and an inclination of π / 4 involved in the present invention.

[0101] Figure 5 It is a schematic diagram of a virtual node domain of any unit involved in the present invention.

[0102] Figure 6 The present invention involves Figure 4 Schematic diagram of a parallelogram unit grid.

[0103] Figure 7a Schematic diagram of a square unit grid (81 nodes) involved in the present invention.

[0104] Figure 7b It is a schematic diagram of 3 to 8 (specific numbers are 3, 4, 5, 8) unit grids (69 nodes) involved in the present invention.

[0105] Figure 7c It is a schematic diagram of the diamond unit grid (108 nodes) involved in the present invention.

[0106] Figure 7d Schematic diagram of the rectangular quadrilateral grid (141 nodes) involved in the present invention.

[0107] Figure 7e Schematic diagram of the parallelogram grid (279 nodes) involved in the present invention.

[0108] Figure 7f It is a schematic diagram of a terrace-type unit grid (353 nodes) involved in the present invention.

[0109] Figure 7g It is a schematic diagram of the irregular quadrilateral grid (84 nodes) involved in the present invention.

[0110] Figure 8 It is a schematic diagram of the simulation error under the forced boundary conditions involved in the present invention.

[0111] Fig. 9 It is a schematic diagram of the simulation error under the mixed boundary conditions involved in the present invention. DETAILED DESCRIPTION

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

[0113] like Figure 1 As shown, the embodiment of this specification discloses a quadrilateral unit transient heat conduction simulation method based on a virtual node domain, including:

[0114] Select a simulation area and divide the simulation area into quadrilateral grids;

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

[0116] Generate a cell matrix:

[0117]

[0118] Where Δt is the time step, t is the time, T n represents the temperature at the nth moment, T n-1 represents the temperature at the n-1th moment, λ represents the thermal conductivity, W P represents the Galerkin weight function, ξ and η are the horizontal and vertical coordinates of the isoparametric unit respectively; represents the integral over isoparametric units; is the Jacobian determinant; the subscript O indicates the value at the geometric center of the unit, Vn P is the area of ​​the virtual node domain of the target node P in the unit, is the time parameter, represents the partial derivative of x, It means to find the partial derivative of y, x and y represent the plane coordinates;

[0119] Generate the overall matrix:

[0120]

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

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

[0123] Substitute the temperature at the previous moment into Substituting the source term Q into Generate equivalent nodal loads and generate a system of simultaneous equations:

[0124]

[0125] Among them, Q n is the source term at the nth moment, Q n-1 is the source term at the n-1th moment;

[0126] Import the boundary conditions and solve the simultaneous equations to obtain the temperature T n ;

[0127] 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.

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

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

[0130]

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

[0132]

[0133] e represents the unit around the node P, ∫ e de represents the integral on the element e, Σ represents the sum of the integrals of each element; for the element integral on the left side of equation (3), the integration by parts theorem and Gauss-Green theorem are applied to obtain equation (4) which is the discretization of the Galerkin finite element method;

[0134]

[0135] For quadrilateral elements, the integral of each element in equation (4) becomes the unified equation (5);

[0136]

[0137] Use equation (14) to replace equation (5) of the Galerkin method;

[0138]

[0139] Vn P is the virtual node domain area of ​​the target node P in the unit;

[0140] Theoretically, it can be proved that the sum of the areas of the virtual node domains of all nodes of the same unit satisfies equation (15), and the conservation law is satisfied at the unit level;

[0141]

[0142] After moving the terms in equation (14), it can be rewritten as equation (18);

[0143]

[0144] Formula (18) is valid at any time, that is, formula (19) is obtained;

[0145]

[0146] Substituting formula (18) into formula (19) yields formula (20);

[0147]

[0148] Substituting the mean value theorem of equation (21) into equation (20), the time derivative term disappears, and the finite element discretization result of the transient heat conduction equation is obtained, that is, a system of simultaneous equations is obtained.

[0149] In some embodiments, for Figure 4 The parallelogram unit, the unit integral on the left side of equation (5) is calculated as equation (6);

[0150]

[0151] The node number of the target node P in the unit is 1, and the node number in the unit is represented by a digital subscript; T1, T2, T3 represent the temperature of nodes 1, 2, and 3 respectively, h is the height, I ij represents the midpoint of the line segment connecting nodes i and j; λ(T1-T2) / h in the first term on the right side of equation (6) is the midpoint I 12 The heat flux density from node 1 to node 2 is approximated by the second order accuracy, and its coefficient h / 2 is the heat flux density from node 1 to node 2. 12 To point I 23 The first term on the right side of equation (6) is the distance of point I. 12 The heat flux at the vertical line segment I is evaluated with second-order accuracy. 12 I 23 The second term on the right side of equation (6) uses the second-order accuracy heat flux density from node 1 to node 3 at the geometric center of the unit O to evaluate the heat flux through the vertical line segment I 23 I 14 The heat flow from node 1 to 3 and the line segment I 23 I 14 The right side of equation (6) indicates that the control line I 12 I 23 and I 23 I 14 Outgoing heat flux; the area enclosed by the cell boundary and the control line is a grayed-out quadrilateral 1I 12 I 23 I 14 , whose area is 3 / 8 of the unit area.

[0152] Figure 4The parallelogram of is a specific type of parallelogram. It is cited among many quadrilaterals to simply illustrate the inconsistency of Galerkin's method. The result of formula (6) is obtained under this specific parallelogram.

[0153] In some embodiments, for Figure 4 The parallelogram element on the right side of equation (5) is calculated by element integration. The contributions from the transient term and the source term are almost equally distributed to the four nodes of the element in the form of equivalent node loads. When and Q are constants inside the unit, the unit integral is as shown in equation (7);

[0154]

[0155] The transient term and source term assigned to node 1 come from only 1 / 4 of the unit area, which is the same as that of quadrilateral 1I. 12 I 23 I 14 The areas are quite different, and there is an inconsistency problem from the perspective of conservation law.

[0156] In some embodiments, for a quadrilateral unit of any shape, the node number of the target node P in the unit is set to 1, and the integral result of the left side of equation (5) is equation (9);

[0157]

[0158] A 1i It is a coefficient that depends only on the coordinates of the element nodes. Its first subscript 1 is due to the node number of the target node P in the unit. Theoretically, it can be proved that formula (10) is valid.

[0159]

[0160] Using formula (10), formula (9) can be rewritten as formula (11);

[0161]

[0162] Introducing L 1i (i=2,3,4), let L 1i =[(x i -x1) 2 +(y i -y1) 2 ] 1 / 2 , that is, it is defined as the distance between node 1 and node i; λ(T1-T i ) / L 1i is the midpoint I of the line segment connecting nodes 1 and i 1i The heat flux density at the second-order accuracy from node 1 to i at the position is calculated based on its coefficient (-A1i L 1i ) is the length (-A 1i L 1i ) that passes through the midpoint I 1i And it is perpendicular to the heat flux, that is, perpendicular to the line segment connecting nodes 1 and i; multiplying the two, the right side of formula (11) becomes the passing line segment (-A) evaluated with second-order accuracy 1i L 1i ) heat flow; in the unit, the line segment (-A 12 L 12 )、(-A 13 L 13 )、(-A 14 L 14 ) as the base, and the three triangles with node 1 as the vertex, their areas are respectively expressed as Vn 12 、Vn 13 、Vn 14 To express, the areas of these three triangles are calculated by formula (12);

[0163]

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

[0165] The area of ​​the virtual node domain Vn1 of node 1 is the sum of the areas of these three triangles. Since node P and node 1 are identical, Vn P It is also calculated by formula (13).

[0166] In some embodiments, for square cells and rectangular cells, Vn P =S / 4, the correction term is zero; for example Figure 4 The unit shown in the figure has Vn1=3S / 8 obtained by equation (13), which is consistent with the quadrilateral 1I 12 I 23 I 14 coincidence; the right side of formula (5), that is, the first term in the brackets on the right side of formula (14) The contribution of the cell integral to the target node is considered to come from the cell area S divided by W P The weighted transient term and source term; the correction term introduced to the Galerkin finite element method, that is, the second term in the bracket on the right side of equation (14) Its area is from the virtual node domain Vn P Subtracted Then multiply it by the transient term and source term at the geometric center O; when the transient phase and source terms are constants in the unit, the square brackets on the right side of equation (14) become That is, the system area considering transient terms and source terms is equal to the system area considering entry and exit.

[0167] In some embodiments, eight Figure 4 The grid composed of parallelograms is integrated using equation (5) to obtain equation (16);

[0168]

[0169] Here, the transient term and source term are assumed to be constant in the grid, and the subscripts -x, +x, -y, and +y represent Figure 6 The corresponding points in (1); the left side of equation (16) is the discretization of the second-order accuracy of the second-order differential term in equation (1), and the transient term and source term on the right side are twice the corresponding terms in equation (1), which obviously does not converge to the original differential equation; the virtual node domain Vn2 = S / 8 of node 2 is obtained from equation (13), and equation (14) after introducing the correction term is used to obtain equation (17) that converges to the original differential equation.

[0170] In some embodiments, when the simulation domain is a square area, 0≤x≤1, 0≤y≤1, set λ=1; for the numerical error in the forced boundary condition problem, the source term is taken as equation (23), the initial condition is taken as equation (24), and the forced boundary condition is taken as equation (25), and the analytical solution is obtained as equation (26);

[0171]

[0172] T 0 = sin(πx)sin(πy) (24);

[0173] T=0(x=0,1, or y=0,1) (25);

[0174]

[0175] The numerical value error is evaluated using formula (32);

[0176]

[0177] Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the numerical 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; the time step is Δt = 0.005, End time t end =1.

[0178] In some embodiments, when the simulation domain is a square area, 0≤x≤1, 0≤y≤1, set λ=1; for the numerical error in the mixed boundary condition problem, take equation (27) as the source term, take equation (28) as the initial condition, take equation (29) in the mixed boundary condition as the forced boundary condition, take equation (30) as the natural boundary condition to obtain the analytical solution as equation (31);

[0179]

[0180] In formula (30), n represents the unit external normal direction;

[0181] The numerical value error is evaluated using formula (32);

[0182]

[0183] Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the numerical 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; the time step is Δt = 0.005, End time t end =1.

[0184] The embodiment of this specification also discloses a quadrilateral unit transient heat conduction simulation system based on a virtual node domain, which is used to implement any of the quadrilateral unit transient heat conduction simulation methods based on a virtual node domain described above. The quadrilateral unit transient heat conduction simulation system based on a virtual node domain includes:

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

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

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

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

[0189] 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;

[0190] The module for generating simultaneous equations is used to substitute the temperature at the previous moment into Substituting the source term Q into Generate equivalent nodal loads and generate simultaneous equations;

[0191] The simultaneous equations solving module is used to import boundary conditions, solve the simultaneous equations, and obtain the temperature T n ;

[0192] 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.

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

[0194] The governing equation for the transient heat conduction problem is shown in equation (1).

[0195] Here, t represents time, is a transient term, x and y are plane coordinates, They represent the partial derivatives of x and y respectively, T represents temperature, λ represents the thermal conductivity, and Q is the source term.

[0196] The Galerkin finite element method is applied to equation (1), and the integral equation of equation (2) is obtained for the target node P considered.

[0197] Here, Ω represents Figure 2 The area around the node P shown is composed of several 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, so equation (2) can be rewritten as equation (3).

[0198] Here, e represents the unit around the node P, ∫ e de represents the integral on the element, and Σ represents the sum of the integrals of each element. Applying the integration by parts theorem and the Gauss-Green theorem to the element integral on the left side of equation (3), equation (4) is obtained as the discretization of the Galerkin finite element method.

[0199] Since this manual does not make any improvements to the natural boundary conditions, they will not be described here.

[0200] For quadrilateral elements, the finite element method introduces a square with a side length of 2, called an isoparametric element. The integral on the x and y planes is converted into the integral on the ξ and η planes of the isoparametric element through coordinate transformation. ξ and η are the horizontal and vertical coordinates of the isoparametric element, respectively. Figure 3a and Figure 3b The transformation of arbitrary shaped cells into isoparametric cells is shown. As a result of the coordinate transformation, the individual cell integrals in equation (4) are transformed into equation (5) in a unified form.

[0201] In the formula, represents the integral over isoparametric units, It is called the Jacobian determinant.

[0202] for Figure 4 For the listed parallelogram units, the unit integration results on the left side of equation (5) are analyzed and sorted to obtain equation (6).

[0203] Here, the node number in the unit is used as usual, represented by a digital subscript, and the node number of the target node P in the unit is set to 1. ij represents the midpoint of the line segment connecting nodes i and j. Therefore, λ(T1-T2) / h in the first term on the right side of equation (6) is the midpoint I 12 The heat flux density from node 1 to node 2 is approximated by the second order accuracy, and its coefficient h / 2 is the heat flux density from node 1 to node 2. 12 To point I 23 The physical meaning of the first term on the right side of equation (6) can be interpreted as the distance from point I 12 The heat flux at the vertical line segment I is evaluated with second-order accuracy. 12 I 23 Similarly, the second term on the right side of equation (6) uses the second-order heat flux density from node 1 to node 3 at the geometric center of the unit O to evaluate the heat flux through the vertical line segment I 23 I 14 Here, the heat flux density from node 1 to 3 is the same as the line segment I 23 I 14 To summarize the above discussion, the right side of equation (6) represents the control line I 12 I 23 and I 23 I 14 The area enclosed by the cell boundary and the control line is a grayed-out quadrilateral 1I 12 I 23 I 14 . , its area is 3 / 8 of the unit area.

[0204] On the other hand, for Figure 4 The parallelogram element is calculated by the element integral on the right side of equation (5), and the contributions from the transient term and the source term are almost equally distributed to the four nodes of the element in the form of equivalent node loads. In particular, when When and Q are constants inside the unit, the unit integral is as shown in equation (7).

[0205] Here, S represents the area of ​​the cell. The transient term and source term assigned to node 1 both come from only 1 / 4 of the cell area.

[0206] When observing the same system from time t to time t+△t, the conservation law of physical quantities becomes equation (8).

[0207]

[0208] In summary, for Figure 4 For the node 1 of the quadrilateral element, the system of the Galerkin finite element method when considering the heat in and out is quadrilateral 1I 12 I 23 I 14 , but when considering transient terms and source terms, another system with nearly 1 / 4 of the unit area is used, 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, and other shapes of units with inconsistency problems can also be listed.

[0209] However, there is no effective improvement measure for the inconsistency problem of Galerkin finite element method. In simulation practice, in order to reduce the error caused by inconsistency, human and financial resources are invested in meshing, and the unit shape is as close to square or rectangle as possible.

[0210] For quadrilateral elements, the present invention aims to solve the inconsistency problem between the system when considering entry and exit and the system when considering transient terms and source terms when examining the Galerkin finite element method with second-order accuracy of heat flux. The present invention proposes a new method for processing transient terms and source terms when using quadrilateral elements for finite element simulation, so that the system used when considering them is equivalent to or close to the system when considering entry and exit, and satisfies the conservation law at the node level with higher accuracy than the traditional Galerkin method, thereby improving the dependency problem of the grid mode and improving the accuracy of the numerical solution obtained by low-quality elements.

[0211] The means for solving the problem in this specification (i.e., the main component of the invention) is to propose a new concept of virtual node domain; for a quadrilateral unit of any shape, the node number of the target node P in the unit is set to 1, and the integral result of the left side of equation (5) is equation (9).

[0212] Here A 1i It is a coefficient that depends only on the coordinates of the element nodes, and its first subscript "1" is due to the node number of the target node P in the unit. Theoretically, it can be proved that formula (10) is valid.

[0213] Using equation (10), equation (9) can be rewritten as equation (11).

[0214] Here we introduce L 1i (i=2,3,4), let L 1i =[(x i -x1) 2 +(y i -y1) 2 ]1 / 2 , that is, it is defined as the distance between node 1 and node i. The λ(T1-T i ) / L 1i is the midpoint I of the line segment connecting nodes 1 and i 1i The heat flux density at the second-order accuracy from node 1 to i at . 1i L 1i ) is the length (-A 1i L 1i ) that passes through the midpoint I 1i And it is perpendicular to the heat flux, that is, perpendicular to the line segment connecting nodes 1 and i. Multiplying the two, the right side of equation (11) becomes the passing line segment (-A) evaluated with second-order accuracy. 1i L 1i ) heat flow. Figure 5 In the figure, the line segments (-A 12 L 12 )、(-A 13 L 13 )、(-A 14 L 14 ) as the base, and the three triangles with node 1 as the vertex, their areas are respectively expressed as Vn 12 、Vn 13 、Vn 14 The areas of the three triangles are calculated by equation (12). Due to the shape of the unit, the areas of the triangles may be negative, but this does not require any special computational processing and does not affect the numerical simulation to be implemented.

[0215] The area of ​​the virtual node domain Vn1 of node 1 is the sum of the areas of these three triangles. Since node P and node 1 are identical, Vn P It is also calculated by formula (13).

[0216] For a quadrilateral of any shape, the heat flux density passing lines cannot be connected to form a broken line (in the case of a rectangle, they can be connected to form a broken line). Therefore, for the virtual node domain, its area is defined by equations (12) and (13), but the shape cannot be defined. In this sense, Figure 5 The areas of the virtual node domains are shown rather than their shapes.

[0217] Here, the treatment of transient and source terms is modified to match the in-out evaluation system with second-order accuracy. We propose to replace Galerkin method (5) with (14).

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

[0219] Regarding the correction term in equation (14), for square cells and rectangular cells, Vn P =S / 4, the correction term is zero. This is because there is no inconsistency between the two types of units. Figure 4 For the unit shown in the figure, Vn1 = 3S / 8 can be obtained from formula (13), which is consistent with the quadrilateral 1I 12 I 23 I 14 The inconsistency is also eliminated for other nodes. For any shape of element, the contribution of the element integral to the target node (called equivalent node load) in the right side of equation (5), that is, the first term in the right side of equation (14) can be considered to come from the unit area S divided by W. P The weighted transient term and source term. The correction term we introduced to the Galerkin finite element method, that is, the second term in the right side of equation (14), is the area of ​​the virtual node domain Vn P Subtracted Then multiply it by the transient term and source term at the geometric center O. In particular, when the transient phase and source terms are constants in the unit, the [] on the right side of equation (14) becomes That is, the system area considering transient terms and source terms is equal to the system area considering entry and exit, so that formula (14) is consistent with the traditional Galerkin method at a higher second-order accuracy level and satisfies the conservation law, which can be expected to improve the simulation accuracy.

[0220] Regarding conservation at the unit level, for units of arbitrary shapes, it can be theoretically proved that the sum of the areas of the virtual node domains of all nodes of the same unit satisfies equation (15).

[0221] This shows that although the introduction of the correction terms (except squares and rectangles) changes the distribution amounts of the transient terms and source terms of the four nodes in the unit, the sum of their distribution amounts does not change, which means that the present invention satisfies the conservation law even at the unit level.

[0222] Regarding convergence, Figure 6 The grid is given an example to illustrate the improvement effect on convergence. Figure 6 Shown are eight Figure 4 The grid composed of parallelograms is integrated using the Galerkin method (5) to obtain (16).

[0223] For simplicity, it is assumed that the transient term and source term are constants in the grid. The left side of equation (16) is the discretization of the second-order differential term in equation (1) with second-order accuracy, while the transient term and source term on the right side are twice the corresponding terms in equation (1), which obviously does not converge to the original differential equation. From equation (13), we get Figure 4 The virtual node domain Vn2 of node 2 is S / 8. By using equation (14) after introducing the correction term, we can obtain equation (17) that converges to the original differential equation.

[0224] Other examples can be listed where the introduction of correction terms turns a non-convergent mesh into a convergent mesh, such as Figure 7f This shows that the introduction of the correction term improves the convergence.

[0225] Regarding the discretization of the time derivative term, equation (14) can be rewritten as equation (18) after shifting the terms.

[0226] Formula (18) holds true at any time, thus we get Formula (19).

[0227] Here, is the time parameter, the superscript n indicates the nth moment, and n-1 indicates the physical quantity at the previous moment. Substituting equation (18) into equation (19) yields equation (20).

[0228] Assuming the time step is Δt, substituting the mean value theorem of equation (21) into equation (20), the time derivative term disappears, and the finite element discretization result of the transient heat conduction equation after applying the present invention is obtained, that is, a set of simultaneous equations.

[0229] All the terms on the right side of the simultaneous equations are known, that is, the temperature value T of each node at time n n is a system of algebraic equations with unknown variables.

[0230] The accuracy improvement effect of the present invention is demonstrated by simulation examples of heat conduction problems under forced boundary conditions and mixed boundary conditions. Referring to the embodiments of the above equations (23) to (32).

[0231] Figure 7a to Figure 7g The seven quadrilateral grid patterns used are shown, with the numbers in brackets indicating the number of summary points. Figures 7a to 7f Due to symmetry, the grid only represents the part where 0≤x≤0.5 and 0≤y≤0.5, and “○” is the position where x=y=0.5. Figure 7c to Figure 7e There are low-quality cells with an internal angle of π in the grid. Figure 7g The grid has density. Figure 8 shows the numerical error in simulation problem 1, which is the problem of enforcing boundary conditions. Fig. 9The numerical error in simulation problem 2, i.e., 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), (c), (d), (e), (f), and (g) represent the numerical errors in the simulation problem 2, i.e., the mixed boundary condition problem. Figure 7a to Figure 7g The error of the vertical axis is defined by equation (32). For a square grid, since the correction term in equation (14) is zero, the two methods obtain the same result. Figure 7d The two methods also obtained the same results under basic boundary conditions. In other cases, it can be confirmed that the analytical accuracy of CSFE is greatly improved. In addition, the mesh dependency problem is also improved, even in the case of internal angles of π( Figure 7d ), the nodes are very close ( Figure 7f ) and other extremely low-quality units that are almost deformed into triangles, the simulation accuracy is still significantly improved.

[0232] In summary, for the finite element simulation of the transient heat conduction equation using quadrilateral units, this specification proposes a new method, which introduces a correction term of the virtual node domain based on the second-order accuracy heat flux density to correct the transient term and the source term to ensure that the conservation law is satisfied at the node level. Theory shows that the introduction of the correction term not only ensures that the conservation law is satisfied at the node level, but also satisfies the conservation law at the unit level, and effectively improves the convergence of the finite element discretization to the original differential equation. The numerical simulation results using a variety of grid modes show that compared with the traditional Galerkin finite element method, this method has significantly improved the accuracy of the simulation in an all-round way. At the same time, the grid mode dependency problem is improved, which greatly improves the credibility of the simulation results.

[0233] The above-described embodiments are used to illustrate the present invention, not to limit the present invention, so changes in the exemplified values ​​or replacement of equivalent elements should still fall within the scope of the present invention.

Claims

1. A transient heat conduction simulation method of quadrilateral elements based on virtual node domains, characterized in that: include: Select a simulation area and divide the simulation area into quadrilateral grids; Set the initial conditions, t 0 =0, n=1; Generate a cell matrix: Where Δt is the time step, t is the time, T n represents the temperature at the nth moment, T n-1 represents the temperature at the n-1th moment, λ represents the thermal conductivity, W P represents the Galerkin weight function, ξ and η are the horizontal and vertical coordinates of the isoparametric unit respectively; represents the integral over isoparametric units; is the Jacobian determinant; the subscript O indicates the value at the geometric center of the unit, Vn P is the area of ​​the virtual node domain of the target node P in the unit, is the time parameter, represents the partial derivative of x, It means to find the partial derivative of y, x and y represent the plane 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; Substitute the temperature at the previous moment into Substituting the source term Q into Generate equivalent nodal loads and generate a system of simultaneous equations: Among them, Q n is the source term at the nth moment, Q n-1 is the source term at the n-1th moment; Import the boundary conditions and solve the simultaneous equations to obtain the temperature T 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 method for simulating transient heat conduction of quadrilateral units based on virtual node domain according to claim 1, 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 target node P considered; is a transient term, Ω represents the area around the node P consisting of several units; ∫ Ω dΩ represents the integral over the region Ω; the integral over Ω is the sum of the integrals of the surrounding units, and equation (2) can be rewritten as equation (3); e represents the unit around the node P, ∫ e de represents the integral on the element e, Σ represents the sum of the integrals of each element; for the element integral on the left side of equation (3), the integration by parts theorem and Gauss-Green theorem are applied to obtain equation (4) which is the discretization of the Galerkin finite element method; For quadrilateral elements, the integral of each element in equation (4) becomes the unified equation (5); Use equation (14) to replace equation (5) of the Galerkin method; Vn P is the virtual node domain area of ​​the target node P in the unit; It is theoretically proved that the sum of the areas of the virtual node domains of all nodes of the same unit satisfies equation (15) and the conservation law is satisfied at the unit level; Where S represents the area of ​​the unit; After moving the terms in equation (14), it can be rewritten as equation (18); Formula (18) is valid at any time, that is, formula (19) is obtained; Substituting formula (18) into formula (19) yields formula (20); Substituting the mean value theorem of equation (21) into equation (20), the time derivative term disappears, and the finite element discretization result of the transient heat conduction equation is obtained, that is, a system of simultaneous equations is obtained.

3. The method for simulating transient heat conduction of quadrilateral units based on virtual node domain according to claim 2, characterized in that: For a parallelogram element with the height equal to the base and tilted π / 4, the element integral on the left side of equation (5) is calculated as equation (6); The node number of the target node P with an internal angle of 3π / 4 in the unit is 1, and the node number in the unit is represented by a numerical subscript; T1, T2, and T3 represent the temperature of nodes 1, 2, and 3 respectively, h is the height, and I ij represents the midpoint of the line segment connecting nodes i and j; λ(T1-T2) / h in the first term on the right side of equation (6) is the midpoint I 12 The heat flux density from node 1 to node 2 is approximated by the second order accuracy, and its coefficient h / 2 is the heat flux density from node 1 to node 2. 12 To point I 23 The first term on the right side of equation (6) is the distance of point I. 12 The heat flux at the vertical line segment I is evaluated with second-order accuracy. 12 I 23 The second term on the right side of equation (6) uses the second-order accuracy heat flux density from node 1 to node 3 at the geometric center of the unit O to evaluate the heat flux through the vertical line segment I 23 I 14 heat flow; Heat flux density from node 1 to 3 and line segment I 23 I 14 The right side of equation (6) indicates that the control line I 12 I 23 and I 23 I 14 Outflow heat flow; the area enclosed by the unit boundary and the control line is quadrilateral 1I 12 I 23 I 14 , whose area is 3 / 8 of the unit area.

4. The method for simulating transient heat conduction of quadrilateral units based on virtual node domain according to claim 2, characterized in that: For a parallelogram element with a height equal to the base and an inclination of π / 4, the element integral on the right side of equation (5) is calculated, and the contributions from the transient term and the source term are almost equally distributed to the four nodes of the element in the form of equivalent node loads. When and Q are constants inside the unit, the unit integral is as shown in equation (7); The transient term and source term assigned to node 1 both come from only 1 / 4 of the cell area.

5. The method for simulating transient heat conduction of quadrilateral units based on virtual node domain according to claim 2, characterized in that: For a quadrilateral element of any shape, the node number of the target node P in the element is set to 1, and the integral result on the left side of equation (5) is equation (9); A 1i It is a coefficient that depends only on the coordinates of the element nodes. Its first subscript 1 is due to the node number of the target node P in the unit. It is theoretically proved that formula (10) is valid. Using formula (10), formula (9) can be rewritten as formula (11); Introducing L 1i (i=2,3,4), let L 1i =[(x i -x1) 2 +(y i -y1) 2 ] 1 / 2 , that is, it is defined as the distance between node 1 and node i; λ(T1-T i ) / L 1i is the midpoint I of the line segment connecting nodes 1 and i 1i The heat flux density at the second-order accuracy from node 1 to i at the position is calculated based on its coefficient (-A 1i L 1i ) is the length (-A 1i L 1i ) that passes through the midpoint I 1i And it is perpendicular to the heat flux, that is, perpendicular to the line segment connecting nodes 1 and i; multiplying the two, the right side of formula (11) becomes the passing line segment (-A) evaluated with second-order accuracy 1i L 1i ) heat flow; in the unit, the line segment (-A 12 L 12 )、(-A 13 L 13 )、(-A 14 L 14 ) as the base, with node 1 as the vertex, and their areas are respectively expressed as Vn 12 、Vn 13 、Vn 14 To express, the areas of these three triangles are calculated by formula (12); VN P =Vn1=Vn 12 +Vn 13 +Vn 14 (13); The area of ​​the virtual node domain Vn1 of node 1 is the sum of the areas of these three triangles. Since node P and node 1 are identical, Vn P It is also calculated by formula (13).

6. The method for simulating transient heat conduction of quadrilateral units based on virtual node domain according to claim 5, characterized in that: For square and rectangular cells, Vn P =S / 4, the correction term is zero; for the parallelogram unit with the height equal to the base and tilted π / 4, Vn1 = 3S / 8 is obtained from formula (13), which is the same as the quadrilateral 1I 12 I 23 I 14 coincidence; the right side of formula (5), that is, the first term in the brackets on the right side of formula (14) The contribution of the cell integral to the target node is considered to come from the cell area S divided by W P The weighted transient term and source term; the correction term introduced to the Galerkin finite element method, that is, the second term in the bracket on the right side of equation (14) Its area is from the virtual node domain Vn P Subtracted Then multiply it by the transient term and source term at the geometric center O; when the transient phase and source terms are constants in the unit, the square brackets on the right side of equation (14) become That is, the system area considering transient terms and source terms is equal to the system area considering entry and exit.

7. The method for simulating transient heat conduction of quadrilateral units based on virtual node domain according to claim 5, characterized in that: The grid consisting of eight parallelograms with a height of h, a height equal to the base and a tilt of π / 4 is integrated using equation (5) to obtain equation (16); Here, the transient term and source term are assumed to be constants in the grid, the node P is the center point of the grid, and the subscripts -x, +x, -y, and +y represent the nodes located to the left, right, bottom, and top of the point P, respectively, at a distance h from the point P. The left side of equation (16) is the discretization of the second-order accuracy of the second-order differential term in equation (1), and the transient term and source term on the right side are twice the corresponding terms in equation (1). The virtual node domain Vn2=S / 8 of node 2 is obtained from equation (13), and equation (14) after introducing the correction term is used to obtain equation (17) that converges to the original differential equation.

8. The method for simulating transient heat conduction of quadrilateral units based on virtual node domain according to claim 2, characterized in that: When the simulation domain is a square area, 0≤x≤1, 0≤y≤1, set λ=1; for the numerical error in the forced boundary condition problem, the source term is taken as equation (23), the initial condition is taken as equation (24), and the forced boundary condition is taken as equation (25), and the analytical solution is obtained as equation (26); T 0 =sin(πx)sin(πy) (24); T=0(x=0,1,ory=0,1)(25); The numerical value error is evaluated using formula (32); Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the numerical 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; the time step is Δt = 0.005, Ending time t end ＝1。 9. The method for simulating transient heat conduction of quadrilateral units based on virtual node domain according to claim 2, characterized in that: When the simulation domain is a square area, 0≤x≤1, 0≤y≤1, set λ=1; for the numerical error in the mixed boundary condition problem, take equation (27) as the source term, take equation (28) as the initial condition, take equation (29) in the mixed boundary condition as the forced boundary condition, and take equation (30) as the natural boundary condition to obtain the analytical solution as equation (31); In formula (30), n represents the unit external normal direction; The numerical value error is evaluated using formula (32); Subscripts Ai, Ni, and Amax represent the analytical solution of node i, the numerical 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; the time step is Δt = 0.005, Ending time t end ＝1。 10. A quadrilateral unit transient heat conduction simulation system based on virtual node domain, characterized in that: The method for realizing the transient heat conduction simulation of quadrilateral units based on virtual node domains according to any one of claims 1 to 9, wherein the transient heat conduction simulation system of quadrilateral units based on virtual node domains comprises: A partitioning module is used to select a simulation area and divide the simulation area into quadrilateral grids; Setting module, used to set initial conditions, t 0 =0, n=1; Unit matrix generation module, used to generate unit matrix: Where Δt is the time step, t is the time, T n represents the temperature at the nth moment, T n-1 represents the temperature at the n-1th moment, λ represents the thermal conductivity, W P represents the Galerkin weight function, ξ and η are the horizontal and vertical coordinates of the isoparametric unit respectively; represents the integral over isoparametric units; is the Jacobian determinant; the subscript O indicates the value at the geometric center of the unit, Vn P is the area of ​​the virtual node domain of the target node P in the unit, is the time parameter, represents the partial derivative of x, It means to find the partial derivative of y, x and y represent the plane coordinates; The overall matrix generation module is used to generate the overall matrix: 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 module for generating simultaneous equations is used to substitute the temperature at the previous moment into Substituting the source term Q into Generate equivalent nodal loads and generate a system of simultaneous equations: Among them, Q n is the source term at the nth moment, Q n-1 is the source term at the n-1th moment; The simultaneous equations solving module is used to import boundary conditions, solve the simultaneous equations, and obtain the temperature T n ; End module, used to determine whether it is finished by time. If it is not finished, n=n+1 and continue to advance the time until it is finished.

Images