## Abstract

The study of oral processing and specifically cutting of the food piece during mastication can lead towards optimization of products for humans or animals. Food materials are complex biocomposites with a highly nonlinear constitutive response. Their fracture properties have not been largely investigated, while the need for models capable of predicting food breakdown increases. In this study, the blade cutting and the essential work of fracture (EWF) methodologies assessed the fracture behaviour of starch-based pet food. Tensile tests revealed rate-dependent stiffness and stress softening effects, attributed to viscoplasticity and micro-cracking, respectively. Cutting data were collected for 5, 10 and 30 mm s^{−1} sample feed rates, whereas the EWF tests were conducted at 1.7, 3.3 and 8.3 mm s^{−1} crosshead speeds corresponding to average crack speeds of 4, 7 and 15 mm s^{−1}, respectively. A reasonable agreement was achieved between cutting and EWF, reporting 1.26, 1.78, 1.76 kJ m^{−2} and 1.52, 1.37, 1.45 kJ m^{−2} values, respectively, for the corresponding crack speeds. These toughness data were used in a novel numerical model simulating the ‘first’ bite mastication process. A viscoplastic material model is adopted for the food piece, combined with a damage law that enabled predicting fracture patterns in the product.

## 1. Introduction

A wide range of bite forces have been reported in mammals [1,2] and specifically in dogs [3]. It is well established that differences in crown morphology are linked to variations in the properties of foods ingested and masticated [4–6]. Skull and tooth shape are principal determinants of how effectively the mastication is performed [5–7], whereas texture plays a key role in the acceptability of the food [6,8,9]. The need for understanding the breakdown of food and quantifying the required chewing effort arises, when designing nutritional treats for various breeds [6]. Because of the multiple parameters involved in canine chewing and the variability in foods, a reproducible sensory test that predicts acceptability and tooth–food interaction is hard to establish. In contrast, a finite-element (FE) mastication model raises greater potential for reconstructing the real process [9]. This is considered in this study.

Such a model requires the stiffness and fracture toughness as input material properties [9–11]. However, the fracture behaviour of food products is difficult to determine owing to their high compliance, nonlinearity as well as their strong time dependence [9,11–13]. Often, their deformation is dissipative, not obeying the requirement of the linear elastic fracture mechanics (LEFM) that the energy for fracture is consumed locally at the crack tip [14]. Alternative test methods have been developed to determine the fracture toughness of soft solids without the requirement for LEFM. Horgan & Smayda [15] demonstrated the use of the trousers test for the tearing of soft biomaterials. Wire cutting was used to determine the fracture toughness of cheese [11], starch gels [13] and gelatin gels [16]. The orthogonal cutting method [9,17–19] has also been employed. The idea of regarding cutting as a fracture process has been proposed [20], but this was interpreted as requiring an observable crack in front of the cutting tool. Recently [21], it has been observed that a cutting tool may touch the crack tip so obviating the need for cracks to precede the tool. When fracture terms are included, positive intercepts of force values at zero machining depths are expected and provide a method of determining toughness, avoiding excessive crack blunting which would otherwise take place in soft solids. Work on machining processes has highlighted the importance of including the fracture energy in the analysis [17,19,22]. Cotterel & Reddel [23] developed the essential work of fracture (EWF) method for obtaining the energy per unit area dissipated locally in forming the fracture surfaces. It was initially applied in thin sheets to give a plane stress toughness, though efforts have been made to derive a plane strain value [24–27]. Recently, the scheme was extended for rate-dependent materials [28–31].

Although the EWF method has been widely used in metal and polymers, the influence of some experimental variables such as speed and specimen dimensions is not fully understood and the EWF applicability in rate-dependent materials may be limited for this reason [29,31,32]. Moreover, no literature is found on the comparison of the EWF with the orthogonal cutting as methods of measuring toughness. The use of toughness as an input material parameter in FE analyses has been discussed [9–13], but this has never been implemented in modelling chewing. In this work, the fracture toughness of starch-based products is evaluated via the EWF and cutting methods. Background theory is provided, followed by details of the food material used, as well as the experimental procedures adopted. Thereafter, the two fracture tests are compared. Finally, the toughness result is used in an FE mastication model where the fracture paths during the food breakdown are discussed.

## 2. Background

### 2.1. Essential work of fracture

The EWF methodology postulates that the total energy input consumed, *W*_{f}, is partitioned into the work involved in creating new crack faces, *W*_{e}, and that associated with plastic/inelastic deformation encountered in the region surrounding the fracture zone, *W*_{p} [27]. The method employs the double edge-notched tensile (DENT) geometry, in which the ligament is the un-notched section that carries the load. The following assumptions are made.

(i) *W*_{e}, is proportional to the ligament sectional area, such that the specific term, *w*_{e}, is introduced to give
2.1

(ii) *W*_{p} is proportional to the volume of the inelastic zone and the shape factor, *β*, is used to yield
2.2
where *B* is the specimen thickness, *l* is the ligament length (described in §3.3), *β* is the inelastic zone shape factor and the parameters *w*_{e} and *w*_{p} denote the energy consumed per unit fractured area (specific EWF) and the energy absorbed per unit inelastic volume (specific plastic work of fracture), respectively. The *β* factor is determined from the relation between the inelastic zone height and *l*, suggesting three different cases for the inelastic zone shape (circular, diamond and elliptical) [33]. Here, a circular zone is assumed such that *β* = *π*/4 [27]. The energy is normalized over the uncracked (ligament) area, to yield
2.3
where *w*_{f} is the specific work of fracture, *P* is the force and *δ* is the displacement applied on the plate.

A linear correlation is therefore formulated between *w*_{f} and *l* such that *w*_{e} is obtained from extrapolation to zero ligament length. The specific work of fracture, *w*_{f}, is computed as the area under the force–displacement graphs. A detailed description of the procedure can be found in [27,33].

### 2.2. Blade cutting

The theory used here originates from the Williams model [34], illustrated in figure 1. The force *F*_{c} is caused by the movement of the tool in the cutting direction that in turn causes a reaction *F*_{t} owing to the conditions at the chip–tool interface. The thickness of the cut strip, *h*_{c}, may be larger than or equal to the cut depth, *h*, depending on the occurrence of shearing, which is a function of forces *S* and *N*, acting parallel and normal to the tool surface, respectively [34]. By using the Tresca criterion and for yield stress, *σ*_{Y}, a shear plane forms at an angle, *ϕ*, along which the stress reads *σ*_{Y}/2.

Patel *et al.* [17] report on several methods by which the fracture toughness in terms of energy release rate, *G*_{c}, can be determined. The one used here requires the measurement of the chip thickness, *h*_{c} (figure 1). Assuming plane strain conditions, *h*_{c} is given by
2.4

Figure 1 shows a crack tip touching condition [35] where the tip of the tool is in contact with the fracture point. This is modelled by assuming a horizontal force acting at the tool tip, equivalent to *G*_{c}, giving rise to a net axial force of (*F*_{c}/*b* − *G*_{c}). Friction is embodied in the force measurements and most machining analyses do use a friction law because they are concerned with predicting *F*_{c} and *F*_{t} from known property values.

Considering equilibrium along the shear plane equation (2.5) is derived
2.5
where *h*/sin *ϕ* is the length of the shear plane. Equation (2.5) is rearranged to give
2.6

Thus, if a set of experiments is performed in which *F*_{c}/*b* and *F*_{t}/*b* are measured for varying *h*, the *G*_{c} and *σ*_{Y} may be found by plotting ((*F*_{c}/*b*) − (*F*_{t}/*b*)tan*ϕ*) versus *h*/2(tan*ϕ* + (1/tan*ϕ*)). The values *σ*_{Y}, *G*_{c} correspond to the slope and the *Y*-intercept of the linear fit, respectively.

## 3. Experimental

### 3.1. Material

Samples were provided in sheet and square section extrudates by Mars Petcare Leicester UK, from which 60 and 10 specimens were shaped for the EWF and cutting tests, respectively. The starch blend included crude fibres (200–500 µm) at less than 5% w/w and proteins at less than 4% w/w, to form a fibre-reinforced biopolymer with an amorphous matrix of gelatinized starch [36]. A degree of alignment of the fibrous phase to the extrusion direction (principal fibre orientation axis) was revealed by *in situ* scanning electron microscope (SEM) tensile tests (10 kV, 56.8 mm, ×60), where different fibre pull-out lengths were observed when the material was stretched to failure in two principal directions (figure 2). However, the results from these tests (not presented) indicated that the effect of the fibre orientation is minimal on the load–displacement trace, so that a homogeneous, isotropic material may be assumed [37]. This was expected due to the low fibre weight fraction and fibre short length [38]. The samples were sealed after production in order to preserve the water content. All the tests were conducted within 3 days, three weeks following production.

### 3.2. Tension

Dumbbell specimens were cut from the extruded sheets using specimen dies (dimensions shown in table 1). The gauge length, *L _{o}*, was defined by marking two dots on the specimens; these were subsequently optically tracked during the test using an HD Pro webcam C920. Point tracking scripts were set up in Matlab 2014a [39] to measure the elongation,

*L*

_{i}. The stress–strain response was obtained assuming incompressibility, such that 3.1 where

*A*,

_{o}*A*are the original and current cross-sectional area, respectively. In order to verify the incompressibility assumption, the Poisson ratio,

_{i}*v*, was determined from an independent test, where the lateral contraction of the sample was recorded using the same set-up (figure 3). For applied force,

*F*, the Hencky strain,

*ε*, and true stress,

*σ*, were calculated as 3.2 and 3.3

An optical microscope (Gigastone S1–100) was used to image the initial cross-sectional area; these optical data were processed in ImageJ v. 1.49 [40]. Five specimens were tested for each of the three strain rates 0.001, 0.01, 0.1 s^{−1}. In addition, two relaxation tests were conducted via loading up to the strains of 0.06 and 0.08 (0.01 s^{−1} rate), which were then held constant, whereas the stress decay was recorded for 30 min. Stress relaxation experiments are widely used in time-dependent materials [9,11–13], and here they served in increasing accuracy in modelling the time dependency.

### 3.3. Essential work of fracture

The DENT specimens were cut from the extruded sheets, using a scalpel. Sample details are shown in table 1. Two dots were marked across the mid-point of the samples, separated by a distance equal to the ligament length of each specimen (figure 4*a*) [27,28,30]. Varying size notches were created by pushing razor blades gripped in a manual lever press, using a reference line drawn across the specimen width, while the specimen was clamped in a drill vice to ensure notch alignment [24,31,41]. Finally, the razor was used to slide across the crack tip to generate a notch tip of sufficient acuity [31,41]. Note that test data from samples where the cracks were judged not to be aligned were not used in the analysis [27,31,32]. The ligament length, ligament sectional area as well as the remote area were measured post fracture, by imaging one of the fractured faces and processing in ImageJ (figure 4*d*). Sixty specimens comprising five ligament lengths of 6, 8, 10, 12, 14 mm were tested under three crosshead speeds of 8.3, 3.3, 1.7 mm s^{−1} and using four replicates per condition [32,42]. To look at the mechanisms occurring at the fracture process zone, *in situ* SEM tensile tests on smaller DENT samples of dimensions 15 × 12 × 3 mm were conducted (figure 4*e*).

Specimen size criteria were addressed towards investigating the method validity, i.e. whether the true fracture toughness, *w*_{e}, is obtained, independent of sample geometry/thickness [24–26]. For the latter to be achieved, the method requires the ligament to be in a state of either plane stress or plane strain. For sample thickness, *B*, and ligament length, *l*, it has been proposed that the plane strain–plane stress transition occurs for values of *l* in the range 3*B* < *l* < 5*B* [25,33]. Other authors claim that for *l* < 3*B* some plane stress condition may still prevail, though they suggest that the intercept value at zero ligament length (equation (2.3)) still gives the true plane strain, *w*_{e} [24,26]. Another point is that edge effects become important when the total specimen width, 2*D*, is not large enough compared with *l*, such that yielding does not remain confined within the ligament area [24,27]. A restriction of has been suggested [24,28,29], although recently it was argued that values of can be adopted [25,31]. In this work, the sample width and thickness were constrained to be 27 ± 0.25 and 4.9 ± 0.42 mm, respectively, from the extrusion process. Thus, ligament lengths were chosen to vary between 6 and 14 mm (*l* < 3*B*), giving rise to the plane strain conditions.

### 3.4. Orthogonal cutting

Figure 5*a* illustrates the experimental set-up. The cutting tool is ground from tool steel which has a rake angle *α* = 10° and a clearance angle *β* = 11° (figure 1). The cutting tip radius was measured to be approximately 5 mm before and after testing using an optical microscope to ensure there was no tip blunting. Tool sharpness has been shown to influence results [43], particularly for low-toughness materials, but, for those considered here, the 5 mm radius is found to be sharp enough (i.e. smaller than the crack-opening displacement) to give *G*_{c} independent of tip radius. To avoid edge effects, the specimen width was *b* = 10 mm, which is smaller than the 12.5 mm tool width. The tool was mounted on a triaxial piezo-electric load cell (PCB 260A01), capable of simultaneously measuring *F*_{c} and *F*_{t}. The tests were conducted at a rate fast enough to negate any signal decay. The steady-state values of *F*_{c} and *F*_{t} during cutting were used in the analysis (equation (2.6)) [17]. The method involved initially cutting the raw surface to give a smooth reference plane and then moving the specimen upwards to the desired depth; an electronic depth gauge (Mitutoyo UKAS 3533733) measured the position of the surface before and after cutting. The specimen was supported into a clamp fitted onto a modified sledge microtome that was powered via an actuator. Figure 5*b*,*c* shows the clamp set-up.

A series of cuts were performed at feed rates of 5, 10 and 30 mm s^{−1}; for each rate, 10 cut depths were chosen, with values uniformly lying in the range 0.01 < *h* < 0.55 mm. Chip thickness measurements were performed using a micrometer (RS Pro UKAS 7051216).

## 4. Results and discussion

Table 1 summarizes specimen dimensions, test parameters and the main results from all the conducted tests. The last are discussed individually in the following sections.

### 4.1. Uniaxial response

The monotonic tests showed a nonlinear stress–strain response where no distinct yield point can be identified but only progressively decreasing stiffness with strain (figure 6*a*). The latter is suggested to be a synergistic effect of time dependency and observed micro-cracking (figure 6*b*) which degrades the load-carrying capacity [44]. The stress level and strain at break both increased with strain rate (figure 6*a*). It should be noted that the mechanical response including failure in a macroscopic manner was of interest in this work and, therefore, the continuum mechanics approach was adopted [44]. Thus, typical composite damage mechanisms, i.e. fibre pull-out or matrix–fibre debonding were not studied. A viscoplastic material model was employed to reproduce the constitutive response of the composite instead [12].

The viscoplastic model from the material library of the commercial FE software Abaqus v. 6.13 [45] was used, by implementing an initial elastic regime followed by a strain rate-dependent yield behaviour [12]. The input includes the elastic modulus, *E*, Poisson's ratio, *v*, values as well as tabular equivalent stress with respect to equivalent plastic strain data at different plastic strain rates. The plastic strain at increment, *i*, was computed based on the total and elastic strain accordingly,
4.1

The initial yield stress and modulus were defined as *σ _{Y}* = 0.5 MPa and

*E*= 40 MPa. This implies a very small initial linear region to fit the stress–strain curves shown in figure 6

*a*[12]. Two additional values of stress and plastic strain were given for zero plastic strain rate (static), such that the long-term stresses obtained from relaxation are captured (figure 6

*c*). Specifically, equation (4.1) was used for the total strains 0.06 and 0.08, where the material relaxes at stresses 0.7 and 0.9 MPa, respectively. The model calibration shown in figure 6

*a*,

*c*leads to a reasonable agreement with the experimental data.

### 4.2. Essential work of fracture

Figure 7 reports the averaged load–elongation curves. The results were repeatable with a typical variability of 6% in the area computed under these curves. The maximum force increases with speed owing to the higher strain rates encountered, leading to a rise in total work. Figure 8 shows that the EWF (intercept of lines) stayed essentially constant with crosshead speed with values of 1.52, 1.37 and 1.45 kJ m^{−2} for the corresponding speeds of 1.7, 3.3 and 8.3 mm s^{−1}. Figure 9 depicts the crack path during the *in situ* SEM test, where micro-cracking is apparent around the crack tip region. This can be considered as a toughening mechanism [35].

### 4.3. Essential work of fracture method applicability

The test results were examined against the self-similarity criterion, which requires the load–deflection traces for samples of different ligament lengths to be similar [25,28,33]. Figure 7 qualitatively illustrates that the requirement is met, though no literature has been found on quantifying this similarity [32].

Because the current material had never been tested via the EWF method, the validity of the data was further investigated. Williams & Rink [27] proposed another check on the validity of the data, which tests the main assumption that the ligament yields fully prior to fracture. It is required that the theoretical value of the maximum stress on the ligament compares well with the tensile yield stress, *σ*_{Y}, such that
4.2

However, a less restrictive version of the scheme postulates that the maximum stresses, *σ*_{max}(*l*), measured for different ligament lengths, *l*, should fall within 10% of the average maximum value, such that 0.9 [27]. For maximum force, *F*_{max}, ligament sectional area, *A*_{l}, and number of tests, *n*, the *σ*_{max}(*l*) and are given by
4.3
and
4.4

The implementation is depicted in figure 10, which shows that the data points fall close to the prescribed boundaries [31].

As the EWF method was developed for elastic–plastic materials, another question regarding the applicability of the method arises [31]. This was addressed by conducting the analysis using the grip displacement instead of the deflection in the gauge length shown in figure 4*a* [28,30]. Markedly greater (three times) fracture toughness values were obtained (not presented here). The overestimation was due to the deformation outside the plastic zone that was likely to be inelastic [29]. The strains were relieved slowly after the test, embodying a viscous dissipation component in the fracture term as found in [28–30]. These authors revealed a strong effect of sample geometry on *w*_{e} when such conditions occur, whereas they showed that the use of a gauge length equal to the ligament length eliminates the error. This study was the first implementation of the procedure on these pet food products; in order to check the method validity, the EWF results were also compared with those derived from the cutting test.

### 4.4. Blade cutting

Figure 11 reports the force data versus cut depth for the three feed rates. The forces increased linearly with cut depth. Figure 12 shows time delayed images of the cutting process where continuous chip formation and curling are observed. The chip thickness was measured larger than the depth of cut, indicating that shearing occurred [17]. Figure 13 shows the data fitted to equation (3.2) via linear regression, extrapolated to the *y*-axis intercept to determine *G*_{c}. The toughness and yield stress were 1.26, 1.78, 1.76 kJ m^{−2} and 7.8, 12.6, 13.8 MPa for the respective speeds of 5, 10 and 30 mm s^{−1}. The question remains as to whether these yield stresses can be correlated to those obtained from the tensile tests, which are lower but of the same magnitude. High yield stress values have also been noted in polymer and metal cutting and this was attributed to work hardening owing to high shear strains imposed on the chips [17,18].

### 4.5. Comparison of essential work of fracture/cutting

Although the EWF and blade-cutting methods both determined the fracture toughness in terms of the energy release rate, caution was taken when comparing their results. Specifically, the controlling variable differs between the tests, being the crosshead speed in EWF and the feed rate in cutting. Therefore, the optical crack length data were used to relate the crosshead speed in EWF to crack tearing speed; note that, in cutting, the feed rate is identical to crack growth rate. The propagation data for the two cracks either side of each DENT specimen were averaged to give a unique crack length–time curve [30]. Figure 14 presents the average curves of the four replicate tests for each ligament length and speed. These plots demonstrate that the crack propagates with increasing rate followed by a steep rise of approximately 100 mm s^{−1}, indicating unstable crack growth [14]. The incremental speeds were determined via time differentiation and the mean speed values over the whole time range were computed. The results are similar between various ligament lengths for constant crosshead speed, and thus the following analysis was conducted only for ℓ = 14 mm.

Figure 15 illustrates crack speed against time for the three crosshead speeds, together with the corresponding mean speed values 4.2, 8.4 and 16.2 mm s^{−1}. The EWF theory establishes *w*_{e} as the dissipation per unit fractured area on average over the entire crack propagation in the ligament. Therefore, instead of crosshead speeds, *w*_{e} is plotted here versus the average crack speeds.

The EWF method assumes that the whole ligament section yields prior to crack propagation, such that the extrapolation procedure as described by equation (2.3) is valid [33]. It is not certain for this material that this assumption is indeed true. Therefore, an additional analysis was performed using the same data. The latter is based on LEFM, where the energy release rate, *G*, can be calculated by assuming elastic material response through [14]
4.5
where *B*, 2*D*, *C*, *U* and *a* are specimen thickness, width, compliance, stored energy and crack length, respectively. Owing to the section of the specimen being not perfectly rectangular, the remote area, *A*_{r}, was used instead of 2*DB*. The energy calibration factor, *Φ*, is introduced to account for deformations arising from the finite width, 2*D* [14],
4.6

Substituting equation (4.6) in equation (4.5) leads to the following expression: 4.7

The derivation of the calibration factor is shown in appendix A. Figure 16 shows an example where force, crack length, *a*, and energy release rate, *G*, are plotted together versus grip displacement. The area under the graph (ABC) represents stored energy, *U*, for an arbitrary crack growth state [14]. The method was applied up to the maximum computed *G*, which corresponds to crack length of 3.3 mm for the 1.7 mm s^{−1} speed. A substantial drop in *G* follows and instability was assumed to dominate the behaviour thereafter. Initiation was visually observed to occur before the maximum load is reached, such that, for crack growth may occur before the inelastic process zone is fully yielded.

Figure 17*a* demonstrates that *G* increases with crosshead speed for constant crack length. However, when the data in figure 17*a* are re-plotted versus crack speed (figure 17*c*), a unique behaviour is obtained which can be deemed material dependent. The raw crack growth data were used to evaluate *G* in contrast to the raw crack speed results that were smoothed via the five-point mean smoothing function using Matlab (figure 17*b*).

Finally, an average curve was obtained to represent the behaviour shown in figure 17*c* and the comparison between cutting, EWF and the LEFM approach is summarized in figure 18. A close agreement is achieved among the cutting and EWF methods, adding credibility to the test data. Specifically, the average *G*_{c} from the cutting and the EWF methods are 1.6 and 1.45 kJ mm^{−2}, respectively, i.e. a difference of only 9%. However, it is argued that the two methods could yield slightly different toughness results when a fibre-reinforced composite is tested. This is because cutting may not allow for fibre bridging effects to take place as the fibres are being cut by the moving blade [46]. Because the toughness results from EWF and cutting do not differ notably, it is suggested that the fibre reinforcement in these materials is rather poor, as discussed in §3.2. Regarding the LEFM approach, it is acknowledged that the method has limited applicability in this material. The assumption that the deformations are entirely reversible does not hold [14]. Consequently, *G* is overestimated (figure 18).

## 5. Modelling

### 5.1. Progressive damage framework

The viscoplastic model presented in §4.1 was combined with the damage initiation and evolution capability available in Abaqus v. 6.13. The scheme is illustrated in figure 19 and consists of four distinct parts [45]:

— the undamaged viscoplastic constitutive response (a–b–c)

— damage initiation criterion according to the plastic equivalent strain

— damage evolution law based on the measured fracture toughness

*G*_{c}— element removal from the calculations once the material stiffness is fully degraded.

The *σ*_{0} and *σ _{y}*

_{0}are the yield stress and stress, respectively, at the onset of damage, whereas is the equivalent plastic strain at failure. The model captures the reduction in the load-carrying capacity of the element, such that during the analysis the stress tensor in the material is given by the isotropic scalar damage equation 5.1 where

*D*is the damage variable and are the stresses that would exist in the absence of damage. At the onset of damage,

*D*= 0, whereas

*D*= 1 at complete failure. The formation of a crack could entail anisotropic/orthotropic softening [47], such that the load-carrying capacity normal to the crack will degrade fast, but the capacity normal to the crack may not. The current model however degrades all the stress values the same, because equation (5.1) is applied to the full stress tensor,

When damage occurs, the use of the stress–strain relation induces strong mesh dependency. To alleviate this, Hillerborg *et al.* [48] proposed a stress–displacement response after damage initiation. A characteristic length, *L*, is introduced which in Abaqus v. 6.13 is related to the element size. The parameter *G*_{f} is introduced to compute the displacement at which *D* = 1, according to
5.2
where is the equivalent plastic displacement, such that at initiation and during progression. For the linear damage evolution law the *D* increment is computed via
5.3
where the equivalent plastic displacement at failure is computed as
5.4

Moreover, the model assumes that the equivalent plastic strain at the onset of damage is a function of the equivalent plastic strain rate and stress triaxiality, *η*, which is defined as
5.5
where *p* is the pressure stress and *q* is the von Mises stress. For uniaxial compression *η* =− 1/3; for pure shear *η* = 0, whereas for uniaxial tension *η* = 1/3. This is a phenomenological model for predicting the onset of damage owing to nucleation, growth and coalescence of voids in metals, on the basis that a lower is expected with increasing stress triaxiality. Recent work on polymers confirmed that fracture occurs at decreasing values as the stress shifts from compressive to tensile states [49].

Experimentally, it was observed that, under compression, no macroscopic failure occurs, but the material excessively deforms (results not shown). Therefore, the criterion served in preventing elements which experience low stress triaxiality and thus large compressive strains to be deleted, via inputting the artificially large value of for *η* =− 1/3. In contrast, under pure tension the criterion was based on the failure strain obtained from the tensile tests (figure 6*a*). The use of the tensile plastic failure strain as the onset of damage for all the stress states would result in significant removal of elements experiencing compression as the food piece is squeezed between the teeth. Thus, the simulations would predict an unrealistic behaviour. The for *η* = 1/3 was 0.22, 0.29, 0.34 for the respective strain rates of 0.001, 0.01 and 0.1 s^{−1}. Figure 20 depicts the criterion for the 0.1 s^{−1} strain rate only. Abaqus 6.13 interpolates linearly between the two data points corresponding to *η* = 1/3 and *η* =−1/3 and assumes plateau behaviour outside the given range as shown (figure 20). Finally, the material model was fully defined by choosing a linear damage evolution law with a fracture energy *G*_{f} = 1.26 kJ m^{−2}; this was the minimum value deduced from cutting.

### 5.2. Chewing model results

The above-described model was used to define the food material in the mastication model. The product was modelled with a rectangular geometry of 20 mm length and cross section identical to the cutting test samples (table 1). A digitized boxer three-dimensional skull, generated from computed tomography, was used to attain the teeth architecture. The animal was a domestic boxer, donated to research with full consent of the owner. The scanning instrument was the Creaform EXAscan, with a resolution of 0.05 mm and maximum accuracy of 40 µm. The data were imported as an STL format in the RapidformXOR software, which converted the file into a Solid works 2014 editable CAD model. Finally, the geometry was imported into ABAQUS v. 6.13. Three upper (mandible) and lower (maxilla) carnassial teeth were extracted from the CT data.

The skull was scanned at the state of centric relation [50]. It was assumed that the occlusal vector (the trajectory of the mandible) remains linear during the occlusal stroke [7], such that the centric occlusion is reached by translating the mandible along the *y*-axis (figure 21), until tooth–tooth contact occurs [50]. The product was positioned between the middle teeth of the jaws (figure 21). This required a 14 mm mandible displacement (jaw opening) from the centric occlusion state. The maxilla remained fixed, whereas the mandible was constrained in all rotations and translations except for the *y*-direction. The *y*-displacement of 14 mm was applied at the constant speed of 50 mm s^{−1} followed by a recession of 8 mm at 145 mm s^{−1}. The bite cycle period was 0.34 s, which was based on biometric data for the boxer breed [51].

The mesh consisted of 132 079 linear rigid quadrilateral elements (R3D4) and 539 400 three-dimensional continuum linear reduced integration hexahedral elements (C3D8R), discretizing the teeth and the product, respectively. Note that the assumption of rigid teeth is on the basis of a much higher stiffness than that of the food piece [52]. Unreasonable penetrations were avoided by adopting the general contact algorithm, which automatically identifies surfaces in an assembly. Specifically, all the exterior and interior elements of the food product were assigned an element-based surface which was included in the contact domain. Thus, all the element faces could potentially be involved in contact, because the algorithm activates and deactivates faces as necessary when elements fail. Penalty contact was enforced with a coefficient of friction, *μ* = 0.7, obtained from sliding tests conducted in similar materials [53].

The dynamic explicit algorithm was used, which is a special-purpose FE analyser that employs an explicit integration scheme to solve highly nonlinear systems involving complex contacts. A consequence of this is that a large number of increments are required, leading to an increase in computation time, which in this work was alleviated by employing the ‘mass scaling’ technique [10]. This enabled the analysis to be performed without artificially reducing the step time but by artificially increasing the material density, which allows for larger time increments. The density was measured as 1410 kg m^{−2}, and the smallest element dimension was 0.094 mm; the software increased the density of each element, such that a target time increment of 5^{−6} s is maintained. The kinetic energy in the model remained below 5% of the internal energy, such that the applied mass scaling had a negligible effect on the results [10,45]. The simulation was run using eight processor cores (Intel i7 CPU at 3.2 GHz).

Figure 22 shows the damage variable contours versus time, *t*, and corresponding displacement, Δ*L*, of the applied chewing cycle. It can be observed that element deletion first occurs at the exterior sample face at *t* = 0.2 s owing to the tensile strains imposed. Therefore, a fracture pattern initiates (*t* = 0.2 s), followed by crack propagation (*t* = 0.25 s) and final sample separation at maximum displacement (*t* = 0.3 s). Permanent deformations are apparent after unloading (*t* = 0.34 s). In the upper surface of the sample, the elements experienced large compressive strains owing to indentation but were not deleted owing to the damage initiation criterion used. If this was not the case, significant material would ‘disappear’ and this further supports the inclusion of the stress triaxiality parameter. A highly non-uniform strain rate field was obtained owing to the indentation loading; this rate range fell within the tested strain rate range of the tensile tests at which the material model was calibrated.

## 6. Conclusion

This work addresses the need for modelling the constitutive and fracture behaviour of complex biocomposites such as starch-based pet foods. Such studies are very important for animal health and nutrition as well as for aiding biologists in evolutionary studies of species. However, very little work has been reported on these materials. A novel numerical model is presented which predicts food breakdown as well as associated forces applied on the canine teeth.

The study enabled the fracture toughness of starch-based pet foods to be determined and confirmed the equivalency of the EWF and cutting methods. Both tests deduced that the fracture toughness was not highly dependent on rate, within the range studied. The EWF result is based on deflection being measured from a gauge length, such that viscoplastic effects not essential to the fracture process are nullified. The analysis of the DENT data using LEFM did not lead to realistic fracture toughness values. Cutting appears to be a more straightforward test and is also capable of measuring toughness for higher crack speeds.

The mechanical behaviour of starch can be accurately represented in an FE model via a viscoplastic law accompanied with a suitable damage initiation and progression scheme. Tensile tests are necessary to obtain the viscoplastic behaviour and fracture strain. The damage evolution is governed by the fracture energy which is equivalent to the toughness derived from the fracture tests. The model's ability in predicting the breakdown of the food is demonstrated. The chewing model can be replicated via three-dimensional printing of the jaws and performing tests using a universal testing machine with a prescribed displacement–time history, while comparing the FE predicted reaction force output with the experimental data. The load–deflection output could be correlated to texture and help in optimizing food without the need to physically produce different products. The method described in this paper is generic and can be applied to any food material and species.

## Data accessibility

Data supporting this publication can be obtained by request from softsolids{at}imperial.ac.uk.

## Competing interests

We declare we have no competing interests.

## Funding

This work was funded by Mars Petcare UK Ltd.

## Acknowledgements

The authors acknowledge Ewan Armstrong, who helped with the cutting test set-up, and Professor Gordon Williams, for useful discussions. In addition, the anonymous reviewers are acknowledged for their useful comments on an earlier draft of this article.

## Appendix A. Energy calibration factor

As defined in equation (4.6), the calibration factor, *Φ*, requires a solution for the compliance, *C*, of a body with respect to the crack length, *a* [14]. The derivation examines the case of a DENT plate in tension, of length *L*, width 2*D*, thickness *B*, and crack length *a*, at either side of the sample. The LEFM approach requires a linear load–deflection response, such that the compliance is
A 1

For a linear elastic material of modulus *E*, it is also assumed that *G* can be expressed in terms of both energy and the stress intensity factor, *K* (plane strain), to yield equations (A 2) and (A 3), respectively [14],
A 2
and
A 3

The load, *P*, is distributed in the remote section of the plate to give the remote stress, *σ*, as
A 4
Owing to finite dimensions, a finite width correction factor, *Y*^{2}, is incorporated to deduce *K* with respect to remote stress as
A 5
where *Y*^{2} is computed via equation (A 6), with coefficients *A _{o}* = 1.98,

*A*

_{1}= 0.36,

*A*

_{2}= −2.12,

*A*

_{3}= 3.42, for the DENT geometry [14], as A 6

An expression for d*C*/d*a* is obtained by substituting *G*, *P* and *K*^{2} from equations (A 3)–(A 5), respectively, in equation (A 2) to yield
A 7

Introducing A 8 and on integrating equation (A 7) the compliance is deduced as A 9

Substituting equation (A 9) into (4.6) and differentiating with respect to *x* gives the final result:
A 10

The calibration factor was computed for each crack state, *x*, from initiation (*x* → 0) to rupture (*x* → 1). Because the specimen length was *L* = 40 mm and the specimen width was 2*D* = 27 mm, the term (*L*/*D*) = 2.96.

## Footnotes

One contribution of 14 to a theme issue ‘Cutting science in biology and engineering’.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.