Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 44, 4, pp. 797-818, Warsaw 2006 NUMERICAL COMPLEXITY OF SELECTED BIOMECHANICAL PROBLEMS Marcin Wierszycki Witold Kąkol Tomasz Łodygowski Institute of Structural Engineering, Poznań University of Technology e-mail: Marcin.Wierszycki@put.poznan.pl; Witold.Kakol@put.poznan.pl; Tomasz.Lodygowski@put.poznan.pl The paper describes various aspects of numericalmodeling of biomecha- nical problems by the finite element method. The authors would like to presentwhat theymean by the numerical complexity ofmodeling of bio- mechanical problems. The attention is focused on numerical simulation of dental implants and human lumbar spine motion segment (L4-L5). In both cases, acquisition and creation of geometry, number of DOFs, combining different types of elements, properties of the material, con- tact definitions, loads and boundary conditions are difficult tasks. The acquisition of geometric data of living bodyparts canbe realized only by using noninvasive techniques like NMR or CT. The processing of these data requires specialized software and methods. The methodologies of defining mechanical parameters of human tissues are usually inaccura- te and have to be used in practice on living people very carefully. The constitutive data in literature are usually grossly inconsistent. In nume- rical simulations, custommaterial formulations andmodeling techniques should be used. It is difficult to describe real-world loads and bounda- ry conditions since both are very complex and changing. Load scheme models are global and force values are very difficult to obtain. Boun- dary conditions are necessarily very simplified but they should reflect specific biological behaviors and conditions. Nevertheless, the numerical simulationbymeans of the finite elementmethod canbehelpful anduse- ful during solving biomechanical problems like fatigue analysis of dental implants or estimating the stiffness of human lumbar spine segment. Keywords:biomechanical,finite elementmethod, fatigueanalysis,dental implant, spine segment 798 M. Wierszycki et al. 1. Introduction Nowadays, the finite elementmethod plays the key role in solving engineering problems inmanyfieldsof scienceand industry,andcanbe successfullyapplied also in simulations of biomechanical systems (Będziński, 1997; Geng et al., 2001; Eberlein et al., 2002; Kąkol et al., 2003a; Sakaguchi and Borgersen, 1993). This method has been a well established one, used in biomechanical simulations for over 20 years so far. It allows taking into consideration the key features like material inhomogeneity and anisotropic mechanical properties of tissues as well as a very complicated geometry of human body parts (Eberlein et al., 2002;, Ogurkowska, 1992; Skaggs and Weidenbaum, 1994; Swartz and Wittenberg, 1991). Furthermore, it is proved (Dierich et al., 1992; Eberlein et al., 2002; Hędzelek et al., 2003; Kąkol et al., 2003c) that function, failure, prediction of changes and remodeling of biostructures are related to stress and strain fields in tissues which may be calculated with the help of FEA. Thus, FEA is an efficient tool for testing biomechanical sets, but it is still often very difficult to obtain useful and valuable results for these kinds of problems. Themain reason for this is the complexity of biostructures and the resulting complexity of numerical simulations. These complexities are meant as difficulties of modeling as well as overcoming them by means of specific simulation techniques. These complexities are significant within each step of a numerical analysis. 1.1. Geometry creation The first step of carrying out virtual simulations of biomechanical systems is correct definition ofmodel geometry. Besides topology, additional data such as volume density, can also be interesting. The acquisition of geometric data of living body parts can only be realized by using a noninvasive technique such asNuclearMagnetic Resonance (NMR) orComputerTomography (CT). Other techniques, such as 3D scanningmay also appear to be helpful but only for validation purposes (Rychlik et al., 2004). The preparation of geometric models with the use of NMR or CT consists of two stages. In the first stage, the rows of two dimensional pictures for subsequent slices of the analyzed parts of body are obtained. In the second step, the datamust be processed with a specialized code that converts it into a 3D model in a format readable with the CAD software. In the case of a complex geometry, the above procedure is not a trivial task. In fact, it enables us to prepare quite accurate geometric descriptions of biomechanical models (Ogurkowska et al., 2002a). A more difficult procedure takes place in the transfer from a geometric (CAD) model to a discrete (FEA)model. The geometric data usually cannot Numerical complexity of selected biomechanical problems 799 be used directly in FEA codes. Advanced preprocessors of FEA can import andmesh the surface of a solid geometricmodel prepared in theCADsoftware. In practice however, this procedure is limited only to objects with not really simple but regular shapes. If the geometry is complex, a lot of errors may occur while it is being imported, and the creation of mesh of finite elements may not succeed. The most significant sources of errors are: loss of geometry creation history (parent-child relations between geometry features), number of details (e.g. numberof triangles) andprecision of geometrical descriptions (e.g. order of the polynomial). Because of complexity of the biostructure geometry and errors occurring during the geometry transfer, it is often not possible to automate the process of FEA mesh generation. Finally, the finite element modelmust be created ”by hand” on the basis of information fromNMR, CT and contact scanners (Ogurkowska et al., 2002b). 1.2. Material definition Mechanical properties of tissues are extremely complex. They have usually nonlinear, anisotropic andviscous characteristics (Eberlein et al., 2002; Skaggs and Weidenbaum, 1994). Moreover, these characteristics are always changing as a result of natural processes such as remodeling. Obtaining a reliable and representative value of these characteristics poses another important problem. Experimental methods of obtaining mechanical parameters of biomechanical materials are often not accurate enough and, first of all, cannot be used in research on living people. There is a body of literature on the assessment of physical parameters of tissues, but using such data is often difficult because the data are sometimes grossly inconsistent (Swartz and Wittenberg, 1991). Among others, the reason of this is the heterogeneity of tissues. This heteroge- neity does not cease to becomemore conspicuouswhile testingmaterials taken from various donors (Skaggs andWeidenbaum, 1994). Practical application of FEA will be possible once the actual basic data describing an individual’s traits are provided. The method meeting the requirements above is the one relying on the existing relation between the structure of osseous tissue and its mechanical properties and radiological density (Ogurkowska, 1992). The selection of the right types of material models and possible simpli- fications depends on prepared simulations, e.g. for a dynamic process, the viscoelastic behaviour will play the key role, and for a long term analysis the remodeling phenomenawill bemost significant.Generally, in the case of static or quasi-static stress-strain analyses for ”hard” tissues, like a bone, linear, iso- and orthotropic elastic materials are used. For ”soft” tissues, like ligaments, a hyperelastic material is used (Eberlein et al., 2002). The next very important aspect is choosing and assigning the appropriate type of a finite element. In 800 M. Wierszycki et al. order to describe the behaviour of very complex tissues, a combination of a few different types of elements is often necessary. 1.3. Load and boundary conditions Load schemes and boundary conditions are very complex and difficult to define for biomechanical numerical models. For simulations of biostructures which always constitute a part of a greater entity – a living body – the loads are applied to the parts of body with muscles, ligaments and by contact be- tween them. Furthermore, directions and values of these loads are constantly changing and are very difficult to obtain. Similarly, characteristics of boun- dary conditions are specified by surrounding tissues and organs. The changes in boundary condition characteristics are caused by the phenomena of tissue remodeling. Because of these reasons, it is not possible to describe in full deta- ils the loads and boundary conditions, and in a numerical model they usually must be globalised and simplified. The theoretical models of global load sche- mes for some biomechanical system can be found in literature (Adams et al., 2002; Będziński, 1997). 1.4. Efficiency of simulation The last but not least group of complexities are technical difficulties of nu- merical calculations. For the sake of the simulation complexity, only the size of numericalmodels, convergence of iteration algorithms and time of calculations maybe a problem.Themodels that describe individual features of living body parts or very complex biomechanical structures must contain a very detailed geometric description. These models can have a few hundred thousands ele- ments and degrees of freedom. It is possible to calculate such large problems with the present-day FEA codes bymeans of computers. Still, if the size of a model is near to onemillion DOFs, it may turn out not to be so easy to carry out such an analysis. In the case of large problems, already common work- stations and also low-class servers have significant limitations. The source of these problems can be computer architecture limitations, memory addressing, CPU time needed or size of temporary and result files. If the main goal of numerical calculations of biomechanical structures is their practical application, the efficient methods of simulations in medical treatment and rehabilitation are needed. For such complex structures and nu- mericalmodels, this efficiency requiresmany simplifications and a hierarchical approach (Kąkol et al., 2003c). In the next sections, more detailed aspects of the complexity of numerical models of biomechanical problems will be discussed. Selected problems, such as fatigue analysis of dental implants and evaluating stiffness of the human lumbar spine segment are considered. Numerical complexity of selected biomechanical problems 801 2. Example I – numerical analysis of dental implants Implants are a commonly applied treatment method of dental restorations. Unfortunately, numerous clinical observations point to the occurrence of both early and late complications. In many cases, these problems are caused by mechanical fractures of the implants themselves (Hędzelek et al., 2003; Kąkol et al., 2002). Themost frequent complications are loosening of the connecting screw, fracture and cracking of dental implant parts. While the loosening of the connecting screw causes mostly patient’s discomfort in the implant usa- ge, the cracking leads to much more serious complications and makes further treatment extremely difficult. To understand the reasons of observed mecha- nical complications, it is necessary to know stress and strain fields in implant components aswell as changes in boundary conditions. Cyclic loads anda cha- racter of fracture indicate material fatigue as the basic cause of this fracture (Kocańda, 1985). The authorswould like to confirm this proposition bymeans of numerical simulations of a dental implant system (Kąkol et al., 2003d). The calculationwas carried outwith the commercial codeABAQUS/Standard.For fatigue calculations, the fe-safeWorks was used. In the analysis of a dental implant, the implant structure is not trivial, however it is a purely mechanical problem. The loads and boundary condi- tions have alreadymore complex, biomechanical reasons. Themost significant complexities of finite element stress-strain and fatigue analyses are geome- try andmesh preparations, fatiguematerial characteristics, implant assembly, physiological changes of loads and, finally, definition of boundary conditions describing it as a bone. 2.1. Geometry Numerical models were created on the basis of technical documentation of the commercial implant system OSTEOPLANT. It is a commonly used sys- tem consisting of an implant and abutment with a non rotational hexagonal connection assembled by a screw. For full simulation of the implant structure behaviour, a geometrically complex three dimensionalmodel is necessary.This model, which includes a spiral thread, enables taking into consideration a few important aspects such as full simulation of kinematics of the implant, descrip- tion of the multiaxial state of stress and, the most important, the possibility of simulation of screw loosening. Themost interesting result will be the rela- tionship between the torque, friction coefficient and loosening or fatigue life of the screw under cyclic loads. Unfortunately, a three-dimensional model is very large (ca. 500000 DOF). Due to the fact that most parts of the implant are axisymmetric, connecting the axisymmetric and fully three-dimensional concept of modeling is a good idea of simplification. The ABAQUS/Standard 802 M. Wierszycki et al. offers axisymmetric solid CAXA elements that are capable of modeling nonli- near asymmetric deformation. The CAXA elements are intended for analysis of hollowbodies such as pipes, but theymay also beused tomodel solid bodies with some limitations. These elements use standard isoparametric interpola- tion in the axi-symmetric plane, combinedwith the Fourier interpolation with respect to the angle of revolution. The asymmetric deformation is assumed to be symmetric with respect to the axi-symmetric plane at an angle 0 or π [1, 2]. Fig. 1. Axi-symmetric model of a dental implant with a part of the jaw bone: CB – cancellous bone, CC – cortical bone, R – implant, S – screw, A – abutment An axi-symmetricmodel of an implantwas createdwith the application of CAXA elements (Fig.1). This approach reduces geometry description of the implantmodel from three dimensional to two dimensional. The threads of the screw and the implant body were simplified to axi-symmetric parallel rings. This assumption results in the size of the problem ca. 75000 DOF and the cost of calculation being significantly reduced (Table 1). Table 1.Dimensions of models used in FEA of the implant Axi-symmetric 3-dimensional model model Elements 5891 76345 Nodes 54160 166125 Degrees of freedom 74820 498376 2.2. Material properties Thepart of the jaw is composed of two kinds of bones, the cancellous bone and the cortical one (Fig.1). The problem of describing mechanical constitu- tive law of bones is very complex. Mechanical characteristics and the internal Numerical complexity of selected biomechanical problems 803 microstructure of cortical and cancellous bones are nonhomogeneous, aniso- tropic and variable in time. Changes of bone characteristics are caused by the phenomenon of tissue remodeling. It is very difficult to take these aspects into consideration in implant models. In finite element analysis, many concepts of description ofmechanical properties of bones could be applied – starting from very simple linearly elastic isotropic, going through more complicated, trans- versely isotropic or orthotropic and ending up with very complex, nonlinear anisotropic ones. The assumed material characteristics of the jaw bones are linearly elastic, homogeneous and isotropic (Table 2). This simplification is ju- stifiable due to the role which the bone plays in fatigue analysis of an implant. The most important here is the influence of bone loss around an implant as well as bone flexibility on implant boundary conditions and implant fatigue life (Hędzelek et al., 2003). The implant is made of a medical titanium alloy themechanical properties of which are nonlinear. Their description was based on the certificate of conformity. Table 2. Stiffness of jaw bone Bone Young’s modulus E [MPa] Cortical all schemes 13000 Cancellous 1st scheme 9500 2nd scheme 5500 3rd scheme 1600 4th scheme 690 2.3. Mechanical assembly The implant system is seemingly simple, but in fact it is quite a com- plex mechanical system (Kąkol et al., 2002; Merz et al., 2000; Sakaguchi and Borgersen, 1993). An important aspect of the implant assembly is the mo- deling of tightening. For this purpose, it is necessary to define the contact surfaces between the root, abutment and screw. The friction characteristic on these surfaces is one of the key parameters influencing the preload axial for- ce, reduction of implant components mobility, resistance to screw loosening but also fatigue life of the whole implant. For the friction coefficient, a value ranging between 0.1 (as in a specially finished surface) and 0.5 (as in dry tita- nium to titanium friction)may be found in literature. In the present analysis, three different friction coefficients (0.1, 0.2 and 0.5) were considered. The fric- tion characteristic is one of the key parameters influencing the preload axial 804 M. Wierszycki et al. forces V , reduction of mobility of implant components and screw loosening resistance. The first step of simulation was the tightening of the screw. To simulate this, the middle part of the screw was subjected to temperature loading. The thermal expansion property of the screw material was orthotropic. It was de- fined in such a way that the shrinking occurred only in the direction of the screw axis. The axial force in the tightened screwwas calculated from an em- pirical equation. It depends on the friction coefficient and torque. The effect of this force (ranging from 80 to 850N) was replaced by a temperature field (Kąkol et al., 2003d). In the case of a three-dimensional model, simulation of the tightening can be defined in such a way that it describes a real physical process. If this analysis is performed in an implicit code, a lot of iterations are required, and nevertheless, it is often very difficult to obtain a convergent solution. For simulations of the screw tightening, quasi-static analysis in an explicit code is recommended. In this case, themass scaling is one of themost useful techniques [1, 2]. In the next stage, the simulation was carried out in an implicit code. 2.4. Loads The external loads of the implant model were applied in the second step of simulation. The values and directions of forces were assumed according to a physiologically proven scheme. To estimate the least favorable distribution of stress, only the maximal realistic occlusal forces were taken into account. The loading of the implant was never axial. The vertical component of it is estimated at 600N and the horizontal at 100N. For fatigue calculations, it was necessary to define the character of load changeability in the shape of a load-time curve, the so-called load signal. In the applied low-cycled scheme of 24-hour loads, the average values were 60N (Hędzelek et al., 2003; Kąkol et al., 2003d). 2.5. Boundary conditions In the first stage of implant analysis, all degrees of freedom at the bottom part of the implant bodywerefixed.This assumption seemed tohave its expla- nation in dental practice where nomovements of implants under physiological loads are acceptable. However, the difference between infinitely stiff fixing and even low flexibility is significant, especially in the cyclic loading and fatigue damage context. In thenext stage of implant analysis, theboundaryconditions of implants are modeled as a small part of the jaw bone. The geometry of a small part of the jaw surrounding the implant is very simplified but it enables us to take into consideration the changes in implant fixing conditions. The Numerical complexity of selected biomechanical problems 805 changing flexibility of the bone and bone loss phenomena are also very impor- tant, especially because the bone loss has significant influence on the implant behaviour, stress distribution and, therefore fatigue damage. The degree of encasement and osseointegration of the implantmay not be 100%. It depends on the bone quality, stresses developed during healing and function as well as the location of the implant in the jaw. This percentagemay decrease to as low as 50%. This is caused by phenomena of bone remodeling. In these analyses, three levels of osseointegration were considered. In the case of the first level, the implant body was fully fixed in the jaw bone. In the next two, the degree of the implant body embedding decreased to 75 and 50%, respectively (Fig.2) (Hędzelek et al., 2003). Fig. 2. Levels of dental implant osseointegration in the jaw bone: (a) 100%, (b) 75%, (c) 50% 2.6. Fatigue analysis Fatigue calculations were carried out by the fe-safeprogramwhich uses al- gorithms incorporating themultiaxial plasticity model to estimate the fatigue life. These algorithms are based on the stress results obtained from the finite element analysis (Fig.3), variations in loading, hysteresis loop cycle closure, andmaterial properties.Elastic stresses fromtheFEAmodel are translated in- to elastic-plastic stresses bymeans of biaxial Neuber’s rule and cyclicmaterial properties. In order to estimate the service life of an implant, a designed life is defined. fe-safe calculates the factor (FOS – Factor Of Strength) by which the stresses at each node can be increased or reduced to give the required life. The above is themost interesting and vivid for our case. During a single ana- lysis concerning each node separately, a 6-stress tensor is used to calculate the principal stresses and strains and their orientation. Stress concentration and scale factors are applied at this stage. A rainflow cycle counting algorithm is used to extract fatigue cycles. For biaxial fatigue methods, the critical plane procedure is used to calculate the orientation of the most damaged planes at nodes (3, Bishop and Sheratt, 2000; Draper, 1999). 806 M. Wierszycki et al. Fig. 3. Huber-Mises equivalent stress distribution after bending for three levels of osseointegration: (a) 100%, (b) 75%, (c) 50% The fatigue life of the implant screw was calculated for nine separate ca- ses of loading, three cases of boundary conditions and three cases of the jaw cancellous bone density. For all of these cases, the same cyclic scheme of lo- ading was assumed. A twenty-four-hour changeability schemewas assumed as the signal (Fig.4), while the number of days corresponding to four years was assumed as the number of cycles. Fig. 4. Cyclic scheme of loading – load signal The FOS distribution analysis for particular cases indicates the axial for- ces in the screw and the changes in the scheme of boundary conditions which have the greatest influence on fatigue changes. For different bone densities and at the same time, divergent stiffness of boundary conditions and signifi- cant differences of stress distributions in the screw are noticeable. Yet, it does not lead to serious fatigue changes. For axial forces above 600N, there is a noticeable increase in the areas endangered by fatigue failure. The degree of required stress reduction reaches ca. 30%. In the most unfavorable load case, the maximum axial force results from a high torque and a very small fric- tion coefficient on the screw thread. In two-part implants, this high tightening force is motivated by biological and medical aspects. However, the increase Numerical complexity of selected biomechanical problems 807 in torsion and decrease in friction coefficients entail reduction of the implant components fatigue life (Hędzelek et al., 2003; Kąkol et al., 2003b,d). 3. Example II – numerical analysis of human spine motion segment Themain goal of this studywas to conductmechanical analysis of the human lumbar spine segment which can be helpful in medical treatment, diagnosis and rehabilitation (Glema et al., 2004b;Kąkol et al., 2003a). Special attention was focused on the role of intervertebral disc in the characteristic of motion segment stiffness (Glema et al., 2004a). The function, failure, prediction of pathological changes and the remodeling of the motion spinal segment are all related to stress and strain fields in its tissues andmay be calculated with the FEA. To calculate the stress and strain fields in the case of such a complex biostructure, computationalmodels shouldbase on experimentally determined material properties, realistic geometry, appropriate boundary conditions and load schemes (Ogurkowska et al., 2002a,b). In the case of a practical application of such a simulation of surgery or analyses of spinal equilibrium and stability, simplified models of the lumbar spinalmotion segment should be used.The simplification of the lumbar spinal motion segment is based on a hierarchical approach. The intervertebral disc is replaced with a connector type element. The kinematic and kinetic rela- tionships describing behaviour of the connector elementmay be very complex (Glema et al., 2004a). This behaviour is definedon thebasis of results obtained fromtheFEAsimulationswhichwere carriedoutwith theABAQUS/Standard commercial code. In the next stage, validation of the simplified model of the lumbar spinal motion segment is performed (Glema et al., 2004a). The anatomy of the motion segment is very complex and is of the utmost importance at each stage of numerical simulation of the motion segment. 3.1. Geometry The shape of all parts of the motion segment is shown in Fig.5. The- ir individual anatomical characteristics may be quite significant. In order to construct a finite element model of the spinal motion segment, geometric da- ta of the real object are required. In the case of acquisition of such complex geometric data, special methods are used to scan the 3D geometry of the real object and transform it into a CAD system. Nowadays in medicine, the most popular techniques of scanning are Computer Tomography (CT) andNuclear 808 M. Wierszycki et al. Magnetic Resonance (NMR) (Ogurkowska et al., 2002a). Both of these me- thods have some limitations and their resultsmay not always be used directly in FEA codes. Fig. 5. Lumbar vertebra (L3-L5) – 3D geometrical model In this study, a CADmodel is constructed with the use of the 3D-Doctor software (obtainable from Able Software Corp.) that converts the data from CT to 3D geometric model. Unfortunately, thismodel cannot be used directly for structural analysis. The geometric model is imported into the SolidWorks CAD system. Following a lot of repairing and rebuilding operations, the geo- metric model is imported into the ABAQUS/CAE preprocessor and meshed in the final step. Special techniques are required to perform themeshing step successfully. Partitioning and smoothing of the geometry are indispensable he- re. At the same time, it is important to maintain the most critical geometric parameters. In fact, it is done in hundreds of small steps, and the intermediate results are compared with some existing available data. The final geometry, which is abasis formeshing, has avery complexdescriptionof topology (Ogur- kowska et al., 2002b). 3.2. Material properties Each part of the motion segment needs a different approach with respect to material modeling. The anatomy and function of individual parts of the motion segmentareof great importance in thecase of selection offinite element formulations,modeling techniques and constitutive laws. Themotion segment of a human spine consists of two vertebral bodies, intervertebral disc, posterior elements and spinal ligaments (Fig.6) (Adams et al., 2002; Będziński, 1997). Basically, the vertebra is composed of the cancellous bone and cortical bo- ne. Formodeling of the vertebra, solid and shell elementswere usedwith linear elastic orthotropic material properties. Two kinds of bone tissue were taken Numerical complexity of selected biomechanical problems 809 Fig. 6. FEAmodel of themotion segment: VB – vertebral bodies, ID – intervertebral disc, L – ligaments, PE – posterior elements, ZJ – zygapophysial joints into consideration: cortical and cancellous. For the cortical bone of the verte- bral body, which is a very thin sheet in fact, shell elements were used. For the cancellous part, in turn, solid 3D hexahedral elements were applied (Table 3). The intervertebral disc is composed of three different parts: annulus fibrosus, nucleus pulposus and vertebral endplates. This complex structure enablesmo- vements between the vertebral bodies. These are lateral and sagittal bending, twisting, and small sliding. The annulus fibrosus consists of 10-20 layers of collagen fibers oriented +30◦ and −30◦ with respect to the circumferential axis. The nucleus pulposus is a hydrated gel located at the center of the disc. Table 3.Types of elements andmaterials in FEAmodel Anatomical part of Element Material model the motion segment type number vertebra body shell 6400 elastic, orthotropic solid 13966 elastic, orthotropic intervertebral disc solid 1250 elastic, isotropic surface 1004 elastic or hyperelastic, with rebars isotropic, no compression fluid 564 incompressible fluid vertebral endplates shell 814 elastic, isotropic ligaments membrane 129 elastic or hyperelastic, isotropic, no compression zygapophysial joint shell 56 elastic or hyperelastic, isotropic solid 94 elastic, isotropic Both components maintain the stiffness of the disc against the compres- sion loading and allow for some degree of movement between the vertebral 810 M. Wierszycki et al. bodies. The annulus fibrosus can bemodeled as an anisotropic material or by structural elements which introduce this anisotropy. In this study, the second approach was adopted. The ground matrix of the disc annulus wasmodeled with 3D solid elements. The collagen fibers we- re modeled with rebar type elements embedded into the 3D solid elements (Table 3). The rebar elements are defined as surface layers of uniaxial rein- forcements in the solid elements with uniformly spaced reinforcing bars. Such layers are treated as smeared layers with a constant thickness which is equal to the area of each reinforcing bar divided by the reinforcing bar spacing. This technique allows for easy alterations in thenumber of layers, section properties and its orientation.What ismore, it is independent of the element re-meshing mentioned above. In our study, we calibrated the 3D embedded elements with data published by Skaggs and Weidenbaum (1994). So, the number of rebars and their cross-sectional areas were calculated on the basis of the published data (Skaggs and Weidenbaum, 1994). The nucleus pulposus, assumed here to be incompressible, was modeled as a fluid-filled cavity using hydrostatic fluid elements with initial pressure 2MPa. The hydrostatic fluid elements co- ver boundaries of the nucleus pulposus. They share the nodes at the cavity boundary with the standard elements of the annulus fibrosus (Table 3). The posterior elements control the position of the vertebral bodies such as a pair of stout pillars of a bone called pedicles, two transverse processes, two superior articular processes and two laminae endings with a spinous process. These se- veral processes serve as sites of attachment formuscles that control the lumbar vertebral column. All parts of the spinal segment are connected by ligaments. The most definite ligament is the ligamentum flavum which consists of an elastin that connects the lower end of the internal surface of one lamina to the upper end of the external surface of the lamina below and closes the gap between the consecutive laminae.Thetransverseprocessesare connectedwith thin sheets of collagen fibers.Theopposingedges of spinousprocesses are connectedwith the collagen fibers referred to as the interspinous ligament and the supraspinous ligament. In addition to the ligaments of the posterior elements, the lumbar vertebral column is reinforced by ligaments that connect the vertebral bodies. These are the posterior longitudinal ligament and the anterior longitudinal ligament. The vertebral bodies, the disc and the posterior elements create a complete structure of the lumbar vertebral column which maintains stability andcontrols ofmovements.Thevertebraeare jointedwithzygapophysial joints whichcreate a lockingmechanismbetween them.Theyblockaxial rotationand forward sliding of the lumbar vertebrae. For seven spinal ligaments,membrane elements describing the material devoid of compression were used. For the zygapophysial joints, shell and solid elements were used (Table 3). Numerical complexity of selected biomechanical problems 811 The superior and inferior vertebral endplates are the cartilage plates that cover the superior and inferior aspects of the discs and bind the disc to their respective vertebral bodies. The cartilaginous endplates were modeled with thin shell elements (Table 3).The intervertebral disc and ligaments play an essential role in motion of the spine and in the contribution to the general stiffness of the spinal segment. As it can be seen, it is not sufficient for such a complex structure as the motion segment to use only material characteristics, even if they are very complex. In this case, it is necessary tomix differentmaterial descriptions and modeling techniques.Material complexitiesmust be consideredwhen choosing material constants for the finite elementmodel. There is rich literature on the assessment of physical parameters of bones and soft tissues, but it is difficult to make use of the basic piece of information because the constitutive data of the same human bones are usually grossly inconsistent. The methods of defining biomechanical parameters are inaccurate and cannot be used in prac- tice on living people. The osseous tissue of the spine is very heterogeneous in its repetition structure and its properties change depending on the sample placement in the vertebra. On the basis of the relationship existing between mechanical properties and radiological density of bones, Young’s modulus of the FEAmodel of the vertebral cancellous bone of themotion segment can be determined (Ogurkowska et al., 2002b). The rest of the data was taken from literature (Adams et al., 2002; Eberlein et al., 2002; Skaggs andWeidenbaum, 1994; Swartz andWittenberg, 1991). 3.3. Load and boundary conditions In reality, loads of themotion segment are applied with neighbouring ver- tebrae, muscles, ligament and other connecting tissues. It is very difficult to take into consideration all of them. Thus, it is necessary to use simplifiedmo- dels of load schemes with global equivalent forces andmoments (Adams et al., 2002; Będziński, 1997). For the first stage of the motion segment, analysis in which the separated disc was taken into consideration, representative load schemes with the unit force/moments were applied. The concentrated force and moments were ap- plied to the reference point. In this analysis, the boundary conditions were also global and reduced to one reference point which was fixed. These refe- rence points couple all nodes on the upper and the bottom surfaces of the intervertebral disc. In the second, validated stage, six load cases with loads and moments typical for the L4-L5 motion segment from literature were taken and applied to the whole spinal motion segment model. The load values were assumed as the half of the damage force. These loads were applied to the reference 812 M. Wierszycki et al. point on the upper surface of the vertebral body. For a singlemotion segment, other motion segments, muscles and connecting tissues, etc., constitute the boundary conditions. Of course, it is not possible to take into account such a realistic description of the boundary conditions. In our model, the bottom surface of the vertebral body was constrained against anymovement. 3.4. Analysis of stiffness of the motion segment The whole complex model of the motion segment was validated for four basic load schemes: axial compression, latteral and sagittal bending and tor- sion. On the basis of the obtained results, it was possible to verify the validity and quality of the model definition. These results are the axial displacement and disc bulge for compression, rotation for bending and axial rotation for torsion. They were compared to values taken from literature (Adams et al., 2002). The numerical model of the whole motion segment had 84385 DOF and 24277 elements (Table 3). The complete calculations lasted 22 hours on a standard PCworkstation and needed 34 increments. Fig. 7. FEAmodel of the intervertebral disc: AF – annulus fibrosus, NP – nucleus pulposus, CF – collagen fibers If we take into consideration the time needed for the preparation of such a complex model, this approach is impossible to apply in future commonme- dical examinations. The basic idea of implification of modeling of the motion segment is the replacement of the complex structure of the intervertebral disc (Fig.7) with one connector type element having a complex definition. The elastic behaviour of this type of element can be described as an equivalent stiffness matrix that is in a general case the relation F = KD, where F is the vector (12 components) of generalized forces that act on the segment, D (12 components) is the vector of mutual displacements between bones and K (144 components) is the stiffness operator of the segment. The values of these matrix components will be obtained from FEA simulations and calcu- Numerical complexity of selected biomechanical problems 813 lated with a self made code in the SciLab (open source software similar to MatLab). For the behaviour of the segment in the range of our interest, i.e. small rotations and small displacements, we can use a tangent stiffness matrix be- cause of the linearity of load displacement relations. In order to follow this approach, we examine a separated disc, which is subjected to twelve loading schemes in the first place. By means of application of the unit loads (axial, shear and moments) and by recording displacements at appropriate points, where relative motion can be estimated, we build the compliance matrix and then the stiffness matrix for the model. Four concepts of connector elements were tested.Thefirst threewere based on twonode, twelveDOF (six transver- se and rotationall) elements the behaviour of which was described by means of a 12 by 12 stiffness matrix. Three definitions of this matrix were prepared for full matrix with 144 non-zero components, reduced matrix with 92 signi- ficant non-zero components and symmetrized with 52 non-zero components. The last one was based on a special type of the connector element. It is a two node element, the stiffness of which is defined by a 6×6matrix. This matrix describes relative motion and rotations of these two nodes (Table 4) (Glema et al., 2004a). The last stagewas comparison of the numericalmodels ofmotion segments containing a one element model of the intervertebral disc withmulti-element, complex disc definition. All four concepts of the equivalent element were stu- died.The six load cases discussed abovewere considered (Glema et al., 2004b). The recorded relative displacements and rotations allow for validations of dif- ferent concepts of intervertebral disc simplifications. Selected results are shown in Fig.8. It is seen that the response of the simplified intervertebral disc mo- del using connector elements is the closest to behaviour of the multi-element model. This element is recommended for the hierarchical modeling technique. Note here that the response results (Fig.8) for other elements employed are nearly the same – the curves coincide with each other. 4. Conclusions On the basis of analysis carried out, it was proved that the finite element me- thod enables one to obtain useful and helpful results for such complex struc- tures as dental implants or human lumbar spine motion segments. On each and every step of these analyses, many obstacles need to be overcome. These obstacles as well as the overcoming them in numerical modeling with speci- fic techniques determine the numerical complexity of biomechanical problems of biostructures. If the main goal of numerical calculations of biomechanical 8 1 4 M . W ie r s z y c k i e t a l . Table 4.Comparison of concepts of disc simplifications Complex finite element One element model One element model One element model Connector element model of disc with full matrix with reducedmatrix with symmetrized matrix model of disc 5 types of elements: 1 type of element 1 type of element 1 type of element 1 type of element solid, membrane, shell 12×12 unsymmetric 12×12 unsymmetric 12×12 symmetrized 6×6 symmetrized surface with rebars stiffness matrix stiffness matrix stiffness matrix stiffness matrix fluid 144 non-zero components 92 non-zero components 52 non-zero components 18 components 3747 elements 1 element 1 element 1 element 1 element 8797 dofs 12 dofs 12 dofs 12 dofs 6 dofs 136 s/iteration∗ 1 s/iteration∗ 1 s/iteration∗ 1 s/iteration∗ 1 s/iteration∗ ∗ CPU time on P-IV, 2GHz, 512RAM Numerical complexity of selected biomechanical problems 815 Fig. 8. 8: Validation of different concepts of intervertebral disc simplifications: deformed shape (a) and rotation plot (b) for sagittal bending, deformed shape (b) and rotation plot (c) for axial torsion (1 – complexmodel of disc, 2 – connector element, 3 – all three others) structures is the practical application of them, efficientmethods of simulations in medical treatment and rehabilitation are needed. For such complex struc- tures and numerical models, this efficiency requires many simplifications and hierarchical approaches. The complexity will still be quite significant, so we can suggest tools andmethodologies to make them viable to be carried out. In the case of dental implants, the use of CAXA elements allowed for a limited 3D analysis with the use of two dimensional geometric models. In the case of a full three dimensional model of an implant, the use of explicit codes enabled effectively simulation of the screw tightening. On the basis of the results of fatigue analysis, it can be claimed that the material fatigue is the basic reason of the observed complications. For simulation of the motion segment, each step needs special techniques. The acquisition of geometry and the creation of the finite element model re- quire the use of three different types of software andmany special procedures to process geometric data. The modeling of soft and hard tissues needs diffe- rent types of elements andmaterial models. Preparation of a numericalmodel 816 M. Wierszycki et al. of one motion segment and its analysis is so complex and needs so much ti- me and computer resources that some simplifications are absolutely necessary. The idea of a simplifiedmodel of the intervertebral discwill enablemechanical analysis of much more complex models of human spine which can be helpful for simulation of surgery or analyses of spinal equilibrium and stability. Acknowledgements This research was financially supported the grant No. 3T11F02628 within the budget resources for 2005 and 2006. References 1. 2005,ABAQUS Analysis User’s Manual, ABAQUS, Inc., Pawtucket 2. 2005,ABAQUS Theory Manual, ABAQUS, Inc., Pawtucket 3. 2005, fe-safe Works User’s Manual, Safe Technology Limited, Sheffield 4. AdamsM.A., BogdukN., BurtonK., DolanP., 2002,The Biomechanics of Back Pain, Elsevier Science Limited, Churchill Livingstone 5. BędzińskiR., 1997,Biomechanika inżynierska, OficynaWydawniczaPolitech- nikiWrocławskiej,Wrocław 6. Bishop N.W.M., Sherratt F., 2000,Finite Element Based Fatigue Calcula- tions, NAFEMS, Glasgow 7. Dietrich M., Kędzior K., Wittek A., Zagrajek T., 1992, Non-linear finite elementanalysis of formationand treatmentof intervertebraldischerniae, Proc. Inst. Mech. Eng., 225-231 8. Draper J., 1999,Modern Metal Fatigue Analysis, HKS, Inc. Pawtucket 9. Eberlein R., HolzapfelG.A., Schulze-BauerC.A.J., 2002, Assessment of a spinal implant by means of advanced FE modeling of intact human inte- rvertebral discs,Fifth World Congress on Computational Mechanics 10. Geng J.P., TanK.B., LiuG.R., 2001,Application of finite element analysis in implant dentistry. A review of literature, J. Prosthet. Dent., 85, 585-598 11. Glema A., Łodygowski A., Kąkol W., Wierszycki M., Ogurkowska M.B., 2004a, A simplified vertebrae disk model for computations of human lumbar spine segment, Sixth World Congress on Computational Mechanics 12. Glema A., Łodygowski A., Kąkol W., Wierszycki M., Ogurkowska M.B., 2004b,Modeling of intervertebraldiscs in thenumerical analysis of spinal segment,ECCOMAS 2004 13. Hędzelek W., Zagalak R., Łodygowski T., Wierszycki M., 2003, The effect of marginal bone loss and bone density on the risk of late implant com- ponents failures, 27th Annual Conf. of the European Prosthodontic Association Numerical complexity of selected biomechanical problems 817 14. KąkolW., Łodygowski T., OgurkowskaM.B.,WierszyckiM., 2003a, Are we able to support medical diagnosis or rehabilitation of human vertebra by numerical simulation, 15th Int. Conf. on Computer Methods in Mechanics 15. KąkolW., Łodygowski T., OgurkowskaM.B.,WierszyckiM., 2003b, Numerical modeling of spinal motion segment,ABAQUS User’s Conf. 16. Kąkol W., Łodygowski T., Ogurkowska M., Wierszycki M., 2003c, The support of medical diagnosis or rehabilitation of the human vertebra by numerical simulation,Acta of Bioengineering and Biomechanics, 5, 1, 235-241 17. Kąkol W., Łodygowski T., Wierszycki M., 2003d, Estimate of tooth implant fatigue under cycling loading, 15th Int. Conf. on Computer Methods in Mechanics 18. KąkolW.,ŁodygowskiT.,WierszyckiM.,HędzelekW.,ZagalakR., 2002, Numerical analysis of the behavior of dental implant, ABAQUS Users’ Conference 2002 19. Kocańda S., 1985, Zmęczeniowe pękanie metali, Wydawnictwo Naukowo- Techniczne,Warszawa 20. Merz B.R., Hunenbart S., Belser U.C., 2000,Mechanics of the implant- abutment connection: an 8-degree taper compared to a butt joint connection, The International Journal of Oral and Maxillofacial Implants, 15, 4, 519-526 21. Ogurkowska M.B., 1992, Application of radiology and rheology method for mechanical testing of the vertebral column, Doctoral Thesis. University School of Physical Education, Poznań 22. Ogurkowska M.B., Rychlik M., Stankiewicz W., Nowak M., Roszak R., Glema A., Wierszycki M., Morzyński M., Łodygowski T., 2002a, The interaction of the L4-L5 spinal segments by FEA analysis. Part 1. Me- thods of geometrical data acquisition and validation, Acta of Bioengineering and Biomechanics, 13th Conference of the European Society of Biomechanics, 4, 1, 98 23. Ogurkowska M.B., Rychlik M., Stankiewicz W., Nowak M., Roszak R., Glema A., Wierszycki M., Morzyński M., Łodygowski T., 2002b, The interaction of the L4-L5 spinal segments by FEM analysis. Part 2. Virtu- al modeling of the structure, Acta of Bioengineering and Biomechanics, 13th Conference of the European Society of Biomechanics, 4, 1, 99 24. Sakaguchi R.L., Borgersen S.E., 1993, Nonlinear finite element contact analysis of dental implant components,The International Journal of Oral and Maxillofacial Implants, 7, 1, 655-661 25. SkaggsD.L.,WeidenbaumM., 1994,Regional variation in tensile properties and biomechanics composition of the human lumbar annulus fibrous, Spine, 9, 120-134 26. Swartz D.E., Wittenberg R.H., 1991, Physical andmechanical properties of calf lumbosacral trabecular bone, J. Biomechanics, 24, 11, 1059-1068 818 M. Wierszycki et al. 27. Rychlik M., MorzyńskiM., Nowak M., Stankiewicz W., Łodygowski T., Ogurkowska M.B., 2004, Acquisition and transformation of biomedical objects to CAD systems, Strojnicky Casopis, 3, 4, 121-135 Numeryczna złożoność wybranych problemów biomechaniki Streszczenie Praca opisuje wybrane aspekty modelowania numerycznego z wykorzystaniem MetodyElementówSkończonych (MES) zagadnień biomechaniki. Autorzy starają się podkreślić złożoność takiej analizy MES. Uwagę skupiono na dwóch przykładach: analizie wszczepu stomatologicznego oraz ruchomości segmentu kręgosłupa ludzkiego (L4-L5). W obu analizowanych i modelowanych przypadkach odtworzenie geometrii, liczba typów elementów i stopni swobodymodeli, przyjęciemodeli i właściwości kon- stytutywnych materiałów biologicznych, definicje problemów kontaktowych czy wa- runków początkowo-brzegowych stanowią o poziomie trudności podejmowanych za- dań. Pokonaniewymienionych trudności, a co za tym idzie zbudowaniemodeli nume- rycznych odtwarzających zachowanie się rzeczywistych elementów biomechanicznych jest przedmiotem podejmowanej w pracy dyskusji. Pomimo tej złożoności zadań, za- stosowanie MES do oceny zachowania się implantów bądź segmentu ruchowego pro- wadzi do wystarczająco zgodnej z eksperymentem oceny trwałości implantów lub sztywności segmentów kręgosłupa. Manuscript received February 2, 2006; accepted for print May 5, 2006