A composite failure simulation method based on extended finite element

Abstract

The application belongs to the technical field of aircraft strength test, and particularly relates to a composite material failure simulation method based on extended finite element. The method comprises the following steps: S1, obtaining stress values of a composite material in a previous increment step during execution of a specified load loading process; S2, determining a fiber failure value and an inter-fiber failure value; S3, when the fiber failure value is greater than 1, calculating a first reduction multiple of material properties according to the fiber failure value, and when the inter-fiber failure value is greater than 1, calculating a second reduction multiple of material properties according to the inter-fiber failure value; S4, when the inter-fiber failure value is greater than 1, further determining a unit normal vector of a potential fracture surface in a global coordinate system for crack propagation; and S5, calculating stress values of a next increment step according to the reduced material properties, and returning to step S1. The application can accurately and efficiently simulate the failure behavior of the composite material.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of aircraft strength test, and particularly relates to a composite material failure simulation method based on extended finite element. BACKGROUND

[0002] Composite materials are widely used in the field of aerospace due to their excellent mechanical properties. The composite material strength theory and failure analysis method has been a hot and difficult point of academic research. Over the past half century, researchers at home and abroad have proposed dozens of theories and analysis methods, but so far, no method can successfully predict all the observed composite material failure behaviors. Due to the weak interlaminar performance, composite materials are prone to delamination damage, and how to introduce initial delamination damage and simulate the composite material failure analysis when the delamination expands is another difficulty.

[0003] In the aspect of composite material failure analysis, relevant researchers have carried out a lot of work. However, in the prior art, although crack prediction is carried out through the matrix failure criterion, the simulation method of matrix failure mostly depends on mesh division, and the calculation amount is large. Or the failure criterion used does not consider the fiber fracture, and the simulation effect is not accurate. SUMMARY

[0004] In order to solve the above technical problems, the application provides a composite material failure simulation method based on extended finite element. The maximum stress criterion is used to determine the fiber failure, and the stiffness of the composite material after fiber failure is exponentially reduced. The inter-fiber failure is determined by using the fracture surface failure criterion, and the matrix crack propagation direction is determined. The inter-fiber failure of the matrix crack propagation behavior is simulated based on the extended finite element method, so as to accurately simulate the progressive failure behavior of the composite material.

[0005] The composite material failure simulation method based on extended finite element provided by the application mainly includes:

[0006] Step S1, obtaining a stress value of a composite material in a specified load loading process in a previous increment step;

[0007] Step S2, determining a fiber failure value F FF At the same time, the potential fracture surface angle of the composite material is determined based on the golden section search method, and the inter-fiber failure value F IFF is determined according to the fracture surface angle.

[0008] Step S3, when the fiber failure value F FF is greater than 1, the material property first reduction multiple d FF is calculated according to the fiber failure value F FF , and when the inter-fiber failure value F IFFgreater than 1, according to the inter-fiber failure value F IFF calculating a second reduction factor d of the material property IFF ;

[0009] Step S4, when the inter-fiber failure value F IFF greater than 1, further determining the unit normal vector n0 of the potential fracture surface under the global coordinate system for crack propagation;

[0010] Step S5, calculating the stress value of the next increment step according to the reduced material property, and returning to step S1.

[0011] Preferably, in step S2, the fiber failure value F FF includes:

[0012]

[0013] wherein X T is the fiber tensile strength, X C is the fiber compressive strength, and σ1 is the fiber direction stress.

[0014] Preferably, in step S2, the inter-fiber failure value F IFF includes:

[0015]

[0016] wherein σ n , τ nl , τ nt are the stresses in each direction on the potential fracture surface, and A1, A2, B2, C2 are coefficients.

[0017] Preferably, the coefficients A1, A2, B2, C2 are determined by the following method:

[0018]

[0019] wherein Y T is the transverse tensile strength, Y C is the transverse compressive strength, S L is the in-plane shear strength, and θ0 is the fracture angle under pure transverse compression.

[0020] Preferably, in step S3, the first reduction factor d of the material property FF includes:

[0021]

[0022] wherein L C is the element characteristic length, G FT , G FC are the longitudinal tensile and compressive fracture toughnesses, respectively.

[0023] Preferably, in step S3, the second reduction factor d IFF comprises:

[0024]

[0025] where L C is the characteristic length of the element, G MT , G MC are the transverse tensile and compressive fracture toughnesses.

[0026] Preferably, in step S4, determining the unit normal vector of the potential fracture surface in the global coordinate system comprises:

[0027] Step S41, determining the unit normal vector n of the potential fracture surface in the material principal axis coordinate system according to the angle θ of the potential fracture surface:

[0028]

[0029] Step S42, determining the unit normal vector n0 in the global coordinate system according to the rotation angle of the material principal axis coordinate system relative to the global coordinate system:

[0030]

[0031] The present application takes into account the progressive failure of fibers and inter-fiber matrix, and can accurately and efficiently simulate the failure behavior of composite materials. BRIEF DESCRIPTION OF DRAWINGS

[0032] Figure 1 is a flowchart of a preferred embodiment of the composite material failure simulation method of the present application based on extended finite elements.

[0033] Figure 2 is a schematic diagram of stress components on a potential fracture surface.

[0034] Figure 3 is a schematic diagram of specimen size. DETAILED DESCRIPTION

[0035] For the purpose, technical solutions and advantages of the embodiments of the present application to be clearer, the technical solutions in the embodiments of the present application will be described in more detail below with reference to the drawings in the embodiments of the present application. In the drawings, the same or similar notations represent the same or similar elements or elements with the same or similar functions throughout. The described embodiments are part of the embodiments of the present application, rather than all the embodiments. The embodiments described below with reference to the drawings are exemplary and are intended to explain the present application, and cannot be understood as limiting the present application. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative labor fall within the scope of protection of the present application. The embodiments of the present application will be described in detail below with reference to the drawings.

[0036] The present application provides a composite material failure simulation method based on extended finite element, as shown in the figure, mainly comprising: Figure 1

[0037] Step S1, obtaining the stress value of the composite material in the execution of the specified load loading process in the last increment step;

[0038] Step S2, determining the fiber failure value F FF based on the stress value, and determining the inter-fiber failure value F IFF based on the potential fracture surface angle of the composite material determined based on the golden search method.

[0039] Step S3, when the fiber failure value F FF is greater than 1, calculating the first reduction factor d FF of the material properties according to the fiber failure value F FF , and when the inter-fiber failure value F IFF is greater than 1, calculating the second reduction factor d IFF of the material properties according to the inter-fiber failure value F IFF .

[0040] Step S4, when the inter-fiber failure value F IFF is greater than 1, further determining the unit normal vector n0 of the potential fracture surface under the global coordinate system for crack propagation.

[0041] Step S5, calculating the stress value of the next increment step according to the reduced material properties, and returning to step S1.

[0042] The present application adopts the maximum stress criterion to determine the fiber failure, and performs exponential reduction on the stiffness of the composite material after the fiber failure. The fiber inter-failure is determined by the fracture surface failure criterion, and the matrix crack propagation direction is determined. The matrix crack propagation behavior after the fiber inter-failure is simulated based on the extended finite element method. ​

[0043] The inter-fiber matrix failure / delamination failure can be regarded as matrix crack, which is simulated by the extended finite element method based on the cohesive force; the fiber failure can be regarded as brittle fracture, and the material properties of the failed element can be reduced by field variables. The overall process includes: in step S1, at the beginning of each time increment step, the stress state of the previous time step is obtained, then in step S2, the fiber failure function and the inter-fiber failure function are determined respectively, then in step S3, the fiber failure value and the inter-fiber failure value are calculated respectively, and whether the material properties are reduced is determined based on the failure value, after the reduction, the material properties change, at the same time, the stress of the composite material changes continuously under the current load loading state, the stress value is recalculated in step S4, and then the cycle is realized by returning to step S1 until the load or displacement reaches the target value, or the composite material reaches the failure condition.

[0044] The present application is based on the ABAQUS platform, and the numerical method for simulating the fiber fracture and matrix crack of the composite material is realized by using the USDFLD and UDMGINI subprograms. As shown in FIG. 1, the USDFLD subprogram first calculates the fiber failure value and the inter-fiber failure value, and the fracture surface angle, and then calls the UDMGINI subprogram to define the crack initiation and crack direction. Figure 1

[0045] In some optional embodiments, in step S2, the fiber failure value F FF is determined by the following formula:

[0046]

[0047] Wherein, X T is the fiber tensile strength, X C is the fiber compression strength, and σ1 is the fiber direction stress.

[0048] F FF is the fiber failure function, and F FF > 1 represents fiber failure. After the fiber failure, the material parameters longitudinal modulus E 11 , transverse modulus E 22 , out-of-plane modulus E 33 , longitudinal shear modulus G 12 , and transverse shear modulus G 23 are reduced to (1-d FF ) times of the undamaged material value. In some optional embodiments, in step S3, the first reduction multiple d FF of the material properties is calculated by the following formula:

[0049]

[0050] Wherein, L C is the element characteristic length, G FT , and G​FC These are longitudinal tensile and compressive fracture toughness, respectively.

[0051] In some alternative implementations, in step S2, the inter-fiber failure value F is determined. IFF include:

[0052]

[0053] Where, σ n τ nl τ nt Let A1, A2, B2, and C2 represent the stresses in each direction on the potential fracture surface, and A1, A2, B2, and C2 be coefficients. It should be noted that the reference... Figure 2 , σ n τ nl τ nt The value of σ is related to the stress value given in step S1 and the fracture surface angle θ. Different θ values ​​correspond to different stresses σ. n τ nl τ nt This corresponds to different inter-fiber failure values ​​F. IFF Therefore, this application first determines the corresponding potential fracture surface angle θ using the golden search method, and then the interfiber failure function F can be determined. IFF The maximum value.

[0054] In some alternative implementations, the coefficients A1, A2, B2, and C2 are determined in the following manner:

[0055]

[0056] Among them, Y T Y represents the transverse tensile strength. C For transverse compressive strength, S L θ is the in-plane shear strength, and θ0 is the fracture angle under pure transverse compression, which is typically 53°.

[0057] After inter-fiber failure, all material parameters are reduced to (1-d) of the undamaged material value. IFF In some alternative embodiments, in step S3, the second reduction factor d of the material property is calculated. IFF include:

[0058]

[0059] Among them, L C G is the characteristic length of the unit cell. MT G MC It exhibits transverse tensile and compressive fracture toughness.

[0060] In some alternative implementations, step S4, determining the unit normal vector of the potential fracture surface in the global coordinate system, includes:

[0061] Step S41, determine the unit normal vector n of the potential fracture surface in the material principal axis coordinate system (such as 1-2-3 coordinate system) according to the potential fracture surface angle θ: Figure 2

[0062]

[0063] Step S42, determine the unit normal vector n0 in the global coordinate system according to the rotation angle of the material principal axis coordinate system relative to the global coordinate system

[0064]

[0065] In step S42, assuming the ply angle The material principal axis coordinate system (1-2-3) is rotated by the global coordinate system (x-y-z) around the z-axis Therefore, the rotation angle That is, the ply angle

[0066] After the inter-fiber matrix fails, the degradation of the fracture surface bonding stiffness is calculated by the built-in program of ABAQUS according to the damage evolution law using the traction-separation response. Considering the complex stress state of the crack tip, the BK criterion of mixed mode can be used to control the stiffness degradation.

[0067] Take the edge open compression virtual test containing impact damage as an example for specific description.

[0068] The specimen size and delamination damage are shown in Figure 3 The ply sequence is [0 / 45 / -45 / 0 / 90 / 0 / 45 / -45 / 0 / 45 / -45 / 0 / 90 / 45 / -45 / 0 / 0 / 45 / -45 / 0]S, a total of 40 layers, and the single-layer thickness is 0.191 mm. The material performance parameters are as follows: longitudinal modulus 166 GPa, transverse modulus 8110 MPa, shear modulus 4140 MPa, in-plane Poisson's ratio 0.319, longitudinal tensile strength 2920 MPa, longitudinal compressive strength 1140 MPa, transverse tensile strength 70.6 MPa, transverse compressive strength 278 MPa, and shear strength 172 MPa.

[0069] The implementation process is as follows:

[0070] 1) According to the provided material performance parameters, determine the USDFLD and UDMGINI subprograms.

[0071] ​​2) Build the geometric model and material model of the test piece: the test piece is divided into 3 layers, the middle layer contains 2 0° layers, and the initial delamination damage is simulated by embedding a circular plane in the middle layer; the composite solid element C3D8R is modeled, and the number of model grids is 166860. The material model is set as anisotropic solid, the field variable is set as 1, and the number of state variables is 3.

[0072] 3) Set the load boundary condition, set the left and right clamping areas as fixed support boundary conditions, and apply a -x axis direction 1.5mm compression displacement to the right side load edge to perform failure analysis.

[0073] 4) During the failure process, the delamination expansion and fiber fracture can be clearly seen, and the final failure load is 221kN, which verifies the feasibility of the method.

[0074] The above is only a specific embodiment of the present application, but the protection scope of the present application is not limited to this. Any person skilled in the art can easily think of changes or replacements within the technical scope disclosed in the present application, which should be covered within the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. An extended finite element-based composite failure simulation method, characterized by, Comprise: Step S1, obtaining the stress value of the composite material in the last incremental step in the execution of the specified load loading process; Step S2, determining a fiber failure value F according to the stress value FF Meanwhile, the potential fracture surface angle of the composite material is determined based on the golden search method, and the inter-fiber failure value F is determined according to the fracture surface angle IFF ; Step S3, when the fiber failure value F FF greater than 1, the fiber failure value F FF calculating a first reduction factor d FF when the inter-fiber failure value F IFF greater than 1, the inter-fiber failure value F IFF calculating a second reduction factor d IFF ; Step S4, when the inter-fiber failure value F IFF when greater than 1, further determining a unit normal vector n0of the potential fracture surface under the global coordinate system that the crack propagates; Step S5, calculating the stress value of the next incremental step according to the reduced material properties, and returning to step S1; wherein in step S2, the fiber failure value F is determined FF comprising: where X T is the fiber tensile strength, X C is the fiber compressive strength, and σ1is the fiber direction stress; In step S2, the inter-fiber failure value F is determined IFF comprising: where σ n , τ nl , τ nt are stresses in each direction on the potential fracture surface, and A1, A2, B2, C2 are coefficients. The coefficients A1, A2, B2 and C2 are determined by the following method: where Y T is the transverse tensile strength, Y C is the transverse compressive strength, S L is the in-plane shear strength, and θ0is the fracture angle for pure transverse compression. In step S3, the first reduction factor d of the material property is calculated FF comprising: wherein L C is the unit characteristic length, G FT , G FC are the longitudinal tensile and compressive fracture toughness, respectively, and E 11 is the longitudinal modulus of the material parameters; In step S3, the second reduction factor d of the material property is calculated IFF comprising: where L C is the unit characteristic length, G MT , G MC is the transverse tensile, compressive fracture toughness, and E2 is the transverse modulus of the material parameter.

2. The extended finite element based composite failure simulation method of claim 1, wherein, In step S4, determining the unit normal vector of the potential fracture surface in the global coordinate system comprises: Step S41, determining the unit normal vector n of the potential fracture surface in the material principal axis coordinate system according to the potential fracture surface angle θ: Step S42, determining the rotation angle of the material principal axis coordinate system relative to the global coordinate system determining a unit normal vector n0 in the global coordinate system:

Images