ARL-CR-0792 ● JAN 2016

US Army Research Laboratory

Quantifying Discretization Effects on Brain 

Trauma Simulations

prepared by Sophia M Haire

American Society for Engineering Education

1818 N St NW, #600

Washington, DC 20036

under contract W911NF-10-2-0076 

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the 

Army position unless so designated by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official 

endorsement or approval of the use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
ARL-CR-0792● JAN 2016

US Army Research Laboratory

Quantifying Discretization Effects on Brain 

Trauma Simulations

prepared by Sophia M Haire

American Society for Engineering Education

1818 N St NW, #600

Washington, DC 20036

under contract W911NF-10-2-0076

Approved for public release; distribution is unlimited.
ii

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing the 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid 

OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

January 2016

2. REPORT TYPE

Final

3. DATES COVERED (From - To)

June 2015–August 2015

4. TITLE AND SUBTITLE

Quantifying Discretization Effects on Brain Trauma Simulations

5a. CONTRACT NUMBER

W911NF-10-2-0076

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Sophia M Haire

5d. PROJECT NUMBER

622618H80

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

American Society for Engineering Education

1818 N St NW, #600

Washington, DC 20036

8. PERFORMING ORGANIZATION REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

US Army Research Laboratory

ATTN: RDRL-WMP-C

Aberdeen Proving Ground, MD 21005-5069

10. SPONSOR/MONITOR'S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT NUMBER(S)

ARL-CR-0792

12. DISTRIBUTION/AVAILABILITY STATEMENT

Approved for public release; distribution is unlimited.

13. SUPPLEMENTARY NOTES

14. ABSTRACT

Numerical models of the brain are becoming an important tool for developing protective headgear for the Soldier. One 

nonphysical parameter that is often taken for granted in these models is the material representation by the mesh; arbitrarily

formed meshes can propagate error when resolving interactions among the skull, cerebrospinal fluid, and brain. We compared 

Lagrangian, pure Eulerian, and embedded Lagrangian/Eulerian schemes under both 2-dimensional cyclic rotation and planar 

blast scenarios. The 3 discretization techniques used differ not only in their finite element implementations but also necessarily 

in their representations of cerebrospinal fluid. The quantities analyzed in each case were the variations in stress magnitude, 

spatial distribution, and wave patterns that arise inside the brain. The effects of various time scales, degrees of rotation, and 

mesh resolutions were compared across methodologies to quantify the solutions and efficiencies of each technique and to 

understand how their formulations contribute to results. Preliminary results show that the wave patterns and spatial distributions 

for the rotational simulations are significantly different among the methodologies and that there is much more agreement among 

the blast loading simulations. Further investigation into specific aspects of several of the models is recommended.

15. SUBJECT TERMS

mesh effect, finite element, brain, interface, sliding 

16. SECURITY CLASSIFICATION OF:

17. LIMITATION

 OF

 ABSTRACT

UU

18. NUMBER

 OF 

 PAGES

48

19a. NAME OF RESPONSIBLE PERSON

Richard C Becker

a. REPORT

Unclassified

b. ABSTRACT

Unclassified

c. THIS PAGE

Unclassified

19b. TELEPHONE NUMBER (Include area code)

410-278-7980

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18
Approved for public release; distribution is unlimited.

iii

Contents

List of Figures iv

List of Tables v

Acknowledgments vi

Student Bio vii

1. Introduction and Background 1

2. Methods 1

3. Rotational Model Analysis 9

3.1 Results across Discretization Methods 9

3.2 Effects of Mesh Resolution on Results 14

3.3 Efficiency Analysis 20

4. Blast Loading Model Analysis 21

4.1 Results across Discretization Methods 21

4.2 Effects of Mesh Resolution on Results 25

4.3 Efficiency Analysis 31

5. Conclusions and Recommendations 31

6. References 33

Appendix. Material Parameters 35

Distribution List 38
Approved for public release; distribution is unlimited.

iv

List of Figures

Fig. 1 Example meshes for each discretization case. The meshes shown use 

Lagrangian, Eulerian, and embedded methods from top to bottom.......3

Fig. 2 Loading node-set for Eulerian rotational problem. The dark shaded 

area around the skull is the area to which the angular velocity loading 

was applied.............................................................................................4

Fig. 3 Load curve for rotational models...........................................................5

Fig. 4 Velocity, energy, and density inflow conditions for planar blast 

loading problem .....................................................................................6

Fig. 5 Example of planar blast loading problem. The dimensions of the air in 

the simulation are 900 mm in the x-direction by 450 mm in the ydirection. The example shown is pure Eulerian.....................................7

Fig. 6 Locations of tracer particles. The particles are located 5, 20, and 

50 mm in from the right and top inner edges of the skull. The example 

shown is a Lagrangian rotational model. The red and green materials 

represent the brain and skull, respectively.............................................8

Fig. 7 Effect of discretization method on pressure distribution. The rotation 

angles shown are 10.1° from equilibrium, 0.9° back toward 

equilibrium from the inflection point (displaced 59.1° from 

equilibrium), and 44.3° back toward equilibrium from the inflection 

point (displaced 15.7° from equilibrium). ...........................................10

Fig. 8 Effect of discretization method on von Mises equivalent stress 

distribution. Like in Fig. 7, the rotation angles shown are 10.1° from 

equilibrium, 0.9° back toward equilibrium from the inflection point 

(displaced 59.1° from equilibrium), and 44.3° back toward 

equilibrium from the inflection point (displaced 15.7° from 

equilibrium)..........................................................................................11

Fig. 9 Effect of discretization method on effective strain distribution. Like 

in Figs. 7 and 8, the rotation angles shown are 10.1° from equilibrium, 

0.9° back toward equilibrium from the inflection point (displaced 

59.1° from equilibrium), and 44.3° back toward equilibrium from the 

inflection point (displaced 15.7° from equilibrium). ...........................12

Fig. 10 Distance between master and slave nodes over time ...........................13

Fig. 11 Pressure over time by mesh resolution at 3 tracer locations, Lagrangian 

model. The graph legends indicate the number of elements used in 

each simulation. The red dashed lines display the inflection points in 

the simulation—that is, the approximate points at which the velocity 

changed directions. Data points were collected every 1.0 ms. ............15
Approved for public release; distribution is unlimited.

v

Fig. 12 Pressure over time by mesh resolution at 3 tracer locations, Eulerian 

model. As in Fig. 11, the graph legends indicate the number of 

elements used in each simulation, and the red dashed lines display the 

inflection points in the simulation—that is, the approximate points at 

which the velocity changed directions. Data points were collected 

every 1.0 ms.........................................................................................17

Fig. 13 Pressure over time by mesh resolution at 3 tracer locations, embedded 

model. As in Figs. 11 and 12, the graph legends indicate the number of 

elements used in each simulation. Data points were collected every 1.0 

ms.........................................................................................................19

Fig. 14 Wall time by mesh resolution of Lagrangian and Eulerian methods, 

rotational models..................................................................................20

Fig. 15 Pressure over time by discretization method at 6 tracer locations, 

highest resolution. The number of elements in the Eulerian, embedded, 

and Lagrangian simulations are 2,531,250; 2,532,018; and 311,994; 

respectively. Data points were collected every 1.0 ms. .......................22

Fig. 16 Pressure over time by mesh resolution at 3 tracer locations, Eulerian 

blast model. The graph legends indicate the number of elements used 

in each simulation. Data points were collected every 1.0 ms. .............26

Fig. 17 Pressure over time by mesh resolution at 3 tracer locations, embedded 

blast model. The graph legends indicate the number of elements used 

in each simulation. Data points were collected every 1.0 ms. .............28

Fig. 18 Pressure over time by mesh resolution at 3 tracer locations, Lagrangian 

blast model. The graph legends indicate the number of elements used 

in each simulation. Data points were collected every 1.0 ms. .............30

Fig. 19 Wall time by mesh resolution of each discretization method, blast 

models. For clarity, a projected linear trend curve is included for the 

Lagrangian model. ...............................................................................31

List of Tables

Table. Average element sizes within the head by simulation...........................9
Approved for public release; distribution is unlimited.

vi

Acknowledgments

The author wishes to acknowledge the mentorship, support, and patience of 

Dr Richard Becker. Funding for this project was provided by the Blast Protection 

for Platforms and Personnel Institute.
Approved for public release; distribution is unlimited.

vii

Student Bio

Sophia M Haire is a junior studying computer science and mathematics at the 

University of Maryland, Baltimore County. The summer of 2015 marked her 

second internship with the US Army Research Laboratory (ARL); her previous 

work consisted of software development and production for ARL shadowgraph 

range data collection. Her primary career and research goals lie in the fields of 

human-centered computing with a focus on assisting patients with brain disorders 

or injuries, including the elderly.
Approved for public release; distribution is unlimited.

viii

INTENTIONALLY LEFT BLANK.
Approved for public release; distribution is unlimited.

1

1. Introduction and Background

The Soldier Protection Sciences Branch of the US Army Research Laboratory 

(ARL) relies heavily on numerical models of the brain and surrounding tissues to 

understand load distributions, which will allow the design of efficient and effective 

protective equipment for the Army. A recent study of interest to the ARL 

community used both 2- and 3-dimensional (2-D) and (3-D) models to simulate the 

pressures and velocities within the brain when it was cyclically rotated about the 

axis of the neck at a targeted rate.

1 The models, meshed using Eulerian (fixed-grid) 

particle-based methods in the code Uintah, were validated via an experiment in 

which human subjects laid down, placed their heads in a device, and rotated their 

heads back and forth; tagged MRI was used to collect relevant data. Although the 

magnitudes and propagation of velocities between the empirical data and simulated 

results in the study were similar, the predicted pressure distribution in the brain was 

shifted in the results of the human experiments. Noticing that the complicated 

interface between the brain, cerebrospinal fluid (CSF), and skull was not well 

resolved in the study (properties were averaged among the 3 materials along the 

boundary), we hypothesized that the discretization scheme may skew such results. 

Our objective in this study is to quantify the effects that discretization methodology 

has on such model results by creating and analyzing a 2-D brain model with 

simplified geometry under 3 different discretization schemes and thus under 3

different methods of resolving the brain/CSF/skull interface. We conducted our 

study using both rotational boundary conditions and a planar blast-loading 

boundary condition for each discretization case.

2. Methods

To create models of the skull, cerebrospinal fluid, and brain, we relied on the 

software ALE3D from Lawrence Livermore National Laboratories.

2 ALE3D 

possesses an internal mesh generator that allows finite element meshes to be created 

and simulations run in Lagrangian, Eulerian, Arbitrary Lagrangian/Eulerian (ALE), 

and embedded Lagrangian/Eulerian schemes. In a Lagrangian simulation, the mesh 

is tied to the material as it moves or deforms; Eulerian simulations keep the mesh 

fixed while allowing material to advect through it; ALE methods move the mesh 

independently from the material motion, often depending on how distorted the 

mesh becomes over time; and the newly developed embedded method 

(implemented via the ALE3D FEusion package) allows coupling of fully 

Lagrangian and fully Eulerian objects within the same simulation. 
Approved for public release; distribution is unlimited.

2

We simplified the geometry of our brain model by reducing its shape to a 2-D 

ellipse with a semi-major axis of approximately 103 mm and a semi-minor axis of 

83 mm, surrounded by a 7-mm-thick skull. For the Lagrangian model, the CSF was 

accounted for using a frictionless laminar slide surface between the brain and skull.

The Eulerian model contained a 1-mm-thick CSF layer that was incorporated via 

mixed material zones in the mesh; in mixed material zones, material properties are 

averaged by volume in each element that contains multiple materials, and the 

boundaries between materials are reconstructed after each time step based on 

volume fractions. The embedded model comprises a Lagrangian skull and a 

Eulerian brain, and the CSF was represented via a frictionless slide surface between 

the two. In the pure Lagrangian and pure Eulerian rotational models, a void 

background was used due to the added simplicity it gave to the problem; in all other 

cases, a Eulerian mesh of air was used as the background. The Lagrangian blast 

loading problem required that embedded methodology be used to resolve the 

boundary between the air and the head. Figure 1 gives examples of each 

discretization type.
Approved for public release; distribution is unlimited.

3

Fig. 1 Example meshes for each discretization case. The meshes shown use Lagrangian, 

Eulerian, and embedded methods from top to bottom.
Approved for public release; distribution is unlimited.

4

The material properties for the air, CSF, brain, and skull were obtained from Dr 

Yolin Huang.

3 The air was modeled as a gamma law gas with coefficient

(γ – 1) = 0.4. The brain was modeled as an incompressible Mooney-Rivlin material 

with the follow equation of state for Cauchy stress: 

σ'= 

μ0

J �ω B�'

-(1- ω)�B�-1�

'

� = 

μ0

J �ω B�'

-(1- ω)((B�2)'- tr (B�) B�') � . (1)

The skull model used was a linear elastic model with a shear modulus of 2,660 MPa 

and a bulk modulus of 3,891 MPa. The units used throughout our simulations for 

material properties, coordinates, and results were based on grams (g), millimeters 

(mm), milliseconds (ms), and megapascals (MPa). These are different than the units 

traditionally used with ALE3D for blast simulations.

The head meshes in the rotational models were rotated using an angular velocity 

load curve applied to a node-set on the skull. In the Lagrangian model, the nodeset chosen was the very outside surface of the skull. In the embedded simulation, 

all sets of nodes inside the skull except for the inner and outer surfaces of the skull 

were used. The Eulerian model used a filled circular mesh object that surrounded 

the skull and was made of the same material as the skull to rotate the head by 

applying the angular velocity load curve to the set of all nodes within a circular 

radius outside of the skull but within the circular object (Fig. 2). This method of 

rotation was used because the Eulerian method prevented rotation of a noncircular 

node-set.

Fig. 2 Loading node-set for Eulerian rotational problem. The dark shaded area around the 

skull is the area to which the angular velocity loading was applied.
Approved for public release; distribution is unlimited.

5

To evaluate our rotational models, we ran each simulation for 520.0 ms with a 

loading curve represented by Fig. 3. In this time, the head rotates 60° to the left, 

120° to the right, and 60° left back to the center position. The rate was selected 

based upon android research performed by Robert Fitzpatrick,

4 which suggested 

that the maximum rate of safe rotation of a human head is 467°/s. For the embedded 

grid simulations, we were only able to run the model for 40.0 ms as opposed to 

520.0 ms because of the long run time of the simulations; so the loading curve, 

although the same shape as the one for the other 2 models, produced much higher 

velocities. Because of the time constraints of the project, no attempt was made to 

optimize the run time of the simulations. The 3-dimensionality of the embedded 

simulations and the Lagrangian blast simulation did not appear to significantly 

affect run time.

Table of Values

ms rad/ms

0 0

6.5 0.00815068

123 0.00815068

129.5 0

130 0

136.5 -0.00815068

383.5 -0.00815068

389.5 0

Fig. 3 Load curve for rotational models

The force applied in the blast loading models was simulated by applying velocity, 

density, and energy inflow conditions to the air node-set on the extreme left edge 

of the model, 340 mm to the left of the outer left-most edge of the skull. The inflow 

conditions were obtained by collecting time history values from a modified version 

of Dr Joel Stewart’s ShockTube model, which simulated a TNT explosion initiated 

at the apex of a rigid right-circular cone that was 14 m long.

5 The time histories 

were collected 2 m radially inward from the widest part of the cone and 

approximately halfway between the top edge of the cone and the midpoint axis of 

the cone in the vertical direction. The ShockTube simulation ran for 5.0 ms, and 

time histories were collected throughout the simulation. The time history data used 

were then restricted to the time in which a blast was actively moving past our 

selected tracer particle location.

Our blast loading simulations were thus run for 2.5 ms. Figure 4 depicts the inflow 

conditions used. 
Approved for public release; distribution is unlimited.

6

Fig. 4 Velocity, energy, and density inflow conditions for planar blast loading problem

Because of the planar symmetry of the blast loading problems, the problem was cut 

in half along the head’s major axis, and only the upper half of the problem was run 

(Fig. 5).
Approved for public release; distribution is unlimited.

7

Fig. 5 Example of planar blast loading problem. The dimensions of the air in the simulation 

are 900 mm in the x-direction by 450 mm in the y-direction. The example shown is pure 

Eulerian.

Data collection was performed using 6 tracer particles located at different locations 

inside the brain in each mesh; the locations of the tracer particles are shown in Fig.

6. For the purposes of our simulations, the right side of the brain is the side of the 

brain closest to the right side of the screen when viewed in a postprocessor during 

the first time cycle. The top side of the brain is the side closest to the top of the 

screen when viewed in a postprocessor during the first time cycle.

In one of the Lagrangian rotational simulations, 2 additional time history tracers 

were used to calculate the positions over time of a pair of adjacent master and slave 

surface nodes on the slide surface to determine how much slippage of the brain 

occurs during rotation.

Pressure (MPa)

0.50

0.25

0.00

200

400

0

-400 -200 0 200 400 INFLOW
Approved for public release; distribution is unlimited.

8

Fig. 6 Locations of tracer particles. The particles are located 5, 20, and 50 mm in from the 

right and top inner edges of the skull. The example shown is a Lagrangian rotational model. 

The red and green materials represent the brain and skull, respectively.

Each type of simulation was run at several different mesh resolutions; the average 

element sizes in the head in each case are given in the following Table. The number 

of elements given for each simulation is the total number of elements, including 

background void or air elements, in that simulation. Average element sizes are 

calculated for the rotation simulations and the Eulerian and embedded blast 

simulations by finding the percentage area of the head within the mesh and dividing 

the area of the head by the total number of mesh elements multiplied by this 

percentage. Average element sizes for the Lagrangian blast simulation were 

determined by finding the exact number of elements inside the head via temporarily 

removing the background air and dividing the area occupied by the head by the 

number of elements; this was done since the mesh resolution of the skull and brain 

in the Lagrangian blast simulation was much higher than the mesh resolution of the 

air. Since the embedded simulations and the Lagrangian blast simulation had to be 

represented in quasi-3-D with a 0.3-mm third-dimensional thickness, the element 

sizes given are the hypothetical sizes of their 2-D counterparts.
Approved for public release; distribution is unlimited.

9

Table. Average element sizes within the head by simulation

Rotation Blast

Lagrange

No. of 

Elements

Avg. Size

(mm2)

Lagrange

No. of 

Elements

Avg. Size

(mm2)

44,460 0.700 22,272 1.969

63,540 0.489 49,230 0.921

200,928 0.155 94,464 0.420

393,120 0.079 152,874 0.247

786,420 0.040 311,994 0.114

Euler

No. of 

Elements Avg. Size

Euler

No. of 

Elements Avg. Size

40,000 1.690 101,250 3.997

90,000 0.751 405,000 0.999

160,000 0.423 911,250 0.444

250,000 0.270 1,620,000 0.250

360,000 0.188 2,531,250 0.160

Embedded

No. of 

Elements Avg. Size

Embedded

No. of 

Elements Avg. Size

27,576 2.942 101,538 3.986

62,190 1.304 405,768 0.444

111,072 0.730 912,018 0.250

... ... 1,620,768 0.160

... ... 2,532,018 0.999

3. Rotational Model Analysis

3.1 Results across Discretization Methods

The distributions of the pressures, von Mises equivalent stresses, and effective 

strains in the brain at fixed angles of rotation were drastically different among the 

discretization schemes, as shown in Figs. 7–9. Again, the embedded simulation was 

run on a smaller time scale than the Lagrangian and Eulerian models. Also, the 

embedded model was rotated through atmospheric pressure to satisfy the 

requirement of the embedded calculations; thus, the range of the pressures and Von 

Mises equivalent stress for the embedded simulation shown in Figs. 7–9 has been 

adjusted for an initial internal pressure of 0.101325 MPa. The effective strain rate 

was given by the following equation:

eff. strain =�2

3 (ε11

2+ ε222+ 2ε12

2

) . (2)
Approved for public release; distribution is unlimited.

10

Fig. 7 Effect of discretization method on pressure distribution. The rotation angles shown 

are 10.1° from equilibrium, 0.9° back toward equilibrium from the inflection point 

(displaced 59.1° from equilibrium), and 44.3° back toward equilibrium from the inflection 

point (displaced 15.7° from equilibrium).

1.0

0.0

-1.0

Pressure (kPa)

1.0

0.0

-1.0

Pressure (kPa)

1.0

0.0

-1.0

Pressure (kPa)

12

0

-12

Pressure (kPa)

12

0

-12

Pressure (kPa)

12

0

-12

Pressure (kPa)

150

100

50

Pressure (kPa)

150

100

50

Pressure (kPa)

150

100

50

Pressure (kPa)
Approved for public release; distribution is unlimited.

11

Fig. 8 Effect of discretization method on von Mises equivalent stress distribution. Like in 

Fig. 7, the rotation angles shown are 10.1° from equilibrium, 0.9° back toward equilibrium 

from the inflection point (displaced 59.1° from equilibrium), and 44.3° back toward 

equilibrium from the inflection point (displaced 15.7° from equilibrium).

5.0

2.5

0.0

5.0

2.5

0.0

von Mises (kPa)

von Mises (kPa)

5.0

2.5

0.0

von Mises (kPa)

5.0

2.5

0.0

von Mises (kPa)

10

5

0

von Mises (kPa) von Mises (kPa)

von Mises (kPa)

5.0

2.5

0.0

10

5

0

von Mises (kPa)

10

5

0

von Mises (kPa)

5.0

2.5

0.0
Approved for public release; distribution is unlimited.

12

Fig. 9 Effect of discretization method on effective strain distribution. Like in Figs. 7 and 8, 

the rotation angles shown are 10.1° from equilibrium, 0.9° back toward equilibrium from 

the inflection point (displaced 59.1° from equilibrium), and 44.3° back toward equilibrium 

from the inflection point (displaced 15.7° from equilibrium).

One immediately visible result is that the pressure, stress, and strain waves 

propagating through the brain take on different shapes and frequencies in each case. 

In the Eulerian and embedded grid models, the pressure waves, Fig. 7, cause peak 

pressures to arise more interior to the brain than in the Lagrangian case. The 

embedded grid model, although unfortunately under-resolved because of the time 

constraints of this study, gives the closest results to the Uintah fixed-grid 

simulation. The Lagrangian model, oddly, appears not to have a clearly defined 

pressure wave but has pressures and stresses that gradually build up over time, as 

would be expected in a problem involving rigid body rotation. Figure 10 plots the 

distance between nodes in a master/slave node pair over time; the distance is 

coherent with what would be expected of a fluid rotating within a solid shell. 

0.4

0.2

0.0

0.4

0.2

0.0

Equivalent strain

Equivalent strain

0.4

0.2

0.0

Equivalent strain

0.4

0.2

0.0

Equivalent strain

0.8

0.4

0.0

Equivalent strain Equivalent strain

Equivalent strain

0.4

0.2

0.0

0.8

0.4

0.0

Equivalent strain

0.8

0.4

0.0

Equivalent strain

0.4

0.2

0.0
Approved for public release; distribution is unlimited.

13

Fig. 10 Distance between master and slave nodes over time

An analysis of what causes pressures to build instead of reflect in the Lagrangian 

model is ongoing, but the preliminary investigation has shown that the pressure 

increase can be completely attributed to volume changes in the brain. The change 

in energy at all time steps in our simulations is zero, so none of the volume changes 

are due to thermal expansion. It may be the case that the volume changes are due 

to roundoff error in our boundary conditions that causes the skull to change shape.

Another immediately visible result is that the pressures in the brain in the 

Lagrangian model are at least an order of magnitude higher than the pressures in 

the Eulerian model. The magnitudes of pressures observed in the embedded grid 

model cannot be subject to comparison with the other 2 methodologies because of

the shorter timescale over which rotation was applied. The scattered pressures 

within the skull in the Eulerian model arise from the exterior skull material used to 

rotate the head.

One notable result is that the highest von Mises stresses and effective strains in all 

cases occur very near to the center of the brain. This differs from the experimental 

results and brings up a number of questions as to why this phenomenon arises in 

our models. The most likely explanation for this is the geometry chosen for our 

simulations. The simplified, elliptical geometry chosen for our simulations neglects 

the folds and lobe patterns within a real brain structure. These folds and lobes 

provide a mechanism similar to hundreds of slide surfaces via which shear stresses 

arising from the rotation of the brain/skull interface are prevented from propagating 

far into the center brain. Without those folds and lobes providing a built-in 

mechanism for slip, the greatest shearing will occur at the center of the brain. 

0

20

40

60

80

100

120

140

160

180

200

1.00E+002.10E+014.10E+016.10E+018.10E+011.01E+021.21E+021.41E+021.61E+021.81E+022.01E+022.21E+022.41E+022.61E+022.81E+023.01E+023.21E+023.41E+023.61E+023.81E+024.01E+024.21E+024.41E+024.61E+024.81E+025.01E+02

Distance, mm

Time, ms

Distance between Master and Slave Nodes over 

Time, 44,460 Mesh Elements
Approved for public release; distribution is unlimited.

14

A concerning phenomenon present in the embedded grid model that is not well 

depicted by Fig. 7 is slight pressure loss and subsequent recovery over the course 

of the simulation. Visual analysis of results from the simulation show that by the 

time the brain reaches an inflection point, the average pressure in the brain has 

dropped by at least 0.02 MPa (the pressure has decreased to at least 80% of initial 

pressure). All pressure is regained by the end of the simulation; however, the 

distortion in results that this generated requires an understanding of what in the 

embedded grid method is causing this pressure flux. One hypothesis is that the 

reconstruction of the brain/skull boundary in the embedded calculations is subject 

to roundoff errors when rebuilding element volumes. Small losses in skull volume 

or stretching of the skull in the radial direction (and thus a miniscule volume 

increase in the brain) would likely lead to unexpected and unwanted pressure 

fluctuations. Another hypothesis is that the boundary conditions applied to the skull 

distort the shape of the skull and cause it to expand slightly; this phenomenon is 

similar to what we are currently investigating in the Lagrangian simulation, which 

is subject to comparable loading conditions.

3.2 Effects of Mesh Resolution on Results

Within discretization methods, it was often the case that the resolution of the mesh 

greatly affected the magnitudes and distributions of the results given by our models. 

An effect of mesh resolution on the Lagrangian simulations was that the magnitudes 

of the pressures observed varied significantly with resolution. Peak pressures at 

some points increased as resolution increased, and others did not. Figure 11 plots 

the pressures observed at 3 of the tracer locations at 3 different mesh resolutions 

over time.
Approved for public release; distribution is unlimited.

15

Fig. 11 Pressure over time by mesh resolution at 3 tracer locations, Lagrangian model. The

graph legends indicate the number of elements used in each simulation. The red dashed lines 

display the inflection points in the simulation—that is, the approximate points at which the 

velocity changed directions. Data points were collected every 1.0 ms.
Approved for public release; distribution is unlimited.

16

The peak pressure patterns and magnitudes, although not entirely congruent, are 

generally consistent across mesh resolutions and show a trend toward a lowerpressure magnitude than is given by our simulations. In all cases, as one nears the 

center of the brain, the order of magnitude of the pressure values decreases slightly, 

as would be expected for a rotating body. The remarkable uniformity in shape of 

the curves shows that not all the pressure effects are mesh related and that our model 

is capturing the physics of rotation of our head model without much interference. 

The resultant pressures in the Lagrangian rotational models are all an order of 

magnitude higher than the corresponding pressures in their Eulerian and embedded 

counterparts. As mentioned previously, the steady increase in pressures is 

completely dependent on changes in the volume of the brain that may in turn result 

from errors in our rotational boundary conditions; investigation into the cause of 

the volume change is ongoing.

In the Eulerian simulations, the peak pressure values were similar to each other, but 

the times at which peak values were reached differed. Also, at locations closer to 

the skull, the pressure patterns observed were inconsistent across resolutions.

However, like the Lagrangian simulations, the pressure curves appear to be 

converging to a single result as mesh resolution increases. Figure 12 displays this 

phenomenon.
Approved for public release; distribution is unlimited.

17

Fig. 12 Pressure over time by mesh resolution at 3 tracer locations, Eulerian model. As in 

Fig. 11, the graph legends indicate the number of elements used in each simulation, and the 

red dashed lines display the inflection points in the simulation—that is, the approximate points 

at which the velocity changed directions. Data points were collected every 1.0 ms.
Approved for public release; distribution is unlimited.

18

For the Eulerian simulations, the primary effect of mesh resolution appears to be 

the phase of the pressure curve as opposed to pressure magnitudes, and the phases 

become more in agreement among the models as one nears the center of the brain. 

Since the Eulerian method averages material properties over mixed-material 

elements (in our model, along the skull/CSF/brain boundary, which the 5-mm tracer 

was very close to), mesh resolution would be expected to be a large factor in the 

results along boundaries between materials and virtually a nonfactor far from 

mixed-material elements. As mesh resolution increases, the area over which 

averaging occurs decreases, potentially improving the accuracy of the Eulerian 

model. One of the most important pieces of Fig. 12 is that the boundary effects 

propagate toward the center of the brain; across resolutions, the highest pressures 

were located 20 mm from the skull as opposed to the normally expected 5 mm from 

the edge of the skull. 

In creating Figs. 11 and 12, sets of results for 5 different mesh resolutions were 

initially plotted, but this number was reduced to 3 to allow for better visualization 

of results. The high, low, and middle resolutions were kept in the graphs.

The unfortunate effects of under-resolution are manifested in the comparison of 

embedded grid results across mesh resolutions as shown in Fig. 13. Note that the 

zero pressure for these simulations is atmospheric pressure instead of 0.0 MPa.
Approved for public release; distribution is unlimited.

19

Fig. 13 Pressure over time by mesh resolution at 3 tracer locations, embedded model. As in 

Figs. 11 and 12, the graph legends indicate the number of elements used in each simulation. 

Data points were collected every 1.0 ms.
Approved for public release; distribution is unlimited.

20

The steepness of the spikes and drops in this figure seen in the meshes with lower 

resolutions can be attributed mainly to the checkerboarding in the simulation that 

resulted from the low mesh resolution. The amplitudes of the peaks decrease with 

increasing mesh resolution, implying that the simulations are converging to a 

gentler pressure curve that hovers very near atmospheric pressure. Like the 

Lagrangian simulations, the highest pressures seen in the embedded grid 

simulations are at the very outside of the brain.

3.3 Efficiency Analysis

The efficiencies of each discretization type with respect to simulation size are 

difficult to compare across all 3 methods; however, some analysis can be 

performed. Figure 14 displays the wall times of the Lagrangian and Eulerian 

simulations with respect to mesh resolution. All simulations were run on Excalibur 

on 64 processors.

Fig. 14 Wall time by mesh resolution of Lagrangian and Eulerian methods, rotational 

models 

In finite element analysis, it is usually expected that Lagrangian simulations run 

faster than Eulerian simulations because Lagrangian models do not have to account 

for advection across element boundaries. This is definitely the case for our models, 

and it appears that both methods provide a nearly linear run time with respect to 

input (mesh) size. Because of the time constraints of our project, we were only able 

to obtain results at 3 different mesh resolutions in the Lagrangian case.

0

0.005

0.01

0.015

0.02

0.025

0.03

0.035

0 100000 200000 300000 400000

Wall Time, seconds/cycle

Number of Mesh Elements

Wall Time per Cycle by Number of Mesh 

Elements, Rotation Simulations

Euler

Lagrange
Approved for public release; distribution is unlimited.

21

Unfortunately, because of the constraints of our simulations and the necessary 

restarts that had to be performed, the exact recorded wall times of the embedded 

simulations were lost. However, we can approximate that the embedded simulation 

with 111,072 elements ran for about 46 h to reach 40.0 ms; the Lagrangian 

simulation with approximately twice that many elements ran for around 14 h to 

reach 520.0 ms. The inefficiency posed by the embedded model was a major 

limiting factor on the quality of results available for this study.

4. Blast Loading Model Analysis

4.1 Results across Discretization Methods

Figure 15 depicts the results of the blast loading simulations using each 

discretization method at the highest resolution we were able to run.
Approved for public release; distribution is unlimited.

22

Fig. 15 Pressure over time by discretization method at 6 tracer locations, highest resolution. 

The number of elements in the Eulerian, embedded, and Lagrangian simulations are 

2,531,250; 2,532,018; and 311,994; respectively. Data points were collected every 1.0 ms.

0

0.05

0.1

0.15

0.2

0.25

Pressure, MPa

Time, ms

Pressure in Brain, 5 mm in from Right Edge of Skull, by 

Discretization Method, Highest Resolution

Euler

Embed

Lagrange

0

0.05

0.1

0.15

0.2

0.25

Pressure, MPa

Time, ms

Pressure in Brain, 20 mm in from Right Edge of Skull, 

by Discretization Method, Highest Resolution

Euler

Embed

Lagrange
Approved for public release; distribution is unlimited.

23

Fig. 15 Pressure over time by discretization method at 6 tracer locations, highest resolution. 

The number of elements in the Eulerian, embedded, and Lagrangian simulations are 

2,531,250, 2,532,018, and 311,994, respectively. Data points were collected every 1.0 ms

(continued).

-0.02

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

Pressure, MPa

Time, ms

Pressure in Brain, 50 mm in from Right Edge of Skull, 

by Discretization Method, Highest Resolution

Euler

Embed

Lagrange

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

Pressure, MPa

Time, ms

Pressure in Brain, 5 mm in from Top Edge of Skull, by 

Discretization Method, Highest Resolution

Euler

Embed

Lagrange
Approved for public release; distribution is unlimited.

24

Fig. 15 Pressure over time by discretization method at 6 tracer locations, highest resolution. 

The number of elements in the Eulerian, embedded, and Lagrangian simulations are 

2,531,250, 2,532,018, and 311,994, respectively. Data points were collected every 1.0 ms

(continued).

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

Pressure, MPa

Time, ms

Pressure in Brain, 20 mm in from Top Edge of Skull, by 

Discretization Method, Highest Resolution

Euler

Embed

Lagrange

0

0.05

0.1

0.15

0.2

0.25

Pressure, MPa

Time, ms

Pressure in Brain, 50 mm in from Top Edge of Skull, by 

Discretization Method, Highest Resolution

Euler

Embed

Lagrange
Approved for public release; distribution is unlimited.

25

Unlike in the rotation problem, the magnitudes and distributions of pressures in the 

blast loading simulations across discretization schemes were remarkably similar. 

Each simulation clearly showed 5 distinct pressure wave reflections within the skull 

over the course of the simulation; these reflections were observed at almost 

identically the same times and in the same places in each case. The magnitudes of 

pressures never deviated from atmospheric pressure by more than 0.1 MPa in any 

of the cases. The primary differences among the sets of results lie at the points 

closest to the brain/skull interface, as expected. The embedded grid methodology 

displays the shallowest wave peaks along the edge of the brain; the air required 

behind the skull for the embedded model may have had the effect of dampening the 

pressure wave. The sharpest and tallest wave peaks, as well as clearest pressure 

modes, are seen in the Lagrangian simulations; since Lagrangian meshing does not 

have mixed material zones, the averaging of properties seen in the other 2

methodologies may be dampening results by providing a gradient of densities over 

which the pressure wave has to pass before reflecting.

4.2 Effects of Mesh Resolution on Results

Mesh resolution of the blast loading meshes was not as large a factor in the pressure 

values recorded over time as it was for the rotational simulations. Figure 16 plots 

the results for the Eulerian blast loading model.
Approved for public release; distribution is unlimited.

26

Fig. 16 Pressure over time by mesh resolution at 3 tracer locations, Eulerian blast model. 

The graph legends indicate the number of elements used in each simulation. Data points were 

collected every 1.0 ms. 
Approved for public release; distribution is unlimited.

27

Pressure magnitude and wave reflection timing are extremely consistent across 

mesh resolutions for the Eulerian simulations. The “shakiness” observed in the 

values across resolutions at the 5-mm tracer particle is likely due to the mixedmaterial elements at the brain/CSF/skull interface; finer resolution of the CSF 

results in mixed elements comprising a smaller portion of the mesh, which should 

improve the accuracy of the simulation.

Mesh resolution appears to be slightly more of a factor in the embedded 

simulations’ results, as shown in Fig. 17.
Approved for public release; distribution is unlimited.

28

Fig. 17 Pressure over time by mesh resolution at 3 tracer locations, embedded blast model. 

The graph legends indicate the number of elements used in each simulation. Data points were 

collected every 1.0 ms.
Approved for public release; distribution is unlimited.

29

However, the convergence among the models is still very tight. Like the Eulerian 

blast simulation, the greatest divergence in pressure values is seen closest to the 

skull. In the case of the embedded simulation, this may be due to accumulated error 

in the FEusion calculation that is performed at the brain/skull boundary.

The Lagrangian blast loading models have the greatest results discrepancies with 

differing mesh resolution (Fig. 18); part of this may be due to the smaller numbers 

of overall elements than their Eulerian and embedded grid counterparts. The results 

appear to converge over time toward a clearly defined pattern of 5 wave reflections 

in the head over the course of the run. Once again, the highest discrepancies are 

closest to the skull. There are no mixed-material zones in the Lagrangian models; 

however, embedded grid methods were used to implement the background air, 

some of which lies behind the skull. It may be that the discretization of this 

calculation is a large factor in results. Another possible method of implementing 

the background air is to use a conformal mesh surrounding the head with a slide 

surface between the air and the skull. Unfortunately, there was not enough time 

available to create a conformal mesh during this study. Further investigation into 

different methods of implementing the background air would be required to 

understand this phenomenon more fully.
Approved for public release; distribution is unlimited.

30

Fig. 18 Pressure over time by mesh resolution at 3 tracer locations, Lagrangian blast model. 

The graph legends indicate the number of elements used in each simulation. Data points were 

collected every 1.0 ms.
Approved for public release; distribution is unlimited.

31

4.3 Efficiency Analysis

An interesting piece of the blast loading simulations is that contrary to what is 

normally expected in modeling, the Lagrangian simulation appears to be the least 

efficient for the number of mesh elements used (Fig. 19). This is probably due to 

the use of embedded grid methods to provide the background air for the run; it is 

recommended that a methodology that does not use an embedded methodology be 

implemented for the Lagrangian run. However, the Lagrangian simulations only 

required one embedded interface (between the skull and air), while the embedded 

simulations required 2 embedded interfaces (between the skull and air and between 

the brain and skull). In this planar situation, the Eulerian model is by far the quickest 

run.

Fig. 19 Wall time by mesh resolution of each discretization method, blast models. For

clarity, a projected linear trend curve is included for the Lagrangian model.

5. Conclusions and Recommendations

We have seen that in under both rotational and blast loadings, brain trauma 

simulation results are significantly affected by the discretization method used in 

creating the mesh. In the rotational problem, results using the same boundary 

conditions were disparate; there were little similarities among the pressure wave 

patterns observed in each case. The embedded grid model appears to have the most

spatially similar pressure patterns to the data from human trials.

1 However, 

although the preliminary results of the method look promising, the pressure loss 

during the simulation and the extremely long run times of the model make it less 

0

0.2

0.4

0.6

0.8

1

1.2

0 500000 1000000 1500000 2000000 2500000 3000000

Wall Time, seconds/cycle

Number of Mesh Elements

Wall Time per Cycle by Number of Mesh 

Elements, Blast Loading Simulations

Euler

Embed

Lagrange
Approved for public release; distribution is unlimited.

32

than desirable for the application we have presented. Investigation into the cause of 

the pressure loss and into techniques that improve the efficiency of the simulation 

are of paramount importance for future use of the embedded method in cases when 

the foreground and background meshes slide against each other. Also, the 

embedded model’s usefulness was limited by the low mesh resolutions required to 

complete the simulations on time; for further investigation, is it essential that a finer 

embedded mesh be used. The Lagrangian rotational simulation requires further 

study as well to determine the cause of the volume change that results in a pressure 

buildup that deviates from expected results. Our hypothesis is that the volume 

change results from error in the boundary conditions, but this hypothesis requires 

verification. 

The most pressing topic for further investigation across all methods is the effect of 

brain structure on the propagations of pressures, stresses, and strains. It is expected 

that a geometry which accurately represents the folds and lobes of the brain will 

have enough internal sliding mechanisms that the shear stresses will remain highest 

on the outside of the brain, which is very different from our results. For the blast 

loading problem, all 3 methods perform similarly and efficiently in areas far from 

the skull. An appropriate next step would be to conduct a study that provides a 

clearer physical explanation for the discrepancies in results near the brain/skull 

interface. Also, methods to implement the Lagrangian model that do not use 

embedded techniques to implement the air are of interest in future work both to 

isolate the Lagrangian method within the simulation and to improve the model’s 

run time. Unlike the rotational problem, it is not anticipated that incorporating an 

accurate brain geometry will affect results since no shear forces are involved in the 

simulation. 
Approved for public release; distribution is unlimited.

33

6. References

1. Daphalapurkar N, Ganpule S. Presentation: Dynamic shearing deformations in 

living human brain: modeling and validation. Baltimore (MD): The Johns 

Hopkins University. Presented to the Soldier Protection Sciences Branch, US 

Army Research Laboratory, Aberdeen Proving Ground, MD; 2015 Feb.

2. Nichols AL, editor. ALE3D users’ manual, arbitrary Lagrangian-Eulerian 3D 

and 2D multi-physics code. Livermore (CA): Lawrence Livermore National 

Laboratory; 2015 Sep. Report No.: LLNL-SM-662355.

3. Huang Y. Simulation of blast on a porcine head. Aberdeen Proving Ground 

(MD): Army Research Laboratory (US); 2015 July. Report No.: ARL-TR7340.

4. Fitzpatrick R. Designing and constructing an animatronic head capable of 

human motion programmed using face-tracking software [master’s thesis]. 

[Worcester (MA)]: Worcester Polytechnic Institute; 2010 May 1.

5. Stewart J. US Army Research Laboratory, Aberdeen Proving Ground, MD. 

Personal communication, July 2015.
Approved for public release; distribution is unlimited.

34

INTENTIONALLY LEFT BLANK.
Approved for public release; distribution is unlimited.

35

Appendix. Material Parameters 

 

This appendix appears in its original form, without editorial change.
Approved for public release; distribution is unlimited.

36

Void

voidinput ss0 1000

Air

matinput

 czero 2.0 qfb 0.0 crq 0.1 pmin -1.0e-01

 eosvmin 1.0e-03 eosvmax 1.0e+03 vhlimit 1.0e+05 

 rho 1.18e-6 cvav {(720.)*(1.18e-6)} p0eosvmax 1

 evfrom tp

 p0 1.01325e-1

 t0 298.

 e0 2.533125e-01

 v0 1.0

 koinput

 isol 0 iform 5

 coef 0.4

Brain

 matinput

 rho 1.04e-3 e0 0.0 v0 1.0

 pmin -2.0e2 epsfail 25 

 czero 2.0 crq .1 qfb 0.15 linq 1 

 cvav 2.926 eosvmin 0.3 eosvmax 2.0 

 t0 300.0 

 msinput 

 ysmodel 146

 shr_mod 13.e-3 w0 0.3 order 2 adv_log 1

 elasmodel 99

 hardmodel 299 

 eosmodel 304

 rhoc2 2200 s1 2.56 s2 -1.986 s3 .2268 g0 0.5 a 0.

 eosvmin 0.3 eosvmax 2.0 

 ec0 -877.8

 em0 800.

 failmodel 499

Cerebrospinal Fluid

 # y0 from Bloomfield et. al., 1998 -> assumed viscosity and 

density of water

 matinput

 rho 1.04e-3 e0 0.0 v0 1.0

 pmin -2.0e2 epsfail 25 

 czero 2.0 crq .1 qfb 0.15 linq 1 

 cvav 2.926 eosvmin 0.3 eosvmax 2.0 

 t0 300.0 
Approved for public release; distribution is unlimited.

37

 msinput 

 ysmodel 100

 elasmodel 0

 shr_mod 2660

 hardmodel 201

 y0 8.9e-7 ybet 0 n 0.1 ymax 5 rate_a 0 rate_b 1 

rate_exp 1 

 eosmodel 304

 rhoc2 2200 s1 2.56 s2 -1.986 s3 .2268 g0 0.5 a 0.

 eosvmin 0.3 eosvmax 2.0 

 ec0 -877.8

 em0 800.

 failmodel 499

Skull

 matinput

 rho 1.412e-3 e0 0.0 v0 1.0 

 cvav 3.772 

 msinput 

 ysmodel 100

 elasmodel 0

 shr_mod 2660

 hardmodel 201

 y0 4 ybet 0 n 0.1 ymax 5 

 eosmodel 304

 rhoc2 3891 s1 0.94 s2 0 s3 0 g0 1.0 a 0.

 ec0 -1131.6

 em0 800.

 failmodel 499
Approved for public release; distribution is unlimited.

38

1 DEFENSE TECHNICAL

(PDF) INFORMATION CTR

DTIC OCA

2 DIRECTOR

(PDF) US ARMY RESEARCH LAB

RDRL CIO LL

IMAL HRA MAIL & RECORDS

MGMT

1 GOVT PRINTG OFC

(PDF) A MALHOTRA

1 UMBC

 (PDF) S HAIRE

10 DIR USARL

(PDF) RDRL WMP A

S R BILYK

RDRL WMP B

C HOPPEL

Y HUANG

S SATAPATHY

S WOZINAK

A SOKOLOW

RDRL WMP C

R BECKER

T BJERKE

RDRL WMP G

S KUKUCK

R EHLERS