EXPERIMENTAL SYSTEMS FOR EX-VIVO AND IN-VITRO STUDIES OF

BLAST-INDUCED TRAUMATIC BRAIN INJURY

By

Suhas Vidhate

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Mechanical Engineering — Doctor of Philosophy

2021

ABSTRACT

EXPERIMENTAL SYSTEMS FOR EX-VIVO AND IN-VITRO STUDIES OF

BLAST-INDUCED TRAUMATIC BRAIN INJURY

By

Suhas Vidhate

Blast-induced traumatic brain injury (bTBI) has been widely accepted as a “signature”

wound affecting military service members in modern conflicts. When a blast-wave gener-

ated by an improvised explosive device explosion propagates through the human head, it

is hypothesized to cause direct mechanical damage to the brain tissue leading to vascular

injury, cerebral edema as well as less detectable but persistent deficits. However, the exact

mechanisms and pathophysiology of bTBI still remain poorly understood. One of the main

reasons for such poor understanding is the technical challenge of reproducing the typical

time-varying loading cycles induced on brain tissue after a blast event under controlled lab-

oratory conditions. Blast events have a sub-millisecond onset of high pressure followed by

complex dynamics resulting from the interaction between the blast wave and the intricate

anatomical structure of the human head. This interaction gives rise to time-varying intracra-

nial pressure profiles, dynamic shear loads, and cavitation events, which are hypothesized to

cause mechanical insults to the brain tissue.

To tackle these experimental challenges, we have developed four experimental devices

that can reproduce the dominant intracranial effects induced by blast loading onto the brain.

These experimental devices are: (1) Multi-material Hopkinson bar (MMHB) actuator, (2)

Generator of broadband pressure cycles, (3) Rig for controlled-cavitation events, and (4)

Generator of dynamic shear cycles. This dissertation details the design and development

of the aforementioned experimental systems along with the preliminary ex-vivo and in-vitro

experiments performed to test their applicability.

The designed apparatuses are comparably inexpensive, compact, easily portable, and

highly controllable, making them well suited for biomedical applications. These devices can

be used to conduct ex-vivo and in-vitro experiments involving animal brain tissue specimens,

cell cultures, and organoids to explore their pathophysiological response to the blast-like

loadings observed during a bTBI event.

Dedicated to my beloved aunt Aruna and my wife Rasika for their unconditional love,

support, and sacrifices . . .

iv

ACKNOWLEDGMENTS

My academic journey to PhD at Michigan State University (MSU) has been possible only

due to the strong support and guidance from many special people at different stages of my

life. I want to take this opportunity to express my sincere gratitude to them all.

First and foremost, I want to thank my advisor Prof. Ricardo Mej´ıa-Alvarez for his

continuous guidance and support throughout my PhD tenure. Dr. Mej´ıa-Alvarez always

encouraged my ideas and gave me the opportunity to work autonomously. He was the first

to applaud me for accomplishments, provided constructive criticism when required, and

steered me in the right direction during setbacks. His immense knowledge and experience

helped bring my ideas into reality and personally motivated me greatly.

I am also deeply grateful to my dissertation committee members Dr. Adam M. Willis

(Maj. USAF), Prof. Thomas J. Pence, and Prof. Ranjan Mukherjee for their invaluable

supervision and assistance. Their inputs and suggestions always persuaded me to explore

beyond the obvious and bring ingenuity to my work.

I want to thank my research collaborators Dr. Joseph Beatty and Dr. Charles Lee Cox

from the Department of Physiology at MSU; Dr. Michael March, Dr. Patrick Sleiman, and

Dr. Hakon Hakonarson from the Center for Applied Genomics at the Children’s Hospital

of Philadelphia (CHOP); Prof. Zane Lybrand from the Department of Biology at Texas

Woman’s University; and Nikolas Merlock from the University of Texas at San Antonio for

their association in conducting the multi-disciplinary research. Their collaborations facili-

tated me to demonstrate the applicability of my doctoral research.

Within MSU, I express my gratitude to Dr. Craig Gunn for all the technical writing

guidance and Mechanical Engineering Department Staff for their organizational support.

v

And I want to thank MSU, U.S. Air Force Research Laboratory (award numbers FA8650-

18-2-6880 and FA865-17-2-6836), and the Qatar National Research Fund (a member of the

Qatar Foundation, NPRP grant No. 8-2424-1-477) for sponsoring my research.

I have had the special support and encouragement from my colleagues Joseph Kerwin,

Atacan Y¨ucesoy, Ankit Gautam, Tushar Kailu, Faezeh Masoomi, Cl´ement Roy, and Paul

Sandherr through several technical contributions and wonderful collaborations. I am grateful

for the cherished time we spent together in laboratory and socially. I also want to thank my

friends Saikat Mondal, Deepak Kumar, Gaurav Chauda, and Dr. Saranraj Karuppuswami

for their companionship which turned my PhD journey into a joyful ride.

I particularly want to thank my schoolteachers Mr. Vivek Supare and Mr. Surendranath

Reddy for cultivating my special interest in Mathematics and Physics. I also want to thank

my professors at the Indian Institute of Technology Bombay, Dr. Salil S. Kulkarni, Dr.

Milind Atrey, and Dr. V. Kartik, for giving me the early research exposure and inspiring me

to pursue doctoral studies.

Finally, a special thanks to my family. I want to remember my late father Mr. Jeevan

Vidhate and my late mother Mrs. Anita Vidhate for all the love and values they deeply

enrooted into me.

I also want to remember my late grandparents Dr. Laxman Vidhate

and Mrs. Kanta Vidhate for their cheers and blessings.

I want to thank my uncles Mr.

Vasudeo Vidhate and Mr. Gajanan Vidhate, and all of my family whom I owe a great deal.

I would not have been here without their support and sacrifices. I want to appreciate my

friends Visat Patel and Shubham Yadav for the motivating conversations about business,

technology, and research which always kept me enthusiastic about my work. Lastly, I want

to thank my dear wife Rasika for her patience, support, and companionship, which helped

me to complete this journey with ease.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi

Chapter 1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1
Introduction and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Experimental Techniques for Impulsive Pressure Generation . . . . .
1.2.2 Experimental Techniques for Generating Cavitation . . . . . . . . . .
1.2.3 Experimental Techniques for Dynamic Shear Loadings
. . . . . . . .
1.3 Aims and Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Devices for Complex Pressure Loadings . . . . . . . . . . . . . . . . .
1.3.2 Device to Generate Cavitation . . . . . . . . . . . . . . . . . . . . . .
1.3.3 Device for Dynamic Shear Cycles . . . . . . . . . . . . . . . . . . . .

Chapter 2 Multi-Material Hopkinson Bar Actuator

. . . . . . . . . . . . .
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1
Stress Wave Propagation in Multi-Material Elastic Bars . . . . . . . .
2.2.2 Multi-material Hopkinson bar (MMHB) Actuator . . . . . . . . . . .
2.2.3 Theoretical Background and Numerical Modeling . . . . . . . . . . .
2.2.3.1 One-dimensional Wave Propagation in Solids
. . . . . . . .
2.2.3.2 Method of Characteristics for Modeling One-dimensional Wave
Propagation in a Multi-Material Bar . . . . . . . . . . . . .
2.2.4 Optimization Methodology . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4.1 Design Parameters . . . . . . . . . . . . . . . . . . . . . . .
2.2.4.2 Optimization Algorithm . . . . . . . . . . . . . . . . . . . .
2.2.4.3 Objective Functions
. . . . . . . . . . . . . . . . . . . . . .
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Case Study-I: Single Frequency Sine Wave . . . . . . . . . . . . . . .
2.3.2 Case Study-II: Single-frequency decaying sine wave . . . . . . . . . .
2.3.3 Case Study-III: Sum of two sinusoidal waves . . . . . . . . . . . . . .
2.3.4 Case Study-IV: Sum of three sinusoidal waves . . . . . . . . . . . . .
2.3.5 Case Study-V: Sum of three decaying sinusoidal waves
. . . . . . . .
2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 3 Generator of Broadband Pressure Cycles . . . . . . . . . . . . .
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Working Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Device Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
7
8
9
10
11
11
12
12

13
13
16
17
20
21
21

24
29
29
31
32
34
34
36
38
39
42
44

47
47
48
48
49

vii

3.3 Control

3.2.3
3.2.4 Numerical Model

3.2.2.1 Design Iteration I . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2.2 Design Iteration II . . . . . . . . . . . . . . . . . . . . . . .
Instrumentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
System Identification . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1
Single Pulse Waveform . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Friedlander Waveform . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 Multi-modal Waveform . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Conclusion and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.3.1
3.3.2 Feedforward Control

Chapter 4 Bench-top Rig for Controlled-Cavitation Experiments . . . . .
4.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Working Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Device Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Instrumentation and Control . . . . . . . . . . . . . . . . . . . . . . .
4.2.3
4.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49
52
53
55
59
59
65
65
66
68
71
73

76
76
77
77
79
80
81
85

5.1
5.2 Design Iteration I

Chapter 5 Generator of Dynamic Shear Cycles . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 Working Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 Device Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.3
Instrumentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.4 Performance Characterization . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.5 Control
System Identification . . . . . . . . . . . . . . . . . . . . . .
5.2.5.1
5.2.5.2 Feedforward Control
. . . . . . . . . . . . . . . . . . . . . .

87
87
88
88
89
90
91
97
97
99
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.3.1 Working Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.3.2 Device Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Instrumentation and Control . . . . . . . . . . . . . . . . . . . . . . . 106
5.3.3
5.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

5.2.6 Limitations
5.3 Design Iteration II

Chapter 6 Preliminary Ex-vivo and In-vitro Experiments

. . . . . . . . . 118
6.1 Ex-vivo Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
6.1.1 Ex-vivo Experiments involving Time-varying Pressure Cycles . . . . . 119
6.1.1.1 Experimental Procedure . . . . . . . . . . . . . . . . . . . . 119
6.1.1.2 Results
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

viii

6.2

6.1.2 Ex-vivo Experiments involving Dynamic Shear Loadings

. . . . . . . 123
6.1.2.1 Experimental Procedure . . . . . . . . . . . . . . . . . . . . 123
6.1.2.2 Results
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
In-vitro Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
In-vitro Experiments involving HeLa cell cultures . . . . . . . . . . . 130
6.2.1
6.2.1.1 Experimental Procedure . . . . . . . . . . . . . . . . . . . . 130
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
6.2.1.2 Results
In-vitro Experiments involving Cerebral Organoids
. . . . . . . . . . 134
6.2.2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 134
6.2.2.2 Experimental Procedure . . . . . . . . . . . . . . . . . . . . 138
6.2.2.3 Results
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

6.2.2

Chapter 7 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
7.1 Future Directions for the MMHB Actuator . . . . . . . . . . . . . . . . . . . 143
7.2 Future Directions for the Pressure Apparatus
. . . . . . . . . . . . . . . . . 144
7.3 Future Directions for the Cavitation Rig . . . . . . . . . . . . . . . . . . . . 145
7.4 Future Directions for the Shear Apparatus . . . . . . . . . . . . . . . . . . . 146

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

ix

LIST OF TABLES

Table 2.1: Material library for multi-material bar components. The material proper-
ties for Material ID 1-5 were taken from [1], and Material ID 6-12 were
obtained from the ANSYS material library [2].
. . . . . . . . . . . . . . .

30

Table 2.2: Design variables for MMHB actuator . . . . . . . . . . . . . . . . . . . . .

30

Table 2.3: RMSE values and the optimized actuator design parameters for the MMHB
. . . . . . . . . . . . . . . . . . . . . . . . . . .

optimization case studies

44

Table 3.1: Material properties for the numerical model of the pressure apparatus

. .

56

Table 3.2: Input parameters of the sine-sweep system identification experiments con-
. . . . . . . . . . . . . . . . . . . . . .

ducted for the pressure apparatus

60

Table 3.3: Input excitation pulses constructed for generating pressure pulses

. . . .

66

Table 3.4: Sinusoidal components of multi-modal target pressure profile

. . . . . . .

71

Table 5.1: Input parameters of the sine-sweep system identification experiments con-
. . . . . . . . . . . .

ducted for the shear apparatus (Design iteration I)

97

Table 6.1: Loading protocols for in-vitro experiments involving HeLa cell cultures . . 132

x

LIST OF FIGURES

Figure 1.1:

Illustration showing a blast wave generated from an explosion and the
resulting injuries from the blast exposure (adapted from [3] ) . . . . . . .

Figure 1.2: A Friedlander waveform having a peak pressure of 450 kPa and a positive
. . . . . . . . . . . . . . . . . . . . . . . . . . .

phase duration of 10 ms

Figure 1.3:

Interaction between a blast wave and the human head (axial-view shown).
The reflection and the transmission of the incident wave and the resulting
skull deformation are also illustrated (adapted from [4]).
. . . . . . . . .

Figure 1.4:

(a) The time-domain and (b) the frequency-domain representation of the
intracranial pressure profile at the contrecoup region (shown in figure 1.3)
resulting from a blast exposure [5, 6]
. . . . . . . . . . . . . . . . . . . .

Figure 1.5:

Illustration of skull flexure resulting from blast wave-skull interaction
(adapted from [7])
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.6: Experimental setup of a split Hopkinson pressure bar (Picture Credits:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

REL Inc.)

Figure 1.7: Experimental setup of a shock-tube (Picture Credits: Center for Injury
Biomechanics, Materials, and Medicine at New Jersey Institute of Tech-
nology)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.1:

(a) Schematic of the three bar configurations used to illustrate the effect
of impedance mismatch. (b) Incident stress pulse (σL) specified as the
boundary condition at the left end x = 0.
. . . . . . . . . . . . . . . . .

Figure 2.2: The stress profiles (σR) at the fixed far end of (a) single material bar, (b)
two-material bar, and (c) three-material bar depicted in figure 2.1a . . .

2

3

3

4

6

8

9

18

19

Figure 2.3: Schematic of the MMHB actuator

. . . . . . . . . . . . . . . . . . . . .

20

Figure 2.4: A spring-dashpot representation of the standard-linear-solid (SLS) mate-
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

rial model

22

Figure 2.5: Stress wave propagation through the interface between materials A and B 22

xi

Figure 2.6: Characteristics solution for (a) Left boundary point, (b) Interior point,
and (c) Right boundary point of the MMHB actuator. The subscript
of the spatial element (x) indicates the location of grid points: 1 and n
denote the left and right boundary, respectively, and i denote an internal
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
point.

25

Figure 2.7: Flowchart of the optimization algorithm . . . . . . . . . . . . . . . . . .

32

Figure 2.8:

(a) Time domain and (b) frequency domain comparison between the target
pressure profile and the output pressure profile generated by the optimized
MMHB actuator for case study-I
. . . . . . . . . . . . . . . . . . . . . .

Figure 2.9:

(a) Time domain and (b) frequency domain comparison between the target
pressure profile and the output pressure profile generated by the optimized
MMHB actuator for case study-II . . . . . . . . . . . . . . . . . . . . . .

Figure 2.10: The set of optimization solutions obtained for the case study-III. These
solutions are plotted in the objective space where the axes denote the
optimization objective functions . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.11: (a) Time domain and (b) frequency domain comparison between the target
pressure profile and the output pressure profile generated by the optimized
. . . . . . . . . . . . . . . . . . . . .
MMHB actuator for case study-III

Figure 2.12: (a) Time domain and (b) frequency domain comparison between the target
pressure profile and the output pressure profile generated by the optimized
MMHB actuator for case study-IV . . . . . . . . . . . . . . . . . . . . .

Figure 2.13: (a) Time domain and (b) frequency domain comparison between the target
pressure profile and the output pressure profile generated by the optimized
MMHB actuator for case study-V . . . . . . . . . . . . . . . . . . . . . .

Figure 3.1: A piston-cylinder assembly (l0: initial length of the water volume inside
the cylinder, up: displacement of the piston) . . . . . . . . . . . . . . . .

Figure 3.2:

(a) Cross-sectional schematic view and (b) Fabricated design of the pres-
. . . . . . . . . . . . . . . . . . . . .
sure apparatus (Design iteration I)

35

37

39

40

41

43

48

50

Figure 3.3: The complete assembly of the pressure apparatus (Design iteration I) . .

51

Figure 3.4:

(a) Cross-sectional schematic view and (b) Fabricated design of the pres-
sure apparatus (Design iteration II) . . . . . . . . . . . . . . . . . . . . .

53

xii

Figure 3.5:

(a) Pressure apparatus fixture consisting of aluminum plates and steel
threaded rods along with the supporting V-block and (b) the complete
assembly of the pressure apparatus (Design iteration II)
. . . . . . . . .

Figure 3.6: One-dimensional domain for the numerical model simulating the pressure
chamber with white dots showing the locations of the nodes used for ex-
tracting pressure profiles
. . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.7: Comparison between the actuator displacement profile (left y-axis) spec-
ified as the input boundary condition to the numerical model and the
output pressure profiles (right y-axis) obtained from the numerical simu-
lation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

55

58

Figure 3.8: A chirp signal

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

Figure 3.9:

(a) Chirp voltage profile given as an input to the pressure apparatus for
a system identification experiment (b) Output pressure profile generated
as a response to the chirp excitation . . . . . . . . . . . . . . . . . . . .

61

Figure 3.10: (a) Frequency-domain representation of the chirp input profile (b) Frequency-

domain representation of the output pressure profile

. . . . . . . . . . .

62

Figure 3.11: (a) Gain curves and (b) phase difference curves for the pressure apparatus
. . . . . . . . . . . .

determined from system identification experiments

Figure 3.12: (a) Sinusoidal excitation pulses given as inputs to the apparatus. (b) Out-
put pressure pulses generated by the apparatus for corresponding input
pulses.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.13: Comparison between the pressure pulse generated by our apparatus and
. . . . . . . . .

the pressure profile obtained from SHPB experiment [8]

Figure 3.14: A Friedlander waveform with 450 kPa maximum pressure and 10 ms pos-
itive phase duration . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.15: Comparison between a Friedlander like input excitation waveform given to
the apparatus (left y-axis) and the output pressure profiles (right y-axis)
generated by the apparatus
. . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.16: Multi-modal target pressure profile constructed for an output tracking
problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

67

68

69

70

71

Figure 3.17: Control input profile determined by the feedforward controller . . . . . .

72

xiii

Figure 3.18: Comparison between the target pressure profile and the experimental pres-
. . . . . . . . . .

sure profiles obtained using the feedforward controller

Figure 4.1: A piston-cylinder assembly filled with water (l0: initial length of the water
volume inside the cylinder, up: displacement of the piston) . . . . . . . .

73

78

Figure 4.2: Experimental setup of the cavitation apparatus

. . . . . . . . . . . . . .

79

Figure 4.3: Comparison between the target piston displacement profile and the mea-
sured piston displacement profile for the cavitation experiment (left y-
axis). The force applied by the motor to achieve the target piston dis-
. . . . . . . . . . . . . . . . . . .
placement is also shown (right y-axis).

Figure 4.4: Displacement of the piston (left y-axis) and the resulting pressure profile
generated inside the chamber (right y-axis). Some of the key instances
of the cavitation event are highlighted on the pressure profile. The corre-
sponding high-speed images of the pressure chamber are shown in figure
4.5.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.5: High-speed images of the pressure chamber showing the expansion and
collapse of the cavitation bubbles. The pressure inside the chamber at
the corresponding instances is highlighted on the pressure profile shown
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
in figure 4.4.

81

83

84

Figure 5.1: Schematic of the shear apparatus (Design iteration I)

. . . . . . . . . .

89

Figure 5.2: Fabricated setup of the shear apparatus (Design iteration I)

. . . . . . .

90

Figure 5.3: A PAA gelatin specimen . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

Figure 5.4: Current profile given as an excitation input to the voice-coil actuator . .

92

Figure 5.5: High-speed images showing the shear deformation of PAA gelatin specimen 93

Figure 5.6: Processed images of the gelatin specimen overlaid with displacement field
vectors at time (A) t = 0 ms, (B) t = 5 ms, (C) t = 10 ms, (D) t = 15 ms,
. . . . . . . . . . . . . . . . . . . . .
(E) t = 20 ms, and (F) t = 25 ms

94

Figure 5.7: Measured top plate displacement for the first two loading cycles . . . . .

95

Figure 5.8: Comparison of the shear strain profiles computed using our in-house de-
veloped PIV script, the PIVLab toolbox, and the measured top plate
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
displacement profile

96

Figure 5.9: Shear strain-rate profile for the first two loading cycles . . . . . . . . . .

96

xiv

Figure 5.10: (a) Chirp voltage profile given as an input to the shear apparatus for a sys-
tem identification experiment. (b) Output displacement profile generated
as a response to the chirp excitation . . . . . . . . . . . . . . . . . . . .

98

Figure 5.11: (a) Gain curves and (b) phase difference curves for the shear apparatus

(Design iteration I) determined from system identification experiments

. 100

Figure 5.12: Multi-modal target displacement profile constructed for an output track-

ing problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

Figure 5.13: Control input profile determined by the feedforward controller . . . . . . 101

Figure 5.14: Comparison between the target displacement profile and the experimental

displacement profiles obtained using the feedforward controller

. . . . . 102

Figure 5.15: (a) Front-view and (b) top-view of the fabricated setup of the shear ap-

paratus (Design iteration II) . . . . . . . . . . . . . . . . . . . . . . . . . 105

Figure 5.16: Comparison between the target displacement profiles and the experimen-
tal displacement profiles generated by the voice-coil actuator for a loading
frequency of (a) 10 Hz, (b) 50 Hz, and (c) 100 Hz . . . . . . . . . . . . . . 108

Figure 5.17: Sinusoidal shear strain profiles generated by the shear apparatus (Design

iteration II) for a loading frequency of (a) 10 Hz, (b) 50 Hz, and (c) 100 Hz 109

Figure 5.18: Shear strain (left y-axis) and shear stress (right y-axis) profiles for PAA

gelatin shear tests at (a) 20 Hz and (b) 50 Hz loading frequency . . . . . 111

Figure 5.19: Variation of the complex shear modulus of PAA gelatin with loading fre-

quency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

Figure 5.20: Variation of the loss factor of PAA gelatin with loading frequency . . . . 114

Figure 5.21: Variation of the storage and loss modulus of PAA gelatin with loading

frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

Figure 6.1: Mouse coronal brain slices in aCSF . . . . . . . . . . . . . . . . . . . . . 119

Figure 6.2: Pressure chamber (Design iteration I) containing the mouse brain slices

in an aCSF environment

. . . . . . . . . . . . . . . . . . . . . . . . . . 120

Figure 6.3: Pressure apparatus (Design iteration I) with the mouse brain slices placed

inside the chamber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

Figure 6.4: Pressure loading profile exerted onto the mouse brain slices . . . . . . . . 121

xv

Figure 6.5: Electrophysiological recordings of the mouse brain slices

. . . . . . . . . 122

Figure 6.6: The shear strain profiles (loading frequency: 50 Hz; shear strain ampli-
tudes: 10% and 20%) exerted onto the brain tissue specimens (the first
10 out of the total 50 loading cycles are shown) . . . . . . . . . . . . . . 124

Figure 6.7: High-speed images of the brain tissue specimen under dynamic shear loads:
(a) undeformed state (0% shear strain) and (b) deformed state (20% shear
strain) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

Figure 6.8: Change in electrophysiological resting membrane potential of neurons re-

sulting from dynamic shear loadings of various strain amplitudes . . . . . 127

Figure 6.9: Change in electrophysiological input resistance of neurons resulting from

dynamic shear loadings of various strain amplitudes . . . . . . . . . . . . 128

Figure 6.10: Change in electrophysiological time constant of neurons resulting from

dynamic shear loadings of various strain amplitudes . . . . . . . . . . . . 128

Figure 6.11: Frequency of neuronal action potential firing versus current injected curves

(F-I curves) for three dynamic shear loading conditions . . . . . . . . . . 129

Figure 6.12: Pressure chamber (Design iteration I) containing a glass coverslip plated

with HeLa cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

Figure 6.13: LDH assay results for in-vitro experiments involving HeLa cell cultures . 132

Figure 6.14: “Dead-Green” staining results for in-vitro experiments involving HeLa cell

cultures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

Figure 6.15: Cerebral organoids reproduce human cortical organization. (A) Timeline
for stem cell protocol to grow cerebral organoid. (B) Compared to the
human cortex, cerebral organoids contain most of the same layered orga-
nization and cells present. (C) Characterization of cortical structure from
deep and upper-layer neurons in developing cortical spheroid (dotted line).
(D) Image of the ventricular zone (VZ), astrocytes, and oligodendrocytes
that form a support niche. . . . . . . . . . . . . . . . . . . . . . . . . . . 136

Figure 6.16: Neurophysiology of the cerebral organoids.

(a) Cerebral organoid dis-
playing a complex network of neurons. (b) Multielectrode array (MEA)
measures 12 channels from each organoid. (c) Neural activity is measured
from the MEA. The bottom image is a magnified time scale of the voltage
signal recorded in the top image. (d) The raw signals can be filtered, and
network oscillations are measured to determine the synchrony of activity
recorded in cerebral organoids.

. . . . . . . . . . . . . . . . . . . . . . . 137

xvi

Figure 6.17: Blasting cerebral organoids.

(a) Timeline of blast experiment (b) Im-
age of 3 cerebral organoids inside of pressure chamber. (c) The pressure
waveform recorded inside the chamber during the exposure (Amplitude
= 250 kPa, Frequency = 3 kHz, and duration = 8 ms). (d) 24 well plate
being loaded with organoids after being blasted. (e) Empty MEA well
showing the layout of 12 electrodes. (f) MEA well filled with media and
cerebral organoid on electrodes.

. . . . . . . . . . . . . . . . . . . . . . 139

Figure 6.18: Blast acutely stimulates cerebral organoids. (a) Raster plot of spike events
recorded from MEA in no blast control organoids. (b) Raster plot of spike
events recorded from blast group at 1 hr post-blast. (c) Quantification of
spike amplitudes in 1 hr and 24 hr post-blast organoids. (d) Quantification
of spike number in cerebral organoids.

. . . . . . . . . . . . . . . . . . . 141

Figure 6.19: Blast acutely increases network synchrony. (A) Two channels with raw
signal filtered using a 5 Hz Morlet Wavelet from control organoids. (B)
The phase angle difference plotted for each time unit between two channels
(black), and the average angle and phase synchronization plotted (green
bar) from control organoids. (C) Filtered channels from blast organoids.
(D) Plot of phase angle difference (black) and average phase synchro-
nization (green) plotted from blast organoids. (E) Quantification of phase
synchronization between the control and blast organoids at 1 hr post-blast.
(n=3 organoids per group; p-values less than 0.05 is considered statisti-
cally significant)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142

xvii

Chapter 1

Introduction

1.1

Introduction and Motivation

Blast-induced traumatic brain injury (bTBI) is one of the major health issues affecting

soldiers on the battlefield. It has been a leading cause of mortality and disability sustained

by armed forces and law enforcement personnel. The statistics of U.S. military casualties in

the operations held in Iraq and Afghanistan reported 253,300 traumatic brain injury cases

between 2000 and 2012 [9]. Consequently, this issue has caused a financial burden of billions

of dollars to the health care system [10].

As stated in [11–13], traumatic brain injuries associated with blast exposure follow four

distinct categories: primary injuries are caused by the interaction of the blast wave directly

with the brain tissue leading to mechanical damage; secondary injuries result from projectile

(shrapnel) penetration into the head; tertiary injuries are caused by blunt impacts to the head

as a result of being propelled away by the explosion; and quaternary injuries include other

forms of injury such as polytrauma, thermal, and chemical burns [14, 15]. The four injury

categories and their causes are illustrated in figure 1.1. Out of all these injuries, primary

bTBI is the least understood and thus have been the main focus of this dissertation.

The blast-wave generated from an improvised explosive device (IED) explosion typically

follows the characteristics of a Friedlander waveform [16,17]. A Friedlander waveform consists

1

Figure 1.1: Illustration showing a blast wave generated from an explosion and the resulting
injuries from the blast exposure (adapted from [3] )

of an instantaneous rise in pressure followed by an exponential decay and, lastly, a negative

pressure phase [16]. The time duration for which the pressure is above the atmospheric

pressure is known as the positive phase duration of the waveform. As an illustration, a

Friedlander waveform having a peak pressure of 450 kPa and a positive phase duration of

10 ms is shown in figure 1.2.

When a human is exposed to a blast wave generated from an IED explosion, the blast

wave interacts with the skull and the intricate anatomical structure inside the head. One

such interaction is illustrated in figure 1.3, where reflection and the transmission of the head-

on incident blast wave through the human head (axial-view) is shown. This interaction alters

the blast wave’s Friedlander wave nature and gives rise to complex and time-varying intracra-

nial pressure (ICP) profiles. The comprehensive numerical studies [5, 6] have simulated a

2

Figure 1.2: A Friedlander waveform having a peak pressure of 450 kPa and a positive phase
duration of 10 ms

Figure 1.3: Interaction between a blast wave and the human head (axial-view shown). The
reflection and the transmission of the incident wave and the resulting skull deformation are
also illustrated (adapted from [4]).

blast wave generated by a detonation of 2.3 kg of C4 at a stand-off distance of 2.3 m, and

analyzed the mechanical loads experienced by brain tissue due to the blast exposure. These

3

05101520253035404550Time (ms)-1000100200300400500Pressure (kPa)studies reported that the ICP profiles are complex, highly dynamic, and exhibit spatial vari-

ation. One such ICP profile experienced by the brain’s rear region (contrecoup region) with

respect to the incident head-on blast wave is shown in figure 1.4a. The frequency-domain

representation of the same ICP profile is shown in figure 1.4b. This ICP profile has a wide

(a)

(b)

Figure 1.4:
(a) The time-domain and (b) the frequency-domain representation of the
intracranial pressure profile at the contrecoup region (shown in figure 1.3) resulting from a
blast exposure [5, 6]

4

01234567Time (ms)-500-400-300-200-1000100200300400500Pressure (kPa)0246810Frequency (kHz)010203040506070Pressure Amplitude (kPa)spectral distribution, which contributes to its complex and highly-dynamic nature. It can be
observed that most of the power is carried by the frequencies in the 0 − 5 kHz range, which
has significantly higher amplitudes compared to the frequencies in the 5 − 10 kHz range.

The multi-material structure of the head offers interfaces with impedance mismatch for

the transmission of the incident blast wave. Because of the impedance mismatch at each

material interface, a part of the incident stress wave reflects, and a part of it is transmit-

ted. This transmission and reflection repeat at every intermediate interface and result in

a complex ICP profile because of the wave superpositions [7, 13]. A blast wave also causes

the skull to deform dynamically (shown in figure 1.3), creating localized regions of high and

low pressure and large pressure gradients that sweep through the brain [18]. This unsteady

deformation of the skull also contributes to the complexity of the ICP profiles [7, 18, 19].

The blast wave interaction with the skull and the intricate anatomical structure inside the

head also results in traveling shear waves. The unsteady deformation of the skull resulting

from reflection and transmission of the incident blast wave at the air-skull interface excites

the skull’s multi-modal response. This skull flexure phenomenon is illustrated in figure 1.5.

The skull flexure resulting from a blast exposure additionally contributes to traveling shear

waves inside the brain soft tissue [18]. As with the case of longitudinal stress waves, the

propagation of shear waves is significantly altered by the impedance mismatch at material

interfaces, which results in energy dissipation at the interfacial zones and causes focal me-

chanical damage [13]. Moreover, the transmission and attenuation of shear waves through

the human brain highly depend on the frequency of the propagating shear waves [20].

The negative ICP resulting from the multiple transmissions-reflections of the incident

blast wave, skull flexure, or the relative motion between the brain and the skull can cause

cavitation in the cerebrospinal fluid (CSF) [7, 13, 21]. Cavitation phenomenon is observed

5

Figure 1.5: Illustration of skull flexure resulting from blast wave-skull interaction (adapted
from [7])

when a liquid is subjected to a tensile or negative pressure higher than a threshold, i.e.,

the liquid’s vapor pressure at that temperature [22]. Cavitation encompasses the formation,

growth, and collapse of the bubbles caused by the rapid changes in pressure [23]. The collapse

of the CSF cavitation bubbles generates high-pressure shock-waves and impinging jets that

can be damaging to the surrounding brain tissue [24–27].

Primary blast injuries are hypothesized to be linked with diffuse axonal injury (DAI),

vascular tears, intracranial hemorrhage, and focal cortical contusion [28]. The mechanical

damage resulting from the negative phase of ICP profiles is often associated with the CSF

cavitation [26,29,30]. The local shearing of the brain tissue resulting from pressure waves and

shear strain loads at the material interfaces is observed to be consistent with brain injuries

6

such as DAI near material interfaces [18,28]. Though brain soft-tissue damage resulting from

shear [29,31], cavitation [26,29,30], traumatic cell deaths, and disrupted cell structures [8,32]

have been identified as possible injury mechanisms involved in bTBI, the exact mechanisms

associated with primary blast exposure and their neuropathological consequences still remain

ambiguous.

One of the major contributing factors for such ambiguity is the difficulty in replicating

the mechanical loads observed during a blast exposure to induce mechanical insults within

brain tissue specimens or cell cultures. These mechanical loads include the complex ICP

profiles, the dynamic shear loadings induced by traveling shear waves, and the high-pressure

shock-waves resulting from the collapse of the CSF cavitation bubbles. To fully understand

the underlined pathological effects of these mechanical insults on brain tissue specimens or

cell cultures, there is a need for experimental devices that can generate such complex loading

cycles under controlled laboratory conditions.

1.2 Literature Review

Although several damage mechanisms associated with bTBI have been reported in the lit-

erature [33], this dissertation mainly focuses on three dominant damage mechanisms: (1)

impulsive, multi-frequency, and fast-varying pressure cycles, (2) cavitation and,(3) dynamic

shear loadings. This section discusses the experimental techniques that have been previously

employed to generate above mentioned mechanical loadings.

7

1.2.1 Experimental Techniques for Impulsive Pressure Generation

Previously, researchers have implemented experimental techniques such as split-Hopkinson

pressure bars (SHPB) [8, 32, 34, 35], shock-tubes [36–42], and blast-tubes [43–46] to generate

pressure profiles mimicking a blast loading. In the SHPB technique, a traveling stress wave

is generated in the incident bar from a striker impact. This wave is then used to compress

water and generate impulsive pressures in the piston-cylinder assembly placed at the end of

the bar. Shock-tubes use a sudden release of compressed gases to generate impulsive pres-

sures. Furthermore, blast-tubes use confined explosive detonations to generate impulsive

pressures. These experimental setups generate pressure profiles that have a single overpres-

sure pulse either sustained or similar to a Friedlander waveform. Other researchers have

used microdetonics to generate blast loading [47,48]. Representative experimental setups for

an SHPB and a shock-tube are shown in figures 1.6 and 1.7.

All of the above-mentioned techniques have limited control over the pressure profiles that

they can generate, are not easily tunable to replicate the ICP dynamics, and are, in general,

very expensive, requiring significant technical expertise and large laboratory spaces.

Figure 1.6: Experimental setup of a split Hopkinson pressure bar (Picture Credits: REL
Inc.)

8

Figure 1.7: Experimental setup of a shock-tube (Picture Credits: Center for Injury Biome-
chanics, Materials, and Medicine at New Jersey Institute of Technology)

1.2.2 Experimental Techniques for Generating Cavitation

There are several experimental techniques reported in the literature employed to gener-

ate cavitation events. Some of the key techniques are a modified SHPB [25, 26, 30, 32, 49,

50], shock-tube [36, 37], drop tower [51, 52], ultrasound [53, 54], laser [27, 55], and a laser-

ultrasound combination [56]. The modified SHPB setups use the reflected stress-pulse to

generate negative pressures and induce cavitation. The shock-tubes and drop towers de-

liver acceleration-induced or inertial cavitation. The impact from a shock-wave or a drop

weight creates a relative motion between the liquid and its surrounding solid boundary,

which generates a negative pressure region at the boundary and induces cavitation. The

SHPB, shock-tube, and drop tower setups are relatively complex and require rigorous tech-

nical training for safe operation. Also, they are extremely bulky and thus, lack portability.

The ultrasound technique uses focused ultrahigh-frequency (typically greater than 1 MHz)

acoustic waves to apply periodic compression-tension cycles. When the acoustic wave am-

plitude is sufficiently high, it generates cavitation bubbles during the negative phase. These

bubbles partially collapse during the subsequent positive phase of the acoustic wave. The

9

laser-based setups use a focused laser beam pulses to generate cavitation bubbles. The ul-

trasound and the laser techniques generate focused cavitation events with extremely short

timescales (order of few microseconds). The timescales and the frequency range of cavitation

events generated by laser and ultrasound techniques significantly differ from the cavitation

associated with bTBI, which typically lasts for several milliseconds.

1.2.3 Experimental Techniques for Dynamic Shear Loadings

To induce dynamic shear loading at high strain-rates, researchers have implemented a mod-

ified SHPB that includes a parallel plate fixture to shear the specimen [57,58]. However, the

shear strain magnitudes generated using the modified SHPB technique are usually very high

and result in specimen tearing. This high strain deformation does not replicate the low shear

strains observed during a primary bTBI event. The shear deformation linked to bTBI has low

strain amplitudes but high strain-rates [10]. Moreover, the SHPB technique cannot generate

oscillatory shear deformation associated with traveling shear waves. Rotational rheometers

are generally used to induce oscillatory shear deformation into a specimen and can be used

to generate shear cycles of low strain amplitudes [59]. However, they are optimized for ma-

terial characterization tasks and cannot be tuned to generate user-specified shear loading

profiles. Additionally, rheometers lack the ability to generate shear deformations at high

oscillating frequencies and strain-rates and are, in general, very expensive. Researchers have

also developed their own custom shear instruments that use linear actuators and a parallel

plate assembly to load soft tissue-like gelatin specimens [60–63]. These reported instruments

could generate either the shear strain amplitudes or loading frequencies typically associated

with a bTBI event, but not both simultaneously.

10

1.3 Aims and Outline

An ideal solution to the experimental challenges involved in reproducing mechanical loads of

a blast event in a controlled manner would be the following: create completely controllable,

portable, reproducible, and inexpensive methods to generate arbitrarily complex pressure

cycles, dynamic shear loadings, and cavitation events onto living tissue and cell cultures.

Building apparatuses that could be easily shared between laboratories would also accelerate

bTBI research. Our research demonstrates the design and development of such experimental

devices that can generate tunable, reproducible, and complicated loading cycles.

1.3.1 Devices for Complex Pressure Loadings

In this dissertation, Chapter 2 discusses the concept of a Multi-Material Hopkinson bar

(MMHB) actuator to generate complicated stress-wave patterns. The MMHB actuator re-

places the incident bar of a conventional split-Hopkinson pressure bar apparatus with a

multi-material incident bar to transform the incident stress pulse into a complex loading

profile. As proof of concept, we have designed the MMHB actuator for generating various

sinusoidal pressure profiles with multi-frequencies and decaying characteristics.

Chapter 3 discusses the design and development of the apparatus for generating broad-

band pressure cycles. This apparatus uses a water-filled piston-cylinder assembly driven by a

piezoelectric actuator to generate complex and fast-varying pressure profiles. The apparatus

also incorporates a feedforward controller to track user-specified target pressure profiles. The

experimental capabilities of this apparatus are demonstrated by generating a single pressure

pulse with various pulse-widths and magnitudes, an approximate Friedlander waveform, and

a multi-modal waveform.

11

1.3.2 Device to Generate Cavitation

Chapter 4 focuses on a bench-top rig designed for controlled-cavitation experiments. This

cavitation apparatus also uses a water-filled piston-cylinder assembly, which is operated

by a high-speed linear motor to generate negative pressures and induce cavitation. The

capabilities of this apparatus are demonstrated by performing a cavitation experiment with

water as the working liquid.

1.3.3 Device for Dynamic Shear Cycles

Chapter 5 discusses the design and development of the apparatus for generating dynamic

shear cycles. This apparatus uses a voice-coil actuator and two parallel plates assembly to

generate oscillatory shear loadings at various strain amplitudes and frequencies. The pre-

liminary experiments conducted to determine the mechanical properties of polyacrylamide

(PAA) gelatin specimens, as well as to test the implementation of the apparatus for loading

samples of living tissue and three-dimensional cell cultures, are also discussed in this chapter.

Chapter 6 discusses the preliminary ex-vivo and in-vitro experiments conducted to explore

the response of biological specimens to complex time-varying pressure and shear loadings.

These experiments served a dual purpose: determining the pathophysiological response of

the animal brain tissue specimens, cell cultures, and cerebral organoids to the blast-like in-

tracranial effects and testing the implementation of our experimental devices for bTBI-related

studies. Lastly, Chapter 7 summarizes the development of the aforementioned experimental

systems for loading biological specimens and provides directions for future developments.

12

Chapter 2

Multi-Material Hopkinson Bar

Actuator

2.1

Introduction

In stress wave analysis and experimental characterization, the assumption of homogeneity

and continuity is mostly inadequate for contemporary engineering practices. Specifically,

there is a need to better understand the formation of complex wave profiles and the response

of the continua to these wave profiles. Conventional testing systems are not fully capa-

ble of generating controlled complex stress waves with predefined spectral characteristics

comparable to the application requirements.

Systems containing multiple materials and material interfaces typically exhibit complex

stress waveforms when subject to impact or blast-like impulsive loads. Structures made up of

layered composite materials are examples of this class of systems. The multiple transmissions

and reflections that the incident wave experiences because of the layered configuration result

in a complicated loading [64], and this phenomenon is utilized to design impact-resistant

protective equipment such as armors and helmets [65].

The response of biological systems to exogenous dynamic loadings is another phenomenon

where complex stress waveforms are observed. One such case is that of blast-induced Trau-

13

matic Brain Injury (bTBI), which is one of the major reasons for the morbidity and mortal-

ity of service-members on the battlefield [9, 66]. While the nature of the incident load has

a significant effect on the magnitude and pattern of injury [67], previous studies show that

the multi-layered morphology of the head, skull, and intracranial contents, with different

impedances for wave propagation, greatly affect the peak amplitude and the spectral char-

acteristic of the intracranial pressure (ICP) profile [7, 18, 68]. Moreover, not only the peak

stress amplitude but also the temporal frequency of the stress profile could contribute to tis-

sue damage [68]. Due to the geometrical and material anisotropy in the brain tissues at the

micron-level, central nervous system (CNS) components such as neuron and glial structures

might produce different responses to the same blast exposure [69, 70]. As a result, bTBI

is diffused in nature due to the heterogeneity of the skull’s contents and the complexity of

intracranial wave interactions. In consequence, the degree to which different spectral com-

ponents participate in brain damage is still an open question. There is also a discrepancy

in blast exposure data obtained from different laboratories because of the variations in the

experimental techniques and protocols implemented to simulate blast [71]. This poses a

challenge while comparing bTBI-related experimental studies from different research groups

and encourages an experimental technique that can consistently replicate the amplitude and

frequency characteristics of stress experienced in human bTBI. Such an experimental tech-

nique might replicate blast-like loading conditions that produce TBI neuropathologies seen

in post-trauma.

The split-Hopkinson pressure bars (SHPB) technique has been extensively used to in-

vestigate the mechanical response of materials at high strain-rate deformation [34, 72]. The

details of the SHPB technique are well-known. SHPB offers a diverse arrangement of spec-

imen, bar, and loading mechanism for various high strain-rate material characterization

14

studies. Additionally, the SHPB has been previously used to load biological tissues with

blast-like pressure profiles [8,32]. However, the conventional SHPB is still challenging to use

on specimens with low strength, stiffness, or impedance, such as biological tissues, foams,

gels, and polymers [72–74]. Another limitation is that the loading profile generated using a

conventional SHPB is typically a single stress wave pulse [34], despite the various available

types of accelerators, e.g., electrostatic, pneumatic, explosive, and pulsed laser [75]. The

impact of the striker generates a stress wave that propagates through the incident bar; part

of the incident stress pulse is transmitted into the specimen and the transmission bar, while

the remaining part reflects into the incident bar. The reflected pulse travels back and forth

along the incident bar as it reflects consecutively on both ends of the bar. But this effect

does not load the sample cyclically in a typical SHPB test, because the sample and trans-

mission bar tend to separate from the incident bar upon the first wave interaction. Thus, a

conventional SHPB is not capable of generating a complex, broadband stress wave akin to

ICP profiles.

In this chapter, we report on a promising methodology to design an actuator that can

generate repeatable cycles of multi-frequency stress waves at strain-rate levels typical of SH-

PBs. Because the loading generated by a conventional SHPB is a single-frequency stress-wave

pulse, we consider the more general concept of a Multi-Material Hopkinson Bar (MMHB)

actuator. Our contention is that such a device is capable of generating loading profiles hav-

ing complicated spectral characteristics as per the user specifications. In this concept, the

incident bar is substituted for multiple bars of different materials arranged in series. The

distribution of materials and lengths of these bars are chosen so as to obtain the characteris-

tics of the targeted loading profile. The underlying idea behind the multi-material incident

bar is to transform the single stress wave impact into a targeted complex wave by leveraging

15

the impedance mismatch between adjacent bars. As the incident wave propagates through

the multi-material bar, each material interface produces a reflected and a transmitted wave.

This process repeats multiple times at each material interface, inducing wave interactions

that result in a more complex output at the end of the multi-material bar. The material

and geometrical features of each bar component constituting the multi-material form the

parameter design space and require a selection that is based on the targeted profile. A

multi-objective optimization algorithm is employed to acquire the multi-material incident

bar configuration that can best replicate the targeted wave profile in both the frequency and

time domains. This design optimization is achieved by coupling the numerical modeling of

the MMHB wave dynamics to the optimization algorithm.

This chapter is organized as follows. Section 2.2 discusses the working principle of the

MMHB actuator along with a numerical model analyzing stress wave propagation in multi-

material bars. The optimization algorithm and the design parameters required to develop

a target profile-specific actuator are also discussed here. The results from a series of case

studies conducted to test the applicability of the MMHB concept are presented in section

2.3. The chapter concludes in section 2.4 with a discussion on the capabilities and limitations

of the MMHB actuator.

2.2 Methods

The stress wave propagation through an interface results in a reflected and a transmit-

ted wave due to the impedance mismatch, as detailed in section 2.2.3.1. A multi-material

bar constructed by joining multiple bars of different materials naturally has interfaces with

impedance mismatch, which create reflections and transmissions of the incident stress wave.

16

The superposition and interference of these wave reflections and transmissions transform a

single localized stress pulse input into a relatively complex and delocalized stress wave at

the end of the bar.

2.2.1 Stress Wave Propagation in Multi-Material Elastic Bars

The principle of generating complex stress waves using a multi-material bar is demonstrated

by analyzing longitudinal stress wave propagation for the three cases shown in figure 2.1a.

These cases are: (1) a single material bar made from aluminum, (2) a two-material bar

made from aluminum and steel, and (3) a three-material bar made from aluminum, steel,

and copper.

The resulting stress history at the end of the multi-material bar is obtained by solving

the governing equations given in section 2.2.3.1. These equations within each bar, along with

the appropriate interface equations between constituent bars, are treated numerically with

MATLAB using the method of characteristics (see section 2.2.3.2). The material properties

used for the numerical modeling presented here are listed in table 2.1.

A stress profile (σL), such that shown in figure 2.1b, is specified on the left end of the

multi-bar assembly as a boundary condition to simulate an impact from a striker. This type

of incident stress pulse is typically observed during an impact from a spherical striker [30].

The velocity at the right end of the multi-bar assembly (V R) is specified as zero to impose

the fixed-end boundary condition. As a sign convention, compressive stress is considered

positive.

The stress profiles (σR) generated at the fixed right-end of the multi-bar assembly is

extracted as output from the MATLAB script. For the case studies depicted in figure 2.1a,

the corresponding computed right-end stress profiles are shown in figures 2.2a, 2.2b, and

17

2.2c. As the number of bar components within the multi-material bar increases, so does

the number of interfaces with impedance mismatch giving rise to an increasing reflection

(a)

(b)

Figure 2.1: (a) Schematic of the three bar configurations used to illustrate the effect of
impedance mismatch. (b) Incident stress pulse (σL) specified as the boundary condition at
the left end x = 0.

18

Three-Material BarTwo-Material BarSingle-Material Barx=1000 mmx=1000 mmx=1000 mm00.20.40.60.81Time (ms)-1-0.500.511.522.53 Stress L (MPa)(a)

(b)

(c)

Figure 2.2: The stress profiles (σR) at the fixed far end of (a) single material bar, (b)
two-material bar, and (c) three-material bar depicted in figure 2.1a

and transmission cascade. Accordingly, the stress profile obtained at the right end of the

bar also increases in complexity (as is evident from figures 2.2a, 2.2b, and 2.2c). The case

studies depicted in figure 2.1a treated the metallic bar materials as purely elastic. The more

19

00.511.522.533.54Time (ms)-6-4-20246Stress R (MPa)00.511.522.533.54Time (ms)-8-6-4-202468Stress R (MPa)00.511.522.533.54Time (ms)-6-4-20246Stress R (MPa)general methodology presented here treats multi-material bars consisting of both elastic and

viscoelastic materials, the latter of which allows explicitly for dissipation.

2.2.2 Multi-material Hopkinson bar (MMHB) Actuator

Motivated by the previous considerations, we seek to replicate a given complex loading

cycle at the right end of the multi-material bar by optimally sequencing bars of different

materials and lengths. Because our primary motivation is the study of bTBI, we focus on

an experimental setup that is appropriate for applying loads on biological samples. More

generally, however, the principles are general ones that can be applied for broader purposes

of material characterization.

The multi-material incident bar shown in figure 2.3 acts as the piston in a piston-cylinder

assembly used for unsteady compression of water. Its movement produces a highly complex

pressure cycle inside the cylinder. This system is basically a pressure chamber for loading a

biological specimen. Water is used as the working fluid inside the pressure chamber due to

its nearly incompressible nature.

Figure 2.3: Schematic of the MMHB actuator

20

Multi-Material Incident BarBar Component-1Bar Component-2Bar Component-NLoading DeviceGas GunStrikerWaterPressure Chamber(Piston-Cylinder Assembly)The actuator design problem involves several design variables such as the lengths of the

various bar components, material selection for those components, and the striker impact

parameters. The design is now optimized in order to generate a pressure profile similar to

the user-specified target pressure profile.

2.2.3 Theoretical Background and Numerical Modeling

2.2.3.1 One-dimensional Wave Propagation in Solids

This section provides a brief review of the essential theoretical background in order to aid

in the understanding of the following sections.

One-dimensional longitudinal wave propagation in a long bar is governed by the momen-

tum balance equation

and the continuity equation

ρ

∂V
∂t

=

∂σ
∂x

,

∂V
∂x

=

∂

∂t

.

(2.1)

(2.2)

The dependent variables, V , σ, and 
 denote the velocity, stress, and strain respectively

at a given location x and time t. The constant ρ denotes the density of the bar material.

Additionally, the relation between the stress and the strain is defined by the constitutive

model of the bar material, which is taken to be viscoelastic.

Wave propagation in a viscoelastic bar has two distinct characteristics: attenuation and

dispersion [1, 76]. Attenuation corresponds to the reduction in stress wave amplitude as the

wave travels along the length of the bar. Dispersion is associated with the change in either

the shape or the spectral character of the stress pulse as it propagates. Here, the standard-

21

linear-solid (SLS) model is used as the constitutive model. It is modeled as a combination

of springs and a dashpot, as shown in figure 2.4.

Figure 2.4: A spring-dashpot representation of the standard-linear-solid (SLS) material
model

The stress-strain relationship for the SLS model is given by

∂σ
∂t

+

σ
θ

= (E1 + E2)

∂

∂t

+ E1



θ

,

(2.3)

where E1 is the long-term (fully relaxed) modulus, E1 +E2 is the short-term (instantaneous)

modulus, and θ is the relaxation time constant [1,77]. The linear elastic model is the special

case of the SLS model with E1 = E and E2 = 0 , where E is the Young’s modulus of the

elastic bar material.

The wave propagation phenomenon through an interface between two bars made up of

materials A and B is shown in figure 2.5. When a stress wave (σI) traveling through material

A reaches the interface where the two bars are bonded, the difference in the material proper-

ties on either side of the interface creates an impedance mismatch for the wave propagation.

Figure 2.5: Stress wave propagation through the interface between materials A and B

22

η=E2E2E1Since the two bars are bonded together, it is assumed that there is no relative motion between

the bars at the interface, implying continuous velocity and stress fields across the interface.

Part of the incident wave is reflected (σR), while the rest of the wave is transmitted into the

next bar (σT).

If both materials A and B are elastic then the amplitudes of the reflected (σR) and

transmitted (σT) stress waves are determined as

σR =

ZB − ZA
ZB + ZA

σI

and σT =

2ZB

ZB + ZA

σI,

where ZA and ZB denote the respective acoustic impedances of materials A and B [78–80].

This acoustic impedance is given by

Z =(cid:112)E × ρ.

As expected, perfect impedance match (ZB = ZA) gives no reflected wave and gives a
transmitted wave that exactly matches the incident wave. More generally (ZB (cid:54)= ZA) the

transmitted wave has the same compressive or tensile character as the incident wave whereas
the reflected wave can change its character depending on the sign of ZB − ZA.

The numerical method described next in section 2.2.3.2 provides a treatment where mate-

rial interfaces are always nodes in the (x, t) discretization. The procedure naturally recovers

the above interface relations, if both materials are elastic. When one or both materials are

truly viscoelastic then the numerical method provides an appropriate extension.

23

2.2.3.2 Method of Characteristics for Modeling One-dimensional Wave Propa-

gation in a Multi-Material Bar

Since an elastic material model can be treated as a special case of the SLS model, a numerical

model is developed to simulate the wave interaction in a multi-material bar made up of

viscoelastic bar components. The governing equations for one-dimensional wave propagation

in a viscoelastic bar are treated using the method of characteristics as described in [1].

In an (x, t) domain governed by equations (2.1), (2.2), and (2.3), there are two sets of

characteristic curves

dx
dt

(cid:115)

= ±

E1 + E2

ρ

≡ ±CV ,

(2.4)

where CV denotes the maximum (leading edge) signal speed in the viscoelastic material. The

(+) sign in equation (2.4) corresponds to the waves propagating in the positive x-direction
while the (−) sign corresponds to the waves propagating in the negative x-direction. The

corresponding compatibility conditions along the characteristics are given by

dV = ± 1
ρCV

dσ ± σ − E1

ρCV θ

dt.

(2.5)

The third set of characteristics is non-propagating and along it the compatibility condition

holds:

dx = 0, and d
 =

dσ

E1 + E2

σ − E1

(E1 + E2)θ

+

dt.

(2.6)

Quantities E1, E2, ρ, θ, and CV depends on location x as determined by the multi-material

bar under consideration.

To solve the compatibility equations along the characteristics, the spatial domain for each

24

of the bar components is discretized to have

∆x = CV × ∆t,

(2.7)

where ∆x denotes the distance between two adjacent nodes, and ∆t denotes the time-step.

Additionally, the length of each of the bar components constituting the multi-material bar

is rounded so that there is an integer number of node points for each bar component, and

the interfaces occur only at the locations of the nodes.

Since the water in the piston-cylinder assembly is modeled as a linear elastic fluid, which

can again be considered a special case of the SLS model, it is numerically treated in the

same way as the rest of the bar components.

The compatibility equations are implemented along the characteristics to determine ve-

locity (V ), stress (σ), and strain (
) fields for a given initial and boundary conditions. The

nodal points are divided into three categories:

left boundary point (figure 2.6a), interior

points (figure 2.6b), and right boundary point (figure 2.6c).

Figure 2.6: Characteristics solution for (a) Left boundary point, (b) Interior point, and
(c) Right boundary point of the MMHB actuator. The subscript of the spatial element (x)
indicates the location of grid points: 1 and n denote the left and right boundary, respectively,
and i denote an internal point.

25

Interface(x1,tj-1)(x1,tj+1)(x2,tj)(xi-1,tj)(xi,tj-1)(xi+1,tj)(xi,tj+1)(xn,tj-1)(xn-1,tj)(xn,tj+1)Left Boundary Pointσ(x1,t) specifiedInterior PointRight Boundary Pointv(xn,t) = 0(a)(b)(c)For the ith interior node point, define ρ+(i), E+

1 (i) , E+

2 (i), θ+(i) , and C+
V

(i) to be the

corresponding material properties of the material just right of the nodal point. Similarly,
define ρ−(i), E−

(i) to be the corresponding material properties

2 (i), θ−(i) , and C−

1 (i) , E−

V

of the material just left of the nodal point. Notice that:

ρ−(i) = ρ+(i − 1),
E−
1 (i − 1),
1 (i) = E+
E−
2 (i − 1),
2 (i) = E+
θ−(i) = θ+(i − 1),
C−
(i − 1).

(i) = C+
V

V

For the left boundary point (i = 1), only ρ+, E+
only ρ−, E−

2 , θ− , and C−

1 , E−

V

are defined for the right boundary point (i = n). Because

1 , E+

2 , θ+ , and C+
V

are defined, whereas

of the way material properties are assigned to the nodal points, in principle, every interior

nodal point can be treated as an interface. This observation, utilized by [81] in the analysis

of stress waves within an elastic multi-material domain, is central to our treatment of the

more general viscoelastic case. In particular, it recovers the impedance conditions of section

2.2.3.1 for the elastic special case.

For the multi-material bar, the intermediate material interfaces between the bar compo-

nents are treated as bonded interfaces implying no relative motion between the bar compo-

nents. This makes the velocity and stress fields continuous across the interfaces. However,
the strain is generally discontinuous across the interfaces, and we define 
− and 
+ as the

strain just to the left and just to the right of the node, respectively.

The velocities and the stresses are determined at the interior points by applying the

26

compatibility conditions given in equation (2.5) along the characteristics given in equation

(2.4). These left-running and right-running characteristics are shown in figure 2.6b. This

gives:

ρ−(i)C−

V

V (xi, tj+1) =

(i)V (xi+1, tj) + σ(xi+1, tj) − σ(xi−1, tj)
(i)V (xi−1, tj) + ρ+(i)C+
V
ρ−(i)C−
(i) + ρ+(i)C+
V

(cid:26) E−
1 (i + 1)
−(xi+1, tj) − σ(xi+1, tj)

(i)

V

+

ρ−(i)C−

V

∆t

(i) + ρ+(i)C+
V

(i)

σ(xi, tj+1) =

ρ+(i)C+
V

(i)V (xi+1, tj) − ρ−(i)C−
2
(i) − ρ−(i)C−
(i)}V (xi, tj+1)

V

−

{ρ+(i)C+
(cid:26) E−
1 (i + 1)
−(xi+1, tj) − σ(xi+1, tj)

2

V

V

+

θ−(i + 1)

+

∆t
2

θ−(i + 1)
1 (i − 1)
+(xi−1, tj) − σ(xi−1, tj)

(cid:27)

− E+

θ+(i − 1)

(i)V (xi−1, tj)

σ(xi+1, tj) + σ(xi−1, tj)

+

1 (i − 1)
+(xi−1, tj) − σ(xi−1, tj)
E+

2

θ+(i − 1)

(2.8)

(cid:27)

(2.9)

If the interior point is not at an interface, the equations (2.8) and (2.9) reduce to the

corresponding ones given in [1].

The strains at the interior points are determined by applying the compatibility conditions

along the characteristics given in equation (2.6). This gives:


+(xi, tj+1) = 
+(xi, tj−1) +

−

2∆t

1 (i) + E+
E+

2 (i)

σ(xi, tj+1) − σ(xi, tj−1)

(cid:26) E+
1 (i)
+(xi, tj−1) − σ(xi, tj−1)

E+
1 (i) + E+

2 (i)

θ+(i)


−(xi, tj+1) = 
−(xi, tj−1) +
(cid:26) E−
1 (i)
−(xi, tj−1) − σ(xi, tj−1)

σ(xi, tj+1) − σ(xi, tj−1)

2 (i)

2∆t

E−
1 (i) + E−
θ−(i)

E−
1 (i) + E−

2 (i)

−

(cid:27)

(cid:27)

(2.10)

(2.11)

If the interior point is not at an interface, it has 
+ = 
− and the equations (2.10) and (2.11)

27

also reduce to the corresponding equation given in [1].

At the left boundary point (i = 1), stress σ(x1, t) is specified as a boundary condition to

simulate an impact from a striker. Velocity at the left boundary is determined by applying

the compatibility equation (2.5) for the wave propagating in the negative x-direction (figure

2.6a). This gives:

V (x1, tj+1) = V (x2, tj) +

+

∆t

ρ+(1)C+
V

(1)

σ(x2, tj) − σ(x1, tj+1)
(cid:27)

(cid:26) E−
1 (2)
−(x2, tj) − σ(x2, tj)

ρ+(1)C+
V

(1)

θ−(2)

(2.12)

.

The strain calculation at the left boundary point is similar to the internal points.


+(x1, tj+1) = 
+(x1, tj−1) +

σ(x1, tj+1) − σ(x1, tj−1)

(cid:26)E+
1 (1)
+(x1, tj−1) − σ(x1, tj−1)

E+
1 (1) + E+

2 (1)

(cid:27)

(2.13)

−

2∆t

E+
1 (1) + E+

2 (1)

θ+(1)

At the right boundary point (i = n), velocity V (xn, t) is set to zero to simulate a fixed

boundary condition. Stress at the right boundary is determined by applying the compati-

bility equation (2.5) for the wave propagating in the positive x direction (figure 2.6c). This

gives:

(cid:27)
σ(xn, tj+1) = σ(xn−1, tj) + ρ−(n)C−
(cid:26) E+
V (xn, tj+1) − V (xn−1, tj)
1 (n − 1)
+(xn−1, tj) − σ(xn−1, tj)

(cid:26)

(n)

V

+ ∆t

θ+(n − 1)

(cid:27)

(2.14)

28

The strain calculation at the right boundary point is also similar to the internal points.


−(xn, tj+1) = 
−(xn, tj−1) +
(cid:26)E−
1 (n)
−(xn, tj−1) − σ(xn, tj−1)

σ(xn, tj+1) − σ(xn, tj−1)

2 (n)

1 (n) + E−
E−
θ−(n)

−

2∆t

E−
1 (n) + E−

2 (n)

(cid:27)

.

(2.15)

A MATLAB script was developed to implement the aforementioned method of charac-

teristics and numerically model the MMHB actuator. For the numerical modeling results

presented in section 2.3, time-step ∆t was chosen to be 1 µs.

2.2.4 Optimization Methodology

2.2.4.1 Design Parameters

A numerical model for the MMHB actuator dynamics is based on the method of characteris-

tics as elaborated in section 2.2.3.2. A theoretical background review for numerical modeling

is presented in section 2.2.3.1. The design of the MMHB is obtained by coupling this dy-

namical model with an optimization algorithm in order to replicate the target stress profile.

The MMHB actuator design problem involves geometric variables, such as the lengths of

the bar components. The impact between the striker and the multi-material bar gives the

stress boundary condition at the left end of the overall bar assembly. This incident stress

pulse is specified to have a single half-wavelength sinusoidal form similar to the one shown

in figure 2.1b. In practice, such an incident stress pulse is typically generated by an impact

from a spherical striker. To account for various magnitudes of impacts, the incident pulse

amplitude and duration are considered as design variables. The material selection for each

bar component is crucial for introducing interfaces with impedance mismatch in the multi-

material bar, which is the key for transforming the incident pulse into a complicated stress

29

wave. Thus, the material assignment for the bar components are significant design variables.

Table 2.1 gives the details of material selection options we have considered for each of

the bar components in terms of the material properties used to simulate these materials. For

an N-bar assembly of the MMHB actuator, there are 2N+2 design variables: N lengths, N

materials, the incident pulse duration and amplitude. The magnitude ranges for these design

variables are given in table 2.2. Additionally, the optimization algorithm allows the same

materials assignment for two adjacent bar components. This reduces the effective number

of bar components while providing additional design flexibility.

Table 2.1: Material library for multi-material bar components. The material properties for
Material ID 1-5 were taken from [1], and Material ID 6-12 were obtained from the ANSYS
material library [2].

Material ID

Material

E1 (GPa) E2 (GPa) θ(µs) Density (kg/m3)

1
2
3
4
5
6
7
8
9
10
11
12

Epoxy

PMMA 1
PMMA 2
PMMA 3

Polycarbonate
Structural Steel
Stainless Steel

Magnesium Alloy
Gray Cast Iron
Copper Alloy

Aluminum Alloy
Titanium Alloy

3.43
2.94
3.14
3.78
2.3
200
193
45
110
110
71
96

3.43
3.07
3.98
5.24
0.73

0
0
0
0
0
0
0

8.57
95.4
67.4
40.5
140

-
-
-
-
-
-
-

1200
1190
1190
1190
1200
7850
7750
1800
7200
8300
2770
4620

Table 2.2: Design variables for MMHB actuator

Design Variable (Unit)

Minimum Value Maximum Value Step/Resolution

Length of each bar component (mm)
Material ID for each bar component

Incident pulse amplitude (MPa)

Incident pulse width (ms)

50
1
0.1
0.1

2000
12
10
2

1
1
0.1
0.01

Water was modeled as a nearly incompressible linear elastic fluid [82]; and since it was

30

fully confined, its bulk modulus was used as the elastic modulus. To avoid any additional

variables, the length of the water column in the cylinder was fixed to 50 mm for all the

optimization problems. For the numerical model, the bulk modulus of water was taken as

2.15 GPa and the density as 998 kg/m3 [83].

2.2.4.2 Optimization Algorithm

An optimization problem with multiple conflicting objectives has a set of Pareto optimal

solutions instead of a single solution [84]. The Non-dominated Sorting Genetic Algorithm

II (NSGA-II) is an evolutionary optimization algorithm that is widely used to find a set

of optimal solutions for such a multi-objective optimization problem [85]. NSGA-II is a

population-based optimization algorithm. First, the objective functions, lower and upper

bounds for the variables, and the population size for the set of solutions are specified. Then,

a random set of solutions termed as the initial population is created by the algorithm by

assigning random values to the variables. Next, the fitness values are determined for all the

initial solutions based on the objective functions. The algorithm further creates a children

population by performing selection, crossover, and mutation operations on the initial pop-

ulation. Here, the selection operator gives more importance to the solutions with higher

fitness values; the crossover is analogous to reproduction, and the mutation operator mimics

genetic mutation. A survivor operator is then used to create a new population by selecting

the fittest and diverse set of solutions from both the initial and children populations. Lastly,

the initial population is replaced with the new population, and the cycle is repeated for a

specified number of generations. After performing a sufficiently large number of generations,

which is specified in advance, the fittest and diverse solutions from the resultant population

can be considered as optimal solutions. Based on the trade-off between the optimal solutions,

31

the user finally needs to choose the solutions for further implementation. The flow chart for

the genetic algorithm is shown in figure 2.7.

Figure 2.7: Flowchart of the optimization algorithm

2.2.4.3 Objective Functions

For every new design of the MMHB actuator, the output pressure profile inside the pressure

chamber is computed from the dynamic simulation. This output pressure profile is then used

to assess the performance of the actuator design based on two objective functions: (1) the

error between the target pressure profile and the pressure profile generated by the actuator

in the time domain, and (2) the error between these two profiles in the frequency domain.

In order to calculate both of these errors, the time-domain signal S(t) is decomposed into

a sum of sinusoids using a conventional Fourier series expansion [86]

S(t) =

Ai sin (ωit + φi),

(2.16)

l(cid:88)

i=1

32

Initialize Population and EvaluateSpecified number ofgenerations completed ?Selection, Crossover andMutation OperationOutput ResultsNOYESEvaluate and SurvivorOperationwhere Ai, wi and φi are the amplitude, frequency and phase shift of the ith mode. By mini-

mizing the error between the two pressure profiles in the frequency domain, the optimization

algorithm seeks to match the respective modal amplitudes of the two pressure profiles. A

MATLAB script computes the Fast Fourier Transform (FFT) of the pressure profiles and

calculates the root-mean-square error (RMSE) between the profiles in the separate time and

frequency domain

RMSE in time domain =

(cid:115)(cid:80)n

j=1(S(tj)Tar − S(tj)Out)2

and

RMSE in frequency domain =

(cid:115)(cid:80)m

n
k=1(ATar
m

k − AOut

k

(2.17)

)2

.

The superscripts Tar and Out here denote the target and output values of the shown quantity.

For the RMSE calculation in the frequency domain, the frequency range was restricted to

0-10 kHz since all the dominant frequencies fell into this range.

The RMSE in the time domain accounts for the mismatch in both the amplitudes and

phases of the sinusoids constituting the profiles. However, the RMSE in the frequency

domain only accounts for the mismatch in the amplitudes, which emphasizes replicating the

spectral characteristics rather than an exact match of the waveforms in the time domain.

These objectives can be conflicting as the optimized MMHB actuator may not achieve the

best match for amplitudes and the phases simultaneously. Thus, this problem is treated as

a bi-objective optimization problem.

For optimizing the MMHB actuator design, numerous designs need to be simulated and

evaluated for their performance. This simulation process becomes time-intensive for a large

number of designs. Our custom-developed script based on the method of characteristics

offers a significantly faster evaluation of performance than finite elements-based commercial

33

software, and it is easier to integrate with the optimization algorithm because it is also

written in MATLAB. This seamless integration of the simulation and optimization scripts

increases the time-efficiency of the process.

2.3 Results

The methodology to design the MMHB actuator is demonstrated by devising the actuator to

generate various target pressure profiles composed of decaying and non-decaying sinusoids.

Solutions to these optimization problems are discussed in this section.

2.3.1 Case Study-I: Single Frequency Sine Wave

For the first test problem, a sinusoidal wave with 2 kHz frequency and 1 MPa amplitude

was chosen as the target profile. Since the target profile had a single frequency, the number

of bar components was set to 1 for the optimization problem. The NSGA-II optimization

algorithm was employed with a population size of 100, and the number of generations was set

to 100. Thus in total, 10000 MMHB actuator designs were simulated during the optimization

process. This contrasts with the 108 number of possibilities in the overall design space, based

on the variables and ranges detailed in table 2.2. The optimum design was chosen to be the

design with the least RMSE in the frequency domain. The time and the frequency domain

comparison between the target pressure profile and the pressure profile generated by the

optimized MMHB actuator are shown in figures 2.8a and 2.8b. The design details for the

optimized actuator are listed in table 2.3.

For this case study, even though the incident bar consists of a single material, the water

effectively adds an additional interface. Thus, wave reflection, transmission, and interference

34

(a)

(b)

Figure 2.8: (a) Time domain and (b) frequency domain comparison between the target
pressure profile and the output pressure profile generated by the optimized MMHB actuator
for case study-I

35

012345Time (ms)-1.5-1-0.500.511.5Pressure (MPa)Target ProfileOutput Profile0246810Frequency (kHz)00.20.40.60.811.2Pressure Amplitude (MPa)Target ProfileOutput Profileremain in the system. This phenomenon makes it difficult to predict the pressure profile

analytically, and so necessitates numerical simulation. However, if the assembly is simply a

single-material bar without water, as shown in figure 2.1a (Case 1), the stress profile at the

end of the bar can be predicted analytically (shown in figure 2.2a).

2.3.2 Case Study-II: Single-frequency decaying sine wave

For the second test problem, a decaying sinusoidal wave with 2 kHz frequency, 2 MPa

amplitude, and 0.1 decay coefficient was chosen as the target profile.

In particular, the

target profile ST ar(t) was taken as

STar(t) = Ae−2πζf tsin(2πf t),

(2.18)

where A denotes the amplitude, ζ is the decay coefficient, and f is the frequency. Again,

the number of bar components was set to 1 for the optimization problem. The optimization

routine was implemented with a population size of 100 and the number of generations was

set to 200. In total, 20000 MMHB actuator designs were simulated during the optimization

process. The optimum design was chosen to be the design with the least RMSE in the time

domain. The time and the frequency domain comparison between the target pressure profile

and the pressure profile generated by the optimized MMHB actuator are shown in figures

2.9a and 2.9b. The design details for the optimized actuator are listed in table 2.3.

The design optimization problem defined in this case study was also solved for a case

where the number of bar components for the multi-material incident bar was set to 2. This

exercise was conducted to test the robustness of the design optimization methodology. No-

tably, the optimum design obtained from the optimization routine resulted in both bars

36

(a)

(b)

Figure 2.9: (a) Time domain and (b) frequency domain comparison between the target
pressure profile and the output pressure profile generated by the optimized MMHB actuator
for case study-II

37

0123456Time (ms)-1.5-1-0.500.511.52Pressure (MPa)Target ProfileOutput Profile0246810Frequency (kHz)00.050.10.150.20.250.30.35Pressure Amplitude (MPa)Target ProfileOutput Profilebeing assigned the same material. Moreover, their combined length was the same as the

previously obtained single-bar optimum design.

2.3.3 Case Study-III: Sum of two sinusoidal waves

For the third test problem, a sum of two sine waves with 0.5 kHz and 2 kHz frequencies,

and 1 MPa and 0.5 MPa amplitudes was used as the target profile. Since the target profile

contained more frequencies compared to previous target profiles, the number of bar elements

was set to 2 for the optimization problem. The optimization routine was implemented with

a population size of 100 and the number of generations was set to 500. Thus in total,

50000 MMHB actuator designs were simulated during the optimization process. This again

contrasts with the huge number of possibilities in the overall design space. Due to the two

bar components, the possibilities in the overall design space are significantly higher than the

108 possibilities available for the case of a single component MMHB actuator. Figure 2.10

shows the set of solutions obtained from the optimization process. The solutions are plotted

in the objective space, i.e., the axes of the graphs are the objective functions. It is clear that

the objective functions are conflicting, and there is a trade-off between the solutions in terms

of the objective values. The time and the frequency domain comparison between the target

pressure profile and the pressure profile generated by the optimized MMHB actuator designs

are shown in figures 2.11a and 2.11b. The optimized designs chosen for the comparison were

the design with the best match in the frequency domain, the design with the best match

in the time domain, and the design with a moderate match in both the frequency and time

domains. The design details for the optimized actuators are listed in table 2.3.

38

Figure 2.10: The set of optimization solutions obtained for the case study-III. These
solutions are plotted in the objective space where the axes denote the optimization objective
functions

2.3.4 Case Study-IV: Sum of three sinusoidal waves

For the fourth test problem, a sum of three sine waves with frequencies of 0.5 kHz, 1kHz,

and 2 kHz, and amplitudes of 2 MPa, 1.5 MPa, and 1 MPa was used as the target profile.

Since the target profile included a higher number of frequencies, the number of bar elements

was set to 5 for the optimization problem allowing the MMHB actuator to generate more

complex and dynamic pressure profiles. The optimization algorithm was employed with a

population size of 100 and the number of generations was set to 1000. Thus in total, 100000

MMHB actuator designs were simulated during the optimization process. The time and the

frequency domain comparison between the target pressure profile and the pressure profile

generated by the optimized MMHB actuator designs are shown in figures 2.12a and 2.12b.

The optimized designs chosen for the comparison were the design with the best match in

the frequency domain, the design with the best match in the time domain, and the design

39

0.20.250.30.350.40.450.5RMSE in time-domain (MPa)0.0150.020.0250.030.0350.040.045RMSE in frequency-domain (MPa)Optimized designsDesign with best frequency domain matchDesign with best time domain matchDesign with a moderate match(a)

(b)

Figure 2.11: (a) Time domain and (b) frequency domain comparison between the target
pressure profile and the output pressure profile generated by the optimized MMHB actuator
for case study-III

40

0123456Time (ms)-2.5-2-1.5-1-0.500.511.522.5Pressure (MPa)Target ProfileDesign with best frequency domain matchDesign with best time domain matchDesign with a moderate match0246810Frequency (kHz)00.20.40.60.811.2Pressure Amplitude (MPa)Target ProfileDesign with best frequency domain matchDesign with best time domain matchDesign with a moderate match(a)

(b)

Figure 2.12: (a) Time domain and (b) frequency domain comparison between the target
pressure profile and the output pressure profile generated by the optimized MMHB actuator
for case study-IV

41

0123456Time (ms)-4-3-2-1012345Pressure (MPa)Target ProfileDesign with best frequency domain matchDesign with best time domain matchDesign with a moderate match0246810Frequency (kHz)00.20.40.60.811.21.41.61.82Pressure Amplitude (MPa)Target ProfileDesign with best frequency domain matchDesign with best time domain matchDesign with a moderate matchwith a moderate match in both the frequency and time domains. The design details for the

optimized actuators are listed in table 2.3.

For this case study, even if the number of bar components for the multi-material bar was

set to 5, the optimization algorithm was allowed to assign the same materials to adjacent

bars which would effectively reduce the number of bar components. That is what exactly

happened here. All of the optimal designs involved two adjacent bar components with the

same material specification. The resulting optimal designs were effectively 4 component bars.

2.3.5 Case Study-V: Sum of three decaying sinusoidal waves

For the fifth test problem, a sum of three decaying sinusoidal waves with frequencies of 1, 2,

and 3 kHz, amplitudes of 1 MPa each, and decay coefficients 0.1, 0.05, and 0.033 was chosen

as the target profile. Like the previous test problem, the number of bar components was

set to 5 for the optimization problem. The optimization routine parameters were also kept

the same as in the previous test problem. The time and the frequency domain comparison

between the target pressure profile and the pressure profile generated by the optimized

MMHB actuator designs are shown in figures 2.13a and 2.13b. The optimized designs chosen

for the comparison were the design with the best match in the frequency domain, the design

with the best match in the time domain, and the design with a moderate match in both the

frequency and time domains. The design details for the optimized actuators are listed in

table 2.3.

42

(a)

(b)

Figure 2.13: (a) Time domain and (b) frequency domain comparison between the target
pressure profile and the output pressure profile generated by the optimized MMHB actuator
for case study-V

43

0123456Time (ms)-1.5-1-0.500.511.522.5Pressure (MPa)Target ProfileDesign with best frequency domain matchDesign with best time domain matchDesign with a moderate match0246810Frequency (kHz)00.050.10.150.20.250.3Pressure Amplitude (MPa)Target ProfileDesign with best frequency domain matchDesign with best time domain matchDesign with a moderate matchTable 2.3: RMSE values and the optimized actuator design parameters for the MMHB
optimization case studies

Case No. of Bar
Study Components Time domain Frequency domain

RMSE in

RMSE in

I
II
III

IV

V

1
1
2

5

5

(MPa)

(MPa)

0.169
0.091
0.467
0.249
0.337
1.137
0.971
1.036
0.281
0.269
0.270

0.0047
0.0119
0.0191
0.04
0.0265
0.1541
0.1647
0.1578
0.0312
0.0345
0.0328

2.4 Discussion

Length of Bar
Components

Material

IDs

(mm)

90
158

[1062, 103]
[1274, 50]
[1192, 50]

11
4

[11, 10]
[11, 11]
[11, 10]

[1600, 1726, 1826, 1998, 1999]
[1434, 1995, 1819, 2000, 2000]
[1537, 1847, 1786, 2000, 2000]

[637, 141, 107, 1930, 704]
[642, 527, 67, 1686, 704]
[642, 527, 64, 1681, 704]

[8, 11, 6, 8, 8]
[8, 11, 6, 8, 8]
[8, 11, 6, 8, 8]
[5, 10, 11, 12, 5]
[5, 12, 7, 12, 5]
[5, 12, 9, 12, 5]

Pulse

Pulse
Amplitude Width
(ms)

(MPa)

1.3
1.3
2
2
2.1
5.5
4.1
4.6
7.2
5.6
6

0.1
0.13
0.24
0.23
0.22
0.35
0.47
0.43
0.12
0.13
0.13

Special purpose actuators are proposed on the basis of modifying the conventional SHPB

design for the purpose of generating complex loading profiles.

Instead of a conventional

single material incident bar, a multi-material incident bar is employed to alter the dynamic

nature of the incident pulse. The design of the MMHB actuator is specifically optimized to

generate a pressure profile similar to a given target pressure profile. The MMHB actuator

design optimization was conducted for various target profile case studies that included a sum

of decaying and non-decaying sinusoids within the target pressure profile.

For the first two case studies, the target pressure profiles contained single frequency

sinusoids, and the pressure profiles generated by the optimized MMHB actuator were in

good agreement with these target pressure profiles. There wasn’t any significant trade-off

between the optimum solutions obtained on the basis of the time domain comparison vs. a

frequency domain comparison. For the third case study, the target pressure profile was a

sum of two sinusoids. The pressure profile generated by the optimized actuator was able to

fairly replicate the dynamics of the target pressure profile. Lastly, for the case studies IV

and V, the target pressure profiles were composed of three sinusoids. The pressure profiles

44

generated by the optimized actuators for these test problems were not able to capture all

the spectral components of the target pressure profiles. As the spectral complexity of the

target pressure profile increases, the spectral characteristics of the pressure profile generated

by the optimized actuator deviates from that of the target pressure profile.

Presumably, a higher number of bar elements constituting the multi-material incident

bar and a wide range of material selection options will give additional flexibility to the

optimization algorithm to design actuators that can replicate even more complicated loading

profiles. However, as the number of rod components used to construct the multi-material

incident bar increases, the total length of the actuator will increase as well as its complexity

for experimental settings. This increased length of the actuator may limit its applicability

because of the laboratory space or experimental constraints.

A strategy that is not considered in this study is to use a parallel combination of bars that

have been each previously optimized as in the above-explored case studies. That is, a design

employing several multi-material bars connected in parallel with each other could possibly

replicate a multi-frequency pressure profile better than that of a single multi-material bar.

While a strategy like this could have more flexibility than single multi-material bars, its

implementation will bring new challenges that should be the subject of future studies.

For the numerical modeling of the MMHB actuator, water was modeled as a nearly

incompressible linear-elastic fluid; and for simplicity, the cavitation phenomenon associated

with water was not incorporated in the numerical modeling. The cavitation phenomenon is

linked to the formation, expansion, and collapse of the bubbles when a liquid is exposed to

large negative pressures [25,30,87]. Consequently, any cavitation of water in the cylinder may

saturate the negative pressure part of the loading profiles. Future studies should consider a

more specific equation of state for the water constituent that accounts for pressure-density

45

variations and cavitation phenomenon [5].

The presented actuator design can be implemented to generate a stress wave that has user-

defined temporal and spectral characteristics. The actuator allows the users to investigate

the response of the specimen to real world loading conditions in a lab environment. Unlike

the other experimental techniques, this actuator might facilitate the repeatable production

of complex stress wave seen in high-impact loadings.

The MMHB actuator suffers from some limitations that are inherent to the SHPB system.

The MMHB actuator can get bulky as the number of bar components used to build the

multi-material bar increases. Moreover, based on the target loading profile, the MMHB

actuator’s design needs to be optimized, and a new multi-material incident bar needs to be

built. These design changes for every new target loading profile make the implementation

of the MMHB actuator costly and time-intensive, which can significantly limit its usability.

These limitations of the MMHB actuator inspired the design of the experimental apparatus

presented in the next chapter.

46

Chapter 3

Generator of Broadband Pressure

Cycles

3.1

Introduction

This chapter introduces the apparatus created for generating broadband pressure cycles.

The principal idea presented here is implementing a piezoelectric actuator to compress water

confined in piston-cylinder assembly to generate complex and fast-varying pressure profiles.

This chapter is organized as follows. Section 3.2 discusses the working principle and the

fabricated design of the apparatus. A simplistic numerical model analyzing the apparatus

is also discussed in this section. Section 3.3 discusses the system identification methodology

used to characterize the input-output relation of the apparatus. A feedforward controller

designed based on the system identification data to generate user-defined pressure profiles

is also detailed in this section. The results from experiments conducted to generate various

pressure profiles using the apparatus are discussed in section 3.4. Lastly, section 3.5 concludes

the discussion on the experimental apparatus.

47

3.2 Design

3.2.1 Working Principle

Pressure is a surface loading that compresses a specimen from all directions. In this study,

a piston-cylinder assembly filled with water (figure 3.1) was used as a variable-volume con-

finement for generating pressure loading. Water is considered a nearly incompressible fluid

because of its very high bulk modulus of 2.15 GPa [83]; hence, its pressure rises sharply as

a response to relatively small volume changes.

Figure 3.1: A piston-cylinder assembly (l0: initial length of the water volume inside the
cylinder, up: displacement of the piston)

The change in pressure inside the cylinder, ∆P , can be estimated using the Hooke’s law

for volumetric loading [82, 88] given by the equation,

∆P = −K

∆V
V0

,

(3.1)

where K is the bulk modulus for water, V0 is the initial volume of water inside the cylinder,

and ∆V denotes the change in volume of water. Assuming the cylinder walls to be perfectly

rigid and the cylinder having a uniform cross-section, equation (3.1) can be further simplified

48

l0upNearly Incompressible Liquid (Water)Piston with O-ringsCylinderto

∆P = K

up
l0

.

(3.2)

Here, l0 and up denote the initial length of the water volume inside the cylinder and the

displacement of the piston respectively as shown in figure 3.1. As the bulk modulus of water

is extremely high, even a piston displacement of a few micrometers can generate a high

pressure spike inside the cylinder, equivalent to a blast event. A piezoelectric stack actuator

(Model P885.95 Physik Instrumente (PI) GmbH & Co. KG, Karlsruhe, Germany) was used

to drive the piston in the piston-cylinder assembly. The piezoelectric actuator expands when

a voltage signal is applied across its leads. By expanding, the actuator displaces the piston,

compressing the water inside the cylinder and raising the pressure as a result. The time scale

of this actuation is of the order of a few microseconds; this is fast enough to generate pressure

loadings similar to those of a blast event (typically lasting for a few milliseconds). A high-

frequency pressure transducer (Model 113B27 PCB Piezotronics Inc., Depew, NY, USA)

measures the pressure inside the cylinder to give feedback about the actual performance of

the piezoelectric actuator.

3.2.2 Device Fabrication

This section discusses the two design iterations of the pressure apparatus. The limitations

observed in the first-iteration design of the pressure apparatus necessitated a design revision.

3.2.2.1 Design Iteration I

A pressure chamber was machined from an acrylic rod (square cross-section 25.4 mm ×

25.4 mm) by drilling a through-hole of 13.1 mm diameter. A custom end-cap was designed

49

and machined from an acrylic rod (circular cross-section with 25.4 mm diameter) to hold the

pressure transducer as well as seal one end of the chamber. The other end of the chamber

was sealed using a piezoelectric actuator and O-rings. A venting channel, i.e. a circular hole

with 1.6 mm diameter drilled perpendicular to the pressure chamber axis, was incorporated

in the pressure chamber to remove air bubbles and excess water while filling the chamber.

The design of the venting channel is similar to the one implemented in [30]. The schematic

and the fabricated design of the apparatus are shown in figures 3.2a and 3.2b. The assembly

was held together using a fixture built from acrylic plates and steel threaded rods. This

(a)

(b)

Figure 3.2: (a) Cross-sectional schematic view and (b) Fabricated design of the pressure
apparatus (Design iteration I)

50

FixedPressure TransducerO-ringsWaterAcrylic BlockEnd Cap Holding Pressure TransducerO-ringFixedSupport for Piezoelectric Actuator25.4 mm13.1 mm50 mmVenting Channel+-PiezoelectricActuatorfixture was also used to pre-compress the assembly to remove any slack that would prevent

effective compression. The complete setup of the apparatus including the fixture is shown

in figure 3.3.

Figure 3.3: The complete assembly of the pressure apparatus (Design iteration I)

This fabricated design was used to conduct preliminary ex-vivo and in-vitro experiments

involving mouse brain slices and HeLa (Henrietta Lacks) cell cultures. The details of these

preliminary experiments are provided in chapter 6 (sections 6.1.1 and 6.2.1). This fabricated

design was observed to have some design flaws which were limiting the performance of the

apparatus. Firstly, for sealing the actuator side of the pressure chamber, O-rings were placed

between the metal bellows of the piezoelectric actuator and the inner wall of the pressure

chamber. Since the metal bellows were easily deformable, a good reliable seal was not

achieved with this arrangement. This improper seal would result in leaks inhibiting the

pressure generation. Secondly, bending of the acrylic plates was observed when the pressure

chamber and actuator assembly were pre-compressed. Thus, this fixture couldn’t rigidly

support the assembly which directly limited the pressure generation. Lastly, the method

51

Acrylic PlateThreaded Rodused for supporting the piezoelectric actuator strained the lead wires of the actuator which

eventually broke the leads after a series of experiments. All of these design flaws were rectified

in the second design iteration by (a) introducing a separate piston part, (b) using aluminum

plates for the fixture, and (c) designing a new holder for the piezoelectric actuator.

3.2.2.2 Design Iteration II

In the new design, a cylindrical geometry was incorporated for the pressure chamber, piezo-

electric actuator holder, and the end cap which held the pressure transducer. The pressure

chamber was machined from an acrylic rod (circular cross-section with 25.4 mm diameter)

by drilling a through-hole of 13.1 mm diameter. A new custom end-cap was designed and

machined from a medical-grade SAE 316L stainless-steel rod (circular cross-section with

25.4 mm diameter) to hold the pressure transducer as well as seal one end of the chamber.

The other end of the chamber was sealed using a stainless steel piston with O-rings. The

same venting channel design was also incorporated. A custom stainless steel holder was

machined to house the piezoelectric actuator. An acrylic tube (inner diameter 25.4 mm and

outer diameter 31.75 mm) was used to align the pressure chamber and the actuator. The

schematic and the fabricated design of the apparatus are shown in figures 3.4a and 3.4b. A

block with a V-shaped groove to support while ensuring alignment was designed and 3D-

printed out of polylactic acid (PLA) filament. The assembly was held together using a new

fixture built from aluminum plates and steel threaded rods as shown in figure 3.5a. The

complete setup of the apparatus including the fixture is shown in figure 3.5b. This updated

design was used for the rest of analysis discussed in this chapter and to conduct the pre-

liminary in-vitro experiments involving cerebral organoids (discussed in chapter 6, section

6.2.2).

52

(a)

(b)

Figure 3.4: (a) Cross-sectional schematic view and (b) Fabricated design of the pressure
apparatus (Design iteration II)

3.2.3 Instrumentation

The input to the apparatus was a 0−10 V signal. This input signal was amplified to 0−100 V

using a power amplifier (Model E504 Physik Instrumente (PI) GmbH & Co.) and then fed

to the piezoelectric actuator. The input voltage signal determined the displacement of the

piston and altered the pressure in the chamber, which was measured using the pressure trans-

ducer. A signal conditioner (Model 482A21 PCB Piezotronics, Inc.) was used to amplify

the pressure signal measured by the transducer, which was finally recorded using a data ac-

quisition device (Model USB-6341 National Instruments, Austin, TX, USA). The same data

53

FixedPressure TransducerPiezoelectricActuatorPiston with O-ringsWaterAcrylic CylinderEnd Cap Holding Pressure TransducerO-ringFixedHolder for Piezoelectric ActuatorAcrylic Tube25.4 mm13.1 mm50 mmVenting Channel+-(a)

(b)

Figure 3.5: (a) Pressure apparatus fixture consisting of aluminum plates and steel threaded
rods along with the supporting V-block and (b) the complete assembly of the pressure
apparatus (Design iteration II)

54

acquisition device was also used to generate the input voltage signal. A custom LabVIEW

(National Instruments) program was developed to perform and synchronize both the signal

generation and the data acquisition tasks at the sampling frequency of 100 k samples per

second.

3.2.4 Numerical Model

A numerical model of the apparatus was developed to further analyze the pressure gen-

eration mechanism and estimate the theoretical limit of maximum pressure that could be

achieved. The pressure chamber, along with the piston and the end-cap, was modeled as a

one-dimensional domain as shown in figure 3.6. This assembly was simulated by solving the

one-dimensional wave propagation equation,

∂2u

∂t2 = α2 ∂2u
∂x2 ,

(3.3)

which has one dependent variable u and two independent variables t and x [81]. For a

longitudinal stress wave traveling through an elastic medium, u denotes the displacement at

location x and time t. The constant α in equation (3.3) denotes the speed of the longitudinal

Figure 3.6: One-dimensional domain for the numerical model simulating the pressure cham-
ber with white dots showing the locations of the nodes used for extracting pressure profiles

55

x30 mm10 mm10 mmSteel PistonWaterSteel End-capFixeduastress wave traveling in the elastic medium. This wave speed was determined from the elastic

modulus (E) and the density (ρ) of the elastic material using the relation

(cid:115)

α =

E
ρ

.

(3.4)

Water was modeled as a linear elastic fluid [82]; and since it was fully confined, its bulk

modulus was used as the elastic modulus. The steel parts were modeled as linear elastic

solids [88]; and since they experienced uniaxial loading, Young’s modulus was used as the

elastic modulus for steel. The material properties (elastic moduli and densities) used for

water and steel are given in table 3.1. The acoustic impedance (Z) for the wave propagation

in each medium was calculated using the relation

Z = α × ρ.

(3.5)

Table 3.1: Material properties for the numerical model of the pressure apparatus

Material Elastic Modulus (GPa) Density (kg/m3)

Water
Steel

2.15 (Bulk Modulus)

200 (Young’s Modulus)

998
7850

The interfaces between water and steel were modeled as bonded interfaces sharing com-

mon nodes, implying no relative motion at the interfaces. This makes the velocity and

stress fields continuous across the interfaces. The stress wave interaction at the interface

was modeled using the equations:

σr =

ZB − ZA
ZB + ZA

σi

and σt =

2ZB

ZB + ZA

σi,

(3.6)

56

where σi denotes the incident stress wave traveling through material A, which gets reflected

from the interface as σr and transmits into material B as σt [78–80].

Displacement ua(t) was specified as a boundary condition at the left end of the piston

(figure 3.6). This boundary condition physically represents the displacement provided by

the piezoelectric actuator. A fixed boundary condition was specified at the right end of

the end-cap (figure 3.6) to model the support by the aluminum plates and threaded rods

fixture. The method of characteristics was implemented using a MATLAB (MathWorks

Inc., Natick, MA, USA) script to numerically solve the one-dimensional wave propagation

equation (equation (3.3)) for each of these media [81].

As a test case, a sinusoidal profile having 7 µm amplitude and 2000 Hz frequency was

specified as a displacement boundary condition. This boundary condition mimics the dis-

placement provided by the piezoelectric actuator at 2000 Hz operating frequency. The appa-

ratus was simulated by numerically solving the governing equations and the pressure profiles

were extracted for two locations of interest: first being the midpoint of the pressure chamber

and second was the location where the pressure transducer was incorporated in the experi-

mental setup. Figure 3.7 shows both the input displacement profile specified as the boundary

condition and the pressure profiles extracted for the specified locations.

It is noticeable that the pressure profiles have the same form as the actuator displacement

profile. Moreover, the two pressure profiles are almost coinciding, suggesting a minimal

variation in the pressure across different locations inside the chamber. This observation

supports the simplification of a hydrostatic pressure distribution inside the chamber as the
piston speed (≈ 0.05 m/s) is several orders of magnitude smaller than the longitudinal wave
speed in water (≈ 1500 m/s). Thus, the pressure generation can be directly linked to the

piston displacement as stated in equation (3.2).

57

Figure 3.7: Comparison between the actuator displacement profile (left y-axis) specified as
the input boundary condition to the numerical model and the output pressure profiles (right
y-axis) obtained from the numerical simulation

This model of the pressure chamber is a simplistic model that doesn’t account for the

components with non-linear response such as deformation of the O-rings and the chamber,

friction between the O-rings and the pressure chamber walls, etc. Still, it is helpful in esti-

mating the shape and the maximum expected magnitude of the generated pressure profiles.

The effect of unaccounted components in the numerical model on the performance of the

apparatus was further experimentally evaluated through system identification, and a control

strategy was developed to generate complex pressure profiles.

58

00.10.20.30.40.50.60.70.80.91Time (ms)0246810Actuator Displacement (m)0100200300400500600Pressure (kPa)Actuator DisplacementPressure at chamber midpointPressure at transducer location3.3 Control

To generate a user-defined, arbitrarily complex pressure loading cycle in a consistent and

reproducible manner, a controller is needed for the pressure apparatus. Based on the user-

specified target pressure profile, the controller needs to compute appropriate control input,

i.e., an excitation voltage profile driving the piezoelectric actuator. This allows the pressure

profile generated by the apparatus to closely follow the target pressure profile. In control

theory, this control problem is termed as an “output tracking problem” [89]. To design

such a controller for the apparatus, first, the apparatus was characterized by correlating the

inputs given to the apparatus and the measured outputs. This performance characterization

is known as “system identification”.

3.3.1 System Identification

The pressure apparatus is a single-input single-output (SISO) system. The input to the

system is an excitation voltage signal provided to the power amplifier driving the piezoelectric

actuator. The output of the system is the pressure profile generated by the apparatus

in response to the actuation of the piezoelectric actuator. For system identification, the

response of the apparatus to various sinusoidal excitations was measured. A chirp signal,

i.e., a constant amplitude sinusoid with continuously varying frequency (illustrated in figure

3.8), was given as input to the system to evaluate its frequency response. This chirp signal

has a starting frequency fo and end frequency fT , and the instantaneous frequency f is

linearly increased from fo to fT over the total duration of the signal (T ) as given in equation

f (t) = f0 + (

fT − f0

T

) × t.

59

(3.7)

The chirp function in MATLAB was used to generate this signal profile.

Figure 3.8: A chirp signal

To evaluate the frequency response of the system, the system was excited with chirp input

signals having f0 = 100 Hz and fT = 5 kHz for a total duration of 100 ms. This frequency

range was specifically chosen as the dominant frequencies in the ICP profiles observed in a

bTBI event fall within this range [5, 6]. Multiple experiments were performed to evaluate

the system’s response to such chirp inputs of various magnitudes. For each of the chirp

signal, an offset equal to the amplitude of the sinusoids was added to make the chirp signal

non-negative. The input parameters for these sine-sweep experiments are given in table 3.2.

Table 3.2: Input parameters of the sine-sweep system identification experiments conducted
for the pressure apparatus

Start Frequency End Frequency Duration Input Amplitude Number of
experiments

(kHz)

(Hz)

(ms)

(V)

100
100
100
100
100

5
5
5
5
5

100
100
100
100
100

2
4
6
8
10

2
2
2
2
2

As a representative result, the chirp input signal having 6 V amplitude is shown in figure

3.9a, and the corresponding output pressure profile generated by the apparatus is shown

60

TimeAmplitudein figure 3.9b. There was a transport delay between the input excitation and the pressure

output profiles, which corresponds to the response time of the system. This transport delay

was measured using the correlation-based finddelay function in MATLAB and was found

to be 80 µs. This transport delay was removed from further analysis by defining a new

start time for the pressure output profile. It can be observed that even though the input

(a)

(b)

Figure 3.9: (a) Chirp voltage profile given as an input to the pressure apparatus for a
system identification experiment (b) Output pressure profile generated as a response to the
chirp excitation

61

020406080100Time (ms)-101234567Input Excitation (V)020406080100Time (ms)-200-1000100200300400Pressure Output (kPa)signal had a fairly constant amplitude for all the excitation frequencies, the output pressure

profile exhibits amplitude variation with input frequencies. This amplitude variation of the

sinusoids constituting the input and the output profiles was quantified by analyzing the

profiles in the frequency domain. The frequency-domain representation of the input and the

output profiles are shown in figures 3.10a and 3.10b.

(a)

(b)

Figure 3.10: (a) Frequency-domain representation of the chirp input profile (b) Frequency-
domain representation of the output pressure profile

62

012345Frequency (kHz)0.040.060.080.10.120.140.160.18Input Excitation Amplitude (V)012345Frequency (kHz)0246810Output Pressure Amplitude (kPa)The amplitude variation was characterized by computing the gain, i.e., the ratio of output

amplitude to input amplitude for a particular frequency. The DC gain was also determined

by computing the ratio between the output and input amplitudes for the zero frequency. Ad-

ditionally, the phase difference between the sinusoids constituting the input and the output

profiles was computed from the frequency-domain data. The gains and the phase differ-

ences for the sinusoids were determined for all the system identification experiments. These

gains and the phase differences computed from the input-output data for each of the 10

experiments are shown in figure 3.11a and figure 3.11b.

The gain and the phase difference curves determined from the experimental data are

consistent over all the experiments. The gain curves clearly capture the local resonances

present in the system even though the exact factors causing the resonances are hard to

isolate or predict analytically. The average gain and the average phase difference curve are

also shown in the figures 3.11a and 3.11b. Thus, by using these average gain and the average

phase difference curves, the system output can be predicted for a sinusoidal input having a

frequency between 100 Hz to 5 kHz.

By knowing how the amplitudes and phases of sinusoids change from input to output

and modeling the system as a linear time-invariant system [90], the response of the system

can be predicted for an arbitrarily complex input profile composed of the sinusoids of these

frequencies. Finally, the system identification curves, the average gain curve and the average

phase difference curve, were used to design a feedforward controller for the output tracking

problem.

63

(a)

(b)

Figure 3.11: (a) Gain curves and (b) phase difference curves for the pressure apparatus
determined from system identification experiments

64

00.511.522.533.544.55Frequency (kHz)0102030405060708090Gain (kPa/V)Experiment01Experiment02Experiment03Experiment04Experiment05Experiment06Experiment07Experiment08Experiment09Experiment10Average00.511.522.533.544.55Frequency (kHz)-2-1.5-1-0.500.5Phass Difference (rad)Experiment01Experiment02Experiment03Experiment04Experiment05Experiment06Experiment07Experiment08Experiment09Experiment10Average3.3.2 Feedforward Control

For the output tracking problem, the target output profile to be tracked is provided as an

input to the controller. A feedforward controller was designed to compute a control input

profile such that the output generated by the system tracks the target output profile. The

feedforward controller used the system identification curves, the average gain curve and the

average phase difference curve, to compute the control input profile.

First, the target output profile was transformed into its frequency-domain representation

to determine the amplitudes and the phases of the sinusoids of which it is composed. Then,

the amplitudes of the sinusoids for the control input signal were determined by dividing

the output amplitudes with corresponding gains. Also, the phases of the sinusoids for the

control input signal were calculated by incorporating the phase shift from the average phase

difference curve into the phases of corresponding sinusoids in the target output profile. Once

the amplitudes and phases of the sinusoids were determined for the control input signal,

its time-domain representation was reconstructed. This control input signal was then pro-

vided to the system and the actual experimental output profile was measured. Lastly, the

target output profile and the experimental output profile were compared to evaluate the

performance of the feedforward controller. The implementation and the performance of this

controller were tested by tracking a multi-modal pressure profile as detailed in section ??.

3.4 Results

To assess the capabilities and the performance of the apparatus, we designed a series of

experiments to generate various pressure waveforms. Initially, we test the capability of the

device to generate two simple waveforms without using the feedforward control. We demon-

65

strate a single pulse waveform in section 3.4.1, and an approximate Friedlander waveform in

section 3.4.2. Subsequently, we demonstrate the implementation of the feedforward control

to generate a multi-modal waveform in section 3.4.3.

3.4.1 Single Pulse Waveform

In order to generate a single pressure pulse, a single sinusoidal voltage pulse was given as an

input to the apparatus. Three such voltage excitation pulses having different pulse widths

and amplitudes were constructed as detailed in table 3.3. The excitation voltage vs time

graph for these pulses is shown in figure 3.12a. The corresponding pressure pulses generated

as a result of these input voltage pulses are shown in figure 3.12b.

Table 3.3: Input excitation pulses constructed for generating pressure pulses

Excitation Pulse No. Pulse Width (ms) Pulse Amplitude (V)

1
2
3

0.50
1.75
6.00

8
10
7

It can be observed that the pressure pulses closely follow the temporal and magnitude

characteristics of the input voltage pulse. Thus, by changing pulse widths and amplitudes

of the input voltage pulses, pressure pulses with various pulse widths and amplitudes can be

easily generated using this apparatus. Important characteristics of the pressure pulses such

as maximum loading rate (kPa/ms) and total impulse, which is the area under the pressure-

time curve, can also be controlled by manipulating input voltage pulse profiles. Molecular

dynamics simulations suggest parameters like peak pressure, pulse duration, maximum load-

ing rate, and total impulse have a significant effect on cell membrane deformation [91]. Thus,

having an experimental apparatus that could generate various pressure loadings spanning

66

(a)

(b)

Figure 3.12: (a) Sinusoidal excitation pulses given as inputs to the apparatus. (b) Output
pressure pulses generated by the apparatus for corresponding input pulses.

this parameter space facilitates to experimentally test the numerical findings.

The comparison between the pressure profile generated using a modified SHPB setup [8]

and our apparatus is also shown in figure 3.13. The input excitation pulse 2 was designed in

such a way that the pressure pulse generated by the apparatus has a similar pulse width and

magnitude compared to the pressure waveform obtained from the SHPB experiment. Our

apparatus could not only replicate the pressure generation capabilities of an SHPB setup

but also exceed them in terms of generating tunable pressure profiles.

67

012345678Time (ms)0246810Input Excitation (V)Excitation Pulse 1Excitation Pulse 2Excitation Pulse 3012345678Time (ms)0100200300400Pressure Output (kPa)Pressure Pulse 1Pressure Pulse 2Pressure Pulse 3Figure 3.13: Comparison between the pressure pulse generated by our apparatus and the
pressure profile obtained from SHPB experiment [8]

At the end of each excitation pulse, the input to the actuator was set to 0 V as shown in

figure 3.12a. Thus, the excitation input to the actuator was stopped, allowing the actuator

to discharge. The actuator and the piston were allowed to come to rest on their own. Part

of the pressure profiles, since the actuation was stopped, shows this transient response of

the apparatus and eventually pressure decays to zero as seen in figure 3.12b. The amplitude

of this transient pressure was always found to be less than 10 % percent of the pressure

amplitude during the loading cycle.

3.4.2 Friedlander Waveform

A Friedlander waveform consists of an instantaneous rise in pressure followed by an ex-

ponential decay and lastly, a negative pressure phase [16]. The time duration for which

the pressure is above the atmospheric pressure is known as the positive phase duration of

68

012345Time (ms)-50050100150200250300350400450Pressure Output (kPa)Pressure Pulse 2SHPB Experiment Pressure Datathe waveform. As an illustration, a Friedlander waveform having a maximum pressure of

450 kPa and positive phase duration of 10 ms is shown in figure 3.14. In order to generate a

Friedlander waveform, the input voltage signal was constructed to have a sharp rise followed

by an exponential decay. The negative phase feature of the Friedlander waveform was not

incorporated in the input voltage signal because of the actuator limitation that it could only
be operated within a 0− 10 V range to generate compression, not tension. The input voltage

signal starts with an instantaneous rise to the top voltage, followed by a hold time of 1 ms,

and an exponential decay that lasts 10 ms. This input excitation and the resulting output

pressure profiles are shown in figure 3.15. The two pressure profiles were obtained when the

same input voltage profile was fed to the actuator. These two pressure profiles replicate each

other very closely, confirming the repeatability of the apparatus.

Figure 3.14: A Friedlander waveform with 450 kPa maximum pressure and 10 ms positive
phase duration

69

05101520253035404550Time (ms)-1000100200300400500Pressure (kPa)Figure 3.15: Comparison between a Friedlander like input excitation waveform given to
the apparatus (left y-axis) and the output pressure profiles (right y-axis) generated by the
apparatus

Even though the input excitation can have an instantaneous rise, the response time of

the system is finite. The response time of the system depends primarily on the operational

frequency-bandwidth of the piezoelectric actuator. If the input pulse had a sharp peak like

a Friedlander wave, instead of the flat-top shown in figure 3.15, the system would not have

had enough time to respond and would not have generated a pressure peak.

After this initial characterization, we determined that the maximum pressure that could

be generated with the current implementation of our apparatus is 450 kPa, which is limited

by the geometry of the confinement and the maximum displacement range of the piezoelectric

actuator. The maximum operating frequency of this particular actuator is 5 kHz, which also

imposes limits to the shortest rise time for pressure generation and high-frequency spectral

components of the loading cycle.

70

0510152025Time (ms)0246810Input Excitation (V)050100150200250300350400450Pressure Output (kPa)Input ExcitationPressure Output Test 1Pressure Output Test 23.4.3 Multi-modal Waveform

The implementation of the feedforward controller was demonstrated by solving the output

tracking problem for a multi-modal target output profile. As a test case, a multi-modal

pressure profile was constructed by adding three sine waves of frequencies 500 Hz, 1000 Hz,

and 2000 Hz. The amplitudes and the phases of the sinusoids used to construct the target

pressure profile are detailed in table 3.4. Also, offsets equal to the amplitudes of sinusoids

were added to make the target profile non-negative. The multi-modal target pressure profile

constructed for the output tracking problem is shown in figure 3.16.

Table 3.4: Sinusoidal components of multi-modal target pressure profile

Frequency (Hz) Amplitude (kPa) Phase (rad) Offset (kPa)

500
1000
2000

75
75
75

1.5π
1.5π
1.5π

75
75
75

Figure 3.16: Multi-modal target pressure profile constructed for an output tracking prob-
lem

The designed feedforward controller was used to determine the control input profile. The

71

0123456Time (ms)050100150200250300350400Pressure (kPa)Target Pressure Profilecontrol input profile generated by the controller for this output tracking problem is shown in

figure 3.17. It is noticeable that the control input profile exhibits the amplitude and phase

modulation characteristics of the system. This control input profile was then fed to the

amplifier that powered the piezoelectric actuator. The output pressure profile generated by

the apparatus in response to the control input profile was recorded. This experiment was

repeated twice and the measured pressure profiles for both the experiments are shown in

figure 3.18 along with the target pressure profile.

Figure 3.17: Control input profile determined by the feedforward controller

The experimental pressure profiles were found to be in close agreement with the target

pressure profile. The two experimental pressure profiles were almost identical, indicating

the repeatability of the experimental results. The experimental pressure profiles provided

an over 96 % fit to the target pressure profile based on the normalized mean squared error

(NMSE) calculation. Thus, this proposed feedforward control strategy can be used to fairly

track a complex pressure profile composed of the sinusoids having frequencies between 100 Hz

to 5 kHz.

72

0123456Time (ms)0246810Input Excitation (V)Feedforward Control InputFigure 3.18: Comparison between the target pressure profile and the experimental pressure
profiles obtained using the feedforward controller

3.5 Conclusion and Discussion

An experimental apparatus was designed to generate complex pressure profiles by using a

piston-cylinder assembly filled with water. A piezoelectric actuator was used to displace the

piston. The displacement of the piston driven by the actuator was controlled through the

voltage excitation input given to the actuator. This displaced piston further compressed the

water inside the cylinder that generated pressure profiles replicating the dynamics of the

input voltage excitation profiles. To test the capabilities of the apparatus, the piezoelectric

actuator was excited with various input voltage profiles that included single sinusoidal pulses

and a Friedlander-like waveform. The pressure profiles obtained for these input excitation

profiles followed the temporal and magnitude characteristics of their respective input pro-

73

0123456Time (ms)-50050100150200250300350400450500Pressure (kPa)Target Pressure ProfileExperiment 1 Pressure DataExperiment 2 Pressure Datafiles. Moreover, a feedforward controller was designed and implemented, which enabled the

apparatus to generate complex user-defined pressure profiles.

Split-Hopkinson Pressure Bars, shock-tubes, and blast-tubes are relatively complex sys-

tems with inherent hazards such as high pressure, shock, impact, and loud noise. In conse-

quence, their safe operation requires rigorous technical training. Also, they lack portability
due to their typically large footprint (∼ 10 m2 or more) and anchoring requirements (e.g.,

vibration isolation); and are in general not tunable in the sense that their ability to produce

arbitrary loading cycles is very limited. Our apparatus overcomes all of these limitations due
to its very small footprint (∼ 4 × 10−2 m2), no special anchoring requirements, no inherent

safety concerns, technical simplicity, and its ability to produce virtually any arbitrary load-

ing cycle within its pressure and frequency bounds. In addition, our apparatus fits inside

any fume hood due to its compactness, increasing its suitability for biomedical applications.

Depending on the choice of piezoelectric actuator, an apparatus like the one presented herein

is suitable for applications that require resonant frequencies in the range 5 kHz – 16 kHz,

and peak pressures of up to 450 kPa.

Typical ex-vivo and in-vitro experiments on living tissue, organoids, or cell cultures

are performed in nearly incompressible media such as artificial cerebrospinal fluid or other

aqueous solutions [8, 32, 35, 92], enabling our apparatus for studies of the effect of dynamical

loading on these kinds of biological systems.

The inner diameter of the pressure chamber is based on the diameter of the piezoelectric

actuator (11.2 mm). Its inner volume was optimized for the operational characteristics of

the piezoelectric actuator.

In principle, the inner diameter of the confinement could be

increased, but it would require some additional design optimization beyond just a change in

scale to offer similar operational characteristics as with the basic design. The same idea goes

74

for any attempt to increase the maximum pressure, or the maximum achievable frequency:

the resonant frequency, stroke, and maximum force of the piezoelectric actuator define the

limitations to each new design.

The current implementation of this apparatus does not include volumetric expansion

of the confinement chamber, limiting its operation to positive gauge pressures. With the

intention to expand our studies to the effect of cavitation [26,29,30,36,37], we have developed

a rig that includes the ability to exert volumetric expansions and generate cavitation. The

design of this new rig is discussed in the next chapter.

75

Chapter 4

Bench-top Rig for

Controlled-Cavitation Experiments

4.1

Introduction

Cavitation phenomenon is observed when a liquid is subjected to a tensile or negative pres-

sure higher than a threshold, i.e., the liquid’s vapor pressure at that temperature [22]. Cav-

itation encompasses the formation, growth, and collapse of the bubbles caused by the rapid

changes in pressure [23]. The collapse of the cavitation bubbles generates high-pressure

shock-waves that can be damaging to the components in the vicinity of the bubbles [24–26].

Impact forces experienced by the human head are hypothesized to cause cavitation in

cerebrospinal fluid (CSF), which could be damaging to the surrounding brain tissue [7,13,21].

The impact forces can arise because of a sudden impact to the head during sports such as

football, boxing, etc., or an automotive accident, or a blast exposure [36,37,93,94]. A person

exposed to such events could sustain a traumatic brain injury (TBI).

To explore the TBI mechanisms associated with cavitation, we have developed a bench-

top apparatus that can generate controlled cavitation in liquids. This apparatus can be used

to expose living brain tissue or cell culture specimens to a cavitation event in a controlled

laboratory setting. Evaluating the biological response of such specimens could provide insight

76

into TBI injury mechanisms.

The presented apparatus is inexpensive, compact, portable, and highly-controllable com-

pared to other experimental techniques used to generate cavitation. The apparatus uses

a high-speed linear motor to generate negative pressures in a cylindrical cavity, a pressure

transducer to measure the pressure profiles, and a control system consisting of a servo drive

and data acquisition hardware to generate a user-defined characteristic for the cavitation

event.

This chapter is organized as follows. Section 4.2 discusses the working principle and

the fabricated design of the cavitation apparatus. The experiment conducted to produce

cavitation in water, and its results are discussed in section 4.3. Lastly, section 4.4 concludes

the discussion on the cavitation apparatus.

4.2 Design

4.2.1 Working Principle

In this study, a piston-cylinder assembly filled with water (shown in figure 4.1) was used as

variable-volume confinement for generating negative pressure and induce cavitation. Water

has a very high bulk modulus of 2.15 GPa [83]; hence, its pressure changes sharply as a

response to relatively small volume changes.

The change in pressure inside the cylinder, ∆P , can be estimated using the Hooke’s law

for volumetric loading [82, 88] given by the equation,

∆P = −K

∆V
V0

,

77

(4.1)

Figure 4.1: A piston-cylinder assembly filled with water (l0:
volume inside the cylinder, up: displacement of the piston)

initial length of the water

where K is the bulk modulus for water, V0 is the initial volume of water inside the cylinder,

and ∆V denotes the change in volume of water. Assuming the cylinder walls to be perfectly

rigid and the cylinder having a uniform cross-section, equation (4.2) can be further simplified

to

∆P = K

up
l0

.

(4.2)

Here, l0 and up denote the initial length of the water volume inside the cylinder and the

piston’s displacement, respectively, as shown in figure 4.1. As a sign convention, the com-

pressive displacement of the piston is considered positive. Since the bulk modulus of water

is extremely high, even a sub-millimeter piston retraction can generate high negative pres-

sures, enough to induce cavitation in water. Equations (4.1) and (4.2) are valid only until

the negative pressure reaches the cavitation threshold pressure. Any piston retraction be-

yond this threshold would just result in the expansion of the bubbles without any further
decrease in the pressure. Here, the piston speed (≈ 1 m/s) is several orders of magnitude
smaller than the longitudinal wave speed in water (≈ 1500 m/s). Thus, the transient wave

propagation effects can be neglected, and the pressure changes can be directly linked to the

piston displacement as stated in equation (4.2). Consequently, the pressure field is uniform

78

l0upWaterPiston with O-ringsCylinderthroughout the cylinder, and it is measured using the pressure transducer.

A linear motor (Model PS01-37x120F-HP-C LinMot USA Inc.) is used to displace the

piston of the piston-cylinder assembly as shown in figure 4.2. The linear motor is used to

retract the piston from the cylinder, increasing the volume of the cylindrical cavity and

decreasing the pressure as a result. The time scale of this actuation is of the order of a

few milliseconds; this is fast enough to generate negative pressure loadings similar to those

of a blast or a blunt TBI event [36]. A high-frequency pressure transducer (Model 113B27

PCB Piezotronics Inc.) measures the pressure inside the cylinder to give feedback about the

generated negative pressure and the cavitation phenomenon.

Figure 4.2: Experimental setup of the cavitation apparatus

4.2.2 Device Fabrication

The pressure chamber was machined from an acrylic block by drilling a through-hole of

13.1 mm diameter. A custom end-cap was designed and machined from a medical-grade

SAE 316L stainless-steel rod (25.4 mm diameter) to hold the pressure transducer as well as

79

Linear motor (slider)Linear motor (stator)L-shaped clampsHigh-speed camera (lens)End-cap housing pressure transducerPressure chamber filled with waterPistonseal one end of the chamber using an O-ring. The other end of the chamber was sealed

using an aluminum piston with O-rings. A venting channel, which is a 1.6 mm diameter

hole drilled perpendicular to the pressure chamber axis, was incorporated in the pressure

chamber to remove air bubbles and excess water while filling the chamber. The design of the

venting channel is similar to the one implemented in [30]. The piston was connected to the
slider of the linear motor using a threaded rod (M8×1.25 threads). The fabricated design

of the apparatus is shown in figure 4.2. Two custom-machined L-shaped clamps were used

to prevent the motion of the pressure chamber. The entire assembly was bolted down to an

optical table.

4.2.3 Instrumentation and Control

The input to the apparatus was a displacement profile specified for the slider. A power supply

unit (Model S01-72/500 LinMot USA Inc.) was used to power the linear motor along with

a servo drive (Model C1250-LU-XC-0S-000 LinMot USA Inc.). The servo drive was used

to control the slider’s displacement using the pre-configured proportional-integral-derivative

(PID) controller [95]. The gains of the PID controller were tuned for the optimal tracking of

the displacement profile. The slider’s motion determined the displacement of the piston and

altered the pressure in the chamber, which was measured using the pressure transducer. A

signal conditioner (Model 482A21 PCB Piezotronics, Inc.) was used to amplify the pressure

signal measured by the transducer, which was finally recorded using a data acquisition device

(Model USB-6341 National Instruments) at a sampling frequency of 100 kHz. A high-speed

camera (Model Phantom V2512) was also added to the system to record the cavitation events

at 75 k frames per second. Lastly, a 5 V transistor-transistor logic (TTL) trigger pulse was

used to synchronize the actuation of the linear motor with the acquisition of pressure and

80

high-speed imaging data.

4.3 Results and Discussion

To assess the capabilities and the performance of the apparatus, an experiment was con-

ducted to expose water to negative pressures and analyze the cavitation phenomenon. In

order to generate the negative pressures in the piston-cylinder assembly, the piston was re-

tracted as per the specified displacement profile. A sinusoidal half-wavelength pulse with an

amplitude of 5 mm, and pulse-width of 10 ms was specified as a target displacement profile

for the piston. The target displacement profile is shown in figure 4.3 along with the mea-

sured displacement of the piston and the force applied by the motor to generate the piston

displacement.

Figure 4.3: Comparison between the target piston displacement profile and the measured
piston displacement profile for the cavitation experiment (left y-axis). The force applied by
the motor to achieve the target piston displacement is also shown (right y-axis).

81

0102030405060Time (ms)-6-4-20246Piston displacement (mm)-300-200-1000100200300Motor force (N)Target piston displacementMeasured piston displacementMeasured motor forceIt can be observed that the motor force magnitude saturated at 255 N, which was the

maximum force that the motor could generate. This force constraint downgraded the PID

controller’s performance and limited tracking of the target displacement profile. However,

the measured displacement profile had a similar magnitude and timescale compared to the

target profile. Also, additional oscillations were observed in the measured displacement

profile compared to that of the target profile. The motor force saturation and the sharp

pressure variations due to the cavitation phenomenon have contributed to the additional

oscillations observed in the measured displacement profile.

The displacement of the piston affected the pressure inside the cylinder. This pressure was

measured using the pressure transducer placed at the other end of the cylinder. The pressure

profile measured by the transducer is shown in figure 4.4 along with the displacement of the

piston. Some of the key instances for the cavitation event are highlighted on the pressure

profile, and the corresponding snapshots of the pressure chamber are shown in figure 4.5.

It is evident from figure 4.4, as the piston was retracted, the pressure inside the cylinder

rapidly decreased and saturated at -110 kPa. The effect of this negative pressure on the

micro-bubbles present inside the chamber can be seen in figure 4.5 (snapshots A and B). As

the piston was retracted further, there wasn’t any significant change in pressure. However,

the additional retraction of the piston resulted in the further expansion of the bubbles, which

can be seen in figure 4.5 (snapshots C, D, and E). The bubbles reached their maximum size

when the piston was at the fully retracted position. As the piston started to return, the

pressure inside the chamber started to increase. This rise in the pressure resulted in the

collapse of the bubbles, which can be seen in figure 4.5 (snapshots F and G). The collapse

of the bubbles suddenly released all the stored energy in the form of shock-waves. The

presence of the shock-waves was confirmed from both the pressure and high-speed imaging

82

Figure 4.4: Displacement of the piston (left y-axis) and the resulting pressure profile
generated inside the chamber (right y-axis). Some of the key instances of the cavitation
event are highlighted on the pressure profile. The corresponding high-speed images of the
pressure chamber are shown in figure 4.5.

data. The shock-waves generated from the collapse of the bubbles sharply increased the

pressure inside the cylinder. This high-pressure exceeded the rated measurement range of

the pressure transducer, and thus, saturation was observed in the measured pressure profile.

This high-pressure cleared all the bubbles from the pressure chamber, which is evident from

4.5 (snapshot H).

When the pressure inside the chamber saturated at -110 kPa for the first time (time =

2 ms in figure 4.4), the piston was at 0.4 mm retracted position. This implied, the piston dis-

placement of only 0.4 mm was sufficient to generate the cavitation pressure. The subsequent

83

0102030405060Time (ms)-5-4-3-2-1012345Piston displacement (mm)-2000200400600800100012001400Pressure (kPa)ABCDEFGH0 kPa-110 kPatime = 2 msPiston displacementPressureFigure 4.5: High-speed images of the pressure chamber showing the expansion and collapse
of the cavitation bubbles. The pressure inside the chamber at the corresponding instances
is highlighted on the pressure profile shown in figure 4.4.

oscillations observed in the piston displacement profile also crossed the 0.4 mm displacement

threshold. Thus, these piston oscillations generated sufficient negative pressures, which re-

peated the cavitation bubble dynamics.

Few micro-bubbles were already present in the pressure chamber before the piston was

84

ABCDEFGH0 ms2 ms4 ms6 ms8 ms10 ms12 ms13 ms2 mmmoved (can be seen in figure 4.5: snapshot A). The negative pressure loading generated by the

apparatus expanded these existing micro-bubble rather than creating new cavitation bubbles

in the water. Thus, if the focus is to generate new cavitation bubbles in a liquid, deliberate

efforts are needed to remove all the pre-existing bubbles from the liquid. The methods that

can be implemented to remove the pre-existing bubbles are detailed in [25, 26, 30].

The pressure sensor can only measure the pressure variations, and the steady-state pres-

sure is always measured as zero pressure. This limitation of the pressure sensor does not

allow us to measure the reference pressure, which is the pressure inside the chamber before

any piston actuation. However, deliberate efforts were made to keep the reference pressure

of the chamber close to the atmospheric pressure, i.e., the gauge pressure would be close

to zero. For the cavitation event, it was observed that the negative pressure got saturated
at -110 kPa. The vapor pressure for water at 25 ◦C is -98 kPa [96]. At the vapor pressure,

the water vapor and water are in thermodynamic equilibrium [22]. Thus, during the cav-

itation event, it is expected to see the negative pressure saturate at the vapor pressure of

water for that particular temperature. The difference in the measured saturation pressure

and the vapor pressure could be because of the deviation of the reference pressure from the

atmospheric pressure.

4.4 Conclusion

We have designed a bench-top apparatus that can generate negative pressures and induce

cavitation in liquids. This apparatus uses a piston-cylinder assembly that would be filled with

the liquid of interest. A linear motor is used to displace the piston. A PID controller is further

used to control the motor and displace the piston as per the user-specified displacement

85

profile. The retraction of the piston increases the volume of the pressure chamber cavity

and generates tensile or negative pressures. As the negative pressure crosses a cavitation

threshold pressure, it induces cavitation in the liquid present in the pressure chamber. This

apparatus was tested by inducing cavitation in water.

Experimental techniques, such as split Hopkinson pressure bars, shock-tubes, blast-tubes,

drop-towers, etc., currently used to generate negative pressures and study cavitation, are

relatively complex and require rigorous technical training for safe operation. Moreover,

these experimental systems are incredibly bulky and thus, lack portability. Our apparatus

overcomes all of these limitations due to its very small footprint and technical simplicity.

Additionally, the pressure profiles generated by our apparatus can be easily tuned based on

the required characteristics of the cavitation events.

The typical pressure profiles generated by the split Hopkinson pressure bars, shock-

tubes, and blast-tubes follow a Friedlander waveform, which has a large positive pressure

peak followed by a negative pressure phase. It is this negative pressure phase that induces

cavitation. However, the pressure profile generated by our apparatus can be controlled

to exclusively include the negative pressure phase without the initial large positive pressure

peak. This feature allows us to load specimens with negative pressure cycles without exposing

them to initial high positive pressure, which is particularly helpful to decouple and isolate the

damage mechanisms associated with cavitation. The shock-waves arising from the collapse

of the cavitation bubbles also generate a high-pressure, but its magnitude and time-scale

significantly differ from the initial Friedlander-like pressure peak.

86

Chapter 5

Generator of Dynamic Shear Cycles

5.1

Introduction

This chapter introduces the apparatus for generating dynamic shear loading cycles. The

principal idea constituting this apparatus is to use a voice-coil actuator [97] and two parallel

plates assembly to generate shear loadings.

This chapter is organized as follows. Section 5.2 discusses the first design iteration of

the shear apparatus. This section details the working principle, fabricated design, and

performance characterization of the shear apparatus. The feedforward controlled for the

shear apparatus is also described in this section. The limitations of the first-iteration shear

apparatus necessitated a second design iteration. The second design iteration of the shear

apparatus is discussed in section 5.3. The preliminary experiments conducted to determine

the mechanical properties of polyacrylamide (PAA) gelatin specimens, as well as to test the

implementation of the apparatus for loading samples of living tissue and three-dimensional

cell cultures, are discussed in section 5.4. Lastly, section 5.5 concludes the discussion on the

shear apparatus.

87

5.2 Design Iteration I

5.2.1 Working Principle

The shear apparatus consists of two parallel plates with the specimen sandwiched between

these plates. The top plate is connected to a voice-coil actuator (Model NCM02-05-005-

4JBMCS, H2W Technologies). Thus, the movement of the top plate is set by controlling

the motion of the voice-coil actuator. A laser-based displacement sensor (Model LK-H057,

Keyence Corporation) is used to measure the displacement of the top plate. By measuring

the displacement of the top plate, the shear strain (γ) is calculated using the equation,

γ =

utp
d

,

(5.1)

where utp denotes the displacement of the top plate and d denotes the specimen thickness.

The bottom plate is attached to a load cell (MDB2.5, Transducer Techniques) that measures

the force (F ) needed to keep the plate fixed. The shear stress (τ ) experienced by the specimen

is determined using the equation,

τ =

F
A

,

(5.2)

where A denotes the contact area between the specimen and the bottom plate. The schematic

of the experimental apparatus is shown in figure 5.1.

The voice-coil actuator generates a displacement of few millimeters at oscillating frequen-

cies up to 100 Hz. Thus, by controlling the displacement of the voice-coil actuator, dynamic

shear strains of variable magnitudes and frequencies can be induced to the specimen. Strain
rates up to 100 s−1 can be achieved using this apparatus. The mechanical properties of

the specimen can be determined by correlating shear stress and strain data obtained under

88

various dynamic loading conditions.

Figure 5.1: Schematic of the shear apparatus (Design iteration I)

5.2.2 Device Fabrication

Each of the parallel plates consisted of two parts: a base plate and a specimen holder plate.

For the top plate, the base plate was attached to the voice-coil actuator using a set screw, and

for the bottom plate, the base plate was attached to the load cell. Slots were incorporated

in the base plates to house the specimen holder plates. This two-part design of the parallel

plates facilitated the easy placement and removal of the specimen from the apparatus. The

base plates were 3D-printed out of VeroWhite resin, and the specimen holder plates were

3D-printed out of PLA filament. A custom holder was designed and 3D-printed out of PLA

to rigidly mount the voice-coil actuator to a support fixture. The support fixture was built

from T-slotted aluminum rods (80/20 Inc.). The completed setup of the shear apparatus is

shown in figure 5.2.

89

Controller & Data AcquisitionComputerLoad CellVoice Coil ActuatorHigh Speed CameraLaser-based Displacement SensorLinear displacement (utp)Top Plate (movable)Bottom Plate (constrained)SpecimendFigure 5.2: Fabricated setup of the shear apparatus (Design iteration I)

5.2.3 Instrumentation

The input to the apparatus was a −10 to 10 V voltage signal. Based on the input voltage
signal, a proportional −3 to 3 A current signal was generated by a power amplifier (Model

LCAM 5/15, H2W Technologies). Then it was fed to the voice-coil actuator. The input

voltage signal determined the force generated by the voice-coil actuator, which moved the

top plate. The top plate’s displacement was measured using the laser-based displacement

sensor, which was driven by a controller (Model LK-G5001P, Keyence Corporation). The

load cell was used along with a signal conditioner (Model NI9237, National Instruments)

to measure the shear force acting on the bottom plate. A data acquisition device (Model

cDAQ9178 with NI9215 and NI9263, National Instruments) was used to generate the control

input voltage signal and acquire the displacement and force measurement signals. A custom

LabVIEW (National Instruments) program was developed to perform and synchronize both

90

the signal generation and the data acquisition tasks at the sampling frequency of 10 k samples

per second. A high-speed camera (Model Phantom V2512) was also used to capture the

deformation of specimens at 10 k frames per second.

5.2.4 Performance Characterization

The performance of the apparatus was characterized by loading a PAA gelatin specimen

under dynamic shear cycles and calculating the shear strain field from the high-speed images.

These shear strain field measurements were then compared to the shear strain calculated from

the top plate displacement data.

The mechanical properties of these gelatin specimens can be tailored to mimic different

regions of the brain tissue e.g., white and gray matter by varying the weight concentrations
of PAA [98–100]. For these experiments, a cuboidal specimen of size 19 mm × 19 mm ×

6.5 mm was molded from PAA gel having 10 % weight concentration. The PAA gel specimen

is shown in figure 5.3. The main reason for using PAA gelating was its easily tunable material

Figure 5.3: A PAA gelatin specimen

91

properties which could be achieved by varying concentrations of the base acrylic monomers.

Moreover, the PAA gel fabrication is fast, has low energy requirements, and can be done

at room temperature. The fabrication technique detailed in [101] was followed to make the

PAA gelatin specimens.

A sinusoidal voltage signal of 20 Hz frequency was given as input signal to the power

amplifier, which powered the voice-coil actuator with a sinusoidal current signal of amplitude

1 A and frequency 20 Hz. The actuator was powered for 10 cycles lasting over 500 ms. The

current profile fed to the actuator is shown in figure 5.4.

Figure 5.4: Current profile given as an excitation input to the voice-coil actuator

The PAA gelatin specimen was sandwiched between the parallel plates and glued to the

specimen holder plates. Sand particles having a size less than 125 µm were attached to the

surface of the gelatin specimen. Once the voice-coil was actuated, it moved the top plate,

whereas the bottom plate was stationary. The PAA gelatin specimen was sheared between

the plates for 10 cycles at 20 Hz frequency. The deformation of the specimen was quantified

by tracking the groups of sand particles in the acquired high-speed images. The particle

92

0100200300400500Time (ms)-1-0.500.51Current Input Signal (A)image velocimetry (PIV) technique was used to track the groups of sand particles. The PIV

technique tracks a group of particles by cross-correlating consecutive images and computing

the displacement field [102].

The sequence of images representing the initial half of the first shearing cycle is shown in
figure 5.5. The raw images (1280 × 800 pixels) acquired from the high-speed camera had a

spatial resolution of 40.33 pixels/mm. These images were cropped for the region of interest
(541 × 230 pixels), and were pre-processed with histogram equalization operation [103].

Figure 5.5: High-speed images showing the shear deformation of PAA gelatin specimen

93

ABCDEF0 ms5 ms10 ms15 ms20 ms25 ms5 mmThese pre-processed images were cross-correlated using PIV to calculate the displacement

field. For the PIV calculations, a window size of 32 pixels and a 50 % window overlap

was used. The processed images overlaid with the displacement field vectors are shown in

figure 5.6. The PIV calculations were performed using two different methods: our in-house

developed PIV code and the open-source PIVLab toolbox [104]. Both of these methods

produced very consistent results validating the implementation of the PIV technique.

Figure 5.6: Processed images of the gelatin specimen overlaid with displacement field
vectors at time (A) t = 0 ms, (B) t = 5 ms, (C) t = 10 ms, (D) t = 15 ms, (E) t = 20 ms,
and (F) t = 25 ms

Once the displacement field between two consecutive frames was obtained using PIV, the

shear strain field was computed using the equation,

γ =

∂u1
∂y

+

∂u2
∂x

,

94

(5.3)

ABCDEFwhere u1 and u2 denote the displacement components in x and y directions [88]. This shear

strain field was averaged to find the mean shear strain experienced by the specimen. Lastly,

the instantaneous mean shear strain measurements were summed cumulatively to obtain

the time evolution of the shear strain profile. The displacement of the top plate was also

measured using the displacement sensor. The top plate displacement measured for the first

two loading cycles is shown in figure 5.7. The shear strain profile was also calculated using

this displacement measurement and equation (5.1). The shear strain profiles computed from

all these methods for the first two loading cycles are compared in figure 5.8.

Figure 5.7: Measured top plate displacement for the first two loading cycles

It is evident that the strain profiles obtained from all three methods closely follow each

other, showing consistent measurements. Thus, by only measuring the top plate’s displace-

ment, the average shear strain experienced by the specimen can be measured accurately.

Hence, only the displacement data obtained from the displacement sensor was used to com-

pute the shear strains for the later experiments.

The shear strain-rate profile was also computed from the shear strain profile. The com-

95

0102030405060708090100Time (ms)-2-101234Displacement of the Top plate (mm)Figure 5.8: Comparison of the shear strain profiles computed using our in-house developed
PIV script, the PIVLab toolbox, and the measured top plate displacement profile

puted shear strain-rate profile for the first two loading cycles is shown in figure 5.9. The

apparatus was able to induce shear deformations of strains up to 50 % at strain-rates up
to 50 s−1 for 20 Hz operating frequency. Similar experiments were conducted at various

frequencies, and the results are detailed in section 5.4.

Figure 5.9: Shear strain-rate profile for the first two loading cycles

96

0102030405060708090100Time (ms)-0.4-0.200.20.40.60.8Shear strainShear Strain calculated using PIV scriptShear Strain calculated using PIVLab toolboxShear Strain calculated from experimental data0102030405060708090100Time (ms)-80-60-40-200204060Shear strain rate (s-1)5.2.5 Control

5.2.5.1 System Identification

For system identification, the shear apparatus is considered as a single-input single-output
(SISO) system. The input to the system is a −10 to 10 V signal provided to the power

amplifier. Based on the input voltage signal, the power amplifier generates a proportional
−3 to 3 A current signal to drive the voice-coil actuator. Since the shear strain induced

in a specimen depends on the displacement of the top plate (equation (5.1)), the resulting

displacement generated by the voice-coil actuator is considered as the output of the system.

The system identification for the shear apparatus was performed similar to the pressure

apparatus (see section 3.3.1), i.e., by exciting the actuator with a chirp (sine-sweep) signal

and evaluating the displacement generated by the actuator. The voice-coil actuator was

excited with a sine-sweep signal having a 10 Hz start frequency, 200 Hz end frequency, and

the total duration of 2 seconds. Also, the amplitude of the sine-sweep input was varied

over the multiple tests. The input parameters for these sine-sweep experiments are given in

table 5.1. The maximum displacement limit for the voice-coil actuator is 3 mm. When the

actuator was excited with sine-sweep inputs of amplitude greater than 2.5 V, the displacement

profile saturated at 3 mm for some of the frequencies. This output saturation could result in

incorrect system identification. Thus, the sine-sweep amplitude for the system identification

Table 5.1: Input parameters of the sine-sweep system identification experiments conducted
for the shear apparatus (Design iteration I)

Start Frequency End Frequency Duration Input Amplitude Number of
experiments

(Hz)

(V)

(Hz)

(s)

10
10
10

200
200
200

2
2
2

97

1.5
2
2.5

5
5
5

experiments was restricted to 2.5 V.

As a representative result, the chirp input signal having 2 V amplitude is shown in figure

5.10a, and the corresponding output displacement generated by the voice-coil actuator is

shown in figure 5.10b. These input and output profiles were analyzed in the frequency

(a)

(b)

Figure 5.10: (a) Chirp voltage profile given as an input to the shear apparatus for a system
identification experiment. (b) Output displacement profile generated as a response to the
chirp excitation

98

00.511.52Time (s)-2-1012Input Excitation (V)00.511.52Time (s)-3-2-10123Displacement Output (mm)domain to determine the gain curves and the phase difference curves. The methodology

used for this analysis is the same as the one used for the system identification of the pressure

apparatus (see section 3.3.1). The gain curve and the phase difference curve were determined

for all of the 15 sine-sweep experiments. The average gain curves and the average phase

difference curves were computed for the tests having the same input excitation amplitudes.

These curves are shown in figure 5.11a and figure 5.11b along with the overall average gain

curve and the overall average phase difference curve. These system identification curves, the

overall average gain curve and the overall average phase difference curve, were used to design

a feedforward controller for the output tracking problem.

5.2.5.2 Feedforward Control

The design and the implementation of the feedforward controller for the shear apparatus is

the same as the feedforward controller for the pressure apparatus (see section 3.3.2). As a test

case for the output tracking problem, a multi-modal displacement profile was constructed

by adding three sine waves with frequencies of 20 Hz, 40 Hz, and 60 Hz, and amplitude of

0.75 mm. The constructed displacement profile is shown in figure 5.12.

The designed feedforward controller was used to determine the control input profile.

The control input profile generated by the controller for this output tracking problem is

shown in figure 5.13. It is evident that the control input profile exhibits the amplitude and

phase modulation characteristics of the system. This control input profile was then fed to

the power amplifier that powered the voice-coil actuator. The output displacement profile

generated by the voice-coil actuator in response to the control input profile was recorded.

This experiment was repeated twice, and the measured displacement profiles for both the

experiments are shown in figure 5.14 along with the target displacement profile.

99

(a)

(b)

Figure 5.11:
(Design iteration I) determined from system identification experiments

(a) Gain curves and (b) phase difference curves for the shear apparatus

100

020406080100120140160180200Frequency (Hz)00.511.5Gain (mm/V)Excitation Amplitude = 1.5 VExcitation Amplitude = 2 VExcitation Amplitude = 2.5 VAverage020406080100120140160180200Frequency (Hz)-2-1.5-1-0.500.511.5Phass Difference (rad)Excitation Amplitude = 1.5 VExcitation Amplitude = 2 VExcitation Amplitude = 2.5 VAverageFigure 5.12: Multi-modal target displacement profile constructed for an output tracking
problem

Figure 5.13: Control input profile determined by the feedforward controller

The displacement profiles generated by the actuator were in moderate agreement with

the target displacement profile. The two experimental displacement profiles were in good

agreement with each other, indicating the repeatability of the experimental results. The

experimental displacement profiles provided over 89 % fit to the target displacement profile

101

00.020.040.060.080.1Time (s)-2-1012Displacement (mm)Target displacement profile00.020.040.060.080.1Time (s)-8-6-4-20246810Input Excitation (V)Feedforward Control InputFigure 5.14: Comparison between the target displacement profile and the experimental
displacement profiles obtained using the feedforward controller

based on the normalized mean squared error (NMSE) calculation. As seen in figure 5.11a,

the gain curves vary with the input excitation amplitude. Using the average gain curve as a

system identification curve does not capture the variation of the gain with the input ampli-

tude. Thus, system identification misses some of the system characteristics and diminishes

the performance of the feedforward controller.

5.2.6 Limitations

This fabricated design of the first-iteration shear apparatus was used to conduct preliminary

experiments on PAA gelatin specimens. The details of these experiments are provided in

section 5.4. This fabricated design was observed to have some design flaws that necessitated

102

00.020.040.060.080.1Time (s)-2-1.5-1-0.500.511.522.53Displacement (mm)Target Displacement ProfileExperiment 1 Displacement DataExperiment 2 Displacement Dataa revision in the design. There is no provision to precisely control the gap between the two

parallel plates or the specimen thickness in the current design. Moreover, the specimen can-

not be precisely pre-compressed between the parallel plates. Consequently, pre-compression

cannot be used to hold the specimen between the plates, requiring the specimen to be glued

to both plates. The requirement of gluing the specimen to the plates reduces this apparatus’s

applicability in studies involving biological specimens.

In the current design, both of the parallel plates are supported at a single location; the

top is connected to the shaft of the voice-coil actuator, and the bottom plate is connected

to the load cell using threaded rods. This arrangement makes the plates susceptible to

transverse vibrations like a cantilever beam. These transverse vibrations can induce ten-

sion or compression in the specimen and prevent loading it in a simple-shear state. These

transverse vibrations were mainly observed when the PAA gelatin specimens were loaded at

high shearing frequencies. Moreover, the voice-coil actuator’s shaft can rotate along its axis,

which makes the top plate susceptible to an out-of-plane motion, and again, can prevent

loading specimen in a simple-shear state.

The designed feedforward controller for the apparatus relies on the average gain curve

and the average phase difference curve for system characteristics. It was observed that the

gain curve for the voice-coil actuator varies with the input amplitude. This variation was

not fully captured in the average gain curve, leading to inaccurate system identification

and limiting the feedforward controller’s performance. Moreover, the feedforward controller

does not account for the variability in specimen stiffness, which could further degrade its

performance while generating user-specified shear strain profiles.

All of the aforementioned limitations of the shear apparatus design were addressed in the

next design iteration. The details of the revised design are provided in the following section.

103

5.3 Design Iteration II

5.3.1 Working Principle

The second-iteration design of the shear apparatus has the same working principle as the

previous design.

It uses a parallel plate assembly to hold the specimen and a voice-coil

actuator to displace the top plate while the bottom plate is held stationary. The major

upgrade in the second iteration is the new voice-coil actuator (Model LAS16-23-000A-P01-

4E, Sensata Technologies, Inc.) that has an integrated position sensor. Since the position

feedback is available from the voice-coil actuator itself, there is no need for an additional

displacement sensor like the one used in the previous design. Moreover, since the role of this

apparatus is to load biological specimens under dynamic shear loadings and not to determine

the mechanical properties on the specimen, the load cell is removed from the apparatus as

shear stress measurement is not required. Rheometers are better suited and optimized for

determining the mechanical properties of a specimen under dynamic shear loads [59]. Thus,

they should be used instead of this apparatus for material characterization applications. A

translation stage and a micrometer were included in the design, and the voice-coil actuator

was mounted on the translation stage using appropriate fixtures. This arrangement enables

to precisely control the gap between the parallel plates or the specimen thickness. Moreover,

a specimen can be precisely pre-compressed and held between the plates without gluing it

to both plates.

5.3.2 Device Fabrication

The fabricated second-iteration design of the shear apparatus is shown in figures 5.15a and

5.15b. The voice-coil actuator is mounted on the translation stage using custom fixtures

104

(a)

(b)

Figure 5.15: (a) Front-view and (b) top-view of the fabricated setup of the shear apparatus
(Design iteration II)

3D printed out of PLA filament. The top plate is machined out of aluminum and coupled

with the shaft of the voice-coil actuator. Two partially threaded studs are attached to the

105

Parallel-plate assembly witha PAA gel specimenVoice-coil actuator with holderAluminum rail fixtureTranslation stage with a micrometerTop plateBottom plate with holderOptical tablevoice-coil as a guide rail for the top plate. These support studs prevent the out-of-plane

motion of the top plate and also minimize the transverse vibrations of the top plate. This

arrangement ensures that the top plate is always parallel to the bottom plate, and the gap

between the plates is constant throughout the actuation cycle. The bottom plate is held

fixed using a custom 3D printed fixture. The voice-coil and the translation stage assembly

is mounted on a T-slotted aluminum rail (80/20 Inc.), bolted to an optical table. Similarly,

the bottom plate fixture is mounted on an aluminum rail (80/20 Inc.), which is again bolted

to the optical table.

5.3.3 Instrumentation and Control

A servo drive (Pluto Servo Drive, Ingenia Motion Control) is used to power as well as control

the voice-coil actuator. The servo drive uses position feedback from the integrated position

sensor of the voice-coil actuator to control it a closed-loop. A proportional-integral-derivative

(PID) controller is implemented using the servo drive to control the displacement of the voice-

coil actuator. The gains of the PID controller were tuned for the optimal tracking of the

user-defined displacement profiles. A 24 V power supply unit (Model MS2-H50, Keyence

Corporation) is used to power the servo drive. The target displacement profile for the
voice-coil actuator is provided in the form of a −10 V to 10 V signal, which corresponds to
−3 mm to 3 mm displacement. A data acquisition device (Model cDAQ9178 with NI9215

and NI9263, National Instruments) was used to generate the target displacement voltage

signal and acquire the actual displacement data measured by the integrated position sensor.

A custom LabVIEW program was developed to perform and synchronize both the signal

generation and the data acquisition tasks at the sampling frequency of 10 k samples per

second. A high-speed camera (Model Phantom V2512) was also added to the system to

106

record the shear loading cycles at 10 k frames per second. Lastly, a 5 V transistor-transistor

logic (TTL) trigger pulse was used to synchronize the actuation of the voice-coil actuator

with the acquisition of displacement and high-speed imaging data.

The performance of the closed-loop PID controller was determined by conducting shear

tests on a PAA gelatin specimen. For these tests, a PAA gelatin specimen of 10 % concen-

tration (weight/volume) and a disk shape with 25.4 mm diameter and 4 mm thickness was

used. The specimen was loaded at three frequencies of 10 Hz, 50 Hz, and 100 Hz, for 10 load-

ing cycles each. For all of these loading cycles, the shear strain amplitude was set to 10 %.

Based on the loading cycles, the target displacement profiles for the voice-coil actuator were

determined and fed to the controller. These target displacement profiles are shown in figure

5.16 along with the actual displacement profiles generated by the actuator. It is evident that

the PID controller could reasonably track the target displacement profile. There is a delay

between the target profile and the measured displacement profile. The inertia of the system

and the response time of the PID controller give rise to this delay. The delay between the

profiles become prominent for the higher frequencies.

The shear strain profiles were determined for the experimental displacement profiles using

the equation (5.1). Since the delay between the target and the measured displacement profiles

does not affect the shear loading experienced by the specimen, the delay was removed while

plotting the strain profiles by defining a new starting time for each test. The shear profiles

are shown in figure 5.17. These profiles exhibit a sinusoidal form with a 10 % shear strain

amplitude. Based on the results, the PID controller’s performance was considered reasonable

for generating sinusoidal shear loading cycles. This second-iteration fabricated design of

the shear apparatus was used to conduct preliminary experiments on mouse brain tissue

specimens. These ex-vivo experiments are discussed in chapter 6 (section 6.1.2).

107

(a)

(b)

(c)

Figure 5.16: Comparison between the target displacement profiles and the experimental
displacement profiles generated by the voice-coil actuator for a loading frequency of (a) 10 Hz,
(b) 50 Hz, and (c) 100 Hz

108

00.10.20.30.40.50.60.70.80.91Time (s)-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5Displacement  (mm)Target profileExperimental data00.020.040.060.080.10.120.140.160.180.2Time (s)-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5Displacement  (mm)Target profileExperimental data00.010.020.030.040.050.060.070.080.090.1Time (s)-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5Displacement  (mm)Target profileExperimental data(a)

(b)

(c)

Figure 5.17: Sinusoidal shear strain profiles generated by the shear apparatus (Design
iteration II) for a loading frequency of (a) 10 Hz, (b) 50 Hz, and (c) 100 Hz

109

00.10.20.30.40.50.60.70.80.91Time (s)-0.15-0.1-0.0500.050.1Shear strain00.020.040.060.080.10.120.140.160.180.2Time (s)-0.15-0.1-0.0500.050.1Shear strain00.010.020.030.040.050.060.070.080.090.1Time (s)-0.15-0.1-0.0500.050.1Shear strain5.4 Results and Discussion

This section presents the results from the material characterization experiments performed

on PAA gelatin specimens. The shear apparatus fabricated in the first design iteration

was used for these shear experiments. In these tests, the PAA gelatin specimens of 10 %

weight by volume concentration were characterized under dynamic shear loading. Each of

the loading cycles consisted of a sinusoidal shear loading of a particular frequency for 10

oscillations. The loading frequency was varied from 10 Hz to 100 Hz in steps of 10 Hz. The

shear test at each frequency was repeated 4 times. Thus, 40 shear tests were performed to

determine the shear properties of PAA gelatin specimens in total. For each test, the shear

strain profile was obtained by measuring the top plate displacement, and the shear stress

profile was computed from the force data measured with the load cell. The representative

shear stress-strain profiles for 20 Hz and 50 Hz loading cycles are shown in figures 5.18a and

5.18b. The shear stress-strain profiles follow a similar temporal character. Also, the stress

profile lags with respect to the strain profile.

For each of these experiments, the current profiles given as excitation inputs to the

voice-coil actuator were sine waves of particular operating frequencies. However, the strain

profiles generated by the apparatus are not perfectly sinusoidal. The system dynamics of

the apparatus, which includes the response of voice-coil to input current, hysteresis, out of

the plane motion or the transverse oscillations of the parallel plates, etc., is responsible for

the non-sinusoidal strain profiles. These limitations of the apparatus were addressed in the

second design iteration.

The shear stress-strain profiles for each experiment were correlated to evaluate the shear

moduli. The amplitude of the shear strain profile (γa) and the amplitude of the shear stress

110

(a)

(b)

Figure 5.18: Shear strain (left y-axis) and shear stress (right y-axis) profiles for PAA
gelatin shear tests at (a) 20 Hz and (b) 50 Hz loading frequency

profile (τa) for the actuation frequency were determined by transforming the profiles in the

frequency domain using the Fast Fourier Transform (FFT). These amplitudes were then

111

00.050.10.150.20.250.30.350.40.450.5Time (s)-0.1-0.0500.050.10.15Shear Strain-1-0.500.511.5Shear Stress (kPa)Shear StrainShear Stress00.020.040.060.080.10.120.140.160.180.2Time (s)-0.15-0.1-0.0500.050.10.150.2Shear Strain-1.5-1-0.500.511.52Shear Stress (kPa)Shear StrainShear Stressused to determine the complex shear modulus (G∗) of the gelatin specimen [60] using the

equation:

G∗ =

τa
γa

.

(5.4)

The complex shear moduli were determined for all the operating frequencies, and the varia-

tion is shown in figure 5.19.

Figure 5.19: Variation of the complex shear modulus of PAA gelatin with loading frequency

The time delay (∆ttotal) between the shear strain and stress profiles consists of two parts:

the transport delay (∆ttransport) and the delay due to the viscoelastic nature of the PAA

gel (∆tviscoelastic). This gives:

∆ttotal = ∆ttransport + ∆tviscoelastic

(5.5)

Since the strain is measured at the top plate, and the stress is measure at the bottom plate,

112

0102030405060708090100110Frequency (Hz)02468101214161820Complex Shear Modulus (kPa)there is a transport delay that corresponds to the time taken by the shear waves to travel

across the specimen. This transport delay was estimated from the elastic shear wave speed

(β) and the specimen thickness (d) using the equation,

∆ttransport =

d
β

.

The elastic shear wave speed is given by the equation [60, 81],

(cid:115)

β =

G∗cos(δ)

ρgel

,

(5.6)

(5.7)

where ρgel denotes the density of the gel specimen and δ is the phase difference between the

shear stress-strain profiles arising from the viscoelastic nature of the specimen. The density

was determined to be 1200 kg/m3. Also, ∆tviscoelastic is related to δ as given in the equation,

∆tviscoelastic =

δ

2πfa

,

(5.8)

where fa is the voice-coil actuation frequency. Combining equations (5.5) to (5.8) gives:

(cid:114) ρgel

∆ttotal = d ×

G∗cos(δ)

+

δ

2πfa

.

(5.9)

The time delay (∆ttotal) between the shear strain and stress profiles was measured by

cross-correlating them using finddelay function in MATLAB. Since all the entities in equation

(5.9) except δ are predetermined, this non-linear equation was solved numerically using

MATLAB to evaluate δ. Consequenctly, the loss factor for the gelatin specimen, i.e., tan(δ),

was determined [77]. The variation of the loss factor with frequency is shown in figure 5.20.

113

Figure 5.20: Variation of the loss factor of PAA gelatin with loading frequency

The loss factor was observed to increase with the increase in loading frequency. Theoret-

ically, since the specimen is in a solid-state at resting condition, the loss factor should always

lie between 0, and 1 [59]. The experimentally evaluated loss factor was consistent with the
theoretical limits. The phase lag δ was further used to determine the storage modulus (G(cid:48)),
and the loss modulus (G(cid:48)(cid:48)) based on the relations [60, 77],

G(cid:48) = G∗cos(δ) and G(cid:48)(cid:48) = G∗sin(δ).

The loss factor represents the ratio of the loss modulus to the storage modulus. The variation

of storage and loss moduli with the loading frequency is shown in figure 5.21.

For the range of loading frequencies analyzed, it can be observed that the loss modulus

increased with the increase in loading frequency. This signifies that more energy is dissipated

114

0102030405060708090100110Frequency (Hz)00.10.20.30.40.50.60.70.80.91Loss FactorFigure 5.21: Variation of the storage and loss modulus of PAA gelatin with loading
frequency

under dynamic shear loading at higher loading frequencies. However, this observation is

limited to the range of frequencies analyzed in these experiments.

As the actuation frequency increases, the relative contribution of the specimen’s inertia

to the force measured at the bottom plate increases. These inertial effects lead to erroneous

measurement of the shear force and, consequently, the shear moduli. Ideally, the contribution

of the inertial effects to the shear force measurement may be neglected if the specimen

thickness is small compared to the elastic shear wavelength (λ) in the specimen [60, 105].

The ratio of specimen thickness to the elastic shear wavelength is given by the equation,

(cid:114) ρgel

G(cid:48) .

= d × fa ×

d
λ

(5.10)

For the performed shear tests, the d/λ ratio was less than 0.1 for actuation frequency up to

115

0102030405060708090100110Frequency (Hz)0246810121416Shear Modulus (kPa)Storage ModulusLoss Modulus50 Hz and increased steadily with the actuation frequency. The d/λ ratio was 0.17 for 100 Hz

actuation frequency. Thus, at higher actuation frequencies, the inertial effects contributed

to the erroneous determination of the shear moduli. The sudden rise in the complex shear

modulus (seen in figure 5.19) at high actuation frequencies results from the overestimation

of the modulus due to the inertial effects. The d/λ ratio can be reduced by decreasing

the specimen thickness. Thus, to reduce any error associated with the inertial effects, the

specimen thickness should be chosen appropriately.

5.5 Conclusion

We have developed an experimental apparatus that can induce shear loadings of various

strain amplitudes and frequencies at the levels relevant to bTBI. This apparatus uses a

voice-coil actuator and two parallel plates assembly to generate dynamic shear loadings. In

this chapter, we have discussed the two design iterations of the shear apparatus. The first-

iteration shear apparatus was used to conduct material characterization experiments for PAA

gelatin specimens. These experiments also tested the implementation of the apparatus for

loading samples of living tissue and three-dimensional cell cultures. The insights from these

material characterization experiments inspired the design improvements implemented in the

second iteration of the shear apparatus. Furthermore, the second-iteration shear apparatus

was used to perform bTBI-related ex-vivo experiments on mouse brain tissue specimens.

Researchers have previously used a modified split Hopkinson pressure bar (SHPB) that

includes a parallel plate fixture to shear the specimen at high strain-rates [57,58]. Usually, the

high shear strain magnitudes generated by SHPB result in specimen tearing. These strain

magnitudes are significantly higher than the ones associated with a primary bTBI event.

116

Moreover, the SHPB technique cannot generate oscillatory shear deformation associated with

traveling shear waves. Rotation rheometers [59], and the other custom shear instruments

developed by researchers [60–63] can generate either the shear strain amplitudes or loading

frequencies typically associated with a bTBI event, but not both simultaneously. Our shear

apparatus (second-iteration design) overcomes these limitations and can induce shear loads

at the levels relevant to bTBI. Moreover, the designed apparatus is inexpensive, compact,

portable, and highly controllable, making it well suited for biomedical applications.

117

Chapter 6

Preliminary Ex-vivo and In-vitro

Experiments

To explore the response of biological specimens to complex time-varying pressure cycles and

dynamic shear loadings, we performed few preliminary ex-vivo and in-vitro experiments.

These experiments served a dual purpose: (1) determining the pathophysiological response

of the animal brain tissue specimens, cell cultures, and cerebral organoids to the blast-like

intracranial effects and (2) testing the implementation of our experimental devices for bTBI-

related studies.

In this chapter, section 6.1 discusses the preliminary ex-vivo experiments involving the

mouse brain tissue specimens. The in-vitro experiments involving cell cultures and cerebral

organoids are discussed in section 6.2.

6.1 Ex-vivo Experiments

This section discusses the preliminary ex-vivo experiments involving the mouse brain tissue

specimens. The experiments were performed in collaboration with the Cox Lab in the

Department of Physiology at Michigan State University. The Institutional Animal Care

and Use Committee at Michigan State University approved all animal procedures performed

here. The tests performed to evaluate the electrophysiological response of neurons to time-

118

varying pressure cycles and dynamic shear loadings are discussed in sections 6.1.1 and 6.1.2,

respectively. The specimen preparation and post-exposure analysis for these tests were

performed by Dr. Joseph Beatty at the Cox Lab.

6.1.1 Ex-vivo Experiments involving Time-varying Pressure Cy-

cles

6.1.1.1 Experimental Procedure

Specimen Preparation: Mouse coronal brain slices containing the primary somatosen-

sory cortex [106] were obtained by slicing the brain shortly after decapitation. The obtained

slices were 300 µm thick. These brain slices were kept submerged in a cold, oxygenated

(95% O2 and 5% CO2) artificial cerebrospinal fluid (aCSF) containing (in mM): 2.5 KCl,

26 NaHCO3, 1.25 NaH2PO4, 10.0 MgCl2, 2.0 CaCl2, 234.0 sucrose, and 11.0 glucose. The

brain slices submerged in aCSF are shown in figure 6.1.

Figure 6.1: Mouse coronal brain slices in aCSF

119

Loading Protocol:

In this study, our first-iteration prototype of the pressure apparatus

(discussed in chapter 3) was used to expose the brain slices to a fast-varying pressure profile,

similar to an intracranial pressure profile induced by a blast exposure [5, 6, 29, 107]. After

obtaining all the brain slices, the slices were transferred to the pressure chamber of the

apparatus (shown in figure 6.2).

Figure 6.2: Pressure chamber (Design iteration I) containing the mouse brain slices in an
aCSF environment

The same aCSF was also used as the pressure chamber fluid. The pressure chamber

filled with the aCSF was then closed and sealed using the end-cap. The actuator and the

pressure chamber assembly was then held together and pre-compressed using the fixture

made of acrylic plates and threaded rods. The experimental apparatus containing the mouse

brain slices is shown in figure 6.3. After pre-compression, the brain slices were exposed to

sinusoidal pressure loading of 250 kPa amplitude and 3 kHz frequency for a 7 ms duration.

The measured pressure loading profile is shown in figure 6.4. This pressure profile has a

similar magnitude and dynamics compared to the intracranial pressure profiles observed

during a bTBI event.

120

Figure 6.3: Pressure apparatus (Design iteration I) with the mouse brain slices placed
inside the chamber

Figure 6.4: Pressure loading profile exerted onto the mouse brain slices

The brain slices were freely floating inside the pressure chamber. As the actuation speed
of the piezoelectric actuator (≈ 0.1 m/s) is several orders of magnitude smaller than the
longitudinal wave speed in water (≈ 1500 m/s), the transient wave propagation effects can

be neglected. The pressure generated inside the chamber can be considered as hydrostatic

pressure. Thus, irrespective of the location and the orientation of the brain slices, each brain

slice was exposed to the same pressure loading.

121

012345678910Time (ms)-100-50050100150200250300Pressure (kPa)The duration of each experiment was minimized as much as possible to reduce the stop-

page of oxygen flow to the aCSF medium. The duration for each test was around 3 minutes.

A total of 4 tests were performed, including a control trial where the brain slices were kept

in the pressure chamber but not exposed to the pressure loading.

6.1.1.2 Results

Following the trauma, the slices were incubated for an hour prior to electrophysiological

recordings. The procedures followed for the electrophysiological recordings are detailed

in [108, 109]. In these recordings, the neurons were excited with current pulses and their

response was recorded through the voltage measurements. A representative result obtained

from these recordings is shown in figure 6.5.

Figure 6.5: Electrophysiological recordings of the mouse brain slices

Initial electrophysiological recordings from layer 5 and layer 2-3 pyramidal neurons re-

vealed depolarized resting membrane potentials (-43 mV compared to -85 mV for neurons

from the control or naive slices), decreased input resistance (calculated from the linear slope

of the voltage-current relationship), and short broad action potentials (extremely small pulse-

widths associated with the neurons’ firings). It was noted that recordings one to two hours

122

after the trauma event appeared to recover to “normal” conditions. These normal conditions

were similar to those slices that were either from naive animals or brain slices placed in the

pressure chamber but not exposed to the pressure wave (control slices). Our early findings

suggest an early trauma to the excitability of the neurons; however, these changes do not

appear to be long-lasting.

6.1.2 Ex-vivo Experiments involving Dynamic Shear Loadings

6.1.2.1 Experimental Procedure

Specimen Preparation: B6 mice of either sex were used in these experiments (postnatal

age: 25 to 38 days). Animals were deeply anaesthetized with 3% isoflurane and decapitated.
The brains were removed and placed into cold (<4 ◦C), oxygenated (95% O2 and 5% CO2)

aCSF physiological solution containing (in mM): 126.0 NaCl, 26.0 NaHCO3, 2.5 KCl, 1.25

NaH2PO4, 2.0 MgCl2, 2.0 CaCl2, and 10.0 glucose. Brains were blocked coronally and glued

caudal side down on the bottom plate of the shear apparatus (Design iteration II).

Loading Protocol: The second-iteration prototype of our shear apparatus (discussed in

chapter 5) was used to exert dynamic shear loads on the brain tissue specimens. The spec-

imens were glued to the bottom-plate of the shear apparatus. The top-plate was precisely

lowered using the translation stage and the micrometer to apply a 10% uniaxial compression

onto the specimens, which ensured the contact between the top-plate and the specimens

throughout the shear loading cycles. The stickiness of the brain tissue specimens and the

applied pre-compression were sufficient to prevent the slippage between the top-plate and

the specimens, eliminating the need of gluing the specimen to both the plates, which sig-

nificantly simplified the tissue handling process. After pre-compression, the specimens were

123

loaded with sinusoidal shear strain profile of 50 Hz frequency for a duration of 1 s (total 50

loading cycles). The specimens were loaded at two amplitude levels: 10% and 20% shear

strain. These dynamic shear loadings mimic the low strain amplitudes and high strain-rates

shear deformation associated with traveling intracranial shear waves [10]. For the control

tests (0% shear strain amplitude), the specimens were pre-compressed but not exposed to

any shear loadings. The measured shear strain profiles for the 10% and 20% shear amplitude

tests are shown in figure 6.6 (the first 10 out of the total 50 loading cycles are shown).

Figure 6.6: The shear strain profiles (loading frequency: 50 Hz; shear strain amplitudes:
10% and 20%) exerted onto the brain tissue specimens (the first 10 out of the total 50 loading
cycles are shown)

A high-speed camera (Model Phantom V2512) was also used to capture brain specimens’

deformation at 4000 frames per second. The high-speed images were used to visually verify

the simple-shear deformation of the brain specimens and examine if there is any slippage

124

00.020.040.060.080.10.120.140.160.180.2Time (s)-0.25-0.2-0.15-0.1-0.0500.050.10.150.20.25Shear strainTarget shear strain amplitude = 10%Target shear strain amplitude = 20%between the parallel plates and the brain tissue specimens. However, high-speed images were

not used for any strain measurements as the top-plate displacement data was sufficient to

determine the shear stains. As an illustration, the undeformed state (0% shear strain) of

the brain tissue specimen and the deformed state with 20% shear strain are shown in figures

6.7a and 6.7b, respectively.

(a)

(b)

Figure 6.7: High-speed images of the brain tissue specimen under dynamic shear loads:
(a) undeformed state (0% shear strain) and (b) deformed state (20% shear strain)

Brain Slice Preparation: Following the shear loads, the brains were placed in cold
(<4 ◦C), oxygenated (95% O2 and 5% CO2) slicing solution containing (in mM): 2.5 KCl, 26

NaHCO3, 1.25 NaH2PO4, 10.0 MgCl2, 2.0 CaCl2, 234.0 sucrose, and 11.0 glucose. Coronal

slices (300 µm thickness) at the level of the primary somatosensory cortex [106] were cut

using a vibrating tissue slicer. The slices were incubated in oxygenated (95% O2 and 5%

CO2) physiological solution containing (in mM): 126.0 NaCl, 26.0 NaHCO3, 2.5 KCl, 1.25
NaH2PO4, 2.0 MgCl2, 2.0 CaCl2, and 10.0 glucose at 36 ◦C for at least 30 minutes, then kept

at room temperature for the remainder of the experiment. Electrophysiological recordings

were performed an hour following decapitation in a recording chamber that was maintained
at 32 ◦C, and slices were continuously superfused (2.5 ml/min) with physiological solution.

125

Electrophysiology: Whole-cell intracellular recordings were obtained from layer 5 pyra-

midal neurons of the primary somatosensory cortex (barrel cortex). The procedures followed

for the electrophysiological recordings are detailed in [108, 109]. Recordings pipettes were

pulled from borosilicate glass and had tip resistances of 4-6 MΩ when filled with a solu-

tion containing (in mM) 117.0 K-gluconate, 13.0 KCl, 1.0 MgCl2, 0.07 CaCl2, 0.1 EGTA,

10.0 HEPES, 2.0 Na-ATP, and 0.4 Na-GTP. The pH and osmolarity of the internal solution

were adjusted to 7.3 and 290 mOsm, respectively. Recordings were obtained using a Multi-

clamp 700B amplifier (Molecular Devices, Sunnyvale, CA). An 8 mV junction potential was

corrected for in all voltage measurements. Signals were digitized at 10 kHz and stored for

subsequent analyses using pCLAMP software (Molecular Devices).

Input resistances were calculated from the linear slope of the voltage-current relationship
obtained by applying constant current pulses ranging from −40 to +40 pA (1 s duration).
The time constant measurements were obtained from a small hyperpolarizing (−10 pA) short

duration (100 ms) current pulse. Action potential output for a single neuron was obtained

by applying constant current pulse ranging from 0 pA to 1600 pA by 40 or 80 pA intervals for

1 second. Action potential frequency versus intensity curves (F-I curves) were then generated
from the collected data. Summarized data is expressed as mean ± standard deviation. P-

values less than 0.05 were considered statistically significant.

6.1.2.2 Results

Intrinsic Properties: Whole-cell current-clamp recordings were used to measure the three

intrinsic neuronal properties: (1) resting membrane potential, (2) input resistance, and (3)

time constant. These three intrinsic properties were used as indicators of cell health following

the three different shear strain conditions: (a) 0% shear strain (control), (b) 10% shear strain,

126

and (c) 20% shear strain.

Resting membrane potential measurements give insight into the neurons’ level of exci-

tation. Resting membrane potential measurements were not significantly different between

the shear strain conditions tested (one-way ANOVA, p>0.05). Neurons exposed to 0% shear
strain possessed an average resting membrane potential of −75.5 ± 4.1 mV (5 cells from 1

animal, shown in figure 6.8 (Black)). While neurons exposed to 10% shear strain possessed
an average resting membrane potential of −69.0 ± 11.1 mV (8 cells from 1 animal, shown

in figure 6.8 (Blue)) and neurons exposed to 20% shear strain possessed an average resting
membrane potential of −71.0 ± 12.5 mV (10 cells from 2 animals, shown in figure 6.8 (Red)).

Figure 6.8: Change in electrophysiological resting membrane potential of neurons resulting
from dynamic shear loadings of various strain amplitudes

Input resistance measurements reflect the neurons’ sensitivity to incoming electrical sig-

nals. Input resistance measurements were not significantly different between the shear strain

conditions tested (one-way ANOVA, p>0.05). Neurons exposed to 0% shear strain displayed
an average input resistance of 64.1 ± 39.1 MΩ (7 cells from 1 animal, shown in figure 6.9

(Black)). While neurons exposed to 10% shear strain displayed an average input resistance
of 111.4 ± 63.9 MΩ (8 cells from 1 animal, shown in figure 6.9 (Blue)) and neurons exposed
to 20% shear strain possessed an average input resistance of 69.5 ± 41.3 MΩ (10 cells from

2 animal, shown in figure 6.9 (Red)).

127

Figure 6.9: Change in electrophysiological input resistance of neurons resulting from dy-
namic shear loadings of various strain amplitudes

Time constant measurements provide insight into the neurons’ membrane capacitance

and membrane resistance. The time constant measurements were not significantly different

between the shear strain conditions tested (one-way ANOVA, p>0.05). Neurons exposed to
0% shear strain had an average time constant of 13.0 ± 4.1 ms (7 cells from 1 animal, shown

in figure 6.10 (Black)). While neurons exposed to 10% shear strain had an average time
constant of 18.7 ± 11.1 ms (8 cells from 1 animal, shown in figure 6.10 (Blue)) and neurons
exposed to 20% shear strain had an average time constant of 17.1 ± 5.2 ms (8 cells from 2

animal, shown in figure 6.10 (Red)).

Figure 6.10: Change in electrophysiological time constant of neurons resulting from dy-
namic shear loadings of various strain amplitudes

Overall, none of the intrinsic properties measured showed any significant difference be-

128

tween the specimens exposed to dynamic shear loads of various strain amplitudes and the

control specimens. These results suggest that these dynamic shear loadings do not produce

an acute (< 12 hours) effect on the intrinsic cell health of layer 5 pyramidal neurons of the

mouse barrel cortex.

Action Potential Output: Whole-cell current-clamp recordings were used to measure

the action potential output of the neurons tested. Frequency of action potential firing versus

current injected curves (F-I curves) were averaged for the three different shear strain con-

ditions: 0% shear strain (control) n=4, 10% shear strain n=5, and 20% shear strain n=8.

These curves are shown in figure 6.11.

Figure 6.11: Frequency of neuronal action potential firing versus current injected curves
(F-I curves) for three dynamic shear loading conditions

On average higher current injections in control (0% shear strain) neurons tended to result

in higher firing rates than the two shear strain conditions (10% shear strain and 20% shear
strain). 0% shear strain neurons had an average maximum frequency of 91 ± 31 Hz, while
10% shear strain was 64 ± 12 Hz and 20% shear strain was 60 ± 25 Hz. These sample

129

sizes are relatively small, and more samples will be needed for statistical analysis. However,

these results suggest there may be an acute reduction of action potential output of layer 5

pyramidal neurons of the mouse barrel cortex following dynamic shear loads.

6.2

In-vitro Experiments

This section discusses the preliminary in-vitro experiments involving HeLa (Henrietta Lacks)

cell cultures and cerebral organoids. The HeLa cell culture experiments were performed

in collaboration with the Center for Applied Genomics (CAG) at Children’s Hospital of

Philadelphia. The cell culture preparation and post-exposure analysis for the tests were

performed by Dr. Michael March at CAG. These tests are discussed in section 6.2.1. The

in-vitro experiments involving cerebral organoids were performed in collaboration with Dr.

Zane Lybrand’s group at Texas Woman’s University. The cerebral organoids preparation

and post-exposure analysis for the tests were performed by Dr. Zane Lybrand and Nikolas

Merlock. The details of these tests are provided in section 6.2.2.

6.2.1 In-vitro Experiments involving HeLa cell cultures

6.2.1.1 Experimental Procedure

Specimen Preparation:

In this study, our first-iteration prototype of the pressure appa-

ratus (discussed in chapter 3) was used to expose HeLa cell cultures to fast-varying pressure

profiles of various amplitudes and frequencies. HeLa cells (immortalized cervical carcinoma

cell line) were grown on 12 mm #2 thickness glass coverslips. For easy handling of the glass

coverslip in and out of the pressure chamber, a semi-circular cylindrical plug made of 2 %

weight/volume concentration agarose gel was used to support the glass coverslip. The glass

130

coverslip placed inside the pressure chamber (Design iteration I) with the agarose gel support

is shown in figure 6.12. The pressure chamber was filled with the cell-culture media buffer

solution.

Figure 6.12: Pressure chamber (Design iteration I) containing a glass coverslip plated with
HeLa cells

Loading Protocol:

In order to test the response of living cells to various frequencies and

pressure levels of a blast exposure, HeLa cells were exposed to high-frequency pressure waves

having magnitudes of 120 kPa and 250 kPa and frequencies of 500 Hz, 1 kHz, and 3 kHz for

a duration of 8 ms. The parametric details of all these loading cycles are given in table 6.1.

A total of 4 sets of loading tests were performed, where each set included the 8 loading

protocols. After loading the coverslips, they were incubated until further analysis.

6.2.1.2 Results

To quantify the cell death, lactate dehydrogenase (LDH) cytotoxicity assays (CyQUANT

LDH Cytotoxicity Assay, Thermo Fisher Scientific) were performed on the coverslips. The

131

Table 6.1: Loading protocols for in-vitro experiments involving HeLa cell cultures

Loading Pressure Amplitude
Protocol

(kPa)

LP1
LP2
LP3
LP4
LP5
LP6
LP7

Control

250
120
250
120
250
120
120

-

Frequency

Loading Time

(Hz)

500
500
1000
1000
3000
3000

500, 1000, and 3000

-

(ms)

8
8
8
8
8
8
10
-

LDH cytotoxicity assay measures the release of a cytoplasmic enzyme into culture super-

natant as a measure of cell death. These assays were performed on two different sets of

coverslips after 6 hrs and 24 hrs from the exposure. The percent specific release obtained

from the LDH cytotoxicity assays as a measure of cell death for all the loading protocols is

shown in figure 6.13. It can be observed that the specific release measured after 6 hrs from

the exposure is significantly higher than the one measured after 24 hrs. This result again

suggests a time-dependent response of the cells to the pressure loadings.

Figure 6.13: LDH assay results for in-vitro experiments involving HeLa cell cultures

132

-10-50510152025303540LP1LP2LP3LP4LP5LP6LP7Control% Specific Release% Specific Release -LDH Cytotoxicity Assay6 Hours24 HoursAnother method called “Dead-Green” staining (Image-IT DEAD Green Viability Stain,

Thermo Fisher Scientific) was also used to quantify cell death. This staining method is

detailed in [110]. The “Dead-Green” staining uses a cell-impermeable dye that cannot enter

live cells but stains cells with damaged membranes green. Thus, by calculating the percentage

of the green stained area, cell death can be quantified. The “Dead-Green” staining was

performed on another set of coverslips 3 hours after the exposure, and results are shown in

figure 6.14.

Figure 6.14: “Dead-Green” staining results for in-vitro experiments involving HeLa cell
cultures

It was observed from both LDH cytotoxicity assay and the “Dead-Green” staining that

LP4 induced maximum damage to cells compared to the other loading protocols. Even

though LP3 loaded the cells with higher pressure amplitude than LP4 at the same loading

frequency, LP4 induced greater damage to the cells than LP3. These results appeared to be

133

Ratio –Area under thresholded Dead Green /  Area under thresholded nucleiLoading Protocol5 10X microscope fields per protocolDead-Green Stainingcounter-intuitive. To gain confidence in these experimental findings, the same experiments

were repeated. However, the results did not support the initial findings. One of the possible

reasons for this discrepancy could be the robustness of the HeLa cells. HeLa cells, which

is an immortalized cervical carcinoma cell line, are lab-grown and very robust compared to

neuronal cell lines. Thus, there could be a very low signal to noise ratio when the mechanical

damage was assessed on the HeLa cells. However, further experimentation is necessary to

test this hypothesis.

6.2.2 In-vitro Experiments involving Cerebral Organoids

6.2.2.1 Introduction

Cerebral organoids are 3D structured cultures of human astrocytes and neurons which have a

similar organization and physiology to human neuronal networks found in the cerebral cortex

[111]. These organoids also exhibit the complex oscillatory electrical currents - analogous to

brain waves [112]. Additionally, these organoids can survive for long durations of time to

follow chronic changes (as well as acute) after a traumatic brain injury (TBI) event.

Cerebral organoids develop brain architecture found in the human cortex: Cere-

bral organoids are grown by dissociating into single cells and plating to form embryoid bodies

(EB). Following a dual SMAD inhibition protocol [111, 113] for neural induction, two types

of spheroids are grown from EBs, cortical and subpallial spheroids. As represented in neural

development, excitatory glutamatergic neurons are generated from the cortical plate, which

is later infiltrated by migrating inhibitory interneurons from the medial ganglionic eminence

of the subpallium. To replicate this process, both cortical and subpallial spheroids are gen-

erated using FGF2 and EGF to generate neural progenitors, followed by BDNF and NT3 to

134

promote neuronal differentiation. To pattern subpallial spheroids, WNT and SHH pathways

are modulated using IWP-2 and SAG. Once cortical and subpallial spheroids are patterned,

the two are fused to generate complete cortical grafts (shown in figure 6.15A). The human

cortex is organized into identifiable layers based on cells found in each layer. The same or-

ganization and cell populations are identified in human cerebral organoids (shown in figure

6.15B). Early-stage cortical spheroids (30 DIV) express distinguishing laminar markers for

deep and superficial cortical layers.

In cortical spheroids, from 70-180 DIV, we observed

TBR1 and CTIP2 expression in deep layer neurons, SATB2 and BRN2 in upper-layer neu-

rons, and HOPX and FAM107, markers from outer radial glial cells (shown in figure 6.15C).

Importantly, cerebral organoids also generate cells essential to maintain the health of the

local brain environment. This includes the ventral zone (VZ) to form internal ventricles,

glial cells like astrocytes, and oligodendrocytes (shown in figure 6.15D). Altogether, these

cerebral organoids can be grown to generate the critical cells and organization found in the

human cortex.

Cerebral organoids generate functional brain networks similar to the human cor-

tex: Using multielectrode arrays (MEA), fused cerebral organoids that display complex

neural networks (shown in figure 6.16a) are plated on top of a 12-electrode grid, and electri-

cal neuronal activity can be measured (shown in figure 6.16b). From this MEA data, specific

neuronal activity is recorded and measured (shown in figure 6.16c). These events represent

extracellular action potentials, or spikes, that originate from neurons of the organoid located

in proximity to an electrode (shown in figure 6.16c bottom). The electoral neuronal activity

can further be filtered in the 9-14 Hz bandwidths to identify the network oscillations (seen

in figure 6.16d) and measure the synchrony across multiple electrodes. These network oscil-

135

Figure 6.15: Cerebral organoids reproduce human cortical organization. (A) Timeline for
stem cell protocol to grow cerebral organoid. (B) Compared to the human cortex, cerebral
organoids contain most of the same layered organization and cells present. (C) Characteriza-
tion of cortical structure from deep and upper-layer neurons in developing cortical spheroid
(dotted line). (D) Image of the ventricular zone (VZ), astrocytes, and oligodendrocytes that
form a support niche.

136

1mmAstrocytesOligodendrocytesDAPIGFAPS100βDAPIBMPVentriclesDeep layerUpper layerCortical NicheCerebral OrganoidABCDlations are reminiscent of cortical brain waves measured using electroencephalogram (EEG)

in human patients (shown in figure 6.16d).

(a) Complex Neural Networks

(b) MEA Analysis

(c) Neural Activity

(d) Network Synchrony

Figure 6.16: Neurophysiology of the cerebral organoids. (a) Cerebral organoid displaying
a complex network of neurons. (b) Multielectrode array (MEA) measures 12 channels from
each organoid.
(c) Neural activity is measured from the MEA. The bottom image is a
magnified time scale of the voltage signal recorded in the top image. (d) The raw signals
can be filtered, and network oscillations are measured to determine the synchrony of activity
recorded in cerebral organoids.

In this study, the second-iteration prototype of our pressure apparatus (discussed in

chapter 3) was used to expose cerebral organoids to a fast-varying pressure profile, similar

to an intracranial pressure profile induced by a blast exposure. The pressure chamber of

the apparatus was filled with culture media along with the organoids, where a piezoelectric

137

actuator compressed the liquid rapidly, generating high-pressure fluctuations with controlled

pressure amplitude and frequency. Following exposure, cerebral organoids were monitored

for neurophysiological properties immediately after exposure (1 hour) and 24 hours later.

We measured an increase in neuronal activity 1 hour after exposure and an increase in the

synchronization of network function. These results demonstrate that in response to the

specific pressure parameters, cerebral organoids’ physiology changes.

6.2.2.2 Experimental Procedure

To determine if cerebral organoids respond to pressure forces similarly observed in a blast

injury, 60 days old fused cerebral organoids were loaded into the pressure apparatus. The

timeline for the tests is shown in figure 6.17a. Briefly, the blast chamber (shown in figure

6.17b) was filled with organoid growth media (orange), and organoids for each group were

placed inside. Later the chamber was closed off without introducing any bubbles within the

chamber. The pressure apparatus was then set and activated to deliver a 3 kHz frequency

waveform that reached 250 kPa pressure (shown in figure 6.17c).

Immediately following

pressure exposure, cerebral organoids were removed and plated on a pre-warmed MEA plate

(shown in figure 6.17d) where they were positioned on the electrode grid (shown in figures

6.17e and 6.17f). For control tests, cerebral organoids were placed inside the pressure cham-

ber, but the piezoelectric actuator was not activated. Once all organoids were plated, the

MEA plate was incubated until 1 hour after the blast exposure. At that time, an initial

2 minutes recording was acquired, and the MEA was returned to incubation over night at
37 ◦C with 5% CO2. The next day at 24 hours post-blast, the MEA plate was removed from

the incubator, and another 2 minutes MEA recording was acquired.

138

(a)

(b)

(c)

(d)

(e)

(f )

Figure 6.17: Blasting cerebral organoids. (a) Timeline of blast experiment (b) Image of
3 cerebral organoids inside of pressure chamber. (c) The pressure waveform recorded inside
the chamber during the exposure (Amplitude = 250 kPa, Frequency = 3 kHz, and duration
= 8 ms).
(e) Empty
MEA well showing the layout of 12 electrodes. (f) MEA well filled with media and cerebral
organoid on electrodes.

(d) 24 well plate being loaded with organoids after being blasted.

139

0246810Time (ms)-50050100150200250300Pressure (kPa)6.2.2.3 Results

Extra-cellular action potentials were quantified between no blast controls (n=3) and blast

organoids (n=3) as shown in figures 6.18a and 6.18b. Over 120 seconds, the total number

of spikes (shown in figure 6.16c bottom) were quantified. Raster plots in figures 6.18a and

6.18b show events over the 120 s recording displayed for each of the 12 channels. One hour

after the exposure, there was a modest increase in the amplitude of the spike (shown in figure

6.18c) and frequency (shown in figure 6.18d) compared to no blast controls. Interestingly,

by 24 hours, this emerging difference was no longer present.

To determine if the increase in spike frequency measured at 1 hr post-blast resulted in a

population-level change of neural network function, the phase synchrony was measured using

a 5 Hz filter (shown in figure 6.19). In control organoids, spontaneous 5 Hz oscillations are

unsynchronized (shown in figure 6.19a). This is seen by a low degree of phase angle difference

measured between two channels of the MEA (shown in figure 6.19b). Following blast expo-

sure, organoid networks were observed to become more synchronized (shown in figure 6.19c)

as measured by the average length of phase synchronization (shown in figures 6.19d and

6.19e). In the phase synchronization chart shown in figure 6.19e, the phase synchronization

length of 1 indicates a highly synchronous network, while the phase synchronization length

of 0 indicates network without any synchrony. The typical phase synchronization length for

a human EEG is between 0.2 and 0.5.

Together these results suggest that the blast parameters testing with the pressure appa-

ratus is capable of acutely altering neuronal function and network synchrony. Further testing

will expand parameter space and be used to investigate the cell and molecular biology of

blast and blunt forces on complex neural networks.

140

(a)

(b)

(c)

(d)

Figure 6.18: Blast acutely stimulates cerebral organoids. (a) Raster plot of spike events
recorded from MEA in no blast control organoids. (b) Raster plot of spike events recorded
from blast group at 1 hr post-blast. (c) Quantification of spike amplitudes in 1 hr and 24 hr
post-blast organoids. (d) Quantification of spike number in cerebral organoids.

141

(a) Control organoids

(b) Control organoids

(c) Blast organoids

(d) Blast organoids

(e) Phase Synchronization

Figure 6.19: Blast acutely increases network synchrony. (A) Two channels with raw signal
filtered using a 5 Hz Morlet Wavelet from control organoids. (B) The phase angle difference
plotted for each time unit between two channels (black), and the average angle and phase
synchronization plotted (green bar) from control organoids. (C) Filtered channels from blast
organoids.
(D) Plot of phase angle difference (black) and average phase synchronization
(green) plotted from blast organoids. (E) Quantification of phase synchronization between
the control and blast organoids at 1 hr post-blast. (n=3 organoids per group; p-values less
than 0.05 is considered statistically significant)

142

Chapter 7

Future Directions

This chapter summarizes the development of the experimental systems designed for bTBI-

related ex-vivo and in-vitro studies, and also provides directions for future advancements.

7.1 Future Directions for the MMHB Actuator

The MMHB actuator discussed in Chapter 2 replaces the incident bar of a conventional

split-Hopkinson pressure bar apparatus with a multi-material incident bar to transform the

incident stress pulse into a complex loading profile. The design parameters for an MMHB

actuator, such as the materials and lengths of bar components, are adjusted so as to closely

replicate the user-defined target loading profile. This is achieved by coupling the numerical

simulation of the MMHB actuator with an optimization algorithm. As proof of concept,

the MMHB actuator was designed for generating various sinusoidal pressure profiles with

multi-frequencies and decaying characteristics.

The designed MMHB actuator could reasonably replicate the dynamics of specified two-

frequency target profiles. Increasing spectral complexity of the target pressure profiles re-

sulted in a relatively higher deviation between the target and generated pressure profiles,

pointing to the necessity of increasing the number of material bars in the design space.

However, as the number of bar components used to construct the multi-material incident

bar increases, the total length of the actuator will increase as well. This increased length

143

of the actuator may limit its applicability because of the laboratory space or experimen-

tal constraints. A strategy to be considered in the future is to use a parallel combination

of several multi-material bars to replicate a multi-frequency pressure profile better than a

single multi-material bar. While a strategy like this could have more flexibility than single

multi-material bars, its implementation will bring new challenges.

The numerical model of the MMHB actuator treats water as a nearly incompressible

linear-elastic fluid. This model does not incorporate the compressibility of water and the

cavitation phenomenon. Cavitation can saturate the negative pressure part of the pressure

loadings generated by the MMHB actuator. Thus, a more specific equation of state for the

water constituent that accounts for pressure-density variations and cavitation phenomenon

[5] should be considered in the future.

7.2 Future Directions for the Pressure Apparatus

The generator of broadband pressure cycles discussed in Chapter 3 uses a water-filled piston-

cylinder assembly driven by a piezoelectric actuator to generate complex and fast-varying

pressure profiles. The apparatus employs a feedforward controller to track user-specified

target pressure profiles. The versatility of this apparatus in producing complex pressure

profiles was demonstrated by generating a single pressure pulse with various pulse-widths

and magnitudes, an approximate Friedlander waveform, and a multi-modal waveform. The

apparatus successfully generated pressures up to 450 kPa at frequencies up to 5 kHz. This

apparatus was used perform the preliminary ex-vivo and in-vitro experiments discussed in

Chapter 6 (sections 6.1.1, 6.2.1, and 6.2.2).

In the future, this pressure apparatus can be scaled to include multiple piezoelectric

144

actuators. Actuating multiple piezoelectric actuators synchronously to compress the water

confined in the pressure chamber should increase the maximum achievable pressure. More-

over, the actuation of the various actuators can be precisely timed to generate even more

intricate pressure profiles by leveraging the superposition of multiple pressure profiles. Incor-

porating multiple actuators can increase the capabilities of the apparatus, but it will further

complicate the system and may reduce its usability.

7.3 Future Directions for the Cavitation Rig

The bench-top cavitation rig discussed in Chapter 4 uses a water-filled piston-cylinder as-

sembly which is driven by a high-speed linear motor. The motor is used to retract the piston,

increasing the pressure chamber’s volume and, consequently, generating a tensile or nega-

tive pressure. This negative pressure gives rise to cavitation. The apparatus can generate

cavitation events of various durations and intensities. The potential of this apparatus is

demonstrated by producing a cavitation event in water.

For the cavitation experiment discussed in Chapter 4, a few micro-bubbles were already

present in the pressure chamber before the piston was actuated. Thus, the negative pressures

generated by the apparatus ended up expanding these existing micro-bubble rather than

creating new cavitation bubbles. In the future, if the focus is to generate new cavitation

bubbles in a liquid, deliberate efforts are needed to remove all the pre-existing bubbles from

the liquid, which could be achieved using the methods reported in [25, 26, 30].

The curved surface of the pressure chamber cavity creates an optical distortion while

imaging the cavitation bubbles due to refraction. This optical distortion creates challenges

in analyzing the bubble dynamics based on high-speed images. One way to circumvent this

145

challenge is to design a new pressure chamber with a flat optical window for high speed-

imaging, similar to the ones reported in [26, 50].

The shock-waves resulting from the collapse of the bubbles sharply increased the pressure

inside the chamber. This pressure exceeded well beyond the pressure transducer’s rated

measurement range, which led to the saturation of the measured pressure profile. Thus, a

pressure transducer with a significantly larger pressure range is required to appropriately

measure the high pressure associated with the shock-waves.

7.4 Future Directions for the Shear Apparatus

The generator of dynamic shear cycles discussed in Chapter 5 uses a voice-coil actuator and

two parallel plates assembly. This apparatus can generate oscillatory shear loadings of various

strain amplitudes and frequencies at the levels relevant to bTBI. The first-iteration shear

apparatus was used to conduct preliminary material characterization experiments for PAA

gelatin specimens. These experiments also tested the implementation of the apparatus for

loading samples of living tissue and three-dimensional cell cultures. The insights from these

preliminary experiments led to the design improvements implemented in the second iteration

of the shear apparatus. The second-iteration shear apparatus employs a PID controller to

control the voice-coil actuator and replicates the user-specified shear strain profiles. It was

also used to perform bTBI-related ex-vivo experiments on mouse brain tissue specimens

(discussed in chapter 6, section 6.1.2).

In the future, this shear apparatus can be scaled to include two voice-coil actuators by

replacing the stationary bottom plate with a movable plate actuated by the second voice-

coil actuator. Actuating both the actuators synchronously to move the plates in opposite

146

directions should increase the maximum achievable shear strain for a given loading frequency

and specimen thickness. Moreover, by precisely timing the motion of the actuators, multi-

frequency dynamic loading cycles can be generated through wave superposition. Again,

incorporating multiple actuators can increase the capabilities of the apparatus, but it will

further complicate the system and may reduce its usability.

As the shear loading frequency increases, the inertial effects inducing the unsteady shear

deformation of the specimen increases. The inertial effects are significant when the speci-

men’s thickness is comparable to the elastic shear wavelength in the specimen. In such cases,

the specimen’s deformation is no longer a simple-shear deformation as the shear strain field

is inhomogeneous. Such an inhomogeneous strain field cannot be reliably estimated just

from the top-plate displacement data. Furthermore, imaging-based techniques are required

to quantify the inhomogeneous deformation of the specimen. Thus, reducing the specimen’s

thickness as much as possible is a definitive strategy to eliminate the undesired inertial

effects.

147

BIBLIOGRAPHY

148

BIBLIOGRAPHY

[1] L. Wang, K. Labibes, Z. Azari, and G. Pluvinage, “Generalization of split Hopkinson
bar technique to use viscoelastic bars,” International Journal of Impact Engineering,
vol. 15, no. 5, pp. 669–686, 1994.

[2] ANSYS, “Ansys academic research mechanical, release 18.1.”

[3] A. Yucesoy, “The role of morphology and residual stress on blast-induced traumatic

brain injury,” Master’s thesis, Michigan State University, 2019.

[4] R. K. Gupta and A. Przekwas, “Mathematical models of blast-induced tbi: current

status, challenges, and prospects,” Frontiers in neurology, vol. 4, p. 59, 2013.

[5] X. Tan, A. Przekwas, and R. Gupta, “Computational modeling of blast wave inter-
action with a human body and assessment of traumatic brain injury,” Shock Waves,
vol. 27, no. 6, pp. 889–904, 2017.

[6] A. Przekwas, H. T. Garimella, X. G. Tan, Z. Chen, Y. Miao, V. Harrand, R. H. Kraft,
and R. K. Gupta, “Biomechanics of blast tbi with time-resolved consecutive primary,
secondary, and tertiary loads,” Military medicine, vol. 184, no. Supplement 1, pp.
195–205, 2019.

[7] E. Fievisohn, Z. Bailey, A. Guettler, and P. VandeVord, “Primary blast brain injury
mechanisms: current knowledge, limitations, and future directions,” Journal of biome-
chanical engineering, vol. 140, no. 2, p. 020806, 2018.

[8] S. Canchi, M. Sarntinoranont, Y. Hong, J. J. Flint, G. Subhash, and M. A. King,
“Simulated blast overpressure induces specific astrocyte injury in an ex vivo brain
slice model,” PloS one, vol. 12, no. 4, p. e0175396, 2017.

[9] H. Fischer, “Us military casualty statistics: operation new dawn, operation iraqi free-
dom, and operation enduring freedom.” Library of Congress Washington Dc Con-
gressional Research service, 2013.

[10] D. F. Meaney, B. Morrison, and C. D. Bass, “The mechanics of traumatic brain injury:
a review of what we know and what we need to know for reducing its societal burden,”
Journal of biomechanical engineering, vol. 136, no. 2, p. 021008, 2014.

[11] C. R. Bass, M. B. Panzer, K. A. Rafaels, G. Wood, J. Shridharani, and B. Capehart,
“Brain injuries from blast,” Annals of biomedical engineering, vol. 40, no. 1, pp. 185–
202, 2012.

149

[12] Y. C. Chen, D. H. Smith, and D. F. Meaney, “In-vitro approaches for studying blast-
induced traumatic brain injury,” Journal of neurotrauma, vol. 26, no. 6, pp. 861–876,
2009.

[13] A. Nakagawa, G. T. Manley, A. D. Gean, K. Ohtani, R. Armonda, A. Tsukamoto,
H. Yamamoto, K. Takayama, and T. Tominaga, “Mechanisms of primary blast-induced
traumatic brain injury: insights from shock-wave research,” Journal of neurotrauma,
vol. 28, no. 6, pp. 1101–1119, 2011.

[14] H. L. Lew, “Rehabilitation needs of an increasing population of patients: Traumatic
brain injury, polytrauma, and blast-related injuries,” Journal of Rehabilitation Re-
search and Development, vol. 42, no. 4, p. xiii, 2005.

[15] S. G. Scott, H. G. Belanger, R. D. Vanderploeg, J. Massengale, and J. Scholten,
“Mechanism-of-injury approach to evaluating patients with blast-related polytrauma,”
Journal of the American Osteopathic Association, vol. 106, no. 5, p. 265, 2006.

[16] J. M. Dewey, “The shape of the blast wave: studies of the friedlander equation,” in
Proceeding of the 21st International Symposium on Military Aspects of Blast and Shock
(MABS), Israel, 2010, pp. 1–9.

[17] N. Chandra, S. Ganpule, N. Kleinschmit, R. Feng, A. Holmberg, A. Sundaramurthy,
V. Selvan, and A. Alai, “Evolution of blast wave profiles in simulated air blasts: ex-
periment and computational modeling,” Shock Waves, vol. 22, no. 5, pp. 403–415,
2012.

[18] W. C. Moss, M. J. King, and E. G. Blackman, “Skull flexure from blast waves: a mech-
anism for brain injury with implications for helmet design,” Physical review letters, vol.
103, no. 10, p. 108702, 2009.

[19] R. Bolander, B. Mathie, C. Bir, D. Ritzel, and P. VandeVord, “Skull flexure as a
contributing factor in the mechanism of injury in the rat when exposed to a shock
wave,” Annals of biomedical engineering, vol. 39, no. 10, p. 2550, 2011.

[20] E. H. Clayton, G. M. Genin, and P. V. Bayly, “Transmission, attenuation and reflection
of shear waves in the human brain,” Journal of The Royal Society Interface, vol. 9,
no. 76, pp. 2899–2910, 2012.

[21] G. Nusholtz, P. Kaiker, and W. Gould, “Two factors critical in the pressure response
of the impacted head.” Aviation, space, and environmental medicine, vol. 58, no. 12,
pp. 1157–1164, 1987.

[22] J.-P. Franc and J.-M. Michel, Fundamentals of cavitation. Springer science & Business

media, 2006, vol. 76.

150

[23] C. E. Brennen, Cavitation and Bubble Dynamics. Cambridge University Press, 2013.

[24] J. R. Blake and D. Gibson, “Cavitation bubbles near boundaries,” Annual review of

fluid mechanics, vol. 19, no. 1, pp. 99–123, 1987.

[25] S. Canchi, K. Kelly, Y. Hong, M. A. King, G. Subhash, and M. Sarntinoranont, “Con-
trolled single bubble cavitation collapse results in jet-induced injury in brain tissue,”
Journal of the mechanical behavior of biomedical materials, vol. 74, pp. 261–273, 2017.

[26] Y. Hong, M. Sarntinoranont, G. Subhash, S. Canchi, and M. King, “Localized tis-
sue surrogate deformation due to controlled single bubble cavitation,” Experimental
Mechanics, vol. 56, no. 1, pp. 97–109, 2016.

[27] E.-A. Brujan, K. Nahen, P. Schmidt, and A. Vogel, “Dynamics of laser-induced cavi-
tation bubbles near an elastic boundary,” Journal of Fluid Mechanics, vol. 433, no. 1,
pp. 251–281, 2001.

[28] K. H. Taber, D. L. Warden, and R. A. Hurley, “Blast-related traumatic brain injury:
what is known?” The Journal of neuropsychiatry and clinical neurosciences, vol. 18,
no. 2, pp. 141–145, 2006.

[29] P. A. Taylor, J. S. Ludwigsen, and C. C. Ford, “Investigation of blast-induced traumatic

brain injury,” Brain injury, vol. 28, no. 7, pp. 879–895, 2014.

[30] M. Bustamante, D. Singh, and D. Cronin, “Polymeric hopkinson bar-confinement
chamber apparatus to evaluate fluid cavitation,” Experimental Mechanics, vol. 58,
no. 1, pp. 55–74, 2018.

[31] M. A. G. Sosa, R. De Gasperi, A. J. Paulino, P. E. Pricop, M. C. Shaughness,
E. Maudlin-Jeronimo, A. A. Hall, W. G. Janssen, F. J. Yuk, N. P. Dorr et al., “Blast
overpressure induces shear-related injuries in the brain of rats exposed to a mild trau-
matic brain injury,” Acta neuropathologica communications, vol. 1, no. 1, p. 51, 2013.

[32] M. Sarntinoranont, S. J. Lee, Y. Hong, M. A. King, G. Subhash, J. Kwon, and D. F.
Moore, “High-strain-rate brain injury model using submerged acute rat brain tissue
slices,” Journal of neurotrauma, vol. 29, no. 2, pp. 418–429, 2012.

[33] J. V. Rosenfeld, A. C. McFarlane, P. Bragge, R. A. Armonda, J. B. Grimes, and G. S.
Ling, “Blast-related traumatic brain injury,” The Lancet Neurology, vol. 12, no. 9, pp.
882–893, 2013.

[34] W. W. Chen and B. Song, Split Hopkinson (Kolsky) bar: design, testing and applica-

tions. Springer Science & Business Media, 2010.

151

[35] M. Nienaber, J. S. Lee, R. Feng, and J. Y. Lim, “Impulsive pressurization of neuronal
cells for traumatic brain injury study,” JoVE (Journal of Visualized Experiments),
no. 56, p. e2723, 2011.

[36] J. Goeller, A. Wardlaw, D. Treichler, J. O’Bruba, and G. Weiss, “Investigation of
cavitation as a possible damage mechanism in blast-induced traumatic brain injury,”
Journal of neurotrauma, vol. 29, no. 10, pp. 1970–1981, 2012.

[37] R. S. Salzar, D. Treichler, A. Wardlaw, G. Weiss, and J. Goeller, “Experimental in-
vestigation of cavitation as a possible damage mechanism in blast-induced traumatic
brain injury in post-mortem human subject heads,” Journal of neurotrauma, vol. 34,
no. 8, pp. 1589–1602, 2017.

[38] M. B. Panzer, K. A. Matthews, A. W. Yu, B. Morrison, D. F. Meaney, and C. R. Bass,
“A multiscale approach to blast neurotrauma modeling: part i–development of novel
test devices for in vivo and in vitro blast injury models,” Frontiers in neurology, vol. 3,
p. 46, 2012.

[39] G. B. Effgen, C. D. Hue, E. Vogel III, M. B. Panzer, D. F. Meaney, C. Bass, and B. Mor-
rison III, “A multiscale approach to blast neurotrauma modeling: part ii: methodology
for inducing blast injury to in vitro models,” Frontiers in neurology, vol. 3, p. 23, 2012.

[40] C. E. Needham, D. Ritzel, G. T. Rule, S. Wiri, and L. Young, “Blast testing issues
and tbi: experimental models that lead to wrong conclusions,” Frontiers in neurology,
vol. 6, p. 72, 2015.

[41] N. J. Logan, H. Arora, and C. A. Higgins, “Evaluating primary blast effects in vitro,”

JoVE (Journal of Visualized Experiments), no. 127, p. e55618, 2017.

[42] R. Campos-Pires, A. Yonis, W. Macdonald, K. Harris, C. J. Edge, P. F. Mahoney, and
R. Dickinson, “A novel in vitro model of blast traumatic brain injury,” JoVE (Journal
of Visualized Experiments), no. 142, p. e58400, 2018.

[43] R. A. Bauman, G. Ling, L. Tong, A. Januszkiewicz, D. Agoston, N. Delanerolle,
Y. Kim, D. Ritzel, R. Bell, J. Ecklund et al., “An introductory characterization of
a combat-casualty-care relevant swine model of closed head injury resulting from ex-
posure to explosive blast,” Journal of neurotrauma, vol. 26, no. 6, pp. 841–860, 2009.

[44] M. D. Alley, B. R. Schimizze, and S. F. Son, “Experimental modeling of explosive

blast-related traumatic brain injuries,” Neuroimage, vol. 54, pp. S45–S54, 2011.

[45] N. C. De Lanerolle, F. Bandak, D. Kang, A. Y. Li, F. Du, P. Swauger, S. Parks,
G. Ling, and J. H. Kim, “Characteristics of an explosive blast-induced brain injury in
an experimental model,” Journal of Neuropathology & Experimental Neurology, vol. 70,
no. 11, pp. 1046–1057, 2011.

152

[46] M. Risling, S. Plantman, M. Angeria, E. Rostami, B.-M. Bellander, M. Kirkegaard,
U. Arborelius, and J. Davidsson, “Mechanisms of blast induced brain injuries, experi-
mental studies in rats,” Neuroimage, vol. 54, pp. S89–S97, 2011.

[47] A. Nakagawa, M. Fujimura, K. Kato, H. Okuyama, T. Hashimoto, K. Takayama, and
T. Tominaga, “Shock wave-induced brain injury in rat: novel traumatic brain injury
animal model,” in Acta Neurochirurgica Supplements. Springer, 2008, pp. 421–424.

[48] N. E. Zander, T. Piehler, R. Banton, and R. Benjamin, “Bioeffects on an in vitro model
by small-scale explosives and shock wave overpressure impacts,” US Army Research
Laboratory Aberdeen Proving Ground United States, Tech. Rep., 2017.

[49] M. Bustamante and D. Cronin, “Cavitation threshold evaluation of porcine cere-
brospinal fluid using a polymeric split hopkinson pressure bar-confinement chamber
apparatus,” Journal of the mechanical behavior of biomedical materials, vol. 100, p.
103400, 2019.

[50] Y. Hong, S. Canchi, M. King, S. J. Lee, M. Sarntinoranont, and G. Subhash, “Devel-
opment of a test system to study brain tissue damage due to cavitation,” Am. Soc.
Biomech. Gainesville, Florida, vol. 2, 2012.

[51] W. Kang, Y. Chen, A. Bagchi, and T. J. OShaughnessy, “Characterization and de-
tection of acceleration-induced cavitation in soft materials using a drop-tower-based
integrated system,” Review of Scientific Instruments, vol. 88, no. 12, p. 125113, 2017.

[52] W. Kang and M. Raphael, “Acceleration-induced pressure gradients and cavitation in

soft biomaterials,” Scientific reports, vol. 8, no. 1, pp. 1–12, 2018.

[53] J. Wu and W. L. Nyborg, “Ultrasound, cavitation bubbles and their interaction with

cells,” Advanced drug delivery reviews, vol. 60, no. 10, pp. 1103–1116, 2008.

[54] W. Zhu, M. Alkhazal, M. Cho, and S. Xiao, “Microbubble generation by piezotrans-
ducer for biological studies,” Review of Scientific Instruments, vol. 86, no. 12, p.
124901, 2015.

[55] J. B. Estrada, C. Barajas, D. L. Henann, E. Johnsen, and C. Franck, “High strain-rate
soft material characterization via inertial cavitation,” Journal of the Mechanics and
Physics of Solids, vol. 112, pp. 291–317, 2018.

[56] B. Gerold, S. Kotopoulis, C. McDougall, D. McGloin, M. Postema, and P. Prentice,
“Laser-nucleated acoustic cavitation in focused ultrasound,” Review of Scientific In-
struments, vol. 82, no. 4, p. 044902, 2011.

153

[57] K. Upadhyay, A. Bhattacharyya, G. Subhash, and D. Spearot, “Quasi-static and high
strain rate simple shear characterization of soft polymers,” Experimental Mechanics,
pp. 1–15, 2019.

[58] M. Trexler, A. Lennon, A. Wickwire, T. Harrigan, Q. Luong, J. Graham, A. Maisano,
J. Roberts, and A. Merkle, “Verification and implementation of a modified split hop-
kinson pressure bar technique for characterizing biological tissue and soft biosimulant
materials under dynamic shear loading,” Journal of the mechanical behavior of biomed-
ical materials, vol. 4, no. 8, pp. 1920–1928, 2011.

[59] T. Mezger, The rheology handbook: for users of rotational and oscillatory rheometers.

European Coatings, 2020.

[60] K. B. Arbogast, K. L. Thibault, B. S. Pinheiro, K. I. Winey, and S. S. Margulies, “A
high-frequency shear device for testing soft biological tissues,” Journal of biomechanics,
vol. 30, no. 7, pp. 757–759, 1997.

[61] R. J. Okamoto, E. H. Clayton, and P. V. Bayly, “Viscoelastic properties of soft gels:
comparison of magnetic resonance elastography and dynamic shear testing in the shear
wave regime,” Physics in Medicine & Biology, vol. 56, no. 19, p. 6379, 2011.

[62] B. Rashid, M. Destrade, and M. D. Gilchrist, “Mechanical characterization of brain
tissue in simple shear at dynamic strain rates,” Journal of the mechanical behavior of
biomedical materials, vol. 28, pp. 71–85, 2013.

[63] M. C. LaPlaca, D. K. Cullen, J. J. McLoughlin, and R. S. Cargill II, “High rate shear
strain of three-dimensional neural cell cultures: a new in vitro traumatic brain injury
model,” Journal of biomechanics, vol. 38, no. 5, pp. 1093–1105, 2005.

[64] F. Perez Jr and M. Al-Haik, “Analysis of elastic stress wave propagation through a
complex composite structure,” International Journal of Theoretical and Applied Mul-
tiscale Mechanics, vol. 1, no. 2, pp. 164–175, 2009.

[65] T. Rahimzadeh, E. M. Arruda, and M. Thouless, “Design of armor for protection
against blast and impact,” Journal of the Mechanics and Physics of Solids, vol. 85,
pp. 98–111, 2015.

[66] C. Morganti-Kossmann, R. Raghupathi, and A. Maas, Traumatic brain and spinal cord

injury: challenges and developments. Cambridge University Press, 2012.

[67] M. Ghajari, P. J. Hellyer, and D. J. Sharp, “Computational modelling of traumatic
brain injury predicts the location of chronic traumatic encephalopathy pathology,”
Brain, vol. 140, no. 2, pp. 333–343, 2017.

154

[68] A. Jean, M. K. Nyein, J. Q. Zheng, D. F. Moore, J. D. Joannopoulos, and
R. Radovitzky, “An animal-to-human scaling law for blast-induced traumatic brain
injury risk assessment,” Proceedings of the National Academy of Sciences, vol. 111,
no. 43, pp. 15 310–15 315, 2014.

[69] S. T. Ahlers, E. Vasserman-Stokes, M. C. Shaughness, A. A. Hall, D. A. Shear,
M. Chavko, R. M. McCarron, and J. R. Stone, “Assessment of the effects of acute
and repeated exposure to blast overpressure in rodents: toward a greater understand-
ing of blast and the potential ramifications for injury in humans exposed to blast,”
Frontiers in neurology, vol. 3, p. 32, 2012.

[70] M. Chavko, T. Watanabe, S. Adeeb, J. Lankasky, S. T. Ahlers, and R. M. McCarron,
“Relationship between orientation to a blast and pressure wave propagation inside the
rat brain,” Journal of neuroscience methods, vol. 195, no. 1, pp. 61–66, 2011.

[71] I. Cernak, “Blast-induced neurotrauma models and their requirements,” Frontiers in

neurology, vol. 5, p. 128, 2014.

[72] B. Song and W. Chen, “Split hopkinson pressure bar techniques for characterizing soft
materials,” Latin American Journal of Solids and Structures, vol. 2, no. 2, pp. 113–152,
2005.

[73] C. Bacon, “An experimental method for considering dispersion and attenuation in a
viscoelastic hopkinson bar,” Experimental mechanics, vol. 38, no. 4, pp. 242–249, 1998.

[74] W. Chen, B. Zhang, and M. Forrestal, “A split hopkinson bar technique for low-

impedance materials,” Experimental mechanics, vol. 39, no. 2, pp. 81–85, 1999.

[75] O. S. Lee, S. S. You, J. H. Chong, and H. S. Kang, “Dynamic deformation under
a modified split hopkinson pressure bar experiment,” KSME International Journal,
vol. 12, no. 6, pp. 1143–1149, 1998.

[76] Q. Liu and G. Subhash, “Characterization of viscoelastic properties of polymer bar
using iterative deconvolution in the time domain,” Mechanics of Materials, vol. 38,
no. 12, pp. 1105–1117, 2006.

[77] A. S. Wineman and K. R. Rajagopal, Mechanical response of polymers: an introduc-

tion. Cambridge university press, 2000.

[78] H. Kolsky, Stress Waves in Solids, ser. Dover Books on Physics. Dover Publications,

1963.

[79] K. F. Graff, Wave Motion in Elastic Solids, ser. Dover Books on Physics Series. Dover

Publications, 1991.

155

[80] M. A. Meyers, Dynamic Behavior of Materials, ser. Wiley-Interscience publication.

Wiley, 1994.

[81] A. Bedford and D. Drumheller, Introduction to elastic wave propagation. John Wiley

& Sons Australia, Limited, 1994.

[82] P. Williams, P. Williams, and S. Brown, “Cavitation phenomena in water involving
the reflection of ultrasound pulses from a free surface, or from flexible membranes,”
Physics in Medicine & Biology, vol. 43, no. 10, p. 3101, 1998.

[83] Engineering Toolbox,

“Bulk Modulus

and Fluid Elasticity,”

https://www.

engineeringtoolbox.com/.

[84] K. Deb, Multi-objective optimization using evolutionary algorithms.

John Wiley &

Sons, 2001, vol. 16.

[85] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective
genetic algorithm: Nsga-ii,” IEEE transactions on evolutionary computation, vol. 6,
no. 2, pp. 182–197, 2002.

[86] K. Kido, Digital Fourier Analysis: Fundamentals. Springer, 2015.

[87] Y. Hong, M. Sarntinoranont, G. Subhash, S. Canchi, and M. King, “Localized tis-
sue surrogate deformation due to controlled single bubble cavitation,” Experimental
Mechanics, vol. 56, no. 1, pp. 97–109, 2016.

[88] M. H. Sadd, Elasticity: theory, applications, and numerics. Academic Press, 2009.

[89] H. K. Khalil, “Nonlinear systems,” Upper Saddle River, 2002.

[90] P. J. Antsaklis and A. N. Michel, Linear systems. Springer Science & Business Media,

2006.

[91] R. Kfoury, B. Marzban, E. Makki, M. L. Greenfield, and H. Yuan, “Effect of pressure
profile of shock waves on lipid membrane deformation,” PloS one, vol. 14, no. 2, p.
e0212566, 2019.

[92] B. Koo, B. Choi, H. Park, and K.-J. Yoon, “Past, Present, and Future of Brain

Organoid Technology,” Molecules and cells, vol. 42, no. 9, p. 617, 2019.

[93] D. St˚alhammar, “Experimental brain damage from fluid pressures due to impact accel-
eration. 1. design of experimental procedure,” Acta Neurologica Scandinavica, vol. 52,
no. 1, pp. 7–26, 1975.

156

[94] J. J. Crisco, B. J. Wilcox, J. T. Machan, T. W. McAllister, A.-C. Duhaime, S. M.
Duma, S. Rowson, J. G. Beckwith, J. J. Chu, and R. M. Greenwald, “Magnitude of
head impact exposures in individual collegiate football players,” Journal of applied
biomechanics, vol. 28, no. 2, pp. 174–183, 2012.

[95] A. Visioli, Practical PID control. Springer Science & Business Media, 2006.

[96] D. R. Lide, CRC handbook of chemistry and physics. CRC press, 2004, vol. 85.

[97] X. Feng, Z. Duan, Y. Fu, A. Sun, and D. Zhang, “The technology and application of
voice coil actuator,” in 2011 Second International Conference on Mechanic Automation
and Control Engineering.

IEEE, 2011, pp. 892–895.

[98] R. Cao, Z. Huang, T. Varghese, and G. Nabi, “Tissue mimicking materials for the de-
tection of prostate cancer using shear wave elastography: A validation study,” Medical
physics, vol. 40, no. 2, 2013.

[99] M.-K. Sun, J. Shieh, C.-W. Lo, C.-S. Chen, B.-T. Chen, C.-W. Huang, and W.-S. Chen,
“Reusable tissue-mimicking hydrogel phantoms for focused ultrasound ablation,” Ul-
trasonics sonochemistry, vol. 23, pp. 399–405, 2015.

[100] J. Kerwin, S. Vidhate, F. Masoomi, M. Tartis, A. M. Willis, and R. Mejia-Alvarez,
“Experimental study of the mechanics of blast-induced traumatic brain injury,” in
Mechanics of Biological Systems & Micro-and Nanomechanics, Volume 4. Springer,
2019, pp. 71–74.

[101] A. Wermer, J. Kerwin, K. Welsh, R. Mejia-Alvarez, M. Tartis, and A. Willis, “Materials
characterization of cranial simulants for blast-induced traumatic brain injury,” Military
medicine, vol. 185, no. Supplement 1, pp. 205–213, 2020.

[102] L. Adrian, R. J. Adrian, and J. Westerweel, Particle image velocimetry. Cambridge

University Press, 2011, no. 30.

[103] R. Gonzalez and R. Woods, Digital Image Processing. Pearson, 2018.

[104] W. Thielicke and E. J. Stamhuis, “PIVlab – towards user-friendly, affordable and
accurate digital particle image velocimetry in MATLAB,” Journal of Open Research
Software, vol. 2, oct 2014.

[105] J. D. Ferry, Viscoelastic properties of polymers. John Wiley & Sons, 1980.

[106] E. S. Lein, M. J. Hawrylycz, N. Ao, M. Ayres, A. Bensinger, A. Bernard, A. F.
Boe, M. S. Boguski, K. S. Brockway, E. J. Byrnes et al., “Genome-wide atlas of gene
expression in the adult mouse brain,” Nature, vol. 445, no. 7124, pp. 168–176, 2007.

157

[107] P. A. Taylor and C. C. Ford, “Simulation of blast-induced early-time intracranial wave
physics leading to traumatic brain injury,” Journal of biomechanical engineering, vol.
131, no. 6, 2009.

[108] S.-H. Lee, G. Govindaiah, and C. L. Cox, “Heterogeneity of firing properties among
rat thalamic reticular nucleus neurons,” The Journal of physiology, vol. 582, no. 1, pp.
195–208, 2007.

[109] J. A. Beatty, E. L. Sylwestrak, and C. L. Cox, “Two distinct populations of projec-
tion neurons in the rat lateral parafascicular thalamic nucleus and their cholinergic
responsiveness,” Neuroscience, vol. 162, no. 1, pp. 155–173, 2009.

[110] T. Delgado, P. A. Carroll, A. S. Punjabi, D. Margineantu, D. M. Hockenbery, and
M. Lagunoff, “Induction of the warburg effect by kaposi’s sarcoma herpesvirus is re-
quired for the maintenance of latently infected endothelial cells,” Proceedings of the
National Academy of Sciences, vol. 107, no. 23, pp. 10 696–10 701, 2010.

[111] F. Birey, J. Andersen, C. D. Makinson, S. Islam, W. Wei, N. Huber, H. C. Fan, K. R. C.
Metzler, G. Panagiotakos, N. Thom et al., “Assembly of functionally integrated human
forebrain spheroids,” Nature, vol. 545, no. 7652, pp. 54–59, 2017.

[112] C. A. Trujillo, R. Gao, P. D. Negraes, J. Gu, J. Buchanan, S. Preissl, A. Wang, W. Wu,
G. G. Haddad, I. A. Chaim et al., “Complex oscillatory waves emerging from cortical
organoids model early human brain network development,” Cell stem cell, vol. 25,
no. 4, pp. 558–569, 2019.

[113] S. P. Paca, “The rise of three-dimensional human brain cultures,” Nature, vol. 553, no.

7689, pp. 437–445, 2018.

158