Capabilities of helmets for preventing head injuries induced by ballistic impacts

The limiting performance of ballistically loaded helmets designed to reduce head injuries is studied analytically. The projectile does not penetrate the helmet. This analysis evaluates the absolute minimum of the peak displacement of the helmet shell relative to the head, provided that criteria measuring the severity of head injuries lie within prescribed limits. Rather than optimize a specific design configuration, e.g. a viscoelastic foam liner, characteristics of a time-dependent force representing the helmet liner are calculated. The formulation reduces the limiting performance analysis to an optimal control problem.


Introduction
Helmets are widely used to prevent impact-induced head injuries in traffic, industry, construction work, sports, and other areas of human activity.Typical examples are helmets for motorcyclists, bicycle riders, soldiers, and ice hockey players and boxers.The design of a helmet, including the selection of materials, depends primarily on the characteristics of the impact disturbance of the helmet's outer shell with a fixed obstacle (for instance, with a roadway surface when a bicycle or motorcycle rider falls) or with a projectile (for example, with a bullet during a military operation or a puck when playing ice hockey).On the one hand, a good helmet must have shock isolation properties that correspond to head injuries lying within tolerable limits.On the other hand, the helmet must be fairly compact and light to be convenient for workers, soldiers or sportsmen to perform their normal functions when no impact has occurred.These requirements contradict each other, thereby giving rise to the problem of optimal design of helmets for head injury prevention.
To calculate the optimal design variables of a helmet, it is necessary to have a mathematical model that adequately describes the mechanical response of a helmeted head to an impact of a given type, as well as a set of functionals characterizing the severity of injuries in terms of this mathematical model.Then the determination of the optimal design variables is reduced to an optimal control problem or a mathematical programming problem.At present, a number of different models are available that enable one to estimate the nature and severity of head injuries depending on the parameters of the impact and shock isolation properties of helmets.For years the debate on the fundamental cause of brain injury has been between linear and angular acceleration.Now both of these are being questioned and other causes, such as strain rate, are being considered [1].The shock isolation properties depend on the mechanical characteristics of the liner that isolates the head from the helmet's shell and serves to cushion impacts.The head/brain injury models tend to be complicated, especially when nonlinear finite-element models are used.The simpler models are represented by a small number of concentrated masses or rigid bodies connected by springs and dampers in various combinations.In these models, the bodies describe the inertial properties of the skull and brain, whereas the springs and dampers represent their viscoelastic properties.Numerical parameters of the models are identified on the basis of experimental data.In these simple models, the severity of injuries is characterized by such quantities as the maximum absolute value of the acceleration of individual bodies of the model, maximum deformation of the springs, and the power or work of elastic and damping forces.The maximum values allowed for these characteristics are determined on the basis of biomechanical and medical data enabling one to establish relationships between the mechanical characteristics of the head's response to an impact and the severity of injuries.A number of mathematical models, which do not take angular motion into account, widely utilized for investigating the head's response to impact loads are presented in [2][3][4].
A helmet for head injury prevention is in essence a shock isolator and, hence, its design fits into the conventional design scheme for such systems.The first stage of this scheme involves a limiting performance analysis to evaluate the absolute minimum of the performance index of the system to be designed irrespective of the specific hardware implementation.For the limiting performance analysis, the isolation device, the shock-isolating liner, and in some cases, the helmet itself, are treated as generators of a control force that acts between the base (a body to which the impact load is directly applied) and the object to be protected.The desired estimate of the performance index is determined as its greatest lower bound over all possible laws of change of the control force.Thus, the limiting performance analysis of shock isolation systems is reduced to an optimal control problem.The current state of the art in the theory of optimal shock and impact isolation is presented in [5].A special section of this book is devoted to the limiting performance analysis of impact isolation systems for the prevention of spinal and head injuries.The book contains a lengthy bibliography on the optimal protection from impact and shock and related optimal control problems.
In the present paper, we study the limiting capabilities of helmets for the prevention of head injuries induced by impacting a helmet with a non-penetrating projectile.As a model of the head's response to an impact, we utilize the translational head injury model proposed by Stalnaker and colleagues [2][3][4].The performance index to be minimized is the peak value of the post-impact displacement of the outer shell of the helmet with respect to the head (the rattlespace of the helmet).The injury criteria subject to constraints are the peak power of the forces caused by the skull bone deformation and the peak acceleration of the brain.The action of the impact-isolating liner is modeled by the control force acting between the helmet shell and the skull.Earlier, a similar study was performed for a helmet to prevent injuries induced by an impact of a helmeted head against a fixed obstacle [6].

Mathematical model of a helmeted head interacting with a non-penetrating projectile
A sketch illustrating a helmeted head impacted by a projectile is shown in Fig. 1.The helmet consists of a rigid shell and a relatively soft liner that isolates the head from the impact-induced shock disturbance thereby reducing the severity of possible injury.To describe the response of a head to a shock we will use a two-mass translational head injury model suggested by Stalnaker, Low, and Lin [2].In the case of a high-speed projectile, there is little reason to believe that angular motion will play an important role.This model involves two masses (m 1 and m 2 ), a damper c 2 , as well as a combination of a spring k and a damper c 1 connected in series.See Fig. 2. The mass m 1 represents the mass of the skull bone interacting with the helmet liner.The mass m 2 represents the mass of the brain and the remaining portion of the head that may move with respect to the skull bone.The stiffness coefficient k and the damping factor c 1 model primarily the elastic and dissipative properties of the skull, whereas the damping factor c 2 models primarily the dissipative properties of the brain.This model accounts only for the translational motion of the head.Its parameters should be adjusted depending on the direction of the impact.For more details, see [2].
The mass m 4 in Fig. 2 represents that of the helmet shell, the force f is due to interaction of the projectile with the shell, and the force u characterizes the interaction of the helmet liner with the skull bone.The force u depends on the design of the helmet.Consider the model of the interaction of the projectile with the helmet shell.Let m p be the mass of the projectile and v p its velocity when hitting the helmet.We assume that the impact is perfectly inelastic, i.e., that the projectile adheres to the shell after the impact.In this case, using the law of conservation of momentum to calculate the velocity of the helmet shell (mass m 4 ) with the projectile embedded in it we obtain This relation is valid if the impulse of the force applied by the liner to the shell during the impact is much less than the momentum of the impacting projectile, m p v p .
An impact-induced force acting on an impacted body is frequently represented as the half-sine pulse where F 0 is the amplitude of the force and T 0 is the duration of the impact.The amplitude F 0 can be expressed in terms of the quantities m p , m 4 , v p , and T 0 .Let x and y be the coordinates of the impacted body (helmet shell) and the projectile, respectively, in a fixed inertial reference frame.The motion of these bodies is governed by the equations with the initial conditions Solving this initial value problem yields Since at the time instant T 0 the projectile adheres to the helmet shell, we have ẋ(T 0 ) = ẏ(T 0 ).Using this condition and Eq. ( 6) we obtain the desired relation If the mass of the helmet shell considerably exceeds that of the projectile, one may use the simplified relation The distance moved by the helmet shell during the impact is measured by If the projectile rebounds from the shell, the right-hand sides of Eqs ( 1), ( 7), ( 8) and ( 9) should be multiplied by 1 + µ, where µ is the restitution coefficient of the impact velocity.To obtain these relations, one should replace the condition ẋ(T 0 ) = ẏ(T 0 ) at the end of the impact interaction process by the condition ẋ(T 0 ) = ẏ(T 0 ) + µv p .
If x(T 0 ) is much smaller than a characteristic length of the system, for example, the thickness of the liner, one can neglect the duration of the impact and utilize the model of an instantaneous impact with where V is the post impact velocity of the helmet shell defined by Eq. ( 1) and δ(t) is Dirac's delta function.The equations of motion of the system have the form where x 1 , x 2 , and x 4 are the coordinates of the bodies representing the skull bone, the brain, and the helmet shell, respectively; x 3 is the coordinate of the point at which the spring k is connected to the damper c 1 .All these coordinates are measured with respect to a fixed inertial reference frame.The coordinate system is chosen such that when all of the coordinates are equal to zero (x i = 0, i = 1, 2, 3, 4), the brain, the skull bone, and the helmet are not stressed.We subject these equations to zero initial conditions The initial value for the velocity ẋ3 is not prescribed in Eq. ( 15), because it is uniquely defined by the initial values of x 1 , x 3 , and ẋ2 in accordance with Eq. ( 13).The initial conditions of Eq. ( 15) correspond to the case where the system is in a state of rest and is unstressed at the time of the impact.

Head injury criteria
To evaluate the severity of head injury in terms of the model described in the previous section, we will use the following two functionals: where T is the time interval on which the motion of the system is considered.
Criterion J 1 measures the peak value of the power developed in the spring k during the interval [0, T ].It represents the power developed by the forces due to elastic deformation of the skull.There is a relationship between the value of J 1 and the probability of skull fracture [3].
Criterion J 2 measures the peak acceleration of the mass m 2 during the interval [0, T ].It can be regarded as the peak acceleration of the head, since m 2 is about 90% of the total mass of the head.This criterion characterizes the severity of brain concussion [3].
To evaluate and optimize the design of the helmet, we supplement the head injury criteria with the criterion measuring the peak displacement of the helmet shell relative to the skull This criterion, the so-called rattlespace, characterizes the thickness of the liner of the helmet.

Limiting performance analysis. Problem formulation
The main idea and purpose of the limiting performance analysis of a shock isolation system is to investigate the theoretical possibility that all performance criteria lie within admissible ranges, thus guaranteeing operational integrity of an object or survivability of a person to be protected.In this analysis, the configuration and engineering implementation of the isolation system is not taken into account, and the control force can be treated as a function of time.Frequently, the limiting performance analysis is treated as an optimal control problem aiming at minimizing one of the performance criteria, provided that the other criteria are subjected to prescribed constraints.This approach is especially appropriate if the constraint to be imposed on one of the criteria is not established reliably or this criterion is critical for the design.Consider the following formulation of the limiting performance problem.Problem 1.For the system governed by Eqs (11)-( 14) with the initial conditions of Eq. (15), find the greatest lower bound J 0 3 for the criterion J 3 , provided that the other criteria are constrained by prescribed constants, i.e., where D i are specified positive numbers.

Transformation of the equations of motion
Some transformations will be introduced to simplify the equations of motion in order to facilitate the limiting performance analysis.First, proceed from the coordinates (x 1 , x 2 , x 3 , x 4 ) to the new coordinates (z, y 1 , y 2 , y 3 , y 4 ) defined by The variable z is the coordinate of the center of mass of the system, and y i are the coordinates of the bodies of this system measured from the center of mass.In terms of the new variables, the equations of motion take the form where The initial conditions for the system of Eqs ( 21)-( 24) corresponding to Eq. ( 15) are given by Note that the system of Eqs ( 21)-( 24) does not involve the variable y 4 .Use Eq. ( 20) to express this variable in terms of y 1 , y 2 , and y 3 as Introduce the notation to reduce Eqs.( 21)-(24) to the system with the initial conditions The performance criteria of Eqs ( 16)-(28) represented in terms of the new variables take the form The structure of Eqs ( 29)-( 33) and ( 35)-(37) implies that for the limiting performance analysis, in which one is interested primarily in the values of the performance criteria rather than in the behavior of the control force, it suffices to use only three equations (Eqs (31)-( 33)) considering w as a new (auxiliary) control variable.Once the control function w(t) is known, the system of Eqs (31)-(33) can be solved to find y 2 (t), ξ(t), and η(t).Then, if necessary, one can determine the control force u(t) from Eq. (30).
To further simplify the governing relations, introduce the dimensionless (primed) variables where In what follows, the primes labelling the dimensionless variables are omitted.The nondimensionalization reduces Eqs (31)-(33) to the form where The performance criteria become For Problem 1 formulated in terms of the dimensionless variables, we have D 1 = 1 and D 3 = 1 in Eq. (19).

Simplified optimal control problem
In addition to Problem 1, consider a simplified problem in which the performance index J 3 of Eq. ( 46) is replaced by Problem 2. For the system of Eqs (40)-(42) subjected to the initial conditions of Eq. (34), find the greatest lower bound of the performance index J3 under the following constraints With reference to Eq. ( 41), the constraint of Eq. ( 49) can be rewritten as Note that when solving Problem 2, Eq. ( 42) can not be taken into account.This is the case because the performance index and constraints of this problem involve only the variables y 2 , ξ and the control function w.For given control function w(t), the variables y 2 and ξ are completely and uniquely defined by Eqs (40) and (41) subject to the respective initial conditions.
The error of using J3 to approximate the criterion J 3 can be estimated as This estimate is based on an upper estimate of the variable η and the familiar inequalities From Eq. ( 52) it follows that and, hence, The solution of the system of Eqs ( 41) and ( 42) with zero initial conditions yields The last relation and the constraint of Eq. (48) imply the inequalities By substituting the notation of Eqs ( 46) and (47) into Eq.( 54) and using the inequality of Eq. ( 57) we obtain the estimate of Eq. ( 51).
Note that this estimate depends on the time T during which the motion is considered.The estimate is simplified and becomes independent of T in the case of α = β.In this case, Solve Problem 2 for the case of the instantaneous impact where the impact force is defined as with V given by Eq. ( 1).The minus sign on the right-hand side was taken for convenience and does not influence the final solution of the problem.In terms of the dimensionless variables, In this case, Eq. ( 40) with the initial conditions y 2 (0) = 0 and ẏ2 (0) = 0 is equivalent to the equation with the initial conditions

Lower estimate for the performance index in Problem 2
To obtain a lower estimate for the performance index we will solve Problem 2 on the time interval 0 t T where T is defined as the first instant at which the velocity of body m 2 vanishes, i.e, ẏ2 (T ) = 0.The time T is not prescribed beforehand but is to be determined during the solution.
The optimal control in this case takes on the maximum value allowed by the constraints of Eqs ( 48) and ( 50) at each instant.The control variable, being represented as a function of ξ, has the form To prove the optimality of this control, represent Eq. ( 41) in the form Let ξ w (t) be the solution of the differential equation of Eq. ( 64) subject to the initial condition ξ(0) = 0 for the admissible control w.Let w 0 (ξ) be the control defined by Eq. (63).By construction, w 0 (ξ) is the maximum value of the control variable admissible for the given ξ.Hence, for any admissible control w, the inequality f w0 (ξ) f w (ξ) holds.Therefore, in accordance with Chaplygin's theorem on differential inequalities [7,8], By substituting the left-hand side of Eq. ( 41) for w into Eq.( 61) we obtain ÿ 2 = − ξ − αξ.Integrate this relation over the interval [0, t] with the initial conditions of Eq. (34) for ξ and Eq.(62) for y 2 to obtain From Eqs (65) and (66) we conclude that ẏ(w0) where y (w0) 2 (t) and y (w) 2 (t) are the solutions of Eq. ( 61) for the control laws w 0 and w, respectively.Integration of the inequality of Eq. (67) from 0 to t with the initial conditions y    To determine the minimal value of the criterion J3 one should define the functions ẏ2 (ξ) and y 2 (ξ) by integrating the system (obtained from the system of Eqs ( 61) and (41) by utilizing the independent variable ξ) with the control of Eq. ( 63) and the initial conditions Then solve the equation ẏ2 (ξ) = 0 to find ξ * .Finally, calculate J3 = y 2 (ξ * ).
The solution of Eq. ( 68) for α 1/4 is given by where For α 1/4, the solution has the form In this case, ξ * and J3 are expressed by the simple relations In the case of α < 1/4, these quantities have to be calculated numerically.
If ξ * α −1/2 , then the value J3 = y 2 (ξ * ) is attained.For example, this is the case for the control coinciding with that of Eq. ( 63) for 0 t T and being identically zero for t > T .It is apparent that such a control does not violate the constraint of Eq. ( 48).It will be shown that this control also does not violate the constraint of Eq. ( 50).On the interval 0 t T , this constraint holds in accordance with the construction of the control of Eq. (63).Integrate Eq. (41) on the interval T t < ∞ for w = 0 and ξ(T ) = ξ * to obtain ξ = ξ * e −α(t−T ) . (75) Then for the left-hand side of the inequality of Eq. (50), Therefore, the inequality of Eq. (50) holds for ξ * α −1/2 .If ξ * > α −1/2 , the constraint of Eq. ( 50) is violated at the time instant T after switching the control of Eq. ( 63) to w = 0.

Upper estimate for the performance index in Problem 2
If we apply the control of Eq. (63) to the system of Eqs ( 40) and (41) on the interval 0 t T and define w ≡ 0 for t > T, the criterion J3 will coincide with the lower bound calculated in the previous section.However, if ξ * > α −1/2 , the constraint of Eq. (50) will be violated.In fact, the left-hand side of Eq. (50) becomes αξ 2 for w = 0.When switching the control of Eq. ( 63) to zero at the instant T, we have αξ 2 = αξ 2 * > 1 immediately after the switch.
In this section, we will modify the control of Eq. (63) to assure that the constraint of Eq. ( 50) is satisfied during the infinite time interval.Denote the time history of the control of Eq. (63) by w(t).To define the function w(t) one should (i) substitute the control of Eq. (63) into Eq.( 41), (ii) integrate this equation with the initial condition ξ(0) = 0 to determine ξ = ξ(t) on the time interval 0 t T , and (iii) substitute the function ξ(t) into Eq.(63).Define the modified control w * (t) as follows: where The function w * (t) is defined by the condition that ξ ξ = −1 during the interval t 1 < t < T 1 (see Eq. ( 49)).Integration of the equation ξ ξ = −1 with the initial condition ξ(t 1 ) = ξ yields By substituting this function into the left-hand side of Eq. ( 41) we obtain Eq. ( 78).
The time instants t 1 and T 1 in Eq. ( 77) are to be defined by the conditions The control of Eqs (77), ( 78) and (80) satisfies the constraints of Eqs ( 48) and (50).Apparently, both of these constraints hold on the interval 0 t t 1 , since the control of Eq. ( 77) coincides with that of Eq. ( 63) on this interval.The constraint of Eq. ( 50) holds on the interval t 1 < t < T 1 , since the control of Eq. (77) on this interval is defined by the relation ξ ξ = −1.For t T 1 , the behavior of the function ξ(t) is governed by Eq. ( 41) with w = 0.

Integration of this equation subject to the initial condition ξ(T
Hence, the control of Eq. ( 77) satisfies the constraint of Eq. ( 50) on the entire nonnegative semiaxis, 0 t < ∞.
Define ẏ2 and y 2 to be functions of the variable ξ on the time interval t 1 < t < T 1 .Since ξ ξ = −1 on this interval, from Eq. (41) w = (αξ 2 − 1)/ξ.Substitute the last expression into Eq.( 68) to obtain Integration of the first equation gives where ẏ2 ( ξ) can be calculated in accordance with the relations of Eqs (70) or (73).Since the relations ξ = α −1/2 and ẏ2 = 0 hold at the time instant T 1 , Eq. (85) implies that One can utilize this relation to determine ξ.By substituting Eq. (86) into Eq.( 85) we finally obtain Substitution of the right-hand side of Eq. ( 87) into the second relation of Eq. ( 84) and integration of the resulting equation yields where y 2 ( ξ) is defined by Eq. (71) or Eq. ( 73).
The value of the criterion J3 for the control of Eq. ( 77) is equal to the value of the variable y 2 at the time instant T 1 , i.e., 4. Calculate the estimate for the absolute error, or the relative error, of the approximation of the minimum rattlespace of the helmet by the quantity Ĵ3 0 .Note that in the expression of Eq. ( 94) for ∆, the parameter T is the dimensionless characteristic time of the control process obtained by dividing the dimensional time by T * defined in Eq. ( 39).
Finally, we can represent the estimate of the minimum rattlespace as The relations of Eq. ( 111) obtained by executing steps 1 to 4 of the above algorithm are represented in terms of dimensionless units.To represent these relations in terms of the original dimensional units, multiply J 0 3 , Ĵ0 3 , and ).The relations of Eq. ( 111) define an interval in which the rattlespace minimum value certainly lies and provide information about the order of magnitude of this quantity.The length of this interval is determined by the quantity ∆ abs .If ∆ abs is sufficiently small as compared with Ĵ0 3 , the value Ĵ0 3 can be regarded as the final quantitative estimate of the minimum rattlespace, which completes the limiting performance analysis.

Numerical results
Two numerical examples are provided to illustrate typical values of the head injury criteria and the absolute minimum of the helmet padding thickness corresponding to the limiting performance characteristics.Consider two types of projectiles hitting a helmeted head -a bullet and an ice hockey puck.Of course, in the case of the bullet, it will be assumed that the helmet shell is strong enough to prevent the bullet from penetrating or substantially deforming the shell.We utilize the translational head injury model of [4] with the following parameters: We assume the mass of the helmet shell, m 4 , to be equal to 0.5 kg.The values D 1 and D 2 measuring the maximum tolerable injury in terms of the criteria J 1 and J 2 , respectively, are defined as For the bullet, define the mass and the impact velocity of the projectile as For the puck, typical values of these quantities are given by The dimensionless parameters α and β defined by Eq. ( 43) and occurring in Eqs (41) and (42) depend only on the parameters of the mechanical model of the response of the head to an impact (head injury model) and the maximum values allowed for the head injury criteria.These parameters are independent of the mass of the helmet, the mass of the projectile, and the impact velocity.For the numerical data of Eqs (112) and (113), α = 5.90 × 10 −2 , β = 5.85 × 10 −2 . (117) The dimensionless parameter λ depends on the characteristics of the head injury model and also on the mass of the helmet shell, but is independent of the mass and velocity of the projectile.In the case under consideration, this parameter is given by λ = 7.16 × 10 3 . (118) The dimensionless impact velocities v 0 of Eq. ( 60) for the bullet and the puck have the values The relative difference between the lower and upper bounds for the criterion J3 calculated in accordance with Sections 6 and 7 is 3 • 10 −5 .Within this error, the value of J3 is defined as According to Eq. ( 97), the relative difference of the optimal value of the criterion J 3 in Problem 2 from the optimal value of the criterion J3 in Problem 1 is estimated to be 40% for the bullet and 70% for the puck.To obtain these estimates, we defined the time T in Eq. ( 51) to be equal to the time of deceleration of the impact velocity v 0 to a complete stop under the action of the (dimensionless) control w = 1.In the dimensionless variables, T = v 0 , in accordance with Eqs (61) and (62).Of course, these estimates indicate a rather high possible error of approximation of the limiting performance characteristics of the helmet by solving Problem 2 instead of solving the original (more complicated) Problem 1.Nevertheless, in the cases under consideration, these estimates are useful, because they restrict the optimal value of the performance index J 3 of Problem 1 to a certain finite interval, while the optimal displacement of the helmet shell relative to the head is small (in comparison with the dimensions of the head and the helmet).These intervals are defined as 1.22mm Jbullet

Conclusions
An analytical technique has been formulated for the limiting performance analysis of a helmet intended to prevent head injuries due to being hit by a non-penetrating projectile.The goal of this analysis is to evaluate the absolute minimum of the peak displacement of the helmet shell relative to the head (the rattlespace of the helmet), provided that quantitative criteria measuring the severity of injuries lie within prescribed admissible limits.The response of the head to an impact was calculated in accordance with the translational head injury model by Rojanavanich and Stalnaker.This model leads to a seventh-order set of ordinary differential equations governing the behavior of the helmet -head system.Among the head injury criteria taken into account were the peak power of elastic forces due to deformation of the skull bone interacting with the helmet liner and the acceleration of the brain.This approach reduces the limiting performance analysis to solving an optimal control problem for a dynamical system governed by three first-order ordinary differential equations, subject to joint constraints on the phase and control variables.This substantially enables the analysis and makes the utilization of numerical methods necessary only for solving relatively simple transcendental equations.This approach was applied to evaluate the minimum rattlespace of helmets preventing the helmeted head from injuries due to being hit by a bullet and an ice hockey puck.The calculations show that a liner of a thickness of several millimeters could be sufficient.

helmet, m 4 v
p , m p m p : mass of projectile

Fig. 2 .
Fig. 2. Model of the head with helmet.
), which completes the proof.