2D macro-mechanical FE simulations for machining unidirectional FRP composite: the influence of damage models

2.1 Assumptions, control volumes and meshes

The geometry and boundary conditions of the FE model are shown in Figure 1. Only a small material area (3 mm ×2 mm) close to the tool tip is modeled and meshed, as 3 mm length of cut is enough to reach a steady machining process.

Figure 1:

Orthogonal machining of fiber-reinforced polymer composites (FRP).

Cutting conditions are set as the same as in the experiment [5]. The geometry of the cutting tool (edge radius, rake and clearance angle), fiber orientation and cutting parameters such as cutting speed, depth of cut, etc, are listed in Table 1.

The machining process is simulated with quasi-static explicit analysis. The explicit dynamic method is used due to the advantages of computational efficiency for applications involving large deformation and high nonlinearities such as machining. 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, offering output of the reaction force by summing the contact forces. For the modeled area, the bottom and both sides of the workpiece are restricted as shown in Figure 1.

The FRP material is modeled as equivalent homogeneous anisotropic and elastic-failure material. For 2D cutting simulation, the workpiece is meshed with plane stress, quadrilateral, linearly interpolated element with reduced integration and hourglass control (CPS4R). The thermal effect is neglected as the cutting is conducted at a very low speed. For this study, the mesh size at the tool tip is around 4 μm. There are, in total, 13,381 elements in the mesh.

Interaction between the cutting tool and the workpiece is modeled as surface-node surface contact available in ABAQUS/Explicit. The contact is assumed to follow Coulomb friction law. The friction is difficult to obtain since the interaction between the tool flank and the workpiece and that between the tool rake and chip is a complex phenomenon. For simplicity, a constant friction coefficient of 0.5 is used during the cutting process.

2.2 Material constitutive law and properties

A global macro-stiffness was used to replace individual properties of fiber and matrix in FRP composites. The material constitutive is shown in Eq. (1). These stress-strain relations are associated with the principal directions (i.e., material coordinate direction in 1, 2, 3 axes). For off-axes load, the stress and strain should be transformed to the fiber direction with the transformation matrix. The properties of the unidirectional FRP were adopted from [8] and are listed in Table 2, which were also used in [17], [18], [19]. The experiment shows that strain rate has a strong effect on strength of FRP [22]. The mechanical properties from the standard test may not reflect exactly the material behavior, since the cutting speed is higher than that used in the tensile test [17]. However, the simulation result is not much influenced as the cutting seed is considerably low when compared to the impact test [22].

2.3 Failure criterion

In this work, four failure theories are studied, the damage initiation criteria of which are listed in Table 3. Hashin and Maximum stress criteria consider different failure modes (such as fiber tension, fiber compression, matrix cracking and matrix crushing) [23]. Note that the stress in the table is also expressed in the material coordinate system.

2.4 Damage evolution and stiffness degradation

The degradation of material stiffness initiates when the failure criterion is met, i.e., the failure index in Table 3 exceeds 1.0. Two stiffness degradation strategies are studied.

(1) Instantaneous damage: The material is assumed to fail (or fail a failure mode) at the onset of damage initiation. It is a common practice to degrade associated elastic properties instantly to zero according to the failure mode. Instantaneous damage is applied to each failure criterion. In the following sections, the FEM models applying instantaneous damage are referred to with corresponding criterion name.

Table 4 shows the degradation rules with respect to stiffness reduction for this damage strategy, where the coefficient & was set as 0. The selective stiffness reduction scheme is applied to Maximum stress and Hashin theories. It is implied that the matrix failure will induce shear failure, while fiber failure will induce both matrix and shear failures. Since Hoffman and Tsai-Hill theories are unable to distinguish failure modes, a non-selective stiffness reduction scheme is used.

In some literatures, the material properties were degraded by setting a small coefficient value, like 0.05 [18]. However, the existence of residual stiffness will lead to a continuously increasing stress calculation and thus a continuously increasing cutting force (Figure 2), as violating experimental observation. So the coefficient of 0 was used in this study.

Figure 2:

Cutting and thrust forces for instantaneous failure with 1% residual stiffness (&=0.01, fiber orientation 45°).

For Maximum stress and Hashin criteria, the status of material is evaluated based on the modes of failure the material has reached. A material point will fail only when fiber failure modes (either tensile or compression) take place, i.e., f_ft or f_fc reaches 1.0.

(2) Progressive degradation: Progressive degradation is applied to Hashin and Maximum stress criteria. The stiffness is degraded gradually rather than instantly to zero. A damage variable is defined for each failure mode to indicate the damage accumulation and evolution [24]. Damage variables for fiber tension/compression, matrix tension/compression and shear failure modes are denoted by d_ft/d_fc, d_mt/d_mc and d_s, respectively.

The constitutive law of damaged material is expressed in the stress-displacement relation shown in Figure 3 for each failure mode. The positive slope corresponds to linear elastic response prior to damage initiation, and the negative slope post-damage behavior. After damage initiation, the damage variable is calculated by Eq. (2) based on equivalent displacement δ_eq, where δeq0 is the equivalent displacement at which the initiation criterion was met and δeqf is the equivalent displacement at which the material is completely damaged.

Figure 3:

Damage evolution for progressive failure.

The values δeqf for various modes depend on the respective energy dissipation G^c value, which corresponds to the area of the triangle OAC in Figure 3. By default, the material will reach complete damage for a failure mode when the corresponding damage variable reaches 1.0.

The equivalent stress and displacement for the Hashin criterion is calculated using the method presented in [24]. It is further adapted and extended in this paper for Maximum stress criterion, as shown in Table 5. L^c represents characteristic length and operator < > is Macauley bracket defined as <a>=(a+|a|)/2.

The energy dissipation for fiber tension, fiber compression, matrix tension and matrix compression modes are specified in Table 6 according to literatures [21], [24]. In addition, the energy for interface shear failure mode must be specified for Maximum stress criterion. A value of 0.05 N/mm is used. It is noted that both criteria use the same energy value while their equivalent stress/displacement are calculated using different expressions. This may also be one of the sources leading to their different simulation results.

The response of the damaged material is governed by Eqs. (3) to (6). Unlike the Maximum stress criterion where its d_s evolves independently by Eq. (2), the Hashin criterion does not include shear failure mode. That explains why the damage variable d_s for the Hashin criterion is calculated based on damage variables of fiber failure and matrix failure as in Eq. (7).

To keep consistency with instantaneous damage models, an additional control logic is also applied in the VUMAT implementation, so that matrix failure will induce shear failure, and fiber failure will induce both matrix and shear failures. A material point will fail only when reaching either fiber tensile or fiber compression failure modes. To avoid numerical difficulties, all damage variables were limited to an upper limit of 0.99, indicating the material point will fail when either d_ft or d_fc reaches 0.99.

In progressive damage models, the stress for damage initiation in Table 3 should be replaced by effective stress tensor σ^ to reflect the influence of damage on the initiation criterion. The relation between σ^ and the nominal stress in expressed in Eq. (8):

2.5 Material separation for chip formation

With the advance of cutting tools, the unwanted material is separated and a new workpiece surface formed. Chip is usually simulated as crack initiation and propagation in the Lagrangian formulation. The material separation technology could be geometrical or physical [25], [26]. In this study, material failure criterion along with element deletion approach [26], [27] was used. When all material points in an element failed, the element was deleted and removed from the mesh. No separation line along the tool tip path is predefined as required in node debond/splitting technology [25], [26].

However, the element deletion will induce loss of volume which violates the law of continuity, so the mesh is designed to be very fine to minimize this volume loss. The arbitrary Lagrangian- Euleiran [27] model is avoided because the state variable has to be interpolated from old mesh to the new mesh frequently. Errors will occur and accumulate during remeshing and the result will be deteriorated.

The material damage models are coded in ABAQUS user material subroutine (VUMAT). Convergence is an important issue for quastistatic cutting simulation. Energy balance is checked for each simulation run to ensure that the ratio of kinetic energy to internal energy meets the requirement.

3.1 Chip shape

Figure 4 demonstrates the influence of damage model on the simulated chip shapes. Hoffman and Tsai-Hill criteria predict powdered chips, because elements at the tool chip contact are deleted throughout the cutting process after internal material points failed. Continuous chips were obtained for Maximum stress and Hashin criteria models since few elements were deleted. The introduction of progressive damage will change the chips slightly, but still continuous chips were obtained.

Figure 4:

Simulated chips (t=5 ms).

A close look shows that the chip predicted by the Maximum stress criterion is thicker than that by the Hashin criterion, either in instantaneous damage or progressive damage models. This is attributed to their different shear damage mechanisms. At the early cutting process when fiber failure does not yet happen, the shear damage for the Hashin criterion is dependent on matrix cracking or crushing damages. The material point will fail shear failure mode when reaching either matrix failure mode. However, for the Maximum stress criterion, the shear damage is evaluated more specifically on the ratio of shear stress against shear strength (in the case of instantaneous damage) or on shear equivalent displacement (in the case of progressive damage). A larger zone of material was affected and finally turned into chips. This can be observed through the damage evolution in chip formation process shown in Figures 5 and 6.

Figure 5:

Damages during chip formation for Hashin-progressive damage (fiber orientation 45°).

Figure 6:

Damages during chip formation for Maximum stress-progressive damage (fiber orientation 45°).

In the case of progressive damage models in Figures 5 and 6, the shear damage is always the first type of damage developed and has the largest area among all failure modes. It initiated and progressed parallel to the fiber until complete debonding was achieved. More material is involved in shear damage for the Maximum stress criterion than for the Hashin criterion. However, the matrix damage is limited to the small zone close to the tool tip in the former, while extended to a larger area along the fiber orientation in the latter.

The chip mechanisms are also illustrated in Figure 4. It is observed that the continuous chips in both progressive damage models bear similarities to the chip formation mechanism. However, according to cutting experiments [4], [5], chips produced for 0° fiber orientation are crushed beams ruptured perpendicular to fibers. For 45° and 90° orientation, the chips are friable ribbon and fine particles, respectively. Neither model can predict, for all fiber orientations, exact chip shapes obtained in experiment, since macro-mechanical FEM models could not simulate the behavior of individual phases of fiber, matrix, interface, as well as the interactions between phases, such as chip slipping along a shear plane [2], [3], [4], [5].

3.2 Cutting and thrust forces

The cutting forces predicted by different damage models are significantly different. A typical variation of a cutting force with time is shown in Figure 7.

Figure 7:

Cutting force and thrust force (fiber orientation 30°).

In Figure 7, the Maximum stress criterion with progressive damage model predicted the largest cutting force magnitude, followed by the Hashin criterion with progressive damage model, then the Maximum stress model, then the Hashin model. Cutting forces in progressive damage models are higher than instantaneous damage model for the same criterion. This is because the material shows a higher average stiffness as the engineering modulus is reduced gradually to zero in the former, a larger reaction force thus being obtained.

The Maximum stress model always predicts a higher force than the Hashin model. This is explained by the difference in damage initiation expressions in Table 3. It is possible that material will remain undamaged and continue elastic deformation when evaluated by the Maximum stress criterion, but become damaged and degraded when evaluated by the Hashin criterion. The difference in equivalent displacement expression also contributes in the case of the progressive degradation model. With identical energy values used for both criteria, the equivalent displacement calculated for Maximum stress is also smaller, meaning that it will take more deformation for material to be completely damaged. A higher simulated cutting force value is thus generated.

Cutting forces predicted by Hoffman and Tsai-Hill criteria have the smallest values. This is because their damage initiation criteria are the easiest to meet. Also, according to the stiffness degradation rule, all the engineering constants are degraded to zero simultaneously, unlike in the case of Maximum stress and Hashin criteria where only the engineering modulus specific to a failure mode is reduced.

Oscillation in the value of cutting forces occurs in each model due to stiffness reduction upon failure. For Hoffman and Tsai-Hill criteria, the cutting forces repeat the pattern of reaching a peak value and then dropping to zero. This is because when elements are deleted, a gap is generated between the tool rake and the workpiece, therefore, no cutting force is recorded before the contact between them is reestablished.

For the cutting conditions in Figure 7, the experimental cutting force and thrust force are about 30 N/unit width [5]. All force predictions are less than experiment values. This is because factors such as tool wear and materials bouncing under tool relief are ignored in FEM. For different fiber orientations, the force predicted by Maximum stress with progressive damage model is always close to experiments.

The cutting force with respect to different fiber orientations is shown in Figure 8. When the trends with different fiber orientations are considered, the force predicted by progressive damage is relatively consistent with experimental observations, for either cutting or thrust force. This result is better than the macro-mechanical models in literature [8], [18] in which the trends of thrust force are far from the experiment.

Figure 8:

Influence of damage models on average cutting force.

3.3 Sub-surface damage

The sub-surface damage for the first complete chip is shown in Figure 9. In instantaneous degradation models, damage for a failure mode occurs where initiation index is ≥1.0, however, in progressive damage models, damage occurs where the corresponding damage variable is >0.

Figure 9:

Sub-surface damage simulation (fiber orientation 90°, t=5 ms).

All damage initiation indexes and damage variables are stored as solution dependent state variables in VUMAT. It is noted that according to the degradation rules for Hashin and Maximum stress criteria, the sub-surface damage for one failure mode could be caused by other modes. For example, in instantaneous damage model, if the matrix crushing initiation index for a material point reaches 1.0, the material will also get damaged in shear failure mode automatically.

The sub-surface damage is measured as the damaged area extended below the plane of the machined surface. For the fiber orientation in Figure 9, Hoffman and Tsai-Hill criteria show similar sub-surface damage, which corresponds to a few layers of elements removed from the mesh. In other damage models, fiber damage is rarely observed, while matrix cracking, matrix crushing and shear damage is obvious. Since both fiber and matrix failure will induce shear failure, sub-surface shear damage prediction is always the deepest. So, in this study, the simulated sub-surface damage is determined based on the maximum shear damage at material points throughout the cutting process.

It is interesting to note that sub-surface matrix crushing/cracking or shear damage contours in progressive degradation models extends to a larger area than in corresponding instantaneous damage models. This is because in the former, damage tends to progress with material deformation until specified fracture energies are dissipated, while in the latter, material is assumed to fail a mode at damage initiation so the damage tends to be more localized.

When a progressive damage model is used, the sub-surface matrix cracking and crushing extend down-left to workpiece for Hashin criterion, but are constrained to the bottom of tool flank and front of tool rake for the Maximum stress criterion. However, the sub-surface shear damage in the former model is limited to a smaller area than in the latter, due to a different shear damage mechanism. A similar phenomenon is observed for other fiber orientations, e.g., 45° orientation, when we compare Figures 5 and 6.

In the experiment, the sub-surface damage [5] is measured by the extent of spread of fluorescent dye along the vertical axis from the trimmed edge, and damage value increases with fiber orientation. The sub-surface damage simulations are compared with experiments in Figure 10. The predicted sub-surface damages in all FEM models are less than experimental observations, just like the case of cutting forces. Maximum stress with progressive damage model predicts the largest sub-surface damage for all fiber orientations, while Hoffman and Tsai-Hill criteria always predict the least. Smaller forces might imply smaller damages, but the simulated sub-surface damage does not strictly follow the order of simulated cutting force magnitude. For example, the Hashin criterion with progressive damage model predicts a larger force than the Maximum stress model, but the sub-surface damage by the former is always less than that by the latter except at 45° fiber orientation.

Figure 10:

Sub-surface damage versus fiber orientation.

In Figure 10, the simulated sub-surface damages in instantaneous damage models do not change apparently with fiber orientations, while they vary with fiber orientations in progressive damage models. The damage values and the trends with respect to fiber orientation agree well with experiments as long as the fiber orientation is ≤45°, but agree poorly when fiber orientations gets larger. This is because the sub-surface damage mechanism changes with fiber orientations in orthogonal cutting. When the orientation angle is small, the sub-surface damage is more related to fiber/matrix cracking or crushing. However, for large fiber orientation angles, the fiber-matrix interface debonding will extend severely below the machined surface and become the dominant factor of sub-surface damage [4], [5]. However, the phase of fiber-matrix interface and its damage process can not be modeled and simulated appropriately in macro-mechanical FEM models.