Verification Process for Finite Element Modeling Techniques Used in Biological Hard Tissue

Patient-specific finite element modeling and simulation have great potential for assisting orthopedic surgeons in planning procedures, predicting risk of fracture, and individualizing surgical implants. However, the field has been severely limited by the use of computationally inaccurate techniques without suitable calculation convergence verification, impractical construction and execution time, and the need for experienced modeling personnel. Standardization of computational verification and validation has evolved greatly in the past decade, maturing in the assessment of finite element model accuracies [1,2]. This work focused solely on the fractional calculation convergence verification of finite element modeling techniques for hard biological tissues.

This study systematically compares calculation convergence behavior of two high-resolution medical image-based finite element modeling techniques in use in modeling biological hard tissues. Model maximum strain and equivalent model stiffness were compared for hard biological tissue in this case for a lumbar vertebral model to investigate hp-convergence and model precision. H-convergence alone refines the subdomain finite element size while p-convergence increases the order of approximations within an element without altering the mesh size [3]. The combination of the two abovementioned convergence techniques is known as hp-convergence in which both mesh size and element order are refined.

Two forerunning methods (voxel-based and geometry-based) exist for generating patient-specific finite element models from medical images of high-resolution computed tomography (CT) or magnetic resonance imaging (MRI) scans, but the convergence of these two modeling techniques have never been compared for a population of vertebral geometries. Advances in CT and MRI scan resolutions and accurate descriptions of the mechanics of biological materials have greatly improved the accuracy of patient-specific finite element modeling. Once a physical system is modeled, the accurate finite element simulation result is found from a convergence study based on refining mesh and/or increasing order, which shows the numerical solution approaching asymptotically to a value. Without a properly conducted convergence study, reported results of the model of the physical system under investigation are not a numerically accurate solution. In a convergence study, the same geometrical model is meshed with varying degrees of refinement and order, by changing mesh seed length (h-convergence) and/or changing element orders (p-convergence). A whole model quantity of interest (QOI), such as strain energy, is examined for each model tested. As the accuracy of spatial discretization increases, the percent change between QOI values should decrease to an acceptable level as the solution for the model asymptotically converges to the mesh and/or order independent solution for the system under investigation. This, referred to as fractional change, is a commonly used method of assessing model convergence [2]. Following model convergence, it is then appropriate to validate the solution through a comparison with experimental results, proving that the predicted model solutions are physiologically realistic.

The geometry-based modeling technique begins with the segmentation of CT or MRI slices of the biological tissue in question and reconstruction into a three-dimensional geometrical model. Smoothing algorithms are used to reduce the effect of slice thickness on the realism of a model [4–6]. The resulting geometry can then be discretized into finite elements using a meshing algorithm. During mesh generation, the elements are not constrained to following the global orientations, resulting in a smooth approximation of the shape. In this method, elements can be any shape and any order. Mesh refinement can occur by decreasing the mesh seed length or by increasing the element order [3,7–9].

Alternatively, the voxel-based modeling technique uses computer algorithms to automatically convert a three-dimensional medical image into a finite element model [10–17]. Spacing between slices generates a voxel, and custom algorithms convert the voxels directly into finite elements. The resulting model always uses hexahedral elements aligned with the global orientations. Mesh seed length, limited to the voxel size and the slice thickness, conventionally falls between 1 and 3 mm [18]. As improving scan resolution requires increased scan time and better scanning equipment; mesh refinement has rarely been investigated and is only possible through an increase in the element order [19].

The accuracy of the geometry-based and voxel-based finite element methods has been compared for simplified geometries, macroscale hard biological tissue, and bone microstructure [5,20,21]. The solution accuracy of the voxel-based method is degraded at domain boundaries due to the jagged boundary requiring an averaging scheme to achieve reasonable accuracy [20]. Therefore, convergence of voxel-based models requires a very refined mesh that can be computationally expensive [21]. Alternatively, an accurate solution can be obtained with a geometry-based model with significantly less elements since the domain boundary is less sensitive to the mesh density. Geometry-based models have also been found to result in more accurate material property characterization in finite element models of bone trabecular structure.

In this study the two aforementioned modeling techniques were systematically compared. All geometry-based numerical test models hp-converged at an acceptable mesh seed length of 0.6 mm, while not all voxel-based models exhibited convergence and no voxel-based models p-converged. Converged geometry-based meshes were found to offer accurate and consistent solutions of the deformed model shape and equivalent vertebral stiffness, while voxel-based models were 6.35% ± 0.84% less stiff (p < 0.0001) and deformed 6.79% ± 0.96% more (p < 0.0001). Element shape and order had no significant effect on model predictions. Based on this study, we propose that convergence-based verification approaches for biological tissue and simulations to be followed to achieve accurate prediction results. In this study, fractional convergence is used to quantify model convergence behavior, but more rigorous convergence checks, such as Roache's grid convergence index [2], are also acceptable convergence criteria metrics. Findings from this study will help as a guide in selecting an appropriate finite element modeling verification technique to be used in biomedical tissue response studies.

Geometrical models of seven L4 vertebral bodies were used in this study. A literature survey was conducted to find relevant vertebral body measurements to scale the models to match patient-specific vertebral dimensions, as summarized in Table 1 and pictured in Fig. 1 [22–25]. Literature values were selected such that each study had both male and female dimensions and no vertebral body abnormalities. Studies were chosen so that oldest test subject was under 35, to minimize the effect of aging. The use of geometrical models equivalent to the anatomical vertebrae was necessary to eliminate potential geometrical discrepancies introduced by fundamental differences between modeling techniques. The models were created with high-fidelity, using accurate material models, appropriate geometrical detail, and realism in the loading and boundary conditions. The experimental setup was designed to simulate ex vivo axial compressive studies of excised cadaveric specimens—an experimental procedure commonly simulated with both voxel- and geometry-based modeling techniques [18,26,27]. To effectively load and bound the models, a cut was generated on the cranial and caudal surface of the control vertebra. All nodes on the caudal surface of each control vertebra were fixed to prevent translational displacement, simulating the common test practice of embedding the vertebra in resin. The entire cranial surface of each control vertebra was subjected to a compressive pressure of 300 N distributed across the cranial surface, applied linearly over 10 s. This resulted in a low strain rate loading condition, with an approximate strain rate of 0.0001/s. Orientation conventions used in this study are included in Fig. 1.

Voxel-based meshes were generated using conventional voxel-model development techniques [10–17]—through conversion of a voxel into a finite element when over 50% of the voxel contains the material in question. An example of the resulting voxel-based models is included in Fig. 2. The geometry-based meshes were generated using the Abaqus default meshing algorithm.

To simulate the relative mechanical contributions of cortical and cancellous material properties, the cortical and cancellous material properties were averaged using the arithmetic mean with respect to the relative thickness of the cortical bone and the cancellous bone in the model's outer elements. This technique is similar to the procedures used to assign material properties in voxel-based modeling [18]. As scan resolution decreases, the voxel size will increase, increasing the volume of tissue that is represented within a voxel. As the brightness of the tissue is correlated with the relative density of the underlying tissue, the mass of the tissue and the volume of the voxel will combine in a sum, making the arithmetic mean ideal for representing the combined mechanical contributions of the tissues. A cortex thickness of 0.5 mm was assumed [30–32].

All simulations were performed with Abaqus v6.12-1 with two Intel^® Xeon^® 2.3 GHz CPUs with 128 GB RAM on 22 cores on a Windows x64 operating system. Numerical accuracy of the installation of the solver was confirmed with Abaqus Verification v6.12-1 for the purposes of software quality assurance. To evaluate convergence Behavior of the models, test subjects were generated from four sets of independent variables: meshing technique (geometry-based, voxel-based), element order (linear, quadratic), element shape (hexahedral, tetrahedral), and mesh seed length (0.6 mm, 0.75 mm, 1 mm, 1.5 mm, 3 mm) [18]. The Cartesian product of these sets was a set of 40 test cases to be investigated, as seen in Fig. 3. The final whole-model internal strain energy was reported for each test subject to evaluate convergence of the models to a monotonic solution. As experimental measurements to the simulations were not available, validity of the solution was not evaluated, and only the convergence characteristics of the solutions were examined. The acceptable convergence limit was defined as a  ± 5% change in internal strain energy between mesh seed lengths and element orders [33,34]. The most converged model was defined as the model that exhibited the smallest percent differences between adjacent mesh seed densities and element orders. The term false convergence was defined as the case when apparent acceptable convergence levels were observed at lower fidelity meshes, but further mesh refinement did not maintain a converged percent change. It is important to note that this definition of convergence is not constant across the literature. This definition of convergence offers a measure of the stability of solution predictions and not validation of the solution. As there is no experimental validation dataset available for this work, the ability for the modeling technique to produce a consistent solution was reported. Fractional change acts as an indicator that the solution has stably converged and does not measure the accuracy of the converged solution. Therefore, the accuracy of the converging solution was reported as the difference between the converging model and the most converged model.

Maximum compressive axial strain $ε$ was calculated to compare final deformed states. The applied load was divided by maximum axial displacement (normalized by patient height) to give the whole model equivalent stiffness, $k$⁠. Axial compressive stress was examined in the final deformed state as an indicator for potential mode of fracture. The midline axial compressive stress $σc$ was defined as the mean axial compressive stress value of all nodes located in the outer shell of elements at the midline of the model, as depicted in Fig. 4. This measurement location was selected as the common fracture location in low strain rate compression of vertebrae is at the midline. Computational time was defined as the elapsed time between submission and completion of the simulation.

As all the models tested in this study were computational, there were no physical test vertebrae to create experimental data to evaluate the predicted solutions of the models tested. As the predicted values by the most converged model (the quadratic hexahedral model at 0.6 mm) most closely resemble the actual system results, the QOIs of each test subject's most converged model were used as the accepted QOIs in calculations. To compare the quality of QOI results, percent changes were compared between the most hp-converged model for each test subject and all other h- or hp-converged models, with a commonly used error value of  ± 5%. Normalized values of the QOI were created by dividing each value by the theoretically most converged value for each test subject.

Student's t-test was used to compare $ε$⁠, $k$⁠, $σc$⁠, and computational time of each test group. Six test groups consisting of four models were compared for each test subject to investigate three independent variables: modeling technique (geometry-based, voxel-based), element order (linear, quadratic), and element shape (tetrahedral, hexahedral). Significance was based on a 95% confidence interval and p-values were reported.

The stability of the solution was investigated by comparing neighboring mesh seed densities. H-convergence was achieved by the smallest mesh seed length investigated (0.6 mm) for the geometry-based models, and the linear hexahedral geometry-based models were converged by the 0.75 mm mesh seed length. An acceptable level of h-convergence was not achieved by the smallest mesh seed investigated (0.6 mm) in all of the voxel-based models. False convergence was observed in the voxel-based model of one test subject and the voxel-based model of one test subject failed to converge. Percent change values between neighboring mesh seed densities are presented in Table 3. The accuracy of the converged solution was investigated by examining the percent change with respect to the most converged geometry-based model with the same element shape and order. The most accurate models were the 0.75 mm geometry-based models. The solution error increased with increasing mesh seed length. The voxel-based models, in general, exhibited a lower strain energy. Comparable strain energies were observed between higher mesh seed lengths and the 0.6 mm geometry-based model equivalent, when the solution had not stably reached convergence.

When investigating stability of the solution, p-convergence was observed when comparing the same mesh seed length of geometry-based models. The hexahedral and tetrahedral geometry-based models reached acceptable convergence levels at 1 mm and 0.75 mm, respectively. P-convergence was not observed for any voxel-based models. The geometry-based models also exhibited a more accurate solution, with stable solutions seen at 0.75 mm for all hexahedral geometry-based models. Again, the solution error increased with the mesh seed length, following the pattern of solution stability. All percent change values between element orders are given in Table 4. Convergence characteristics of different element types are depicted in Fig. 5 by plotting the internal strain energy by the total number of elements for each test subject.

Significant differences were observed between computation time for linear and quadratic models in all seven test subjects. Additionally, significant differences between voxel- and geometry-based models were observed in all seven test subjects in $ε$ and in six test subjects in $k$⁠. No significant difference in computational time was observed between voxel- and geometry-based models. Significance was not observed in $σc$⁠. p-values are included in Table 5.

For each test subject, the converged models were isolated and compared to determine the difference in QOI outputs. All geometry-based models that were hp-converged acceptably and all voxel-based models that were h-converged were isolated. As p-convergence was not reached by voxel-based models, h-convergence was deemed acceptable. This resulted in a subset of 24 test cases, detailed in Table 6. The most converged model was closest to the mesh size and element order independent result and used as a base of comparison. The percent changes in $ε$⁠, $k$⁠, and $σc$ were compared.

Voxel-based models over-predicted $ε$ by 6.79% ± 0.96%, while geometry-based models under-predicted $ε$ by an acceptable 0.77% ± 0.66% (p < 0.0001). Voxel-based models over-predicted $k$ by 6.35% ± 0.87%, while geometry-based models acceptably under-predicted $k$ by 0.77% ± 0.67% (p < 0.0001). The $σc$ varied with mesh seed density. When comparing normalized mean midline axial stress values $σc_norm$⁠, significant differences were observed between all mesh seed densities (p < 0.0001). Therefore, only mesh seed lengths of 0.6 mm were compared to determine percent change in $σc$⁠, resulting in 16 test cases. In this smaller subset, voxel-based models unacceptably under-predicted $σc$ by 5.01% ± 0.42% and geometry-based models acceptably under-predicted QOI results by 0.02% ± 0.10% (p < 0.0001). Normalized QOI values for maximum compressive axial strain $εnorm$⁠, whole model equivalent stiffness $knorm$⁠, and mean midline axial stress $σc_norm$ were compared for voxel- and geometry-based models in Fig. 6.

Although this study mainly focused on fractional calculation verification, we discuss validation for completeness. As the control vertebrae are theoretical, no experimental validation is possible for this study, rendering a controlled computational comparison necessary. However, the realistic and physiologically accurate generation of the geometrical vertebral models and the simulation of common experimental test conditions allowed for comparison between QOIs and experimental results. The average computed stiffness $k$ of 18.6  ±  3.43 kN/mm from this study is consistent with finite element predictions for $k$ in lumbar vertebrae of 5–27.5 kN/m [35,36]. Physiological values of other investigated QOIs are not reported in literature with enough accuracy to enable further validation of the model.

In most test cases, the h-converged mesh seed length was smaller than the conventional mesh seed lengths of 1 mm, 1.5 mm, or 3 mm found in literature in the modeling of vertebrae [18,37,38]. The only exception is convergence being achieved by 1 mm in the voxel-based model of one test subject. This convergence trend at literature observed mesh seed lengths was not observed in other test subjects. In most finite element studies of biological hard tissue, a convergence study is not presented in the published results. However, without the presentation of proper convergence verification for these models, the results reported by studies with larger mesh seed lengths offer no proof that reported results are comparable to physiological responses of biological tissues. It is strongly recommended that future studies properly confirm the numerical accuracy of the model predictions when presenting their results.

The absence of acceptable levels of p-convergence in voxel-based models, regardless of element type, implies that this technique may require further refinement and extensive investigation to ensure that reported conclusions are accurate. This observation, coupled with the observed p-convergence by geometry-based models, indicates that nonconvergence is likely due to incongruities in model shape. Mesh refinement in geometry-based modeling preserved the model's shape much better than voxel-based modeling, resulting in smooth convergence. Additionally, quadratic element order resulted in the best approximation of the vertebral shape, resulting in a smoother convergence than linear models. As expected, hexahedral elements resulted in a solution with, theoretically, higher numerical accuracy than tetrahedral elements. This allowed for smooth and swift convergence. Voxel-based mesh refinement occurs through a decrease in CT or MRI image spacing, resulting in small changes in overall model shape. Model geometrical inconsistencies, such as those created by voxel-based modeling, have been shown to inhibit the strain energy convergence of finite element models of bone [39]. Therefore, convergence of voxel-based models must be confirmed with local stress values and validation of predicted QOIs before finite element predictions are to be trusted. In addition, voxel-based models introduce artificial corners on domain boundaries causing pseudo-stress/strain singularities that cannot be eliminated with mesh size refinement, element shape, or element order increase.

Significant differences were seen between voxel- and geometry-based models in both $ε$ and $k$⁠. The difference in overall model shapes was anticipated to contribute to differences in the deformed shape of the model. Theoretically, as the characteristic length of voxel elements approaches zero, the shape of the model will more closely reflect the shape of the vertebral body. As the number of elements increases and mesh seed length decreases, convergence of the solution indicates that the error associated with deformation of the model decreases to an acceptable level. The significant difference in $ε$ and $k$ highlights the importance of convergence when comparing solutions.

Quantity of interest solution values were compared amongst converged models. Unacceptable percent differences in solutions between voxel- and geometry-based models were observed in all QOIs investigated. There were significant differences in terms of $σc$ values between geometry- and voxel-based models, despite level of convergence. When comparing the $σc$ of each test subject, all test subjects with acceptably converged voxel- and geometry-models, all geometry-based models predicted greater axial stress at the midline of the model. This is an anticipated result, as these changes are due to the fundamental differences in shapes of the voxel- and geometry-based models. The low-resolution of voxel-based modeling results in the midsection of the vertebrae being represented as a column rather than approximating the natural curvature of a healthy vertebra, acting like a short compressive member and not recognizing the stress concentrations that exist at the vertebral midline. Although beyond the scope of this study, it is likely that the compressive stresses within the voxel-based models show greater deviation from the geometry-based models away from the midline. Due to the plane of symmetry that exists at the midline, the singularities in the voxel-based models due to reentrant trihedral and dihedral angles vanish. Therefore, it is probable that deviations between the voxel and geometry-based models will be smallest at the midline and increase away from the plane of symmetry. Additional investigation into the differences in compressive stress distributions between appropriately converged voxel- and geometry-based models is recommended before investigation of alternative fracture mechanisms.

This has definite implications in the use of voxel-based modeling in the simulation of fracture mechanics. The under-prediction of stress values by voxel-based models imply that they are not as adept at modeling compressive Behavior or predicting fracture sites. Therefore, it is recommended that geometry-based models be favored for predicting fracture mode or vertebral strength. The voxel-based modeling technique, if used, should be created with high-resolution and small mesh seed lengths.

Significance was observed in differences between mesh seed density and the $σc$ of all models, regardless of the level of convergence. This phenomenon, observed in this study, is attributed to the method of simulating cortical bone and highlights the current weaknesses of artificially imposing a continuous cortical shell. The cortical bone in a vertebral body is thin, discontinuous, and seamlessly transitions into peripheral trabeculae of the cancellous core. This contributes to how vertebrae carry and transmit load. The cortical shell when excised from the cancellous core is thin and discontinuous, relying on peripheral trabeculae to transmit the load to the cancellous core and to bridge gaps in the shell [31]. However, traditional CT-scan voxel size is not fine enough to capture these discontinuities with enough fidelity to reflect their inclusion in finite element models. The work presented here has an idealized perfectly continuous cortical shell, resulting in the cortical bone stress shielding the cancellous bone and an inaccurate stress distribution. As the thickness of the outer shell decreased and the density of the outer shell of elements approached pure cortical bone, the outer shell of elements was responsible for carrying more of the load.

Limitations of this study include the lack of solution convergence seen or verified for voxel-based models, making all predicted values and trends in this paper numerically less accurate solutions—a converged model would more reliably predict model behavior. Finally, using the most converged model in determining percent errors introduced inaccuracies into the results. A finite element model is an approximation of real materials and geometry and, therefore, some degree of inaccuracy is inevitable during modeling. Mesh- or/and order-based convergence studies assess the numerical accuracy of the solution and must always be conducted, but biological experimental studies should be used for validation whenever possible.

Results from this study may help researchers as a guide in fractional calculation verification process and selection of appropriate finite element modeling techniques to be used in biomedical tissue response studies. This study reasons for the use of voxel-based modeling with care and precision, with multiple QOIs used to confirm both h- and p-convergence. Maximum compressive axial strain and whole model equivalent stiffness can be reported with any appropriately converged geometry-based model with confidence. Modeling of fracture modes is best accomplished through the use of geometry-based modeling, and when voxel-based models are used, it should be verified that the mesh is fine enough that they adequately simulate the natural vertebral curvature. Additionally, the artificial superposition of an ideal cortical shell stress shields the cancellous core and patient-specific bone mineral density is the best method of capturing the micromechanical relationship between cortical and cancellous bone in vertebrae.

• U.S. National Science Foundation (NSF) (Grant No. 1148897; Funder ID: 10.13039/100000001).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

$Cijkl$ =

elastic stiffness matrix

CT =

computed tomography

$Ep$ =

in-plane elastic modulus

$Et$ =

transverse elastic modulus

$Gtp$ =

transverse shear modulus

$k$ =

whole model equivalent stiffness

$knorm$ =

normalized whole model equivalent stiffness

MRI =

magnetic resonance imaging

Greek Symbols

$γ$ =

parameter gamma

$εkl$ =

strain

$εnorm$ =

normalized strain

$υp$ =

in-plane Poisson's ratio

$υpt$ =

transverse Poisson's ratio

$σc$ =

midline axial compressive stress

$σc_norm$ =

normalized midline axial compressive stress

$σij$ =

stress