Dynamic characteristics of finite element (FE) model of the human head-neck

2.1. Finite element method and governing equation

In order to determining modal responses, modal analysis using finite element method (FEM) is performed using an implicit finite element code (Abaqus). The problem analyzed is that of looking for eigenvalues and eigenvectors of the rigidity and mass matrices of the model. After system discretization by finite elements, the equation governing the dynamic response is given as follows:

where $\left[M\right]$, $\left[C\right]$ and $\left[K\right]$ are the global mass, damping and stiffness matrices of the model; $\left\{\ddot{u}\right\}$, $\left\{\dot{u}\right\}$, and $\left\{u\right\}$ are the nodal acceleration, velocity and displacement vectors; $\left\{F\right\}$ represents the sum of the external loads applied to the structure.

For undamped free vibration (i.e. $\left\{F\right\}=0\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\left[C\right]=0)$), the solution of the above equation can be written as follows:

where $\left\{U\right\}$ represents the amplitudes of all the masses (mode shapes or eigenvectors) and $\omega (=2\pi f)$ represents each eigenvector’s corresponding eigenfrequency in rad/s, while $f$ represents the natural frequency in hertz. Thus Eq. (1) reduces to:

The above equation is known as eigenvalue problem in matrix algebra and is considered as linear by replacing ${\omega }^2$ by $\lambda$. The system solution, which relies on determining each eigenvector $\left\{X^i\right\}$, with its corresponding eigenvalues ${\omega }_i^2$, is solved by the natural frequency extraction using Lanczos eigensolver [12].

For the case of damped free vibration (i.e. $\left\{F\right\}=0\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\left[C\right]\ne 0$), the eigenvalue equation is quadratic as it contains an additional term of $i\omega \left[C\right]$:

The complex eigenvalue problem is solved by the complex eigenvalue extraction procedure using a projection method [12]. In this present study, four global material damping factors of 0.1, 0.2, 0.3 and 0.4, which are within the range of damping factors of head and neck found in the literature [7, 9, 10, 13], are included in this complex eigenvalue analysis of the multi-degree of freedom problem. In addition, the biomechanical parameters such as skull stresses and intracranial pressure (ICP) are evaluated at each individual mode for both damped and undamped cases. It shall be noted that the stresses evaluated are not due to any external applied force, in fact, they are element stresses ${\left\{\sigma \right\}}_e$ which are computed based on the following equation in FEM [14]:

where ${\left[D\right]}_e$ is the constitutive matrix of the element; ${\left\{\epsilon \right\}}_e$ is the elemental strain.

2.2. Mesh development and interface conditions

A male volunteer with anthropometric characteristics close to the 50th percentile Singapore Chinese male was recruited to develop an extensive image data set. 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. The geometries of the bony structures and soft tissues of the volunteer head region were reconstructed using the CT and MRI scanned images, using the segmentation method developed by Dale et al. [15] and later on by Fischl et al. [16]. With minimum manual edition, we sought to align the MRI to the CT, and registration accuracy was evaluated by performing analysis of the coordinate differences between CT and MR anatomical landmarks along the $x$-, $y$- and $z$-axes. The human brain was segmented into cerebellum, gray and white matters, the entire ventricular system of the brain.

Despite recent advancements in segmentation methods for brain tissue with Magnetic Resonance Images (MRI) [17], there is no automatic segmentation tool available for non-brain tissues such as extracranial tissues like cartilages, fats, and neck muscles. This was owing to the fact that segmentation of these tissue types was often ignored since these tissues were regarded as less important as compared the skull-brain tissue and were not usually considered in the FE head model. Based on the reference to available atlas of head anatomy [18], the geometry of the cartilages, namely the cartilage of septum, the lower and upper lateral cartilages of the human nose, are reconstructed semi-automatically using an adaptive moving mesh technique and shape preserving parameterization.

Multi-tissue mesh generation on medical images is a fundamental step for building a realistic biomechanical model. Mesh elements with large or low dihedral angles are undesirable. In the literatures, there have been studies on multi-tissue meshing based on Delaunay refinement [19, 20].However, elements with small dihedral angles are likely to occur in Delaunay meshes, because elements can be removed only when their radius-edge ratio is large.

Fig. 1Fixed boundary condition is applied at the surface nodes of the base of the neck

Unlike above Delaunay-based methods, Zhang et al. [21] presented a new method to generate a hexahedral and tetrahedral mesh. Firstly, this method identified the interface between different tissues and non-manifold nodes on the boundary. Then, all tissue regions were meshed with conforming boundaries cooperatively. Finally, geometric flow schemes and edge-contraction were used to improve the quality of the hexahedral mesh. In our work, we incorporated mesh quality, fidelity and smoothing into one point based registration framework. The existing hexahedral meshes were optimized on combining Laplacian and optimization-based mesh smoothing, nodal points deletion and insertion as well as local remeshing. The resulting meshed head model is composed of 403, 176 hexahedral elements and 483, 711 nodes (Fig. 1).

The surface nodes at the base of the head-neck model are constrained in all the six degrees of freedom (Fig. 1) while all the interfaces between components are prescribed with tie conditions since interaction properties are not available in modal analyses.

4.1. Comparison of fundamental frequency

In the present work, Meyer et al. [22]’s FE study is chosen for comparison of the modal responses as it is the only study involving the head-neck system which is similar to the present FE model. Our computed fundamental frequency in the free vibration is 35.57 Hz while Meyer et al. [22] reported that the fundamental frequency is 3.01 Hz. This discrepancy in the fundamental frequency arises between the two FE models, probably due to the difference in the way in modeling and material properties, with the former factor being more likely to contribute to the discrepancy. With the prioritized focus on the neck injury, the head and all the cervical vertebrae of Meyer et al. [22]’s FE model were modeled as rigid bodies, with all the individuals’ masses and inertial moment taken into account. Only the intervertebral discs were modeled as deformable bodies. On the contrary, all the components in our FE head-neck model are modeled as deformable bodies with the ability to dissipate external energy through deformation. In addition, unlike the rigid head of Meyer et al. [22]’s FE model, our FE head model contains skull, subarachnoid space (including CSF and membranes), cartilages as well as brain tissues. By having the various distinct components in the head model, especially with the lumped mass of viscous intracranial content, the fundamental frequency of the multi-components system is expected to be lower than that of a one-component system as additional natural frequencies appear on the frequency spectrum of the one-component system. This phenomenon is consistent with observations by Guarino and Elger [23] and Chu et al. [24].

In summary, our computed fundamental frequency of 35.57 Hz is reasonably close to the reported range of 20 to 30 Hz in book reviews [25-27] and handbook [28]. Although the present value is lower than the few hundred hertz of the analytical skull-brain-neck system by Charalambopoulos et al. [2] and higher than the few hertz as predicted by Meyer et al. [22]’s FE study, it is still considered to fall within a reasonable range since the consensus on this dynamic characteristic remains debatable.

Fig. 2Mid-sagittal view of the displacement magnitude (in mm) contour plots of the head-neck model, showing various mode shapes and their corresponding frequencies for the undamped vibration case

Fig. 3Comparison of the mode shapes with Meyer et al. [22]’s FE head-neck model

4.2. Comparison of mode shapes

Most of the previous studies compared their resonant frequencies with those found in the literature whilst modal comparison in terms of mode shapes was often ignored. A limited number of published works [7, 24] on the vibration patterns of the human head can be found in literature and some of them are unsuitable for comparing human head’s mode shapes. Khalil et al. [7]’s simplified representation of mode shapes in terms of stationary nodal lines were only indicated on the top and lateral sides of the skull, due to extremely complex geometry of the skull base and the facial region, whilst Chu et al. [24]’s modal analysis using a 2D FE skull-brain lacked detailed vibration patterns. To the best knowledge of the author, the only work that had recorded the detailed description of mode shapes of human head and neck, was the FE study by Meyer et al. [22]. In order to further validate our newly developed 3D FE model of head-neck, the simulated fundamental frequency is compared with those in the literature while the mode shapes obtained from the simulations are also compared with the mode shapes as predicted by Meyer et al. [22]’s FE head-neck model.

Fig. 4(Still Instance No. 1) Comparison of the mode shapes with Meyer et al. [22]’s FE head-neck model

It is demonstrated in Fig. 3-5 that the first three mode shapes in our model are similar to that of Meyer et al. [22]’s FE model, which correspond to extension-flexion mode, lateral flexion mode as well as axial rotation mode of the head respectively. It is then followed by the fourth mode shape resembling the resonance of the mandible, which were not included in Meyer et al. [22]’s FE model. Our fifth mode shape matches well with the S-shaped anterior-posterior retraction mode found in Meyer et al. [22]’s fourth mode. However, the fifth mode shape of lateral retraction is not captured in the simulations of the present model. Instead, the eighth mode of the current study is considerably similar to the vertical translation mode (sixth mode) of the neck in Meyer et al. [22]’s model. The current model, with detailed facial features, is able to simulate or capture the additional modes of the lateral flexion of the nasal lateral cartilages and the “mastication” mode of the mandible which have never been shown in any vibration analysis. The presence of these additional modes also indicates the importance of detailed modeling in identifying all the resonant modes of the individual components.

Fig. 5(Still Instance No. 2) Comparison of the mode shapes with Meyer et al. [22]’s FE head-neck model