Verification of Bone Strength Pipeline in Abaqus

Final Year Dissertation Project

Introduction

Performed a thorough literature review, assessing the impact of Osteoporosis (especially in elderly women), the use of computational (FEA) models to predict femur strength, and the format and conversion of VTK files to Abaqus INP files.

Successfully developed an Abaqus-based femur strength prediction Finite Element model, based on the Linear model obtained from existing literature, that achieves within 1% Minimum Femur Strength of the existing Ansys-based model for 15 loading orientations.

Developed the “VTK to Abaqus INP Converter and Reader” and “Results Post Processor” applications, using Python and relevant libraries, to automate and streamline the femur strength pre and post-processing to predict the MFS. Lays a firm foundation for future automation and integration with other applications.

The “VTK to Abaqus INP Converter and Reader” application: converts surface, volume, and volume mesh with material properties VTK files to corresponding Abaqus INP files, identifies the element type of input VTK file, and filters nodes based on various input parameters.

The “Results Post Processor” application: calculates the average minimum principal strain for each node from values yielded from adjacent elements and determines the average minimum principal strain and Minimum Femur Strength when considering each node as the center node of a 3mm radius using the linear force-strain relationship.

Brief:

Osteoporotic hip fractures are a major health concern, especially for women over 45, with side falls being the cause of 95% of such fractures. Numerous Finite-Element femur strength studies exist simulating a side-fall configuration for an older woman’s femur. This research project set out to verify an existing Ansys-based femur strength analysis in Abaqus, focusing on an 80-year-old woman’s femur under a side-fall condition. The project sought to ensure software enhancements, cross-platform accessibility, and alignment with industry requirements.

To achieve these goals, ASCII VTK files were converted to corresponding Abaqus INP files using Python scripts, and boundary conditions were applied. A Linear model was selected due to its extensive usage in previous studies and its ability to minimize potential discrepancies from other models. The primary objective was to ensure that the Abaqus-predicted Minimum Femur Strength (MFS) for the neutral orientation was within 1% of the Ansys-based counterpart.

The Abaqus-based femur strength model predicted the MFS as 2162.7N, which varied by only 0.4% from the Ansys-based counterpart in the neutral loading orientation. Fourteen other loading orientations achieved an average MFS of 2394N with an average percentage difference of nearly 0.1. However, thirteen loading orientations failed to reach the 1% MFS verification threshold, mainly due to inaccuracies in the boundary condition at the greater trochanter.

To streamline the femur strength analysis process, two Python-based applications, “VTK to Abaqus INP Converter and Reader” and “Results Post Processor“, were developed. These tools automated both the pre-processing and results post-processing, thus enabling efficient MFS predictions. The project’s findings lay a firm foundation for future work, with the potential for further automation or integration with other applications.

Aim:

Develop a femur strength analysis of an 80-year-old woman under side-fall conditions in Abaqus and verify based on the existing Ansys pipeline.

Objectives:

Covert an existing Ansys-based femur strength finite element analysis (FEA) pipeline for an 80- year-old woman’s femur under side-fall configuration to an Abaqus-based pipeline through the following process:

L1. Convert surface mesh (SM), volume mesh (VM), and volume mesh with material properties (VMMP) VTK files to Abaqus format INP files.

L2. Apply relative boundary and loading conditions on the respective node(s) for the chosen FE modelling approach.

Perform sanity checks for L1 and L2 to ensure the methods’ accuracy.

 L3. Select an appropriate failure criterion and implement results post-processing methods to predict the MFS.

Ensure that the newly developed Abaqus-based model achieves an MFS within 1% of the existing Ansys-based model.

Hip fracture is a severe and costly injury affecting primarily older adults over 65, especially women. Approximately 76,000 hip fractures occur in the UK each year. It is typically a consequence of a weak bone structure (known as Osteoporosis) and a side fall and poses a greater risk to people of Asian or Caucasian descent. One in five men and one in three women over fifty will suffer fractures related to osteoporotic conditions.

Osteoporosis is a health condition that significantly weakens bones making them brittle and prone to fracture due to a decline in bone mass (the total amount of bone substance), hence becoming less dense and more vulnerable to fracture. The condition develops slowly over several years. Osteopenia is the stage before Osteoporosis, where the bone density obtained from bone density scans is lower than average for that age, gender, and ethnicity but not low enough to be classed as Osteoporosis. The condition develops slowly over several years and is usually identified upon sudden impacts such as side falls.

Osteoporosis is currently diagnosed by determining bone mineral density (BMD) using dual-energy X-ray (DXA) scans. The determined BMD value is compared with the average value for healthy females. Osteoporosis is defined when the BMD value equals or exceeds 2.5 standard deviations. Once diagnosed, the fracture risk is estimated using the BMD value and other parameters involving complex instruments. Although commonly used, this method (using the BMD value) predicts fractures at a rate of around 40% with a false positive of around 20%.

Computed tomography-based finite element methods have been increasingly used over the past four decades to simulate the mechanical behaviour of bones (predominantly adult bones). Given that 95% of Osteoporotic hip fractures are due to side falls, numerous femur strength FE studies exist that simulate an elderly woman’s femur under a side-fall configuration.

An existing Ansys model simulates the side-fall configuration of an 80-year-old woman’s femur. Converting the Ansys femur strength pipeline to an Abaqus-based pipeline offers several benefits. It ensures compatibility with software advancements and enhances accessibility for devices, users, and platforms where Ansys is not the primary FEA software. Developing and verifying an Abaqus-based femur strength pipeline will also enable SIMULIA to introduce software updates that align with industry requirements.

Literature Review

  • Reviewed the classifications of hip fractures, focusing on common fracture locations such as the femoral neck and intertrochanteric region.
  • Analyzed various risk factors associated with rapid bone density decrease leading to Osteoporosis. Factors include long-term use of high-dose steroid tablets, a family history of Osteoporosis, low body-mass index due to malnutrition or chronic illnesses, deficiency of essential nutrients like calcium and vitamin D, smoking and excessive alcohol consumption, and an inactive lifestyle. Also discussed hormone-related disorders as potential causes.
  • Investigated the demographics most affected by Osteoporosis, noting the high risk in women, especially those of Asian or Caucasian descent, who have undergone surgical removal of ovaries or experienced early menopause.
  • Explored the methods of diagnosing Osteoporosis, with a primary focus on the use of dual-energy X-ray (DXA) scans to measure bone density. Discussed the interpretation of results using ‘T Score’ and ‘Z Score’.
  • Examined the significant decrease in bone strength due to Osteoporosis, highlighting the correlation between reduced bone density and increased fracture risk. Reviewed the contrast between a healthy vertebra and an osteoporotic vertebra in terms of bone density.
  • Examined the impact of loading conditions on bone deformation and stress, highlighting the necessity of numerical methods, such as finite element analysis, for understanding this complex, heterogenous system.
  • Described the process of finite element analysis, emphasizing its basis in the generalized Galerkin methods, which convert continuous operator problems into discrete models.
  • Analyzed and compared three finite element models (Linear, Multi Point Constraints, Contact) used for simulating side-fall loading conditions on the femur. Highlighted the choice of the Linear model in this study due to its widespread use and lower potential for discrepancies.
  • Detailed the significant impact of side-fall orientation on maximum femoral strength (MFS), noting that a posterolateral fall is considered the weakest orientation.
  • Reviewed the range of loading orientations applied in the Ansys model, focusing on the importance of the neutral axis for verifying MFS in this study.
  • Overview of the CT2S (computer tomography to strength) pipeline, which uses non-invasive medical imaging to estimate bone strength, applying varying density material properties to a 3D geometry obtained from a CT scan.
  • Quantitative Computed Tomography (QCT)
    • Explanation of how CT scans are used to produce detailed cross-sectional images of the body, providing more detailed 3D images compared to 2D images from DXA scans.
  • Densitometric Calibration
    • Description of inline and external densitometric calibration methods to translate the grayscale image from the CT scan into accurate density values.
  • Segmentation
    • Introduction to image segmentation, the process of partitioning CT-scanned images into regions of interest corresponding to anatomical structures.
  • Meshing
    • Outline of the meshing process, which involves defining continuous geometric grids on the 3D geometry. Discussed the common mesh types: tetrahedral and hexahedral, and the necessity of conducting a mesh convergence/refinement study.
  • Material Properties
    • Explanation of how the heterogeneous material properties of bone are estimated through a series of calculations involving CT density, ash density, apparent density, and Young’s Modulus.
  • Reference System & Boundary Conditions
    • Description of how a reference system is selected based on anatomical landmarks identified in the CT images, and how boundary conditions are assigned based on the type of loading simulation required.
  • Simulation Results & Predicted Bone Strength
    • Overview of how simulation results are interpreted to identify potential femur fractures and predict maximum bone strength through the linear relationship between force and strain.
  • Introduction of Visualization Toolkit (VTK) file format used for visualizing scientific data and the INP format used for finite element analysis (FEA) in the Abaqus software.
  • Overview of VTK File Format
    • Explanation of the VTK file format used by the Visualization Toolkit for storing and exchanging visualization data.
    • Explanation of the VTK file’s role in geometry and mesh generation, material property assignment, and issue identification in model visualization.
    • Introduction of key components of VTK format files, including the header, file type, dataset structure, field data, points, cells, cell type, and associated cell and point data.
    • Explanation of the two file formats available in VTK, the legacy VTK (.vtk) file format, and the implications of using ASCII or binary formats.
  • Conversion of Binary Legacy VTK to ASCII Legacy VTK
    • Discussion on the process and benefits of converting binary legacy VTK files to ASCII legacy VTK files using ParaView.
    • Mention of the slight decrease in precision due to the conversion from binary to ASCII legacy VTK formats.
  • Overview of INP File Format
    • Explanation of the INP file format used by the Abaqus finite element analysis software.
    • Discussion of key sections in an Abaqus INP file, including the header, nodes, elements, materials, boundary and loading conditions, and output requests.
    • Highlight of the INP file’s readability and compatibility with Abaqus or other FEA software, preprocessors, or Python scripts.
  • General Legacy ASCII VTK to Abaqus INP File Conversion Process
    • Description of the process for converting Legacy ASCII VTK files to Abaqus INP files to perform simulations using the FEA software.
    • Discussion on the flowchart for the conversion process, including the analysis of ASCII legacy VTK file, extraction of required data, data preparation for conversion, and saving as an INP file.
    • Explanation of the sanity check or verification mode to ensure the accuracy of the conversion.

Methodology

Identifying Finite Element Model Loading & Boundary Conditions:

  • Explanation of the choice for the linear model from the most common boundary and loading condition models based on literature review, for application to the femur in side-fall configuration.
  • Description of the implementation of the linear model, involving constraints on specific nodes of the femur.

Modelling Various Loading Orientations:

  • Discussion on the simulation of 28 loading orientations based on the literature review, to ensure the accuracy of the finite element modelling approach, methods, and results post-processing.
  • Introduction of the anatomical coordinate system used for the hip, with the origin at the center of the femoral head.
  • Explanation of the determination of unit vectors 𝐹𝑥, 𝐹𝑦, and 𝐹𝑧 from the anatomical coordinate system and their use in expressing forces along various orientations.
  • Discussion on the normalization and application of force magnitude to the 𝐹(𝑗, 𝑘) values for each orientation.
  • Mention of the creation of a consistent set of load cases with the same force magnitude applied in different directions.
  • Explanation of the neutral loading direction, represented by 𝐹(3, 0), equalling the 𝐹𝑦 unit vector after normalization.
  • Accurately converted surface mesh, volume mesh, and volume mesh with material properties from VTK files to corresponding Abaqus INP files, employing Python-based approaches (the widely used Meshio library and a custom Python code).
  • Conducted thorough sanity checks of the converted INP files, which involved:
    • Comparing INP files produced using both the Meshio library and the custom codes to verify node and element connectivity data.
    • Performing a quantitative comparison of node and element connectivity data between the original VTK files and the converted INP files.
    • Conducting a visual inspection of the imported INP files in Abaqus to confirm the accurate representation of the original VTK files.
    • Verifying the material properties (Young’s Modulus, Poisson’s ratio, and density) in the converted INP files against the original VTK files.
    • Cross-checking material properties assignment between the converted INP file and Abaqus.
  • Identified the center node on the femoral head for load application using a node filtering and elimination technique based on x, y, and z coordinates. Performed a sanity check through Euclidean distance calculations using filtered and non-filtered node datasets to confirm accurate center femoral head node identification.
  • Used Python scripts and the linear modelling approach to identify nodes for constraint application on the greater trochanter and the femoral shaft.
  • Developed the “VTK to Abaqus INP Converter and Reader” application using Python and relevant libraries to automate and streamline:
    • Identification of VTK file and cell (Element) types.
    • Conversion of VTK files of different Element Types (Surface Mesh – CPE3, cell type 4; Volume Mesh – C3D10, cell type 24; Volume Mesh with Material Properties – C3D10, cell type 24) to corresponding Abaqus INP files.
    • Node identification for boundary and loading conditions, which includes the determination of the closest node to a reference point (based on a reference coordinate, not a node, for greater flexibility), identification of the least and greatest nodes along a coordinate axis (x, y, or z), and filtering of nodes based on the coordinate axis, direction (greatest or least), and a distance threshold from the extreme node on the selected coordinate system.
  • Adapted the INP files to apply specific boundary and loading conditions on desired nodes, as Abaqus doesn’t support node selection by node number and only offers manual node selection.
  • Loading Assignment on Femoral Head Node:
    • Defined the determined loading node (at the femoral head center) in the “Assembly” section of the INP file using a specific heading.
    • Applied the load to each INP file (for the 28 loading orientations) using components obtained from multiplying the normalized 𝐹(𝑗, 𝑘) component (for each specific loading case) by 500.
  • Constraining Node on Greater Trochanter:
    • Defined the determined most lateral greater trochanter node in the “Assembly” section of the INP file using a specific heading.
    • Applied the XSYMM boundary condition (symmetric boundary condition constraining the node’s movement in the X direction while allowing free translation in the Y and Z directions) to the node under the “BOUNDARY CONDITIONS” heading and within the initial step.
  • Constraint Assignment on Femoral Shaft Nodes:
    • Defined the femoral shaft nodes (distal end) to be constrained in all directions in the “Assembly” section of the INP file using a specific heading. All nodes (around 850) were defined using a list format to prevent Abaqus INP file processing errors.
    • Applied the “PINNED” boundary condition (constrains node set in all three translational degrees of freedom, U1 = X, U2 = Y, and U3 = Z, but is still free to rotate) to the node set under the “BOUNDARY CONDITIONS” heading and within the initial step.
  • Verification/Sanity Check Using Abaqus CAE:
    • Loaded the modified INP file, with the appropriate loading and boundary conditions included, into Abaqus to confirm whether the correct nodes were being loaded/constrained appropriately. The applied loading and boundary conditions perfectly resembled the Linear model reported in the literature and outlined in the modelling approach.

Fracture Sites:

  • Anticipated a femoral neck fracture as per the Ansys model.
  • Selected the neck surface area of the sub-capital and transcervical fractures’ locations as the region of interest (ROI) for the femur strength analysis.
  • Investigated only compressive (minimum principal) strains as the fracture was identified due to compression by the existing Ansys-based femur strength model.

Failure Criteria:

  • Implemented the maximum principal strain criterion to predict the femur strength.
  • Only considered minimum (compressive) strains as it was established that compressive strains caused the fractures. The maximum or tensile strains were not taken into account.
  • Averaged the compressive strains on a 3mm radius at the surface nodes to adhere to the continuum hypothesis and avoid the local effects of nodes.
  • Adopted the threshold compressive strain value of 1.04% as used in the study by Harun H. Bayratkar and in the existing Ansys-based femur strength model.

Results and Analysis

  • Implemented a post-processing approach to calculate minimum femur strength instead of defining the failure criterion in the simulation input file.
  • This approach provided more flexibility in results analysis, eliminated local effects, and was more resource-efficient.
  • Applied the maximum principal strain criterion to calculate minimum femur strength, using the linear relationship between the applied force and produced strain.
  • Minimum principal strains were averaged over a 3mm radius at surface nodes.

Results Post Processor Application

  • The post-processor tool was designed to analyse results, particularly nodes, coordinates, minimum principal strains, and MFS for each loading case.
  • The application integrated two key Python scripts to calculate the average minimum principal strain for each node and determine the average minimum principal strain and MFS for each node as the center node of the 3mm radius.
  • The tool enabled efficient and accurate processing of all results simultaneously through a unified process.
  • Two output files generated: One with the average minimum principal strain and MFS values and the other with details of nodes within a 3mm radius when each node is considered the center node.
  • Nodes within the 3mm radius of each node are identified, with checks performed using Abaqus Visualization to ensure accuracy and consistency.
  • MFS values sorted in ascending or descending order to identify the center node with the desired MFS.

Comparing Abaqus-Predicted MFS to Femur Strengths from the Ansys Model

  • Compared the failure load predicted by the Ansys model to the minimum femur strength predicted by the current model.
  • The predicted minimum femur strength for different loading orientations was obtained, with the MFS for the neutral loading orientation being 2162.7N.
  • The average minimum femur strength predicted by the current model was 2394N (±289) for 15 loading orientations and 2233N (±277) for all 28 loading orientations.
  • Results showed that the Abaqus model developed in this study performs closely with the Ansys-based model, with some discrepancies for posterior loading orientations with medial angles greater than 0°.

Comparing Abaqus-Predicted MFS to Femur Strengths from Clinical Trials

  • The femur strengths of different patients from clinical trials were compared to the femur verified in Abaqus.
  • A normalized histogram showed that the determined MFS from the Abaqus model was specific to the patient and not general to all patients, thus confirming the accuracy of the femur strength analysis pipeline.
  • Slight discrepancy noted between the Abaqus model’s MFS and patient 12’s femur strengths due to error from the inaccurate application of the boundary condition at the greater trochanter.

Impact of the Constrained Distal Femoral Shaft Nodes (𝑽 Value) on the MFS – Sensitivity Analysis

  • The impact of the number of nodes pinned at the distal end of the femoral shaft on the minimum femur strength was analyzed.
  • The variation in minimum femur strengths for different V values reached a maximum of around one-tenth of a percent from the average MFS.
  • Though the impact on the minimum femur strength was not substantial, a 𝑉 value of 0.75 was selected for this study because of it being the average value with minimal deviation and most reasonable for finite element modelling.

Job Summary:

  • The model consists of 198,489 elements and 872,589 nodes, with 277,122 nodes defined by the user and 595,467 generated internally by Abaqus.
  • Each node was allowed 3 degrees of freedom (DOF), yielding a total of 2,220,789 variables in the model.
  • The minimum memory required to run the analysis was approximately 2.49 GB, while the memory to minimize input/output operations was 20.12 GB.
  • The sparse, direct linear equation solver was used due to the sparse nature of the stiffness matrices in the model.

Computational Resource:

  • Allocating more cores and memory significantly reduced computational time.
  • Utilizing Hyper-Threading significantly reduced (by a factor of 10 times) the computation time.
  • The use of High-Performance Computing (HPC) at The University of Sheffield (specifically the Bessemer), further reduced the wallclock time to just 3 minutes, with eight cores and 32 GB memory allocated.
  • The time for computation did not vary significantly across different loading orientations.

Limitations and Future Work

  • The primary limitation was the inaccurate application of the boundary condition at the greater trochanter, confining the verification to 15 out of 28 simulated loading orientations.
  • Soft tissues surrounding the femur, particularly the greater trochanter, were not considered, potentially influencing the predicted minimum femur strength due to damping effects.
  • The impact of clothing on femur strength was not taken into account, which may contribute to more conservative minimum femur strength predictions.
  • The study did not consider multiple contact points during a fall, suggesting a dynamic model that considers this factor would be more appropriate.
  • In future work, the node to be constrained on the greater trochanter should be the one farthest along the loading direction, determined by taking the dot product of the coordinate vector of all the nodes on the surface of the greater trochanter with the 𝐹(𝑗, 𝑘) vector.
  • Implementing multi-point constraint (MPC) and contact models could yield results more aligned with clinical trials.
  • An advanced dynamic model considering patient’s weight, height, soft tissues/clothing-imposed damping effects, and multiple contacts upon impact would offer a more accurate assessment of hip fracture risk.
  • The current “VTK to Abaqus INP Converter & Reader” application could be extended to allow conversions of more than two element types and to include all necessary conversions and assignments before importing to Abaqus.
  • The “Results Post Processor” application could be improved for further automation and direct results processing.
  • These improvements might increase the adoption of Abaqus in the biomedical FEA industry, bridge the gap in user preference between Ansys and Abaqus, and offer a more comprehensive and seamless simulation process.
  • The femur strength pipeline could provide valuable insights for targeted rehabilitation strategies, ergonomics/biomechanical studies, bone health research, and orthopedic implant design and development.
  • The femur strength pipeline and the two applications could also be adapted for other fracture-prone sites such as the spine or pelvis, relevant for geriatrics with osteoporosis.