Eindhoven University of Technology MASTER On the multi-axial testing of membranes Brouwers, H.

(1)

Eindhoven University of Technology

MASTER

On the multi-axial testing of membranes

Brouwers, H.

Award date:

1997

Link to publication

Disclaimer

This document contains a student thesis (bachelor's or master's), as authored by a student at Eindhoven University of Technology. Student theses are made available in the TU/e repository upon obtaining the required degree. The grade received is not published on the document as presented in the repository. The required complexity or quality of research of student theses may vary by program, and the required minimum study period may vary in duration.

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.

• Users may download and print one copy of any publication from the public portal for the purpose of private study or research.

• You may not further distribute the material or use it for any profit-making activity or commercial gain

(2)

On the multi-axial testing of membranes

H. Brouwers

WFW reportnr. 97.062 Master’s thesis

Supervisors: Prof.Dr.Ir. J.D. Janssen Dr.Ir. C.W.J. Oomens Prof. P.J. Hunter

t$ Eindhoven University of Technology, Department of Mechanical Engineering.

Department of Engineering Science.

University of Auckland,

September 1997

(3)

c

(4)

This study has been performed as a final assignment to obtain the Master degree in Biome- chanical Engineering at the Eindhoven University of Technology.

The experimental part of this study has been carried out at the University of Auckland, New Zealand, whereas the preparations, analyses and completion were fulfilled in Eindhoven.

First of all, I want to thank the Dutch Heart Association for the financial support they gave me. Without their contribution, the trip to New Zealand was probably not possible.

Furthermore, I want to thank:

o From the Eindhoven University of Technology:

- Dr.Ir. Cees Oomens and Prof.Dr.Ir. Jan Janssen for their coaching and the chance they gave me to visit New Zealand;

- Ir. Marcel Meuwissen and Ing. Rob Petterson for giving me useful suggestions during my project.

- The University of Eindhoven for giving a financial contribution for my trip to New Zealand.

o From the University of Auckland:

- Prof. Peter Hunter for his invitation, hospitality and supervision during my stay in New Zealand, and his wife Karen for the lovely "Kiwi" meals she sometimes prepared;

- Dr. Paul Charette for solving problems with X V G , even after his departure to Boston, and for his contribution to the experimental setup;

- Dr. Poul Nielsen for his help with setting up the rig during the first six weeks of my stay in Auckland;

- Phd-student Chris Bradley for solving probiems with the software and Phd-student Rob Kirton for his useful ideas about the experimental setup.

Last but not least,

I

want to thank Melissa for the love and support she gave me during this project and for the wonderful time we had during our stay in New Zealand.

Henri Brouwers

i

(5)

This project was performed as a cooperation between the Eindhoven University of Tech- nology, The Netherlands and the University of Auckland, New Zealand.

Eindhoven University of Technology

The University of Auckland

It was financially supported by the Dutch Heart Association (Nederlandse Hartstichting) and the Eindhoven University of Technology.

(6)

Biologische materialen vertonen mechanisch gezien een zeer complex gedrag. Het is wenselijk om dit mechanische gedrag te vangen in een wiskundig model, met behulp waarvan uitspraken gedaan kunnen worden over hoe een bepaald materiaal zich gedraagt onder invloed van een mechanische belasting. Deze modellen kunnen, behalve de geometrie van het materiaal, ook de relatie beschrijven tussen spanningen en rekken in het materiaal, de zogenaamde constitutieve wet. Aangenomen dat er voor elk type materiaal een voldoende nauwkeurige constitutieve wet bestaat, kan deze wet gekwantificeerd worden met een set materiaalparameters. Het is de uitdaging om de waarden van deze materiaalparameters zo nauwkeurig mogelijk te bepalen.

Het is bekend dat biologische weefsels over het algemeen een inhomogeen gedrag vertonen:

de materiaaleigenschappen veranderen van plaats tot plaats. De identificatiemethode (een numeriek-experimentele methode), is uitermate geschikt voor de karakerisering van inhomo- gene materialen. Een proefstuk van een bepaald materiaal en met een bepaalde geometrie, en een numeriek model van dit proefstuk worden onderworpen aan dezelfde randvoorwaar- den, de input, zoals bijvoorbeeld een set voorgeschreven verplaatsingen. Outputvariabelen, zoals reaktiekrachten en het verplaatsingsveld, kunnen gemeten worden. Met behulp van een schattingsalgoritme wordt het verschil tussen de outputvariabelen van het experiment en het numeriek model geminimaliseerd om tot een optimale schatting van parameters te komen. In dit onderzoek is gebruik gemaakt van een multi-axiale trekbank, waarmee de experimentele vrijheid drastisch is toegenomen ten opzichte van een bi-axiale trekbank. Het doel van dit onderzoek is om enig inzicht te verkrijgen in de voor- en nadelen van de toegenomen experimentele vrijheid.

Met behulp van numerieke simulaties, waarbij een experiment als het ware wordt nagebootst, is aangetoond dat het schatten van de materiaalparameters van een lineair elastisch orthotroop materiaal mogelijk is bij gebruik van een multi-axiale trekbank. Daar het accent ligt op de inhomogeniteit van het materiaal, is met opzet gekozen voor een lineair elastisch materiaalge- drag om de rekentijden zo klein mogelijk te houden. Door een beperking in de software tijdens deze simulaties kan alleen de informatie van het verplaatsingsveld worden meegenomen in het schat t ingspro ces.

Experimenten zijn uitgevoerd met twee verschillende materialen, namelijk rubber en pericardium (het hartzakje). Deze proefstukken kunnen met een grote nauwkeurigheid worden opgerekt door middel van zestien actuatoren die zich als een cirkel om het proefstukje heen bevinden. Op deze manier wordt het proefstuk inhomogeen belast. Na het voor-rekken van het proefstukje wordt een referentietoestand bepaald, ten opzichte waarvan alle verplaatsingen (rekken) en krachten gedefinieerd zullen worden. De reaktiekrachten ten gevolge van de

111

...

(7)

extensie van het membraan worden direct gemeten met behulp van zestien krachtopnemers.

Het verplaatsingsveld over het proefstuk wordt met een zeer grote dichtheid gesampled. In tegenstelling tot bij de numerieke simulaties kan nu zowel kracht als verplaatsingsdata worden meegenomen in het schattingsproces. Voor het rubber wordt de neo-Hookean rekenergiefunctie gebruikt als de constitutieve wet en voor het pericardium wordt een extra term toegevoegd aan deze Neo-Hookean rekenergiefunctie.

Het blijkt dat de verplaatsingsdata vrijwel geen informatie bevat omtrend materiaalparameters; gemeten reaktiekrachten bevatten hierentegen wel essentiële informatie. Goede schat- tingsresultaten worden verkregen voor het rubber proefstukje; de resultaten van het pericardium zijn sterk afhankelijk van de gekozen referentietoestand en bijbehorende voorspan- ningen. Verder blijkt dat het schatten beter gaat als de informatie van meerdere belast- ingss t appen wordt samengevoegd.

(8)

Biological materials exhibit a very complex behaviour from the mechanical point of view.

It is preferable to catch this mechanical behaviour into a mathematical model, from which statements can be posed about the behaviour of a certain material under the influence of a mechanical load. These models are able to describe the geometry of the material as well as the relation between strains and stresses in the material, i.e., the constitutive law. Under the assumption that for all material types a sufficiently accurate constitutive law exists, this law can be quantified by means of a set of material parameters. The challenge is to determine the values of these material parameters as accurate as possible.

It is known that biological tissues show in general an inhomogeneous behaviour: the material properties change from place to place. The identi6cation method (a numerical-experimental method), is exceedingly suited for the characterization of inhomogeneous materials. A specimen of a specific material with a certain geometry, and a numerical model of this specimen are subjected to the same boundary conditions, the i n p u t , such as a set of prescribed displacements. O u t p u t variables, such as reaction forces and a displacement field, can be measured.

With an estimation algorithm, the difference between the output variables of the experiment and the numerical model is minimized to obtain an optimal estimate of the parameters. In this study, a multi-axial rig is used with the result that the experimental freedom has increased dramatically with respect to a bi-axial rig. The aim of this study is to obtain insight in the advantages and disadvantages of the increased experimental freedom.

With numerical simulations, which are, so to speak, imitated experiments, is demonstrated that the estimation of material parameters of a linear elastic orthotropic material is possible when a multi-axial rig is used. As the point of interest is formed by the inhomogeneity of the material, linear elastic material behaviour was deliberately chosen to reduce calculation times. Due to a restriction in the software during the numerical simulations, only information of the displacement field can be used in the estimation process.

Experiments are performed with two different materials, namely rubber and pericardium (the heart sack). These specimens can be extended with a high accuracy with sixteen actuators arrayed circularly around the specimen. In this way the specimen will be inhomogeneously loaded. After the pre-stretching of the specimen, a reference state is defined, with respect to which all displacements (strains) and forces will be defined. The reaction forces which are a result of the extension of the membrane are directly measured with sixteen force transducers.

The displacement field over the specimen is sampled with a high density. In contrast to the numerical simulations, both reaction forces and displacements can be used in the estimation process. For the rubber, the neo-Hookean strain energy function is used as the constitutive

V

(9)

model and for the pericardium, an extra term is added to this neo-Hookean strain energy function.

Displacement data seems to contain almost no information about material parameters; on the contrary, measured reaction forces contain essential information. Good estimation results are obtained for the rubber specimen; the results of the pericardium are strongly dependent on the chosen reference state and corresponding pre-stresses. Furthermore, it seems as if the estimation goes better if the information of multiple load cases is combined.

(10)

Preface

i

Samenvatting iii

Summary

V

1

Introduction

1

2

3

Theory

4

2.1 The composition of pericardium

. . .

4

2.4 Constitutive laws in general . . . 6

2.4.1 Nonlinear constitutive laws for soft elastic tissue

. . .

6

2.5 Identification method

. . .

7

2.5.1 Parameter estimation

. . .

8

2.2 Kinematics

. . .

2.3 Stresses

. . .

5 6

Numerical simulations

11 3.2 Methods

. . .

12

3.2.1 Finite Element mesh . . . 12

3.2.2 Boundary conditions and the positions of the observed displacements

.

i 3 3.2.3 Constitutive behaviour of the sample material

. . .

14

3.2.4 Initial parameter values and weighting factors . . . 14

. . .

3.1 Introduction 11

. . .

3.3 Results 15 3.4 Discussion

. . .

17

4

Experiments

19 4.1 Introduction

. . .

19

4.2 Experimental setup

. . .

20

4.2.1 The force transducers . . . 20

4.2.2 Data acquisition and control hardware subsystem . . . 21

4.2.3 The mounting of the specimen on the rig . . . 23

4.3 Determination of the experimental displacement field . . . 26

4.3.1 Selection of markers . . . 26

4.3.2 Transition from an Eulerian to a Lagrangian description . . . 4.4 Finite element modeling and analysis

. . .

28

27

vii

(11)

v111

...

C o n t e n t s

4.5 Performance of the experiments

. . .

30

4.5.1 Rubber membrane . . . 30

4.5.2 Pericardium membrane

. . .

30

4.6 Results

. . .

31

4.6.1 Rubber membrane . . . 31

4.6.2 Pericardium tissue . . . 37

5 Discussion. conclusions and recommendations 42 5.1 Discussion

. . .

42

5.2 Conclusions

. . .

43

5.3 Recommendations

. . .

44

References 46 A Calibration of the sixteen force transducers 48 B Correction of the measured forces 51 C One dimensional analysis of two different strain energy functions C.1 Rubber . . . 53

C.2 Pericardium

. . .

54

53 D Strain estimation from measured marker displacements 56 D.l Introduction

. . .

56

D.2 Estimation of deformed quantities

. . .

57

E Personal experiences in New Zealand 60

(12)

Introduction

Inverse problems can be described as problems from which the causes of certain phenomena are tried to be explained, whereas the results of these phenomena are known by observations.

An example of an inverse problem is the characterization of the mechanical behaviour of materials: What is the cause of a certain deformation of a material when a force is applied to it?

This study is fixed on the aspect of material characterization, with the accent on biological materials.

The knowledge of the mechanical behaviour of biological tissues is important for many reasons. Knowledge of the mechanics of cardiac tissue can, for example, contribute to a better understanding of the coronary blood supply and chemical energetics. It is important to iden- tify the structure and mechanical behaviour of cardiac tissue in order to determine the causes of the rhythm disturbances that lead to failure. In this study, the pericardium (Figure 1.1B) is considered, which is the membraneous sack around the heart. It is known that its role in cardiac mechanics can not be neglected.

In order to be able to obtain predictions of certain phenomena in the human body, numerical models are a useful tool. An example of a numerical model (here: a finite element model) of the human heart is shown in Figure l.ik, where the surface images of a reai heart were

”texture-mapped” onto the model to give it a realistic appearance. The mechanical behaviour of biological tissues is in general non-linear, visco-elastic, anisotropic and inhomogeneous, the latter means that the material properties vary with position. Also, intra- and interindividual variances in mechanical behaviour exist, as well as an age-related dependence of biomechan- ical behaviour, which makes the modeling of biological tissues very compiex. The challenge is to find mathematical descriptions of the stress-strain relations (i. e. constitutive laws) for these types of materials together with the possibility to estimate the parameters of these constitutive laws. In this study, only the non-linear elastic part of materials is taken into account.

The mechanical behaviour of biological tissues can only be determined by means of experiments. In history, these materials were often characterized by means of standard testing methods. These traditional methods have a few disadvantages, especially for biological materials: a large number of experiments must be performed to get a complete description of the mechanical behaviour of the considered piece of tissue. The extraction of test specimen from it’s natural environment may deteriorate the internal structure of the material. Further

1

(13)

2 C h a p t e r i

A B

Figure 1.1: A: finite element model of the heart, developed at the Department of Engineering Science, Auckland; B: schematic cross-section of the heart where the white band around the outer side of the heart represents the pericardium

more, the majority of these experimental methods require a uniform (homogeneous) stress and strain distribution over the specimen, which can not be realized for complex materials.

The identification method [ 5 ] , based on the confrontation between experimental data and a numerical analysis with a resulting adaptation of material parameters, seems to be partic- ularly suitable for the characterization of inhomogeneous materials. It appears to be possible to estimate the material parameters of an inhomogeneous linear elastic orthotropic solid. Van Ratingen [lo] used the identification method together with a linear elastic orthotropic model to estimate the parameters of dog skin with a biaxial testing rig.

Van Kemenade [7] proved by means of numerical simulations that it is possible to estimate the parameters of a non-linear elastic constitutive law with the identification method. Although, when the assumed errors in the observed field data exceeded specific borders, the estimation process diverged. In [14], the identification method was used together with a non-linear strain energy density function to determine the parameters of skin in vivo. Van der Voorden used Lanir’s skin model to estimate the parameters of human skin in vivo [13]. I t seemed to be difficult to obtain reproducible estimates of the parameter values.

In the above named studies, the assumption existed that there is some sort of correlation between the inhomogeneity of the strain field of the deformed sample and the amount of information that this field contains: the more inhomogeneous the strain field is, the more information it might contain, which results in a better convergence of the parameters [5, 101.

It was recommended to use a multi-axial testing machine to obtain the possibility of creating a more inhomogeneous strain field. In [li], a multi-axial tensile device was used to restore the original in vivo geometry of a skin sample after excision and this geometry was used as

(14)

Introduction 3

a reference for the determination of the two-dimensional elastic properties of human skin in terms of a n incremental model at the in vivo configuration.

In this study, the identification method is used together with a multi-axial testing rig, which offers more experimental freedom compared to biaxial testing devices. In Chapter 2, a short outline is given about the structure of the pericardium and the theoretical aspects of constitutive laws and the identification method. Chapter 3 deals with the numerical simulations of the multi-axial experiments, which were performed to obtain insight about experimental conditions. Experiments with membranes of two different materials, rubber and pericardium, together with the results, are elucidated in Chapter 4. Chapter 5 will give a discussion, conclusions and recommendations.

(15)

Chapter 2 Theory

Theoretical aspects which are used in this study will be outlined in this chapter. First, a short description of the structure of the pericardium and its accompanying mechanical behaviour is given in section 2.1. Relevant quantities which can describe the kinematics of and the stresses in a deforming material are introduced in section 2.2 and section 2.3. A general outline about constitutive laws is given in section 2.4. Section 2.5 of this chapter elucidates the identification method in which the parameter estimation process is emphasized.

2.1 The composition of pericardium

In this study, bovine pericardium is taken as the biological material under investigation. P e r i - cardium is the thin membraneous sack which encloses the heart in most mammals. It is made out of two different tissue types: serous pericardium, which is composed of a single smooth layer of epithelial cells and the much thicker fibrous pericardium which is composed of inter- laced collagen and elastin fibers. Pericardium consists of two membranes made out of these two different tissue types, see Figure 2.1. The inner membrane, the visceral pericardium or epicardium, consists of a single layer of serous pericardium which covers the outside of the heart wall. The outer membrane, the parietal pericardium or pericardial sack, consists of an outer layer of fibrous pericardium and an inner layer of serous pericardium. This structure causes the non-linear behaviour of the pericardium. The sack is generally attached to the great vessels, the sternum the vertebral column and the diaphragm. The visceral and parietal paricardia are separated by the pericardial cavity, which is filled with the pericardial fluid.

Of the two pericardial membranes, it is the parietal pericardium which has been studied mostly because it plays a mayor role in cardiac mechanics. In this study, it will simply be referred to as ”pericardium”. The visceral pericardium is not studied on it’s own because it is considered as a part of the heart itself. Pericardium is not essential to life and is often removed during cardiac surgery. However, constrictive pericardial function can be fatal in pathological situations. An example of this is cardiac tamponade, which is an impairment of diastolic filling of the heart caused by an unchecked rise in the inter-pericardial pressure. Also, constrictive pericarditis may occur, which is an inflammation of the pericardium, resulting in constriction of the heart [la]. It it thought that the pericardium has several functions: it

(16)

Theorg 5

Parietal paricardium (pericardial sack)

Fibrous paricardium Serous pericardium

Pericardial cavity

Heart wall

Visceral pericardium (epicardium) Myocardium

(heart muscle) Endocardium Trabeculai!

Figure 2.1: A cross-section of the pericardium, together with a p a r t of the heart wall

holds the heart in a fixed position, prevents over-distension of the heart and protects it from neightbouring organs. Furthermore it enhances the interaction between the ventricles and affects the ventricular pressure-volume curve. Its restriction of free myocardial expansion can effect myocardial wall stresses.

It is obvious that a detailed knowledge of the mechanical behaviour of the pericardium is required to determine its influence on the mechanics of the pumping heart.

2.2 Kinematics

When the mechanical behaviour of highly deformable materials is investigated, it is necessary to have excellent mathematical tools which are able to describe the deformation quantities of the particular materials, the kinematics. The kinematics of a deforming body can be described with the deformation tensor I?, which transforms an infinitesimal material vector dzo in a reference or undeforrned situation' into a vector d z in a deformed situation

This tensor contains information about the rotation in addition to the stretching of material vectors. Another useful tensor is the Green deformation t e n s o r

C

=

@F,

which is independent of rigid body motion and contains information on squared materia1 vector lengths;

*

^indicates

that the tensor is transposed.

The matrix

C

is a second order symmetric tensor which has three i n v a r i a n t s 11, 12 and 13, which remain unchanged under coordinate rotation at a given deformation state. These

'the idea "reference state" or "reference situation" will be elucidated further on in this thesis

(17)

6 Chapter 2

invariants are given by:

where tr(C) is the trace of

C

and d e t ( C ) is the determinant of

C,

which is a measure of volume change. For incompressible materials, det(C) = 1.

To describe strain in a material, the Lagrangian strain tensor

E

is introduced, which is defined by E =

%(C

-

I).

When a material is assumed to be undeformed, the components of E will be zero, since

C

=

1,

so E is defined with respect to the reference state. The eigenvalues of

E

have a physical meaning: They represent the values of the so-called principle strains, and the eigenvectors of E represent their direction, that are the directions in which no shear strains are present. In general, the directions and values of the principle strains will vary over the specimens surface. In a two dimensional situation, two principle strains ( E ~ , E Z ) can be determined in each point on the material. The principle strain domain, which is the collection of a large number of points on the considered material surface, can give information about the inhomogeneity of the strain field over this surface.

2.3 Stresses

Stress is not a uniquely defined quantity as the components of stress can be referred to any of the various material and reference coordinate systems. Furthermore, it can be measured per unit area in either the deformed tissue or the undeformed tissue.

As will be shown in Section 2.4, a constitutive law is a relationship between stress and strain which must be independent of rigid body motion. By introducing the second Piola Kirchhofl stress t e n s o r Y, one is able to describe force per unit undeformed area, which is independent of rigid body motion. Therefore,

Y

is well suited to be used in a constitutive law.

2.4 Constitutive laws in general

Constitutive laws are the mathematical relations between dependent and independent quantities, which are used to describe material behaviour. The relation between stresses (the force per unit area acting on an infinitesimally small plane within the material) and strains (a measure of length change or displacement gradient) in a material is an example of a constitutive law, where the stress is often assumed to be the dependent quantity, caused by the independent quantity strain. A constitutive law must provide a reliable representation of the material behaviour over the range of stresses and strains likely to be encountered in the application of the law. To describe complex materials, for example soft tissues, constitutive laws are often nonlinear functions.

2.4.1

Nonlinear constitutive laws for soft elastic tissue

In formulating a constitutive law for a nonlinearly elastic material certain theoretical guide- lines must be observed [4]. The two axioms of most importance here are:

(1)

the axiom of objectivity, which requires the constitutive law to be invariant to rigid motion of the spatial

(18)

T h e o r y 7

frame of reference; and (2) the axiom of material invariance, which implies certain symmetry conditions dependent on the type of anisotropy exhibited by the material. The first axiom can be satisfied by postulating the existence of a strain energy function W dependent only on the Green deformation tensor

C

or Lagrangian strain tensor

E.

Since

C

(or E) depends only on the material coordinate system, rigid body movement has no influence on the strain energy.

The components of the second Piola-Kirchhoff stress tensor

T

are given by the derivatives of W ( C ) or W(E) with respect to the components of

C

or E, respectively:

For isotropic materials, i.e. the material behaviour in all points of the material is identical in all directions, W is restricted by the second axiom above t o a dependence on the invariants of

C

(11, 1 2 , 13), since rotations of the material coordinate system then have no influence on W . If the material is also incompressible the kinematic constraint 13 = 1 restricts the strain energy to a function of 11 and 1 2 only and the hydrostatic pressure p is introduced into the stress tensor as a Lagrange multiplier to enforce this constraint.

W' = W ( L 1 2 ) + p ( & - 1)

A suitable form of W(Il,12) must then be chosen, based on experimental testing of the (assumed) isotropic, incompressible material. Choosing

W = q ( I 1 - 3)

where e1 is a m a t e r i a l parameter, is suited for certain types of rubber for relatively small strains.

2.5 Identification method

Mixed numerical-experimental methods are suited for the characterization of complex materials. This methods abandon the requirement of the presence of a homogeneous stress and strain field, which leads to more experimental freedom. Hendriks developed a numericai- experimental method for the characterization of the mechanical behaviour of complex materials [5]. This method is based on the assumption that a sufficiently accurate constitutive law is available for the material under investigation. The goal of this method, which will further be referred to as the identification method, is to estimate the parameters in this constitutive law. The identification method, schematically represented by Figure 2.2, is based on the combination of three elements:

o Experiments, which are the measurements of field data (for example displacements and strains) and boundary forces on a multi-axially loaded specimen with arbitrary geometry;

o Finite element modeling of the experiment, resulting in calculated field data and boundary forces;

o Parameter estimation by means of a comparison between experimental and finite element model output, to obtain estimates for the material parameters.

(19)

8 Chapter 2

Figure 2.2: Schematic representation of the identijîcation method

The first two elements (experiments and finite element modeling) will be dealt with further on in this report. The last element (parameter estimation) will now be discussed in more detail.

2.5.1

Parameter estimation

The essential part of the identification method is the parameter estimation process that is performed by means of an iterative procedure. In this iterative procedure, the differences between measured and calculated field data (i.e. residuals) is used to make an improved update of the parameters in the constitutive model. The iteration procedure that was used in [5] is now improved and the improved version will be used in this study.

The unknown material parameters are stored in column

e.

The number of unknown parameters is N . All quantities in the experiment that are s e t by the experimenter, e.g. prescribed displacements, are stored in column ZA, the experimental and model i n p u t . In other words, are the boundary conditions of the experiment as well as the model and 9 is assumed to be exactly known.

When 9 is applied to a specimen, the specimen will undergo a change, e.g., when 9 is a column of prescribed displacements, the specimen will undergo a deformation. From the experiment, certain quantities can be measured, which are related to this change. These quantities are stored in column n ~ , the experimental output. In this study, n~ contains reaction forces and

(20)

Theorv 9

field quantities (displacements) and can therefore be written as:

F m = [ g D ]

where

m”

is a column with measured forces and

mD

is a column with measured displacements.

The lengths of

m”

and m” are respectiveìy the totai number of reaction forces i“vi‘” and the number of displacements M D that are measured in the experiment. i n case of the model, the computed values of the quantities that are measured in the experiments are stored in column y, the model output. The column y is computed via the known constitutive model and an estimate for the parameter values, i e , the model output is dependent on the material parameters: y

-

= y(@).

-

Similar to Equation 2.6, y(@) can be subdivided:

-

y F (O)

[ 1

where yF(@ is a column with calculated forces and yD((3) is a column with calculated displacements. In general, % differs from y:

-

I

m

=

y(@) + 2

(2.8)

where

6

⁼

[ ]

is a column which represents the difference between experimental and model output, the residuals. When the model is assumed to be perfect, ( ^Iwill represent the measurement errors in the experiment.

A weighted least squares estimate of O is obtained by minimizing the following scalar function J , the so-called object function:

+’

²

^po ^-e]-* ^po ^-e]

where 6o represents an initial estimate of the parameter column

e

^and

^y

^contains

^{y F}

^and

- V D . ^{1 -}V F ^{7 -}V D and

w

are symmetrical positive definite weighting matrices. The matrices

y F ,

- V D and are introduced to create the possibility of giving more or less confidence to the measured reaction forces

mF,

the measured displacement field mD and the initial parameter estimate

go.

So J consists of the weighted difference between measured and computed output and the weighted difference between the initial parameter value and the true parameter value. It is wisely to base the weighting matrices on knowledge about the experiment. For example,

y

can be based on information about the inaccuracy of the measurements; in this study,

y F

and

y D

will be taken inverse proportional to the quadratic of the RMS-value of the measurement inaccuracy in mF and

mD

respectively; the ratio of M F and M D (which represents the lengths of

mF

and g D respectively), is also settled in either

y F

or

Y D .

(21)

10 C h a p t e r 2

- .

H F = & - a y F

-

a@

Now J is minimized by setting CY = O, which leads to:

ae

(HF>T(e> Y F [mF

^-

y F ( 8 ) ] + ED [mD

^-

y D ( e ) ] + w [e, -e]

⁼

⁰

(2.10)

-

^-ayF ^I

. . . ayl- -

^-ayD ¹

... dYf

ae, ae, a y D

a!?, aQ,

ay

^{M D}

. . .

-

aY;F

&

. . .

-

ae1 ae,

- -

ae, ae,

-

- . .

H D = -

7 -

a8

. ay -MF

(2.11) Equation 2.10 is non-linear in

e

and can be solved by means of an iteïative procedure. If

e

⁼

^ei ⁺

^{S e i ,}^where^ûiis the i-th parameter estimate and Sei is the difference between estimation and real value, it can be proved that an (i+l)-th estimate û,,, of

e

can be obtained from the i-th estimate ûi as follows:

where

H

contains both

H F

and

H D .

This iterative scheme can only be applied in case of a single set of measured quantities, i.e.,

vF

and

vD

can only contain measurements at one single point in time. Other iteration schemes exist which make it possible to obtain a combined parameter estimate with different sets of measured data, which will not be further discussed. Parameters are updated each iteration step and the iteration process stops when a specified convergence criterion is satisfied. The final p a r a m e t e r estimate û, is then obtained.

(22)

Chapter 3 Numerical simulationis

Numerical simulations of a multi-axial tensile experiment on a membrane are performed with the finite element package DIANA. After a short introduction in the first section, the numer- ical model and the boundary conditions in the simulations are described in the section 3.2, together with the constitutive model that was used to describe the material behaviour of the specimen during the simulations. Also, a choice for initial parameters values and weighting factors are dealt with in this section. In section 3.3, the results of the numerical simulations are shown. This chapter will be ended with a short discussion in section 3.4.

3.1 Introduction

Numerical simulations can be carried out to design a real experiment. The experimental restrictions, such as the size and shape of the specimen and the positions of boundary conditions are known beforehand. With this information, a numerical model with a specified geometry that satisfies these restrictions can be created. Numerical simulations distinguish themselves from real experiments by the fact that the whole experiment is repiaced by a straight forward numerical analysis. During this analysis an assumed constitutive law together with a set of material parameters values are used to characterize the material behaviour of the nu- merical specimen.

The undeformed (unloaded) numerical specimen is taken as the reference state. One or more deformed states of the specimen are created by stretching the specimen by means of one or more sets of boundary conditions

u,

where every set of boundary conditions result in a deformation called a load case. The deformed states can be observed: field quantities and reaction forces, which follow from the numerical analysis, can be used as observed data for a parameter estimation process. Now, the goal of this estimation process is to find the parameters back with the observed data.

In a real experiment, the observational data m =

meZP

is usually not correct: it always contains a m e a s u r e m e n t inaccuracy. For example, when a displacement field is registered by means of a camera, the displacement fieId will contain an error due to discretization of the

11

(23)

12 C h a p t e r 3

continuous displacement field into pixels. Therefore, the expected measurement inaccuracy in the observed quantities meZP in the real experiment have to be added to the observation variables in the deformed numerical model to obtain one numerical observation column

m,,,

per load case;

m,,,

can be used as an input for the parameter estimation process, together with an initial estimate ûi for the parameter values. In theory, a setup of a numerical experiment is successful for material characterization, if the material parameters

ereal

can be traced back from one or more observation columns

m,,,.

In this study, multi-axial tensile experiments with a inhomogeneous orthotropic linear elastic material are simulated in order to obtain knowledge about the influence of the boundary conditions (prescribed displacements) on the parameter estimation process. Due to a restriction in the software at the time of the numerical simulations, only observations of the displacement field y$!!, are used as an input for the parameter estimation process; boundary forces can not be used.

3.2 Methods

3.2.1 Finite Element mesh

The experimental data reflects the distributed characteristics of the material under investigation, e.g., the magnitude, direction and spatial distribution of the displacement field. The data obtained from a number of different stretch steps in an experiment should be linked into one consistent sequence in time.

An appropriate method for the analysis of the large volume of data generated by the experiments is the finite element method. During numerical simulations, the experimental data is generated by the finite element method itself. The finite element method will be used for the creation of the geometry of the specimen and for the characterization and quantification of the constitutive laws of the material under consideration, i.e., parameter estimation. The geometry of the n u m e r i c a l specimen together with the characterization and quantification of its material properties is from now on called the numerical model.

The basic idea behind representing the experimentally measured field quantities (here: displacements) with finite elements is to track the positions of a finite set of predefined material points on the material (the nodes), throughout the deformations. The displacements of any other point on the membrane can then be interpolated from this finite set of tracked points.

The order of the interpolation and the number of tracked material points are selected so that the displacements of all other membrane points can be interpolated to a satisfactory degree o f accuracy.

The finite element mesh that is used for the numerical simulations of the multi-axial experiment is depicted in Figure 3.1. This mesh represents a flat circular specimen shape with a radius

R

= 10 [mm] and a thickness of 1.0 [mm]; these are only proposed dimensions, the dimensions in case of the real experiment were not exactly known during the numerical simulations. The mesh consists of 1184 8-noded plane stress elements (DIANA membrane element CQ16M with eight integration points).

(24)

N u m e r i c a l s i m u l a t i o n s 13

Figure 3.1: Finite element mesh where O represent the points of affection of the boundary conditions, i.e. prescribed displacements

3.2.2

Boundary conditions and the positions of the observed displacements

The positions of the prescribed displacements 21 that are used as an input for the numerical model are given in Figure 3.2A. These displacements are prescribed on the circumference of

0 0 0 0 0 0 0 0 0 O O O O Q O O O O 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Q

0 0 0 O O O O O Q . h

-1.

1 0

o o o o o o ~ o a 0 0 0 0 0 0 0 0 ~ 0 0 0 0 0 0 0 0 0

B

Figure 3.2: A: Position of the prescribed displacements o n the sample and their index; B: Position of the markers in the reference state and the dimensions of the specimen and marker spacing in [mm]

the specimen, perpendicular to this circumference and the mutual distance between their point

(25)

14 Chawter 3

of affection is uniform; All these boundary conditions can be described independently and each individual boundary condition is prescribed in one single node of the finite element mesh.

As a result of the apply of the boundary conditions, fixed points on the surface of the specimen will undergo a displacement between the reference state and a deformed state. A number of fixed points on the material (markers) are followed accurately between two load cases with respect to a fixed coordinate system (z,y) to obtain an observation column The local strains are inferred from the relative displacements of the markers in a certain strain group (see Appendix D). The position of the markers in the reference state are given in Figure 3.2B.

3.2.3

The material that is used in the simulations is assumed to be a physically linear elastic, orthotropic material with varying fibre direction. It is a strong simplification of the real mechanical behaviour of biological tissues. The mechanical behaviour of the material is quantified by means of the following set of independent parameters [lo]:

e,,,,

⁼

e

⁼

(Ei,

E 2 , ~ 1 2 , G 1 2 , bo, b i , b2IT with:

El ^:Young’s modulus in the direction of the fibres;

E2 : Young’s modulus perpendicular to the fibres direction;

v12 : Contraction or Poisson’s ratio;

G 1 2 : Shear modulus;

bo, b l , b2 : Three parameters of a bilinear function, which is used for the description of the fibres direction over the sample domain; The bilinear function is written as (see [ 5 ] ) :

Constitutive behaviour of the sample material

d(x, y) = bo

+

^bi

^* +

^b2

*

y

where

d(z,

y) is the angle between the fibre direction at position (z,y) and the positive z-axis and x and y are the position in [mm] along the z- and y-axis respectively.

An example of a realization of a varying fibre direction over the specimen is shown in Fig- ure 3.3 by choosing specific values for bo, bl and b2. In reality most biological tissues show a less varying fibre direction. In the simulations, the values for bi and b2 have deliberately an unrealistic value. In this manner it can be investigated if the parameters of B strongly inhomogeneous material (in this case a varying fibre direction) can be estimated back from observational data.

3.2.4

Before the parameter estimation process is started with input data

mfUm,

initial values

Bo

for the parameters have to be chosen.

Initial parameter values and weighting factors

As no forces are taken into account, material parameter El will not be estimated and its value will be fixed during the estimation process. Further, the weighting matrix

W

is considered as a diagonal matrix and the main diagonal elements are equal to the squared errors in the initial parameter guess. The weighting matrix

rD

of the measured displacements is set to the inverse of the squared value of the RMS-value of the expected measurement inaccuracy.

This is in fact the RMS-value of the measurement inaccuracies which is also contained in

%,U,.

D

(26)

Numerical simulations 15

Figure 3.3: Fibre direction over the sample geometry

3.3 Results

As stated before, different sets of observations i.e. the displacement data of different deformed states of the specimen, are used to estimate the parameters

ereal

back from an initial parameter estimate

go.

Different deformation states are obtained from different sets of boundary conditions.

Initially, a symmetrical deformation is applied to the membrane. All sixteen prescribed dispiacements (25% of the initia1 radius of the specimen) have the same magnitude and the global deformed shape of the specimen remains circular. In Figure 3.4, the principle strain domain ( E I , E ~ ) , determined via the observed 81 marker displacements m:!, (appendix D), is illustrated. Furthermore, the values of the six parameters ûi are presented as a function of the iteration step i (solid line) and the real parameter values

ereal

are illustrated by the dashed lines. The residual displacement field

C D

after the final iteration step (here: 10) is also shown. The principle strain domain showsthe fact that the strain field is quite inhomogeneous. The parameters v12 and

G12

do not converge nicely to their real values, where the other four parameters do. Due to the absence of modeling errors, the displacement residuals

2”

are randomly distributed and the RMS-value of these residuals equal the EMS-value of the assumed measurement inaccuracy, that is included in TTJ:~,.

Other sets of boundary conditions are applied to the specimen as a trial to improve the estimation results for v12 and

G12.

Asymmetrical prescribed displacements are applied and also cases are simulated where several affection points are released. Below, the most important outcomes of these numerical simulations will be highlighted:

(27)

16 Chawter 3

o 1

008.. ... ....

O 06-

Principle strain domain Residuals

5 0.1 0.15

Li-

E2 ...

'E

^O0 ^Iteration ¹⁰

0 6

' i i

0 3

I

5 10

0.2' O

Iteration

v12

bl

o 3 ...

o 2

o 1

~~~ O0 5 10

Iteration

G12 Z5 ...

20

ii

¹⁵^O0 ¹⁰

Iteration

b2

Iteration

Tc

Figure 3.4: A: Principle strain domain (Q,Q); B: residual displacement field (x 100); G: the course of the parameter values as a function of the iteration step

o In general, the parameters bo, bi and bz converge to their real values bo,real, bi,real and

b2,real, together with Ez. The parameters that give most of the problems are v12 and

GI2 for almost all load cases.

o Furthermore it is obvious that the parameter convergence is better, i.e. final parameter are larger. This fact is estimates

e,

are closer to

ereal,

if the prescribed displacements

(28)

Numerical simulations 17

quite logical since the relative measurement inaccuracy of the marker displacements de- creases when the values of the boundary conditions increase (which results in increasing marker displacements).

o If the assumed measurement inaccuracy exceeds certain borders, the value of the final parameter estimate û, is strongly affected by the realization of the assumed measurement inaccuracy.

o After the last iteration step of the parameter estimation process, displacement residuals are randomly distributed and the RMS-value of these residuals is approximately equal to the RMS-value of the assumed measurement inaccuracy.

o The realization of the boundary conditions 9 has a minor influence on the inhomogeneity of the strain field near the markers and also on the values of the final parameter estimates

8,.

3.4 Discussion

Ideally, it is possible to estimate the parameters of an inhomogeneous orthotropic linear elastic material with the used experimental setup under the assumption that the (assumed) measurement inaccuracy does not exceed certain borders. The assumed measurement inaccuracy in this chapter is a large over-estimation with respect to real experimental noise, so the estimates that follow from the real experimental observations are less influenced by the measurement inaccuracy.

Only displacements of the center part of the specimen are observed. Near the edge of the specimen, the strain field becomes more inhomogeneous. It may be worthwhile to observe the whole surface of the specimen.

When real experiments are carried out, modeling errors will play an important role in the estimation process. A profound inspection of the residuals can provide information about the severity of these errors. In general, modeling errors will have a huge effect on the outcome of the parameter estimation process.

The estimation process that is used to perform the numerical simulations is an old version of the one described in chapter 2. To estimate the parameters of the real experiment, the iterative scheme of chapter 2 will be used, which seems to work faster. The estimation process as described in chapter 2 has two additional advantages: measured forces

mF

can also be used as observed data and the observations of different load cases can be combined. As a consequence, parameters related to absolute stiffnesses in a material can all be estimated.

Also, the resulting parameter estimates are not based on one strain field, but a collection of strain fields. This is expected to give a more fair parameter estimate, because a larger area of stresses and strains will be fitted into the model.

In case fibres are apparent in the material under consideration, it is preferable to measure the course of the fibre direction over the samples surface directly. In this way, fairly good initial estimates of the parameters that describe this course can be used.

(29)

18 Chapter 3

The inhomogeneity of the material here consists of a varying fibre direction. It may be well possible that also stiffness parameters will vary locally over materials surface, especially in case of complex materials (e.g. biological materials). Although this possibility exists in a complex material, it is not taken into account in the numerical simulations in this study.

(30)

Chapter 4 Experiments

In this chapter, the performance and the results of experiments on rubber and pericardium will be elucidated. After an introduction in section 4.1, the experimental setup will be discussed in section 4.2. The determination of displacements in the experiment are explained in section 4.3, and the numerical model of the specimen is dealt with in section 4.4. The performance of the experiments is clarified in section 4.5, where the results of these experiments are presented in section 4.6.

4.1 Introduction

In order to model the mechanical behaviour of a material, the relation between the stresses and the strains needs to be discovered for the material. In fact, these two physical quantities can be derived from two directly measurable variables: forces and material displacements respectively. When complex materials are mechanically tested, it is therefore strongly recom- mendable to fulfill the following conditions:

o The boundary forces, ie., the reaction forces

mF

from the specimen due to the fact that it is extended, should be measured directly.

o The spatial sampling rate of the displacement field over the specimen region must be high enough to ensure that the measurements closely reflect the total displacement field.

In [ 3 ] , a novel method for mechanical testing of complex materials is presented which fulfills the above conditions. An apparatus was developed which could apply any arbitrary planar stress configuration to the membrane under investigation by means of thirteen independently moving actuators which were positioned at the boundaries of the membrane. The displacement measurements were performed using electronic speckle pattern interferometry. Despite the fact that the different components of the apparatus were not implemented to their fullest level of complexity, the technique seemed to be successful.

The apparatus was rebuilt to exploit its full potential. The number of actuators was extended to sixteen and the displacement and strain fields were determined via FT-correlation

19

(31)

20 Chapter

4

techniques, allowing larger deformation steps. Reaction forces of the membrane were measured by strain gauges. Only one experiment was performed with this setup due to a breakage of the force transducers [i]. More robust force transducers are needed to be built, which is done during this study and the final version of these transducers will be described in the next section.

4.2 Experiment al set up

As stated before, the apparatus must be able to apply an arbitrary state of plane stress to the membrane under test. This is achieved by using sixteen actuators arrayed circularly around the membrane, each moving independently in radial direction. The experimental setup roughly consists of a frame and sixteen computer-controlled actuators which are attached to sixteen corresponding force transducers. The frame has a shape of a cylinder with an outer radius of 98 [mm] and a height of 78 [mm] and was designed in a way that it can be used as a bath for biological tissues in the future'. In figure 4.1, the mechanical subsystem is shown; in the schematic drawing (right), eleven actuators are not shown since the system is completely rot at ion-symmet ric.

The actuators consist of a massive cylindrical inner part and a hollow outer part. One tip of the outer part is caught in a hole in the frame and is sticked by screwing a bolt. The inner part of an actuator, which is directly connected to a force transducer, is able to move with respect to the outer part of the actuator in a direction along its axis. The actuators are equally spaced along the edge of the frame. The position of the inner part of the actuators could be read directly on a recess in the outer part of the actuator with an accuracy of about 0.1 [mm].

Two halogen light sources are placed above the specimen and their angles of incidence are optimized to obtain a light distribution over the specimens surface which is as homogeneous as possible. A specimen, which can directly be attached to the force transducers, is extended by a prescribed movement of the actuators and the reaction forces of the specimen can be directly measured. In subsection 4.2.1, a more detailed description of the force transducers is given.

4.2.1

The force transducers

One of the sixteen force transducers that are used during the experiments is shown in Fig- ure 4.2. The force transducer consists of a bronze cantilever beam, which is fixed in a frame at one tip and has one degree of freedom at the other tip, where the force is applied on the membrane via a metal hook. When a force is applied on the hook, the cantilever beam un- dergoes a small deformation in the plane of symmetry of the transducer (the plane of drawing in Figure 4.2), which is almost uniquely related to the applied force. The displacement of the beam can then be measured with a Hall-effect sensor.

'Before this becomes reality, first the size and shape of the force transducers has to be adapted to create this possibility

(32)

Exveriments 21

A

m

l

actuator outer p d i l

// 1

recess

I

B

Figure 4.1: Top view of a part of the mechanical subsystem; A: picture; B: schematic drawing, where the outer edge of the specimen (the membrane) i s depicted b y a dashed line

Two roughly identical magnets (the shaded objects in Figure 4.2) are glued on the frame of the transducer with the N-poles directed towards each other, they will cause a static position- dependent magnetic field in the space between the magnets. This magnetic field strength can then be measured as a function of position by means of a Hall effect sensor (Honeywell, type SS495A), The sensor produces an output voltage VOut which is related to

11?

(Figure 4.3).

When the Hall-effect sensor is mounted on the cantilever beam close to the hook, the position of that part of the beam near the sensor can be accurately determined. Advantages of the sensor is that it is relatively cheap and insensitive for environment influences such as light or temperature. As the cantilever has only one degree of freedom, only normal forces of the membrane can be measured. The sixteen force transducers are calibrated (see Appendix A) and after calibration, an inaccuracy in the measurement of the forces of about 3-4% was apparent.

4.2.2

The entire experiment is controlled in real time by a workstation (IBM RISC System/6000) and a UNIMA (UNiversal Interface Modules and Adapters) box. A flow diagram for control and data acquisition is shown in Figure 4.4. The computer interfaces to external devices via the UNIMA 1/0 system, which is a real time analog and digital 1/0 interface system.

A

graphic user (biaxial control) interface provides the user an easy way of adjusting system parameters, moving the actuators and running the experiments. The biaxial control interface sends commands to the UNIMA 1/0 card which sends digital signals via a 24 bit custom bus to four ”Quad Controller” printed circuit boards (PCBs).

Data acquisition and control hardware subsystem

(33)

22 C h a p t e r

4

A

Intersection A

beam

33

34

J

B

Figure 4.2: T h e force transducer; A : picture; B: drawing, where the connector is left out of consid- eration; dimensions in [mm]

Each quad controller PCB contains four HCTL-1100 chips, which individually control a single actuator. Signals from the PCBs are passed through small actuator control cards connected to the actuators. These consist of a micrometer with a non-rotating spindle coupled without play to a DC-motor. As they are high precision drives they allow displacement steps of less than 0.1 [ p m ] and have a travel range of 50 [mm].

The analog output voltages from the force transducers are sent to voltage dividers, which are positioned on the same board as the actuator control cards. These transform the original output voltage span, which was [-10;10]

[u,

^to^a^{span of}[-350;350] [ m u ; These analog voltage outputs from the board are acquired by EMAP (Electro physiological MAPping system [ S I ) , converted to digital signals and sent via UNIMA to the workstation and the controlling interface where they were timed and recorded to disk.

A black and white camera (640 x 480 pixels, Pulnix TM-9701 Progressive Scanning) was mounted above the test rig. The camera signal is monitored and send to a Silicon Graphics

0 2 workstation. The 0 2 contains a frame-grabber which communicates with the biaxial control interface on the R6000 via the Internet. The biaxial control diagram sends a message to the 0 2 which would save the incoming video image (the frame) to disk with a reference number

(34)

Experiments 23

_ -

Magnetic field strength Output voltage of sensor along 21-axis along 21-axis

Figure 4.3: Relationships of the magnetic field strength

G(ë‘1)

and the output voltage of the sensor VOut(&) to its position along the ë‘l-axis (approximated). Also, the magnetic field lines are drawn

(left); The characters N and S illustrate the North and South pole of the magnets respectively

attached. The IBM/6000 workstation is controlled by a custom application program called X V G [3], running under the UNIX operating system on the Silicon Graphics 0 2 workstation.

4.2.3

In case of the rubber and the pericardium, a circular shaped specimen with a radius of approximately 30 [mm] is made up. For the pericardium, this is the largest possible size. In case of the rubber, a random square black and white dot pattern is generated and directly printed on the specimen by a inkjet printer. With the pericardium, one has to be more careful and the random pattern was hand-painted with Indian ink. The reason for putting on this dot pattern will be explained in section 4.3. Also, pieces of paper with a black and white dot pattern are glued on top of the actuators (see Figure 4.5A). The dot size is chosen approximately three times as large as the pixel size of the camera. Also, the magnetic parts of the force transducers are covered with a black paper to prevent reflections of the magnets.

It is decided not to attach the specimen directly on the hooks of the force transducers but in an indirect way by means of strong surgical wires. Two important reasons for this are:

The mounting of the specimen on the

rig

o The size of the specimens are too small; the minimum radius described by the hooks of the force transducers is about 40 [mrn]; when they move further inwards they will hit each other.