Determining rigid body motion from accelerometer data through the square-root of a negative semi-definite tensor, with applications in mild traumatic brain injury
DDetermining rigid body motion from accelerometer datathrough the square-root of a negative semi-definite tensor,with applications in mild traumatic brain injury
Yang Wan , Alice Lux Fawzi , and Haneesh Kesari Brown University School of Engineering, 184 Hope St., Providence, RI, USA * Corresponding author, [email protected] 27, 2021
Abstract
Mild Traumatic Brain Injuries (mTBI) are caused by violent head motions or impacts. Most mTBIprevention strategies explicitly or implicitly rely on a “brain injury criterion”. A brain injury criteriontakes some descriptor of the head’s motion as input and yields a prediction for that motion’s potential forcausing mTBI as the output. The inputs are descriptors of the head’s motion that are usually synthesizedfrom accelerometer and gyroscope data. In the context of brain injury criterion the head is modeledas a rigid body. We present an algorithm for determining the complete motion of the head using datafrom only four head mounted tri-axial accelerometers. In contrast to inertial measurement unit basedalgorithms for determining rigid body motion the presented algorithm does not depend on data fromgyroscopes; which consume much more power than accelerometers. Several algorithms that also makeuse of data from only accelerometers already exist. However, those algorithms, except for the recentlypresented AO-algorithm [Rahaman MM, Fang W, Fawzi AL, Wan Y, Kesari H (2020):
J Mech PhysSolids
Keywords: mTBI, Tensor square root, Rigid body motion, Accelerometers, Inertial navigation, Con-tinuum mechanics
Mild Traumatic Brian Injury (mTBI) is the most common injury among military personnel and it isestimated that as many as 600 per 100,000 people experience mTBIs each year across the world [1, 2].Mild Traumatic Brain Injuries are caused by violent head motions, that may occur from intense bluntimpacts to the head in contact sports, motor vehicle accidents, falls following blasts, etc. In mTBI, themotion of the head causes the soft tissue of the brain to deform. The magnitude and time rate of braindeformation can cause brain cells to die [3, 4, 5, 6, 7].There have been many strategies aimed at preventing mTBI. In sports, new rules aim to modify playerbehavior in order to decrease or eliminate exposure to blunt impacts [8]. Helmets and neck collars areexamples of equipment that can alter the motion experienced by the head. A jugular vein compressioncollar aims to change the stiffness of the brain making it less susceptible to injury [9].Most mTBI prevention strategies explicitly or implicitly rely on a “brain injury criterion” for theireffective synthesis, implementation, and evaluation. A brain injury criterion takes some descriptor of thehead’s motion as input and yields a prediction for that motion’s potential for causing mTBI as the output.When we refer to any aspect of the head’s motion, we, in fact, are referring to that aspect as it pertainsto the skull’s motion; since it is the skull’s motion that is, at least currently, observable and quantifiablein the field, either using video recording equipment or inertial sensor systems. The Young’s modulus1 a r X i v : . [ phy s i c s . m e d - ph ] J a n f bone from human skulls generally lies in the 2–13 GPa range [10, 11, 12]. In comparision, braintissue is extremely compliant. Recent indentation tests on brain slices that were kept hydrated showthat the Young’s modulus of brain tissue lies in the 1–2 kPa range (white matter . ˘ . kPa, and graymatter . ˘ . kPa) [13]. Due to this large disparity between the skull’s and the brain’s stiffnesses,in biomechanical investigations of mTBI the skull is usually modeled as a rigid body [14, 15]. Thus,inputs to the brain injury criteria are rigid body motion descriptors, such as angular velocity time series,translational acceleration times series, etc., or a combination of such time series.Rigid body motion can be thought of as a composition of translatory and rotatory motions. In initialbrain injury criteria the focus was on the head’s translatory motion. Two of the first published injurycriteria are the Gadd Severity Index (SI) and the Head Injury Criterion (HIC) [16, 17]. Both SI and HICignore the head’s rotations and take the head’s translational acceleration as their input. Later, however,it was realized that in the context of mTBI the head’s rotations play an even more important role incausing injury than its translations. The first brain injury criterion to take the rotational aspect of thehead’s motion into consideration was GAMBIT [18]. The input to GAMBIT is the tuple of center-of-mass-acceleration and angular-acceleration time series. Following the development of GAMBIT, brain injurycriteria that use descriptors that only depend on the head’s rotational motion as inputs have also beenput forward. One such criteria is the Brain Injury Criteria (BrIC) [19]. Aiming to compliment HIC, BrIConly uses the head’s angular velocity time series as input. We also note that there is currently significantactivity in applying finite element modeling using 2D/3D anatomically consistent discrete geometry headmodels to evaluate or develop new brain injury criteria [20, 21, 22].Irrespective of which existing, or yet to be developed, brain injury criterion will be used in the future,its successful application will hinge on the availability of a robust algorithm for constructing the motiondescriptor that the criterion takes as input from easily measurable data. Currently, different algorithmsare used to obtain the descriptors taken by the injury criteria as inputs. The inputs to GAMBIT can beobtained from the measurements of one tri-axial accelerometer and one tri-axial gyroscope mounted in amouthguard [23] if the center-of-mass-acceleration and angular-acceleration are obtained by processingthe data using the algorithm in [24]. In another example the input to BrIC (i.e., angular velocity) isprepared by numerically integrating the angular acceleration, which is determined by applying the 6DOFalgorithm [25] to the data from 12 single-axis accelerometers mounted in a helmet [26]. Interestinglythe inputs to most of the currently employed brain injury criteria can be prepared from the knowledgeof a few key rigid body motion descriptors. To make this idea more concrete, consider the followingequation, which is often used to describe rigid body motion, x p τ q “ Q p τ q X ` c p τ q . (1.1)In (1.1) τ is a real number that denotes a non-dimensional time instant; X is a column matrix ofreal numbers that denotes the initial position vector of a rigid body material particle; x p τ q is the columnmatrix of real numbers that denotes that material particle’s position vector at the time instance τ ; Q p τ q is a time dependent square matrix of real numbers with positive determinant whose transpose equalsits inverse; and c p τ q is a time dependent column matrix of real numbers. The matrix Q p τ q quantifiesthe rotation or orientation of the rigid body at the time instance τ , while c p τ q quantifies the rigid body’stranslation at that time instance. The inputs to most current brain injury criteria can be computed froma knowledge of the maps Q and c and their first and second-order time derivatives, i.e., Q , c , Q , c .In this manuscript we present an algorithm for determining these maps and their derivatives using datafrom only four tri-axial accelerometers.The presented algorithm has some similarities to the one recently presented by Rahaman et al. [27] ,which is referred to as the AO (accelerometer-only) algorithm. For reasons that will become clear shortly,we refer to the algorithm that we present in this manuscript as the ? AO-algorithm. The AO-algorithmalso presents a framework for completely determining the rigid body’s motion, i.e., for constructingthe maps Q and c and their time derivatives, using data only from four tri-axial accelerometers. The ? AO-algorithm has all the advantages of the AO-algorithm.There are existing algorithms for completely determining a rigid body’s motion using sensor data.However, these algorithms take data from sensor systems called inertial measurement units (IMUs).These units contain one or more gyroscopes. One of the primary advantages of the AO and ? AO al- A graphical user interface for applying the AO algorithm to different types of data sets and visualizing its results is freely available [28] § et al. [32], Genin et al. [33] and Naunheim et al. [34]). Thesealgorithms, however, give much more limited information than is given by the AO and the ? AO algo-rithms. For example, all these algorithms give the rigid body’s acceleration field in terms of the bodyframe, which is a set of vectors that are attached to the rigid body, and hence move with it. Thesealgorithms do not provide any information of how the body frame is oriented in space. However, thatinformation is critical for constructing inputs for the upcoming finite element based brain injury criteria.The AO and the ? AO algorithms provide complete information of how the body frame is oriented inspace. See § ? AO algorithms overother algorithms that also make use of only accelerometer data.Despite its many advantages we note that the AO-algorithm has one critical limitation. It is quitesensitive to bias type errors in the accelerometer data. Bias type errors are distinct from random errorsin that they do not arise as a consequence of stochastic processes. For accelerometers, bias type errorscan arise as a consequence of inaccurately defining sensor position and orientation (see Fig. 1). As weexplain below, the advantage of the ? AO-algorithm over the AO-algorithm is that it is far less sensitiveto bias type errors than the AO-algorithm.One of the critical steps in the ? AO and the AO algorithms is the determination of the map τ ÞÑ W p τ q .Here W p τ q is time dependent skew-symmetric matrix of real numbers that is related to the rigid body’sangular velocity. In the AO-algorithm W is determined by numerically integrating the equation ([27,§3.12]) W p τ q “ skew part of P p τ q . (1.2)Here P p τ q is a square matrix of real numbers that is to be computed from the accelerometers’ data,relative locations, and orientations. Due to numerical integration any bias type errors in P will give riseto errors in W that grow with time. In the ? AO-algorithm we alternatively determine W by taking thesquare-root of the equation W p τ q W p τ q “ symmetric part of P p τ q . (1.3)We derive (1.3) in B. Due to the elimination of the numerical integration step associated with the solutionof (1.2), the ? AO-algorithm gives much better persistent accuracy over time when applied to the datacontaining bias type errors, compared to the AO-algorithm.In § ? AO-algorithm. In § ? AO-algorithm. In § ? AO-algorithm and present a procedure fortaking the square root of (1.3). In § ? AO-algorithm.We do so by feeding in virtual accelerometer data, to which differing amounts of bias and noise typeerrors have been added, to both the ? AO and AO algorithms and comparing their predictions. Usingthose predictions in § ? AO-algorithm is less sensitive to bias type errors than theAO-algorithm. We make a few concluding remarks in § In this section we briefly recapitulate the mathematics and kinematics of rigid body motion from [27,§2] that are needed for the development of the proposed ? AO-algorithm.
Let E be a finite dimensional, oriented, Hilbert space, i.e., a Euclidean vector space. The Euclidean pointspace (cid:69) has E as its associated vector space. Let o P (cid:69) be (cid:69) ’s origin. The spaces E and (cid:69) are related toeach other such that for any point x P (cid:69) there exists a vector x P E such that o ` x “ x . The topological3 a Figure 1: Bias and noise type errors in acceleration component measurements. When there are errors δ X and δθ in accurately defining an accelerometer’s position and orientation, respectively (see (a)), biastype errors can occur in the measurement of acceleration components (see (b)). Noise type errors in theacceleration component measurements are usually a consequence of seismic, electrical, and other typesof noise.space (cid:66) serves as our model for a rigid body that executes its motion in (cid:69) . For that reason, we referto E and (cid:69) as the physical Euclidean vector space and point space, respectively. The spaces E R and (cid:69) R are another pair of Euclidean vector and point spaces, respectively, that are related to each other in thesame way that E and (cid:69) are related to each other. We refer to E R and (cid:69) R as the reference Euclideanvector and point spaces, respectively. The spaces E , (cid:69) , E R , and (cid:69) R have the same dimension, which wedenote as n sd . The dimension of (cid:66) is less than or equal to n sd . We call a select continuous, injective mapfrom (cid:66) into E R the reference configuration and denote it as κ R . The elements of (cid:66) are called materialparticles. We call X “ κ R p (cid:88) q , where (cid:88) P (cid:66) , the reference position vector of the material particle (cid:88) ,and we call the set κ R p (cid:66) q “ (cid:32) κ R p (cid:88) q P E R ˇˇ (cid:88) P (cid:66) ( the reference body (see Fig. 2). When we referto X as a material particle we in fact mean the material particle κ ´ p X q P (cid:66) . We model time as aone-dimensional normed vector space T and denote a typical element in it as τ “ τ s , where τ P R and s is a fixed vector in T of unit norm. We model the rigid body’s motion using the one-parameter familyof maps x τ : E R Ñ E (see Fig. 2). We call x τ the deformation map and x “ x τ p X q the material particle X ’s position vector at the time instance τ . The set κ τ p (cid:66) q = (cid:32) x τ p X q P E ˇˇ X P κ R p (cid:66) q ( (see Fig. 2) iscalled the current body. The sets p E i q i P (cid:73) and p e i q i P (cid:73) , where (cid:73) “ p , . . . , n sd q , are orthonormal sets of basis vectors for E R and E , respectively. By orthonormal we mean that the inner product between E i and E j , or e i and e j , where i, j P (cid:73) , equals δ ij , the Kronecker delta symbol, which equals unity iff i “ j and zero otherwise. Wecall X i the component of X w.r.t. E i iff X i “ X ¨ E i , where the dot denotes the inner product in E R .The dot in other expressions is to be similarly interpreted noting the space to which the vectors belong.We call the ordered set p X i q i P (cid:73) the component form of X w.r.t. p E i q i P (cid:73) and denote it as X or (cid:77) X . Wedenote the space of all m ˆ n real matrices, where m, n P N , (cid:77) m,n p R q ; here N and R denote the set ofnatural numbers and the space of real numbers, respectively. Thus, X P (cid:77) n sd , p R q . We access the i th component, where i P (cid:73) , of X , which of course is X i , as p X q i . Similarly, we denote the component of x w.r.t. e i as x i and call x “ p x i q i P (cid:73) P (cid:77) n sd , p R q the component form of x w.r.t. p e i q i P (cid:73) .Say W and U are two arbitrary, oriented, finite dimensional Hilbert spaces; for instance, they can be E R and E . We denote the space of all linear maps (transformations/operators) from W to U as (cid:76) p W , U q .We denote the norm of a vector w in W that is induced by W ’s inner product, i.e., p w ¨ w q { , as (cid:107) w (cid:107) .For u P U , the expression u b w denotes the linear map from W to U defined as p u b w q w “ u p w ¨ w q , (2.1) In our previous paper [27], we denoted the set of bounded linear operators from U to W as B p U , W q . As a linear operator on a finitedimensional normed space is automatically bounded, here we use (cid:76) p U , W q instead of B p U , W q to denote the set of all linear operators from U to W . w P W . If the sets p u i q i P (cid:73) and p w i q i P (cid:73) provide bases for U and W , respectively, then it can beshown that ´ p u i b w j q j P (cid:73) ¯ i P (cid:73) , which we will henceforth abbreviate as p u i b w j q i,j P (cid:73) , provides a basisfor (cid:76) p W , U q . The number T ij , where i, j P (cid:73) , is called the component of T P (cid:76) p W , U q w.r.t. u i b w j iff T ij “ u i ¨p T w j q . We call the nested ordered set p T ij q i,j P (cid:73) the component form of T w.r.t. p u i b w j q i,j P (cid:73) ,and denote it as (cid:77) T , or, when possible, briefly as T . We sometimes access the i th , j th component of T ,where i, j P (cid:73) , as p T q ij .Figure 2: Some mathematical quantities used in the description of motion. Illustration of the referenceEuclidean vector space E R , reference body κ R p (cid:66) q , a material particle X , the deformation map x τ ,current body κ τ p (cid:66) q , the (physical) Euclidean vector space E , and the location of the material particle X in E , i.e., the material particle X ’s spatial position vector x . See § (cid:73) . For example, the expression X i E i represents the sum ř i P (cid:73) X i E i . And anunrepeated index in a term will signify a set of n sd terms. For example, the term E i represents the set (cid:32) E i ˇˇ i P (cid:73) ( . For the case of rigid body motion x τ takes the form x τ p X q “ Q τ X ` c p τ q , (2.2)where Q τ is a proper (orientation preserving), linear isometry from E R into E and c p τ q “ c i p τ q e i ,where c i belongs to the space of twice continuously differentiable real valued functions over R , i.e.,to C p R , R q . The operator Q τ can be written as Q ij p τ q e i b E j , where Q ij P C p R , R q and satisfy Q ki p τ q Q kj p τ q “ δ ij for all τ P R . We abbreviate p Q ij p τ qq i,j P (cid:73) P (cid:77) n sd , n sd p R q , p c i p τ qq i P (cid:73) P (cid:77) n sd , p R q ,and p δ ij q i,j P (cid:73) P (cid:77) n sd , n sd p R q as Q p τ q , c p τ q , and I , respectively. The component or non-dimensional formof (2.2) is (1.1). Since Q τ is a proper isometry, it follows that Q p τ q , which we refer to as the rotationmatrix, belongs to the special orthogonal group SO p n sd q . As a consequence of belonging to SO p n sd q thematrix Q p τ q satisfies the equations Q T p τ q Q p τ q “ I , (2.3a)5nd Q p τ q Q T p τ q “ I , (2.3b)where Q T p τ q is the transpose of Q p τ q , i.e., Q T p τ q “ p Q p τ qq T .We call (cid:76) p T , E q the physical velocity vector space and denote it as V . It can be shown that the set p v i q i P (cid:73) , where v i P V and are defined such that v i τ “ τ e i , provides an orthonormal basis for V . Thevelocity of a material particle X executing its motion in E lies in V . The velocity of the material particle X at the instant τ , which we denote as V τ p X q , equals the value of the Fréchet derivative of the map T Q τ ÞÑ x X p τ q P E , where x X p τ q “ x τ p X q , at the time instance τ . Thus, it follows from (2.2) that V τ p X q “ L τ X ` c p τ q , (2.4)where L τ : “ Q ij p τ q v i b E j and c p τ q : “ c i p τ q v i , and Q ij and c i are the derivatives of Q ij and c i ,respectively. We abbreviate ` Q ij p τ q ˘ i,j P (cid:73) P (cid:77) n sd , n sd p R q and p c i p τ qq i P (cid:73) P (cid:77) n sd , p R q as Q p τ q and c p τ q ,respectively. Using (2.4) and (2.2) it can be shown that the velocity at the time instance τ of the materialparticle occupying the spatial position x P E at the time instance τ is W τ p x ´ c p τ qq ` c p τ q , where thelinear map W τ : E Ñ V is defined by the formula W τ x “ L τ Q ˚ τ x , (2.5)for all x P E . The operator Q ˚ τ is the Hilbert-adjoint of Q τ and is equal to Q ji p τ q E i b e j . Let thecomponent form of W τ w.r.t. p v i b e j q i,j P (cid:73) be p W ij p τ qq i,j P (cid:73) P (cid:77) n sd , n sd p R q , which we abbreviate as W p τ q . It follows from (2.5) that W ij p τ q “ Q ik p τ q Q jk p τ q , or equivalently, W p τ q “ Q p τ q Q T p τ q . (2.6)We call (cid:76) p T , V q the physical acceleration vector space and denote it as A . It can be shown that theset p a i q i P (cid:73) , where a i P A and are defined such that a i τ “ τ v i , provides an orthonormal basis for A .The acceleration of a material particle X executing its motion in E lies in A . The acceleration of X atthe time instance τ equals the value of the Fréchet derivative of the map T Q τ ÞÑ V X p τ q P V , where V X p τ q “ V τ p X q , at the time instance τ . Thus, it follows from (2.4) that A τ p X q “ M τ X ` c p τ q , (2.7)where the map M τ : E R Ñ A is defined by the equation M τ : “ Q ij p τ q a i b E j (2.8)and c p τ q : “ c i p τ q a i , where Q ij and c i are the derivatives of Q ij and c i , respectively. Let A τi p X q be thecomponent of A τ p X q w.r.t. a i . We abbreviate the ordered sets p A τi p X qq i P (cid:73) P (cid:77) n sd , p R q , ` Q ij p τ q ˘ i,j P (cid:73) P (cid:77) n sd , n sd p R q , and p c i p τ qq i P (cid:73) P (cid:77) n sd , p R q as A τ p X q , Q p τ q , and c p τ q , respectively.We will predominantly be presenting the ensuing results in component form. The component formcan be converted into physical or dimensional form. Therefore, from here on we will be often omit explic-itly using the qualification “is the component form of” when referring to the component form of a physicalquantity. For example, instead of saying “ A τ p X q as the component form of the acceleration of the materialparticle X at the time instance τ ”, we will often write “ A τ p X q is the acceleration of the material particle X at the time instance τ ”. The acceleration A τ p X q can be interpreted as the value of the (non-dimensional)acceleration field A τ : B R p (cid:66) q Ñ R , where we call B R p (cid:66) q : “ (cid:32) p X , X , X q P R ˇˇ X i E i P κ R p (cid:66) q ( thenon-dimensional reference body. For the definition of Fréchet derivative in the context of the current work see [27, §2.1] Review of the AO-algorithm
Let Q τ : A Ñ E R be defined by the equation Q τ “ Q ji p τ q E i b a j , (3.1)then we call the map A τ : κ R p (cid:66) q Ñ E R defined by the equation A τ p X q “ Q τ A τ p X q (3.2)the “Pseudo-acceleration field”. Say A τi p X q is the component of A τ p X q w.r.t. E i , then we abbreviate ` A τi p X q ˘ i P (cid:73) P (cid:77) n sd , p R q , the component form of A τ p X q w.r.t. E i , as A τ p X q . From the definitions ofthe pseudo acceleration field A τ (3.2), and A τ p X q , and the definitions of A τ p X q , and Q p τ q , which aregiven in § A τ p X q “ Q p τ q A τ p X q . (3.3)In [27, §2.1.1] it was shown that A τ p X q “ P p τ q X ` q p τ q , (3.4)where P p τ q “ Q T p τ q Q p τ q (3.5)is the component form of the linear map P τ : “ Q τ ˝ M τ w.r.t. p E i b E j q i,j P (cid:73) , and q p τ q P (cid:77) n sd , p R q isthe component form of q p τ q : “ Q τ c p τ q (3.6)w.r.t. p E i q i P (cid:73) . Thus, the acceleration field A τ is taken to be fully determined once Q p τ q , P p τ q , and q p τ q have been computed.Both the AO and the ? AO algorithms can be described as consisting of three primary steps. TheAO-algorithm’s three steps can briefly be described as follows:
AO-Step 1
Compute (time discrete versions of) the maps τ ÞÑ P p τ q and τ ÞÑ q p τ q using the measure-ments and the geometry of the arrangement of the four tri-axial accelerometers. AO-Step 2
Compute the map τ ÞÑ W p τ q , where W p τ q : “ Q T p τ q W p τ q Q p τ q , (3.7)using the map P computed in AO-Step 1 and numerical integrating (1.2). From Lemma B.1we have that the matrix W p τ q belongs to the space of n sd ˆ n sd real skew-symmetric matrices,which we denote as so p R , n sd q . AO-Step 3
Compute the map τ ÞÑ Q p τ q using the W map computed in AO-Step 2 and numericallyintegrating the equation Q p τ q “ Q p τ q W p τ q . (3.8)Equation (3.8) is from [27], where it appears as equation . .Step one of the ? AO-algorithm has two sub-steps: the predictor step and the corrector step. Thepredictor step is the same as
AO-Step 1 of the AO-algorithm. The corrector step is necessary for carryingout step two of the ? AO-algorithm. In step two of the ? AO-algorithm instead of obtaining W from(1.2), as is done in the AO-algorithm, we obtain it from (1.3). More precisely, in the ? AO-algorithm W p τ q is obtained as the square root of the symmetric part of P p τ q . We use sym p P p τ qq to denote thesymmetric part of P p τ q . The derivation of (1.3) is presented in B. A procedure for determining W p τ q asthe square root of sym p P p τ qq , i.e., for solving (1.3) for W p τ q with given P p τ q , is presented in § ? AO-algorithm is to compute Q using the W computed in step two. It involvesusing a slightly modified version of the numerical integration scheme described by equations . , . ,and . in [27, §3.2] to solve (3.8). We discuss it in § The ? AO-algorithm
As we mentioned in § ? AO algorithm consists of three primary steps. Those steps are as follows:Step 1 Compute (time discrete versions of) the maps τ ÞÑ P p τ q and τ ÞÑ q p τ q using the measurementsand the geometry of the arrangement of the four tri-axial accelerometers (see § . for details).Step 2 Use (1.3) and the P map obtained from Step 1 to solve for τ ÞÑ W p τ q . That is, for each τ in adiscrete sequence of time instances, compute W p τ q as the square root of the symmetric part of P p τ q (for details see § . ).Step 3 Compute (a time discrete version of) the map τ ÞÑ Q p τ q using the W map computed in Step 2and numerically integrating (3.8) (details in § ? AO-algorithm, Step 1 of In § . of [27] a method was presented to estimate P p τ q and q p τ q from the accelerometer measurementscorresponding to the time instance τ . Applying that method for each τ in a discrete time sequence yieldsa numerical approximation for the maps τ ÞÑ P p τ q and τ ÞÑ q p τ q . We present here an augmented versionof that method for computing similar numerical approximations. The primary difference between ourmethod and that presented in [27] is that the estimate for P p τ q yielded by our method is certain toretain some of the mathematical properties that are expected of P p τ q based on our theoretical analysis.Specifically, it follows from Lemmas C.1 and C.2 that sym p P p τ qq is a negative semidefinite matrix with itsnegative eigenvalues, if any, being of even algebraic multiplicities. These mathematical properties of P p τ q are critical for carrying out Step 2 of the ? AO-algorithm. We found that experimental noise and errorscan cause the estimate for P p τ q provided by the method presented in [27] to lose the aforementionedmathematical properties. Our method, on the contrary, ensures that the symmetric part of the estimated P p τ q is negative semidefinite and that its negative eigenvalues, when they exist, are of even algebraicmultiplicities. Once P p τ q is estimated, our method to estimate q p τ q is exactly the same as that in [27].We review it in 4.1.2. P p τ q Our method for estimating P p τ q can be described as consisting of two steps: a predictor step and a corrector step . In the predictor step we use the method presented in [27, §3.1] for estimating P p τ q tocompute a prediction for P p τ q . We denote this prediction as P p τ q p . In the corrector step we estimate P p τ q as the sum of P p τ q p and a correction term, which we construct using P p τ q p . The correction termis constructed such that the estimated P p τ q is as close as possible to P p τ q p under the constraint that theestimated P p τ q ’s symmetric part is negative semidefinite and its negative eigenvalues (if they exist) areof even algebraic multiplicities. Predictor step
Say the four tri-axial accelerometers are attached to the rigid body (cid:66) at the materialparticles ` (cid:108) (cid:88) ˘ (cid:96) P (cid:74) , where (cid:74) : “ p , . . . , q , and let the position vectors of those particles in E R , respec-tively, be ` (cid:108) X ˘ (cid:96) P (cid:74) (see Fig. 3). A tri-axial accelerometer is capable of measuring the components of itsacceleration in three mutually perpendicular directions. We refer to those directions as the accelerom-eter’s measurement directions. The measurement directions are usually marked on the accelerometerpackage by the manufacturer as arrows that are labeled x , y , and z . As (cid:66) moves in E , the attachedaccelerometers move with it, and, therefore, the measurement directions (in E ) can change with time.For an accelerometer (cid:96) , where (cid:96) P (cid:74) , we denote its time varying measurement directions in E using theorthonormal set ` (cid:108) e τi ˘ i P (cid:73) . Assuming that the accelerometers remain rigidly attached to (cid:66) , i.e., theirpositions and orientations w.r.t. (cid:66) do not change as (cid:66) moves in E , it can be shown that Q ˚ τ (cid:108) e τi , where (cid:96) P (cid:74) , i P (cid:73) , is a constant vector in E R , which we denote as (cid:108) E i . The position vectors ` (cid:108) X ˘ (cid:96) P (cid:74) and thedirections ` (cid:108) E i ˘ (cid:96) P (cid:74) ,i P (cid:73) are known from the arrangement and orientation of the accelerometers at theexperiment’s beginning. 8or (cid:96) P (cid:74) , let (cid:108) A p τ q : “ ` (cid:108) α j p τ q (cid:108) E j ¨ E i ˘ i P (cid:73) (no sum over (cid:96) ), where (cid:108) α i p τ q , i P (cid:73) , is the measurementreported by accelerometer (cid:108) X for the (non-dimensional) component of its acceleration in the (cid:108) e τi direc-tion at the time instance τ . And let (cid:108) X : “ ` (cid:108) X ¨ E i ˘ i P (cid:73) . Then, we compute P p τ q p as P p τ q is estimatedin [27] using the equation P p τ q p “ ´` i ∆ A p τ q ˘ T ` j ∆ X ˘¯ ´` i ∆ Y ˘ ` j ∆ Y ˘ T ¯ , (4.1)where i ∆ A p τ q : “ i ` A p τ q ´ A p τ q , i ∆ X : “ i ` X ´ X . The ordered sets i ∆ Y belong to (cid:77) n sd , p R q and aredefined by the equation ` ∆ Y , . . . , n sd ∆ Y ˘ “ ` ∆ X , . . . , n sd ∆ X ˘ ´ T , (4.2)where p¨q ´ T is the operator that acts on an invertable element of (cid:77) n sd , n sd p R q and returns the transposeof its inverse.Figure 3: Schematic of the locations and orientations of four tri-axial accelerometers (left) and theirmotion (right) (modified from [27], copyright 2020, Elsevier). Corrector step
In D.1 we show that sym p P p τ qq allows itself to be decomposed as N p τ q D p τ q N T p τ q , (4.3a)where N p τ q P (cid:77) n sd , n sd p R q is an orthogonal matrix, i.e., N T p τ q N p τ q “ I , (4.3b)and D p τ q P (cid:77) n sd , n sd p R q is a diagonal matrix that for n sd “ and , respectively, has the form D p τ q “ diag ` ´ λ p τ q , ´ λ p τ q ˘ and diag ` , ´ λ p τ q , ´ λ p τ q ˘ , (4.3c)where λ p τ q P R and the function diag p¨q : F n sd Ñ (cid:77) n sd , n sd p F q , where F is either R or C , is defined suchthat diag p a , . . . , a n sd q is a diagonal matrix with diagonal entries a , . . . , a n sd . Or to be mathematically precise, in the (cid:108) a τi P A direction that is defined such that ´ (cid:108) a τi s ¯ s “ (cid:108) e τi . p P p τ qq allowing the decomposition (4.3) is critical for carrying out Step 2 of the ? AO-algorithm. In an ideal scenario, in which there are no experimental errors or noise in the accelerometermeasurements, P p τ q p would be the same as P p τ q . However, due to the experimental noise and othererrors P p τ q p will generally be different from P p τ q . In general, such a deviation would not be of muchconsequence, since, experimental measurements of physical quantities, more often than not, are differentfrom the true values of those quantities. Thus, generally, we would, as done by [27], take P p τ q p tobe the final estimate for P p τ q and no longer distinguish between P p τ q p and P p τ q . However, in thepresent case the deviation of P p τ q p from P p τ q has an important consequence which requires us to nottake P p τ q p as P p τ q ’s final estimate. The important consequence is that in general sym p P p τ q p q will notallow a decomposition of the form (4.3). In general, it will only allow itself to be decomposed as N (cid:112) p τ q diag p λ i p τ q , . . . , λ n sd p τ qq N (cid:112) T p τ q , where N (cid:112) p τ q ’s columns are the eigenvectors of sym p P p τ q p q thatare chosen such that N (cid:112) p τ q is orthogonal and their corresponding eigenvalues λ i p τ q P R form a non-increasing sequence, i.e., λ p τ q ě . . . ě λ n sd p τ q . Therefore, instead of taking P p τ q p as the final estimateof P p τ q we derive the final estimate for P p τ q , as we detail next, using P p τ q p so that its symmetric partdoes allow a decomposition of the form (4.3).We take the skew-symmetric part of our final estimate for P p τ q to be the same as that of P p τ q p . Wetake its symmetric part to be N (cid:112) p τ q ˇ D p τ q N T (cid:112) p τ q , (4.4a)where ˇ D p τ q : “ diag ` , ´ ˇ λ p τ q , ´ ˇ λ p τ q ˘ , (4.4b)with ˇ λ p τ q : “ ´ λ p τ q` λ p τ q , λ p τ q ` λ p τ q ď , , λ p τ q ` λ p τ q ą , (4.4c)for n sd “ , and ˇ D p τ q : “ diag ` ´ ˇ λ p τ q , ´ ˇ λ p τ q ˘ , (4.4d)with ˇ λ p τ q : “ ´ λ p τ q` λ p τ q , λ p τ q ` λ p τ q ď , , λ p τ q ` λ p τ q ą , (4.4e)for n sd “ .The orthogonal matrix N (cid:112) p τ q and the eigenvalues λ i p τ q can be obtained from the spectral or symmetric-Schur [35, §8] decomposition of sym p P p τ q p q . Since sym p P p τ q p q is a real symmetric matrix, it is alwayspossible to carry out sym p P p τ q p q ’s spectral or symmetric Schur decomposition.To summarize, we take P p τ q ’s final estimate to be P p p τ q ` ∆ P p τ q , (4.5a)where ∆ P p τ q : “ sym ` ˇ P p τ q ´ P p p τ q ˘ , (4.5b) ˇ P p τ q : “ N (cid:112) p τ q ˇ D p τ q N T (cid:112) p τ q . (4.5c)10t can be ascertained that the symmetric part of our final estimate for P p τ q , namely ˇ P p τ q , allows adecomposition of the form (4.3). In fact, that decomposition is precisely the one given by (4.4). For n sd “ or it can be shown that ˇ P p τ q is the best approximation , in the Frobenius norm, to sym p P p p τ qq in the set of n sd ˆ n sd real symmetric negative-semidefinite matrices whose negative eigenvalues (whenthey exist) are of even algebraic multiplicities. q p τ q After estimating P p τ q as described in § q p τ q as (cid:108) ¯ A p τ q ´ P p τ q (cid:108) X , (4.6)where (cid:108) is some particular integer in (cid:74) . ? AO-algorithm, Step 2 of In D.2 we show using P p τ q ’s decomposition (4.3) that W p τ q can be computed from (1.3) as W p τ q “ ˘ N p τ q ‹ p λ p τ qq N T p τ q , n sd “ , N p τ q ‹ pp λ p τ q , , qq N T p τ q , n sd “ , (4.7)where the map ‹ p¨q is defined in A.Using the decomposition (4.4) for the symmetric part of our final estimate for P p τ q and similarcalculations as those used in § D.2, it can be shown that if we compute our final estimate for W p τ q as ˘ N (cid:112) p τ q ‹ ` ˇ λ p τ q ˘ N T (cid:112) p τ q , n sd “ , N (cid:112) p τ q ‹ `` ˇ λ p τ q , , ˘˘ N T (cid:112) p τ q , n sd “ , (4.8)then it and our final estimate for P p τ q will satisfy (1.3).We take the time discrete versions of the W and P maps to be constant over each time interval ∆ τ n : “ r n ∆ τ, p n ` q ∆ τ q , where n P p , , ¨ ¨ ¨ q and ∆ τ P R . We denote the values of these two mapsover ∆ τ n as W p n q and P p n q , respectively. The quantity W p q is known from initial conditions. For n ą we compute W p n q using (4.8). Using the positive and negative signs in (4.8) will give us two differentestimates for W p n q . Among those two estimates we choose the one that is closer to W ’s value over theprevious time interval. To be precise, we choose the estimate that gives a lower value for the metric m p W p n q , W p n ´ qq , where m : so p R , n sd q Ñ R , m p W p n q , W p n ´ qq “ arccos ˜ ‹ ` W p n q ˘ ¨ ‹ ` W p n ´ q ˘ (cid:107) ‹ ` W p n q ˘ (cid:107)(cid:107) ‹ ` W p n ´ q ˘ (cid:107) ¸ . (4.9)The metric (4.9) is a measure of the difference between between W p n q and W p n ´ q . Our criterionfor choosing between the two estimates given by (4.8) is essentially based on the assumption that inmost practical scenarios W should be a continuous function of time, and on the ansatz that due to thecontinuity of W the true W p n q would be the one that is closer to W p n ´ q when (cid:107) ‹ ` W p n ´ q ˘ (cid:107) is large.If W p n ´ q “ or when (cid:107) ‹ ` W p n ´ q ˘ (cid:107) is small, then our above criterion for choosing betweenthe two estimates for W p n q cannot be used. In such cases we first derive a prediction for W p n q byapplying the AO-algorithm to the previous time interval and then choose the estimate that is closer tothat prediction. ? AO-algorithm, Step 3 of We use a slightly modified version of the numerical integration scheme described by equations . , . , and . in [27, §3.2] to solve (3.8). As we did with W and P , we assume that the discrete We plan to publish this result along with its proof elsewhere. Q remains constant over each time interval ∆ τ n and denote its constant values as Q p n q , where n P p , , ¨ ¨ ¨ q . The matrix Q p q is taken to be known from the initial conditions of the experiment. For n ą the matrix Q p n q is computed as Q p n q “ Q p n ´ q e ∆ τ W n ´ , (4.10a)where the map e p¨q : so p R , n sd q Ñ SO p n sd q is defined by the equation e p¨q “ I ` sinc p (cid:107) ‹ p¨q (cid:107) q p¨q ` ˆ sinc ˆ (cid:107) ‹ p¨q (cid:107) ˙˙ p¨q , (4.10b)and W n ´ : “ ` W p n q ` W p n ´ q ˘ . (4.10c)The difference between the integration scheme (4.10) and that given by . , . , and . in[27, §3.2] is the manner in which W n ´ is computed. In Rahaman et al. ’s integration scheme, W n ´ is computed as W p n ´ q ` ∆ τ skew p P p n ´ qq , whereas we compute it using (4.10c). Here we useskew p P p n ´ qq to denote the skew-symmetric part of the matrix P p n ´ q . In silico validation, evaluation and comparison of the ? AO-algorithm
In this section, we check the validity and robustness of the ? AO-algorithm. We do that by feeding invirtual accelerometer data, to which differing amounts of bias and noise type errors have been added, tothe ? AO and AO algorithms, and comparing their resulting predictions. We discuss the creation of thevirtual accelerometer data in § § § The virtual accelerometer data we use for the comparison is from the numerical simulation of a rigidellipsoid impacting an elastic half-space. This simulation is presented and discussed in detail in [27],starting in § . However, for the readers convenience we give a very brief description of that simulationhere.In the simulation an ellipsoid, (cid:66) , is dropped onto an elastic half-space, H , under the action ofgravity with the initial angular and translational velocities prescribed (see Fig. 5). In the simulation theEuclidean point space (cid:69) , in which the ellipsoid and the half-space, respectively, execute their motionand deformation, is taken to be three dimensional, i.e., n sd “ . The vectors E i and e i , i P (cid:73) , aretaken to have units of meters and s to have units of seconds. Hence v i and a i , i P (cid:73) , have units ofmeters-per-second and meters-per-second-squared, respectively.The reference configuration of the ellipsoid is given in Fig. 4. In (cid:69) R the ellipsoid occupies the region ! p X , X , X q ˇˇ p X { a q ` p X { b q ` p X { c q ď ) , where p a, b, c q “ p . , . , . q . The half-spacewhen it is undeformed in (cid:69) occupies the region x ă . The initial location and orientation of (cid:66) w.r.t. H in (cid:69) are shown in Fig. 5. They correspond to the initial conditions Q p q “ diag p , , q and c p q “ p , , . q .The mechanics of H is modeled using the theory of small deformation linear elasto-statics and taking H ’s Young’s modulus and Poisson’s ratio to be Pa and . , respectively. The ellipsoid is rigid andhomogeneous. Its density and total mass are .
44 Kg { m and
10 Kg , respectively. This equation is the corrected version of equation . in [27, §3.2], which has two typos in it. (cid:66) and (cid:66) ’s interactionwith H ; and the torque exclusively from (cid:66) ’s interaction with H . The interaction between (cid:66) and H ismodeled using the Hertz contact theory, e.g., [36, 37]. For details on effecting a numerical solution tothe balance equations see [27, §B1.1]. For more details regarding the contact modeling see [27, §B1.2].Four virtual accelerometers are taken to be rigidly attached to the ellipsoid’s material particles (cid:108) (cid:88) , (cid:96) P (cid:74) . The locations and orientations of those accelerometers w.r.t. (cid:66) in (cid:69) R are shown in Fig. 4. Theirposition vectors (cid:108) X and orientations ` (cid:108) E i ˘ i P (cid:73) , (cid:96) P (cid:74) , are given in the caption of Fig. 4. The accelerationof any of the ellipsoid’s material particles can be obtained from the simulation results using the procedureoutlined in [27, §B2]. For i P (cid:73) , the values of α i , which is the component of (cid:88) ’s acceleration in the e τi direction, or to be more precise the a τi direction (see Fig. 5), at a sequence of time instances areshown in Fig. 6(a). The acceleration components (cid:108) α i , (cid:96) P (cid:74) , from the simulation do not contain any errors; other than, ofcourse, the errors that arise due to numerical discretization of the balance equations, numerical round-off, etc. However, those type of errors are of insignificant magnitude. Using the error free virtualaccelerometer data (cid:108) α i , (cid:96) P (cid:74) , from the simulation we generate virtual error-inclusive accelerometerdata (cid:108) α Error i , (cid:96) P (cid:74) , as (cid:108) α Error i p τ q “ (cid:108) α i p τ q ` η τ . (5.1)In equation (5.1) η τ denotes a particular realization of the Ornstein-Uhlenbeck (OU) process [38].We will describe shortly what we mean by a “realization”. The OU process is a continuous time and statestochastic process that is defined by the integral equation η τ ` τ ´ η τ “ β ż τ ` τ τ p µ ´ η τ q dτ ` σ ż τ ` τ τ dW τ , (5.2)where the second integral on the right is an Itô integral and W τ is the Wiener process [39]. The realnumber µ is called the mean value, σ ě the diffusion coefficient, and β ą the drift coefficient. Thesymbols τ , τ denote any two (non-dimensional) time instances. Since the OU process is a stochasticprocess, a given set of OU parameters, i.e., a particular set of µ , σ , and β values, define an entire familyor population of real valued functions on R . For a given OU parameter set, a particular realization of theOU process is obtained by drawing η from a Gaussian distribution of mean µ and variance σ { p β q andsolving (5.2). As a consequence of (5.1), (cid:108) α Error i , (cid:96) P (cid:74) , too are stochastic processes.For i P (cid:73) and (cid:96) P (cid:74) , when µ ‰ and σ “ , any particular realization of (cid:108) α Error i will contain only biastype errors. A representative realization of α Error2 for µ “ , σ “ , and β “ is shown in Fig. 6(b).Alternatively, when µ “ and σ ‰ any particular realization of (cid:108) α Error i will only contain noise typeerrors. A representative realization of α Error2 for µ “ , σ “ , and β “ is shown in Fig. 6(c).In general when µ and σ are both non-zero, realizations of (cid:108) α Error i will contain both bias and noise typeerrors. A representative realization of α Error2 for µ “ , σ “ , and β “ is shown in Fig. 6(d).From here on unless otherwise specified the value of β will always be equal to . ? AO and AO algorithms using virtual error-inclusive accelerom-eter data
We compare the predictions of the ? AO and AO algorithms for the following categories of OU parametersets.
Category I
Exclusively bias type errors: µ “ , . , . , . , and , and σ “ (see Table. 1). Category II
Exclusively noise type errors: µ “ , and σ “ , , , , and (see Table. 2). Category III
Both bias and noise type errors: µ “ , . , . , . , and , and σ “ (see Table. 3).13or a given OU parameter set we generate a large number of (cid:108) α Error i , (cid:96) P (cid:74) , realizations. We applythe ? AO and AO algorithms to each of those realizations and derive a population of predictions for theacceleration of the material particle (cid:88) (see Fig. 4). We denote the error-free (non-dimensional) accel-eration of (cid:88) , which we know from the rigid-ellipsoid-impact-simulation’s results, at the time instance τ as A p τ q P (cid:77) , p R q . The components of A p τ q , i.e., ` A p τ q ˘ i , i P (cid:73) , for a sequence of time instances are,respectively, shown in subfigures (a), (b), and (c) in each of Figs. 7–9. They are shown using thick graycurves.Since the predictions of the ? AO and AO algorithms are derived by, respectively, feeding the ? AOand AO algorithms the stochastic processes (cid:108) α Error i , (cid:96) P (cid:74) , they too, in fact, are stochastic processes.Representative realizations of the predictions from the ? AO (resp. AO) algorithm for different OUparameter sets are, respectively, shown in Figs. 7–9 in green (resp. red).
Among the OU parameter sets belonging to
Category I the set corresponding to the most amount of erroris p µ, σ q “ p . , . q . Representative realizations of the predictions from the ? AO and AO algorithmsfor this parameter set are, respectively, shown in green and red in Fig. 7. In Fig. 7 the realization ofthe ? AO-algorithm’s prediction appears to be more accurate than that of the AO-algorithm’s prediction,especially with increasing time. In order to make a more quantitative comparison between the ? AO andAO algorithms’ predictions, we focus on the time interval r , s and make use of the the metrics (cid:15) ´ ? AO ¯ : “ (cid:107) ? AO ` A ˘ ´ A (cid:107) (cid:107) A (cid:107) , (6.1) (cid:15) p AO q : “ (cid:107) AO ` A ˘ ´ A (cid:107) (cid:107) A (cid:107) , (6.2)where (cid:107) f (cid:107) : “ bş (cid:107) f p τ q (cid:107) dτ ; and R Q τ ÞÑ ? AO ` A ˘ p τ q P (cid:77) , p R q , and R Q τ ÞÑ AO ` A ˘ p τ q P (cid:77) , p R q are, respectively, particular realizations of the ? AO and AO algorithms’ predictions for R Q τ ÞÑ A p τ q P (cid:77) , p R q . The metric (cid:15) ` ? AO ˘ (resp. (cid:15) p AO q ) is constructed such that the smaller its valuethe more accurate the realization used in computing it. The values of (cid:15) ` ? AO ˘ and (cid:15) p AO q for therealizations shown in Fig. 7 are, respectively, . and . . The metric (cid:15) ` ? AO ˘ ’s smaller value incomparison to that of (cid:15) p AO q corroborates our earlier assertion that among the ? AO ` A ˘ and ? AO ` A ˘ shown in Fig. 7 the realization ? AO ` A ˘ is more accurate. This comparison between the ? AO and AOalgorithms’ predictions’ particular realizations prompts us to hypothesize that the ? AO algorithm is moreaccurate than the AO algorithm.In order to compare the ? AO and AO algorithms’ predictions in a more well-balanced and compre-hensive manner we calculated (cid:15) ` ? AO ˘ and (cid:15) p AO q , respectively, for a large number of realizations(population size N “ ) of the predictions from the ? AO and AO algorithms. The mean values ofthe thus generated populations of (cid:15) ` ? AO ˘ and (cid:15) p AO q are . and . , respectively (see rownumber 5 of Table. 1). The mean value of the population of (cid:15) ` ? AO ˘ being lower than the mean valueof the population of (cid:15) p AO q further supports our earlier hypothesis that ? AO-algorithm is more accuratethan the AO-algorithm.To recall, the discussion so far in this section exclusively relates to the p µ, σ q parameter set p . , . q .We performed analysis similar to the one discussed in the previous paragraph for the parameter sets p . , . q , p . , . q , p . , . q , and p . , . q as well. The means of (cid:15) ` ? AO ˘ ’s and (cid:15) p AO q ’s populationsfor these other parameter sets are, respectively, given in the first and second columns of Table. 1. It canbe seen from Table. 1 that the means of the (cid:15) ` ? AO ˘ populations are consistently smaller than thoseof (cid:15) p AO q populations across all the parameter sets considered. Furthermore, in Table. 1 the differencebetween the means of a (cid:15) ` ? AO ˘ population and a (cid:15) p AO q population corresponding to the same pa-rameter set increases with the amount of error, i.e., with the magnitude of µ in the present category.14hus for the category of exclusively bias type errors, in addition to the ? AO-algorithm appearing to bemore accurate than the AO-algorithm, it further appears that the ? AO-algorithm’s performance over theAO-algorithm increases with increasing amount of error. In Category II we consider the p µ, σ q parameter sets p . , . q , p . , . q , p . , . q , p . , . q , and p . , . q . The means of (cid:15) ` ? AO ˘ ’s and (cid:15) p AO q ’s populations for these parameter sets are, respec-tively, given in the first and second columns of Table. 2. It can be seen from Table. 2 that the meansof the (cid:15) ` ? AO ˘ populations are approximately the same as those of (cid:15) p AO q populations across all theparameter sets considered. Thus, for the category of exclusively noise type errors the ? AO-algorithmappears to perform on par with the AO-algorithm.Among the OU parameter sets belonging to
Category II the set corresponding to the most amountof error is p µ, σ q “ p . , . q . Representative realizations of the predictions from the ? AO and AOalgorithms for this parameter set are, respectively, shown in green and red in Fig. 8. In Fig. 8, atleast at the earlier time instances, the ? AO and AO algorithms’ predictions’ realizations are almostindistinguishable from one another. However, at later time instances the ? AO-algorithm seems to beperforming better than the AO-algorithm (This feature is likely not reflected in the results presentedin Table. 2 because they are calculated only using data from the initial time instances, or, to be moreprecise, from the r , s time interval.). Based on this observation we venture to speculate that even whenthe errors are predominantly of the noise type, the ? AO-algorithm will eventually begin to outperformthe AO-algorithm. In Category III we consider the p µ, σ q parameter sets p . , . q , p . , . q , p . , . q , p . , . q , and p , . q . The means of (cid:15) ` ? AO ˘ ’s and (cid:15) p AO q ’s populations for these parameter sets are, respectively,given in the first and second columns of Table. 3. Among the OU parameter sets belonging to Category III the set corresponding to the most amount of error is p µ, σ q “ p . , . q . Representative realizations ofthe predictions from the ? AO and AO algorithms for this parameter set are, respectively, shown in greenand red in Fig. 9.It can be seen from Table. 3 that the means of the (cid:15) ` ? AO ˘ populations are consistently smallerthan those of (cid:15) p AO q populations across all the parameter sets considered. Furthermore, in Table. 3 thedifference between the means of a (cid:15) ` ? AO ˘ population and a (cid:15) p AO q population corresponding to thesame parameter set increases with the amount of error, i.e., with the magnitudes of σ and µ . Thus, in Category III the relative performance of the ? AO and AO algorithms is very similar to that in
Category I .From the discussion in § ? AO-algorithm and from the discussion in § ? AO and AO algorithms are,approximately, equally sensitive to noise type errors. From the results in this section it appears that the ? AO-algorithm outperforms the AO-algorithm as long as the errors have some bias type component inthem, irrespective of what the amount of the noise type component in them is.15able 1: The mean and standard deviation of the error measure (cid:15) for 200 realizations of the accelerom-eter data only containing bias type error with σ “ . In this case, as the value of standard deviation isquite small compared to the value of mean, we do not show the value in the table. (cid:15) ˆ (mean ˘ std) µ ? AO-algorithm AO-algorithm .
01 0 . . .
87 41 . . .
73 83 . . .
22 220 .
281 87 .
88 471 . Table 2: The mean and standard deviation of the error measure (cid:15) for 200 realizations of the accelerom-eter data only containing noise type error with µ “ (cid:15) ˆ (mean ˘ std) σ ? AO-algorithm AO-algorithm0 .
01 0 .
021 1 . ˘ .
09 1 . ˘ . . ˘ .
64 11 . ˘ . . ˘ .
33 56 . ˘ . . ˘ .
11 112 . ˘ . Table 3: The mean and standard deviation of the error measure (cid:15) for 200 realizations of the accelerom-eter data containing bias and noise type errors with σ “ (cid:15) ˆ (mean ˘ std) µ ? AO-algorithm AO-algorithm . ˘ .
64 11 . ˘ . . . ˘ .
14 42 . ˘ . . . ˘ .
43 84 . ˘ . . . ˘ .
47 221 . ˘ .
021 88 . ˘ .
44 472 . ˘ . (cid:69) R the ellipsoid, (cid:66) , occupies the region ! p X , X , X q ˇˇ p X { a q ` p X { b q ` p X { c q ď ) , where p a, b, c q “ p . , . , . q . Four virtual accelerometers are, respectively, attached to the ellipsoid’smaterial particles (cid:108) (cid:88) , (cid:96) P (cid:74) . The reference position vectors of (cid:108) X , (cid:96) P (cid:74) , are, respectively, c E , b E , a E , and ´ a E , where p E i q i P (cid:73) are shown in the figure as well. The accelerometers’ orien-tations are given by ` (cid:108) E i ˘ i P (cid:73) , (cid:96) P (cid:74) . The component representation of ` E i ˘ i P (cid:73) w.r.t. p E i q i P (cid:73) is pp , , q , p , , q , p , , qq ; of ` E i ˘ i P (cid:73) is ´´ ´ ? , , ´ ¯ , ´ ? , , ´ ¯ , ´ , ? , ? ¯¯ ;of ` E i ˘ i P (cid:73) is ´ p , , q , ´ ? , , ´ ? ¯ , ´ ? , , ? ¯¯ ; and of ` E i ˘ i P (cid:73) is ´ p , , q , ´ ´ ? , , ´ ? ¯ , ´ ´ ? , , ? ¯¯ . We apply the ? AO-algorithm to the accelerometerdata from the four virtual accelerometers (cid:108) (cid:88) , (cid:96) P (cid:74) to predict the acceleration of the material particle (cid:88) . The reference position vector of (cid:88) is ´ c E . (modified from [27], copyright 2020, Elsevier)17 urface of the elastic half-space, Figure 5: Configuration of the rigid ellipsoid at different time instances in the simulation of it impactingan elastic half space (see §5.1 for details). In the simulation, the ellipsoid, (cid:66) , is dropped onto an elastichalf-space, H , under the action of gravity with the initial angular and translational velocities prescribed.The ellipsoid’s initial position in (cid:69) is prescribed by taking c p q “ p , , . q , and Q p q “ diag p , , q .Its initial velocities are prescribed by setting ‹p W p qq “ p , , q , and c p q “ p . , , q . (modified from[27], copyright 2020, Elsevier) 18 bad Figure 6: The acceleration components α i p τ q , i P (cid:73) , of the virtual accelerometer (cid:88) before and afteraddition of synthetic errors (see §5.2 for details). (a) shows the acceleration components before theaddition of synthetic errors. (b)–(d) show the error-inclusive acceleration component α Error2 , which isgenerated by adding different errors to the acceleration component α . In (b), (c), and (d) the errortime signals are particular realizations of the OU process for the OU parameter sets p µ, σ, β q “ p , , q , p , , q , and p , , q , respectively. The error in (b) corresponds to Category I (exclusively biastype errors); in (c) to
Category II (exclusively noise type errors); and in (d) to
Category III (a combinationof bias and noise type errors). 19 .2 0.4 0.6 0.8 1.0 1.2 1.4 - - - - - a bc Figure 7: Comparison of the predictions from the ? AO and AO algorithms for the acceleration of thematerial particle (cid:88) (see Fig. 4) in the rigid ellipsoid impact simulation (see §5.1 for details). Boththe ? AO and AO algorithms were fed the same virtual error-inclusive accelerometer data. The datawas generated by adding a particular realization of the OU process to the virtual accelerometer datafrom the rigid ellipsoid impact simulation. The OU realization corresponded to the OU parameter set p µ, σ, β q “ p , , q . Subfigures (a), (b), and (c), respectively, show the comparison for the componentof (cid:88) ’s acceleration in the e i , i P (cid:73) , directions. 20 .2 0.4 0.6 0.8 1.0 1.2 1.4 - - - - - a bc Figure 8: Comparison of the predictions from the ? AO and AO algorithms for the acceleration of thematerial particle (cid:88) (see Fig. 4) in the rigid ellipsoid impact simulation (see §5.1 for details). Both the ? AO and AO algorithms were fed the same virtual error-inclusive accelerometer data. The data wasgenerated by adding a particular realization of the OU process to the virtual accelerometer data from therigid ellipsoid impact simulation. The OU realization corresponded to the OU parameter set p µ, σ, β q “p , , q . Subfigures (a), (b), and (c), respectively, show the comparison for the component of (cid:88) ’sacceleration in the e i , i P (cid:73) , directions. 21 .2 0.4 0.6 0.8 1.0 1.2 1.4 - - - - - a bc Figure 9: Comparison of the predictions from the ? AO and AO algorithms for the acceleration of thematerial particle (cid:88) (see Fig. 4) in the rigid ellipsoid impact simulation (see §5.1 for details). Boththe ? AO and AO algorithms were fed the same virtual error-inclusive accelerometer data. The datawas generated by adding a particular realization of the OU process to the virtual accelerometer datafrom the rigid ellipsoid impact simulation. The OU realization corresponded to the OU parameter set p µ, σ, β q “ p , , q . Subfigures (a), (b), and (c), respectively, show the comparison for the componentof (cid:88) ’s acceleration in the e i , i P (cid:73) , directions. The predictions of the ? AO algorithm when fed just thevirtual accelerometer data, i.e., with no added errors, is also shown in (a), (b), and (c) using black opencircles. 22
Concluding remarks
1. The results discussed in § ? AO-algorithm provides a valid approach to determinethe complete motion of a rigid body using only data from four tri-axial accelerometers. How-ever, the ? AO-algorithm’s practical validity in the field still remains to be explored. In the future,we plan to conduct an experimental evaluation of the ? AO-algorithm to compliment its in silico validation that we presented in this paper.2. The comparison in § ? AO-algorithm is less sensitiveto bias type errors compared to the AO-Algorithm. However, we have not provided a mathematicalproof that the ? AO-algorithm is better than the AO-algorithm with regard to bias type errors.Thus, though the comparison presented in § ? AO-algorithm is less sensitive to bias type errors than the AO-algorithm, it by no means providesa proof for the hypothesis. A definitive resolution to the question of whether the hypothesis is truerequires an error analysis of both the ? AO-algorithm as well as the AO-algorithm. We currentlyhave not carried out such analyses. Nevertheless, irrespective of the relative merit of the ? AO overthe AO algorithm, it is quite clear from its derivation and the results discussed in § ? AO-algorithm retains all the benefits of the AO-algorithm. Both algorithms provide the com-plete motion of the rigid body in the fixed laboratory frame. Without integration or differentiation,both algorithms are able to determine the pseudo acceleration field, providing the magnitude ofacceleration for all material particles. Both algorithms can be applied to any arrangement of fourtri-axial accelerometers as long as they do not lie in the same plane. There is no restriction on theorientation of the tri-axial accelerometers.4. In the in silico validation, evaluation, and comparison of the ? AO-algorithm that we set up in § µ in the OU process as a measure of the bias type errors in the OUprocess’ realizations. There of course exist bias type errors that cannot be modeled in this manner.Thus, our evaluation of the relative sensitivities of the ? AO and the AO algorithms to bias typeerrors was carried out using a limited form of bias type errors. A more general method to representbias type errors in the virtual accelerometer data would provide a more comprehensive comparisonof the relative sensitivities of the ? AO and the AO algorithms to bias type errors.5. We expect the ? AO-algorithm to be especially useful for constructing inputs to the upcoming finiteelement based brain injury criteria. The finite element based brain injury criteria are based on themechanics of head motion and brain deformation, while the traditional brain injury criteria havemostly been developed empirically. Therefore, we expect the finite element based brain injurycriteria to find increased use in the future.
Funding information
The authors gratefully acknowledge support from the Panther Program and the Office of Naval Research(Dr. Timothy Bentley) under grants N000141812494 and N000142112044.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships thatcould have appeared to influence the work reported in this paper.
Acknowledgments
The authors thank Sayaka Kochiyama for her help in preparing some of the figures in the manuscript.23
Definition of the map ‹ p¨q
For n sd “ , the map ‹ p¨q : so p R , q Ñ R is defined by the equation ‹ p¨q “ p¨q . The inverse of ‹ p¨q isthe map ‹ ´ p¨q : R Ñ so p R , q defined by the equation ‹ p α q “ ˆ ´ αα ˙ . (A.1)For n sd “ , the map ‹ p¨q : so p R , q Ñ (cid:77) , p R q is defined by the equation ‹ p¨q “ pp¨q , p¨q , p¨q q .The inverse of ‹ p¨q is the map ‹ ´ p¨q : (cid:77) , p R q Ñ so p R , q defined by the equation ‹ pp α , α , α qq “ ¨˝ ´ α α α ´ α ´ α α ˛‚ . (A.2)To make our notation appear less cumbersome we denote ‹ ´ p¨q too as ‹ p¨q . Whether we mean ‹ p¨q or ‹ ´ p¨q will be clear from the argument of ‹ p¨q . B Derivation of (1.3) , i.e., proof of the statement that square of W p τ q is equal to the symmetric part of P p τ q The following lemmas can be shown to be equivalent to some of the standard results in the mechanics ofrigid solids, see the work in [40, §2.5.2] and [41, §6.4], which treats the rigid body motion in a moderncontinuum mechanics style; or see the work in [42, §9.4] and [43, §15], which treats the rigid bodymotion from a perspective of geometric mechanics. However, at a cursory level, due to our notation andformalism, those results might appear to be different from the below lemmas. The differences in notationand formalism are primarily due to the fact that in our work we distinguish between the vector spaces towhich the various physical quantities, e.g., the rotation operation, belong to and the (non-dimensional)matrix vector spaces to which the component representations of those quantities belong to. For thatreason, we believed that it would be helpful to the reader if we presented the following lemmas usingthe notation and formalism that we use in the current work.
B.1 Skew symmetry of W p τ q Lemma B.1.
The matrix W p τ q , defined in (3.7) , is skew-symmetric.Proof . It can be shown using W p τ q ’s definition (3.7) and equations (2.6) and (2.3a) that W p τ q “ Q T p τ q Q p τ q . (B.1)Differentiating (2.3a) we get ´ Q T ¯ p τ q Q p τ q ` Q T p τ q Q p τ q “ . (B.2)Noting that ´ Q T ¯ p τ q “ ` Q ˘ T p τ q we see that the first term on the left hand side of (B.2) is equal to ` Q ˘ T p τ q Q p τ q , which, in fact, is equal to the transpose of the second term on the left hand side of (B.2).Thus, it follows from (B.2) that sym ´ Q T p τ q Q p τ q ¯ “ . That is, that Q T p τ q Q p τ q is skew-symmetric.The result that W p τ q too is skew-symmetric now immediately follows from (B.1). B.2 Derivation of equation (1.3)
Lemma B.2.
The symmetric part of P p τ q is equal to the square of W p τ q . roof . Differentiating (2.3a) twice and rearranging we get that ´ Q T ¯ p τ q Q p τ q ` Q T p τ q Q p τ q “ ´ ´ Q T ¯ p τ q Q p τ q . (B.3)The equation (B.3) on noting that the first and second terms on its left hand side are in fact transposesof each other simplifies to sym ´ Q T p τ q Q p τ q ¯ “ ´ ´ Q T ¯ p τ q Q p τ q . (B.4)Writing the term on the right hand side of (B.4) as ´ ´ Q T ¯ p τ q I Q p τ q , and then using (2.3b) andreplacing the I in the resulting equation with Q p τ q Q T p τ q , we getsym ´ Q T p τ q Q p τ q ¯ “ ´ ˆ´ Q T ¯ p τ q Q p τ q ˙ ´ Q T p τ q Q p τ q ¯ . (B.5)Noting that ´ Q T ¯ p τ q “ ` Q ˘ T p τ q we see that the first factor on the right hand side of (B.5) is equalto ` Q ˘ T p τ q Q p τ q , which is the transpose of the second factor on the right hand side of (B.5), namely Q T p τ q Q p τ q . We, however, know from (B.1) that this second factor is equal to W p τ q . Thus, we get from(B.5) that sym ´ Q T p τ q Q p τ q ¯ “ ´ W T p τ q W p τ q , (B.6a)which simplifies on using Lemma B.1 tosym ´ Q T p τ q Q p τ q ¯ “ W p τ q . (B.6b)Using (3.5) and replacing the quantity Q T p τ q Q p τ q appearing on the left hand side of (B.6b) with P p τ q , we get sym p P p τ qq “ W p τ q . (B.7) C The matrix sym p P p τ qq is negative semidefinite and its negativeeigenvalues, if they exist, have even algebraic multiplicities The entries of W p τ q and sym p P p τ qq are all real numbers. However, in this section we consider W p τ q andsym p P p τ qq to be elements of (cid:77) n sd , n sd p C q , where (cid:77) n sd , n sd p C q is the space of all n sd ˆ n sd matrices whoseentries belong to C , the space of complex numbers. C.1 Negative semi-definiteness of sym p P p τ qq Lemma C.1.
The matrix sym p P p τ qq is negative semidefinite.Proof . The matrix sym p P p τ qq is self-adjoint since it is equal to its transpose, which, as all of sym p P p τ qq ’sentries are real, is equal to its conjugate-transpose, i.e., to its adjoint. Hence it follows from [44,7.31] that the matrix sym p P p τ qq is negative semidefinite iff the inner product x sym p P p τ qq X , X y , where X P (cid:77) n sd , p C q but is otherwise arbitrary, is always non-positive.It follows from (1.3) that x sym p P p τ qq X , X y “ x W p τ q X , X y . (C.1)25ince we know from Lemma B.1 that W p τ q is skew-symmetric, we can write W p τ q on the right handside of (C.1) as ´ W T p τ q W p τ q . On doing so and using the properties of the inner product we get x sym p P p τ qq X , X y “ ´x W p τ q X , W p τ q X y . (C.2)It also follows from the properties of the inner product that x W p τ q X , W p τ q X y is always non-negative.Therefore we get from (C.2) that x sym p P p τ qq X , X y is always non-positive, or equivalently that sym p P p τ qq is negative semidefinite. C.2 Form of the eigenvalues of sym p P p τ qq Lemma C.2.
The matrix sym p P p τ qq ’s negative eigenvalues, if they exist, have even algebraic multiplicities.Proof. Say S P (cid:77) n sd , n sd p C q then S is said to be normal when it commutes with its conjugate-transpose S H , i.e., when S S H “ S H S . Note that W p τ q W T p τ q “ W p τ q ` ´ W p τ q ˘ “ ` ´ W p τ q ˘ W p τ q “ W T p τ q W p τ q . (C.3)The first and the third equalities in (C.3) follow from the fact that W p τ q is skew-symmetric (Lemma B.1).Since all the entries of W p τ q are real, the transpose of W p τ q is equal to its conjugate-transpose. For thisreason it follows from (C.3) that W p τ q W H p τ q “ W H p τ q W p τ q , or equivalently that W p τ q is normal.Since W p τ q is normal it follows from the Complex Spectral Theorem e.g., see, [44, 7.24] that thereexists an unitary matrix U p τ q P (cid:77) n sd , n sd p C q such that U H p τ q W p τ q U p τ q “ diag p µ i p τ qq i P (cid:73) , where µ i p τ q P C and diag p µ i p τ qq i P (cid:73) is a diagonal matrix whose diagonal entries are µ p τ q , µ p τ q . . . , µ n sd p τ q . That is, U H p τ q W p τ q U p τ q “ ¨˚˚˚˝ µ p τ q . . . µ p τ q . . . ... ... . . . ... . . . µ n sd p τ q ˛‹‹‹‚ . (C.4)The complex numbers µ i p τ q , not necessarily distinct, are the eigenvalues of W p τ q and the columns of U p τ q are the eigenvectors of W p τ q . To be more specific, W p τ q u i p τ q “ µ i p τ q u i p τ q (no sum over i ) , (C.5)where u i p τ q “ ´ p U p τ qq ji ¯ j P (cid:73) . Applying the operation of complex-conjugation to both sides of (C.5) andnoting that the entries of W p τ q are all real we get that W p τ q u ˚ i p τ q “ µ ˚ i p τ q u ˚ i p τ q (no sum over i ) , (C.6)where µ ˚ i p τ q and u ˚ i p τ q are, respectively, the complex-conjugates of µ i p τ q and u i p τ q . It followsfrom (C.6) that if µ i p τ q is an eigenvalue of W p τ q then so is µ ˚ i p τ q . Thus p µ i p τ qq i P (cid:73) has the form ς p z p τ q , z ˚ p τ q , z p τ q , z ˚ p τ q , . . . , z k p τ q , z ˚ k p τ q , α p τ q , α p τ q , . . . , α l p τ qq , where ς p¨q is the permutation op-eration, z i p τ q P C with Im p z i p τ qq ‰ , z ˚ i p τ q is the complex-conjugate of z i p τ q , ď k ď t n sd { u , α i p τ q P R , and l “ n sd ´ k . It is not necessary that the complex numbers z i p τ q be distinct from oneanother. The same is the case with the real numbers α i p τ q .Taking the square on both sides of (C.4) and using our knowledge about the form of p µ i p τ qq i P (cid:73) weget that ´ U H p τ q W p τ q U p τ q ¯ “ diag ς ´ z p τ q , z ˚ p τ q , . . . , z k p τ q , z ˚ k p τ q , α p τ q , . . . , α l p τ q ¯ . (C.7)The expression on the left hand side of (C.7) can be simplified as ´ U H p τ q W p τ q U p τ q ¯ “ U H p τ q W p τ q U p τ q U H p τ q W p τ q U p τ q “ U H p τ q W p τ q U p τ q “ U H p τ q sym p P p τ qq U p τ q , where the second Here t ¨ u is the floor function. U p τ q is an unitary matrix, and the third equality from (1.3). Thus wecan get from(C.7) that U H p τ q sym p P p τ qq U p τ q “ diag ς ´ z p τ q , z ˚ p τ q , . . . , z k p τ q , z ˚ k p τ q , α p τ q , . . . , α l p τ q ¯ , (C.8a)which implies that u i p τ q are the eigenvectors of sym p P p τ qq as well, but with their corresponding eigen-values being ς ´ z p τ q , z ˚ p τ q , . . . , z k p τ q , z ˚ k p τ q , α p τ q , . . . , α l p τ q ¯ . We know from Lemma C.1 thatall of sym p P p τ qq ’s eigenvalues are non-positive. Therefore, each z i p τ q must be of the form λ i p τ q?´ ,where λ i p τ q P R and λ i p τ q ‰ , and α i p τ q “ . Hence, we get from (C.8a) that U H p τ q sym p P p τ qq U p τ q “ diag ς ` ´ λ p τ q , ´ λ p τ q , . . . , ´ λ k p τ q , ´ λ k p τ q , , . . . , ˘ , (C.9)where, to reiterate, ď k ď t n sd { u and λ i p τ q , when they exist, are non-zero and not necessarily distinct.It can be noted from this last assertion that all of sym p P p τ qq ’s negative eigenvalues, specifically thosecorresponding to λ i p τ q , are of even geometric multiplicities. For the case of symmetric matrices, algebraicand geometric multiplicities are one and the same. Therefore, sym p P p τ qq ’s negative eigenvalues, whenthey exist, are also of even algebraic multiplicities. D Calculating W p τ q as the square-root of sym p P p τ qq D.1 A spectral decomposition of sym p P p τ qq Since sym p P p τ qq is a real symmetric matrix it follows from the Real Spectral Theorem [44, 7.29] that itcan be decomposed as N p τ q D p τ q N T p τ q , (D.1)where D p τ q and N p τ q belong to (cid:77) n sd , n sd p R q . We first describe N p τ q and then D p τ q , in the next paragraph.The matrix N p τ q : “ p n p τ q , . . . , n n sd p τ qq T where n i p τ q P (cid:77) n sd , p R q are sym p P p τ qq ’s eigenvectors that areconstructed such that n i p τ q ¨ n j p τ q “ δ ij , or equivalently N p τ q N T p τ q “ I . (D.2)Using (C.9) it can be shown that for n sd “ , D p τ q “ p q ; for n sd “ , D p τ q “ diag p , q or diag ` ´ λ p τ q , ´ λ p τ q ˘ , where λ p τ q ‰ ; and for n sd “ , D p τ q “ diag p , , q ordiag ς ` ´ λ p τ q , ´ λ p τ q , ˘ . The last two results can be summarized by saying that when n sd “ , D p τ q “ diag ` ´ λ p τ q , ´ λ p τ q ˘ , (D.3)where λ p τ q P R , and when n sd “ , D p τ q “ diag ς p´ λ p τ q , ´ λ p τ q , q . Without loss of generality, wecan choose the order of n i p τ q so that their respective eigenvalues form a non-increasing sequence .Therefore, for concreteness in the case of n sd “ we take D p τ q “ ¨˝ ´ λ p τ q
00 0 ´ λ p τ q ˛‚ . (D.4) D.2 Calculation of W p τ q from sym p P p τ qq using (1.3) Let F p τ q : “ N T p τ q W p τ q N p τ q . (D.5)It can be shown using F p τ q ’s definition, equations (D.2), and (1.3), and sym p P p τ qq ’s decomposition thatis derived in D.1 and summarized in (4.3) that The matrix W p τ q as the square root of sym p P p τ qq does not depend on the order of n i p τ q . Different orders will lead to the same W p τ q . p τ q “ D p τ q . (D.6)Substituting D p τ q in (D.6) from (D.3) and (D.4) for then noting from Lemma B.1 and F p τ q ’s definitionthat F p τ q is skew-symmetric, it can be shown that for n sd “ and F p τ q “ ˘ ‹ p λ p τ qq (D.7)and F p τ q “ ˘ ‹ pp λ p τ q , , qq , (D.8)respectively. Equations . follow from (D.5), (D.7), and (D.8). ReferencesReferences [1] CDC, NIH, DoD, and VA Leadership Panel, Report to congress on traumatic brain injury in the unitedstates: Understanding the public health problem among current and former military personnel.,Tech. rep., Centers for Disease Control and Prevention (CDC), the National Institutes of Health(NIH), the Department of Defense (DoD), and the Department of Veterans Affairs (VA) (2013).[2] J. D. Cassidy, L. Carroll, P. Peloso, J. Borg, H. Von Holst, L. Holm, J. Kraus, V. Coronado, Incidence,risk factors and prevention of mild traumatic brain injury: results of the WHO collaborating centretask force on mild traumatic brain injury, Journal of Rehabilitation Medicine 36 (0) (2004) 28–60.[3] S. Ganpule, N. P. Daphalapurkar, K. T. Ramesh, A. K. Knutsen, D. L. Pham, P. V. Bayly, J. L. Prince,A three-dimensional computational human head model that captures live human brains dynamics,Journal of Neurotrauma 34 (13) (2017) 2154–2166.[4] A. C. Bain, D. F. Meaney, Tissue-level thresholds for axonal damage in an experimental model ofcentral nervous system white matter injury, Journal of Biomechanical Engineering 122 (6) (2000)615–622.[5] D. F. Meaney, B. Morrison, C. D. Bass, The mechanics of traumatic brain injury: a review of whatwe know and what we need to know for reducing its societal burden, Journal of BiomechanicalEngineering 136 (2) (2014) 021008.[6] S. Kleiven, Why most traumatic brain injuries are not caused by linear acceleration but skull frac-tures are, Frontiers in Bioengineering and Biotechnology 1 (2013) 15.[7] E. Bar-Kochba, M. T. Scimone, J. B. Estrada, C. Franck, Strain and rate-dependent neuronal injuryin a 3D in vitro compression model of traumatic brain injury, Scientific Reports 6 (2016) 30550.[8] P. K. Kriz, S. J. Staffa, D. Zurakowski, M. MacAskill, T. Kirchberg, K. Robert, J. Baird, G. Lockhart,Effect of penalty minute rule change on injuries and game disqualification penalties in high schoolice hockey, The American Journal of Sports Medicine 47 (2) (2019) 438–443.[9] R. Mannix, N. J. Morriss, G. M. Conley, W. P. Meehan III, A. Nedder, J. Qiu, J. Float, C. A. DiCesare,G. D. Myer, Internal jugular vein compression collar mitigates histopathological alterations afterclosed head rotational head impact in swine: A pilot study, Neuroscience 437 (2020) 132–144.[10] J. H. McElhaney, J. L. Fogle, J. W. Melvin, R. R. Haynes, V. L. Roberts, N. M. Alem, Mechanicalproperties of cranial bone, Journal of Biomechanics 3 (5) (1970) 495 – 511.[11] R. Delille, D. Lesueur, P. Potier, P. Drazetic, E. Markiewicz, Experimental study of the bone be-haviour of the human skull bone for the development of a physical head model, InternationalJournal of Crashworthiness 12 (2) (2007) 101–108.2812] J. A. Motherway, P. Verschueren, G. Van der Perre, J. Vander Sloten, M. D. Gilchrist, The mechanicalproperties of cranial bone: The effect of loading rate and cranial sampling position, Journal ofBiomechanics 42 (13) (2009) 2129 – 2135.[13] S. Budday, R. Nay, R. de Rooij, P. Steinmann, T. Wyrobek, T. C. Ovaert, E. Kuhl, Mechanical prop-erties of gray and white matter brain tissue by indentation, Journal of the Mechanical Behavior ofBiomedical Materials 46 (2015) 318 – 330.[14] R. Willinger, H. Kang, B. Diaw, Three-dimensional human head finite-element model validationagainst two experimental impacts, Annals of Biomedical Engineering 27 (3) (1999) 403–410.[15] R. M. Wright, A. Post, B. Hoshizaki, K. Ramesh, A multiscale computational approach to estimatingaxonal damage under inertial loading of the head, Journal of Neurotrauma 30 2 (2013) 102–18.[16] C. W. Gadd, Use of a weighted-impulse criterion for estimating injury hazard, Tech. rep., SAETechnical Paper (1966).[17] C. C. Chou, G. W. Nyquist, Analytical studies of the head injury criterion (HIC), SAE Transactions(1974) 398–410.[18] J. A. Newman, A generalized acceleration model for brain injury threshold (GAMBIT), in: Proceed-ings of the 1986 International IRCOBI Conference on the Biomechanics of Impact, 1986.[19] E. G. Takhounts, M. J. Craig, K. Moorhouse, J. McFadden, V. Hasija, Development of brain injurycriteria (BrIC), Tech. rep., SAE Technical Paper (2013).[20] K. Laksari, M. Fanton, L. C. Wu, T. H. Nguyen, M. Kurt, C. Giordano, E. Kelly, E. O’Keeffe, E. Wallace,C. Doherty, M. Campbell, S. Tiernan, G. Grant, J. Ruan, S. Barbat, D. B. Camarillo, Multi-directionaldynamic model for traumatic brain injury detection, Journal of Neurotrauma 37 (7) (2020) 982–993.[21] L. F. Gabler, J. R. Crandall, M. B. Panzer, Development of a metric for predicting brain strainresponses using head kinematics, Annals of Biomedical Engineering 46 (7) (2018) 972–985.[22] L. F. Gabler, J. R. Crandall, M. B. Panzer, Development of a second-order system for rapid estimationof maximum brain strain, Annals of Biomedical Engineering 47 (2019) 1971–1981.[23] F. Hernandez, L. C. Wu, M. C. Yip, K. Laksari, A. R. Hoffman, J. R. Lopez, G. A. Grant, S. Kleiven,D. B. Camarillo, Six degree-of-freedom measurements of human mild traumatic brain injury, Annalsof Biomedical Engineering 43 (8) (2015) 1918–1934.[24] D. B. Camarillo, P. B. Shull, J. E. Mattson, R. Shultz, D. Garza, An instrumented mouthguard formeasuring linear and angular head impact kinematics in american football, Annals of BiomedicalEngineering 41 (2013) 1939–1949.[25] J. Chu, J. Beckwith, J. Crisco, R. Greenwald, A novel algorithm to measure linear and rotationalhead acceleration using single-axis accelerometers, Journal of Biomechanics 39 (Supplement 1)(2006) S534.[26] S. Rowson, G. Brolinson, M. Goforth, D. Dietter, S. Duma, Linear and angular head accelerationmeasurements in collegiate football, Journal of Biomechanical Engineering 131 (6) (2009) 061016.[27] M. M. Rahaman, W. Fang, A. L. Fawzi, Y. Wan, H. Kesari, An accelerometer-only algorithm fordetermining the acceleration field of a rigid body, with application in studying the mechanics ofmild traumatic brain injury, Journal of the Mechanics and Physics of Solids 143 (2020) 104014.[28] W. Fang, Y. Wan, H. Kesari, AO-Desktop-App, https://github.com/AppliedMechanicsLab/AO-Desktop-App.githttps://github.com/AppliedMechanicsLab/AO-Desktop-App.git