Transient heat conduction and thermal stress simulation method for composite material based on two-dimensional arbitrary polygon hybrid finite element

Abstract

The invention relates to the technical field of transient heat conduction, in particular to a composite material transient heat conduction and thermal stress simulation method based on a two-dimensional arbitrary polygon hybrid finite element, which comprises the following steps: S1, determining an analysis model and material attributes, and establishing a transient heat conduction basic equation including boundary and initial conditions; s2, constructing and expanding a generalized variational functional, processing boundary conditions in combination with a Lagrange multiplier method, and discretizing a solution domain into a plurality of independent units; s3, adopting hypothetical temperature field, heat flux field and thermal stress field functions, introducing a discrete finite element format, and deducing a corresponding finite element equation from a functional; and S4, solving a transient temperature field through programming, introducing the transient temperature field into stress analysis, and finally obtaining temperature and stress distribution of the material at different time points. According to the method, temperature distribution and thermal stress of transient heat conduction are solved by hybridizing heat flux finite elements, and a new way is provided for transient heat conduction and thermal stress analysis of the composite material.

Description

Technical Field

[0001] The invention relates to the technical field of transient heat conduction, and in particular to a composite material transient heat conduction and thermal stress simulation method based on two-dimensional arbitrary polygon hybrid finite element. Background Art

[0002] Composite materials and composite structures used in mechanical manufacturing, aerospace, and medical fields will produce non-uniform transient temperature fields due to the instantaneous changes in service temperature, which will in turn produce complex thermal stresses and affect the service life of the structure. Therefore, in order to meet the design and application needs of new multifunctional composite materials in many fields, the study of their transient heat conduction and thermal stress is an extremely important link. In recent years, since thermodynamic experiments on composite materials generally require high costs, numerical simulation methods are currently one of the main research methods.

[0003] In the field of numerical simulation of transient heat conduction and thermal stress, Faruk et al. used the traditional finite element method to solve the thermal stress distribution of steel fiber reinforced curved composite disks. Bonito et al. developed an adaptive finite element method based on the finite element method to solve the diffusion problem in porous media. Compared with the finite element method, the meshless method can partially or completely get rid of the dependence on the computational grid. Ren et al. proposed a thermal-mechanical coupling algorithm based on the meshless method to simulate the ductile failure problem of materials. However, the meshless method is often only applicable to large gradient and large deformation problems. In addition, the multiscale analysis method is also a numerical method for thermodynamic analysis. It can capture the microscopic scale behavior of materials with less computational effort, but the solution of its homogenization equation still relies on the finite element method. For example, Yang et al. combined the multiscale analysis method with the finite element method and proposed a method for predicting the thermodynamic properties of materials. Therefore, it has to be admitted that the finite element method is one of the most effective numerical methods for solving thermodynamic problems. In recent years, in order to improve the disadvantages of the traditional finite element method such as sensitivity to mesh distortion and dependence on shape functions, the hybrid element method with advantages such as efficient solution has attracted much attention. However, hybrid elements seem to still be limited to steady-state conditions in the calculation of heat conduction and thermal stress. For example, Rao, Hao, and others used the stress hybrid element method to solve the thermal stress field with simple temperature boundary conditions, that is, the heat conduction of the material was not considered, and the entire model had the same temperature change value. Wang and others further solved the temperature in the domain through the node temperature, calculated the steady-state heat conduction and thermal stress using hybrid elements, and achieved good results. Therefore, proposing a hybrid element method that can calculate transient heat conduction and thermal stress will improve the efficiency of solving thermodynamic problems, and can more quickly study the thermal and mechanical properties of composite materials, which has a strong practical significance for the research and development and practical application of new composite materials.

[0004] The present invention is based on the generalized variational principle, and proposes a multivariate variational functional including two variables, heat flux and temperature, by introducing independent temperature fields and heat flux fields at the boundary and inside of the unit. In this method, the total heat flux of the system is divided into the sum of the steady-state heat flux that satisfies the heat balance equation in advance and the transient heat flux that balances the dynamic thermal energy of the system. The independent temperature field is introduced into the unit, and the temperature field obtained by the node temperature interpolation is assumed to be on the boundary of the unit, and the weak continuous form between the two parts of the temperature is introduced into the functional by the Lagrange multiplier method. Then, based on the functional, an arbitrary polygon hybrid heat flux finite element (PT-FEM) for solving two-dimensional transient heat conduction problems is derived. Next, the transient temperature change calculated in the previous step is introduced into the modified minimum complementary energy principle, and then the arbitrary polygon hybrid stress element for calculating transient thermal stress is derived. The calculation of transient temperature and transient thermal stress maintains the same grid form, which means that the advantage of the high efficiency of the coarse grid of the hybrid stress element can still be reflected in the calculation of transient thermal stress. Based on the above theory, the present invention realizes the solution of transient temperature field and transient thermal stress of two-dimensional composite materials through a self-compiled finite element program, and verifies the proposed method through some numerical examples. The verification results show that the method of the present invention has a flexible grid division method and satisfactory accuracy, which greatly reduces the number of units required by the finite element method in solving the transient temperature field and thermal stress field of composite materials, and provides a new idea for the efficient calculation of thermal stress of composite materials. Summary of the invention

[0005] The purpose of the present invention is to provide a method for simulating transient heat conduction and thermal stress of composite materials based on two-dimensional arbitrary polygon hybrid finite element. A hybrid heat flux finite element with high-order temperature field and heat flux field is developed to solve the temperature distribution of transient heat conduction, and then the calculated temperature field is transferred to the hybrid stress finite element in real time to solve the thermal stress. The feasibility and effectiveness of the method are verified through numerical examples, and the flexibility of arbitrary polygon units in meshing composite materials is demonstrated, providing a new way for transient heat conduction and thermal stress analysis of composite materials.

[0006] In order to achieve the above technical objectives and the above technical effects, the present invention is implemented through the following technical solutions:

[0007] A method for simulating transient heat conduction and thermal stress of composite materials based on two-dimensional arbitrary polygon hybrid finite element, comprising the following steps:

[0008] S1: Determine the analysis model and material properties, and establish the basic transient heat conduction equation including boundary and initial conditions;

[0009] S2: Construct and expand the generalized variational functional, combine the Lagrange multiplier method to deal with boundary conditions, and discretize the solution domain into multiple independent units;

[0010] S3: Using the assumed temperature field, heat flux field and thermal stress field functions, the discretized finite element format is introduced, and the corresponding finite element equations are derived from the functional;

[0011] S4: Solve the transient temperature field through programming and introduce it into the stress analysis, and finally obtain the temperature and stress distribution of the material at different time points.

[0012] Beneficial effects of the present invention:

[0013] This method overcomes the strict requirements of traditional finite element methods on mesh shape by using arbitrary polygons for meshing. This strategy can not only adapt to complex geometric boundaries, but also improve the flexibility of meshing. This method introduces independent temperature fields and high-order heat flux fields inside the unit, so even with a coarser grid, it can capture the local thermal conduction characteristics in detail. Through this flexible meshing, the introduction of redundant units is reduced, and the overall computational efficiency is significantly improved. It is a fast calculation method for predicting and testing the thermodynamic properties of new materials.

[0014] When dealing with transient heat conduction problems, the heat flux field and stress field in the units (triangles and quadrilaterals) used by traditional finite elements are determined by the temperature field and displacement field, respectively. Due to the limitation of shape functions, these units are low-order units, while the present invention provides a more accurate transient analysis framework by decomposing the total heat flux of the system into two parts: a steady-state heat flux that satisfies the heat balance equation and a dynamic equilibrium transient heat flux, while introducing a high-order stress field in the unit. Such a heat flux decomposition can ensure the accurate solution of the heat balance equation. Even when the temperature gradient inside the material changes dramatically, this method can obtain a more accurate transient temperature field evolution process with a smaller number of grids.

[0015] The present invention constructs a multivariate variational functional including a steady-state heat flux field, a transient heat flux field, and a thermal stress field, and uses the Lagrange multiplier method to introduce a weakly continuous form of boundary conditions, thereby effectively reducing the computational complexity and cost in the numerical solution process. This method can achieve efficient solutions to transient temperature fields and their induced thermal stresses within a limited number of units, avoiding the high computational costs of traditional methods due to multiple grid reconstructions. This advantage is particularly significant in the analysis of large-scale, complex structure composite materials.

[0016] The hybrid finite element method created by the present invention is not limited to solving heat conduction problems. Through appropriate deformation and extension, it can also be applied to other physical fields such as electromagnetic field analysis and fluid mechanics. Such wide applicability stems from its unique processing mechanism for the coupling of boundary conditions and field variables in multi-physics field problems, which makes this method have strong theoretical scalability. In addition, combined with multi-scale analysis methods, it can further enhance its precision in capturing the microscopic behavior of materials, providing important support for the design and performance optimization of new composite materials.

[0017] Of course, any product implementing the present invention does not necessarily need to achieve all of the advantages described above at the same time. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0019] Figure 1 Schematic diagram of the solution domain and boundary conditions;

[0020] Figure 2 Schematic diagram of the partitioning scheme of the integral domain;

[0021] Figure 3 is the calculation flow chart;

[0022] Figure 4 Schematic diagram of the calculation model and grid division scheme;

[0023] Figure 5 Computational cloud diagrams and numerical convergence verification for different grids, at t = 0.5s;

[0024] Figure 6 It is the temperature cloud map at different times and the comparison chart of the path sampling point values;

[0025] Figure 7 is the stress cloud diagram in the x direction at different times;

[0026] Figure 8 is the stress cloud diagram in the y direction at different times;

[0027] Fig. 9 It is the stress cloud diagram in the xy direction at different times;

[0028] Fig.10 This is a numerical comparison chart of the path sampling points at different times;

[0029] Fig.11 Schematic diagram of the calculation model and grid division scheme;

[0030] Fig.12 It is the temperature cloud map at different times and the comparison chart of the path sampling point values;

[0031] Fig.13 is the stress cloud diagram in the x direction at different times;

[0032] Fig.14 is the stress cloud diagram in the y direction at different times;

[0033] Fig.15 It is the stress cloud diagram in the xy direction at different times;

[0034] Fig.16 This is a numerical comparison chart of the path sampling points at different times;

[0035] Fig.17 Schematic diagram of the calculation model and grid division scheme;

[0036] Fig.18 Temperature cloud diagram at different times;

[0037] Fig.19 is the stress cloud diagram in the x direction at different times;

[0038] Fig. 20 is the stress cloud diagram in the y direction at different times;

[0039] Fig.21 It is the stress cloud diagram in the xy direction at different times. DETAILED DESCRIPTION

[0040] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. 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.

[0041] Example 1

[0042] The method for simulating transient heat conduction and thermal stress of composite materials based on two-dimensional arbitrary polygon hybrid finite element described in this embodiment includes the following steps:

[0043] S1: Determine the analysis model and material properties, and establish the basic transient heat conduction equation including boundary and initial conditions;

[0044] S2: Construct and expand the generalized variational functional, combine the Lagrange multiplier method to deal with boundary conditions, and discretize the solution domain into multiple independent units;

[0045] S3: Using the assumed temperature field, heat flux field and thermal stress field functions, the discretized finite element format is introduced, and the corresponding finite element equations are derived from the functional;

[0046] S4: Solve the transient temperature field through programming and introduce it into the stress analysis, and finally obtain the temperature and stress distribution of the material at different time points.

[0047] Example 2

[0048] Derivation of two-dimensional transient temperature-thermal stress equation

[0049] Hybrid Heat Flux Finite Element

[0050] The analysis domain is defined as a bounded closed region, with the inner part and the boundary of the domain denoted by Ω and The heat flux vector q(x,t) and the temperature T(x,t) are defined in Ω. The vector x represents a point in a given Cartesian coordinate system, and t∈I is the time variable, where I=(0,t end ).temperature and heat flux At the temperature boundary and heat flow boundary Applied on, where

[0051] When the heat source is not considered, the basic equation of the transient heat conduction problem can be written as:

[0052]

[0053] The boundary conditions are:

[0054]

[0055] The initial conditions at time t = 0 are:

[0056] T|0=T0,inΩ(4)

[0057] In the above equation, is the Hamiltonian operator, the superimposed points represent the time integral, ·0 represents the value corresponding to t = 0. k is the thermal conductivity of the material, ρ and c represent the density and specific heat respectively. n is the boundary The external normal vector on , T0 is the given initial temperature.

[0058] The method proposed in this embodiment includes both the temperature field and the heat flux field in the unit, and is a generalized variational functional. First, the multivariate variational functional of steady-state heat conduction can be generalized as:

[0059]

[0060] This formula can be extended to transient heat conduction by the following expression

[0061]

[0062] Among them, δΠ C Item and I D The terms are used to introduce transient thermal energy and relax the initial conditions, respectively.

[0063]

[0064] Discretization of the solution domain

[0065] like Figure 1 As shown in a, the solution domain Ω is discretized into multiple subdomains Ω with independent boundaries. e , each subdomain is a unit. In the heat conduction problem, if the subdomain boundary With temperature boundary or heat flow boundary intersect, then they are recorded as and In thermodynamic problems, if the subdomain boundaries With displacement boundary or traction boundary intersect, then they are recorded as and In addition, if Figure 1 As shown in b and 1c, since independent temperature fields are introduced at the unit boundary and inside the unit, it is necessary to establish formula (12) to relax the unit internal temperature T and the unit boundary temperature For transient heat conduction problems, the heat flux is given by the steady-state heat flux q s and transient heat flux q d In addition, the heat flux density q along the normal direction outside the unit boundary n Steady state and transient In the above discretization method, the following conditions exist: (i) the boundary temperature is Figure 1 d, and must be continuous between cells. Figure 1 The temperature and heat flux assumed in the cell shown in b are Figure 1 The unit boundary shown in c is weakly continuous. (iii) The heat flux density on the boundary is Figure 1 d) and is weakly continuous on the cell boundaries.

[0066] Based on the above assumptions, the heat flux of transient heat conduction consists of two parts:

[0067] q=q s +q d , inΩ e ×I (8)

[0068] Among them, the steady-state term satisfies the thermal balance equation in advance, while the transient term is an uncertain field and needs to ensure that it is in balance with the thermal energy in formula (2). That is, after substituting formula (8) into formula (2), the following equation is obtained:

[0069]

[0070] The heat flux in the normal direction is divided into a steady-state part and a transient part

[0071]

[0072] The steady-state part is defined by the steady-state heat flux and the external normal vector, while the transient term needs to be assumed separately.

[0073] Based on the above discretization method, the following conditions need to be introduced into the functional using the Lagrange multiplier method:

[0074]

[0075] and are the heat flux and temperature on the unit boundary, respectively, where It is obtained by interpolating the node temperature and is forced to be continuous on the unit boundary; considering the continuity of heat flow between adjacent units, the heat flux density The following conditions should be met

[0076]

[0077] in, is the boundary of adjacent cells, such as Figure 1 As shown in d.

[0078] The hybridization heat flux functional within each cell is modified as follows,

[0079]

[0080] Where K is the thermal conductivity matrix, ▽ is the Laplace operator,

[0081] So the discretized formula (6) becomes

[0082]

[0083] in, and As shown in formula (7). So far, we have obtained the modified multivariate variational functional containing three variables: heat flux in the domain, temperature in the domain, and boundary temperature. Next, we will derive the finite element format for programming calculation based on this functional.

[0084] Derivation of the finite element format

[0085] Based on the ideas of formulas (8) and (9), the heat flux field q and temperature field T in the unit are introduced independently. In addition, for the method proposed in this embodiment, it is also necessary to introduce the transient heat flux density in the normal direction of the unit boundary: For two-dimensional problems, the heat flux is a column vector containing two components. This embodiment introduces the assumed field by combining the known basis function P and the unknown coefficient β. The heat flux field and temperature field in the domain are shown below:

[0086]

[0087] Here, the subscripts h and θ are used to distinguish the heat flux and temperature, and s and d represent steady state and transient state, respectively.

[0088] The transient heat flux density of the normal component of the unit boundary,

[0089]

[0090] The boundary temperatures are interpolated from the node temperatures.

[0091]

[0092] Where Θ(t) is the node temperature at time t, is the temperature interpolation function of the ith edge, l i and a are the length of the ith edge and the distance from the ith node, respectively.

[0093] Substituting the above assumed fields into equations (7) and (15), we obtain

[0094]

[0095] in,

[0096]

[0097] Where, the subscripts s, d, and θ represent the parameters related to steady-state heat flux, transient heat flux, and temperature, respectively, and θ is the node temperature. P is a coordinate-related matrix whose form is known and independent of time. β is a coefficient matrix, which is a time-dependent quantity.

[0098] Partitioning the integration domain and solving equations

[0099] Formulas (23)-(26) involve a large number of domain integrals and boundary integrals. Choosing a reasonable integral domain partition scheme is an important step in the calculation. Figure 2 Taking the pentagonal unit shown in the figure as an example, the integral domain partitioning scheme used in this embodiment is introduced. Figure 2 Take node 1 in the example, insert midpoints a1 and b1 between adjacent nodes 2 and 5, and then insert points m1, m2, and m3 near the boundary. In this way, one node is divided into three integral domains, pentagonal elements have a total of 5×3 integral domains, and n-gonal elements are divided into n×3 integral domains. This division scheme increases the number of integral points at the boundary of the element, so that the calculation accuracy between adjacent elements is guaranteed to a certain extent, while also reducing the number of integral points required for the entire element as much as possible, thereby ensuring calculation efficiency.

[0100] Next, combining formula (16) and formula (20)-(22), we get

[0101]

[0102] in,

[0103]

[0104] From the variational condition of formula (16), we can get

[0105]

[0106] After sorting, we get

[0107]

[0108] in,

[0109]

[0110] Therefore, the solution form of the functional is transformed into solving the following differential equation,

[0111]

[0112] From the above derivation, it can be seen that the finite element established according to the functional of formula (15) is ultimately reduced to the matrix temperature method, that is, the temperature field and heat flux field are solved by solving the node displacement, which is a hybrid heat flux model.

[0113] According to the self-compiled finite element program, the above equations can be assembled and the temperature boundary and heat flow boundary can be applied. Then, it can be combined with the classical difference theory to achieve the solution of the time domain integral. Of course, a new time domain solution can also be developed based on formula (32). The time domain solution of this embodiment is as follows.

[0114]

[0115] Introduce the following equivalent mass matrix and equivalent load matrix,

[0116]

[0117] So we can get a relatively unified solution form,

[0118]

[0119] Among them, λ represents different difference methods, as shown in Table 1:

[0120] Table 1 Different difference formats

[0121]

[0122] Hybrid stress finite element method considering transient temperature field

[0123] The temperature T at time t solved above is now introduced into the hybrid stress functional established by Pian et al. for correction, and a hybrid thermal stress functional used in combination with the above-mentioned method is obtained.

[0124]

[0125] The functional contains two field variables, which act on Figure 1 The stress σ in the domain shown in e and the stress acting on Figure 1 f is the boundary displacement u. Among them, It acts on Figure 1 The given traction force on the traction force boundary shown in a, ε th It is the thermal strain due to temperature change and its form is as follows:

[0126] ε th =(T-T0)αA(38)

[0127] Where T and T0 are the temperature at time t and the initial temperature respectively, α is the thermal expansion coefficient, and A =

[110] T .

[0128] Next, we assume the stress field in the domain and the displacement field on the boundary respectively.

[0129] σ=P th β th , u=Ld(39)

[0130] Among them, L is the interpolation function on the boundary, so the functional can be transformed into:

[0131]

[0132] in,

[0133]

[0134] The integration method of the above equation is still the same Figure 2 As shown in the figure, because the two methods are unified in terms of grid division and integration, fast combined calculations can be performed.

[0135] Let the total energy be the stress coefficient β th and the first variation of the node displacement d is zero, we obtain:

[0136]

[0137] By combining the above equations, we can get a solution format similar to that of displacement finite element method.

[0138]

[0139] in, is the element stiffness matrix, It is called equivalent nodal load.

[0140]

[0141] At this point, the corresponding thermal stress can be quickly calculated by solving the temperature at each instant. It is worth noting that the current work has not established the coupling relationship between temperature and displacement. Therefore, this method is a transient-static combined hybrid finite element method.

[0142] Selection of field functions

[0143] Temperature field function and heat flux field function

[0144] For plane problems, the heat flux of the steady-state part is pre-satisfied with the heat balance equation. The steady-state heat flux field can be assumed by the heat flux function Φ.

[0145]

[0146] According to formula (9), the heat flux field in the transient part can be assumed by the temperature field and can be expressed as:

[0147]

[0148] The temperature field in equation (17) is approximated by a complete polynomial, and the first three orders are as follows:

[0149]

[0150] Taking the temperature field as a 2nd order complete polynomial as an example, the form of each field is as follows:

[0151]

[0152] Among them, in order to ensure the reversibility of the H matrix, it is necessary to select unrelated P s and P d Similarly, in order to ensure G d The reversibility of matrices, and P θ It can be chosen to be consistent. The number of terms in the basis function can be varied.

[0153] Stress Field Function

[0154] For plane problems, the thermal stress field is pre-satisfied with the equilibrium equation, which can be represented by the Airy stress function polynomial. The following is the stress function form containing 12 stress coefficients β.

[0155]

[0156] The above stress field is a field with 2nd order completeness, which is applicable to units with 7 or less nodes. In order to make the arbitrary polygonal unit proposed in this embodiment universal, the selection of the number of stress terms is discussed above. So far, the theoretical framework of the method proposed in this embodiment has been explained. Figure 3 The corresponding calculation flow chart is shown.

[0157] Numerical Examples of Arbitrary Polygonal Elements

[0158] Model with a single inclusion

[0159] This embodiment will use the commercial finite element software MSC.MARC and the method proposed in this embodiment (PT-FEM) to analyze the same typical example containing only a single inclusion. It is generally believed that the calculation results of MARC tend to get closer and closer to the true solution as the mesh is refined. Figure 4 a is a schematic diagram of the calculation model. Its size is 1×1mm, the radius of the circular inclusion is 0.15mm, the initial temperature of the transient heat conduction analysis is 25℃, and the value of the bottom of the model is 100℃ / mm 2 The other three sides are adiabatic boundaries, and the boundary conditions for thermal stress analysis are fixed on all sides. Figure 4As shown in b, compared with Marc's quadrilateral mesh, the polygonal mesh only needs to increase the number of edges on the inclusion boundary to well describe the shape of the inclusion. The material parameters used in the calculation are shown in Table 2. It is worth noting that in order to obtain a more obvious heat conduction trend in a shorter time step, the parameters used do not represent the actual situation and are only used for method verification. The time step of this embodiment is 0.05s, and the total step is 0.5s. In order to take into account both calculation accuracy and calculation efficiency, it is necessary to select appropriate temperature field functions and thermal stress field functions respectively. Therefore, while performing the MARC mesh convergence verification, PT-FEM with different numbers of items is also selected to solve the problem. The results are shown in Figure 5 shown.

[0160] Table 2 Material parameters

[0161]

[0162] from Figure 5 As can be seen from Figure 5, for MARC, as the number of units increases, the calculation results of the temperature field do not change significantly. Similarly, for PT-FEM, increasing the number of terms in the temperature field function also yields similar results. This is because in this embodiment, the form of the temperature field is relatively simple, and only a few units (MARC) or a lower-order temperature field function (PT-FEM) can be used to describe the temperature field of the problem. However, Figure 5 bd shows that it is far from enough to calculate the thermal stress field with only a few units, and the result gradually converges with the increase of the number of units. Therefore, for MARC, even if fewer units can complete the temperature calculation in heat conduction, in order to ensure the calculation accuracy of thermal stress, a relatively fine grid is still required, which reduces the efficiency of solving the problem. Figure 5 a shows that PT-FEM can not only calculate the temperature field with high accuracy, but also capture the stress concentration phenomenon of the thermal stress field under the premise of using only 68 arbitrary polygonal units. This is because the thermal stress field assumed by PT-FEM in each unit is a high-order field, so PT-FEM can still maintain satisfactory accuracy while greatly reducing the number of units to solve the problem. The comparative results of this example prove the computational efficiency and effectiveness of the method proposed in this embodiment.

[0163] In addition, in order to ensure both calculation efficiency and calculation accuracy, all calculation examples in this embodiment use the third-order (10β) temperature field function and the eighth-order (52α) thermal stress field function for calculation. Figure 6 is the comparison result of the temperature cloud map calculated by heat conduction and the numerical value of the sampling point, where Figure 6 jThe temperature change process of the material during transient heat conduction is reflected through the three points A, B and C. Figure 7-9The stress cloud diagrams of PT-FEM and MARC with 3168 elements at different times are shown. At the same time, the numerical comparison results of the calculation results of the two methods on different paths are shown in Figure 2. Fig.10 As shown, the numerical errors of all sampling points are within 7%, which once again proves the effectiveness of the method proposed in this embodiment.

[0164] Uniformly distributed multi-inclusion model

[0165] Generally speaking, particles in particle-reinforced composite materials often interact with each other, further changing the thermal and mechanical properties of the material. Therefore, whether the stress distribution between inclusions can be captured is also one of the indicators for measuring the calculation method. This embodiment uses the above model as a cell to array a 4×4mm model with 16 inclusions. The inclusions inside are evenly distributed. It is an idealized model and still uses the material parameters in Table 2. The calculation model and boundary conditions are as follows: Fig.11 As shown in a, Fig.11 b and 11c show the PT-FEM model of arbitrary polygonal elements and the MARC model of quadrilateral four-node elements, respectively. The total step length of this embodiment is 1 s and the time step is 0.1 s. Fig.12 The temperature cloud diagrams of the two methods at different times and the numerical comparison results on the x=0 path are shown. The calculated cloud diagram of thermal stress is shown in Figure 12-15 middle, Fig.16 It shows the numerical comparison results of stresses in different paths. It can be seen that in the early stage of the calculation, the temperature has not been transmitted to the entire model, but stress is still generated at locations far away from the heating boundary. This is because the model is fixed on all sides. Even if the temperature away from the heating boundary has not changed, the expansion of the bottom of the model causes the external displacement constraint to generate a reaction force, which in turn affects the stress distribution of the entire model. Therefore, the stress generated at a location far from the boundary is not caused by the temperature change at that point, but by the temperature change at the bottom. This example not only proves that HT-FEM can capture stress concentration near a single inclusion, but also describe the complex stress distribution of the entire composite material.

[0166] Randomly distributed multi-inclusion model

[0167] Homogenized particle reinforced composites are an ideal model. Some special methods may be able to design such materials, but the particle distribution in common composites is often random. Since the effectiveness of the calculation method has been verified, in order to simulate a more realistic situation, some random inclusion particles will be generated next to demonstrate the potential of the method proposed in this embodiment in practical applications. Fig.17As shown in a, the 2×2 mm model generated by a computer program has five inclusions of different shapes and sizes, and the position distribution is also random. This example still uses the material parameters in Table 2, the calculation time step is 0.05 s, and the total time is 0.5 s.

[0168] Since PT-FEM is an arbitrary polygonal unit, it has a very flexible meshing scheme, so it can be combined with flexible methods such as quadtree for pre-processing. For example, the developed polygon-quadtree hybrid method can be used to achieve fast meshing, which greatly reduces the initial modeling work. Fig.17 As shown in b, this embodiment uses a hybrid quadtree grid containing arbitrary polygons for solution. The calculated cloud map is shown in Figure 17-20 As shown in the figure, it can be seen that HT-FEM not only captures the stress concentration phenomenon around the inclusions very well, but also successfully describes the influence of the interaction between particles on the stress field distribution.

[0169] In summary, the present invention proposes an arbitrary polygon hybrid finite element with high computational efficiency and numerical accuracy for studying the two-dimensional transient heat conduction and thermal stress problems of composite materials. Since the method assumes high-order temperature field, heat flux field and stress field separately, the number of grids required for finite element calculation is greatly reduced under the premise of ensuring accuracy. It is a coarse grid high-performance finite element method. In addition, because the internal physical field no longer depends on the interpolation calculation of the shape function, the arbitrary polygon unit improves the flexibility and convenience of pre-processing work such as mesh division. The present invention successfully predicts the distribution of transient temperature field and thermal stress field of some particle reinforced composite examples, and compares them with the calculation results of traditional finite elements. The results show that the arbitrary polygon has a high adaptability to the inclusion particles in mesh division, and the calculation accuracy is also satisfactory. It not only captures the stress concentration phenomenon around the inclusion well, but also successfully describes the influence of the interaction between particles on the stress field distribution. This method provides an effective new way for the analysis of transient heat conduction and thermal stress of composite materials, and has an optimistic application prospect in the field of thermodynamic property prediction and development of composite materials. In addition, the present invention proposes a transient-static combined calculation method. The development of efficient hybrid elements for solving the bidirectional coupling of temperature and stress will be one of the difficult problems that need to be solved in the future.

[0170] The preferred embodiments of the present invention disclosed above are only used to help illustrate the present invention. The preferred embodiments do not describe all the details in detail, nor do they limit the invention to the specific implementation methods described. Obviously, many modifications and changes can be made according to the content of this specification. This specification selects and specifically describes these embodiments in order to better explain the principles and practical applications of the present invention, so that those skilled in the art can understand and use the present invention well. The present invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for simulating transient heat conduction and thermal stress of composite materials based on two-dimensional arbitrary polygon hybrid finite element, characterized in that: The following steps are involved: S1: Determine the analysis model and material properties, and establish the basic transient heat conduction equation including boundary and initial conditions; S2: Construct and expand the generalized variational functional, combine the Lagrange multiplier method to deal with boundary conditions, and discretize the solution domain into multiple independent units; S3: Using the assumed temperature field, heat flux field and thermal stress field functions, the discretized finite element format is introduced, and the corresponding finite element equations are derived from the functional; S4: Solve the transient temperature field through programming and introduce it into the stress analysis, and finally obtain the temperature and stress distribution of the material at different time points.

2. The composite material transient heat conduction and thermal stress simulation method based on two-dimensional arbitrary polygon hybrid finite element as claimed in claim 1, characterized in that: The step S1 specifically includes: The analysis domain is defined as a bounded closed region, with the inner part and the boundary of the domain denoted by Ω and = ... end ), temperature T(x,t) and heat flux At the temperature boundary and heat flow boundary Applied on, where The basic equation of transient heat conduction problem without considering the heat source is: q=-k▽T inΩ×I (1) The boundary conditions are: The initial conditions at time t = 0 are: T|0=T0,inΩ (4) ▽ is the Nabla operator, the superimposed points represent the time integral, ·|0 represents the value corresponding to t=0; k is the thermal conductivity of the material, ρ and c represent the density and specific heat respectively; n is the boundary The external normal vector on , T0 is the given initial temperature; The multivariate variational functional of steady-state heat conduction can be generalized as: This formula can be extended to transient heat conduction by the following expression Among them, δΠ C Item and I D The terms are used to introduce transient thermal energy and relax initial conditions respectively; 3. The composite material transient heat conduction and thermal stress simulation method based on two-dimensional arbitrary polygon hybrid finite element as claimed in claim 1, characterized in that: The step S2 comprises the following sub-steps: Discretize the solution domain Ω into multiple subdomains Ω with independent boundaries e , each subdomain is a unit; in heat conduction, if the subdomain boundary With temperature boundary or heat flow boundary intersect, then they are recorded as and In thermodynamics, if the subdomain boundaries With displacement boundary or traction boundary intersect, then they are recorded as and By establishing formula (12), the internal temperature T of the unit and the boundary temperature of the unit are relaxed. For transient heat conduction problems, the heat flux is given by the steady-state heat flux q s and transient heat flux q d In addition, the heat flux density q along the normal direction outside the unit boundary is n Steady state and transient Two parts; The heat flux of transient heat conduction consists of two parts q=q s +q d ,inΩ e ×I (8) Substituting equation (8) into equation (2), we get the following equation: The heat flux in the normal direction is divided into a steady-state part and a transient part The steady-state part is defined by the steady-state heat flux and the external normal vector, while the transient term needs to be set separately; Based on the above discretization method, the following conditions are introduced into the functional using the Lagrange multiplier method: and are the heat flux and temperature on the unit boundary, respectively, where It is obtained by interpolating the node temperature and is forced to be continuous on the unit boundary; considering the continuity of heat flow between adjacent units, the heat flux density Meet the following conditions in, is the adjacent cell boundary; The hybridization heat flux functional within each cell is modified as follows, Where K is the thermal conductivity matrix, ▽ is the Laplace operator, So the discretized formula (6) becomes in, and As shown in formula (7), the modified multivariate variational functional containing three variables, namely, the heat flux within the domain, the temperature within the domain, and the boundary temperature, is obtained. Based on the above functional, the finite element format for programming calculation is derived.

4. The composite material transient heat conduction and thermal stress simulation method based on two-dimensional arbitrary polygon hybrid finite element as claimed in claim 1, characterized in that: The step S3 specifically includes: Introduce the heat flux field q and temperature field T in the unit; introduce the transient heat flux density in the normal direction of the unit boundary The set field is introduced by combining the known basis function P and the unknown coefficient β; the heat flux field and temperature field in the domain are shown below: Wherein, the subscripts h and θ are symbols used to distinguish heat flux and temperature, s and d represent steady state and transient state, respectively; The transient heat flux density of the normal component of the unit boundary, The boundary temperatures are interpolated from the node temperatures. Where Θ(t) is the node temperature at time t, is the temperature interpolation function of the ith edge, l i and a are the length of the ith edge and the distance from the ith node respectively; Substituting the above set fields into equations (7) and (15), we obtain in, Among them, the subscripts s, d and θ represent the parameters related to steady-state heat flux, transient heat flux and temperature respectively, θ is the node temperature; P is the coordinate-related matrix, its form is known and has nothing to do with time; β is the coefficient matrix, which is a time-related quantity; Combining formula (16) and formula (20)-(22), we get in, From the variational condition of formula (16), we can get After sorting, in, Therefore, the solution form of the functional is transformed into solving the following differential equation, 5. The composite material transient heat conduction and thermal stress simulation method based on two-dimensional arbitrary polygon hybrid finite element as claimed in claim 1, characterized in that: The step S4 comprises the following sub-steps: Apply temperature boundaries and heat flow boundaries; combine with classical difference theory to solve the time domain integral; Introduce the following equivalent mass matrix and equivalent load matrix, We get a relatively unified solution form. Among them, λ represents different difference methods. In the forward difference format, λ=0; in the central difference format, λ=1 / 2; in the backward difference format, λ=1; The temperature T at time t solved above is introduced into the hybrid stress functional for correction to obtain the hybrid thermal stress functional; The functional contains two field variables, namely the stress σ in the domain and the boundary displacement u; where F is the given traction on the traction boundary, ε th It is the thermal strain due to temperature change and its form is as follows: ε th =(T-T0)αA(38) Where T and T0 are the temperature at time t and the initial temperature respectively, α is the thermal expansion coefficient, and A = [110] T ; Set the stress field in the domain and the displacement field on the boundary separately; σ=P th b th ,u=Ld(39) Among them, L is the interpolation function on the boundary, so the functional can be transformed into: in, Let the total energy be the stress coefficient β th The first variation of the node displacement d is zero, and we get: By combining the above equations, we can get a solution format similar to that of displacement finite element. in, is the element stiffness matrix, It is called equivalent nodal load; For plane problems, the heat flux of the steady-state part is pre-satisfied with the heat balance equation. The steady-state heat flux field is set by the heat flux function Φ. According to formula (9), the heat flux field of the transient part is set by the temperature field, which is expressed as: The temperature field in equation (17) is approximated by a complete polynomial, and the first three orders are as follows: Taking the temperature field as a 2nd order complete polynomial as an example, the form of each field is as follows:

Images