[PDF] Continuum Mechanical Investigations of the Intervertebral Disc

Download Continuum Mechanical Investigations of the Intervertebral Disc...

Diss. ETH Nr. 19545

Continuum Mechanical Investigations of the Intervertebral Disc A dissertation submitted to ETH ZURICH for the degree of Doctor of Sciences

presented by ¨ JORG HELFENSTEIN Dipl. Masch.-Ing. ETH born September 15, 1981 citizen of Sempach, Lucerne

accepted on the recommendation of Prof. Dr. E. Mazza, examiner Prof. Dr. M. B. Rubin, co-examiner Prof. Dr. J. G. Snedeker, co-examiner

2011

Acknowledgments

The present work was part of the CO-ME (COmputer aided and image guided MEdical interventions) Project 9, financed by the Swiss National Science Foundation. This thesis was developed during my employment as a research assistant at the Institute of Mechanical Systems, ETH Zurich, under the guidance of Prof. Dr. Edoardo Mazza, to whom I owe gratitude for giving me the opportunity to work in such a challenging research field. I would like to thank Prof. Dr. Miles B. Rubin and Prof. Dr. Jess G. Snedeker, my co-examiners, for providing insightful comments and carefully reviewing this thesis. I am grateful to all the people who gave their contribution to this work, in particular to: Prof. Dr. Sanjay Govindjee and Dr. Mahmood Jabareen, who helped to advance the subject about a non-physical response in constitutive models, Giuseppe Barbarino, Marc Hollenstein, Ondrej Papes and Andreas Schifferle, who contributed to this thesis with their expertise in discussions and collaborations, Prof. Dr. Brigitte von Rechenberg, who gave permission to collect ovine spines from the Veterinary Clinic of Zurich, her students Dr. Karina Klein, Dr. Mich`ele Sidler and Dr. Anja-Christina Waselau, and the butcher Mr. Bruno Gerzner, who helped with slaughtering and sample taking, Dr. Thomas M¨ uller, Dr. Stephan Weiss and Michael Wyss, who provided CT- and MRscans of the functional spinal units, Michael M¨oller and Adrian Schweizer, who helped me with their programming skills, Christian Maier, whose knowledge about design and construction is far better than mine, Nicolas Jochum, to whom I still owe a beer for milling aluminum casting boxes for the cruciform sample preparation, Dr. Philippe B¨ uchler, who provided access to the software package AMIRA, Max Gay, Jan Hunger, Evelin Mattle, Joris Moerkerken, Roman Schneider and

Michael Windler, who contributed with their bachelor or master theses to this research, Beni Cadonau, Traude Junker, Dr. Stephan Kaufmann, Urs Loy, Ueli Marti, Gabi Squindo, Daniela Tanner and Jean-Claude Tomasina for their technical support and for providing the infrastructure at the Center of Mechanics, all professors and colleagues at the Center of Mechanics for the excellent working atmosphere and the good times spent together. Finally, I am most grateful to all members of my family and my wife Ursula, for their constant support, patience and encouragement throughout my dissertation: this work is dedicated to them.

J¨org Helfenstein Zurich, March 2011

Contents

Abstract

iii

Zusammenfassung

Nomenclature

vii

1 Introduction

2 Spine Anatomy 2.1 The spinal column . . . . . . . . . . . . . . . 2.2 The intervertebral disc . . . . . . . . . . . . . 2.3 Loading of the intervertebral disc . . . . . . . 2.4 Comparison of human and ovine vertebrae and

. . . .

7 8 10 14 16

. . . .

19 19 21 23 33

Response in Constitutive Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35 36 41 44

3 Nonlinear Continuum Mechanics 3.1 Kinematics . . . . . . . . . . . 3.2 Stress measures . . . . . . . . . 3.3 Constitutive relations . . . . . . 3.4 Elasticity tensors . . . . . . . . 4 A Non-Physical 4.1 Methods . . 4.2 Results . . . 4.3 Discussion .

5 Specimen Design 5.1 Motivation . . 5.2 Methods . . . 5.3 Results . . . . 5.4 Discussion . . 5.5 Conclusion . .

. . . .

for Planar-Biaxial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . . . . . . . . . intervertebral

. . . .

Materials Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . . . . . . . discs

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

47 48 49 52 53 56 i

6 Uniaxial Experiments with Anulus Fibrosus 6.1 Deformation modes in the intervertebral disc 6.2 Pre-studies . . . . . . . . . . . . . . . . . . . 6.3 Sample preparation . . . . . . . . . . . . . . 6.4 Experiments . . . . . . . . . . . . . . . . . . 6.5 Results . . . . . . . . . . . . . . . . . . . . . 6.6 Discussion . . . . . . . . . . . . . . . . . . . 7 Experiments with Functional Spinal 7.1 Preliminaries . . . . . . . . . . . . 7.2 Protocols and results . . . . . . . . 7.3 Discussion . . . . . . . . . . . . . .

Samples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . .

57 58 61 67 69 71 72

Units 75 . . . . . . . . . . . . . . . . . . . . . 76 . . . . . . . . . . . . . . . . . . . . . 88 . . . . . . . . . . . . . . . . . . . . . 100

8 Calibration of Elasto-Viscoplastic Equations 8.1 The finite element model . . . . . . . . . . . . 8.2 Calibrations . . . . . . . . . . . . . . . . . . . 8.3 Results . . . . . . . . . . . . . . . . . . . . . . 8.4 Discussion . . . . . . . . . . . . . . . . . . . . 9 Simulations of Axial Loading and Bending 9.1 Reconstruction of the geometry . . . . . . . 9.2 The finite element model . . . . . . . . . . . 9.3 Results . . . . . . . . . . . . . . . . . . . . . 9.4 Discussion . . . . . . . . . . . . . . . . . . .

. . . .

. . . . . . . .

. . . .

103 104 105 108 110

. . . .

115 117 120 126 131

10 Conclusions and Outlook

135

A Polyconvexity

139

B Constitutive Equations in MATHEMATICA

143

C Routines and ABAQUS Input Files C.1 Customization of the optimization environment C.2 ABAQUS job generation: PYTHON script . . . C.3 Deformation modes . . . . . . . . . . . . . . . . C.4 MATLAB code for mesh generation . . . . . . . C.5 Fluid-cavity elements . . . . . . . . . . . . . . . C.6 User defined anisotropic materials in ABAQUS . C.7 Writing fields to .odb-files . . . . . . . . . . . . C.8 Reading .odb-files with PYTHON . . . . . . . .

. . . . . . . .

147 147 150 156 162 175 178 181 191

D Data Sheet of Commercial Products

193

References

195

Curriculum Vitae

211

Abstract

This thesis aims at providing novel insights in the whole procedure of model generation for biological systems. To this end, aspects as the selection of suitable material models, the performance of appropriate experiments for parameter determination of constitutive equations and the model validation are analyzed. It is not the objective to provide a model for biomedical applications, but rather to explore continuum mechanical related aspects of the modeling procedure. The vertebrate spinal column is chosen as example structure for its relevance in medicine and the amount of available literature data. An ovine functional spinal unit is modeled with the finite element method: while the vertebrae are represented with rigid bodies, the intervertebral disc is deformable. The latter is subdivided into the nucleus pulposus, which is described by an inviscid, incompressible fluid, and the anulus fibrosus (AF), which is modeled with anisotropic elasto-viscoplastic equations [Rubin and Bodner, 2002]. For first simulations of the AF tissue anisotropic hyperelastic equations [Holzapfel et al., 2000] were used as for these equations parameters are available. The results of computations with this equations did not match the expectations. In fact, simulations involving uniaxial stress configurations in transversely isotropic materials reveal volume growth at rather small stretches. The same phenomenon is also found with other constitutive equations, too. It is found that the additive split of the Helmholz free-energy into a volumetric and an isochoric part, applied to the matrix and the fiber contribution, is the origin of this shortcoming. A solution is presented how to avoid this problem on the constitutive level rather than to use numerical methods such as the Augmented Lagrangian method. The calibration of the constitutive equations used to describe the AF tissue was done with data from (quasi-) uniaxial tensile experiments. Also biaxial experiments were considered but finally excluded for the reasons described hereafter. Our own studies revealed that results from planar biaxial tests on regular cruciform specimen do not predict the true equibiaxial material behavior and are specimen geometry dependent. Only a small central area deforms equibiaxially, the measured global displacement and force values do not well represent the biaxial materials response. It is shown that already a low number of slots in the limbs of the cruciform samples helps to increase the area of hom*ogeneous equibiaxial loading significantly. Only, the outcomes highlight, with regard to the small size of the possible AF specimens and their anisotropy, the difficulties to iii

generate adequate biaxial stress fields for a faithful material parameters determination. Consequently uniaxial tensile experiments only were used to characterize the mechanical behavior of the AF tissue. The results of simulations of a functional spinal unit under cyclic axial loading and bending are compared with our own experimental data. The realization of appropriate experiments led to several challenges, including the design of a reproducible solution to the clamping problem by means of suitable mountings and the verification of the testing machines accuracy. Preliminary experiments with steel springs had the objective to demonstrate the performance of the constructed mountings in different load cases and to verify the suitability of the used materials testing machine. These experiments yielded improvements of the setup that would not have been realized without them. Not enough experiments were made in order to make them statistically representative, but the applied procedures for validation and measurement are thought to provide interesting insights in how future experiments should be designed. The comparison of the simulations with the experiments indicates that the finite element model behaves too stiff, especially at larger deformations. It is shown that the fibers in the AF contribute only to a minor extend to the high stiffness. Instead, two other sources are discussed which could lead to an overestimation of the AF stiffness: the low initial stiffness of the tissue which makes it difficult to define the unstressed configuration and the exponential form of the used constitutive equations.

Zusammenfassung

Diese Arbeit liefert neue Einsichten in den gesamten Prozess der Modellgenerierung f¨ ur biologische Systeme. Zu diesem Zweck werden Aspekte wie die Auswahl von geeigneten Materialmodellen, die Durchf¨ uhrung geeigneter Experimente f¨ ur die Parameterbestimmung von konstitutiven Gleichungen und die Modellvalidierung analysiert. Es ist nicht das Ziel ein Modell f¨ ur biomedizinische Anwendungen zu entwickeln, vielmehr werden kontinuumsmechanische Aspekte der Modellierung untersucht. Wegen ihrer Relevanz f¨ ur die Medizin und der Menge an verf¨ ugbarer Literatur wurde die Wirbels¨aule als Beispiel gew¨ahlt. Eine funktionelle Einheit einer Schafwirbels¨aule wurde mit der Methode der Finiten Elemente modelliert. Die Wirbelk¨orper werden als starr angenommen, die Bandscheibe hingegen als ein deformierbarer K¨orper. In der Bandscheibe wird zudem zwischen dem Nucleus pulposus, welcher durch eine reibungsfreie, inkompressible Fl¨ ussigkeit beschrieben wird, und dem Anulus fibrosus (AF), modelliert mit anisotropen elasto-viskoplastischen Gleichungen [Rubin and Bodner, 2002], unterschieden. F¨ ur erste Simulationen des AF-Gewebes wurde ein anisotropes hyperelastisches Stoffgesetz [Holzapfel et al., 2000] verwendet, da f¨ ur dieses Parameter vorhanden sind. Die Ergebnisse von Berechnungen mit diesen Gleichungen entsprachen nicht den Erwartungen. In der Tat zeigen Simulationen von einachsigen Spannungszust¨anden in transversalisotropen Materialien schon bei kleinen Dehnungen Volumenwachstum. Das gleiche Ph¨anomen wurde auch mit anderen Stoffgesetzen gefunden. Als Ursache wurde die additive Trennung der freien Helmholz-Energie in einen volumetrischen und einen deviatorischen Teil, angewandt auf die Matrix und den Faserbeitrag, identifiziert. Es wird ein Vorgehen pr¨asentiert, wie dieses Problem auf Ebene der Stoffgesetze gel¨ost werden kann, anstatt numerische Methoden wie die Augmented-Lagrangian-Methode zu verwenden. F¨ ur die Kalibrierung der Stoffgesetze zur Beschreibung des AF-Gewebes wurden (quasi-) einachsige Experimente durchgef¨ uhrt. Biaxiale Versuche wurden ebenfalls in Betracht gezogen, aber schliesslich aus den nachfolgend beschriebenen Gr¨ unden ausgeschlossen. Eigene Studien ergaben, dass Ergebnisse aus planaren biaxialen Versuchen mit regelm¨assigen kreuzf¨ormigen Proben nicht das wahre equibiaxial Werkstoffverhalten beschreiben und zudem Geometrie abh¨angig sind. Nur ein kleiner zentraler Bereich verformt sich equibiaxial, die gemessenen globalen Verschiebungs- und Kraftwerte repr¨asentieren nur ungen¨ ugend das biaxiale Materialverhalten. Bereits eine geringe Anzahl von Langl¨ochern in den Armen der kreuzf¨ormigen Proben vergr¨ossert den Bereich v

unter hom*ogener equibiaxialer Belastung deutlich. Die Resultate zeigen ebenfalls, in Hinblick auf die geringe Gr¨osse der m¨oglichen AF-Proben und deren Anisotropie, die Schwierigkeiten auf ad¨aquate zweiachsige Spannungsfelder f¨ ur die zuverl¨assige Kalibration von Stoffgesetzen zu generieren. Deshalb wurden nur einachsige Zugversuche durchgef¨ uhrt um das mechanische Verhalten des AF-Gewebes zu charakterisieren. Die Ergebnisse der Simulationen einer funktionelle Einheit der Wirbels¨aule unter zyklischer axialer Belastung und Biegung werden mit eigenen experimentellen Daten verglichen. Die Realisierung der entsprechenden Experimente f¨ uhrte zu mehreren Herausforderungen, einschliesslich der Gestaltung geeigneter Halterungen und der Verifikation der Genauigkeit der verwendeten Materialpr¨ ufmaschine. Vorg¨angige Experimente mit Stahlfedern hatten das Ziel, die Einsetzbarkeit der konstruierten Halterungen in verschiedenen Lastf¨allen zu demonstrieren und die Eignung der verwendeten Pr¨ ufmaschine zu testen. Diese Experimente haben Verbesserungen des Setups hervorgebracht welche ohne sie nicht realisiert worden w¨aren. Es wurden nicht gen¨ ugend Versuche gemacht damit die Resultate statistisch repr¨asentativ w¨aren, allerdings geben die beschriebenen Prozeduren zur Validierung und Messung interessante Einblicke wie zuk¨ unftige Experimente gestaltet werden k¨onnen. Der Vergleich der Simulationen mit den Experimenten zeigt, dass sich das Finite Elemente Modell zu steif verh¨alt, vor allem bei gr¨osseren Verformungen. Es wird gezeigt, dass die Fasern im AF nur geringf¨ ugig zu der hohen Steifigkeit beitragen. Stattdessen ¨ werden zwei andere Quellen diskutiert welche zu einer Ubersch¨ atzung des AF-Steifigkeit f¨ uhren k¨onnen: die anf¨angliche Weichheit des Gewebes, die es schwierig macht die undeformierte (Ausgangs-) Konfiguration zu bestimmen und die exponentielle Form des verwendeten Stoffgesetzes.

Nomenclature

The notation in this thesis follows mainly the one used in Holzapfel [2000] which gives in its first chapter a good introduction into the algebra of vectors and tensors. In the following only a short recapitulation of the most important notions. Scalars are typeset in standard letters (e.g. α, a), vectors are printed bold (e.g. ζ, x), second order tensors are set bold and sanserif (e.g. σ, F). Higher order tensors are set in blackboard bold letters (e.g. C, P). In chapters directly related with continuum mechanics it is differentiated between quantities in the reference configuration (upper case letters) and quantities in the actual configuration (lower case letters). Indices designate components of tensor quantities (e.g. σij = ei ·σej is the {i, j}-component of the Cauchy stress σ). The Kronecker delta δij designates the components of the second order identity tensor 1 (δij = 1 [-] if i = j and 0 elsewhere). Square brackets define function arguments (e.g. f [x] is a function f with argument x). On the next page all symbols of mathematical operators used are defined and the most important abbreviations are listed. The nomenclature concludes on Page ix with a presentation of the anatomical coordinate system.

vii

Operators h•i . . . . . . . . . . . . Macaulay brackets, h•i = ¯ . . . . . . . . . . . . isochoric part (•) ˙ . . . . . . . . . . . . total time derivative (•)

1 2

(• + abs[•])

(•)0 . . . . . . . . . . . deviator (•)T . . . . . . . . . . transpose (•)vol . . . . . . . . . volumetric part (•) : (N) . . . . . . double contraction (•) · (N). . . . . . . contraction

(•) ⊗ (N) . . . . . tensor product

(•) × (N) . . . . . vector cross product ∇[•] . . . . . . . . . . gradient

abs[•] . . . . . . . . . absolute value adj[•] . . . . . . . . . adjugate d[•] . . . . . . . . . . . total derivative det[•] . . . . . . . . . determinant eig[•] . . . . . . . . . eigenvalues grad[•] . . . . . . . . gradient ln[•] . . . . . . . . . . natural logarithm max[•] . . . . . . . . maximum tr[•] . . . . . . . . . . trace

Abbreviations AF . . . . . . . . . . . Anulus Fibrosus CT . . . . . . . . . . . Computer Tomography FCE [∼s] . . . . . Fluid Cavity Element FE . . . . . . . . . . . Finite Element FSU [∼s] . . . . . Functional Spinal Unit IVD [∼s]. . . . . . InterVertebral Disc NP . . . . . . . . . . . Nucleus Pulposus PMMA . . . . . . . PolyMethylMethAcrylate; a thermoplastic SC . . . . . . . . . . . . Spinal Column STATEV [∼s] . solution-dependent STATE Variable; variables used in a UMAT UMAT. . . . . . . . User MATerial; ABAQUS user subroutine to define constitutive equations V [∼e] . . . . . . . . Vertebra

viii

The anatomical coordinate system

z compression

tension

left axial torsion right axial torsion

right lateral bending posterior shear x

extension left lateral bending anterior shear

flexion left lateral shear

right lateral shear y

Figure I.1: Illustration of the anatomical coordinate system and the definitions of different motions [after Panjabi et al., 1992].

frontal plane: y-z-plane median plane: x-z-plane sagittal planes: all planes coplanar to the median plane transversal plane: x-y-plane

Few tests have been made and the data are insufficient. However, the findings to date are of interest and encourage us to continue these investigations. Brown et al. [1957]

Chapter

Introduction In 2001 the Swiss National Science Foundation (SNSF) started a twelve years project called CO-ME (COmputer aided and image guided MEdical interventions)1 . The aim of this project is the development and evaluation of basics and applied technologies in the field of diagnosis, surgical planning and therapeutic intervention. The present work is part of phase II (years 2005 to 2009) and is included in a project cluster with the title “soft tissue modeling: from mechano-biology to real-time simulation”. The target of this cluster is an improved understanding of soft-tissue mechanics including the structural multiscale relationships, the modeling of tissue behavior and the enhancement of simulation and visualization tools. This thesis aims at providing novel insights in the whole procedure, from the selection of suitable model formulations to the performance of appropriate experiments, for parameter determination of nonlinear, anisotropic, viscoelastic constitutive equations and finally to the model creation of a biological system. It is not its objective to provide a model for biomedical applications, but rather to explore continuum mechanical related aspects of the modeling procedure. For its relevance in medicine and the amount of available literature data, the vertebrate spine was chosen as example structure. In industrialized Western countries back pain is a major condition, from which almost everybody suffers once or more in his lifetime.i The number of surgical interventions and other medical treatments is constantly increasing. Back pain is not only a serious health problem but has an important economic impact by its enormous costs for medical treatments and work compensations [Frank et al., 1996] as well. A better understanding of the complex, coupled behavior of the intervertebral disc (IVD) mechanics is needed to identify the causes of back pain and to give best treatment to patients. Realistic simulations of the spine might serve as numerical laboratories [e.g. Magnier et al., 2009] and are an aid for the design of novel implants [e.g. Ferguson et al., 2006]. The key hereto is the qualitative and quantitative understanding of the mechanical properties and their role and importance in simulations. 1 i

http://co-me.ch Andersson [1999]; Bao et al. [1996]; Parent-Thirion et al. [2007]; Raspe et al. [2008]

CHAPTER 1. INTRODUCTION

In the last decades many constitutive models for soft biological tissues were proposed, some of them are rather complex and use a large number of parameters.i The precise determination of material parameters with respect to sound experimental data is necessary to assess the predictive power of those models. Experiments should be designed in a way to address individual features of the models, e.g. elastic properties, viscous properties, etc. In such experiments it is essential to minimize any undesired influence from the testing environment on the measurements [e.g. Helfenstein et al., 2009b]. In the remaining part of this chapter, the main subjects of investigation on theory, experiment and simulation related aspects of the present work are introduced. At the same time the structure of the document is presented.

Outline of the thesis The first two chapters familiarize the reader with the basics of the spine anatomy and the mechanics of continuous deformable bodies. Chapter 2 gives an overview on the anatomy of the (human) spinal column and, in particular, of the IVD. Also, the physiological loading of the IVD is discussed. As human tissues for in-vitro studies are difficult to obtain the question arises if experiments cannot be with using animal material. Except for the simpler supply, other positive aspects of animal subjects are the smaller variability and the availability of young subjects. In this study samples from sheep were used and consequently the chapter concludes with a comparison of the human and the ovine spine anatomy. Chapter 3 presents the basic aspects of nonlinear continuum mechanics. It starts with the introduction of kinematic quantities and continues with the stress measures. Thereafter constitutive equations used throughout the theses are explained and eventually the elasticity tensors are defined.

Theoretical aspects Several anisotropic elastic constitutive equations exist that try to capture the essential behavior of anisotropic tissues, e.g. the anulus fibrosus (AF).ii As for the equations of Holzapfel et al. [2000] parameters for the AF-tissue exist [Eberlein et al., 2001], these were used for first in-silico studies on the disc mechanics. With the high fluid content of the nucleus pulposus (NP) and AF-tissues comes a pronounced viscous behavior [Ambard and Cherblanc, 2009], such that elastic equations fail to represent more than the (quasi-) static response. Rubin and Bodner [2002, 2004] presented anisotropic elasto-viscoplastic equations that were successfully applied to the superficial musculo aponeurotic system (SMAS) and to facial skin. For their ability to i

Arruda and Boyce [1993]; Ehret et al. [2010]; Elliott and Setton [2000]; Holzapfel et al. [2000]; Merodio and Ogden [2005]; Rubin and Bodner [2002]; Weiss et al. [1996] ii Elliott and Setton [2000]; Holzapfel et al. [2000]; Weiss et al. [1996]

represent time dependent behavior these equations were finally chosen for the modeling of the AF tissue. Chapter 4 explores a non physical effect observed in simulations with an implementation of the equations proposed by Holzapfel et al. [2000]. In fact, simulations involving uniaxial stress configurations in transversely isotropic materials reveal volume growth at rather small stretches. The same phenomenon is also found with other constitutive equations, too. Investigations led to the conclusion that the additive split of the Helmholz free-energy into a volumetric and an isochoric part, applied to the matrix and the fiber contribution, is the origin of this shortcoming. It is shown that this shortcoming exists independently of the polyconvexity of the models in concern. A solution was found to avoid this problem on the constitutive level rather than to use numerical methods such as the Augmented Lagrangian method [Zienkiewicz et al., 2005]. These findings were presented at the Society of Engineering Science (SES) conference (2008) and published in Helfenstein et al. [2010]

Experimental work A characterization of the overall mechanical properties of the spinal column might be done on the level of the functional spinal unit (FSU), the unit element of the spinal column. A more detailed characterization will try to capture not only the collective mechanical behavior but the properties of the single components that are present in the FSU. One well studied component is the AF and its lamellae.i Under axial load the two vertebrae of an FSU approach, causing the NP to extend radially, what causes the AF to bulge. As this happens, the AF circumference enlarges and the tissue gets stretched in hoop direction. In experiments, also for other motions extensional strains in circumferential direction and in the direction of the collagen fibers were found [Shah et al., 1978; Stokes, 1987]. Though, tension is an important loading mode of the AF-tissue. Nevertheless, its complete mechanical characterization requires data not only from experiments in fiber direction but should include biaxial experiments as well [Bass et al., 2004]. Consequently the practicability of biaxial experiments with AF tissue was studied, especially the aspect of the samples geometry came under scrutiny. Chapter 5 reports the results of a collaboration with M. Hollenstein where a numerical method was developed which maximizes the area of hom*ogeneous equibiaxial loading in isotropic test specimens by following the idea of M¨onch and Galster [1963]. For their ease of handling cruciform specimen geometries are often used but so far little research was done in the soft tissue community to improve the specimen geometry with respect to the materials parameter determination. The findings suggest that already a low number of slots in the limbs of the cruciform samples helps to increase the area of i

Adams and Green [1993]; Ebara et al. [1996]; Elliott and Setton [2001]; Fujita et al. [1997]; Galante [1967]; Green et al. [1993]; Hirsch and Galante [1967]; Holzapfel et al. [2005]; Panagiotacopulos et al. [1979, 1987]; Skaggs et al. [1994]; Wagner and Lotz [2004]; Wu and Yao [1976]

CHAPTER 1. INTRODUCTION

hom*ogeneous equibiaxial loading significantly. The in-silico results are in good agreement with experimental findings obtained with casted silicone samples. This subject was presented at the European Conference on Constitutive Models for Rubber (ECCMR, 2009) [Helfenstein et al., 2009a]. The above mentioned outcomes highlight, with regard to the small size of the possible AF specimens and their anisotropy, the difficulties to generate adequate biaxial stress fields for a faithful material parameters determination. In consequence only (quasi-) uniaxial experiments with AF samples were carried out. Chapter 6 presents the uniaxial tensile experiments with AF-strips that were executed with the aim to calibrate anisotropic, elasto-viscoplastic constitutive equations. Several prestudies were made before the actual experiments started. The first analyzes to which extent, in dependence of the samples aspect ratio, the quasi-uniaxial experiments can be approximated by a purely uniaxial model. A second estimates the error in the strain measurement due to a misalignment of the sample. The compliance of the used force transducer was quantified and the effect of freeze storage on the mechanical properties was qualified. Eventually, circumferential AF-strips were cut from fresh ovine lumbar spines. The experiments include displacement driven load ramps with intermediate relaxation phases at defined strain levels. Results show the presence of an initial toe-region, a distinct stiffening and give insights into the time-dependent material behavior. The data from these uniaxial tensile experiments are used to calibrate elasto-viscoplastic equations [Rubin and Bodner, 2002] for the use in a finite element (FE) model of a FSU. This model would need to be validated with experimental data of the overall responses of (ovine) FSUs to axial loading and bending cycles. Only, in literature few data is available on the response of FSUs subjected to cyclic loadingi , none of them covering ovine samples. Therefore own cyclic axial loading and bending experiments with ovine FSUs were designed and executed, which has the additional benefit to have access on the exact geometries of the characterized samples for the FE model generation. Chapter 7 depicts the challenges met during the realization of the FSU experiments and the present results. Preliminary experiments with steel springs had the objective to demonstrate the performance of the constructed mountings in different load cases and to verify the suitability of the used materials testing machine. These experiments yielded improvements of the setup that would not have been realized without them. The results of this design and verification process were presented in a young researcher minisymposium at the Annual Meeting of the International Association of Applied Mathematics and Mechanics (GAMM) [Helfenstein et al., 2009b]. i

Boxberger et al. [2009]; Chow et al. [2004]; Crisco et al. [2007]; Hansson et al. [1987]; Izambert et al. [2003]; Koeller et al. [1984]; Liu et al. [1983]

For the experiments whole lumbar spines were dissected into FSUs, freeze stored until testing and finally characterized using a standard uniaxial materials testing machine. The experimental protocols comprised axial loading and bending cycles with different load rates to give insight into the viscoelastic properties of the FSUs. Cycles were rerun several times in order to check for repeatability. The results of these experiments reveal typical characteristics of viscoelastic materials as are nonlinearity, distinct dissipation and ratcheting.

Simulations Insights into the mechanics of biological systems are mainly obtained by appropriate experiments. But sometimes experimentation is not possible due to prohibitive ethical, technical or other reasons; if this is the case in-silico studies might present an alternative. Modeling techniques, and in particular the FE method, provide powerful tools for the study of the IVD mechanics. Belytschko et al. [1974] were the first to present a FE model of the disc. In their axisymmetric model six tissues with different linear orthotropic properties are distinguished. Ever since models gained complexity, more and more geometric details were included and the tissues are represented with more elaborate constitutive equations. Anyhow, quality and accuracy of results obtained from such models depend on the underlying assumptions. Model validations are therefore a necessary condition before results can be trusted or even model based predictions are made. Purely phenomenological models [e.g. Groth and Granata, 2008; Li et al., 1995] may well reproduce motions for which they were calibrated but fail for any others. Bottom-up models capture not only the systems collective mechanical behavior but, in contrast, the properties of the single components present in the FSU. Such models have far better predictive capabilities and should, whenever possible, be favored. In the present work a component based FE model was created where the constitutive equations that describe the AF are calibrated with experimental data. Chapter 8 specifies the calibration process of elasto-viscoplastic equations [Rubin and Bodner, 2002] with the data obtained from uniaxial tensile experiments with AF-strips. A stepwise calibration procedure was chosen in order to structure the process and to analyze the importance of the different aspects of the equations. In a first calibration step only the elastic part of the model is used, the second includes rate-independent viscoplasticity and the third adds rate-dependence. As no suitable measurement data for the first two steps was available, reference curves were extracted best possible from the experimental data. Results partially negative the arguments for the stepwise procedure: parameter sets from preceding steps needed to be recalibrated in next steps. The parameter sets found for the elastic and rate-independent part are quite promising while for the rate-dependent model the viscous behavior is not reproduced to full satisfaction. Chapter 9 describes the construction of the FE model of the FSU and presents the results of a validation with data from axial loading and bending experiments. 5

CHAPTER 1. INTRODUCTION

Serial computer tomography scans of an ovine FSU were used to reconstruct the original geometry. Vertebrae are represented in the model as rigid bodies and only the IVD is deformable. The disc is further subdivided into the NP which is, according to the literature, modeled as an incompressible, inviscid fluid [e.g. Eberlein et al., 2001] and the AF. This latter is described by elasto-viscoplastic constitutive equations [Rubin and Bodner, 2002] that are calibrated with data from uniaxial tensile experiments with AF strips. A validation indicates that the model behaves too stiff, especially at larger deformations. It is shown that the NP plays an important role in the activation of the fibers in the AF. The encountered discrepancies between the experimental and the in-silico results could be numerical artifacts but it rather seems they originate from the exponential form of the constitutive equations used for the AF. Chapter 10 summarizes eventually the achievements of the present work and raises open questions that need future investigation.

Chapter

Spine Anatomy This chapter presents the anatomy of the (human) spinal column and its constituents. As the main skeletal structure, the spinal column is responsible for the stabilization of the torso and the protection of the spinal cord. Its double-S-shape and the alternating structure of vertebrae and intervening discs renders it well suited to damp impacts on the body. The intervertebral discs are viscoelastic structures with a high water content. Their inner part, the nucleus pulposus, is made up by an incompressible isotropic fluid like material. Under compression it transforms the axial load into circumferential loads in the anulus fibrosus. This latter structure surrounds the nucleus and is, due to its fiber reinforcement, adapted to carry such loads. In-vivo , muscles and ligaments prestress the spinal column and assure therewith its stability. Studies on the structure of quadruped vertebrae suggest even higher preloads in the animal spine. Despite some differences in the exact geometry and composition, compared to the human spine, ovine models seem to be valid alternatives for (qualitative) mechanical studies.

As already their name vertebrate indicates, the main skeletal structure of their bodies is the spinal column (SC) with its vertebrae (Ve) and intervertebral discs (IVDs). Supported by muscles and ligaments the SC stabilizes the torso while not preventing it from being flexible. Besides stabilization, the SC absorbs impacts that occur under physiological conditions as running or walking. The following sections detail the anatomical and functional descriptions of the human SC and its constituents. Although some anatomical drawings are shown, the interested reader might take a look at the excellent drawings in Sch¨ unke et al. [2005]. In the present work experiments with ovine spinal parts were done and consequently the last section of this chapter will provide a comparison between the human and ovine spine. 7

CHAPTER 2. SPINE ANATOMY

2.1

The spinal column

The basic functions of the SC are to stabilize the vertebrates torso, to absorb mechanical vibrations and shocks and to protect the spinal cord [Keller et al., 1987]. The human SC is divided into five main levels, which are, from top to bottom, (cf. Fig. 2.1a): -

cervical spine (pars cervicalis) CI-CVII, thoracic spine (pars thoracalis) ThI-ThXII, lumbar spine (pars lumbalis) LI-LV, sacral spine (os sacrum) SI-SV and tail bone (os coccygis) with four to five Ve.

In the first three spinal levels Ve alternate with IVDs (except for the first and second cranial Ve), while the Ve of the sacral spine and the tail bone are grown together and no IVDs are present. Even though the individual Ve of the different levels are single in their shape and size, the general layout is the same [Sch¨ unke et al., 2005]. The size of the Ve increases from top to bottom since the body weight that lasts on them increases. Each V consists of a body (corpus vertebrae) to which a bony neural arch (arcus vertebrae) is attached, see Figure 2.2. Diverse processes are attached on the outer side of the neural arch that act as point of load incidence of muscles and ligaments. The vertebral body is not a compact bone but a bony shell reinforced by thin bony rods, the so called trabeculae. This inner part is vascularized and filled with blood, such that it resembles a sponge; hence it is sometimes called vertebral spongiosa [Bogduk, 2005]. Designed to support axial forces the vertebral body has nearly flat top and bottom surfaces (cf. Fig. 2.2). The series of neural arches in the SC form the so called vertebral canal (canalis vertebralis) that surrounds and protects the spinal cord. The segmental structure of the SC allows for complex movements of the torso without interference or even damage to the spinal cord. In contrast, a common impairment develops if a part of the IVD bulges pathologically and compresses the spinal cord, a so called disc hernia. Since all nerves from extremities end in the spinal cord, a compression may end in pain that can radiate all-over the body. Due to their flat surfaces two vertebral bodies cannot withstand any force other than axial. To stabilize the relative positions of two neighboring Ve they are guided by the zygapophysial joints (articulationes zygapophysiales) formed by the superior and inferior articular processes (processus articularis sup. and inf., respectively) of two neighboring Ve. In the lumbar spine these synovial joints resist forward sliding and twisting of the Ve. While the zygapophysial joints mainly determine the possible directions of motion, the IVDs and ligaments limit the motion in its range (more details about the range of motion is given in Sec. 2.4). The four basic movements that the SC allows are flexion (forward bending), extension (backward bending), lateral flexion and axial torsion (cf. as well Fig. I.1 in the nomenclature).

2.1. THE SPINAL COLUMN CI anulus fibrosus with crossbred fibers

CVII ThI vertebra

vertebra

@ @ @ @ @ @ @

intervertebral disc

(b) frontal view on a functional spinal unit

ThXII

anulus fibrosus

vertebra @ @ @ @ @

nucleus pulposus

PP PP PP endplate

LV (c) cut through a functional spinal unit

sacral spine

tail bone (a) human spinal column

Figure 2.1: View on an entire human spinal column with cervical (CI-CVII), thoracic (ThI-ThXII) and lumbar (LI-LV) vertebrae (a). Frontal view on a lumbar functional spinal unit (b). Cut through a functional spinal unit (c). Illustrations from Sch¨ unke et al. [2005] 9

CHAPTER 2. SPINE ANATOMY arcus vertebrae proc. articularis superior proc. mammillaris HH HH proc. accessorius HH @ proc. costalis @ HH @ H proc. spinosus @ @ HH @ @ H proc. articularis inferior @ @ @ @ @ @ corpus vertebrae @ foramen vertebrale @

Figure 2.2: Fourth lumbar vertebra (LIV) and its processes. Illustrations from Sch¨ unke et al. [2005]

The capacity of the SC to absorb mechanical vibrations and shocks bases on two of its characteristics. The first is the double-S-shape of the SC in the sagittal plane (cf. Fig. 2.1a) that eases its bending under axial load [Sch¨ unke et al., 2005]. The second is the constitution of the SC with its alternating structure of Ve and IVDs. The latter have, due to their viscous nature (cf. Sec. 2.2) the ability to dissipate mechanical energy [Hirsch, 1955]. Vibrations and shocks, as they arise e.g. during walking or running, are effectively damped. Two consecutive Ve, the intervening IVD, the zygapophysial joints and the muscles and ligaments that run between the two Ve form a functional spinal unit. For what follows the term functional spinal unit is defined slightly different, muscles and ligaments will not be considered to be part of it (cf. Fig. 2.1b). In an functional spinal unit three joints are present: the posterior zygapophysial joints (two) and the the interbody joint formed by the IVD.

2.2

The intervertebral disc

The IVD forms an interbody joint that allows for relative motions of two neighboring Ve (even though the effective mobility of the two Ve is limited by the zygapophysial joints). One third to nearly a half of the SC height is made up by the IVDs [De Puky, 1935; Virgin, 1951]. The structure of the IVD is designed to respect three basic demands [Johannessen et al., 2004]: - load support, - flexibility and - energy dissipation. On top and bottom each IVD is covered by 0.1 – 1.6 mm thick cartilagineous vertebral endplates [Roberts et al., 1989] that separate the IVD from the adjacent vertebral bodies (cf. Fig. 2.1c). Between the two endplates lays the nucleus pulposus (NP) that is 10

2.2. THE INTERVERTEBRAL DISC

circumferentially surrounded by the anulus fibrosus (AF). Between NP and AF no clear boundary exists [Markolf and Morris, 1974], they rather merge continuously from one to the other. All three structures, vertebral endplates, NP and AF have basically the same constituents and differ mainly in their composition. The basic constituent of the three is water with concentrations ranging from 60 – 90 % (AF and NP, respectively). The high water content influences the elastic properties of the IVD [Andersson and Schultz, 1979] and its viscoelastic properties are largely dependant on the ability to absorb and lose fluid [Virgin, 1951]. In fact, the fluid content is not constant but depends on the mechanical load on the IVD (e.g. muscle forces, body weight, external loads). The fluid content is always defined by an osmotic equilibrium between these loads and the swelling pressure in the tissue [Johnstone et al., 1992]. The solid constituents are proteoglycans and collagen. These two components stabilize the tissues by their interconnections via electrostatic and covalent bindings. Proteoglycans are very large molecules made up by a lot of glycosaminoclycans interconnected by proteins. As in articular cartilage, the essential proteoglycan of the IVD is aggrekan. In vivo this proteoglycans convolute around each other to form extensive three-dimensional complexes that own hydrophile and hydrophobic regions maintaining the high fluid content of the IVD. Seven types of collagen are known to be present in the IVD (I, II, III, V, VI, IX and XI) of which types-I and -II account for about 80 % by weight [Eyre, 1988]. Type-I collagen is very expansible and mostly found in tissues that are stretched or compressed. Type-II collagen is elastic and often found in tissues that withstands compressive loads. Old proteoglycan and collagen are systematically removed and replaced. Since the IVD is an avascular structure (at least after the age of 25 years [Nachemson et al., 1970]), nutrition of cells that produce new components is done by dint of diffusioni . Nutrients are thought to be transported through the porous vertebral endplates (and the AF) to the inner parts of the IVD. The endplates are in contact with the vascular supply [Roberts et al., 1989], by their small thickness they should not act as a significant diffusion barrier. Far more important are the dimensions of the IVD that has a height of some millimeters such that nutrients must diffuse a long way to reach the center part. Constant changes in posture are important to hydrate and dehydrate the IVD [Wilke et al., 2001], what might be positive for nutrition [Albro et al., 2008]. The nucleus pulposus The NP is an oval gelatinous part that occupies the center part of the IVD. Located closer to the posterior than to the anterior border of the disc, it occupies 50 – 60 % of its cross-sectional area [Markolf and Morris, 1974]. The NP is, for young and healthy adults, a viscous material made up by some cartilage cells and randomly arranged collagen fibers embedded in a base substance. Even though the content varies with age, the NP embodies 70 – 90 % (weight) of waterii . With 65 % dry weight proteoglycans are the second common components i ii

Brodin [1954]; Huber et al. [2007]; Nachemson et al. [1970] Gower and Pedrini [1969]; Markolf and Morris [1974]; Roberts et al. [1989]

CHAPTER 2. SPINE ANATOMY

[Gower and Pedrini, 1969] followed by collagen (especially type-II) which adds 6 – 25 % to the dry weight of the NP [Eyre, 1988]. Its high fluid content renders the NP fluid like and incompressible (at short time scales) while proteoglycans and collagen add viscosity and density (cf. as well Sec. 2.3). This high fluid content is coupled with the content of proteoglycan and therefore a reduction of proteoglycan, as with increasing age [Gower and Pedrini, 1969], lowers the load bearing capability of the IVD. It was shown experimentally that a disc with removed NP has the same short time strength under axial load as an intact one [Markolf and Morris, 1974]. What is different is the long term response: liquid is squeezed out of the AF and it buckles, therefore the importance of the NP. As the NP is incompressible by its high water content it will, under axial load, reduce its height by an increase of its lateral dimensions. This radial extension bulges the AF radially and induces therein circumferential stresses that withstand the axial load. Additionally, the radial pressure from the NP prevents the AF from buckling such that it retains its axial strength. An IVD under an axial load of 40 kg will become 1 mm thinner and extends radially by 0.5 mm [Hirsch and Nachemson, 1954]. The anulus fibrosus The AF is a stiff structure that runs “concentrically” around the NP (cf. Fig. 2.1c). It is build up by about 20 – 90 layers [Marchand and Ahmed, 1990; Panagiotacopulos et al., 1987], called lamellae, where each consists of a transversely isotropic material. The ground substance contains mainly water (around 70 – 80 % of its weighti ), proteoglycans (around 20 % of the dry weight) and non-collagenous proteins. The concentration of proteoglycans increases from outer to inner regions [Best et al., 1994]. This matrix is reinforced by collagen fibers of types I and II. In total the collagen makes 50 – 60 % of the dry weight or about 16 % of the volume of the AF [Galante, 1967]. The type-I collagen content increases in the AF from 0 % in the transition zone of NP and AF up to 80 % in the outer AF. Collagen type-II is found all over the AF with increasing content from 0 % in the outer rim to 80 % in the transition zone [Eyre, 1988; Roberts et al., 1991]. The AF contains as well elastic fibers (elastin) that run circular, diagonal or vertical inside the lamellae [Johnson et al., 1982]. In each layer collagenous fibers run in parallel from the lower V to the upper giving the lamellae their transversal isotropic properties. According to Cassidy et al. [1990] the fibers are not fully stretched but possess a planar crimped waveform. With the vertebral endplates the fibers make an angle of ±30 – 35 ◦ [Hickey and Hukins, 1980a,b; Holzapfel et al., 2005]. The sign of the angle alternates from one to the next lamellae (cf. Fig. 2.3). The angle is not constant, but varies in radial and circumferential directions. Cassidy et al. [1989] find that the lamellar fiber angle increases from 28◦ at the outer rim to 45◦ beside the NP. Marchand and Ahmed [1990] report varying fiber angles from 20◦ at the ventral position to 70◦ at the dorsal position. According to Holzapfel et al. [2005] the fiber angle θ is described approximately by the function |θ| = 23.2 + 0.13α [◦ ], where α i

Eyre [1988]; Gower and Pedrini [1969]; Roberts et al. [1989]

2.2. THE INTERVERTEBRAL DISC

is the polar angle associated with the circumferential position (α = 0 [◦ ] for a medial ventral position and α = 180 [◦ ] for a medial dorsal position). In contrast to Cassidy et al. [1989], Holzapfel et al. [2005] don’t find a strong dependence of the fiber angle on the radial position. By the same authors it is shown that the absolute values of positive and negative oriented fibers do not differ. The collagen fibers of the inner lamella enter the endplates and run to their center, a fact that might help the endplates to carry the NP pressure [Inoue and Takeda, 1975; Roberts et al., 1989]. Though the NP is more or less encapsulated by collagen fibers. The fibers of the outer lamella, that are not covered by the endplates, are anchored in the bone of the vertebral body. Although it is suggested by the above description (and shown in Fig. 2.3) that each lamella runs completely around the IVD, this is not the case. Around 40 – 50 % of all lamellae are incomplete and not closed [Marchand and Ahmed, 1990]. Neither all lamellae have the same thickness nor a single lamella has a constant thickness. The lamellae get thicker to the center [Marchand and Ahmed, 1990] and lamellae are around two times thinner at posterior sites than at anterior or lateral [Markolf and Morris, 1974]. Holzapfel et al. [2005] report lamellae thicknesses of 0.69±0.07 mm and 0.76±0.10 mm at outer and inner ventro-lateral sites, respectively. At dorsal sites the lamellae thickness varies from 0.38±0.04 mm to 0.40±0.06 mm from outer to inner positions. As already mentioned in Section 2.2, the AF bulges under an axial load causing circumferential stretches in the lamellae. Shah et al. [1978] and Stokes [1987] showed experimentally for various loads the presence of extensional strains in circumferential direction and in direction of the fiber bundles, respectively. Osmotic effects result in additional tensile stresses in the tissue [Laible et al., 1993]. Though it seems that tension is an important loading mode of the AF. As a conclusion Holzapfel et al. [2005] suggest the single lamella as the elementary structural unit of the AF. It is worth to think more about the function of the collagen fibers in the lamellae and why they are oriented as they are. If the fibers would be oriented perpendicular to the vertebral endplates they would provide maximum strength in tension. But on the other hand they would offer minimal strength against relative sliding of two Ve. Mathematical calculations showed that a 35◦ orientation is optimal for the AF to withstand different deformations [Hickey and Hukins, 1980b]. The alternating orientation of the collagen fibers allows the IVD to resist axial rotations in both directions.

Figure 2.3: The laminate structure of the AF: fibers include alternatingly the angle ±θ with the transversal plane. After Bogduk [2005]

CHAPTER 2. SPINE ANATOMY

2.3

Loading of the intervertebral disc

The SC has to withstand most of the forces acting on the vertebrate body but still needs flexibility to allow for changes in posture. Flexibility of the SC is given by the deformability of the IVDs, the mobility of the facet joints and compressibility of the motion segments [Brinckmann et al., 1983]. Build up by an alternating sequence of vertebral bodies interconnected by ligaments and discs, the intrinsic stability of the SC is low. Nachemson [1966b] reports a buckling load of about 20 N if applied to the first thoracic V of an erect standing dissected spine. Though it is clear that the SC alone would not have the capability to withstand the forces acting on it in daily live and therefore has to be stabilized by external means as are the paraspinal and trunk muscles. In-vivo, ligaments, muscles and external loading keep the discs under permanent load [Urban and Maroudas, 1981] that counterbalances the swelling pressure of the discs [Pflaster et al., 1997]. The most obvious load on the SC is the body weight. Ruff [1950] determined the percentage of the body weight acting on the lumbar spine to be 55.4 % (average over the five lumbar Ve). This (gravitational) load cannot entirely explain the intradiscal pressure found by Nachemson [1966b]. Electromyographic studies revealed that the psoas muscle helps maintaining an upright posture and adds a compressive force to the lumbar spine [Nachemson, 1966b]. A minor role have the interlaminar ligaments (ligamenta flava) (with around 0.15 MPa) which are the only posterior elements that prestress the discs [Nachemson, 1966a; Nachemson and Evans, 1968]. In all postures the IVDs carry most of the load acting on the SC. In erect standing only 16 % of the total load is carried by the apophysial joints; most of the force is taken by the IVDsi . Nachemson was the first measuring the load applied to the disc in-vivo. A needle connected to a pressure transducer was introduced to the inside of lumbar discs [e.g. Nachemson, 1992]. With this procedure the fluid like behavior of the NP could be shown and the load in different postures was measured. Wilke et al. [2001] found NP pressures being 0.1 MPa in a lying and 0.5 MPa in a relaxed standing person; especially they found the nuclear pressure being smaller in sitting than in standing position. Nachemson [1966a] found that the pressure in the NP is about 50 % greater than the externally applied vertical load and does increase linearly for external loads up to 2 kN. Furthermore the load on the lumbar spine is (in sitting position) approximately three times the body weight above the given lumbar level. From this global load values Nachemson [1966a] estimates the load on the AF by means of simple deformable body mechanics. Only half of the external vertical load applied on the SC is supported by the AF while the other is carried by the NP. Tensile stresses in the AF depend on its thickness but seem to be 4 – 5 times the applied external load per unit area. Only few in-vivo studies are available that determine the real-time loading of the spine. Ledet et al. [2000] found in an in-vivo study with baboons peak forces higher i

Adams and Hutton [1980]; Lin et al. [1978]; Nachemson [1963]

2.3. LOADING OF THE INTERVERTEBRAL DISC

than 2.8 times the body weight (of about 40 kg). This and results from other studies [Nachemson, 1966a; Wilke et al., 2001], indicate that peak forces being 1.5 times the body weight are reasonable, not inducing irreversible damage to the tissue [Huber et al., 2007]. While the NP is considered to behave mechanically as an incompressible fluid [Nachemson, 1966a] or a perfectly plastic material [Keller et al., 1987] the AF might be thought as a pressure vessel (cf. as well Sec. 6.1). As the collagen fibers that reinforce the AF are grown into the cartilage endplates they must stretch when the disc bulges radially [Cassidy et al., 1990]. Calculations show that if the fiber tilt is at most 35.3◦ with the transverse plane fibers get stretched when the AF bulges. Though the fibers of the AF are loaded in tension and their mode of deformation is comparable to that of fibers in tendons or ligaments. The fiber stretch balances pressure in the NP and therefore the axial load [Hickey and Hukins, 1980b; Klein et al., 1982]. By X-ray diffraction it was found that a compression of the IVD lowers the tilt of the collagen fibers, indeed [Klein and Hukins, 1982]. Axial stress relaxation experiments on IVDs revealed that the initial disc bulge almost completely recovers with time [Cassidy et al., 1990]. Though the volume of the disc reduces; a phenomenon that is explained by a fluid loss of the disc. The higher the axial strain the higher the overall decrease in fluid content. Adams et al. [1990] report an average fluid loss of 12 % in the AF and 5 % in the NP, and a decrease in disc height by 1.5 mm, after 4 h creep with 700 N load. Cassidy et al. [1990] did not register any condensed fluid on the disc surface and conclude that the fluid escapes through the cartilage endplates into the vertebral bodies [also Ayotte et al., 2000]. The same authors make the hypothesis that “the mechanism responsible for the viscoelastic character of the disc response in compression, . . . , is the transport of water out of the disc through the cartilage end-plates into the vertebral bone”. During standing the human SC is under permanent compressive load that induces a hydrostatic pressure in the IVD. This pressure is 3 – 5 times larger than the discs osmotic swelling pressure, which is about 0.15 – 0.2 MPa [Adams et al., 1996]. Ayotte et al. [2001] found that discs lose up to 20 % of their volume during the day. This volume loss is compensated during bed time when the load is reduced and the IVDs absorb fluid due to their osmotic swelling pressure [Adams et al., 1990]. A reduction in disc volume results in a reduction in body length and, in fact, the human body length oscillates in a diurnal rhythm. De Puky [1935] measured an average change in body height of 1.022 % (that was not only assigned to the contribution of the IVD volume change but as well to the curvature of the SC that becomes more pronounced during day). The change in body height is especially pronounced the first minutes after out of bed [De Puky, 1935]. It is to be expected that this diurnal changes in fluid content and disc height have an influence on the mechanical properties of the spine. Adams et al. [1990] found that with creep loading the IVDs become more elastic, stiffer in compression and more compliant in bending and their affinity for water increases.

CHAPTER 2. SPINE ANATOMY

Under an axial load the discs bulge not only radially but axially as well, i.e. the cartilagineous endplates bulge into the spongiosa of the neighboring Ve. Brinckmann et al. [1983] report same orders of magnitude for change in disc height, radial and axial bulge where Adams et al. [1990] found the radial bulging being one-third of the axial displacement. Irreversible axial bulging is the result of an axial overload of the SC, in fact, the vertebral end-plates are the first structures to break [Nachemson, 1992]. Ruff [1950] observed in clinical observations that the lumbar spine can withstand around 10 kN of axial load before failing. In pure compression Brinckmann et al. [1983] found the end-plates rupturing at loads of 7.5 kN after a deflection of around 0.5 mm. Cassidy et al. [1990] found two different failure modes of spine segments under axial load: 1) the stress strain curve suddenly decreases and a “popping or cracking sound” can be heard and 2) the stress strain curve decreases with increasing strain. While the segments failed with mode 1) revealed ruptured end-plates, segments with mode 2) failure showed fissures in the AF through which NP material leaked.

2.4

Comparison of human and ovine vertebrae and intervertebral discs

As human tissue for in-vitro studies is difficult to obtain and in-vivo studies on human subjects are hard to get approved, the question arises if such studies cannot be done using animals and animal tissues. Except for the simpler supply, positive aspects of animal subjects are the smaller variability if taken from the same breed and age and the availability of young subjects. And, last but not least, bacterial or viral infections (e.g. hepatitis or AIDS) are ever-present risks when working with human material, even if the cadavers were screened. It seems obvious that a quadruped spine is under a different load that an upright human spine (cf. Fig. 2.4) and therefore one would expect major differences in the anatomy. Especially lower axial loads would be to expect. But in contrary, Smit [2002] found caprine Ve to have a two times denser trabecular bone than human Ve and argues that, according to Wolff’s law [Wolff, 1986], the axial load in the Ve is higher in quadrupeds. It is suggested that the high axial load is necessary to control the posture of the spine as the spine itself cannot withstand the bending moments exerted by the body weight. If the quadruped spine is mainly loaded by axial loads, its loading conditions are similar to the ones in the upright human spine. In the following special attention is laid on the comparison of the human and ovine lumbar spine since the latter is used in the mechanical experiments effectuated during this thesis. Smeathers [1984] finds “. . . sufficient similarities in the [ovine and human] disc thickness, types and distribution of tissues to justify using them [the ovine spines] to demonstrate the techniques prior to experimenting on the much less available human cadaveric material”.

2.4. COMPARISON OF HUMAN AND OVINE VERTEBRAE AND INTERVERTEBRAL DISCS

Figure 2.4: Human and ovine skeleton. After Nickel et al. [1992]

The ovine spine has some more Ve than the human spine; it consists of 7 cervical, 12 – 14 thoracic and 6 – 7 lumbar Ve [Nickel et al., 1992] (vs. 7, 12 and 5 in the human spine). Human and ovine lumbar Ve have an oval shape where the width is greater than the depth. But while human Ve are approximately two times wider than tall the sheep Ve are taller than wide [Wilke et al., 1997b]. Another difference concerns the facet joints for which in quadrupeds the articular surfaces and the frontal plane include an angle of about 28◦ while in humans this angle is about 58◦ [for cattle, Cotterill et al., 1986]. Easley et al. [2008] found that in the sheep spine the posterior elements carry more load (relative to the total load) than this is the case for young healthy humans. Ovine as well as human lumbar IVDs have a wedge like shape [Easley et al., 2008]. The ovine IVD height is with 5 – 6.3 mm [Easley et al., 2008; Wilke et al., 1997b] thinner than in the human spine where disc heights of 8 – 16 mm were measuredi . This smaller disc height in sheep is in accordance with the smaller anterior-posterior and mediallateral dimensions of the vertebral body [Easley et al., 2008]. Wilke et al. [1997a] compare the ranges of motion of ovine and human spine levels and conclude that the sheep spine might be used for the evaluation of spinal implants. Table 2.1 summarizes the range of motion values found in the cited study for the lumbar level of the ovine spine. The range of motion for flexion-extension, lateral bending and axial rotation at the LIV-LV level are comparable in sheep and human. Easley et al. [2008] find good agreement in flexion-extension, lateral bending and axial rotation for the ovine spine (levels LIII-LIV and LIV-LV) with in-vitro measurements on human lumbar spines but not with in-vivo measurements. This differences are attributed to a strong in-vivo muscle activity that generates larger ranges of motion. i

Campbell-Kyureghyan et al. [2005]; Costi et al. [2002]; Nachemson et al. [1979]; Pfirrmann et al. [2006]

CHAPTER 2. SPINE ANATOMY Table 2.1: Range of motion [◦ ] of the ovine lumbar spine [after Wilke et al., 1997a]. In all four tested cases a pure moment of 7.5 N m was applied. flexion L1 L2 L3 L4 L5 L6

– – – – – –

L2 L3 L4 L5 L6 L7

3.99 3.97 4.02 4.07 3.24 5.29

±0.64 ±0.73 ±0.91 ±1.40 ±0.49 ±0.82

extension 5.70 5.60 5.08 4.57 3.89 5.72

±0.66 ±0.55 ±0.71 ±0.94 ±0.47 ±0.47

axial rotation 1.31 0.95 0.73 0.73 0.59 0.62

±0.18 ±0.19 ±0.12 ±0.20 ±0.18 ±0.15

lateral bending 6.20 5.02 4.40 4.29 4.31 4.14

±0.91 ±0.77 ±0.43 ±0.80 ±0.69 ±1.07

Along with this comparisons of the overall anatomy and mechanics only few data is available on the biochemistry and structure of the ovine IVD in special and the IVD of other quadrupeds in general. Meakin and Hukins [2000] report a water content of about 81±2 % in the ovine NP which is comparable with the range of 70 – 90 % reported for the human NP (cf. Sec. 2.2). Cassidy et al. [1990] find a similar hierarchical structure in human and canine lumbar discs but no details are given. Huber et al. [2007] believe that their qualitative findings might apply to human tissues as well, due to the “structural similarities between human and sheep discs”.

Chapter

Nonlinear Continuum Mechanics: a Short Introduction The theory of continuum mechanics presents a framework to treat the interplay of loads and deformations acting on a body. In most processes in nature and, especially, in biological systems, this interplay is highly nonlinear. Nonlinearities might arise from different sources which are 1) large strains, i.e. large relative displacements, 2) a nonlinear relationship of strains and forces and 3) boundary conditions. Nonlinear continuum mechanics provides and establishes a theory which allows a matching mathematical treatment of such systems.

This chapter gives a brief overview into the theory of nonlinear continuum mechanics. It is designed to present the necessary quantities and objects that will be used in the subsequent chapters. Far from being a full-fledged introduction into the depth of this field, the reader is referred to the classic text-books of (nonlinear) continuum mechanics [e.g. Ogden, 1997]. The notation used in the present work is inspired by the one used in Holzapfel [2000], a textbook which gives a good introduction into tensor analysis and solid mechanics.

3.1

Kinematics

Consider a body Ω0 in the (undeformed) reference configuration with material points P ∈ Ω0 . In a Cartesian coordinate system with basisvectors ei , i = 1, . . . , 3, each point P is described by its coordinates X (cf. Fig. 3.1). Under an imposed load Ω0 experiences a deformation and maps into the deformed (current) configuration Ω where the material points p are described by vectors x[X; t]. 19

CHAPTER 3. NONLINEAR CONTINUUM MECHANICS ϕ[X; t]

Ω

Ω0 U [X; t]

e3 x[X; t]

X e2

e1 Figure 3.1: Reference configuration of a body Ω0 and its mapping into a deformed configuration Ω.

Displacement The transition from the reference to the current configuration is described by the mapping ϕ[X; t], called motion, such that x = ϕ[X; t]; x is a function of X and the time t. The displacement U [X; t] of a point X is therefore defined as U [X; t] = x[X; t] − X. In the following the function arguments will be omitted if not necessairy for clarity. Deformation gradient We now define the mapping of an infinitesimal line-element dX in the reference configuration into an infinitesimal line-element dx in the deformed configuration ∂ϕ ∂U dx = dX = 1 + dX = FdX (3.1) ∂X ∂X where 1 is the second order identity and F is called the deformation gradient. F can be decomposed in F = RU = vR (3.2) where R is a proper orthogonal tensor describing a rotation and U and v are the right and left stretch tensors, respectively. Thus, F describes the stretch and rotation of an infinitesimal line element dX. Volume ratio The determinant of the deformation gradient F is the volume ratio J, i.e. the ratio between the volume of a deformed infinitesimal volume element and its volume in the reference configuration: J = det[F] = dV /dV0 . (3.3) An (near) incompressible material requires that J stays (close to) one for all deformations. 20

3.2. STRESS MEASURES

Stretch measures Since only stretches (but not pure rotations) contribute to the deformation energy a body takes under a superimposed load, it is customary to define a right and left Cauchy-Green deformation tensor C and b, respectively, that depend only on the stretches (cf. Eq. 3.2) C := FT F = UT U,

b := FFT = vvT .

(3.4a,b)

Principal stretches and directions The deformation gradient F can be written in terms of its eigenvalues λi , the principal stretches, and -vectors ni and N i , i = 1, . . . , 3 F=

3 X i=1

λi ni ⊗ N i .

(3.5)

The N i are vectors defined in the reference configuration while the ni are vectors defined in the deformed configuration. It follows directly for the right and left Cauchy-Green deformation tensors

3 X i=1

λ2i N i

⊗ N i,

3 X i=1

λ2i ni ⊗ ni

(3.6a,b)

from where can be seen that the eigenvalues of C and b are the squares of the eigenvalues λi of the deformation gradient F. Often materials tests, e.g. uniaxial tensile tests, result in strain/stress states where the directions of the principal axes are knows. It is customary to describe such tests in terms of the principal strain/stress values. The deformation gradient F, the right and left Cauchy-Green deformation tensor C and b, respectively, might be split into an isochoric (distortional) and a volumetric part, ¯ and (•)vol , respectively [Flory, 1961]: (•) ¯ = J 1/3 1F, ¯ F = Fvol F

¯ = J 2/3 1C, ¯ C = Cvol C

¯ b = bvol b¯ = J 2/3 1b.

(3.7a–c)

It can easily be demonstrated that the determinants of the deviators equal unity for all deformations ¯ = det[C] ¯ = det[b] ¯ ≡ 1. det[F] (3.8)

3.2

Stress measures (i)

Imagine a body Ω in a (deformed) configuration where external loads F ext are applied (cf. Fig. 3.2a). If one virtually cuts that body into two parts, Ω(1) and Ω(2) , two cut(1) (2) ting surfaces ∂Ω(1) and ∂Ω(2) can be defined. Internal forces F int and F int such that (1) (2) F int = −F int are necessary to hold the two parts together (Fig. 3.2a). 21

CHAPTER 3. NONLINEAR CONTINUUM MECHANICS

(3)

F ext

F ext Ω(2) (2) ds

Ω(2)

t(2) ∂Ω(2)

∂Ω(2)

(2) F int

n(2)

(1)

F int

n(1) t(1) ds(1)

∂Ω(1)

(2)

Ω(1)

F ext

(1)

(2)

Ω(1)

F ext

(1)

F ext

(a)

(b)

Figure 3.2: Virtual cut of a body Ω into two parts Ω(1) and Ω(2) . (a) Introduction of internal (1) (2) forces F int and F int . (b) Introduction of traction vectors t(1) and t(2) .

On the cutting surfaces ∂Ω(i) infinitesimal surface areas ds(i) with unit outward normals n(i) can be defined (Fig. 3.2b). On each surface area a traction t(i) acts that has the physical dimension of force per area. Here again holds n(1) = −n(2) and t(1) = −t(2) . The integration of all traction vectors over the cutting surface gives the internal forces: Z (i) F int = t(i) ds(i) . (3.9) ∂Ω(i)

The tractions t defined in the current configuration are called Cauchy traction vectors. Pseudo traction vectors T , called first Piola-Kirchhoff traction vectors, measure the current force per unit area dS in the reference configuration. The tractions depend on the position, outward normal and time t = t[x, n, t],

T = T [X, N , t].

(3.10a,b)

Cauchy’s stress theorem says that unique second-order tenors σ and P exist such that t = σn,

T = PN .

(3.11a,b)

σ is called the Cauchy (or true) stress and P the first Piola-Kirchhoff (or nominal) stress. By using the balance law of the angular momentum one can show that σ is symmetric (while P is in general asymmetric but satisfies the condition PFT = FPT ). Nanson’s formula, ds = JF−T dS, which connects surface elements in different configurations (and especially surface elements in the current and the reference configuration) yields a relation between the two stress measures: P = JσF−T . 22

(3.12)

3.3. CONSTITUTIVE RELATIONS

The so called second Piola-Kirchhoff stress S has no direct interpretation but proves to be very useful for many applications. Its relation to the previously introduced stress measures is given by S = JF−1 σF−T ,

3.3

S = F−1 P.

(3.13a,b)

Constitutive relations

Constitutive equations describe the relation between the strains undergone by a certain material and the resulting stress. Linear elasticity The most simple and best known constitutive relation is the one for isotropic, linear elastic materials. Hook’s law describes a linear relationship between the deformations (strains) ε and the true stress σ: σ = 2µε + λtr[ε]1. The strain ε is defined as ε=

1 2

∇u + (∇u)

(3.14)

(3.15)

The Lam´e constants µ [MPa] and λ [MPa] are given in terms of the (more classic) Young’s modulus E [MPa] and Poisson’s ratio ν [-]: µ=

3.3.1

1 E , 21+ν

λ=

Eν . (1 − 2ν) (1 + ν)

(3.16a,b)

Green-elasticity for isotropic materials

Elastic materials that deform with negligible energy dissipation are called Green- or hyperelastic. For such materials the stress can be deduced from a strain energy density Ψ (the Helmholz free-energy per unit reference volume) by taking its derivative with respect to the deformation (for notational ease, here and in the following, different strain energy functions will be written with the same symbol Ψ): P=

∂Ψ[F] , ∂F

S=2

∂Ψ[C] . ∂C

(3.17a,b)

By the fact that Ψ should not depend on superimposed rigid body motions (material frame indifference) it is possible to show that it must depend on the right stretch tensor U (cf. Eq. (3.2)), only. Typically, Ψ is written as a function of the right Cauchy-Green deformation tensor C which depends on U only (cf. Eq. (3.4a)): Ψ = Ψ[U[F]] = Ψ[C[F]].

(3.18) 23

CHAPTER 3. NONLINEAR CONTINUUM MECHANICS

The invariance of Ψ under transformations that respect the material symmetries (in case of isotropic materials this are arbitrary rotations) leads to a representation that depends only on the so called invariants of C or b: Ψ = Ψ[I1 [C], I2 [C], I3 [C]] = Ψ[I1 [b], I2 [b], I3 [b]]

(3.19)

with = λ21 + λ22 + λ23

I1 = tr[C] 2

I2 = 21 (tr[C]) − tr[C2 ] = I3 = J 2 = det[C] =

λ21 λ22 + λ21 λ22 λ23

λ21 λ23

(3.20a) +

λ22 λ23

(3.20b) (3.20c)

From Ψ[I1 , I2 , I3 ] the first Piola-Kirchhoff stress P can be obtained by taking the derivative with respect to the deformation gradient P = 2Ψ,1 F + 2Ψ,2 (I1 F − 2FC) + JΨ,3 F−T

(3.21)

where Ψ,i depicts the partial derivative of Ψ with respect to the i’th invariant Ii . By the multiplicative split of the right Cauchy-Green tensor C (cf. Eq. (3.7b)) it is possible to write the strain energy Ψ in an uncoupled form ¯ I¯1 , I¯2 ]. Ψ = U [J] + Ψ[

(3.22)

¯ depends only on Here U is the response of the material to volume changes and Ψ the isochoric part of the deformation. I¯1 and I¯2 are the invariants of the unimodular ¯ of the right Cauchy-Green deformation tensor C: (isochoric) part C ¯ ¯ 2 − tr[C ¯ 2] . I¯1 = tr[C], I¯2 = 1 tr[C] (3.23a,b) 2

The first Piola-Kirchhoff stress reads in this uncoupled form1 ¯ ,1 F − 1 I1 F−T P = JU,J F−T + 2J −2/3 Ψ 3 −4/3 ¯ 2 −T + 2J Ψ,2 I1 F − I2 F − FC

(3.24)

Sansour [2008] shows the additive split to be a consequence of the assumption that the pressure is solely a function of J, i.e. of the volume changes. The split (3.22) enables distinct modeling of high resistance to volumetric deformation and low resistance to isochoric deformation as in the case of nearly incompressible materials. For fully incompressible materials the constitutive pressure cannot be deduced from the constitutive equations but has to be determined using the boundary conditions. In Equation (3.24) the term including U,J has to be replaced by the unknown pressure p. 1

If the second Piola-Kirchhoff stress S is considered the extra terms reveal their nature: ¯ ∂ I¯1 ∂ I¯2 ∂C −1 ¯ ¯ S = JU,J C + Ψ,1 ¯ + Ψ,2 ¯ : ∂C ∂C ∂C ¯1 ∂ I ∂ I¯2 −1 −2/3 ¯ ¯ = JU,J C + J P : Ψ,1 ¯ + Ψ,2 ¯ . ∂C ∂C

P is the projection tensor in the reference configuration which makes the respective stress components deviatoric. This means that the stress components originating from a dependency on I¯1 or I¯2 are independent on a volume change.

3.3. CONSTITUTIVE RELATIONS

Polynomial models The polynomial models are truncated Taylor-series of the isochoric strain energy den¯ around the initial state where I¯1 = 3 [-] and I¯2 = 3 [-] [Mooney, 1940; Rivlin, sity Ψ 1949]: =0

}| { z ¯ ¯ ¯ ¯ ¯ ¯ Ψ[I1 , I2 ] = Ψ[I1 = 3, I2 = 3] ∞ X ¯ 1 ∂ i+j Ψ ¯1 − 3 i I¯2 − 3 j . + I i!j! ∂ I¯1i ∂ I¯2j I¯1 =3,I¯2 =3 i+j=1 {z } |

(3.25)

=cij

¯ at the reference state together with the preceding constants The partial derivatives of Ψ are treated as parameters cij [MPa]: ¯P = Ψ

N X i+j=1

j i cij I¯1 − 3 I¯2 − 3 + U [J].

(3.26)

Here again the volumetric energy U [J] was used. The small strain shear modulus µ equals 2 (c10 + c01 ) in this model. The expression for the first Piola-Kirchhoff stress of the polynomial model PP is given in Table 3.1. As shown by Yeoh [1993] the dependency on the second invariant is generally much smaller than on the first, such that often only the latter is taken into account, the reduced polynomial model: N X i ¯ ΨRP = ci I¯1 − 3 + U [J]. (3.27) i=1

Some of the most often used models for isotropic materials are obtained by a limitation of Equation (3.26) to only the first terms. The Mooney-Rivlin model restricts N to one, such that ¯ MR = c10 I¯1 − 3 + c01 I¯2 − 3 + U [J]. Ψ (3.28) An even simpler model is the Neo-Hookean where only the first term of the Reduced polynomial model is considered: ¯ NH = c1 I¯1 − 3 + U [J]. Ψ

(3.29)

Remark Equations (3.28) and (3.29) are generalizations of the Mooney-Rivlin and Neo-Hookean materials, respectively, for the compressible response. The original forms of these materials were for the incompressible response. 25

CHAPTER 3. NONLINEAR CONTINUUM MECHANICS

The Ogden model Ogden [1972] proposed a different model that does not use the invariants I1 and I2 but the principal stretches λi as arguments: ΨO =

N X 2µi i=1

αi

¯ αi + λ ¯ αi + λ ¯ αi − 3 + U [J]. λ 1 2 3

(3.30)

The parameters µi [MPa] have the dimension of a pressure, the αi [-], i = 1, . .P . , N , are dimensionless. The small strain shear modulus of this model is given by µ = N i=1 µi . Note that for the special case N = 1 [-] and α1 = 2 [-] the Neo-Hookean model (cf. Eq. (3.29)) is regained. Table 3.1 shows the expression for the first Piola-Kirchhoff stress of the Ogden model PO .

3.3.2

Green-elasticity for fiber-reinforced materials

Anisotropic materials have, in contrast to isotropic ones, direction dependent mechanical properties. In biological tissues an often encountered reason for anisotropy are thin and stiff (collagen) fibers as there are in ligaments and tendons, arterial walls or the anulus fibrosus. Ligaments and tendons behave transversally isotropic due to the presence of one parallel oriented fiber family (cf. Fig. 3.3a) while tissue from arterial walls is characterized by two crossed families of fibers (Fig. 3.3b). Anisotropy due to fibers can be included into the constitutive formulations by use of structural tensors A(i) = ζ (i) ⊗ ζ (i) (3.31)

where ζ (i) is a normalized vector describing the orientation of the i’th fiber family in the reference configuration. To the invariants defined in Eqs. (3.20) a new one can be added that depends on C and A(i) : (i)

I4 = ζ (i) · Cζ (i) = C : A(i) = (ζ (i) )2 .

(3.32)

ζ is exactly the fiber stretch. Note that the set of invariants Ii , i = 1, . . . , 4 is not the full set of invariants of the tensors C and A(i) [Antman, 2005, Chap. 12]. ¯ and a volumetric Splitting the right Cauchy-Green tensor C into an isochoric part C (i) ¯ part Cvol (cf. Eq. (3.7b)) the modified invariant I4 is obtained: (i) ¯ = J −2/3 (ζ (i) )2 = (ζ¯(i) )2 . I¯4 = A(i) : C

It should be noted that, in general, ζ¯ does not represents the fiber stretch.

(3.33)

3.3. CONSTITUTIVE RELATIONS

e3 e2 (b)

e2 (a)

Figure 3.3: Schematic of two often observed types of anisotropy in biological tissues: transverse isotropy due to parallel fibers (a) and anisotropy due to two crossed fiber families in the 1-2plane (b). One member of each fiber family is highlighted in black.

¯ is to divided it into Ψ ¯ gs and Ψ ¯ f for In case of fiber-reinforced materials an ansatz for Ψ the energy contributions of the matrix (ground substance) and the fibers, respectively: (i)

¯ gs [I¯1 , I¯2 ] + Ψ ¯f [I¯4 ]. Ψ[C, Ai ] = U [J] + Ψ

(3.34)

A fourth term adds to the first Piola-Kirchhoff stress from the isotropic case (cf. Eq. (3.24)): ¯ gs,1 F − 1 I1 F−T P = JpF−T + 2J −2/3 Ψ 3 ¯ gs,2 I1 F − 2 I2 F−T − FC + 2J −4/3 Ψ 3 ¯ f,4 FA − 1 I4 F−T . + 2J −2/3 Ψ (3.35) 3 The upcoming sections present two Green-elastic models that are used in this work by means of their strain energy functions. Their respective first Piola-Kirchhoff stresses are given in Table 3.1. A model proposed by Weiss et al. [1996] Weiss et al. suggest a constitutive model that uses the additive split of the strain energy function as presented in Equation (3.34). The individual contributions are given by κ (ln[J])2 2 ¯ gs [I¯1 , I¯2 ] = c1 I¯1 − 3 + c2 I¯2 − 3 Ψ N X c4 ¯(i) −1 c c I 3 (i) (i) 4 4 ¯f [I¯4 ] = Ψ exp − I¯4 c4 i=1 U [J] =

(3.36a) (3.36b) (3.36c)

where κ [MPa] is the small strain bulk modulus, c1 [MPa] and c2 [MPa] are parameters of the isotropic matrix and c3 [MPa] and c4 [-] are parameters of the fiber contribution. 27

CHAPTER 3. NONLINEAR CONTINUUM MECHANICS

c4 was introduced in the present work as an additional parameter that allows to shift the onset of the fiber response (c4 = 1 [-] retrieves the original model). N [-] is the number of fiber families that are modeled. Note that the contribution of the ground substance is described by a Mooney-Rivlin model with the small strain shear modulus µ = 2 (c1 + c2 ). A model proposed by Holzapfel et al. [2000] Holzapfel et al. propose a model according to Equation (3.34) for modeling the arterial wall. The explicit forms of the energy contributions are given by κ (J − 1)2 2 µ ¯ ¯ ¯ Ψgs [I1 ] = I1 − 3 2 N X 2 ¯(i) ¯f [I¯4(i) ] = k1 expk2 (I4 −1) −1 Ψ 2k2 i=1 U [J] =

(3.37a) (3.37b) (3.37c)

where κ [MPa], µ [MPa], k1 [MPa] and k2 [-] are the small strain bulk and shear moduli and two more material parameters, respectively. The ground substance is modeled as a Neo-Hookean material. Remark 1 In both previous models the fibers are thought to be active in tension only and therefore all contributions (energy, stress and stiffness) of fiber family i are set to (i) zero if I¯4 < 1 [-]. Remark 2 The model proposed by Weiss et al. can be shown to be polyconvex if c2 = 0 [MPa] and the volume ratio J does not exeed the natural number e [-]. The model proposed by Holzapfel et al. is polyconvex for all parameters and deformations. More details of the polyconvexity of these models are given in Appendix A.

3.3.3

An elasto-viscoplastic model

Biological tissues often exhibit a time dependent behavior that cannot be described by Green elastic constitutive equations. A model presented by Rubin and Bodner [2002] uses a strain energy function that includes an energy dissipating part such that viscoplastic effects can be modeled: µ0 qΨˆ exp −1 (3.38) ΨRB = 2q ˆ is the sum of four contributions where µ0 [MPa] and q [-] are material constants. Ψ ˆ =Ψ ˆ 1 [J] + Ψ ˆ 2 [I¯1 ] + Ψ ˆ 3 [ζi ] + Ψ ˆ 4 [l]. Ψ

(3.39)

ˆ i , i = 1, . . . , 4, characterize the response to total dilatation, total distorThe functions Ψ tion, fiber stretching and distortional deformation of a dissipative component, respec28

3.3. CONSTITUTIVE RELATIONS

ˆ determine the elastic response of the model: tively. The first three Ψ ˆ 1 [J] = 2m1 ((J − 1) − ln[J]) Ψ ˆ 2 [I¯1 ] = m2 I¯1 − 3 Ψ ˆ 3 [ζi ] = Ψ

(3.40a) (3.40b)

N X 2m4 m3 (i) ζ −1 . m4 i=1

(3.40c)

The mi [-], i = 1, . . . , 4, are material parameters and N [-] is the number of fiber families. The small strain bulk and shear moduli are given by µ = µ0 m1 and κ = µ0 m2 , respectively. h•i are the Macaulay brackets defined by hxi = 12 (x + |x|) and eliminate the energy contribution of the fiber family i if ζ (i) < 1 [-]. The following description of the inelastic part follows the description given in Papes and Mazza [2008]. Therein, and in the original paper by Rubin and Bodner [2002], more concise descriptions of the model can be found. ˆ 4 is a function of the distortional deformation of a dissipative component: Ψ ˆ 4 [l] = tr[b¯de ] − 3 Ψ

(3.41)

where b¯de is the state of isotropic distortion. To describe its time dependent behavior an evolution equation is set up b¯˙ de = lb¯de + b¯de lT − 32 tr[d]b¯de −Γad , | {z }

b¯˙ 0de = 1.

(3.42a,b)

¯˙ de b

˙ −1 is the spatial velocity gradient, its symmetric part is the rate of deforHere l = FF mation d = 12 l + lT . Γ is a soft switch that enables or disables the relaxation of b¯de in direction of ad such that b¯de evolves toward 1. The initial value of b¯de is (for a virgin material) given by Equation (3.42b). Γ and ad are defined by 2n

Γ = (Γ1 + Γ2 ε) ˙ exp−1/2(β/βde ) ,

ad = b¯de −

3 1. ¯ tr[b−1 de ]

(3.43a,b)

The material parameter n [-] defines the sharpness of the elastoplastic transition function Γ. Γ1 [Hz] and Γ2 [-] are additional material parameters that define the time dependent and time independent plastification, respectively. ε is the effective total distortional deformation and defined as q ε˙ = 23 d0 : d0 , d0 = d − 13 tr[d]1. (3.44a,b) In the last equation the deviator d0 of the rate of deformation tensor d is introduced. For the strain measure βde a formulation equivalent to the von-Mises stress is chosen: q b¯0de = b¯de − 13 tr[b¯de ]1. (3.45a,b) βde = 23 b¯0de : b¯0de , 29

CHAPTER 3. NONLINEAR CONTINUUM MECHANICS

The variable yield strain β introduces hardening and recovery effects: r1 r3 + r2 ε˙ β˙ = Γβde − r4 β r5 , r3 + ε˙

β0 = 0

(3.46a,b)

where the ri , i = 1, . . . , 5, are parameters (r3 and r4 have the dimension of Hertz, all others are dimensionless). While the first term raises the point of elasto-plastic transition the second lowers it such that the material returns to its initial state. For a virgin material β initially equals zero (3.46b). Table 3.1 gives the formula for the first Piola-Kirchhoff stress PRB of the elastic part. Remark 1 Note that if only the elastic part is considered (cf. Eqs. (3.38) – (3.40)) the model is overparameterized and can be reformulated with only five parameters: µ0 /q → p,

µ0 m1 → κ,

µ0 m2 → µ,

µ0 m3 → m ˜ 3,

m4 → m4 .

(3.47a–e)

Here κ [MPa] and µ [MPa] are the small strain bulk and shear moduli, respectively, and p [-] is an alternative parameter. With this new parameters the (elastic) strain energy functions reads 1 ˆ p µ0 Ψ p exp −1 (3.48) ΨRB = 2 ˆ is defined by Equation (3.39) (with Ψ ˆ 4 = 0). where Ψ In the limit case p → ∞ the elastic part simplifies to a model where the ground substance is modeled as a Neo-Hookean material: 1 ˆ lim ΨRB = µ0 Ψ. p→∞ 2 Remark 2 The elastic part of the above model is polyconvex (cf. App. A).

(3.49)

a=1

i+j=1 i≥1

N X

i=1

1 ¯ αi ¯ αi ¯ αi −1/3 ¯ αi −1 J λa − λ + λ2 + λ3 3λa 1

j≥1

N X i j−1 ∂ I¯2 j ∂ I¯1 + jcij I¯1 − 3 I¯2 − 3 I¯2 − 3 ∂F ∂F i+j=1

∂ΨO X 2µi = ∂λa αi

i−1

with

icij I¯1 − 3

3 X 1 ∂ΨO na ⊗ Na PO = λa ∂λa

PP = JpF−T +

Model by Holzapfel et al. [2000]

−T

i=1

X (i) µ ∂ I¯1 k2 + k1 I¯4 − 1 exp 2 ∂F

i=1

2 (i) I¯4 −1

(i)

∂ I¯4 ∂F

¯(i) N (i) X ∂ I¯1 ∂ I¯2 c4 I¯4 −1 (i) c4 −1 ∂ I4 ¯ + c1 + c2 + c3 exp − I4 ∂F ∂F ∂F

PH = JpF−T +

PW = JpF

Model by Weiss et al. [1996]

(3.53)

(3.52)

(3.51)

(3.50)

continued on the next page

where λa are the principal stretches and na and Na are the respective eigenvectors of the left and right Cauchy-Green deformation tensors b and C, respectively (cf. Eqs. (3.6a,b)).

Ogden model

Polynomial model

Table 3.1: 1st Piola-Kirchhoff stresses for some selected constitutive models

3.3. CONSTITUTIVE RELATIONS

32 ∂I1 ∂F ∂I2 ∂F ∂J ∂F ∂I4 ∂F ∂ζ ∂F = ζ −1 FA

= 2FA

= JF−T

= 2 (I1 F − FC)

= 2F

∂ I¯4 = 2J −2/3 FA − 31 I4 F−T ∂F ∂ ζ¯ = J −1/3 ζ −1 FA − 31 ζF−T ∂F

∂ I¯1 = 2J −2/3 F − 31 I1 F−T ∂F ∂ I¯2 = 2J −4/3 I1 F − 23 I2 F−T − FC ∂F

Derivatives of the invariants with respect to the deformation gradient F

i=1

Model by Rubin and Bodner [2002] Only the stress of the elastic part (cf. Eqs. (3.38) – (3.40)) is given using the parameters defined in Equation (3.47). ! n D E2m4 −1 ∂ζ (i) X ¯ 1 ˆ 1 1 ∂J ∂ I 1 µ Ψ PRB = exp p 0 2κ 1 − +µ +2 m ˜ 3 ζ (i) − 1 2 J ∂F ∂F ∂F

Table 3.1: 1st Piola-Kirchhoff stresses for some selected constitutive models, continued

(3.54)

CHAPTER 3. NONLINEAR CONTINUUM MECHANICS

3.4. ELASTICITY TENSORS

3.4

Elasticity tensors

The solution of finite (in-) elasticity problems requires often the use of numerical methods, e.g. the finite element (FE) method. It achieves solutions that satisfy the weak form of static or dynamic equilibrium by solving iteratively the nonlinear problem. Trial solutions are checked if they fulfill the equilibrium condition and if not a Newton-Raphson algorithm is used to find a correction. This strategy relies on the treatment of the linearized constitutive equations, which involves the derivative of the stress with respect to the deformation measures, the so called elasticity (or stiffness) tensors. The material elasticity tensor C, which is a fourth-order tensor, is defined by C=2

∂ 2 Ψ[C] ∂S[C] =4 . ∂C ∂C ⊗ ∂C

(3.55)

The last equality in (3.55) holds only true for the case of a hyperelastic material. Numerical methods often use the spatial elasticity tensor applying a Piola transformation on each slot of C: C ijkl

= J −1 FiI FjJ FkK FlL CIJKL .

which is obtained by (3.56)

A user implementation of constitutive equations for use with a FE program needs good knowledge on the required quantities. For the derivation of stiffness quantities for a certain material the reader is pointed to textbooks (e.g. Holzapfel [2000] and others). Also for experienced users it is advantageous to test an implementation with a benchmark program. MATHEMATICA (cf. App. D) allows in a rather simple way to implement a material law by just programming the strain energy density while all derivations are done numerically; a manual and hints are given in Appendix B.

Chapter

Investigation of a Non-Physical Response in Transversely Isotropic Constitutive Models Anulus fibrosus tissue, as many other soft biological tissues, is composed of thin and stiff collagen fibers in a soft matrix leading to a strong anisotropy. Commonly, constitutive models for quasi-incompressible materials, as for soft biological tissues, make use of an additive split of the Helmholz free-energy into a volumetric and an isochoric part that is applied to the matrix and fiber contribution. This split offers conceptual and numerical advantages. Working with a user implementation of a model proposed in Holzapfel et al. [2000] a non physical effect was observed. In fact, simulations involving uniaxial stress configurations reveal volume growth at rather small stretches. Numerical methods such as the Augmented Lagrangian method might be used to suppress this behavior. Investigations on this effect led to an alternative approach which solves this problem on the constitutive level.

Soft matter often exhibits complex mechanical behavior related to inhom*ogeneity, anisotropy and (near) incompressibility. This is the case for soft biological tissues that are sometimes characterized by strong anisotropy due to the presence of thin and stiff collagen fibers in a soft extracellular matrix. Examples for this are ligaments and tendons, arterial walls or the anulus fibrosus. Under physiological conditions these tissues might undergo large deformations, necessitating a mechanical description in the framework of nonlinear continuum mechanics. Numerical methods are often used to simulate and understand their mechanical behavior. This has motivated a number of recent studies on constitutive modeling approaches for fiber reinforced hyperelastic materials. Classical hyperelastic models, such as the Polynomial forms [Rivlin and Saunders, 1951] or the Ogden model [Ogden, 1972] were designed for the modeling of isotropic materials. Later these models were modified to allow for an additive split of the Helmholz free-energy into a volumetric and an isochoric part, to some extent motivated by the multiplicative decomposition of the deformation gradient introduced by Flory [1961]. 35

CHAPTER 4. A NON-PHYSICAL RESPONSE IN CONSTITUTIVE MODELS

It was shown by Ehlers and Eipper [1998] that, if this split is applied to materials not restricted to (nearly) incompressible behavior, it might lead to unphysical responses in uniaxial tension. Anisotropy itself can be modeled by explicitly including a fiber contribution in the strain energy formulation (cf. Sec. 3.3.2). Thereby the idea of the additive split might also be applied to the fiber contribution. A first fiber-reinforced model using this additive split for (nearly) incompressible materials was proposed by Weiss et al. [1996]. In 2000, Holzapfel et al. published a model that has since then become very popular in simulations involving anisotropic biological tissues. The fact that this latter model was already successfully applied to simulations of intervertebral discs [Eberlein et al., 2001] motivated an implementation in a usersubroutine for the finite element program ABAQUS (cf. App. D). Simple computations with this model produced results that did not match the expectations. A closer look revealed a non physical effect in the application of the additive split to fiber-reinforced (nearly) incompressible materials. In fact, numerical simulations involving uniaxial stress configurations lead to volume growth when using the models by Weiss et al. or Holzapfel et al. To prevent this behavior, numerical costly methods such as the Augmented Lagrangian method [Simo and Taylor, 1991] might be used. An alternative approach was found, it is to avoid the multiplicative decomposition of the deformation gradient in the fiber contribution. Although several authors [Holzapfel et al., 2004; Schr¨oder et al., 2005] have already proposed model formulations without this decomposition, no statements were found in the papers that this was motivated by the here reported problems. In general no criterion is proposed for or against this decomposition being applied to the fiber contributions to the free energy. Rather it seems the invariants and their modified equivalents are considered “to be equivalent” in case of incompressibility. Though it will be shown that the choice of either one has an influence on the stress and stiffness terms of the fiber contribution. The results are also compared with the predictions from a model proposed by Rubin and Bodner [2002], which avoids the volumetric split for the fiber contribution and uses an additive split in the log-energies. The findings reported in this chapter were presented at the Society of Engineering Science (SES) conference (2008) and published in Helfenstein et al. [2010].

4.1 4.1.1

Methods Constitutive models

The phenomenon under investigation was first observed in a model proposed by Holzapfel et al. [2000] but was found as well in another model proposed by Weiss et al. [1996]. Both models have in common that they use an additive split of the strain energy in an isochoric and a volumetric part. To complete the study another model was chosen that uses an additive split in the log-energies [Rubin and Bodner, 2002]. In Sections 3.3.2 36

4.1. METHODS

and 3.3.3 these three models are presented in some details. Only the elastic part of the model by Rubin and Bodner [2002] was used and the reparametrization described in Remark 1 (Sec. 3.3.3) was applied. For the purpose of the present analysis, the three models were considered as well in a modified variant where the measure of fiber stretch I¯4 is exchanged with the unmodified invariant I4 , for the first two models, and ζ with ζ¯ for the last model. Table 4.1 summarizes the different model formulations analyzed in this work. All models can be shown to be polyconvex, details of the proofs are given in Appendix A. The model by Weiss et al. is polyconvex only if c2 = 0 [MPa] and the volume ratio J does not exceed the natural number e [-].

4.1.2

Numerical experiment

To investigate the behavior of the highlighted models, a simple example is selected that can easily be implemented numerically. It consists of a cube of hom*ogeneous material that is aligned with a Cartesian coordinate system in 3-dimensional space (cf. Fig. 4.1). One fiber family aligned with the 1-direction reinforces the material. The load consists of an uniaxial stress state in the fiber direction (the lateral surfaces of the cube are traction free). The numerical experiment and all constitutive models were implemented in MATLAB (cf. App. D). Axial stretches are prescribed as boundary conditions, the lateral stretches are found by a numerical minimization of the strain energy using the MATLAB function fminsearch(). Calibration of the models Even though the predictive capabilities of the six models are not to be judged, they were calibrated for their responses to uniaxial tension in the fiber direction and in the transverse direction (where only the matrix is load bearing). This is done in order to activate the models similarly in the numerical experiment described above. Table 4.1: Models analyzed in this work and their references. reference

name

measure of fiber stretch

Weiss et al. [1996]

¯ W W

I¯4 I4

Holzapfel et al. [2000]

¯ H H

I¯4 I4

Rubin and Bodner [2002] (elastic part only)

¯ R R

ζ¯4 ζ4 37

CHAPTER 4. A NON-PHYSICAL RESPONSE IN CONSTITUTIVE MODELS

Figure 4.1: Example configuration consisting of a cube that is reinforced with fibers in the 1-direction and aligned with a Cartesian coordinate system. A stretch is applied in the fiber direction, whereas the lateral sides are traction free (uniaxial stress state).

¯ are selected according to the results Numerical values for the parameters of model H of an experimental program, Eberlein et al. [2001]. All other models were calibrated to provide the same uniaxial response in the axial and transverse directions with respect to the fibers for stretches up to λ1 = 1.2 [-]. The parameter set from Eberlein et al. [2001] is based on experimental data obtained for measurements with a stretch range from λ1 = 1 to 1.08 [-]. Note that we extrapolated here the validity of these parameters for stretches up to 1.2 [-] but no experimental evidence exists of the practical relevance of such a model for 1.08 < λ1 < 1.2 [-]. For the calibration the small strain shear modulus is fixed for all models such that it equals the small strain shear modulus µ∗ of the reference material. Furthermore, full incompressibility is assumed as it was done in the original parameter calibration [Eberlein et al., 2001] and therefore the parameters describing the volumetric response can not be obtained by such a procedure. A reference small strain bulk modulus κ∗ equal 2200 [MPa] is chosen, which corresponds to the bulk modulus of water [Rubin and Bodner, 2002] and is thus reasonably accurate for the materials of interest. The parameters describing the volumetric response of all models are chosen to give the same ratio of small strain bulk to shear modulus. ¯ (if J does not In order to guarantee the polyconvexity of the models W and W exceed the natural number e [-]) the parameter c2 is set equal 0.0 [MPa] a priori. Doing ¯ , H and H ¯ have the same ground substance contribution to the so, the models W , W Helmholz free energy and their response to uniaxial tension in the direction transverse to the fibers is the same. As remarked at the end of Section 3.3.3 the elastic part of the model proposed by Rubin and Bodner reduces in the limit case p → ∞ [MPa] (cf. Eq. (3.47a)) to a NeoHookean material similar to the models proposed by Weiss et al. (with c2 = 0 [MPa]) and Holzapfel et al. In order to preserve the special structure of this model, p is set equal 1 [MPa]. Unfortunately this goes in hand with a non perfect calibration of the ¯ in the direction transverse to the fibers. models R and R Table 4.2 reports the parameters for all models that are used. 38

4.1. METHODS Table 4.2: Material parameters of the models that were analyzed. The parameter set of model ¯ was used as reference for the calibration of the other models. H

¯ W W

¯ H H

¯ R R

κ [MPa]

c1 [MPa]

c2 [MPa]

c3 [MPa]

c4 [-]

κ∗ κ∗

1 ∗ 2µ 1 ∗ 2µ

0.000 0.000

3.362 × 10−4 6.689 × 10−4

3.864 × 101 3.865 × 101

κ [MPa]

µ [MPa]

k1 [MPa]

k2 [-]

κ∗ κ∗

µ∗ µ∗

3.000 6.000

4.500 × 101 4.500 × 101

p [MPa]

κ [MPa]

µ [MPa]

m ˜ 3 [MPa]

m4 [-]

1.000 1.000

κ∗ κ∗

µ∗ µ∗

1.539 × 103 1.539 × 103

1.583 1.583

− c2 − c2

The reference bulk and shear moduli are given by κ∗ = 2200 [MPa] and µ∗ = 0.5 [MPa].

4.1.3

Tangent Poisson’s ratio

The kinematic response of a cube under uniaxial tension is (normally) characterized by contractions in the lateral directions.1 This should be true as well for a fiber reinforced cube under uniaxial tension in the fiber direction. In order to get information about the lateral behavior of the example configuration (Sec. 4.1.2) the tangent Poisson’s ratio νtan is calculated. To begin with, the time rate of the Cauchy stress σ of a hyperelastic material can be written as ˙ ˙ T, σ˙ = J −1 FSFT = −tr[l]σ + lσ + σlT + J −1 FSF where S˙ is the time rate of the second Piola-Kirchhoff stress: S˙ = 12 C : C˙ = 12 C : 2FT lF .

(4.1)

(4.2)

Here C is the material stiffness (cf. Eq. 3.55). The spatial velocity gradient l is defined ˙ −T . by l = FF Insertion of Equation (4.2) into (4.1) yields the expression: σ˙ = −tr[l]σ + lσ + σlT + C : l,

(4.3)

where C is the spatial material stiffness (cf. Eq. 3.56). Now the example configuration described in Section 4.1.2 is considered. Whether the material is fiber-reinforced or not is not important for the following considerations. Let the cube be under an uniaxial load in 1-direction, i.e. it is elongated in 1-direction with a given stretch λ1 ≥ 1 [-]. Since only the axial stress is nonzero and the transverse 1

Except for some rare materials [e.g. Evans, 1991; Lubarda and Meyers, 1999; Milton, 1992] the Poisson’s ratio is positive (0 ≤ ν ≤ 0.5). No evidence was found in literature that there exist natural biological tissues that exhibit a negative Poisson’s ratio.

CHAPTER 4. A NON-PHYSICAL RESPONSE IN CONSTITUTIVE MODELS

stresses equal zero for all times, it holds σ˙ 11 = −tr[l]σ11 + 2l11 σ11 + [C : l]11 ≥ 0 σ˙ 22 = [C : l]22 ≡ 0 σ˙ 33 = [C : l]33 ≡ 0

(4.4a) (4.4b) (4.4c)

For the given uniaxial stress state the deformation and spatial velocity gradients, F and l, respectively, are given by   λ1 0 0 F =  0 λ2 0  (4.5) 0 0 λ3 

˙ ln[λ1 ]

 ˙ −1 =  l = FF 

0 ˙ ln[λ2 ]



 0  . ˙ ln[λ3 ]

(4.6)

With the above relation for the spatial velocity gradient l it is possible to write the system of Equations (4.4) in matrix form. If an isotropic system or a transversely isotropic system is considered with its principal direction aligned with the 1-direction, the 2- and 3-directions are equivalent.        ˙ ˙ σ˙ σ + C 1111 −2σ11 + 2C 1122 ln[λ1 ] ln[λ1 ]   11  =  11  = k  (4.7) ˙ ˙ 0 C 1122 C 2222 + C 2233 ln[λ2 ] ln[λ2 ] At this step it is important to remind that the two stiffness components C 2222 and C 2233 are not the same even though the directions 2 and 3 are equivalent. Inversion of the above Equation (4.7) leads to the expressions for the rates of the logarithmic stretches: ˙ ln[λ1 ] =

+ C 2233 σ˙ 11 det[k]

C 2222

C 1122 ˙ ln[λ2 ] = − σ˙ 11 . det[k]

(4.8a) (4.8b)

According to Bahuaud and Boivin [1968] a tangent Poisson’s ratio is defined by νtan

˙ ln[λ2 ] λ1 λ˙ 2 := − =− . ˙ λ2 λ˙ 1 ln[λ ]

(4.9)

In contrast to the definition of a secant Poisson’s ratio νsec 2 the above Equation (4.9) is sensitive to the rate of change of the volume. 2

νsec = − (λ2 − 1) / (λ1 − 1)

4.2. RESULTS

Using Equations (4.8) together with (4.9) the tangent Poisson’s ratio can be written in terms of the spatial stiffness only: νtan =

C 1122 C 2222

+ C 2233

(4.10)

By replacing the stiffness terms by their corresponding linear elastic components it is straightforward to see that for linear elasticity νtan corresponds with the classical Poisson’s ratio ν. Remark An increasing stress (σ˙ > 0 [MPa/s]) will result in an increasing axial stretch (λ˙ 1 > 0 [MPa]). According to Equation (4.8a) and the monotonicity of the logarithm, this is only the case if C 2222 + C 2233 > 0 [MPa] (given that k is positive definite, what still needs to be proven). Therefore the only way for νtan to become negative is a negative stiffness component C 1122 (cf. Eq. 4.10). Therefore if C 1122 < 0 [MPa] the rate of the transverse stretch becomes positive, and therefore the lateral stretch begins to increase.

4.1.4

Condition number

A condition number c can be defined [cf. e.g. Stoffer and Nipp, 1998] as the ratio of the maximum and the minimum eigenvalue of the spatial stiffness C: c :=

max[eig[C]] . min[eig[C]]

(4.11)

A low condition number is favorable in computations.

4.2

Results

Simulations of the numerical example (cf. Sec. 4.1.2), using the quasi-incompressible formulations of the constitutive equations, reveal that only the models W , H and R reproduce the high stresses they showed while using the assumption of full incompressibility ¯ shows a slightly softer behavior while the two remaining (cf. Fig. 4.2). The model R ¯ and H ¯ are significantly too soft. models W This goes in hand with the observation that the three models using the modified ¯ respectively, already violate the incompressibility measures of the fiber stretch I¯4 and ζ, constraint at rather small stretches (cf. Fig. 4.3). An alternative and instructive representation of this effect is given in Figure 4.4. ¯ and H ¯ the tangent Poisson’s ratio νtan (cf. Eq. 4.10) decreases, after For the models W being constant in the beginning, and even becomes negative at a critical stretch λcrit , ¯ does such that with increasing axial stretch the lateral stretch grows as well. Model R not keep the Poisson’s ratio constant equal 0.5 [-] either, but shows a smaller deviation ¯ and H. ¯ The models W , H and R keep a tangent Poisson’s ratio (close than models W to) equal 0.5 [-], which agrees well with an (quasi) incompressible material. ¯ Figure 4.5 shows the condition number (cf. Eq. 4.11) of all models. The models W ¯ have an only slightly increasing condition number until a critical stretch where and H 41

CHAPTER 4. A NON-PHYSICAL RESPONSE IN CONSTITUTIVE MODELS

25 ¯ W

¯ H

¯ R

σ11 [GPa]

15 R H W

¯ R

¯ W 1

1.02

1.04

1.06

1.08

1.1 λ1 [-]

1.12

1.14

1.16

¯ H

1.18

1.2

Figure 4.2: Cauchy-stress σ11 in axial direction plotted against the axial stretch λ1 . The ¯ and H ¯ do not reproduce the high stresses observed while using the assumption of models W ¯ behaves slightly softer. incompressibility during the parameter calibration; the model R

1.12 ¯ W

¯ H

¯ R

1.1

J [-]

1.08 ¯ H

1.06 1.04

¯ W

1.02

¯ R

1 W H R 0.98

1.02

1.04

1.06

1.08

1.1 λ1 [-]

1.12

1.14

1.16

1.18

1.2

Figure 4.3: Volume ratio J plotted against the axial stretch λ1 . Models W , H and R preserve ¯ and H. ¯ the volume. The largest deviations from J = 1 [-] are obtained by the models W

4.2. RESULTS

W H R 0.5 ¯ R

νtan [-]

0.25 0 −0.25

¯ H λcrit

−0.5 −0.75

¯ W

¯ H

¯ R

1.02

1.04

1.06

1.08

1.1 λ1 [-]

1.12

¯ W 1.14

1.16

1.18

1.2

Figure 4.4: Tangent Poisson’s ratio νtan (cf. Eq. (4.10)) plotted against the axial stretch λ1 . ¯ and H ¯ the ratio becomes negative (at λcrit , shown for W ¯ ) indicating growing For models W lateral stretches.

108 ¯ W

¯ H

¯ R

107 λcrit. c [-]

H 106 W

¯ W

105

¯ H 104

1.02

1.04

1.06

1.08

1.1 λ1 [-]

1.12

1.14

1.16

1.18

R ¯ R 1.2

Figure 4.5: Condition number c (cf. Eq. (4.11)) plotted against the axial stretch λ1 . Up to the ¯ ) the condition numbers of models W ¯ and H ¯ increase only critical stretch λcrit (shown for W slightly, while the condition numbers for models W and H increase rapidly. 43

CHAPTER 4. A NON-PHYSICAL RESPONSE IN CONSTITUTIVE MODELS

the condition number begins to increase rapidly (in Fig. 4.5 the critical stretch λcrit ¯ ). This critical stretch is exactly at the value for which the tangent is shown for W Poisson’s ratio νtan passes through zero (c1122 = 0 [MPa], cf. Fig. 4.4). An analysis of the eigenvectors of the spatial stiffness tensor C reveals that for λ1 > λcrit , the deformation with the largest increase in the deformation energy is no longer volume change but rather axial stretch in the fiber direction. The condition numbers for the models W and H (using I4 ) increase much more than for all other formulations and are, at λ1 = 1.2 [-], two orders of magnitude higher than initially. ¯ and R show only slightly increasing condition numbers. Both models R

4.3

Discussion

¯, The calculations show that for the quasi-incompressible formulation of the models W ¯ and R ¯ the volume ratio J does not stay constant a priori (cf. Fig. 4.3). The first H two models have in common that the isochoric contribution of the ground substance and the volumetric contribution have polynomial forms while the fiber contributions are of exponential type. The latter model has a polynomial log-energy contribution of ¯f and Ψ ˆ 3 (cf. Eqs. (3.36c), (3.37c) the fibers. This means that the fiber contributions Ψ and (3.40c)), respectively, grow much faster than the ground substance and volumetric contributions as their respective arguments increase. ¯ respectively, decrease with an increasing ¯f and Ψ ˆ 3 , i.e. I¯4 and ζ, The arguments of Ψ volumetric deformation (cf. Eq. (3.33)). Since the axial stretch is given by boundary conditions, the strain energy might be lowered through lateral expansion leading to an increasing volume ratio J and to a reduction in stress when compared to the incompressible case (cf. Fig. 4.2). Thus the material tends, under uniaxial load, to a purely volumetric deformation state in order to lower the fiber energy at the expense of an increase in the ground substance and volumetric energy. The level of deformation at which this effect will take place depends on the specific formulation and the material parameters, in particular on the selected value of bulk modulus. In Section 4.1.3 is shown that, in analogy to linear elasticity, a tangent Poisson’s ratio (Eq. (4.10)) can be defined. This Poisson’s ratio needs to be (close to) 0.5 [-] in ¯ and H ¯ the order to describe an (nearly) incompressible material. For both models W tangent Poisson’s ratio deviates from the incompressible limit and becomes negative (cf. Fig. 4.4) even at seemingly innocuous strain levels. While the tangent Poisson’s ¯ departs as well from 0.5 [-] it does not do so as much as the other two models ratio of R using the modified invariant I¯4 . The model proposed by Rubin and Bodner seems to be less sensitive to a use of ¯ The reason for this might be the fact that all contributions the modified invariant ζ. to the strain energy are in the exponent and are therefore more equilibrated in their contribution to the overall strain energy. Changing the argument of the fiber contribution to I4 disables the ability of the material to reduce its total energy by a more spherical deformation, due to the changed energetic couplings. Simulations indicate that with this change the volume ratio stays 44

4.3. DISCUSSION

constant, equal to unity, and the tangent Poisson’s ratio is kept close to 0.5 [-]. The unphysical response of growing lateral stretches in uniaxial tension has previously been reported in a different context by Ehlers and Eipper [1998]. Their setting differs in two major points from the here used one: 1) they study isotropic materials and 2) in contrast to the here made assumption of (near) incompressible materials, they consider compressible materials with κ/µ = 5/3 [-]. Ehlers and Eipper report growing lateral stretches when using strain-energy forms that use an additive split into deviatoric and volumetric part equivalent to Equation (3.34) without a fiber contribution. For axial stretches of up to λ1 = 10 [-] none of the three models W , H and R of this study shows a growing lateral stretch when used without fibers (c3 = 0 [MPa], k1 = 0 [MPa] and m1 = 0 [MPa], respectively). This was to be expected since here a ratio κ/µ as high as 4400 [MPa] is used. The here reported unphysical response originates from the fiberreinforcement and is distinct from that studied by Ehlers and Eipper. The underlying mechanism leading to this behavior, however, is the same: the total energy of the system can be lowered with a more spherical deformation in case of a large tangent stiffness related to the isochoric energy contribution. The models using I¯4 and ζ¯ can reduce the fiber contribution to the stiffness by volume growth and thus keep the condition number “moderate” (cf. Fig. 4.5). Using I4 instead, the fibers contribute according to their true stretch and generate large entries in the stiffness matrix, thus leading to a dramatic enhancement of the condition number. This is not the case for the model proposed by Rubin and Bodner, where an additive split is applied to the log-energies (cf. Eq. 3.38); this formulation keeps the individual contributions closer together. The price for this is a more involved physical interpretation of each single term in the energy function as compared with the formulations using the additive split of the single contributions of volumetric and isochoric deformation of the isotropic ground substance and fiber stretch. It has to be mentioned that if the true nature of a material is such that it combines fibers of exponentially increasing stiffness with a soft matrix the increase of the condition number corresponds to the physical reality and is therefore not bad or wrong a priori. Its negative implications for numerics have to be solved on the computational side. ¯f also has consequences for the contributions of the The choice of the argument of Ψ fibers to the stress and stiffness. The Cauchy stresses due to the fibers can be written as ¯f ∼ F σ

1 ∂ I¯4 T F = J −2/3 a − I¯4 1, ∂C 3

σf ∼ F

∂I4 T F = a. ∂C

(4.12a,b)

Here a is the structural tensor in the actual, deformed, configuration. It can clearly be seen that the use of I4 leads only to stresses along the fibers while using I¯4 induces stresses with components perpendicular to the fiber direction. The latter does not correspond to the idea of fibers behaving like uniaxial springs where 1) the resistance only depends on the stretch and 2) the corresponding force is in the direction of the spring. Also in the stiffness terms the use of I¯4 induces fiber contributions that are orthogonal to the fiber direction. A more detailed discussion about the implications of either using I¯4 or I4 on the 45

CHAPTER 4. A NON-PHYSICAL RESPONSE IN CONSTITUTIVE MODELS

coupling between isochoric and volumetric parts of deformation and stress can be found in Sansour [2008]. The reported effect of undesired volume growth can be avoided by introducing specific constraints in the numerical scheme such as Augmented Lagrangians, which effectively entails a change in the actual model to a strictly incompressible one. The alternative approach, here suggested, which preserves the near incompressible character of the material, is to use I4 in place of I¯4 for the fiber contribution to the strain energy function. Weakening the principle of the additive split and allowing the fibers to contribute according to the total deformation solves the problem on a constitutive level. Remark From version 6.8 on, the material proposed by Holzapfel et al. is implemented in ABAQUS (in a version that also accounts for fiber dispersion [Gasser et al., 2006]). Considering the results obtained in this study it is interesting to see how the implementation in ABAQUS behaves when simulating the numerical example (cf. Sec. 4.1.2). Results are sobering. Using a fully incompressible material (1/κ = 0 [1/MPa]) no transverse growth is found, simulations crash when the Cauchy stress in axial direction is in the order of 1019 MPa. If a nearly incompressible material is used (κ = κ∗ ) things become worse. Lateral growth appears at rather the same critical stretch λcrit as found in this study. Differences in the critical stretch might arise due to a different formulation of the isochoric strain energy function. But still the non-physical effect is reproducible.

Chapter

Optimal Specimen Design for Planar-Biaxial Materials Testing of Soft Materials Planar biaxial tests have become popular to complement classic uniaxial tests on soft materials. For their ease of handling cruciform shaped specimen geometries are often used. So far little research was done in the soft tissue community to improve the specimen geometry with respect to materials parameter determination. Following the idea of M¨ onch and Galster [1963] a numerical method to maximize the area of hom*ogeneous equibiaxial loading is presented. The findings of the study suggest that already a low number of slots in the sample limbs help to increase this area significantly.

In order to characterize soft materials like rubbers or biological tissues, planar-biaxial tests have become popular to complement classic uniaxial materials tests. It is known that the latter only are not sufficient to fully characterize the materials response [Bass et al., 2004]. A first planar-biaxial test setup for soft biological tissues was presented by Lanir and Fung [1974]. Ever since the interest in biaxial materials testing was growing in the soft tissue community [Sacks, 2000]. An especially appealing feature of planar-biaxial testing is that a typical setup (i.e. independent control of the two principal test axes) allows for an elastic, isotropic and (nearly) incompressible material to assess its complete specific strain energy potential, i.e. to test over the complete strain-space. Starting in the field of metals testing, planar biaxial materials tests often use cruciform shaped specimen geometries for their ease of handling (e.g. clamping). Unlike in other scientific communitiesi , so far little research was done in the soft tissue community to improve the specimen geometry with respect to the model parameter determination [Waldman and Lee, 2005]. Results from the field of metals testing cannot be applied directly to soft tissue testing as here much larger deformations are considered. Additionally to the challenge of large displacements that come with soft biological i

Abdul-Aziz and Krause [2006]; Demmerle and Boehler [1993]; Hardacker [1981]; Smits et al. [2006]; and others

CHAPTER 5. SPECIMEN DESIGN FOR PLANAR-BIAXIAL MATERIALS TESTING

tissues, these often present an anisotropic behavior due to the presence of fibers. Waldman et al. [2002] studied the effect of different holding methods at square samples of biological tissues. They found that suturing sample edges weakens the load response compared to clamping the sample edges and relate this to the localized force application in the first case [as well Waldman and Lee, 2002]. Waldman and Lee [2005] studied the influence of the sample arm length for cruciform shape specimens and conclude that this geometry is non appropriate for connective tissue testing. What was not done so far is to optimize the cruciform sample geometry with respect to the size of the biaxial area. This study focusses on isotropic soft elastic materials that undergo large deformations. To maximize the biaxially loaded area it helps to slot the limbs of the cruciform specimen [M¨onch and Galster, 1963]. In this light, the specimen geometry along with the force transmission into the specimen was optimized towards uniformity of the induced strain- and stress-field in the test region by use of a numerical optimization scheme. The specimen shapes obtained by the in-silico study were realized as silicone samples. With these samples the predictions were verified by means of the displacement fields obtained under equibiaxial stretch. Most of the results presented in this chapter were presented at the 6th European Conference on Constitutive Models for Rubber (ECCMR) Helfenstein et al. [2009a].

5.1

Motivation

In classic materials testing the forces and displacements at the clamp interfaces to the specimen are measured. Based on these global measurements the local stress and strain state in the specimens test region must be assessed in order to calibrate constitutive models by means of regression analyses. This in turn implies that the states of strain and stress in the test region must be compatible with the nominal stress and strain values obtained from the global measurements, i.e. the clamp force divided by the corresponding cross-sectional area and the displacement divided through the clamp to clamp initial free gauge length, respectively. Therefore, appropriate design of the test setup and specimen geometry is a fundamental prerequisite to materials testing in order to establish the required well-defined state of strain and stress in the test region. In contrast, one could argue that once the interface forces, displacements and the probe geometry are known, the setup can be modeled within a finite element (FE) framework and the material parameters determined from solving the inverse problem. What follows demonstrates that this generally is not a suitable approach. To this end an in-silico planar biaxial test is performed on two non-slotted cruciform specimen geometries with different limb lengths, where the material behavior is prescribed as an Ogden-type (cf. Sec. 3.3.1). A nominal stretch λnom = 1.3 [-] is applied equally to both cruciform axes (equibiaxial tension) and the resulting reaction forces are computed at the virtual clamp interfaces. The prescribed clamp displacements and the computed reaction forces are then used to calibrate a third order Reduced polynomial (cf. Sec. 3.3.1) to the test data by solving the inverse problem. The used material parameters and the results 48

5.2. METHODS

are given in Table 5.1 and figure 5.1, respectively. When plotting the equibiaxial characteristic of the in-silico material against the reduced polynomial fit for the short and long limbs, two observations become apparent: 1) the results obtained from the planar biaxial tests on an regular cruciform specimen do not predict the true equibiaxial material behavior, i.e. the obtained material parameters are not associated with an equibiaxial state of deformation, and 2) the results are specimen geometry dependent. The observations can be explained, as the specimen test region is in a mixed deformation mode; there exits no clear relation between the nominal stress and strain values and the prevailing stress and strain fields. Thus a similar error is to be expected if the constitutive parameters were estimated from the nominal measurements. Concluding, this demonstrates the importance of the specimen geometry and the requirement for explicit assignment of the model parameters to the mode of deformation in nonlinear materials testing. These observations motivate appropriate design of specimens and the general test setup such that the state of strain and the associated state of stress in the test region become accessible by means of the available measuring data. Typically, in planar biaxial test setups in-plane deformations are measured by means of optical methods; the full kinematics is well-defined if incompressible material behavior can be assumed. In order to obtain a significant relation between the global force measurements and the local stress state, the hom*ogeneous biaxial test region in the mid-section of the cruciform needs to be maximized.

5.2

Methods

The problem of biaxial stretch in a hom*ogeneous isotropic probe can be simplified by several means: 1) the symmetries of the problem can be used such that only one quarter of the specimen has to be considered1 , 2) far away from the clamping the loading case can 1

In fact, even only one eighth of the specimen would be sufficient. The problem herewith is that boundary conditions defined in two different coordinate systems are necessary. These BCs meet in one point, what causes problems to ABAQUS. In later studies a solution to this problem was found, it

Table 5.1: Material parameters found by calibration of a Reduced polynomial model on data generated with an Ogden type material. Ogdena

µ1 = 4.62 [MPa],

α1 = 3.95 [-]

µ2 =−2.65 [MPa],

α2 =−0.8 [-]

RP, short limbs

c1 = 1.02 [MPa],

c2 =−0.51 [MPa],

c3 =0.55 [MPa]

RP, long limbs

c1 = 1.12 [MPa],

c2 =−0.35 [MPa],

c3 =0.46 [MPa]

The parameter set was obtained by calibration of the Ogden model with experimental data on vulcanized rubber [Treloar, 1944] (the experimental data was rescaled in order to resemble data obtained from biological tissues). Later use of the parameter set revealed its invalidity, it was found that under quasi -uniaxial load the material starts to grow laterally; it does not grow laterally under pure-uniaxial load.

CHAPTER 5. SPECIMEN DESIGN FOR PLANAR-BIAXIAL MATERIALS TESTING

reference material RP short limbs RP long limbs

Cauchy stress [MPa]

2.5 2 1.5 1 0.5 0 1

1.05

1.1

1.15 λnom [-]

1.2

1.25

1.3

Figure 5.1: Comparison between the equibiaxial characteristics of the in-silico reference material and the calibrated reduced polynomials (RP).

be considered as a plane stress situation. In addition to these two general assumptions the study is restricted to the case of equi -biaxial loading. Definition of the biaxial area A point is assumed to be under equibiaxial and hom*ogeneous load if - its principal stresses do not differ for more than τ1 from each other and - the mean principal stress does not differ for more than τ2 from the mean principal stress at the center point. The side length L of the maximum square that is inscribed to this set of points, without including any other points, is called the size of the biaxial area (Fig. 5.2b). Simulations For the maximization of the biaxial areas size MATLAB (cf. App. D) optimization algorithms were used in a framework developed by Hollenstein [2011]. As the cooperation of MATLAB and the FE program ABAQUS (cf. App. D) is sometimes tricky some of the adjustments made are explained in Appendix C.1. Two algorithms from the optimization toolbox of MATLAB are used: an unconstrained nonlinear optimization (fminsearch) that uses a Nelder-Mead simplex method [Lagarias et al., 1998] and a constrained genetic algorithm (ga) [Goldberg, 1989]. The consists of disabling the consistency checking when the job is started (consistencyChecking=OFF).

5.2. METHODS y

b ∆

a li di

ui wi

∆ x

(a)

L 2

(b)

Figure 5.2: Quarter of a cruciform specimen: undeformed with parametrization of the geometry and boundary conditions (a) and deformed with area of points under equibiaxial and hom*ogeneous load (hatched vertically) and the biaxial area (hatched horizontally) (b).

genetic algorithm allows to define an initial population consisting of parameter sets and the corresponding fitness values (errors). For parameter sets that cover a wide range of the design space the size of the biaxial area L (and the respective error 1/L) were computed at the deformed state with a global stretch of 1.5 [-]. The two tolerances τ1 and τ2 were set to 5 % and 1 %, respectively. The twenty best parameter sets (corresponding to the size of the initial population) were taken as initial population. Both algorithms call a function that starts the FE software ABAQUS and that, after termination of the simulation, computes the error defined as 1/L. After being called ABAQUS starts a PYTHON script (cf. App. C.2) that generates the geometry, assigns the material properties, meshes the geometry and initiates the FE computation. In order to transfer the trial parameter set to ABAQUS, the MATLAB code writes them in a text file that is read by the PYTHON script. Vice versa, the results from ABAQUS (coordinates and principal stresses of all nodes) are written to a text file that is read by the MATLAB routine. Figure 5.2a shows the parametrization of the geometry with slots in the specimen limbs. The symmetries with respect to the principal axes of the specimen are exploited and only one quarter is simulated. The end of the limbs are clamped, displacements (∆) are applied such that global stretches of 1.5 [-] are obtained. The fixed dimensions a = 40 [mm] and b = 20 [mm] are half of the total size of the cruciform specimen and half of the limb width, respectively. The set {wi , li , ui , di }, i = 1, . . . , n, parameterizes one single slot. n is the total number of slots in a half limb. wi designates the width, li the length, ui the distance from the slot to the limbs end and di the distance from the centerline of the i − 1 slot to the centerline of the ith slot. d1 is the distance from the limbs centerline to the centerline of the first slot. Table 5.2 summarizes the constraints that are used for the parameters in the case of the genetic algorithm. Condition C7 ensures that there is a ligament of at least 1.0 mm 51

CHAPTER 5. SPECIMEN DESIGN FOR PLANAR-BIAXIAL MATERIALS TESTING

between two slots. Simulations are done with quadrilateral two dimensional eight node biquadratic plane stress elements (CPS8). The typical element size is set to 0.5 mm. Incompressibility is assumed for the material, a typical assumption often made for soft biological tissues and rubbers. The elastic properties are modeled as Neo-Hookean with an initial shear modulus µ∗ = 2 [MPa]. Experimental validation To validate the numerical results two silicone samples that correspond to the situation without any slot and to the best solution with four slots per limb were fabricated. The outlines follow the numerically obtained geometries, the thickness is 2 mm. In order to clamp the samples correctly the limbs are made 15 mm longer. Aluminum casting boxes were milled using a computerized numerical control machine. With this molds silicone samples were cast (cf. App. D). The samples are fixed with custom made clampings in our biaxial materials testing machine (cf. App. D) and tested with global stretches up to 1.5 [-]. The displacements are captured by a digital camera installed over the test area. Red chalk is used to produce a stochastic pattern on the surface such that the displacement field can be determined using a digital image correlation software (cf. App. D).

5.3

Results

First optimizations were done using the Nelder-Mead simplex method but it turns out that the solutions of this optimization procedure stay in the very same region as defined by the initial parameter set. Therefore we changed to a genetic algorithm. Table 5.3 summarizes the best solutions with no, two, four and six slots per limb. When no slots are present the size of the biaxial area is equal to 9.966 mm. Already two slots give a significant increase of 242 %, L = 34.058 [mm]. Including more and more slots the size of the biaxial area increases as well. In the solution with four cuts the innermost slots do coincide (d1 = 0 [mm]) such that in fact only three cuts are present. The best solution obtained with six slots raises L up to 38 mm. The second slot in this solution is very short compared to the others (l2 = 1.815 [mm]) but neglecting this slots degrades the length L to 23.292 mm (−39 %). The left of Figure 5.3 shows the simulated displacement fields at the central part of two configurations (only one quarter is shown). On top the displacement field where Table 5.2: Constraints that are used with the genetic algorithm.

0.2

(C1)

(C2)

0.2

(C3)

(C4)

0.2

(C5)

≥

(C6)

0.5 × (wi−1 + wi ) + 1.0 for i > 2

(C7)

5.4. DISCUSSION Table 5.3: Results of the optimizations. The first column (]) indicates the number of slots per limb. ]

L [mm]

0 2 4 6

9.966 34.058 36.066 38.000

wi [mm]

li [mm]

ui [mm]

di [mm]

2.938 1.133, 3.000 1.214, 1.402, 1.862

15.569 15.750, 18.000 18.427, 1.815, 18.013

5.306 4.000, 2.000 1.977, 18.250, 2.328

5.722 0.000, 7.000 2.488, 3.326, 4.088

no slots are used (a): it can be seen, that only a rather small area with a side length of approximately the above mentioned 5 mm can be considered as equibiaxially loaded. In contrast to the case where four (three) slots per limb are simulated (b), here the hom*ogeneity of the deformation increases significantly. These in-silico results are in good agreement with the experimental results shown on the right of Figure 5.3. Since the stiffness of the limbs is smaller than the stiffness of the central part, the global stretch of 1.5 [-] is not fully regained at the center part where local stretches in the order of 1.3 [-] are obtained. Figure 5.4 shows the results of a study of the force flux in the specimen without and with four slots per limb (the forces are normalized with the maximum value of the respective applied total force). The imposed force at the clamping is partially transferred to the biaxial area. Without any slots 14 % of the total force are transferred, using four slots increases the ratio up to 62 %. This ratios are not constant and change, in the case with four slots, for approximately 4 % over the whole deformation process (1 % in the case without any slots). The standard deviation of the force acting on different sections in the biaxial area is in the range of 0.7 % of the mean force amplitude.

5.4

Discussion

A point is considered under equibiaxial and hom*ogeneous load if its principal stresses do not differ for more than τ1 = 5 % and the mean stress at that point do not differ for more than τ2 = 1 % from the mean stress at the center point of the probe. The size of the biaxial area is measured by the side length L of the largest square inscribed to that area that contains no other points. The hom*ogeneous equibiaxial deformation of a square can be described analytically such that an analytical regression analysis can be envisaged. A cruciform specimen geometry with a total width of 80 mm and a limb length of 20 mm each is used, such that the central region of 40 × 40 mm2 results. The arms ends are clamped and displacements of 20 mm are applied to each (such that global stretches up to 1.5 result). For a non-slotted specimen the biaxial area L has a size of 9.966 mm at that stretch. A consequence of this small biaxial area is that only little of the specimen is under 53

CHAPTER 5. SPECIMEN DESIGN FOR PLANAR-BIAXIAL MATERIALS TESTING

20 [mm]

[mm]

10 15 [mm]

15 20 [mm]

(a) sample without slots

20 [mm]

[mm]

10 15 [mm]

(b) sample with four slots reference configuration

deformed configuration

Figure 5.3: Displacement fields observed in the simulations (left) and experiments (right): sample without any slots (a), sample with 4 slots per limb (b).

5.4. DISCUSSION

1.2 imposed force

mean force at center

ratio

0.8

0.6

0.8 0.6

0.4

ratio [-]

normalized force [-]

0.4 0.2 0.2 0

0 1

1.1

1.2 1.3 global stretch [-]

1.4

1.5

Figure 5.4: Evolution of the imposed force, the mean force that goes through the biaxial area and the ratio thereof with the global stretch. (◦: specimen without slots, ×: specimen with 4 slots per limb)

biaxial load and the global response is dominated by other load cases. Only 14 % of the global force goes through the sections of the biaxial area. In order to increase the size of the biaxially loaded area M¨onch and Galster [1963] proposed to cut slots in the specimen limbs. Adapting their idea to numerical optimization methods allows to search for non intuitive optimal solutions for the positioning and size of the slots. The best solution found with two slots per limb increases the size of the biaxial area by 242 %. This means that a larger portion of the specimen is effectively under biaxial load. Including up to four more slots increases the size of the biaxial area by another 11 %. The solution with six slots shows that rather small changes influence the solution significantly. Omission of the two smallest slots with a length of only 10 % of the next longer slots degrades the solution by 39 %. The left of Figure 5.3 shows the computed displacement fields for the two configurations with no and with four slots per limb. As expected from the computed sizes L the deformation field for the latter configuration is much more hom*ogeneous for a larger area. A comparison with the experimentally determined deformation fields (on the right of Fig. 5.3) shows a good (qualitative) agreement and verifies the in-silico study. In case of the geometry with four slots the ratio between the global force and the force acting on the sections of the biaxial area is significantly increased compared with the situation without any slots (62 % vs. 14 %). This finding suggests that most of the global force is held by the biaxial area and that the influence of modes other than 55

CHAPTER 5. SPECIMEN DESIGN FOR PLANAR-BIAXIAL MATERIALS TESTING

equibiaxial is reduced. With ongoing deformations the force flux varies only negligibly considering other imprecisions that occur during soft tissue testing. Unfortunately this small variation is due to the material chosen for the simulations and not an overall property.

5.5

Conclusion

Material parameters determination from planar biaxial materials testing needs knowledge of the local strains and stresses at the biaxially loaded part in the specimen. As shown in Section 5.1 fitting global displacements and forces acting on cruciform specimens only can lead to significant errors. Local strains can be measured by optical methods and are therefore directly accessible. This is different for the local stresses, their value has to be estimated from the global forces measured at the clampings. The ratio between global forces and local stresses depends on the geometry of the specimens and, for the case of nonlinear material behavior as well, on the material parameters. Linear elastic materials are a special case, where the latter dependency does not hold. For such materials Demmerle and Boehler [1993] presented a specimen shape for which the local stress equals the nominal stress. For soft biological materials that have generally a nonlinear mechanical behavior this dependency of the ratio between global and local forces can be minimized by enlarging the biaxial area. A large area under biaxial load ensures that most of the force flux passes it and therefore the dependency of the ratio becomes smaller. If, in addition, the stress field in the biaxial area is hom*ogeneous the stresses can be determined from the globally acting forces. Following the idea of M¨onch and Galster [1963] a numerical method to maximize the area of hom*ogeneous (equi-) biaxial loading in cruciform shaped specimens undergoing biaxial materials testing is presented. The findings suggest that already a low number of slots helps to increase the biaxial area significantly.

Chapter

Uniaxial Experiments with Anulus Fibrosus Samples Mostly shear and mixed deformation modes are present in the anulus fibrosus, but uniaxial and biaxial deformed zones are present as well. Therefore the need to characterize this tissue under several different loads. Eventually, feasibility reasons led to an uniaxial characterization only. Prestudies explored the influence of the aspect ratio and the sample alinement on the outcome of uniaxial experiments. The compliance of the used force transducers were quantified and the influence of freeze storage on the mechanical properties was qualified. Circumferential anulus fibrosus samples were prepared from fresh ovine lumbar spines. Experiments included displacement driven load ramps with intermediate relaxation phases at defined strain levels. Interspecimen results are qualitatively but not quantitatively very similar and show the presence of an initial toe-region, a distinct hysteresis and stiffening.

Although it is know that uniaxial materials tests only are not sufficient to fully characterize the materials response [Bass et al., 2004] they are widely used for materials characterization. Starting with Wertheim [1847] the experience in characterizing biological tissues in tension grew steadily leading to a well-understanding of this type of test. Therefore it is not surprising that tensile experiments are also nowadays used in biomechanics where a variety of tissues, including anulus fibrosus (AF) tissue, are characterized uniaxially. Section 6.1 gives an overview on the modes of deformation present in the AF and motivates the experiments chosen to achieve. Also the reasons not to perform biaxial or pure shear experiments are discussed. The following sections describe prestudies that were conducted, the sample preparation, the experiments and eventually the results. 57

CHAPTER 6. UNIAXIAL EXPERIMENTS WITH ANULUS FIBROSUS SAMPLES

6.1

Deformation modes in the intervertebral disc

As described in Section 2.3, the intervertebral disc (IVD) can be considered, in a first approximation, as a vessel (the AF) filled with an incompressible medium (the nucleus pulposus (NP)). Under an axial load the height of the IVD decreases, the incompressibility of the NP forces the AF to bulge radially outward in order to maintain its volume (cf. Fig. 6.1). In practice things are more complicated, just to mention a little of the complexity: the IVD is not cylindric nor has it a constant height, the AF has a varying radial thickness and its tissue is strongly anisotropic. Though the question is what the stresses and strains look like in a real geometry with (approximately) the correct mechanical properties. In anticipation of Chapter 9, where the finite element model is presented in more details, the results of some simulations shall be presented. For different anatomical motions of a (ovine) functional spinal unit (FSU) the strain and stress fields were computed. The geometry of the model originates from a computer tomography scan of a real FSUs. The vertebrae (Ve) are modeled as rigid bodies. The AF is divided into ten circumferential layers for which the elastic modulus of the collagen fibers increases linearly from zero at the innermost layer to their largest values at the outermost layer. The fibers include with the horizontal, circumferential direction angles of ±30◦ . The constitutive model and the according properties of the AF fibers and the ground substance are taken from literature [Eberlein et al., 2001]. In short, the AF is modeled as an incompressible, anisotropic, hyperelastic material and the NP as an incompressible fluid. For more details please refer to Chapter 9. Figure 6.2 displays the modes of deformation in terms of the invariant K3 (see App. C.3 for more details about K3 ) for each of the motions axial compression (b), axial torsion (c), flexion (d), extension (e), and lateral bending (f). The geometry was taken from an ovine FSU, therefore the amplitudes of deformation are adjusted to the range of motion of the sheep spine (cf. Tab. 2.1). As it can be seen, the situation is totally different from the ideal case of a circular AF. For a such under an axial load one would expect a circumferential zone of uniaxial extension at the outer AF where the material is stretched in peripheral direction while being compressed in axial and radial direction. Figure 6.2b clearly shows that this is not true in an actual AF where a multitude of different deformation modes, ranging from uniaxial to biaxial extension, applies. Though it is clear that a characterization of the mechanical properties of the AF tissue should include all those deformation modes. For experiments with AF samples three setups that activate different deformation modes were evaluated. This are i) biaxial experiments with extension in circumferential and axial direction, ii) pure shear experiments with a biaxial extension in circumferential and radial direction and iii) uniaxial experiments in circumferential direction. Only uniaxial experiments were eventually realized, the first two methods all showed severe practical limitations that restrained their use, as explained hereafter:

6.1. DEFORMATION MODES IN THE INTERVERTEBRAL DISC

unloaded IVD A

IVD under axial load

σa

A p

σθ

Section A – A

σθ Figure 6.1: Schematic drawing of a quarter section of an IVD: unloaded (left) and under axial load (right). The external axial stress σa is counterbalanced by an internal pressure p (top right). Under the pressure p the AF bulges radially inducing circumferential stresses σθ (bottom right). (The proportions of the stress values are taken from Nachemson [1966a].)

Biaxial experiments are a popular complement to classic uniaxial materials tests from which is known that they only are not sufficient to fully characterize the materials response [Bass et al., 2004]. Bass et al. used biaxial materials experiments to characterize the constitutive response of human AF tissue with some success. They used cruciform shaped specimen where in (anatomical) axial direction the vertebral bone was retained such that the bone could be clamped. The acquisition of a biaxial materials testing machine for our laboratory led to a study about the optimal specimen shape (cf. Chap. 5). From this study it got clear that simple cruciform shaped specimens are not sufficient to get meaningful results. More engineered shapes should be used in order to get appropriate responses that can clearly be attributed to the biaxial materials response. This results and the smallness of the possible AF samples led to the conclusion that biaxial experiments were not feasible within this work. Pure shear experiments require a much simpler equipment than the classic biaxial tests to produce a multiaxial stress response as well [cf. Hollenstein et al., online first]. 59

CHAPTER 6. UNIAXIAL EXPERIMENTS WITH ANULUS FIBROSUS SAMPLES

uniaxial extension

0.5 0

undeformed / pure shear

−0.5 (a) undeformed

−1

biaxial extension

(b) axial compression

(d) flexion

(e) extension

(f) lateral bending Figure 6.2: Modes of deformation for different postures, predicted by finite element simulations. Since the AF is (almost) symmetric with respect to the median plane only half of it is shown. Undeformed geometry (a), 10 % axial compression (b), 1◦ axial torsion (c), 4◦ flexion (d), 4◦ extension (e) and 4◦ lateral bending (f ). 60

6.2. PRE-STUDIES

A materials testing machine with only one axis of actuation, as for uniaxial materials testing, suffices. In contrast to uniaxial experiments where, in order to allow free contraction in the sample center, long and thin samples are clamped at their long ends, in pure shear tests short and thin but broad samples are clamped along their long side. This is done in order to prohibit lateral contraction as much as possible such that a biaxial stress state develops. Long samples of the AF can only be obtained with the long side along the circumferential direction. Though, for the short and thin dimension only the axial and radial directions are left. It is not appropriate to pull along the radial direction since laminae from the AF would separate. Therefore the (anatomical) axial direction had to be the direction of load application. But as the height of the AF is in the order of some millimeters only this dimension would be very small. Clamping problems, especially slippage, would influence the experimental results. This major problem would be accompanied by problems with sample preparation and, as the samples would need to be extremely thin, the difficulty to prevent dehydration.

Uniaxial tensile tests on AF strips are motivated by findings of Shah et al. [1978] who demonstrated experimentally the existence of tensile circumferential strains on the AF surface. An argument that justifies the number of mechanical characterizations of the AF-tissue in this direction.i The AF is, as described in Chapter 2, build up by several concentric layers (lamellae) of a transversely isotropic material. Holzapfel et al. [2005] find “that the single lamella is the elementary structural unit of the anulus fibrosus” and study consequently the unilamellar mechanical properties in fiber direction, as did Skaggs et al. [1994]. The low number of studies that used unilamellar samples (with Panagiotacopulos et al. [1979] only three were found) might be explained by the difficulty to separate the individual layers; most studies are done with multilamellar samples. Due to the problems related with the realization of biaxial or pure shear experiments it was decided to relinquish these and to concentrate on uniaxial experiments with multilamellar AF strips cut in circumferential direction.

6.2

Pre-studies

Before the actual experiments were realized four pre-studies were carried out with the aims to: i) quali- and quantify the errors admitted if the true (quasi-uniaxial) experiments are approximated by purely uniaxial stress states, ii) estimate the error in the strain calculation if the samples are not perfectly aligned with the tensile axis, iii) measure and verify the compliance of the force transducer and iv) qualify the changes induced by the freezing process used for storage and sample preparation. i

Ebara et al. [1996]; Elliott and Setton [2001]; Galante [1967]; Hirsch and Galante [1967]; Panagiotacopulos et al. [1987]; Wagner and Lotz [2004]; Wu and Yao [1976]

CHAPTER 6. UNIAXIAL EXPERIMENTS WITH ANULUS FIBROSUS SAMPLES

6.2.1

The influence of the samples aspect ratio

The question is what the error will be if a real experiment, where the sample is clamped and lateral contractions at the grips are prohibited, is interpreted as a perfect uniaxial experiment where only axial stresses occur. This was analyzed by taking into consideration a hom*ogeneous, isotropic sample with a given initial length L0 and cross-section A0 (cf. Fig. 6.3). From the experiment only the initial free clamp length L0 , the clamp displacement ∆L and the clamp force F ∗ are known. If a uniaxial stress state is assumed, the strain ε = ∆L/L0 is constant over the length of the sample. In the case of a quasiuniaxial stress state the strain is no longer constant but depends on the position. Near the clamps the axial strain, due to the constraint lateral contraction, will be smaller than ε. In order to get a total change in length of ∆L the strain at the center ε∗ has to be larger than ε. This means that, if the real quasi-uniaxial test is assumed to be a uniaxial one, the true strain ε∗ at the center will be underestimated. If, as in the present work, calibrations of constitutive models are done by means of displacements and forces, it is of interest to know what the forces are in comparison with the uniaxial case. The total energy of the quasi-uniaxial system is higher than of the uniaxial system. This can easily be seen by the fact that the uniaxial system can be transformed into the quasi-uniaxial system by additional lateral forces. As the total energy of the systems is given by the outer forces F and F ∗ , respectively, times the displacement ∆L the force F ∗ has to be higher than F . Therefore, if the experiments are calibrated with an uniaxial model, the material stiffness will be overestimated. Given a constant cross-section, the influence of the clamping will become smaller and smaller if the sample becomes longer. Though, as the aspect ratio (length:width) grows the closer the two systems get and the smaller the error made by the assumption of an uniaxial stress state becomes. It is clear that the precise errors will depend on the material itself. To get an idea of the order of magnitude, finite element simulations were run assuming an incompressible Neo-Hookean material. The aspect ratio of the samples was varied between 1:1 and 10:1. As can be seen in Figure 6.4a the true stretch at the center λc can reach much higher values than would be assumed by the stretch derived from the grip to grip separation λg . For higher aspect ratios the error becomes smaller, reaching 1 % in maximum for the aspect ratio 10:1. The same holds true for the forces (cf. Fig. 6.4b), where at an aspect ratio of 10:1 the total force in a quasi-uniaxial experiments is only 1 % higher than the force would be for an uniaxial experiment at the same stretch level.

6.2.2

The alignment of the sample

The samples alignment with the tensile axis has an influence on the true stretch “felt” by the sample. Lets assume that one end of the sample is clamped with an offset d from 62

6.2. PRE-STUDIES l0 A0

undeformed L0

l uniaxial

A L0

F ∆L

l∗ quasi-uniaxial F ∗

A∗ L0

F∗ ∆L

Figure 6.3: Comparison between two deformations where in one acts a uniaxial stress state while in the other a quasi-uniaxial stress state is present.

the tensile axis that is defined as going through the clamped area at the other side of the sample (Fig. 6.5). If the clamps have an initial grip to grip separation of L0 and after displacement a separation of L then the computed (nominal) sample stretch would be λn = L/L0 . The (true) stretch experienced by the sample λt = l/l0 is smaller, depending on d or α0 , respectively. The relation between the nominal and the true stretch is given by the formula 2 1/2 1/2 2 L + d2 λn + tan2 [α0 ] λt = . (6.1) = L20 + d2 1 + tan2 [α0 ] As can be seen from Figure 6.5b the error is small, even for unrealistic large misalignments of α0 = 20◦ . In real experiments offsets of 1 mm might be possible, with a free sample length of about 15 mm, as in the later experiments, the error is smaller than 1 %.

6.2.3

The compliance of the force transducer

The compliances of the force transducers were shown to be non negligible for two of the three used force transducers. The 10 N force transducer (AST KAW-S 10N, cf. App. D) that was used for the qualification of the effects of freeze storage has a specified nominal measurement displacement of 0.3 mm. A check with a laser interferometer (opto NCDT 2000-10, cf. App. D) emphasized that this value is about 2.5 times too low, the true nominal measurement displacement was found to be around 0.865 mm, i.e. a compliance of 0.0865 mm/N. The values of the 50 N force transducer (AST KAP-S 50N, cf. App. D) instead, for which the data sheet declares a nominal measurement displacement of 0.25 mm, i.e. a compliance of 0.05 mm/N, was verified. 63

CHAPTER 6. UNIAXIAL EXPERIMENTS WITH ANULUS FIBROSUS SAMPLES

1.16

1:1 2:1 3:1 5:1 10 : 1

1.14 1.12

λc /λg [-]

1.1 1.08 1.06 1.04 1.02 1 1

1.1

1.2

1.3

1.4

1.5 λg [-]

1.6

1.7

1.8

1.9

(a) Ratio λc /λg

1.25

1:1 2:1 3:1 5:1 10 : 1

Fq /Fu [-]

1.2

1.15

1.1

1.05

1 1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

λg [-] (b) Ratio Fq /Fu Figure 6.4: Visualization of the errors made if a quasi-uniaxial tensile experiment is approximated by a purely uniaxial one: ratio of the true stretch at the center λc and the stretch derived from the grip to grip separation λg (a) and ratio of the total force in a quasi-uniaxial tensile test Fq and the total force in a purely uniaxial one Fu (b).

6.2. PRE-STUDIES

1 α0

0.99 λt /λn [-]

0.98

0.97 L L0

α0 α0 α0 α0 α0

0.96

= 0◦ = 5◦ = 10◦ = 15◦ = 20◦ 1.2

1.4

1.6

1.8

λn [-] (a) Sketch

(b) Ratio λt /λn

Figure 6.5: Error in the calculation of the (true) stretch if the sample is not properly aligned with the tensile axis: sketch of the situation (a) and graph with the ratio of the true stretch λt and the nominal stretch λn (b).

Furthermore, these measurements showed that the compliance is a constant in good approximation and can be corrected. The grip to grip separation from which the sample stretches were computed are corrected by the force times the corresponding compliance. This can clearly be seen in Figure 6.6a where at the relaxation phases, not only the force decreases but the stretch increases. This increase in stretch comes from the fact that the samples relax, and therefore smaller forces act on the force transducer such that it retracts a little bit. This means, on the other hand, that the relaxation experiments were done with a (slightly) changing stretch.

6.2.4

Effects of freeze storage

Freeze storage that lasts for some days or months is common practice and helps to maintain the in-situ water content during storage [Elliott and Setton, 2001]. A second reason to congeal samples is their preparation. The softness of many biological tissues is a problem when well-defined sample geometries have to be cut; solid tissues are much easier to shape. Only one study was found where fresh specimens were used [Holzapfel et al., 2005]. Freezing might damage tissues due to ice formation as shown by Rubinsky et al. [1990] at the example of liver tissue. Effects of freeze storage for IVD tissues and FSUs are 65

CHAPTER 6. UNIAXIAL EXPERIMENTS WITH ANULUS FIBROSUS SAMPLES

controversially discussed, some studies report no influence of freezing on the measured mechanical properties while others found differences. Hirsch and Galante [1967] tested AF samples in uniaxial tension with the conclusion that freezing under dry ice snow (for short periods) and thawing at 100 % humidity preserves the mechanical properties. Studies that analyzed FSUs that were freeze stored for up to seven months did not find any changes in the mechanical responses [Gleizes et al., 1998; Panjabi et al., 1985]. Hickey and Hukins [1979] investigated with X-ray diffraction the effect of different types of freezing (immersing specimens in liquid nitrogen and putting them in a freezer) on the arrangement of collagen fibrils in the AF. They could not detect any disruptions in the fiber arrangement. In contrast, Evans and Lissner [1959] found that freezing might stiffen tissues. Bass et al. [1997] tested the effect of freezing on the long term creep response of FSUs and found permanent and not recoverable changes for the frozen samples. Nonetheless it was found that the elastic behavior was not significantly affected. The measured differences to previous studies are attributed to the fact that in their study the long term response was measured and not responses for short terms. Due to this controversial discussion about the effects of freezing it was decided to qualify the changes by own means. This was done by measuring AF tissue samples in uniaxial tensile experiments before and after freezing such that any existent differences should be visible in the force-displacement diagrams. The same samples were characterized before and after freezing. As cutting in the fresh state is unhandy, samples were tested as they were dissected from the spine (cf. Fig. 6.7b). Though, their geometry was not well defined and results between different samples cannot be compared, but still, a comparison of the results of one sample before and after freezing is possible. In order to minimize the influence of clamping and unclamping samples were gripped only once and frozen with the clamps. The preparation and experimental protocol was the following: Dissection of the samples from the spine and subsequent cutting in a roughly oblong shape. This samples were gripped with custom made clamps. To ensure a good positive locking the inner sides of the clamps were faced with rough sandpaper. Clamping forces such that the samples were slightly squeezed were applied. Samples with clamps were mounted in the materials testing machine and measured twice. Displacement ramps with 0.1 mm/s were applied until a maximum force of 5 N, once this force was reached a relaxation phase of about 60 s followed. For a 15 mm long sample the speed of 0.1 mm/s results in a stretch rate of approximately 0.7 %/s what is within the range of strain rates reported in literature where rates from 0.009 %/s up to 2 %/s are reported [Holzapfel et al., 2005; Skaggs et al., 1994]. After the first run samples were allowed to relax for about 4 minutes at an undeformed state before the second run was started. Samples were wetted with saline, using a pipet, before and after each run. After two consecutive runs the samples and the clamps were deep frozen at −30 ◦C for about 30 minutes (one sample was frozen for even 5 days). In order to prevent dehydration, specimens were sealed in plastic bags. 30 minutes were more than enough to freeze the samples to the core. Following, the samples were thawed in a saline bath at room temperature and tested again twice. 66

6.3. SAMPLE PREPARATION

As all samples had different lengths, and therefore the force-displacement characteristics could not be compared directly, the sample stretch was computed. In order to do so, the initial lengths of the samples were determined in the postprocessing. A limit force of 0.01 N, that is thought to be negligible, was defined, the sample length at this load was taken as initial length. The results of the five specimens that were characterized with this protocol are shown in Figure 6.6. Since the geometry of the specimens are neither constant in shape nor in geometry results needed to be somehow normalized to be comparable between the specimens. The result of the first run of each sample was scaled such that the maximum force and the related stretch or time points coincide with the mean maximum force and related mean stretch and mean time points. The three subsequent results were scaled with the same factors. In this way it is possible to compare qualitatively relative changes due to the freezing and unfreezing process. As can be seen in both sub-diagrams (a) and (b) of Figure 6.6 no clear difference in the responses of fresh and unfrozen specimens can be detected. This is true even for the sample that was frozen for 5 days. If differences exist, the are hidden in the general scatter. Therefore it was concluded that freezing for storage and preparation purposes does not influence the mechanical behavior in an unacceptable way.

6.3

Sample preparation

Fresh lumbar sheep spines were obtained from the Veterinary Clinic of Zurich, Switzerland, few hours after euthanasia. Tissue harvesting, preparation and experimentation was done as soon as possible after exitus in order to minimize the post mortem changes as it is known that biological tissues and especially soft biological tissues lose their in-vivo characteristics and degenerate after exitus [Fitzgerald, 1975]. For IVDs, that contain mainly water and proteoglycans, sources of degeneration are dehydration and the breakdown of the proteins. Both processes can be decelerated or even stopped by keeping the samples moist and cool. The spines were often in their surrounding muscular tissues, rarely the butcher already dissected the spinal column. The lumbar spines were brought to our laboratory and, if necessary, dissected by use of a scalpel. Soft tissue adhering to the IVD was removed until the diagonal orientation of the AF collagen fibers got visible. Afterwards, anterior portions of the IVD were excised with a scalpel (cf. Fig. 6.7a). To do so, cuts were made in parallel to the caudal and cranial endplates of the neighboring Ve, this procedure mostly resulted in wedge like pieces with AF tissue at the thicker (outer) side and some NP tissue at the thinner (inner) side. The dispensable NP tissue was cut away. Due to the kidney shape of the IVD cross-section these wedges were rounded at their outer side (cf. Fig. 6.7b). In the following the samples should be die-cut to rectangular samples with the long side in parallel with the lamellae. In order to facilitate the stamping samples were stretched and pined with needles onto a piece of styrofoam. The styrofoam was packed in a sealed plastic bag and frozen at −30 ◦C until experimenting. A custom built stamping instrument was used to die-cut the samples (cf. Fig. 6.7c). It was designed such that two razor blades can be arranged in parallel with a given 67

CHAPTER 6. UNIAXIAL EXPERIMENTS WITH ANULUS FIBROSUS SAMPLES

fresh unfrozen

force [N]

4 3 2 1 0 0.9

0.95

1.05 stretch [-]

1.1

1.15

(a) force vs. stretch

fresh unfrozen

force [N]

4 3 2 1 0 −20 −10

100

110

time [s]

(b) force vs. time Figure 6.6: Comparison of the mechanical response of samples that were tested fresh and unfrozen: force vs. stretch (a) and force vs. time (b). A • marks the results of a sample that was frozen for 5 days.

6.4. EXPERIMENTS

distance in between that can be adjusted by means of intermediate flat blanks. Diecutting was done minutes before measuring when the samples were still in a (semi-) frozen state. The samples were first stamped in the radial-circumferential plane, turned by 90◦ , and then die-cut in the axial-circumferential direction. In this manner specimens with square cross-sections were produced. If the samples were still deep-frozen during the stamping process, sometimes it happened that the razor blades bent outward, resulting in unprecise cross-sections. Due to the small size of the samples (the samples had cross-sections of 1.5 × 1.5 mm2 ) unfreezing was very fast and minutes after removal from the deep-freezer they could be characterized.

6.4

Experiments

Tensile experiments were done at the Balgrist University Hospital (Zurich, Switzerland) on a standard uniaxial materials testing machine (Zwick 1456, cf. App. D) equipped with a 50 N load transducer (AST KAP-E 50N, cf. App. D). Displacements were measured by means of an external extensometer (Zwick Multisens, cf. App. D) for reasons of accuracy (cf. Fig. 6.7d). The lower arm of the extensometer was attached to the lower clamp and the upper to the upper part of the force transducer. It could not be attached to the upper clamp directly since the extensometer would disturb the force reading. Though it was necessary to correct the recorded grip-to-grip distance by the compliance of the force transducer. The raw samples were taken out of the freezer and cut into strips with 1.5 × 1.5 mm2 cross-section,1 as discussed in the previous section. These strips were mounted in custom made clamps, the compressive force was chosen such that the samples were slightly squeezed; sandpaper glued to the clamp surfaces increased the force locking. Subsequently the sample and clamps were mounted in the materials testing machine. Right after, and at all the following times, when no measurements were made, the sample was dipped with saline such that it was always covered by a fluid film. The experimental protocol was designed to include three displacement driven load ramps with a specific strain rate up to given strains and 90s relaxation just afterwards. After the last relaxation the jaws were closed such that the sample could relax for 2 minutes. Due to its viscous nature the tissue will regain its initial properties with time [Galante, 1967]. This protocol was driven five times for each sample. After the experiments the samples were freeze stored once more until the determination of their cross-sections. The five cycles were taken as preconditioning phase. After the first four cycles the 1

For hom*ogeneous and isotropic materials exist standards that prescribe the exact shape of the specimens [e.g. EN ISO 527, Europ¨ aisches Komitee f¨ ur Normung, 1995 – 1997]. For several advantages, these material classes are often tested with dumbbell shaped specimens: 1) clamping forces are well distributed due to the wider ends, 2) stresses are concentrated at the central part where a uniform stress state dominates and 3) samples fail at the center part. Volumes four and five of standard EN ISO 527 cover isotropic, anisotropic and unidirectional fiber-reinforced synthetic composite materials. For those classes of materials rectangular shaped specimens are prescribed.

CHAPTER 6. UNIAXIAL EXPERIMENTS WITH ANULUS FIBROSUS SAMPLES

(b) Excised sample from the AF

(a) Spine where the samples were excised from

Figure 6.7: Detail of a spine section where the samples were excised from (a). A typical sample before it was pinned on styrofoam (b). Custom made stamping device (c). Sample clamped in the materials testing machine: custom made clamps, load cell and extensometer (d). (d) AF sample mounted in the materials testing machine (photo taken by J. Hunger) 70

6.5. RESULTS

material was said to be preconditioned such that the last cycle would represent the “trained” material response. In order to determine the initial length of the samples, they were, after being mounted in the materials testing machine, stretched up to a force of 0.01 N. The clamp to clamp distance at this force level was taken as the initial sample length L0 . As the control software of the testing machine does not support strain rates as driving variables but only velocities, the initial sample length was used to calculate the needed grip displacement velocity such as to obtain the desired strain rate of 1 %/s: v = ε˙ × L0 . The initial sample length determined as well the three holding points (given as absolute positions) at 12.5, 15 and 17.5 % strain, respectively. First tests runs with this protocol showed that the samples did, after one cycle, not return to their initial length but stayed longer. Though for the final experiments it was decided to account for this changes in length by measuring the sample length for each cycle a new, including a change of the velocities and the holding points. The strain rate of 1 %/s was chosen according to the literature in order to render the the experiments quasi-static. Fujita et al. [1997] performed pilot studies with strain rates ranging from 0.005 %/s to 50 %/s and did not find any dependency of the stressstrain behavior on the strain rate. This findings are second by Holzapfel et al. [2005] who reports only a small dependency of the tensile behavior on the strain rate. Cross-section areas were determined under a microscope which can take digital images (cf. App. D). From the frozen samples five cross-sections were cut at different locations that were stored in saline. Just before taking a picture under the microscope the sections were taken out of the saline and their surface dried. This had to be done quickly since the small portions dehydrated very fast. The software AxioVision allows the measurement areas in the pictures taken. The five measurements were averaged such that a mean cross-section area was determined.

6.5

Results

The mean free sample length (±standard deviation) averaged over all measured samples and all cycles was 15.6±2.3 mm. The cross-sectional areas measured by use of an optical microscope were much larger than envisaged (cf. Fig. 6.8). More importantly the areas were neither constant nor were they square shaped. Nevertheless, the variabilities in sample length and cross-section are comparable with the ones reported in literature [Elliott and Setton, 2001; Fujita et al., 1997; Holzapfel et al., 2005]. Aspect ratios greater than 10:1 were realized and thus the approximation of the tests with an uniaxial stress state are legitimate. In total six samples were characterized with the protocol described in Section 6.4. One of the experiments had to be discarded due to loose clamping which resulted in slippage of the sample. Figure 6.9 shows the load-elongation and load-time curves for one selected sample. 71

CHAPTER 6. UNIAXIAL EXPERIMENTS WITH ANULUS FIBROSUS SAMPLES

16 12.65 (±1.35)

cross-sectional area [mm2 ]

14 12

9.20 (±1.50)

7.97 (±0.85) 6.66 (±0.97)

7.14 (±0.66)

7.18 (±0.28)

6 4

3 4 sample number

Figure 6.8: Cross-sectional areas of the characterized samples; for each sample five crosssections were sized. The values on top indicate the average (±standard deviation) areas.

6.6

Discussion

The experiments were effectuated on AF samples with approximately 15 mm length and nominal cross-sections of 1.5 × 1.5 mm2 . In a step after the experiments the actual cross-sections were determined under an optical microscope. The areas found were all much larger than expected, ranging from 6.66 mm2 up to 12.65 mm2 , and not square shaped. Several sources were identified that could be part of the problem: 1) the two coplanar razor-blades of the stamping tool bent outward while cutting, especially if the samples were still deep frozen (observed at the sample preparation), 2) the film of saline covering the samples during experimentation might have led to swelling [Hirsch and Galante, 1967] and 3) swelling might have taken place under the microscope. Although the cross-sectional areas were larger than planned, aspect ratios of more than 10:1 were obtained. Consequently, errors in the interpretation of the data by assuming an underlying uniaxial stress state are small. In other studies five [Holzapfel et al., 2005] to ten [Wagner and Lotz, 2004] preconditioning cycles were used to assure that the specimens are well gripped and different samples have the same load history (at least for a short time in the past). In the present study five cyclic loadings were applied to all samples. Despite some convergence can be observed, no steady state reactions were found at all. It seems that more load cycles would have been needed to fully precondition the samples. 72

6.6. DISCUSSION

force [N]

cycle cycle cycle cycle cycle

1 2 3 4 5

6 4 2 0 1

1.02

1.04

1.06

1.08 1.1 stretch [-]

1.12

1.14

1.16

1.18

(a) force vs. stretch

force [N]

cycle cycle cycle cycle cycle

1 2 3 4 5

6 4 2 0 0

100

150

200 time [s]

250

300

350

400

(b) force vs. time Figure 6.9: Results from cyclic uniaxial tensile experiments on one sample: force vs. stretch (a) and force vs. time (b).

CHAPTER 6. UNIAXIAL EXPERIMENTS WITH ANULUS FIBROSUS SAMPLES

While the recorded load curves are qualitatively very similar they differ a lot in absolute values. Not only the forces are different but the nominal stresses obtained by dividing the forces by the respective cross-sectional areas as well. Finally the small number of characterized samples makes it impossible to determine representative values, e.g. for the stiffness. Therefore a measurement that is thought to represent a typical behavior was chosen. The parameter determination explained in Chapter 9 was made based on the data from this one sample. Some observations that can be made concern the presence of - a pronounced initial toe-region, - a distinct hysteresis and - an increase in stiffness from one cycle to the next. At small deformations the load-elongation curves are characterized by an initial nonlinear toe-region that is followed by a quasi-linear region. Such traits were already described by many others for AFi and other collagenous tissues. This behavior can be attributed to two different load mechanisms: in an initial phase the fibers are crimped and mainly the matrix withstands the load while at higher strains the fibers are stretched and in tension [Elliott and Setton, 2001]. As Wagner and Lotz [2004] point out, the toe region might be important for smooth passage from the compliant compression regime to the stiffer behavior in tension. Integration of the area under the load curves showed that the energy dissipation increases with the cycle number (only for one sample the dissipation was lower during the fifth cycle than for the two previous ones). The increase in stiffness from one cycle to the next can also be seen in the forcetime curves shown in Figure 6.9b. As well recognizable in this graph is the relaxation behavior: while for the first cycles and for smaller loads 90 s are sufficient for the material to relax, longer relaxation phases would be needed to reach equilibrium for higher load at later cycles.

Ebara et al. [1996]; Elliott and Setton [2001]; Fujita et al. [1997]; Skaggs et al. [1994]; Wagner and Lotz [2004]

Chapter

Experiments with Functional Spinal Units Ovine functional spinal units were characterized in longitudinal loading and bending. The complex geometry of the vertebrae necessitated the design of special mountings in which the vertebrae are ingrained. Validation of the mounting method and the used materials testing machine provided valuable information for improvements of the experimental equipment and procedures. The final experimental protocols comprised cyclic tests with different load rates to give insight into the viscoelastic properties of the tissue. Cyclic tests were repeated several times in order to check for repeatability. Results reveal typical characteristics for viscoelastic materials as are nonlinearity, distinct dissipation and ratcheting.

Realistic simulations of organs and other biological tissues and structures are of increasing importance for medical therapies and training. The key hereto is the qualitative and quantitative understanding of the mechanical properties. In the last decades many constitutive models for soft biological tissues were proposedi , some of them are rather complex and use a large number of parameters. A way to assess the predictive power of the various models is the comparison with sound experimental data. In this work the viscoelastic properties of ovine lumbar intervertebral discs (IVDs) in axial loading and in bending was determined. The complex geometry as well as the small imposed displacements combined with rather high forces and moments are to be treated. The two vertebrae (Ve) adjacent to the IVD are mounted using custom made clampings where attention is paid not to prestress the tissue. Extensive preliminary experiments were run with steel springs in order to check the suitability of the used materials testing machine. No significant compliance of the machine could be detected but a negative hysteresis in the response to cyclic axial exitation. First experiments with IVDs gave an in-depth insight into the underlying mechanics such that appropriate experimental protocols could be defined. i

Arruda and Boyce [1993]; Ehret et al. [2010]; Elliott and Setton [2000]; Holzapfel et al. [2000]; Merodio and Ogden [2005]; Rubin and Bodner [2002]; Weiss et al. [1996]

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

Eventually six functional spinal units (FSUs) were characterized in longitudinal loading (compression and tension) and one in bending (flexion and extension). The obtained results give insights into the viscoelastic properties; cyclic experiments with different load rates reveal a distinct dependency on the strain rate. In Chapter 9 the data sets gained by these experiments are used as test cases for a finite element (FE) model calibrated with data from uniaxial experiments on anulus fibrosus tissue. Parts of what is presented in this chapter originate from a collaboration with M. Hollenstein and O. Papes from the same laboratory and were presented at the 80th Annual Meeting of the International Association of Applied Mathematics and Mechanics (GAMM) (Gdansk, Poland, 2009). Thereof resulted a publication in the Proceedings in Applied Mathematics and Mechanics (PAMM) [Helfenstein et al., 2009b].

7.1

Preliminaries

Several steps had to be taken before starting with the experiments with FSUs. The first was to think of an appropriate mounting system that - assures a good alignment of the samples, - has a higher stiffness than the samples to be characterized and - is free of any backlash. To validate the designed mountings experiments with steel springs and on FSUs were carried out. The resulting experiences led in several iterations to the final setup and protocols including improvements of the mountings, the technique for the displacement measurement and the protocol for sample preparation and experimentation.

7.1.1

The mounting

The experiments were performed on a standard uniaxial materials testing machine (Zwick 1456, cf. App. D) which offers plug and socked connections that can be used for fixation of (custom made) mountings (cf. Fig. 7.1a). This connections consist of a plug on which the mounting can be put, a lateral fixation pin and a knurled nut to exclude any backlash from the system. Depending on the force range to be measured different interface sizes exist. In several studies in which FSUs were tested in compression only, the upper and lower Ve were cut to form two flat surfaces that were be used for the load application [e.g. Virgin, 1951]. As in the present experiments tension is considered as well another solution had to be found. Inspired by the literaturei , the mounting of the FSU is based on the idea of ingraining the Ve in polymethylmethacrylate (PMMA) (cf. App. D). Hollow cylinders that can be mounted on the testing machine (Figs. 7.1 and 7.2) are used as pick-up. In order to increase the form closure between V and PMMA three wood screws are screwed into the free vertebral endplates prior ingraining. i

e.g. Bass et al. [1997]; Costi et al. [2002]; Heuer et al. [2007]; Pitzen et al. [2002]

7.1. PRELIMINARIES

fixation pin plug knurled nut

(a) Version 1

(b) Version 4

Figure 7.1: Drawings of the first (a) and the fourth (b) version of the mountings. Figure (a) shows as well the (small) plug and socket connection of the materials testing machine (including plug, fixation pin and knurled nut).

vertebra steel ring set screw wood screw PMMA

IVD

mounting

vertebra steel ring

machine interface

(a) Sketch

set screw mounting (b) Detail view

Figure 7.2: Custom made clampings are used for the experiments with FSUs. Sketch of the mounting principle (a): three wood screws are screwed into the endplates of the Ve, screws and bone are ingrained using PMMA. A steel ring is bolt from the top to pretension the polymer. On the right a detail view of a mounted FSU (b).

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

A drawing of the first version of these mountings is shown in Figure 7.1a. Basically it consists of an aluminum cylinder that has on its one end a cavity in which the V is ingrained (cf. Fig. 7.2, which shows the final mountings version). Four lateral threads allow the fixation of set screws that hold the PMMA in place, which was thought to be important for tension experiments. The bottom of this cavity has three threads that are used to eject the sample after ingraining or experimentation.

7.1.2

Validation experiments

With steel springs The compliance of the testing system was checked using steel springs (cf. App. D) in axial loading experiments. Their stiffness (709 mm/N) was chosen such that it is similar to the stiffness of an ovine FSU (determined in a first experiment to be around 955 mm/N).1 Two types of testpieces were used for those validation experiments. To measure the compliance of the materials testing machine only, a spring was welded on steel adapters (Fig. 7.3a). To test the whole system with the PMMA a spring was ingrained using intermediate steel parts welded to it (Fig. 7.3b). With this arrangement the compliance of the PMMA in series with the materials testing machine was measured. With both testpieces cyclic uniaxial experiments were run with cycles centered at −450 N (compression) and peak to peak amplitudes of 900 N in a displacement controlled mode at a rate of 2 mm/s. The validation experiment with a spring mounted on steel adapters allows to attribute any measured compliance to the testing machine. Linear interpolation of cyclic experiments reveals a very good accordance of the measured spring constant with its declared value (708 mm/N measured vs. 709 mm/N declared). The materials testing machine is therefore supposed to have no significant compliance in the tested force range from 0 to −900 N. A comparison of the measured minimum/maximum reversal points with the targeted values (−903 N/3 N measured vs. −900 N/0 N defined) shows the ability of the controller to accurately drive the system with displacement rates up to 2 mm/s at the given stiffness. The experiments with the steel spring mounted on steel adapters and the ingrained steel spring gave almost identical results. The conclusion thereof is the the PMMA does not add a significant compliance to the testing system nor does it have any other influence on the measured results. Figure 7.4 shows the result of the cyclic experiment with the molded spring. The ten cycles reveal a hysteresis with a width of 6 – 10 µm (corresponding to 4.25 – 7.08 N). A closer look at the force-displacement data reveals that the force signal lags behind the displacement signal by 0.15 s, indicating a negative hysteresis. A negative hysteresis suggests that the sample generates energy every loading-unloading cycle, something that is known to be non-physical (in contrast to a positive hysteresis that indicates a 1

In fact, the true stiffness of the ovine FSU might even be higher. Smeathers and Joanes [1988] report a reduced compressive stiffness of FSUs due to freezing that might relate with a change of the osmotic pressure of the tissue.

7.1. PRELIMINARIES

(a) steel spring welded to steel adapters

(b) steel spring ingrained in PMMA Figure 7.3: Samples for the validation of the materials testing machine: steel spring welded to steel adapters (a) and steel spring ingrained in PMMA (b).

dissipative process). The consequence of this finding is discussed in the second paragraph of Subsection 7.1.3. With a functional spinal unit For the experiments ovine lumbar spines were obtained within 2 hours post mortem from the veterinarian hospital (University of Zurich, Vetsuisse Faculty). In our laboratory the Ve LV and LVI with intervenient IVD were dissected, surrounding muscular and ligamentous tissues and posterior elements were removed. One V was ingrained in our laboratory. To do so the cranial vertebra was held by a burette clamp such that the transversal plane of the IVD was as horizontal as possible. The PMMA was filled into the cylinder and let harden. Once the polymer hardened, the sample was ejected from the mounting, sealed in a plastic bag and stored at −30 ◦C until testing. Four days before the experiments the sample was taken out from the freezer and allowed to unfreeze in a refrigerator. The evening prior measuring the sample was mounted at the upper interface of the materials testing machine and the lower V was ingrained. The moulding is done directly in the materials testing machine such that any misalignment is excluded and the induced pre-loads are minimized. After the time to harden the polymer, the position of the machine traverse was adjusted such that the displayed force equaled zero. This position was kept overnight until the experiments started the following morning. In order to keep the samples wet, they were loosely wrapped in saline soaked gauze covered by cling film (a measure that is also used by others, e.g. Huber et al. [2003]). All manipulations and experiments were done at room temperature. Studies on humans and baboons showed that peak forces up to 1.5 times the body79

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

weight are supported by the FSU without any damages;i for sheep this suggests peak forces down to −900 N. Therefore test protocols were designed such that this limit is not undershot. With FSUs from different sheep the following types of experiments were performed: - cyclic experiments centered at −450 N with peak to peak amplitudes of 900 N and a constant displacement rate of 2 mm/s, - creep experiments at a force of −450 N, with a duration of 2 h, - relaxation experiments with a peak force of −900 N and a duration of 1 h (beginning after the creep experiments at −450 N). Experiments are done at room temperature. In order to keep them moist they are wrapped in a saline soaked gauze. Figure 7.5 shows the results of a cyclic experiment with a FSU. It reveals the nonlinear material behavior including nonlinear elasticity, viscosity and plasticity. Above −200 N approximately, the stiffness decreases and approaches zero at 0 N almost. A maximum (positive) hysteresis with a width of 14 – 16 µm occurs in the ten cycles. This width includes the measurement induced (negative) hysteresis experienced in the experiments with the steel springs. Therefore the material induced hysteresis might be 20 – 26 µm in reality. From cycle to cycle the FSU accumulates permanent deformation which leads to a racheting towards a thinner IVD. This might be due to the high fluid content of the IVD. The fluid is squeezed out during testing either through the outer surface of the anulus fibrosus or through the vertebral end platens into the Ve [Huber et al., 2007]. The fluid loss is possibly the reason for the observed highly viscous response in the creep experiment (cf. Fig. 7.6), too. Even after 2 hours no significant deceleration of the strain rate is observed. The extreme deflections in the signal occurring after 2 hours of creep (cf. Fig. 7.6) are due to a spontaneous destabilization of the controller, upon which the traverse of the materials testing machine began to oscillate. This problem occurred unpredictably several times after different test durations. Consultations with the manufacturer of the materials testing machine revealed that the most possible cause for this behavior are electro-magnetic waves due to a freezer operating nearby. The force-transducer cables are not well shielded and act like antennas.2 As the creep experiments are done in a force controlled mode, this unexpected changes in signal disturb the control loop of the testing machine and the system begins to oscillate. In a displacement controlled mode fluctuations in the force signal do not disturb the controller and have therefore no influence on the control stability. Figure 7.7 presents the outcome of 1 hour of relaxation: it is almost sufficient for the sample to attain equilibrium. A comparison with a calibrated exponential function shows that the force reached almost 10 % of the asymptotic force value that can be expected. Within 1 hour of relaxation the force drops by 53 % of the peak force (−900 N). Relaxation experiments reach equilibrium much faster than creep experiments. Therefore they could be used to determine the long term elastic properties of the FSU by i 2

Ledet et al. [2000]; Nachemson [1966a]; Wilke et al. [2001] Newer machine versions are better shielded and this type of phenomenon should no longer occur.

7.1. PRELIMINARIES

force [kN]

−0.2 −0.4 4.25 – 7.08 N −0.6

6 – 10 µm

lo un

−0.8 −1

−1.2

−1

−0.8

−0.6

−0.4

displacement [mm]

loa

−0.2

Figure 7.4: Result from a cyclic experiment with the spring molded in PMMA. By linear interpolation a stiffness of 708 mm/N is measured. The detail plot shows the sign of the hysteresis (arrows) and its magnitude.

0.2 0 softening

force [kN]

−0.2 −0.4

hysteresis

−0.6 −0.8 −1 −0.4

racheting −0.35

−0.3

−0.25

−0.2

−0.15

displacement [mm]

−0.1

−0.05

0.0

Figure 7.5: Result from a cyclic experiment with a FSU. The last cycle is highlighted with a dark line. Clearly visible: the softening at lower compressive loads, the hysteresis and the racheting. 81

displacement [mm]

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

0 −0.2 −0.4 −0.6 −0.8 0

4 time [s]

7 ·103

0 force [kN]

−0.5

8 ·103

displacement [mm]

Figure 7.6: Creep experiment with a duration of more than 2h; no equilibrium was attained.

0 −0.1 −0.2 −0.3

0.5

1.5

2 time [s]

2.5

2 time [s]

2.5

·103

0 force [kN]

3.5

−0.2 −0.4 −0.6 −0.8 −1

0.5

1.5

3.5 ·103

Figure 7.7: Relaxation experiment lasting for 1h. The force dropped by 53 % of the peak force, equilibrium is almost attained. indicates the asymptotic force value.

7.1. PRELIMINARIES

fitting series of equilibrium states at different deformations. After this first test runs it got obvious that the displacement measurement had to be improved. Besides this result from the obtained test data, the experiences with the sample preparation and handling led to some changes.

7.1.3

Improvements

Mountings A critical revision of the mountings design after first test runs led to improvements that culminated in the last and final version of the mountings (version four, cf. Fig. 7.1b). The main idea of a hollow cylinder remained but much else changed. The aluminum from the first design was replaced by stainless steel. The interface size was found to be insufficient for the high forces that were to be expected (in the order of 1 kN). Therefore it was redesigned such that an 18 mm pin can be used to lock the interface. With this interface the mountings can be attached directly to the 20 kN force transducer (GTM-K 20kN, cf. App. D). The cavity in which the Ve are ingrained was made less deep such that less PMMA is used. A steel ring was introduced that is bolt from top onto the hollow cylinder. The aim is to pretension the PMMA in order to avoid any backlash in the mountings. Displacement measurement The negative hysteresis observed when measuring steel springs cyclically highlighted the need for an improvement. A first step was the full understanding of the sources for this phenomenon. A negative hysteresis means that more energy is freed during unloading than was put into the system during loading. Therefore it is certainly not an attribute of the material response. An alternative description of hysteresis, besides its characterization by its width, brings some insight. Hysteresis is nothing else than a temporal shift of the displacement and force signal (cf. Fig. 7.8). If the displacement lags behind the force a positive hysteresis results, if the contrary is true the result is a negative hysteresis. In the above mentioned experiments the time by which the force signal lags behind is approximately 0.15s. Two possible sources for this behavior were identified and one was shown to be the true one. The first mechanism could be that the force signal at a given time point is correlated with the displacement signal of an other time point. This might happen if data acquisition is not made correctly, something that can be excluded (in the observed order of magnitude) for a commercial materials testing machine. The second mechanism is some backlash in the drive mechanism for which is not accounted in the displacement acquisition. And this is the case with the used materials testing machine which uses a spindle to drive the traverse. The displacement of the traverse on the other hand is not measured locally but is computed from the turns of the motor driving the gear. So, what happens is the following: the motor begins 83

d | f [arb.]

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

f [arb.]

t [arb.]

d [arb.] (a)

d [arb.] (b)

d [arb.] (c)

Figure 7.8: Visualization of hysteresis (t: time; d: displacement [ ]; f : force [ ]). (a): no hysteresis, force and displacement are in phase; (b): positive hysteresis, the force advances the displacement; (c): negative hysteresis, the force lags behind the displacement.

to turn (and therefore some displacement is calculated for the traverse) but, due to backlash, the traverse does not move until the backlash is not overcome. Since no real displacement happens the first instant the force does not change either, therefore the force lags behind the (computed) displacement. This mechanism holds especially for load reversal points. The conclusion is to use alternative displacement transducers that measure the displacements as close to the sample as possible. The use of an external extensometer (Multisens, cf. App. D) brought the desired improvements. The extensometer was placed such that it measures the change in the free length between the two steel mountings.

Protocol for the sample preparation One major change in the sample preparation raised from the first experiments with FSUs. It was found that 24 hours of defrosting in a refrigerator is enough to unfreeze the samples. The overnight stay mounted in the materials testing machine at room temperature would remove any residual coldness inside. In contrast to what is done in other studies [e.g. Goel et al., 1985], here the samples are made ready as much as possible (dissection, casting of one side) before freeze storage. 84

7.1. PRELIMINARIES

Though the final protocol for the sample preparation was the following: - acquisition of fresh ovine lumbar spines, - dissection of the FSUs (including LV-LVI, removal of surrounding muscular and ligamentous tissues), - detach the posterior elements, - ingrain LVI in the lower mounting, taking care that the transversal plane of the IVD is as horizontal as possible, - ejection of the sample after the polymer hardened, - wrap the sample with saline soaked gauze, - seal the sample in a plastic bag and freeze storage until the experiments start, - defrost the sample in a refrigerator 24 hours prior measuring, - the evening before the experiments: -

insert the sample in the upper mounting, fix the mounting at the top interface of the materials testing machine, ingrain of the second V, wrapping the sample loosely with saline soaked gauze, covering all with cling film, - wait until the polymer hardened, - position the traverse such that zero force results, - start of the experimental protocol. During all this steps care was taken that the sample does not dehydrate. The duration for which the sample were in direct contact with ambient air was held as short as possible and samples were repeatedly sprayed with saline. Shortly after the mounting of the sample in the materials testing machine it was loosely wrapped with saline soaked gauze and covered with cling film. Checks after the experiments found the gauze still wet, therefore any dehydration of the sample seems excludable.

7.1.4

Extension for bending experiments

Principle The experiments should be realized with the same uniaxial materials testing machine and the same mountings as the axial loading experiments. It was therefore necessary to extend the setup such that moments can be induced. To this aim two lever arms were used, see Figure 7.9. The lower mounting is directly attached to the traverse of the testing machine. A lever arm is attached to the upper mounting such that a moment can be applied to the disc by applying a force on it. The vertical force is applied via a 50 N load transducer (cf. App. D) attached to a second lever arm fixed to the upper machine interface. The front (defined by the ventral direction of the FSU) of the lower lever arm is designed such that only a point contact exists between the lever arm (which has a convex surface) and the force transmitting cylindrical pin (cf. detail view in Fig. 7.9); besides some friction only vertical forces are transmitted. 85

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

With this setup two quantities are measured: the force applied to the lower lever arm and by means of an additional extensometer the displacement of a point on the rear end of the lever arm. For well definiteness of the contact points of the extensometer with the lever arms flats are machined on the rear ends of the latter. It is clear that with such a setup no pure moment can be applied to the disc. The FSU is preloaded by the weight of the upper mounting and the lower lever arm (approx. 2.5 kg), and additionally, it must support not only the induced moment but the vertical load as well (up to 20 N). Taking the FSU stiffness in the order of 900 mm/N into account, the vertical load of maximal 45 N is neglected. Data evaluation is done as follows (Fig. 7.10 shows schematically all necessary geometric quantities): It is assumed that all parts behave as rigid bodies except for the IVD which can bend and change length and the force transducer which has a known compliance cs (cf. App. D). As the materials testing machine moves the traverse by δt the lower lever arm will be translated, too. As this movement is restrained by the force transducer, which complies with δs = F cs , where F is the exerted force, the lower lever arm will rotate around the contact point of the force transducer, though an angle α forms with the horizontal. This rotation induces a vertical movement of the contact point of the extensometer on the rear end of the lever arm such that the extensometer measures a change in length δe . Therefore the relative vertical positions of three points are known: 1) a point on the rear end of the lever arm (measured by the extensometer), 2) the point of load incidence of the force transducer (known via its known compliance and the measured force) and 3) a point on the traverse (measured by the built-in displacement transducer). With the above mentioned assumptions about the possible displacements the rotation angle α can be calculated: α = arctan[

δe − F cs δe − δs ] = arctan[ ]. le0 + ls0 le0 + ls0

(7.1)

Experiments with a steel spring First experiments were done with the steel spring ingrained in PMMA described in Section 7.1.2. Figure 7.11 shows one example result that was obtained for cycles between ±50 N, resulting in ±10 N m bending moment. A first observation is the small region around zero angle with constant moment. In terms of displacements this are roughly 0.2 mm, which can be attributed to the backlash present between the lower and upper sliding surfaces and the cylinder pin with which the force to the lever arm is applied (see detail view in Fig. 7.9). As this “gap” is small and it takes not much time to traverse it might be appropriate to neglect it or to correct the experimental data for it. The second remark is the loop present on top of the loading curve in flexion. In some experiments with the spring it was there, in others not. It is assumed that the existence of this loop is a control issue and depends on the speed with which the traverse is driven and on the stiffness of the sample. Eventually, it was not present in the experiments with the FSUs. 86

7.1. PRELIMINARIES

machine interface upper lever arm flats for extensometer lower lever arm

force transducer @ @ @ @

upper mounting

@ @ @ @

FSU lower mounting

B B B

sliding surface

flats for attachment of a fixation

18mm pin machine interface detail view Figure 7.9: Drawing of the setup used for the bending experiments. The lower machine interface is mobile: an upward translation results in flexion, a downward translation in extension of the disc.

α ls

δe

δs le0 h0

ls0 h

δt

Figure 7.10: Scheme defining all variables that are used for the interpretation of the bending experiments. The traverse is displaced by δt . The sample is, by action of the lever arm, compressed by h0 − h and bent by α. Due to this bending the extensometer changes its length by δe . The compliance of the force transducer cs is taken into account by δs = F cs . 87

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

12 10 8

loop

d loa

6 moment [N m]

4 2

backlash

0 −2 −4 −6

loa

−8 −10 −12

extension −2.5

−2

−1.5

−1

flexion −0.5

0 angle

0.5

1.5

2.5

[◦ ]

Figure 7.11: Results of a bending experiments with a steel spring. Remarkable are i) the horizontal segment around zero angle, ii) the loop at the top end of flexion, iii) the presence of three different stiffness in extension and iv) the saw tooth pattern at higher moments.

Third, in extension (and in other experiments as well in flexion) three distinct stiffness exist and therefore a hysteresis is present. As the steel spring is elastic the hysteresis must have its origin somewhere in the test setup. The sharpness of the transitions from one stiffness to the next indicates some kind of backlash problem. Up to date no solution for this problem was found. For this reason it was decided not to use the full cyclic data but just the loading curves. Especially at the loading curves a saw tooth pattern can be observed. This was found to be the effect of friction that acts between the sliding surfaces and the cylinder pin with which the force is transmitted.

7.2

Protocols and results

The ambition of the experiments with the FSUs was to obtain a reliable data basis for a comparison with predictions obtained by FE simulations. For this reason FSUs were characterized in two different types of motions: 1) axial compression/tension and 2) flexion/extension. Both experiments include cycles with different deformation rates, repetitions of the cyclic experiments give insights into the repeatability of the measurements (plasticity). For the final experiments, samples were prepared as described in Section 7.1.2. 88

7.2. PROTOCOLS AND RESULTS

7.2.1

Axial loading

Experimental protocol The detailed experimental protocol is given in Table 7.1, durations of each run were estimated from a first experiment, the exact durations depend on the sample stiffness. The setup including the sample mounting were prepared the evening prior measuring, as explained in Section 7.1.2. This was done to allow the FSU to relax overnight at its initial length (#0) and at room temperature. The initial length of the FSU was determined by zeroing the load manually. In order to prevent dehydration the samples were covered with loosely wrapped saline soaked gauze whenever possible and especially during the experiments. The first experiments were done in a compressive regime (which is the state in which the IVD is mostly loaded [Koeller et al., 1984]). Cycles with three different displacement rates (0.05, 0.5 and 5 mm/min) at forces between −900 N and −100 N were run (#1, #3 and #5). The limit forces were chosen according to findings on humans and baboons who showed that peak forces up to 1.5 times the bodyweight are supported by the FSU without any damages.i For sheep this suggests peak forces down to −900 N. Due to the small speeds at the lowest rate only six cycles (otherwise ten) were run. Between two cyclic runs the FSU was allowed to relax for 1 hour at its initial length (#2, #4 and #6). After completion this protocol was rerun a second time (#7 to #12). After completion of the experiment in the compressive regime, six cycles were run that cover changes from the compressive to the tensile regime (and vice versa, #13). Experiments were run with a displacement rate of 0.5 mm/min between −450 and 450 N. The experimental protocol ends with six cycles in the tensile regime from 100 – 900 N (#15) and a final relaxation phase. The total duration (including the relaxation over night) was estimated to be about 27 hours; the effective durations were even shorter (about 23 hours). This duration is comparable to other experiments that lasted for 28 hours [Smeathers, 1984]. In total six specimens were characterized with the above protocol. After experiments specimens were removed from the materials testing machine and freeze stored for later examinations. Results The absolute values of the measured quantities vary considerably from one sample to the next (not shown). This might be explained by the variability of the FSU sizes. In fact, the variability of relative measures is much smaller. If, for example, the maximum displacements measured in the first run of the cycles with 0.05 mm/min (#1) are considered, the standard deviation is rather high: 0.71±0.3 mm (average ±standard deviation). But if the changes of the maximum displacements from the first to the second run of the cycles (#1 and #7, respectively) are considered, the variability is quite small: 1.08±0.02 –. i

Ledet et al. [2000]; Nachemson [1966a]; Wilke et al. [2001]

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

Table 7.1: Protocol that was used for the uniaxial experiments with the FSUs. Durations for each run were estimated from a first measurement. #

description

duration [min]

relaxation over night 0

900’∗

relaxation at the initial length

2nd run

1st run

experiments in compression 100’∗

6 cycles with 0.05 mm/min between −900 N and −100 N

relaxation at the initial length

60’

10 cycles with 0.5 mm/min between −900 N and −100 N

20’∗

relaxation at the initial length

60’

10 cycles with 5 mm/min between −900 N and −100 N

relaxation at the initial length

6 cycles with 0.05 mm/min between −900 N and −100 N

relaxation at the initial length

60’

10 cycles with 0.5 mm/min between −900 N and −100 N

20’∗

relaxation at the initial length

60’

10 cycles with 5 mm/min between −900 N and −100 N

relaxation at the initial length

2’∗ 60’

100’∗

2’∗ 60’

experiments in compression and tension 13

6 cycles with 0.5 mm/min between −450 N and 450 N

30’∗

relaxation at the initial length

60’

experiments in tension 15

6 cycles with 0.5 mm/min between 100 N and 900 N

30’∗

relaxation at the initial length

15’ total duration: 1639’∗ 27.3h∗

∗

estimated durations

7.2. PROTOCOLS AND RESULTS

Figure 7.12 gives an overview of the measured behavior of one representative FSU. The total measurement duration was about 23 hours. Already here, a remarkable repeatability in the force signal can be observed. Even after the latest cyclic experiments there is no considerable deviation from zero force once the sample is set back to its initial length. More details are shown in Figure 7.13 that shows the hysteresis curves for all the eight cyclic experiments that were done. Cyclic experiments with same displacement rate and same force range are shown together such that the repeatability can be compared. Looking at Figures 7.13a – c reveals in fact a remarkable repeatability in the initial force that is almost zero for all six measurements. At all cyclic experiments hysteresis and ratcheting can be observed. Both phenomenon become smaller with increasing displacement rates. Another thing worthy to note is that with increasing displacement rates the first and second runs become more and more distinct. At a rate of 0.05 mm/min “loops” can be observed at the upper reversal point: at the start of a next loading the material behaves stiffer than at the previous unloading such that the loading path crosses the unloading path. At higher rates this “loops” are present as well, the crossing point moves to higher (absolute) force values. While the results in compression are nice and smooth this is not absolutely true for experiments in the tensile regime (cf. Figs. 7.13d – e). It seems there is still some backlash combined with friction in the test setup: first loading curves in the tensile regime show a discontinuity in the force-displacement characteristics. Once the first loading in tension is finished no other discontinuities happen; especially for the experiments cycling between compression and tension (cf. Fig. 7.13d) this is astonishing. Besides this shortcomings in the tensile experiments some statements can be made. The first is that the FSU is softer in tension than in compression (cf. Fig. 7.13d): at 400 N in extension the stiffness is about 1840 N/mm while for the same compressive force a stiffness of about 3430 N/mm is measured. A further observation is that the FSU has a non vanishing stiffness at the neutral position which is about 450 N/mm.

displacement [mm]

−1

−0.8

−0.6

−0.4

−0.2

0.2

0.4

0.6

0.8

−0.8

−0.6

−0.4

−0.2

0.2

0.4

0.6

time [s]

#1 #2

#6 #7 #8

#10 #12

#14

·10

Figure 7.12: Visualization of the results of one representative experiment. The total duration was approximately 23 hours.

force [kN]

0.8

#16

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

7.2. PROTOCOLS AND RESULTS

force [kN]

−0.2

−0.4

−0.6

−0.8

−1

−0.7

−0.6

−0.5

−0.4

−0.3

displacement [mm]

−0.2

−0.1

(a) cycles with 0.05 mm/min between −900 N and −100 N (#1 and #7) 0

force [kN]

−0.2

−0.4

−0.6

−0.8

−1

−0.7

−0.6

−0.5

−0.4

−0.3

displacement [mm]

−0.2

(b) cycles with 0.5 mm/min between −900 N and −100 N (#3 and #9) Figure 7.13: Details of the cycles: first run ( pages 94 and 95.

) and second run (

). Continued on

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

force [kN]

−0.2

−0.4

−0.6

−0.8

−1

−0.7

−0.6

−0.5

−0.4

−0.3

displacement [mm]

−0.2

−0.1

force [kN]

0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5

−0.4

−0.3

−0.2

−0.1

0.1

0.2

0.3

0.4

displacement [mm]

(d) cycles with 0.5 mm/min between −450 N and 450 N (#13) Figure 7.13: Details of the cycles: first run ( page 93 and continued on page 95. 94

) and second run (

). Continued from

7.2. PROTOCOLS AND RESULTS

0.8

force [kN]

0.6

0.4

0.2

0 0

0.1

0.2

0.3 0.4 displacement [mm]

0.5

0.6

0.7

(e) cycles with 0.5 mm/min between 100 N and 900 N (#15) Figure 7.13: Details of the cycles. Continued from pages 93 and 94.

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

7.2.2

Bending

Experimental protocol FSUs were ingrained into the lower mounting such that the IVD was horizontal as perfect as possible and such that its ventral direction pointed in the direction of the lever arm (such that a downward directed vertical force on the lever arm results in flexion). To ingrain the second V the other mounting was fixed rigidly to the first one by use of a custom made fixation plate attached to the flats machined on the sides of the mountings (cf. Fig. 7.9). The FSUs were ingrained such that their IVDs were vertically centered between the two mountings. Whenever manipulations it allowed, samples were covered with loosely wrapped saline soaked gauze which aimed at preventing dehydration. Also during the experiments samples were always covered. The protocol consisted of four cyclic experiments with two different displacement rates (of the traverse) and intermediate relaxation phases of 15 minutes (cf. Tab. 7.2). The experiment started with the fixation of the sample in the testing machine and removal of the fixation plate such that the FSU was free to deform. The position of the traverse was adapted such that the force transmitting cylinder pin had no contact with the lower and upper sliding surfaces of the lower lever arm. Once the setup was ready the sample relaxed for 15 minutes (#0). Afterwards six cycles between ±20 N (corresponding to ±3 N m) were driven with 17 mm/min displacement rate (#1 and #5) followed by 15 min relaxation. Ten cycles were driven between the same force limits but with a ten times higher displacement rate (#3 and #7). This full sequence was repeated once more such in total four cyclic experiments were conducted on the same sample. Results Figure 7.14 visualizes the results of the one FSU characterized in flexion/extension. The visible gaps in the data result from the time needed to switch between the experimental protocols for relaxation and cyclic experiments and, after the cyclic experiments, to set the traverse to its initial position. A closer look at the angle vs. time graph reveals that the angles change almost linear with time. In fact, the time rates of the angle are about 6 ◦ /min for a displacement rate of 17 mm/min and 60 ◦ /min for 170 mm/min. Figure 7.15 shows the results of the cyclic experiments in more details, with moments reported against the bending angle α (calculated according to Eq. (7.1)). As in the results of the experiments with the steel spring (cf. Fig. 7.11) three different types of stiffness regimes can be seen: the two nonlinear regimes for loading and unloading and a third, linear, between loading and unloading. Therefore, as already mentioned above, only the loading paths will be used for a comparison with the FE simulations (cf. Chap. 9). Obvious is the high repeatability of the measurements. The two cyclic experiment with the higher displacement rates are congruent. The first cycles at lower rates are a little stiffer than the second ones, that are congruent with the results with the higher rates as well. 96

angle [◦ ]

−4

−3

−2

−1

−8

−6

−4

−2

time [s]

6 ·103

Figure 7.14: Visualization of the results of one representative experimental series. The total duration was approximately 2 hours.

moment [N m]

7.2. PROTOCOLS AND RESULTS

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

moment [N m]

2 1 0 −1 −2 −3

extension −8

−6

flexion −4

−2

angle

[◦ ]

(a) cycles with 17 mm/min between ±3 N m (nominal). 3

moment [N m]

2 1 0 −1 −2 −3

extension −8

−6

flexion −4

−2

angle

[◦ ]

(b) cycles with 170 mm/min between ±3 N m (nominal). Figure 7.15: Details of the cycles: first run (

) and second run (

7.2. PROTOCOLS AND RESULTS Table 7.2: Protocol that was used for the bending experiments with the FSUs. #

description

duration [min]

initial relaxation 0

relaxation at the initial length

15’

2nd run

1st run

bending experiments 1

6 cycles with 17 mm/min between −20 N and 20 N

22’∗

relaxation at the initial length

15’

6 cycles with 170 mm/min between −20 N and 20 N

relaxation at the initial length

15’

6 cycles with 17 mm/min between −20 N and 20 N

22’∗

relaxation at the initial length

15’

6 cycles with 170 mm/min between −20 N and 20 N

5’∗

total duration: 114’∗

∗

estimated durations

The graphs in Figure 7.15 indicate an asymmetric behavior of the disc with higher stiffness in flexion than in extension. Visual inspection at the start of the first relaxation phase showed that the lever arm is not in a perfect horizontal position once the fixation is removed from the two mountings but the disc is already in flexion. This can be understood by the fact that the lever arms are not well balanced and have more weight on the ventral side, though the FSU will flex. Digital pictures were taken at a relaxation phase from which the initial flexion angle was computed to be approximately 1.7◦ . As the angle and force at the position at which the FSU relaxed did not change much over time, the determined angle can be considered as initial angle for all cyclic experiments. The vertical and horizontal lines in Figure 7.15 indicate the corrections for the angle and the moment that have to be done in order to get the unloaded (except for the axial load) position for the FSU. The initial moment of 0.6 N m acting on the FSU was calculated with the computer aided design (CAD) drawings of the setup. It follows that the FSU was loaded more in flexion (approx. 3.6 N m) than in extension (approx. 2.4 N m). With respect to this two lines the response is almost point symmetric. 99

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

7.3

Discussion

One of the most important measures that resulted from the preliminary experiments is the abandonment of the built-in deformation transducer for the axial loading experiments. It was shown that the reliance on this transducer is deceiving since it induces non-physical phenomenons. As the displacement is not measured where it is applied to the sample it is measured at the gear. Backlash in the mechanical construction of the materials testing machine (which is almost inevitable) leads to negative hysteresis. In immediate experiments with a viscoelastic material this phenomenon would probably not have been recognized because the (positive) constitutive hysteresis would outnumber the negative, measurement induced one. The use of external transducers that measure the change in free length between the two mountings solves the problem. The same mounting principle as for the longitudinal loading was to be used for bending experiments as well. Though the setup was extended by two lever arms that are used to induce a moment in the FSU. A phenomenon that is present in bending experiments with the steel spring and FSUs, is the existence of three distinct stiffness regimes. Unfortunately no solution for this phenomenon was found and therefore only results for the loading curves are used for a comparison with FE simulations (cf. Chap. 9). The fact the the samples were freeze stored prior measuring might have an influence on the results of the experiments. This subject is discussed controversially in literature, nevertheless many studies use freeze storage for sample conservation. Smeathers and Joanes [1988] report a reduced compressive stiffness of FSUs due to freezing that might relate with a change of the osmotic pressure of the tissue. As the osmotic pressure is equivalent to a mechanical stress [Schneiderman et al., 1986], its change alternates the load bearing capabilities of the IVD. In fact, Bass et al. [1997] tested the effect of freezing on the long term creep response of FSUs and found permanent and not recoverable changes for the frozen samples: a decreased osmotic swelling pressure of the nucleus pulposus (−25 %) and a higher permeability. Nonetheless it was found that the elastic behavior was not significantly affected. The measured differences to previous studies are attributed to the fact that in this study the long term response was characterized and not responses for short times. Other studies report no changes of the mechanical responses in flexion, lateral bending and axial torsion for vertebral segments that were frozen for up to 7 months [Gleizes et al., 1998; Panjabi et al., 1985]. The protocols for axial loading and bending experiments were designed to highlight viscous properties of the IVD. The fact that the experiments were run at room temperature is expected to affect the viscoelastic behavior of the IVDs only quantitatively but not qualitatively [Koeller et al., 1986]. Cyclic experiments at different speeds allow an insight into rate dependency and repetitions of the load cycles give information about repeatability and plastic effects. Experimental protocols, and especially for the axial loading case, are long lasting; care 100

7.3. DISCUSSION

was taken to prevent dehydration. Samples were wrapped in saline soaked gauze covered by cling film. Even after 23 hours of experimentation the gauze was found humid and therefore dehydration of the disc is excluded. Another aspect related to the hydration of the samples should be mentioned. Under pressure fluid is squeezed out of the IVD and enters, via the endplates, the vertebral vascular and trabecular space [Ayotte et al., 2000; Beutler et al., 2002]. Casting the Ve, as in the present study with PMMA, makes it impossible for the fluid to escape via the vertebral surfaces. Huber et al. [2007] studied the influence of the prohibited fluid flow and suggest to use a fluid reservoir that is connected to the cancellus vertebral bone part such that fluid can leave and reenter the Ve. They found that with such a reservoir peak forces in creep experiments are more consistent than in the case of sprayed samples; the long-term response is less affected. Results of the axial loading experiments show a remarkable repeatability of the force at the initial position even after the last experiment. Also this is taken as evidence that no structural changes or dehydration took place. Other phenomenon that can be observed are ratcheting and a distinct hysteresis. At all displacement rates “loops” at the upper reversal points of the force-displacement diagrams are present. This loops arise from a stiffer characteristics in loading than in unloading and can also be observed in simple spring-dashpot models (e.g. a spring in parallel with a spring and dashpot in series). Bending experiments with FSUs are influenced by the unbalance of the test equipment. It was found by visual inspection and proved by some calculations that the FSUs are not loaded from their neutral position on. Instead, due to the unbalance in weight of the lever arm, the disc is prestressed by a 0.6 N m bending moment in flexion. The FSUs are not cycled between ±3 N m but between 3.6 N m in flexion and 2.4 N m in extension. Next experiments could account for this by defining the neutral point differently or some additional weight on the rear end of the lower lever arm could balance the setup. Finally, no pure moment is induced by the setup as it was used, but a moment superimposed by some axial load due to the weight of the setup and the vertical force component applied to the lever arm. The total axial load (25 – 45 N) is small compared with the stiffness of the IVD (in the order of 900 N/mm) and might be neglected. In fact, pure bending is seldom present in the spine, most often bending moments and axial forces superimpose (e.g. muscle forces, body weight) [Costi et al., 2002]. Pure bending moments might be realized with a setup proposed by Panjabi et al. [1994]. With the final experimental protocols six FSUs were characterized under axial loading and one in bending. This low number does not allow to make any valid statistical statements and this was not the aim of this study. One data set from the axial loading case that is felt to be “representative” and the one result from the bending load case are taken as benchmarks for FE simulations (cf. Chap. 9). The mountings and additional setup that was used for the axial loading and bending experiments are a custom solution. Validation experiments with steel springs revealed 101

CHAPTER 7. EXPERIMENTS WITH FUNCTIONAL SPINAL UNITS

some of its major shortcomings, which were fixed consequently. Still, it is not perfect and with the increasing number of experiments the knowledge about further improvements grows. Despite the remaining uncertainties and the low number of experiments, the results obtained are thought to be relevant for the modeling process.

102

Chapter

Calibration of Elasto-Viscoplastic Equations Elasto-viscoplastic constitutive equations are calibrated with results from uniaxial experiments on strips of anulus fibrosus tissue. A stepwise procedure is applied where first only the elastic part is calibrated, succeeded by the rate-independent parts and finally the rate-dependent parts. As no suitable measurement data for the first two steps was available, reference curves were extracted best possible from the experimental data. Results partially negative the arguments for the stepwise procedure: parameter sets from preceding steps needed to be recalibrated in next steps in order to get acceptable outcomes. The parameter sets found for the elastic and rate independent part are quite promising while for the model including also the rate dependent parts the viscous behavior is not reproduced to full satisfaction.

The results of the uniaxial experiments on the anulus fibrosus (AF) strips (cf. Sec. 6.5) are used as reference data for a parameter calibration of an elasto-viscoplastic model. This model, proposed in 2002 by Rubin and Bodner, is briefly presented in Section 3.3.3. In order to structure the calibration process and to analyze the importance of the different aspects of the model, which are elasticity, viscoplasticity and rate dependance, three calibrations are done, each one adding some more degrees of freedom to the model. The first calibration uses only the elastic part of the model, the second includes rate-independent viscoplasticity and the third adds rate-dependence. The calibration is done using the MATLAB (cf. App. D) optimization algorithm patternsearch in a routine developed and described by Hollenstein [2011]. In the present work the finite element (FE) program ABAQUS (cf. App. D) was used for the evaluation of a given parameter set, i.e. given material parameters. As the cooperation of MATLAB and ABAQUS is sometimes tricky some of the adjustments made to the optimization environment are explained in Appendix C.1. In the following the FE model, the reference data sets and the calibrations with results are presented and discussed. 103

CHAPTER 8. CALIBRATION OF ELASTO-VISCOPLASTIC EQUATIONS

8.1

The finite element model

As the aspect ratio of the measured samples is high (cf. Sec. 6.5) an uniaxial stress state was assumed for the calibration. Therefore only one brick element (type C3D8H) is used for the calibration (cf. Fig. 8.1). The element was defined to have the same length L0 as the sample such that it can be loaded with the absolute displacements known from the experiments (L0 = 17.1992 [mm]). In the experiments it was found that the samples get longer after one cycle. I.e. when the initial gauge length, at which the samples are stretched before the cycle, is regained, the samples are slack. The measured force, indeed, vanishes before the initial gauge length is regained. But as the samples have a low bending stiffness they cannot withstand any compressive load and, therefore, no compressive forces were measured. How to simulate this with FE? In a simulation of a perfect geometry, as present with a single prismatic element, bending cannot be simulated and another solution has to be found. An option would be to change the constitutive model such that it does not produce compressive forces. As the user subroutine that defines the constitutive behavior is given [Papes and Mazza, 2008] also this is no solution. The remedy found is to connect one end of the FE model with the (virtual) ground by an axial slider which has an attributed “stop”-behavior. This means that the slider disables movements in one direction (the tensile direction) but allows for displacements in the other (compressive) direction (cf. Fig. 8.2 for a schematic sketch). This axial slider was applied to only one point P1 of the element (cf. Fig. 8.1), points Pi , i = 2, . . . , 4 could not be connected to the ground since this would prohibit them from their movements induced by the lateral contraction. In order to keep them in-plane with P1 they are connected to the latter by slide-planes that allow for in-plane movements (lateral contraction) but not for out-of-plane movements. With this kinematic model it is possible to simulate adequately the slacking of the material and therefore to reproduce the experiments. In reality two families of fibers are present in the AF. Both are alternatingly present in the different layers that make up the structure of the AF. The fibers include angles of ±30◦ with the transversal body plane (see as well Sec. 2.2). For the simulations two assumptions were made on the distribution of the fibers:

1. the strips were cut in circumferential direction 2. each strip contains several layers of the AF such that both fiber families are equally present and therefore the fiber distribution might be hom*ogenized.

For the simulations this implies that, other than in reality, both fiber families are present throughout the model. The fibers are distributed symmetrically around the y-z-plane and make an angle ϕ = ±30◦ with it (cf. Fig. 8.1). 104

8.2. CALIBRATIONS y L0 P6 P2 P7

P3 2ϕ P5

P1 P8 x

Figure 8.1: Visualization of the FE model used for the calibration and its boundary conditions. The prescribed displacement is applied to points Pi , i = 5, . . . , 8. Point P1 is fixed to the (virtual) ground by an axial slider (not shown). To keep points Pi , i = 2, . . . , 4 in-plane with point P1 they are connected to the latter by slide-planes. The two fiber families are shown exemplarily, they include an angle of 2ϕ = 60◦ .

8.2

Calibrations

One experiment with the AF strips that is felt to be representative is taken as reference data set for the calibration. This is the set already presented in Figure 6.9a in terms of a force vs. stretch diagram. For this plot, as explained in the corresponding Section 6.5, the initial length of the sample was (re-)determined at the beginning of every cycle. For the calibration absolute changes in the clamp-to-clamp distance (the displacement) are used. Such a force vs. displacement diagram is shown in Figure 8.3a. Two text files were prepared for the calibration. The first contains a list with times and corresponding reference displacements. This file is read in by ABAQUS and determines the displacements of the nodes Pi , i = 5, . . . , 8 (cf. Fig. 8.1). After the simulation the results (forces) are compared with the reference forces read from the second text file listing times, displacements and forces and an error value f is computed. This value is defined as the sum of the absolute difference between the measured (reference) and the computed forces at all time increments: f=

N X i=1

abs[Firef − Fisim ].

(8.1)

Here F ref and F sim are lists (of same length) containing the reference and the simulated force values at N time points ti , respectively. The optimization algorithm decides, based on f and its history, which parameters to change in order to further minimize the difference between reference and simulations. Elasto-viscoplastic equations [Rubin and Bodner, 2002] were calibrated; a short description of the equations is given in Section 3.3.3. O. Papes implemented this model into a FORTRAN subroutine (UMAT) such that it can be used with ABAQUS [Papes and Mazza, 2008]. 105

CHAPTER 8. CALIBRATION OF ELASTO-VISCOPLASTIC EQUATIONS

∆

(a)

∆

L∆ L∗0

(b)

(c)

Figure 8.2: Schematic of the model that was chosen for the calibration of the uniaxial experiments on AF strips: initial configuration (a), deformed configuration (b) and configuration after deformation (c). A Kelvin-Voigt model is shown as illustration, for the calibration a different model was used.

In total the model has fourteen parameters of which six describe the elastic response of the material (cf. Tab. 8.1). Though a stepwise calibration procedure was designed in order to improve accuracy and to study the influence of the different aspects of the model. In a first step only the elastic part of the model was fitted to the elastic response of the material. As the elastic response of the AF strips was not measured explicitly, it had to be approximated the best way possible. This was done taking the last cycle of the reference data, which is the most preconditioned materials response. It still inherits a large amount of viscosity, as can be seen in Figure 8.3b. To get a reference curve which describes the elastic behavior the best, only the stress-strain data points after each relaxation phase were taken. In a second step rate independent visco-plastic effects were included, adding another three parameters to the model (cf. Tab. 8.1). Now the parameters describing the elastic behavior were already calibrated, such that, in theory, only the new parameters had to be fit. In practice it was found that the four reference data points for the elastic behavior are not sufficient to fully characterize the elastic part of the model. Nevertheless the first step gave a hint what the elastic parameters are like. To fit the rate independent behavior, which as well was not measured explicitly, a composite curve was designed (cf. Fig. 8.3b). This curve follows the original curve until the first relaxation pause. The next two parts are the shifted respective loading paths after the relaxation pauses. The unloading path is defined by the shifted measured unloading path. In the third, and final step, also the rate dependency of the experimental data was considered. This was done by calibrating an additional parameter Γ1 (cf. Tab. 8.1) to the full reference data set, including the relaxation phases.

106

8.2. CALIBRATIONS

last cycle 10

force [N]

8 6 4 2 0 0

0.5

1.5

2 2.5 3 displacement [mm]

3.5

4.5

(a) full data set

last cycle points of elastic curve curve for rate indep. calibration

force [N]

8 6 4 2 0 0

0.5

1.5 2 displacement [mm]

2.5

(b) reference data set used for the calibration of the elastic and rate independent model Figure 8.3: Reference data for the calibration of the Rubin and Bodner model. Calibrations were done with an elastic reference curve (b), a reference curve accounting for rate independent viscosity (b) and the full reference data (a).

107

CHAPTER 8. CALIBRATION OF ELASTO-VISCOPLASTIC EQUATIONS

Remarks: 1. Also in the last calibration step only a reduced model with ten parameters (of which eight are free) was used. In fact, setting r1 = r5 = 1 [-] and r3 = r4 = 0 [Hz] reduces the number of free parameters, the evolution equation for the variable yield strain β (3.46a) reduces to β˙ = r2 Γβde . Analytically being simple, the numerical implementation might get into troubles, therefore in the UMAT the evolution equation (3.46a) was replaced by the new one (in Tab. 8.1 indicated by ×). 2. As the model scales linearly with µ0 this parameter was, for the calibration, set arbitrarily to µ∗0 = 2 [MPa]. The (small strain) bulk modulus of the material (µ0 m1 ) was assumed to equal the bulk modulus of water (κ∗ = 2200 [MPa]). Therefore the parameter m1 was set equal to 1100 [MPa]. 3. If only the elastic part of the model is used by setting parameters Γi , n and ri appropriately (in the Tab. the column “elastic”), the effective parameter meff 2 ˆ 1, Ψ ˆ 2 and Ψ ˆ 3 (cf. Eq. 3.39) is given by of an elastic model considering only Ψ meff 2 = m2 + 1. Remark 2 necessitates a change in the routine that computes the error value for the optimization algorithm. As the scaling constant µ0 is defined arbitrarily, absolute force values are meaningless and therefore both, reference and computed forces, were normalized with their maximum values. The error value was calculated based on these normalized data sets. The effective values of µ0 and m1 were determined after the optimization is terminated by µ0 =

8.3

µ∗0

max[F ref ] , max[F sim ]

κ∗ max[F sim ] m1 = ∗ . µ0 max[F ref ]

(8.2a,b)

Results

Elastic model First calibrations were done on the quasi-elastic response determined from the last cycle of the reference data. As only the data points after each relaxation phase were taken only four points were available to describe the elastic behavior (cf. Fig. 8.3b). First results looked very promising, see Figure 8.4, unfortunately the solution is not unique and there exist a multitude of parameter sets {µ0 , q, m1 , m2 , m3 , m4 } that give almost identical results. The five parameter sets given in Table 8.2 cover a wide range in the parameter space and no preference is observed. Rate independent model In a next step the rate independent part of the constitutive model should be calibrated. At the beginning the elastic parameter set was fixed and only the three additional parameters Γ2 , n and r2 were free to vary. Calibrations with all five elastic parameter sets known from the first (elastic) calibration were not successful and therefore also these 108

visco-plastic

q q

µ∗0

rate indep.

rate dep.

m∗1 m∗1

m∗1

m3 m4

m4 Γ1

0 Γ2

Γ2 n

1 (×) r2

1 (×)

0 (×)

0(×) 1 (×)

DOFs

8.189·10−1 8.063 1.479·101 1.481·101 1.481·101 1.800·101 1.412·101

q [-] 8.255·103 8.377·102 4.572·102 4.571·102 4.569·102 2.758·103 2.838·103

m1 [-] 1.719·101 * 8.359·10−1 * 9.000·10−4 * * 0.000 −6 * 2.000·10 2.000 7.971·10−1

m2 [-] 1.000·101 1.000 0.000 4.860·10−4 1.000 3.726·101 3.151·101

m3 [-]

1.000·101 9.000 3.000·101 1.000 9.000 6.998 8.041

m4 [-]

0 0 0 0 0 0 9.172·10−4

Γ1 [Hz]

0 0 0 0 0 1.211·101 0.000

Γ2 [-]

ˆ 1, Ψ ˆ 2 and Ψ ˆ 3 , cf. Eq. 3.39) is used, the effective parameter equals meff = m2 + 1. If an elastic model (considering only Ψ 2

elastic

2.665·10−1 2.626 4.812 4.813 4.815 rate indep. 7.978·10−1 rate dep. 7.752·10−1

µ0 [MPa]

1 1 1 1 1 1.700·101 1.710·101

n [-]

0 0 0 0 0 2.000·10−6 2.000·10−6

r2 [-]

Table 8.2: Parameters found by calibration of the elastic, rate independent and rate dependent part of the respective uniaxial experimental data of AF strips. The omitted parameters (r1 , r3 , r4 and r5 ) are not shown.

µ∗0

elastic

full model µ0 [MPa] q [-] m1 [-] m2 [-] m3 [-] m4 [-] Γ1 [Hz] Γ2 [-] n [-] r1 [-] r2 [-] r3 [Hz] r4 [Hz] r5 [-] adm. range > 0 >0 >0 ≥0 ≥0 ≥0 ≥0 ≥0 >0 ≥0 ≥0 >0 ≥0 >0

elastic

Table 8.1: Summary of the degrees of freedom (DOFs) that were used for the three different calibrations done on the uniaxial data. µ∗0 is arbitrarily set equal 2 [MPa] and m∗1 = κ∗ /µ∗0 , where κ∗ = 2200 [MPa] is the bulk modulus of water.

8.3. RESULTS

109

CHAPTER 8. CALIBRATION OF ELASTO-VISCOPLASTIC EQUATIONS

8 7

reference points simulation

force [N]

5 4 3 2 1 0 1

1.02

1.04

1.06

1.08 1.1 stretch [-]

1.12

1.14

1.16

1.18

Figure 8.4: Result of the calibration with the quasi elastic experimental data (using the first elastic parameter set from Tab. 8.2).

four parameters were calibrated once more. The resulting parameter set, which gives a quite good agreement with the experimental data (see Fig. 8.5) is shown in Table 8.2. Only one parameter set was found that agrees as good with the experimental data. Rate dependent model The third and last stage in the process was the calibration of the rate dependent model with eight degrees of freedom. Also here first trials were made with a fixed parameter set for the elastic and rate independent part such that only the parameter Γ1 could vary. This proceeding did not give any valuable results, only when also the other seven parameters were freed better fits were obtained. The best parameter set found is given in Table 8.2 and the respective graphs are shown in Figure 8.6a. The hardening of the material is to a greater or lesser extend correctly reproduced but not the relaxation behavior.

8.4

Discussion

The FE model used to represent the AF tissue experimented in uniaxial tension consists of one single brick element under uniaxial tension. This representation is justified by the high aspect ratio of the samples (>10:1, see Sec. 6.5) and the resulting good agreement between the real experiment and a uniaxial stress situation. The sample was assumed to contain several lamellae and an equal amount of fibers of the two families. This assumption is legitimated by the fact that the samples were thick 110

8.4. DISCUSSION

8 7

reference curve simulation

force [N]

5 4 3 2 1 0 1

1.02

1.04

1.06 1.08 stretch [-]

1.1

1.12

1.14

Figure 8.5: Result of the calibration with the quasi rate independent experimental data.

compared to the thickness of the AF lamellae such that they contained multiple layers. At the level of the model it is assumed that the two fiber families can be hom*ogenized and represented by a continuous distribution of crossed fibers. In experiments samples were found to accumulate viscoplastic strain that keeps the samples longer after unloading. As the samples have a low bending stiffness they buckled and a negligible amount of compressive stress was measured. In the FE simulation this behavior is modeled by use of an axial slider with a one-sided contact attached to the sample and the ground (cf. Fig. 8.1). If the sample accumulates viscoplastic strain it elongates and, at the unloading, detaches from the dead stop and does no longer transmit any force. With this solution no modifications had to be made to the implementation of the constitutive model in FORTRAN. An elasto-viscoplastic model should be calibrated with data from uniaxial experiments with AF tissue. To this end one data set, representing the results of a repeated load protocol with intermediate relaxation times, was used to represent three characteristics of the material: 1) the elastic response, 2) a rate-independent response and 3) a rate dependent response. This procedure was chosen to simplify the calibration process by addressing specific features of the model and, second, to gain some knowledge about the importance of the specific features. A modified version of the constitutive equations with a reduced viscoplastic part was used. The evolution equation for the variable yield strain β was simplified and four 111

CHAPTER 8. CALIBRATION OF ELASTO-VISCOPLASTIC EQUATIONS

experiment simulation

force [N]

8 6 4 2 0 0

0.5

1.5

time [s]

2.5 ·103

(a) force vs. time

experiment simulation

1.25

stretch [-]

1.2 1.15 1.1 1.05 1 0

0.5

1.5 time [s]

(b) stretch vs. time

2.5 ·103

Figure 8.6: Results of the calibration with the full experimental data (a). Figure (b) shows the stretch characteristics determined in the experiment and the idealized one that was used for the calibration.

112

8.4. DISCUSSION

parameters were suppressed (see remark 1, Sec. 8.2). As long as not more experimental data is available the full complexity of the equations cannot be explored. In order to further simplify the calibration the number of degrees of freedom was reduced to the minimum possible. To this end the linear scaling factor µ0 was arbitrarily set equal 2 MPa. Therefore the absolute force values are meaningless and both, the reference and computed force, were normalized with their maximum values. The (small strain) bulk modulus κ was set equal the bulk modulus of water κ∗ = 2200 [MPa]. As κ is given by µ0 m1 this determined as well the value of m1 = κ∗ /µ0 . After calibration the two parameters were determined using Equations (8.2a,b). The elastic characteristics is represented by only four force-stretch points of the last cycle in the uniaxial experiment, that are the initial position and three points at the respective ends of each relaxation phase (see Fig. 8.3b). These last three points largely overestimate the elastic response of the material, as can be seen from the curve representing the experimental data in Figure 8.6: at the end of the relaxation phases the force is far from equilibrium. Nevertheless this three points are the best approximation of the elastic response that could be made from the available experimental data. Calibration with this reference data gave qualitatively good results, but a multitude of possible parameter sets (with almost equal qualities) was found. These parameter sets are very different and show no correlation of the material parameters (cf. Tab. 8.2). The reason for this plurality of parameter sets is attributed to the fact that the reference data is sparse and therefore no detailed knowledge about the shape of the elastic curve exists. As all experiments were interrupted for relaxation phases no continuous cyclic data is available from the experimental data. Therefore, to calibrate the rate-independent part of the model, a composite curve was constructed. It consists of the appropriately shifted four loading and unloading path segments and is thought to represent the continuous experiments. Knowledge from the previous calibrations of the elastic part were thought to be a good starting point for this next calibration. Consequently, for first runs, the elastic parameters were fixed and only the three new viscoplastic parameters were varied. Unfortunately this stepwise procedure did not result in any good parameter set and only allowing the elastic parameters to vary, too, gave better fits. For the rateindependent model only one single parameter set was found that gives good results. As can be seen in Figure 8.5 hysteresis and nonlinear stiffening is represented quite well. The last calibration was done with respect to the full experimental data including all five cycles and the relaxation phases. Also here first runs were made with fixed parameters of the rate independent model (only Γ1 was free to vary) and also this time the stepwise procedure did not give any valuable results. Consequently all eight parameters were optimized, anew. The best parameter set found represents the hardening of the material quite well but fails to represent the exponential like decay of the relaxation behavior (cf. Fig. 8.6). A similar phenomenon can be found in Nava [2007, Fig. 4.24] where the constitutive equations fail to reproduce the creep behavior of liver tissue. In summary it can be said that the enormous complexity of the used constitutive 113

CHAPTER 8. CALIBRATION OF ELASTO-VISCOPLASTIC EQUATIONS

model makes a structured calibration necessary in this case. More experimental data would be needed for a stepwise calibration, and it would be interesting to see if data sets better representing the elastic and rate-independent behavior of the material, would help to take the steps. Nevertheless the parameter sets represent to some extent the complex behavior of the AF tissue and it will be interesting to see how good they capture the axial loading and bending properties of the functional spinal units discussed in the following chapter.

114

Chapter

Simulations of Axial Loading and Bending One ovine functional spinal unit that gave representative results in axial loading experiments was digitized by means of serial computer tomography scans. The acquired data was used to create a finite element model of the functional spinal unit composed of the three main components vertebrae, nucleus pulposus and anulus fibrosus. The bony vertebrae are represented by rigid bodies and the nucleus by an incompressible fluid. The anulus fibrosus consists of “concentric” layers of an anisotropic elasto-viscoplastic material that was calibrated with data from uniaxial tensile experiments on anulus fibrosus strips. The model is validated with data from axial loading and bending experiments with ovine functional spinal units. Results indicate that the model behaves too stiff, especially at larger deformations. It is shown that the nucleus pulposus plays an important role in the activation of the fibers in the anulus fibrosus.

Insights into the mechanics of biological systems as the spine are mainly gained by appropriate experiments. But sometimes experimentation is not possible due to prohibitive ethical, technical or other reasons, or quantities that are not measurable experimentally, e.g. stress distributions inside a body, are of interest. In such cases in-silico studies might present an alternative. Modeling techniques, and in particular the finite element (FE) method, provide powerful tools for the study of the functional spinal unit (FSU) mechanics. The golden goal is to have a model with good predictive capabilities that does not need separate calibrations for each different loading pattern. Purely phenomenological models [e.g. Groth and Granata, 2008; Li et al., 1995] may well reproduce motions for which they were calibrated but fail for any others. More physically based models capture not only the systems collective mechanical behavior but, in contrast, the properties of the single components present in the FSU. Such models have far better predictive capabilities and should, whenever possible, be favored. Realistic models include all important aspects related to the problem of investigation. For FSU models Shirazi-Adl and Parnianpour [2001] list some of them: 115

CHAPTER 9. SIMULATIONS OF AXIAL LOADING AND BENDING

three-dimensional geometry, material nonhom*ogeneities, facet articulation, material and geometrical nonlinearities, time dependency, boundary conditions and modifications in the structure (e.g. failure or tissue remodeling).

Belytschko et al. [1974] were the first to present a FE model of the disc. In their two-dimensional axisymmetric model six tissues with different linear orthotropic properties are distinguished (cf. Fig. 9.1). Ever since, models gained complexity, more and more geometric details were included and the tissues are represented with more elaborate constitutive equations, accounting for time dependent effects, too. Entire motion segments [Polikeit et al., 2003a] or even multimotion segments [Eberlein et al., 2004] are modeled in three dimensions. Nowadays, FE simulations are used to study a wide variety of phenomenons in the disc/spine: disc bulge [Wang et al., 1998], fiber strains [Schmidt et al., 2010], intradiscal pressure [Shirazi-Adl, 2006], contact forces in the posterior joints [Shirazi-Adl, 1994], stress distributions [Polikeit et al., 2003a,b] and many more. FE models can provide new insights in the mechanics of the intervertebral disc (IVD), but the quality and accuracy of results obtained from such models depend on the underlying assumptions. Model validations are therefore a necessary condition before results can be trusted or even model based predictions are made. In the present work, a component based model of a FSU is tested for its (predictive) qualities. One ovine FSU was digitized by means of serial computer tomography (CT) scans and a three-dimensional FE model was created. In this model vertebrae (Ve), anulus fibrosus (AF) and nucleus pulposus (NP) are distinguished and obey individual constitutive equations. The vertebral endplates are not represented as their location and low moduli makes their mechanical effect small [Shirazi-Adl et al., 1986]. The high stiffness of the Ve, compared to the disc tissue, let to their representation z vertebral core bony endplate PP

PP PP PP

outer shell

cartilagineous endplate nucleus pulposus

anulus fibrosus horizontal axis of symmetry

Figure 9.1: Illustration of the axisymmetric model of Belytschko et al. [1974] with the six regions that were differentiated. [after Belytschko et al., 1974]. 116

9.1. RECONSTRUCTION OF THE GEOMETRY

as rigid bodies [Breau et al., 1991]. The NP was shown by Nachemson [1966a] to behave fluid-like and due to its high fluid content (cf. Sec. 2.2) it is assumed to be incompressible. Consequently it is modeled as an inviscid, incompressible fluid that is circumferentially surrounded by the AF. The AF itself is represented by a layered structure of an anisotropic elasto-viscoplastic material. Ground substance and collagen fiber contributions are hom*ogenized in each layer, fibers are not represented explicitly. This choice is justified by a study which shows that a hom*ogenized material model yields smoother strains distributions that if fibers are modeled explicitly (e.g. with one-dimensional nonlinear spring elements) [Eberlein et al., 2001]. The parameters of the constitutive equations [Rubin and Bodner, 2002] were calibrated with data from uniaxial experiments with AF strips (cf. Sec. 8). The resulting model is validated with data from cyclic axial loading and bending experiments (cf. Chap. 7).

9.1

Reconstruction of the geometry

For the simulations of the axial loading experiments the exact geometry of the FSU that gave representative results (cf. Sec. 7.2.1) was reconstructed. To obtain the threedimensional geometry, including not only the surface but as well information on the internal structure (position and height of the IVD, shape and area of the NP, etc.) medical imaging technologies were used. Trials were made with magnetic resonance imaging (MRI) and CT. The image data obtained by those techniques were semi-automatically segmented and the volumes of the Ve and the IVD reconstructed. The resulting polygon models were further processed and finally imported in the FE software ABAQUS (cf. App. D). Image data acquisition One requirement that was identified for the acquisition of volumetric data is the differentiation between variable soft tissues. This is important for a distinction between AF and NP tissue in the disc. Consequently the first imaging technique used was magnetic resonance imaging which gives excellent soft tissue contrasts. Unfortunately this technique was found to be non applicable to the samples as they were used for the mechanical experiments. The brass screws inserted into the lower and upper vertebral endplates to increase the form closure with the polymethylmethacrylate (PMMA) (cf. Sec. 7.1.1) induced artifacts in the acquired images (cf. Fig. 9.2). Brass was chosen on purpose for its first order compatibility with the magnetic resonance imaging technique [Wintermantel and Ha, 2002]. Nevertheless the screws disturb the magnetic field and induce black spots near the position of the screws and distort the images what makes it impossible to reconstruct reliably the geometry. Though the switch to another imaging technique, i.e. CT. One sample was scanned using CT technique with an isotropic resolution of 82 µm (cf. App. D) what would even allow to model the trabecular space in the vertebral bodies (cf. Fig. 9.3). The resulting image data lacks a good differentiation between 117

CHAPTER 9. SIMULATIONS OF AXIAL LOADING AND BENDING

Figure 9.2: Sequential series of four lateral FSU images taken in a magnetic resonance scanner. The crinkly structure at the outer surface is the saline soaked gauze used to keep the sample moist. Brass screws in the outer vertebral endplates disturb the magnetic field and induce artifacts. These are black spots near the screws (at the lower and upper bottoms of all images) and geometric distortion (images 3 and 4, on the right at the level of the IVD).

118

9.1. RECONSTRUCTION OF THE GEOMETRY

brass screws PMMA cylinder upper vertebra IVD lower vertebra

Figure 9.3: Three example image layers from the CT data set. The detail resolution of this scan is extremely high but the disc tissue is not as well differentiated as in the magnetic resonance scans: AF and NP cannot be recognized.

the various soft tissues and AF and NP cannot be distinguished; in the final model their respective volumes and positions were determined with literature data. The brass screws result in shadowlike artifacts in the single images. As the scans were done with layers approximately perpendicular to the longitudinal axis of the sample the artifacts are restricted to the levels where the screws are present, no geometric distortions were introduced. Segmentation and volume generation The data obtained by CT were imported in AMIRA (cf. App. D) a visualization software that allows to segment three dimensional image data and eventually to generate polygonized volume objects. As the total amount of data (approx. 1.3 gigabyte) was too large for the computer system on which AMIRA was running (an office computer with two gigabyte RAM) the two Ve were segmented separately. To this end an automatic segmentation routine based on a threshold gray level value was used. As the resolution of the CT data is so high that the single trabeculae are visible, this procedure resulted in a segmentation with lots of included holes which were removed manually. In the next step the surface of the Ve were computed and the number of elements reduced by a smoothing algorithm. The geometry of the IVD and the relative positions of the AF and NP could not exactly be detected in the CT data. Therefore they were not segmented in AMIRA but approximated in later steps. 119

CHAPTER 9. SIMULATIONS OF AXIAL LOADING AND BENDING

The surface of the PMMA plugs were also reconstructed from which the loading axis was determined. Here the cylindrical part of the plugs was of special interest (cf. Fig. 7.2), unfortunately this is exactly the part that is the most deteriorated by the artifacts induced by the screws. Though, the surfaces were reconstructed only fragmentary (cf. Fig. 9.4), but still this was sufficient. Once the four surfaces (upper and lower Ve and upper and lower PMMA plugs) were reconstructed they were further treated in GEOMAGIC (cf. App. D). In GEOMAGIC important steps were done: 1) the loading axis was determined and the model rotated such that this axis corresponds to the z-axis, 2) the polygonized models from AMIRA were transformed in spline based geometries that can be imported in the FE software ABAQUS and 3) the IVD was created. The alignment of the model with the z-axis is based on the reconstructed surfaces of the PMMA plugs. As the final ingraining of the FSU was done in the materials testing machine, the plugs define the loading axis with a high precision. In order to determine the orientation of the surface reconstructions two coaxial cylinders with the same radius as the PMMA plugs were aligned with them. From those cylinders the orientation was determined and the whole model reoriented such that its loading axis is aligned with the z-axis (cf. Fig. 9.4). The polygonized Ve geometries from AMIRA were further simplified by removing remaining inclusions, holes and other distortions. Afterwards, the spline based geometry was generated by the intermediate steps of contour determination and by patch and grid construction. The third and last step in GEOMAGIC was the construction of the IVD surface. To this end the two Ve were imported and deleted except for the two median endplates. Those two were connected by straight bridges and the intermediary spaces filled. This procedure ensures the good fit of the discs top and bottom surfaces and the adjacent vertebral endplates. The two Ve and the disc were saved in the STEP (STandard for the Exchange of Product model data) format that can be read by ABAQUS where, with help of some MATLAB (cf. App. D) routines, the FE mesh and model was built.

9.2

The finite element model

The main idea behind the FE model is to represent the Ve as rigid bodies and only the IVD as deformable. The disc is further subdivided into the NP and the AF (cf. Fig. 9.5). The position and size of the NP could not be obtained by the CT scans and therefore it was constructed using literature data. According to literature (cf. Sec. 2.3) the NP is modeled as an incompressible fluid (no fluid leaves the NP). The AF is modeled with ten layers for which the fibers have an increasing stiffness from the inner- to the outermost layer. This structure is motivated by the fact that the collagen type-I content in the AF was found to increase from 0 % in the transition zone of NP and AF to 80 % in the outer AF (cf. Sec. 2.2). 120

9.2. THE FINITE ELEMENT MODEL z

(a)

(b)

Figure 9.4: Polygonized reconstruction of the FSU in GEOMAGIC (a). The (common) axis of the two cylinders on top and bottom of the model defines the loading axis (z-axis). The exact orientation of the two cylinders is defined by the best fit on the fragmentary surface reconstruction of the PMMA cylinders in which the Ve were ingrained. In GEOMAGIC the polygon models of the two Ve and the IVD (a) were transformed into spline based models that can be used in ABAQUS (b)

The vertebrae The Ve models from GEOMAGIC were successively imported in ABAQUS as rigid shell objects and meshed with rigid triangular facet (R3D3) and bilinear quadrilateral (R3D4) elements (cf. Fig. 9.5).

9.2.1

The intervertebral disc

The IVD should be meshed with brick elements, a task that is not supported by ABAQUS due to its irregular geometry. Therefore a user routine was written that allows to adapt a meshed template geometry to the geometry of the IVD. The code and a detailed description of the routine is given in Appendix C.4. The spline based IVD model from GEOMAGIC was imported as rigid shell ob121

CHAPTER 9. SIMULATIONS OF AXIAL LOADING AND BENDING

upper vertebra discrete rigid elements (R3D3 / R3D4)

anulus fibrosus solid continuum elements (C3D8H)

nucleus pulposus hydrostatic fluid elements (F3D4)

lower vertebra discrete rigid elements (R3D3 / R3D4)

Figure 9.5: Exploded view drawing of the FSU model. The NP and the upper V are shown partially open to illustrate their shell like (model) character. The left side of the AF presents the orientation of one fiber family1 , the right side shows its ten layers. The “circumferential” nodes of the NP are collocated with the respective AF nodes. Top and bottom nodes of the AF and NP are tied to the “endplates” of the respective Ve. 1

122

The visualization of the fibers was realized with a routine described in Appendix C.7.

9.2. THE FINITE ELEMENT MODEL

ject into the FE software and meshed with rigid triangular facet elements (R3D3) (cf. Fig. C.1a). This mesh served as reference geometry for the mesh adaption routine and was therefore exported via a dummy job that creates an input file from where the nodal coordinates and element connectivities can be extracted. To generate the meshed template geometry the IVD model was imported once more from GEOMAGIC, this time as solid object. A datum plane normal to the z-axis was created on which a sketch of the disc contour was drawn (cf. Fig. 9.6a). This sketch was saved and exported as step file (“.stp”) and served as basis for the computation of the NP and AF layers positions and sizes. A MATLAB routine was written that reads the positions of the control points in the step file, interprets them as vertexes of a polygon and computes an approximate area of the disc. For the definition of the AF layers and NP the control points were shifted along the inward normal at each vertex such that the resulting new areas are 55 % for the NP and equally spaced values for the AF layers up to 100 % for the outermost. According to the literature (cf. Sec. 2.2) the lamellae thicknesses are not constant but depend on the circumferential position. As the modeled AF layers should resemble the lamellae, the relative shift of the control points was computed based on their circumferential position (δrel = (π + ϕ) /π, where ϕ = 0 [rad] at the dorsal and ϕ = π [rad] at the ventral position). For each new set of control points a step file was generated. The step files (in total eleven) were imported in ABAQUS as sketches and the template mesh was generated according to those (cf. Fig. 9.6b). The interior sketches defined the partitions for the ten AF layers and the NP. The geometry was further partitioned in order to get a regular mesh with hybrid brick elements (C3D8H) (cf. Fig. C.1b). This mesh defined the template mesh for the mesh adaption routine. A dummy job created an input file from where the nodal coordinates and element connectivities of the mesh can be extracted. The final mesh for the IVD was obtained by applying the MATLAB routine described in Appendix C.4 to the template mesh. The output of this routine is a mesh where the nodal coordinates of the template mesh are adapted such that it approximates best possible the original disc geometry (cf. Fig. C.1d).

(a)

(b)

Figure 9.6: Topview on spline based model imported from GEOMAGIC into ABAQUS with superimposed approximate contour (a). This contour is used to create a template geometry with NP and ten AF layers (b). 123

CHAPTER 9. SIMULATIONS OF AXIAL LOADING AND BENDING

The nucleus pulposus The NP behaves mechanically like an incompressible fluid (cf. Sec. 2.3). In few studies it is modeled by a (quasi-) incompressible solid [Polikeit et al., 2004] but most often as an inviscid, incompressible fluid [Eberlein et al., 2004]. In the present work it is represented with fluid cavity elements (FCEs). These FCEs are surface elements (F3D4) that define an (incompressible) fluid volume. Here they were applied to the faces of the innermost elements of the AF. Unfortunately this element type cannot be defined in the CAEmodule of ABAQUS — a MATLAB routine was written for this job (cf. App. C.5). This routine deactivates the solid elements that constitute the NP until to this point and defines the new FCEs plus an associated reference points at which the cavity pressure is defined. In the job input file the following lines had to be added to the part definitions in order to use the FCEs (additionally to the element definitions): ∗ F l u i d p r o p e r t y , e l s e t=C a v i t y S u r f a c e , r e f node=Wall−RefPt , type=h y d r a u l i c ∗FLUID DENSITY 1

These lines define a hydraulic (incompressible) fluid on the element set CavitySurface with the reference point Wall−RefPt . The required keyword ∗FLUID DENSITY defines the reference fluid density but is meaningless in the simulations done. The anulus fibrosus The AF mesh was already obtained in the previous step when the template mesh was adapted to the effective IVD geometry. It consists of seven elements in axial and 58 elements in circumferential direction. Ten “concentric” layers constitute a structure similar to the lamellae in the AF. All layers have the same constitutive behavior, which is described by a model presented by Rubin and Bodner [2002, cf. as well Sec. 3.3.3]. An implementation of this model is available as a FORTRAN subroutine for ABAQUS (UMAT) [Papes and Mazza, 2008]. The parameters were taken from a calibration on uniaxial data of AF strips (Tab. 8.2, row “rate dep.”). In contrast to the lamellae, where in each exists only one fiber family, two (hom*ogenized) fiber families are modeled in each AF layer. The only difference between the layers is the stiffness of the (collagen) fibers. This structure is motivated by the fact that the collagen type-I fiber content in the AF was found to increase from 0 % in the transition zone of NP and AF to 80 % in the outer AF (cf. Sec. 2.2). Accordingly the linear parameter of the fiber contribution m3 is linearly decreased from 100 % of its nominal value in the outermost layer to 10 % in the innermost layer. A similar model was already used by Kulak et al. [1976] What still lacks are the orientations of the (collagen) fibers. According to the literature the angle they make with the transversal plane was fixed to ±30◦ (cf. Sec. 2.2). In contrast to some other studies where the fiber angles are varied in circumferential direction [Eberlein et al., 2001], from outer to inner layers [Polikeit et al., 2004], or both 124

9.2. THE FINITE ELEMENT MODEL

[Schmidt et al., 2006], fiber angles are constant throughout the AF. The reason was that no data on the local variation of the fiber angle in the ovine AF was found. Still, the in-plane fiber directions are unknown and need to be defined. To do so, it is assumed that the fibers stay inside their respective layers, i.e. they run “in parallel” with the contours of the AF. This can be achieved if the fibers are interpreted as tangents to streamlines of some flow that is enclosed by the AF. A simple way to achieve such streamlines in a FE software is to run a steady state heat transfer simulation, where the heat is the quantity that flows. Next paragraphs describe in details how the direction vectors were determined. Steady state heat transfer simulations were run with the final IVD mesh and the direction of the heat flux in the AF was taken as the in-plane directions of the fibers. To this end the AF was separated along the median plane into two parts. To the one half and the NP zero conductivity and to the other half a non vanishing conductivity was assigned. At one of the two cutting surfaces through the AF a temperature of 0 K and at the other a temperature of 100 K was defined. The resulting heat fluxes (ABAQUS output variable HFL) at all nodes of the one AF half with non-vanishing conductivity were written to a text file. The conductivities of the two AF halves were switched, the simulation repeated and the heat fluxes written to a second text file. The two text files containing the heat fluxes were read by a MATLAB routine and for each element an average heat flux vector was computed. This vector defines the second basis vector e2 of a local (element wise) coordinate system in which the fiber direction is defined. The third basis vector e3 is defined by the (global) z-direction and is identical to ez . The remaining basis vector e1 is obtained by the crossproduct e1 = e2 × e3 . In this local coordinate system the directions of the two fiber families are given by h i i h ζ (1) = 0 cos[30◦ ] sin[30◦ ] , ζ (2) = 0 cos[30◦ ] − sin[30◦ ] , (9.1a,b) 123

123

respectively. As the two fiber families differ only in the sign of their z-components only one fiber family was defined at this step. Now the direction vectors defined in the local coordinate system were transformed into the global coordinate system with basis vectors ex , ey and ez that is used in the FE simulations. In the implementation of the constitutive model the fiber directions are realized in terms of solution-dependent state variables that are initialized in the FORTRAN subroutine SDVINI (cf. App. C.6). To that end quadruples of element numbers and corresponding fiber directions were written to a text file that is read from the SDVINI. The final FSU model with its different parts, elements and fiber directions is shown in Figure 9.5 as an exploded view drawing.

9.2.2

Boundary conditions

The four instances of the model (upper and lower Ve, AF and NP) were tied together at the intended locations. The AF and the NP share common nodes at their interface and needed therefore no special attention. Ve and IVD were connected in the Assembly module of ABAQUS by (node-to-surface) tie-contacts. 125

CHAPTER 9. SIMULATIONS OF AXIAL LOADING AND BENDING

Axial loading Boundary conditions fix the lower V rigidly to the ground. Displacements in z-direction were applied to the upper V while all other displacements and rotations were restrained. The prescribed displacements of the upper V correspond to the experimental ones measured with the extensometer (cf. Sec. 7.1.3). Bending The kinematics in the bending experiments is more involved than in the axial loading case. Figure 9.7 shows the whole assembly as it was realized in the FE software (compare with Fig. 7.9 which displays the experimental setup). The reference point RP of the upper V was shifted 100 mm above the level of the IVD, what corresponds to the distance of the disc and the lower lever arm. To simulate the lower lever arm, RP was connected to a point P1 shifted horizontally by 150 mm (corresponding to the experimental lever arm length). This connection was realized with a translator that allows only relative displacements along the connecting line and no rotations. The influence of the sensors compliance is taken into account by an additional connector that joins P1 with a second point P2, situated above P1. This second connector restricts relative movements other than along the connecting line but does not prohibit rotations (slot). An elastic modulus mimicking the sensors compliance was assigned to this connector. In a first step of the simulations an initial moment is applied that corresponds to the initial moment applied by the experimental setup (cf. Sec. 7.2.2). To that end the lower V is rigidly fixed to the ground and a force of 4 N charges the point P1, resulting in flexion of the disc. P2 is free to move vertically and follows the displacement of P1. At the end of this step P2 is fixed at its momentary position and the system is allowed to relax for 15 minutes. Then the actual load protocol starts and vertical displacements are applied to the lower V. The imposed displacements result from the recordings of the standard displacement transducer in the experiments.

9.3

Results

Axial loading case simulations failed to converge in the third loading path of the first cyclic compression. The reaction force reached, in the increment before the simulation crashed, a value of −1.1 · 108 N what is five orders of magnitude higher than in the experiment (min. −900 N). The NP pressure surpassed 3 · 105 MPa what is an unnatural high pressure. The reaction force significantly decreased in simulations where the NP was omitted and the full experimental protocol could be simulated. In Figure 9.8 the computed force characteristics for the six compressive cycles is shown. The minimum compressive force attained is about −19 kN, what is still considerably lower than the experimental values. Figure 9.9 presents a detail view on the first compressive cycle and compares the simulations with the experimental results. For the first two or even three loading paths the forces can be considered comparable while for the subsequent loading paths results diverge more and more. 126

9.3. RESULTS

50 mm

slot

150 mm RP

100 mm

translator

Figure 9.7: Sketch of the ABAQUS model used for the bending simulations. The reference point (RP) of the upper V is positioned 100 mm above the IVD center. Points P1 and P2 are the points of force application to the lower lever arm and the fixation of the force transducer, respectively. Connectors (slot and translator) restrict the relative displacements between the three points.

In order to understand why this happens the fiber stretches λf were analyzed (cf. Fig. 9.10). What can be observed is that the fiber stretches are higher for the model with incompressible NP than without. At the last loading path that was simulated for the IVD with NP the median fiber stretch is around 1.0687 [-] and the maximum equals 1.4116 [-]. In case of the IVD without NP the median fiber stretches are for all loading paths smaller one and also the maximum values are smaller than for the IVD with NP. In both simulations fiber stretches reach values that exceed the fiber stretches attained in the experiment with which the constitutive model was calibrated (cf. Sec. 6.4). As a consequence of the above observations the simulation with NP was repeated with more compliant fibers that (the parameter c3 was reduced by 10 % for all AF layers). Neither the reaction force nor the NP pressure values show any significant difference and therefore fibers seem not to be the principal source of the high forces. In case of the bending simulations both models, with and without NP, give similar results and consequently only the results of the model with NP are discussed (cf. Fig. 9.11). Also in this simulation forces are too high, the maximally attained moment is 162 N m 127

128

force [kN]

0.2

0.4

0.6

0.8

1.2

1.4 time [s]

1.6

1.8

2.2

Figure 9.8: Results of the axial loading simulations (IVD without NP). Only the compressive cycles are shown.

−20

−18

−16

−14

−12

−10

−8

−6

−4

−2

·104

2.4

2.6

CHAPTER 9. SIMULATIONS OF AXIAL LOADING AND BENDING

9.3. RESULTS

force [kN]

−1

−0.1

−2

−0.2 −0.3

−3

−0.4

−4

−0.5

−5

force experiment force simulation displacement

−6 0

0.5

displacement [mm]

−0.6 1.5 2 time [s]

2.5

3.5

−0.7

·103

Figure 9.9: Comparison of the first cycle of the axial loading simulations (IVD without NP) with the experiments.

fiber stretch λf [-]

1.4

1.2

1.0

0.8

w/o 1st

w/o 2nd

w/o

3rd 4th # of loading path

w/o 5th

w/o 6th

Figure 9.10: Distribution of the fiber stretches for the simulations of the first axial loading cycle. The distributions are shown at the end of the loading paths for the model with NP (w) and without NP (w/o).

129

130

−50

100

150

Figure 9.11: Results of the bending simulations (IVD with NP).

−150

−100

angle [◦ ]

moment [N m]

6 5 4 3 2 1 0 −1 −2 −3 −4 −5 −6

time [s]

6 ·10

CHAPTER 9. SIMULATIONS OF AXIAL LOADING AND BENDING

9.4. DISCUSSION

(absolute value), compared to 3 N m in the experiment. The initial flexion angle, in the experiments result from of the unbalanced setup and in the simulation obtained by application of an equivalent flexion moment of 0.6 N m, equals 1.6◦ , what is in good accord with the value obtained in the experiment (1.7◦ ). In the four bending cycles the angles were reproduced perfectly for flexion (5.4◦ in-silico and in the experiment) but not as good for extension (−4.6◦ computed vs. −7.8◦ in the experiment). In Figure 9.12 the computed hysteresis curves of the first cyclic tests are compared with the experimental ones. The in-silico constitutive response is characterized by an exponential stiffening and already at rather small angles the bending moment surpasses the maximum experimental bending moment. The two stiffness values at the neutral position (at approximately −1.7◦ ) are about 0.18 N m/◦ for the experiment and 0.32 N m/◦ for the simulation. At higher bending angles the two stiffness values differ much more and read 0.6 N m/◦ vs. 118.0 N m/◦ in extension and 0.9 N m/◦ vs. 132.9 N m/◦ in flexion (experiment vs. simulation, respectively).

9.4

Discussion

Simulations of the axial compression experiments revealed that the FE model as presented in Section 9.2 fails to reproduce the experimentally experienced forces and predicts much higher values. Removal of the NP, which is modeled as an incompressible inviscid fluid, reduces the differences of the reaction forces from five to one orders of magnitude. Under axial load, a perfectly incompressible NP can reduce its height only by expanding laterally, though by bulging the AF. If the AF excels with a high stiffness the whole disc becomes very stiff. In bending where the global deformation of the IVD is wedgelike, with a zone in compression and another in tension, the necessary volume change of the NP is much smaller and therefore less AF distortion is present. This might explain the very similar results of the bending simulations for the FSU with and without NP. The incompressible NP stiffens the FSU in axial loading only if it is surrounded by a material that withstands its high fluid pressure. In the constitutive model of the AF the fibers present the stiffest constituents and therefore their stretch was analyzed. In both simulations, with and without NP, fibers are present that are much more stretched than in the experiments with which the constitutive model was calibrated. A simple numerical experiment proves this not to be the correct explanation for the stiff axial behavior of the FSU with incompressible NP. Fibers were weakened by 10 % in an additional simulation — if the fibers would be the main source of stiffness major changes should occur. This is not the case and therefore it is concluded that the fiber contribution to the total stiffness is only minor. Also the bending experiments could not be reproduced in-silico. Experimentally moments of ±3 N m were applied to the FSU, more than 50 times higher moments were obtained in the FE simulations. The hysteresis curves of the first cyclic bending test 131

CHAPTER 9. SIMULATIONS OF AXIAL LOADING AND BENDING

4 experiment simulation

moment [N m]

2 1 0 −1 −2 −3 −4

−8

−7

−6

−5

−4

−3

−2

−1

angle

[◦ ]

Figure 9.12: Comparison of the hysteresis curves of the first cycle of the bending simulations (IVD with NP) with the experiments.

(cf. Fig. 9.12) show that the simulated structure stiffens at much smaller bending angles than the real structure. While at the neutral position the measured and the modeled stiffness differ only by a factor of about two this factor becomes almost 200 at the maximum angle in extension. This might be due to the constitutive behavior of the AF material or could be a numerical artifact as well. Numerical stiffening can arise from an insufficient discretization and therefore a convergence study is necessary. Eberlein et al. [2004] made a such and found that for the AF four elements in radial, eight in axial and 24 in hoop direction are sufficient to give converged results for spine flexion. Goel et al. [1995] found eight elements in axial and six in radial direction enough to yield good results for the disc bulge under axial load. With regard to this findings it is believed that the present model (with ten elements in radial, seven in axial and 58 in circumferential direction) does not suffer from an insufficient discretization. Another crucial point when it comes to FE simulations is the choice of the element formulation. In the present model, hybrid hexahedral elements with linear shape functions were selected for the AF. These elements should not suffer from numerical stiffening, nevertheless, in a future study their suitability should be checked. A reason for the larger differences in bending than in the axial loading might be the presence of more pronounced displacement gradients. Under axial load the whole disc is axially compressed and bulges radially while for bending one part is in compression, the opposed part in tension and lateral parts do not deform a lot. 132

9.4. DISCUSSION

Besides this numerical issues the fast stiffening could be due to the constitutive model used for the AF. First to be mentioned is the exponential functional form of the equations (cf. Eq. (3.38)) and though their sensitivity to small changes of the arguments. Especially noteworthy are predictions for deformations that lay outside the range for which the model was calibrated. Already small deviations from the calibrated range may result in large errors in the forces/stresses. Furthermore, the stress-strain behavior of soft biological tissues is often prone to errors for calibration purposes. Many soft tissues exhibit a distinct J -shaped stress-strain relation with a very soft initial behavior (cf. Figs. 6.6a and 6.9a), which makes it difficult to determine the initial length of such samples. In the uniaxial experiments with AF strips described in Section 6.4 a preload of 0.01 N was applied that should guarantee non-slack samples at the beginning of the tests. This preload was said to be negligible and consequently it was omitted for the subsequent calibration (cf. Chap. 8). As the material is very soft at small deformations this procedure may lead to a severe underestimation of the stretch and therefore to an overestimation of the samples stiffness. Together with the exponential character of the constitutive equations the increased bending stiffness of the FSU might come into existence. In the axial loading simulations the incompressible NP makes the behavior of the FSU much too stiff. Neglecting the NP changes dramatically the kinematics of the disc as nothing withstands the inward bulging of the AF and seems to reduce the load carrying function of the collagen fibers. While in the simulation with NP almost 75 % of the fibers are in tension, only a little more than 25 % are for the simulations without NP. And, with increasing load, the median fiber stretch becomes larger for the FSU with NP and smaller for the FSU without NP. Though, a realistic FSU model has to include the NP. In reality the NP is a gelatinous part with a high fluid content (cf. Sec. 2.2) which, in fact, renders its behavior incompressible for short time scales. Under long lasting loads fluid leaves the proteoglycan matrix and the NP volume reduces (cf. Sec. 2.3). The slowest cycles in the axial loading lasted for approximately 1 hour, it is questionable if, for this duration, the NP can be assumed incompressible or if not the changing fluid content needs to be incorporated. This could be done either in a phenomenological way, e.g. by using a time dependent bulk modulus, or by a physically based approach as for example by use of a poroelastic model [e.g. Cheung et al., 2003]. The model of a time independent incompressible NP might be a part of the explanation of the unnatural high stiffness in axial compression. Another share could come from the rigid Ve which do not allow for any axial bulging of the NP in the simulation. It is true that the stiffness of the FSU is determined by the sagittal bulging of the IVD at low loads [Hulme et al., 2008] and therefore the endplate stiffness has a minimal influence. But at higher loads it is the deformation of the endplates that defines the compliance. Hence a model that allows the NP to bulge axially too, might provide more realistic results.

133

Chapter

Conclusions and Outlook Summary of achievements This thesis aimed at providing novel insights in the whole procedure of model generation for biological systems. To this end, aspects as the selection of suitable model formulations, the performance of appropriate experiments for parameter determination of constitutive equations and the model validation were analyzed. It was not the objective to provide a model for biomedical applications, but rather to explore continuum mechanics related aspects of the modeling procedure for a complex system as the intervertebral disc. Insights into several aspects were gained and contributions could be made. In Chapter 4 a non physical effect is discussed that might take place with anisotropic constitutive equations that use an additive split of the strain energy into a volumetric and an isochoric part. In fact, simulations involving uniaxial stress configurations in anisotropic fiber-reinforced materials reveal volume growth if this split is applied to the fiber contribution, too. A solution is proposed that solves the problem on the constitutive level instead of using numerical costly methods as the Augmented Lagrangian method. Chapter 5 deals with the calibration of constitutive equations for soft (biological) tissues using planar biaxial materials tests and the herewith connected question of the optimal specimen design. It is shown that results obtained from planar biaxial tests on regular cruciform specimen do not predict the true equibiaxial material behavior and the results are specimen geometry dependent. A numerical optimization scheme is used to maximize the biaxially loaded area by slotting the limbs of the specimen. The findings suggest that already a low number of slots helps to increase the biaxial area significantly. The above results and the smallness of the possible anulus fibrosus (AF) samples led to the conclusion that biaxial experiments are not feasible within the present work. Consequently, only uniaxial experiments with AF-strips were realized, these are described in Chapter 6. While the recorded load curves are qualitatively very similar for different samples they vary a lot in absolute values. Nonetheless, results highlight the presence of a pronounced initial toe-region, a distinct hysteresis and stiffening from one cycle to the next. Chapter 7 covers the axial loading and bending experiments that were done with 135

CHAPTER 10. CONCLUSIONS AND OUTLOOK

ovine functional spinal units (FSUs). Validation of the setup by use of (elastic) steel springs as test pieces revealed several shortcomings that would not have been recognized if experiments with FSUs were directly started. In the following these deficiencies in the mountings, the data acquisition system and the test protocol could be corrected. Repeated cyclic experiments with FSUs at different load rates allow an insight into rate dependency, give information about repeatability and plastic effects. Other phenomenons that can be observed are ratcheting and a distinct hysteresis. The calibration of elasto-viscoplastic constitutive equations with measurement data from uniaxial tensile experiments on AF-strips is exemplified in Chapter 8. The arguments for the applied stepwise procedure, where successively the elastic, the rateindependent and the rate-dependent response was calibrated, were partially negatived by the need for recalibrations in every new step. The final parameter set does not reproduce the experimental data to full satisfaction: it reflects the hardening of the material quite well but not the relaxation behavior. Nevertheless, it represents quite a bit of the complex behavior of the AF tissue. In Chapter 9 the predictions of a finite element model of a FSU are compared with the data gained by axial loading and bending experiments. For both motions the numerical models behave too stiff and fail to reproduce the measured forces. It is shown that the AF fiber contribution to the stiffness is minor. Instead, two other sources are discussed which could lead to an overestimation of the AF stiffness: the low initial stiffness of the tissue which makes it difficult to define the unstressed configuration and the exponential form of the used constitutive equations.

Final remarks The experience gained in the modeling process lead to important considerations that are summarized here: i. Non-physical responses of constitutive models might happen despite those were checked for their robustness. The question of existence and uniqueness of solutions provoked the introduction of mathematical concepts that prove or disprove them for given equations. Also for the three analyzed constitutive models existence and uniqueness of solutions was verified by the notion of polyconvexity [Ball, 1977a]. But those mathematical concepts do not guarantee the physical relevance of the solutions. In fact, all models that use I¯4 as measure of the fiber stretch show a nonphysical response if used to simulate uniaxial tension in a transversely isotropic, fiber-reinforced material. ii. The calibration of constitutive equations requires an explicit assignment of the model parameters to the mode of deformation in nonlinear materials testing. In planar biaxial tests on regular cruciform specimens this is not the case; results do not predict the true equibiaxial material behavior and results are specimen geometry dependent. Only a small central area deforms equibiaxially, the measured global displacement and force values do not well represent the biaxial materials response. Local 136

measures of the deformation field are obtained by optical methods, just, the local stress field is unknown and cannot be measured. A larger area under hom*ogeneous biaxial load would allow to better estimate the stress acting in this area. If biaxial stress and stretch are known, they can be used for the calibration of constitutive equations. iii. A mounting system for stiff (biological) samples was designed and validated by means of experiments with steel springs. These measurements unveiled many shortcomings in the mountings and the experimental setup that would not have been recognized (and corrected) without them. Not enough experiments were made in order to make them statistically representative, but besides the quantitative aspects of the measurements, the applied procedures for validation and measurement are thought to provide interesting insights in how future experiments should be designed.

Outlook A future study that concentrates more on the aspect of a realistic finite element model would need further improvements of the experiments for calibration of constitutive equations and the model itself. If for the representation of the AF such complex constitutive equations as the one of Rubin and Bodner [2002] are used, more experiments that address selected aspects of the model should be realized. Quasi-static experiments could be used to determine the elastic response of the AF tissue. Visco-elastic effects could be characterized by relaxation and cyclic experiments with different load rates. Repetitions of cyclic experiments would give insights in plastic effects that might take place. A higher number of experiments would allow to make statistically representative conclusions about the mechanical behavior of the AF tissue. In the present work, the experiments with AF-strips suffered from the unprecise sample geometries. It was found that the used, custom made, cutting device lacks the needed stability to guarantee well defined and constant sample dimensions. Other techniques could be used to cut specimens of constant thickness as for example a freezing microtome [Hirsch and Galante, 1967]. Few studies were done with unilamellar samplesi but as the single lamella constitutes the unit element of the AF [Holzapfel et al., 2005], experiments with such samples could improve the understanding of the disc mechanics. With more experimental data about the different aspects of the mechanical behavior of the AF tissue, a stepwise calibration of the constitutive equations might unveil its potential. In fact, if the elastic components of the constitutive model are well calibrated an additional calibration of the visco-plastic parts only should give good results. An improvement of the finite element model of the FSU requires, in a first step, a better understanding which parameters have an influence on the predictions of the model. A parametric study on the existing model would help to clarify to which extend single i

Holzapfel et al. [2005]; Panagiotacopulos et al. [1979]; Skaggs et al. [1994]

137

CHAPTER 10. CONCLUSIONS AND OUTLOOK

parameters are responsible for the discrepancies between experiments and simulations. Once the most important parameters are identified they could be optimized. The finite element model of the FSU, as it was presented in this thesis, bases on serial computer tomography scans of an ovine sample. With the applied technique it was not possible to distinguish the different soft tissues that are present in the intervertebral disc, i.e. AF, nucleus pulposus and vertebral endplates. A future model could rely on magnetic resonance image data in which soft tissues are better differentiated. More importantly, a next model should account for the stiffness of the vertebrae and if possible include the vertebral endplates. Such a model is expected to behave more compliant in axial loading and therefore to predict more realistic reaction forces. In a further step a time dependent bulk modulus could be assigned to the nucleus pulposus and, if required, to the AF, too. For the verification of this new model the axial experiments which were realized within this work could be used. The bending experiments should be repeated with an improved setup that is well balanced and does not induce initial moments. Or a new setup could be designed, similar to the one presented by Panjabi et al. [1994], that applies a pure moment to the FSU. Despite the fact that pure bending is seldom present in-vivo, it would help to understand the basic properties of the FSU in this important motion. Additional verification experiments could be added that account for the other motions that are present in daily life: axial torsion and lateral bending. A model that could reliably predict the mechanical response in these four motions would be of high interest.

138

Appendix

Polyconvexity The mechanical problem (A.1a) for a hyperelastic material is to find a displacement field u out of the set of all admissible displacement fields u∗ that fulfills the boundary conditions on ∂Ω and that minimizes the total energy Π of the system. u = arg min Π[ν] ∗ Z ν∈u Z Π[u] = Ψ[grad[u]]dΩ0 − Ω0

Ω0

(A.1a) ρ0 b · udΩ0 −

Z ∂Ω0

t · ud∂Ω0 .

(A.1b)

Here Ω0 and ∂Ω0 are a body and its boundary, respectively, Ψ is the Helmholz free energy, ρ0 the initial density of the body, b and t are volume and surface tractions, respectively and u is the displacement. The question is how the functional form of the strain energy density Ψ[grad[u]] influences the existence and uniqueness of solutions for the minimization problem (A.1a). Several concepts, dealing with the functional dependence of Ψ on the deformation gradient F = grad[u], exist that guarantee the existence and uniqueness of a solution of problem (A.1a): convexity ⇒ polyconvexity ⇒ quasiconvexity ⇒ rank-one-convexity .

Each concept from left to right implicates the consecutive, while the inverse is shown not to be true in full generality [Dacorogna, 2007]. The requirement of convexity with respect to F can be shown to violate fundamental principles of solid continuum mechanics [Ball, 1977b], such that this concept cannot be used. Quasiconvexity (introduced by Morrey Jr. [1952]) is an integral inequality and therefore rather complicated to handle. A more practical notion is the one of polyconvexity [Ball, 1977a,b]. Schr¨oder and Neff [2003, Definition 3.1] say that Definition A.1 (Polyconvexity) “F 7→ W [F] is polyconvex if and only if there exists a function P : M3×3 ×M3×3 ×R 7→ R (in general non-unique) such that W [F] = P [F, adj[F], det[F]] and the function R19 7→ R, ˜ Y˜ , Z) ˜ 7→ P [X, ˜ Y˜ , Z] ˜ is convex for all points X ∈ R3 .” (X, 139

APPENDIX A. POLYCONVEXITY

An immediate consequence of the above Definition A.1 is given by Schr¨oder and Neff [2003, Corollary 3.2]. Corollary A.1 (Additive polyconvex functions) “Let W [F] = W1 [F] + W2 [adj[F]] + W [det[F]]. If Wi , i = 1, 2 are convex in the associated variable respectively and W3 : R+ 7→ R is convex in the associated variable as well, then W is altogether polyconvex.” It can be proven that the following Lemma holds [Schr¨oder and Neff, 2003, Lemma B.9]: Lemma A.1 (Convexity and monotone composition) “Let P : Rn 7→ R be convex and let m : R 7→ R be convex and monotone increasing. Then the function Rn 7→ R, X 7→ m(P (X)) is convex.” Corollary A.1 and Lemma A.1 are very usefull for the following proves of polyconvexity of the constitutive models by Weiss et al. [1996], Holzapfel et al. [2000] and Rubin and Bodner [2002]. Model by Weiss et al. [1996] According to Corollary A.1 it is sufficient for this equations (cf. Eq. (3.34)) to prove the polyconvexity of each of the three energy contributions. The volumetric contribution U (Eq. (3.36a)) is convex only for J smaller equal the natural number e [Hartmann and Neff, 2003, Tab. 4 and text]. ¯ gs (Eq. (3.36b)) involving the first modThe contribution of the ground substance Ψ ified first invariant I¯1 is proven to be polyconvex in Hartmann and Neff [2003, Lemma 2.2]. The term involving the second invariant I¯2 is shown by the same authors to be non-polyconvex hence the parameter c2 has to equal zero in order to assure the overall model being polyconvex. ¯ f (Eq. (3.36c)): Schr¨oder and Neff [2003, Eq. 3.44] prove Polyconvexity of Ψf or Ψ the polyconvexity of the fourth invariants I4 and I¯4 . Though, according to Lemma A.1 ¯ f is a monotone increasing function. This can be done by it is sufficient to prove that Ψ the following procedure: Lets rewrite (3.36c) as ( expc(x−1) −xc x ≥ 1 g[x] = with c ≥ 1. (A.2) 0 x 1 √ √ 2m−2 (2m − 1) ( x − 1) x ≥ 1 > 0. = (2m − 1) √ √ 2m−1 −1/2 x−1 x ( x − 1) | {z } | {z } ≥1

Therefore g 00 ≥ 0 and g is convex.

142

≥1

(A.9)

Appendix

Constitutive Equations in MATHEMATICA A quick but reliable implementation of new constitutive equations is essential when it comes to a comparison of different models or an evaluation of their ability to reflect experimental data. In the present work this was the case in Chapter 4 where different equations were calibrated against each other. An implementation should be fast, user friendly and easy to adapt to new (simple) load cases. MATHEMATICA (cf. App. D) offers a framework that can be used for analytical derivatives, e.g. of the strain energy density Ψ with respect to the right Cauchy-Green deformation tensor. The main difficulty is the definition of the arguments of Ψ to which the differentiation should be made. How to define the deformation tensor C The arguments of Ψ are often the (mixed) invariants Ii of the right Cauchy-Green deformation tensor C (and some structural tensors). The calculation of the stress and stiffness components bases on the computation of the first and second derivative of Ii with respect to C, respectively. In MATHEMATICA this derivatives can easily be obtained as follows: 1. define C in terms of its components, 2. define Ψ and its arguments as functions of the components of C, 3. compute componentwise the derivatives using the Table function. The question is how to define C in order to get the correct result, for example for the derivative of tr[C2 ], as in I2 , where the result is 2C. In the calculation of the stiffness, derivatives of C with respect to itself appear, here the result is known as well: 1 ∂Cij = (δik δjl + δil δjk ) . ∂Ckl 2

(B.1)

This result has all the essential minor and major symmetries that the stiffness tensor C has (i ↔ j, k ↔ l and {i, j} ↔ {k, l}). 143

APPENDIX B. CONSTITUTIVE EQUATIONS IN MATHEMATICA

Since C is symmetric, one could define the right Cauchy-Green tensor with symmetric entries from beginning on:     2C˜11 4C˜12 4C˜13 C˜11 C˜12 C˜13 2 ∂tr[Cs ] ei ⊗ ej =  4C˜12 2C˜22 4C˜23  (B.2a,b) Cs =  C˜12 C˜22 C˜23  , ⇒ ∂ (C ) s ij 4C˜13 4C˜23 2C˜33 C˜13 C˜23 C˜33 As can be seen from Equation (B.2b) this leads to the wrong result. What happens if C is defined without symmetric entries?   C˜11 C˜12 C˜13 Cg =  C˜21 C˜22 C˜23  (B.3) ˜ ˜ ˜ C31 C32 C33 ⇒

∂tr[C2g ] = 2 (Cg )ji , ∂ (Cg )ij

∂ (Cg )ij ∂ (Cg )kl

= δik δjl .

(B.4a,b)

The Equation (B.4a) gives the correct result (remember that eventually C˜ij = C˜ji ), but (B.4b) shows that the stiffness would be wrong. The correct way to define the entries of the right Cauchy-Green tensor is 1) to define a non-symmetric tensor Cg and 2) to compute its symmetric part Csym = 21 Cg + CTg and define all subsequent quantities with respect to Csym . The derivatives are done with respect to the components of Cg . (B.5) Csym = 12 Cg + CTg ⇒

∂tr[C2sym ] = 2 (Csym )ij , ∂ (Cg )ij

∂ (Csym )ij ∂ (Cg )kl

= 12 (δik δjl + δil δjk ) .

(B.6a,b)

Both derivatives are correct. This procedure was used in the MATHEMATICA implementation that is described on the following pages.

144

145

RCG33=RCGF[[3,3]]; σ=J^-1 F.S.Transpose[F]; stiff=J^-1 Table[Sum[ F[[i,w]]F[[j,x]]F[[k,y]]F[[l,z]]Stiff[[w,x,y,z]], {w,3},{x,3},{y,3},{z,3}],{i,3},{j,3},{k,3},{l,3}]; λ1num=Table[1+(λmax-1)/n*i, {i,0,n}]; λ2num=Table[0, i,0,n]; λ2num[[1]] = 1.0; For[i=2,i≤n+1,i++, sol=FindRoot[σ[2,2]==0/.λ1→ λ1num[[i]],{λ2,λ2num[[i-1]]}]; λ2num[[i]]=λ2/.sol; ] σ11num=σ[[1,1]]/.{λ1→ λ1num,λ2→ λ2num};

RCGF = F.Transpose[F]; RCG11=RCGF[[1,1]]; .. .

S=Table[∂RCG[[i,j]] Ψ,{i,3},{j,3}]; Stiff=Table[∂RCG[[k,l]] S[[i,j]],{i,3},{j,3},{k,3},{l,3}]; F={{λ1,0,0},{0,λ2,0},{0,0,λ2}};

.. . Compute here the arguments of Ψ as functions of Csym .. .

Ψ = ... RCG={{RCG11,RCG12,RCG13},{RCG21,RCG22,RCG23}, {RCG31,RCG32,RCG33}}; RCGSym=1/2(RCG + Transpose[RCG]);

MATHEMATICA pseudo code

Compute the numerical values of the Cauchy stress component σ11

Determine λ2 such that σ22 = 0

Define vector λ1 in the range of 1 to λmax with n increments Define vector λ2 with n entries and λ(1)2= 1.0

Compute the spatial stiffness

Compute the Cauchy stress σ

Define the components of C as functions of λi

Compute the second Piola-Kirchhoff stress S Compute the material stiffness C Define the deformation gradient F (here shown for the case of uniaxial tension) Compute C as a function of the stretches λi

Define a symmetric tensor Csym that is used to compute the invariants

Define C as a general tensor, not using its intrinsic symmetry

Define Ψ

Appendix

Routines and ABAQUS Input Files C.1

Customization of the optimization environment

Several times the finite element software ABAQUS (cf. App. D) was used in cooperation with optimization algorithms of MATLAB (cf. App. D) (cf. Chaps. 5 and 8). Whenever this was the case a MATLAB environment developed and described by Hollenstein [2011] was used. This environment allows a simple parametrization of the optimization by splitting it up into the four parts 1) global settings (main.m), 2) settings for the optimization algorithm (kernel.m), 3) control of external programs (COREsearch.m) and 4) error computation (USERerr.m). First optimizations with ABAQUS struggled with several major problems, this were non-existent or incorrect input files (if they were automatically regenerated for each run), lack of free network licenses and old results (.odb-) files. To achieve the required stability for optimizations with several ten thousand iterations, the MATLAB environment was customized by an error handling system. This and other customizations resulting from the experience with the optimization routines are presented in the following. main.m A manual stop of the optimization lets ABAQUS forget some files — at the beginning of each optimization these files are deleted (if they exist). This is especially important in the case of .lck-files, if a such exists no new job can start. delete ( ’ ∗ . rpy ∗ ’ , ’ ∗ . r e c ’ , ’ ∗ . odb ’ , ’ ∗ . l c k ’ ) ;

kernel.m The MATLAB optimization routines do not save intermediate parameter sets and results; if the optimization routine stops unexpectedly (by a crash or manual stop) all results are lost. Therefore a text file is created in which at each iteration the error 147

APPENDIX C. ROUTINES AND ABAQUS INPUT FILES

and the respective parameters are written (in COREsearch to which the file-handle is passed). f i d h i s t = fopen ( ’ h i s t o r y . t x t ’ , ’w ’ ) ;

.. . fclose ( f i d h i s t ) ;

COREsearch.m The error f is initialized with an infinite value, this is done to return an infinite error if the ABAQUS job did not finish successfully. f=i n f ;

.. . The status variable abaqusSTS is set equal −1, this indicates that no job was started yet. As long as abaqusSTS equals −1 it is tried to start the ABAQUS job. This can fail if no ABAQUS licences are available1 . If a result (.odb-) file exists (what means that it could not be deleted previously because it is still opened) ABAQUS will not start the new job and the while-loop is interrupted. After job completion the data of interest is read by the Python script readOdb.py (cf. App. C.8). This script writes the job status to the text file abaqusSTS.sts (1 for job completed successfully and 0 else). abaqusSTS=−1; while abaqusSTS==−1 i f e x i s t ( ’ Uniax2 . odb ’ , ’ f i l e ’ ) break end dos ( ’ abq691 j o b=Uniax2 i n t e r a c t i v e ’ ) ; dos ( ’ abq691 c a e noGUI=readOdb . py ’ ) ; load abaqusSTS . s t s ; end

.. . MATLAB determines the job status by reading the file abaqusSTS.sts. If the status variable abaqusSTS equals 1 the error f is computed. .. . The error and the respective parameter set are written to the text file opened in kernel.m. f p r i n t f ( f i d h i s t , ’%f , ’ , f ) ; 1

If no licences are available, it sometimes happens that ABAQUS gets into the queue and sticks therein. To avoid this unhandiness the following line can be added to the ABAQUS environment file abaqus v6.env: l m l i c e n s e q u e u i n g=OFF With this line included, ABAQUS does not queue and the job aborts — it has to be restarted manually.

148

C.1. CUSTOMIZATION OF THE OPTIMIZATION ENVIRONMENT f p r i n t f ( f i d h i s t , ’%e , ’ , param ) ; f p r i n t f ( f i d h i s t , ’ \n ’ ) ;

.. . The next iteration is prepared by resetting the value in the status file abaqusSTS.sts to −1 and by deletion of remaining ABAQUS files. abaqusSTS=−1; save ( ’ abaqusSTS . s t s ’ , ’ abaqusSTS ’ , ’− a s c i i ’ ) delete ( ’ ∗ . rpy ∗ ’ , ’ ∗ . r e c ’ , ’ ∗ . odb ’ )

149

APPENDIX C. ROUTINES AND ABAQUS INPUT FILES

C.2

ABAQUS job generation: PYTHON script

The PYTHON interface of ABAQUS (cf. App. D) allows the definition of finite element models by use of scripts. This is especially useful in parametric studies where, as in the present case, a parameterized geometry is used that changes according to an external algorithm. The first step of such a PYTHON script is to import the ABAQUS objects and modules that give access to the proprietary functions and operations. The structure is very much based on the module structure in the graphical user interface of ABAQUS (CAE). A clever way to generate such a script is to make use of the CAE to define the model. Doing so, ABAQUS creates a .rpy-file that already contains the structure and most of the commands needed in the script. In the following the (commented) script that was used for the optimization of the cruciform shaped geometry (cf. Chap. 5) is presented. The code # Import Abaqus o b j e c t s from abaqus import ∗ from abaqusConstants import ∗ from odbAccess import ∗ from jobMessage import ∗ # Import modules import part import regionToolset import assembly import step import l o a d import mesh import j o b import odbAccess import v i s u a l i z a t i o n ########################################################## # Input s e c t i o n # Geometry o f t h e c r u c i f o r m T o t l e n g t h =40.0 Armlength =20.0 # Kinematic boundary c o n d i t i o n s BC U1=20.0 BC U2=20.0 # S i z e o f t h e s e e d s f o r t h e mesh S e e d S i z e =0.5 # Element type 150

C.2. ABAQUS JOB GENERATION: PYTHON SCRIPT ElementType=CPS4 ########################################################## # C r e a t e a model . myModel = mdb . Model ( name= ’ Cruciform ’ ) # C r e a t e a new v i e w p o r t i n which t o d i s p l a y t h e model and # the r e s u l t s of the a n a l y s i s . myViewport = s e s s i o n . Viewport ( ’ Cruciform Geometry ’ ) # Create a sketch f o r the base f e a t u r e . mySketch = myModel . S k e t c h ( name= ’ Cruciform shape ’ , s h e e t S i z e =60.) # C r e a t e t h e geometry . mySketch . Line ( p o i n t 1 = ( 0 . 0 , 0 . 0 ) , p o i n t 2 = ( 0 . 0 , T o t l e n g t h ) ) mySketch . Line ( p o i n t 1 = ( 0 . 0 , T o t l e n g t h ) , p o i n t 2 =( Totlength −Armlength , T o t l e n g t h ) ) mySketch . Line ( p o i n t 1 =( Totlength −Armlength , T o t l e n g t h ) , p o i n t 2 =( Totlength −Armlength , Totlength −Armlength ) ) mySketch . Line ( p o i n t 1 =( Totlength −Armlength , Totlength −Armlength ) , p o i n t 2 =( Totlength , Totlength −Armlength ) ) mySketch . Line ( p o i n t 1 =( Totlength , Totlength −Armlength ) , p o i n t 2 =( Totlength , 0 . 0 ) ) mySketch . Line ( p o i n t 1 =( Totlength , 0 . 0 ) , p o i n t 2 = ( 0 . 0 , 0 . 0 ) ) # C r e a t e a two−d i m e n s i o n a l , d e f o r m a b l e p a r t . myCruciform = myModel . Part ( name= ’ Cruciform ’ , d i m e n s i o n a l i t y=TWO D PLANAR, t y p e=DEFORMABLE BODY) # Create the part ’ s base f e a t u r e myCruciform . B a s e S h e l l ( s k e t c h=mySketch ) # Prepare t h e s k e t c h f o r t h e c u t s s = myModel . C o n s t r a i n e d S k e t c h ( name= ’ Cuts ’ , s h e e t S i z e =56.56 , g r i d S p a c i n g =1.41) s . s e t P r i m a r y O b j e c t ( o p t i o n=SUPERIMPOSE) p = myModel . p a r t s [ ’ Cruciform ’ ] p . p r o j e c t R e f e r e n c e s O n t o S k e t c h ( s k e t c h=s , f i l t e r =COPLANAR EDGES) # Read i n t h e p a r a m e t e r s f o r t h e g e n e r a t i o n o f t h e c u t s f = open ( ’ param . param ’ , ’ r ’ ) a=f . r e a d l i n e s ( ) f . close () # Width o f t h e c u t s a0=a [ 0 ] . l s t r i p ( ’ [ ’ ) . r s t r i p ( ’ ] \ n ’ ) . s p l i t ( ’ , ’ ) # Length o f t h e c u t s 151

APPENDIX C. ROUTINES AND ABAQUS INPUT FILES a1=a [ 1 ] . l s t r i p ( ’ [ ’ ) . r s t r i p ( ’ ] \ n ’ ) . s p l i t ( ’ , ’ ) # D i s t a n c e between t h e end o f t h e arms and t h e b e g i n n i n g # of the cuts a2=a [ 2 ] . l s t r i p ( ’ [ ’ ) . r s t r i p ( ’ ] \ n ’ ) . s p l i t ( ’ , ’ ) # D i s t a n c e between t h e c e n t e r l i n e and t h e f i r s t c u t and # between t h e f o l l o w i n g c u t s a3=a [ 3 ] . l s t r i p ( ’ [ ’ ) . r s t r i p ( ’ ] \ n ’ ) . s p l i t ( ’ , ’ ) width=map( f l o a t , a0 ) l e n g t h=map( f l o a t , a1 ) u p p e r l i m i t=map( f l o a t , a2 ) d i s t a n c e=map( f l o a t , a3 ) # D e f i n e t h e geometry o f t h e c u t s d i s t t o t =0.0 f o r i in range ( l e n ( width ) ) : d i s t t o t+=d i s t a n c e [ i ] x p o s l=d i s t t o t −width [ i ] / 2 . 0 x p o s r=d i s t t o t+width [ i ] / 2 . 0 yposu=Totlength −u p p e r l i m i t [ i ]− width [ i ] / 2 . 0 y p o s l=Totlength −u p p e r l i m i t [ i ]− l e n g t h [ i ]+ width [ i ] / 2 . 0 s . Line ( p o i n t 1 =( x p o s l , y p o s l ) , p o i n t 2 =( x p o s l , yposu ) ) s . Line ( p o i n t 1 =( xposr , y p o s l ) , p o i n t 2 =( xposr , yposu ) ) s . ArcByCenterEnds ( c e n t e r =( d i s t t o t , yposu ) , p o i n t 1 =( xposr , yposu ) , p o i n t 2 =( x p o s l , yposu ) ) s . ArcByCenterEnds ( c e n t e r =( d i s t t o t , y p o s l ) , p o i n t 1 =( x p o s l , y p o s l ) , p o i n t 2 =( xposr , y p o s l ) ) s . Line ( p o i n t 1 =( y p o s l , x p o s l ) , p o i n t 2 =(yposu , x p o s l ) ) s . Line ( p o i n t 1 =( y p o s l , x p o s r ) , p o i n t 2 =(yposu , x p o s r ) ) s . ArcByCenterEnds ( c e n t e r =(yposu , d i s t t o t ) , p o i n t 1 =(yposu , x p o s l ) , p o i n t 2 =(yposu , x p o s r ) ) s . ArcByCenterEnds ( c e n t e r =( y p o s l , d i s t t o t ) , p o i n t 1 =( y p o s l , x p o s r ) , p o i n t 2 =( y p o s l , x p o s l ) ) # Make t h e c u t s # I t might happen t h a t t h e p a ramters o f t h e c u t s a r e such # t h a t c u t t i n g i s not p o s s i b l e . In t h i s c a s e w r i t e ’ 0 ’ # t o ’ abaqusSTS . s t s ’ and e x i t . try : p . Cut ( s k e t c h=s ) except AbaqusException , message : output = open ( ’ abaqusSTS . s t s ’ , ’w ’ ) output . w r i t e ( ’ 0 ’ ) output . c l o s e ( ) sys . e x i t () s . unsetPrimaryObject () # Create a ma t e r ia l 152

C.2. ABAQUS JOB GENERATION: PYTHON SCRIPT myModel . M a t e r i a l ( name= ’ NeoHookean Material ’ ) myModel . m a t e r i a l s [ ’ NeoHookean Material ’ ] . H y p e r e l a s t i c ( t e s t D a t a=OFF, t y p e=NEO HOOKE, v o l u m e t r i c R e s p o n s e=VOLUMETRIC DATA, t a b l e = ( ( 1 . 0 , 0 . 0 ) , ) ) mySection = myModel . H o m o g e n e o u s S o l i d S e c t i o n ( name= ’ NeoHookean Section ’ , m a t e r i a l= ’ NeoHookean Material ’ , t h i c k n e s s =1.0) # Assign the s e c t i o n myregion = r e g i o n T o o l s e t . Region ( f a c e s=p . f a c e s ) myCruciform . S e c t i o n A s s i g n m e n t ( r e g i o n=myregion , sectionName= ’ NeoHookean Section ’ , o f f s e t =0.0) # Import t h e p a r t t o t h e assembly myAssembly = myModel . r o o t A s s e m b l y myInstance = myAssembly . I n s t a n c e ( name= ’ C r u c i f o r m I n s t a n c e ’ , p a r t=myCruciform ) # Create a step myModel . S t a t i c S t e p ( name= ’ Step −1 ’ , p r e v i o u s= ’ I n i t i a l ’ , t i m e P e r i o d =1.0 , i n i t i a l I n c =0.1 , d e s c r i p t i o n= ’ ’ , nlgeom=ON, maxNumInc=1000 , e x t r a p o l a t i o n=PARABOLIC) # Define the f i e l d outputs mdb . models [ ’ Cruciform ’ ] . f i e l d O u t p u t R e q u e s t s [ ’F−Output−1 ’ ] . s e t V a l u e s ( v a r i a b l e s =( ’U ’ , ’ S ’ , ’COORD’ ) ) # C r e a t e t h e Boundary C o n d i t i o n s and Loads point =(0.0 ,0.1 ,0.0) edge1=myAssembly . i n s t a n c e s [ ’ C r u c i f o r m I n s t a n c e ’ ] . edges . findAt ( ( point , ) ) p o i n t = ( 0 . 0 , Totlength − 0 . 1 , 0 . 0 ) edge2=myAssembly . i n s t a n c e s [ ’ C r u c i f o r m I n s t a n c e ’ ] . edges . findAt ( ( point , ) ) p o i n t = ( 0 . 1 , Totlength , 0 . 0 ) edge3=myAssembly . i n s t a n c e s [ ’ C r u c i f o r m I n s t a n c e ’ ] . edges . findAt ( ( point , ) ) myModel . DisplacementBC ( name= ’BC−1 ’ , createStepName= ’ I n i t i a l ’ , r e g i o n =[ edge1 , edge2 , edge3 ] , u1=SET) myModel . DisplacementBC ( name= ’BC−2 ’ , createStepName= ’ Step −1 ’ , r e g i o n =[ edge3 ] , u2=BC U2) point =(0.1 ,0.0 ,0.0) edge1=myAssembly . i n s t a n c e s [ ’ C r u c i f o r m I n s t a n c e ’ ] . e d g e s . findAt ( ( point , ) ) p o i n t =( Totlength − 0 . 1 , 0 . 0 , 0 . 0 ) edge2=myAssembly . i n s t a n c e s [ ’ C r u c i f o r m I n s t a n c e ’ ] . e d g e s . findAt ( ( point , ) ) p o i n t =( Totlength , 0 . 1 , 0 . 0 ) 153

APPENDIX C. ROUTINES AND ABAQUS INPUT FILES edge3=myAssembly . i n s t a n c e s [ ’ C r u c i f o r m I n s t a n c e ’ ] . e d g e s . findAt ( ( point , ) ) myModel . DisplacementBC ( name= ’BC−3 ’ , createStepName= ’ I n i t i a l ’ , r e g i o n =[ edge1 , edge2 , edge3 ] , u2=SET) myModel . DisplacementBC ( name= ’BC−4 ’ , createStepName= ’ Step −1 ’ , r e g i o n =[ edge3 ] , u1=BC U1) # Seed t h e p a r t i n s t a n c e . myAssembly . s e e d P a r t I n s t a n c e ( r e g i o n s =(myInstance , ) , s i z e=S e e d S i z e ) # A s s i g n an e l e m e n t type t o t h e p a r t i n s t a n c e . elemType = mesh . ElemType ( elemCode=ElementType , e l e m L i b r a r y=STANDARD) f 1 = myAssembly . i n s t a n c e s [ ’ C r u c i f o r m I n s t a n c e ’ ] . f a c e s f a c e s 1 = f 1 . getSequenceFromMask ( mask=( ’ [#1 ] ’ , ) , ) p i c k e d R e g i o n s =( f a c e s 1 , ) myAssembly . setElementType ( r e g i o n s=p i c k e d R e g i o n s , elemTypes=(elemType , ) ) p a r t I n s t a n c e s =(myAssembly . i n s t a n c e s [ ’ C r u c i f o r m I n s t a n c e ’ ] , ) myAssembly . generateMesh ( r e g i o n s=p a r t I n s t a n c e s ) # Create a job myjobName= ’ C r u c i f o r m O p t i m i z a t i o n ’ myjob=mdb . Job ( name=myjobName , model= ’ Cruciform ’ , t y p e=ANALYSIS , n o d a l O u t p u t P r e c i s i o n=FULL) # Run t h e j o b and h a l t u n t i l j o b has f i n i s h e d myjob . s u b m i t ( ) myjob . w a i t F o r C o m p l e t i o n ( ) # Write t h e o u t p u t s t o ’ abaqus . t x t ’ o1 = s e s s i o n . openOdb ( name= ’ C r u c i f o r m O p t i m i z a t i o n . odb ’ ) s e s s i o n . v i e w p o r t s [ ’ Viewport : 1 ’ ] . s e t V a l u e s ( d i s p l a y e d O b j e c t=o1 ) myodb = s e s s i o n . o d b s [ ’ C r u c i f o r m O p t i m i z a t i o n . odb ’ ] l f i d x=l e n ( myodb . s t e p s [ ’ Step −1 ’ ] . frames )−1 s e s s i o n . f i e l d R e p o r t O p t i o n s . s e t V a l u e s ( p r i n t T o t a l=OFF, printMinMax=OFF) s e s s i o n . w r i t e F i e l d R e p o r t ( f i l e N a m e= ’ abaqus . t x t ’ , append=OFF, s o r t I t e m= ’ Node La be l ’ , odb=myodb , s t e p =0, frame=l f i d x , o u t p u t P o s i t i o n=NODAL, v a r i a b l e =(( ’COORD’ , NODAL, ( (COMPONENT, ’COOR1 ’ ) , (COMPONENT, ’COOR2 ’ ) , ) ) , ( ’ S ’ , INTEGRATION POINT, ( (INVARIANT, ’Max . P r i n c i p a l ’ ) , (INVARIANT, ’ Mid . P r i n c i p a l ’ ) , (INVARIANT, ’ Min . P r i n c i p a l ’ ) , ) ) , ) )

# Get t h e s t a t u s o f t h e computation and w r i t e i t t o 154

C.2. ABAQUS JOB GENERATION: PYTHON SCRIPT # ’ abaqusSTS . s t s ’ : 1 f o r ok , 0 f o r e r r o r . s t a t u s = s t r ( myodb . d i a g n o s t i c D a t a . j o b S t a t u s ) i f l e n ( s t a t u s ) >20: i f ( s t a t u s [ 2 1 ] == ”S” ) : statusOut = ’ 1 ’ else : statusOut = ’ 0 ’ else : statusOut = ’ 0 ’ output = open ( ’ abaqusSTS . s t s ’ , ’w ’ ) output . w r i t e ( s t a t u s O u t ) output . c l o s e ( )

155

APPENDIX C. ROUTINES AND ABAQUS INPUT FILES

C.3

Deformation modes

For soft biological tissues most constitutive models used by now are phenomenological ones. And even those that use some physically based assumptions do not represent the complex structure of the tissues in all details. This means that this models do have (if even) only restricted predictive capabilities when it comes to deformations that are 1) larger and/or 2) different in their mode than those used for the model calibration. The converse argument is that the models should be calibrated in 1) the range of motion and 2) the mode of the deformations that are encountered in the simulations. Therefore the need for knowledge about the modes of deformation that occur in a structure undergoing a given load case. An experienced user might “feel intuitively” what mode of deformation is predominant or which modes of deformation are existent in a certain load case while a less experienced might not. Though it would be favorable to have a qualitative or even a quantitative measure for the modes of deformation. In the following one such is proposed and its implementation for the use in ABAQUS (cf. App. D) discussed. Theory Criscione et al. [2000] presented three invariants Ki , i = 1, . . . , 3 of the natural strain η for the definition of a constitutive model. These three invariants have a physical meaning: K1 : amount of dilatation, K2 : magnitude of distortion, K3 : mode of distortion. The natural strain is defined as η = ln[v] where v is the left stretch tensor according to Equation (3.2).1 The invariants can be written as K1 = tr[η] = ln[J] p K2 = η 0 : η 0 √ K3 = .3 6 det[η 0 /K2 ]

(C.1a) (C.1b) (C.1c)

It can be shown that K1 ∈ (− inf, inf), K2 ∈ [0, inf) and K3 ∈ [−1, 1]. K3 can be used to distinguish the mode of distortion in different regions of a deformed body, where    1 uniaxial extension K3 = 0 pure shear   −1 biaxial extension 1

For two different positive definite but coaxial tensors A and B one can show (using their respective spectral decompositions) that ln [AB] = ln [A] + ln [B]. v is known to be a positive definite and symmetric tensor and therefore 1 1 ln [v] = ln [vv] = ln [b] 2 2 where for the last equality Equation (3.4b) was used. This shows that for the calculation of the logarithmic strains no polar decomposition is necessary.

156

C.3. DEFORMATION MODES

If ηi , i = 1, . . . , 3 are the principal values of the natural strain η, K3 is defined by K3 = −

(η1 + η2 − 2η3 ) (−2η1 + η2 + η3 ) (η1 − 2η2 + η3 ) 3/2

2 (η12 + η22 + η32 − η1 η2 − η2 η3 − η3 η1 )

(C.2)

Now assume a deformation case with an axial stretch λ and two identical lateral stretches λT . In this case the third invariant K3 is given by K3,uniax = sign[log[

λ ]] λT

(C.3)

As can be seen by the above formula K3 = 1 if λ > λT (uniaxial extension) and K3 = −1 if λ < λT (biaxial extension). To see that K3 in fact equals zero for the case of pure shear lets write it in terms of the principal stretches λi : √ 3 6 (ln[λ1 ] ln[λ2 ] ln[λ3 ] − ln[λ1 λ2 λ3 ]) . (C.4) K3 = K2 3 While the first term in brackets equals zero because of the fact that one of the principal stretches is hold constant equal one, the second vanishes because the volume is kept constant during a pure shear deformation. Implementation In a first version Equation (C.2) was implemented in ABAQUS using the subroutine UVARM. This subroutine allows to define user defined variables UVAR. A problem that arises with this implementation is that the variables UVAR are defined at the integration points of an element and not at its nodes. Since the colors in the visualization module of ABAQUS represent the extrapolated (and interpolated) values of a quantity defined at the integration points it happens that the absolute value of UVAR becomes greater than 1 [-]. A second implementation was done using the C++ interface of ABAQUS2 (see the following pages). With C++ it is possible to manipulate .odb-files, i.e. it allows to read and write field variables and to create new ones — also at the element nodes. This is exactly what this implementation does: it reads the principal values of the logarithmic strains “LE” at the element nodes, computes K3 and writes it in a new field “K3” defined at the element nodes. Figure 6.2 gives an example of the use of K3 .

ABAQUS also has a PYTHON interface that offers more comfort with ABAQUS objects. Unfortunately it lacks in efficiency, especially when it comes to loops, e.g. over steps or frames.

157

158

// Reading t h e name o f t h e i n s t a n c e t h a t s h o u l d be used o d b S t r i n g instanceName ( argv [ 2 ] ) ; // Opening t h e odb f i l e odb Odb& Odb = openOdb ( odbName . CStr ( ) ) ; // Opening t h e rootAssembly

// Main program int ABQmain( int argc , char ∗∗ argv ) { // D e f i n i t i o n s f o r t h e new f i e l d o d b S t r i n g FieldName=”K3” ; o d b S t r i n g F i e l d D e s c r i p t i o n=” Deformation modes” ; odb Enum : : odb DataTypeEnum FieldType=odb Enum : : SCALAR; o d b S e q u e n c e S t r i n g FieldComponentName ; FieldComponentName . append ( ”K3” ) ; odb Enum : : o d b R e s u l t P o s i t i o n E n u m F i e l d P o i n t o f D e f i n i t i o n=odb Enum : : NODAL; // Reading t h e name o f t h e odb− f i l e t h a t s h o u l d be used o d b S t r i n g odbName ( argv [ 1 ] ) ;

// Function t o compute K3 a t a node double compK3 ( v a l a r r a y pv ) ;

struct n e w f i e l d s t y p e { int Counter ; v a l a r r a y P r i n c i p a l V a l u e s ; };

#include #include #include #include using namespace s t d ;

The code

APPENDIX C. ROUTINES AND ABAQUS INPUT FILES

o d b A s s e m b l y& r o o t A s s y = Odb . r o o t A s s e m b l y ( ) ; // I n s t a n c e a t which t h e K3 s h o u l d be w r i t t e n o d b I n s t a n c e& i n s t a n c e = r o o t A s s y . i n s t a n c e s ( ) [ instanceName ] ; // Loop o v e r a l l s t e p s i n t h e odb o d b S t e p R e p o s i t o r y I T s t e p I t e r (Odb . s t e p s ( ) ) ; for ( s t e p I t e r . f i r s t ( ) ; ! s t e p I t e r . isDone ( ) ; s t e p I t e r . n e x t ( ) ) { c o u t

[PDF] Continuum Mechanical Investigations of the Intervertebral Disc - Free Download PDF (2024)

References