Three-dimensional VCFEM numerical simulation method considering temperature and plasticity for porous material

Abstract

The invention discloses a three-dimensional VCFEM numerical simulation method considering temperature and plasticity for a porous material. The plasticity and temperature-dependent thermal expansion at high temperature are considered. On the basis of the principle of minimum complementary energy, a corrected complementary energy functional in an increment form is deduced under the condition of considering thermal strain. And after the units are zoomed, each unit is subdivided into a plurality of Delaunay tetrahedrons, integration is carried out by using a Hammer numerical integration method, and a three-dimensional stress function considering an ellipsoid is constructed. Compared with the traditional finite element analysis, the three-dimensional VCFEM has the advantages that the calculation time required by each increment step in the nonlinear simulation is shorter, only less than 600 seconds are required for the same model and load condition, and the wall clock time is reduced by nearly 54 times. And the efficiency is greatly improved while the precision is kept.

Description

Technical Field

[0001] The present invention belongs to the technical field of material simulation and analysis, and in particular relates to a three-dimensional VCFEM numerical simulation method for porous materials taking temperature and plasticity into consideration. Background Art

[0002] Porous materials are widely found in nature and industrial applications, and can have a wide range of topological structures. From the nearly perfectly ordered grid structure found in honeycombs to the disordered three-dimensional networks observed in coral reefs. Due to their unique properties, porous materials have been widely used in various fields. For example, in the field of thermodynamics, porous insulation materials can effectively prevent heat transfer; in the field of aerospace, lightweight porous materials with high porosity can achieve aircraft structure weight reduction and performance improvement.

[0003] Porous materials are often exposed to complex environments during application, and in complex environments, porous materials will produce nonlinear strains, such as plastic strain and temperature-dependent thermal strain. When predicting the mechanical properties of porous materials in complex environments, the special pore microstructure will have a significant impact on the results. Therefore, it is generally necessary to predict the performance of porous materials. The main means at present is to analyze them through finite element analysis. The current finite element analysis mainly includes extended finite element method (XFEM), embedded finite element method (EFEM), hybrid finite element method (HFEM) and adaptive finite element method (AFEM). XFEM enhances the traditional finite element method by extending shape functions and using improved displacement interpolation functions to handle models with special microstructures (such as cracks, inclusions and voids). The embedded finite element method considers the strong discontinuity of displacement and stress fields in porous media, but its accuracy depends on the accuracy of embedded local details. The traditional displacement finite element method (FEM) usually requires a large number of units and multiple incremental steps to achieve a certain accuracy, resulting in high computational cost and low efficiency. Specifically, the traditional finite element method usually requires a large number of elements and incremental steps. For example, a 3D model with multiple pores may require more than 1,000 increments for nonlinear strain simulation, which requires finite element software (such as MSC MARC) to spend a lot of computing time, thus reducing the computational efficiency; there are still defects to a certain extent. Summary of the invention

[0004] In order to solve the above technical problems, the present invention designs a three-dimensional VCFEM numerical simulation method for porous materials considering temperature and plasticity. This method takes into account plasticity and temperature-dependent thermal expansion at high temperatures. Based on the principle of minimum complementary energy, a modified incremental complementary energy functional is derived for the case of considering thermal strain. After scaling the unit, each unit is subdivided into multiple Delaunay tetrahedrons, and the Hammer numerical integration method is used for integration to construct a three-dimensional stress function considering the ellipsoid. Compared with traditional finite element analysis, three-dimensional VCFEM requires less calculation time for each incremental step in nonlinear simulation. The same model and load conditions only require less than 600 seconds, and the wall clock time is reduced by nearly 54 times. While maintaining accuracy, efficiency is greatly improved.

[0005] In order to achieve the above technical effects, on the one hand, the present invention discloses a three-dimensional VCFEM numerical simulation method for porous materials considering temperature and plasticity, comprising the following steps:

[0006] S1: According to the existing initial complementary energy functional expression of the stress hybrid element model, according to the relationship between temperature, thermal strain and strain, and considering the boundary adjustment according to the plasticity increment theory, the complementary energy functional can be rewritten as follows:

[0007]

[0008] In the formula, Ω m is the area of ​​the matrix material, Γ e is the displacement boundary of the matrix material, Γ t is the force boundary of the whole unit, Γ i is the interface between the matrix material and the pores;

[0009] S2: Each increment is equal to the sum of the total increment of the previous step plus the increment of the current step, and the total energy of the system is obtained; at the same time, the motion relationship within the unit is obtained according to the Euler equation, which can be expressed as:

[0010]

[0011] Where, Δσ is the stress increment in the unit domain; Δε is the strain increment in the unit domain; Δu is the node displacement increment on the unit boundary; Δu1 is the node displacement increment on the matrix and hole interface close to the matrix side; Δu2 is the node displacement increment on the matrix and hole side close to the hole side; Δε is e —elastic strain increment within the unit; Δε p —Plastic strain increment within the element; Δε th —thermal strain increment within the element;

[0012] S3: The total energy of the system Π mcThe first-order variation of the displacement increment Δu can obtain the corresponding stress coordination relationship. At the same time, according to the stress function of the matrix material in the three-dimensional porous unit, the Maxwell stress function of the stress field and the interaction function; and the displacement increment of the unit boundary is expressed as the interpolation of the generalized displacement increment of each boundary segment node; the weak form of the system complementary energy after discretization can be expressed as the following formula:

[0013]

[0014] Δβ m Performing first-order variational calculus can obtain the stationary value condition of the modified complementary energy functional: And the weak form of the motion relationship within each unit can be obtained; at the same time, in the iterative solution process, the iterative expression of the i-th iteration step can be expressed as:

[0015]

[0016] [dβ mi ]=-[H m ] -1 [H m ] th +[H m ] -1 [G][dq] (5)

[0017] in,

[0018] S4: The first-order variation of the total capacity of the system with the node displacement increment obtains the weak form of the surface condition to be solved:

[0019]

[0020] Then substitute formula 5 into formula 6 to obtain:

[0021]

[0022] The solution of the matrix is:

[0023]

[0024] in, According to equations (8) and (9), we can get the solution relationship as follows:

[0025]

[0026] in:

[0027]

[0028] According to the above formula, the solution relationship can be obtained again, and the internal node displacement increment, the unit solution relationship, and the equivalent node force increment can be obtained; finally, it can be obtained:

[0029]

[0030] S5: Then the elastic-plastic constitutive equation of thermal strain is considered, and further derivation is carried out according to the basic laws of plastic mechanics. At the same time, this method is based on the small deformation thermal elastic-plastic constitutive relationship, considering the Mises yield condition and isotropic strengthening material; the specific basic laws of plastic mechanics include initial yield condition, plastic flow law, and plastic hardening law;

[0031] The specific steps include:

[0032] S5.1: Using the Mises yield condition, when the second stress invariant J2 = C (C is the material constant), the material will enter the yield state. The value of the material constant C can be obtained through uniaxial tension test. According to the expression of J2, the Mises yield condition can be rewritten as:

[0033]

[0034] S5.2: According to the Mises plastic potential concept expression, it is transformed into a flow law expression:

[0035]

[0036] in, The defined vector is along the normal direction of the subsequent yield surface F = 0 in space;

[0037] When the yield condition is the Mises yield condition, substituting the Mises yield function equation (the expression of J2) into equation (17) yields:

[0038]

[0039] Then rewrite the tensor expression (17) into a general expression, transform the equivalent stress, and then substitute (18) into the general equation to obtain:

[0040]

[0041] S5.3: The plastic hardening law considers the process as follows: The hardening law stipulates that the material enters the loading surface determined by the subsequent yield function after plastic deformation, so the general loading function is used:

[0042]

[0043] Among them, C is the hardening parameter. For the 2024-T6 aluminum alloy material used in this paper, since there is no reverse loading and unloading process, it will not reflect the Bauschinger effect, so the linear isotropic hardening law is adopted; then the Mises yield condition is adopted, and the subsequent yield function of isotropic hardening can be expressed as:

[0044]

[0045] Among them, σ ij is the elastic-plastic stress, which is the equivalent plastic strain The function of It can be obtained from the σ-ε curve of the uniaxial tensile test of the material; definition: Among them, H is the plastic modulus, also known as the hardening modulus, then H is related to the elastic modulus E and the tangent modulus The relationship is:

[0046]

[0047] S6: According to the Parandtl-Reuss theory, the strain increment (elastic strain increment and plastic strain increment) is considered and further combined with the Mises yield law and hardening law to derive the expression:

[0048]

[0049] At the same time, the unit flexibility matrix is ​​established for the three-dimensional problem, and the final expression is:

[0050]

[0051] S7: In order to determine whether plastic loading or elastic unloading should be performed in 3D VCFEM, and thus decide whether to use the plastic constitutive relation or the elastic constitutive relation, the following formula is established for judgment:

[0052]

[0053] in It is a material parameter determined by the yield function. When condition (1) in the equation is met, plastic loading is continued. When condition (2) is met, elastic loading is selected. When condition (3) is met, for ideal elastic-plastic materials, new plastic flow is generated. For hardened materials, no new plastic flow occurs.

[0054] Furthermore, in S1, the initial complementary energy functional expression of the existing stress hybridization metamodel is:

[0055]

[0056] Among them, S is the elastic-plastic flexibility matrix, u is the given boundary condition, T is the surface force of the current boundary, Γt e is the displacement boundary of the unit;

[0057] The expression of thermal stress is:

[0058] Among them, α(T) is the temperature-dependent thermal expansion coefficient, and α(T) can be defined as a function related to the temperature field variables.

[0059] Furthermore, in S2, each increment is equal to the sum of the total increment of the previous step plus the increment of the current step:

[0060]

[0061] The total energy of the system is expressed as:

[0062]

[0063] The beneficial effects of the present invention are:

[0064] The present invention discloses a modified complementary energy functional representation method based on a three-dimensional Voronoi unit finite element method (3D VCFEM), which takes into account plasticity and temperature-dependent thermal expansion at high temperatures. Based on the minimum complementary energy principle, a modified incremental complementary energy functional is derived for the case of thermal strain. After scaling the unit, each unit is subdivided into multiple Delaunay tetrahedrons, and the Hammer numerical integration method is used for integration to construct a three-dimensional stress function considering the ellipsoid. Compared with traditional finite element analysis, three-dimensional VCFEM requires less calculation time for each incremental step in nonlinear simulation. The same model and load condition only requires less than 600 seconds, and the wall clock time is reduced by nearly 54 times. While maintaining accuracy, efficiency is greatly improved.

[0065] In addition, through further optimization of this method, the three-dimensional VCFEM not only achieves the same accuracy as the traditional finite element method (FEM) while using fewer nodes and elements, but also highlights that the correction analysis performed by this method has a higher computational efficiency than the existing finite element analysis and has outstanding effects; it can simulate the nonlinear behavior of porous materials in complex environments, which is a significant improvement. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0067] Figure 1 is the average temperature coefficient diagram of the present invention;

[0068] Figure 2 is a tetrahedron diagram in three-dimensional space in the present invention, in which a is a tetrahedron in three-dimensional space, b is a triangular face of the tetrahedron, and c is a schematic diagram of a point P in the triangle;

[0069] Figure 3 It is a comparison diagram of the results of MSC MARC for different POL values ​​in the present invention;

[0070] Figure 4 Schematic diagram of the single hole displacement boundary condition of the present invention; in the figure, a is the YZ plane of the single hole displacement boundary condition, b is the YX plane of the single hole displacement boundary condition, c is the single hole working condition diagram, d is the YZ plane of the regular distribution displacement boundary condition, e is the YX plane of the regular distribution displacement boundary condition, and f is the regular distribution working condition diagram;

[0071] Figure 5 is the stress distribution diagram of the present invention, where a is σ y Stress distribution diagram (25th increment in 3D VC FEMC), b is σ y Stress distribution diagram (25th increment in MSC MARC);

[0072] Figure 6 This is a slice diagram of the stress result of the temperature-dependent thermal expansion of a single hole in the present invention, where a is σ y Stress distribution diagram (10th increment in 3D VC FEMC), b is σ y Stress distribution diagram (10th increment in MSC MARC), c is σ y Path along the Y axis; d is σ y Along the Y axis;

[0073] Figure 7 This is a slice diagram of the stress result when the temperature depends on thermal expansion of the present invention, where a is σ y Stress distribution diagram (10th increment in 3D VC FEMC), b is σ y Stress distribution diagram (10th increment in MSC MARC) c is σ y Path along the Y axis; d is σ y Along the Y axis

[0074] Figure 8 This is a slice diagram of the stress result during plastic strain of the present invention, where a is σ y Stress distribution diagram (25th increment in 3D VC FEMC), b is σ y Stress distribution diagram (25th increment in MSC MARC), c is σ y Path along the Y axis; d is σ y Along the Y axis;

[0075] Fig. 9This is an analysis diagram of randomly distributed holes of the present invention; a in the figure is a VCFEM grid; b is σ y Stress distribution diagram (25th increment), c is a is VCFEM mesh (10th increment), d is σ y Stress distribution diagram (10th increment). DETAILED DESCRIPTION

[0076] Example 1

[0077] The present invention discloses a three-dimensional VCFEM numerical simulation method for porous materials considering temperature and plasticity, comprising the following steps:

[0078] According to the hypothetical stress hybrid element model proposed by Pian's [10-12], the initial complementary energy functional expression (1) is given:

[0079]

[0080] Where S is the elastic-plastic flexibility matrix, is the given boundary condition, T is the surface force of the current boundary, Γt e is the displacement boundary of the unit. When the thermal deformation (expansion or contraction) caused by temperature change is constrained, stress will be generated; this type of stress is usually called thermal stress / temperature stress. Assuming that the temperature change of the porous material is ΔT, it can be known that its thermal strain along a certain direction is:

[0081]

[0082] Among them, α(T) is the temperature-dependent thermal expansion coefficient, and α(T) can be defined as a function related to the temperature field variable. The thermal strain caused by temperature change can be defined as ε th , then the total strain of the porous material can be defined as:

[0083] ε=σS+ε th (3)

[0084] Assuming that the traction between the elements is balanced, the following equation is satisfied at a given traction boundary:

[0085]

[0086] At a given cell boundary On the surface, the surface force satisfies the equilibrium condition of equation (5):

[0087] n + σ=n - ·σ (5)

[0088] According to the plastic increment theory, the complementary energy functional can be rewritten as:

[0089]

[0090] Among them, the area of ​​the matrix material Ω m , the displacement boundary of the matrix material Γ e , the force boundary Γ of the whole unit t , the interface between the matrix material and the pores Γ i Each increment is equal to the sum of the total increment of the previous step plus the increment of the current step:

[0091]

[0092] The total energy of the system is:

[0093]

[0094] The unit energy The first-order variation of the stress increment Δσ can be used to obtain the motion relationship within the unit (Euler equation):

[0095]

[0096] Among them, the unknown field variables are: Δσ—stress increment in the unit domain, Δε—strain increment in the unit domain, Δu—node displacement increment on the unit boundary, Δu1—node displacement increment on the matrix and hole interface close to the matrix side, Δu2—node displacement increment on the matrix and hole side close to the hole side. Δε e —Elastic strain increment within the unit, Δε p —Plastic strain increment within the element, Δε th —The thermal strain increment within the unit. mc For the first-order variation of the displacement increment Δu, the corresponding stress coordination relationship can be obtained:

[0097]

[0098] n i (σ i +Δσ i )=0 (15)

[0099] The stress function of the matrix material in the three-dimensional porous element should be selected as:

[0100] [Δσ m ]=[P m ][Δβ m ] (16)

[0101] The construction of the stress field can be composed of two parts: the Maxwell stress function and the interaction function, as shown in equation (17):

[0102]

[0103] The displacement increment of the element boundary can be expressed as the interpolation of the generalized displacement increments of each boundary segment node:

[0104] [Δu1]=[L1][Δq e ] (18) [Δu2]=[L2][Δq i ] (19) The weak form of the system complementary energy after discretization can be obtained:

[0105]

[0106] Equation (20) can be rewritten as:

[0107]

[0108] Δβ m Performing first-order variation can obtain the stationary value condition of the modified complementary energy functional:

[0109]

[0110] The weak form of the motion relationship within each unit can be obtained:

[0111]

[0112] In the iterative solution process, the stress parameter increment of the i-th iteration step is expressed as:

[0113] [Δβ mi ]=[Δβ mi-1 ]+[dβ mi ] (twenty four)

[0114] The displacement increment of the i-th step can be expressed as:

[0115]

[0116] The iteration expression is:

[0117]

[0118] [dβ mi ]=-[H m ]- 1 [H m ] th +[H m ]- 1 [G][dq] (27) where:

[0119]

[0120] The first-order variation of the total capacity of the system with the node displacement increment obtains the weak form of the surface condition to be solved:

[0121]

[0122] (-[β m ] T +[Δβ mi ])[G m 2 ]=0 (31)Wherein:

[0123]

[0124]

[0125] The above equation can be rewritten as:

[0126]

[0127] Substituting equation (27) into equation (34), we have:

[0128]

[0129] At this time, the solution of the matrix can be written as:

[0130]

[0131] Equation (36) can be rewritten as:

[0132]

[0133] in:

[0134]

[0135] The solution relationship can be rewritten as:

[0136]

[0137] in:

[0138]

[0139] The right side of equation (39) can be rewritten as:

[0140]

[0141] Equation (45) and Equation (46) can be rewritten as:

[0142]

[0143] The internal node displacement increment can be expressed as:

[0144] dq in i =K 22 -1 (dF in -K 12 T dq ei ) (49)

[0145] Substituting equation (48) into equation (49), we obtain:

[0146]

[0147] The solution relation of the unit can be re-expressed as:

[0148]

[0149] in:

[0150]

[0151] The equivalent nodal force increment can be obtained:

[0152]

[0153] Finally we get:

[0154]

[0155] Then, the Mises yield condition is used. When the second stress invariant J2=C (C is the material constant), the material will enter the yield state. The value of the material constant C can be obtained through uniaxial tension test. According to the expression of J2, the Mises yield condition can be rewritten as:

[0156]

[0157] According to the Mises plastic potential concept expression, it is transformed into the flow law expression:

[0158]

[0159] in, The defined vector is along the normal direction of the subsequent yield surface F = 0 in space;

[0160] When the yield condition is the Mises yield condition, substituting the Mises yield function equation (the expression of J2) into equation (56) yields:

[0161]

[0162] Then rewrite the tensor expression (17) into a general expression, transform the equivalent stress, and then substitute (18) into the general equation to obtain:

[0163]

[0164] The plastic hardening law considers the process of insertion. The hardening law stipulates that the material enters the loading surface determined by the subsequent yield function after plastic deformation, so the general loading function is used:

[0165]

[0166] Among them, C is the hardening parameter. For the 2024-T6 aluminum alloy material used in this paper, since there is no reverse loading and unloading process, it will not reflect the Bauschinger effect, so the linear isotropic hardening law is adopted; then the Mises yield condition is adopted, and the subsequent yield function of isotropic hardening can be expressed as:

[0167]

[0168] Among them, σ ij is the elastic-plastic stress, which is the equivalent plastic strain The function of It can be obtained from the σ-ε curve of the uniaxial tensile test of the material; definition: Among them, H is the plastic modulus, also known as the hardening modulus, then H is related to the elastic modulus E and the tangent modulus The relationship is:

[0169]

[0170] S6: According to the Parandtl-Reuss theory, the strain increment (elastic strain increment and plastic strain increment) is considered and further combined with the Mises yield law and hardening law to derive the expression:

[0171]

[0172] At the same time, the unit flexibility matrix is ​​established for the three-dimensional problem, and the final expression is:

[0173]

[0174] It is determined whether to use the plastic constitutive relation or the elastic constitutive relation, so the following formula is used to make a judgment:

[0175]

[0176] in It is a material parameter determined by the yield function. When condition (1) in the equation is met, plastic loading is continued. When condition (2) is met, elastic loading is selected. When condition (3) is met, for ideal elastic-plastic materials, new plastic flow is generated. For hardened materials, no new plastic flow occurs.

[0177] The calculation of the present invention uses 2024-T6 aluminum alloy to simulate its average thermal expansion coefficient in a temperature environment of 20°C-1000°C using Jamtpro, such as Figure 1 As shown by Figure 1 The fitting equation can be obtained as:

[0178]

[0179] Using the 3D VCFEM method, the element is first divided into multiple Delaunay tetrahedrons, each face of which is a triangular surface. In order to solve the final displacement, it is necessary to construct an interpolation function on the surface of these tetrahedrons. In the rectangular coordinate system, the coordinates of vertices i, j, and k are given, and the displacement of the vertex is as follows Figure 2 (b) Then, the boundary displacement of the element surface can be obtained by interpolating the node displacement:

[0180]

[0181] in Δ represents the area of ​​a triangle, which can be calculated based on the coordinates of the triangle vertices in space, and Δ m represents the area of ​​the smaller triangle consisting of two vertices and an internal point p(x, y, z), so the above formula can be rewritten as: u = Ld (67)

[0182] Equation () can also be expressed as the total amount of equation () and equation ()

[0183] in:

[0184]

[0185] Where L m (m=ijk) represents the displacement interpolation function of the triangle. Connecting any point p(x,y,z) inside the triangle to the vertex of the triangle can subdivide the large triangle into three small triangles, and the area of ​​each small triangle is Δ m (m=i,j,k), the total area of ​​the triangle can be calculated based on the coordinates of the triangle vertices Δ=Δ i +Δ j +Δ k And the area of ​​each small triangle:

[0186]

[0187] According to the equation ()L m , Δ and Δ m The relationship between can be expressed as equation ():

[0188]

[0189] In addition, according to the above conditions, the coordinates of any point p(x,y,z) on each triangular face can be determined:

[0190]

[0191] When solving the G matrix in the formula, it is also necessary to calculate the cosine of the normal vector of each face of the tetrahedron. In this way, the normal vector matrix can be constructed, where the normal vector of each triangular face is given by the formula:

[0192]

[0193] Therefore, the triangular cosine matrix is ​​expressed as:

[0194]

[0195] in:

[0196]

[0197] In three-dimensional space, the stress function tensor can be expressed as follows:

[0198]

[0199] The stress form that satisfies both the equilibrium equation and the compatibility equation is given by Equation (101) and Equation (102):

[0200]

[0201]

[0202] Equation (78) represents Maxwell's formula for solving the stress components, and equation (79) corresponds to Morella's formula for solving the stress components.

[0203] In three-dimensional problems, it is very important to choose an appropriate stress function to accurately represent the stress field inside the element. The stress function of the matrix is ​​expressed as , and the stress function of the interaction force in the matrix is ​​expressed as The construction of the polynomial stress function and the interaction force function must follow certain criteria: the columns in the H matrix must be linearly independent to ensure the invertibility of the H matrix. The choice of the polynomial stress function should take into account the shape of the pores, ensure that the interaction stress terms approach zero far away from the boundaries where different phases contact, and maintain the continuity of the surface forces on these boundaries.

[0204]

[0205] In this study, the interface interaction stress function derived by Ghosh was further simplified and expressed as:

[0206]

[0207] For each Voronoi element containing a void, the stress function of the matrix can be expressed as:

[0208]

[0209]

[0210] This method uses the Maxwell stress function commonly used in hybrid stress elements as the stress field introduced in the hybrid stress element. Take the third-order perfect polynomial as an example:

[0211]

[0212] The present invention adopts Maxwell stress function. Therefore, the stress components can be obtained as follows:

[0213]

[0214] From the construction method of the element stiffness matrix in the stress hybrid element, to ensure that the stiffness matrix is ​​non-singular, the following equation must be satisfied

[0215] n β ≥n n ·DD F (86)

[0216] Among them, n β is the number of stress parameters, D is the number of nodes in the element, and D F is the number of degrees of freedom per node. For 3D elements, each node usually has three degrees of freedom. In addition, represents the rigid body displacement degrees of freedom of the element. For example, a hexahedral Voronoi element with 8 nodes has 24 nodal degrees of freedom and 6 rigid body displacements, so the minimum value is 18.

[0217] Example 2

[0218] In this embodiment, the final expression derived and corrected in the above embodiment 1 is further verified, and the specific verification process is as follows:

[0219] As can be seen from Example 1, to analyze the high-order stress field in three-dimensional space, [Δβ m ]The minimum order required is 6. However, to determine the optimal stress order, it must be verified that [Δβ] pol and [Δβ] rec Convergence of the selected order. , for and The convergence of the selected options was independently verified. The order of is selected between 16 and 161 (too high an order may result in excessive stiffness, thus distorting the calculated value), and for The order of is chosen between 16 and 52. Taking the plasticity of a single hole as an example, the values ​​calculated by MSC MARC are used as reference to verify the effects of different orders in 3D VCFEM. and The convergence of Figure 3 As shown in (a), when REC = 31, the values ​​of POL are shown in Table 1. It can be seen that when POL = 31 and POL = 52, there is a big difference compared with the MSC MARC results, indicating that the stress results obviously lack convergence under these orders. Therefore, after eliminating these two orders, the stress comparison results are shown in Figure (3).

[0220] Table 1 Relative error of the maximum stress concentration point

[0221]

[0222]

[0223] *Among them, the “-” sign indicates that the relative error is too large

[0224] To quantify the specific difference in stress values, the relative error is defined as follows:

[0225]

[0226] Among them, σ m represents the stress value calculated by MSC MARC, σ ij Indicates the stress value calculated by 3D VCFEM at POL = i and REC = j. In order to accurately capture the difference between 3D VCFEM and MSC MARC at the stress concentration point, only the two symmetrical points closest to the void on the Y-axis path are selected as sampling points, and the points that are far away on the Y-axis path and cannot reflect the difference between MSC MARC and 3D VCFEM are not selected as sampling points.

[0227] from Figure 3From the results, it can be seen that at the element boundary far away from the cavity, POL=116 and POL=80 have little difference compared with the results of MSCMARC. However, at the maximum stress concentration point near the cavity, POL=116 shows the smallest relative error compared with MSC MARC. The detailed comparison is shown in Table 2;

[0228] Table 2 Material properties

[0229]

[0230] “-” indicates that in the 3D VCFEM calculation, the relative error is too large compared with the MSC MARC calculation value, resulting in no convergence. Figure 3 As shown in (b), the reason why two relative error values ​​appear at POL=80 is that the stress values ​​of the two symmetrical points near the hole are asymmetric, resulting in two different relative error values. According to Table 3, at the location of maximum stress concentration, the relative error between 3D VCFEM and MSC MARC does not exceed 10%, while at the location far away from the void, the relative error between 3D VCFEM and MSC MARC does not exceed 3%. This fully proves the convergence of the high-order stress function constructed under the conditions of POL=116 and REC=31.

[0231] Example 3

[0232] Regarding the plastic strain of a single hole in a cube, the model structure is as follows Figure 4 (c) The model is centered at point (0, 0, 0). The model size is Figure 4 (a) and Figure 4 (b). The displacement boundary conditions of each increment are shown in Table 2, and the elastic-plastic parameters of the material are shown in Table 1. The number of increments is set to 25, and thermal strain is not considered in this numerical example. Since the values ​​of POL and REC have been determined in Example 1 and the relative error of the stress value has been given, in this numerical example, POL is set to 116 and REC is set to 31, and only the stress distribution diagram is displayed.

[0233] The stress distribution slices calculated by MARC and VCFEM are shown in Figure 2. Figure 5 (a) and Figure 5 (b) is shown. The normal vector of the slice plane is (1, 0, 0) and passes through the point (0, 0, 0). By observing the stress concentration and the distribution of stress contours at the 25th increment, it can be concluded that the results of VCFEM and MARC are basically consistent, thus verifying the effectiveness of VCFEM. The stress values ​​of MARC and VCFEM along the Y-axis path are extracted and compared. The results are shown in Figure 3 (b) as shown.

[0234] like Figure 3As shown in (b), the stress variation trends of VCFEM and MARC are basically the same, proving the accuracy of VCFEM. This confirms that VCFEM can achieve the same level of accuracy as the finite element method (FEM) under the condition of material plastic strain.

[0235] Example 4

[0236] There is a temperature-dependent thermal expansion of a single hole in a cube, the dimensions and boundary conditions of the model are as follows Figure 4 (a) Figure 4 (b) and Figure 4 As shown in (c), the center of the model is at (0,0,0). Table 1 lists the detailed elastic-plastic parameters of the matrix material. This example only considers the thermal expansion that changes with temperature in a heating environment, and the final steady-state temperature is 530°C. The specific average thermal expansion coefficient is given in Example 1. The specific displacement boundary conditions and increments of the model are shown in Table 2. Since the influence of the selected polynomial order on the results in the three-dimensional VCFEM has been verified, and the relative errors have been compared, and the element structure remains unchanged, the polynomial order used in the calculation is the same as that verified in the previous section. Figure 6 (a) and Figure 6 As shown in (b), the stress results of the last incremental step are sliced. The normal vector of the slice is (1,0,0) and passes through the center of the model (0,0,0).

[0237] By observation Figure 6 (a) and Figure 6 From the results in (b), it can be concluded that when considering temperature-dependent thermal expansion, the stress concentration trend in 3D VCFEM is basically consistent with that calculated by MSC MARC, which verifies that 3D VCFEM can accurately capture the stress behavior of porous materials under temperature-dependent thermal expansion. Figure 6 The stress values ​​extracted by the red path in (c) also generate a stress value comparison graph, such as Figure 6 (d) As shown in Figure 96, the stress values ​​near the hole obtained by 3D VCFEM are slightly higher than those obtained by MSC MARC. This is because more polynomial terms are used in the interpolation function to improve the calculation accuracy. In addition, under the same model and working conditions, the number of elements and nodes used by 3D VCFEM is significantly less than that of MSC MARC, so the wall clock time required for 3D VCFEM to complete the numerical simulation of 10 increments is much shorter.

[0238] In summary, under the condition that the thermal expansion coefficient is temperature-dependent, 3D VCFEM significantly improves the computational efficiency while maintaining a certain accuracy, making it a promising application tool in certain industrial fields that may involve extreme high temperature environments in the future.

[0239] Example 5

[0240] In this embodiment, there are regularly distributed holes in the cube. Regarding the temperature-dependent thermal expansion, this study constructed a porous model containing 27 voids to study the temperature-dependent thermal expansion effect at high temperatures. The high-temperature thermal expansion coefficient of the porous material is given by equation (87). The displacement boundary conditions of each increment are shown in Table 2, 10 increment steps.

[0241] The calculation results are obtained through the points (0, 0, 0), (20, 0, 0) and (-20, 0, 0), and multiple slices are made with the vector (1, 0, 0) as the normal. The stress contour diagram of the slice is shown in Figure 6 (a) and 6(b). Figure 7 The red line path in (c) extracts stress data to form a stress path. The results are compared in Figure 7 (d).

[0242] The stress calculation results of multiple elements show that the stress results obtained by 3D VCFEM and MSCMARC are very consistent under heating conditions from 20°C to 530°C. This shows that VCFEM not only achieves good agreement with MSC MARC under low-temperature linear expansion conditions, but also has high accuracy in stress calculation results under high-temperature conditions.

[0243] Regarding plastic strain, Figure 4 As shown in (d), a porous model with 27 pores is constructed under the consideration of plastic strain. The elastic-plastic material parameters are listed in Table 2. 2. The displacement boundary conditions of each increment are shown in Table 2, and there are 25 increments in total.

[0244] The calculation results are obtained through the points (0, 0, 0), (20, 0, 0) and (40, 0, 0), and the vector (1, 0, 0) is used as the normal line to perform multi-slice calculations. The stress contours of the slices are shown as follows: Figure 8 As shown in (a) and (b). Figure 8 The red line path in (c) extracts the stress contour line slice results to form the stress path. The results are compared in Figure 8 (d).

[0245] The calculation results of multiple elements show that the stress calculation results of 3D VCFEM and MSC MARC are very consistent under the Mises yield criterion and linear isotropic hardening conditions. This shows that VCFEM can achieve high accuracy not only under thermoelastic conditions but also under plastic conditions. In addition, since the number of nodes and elements used by 3D VCFEM is much smaller than that of MSCMARC, the calculation efficiency is significantly improved. This is also one of the advantages of VCFEM in analyzing nonlinear problems of models containing multiple microstructures.

[0246] In this example, there are randomly distributed holes in the cube. The above numerical examples verify the effectiveness and efficiency of 3D VCFEM in solving nonlinear strain problems with regularly distributed pores. However, under realistic conditions, the pores may be randomly distributed in the model. In order to simulate the mechanical behavior of porous materials under realistic conditions (random distribution), we used 3D VCFEM to construct models with 40 randomly distributed pores and 90 randomly distributed pores, respectively. Fig. 9 (a) and Fig. 9 (c) is shown. Fig. 9 In the model (a), 25 increments of plasticity are considered, while Fig. 9 In the model of (c), only the temperature-dependent thermal expansion at high temperatures is considered, with an increment of 10. Fig. 9 (b) and Fig. 9 (d) shows the stress distribution diagrams of the 25th and 10th increments, indicating that the 3D VCFEM results have good regularity for the model with randomly distributed pores.

[0247] In summary, the present invention discloses a modified complementary energy functional representation method based on the three-dimensional Voronoi unit finite element method (3D VCFEM), which takes into account plasticity and temperature-dependent thermal expansion at high temperatures. Based on the principle of minimum complementary energy, a modified incremental complementary energy functional is derived for the case of thermal strain. After scaling the unit, each unit is subdivided into multiple Delaunay tetrahedrons, and the Hammer numerical integration method is used for integration to construct a three-dimensional stress function considering the ellipsoid. Compared with traditional finite element analysis, three-dimensional VCFEM requires less calculation time for each incremental step in nonlinear simulation. The same model and load condition only requires less than 600 seconds, and the wall clock time is reduced by nearly 54 times. While maintaining accuracy, efficiency is greatly improved.

[0248] In addition, through further optimization of this method, the three-dimensional VCFEM not only achieves the same accuracy as the traditional finite element method (FEM) while using fewer nodes and elements, but also highlights that the correction analysis performed by this method has a higher computational efficiency than the existing finite element analysis and has outstanding effects; it can simulate the nonlinear behavior of porous materials in complex environments, which is a significant improvement.

[0249] The preferred embodiments of the present invention disclosed above are only used to help explain the present invention. The preferred embodiments do not describe all the details in detail, nor do they limit the invention to only the specific implementation methods described.

Claims

1. A three-dimensional VCFEM numerical simulation method for porous materials considering temperature and plasticity, characterized in that: The following steps are involved: S1: According to the existing initial complementary energy functional expression of the stress hybrid element model, according to the relationship between temperature, thermal strain and strain, and considering the boundary adjustment according to the plasticity increment theory, the complementary energy functional can be rewritten as follows: In the formula, Ω m is the area of ​​the matrix material, Γ e is the displacement boundary of the matrix material, Γ t is the force boundary of the whole unit, Γ i is the interface between the matrix material and the pores; S2: Each increment is equal to the sum of the total increment of the previous step plus the increment of the current step, and the total energy of the system is obtained; at the same time, the motion relationship within the unit is obtained according to the Euler equation, which can be expressed as: Where, Δσ is the stress increment in the unit domain; Δε is the strain increment in the unit domain; Δu is the node displacement increment on the unit boundary; Δu1 is the node displacement increment on the matrix and hole interface close to the matrix side; Δu2 is the node displacement increment on the matrix and hole side close to the hole side; Δε is e —elastic strain increment within the unit; Δε p —Plastic strain increment within the element; Δε th —thermal strain increment within the element; S3: The total energy of the system Π mc The first-order variation of the displacement increment Δu can obtain the corresponding stress coordination relationship. At the same time, according to the stress function of the matrix material in the three-dimensional porous unit, the Maxwell stress function of the stress field and the interaction function; and the displacement increment of the unit boundary is expressed as the interpolation of the generalized displacement increment of each boundary segment node; the weak form of the system complementary energy after discretization can be expressed as the following formula: For Δβ m Performing first-order variational calculus can obtain the stationary value condition of the modified complementary energy functional: And can get To the weak form of the motion relationship within each unit; at the same time, in the iterative solution process, the iterative expression of the i-th iteration step can be expressed as: [dβ mi ]=-[H m ] -1 [H m ] th +[H m ] -1 [G][dq](5)where: S4: The first-order variation of the total capacity of the system with the node displacement increment obtains the weak form of the surface condition to be solved: Then substitute formula 5 into formula 6 to obtain: The solution of the matrix is: in, According to equations (8) and (9), we can get the solution relationship as follows: in: According to the above formula, the solution relationship can be obtained again, and the internal node displacement increment, the unit solution relationship, and the equivalent node force increment can be obtained; finally, it can be obtained: S5: Then the elastic-plastic constitutive equation of thermal strain is considered, and further derivation is carried out according to the basic laws of plastic mechanics. At the same time, this method is based on the small deformation thermal elastic-plastic constitutive relationship, considering the Mises yield condition and isotropic strengthening material; the specific basic laws of plastic mechanics include initial yield condition, plastic flow law, and plastic hardening law; The specific steps include: S5.1: Using the Mises yield condition, when the second stress invariant J2 = C (C is the material constant), the material will enter the yield state. The value of the material constant C can be obtained through uniaxial tension test. According to the expression of J2, the Mises yield condition can be rewritten as: S5.2: According to the Mises plastic potential concept expression, it is transformed into a flow law expression: in, The defined vector is along the normal direction of the subsequent yield surface F = 0 in space; When the yield condition is the Mises yield condition, substituting the Mises yield function equation (the expression of J2) into equation (17) yields: Then rewrite the tensor expression (17) into a general expression, transform the equivalent stress, and then substitute (18) into the general equation to obtain: S5.3: The plastic hardening law considers the process as follows: The hardening law stipulates that the material enters the loading surface determined by the subsequent yield function after plastic deformation, so the general loading function is used: Among them, C is the hardening parameter. For the 2024-T6 aluminum alloy material used in this paper, since there is no reverse loading and unloading process, it will not reflect the Bauschinger effect, so the linear isotropic hardening law is adopted; then the Mises yield condition is adopted, and the subsequent yield function of isotropic hardening can be expressed as: Among them, σ ij is the elastic-plastic stress, which is the equivalent plastic strain The function of It can be obtained from the σ-ε curve of the uniaxial tensile test of the material; definition: Among them, H is the plastic modulus, also known as the hardening modulus, then H is related to the elastic modulus E and the tangent modulus The relationship is: S6: According to the Parandtl-Reuss theory, the strain increment (elastic strain increment and plastic strain increment) is considered and further combined with the Mises yield law and hardening law to derive the expression: At the same time, the unit flexibility matrix is ​​established for the three-dimensional problem, and the final expression is established as follows: S7: In order to determine whether plastic loading or elastic unloading should be performed in 3D VCFEM, and thus decide whether to use the plastic constitutive relation or the elastic constitutive relation, the following formula is established for judgment: in It is a material parameter determined by the yield function. When condition (1) in the equation is met, plastic loading is continued. When condition (2) is met, elastic loading is selected. When condition (3) is met, for ideal elastic-plastic materials, new plastic flow is generated. For hardened materials, no new plastic flow occurs.

2. According to claim 1, a three-dimensional VCFEM numerical simulation method for porous materials considering temperature and plasticity is characterized in that: In S1, the initial complementary energy functional expression of the existing stress hybridization metamodel is: Where S is the elastic-plastic flexibility matrix, is the given boundary condition, T is the surface force of the current boundary, Γt e is the displacement boundary of the unit; The expression of thermal stress is: Among them, α(T) is the temperature-dependent thermal expansion coefficient, and α(T) can be defined as a function related to the temperature field variables.

3. According to claim 2, a three-dimensional VCFEM numerical simulation method for porous materials considering temperature and plasticity is characterized in that: Furthermore, in S2, each increment is equal to the sum of the total increment of the previous step plus the increment of the current step: The total energy of the system is expressed as:

Images