Geometrical Simulation Model for Milling of Carbon Fiber Reinforced Polymer (CFRP) Composites

This paper proposes a geometrical 3D milling simulation algorithm for Carbon Fiber Reinforced Polymer (CFRP) milling. In this simulation model, milling tools are simplified into layers of circles while CFRP laminates are simplified into layers of Dexel lines, which can realize simulations for various complex milling conditions. Significant geometrical parameters, for example, cutting angle and cutting length, can be computed with high efficiency. With some geometry-related physical models, the machining results can be predicted for the entire milling process. The effectiveness of this simulation model has been validated by the milling force prediction and the delamination prediction. The performance of this simulation model benefits industrial CFRP manufacturing and provides a new method for online or long-time-interval simulation of CFRP machining.


Introduction
The increasing demand for high-performance materials has never ceased according to the history of science and technology. Today, the harsh operating environment from the aerospace industry is pushing people in search of materials that have Lightweight, hard strength, low thermal expansion coefficient, and high corrosion resistance. The Carbon Fiber Reinforcement Polymer (CFRP) satisfies all these demanding requirements and is therefore dominated in aerospace fields [1]- [3].
However, as one type of difficult-to-machine material, CFRP machining is a major challenge for the manufacturing industry. Milling process with defects including delamination, pull-out fibers and uncut fibers ( [1], [4], [5]) always happens, which can lead to ruinous cost in massive production. To avoid potential losses and to evaluate machining quality, simulations are suggested to be conducted before the milling process to direct processing plans. Many researchers contribute their works to develop the milling simulation for CFRP laminates, and most of these works focus on the Finite Element (FE) model. Kerrigan et al. [6] employed a wireless temperature sensor to measure the heat flux in milling processes and built a FE model to assist calibrating the sensor. Pecat et al. [7], [8] compared implicit and explicit fiber representations in FE model, and found the correlation between the cutting forces and the surface quality. Gao et al. [9] analyzed the influence of milling parameters on machined surface quality and milling force, and created one FE model whose result shows very good agreement with the experiment. In addition, many other researchers have been conducted on CFRP orthogonal cutting and drilling simulations based on FE model, which helps to understand the cutting mechanism of CFRP milling [10]- [12]. From these works, it can be seen that most of the simulations have done predictions under only a few cutting conditions. However, in industrial applications, milling parameters are always changing with some complex toolpath, which makes most of the current simulations inapplicable. Also, the simulation scale is relatively small due to massive computations required by the FE model, which cannot be applied to a largescale condition in real industry. These simulation models are not applicable to the industrial condition.
A couple of studies focused on finding the cutting mechanism of CFRP machining. From these works, the researchers found that geometrical parameters in CFRP milling (e.g., fiber cutting angle and fiber cutting length) played a key role in determining the machined quality and cutting force. Koplev et al. [13] disclosed different cutting mechanisms under CFRP machining according to different orientation of fibers and analyzed the machined quality and milling force. Karpat et al. [14] found the relationship between cutting coefficients and the fiber cutting angle, which leads to a good prediction for CFRP milling force. Hintze et al. [4] investigated the occurrence of delamination and provided delamination criteria based only on fiber cutting angles. Voss et al. [15] studied the influence of the tool geometry and analyze the effect of the tool wear on milling force. Numerous researches bridged relationships between milling geometrical parameters and machining results [16], [17]. These relationships validated the feasibility of milling simulations based on pure geometry calculations.
Therefore, to realize a milling simulation for the entire milling process, a new CFRP contour-milling simulation model based on geometry calculations is provided in this paper. In this model, geometrical parameters can be generated with the given information about the geometry of milling tools, CFRP laminates, and the toolpath. The change of these geometrical parameters, such as cutting angles and cutting lengths, can be captured in the milling process. The record of changing geometrical parameters can not only help to predict milling results with variable milling parameters under complex tool path in time domain, but also assist further studies of cutting mechanisms. This paper is organized as follows: Section 2 introduces simulation model and elaborates simulation algorithm for calculating geometrical parameters. Section 3 shows the details of the simulation settings and validates the geometrical calculation results. Two case studies -milling force prediction and delamination prediction, are presented in Section 4 verify the feasibility of this model. Conclusion follows.

CFRP milling simulation model
Milling is a complex machining process that involves complex physics. Massive factors, including milling parameters, coolant condition, and material properties can influence the cutting mechanism, which can lead to distinct milling results. However, from the pure geometrical aspect, the cutting process between the milling tool and the workpiece is relatively simple. Cutters of the milling tool run in a cycloidal motion and cut into the material, which leads to the Boolean difference between the milling tool and the workpiece. The removed volume, as well as other geometrical parameters like the cutting angle, can be exactly determined by overlapping space between milling tool and workpiece if other cutting defects can be ignored. Also, with these geometrical parameters, milling results can be predicted by using geometry-based models provided by other studies. This paper provides a Dexel based milling simulation algorithm to validate this idea and uses the CFRP material as one example to elaborate this algorithm. In this model, shown in Figure 2, Dexel lines are employed to present CFRP laminates; layers of circles (shown in Figure 3) are adopted as the simulation representation for the milling tool.
For the convenience of elaborating geometrical calculation in this simulation, this paper denotes all the points by P and uses subscript to specify certain point (e.g., for point T), whose coordinate contains position along x axis and y axis according to Figure 1. And use { , ′ } to represent line segments. This simulation employs the Dexel model as the simulation representation for CFRP laminates. Dexel lines, as one of fundamental computer representations, have been widely used in computation graphics due to its low computational complexity and efficient operating performance [18]. Numerous researches have conducted milling simulation with this model for isotropic and continuum materials [19]- [22]. Besides these, Dexel model has its unique advantages for CFRP simulation. With the use of this model, arrays of oriented Dexel lines with the uniform direction and spacing distance can comprise one layer of CFRP laminate, which exactly matches CFRP's property of discontinuity and anisotropy. Also, each fiber in laminates corresponds its own Dexel line entity during cutting simulation, which can assist to analyze and capture the change of geometrical parameters for every single fiber. Further, layers of Dexel lines with variable orientation can represent both unidirectional and multi-directional CFRP laminates, and are able to form various 2D shapes required in industrial applications. All those properties make Dexel model a primarily considered representation for CFRP laminates.

Figure 2 Different CFRP workpieces
This simulation represents CFRP workpiece by layers of fibers. In each layer, fibers align homogeneously with uniform orientation and spacing distance , shown in Figure 2. Different orientation of fibers and different fiber density can be changed by adjusting and , and the shape information of CFRP laminates should be provided to fully determine the entire geometry of workpiece. CFRP laminates with various geometrical properties are able to be built in this model, shown in Figure 2. All geometrical parameters to fully define one layer of CFRP laminates are (1) fiber orientation, ; (2) fiber spacing distance, ; (3) height of fiber layer, ; (4) Shape information (e.g., the rectangular laminate should provide origin position, width and length, shown in Figure 2). Parameters for each fiber should also be given in this model for future calculations, which are fiber line equation (2D line equation in vector form) and existent fiber segments. The existent fiber segment works when fibers are cut into several pieces during the simulation. A list of line segments should be given to record the change of geometrical parameters for each fiber. Existent fiber segments can be expressed by position of its end points, and ′ , denoted by { , ′ }.

Simulation model for milling tool
Choosing the simulation model for the milling tool also needs the consider of its trajectory. For contour milling according to Figure 1, one cutter's trajectory for a linear milling path can be represented as: Where , represent the position of the tooth at time t in cartesian coordinate (t is the variable); v x , v y are instant federate along and direction at time ; is the radius of the milling tool; ω is the rotation speed of the milling tool.
As CFRP fibers have been simplified into Dexel lines, once the intersection points between the cutter's trajectory and fiber lines are calculated, some geometrical parameters, for example, the cutting length per tooth, the immersion angle at one instance, can be conducted. However, the intersection calculation between the cutter trajectory and the fiber lines is too complex, which can generate a transcendental equation without any analytic solutions. Therefore, to increase the computing efficiency, and to derive precise scalar solution for intersection points, proper approximations for the milling tooth trajectory should be made.
In order to make this approximation, milling tool in one-turn rotation is considered. When the tool rotates from 0° to 360°, the time runs Δt = 2 according to the initial time 0 . The equation for the new position (x ′ , y′) of the milling tooth after time Δt can be expressed by: = v y 0 + 2πv y + sin(ω(t + Δt)) x′ = v x ( 0 + Δt) + Rcos(ω(t + Δt)) = v x 0 + 2πv x + cos(ω(t + Δt)) (2) Generally, the feedrate along x and y direction is much smaller than the rotation speed . Therefore, the milling tooth's trajectory in time interval t ∈ [0, 2 ] (in one rotation) can be expressed by: Which is a circular path. This result also works when the milling tool runs in a free-curve path. One secondary infinite-small term is involved, but can be ignored and still lead to a circular-path result. From above, the trajectory for the milling tooth can be approximately expressed by one circular path which moves discretely along the feed direction.
According to this approximation, the milling tool can be simplified into layers of circles which can move along the feed direction. To match the layers of fibers, each circle layer has the same height as its corresponding layer of the CFRP laminate comprised by Dexel lines, which implies the cross-section of the milling circle and that of the CFRP workpiece are in the same height along Z direction. Each layer of the circle can be set to have different geometries, including different radius, different flute position, and different flute numbers. All geometrical parameters to fully define one layer of tool are (1) circle radius, ; (2) tool/circle center position, at time ; (3) height of the circle layer, (the same as its corresponding CFRP laminate layer); (4) helix angle at the current layer, ; (5) rotation direction, (equals to -1 when tool rotates clock-wisely; 1 when counter-clock-wisely); (6) tip position for th cutter , ( , ) at time ; (7) rake angle for th cutter, ; (8) clearance angle for th cutter, , where = 1, 2, 3, ⋯ , .
These parameters allow to simulate CFRP milling with various types of milling tools shown in Figure 3. To reduce the computation load, only the position of cutters is chosen to be displayed by showing lines in the simulation. Angle properties for each cutter (e.g., rake angle) are omitted in the simulation display, but can be considered during the calculation for geometrical parameters.

Figure 3 Tool parameters
With discretizing contour milling by layers, cutting simulation between the most general commercially available milling tool and most common CFRP laminates can be built and defined in this model.

Simulation algorithm
With comprehensive definition of CFRP laminates and milling tools in geometry, proper algorithm can be conducted to get geometrical parameters. In this model, processing time is discretized into steps according to the approximation shown in Section 2.2. In one step, the milling tool rotates from 0° to 360° in reality, while the tool circle moves along the feed direction with the distance of feed per rotation in the simulation. The cutting calculation for one layer of fibers at two adjacent time, t − Δt and t, is shown in in Figure 4, where in time interval Δt, the milling tool rotates only one round as mentioned.

Figure 4 Simulation tasks
In Figure 4, the milling tool initially locates at point −Δt and engages with the CFRP laminate at time t. When time Δt passes, the milling tool moves along the feedrate direction and translates to the new center point . The tool circle's equation can be built up by its center position (x , y ) and its tool radius . Since the tool moves one round between two steps, all the cutters do the same translation without any rotation and keep the same relative position with accord to current circle center +Δt . After the tool's translation, the geometrical calculation can be operated. Four types of basic geometrical parameters (shown in Figure 5) will be derived: (2) Lines Two types of lines will be derived in this part: one is cutting chips and remaining fibers, the other type is reference lines, for example, tangent lines on cutting points.
Cutting chips are represented by overlapped line segments inside the tool circle shown in Figure 6. Also, with the trim of chips, the remaining line segments represent the remaining part of fibers after milling in simulation.
Here, the trimming function, which is the 2D Boolean difference operation, is largely simplified with the use of the Dexel model, for the simulation entities are only circle and lines. Given the segment of the fiber { , ′ } before trimming, and intersection points and ′ (where ≤ ′ , ≤ ′ ), remaining fibers and chips can be calculated as shown in Table 1. Meanwhile, the cutting length for each chip can be derived according to their coordinates, and so does remaining fibers.

No Conditions
Remaining fibers Chips 1 , Table 1 Remaining fibers and chips Reference lines can also be conducted. Reference lines are shown to clarify the geometrical relationship between the tool and fibers, which can contribute to the calculation of angles for the next step. One representative type of reference lines is the tangent line on cutting points, which is the tangent line of the tool circle at the cutting point. Given as the length of reference line, the end of tangent line segment, , can be defined as: Where is tool's center at current step; is the calculated cutting point, and also the beginning point of the tangent line; shows the rotation direction of the milling tool, assigned -1 when tool rotates clockwisely; 1 when counter-clockwisely.
Besides the tangent line, other reference lines can be added according to processing requirements and further calculations.

Figure 7 Angle results
Cutting angle, which is the angle between the tool's edge and the tangent cutting direction, is one of the most dominant factors determining the cutting mechanism during the milling process. In this simulation, the cutting angle is defined by the angle shown in Figure 7. As line segments besides the cutting angle (tangent line and remaining fiber's line) have been derived, various methods can be operated to find cutting angles on each cutting point. Equation (5 is applied to calculate the cutting angle in this model, and cross product is used to involve the influence of the rotation direction. Denote as vectors originating from the cutting angle's vertex in Figure 7, the angle ′ between vectors can be expressed by: With consider of the rake angle, the cutting angle for th cutter should be: Rotation angle is the angle indicating time order, which is defined by the angle between the cutter and the cutting point. For th cutter, vectors defined by ) can be applied in Equation (5 to calculate the rotation angle. Immersion angle is the angle related to the chip thickness which influences the milling force. The immersion angle on certain cutting point depends on the feedrate direction, which can be calculated by the vector { , } and { , } by using Equation- (7, where { , } is the vector which is perpendicular to the feedrate direction, and the coordinates of can be expressed by: Other angles, for instance, engagement angle (originates from tool center between engage-in point and engage-out point), can be calculated with the same method according to Equation (5. Table 2 shows calculation parameters for each type of angle with details.
Angle Table 2 Angle results

(4) Arrays
With respect to the ascending order of the rotation angle, all geometrical results derived in one step can be aligned by sorting algorithms, and their change can be deduced according to time domain. The final geometrical parameters output can be the array whose format is showing below:

Figure 8 Output for geometrical parameters
From this array, some key geometrical parameters, such as rotation angles, immersion angles and cutting angles at the engage-in and engage-out points, can be derived according to their order.
After the derivation of all basic geometrical parameters for points, lines, angles and arrays, chips can be trimmed from fibers. And remaining fibers at current step are the input fibers for the next step. For each step, all layers experience the same simulation process as described above, and keep this operation until the end of the entire simulation. The change of geometrical parameters for the entire milling process can be captured by this simulation algorithm for various milling conditions by defining different tools and CFRP laminates geometries. Reliability of this simulation model is verified in the next section.

Simulation results and model realization
One simulation software has been developed according to the given simulation. The panel for the simulation software is shown in Figure 9 Simulation software interface. Calculated geometrical parameters can be presented by both online chart and offline XML-style file. Online charts include curves of cutting lengths, cutting angles, engagement angles and so on with respect to time, which assists users supervise the trend of the change for output geometrical results during simulation. The offline file offers all calculated geometrical information aligned in arrays mentioned in Section 2.3 for the entire milling process, which serves for users' further studies.

Fundamental simulation experiments
In order to validate the simulation capability of this model, four fundamental simulation experiments have been conducted, as shown in Table 3. In this 2D simulation, a one-flute milling tool with 10mm diameter, 0° helix angle, and 0° rake angle, conducts full-immersion milling on CFRP laminates with different fiber orientations ( = 0°, 45°, 90°, 135°).In this fundamental experiment, we use the cutting angle and the cutting length outputs as the example to show the reliability of this simulation model.
The spacing distance between fibers has set to be 10μm, and the feedrate of the milling tool is 30μm per revolution. Two key geometrical parametersthe cutting length and the cutting angle -are chosen to validate the correctness of this simulation model. The change of these two geometrical parameters in one simulation step can be drawn using offline files. Table 3 shows the change of the cutting angle according to the increase of the tool rotation angle in time domain. It can be seen that the increase of the cutting angle is proportional to the increase of rotation angle. When cutting angle reaches 180°, the milling tool circle is tangent to the fiber line at this cutting point, and the cutting angle suddenly drops to 0°.

Figure 10 Illustration of angles
Analytical equation for the change of the cutting angle has been derived to verify the reliable of the output curve provided by this simulation model. Figure 10 Illustration of angles illustrates the geometrical relationship between the cutting angle and the rotation angle. In Figure 10 It can be seen that curves in Table 3 follows the Equation (8 and hence validate the reliability of the calculation results from the simulation. Table 3 shows the change of the cutting length according to the increase of the tool rotation angle in time domain. It should be noticed that two types of fibers can be cut in one tool rotation, shown in Figure 11. The first type has no free end and its chip forms when both ends of the fiber is cut by the milling tool; The second type has one free end due to the previous cutting and its chip forms when the other end is cut by the milling tool in the current rotation step. In this figure, the cutting length for the first type of fibers is calculated by the half of the total chip length. From Figure 11, it can be seen that, except when the orientation of the fiber is 90°, the maximum fiber cutting length always happens around the tangent point between the fiber line and the tool circle. There is a suddenly drop at the tangent point, since no fiber can be cut when the tool is tangent to the fiber line. For 90° fiber orientation, the fiber cutting length is constant according to the increase of the rotation angle, which matches the real cutting condition since the milling tool is cutting along the fiber direction.
In order to verify the correctness of cutting-length curves, the ideal curve for the change of the cutting length has been derived and its result is shown in Table 3 by dash lines. It can be seen that the cutting length curve fits the ideal curve, except when the tool circle is tangent to the fiber line. In simulation, milling tool cannot always reach tangent points since fibers have spacing distance, while the ideal curve counts the tangent point and therefore gives zero output. Overall, for the fundamental simulation experiment, it shows that the simulation algorithm provides the correct geometrical calculation results. Generally, the diameter of a single carbon fiber with ambient matrix can be approximately 10μm scale. The spacing distance can be set as 10μm to provide a precise simulation.

Condition
Cutting Length Cutting Angle Table 3 Full immersion milling with different fiber orientation (φ = 0°, 45°, 90°, 135°) For micro-milling, or common milling when the workpiece is not too large, setting =10μm shows sufficient computation efficiency. In this condition, this simulation model can also show the micro-scale effect mentioned in Table 3, where the cutting angle and the cutting length cannot reach zero at the tangent point. This micro-scale effect can be enlarged when the micromilling is applied. For instance, a one-flute milling tool with 1mm diameter, 0° helix angle, and 0° rake angle, has been used to conduct full-immersion milling on CFRP laminates when the orientation of the fiber is 45°. When the feedrate of the milling tool is 3010μm per revolution, the change of the cutting angle and the cutting length is shown in Figure 12.
Comparing with Table 3, it can be seen that, though the only parameter changed in this simulation is the radius of the milling tool, the trend of this curve is not the same as Table 3. This is because the fibers' distribution in micro-scale is too sparse, which decreases the number of cutting points and provides less data points in Figure 12. This leads to the difference between two figures. This micro-scale effect can happen in the real milling process according to the discontinuous property of CFRP laminates. Figure 12 The change of the cutting length and the cutting angle when =10μm However, for macro-milling, or common milling with large workpiece, the setting for the spacing distance should be changed. Small spacing distance leads to a large number of fiber lines which can decrease the calculation efficiency or even go beyond the computational power. Therefore, large spacing distance should be applied in this condition to accelerate the calculation in one simulation step. But large spacing distance cannot provide simulation results with high resolution. For instance, four fundamental simulation experiments have been conducted according to Section 3.1. All milling parameters keep the same in these four simulations, except the spacing distance between fibers is 100μm. The change of geometrical parameters in one step is shown in Figure 13. Figure 13 The change of the cutting length and the cutting angle when =100μm Comparing with Table 3, Figure 13 shows large distortion for the curve of the cutting angle and cutting length. This distortion happens mainly near the tangent point between the tool circle and fiber lines. The cutting angle cannot reach 0° and 180° according to this curve. Also, cutting length curves keep constant near the tangent point. This phenomenon is caused by the scarce cutting points near the tangent point.
Therefore, for CFRP macro-milling simulation, proper spacing distance should be calculated before the simulation to guarantee that simulation results have enough resolution, especially near the area when the fiber is tangent to the tool circle. The required spacing distance (whose geometrical relationship is shown in Figure 14) can be calculated by the corresponding central angle of the secant. When the required angle resolution is rq , milling along the direction that is perpendicular to fiber's orientation requires the minimum spacing distance to provide sufficient cutting points. Therefore, the maximum allowable spacing distance max can be expressed by: Which implies the required spacing distance should be smaller than max . Figure 14 Require resolution for the simulation

Tool geometry effect
A simple milling tool with one tooth, zero helix angle and rake angle has been employed in Section 3.1 to conduct four fundamental simulation experiments. However, in real industrial applications, milling tools with complex geometry (e.g., multiple teeth, complex profile, non-zero or changeable rake angle and helix angle) are widely used. Therefore, the simulation model should verify its capability of computing geometrical parameters influenced by these tool geometry parameters. This section shows the influence of these tool geometry parameters on output geometrical results, and verifies the reliability of this simulation model under milling condition with complex tools.

(1) Rake angle
The change of the rake angle does not change the point or the line type of calculated geometrical results. It only changes the angles from the output. For example, two one-flute common milling tools with = 0° and ′ = 15°r ake angle respectively, conduct slot milling while their cutting angle outputs are shown in Figure 15. The feedrate is 30μm per revolution, and the helix angles for two tools are both 0°. Figure 15 The change of cutting angle with different rake angle Compared with Table 3, it can be seen that the rake angle changes the phase of the cutting angle curve. Figure 15 shifts 15° to the right in rotation angle. This phase delay can be explained by the geometry relationship between tooth and fibers shown in Figure 10, which verifies the reliability of this simulation model.

(2) Helix angle
Milling simulation for zero-helix-angle milling tool can be realized by onelayer milling simulation model since the milling condition does not change along z axis. However, the helix angle changes cutter's tip position in different layers, and the milling simulation for a helix milling tool cannot be realized by a one-layer simulation model. The change of the tip position does not influence the value of point-type or line-type geometrical output. It changes the time-dependent angle-type output results, for example, immersion angle and rotation angle, so that the order of the output geometrical parameters are changed. The tip position angle is defined to illustrate the effect of the tip-position change.
is the angle between positive direction of the X axis and the vector of the tooth, expressed by: One simulation example has been conducted to illustrate the change of the output array for one layer caused by the change of the tip position. Here, the simulation condition is the same as the condition mentioned in Table 3 except changing the tip-position angle for the flute to = 45° ( = 0° in Table 3).

Figure 16 The influence of the helix angle
Simulation shown in Figure 16 can be treated as one of the upper layers of a milling simulation with a non-zero helix angle milling tool, where the simulation shown in Table 3 can be treated as the bottom layer with the tip position angle = 0°. With the comparison between Table 3 and Figure 16, it is obvious that varying tip-position angle provides the same magnitude delay to the rotation angle. This delay shifts the cutting angle curve and the cutting length curve to the right. Figure 16 only shows the tip position change effect for one layer, and this effect can be applied to multiple layers. For a non-zero helix angle milling tool, the difference of the tip-position angle between two neighboring layers, ∆ , can be approached by: Where is the radius of bottom layer, ∆ and ∆ are the height and radius difference between these two layers. This approximation shows sufficient precision when the simulation resolution along z-axis is high enough. According to this equation, a 15° helix angle milling tool with 0° rake angle and 10mm diameter has been used for slot milling simulation. Resolution along z-axis has been set as 0.1mm and the radial depth of the cut is 10mm. The orientation of the fiber is 135° with spacing distance =10μm, and the feedrate is 30μm/rev. The result is shown in Figure 17.
It can be seen that by changing tip-position angle, milling tools with nonzero helix angle can be simulated by this simulation model and provides reliable results. For the cutting length curve, adding one flute to the tool does not only provide a rotation angle delay, but also change the magnitude of the cutting length since the teeth cut fibers many times in one simulation step. The precise value of the cutting length for each tool cannot be computed directly due to the approximation of the milling process. The micro-scale effect mentioned in Section 3.2 can also happen on each cutting tooth in one simulation step. However, the magnitude of the cutting length can be calculated approximately by dividing flute number comparing with one-flute milling simulation. This approximation works precisely for the macromilling.
A one-layer slot milling simulation has been conducted with the same condition as that in Table 3 when =135°, except the number of the flute is four. The result is shown in Figure 18. It can be seen that simulation results matches the analysis and provides reliable outputs.
Tool geometrical effects have been fully discussed and verified in this chapter. Milling simulation with complex tool geometry can be conducted with the combination of these effects.

Other effects
Besides the tool geometry, there are some other effects that influence the output geometrical parameters, including settings for multi-layer simulation and the toolpath.
Multi-layer milling simulation occurs when milling tool has helix angle, or when CFRP laminates have different properties for each layer. Generally, the multi-layer simulation can be treated as the superposition of multiple one-layer milling simulations with different tool geometry and CFRP lamina geometry. However, it should be noted that the spacing distance for each layer should be adjusted, especially when micro-scale effect happens. For example, the slot milling with a ball-end milling tool has been simulated by this model. The milling tool has 15° helix angle, 0° rake angle and 10mm diameter. Resolution along z-axis has been set as 0.1mm and the thickness of the CFRP laminates is 10mm. The orientation of the fiber is 45°, and the feedrate is 30μm/rev.

Figure 19 Full immersion milling with the ball-end milling tool
Due to the variation of the tool diameter along Z axis, the milling tool circle at the bottom layer can be very small, where the micro-scale effect should be considered. Therefore, the spacing distance for each layer should be different to provide sufficient simulation precision. According to Equation (9, set the required resolution as rq = 1° and keep the minimum spacing distance as max = 10μm. Then, the spacing distance and the diameter of tool circle for each layer have been plotted in Figure 19. This setting not only guarantees computational efficiency, but also provides sufficient computing accuracy for the outputs. The change of the cutting angle and the cutting length has been presented in Figure 19, it can be seen that no distortion appears from output curves, which shows the precision of simulation results.
The change of milling direction can also influence geometrical outputs of this simulation model. However, complex tool path can be simplified into the combination of short linear toolpath, which enables this simulation model provide prediction for milling with complex toolpath. One example has been conducted to show a simulation under complex milling condition.
With the same ball-end milling tool mentioned in Figure 19, a complex toolpath has been applied on 45°-fiber-orientation unidirectional CFRP laminates. The change of the average cutting length and the average cutting angle has been presented in Figure 20. From analytical results and from output geometrical parameters, this simulation model shows its reliability and efficiency for simulating complex milling processes. With these parameters, milling quality, cutting defects and path optimization can be derived for the entire milling process with some geometry-based physical models. Also, this assists researchers to disclose more cutting mechanisms based on geometrical relationships between the cutter and fibers.

4Physical milling simulation results
Based on geometrical parameters generated from the simulation, the CFRP milling results can be predicted based on some geometry-based models. This section uses two geometry-based models as examples to presents the details for the milling force prediction and delamination prediction.

Milling force simulation
Milling force is the main factor that influences the processing results. The change of the milling force can lead to various stress states on CFRP laminates, which can cause different fracture behaviors. The uncontrolled cracks in CFRP machining does not only generate machining defects, like pull-out and delamination on machined surfaces, but also lead to future failures inside the material [5]. Therefore, a practical and productive simulation model should have the ability to predict the milling force for CFRP material. Therefore, in this simulation, chip thickness model is employed for the force prediction.
Chip thickness model has been originated from [23], which is a commonly used model for milling force predictions. In this model, the chip thickness ℎ at certain cutting point can be approximated by the product of the feedrate per tooth and the instantaneous immersion angle at that point: , and φ and φ are the immersion angle at engage-in and engage-out points, shown in Figure 21.
We define one infinite small value, d , as the height for milling tool's cross section on the chosen milling plane. Forces generated by this cross section (tangential force ( ) , radial force ( ) and axial force ( ) , whose direction is shown in Figure 21) can be derived by multiplying cutting coefficients: [ Since cutting coefficients can be deduced by various ways, including orthogonal cutting [24] and milling force experiments, this method has been widely employed for industrial and research use. For most of the continuous isotropic material, cutting coefficients are constant. However, composite materials like CFRP have complex cutting mechanisms, which leads to the change of cutting coefficients under different milling conditions. Therefore, many researchers have been focusing on finding proper ways to fit the variation of cutting coefficients. Sheikh-Ahmad et al. [25] have employed Neural Network method to build up the function for calculating cutting coefficients. Haiyan et al. [26] have used response surface methodology to fit the curve of cutting coefficients for helical milling with the consider of the cutting depth, spindle speed and feedrate.

Figure 21 Milling force calculation
Besides these methods, Karpat et al. [14] has provided one simple model, which is only dependent on geometrical parameters and can be applied to this simulation. A zero-helix-angle milling tool has been employed to conduct full-immersion milling experiments. According to [14], the cutting coefficients is dependent on the cutting angle and follows the sinusoidal change. Cutting coefficients for tangent and radial directions can be expressed as: = 830 + 410sin(2 + 215) = 0 = 3000 + 1810sin(2 + 175) = 0 (14) Only 2D milling toolpath is considered in the following experiment. Since the helix angle for the tool is 0°, and no axial force is generated, we have =0, =0.
According to Equation (13, following angles should be generated to predict milling force based on chip thickness model: cutting angle , rotation angle , absolute rotation angle abs and immersion angle (and that at engage-in and engage-out point: st , ed ). The angles need to be generated homogeneously with respect to time for each layer, which requires uniform discretization of the rotation angle. Assume ∆ as the minimum increment of the rotation angle in one rotation, after time t = •Δ (where n is the integer which satisfies = 0,1,2, … 2 Δ ), the immersion angle can be expressed by: Where st is the rotation angle at the engage-in point. The value of st , and st can be provided by arrays in Section 2.3. The derivation for the cutting angle has been discussed in Equation (5. Then, the cutting angle and the immersion angle can be generated at any time from one tool's rotation. With Equation (13, , and can be derived with respect to time in one rotation.
The absolute rotation angle γ abs expresses the absolute direction of tangent force in global X, Y, Z coordinate, shown in Figure 21. With γ abs , forces for one layer in local coordinate system (tangential force ( ), radial force ( ) and axial force ( )) can be converted into forces along X, Y, Z direction (denoted as ( ), ( ), ( )). We define ⃗ 1 = (1,0) as the vector along positive direction of X axis; vector ⃗ 2 which is ( − , − ) along tangent line direction; rotation direction as the counter-clockwise direction that is = 1. Then, γ abs can be calculated by Equation (5. The coordinate conversion can be expressed by: cos γ abs −sin γ abs 0 sin γ abs cos γ abs 0 0 Then, the sum of forces ( ), ( ), ( ) for each cutter and each layer can be operated, and the total milling force at one step can be deduced.
One simulation software for milling force prediction has been developed based on the given algorithm to predict the change of milling force for the entire milling process. Full-immersion slot milling simulation has been conducted when the fiber orientation is 0/45/90/135 degrees according to [14]. The diameter of the milling tool is 10mm, and the helix angle and the rake angle are both 0°. Spacing distance of fibers are 0.1mm which satisfies the required resolution rq =1°. The thickness of CFRP laminates is 3mm, and the milling feedrate is 0.03mm/rev. The cutting force along local coordinates has been output shown in Figure 22. Compared with the experiment results from [14], it can be seen that the predicted milling force follows the trend of the real milling force as investigated in [14].
Since superposition principle [14] can be applied when the orientation of fibers varies in different layers, the application of this simulation software can be broadened to more complex milling conditions. This prediction algorithm can be applied for computing milling force with different milling tools shape and complex milling toolpath. Two examples for milling force prediction with complex milling toolpath are shown in Figure 23.

Delamination simulation
Delamination happens when layers of CFRP laminates separate from each other [4]. Delamination largely damages craft's surface integrity and lowers machining precision. CFRP laminates with delamination defects usually require additional trimming process to achieve the required surface quality, which leads to extra expense and time. Therefore, it's necessary to predict the occurrence of the delamination. Hintze et al. [4] have investigated occurrence and propagation of delamination on top layers of CFRP laminates. The delamination prediction criteria from [4] are only based on cutting angles, which can be summarized as below: (1) Delamination occurs where fibers are cut within the critical cutting angle range.
(2) When fibers' initial cutting angle is in the critical cutting angle range, delamination propagates regardless how large the cutting angle is in the following cut.
In [4], a very worn tool has been used with 0° helix angle and 0° rake angle, whose critical cutting angle range is [90,180] degrees. In terms of these criteria, four milling conditions can happen. Delamination occurrence with different fibers' orientations is shown in Figure 24. From Figure 24, it's obvious that for most of the cases (e.g., when fiber's orientation is 0°, 45°, 90°), delamination happens when the cutting angle is within range of 90°-180° according to the criterion (1). However, when the orientation of the fiber is 135°, delamination also occurs, though the cutting angle can be less than 90°. This extra area of defects, is caused by delamination propagation, shown in Figure 25. In Figure 25, one fiber with 135° orientation has been cut 4 times. It can be seen that for the first 3 cuts, the cutting angle for one end of the fiber decreases from 180° to 90°, which means delamination can happen according to the criterion (1). However, for the fourth cut, the cutting angle is below 90° and reaches 45° in the end, which is below 90°. But the delamination still happens since the delamination from the previous cut can propagate to the new cut. This phenomenon can be described by the criterion (2). The fiber's initial cutting angle is 180°, which is over 90°, so that the delamination can propagate and leaves the defect.
Since these criteria are founded purely on geometrical conditions, this can be adopted by our simulation model. We only consider the cases when the simulation runs with the critical cutting angle ranging from 90° to 180°. Once the tool has been changed, the critical cutting angle range can also change. But the same simulation process can be followed to simulate occurrence of delamination for other critical cutting angle ranges.

Prediction algorithm for delamination points
For simulation, criterion (1) can be directly used by comparing the cutting angle provided by the geometrical output. Delamination status function 1 ( ) for criterion (1) can be: Equation (17 indicates that, when delamination happens, 1 ( ) =0; no delamination happens when 1 ( )=1. However, two special cases need to be considered.
One special case happens when the tool circle is tangential to certain fiber, shown in Figure 26 (a), two cutting angles exist at the same cutting point: one is 180° opposite to tool's rotation direction, the other is 0° along tool's rotation direction. For the most of the cases, delamination occurs at the tangential point since the cutting angle is within the critical range. According to Table 1, when the tool circle has tangent points with fibers, condition 1 happens and = ′ . Therefore, the additional criterion should be added into 1 ( ), which can be expressed as: The other special case happens at the engaging-in and engaging-out point. As described above, delamination occurs usually at the tangential point, but this cannot be applied at engaging-in point when immersion angle is 0°, and the tool's moving direction is parallel to fiber's direction, shown in Figure 26 (b). Delamination does not occur since the milling tool does not begin to cut the fiber at this point. The same condition can happen at the engaging-out point, when the immersion angle equals to the engagement angle and milling tool's moving direction is vertical to fiber's direction. Denote the rotation angle at engage-in point and engage-out points as st and ed respectively, and one extra criterion is added into 1 ( ): Criterion (2) is based on Criterion (1), which is more complex since the cutting angle ′ from former steps on the same fiber should be considered.
For instance, in Figure 25 Delamination occurrence and propagation when fiber's orientation is =135°, is the cutting point calculated in the fourth cut, while ′ is the cutting point from the third cut. Line { , ′ } forms the chip in the fourth cut step, and the delamination state of can be influenced by that of ′ . When delamination occurs at ′ , is a delaminationpropagation point.
In summary, the delamination status function for criterion (2), 2 ( ), can be written as: Where 1j represents the delamination state 1 of the cutting point in step by using criteria (1), which can be found from chips linking the cutting point. The entire machining process has n steps in total (where n = 2, 3…). When n = 1, fibers are cut for the first time so no delamination can propagate. When 2 =0, delamination occurs because of its propagation; when 2 =1, no propagation happens.

Figure 26 Special cases for delamination occurrence
The final delamination state for the machined end of the fiber can be described as function : When D=0, delamination happens.

Prediction algorithm for hangout fibers
In paper [4], criteria to predict the maximum length of hangout fibers are also provided. Since fibers cannot be cut normally because of delamination, the longest hangout fiber starts where the delamination happens, and ends at the machined edge. For one fiber, delamination happens when its cutting angle initially reach the critical angle range according to criterion (1). Based on Equation (17, the beginning point of the maximum hangout fiber can be judged and recorded when 1 ( ) = 0 for the first time on one fiber, and denoted by which is a fixed point on this fiber. Denote cutting point on this fiber at current step as , then hangout fiber segment for this fiber can be expressed by { , }.
Though the hangout fiber cannot be cut normally by milling, the fiber still cannot keep its maximum length since the cutting can be conducted even if the separation from its original layer leads to poor contact between tool's edge and fibers. However, after milling, the remaining length of hangout fiber still cannot be predicted due to the lack of full understanding of its mechanism. Therefore, to fit simulation results to the real processing results, two factors, 1 and 2 , are involved to control the existence of the hangout fiber and the length of the remaining hangout fiber. Denote the beginning point of the remaining hangout fiber as ′ , then its coordinate is: Where 1 and 2 are both uniformly random number between 0 and 1, denoted by U[0,1]. When 1 is larger than a constant number , 1 equals to 1 so that the hangout fiber exists, otherwise the fiber can be totally cut. 2 controls the length of hangout fiber.

Realization for delamination simulation
One simulation software is provided according to the algorithm for delamination prediction, shown in Figure 27. One data struct for chips has been defined to record the cutting angle history and their corresponding cutting points. The entire simulation flowchart is given in Figure 27. In this simulation, delamination caused by criteria (1) has been marked by red on its cutting points and hangout fibers, while that by criteria (2) are colored as green. Choosing one of the delamination points can pop out the chart showing the change of the cutting angle history, which assists to analyze delamination states in each step.

Figure 27 Simulation algorithm and simulation software
Four fundamental experiments were simulated and shown in Figure 28 when the orientation of fibers is 90/135/180/0 degrees, and =0.5. The tool runs in a linear path in these simulations, and rotates clock-wisely. It can be seen that, though the model for hangout fiber length prediction is relatively simple, the simulation result is very close to the experiment results from [4]. More researches should be conducted to precisely predict the length of hangout fiber in the future. The delamination prediction with a complex toolpath is also shown in Figure 29.
Overall, it can be seen that this simulation model is capable to predict the milling force and the delamination area for the entire milling process under various complex milling conditions. More physical models can be applied to predict machined quality for long processing intervals with high accuracy and efficiency.

Conclusions
Nowadays, in order to satisfy the strict demand for machining accuracy, various simulation methods have been developed to avoid machining defects for CFRP manufacturing. However, due to the limitation of simulation scale and the lack of computational efficiency, no simulation method is capable to simulate the entire milling process under complex machining conditions. This paper presents a geometry-based 3D simulation model for the CFRP routing process. Simplifying milling tools as layers of circles, and using layers of Dexel lines to represent CFRP laminates, this simulation model can simulate massive milling conditions with various tools and CFRP geometries. With the use of simple and efficient Boolean operations between fiber lines and milling circle tools, the simulation model provides four types of output geometrical parameters, including points (e.g., cutting points), lines (e.g., cutting length), angles (e.g. cutting angle), and arrays, which provides comprehensive key geometrical parameters for entire milling process.
Fundamental simulation experiments compared with analytical equations have been conducted to verify the reliability of this simulation model. Effects caused by micro-scale fiber distribution, tool geometry, CFRP geometry, and tool path have been fully discussed to testify this simulation model's capability and feasibility to compute the change of geometrical outputs under various complex milling conditions. Combining with geometry-related physical models, this simulation model can predict milling process and machined defects for long processing intervals with high efficiency. Two case studies have been provided to validate the effectiveness of the simulator. One case study uses the chip thickness model to predict the milling force with a complex toolpath. The prediction results for the slot milling test are similar to the experiment which validates its reliability for computing cutting force results. The other case study employs the cutting angle-based criteria to predict the delamination area. The simulation results for full-immersion milling also match the real experiments. Both case studies can realize the prediction on a large scale with a complex toolpath. More researches need to be conducted to investigate the hangout fiber and delamination in the future. Other CFRP-related machining processes, for example, drilling and turning, can also be realized based on this simulation model with a minor change.