3D Progressive Damage Based Macro-Mechanical FE Simulation of Machining Unidirectional FRP Composite

Finite element (FE) simulation is a powerful tool for investigating the mechanism of machining fiber‐reinforced polymer composite (FRP). However in existing FE machining simulation works, the two‐dimensional (2D) progressive damage models only describe material behavior in plane stress, while the three‐dimensional (3D) damage models always assume an instantaneous stiffness reduction pattern. So the chip formation mechanism of FRP under machin‐ ing is not fully analyzed in general stress state. A 3D macro‐mechanical based FE simulation model was developed for the machining of unidirectional glass fiber reinforced plastic. An energy based 3D progressive damage model was proposed for damage evolution and continuous stiffness degradation. The damage model was implemented for the Hashin‐type criterion and Maximum stress criterion. The influences of the failure criterion and fracture energy dissipa‐ tion on the simulation results were studied. The simulated chip shapes, cutting forces and sub‐surface damages were verified by those obtained in the reference experiment. The simulation results also show consistency with previous 2D FE models in the reference. The proposed research provides a model for simulating FRP material behavior and the machining process in 3D stress state.


Introduction
Fiber reinforced polymer composite (FRP) material is finding increasing applications in modern aerospace, automobile and sports industries. A significant amount of machining work is necessary as material removal is required for the FRP component to meet the dimensional requirements before assembly. To understand the machining process and to avoid machining damages such as fiber breakage, debonding, or delamination, many experiments have been performed in past years [1][2][3][4].
On the other hand, finite element (FE) simulation has been widely used for investigating the machining or forming process of homogeneous and isotropic materials [5,6]. This technology is also extended to inhomogeneous and anisotropic FRP composite when material behavior is defined appropriately [7]. Most existing work in this area focused on orthogonal cutting of unidirectional FRP. Generally, these FE models are based on either macromechanical or micro-mechanical approaches. In macromechanical FE machining models, the FRP material is modeled as an equivalent homogeneous anisotropic material (EHM). In micro-mechanical models, multiphase and different constituents (matrix and fiber) in FRP are modeled separately [8][9][10]. There are also studies which combine these two approaches in one model to leverage the advantages of both methods [11,12].
Many macro-mechanical FE models are built as twodimensional, using various material damage models and failure criteria. The Tsai-hill criterion was used in the 2D FE machining models developed by Arola et al. [13,14], Mahdi et al. [15], and Mkaddem et al. [16]. In the 2D FE model developed by Lasri et al. [17], damage analysis was conducted using the Hashin, Maximum stress and Hoffman failure criteria together with the instantaneous stiffness degradation strategy. Other 2D plane stress FE He et al. Chin. J. Mech. Eng. (2018) 31:51 models [18][19][20] apply the Hashin criterion together with the progressive damage model built in Abaqus software. In these models, the predicted principle cutting force generally agreed with experimental records, but a significant difference was observed between the predicted and experimental thrust forces [14,17,18].
3D FE models for FRP machining simulation are not frequently found in the literature. Mahdi and Zhang [21] presented a 3D FE model in which the maximum shear stress was used as the material separation criterion. Rao et al. [22] developed a 3D FE model with the Tsai-hill failure criterion to simulate the cutting force and chip formation for orthogonal machining of UD-CFRP composite. Santiuste et al. [23] developed a 3D laminate model based on the three-dimensional Hoffman criterion. In these limited research models, the FRP is modeled as an elastic material with instantaneous failure behavior, and neither the chip shape nor the sub-surface damage was studied in a quantitative manner.
Recently, finite element modeling was also used in simulation of drilling or milling of FRP composite. In the meso-scale CFRP drilling FE models developed by Isbilir et al. [24] and Feito et al. [25], sudden stiffness degradation was assumed in intra-lamina failure analysis with the 3D Hashin and Hou criteria, respectively. Phadnis et al. [26] developed a meso-scale FE model for CFRP drilling simulation. The 3D Hashin criterion was used in fiber failure prediction, and the Puck criterion was used to model matrix failure. Still, a sudden degradation model was applied in intra-ply damage analysis such that the material point offers no resistance to deformation if any damage initiation condition was met.
From the existing research studies, there lacks a FE machining simulation work that analyses the 3D progressive damage model of FRP. As such, the FRP material behavior during cutting is not fully investigated in a general stress state, because the 2D damage model is only applicable with plane stress assumption. Additionally, different failure criteria were used in different FE models, but the influence of failure criterion on the simulation results was not discussed. To reveal the 3D material behavior that explains the phenomena in FRP cutting, a 3D macro-mechanical FE model was developed according to the orthogonal cutting experiment by Nayak et al. [3]. An energy based progressive damage model was proposed with the Hashin-type and Maximum stress criteria for continuous stiffness degradation of UD-GFRP. The material damage model was implemented in a user subroutine (VUMAT) and incorporated into the FE models for machining simulation.
The proposed 3D macro-mechanical FE model was verified by the reference experiment [3] in terms of chip formation, cutting forces, and sub-surface damages. In addition, the simulation observations were also consistent with the results from existing FE models in the literature. The influences of failure criterion and fracture energy dissipation on the simulation results were also studied. The proposed damage model can also be used in FRP structure analysis (e.g., under low-velocity impact load).

Assumptions, Control Volume and Boundary Conditions
The FE model is developed based on the commercially available FEM software Abaqus v6.11. The geometry and boundary conditions of the model are shown in Figure 1.
The workpiece material is UD-GFRP. Only a small region of the workpiece (2 mm × 1 mm) close to the tool tip is modeled and meshed, while enough for cutting process to reach steady state conditions. For the modeled area, the bottom is restrained for displacement in the x, y, and z directions. The displacement of the entire left side and partial right side (corresponding to the area fixed by the fixture) were also restrained in the cutting direction. For consistency with the experiment [3], a tool with a 0.05 mm edge radius, 6° clearance angle and 10° rake angle, a cutting speed of 0.5 m/min, and a depth of cut of 0.2 mm were used in the simulation. Different fiber orientations (15°, 30°, 45°, 60°, 75°, and 90°) were tested. The machining process is simulated as a quasi-static explicit analysis. The cutting tool is modeled as a rigid body with a predefined velocity v in the negative x-direction. A reference point at the top right corner is defined to control the movement of the tool and offer output of the reaction force. The UD-GFRP material is modeled as an equivalent homogeneous and orthotropic material with elastic-failure behavior, and is meshed with eightnode brick elements (C3D8R). The mesh contains 20 layers of elements through the thickness direction, and was refined in the zone surrounding tool tip with a size of approximately 5 microns. The total element number is 160000. The mesh will influence the accuracy and step increment of the simulation, but the sensitivity analysis of the numerical parameter is outside the scope of this study.
The interaction between the cutting tool and workpiece is modeled as surface-node surface contact available in Abaqus/Explicit. The contact is assumed to follow the Coulomb friction law. The experiment shows that friction increases with the fiber orientation [9], so six different friction coefficients, i.e., 0.4, 0.5, 0.6, 0.7, 0.8, and 0.9, were set for fiber orientations of 15°, 30°, 45°, 60°, 75°, and 90°, respectively in the simulation. It is noted that these coefficients are obtained from a pin-on-disc experimental test with a normal load of 120 N, and thus they are used in machining simulation as an approximate estimation that reflects the influence of fiber orientation on the friction. The thermal effect is neglected since the cutting process is conducted at a rather low speed.

Material Separation for Chip Formation
In the literature, the chip can be simulated by means of material separation criterion (node/element splitting), element deletion, or a pure deformation process based on adaptive remeshing [27][28][29]. For remeshing, the state variables have to be frequently interpolated from the old mesh to the new mesh, leading to error accumulation and deterioration of results. To avoid this problem, chip formation is achieved based on the material failure criterion and element deletion approaches, which were also applied in previous machining simulation studies [17][18][19]. When all material points in an element fail, the element is deleted and removed from the mesh. No separation line along the tool tip path is predefined as required in the node splitting technology.
It should be noted that element deletion induces loss of volume, which violates the law of continuity, and thus the mesh design is notably fine to minimize this volume loss. The issue of convergence is important for quasi-static simulation, as indicated by the energy approximation balance in cutting simulation. The energy balance is checked for each simulation run to ensure that the ratio of kinetic energy to internal energy meets the requirement.

Material Constitutive Law and Properties
The global stiffness is used to replace the individual properties of the fiber and matrix in FRP in macro-mechanical FE model. The material constitutive for undamaged material is shown Eqs. (1)−(3): Key properties of UD-GFRP are adopted from Ref. [4] and extended as shown in Table 1. These properties are

Failure Criterion
Numerous failure criteria are available to describe FRP damage, among which the 3D Hashin criterion is most widely used. Varieties of modified criteria based on the Hashin criterion have been proposed [30]. In this study, the 3D Hashin-type and Maximum stress failure criteria are used in failure prediction, both of which distinguish various failure modes. The corresponding damage initiation criteria are shown in Table 2, where the stress is also expressed in the material coordinate system. The failure mode is generally named with the word "fiber" or "matrix" although homogeneous anisotropic material is assumed.
The fiber failure mode is associated with the longitudinal direction (Dir. 1), with f 1T and f 1C denoting tension and compression failure indexes respectively [17,31]. The matrix failure is associated with the transverse (Dir. 2) and thickness (Dir. 3) directions, with failure indexes represented as f 2T , f 2C and f 3T , f 3C . The mechanism of delamination damage is ignored. This assumption is reasonable since the out-of-plane stress is not significant and plane stress was used in most orthogonal cutting FE models [14][15][16][17][18][19][20].
When the failure criterion is met (i.e., the failure index in Table 2 exceeds 1.0), the material point reaches the onset of damage and stiffness degradation is initiated. Two types of material degradation strategies, i.e. instantaneous failure and progressive failure, have been widely used [32].
For instantaneous failure, the material is assumed to fail immediately in a mode at damage initiation, and thus the material properties associated with that failure mode are degraded instantly. It's a common practice to reduce selected stiffness with respect to the failure mode to zero or a near zero value. Table 3 lists two material degradation rules that are commonly used [31][32][33]. The difference between the two rules lies in how the fiber failure influences the matrix failure mode. In Rule 1, all of the moduli are set to 0 when fiber failure is detected, whereas in Rule 2, the transverse modulus remains unchanged at fiber failure. In both rules, the transverse  Fiber failure For the progressive damage model, the material property is gradually degraded. Various approaches address progressive degradation using continuum damage mechanics (CDM) theory [34][35][36]. The CDM replaces the mechanical properties of the damaged materials with those of the homogenous materials by associating the damage mechanisms with their effects on the elastic constants of the materials. Several energy based stiffness degradation models are presented by Fang et al. [37], and Lapczyk et al. [38], while this study extends previous studies [38,39]. For each failure mode, the damage evolution in the post-damage initiation phase is indicated by a damage variable.

Damaged Material Response
The stiffness coefficient of the material is reduced when the effect of damage is considered. The principal damage variables d i ∈ [0, 1] is used to express the damage states in direction i, where 0 means undamaged and 1 means completely damaged. Because two failure modes are associated with one direction, the final damage variable d i in each direction is expressed as: where d iT and d iC indicate the damage variables for the tension and compression modes, respectively, in direction i.
The effective stress σ was introduced in the progressive damage model. The relationship between the effective stress σ and nominal stress σ is where in which d 1 , d 2 , and d 3 are the damage variables for direction 1, 2, and 3 as defined in Eq. (4), and d 12 , d 13 , and d 23 are the damage variables for shear failure modes in the 1-2, 1-3, and 2-3 planes. The damaged stiffness matrix C d is expressed in matrix form using the damage variables and the undamaged stiffness matrix, as shown in Eqs. (6)−(8): The constitutive relationships for the damaged material can be obtained as follows:

Damage Evolution
The evolution of the damage variable is based on the fracture energy dissipated during the damage process. The crack band model is adopted to alleviate the mesh dependence of the energy dissipated in the finite element result [37,38]. The fracture energy is preserved and is the function of stress, the characteristic element length and the strain at strength limit.
The damage evolution is expressed as the stress-displacement relation in Figure 2 for each failure mode. The positive slope corresponds to the linear elastic material behavior prior to damage initiation, and the negative slope corresponds to damage accumulation till failure. After damage initiation, the damage variable evolves from 0 (damage initiation) to 1 (completely damaged) based on equivalent displacement δ eq . In Figure 2, δ 0 eq is the equivalent displacement at which the damage initiation is met and δ f eq is the equivalent displacement at which the material is completely damaged. The area of the triangle 0AC corresponds to the energy dissipation G for that mode.
(9) σ = C d ε. The equivalent displacement and equivalent stress are calculated by the following expressions, where L c represents the characteristic element length, and operator �x� = (x + |x|)/2.

For shear modes:
The damage variable is subsequently expressed as: For the Hashin-type criterion, it is assumed that the damage variables for shear failure, i.e., d 12 , d 13  For the Maximum stress criterion which includes individual shear failure modes, the shear damage variable evolves on its own and is calculated in a manner similar to that of the fiber or matrix failure modes. However another consideration is that shear damage is also influenced by fiber failure and matrix failure, as indicated in Table 3. Therefore, the final shear damage variables d 12 , d 13 , and d 23 for the Maximum stress criterion are computed as follows: In Eq. (20), δ If eq is the completely damaged equivalent displacement of failure mode I: σ k eq = L c (�σ k ��ε k �) δ k eq , k = 12, 13, 23.  The equivalent displacement and stress at damage initiation, δ I0 eq and σ I0 eq , are computed as: where f I is the failure index for mode I in Table 2.
The exact fracture energy for progressive damage is difficult to obtain. Due to the absence of experimental data, the energies for various modes were assumed for the numerical computations. The energies for fiber tension, fiber compression, matrix tension and matrix compression modes were specified as 8 N/mm, 8 N/mm, 1 N/ mm, and 2 N/mm, respectively. For the Maximum stress criterion, an additional value of 0.5 N/mm was assumed for each shear damage mode. Moreover, other energies were also specified to investigate the influence of the fracture energy level on the simulation results. In future, more accurate degradation parameters could be obtained through experiments.

Subroutine Implementation
The Abaqus user material subroutine (VUMAT) implementation of progressive degradation is illustrated in Figure 3. The program first reads the input values of material modulus and strength. For each material point, the program obtains the damage variable of each failure mode, and calculates the damaged stiffness matrix, strain and stress using Eq. (6). The equivalent displacement and stress are calculated using Eqs. (10)−(19) before the damage initiation criterion in Table 2 is evaluated for each failure mode. If damage is initiated for one failure mode, the initiation and complete damage equivalent stress and equivalent displacement corresponding to that mode are computed. Finally, the damage variables for the current increment are calculated by Eq. (20), and new stress components are modified if necessary.
All damage initiation indexes and damage variables are stored as solution dependent state variables (SDSV) in VUMAT. To avoid numerical difficulties, the damage variable was limited to a maximum value of 0.999 for each failure mode. In explicit analysis, a material point fails when it reaches either the fiber tensile or fiber compression mode. As such, the material point is flagged as failed and is removed from the mesh when d 1 reaches the critical value of 0.999. 5 show an example of progressive damage analysis of chip formation using the Hashin-type and Maximum stress criteria for 45° fiber orientation. In Figure 4, the matrix damage contour in Dir. 2 is severe, whereas longitudinal (fiber) damage and the matrix damage in Dir. 3 are minor. Transversal damages in Dir. 2 initiate at the tool-workpiece contact point and progress in the direction parallel to fiber till chip is formed. Matrix cracking in Dir. 2 also extends vertically to form sub-surface damage. The chip formation simulation using the Maximum stress criterion is demonstrated in Figure 5. For the same   failure mode, the damage generally initiates later because the material is less likely to reach the damaged state when evaluated with the Maximum stress criterion. Among all failure modes, the fiber damage in Dir. 1, the matrix damage in Dir. 3, and shear damage in the 1-3 plane are rarely observed. The transversal matrix damages in Dir. 2 are more obvious, but are limited to the vicinity of the cutting tool. Matrix cracking damage exists below the tool flank and extends vertically. The matrix crushing damage initiates near the tool tip-workpiece contact area and propagates in the vicinity of the tool edge. Both matrix damage patterns are different from those in the Hashintype model in which the damage progresses along fiber direction, because the shear stress does not contribute to transversal damage initiation in the Maximum stress criterion. The most obvious damage is shear damage in the 1-2 plane which propagates along fiber direction to the top surface of the workpiece due to the nature of the shear damage initiation criterion and its evolution mechanism. The damage extension pattern implies that the shear damage in the 2-3 plane is primarily induced by transversal matrix cracking and crushing.

Figures 4 and
The final simulated chip shapes are shown in Figure 6. The two chips look similar, while close examination shows that the chip predicted by the Maximum stress criterion has a larger size than that predicted by the Hashintype criterion. This result is attributed to the difference in the damage initiation expression and in the equivalent displacement definition. It usually requires larger deformations for the material to reach damage initiation and full damage state in the Maximum stress model such that a larger chip is produced in the numerical simulation. The chip mechanism [1] and the simulated chip shapes from 2D models by Lasri et al. [17], and Santiuste et al. [18] are also demonstrated in Figure 6. Although a detailed figure of the chip formation in referred experiment is not available, both chips bear similarities to the chip mechanism and also to those in previous studies.
The simulated chips for different fiber orientations are shown in Figure 7. For the same fiber orientation, two models predict similar chip shapes, and those predicted by the Maximum stress criterion always display a larger chip thickness or chip size. In both models, it is also obvious that the simulated chip size decreases with the increase in fiber orientation. This result is consistent with the experiment [1,3] and also with the previous simulation results [17]. The reason for this result is that the length of the uncut fiber above the machined surface is reduced for increased fiber orientation [1].
Different levels of the fracture energy value were specified in the simulation to investigate the influence of fracture energy on the simulated chip shape. The results are shown in Figure 8. For both criteria, the chip shape does not display obvious changes at the three energy levels. Only the chip predicted using 0.1 times the energy level is slightly smaller. A similar phenomenon was observed in the 2D FE simulation model by Soldani et al. [19]. This is because the damage progresses with material deformation till the specified fracture energy is dissipated. When the energy is small, the material is assumed to reach complete damage at a small deformation for a failure mode, and thus a slim chip is obtained. For the same reason, it could also be imagined that the damage tends to remain more localized when a smaller energy value is specified, as demonstrated in Section 4.3. However, the chips predicted using the default energy are almost identical to those predicted with 5 times the energy. This result might indicate that when the fracture energy reaches a critical value, further increases in the energy value do not lead to obvious change in the simulated chip shapes.

Cutting Forces
An example of the cutting force is shown in Figure 9 for the 30° fiber orientation. The cutting forces predicted by the two criteria are significantly different at the same energy value. To maintain consistency with the experiment, the cutting and thrust forces were found by averaging the force per unit width (mm) during the quasi steady-state region of the force vs. time plot [3]. For the 30° fiber orientation, the Maximum stress criterion predicts a higher principal cutting force (33 N) than the Hashin-type model (13 N). This result can be explained by their different failure initiation expressions in Table 2. There is a chance that the damage initiation criterion is not yet met when evaluated by the maximum stress criterion, but is met when evaluated by the Hashin-type criterion for the same stress state. In the former case, the material continues elastic deformation such that the reaction force continues to improve. The difference in the equivalent displacement expression also contributes. With the same energy value, the material takes on additional deformation before reaching complete damage since the equivalent displacement calculated in the Maximum stress model is always smaller. A higher simulated cutting force value is thus obtained. The experimental cutting force for this orientation is approximately 30 N per unit width [3], which is close to the value predicted by the maximum stress criterion.
Oscillation in the value of the cutting forces occurs in each model due to material stiffness reduction upon failure. Additionally, it is observed that the cutting force simulation reaches steady state earlier in the Hashin-type model than that in the maximum stress model. The thrust forces predicted by the two criteria are also different. Both thrust forces remain negative at the early stage but increase gradually and turn into positive, reaching nearly 0.5 N and 3.5 N, respectively, at a later stage. In the beginning, the thrust force is negative because at the initial stage of penetration of the tool edge into workpiece, the workpiece contacts the tool edge and the rake area. Besides the reaction forces ahead in the horizontal direction, the tool also encounters a reaction force from upper left perpendicular to the tool rake. The resultant force has a negative (downwards) thrust force component. However, with additional advances of the tool, when the contacts between the workpiece and the bottom of tool edge and the tool flank are established, the reaction force from the tool bottom becomes dominant, which results in a positive thrust force.
The experimental thrust force is approximately 32 N per unit width for this condition. Both models predict far smaller thrust forces than the experiment, as is consistent with the 2D FE models in Refs. [9,17,18] in which significant differences are observed between the predicted and experimental thrust force results. According to Wang and Zhang [2], this difference can be explained by the "bouncing back" mechanism of the material under the tool flank. During cutting, the material under the tool tip is pushed down and springs back partially elastically after the tool has passed. The real depth of cut is different from the nominal one, and thus the experimental thrust force is strongly influenced.
The influence of the fracture energy value on the cutting forces simulation is demonstrated in Figure 10 and Figure 11. In both models, the increase in energy from the 0.1 times value to the default value leads to an increase in both the cutting force and thrust force, because higher energy allows larger deformation before the material point reaches complete damage. Once more, when the energy reaches a critical value, the simulated cutting and thrust force magnitude and the evolving pattern do not show an obvious change if the fracture energy increases further thereafter.
For different fiber orientations, the simulated cutting forces are compared with the experimental data in Figure 12. From Figure 12(a), for all fiber orientations, the Maximum stress model predicts a higher principal cutting force than the Hashin-type model and is closer to the experimental data. This result is reasonable considering the nature of these two criteria. The experimental cutting force increases with increasing fiber orientation. Beyond 30° orientation, both models predict an increasing tendency with increasing fiber orientation. However, for fiber orientations smaller than 30°, the predicted trends are not strictly consistent with the experiment, as is similar to the result of 2D FE model by Lasri et al. [17] in which the minimum simulated cutting forces occur in the 30° fiber orientation. The discrepancy is related to the limitation of the macro mechanical FE model in which the reaction force is obtained based on assumed material damage behavior throughout the chip formation process. However in real cutting, the chip formation process is composed of events of material fracture and chip slipping due to the pushing effect of the tool rake. The smaller the fiber orientation, the easier for chips to slip along shear plane or the fiber/matrix bonding interface [1], and thus a smaller experimental principal cutting force is recorded for lower fiber orientations.
As shown in Figure 12(b), the thrust forces predicted by both models are smaller than experiment for all fiber orientations, due to the bouncing mechanism explained previously. The same phenomenon was observed in the 2D FE model in Refs. [9,17,18]. However, the trends in both models agree with the experimental observations.

Sub-surface Damages
In VUMAT, the damage variables for different failure modes are stored in different SDSVs and updated in each increment. The progression of various sub-surface damage modes can also be observed in Figure 4 and Figure 5. Damage occurs where the corresponding damage variable is larger than 0. For 45° fiber orientations, the subsurface damage in the Hashin-type criterion model is characterized by transversal matrix damage in Dir. 2 and subsequent shear damages in 1-2 and 2-3 planes. In the Maximum stress model, the sub-surface damage is attributed to shear damage in the 1-2 plane, the transversal matrix damage and subsequent shear damage in 2-3 plane.
Since a material point is associated with nine possible failure modes (d 1T , d 1C , d 2T , d 2C , d 3T , d 3C , d 12 , d 13 , and d 23 ) in the simulation, the maximum value among these nine damage variables throughout the cutting process was recorded as simulated sub-surface damage value. The damage variable only increases during the cutting process since the damage is irreversible. With this method, the predicted sub-surface damage contours for the final chip are shown in Figure 13, with the values of approximately 178 microns in the Hashin-type model and 330 microns in the Maximum stress model when the default energy is specified. Larger sub-surface damage was predicted by the Maximum stress criterion due to exactly the same reasons of its thicker chip and higher cutting force. This result also agrees with the fact that higher cutting forces always imply more damage. Figure 13 also shows the influence of the energy level on sub-surface damages, with the values compared in Figure 14. An increase in the energy value leads to increasing simulated sub-surface damages for both criteria. As discussed in Section 4.1, with a higher energy value, the material point experiences larger deformation before reaching the specified fracture energy value, and thus the damage tends to extend to a larger (or deeper) area. Again, when the energy reaches a critical value, the sub-surface damage tends to remain stable.
The variation of simulated sub-surface damage with fiber orientation is compared with that obtained in experiment in Figure 15. The experimental sub-surface damage is measured by the extent of the spread of fluorescent dye along the vertical axis from the trimmed edge, with the values showing an increasing tendency with fiber orientation [3]. For each fiber orientation, the maximum stress model always predicts a larger sub-surface damage.
In both models, the simulated sub-surface damages increase with increasing fiber orientation, and the values and trends are close to the experiment for fiber orientations below 75°. However, the predicted damage drops when the fiber orientation increases further to 90°, leading to a large discrepancy between the simulated and experimental sub-surface damage at larger fiber orientations. This result occurs because the simulated sub-surface damage is only related to the material failure (fiber/ matrix cracking or crushing). However during the cutting experiment, fiber-matrix interface debonding also contributes to sub-surface damages. This mechanism plays a minor role for small fiber orientations but a significant role at 90° orientation, where the interface debonding can extend severely below the machined surface [1][2][3]. The fiber-matrix interface and its damage process are not modeled and simulated appropriately in the macromechanical FE models, and a multiphase micro-mechanical model is thus required.

Conclusions
(1) The failure criterion has a strong influence on the simulation results, especially on the cutting force. The Maximum stress criterion always predicts thicker chip shapes, higher cutting forces, and larger sub-surface damages due to the difference in expression for failure initiation and damage evolution.
(2) For both criteria, an increase in the fracture energy level results in increases in chip thickness, cutting and thrust forces, as well as sub-surface damage predictions. However, when the energy reaches a critical value, the simulation results tend to remain stable, and no obvious change is observed if the energy is increased further. (3) For the assumed fracture energy value, the chip shapes simulated by both FE models are similar to the mechanisms. The simulated chip length decreases with increasing fiber orientation. The trends of cutting force and thrust forces with respect to fiber orientations agree reasonably with the experiment, with the cutting force increasing with fiber orientation, and the thrust force increasing initially and subsequently decreasing. Consistent with other 2D macro-mechanical FE models in the literature, the simulated thrust force has an order of magnitude smaller than the experiment. (4) The sub-surface damages predicted by both models show an increasing tendency with fiber orientation and agree with the experimental results for fiber orientations up to 75°. A significant difference exists between the predicted and experimental results at higher fiber orientation (above 75°) since the macro-mechanical model does not simulate the fiber/matrix interface debonding behavior, revealing one of the limitations of the macro-mechanical FE model.