Abstract
Traumatic brain injury (TBI) caused mainly by external head impacts poses significant health concerns globally. Understanding the mechanics behind TBI during head impacts is crucial for developing effective preventive and therapeutic strategies. In this study, the fluid–structure interaction within the system comprising the skull, brain, and the cerebrospinal fluid (CSF)-filled gap is investigated using finite element methods (FEM). Unlike most studies that model CSF as a solid, this research models CSF as a fluid, focusing on its fluid dynamics and their impact on brain tissue response to external head impacts. Additionally, this study is the first to model CSF as a non-Newtonian fluid, exploring its influence on injury metrics compared to a Newtonian CSF model. The results demonstrate significant pressure build-up and shear rate variations within the CSF due to impact. The model shows that maximum strain values are concentrated in the central regions of the brain tissue rather than at its interface with the CSF. Comparative analysis of the first and third principal strains shows that the tissue experiences twice as much compressive strain compared to tensile strain. Further, the comparison between Newtonian and non-Newtonian CSF models shows that the non-Newtonian model results in lower shear rates. This leads to a decrease in tissue strain, with a 4.3% reduction in the first principal strain and a 6.7% reduction in the third principal strain for the non-Newtonian CSF model. These findings underscore the importance of accurately modeling CSF properties to better understand TBI mechanisms.
Introduction
Traumatic brain injury (TBI), causing a substantial number of deaths and disabilities annually, is a significant health concern [1]. External head impacts are identified as the leading cause of TBI, underscoring the importance of studying the mechanical impacts on the head [2]. Cerebrospinal fluid (CSF) as an integral part of the glymphatic system plays various roles, including waste clearance from the central nervous system [3,4]. This fluid, surrounding the brain tissue inside the skull, plays an important role during external head impacts. The system including the brain, CSF, and skull, represents a complex fluid–structure interaction [5]. A better understanding of CSF dynamics and its interaction with brain tissue during an impact is essential for understanding the mechanisms behind TBI and developing more effective protective gear and treatments.
Numerical simulations and the finite element method (FEM) are powerful tools for studying brain tissue under head impacts [6]. To accurately assess injury and damage using numerical models, it is important to select appropriate injury metrics that quantitatively represent the location and magnitude of damage. Research has shown a strong correlation between tissue damage and the strain and stress applied to the tissue [7,8]. Zhan et al. [9] highlighted that strain and stress are reliable indicators of injury in impact-induced studies. Bar-Kochba et al. conducted a study using a 3D in vitro compression model to explore the dependency of neuronal injury on strain magnitude and rate. They reported that while the time of neuronal death is affected by strain magnitude, the pathomorphology and extent of population injury are influenced by strain rate [10]. These findings underscore the importance of stress and strain from FEM simulations as injury metrics.
Most studies to date have modeled CSF as a viscoelastic solid, overlooking the significant effect that its fluid dynamics have on brain tissue deformation and damage [11,12]. Yang et al. conducted a study on different modeling approaches and their influence on injury metrics. They found that material properties and brain skull interface significantly affect the stress and strain in brain tissue, highlighting the importance of selecting appropriate material models for simulations [13]. There are, however, very few studies that account for the fluidic behavior of CSF. Lang and Wu studied CSF fluidic behavior under external head impacts by simplifying the geometry to two concentric cylinders with a gap filled with fluid. They provided analytical solutions for pressure and velocity variations under loading applied to the outer cylinder, showing that the fluidic behavior of CSF significantly influences brain response and reporting the importance of the subarachnoid space gap size [14]. Experimental studies have shown the importance of CSF pressure variation during external head impacts as it may result in microbubble generation and further tissue damage due to bubble burst [15–17]. To the best of our knowledge, no studies have comprehensively addressed the non-Newtonian behavior of CSF in the context of head impacts. Modeling CSF as a fluid, especially as a non-Newtonian fluid, is expected to provide a more accurate representation of its effects on brain tissue during external impacts.
This study aims to fill this gap by using numerical simulations to analyze the interaction between CSF and brain tissue for both Newtonian and non-Newtonian CSF. In this study, we focus on CSF and fluid–structure interaction, reporting injury metrics and how Newtonian and non-Newtonian models of CSF are going to affect strain and stress as injury metrics, hence brain response to an external impact. Shear rate and pressure variation of CSF as a fluid are studied and reported alongside the response of the brain tissue including the stress and strain distribution. Additionally, the non-Newtonian model of CSF is studied and its effect on brain tissue response is reported. This comprehensive study of CSF fluidic behavior for Newtonian and non-Newtonian models provides a better understanding of the mechanics of brain injury during impacts and contributes to practical implications in protective, preventive, and treatment measures.
Methods
In this study, we aim to investigate the fluid dynamics of CSF and its effect on brain tissue during external head impacts, using fluid–structure interaction simulations performed with comsol multiphysics, and compare the influence of Newtonian and non-Newtonian CSF models on tissue. Figure 1(a) illustrates the geometric model that consists of two concentric spheres representing the brain (inner sphere, radius ) and the skull (outer sphere, radius ), with the gap between them filled with CSF. The computational domain, as shown in Fig. 1(b), was discretized using a finite element mesh. The mesh was refined in regions of high deformation and fluid–structure interaction, particularly within the cerebrospinal fluid gap, to accurately capture the dynamic response during head impacts. Additionally, a moving and deforming mesh approach was used to account for the motion of the skull and the resulting changes in the gap size.
![Model setup and validation for CSF dynamics during head impact. (a) Schematic representation of the geometric model used in the simulation. (b) Mesh of the simulation domain showing initial mesh configuration and deformed mesh at a later time, highlighting mesh adaptation during simulation. (c) Time history of the applied acceleration loading. (d) Experimental data of CSF non-Newtonian behavior and the proposed Casson–Papanastasiou model. (e) Validation of ICP dynamics against experimental data from cadaver studies by Nahum et al. [18].](https://asmedc.silverchair-cdn.com/asmedc/content_public/journal/openengineering/4/10.1115_1.4067643/1/m_aoje_4_041003_f001.png?Expires=1744125227&Signature=styQTDRt1kCjySO52xJLoTPgKVR4ekKuKAQ3b9CoAFgTetQlkqJXkHo~xVjLckCxgFXH3-oNu1xvsAA-xhBrgJiHF6sWbrMIF4GPBDV~zV4RYtqS1HLNYV~0g1lqwuPbKYM5HBXNP6oteTC8RHaRRLu9vvwjSJbIJcc3fUCRI0UoBn0SBS6Z4ZrIvEmRPPqDcxxFx8MK33jkFRip2J3x-GRTLnfpt0g64RAffODB5sjzKJZunTqqPitrrsURqC6jF-397gwd7o-HwOjgyQEJAdlHPnswPOzUGcLCIHyANR3klDq8bKUBakHBJfcHhbkp5U4KJkvrlK-5PXtdYGoAOw__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Model setup and validation for CSF dynamics during head impact. (a) Schematic representation of the geometric model used in the simulation. (b) Mesh of the simulation domain showing initial mesh configuration and deformed mesh at a later time, highlighting mesh adaptation during simulation. (c) Time history of the applied acceleration loading. (d) Experimental data of CSF non-Newtonian behavior and the proposed Casson–Papanastasiou model. (e) Validation of ICP dynamics against experimental data from cadaver studies by Nahum et al. [18].
![Model setup and validation for CSF dynamics during head impact. (a) Schematic representation of the geometric model used in the simulation. (b) Mesh of the simulation domain showing initial mesh configuration and deformed mesh at a later time, highlighting mesh adaptation during simulation. (c) Time history of the applied acceleration loading. (d) Experimental data of CSF non-Newtonian behavior and the proposed Casson–Papanastasiou model. (e) Validation of ICP dynamics against experimental data from cadaver studies by Nahum et al. [18].](https://asmedc.silverchair-cdn.com/asmedc/content_public/journal/openengineering/4/10.1115_1.4067643/1/m_aoje_4_041003_f001.png?Expires=1744125227&Signature=styQTDRt1kCjySO52xJLoTPgKVR4ekKuKAQ3b9CoAFgTetQlkqJXkHo~xVjLckCxgFXH3-oNu1xvsAA-xhBrgJiHF6sWbrMIF4GPBDV~zV4RYtqS1HLNYV~0g1lqwuPbKYM5HBXNP6oteTC8RHaRRLu9vvwjSJbIJcc3fUCRI0UoBn0SBS6Z4ZrIvEmRPPqDcxxFx8MK33jkFRip2J3x-GRTLnfpt0g64RAffODB5sjzKJZunTqqPitrrsURqC6jF-397gwd7o-HwOjgyQEJAdlHPnswPOzUGcLCIHyANR3klDq8bKUBakHBJfcHhbkp5U4KJkvrlK-5PXtdYGoAOw__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Model setup and validation for CSF dynamics during head impact. (a) Schematic representation of the geometric model used in the simulation. (b) Mesh of the simulation domain showing initial mesh configuration and deformed mesh at a later time, highlighting mesh adaptation during simulation. (c) Time history of the applied acceleration loading. (d) Experimental data of CSF non-Newtonian behavior and the proposed Casson–Papanastasiou model. (e) Validation of ICP dynamics against experimental data from cadaver studies by Nahum et al. [18].
The fluid and solid domains are coupled, ensuring the continuity of stress and velocity at the interface of the domains.
The simulation results were validated using experimental data from the study by Nahum et al., which investigated intracranial pressure (ICP) dynamics during head impacts using cadaver experiments [18]. Given the simplified spherical geometry used in the simulations, ICP was monitored at two points: one at the frontal region and the other at the occipital region. The comparison between the normalized ICP from the simulation and the experimental data is shown in Fig. 1(e). The results demonstrate a good agreement between the simulated and experimental ICP variation over time.
Results and Discussion
The behavior of CSF and brain tissue during an external impact is reported and discussed for both Newtonian and non-Newtonian CSF models. First, the Newtonian CSF case is investigated, focusing on the fluid dynamics of CSF, including its time-dependent pressure and shear rate, and the tissue response to head impacts, such as stress and strain variations. During an external head impact, the variation in CSF pressure is important as it influences the tissue response to the impact. The pressure response of CSF to the external impact is shown in Fig. 2. The pressure contours in Fig. 2(a) illustrate that, initially, the CSF pressure is almost uniform everywhere. As the force magnitude increases, there is a pressure build-up in the coup region and a pressure drop in the counter-coup region inside the CSF gap. At the impact peak, the pressure difference between the coup and counter-coup regions reaches its maximum. After this point, the pressure difference starts to decrease, and eventually, at the end of the impact duration, the pressure returns to being nearly uniform again. This shows the instantaneous response of pressure to the impact. Figure 2(c) shows pressure variation at the interface of CSF and brain, along the arc indicated in Fig. 2(b). Initially, pressure along the arc is uniform, however, as the impact reaches its maximum, pressure along the arc changes with its maximum at the coup and minimum at the counter-coup part. By the end of the impact duration, pressure along the arc goes back to uniform. Figure 2(d) illustrates pressure evolution by time at the selected points, indicating instantaneous pressure variation at the impact peak. CSF shear rate at different time steps is shown in Fig. 3(a). Initially, the shear rate is uniformly low throughout the CSF. At later time steps, the shear rate increases as the impact magnitude increases. However, unlike pressure, the shear rate reaches its maximum with a delay relative to the impact peak. This figure also shows how the location of the maximum shear rate is noticeably affected by the changes in gap size over time. Figure 3(b) presents the average shear rate in the gap, illustrating the increase in shear rate caused by the impact. Figure 3(c) shows the shear rate along the CSF and tissue interface. In contrast to pressure, the maximum shear rate occurs in the central parts of the arc caused by the impact. Figure 3(c) shows the shear rate along the CSF and tissue interface. In contrast to pressure, the maximum shear rate occurs in the central parts of the arc rather than at the ends. von Mises stress distribution inside the brain is shown in Fig. 4. At the initial time steps, stress inside the tissue is uniform. The stress increases by approaching the impact peak. The location of maximum stress is more central rather than peripheral, and during the impact its magnitude increases. Maximum stress moves toward the coup region, but throughout the impact, it keeps being in the central regions rather than the interface of the brain and CSF.

Pressure distribution and dynamics in CSF during head impact. (a) Pressure contours in the CSF at different normalized times. (c) Pressure distribution along the arc length of the CSF at different normalized times. (d) Time history of pressure at the four specific points indicated in figure (b).

Shear rate dynamics in CSF during head impact. (a) Shear rate distribution within the CSF at different normalized times, the applied load direction is indicated. Insets provide magnified views of specific regions to highlight detailed shear rate variations and changes in gap size. (b) Average shear rate over t*, the graph shows a rapid increase in shear rate as the impact force is applied. (c) Shear rate along the arc length of the CSF at various normalized times.

Shear rate dynamics in CSF during head impact. (a) Shear rate distribution within the CSF at different normalized times, the applied load direction is indicated. Insets provide magnified views of specific regions to highlight detailed shear rate variations and changes in gap size. (b) Average shear rate over t*, the graph shows a rapid increase in shear rate as the impact force is applied. (c) Shear rate along the arc length of the CSF at various normalized times.

Time-dependent stress distribution in brain tissue under applied load. The dynamic response of brain tissue to an applied load over time, highlighting the development and propagation of stress through the tissue.
Principal strain as a significant injury metric is extracted for different time steps. Following a trend close to von Mises stress, it is maximum in the central regions of the brain. Strain contours are shown in Fig. 5(a). It starts initially at zero and increases during the impact. The maximum strain value is not at the interface of the brain and tissue. The location and magnitude of maximum strain changes throughout the impact. Figure 5(b) shows principal strains along the center line of the brain. The first principal strain, being positive, represents the maximum tensile strain of the tissue, whereas the third principal strain, being negative, represents the maximum compressive strain of the tissue. The third principal strain, being twice the first principal strain, represents that compressive strain is higher than tensile strain. It also shows how the maximum strain remains at the central parts of the tissue, and it moves toward coup regions. Figure 5(c) shows a snapshot of brain tissue and CSF, illustrating the CSF shear rate and maximum strain inside the tissue. Shear rate along the interface of CSF and brain tissue correlates to the region of maximum strain inside the tissue. The non-Newtonian CSF model is investigated and its dynamics along with tissue response are reported in this section. The results for Newtonian and non-Newtonian CSF are studied and compared. For the sake of comparison, since the viscosity of non-Newtonian CSF varies between the viscosity of water and double the viscosity of water, an additional model with a fluid viscosity of two times that of water is included.

Strain distribution in brain tissue under applied load. (a) First principal strain inside tissue for selected time steps. (b) First and third maximum principal strain inside the brain tissue along the midline for Newtonian CSF. (c) Brain first principal strain and CSF shear rate correlation. Maximum first principal strain happens at locations with maximum CSF shear rate at the fluid–structure interface.

Strain distribution in brain tissue under applied load. (a) First principal strain inside tissue for selected time steps. (b) First and third maximum principal strain inside the brain tissue along the midline for Newtonian CSF. (c) Brain first principal strain and CSF shear rate correlation. Maximum first principal strain happens at locations with maximum CSF shear rate at the fluid–structure interface.
Figure 6(a) compares the time-dependent maximum shear rate in the CSF for three different cases: Newtonian CSF with the viscosity of water, non-Newtonian CSF, and Newtonian CSF with twice the viscosity of water. At the peak of the impact, all three cases exhibit an initial peak in shear rate. The Newtonian CSF with lower viscosity shows the highest initial peak, followed closely by the non-Newtonian CSF and then the Newtonian CSF with the higher viscosity. The Newtonian CSF with the lowest viscosity consistently shows a higher shear rate compared to the other two cases. The non-Newtonian CSF shows a slightly lower shear rate but follows a similar trend to the Newtonian CSF. The Newtonian CSF with higher viscosity has the lowest shear rate. Over time, the maximum shear rate in all cases decreases as the impact force decreases. However, the Newtonian CSF with lower viscosity maintains higher shear rates for a longer duration compared to the other two models. The comparison of shear rates among the three CSF models shows noticeable differences in their response to head impact, due to their distinct fluid properties. The Newtonian CSF with the smallest viscosity has the highest shear rate, suggesting that lower viscosity results in a higher shear rate, potentially leading to more significant strain and stress on brain tissue. The non-Newtonian CSF shows a moderated shear rate response. Understanding these differences is crucial for accurately modeling brain biomechanics.

Comparison of shear rate, pressure, and strain distribution for Newtonian and non-Newtonian CSF models. (a) Maximum shear rate over time for Newtonian CSF with water viscosity, non-Newtonian CSF, and Newtonian CSF with double the water viscosity. (b) Maximum pressure over time for Newtonian and non-Newtonian CSF models. (c) Normalized strain difference between Newtonian and non-Newtonian CSF models at normalized times. The contour plots illustrate the spatial distribution of strain differences.

Comparison of shear rate, pressure, and strain distribution for Newtonian and non-Newtonian CSF models. (a) Maximum shear rate over time for Newtonian CSF with water viscosity, non-Newtonian CSF, and Newtonian CSF with double the water viscosity. (b) Maximum pressure over time for Newtonian and non-Newtonian CSF models. (c) Normalized strain difference between Newtonian and non-Newtonian CSF models at normalized times. The contour plots illustrate the spatial distribution of strain differences.
Time-dependent maximum pressure within the CSF for Newtonian and non-Newtonian CSF models is shown in Fig. 6(b). Initially, the pressure in both models is zero. As the impact force is applied, both models show a rapid increase in pressure, reaching the maximum pressure at the peak of impact. The peak pressure value for both models is the same, . Following the impact peak, the pressure in both models decreases rapidly, returning to near-zero levels by . Similar peak pressure responses suggest that the non-Newtonian properties of CSF do not significantly alter the maximum pressure values during the impact.
Figure 6(c) shows the normalized strain difference between the Newtonian and non-Newtonian CSF models. The contour plots illustrate variation of maximum principal strain at different locations of the brain, providing a detailed comparison of the influence of the fluid properties on the strain within the brain tissue during an external head impact. At the initial time steps, there are minor differences between the Newtonian and non-Newtonian models. At the peak of impact, strain in the counter-coup region is around more for the non-Newtonian CSF model. At corresponding to the end of impact duration, both models result in almost similar strain in the tissue, indicating no noticeable variations. However, postimpact at , the differences are most significant. Newtonian model shows a variation from the non-Newtonian model at the interface of CSF and brain tissue. The Newtonian CSF model shows higher strain at the top region, whereas the non-Newtonian model has a more evenly distributed strain.
In Fig. 7, the first and third principal strain along the midline of the brain tissue for both Newtonian and non-Newtonian CSF models at various time points during an external head impact is presented. At the initial time point, both models show negligible strain along the midline, indicating no deformation before the impact. At , the strain begins to develop along the midline, with both models showing a similar pattern, although the strain magnitude is higher in the Newtonian model. At , the strain values reach their maximum which corresponds to the second peak of maximum shear rate. At and , the differences between the two models become more noticeable, with the Newtonian CSF model exhibiting higher first principal strain and higher third principal strain compared to the non-Newtonian model. At later times, the strain values gradually decrease for both models, but the Newtonian model continues to show higher strain values, although the difference between the models becomes less distinct over time.

Comparison of principal strains along the midline for Newtonian and non-Newtonian CSF models. First and third principal strains along the midline at various time points illustrated in plots (a) and (b), respectively.
The results highlight the critical role of CSF fluid dynamics during an impact. There is a sudden pressure variation within the CSF, with pressure building up at the coup region and dropping at the counter-coup region. The shear rate of CSF also varies during the impact, influenced by changes in the gap size. Maximum shear rate is observed along the CSF–brain interface and is time-dependent, showing an inverse relationship with viscosity, as viscosity decreases, shear rate increases. Stress and strain within the brain tissue are also reported, showing that strain increases with the impact and peaks in the central regions rather than the periphery. The non-Newtonian CSF model shows noticeable differences in shear rate and tissue strain compared to the Newtonian model. Specifically, the lower viscosity CSF model results in a higher shear rate and consequently higher strain values within the tissue. Therefore, the non-Newtonian CSF model exhibits lower maximum principal strain compared to the Newtonian CSF model, underscoring the importance of cerebrospinal fluid behavior and properties in brain biomechanics during impacts.
Conclusion
This study provides a comprehensive analysis of the fluid–structure interaction between the brain, CSF, and skull during external head impacts using the finite element method. By modeling CSF as a fluid, we focused on its fluid dynamics and their impact on brain tissue response. Our findings show the critical role of CSF in brain strain during external head impacts. The results show significant pressure variations and shear rates within the CSF, with the pressure build-up and drop of CSF during the impact. The maximum principal strains are concentrated in the central regions of the brain tissue rather than at the CSF–brain interface. Furthermore, this study investigates the non-Newtonian CSF behavior, which shows noticeable differences compared to the Newtonian model. The non-Newtonian CSF model results in lower shear rates and principal strains within the brain tissue. Specifically, the first principal strain decreased by , and the third principal strain decreased by , demonstrating the potential protective effect of non-Newtonian fluid behavior during head impacts. These insights emphasize the necessity of accurate CSF modeling to understand the biomechanical behavior of the brain during external impacts. The observed differences between Newtonian and non-Newtonian CSF models suggest that incorporating CSF fluid dynamics into FEM simulations can enhance the predictive capability of brain injury models. This could be applied to the design of more effective protective gear and therapeutic strategies, contributing to better outcomes in the prevention and treatment of traumatic brain injuries.
Acknowledgment
This work was partially supported by the UCLA Brain Injury Research Center (BIRC). The help of John Hollister from UCLA in providing the data on CSF rheology is gratefully acknowledged.
Conflict of Interest
There are no conflicts of interest. This article does not include research in which human participants were involved. Informed consent is not applicable. This article does not include any research in which animal participants were involved.
Data Availability Statement
The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.