Finite element implementation of a generalized Fung-elastic constitutive model for planar soft tissues

Numerical simulations of the anisotropic mechanical properties of soft tissues and tissue-derived biomaterials using accurate constitutive models remain an important and challenging research area in biomechanics. While most constitutive modeling efforts have focused on the characterization of experimental data, only limited studies are available on the feasibility of utilizing those models in complex computational applications. An example is the widely utilized exponential constitutive model proposed by Fung. Although present in the biomechanics literature for several decades, implementation of this model into finite element (FE) simulations has been limited. A major reason for limited numerical implementations are problems associated with inherent numerical instability and convergence. To address this issue, we developed and applied two restrictions for a generalized Fung-elastic constitutive model necessary to achieve numerical stability. These are (1) convexity of the strain energy function, and (2) the condition number of material stiffness matrix set lower than a prescribed value. These constraints were implemented in the nonlinear regression used for constitutive model parameter estimation to the experimental biaxial mechanical data. We then implemented the generalized Fung-elastic model into a commercial FE code (ABAQUS, Pawtucket, RI, USA). Single element and multi-element planar biaxial test simulations were conducted to verify the accuracy and robustness of the implementation. Results indicated that numerical convergence and accurate FE implementation were consistently obtained. The present study thus presents an integrated framework for accurate and robust implementation of pseudo-elastic constitutive models for planar soft tissues. Moreover, since our approach is formulated within a general FE code, it can be straightforwardly adopted across multiple software platforms.     

Effect of glenoid prosthesis design on glenoid bone remodeling: Adaptive finite element based simulation

Glenoid prosthesis loosening is the most common cause for revision total shoulder arthroplasty. Improved glenoid prosthesis design requires looking beyond initial post-implantation static stress analyses. Adaptive bone remodeling simulations based on Wolff’s law are needed for predicting long-term glenoid prosthesis results. This study demonstrates the capability of predicting glenoid bone remodeling produced by changing prosthesis design features. Twelve glenoid prostheses were designed to fit each of six donor human glenoids, using combinations of three peg types and four backing-peg material combinations (polyethylene and or metal). The twelve FE prosthesis models were individually combined, simulating surgical implantation, with the glenoid models. Remodeling simulations, using a validated adaptive bone remodeling simulation, commenced with homogeneous glenoid bone density. To produce bone remodeling, center, posterior-offset, and anterior-offset physiologic loads were consecutively applied to the bone–prosthesis FE models for 300 iterations. Upon completion, region-specific mean predicted bone apparent densities were compared between bone–prosthesis and intact glenoid FE models. Metal fixations significantly increased proximal-center bone density. Polyethylene fixations resulted in bone density approximately equal to intact. Two angled polyethylene peg designs with longer-anterior and shorter-posterior pegs, reflecting natural glenoid shape, best maintained mid and distal glenoid bone density. While these initial results were not validated, they demonstrate the capability and potential of adaptive glenoid bone remodeling simulation. We expect FE glenoid bone remodeling simulations to become powerful and robust tools in the design and evaluation of glenoid prostheses.     

A patient-specific finite element methodology to predict damage accumulation in vertebral bodies under axial compression, sagittal flexion and combined loads

Due to the inherent limitations of DXA, assessment of the biomechanical properties of vertebral bodies relies increasingly on CT-based finite element (FE) models, but these often use simplistic material behaviour and/or single loading cases. In this study, we applied a novel constitutive law for bone elasticity, plasticity and damage to FE models created from coarsened pQCT images of human vertebrae, and compared vertebral stiffness, strength and damage accumulation for axial compression, anterior flexion and a combination of these two cases. FE axial stiffness and strength correlated with experiments and were linearly related to flexion properties. In all loading modes, damage localised preferentially in the trabecular compartment. Damage for the combined loading was higher than cumulated damage produced by individual compression and flexion. In conclusion, this FE method predicts stiffness and strength of vertebral bodies from CT images with clinical resolution and provides insight into damage accumulation in various loading modes.     

Adaptive glenoid bone remodeling simulation

Glenoid prosthesis loosening is the most common cause for revision total shoulder arthroplasty. Stress-induced bone remodeling may compromise long-term prosthesis fixation and significantly contribute to loosening. Realistic, robust analysis of bone-prosthesis constructs need to look beyond initial post-implantation mechanics provided by static finite element (FE) simulation. Adaptive bone remodeling simulations based on Wolff's law are needed for evaluating long-term glenoid prostheses fixation. The purpose of this study was to take a first step towards this goal and create and validate two-dimensional FE simulations, using the intact glenoid, for computing subject-specific adaptive glenoid remodeling. Two-dimensional glenoid FE models were created from scapulae computed tomography images. Two distinct processes, “element” and “node” simulations, used the forward-Euler method to compute bone remodeling. Initial bone density was homogeneous. Center and offset load combinations were iteratively applied. To validate the simulations we performed location-specific statistical comparisons between predicted and actual bone density, load combinations, and “element” and “node” processes. Visually and quantitatively “element” simulations produced better results (p>0.22), and correlation coefficients ranged 0.51–0.69 (p<0.001). Having met this initial work's goals, we expect subject-specific FE glenoid bone remodeling simulations together with static FE stress analyses to be effective tools for designing and evaluating glenoid prostheses.     

Concept and development of an orthotropic FE model of the proximal femur

PURPOSE: In contrast to many isotropic finite-element (FE) models of the femur in literature, it was the object of our study to develop an orthotropic FE "model femur" to realistically simulate three-dimensional bone remodelling. METHODS: The three-dimensional geometry of the proximal femur was reconstructed by CT scans of a pair of cadaveric femurs at equal distances of 2mm. These three-dimensional CT models were implemented into an FE simulation tool. Well-known "density-determined" bony material properties (Young's modulus; Poisson's ratio; ultimate strength in pressure, tension and torsion; shear modulus) were assigned to each FE of the same "CT-density-characterized" volumetric group.In order to fix the principal directions of stiffness in FE areas with the same "density characterization", the cadaveric femurs were cut in 2mm slices in frontal (left femur) and sagittal plane (right femur). Each femoral slice was scanned into a computer-based image processing system. On these images, the principal directions of stiffness of cancellous and cortical bone were determined manually using the orientation of the trabecular structures and the Haversian system. Finally, these geometric data were matched with the "CT-density characterized" three-dimensional femur model. In addition, the time and density-dependent adaptive behaviour of bone remodelling was taken into account by implementation of Carter's criterion. RESULTS: In the constructed "model femur", each FE is characterized by the principal directions of the stiffness and the "CT-density-determined" material properties of cortical and cancellous bone. Thus, on the basis of anatomic data a three-dimensional FE simulation reference model of the proximal femur was realized considering orthotropic conditions of bone behaviour. CONCLUSIONS: With the orthotropic "model femur", the fundamental basis has been formed to realize realistic simulations of the dynamical processes of bone remodelling under different loading conditions or operative procedures (osteotomies, total hip replacements, etc).     

Computational efficiency of numerical approximations of tangent moduli for finite element implementation of a fiber-reinforced hyperelastic material model

In this study, we evaluated computational efficiency of finite element (FE) simulations when a numerical approximation method was used to obtain the tangent moduli. A fiber-reinforced hyperelastic material model for nearly incompressible soft tissues was implemented for 3D solid elements using both the approximation method and the closed-form analytical method, and validated by comparing the components of the tangent modulus tensor (also referred to as the material Jacobian) between the two methods. The computational efficiency of the approximation method was evaluated with different perturbation parameters and approximation schemes, and quantified by the number of iteration steps and CPU time required to complete these simulations. From the simulation results, it can be seen that the overall accuracy of the approximation method is improved by adopting the central difference approximation scheme compared to the forward Euler approximation scheme. For small-scale simulations with about 10,000 DOFs, the approximation schemes could reduce the CPU time substantially compared to the closed-form solution, due to the fact that fewer calculation steps are needed at each integration point. However, for a large-scale simulation with about 300,000 DOFs, the advantages of the approximation schemes diminish because the factorization of the stiffness matrix will dominate the solution time. Overall, as it is material model independent, the approximation method simplifies the FE implementation of a complex constitutive model with comparable accuracy and computational efficiency to the closed-form solution, which makes it attractive in FE simulations with complex material models.     

A constitutive model for the time-dependent,nonlinear stress response of fibrin networks

Blood clot formation is important to prevent blood loss in case of a vascular injury but disastrous when it occludes the vessel. As the mechanical properties of the clot are reported to be related to many diseases, it is important to have a good understanding of their characteristics. In this study, a constitutive model is presented that describes the nonlinear viscoelastic properties of the fibrin network, the main structural component of blood clots. The model is developed using results of experiments in which the fibrin network is subjected to a large amplitude oscillatory shear (LAOS) deformation. The results show three dominating nonlinear features: softening over multiple deformation cycles, strain stiffening and increasing viscous dissipation during a deformation cycle. These features are incorporated in a constitutive model based on the Kelvin–Voigt model. A network state parameter is introduced that takes into account the influence of the deformation history of the network. Furthermore, in the period following the LAOS deformation, the stiffness of the networks increases which is also incorporated in the model. The influence of cross-links created by factor XIII is investigated by comparing fibrin networks that have polymerized for 1 and 2 h. A sensitivity analysis provides insights into the influence of the eight fit parameters. The model developed is able to describe the rich, time-dependent, nonlinear behavior of the fibrin network. The model is relatively simple which makes it suitable for computational simulations of blood clot formation and is general enough to be used for other materials showing similar behavior.     

The biomechanics of human femurs in axial and torsional loading: comparison of finite element analysis, human cadaveric femurs, and synthetic femurs

To assess the performance of femoral orthopedic implants, they are often attached to cadaveric femurs, and biomechanical testing is performed. To identify areas of high stress, stress shielding, and to facilitate implant redesign, these tests are often accompanied by finite element (FE) models of the bone/implant system. However, cadaveric bone suffers from wide specimen to specimen variability both in terms of bone geometry and mechanical properties, making it virtually impossible for experimental results to be reproduced. An alternative approach is to utilize synthetic femurs of standardized geometry, having material behavior approximating that of human bone, but with very small specimen to specimen variability. This approach allows for repeatable experimental results and a standard geometry for use in accompanying FE models. While the synthetic bones appear to be of appropriate geometry to simulate bone mechanical behavior, it has not, however, been established what bone quality they most resemble, i.e., osteoporotic or osteopenic versus healthy bone. Furthermore, it is also of interest to determine whether FE models of synthetic bones, with appropriate adjustments in input material properties or geometric size, could be used to simulate the mechanical behavior of a wider range of bone quality and size. To shed light on these questions, the axial and torsional stiffness of cadaveric femurs were compared to those measured on synthetic femurs. A FE model, previously validated by the authors to represent the geometry of a synthetic femur, was then used with a range of input material properties and change in geometric size, to establish whether cadaveric results could be simulated. Axial and torsional stiffnesses and rigidities were measured for 25 human cadaveric femurs (simulating poor bone stock) and three synthetic "third generation composite" femurs (3GCF) (simulating normal healthy bone stock) in the midstance orientation. The measured results were compared, under identical loading conditions, to those predicted by a previously validated three-dimensional finite element model of the 3GCF at a variety of Young's modulus values. A smaller FE model of the 3GCF was also created to examine the effects of a simple change in bone size. The 3GCF was found to be significantly stiffer (2.3 times in torsional loading, 1.7 times in axial loading) than the presently utilized cadaveric samples. Nevertheless, the FE model was able to successfully simulate both the behavior of the 3GCF, and a wide range of cadaveric bone data scatter by an appropriate adjustment of Young's modulus or geometric size. The synthetic femur had a significantly higher stiffness than the cadaveric bone samples. The finite element model provided a good estimate of upper and lower bounds for the axial and torsional stiffness of human femurs because it was effective at reproducing the geometric properties of a femur. Cadaveric bone experiments can be used to calibrate FE models' input material properties so that bones of varying quality can be simulated.     

An anatomically based patient-specific finite element model of patella articulation: towards a diagnostic tool

A 3D anatomically based patient-specific finite element (FE) model of patello-femoral (PF) articulation is presented to analyse the main features of patella biomechanics, namely, patella tracking (kinematics), quadriceps extensor forces, surface contact and internal patella stresses. The generic geometries are a subset from the model database of the International Union of Physiological Sciences (IUPS) (http://www.physiome.org.nz) Physiome Project with soft tissue derived from the widely used visible human dataset, and the bones digitised from an anatomically accurate physical model with muscle attachment information. The models are customised to patient magnetic resonance images using a variant of free-form deformation, called 'host-mesh' fitting. The continuum was solved using the governing equation of finite elasticity, with the multibody problem coupled through contact mechanics. Additional constraints such as tissue incompressibility are also imposed. Passive material properties are taken from the literature and implemented for deformable tissue with a non-linear micro-structurally based constitutive law. Bone and cartilage are implemented using a 'St-Venant Kirchoff' model suitable for rigid body rotations. The surface fibre directions have been estimated from anatomy images of cadaver muscle dissections and active muscle contraction was based on a steady-state calcium-tension relation. The 3D continuum model of muscle, tendon and bone is compared with experimental results from the literature, and surgical simulations performed to illustrate its clinical assessment capabilities (a Maquet procedure for reducing patella stresses and a vastus lateralis release for a bipartite patella). Finally, the model limitations, issues and future improvements are discussed.     

A mechanistic model for internal bone remodeling exhibits different dynamic responses in disuse and overload

Bone is a dynamic tissue which, through the process of bone remodeling in the mature skeleton, renews itself during normal function and adapts to mechanical loads. It is, therefore, important to understand the effect of remodeling on the mechanical function of bone, as well as the effect of the inherent time lag in the remodeling process. In this study, we develop a constitutive model for bone remodeling which includes a number of relevant mechanical and biological processes and use this model to address differences in the remodeling behavior as a volume element of bone is placed in disuse or overload. The remodeling parameters exhibited damped oscillatory behavior as the element was placed in disuse, with the amplitude of the oscillations increasing as the severity of disuse increased. In overload situations, the remodeling parameters exhibited critically sensitive behavior for loads beyond a threshold value. These results bear some correspondence to experimental findings, suggesting that the model may be useful when examining the importance of transient responses for bone in disuse, and for investigating the role fatigue damage removal plays in preventing or causing stress fractures. In addition, the constitutive algorithm is currently being employed in finite element simulations of bone adaptation to predict important features of the internal structure of the normal femur, as well as to study bone diseases and their treatment.     

Differences in 3D vs. 2D analysis in lumbar spinal fusion simulations

Lumbar interbody fusion is currently the gold standard in treating patients with disc degeneration or segmental instability. Despite it having been used for several decades, the non-union rate remains high. A failed fusion is frequently attributed to an inadequate mechanical environment after instrumentation. Finite element (FE) models can provide insights into the mechanics of the fusion process. Previous fusion simulations using FE models showed that the geometries and material of the cage can greatly influence the fusion outcome. However, these studies used axisymmetric models which lacked realistic spinal geometries. Therefore, different modeling approaches were evaluated to understand the bone-formation process.Three FE models of the lumbar motion segment (L4–L5) were developed: 2D, Sym-3D and Nonsym-3D. The fusion process based on existing mechano-regulation algorithms using the FE simulations to evaluate the mechanical environment was then integrated into these models. In addition, the influence of different lordotic angles (5, 10 and 15°) was investigated. The volume of newly formed bone, the axial stiffness of the whole segment and bone distribution inside and surrounding the cage were evaluated.In contrast to the Nonsym-3D, the 2D and Sym-3D models predicted excessive bone formation prior to bridging (peak values with 36 and 9% higher than in equilibrium, respectively). The 3D models predicted a more uniform bone distribution compared to the 2D model.The current results demonstrate the crucial role of the realistic 3D geometry of the lumbar motion segment in predicting bone formation after lumbar spinal fusion.     

Modified Bilston nonlinear viscoelastic model for finite element head injury studies

This paper proposes a modified nonlinear viscoelastic Bilston model (Bilston et al., 2001, Biorheol., 38, pp. 335-345). for the modeling of brain tissue constitutive properties. The modified model can be readily implemented in a commercial explicit finite element (FE) code, PamCrash. Critical parameters of the model have been determined through a series of rheological tests on porcine brain tissue samples and the time-temperature superposition (TTS) principle has been used to extend the frequency to a high region. Simulations by using PamCrash are compared with the test results. Through the use of the TTS principle, the mechanical and rheological behavior at high frequencies up to 10(4) rads may be obtained. This is important because the properties of the brain tissue at high frequencies and impact rates are especially relevant to studies of traumatic head injury. The averaged dynamic modulus ranges from 130 Pa to 1500 Pa and loss modulus ranges from 35 Pa to 800 Pa in the frequency regime studied (0.01 rads to 3700 rads). The errors between theoretical predictions and averaged relaxation test results are within 20% for strains up to 20%. The FEM simulation results are in good agreement with experimental results. The proposed model will be especially useful for application to FE analysis of the head under impact loads. More realistic analysis of head injury can be carried out by incorporating the nonlinear viscoelastic constitutive law for brain tissue into a commercial FE code.     

Enamel dictates whole tooth deformation: A finite element model study validated by a metrology method

In order to understand whole tooth behavior under load the biomechanical role of enamel and dentin has to be determined. We approach this question by comparing the deformation pattern and stiffness of intact teeth under load with the deformation pattern and stiffness of the same teeth after the enamel has been mechanically compromised by introducing a defect. FE models of intact human premolars, based on high resolution micro-CT scans, were generated and validated by in vitro electronic speckle pattern interferometry (ESPI) experiments. Once a valid FE model was established, we exploit the flexibility of the FE model to gain more insight into whole tooth function. Results show that the enamel cap is an intrinsically stiff biological structure and its morphology dictates the way a whole tooth will mechanically behave under load. The mechanical properties of the enamel cap were sufficient to mechanically maintain almost its entire stiffness function under load even when a small defect (cavity simulating caries) was introduced into its structure and breached the crown integrity. We conclude that for the most part, that enamel and not dentin dictates the mechanical behavior of the whole tooth.     

Modeling of long-term fatigue damage of soft tissue with stress softening and permanent set effects

One of the major failure modes of bioprosthetic heart valves is non-calcific structural deterioration due to fatigue of the tissue leaflets. Experimental methods to characterize tissue fatigue properties are complex and time-consuming. A constitutive fatigue model that could be calibrated by isolated material tests would be ideal for investigating the effects of more complex loading conditions. However, there is a lack of tissue fatigue damage models in the literature. To address these limitations, in this study, a phenomenological constitutive model was developed to describe the stress softening and permanent set effects of tissue subjected to long-term cyclic loading. The model was used to capture characteristic uniaxial fatigue data for glutaraldehyde-treated bovine pericardium and was then implemented into finite element software. The simulated fatigue response agreed well with the experimental data and thus demonstrates feasibility of this approach.     

A study of the viscoelastic effect in a bone remodeling model

The present paper addresses the following question can a simple regulatory bone remodeling model predict effects of viscosity on the trabecular morphology? For that, we propose an extension of a previous bone remodeling model by taking into account the viscosity properties of the tissue. Zener’s law is used to describe the mechanical behavior of the bone and a specific law of the apparent bone density rate is proposed. Based on stability analysis, numerical simulations are then performed to investigate the viscosity role on simulations of the bone remodeling process. We show that the viscous contribution affects the evolution of the apparent bone density, by slowing down the adaptation process, which seems to be confirmed by simulations with real data obtained from rat tibia.     

Intervertebral disc creep behavior assessment through an open source finite element solver

Degenerative Disc Disease (DDD) is one of the largest health problems faced worldwide, based on lost working time and associated costs. By means of this motivation, this work aims to evaluate a biomimetic Finite Element (FE) model of the Intervertebral Disc (IVD). Recent studies have emphasized the importance of an accurate biomechanical modeling of the IVD, as it is a highly complex multiphasic medium. Poroelastic models of the disc are mostly implemented in commercial finite element packages with limited access to the algorithms. Therefore, a novel poroelastic formulation implemented on a home-developed open source FE solver is briefly addressed throughout this paper. The combination of this formulation with biphasic osmotic swelling behavior is also taken into account. Numerical simulations were devoted to the analysis of the non-degenerated human lumbar IVD time-dependent behavior. The results of the tests performed for creep assessment were inside the scope of the experimental data, with a remarkable improvement of the numerical accuracy when compared with previously published results obtained with ABAQUS^®. In brief, this in-development open-source FE solver was validated with literature experimental data and aims to be a valuable tool to study the IVD biomechanics and DDD mechanisms.     

Statistical finite element method for real-time tissue mechanics analysis

The finite element (FE) method can accurately calculate tissue deformation. However, its low speed renders it ineffective for many biomedical applications involving real-time data processing. To accelerate FE analysis, we introduce a novel tissue mechanics simulation technique. This technique is suitable for real-time estimation of tissue deformation of specific organs, which is required in computer-aided diagnostic or therapeutic procedures. In this method, principal component analysis is used to describe each organ shape and its corresponding FE field for a pool of patients by a small number of weight factors. A mapping function is developed to relate the parameters of organ shape to their FE field counterpart. We show that irrespective of the complexity of the tissue's constitutive law or its loading conditions, the proposed technique is highly accurate and fast in estimating the FE field. Average deformation errors of less than 2% demonstrate the accuracy of the proposed technique.     

A new approach for assigning bone material properties from CT images into finite element models

Generation of subject-specific finite element (FE) models from computed tomography (CT) datasets is of significance for application of the FE analysis to bone structures. A great challenge that remains is the automatic assignment of bone material properties from CT Hounsfield Units into finite element models. This paper proposes a new assignment approach, in which material properties are directly assigned to each integration point. Instead of modifying the dataset of FE models, the proposed approach divides the assignment procedure into two steps: generating the data file of the image intensity of a bone in a MATLAB program and reading the file into ABAQUS via user subroutines. Its accuracy has been validated by assigning the density of a bone phantom into a FE model. The proposed approach has been applied to the FE model of a sheep tibia and its applicability tested on a variety of element types. The proposed assignment approach is simple and illustrative. It can be easily modified to fit users’ situations.