A Study on the Finite Element Model for Head Injury in Facial Collision Accident

: In order to predict and evaluate injury mechanism and biomechanical response of the facial impact on head injury in a crash accident. With the combined modern medical imaging technologies, namely computed tomography (CT) and magnetic resonance imaging (MRI), both geometric and finite element (FE) models for human head-neck with detailed cranio-facial structure were developed. The cadaveric head impact tests were conducted to validate the head-neck finite element model. The intracranial pressure, skull dynamic response and skull-brain relative displacement of the whole head-neck model were compared with experimental data. Nine typical cases of facial traffic accidents were simulated, with the individual stress wave propagation paths to the intracranial contents through the facial and cranial skeleton being discussed thoroughly. Intracranial pressure, von Mises stress and shear stress distribution were achieved. It is proved that facial structure dissipates a large amount of impact energy to protect the brain in its most natural way. The propagation path and distribution of stress wave in the skull and brain determine the mechanism of brain impact injury, which provides a theoretic basis for the diagnosis, treatment and protection of craniocerebral injury caused by facial impact.


Introduction
In traffic accidents, head injury has become the main cause of serious injuries or death. According to WHO statistics, vulnerable road users account for 50% of the deaths in road traffic accidents, and pedestrians account for 22% of the deaths [1]. In China, the incidence of craniocerebral injury with the number of disabilities and deaths is increasing annually. Therefore, it is of great practical significance to study the head injury of pedestrians.
Facial injury and concomitant traumatic brain injury (TBI) have been the focus of numerous investigations over the past few decades. On account of the close anatomical proximity of the facial skeleton and cranium, it is not surprising that patients with facial trauma are at higher risk for suffering brain injuries. Several earlier studies had been conducted in evaluating the incidence of facial injuries and associated injuries [2][3][4]. Prior statistical findings from retrospective clinical cases reported a wide range of incidence rates of brain injuries associated with facial fractures; with the lowest rate as 5.4% [5] whereas some rates as high as 80% [6][7]. With the fleeting progress of computer technology, the finite element model has become an important tool to learn about the biomechanics of craniocerebral injury [8][9][10]. Santos et al. [11] used finite element method to analyze stress distribution of edentulous mandibular median commissure, mandibular foramen and mandibular body under traumatic load. Tuchtan et al. [12] conducted a study regarding direct mandible impact on the level of energy required to create a mandible fracture and on the energy dispersion phenomenon to the skull and the brain. In 2015, Huempfner-Hierl et al. [13] established a finite element model of human head without mandible, then applied this model to simulate three collision tests of the orbital margin, the thyroid gland of the nose and the upper rim of the orbit.
In summary, these previous studies provide effective information about the propagation of stress waves in the skull, but many of these finite element head models are too simplified or ignore the mandible and intracranial tissue, and only aim at a small number of collisions locations, which has certain limitations for the correlation study between facial impact and traumatic brain injury.
In this research, a 50th percentile finite element model of head with detailed craniofacial structure is established. The biofidelity of the model is validated by the corpse experimental data of Nahum et al. [14] and Trosseille et al. [15]. Through the analysis of the collected face collision accidents, nine kinds of common facial collision scenes are simulated to clarify the relationship between facial collision and human craniocerebral injury, which provides the basis for the study of human head injury mechanism and injury evaluation criteria.

Geometric Model
Through tomography and MRI scanning in line with the body features of the 50th percentile adult Chinese male head, a tomography image has been achieved. The resolution/thickness of the computed tomography (CT) and magnetic resonance imaging (MRI) scans were 0.488/1.0 mm and 0.500/4.0 mm, respectively. Then this study imports the image into Mimics for pre-processing, which includes filtering, image segmentation and contour extraction, and then obtains the point cloud data and imports it into Geomagic Studio. The point cloud data is parameterized, and then generates surface, detects the surface, removes the noise points and repeated points on the surface, and then three times Bessel curve is used to fit each surface to get the 3D geometric model. The final three-dimensional solid model of the head consists of skull, facial skeleton, cartilage, teeth, brain, cerebellum, meninges, cerebrospinal fluid, etc., as shown in Fig. 1.

Finite Element Model
Three-dimensional geometric model is imported into Hypermesh for grid division, and then solid three-dimensional grids are generated. It is difficult to use hexahedral elements to model because of the high folds of the meninges in the model. Tetrahedron is adaptable to highly complex surface structure, so tetrahedron is more suitable for the discretization of the mesh. Thus the unit type is set as a linear tetrahedral element in this study. The divided tetrahedral meshes are inserted and deleted through nodes, and optimizes the local grid. Finally, the three-dimensional entity model is generated and imported into Abaqus, the material properties of each tissue are shown in Tab. 1.
In the model, the dura mater is attached to the internal cranium, and the pial membrane is adhered to the brain upper surface. And the CSF (cerebrospinal fluid) layer is between the dura mater and pia mater.
Finally, the finite element model of the head is achieved, in which the number of single element is 1,337,903, the number of nodes are 327,536, the average length of unit is 1.57 mm, and a unit average aspect ratio of 1.61, as shown in Fig. 2.

Material Properties
After meshing the model, it is necessary to assign material parameters corresponding to the actual situation to each part of head finite element model. In all the latest finite element head models, skeletal tissues such as cartilage, skull and facial bone are recognized as linear elastic and isotropic materials. Like other finite element models of the head, the linear viscoelastic material combined with large deformation theory is used to simulate the time-varying relative movement of brain tissue. The material properties used in the model are referenced in the literature [16]. The material properties of these skeletal and brain tissues are summarized in Tab. 1.

Experimental Verification of Nahum Intracranial Pressure
The model is first validated against the ICP-time histories produced by Nahum et al. [14]. In order to test the feasibility of the cadavers for collision damage research, Nahum et al. [14] used a human cadaver without antiseptic treatment as samples. Special experimental tools were used to clamp the head of the cadaver to avoid the rotation acceleration of the head after impact. In this study, a cylindrical rigid body is used to hit the frontal bone of human head finite element model at the speed of 9.94 m/s, the horizontal plane is 45 degrees to the Frankfurt plane.
Nahum et al. [14] only gave the various curves of experimental samples numbered thirty-seventh groups, for other experiments, he did not give a detailed description. In the next 30 years, the experimental data was verified by most of the major human head finite element models. Therefore, this research compares the simulation results with the experimental curves of the thirty-seventh groups of Nahum et al. [14].
The pressure curves extracted from the impact side, impact opposite side, parietal lobe and occipital bone are shown in Fig. 3. It can be seen in the four parts (a-d) of the Fig. 4 that the changing trends of the simulation curve and the experimental curve are basically the same, the curve anastomosis is good on the whole. The peak pressure in four collision regions is about 1 ms ahead of the experimental results of Nahum et al. [14], the maximum peak is slightly different. This is mainly because there are some distinctions between the simulation model and the experimental model in the geometry size and structure, but the peak duration is roughly equal to the experimental results.

Experimental Verification of Trosseille Intracranial Dynamic Response
In 1992, Trosseille et al. [15] reconstructed the impact process between the human head and the car parts in a vehicle accident. The impact positions included the face, forehead and chest, and the intracranial pressure of the corresponding points was measured. This research uses steering wheel, which is modelled as a rigid cylinder with an inertial mass of 23.4 kg, to impact at the nasal of the head finite element model. In this experiment, the bottom surface nodes of C7 vertebrae are set as a fixed boundary condition. The translation acceleration of the Y axis is ignored in the simulation as the impact occurs in the head median sagittal plane. The acceleration of the X axis and Z axis in the head centre and the rotation acceleration of the Y axis are loaded on the rigid skull as the initial boundary conditions of the simulation. The pressure of frontal lobe, temporal lobe and occipital lobe are compared with experiment values, as shown in Fig. 4.
As can be seen from Fig. 4, the simulation curve is anastomotic well to the experimental curve, and the pressure distribution trend in three positions is similar to the experimental results. For the peak pressure, temporal and occipital peak values are more close to the experimental results, the difference is less than 10%, and the frontal pressure value is high. For the curve trend, frontal and occipital are consistent with experimental results and the simulation results of the temporal pressure is negative firstly, and the experimental result is almost 0. For the peak time, the three positions are different, the duration of peak pressure is also slightly different from the experimental results.
Similar deviations have also been observed in literature [17]. It may be due to the large differences between the model and experimental samples. In the experiment, the impactor is used to impact the head, and the acceleration is directly loaded in the simulation to avoid the influence of neck.

Simulation Experiment of Face Collision
As shown in Fig. 5, nine impact directions, which are believed to be commonly encountered in the real-world impact scenario, are simulated in the present study [18], it contributes to a deeper understanding of the mechanism of maxillofacial and craniocerebral injury. Similar to the impactor used by Cormier et al. [19] in the cadaver experiment, a cylindrical rigid body is used to impact the human head finite element model at the speed of 2.5 m/s. However, unlike the rigid constraints generated by mounting brackets in the Cormier et al.'s experiment, a fixed boundary condition is imposed on the nodes of the subvertebral surface of C7 in the experiment, the surface nodes at the base of the model are constrained in all the six degrees of freedom.

Tolerance Limit of Craniocerebral Injury
In order to analyze potential fractures and distinguish the region of the facial and cranial bones fracture, von Mises criterion is used in the current work because bone exhibits both ductile and brittle material properties. Keyak and Rossi [20] found that the "von Mises criterion is the most accurate for predicting fracture location, when isotropic material properties are used". Therefore, an average value of von Mises stresses of 75 MPa is used as the yielding limit of the skull bone according to the reported range of McElhaney et al. [21]. We compare the biomechanical parameters after trauma with the brain tissue tolerance limit in the literature, which can be used to evaluate traumatic brain injury.

Simulation Results Discussion
Through simulation, the propagation process of head stress wave in nine cases is obtained. In Scene 1, the stress first concentrates on the impact position, then spreads to the bilateral frontal lobes of the maxilla and the posterior end of the nasal bone. As the base of the model is fixed, the bending of the head and neck contributes to local stress at the lower neck, it also causes the relative movement of the skull and brain, which results in stress concentration in frontal lobe and temporal lobe. Recoverable injury and axonal injury may occur at the inferior frontal lobe. At t = 2.0 ms, there is a stress of 65 Mpa at the anterior nasal spines, nasal septum cartilage and lateral cartilage junction. According to the von Mises yield criterion, fractures occur easily around the nasal bone.
In Scene 2, starting from the impact position, stress waves propagate in two ways: One is backward through the frontal and ethmoid plate, and then reaches the inferior frontal lobe of the brain. The other propagates downward through the perpendicular plate and nasal septum. Afterwards, the stress wave reaches the palatal floor of maxilla. On the basis of von Mises yield criterion, there is a risk of fractures at the nasal bone and bilateral maxilla frontal processes. In addition, the fracture regions also include orbital floor of maxilla, bilateral faces nasalis maxillae, cribriform plate of ethmoid bone and inferior nasal septum. In addition, mild TBI (traumatic brain injury) will appear at the temporal pole and inferior frontal lobe, and severe injury or mild DAI (diffuse axonal injury) will occur at the lateral occipito-temporal lobe and cerebellum.
In Scene 3, the stress wave first passes through the upper lateral nasal cartilage and the posterior nasal bone, and then spreads to inferior frontal lobe from the lower lateral cartilage, then through the nasal septal cartilage and the sieve disk of plow bone, and then spreads to the sphenoid from the frontal bone and ethmoidal cribriform plate. During the impact, the nasal septum cartilage and the lateral nasal cartilage are collapsed by the impactor, which causes that stress is raised and concentrates around the tip of nose and the nasal septal cartilage. According to the von Mises yield criterion, there is a risk of fracture at the anterior nasal bone and lateral nasal cartilage. Similar to Scene 1, the stress is concentrated on the frontal lobe and temporal lobe, thus minor recoverable injury and mild injury may occur at the anteriorinferior frontal lobe and inferior frontal lobe.
In Scene 4, starting from the dense incisors, a part of the stress wave passes via the alveolar processes of mandible, vomer and maxillary palatal process before reaching the sphenoid region. Another part of the stress wave reaches the bilateral maxilla frontal process through the alveolar process of mandible along the pyriform aperture, it results in local stress concentration in the temporal lobe, foramen magnum and the occipital lobe. After about 4.4 ms, a part of the stress spreads to the exterior upper and bilateral zygomatic bone, this leads to the peak von Mises stress at the maxilla, the value of which is 203.66 MPa. Therefore, Le Fort I fracture is easy to occur at the alveolar process of mandible, at the same time, high stress is easily produced at the foramen magnum and bent lower neck.
Scene 5 is a vertical impact on the base of mandible, and Scene 9 is an oblique impact on the base of the mandible. In the process of impact, the stress wave passes through the bilateral mandibular ramus, mandibular arch and condyles, and then spreads to the zygoma along the zygomatic arch. The head is elastic during the impact, which leads to the relative rotation of the skull and brain, it makes the stress concentrate at the frontal lobe. In Scene 5, at t = 3.0 ms, peak stress concentrations of 209.55 MPa appear at the mandibular ramus and mandibular body, which is likely to lead to Le Fort III fractures in the ethmoidal cribriform plate. Moreover, Mild injury will appear at the frontal lobe, local contusion may appear at the anterior upper brainstem. At the anterior-inferior frontal lobe, it may produce mild DAI. In Scene 9, peak stress concentrations of 267.52 MPa appears at t = 3.6 ms, mild brain damage is likely to occur at the anterior inferior frontal lobe and temporal lobe, and concussion may produce at the frontaltemporal lobe.
In Scene 6, the stress quickly propagates to the condyle of mandible and mandible body. Then, from the right condyle of mandible, to the temporal bone and zygomatic bone, to the maxilla nasal surface. The stress wave propagates from the mandibular body to the contralateral ramus of mandible. The peak stress appears at the mandibular condyle and mandibular incisor, the value of which is 328.89 MPa. The stress then spreads to the sphenoid bone, foramen magnum and vertebrae. After about 1.8 ms, the lateral bending of the head causes the stress concentrate at the right frontal lobe, right parietal lobe and right temporal lobe, mild brain damage is likely to occur in these regions.
In Scene 7 and Scene 8, the stress wave propagates mainly through the maxillary palatine process, the maxilla and the sphenoidal orbit surfaces to the sphenoid and middle surface structures. In Scene 7, peak stress of 244.69 MPa appears at the zygomatic region after about 1 ms, and there is local stress concentration on the maxilla, zygomatic orbit surfaces and frontal process of the zygoma. In Scene 8, peak stress of 537.94 MPa appears at the zygoma after about 1.2 ms. In the two cases, there are a large number of stress concentration at the zygomatic and maxillary orbit surfaces, particularly in Scene 8, because the collision position is proximal to the orbital area. Although the stress concentrations in Scene 7 and Scene 8 are lower than that in Scene 6, there is still a higher stress concentration at the anteriorinferior frontal lobe and right temporal lobe, which is due to the large relative displacement between the frontal lobe and frontal area of skull.

Intracranial Biomechanical Parameters
In nine impact scenarios, when the impact force reaches the peak value, the maximum values of the intracranial biomechanical parameters of these locations are achieved, as shown in Tab. 2. Meanwhile, the contour plots of median sagittal plane as shown in Fig. 6 to Fig. 11 are obtained.  Fig. 6 indicates that the maximum ICP of frontal collision and mandibular base impact occurs at the inferior frontal lobe, which is close to the mid-sagittal plane, while the maximum ICP appears near the impact side in lateral impact. Among all the scenes, the frontal collisions at the nasal bone (Scene 1 and Scene 2) are the most serious, and their peak pressure are close to the brain tolerance limit of 235 KPa, which proposed by Ward et al. [22]. The second are the oblique impacts at the zygomatic maxillary area (Scene 7 and Scene 8), the peak pressure appears at the frontal-temporal lobe. In the various impacts of the mandible (Scene 5, Scene 6 and Scene 9), the ICP value has a moderate increase. In frontal impacts (Scene 3 and Scene 4), the intracranial pressure in the lateral cartilage and teeth are relatively low.

Von Mises Stress
As shown in Fig. 7, the maximum von Mises stress appears at the frontal-temporal lobe and frontal lobe. While in Scene 6, the peak von Mises stress is situated at the posterior cerebellum. In all cases, Scene 1 and Scene 2 are the nasal bone impacts, which are the most serious, came along by lateral impacts (Scene 7 and Scene 8), which has the risk of mild TBI (traumatic brain injury). When the maximum von Mises stress is produced by the frontal collision on teeth (Scene 4) and the base impact on the mandible (Scene 9), it may develop into moderate neurological lesions [23]. The injury caused by frontal collision t on the lateral nasal cartilage (Scene 3) is the lowest, with the maximum stress of 5.72 KPa, there is a risk of producing local contusion [15].

Shear Stress
From Fig. 8 to Fig. 10, we can see that in all cases, the critical shear stress positions are mainly in three directions (i.e., fore-and-aft, left-right, up-down). In the lateral collisions on the zygomatic maxillary bone (Scene 7 and Scene 8), the maximum shear stress is generated in the left-right direction, whereas the base collisions at the mandible lead to the maximum shear deformation of the brain tissue in the up-down direction. According to the tolerance limit of shear stress [24], except Scene 3, there is a possibility of mild TBI or mild DAI in most scenes.

Strain
You can see from Fig. 11, the critical strains appear at the frontal-temporal lobe region nearby the anterior-superior brainstem in most scenes, except Scene 3, Scene 4 and Scene 6. In Scene 2, the peak strain is about 0.13, it is the conservative threshold for 25% probability of developing mild TBI. In other scenes, the peak strains are less than 0.10, and the injuries is reversible according to the strain injury limit of Galbraith et al. [25]. The strains of the zygomaticomaxillary bone (Scene 7 and Scene 8) are more than 0.09, and the strain is the lowest in the frontal impact on nasal tip (Scene 3), and the value is 0.023.

Discussion
In frontal impact on nasal bone, the general trend of maximal intracranial parameters reveals that the distance from the impact position to the brain determines the severity of brain injury, as well as in the lateral impacts at the zygomatic maxillary region. It verifies the conclusion of Lee: "The risk of TBI caused by facial fractures adjacent to the brain is very high." In facial trauma, the fractures of nasal bone are the most common fractures, which are normally accompanied with the maxilla orbital floor fracture as well as the fractures of the ethmoidal cribriform plate and bilateral maxillary nasal surface. Therefore, this kind of injury is called the nasal-orbit-ethmoidal complex fractures.
In Scenes 1, 2 and 4, they are likely to result in maxillary fractures, which are normally accompanied by nasoethmoidal fractures, it will further lead to the occurrence of mild DAI at the medulla oblongata. Meanwhile, moderate neurological lesions will appear at the frontal-temporal lobe. Just like Chang's hypothesis: "A structure is formed by maxilla and other facial bones in midface, which can absorb considerable energy to protect the brain from impact". When the face is impacted, a portion of collision energy is assimilated by the midfacial structure, and the other energy is transferred to the intracranium.
Nasal bone fractures are found in Scenes 1, 2 and 3, but in Scene 3, the risk of skull fractures and craniocerebral injuries is the lowest, because nasal lateral cartilage has less rigid structure, so it can absorb impact energy by means of deformation and fractures. In accordance with the clinical review of Weller et al. [26], this research also illustrates that the nasal cartilages can behave as a buffer region to decrease the impact force that will transmit to the skull and minimize the risk of sustaining severe traumatic brain injuries.
Mandibular fractures may exist in Scenes 5, 6 and 9, the associated facial fractures also include rami mandibular, coronoid process of mandible and bilateral facies nasalis maxillae. Thus, the incidental brain traumas range from mild recoverable injury to mild DAI in the brainstem, cerebellum as well as inferior frontal lobe. Mandibular fracture is usually accompanied by midfacial fracture, though the mandible is not close to the brain. Therefore, it still causes mild to moderate traumatic brain injury, especially in inferior frontal lobe.
The propagation of stress waves in cranial bones decides the onset and progression of various facial fracture patterns. Moreover, it has important clinical significance for the mechanism research of intracranial injury. This study shows that a large number of impact energy is transmitted to the brain through the facial structure, that is, the collisions do not affect the brain directly. This may lead to speculation that the impact is passed through the midfacial skeleton structure to the brain. Nevertheless, through observing the propagation of stress, it can be found that stress waves are mainly concentrated in the prefrontal lobe before stress wave passes through the frontal sinuses. It indicates that the relative motion of skull and brain leads to the elevation of frontal lobe stress, in which it strikes the protuberance of the skull internal surface, rather than the stress wave propagation. Nevertheless, localized stress concentration, which causes by the propagation of stress waves, is still a relevant factor leading to injury in other locations.
The results show that in most scenes, the von Mises stress at the VE region is higher, it is consistent with clinical diagnosis. According to the study of Aykan et al. [27], the plowbone ethmoid located at the posterior part of nasal septum cartilage is the most common fracture area in nasal trauma. This suggests that the skeletal structure in the middle of the face may be the focus of stress transmission, where the collision energy is further dissipated to "buffer zone" of the ethmoid sinus, reducing the risk of craniocerebral injury.

Conclusions
By establishing and verifying a fifty percentile head finite element model with detailed craniofacial characteristics, and simulating nine kinds of common facial impact scenes in vehicle traffic accidents, we draw the following conclusions: (1) When the face is impacted, the facial structure can protect the brain by assimilating a lot of impact energy through deformation and fracture, thereby decreasing the risk of craniocerebral injury.
(2) According to the shear stress tolerance, facial impact may lead to TBI or DAI. Even patients with minor facial injuries may lead to craniocerebral injury.
(3) The propagation path and distribution rule of stress waves in the skull and brain determine the mechanism of craniocerebral impact injury, which provide a theoretical basis for the diagnosis, treatment and protection of craniocerebral injury caused by facial impact.