Modeling method of piezoelectric particle composite material interface bonding unit model

Abstract

The invention relates to the technical field of mechanical-electric coupling, in particular to a modeling method of a piezoelectric particle composite material interface bonding unit model, which comprises the following steps of: S1, eliminating strain and electric field intensity in a functional by introducing a mixed complementary energy functional to obtain a simplified equation; s2, selecting a Voronoi grid to carry out domain division, dividing the material into a plurality of Voronoi units, and enabling the boundary of each unit to be composed of inclusions and central bisectors of adjacent inclusions; s3, adopting a Lagrange multiplier method to realize the introduction of an interface condition; and S4, through comparison with ABAQUS commercial finite element software, verifying the accuracy and reliability of the model in practical application. According to the method, the material is divided into Voronoi units, so that the generation amount of grids is greatly reduced. In addition, the Voronoi unit can naturally adapt to irregular distribution of inclusions and matrixes in the material and accurately describe the interface position and complex geometric characteristics of the interface position.

Description

Technical Field

[0001] The present invention relates to the technical field of force - electricity coupling, and particularly relates to a modeling method for an interfacial bonding unit model of a piezoelectric particle composite material. Background Art

[0002] Piezoelectric material structures have a unique force - electricity coupling effect, and are widely used in intelligent structures due to their high precision and high sensitivity. They can be used for the production of both actuators and sensors, and are continuously applied in fields such as manufacturing, automotive, aerospace engineering, medical instruments, information communication, etc.

[0003] Monolithic piezoelectric ceramics and piezoelectric polymers have their own limitations. Therefore, the idea of combining the two to make up for each other has led to the emergence of piezoelectric composites. Piezoelectric ceramic particle / polymer composites are piezoelectric composites formed by filling piezoelectric ceramic particles in a polymer matrix. Its dielectric constant is extremely low, but its force - electricity coupling coefficient is very high, several times higher than that of pure piezoelectric ceramics, and its flexibility is much better than that of piezoelectric ceramics. Therefore, the comprehensive performance of piezoelectric composites is superior to that of average piezoelectric materials and has great application value.

[0004] Under the same deformation, piezoelectric materials may exhibit higher stress due to the generation of charges and the interaction of electric fields. This is because the piezoelectric effect causes the material to generate additional charges when subjected to stress, and these charges interact with the external electric field, further increasing the internal stress of the material. People use various methods, including analytical, experimental, and numerical methods, to study the mechanical properties of piezoelectric structures. Researchers have widely explored the static, dynamic, and electromechanical coupling characteristics of piezoelectric structures using various finite - element methods.

[0005] It is also very important to understand the local behavior of piezoelectric composites and structures. Local stress concentration is one of the main factors that may cause material damage. In particular, brittle piezoelectric components are mostly used as particles embedded in a polymer matrix. In many cases, fractures occur at the interfaces, such as debonding and cracks. The appearance of interface cracks will reduce the electromechanical conversion rate of the piezoelectric structure, and piezoelectric composites are extremely prone to interface debonding and premature fracture failure. Therefore, further understanding of the interface fracture mechanism will contribute to the design and safe application of related piezoelectric composites and intelligent devices. In the study of interface fracture in piezoelectric composites, Parton provided a reasonable method for determining the critical load conditions that lead to crack propagation. Govorukha and Kamlah described the complex singular fields at the tips of interface cracks between various piezoelectric bodies by combining finite element and asymptotic techniques. Wang et al. used the finite element method to capture the initiation and propagation of interface fracture behavior in the structure of piezoelectric composites. However, although simple finite elements can model the global and local behavior of piezoelectric composites and structures, it is very inefficient and requires a large amount of computational time, resources, and manpower when generating complex meshes. Therefore, most simple finite element analyses in the literature use simple elements with only one inclusion, which oversimplifies the complex material microstructure.

[0006] For the case of randomly distributed large-scale inclusions, traditional displacement-based finite elements are difficult to simulate. A large number of fine meshes need to be divided at the interfaces between inclusions and the matrix, restricting the computational region to a very small area. The present invention proposes a method VCFEM to improve the computational efficiency for simulating piezoelectric composites. Compared with traditional finite element models, this method has a small amount of mesh generation and high computational accuracy. The Voronoi cell finite element method proposed by Ghosh and Moorthy has been used as an efficient and accurate tool. The Voronoi mesh diagram can accurately describe the size, shape, and randomness of the spatial distribution of the composite material microstructure, enabling users to better simulate the shape and randomness of the spatial distribution of inclusions, voids, and / or microcracks in the composite material. Summary of the Invention

[0007] The purpose of the present invention is to provide a modeling method for an interface bonding element model of a piezoelectric particle composite material. By dividing the material into Voronoi cells, the amount of mesh generation is significantly reduced. In addition, Voronoi cells can naturally adapt to the irregular distribution of inclusions and the matrix in the material and accurately describe the interface position and its complex geometric characteristics.

[0008] To achieve the above technical purposes and reach the above technical effects, the present invention is realized through the following technical solutions:

[0009] A modeling method for an interface bonding element model of a piezoelectric particle composite material, comprising the following steps:

[0010] S1: Construct a hybrid energy functional, which is described by the generalized variational principle, dealing with six independent variables including stress, strain, displacement, electric displacement, electric field strength and electric potential. By introducing a hybrid complementary energy functional, strain and electric field strength are eliminated in the functional to obtain a simplified equation;

[0011] S2: Select Voronoi meshes for domain division, dividing the material into multiple Voronoi cells, such that the boundary of each cell is composed of the bisector of the inclusion and the center of adjacent inclusions;

[0012] S3: Adopt the Lagrange multiplier method to introduce interface conditions, ensuring the balance of stress and electric displacement; handle rigid body motion by adding four constraint conditions, thereby ensuring the non - singularity and numerical stability of the column matrix.

[0013] S4: Verify the accuracy and reliability of the model in practical applications by comparing it with the commercial finite element software ABAQUS. In different examples, test the single - inclusion and multi - inclusion cases to ensure that the model can correctly simulate the multi - field coupling behavior of piezoelectric composites.

[0014] Advantages of the present invention:

[0015] The present invention adopts the Voronoi cell finite element method (VCFEM) to model the microstructure of materials through the characteristics of Voronoi diagrams. Traditional finite element methods usually require generating a large number of fine meshes when simulating complex multiphase materials, which not only increases the computational amount, but also easily leads to accuracy loss when dealing with heterogeneous materials at different scales. VCFEM divides the material into Voronoi cells, significantly reducing the amount of mesh generation. In addition, Voronoi cells can naturally adapt to the irregular distribution of inclusions and matrix in the material and accurately describe the interface position and its complex geometric characteristics. Since this modeling method can directly reflect the microstructure characteristics of materials, it shows high accuracy in capturing local stress concentration phenomena. This method solves the bottleneck of traditional methods in terms of computational efficiency and accuracy, especially having significant advantages in large - scale calculations.

[0016] The characteristics of piezoelectric materials result in a complex coupling effect between the force field and the electric field, making it difficult for traditional finite element methods to effectively capture the characteristics of these multi-physical fields simultaneously. The method of the present invention establishes a hybrid energy functional, taking stress, strain, displacement, electric displacement, electric field strength, and electric potential as independent variables, and by introducing the Lagrange multiplier method, ensuring the continuity and balance of interface stress and electric displacement. This processing method not only solves the computational instability caused by the complexity of the material interface but also enhances the accurate simulation ability of the force-electric coupling effect. In other words, the model can finely describe the local stress and charge changes in the material due to the piezoelectric effect, providing a solid theoretical basis for the study of multi-field coupling problems.

[0017] In practical applications, piezoelectric composites often face the challenge of complex structures with multiple inclusions. The modeling method of the present invention not only exhibits superior performance in simple structures with a single inclusion but also can handle complex structures with multiple inclusions of different shapes and sizes coexisting. This benefits from the flexibility of Voronoi cells, which can adapt to different boundary conditions and material inhomogeneities. Through comparison and verification with commercial finite element software such as ABAQUS, the simulation accuracy and stability of this method under various complex working conditions are confirmed. This verification method fully demonstrates the reliability of the method, enabling it to be promoted and applied in engineering practice, especially in fields such as aerospace and automotive engineering that require precise structural design and verification, showing extremely high application value.

[0018] Of course, it is not necessary for any product implementing the present invention to simultaneously achieve all the above-mentioned advantages. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for describing the embodiments will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on these drawings.

[0020] Figure 1 Schematic diagram of a two-dimensional Voronoi piezoelectric unit model with one inclusion;

[0021] Figure 2 Schematic diagram of a two-dimensional Voronoi unit model with one inclusion considering interface delamination;

[0022] Figure 3 Schematic diagram of the model loading condition;

[0023] Figure 4 VCFEM mesh diagram and ABAQUS mesh diagram;

[0024] Figure 5 is σ x Stress spectrogram;

[0025] Figure 6 is σ y Stress spectrogram and D y Electric displacement spectrogram;

[0026] Figure 7 is the stress path selection diagram;

[0027] Figure 8 is the grid diagram for verifying VCFEM and ABAQUS with multi-inclusion elements;

[0028] Figure 9 is σ x Stress and D y Electric displacement nephogram;

[0029] Figure 10 is σ x Path diagram;

[0030] Figure 11 is the simulation grid diagram of multi-inclusion piezoelectric elements;

[0031] Figure 12 is σ x Stress spectrogram and D y Electric displacement nephogram. Specific implementation manners

[0032] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts belong to the scope of protection of the present invention.

[0033] Embodiment 1

[0034] A modeling method for an interface bonding unit model of a piezoelectric particulate composite material described in this embodiment includes the following steps:

[0035] S1: Construct a mixed energy functional, which is described by the generalized variational principle, and processes six independent variables of stress, strain, displacement, electric displacement, electric field strength, and electric potential. By introducing a mixed complementary energy functional, strain and electric field strength are eliminated in the functional to obtain a simplified equation;

[0036] S2: Select Voronoi grids for domain division, divide the material into multiple Voronoi units, and make the boundary of each unit consist of the inclusion and the bisector of the centers of adjacent inclusions;

[0037] S3: The Lagrange multiplier method is adopted to introduce the interface conditions and ensure the balance of stress and electric displacement. Four constraint conditions are added to handle the rigid body motion, thus ensuring the non-singularity and numerical stability of the column matrix.

[0038] S4: By comparing with the commercial finite element software ABAQUS, the accuracy and reliability of the model in practical applications are verified. In different examples, the single inclusion and multi-inclusion cases are tested to ensure that the model can correctly simulate the multi-field coupling behavior of piezoelectric composites.

[0039] Example 2

[0040] Hybrid complementary energy functional of piezoelectric particulate composites

[0041] The most general variational principle is the generalized variational principle including six independent variables of stress, strain, displacement, electric displacement, electric field strength and electric potential. The fixed conditions of this functional are equivalent to all the governing equations and boundary conditions of the piezoelectric elastomer. Existing hybrid energy functionals:

[0042]

[0043] In this example, are used to represent the specified displacement at the displacement boundary and the specified traction force at the external force boundary respectively. is used to represent the specified electromotive force at the electric potential boundary and the specified surface charge density at the electric charge density boundary. In the finite element analysis method, both the mechanical and electrical boundary conditions can satisfy the displacement and electric potential as nodal variables. Through the two constraint conditions u and t it can be simplified to: represent the specified electromotive force at the electric potential boundary and ω the specified surface charge density at the electric charge density boundary. In the finite element analysis method, both the mechanical and electrical boundary conditions can satisfy the displacement and electric potential as nodal variables. Through the two constraint conditions and it can be simplified to:

[0044]

[0045] Assumptions are made for the field function variables of stress σ, displacement u, electric displacement D, and electric potential φ. Using the piezoelectric constitutive relations σ = C E γ - e T E, D = eγ - ε γ E, the strain γ and the electric field strength E are eliminated in the functional to obtain a simplified equation, denoted as Π σD :

[0046]

[0047] Assume that the piezoelectric material is two-dimensional isotropic. Under the plane stress state, the constitutive relation equation can be obtained:

[0048]

[0049] Wherein, c, e, and d are the elastic constant, piezoelectric coefficient, and dielectric constant, respectively.

[0050] If we want to model the microstructure containing inclusions, a simple finite element program requires a large amount of time and resources for calculation and mesh generation. Therefore, it is decided to develop a new VCFEM for more efficient simulation of piezoelectric composites.

[0051] Inside the piezoelectric composite, the topological structure is quite complex, and the inclusions are of various sizes and randomly distributed. The deformation and stress at each point in the matrix are affected by these complex topological structures of the material. The greatest influence comes from the inclusion closest to that point. Therefore, in this embodiment, all the points closest to the center of the inclusion form an aggregate to constitute a unit. In this way, the edges of the unit are formed by the central bisectors of this inclusion and adjacent inclusions. Mathematically, the mesh formed by such edges is called a Voronoi mesh diagram. Therefore, the unit composed of the points closest to the inclusion is called a Voronoi unit.

[0052] Now, considering the characteristics of the Voronoi unit, the domain Ω is discretized into sub-domains Ω e , so that Ω = Σ e Ω e . At the same time, there is an inclusion Ω e in each Ω c , which satisfies Use Ω m to represent the matrix material, that is, Ω e = Ω m + Ω c . In order to develop a new VCFEM, consider the complete functions σ e and D i independently assumed within each sub-domain Ω i . In each discretized Voronoi unit domain, there are nodes on the outer and inner boundaries. Figure 1 shows a Voronoi unit containing a ten-node inclusion. Use D c to represent the fields in and respectively. Use Φ m , Φ c to represent the displacements and electric potentials on the boundaries and respectively.

[0053] The following will pre-satisfy the essential boundary conditions and the The interface between the matrix and the inclusion The force continuity and electric displacement continuity conditions

[0054] Voronoi cell construction for particle delamination under the force - electricity coupling effect

[0055] Refinement of the derivation of the piezoelectric element formula

[0056] In order to analyze the influence of the delamination between the inclusion and the matrix in the piezoelectric particle - reinforced composite on the structural micro - evolution and macroscopic properties, in this embodiment, based on the stress hybrid element, the Voronoi cell of the interface delamination in the force - electricity coupling field is derived

[0057] Consider a two - dimensional Voronoi cell model of an inclusion delamination as Figure 2 shown. In the force - electricity coupling region and on the boundary, the following should be satisfied

[0058]

[0059] Based on the above - derived mixed complementary energy functional formula, applying the Lagrange multiplier method to implement the above - mentioned constraint conditions, the modified mixed complementary energy functional is obtained

[0060]

[0061] Among them,

[0062] The boundary satisfies

[0063] The displacements and electric potentials of each boundary are

[0064]

[0065] The boundary displacement / electric potential d is the interpolation of the generalized nodal displacement / electric potential q

[0066] d = Lq (18)

[0067] Assumed independent stress fields / electric displacement fields are introduced into the Voronoi cell. Their construction should satisfy the equilibrium conditions, and the equilibrium conditions can be obtained through the stress function / electric displacement function. In this embodiment, the stress (electric displacement) functions of the matrix and the inclusion are respectively

[0068]

[0069] Among them, and are the Airy stress functions of the matrix and inclusion parts. Considering the existence of the electric displacement field, and are defined as the electric displacement functions corresponding to the stress functions and are the interaction functions of stress and electric displacement, respectively, which consider the complete polynomials based on the reciprocal of the shape function at the matrix-inclusion interface. The detailed construction of these functions can be found in reference

[29] . In this embodiment, by differentiating the stress function and the electric displacement function, the stress and electric displacement at each point of the element can be obtained:

[0070]

[0071] where is the stress matrix of the matrix and the inclusion, is the electric displacement matrix of the matrix and the inclusion, is the stress parameter of the matrix and the inclusion, is the electric displacement parameter of the matrix and the inclusion.

[0072]

[0073] After finite element discretization and substitution and arrangement, we get

[0074]

[0075] where

[0076]

[0077]

[0078] Performing variational on equation (24), according to the stationary value principle, from we can obtain

[0079]

[0080] can be simplified to β = H -1 Gq (34)

[0081] Substituting equation (34) into equation (24), according to the stationary value principle, from we can obtain the system of equations for solving the generalized displacement / electric potential

[0082]

[0083] where the element stiffness matrix is equal to:

[0084]

[0085] Elimination of the rigid body displacement / electric potential of the inclusion

[0086] Since there is no connection between the internal and external nodes of the element, it is not possible to fully establish a connection between the interface node displacements / potentials and the displacements / potentials outside the element. Four additional displacement constraints must be implemented to ensure the consistency of the rigid body displacements / potentials of the inclusion and the matrix. The rigid body displacements / potentials of the element boundary nodes can be expressed as

[0087]

[0088] where x i , y i are the coordinates of node i, α 1 , α 2 are the rigid body translational momenta and α 3 represents the corresponding potential constant, and α 4 is the degree of freedom in the rigid body rotation direction. The interface rigid body node displacements / potentials can be represented by the matrix φ'. Then the node displacements / potentials of the interface and the element boundary can be expressed as

[0089] dq' = φ'α + dq' def (38)

[0090] dq = φα + dq def (39)

[0091] where dq' def and dq def are the pure deformation displacements / potentials of the interface and element boundary nodes respectively. Since the pure deformation mode and the rigid body mode spaces are orthogonal, multiplying both sides of Eqs. (38) and (39) by φ' T and φ T gives:

[0092] φ' T dq' = φ' T φ'α (40)

[0093] φ T dq = φ T φα (41)

[0094] In the iterative solution process, for the same element, the iterative increment of the element boundary node displacements / potentials is dq ei and the iterative increment of the node displacements / potentials on the matrix side of the matrix-inclusion interface is equal to

[0095]

[0096] The rigid body displacements of the two parts are the same, so it satisfies:

[0097]

[0098] where, Φ = [{φ T φ} -1 φ T -{φ' T φ'} -1 φ' T

[0099] Equation (43) provides four additional constraint conditions to exactly eliminate the rigid body displacements / electric potentials of the inclusions within the element, ensuring that the element stiffness matrix K e is non-singular. Introducing the above constraint conditions into the modified complementary energy functional by the Lagrange multiplier method gives

[0100]

[0101] According to the principle of stationary complementary energy functional, we obtain

[0102]

[0103] Similarly, according to we obtain

[0104]

[0105] Combining equations (45) and (46), we obtain the solution equation

[0106]

[0107] where the stiffness matrix is equal to

[0108]

[0109] Condensation of internal degrees of freedom

[0110] From the structure of the element, it can be known that the nodal displacements / electric potentials on the internal interfaces of the element are only related to the external nodal displacements / electric potentials of the same element, and have no direct relationship with the nodes of other elements. To establish the connection, rewrite equation (47) as:

[0111]

[0112] where

[0113]

[0114]

[0115] Divide equation (49) into two parts. The first group only involves the external nodal displacements / electric potentials of this element and the external nodal displacements / electric potentials of other elements. Thus, we have

[0116] ​

[0117] The second group only involves the relationship between the displacements / potentials of the internal nodes and those of the external nodes within each element, satisfying:

[0118]

[0119] That is

[0120]

[0121] Substituting (58) into (56) gives

[0122]

[0123] At this time, the expression of the element stiffness matrix is

[0124]

[0125] In this way, solving the relation (60) only involves the nodal displacements / potentials outside the element, and the condensation of the internal degrees of freedom greatly reduces the global stiffness matrix.

[0126] Embodiment 3

[0127] Verification of single / multi-inclusion piezoelectric elements

[0128] To verify the effectiveness and accuracy of the VCFEM considering the force-electric coupling, models with the same structure and boundary conditions were established using the commercial finite element analysis software ABAQUS and VCFEM. It is assumed that the 10th-order Airy stress function and its interaction function are assumed to be of the 8th order. It is assumed that the 8th-order electric displacement function and its interaction function are assumed to be of the 4th order. And the three-node piezoelectric triangular element (CPS3E) is used in ABAQUS.

[0129] In this embodiment, single-inclusion and multi-inclusion models are used to simulate the stress and inclusion delamination phenomena of composite piezoelectric materials under the force-electric coupling field. For more concise and intuitive analysis in the analysis, the physical quantities are dimensionless. Considering that the x-y plane is the isotropic plane of the piezoelectric material, both the matrix and the inclusion materials have strong piezoelectric properties, and the material parameters of the matrix and the inclusion are shown in Table 1.

[0130] Table 1 Dimensionless material parameter table

[0131]

[0132] Verification of a completely bonded single piezoelectric element

[0133] Select the model as a circular inclusion with a radius of 0.2 in the middle of a 2×2 square matrix. For easy comparative analysis, the working conditions of the VCFEM model and the ABAUS model are exactly the same. The model diagram is asFigure 3 As shown, the displacements of the top and bottom boundaries are 0, the horizontal displacement of the left boundary is 0, and the horizontal displacement of the right boundary is 0.0002.

[0134] In this example, the VCFEM model contains a simple Voronoi cell, while in ABAQUS, the same model is divided into 19378 cells. The mesh diagrams of VCFEM and ABAQUS are as Figure 4 shown.

[0135] By observing Figure 5 and Figure 6 the stress curve diagrams, it can be seen that the results of VCFEM and ABAQUS are highly consistent, which proves the feasibility of this theory and verifies the effectiveness of the single inclusion model under the electro-mechanical coupling effect. Two special paths (the horizontal midline and the vertical midline) are selected to extract the stress in the x direction and the stress in the y direction as Figure 7 shown, and the fitting curve is compared with the ABAQUS result.

[0136] As Figure 7 shown, the change trend of VCFEM is basically consistent with the result of ABAQUS, which proves the accuracy of VCFEM in the single inclusion model under the electro-mechanical coupling effect.

[0137] Verification of multi-inclusion cells

[0138] Since the number of sides of the Voronoi cell is not fixed. To verify the effectiveness of multiple Voronoi cells with different numbers of sides in the electro-mechanical coupling field, a model with 5 Voronoi cells is established in this embodiment for comparison with the ABAQUS model.

[0139] Five elliptical inclusions with different shapes are randomly added to a 2×2 square matrix, and the range of the semi-major axis length of the inclusions is 0.070 - 0.098. The critical normal stress value on the cell interface is set to 13631608. A displacement load is continuously applied to the right boundary of this model, and the other boundaries are fixed. After calculation by VCFEM, the inclusion cell will experience the first delamination when the displacement load is 0.0004.

[0140] In this example, the VCFEM model contains 5 simple Voronoi cells, while in ABAQUS, the same model is divided into 29842 cells. The mesh diagrams of VCFEM and ABAQUS are as Figure 8 shown.

[0141] Figure 9The stress nephogram of 5 inclusion units, with the displacement load being 0.0001 at this time. It can be seen that when simulating multiple Voronoi units with different numbers of sides, the calculation results are in good agreement with ABAQUS. This proves that VCFEM can still ensure the calculation accuracy when calculating complex cases with multiple inclusions.

[0142] Φ poly and Φ rec Influence on calculation accuracy

[0143] The correct selection of the number of terms of the stress function and the number of terms of the electric displacement function has an important impact on the convergence and efficiency of VCFEM. Since stress concentration generally occurs near the matrix-inclusion interface, in this embodiment, this embodiment will discuss Φ poly and Φ rec The influence of the increasing order of different combinations.

[0144] The selected model is a circular inclusion with a radius of 0.2 in the middle of a 2×2 square matrix. The material parameters and loading conditions are the same as above, and the stress on the horizontal center line is extracted. Different numbers of terms in formulas (19) and (20) are used in the calculation. In this embodiment, the calculation results of VCFEM are compared with the calculation results of the finite element software ABAQUS to ensure the correctness of the calculation results.

[0145] If the total number of system degrees of freedom is assumed to be n and the rigid body degrees of freedom is l. From the perspective of constructing the stiffness matrix, the necessary condition to ensure that the stiffness matrix is non-singular is m≥n-l, where m is the total number of independent stress coefficients. In this embodiment, according to this condition, considering the interaction term Φ rec When the orders are the same, the polynomial term Φ poly The influence brought by different orders.

[0146] It should be noted that as the orders of the stress function and the electric displacement function Φ poly of the inclusion are selected differently, it is difficult to meet the condition of m≥n-l, resulting in the singularity of the stiffness matrix. Therefore, when discussing the order of the polynomial term number Φ poly in this embodiment, it is necessary to ensure that the order of the interaction term number Φ poly meets the necessary condition for the matrix to be non-singular.

[0147] Discuss the order of the polynomial term number Φ poly of the stress function and the electric displacement function of the matrix and the inclusion, and observe the convergence of VCFEM. Extract the stress on the horizontal center line as Figure 10 shown. It can be seen that when the order of the interaction term number Φ rec of the stress function and the electric displacement function is 3, 4, 5, 6 orders, the polynomial term number Φ polyNon - convergence occurs when the order is 8. When the number of interaction terms Φ of the stress function and the electric displacement function rec is 7 or 8, the calculation results converge. And when the number of polynomial terms Φ poly increases from 3 to 10, the calculation accuracy becomes more precise. Therefore, in this embodiment, by selecting an appropriate order, both calculation accuracy and calculation efficiency can be taken into account.

[0148] Simulation of multi - inclusion piezoelectric elements

[0149] On the premise that the accuracy of VCFEM has been verified previously, this example combines the actual simulation of a relatively complex piezoelectric particle - reinforced composite material. In the target model, 20 elliptical inclusions of arbitrary sizes and directions are randomly distributed on a 1×1 mm² square matrix, and the semi - major axis lengths of the inclusions are about 0.024 mm - 0.033 mm. In this example, the VCFEM model contains 20 simple Voronoi elements, while in ABAQUS, the same model is divided into 37406 elements, as Figure 11 shown. The material parameters of the matrix and the inclusions are shown in Table 2.

[0150] Table 2 Material parameters of particle - reinforced composite piezoelectric materials

[0151]

[0152]

[0153] Set the critical normal stress value on the interface to 2 N. Apply a displacement load of 0.00001 m to each increment step of the right boundary of this model, and fix the other boundaries.

[0154] After Figure 12 comparison, it effectively proves the accuracy of 20 inclusion elements. Thus, this method demonstrates its effectiveness in simulating complex piezoelectric composite material structures, has high calculation accuracy, and can greatly simplify the calculation cost and improve the calculation efficiency compared with the finite - element software ABAQUS.

[0155] The present invention proposes a force - electricity coupled VCFEM with perfect bonding at the interface of piezoelectric particle composites and a force - electricity coupled VCFEM with interface debonding. By comparison with the finite element software ABAQUS, the accuracy and effectiveness of the VCFEM are verified. Compared with the ordinary displacement finite element, the number of elements of the VCFEM is less, which is particularly obvious in large - scale calculations. Through examples, it can be found in the present invention that even if only one element is used in the VCFEM, the stress concentration phenomenon at the matrix and the interface under the force - electricity coupling effect can be well captured. This shows that the Voronoi element has higher precision and faster calculation efficiency in micro - structure analysis. By introducing the electric potential field and the electric displacement field into the Voronoi element, the loading of the piezoelectric composite material is accurately simulated. The present invention lays a foundation for the VCFEM in future multi - physical - field coupling numerical calculations.

[0156] In addition, through the comparison of different orders of Φ poly and Φ rec in the present invention, it is found that in the VCFEM calculation, only when Φ poly and Φ rec reach a certain order, the calculation results will converge. The order of Φ poly and Φ rec will affect the accuracy of the VCFEM calculation results, and the order size is positively correlated with the accuracy. According to this, different orders can be selected in the present invention to ensure a balance between accuracy and calculation efficiency and achieve the most efficient numerical simulation.

[0157] 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 embodiments described. Obviously, many modifications and variations can be made according to the content of this specification. These embodiments are selected and specifically described in this specification to better explain the principle and practical application of the present invention, so that those skilled in the art in the relevant technical field can understand and utilize the present invention well. The present invention is only limited by the claims and their full scope and equivalents.

Claims

1. A modeling method for an interface bonding unit model of a piezoelectric particle composite material, characterized in that: The following steps are involved: S1: Construct a hybrid energy functional, describe it by generalized variational principle, deal with six independent variables: stress, strain, displacement, electric displacement, electric field intensity and electric potential. By introducing the hybrid complementary energy functional, eliminate the strain and electric field intensity in the functional and obtain a simplified equation; S2: Select Voronoi grid for domain partitioning, divide the material into multiple Voronoi cells, and make each cell boundary consist of the center bisector of the inclusion and the adjacent inclusion; S3: The Lagrange multiplier method is used to introduce interface conditions to ensure the balance between stress and electric displacement; four constraints are added to handle rigid body motion to ensure the non-singularity and numerical stability of the column matrix; S4: The accuracy and reliability of the model in practical applications are verified by comparing with the commercial finite element software ABAQUS. In different examples, single inclusion and multi-inclusion cases are tested to ensure that the model can correctly simulate the multi-field coupling behavior of piezoelectric composites.

2. The modeling method of the piezoelectric particle composite material interface bonding unit model according to claim 1, characterized in that: The step S1 specifically includes: using Respectively represent S u The prescribed displacement and S at the displacement boundary t The prescribed traction at the external force boundary; express The prescribed electromotive force and S at the potential boundary ω The prescribed surface charge density at the electric density boundary; in the finite element analysis method, both mechanical and electrical boundary conditions can satisfy displacement and potential as node variables; through two constraints and can be simplified to The stress σ, displacement u, electric displacement D, and electric potential φ field function variables are set, and the piezoelectric constitutive relation σ=C E γ-e T E、D=eγ-ε γ E, eliminating the strain γ and the electric field strength E in the functional, we get a simplified equation, denoted as Π σD : Assuming that the piezoelectric material is two-dimensionally isotropic, under plane stress state, its constitutive relation equation can be obtained: Where c, e, and d are the elastic constant, piezoelectric coefficient, and dielectric constant, respectively.

3. The modeling method of the piezoelectric particle composite material interface bonding unit model according to claim 1, characterized in that: The step S2 specifically includes: discretizing the domain Ω into subdomains Ω e , so Ω=Σ e Ω e ; At the same time, each Ω e There is an inclusion Ω c , which satisfies Using Ω m Represents the matrix material, i.e. Ω e =Ω m +Ω c ; In each subdomain Ω e The complete function σ of the internal independence assumption i and D i ; After discretization, each Voronoi unit domain has nodes on the outer and inner boundaries. D c Respectively and field in Φ m , Φ c Respectively represent the boundaries and The displacement and electric potential on .

4. The modeling method of the piezoelectric particle composite material interface bonding unit model according to claim 1, characterized in that: The step S3 specifically includes: pre-satisfying the basic boundary conditions and the constraints between the units by applying the constraints implemented by the Lagrange multiplier. Interface between matrix and inclusion The force continuity and electric displacement continuity conditions; The following conditions should be met in the electromechanical coupling area and on the boundary: At a given traction boundary S t Previous (5) At a given electric density boundary S w Previous (6) At the border between cells Previous (7) At the border between cells Previous (8) At the delamination interface Previous (9) At the delamination interface Previous (10) At the delamination interface Previous (11) At the delamination interface Previous (12) At the undelaminated interface Previous (13) At the undelaminated interface Previous (14) Based on the above derived hybrid complementary energy functional formula, the Lagrange multiplier method is applied to implement the above constraints and the modified hybrid complementary energy functional is obtained: in, Boundaries meet: The displacement and potential of each boundary are Boundary displacement / potential d is the interpolation of generalized nodal displacement / potential q d=Lq (18) The Voronoi unit introduces an assumed independent stress field / electric displacement field, which should satisfy the equilibrium condition, and the equilibrium condition can be obtained through the stress function / electric displacement function; the stress functions of the matrix and inclusion are in, and is the Airy stress function of the matrix and the inclusions; considering the existence of the electric displacement field, and It is defined as the electric displacement function corresponding to the stress function; and are the interaction functions of stress and electric displacement, respectively, which consider the complete polynomial based on the inverse of the shape function at the matrix-inclusion interface; by differentiating the stress function and the electric displacement function, the stress and electric displacement at each location of the unit can be obtained: in, is the stress matrix of matrix and inclusion, is the electric displacement matrix of the matrix and inclusions, are the stress parameters of matrix and inclusion, is the electric displacement parameter of matrix and inclusion; After finite element discretization, we can get in By varying equation (24), according to the stationary value principle, we can obtain Can get It can be simplified to β=H -1 Gq(34) Substituting equation (34) into equation (24), according to the stationary value principle, we have The system of equations for solving the generalized displacement / potential can be obtained Here, the element stiffness matrix is ​​equal to: The rigid body displacement / potential for the unit boundary nodes can be expressed as Among them, x i ,y i , is the coordinate of node i, α1, α2 are the rigid body translation moments and α3 represents the corresponding potential constant, α4 is the rigid body rotational freedom; interface rigid body node displacement / potential It can be represented by the matrix φ'; then the node displacement / potential at the interface and unit boundary can be expressed as dq'=φ'α+dq' def (38) dq=φα+dq def (39) where dq' def and dq def are the pure deformation displacement / potential of the interface and unit boundary nodes respectively; because the pure deformation mode and rigid body mode space are orthogonal, both sides of equations (38) and (39) are multiplied by φ' T and φ T have: f' T dq'=φ' T f'a (40) f T dq=φ T fa (41) During the iterative solution process, for the same unit, the unit boundary node displacement / potential iteration increment is dq ei The nodal displacement / potential iteration increment on the matrix side of the matrix-inclusion interface is equal to The rigid body displacements of the two parts are the same, so: where Φ = [{φ T φ} -1 φ T -{φ' T φ'} -1 φ' T ​ Formula (43) provides four additional constraints to eliminate the rigid body displacement / potential in the unit, ensuring that the unit stiffness matrix K e It is not singular. The constraints of the above formula are introduced into the modified complementary energy functional through the Lagrange multiplier method to obtain According to the stationary value principle of the co-energy functional, get Also according to get Combining equations (45) and (46), we can obtain the solution equation The stiffness matrix is ​​equal to

Images