*2.1. Model Configuration*

A two-dimensional orthogonal, plane strain cutting model is configured in ABAQUS (version 6.14-2), as illustrated in Figure 1. There are two main sections in this model. The top section (named the chip layer) is the ECZ domain where CZ elements are embedded all around the main elements. The bottom section is the regular finite element domain without CZ elements. This configuration saves computational time compared to a fully embedded CZ model since the bottom layer is not directly involved in the tool–workpiece interaction. Instead, the bottom layer provides compliance to the material during cutting. To ensure the stress continuity, the nodes on both sides of the interface must be merged or tied with all degrees of freedom. For this reason, the mesh sizes on both sections must match to have a perfect node-to-node alignment.

**Figure 1.** Schematic of embedded cohesive zone–finite element method (ECZ–FEM) model configuration, boundary conditions, and element arrangements.

Table 1 shows the actual model dimensions used for two depths of cut (DOC), 0.1 mm and 0.3 mm. The boundary of the bottom layer is fixed in both translational directions (X and Y). The element size, *d*, is set at 0.01 mm. The bottom layer is meshed structurally with brick elements, while for the top layer, the elements are tilted by 45◦. The inclined elements are necessary because the maximum shear stress to fracture the material is expected to be around 45◦ based on the Merchant's Circle [15]. This configuration can avoid numerical instability when no CZ mesh aligns with the preferred fracture direction. The main elements are the four-node plane strain elements CPE4R, and the CZ elements are the four-node two-dimensional cohesive elements COH2D4. To embed zero-thickness CZ elements, all elements and nodes of the chip layer need to be assigned through the input file directly because each CZ element shares nodes with adjacent two main elements, as shown in Figure 1. The CZ element is defined by nodes ABCD, in which A and D belong to the element on the left side (identical to Nodes 1 and 4), while B and C belong to the right side (identical to Nodes 2 and 3). Since these two pairs of nodes are overlaid geometrically, they cannot be identified from the graphic user interface. The meshing process is automatized by a separate MATLAB code.


**Table 1.** The model dimensions and depths of cut (DOC) used.

A complete mesh is imported to ABAQUS/EXPLICIT to set up other boundary conditions. The plane strain thickness of 3 mm is also applied to the model to be consistent with the thickness of the actual sample. The cutting tool is modeled as a rigid body with a constant speed at 10 m/min to match with the experiment. The tool has a rake angle of zero, a clearance angle of 7◦, and an edge radius of 11 μm.

#### *2.2. Damage Criteria*

To apply the ECZ–FEM to a brittle cutting process, the material properties of the main and CZ elements and their damage criteria are defined separately despite being within the same entity. Assuming an isotropic, brittle material, the main element is defined by the modulus of elasticity (*E*), Poisson's ratio (μ), the ultimate strength (σu), and damage criteria of the material. Although the model

is meant to impart fracture-based failure on the CZ mesh, the main element should still allow failing to avoid excessive element distortion when no fracture occurs. For this reason, the damage to the main elements is defined by an initiation at the ultimate strength followed by a progressive damage evolution by the Hillerborg's fracture energy theory. The total energy required to completely degrade the element after the damage initiation is *Gf*, which can be calculated from the material's fracture toughness *Kc* by Equation (1):

$$
\mathcal{G}\_f = \left(\frac{1-\nu^2}{E}\right) \mathcal{K}\_c^2. \tag{1}
$$

The degradation is in a linear manner [16], such that

$$D = \frac{\overline{u}}{\overline{u}\_f} \, \tag{2}$$

where *u* is the equivalent element displacement after the damage initiation; *uf* represents the equivalent displacement at failure. The displacement at failure is calculated by

$$
\overline{u}\_f = \frac{2G\_f}{\sigma\_u},
\tag{3}
$$

where σ*<sup>u</sup>* represents the ultimate stress. These are standard steps to simulate material failure for metal cutting [16]. It should be emphasized that this damage definition for the main element is to ensure the model stability by avoiding excessive element distortion.

The properties associated with the CZ elements embedded in the chip layer are defined differently. The cohesive zone is a mathematical approach in which the work is done to overcome the energy needed to open a crack. This work can be described by a traction–displacement relationship, *t-*δ*,* as seen in Figure 2. Damage initiation is related to the interfacial strength (i.e., the maximum traction *tc*) on the traction–displacement relation, and the area under the relation represents the fracture energy, *Gf*, as defined in Equation (1).

**Figure 2.** Bilinear traction–displacement (*t*-δ) model for the cohesive element.

In this study, a bilinear traction–separation law is adopted along with the mixed-mode progressive damage. The maximum traction *tc* should be equal or less than the strength of the material to be able to fail, while a lower strength can improve the convergence rate of the solution. In general, the variations of the maximum strength do not have a strong influence on the results [12]. Hence, the 80% ultimate stress is selected here. The initial stiffness *k* should be large enough to ensure the continuum between the two adjacent bulk elements, but small enough to avoid numerical issues such as spurious oscillations of the tractions in an element. Studies suggest that the initial stiffness of CZ elements can be calculated from Equation (4), which balances accuracy and simulation stability [12,17,18].

$$k = a \frac{E}{d'} \tag{4}$$

where *E* is the bulk elasticity, *d* is the maximum element size, and α is taken as 1.

The maximum deflection of a CZ element δ*<sup>c</sup>* is determined by given *Gf* and *t*c, as shown in Figure 2. Therefore, the deflection can become relatively large compared to the element size when a fine mesh is used. A large deflection is infeasible since it increases the material ductility when a CZ mesh is embedded in the material, as shown in Figure 3. When the material is subject to stresses to deform, the original element size (*d*) will increase to (*d*' + δ), which adds additional elongation δ to the material. Because of this limitation, a scaling factor (denoted as *f*) is introduced here to limit the maximum deflection of CZ elements, as shown by *f*δ*<sup>c</sup>* in Figure 2, and thus to control the chip behavior. Chip behavior is a critical indicator as the cutting force can be affected by the number of cracks during cutting (i.e., work done vs. total fracture energy).

**Figure 3.** A schematic drawing to show unrealistic deformation due to the deflection of cohesive zone (CZ) elements.

When the deflection is scaled to control chip behavior, the fracture energy and, therefore, the cutting force will be scaled accordingly. Thus, the cutting force must be inversely scaled to represent the actual force. To properly apply this model with the scaling factor, the following assumptions are made. First, beyond the elastic deformation, no plastic deformation occurs in the material and all the force contributes to material removal. Second, the specific cutting energy (energy required to remove a unit volume of material) is based solely on the fracture energy. Given constant cutting velocity *vf* and material removal rate (MRR), the cutting force (*Fc*) will be linearly proportional to the specific cutting energy (*p*), as described in Equation (5),

$$F\_{\mathfrak{C}}v\_f = MR\mathbb{R} \cdot p.\tag{5}$$

This implies that the cutting force is scaled linearly with the CZ element's fracture energy. This concept will be validated in the experimental study.

## *2.3. Other Material Properties*

The brittle materials used for the experiment are two types of solid bone-mimetic materials made of high-density polyurethane (PU) foam (Sawbones, Vashon, WA, USA). This material provides consistent and uniform material properties; it is isotropic and does not require a large force to cut. It is ideal for the modeling purpose and experimental validation without extraneous variables such as vibration, impact shock, and heat. These two foams are named based on their densities, 30 and 40 pcf (pound per cubic foot), which equates to 480 kg/m<sup>3</sup> and 640 kg/m<sup>3</sup> , respectively. The 30 pcf has the ultimate strength of 12 MPa and the elasticity modulus of 592 MPa; the 40 pcf has the strength of 19 MPa, and the modulus of 1000 MPa, respectively, based on the manufacturer provided data [19]. The fracture toughness, *Kc*, of these foams is obtained from a separate three-point bending experiment following ASTM D5045–93. The averaged *Kc* for the 30 pcf is 0.46 MPa.m1/<sup>2</sup> and that of 40 pcf is 1.13 MPa.m1/2. The 40 pcf is stiffer and also tougher than the 30 pcf. Based on these properties, the original CZ element properties are calculated in Table 2 below. As seen, the allowed cohesive element deformations are both larger than the element itself (0.01 mm).


**Table 2.** The CZ element properties for the testing materials 30 pound per cubic foot (pcf) and 40 pcf.

#### *2.4. Scaling Factor*

The scaling factor is necessary to control the maximum deflection of CZ elements and thus the material ductility. In the case of 30 pcf, the original CZ deflection goes up to 0.064 mm. With the adjacent element size being 0.01 mm, this allowable deflection is equivalent to a 600% additional elongation (0.064/0.01), which is unrealistic. Figure 4 shows four different scenarios when using the original *Gf* and scaled *Gf* that limits the deflection to be 0.00512 mm (51.2% elongation), 0.00128 mm (12.8% elongation), and 0.00032 mm (3.2% elongation), respectively. As seen in Figure 4a with the original *Gf*, the workpiece and elements experience excessive deformation. Many stretched CZ elements remain alive though the chip has been distorted significantly. Figure 4b shows small but consistent chips generated from the shear plane, which is similar to cutting of brittle metals like high carbon steels. Figure 4c begins to generate fragmented, irregular debris accompanied by dusty pieces, which can be similar to ceramic materials. Figure 4d shows a more extreme case, where the workpiece shatters upon the tool contact. From these simple examples, it can be seen that a fairly small scaling factor is needed in order to force the material to behave as brittle. Note that in the simulation no self-contact is employed because the elements are supposed to support each other via CZ elements. Self-contact is possible but will exponentially increase the computational load due to the larger number of surfaces involved in the contact algorithm.

**Figure 4.** Material responses to the cutting tool with different scaling factors: (**a**) *f* = 1, (**b**) *f* = 0.08, (**c**) *f* = 0.02, and (**d**) *f* = 0.005. The stress is based on the testing material 30 pcf.
