A selective smoothed finite element method for 3D explicit dynamic analysis of the human annulus fibrosus with modified composite-based constitutive model

Highlights

• A selective smoothed FEM is developed for dynamic analysis of human annulus fibrosus.

• A modified composite-based constitutive model is presented for material description.

• Compressibility of the matrix and the fiber-matrix shear interaction are considered.

• Nonlinear response of annulus fibrosus and fiber angle change are predicted accurately.

• The present method shows good robustness in large deformation of annulus fibrosus.

Abstract

The human annulus fibrosus (HAF) is the major component in response to external forces for the intervertebral disk (IVD), which maintains the stability and flexibility of human spine. It can be assumed to be an anisotropic nearly incompressible hyperelastic composite consisting of collagen fibers and matrix in the numerical simulations of biomechanics. However, due to the geometric complexity and material nonlinearity of HAF, the conventional Finite Element Method (FEM) often gets into difficulties in mesh generation and uncertainty of accuracy control. In this paper, a modified composite-based constitutive model, which considers the slight compressibility of ground substance and the shear interaction between collagen fibers and matrix, is developed to describe the mechanical behavior of HAF. In addition, based on the gradient smoothing techniques, the selective 3D-edge-based and node-based smoothed finite element method (Selective 3D-ES/NS-FEM) is developed to alleviate volume locking and improve the accuracy of linear four-node tetrahedral (TET4) elements. Combined with the modified constitutive model, the Selective 3D-ES/NS-FEM is applied into the explicit dynamic analysis of HAF undergoing large deformation. By comparing with the experiment data in the literatures and the numerical results produced by conventional FEM, the presented approach is proved to possess excellent accuracy and efficiency in predicting the nonlinear mechanical behavior of HAF, as well as the orientation change of the collagen fibers. Moreover, the Selective 3D-ES/NS-FEM is demonstrated to have robust capability in handling element distortion, even with the simplest TET4 mesh. This study is significant to the biomechanical research of HAF, and has potential value for guiding the prevention and treatment of low back pain.

Introduction

Low back pain (LBP) is a very common symptom that occurs in all age groups from children to the elderly population, and is now the leading cause of disability worldwide [1]. It is estimated that 70% to 80% of adults experience at least one episode of LBP during their lifetime [2]. Most people with LBP can recover, but reoccurrence is common and for small percentage of people the condition will become chronic and disabling [1]. Globally, the disability and costs attributed to the LBP will become an enormous social, psychological, and economic burden in coming decades. The causes of LBP are multifactorial, among which the degeneration of spine is one of the main etiologic factor. The common symptoms of spinal problems are bulging, herniated and ruptured intervertebral discs (IVD), which are closely related to the mechanical damage or degeneration of human annulus fibrosus (HAF). Therefore, the study of mechanical properties of HAF is of great significance to understand the intrinsic sources of related diseases and to guide the prevention and treatment of LBP.

In normal human spine, HAF is a tough outer wall, surrounding and protecting the nucleus pulposus in the center of the IVD. As shown in Fig. 1, the HAF is a concentric lamellar structure composed of collagen fibers and ground substance [3]. It enhances the stability and flexibility of spine during rotation and flexion, and helps the spine resist pressure. In order to have a better understanding of HAF, scholars have conducted extensive research on the mechanical behavior of HAF through experiments, theoretical analysis and numerical simulations.

Over the past several decades, numerous experiments have been carried out on the HAF specimens to explore the intrinsic mechanical properties of HAF. Wu and Yao [4] not only investigated the tensile properties of HAF specimens, but also presented the constitutive equations for HAF material. Adams and Green [5] explored the effect of sample size to the tensile behavior of HAF and found that fiber-matrix interactions make a remarkable contribution to the tensile stiffness and strength of HAF. Skaggs et al. [6] observed the tensile properties of single lamella specimens isolated from different regions of normal HAF. Fujita et al. [7] investigated the time-independent, anisotropic shear behavior of cube and sheet specimens. Their experiments demonstrated that collagen fiber interactions significantly increase the shear stiffness in the plane of the annular lamellae. Besides, Holzapfel et al. [8] performed in vitro experiments on single annulus lamellae from fresh non-degenerate human lumbar spines. They investigated the tensile response of single lamellar HAF specimens and determined the orientation of lamellar collagen fibers and their regional variation. Through these substantial experimental studies, HAF is proved to be a nonlinear anisotropic hyperelastic material, and its mechanical properties are related to the anatomical region and degeneration degree.

For the sake of theoretically describing the mechanical properties of HAF, various material constitutive models have been proposed according to the experimental results. A common approach is to treat HAF as a hyperelastic matrix embedded with aligned collagen fibers. Based on the constitutive theory for fiber-reinforced composites developed by Spencer [9], the material constitutive model is represented by a strain energy function which depends on the invariants of strain tensor and fiber directions. Elliott and Setton [10] utilized a linear fiber-induced anisotropic constitutive model allowing for compressibility to describe HAF. The results indicated that the physical effects of equivalent "matrix", "fiber", "fiber-matrix" and "fiber-fiber" interactions may be important factors affecting the mechanical properties of HAF. Besides, Elliott and Setton [11] also measured the tensile moduli and Poisson's ratios in samples of HAF along circumferential, axial, and radial directions at different regions. Furthermore, by applying the results to a linear anisotropic constitutive model, they determined a complete set of model properties and predicted the material behavior of HAF under idealized kinematic state. Bass et al. [12] performed uniaxial and biaxial tests on HAF specimens, and then developed an objective methodology for determination of strain energy function through the mechanical behavior of HAF in biaxial tension. Wagner and Lotz [13] developed a constitutive formulation to predict the response of annulus fibrosus under multi-axial tension and compression. They decomposed the strain energy function into three distinct components, including an isotropic matrix, two bundles of collagen fibers, and the interaction of fiber crosslinks. Peng et al. [14] proposed an anisotropic hyperelastic constitutive model with fiber-matrix shear interaction for HAF, and determined the material parameters by matching uniaxial tensile experimental data. Moreover, Guo et al. [15] developed a composites-based hyperelastic constitutive model for soft tissue and applied it to the material description of HAF. They assumed the matrix and fibers are both incompressible materials, and introduced fiber-matrix interaction into the strain energy function. This constitutive model has been demonstrated to be in good agreement with experimental results including fiber angle changes.

Nevertheless, due to the experimental measurement of stress in vivo HAF are difficult to operate and ethically controversial, most experiments are limited to in vitro and are performed on simple sliced specimens. In the meantime, the numerical methods are known for their flexible application and detailed observation of the object. Hence, it is necessary to introduce the numerical methods into the field of biomechanics in order to obtain more specific information of the internal force distribution. What should be noted is that an effective numerical result depends on the reliable discrete models, the appropriate material properties with correct constitutive laws, and the accurate and stable solvers. Thanks to the rapid development of computer technology, as well as the solid foundation of the previous experimental and theoretical studies on HAF materials, the numerical methods can fulfill the potential in the biomechanical simulation of HAF and further provide more detailed mechanical guidance for LBP.

As is known, Finite Element Method (FEM) is the most commonly used numerical method for simulation and analysis of complex biological systems. In recent decades, scholars have done a lot of research on biomechanical analysis of spine based on the conventional FEM. For example, Shirazi-Adl et al. [16], [17], [18] constructed a nonlinear 3D FE model considering the annulus as a composite of collagenous fibers embedded in a matrix of ground substance. Yin et al. [19] proposed a homogenization model of HAF, and then provided a comparison with the structural FE models using ABAQUS software and the fiber-reinforced strain energy models. Schroeder et al. [20] developed a 3D fiber-reinforced poroviscoelastic swelling FE model to compute the interplay of osmotic, viscous and elastic forces in IVD under uniaxial compressive load. Schmidt et al. [21,22] developed a calibration method for 3D nonlinear FE model of annulus fibrosus ring, in which collagen fibers were simulated by unidirectional spring elements. Malandrino et al. [23] regarded the lumbar annulus model as a porous hyperelastic material and considered the strain energy density terms related to the direction of fibers. Jacobs et al. [24] developed a FE model of IVD utilizing analytic geometry and modified material properties by combing the biphasic-swelling theory. Finley et al. [25] presented open-access FE models of healthy and degenerated human lumbar spine using the latest FE package FEBio proposed by their team. However, due to complex geometry, inhomogeneity and incompressibility of HAF, the conventional FEM has some common drawbacks such as difficulty in meshing with structured elements, as well as element distortion, volume locking and unstable accuracy in irregular elements when dealing with large deformation problems.

In order to address the deficiencies of conventional FEM, the Smoothed Finite Element Method (S-FEM) based on G space theory [26,27] and weakened weak (W²) form [28] was proposed by Liu and coworkers [29]. This method introduces the strain smoothing technique [30] of Galerkin meshfree method [31], [32], [33] into the framework of conventional FEM, which not only retains the advantages of FEM but also remedies most of its deficiencies [34]. So far, the family of S-FEM has been established through various constructions of smoothing domains, in which the most representative ones are cell-based smoothed FEM (CS-FEM) [35,36], node-based smoothed FEM (NS-FEM) [37,38], edge-based smoothed FEM (ES-FEM) [39,40] and face-based smoothed FEM (FS-FEM) [41,42]. They have been widely applied in the mechanical problems of solids and structures, and even in the multi-physical coupling problems, such as thermo-mechanical problems [43,44], structural-acoustic coupled problems [45,46], piezoelectric structures [47,48] and magneto-electro-elastic (MEE) structures [49,50]. Moreover, based on the hybrid types of smoothing domains, αS-FEM [51], βS-FEM [52], H-SFEM [53], ES-XFEM [54], etc., have been developed for some specific complex problems in recent years.

Substantial numerical experiments have demonstrated that NS-FEM is a volumetric locking free method due to its overly-underestimated stiffness [55,56], which is extremely valuable for the simulation of nonlinear problems like hyperelastic bio-tissues. However, because of the temporal instability of NS-FEM, it is not suitable for dynamic or static large deformation problems [57,58]. Fortunately, ES-FEM has superior h-convergence properties, computational accuracy and efficiency, as well as spatial and temporal stability [39,[59], [60], [61]]. It can greatly improve the performance of tetrahedral elements and solve the dynamic problems well, but it also suffers from the volumetric locking for simulations of nearly incompressible materials. Therefore, combining the advantages of ES-FEM and NS-FEM, the Selective 3D-ES/NS-FEM was developed and applied in analysis of nearly incompressible bio-tissues under large deformation, such as adventitial strips of artery [62], passive rabbit ventricles [63] and brain sulcus tissues [64].

In this work, the Selective 3D-ES/NS-FEM is extended for explicit dynamics analysis of nearly incompressible HAF undergoing large deformation. In Section 2, we introduce the 3D-edge-based and node-based gradient smoothing techniques for continuous solid. Then, a modified constitutive model of nearly incompressible HAF considering the fiber-matrix interaction is presented in Section 3. Employing the Selective 3D-ES/NS-FEM as solver, the smoothed strain energy and smoothed stress tensor are decoupled subsequently. In Section 4, an effective explicit dynamic algorithm and programs are developed for analysis of the HAF in nonlinear large deformation. In Section 5, the proposed approach is accordingly implemented to 3D numerical simulations of HAF sliced specimens to validate its effectiveness, accuracy and robustness. Furthermore, the intact annulus fibrosus is modelled and analyzed under uniaxial compression. In the end, some conclusions are given in Section 6.