Biomechanical Computer Models

1. Introduction

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:

How can we optimise movements (e.g. in sports or rehabilitation) in order to minimise the loads on the joints? How can we better understand mechanisms of injury and thus improve prevention? These questions are related to the human body.

How do we have to design the equipment to optimally suit the patient’s or the athlete’s requirements in terms of mechanical behaviour? These questions are related to the medical or technical equipment. This optimisation process has three aspects:

• 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).

• 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: Multi body systems (MBS) have been set up yielding important results in prevention and sports (Gruber et al.,1998; Blajer & Czaplicki, 2001; Wojtyra, 2003); 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 (Wangerin et al., 2007).

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 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.

$$F=ma$$E1

$$M=J*a$$E2

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

Figure 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.

Figure 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).

Figure 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 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:

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

2. Defining the physical behaviour

3. Setup of the stiffness matrix

4. 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

$$KU=F$$E3

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

$$KU={\omega }^{2}MU$$E4

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.

$$M\underset{}{\overset{\mathrm{..}}{{\displaystyle U}}}(t)+M\underset{}{\overset{.}{{\displaystyle U}}}(t)+KU(t)=0$$E5

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).

Figure 4.

From medical scan data to finite elements.

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.

Figure 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.

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 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.

Figure 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).

4.2.3. Segmentation process

The segmentation of the images is done to define anatomical structures in the three-dimensional 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 semi-automatic 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 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.

Figure 7.

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

Figure 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.

Figure 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 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.

Figure 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 reconstruction of the human knee joint is shown in Figure 10.

4.2.5. Implementation of the model

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 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/m³. 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 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.

Figure 11.

Calculated contact forces of the simulation of a knee bending.

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

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

• predictive simulations developed from these collections,

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

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

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.

Acknowledgments

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.