Using geometric algebra to represent curvature in shell theory with applications to Starling resistors

We present a novel application of rotors in geometric algebra to represent the change of curvature tensor that is used in shell theory as part of the constitutive law. We introduce a new decomposition of the change of curvature tensor, which has explicit terms for changes of curvature due to initial curvature combined with strain, and changes in rotation over the surface. We use this decomposition to perform a scaling analysis of the relative importance of bending and stretching in flexible tubes undergoing self-excited oscillations. These oscillations have relevance to the lung, in which it is believed that they are responsible for wheezing. The new analysis is necessitated by the fact that the working fluid is air, compared to water in most previous work. We use stereographic imaging to empirically measure the relative importance of bending and stretching energy in observed self-excited oscillations. This enables us to validate our scaling analysis. We show that bending energy is dominated by stretching energy, and the scaling analysis makes clear that this will remain true for tubes in the airways of the lung.


Introduction
Self-excited oscillations of flexible tubes driven by fluid flow have been a subject of interest for some time, and there is a considerable literature on the subject, which is reviewed in [1][2][3][4][5][6]. Experimental rigs designed to study this phenomenon are often called Starling resistors. We are interested in this phenomenon because of its possible relevance to wheezing in the lung [7], which is one of the most commonly heard lung sounds used for diagnosis [8,9]. Previous work on Starling resistors has largely used water as the working fluid. In the lung, the working fluid is air. This means that the density ratio between the working fluid and the tube material in the lung is significantly different from almost all of previously completed work on Starling resistors. Previous modelling work 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited. usually neglects wall inertia, using instead a 'tube law' [10], but there is strong evidence that wall inertia is significant in the lung from [8,11], where it was shown that when the density of the fluid breathed in is changed, the frequency of the wheezes is not affected significantly. It is clear, therefore, that the change in density ratio results in a qualitatively different mechanism. For this reason, we have been conducting our own experiments, and creating models to understand the onset of oscillations.
The flexible tube itself is generally modelled as an elastic shell. Traditional shell theories [12][13][14][15][16][17] are well developed but difficult to implement. We have found that linearized shell theories do not provide good predictions of the frequencies of oscillation, and we believe that this is due in part to the fact that we have observed that oscillations start from a collapsed or partially collapsed state. To use current geometrically nonlinear shell theories would require a numerical simulation of a very complex fluid structure interaction problem, which would be of similar value to experimental results (though arguably harder to implement), and would provide the same problem of being difficult to physically interpret due to the complex nature of the oscillations observed. Instead, we would like to gain a physical understanding of the important mechanisms behind the oscillations, and this is difficult with shell theories based in differential geometry [14,15], in particular, due to the lack of physical interpretations of the change of curvature tensor in general situations. We recently introduced geometric algebra to shell theory [18], which allowed us to express the fundamental laws in a component-free form and clarify the role of angular velocity and moments through the use of bivector representation. For an introduction to the basics of geometric algebra see [19]. One of the most powerful aspects of geometric algebra lies in the use of rotors to represent rotations. In [20], these have been used to simplify Simo and Vu Quoc's numerical algorithm [21] for modelling the nonlinear behaviour of rods. In projective and conformal geometry [19, §10] rotors have allowed geometric primitives to be represented in a more simple and lucid manner, and in relativity [22] and relativistic analogies [23] rotors can simplify transformation between frames of reference. In this paper, we make use of rotors to better understand the change of curvature of a shell, which is of prime importance to the constitutive law of the shell, but whose representation has long caused controversy. In [12] at least 10 different linearized shell theories are presented, and the differences are primarily caused by disagreements over how to represent changes of curvature. The author in [14] has provided a tensor definition of the change of curvature that has become accepted; however, the utility of this expression is limited by its complexity. We have been able to simplify the representation of this tensor using rotors, allowing a more lucid and physical interpretation of changes of curvature. We take advantage of this to allow us to understand the importance of the change of curvature in the context of our Starling resistor experiments.
To compare results from the shell theory to our experimental results, we need to be able to calculate the kinematic parameters associated with the deformation of the flexible tube. To enable this, we use stereoscopic imaging, which, to our knowledge, is the first use in the study of Starling resistors. We take high-speed video of the tube at the onset of oscillation, and are able to track the motion of the surface, and consequently compare the predictions of shell theory with empirical calculation.

Understanding changes in curvature
There is an energy associated with any deformation of a shell, and Koiter [17] proposes the following form for this energy: Here, U is the internal energy per unit mass of the shell, defined on the reference configuration, ρ 0 is the time-independent area density of the shell on the reference configuration, E is Young's modulus, ν is Poisson's ratio, h is the shell thickness (which in Koiter's theory is assumed constant), E is the twodimensional Green-Lagrange strain tensor defined on the reference configuration, H is the change of curvature tensor and tr is the trace operator. From (2.1), we can derive the governing equations of the shell (for more details see [18]). The first term on the right-hand side of (2.1) represents the stretching energy, and the second term represents the bending energy. In general, a shell is a body in which the thickness is smaller than the other relevant defining length scales. The use of (2.1) for the energy of deformation implies that we additionally assume that the midsurface of the body remains the mid-surface under deformation, a material line that is normal to the mid-surface remains normal to it under deformation, the shell thickness remains constant with time, the first and second moments of density relative to the mid-surface are zero, and strains within the shell are small and so is the normal stress (see [18] for further discussion). Following the notation of [18], we take B and S to be the reference and spatial configurations of the shell, and X ∈ B and x ∈ S to be locations on these configurations. By φ t we denote the motion of the shell, meaning that, at time t, the point X ∈ B is at the position φ t (X) ∈ S; G and g are the identity functions on the reference and spatial configurations. Y and y are vectors within the tangent spaces of B and S, respectively; {X i }, i = 1, 2 is a coordinate system on the reference configuration B, which we can then use to define the frame on the reference configuration {E i = ∂X/∂X i }. By {E i } we denote the reciprocal frame that satisfies E i · E j = δ i j , and {x i }, {e i } and {e i } are the similarly defined coordinate system, frame and reciprocal frame on the spatial configuration, respectively. The shell undergoing deformation is embedded within a flat three-dimensional Euclidean space E 3 .
If A and B are general multivectors, then AB is the geometric product between them, A · B is the inner product, A ∧ B is the outer product and A × B is the commutator product, defined by A × B = 1 2 (AB − BA) (see [19, §4.1.3]). We also take × (compared to ×) to be the cross product between two vectors. If I is the pseudoscalar of a three-dimensional space, and a and b are vectors, then a×b = −Ia ∧ b.
We are particularly interested in the change of curvature that is encoded in H. To understand this, we must understand the curvature tensors on the reference and spatial configurations B and b. If E 3 and e 3 are the normal vectors to the reference and spatial configurations, respectively, then B(Y) and b(y) are given by where ∂ is the intrinsic vector derivative to any surface. The relationship between ∂ and the vector derivative of E 3 is explained in [18].
where F is the deformation gradient, defined by F(Y) = Y · ∂φ t (X). We see that F maps from the tangent space of B to the tangent space of S, providing information about the local deformation of the surface, andF is the adjoint of F, i.e. F(Y) · y = Y ·F(y). The strain tensor E, used in the constitutive law (2.1), is given by This much is well known, though in other treatments coordinate-dependent definitions of H are used (e.g. in [14]). To make further progress, we will now use rotors to better understand what will produce changes in H.
To begin, we note that we can perform a polar decomposition on F such that F(Y) = RU(Y), whereR = R −1 , det R = 1 andŪ = U. We see that R encodes rotation, and U encodes stretching. We can choose to consider a frame {E i } on the reference configuration that is locally orthonormal. 1 In this case, To find e 3 , we need two unit vectors in the tangent space of the spatial configuration that are oriented in the same way as the pair F(E 1 ), F(E 2 ), and are orthonormal. This pair of unit vectors is given by R(E 1 ), R(E 2 ) and so e 3 is given by ). The fact that the function R is a rotation means that it has an associated rotor R such that R(Y) = RYR, where R is an even multivector that satisfiesRR = RR = 1, andR is the reverse of R. This allows us to write e 3 as where we have used the fact that any rotor will commute with I 3 . Thus we have shown that the rotation associated with the deformation is also the rotation between the normal vectors E 3 and e 3 , which makes intuitive sense. We have also extended the range and domain of R to E 3 , while the range and domain of U is still constrained to the tangent space of the reference configuration, and the range and domain of F are constrained to the tangent spaces of the spatial and reference configuration, respectively. Two results that we will find useful are We see that (2.6a) follows fromRR = 1; (2.6b) has implicit assumptions that require explanation. The expression on the left of (2.6b) tells us how e 3 varies over the spatial configuration in the direction defined by F(Y), which lies in the tangent space of the spatial configuration. On the right of (2.6b), e 3 = e 3 (x) has been mapped to a vector field on the reference configuration such that e 3 (X) = e 3 (φ t (X))∀X ∈ B. This allows the expression on the right to tell us how e 3 varies over the reference configuration in the direction defined by Y, which is tangent to the reference configuration. The equality of these expressions is a standard result when mapping derivatives between manifolds, which can be proved by considering derivatives with respect to convected coordinates that satisfy . From this point we will assume that {x i } are convected coordinates.
Using these results, we can writeFbF(Y) asFbF(Y) = −F(Y · ∂e 3 ), and the argument ofF can be expressed as where × is the commutator product. Hence, we can expressFbF(Y) as and finally we obtain an expression for H, This shows that there are two contributions to H. Firstly, if the reference configuration is at all curved (i.e. B(Y) is non-zero), then the strain of the shell, encoded in U − G, will result in a change of curvature. The second contribution is due to variation of the rotor R over the shell. These two kinds of change of curvature are illustrated well by an inflating sphere and deformation of a flat plate. As a sphere is inflated to become a larger sphere, the normal vector is unchanged, i.e. e 3 = E 3 , hence R = 1 everywhere. This means that the second term in our expression for H will be zero. However, the surface of the sphere will stretch, meaning that U − G will be non-zero. In addition, B will be non-zero for a sphere, which tells us that the first term in our expression for H will be non-zero. By contrast, for a flat plate B(Y) will be zero, meaning that R must vary over the plate in order for there to be any change of curvature. A variant on this expression for H can be obtained if we express the rotor R as R = exp(−A/2) = exp(−Âθ/2), where A is a bivector aligned with the plane of rotation whose magnitude is equal to the angle of rotation θ andÂ is a unit bivector (Â 2 = −1). Note that the direction of rotation is defined by the sign of θ and the orientation ofÂ together. We can set the convention that θ ≥ 0, in which case θ andÂ are uniquely defined if A is known. Given this definition, we can express Y · ∂R as, Using this, we can express H as The term thatF operates on is the commutator product of a vector and bivector, so we can replace the commutator product with a dot product, Taking the inner product of a bivector with RE 3R = e 3 means that only vectors tangential to the spatial configuration are retained, which then means thatF can operate and return vectors tangential to the reference configuration. Hence, our description confirms that the range and domain of H are both the tangent space of the reference configuration. The explicit expression for the two possible contributions to change of curvature, shown in (2.12), gives a new decomposition of the change of curvature tensor which will be of use when we try to understand the importance of bending in the Starling resistor.  (1), and then through a rotameter (2) used to monitor flowrate. The noise that the rotameter introduces into the flow, and any other noise, is isolated from the flexible tube by the upstream settling chamber (3). Air flows into the upstream clean flow tube (5) section via a shaped inlet (4) that reduces separation. A contraction (6) leads to the flexible tube (7), before an expansion (6 ) leads to the downstream clean flow tube (5 ) that exits into the downstream settling chamber (3 ). Suction is provided by a fan (8). The downstream settling chamber (3 ) isolates the flexible tube from the noise from this fan. Experiments were performed in the Acoustics Laboratory in the Department of Engineering at the University of Cambridge.

Experimental set-up
In the experiments relevant to this paper, the suction at (8) is gradually increased until the flexible tube just starts oscillating. With the tube oscillating in this way, high-speed video is recorded from 2 FASTCAM-ultima APX cameras (produced by Photron, https://photron.com) with a frame rate of 12 500 fps and a resolution of 512 × 256 pixels (greyscale). Our experiments require us to focus on a small flexible tube at reasonably close range. The Photron camera has an adaptor for Nikon lenses. We use a 50 mm lens combined with a 7 mm extension tube to allow us to focus on the tube and have it fill most of the frame. An aperture of f/2.8 is used.
The flexible tubes used are made out of rubber latex for which E = 1 MPa and ν = 0.5. The tube diameter is 6 mm, the wall thickness is 0.3 mm and the unstrained length is 19 mm. The tubes are held in an axially strained state, so the length of the tubes during the experiment is 25 mm.

Image processing
The high-speed cameras record at 12 500 fps, and are triggered together, so that every t = 80 µs, two images of the flexible tube are taken. A schematic of the two cameras and the flexible tube is shown in figure 2. Dots are drawn on the flexible tube (shown in white in figure 2), which indicate a set of material points we aim to track over time in three dimensions.
It is possible to find the characteristics 2 of two cameras such that if a point appears in simultaneous images from both cameras, the point's position in three-dimensional space can be triangulated [24,25]. Calibration involves taking at least 3, and in general between 10 and 20, simultaneous images of a chequerboard pattern in various orientations. From this the position of the two cameras relative to each other, and their internal parameters, can be calculated. In this method, cameras are modelled as pinhole cameras, meaning that the focal length, pixel size and skew are the important internal parameters. In addition, it is possible to account for radial distortion of the image by the camera lens, and tangential distortion, which occurs when the image sensor is not perfectly perpendicular to the line of sight of the camera. This calibration is performed using the computer vision system toolbox of Matlab [26]. The images produced by these pinhole camera models is what is illustrated in figure 2.
To find the three-dimensional tracks of the material points, we must first find the locations of the dots within each image. We refer to these coordinate pairs as points. The dots on the tube surface are drawn in white, are spaced by approximately 1 mm and have a diameter of approximately 0.7 mm. We take the material points to be the centres of the dots, and we find them by first taking a two-dimensional convolution of the image with a 'mexican hat' function of the form 1 πσ 4 where σ is the expected radius of the dot in the image. The convolution effectively smooths the image, removing any artefacts from drawing the dots, leaving peaks in the centroid of each dot. These peaks are then used as the locations of each point. Hence at each time instance, we have two collections of points, representing the material points as seen from each camera. If there are n points on the tube, and m frames in our video, then in total we will have found the locations of 2nm points. To make use of this data, we need to identify each unique material point in each camera and over time.
To associate points across time for a single camera's set of images, the points at t and t + t are compared, and if two points are within a certain distance of each other, then it is assumed that these represent the same point. This works because the frame rate (12 500 fps) is much larger than the frequency of the observed vibrations (approx. 500 Hz), so motions between frames are small.
A pair of corresponding points in the two camera images are illustrated in figure 2, but finding these pairings at each instant in time is more complex. First, we consider the line drawn from the focal point of camera 2 to dot i in the image, which we will call a ray. Anything on this ray in three-dimensional space will appear at the same highlighted location in camera 2. However, from camera 1, the ray will appear as a line. Therefore, if a point's location is known in one camera image, then it must lie on a specific line in the other image. This line is known as an epipolar line [25]. Hence, for a pair of points, one in each camera image, to correspond to the same material point, they must each lie on the epipolar line of the other. However, because of the specific arrangement of the cameras and dots, this does not usually provide a unique set of pairs. The relative positioning of the cameras means that the epipolar lines are all approximately horizontal, and the dots drawn on the tube are arranged in horizontal rows, so multiple dots can be very close to a given epipolar line. To overcome this, we specify 10 corresponding point pairs between the two images (20 points in total), and these 10 pairs are then used to find the best fitting projective transformation from camera 1 to camera 2. Applying this transformation to the image from camera 1 places each point close to the corresponding point in the image from camera 2, allowing all the remaining point pairs to be found. This result is then checked for consistency with the epipolar line condition.
In addition to tracing material points over time as the self-excited oscillations occur, we must also find the locations of the material points on the tube when it is unstrained. This is necessary for the calculation of the kinematic variables used in the expression for energy of deformation (2.1). To achieve this, a single image is taken from each camera when the tube is held in its unstrained state. Pairing of points must then also be completed between the two images of the tube in its unstrained state and images from the high-speed video of self-excited oscillations. This pairing is done using the methods described in the previous paragraph.
Once point pairs are known over time, the camera calibration can be used to find three-dimensional point traces over time. The spatial resolution of this trace is limited by the size of the pixels in the highspeed video. This results in point traces with distinct jumps in position. These jumps are by no more than 0.1 mm in three-dimensional space, compared to variations in position of the order of 2 mm over the course of the self-excited oscillations. For this reason, we smooth the three-dimensional point traces by fitting functions of the form (3.2) to the three position components, where A i and ω i are chosen to fit the empirical data. These fits work well because the videos are of quasi-steady behaviour at the onset of oscillation, and the observed motions are close to sinusoidal. Eight terms have been found to be sufficient to match the experimental data. This smoothing is necessary in order for derivatives of the point traces to give meaningful results.

Kinematic calculations
In figure  The surface fit and its derivatives are used to calculate all the kinematic properties of the undeformed and deformed surfaces.
In figure 3, we show a typical image of the principal strains of the flexible tube, i.e. the eigenvectors and values of E. An eigenvalue of 1 corresponds to no strain, a value less than 1 corresponds to compression and a value greater than 1 corresponds to tension. The eigenvectors give the direction in which the strain is occurring. We can see that the first principal strain is approximately aligned with the longitudinal direction and is tensile. This is due to the dominant pre-strain of the elastic tube. By contrast, the second principal strain, which is primarily in the azimuthal direction, is a mixture of compression and tension, and is generally closer to 1. In figure 4, we show the principal curvatures of the deformed surface, i.e. the eigenvectors and values of b. For a cylinder, which the tube is in its undeformed state, the principal curvatures would be 0 in the longitudinal direction and 1/a = 0.33 mm −1 in the azimuthal direction. We see that, in the deformed tube, the curvatures are still closely aligned to the longitudinal and azimuthal  directions, and can see that the squashing of the tube results in a slightly negative longitudinal curvature and a slight reduction in the azimuthal curvature towards the centre of the tube. But these effects are fairly small.

Calculation of bending and stretching energies
We now aim to gain more understanding of how the tube deforms. To do this, we will use a mixture of empirical and analytical techniques. More specifically, we can use the high-speed video reconstructions combined with the mathematical framework for shells already developed in §2.

Scaling analysis
We start by estimating the bending and stretching energies analytically, which requires us to estimate the values of tr(E 2 ), tr(E) 2 , tr(H 2 ) and tr(H) 2 .
To understand bending in the tube, we consider the representation of H derived in §2, and given in (2.12). We can write the action of U(Y) as U(Y) = λ 1 (Y ·Ŵ 1 )Ŵ 1 + λ 2 (Y ·Ŵ 2 )Ŵ 2 , whereŴ i are the unit eigenvectors of U. Using the approximation λ 1 ≈ λ and λ 2 ≈ 1, this becomes U(Y) ≈ λ(Y ·Ŵ 1 )Ŵ 1 + (Y ·Ŵ 2 )Ŵ 2 . We can write B(Y) as B(Y) = C(Y ·Ê 2 )Ê 2 , whereÊ i are the unit eigenvectors of B, and C is the principal curvature of the undeformed tube in the azimuthal direction. Combining these we have We know from figure 3 thatŴ 1 andŴ 2 are approximately aligned with the longitudinal and circumferential directions, so we can writeÊ 2 ·Ŵ 1 ≈ 0 andÊ 2 ·Ŵ 2 ≈ 1. Using this, we obtain This is saying that because the directions of principal strain and principal curvature are approximately perpendicular, the influence of strain on the change of curvature is removed.
We are now in a position to consider the rotor R, because it is changes in this multivector over the surface of the shell that are responsible for the change of curvature. By R we represent the rotation, and as is shown in §2 it is characterized by the bivector A = θÂ, whose magnitude gives the rotation angle in radians, and whose plane gives the plane of rotation. In figure 5, we visualize the angle and axis of rotation encoded in the rotation tensor R, corresponding to a typical deformation. The axes of rotation shown in figure 5 are primarily tangential to the surface, so the bivector A will be dominated by the components e 1 ∧ e 3 and e 2 ∧ e 3 , with little rotation in the e 1 ∧ e 2 plane, i.e. about the normal vector e 3 . Hence, we can write A as We have used the reciprocal frame {e i } instead of {e i } because it will allow us to use the propertyF(e i ) = E i . 4 We can extend the frames {e i } and {e i } to span E 3 by using the normal vector e 3 . Because e 3 is a unit vector and perpendicular to all of e i , e i , we can also write e 3 = e 3 , and we have the frame {e a }, a = 1, 2, 3 and {e a }.
Using this, we define the Christoffel coefficients γ a ib = e a · ∂e b /∂x i , i = 1, 2; a, b = 1, 2, 3. These also satisfy e b · ∂e a /∂x i = −γ a ib . Substituting A into the second part of the change of curvature tensor given in (2.12), using the fact that e 3 is normal to e 1 and e 2 , and γ 3 i3 = 0, we obtain Therefore, given that the first part of H in (2.12) is zero, E i · H(E j ) = H ij is given by (4.5) 3 The strained length of the tube is l = 25 mm and the unstrained length is l 0 = 19 mm, giving λ = l/l 0 ≈ 1.3. 4 The coordinate system is convected, so e i = F(E i ), and as always  We know that H is symmetric, so from this we see that our earlier assumption on the form of A must be joined by the condition ∂ i θ j = ∂ j θ i to produce consistent results. The frame {E i } can be chosen by us to be orthonormal. More specifically, we can align E 1 with the longitudinal direction and E 2 with the azimuthal direction on the unstrained cylindrical tube. From figure 3, we can see that under a typical deformation these basis vectors remain close to the axial and azimuthal directions. In figure 6, we give a schematic illustration of how E i maps to e i . In figure 6, we have also labelled the values of θ i where they are obvious. The regions where |θ 2 | ≈ π/4 can be seen empirically in figure 6. There are two unknown values of θ 1 shown as question marks. We can estimate the largest value of these rotations by assuming a straight line from the clamped tube end and the centre of the tube when the tube collapses completely at the centre. In this case, θ 1 ∼ arctan(a/(l/2))λ, where a is the tube radius and l is the tube length in its deformed state. The multiplication by λ is necessary because θ 1 is the e 1 ∧ e 3 component and e 1 is shortened by a factor of λ compared to the unit vector E 1 . Up to angles of 30 • , tan θ is within 10% of θ , so we will take θ 1 ∼ λa/(l/2) = a/(l 0 /2) at the point in question, where l 0 is the unstrained length of the tube. This, and the values of θ i shown in figure 6, allow us to make the following estimates: Because of replacement of arctan with the identity function, our estimate for ∂ 1 θ 1 will be an overestimate when the tube is very short. We can also estimate the values of the coefficients γ i jk using figure 6,  From this we see that the changes in the basis vectors are primarily in the e 3 direction, meaning that they do not contribute to H ij . If we take l 0 to be much larger than a, then the dominant term in H ij will be 1/a, but even if l 0 and a are of a similar order of magnitude, all of the H ij terms will be of the order 1/a. Hence, we expect tr(H 2 ) and tr(H) 2 will scale as (1/a) 2 = (1/3 mm) 2 = 0.1 mm −2 .

Direct calculation from data
We can use the high-speed video data to calculate tr(E 2 ), tr(E) 2 , tr(H 2 ) and tr(H) 2 from (2.1). Typical plots are shown in figure 7, from which we see that in the units chosen these have similar orders of magnitude. We also see that the scalings obtained in the previous section agree with these plots well, providing support for the assumptions made. We can also calculate the bending and stretching energy, and this is shown in figure 8. This agrees very well with the scaling values of the previous section, again supporting our conclusions.

Conclusion
We have developed a new method of representing the change of curvature tensor using rotors (see (2.12)), increasing our understanding of bending in shells. We have used this representation to explain results from stereographic imaging of Starling resistors that demonstrate that the bending energy in these deformations is approximately 2 orders of magnitude lower than the stretching energy. We have been able to show that this relies on the fact that the strain energy is dominated by the effects of pre-strain, the axis aligned with the largest strain remains close to perpendicular to the axis aligned with the largest curvature, rotations are mostly about axes tangential to the shell and changes in the rotation scale with the change of rotation about the longitudinal axis in the azimuthal direction. Further to this, our scaling analysis remains valid even when the tube length gets close to the tube diameter. This is of significance to our work in understanding wheezing, because the length to diameter ratio of tubes in the lung typically varies from 1 to 6. Hence we have provided a scaling analysis, confirmed by experiment, that allows us to say that bending energy is dominated by stretching energy during self-excited oscillations in the airways of the lung. This should allow the use of membrane theory to model the tube, which reduces the order of the equations of motion from 4 to 2.
Ethics. This paper adheres to the Royal Society publishing and research ethics policy. Data accessibility. Data for this paper can be accessed at University of Cambridge Repository: https://doi.org/10.17863/ CAM.10363 [28].
Authors' contributions. The work was completed as part of the PhD project of A.L.G, with A.A. supervising and J.L.
advising, in particular, on image processing and the uses of geometric algebra.