evaluation of the mechanical behavior and material properties of native and tissue-engineered

155 Pages·2014·7.31 MB·English

EVALUATION OF THE MECHANICAL BEHAVIOR AND
MATERIAL PROPERTIES OF NATIVE AND
TISSUE-ENGINEERED CARTILAGE USING FINITE
ELEMENT ANALYSIS AND ULTRASONIC
ELASTOGRAPHY MEASUREMENT

by

CHEN-YUAN CHUNG

Submitted in partial fulfillment of the requirements
for the degree of Doctor of Philosophy

Dissertation Advisor: Professor Joseph M. Mansour

Department of Mechanical and Aerospace Engineering
CASE WESTERN RESERVE UNIVERSITY

January, 2015
CASE WESTERN RESERVE UNIVERSITY
SCHOOL OF GRADUATE STUDIES

We hereby approve the thesis/dissertation of
Chen-Yuan Chung
candidate for the Ph.D. degree*

(signed) Dr. Joseph M. Mansour
(Chair of the committee)
Dr. Jean F. Welter
(Member of the committee)
Dr. Ozan Akkus
(Member of the committee)
Dr. Umut Gurkan
(Member of the committee)

(date of defense) September 24, 2014

*We also certify that written approval has been obtained for any
proprietary material contained therein.
Table of Contents

Table of Contents………………………………………………………………............i
List of Tables………………………………………………………………………….iv
List of Figures…………………………………………………………………………v
Acknowledgements…………………………………………………………………...xi
Abstract……………………………………………………………………………....xii

Chapter 1: Introduction………………………………………………………………..1
1.1 Background………………………………………………………………...1
1.2 Aim of Research……………………………………………………………2

Chapter 2: Poroelastic Finite Element Analyses of Stress-relaxation in Cartilage……6
2.1 Abstract…………………………………………………………………….6
2.2 Introduction………………………………………………………………...6
2.3 Theory of Porous Media……………………………………………………8
2.4 Isotropic Model under Unconfined Compression………………………...10
2.4.1 Modelling and Boundary Setting…………………………………..11
2.4.2 Results and Discussion…………………………………………….13
2.5 Depth- and Strain-dependent Model under Unconfined Compression…...24
2.5.1 Equation-based Modelling………………………………………...24
2.5.2 Results and Discussion…………………………………………….26
2.6 Stress Relaxation of Cartilage under Simple Shear and Compression:
Experiments and Finite Element Analyses………………………………..30
2.6.1 Experimental Procedures…………………………………………..30
i
2.6.2 Finite Element Model……………………………………………...31
2.6.3 Results and Discussion…………………………………………….31

Chapter 3: Using Regression Models to Determine the Poroelastic Properties of
Cartilage……………………………………………………………………………...37
3.1 Abstract…………………………………………………………………...37
3.2 Introduction……………………………………………………………….37
3.3 Methods…………………………………………………………………...38
3.4 Results…………………………………………………………………….41
3.5 Discussion………………………………………………………………...43

Chapter 4: Determination of Poroelastic Properties of Cartilage Using Constrained
Optimization Coupled with Finite Element Analysis………………………………...55
4.1 Abstract…………………………………………………………………...55
4.2 Introduction……………………………………………………………….55
4.3 Methods…………………………………………………………………...58
4.4 Results…………………………………………………………………….61
4.5 Discussion………………………………………………………………...62

Chapter 5: Prediction of Hydrogel Material Properties from Ultrasonic Speed
Measurement for Functional Tissue Engineering Applications……………………...73
5.1 Abstract…………………………………………………………………...73
5.2 Introduction……………………………………………………………….73
5.3 Materials and Methods……………………………………………………74
5.4 Results and Discussion……………………………………………………77

ii
Chapter 6: Ultrasound Elastography for Estimation of Regional Strain of Multilayered
Hydrogels and Tissue-Engineered Cartilage…………………………………………84
6.1 Abstract…………………………………………………………………...84
6.2 Introduction……………………………………………………………….84
6.3 Methods…………………………………………………………………...86
6.3.1 Estimate the Depth-dependent Axial Strain of a Three-layered
Hydrogel Construct Using Ultrasound Elastography……………...86
6.3.2 Poroelastic Contact Finite Element Model for the Elastography
Measurement of a Three-layered Construct……………………….92
6.3.3 Poroelastic Contact Finite Element Model for a Verification
Example of Indentation……………………………………………93
6.4 Results…………………………………………………………………….96
6.5 Discussion………………………………………………………………...97
6.6 The Feasibility of Using Time-delay Estimation in Elastography to
Evaluate Mechanical Behavior of Developing Engineered Cartilage in a
Sterile Bioreactor………………………………………………………….99

Chapter 7: Conclusions and Future Studies………………………………………...124
7.1 Conclusions……………………………………………………………...124
7.2 Future Studies……………………………………………………………126

Appendix A: MATLAB Code for Time-delay Estimation from Pre- and
Post-compression Ultrasonic Signals……………………………………………….128

References…………………………………………………………………………..130
iii
List of Tables

Table 2.1 The major types of constitutive models proposed for articular cartilage..23
Table 3.1 Values of material properties used to simulate stress relaxation in the
transversely isotropic finite element model………………………………...53
Table 3.2 Predicted percentage error of material properties for combinations of five
regression equations at different time points during stress relaxation and the
corresponding value of  ………………………………………………...54
at
Table 4.1 Initial guesses of material properties (E ,E , , ,k) are set relatively
t a t at
far from one another in the space between their lower and upper bounds to
solve for the corresponding optimized material parameters……………….72
Table 4.2 Material properties (mean and standard deviation) were obtained from the
zero-order method (Case a) and gradient-based algorithm (Case b)……….72
Table 6.1 Material properties for the model of cartilage under a porous indenter
[100] where the value of Young’s modulus was also computed from the
aggregate modulus and Poisson’s ratio found in Spilker et al. [99]………120
Table 6.2 Material properties of three hydrogel layers were measured by the
indentation apparatus (Figure 6.6). These are the values of intrinsic material
properties that were used in the finite element model of the three-layered
hydrogel construct………………………………………………………...121
Table 6.3 Local strains predicted by the finite element model (FEA) and from
ultrasound elastography (US) for the last three steps of applied compression
on the top of the three-layered gel construct……………………………...121
Table 6.4 Global strains predicted by the definition in strength of materials () and
from ultrasound elastography (US) for the last three steps of applied
compression on the top of the three-layered gel construct………………..122
Table 6.5 Local strains predicted from ultrasound elastography using the equations
(6.8) to (6.10) where the estimated time delay was plus or minus the
sampling period, i.e. Ti 1108, i 1,2,3 at each step……………….122
d
Table 6.6 Internal displacement of TE cartilage for each region in Figure 6.18 under
three steps of compression………………………………………………..123
Table 6.7 Local strain of TE cartilage for each region in Figure 6.18 under three
steps of compression……………………………………………………...123

iv
List of Figures

Figure 1.1 The structure of articular cartilage showing four distinct regions:
superficial, middle, deep, and calcified cartilage zones……………………..2
Figure 2.1 Schematic of the unconfined compression test of a cylindrical disk of
hydrated cartilage…………………………………………………………..17
Figure 2.2 Specified axial displacement history for the numerical tests…………..17
Figure 2.3 Axisymmetric finite element model of the cartilage specimen………...18
Figure 2.4 Comparison of the normalized radial displacement of the cylindrical
boundary surface predicted using ANSYS with that predicted by the
analytical solution of Armstrong et al. [14] for a stress relaxation subsequent
to a Heaviside displacement function of axial compression……………….18
Figure 2.5 Plot of the normalized fluid pressure ([1]p E ) along the radius in
s 0
dimensionless form (r a) for different times after the application of a
Heaviside displacement……………………………………………….……19
Figure 2.6 Expanded results of fluid pressure ([1]p) to 90 about Y axis when
the time is 449.96 sec (t 0.07)…………………………..……………...19
n
Figure 2.7 Distributions of the normalized effective stress ({t  p } E ) along
rr 0 s 0
the dimensionless radius (r a) for different times………………………...20
Figure 2.8 Expanded results of effective radial stress (t  p ) to 90 about Y
rr 0
axis when the time is 89.992 sec (t 0.014)…………..…………………20
n
Figure 2.9 Comparison of the axial load intensity predicted using ANSYS with that
predicted by the analytical solution of Armstrong et al. [14] during a ramp
loading displacement and the subsequent relaxation process……………...21
Figure 2.10 Mesh convergence test of the finite element model by nodes,
corresponded to the case of  0 in Figure 2.9………….……………..21
s
Figure 2.11 The fluid velocity field (μm sec) is taken at nondimensional times
 t  (a) 0.04, (b) 0.1, (c) 0.13 and (d) 0.25. For the case shown  0.15
0 s
and results were plotted on the undeformed shape…………….…………...22
Figure 2.12 Simulation of experimental and computational results of Brown and
Singerman [21] using the proposed isotropic model……………………….23
Figure 2.13 Plot of the depth-dependent permeability profile of articular cartilage.
The permeability is seen to increase from the superficial zone and reach a
peak at one third zone before decreasing in the deep zone………………...27
v
Figure 2.14 Comparison of the axial load intensity predicted using COMSOL with
that predicted by the analytical solution of Armstrong et al. [14] during a
ramp loading displacement and the subsequent relaxation process………..28
Figure 2.15 Comparison of the axial load intensity predicted using COMSOL with
that predicted using ANSYS [45]. The permeability stays constant in the
linear isotropic model while it varies with depth and strain in the nonlinear
isotropic model..............................................................................................28
Figure 2.16 The fluid velocity field plotted on the deformed shape at  t 0.13
0
for the linear isotropic model (a) and the nonlinear isotropic model (b),
respectively: the contours and arrows represent fluid velocity quantified by
the left color legend (μm/s), the surface plot represents hydrostatic pressure
quantified by the right color legend (Pa). The results of the linear isotropic
model are in agreement with Figure 2.11(c) that was developed by
ANSYS…………………………………………………………..…………29
Figure 2.17 The designed device sat on the confocal microscope (Olympus
FluoViewFV1000), allowing the visualization of deformed and undeformed
photobleached lines. Displacements were applied in the directions of the
arrows shown…………………………………………………..…………...33
Figure 2.18 Deformed photobleached gridlines through the thickness of the sample
after the third ramp of shear displacement…………………………………33
Figure 2.19 Shear stress (a) and compressive stress (b) in the model with
E 0.6 MPa, 0.2 at the end of the third ramp (440 sec) in Figure
2.20...........................................…………………………………………….34
Figure 2.20 Stress relaxation for three steps of shear displacements………………35
Figure 2.21 History of compressive stress relaxation…………………….………..35
Figure 2.22 Compressive stress predicted by the model with E 0.6 MPa and
0.2 at the equilibrium (3000 sec) in Figure 2.21……………….…….36
Figure 3.1 Schematic of the unconfined compression test of a cylindrical disk of
hydrated cartilage……………………………………….………………….46
Figure 3.2 Comparison of the present isotropic and transversely isotropic models
using ANSYS with the analytical solution of those models and experimental
results from Cohen et al. [20]………………………………………………47
Figure 3.3 Parametric dependence of the non-dimensional load intensity on the
in-plane Poisson’s ratio  between 0 and 0.9; for the case shown, the
t
out-of-plane Poisson’s ratio  0………………………………………47
at
Figure 3.4 Parametric dependence of the ratio of the load intensity at peak to that at
equilibrium for various values of E E and  . All other material
t a t
properties are the same as used in Figure 3.2………………………………48
vi
Figure 3.5 Plots of load intensity at (a) 131 sec and (b) 1000sec as a function of
Young’s modulus in plane and out of plane where  0.49,  0 and
t at
k 5000 μm4 (μNsec)………………………….……..………………...48
Figure 3.6 Plots of load intensity at (a) 131sec and (b) 1000sec as a function of
Poisson’s ratio in plane and out of plane where E 4.3 MPa ,
t
E 0.64 MPa and k 5000 μm4 (μNsec)……………………………49
a
Figure 3.7 Plots of load intensity at (a) 131sec and (b) 1000sec as a function of
Poisson’s ratio in plane and Young’s modulus in plane where
E 0.64 MPa,  0 and k 5000 μm4 (μNsec)…………………49
a at
Figure 3.8 Plots of load intensity at (a) 131sec and (b) 1000sec as a function of
Young’s modulus out of plane and Poisson’s ratio out of plane where
E 4.3 MPa,  0.49 and k 5000 μm4 (μNsec)………………..50
t t
Figure 3.9 Plots of load intensity at (a) 131sec and (b) 1000sec as a function of
Young’s modulus in plane and permeability where E 0.64 MPa ,
a
 0.49 and   0…………………………………………….………50
t at
Figure 3.10 Plots of load intensity at (a) 131sec and (b) 1000sec as a function of
Poisson’s ratio in plane and permeability where E 4.3 MPa ,
t
E 0.64 MPa and  0………………………………………………51
a at
Figure 3.11 Plots of load intensity at (a) 131sec and (b) 1000sec as a function of
permeability and Young’s modulus out of plane where E 4.3 MPa,
t
 0.49 and  0…………………………………………………….51
t at
Figure 3.12 Plots of load intensity at (a) 131sec and (b) 1000sec as a function of
Young’s modulus in plane and Poisson’s ratio out of plane where
E 0.64 MPa,  0.49 and k 5000 μm4 (μNsec)……………..52
a t
Figure 3.13 Plots of load intensity at (a) 131sec and (b) 1000sec as a function of
permeability and Poisson’s ratio out of plane where E 4.3 MPa ,
t
E 0.64 MPa and  0.49……………………………………………52
a t
Figure 3.14 Plots of load intensity at (a) 131sec and (b) 1000sec as a function of
Young’s modulus out of plane and Poisson’s ratio in plane where
E 4.3 MPa ,  0 and k 5000 μm4 (μNsec)…………………53
t at
Figure 4.1 The flow of the design optimization using the zero-order method
(subproblem approximation method) established within the ANSYS
parametric design language………………………………………………...66
vii
Figure 4.2 Flow chart of finite element analysis integrated with the gradient-based
optimization algorithm for the determination of material properties from
experimentally measured data……………………………………..……….67
Figure 4.3 History of the first design variable E evolved from (a) zero-order
t
method and (b) gradient-based algorithm during optimization…………….68
Figure 4.4 History of the second design variable E evolved from (a) zero-order
a
method and (b) gradient-based algorithm during optimization……….……68
Figure 4.5 History of the third design variable  evolved from (a) zero-order
t
method and (b) gradient-based algorithm during optimization…………….69
Figure 4.6 History of the forth design variable  evolved from (a) zero-order
at
method and (b) gradient-based algorithm during optimization…………….69
Figure 4.7 History of the fifth design variable k evolved from (a) zero-order
method and (b) gradient-based algorithm during optimization…………….70
Figure 4.8 History of the objective function evolved from (a) zero-order method and
(b) gradient-based algorithm during optimization………………………….70
Figure 4.9 Experimental unconfined compression stress-relaxation data, analytical
solution with the least-squares optimization [20] and corresponding finite
element simulations with the zero-order and gradient-based optimizations
for the transversely isotropic model………………………………………..71
Figure 5.1 Schematics of apparatus included an ultrasound transducer in series with
a load cell. The face of the transducer was acoustically coupled to an open
bioreactor containing the gel sample……………………………………….79
Figure 5.2 Young’s and storage moduli were mechanically derived using the
Rheometrics RSA-II Solids Analyzer (Rheometrics Inc., Piscataway, NJ)..79
Figure 5.3 Box plot for SOS through the hydrogels grouped according to various
concentrations………………………………………………………………80
Figure 5.4 Box plot of population data showing the variations in (a) Young’s
modulus and (b) storage modulus between the concentrations…………….81
Figure 5.5 Linear correlation between SOS and agarose concentration (volume
fraction)………………………………………………………………..…...82
Figure 5.6 Linear correlation between (a) Young’s modulus as well as (b) storage
modulus derived from conventional mechanical tests and the aggregate
modulus derived from ultrasound………………………………………….83
Figure 6.1 Schematic illustration of the custom designed elastography test rig….104
Figure 6.2 Illustration of the pre- and post-compression echo signals and the two
observation windows used for computing the longitudinal strain in tissue.
The pre-compression and post-compression signals in the indexed window
pair are compared by cross correlation……………………………………105
viii