An improved arc-length method for efficient solution of nonlinear thermal boundary problems

Abstract

The present invention discloses an improved arc length method for efficiently solving nonlinear thermal boundary problems, and belongs to the field of high-performance computing simulation technology. The method first needs to establish a geometric model to be analyzed, and then the geometric model to be analyzed is meshed with any unit type, and then the boundary conditions and load conditions of the model to be analyzed are set, and finally the above conditions are input into the algorithm solver. In the solution, the size of the second norm of the tangent increment obtained in the prediction stage of adjacent incremental steps is determined. The arc length calculation method of the algorithm prediction stage is selected according to the determination result. Thereby, the algorithm solution process is accelerated under the premise of realizing reliable solution of heat transfer nonlinear problems. The present invention combines heat transfer theory, finite element method, and nonlinear algorithm control strategy, improves on the theoretical basis of the existing arc length method, and uses the excellent numerical calculation ability of computers as a carrier to solve thermal nonlinear problems, which not only improves the numerical accuracy of the calculation results, but also improves the speed of solving thermal nonlinear problems.

Description

Technical Field

[0001] The invention belongs to the technical field of high performance computing simulation, and in particular is an improved arc length method for efficiently solving nonlinear thermal boundary problems. Background Art

[0002] Heat transfer is everywhere. Heat transfer is not only a common natural phenomenon, but also widely exists in the field of engineering technology. In the energy and power, chemical and pharmaceutical, material metallurgy, machinery manufacturing, and other sectors, there are prominent problems of heat transfer.

[0003] At present, although the theory of heat transfer has been gradually improved, since most of the problems dealt with in heat transfer are nonlinear practical problems, they are generally solved by breaking down the nonlinear problem into multiple load increments, applying the traditional solution idea of ​​replacing the curve with a straight line to each single load increment, and approximating the nonlinear problem with multiple linear solution iterations. However, with the rapid development of electronic computer hardware and software technology and the continuous optimization of numerical simulation theory and algorithms. In terms of nonlinear thermal analysis, the requirements for numerical accuracy, reliability of algorithm convergence, and solution speed are constantly increasing. The solution algorithm based on the Newton iteration method has been continuously improved and upgraded, and a number of highly applicable algorithms for solving nonlinear thermal problems have emerged. In addition, by improving the control strategy of loads and temperatures, refining the load increments, and adjusting the solution speed during the solution process, the speed has been improved at the algorithm level while ensuring the numerical accuracy of the calculation results.

[0004] Therefore, the present invention establishes an improved arc length method for efficiently solving nonlinear thermal boundary problems, and utilizes the well-developed finite element numerical calculation theory, heat transfer theory, and nonlinear solution strategies to focus on solving thermal nonlinear problems. While accurately solving nonlinear thermal problems, the solution speed is improved. It can not only be applied to large-area heat transfer analysis of high-end equipment such as aerospace, but also provide a more accurate solution method for nonlinear problems in high-precision fields such as chip design and micro-nano manufacturing, which can greatly promote the optimization of thermal and thermodynamic properties of material and structural design. Summary of the invention

[0005] The present invention discloses an improved arc length method for efficiently solving nonlinear thermal boundary problems. The method combines heat transfer theory, finite element method and nonlinear algorithm control strategy, improves the existing arc length method on the theoretical basis, and uses the excellent numerical calculation ability of computers as a carrier to solve thermal nonlinear problems, which not only improves the numerical accuracy of the calculation results, but also improves the speed of solving thermal nonlinear problems.

[0006] The technical solution provided by the present invention to solve the above technical problems is: an improved arc length method for efficiently solving nonlinear thermal boundary problems, comprising the following steps:

[0007] S1. Establishing a geometric model or geometric assembly model of a thermal analysis object to be determined;

[0008] S2, spatially discretize the geometric model to be analyzed, i.e., mesh division, to generate the mesh data required for calculation;

[0009] S3, setting parameters such as boundary conditions, loads, and material thermal conductivity for the divided grid model;

[0010] S4. Set the value of the initial load factor of the first incremental step of the modified arc length method λ initial , iteration tolerance value t, load factor symbol S, total load factor value λ, number of Hammer integration points N;

[0011] S5. To solve the problem under given load boundary conditions using the modified arc length method, it is necessary to calculate the current reference load heat flow rate vector {Q ref}. Enter the load increment prediction stage to solve. Calculate the tangent conduction matrix of the unit, the thermal conductivity of each unit under the current configuration, and the tangent prediction increment solution. Determine the size of the bi-norm of the tangent increment of the adjacent load increment step, and determine the acceleration method based on the judgment result;

[0012] S6, enter the correction solution stage, first update the current external load vector {Q e}、Internal heat flow rate vector {Q a}, residual vector Secondly, determine whether the current load increment step has reached the convergence standard. After the final iteration phase is completed, determine whether the total load factor exceeds the preset maximum load factor value;

[0013] S7, complete the thermal analysis solution and collect the temperature vector results of the final temperature field;

[0014] S8, calculate the unit conduction matrix, and calculate the temperature gradient and heat flow to achieve result post-processing;

[0015] S9. Calculate the average heat flux of the node to complete the solution of the nonlinear problem of thermal analysis under the given load boundary conditions; record the total time consumed by the solver.

[0016] A further technical solution is that the specific implementation method of step S1 is:

[0017] Establish a geometric model or geometric assembly model of any size at the macro scale and establish a continuous solution domain;

[0018] A further technical solution is that the specific implementation method of step S2 is:

[0019] S21, setting of discrete parameters;

[0020] Combined with the geometric features of large-scale geometric continuous models, the global mesh seed of the geometric model and the local seeds of some characteristic geometric edges are set, and local refinement parameters are set at the characteristic edges with complex geometric features.

[0021] S22. Mathematical modeling of heat transfer problems will produce partial differential equations (PDEs). Numerical analysis, i.e., finite element method calculation, enables the partial differential equations to be approximately solved by numerical solutions of the original model domain.

[0022] K(T)=Q,Q∈R n

[0023] Where n is the total degree of freedom of the discretized system, K is the coefficient matrix of the linear equations of the system under the current boundary load problem, Q is the load vector, and the unknown quantity u can be solved using linear algebra methods;

[0024] S23, establishing a discretized geometric model;

[0025] Set meshing parameters and use meshing tools to mesh geometric models. This can realize spatial discretization of large-scale geometric models of different shapes, discretize the model into spatial continuous discrete bodies with tetrahedrons as units, and write the discretized data information of the model to local files.

[0026] A further technical solution is that the specific implementation method of step S3 is:

[0027] S31, set the following material parameters, material thermal conductivity K, material thermal conductivity and temperature equation K(T), ambient temperature Ta, etc., and output the material parameters required for solving to a local file;

[0028] S32, applying boundary conditions and loads to the discretized model, specifying the geometric position of the boundary constraints, the load size, and the geometric position of the load application. And outputting the geometric position set information of the constrained load to a local file;

[0029] A further technical solution is that the specific implementation method of step S4 is:

[0030] S41. When solving with the modified arc length method, the initial load factor value needs to be set. The initial value of the load factor will affect the convergence of the solution. The iteration tolerance value is the convergence standard of the iteration step during the correction phase. The load factor sign S is positive by default. The maximum value of the load factor λ max =1, i.e. full load condition - reference load vector {Q ref} is equal to {Q e}External heat flux vector. The default initial load factor value is 0.1 and the iteration tolerance value is 10 -3According to the accuracy requirements of tetrahedron numerical integration, the number of integration points N of Hammer integration is preset to 4.

[0031] A further technical solution is that the specific implementation method of step S5 is:

[0032] S51. Before entering the improved arc length method solver, it is necessary to read the model information, including the unit node information, material information, boundary and initial load factor of the discretized model and store them. First, according to the given load and boundary conditions, the reference load heat flux vector {Q ref}. Secondly, read the initial load factor parameter λ from the local file initial , iteration tolerance t, maximum load factor λ max , the number of numerical integration points N.

[0033] The specific calculation steps are as follows:

[0034] First, in the jth iteration step under the i-th load increment step, the established load and boundary conditions are read, and the boundary conditions and load magnitudes of the constrained degrees of freedom are stored;

[0035] Secondly, the load value is calculated in the constrained area based on the finite element theory.

[0036]

[0037]

[0038] Where: is the overall conduction matrix of the jth iteration step under the i-th load increment; Ω is the unit domain; h is the convective heat transfer coefficient; N is the unit shape function; Γ is the load loading area after the model is discretized. G is the volume heat flux density, q is the surface heat flux density, T a is the ambient temperature.

[0039] S52, enter the main loop of the load increment step of the improved arc length method solver, first traverse the model unit, extract the node temperature of each unit of the discretized model, and calculate the thermal conductivity K at the Gaussian point of the current unit based on the material thermal conductivity and temperature relationship K(T). Secondly, the numerical integration method is used to solve the unit tangent conduction matrix under the jth iteration step of the i-th load increment step.

[0040] The specific calculation steps are as follows:

[0041] First, calculate the shape function N of the tetrahedral element of the current element in the local coordinate system (r, s, t) of the jth iteration step under the i-th load increment step;

[0042] Secondly, the partial derivative of the shape function for each coordinate component in the local coordinate system of the parent unit is obtained to obtain the partial derivative of the shape function of the parent unit for each coordinate axis in the local coordinate system.

[0043] Then the shape function in the local coordinate system (r, s, t) is connected with the shape function in the global coordinate system (x, y, z) through the shape function offset vector through Jacobian [J], and the calculation formula is:

[0044]

[0045]

[0046] N i～l is the element shape function, the subscripts i, j, k, l are the node numbers in the tetrahedral element; x i～l is the coordinate component of the x-axis of the unit node in the global coordinate system;

[0047] Finally, the gradient matrix [B] is obtained as:

[0048]

[0049] The thermal conductivity matrix [D] of isotropic materials is:

[0050]

[0051] Where: K x , K y , K z It is an intrinsic property of the material and its value is specified by the thermal conductivity of the material.

[0052] Unit finite element linear heat transfer matrix The calculation formula is:

[0053]

[0054] Where: is the unit thermal conductivity matrix under the jth iteration step of the i-th load increment step, and Ω is the unit domain.

[0055] Taking a single tetrahedron unit as an example, the finite element equation The specific form is as follows:

[0056]

[0057] For nonlinear problems, the specific steps for calculating the element tangent conduction matrix at the jth iteration step under the i-th load increment step are as follows:

[0058] When the model element type is a first-order tetrahedral element, the element Gauss point temperature calculation formula is:

[0059] T gauss =N i T i +N j T j +N k T k +N l T l

[0060] According to the relationship between thermal conductivity and temperature K(T) in the material parameter setting step, the temperature of the Gauss point is substituted into K(T) to obtain the unit thermal conductivity under the current configuration.

[0061] The linear term [D] of the thermal conductivity matrix of isotropic materials is:

[0062]

[0063] Where: K(T) x 、K(T) y , K(T)z is calculated based on the temperature value of the Gaussian integration point of the unit where the current material is located through the relationship between thermal conductivity and temperature K(T).

[0064] The nonlinear term [Dn] of the thermal conductivity matrix of isotropic materials is:

[0065]

[0066] The gradient matrix [B] is:

[0067]

[0068] At this point, the element tangent conduction matrix at the jth iteration step under the i-th load increment step during the nonlinear solution is The linear and nonlinear terms are calculated as follows:

[0069] The linear term is calculated as [K l ] e =∫ Ω [B] T [D][B]dΩ, the specific form is as follows:

[0070]

[0071] The calculation formula of nonlinear term is [K nl ] e =∫ Ω [B] T [D n ][B][T][N]dΩ, the specific form is as follows:

[0072]

[0073] The tangent conduction matrix of the cell is thus And store each matrix;

[0074] S53. The unit tangent conduction matrix in the above-mentioned global coordinate system saves the position index of its non-zero elements and the value of the corresponding position matrix element, and uses the COO storage method in the sparse matrix CSR method to store the sparse matrix; completes the global tangent conduction matrix of the solution domain at the jth iteration step under the i-th load increment step. Assembly.

[0075] S54, the overall tangent conduction matrix of the jth iteration step under the i-th load increment step is calculated by the large number method. External load heat flow rate vector {Q e} to process. Based on the boundary conditions, place a large value such as 10 on the degrees of freedom of the constrained point. 8 Even larger, which means that the constrained degrees of freedom are equal to the set boundary value dofvalue when solving. dofvalue can be obtained by reading the constrained degrees of freedom and boundary values ​​from the text file after the model is discretized.

[0076] The specific form of the large number method is as follows:

[0077]

[0078] S55. Configure the third-party library Eigen and set the linear equation solver to the conjugate gradient solver <bicgstab>, initialize the conjugate gradient solver and set the numerical precision of the solver. The overall tangent conduction matrix after the jth iteration step under the i-th load increment step is The external heat flow rate vector {Q e }Incoming, calling statement Solve the linear equations to obtain the tangent increment in the prediction phase The equation form is:

[0079]

[0080] S56. Determine whether the current load increment step is the first load increment step. If the current increment step is the first load increment step, there is no need to calculate the arc length. If the current increment step is not the first load increment step, it is necessary to determine the tangent increment of the previous load increment step prediction stage. Tangent increment to the current load increment The size of the second norm of is used to accelerate the algorithm according to the judgment result. The superscript i is the current load increment step, and the subscript 1 represents the prediction stage iteration step.

[0081] The arc length calculation formula is: ΔS i ={ΔT i-1 } T {ΔT i-1 }where superscript i represents the current load increment;

[0082] For heat conduction analysis with non-negative ambient temperature Ta, when the tangent increment norm in the incremental prediction stage meets the above requirements, the arc length calculation formula can be modified as follows:

[0083] (n is the system degree of freedom)

[0084] The specific calculation expression of arc length ΔS is as follows:

[0085] (n is the system degree of freedom)

[0086] Where: {ΔT} T The superscript T means transposition, ΔT is the temperature increment of the load increment step; The superscript i represents the number of load increment steps, and the subscript 1 represents the iteration step in the prediction phase;

[0087] S57. If the current load increment is the first load increment, the load factor The size is the initial load factor size λ initial , record the tangent increment vector of the load increment prediction stage. If the current load increment is not the first load increment, the load factor The specific calculation steps are as follows:

[0088]

[0089] Where: is the tangent increment obtained in the prediction stage; i is the load increment step; T is the transposition sign.

[0090] S58, Prediction stage load factor The calculation is completed, the tangent increment vector of the current load increment prediction phase is saved, the temperature vector and the load factor λ are updated, and the load increment prediction phase solution is ended.

[0091] The calculation formula is:

[0092]

[0093]

[0094] Where: is the tangent increment obtained in the prediction phase; i is the load increment step; T is the transposition sign

[0095] A further technical solution is that the specific implementation method of step S6 is:

[0096] S61. First, based on the updated load factor of the jth iteration step under the i-th load increment step, Update external load heat flow rate vector And the total load factor λ. The specific calculation steps are as follows;

[0097]

[0098]

[0099] Where: i is the load increment step; j is the iteration step; is the load factor increment step increment; is the load factor iteration step increment

[0100] S62, then traverse the unit and calculate the linear conduction matrix of the unit at the jth iteration step under the i-th load increment step And assembled into an overall linear conduction matrix The overall linear conduction matrix Multiplying by the node temperature vector {T} gives the element internal heat flow rate vector

[0101] The specific calculation steps are as follows:

[0102]

[0103] Where: The calculation of the overall linear conduction matrix is ​​the same as described in S52 and will not be repeated here. The subscript j-1 is because no calculation is performed in the prediction stage.

[0104] S63. Update the residual vector under the current configuration Determine the residual vector Is it less than the given residual? The residual vector here is The constrained degree of freedom position should be processed according to the boundary conditions, and the internal force at the constrained degree of freedom position is 0;

[0105] The calculation steps of the residual vector of the jth iteration step under the i-th load increment step are as follows:

[0106]

[0107] like The program will exit the calibration phase loop if the load factor λ is less than λ max If , the program will execute the next load increment calculation. Otherwise, the solution of the nonlinear thermal problem is completed. Then it enters the iterative correction stage.

[0108] S64, the iterative correction stage is mainly to eliminate the imbalance of heat flow rate inside and outside the system. First, the overall tangent conduction matrix of the jth iteration step under the i-th load increment step is calculated. Next is the temperature iteration correction vector and the residual iterative correction vector Solve, then update the load factor size according to the arc length constraint equation, and finally update the total displacement vector and the total load factor size.

[0109] The main calculation steps are as follows:

[0110] The constraint equation used in the algorithm is the cylindrical constraint equation, and its equation form is:

[0111] {ΔT i-1 } T {ΔT i-1 }=ΔS 2

[0112] Where: i represents the load increment step; ΔT i-1 represents the temperature increment of the previous load increment step;

[0113] The finite element equation for the jth iteration step under the i-th load increment step is as follows:

[0114]

[0115] The external heat flux load for the jth iteration step under the i-th load increment step is as follows:

[0116]

[0117] The temperature increment in the jth iteration step under the i-th load increment step is as follows:

[0118]

[0119] Therefore, the final form of the equation for the jth iteration step under the ith load increment step is as follows:

[0120]

[0121] Where: i represents the number of load increment steps; j represents the number of iteration steps; represents the temperature iteration correction vector; represents the residual iterative correction vector;

[0122] Therefore, the solution formula in the correction stage is:

[0123]

[0124]

[0125] The calculation steps of the unit and overall tangent conduction matrices are the same as S52 and will not be repeated here.

[0126] The external heat flow rate reference vector {Q ref }, the updated residual vector The processing is performed according to the boundary conditions, and the position where the degree of freedom is constrained is processed by the method of setting large numbers. The processing method is the same as S54, and will not be repeated here.

[0127] Configure the third-party library Eigen and set the linear equation solver to the conjugate gradient solver <bicgstab>, initialize the conjugate gradient solver and set the numerical precision of the solver. The temperature iteration correction vector can be obtained Similarly, the residual iterative correction vector can be obtained

[0128]

[0129] S65, Load factor solution in iterative correction stage

[0130] Temperature increment at the jth iteration under load increment i It can be expressed as follows:

[0131]

[0132] Where: i represents the number of load increment steps; j represents the number of iteration steps; represents the temperature iteration correction vector; represents the residual iterative correction vector;

[0133] Temperature increment Substituting this into the cylindrical constraint equation gives the load factor for the current iteration step. The quadratic equation for the variable is:

[0134]

[0135] Then the load factor of the correction phase of the jth iteration step under the i-th load increment step is The calculation formula is:

[0136]

[0137] Where: i represents the number of load increment steps; j represents the number of iteration steps; represents the tangent increment in the prediction phase; represents the temperature iteration correction vector; represents the residual iterative correction vector;

[0138] S66: Select the calculation method for the calibration phase according to the determination result in S56 and calculate the load factor Update the temperature vector and the total load factor λ. The iteration stops when the residual is less than the preset iteration tolerance.

[0139] A further technical solution is that the specific implementation method of step S7 is:

[0140] S71. When the total load factor value is greater than the preset maximum load factor value, the load increment step main loop is exited, and the total temperature vector under the current configuration is recorded. The solver completes the solution and stores the temperature vector of each node of the current model.

[0141] A further technical solution is that the specific implementation method of step S8 is:

[0142] S81. After the temperature field calculation of the entire solution domain is completed, the temperature field vector is obtained, and the temperature gradient matrix {g} is calculated and solved by traversing the units. The calculation formula is as follows:

[0143]

[0144] Where: T i , T j , T k , T l is the unit node temperature, the intermediate matrix is ​​the gradient matrix [B] matrix, N i 、N j 、N k 、N l is the unit shape function;

[0145] S82, calculation unit heat flow gradient {q} e After the temperature gradient calculation of the entire solution domain is completed, the heat flow of the calculation unit and the heat conduction satisfy Fourier's law. The formula is as follows:

[0146]

[0147] Where: q x ,q y ,q z is the heat flux in the x, y, and z directions; k is the thermal conductivity, dT / dx, dT / dy, and dT / dz are the temperature gradients;

[0148] A further technical solution is that the specific implementation method of step S9 is:

[0149] S91, according to the unit heat flux {q} e Calculate the average heat flux of the node. The average heat flux vector calculation formula is:

[0150]

[0151] Where: is the node average heat flow vector; {q} e is the heat flux size of the unit containing the node unit, and n is the number of units containing the node.

[0152] The beneficial effects of the present invention are as follows: the present invention discretizes the geometric model based on LiToMesh to establish a finite element model, and discretizes the geometric model with a tetrahedral unit type. The model information is then imported into a self-developed solver. By determining the size of the second norm of the tangent increment obtained in the prediction stage of adjacent incremental steps, the algorithm is accelerated according to the determination result, thereby achieving a reliable solution to the heat transfer nonlinear problem without using parallel technology to accelerate the algorithm solution process. Not only is the numerical accuracy of the calculation results improved, but also the speed of solving thermal nonlinear problems is improved. A solution is provided for heat transfer calculation scenarios that have high requirements for solution speed. BRIEF DESCRIPTION OF THE DRAWINGS

[0153] Figure 1 It is a flow chart of an improved arc length method for efficiently solving nonlinear thermal boundary problems of the present invention;

[0154] Figure 2 It is a schematic diagram of the geometric calculation model of the present invention;

[0155] Figure 3 is a schematic diagram of discrete unit types;

[0156] Figure 4 It is the finite element mesh diagram after the geometric model is discretized;

[0157] Figure 5 is a graph of the relationship between nonlinear heat flux load and temperature;

[0158] Figure 6 is the temperature field result obtained by solution;

[0159] Figure 7 It is the heat flow of nodes in all directions in three-dimensional space. DETAILED DESCRIPTION

[0160] The technical solution of the present invention will be described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0161] like Figure 1 As shown, the present invention provides an improved arc length method for efficiently solving nonlinear thermal boundary problems, comprising the following steps:

[0162] S1. Establishing a geometric model or geometric assembly model of a thermal analysis object to be determined;

[0163] The specific implementation steps are:

[0164] Establish a geometric model or geometric assembly model of any size at the macro scale and establish a continuous solution domain; the schematic diagram of the geometric model is as follows Figure 2 As shown;

[0165] S2, spatially discretize the geometric model to be analyzed, i.e., mesh division, to generate the mesh data required for calculation;

[0166] The specific implementation steps are:

[0167] S21, combining the geometric features of the large-scale geometric continuous model, setting the global mesh seed of the geometric model, the local seeds of some characteristic geometric edges, and setting local refinement parameters at the characteristic edges of some geometric features with complex features;

[0168] S22. Mathematical modeling of heat transfer problems will produce partial differential equations (PDEs). Numerical analysis, i.e., finite element method calculation, enables the partial differential equations to be approximately solved by numerical solutions of the original model domain.

[0169] K(T)=Q,Q∈R n

[0170] Where n is the total degree of freedom of the discretized system, K is the coefficient matrix of the linear equations of the system under the current boundary load problem, Q is the load vector, and the unknown quantity u can be solved using linear algebra methods;

[0171] S23, establishing a discretized geometric model;

[0172] By setting the meshing parameters and using the meshing tool to mesh the geometric model, the spatial discretization of large-scale geometric models of different shapes can be realized. The model can be discretized into a spatial continuous discrete body with tetrahedrons as units, and the data information of the discretized model can be written to a local file. The unit model is as follows: Figure 3 As shown, the discretized model is Figure 4 Shown

[0173] S3, setting parameters such as boundary conditions, loads, and material thermal conductivity for the divided grid model;

[0174] The specific steps include:

[0175] S31, set the following material parameters, material thermal conductivity K, material thermal conductivity and temperature equation K(T), ambient temperature Ta, etc., and output the material parameters required for solving to a local file;

[0176] S32, apply boundary conditions and loads to the discretized model, specify the geometric position of the boundary constraints, the load size, and the geometric position of the load application. Output the geometric position set information of the constrained load to a local file; the nonlinear function of the load size and the Gauss point temperature is as follows: Figure 5 As shown;

[0177] S4. Set the value of the initial load factor of the first incremental step of the modified arc length method λ initial , iteration tolerance value t, load factor symbol S, total load factor value λ, number of Hammer integration points N;

[0178] The specific steps include:

[0179] S41. When solving with the modified arc length method, the initial load factor value needs to be set. The initial value of the load factor will affect the convergence of the solution. The iteration tolerance value is the convergence standard of the iteration step during the correction phase. The load factor sign S is positive by default. The maximum value of the load factor λ max =1, i.e. full load condition - reference load vector {Q ref } is equal to {Q e }External heat flux vector. The default initial load factor value is 0.1 and the iteration tolerance value is 10 -3 According to the accuracy requirements of tetrahedron numerical integration, the number of integration points N of Hammer integration is preset to 4.

[0180] S5. To solve the problem under given load boundary conditions using the modified arc length method, it is necessary to calculate the current reference load heat flow rate vector {Q ref }. Enter the load increment prediction stage to solve. Calculate the tangent conduction matrix of the unit, the thermal conductivity of each unit under the current configuration, and the tangent prediction increment solution. Determine the size of the bi-norm of the tangent increment of the adjacent load increment step, and determine the acceleration method based on the judgment result;

[0181] The specific steps include:

[0182] S51. Before entering the improved arc length method solver, it is necessary to read the model information, including the unit node information, material information, boundary and initial load factor of the discretized model and store them. First, according to the given load and boundary conditions, the reference load heat flux vector {Q ref }. Secondly, read the initial load factor parameter λ from the local file initial , iteration tolerance t, maximum load factor λ max , the number of numerical integration points N.

[0183] The specific calculation steps are as follows:

[0184] First, in the jth iteration step under the i-th load increment step, the established load and boundary conditions are read, and the boundary conditions and load magnitudes of the constrained degrees of freedom are stored;

[0185] Secondly, the load value is calculated in the constrained area based on the finite element theory.

[0186]

[0187]

[0188] Where: is the overall conduction matrix of the jth iteration step under the i-th load increment; Ω is the unit domain; h is the convective heat transfer coefficient; N is the unit shape function; Γ is the load loading area after the model is discretized. G is the volume heat flux density, q is the surface heat flux density, T a is the ambient temperature.

[0189] S52, enter the main loop of the load increment step of the improved arc length method solver, first traverse the model unit, extract the node temperature of each unit of the discretized model, and calculate the thermal conductivity K at the Gaussian point of the current unit based on the material thermal conductivity and temperature relationship K(T). Secondly, the numerical integration method is used to solve the unit tangent conduction matrix under the jth iteration step of the i-th load increment step.

[0190] The specific calculation steps are as follows:

[0191] First, calculate the shape function N of the tetrahedral element of the current element in the local coordinate system (r, s, t) of the jth iteration step under the i-th load increment step;

[0192] Secondly, the partial derivative of the shape function for each coordinate component in the local coordinate system of the parent unit is obtained to obtain the partial derivative of the shape function of the parent unit for each coordinate axis in the local coordinate system.

[0193] Then the shape function in the local coordinate system (r, s, t) is connected with the shape function in the global coordinate system (x, y, z) through the shape function offset vector through Jacobian [J], and the calculation formula is:

[0194]

[0195]

[0196] N i～l is the element shape function, the subscripts i, j, k, l are the node numbers in the tetrahedral element; x i～l is the coordinate component of the x-axis of the unit node in the global coordinate system;

[0197] Finally, the gradient matrix [B] is obtained as:

[0198]

[0199] The thermal conductivity matrix [D] of isotropic materials is:

[0200]

[0201] Where: K x , K y , K z It is an intrinsic property of the material and its value is specified by the thermal conductivity of the material.

[0202] Unit finite element linear heat transfer matrix The calculation formula is:

[0203]

[0204] Where: is the unit thermal conductivity matrix under the jth iteration step of the i-th load increment step, and Ω is the unit domain.

[0205] Taking a single tetrahedron unit as an example, the finite element equation The specific form is as follows:

[0206]

[0207] For nonlinear problems, the specific steps for calculating the element tangent conduction matrix at the jth iteration step under the i-th load increment step are as follows:

[0208] When the model element type is a first-order tetrahedral element, the element Gauss point temperature calculation formula is:

[0209] T gauss =N i T i +N j T j +N k T k +N l T l

[0210] According to the relationship between thermal conductivity and temperature K(T) in the material parameter setting step, the temperature of the Gauss point is substituted into K(T) to obtain the unit thermal conductivity under the current configuration.

[0211] The linear term [D] of the thermal conductivity matrix of isotropic materials is:

[0212]

[0213] Where: K(T) x 、K(T) y 、K(T) z It is calculated based on the temperature value of the Gaussian integration point of the unit where the current material is located through the relationship between thermal conductivity and temperature K(T).

[0214] The nonlinear term [Dn] of the thermal conductivity matrix of isotropic materials is:

[0215]

[0216] The gradient matrix [B] is:

[0217]

[0218] At this point, the element tangent conduction matrix at the jth iteration step under the i-th load increment step during the nonlinear solution is The linear and nonlinear terms are calculated as follows:

[0219] The linear term is calculated as [K l ] e =∫ Ω [B] T [D][B]dΩ, the specific form is as follows:

[0220]

[0221] The calculation formula of nonlinear term is [K nl ] e =∫ Ω [B] T [D n ][B][T][N]dΩ, the specific form is as follows:

[0222]

[0223] The tangent conduction matrix of the cell is thus And store each matrix;

[0224] S53. The unit tangent conduction matrix in the above-mentioned global coordinate system saves the position index of its non-zero elements and the value of the corresponding position matrix element, and uses the COO storage method in the sparse matrix CSR method to store the sparse matrix; completes the global tangent conduction matrix of the solution domain at the jth iteration step under the i-th load increment step. Assembly.

[0225] S54, the overall tangent conduction matrix of the jth iteration step under the i-th load increment step is calculated by the large number method. External load heat flow rate vector {Q e } to process. Based on the boundary conditions, place a large value such as 10 on the degrees of freedom of the constrained point. 8 Even larger, which means that the constrained degrees of freedom are equal to the set boundary value dofvalue when solving. dofvalue can be obtained by reading the constrained degrees of freedom and boundary values ​​from the text file after the model is discretized.

[0226] The specific form of the large number method is as follows:

[0227]

[0228] S55. Configure the third-party library Eigen and set the linear equation solver to the conjugate gradient solver <bicgstab>, initialize the conjugate gradient solver and set the numerical precision of the solver. The overall tangent conduction matrix after the jth iteration step under the i-th load increment step is The external heat flow rate vector {Q e }Incoming, calling statement Solve the linear equations to obtain the tangent increment in the prediction phase The equation form is:

[0229]

[0230] S56. Determine whether the current load increment step is the first load increment step. If the current increment step is the first load increment step, there is no need to calculate the arc length. If the current increment step is not the first load increment step, it is necessary to determine the tangent increment of the previous load increment step prediction stage. Tangent increment to the current load increment The size of the second norm of is used to accelerate the algorithm according to the judgment result. The superscript i is the current load increment step, and the subscript 1 represents the prediction stage iteration step.

[0231] The arc length calculation formula is: ΔS i ={ΔT i-1 } T {ΔT i-1 }where superscript i represents the current load increment;

[0232] For heat conduction analysis with a non-negative ambient temperature Ta, when the second norm of the tangent increment in the incremental prediction stage meets the above requirements, the arc length calculation formula can be enlarged as follows:

[0233] (n is the system degree of freedom)

[0234] The specific calculation expression of arc length ΔS is as follows:

[0235] (n is the system degree of freedom)

[0236] Where: {ΔT} T The superscript T means transposition, ΔT is the temperature increment of the load increment step; The superscript i represents the number of load increment steps, and the subscript 1 represents the iteration step in the prediction phase;

[0237] S57. If the current load increment is the first load increment, the load factor The size is the initial load factor size λ initial , record the tangent increment vector of the load increment prediction stage. If the current load increment is not the first load increment, the load factor The specific calculation steps are as follows:

[0238]

[0239] Where: is the tangent increment obtained in the prediction stage; i is the load increment step; T is the transposition sign.

[0240] S58, Prediction stage load factor The calculation is completed, the tangent increment vector of the current load increment prediction phase is saved, the temperature vector and the load factor λ are updated, and the load increment prediction phase solution is ended.

[0241] The calculation formula is:

[0242]

[0243]

[0244] Where: is the tangent increment obtained in the prediction phase; i is the load increment step; T is the transposition sign

[0245] S6, enter the correction solution stage, first update the current external load vector {Q e }、Internal heat flow rate vector {Q a }, residual vector Secondly, determine whether the current load increment step has reached the convergence standard. After the final iteration phase is completed, determine whether the total load factor exceeds the preset maximum load factor value;

[0246] The specific steps include:

[0247] S61. First, based on the updated load factor of the jth iteration step under the i-th load increment step, Update external load heat flow rate vector And the total load factor λ. The specific calculation steps are as follows;

[0248]

[0249]

[0250] Where: i is the load increment step; j is the iteration step; is the load factor increment step increment; is the load factor iteration step increment

[0251] S62, then traverse the unit and calculate the linear conduction matrix of the unit at the jth iteration step under the i-th load increment step And assembled into an overall linear conduction matrix The overall linear conduction matrix Multiplying by the node temperature vector {T} gives the element internal heat flow rate vector

[0252] The specific calculation steps are as follows:

[0253]

[0254] Where: The calculation of the overall linear conduction matrix is ​​the same as described in S52 and will not be repeated here. The subscript j-1 is because no calculation is performed in the prediction stage.

[0255] S63. Update the residual vector under the current configuration Determine the residual vector Is it less than the given residual? The residual vector here is The constrained degree of freedom position should be processed according to the boundary conditions, and the internal force at the constrained degree of freedom position is 0;

[0256] The calculation steps of the residual vector of the jth iteration step under the i-th load increment step are as follows:

[0257]

[0258] like The program will exit the calibration phase loop if the load factor λ is less than λ max If , the program will execute the next load increment calculation. Otherwise, the solution of the nonlinear thermal problem is completed. Then it enters the iterative correction stage.

[0259] S64, the iterative correction stage is mainly to eliminate the imbalance of heat flow rate inside and outside the system. First, the overall tangent conduction matrix of the jth iteration step under the i-th load increment step is calculated. Next is the temperature iteration correction vector and the residual iterative correction vector Solve, then update the load factor size according to the arc length constraint equation, and finally update the total displacement vector and the total load factor size.

[0260] The main calculation steps are as follows:

[0261] The constraint equation used in the algorithm is the cylindrical constraint equation, and its equation form is:

[0262] {ΔT i-1 } T {ΔT i-1 }=ΔS 2

[0263] Where: i represents the load increment step; ΔT i-1 represents the temperature increment of the previous load increment step;

[0264] The finite element equation for the jth iteration step under the i-th load increment step is as follows:

[0265]

[0266] The external heat flux load for the jth iteration step under the i-th load increment step is as follows:

[0267]

[0268] The temperature increment in the jth iteration step under the i-th load increment step is as follows:

[0269]

[0270] Therefore, the final form of the equation for the jth iteration step under the ith load increment step is as follows:

[0271]

[0272] Where: i represents the number of load increment steps; j represents the number of iteration steps; represents the temperature iteration correction vector; represents the residual iterative correction vector;

[0273] Therefore, the solution formula in the correction stage is:

[0274]

[0275]

[0276] The calculation steps of the unit and overall tangent conduction matrices are the same as S52 and will not be repeated here.

[0277] The external heat flow rate reference vector {Q ref }, the updated residual vector The processing is performed according to the boundary conditions, and the position where the degree of freedom is constrained is processed by the method of setting large numbers. The processing method is the same as S54, and will not be repeated here.

[0278] Configure the third-party library Eigen and set the linear equation solver to the conjugate gradient solver <bicgstab>, initialize the conjugate gradient solver and set the numerical precision of the solver. The temperature iteration correction vector can be obtained Similarly, the residual iterative correction vector can be obtained

[0279] S65, Load factor solution in iterative correction stage

[0280] Temperature increment at the jth iteration under load increment i It can be expressed as follows:

[0281]

[0282] Where: i represents the number of load increment steps; j represents the number of iteration steps; represents the temperature iteration correction vector; represents the residual iterative correction vector;

[0283] Temperature increment Substituting this into the cylindrical constraint equation gives the load factor for the current iteration step. The quadratic equation for the variable is:

[0284]

[0285] Then the load factor of the correction phase of the jth iteration step under the i-th load increment step is The calculation formula is:

[0286]

[0287] Where: i represents the number of load increment steps; j represents the number of iteration steps; represents the tangent increment in the prediction phase; represents the temperature iteration correction vector; represents the residual iterative correction vector;

[0288] S66: Select the calculation method for the calibration phase according to the determination result in S56 and calculate the load factor Update the temperature vector and the total load factor λ. The iteration stops when the residual is less than the preset iteration tolerance.

[0289] S7, complete the thermal analysis solution and collect the temperature vector results of the final temperature field;

[0290] The specific steps include:

[0291] S71. When the total load factor value is greater than the preset maximum load factor value, the load increment step main loop is exited and the total temperature vector under the current configuration is recorded. The solver completes the solution and stores the temperature vector of each node of the current model. The temperature field result is as follows: Figure 6 As shown;

[0292] S8, calculate the unit conduction matrix, and calculate the temperature gradient and heat flow to achieve result post-processing;

[0293] The specific steps include:

[0294] S81. After the temperature field calculation of the entire solution domain is completed, the temperature field vector is obtained, and the temperature gradient matrix {g} is calculated and solved by traversing the units. The calculation formula is as follows:

[0295]

[0296] Where: T i , T j , T k , T l is the unit node temperature, the intermediate matrix is ​​the gradient matrix [B] matrix, N i 、N j 、N k 、N l is the unit shape function;

[0297] S82, calculation unit heat flow gradient {q} e After the temperature gradient calculation of the entire solution domain is completed, the heat flow of the calculation unit and the heat conduction satisfy Fourier's law. The formula is as follows:

[0298]

[0299] Where: q x ,q y ,q z is the heat flux in the x, y, and z directions; k is the thermal conductivity, dT / dx, dT / dy, and dT / dz are the temperature gradients;

[0300] S9. Calculate the average heat flux of the node to complete the solution of the nonlinear problem of thermal analysis under the given load boundary conditions; record the total time consumed by the solver.

[0301] The specific steps include:

[0302] S91, according to the unit heat flux {q} e Calculate the average heat flux of the node. The average heat flux vector calculation formula is:

[0303]

[0304] Where: is the node average heat flow vector; {q} e is the unit heat flux size of the node unit, and n is the number of units containing the node. The node heat flux calculation results are as follows: Figure 7 The time statistics of each method are shown in Table 1.

[0305] Table 1 Time consumption of the improved arc length method for solving nonlinear heat transfer problems

[0306]

[0307] The present invention discloses an improved arc length method for efficiently solving nonlinear thermal boundary problems. The method combines heat transfer theory, finite element method and nonlinear algorithm control strategy, improves the existing arc length method on the theoretical basis, and uses the excellent numerical calculation ability of computers as a carrier to solve thermal nonlinear problems, which not only improves the numerical accuracy of the calculation results, but also improves the speed of solving thermal nonlinear problems.

[0308] The above description does not limit the present invention in any form. Although the present invention has been disclosed by the above embodiments, it is not used to limit the present invention. Any technician familiar with the profession can make changes or modifications to equivalent embodiments of equivalent changes by using the technical contents disclosed above without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent changes and modifications made to the above embodiments according to the technical essence of the present invention without departing from the content of the technical solution of the present invention still fall within the scope of the technical solution of the present invention.< / bicgstab> < / bicgstab> < / bicgstab> < / bicgstab>

Claims

1. An improved arc length method for efficiently solving nonlinear thermal boundary problems, characterized in that: The following steps are involved: S1. Establishing a geometric model or geometric assembly model of a thermal analysis object to be determined; S2, spatially discretize the geometric model to be analyzed, i.e., mesh division, to generate the mesh data required for calculation; S3, setting boundary conditions, loads, and material thermal conductivity parameters for the divided grid model; S4. Set the value of the initial load factor of the first incremental step of the modified arc length method λ initial , iteration tolerance value t, load factor symbol S, total load factor value λ, number of Hammer integration points S5. To solve the problem under given load boundary conditions using the modified arc length method, it is necessary to calculate the current reference load heat flow rate vector {Q ref }; Enter the load increment prediction stage for solution; Calculate the tangent conduction matrix of the unit, the thermal conductivity of each unit under the current configuration, and solve the tangent prediction increment; determine the size of the second norm of the tangent increment of adjacent load increment steps, and determine the acceleration method based on the judgment result; S6, enter the correction solution stage, first update the current external load heat flow rate vector {Q e }、Internal heat flow rate vector {Q a }, residual vector Secondly, determine whether the current load increment step has reached the convergence standard; then enter the iterative correction solution process based on the judgment result; after the final iteration stage is completed, determine whether the total load factor exceeds the preset maximum load factor value; S7, complete the thermal analysis solution and collect the temperature vector results of the final temperature field; S8, calculate the unit conduction matrix, and calculate the temperature gradient and heat flow to achieve result post-processing; S9, calculate the average heat flux of the node to complete the solution of the nonlinear problem of thermal analysis under the given load boundary conditions; record the total time consumed by the solver; The specific implementation method of step S5 is: S51. Before entering the improved arc length method solver, it is necessary to read the model information, including the unit node information, material information, boundary and initial load factor of the discretized model and store them; first, according to the given load and boundary conditions, calculate the reference load heat flux vector {Q ref }; Secondly, read the initial load factor parameter λ from the local file initial , iteration tolerance t, maximum load factor λ max 、Number of Hammer points The specific calculation steps are as follows: First, in the jth iteration step under the i-th load increment step, the established load and boundary conditions are read, and the boundary conditions and load magnitudes of the constrained degrees of freedom are stored; Secondly, the load value is calculated on the constrained area based on the finite element theory; Where: is the global tangent conduction matrix of the jth iteration step under the i-th load increment; is the heat flux vector of the jth iteration step under the i-th load increment; Ω is the unit domain; h is the convective heat transfer coefficient; N is the shape function of the tetrahedral unit; Γ is the load loading area after the model is discretized; G is the volume heat flux density, q is the surface heat flux density, T a is the ambient temperature; S52, enter the main loop of the load increment step of the improved arc length method solver, first traverse the model unit, extract the node temperature of each unit of the discretized model, and calculate the thermal conductivity K at the Gaussian point of the current unit according to the material thermal conductivity and temperature relationship K(T); secondly, use the numerical integration method to solve the unit tangent conduction matrix under the jth iteration step of the i-th load increment step; The specific calculation steps are as follows: First, calculate the shape function N of the tetrahedral element of the current element in the local coordinate system (r, s, t) of the jth iteration step under the i-th load increment step; Secondly, the partial derivative of the shape function for each coordinate component in the local coordinate system of the parent unit is obtained to obtain the partial derivative of the shape function of the parent unit for each coordinate axis in the local coordinate system. Then the shape function in the local coordinate system (r, s, t) is connected with the shape function in the global coordinate system (x, y, z) through the shape function offset vector through Jacobian [J], and the calculation formula is: N i～l is the shape function of the tetrahedral element, and the subscripts i, j, k, and l are the node numbers in the tetrahedral element; x i～l is the coordinate component of the x-axis of the unit node in the global coordinate system; Finally, the gradient matrix [B] is obtained as: The thermal conductivity matrix [D] of isotropic materials is: Where: K x , K y , K z It is an intrinsic property of the material, and its value is specified by the thermal conductivity of the material; Unit finite element linear heat transfer matrix The calculation formula is: Where: is the unit finite element linear heat transfer matrix under the jth iteration step of the i-th load increment step, Ω is the unit domain; Finite element equations for a single tetrahedron element The specific form is as follows: Where: Q 1～4 is the heat source item corresponding to the unit node; For nonlinear problems, the specific steps for calculating the element tangent conduction matrix at the jth iteration step under the i-th load increment step are as follows: When the model element type is a first-order tetrahedral element, the element Gauss point temperature calculation formula is: T gauss =N i T i +N j T j +N k T k +N l T l Where, T i～l is the unit node temperature, N i～l is the shape function of the tetrahedral element, and the subscripts i, j, k, and l are the node numbers in the tetrahedral element; According to the relationship between thermal conductivity and temperature K(T) in the material parameter setting step, the temperature of the Gauss point is substituted into K(T) to obtain the unit thermal conductivity under the current configuration; The linear term [D] of the thermal conductivity matrix of isotropic materials is: Where: K(T) x 、K(T) y 、K(T) z It is calculated based on the temperature value of the Gaussian integral point of the unit where the current material is located through the relationship between thermal conductivity and temperature K(T); Nonlinear terms of thermal conductivity matrix of isotropic materials [D n ]for: Where: It is defined as the differential of thermal conductivity K(T) with respect to temperature T in the x, y, and z directions; The gradient matrix [B] is: At this point, the element tangent conduction matrix at the jth iteration step under the i-th load increment step during the nonlinear solution is The linear and nonlinear terms are calculated as follows: The linear term is calculated as [K l ] e =∫ Ω [B] T [D][B]dΩ, the specific form is as follows: The calculation formula of nonlinear term is [K nl ] e =∫ Ω [B] T [D n ][B][T][N]dΩ, the specific form is as follows: The tangent conduction matrix of the cell is thus And store each matrix; S53. The unit tangent conduction matrix in the above-mentioned global coordinate system saves the position index of its non-zero elements and the value of the corresponding position matrix element, and uses the COO storage method in the sparse matrix CSR method to store the sparse matrix; completes the global tangent conduction matrix of the solution domain at the jth iteration step under the i-th load increment step. Assembly; S54, the overall tangent conduction matrix of the jth iteration step under the i-th load increment step is calculated by the large number method. External load heat flow rate vector {Q e } to process; according to the boundary conditions, place a greater than 10 degree of freedom on the constrained point 8 The value of , which means that the constrained degrees of freedom are at the set boundary value dofvalue when solving; dofvalue can be obtained by reading the constrained degrees of freedom and boundary values ​​from the text file after the model is discretized; The specific form of the large number method is as follows: Where: K ij is the heat conduction from node i to j or the influence of structural stiffness, T1-Tn is the temperature vector of the node, Qe1-Qe n is the external load heat flow rate vector; S55. Configure the third-party library Eigen and set the linear equation solver to the conjugate gradient solver <bicgstab>, initialize the conjugate gradient solver, set the numerical precision of the solver; convert the overall tangent conduction matrix after the jth iteration step under the i-th load increment step External load heat flow rate vector {Q e }Incoming, calling statement compute(Q e ) to solve the linear equations and obtain the tangent increment in the prediction stage The equation form is:< / bicgstab> S56. Determine whether the current load increment step is the first load increment step. If the current increment step is the first load increment step, there is no need to calculate the arc length. If the current increment step is not the first load increment step, it is necessary to determine the tangent increment of the previous load increment step prediction stage. Tangent increment to the current load increment The size of the second norm of is used to accelerate the algorithm according to the judgment result; the superscript i is the current load increment step, and the subscript 1 represents the prediction stage iteration step; The arc length calculation formula is: ΔS i ={ΔT i-1 } T {ΔT i-1 }where superscript i represents the current load increment, ΔT i-1 represents the temperature increment of the previous load increment step; For heat conduction analysis with non-negative ambient temperature Ta, when the tangent increment norm of the incremental prediction stage meets the specified convergence criterion, the arc length calculation formula can be enlarged as follows: n2 is the system degree of freedom, The specific calculation expression of arc length ΔS is as follows: n is the system degree of freedom Where: {ΔT} T The superscript T means transposition, ΔT is the temperature increment of the load increment step; The superscript i represents the number of load increment steps, and the subscript 1 represents the iteration step in the prediction phase; S57. If the current load increment is the first load increment, the load factor The size is the initial load factor size λ initial , record the tangent increment vector of the load increment prediction stage; if the current load increment is not the first load increment load factor The specific calculation steps are as follows: Where: is the tangent increment obtained in the prediction stage; where i is the load increment step; T is the transposition symbol; S58, Prediction stage load factor After the calculation is completed, the tangent increment vector of the current load increment step prediction phase is saved, the temperature vector and the total load factor λ are updated; the load increment step prediction phase solution is ended; The calculation formula is: Where: is the tangent increment obtained in the prediction stage; i is the load increment step; T is the transposition sign; Δλ i-1 Represents the increment of the load factor completed in the load increment step i-1.

2. The improved arc length method for efficiently solving nonlinear thermal boundary problems according to claim 1, characterized in that: The specific implementation method of step S1 is: establishing a geometric model or a geometric assembly model of any scale at a macroscopic scale, and establishing a continuous solution domain.

3. The improved arc length method for efficiently solving nonlinear thermal boundary problems according to claim 1, characterized in that: The specific implementation method of step S2 is: S21, setting of discrete parameters; Combined with the geometric features of large-scale geometric continuous models, the global mesh seed of the geometric model and the local seeds of some characteristic geometric edges are set, and local refinement parameters are set at the characteristic edges with complex geometric features. S22. Mathematical modeling of heat transfer problems will produce partial differential equations (PDEs). Numerical analysis, i.e., finite element method calculation, enables the partial differential equations to be approximately solved by numerical solutions to the original model domain. Where: n1 is the total degree of freedom of the discretized system, K is the coefficient matrix of the linear equations of the system under the current boundary load problem, Q is the load vector, and the unknown quantity T can be solved using linear algebra methods; S23, establishing a discretized geometric model; By setting meshing parameters and using meshing tools to perform meshing of geometric models, the spatial discretization of large-scale geometric models of different shapes can be realized. The model can be discretized into spatial continuous discrete bodies with tetrahedrons as units, and the discretized data information of the model can be written to a local file.

4. The improved arc length method for efficiently solving nonlinear thermal boundary problems according to claim 1, characterized in that: The specific implementation method of step S3 is: S31, set the following material parameters, material thermal conductivity K1, material thermal conductivity and temperature related equation K1(T), ambient temperature Ta, and output the material parameters required for solving to a local file; S32, applying boundary conditions and loads to the discretized model, specifying the geometric position of the boundary constraints, the load magnitude, and the geometric position where the load is applied; And output the geometric position set information of the constrained load to a local file.

5. The improved arc length method for efficiently solving nonlinear thermal boundary problems according to claim 1, characterized in that: The specific implementation method of step S4 is: S41. When solving with the modified arc length method, the initial load factor value needs to be set. The initial value of the load factor will affect the convergence of the solution. The iteration tolerance value is the convergence standard of the iteration step during the correction phase. The load factor symbol S is positive by default. The maximum value of the load factor λ max =1, i.e., the full load reference load heat flow rate vector {Q ref } is equal to {Q e External load heat flux vector; preset initial load factor value is 0.1 and iteration tolerance value is 10 -3 ; According to the accuracy requirements of tetrahedron numerical integration, the number of integration points N of Hammer integration is preset to 4.

6. The improved arc length method for efficiently solving nonlinear thermal boundary problems according to claim 1, characterized in that: The specific implementation method of step S6 is: S61. First, based on the updated load factor of the jth iteration step under the i-th load increment step, Update external load heat flow rate vector And the total load factor λ; the specific calculation steps are as follows; Where: i is the load increment step; j is the iteration step; is the load factor increment step increment; is the load factor of the jth iteration step under the i-th load increment step; S62, then traverse the unit and calculate the linear conduction matrix of the unit at the jth iteration step under the i-th load increment step And assembled into an overall linear conduction matrix The overall linear conduction matrix Multiplying by the node temperature vector {T} gives the element internal heat flow rate vector The specific calculation steps are as follows: Where: The subscript j-1 is because no calculation is performed in the prediction stage; S63. Update the residual vector under the current configuration Determine the residual vector Is it less than the given residual? The residual vector here is The constrained degree of freedom position should be processed according to the boundary conditions, and the internal force at the constrained degree of freedom position is 0; The calculation steps of the residual vector of the jth iteration step under the i-th load increment step are as follows: like The program will exit the calibration phase loop if the total load factor λ is less than λ max When , the program executes the next load increment calculation; otherwise, the solution of the nonlinear thermal problem is completed; like Then it enters the iterative correction stage; S64, the iterative correction stage is to eliminate the imbalance of heat flow rate inside and outside the system; first, the overall tangent conduction matrix of the jth iteration step under the i-th load increment step is calculated Next is the temperature iteration correction vector and the residual iterative correction vector Solve, then update the load factor size according to the arc length constraint equation, and finally update the total displacement vector and the total load factor size; The calculation steps are as follows: The constraint equation used in the algorithm is the cylindrical constraint equation, and its equation form is: {ΔT i-1 } T {ΔT i-1 }=ΔS 2 Where: i represents the load increment step; ΔT i-1 represents the temperature increment of the previous load increment step; The finite element equation for the jth iteration step under the i-th load increment step is as follows: The external heat flux load for the jth iteration step under the i-th load increment step is as follows: The temperature increment in the jth iteration step under the i-th load increment step is as follows: Therefore, the final form of the equation for the jth iteration step under the ith load increment step is as follows: Where: i represents the number of load increment steps; j represents the number of iteration steps; represents the temperature iteration correction vector; represents the residual iterative correction vector; Represents the increment of the load factor in the jth iteration step under the ith load increment step With the reference load heat flow rate vector {Q ref }; Therefore, the solution formula in the correction stage is: The reference load heat flow rate vector {Q ref }, the updated residual vector According to the boundary conditions, the position where the degree of freedom is constrained is processed by the method of setting large numbers; the processing method adopts the method of step S54; Configure the third-party library Eigen and set the linear equation solver to the conjugate gradient solver <bicgstab>, initialize the conjugate gradient solver and set the numerical precision of the solver; pass in The temperature iteration correction vector can be obtained Similarly, the residual iterative correction vector can be obtained < / bicgstab> S65, Load factor solution in iterative correction stage Temperature increment at the jth iteration under load increment i It can be expressed as follows: Where: i represents the number of load increment steps; j represents the number of iteration steps; represents the temperature iteration correction vector; represents the residual iterative correction vector; Temperature increment Substituting this into the cylindrical constraint equation gives the load factor for the current iteration step. is a quadratic equation with one variable as the variable; therefore, during the calibration phase, the quadratic equation for the load factor is: Where: represents the temperature change increment from the j-1th iteration step to the jth iteration step under the i-th load increment step; Then the load factor of the correction phase of the jth iteration step under the i-th load increment step is The calculation formula is: Where: i represents the number of load increment steps; j represents the number of iteration steps; represents the tangent increment in the prediction phase; represents the transpose of the temperature iteration correction vector; represents the residual iterative correction vector; S66: Select the calculation method for the calibration phase according to the determination result in S56 and calculate the load factor Update the temperature vector and the total load factor λ; stop the iteration when the residual is less than the preset iteration tolerance.

7. The improved arc length method for efficiently solving nonlinear thermal boundary problems according to claim 1, characterized in that: The specific implementation method of step S7 is: S71. When the total load factor value is greater than the preset maximum load factor value, the load increment step main loop is exited, and the total temperature vector under the current configuration is recorded; the solver completes the solution and stores the temperature vector of each node of the current model.

8. The improved arc length method for efficiently solving nonlinear thermal boundary problems according to claim 1, characterized in that: The specific implementation method of step S8 is: S81. After the temperature field calculation of the entire solution domain is completed, the temperature field vector is obtained, and the temperature gradient matrix {g} is calculated and solved by traversing the units. The calculation formula is as follows: Where: T i , T j , T k , T l is the unit node temperature, the intermediate matrix is ​​the gradient matrix [B] matrix, N i 、N j 、N k 、N l is the unit shape function; S82, calculate the unit heat flux vector {q} e After the temperature gradient calculation of the entire solution domain is completed, the heat flow of the calculation unit and the heat conduction satisfy Fourier's law. The formula is as follows: Where: q x ,q y ,q z is the heat flux in the x, y, and z directions; k is the thermal conductivity, and dT / dx, dT / dy, and dT / dz are the temperature gradients.

9. The improved arc length method for efficiently solving nonlinear thermal boundary problems according to claim 1, characterized in that: The specific implementation method of step S9 is: S91, according to the unit heat flux vector {q} e Calculate the average heat flux of the node. The average heat flux vector calculation formula is: Where: is the node average heat flow vector; {q} e is the element heat flux vector, and n3 is the number of elements containing the node.

Images