Introduction

Over the course of human history, the details surrounding many seminal events have been lost to time. The invention of the wheel is an example of one such episode about which we know very little. The wheel has been described as the most important mechanical invention of all time [1] and has been credited with creating seismic social and economic shifts in the trajectory of out of alignment, which necessitated periodic manual readjustment, as has been reported with the use of bilateral rollers [12]. So while the adoption of unilateral rollers imposed some cost in the form of increased friction, this strategy conferred significant advantages over free rollers, since it eliminated the need for periodic realignment and replenishing of rollers. Note also that when operating inside a tunnel, spent rollers could not simply be rolled around to the front of the cart but rather they would have to be lifted and carried or dragged, which would require additional energy expenditure.

Since unilateral rolling reintroduces a friction surface, it may appear as though this would be the mechanical equivalent of simply dragging the cart on skids, which was a commonly used method of transport throughout the ancient world [10,13]. However, bilateral rolling offers two significant advantages over dragging. First, by transferring the friction surface from the roller-ground interface to the roller-cart interface, one would have control over both frictional surfaces. While cart operators would have been limited in their ability to process the surface of the ground to reduce roughness and friction, they could easily process the surfaces of both the roller and the cart in order to reduce the coefficient of friction. The second advantage of limiting friction to the wheel-cart interface is that the cart operators could now apply lubricant to the friction surface, thereby significantly reducing the friction force. For this reason, the presence of grooves on the roller surface is highly significant. Grooves would have provided a channel in which to apply lubricant where it could remain free from contact with the ground. Whereas grooveless rollers would summarily deposit most of their lubricant on the ground after a few revolutions, a grooved roller could retain its lubricant for a much longer duration.

Wheel historians have long theorized about the use of grooved rollers [1,14]. In his 2011 treatise on the history of rotating machinery, Rao surmises that grooves initially formed inadvertently from the indentation of the rollers caused by narrow sledges that sat between the rollers and the heavy cargo being transported [1]. In his 2016 book on the history of the wheel, Bulliet observes that this theory is ubiquitous, but as yet unsubstantiated by archaeological evidence [14]. Whether their introduction was incidental or deliberate, these grooves set in motion the next stage of the evolution, which ultimately led to the familiar wheel-and-axle configuration. We propose that following the introduction of grooved unilateral rolling, there was an expansion and coalescing of the roller's grooves to form a single channel in the centre of the roller (see figure 1c). It is possible that this shape change was initially motivated by the desire for clearance so that the cart could easily pass over small objects obstructing its path.

By creating this large channel in the centre of the roller, designers effectively forged a central axle. In addition to providing clearance, this modification yielded a lighter more portable roller, while preserving structural stiffness. But the most significant benefit of the narrow axle was mechanical advantage. Because of this principle, an axle whose diameter was one-tenth that of the wheel would require roughly one-tenth the pushing force needed for a prismatic roller with no axle. Our theory differs from earlier theories in that previous authors did not address the transition to unilateral rolling, nor did they relate the evolution of the axle to mechanical advantage.

In §4, we derive an equation that quantifies mechanical advantage and describes the required pushing force as a function of the axle radius. Next, we simulate the evolution of the wheel-and-axle structure from a unilateral roller using a structural optimization algorithm in which the objective function (also known as the cost function) is the mechanical advantage formula we previously derived. In keeping with the forensic nature of our study, the investigation begins with the wheel-and-axle geometry as the starting point. We then attempt to identify a set of boundary conditions, a design objective and an initial baseline design (i.e. the initial roller design from which the wheel evolved) that would cause the algorithm to independently arrive at the wheel-and-axle geometry. Therefore, the algorithm, which comprises just one part within our broader investigation, has no explicit prior knowledge of the wheel-and-axle concept. Instead, it is provided only with a mathematical model with which to calculate how each candidate design would perform in terms of effort required to push the cart. From this information, the algorithm converges upon the familiar wheel-and-axle design, as shown in figure 2. The algorithm is iterative in that it makes a series of incremental updates to the design with each iteration. This result illustrates how the wheel-and-axle could have gradually evolved through a series of successive incremental improvements.

The output of the optimization algorithm is a monolithic structure in which the axle and wheels turn in unison. This version of the wheel-and-axle is referred to as a wheelset, and it is believed that this design preceded multi-body wheel-and-axle systems in which the wheels turn independently of the axle (see figure 1d) [15]. The monolithic nature of wheelsets meant that, compared with their multi-body counterparts, they performed poorly during turns, making it difficult to change direction. Whereas a multi-body design would allow one wheel to rotate faster than the other to cover more ground during the turn, turning on a wheelset would cause one wheel to drag along the ground, producing additional friction. A four-wheeled wagon using two wheelsets would have an even harder time negotiating turns. Here again, the context of the mine plays a crucial role. Since the mine passages were not natural geographic features, they could be made straight, thereby significantly reducing the need to turn and accommodating the use of four-wheeled carts.

Another disadvantage of the wheelset is its increased susceptibility to damage and failure. Using computer simulations, we show that in a wheelset, the wheel-axle junction experiences high stresses due in part to the torsional loads transmitted from the wheel to the axle. These loads are reduced in the multi-body wheel-axle system, which can make it more robust to damage. This factor combined with its improved manoeuverability likely led to the invention of the multi-body wheel-and-axle some 500 years after the initial invention of the wheelset [14]. This development marks the third innovation in the sequence that defines the early evolution of the wheel.

Historical and archaeological evidence for the Carpathian-roller hypothesis

Archaeologists excavating sites in the Carpathian Mountain region of eastern Europe have unearthed more than 150 clay models of four-wheeled wagons, all of them, judging from a big loop handle at one end, designed for use as drinking mugs (see figure 3a). Carbon-14 analysis dates the formation of the Boleráz culture that produced the mugs to no later than 3600 BC [15]. This makes the mugs the world's earliest known representations of wheeled transport. But it leaves a key question unanswered: why did the Boleráz people enjoy drinking out of square mugs on wheels? Ceremonial, even iconic, the mugs must reflect an important aspect of Boleráz culture, though one that had no apparent precursor or successor, at least in the area of tableware. Several particulars point to the inspiration for the mugs being small, wheeled baskets used to carry ore in the trenches or tunnels of copper mines. First, the mugs have wheelsets rather than wheels that revolve independently at either end of a non-revolving axle. This is consistent with a mine environment which, unlike a farmer's field, facilitates digging and smoothing pathways that are straight and level. Second, pre-modern miners ordinarily pushed minecarts along pathways that were too constricted for draft animals, and the mugs do not show traces of yokes or other harnessing. Third, some mugs feature grooved side panels (see figure 3b) suggestive of the basketry shown in depictions of ancient mining [14].

The Copper Age began in the Balkan Mountains, south of the Carpathians, before the formation of the Boleráz culture. The Boleráz people had less of the metal than their southern precursors, indicating perhaps, poorer ores [14]. However, if poorer ores meant that a greater weight had to be transported from the mine to the smelter, then the invention of a wheeled basket could have been a technological breakthrough well worth celebrating with a round of drinks, at least until the wheel idea found copycats outside the mine. Notably, the earliest known transport wheels come from slightly later sites bordering the Carpathians (see figure 3c) where larger, more steerable carts equipped with a single wheelset were pulled by draft animals [14].

While there is no direct evidence for the use of rollers in the Carpathian region, the historical record contains multiple reports of roller-based transport being used by various civilizations that span vast distances of space and time from pre-colonial Fiji [17] to Assyria [8]. Furthermore, experimental archaeologists have successfully replicated the process of transporting megaliths using rollers in order to test their viability as an ancient method of transport [18][19][20][21]. Included among these studies is the 2001 investigation by Osenton, who showed that rollers could be fabricated using axe-based technology, and the resulting contraption enabled small groups of skilled workers to transport massive objects [21] (figure 4).

Computational analysis and design methodology

Here, we present the methods used for the analysis and generation of the optimized wheelset design. We developed several novel topology optimization methods for this purpose, the first of which is an original formulation for design-dependent variable loads in three dimensions. A variable void region, i.e. a region of empty space that is movable by the optimization algorithm, prevents the load from being applied inside the three-dimensional structure while also allowing the optimizer to place it at any location in the design domain. This enabled us to solve the wheelset design problem, which was not previously possible using any currently available topology optimization methods. Another novelty of the work is our method for achieving an axisymmetric geometry on a fixed, rectangular grid of finite elements. We use a 'turning' filter method which replaces the density filter matrix of standard topology optimization. Additionally, we model the orthotropic wood material in the same cylindrical coordinate system as the axisymmetric geometry filter. A wickerwork pattern can be seen on the side panel of the cart, suggesting basketry as the likely manufacturing technique [14]. (c) The Ljubljana Marshes Wheel, the oldest known transport wheel, discovered in Slovenia and estimated to be between 5100 and 5350 years old [3]. The square hole in the middle of the wheel indicates that this was a wheelset in which the wheels were fixed to the axle. (d) A wheeled dog figurine from pre-Columbian Mesoamerica created in the eighth century AD [16].

Physics of a wheelset

Before we can perform topology optimization, we need an objective function to minimize. For this, the expression for the force required to achieve the rolling of a wheelset in a plain sliding contact bearing is derived using statics. A free-body diagram of a cart with two wheelsets is shown in figure 4a, where P is the applied pushing force, m c is the mass of the cart excluding the wheelsets, m w is the mass of one wheelset, g is the acceleration due to gravity and F gx and F gy are the reaction forces with the ground on each wheelset. Summing the horizontal forces gives (4.1)

Summing the moments about the front wheelset's contact point with the ground leads to (4.2)

and summing moments about the rear wheelset's contact point yields

Now, we consider one wheelset in its plain bearing by itself. The diameter of the bearing is slightly larger than the axle so that they only come into contact at a single point. Initially, with no external load applied, the bearing's surface rests on the top of the axle. Then, as the pushing force is applied slowly, the contact point moves to the left along the axle's surface until the friction between the bearing and axle is overcome and rolling begins [22]. In the analysis that follows, we assume the wheelsets have constant angular velocity, and therefore, they are in rotational and translational equilibrium such that all external forces and moments sum to zero. The free body diagram for a wheelset at equilibrium is shown in figure 4b, where the overall force applied to the axle, P ‾, consists of an unknown horizontal and vertical component written in terms of the mass of the wheelset, m w , and the reaction forces with the ground, F gx and F gy . P ‾ can also be decomposed into the components normal and tangent to the axle surface. The tangent component is the force of friction, F f , which depends on the normal component N as F f = μN, where μ is the coefficient of friction. We sum the moments about the centre of the axle, giving (4.4) F gx = F f r axle r wheel .

Using Pythagoras' theorem to find an expression for N, we write the force of friction as (4.5)

1 and substitute it into equation (4.4). Solving for F gx leads to (4.6) F gx = F gym w g μ r axle r wheel 1 + μ 2 1 -r axle r wheel 2 . Now, with an expression for the horizontal reaction force, we insert this identity into equation (4.1), eliminating F gx and resulting in (4.7) P = F gy1 + F gy2 -2m w g μ r axle r wheel 1 + μ 2 1 -r axle r wheel 2 . Substituting equations (4.2) and (4.3) into equation (4.7) for the vertical reaction forces leads to our final expression for the force required to push the cart, (4.8) P = m c gμ r axle r wheel 1 + μ 2 1 -r axle r wheel 2 . royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 11: 240373 Downloaded from https://royalsocietypublishing.org/ on 29 January 2025

The equation is plotted in figure 5. Note that this equation holds for both the wheelset and the multi-body wheel-and-axle, provided the sliding surface contains a plain bearing, as was the case for early wheels.

Topology optimization with a variable load

We use a three-dimensional density-based topology optimization method with a variable load location based on our previous work [23]. The variable load location is used in this problem to model the contact force between the wheelset and the cart that rests upon it. In contrast, other works modelling contact between components of a structural assembly use techniques such as spring, bar or beam element connections between coupled design domains [24][25][26].

A rectangular design domain of dimensions L x × L y × L z is discretized into a regular grid of voxels, which each contain a continuously variable density of material that can range from nearly empty space to completely solid. These densities are controlled indirectly through a vector of density design variables, ρ, while the location of the applied force representing the weight of the cart is controlled directly by three variables x f , y f and z f representing its coordinates in space. All quantities used for the topology optimization are formulated to be continuously differentiable with respect to the design variables in order to allow for gradient-based optimization using the method of moving asymptotes (MMA) [27].

The finite element method is used to evaluate the structure defined by the design variables, treating each voxel as an eight-node hexahedral finite element. Using linear elastic physics, the finite element problem consists of solving the following linear system for the vector of nodal displacements, U:

where K is the global stiffness matrix and F is the vector of nodal forces. We use the standard solid isotropic material with penalization (SIMP) method, which interpolates the element's stiffness between solid and void while penalizing the stiffness-to-weight ratio of intermediate densities. The assembly of the global stiffness matrix is then given by (4.10)

where N e is the number of elements, k e 0 is the element stiffness matrix of a fully solid element, ρ ‾e is the physical density of the element (determined from ρ through filtering, described in §4.3), p is the penalization factor, and ρ min is the minimum allowable value for the densities (here chosen to be 1 × 10 -6 ) to avoid singular global stiffness matrices. The element stiffness matrix is given by where V e is the volume of the element, B e is the strain-displacement matrix and C e 0 is the constitutive matrix for the solid material.

A continuously differentiable variable load is achieved by applying forces to every element centroid, as illustrated in figure 6, and controlling the magnitudes of those forces through the design variables x f , y f and z f to change the effective location of a single applied load. The shape of the applied load is defined using a super-Gaussian projection function [23], which takes a distance function as input, and outputs a geometric shape of a given radius. In this case, we use the distance function for a single point to model a concentrated load, where the distance from each element centroid to the load coordinate is given by (4.12)

The distance function is passed to the super-Gaussian function, which projects it into a spherical shape of radius r defining the distribution of force magnitudes, (4. 13)

Here, f e is the force applied to the element centroid, A f is the maximum magnitude of the distributed force and S controls the sharpness of the transition of the force magnitude from A f to 0. The radius parameter r determines the length from the load coordinate to a point in the transition region with magnitude A f divided by the base parameter b. The maximum distribution magnitude, A f , is chosen to give a total load equivalent to a specified number, (4.14)

This is assumed to be a constant and is calculated only once before the optimization begins to make it easy for the user to directly specify how much load is applied.

With the spatial distribution of element centroid forces defined, the global vector of nodal forces in the finite element problem is assembled as follows:

(4.15)

where the force at the centroid is evenly distributed to the nodes by dividing it by the number of nodes in the element, n n e . The vector Θ assigns a direction to the nodal forces by applying them to the appropriate degrees of freedom, which in this case is straight downwards, Figure 6. A finite element mesh with forces applied to the element centroids, which are evenly distributed to the nodes. The magnitudes (and orientations [23]) of the element forces can be varied to smoothly control the effective location of an applied load in a differentiable way.

Axisymmetric geometry

To limit the design space to only axisymmetric structures, a two-stage density filter is used to constrain the possible geometries to only those that can be made using a turning manufacturing process, such as a lathe. Vatanabe et al. have previously implemented a turning manufacturing constraint using a mapping approach [28]; however, our turning filter is an extension and modification of the casting and milling restriction filter proposed by Guest & Zhu [29]. Filtering-based manufacturing constraints have also been used for overhang limits in additive manufacturing [30,31]. In the first stage of our method, the base density design variables are mapped to an axisymmetric configuration using an N e × N e filter matrix W, (4.17 N e w ei ) .

The weight factor w ei is non-zero only when element e is closer to the axis of symmetry than element i, and its magnitude is inversely proportional to the area of the circle centred on the axis of symmetry and intersecting the centroid of element e, (4. 19)

These circular regions are illustrated graphically in figure 7. The effect of applying the filter matrix is to project each unfiltered density ρ i to every element within the circular projection domain of that element. Each unfiltered density adds to the densities of the filtered variables, so the weight is increased closer to the axis of symmetry to compensate for fewer unfiltered variables contributing to the density of the filtered ones.

The second stage of the turning filter consists of a smoothed Heaviside step function.

(4.20)

where the parameter β affects the sharpness of the step function.

To avoid having the variable load penetrate the solid structure, we create an axisymmetric region of empty space that follows the location of the load. The void region begins just beyond the radius of the load application zone and extends radially outward, such that the edge of the loading zone is always Filtered element e is inside the unfiltered element i's region of influence (blue), and therefore, the filter matrix is assigned a weight at entry (e, i) that is inversely proportional to the area of the circle intersecting element e (red).

applied at the solid-void interface. We use the super-Gaussian projection function to create this void space, using the minimum distance function for a yz plane with a circular hole centred on the axis of symmetry, as shown in figure 8. Defining a vector a, which points from the centre of the void region to the point that is a distance 2r above the load coordinate,

and a vector h = ℎ x , ℎ y , ℎ z , which points from the centre of the void region to an element e's centroid,

the minimum distance from the centre of the void region geometry to that element's centroid is given by (4. 23)

The distance function is then projected to a density field of the same radius as the load application region before being subtracted from a uniform field of 100% density, (4.24)

S

, where the coefficient (1ρ min ) is used to maintain the minimum density value. Finally, using a smooth minimum function, the void region is combined with the filtered densities to give the physical density field that represents the actual design on which the finite element simulations are based, (4.25)

Here, Q is a negative number such that the smooth minimum function approaches the true minimum as the magnitude of Q increases.

Orthotropic material properties

We use orthotropic material properties in a cylindrical coordinate system about the axis of symmetry. Topology optimization using orthotropic materials is not new and has been used for applications such as continuous fibre composites [32]. In our application of a structure made from a single piece of wood, the material coordinates only need to be defined once before the optimization begins. Wood's three mutually perpendicular axes are the longitudinal, radial and tangential directions. The longitudinal direction runs parallel to the axis of symmetry, the radial direction is normal to the log's growth rings, and the tangential direction is parallel to the rings [33]. We use the approximate material properties of spruce trees obtained from data in [33] and shown in table 1. The remaining Poisson's ratios ν 23 , ν 31 and ν 21 are calculated as (4.26)

The constitutive matrix for a solid element in a local coordinate system is (4.27)

, where the subscript 1 is the longitudinal direction, 2 is tangential and 3 is radial. Now, for each element in the mesh, the constitutive matrix is rotated an angle θ so that the radial direction points away from the axis of symmetry, (4.28)

The rotation matrix R(θ) is

Numerical implementation

The optimization task is formulated as a multi-objective problem where we seek to minimize the force required to move the cart as well as the structural compliance, given a fixed amount of material to work with. The pushing force is given by equation (4.8), and the compliance, C, is the dot product of the load and displacement vectors, (4.30)

The amount of material is constrained using the volume fraction, (4.31)

Our optimization problem is then written as aC

where P is the pushing force evaluated using equation (4.8) with the axle radius calculated as the distance from the point of load application to the axis of symmetry,

and the coefficient a = 1 m -1 , which ensures dimensional consistency across terms within the objective function. Bounds on the design variables are set such that the densities can vary between nearly void and solid, and the applied load can move anywhere in the upper half of the domain.

Taking advantage of the symmetry of the problem, only one-half of the wheelset is modelled and optimized. The symmetry is modelled using roller boundary conditions by fixing the displacements in the x direction on the face of the domain at x = 0. Boundary conditions representing the ground that the wheel rests on are placed at 10% of the domain's width above the lower face and away from the side faces so that the wheel can be supported by more than a single point by sinking into the ground slightly. The ground extends 20% the width of the domain from the face at x = L x , and the x and z displacements are fixed on this area. A single node at the middle of the inner edge is also fixed in the y direction to fully constrain the problem by preventing rigid body translation. The variable load is initialized at the top face and near the x = 0 face, with a magnitude equivalent to one-quarter of the weight of a 100 kg cart. The domain and boundary condition set-up described here is illustrated in figure 10.

The program is implemented in C++ using the PETSc library [34,35], based on the code originally written by Aage et al. [36]. Using a mesh of 104 × 104 × 104 (1,124,864) eight-node hexahedral (cube) elements and the parameters shown in table 2, we ran the problem on the Illinois Campus Cluster using 32 nodes and 640 processes. The optimization search was considered to have converged when the maximum change in the density design variables was less than 0.001. The optimization completed in 9 h and 38 min after 776 iterations, with the final design shown in figure 2.

Failure analysis

Failure of orthotropic materials, like wood, can be predicted by performing a stress analysis and computing a failure index. Many different failure criteria for orthotropic materials exist, such as the Tsai-Hill [37], Hoffman [38] or Tsai-Wu [39] failure criteria. The Tsai-Hill index has been found to work adequately well for predicting wood failure [40,41]. We use the approximate strength properties Move Limits: ρ 0.05

Move Limits: x f , y f , z f 0.9375 cm of spruce wood [33], shown in table 1, to analyse the wheelset design. The Tsai-Hill failure index is given by (4.33)

where the superscript f denotes a failure stress material property. Violation of the inequality indicates that failure has occured. Linear elastic finite element analysis was performed in Ansys using the orthotropic material properties from tables 1 and 3 in a cylindrical coordinate system about the axis of symmetry. Two models were created: one as a single solid component, and a second where the axle is separated from the wheel to create a wheel-and-axle assembly. Only half of the wheelset was modelled to save on computational cost, and a portion of the bottom of wheel was cut off to apply a fixed boundary condition simulating a contact patch. A vertical downward load of 2000 N was applied to the face of the axle next to the wheel where the load application zone was applied in the topology-optimized design. Using equation (4.8) with the dimensions of the model (r axle = 0.037 m and r wheel = 0.17 m) gives a corresponding pushing force of 986.5 N, which was applied as a horizontal load on the same face of the axle. In the wheel-and-axle assembly model, a frictionless contact condition was applied to the bearing surface. Each model was meshed with approximately 60,000 quadratic tetrahedral elements, with increased refinement in the axle.

Results

Figure 11 shows the progression of the wheelset structure during the execution of the topology optimization algorithm. The design search begins with a solid prismatic roller. From this starting point, two grooves appear and grow progressively deeper. This is followed be the narrowing of the central portion of the roller to form an axle, and the expansion of the outer discs to form wheels. The resulting wheelset design maximizes mechanical advantage and structural stiffness.

Figure 12 shows the optimization convergence history for the result shown in figure 11. The first plot shows the convergence of the combined objective function, which is defined as the sum of the pushing force plus the structural compliance. The second plot contains the convergence history of the spatial coordinates of the point of load application. There is a noticeable dip in the objective function beginning around iteration 300, when the optimizer pushes the contact point outward towards the wheels, causing a reduction in structural compliance. Note also that both convergence plots contain regions where the objective and constraint function values oscillate. These regions reflect the inherent difficulty of numerically solving a mixed topology optimization problem that combines density variables with explicit variables representing the load coordinates. This combination negatively impacts the convergence of the algorithm. Here, we mitigate this effect by imposing move limits on the design variables.

The pushing force and structural stiffness objectives are in conflict. This can be shown by plotting a Pareto front [42,43] of the optimization problem described in equation (5.1). The objective function is modified by multiplying the pushing force (P) and the compliance (C) terms by a variable weighting Table 3. Strength properties.

Using a resolution of 80 × 80 × 80 elements and a volume fraction limit of 10%, we solved the optimization for several values of α ranging from 0 to 0.999. When α = 0, the optimization is purely a compliance minimization problem. As α approaches 1, the pushing force becomes more heavily weighted. In this case, the optimizer drives the pushing force to its minimum possible value of zero. The values of the pushing force P and the structural compliance C are recorded for each solution, and the Pareto front is plotted in figure 13 along with images of some of the optimal designs. It can be seen from the Pareto plot that if one objective decreases, the other increases, hence the trade-off. 0.5 0.475 0.45 0.425 0.4 0.375 0.35 0.325 0.3 0.275 0.25 0.225 0.2 0.175 0.15 0.125 0.1 0.075 0.05 0.025 0 Figure 14. Finite element simulation of the wheelset (top row) versus the multi-body wheel-and-axle (bottom row). The contours show the Tsai-Hill stress distribution within the volume of the structure. The results show that the multi-body wheel-and-axle system experiences reduced stress and is therefore less likely to sustain material damage and structural failure.

Figure 14 shows a comparison of the stress distributions in wheelset and multi-body mechanisms with identical geometries and loading conditions. In the monolithic wheelset, the peak Tsai-Hill failure index on the axle is 0.56, while on the wheel-and-axle assembly, it only reaches 0.43 because of the frictionless bearing preventing the transfer of a torsional stress between the wheel and axle. This result supports the hypothesis that a multi-body wheel-and-axle system would have been more robust to damage than a wheelset structure. Note that the optimization results shown in figures 11-13 did not account for the Tsai-Hill failure index and that failure analysis was performed as a post-processing step to analyse the optimized designs.

Significance and implications of our findings

These results indicate that our proposed sequence for the progression of the wheel and axle offers a plausible descent path along which the wheel could have evolved as its users sought more energyefficient designs. Our findings also demonstrate the critical role that environmental factors played in the creation of wheeled technology. The unique features of the mine environment accentuated the advantages of the wheelset over its predecessor while negating its most significant disadvantage: the inability to turn. Despite its shortcomings, the wheelset played a crucial role in the evolution of the wheel due to its conceptual proximity to the roller, which served as a technological bridge.

It should be noted that our investigation does not span the full breadth of the wheel's historical development, and the wheel did not stop evolving with the advent of the multi-body wheel-and-axle. The technology came full circle in 1869 with the invention of radial ball bearings [44], which, ironically, reintroduced bilateral rolling but solved the issue of spent rollers by enclosing the rollers (i.e. bearings) within a cage. Nor do we wish to imply that this was the only time in history that the wheel was independently discovered. The indigenous peoples of the Americas also had knowledge of wheeled locomotion as evidenced by the presence of wheeled figurines from the eighth century BC in what is now southern Mexico (see figure 3b) [16,45]. However, the Boleráz wheels represent the earliest known instance of wheeled transport, and they are a likely source of the wheel's early proliferation throughout the Old World.

This chapter in human history runs counter to the popular belief that technologies arise abruptly from the sudden epiphany of a lone inventor. Consequently, some scholars have embraced the other extreme, claiming that the wheel had no origin point and no inventor, but rather it developed gradually across a broad geographical area [5]. A similar theory has also been applied broadly to technology in general [46]. Our investigation modifies this theory and adds several new insights. The wheel evolved over a period spanning many centuries, but this evolution was punctuated by discrete innovations. We also find that, much like biological organisms, the wheel evolved from an earlier technology from which it inherited certain advantageous traits, and this evolution was fundamentally linked to the local environment. Lastly, our study is an example of how we can leverage the techniques of design science and computational mechanics to uncover new knowledge about episodes from our distant past for which there is no written record, and our only clues are the artefacts left behind by the designers who came before us.

Ethics. This work did not require ethical approval from a human subject or animal welfare committee.