The Response of Cranial Biomechanical Finite Element Models to Variations in Mesh Density.код для вставкиСкачать
THE ANATOMICAL RECORD 294:610–620 (2011) The Response of Cranial Biomechanical Finite Element Models to Variations in Mesh Density JEN A. BRIGHT* AND EMILY J. RAYFIELD Department of Earth Sciences, University of Bristol, Wills Memorial Building, Queens Road, Bristol, United Kingdom ABSTRACT Finite element (FE) models provide discrete solutions to continuous problems. Therefore, to arrive at the correct solution, it is vital to ensure that FE models contain a sufﬁcient number of elements to fully resolve all the detail encountered in a continuum structure. Mesh convergence testing is the process of comparing successively ﬁner meshes to identify the point of diminishing returns; where increasing resolution has marginal effects on results and further detail would become costly and unnecessary. Historically, convergence has not been considered in most CT-based biomechanical reconstructions involving complex geometries like the skull, as generating such models has been prohibitively time-consuming. To assess how mesh convergence inﬂuences results, 18 increasingly reﬁned CT-based models of a domestic pig skull were compared to identify the point of convergence for strain and displacement, using both linear and quadratic tetrahedral elements. Not all regions of the skull converged at the same rate, and unexpectedly, areas of high strain converged faster than low-strain regions. Linear models were slightly stiffer than their quadratic counterparts, but did not converge less rapidly. As expected, insufﬁciently dense models underestimated strain and displacement, and failed to resolve strain ‘‘hot-spots’’ notable in contour plots. In addition to quantitative differences, visual assessments of such plots often inform conclusions drawn in many comparative studies, highlighting that mesh convergence should be performed on all ﬁnite element models before further analysis takes place. Anat Rec, 294:610–620, C 2011 Wiley-Liss, Inc. 2011. V Key words: ﬁnite element analysis; biomechanics; feeding The use of ﬁnite element analysis (FEA) in biomechanics has been steadily increasing over recent decades. It is now a powerful investigative tool in applied medical and academic research where quantiﬁcation of skeletal stresses and strains is desired. As a computational technique, FEA offers a repetitive means of mechanically testing complex 3D geometries of multiple materials under a variety of physiological or destructive loading conditions. Orthopedic models are often developed from CT scans, allowing for generalized reconstructions to predict mechanical behaviour (Roth et al., 2010), or subject-speciﬁc representations of an individual’s bones for precise reconstructions, which can non-invasively predict how a patient may respond to complex surgical C 2011 WILEY-LISS, INC. V convergence; skull; Additional Supporting Information may be found in the online version of this article. Grant sponsor: Natural Environment Research Council (NERC) studentship; Grant number: NE/F007310/1. *Correspondence to: Jen A. Bright, Department of Earth Sciences University of Bristol, Wills Memorial Building, Queens Road, Bristol BS8 1RJ, United Kingdom. Fax: 0117-925-3385. E-mail: firstname.lastname@example.org Received 18 October 2010; Accepted 22 December 2010 DOI 10.1002/ar.21358 Published online 2 March 2011 in Wiley Online Library (wileyonlinelibrary.com). CRANIAL FEA MESH DENSITY treatment (Remmler et al., 1998). In comparative studies, FEA is being used to investigate the wider association between function and form, by testing the mechanical response to certain loads (Dumont et al., 2005; Bourke et al., 2008; Moazen et al., 2008; Moreno et al., 2008), and by using the fact that morphologies can be manipulated within the computer to isolate the effects of speciﬁc features on the strain environment (Rayﬁeld et al., 2007; Strait et al., 2007; Farke, 2008; Tanner et al., 2008). FEA is also gaining popularity with palaeontologists, for whom it offers a rare opportunity to elucidate mechanically feasible behaviours, and ideally, to address broader ecological and evolutionary questions in fossil organisms (Rayﬁeld, 2005, 2007; McHenry et al., 2007; Wroe et al., 2007; Wroe, 2008; Tseng, 2009). The decreasing cost of computer resources and wide availability of software have made ﬁnite element modelling easily accessible. However, as with all modelling, FEA describes a version of reality based on several simplifying assumptions. In skeletal biomechanics, this usually means the exclusion of soft tissues, or using bone material properties from a related species in the absence of data for the animal under investigation. An incomplete understanding of how these assumptions affect the results can make it possible to build inaccurate models, or models laden with inappropriate assumptions. Thus, while FE models can potentially yield vast amounts of important information, it is vital to ensure that they provide a precise and accurate approximation of the scenario they have been designed to address. Practices for the construction of successful biomechanical models have been outlined, highlighting the importance of mesh veriﬁcation and validation (Viceconti et al., 2005; Anderson et al., 2007). Veriﬁcation, often referred to as precision, is concerned with the underlying mathematics of the model, ensuring that the correct equations are being solved, and that the model contains sufﬁcient elements to approximate a continuum solution (known as mesh convergence). Validation, often referred to as accuracy, is concerned with ensuring that the model is providing a realistic simulation by assigning appropriate boundary conditions. The distinction between the two is important: a model may be mathematically precise, but if it has made incorrect assumptions about the boundary conditions, then it will only be providing a precise answer to the wrong question. Unfortunately, although crucial, both steps are frequently overlooked by FE users, particularly in studies of animal morphology. Commercial software veriﬁes the equations it uses, but the responsibility of performing a mesh convergence test lies with the modeler. Convergence tests require the production of meshes with successively smaller element size, reﬁning the model until enough elements are used to demonstrate that increasing the resolution does not alter stress, strain, or displacement values calculated at the same region of the structure. Biomechanical FE meshes are often based off stereolithography (.stl) ﬁles that reconstruct internal and external surface geometries using triangular facets. The production of even one complex 3D mesh by this method has historically been an incredibly time-consuming exercise. Stereolithography surfaces are often composed of facets whose shapes are of poor ﬁnite element quality (e.g., high aspect ratio), or too ﬁne to solve on a normal desktop computer. As either problem requires that the facets be completely remeshed, rather than 611 just simply subdivided, and the boundary conditions reapplied to the new mesh, convergence testing is rarely performed. When this is coupled with the fact that more reﬁned meshes frequently require prohibitively high computer power, convergence tests often only compare two or three successive meshes (Jones and Wilcox, 2007). These may be misleading as they only capture a small portion of the convergence curve, which would ideally show a rapid initial increase in strain or displacement followed by a plateau once convergence was achieved. In reality, convergence curves may not behave smoothly and important details could be lost (Schmidt et al., 2009). The vertebrate skull has a highly curved and complicated geometry, with numerous, interconnected internal cavities and multidirectional loading regimes. It is of particular interest to zoologists and palaeontologists as it houses the sensory organs and feeding apparatus, thus characterising many of the interactions an animal will have with its environment. Although a number of papers have been published on the convergence of biomechanical models, these are often concerned with medical models of the interaction between bones and prosthesis (Schmidt et al., 2009), or with the relatively simple geometries of the long bones (Ramos and Simões, 2006), phalanges (Richmond, 2007), pelvis (Anderson et al., 2005), or vertebrae (Crawford et al., 2003; Jones and Wilcox, 2007). These studies are not directly relevant to those interested in the behaviour of the whole skull, such as is often the case in vertebrate morphology. Thankfully, the falling cost of computer power and the continuing improvement of mesh generating software are now making convergence testing on highly complex geometries like the skull a reasonable step in model generation. The aim of this study was to perform convergence testing on a complex skull geometry, to gauge what resolution of FE-mesh was required to achieve convergence and at what rate different regions of the skull converged, and why. Convergence tests were performed on specimen-speciﬁc ﬁnite element models of a modern domestic pig skull by comparing strains and displacements in increasingly reﬁned meshes. In particular it was predicted that areas of the skull that experience high strain gradients would converge at higher mesh densities than regions of the skull with low strain gradients, as the position of high strain gradients would be better resolved in meshes with greater precision. These locations are not discernable a priori, but were expected to occur in close proximity to loads and constraints. It was also expected that areas of complex geometric detail would require higher resolutions to appropriately resolve strain than simple geometries. MATERIALS AND METHODS Mesh Generation The skull of a modern domestic pig (Large White breed, age 6 months; skull dimensions 247 141 133 mm) was CT-scanned at the Royal Veterinary College on a Picker PQ5000 medical scanner (0.55 mm pixel size, 1 mm slice thickness, 120 kV, 200 mA). The CT slices were imported into Amira 4.1 (Mercury Computer Systems), and the bony and dental structures were segmented out and used to reconstruct a stereolithography 612 BRIGHT AND RAYFIELD Fig. 1. Examples of (A) low (4 mm); (B) medium (2 mm); and (C) high (0.92 mm) density meshes, illustrating the loss of geometric details at lower mesh resolution. D: Locations of strain gauges used as convergence points in the study, and details of boundary conditions. Note that the applied loads are asymmetrical. TMJ ¼ temporomandibular joint; M1 = 1st molar tooth; G ¼ gauge. (.stl) ﬁle. This was then imported into HyperMesh 10.0 (Altair Engineering) for ﬁnite element pre-processing. A series of successively reﬁned 2D shell meshes were generated using the HyperMesh ‘‘shrink-wrap’’ function, which ﬁts a 2D shell as tightly as possible over an existing geometry, with a user-controlled modal average element size. This results in meshes of homogeneously sized elements with good aspect ratios, and also retains internal geometries, such as the brain and sinus cavities. In instances where the desired element size is too large to ﬁt to the local curvature, the shrink-wrap function may adopt a slightly smaller element size to capture ﬁne details (e.g., cusps on the teeth), however, with the very coarsest meshes, a certain degree of ‘‘smoothing out’’ of the ﬁner geometric details is unavoidable (Fig. 1A–C). These 2D shells then formed the starting point for an automated, surface-based tetrahedral meshing algorithm to generate a series of 3D FE meshes. Initially, ﬁve meshes of element sizes 1, 2, 3, 4, and 5 mm were generated, then 13 subsequent ﬁner meshes were created post hoc, 11 with element sizes between 1 and 2 mm. The number of elements in a model increases rapidly with decreasing element size, therefore a large increase in the number of elements was encountered between 1 and 2 mm, necessitating additional models to appropriately capture convergence detail. This resulted in a total of 18 meshes, with modal average element sizes between 5 and 0.83 mm (Supporting Information 1). Boundary Conditions Once the 3D meshes had been constructed, we applied boundary conditions that mimicked the set-up of an in vitro strain experiment that was being used to validate the model as part of a wider study. This set-up utilized loads totaling 755 N applied to the masseter and temporalis muscle attachment sites, which pulled the 613 CRANIAL FEA MESH DENSITY specimen down on to supports located bilaterally at the ﬁrst molars and temporomandibular joints (TMJ) on a custom built-testing rig. This provided a reasonable approximation of a biting scenario that was repeatable and easy to model (Fig. 1D; Supporting Information 2). Loads were applied to the FE-model via rigid body elements (RBE3 in HyperMesh) and constraints were deﬁned to prevent translation in the Y axis at the teeth (dorso-ventral motion) and the XYZ axes at the TMJ, to mimic the experimental loading rig. Loads were initially applied to the original .stl import and then, using this as a guide, were used to place the loads on to all the subsequent meshes to ensure repeatability of load positions. Although every effort was taken to ensure that the boundary conditions applied were as consistent as possible amongst models, the changing resolution of the meshes meant that the exact locations on the original mesh did not always match with a node on the shrinkwrapped meshes. In these instances, the nearest node was used to apply the loads, however this meant that, in the coarsest meshes, the positioning of loads and constraints may have been in error by up to 3.5 mm. In the ﬁner meshes, this error was <1 mm. Given that the shortest dimension (height) of the specimen measured 133 mm, the maximum placement error of <3% was deemed to be acceptable. At this stage, node sets corresponding with the location of sixteen 45 mm2 planar rosette strain gauges (C2A-06-062LR-350; Vishay Measurements Group UK, Basingstoke, UK) applied to the experimental pig were deﬁned, again initially on the original .stl and then on the other meshes (Gauges 1–16, Fig. 1D). Deﬁning the node sets at this stage ensured consistency amongst the models, and allowed for strain data to be queried in these exact locations once the analyses were complete. Inevitably, as the mesh resolution increased, so did the number of data points queried per gauge in each subsequent analysis (Supporting Information 3). In the absence of material properties data on pig bone, all models were assigned the properties of adult human cranial bone (E ¼ 12.5 GPa, m ¼ 0.35, Peterson and Dechow, 2003). Bone was assumed to be isotropic and homogeneous, and the teeth were also assigned the material properties of bone. Cranial sutures were not considered in these models. Although investigation continues into the sensitivity of biomechanical FE models to these assumptions, and their effects on model validity (Ross et al., 2005; Strait et al., 2005; Kupczik et al., 2007; Panagiotopoulou et al., 2010, 2011; Rayﬁeld, 2011), the convergence test seeks only to verify, not validate. All models were exported with both linear 4-noded (TET4), and quadratic 10-noded (TET10) tetrahedral elements, resulting in a total of 36 models: 18 with TET4 and 18 with TET10. The size of the elements in a given model remains the same, but TET10 quadratic elements have additional nodes located at the mid-point of the edges of each element, allowing more complex deformation to be modeled. TET4 elements are traditionally considered to be overly stiff (Mac Donald, 2007), and so observing both types of element would determine any difference in ‘‘convergence rates’’ (i.e., whether convergence would be achieved earlier in lower resolution meshes in one element type compared to the other). Finite Element Analysis The 36 models were solved in Abaqus 6.8.2 (Dassault Systèmes Simulia, Providence RI) on a desktop PC (Windows 64-bit Vista Business, Intel Xeon x5450 3.00 GHz CPU, 64 GB RAM). At this speciﬁcation, the lowest resolution model (5 mm TET4) ran in 7 sec, and the highest resolution model (0.83 mm TET10) ran in 10.6 hr. For each model, maximum (emax) and minimum (emin) principal strain and displacement (U) values were exported from each of the 16 node sets representing the strain gauges, and contour plots of the models were visually inspected. Strain and displacement values from nodes under each gauge site were averaged, and plotted against the number of elements in the mesh. Nodes selected only considered element apex nodes, and did not include midpoint nodes from quadratic models, to ensure consistency when comparing with linear models. The error between successive models was calculated as the percentage difference in strain or displacement magnitudes between the current mesh and the previous (less reﬁned) mesh (Schmidt et al., 2009) at each gauge location. We chose to identify the convergence where <10% or <5% change in strain or displacement magnitudes was found between successive meshes. To our knowledge, there appears to be no standard criterion of convergence, although <5% is often quoted (Anderson et al., 2005; Schmidt et al., 2009). Ultimately, the target level of convergence that is considered appropriate, be it 10% or 0.1% change, is left to the user, who must decide whether these errors are acceptable within the remit of their analysis. The highest density model (0.83 mm, TET10) was the most precise, and was therefore used to determine the strain and displacement ranges (max.–min. value) under each gauge site, thus indicating the gradients encountered. Distorted elements (deﬁned by Abaqus as having Tetrahedral Quality Measure <0.02; Supporting Information 1) were identiﬁed for each model, and made up <0.01% of the total elements per model on average [greatest value at 5 mm mesh (32,000 elements) ¼ 0.025% distorted elements, see Supporting Information 1]. RESULTS TET4 Elements All element numbers quoted in the text are rounded to the nearest 10,000. For actual values please see Supporting Information 1. The results displayed in Fig. 2A– C (see Supporting Information 4 for numerical values) show that, by 1.05 mm (1,250,000 elements), all TET4 gauge sites have converged to within 10% change in values from meshes with successively smaller element size, with two exceptions: emax G7 (12%) and displacement G2 (11%). By 0.92 mm (1,750,000 elements), almost all variables are converged to within 5%, with the following exceptions: emax G4, G7, G9; emin G2, G3, G7; displacement G1, G2. This indicates that different locations in the skull are converging at different rates, with some reaching the target of 10% convergence at far lower mesh densities than others: emin G8, for instance, achieves 5% convergence in a mesh as coarse as 2.1 mm (250,000 elements) (Fig. 2B). 614 BRIGHT AND RAYFIELD Fig. 2. Distribution of convergence in the series of models for (A) emax; (B) emin; and (C) displacement. Dark gray indicates convergence within 5%, light gray within 10%, and white >10%. Anomalous points are marked as ‘‘x.’’ TET10 Elements The TET10 models, as expected, are less stiff than their 4-noded counterparts (Fig. 3). Although both types of element generally follow the same convergence curve, with the TET10 results plotting with slightly higher values of strain (except G5; see Fig. 3) and displacement, the differences between TET4 and TET10 models are not consistent. In extreme examples, the TET10 models may give anomalously high peaks of strain, evident as clear outliers (marked by ‘‘x’’ in Fig. 2A,B; Fig. 3: G2, G12, G13, G14, G16; bold numbers in Supporting Information 4), which are absent in the TET4 results. We deﬁned anomalies as wherever convergence could be demonstrated to within 10% between two or more successive models on either side of the anomalous point, or wherever the change in strain between two models was over 100% more than the change between the models preceding it. The positions of these anomalous peaks are inconsistent: for instance, at 1 mm (1,500,000 elements) peaks are present in G16, but not G15, to which it is directly adjacent (Fig. 1D). Nor are all anomalous peaks in the same model (e.g., G2 shows peaks at 1.05 mm, rather than 1 mm as in the previous example). Removal of the anomalous points and recalculation of the percentage change based on the previous model (along with qualitative observation of the convergence plots) suggests that, without these anomalies, 5–10% strain convergence would be obtained by 1.05 mm or earlier (Fig. 4A,B). In all models without clear anomalies, 10% convergence is achieved by 0.92 mm, with the following exceptions: emax G7; displacement G1, G2, G4, G6, G7. Displacement results therefore appear to be less stable than those of strain using TET10 elements, with no gauges achieving 5% convergence by the highest resolution model (0.83 mm), whereas 14 out of 16 gauges reach this by 1.05 mm using TET4 (Table 1). Relationship with Strain Gradients We predicted that regions with the highest strain gradients would take longer to converge than low strain CRANIAL FEA MESH DENSITY Fig. 3. Convergence plots of strain (solid lines) and displacement (dashed lines) at all gauge sites (G1– G16) for TET4 (closed markers) and TET10 (open markers) models. 615 616 BRIGHT AND RAYFIELD Fig. 3. (Continued) 617 CRANIAL FEA MESH DENSITY Fig. 4. Distribution of convergence in the series of models following the removal of anomalous points (marked by ‘‘x’’), for (A) emax and (B) emin. Dark gray indicates convergence within 5%, light gray within 10%, and white >10%. TABLE 1. Number of gauges out of a total of 16 achieving convergence by 0.92 mm (1,749,149 elements) emin emax TET4 TET10 Displacement 5% 10% 5% 10% 5% 10% 13 8 15 11 13 7 16 13 13 0 15 11 gradient regions. Strain ranges predicted by the highest resolution model (0.83 mm, TET10) were obtained from each gauge site, and are presented in Table 2. G5 and G8 have the highest strain ranges and G2 has the lowest. This was expected, given their proximities to loading on the zygomatic arch (Fig. 1D). Contrary to predictions, however, G5 and G8 converged very rapidly, whereas G2 was one of the slowest gauges to converge (Fig. 2A–C). Association with Geometry The convergence locations we chose were equivalent to the positions of strain gauges in an associated in vitro validation study. The placement of strain gauges requires that the bones be relatively ﬂat and featureless, and so curvature of the bones should not be affecting the convergence results. However, G7, which fails to achieve convergence within 10%, is directly adjacent to a prominent groove in the thin anterior frontal bone, associated with a large blood vessel (Fig. 1D). No other gauges are in such close proximity to notable geometric features. Consideration of Contour Plots Most paleontological use of ﬁnite element analysis is concerned not with the absolute values generated by the models, but a more qualitative comparison of the contour plots they produce. Such comparisons allow conclusions to be drawn about the relative performances of two different structures (Rayﬁeld et al, 2007; Jasinoski et al., 2009). Qualitative observations of the contour plots produced by the convergence test show how these change with increasing model resolution. Low resolution models tend to underestimate values and miss ‘‘hot-spots’’ of higher strain (Fig. 5). By 1 mm element size, minor differences in the location of high strains are still discernable, but contour plots have generally settled on a stable pattern (Fig. 5P,H0 ). Contour plots for displacement are more homogeneous, but differences in the rostrum are still apparent even in the highest resolution models (compare Fig. 5A00 –B00 ). DISCUSSION The results show that in most areas of the skull, TET4 models converge to within 5% and TET10 to within 10% by 0.92 mm (1,750,000 elements). The relatively poor performance of the TET10 models was surprising, as they often failed to converge as rapidly as their TET4 counterparts. When strain is the variable of interest, this apparently poor performance can be attributed to the 618 BRIGHT AND RAYFIELD TABLE 2. Predicted strain ranges at each gauge site in the highest resolution model (0.83 mm TET10) Range Gauge site G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 G11 G12 G13 G14 G15 G16 Displacement (mm) emax emin 3.9E-05 9.2E-06 9.8E-05 7.0E-05 3.3E-03 1.3E-04 5.5E-04 5.0E-03 1.6E-04 4.6E-05 6.4E-05 1.1E-04 9.1E-05 1.3E-04 2.9E-05 6.6E-05 3.0E-06 1.0E-06 1.0E-05 8.1E-06 2.0E-04 3.1E-06 5.7E-05 1.4E-04 6.7E-06 5.3E-06 2.3E-06 4.3E-06 2.7E-06 2.5E-06 1.8E-06 2.8E-06 2.0E-06 1.4E-06 6.1E-06 1.3E-05 1.2E-04 4.6E-06 2.7E-05 7.6E-05 6.0E-06 5.1E-06 3.9E-06 7.5E-06 5.6E-06 5.5E-06 4.2E-06 4.3E-06 Bold normal type numbers indicate the highest values and bold italics indicate the lowest values for each variable. presence of anomalously high strains at certain locations (‘‘x’’ in Fig. 2A,B; anomalous peaks in Fig. 3; bold numbers in Supporting Information 4). The reason for the anomalous peaks is unclear, given that the only difference between TET4 and TET10 models is whether the elements are linear or quadratic, and by the stage of mesh reﬁnement that most anomalies appear (1 mm), the maximum possible discrepancy in boundary condition application point between models was < 1 mm. As the percentage of distorted elements was low (< 0.01%) and bad elements were located away from the gauge sites, element distortion is unlikely to be responsible. If the anomalous data points are disregarded, TET10 models converge to within 5% relative difference at element sizes equal to or larger than the TET4 models (compare Fig. 2A,B with Fig. 4A,B). However, it can be reasonably demonstrated that, by 0.92 mm element size, both the TET4 and TET10 models are converged in regards to principal strain to within 5% at all locations except G7 (emax 12–13%, emin 7–8%). As noted above, the G7 site is close to an area of complex geometry, which may be responsible for this failure to converge satisfactorily. Fig. 5. TET4 contour plots of emax (A–R), emin (S–J0 ), and displacement (K0 –B00 ) in order of increasing resolution. 619 CRANIAL FEA MESH DENSITY When displacement is the variable of interest, most TET10 model locations converge to within 10% by 0.92 mm element size (1,750,000 elements), but none to within 5%. This is an unexpected result, as displacement is often demonstrated to converge at lower mesh densities than strain (Mac Donald, 2007). Conversely, in TET4 models 13 gauges achieve 5% convergence by this stage (Table 1). One would expect that the increased resolution offered by TET10 elements should improve the results of these models, but perhaps the greater degree of deformation permitted by TET10 elements slows their rate of displacement convergence. We note that the convergence curves, and consequently the magnitudes of relative error, are not smooth. In all gauges, even those that ultimately achieve convergence, successive models may report a relative error of <5%, only to jump up to higher errors in the next model in the sequence (Supporting Information 4). This is a similar observation to that of Schmidt et al. (2009), who found that Von Mises stress sometimes increased in their models after a period of apparent convergence. These periods of apparent convergence followed by jumps in strain present a dilemma: how can we be certain that convergence has been achieved in our ﬁnest models when further resolution increases may show another strain jump? In the two ﬁnest meshes modelled (0.92–0.83 mm), a decrease of <0.1 mm element size resulted in a model with nearly 450,000 more elements. Conversely, a size decrease of 0.125 mm for models earlier in the analysis (1.500–1.375 mm) gives element gains of only 130,000. Because of the greater range of mesh densities represented towards the end of the analysis, and that stable plateaus are observed in the length of the convergence curve that we have covered, we are conﬁdent that, in gauges we have deemed converged, further jumps are unlikely. Of course, as ﬁnite element models represent discretizations of continuous problems, the only way to be certain that smaller element sizes would not see jumps in their results would be to increase mesh resolution ad inﬁnitum, although in practice this is clearly not possible. Our results therefore support Schmidt et al.’s (2009) assertion that mesh convergence studies comprising only two or three mesh densities may be insufﬁcient. For instance, the peaks of anomalously high strain identiﬁed in this study would not have been apparent without the wider assessment of error that identiﬁed them as outliers, nor could we conﬁdently predict that further element increases will not result in strain jumps. Predictions that gauge sites subjected to higher strain gradients would be slower to converge were shown to be incorrect. The gauge sites with the highest gradients (G5 and G8) in fact showed the fastest convergence rates. These sites were located on the zygomatic arch, close to the application of the loading. On further consideration, this result is perhaps not that surprising. The sites closest to the loading experience the highest absolute strains in the model. Consequently, to produce a variation of 5% requires more absolute change in strain values at these locations than at sites with lower average strain values. Similarly, sites with lower average strains will be very sensitive to small changes in the model, and thus require higher resolution meshing to achieve convergence. However, depending on the nature of the question being addressed, it is possible that regions of low strain are not of interest to the investigator. If this is the case and absolute values are unimportant, then perhaps lows strains can be ignored. Then, it is sufﬁcient only to resolve the model so that (i) areas of high strain converge and (ii) areas identiﬁed as low strain are indeed areas of low strain (rather than an artefact of the model being under-converged; Fig. 5). The meshes used in this study can be thought of as homogeneous in terms of element size. Although this approach is ideal for a convergence study of the whole skull, the observation that different locations converge at different rates indicates that element size homogeneity is not the most efﬁcient meshing method: some gauges converge at lower mesh densities than others and do not require the computationally expensive high resolution forced upon them by homogeneous high-density meshes. Similarly, even regions that require a high mesh density to converge, such as some of the low strain areas identiﬁed in this study, may be of little interest to the investigation. As these regions would not require the most precise solution, a coarser mesh may sufﬁce in such locations. Gauge 7 was unable to show 10% emax convergence even at the highest mesh density considered here. Adaptive re-meshing, where troublesome areas are re-meshed with ﬁner elements after identifying them by the homogeneous approach described above, is the solution to this efﬁciency dilemma. However, this is currently difﬁcult to put into practice using the automated meshing software utilized by most morphologists attempting to capture complicated and highly curved 3D structures. Surfaces generated from .stl ﬁles automatically apply smaller triangular facets to capture highly curved geometry, but conversely, use larger facets for ﬂat areas. As these local variations in facet density are dictated by geometry and not strain, it cannot be assumed that a direct conversion of these facets into elements will give a model that is appropriately converged: for instance, strains developed in ﬂatter regions may not be fully resolved by large elements the size of the original facet (all locations in this analysis were relatively ﬂat, due to the need to attach strain gauges). In addition, errors may arise from the close proximity of very large to very small elements that may be generated in an .stl. Despite its qualitative nature, the fact that contour plots also show sensitivity to mesh density is important. Low resolution meshes underestimated strain and failed to resolve strain ‘‘hot-spots’’ in the contour plots. This may signiﬁcantly alter an interpretation of structural behaviour, or mislead comparisons with a validation dataset. In palaeontology, comparative studies often focus on the response of morphologies to idealised loads and interpret differences in contour plots because absolute values cannot be validated. In such studies, different morphologies will likely converge at different rates and as such, veriﬁcation of all morphologies being considered is particularly crucial so that inappropriate mesh resolutions do not confound the true effects of differences in morphology. CONCLUSIONS For most locations, convergence in this study can be said to have been successful by element sizes of 1 mm. It is important to emphasize that the speciﬁcs of this convergence study are unique to this particular model 620 BRIGHT AND RAYFIELD set-up. Convergence tests should ideally be performed in all ﬁnite element analyses, and the current study highlights mesh veriﬁcation as an important step before any validation or interpretation of models of the skull can occur. If validation cannot be performed and only contour patterns are of interest, convergence is still vital to ensure that the correct patterns have been resolved. ACKNOWLEDGMENTS The authors thank Daniel Nieto (Altair UK) for his extensive help and advice on mesh generation, Chris Lamb (Royal Veterinary College) for scanning the specimen, and Michael Fagan (University of Hull) for his assistance in resolving early mesh errors, as well as two anonymous reviewers for their constructive comments. LITERATURE CITED Anderson AE, Ellis BJ, Weiss JA. 2007. Veriﬁcation, validation and sensitivity studies in computational biomechanics. Comput Meth Biomech 10:171–184. Anderson AE, Peters CL, Tuttle BD, Weiss JA. 2005. A subject-speciﬁc ﬁnite element model of the pelvis: development, validation and sensitivity studies. J Biomech Eng 127:364–373. Bourke J, Wroe S, Moreno K, McHenry C, Clausen P. 2008. Effects of gape and tooth position on bite force and skull stress in the dingo (Canis lupus dingo) using a 3-Dimensional Finite Element approach. PLoS ONE 3:e2200. Crawford RP, Rosenburg WS, Keaveny TM. 2003. Quantitative computed tomography-based ﬁnite element models of the human lumbar vertebral body: effect of element size on stiffness, damage, and fracture strength. J Biomech Eng 125:434–438. Dumont ER, Piccirillo J, Grosse IR. 2005. Finite-element analysis of biting behaviour and bone stress in the facial skeletons of bats. Anat Rec 283A:319–330. Farke AA. 2008. Frontal sinuses and head-butting in goats: a ﬁnite element analysis. J Exp Biol 211:3085–3094. Jasinoski SC, Rayﬁeld EJ, Chinsamy A. 2009. Comparative feeding biomechanics of Lystrosaurus and the generalised dicynodont Oudenodon. Anat Rec 292:862–874. Jones AC, Wilcox RK. 2007. Assessment of factors inﬂuencing ﬁnite element vertebral model predictions. J Biomech Eng 129:898–903. Mac Donald BJ. 2007. Practical stress analysis with ﬁnite elements. Dublin: Glasnevin Publishing. McHenry CR, Wroe S, Clausen PD, Moreno K, Cunningham E. 2007. Supermodeled sabercat, predatory behaviour in Smilodon fatalis revealed by high-resolution 3D computer simulation. PNAS 104:16010–16015. Moazen M, Curtis N, Evans SE, O’Higgins P, Fagan MJ. 2008. Combined ﬁnite element and multibody dynamics analysis of biting in a Uromastyx hardwickii lizard skull. J Anat 213:499–508. Moreno K, Wroe S, Clausen P, McHenry C, D’Amore DC, Rayﬁeld EJ, Cunningham E. 2008. Cranial performance in the Komodo dragon (Varanus komodoensis) as revealed by high-resolution 3-D ﬁnite element analysis. J Anat 212:736–746. Panagiotopoulou O, Curtis N, O’Higgins P, Cobb SN. 2010. Modelling subcortical bone in ﬁnite element analyses: a validation and sensitivity study in the macaque mandible. J Biomech 43:1603–1611. Panagiotopoulou O, Kupczik K, Cobb SN. 2011. The mechanical function of the periodontal ligament in the macaque mandible: a validation and sensitivity study using ﬁnite element analysis. J Anat 218:75–86. Peterson J, Dechow PC. 2003. Material properties of the human cranial vault and zygoma. Anat Rec 274A:785–797. Ramos A, Simões JA. 2006. Tetrahedral versus hexahedral ﬁnite elements in numerical modelling of the proximal femur. Med Eng Phys 28:916–924. Rayﬁeld EJ. 2005. Using ﬁnite-element analysis to investigate suture morphology: a case study using large carnivorous dinosaurs. Anat Rec 283A:349–365. Rayﬁeld EJ. 2007. Finite element analysis and understanding the biomechanics and evolution of living and fossil organisms. Annu Rev Earth Pl Sc 35:541–576. Rayﬁeld EJ. 2011. Strain in the ostrich mandible during simulated pecking and validation of specimen-speciﬁc ﬁnite element models. J Anat. 218:47–58. Rayﬁeld EJ, Milner AC, Xuan VB, Young PG. 2007. Functional morphology of spinosaur ‘‘crocodile-mimic’’ dinosaurs. J Vertebr Paleontol 27:892–901. Remmler D, Olson L, Ekstrom R, Duke D, Matamoros A, Matthews D, Ullrich CG. 1998. Pre-surgical CT/FEA for craniofacial distraction: I. Methodology, development, and validaton of the cranial ﬁnite element model. Med Eng Phys 20:607–619. Richmond BG. 2007. Biomechanics of phalangeal curvature. J Hum Evol 53:678–690. Ross CF, Patel BA, Slice DE, Strait DS, Dechow PC, Richmond BG, Spencer MA. 2005. Modeling masticatory muscle force in ﬁnite element analysis: sensitivity analysis using principle coordinates analysis. Anat Rec 283A:288–299. Roth S, Raul J-S, Willinger R. 2010. Finite element modelling of paediatric head impact: global validation against experimental data. Comput Meth Prog Biol 99:25–33. Schmidt H, Alber T, Wehner T, Blakytny R, Wilke H-J. 2009. Discretization error when using ﬁnite element models: analysis and evaluation of an underestimated problem. J Biomech 42:1926–1934. Strait DS, Richmond BG, Spencer MA, Ross CF, Dechow PC, Wood BA. 2007. Masticatory biomechanics and its relevance to early hominid phylogeny: an examination of palatal thickness using ﬁnite element analysis. J Hum Evol 52:585–599. Strait DS, Wang Q, Dechow PC, Ross CF, Richmond BG, Spencer MA, Patel BA. 2005. Modeling elastic properties in ﬁnite-element analysis: how much precision is needed to produce an accurate model? Anat Rec 283A:275–287. Tanner JB, Dumont ER, Sakai ST, Lundrigan BL, Holekamp KE. 2008. Of arcs and vaults: the biomechanics of bone-cracking in spotted hyenas (Crocuta crocuta). Biol J Linn Soc 95: 246–255. Tseng ZJ. 2009. Cranial function in a late Miocene Dinocrocuta gigantea (Mammalia: Carnivora) revealed by comparative ﬁnite element analysis. Biol J Linn Soc 96:51–67. Viceconti M, Olsen S, Nolte L-P, Burton K. 2005. Extracting clinically relevant data from ﬁnite element simulations. Clin Biomech 20:451–454. Wroe S. 2008. Cranial mechanics compared in extinct marsupial and extant African lions using a ﬁnite-element approach. J Zool 274:332–339. Wroe S, Clausen P, McHenry C, Moreno K, Cunningham E. 2007. Computer simulation of feeding behaviour in the thylacine and dingo as a novel test for convergence and niche overlap. Proc R Soc B 274:2819–2828.