Computational fluid dynamics as a tool for testing functional and ecological hypotheses in fossil taxa

Computational fluid dynamics is a method for simulating fluid flows that has been widely used in engineering for decades, and which also has applications for studying function and ecology in fossil taxa. However, despite the possible benefits of this approach, computational fluid dynamics has been used only rarely in palaeontology to date. The theoretical basis underlying the technique is outlined and the main steps involved in carrying out computer simulations of fluid flows are detailed. I also describe previous studies that have applied the method to fossils and discuss their potential for informing future research directions in palaeontology. Computational fluid dynamics can enable large‐scale comparative analyses, as well as exacting tests of hypotheses related to the function and ecology of ancient organisms. In this way, it could transform our understanding of many extinct fossil groups.

Abstract: Computational fluid dynamics is a method for simulating fluid flows that has been widely used in engineering for decades, and which also has applications for studying function and ecology in fossil taxa. However, despite the possible benefits of this approach, computational fluid dynamics has been used only rarely in palaeontology to date. The theoretical basis underlying the technique is outlined and the main steps involved in carrying out computer simulations of fluid flows are detailed. I also describe previous studies that have applied the method to fossils and discuss their potential for informing future research directions in palaeontology. Computational fluid dynamics can enable large-scale comparative analyses, as well as exacting tests of hypotheses related to the function and ecology of ancient organisms. In this way, it could transform our understanding of many extinct fossil groups. R E C O N S T R U C T I N G how ancient organisms moved and fed is key to understanding past ecosystems, but is often dismissed as speculative because such inferences can be hard to test. In particular, this work has long been hindered by a lack of objective data on the function of extinct species, especially those without extant analogues. However, the recent development and increasing availability of techniques for visualizing and analysing specimens digitally and in three dimensions, including virtual modelling approaches, provides a quantitative framework for testing specific hypotheses even in problematic fossil taxa. Thus, it is now more feasible than ever before to study the function and ecology of ancient organisms.
One such method for palaeontological functional analysis is computational fluid dynamics (CFD). CFD is a tool for simulating flows of fluids (liquids or gases) and their interaction with solid surfaces. Complex equations describing the motion of fluids are solved numerically using a computer, and the results can be visualized as plots of fluid properties (e.g. velocity and pressure) within the flow domain. The first computational simulations of fluid flow were undertaken in the 1950s and 1960s (e.g. Evans & Harlow 1957;Harlow & Welch 1965;Hess & Smith 1967), fuelled by the development of increasingly powerful computers. CFD has subsequently been used to analyse design and optimize performance for a wide variety of structures and machines, ranging from nuclear reactors (e.g. Yadigaroglu 2005;Mahaffy et al. 2007) to aircraft (e.g. Agarwal 1999;Spalart & McLean 2011) and it is well established in engineering. CFD is also becoming increasingly important in biology and medicine; for example, to study blood flow in the human cardiovascular system (e.g. Migliavacca et al. 2000;Hoi et al. 2004;Young et al. 2014) and the functional performance of flying and swimming animals (e.g. Liu 2002;Kato & Kamimura 2008;Deng et al. 2013). To date, however, it has been employed only rarely in palaeontology (e.g. Rigby & Tabor 2006;Shiino et al. 2012;Rahman et al. 2015a).
In this paper, I introduce the fundamentals of fluid dynamics and outline the main steps involved in conducting a CFD analysis. I also discuss the applications of CFD to the study of fossil material, encompassing work on extinct vertebrates and invertebrates, as well as different types of flows. This approach has great potential for informing rigorous tests of long-standing functional and ecological hypotheses, and could thus impact on a variety of research directions in palaeobiology.

FLUID DYNAMICS
CFD is an approach in fluid dynamics, which is the branch of fluid mechanics concerned with fluids in motion (as opposed to fluid statics, which deals with fluids at rest). Fluid dynamics relies on three conservation laws of physics: the conservation of mass; the conservation of momentum (Newton's second law); and the conservation of energy. These can be formulated as mathematical equations (such as the Navier-Stokes equations that represent the conservation of momentum and the continuity equation that represents the conservation of mass) which are solved to predict fluid flow. The specific equations used depend on the properties of the fluid; for example, the Navier-Stokes equations are different for compressible (changes in pressure cause changes in density) and incompressible (changes in pressure cause no changes in density) fluids. Liquids are usually treated as incompressible because any density changes are so small as to be negligible; gases are considered as compressible or incompressible according to the flow velocity (compressible for Mach numbers greater than about 0.3).
The behaviour of the governing equations depends on the contribution of the forces acting on the fluid. These include the inertial forces, which are related to the momentum of the fluid, and the viscous forces, which are frictional forces (shear stress) due to the relative motion of adjacent layers of the fluid. Viscosity is the resistance of the fluid to deformation by these frictional forces. The fluid is said to be Newtonian when the viscosity is constant and non-Newtonian when the viscosity is dependent on the shear rate; air and water are examples of Newtonian fluids, while blood is a non-Newtonian fluid. The velocity of the flow decreases in the immediate vicinity of a solid surface due to increased frictional forces, producing a velocity gradient that reaches zero at the surface. This thin layer of fluid is called the boundary layer.
The dimensionless Reynolds number (Re) is used to characterize the flow and is defined as the ratio between inertial forces and viscous forces: where q is the density of the fluid (kg/m 3 ), U is the characteristic velocity (m/s), L is the characteristic dimension (m) and l is the dynamic viscosity of the fluid (kg/ (m·s)). When viscous forces are important (i.e. the viscosity is high, the velocity is low and/or the dimension is small), the Reynolds number will be low and the flow laminar, meaning that the fluid flows orderly in parallel layers (Fig. 1A). As inertial forces begin to dominate (e.g. with increasing velocity), the Reynolds number increases and the flow becomes turbulent, characterized by the formation of eddies and chaotic motion (Fig. 1B) and a steeper velocity gradient in the boundary layer. In laminar flow, there is little mixing between fluid layers, while in turbulent flow, considerable mixing occurs. The exact values, or critical Reynolds numbers, over which the flow regime changes from laminar to turbulent vary on a case-by-case basis, and there is a broad transitional zone in which both types of flow can occur. Laminar flows can be described through the Navier-Stokes equations, but it is computationally infeasible to resolve the spatial and temporal variation of most turbulent flows using these equations. Typically, an alternative approach is required using, for example, the Reynolds-averaged Navier-Stokes (RANS) equations-time-averaged equations for fluid flow-and a mathematical model to predict the effects of turbulence. Lastly, the effects of time on the flow fluid should be considered. Where the fluid properties are invariant with time, the flow is said to be steady (or stationary). However, if the fluid properties vary with time, the flow is unsteady (or transient). Laminar flows are often steady, whereas turbulent flows are always unsteady. Nevertheless, turbulent flows can be treated as statistically stationary if the mean flow properties do not change with time (Pope 2000). Steady flows are more tractable than unsteady ones, which must incorporate time into calculations and, hence, take significantly longer to solve.

HOW DOES CFD WORK?
There are a variety of commercial (e.g. ANSYS Fluent, Autodesk CFD and COMSOL Multiphysics) and open- source (e.g. Fluidity, OpenFOAM and SU2) software packages for CFD. All these programs share a number of common steps which are outlined below.

Computational domain
The first step is to create the computational domain, which consists of any objects around/through which the fluid will flow ( Fig. 2A). In many cases, these objects are imported from separate computer-aided design (CAD) programs, but simple shapes can be built in the CFD software itself. The dimensionality of the geometry will depend on the objective of the analysis and the available computer power; one-and two-dimensional models are much less computationally expensive than three-dimensional ones, but might not be appropriate for complex shapes that are hard to define with a single line or crosssection, as is often the case for biological and palaeontological specimens. Three-dimensional digital reconstructions are now routinely produced in studies of fossil morphology (Abel et al. 2012;Cunningham et al. 2014;Sutton et al. 2014Sutton et al. , 2017, and in many cases these 'virtual fossils' will be suitable for use in CFD analyses, albeit after correcting for taphonomic and diagenetic distortion (e.g. Lautenschlager 2016), removing noise and other artefacts and conversion into an appropriate file format for the CFD software. When flow is to be simulated around the exterior of the fossil, internal features are unimportant, and thus reconstructions could be derived from tomographic (e.g. X-ray computed tomography) or surface-based methods (e.g. laser scanning or photogrammetry); however, internal flows (e.g. Bourke et al. 2014) require reconstructions that incorporate interior features, meaning that surfacebased techniques will be inappropriate. Alternatively, three-dimensional models of fossils can be created through box modelling, which is particularly useful for taxa lacking complete, three-dimensionally preserved specimens (Rahman & Lautenschlager 2017). A plethora of software is available for digitally reconstructing or modelling fossils in three dimensions (Abel et al. 2012;Cunningham et al. 2014;Sutton et al. 2014Sutton et al. , 2017. In simulations of external flows, the domain surrounding the fossil must also be defined, and this will typically comprise a cuboid or cylinder (for three-dimensional models) ( Fig. 2A). The size of the domain should be sufficiently large to allow the flow to fully develop on all sides of the fossil. For example, previous studies (e.g. Shiino et al. 2009Shiino et al. , 2012Shiino et al. , 2014 have used a computational domain that extended three times the length of the fossil upstream, ten times the length of the fossil downstream and five times the size of the fossil in all other directions; however, in some cases it might be possible to accurately resolve the flow using a smaller domain. Sensitivity analyses can be conducted to determine the optimal domain size (see below).

Material properties and flow model
After creating the computational domain, the physical properties (e.g. density and viscosity) of the fluid can be assigned to the model using a built-in materials library or based on values taken from the literature (e.g. Poling et al. 2000). An appropriate flow model is then selected depending on whether the flow regime is characterized as laminar (low Re) or turbulent (high Re). The available flow models (i.e. laminar and different turbulence models, such as the Spalart-Allmaras, k-epsilon and SST models) will differ according to the software used; models vary in terms of the calculations performed to simulate the flow, and choosing between them is partly a compromise between accuracy and computational time/memory. See Wilcox (2006) for an in-depth discussion of alternative turbulence models.

Boundary conditions
Boundary conditions representing flow variables at the boundaries of the domain must be added to the model ( Fig. 2A). Typically, an inlet boundary condition is defined at one end of the computational domain, representing net flow into the domain, with an outlet boundary condition at the opposing end of the domain (net outflow). A common configuration for incompressible flows has the velocity specified at the inlet and the pressure fixed at zero at the outlet. The inlet velocity will depend on the palaeobiology and inferred environment of the fossil taxon; in simulations of organisms moving through the flow, comparisons with theoretical or experimentally-derived flying or swimming speeds of extant animals (e.g. Azuma 2006) could provide a realistic range, whereas for stationary organisms, data on water or wind currents in modern environments (e.g. Siedler et al. 2013;Randall 2015) should be consulted. A wall boundary condition describing the fluid velocity at the fluid-solid interface is assigned to all solid surfaces; in most cases (laminar and many turbulence models), a no-slip boundary condition is used, which constrains the fluid velocity at zero relative to the object. For external flows, a slip boundary condition can be assigned to the edges of the domain surrounding the fossil, allowing the flow to pass along the walls without friction.

Discretization
In most cases, the next step is to divide the computational domain into a finite number of discrete cells (i.e. a mesh; Fig. 2B, C). Mesh generation can be automated in many CFD programs, but manual refinement might be necessary for complicated structures. For three-dimensional models, meshes are made up of simple shapes like tetrahedra or hexahedra, which vary in number and size (including within the same mesh) according to the complexity of the geometry and the flow model; thin layers of prismatic elements can be inserted along the interface between the fluid F I G . 2 . Main steps in a CFD study using the cinctan Protocinctus mansillaensis as an example. A, threedimensional plot of computational domain with boundary conditions marked. B-C, two-dimensional plots (vertical cross-sections) of mesh used in analysis (red dashed box in B marks the position of C). D, twodimensional plot (vertical cross-section) of velocity magnitude with streamlines from simulation of water flow. and the solid (Fig. 2C) to better capture flow in the boundary layer. The governing equations of fluid flow are then discretized so they can be solved on each of the mesh cells. There are a number of different techniques for discretization, the most common of which are the finite-difference, finite-element and finite-volume methods. The CFD software selected governs the discretization method; for example, ANSYS Fluent and OpenFOAM use the finite-volume method, while Autodesk CFD and COMSOL Multiphysics use the finite-element method. Meshless approaches also exist (e.g. smoothed particle hydrodynamics; Monaghan 2005), but are not as widely used at present.

Solver
After completing the steps outlined above, the discretized equations can be solved, and hence fluid flow simulated, using the chosen CFD software. The time required for this will depend on a number of factors, including the complexity of the geometry, the flow model and the mesh, as well as the available computer power. The solver type is also important. When the mean properties of the flow do not change over time (i.e. steady flow), a stationary solver can be used to compute the steady-state solution. In contrast, when the flow properties vary with time (i.e. unsteady flow), a time-dependent solver should be used, which computes the solution at select time steps and therefore requires considerably longer than a stationary solver.

Visualization and post-processing
The results of CFD simulations can be visualized in a variety of different ways. Commonly, plots of flow velocity, pressure or vorticity are created; for time-dependent simulations, these can be generated for different time steps and combined to produce animations. Plots can be supplemented with, for example, streamlines tracing the movement of particles (Fig. 2D). It is also possible to evaluate the forces exerted by the fluid on the fossil, such as the drag (parallel to flow direction) and lift (perpendicular to flow direction) forces. Dimensionless coefficients of drag (C D ) and lift (C L ) can be calculated to facilitate comparisons between models. For three-dimensional models, these are defined as: where F D is the drag force exerted by the fluid (N), F L is the lift force exerted by the fluid (N), q is the density of the fluid (kg/m 3 ), U is the characteristic velocity (m/s) and A is the characteristic area (m 2 ). Depending on the geometry and the flow conditions, the characteristic area is normally the projected frontal area or the surface area.

EXPERIMENTAL COMPARISON AND SENSITIVITY ANALYSIS
CFD is typically cheaper, faster and more flexible than laboratory experiments in flumes or wind tunnels, making it an attractive approach for studying fluid flows. However, comparisons with experimental data are necessary to establish how accurately the computer simulations replicate real-world flow patterns. Studies of this sort are widespread in engineering, especially in the aerospace and automotive industries. They often entail comparing the results of experiments and CFD for benchmark models (e.g. air flow over the Ahmed body, see Ahmed et al. 1984;Lienhart & Becker 2003), but models of more complex geometries can also be utilized (e.g. 3-D printed fossils; Rahman et al. 2015a). The results can be compared qualitatively, or a metric such as the drag coefficient can be used to quantify the extent to which the two methods agree.
Testing the sensitivity of results to changes in the simulation parameters is another important part of a CFD study. The mesh size, in particular, can have a substantial impact on the accuracy of the solution; meshes composed of numerous smaller hexahedral cells with multiple thin prismatic layers at the fluid-solid interface tend to more accurately represent the flow, but require longer simulation times and more computational memory than coarser meshes. Consequently, a sensitivity analysis should be carried out in which simulations are repeated using successively finer meshes in order to establish the coarsest possible mesh size that guarantees the computational accuracy of the CFD simulation, while also minimizing the simulation duration. In certain cases, it will also be necessary to carry out sensitivity analyses to test the influence of other parameters, such as the domain size, flow model or solver type, on the CFD results.

EXAMPLES IN PALAEONTOLOGY
CFD has much potential for testing hypotheses pertaining to an organism's functional performance, but compared with computer modelling methods like finite-element analysis (e.g. Rayfield 2007; Bright 2014) it has been used infrequently in palaeontology thus far. Rigby & Tabor (2006) were among the first to apply CFD to fossils, simulating laminar flow around two-and three-dimensional models of graptolites (Fig. 3A). Subsequent studies have instance, Shiino et al. (2012) tested the hypothesis that the long, forked hypostome of trilobites aided swimming performance. In this study, they ran simulations of flow around three-dimensional models of the trilobite Hypodicranotus striatulus with and without a hypostome (Fig. 3B). The results showed that the model with a hypostome generated generally lower drag coefficients and more stable lift coefficients than the model without a hypostome, and this suggests that the hypostome played an important role in active swimming. Rahman et al. (2015a) used CFD to test between competing hypotheses for the feeding mode of cinctan echinoderms. They ran simulations for a three-dimensional model of the cinctan Protocinctus mansillaensis using different boundary conditions at the mouth to represent either passive or active suspension feeding (Fig. 3C). These simulations showed that there was almost no flow to the mouth in simulations of passive feeding, whereas considerable flow towards the animal was observed for active feeding, supporting the idea that cinctans were active suspension feeders. Liu et al. (2015) tested between alternative theories for the locomotory mode of plesiosaurs. Using a three-dimensional, articulated model of the plesiosaur Meyerasaurus victor, they conducted simulations of swimming assuming different ranges of motion for the forelimbs and hindlimbs (Fig. 3D). The results suggested that efficient forward motion entailed using the forelimbs to provide the majority of thrust, with the hindlimbs providing only weak thrust.

FUTURE PROSPECTS
CFD has long been an established method in engineering and has developed as a valuable tool in medicine and biology over the past two decades, but it is only in the last 5-10 years that palaeontologists have begun to make use of the approach. In part, this lack of uptake probably reflects the perceived difficulty in applying CFD to fossils; in this paper, I have sought to address this by outlining the main considerations and key steps involved in performing such analyses. Another obstacle is the cost of the software and hardware needed for CFD, but this can be offset through collaborations with colleagues already working in the field, as well as by making use of opensource software and high-performance workstations optimized for visualizing tomographic data (which are increasingly common in palaeontology research groups). Finally, some researchers might be concerned that simulations do not accurately capture real flow patterns around living, moving organisms. This could be addressed through comparisons with laboratory experiments (see above) and by carrying out computer simulations of articulated models (e.g. Liu et al. 2015) to evaluate the interaction between solid objects in motion and the surrounding fluid. Consequently, the major barriers to using CFD in palaeontology can now be overcome.
The overwhelming majority of palaeontological CFD studies have simulated flows of water around the exterior of fossil organisms (although see Bourke et al. 2014 for one exception). Thus, there is great potential to use CFD to analyse internal fluid flows (e.g. air or blood), as well as for studying the performance of extinct flying organisms, such as pterosaurs. Another area of considerable promise is in using articulated models to explore how flow patterns can vary during locomotion, as pioneered by Liu et al. (2015) in their study of plesiosaur swimming. Such work will enable much more realistic simulations of swimming and flight in fossil taxa, and might even have applications in robotics or the aerospace and automotive industries (i.e. palaeobiomimetics). Lastly, comparative analyses of different fossil taxa, including hypothetical morphologies (e.g. Shiino & Kuwazuru 2010;Shiino et al. 2012;Rahman et al. 2015b), should prove especially informative for addressing outstanding questions in palaeobiology, and will enable thorough tests of evolutionary scenarios. This work is increasingly feasible now that methods for the three-dimensional visualization of fossils have become a mainstay of palaeontology (Cunningham et al. 2014;Sutton et al. 2014Sutton et al. , 2017. One outstanding challenge for researchers using CFD is establishing how best to make the raw data underlying their work available. Most CFD software generates simulation results in proprietary file formats, which will not be accessible to all potential users. This represents a barrier to the verification and reproducibility of CFD-based studies, while also limiting the possible reuse of datasets arising from this work (e.g. in large-scale comparative analyses). Ideally, software developers will devise ways to make full details of the computational domain, fluid properties, flow model, boundary conditions and meshes available in accessible file formats. In the meantime, researchers should ensure that this information is reported in the methods sections of all relevant publications (Davies et al. 2017).
CFD is an extremely powerful approach that can be used to rigorously evaluate functional and ecological hypotheses in ancient organisms, opening up an exciting new avenue in palaeontological functional analysis. Thus far, studies have focused predominantly on analysing water flows around static models of Palaeozoic marine invertebrates, but there is enormous potential for investigating feeding and movement in a wide range of groups in both aerial and aquatic environments. Further work comparing the results of CFD to laboratory experiments will be important for establishing the accuracy of these computational analyses. In the coming years, I anticipate that CFD will become an important part of the palaeontologist's toolkit for addressing questions related to the function of fossil species, transforming our understanding of how ancient organisms moved and fed.
Acknowledgements. I thank Andrew Smith for inviting me to write this review article for the Frontiers in Palaeontology series. I am grateful to Stephan Lautenschlager, Susana Gutarra D ıaz, Sally Thomas and two anonymous reviewers, whose comments greatly improved this article. This work was funded by the Royal Commission for the Exhibition of 1851 and the Oxford University Museum of Natural History.