3.1. Weak Form, Linearization and Finite Element Discretization.

Generally, a finite element formulation is established by a variational statement, which represents the weak form of a physical law. In mechanics, this can either be written as the minimization of an energy function or, alternatively, the virtual work equation can serve as the starting point [22]. The spatial form of the virtual work equation, representing the weak form of conservation of linear momentum for a deformable body, can be written as

Here, $\mathbf{\sigma }$ is the second-order Cauchy stress tensor, $\delta \mathbf{d}$ is the second-order virtual rate of deformation tensor, $\delta \mathbf{v}$ is a virtual velocity, and $v$ and $\partial v$ represent the volume and surface in the deformed configuration, respectively. Equation (1) is a highly nonlinear function of the deformation, and an iterative method is necessary to solve for the deformation. In FEBio, this equation is solved using an incremental-iterative strategy based on Newton's method, which requires the linearization of Eq. (1) [23].

The directional derivative of Eq. (1) is needed for the linearization. In an iterative procedure, the deformation $\mathbf{\varphi }$ will be approximated by a trial solution ${\mathbf{\varphi }}_k$. Linearization of the virtual work equation around this trial state gives

where $D\delta W({\mathbf{\varphi }}_k,\delta \mathbf{v})·\mathbf{u}$ is the directional derivative of the virtual work in the direction u.

In the finite element method, the domain is divided into connected subunits called finite elements. The discretization process is established by interpolating the position X of a point within a finite element in terms of the coordinates ${\mathbf{X}}_a$ of the nodes a that define the finite element and the shape functions, N _a:

where n is the number of nodes in a finite element. Quantities such as spatial coordinates, displacement, velocity, and their virtual equivalents can be interpolated similarly.

The discretization of Eq. (2) with Eq. (3) leads to the discrete form of the nonlinear finite element equations [24]:

Here, K is the stiffness matrix, which is mostly defined by the constitutive model, u is the vector of nodal displacements, and R is the residual vector, which measures the difference between internal and external forces. As the virtual velocities $\delta \mathbf{v}$ are arbitrary, they can be eliminated. An iterative solution scheme based on Newton's method can be formulated as follows:

Here, ${\mathbf{x}}_k={\varphi }_k\left(\mathbf{X}\right)$ and ${\mathbf{x}}_{k+1}$ are the current nodal coordinates at iterations k and k+1, respectively. Ideally, Newton's method achieves a quadratic convergence rate in the neighborhood of a solution. However, it also requires the formation and factorization of the stiffness matrix at every iteration, which is the most costly part of nonlinear finite element analysis from a computational standpoint. Quasi-Newton methods offer an attractive alternative. For these methods, the true stiffness matrix K is replaced by an approximation that is relatively easy to calculate. FEBio uses the BFGS method [23], which is one of the more effective quasi-Newton methods for computational solid mechanics.

3.2. Element Technology.

The 3D solid elements in FEBio are isoparametric elements, meaning that the element geometry and displacement fields on the element are interpolated with the same shape functions. Currently available elements include the linear hexahedron, pentahedron, and tetrahedron. It is well-known that these elements tend to “lock” for nearly and fully incompressible material response when using a displacement-only formulation. FEBio deals with this problem by using a three-field element formulation (the mean dilatation method) for the hexahedral and pentahedral elements [25]. For tetrahedral elements, a nodally integrated tetrahedron provides enhanced performance for finite deformation and near-incompressibility compared to the standard constant strain tetrahedron [26].

Quadrilateral and triangular shell elements in FEBio use an extensible director formulation [27]. Six degrees of freedom are assigned to each shell node: three displacement degrees of freedom and three director degrees of freedom. A three-point Gaussian quadrature rule is used for integration through the thickness of the shell.

FEBio also solves porous media problems using the mixture framework known as the biphasic theory [28]. For such problems, a fluid phase exists along with the solid phase and therefore requires the solution of a separate fluid field problem. In FEBio, a u–p formulation is used to solve the coupled solid–fluid problem for the unknown displacements and fluid pressures [29]. Like the displacements, the fluid pressures are defined at the nodes and the biphasic elements are integrated using the same quadrature rule as the structural elements. In the formulation an additional equation is needed to satisfy the conservation of mass for the mixture. The weak formulation now requires the solution of a coupled displacement-pressure problem.

and

Equation (6) is recognized as the virtual work Equation (1). In Eq. (7), $\delta p$ is a virtual fluid pressure and $\mathbf{w}$ is the flux of fluid relative to the solid, which can be calculated for instance by Darcy's law $\mathbf{w}=-\mathbf{k}\nabla p$, where $\mathbf{k}$ is the permeability tensor, and $\mathbf{n}$ is the unit outward normal to $\partial v$.

3.3. Constitutive Models.

The material representations in FEBio span the range of elastic solids, viscoelastic solids, biphasic materials, biphasic-solid materials, and rigid bodies. Many of the constitutive models are based on the concept of hyperelasticity. Constitutive models with isotropic, transversely isotropic [30], and orthotropic material symmetries are available. For most materials, uncoupled formulations are used so that these materials can be used effectively in nearly incompressible analyses. Detailed descriptions of all of the materials in FEBio can be found in the FEBio Theory Manual [21].

3.4. Rigid Bodies.

Rigid bodies are supported as a separate material type. For groups of elements, faces and/or nodes that are assigned the “rigid body” material type in FEBio, nodal degrees of freedom are condensed down to six degrees of freedom. The result is a vast reduction in the number of the total degrees of freedom and leads to a very efficient implementation of the rigid body constraint [33]. Rigid bodies can also be connected by rigid joints. The rigid joint constraints are enforced using an augmented Lagrangian method.

3.5. Linear Solvers.

As over 90% of the execution time of large problems in an implicitly integrated nonlinear FE code is spent in the linear solver (which is called repeatedly for the Newton or quasi-Newton solution method), the linear solver is typically the bottleneck for large computations. Consequently, an efficient and robust linear solver was needed. The linear solver must run on all supported platforms (Linux, Windows, Mac), take optimal advantage of the sparseness of the stiffness matrix, handle both symmetric and nonsymmetric matrices (necessary for biphasic problems), and support parallel execution. For this reason, FEBio includes support for PARDISO [34,35], SUPERLU [36], and WSMP [37]. These solvers take optimal advantage of the sparse stiffness matrix by storing only the nonzero matrix elements using the Harwell–Boeing matrix storage format. Further, they can handle nonsymmetric matrices and support parallel execution on multicore shared-memory architectures.

4.1. Uniaxial Tension-Compression Material Tests.

The single-element uniaxial stress test is a standard verification test for material models. As the deformation is homogeneous, the analytical solution can often be obtained in closed form. This allows for direct comparison of the FEBio results with the analytical solution. In this case study, results for the Mooney–Rivlin and Ogden hyperelastic constitutive models are presented. Both materials are often used to represent the matrix component of biological tissues. The Mooney–Rivlin constitutive model in FEBio has the following uncoupled strain-energy function:

where $C_1$ and $C_2$ are material parameters, K is the bulk modulus, ${\tilde{I}}_1$ and ${\tilde{I}}_2$ are the first and second deviatoric invariants of the right Cauchy–Green deformation tensor $\mathbf{C}={\mathbf{F}}^T\mathbf{F}$, F is the deformation gradient, and $J=det\left(\mathbf{F}\right)$ is the Jacobian of the deformation. The material parameters were chosen to be $C_1=6.8\mathrm{\hspace{1em}\hspace{1em}}\mathrm{MPa}$, $C_2=0$, which approximates the material parameters of articular cartilage during fast loading [40]. The value of K was chosen sufficiently high to satisfy near-incompressibility ($K~{10}^5\mathrm{\hspace{1em}\hspace{1em}}\mathrm{MPa}$).

The Ogden constitutive model with uncoupled deviatoric and dilatational response has the following strain-energy function [41]:

Here, ${\tilde{\lambda }}_i$ are the deviatoric principal stretches and $c_i$ and $m_i$ are material coefficients. Material incompressibility was enforced using the augmented Lagrangian method. Material coefficients were chosen to represent the bulk soft tissue properties of the heel-pad (N=1, $c_1=0.0329\mathrm{\hspace{1em}\hspace{1em}}\mathrm{MPa}$, $m_1=6.82$) [42]. Again, the bulk modulus K was chosen sufficiently high to enforce near-incompressibility. Results from FEBio were compared to analytical solutions. Excellent agreement was achieved between the analytical solutions and the predicted responses from FEBio for both constitutive equations (Fig. ​(Fig.11).

Open in a separate window

Fig. 1

Simulation of uniaxial tension-compression material tests for the Mooney–Rivlin (top) and Ogden (bottom) constitutive models. There was excellent agreement between the analytical solutions and the predictions from FEBio for both constitutive equations.

4.2. Confined Compression Creep.

The confined compression creep test is often used for material characterization of deformable porous media [28,43,44]. A cylindrical tissue sample is placed in a confining chamber with rigid impermeable bottom and side walls. The top of the tissue is compressed with a free-draining rigid porous indenter. If a constant applied load is used, a creep deformation is produced, increasing over time (Fig. ​(Fig.22).

Open in a separate window

Fig. 2

Confined compression creep testing of deformable porous media. (Top) Loading and boundary conditions for the confined compression creep simulation. (Bottom) Creep displacement as a function of time, calculated from the analytical solution and predicted by FEBio, showing the excellent agreement between the two results.

A unit cube was created and the nodes were constrained so that only the top surface could move along the normal direction. A pressure of 0.001 MPa was applied to the top surface. The material was assumed to be biphasic. The solid phase was modeled as a neo-Hookean material with strain energy given by

Here, $\mu$ and $\lambda$ are the Lamé parameters, $I_1=\mathrm{tr}\mathbf{C}$ is the first invariant of the right Cauchy–Green deformation tensor. For this problem, the material parameters were chosen to be $\lambda =1.43\mathrm{\hspace{1em}\hspace{1em}}\mathrm{MPa}$ and $\mu =0.357\mathrm{\hspace{1em}\hspace{1em}}\mathrm{MPa}$ and the permeability of the mixture was chosen to be ${10}^{-3}\hspace{1em}{\mathrm{mm}}^4/\mathrm{N}\hspace{1em}\mathrm{s}$. Results from FEBio were compared to the analytical solution in Mow et al. [28] for verification. The creep response predicted by FEBio was in excellent agreement with the analytical solution for the confined compression problem (Fig. ​(Fig.22).

4.3. Unconfined Compression Stress Relaxation.

In this test, a cylindrical tissue sample is exposed to a prescribed displacement in the axial direction while left free to expand radially. After loading the tissue, the displacement is held constant while the tissue under the displacement relaxes in the radial direction due to interstitial fluid flow through the material. For porous media, the outer radial boundary is free-draining (Fig. ​(Fig.3,3, arrows).

Open in a separate window

Fig. 3

Unconfined compression stress relaxation. (Top) Quarter-symmetry model used for the simulation with the free fluid flow through the outer boundaries indicated with arrows. (Bottom) Axial stress versus time calculated from the analytical solution and predicted by FEBio, showing the excellent agreement between the two results.

A quarter-symmetry mesh was used to model the cylindrical tissue sample. The outer radius and height of the cylinder were 1 mm, whereas the axial compression was 0.01 mm. The bottom of the tissue was constrained and quarter-symmetry boundary conditions were applied. The fluid pressure was constrained to zero at the outer radial surface, whereas the fluid pressure throughout the rest of the body was determined by solution of the equilibrium equations for porous media. The solid phase was represented as a neo-Hookean material (Eq. (11), $\lambda =0.186\mathrm{\hspace{1em}\hspace{1em}}\mathrm{MPa}$, $\mu =0.435\mathrm{\hspace{1em}\hspace{1em}}\mathrm{MPa}$) and the permeability of the mixture was chosen to be ${10}^{-3}{\mathrm{mm}}^4/\mathrm{N}\hspace{1em}\mathrm{s}$.

The unconfined compression response can be modeled with the biphasic theory [28]. For the special case of a cylindrical geometry and assumptions regarding the direction of the fluid flow, Armstrong et al. [45] found a closed-form analytical solution for the average axial stress on the sample in response to a step loading function. The axial stress predicted by FEBio was nearly identical to the analytical solution provided by Armstrong et al. (Fig. ​(Fig.33).

4.4. Strip Biaxial Stretching of an Elastic Sheet with a Circular Hole.

This example demonstrates the use and verifies the results of a hyperelastic material in plane stress under large deformation. A thin, initially square sheet containing a centrally located circular hole was subjected to strip biaxial stretch. The undeformed square sheet was 165×165×2 mm, with a centrally located internal hole of radius 6.35 mm (Fig. ​(Fig.4).4). The sheet was meshed with 32 hexahedral elements. This example is identical to the problem in Sec. 1.1.8 of the ABAQUS Benchmark Manual (Version 6.7).

Open in a separate window

Fig. 4

Strip biaxial stretching of an elastic sheet with a circular hole. (Top) Quarter-symmetry model used for this simulation (A) and deformed configuration after applied strain (B). (Bottom) Total nodal reaction force versus percent elongation for FEBio, NIKE3D, and ABAQUS. FEBio and NIKE3D predicted identical results that were slightly different than ABAQUS due to the different algorithms used to enforce the material incompressibility.

The sheet was stretched to a length of 1181 mm (615.76% strain), whereas the edges parallel to the displacement were restrained from contracting. The sheet was represented as a hyperelastic Mooney–Rivlin material with uncoupled deviatoric and volumetric behavior (Eq. (9)). The material coefficients were $C_1=0.1863\mathrm{\hspace{1em}\hspace{1em}}\mathrm{MPa}$, $C_2=0.00979\mathrm{\hspace{1em}\hspace{1em}}\mathrm{MPa}$, and K=100 MPa, as used by Oden [46] to match the experimental results for a rubber sheet first reported by Treloar [47]. The reaction force predicted by FEBio was compared to results from ABAQUS and NIKE3D as a function of percent elongation.

The total nodal reaction force predicted by FEBio was identical to the results produced by NIKE3D, but was slightly higher than predictions from ABAQUS (Fig. ​(Fig.4).4). The small difference is likely due to the slightly different implementation of the Mooney–Rivlin material in ABAQUS. More specifically, ABAQUS uses a different function for the dilatational strain energy (that is, the last term in Eq. (9)) than FEBio and NIKE3D.

4.5. Geometrically Nonlinear Analysis of a Cantilever Beam.

A 10 m long cantilever beam was subjected to a tip load in the transverse direction to produce over 8 m of lateral deflection (Fig. ​(Fig.5).5). The beam had a solid rectangular cross section (100 mm thick×150 mm height) and three different mesh densities were used (100, 200, and 400 elements along the length of the beam). One element was used through the thickness in all cases. One end of the beam was constrained using nodal constraints, whereas a vertical load of 269.35 N was applied at the opposite tip of the beam. The beam was represented as a hyperelastic St. Venant–Kirchhoff elastic material:

Open in a separate window

Fig. 5

Geometrically nonlinear analysis of a cantilever beam. (Top) Cantilever beam in the undeformed and deformed configurations showing the beam length and deflection. (Bottom) End displacement versus percentage of applied load for the analytical solution and FEBio simulations using 100, 200, and 400 elements along the length of the beam. The beam deflection predicted by FEBio was nearly identical to the analytical solution when 400 elements were used to discretize the beam.

where $\lambda$ and $\mu$ are the Lamé parameters (λ=0 MPa, μ=50 MPa) and E is the Green–Lagrange strain tensor. The use of the Green–Lagrange strain tensor in the formulation of the St. Venant–Kirchhoff constitutive equation makes it appropriate for the large rotations and small strains in this problem. The cantilever tip deflection predicted by FEBio was compared to the analytical solution for the large deflection of a cantilever beam [48].

The beam deflection predicted by FEBio was nearly identical to the analytical solution for the geometrically nonlinear analysis of a cantilever beam for a model with 400 elements along the length of the beam (Fig. ​(Fig.5).5). Models with 100 and 200 elements produced a slightly stiffer response.

4.6. Twisted Ribbon Test for Shells.

This is one of a number of verification problems that are used to verify the shell element formulation in FEBio. A plate discretized with shell elements was twisted by constraining one end and applying equal but opposite forces to the corners of the nonconstrained end (Fig. ​(Fig.6).6). The problem description and analytical solution were originally presented by Batoz [49]. The plate dimensions were 1.00×2.00×0.05 mm and five models with mesh densities from 32 to 2048 elements were used for a convergence study. Nodes along one short edge were clamped using nodal constraints and the plate was twisted by applying equal but opposite forces of 1 N to the corners nodes of the opposite short edge. The St. Venant–Kirchhoff material was used (Eq. (12), λ=2.5×10⁷ MPa, μ=2×10⁶ MPa). The displacement of the corner node predicted by FEBio was compared to the analytical solution.

Open in a separate window

Fig. 6

Twisted ribbon test for shells. (Top) Model in the undeformed configuration with loading and boundary conditions indicated. (Bottom) Displacement of the outside corner node versus number of elements used to discretize the model. 2048 elements were needed to exactly match the analytical solution, but there was <2% difference when only half that many elements were used to discretize the model.

The corner displacement predicted by FEBio was nearly identical to the analytical solution for the twisted ribbon test for shells when 2048 elements were used to mesh the plate (Fig. ​(Fig.6).6). To put this result in perspective, the corner displacement predicted by FEBio was <2% different when half that many elements were used and ~30% different when only 32 elements were used.

4.7. Upsetting of an Elastic Billet.

This problem is commonly used to test finite element approaches for enforcing near- and full-incompressibility, as the domain is highly constrained. A billet of material was compressed between two rigid flat surfaces (Fig. 7(A)). As the material was only compressed to the point when the material bulging from the middle of the billet makes contact with the flat surfaces (Fig. 7(B)), this problem is a simplification of the contact problem presented in the next example (Fig. 7(C)).

Open in a separate window

Fig. 7

Upsetting of an elastic billet. (Top) Quarter-symmetry model in the initial configuration (A); billet deformed until there is contact between the inner material of the billet and the rigid planes (B); and deformed billet with contact between the inner material and the rigid planes (C). (Bottom) Lateral displacement (bulge) versus percent compression for FEBio, NIKE3D, and ABAQUS. There was excellent agreement in the predicted results from all three FE codes for both the contact and noncontact versions of this problem.

A quarter-symmetry, plane strain model was produced (Fig. 7(A)). The model dimensions were 1.0×1.0×0.1 mm, and a mesh of 10×10×1 elements produced converged results. Boundary and loading conditions were prescribed using nodal constraints and nodal loads, respectively. The billet was represented as an uncoupled Mooney–Rivlin material (Eq. (9), C ₁=1 MPa, C ₂=10MPa, K=10,000MPa). The augmented Lagrangian method was used to enforce incompressibility [32]. The displacement of the top right node of the mesh (representing the amount of lateral bulging) was plotted versus the percent compression and the results from FEBio were compared to results from ABAQUS and NIKE3D.

The lateral displacement of the billet predicted by FEBio was nearly identical to the results produced by ABAQUS and NIKE3D (Fig. ​(Fig.7,7, right panel).

To verify the contact algorithm implemented in FEBio, contact between the incompressible material of the billet and the rigid plane was added to the problem. The billet was compressed 60% of its original height, causing extensive contact between the rigid plane and the material originally in the middle of the billet. The boundary conditions for this problem were the same as the previous problem, but here the compression of the billet was caused by frictionless contact between the billet and rigid plane. This contact was enforced using an augmented-Lagrangian method. Again, the lateral displacement of the billet predicted by FEBio was nearly identical to the results produced by ABAQUS and NIKE3D. The results from the noncontact and contact problems lined up through the smaller range of compression (Fig. ​(Fig.7,7, right panel).

4.8. Crushing of a Pipe.

This problem simulated the crushing of a long, straight pipe of cylindrical cross section between two flat, frictionless anvils. This benchmark problem includes large strains/rotations and contact between deformable and rigid materials. A quarter-symmetry, plane strain model of a 114.3 mm radius pipe with an 8.87 mm thickness was created for this problem (Fig. ​(Fig.8).8). The model was meshed with 24 elements along the arc, four elements through the thickness, and one element through the length. Nodal constraints were used for the plane-strain and quarter-symmetry boundary conditions, and the problem was driven by prescribing a 50 mm displacement to the rigid body. Contact between the anvil and pipe was enforced using the augmented Lagrangian method. The St. Venant–Kirchhoff material was used for this problem ($\lambda =107\mathrm{\hspace{1em}\hspace{1em}}\mathrm{GPa}$, $\mu =71.5\mathrm{\hspace{1em}\hspace{1em}}\mathrm{GPa}$). The rigid body reaction force predicted by FEBio was compared to the results from NIKE3D and ABAQUS as a function of anvil displacement.

Open in a separate window

Fig. 8

Crushing of a pipe. (Top, left) Quarter-symmetry model with rigid crushing planes and key dimensions shown. (Top, right) Pipe in the final crushed configuration. (Bottom) Rigid body reaction force versus percent applied displacement for the crushing planes. There was less than a 3% difference in the peak forces predicted by the three FE codes.

The results predicted by the three FE codes are shown in Fig. ​Fig.8,8, demonstrating good agreement among all three codes. Small observed differences were most likely due to the different algorithms to enforce contact. Still, at peak displacement there was less than 3% difference in the predicted forces.

4.9. Cartilage Layer Compressed by a Flat, Rigid, Impermeable Surface.

A half-symmetry finite element frictionless contact analysis was performed between a spherical elastic layer anchored to a rigid impermeable substrate and a flat impermeable rigid surface (Fig. ​(Fig.9).9). The geometry was representative of the articular layer of an immature bovine humeral head, with a cartilage surface radius of 46.3 mm and a cartilage layer thickness of 0.8 mm. The converged mesh consisted of 20 tri-linear hexahedral elements through the thickness, 50 along the radial direction, and 14 along the circumferential direction. The top of the cartilage layer was constrained by defining the top layer of nodes as a rigid body with all 6 degrees of freedom constrained. The deformation at the center of the articular layer was set to 0.095 mm (~12% of the thickness) by displacing the bottom nodes using nodal prescribed displacements. The cartilage layer was represented by the Mooney–Rivlin constitutive model in Eq. (9) to approximate the short-time response of the tissue, again using $C_1=6.8\mathrm{\hspace{1em}\hspace{1em}}\mathrm{MPa}$ and $C_2=0$. The radial and circumferential stresses predicted by FEBio through the thickness of the middle of the cartilage layer were compared to the results from ABAQUS and NIKE3D.

Open in a separate window

Fig. 9

Cartilage layer compressed by a flat, rigid, impermeable surface. (Top) Schematic of the model illustrating the boundary and loading conditions. (Middle) Radial stress versus percent thickness through the middle of the cartilage layer predicted by FEBio, NIKE3D, and ABAQUS. (Bottom) Circumferential stress versus percent thickness predicted by the three FE codes. The stresses predicted by FEBio and NIKE3D were nearly identical and the stresses predicted by ABAQUS were less than 3% different than the other codes.

The predicted radial and circumferential stresses through the thickness of the middle of the cartilage layer were the same for FEBio and NIKE3D, but the stresses predicted by ABAQUS were consistently slightly larger (Fig. ​(Fig.9).9). The small differences (<3%) are most likely due to differences in the algorithms used to enforce the contact constraint: FEBio and NIKE3D use an augmented Lagrangian method to enforce contact, whereas ABAQUS uses a Lagrange multiplier method.

Financial support from the National Institutes of Health (1R01GM083925) is gratefully acknowledged. The authors thank Dr. Ahmet Erdemir (Cleveland Clinic Foundation) and Dr. Alexander Veress (University of Washington) for their feedback related to the FEBio software development.