**Biomechanical Computer Models**

K. Engel1, R. Herpers2 and U. Hartmann3

*1German Sport University Cologne 2Bonn-Rhine-Sieg University of Applied Sciences 3University of Applied Sciences Koblenz Germany* 

### **1. Introduction**

16 Will-be-set-by-IN-TECH

92 Theoretical Biomechanics

Ramsay, J. O., Hooker, G. & Graves, S. (2009). *Functional Data Analysis with R and MATLAB*,

Ramsay, J. O. & Silverman, B. W. (2002). *Applied functional data analysis: methods and case studies*,

Ramsay, J. O. & Silverman, B. W. (2005). *Functional data analysis*, Springer series in statistics, 2

Ramsay, J. O., Wickham, H., Graves, S. & Hooker, G. (2010). *fda: Functional Data Analysis*. R

Røislien, J., Skare, O., Gustavsen, M., Broch, N. L., Rennie, L. & Opheim, A. (2009).

Sadeghi, H., Allard, P., Barbier, F., Sadeghi, S., Hinse, S., Perrault, R. & Labelle, H. (2002).

Sadeghi, H., Prince, F., Sadeghi, S. & Labelle, H. (2000). Principal component analysis of the

Simultaneous estimation of effects of gender, age and walking speed on kinematic

Main functional roles of knee flexors/extensors in able-bodied gait using principal

power developed in the flexion/extension muscles of the hip in able-bodied gait,

Use R!, Spinger, New York.

edn, Springer, New York.

package version 2.2.5.

Springer series in statistics, Springer, New York.

URL: *http://CRAN.R-project.org/package=fda*

gait data, *Gait & Posture* 30(4): 441–445.

component analysis (i), *The Knee* 9(1): 47–53.

*Medical Engineering & Physics* 22(10): 703–710.

In the past decade computer models have become very popular in the field of biomechanics due to exponentially increasing computer power. Biomechanical computer models can roughly be subdivided into two groups: multi-body models and numerical models. The theoretical aspects of both modelling strategies will be introduced. However, the focus of this chapter lies on demonstrating the power and versatility of computer models in the field of biomechanics by presenting sophisticated finite element models of human body parts. Special attention is paid to explain the setup of individual models using medical scan data. In order to reach the goal of individualising the model a chain of tools including medical imaging, image acquisition and processing, mesh generation, material modelling and finite element simulation –possibly on parallel computer architectures- becomes necessary. The basic concepts of these tools are described and application results are presented. The chapter ends with a short outlook into the future of computer biomechanics.

The field of biomechanics suffers from one very severe restriction; in general it is not possible for ethical reasons to measure forces and pressure inside the human body. Thus, typical measurement technology in biomechanics works on the interface between body and environment. Force platforms dynamically quantify reaction forces when a person is walking or running across the sensor, electromyography (EMG) monitors action potentials of contracting muscles with electrodes attached to the human skin. The information provided by the measurement technology is very important to investigate the mechanics of movements, but is not sufficient to answer questions like:

	- Ineffective and destructive loading has to be minimised (e.g. by damping through innovative cushions in the shoe sole or improved orthotic devices).
	- Body protection has to be maximised (e.g. through better design of cyclists' helmets or the interior of automobiles).

Biomechanical Computer Models 95

on forward models. The physical forces and moments (e.g. gravitation and external forces, see Figure 2) are the given quantities. Together with geometric data (e.g. length of body segment) and mechanical parameters (e.g. moments of inertia J, centre of mass) the resulting body movement can be computed. The movement is fully described by the laws of

In general, **F**, **a**, **M** and **α** are three-dimensional vectors denoting force, linear acceleration, moments and angular acceleration respectively. The mass m is a scalar and the moment of inertia J is in general a tensor represented by a 3x3 matrix. To be able to solve the underlying differential equations a kinematical chain of simple geometric entities (e.g. an ellipsoid representing the thigh) has to be set up (see Figures 1 and 2). In order to arrive at the quantities of interest (e.g. the load in the hip joint while landing) segment after segment has to be processed. The segments are linked by joints whose degrees of freedom are defined according to those of the anatomical articulations. The rigid body motion (i.e. translation and rotation) of the current segment (e.g. the foot) can be calculated when the external forces (e.g. measured by a force platform) and the body weight of the segment are given. The mechanical influence of the next segment in the chain (e.g. the shank) on the current segment is represented by a reaction force. These reaction forces are most interesting as they give a measure of the load present in the joint during sports exercise. The reaction force determined for one segment is then input (only switching the sign) for the calculation of the rigid body motion of the next segment. After having processed all body segments the complete dynamics for the motion of interest is computable if a time integration scheme is applied. Figure 2 shows the results of such a computation for different kinds of landing after a jump. The reaction forces in the hip joint are depicted for a soft landing (with the use of the body damping system) and a very hard landing (with the body kept stiff). The analysis shows that peak force in the hip joint is more than doubled for the hard landing. The forward model has to be validated against experimental data. Here, the measurement of the landing phase for the different types of landing are monitored and subsequently compared with the body reaction forces calculated with the multi-body model. The simulation results are in good accordance with the experimental data. Another example application of MB models in the field of sports is briefly presented in Figure 3. A back flip jump has been investigated with three different analysis tools. The last row shows the results of a forward solution using a multi-body model of the human body. The computed motion almost perfectly matches the real body movement that is captured by a high-speed camera depicted in the second row. The first

The finite element method (FEM) is a very popular tool in the field of (fluid)-mechanical and electrical engineering to solve partial differential equations (PDEs) on a computer. The first application of the FEM in biomechanics was presented in 1972. FEM was used to build a model of the human femur for orthopaedic applications (Brekelmans et al., 1972). The authors

**F** = *m* **a** (1)

**M** = *J* **α** (2)

Newtonian dynamics that is mathematically described by the following equations.

row shows the typical result of a motion analysis.

**3. The finite element method** 

**3.1 A short history** 

• The patient's outcome and the athlete's output have to be maximised (e.g. by a novel concept for the design a prosthesis or the aerodynamics of a bob sleigh).

The only possible way to answer the questions related to the human body is to set up mechanical models of the human body or at least of human body parts. The setup of a biomechanical model is the only possibility to gain insight into the mechanical behaviour of the inner human body. Thus, modelling has become extremely popular because it –to some extent- heals the fundamental dilemma of biomechanics that has been described at the beginning of the chapter. Therefore, a software tool chain for the automatic generation of *individual* anatomical structures is presented in this book chapter. The setup of models representing medical or sports equipment is also an important task in computational biomechanics. Here, computer analyses are applied to better understand the mechanics of equipment during exercise. Despite the success of modelling techniques in the field of biomechanics, measurement technology still plays an important role: the information provided by the measurements is either input for a model (e.g. reaction forces acting on the foot) or used to validate the model (e.g. comparing the simulation results with EMG data). Thus, it is almost trivial to mention that models are of no scientific use without validation. When looking at the progress being made in the field of biomechanical modelling one can observe two major directions of development:


Both types of models have their advantages and disadvantages: numerical models enable the computation of the whole body's deformation, whereas multi-body models only provide forces at a limited number of body points. But, the price to pay for taking into account the body's deformation is a much higher demand in terms of computer power for the finite element method. This fact partially explains why sometimes parallel computers enter the stage in order to arrive at accurate FE results. Meanwhile, software packages are available that enable the user to use both methods, MBS and FEM simultaneously. The remainder of this text is organised as follows: the next chapter deals with the basic concept of multi body models to then discuss one MBS model in detail. The subsequent chapter gives a short introduction to the finite element method, describes a software tool chain to set up FE models of human body parts starting from medical image data. Chapter 4 provides insight into two example applications. Finally, a short outlook into the future of computer modelling in the field of biomechanics is given.

### **2. Multi body systems in biomechanics**

Within the last two decades multi body models have frequently been applied to solve biomechanical problems (see Figure 1). Some selected MB applications are for instance simulations for impact analysis (Gruber et al.,1998), investigations of the trampoline jumping (Blajer & Czaplicki, 2001), studies about walking and running (Wojtyra, 2003) or the discussion of forces encountered in bicycle sports (Wangerin et al., 2007). Multi body systems can be used in two different ways: as forward or as inverse models. Here, we focus

• Multi body systems (MBS) have been set up yielding important results in prevention

• Numerical models using finite elements (FEM) or computational fluid dynamics (CFD) have been successfully applied to a variety of biomechanical problems such as improving injury prevention (Penrose & Hose, 1998), ameliorating the design of equipment (Dabnichki & Avital, 2006) and for optimising movement techniques

Both types of models have their advantages and disadvantages: numerical models enable the computation of the whole body's deformation, whereas multi-body models only provide forces at a limited number of body points. But, the price to pay for taking into account the body's deformation is a much higher demand in terms of computer power for the finite element method. This fact partially explains why sometimes parallel computers enter the stage in order to arrive at accurate FE results. Meanwhile, software packages are available that enable the user to use both methods, MBS and FEM simultaneously. The remainder of this text is organised as follows: the next chapter deals with the basic concept of multi body models to then discuss one MBS model in detail. The subsequent chapter gives a short introduction to the finite element method, describes a software tool chain to set up FE models of human body parts starting from medical image data. Chapter 4 provides insight into two example applications. Finally, a short outlook into the future of computer

Within the last two decades multi body models have frequently been applied to solve biomechanical problems (see Figure 1). Some selected MB applications are for instance simulations for impact analysis (Gruber et al.,1998), investigations of the trampoline jumping (Blajer & Czaplicki, 2001), studies about walking and running (Wojtyra, 2003) or the discussion of forces encountered in bicycle sports (Wangerin et al., 2007). Multi body systems can be used in two different ways: as forward or as inverse models. Here, we focus

and sports (Gruber et al.,1998; Blajer & Czaplicki, 2001; Wojtyra, 2003).

observe two major directions of development:

modelling in the field of biomechanics is given.

**2. Multi body systems in biomechanics** 

(Wangerin et al., 2007).

• The patient's outcome and the athlete's output have to be maximised (e.g. by a novel concept for the design a prosthesis or the aerodynamics of a bob sleigh). The only possible way to answer the questions related to the human body is to set up mechanical models of the human body or at least of human body parts. The setup of a biomechanical model is the only possibility to gain insight into the mechanical behaviour of the inner human body. Thus, modelling has become extremely popular because it –to some extent- heals the fundamental dilemma of biomechanics that has been described at the beginning of the chapter. Therefore, a software tool chain for the automatic generation of *individual* anatomical structures is presented in this book chapter. The setup of models representing medical or sports equipment is also an important task in computational biomechanics. Here, computer analyses are applied to better understand the mechanics of equipment during exercise. Despite the success of modelling techniques in the field of biomechanics, measurement technology still plays an important role: the information provided by the measurements is either input for a model (e.g. reaction forces acting on the foot) or used to validate the model (e.g. comparing the simulation results with EMG data). Thus, it is almost trivial to mention that models are of no scientific use without validation. When looking at the progress being made in the field of biomechanical modelling one can on forward models. The physical forces and moments (e.g. gravitation and external forces, see Figure 2) are the given quantities. Together with geometric data (e.g. length of body segment) and mechanical parameters (e.g. moments of inertia J, centre of mass) the resulting body movement can be computed. The movement is fully described by the laws of Newtonian dynamics that is mathematically described by the following equations.

$$\mathbf{F} = m \,\,\mathbf{a} \tag{1}$$

$$\mathbf{M} = J \text{ } \mathbf{a} \tag{2}$$

In general, **F**, **a**, **M** and **α** are three-dimensional vectors denoting force, linear acceleration, moments and angular acceleration respectively. The mass m is a scalar and the moment of inertia J is in general a tensor represented by a 3x3 matrix. To be able to solve the underlying differential equations a kinematical chain of simple geometric entities (e.g. an ellipsoid representing the thigh) has to be set up (see Figures 1 and 2). In order to arrive at the quantities of interest (e.g. the load in the hip joint while landing) segment after segment has to be processed. The segments are linked by joints whose degrees of freedom are defined according to those of the anatomical articulations. The rigid body motion (i.e. translation and rotation) of the current segment (e.g. the foot) can be calculated when the external forces (e.g. measured by a force platform) and the body weight of the segment are given. The mechanical influence of the next segment in the chain (e.g. the shank) on the current segment is represented by a reaction force. These reaction forces are most interesting as they give a measure of the load present in the joint during sports exercise. The reaction force determined for one segment is then input (only switching the sign) for the calculation of the rigid body motion of the next segment. After having processed all body segments the complete dynamics for the motion of interest is computable if a time integration scheme is applied. Figure 2 shows the results of such a computation for different kinds of landing after a jump. The reaction forces in the hip joint are depicted for a soft landing (with the use of the body damping system) and a very hard landing (with the body kept stiff). The analysis shows that peak force in the hip joint is more than doubled for the hard landing. The forward model has to be validated against experimental data. Here, the measurement of the landing phase for the different types of landing are monitored and subsequently compared with the body reaction forces calculated with the multi-body model. The simulation results are in good accordance with the experimental data. Another example application of MB models in the field of sports is briefly presented in Figure 3. A back flip jump has been investigated with three different analysis tools. The last row shows the results of a forward solution using a multi-body model of the human body. The computed motion almost perfectly matches the real body movement that is captured by a high-speed camera depicted in the second row. The first row shows the typical result of a motion analysis.

### **3. The finite element method**

### **3.1 A short history**

The finite element method (FEM) is a very popular tool in the field of (fluid)-mechanical and electrical engineering to solve partial differential equations (PDEs) on a computer. The first application of the FEM in biomechanics was presented in 1972. FEM was used to build a model of the human femur for orthopaedic applications (Brekelmans et al., 1972). The authors

Biomechanical Computer Models 97

Fig. 2. Link segment model of the human leg and calculated forces in the hip joint for different kinds of landing (dotted line=soft landing, straight line =hard landing).

**Time in ms**

0 5 10

**FR**

**FG** 

**FG** 

**FG** 

**R**

Fig. 3. The analysis of a back flip by motion analysis, high speed video and multi body

The basic idea of FEM is quite simple: let us assume that a physical problem that can be described by a partial differential equation which is defined on a complex geometry (e.g. the human head or knee). Due to the complex geometry an analytical (i.e. a mathematically correct) solution is impossible. Perfect mathematical solutions can only be determined if the system under consideration has a simple geometry (e.g. a sphere or an ellipsoid). At this point FEM enters the stage: the idea is to approximate the complex geometry by a multitude

**3.2 The basic concept of the finite element method** 

models.

**Multiples of body weight**

Fig. 1. Multi-body models in different visualisations.

set up a two-dimensional isotropic model with roughly 1000 triangular elements and investigated different loading conditions for the human leg. Due to increasing computer power it soon became possible to realise three-dimensional models enabling an improved representation of the complex anatomy of human body parts. Huiskes and Chao, 1983 gave an overview of the first decade of Finite Element modelling in the field of orthopaedic biomechanics. In the last three decades the number of FE models grew exponentially and the application areas cover a large variety of biomechanical problems. Currently, FE models are widely used in research and industry. The application spectrum ranges from trauma to sports-medical biomechanics. This is not only explained by the rapidly increasing availability of fast computers (today with multi-core processors) but also by the advent of efficient algorithms exploiting the advantages of parallel hardware architectures.

Fig. 1. Multi-body models in different visualisations.

set up a two-dimensional isotropic model with roughly 1000 triangular elements and investigated different loading conditions for the human leg. Due to increasing computer power it soon became possible to realise three-dimensional models enabling an improved representation of the complex anatomy of human body parts. Huiskes and Chao, 1983 gave an overview of the first decade of Finite Element modelling in the field of orthopaedic biomechanics. In the last three decades the number of FE models grew exponentially and the application areas cover a large variety of biomechanical problems. Currently, FE models are widely used in research and industry. The application spectrum ranges from trauma to sports-medical biomechanics. This is not only explained by the rapidly increasing availability of fast computers (today with multi-core processors) but also by the advent of

efficient algorithms exploiting the advantages of parallel hardware architectures.

Fig. 2. Link segment model of the human leg and calculated forces in the hip joint for different kinds of landing (dotted line=soft landing, straight line =hard landing).

Fig. 3. The analysis of a back flip by motion analysis, high speed video and multi body models.

### **3.2 The basic concept of the finite element method**

The basic idea of FEM is quite simple: let us assume that a physical problem that can be described by a partial differential equation which is defined on a complex geometry (e.g. the human head or knee). Due to the complex geometry an analytical (i.e. a mathematically correct) solution is impossible. Perfect mathematical solutions can only be determined if the system under consideration has a simple geometry (e.g. a sphere or an ellipsoid). At this point FEM enters the stage: the idea is to approximate the complex geometry by a multitude

Biomechanical Computer Models 99

In general, only a few of the eigenmodes U are of interest. To solve the above equation system sophisticated software is needed. For a dynamical analysis the velocities and accelerations of the mesh nodes are calculated as well. The integration scheme is then

The purpose of the next chapter is to give a detailed description of the workflow for creating

In principle, accurate FE calculations can be carried out on personal computers. However, very complex simulations are candidates for the use of high performance computing. In particular, applications stemming from the field of biomechanics can profit from the use of parallel computers (Hartmann et al., 2001). The complex geometry of anatomical structures often requires a very high spatial resolution of the model leading to possibly huge demand of computer power. In Figure 4 an automatic tool chain to generate *individual* models of human body parts is presented that takes advantage of the application of parallel computers. The software environment has been tested for impact analyses on the human head. It includes the complete chain of tools necessary for the geometric model generation from medical scan data (segmentation, mesh generation and mesh manipulation), computer simulation and visualisation. The software tools are public domain and are available on the

finite element models of human body parts starting with medical image acquisition.

**MU**(*t*) + **DU**(*t*) + **KU**(*t*) = 0 (5)

solving the following system of linear equations.

**4. FE Applications in biomechanics** 

web (Berti,2010).

**4.1 A high resolution model for impact studies** 

Fig. 4. From medical scan data to finite elements.

of simple geometric entities (e.g. triangular, tetrahedral or hexahedral elements). Each element is then individually processed. Within an element adequate functions (so-called basis functions) are chosen that are used to interpolate the *unknown* function. The general procedure of FEM is summarised in six steps:

a. Generation of a discrete representation of the object of interest.

The mesh generation step is crucial for the rest of the analysis as the spatial resolution and the element quality to a large extent affect the robustness and the accuracy of the simulation result. In order to achieve high accuracy in a physically interesting region (e.g. the area of large deformation in an impact scenario) a mesh refinement in such regions of special interest is indicated. Mesh refinement means to increase the number of elements and likewise the number of mesh nodes (in the simplest case nodes can be regarded as the corner points of an element). Thus, FE meshes can have millions of mesh nodes. The result of 3D mechanical FE study is the displacement function of the nodes (i.e. the body's deformation) in three dimensions of space.

b. Defining the physical behaviour

In order to set up a realistic model the physical behaviour of the object has to be integrated. At this stage two important decisions have to be made. An appropriate material model (e.g. elastic or viscoelastic behaviour) has to be chosen and the boundary conditions (e.g. fixed or free) have to be selected. The simulation results are very sensitive to these settings. Human soft tissue can be modelled in the simplest case by assuming linear elastic behaviour. Here, two mechanical parameters (Young's modulus and Poisson's ratio) suffice. However, the results provided by such a simple model might not be satisfactory.

c. Setup of the stiffness matrix

For each finite element a stiffness matrix has to be set up. A prerequisite for this step is the choice of adequate interpolation or basis functions within an element. Depending on the number of nodes of the element (e.g. a tetrahedral element can be defined to have four nodes on each of its corners or ten if nodes on the connecting lines are added) the interpolation functions are chosen to be linear or quadratic. Increasing the polynomial order improves the accuracy of the simulation results, but the price to pay is a higher demand in terms of computer resources. The element stiffness matrices are then summed up (taking into account the element connectivity) to form the global stiffness matrix K that represents the complete system. If a transient or dynamic FE analysis is desired a mass matrix M and a damping matrix D have to be assembled to describe the inertial and the damping behaviour of the system respectively.

d. Numerical solution of the problem

After the system matrices M, D and K have been assembled a numerical equation solver is started to calculate the vector of nodal displacements U of the body when an external force defined by a force F is applied. The simplest FE analysis is the computation of the new static equilibrium for a body under load. The equation system to solve is then given by

$$\mathbf{K}\mathbf{U}=\mathbf{F}\,.\tag{3}$$

Sometimes it is of great interest for the developer what the eigenfrequencies ω of an object are and how the corresponding eigenmodes U look like. Then the eigenvalue equation system looks like

$$\mathbf{K}\mathbf{U} = \boldsymbol{\phi}^2 \mathbf{M}\mathbf{U} \tag{4}$$

of simple geometric entities (e.g. triangular, tetrahedral or hexahedral elements). Each element is then individually processed. Within an element adequate functions (so-called basis functions) are chosen that are used to interpolate the *unknown* function. The general

The mesh generation step is crucial for the rest of the analysis as the spatial resolution and the element quality to a large extent affect the robustness and the accuracy of the simulation result. In order to achieve high accuracy in a physically interesting region (e.g. the area of large deformation in an impact scenario) a mesh refinement in such regions of special interest is indicated. Mesh refinement means to increase the number of elements and likewise the number of mesh nodes (in the simplest case nodes can be regarded as the corner points of an element). Thus, FE meshes can have millions of mesh nodes. The result of 3D mechanical FE study is the displacement function of the nodes (i.e. the body's

In order to set up a realistic model the physical behaviour of the object has to be integrated. At this stage two important decisions have to be made. An appropriate material model (e.g. elastic or viscoelastic behaviour) has to be chosen and the boundary conditions (e.g. fixed or free) have to be selected. The simulation results are very sensitive to these settings. Human soft tissue can be modelled in the simplest case by assuming linear elastic behaviour. Here, two mechanical parameters (Young's modulus and Poisson's ratio) suffice. However, the

For each finite element a stiffness matrix has to be set up. A prerequisite for this step is the choice of adequate interpolation or basis functions within an element. Depending on the number of nodes of the element (e.g. a tetrahedral element can be defined to have four nodes on each of its corners or ten if nodes on the connecting lines are added) the interpolation functions are chosen to be linear or quadratic. Increasing the polynomial order improves the accuracy of the simulation results, but the price to pay is a higher demand in terms of computer resources. The element stiffness matrices are then summed up (taking into account the element connectivity) to form the global stiffness matrix K that represents the complete system. If a transient or dynamic FE analysis is desired a mass matrix M and a damping matrix D have to be assembled to describe the inertial and the damping behaviour

After the system matrices M, D and K have been assembled a numerical equation solver is started to calculate the vector of nodal displacements U of the body when an external force defined by a force F is applied. The simplest FE analysis is the computation of the new static

Sometimes it is of great interest for the developer what the eigenfrequencies ω of an object are and how the corresponding eigenmodes U look like. Then the eigenvalue equation

> **KU MU** <sup>2</sup> = ω

**KU** = **F** . (3)

(4)

equilibrium for a body under load. The equation system to solve is then given by

procedure of FEM is summarised in six steps:

deformation) in three dimensions of space. b. Defining the physical behaviour

c. Setup of the stiffness matrix

of the system respectively.

system looks like

d. Numerical solution of the problem

a. Generation of a discrete representation of the object of interest.

results provided by such a simple model might not be satisfactory.

In general, only a few of the eigenmodes U are of interest. To solve the above equation system sophisticated software is needed. For a dynamical analysis the velocities and accelerations of the mesh nodes are calculated as well. The integration scheme is then solving the following system of linear equations.

$$\mathbf{M}\ddot{\mathbf{U}}(t) + \mathbf{D}\dot{\mathbf{U}}(t) + \mathbf{K}\mathbf{U}(t) = 0\tag{5}$$

The purpose of the next chapter is to give a detailed description of the workflow for creating finite element models of human body parts starting with medical image acquisition.

## **4. FE Applications in biomechanics**

### **4.1 A high resolution model for impact studies**

In principle, accurate FE calculations can be carried out on personal computers. However, very complex simulations are candidates for the use of high performance computing. In particular, applications stemming from the field of biomechanics can profit from the use of parallel computers (Hartmann et al., 2001). The complex geometry of anatomical structures often requires a very high spatial resolution of the model leading to possibly huge demand of computer power. In Figure 4 an automatic tool chain to generate *individual* models of human body parts is presented that takes advantage of the application of parallel computers. The software environment has been tested for impact analyses on the human head. It includes the complete chain of tools necessary for the geometric model generation from medical scan data (segmentation, mesh generation and mesh manipulation), computer simulation and visualisation. The software tools are public domain and are available on the web (Berti,2010).

Fig. 4. From medical scan data to finite elements.

Biomechanical Computer Models 101

Due to the complex anatomy and the high loading the knee joint is one of the most frequently injured joints of the human body. In addition to injuries like ruptures of the ligaments caused by accidents, chronic diseases as osteoarthritis affect a high amount of the population. Caused by the estimated demographic progression an increase of total knee joint replacements will be observed in the near future. FE models help to understand the behaviour of the structures of the human knee joint under dynamic loading. The strain on the cartilage and the load distribution effects of the menisci can show how the load is transferred through the different anatomical parts. Another line of application is the investigation of the mechanics of the main ligaments of the knee joint. The anterior and posterior cruciate ligaments (ACL, PCL) and the medial and lateral collateral ligaments (MCL, PCL) play an important role in stabilising the knee joint during daily activities and restrain the joint motion in nearly all degrees of freedom. The understanding of the combined role of the knee joint structures in stabilising and restraining the motion of the joint is fundamental for diagnostic procedures (Blankevoort & Huiskes, 1991). This knowledge can help to develop interventions such as training strategies or surgical procedures for the prevention of a knee joint disease. Also an optimization of orthopaedic devices is possible by the use of FE simulations. Even though many studies have been done, the kinematic and kinetic behaviour of the human knee joint is still not fully understood. Furthermore, the complex injury mechanisms and the development of degenerative joint injuries need additional research. Experimental studies often cannot satisfy this demand due to several limitations like for instance their high costs, the difficult measuring situation and the fact that the experiments can hardly be reproduced. As shown at the beginning of this section accurate FE models have proven to be a good alternative in determining knee joint mechanics. In order to get reliable and accurate results a good segmentation of the anatomical structures is highly important. Other quality influencing parameters are the mesh resolution and the choice of the finite elements. Most of the commercial meshing software products offer the possibility to create surface or volume meshes. The quality of the models has become much better within the last decade because of the high-performance computer chips that enable more complex and detailed models with smaller elements. In order to get an idea of the improvement of FE modelling related to biomechanics of the knee joint we will briefly

**4.2 A model for contact analysis in the tibio-femoral joint** 

describe some interesting models developed in the last twenty years.

A three dimensional model of the human knee joint was presented by Blankevoort et al., 1991, to investigate the articular contact of femur and tibia. The paper gave a good mathematical description of the modelling process and the material properties for the different structures. A more detailed model with 798 eight-node solid elements for the cartilage and the menisci, 1212 truss elements for the reinforcing of the solid elements was built to simulate the compression of the knee joint. Only 39 uniaxial elements represented the ligaments (Bendjaballah et al., 1995). Two years later the same authors restructured the model to study the biomechanical influence of varus-valgus (Bendjaballah et al., 1997) followed by an analysis of the tibio-femoral joint in axial rotation (Jilani et al., 1997). These studies are good examples for the advantages of FE models. After a detailed geometry reconstruction and a reliable validation process a lot of biomechanical problems can be explored with only few variations of the same model saving time and money. The availability of commercial software tools and faster multi-core computers lead to an

**4.2.1 Introduction** 

The first step is the identification of anatomical structures of interest in a 3D voxel dataset derived from medical scan data (CT or MRI). After this image processing step follows the geometric modelling of the structures identified in the previous step. The important difference to standard engineering problems is the lack of continuum geometry. The starting point is a discrete voxel dataset with a typical resolution of 1 mm or less. The chosen algorithm to generate a discrete representation of the object suitable for FE simulations is based on an octree data structure. It allows the fast generation of very large tetrahedral and hexahedral meshes. The scan data provide information about geometry only. The computer model is not complete without initial and boundary conditions and reliable information about material behaviour. The FE code for solving the underlying PDE is fully parallel and runs on a Linux cluster with distributed memory architecture. This enables the calculation of millions of unknowns.

Fig. 5. Evolution of a pressure wave in the brain after frontal impact*.*

Figure 5 shows an example of a high-resolution impact analysis based on the above described tool chain. This analysis investigated the mechanics of the *coup-contrecoup*, a phenomenon that is well known in neurology. After a heavy (often deadly) impact to the forehead the neurologist observes a brain contusion at the frontal brain and surprisingly enough a second contusion at the occipital brain lobe. The calculation that provides an FE solution for equation (5) shows that there are two distinct regions of pressure developing in the brain. In the frontal region an area of increased pressure is spreading out whilst on the opposite side an area of low pressure is starting to grow. Low and high-pressure zones are related to an increased risk of brain damage. Thus, the simulation confirms the clinical findings. Furthermore, the results have been validated against crash experiments proving the reliability of the FE model (Hartmann & Kruggel, 1998). However, for many applications the use of parallel computers is not necessary. But the pre-processing becomes much more sophisticated because the number of finite elements has to be controlled during the meshing process. Mesh regions of high physical interest have to be refined whereas other mesh areas might show a lower spatial resolution. This explains the demand for commercial meshing and simulation tools. In the following chapter a very complex model of the human knee joint is described. The model is developed on a typical desktop PC making use of commercial software.

The first step is the identification of anatomical structures of interest in a 3D voxel dataset derived from medical scan data (CT or MRI). After this image processing step follows the geometric modelling of the structures identified in the previous step. The important difference to standard engineering problems is the lack of continuum geometry. The starting point is a discrete voxel dataset with a typical resolution of 1 mm or less. The chosen algorithm to generate a discrete representation of the object suitable for FE simulations is based on an octree data structure. It allows the fast generation of very large tetrahedral and hexahedral meshes. The scan data provide information about geometry only. The computer model is not complete without initial and boundary conditions and reliable information about material behaviour. The FE code for solving the underlying PDE is fully parallel and runs on a Linux cluster with distributed memory architecture. This enables the calculation of

Fig. 5. Evolution of a pressure wave in the brain after frontal impact*.*

Figure 5 shows an example of a high-resolution impact analysis based on the above described tool chain. This analysis investigated the mechanics of the *coup-contrecoup*, a phenomenon that is well known in neurology. After a heavy (often deadly) impact to the forehead the neurologist observes a brain contusion at the frontal brain and surprisingly enough a second contusion at the occipital brain lobe. The calculation that provides an FE solution for equation (5) shows that there are two distinct regions of pressure developing in the brain. In the frontal region an area of increased pressure is spreading out whilst on the opposite side an area of low pressure is starting to grow. Low and high-pressure zones are related to an increased risk of brain damage. Thus, the simulation confirms the clinical findings. Furthermore, the results have been validated against crash experiments proving the reliability of the FE model (Hartmann & Kruggel, 1998). However, for many applications the use of parallel computers is not necessary. But the pre-processing becomes much more sophisticated because the number of finite elements has to be controlled during the meshing process. Mesh regions of high physical interest have to be refined whereas other mesh areas might show a lower spatial resolution. This explains the demand for commercial meshing and simulation tools. In the following chapter a very complex model of the human knee joint is described. The model is developed on a typical desktop PC making use of

millions of unknowns.

commercial software.

### **4.2 A model for contact analysis in the tibio-femoral joint 4.2.1 Introduction**

Due to the complex anatomy and the high loading the knee joint is one of the most frequently injured joints of the human body. In addition to injuries like ruptures of the ligaments caused by accidents, chronic diseases as osteoarthritis affect a high amount of the population. Caused by the estimated demographic progression an increase of total knee joint replacements will be observed in the near future. FE models help to understand the behaviour of the structures of the human knee joint under dynamic loading. The strain on the cartilage and the load distribution effects of the menisci can show how the load is transferred through the different anatomical parts. Another line of application is the investigation of the mechanics of the main ligaments of the knee joint. The anterior and posterior cruciate ligaments (ACL, PCL) and the medial and lateral collateral ligaments (MCL, PCL) play an important role in stabilising the knee joint during daily activities and restrain the joint motion in nearly all degrees of freedom. The understanding of the combined role of the knee joint structures in stabilising and restraining the motion of the joint is fundamental for diagnostic procedures (Blankevoort & Huiskes, 1991). This knowledge can help to develop interventions such as training strategies or surgical procedures for the prevention of a knee joint disease. Also an optimization of orthopaedic devices is possible by the use of FE simulations. Even though many studies have been done, the kinematic and kinetic behaviour of the human knee joint is still not fully understood. Furthermore, the complex injury mechanisms and the development of degenerative joint injuries need additional research. Experimental studies often cannot satisfy this demand due to several limitations like for instance their high costs, the difficult measuring situation and the fact that the experiments can hardly be reproduced. As shown at the beginning of this section accurate FE models have proven to be a good alternative in determining knee joint mechanics. In order to get reliable and accurate results a good segmentation of the anatomical structures is highly important. Other quality influencing parameters are the mesh resolution and the choice of the finite elements. Most of the commercial meshing software products offer the possibility to create surface or volume meshes. The quality of the models has become much better within the last decade because of the high-performance computer chips that enable more complex and detailed models with smaller elements. In order to get an idea of the improvement of FE modelling related to biomechanics of the knee joint we will briefly describe some interesting models developed in the last twenty years.

A three dimensional model of the human knee joint was presented by Blankevoort et al., 1991, to investigate the articular contact of femur and tibia. The paper gave a good mathematical description of the modelling process and the material properties for the different structures. A more detailed model with 798 eight-node solid elements for the cartilage and the menisci, 1212 truss elements for the reinforcing of the solid elements was built to simulate the compression of the knee joint. Only 39 uniaxial elements represented the ligaments (Bendjaballah et al., 1995). Two years later the same authors restructured the model to study the biomechanical influence of varus-valgus (Bendjaballah et al., 1997) followed by an analysis of the tibio-femoral joint in axial rotation (Jilani et al., 1997). These studies are good examples for the advantages of FE models. After a detailed geometry reconstruction and a reliable validation process a lot of biomechanical problems can be explored with only few variations of the same model saving time and money. The availability of commercial software tools and faster multi-core computers lead to an

Biomechanical Computer Models 103

The segmentation of the images is done to define anatomical structures in the threedimensional image stack. This image processing is an important step in the workflow and crucial for the quality of the simulation results. Segmentation can be done with respect to different kinds of concepts; but the output remains the same. A range of grey values characterise the different structures in an image. Different automatic segmentation principles are based on pixel grey values, region or edge detection. These can be used for images in which the anatomical structures show a huge difference in the grey values. For images of the knee joint obtained by MRI a full automatic segmentation is *not* possible due to the fact that the differences of grey values are not high enough. Therefore it is recommended to use a semiautomatic segmentation, which starts with one of the automatic segmentation principles for detection of large connected areas. After this pre- segmentation an interactive improvement of the segmentation by hand leads to a clearly defined differentiation of the structures. A short overview of the numerous automatic methods and the workflow of the process will be shown in the following. The threshold operation is one of the best-known automatic detection procedures. It belongs to the pixel-based procedures, is fast but not very accurate for homogeneous images. It is often used for fast overview segmentations, if no detailed information is necessary. The algorithm compares each pixel grey value with a user-defined threshold; the pixel is assigned to the structure if the pixel value is in the range of the threshold. The application on the whole image is called global thresholding which can often lead to an excluding of pixels that belong to the same structure. Better results can be obtained by using a local thresholding, where the image is separated into different areas with varying thresholds. As mentioned in the beginning it is not sensible enough to detect small changes in an image, but can be applied on regions with nearly the same grey values like the bones. Another method to identify different regions in an image is based on edge detection**.** The principle of operation can be compared to the threshold method. Instead of defining a specific value the algorithm checks for grey value changes in an image by comparing two neighbouring pixels. Is there a clear difference of their grey values a part of the edge is detected and two structures are found. The combining of all grey values differences provides a line between the detected parts and all pixels of one region are then merged. In the field of image processing the Sobel-operator, the Laplace-operator or the gradient search are well known filters to find edges in an image. The detected edges have to be combined by special edge tracing algorithms. Like with the threshold method the edge detection algorithm is often not sensitive enough for a satisfying segmentation, so that a manual improvement is applied. Figure 7 illustrates the mentioned problems occurring with the automatic segmentation algorithms. As a simple example a pure automatic segmentation of the femur is shown on the left (a). It can be seen that it was not possible to define the threshold such that only the bony parts are detected. The expected edge is not straight-lined; in regions with similar grey values protuberances can be found and the inner part is not completely detected. Without a manual correction of the border it is not possible to detect a clear defined region by only

Therefore different tools are available to mark areas, which are incorrectly detected or

In this example the automatic detection can be applied to all 100 slices at the same time, but this also increases the number of errors because the grey values of a specified structure differ from image to image in the image stack. For the above-mentioned reasons it is only applicable for the large bony structures. The soft tissues like the meniscus or the ligaments

omitted. After the whole structure is marked it is added to the specified part.

**4.2.3 Segmentation process** 

defining a threshold.

increased number of models and application fields such as prosthesis development (Godest et al., 2002; Soncini et al., 2002).

### **4.2.2 Image acquisition**

The reconstruction process of the individual geometry of the investigated object is based on medical image data obtained by Computed Tomography (CT) or Magnetic Resonance Imaging (MRI). Both methods create two-dimensional images with a predefined slice thickness and slice distance of a specified region. A set of these slices is called image stack and represents a three-dimensional recording of the knee joint. The advantage of using MRI is the possibility to visualise muscles and ligaments in addition to the bony structures whereas CT gives a better representation of the bony structures. It is obvious that the choice of the imaging procedure is depending on the specific research question.

Fig. 6. Image stack of the knee joint from MRI, b) different orientations of the images

In the study presented it is important to have a nearly full representation of the knee joint including the articulating bones (femur, tibia and patella) and also the stabilizing ligaments and the meniscus. The images of a healthy female knee joint are taken in a Philips Gyroscan Intera MRI with a field intensity of 1.0 Tesla.

We use the software Amira (Version 4.0, Mercury Systems) for the data visualisation and segmentation. The program offers the option to apply several image processing tools such as filters for a better differentiation of the structural edges resulting in a more accurate segmentation. In order to get a first and fast overview of the joint geometry a volume rendering depending on the grey values can be done. From a clinical point of view the deviations of the grey values within a structure can show possible pathological changes (for example in the cartilage or the bones).

increased number of models and application fields such as prosthesis development (Godest

The reconstruction process of the individual geometry of the investigated object is based on medical image data obtained by Computed Tomography (CT) or Magnetic Resonance Imaging (MRI). Both methods create two-dimensional images with a predefined slice thickness and slice distance of a specified region. A set of these slices is called image stack and represents a three-dimensional recording of the knee joint. The advantage of using MRI is the possibility to visualise muscles and ligaments in addition to the bony structures whereas CT gives a better representation of the bony structures. It is obvious that the choice of the imaging procedure is depending on the specific research

Fig. 6. Image stack of the knee joint from MRI, b) different orientations of the images

Intera MRI with a field intensity of 1.0 Tesla.

**(a) (b)**

example in the cartilage or the bones).

In the study presented it is important to have a nearly full representation of the knee joint including the articulating bones (femur, tibia and patella) and also the stabilizing ligaments and the meniscus. The images of a healthy female knee joint are taken in a Philips Gyroscan

We use the software Amira (Version 4.0, Mercury Systems) for the data visualisation and segmentation. The program offers the option to apply several image processing tools such as filters for a better differentiation of the structural edges resulting in a more accurate segmentation. In order to get a first and fast overview of the joint geometry a volume rendering depending on the grey values can be done. From a clinical point of view the deviations of the grey values within a structure can show possible pathological changes (for

et al., 2002; Soncini et al., 2002).

**4.2.2 Image acquisition** 

question.

### **4.2.3 Segmentation process**

The segmentation of the images is done to define anatomical structures in the threedimensional image stack. This image processing is an important step in the workflow and crucial for the quality of the simulation results. Segmentation can be done with respect to different kinds of concepts; but the output remains the same. A range of grey values characterise the different structures in an image. Different automatic segmentation principles are based on pixel grey values, region or edge detection. These can be used for images in which the anatomical structures show a huge difference in the grey values. For images of the knee joint obtained by MRI a full automatic segmentation is *not* possible due to the fact that the differences of grey values are not high enough. Therefore it is recommended to use a semiautomatic segmentation, which starts with one of the automatic segmentation principles for detection of large connected areas. After this pre- segmentation an interactive improvement of the segmentation by hand leads to a clearly defined differentiation of the structures. A short overview of the numerous automatic methods and the workflow of the process will be shown in the following. The threshold operation is one of the best-known automatic detection procedures. It belongs to the pixel-based procedures, is fast but not very accurate for homogeneous images. It is often used for fast overview segmentations, if no detailed information is necessary. The algorithm compares each pixel grey value with a user-defined threshold; the pixel is assigned to the structure if the pixel value is in the range of the threshold. The application on the whole image is called global thresholding which can often lead to an excluding of pixels that belong to the same structure. Better results can be obtained by using a local thresholding, where the image is separated into different areas with varying thresholds. As mentioned in the beginning it is not sensible enough to detect small changes in an image, but can be applied on regions with nearly the same grey values like the bones.

Another method to identify different regions in an image is based on edge detection**.** The principle of operation can be compared to the threshold method. Instead of defining a specific value the algorithm checks for grey value changes in an image by comparing two neighbouring pixels. Is there a clear difference of their grey values a part of the edge is detected and two structures are found. The combining of all grey values differences provides a line between the detected parts and all pixels of one region are then merged. In the field of image processing the Sobel-operator, the Laplace-operator or the gradient search are well known filters to find edges in an image. The detected edges have to be combined by special edge tracing algorithms. Like with the threshold method the edge detection algorithm is often not sensitive enough for a satisfying segmentation, so that a manual improvement is applied.

Figure 7 illustrates the mentioned problems occurring with the automatic segmentation algorithms. As a simple example a pure automatic segmentation of the femur is shown on the left (a). It can be seen that it was not possible to define the threshold such that only the bony parts are detected. The expected edge is not straight-lined; in regions with similar grey values protuberances can be found and the inner part is not completely detected. Without a manual correction of the border it is not possible to detect a clear defined region by only defining a threshold.

Therefore different tools are available to mark areas, which are incorrectly detected or omitted. After the whole structure is marked it is added to the specified part.

In this example the automatic detection can be applied to all 100 slices at the same time, but this also increases the number of errors because the grey values of a specified structure differ from image to image in the image stack. For the above-mentioned reasons it is only applicable for the large bony structures. The soft tissues like the meniscus or the ligaments

Biomechanical Computer Models 105

Due to the image acquisition procedure a specified slice thickness and distance is given causing errors in the reconstruction process. The gap between each slice does not provide any information about the shape of the anatomic structures. This leads to stairs in the 3D reconstruction in areas where the alteration of the structures is smaller then the slice distance. It is obvious that these stairs adulterate not only the design but also the simulation results. Therefore it is mandatory to make the stairs as small as possible by interpolating the

Fig. 9. a) Reconstructed femur with existing stairs due to the slice distance, b) surface of the

In the interpolation process imaginary points, calculated with the existing points of the segmentation are included between two slices thus minimising the effect of the stairs. This can either be done by a linear or a cubic interpolation. Following the interpolation step the surface is reconstructed from the segmented structures. This is a pre-processing step before the real meshing procedure where the surface will be divided into the small tetrahedral elements. Directly after the meshing process a smoothing filter is applied. A comparison between the femur condyles after the surface generation with existing stairs and after the

The reconstructed three-dimensional surface of the anatomic structures is transformed into a surface mesh consisting of triangles using a Marching Cube algorithm before an Advancing Front algorithm constructs the volume mesh. For the surface triangulation Amira uses a modified type of the Marching Cube algorithm developed by Lorensen and Cline [20]. The algorithm assumes a lattice model of a volumetric dataset with a specific value at each intercept point of the lattice. The space is subdivided into a constant grid consisting of cubes, where each cube has eight pixels and eight grey values from two successive slices. Four pixels belong to a slice. In the case that the vertices of a cube's edge are lying on

femur after the interpolation and smoothing process

(a) (b)

interpolation and smoothing is shown in Figure 9.

**4.2.4 Meshing of the segmented images** 

segmented structures and smoothing the surface, before the meshing begins.

cannot be detected automatically. The homogeneous grey value parts of the soft tissue in the middle of the knee joint and the partly thin cartilage are completely manually segmented with a live wire procedure. In total 11 structures in the 100 slices are segmented. This is needed in the following meshing procedure and for the subsequent calculations of the FE solver. An example of one segmented slice is shown in Figure 8.

Fig. 7. a) Automatic segmentation of the femur with protuberances, b) segmented femur after improvement

Fig. 8. Fully segmented slice of the knee joint.

cannot be detected automatically. The homogeneous grey value parts of the soft tissue in the middle of the knee joint and the partly thin cartilage are completely manually segmented with a live wire procedure. In total 11 structures in the 100 slices are segmented. This is needed in the following meshing procedure and for the subsequent calculations of the FE

Fig. 7. a) Automatic segmentation of the femur with protuberances, b) segmented femur

(b )

solver. An example of one segmented slice is shown in Figure 8.

after improvement

(a )

Fig. 8. Fully segmented slice of the knee joint.

Due to the image acquisition procedure a specified slice thickness and distance is given causing errors in the reconstruction process. The gap between each slice does not provide any information about the shape of the anatomic structures. This leads to stairs in the 3D reconstruction in areas where the alteration of the structures is smaller then the slice distance. It is obvious that these stairs adulterate not only the design but also the simulation results. Therefore it is mandatory to make the stairs as small as possible by interpolating the segmented structures and smoothing the surface, before the meshing begins.

Fig. 9. a) Reconstructed femur with existing stairs due to the slice distance, b) surface of the femur after the interpolation and smoothing process

In the interpolation process imaginary points, calculated with the existing points of the segmentation are included between two slices thus minimising the effect of the stairs. This can either be done by a linear or a cubic interpolation. Following the interpolation step the surface is reconstructed from the segmented structures. This is a pre-processing step before the real meshing procedure where the surface will be divided into the small tetrahedral elements. Directly after the meshing process a smoothing filter is applied. A comparison between the femur condyles after the surface generation with existing stairs and after the interpolation and smoothing is shown in Figure 9.

### **4.2.4 Meshing of the segmented images**

The reconstructed three-dimensional surface of the anatomic structures is transformed into a surface mesh consisting of triangles using a Marching Cube algorithm before an Advancing Front algorithm constructs the volume mesh. For the surface triangulation Amira uses a modified type of the Marching Cube algorithm developed by Lorensen and Cline [20]. The algorithm assumes a lattice model of a volumetric dataset with a specific value at each intercept point of the lattice. The space is subdivided into a constant grid consisting of cubes, where each cube has eight pixels and eight grey values from two successive slices. Four pixels belong to a slice. In the case that the vertices of a cube's edge are lying on

Biomechanical Computer Models 107

of the results. It is obvious that a compromise between size and amount of elements and the computation time has to be found. After the simplification step the mesh is smoothed to eliminate the mentioned stairs in case they still exist. During this second smoothing the edge points of the elements are moved. For the calculation of the new position, the middle point of the neighbour elements is determined, and the point to be moved is shifted into the new direction. This is done for all elements, which do not have a uniform changeover from one element to another. Is the form of the mesh satisfying and close to the original structures, its quality has to be determined through several tests. If the mesh quality is not satisfactory, an additional smoothing can be helpful. Checking the dihedral angle, the orientation and the aspect ratio of the elements provides a good measure of the mesh quality. These tests should be done consecutively, as well as the intersection and closeness test. If a problem or a bad triangle is found, the element is highlighted and can be improved using the editor by shifting, subdividing or drawing. If no problems are present the generation of the volume mesh can be started. Based on the surface the mesh of the inner volume is constructed. The described tests for the surface should also be done for the volume mesh. For special regions of interest, e.g. the meniscus or the articulating surfaces of femur and tibia, a refinement of the mesh is taken into account. The existing elements are bisected, so that a finer and smoother structure is achieved. This increases not only the number of nodes and elements but also the quality of the results of the simulation. The complete mesh of the 3D

Madymo (Mathematic dynamic modelling, Version 6.3, Tass–safe, Netherlands) is used for the FE analyses based on the meshes depicted in Figure 10. Madymo enables a combination of a multi-body and finite element parts. Models are built in a kinematic chain where the rigid multi-body elements are connected by joints. The mesh obtained by Amira is attached to those bodies thus ensuring that a prescribed joint motion results in a movement of the bones. The several steps of the implementation in Madymo are described in the following. After having chosen a reference frame three bodies are defined; one for the femur, one for the tibia and one for the patella. To each of these bodies, the elements and nodes belonging to the specified bony structure are assigned. The described x, y and z coordinates of the nodes, the part numbers and the elements for each part, which were exported from Amira, are imported in two different tables. The part numbers are crucial for the assignment to the bodies and for the differentiation of the material properties. Free joints, which have three translational and three rotational degrees of freedom are used for the connection of the bodies. To guarantee a realistic kinematic and kinetic behaviour of the model, the physical properties like the moment of inertia and the mass with its centre are essential. The determination of these parameters is a complex physical procedure and individual for each subject. For this reasons it is common to take values from literature, which were gained, in cross-sectional studies. The inertial characteristics of several human segments can be found

The material properties are important for a realistic model and essential for the stress calculation. In the material sciences there are different typical laws known to describe a material's behaviour with mathematical formulae. Such mathematical material models are also used in finite element simulations to compute the resulting stress from the displacement field of the nodes. The properties of a biological material are always an approximation

reconstruction of the human knee joint is shown in Figure 10.

**4.2.5 Implementation of the model** 

in Zatsiorsky, 2002.

different sides of the isosurface, it is given that this edge is cut by the surface. The detected intersection points of the cube and the isosurface are combined so that a polygon and a triangle are created. Theoretically there are 256 different possible combinations given how the cut goes through the surface. Due to symmetries these can be reduced to 14 classes. The summary of all polygons leads to the 3D reconstruction of the surface. This surface is used to create a volumetric mesh using the above-mentioned Advancing Front algorithm. With this it is possible to create elements and nodes at the same time. The edge is set by nodes, which are divided into connections lines. The length of such a line is depending on the specified size of the elements and results in a closed polygon, the "advancing front". Along this front the elements are created, which should go into the inner part of the area to be meshed. An update of this front includes only the connection lines, which are needed for subsequent construction of the mesh. This method is properly suited for complex geometries with inner edges because the boundary of the area is done first. Both algorithms are often used in commercial meshing software and described in detail by Nigel et al., 1999.

Fig. 10. Mesh of the human knee joint consisting of tetrahedral elements.

Before starting the surface generation the minimal edge length of the triangles and the choice for a first smoothing of the surface have to be set. The result is a closed surface mesh with rather small elements, which leads to a high amount of nodes and elements. The more elements exist, the better is the reconstruction of the surface but also the computation times increases extremely. Therefore it is recommended to reduce the amount of elements using the simplification editor. The positions of the elements are recalculated, so that the same surface is constructed using fewer elements. The meshing procedure is highly important for the whole FE analysis as the quality of the mesh influences the robustness and the accuracy

different sides of the isosurface, it is given that this edge is cut by the surface. The detected intersection points of the cube and the isosurface are combined so that a polygon and a triangle are created. Theoretically there are 256 different possible combinations given how the cut goes through the surface. Due to symmetries these can be reduced to 14 classes. The summary of all polygons leads to the 3D reconstruction of the surface. This surface is used to create a volumetric mesh using the above-mentioned Advancing Front algorithm. With this it is possible to create elements and nodes at the same time. The edge is set by nodes, which are divided into connections lines. The length of such a line is depending on the specified size of the elements and results in a closed polygon, the "advancing front". Along this front the elements are created, which should go into the inner part of the area to be meshed. An update of this front includes only the connection lines, which are needed for subsequent construction of the mesh. This method is properly suited for complex geometries with inner edges because the boundary of the area is done first. Both algorithms are often used in commercial meshing software and described in detail by Nigel et al., 1999.

Fig. 10. Mesh of the human knee joint consisting of tetrahedral elements.

Before starting the surface generation the minimal edge length of the triangles and the choice for a first smoothing of the surface have to be set. The result is a closed surface mesh with rather small elements, which leads to a high amount of nodes and elements. The more elements exist, the better is the reconstruction of the surface but also the computation times increases extremely. Therefore it is recommended to reduce the amount of elements using the simplification editor. The positions of the elements are recalculated, so that the same surface is constructed using fewer elements. The meshing procedure is highly important for the whole FE analysis as the quality of the mesh influences the robustness and the accuracy of the results. It is obvious that a compromise between size and amount of elements and the computation time has to be found. After the simplification step the mesh is smoothed to eliminate the mentioned stairs in case they still exist. During this second smoothing the edge points of the elements are moved. For the calculation of the new position, the middle point of the neighbour elements is determined, and the point to be moved is shifted into the new direction. This is done for all elements, which do not have a uniform changeover from one element to another. Is the form of the mesh satisfying and close to the original structures, its quality has to be determined through several tests. If the mesh quality is not satisfactory, an additional smoothing can be helpful. Checking the dihedral angle, the orientation and the aspect ratio of the elements provides a good measure of the mesh quality. These tests should be done consecutively, as well as the intersection and closeness test. If a problem or a bad triangle is found, the element is highlighted and can be improved using the editor by shifting, subdividing or drawing. If no problems are present the generation of the volume mesh can be started. Based on the surface the mesh of the inner volume is constructed. The described tests for the surface should also be done for the volume mesh. For special regions of interest, e.g. the meniscus or the articulating surfaces of femur and tibia, a refinement of the mesh is taken into account. The existing elements are bisected, so that a finer and smoother structure is achieved. This increases not only the number of nodes and elements but also the quality of the results of the simulation. The complete mesh of the 3D

### **4.2.5 Implementation of the model**

reconstruction of the human knee joint is shown in Figure 10.

Madymo (Mathematic dynamic modelling, Version 6.3, Tass–safe, Netherlands) is used for the FE analyses based on the meshes depicted in Figure 10. Madymo enables a combination of a multi-body and finite element parts. Models are built in a kinematic chain where the rigid multi-body elements are connected by joints. The mesh obtained by Amira is attached to those bodies thus ensuring that a prescribed joint motion results in a movement of the bones. The several steps of the implementation in Madymo are described in the following.

After having chosen a reference frame three bodies are defined; one for the femur, one for the tibia and one for the patella. To each of these bodies, the elements and nodes belonging to the specified bony structure are assigned. The described x, y and z coordinates of the nodes, the part numbers and the elements for each part, which were exported from Amira, are imported in two different tables. The part numbers are crucial for the assignment to the bodies and for the differentiation of the material properties. Free joints, which have three translational and three rotational degrees of freedom are used for the connection of the bodies. To guarantee a realistic kinematic and kinetic behaviour of the model, the physical properties like the moment of inertia and the mass with its centre are essential. The determination of these parameters is a complex physical procedure and individual for each subject. For this reasons it is common to take values from literature, which were gained, in cross-sectional studies. The inertial characteristics of several human segments can be found in Zatsiorsky, 2002.

The material properties are important for a realistic model and essential for the stress calculation. In the material sciences there are different typical laws known to describe a material's behaviour with mathematical formulae. Such mathematical material models are also used in finite element simulations to compute the resulting stress from the displacement field of the nodes. The properties of a biological material are always an approximation

Biomechanical Computer Models 109

shapes of the femur and the menisci. The gap can either be set as a function or is defined by program. The contact edge search algorithm can find additional contacts, which could be lost in presence of the curved shapes. The penalty method is used for the calculation of contact forces, with a damping coefficient of 0.05 and a variable FE time step, given that the

For the kinematic boundary conditions a simple knee bending is used for the prescribed motion, which is obtained in a three-dimensional movement analysis using retro reflective markers and infrared cameras (Vicon Nexus, Vicon Motion System, Oxford, UK). The calculated joint angles in three directions are given as an angle vs. time function. The kinetic boundary conditions are applied by a gravity function which includes the mass of the

The goal of the simulation is the calculation of the contact forces and stresses in three directions as well as the resultants. A distribution of the forces and stresses acting on the articulating surfaces is shown with a specified colour bar referring to the calculated values. The results provide a better understanding of the transduction of applied outer forces in the inner structures of the human knee joint. In addition different scenarios like meniscal tears or ligament ruptures can easily be simulated and the biomechanical effects on the joint investigated. The influence of joint alignment like a varus or valgus condition is of great interest in these simulations. In Figure 11 the results of the knee bending simulation under body weight condition are shown. The menisci and the cruciate ligaments are removed for a better illustration of the forces on the articulating cartilage. The distribution shows a peak in

contact force is dependent on the time step, to keep the simulation stable.

the direct contact area and a decrease from the middle to the border.

Fig. 11. Calculated contact forces of the simulation of a knee bending.

investigated subject.

**4.2.6 First results** 

because they vary due to different influencing factors like age or training. For this reason simplifications in expressing the material laws are often made. The isotropic elastic material is one of the easiest to describe because only two parameters are needed. Elastic means that a force-loaded material returns completely to its initial shape after it is unloaded. Isotropic indicates, the loading direction is of no matter, so there is no distinct direction. The describing parameters for an isotropic elastic material are Young's modulus **E** (in Mega Pascal, MPa) and the Poisson ratio ν (without a unit). When expecting small strains, it is common to use these two parameters for the material description because of the low computation costs. Other biomechanical models use a hyperelastic material law for the description, due to the fact, that the results obtained by a linearised model are not applicable for finite strain calculations.

The response of biological soft tissues of the human body is often time or load dependent. This means, that these viscoelastic materials show creep or relaxation effects. Applying a constant force to a material produces a lengthening. Creeping is given if a material shows an enlarging displacement in addition to the normal response to the outer force. Relaxation of a material stands for the reduction of the inner stresses if a constant force is applied. For the human knee joint the different material responses mentioned above are used to specify the biological tissues. All properties are taken from the literature. The femur, the tibia and the patella are assumed to be rigid with a density of 1000 kg/m3. For the cartilage of the bones viscoelastic behaviour was neglected due to a short loading time. Therefore the assumption of a linear elastic and isotropic material was made and a Young's modulus of E = 5 MPa and a Poisson ratio of ν = 0.46 were chosen (Li et al., 2001). The menisci are hydrated tissues. They are modelled as linear elastic for the same reason that holds true for the cartilage but with a much higher Young's modulus of E = 59 MPa and a Poisson ratio of ν = 0.49 (LeRoux & Setton, 2002). The ligament properties are defined by a stress strain function (Lee & Hyman, 2002).

The articulating surfaces of the human knee joint show a complex geometry with curved shapes in all directions. Due to this fact it is important to have a look on the contact modelling that will influence the computed strains. During a movement of the femur, the menisci and the cartilage parts of the bones are in contact with each other. Therefore a contact definition for each contact pair is defined, the result of a contact between two surfaces is a deformation of at least one of these. Also a self-contact of one structure is possible. Three algorithms are well known for the description of a finite element contact: Node-to surface, point-to surface and surface-to surface.

In the node-to surface contact the nodes of structure 1 are in contact to the surface of the elements of structure 2. The time consuming tracking of contacts between nodes and the surface is simplified with special search strategies. For a three-node element, all nodes can have contact to the elements of the penetrated surface. If elements have an additional point in the middle of the element, the point-to-surface method should be used, which also looks for other points that can apply forces due to penetration. The surface-to-surface method is mainly used for the contact between rigid bodies with a special penetration characteristic and for heavily curved shapes.

The first step is to make a decision, which part is *Master* and which one is *Slave.* For the case that the mesh has different element sizes due to a refinement of a structure, the finer mesh should be chosen as *Slave*. Is this not the case the surface with the more curved shape should be the *Slave*. In this model, the surface-to-surface method is chosen, because of the curved

because they vary due to different influencing factors like age or training. For this reason simplifications in expressing the material laws are often made. The isotropic elastic material is one of the easiest to describe because only two parameters are needed. Elastic means that a force-loaded material returns completely to its initial shape after it is unloaded. Isotropic indicates, the loading direction is of no matter, so there is no distinct direction. The describing parameters for an isotropic elastic material are Young's modulus **E** (in Mega Pascal, MPa) and the Poisson ratio ν (without a unit). When expecting small strains, it is common to use these two parameters for the material description because of the low computation costs. Other biomechanical models use a hyperelastic material law for the description, due to the fact, that the results obtained by a linearised model are not applicable

The response of biological soft tissues of the human body is often time or load dependent. This means, that these viscoelastic materials show creep or relaxation effects. Applying a constant force to a material produces a lengthening. Creeping is given if a material shows an enlarging displacement in addition to the normal response to the outer force. Relaxation of a material stands for the reduction of the inner stresses if a constant force is applied. For the human knee joint the different material responses mentioned above are used to specify the biological tissues. All properties are taken from the literature. The femur, the tibia and the patella are assumed to be rigid with a density of 1000 kg/m3. For the cartilage of the bones viscoelastic behaviour was neglected due to a short loading time. Therefore the assumption of a linear elastic and isotropic material was made and a Young's modulus of E = 5 MPa and a Poisson ratio of ν = 0.46 were chosen (Li et al., 2001). The menisci are hydrated tissues. They are modelled as linear elastic for the same reason that holds true for the cartilage but with a much higher Young's modulus of E = 59 MPa and a Poisson ratio of ν = 0.49 (LeRoux & Setton, 2002). The ligament properties are defined by a stress strain function (Lee &

The articulating surfaces of the human knee joint show a complex geometry with curved shapes in all directions. Due to this fact it is important to have a look on the contact modelling that will influence the computed strains. During a movement of the femur, the menisci and the cartilage parts of the bones are in contact with each other. Therefore a contact definition for each contact pair is defined, the result of a contact between two surfaces is a deformation of at least one of these. Also a self-contact of one structure is possible. Three algorithms are well known for the description of a finite element contact:

In the node-to surface contact the nodes of structure 1 are in contact to the surface of the elements of structure 2. The time consuming tracking of contacts between nodes and the surface is simplified with special search strategies. For a three-node element, all nodes can have contact to the elements of the penetrated surface. If elements have an additional point in the middle of the element, the point-to-surface method should be used, which also looks for other points that can apply forces due to penetration. The surface-to-surface method is mainly used for the contact between rigid bodies with a special penetration characteristic

The first step is to make a decision, which part is *Master* and which one is *Slave.* For the case that the mesh has different element sizes due to a refinement of a structure, the finer mesh should be chosen as *Slave*. Is this not the case the surface with the more curved shape should be the *Slave*. In this model, the surface-to-surface method is chosen, because of the curved

Node-to surface, point-to surface and surface-to surface.

and for heavily curved shapes.

for finite strain calculations.

Hyman, 2002).

shapes of the femur and the menisci. The gap can either be set as a function or is defined by program. The contact edge search algorithm can find additional contacts, which could be lost in presence of the curved shapes. The penalty method is used for the calculation of contact forces, with a damping coefficient of 0.05 and a variable FE time step, given that the contact force is dependent on the time step, to keep the simulation stable.

For the kinematic boundary conditions a simple knee bending is used for the prescribed motion, which is obtained in a three-dimensional movement analysis using retro reflective markers and infrared cameras (Vicon Nexus, Vicon Motion System, Oxford, UK). The calculated joint angles in three directions are given as an angle vs. time function. The kinetic boundary conditions are applied by a gravity function which includes the mass of the investigated subject.

### **4.2.6 First results**

The goal of the simulation is the calculation of the contact forces and stresses in three directions as well as the resultants. A distribution of the forces and stresses acting on the articulating surfaces is shown with a specified colour bar referring to the calculated values. The results provide a better understanding of the transduction of applied outer forces in the inner structures of the human knee joint. In addition different scenarios like meniscal tears or ligament ruptures can easily be simulated and the biomechanical effects on the joint investigated. The influence of joint alignment like a varus or valgus condition is of great interest in these simulations. In Figure 11 the results of the knee bending simulation under body weight condition are shown. The menisci and the cruciate ligaments are removed for a better illustration of the forces on the articulating cartilage. The distribution shows a peak in the direct contact area and a decrease from the middle to the border.

Fig. 11. Calculated contact forces of the simulation of a knee bending.

Biomechanical Computer Models 111

Blajer, W. & Czaplicki, A. (2001). Modelling and inverse simulation of somersaults on the

Blankevoort, L. & Huiskes, R. (1991). Ligament-bone interaction in a three-dimensional model of the knee. *Journal of Biomechanical Engineering* 113 (3), pp. 263-269 Blankevoort, L., Kuiper, J. H., Huiskes, R. & Grootenboer, H. J. (1991) Articular contact in a three-dimensional model of the knee. *Journal of Biomechanics* 24 (11), pp. 1019-1031 Brekelmans, W., Poort, H. & Slooff, T. (1972). A new method to analyse the mechanical behaviour of skeletal parts. *Acta orthopaedica Scandinavica*, 43(5), pp. 301-317 Coveney P. (2010) Virtual Physiological Human Network of Excellence, 04/09/2011,

Dabnichki, P. & Avital, E. (2006). Influence of the position of crewmembers on

Godest, A., Beaugonin, M., Haug, E., Taylor, M., & Gregson, P. (2002). Simulation of a knee

Gruber, K., Ruder, H., Denoth, J. & Schneider, K. (1998). A comparative study of impact

Hartmann, U. & Kruggel, F. (1998). Transient Analysis of the Biomechanics of the Human

Hartmann, U., Kruggel, F., Hierl, T., Lonsdale, G. & Kloeppel, R. (2001) Skull mechanics

Lee, M., & Hyman, W. (2002). Modeling of failure mode in knee ligaments depending on the

LeRoux, M. A., Setton. L. A. (2002). Experimental and biphasic FEM determinations of the

Li, G., Lopez, O., Rubash, H. (2001). Variability of a three-dimensional finite element model

Lorensen, W. & Cline, H. (1987) Marching Cubes: A high resolution 3D surface

Penrose, J.M.T. & Hose, D.R. (1998). Finite Element Impact Analysis of a Flexible Cricket Bat

Soncini, M., Redaelli, A., & Montevecchi, F., (2002). FE Dynamic Analysis Of A Knee Joint

analysis. *ASME Journal of Biomechanical Engineering* 123, pp. 341-346

reconstruction algorithm. *Computer Graphics*, 21 (4), pp 163-169

material properties and hydraulic perme- ability of the meniscus in tension. *ASME* 

constructed using magnetic resonance images of a knee for joint contact stress

for Design Optimisation. In: The Engineering of Sport - Design and Development,

Prosthesis During The Walking Cycle. *Proceeding of 2003 Summer Bioengineering* 

*Conference on Computational Solid and Fluid Mechanics* Boston May 2001 Huiskes, R. & Chao, E. (1983). A survey of finite element analysis in orthopaedic biomechanics: the first decade. *Journal of Biomechanics* 16 (6), pp. 385-409 Jilani, A., Shirazi-Adl, A. & Bendjaballah, M. (1997). Biomechanics of human tibio-femoral

*Biomechanics and Biomedical Engineering* 2(1), pp. 49-64

joint in axial rotation. *The Knee* 4, pp. 203-213

strain rate. *BMC musculoskeletal disorders* 3, 3, pp.

*Journal of Biomechanical Engineering* 124, pp. 315-321

S. Haake (ed.), pp. 531-539, Blackwell Science Ltd.

*Conference*, Key Biscayne, Florida, June 2003

aerodynamics performance of two-man bobsleigh. *Journal of Biomechanics* 39(15) ,

joint replacement during a gait cycle using explicit finite element analysis. *Journal of* 

dynamics: wobbling mass model versus rigid body models. *Journal of Biomechanics*

Head with a High-Resolution 3D Finite Element Model. *Computer Methods in* 

simulations with the prototype SimBio environment. *In: Proceedings of the M.I.T.* 

trampoline. *Journal of Biomechanics* 34, pp. 1619-1629

Available from http://www.vph-noe.eu/

pp. 2733-2742

31, pp. 439-444

*Biomechanics* 35 (2), 267-275

Having the option to prescribe motion enables the determination of the influence on the motion of the inner structures, like shifting of the menisci or the stretching of the ligaments. Using simulated muscle forces as input for the simulation the resulting locomotion and the acting forces on the articulating surfaces can be studied.

## **5. Conclusion**

The increasing availability of computational resources will lead to more detailed and complex models of the human body. In the future, efforts like the Virtual Physiological Human (Coveney, 2010) will provide a framework for developing and coupling biomechanical and bioelectric models of different scales and body parts. The term Virtual Physiological Human (VPH) indicates a shared resource formed by a federation of disparate but integrated computer models of the mechanical, physical, and biochemical functions of living human body in both physiological and pathological conditions. VPH models will be both descriptive and predictive. They comprise


The Virtual Physiological Human is a methodological and technological framework that once established, will enable collaborative investigation of the human body as a single complex system. VPH is not 'the supermodel' that will explain all possible aspects of human physiology or pathology. It is away to share observations, to derive predictive hypotheses from them, and to integrate them into a constantly improving understanding of human physiology/pathology, by regarding it as an integrated system.

## **6. Acknowledgements**

The authors would like to thank Prof. Dr. Karin Gruber (University of Koblenz, Germany) for providing valuable input related to MB models and motion analysis. Many thanks go as well to J. Seybold (University of Stuttgart) for making available his computer studies on the reaction forces for soft and hard landing. Special thanks go to Dr. Torsten Hans (University of Tübingen, Germany) for the visualisation of the jumping gymnast.

## **7. References**


Having the option to prescribe motion enables the determination of the influence on the motion of the inner structures, like shifting of the menisci or the stretching of the ligaments. Using simulated muscle forces as input for the simulation the resulting locomotion and the

The increasing availability of computational resources will lead to more detailed and complex models of the human body. In the future, efforts like the Virtual Physiological Human (Coveney, 2010) will provide a framework for developing and coupling biomechanical and bioelectric models of different scales and body parts. The term Virtual Physiological Human (VPH) indicates a shared resource formed by a federation of disparate but integrated computer models of the mechanical, physical, and biochemical functions of living human body in both physiological and pathological conditions. VPH models will be

• large collections of anatomical, physiological, and pathological data stored in digital

• services aimed to support the researchers in the creation and maintenance of these

• services aimed to empower clinical, industrial and societal users in the use of the VPH

The Virtual Physiological Human is a methodological and technological framework that once established, will enable collaborative investigation of the human body as a single complex system. VPH is not 'the supermodel' that will explain all possible aspects of human physiology or pathology. It is away to share observations, to derive predictive hypotheses from them, and to integrate them into a constantly improving understanding of human

The authors would like to thank Prof. Dr. Karin Gruber (University of Koblenz, Germany) for providing valuable input related to MB models and motion analysis. Many thanks go as well to J. Seybold (University of Stuttgart) for making available his computer studies on the reaction forces for soft and hard landing. Special thanks go to Dr. Torsten Hans (University

Bendjaballah, M., Shirazi-AdI, A. & Zukor, D. (1995) Biomechanics of the human knee joint

Bendjaballah, M., Shirazi-Adl, A. & Zukor, D. (1997). Finite element analysis of human knee

Berti, G. (Nov 2010). Open Source Software for Medical Simulation, 04/12/2011, Available

joint in varus-valgus, *Clinical Biomechanics* 12 (3), pp. 139-148

from http://www.rheinahrcampus.de/~hartmann/

in comparison reconstruction, mesh generation and finite element analysis. *The* 

acting forces on the articulating surfaces can be studied.

both descriptive and predictive. They comprise

• predictive simulations developed from these collections,

physiology/pathology, by regarding it as an integrated system.

of Tübingen, Germany) for the visualisation of the jumping gymnast.

**5. Conclusion** 

format,

models and

**6. Acknowledgements** 

*Knee* 2 , pp. 69-79

**7. References** 

resource.


**6** 

*Finland* 

**Biomechanics and Modeling** 

Rami K. Korhonen1 and Simo Saarakkala2,3

*1Department of Applied Physics, University of Eastern Finland, Kuopio* 

Articular cartilage is a specialized connective tissue that covers the ends of the bones in the diarthrodial joints. The thickness of human articular cartilage is typically between 1-6 mm. The main functions of articular cartilage are to dissipate and distribute contact stresses during joint loading, and to provide almost frictionless articulation in diarthrodial joints. In order to accomplish these demanding tasks, articular cartilage has unique mechanical properties. The tissue is a biphasic material with an anisotropic and nonlinear mechanical

Articular cartilage is composed of two distinct phases. Fluid phase of the cartilage tissue consists of interstitial water and mobile ions. The water phase constitutes 68-85 % of the cartilage total weight and is an important determinant of the biomechanical properties of the tissue. Solid phase (or solid matrix) of the cartilage tissue consists mainly of collagen fibrils and negatively charged proteoglycans. The cell density is relatively small – in human adult tissue only ~2% of the total cartilage volume is occupied by the chondrocytes. Collagen molecules constitute 60-80% of the cartilage dry weight or approximately 10-20% of the wet weight. The collagen molecules assemble to form small fibrils and larger fibers that vary in organization and dimensions as a function of cartilage depth. The diameter of collagen fibers is approximately 20 nm in the superficial zone and 70-120 nm in the deep zone, and it varies between different collagen types. The collagen fibrils of the cartilage tissue consist mainly of type II collagen, although small amounts of other collagen types can be also found in cartilage, e.g. collagen type VI is common form in the vicinity of cells (pericellular matrix). In addition to the collagen fibrils, proteoglycan macromolecules constitute 20-40% of the cartilage dry weight or approximately 5-10% of the wet weight. The proteoglycan aggrecan is composed of a protein core and numerous glycosaminoglycan (GAG) chains attached to the core. Many aggrecan molecules are further bound to a single hyaluronan chain to form a proteoglycan

The basic structure of the articular cartilage can be divided into four zones based on the arrangement of collagen fibril network (Benninghoff, 1925): 1) *Superficial zone*: here the

**1. Introduction** 

behaviour.

aggregate.

**1.1 Articular cartilage** 

*2Department of Medical Technology, University of Oulu, Oulu* 

*University of Oulu and Oulu University Hospital, Oulu* 

**of Skeletal Soft Tissues** 

*3Department of Diagnostic Radiology* 

