A damage analysis method for multi-layer composite materials based on numerical simulation calculations

Abstract

The present invention discloses a method for damage analysis of multi-layer composite materials based on numerical simulation calculations, which relates to the technical field of composite material damage analysis and includes the following steps: S1: Construct the initial single-layer stiffness matrix C' of each layer; S2: Construct the single-layer stiffness matrix A layer and the single-layer mass matrix M layer ; S3: Set free boundary conditions; S4: Construct the global stiffness matrix A global and the global mass matrix M global ; S5: Perform regularization processing on the global stiffness matrix A global to obtain #imgabs0#; S6: Construct a dispersion equation, solve the dispersion equation to obtain eigenvalues λ and eigenvectors Ф, and calculate the phase velocity v of Lamb waves p and the group velocity v g . The method of the present invention is used for damage analysis of orthotropic multi-layer composite materials and improves the accuracy of calculations.

Description

Technical Field

[0001] The present invention relates to the technical field of composite material damage analysis, and particularly to a method for analyzing damage of multi-layer composite materials based on numerical simulation calculation. Background Art

[0002] As a non-destructive testing method, Lamb wave is widely used in the field of health monitoring of composite materials due to its sensitivity to internal damage of materials. For example, non-destructive monitoring is carried out through Lamb wave, and damage analysis of composite materials is carried out through numerical simulation calculation.

[0003] Analyzing the damage of composite materials through numerical simulation calculation means that first, a model of the composite material is established, and then a wave number k is generated through numerical simulation (instead of actually applying a Lamb wave to the composite material), and on this basis, the propagation characteristic parameters of the Lamb wave in the composite material are calculated for the damage analysis of the composite material.

[0004] Currently, the damage analysis of composite materials through numerical simulation calculation mainly focuses on single-layer materials. However, for multi-layer composite materials, the case of interlayer coupling effect is often ignored. In multi-layer composite materials, the interlayer coupling effect has a significant impact on the propagation of Lamb wave. Ignoring this point may lead to errors.

[0005] In view of this, this invention is specifically proposed. Summary of the Invention

[0006] The purpose of the present invention is to solve the disadvantages existing in the prior art, and a method for analyzing damage of multi-layer composite materials based on numerical simulation calculation is proposed, which is used for the damage analysis of orthotropic multi-layer composite materials. The propagation characteristic parameters of Lamb wave in multi-layer composite materials are calculated through numerical simulation calculation, and the calculation accuracy is improved.

[0007] In order to achieve the above purpose, the present invention adopts the following technical scheme:

[0008] A method for analyzing damage of multi-layer composite materials based on numerical simulation calculation includes the following steps:

[0009] S1: Construct the initial single-layer stiffness matrix C' of each layer;

[0010] S2: According to the initial single-layer stiffness matrix C' of each layer, construct the single-layer stiffness matrix A layer and the single-layer mass matrix M layer ;

[0011] S3: Set free boundary conditions;

[0012] S4: Cumulatively construct layer by layer to obtain the global stiffness matrix A globaland the global mass matrix M global During the layer-by-layer cumulative construction process, perform continuous processing calculations between layers according to the interlayer coupling effect;

[0013] S5: Regularize the global stiffness matrix A global to obtain ;

[0014] S6: Based on and M global construct a dispersion equation, solve the dispersion equation to obtain eigenvalues λ and eigenvectors Ф, and calculate the group velocity v of Lamb waves p and the phase velocity v g .

[0015] Furthermore, S1 includes the following steps:

[0016] S11: Construct the original single-layer stiffness matrix C, in the form of:

[0017] ,

[0018] where

[0019] ,

[0020] ,

[0021] ,

[0022] ,

[0023] ,

[0024] ,

[0025] ,

[0026] ,

[0027] ,

[0028] where , , , are the Poisson ratios in the x-y, y-x, x-z, and y-z planes respectively, , , are the elastic moduli in the x, y, and z directions respectively, , , Represent the shear moduli in the x-y, y-z, and x-z planes, respectively;

[0029] S12: Construct a rotation matrix T in the following form:

[0030] ,

[0031] where,

[0032] ,

[0033] ,

[0034] where θ is the rotation angle;

[0035] S13: Calculate the initial single-layer stiffness matrix C' according to the following formula

[0036] .

[0037] Furthermore, S2 includes the following steps:

[0038] S21: Use the Gauss-Legendre integration method to calculate the elements of the single-layer stiffness matrix A layer and the single-layer mass matrix M layer for each layer through the following formula:

[0039] ,

[0040] ,

[0041] where N is the number of integration points, i, j = 0, 1,..., N-1, a and b are the indices of the Gauss-Legendre integration points, is the Legendre polynomial basis function of order N-1, is the weight, h is the layer thickness, is the density.

[0042] Furthermore, S3 includes the following steps:

[0043] S31: For the top layer of the multi-layer composite material, set the free boundary conditions through the following formula:

[0044] ,

[0045] ,

[0046] S32: For the bottom layer of the multi-layer composite material, set the free boundary conditions through the following formula:

[0047] ,

[0048] ,

[0049] Among them, M is the number of layers of the multi-layer composite material.

[0050] Furthermore, the said S4 includes the following steps:

[0051] S41: Cumulatively construct the global stiffness matrix of the first layer according to the following formula and the global mass matrix ,

[0052] ,

[0053] ,

[0054] S42: Cumulatively construct the global stiffness matrix of the next layer according to the following formula and the global mass matrix ,

[0055] ,

[0056] ,

[0057] where c = 1, 2,... M, and M is the number of layers of the multi-layer composite material;

[0058] S43: Update the elements of the global stiffness matrix upper and the global mass matrix lower by superimposing the adjacent layer matrix blocks A and to satisfy the displacement and stress continuity conditions;

[0059] S44: Loop and execute S42 and S43 until the cumulative construction is completed to obtain the global stiffness matrix A global and the global mass matrix M global .

[0060] Furthermore, the said S43 includes the following steps:

[0061] S431: Construct the matrix block A upper at the bottom of the upper layer, and calculate the elements in A upper according to the following formula:

[0062] ,

[0063] where F is the number of basis functions used in discretization;

[0064] S432: Construct the matrix block A lower at the top of the current layer, and calculate the elements in A lower according to the following formula:

[0065] ,

[0066] S433: Update the elements in according to the following formula:

[0067] .

[0068] Further, the S5 includes the following steps:

[0069] S51: Calculate and obtain ,

[0070] ,

[0071] where = 10 -6 , I is the identity matrix of the same order as A global , and is the condition number of A global .

[0072] Further, the S6 includes the following steps:

[0073] S61: Construct a dispersion equation in the following form:

[0074] ,

[0075] S62: Solve the dispersion equation to obtain a series of eigenvalues λ = , and the corresponding eigenvectors Ф, where ω is the angular frequency;

[0076] S63: Delete the imaginary eigenvalues λ to obtain the effective eigenvalues λ and the effective eigenvectors Ф;

[0077] S64: Calculate the phase velocity v of the Lamb wave according to the following formula p and the phase velocity v g，

[0078] ,

[0079] ,

[0080] where k is the wave number of the multi-layer composite material and f is the frequency.

[0081] Further, it also includes the following steps:

[0082] S7: Construct a visualization image with the frequency f as the x-axis, the phase velocity v p or the phase velocity v g as the y-axis, and the effective eigenvector Ф as the z-axis.

[0083] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0084] It can calculate the phase velocity and group velocity of Lamb waves in different frequency ranges of orthotropic multi-layer composite materials; by accurately modeling the global stiffness matrix and global mass matrix of the multi-layer composite material, and performing interlayer continuity processing calculations according to the interlayer coupling effect, it can accurately reflect the stress-strain relationship of each layer. Compared with the prior art, the accuracy of Lamb wave propagation characteristic analysis is improved, making the damage analysis and research of multi-layer composite materials more accurate and reliable, and providing accurate theoretical support for the damage detection of multi-layer composite materials. BRIEF DESCRIPTION OF THE DRAWINGS

[0085] Figure 1 is a flowchart of a damage analysis method for multi-layer composite materials based on numerical simulation calculations;

[0086] Figure 2 is a schematic diagram of a visualization image. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0087] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to 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.

[0088] Embodiment 1:

[0089] A damage analysis method for multi-layer composite materials based on numerical simulation calculations, as Figure 1 shown, includes the following steps:

[0090] S1: Construct the initial single-layer stiffness matrix C' of each layer,

[0091] S2: Construct the single-layer stiffness matrix A layer and the single-layer mass matrix M layer，

[0092] S3: Set the free boundary conditions,

[0093] S4: Cumulatively construct layer by layer to obtain the global stiffness matrix A global and the global mass matrix M global , and during the cumulative construction process layer by layer, perform interlayer continuity processing calculations according to the interlayer coupling effect,

[0094] S5: Perform regularization processing on the global stiffness matrix A global to obtain ,

[0095] S6: Based on and Mglobal Construct a dispersion equation, solve the dispersion equation to obtain the eigenvalue λ and the eigenvector Ф, and calculate the group velocity v of the Lamb wave p and the phase velocity v g .

[0096] The method for analyzing damage of multi-layer composite materials based on numerical simulation calculation in this embodiment can calculate the phase velocity and group velocity of Lamb waves in different frequency ranges of orthotropic multi-layer composite materials; by accurately modeling the global stiffness matrix and global mass matrix of the multi-layer composite material, and performing interlayer continuity processing calculation according to the interlayer coupling effect, the stress-strain relationship of each layer can be accurately reflected. Compared with the prior art, the accuracy of the analysis of Lamb wave propagation characteristics is improved, making the damage analysis and research of multi-layer composite materials more accurate and reliable, and providing accurate theoretical support for the damage detection of multi-layer composite materials.

[0097] In an optional embodiment, S1 includes the following steps:

[0098] S11: Construct an original single-layer stiffness matrix C, in the following form:

[0099] ,

[0100] where

[0101] ,

[0102] ,

[0103] ,

[0104] ,

[0105] ,

[0106] ,

[0107] ,

[0108] ,

[0109] ,

[0110] where , , , are the Poisson's ratios in the x-y, y-x, x-z, and y-z planes respectively, , , The elastic moduli in the x, y, and z directions respectively, which respectively reflect the ability of the material to resist elastic deformation in the x, y, and z directions, and and respectively represent the shear moduli in the x-y, y-z, and x-z planes, which respectively reflect the ability of the material to resist shear deformation in the x-y, y-z, and x-z planes; in addition, the Poisson's ratio, elastic modulus, and shear modulus are material properties, where can be obtained by calculation through the following formula,

[0111] .

[0112] S12: Construct a rotation matrix T in the following form:

[0113] ,

[0114] where

[0115] ,

[0116] ,

[0117] where θ is the rotation angle, which is used to describe the rotation degree of the material ply relative to the reference coordinate system.

[0118] S13: Calculate the initial single-layer stiffness matrix C' according to the following formula

[0119] .

[0120] In this alternative embodiment, the original single-layer stiffness matrix C is constructed based on the material orthotropic property, and the initial single-layer stiffness matrix C' is constructed by combining the rotation angle θ of each layer for rotation transformation, thereby ensuring that the matrix C' can accurately reflect the stress-strain relationship of each layer.

[0121] In an alternative embodiment, S2 includes the following steps:

[0122] S21: Use the Gaussian-Legendre integration method to calculate each element in the single-layer stiffness matrix A layer and the single-layer mass matrix M layer by the following formula:

[0123] ,

[0124] ,

[0125] where N is the number of integration points, i, j = 0, 1,..., N - 1, a and b are the index values of the Gaussian-Legendre integration points, is the Legendre polynomial basis function of order N - 1, is the weight, h is the layer thickness, is the density; in addition, the density is a material property.

[0126] Preferably, according to the convergence analysis, the number of integration points N can be taken as 5, and at this time the relative error is less than 1%.

[0127] In an alternative embodiment, the S3 includes the following steps:

[0128] S31: For the top layer of the multi-layer composite material, set the free boundary conditions through the following formula:

[0129] ,

[0130] .

[0131] S32: For the bottom layer of the multi-layer composite material, set the free boundary conditions through the following formula:

[0132] ,

[0133] ,

[0134] where M is the number of layers of the multi-layer composite material.

[0135] In an alternative embodiment, the S4 includes the following steps:

[0136] S41: Cumulatively construct the global stiffness matrix of the first layer according to the following formula and the global mass matrix ,

[0137] ,

[0138] .

[0139] S42: Cumulatively construct the global stiffness matrix of the next layer according to the following formula and the global mass matrix ,

[0140] ,

[0141] ,

[0142] where c = 1, 2,... M, and M is the number of layers of the multi-layer composite material.

[0143] S43: Update the global stiffness matrix by superimposing the adjacent layer matrix blocks A upper and A lower ​ and the global mass matrix whose elements satisfy the displacement and stress continuity conditions.

[0144] S44: Loop through S42 and S43 until the cumulative construction is completed, i.e., the cumulative calculation reaches and to obtain the global stiffness matrix A global and the global mass matrix M global .

[0145] In this alternative embodiment, first calculate and to calculate and , and then perform layer-by-layer cumulative calculation until the cumulative calculation of and is completed.

[0146] In an alternative embodiment, S43 includes the following steps:

[0147] S431: Construct the matrix block A upper at the bottom of the upper layer, and calculate the elements in A upper according to the following formula:

[0148] ,

[0149] where F is the number of basis functions used in discretization.

[0150] S432: Construct the matrix block A lower at the top of the current layer, and calculate the elements in A lower according to the following formula:

[0151] .

[0152] S433: Update the elements in according to the following formula:

[0153] .

[0154] Starting from the second layer of the multi-layer composite material, it is necessary to ensure the stress continuity between adjacent layers. In this alternative embodiment, by performing interlayer continuity processing calculation based on the interlayer coupling effect, the stress continuity between adjacent layers can be ensured.

[0155] Specifically, through the method of matrix block diagonal splicing, combined with interlayer continuity processing, the global stiffness matrix A global and the global mass matrix M globalThe accuracy and integrity are achieved, enabling the modeling of the overall characteristics of multi-layer composite materials and providing a solid foundation for solving the dispersion equation.

[0156] In an alternative embodiment, S5 includes the following steps:

[0157] S51: Calculate and obtain through the following formula ,

[0158] ,

[0159] where = 10 -6 , I is the identity matrix of the same order as A global ; in addition, it should be noted that when 10 is the condition number of A global ; in addition, it should be noted that when 10 −8 ≤ ≤ 10 −4 , the deviation of the phase velocity calculation result is less than 0.5%.

[0160] In an alternative embodiment, S6 includes the following steps:

[0161] S61: Construct the dispersion equation in the following form:

[0162] ,

[0163] S62: Solve the dispersion equation to obtain a series of eigenvalues λ = , and the corresponding eigenvectors Ф for each, where ω is the angular frequency.

[0164] In this alternative embodiment, the dispersion equation can be solved by numerical methods (such as the eigh function).

[0165] Different angular frequencies ω correspond to different propagation modes of Lamb waves in multi-layer composite materials, such as the symmetric mode S n and the anti-symmetric mode A n ; the eigenvector Ф describes the vibration mode of Lamb waves at the corresponding angular frequency, and the eigenvector Ф includes the displacement components u 1 , u 2 and u 3 , and can be expressed as follows:

[0166] ,

[0167] The formula for calculating its amplitude is as follows:

[0168] .

[0169] S63: Delete the imaginary eigenvalue λ to obtain the effective eigenvalue λ and the effective eigenvector Ф.

[0170] Since the imaginary eigenvalue λ corresponds to an evanescent wave or a non-physical solution, in this optional embodiment, the imaginary eigenvalue λ is deleted, that is, the real eigenvalues λ are screened out from a series of eigenvalues λ obtained by calculation in S62, and then the effective eigenvalue λ and the effective eigenvector Ф are obtained to characterize the propagating Lamb wave mode.

[0171] S64: Calculate the phase velocity v of the Lamb wave according to the following formula p and the phase velocity v g，

[0172] ,

[0173] ,

[0174] where k is the wave number of the multi-layer composite material and f is the frequency.

[0175] For the wave number k, an initial guessed value of the wave number k 0 can be given first. This initial value can be selected according to the prior knowledge or experience of the physical problem, and then the Newton-Raphson method is used for iterative solution to obtain the wave number k. In addition, the angular frequency ω and the frequency f can be converted according to the following formula.

[0176] In this optional embodiment, by solving the generalized eigenvalue problem, combining eigenvalue deletion and screening, the dispersion equation can be accurately solved, and the phase velocity v p and the phase velocity v g of the Lamb wave can be calculated, so as to comprehensively describe the propagation characteristics of the Lamb wave in the multi-layer composite material.

[0177] In an optional embodiment, the following steps are further included:

[0178] S7: Construct a visualization image with the frequency f as the x-axis, the phase velocity v p or the phase velocity v g as the y-axis, and the amplitude of the effective eigenvector Ф as the z-axis.

[0179] In this optional embodiment, by constructing a visualization image with the frequency f as the x-axis, the phase velocity v p or the phase velocity v g as the y-axis, and the amplitude of the effective eigenvector Ф as the z-axis, the changing trends of the phase velocity and the group velocity with the change of the frequency can be clearly presented, making the analysis results of the propagation characteristics of the Lamb wave more intuitive and easy to understand, facilitating technicians to observe and analyze, and improving the efficiency and accuracy of damage analysis and research.

[0180] Specifically, the phase velocity refers to the propagation velocity of the phase front of the Lamb wave. Therefore, the curve constructed with the phase velocity can show the phase propagation velocity of the wave at different frequencies; the group velocity refers to the energy propagation velocity of the Lamb wave. Therefore, the dispersion curve constructed with the group velocity can show the energy propagation velocity of the wave at different frequencies; by observing the dispersion curve graph, the Lamb wave mode at different frequencies can be determined; the dispersion phenomenon can be identified through the bending and change of the dispersion curve, which is of great significance in both non-destructive testing and structural health monitoring.

[0181] Specifically,

[0182] To further demonstrate the calculation results of this embodiment, carbon fiber reinforced material (CFRP) T700M913 with a size of 500*500 mm, 8 layers, and the rotation angle of the upper four layers being [0 / 90] 2s , and the rotation angle of the lower four layers being [90 / 0] 2s , with each layer having a thickness of 0.15 mm, is taken as an example for illustrative calculation. The material properties of this material are shown in Table 1.

[0183] Table 1 Material properties of carbon fiber reinforced material

[0184] <![CDATA[ρ (kg / m 3 )]]> <![CDATA E x ( GPα )]]> <![CDATA E y ( GPα )]]> <![CDATA E z ( GPα )]]> <![CDATA G xy ( GPα )]]> 1.56 135 10 10 5 <![CDATA G xz ( GPα )]]> <![CDATA G yz ( GPα )]]> <![CDATA v xy > <![CDATA v xz > <![CDATA v yz > 5 5 0.3 0.25 0.25

[0185] The visualized image is as Figure 2 shown, where the color bar on the right represents the amplitude of the effective eigenvector Ф. In this embodiment, the ordinate is taken as the phase velocity for demonstration.

[0186] As mentioned above, the above are only the preferred specific embodiments of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention, according to the technical solution and inventive concept of the present invention, makes equivalent substitutions or changes, and all should be covered within the protection scope of the present invention.

Claims

1. A damage analysis method for multilayer composite materials based on numerical simulation calculation, characterized in that: The steps include: S1: Construct the initial single-layer stiffness matrix C' of each layer; S2: Construct the single-layer stiffness matrix A of each layer according to the initial single-layer stiffness matrix C' of each layer layer and the single-layer mass matrix M layer ; S3: Set free boundary conditions; S4: Accumulate layer by layer to obtain the global stiffness matrix A global and the global mass matrix M global , in the process of layer-by-layer accumulation and construction, inter-layer continuous processing and calculation are performed according to the inter-layer coupling effect; S5: For the global stiffness matrix A global After regularization, we can obtain ; S6: Based on and M global Construct the dispersion equation, solve the dispersion equation to obtain the eigenvalue λ and eigenvector Ф, and calculate the phase velocity v of the Lamb wave p and group velocity v g .

2. A multilayer composite material damage analysis method based on numerical simulation calculation according to claim 1, characterized in that: The S1 comprises the following steps: S11: Construct the original single-layer stiffness matrix C in the following form: ， in, ， ， ， ， ， ， ， ， ， in, , , , are the Poisson's ratios of the xy, yx, xz, and yz planes, respectively. , , are the elastic moduli in the x, y, and z directions, respectively, , , denote the shear modulus in the xy, yz, and xz planes, respectively; S12: Construct the rotation matrix T in the following form: ， in, ， ， Where θ is the rotation angle; S13: Calculate the initial single-layer stiffness matrix C' according to the following formula: 。 3. A multilayer composite material damage analysis method based on numerical simulation calculation according to claim 1, characterized in that: The S2 comprises the following steps: S21: Using the Gauss-Legendre integral method, the single-layer stiffness matrix A of each layer is calculated by the following formula layer and the single-layer mass matrix M layer The elements in: ， ， Where N is the number of integration points, i, j = 0, 1, ..., N-1, a and b are the index values ​​of the Gauss-Legendre integration points, is the N-1 order Legendre polynomial basis function, is the weight, h is the layer thickness, is the density.

4. A multilayer composite material damage analysis method based on numerical simulation calculation according to claim 3, characterized in that: The S3 comprises the following steps: S31: For the top layer of a multilayer composite, the free boundary condition is set by the following formula: ， ， S32: For the bottom layer of a multilayer composite material, the free boundary condition is set by the following formula: ， ， Wherein, M is the number of layers of the multilayer composite material.

5. The damage analysis method of multilayer composite materials based on numerical simulation calculation according to claim 1, characterized in that: The S4 comprises the following steps: S41: The global stiffness matrix of the first layer is accumulated according to the following formula: and the global mass matrix , ， ， S42: The global stiffness matrix of the next layer is accumulated according to the following formula: and the global mass matrix , ， ， Where, c=1,2,...M, M is the number of layers of the multilayer composite material; S43: By superimposing adjacent layer matrix blocks A upper and A lower , update the global stiffness matrix and the global mass matrix The elements satisfy the displacement and stress continuity conditions; S44: Execute S42 and S43 repeatedly until the cumulative construction is completed and the global stiffness matrix A is obtained. global and the global mass matrix M global .

6. A multilayer composite material damage analysis method based on numerical simulation calculation according to claim 5, characterized in that: The S43 comprises the following steps: S431: Construct the matrix block A at the bottom of the previous layer upper , and calculate A according to the following formula upper Elements in: ， Where F is the number of basis functions used in discretization; S432: Build the matrix block A at the top of the current layer lower , and calculate A according to the following formula lower Elements in: ， S433: Update according to the following formula Elements in: 。 7. The damage analysis method of multilayer composite materials based on numerical simulation calculation according to claim 1, characterized in that: The S5 comprises the following steps: S51: calculated by the following formula , ， in, =10 -6 , I is in A global The identity matrix of the same order, A global condition number.

8. The damage analysis method of multilayer composite materials based on numerical simulation calculation according to claim 1, characterized in that: The S6 comprises the following steps: S61: Construct the dispersion equation in the following form: ， S62: Solve the dispersion equation and obtain a series of eigenvalues ​​λ= , and the eigenvector Ф corresponding to each one, where ω is the angular frequency; S63: Delete the imaginary eigenvalue λ, and obtain the effective eigenvalue λ and the effective eigenvector Φ; S64: Calculate the phase velocity v of the Lamb wave according to the following formula p and phase velocity v g， ， ， Where k is the wave number of the multilayer composite material and f is the frequency.

9. A multilayer composite material damage analysis method based on numerical simulation calculation according to claim 8, characterized in that: The following steps are also included: S7: With frequency f as x-axis and phase velocity v as p Or phase velocity v g As the y-axis, a visualization image is constructed with the effective eigenvector Ф as the z-axis.

Images