A Fast Algorithm for Two-Dimensional Nonlinear Steady-State Temperature Fields Based on an Improved Domain Decomposition Method

Abstract

This invention discloses a fast algorithm for solving two-dimensional nonlinear steady-state temperature fields using an improved domain decomposition method. This algorithm divides the solution domain into multiple non-overlapping sub-regions, enabling parallel computation of the nonlinear physical field distribution in each sub-region, thus avoiding the need to construct and store large global matrices in traditional finite element methods. Compared to overlapping domain decomposition methods, this invention simplifies the sub-domain partitioning process and only considers the sub-domain boundary nodes as overlapping parts, reducing the computational scale. Compared to non-overlapping domain decomposition methods, this invention relies on Dirichlet boundary condition interaction information, resulting in faster iterative convergence. This invention uses the improved domain decomposition method to solve the finite element model of the nonlinear steady-state temperature field. The steady-state temperature field distribution in each sub-region can be computed in parallel using the multi-threaded processor of the computing platform. Compared to finite element software, this invention can improve the computational speed for solving the nonlinear steady-state temperature field distribution of a contactor by up to 3.38 times.

Description

Technical Field

[0001] This invention belongs to the field of electrical numerical calculation and relates to a solution algorithm for a two-dimensional nonlinear steady-state temperature field. Specifically, it relates to a method for accelerating the parallel finite element solution of a two-dimensional nonlinear steady-state temperature field by dividing the solution domain into multiple sub-regions for parallel computation using an improved domain decomposition method. Background Technology

[0002] In modern engineering, the coupling of electromagnetic and temperature fields is crucial for electromagnetic compatibility analysis and reliability optimization design of electrical appliances. The key to solving these problems lies in accurately calculating the temperature and electromagnetic field distributions of the appliances. Furthermore, the complexity of these problems demands that numerical analysis methods possess not only high precision but also efficient computational capabilities. The finite element method (FEM) is the most commonly used numerical calculation method in industrial design, widely adopted by numerous commercial simulation software programs. However, due to the complex internal structure of electrical appliances, traditional FEM methods face computational efficiency bottlenecks when calculating their nonlinear temperature field distributions.

[0003] Therefore, it is necessary to study a parallel computing method to solve the problems of long solution time and low efficiency caused by the traditional finite element method in solving nonlinear temperature field problems. Summary of the Invention

[0004] This invention provides a fast algorithm for two-dimensional nonlinear steady-state temperature fields using an improved domain decomposition method. This algorithm can solve the problem in parallel, addressing the issues of long solution time and low efficiency inherent in Newton's iteration method for solving finite element nonlinear problems.

[0005] The objective of this invention is achieved through the following technical solution:

[0006] A fast algorithm for two-dimensional nonlinear steady-state temperature fields using an improved domain decomposition method includes the following steps:

[0007] Step 1: Abstract the problem to be solved and establish a geometric model in a plane coordinate system or cylindrical coordinate system according to the cross-sectional dimensions of the object under study;

[0008] Step 2: Discretize and mesh the solution domain using triangular elements to obtain the finite element network of the model to be solved. The total number of triangular elements in the finite element network is N, and the total number of nodes is M. The triangular elements and nodes are numbered from 1 to N and from 1 to M, respectively.

[0009] Step 3: Select the number of sub-regions according to the requirements, and use METIS library functions to implement multi-region division of the finite element model to obtain the mesh information after sub-region division;

[0010] Step 4: List the governing equations and boundary conditions of the two-dimensional nonlinear steady-state temperature field, resulting in a set of differential equations, where:

[0011] The governing equations are:

[0012]

[0013] Among them, P r λ is the heating power of the heat source, and λ is the thermal conductivity of the material. For the temperature gradient, It is a vector differential operator;

[0014] The boundary conditions are:

[0015] Γ1:T=298.15K

[0016]

[0017] Where Γ1 represents the temperature at the model boundary, and Γ2 represents the rate of change of temperature T along the outward normal direction of the boundary;

[0018] Step 5: Calculate the element equation for each triangular element:

[0019]

[0020] Where T1, T2, and T3 are the temperatures at the three vertices of the triangular unit, Δ is the area of ​​the triangle, and r i z i Let i be the x and y coordinates of the three vertices of the triangular unit, i = 1, 2, 3, b i =z m -z n c i =r n -r m The relationship between i and m and n is as follows: i = 1, m = 2, n = 3; i = 2, m = 3, n = 1; i = 3, m = 1, n = 2;

[0021] Step 6: Based on the partitioning information obtained in Step 3, assemble all the triangular unit equations into a sub-region global matrix to obtain the global matrix PartS and PartF for each sub-region. PartS is a square matrix with dimension equal to the number of nodes in the sub-region, and PartF is a vector with dimension equal to the number of nodes in the sub-region.

[0022] Step 7: Based on the partition information obtained in Step 3, retrieve all boundary nodes of different sub-regions, record these nodes as additional nodes, number them starting from 1, and record the unit number and sub-region number to which these nodes belong.

[0023] Step 8: Assign Dirichlet boundary conditions to the sub-region boundary nodes retrieved in Step 7 according to the following formula:

[0024] T q =T1-T r

[0025] Among them, T q Let q be the temperature value of the node within the subdomain during the iteration process, and q be the number of iterations for the subdomain boundary conditions. Figure 2 The number of times steps 8 to 13 are executed in the loop, where T1 is the additional node temperature value, and T... r The temperature of the boundary node at the end of the previous iteration step;

[0026] Step 9: Update the element thermal conductivity λ based on the element temperature at the end of the previous iteration and the material properties of the region to which the element belongs. e and the derivative of the unit thermal conductivity with respect to the unit temperature Wherein: the cell temperature during the first iteration is the ambient temperature;

[0027] Step 10: Calculate the Jacobian matrix PartJ of the sub-region coefficient matrix PartS according to the following formula:

[0028]

[0029] in, i,j = 1,2,3,PartT1,PartT 2、 PartT3 represents the temperature of nodes numbered 1, 2, and 3 in the sub-region. e For the temperature of the sub-region unit, Partλ e The thermal conductivity of the sub-region unit;

[0030] Step 11: Solve the nonlinear equations for each subregion using the Newton-Raphson method:

[0031] [PartS][PartT] = [PartF]

[0032] Part J k ][PartT k+1 ] = [PartF] + ([PartJ] k Parts k Part T k ]

[0033] Where PartT is the temperature vector of the sub-region nodes, and its dimension is the number of nodes in the sub-region; k represents the number of iterations when solving the nonlinear equation system, i.e. Figure 2 The number of times steps 9 through 11 are executed in a loop;

[0034] Step 12: After solving the system of equations, update the temperature of the boundary nodes between subdomains according to the following formula:

[0035] T r =T1′-T q

[0036] Where T1' is the temperature of the boundary node of the subdomain in the current iteration step;

[0037] Step 13: Assemble the node temperature vector PartT of each sub-region into the global node temperature vector T, and calculate the convergence error of the current iteration step according to the following formula:

[0038]

[0039] Where T is the global node temperature vector calculated in the current iteration step, T old Let T be the global node temperature vector from the previous iteration. In the first iteration... old The size of the element in the image is the ambient temperature.

[0040] If the convergence error ε is less than the expected value, the temperature field calculation is completed, and a temperature field distribution cloud map is drawn based on the node temperature T; if the convergence error ε is not less than the expected value, the next iteration is initiated, and steps 8 to 13 are repeated.

[0041] Compared with the prior art, the present invention has the following advantages:

[0042] 1. This invention divides the solution domain into multiple non-overlapping sub-regions, enabling parallel computation of the nonlinear physical field distribution in each sub-region, thus avoiding the problem of constructing and storing large global matrices in traditional finite element methods.

[0043] 2. Compared with the overlapping region decomposition method, the subdomain partitioning process of this invention is simple and only the boundary nodes of the subdomains are considered as overlapping parts, which reduces the computational scale.

[0044] 3. Compared with the non-overlapping domain decomposition method, the present invention relies on the interaction information of Dirichlet boundary conditions, resulting in faster iterative convergence speed.

[0045] 4. This invention uses an improved domain decomposition method to solve the finite element model of the nonlinear steady-state temperature field. The steady-state temperature field distribution of each sub-region can be calculated in parallel by the multi-threaded processor of the computing platform. Compared with finite element software, the calculation speed of this invention to solve the nonlinear steady-state temperature field distribution of the contactor can be increased by up to 3.38 times. Attached Figure Description

[0046] Figure 1 This is a schematic diagram of the improved domain decomposition method. (a) The domain to be solved, and (b) the decomposition boundary points.

[0047] Figure 2 Flowchart for solving the nonlinear steady-state temperature field using the improved domain decomposition method;

[0048] Figure 3 The electromagnetic mechanism of a certain type of contactor and its finite element model;

[0049] Figure 4 A comparison of temperature field distribution between the improved zone decomposition method and finite element software;

[0050] Figure 5 For contactor temperature rise measurement platform;

[0051] Figure 6 This study compares the computation time of the improved domain decomposition method with that of finite element software. Detailed Implementation

[0052] The technical solution of the present invention will be further described below with reference to the accompanying drawings, but it is not limited thereto. Any modifications or equivalent substitutions to the technical solution of the present invention that do not depart from the spirit and scope of the technical solution of the present invention should be covered within the protection scope of the present invention.

[0053] This invention provides a fast algorithm for two-dimensional nonlinear steady-state temperature fields using an improved domain decomposition method. This algorithm is suitable for calculating the nonlinear steady-state temperature field distribution of electromechanical components modeled in planar or cylindrical coordinate systems, such as relays and contactors. Figure 1 and Figure 2 As shown, the method includes the following steps:

[0054] Step 1: Abstract the problem to be solved and establish a geometric model in a plane coordinate system or cylindrical coordinate system according to the cross-sectional dimensions of the object under study.

[0055] Step 2: Discretize the solution domain using triangular elements to obtain the finite element network of the model to be solved. The total number of triangular elements in the finite element network is N, and the total number of nodes is M. The triangular elements and nodes are numbered from 1 to N and from 1 to M, respectively.

[0056] Step 3: Select the number of sub-regions to be divided according to the requirements, and use METIS library functions to implement multi-region division of the finite element model to obtain the mesh information after dividing the sub-regions.

[0057] Step 4: List the governing equations and boundary conditions of the two-dimensional nonlinear steady-state temperature field, resulting in a set of differential equations. The governing equations are as follows:

[0058]

[0059] Among them, P r λ is the heating power of the heat source, and λ is the thermal conductivity of the material. For the temperature gradient, It is a vector differential operator.

[0060] The boundary conditions are:

[0061] Γ1:T=298.15K

[0062]

[0063] Where Γ1 represents the temperature of the model boundary, and Γ2 represents the rate of change of temperature T along the outward normal direction of the boundary.

[0064] Step 5: Calculate the element equation for each triangular element:

[0065]

[0066] Where T1, T2, and T3 are the temperatures at the three vertices of the triangular unit, Δ is the area of ​​the triangle, and r i z i Let i be the x and y coordinates of the three vertices of the triangular unit, i = 1, 2, 3, b i =z m -z n c i =r n -r m The relationship between i and m and n is shown in Table 1:

[0067] Table 1. Relationship between coefficient numbers and vertex numbers in unit equations

[0068]

[0069] Step 6: Based on the partitioning information obtained in Step 3, assemble all the triangular unit equations into a sub-region global matrix to obtain the global matrices PartS and PartF for each sub-region. PartS is a square matrix with dimension equal to the number of nodes in the sub-region, and PartF is a vector with dimension equal to the number of nodes in the sub-region.

[0070] Step 7: Based on the partition information obtained in Step 3, retrieve all boundary nodes of different sub-regions, record these nodes as additional nodes, number them starting from 1, and record the unit number and sub-region number to which these nodes belong.

[0071] Step 8: Assign Dirichlet boundary conditions to the sub-region boundary nodes retrieved in Step 7 according to the following formula:

[0072] T q =T1-T r

[0073] Among them, T qLet T1 be the temperature value of the node within the subdomain during the iteration process, q be the number of iterations for the subdomain boundary conditions, T1 be the temperature value of the additional node, and T be the temperature value of the node within the subdomain. r This represents the temperature of the boundary node at the end of the previous iteration step.

[0074] Step 9: Update the element thermal conductivity λ based on the element temperature at the end of the previous iteration (ambient temperature during the first iteration) and the material properties of the region to which the element belongs. e and the derivative of the unit thermal conductivity with respect to the unit temperature

[0075]

[0076] Step 10: Calculate the Jacobian matrix PartJ of the sub-region coefficient matrix PartS according to the following formula:

[0077]

[0078] in, i,j = 1,2,3,PartT1,PartT 2、 PartT3 represents the temperature of nodes numbered 1, 2, and 3 in the sub-region. e For the temperature of the sub-region unit, Partλ e denoted as the thermal conductivity of the sub-region unit.

[0079] Step 11: Solve the nonlinear equations for each subregion using the Newton-Raphson method:

[0080] [PartS][PartT] = [PartF]

[0081] Part J k ][PartT k+1 ] = [PartF] + ([PartJ] k Parts k Part T k ]

[0082] Where PartT is the temperature vector of the sub-region nodes, and its dimension is the number of sub-region nodes; k represents the number of iterations when solving the nonlinear equation system.

[0083] Step 12: After solving the system of equations, update the temperature of the boundary nodes between subdomains according to the following formula:

[0084] T r =T1′-T q

[0085] Where T1' is the temperature of the boundary node of the subdomain in the current iteration step.

[0086] Step 13: Assemble the node temperature vector PartT of each sub-region into the global node temperature vector T, and calculate the convergence error of the current iteration step according to the following formula:

[0087]

[0088] Where T is the global node temperature vector calculated in the current iteration step, T old Let T be the global node temperature vector from the previous iteration. In the first iteration... old The size of the element in the equation is the ambient temperature.

[0089] If the convergence error ε is less than the expected value (set according to the calculation accuracy requirements), the temperature field calculation is completed, and the temperature field distribution cloud map is drawn according to the node temperature T; if the convergence error ε is not less than the expected value, the next iteration is entered, and steps 8 to 13 are repeated.

[0090] Example:

[0091] This embodiment is for Figure 3 The nonlinear steady-state temperature field of a certain type of axisymmetric contactor is calculated. The method proposed in this invention and finite element software are used to solve the steady-state temperature field distribution of the two-dimensional axisymmetric finite element model of the contactor, and the temperature field contour map is plotted. Figure 4 After prolonged operation, the heat generation and dissipation inside the contactor reach equilibrium. At this point, the armature temperature rise calculated by the improved domain decomposition method (DDM) and the finite element method (FEA) is 60.24 K and 60.62 K, respectively, with a calculation error of 0.76%. To further verify the accuracy of DDM in calculating the temperature field distribution, a system was constructed as follows: Figure 5 The contactor temperature rise measurement platform is shown. The temperature sensor of the temperature acquisition device is fixed at the armature test point, and the temperature change at the test point over time is recorded after the coil is energized. After the coil is energized, the internal temperature of the contactor rises rapidly. As the energizing time continues to increase, the rate of temperature rise gradually slows down until a steady-state value is reached. The difference between this temperature and the initial temperature is the steady-state temperature rise of the contactor. The experimentally measured steady-state temperature rise of the contactor is 58.1 K. The error between the improved domain decomposition method and the experiment is 3.68%, proving that the DDM method is effective in calculating the steady-state temperature rise of the contactor.

[0092] pass Figure 6As shown in Table 2, under the same mesh and convergence error, DDM achieves a speedup of 1.89 compared to finite element software, even when the solution domain is divided into only two sub-regions, in solving the nonlinear temperature and static magnetic fields of the contactor. When the solution domain is divided into 16 sub-regions, the computational speed of DDM in solving the nonlinear temperature field distribution of the contactor is up to 3.38 times faster than that of finite element software. With the increase in the scale of the problem and the improvement of the parallel processing capabilities of the computing platform, the speed advantage of DDM over finite element software will further increase.

[0093] Table 2 Speedup Ratio of Improved Domain Decomposition Method

[0094]

Claims

1. A two-dimensional nonlinear steady-state temperature field fast calculation method of an improved domain decomposition method, characterized in that The method includes the following steps: Step 1: Abstract the problem to be solved and establish a geometric model in a plane coordinate system or cylindrical coordinate system according to the cross-sectional dimensions of the object under study; Step 2: Discretize and mesh the solution domain using triangular elements to obtain the finite element network of the model to be solved. The total number of triangular elements in the finite element network is N, and the total number of nodes is M. The triangular elements and nodes are numbered from 1 to N and from 1 to M, respectively. Step 3: Select the number of sub-regions according to the requirements, and use METIS library functions to implement multi-region division of the finite element model to obtain the mesh information after sub-region division; Step 4: List the governing equations and boundary conditions of the two-dimensional nonlinear steady-state temperature field to obtain a set of differential equations; Step 5: Calculate the element equation for each triangular element: wherein , , is the temperature of the three vertices of the triangle element, is the area of the triangle, , is the horizontal and vertical coordinates of the three vertices of the triangle element, = 1, 2, 3, , , and and are related by: = 1, = 2, = 3; = 2, = 3, = 1; = 3, = 1, = 2; Step 6: Based on the partitioning information obtained in Step 3, assemble all the triangular unit equations into a sub-region global matrix to obtain the global matrix for each sub-region. , ,in, It is a square matrix with dimension equal to the number of nodes in the sub-region. It is a vector whose dimension is the number of nodes in the sub-region; Step 7: Based on the partition information obtained in Step 3, retrieve all boundary nodes of different sub-regions, record these nodes as additional nodes, number them starting from 1, and record the unit number and sub-region number to which these nodes belong. Step 8: Assign Dirichlet boundary conditions to the sub-region boundary nodes retrieved in Step 7 according to the following formula: in, The temperature value of the node within the subdomain during the iteration process. The number of iterations for the subdomain boundary conditions. For additional node temperature values, The temperature of the boundary node at the end of the previous iteration step; Step 9: Update the element thermal conductivity based on the element temperature at the end of the previous iteration and the material properties of the region to which the element belongs. and the derivative of the unit thermal conductivity with respect to the unit temperature Where: the unit temperature during the first iteration is the ambient temperature; Step 10: Calculate the sub-region coefficient matrix according to the following formula. Jacobi matrix : in, i,j=1,2,3 , , The temperature of nodes numbered 1, 2, and 3 in the sub-region. For the temperature of the sub-region unit, The thermal conductivity of the sub-region unit; Step 11: Solve the nonlinear equations for each subregion using the Newton-Raphson method: in, This is the temperature vector of the sub-region nodes, with the dimension being the number of nodes in the sub-region. This represents the number of iterations when solving a system of nonlinear equations. Step 12: After solving the system of equations, update the temperature of the boundary nodes between subdomains according to the following formula: in, The temperature of the boundary node of the subdomain in the current iteration step; Step 13: Convert the node temperature vector of each sub-region Assembled as a global node temperature vector The convergence error of the current iteration step is calculated using the following formula: in, This is the global node temperature vector calculated in the current iteration step. The global node temperature vector from the previous iteration, during the first iteration. The size of the element in the image is the ambient temperature. If convergence error If the value is less than the expected value, the temperature field calculation is completed based on the nodal temperatures. Plot the temperature field distribution contour map; if the convergence error If the result is not less than the expected value, proceed to the next iteration and repeat steps 8 to 13.

2. The method for rapid calculation of two-dimensional nonlinear steady-state temperature field using the improved domain decomposition method according to claim 1, characterized in that... The governing equation is: in, It is the heating power of the heat source. The thermal conductivity of the material. For the temperature gradient, It is a vector differential operator.

3. The method for rapid calculation of two-dimensional nonlinear steady-state temperature field using the improved domain decomposition method according to claim 1, characterized in that... The boundary conditions are as follows: in, The temperature at the model boundary is represented. Indicates temperature The rate of change along the direction of the outward normal to the boundary.

Images