One step further in biomechanical models in palaeontology: a nonlinear finite element analysis review

Non-linear materials: hyperelasticity and plasticity

Non-linear materials are materials in which the constitutive equation that defines their behaviours establishes a relationship between stress and strain that is not proportional to a constant. Typical material non-linearities can be found in phenomena such as plasticity and hyperelasticity. Plasticity describes the deformation of a material undergoing non-reversible changes of shape in response to applied forces. In a typical stress-strain curve for a plastic material there is a linear elastic region which satisfies Hooke’s law and a plastic region before fracture that can also follow a linear law or can be defined using different linear sections (Fig. 3). The transition from elastic behaviour to plastic behaviour is called yield and a non-linear solution is required because the solver needs to discover if the body is in the plastic region or not. The total strain is defined by ${\epsilon }_{total}={\epsilon }_{elastic}+{\epsilon }_{plastic}$ and the value of stress will depend on the value of this total strain. In a biological context, plasticity can be found in trabecular bone formulations to capture tension-compression asymmetry in the yield strength (Gupta et al., 2007) or, more generally, in studies where a permanent deformation or plasticity is assumed in cortical bone or other biological materials such as dentin, enamel or nacre (An, 2016). In materials such as dentine or enamel, plasticity represents the extremes of the loading environment rather than everyday behaviour.

A hyperelastic material is one that shows extreme elastic behaviour, in that it can return to its original shape even after experiencing very high strains. They are ideal elastic materials in which the stress-strain relationship is non-linear because derives from a strain energy function instead of Hooke’s law (Fig. 3). Moreover, these materials use the large deformation theory already mentioned above. However, the response of the material is not plastic because deformations are fully recoverable. Typical formulations of hyperelastic materials are, among others, phenomenological descriptions of observed behaviour in Mooney–Rivlin and Ogden formulations or equations describing the underlying structure of the material in the Neo–Hookean model (Ogden, 1984). Hyperelastic formulations are common in soft tissues such as ligaments or tendons (Shearer, 2015). Specifically, they can be found in the periodontal ligament (PDL) (Bucchi et al., 2019), muscles such as the pelvic floor (Stansfield et al., 2021), the abdominal muscle (Tuset et al., 2019) or generic muscular tissues (Hedenstierna, Halldin & Brolin, 2008), skin (Ito et al., 2022), corneas (Shan et al., 2010), cartilage (Pataky, Koseki & Cox, 2016), the temporomandibular joint (Sagl et al., 2019) or when modelling blood vessels (Vorp, 2007).

Sometimes the equations that are defined to control soft tissue behaviour include a viscous term (Huang et al., 2017). Viscoelasticity describes the variation of material response within time containing an elastic and a viscous part. The viscous part can describe creep, when stress remains constant and the deformation increases with time, or relaxation, when the deformation remains constant and stress decreases over time. On the other hand, the elastic response is instantaneous and can be defined using a linear material (Booker & Small, 1977) or a nonlinear hyperelastic material (Kulkarni et al., 2016).

More complex models, including fibres in their formulation, exist for the arterial vessels (Gasser, Ogden & Holzapfel, 2006) or the intervertebral discs (Noailly, Planell & Lacroix, 2011) among others. Despite the complexity of these formulae, which combines the overlay of the stiffness in the preferred directions of the fibres with the hyperelastic formulation of the matrix, the constitutive equation is also nonlinear, and it must be solved following an iterative procedure.

Non-linearities in geometry: buckling

In a linear problem, the equations of equilibrium are formulated in the original undeformed state and are not updated with the deformation. This is common in most engineering problems because the deformations are small enough to avoid differentiating the original geometry and the deformed one. However, there are cases where the deformation cannot be ignored, and we need to include large displacement non-linearities due to the geometrical update during the application of forces: This is the case of buckling.

Buckling implies a sudden change in shape of a body under load because the loss of stability when this load reaches certain critical value (Fig. 4). If a body- such as a column under compression or a plate under shear, for example- is subjected to a gradually increasing load, when the load reaches the critical value, the body may suddenly change shape. Although buckling appears before failure, it can be decisive in the ergonomics of certain biological bodies, limiting the range of forces under which they are able to remain functional. Buckling is caused by nonlinearities in the geometry and can be approached by a linearisation that drives to a bifurcation problem of eigenvalues. Therefore, the linear buckling analysis is done in parallel to a linear elastic analysis. Otherwise, the full nonlinear solution of the point of collapse can be obtained by increasing the load in smaller steps with an iterative method while the geometry is updated to its deformed state. This latter is significantly more computationally expensive but might be more accurate than the linear buckling. In a biomechanical context, buckling can be found when study slender bodies such as the swordfish rostrum (Habegger et al., 2020), the weevil rostrum (Matsumura et al., 2021) or even in bones under compression such as the vertebrae (Williams et al., 2021). Buckling is also considered in humans as a cause of fractures of postcranial bones (Lee et al., 2009).

Non-linearities in contacts

Contacts between two bodies are divided between linear contacts and non-linear contacts. Linear contacts can be included in a linear elastic model without modifying the solving mode and keeping the direct solution. It also involves a low computational cost that is simply a function of the size of the finite element mesh (namely, the number of elements and nodes). However, the inclusion of non-linear contacts changes the solving mode to a non-linear solution with an iterative solver, increasing the computational cost of the analysis. Contacts can be described according to the relationship between the two separate surfaces of each body that become mutually tangential in five general different types according to how they can move perpendicularly to each other and how they can move in the tangential plane. In other words, if they are allowed to separate and slide (Fig. 5).

1. Bonded contacts: when separation and sliding is not allowed. It is a linear contact.

2. No-separation contact: when separation is not allowed but sliding in the tangential plane is allowed. It is a linear contact.

3. Frictionless contact: when separation and sliding is allowed. It is a non-linear contact.

4. Rough contact: when separation is allowed but sliding in the tangential plane is not allowed. It is a non-linear contact.

5. Frictional contact: when separation is allowed but sliding in the tangential plane is controlled by a friction coefficient. It is a non-linear contact.

Frictional contact can be understood as an intermediate status, where sliding in the tangential plane is not free but is allowed and bonded contact is used when we have two bodies that are perfectly joined but they are created or defined as separate bodies during the FEA modelling. For example can be used for defining all the pieces involving a teeth such as the cortical bone, dentine, enamel, pulp and the PDL (Benazzi et al., 2013; Bucchi et al., 2019)

In general, contacts are found in FEA models involving more than one body and the definition of each contact depends on the nature of its behaviour. It can be found in models when studying the carpal bones of the wrist (Gíslason et al., 2017; Püschel et al., 2020a) or the feet (Ito et al., 2022), the ossicles of the auditory system (Marcé-Nogué & Liu, 2020), the intervertebral discs and the vertebrae of the spine (Guan et al., 2019), all the tissues in the hip (Fleps et al., 2018) or the patella (Fitzpatrick & Rullkoetter, 2012), the mandible, the tempomandibular joint and the skull (Sagl et al., 2019) or the interaction between the bodies in the wings of dragonflies (Rajabi et al., 2016b), bees and wasps (Eraghi et al., 2021), among many others. Therefore, contacts can be used to establish relationship between bones or soft tissues. Contacts are also useful when studying occlusal forces during mastication to model the interaction between teeth and food (Skamniotis, Elliott & Charalambides, 2019) or even the impact of eggshells with the floor (Sellés et al., 2019).

Non-linear soft tissues

Little is known about soft tissue properties in fossils. The direct examination of fossil soft tissues and preserved blood cells is of little value when studying palaeontological remains due to the degradation or the contamination from modern remains (van Dongen et al., 2017). The reconstruction of soft tissues from fossils is an issue that it is unresolved but can be approached through investigating extant relatives to infer the palaeo-physiology of extinct taxa, e.g., via the phylogenetic bracketing approach (Witmer, 1995). Therefore, any FEA models can potentially include an inference regarding soft tissues properties. As an example, cranial sutures are deformable joints between adjacent bones bridged by collagen fibres and there are several works on fossil taxa that have include soft tissues modelling sutures in Tyrannosaurus rex skull (Cost et al., 2020), australopithecines (Dzialo et al., 2014) or in dicynodonts (Jasinoski, Rayfield & Chinsamy, 2009), as well as FEA models of current lizard species (Dutel et al., 2021), Sphenodon (Curtis et al., 2013) or even some mammals (Bright & Gröning, 2011). All of these examples used linear material properties to characterize the elastic behaviour of soft tissues which can be an appropriate simplification if this is validated experimentally (Bright & Gröning, 2011). However a recent diagnosis suggested that the lack of sutures or and inappropriate modelling can result in inaccurate results of stress, strain or deformation (Rayfield, 2019) although it is not clear how the soft tissue can be accurately predicted in fossils (Broyde et al., 2021). It is at this point that the researcher needs, at least, to be aware that a more accurate modelling of these soft tissues should be done using nonlinear material properties, which in turn involve an increase in the computational cost of the model. Unfortunately, soft tissue models in living animals have not been extensively documented, with the exception of biomedical contexts (Tanaka et al., 2000). Hopefully, studies focused on experimentally testing non-linear models in diverse biological taxa will be carried out in the near future.

Plasticity in retrodeformations

Retrodeformation is very common in fossil taxa as the process that produces the original form of the taxon prior to fossil diagenesis when this has been recovered in any deformed way. Deformation in fossils is produced due to a multiple array of taphonomic and tectonic processes. Overburden stress due to the weight of the overlying sediments linearly compacts the fossil from above causing the fossil to break and/or warp. Other causes of fossil deformation include tectonic stresses and sediment cracking. Under the action of these loads, the fossil can break in a brittle manner or can be distorted plastically, preserving the structure of the fossil due to the lack of breakage. Fossils under plastic deformation, where forces applied during time modify the original shape of the bone structure may be restored. Although there are several techniques to virtually restore deformed specimens available without using mechanical equations (Lautenschlager, 2016), it has sense to use methods from mechanics such as FEA that involve forces if one want to infer which was the actual process that drove the fossil to be deformed (Arbour & Currie, 2012; Di Vincenzo et al., 2017). Modelling this phenomenon would require including the nonlinear plastic behaviour of bone, because retrodeformation is a permanent deformation in cortical bone. In this case FEA could be inversely applied by defining the plastic behaviour of the fossil material and then, setting the forces applied at the fossil as the unknowns of the problem. It would allow to answer the question of which forces do we need to apply deformed bodies to recover its original form.

Buckling in slender bones

In palaeontology there are a lot of slender structures that are susceptible to be analysed using buckling. Probably the most common and useful case would be in bones under compression such as the leg bones of large, heavy dinosaurs and mammals. This is because mass is considered as one of the main factors affecting the morphology and osteological adaptation of these bones (Etienne et al., 2020). To understand how these bones are adapted to the heavy weight that they needed to support, evaluation of the maximal stress as a measure of bone strength is not the only informative metric (Hutchinson, 2021). This is because bone may fail without involving fracture. Alternatively, bone could fail through buckling if it is not stiff enough (Currey, 2012). In this case, buckling needs to be considered, because it can cause the collapse of the legs before the fracture of the bone. Usually, buckling reduces the capacity of the strength of the structure because it appears in a lower value than the yield stress and the fracture stress that defines the strength of the material.

If we assume that the leg bones in heavy dinosaurs as slender columns like those from a building, Euler’s critical load is defined as the compressive load at which the column will suddenly buckle (Timoshenko, 1955). This equation can give clues about the relationship between geometrical factors such as the length of the bones or how they are joined to the articulations. Given that the length, material, or boundary conditions cannot be modified from the original model, Euler’s critical force will depend on the second moment of area or moment of inertia. Increasing the value of the critical force implies a modification of the cross-section of the bone through more inertial geometries. Therefore, if we assume the cross-section of leg bones as an annulus, thicker annulus will increase the inertia. But also, if the thickness is kept constant, a broader annulus will increase the inertia of the cross section. This simple consequence can be obtained assuming leg bones with a straight morphology not close enough to the reality, but very useful for the purpose of study. However, in case of analysis of the real and irregular geometry of the bones, the simple formula of Euler cannot be used but the problem of buckling can be solved via computational methods by means of FEA solving an eigenvalue problem. Few works are paying attention to it, discarding the effect of buckling in the morphology of the long bones in living mammals (Brassey et al., 2013). Considering than an eigenvalue problem in a FEA model is not increasing the computational cost of the analysis too much, it would be worth to more exhaustively test if the leg bones of heavy dinosaurs or mastodontic extant mammals are affected by buckling, as was suggested in horses (Currey, 2003). Another case where buckling could be a concern is in bones with a high aspect ratio where the walls are substantially thinner (De Margerie et al., 2005). This could be the case of wing bones (Palmer & Dyke, 2012), and hence buckling should be also explored.

Bone grouping using contacts

Functional implications of fossil bones have been widely studied in fossil taxa using FEA models (Richmond et al., 2005). Depending on the purpose, bones can be studied alone or as a group and the main difference between these two cases is the absence or presence of contacts. When separation between bones is not desired, for example in the analysis of teeth, considering the bonding of the cortical bone, dentine, enamel, pulp and the PDL (Benazzi et al., 2013), the contacts used are linear and it does not require an increase in the computational cost of the solving process. This is something that can be considered when creating FEA models because it allows the inclusion of several bones in the model without nonlinearities.

On the other side, nonlinear contacts allow separation between the bones. Although the inclusion of this contacts implies an iterative solution through convergence, it may be necessary to implement when a group of bones need to be studied together such as the carpal bones of the wrist (Gíslason et al., 2017) or the bones of the foot (Ito et al., 2022). This has been done extensively in biomechanical models of living primates; therefore, it should be considered in other FEA works in the field of the palaeontology and anthropology. In fact, the literature is full of biomechanical analysis of kinematics and dynamics of solid bodies where bones of fossils are grouped to study its performance (Sellers et al., 2017; Bishop et al., 2021). Therefore, it makes sense when creating FEA models, to include more than one bone in the model if it can be useful for the desired analysis despite increasing the computational cost of the solution. Also when the contact between bones is through articular cartilage, the contact can be defined between cartilages that are also in contact with the bone (Püschel et al., 2020a).

Models with shells, plates, beams, springs, and trusses

Finally, although this is not related with the use of a nonlinear iterative solving, the use of other kind of elements other than solid elements can have a great advantage when dealing with nonlinear models. This is because they provide a useful way to reduce the number of elements and nodes of the FEA mesh and, consequently, a reduction of the time spent solving the equation in each iteration. Hence, a nonlinear model will particularly benefit from the use these elements.

The use of shell elements to model cortical bone in morphologies that can be assumed as thin and with a constant thickness, such as carpal bones or talar morphologies, require a lower number of elements and nodes because there is only one mesh layer. Using solid elements in the same morphology would need at least four or five layers of elements would along the thickness to properly build an adequate mesh to accurately capture the results. This was used in an analysis of carpal bones (Gíslason et al., 2017) to model both the cortical bone and the articular cartilage, reducing significantly the number of elements to allow a smooth non-linear solution, due to the presence of non-linear contacts. The same example uses non-linear spring elements to model the behaviour of the ligaments. This decision is also in favour of not increasing the number of nodes and elements of the model, because spring or truss elements can be defined using only one element with the origin and final nodes. In this manner, the model avoids the inclusion of a three-dimensional geometry modelled with solid elements for each ligament, which would exponentially increase the number of nodes and elements in the mesh and consequently, increase the computational cost of the solution.

When creating FEA models of fossils and considering the inclusion of some of the non-linearities it is a good option to evaluate if the use of simpler elements can reduce computational cost. Although researchers should be aware of the potential ramifications of simplifying their models, it is also true that any model will necessarily not represent a literal representation of reality. Instead, the requirements necessary to answer the research question of interest should always be kept in mind when making decisions about model complexity.

Summarizing, this text highlighted the usefulness of non-linearities in FEA palaeontological and anthropological models in spite of increasing their complexity and the computational costs. Nowadays, most of the commercial and non-commercial FEA packages include the resolution of non-linear problems in their capabilities, and they also documented with tutorials and examples how to deal with them. Therefore, the main aim of this review is to provide a road map for the next generation of palaeontologists, anthropologists, and functional morphologists by showing them unexplored ways that could a profound impact in finite element analysis and how they can explore these methods further.