Two-dimensional nonlinear steady-state temperature field fast algorithm of improved region decomposition method

Abstract

The invention discloses a two-dimensional nonlinear steady-state temperature field fast algorithm of an improved regional decomposition method, which is characterized in that a solution domain is divided into a plurality of non-overlapped sub-regions, the nonlinear physical field distribution of each sub-region can be calculated in parallel, and the problem that a large global matrix needs to be constructed and stored in a traditional finite element method is avoided. Compared with an overlapping type region decomposition method, the method has the advantages that the sub-domain division process is simple, only sub-domain boundary nodes serve as overlapping parts, and the calculation scale is reduced. Compared with a non-overlapping region decomposition method, the method depends on Dirichlet boundary condition interaction information, and the iteration convergence speed is higher. According to the method, the finite element model of the nonlinear steady-state temperature field is solved by using an improved region decomposition method, the steady-state temperature field distribution of each sub-region can be calculated in parallel through a multi-thread processor of a calculation platform, and compared with finite element software, the calculation speed for solving the nonlinear steady-state temperature field distribution of the contactor can be increased by 3.38 times at most.

Description

Technical Field

[0001] The present invention belongs to the field of electrical numerical calculation, and relates to a two-dimensional nonlinear steady-state temperature field solution algorithm, and specifically to a method for implementing finite element parallel accelerated solution of a two-dimensional nonlinear steady-state temperature field by using an improved regional decomposition method to divide a solution domain into multiple sub-regions for parallel calculation. Background Art

[0002] In the field of modern engineering, the coupling problem of electromagnetic field and temperature field 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 field and electromagnetic field distribution of electrical appliances. At the same time, the complexity of these problems requires that the numerical analysis method must not only have high precision but also high efficiency. The finite element method is the most commonly used numerical calculation method in industrial design. It is adopted by many commercial simulation software and is widely used. However, the internal structure of electrical appliances is complex, and when calculating their nonlinear temperature field distribution, the traditional finite element method has a technical bottleneck in calculation efficiency.

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

[0004] The present invention provides a fast algorithm for two-dimensional nonlinear steady-state temperature field of an improved regional decomposition method. The algorithm can perform parallel calculations on the required problems, and is used to solve the problems of long solution time and low efficiency when solving finite element nonlinear problems using the Newton iteration method.

[0005] The objective of the present invention is achieved through the following technical solutions:

[0006] A fast algorithm for two-dimensional nonlinear steady-state temperature field based on 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 a cylindrical coordinate system according to the cross-sectional dimensions of the research object;

[0008] Step 2: Use triangular elements to discretize the solution domain to obtain a 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 1 to M respectively.

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

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

[0011] The governing equation is:

[0012]

[0013] Among them, P r is the heat source's heating power, λ is the material's thermal conductivity, is the temperature gradient, is the vector differential operator;

[0014] The boundary conditions are:

[0015] Γ 1 :T=298.15K

[0016]

[0017] Among them, Γ 1 represents the temperature of the model boundary, Γ 2 It represents the rate of change of temperature T along the outer normal direction of the boundary;

[0018] Step 5: Calculate the unit equation of each triangular unit:

[0019]

[0020] Among them, T 1 、T 2 、T 3 is the temperature of the three vertices of the triangular unit, Δ is the area of ​​the triangle, r i 、z i are the horizontal and vertical coordinates of the three vertices of the triangle 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: i = 1, m = 2, n = 3; i = 2, m = 3, n = 1; i = 3, m = 1, n = 2;

[0021] Step 6: According to the partition information obtained in step 3, assemble all triangular unit equations into a sub-region global matrix to obtain the global matrices PartS and PartF of each sub-region, where PartS is a square matrix with a dimension equal to the number of nodes in the sub-region, and PartF is a vector with a dimension equal to the number of nodes in the sub-region;

[0022] Step 7: According to the partition information obtained in step 3, retrieve all the boundary nodes of different sub-areas, record these nodes as additional nodes, number them starting from 1 and record the unit number and sub-area 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 =T 1 -T r

[0025] Among them, T q is the temperature value of the node in the subdomain during the iteration process, and q is the number of iterations of the subdomain boundary conditions, that is, Figure 2 The number of times steps 8 to 13 are executed in a loop, T 1 is the additional node temperature, T r is the temperature of the boundary node at the end of the previous iteration;

[0026] Step 9: Update the unit thermal conductivity λ according to the unit temperature at the end of the previous iteration and the material properties of the area to which the unit belongs e and the derivative of the element thermal conductivity with respect to the element temperature Among them: the unit temperature in 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,PartT 1 、PartT 2、 Part T 3 is the temperature of the nodes numbered 1, 2, and 3 in the sub-region element, PartT e is the temperature of the sub-region unit, Partλ e is the thermal conductivity of the sub-region unit;

[0030] Step 11: Use the Newton-Raphson method to solve the nonlinear equations in each sub-region:

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

[0032] [Part J k ][Part T k+1 ]=[PartF]+([PartJ k ]-[PartS k])[Part T k ]

[0033] Among them, PartT is the temperature vector of the sub-region node, and its dimension is the number of sub-region nodes; k represents the number of iterations when solving the nonlinear equation system, that is, Figure 2 The number of times steps 9 to 11 are executed in a loop;

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

[0035] T r =T 1 ′-T q

[0036] Among them, T 1 ' 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 is the global node temperature vector of the previous iteration step. In the first iteration, T old The element size in is the ambient temperature;

[0040] If the convergence error ε is less than the expected value, the temperature field calculation is completed, and the temperature field distribution cloud diagram 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.

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

[0042] 1. The present invention divides the solution domain into multiple non-overlapping sub-regions, and can calculate the nonlinear physical field distribution of each sub-region in parallel, thus avoiding the problem of needing to construct and store large global matrices in traditional finite element methods.

[0043] 2. Compared with the overlapping region decomposition method, the subdomain division process of the present invention is simple and only the subdomain boundary nodes are regarded as overlapping parts, which reduces the calculation scale.

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

[0045] 4. The present invention uses an improved regional 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 the finite element software, the calculation speed of the present invention for solving the nonlinear steady-state temperature field distribution of the contactor can be increased by up to 3.38 times. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 Schematic diagram of the improved domain decomposition method, (a) the region to be solved, (b) the split intersection point;

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

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

[0049] Figure 4 This is a comparison diagram of temperature field distribution between the improved regional decomposition method and the finite element software;

[0050] Figure 5 It is a contactor temperature rise measurement platform;

[0051] Figure 6 Comparison of calculation time between improved domain decomposition method and finite element software. DETAILED DESCRIPTION

[0052] The technical solution of the present invention is further described below in conjunction with the accompanying drawings, but is not limited thereto. Any modification or equivalent replacement of the technical solution of the present invention without departing from the spirit and scope of the technical solution of the present invention should be included in the protection scope of the present invention.

[0053] The present invention provides a fast algorithm for two-dimensional nonlinear steady-state temperature field of an improved domain decomposition method, which is suitable for calculating the nonlinear steady-state temperature field distribution of electromechanical devices modeled in a plane coordinate system or a cylindrical coordinate system, such as relays, contactors, etc. Figure 1 and Figure 2 As shown, the method comprises the following steps:

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

[0055] Step 2: Use triangular elements to discretize the solution domain to obtain a 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 1 to N and 1 to M respectively.

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

[0057] Step 4: List the control equations and boundary conditions of the two-dimensional nonlinear steady-state temperature field to obtain a set of differential equations, whose control equations are:

[0058]

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

[0060] The boundary conditions are:

[0061] Γ 1 :T=298.15K

[0062]

[0063] Among them, Γ 1 represents the temperature of the model boundary, Γ 2 It represents the rate of change of temperature T along the outer normal direction of the boundary.

[0064] Step 5: Calculate the unit equation of each triangular unit:

[0065]

[0066] Among them, T 1 、T 2 、T 3 is the temperature of the three vertices of the triangular unit, Δ is the area of ​​the triangle, r i 、z i are the horizontal and vertical coordinates of the three vertices of the triangle unit, i = 1, 2, 3, b i =z m -z n , c i =r n -r m , the relationship between i, m and n is shown in Table 1:

[0067] Table 1 Relationship between unit equation coefficient numbers and vertex numbers

[0068]

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

[0070] Step 7: According to the partition information obtained in step 3, retrieve all the boundary nodes of different sub-areas, record these nodes as additional nodes, number them starting from 1 and record the unit number and sub-area 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 =T 1 -T r

[0073] Among them, T q is the temperature value of the node in the subdomain during the iteration process, q is the number of iterations of the subdomain boundary condition, T 1 is the additional node temperature, T r is the temperature of the boundary node at the end of the previous iteration.

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

[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,PartT 1 、PartT 2、 Part T 3 is the temperature of the nodes numbered 1, 2, and 3 in the sub-region element, PartT e is the temperature of the sub-region unit, Partλ e is the thermal conductivity of the sub-region unit.

[0079] Step 11: Use the Newton-Raphson method to solve the nonlinear equations in each sub-region:

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

[0081] [Part J k ][Part T k+1 ]=[PartF]+([PartJ k ]-[PartS k ])[Part T k ]

[0082] Among them, PartT is the temperature vector of the sub-region node, 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 equations, update the temperature of the boundary nodes between subdomains according to the following formula:

[0084] T r =T 1 ′-T q

[0085] Among them, T 1 ' is the temperature of the subdomain boundary node 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 is the global node temperature vector of the previous iteration step. In the first iteration, T old The element size in is the ambient temperature.

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

[0090] Example:

[0091] This embodiment is aimed at Figure 3 The nonlinear steady-state temperature field of a certain type of axisymmetric contactor is calculated. The method proposed in the present invention and the 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 cloud map is drawn on Figure 4After working for a long time, the heat generation and heat dissipation inside the contactor reach a balance. At this time, the armature temperature rise calculated by the improved domain decomposition method and the finite element software are 60.24K and 60.62K respectively. The calculation error of the temperature rise by the improved domain decomposition method (DDM) and the finite element software is 0.76%. To further verify the accuracy of DDM in calculating the temperature field distribution, the following Figure 5 The contactor temperature rise measurement platform shown in the figure. The temperature collector sensor is fixed on the armature test point, and the change of the test point temperature over time after the coil is powered on is recorded. After the coil is powered on, the internal temperature of the contactor rises rapidly. As the power-on time continues to increase, the internal temperature rise rate of the contactor gradually slows down until it reaches a stable value. The difference between the temperature at this time and the initial temperature is the steady-state temperature rise of the contactor. The steady-state temperature rise of the contactor measured in the experiment is 58.1K. The error between the improved regional decomposition method and the experiment is 3.68%, which proves that the DDM method is effective in calculating the steady-state temperature rise of the contactor.

[0092] pass Figure 6 As can be seen from Table 2, when the mesh is the same and the convergence error is the same, DDM can achieve a speedup ratio of 1.89 when solving the nonlinear temperature field and static magnetic field of the contactor compared with the finite element software, even if it is divided into only two sub-areas. When the solution domain is divided into 16 sub-areas, the calculation speed of DDM solving the nonlinear temperature field distribution of the contactor is increased by up to 3.38 times compared with the finite element software. With the increase in the scale of the problem to be solved and the improvement of the parallel capability of the computing platform, the advantage of DDM in terms of calculation speed over the finite element software will be further increased.

[0093] Table 2 Speedup ratio of improved domain decomposition method

[0094]

Claims

1. A fast algorithm for two-dimensional nonlinear steady-state temperature field based on an improved domain decomposition method, characterized in that The method comprises the following steps: Step 1: Abstract the problem to be solved and establish a geometric model in a plane coordinate system or a cylindrical coordinate system according to the cross-sectional dimensions of the research object; Step 2: Use triangular elements to discretize the solution domain to obtain a 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 1 to M respectively. Step 3: Select the number of sub-regions to be divided according to the requirements, use the METIS library function to implement the multi-region division of the finite element model, and obtain the mesh information after the sub-region division; Step 4, list the control equations and boundary conditions of the two-dimensional nonlinear steady-state temperature field to obtain a set of differential equations; Step 5: Calculate the unit equation of each triangular unit: Among them, T1, T2, and T3 are the temperatures of the three vertices of the triangular unit, Δ is the area of ​​the triangle, and r i 、z i are the horizontal and vertical coordinates of the three vertices of the triangle 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: i = 1, m = 2, n = 3; i = 2, m = 3, n = 1; i = 3, m = 1, n = 2; Step 6: According to the partition information obtained in step 3, assemble all triangular unit equations into a sub-region global matrix to obtain the global matrices PartS and PartF of each sub-region, where PartS is a square matrix with a dimension equal to the number of nodes in the sub-region, and PartF is a vector with a dimension equal to the number of nodes in the sub-region; Step 7: According to the partition information obtained in step 3, retrieve all the boundary nodes of different sub-areas, record these nodes as additional nodes, number them starting from 1 and record the unit number and sub-area 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: T q =T1-T r Among them, T q is the temperature value of the node in the subdomain during the iteration process, q is the number of iterations of the subdomain boundary condition, T1 is the temperature value of the additional node, T r is the temperature of the boundary node at the end of the previous iteration; Step 9: Update the unit thermal conductivity λ according to the unit temperature at the end of the previous iteration and the material properties of the area to which the unit belongs e and the derivative of the element thermal conductivity with respect to the element temperature Among them: the unit temperature in the first iteration is the ambient temperature; Step 10: Calculate the Jacobian matrix PartJ of the sub-region coefficient matrix PartS according to the following formula: in, i,j=1,2,3,PartT1、PartT 2、 PartT3 is the temperature of the nodes numbered 1, 2, and 3 in the sub-region element. e is the temperature of the sub-region unit, Partλ e is the thermal conductivity of the sub-region unit; Step 11: Use the Newton-Raphson method to solve the nonlinear equations in each sub-region: [PartS][PartT]=[PartF] [PartJ k ][PartT k+1 ]=[PartF]+([PartJ k ]-[PartS k ])[PartT k ] Where PartT is the temperature vector of the sub-region node, and its dimension is the number of sub-region nodes; k represents the number of iterations when solving the nonlinear equation system; Step 12: After solving the equations, update the temperature of the boundary nodes between subdomains according to the following formula: T r =T1′-T q Among them, T1' is the temperature of the boundary node of the subdomain in the current iteration step; 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: Where T is the global node temperature vector calculated in the current iteration step, T old is the global node temperature vector of the previous iteration step. In the first iteration, T old The element size in is the ambient temperature; If the convergence error ε is less than the expected value, the temperature field calculation is completed, and the temperature field distribution cloud diagram 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.

2. The fast algorithm of two-dimensional nonlinear steady-state temperature field of the improved domain decomposition method according to claim 1 is characterized in that The control equation is: Among them, P r is the heat source's heating power, λ is the material's thermal conductivity, is the temperature gradient, is a vector differential operator.

3. The fast algorithm of two-dimensional nonlinear steady-state temperature field of the improved domain decomposition method according to claim 1 is characterized in that The boundary conditions are: Γ1:T=298.15K C2: Among them, Γ1 represents the temperature of the model boundary, and Γ2 represents the rate of change of temperature T along the outer normal direction of the boundary.

Images