Dynamic Progressive Failure Analysis Method For Composite Multi-Scale Model

Abstract

This patent studies a scale-span modeling method to simulate the structural mechanical responses and dynamic progressive failure behaviors of carbon fiber reinforced plastics (CFRPs) in drilling. Firstly, considering the different mechanical behaviors of fiber and matrix in micro state, a three-dimensional multi-scale dynamic progressive damage evolution model based on micro failure theory is proposed. Based on the degradation elastic parameters of microcomponent in typical volume element model, a new damage evolution model of fiber and resin matrix and an auxiliary deletion criterion of failure element are proposed. Secondly, the relationship between the macro stress and the micro stress of representative volume element in the composite model is established by using the stress amplification factor. Combined with the bilinear cohesion element model, the damage behavior of the composite in and between layers under the cutting action of dagger drill is simulated.

Claims

1. A dynamic progressive failure analysis method based on multi-scale model of composite material is proposed as follows: In the first step, ignoring the initial damage of composite laminates in the process of fabrication, the composite is idealized as the superposition of the preset stacking sequence of unidirectional laminates, and the macro structure of the composite is analyzed, considering the influence of macro structure design, material properties, stacking method and loading form of the whole structure on the internal stress and strain distribution of the structure; In the second step, according to the mechanical analysis of unidirectional laminates, the fiber distribution of unidirectional laminates is reasonably simplified; Based on the fiber volume fraction, a representative volume element of micromechanics oriented to the apparent mechanical properties of composite structures is established, and the multi-scale analysis method is used to transfer the load of macro structural elements or integral points to the micro model, Under this load, the stress distribution of each element integral point of fiber and matrix component material is simulated; In addition, the micro failure theory is used to judge the failure of the matrix and fiber respectively; Assuming that each corresponding element is damaged, the stiffness of the damaged element is reduced and the failure element is deleted until it is deleted; If the element is not deleted, the whole cell model is homogenized to obtain the macro mechanical properties of the damaged material; In the third step, the macro properties of the damaged material are assigned to the corresponding elements in the macro structure, and then the macro analysis of the whole structure is carried out in the next step, which is repeated until the set time analysis step is completed; Finally, the multi-scale numerical simulation analysis of the composite material from the component to the structure under the cutting action of the cutting tool is realized.

2. The method according to claim 1, which is characterized in that the macro stress of the composite laminate structure and the micro stress of each component structure are bridged by the stress amplification factor, and the stress amplification factor can be obtained by the finite element analysis results of representative volume elements including fiber and matrix.

3. The method according to claim 1, which is characterized by comprising the following finite element algorithm: In the first step: At the beginning of the n th step, calculate the total global strain ε.sup.n via adding the previous time strain ε.sup.n−1 and the global strain increment Δε.sup.n; The macroscopic stress for each integration point σ.sup.n is computed by the UD-CFRPs laminate stiffness matrix C.sub.ij.sup.n−1,d stored at the (n−1) th step; In the second step: The microscopic stresses in fiber and matrix are calculated by multiplying the SAFs of the fiber M.sub.fσ.sup.n,pk and matrix M.sub.mσ.sup.n,pk for all the pre-set reference points, respectively; Then apply the modified MMF criterion for the fiber and matrix to judge their failure factor, respectively; The 3D dynamic progressive damage model is employed to gain the damage variables for fiber and matrix; The maximum values of d.sub.f.sup.n,T(C) and d.sub.m.sup.n,T(C) for all pre-set reference points are defined for the fiber and matrix damage coefficient of the RVE at microscopic, respectively; In the third step: Based on the established damage evolution model, the damage variables at macroscale for each integration point could be computed via the maximum values of d.sub.f.sup.n,T(C) and d.sub.m.sup.n,T(C) to realize the stiffness degradation of elements; Larger value of the macroscopic damage variables at the nth step are chosen by comparing with that at the n−1 step since they are irreversible; After that, the failure element deletion auxiliary criterion is added to implement the serious distortion element deletion in order to prevent nonconvergence in the calculation; Then the damaged UD-CFRPs laminate stiffness matrix C.sub.ij.sup.n,d is re-calculated and updated to determine whether the element is required to be deleted or not; The procedure will end if the pre-set time increments is reached.

4. the method according to claim 1, which is characterized in that the expression of macro micro transformation equation is shown in formula (1), (2) and (3): $\begin{array}{cc}{\sigma }_i=M_{\sigma }^{\mathrm{pk}\left(n\right)}{\stackrel{\_}{\sigma }}_j& \left(1\right)\\ \{\begin{array}{c}{\sigma }_i=[\begin{array}{cccccc}{\sigma }_1& {\sigma }_2& {\sigma }_3& {\sigma }_4& {\sigma }_5& {{\sigma }_6]}^T\end{array}\\ {\stackrel{\_}{\sigma }}_j=[\begin{array}{cccccc}{\stackrel{\_}{\sigma }}_1& {\stackrel{\_}{\sigma }}_2& {\stackrel{\_}{\sigma }}_3& {\stackrel{\_}{\sigma }}_4& {\stackrel{\_}{\sigma }}_5& {{\stackrel{\_}{\sigma }}_6]}^T\end{array}\end{array}& \left(2\right)\\ M_{\sigma }^{\mathrm{pk}\left(n\right)}={\left[\begin{array}{cccccc}{M}_{\sigma }^{11}& M_{\sigma }^{12}& M_{\sigma }^{13}& M_{\sigma }^{14}& M_{\sigma }^{15}& M_{\sigma }^{16}\\ M_{\sigma }^{21}& M_{\sigma }^{22}& M_{\sigma }^{23}& M_{\sigma }^{24}& M_{\sigma }^{25}& M_{\sigma }^{26}\\ M_{\sigma }^{31}& M_{\sigma }^{32}& M_{\sigma }^{33}& M_{\sigma }^{33}& M_{\sigma }^{34}& M_{\sigma }^{35}\\ M_{\sigma }^{41}& M_{\sigma }^{42}& M_{\sigma }^{43}& M_{\sigma }^{44}& M_{\sigma }^{45}& M_{\sigma }^{46}\\ M_{\sigma }^{51}& M_{\sigma }^{52}& M_{\sigma }^{53}& M_{\sigma }^{54}& M_{\sigma }^{55}& M_{\sigma }^{56}\\ M_{\sigma }^{61}& M_{\sigma }^{62}& M_{\sigma }^{63}& M_{\sigma }^{64}& M_{\sigma }^{65}& M_{\sigma }^{66}\end{array}\right]}^{\left(i\right)}& \left(3\right)\end{array}$ Where in Eq. (1)-(3), “σ.sub.i” denotes the microscopic stress at point n of fiber or matrix, “n” denotes the number of the pre-set reference point, i=1˜6; “σ.sub.j” denotes the macroscopic stress of CFRPs laminate, j=1˜6; “M.sub.σ.sup.pk(n)” denotes the SAFs at point n under the mechanical loading, which establish the relationship between the macroscopic and microscopic levels, can be computed by RVE, p=1˜6, k=1˜6.

5. The method according to claim 1, which is characterized in that the equation for calculating the element characteristic length can be expressed as follows: $\begin{array}{cc}L=\sqrt{\frac{{\left(L_{\mathrm{initial}}\right)}^3}{i_z}}& \left(4\right)\end{array}$ Where “L.sub.initial” denotes the characteristic length of an element, In the material subroutine, “charlength” function is embedded in VUMAT; “l.sub.z” denotes the dimension in the thickness direction of an element.

6. The dynamic progressive failure analysis method of composite multi-scale model according to claim 1, which is characterized in that in order to prevent non convergence in calculation due to severe distortion elements in the finite element software, the application introduces the maximum principal strain and minimum principal strain criteria to realize auxiliary deletion of excessive distortion elements, and is based on the polar decomposition theorem, The motion of any object can be decomposed into pure stretching and pure rigid body rotation in three directions: $\begin{array}{cc}\{\begin{array}{c}\tilde{F}=\tilde{V}.Math.\tilde{R}\\ \tilde{F}=\frac{\partial x}{\partial \stackrel{\_}{X}}\end{array}& \left(5\right)\end{array}$ “{tilde over (F)}” represents the strain gradient matrix; “{tilde over (V)}” represents the pure tensile deformation of the material; “{tilde over (R)}” represents the pure rigid body rotation of the material; “{tilde over (X)}” and “x” represent the initial position and the deformed position of the micro-element, respectively; The macroscopic damage-constitutive stiffness matrix for CFRPs is updated according to reference when damage occurs, which is written as $\begin{array}{cc}{C}_{\mathrm{ij}}^d=\left[\begin{array}{ccc}\left(1-d_L\right)C_{11}^0& \left(1-d_L\right)\left(1-d_T\right)C_{12}^0& \left(1-d_L\right)C_{13}^0\\ \left(1-d_L\right)\left(1-d_T\right)C_{12}^0& \left(1-d_T\right)C_{22}^0& \left(1-d_T\right)C_{23}^0\\ \left(1-d_L\right)C_{13}^0& \left(1-d_T\right)C_{23}^0& C_{33}^0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ \left(1-d_t\right)C_{44}^0& 0& 0\\ 0& \left(1-d_L\right)C_{55}^0& 0\\ 0& 0& \left(1-d_L\right)\left(1-d_T\right)C_{66}^0\end{array}\right]& \left(6\right)\end{array}$ Where “d.sub.L” and “d.sub.T” denote the longitudinal and transverse macroscopic damage variables of CFRPs laminate, respectively, “C.sub.ij.sup.0” (i, j=1˜6) denotes the undamaged macroscopic stiffness coefficient for CFRPs laminate, and the variable equations are written as $\left(7\right)$ $\{\begin{array}{c}\begin{array}{ccc}{C}_{11}^0=\frac{E_1\left(1-v_{23}{v}_{32}\right)}{\Lambda }& C_{12}^0=\frac{E_2\left(1-v_{13}{v}_{32}\right)}{\Lambda }& C_{13}^0=\frac{E_3\left(v_{13}-v_{12}{v}_{23}\right)}{\Lambda }\\ C_{22}^0=\frac{E_2\left(1-v_{13}{v}_{31}\right)}{\Lambda }& C_{23}^0=\frac{E_3\left(v_{23}-v_{21}{v}_{31}\right)}{\Lambda }& C_{33}^0=\frac{E_3\left(1-v_{12}{v}_{21}\right)}{\Lambda }\\ C_{44}^0\Lambda G_{23}& C_{55}^0\Lambda G_{13}& C_{66}^0\Lambda G_{12}\\ & & \end{array}\\ \Lambda =1-v_{12}{v}_{21}-v_{23}{v}_{32}-v_{31}{v}_{13}-2v_{12}{v}_{23}{v}_{31}\end{array}$ Where “E.sub.i”, “v.sub.ij” and “G.sub.ij” (i,j=1,2,3) denote the macroscopic elastic modulus, Poisson's ratio, macroscopic shear modulus of CFRPs, respectively; The macroscopic damage variables are directly evaluated by the degraded elastic parameters of the microscopic constituents computed from RVE via the mixture rule, which is given as $\begin{array}{cc}\{\begin{array}{c}{d}_L^{T\left(C\right)}=1-\frac{\left(E_{f1}^{d,T\left(C\right)}{V}_f+E_m^{d,T\left(C\right)}{V}_m\right)}{E_1}\\ d_T^{T\left(C\right)}=1-\frac{\left[\begin{array}{c}\left(1-C\right)\frac{E_{f2}^{d,T\left(C\right)}{E}_m}{E_{f2}^{d,T\left(C\right)}{V}_m+E_m^{d,T\left(C\right)}{V}_f}\\ +C\left(E_{f2}^{d,T\left(C\right)}{V}_f+E_m^{d,T\left(C\right)}{V}_m\right)\end{array}\right]}{E_2}\\ C=0.2\left(V_f-V_m\right)+0.175\end{array}& \left(8\right)\\ \{\begin{array}{c}{d}_L=1-\left(1-d_L^T\right)\left(1-d_L^C\right)\\ d_T=1-\left(1-d_T^T\right)\left(1-d_T^C\right)\end{array}& \left(9\right)\end{array}$ Where “V.sub.f” and “V.sub.m” denote the fiber volume fraction and matrix volume fraction of CFRPs, respectively; “E.sub.f1.sup.d,T(C)” and “E.sub.f2.sup.d,T(C)” denote the longitudinal and transverse damaged fiber elastic modulus at microscopic, respectively; “E.sub.f1” and “E.sub.f2” denote the initial fiber elastic modulus at microscopic, respectively, “E.sub.m.sup.d,T(C)” indicates the damage state variable of the matrix under tensile and compression conditions.

Description

BRIEF DESCRIPTION OF THE DRAWINGS:

[0008] FIG. 1 Scale-span FE drilling schematic diagram.

[0009] FIG. 2 RVE model and the responding reference points in fiber and matrix constituents.

[0010] FIG. 3 The stress-strain behavior in fiber and matrix tension and compression cases.

[0011] FIG. 4 Schematic diagram of element characteristic length calculation and fiber failure state.

[0012] FIG. 5 Numerical algorithm flowchart of VUMAT for Span-scale model in ABAQUS/Explicit.

[0013] FIG. 6 Schematic diagram of bilinear cohesive model .

[0014] FIG. 7 ABAQUS-Python code script flowchart for calculating the SAFs in RVE.

[0015] FIG. 8 Schematic diagram of the layered failure model of the cohesion unit.

[0016] FIG. 9 Schematic diagram of periodic boundary condition loading of representative volume element model.

[0017] FIG. 10 SAFs of the reference points in fiber.

[0018] FIG. 11 SAFs of the reference points in matrix.

[0019] FIG. 12 Schematic of the experimental setup for thrust force in the drilling of CFRPs.

[0020] FIG. 13 Setup showing the Span-scale model of drilling CFRPs.

[0021] FIG. 14 Comparison of the thrust force and torque.

[0022] FIG. 15 Comparison of hole's surface morphology.

[0023] FIG. 16 Comparison of damage morphology at drill entry and exit.

[0024] FIG. 17 Comparison of the delamination factor and deviation.

DETAILED DESCRIPTION:

[0025] The macroscopic stress of the composite laminate structure and the microscopic stress of each component structure are bridged by the stress amplification factor, and the stress amplification factor can be obtained by the finite element analysis result of the representative volume element containing the fiber and the matrix. After the stress amplification factor is obtained, the macro-micro conversion equation is used to realize the conversion of macro-stress to micro-stress. The expression of the macro-micro conversion equation is shown in formulas (1), (2), (3).

[00001]$\begin{array}{cc}{\sigma }_i=M_{\sigma }^{pk\left(n\right)}{\overline{\sigma }}_j& \left(10\right)\end{array}$$\begin{array}{cc}\{\begin{array}{c}{\sigma }_i={\left[{\sigma }_1{\sigma }_2{\sigma }_3{\sigma }_4{\sigma }_5{\sigma }_6\right]}^T\\ {\overline{\sigma }}_j={\left[{\overline{\sigma }}_1{\overline{\sigma }}_2{\overline{\sigma }}_3{\overline{\sigma }}_4{\overline{\sigma }}_5{\overline{\sigma }}_6\right]}^T\end{array}& \left(11\right)\end{array}$$\begin{array}{cc}{M}_{\sigma }^{pk\left(n\right)}={\left[\begin{array}{cccccc}{M}_{\sigma }^{11}& M_{\sigma }^{12}& M_{\sigma }^{13}& M_{\sigma }^{14}& M_{\sigma }^{15}& M_{\sigma }^{16}\\ M_{\sigma }^{21}& M_{\sigma }^{22}& M_{\sigma }^{23}& M_{\sigma }^{24}& M_{\sigma }^{25}& M_{\sigma }^{26}\\ M_{\sigma }^{31}& M_{\sigma }^{32}& M_{\sigma }^{33}& M_{\sigma }^{33}& M_{\sigma }^{34}& M_{\sigma }^{35}\\ M_{\sigma }^{41}& M_{\sigma }^{42}& M_{\sigma }^{43}& M_{\sigma }^{44}& M_{\sigma }^{45}& M_{\sigma }^{46}\\ M_{\sigma }^{51}& M_{\sigma }^{52}& M_{\sigma }^{53}& M_{\sigma }^{54}& M_{\sigma }^{55}& M_{\sigma }^{56}\\ M_{\sigma }^{61}& M_{\sigma }^{62}& M_{\sigma }^{63}& M_{\sigma }^{64}& M_{\sigma }^{65}& M_{\sigma }^{66}\end{array}\right]}^{\left(i\right)}& \left(12\right)\end{array}$

[0026] Where in Eq. (1)-(3), σ.sub.i denotes the microscopic stress at point n of fiber or matrix, n denotes the number of the pre-set reference point, i=1˜6. σ.sub.j denotes the macroscopic stress of CFRPs laminate, j=1˜6. M.sub.σ.sup.pk(n) denotes the SAFs at point n under the mechanical loading, which establish the relationship between the macroscopic and microscopic levels, can be computed by RVE. p=1˜6, k=1˜6.

[0027] Based on the failure theory of micromechanics, by setting some reference points to the representative volume element model, the stress amplification factor method is used to extract the microscopic stress of the fiber and the matrix. As shown in FIG. 2, according to the observation of the interface of the unidirectional composite laminate by the ultra-depth-of-field electron microscope, the staggered arrangement of the representative volume unit model of the rectangular structure is closer to the random distribution of fibers, and the rectangular model (FIG. 2c) Compared with the square representative volume unit model (FIG. 2b), it is more accurate, so this application will use the rectangular representative volume unit model for analysis.

[0028] Load periodic boundary conditions on the representative volume element model, where the displacement relationship of each pair of nodes on the opposite surface in the representative volume element model satisfies formulas (4):

u.sub.i.sup.j+(x, y, z)−u.sub.i.sup.j−(x, y, z)=c.sub.i.sup.j(i, j=1,2,3)   (13)

[0029] In formula (4): u.sub.i.sup.j+(x, y, z), u.sub.i.sup.j−(x, y, z) respectively represent the displacement of nodes in the x, y, and z planes along the j-direction normal vector boundary surface in the i direction of the representative volume element model; the superscripts “+” and “−” represent opposite edges Node pairs on the interface; c.sub.i.sup.j is the difference in the displacement of the node pairs in the i direction relative to the boundary surface, which is mainly determined by the average strain of the representative volume element; assuming that the average strain of the representative volume element has been determined, c.sub.i.sup.j changes Constant.

[0030] The representative key points for selecting the structure of each component in the representative volume element model are shown in FIG. 2. Among them, the fiber selects 10 reference points (F1˜F10), and the matrix selects 15 reference points (M1-M15).

[0031] M.sub.σ.sup.pk(n) is obtained through the key reference points set in the linear finite element analysis in the representative volume element model, and the obtained equation is:

[00002]$\begin{array}{cc}{M}_{\sigma }^{pk\left(i\right)}=\psi /\overline{\psi }& \left(5\right)\end{array}$$\begin{array}{cc}\{\begin{array}{c}\psi ={\left[\begin{array}{cccccc}{\zeta }_1& {\zeta }_2& {\zeta }_3& {\zeta }_4& {\zeta }_5& {\zeta }_6\end{array}\right]}_{6\times 6}\\ \stackrel{\_}{\psi }={\left[\begin{array}{cccccc}{\stackrel{\_}{\zeta }}_1& {\stackrel{\_}{\zeta }}_2& {\stackrel{\_}{\zeta }}_3& {\stackrel{\_}{\zeta }}_4& {\stackrel{\_}{\zeta }}_5& {\stackrel{\_}{\zeta }}_6\end{array}\right]}_{6\times 6}\end{array}& \left(6\right)\end{array}$

[0032] Where ξ.sub.i (i=1,2,3,4,5,6) denote the microscopic stress arrays of each reference points under six macroscopic stress load cases, and ξ.sub.i (i=1,2,3,4,5,6) denote the six macroscopic uniform stress load cases.

[0033] Meanwhile, ξ.sub.i and ξ.sub.i are assigned as follows:

[00003]$\begin{array}{cc}{\zeta }_i=\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{12}\\ {\sigma }_{23}\\ {\sigma }_{13}\end{array}\right]{\overline{\zeta }}_i=\left[\begin{array}{c}{\stackrel{\_}{\sigma }}_{11}\\ {\stackrel{\_}{\sigma }}_{22}\\ {\stackrel{\_}{\sigma }}_{33}\\ {\stackrel{\_}{\sigma }}_{12}\\ {\stackrel{\_}{\sigma }}_{23}\\ {\stackrel{\_}{\sigma }}_{13}\end{array}\right]i=1,2,3,4,5,6& \left(7\right)\end{array}$

[0034] Where σ.sub.i (i=11, 22, 33, 12, 23, 13) denotes the obtained microscopic stress of the RVE model under six macroscopic stress load cases, and σ.sub.i (i=11, 22, 33, 12, 23, 13) denotes the loaded six macroscopic uniform stress load cases. 11, 22, 33, 12, 23, 13 denote three uniaxial tensile along X, Y, Z direction and the macroscopic shear along Z, X, Y plane, respectively, which is shown in Error! Reference source not found. Taking X direction tensile for example, macroscopic stress ξ.sub.i equals to [1 0 0 0 0 0].sup.T. The microscopic stress could be extracted from the reference points in RVE model under the unit stress in six load cases, and the SAFs in X direction are easy to be obtained from the microscopic stresses with Eq. (5). Similarly, the analogous approach is adopted to obtain the SAFs for other load cases.

[0035] Due to it is 3D solid element, the equation involves six stress components, which is written as follows:

f.sub.1σ.sub.1+f.sub.2σ.sub.2+f.sub.3σ.sub.3+f.sub.11σ.sub.1.sup.2+f.sub.22σ.sub.2.sup.2+f.sub.33σ.sub.3.sup.2+2 f.sub.12σ.sub.1σ.sub.2+2 f.sub.23σ.sub.2σ.sub.3+2 f.sub.13σ.sub.3σ.sub.1+f.sub.44σ.sub.4.sup.2+f.sub.55σ.sub.5.sup.2+f.sub.66σ.sub.6.sup.2=1   (8)

[0036] Where f.sub.i and f.sub.ij denote the strength tensor in i direction and ij plane, respectively. i,j=1,2,3 (i denotes fiber longitudinal direction). σ.sub.i denote the microscopic normal stress of fiber. In addition, the equation of each strength tensor could be written by

[00004]$\begin{array}{cc}\{\begin{array}{c}{f}_1=\frac{1}{X_f^T}-\frac{1}{X_f^C};f_2=\frac{1}{Y_f^T}-\frac{1}{Y_f^C};f_3=\frac{1}{Z_f^T}-\frac{1}{Z_f^C}\\ f_{11}=\frac{1}{X_f^T{X}_f^C};f_{22}=\frac{1}{Y_f^T{Y}_f^C};f_{33}=\frac{1}{Z_f^T{Z}_f^C}\\ f_{44}=\frac{1}{{\left(S_f^{YZ}\right)}^2};f_{55}=\frac{1}{{\left(S_f^{XZ}\right)}^2};f_{66}=\frac{1}{{\left(S_f^{XY}\right)}^2}\\ \begin{array}{c}{f}_{12}=\frac{-1}{2\sqrt{X_f^T{X}_f^C{X}_f^T{X}_f^C}};f_{23}=\frac{-1}{2\sqrt{Y_f^T{Y}_f^C{Y}_f^T{Y}_f^C}};\\ f_{13}=\frac{-1}{2\sqrt{Z_f^T{Z}_f^C{X}_f^T{X}_f^C}}\end{array}\end{array}& \left(9\right)\end{array}$

[0037] Where X, Y and Z denote the failure strength in each direction, respectively. The superscripts T and C denote tensile and compressible conditions, respectively. The subscripts f denotes the microscopic fiber. S denotes the shear failure strength in each direction, and i, j=X, Y, Z.

[0038] To simulate the real failure behavior of fiber in dynamic drilling condition more accurately, the continuous degradation law for fiber tension and compression are considered by

[00005]$\begin{array}{cc}{d}_f^{T\left(C\right)}=\frac{{\epsilon }_{f,1}^{T\left(C\right)}}{{\epsilon }_{f,1}^{T\left(C\right)}-{\epsilon }_{0,1}^{T\left(C\right)}}\left(1-\frac{{\epsilon }_{0,1}^{T\left(C\right)}}{{\epsilon }_{11}}\right)& \left(10\right)\end{array}$

[0039] Where d.sub.f.sup.T(C) denotes the damage coefficient of each reference points of fiber in RVE model. ε.sub.f,1.sup.T(C) denotes the critical failure strains when the damage variables reach one. ε.sub.11 denotes the real strain. ε.sub.0,1.sup.T(C) denotes the initial strain.

[0040] The equations which is used to obtain the damage variable are written by

[00006]$\begin{array}{cc}{\epsilon }_{f,1}^{T\left(C\right)}=\frac{2{\Gamma }_f^{T\left(C\right)}}{X_f^{T\left(C\right)}L}& \left(14\right)\end{array}$$\begin{array}{cc}{\epsilon }_{0,1}^{T\left(C\right)}=\frac{X_f^{T\left(C\right)}}{E_{f1}}& \left(15\right)\end{array}$

[0041] Where in Eq. (14) and Eq. (15), Γ.sub.f.sup.T denotes the critical failure strain by regularizing the fiber fracture energy, and L denotes the characteristic length of microscopic element.

[0042] Due to the hexahedral element type is adopted in this study, and elements that may be damaged have the same in-plane length, thus the equation is written by

[00007]$\begin{array}{cc}L=\sqrt{\frac{{\left(L_{\mathrm{initial}}\right)}^3}{l_z}}& \left(16\right)\end{array}$

[0043] Where L.sub.initial denotes the characteristic length of an element, l.sub.z denotes the dimension in the thickness direction of an element as shown in Error! Reference source not found. (a).

[0044] For the tension case, fibers are completely failure in direction X if the damage coefficient d.sub.f.sup.T reach 1. For the compression case, the elements are considered to exist some residual loading capacity.

[0045] The fiber constitutive relationship at microscale is expressed as follows:

[00008]$\begin{array}{cc}\{\begin{array}{c}{\sigma }_f=\left(1-d_f^{T\left(C\right)}\right)C_f{\epsilon }_f\\ E_{f1}^{d,T\left(C\right)}=\left(1-d_f^{T\left(C\right)}\right)E_{f1}\\ E_{f2}^{d,T\left(C\right)}=\left(1-d_f^{T\left(C\right)}\right)E_{f2}\end{array}& \left(17\right)\end{array}$

[0046] Where σ.sub.f, C.sub.f and ε.sub.f denote the microscopic stresses, intact stiffness and microscopic strains for fiber, respectively. E.sub.f1.sup.d,T(C) and E.sub.f2.sup.d,T(C) denote the longitudinal and transverse damaged fiber elastic modulus at microscopic, respectively. E.sub.f1 and E.sub.f2 denote the initial fiber elastic modulus at microscopic, respectively.

[0047] A modified Von Mises yield criterion is used to predict the tensile and compressive failure criterion of matrix in this study, which is shown in Error! Reference source not found. (b), and the equation is written as follows:

[00009]$\begin{array}{cc}\frac{{\sigma }_{VM}^2+I_1\left(C_{mi}-T_{mi}\right)}{C_{mi}{T}_{mi}}=1& \left(18\right)\end{array}$

[0048] Where T.sub.mi and C.sub.mi denote the initial tensile strength and compression strength, respectively. σ.sub.VM and I.sub.1 denote equivalent stress and first stress invariants, respectively, these equations could be written by

[00010]$\begin{array}{cc}\{\begin{array}{c}{\sigma }_{VM}=\sqrt{\begin{array}{c}\frac{1}{2}[{\left({\sigma }_{m1}-{\sigma }_{m2}\right)}^2+{\left({\sigma }_{m2}-{\sigma }_{m3}\right)}^2+\\ {\left({\sigma }_{m1}-{\sigma }_{m3}\right)}^2+3\left({\sigma }_{m4}^2+{\sigma }_{m5}^2+{\sigma }_{m6}^2\right)]\end{array}}\\ I_1={\sigma }_{m1}+{\sigma }_{m2}+{\sigma }_{m3}\end{array}& \left(19\right)\end{array}$

[0049] Where σ.sub.mi denotes the principal stress (i=1˜3) and shear stress (i=4˜6) of matrix at microscopic.

[0050] Due to there is a great difference between tensile strength and compression strength, an equivalent stress σ.sub.eq called the Stassi's stress is adopted to define the matrix damage evolution.

[00011]$\begin{array}{cc}{\sigma }_{eq}=\frac{\left(\beta -1\right)I_{1^{+}}\sqrt{{\left(\beta -1\right)}^2{I}_1^2+4\beta {\sigma }_{VM}^2}}{2\beta }& \left(20\right)\end{array}$

[0051] Where β is denotes as the ratio of the matrix initial compressive to initial tensile strength.

[0052] The value of matrix damage coefficient is determined by the stress state, which varies between 0 and 1. In RVE, the damage evolution law and the equivalent stress of the matrix could be expressed as follows:

[00012]$\begin{array}{cc}{\overline{d}}_m^{T\left(C\right)}=\frac{\overline{\kappa }\gamma }{T_{mi}}\left(1-d_m^{T\left(C\right)}\right)& \left(21\right)\end{array}$$\begin{array}{cc}f\left({\sigma }_{eq},{\kappa }_m\right)={\sigma }_{eq}-{\kappa }_m& \left(22\right)\end{array}$

[0053] Where d.sub.m.sup.T(C) denotes the matrix damage coefficient which is obtained by the reference points in RVE. γ denotes a parameter used to calibrate the uniaxial stress-strain curves based on test data. κ.sub.m denotes a scalar history variable, defined as the ratio of σ.sub.eq to T.sub.mi. The determination and update of the initial damage are required by the Kuhn-Tucker loading-unload condition, namely:

f≤0, {dot over (κ)}.sub.m≥0, f{dot over (κ)}.sub.m=0   (20)

[0054] Then, the matrix damage coefficient in tensile and compression are expressed as

[00013]$\begin{array}{cc}{d}_m^{T\left(C\right)}=1-\exp \left[\gamma \left(1-\frac{{\sigma }_{eq}}{T_{mi}}\right)\right]& \left(23\right)\end{array}$

[0055] Under the consideration that the damage initial and evolution laws are isotropic, the constitutive relationship of the matrix by considering the degradation of the stiffness at microscopic is shown in Error! Reference source not found, and the entire equation is given as

[00014]$\begin{array}{cc}\{\begin{array}{c}{\sigma }_m=\left(1-d_m^{T\left(C\right)}\right)C_m{\epsilon }_m\\ E_m^{T\left(C\right)}=\left(1-d_m^{T\left(C\right)}\right)E_m\end{array}& \left(24\right)\end{array}$

[0056] Where σ.sub.m, C.sub.m and ε.sub.m denote the stress, stiffness and strain of matrix at microscopic, respectively. E.sub.m.sup.T(C) and E.sub.m denote the damaged and initial elastic modulus of matrix at microscopic, respectively.

[0057] The macroscopic damage-constitutive stiffness matrix for CFRPs is updated according to reference when damage occurs, which is written as

[00015]$\begin{array}{cc}{C}_{\mathrm{ij}}^d=\left[\begin{array}{ccc}\left(1-d_L\right)C_{11}^0& \left(1-d_L\right)\left(1-d_T\right)C_{12}^0& \left(1-d_L\right)C_{13}^0\\ \left(1-d_L\right)\left(1-d_T\right)C_{12}^0& \left(1-d_T\right)C_{22}^0& \left(1-d_T\right)C_{23}^0\\ \left(1-d_L\right)C_{13}^0& \left(1-d_T\right)C_{23}^0& C_{33}^0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ \left(1-d_t\right)C_{44}^0& 0& 0\\ 0& \left(1-d_L\right)C_{55}^0& 0\\ 0& 0& \left(1-d_L\right)\left(1-d_T\right)C_{66}^0\end{array}\right]& \left(25\right)\end{array}$

[0058] Where d.sub.L and d.sub.T denote the longitudinal and transverse macroscopic damage variables of CFRPs laminate, respectively. C.sub.ij.sup.0 (i, j=1˜6) denotes the undamaged macroscopic stiffness coefficient for CFRPs laminate, and the variable equations are written as

[00016]$\left(26\right)$$\{\begin{array}{c}\begin{array}{ccc}{C}_{11}^0=\frac{E_1\left(1-v_{23}{v}_{32}\right)}{\Lambda }& C_{12}^0=\frac{E_2\left(1-v_{13}{v}_{32}\right)}{\Lambda }& C_{11}^0=\frac{E_3\left(v_{13}-v_{12}{v}_{23}\right)}{\Lambda }\\ C_{22}^0=\frac{E_2\left(1-v_{13}{v}_{31}\right)}{\Lambda }& C_{23}^0=\frac{E_3\left(v_{23}-v_{21}{v}_{31}\right)}{\Lambda }& C_{33}^0=\frac{E_3\left(1-v_{12}{v}_{21}\right)}{\Lambda }\\ C_{44}^0\Lambda G_{23}& C_{55}^0\Lambda G_{13}& C_{66}^0\Lambda G_{12}\\ & & \end{array}\\ \Lambda =1-v_{12}{v}_{21}-v_{23}{v}_{32}-v_{31}{v}_{13}-2v_{12}{v}_{23}{v}_{31}\end{array}$

[0059] Where E.sub.i, v.sub.ij and G.sub.ij (i,j=1,2,3) denote the macroscopic elastic modulus, Poisson's ratio, macroscopic shear modulus of CFRPs, respectively.

[0060] The macroscopic damage variables are directly evaluated by the degraded elastic parameters of the microscopic constituents computed from RVE via the mixture rule, which is given as

[00017]$\begin{array}{cc}\{\begin{array}{c}{d}_L^{T\left(C\right)}=1-\frac{\left(E_{f1}^{d,T\left(C\right)}{V}_f+E_m^{d,T\left(C\right)}{V}_m\right)}{E_1}\\ d_T^{T\left(C\right)}=1-\frac{\left[\begin{array}{c}\left(1-C\right)\frac{E_{f2}^{d,T\left(C\right)}{E}_m}{E_{f2}^{d,T\left(C\right)}{V}_m+E_m^{d,T\left(C\right)}{V}_f}\\ +C\left(E_{f2}^{d,T\left(C\right)}{V}_f+E_m^{d,T\left(C\right)}{V}_m\right)\end{array}\right]}{E_2}\\ C=0.2\left(V_f-V_m\right)+0.175\end{array}& \left(27\right)\\ \{\begin{array}{c}{d}_L=1-\left(1-d_L^T\right)\left(1-d_L^C\right)\\ d_T=1-\left(1-d_T^T\right)\left(1-d_T^C\right)\end{array}& \left(28\right)\end{array}$

[0061] Where V.sub.f and V.sub.m denote the fiber volume fraction and matrix volume fraction of CFRPs, respectively.

[0062] In order to prevent non-convergence in the calculation due to the occurrence of severely distorted elements in the finite element software, this paper introduces the maximum principal strain and minimum principal strain criteria to realize the auxiliary deletion of excessively distorted elements. Based on the polar decomposition theorem, the motion of any object is It can be decomposed into pure stretching and pure rigid body rotation in three directions, namely:

[00018]$\begin{array}{cc}\{\begin{array}{c}\tilde{F}=\tilde{V}.Math.\tilde{R}\\ \tilde{F}=\frac{\partial x}{\partial \stackrel{\_}{X}}\end{array}& \left(29\right)\end{array}$

[0063] {tilde over (F)} represents the strain gradient matrix; {tilde over (V)} represents the pure tensile deformation of the material; {tilde over (R)} represents the pure rigid body rotation of the material; {tilde over (X)} and {tilde over (x)} represent the initial position and the deformed position of the micro-element, respectively.

[0064] The laminated plate models of the composite materials in this application all adopt the hexahedral element model, and the characteristic roots of the {tilde over (V)} matrix are solved, and the main stretching ratios of the three directions of the element λ.sub.1, λ.sub.2, λ.sub.3 can be obtained, which are respectively calculated according to the definition of Biot strain (nominal strain) 3 principal strains.

ε.sub.i=λ.sub.i−1   (30)

[0065] ε.sub.i (i=1,2,3) represents the principal strain of the element; λ.sub.i (i=1,2,3) represents the principal stretch ratio of the element.

[0066] According to the computed principal strain value of each element, the auxiliary element deletion criterion can be expressed as

[00019]$\begin{array}{cc}\{\begin{array}{cc}{\mathrm{.Math.}}_i^{\max }-{\mathrm{.Math.}}_i^{\max ,s}\ge 0& \mathrm{failure}\\ {\mathrm{.Math.}}_i^{\max }-{\mathrm{.Math.}}_i^{\max ,s}<0& \mathrm{no}\mathrm{failure}\\ {\mathrm{.Math.}}_i^{\min }-{\mathrm{.Math.}}_i^{\min ,s}\ge 0& \mathrm{failure}\\ {\mathrm{.Math.}}_i^{\min }-{\mathrm{.Math.}}_i^{\min ,s}<0& \mathrm{no}\mathrm{failure}\end{array}& \left(31\right)\end{array}$

[0067] In Eq. (31), ε.sub.p.sup.max and ε.sub.p.sup.min denote the maximum and minimum principal strain of the failure elements, respectively, ε.sub.p.sup.max,s and ε.sub.p.sup.min,s denote the pre-set value of the maximum and minimum principal strain according to educated estimation, respectively. ε.sub.o.sup.max,s=1.0 and ε.sub.o.sup.min,s=0.8 are adopted in this study.

[0068] A detailed algorithm flowchart of the logic in the implementation of VUMAT is shown in Error! Reference source not found. The procedure of span-scale numerical implementation is mainly divided into three steps.

[0069] Step (1): At the beginning of the n th step, calculate the total global strain ε.sup.n via adding the previous time strain ε.sup.n−1 and the global strain increment Δε.sup.n. The macroscopic stress for each integration point σ.sup.n is computed by the UD-CFRPs laminate stiffness matrix C.sub.ij.sup.n−1,d stored at the (n−1) th step.

[0070] Step (2): The microscopic stresses in fiber and matrix are calculated by multiplying the SAFs of the fiber M.sub.fσ.sup.n,pk and matrix M.sub.mσ.sup.n,pk for all the pre-set reference points, respectively. Then apply the modified MMF criterion for the fiber and matrix to judge their failure factor, respectively. The 3D dynamic progressive damage model is employed to gain the damage variables for fiber and matrix. The maximum values of d.sub.f.sup.n,T(C) and d.sub.m.sup.n,T(C) for all pre-set reference points are defined for the fiber and matrix damage coefficient of the RVE at microscopic, respectively.

[0071] Step (3): Based on the established damage evolution model, the damage variables at macroscale for each integration point could be computed via the maximum values of d.sub.f.sup.n,T(C) and d.sub.m.sup.n,T(C) to realize the stiffness degradation of elements. Larger value of the macroscopic damage variables at the nth step are chosen by comparing with that at the (n−1) step since they are irreversible. After that, the failure element deletion auxiliary criterion is added to implement the serious distortion element deletion in order to prevent nonconvergence in the calculation. Then the damaged UD-CFRPs laminate stiffness matrix C.sub.ij.sup.n,d is re-calculated and updated to determine whether the element is required to be deleted or not. The procedure will end if the pre-set time increments is reached.

[0072] The initial and damage constitutive expression can be written as

[00020]$\begin{array}{cc}\{\begin{array}{cc}{T}_m=K_i{\delta }_m,& {\delta }_m^{\max }\le {\delta }_m^0\\ T_m=\left(1-D_{\mathrm{mix}}^s\right)K_i{\delta }_m,& {\delta }_m^0\le {\delta }_m^{\max }<{\delta }_m^f\\ T_m=0,& {\delta }_m^{\max }\ge {\delta }_m^f\end{array}& \left(30\right)\end{array}$

[0073] Where T.sub.m denotes the mixed-mode nominal stress. K.sub.i (i=n, s, t denotes the mode N, model S and model T delamination, respectively) denotes the constitutive stiffness of CEs. D.sub.mix.sup.s denotes the mixed-mode damage variables of the CEs. δ.sub.i denotes the mixed-mode displacement. δ.sub.m.sup.max denotes the maximum value of the mixed-mode displacement. δ.sub.m.sup.0 denotes the effective displacement at the damage initiation. δ.sub.m.sup.f denotes the mixed-mode displacement at complete failure. Each expression can be written as follows:

[00021]$\begin{array}{cc}{D}_{\mathrm{mix}}^s=\frac{{\delta }_m^f\left({\delta }_m^{\max }-{\delta }_m^0\right)}{{\delta }_m^{\max }\left({\delta }_m^f-{\delta }_m^0\right)}& \left(31\right)\\ \{\begin{array}{c}{\delta }_m^{\max }=\max \left({\delta }_m^{\max },{\delta }_m\right)\\ {\delta }_m=\sqrt{{\delta }_n^2+{\delta }_s^2+.Math.{\delta }_t^2.Math.}\end{array}& \left(32\right)\\ {\delta }_m^0=\{\begin{array}{cc}{\delta }_t^0{\delta }_n^0\sqrt{\frac{1+{\beta }^2}{{\left({\delta }_2^0\right)}^2+{\left({\mathrm{\beta \delta }}_t^0\right)}^2}}& {\delta }_t>0\\ \sqrt{{\left({\delta }_n^0\right)}^2+{\left({\delta }_s^2\right)}^2}& {\delta }_t\le 0\end{array}& \left(33\right)\\ {\delta }_m^f=\{\begin{array}{cc}\frac{2}{K_i{\delta }_m^0}\left[G_{\mathrm{Tn}}+\left(G_{\mathrm{Ts}}-G_{\mathrm{Tn}}\right){\left(\frac{{\beta }^2}{1+{\beta }^2}\right)}^{\eta }\right]& {\delta }_t>0\\ \sqrt{{\left({\delta }_n^0\right)}^2+{\left({\delta }_s^0\right)}^2}& {\delta }_t\le 0\end{array}& \left(32\right)\end{array}$

[0074] Where in Eq. (32), Eq. (33) and Eq. (32), <.circle-solid.> denotes the MacAuley bracket. δ.sub.n, δ.sub.s, δ.sub.t denote the displacement in different directions (n denotes normal direction, s denotes the first shear direction and t denotes the second shear direction), respectively.) β denotes the mode mixing ratio, and it can be written as

[00022]$\begin{array}{cc}\beta =\sqrt{\frac{\left({\delta }_n^2+{\delta }_s^2\right)}{{\delta }_t^2}}{\delta }_t>0& \left(33\right)\end{array}$

[0075] Where η denotes the parameter defined in the Benzeggagh-Kenane fracture energy law, and η=1.45 is adopted due to CFRPs being the research object. G.sub.Tn and G.sub.Ts denote the interlaminar fracture energy in mode N and mode S, respectively. (Mode N and mode S denote crack when the upper and lower surfaces of elements are subject to opposite normal and tangential force displacement, respectively. In addition, Turon et al suggested that the CEs zone length is chosen as approximately 2-3 element side length, which is apply for regularizing the delamination fracture toughness for the sake of eliminating the mesh sensitivity.

[0076] The relationship between actual stress and maximum stress can be expressed as

σ=(1−D.sup.s)σ  (36)

[0077] In the calculation process of the finite element software, if the damage variable D.sup.s at the cross-section integration point of a cohesive element reaches D.sub.max.sup.s, the element will be deleted. The dimensionless stiffness degradation coefficient SDEG in the ABAQUS software and the damage variable D.sup.s in the model have the same meaning, that is, when SDEG=1, the unit is judged to be invalid and deleted.

[0078] The material properties of fiber and matrix are shown in Tab. 1.

TABLE-US-00001 TABLE 1 Material parameters of fiber and matrix Fiber Matrix E.sub.f1(GPa) 230 X.sub.f.sup.T/X.sub.f.sup.C 4.9/4.5 E.sub.m(GPa) 2.9 (MPa) E.sub.f2/E.sub.f3(GPa) 15 Y.sub.f.sup.T/Y.sub.f.sup.C 40.0/70.0 v.sub.m 0.34 (MPa) v.sub.f12/v.sub.f13 0.21 Z.sub.f.sup.T/Z.sub.f.sup.C 40.0/70.0 G.sub.m(GPa) 1.31 (MPa) v.sub.f23 0.307 S.sub.f.sup.XY/S.sub.f.sup.XZ 100.0/100.0 T.sub.mi(MPa) 80.0 (MPa) G.sub.f12/G.sub.f13(GPa) 9 S.sub.f.sup.YZ  58.0 C.sub.mi(MPa) 120.0 (MPa) G.sub.f23(GPa) 5.03 Γ.sub.f.sup.T/Γ.sub.f.sup.C 100.0 γ 1.5 (N/mm)

[0079] When the representative volume element model is meshed, the position and number of the element nodes set on each symmetry plane should be kept completely consistent, and the mesh should be generated by mapping, as shown in FIG. 9(a). In order to avoid over-constrained mesh nodes on the vertices and edges in the representative volume element model, when applying periodic boundary conditions to the node set, the vertices and nodes on the edges must be removed from the surface, and the vertices must be removed at the same time. The node on is removed from the edge, a vertex contains only one node, but each node set on all faces and edges can contain multiple nodes at the same time.

[0080] In the global coordinate system, the side lengths of the representative volume element model established are defined as W.sub.x, W.sub.y, and h. The coordinate origin is located at D. Under six kinds of macro stress loads (ε.sub.x.sup.0, ε.sub.y.sup.0, ε.sub.z.sup.0, ε.sub.xy.sup.0, ε.sub.xz.sup.0, ε.sub.yz.sup.0), the period established according to formula (4) Under nonlinear boundary conditions, the following linear constraint equations are used to realize the loading of vertices, edges and planes in response to load conditions.

[0081] In the ABAQUS/Standard analysis step, the representative volume element model is respectively applied to each unit stress load under six working conditions (triaxial macroscopic tension and shear) to obtain the microscopic stress distribution, as shown in FIG. 10. 10 red solid circles are used to represent the key reference points of micro-fibers and 15 purple solid circles are used to represent the stress magnification factor calculated from the reference points of the micro-matrix of the matrix, as shown in FIG. 2. The F10 reference point and the M10˜M15 points are located at the center of the fiber and the center of the matrix, respectively. F1˜F9 reference points and M1˜M9 consider the key reference points of the position where the fiber and the matrix join. Since the external stress load is 1 MPa, the microscopic stress distribution of each key reference point is the stress amplification factor. The stress of each key reference point on the fiber and matrix under different working conditions is shown in FIG. 10 and FIG. 11. Show. Since the inherent properties of the material will not change, the representative volume element model only needs to be solved once and the key reference point results are extracted. After obtaining the stress amplification factor, store it in a separate file in the Fortran code. In the macro composite drilling finite element model, use the “#include” function to call the file based on the ABAQUS VUMAT code to realize the dynamics of the composite material. Multi-scale analysis of progressive failure.

[0082] High temperature-resistant CFRPs (T700S-12K/YP-H26) are used as the research object in the FE modeling. The fiber volume fraction is approximately 59%, and the layup sequence of the CFRPs is [(0°/90°/45°/−45°)s].sub.2. Orthotropic material property was assigned to each UD-CFRPs laminate according to the fiber orientation by using a predefined local coordinate system. The CEs, which is modeled as having a thickness 0 mm, was used for simulating the delamination phenomenon. The materials parameters of UD-CFRPs laminate elements and cohesive-zone elements are reported in Tab. 2 and Tab. 3, respectively.

TABLE-US-00002 TABLE 2 Material parameters used to model unidirectional CFRPs laminate elements Value Strength Value Elastic parameters (GPa) parameters (GPa) E.sub.1 138.7 X.sup.T/X.sup.C 1870/1026  E.sub.2 = E.sub.3 7.04 Y.sup.T/Y.sup.C 45/156 v.sub.12 = v.sub.13 0.25 Z.sup.T/Z.sup.C 40/145 v.sub.23 0.31 S.sub.XY 87 G.sub.12 = G.sub.13 2.959 S.sub.XZ 87 G.sub.23 2.505 S.sub.YZ 58

TABLE-US-00003 TABLE 3 Material parameters used to model interface CEs Value Strength Value Fraction Value Stiffness (N/mm.sup.3) parameters (MPa) energy (N/mm) K.sub.n 4 × 10.sup.6 δ.sub.n 60 G.sub.n 0.2 K.sub.s = K.sub.t 1 × 10.sup.6 δ.sub.s = δ.sub.t 90 G.sub.s = G.sub.t 1

[0083] According to the actual working conditions of the drill bit feeding the composite material along the axial direction, boundary conditions such as speed and feed rate are applied to the entire drilling finite element model. Since the motion state of the dagger drill model is constrained at the top reference point, the reference point is The displacement in the X and Y directions is limited, and the feed speed is applied in the Z direction. Similarly, the rotation speed in the X and Y directions is limited, and the clockwise rotation speed is applied in the Z direction, and the composite material layer The four vertical surfaces of the plywood are fixed. After the simulation calculation is completed, this application verifies the model through an experimental platform built. The experimental platform mainly includes three systems: CFRP drilling, data acquisition, and experimental observation. The schematic diagram of the experimental plan and the connection of the test device are shown in FIG. 13.

[0084] According to the comparative analysis results of experiment and simulation shown in FIG. 14, the axial force and torque of the drilling multi-scale finite element model and experiment have the same value and change trend at each stage of drilling, only in the T.sub.IV and T.sub.V stages. Compared with the experiment, the reduction speed of axial force and torque is slower. The main reason is that the fiber adopts the degradation method based on fracture toughness in the multi-scale finite element model. The set dissipation energy value is obtained based on empirical evaluation. Due to the effect of the CEs in the bonding area in the simulation of the interlayer delamination phenomenon, some macro-units in the model are not deleted, which leads to too long contact time between the unit and the drill bit, which in turn generates some small forces and torques.

[0085] It can be clearly observed from FIG. 15(a) that the CFRP pore wall is mainly composed of resin-coated surface and fiber fracture surface, and the pore wall coating surface near the entrance is more obvious. The damage of CFRP from the entrance to the exit of drilling is mainly manifested as damages such as splitting at the entrance, tearing at the exit, burrs, separation between layers, radial crushing, and microcracks. According to the comparative analysis of the damage position of each layer in FIGS. 15(a) and (b), it is found that the pits and the hole wall damage positions are almost the same. Therefore, the multi-scale finite element model established in this paper can simulate the damage state of the CFRP drilling hole wall more realistically.

[0086] CFRP is prone to burrs, tearing, delamination and other processing damages at the entrance and exit during the process of making holes. Because the dagger drill is a special CFRP drill with integrated drilling and reaming, only the fiber fracture at the entrance of the hole when drilling the CFRP laminate It is relatively flat, the entrance tearing phenomenon is not obvious and there are fewer burrs, but through the finite element analysis model and test results, there are more obvious burrs, tears, and delamination damage at the exit, as shown in FIG. 16(a) . However, because the multi-scale finite element model is based on the deletion of macroscopic elements to realize the simulation of material removal, the scale effect of the element can only simulate the damage phenomenon at the macro level, and cannot achieve the quantification of burrs and tear damage. The delamination damage in CFRP holes generally only exists between a few layers of material on the side of the exit. The size of the internal damage needs to be observed with observation equipment such as ultrasonic scanning, but the delamination at the exit will follow the direction of the dagger drilling Extending to the outside of the material, there is a part of the height of the bulge in the damaged area, which only needs to be observed with a super-depth microscope, as shown in FIG. 16(a).

[0087] The delamination damage system is based on the tool diameter. The area calculated by the number of deleted elements around the hole is divided by the area of the reference hole to calculate the delamination damage coefficient. The measurement method is to export the picture in the ABAQUS software, and then import it into the CAD software. The belt tool adopts the equal proportion method to realize the measurement respectively. Since the rotation speed of the drill bit is too fast under actual working conditions, and the deletion of the multi-scale finite element model unit is irrecoverable, the area of the damaged unit in the center of the layered area is used to calculate the layered damage coefficient, as shown in FIG. 24(b) Shown. Therefore, the delamination damage coefficient can be expressed as

[00023]$\begin{array}{cc}{D}_d=\left(\frac{A_d}{A_{\mathrm{nom}}}\right)\%& \left(34\right)\end{array}$

[0088] Where D.sub.d denotes the delamination factor, A.sub.d denotes nominal diameter, and A.sub.max denotes the delamination areas. The delamination damage coefficient at the exit of the drilling CFRP under different processing parameters in the multi-scale finite element model and test conditions is calculated. It can be seen from FIG. 17 that the finite element model and the test results are relatively close to the results of the damage coefficient analysis. At the same time, the change trend of the machining parameters also tends to be consistent with the change of the machining parameters.

[0089] In summary, the multi-scale finite element model established in this paper can truly simulate the CFRP damage state under actual working conditions. The maximum errors of axial force, torque, and outlet delamination damage coefficient are 3.37%, 7.69%, and 4.28%, respectively. At the same time, it can also simulate the phenomenon of hole wall surface damage and exit delamination more realistically.