Finite Element Analysis of Bone Remodeling -

Implementation of a Remodeling Algorithm in

MATLAB and ANSYS

Martin Groß

Technical University of Munich1

29th March 2006

1Munich, Germany

---

---

Abstract

English: The process of adaptive bone remodeling can be described mathemati-

cally and simulated in a computer model, integrated with the finite element

method. The main focus of this thesis is the implementation of a bone remod-

eling algorithm in MATLAB and ANSYS on the basis of FEM. The strain

energy density is used as mechanical stimulus. The cortical and trabecular

bone are described as continuous materials with variable density.

This thesis can be divided into four main parts. The first part is due to the ma-

terial properties of cortical and trabecular bone. The second part is about the

remodeling theory and gives an historical review of the developed numerical

approaches up to now. The implementation of the remodeling algorithm in

ANSYS and MATLAB as well as its validation is topic of part three. In last

main part, the algorithm is applied to a 2D FE-model of a human proximal

femur.

Deutsch: Der Vorgang des adaptiven Knochenumbaus kann mathematisch

beschrieben und mit Hilfe der Finiten Elemente Methode an einem Comput-

ermodell simuliert werden. Das Hauptaugenmerk dieser Arbeit liegt in der

Implementierung eines Algorithmus in MATLAB und ANSYS, mit dem der

Knochenumbau auf Basis der FEM simuliert werden kann. Die Dehnungsen-

ergiedichte dient hierbei als mechanischer Stimulus. Die Spongiosa und Ko-

rtikalis werden dabei als kontinuierliches Material mit variierender Dichte

beschrieben.

Die Arbeit kann in vier Hauptteile gegliedert werden. Im ersten Teil wer-

den die Materialeigenschaften von kortikalen und spongiösen Knochen be-

trachtet. Der zweite Teil geht über die Theorie des Knochenumbaus und

gibt einen historischen Rückblick über bisher entwickelte numerische An-

sätze. Die Implementierung des Algorithmus in ANSYS und MATLAB,

sowie dessen Validierung ist Thema des dritten Teils. Im letzten Teil wird der

Algorithmus an einem 2D FE-Modell eines menschlichen proximalen Femurs

angewendet.

Keywords Bone remodeling, finite element method, strain energy density, mechan-

ical stimulus

ii

---

Preface

This Master thesis project was carried out between January 2005 and February 2006

at the Lehrstuhl für Statik at the Technische Universität at München in coopera-

tion with the Labor für Werkstoffkunde und Metallographie (LWM) at the Fach-

hochschule Regensburg.

First of all I like to thank the first supervisor of this thesis, Sebastian Dendorfer, for

his excellent support in every respect. His deep knowledge regarding computational

mechanics and biomechanics in general and especially regarding the remodeling

theory of bone proofed vital for this thesis.

I would also like to thank all my colleagues of the LWM, especially Prof. Dr.

Joachim Hammer, for providing the necessary software and hardware, and for giv-

ing me the opportunity to work at the LWM.

And finally I would like to thank Dr. habil. Manfred Bischoff, who also supervised

this thesis from the Lehrstuhl für Statik.

Martin Groß1, München March 2006

1please feel free to contact me:

in case of any questions concerning

this thesis or the implementation

iii

---

To K.G.

and

to R. & J.G.

iv

---

Contents

Contents

v

List of Figures

viii

List of Tables

xi

1

Introduction

2

1.1

Motivation .

2

1.2

Aim .

3

2

Bone as Material

4

2.1

Composition of Bone .

4

2.2

Difference between Cortical and Cancellous Bone .

5

2.3

Material Properties of Bone Tissue .

8

2.3.1

Cortical Bone Tissue .

8

2.3.2

Cancellous Bone Tissue .

10

3

Bone Remodeling - From Nature to Model

14

3.1

Remodeling Theory .

14

3.1.1

Difference between Modeling and Remodeling .

14

3.1.2

The Remodeling Process .

15

Basic Multicellular Unit .

16

Remodeling Cycle Duration .

20

v

---

CONTENTS

3.1.3

Types of Remodeling .

21

Osteonal Remodeling .

21

Trabecular Remodeling .

22

Endosteal and Periosteal Remodeling .

22

3.2

Theories from the Beginning up to now

.

24

3.2.1

Mechanically Excited Bone Adaption Theories (1865 -1920) 24

3.2.2

Bone Adaptation: General Relationships of Mechanics to

Bone Physiology (1920 - 1970) .

27

3.2.3

Bone Adaptation: Experimental Study of Mechanically

Mediated Bone (1970 - 1984) .

28

3.2.4

Theories of Bone Adaptation:

Numerical Simulations

(1985 to present) .

31

4

Simulation of Remodeling

35

4.1

Implementation of Optimization Algorithm .

35

4.1.1

Reference System .

36

Material Properties .

38

Load

.

38

4.1.2

Optimization Part .

39

Stimulus

.

39

Adaptation Functions .

41

Convergence Criteria .

43

4.2

Results .

44

4.2.1

Reference System 1 .

44

4.2.2

Reference System 2 .

46

4.2.3

Convergence behavior .

47

vi

---

CONTENTS

5

Applications

48

5.1

Modeling of Human Proximal Femur .

49

5.1.1

Load Definitions .

50

5.1.2

Initial Configuration and Boundary Conditions .

52

5.2

Results .

55

5.2.1

Proximal Femur with Stepwise Adaptation Function . . . .

55

Initial Homogeneous Density Distribution .

55

Initial Stochastic Density Distribution .

58

5.2.2

Proximal Femur with Linear Adaptation Function .

61

Initial Homogeneous Density Distribution .

61

Initial Stochastic Density Distribution .

62

5.2.3

Convergence .

63

6

Discussion

65

6.1

Approach .

65

6.2

Details of Algorithm .

65

6.2.1

Building the Model .

65

6.2.2

Development of Bone Structure .

66

6.2.3

Convergence Behavior .

68

6.3

Conclusion and Outlook

.

68

Index

70

Bibliography

72

A Implementation in MATLAB and ANSYS

A

A.1 MATLAB main file run.m

.

A

A.2 ANSYS macro

.

E

A.3 MATLAB function dr.m .

K

B Basic Anatomic Terminology

L

vii

---

List of Figures

2.1

Sketch of some important features of typical long bone, from [93] .

6

2.2

SEM micrograph of ground trabecular vertebra bone, from [76] . . .

6

2.3

3D reconstruction of trabecular bone (4 × 4 × 4mm3 cube) from [4] .

6

2.4

The influence of loading rate on the tensile strength and modulus of

cortical bone from [92] .

8

2.5

Typical stress strain curves for trabecular bone of different densities,

from [43] .

11

2.6

Stress against strain of trabecular bone specimen under compression

from [76] .

12

2.7

Experimental determination of Young′s modulus E against density

. From [85, 14, 2, 57, 64, 79] .

13

3.1

Photomicrograph of an osteonal basic multicellular unit .

16

3.2

Schematic sketch of an osteonal BMU. Cross-sectional view at the

bottom right.

.

17

3.3

The Six Phases of an Osteon′s Lifetime. a) Activation AC, b) Re-

sorption RES, c) Reversal REV, d) Formation FO, e) Mineralization

MI, f) Quiescence QU. from Ott with permission [75] .

19

3.4

BMU activation rate vs. age for human ribs. (From data by [34]) . .

22

3.5

Stress trajectories in curved Culmann crane (left) compared with a

schematic representation of the trabecular pattern in the proximal

femur, from Wolff, 1870 .

24

3.6

Change of trabecular structure in post-fracture , from Wolff, 1870

.

25

3.7

Frost′s description of the different adaptive responses for the ado-

lescent and the adult skeleton.

.

29

viii

---

LIST OF FIGURES

3.8

Density distribution in the femoral head by Fyhrie .

33

3.9

Adaptation function according to equation (3.6) .

34

4.1

Simulation algorithm .

35

4.2

Reference systems .

37

4.3

Higher order 2-D element, from Zienkiewicz [94].

.

37

4.4

Young′s modulus against discretized density according to equa-

tion (4.2)

.

39

4.5

Defined load steps i against time steps .

40

4.6

Overview of implemented adaptation functions .

41

4.7

Resulting material (density) distribution of reference system 1. Fig-

ure 4.7(o) is done with Ole Sigmud′s code from [86]

.

44

4.8

Resulting material distribution of reference system 2. Figure 4.8(o)

is done with Ole Sigmud′s code from [86] .

46

4.9

Convergence plot of reference system 1 .

47

4.10 Convergence plot of reference system 2 .

47

5.1

Anatomy of the human proximal femur. From [89] .

48

5.2

3-D femur model with section plane .

49

5.3

2-D Finite element mesh of the proximal femur with 7124 elements.

49

5.4

Element FLUID79 from [1] .

50

5.5

Proximal femur with muscles.

.

51

5.6

Hip contact force against time for human normal walking from [7]. .

51

5.7

Overview load cases for normal walking .

52

5.8

Initial configurations with two different density distributions .

53

5.9

Remodeling ratio coefficient B(n)in ( g )2MPa-1.

.

53

cm3

5.10 Remodeling progress in human proximal femur with a stepwise

adaptation function and initial homogeneous density distribution. . .

55

5.11 Comparison of v. Mises stresses at initial and converged state. . . .

56

5.12 Comparison of principal stresses at initial and converged state. . . .

57

ix

---