Computer simulations of food oral processing to engineer teeth cleaning

Oral biofilm accumulation in pets is a growing concern. It is desirable to address this problem via non-invasive teeth cleaning techniques, such as through friction between teeth and food during chewing. Therefore, pet food design tools are needed towards optimising cleaning efficacy. Developing such tools is challenging, as several parameters affecting teeth cleaning should be considered: the food’s complex mechanical response, the contacting surfaces topology as well as the wide range of masticatory and anatomical characteristics amongst breeds. We show that Finite Element (FE) models can efficiently account for all these parameters, through the simulation of food deformation and fracture during the first bite. This reduces the need for time consuming and costly in-vivo or in-vitro trials. Our in-silico model is validated through in-vitro tests, demonstrating that the initial oral processing stage can be engineered through computers with high fidelity.


Supplementary discussion 2
This discussion concerns the sensitivity of the numerical − results for ̇= 16.6 mm s −1 and = 0.3, upon the characteristics of the FE mesh. Mesh dependency in the current mastication model can be induced by three main mechanisms. One is through variation in the characteristic length of each element, ch (≈ ), based on which the damage computations are mainly performed; the latter principally influences the crack propagation results, such that these may not converge into a true response with mesh refinement [2]. The second mechanism of mesh dependency relates to insufficient discretisation level i.e. when the elements are too large in size (coarse mesh). The third concerns the element type-shape and specifically for the hexahedral elements used here, the element aspect ratio, , with respect to the direction of maximum absolute strain (here it is the direction of teeth indentation).
Regarding element type-shape, the superiority of hexahedral elements versus tetrahedral elements for this application, is already established [1]; also note that multiple integration points within elements were not employed as this would have led to a prohibitive computational cost. Furthermore, in order to reduce element distortion at large indentations, an optimised element aspect ratio value of three has also been previously proposed [4]; a fixed = / ≈ 3 is therefore used throughout the following mesh sensitivity analysis.
The effect of ch and element size, is investigated (for = 0.3) in two steps. Firstly, a fixed element thickness, = 0.05 mm, is uniformly used along the length of fine mesh region in the food item, while the average element size is scaled by varying the element dimensions , to satisfy the following / ratios: 0.15/0.05, 0.3/0.1, 0.45/0.15, 0.6/0.2. Thus, in this step, the effect of element size i.e. discretisation level is isolated, since food separation is not affected significantly due to a constant ch ≈ = 0.05 mm. The corresponding results are illustrated in supplementary Figure 2(a), which shows that within the studied / range, the − response is insensitive to mesh density up to the onset of breakdown (at = o ). For > o , however, a / = 0.6/0.2 leads to an under-predicted breakdown resistance i.e. force drops steeply, while the other three / ratios continue to give similar predictions, adding credibility to the mesh design employed in the chewing model (corresponding to / = 0.3/0.1). The fact that a / = 0.3/0.1 reduced the analysis time by approximately 50% compared to / = 0.15/0.05, explains the reason for using / = 0.3/0.1 throughout the study. Now for fixed / = 0.3/0.1 (and = 0.3), the individual effect of ch on the results is investigated by varying between: 0.02 mm, 0.05 mm, 0.1 mm, 0.2 mm. In the regime > o supplementary Figure 2(b) shows that the food breakdown resistance increases (smoother drop in force) with decreasing ch . Similar phenomena have been reported in a recent study [2], and have been attributed to the fact that very small elements must deform exceedingly in order to reach the required displacement, f , postulated by the f criterion (see Equation (7)); this affects the response of elements adjacent to the ones undergoing damage, leading to crack tip mesh distortion effects of which the severity varies depending on the mesh density. These effects however can be reduced by regulating the crack tip deformation field between different mesh densities, through the numerical correction strategy described extensively in [2]. Therefore, in order to reduce the mesh dependency shown here in supplementary Figure 2(b), the same strategy is applied. It involves using a strain energy density dissipation parameter, c (kJ m −3 ), in Equation (11), such that f is now defined by: where c is defined as [2]: where ch * is the specific characteristic element length which leads to the best agreement with experimental crack propagation data when f is defined by Equation (11), i.e. when the above strategy is not used. The crack propagation data are obtained from an independent tensile fracture experiment; for this material, this was performed in [2] and gave ch * = 0.05 mm (and thus c = 9.3 kJ m −3 since = 0.93 kJ m −2 ). Therefore, note that the FE results presented in Figures 4−6 remain practically the same when the c parameter is used, since these correspond to ch = ch * = 0.05 mm; for this reason the above strategy was not necessary for studying food breakdown and cleaning efficacy. Nevertheless, when ch ≠ 0.05 mm, supplementary Figure 2(c) shows that using the c parameter has merits in providing more accurate results compared to supplementary Figure 2 Figure 2(c)), thus adding credibility to the model. The same ch range of 0.05−0.1 mm was also found to give accurate fracture predictions for this material both in the single-edge-notched-tensile geometry studied in [2], as well as in the cylindrical indentation geometry reported in [4]. This agreement between different fracture geometries adds further credibility to our numerical modelling methodology.
Regarding potential mesh orientation effects on the results, these are assumed here to be negligible. Specifically, although in Figures 4(c)&(e) the major crack faces in the food item have been shown to approximately coincide with the original plane of element facets, this is triggered by the fact that the original plane of element facets is normal to the direction of the applied maximum tensile strain, 1 (which causes separation into two pieces). It has been generally observed that the morphology of the main crack faces may be influenced by several factors, one of these being the ratio between c and the overall stress levels applied in the food (which depends on ̇) . As a result, although the major crack faces aligned significantly with the original plane of element facets for the conditions ch = 0.05 mm, / = 0.3/0.1, ̇= 16.6 mm/s, this is not true when other conditions are simulated e.g. when ch = 0.05 mm, / = 0.3/0.1, ̇= 166.6 mm/s. Stronger evidence of arbitrary crack predictions by using the current damage scheme can be found in [4]; the former study demonstrated cracking inclined at 45° with respect to the original plane of the element facets. The above lead to the conclusion that the orientation of the mesh does not necessarily predefine the crack path, indicating that our model is a powerful tool for predicting complex fracture patterns during food oral processing.
Supplementary Figure 2. Summary of mesh sensitivity analysis results. (a) shows mandible force-displacement predictions between using four different element average dimensions, / = 0.15/0.05, 0.3/0.1, 0.45/0.15, 0.6/0.2, in the food item, while keeping the thickness (corresponding to the direction of the food specimen length) of all the elements in the fine mesh region fixed at 0.05 mm (for the given direction of cracking here the thickness is practically equal to ℎ ); the parameter was used in Equation (11). (b) shows corresponding results for consistent average element dimensions, / = 3, while varying the thickness of the all the elements in the fine mesh region between: 0.02 mm, 0.05 mm, 0.1 mm, 0.2 mm, and using in Equation (11). (c) shows the same results as in (b) when the numerical strategy against mesh dependency described in [2] is applied i.e. when a term ℎ is used in Equation (11) instead of the parameter . For each condition the corresponding mesh size (number of elements) is indicated, while a friction coefficient of 0.3 is used throughout.