Advanced finite element analysis of L4–L5 implanted spine segment

In the paper finite element (FE) analysis of implanted lumbar spine segment is presented. The segment model consists of two lumbar vertebrae L4 and L5 and the prosthesis. The model of the intervertebral disc prosthesis consists of two metallic plates and a polyurethane core. Bone tissue is modelled as a linear viscoelastic material. The prosthesis core is made of a polyurethane nanocomposite. It is modelled as a non-linear viscoelastic material. The constitutive law of the core, derived in one of the previous papers, is implemented into the FE software Abaqus®. It was done by means of the User-supplied procedure UMAT. The metallic plates are elastic. The most important parts of the paper include: description of the prosthesis geometrical and numerical modelling, mathematical derivation of stiffness tensor and Kirchhoff stress and implementation of the constitutive model of the polyurethane core into Abaqus® software. Two load cases were considered, i.e. compression and stress relaxation under constant displacement. The goal of the paper is to numerically validate the constitutive law, which was previously formulated, and to perform advanced FE analyses of the implanted L4–L5 spine segment in which non-standard constitutive law for one of the model materials, i.e. the prosthesis core, is implemented.

ankylosis after an average follow-up of 17 years. The other studies [9] indicated that only in 50 % of the cases, where two adjacent prostheses were applied, a normal kinematics of the operated spine region was achieved after the average follow-up 12.4 years. Since the age of the patients that undergo total disc replacement is usually between 30 and 50 years a prosthesis used during the surgery should last for at least 40 years. The authors of the research related to total disc replacement [10] stated that it is doubtful that the lifetime of the disc prosthesis designs now available is that long. Therefore, studies on long-term complications after the prostheses implantations and on the constructions of those prostheses are highly recommended.
Clinical result of disc prosthesis application may be enhanced by implementation of a new material, especially with regard to the prosthesis core. Recently, more and more often the intervertebral disc prosthesis core has been made of polymers, in particular polyurethanes, see e.g. [11,12]. Polyurethanes have been used for many years for production of implant components [13,14]. They are characterised by high-degree biocompatibility, high resistance to degradation in human body and high wet angle, which is very important to avoid cell growth in the implantation region.
Introduction of a new material in the prosthesis requires many mechanical and biological tests. Computer simulation seems to be a very helpful and useful tool in the area of a new material study in the context of its application in orthopaedics. Usually such simulations are performed by means of finite element (FE) analyses. However, in order to correctly simulate the material behaviour one has to provide an adequate mathematical model determining the relation between stress and strain, i.e. a constitutive law.
The main objective of the paper is to validate the previously formulated non-linear viscoelastic constitutive law for the polyurethane nanocomposite by making FE analyses of the L4-L5 implanted segment. The intervertebral disc prosthesis consists of two metal endplates and the core made of the polyurethane nanocomposite material. The implementation of the constitutive law for the core into the FE commercial software Abaqus is a novelty aspect of the presented research.
In the present paper we are not considering the possible reconstruction of parts of the spine by means of the addition of some bioresorbable materials. The biomechanical performances of the whole system implanted prosthesis-growing bone may be greatly improved by the application of such a technique, even if serious theoretical and experimental investigations are still needed [15][16][17][18]. In addition to this, another problem related to an implant application, i.e. prosthesis micromotions, is not taken into account. This is an important problem, therefore, the present investigations may need to be completed by paralleling those presented in [19].

Derivation of fourth order stiffness tensor to be applied in Abaqus
The constitutive law, which was implemented into Abaqus , was formulated for the polyurethane nanocomposite using the strain energy function approach. Studies conducted in this regard showed that the Mooney-Rivlin hyperelastic function and the Ogden function (N = 4) were the most adequate ones for modelling pure elastic response of the material. In addition, strain history was taken into consideration by means of hereditary integral. The constitutive model was formulated on the basis of experimental tests. It is capable to simulate stress relaxation and hysteresis loop. The detailed mathematical derivation of the constitutive model for the polyurethane material was presented in [20].
The fact that the constitutive model is to be implemented in FE analyses requires that the potential function Ψ e be split into volumetric part Θ and isochoric part Ψ : The isochoric-volumetric decoupling of the energy function is frequently applied in the context of its application in displacement-based finite element analyses. An advantage of the decoupling is that the isochoric and the volumetric material behaviour can be treated as completely independent. This permits their decoupled treatment in the FE simulations and in derivation of stiffness tensor to be implemented in FE analyses. The general constitutive equation for the non-linear visco-hyperelastic material presented in [20] is of the form: where S-second Piola-Kirchhoff stress tensor, g i -viscoelastic constants, g ∞ = 1− n i=1 g i , τ i -relaxation times.
The elastic part of the stress S e in (5) can be calculated: where the following identities were utilised: In Eqs. (3) and (4) l is a fourth-order identity tensor, the sign ":" denotes double contraction of two tensors, "⊗" represents the dyadic product of two tensors. Tensor C is calculated using equation: where C is the right Cauchy deformation tensor. The formula for the material stiffness fourth-order tensor C is as follows: After implementation of the chain rule one obtains the following form of the decoupled stiffness tensor C: In the contemporary finite element softwares, including Abaqus , it is required that the Zaremba-Jaumann objective stress rate, ∇ T, of Kirchhoff stress T be introduced: where: W and D are the symmetric and anti-symmetric parts of the spatial velocity gradient L, respectively; C T ZJ is the tangent modulus tensor for the Zaremba-Jaumann rate of the Kirchhoff stress. The components of the tensor C T ZJ can be calculated by means of Eq. (10): where δ i j , δ jl are the Kronecker deltas. The tensor C Tc i jkl is associated with the convected rate of the Kirchhoff stress and is obtained by the use of the following transformation rule: The components of C T ZJ i jkl were introduced in Abaqus software by means of UMAT user-subroutine code. As the material in interest is modelled as non-linearly viscoelastic, those components were multiplied by the viscous constituent g v : where t is a time step of numerical calculations. This results from the theory described in [20]. Another interesting example of the stiffness tensor derivation for a two-dimensional case is presented in [22]. The authors showed two-dimensional constitutive equations for a plate starting from the equations of the nonlinear elastic body and describing the small deformations superposed on the finite deformation.
In addition, Cauchy stress was calculated using the formula: To make the set of needed equations complete an expression for Kirchhoff stress has to be introduced: In Eq. (13) the second Piola-Kirchhoff stress tensor S (total stress) was calculated by means of relation (2) and the hyperelastic part of the stress. S e by means of Eq. (3). It has to emphasised, that the hereditary integral in (2) was determined using the theory presented in [23]. The number n in (2) was determined using data from stress relaxation test and denotes the number of Maxwell elements in the viscoelastic model of the material. The values of the hyperelastic constants c 10 , c 01 and viscoelastic constants g i , τ i (i = 1, . . . , n = 7) are listed in Table 1.

FE model of the lumbar spine segment
The geometrical models of the L4 and L5 vertebrae were created on the basis of computer tomography (CT) projections of the segment. The CT data, which was acquired in DICOM format, was first properly processed, i.e. bone tissue thresholding and segmentation were performed, and then transferred into ProENGINEER CAD system. The solid models of the vertebrae were created, as well as the model of the intervertebral disc prosthesis. The geometrical model of the whole implanted lumbar segment was transferred into Abaqus , in which the numerical model was created. The segment model consisted of approx. 300 000 tetrahedral finite elements. The model is presented in Fig. 1, which also shows symbolically the boundary and load conditions. All the degrees of freedom were fixed in the nodes of the lower surface of the L5 vertebral body. The load was subjected onto the upper surface of the L4 vertebral body. Between the metallic endplates and the vertebra bodies contact without the possibility of sliding was defined. On the other hand, the defined interaction between the polyurethane core and the endplates allowed for free adjustment of the core position according to the applied force.
The model of the prosthesis, which is based on the other artificial intervertebral disc described in [24], is shown in Fig. 2. Both upper and lower endplates have flanges that are to prevent the core from slipping from between the plates. The spherical surfaces facilitate positioning of the core with respect to the metallic endplates. The core is modelled as a non-linear viscoelastic material by means of the constitutive law implemented into Abaqus using UMAT subroutine. The endplates are modelled as an elastic material. Their Young modulus is equal E p = 200 000 MPa and Poisson's ratio ν p = 0.3.  Bone tissue is modelled as a linear viscoelastic material taking advantage of the Abaqus material library. The mechanical properties change with time is defined by means of the Prony series [25]: where: K (t) , G(t) are bulk modulus and shear modulus of bone, respectively; G 0 -instantaneous shear modulus, K 0 -instantaneous bulk modulus. The values of viscoelastic constants g i and relaxation times τ i were determined on the basis of relaxation test. The exact methodology can be found in [20]. The values of those quantities (for bone i = 1, . . . , n = 8) are listed in Table 2. In the FE analyses the young modulus of bone is equal E = 12 000 MPa, and Poisson's ratio ν = 0.2 [26]. The instantaneous shear and bulk moduli were calculated using formulae: In the first analysis the segment was subjected monotonically within 30 s to the force 3,000 N and then unloaded at the same rate to 0 N. The maximal force corresponded to the load acting on the disc during sitting position of a man weighing 70 kg and holding a weight 200 N in his hands [27].
In the second analysis stress relaxation test was simulated. The segment was loaded kinematically to the same value of corresponding force within 1 s. The exerted displacement was then held unchanged for approx. 15 min.

Results
The volume-changing part was of the following form: where 1/D represents bulk modulus of the infinitesimal theory, J is the determinant of the deformation gradient tensor F. In numerical simulations of incompressible materials the constant D is assumed to be very small making, thus, the bulk modulus tends to infinity. In the presented FE simulation D was 1 · 10 −10 . In paper [21] the authors discuss various forms of volumetric potential functions. In the presented research the simplest possible form of the function Θ was assumed. The volume-preserving part Ψ was assumed to be of the Mooney-Rivlin form, i.e.: where I 1 and I 2 are the first and second invariants of the distortional part of the right Cauchy tensor C, respectively, c 10 , c 01 are hyperelastic constants. For the assumed form of volumetric potential function (19) and isochoric strain energy function (20) the decoupled stiffness tensor C can be derived using Eq. (8): where: In Eqs. (21)-(24) trC 2 denotes the trace of C 2 , I-second-order identity tensor, the components of I C −1 are: Using Eq. (3) the elastic part of the second Piola-Kirchhoff stress can be calculated: Kirchhoff stress tensor is then derived using Eqs. (13) and (14): Now implementing Eqs. (10) and (11) the components of the fourth-order tensor C Tc i jkl can be calculated: where: In Eqs. (27) The components C Tc i jkl , multiplied by g v (see Eq. (12)), were implemented into the analyses by means of UMAT subroutine. Those components define the Jacobian matrix of the constitutive model DDSDDE in the code. The matrix defines the change in the ith stress component at the end of the time increment caused by an infinitesimal perturbation of the jth component of the strain increment array [28]. Before the constitutive model for the polyurethane material was used in the FE analyses of the implanted lumbar segment, correctness of the model implementation by means of UMAT subroutine was verified. This was done in a simple analysis of a cube 10 × 10 × 10 cm (Fig. 3a) which was compressed by 10 % within 60 s and decompressed to strain 0 % at the same rate. Also a relaxation test was simulated. The cube was compressed to strain 10 % within 1 s and the change in stress was calculated for 120 s. The results of the model numerical validation in Abaqus are shown in Fig. 3b, c. In those figures solid lines represent stress (the component along the load direction) calculated in Abaqus and circles represent stress in the same direction calculated by means of Eq. (2).
In Fig. 4 nominal strain and viscoelastic stress distributions for the first case of the segment load are shown. Figure 5 presents distribution of Huber-von Mises stress in the whole segment. Those distributions were plotted at the moment when the highest value of force was acting on the lumbar segment. The stress-strain characteristic of the polyurethane material determined in the FE analysis is shown in Fig. 6. In Figs. 3b and 6 horizontal axes represent the stretch ratio in the direction of the load while vertical axes-viscoelastic stress (2) in the same direction.
The best way to present the results of the second case of the segment load, i.e. the stress relaxation test, is to show the curve representing change of the vertical component of stress (2) with time (Fig. 7). This curve, as well as the curve in Fig. 6, was drawn on the basis of stress components calculated in a chosen finite element. The element was selected from the most loaded region of the polyurethane material.

Discussion
In the paper the numerical application of the constitutive law is presented. The constitutive model was formulated in [20] for a polyurethane nanocomposite, which is to be applied as the core of intervertebral disc prosthesis. The material was treated as non-linearly viscoelastic.
In the considered prosthesis structure the layer of the polyurethane inlay is rather thin relatively with the whole thickness of the prosthesis. In connection with this it is interesting to derive some 2D finite elements or use some 2D formulations as in [29,30].
In order to implement the constitutive model into FE software Abaqus the strain energy function was split in isochoric and volumetric parts. The expression for the second Piola-Kirchhoff stress components was derived and the stiffness fourth-order tensor was determined. Push-forward operation of those tensors was then performed and their components were introduced in the UMAT subroutine.
The polyurethane material applied as one of the intervertebral disc prosthesis components was manufactured at Warsaw University of Technology, Warsaw, Poland. It was, thus, a new material of unknown properties. Its possible clinical application has to be preceded by computer simulations by means of finite element method [31,32]. In order to simulate the mechanical behaviour of the new material a constitutive law was formulated for it. The experimental studies performed on the polyurethane samples and described in [20] showed that the material indicated rheological effects, i.e. stress relaxation, hysteresis and rate-dependence of stress-strain characteristics. Those effects could not be simulated by standard materials available in the material library of any FE software including Abaqus . In addition, the effects could not be simulated by means of a linear viscoelastic constitutive model. A new non-linear constitutive law had to be formulated which would be capable to model the rheological properties of the material. The model was based on stress relaxation and monotonic compression tests. The number of the relaxation times τ i was determined using stress relaxation data. The hyperelastic and viscoelastic constants were calibrated on the basis of monotonic compression tests results. As the compression tests were conducted at three strain rates and the calibration was performed for the strain rates simultaneously the constitutive model was rate dependent. During the compression tests the polyurethane samples were loaded to a certain strain level and unloaded. The created constitutive model is, then, capable to capture the hysteresis phenomenon.
The novelty aspect of the paper is a FE application of the new constitutive law for a new polyurethane material in FE analyses. To complete the task the fourth-order stiffness tensor was derived and implemented into FE analyses in Abaqus by means of UMAT subroutine. Definition of the new stiffness tensor allows for simulation the rheological effects that the real polyurethane material indicates, i.e. stress relaxation and hysteresis.
The numerical validation of the non-linear viscoelastic constitutive law was completed by means of simple FE analyses. The results of those analyses proved that the implementation of the derived fourth-order stiffness tensor by UMAT subroutine was correct. One can see a perfect match of stress calculated in Abaqus with that calculated with Eq. (2). The positive numerical validation of model (2) allowed us to utilise it in FE analyses of the implanted L4-L5 spine segment. The bone tissue of the vertebrae was simulated as a linear viscoelastic material. The mechanical properties were determined on the basis of the experimental tests similar to those described in details in [20] that are in good agreement with the results presented in [26].
The strain and stress distributions presented in Fig. 4 indicate that the spherical region of the polyurethane core, which was in the direct contact with the metal plates, was heavily compressed. The core adjusted its position according to the direction and value of the load acting on the upper vertebra. Stress concentrations in bone tissue are placed in the vicinity of the spikes of the prosthesis metal plates and in the vertebral processes. Figures 6 and 7 show the typical response of a rheological material to loading-unloading simulated test and relaxation test.
In the studies presented in [33] the authors made FE analyses of an implanted L4-L5 segment. The prosthesis core was made of a polyurethane material. The mechanical behaviour of the core was simulated by means of a non-linear hyperelastic Gent model. The authors performed the analyses in order to evaluate the characteristics of force-displacement and moment-angle for the new design of the prosthesis with a separate inlay under various loads. They concluded that the presented prosthesis requires some modifications in the shape of its elements. Their studies, however, proved that prostheses with polyurethane cores are able to mimic the behaviour of the human intervertebral disc.
The aim of the study presented in [34] was to determine the range of motion and facet joint forces in L4-L5 segment. The authors made FE analyses of the Charité prosthesis and two types of Slide-Disc prosthesis, i.e. unconstrained design and constrained one. The cores of those prostheses were made of the polyethylene (UHMWP). In spite of that fact it is a viscoelastic solid [35] it was modelled as an elastic material, which is a kind of a simplification. The authors indicated that facet forces were strongly dependent on the implant design. However, they finally concluded that the forces may be more dependent on the individual spine geometry rather than a specific implant design.
Also in other FE analyses [36,37] the cores of artificial replacements of the intervertebral disc were modelled as elastic or hyperelastic materials. The current studies reveal the possibility to model polymer elements of the disc prostheses by means of non-linear viscoelastic constitutive law. The methodology is very useful especially in the case of the application of a new polyurethane material that exhibits such mechanical behaviour. The results of such advanced FE analyses are, thus, more reliable because the material of interest is not modelled by a standard constitutive law but by a law that captures rheological properties the material shows. The analyses presented in the paper can be utilised in prosthesis shape design and/or material optimisation. Also they may be further extended for damage and poromechanical effects [38,39]. The limitation of the constitutive model presented in the paper is that it was formulated on the basis of relaxation tests and monotonic compression tests. Inclusion of, for instance, shear tests would make the model more universal. However, in that case the material constants calibration would be more difficult or even impossible within a reasonable error.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.