Computation for biomechanical analysis of aortic aneurysms: the importance of computational grid

Abstract

Aortic wall stress is the most common variable of interest in abdominal aortic aneurysm (AAA) rupture risk assessment. Computation of such stress has been dominated by finite element analysis. However, the effects of finite element (FE) formulation, element quality, and methods of FE mesh construction on the efficiency, robustness, and accuracy of such computation have attracted little attention. In this study, we fill this knowledge gap by comparing the results of the calculated aortic wall stress for ten AAA patients using tetrahedral and hexahedral meshes when varying the FE formulation (displacement-based and hybrid), FE shape functions, spatial integration scheme, and number of elements through the wall thickness.

Keywords:

• Abdominal aortic aneurysm
• computational biomechanics
• stress analysis
• finite element analysis

1. Introduction

An abdominal aortic aneurysm (AAA) is an abnormal dilatation of the lower part of the aorta by 1.5 times its normal diameter (Johnston et al. 1991). The most severe complication of AAA, aortic rupture, is a catastrophic clinical event with a mortality rate of up to 90% (Upchurch and Schaub 2006). Currently, AAA clinical management is based on the maximum aortic diameter. However, 60% of AAAs with aortic diameter over 55 mm remain stable and do not rupture in the patient’s lifetime (Darling et al. 1977). Multiple studies have proposed that patient-specific biomechanical analysis can be potentially a more accurate and reliable solution for predicting AAA rupture than the AAA diameter measurement (Gasser et al. 2010; Polzer et al. 2013; Gasser et al. 2014; Joldes et al. 2017; Polzer et al. 2020).

Biomechanical analysis of AAA strongly relies on the computation of stress within the aortic wall (Fillinger et al. 2002). Aortic wall stress has been the main variable of interest in almost all proposed biomechanical rupture risk indicators of AAA (Fillinger et al. 2003; Gasser et al. 2010; Gasser et al. 2014). However, some studies found no evidence of a relationship between AAA symptoms and stress in the AAA wall (Miller et al. 2020). Evaluation of the extent to which AAA stress can be regarded as a reliable and robust indicator AAA rupture risk is beyond the scope of this study.

Finite Element Analysis (FEA) has been the method of choice for computing stress in AAA walls (Raghavan et al. 2006; Joldes et al. 2016; Wittek et al. 2022). This includes commercial and open-source software packages for AAA biomechanical analysis, such as A4Clinics VASCOPS (http://www.vascops.com/) and BioPARR (https://bioparr.mech.uwa.edu.au/). Progress in the application of FEA in computing stress in the AAA wall has been underpinned by a substantial research effort in methods for incorporating artery residual stress (Polzer et al. 2013; Joldes et al. 2018), determining and formulating material models for aortic tissue (Holzapfel et al. 2000), determining the relationship between the aortic wall thickness and predicted stress (Conlisk et al. 2016; Miller et al. 2019), and investigating the effects of intraluminal thrombus (ILT) on the AAA wall stress (Joldes et al. 2017). However, the effects of finite element (FE) formulation (displacement-based and hybrid), element topology (tetrahedra and hexahedra), order of polynomial shape functions, spatial integration scheme, and the FE mesh quality on the results of AAA stress computation have attracted little attention of the AAA biomechanics community. This is despite the fact that the accuracy of stress predicted using FEA and the computation time needed to obtain the results depend on the element topology, formulation, and quality of FE mesh used (Cifuentes and Kalbag 1992; Tadepalli et al. 2011). In particular for the peak stress, a variable often used in biomechanical indicators of the AAA rupture risk (Fillinger et al. 2002; Singh et al. 2020).

Discussion of whether tetrahedral or hexahedral elements are most suitable for continua with complex geometry has continued for more than 30 years (Cifuentes and Kalbag 1992; Benzley et al. 1995; Knupp 2001; Wittek et al. 2016). Tetrahedral elements tend to be used in AAA computational models as they allow automated mesh generation (Raut et al. 2015; Joldes et al. 2017). However, a high number of tetrahedral elements may be required to obtain a converged solution, with AAA computational biomechanics models typically consisting of 500,000 to over a million tetrahedra (Miller et al. 2020; Wittek et al. 2022). For the same element edge size and the same number of nodes, tetrahedral meshes require around seven times more elements than hexahedral meshes (Cubit 2021). Furthermore, for hexahedral meshes, the solution converges faster compared to tetrahedral meshes (Alkhatib et al. 2023). Meshes consisting of 30,000 to 50,000 hexahedral elements have been used in AAA computational biomechanical analysis (Wittek et al. 2022).

In this study, we compare the results of AAA wall stress computation for ten AAA patients using FEA when varying the element topology (i.e. tetrahedra and hexahedra), FE shape functions, element quality, spatial integration scheme, and number of elements across the AAA wall thickness.

2. Materials and methods

2.1. Computational biomechanical analysis of AAA

We used the method embedded in the freely-available open-source software platform BioPARR – Biomechanics-based Prediction of Aneurysm Rupture Risk (https://bioparr.mech.uwa.edu.au/) (Joldes et al. 2016; Joldes et al. 2017) to perform the biomechanical analysis of AAA in this study. BioPARR relies on the observation that as patient images show the deformed AAA geometry (due to loading by the blood pressure), the internal forces within the blood vessel wall that balance the applied pressure can be calculated from the principles of statics (Joldes et al. 2016; Joldes et al. 2017). This implies that if the computational biomechanics simulations are set-up in such a way that the deformed AAA geometry remains unchanged when loaded by blood pressure, one should obtain the stress field in the AAA wall that balances the pressure (Joldes et al. 2017). Such stress field is for practical purposes independent of the material properties of aortic tissue (Zelaya et al. 2014; Biehler et al. 2015).

The residual stresses are incorporated according to Fung’s Uniform Stress Hypothesis (USH) (Fung 1991) by averaging the stress across the vessel wall thickness (Joldes et al. 2018). This makes the computation of residual stress in the AAA wall a simple post-processing step during AAA biomechanical stress analysis (Joldes et al. 2016). We previously verified (Joldes et al. 2018) that for a given material model and material properties, the residual stress fields obtained using this method are very close to the ones obtained by the more complex non-linear method by Polzer et al. (Polzer et al. 2013).

For AAA stress computation, BioPARR can use any verified linear FE solver. Following our previous studies of AAA biomechanics (Joldes et al. 2016; Joldes et al. 2017), we used ABAQUS/Standard FE solver (Abaqus 2018) (https://www.3ds.com/products-services/simulia/products/abaqus/) called from within BioPARR.

2.1.1. Boundary conditions and load

The top and bottom nodes of the aneurysm were rigidly constrained. As the patients were asked to hold their breath during the image acquisition, there was no need to account for the respiration-induced aorta movement. The spine and tissues surrounding the aorta were not included in the AAA computational biomechanics models used in this study. Such an approach is frequently used in the literature (Gasser et al. 2010; Joldes et al. 2016), and its validity has been confirmed by our recent study (Liddelow et al. 2023).

The mean arterial blood pressure (Cywinski 1980), measured five minutes prior to the image acquisition, was used as patient-specific load. Following (Vorp et al. 1996), the pressure load was uniformly distributed over the ILT luminal surface.

2.1.2. Material models and material properties

We used a linear elastic material model for the aortic wall and ILT tissues, which is an approach used in BioPARR. Validity of application of linear material model in computing stress in AAA wall has been confirmed in our previous studies (Joldes et al. 2016; Joldes et al. 2017; Wittek et al. 2022).

The aorta was defined as a very stiff material with Young’s modulus of 100 MPa and Poisson’s ratio of 0.49 (Joldes et al. 2016; Joldes et al. 2017; Miller et al. 2020). This ensures that the deformed geometry of AAA represented by the model remains unchanged when loaded by blood pressure. The ILT was assumed to be 20 times more compliant than the AAA wall with Young’s modulus of 5 MPa and Poisson’s ratio of 0.49 (Miller et al. 2020).

2.2. Patient-specific AAA geometry from medical images

We used ten anonymised contrast-enhanced computed tomography angiography (CTA) image data sets of AAA patients. Patients were recruited at Fiona Stanley Hospital (Western Australia, Australia) and gave their informed consent before their participation in the study. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by Human Research Ethics and Governance at South Metropolitan Health Service (HREC-SMHS) (approval code RGS3501), and by Human Research Ethics Office at the University of Western Australia (approval code RA/4/20/5913).

We extracted patient-specific AAA geometries through CTA image segmentation using 3D Slicer (https://www.slicer.org/), an open-source software platform for medical image informatics, image processing, and three-dimensional visualization (Fedorov et al. 2012). Region of interest for the AAAs was defined as distal to the renal arteries and proximal to the common iliac bifurcations (Figure 1(a)).

Figure 1. Example of patient-specific abdominal aortic aneurysm (AAA) geometry (Patient 1): (a) AAA model segmented from CTA image – sagittal view, and (b) a cross-section of AAA geometry showing the intraluminal thrombus (ILT), lumen, and AAA external wall (blue color) created automatically using BioPARR (Joldes et al. 2017).

As contrast enhancement was used for all patients, we were able to segment the lumen automatically using the threshold algorithm (Lesage et al. 2009) in 3D Slicer. We manually segmented the external surface of the AAA wall as the AAA wall intensity in the CTA images is notably similar to that of the surrounding tissues.

Surfaces defining the AAA wall and the ILT were automatically created from the AAA segmentations using BioPARR software (Figure 1(b)). We assumed a constant AAA wall thickness of 1.5 mm (Miller et al. 2020; Polzer et al. 2020).

2.3. Patient-specific finite element meshes for biomechanical analysis of AAA

2.3.1. Generation of tetrahedral and hexahedral finite element meshes of AAA

Patient-specific tetrahedral FE meshes were created automatically from the extracted AAA geometries using the open-source Gmsh mesh generator (https://gmsh.info/) (Geuzaine and Remacle 2009), called from within BioPARR software. To generate hexahedral meshes of the aneurysm wall, we followed the semi-automated procedure described in detail in our previous work (Alkhatib et al. 2023), using the HyperMesh FE mesh generator (Altair, MI, USA; https://www.altair.com/hypermesh).

2.3.2. Finite element mesh quality measures

For tetrahedral meshes, we used two element quality measures (Parthasarathy et al. 1994):

1. The minimum/maximum allowable interior angles, which are determined between any two edges of the tetrahedron face. The acceptable interior angle ranges between 30° and 120° (Yang 2018).

2. The volumetric skew (AltairEngineering 2022) that determines how close a given element is to the ideal reference shape (in the sense of isoparametric FE analysis) (Ruggiero et al. 2019).

(1) $$\text{Volumetric skew}=1-\frac{\text{Actual tetrahedron volume}}{\text{Ideal tetrahedron volume}},$$(1) where the ideal tetrahedron is an equilateral tetrahedron with the same circumradius as the actual (analysed) tetrahedron. The circumradius is the radius of a sphere passing through all four vertices of the tetrahedron.

The volumetric skew is zero for an ideal tetrahedron, and one for a degenerated (co-planar/‘flat’) tetrahedron. Values recommended in practice vary between zero and 0.98 (recommendation from HyperMesh (AltairEngineering 2022) and ABAQUS (Abaqus 2018)). In this study, we regarded elements with volumetric skewness > 0.98 exhibiting unacceptable poor-quality (slivers).

For hexahedral meshes, we used two element quality measures (Zhang 2018):

1. Jacobian ratio (normalised Jacobian), defined as the ratio between the smallest and largest Jacobian matrix determinant for each Gauss point in all elements (Yang 2018). It ranges from zero to one, where one indicates an ‘ideal’ hexahedral element. We follow the recommendations by Yang (Yang 2018) who suggested that a Jacobian ratio should be at least 0.6.

2. The minimum and maximum allowable interior angles of a quadrilateral element face. The acceptable range of interior angles is between 45° and 135° (Yang 2018).

2.3.3. Finite element formulation and spatial integration

FEA of incompressible and nearly incompressible continua, such as aortic wall tissue, requires the application of finite elements that are not susceptible to volumetric locking (artificial stiffening for Poisson’s ratio close to 0.5). For linear 4-noded tetrahedral elements, a range of volumetric locking countermeasures, such as mixed (hybrid) FE formulation, need to be applied (Bonet and Burton 1998; Bathe 2006; Joldes et al. 2009). Quadratic 10-noded tetrahedral and 20-noded hexahedral elements are typically regarded as not susceptible to volumetric locking. However, accurate computation of the stress field for the Poisson’s ratio close to 0.5 may still require application of the hybrid FE formulation (Bathe 2006).

Reduced spatial integration has been reported in the literature to improve the accuracy of FE approximations (Zienkiewicz et al. 1971; Sandhu and Singh 1978; Bathe 2006) and reduce the computation cost by 30% compared to fully integrated quadratic hexahedral elements (Abaqus 2018). It has also been indicated that quadratic elements using reduced spatial integration provide accurate results for nearly incompressible materials (Naylor 1974; Malkus and Hughes 1978). This makes such elements particularly suitable for modelling the aorta and other soft tissues. Although quadratic tetrahedral elements with reduced integration have been described in the literature (Malkus and Hughes 1978), they are rarely used in practice.

We used 10-noded quadratic tetrahedral, and 8-noded linear and 20-noded quadratic hexahedral elements. We varied the FE formulation for each element between displacement-based and hybrid formulations. We also varied the spatial integration using full and reduced integration (Gauss) points, Table 1 summarises the finite elements analysed in this study.

Table 1. Tetrahedral and hexahedral finite elements from ABAQUS finite element code library (Abaqus 2018) analysed in this study.

Element name in ABAQUS	Element description	No. of nodes	No. of integration (Gauss) points	Finite element formulation
C3D10	10-node quadratic tetrahedron	10	4	Displacement-based
C3D10H	10-node quadratic tetrahedron, hybrid with constant pressure	10	4	Hybrid (Displacement/pressure-based)
C3D20	20-node quadratic hexahedron	20	27	Displacement-based
C3D20H	20-node quadratic hexahedron, hybrid with linear pressure	20	27	Hybrid (Displacement/pressure-based)
C3D20R	20-node quadratic hexahedron, reduced integration	20	8	Displacement-based
C3D20RH	20-node quadratic hexahedron, reduced integration, hybrid with linear pressure	20	8	Hybrid (Displacement/pressure-based)
C3D8	8-node linear hexahedron	8	8	Displacement-based
C3D8H	8-node linear hexahedron, hybrid with linear pressure	8	8	Hybrid (Displacement/pressure-based)
C3D8R	8-node linear hexahedron, reduced integration with hourglass control	8	1	Displacement-based
C3D8RH	8-node linear hexahedron, reduced integration with hourglass control, hybrid with constant pressure	8	1	Hybrid (Displacement/pressure-based)

2.3.4. AAA mesh convergence analysis

In this study, we focused on analysing the effects of mesh refinement for linear and quadratic hexahedral elements (Babuska and Szabo 1982) on the results of peak maximum principal stress (MPS). For linear hexahedral elements, we used three meshes with two, four, and eight hexahedral elements through the AAA wall thickness. For quadratic hexahedral elements, we used three meshes with one, two, and four hexahedral elements through the wall thickness.

We used the Grid Convergence Index (GCI) (Roache 1994) to quantify how far the obtained solution is from the estimated ‘asymptotic’ value. We selected GCI as it has been used successfully in solid and structural mechanics (Schwer 2008; Kwaśniewski 2013). GCI is calculated based on Richardson’s extrapolation (Slater 2021): (2) $$\mathrm{\text{GCI}}=F_{\mathrm{s}}\frac{\epsilon }{r^p-1}\times 100\mathrm{\%},$$(2) where F_s is a safety factor based on the experience from computational fluid dynamics calculations (F_s equals one in our study), ε is the relative difference between subsequent solutions, r is a mesh refinement ratio (r equals two in our study), and p is the order of convergence given by: (3) $$p=\frac{\mathrm{\text{ln}}\hspace{0.17em}\left(\frac{f_3-f_2}{f_2-f_1}\right)}{\mathrm{\text{ln}}\left(r\right)},$$(3) where f₁, f₂, and f₃ are the results of the variable of interest (peak MPS in this study) for three successive meshes refined by the factor r.

In-depth mesh convergence analysis was conducted for Patient 1 only. As residual stress affects the stress distribution (makes it uniform) through the AAA wall, we did not include it in the analysis of mesh convergence and effects of FE formulation on the computed AAA stress.

2.4. Analysis of the results and simulation matrix

We analysed the MPS at the nodes of the studied computational biomechanical FE models of AAA (Table A1). Stress at each node was computed by averaging the stress tensor components from all the elements with a vertex located at a given node (Abaqus 2018). When comparing the results, we analysed the MPS magnitude and distribution throughout the AAA and constructed MPS histograms and percentile plots with a percentile rank (in the population of all nodes in the AAA FE model of a given patient) of the MPS at a given node on the horizontal axis and MPS on the vertical axis (Wittek et al. 2022).

3. Results

3.1. Generation of tetrahedral and hexahedral finite element meshes of AAA

To compare the time required to generate patient-specific tetrahedral and hexahedral meshes of AAA and analyse the quality of elements in tetrahedral and hexahedral meshes, we built meshes with two elements through the AAA wall thickness (Figure 2). The automated tetrahedral mesh generation (with optimisation to ensure the element quality) took between 40–50 min per patient using Intel(R) Core (TM) i7-5930K CPU @3.50 GHz with 64GB of RAM running Windows 8 operating system. The number of quadratic tetrahedral elements in the AAA wall mesh varied between around 620,000–860,000 (Table A2), and the number of nodes between 820,000 and 1,300,000.

Figure 2. Example of finite element meshes generated in this study (Patient 3). Blue color: aneurysm wall (with a constant wall thickness of 1.5 mm). Red color: intraluminal thrombus (ILT) meshed using tetrahedral elements. (a) Tetrahedral meshes of aneurysm wall, and (b) hexahedral meshes of aneurysm wall. Both meshes have two layers of elements through the wall thickness.

AAA hexahedral mesh generation required 20–30 min of work by the analyst (F.A., a PhD student with over five years of experience in FEA) for each patient. Hexahedral meshes of the AAA walls created in this study consisted of 21,000 to 53,000 hexahedral elements (Table A2) and between 116,000 and 290,000 nodes.

For each FE mesh, the number of poor-quality elements was less than 0.5% of the total elements (Table A2). Most of the poor-quality elements were located at the top and bottom edges of the AAAs. These areas are of no interest from the perspective of stress analysis in the AAA wall.

3.2. AAA mesh convergence analysis and effects of element formulation

3.2.1. Quadratic tetrahedral elements

We observed differences in peak values of MPS between the displacement-based (C3D10) and hybrid (C3D10H) FE formulations (Table A3 and Figure 3). We could not calculate the convergence parameters for the element type C3D10 because of the inconsistency of the MPS results and observed clear symptoms of volumetric locking (Figure A1). This indicates that the hybrid formulation is essential when using the quadratic tetrahedral elements for AAA stress calculations. From very small (practically negligible) differences between the MPS histograms for two and four quadratic tetrahedral elements through the AAA wall thickness (Table A3), we conclude that two layers of quadratic tetrahedral elements through the wall thickness are sufficient to achieve converged results.

3.2.2. Linear hexahedral elements

The differences between the peak MPS obtained using displacement-based (C3D8 and C3D8R) and hybrid (C3D8H and C3D8RH) FE formulations were negligible. This observation applies to full (C3D8 and C3D8H) and reduced (C3D8R and C3D8RH) spatial integration elements. The results show only a moderate (order of 25%) increase in the computation time due to the application of hybrid formulation and a small (order of -15%) decrease in the computation time for reduced spatial integration (Table A3).

Convergence analysis indicates virtually no difference in the grid convergence parameters and the estimated asymptotic solution (f_estimated) of peak MPS between linear hexahedral elements using displacement-based and hybrid formulations. This observation is consistent with the results reported when comparing the peak MPS for displacement-based and hybrid FE formulations (Table A3).

Elements with reduced spatial integration (C3D8R and C3D8RH) exhibited a lower convergence order and a higher GCI than elements with full integration (C3D8 and C3D8H).

Noticeable differences between the peak MPS obtained using the meshes constructed with two elements through the AAA wall and the f_estimated of peak MPS of 0.3192 MPa (Table A3), as well as the differences in the MPS magnitude obtained using the meshes with two and four hexahedral elements through the AAA wall (Figure 3) suggest that four layers of linear hexahedral elements through the AAA wall may be needed to achieve converged results.

Figure 3. Histograms of the maximum principal stress in the aneurysm wall (Patient 1) calculated using tetrahedral and hexahedral finite element meshes with different number of elements through the wall thickness (1 layer, 2 layers, 4 layers and 8 layers of elements) and for different element type/formulation (a) Quadratic tetrahedra with displacement formulation (C3D10); (b) Quadratic tetrahedra with hybrid formulation (C3D10H); (c) Linear hexahedra with displacement (C3D8) and hybrid (C3D8H) formulations. Only three histograms are visible as the results for C3D8 and C3D8H overlap; (d) Linear hexahedra with reduced integration, and displacement (C3D8R) and hybrid (C3D8RH) formulations. Only three histograms are visible as the results for C3D8R and C3D8RH overlap; (e) Quadratic hexahedra with displacement formulation (C3D20); (f) Quadratic hexahedra with hybrid formulation (C3D20H), (g) Quadratic hexahedra with reduced integration and displacement formulation (C3D20R), and (h) Quadratic hexahedra with reduced integration and hybrid formulation (C3D20RH).

3.2.3. Quadratic hexahedral elements

Noticeable differences between the peak MPS obtained using fully integrated displacement-based (C3D20 and C3D20R) and hybrid (C3D20H and C3D20RH) FE formulations were observed. The differences for reduced spatial integration were much smaller and can be regarded as negligible (Table A3 and Figure 3). As the constitutive behaviour of the AAA wall, tissues were described using linear elastic material with Poisson’s ratio of 0.49, the observed differences between the peak MPS for displacement-based and hybrid elements may indicate volumetric locking in very few (located at the top or bottom edges of the wall) quadratic hexahedral elements if full spatial integration and displacement-based formulation are used. This is further confirmed by the differences in the pressure (stress tensor trace) magnitude within the AAA wall obtained using such elements (peak pressure of 0.9178 kPa) and the hybrid elements (peak pressure of 0.6659 kPa).

The results for quadratic hexahedral elements using hybrid formulation indicate a higher convergence rate and lower CGI for reduced (C3D20RH) than full (C3D20H) integration. They also point to around a 10% decrease in the computation time for reduced integration (Table A3).

Negligible differences between the peak MPS obtained using the meshes constructed with two elements through the AAA wall and the f_estimated of the peak MPS of 0.3129 MPa (Table A3), as well as very small differences in the MPS magnitude obtained using the meshes with two and four quadratic hexahedral elements through the AAA wall, suggest that two layers of quadratic hexahedral elements through the AAA wall are sufficient to achieve converged results in terms of computation of the AAA wall stress. Detailed analysis indicates that for the discretisation using two layers of quadratic hybrid hexahedral elements through the AAA wall, the MPS magnitude and distribution obtained using full (C3D20H) and reduced (C3D20RH) spatial integration are practically indistinguishable (Figure 3).

3.3. Comparison of aneurysm wall stress computed using tetrahedral and hexahedral meshes

In this section, we compare the computation time, peak values, and distribution of the MPS in the AAA walls obtained using the 10-noded quadratic tetrahedral (C3D10H) and the 20-noded quadratic reduced integration hexahedral (C3D20RH) meshes with two elements through the AAA wall thickness. The comparison was conducted for all ten patients analysed in this study. The residual stress was included when computing stress in the AAA walls according to Fung’s Uniform Stress Hypothesis (Fung 1991).

Distributions of the MPS in the AAA wall predicted using the quadratic tetrahedral and hexahedral elements were very similar. Contour plots of the MPS for the tetrahedral meshes were practically indistinguishable from those for hexahedral meshes (Figure 4). This observation is consistent with the fact that no appreciable differences in the general behaviour of the percentile plots of the MPS (Figure 4) and peak MPS values (Table A4) could be observed between the tetrahedral and hexahedral meshes, although differences in the stress magnitude were noticeable for some patients (Patients 3, 4, 5, and 7). The relative absolute difference between the MPS in AAA wall obtained using tetrahedral and hexahedral meshes varied between less than 1% (Patients 3, 4, and 9) to up to 4% (Patient 10) (Table A4 in Appendix A).

Figure 4. (a) Contour plots of the maximum principal stress (MPS) including residual stress (RS) obtained using C3D10H 10-noded quadratic tetrahedral hybrid elements (left-hand-side figures) and C3D20RH 20-noded quadratic hexahedral hybrid elements with reduced integration (right-hand-side figures) for the aneurysm wall; and (b) percentile plots of MPS for the studied AAA patients.

To provide quantitative measure of the differences between the results obtained using tetrahedral and hexahedral meshes, we conducted statistical analysis of the peak values of MPS in the AAA walls discretised using tetrahedral and hexahedral elements. We used the non-parametric Mann-Whitney U test (null hypothesis that the difference between the median peak MPS obtained using tetrahedral and hexahedral meshes of AAA wall equals zero) (Mann and Whitney 1947), and two sample t-test (null hypothesis that the difference between the mean peak MPS obtained using tetrahedral and hexahedral meshes of AAA wall equals zero) (Rice 2006). The application of Mann-Whitney U test gave a pvalue of 0.5708 and t-test a pvalue of 0.8334. When treating patient-specific computation of the peak value of MPS in the aneurysm wall using AAA tetrahedral and hexahedral meshes as independent observations, these results should be interpreted as an indication that the differences between the peak MPS obtained using tetrahedral and hexahedral meshes may be noticeable but are not significant.

4. Discussion

Our results indicate that FE meshes constructed using quadratic rather than linear elements are preferred for computing stress in the AAA walls. This conclusion applies to both tetrahedral and hexahedral elements. For quadratic tetrahedral and quadratic hexahedral elements, the converged solution was obtained for two elements through the AAA wall thickness. For linear hexahedral elements, four elements through the AAA wall thickness were needed, i.e. quadratic hexahedral elements facilitated the converged solution for around two times smaller number of nodes than linear hexahedral elements. Furthermore, linear tetrahedral elements are prone to volumetric and shear locking (Bonet and Burton 1998; Bathe 2006; Joldes et al. 2009), a challenge that needs to be addressed when analysing incompressible (or nearly incompressible) soft continua such as aortic wall tissues. The results obtained in this study suggest that even for quadratic hexahedral elements with full spatial integration, hybrid FE formulation may be needed to ensure that volumetric locking does not affect the results of stress computation in the AAA wall.

The distribution and peak value of the MPS in the AAA wall obtained using converged meshes (two elements through the AAA wall thickness) constructed from quadratic hexahedral and quadratic tetrahedral elements with hybrid FE formulation were for practical purposes the same. Hence, the key factors to consider when deciding whether to use tetrahedral or hexahedral meshes are the time and effort required to build patient-specific meshes of AAA and the computation time associated with a given mesh type.

As AAA hexahedral meshes contain substantially fewer elements and nodes (and degrees of freedom) than tetrahedral meshes, their computation time was more than five times shorter than for the tetrahedral meshes. However, it should be noted that the application of more advanced computing hardware would reduce the computation time, making the computational efficiency of hexahedral meshes less important. The need for 20–30 min of experienced FE analyst time when generating AAA hexahedral meshes can be seen as important disadvantage in clinical environment (hospitals are unlikely to hire experts in computational mechanics).

5. Conclusion

The results obtained in this study indicate that quadratic finite elements with a hybrid (displacement/pressure) formulation are needed to ensure fast (for a smaller number of elements and nodes) convergence and to prevent volumetric locking when applying finite element models in patient-specific biomechanical analysis of Abdominal Aortic Aneurysms. In contrast, automation of discretization of 3D objects with complex geometry using hexahedra remains an unsolved problem that is a subject of ongoing research effort. Therefore, we conclude that tetrahedra rather than hexahedra are the element of choice for AAA patient-specific biomechanical models.

Acknowledgements

Authors thank Fiona Stanley Hospital (Murdoch, Western Australia) for providing patient images. Author F.A. thanks Ms. Giuliana D'Aulerio of the University of Western Australia Medical School/Division of Surgery and Ms. Kiara Beinart of the University of Western Australia Intelligent System for Medicine Lab for their support in preparing the ethics approval application. Contributions of Christopher Wood and Jane Polce, radiology technicians at Medical Imaging Department, Fiona Stanley Hospital (Murdoch, Western Australia) to patient image acquisition are gratefully acknowledged.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

Author F. A. acknowledges support of the Australian Government Research Training Program Scholarship at The University of Western Australia. This work was supported by the National Health and Medical Research Council NHMRC Ideas grant no. APP2001689.

Appendix A

Figure A1. AAA biomechanical analysis (the figure shows the pressure distribution in MPa) for Patient 1 using finite element mesh constructed with 10-noded quadratic displacement-based tetrahedral elements (element type C3D10 in ABAQUS finite element code). Volumetric locking, as indicated by nonphysical chess-board-like pattern of the pressure distribution, is clearly visible.

Table A1. Summary of the patient-specific computational biomechanics finite element (FE) models of AAA created and used in this study; a total of 60 models. The number of elements through the wall thickness determines the AAA mesh density. ILT is the intraluminal thrombus. The symbol ✗ indicates that ILT was not included in the model.

Patient no.	Applied MAP (kPa)	Model no.	Finite Element type (from ABAQUS)		No. of elements across AAA wall thickness
			AAA Wall	ILT
1	12	1	C3D10H	C3D10H	2
		2	C3D20H	C3D10H	2
		3	C3D20RH	C3D10H	2
		4	C3D10	✗	1
		5	C3D10H	✗	1
		6	C3D10	✗	2
		7	C3D10H	✗	2
		8	C3D10	✗	4
		9	C3D10H	✗	4
		10	C3D20	✗	1
		11	C3D20H	✗	1
		12	C3D20R	✗	1
		13	C3D20RH	✗	1
		14	C3D8	✗	2
		15	C3D8H	✗	2
		16	C3D8R	✗	2
		17	C3D8RH	✗	2
		18	C3D20	✗	2
		19	C3D20H	✗	2
		20	C3D20R	✗	2
		21	C3D20RH	✗	2
		22	C3D8	✗	4
		23	C3D8H	✗	4
		24	C3D8R	✗	4
		25	C3D8RH	✗	4
		26	C3D20	✗	4
		27	C3D20H	✗	4
		28	C3D20R	✗	4
		29	C3D20RH	✗	4
		30	C3D8	✗	8
		31	C3D8H	✗	8
		32	C3D8R	✗	8
		33	C3D8RH	✗	8
2	13	1	C3D10H	C3D10H	2
		2	C3D20H	C3D10H	2
		3	C3D20RH	C3D10H	2
3	13	1	C3D10H	C3D10H	2
		2	C3D20H	C3D10H	2
		3	C3D20RH	C3D10H	2
4	14	1	C3D10H	C3D10H	2
		2	C3D20H	C3D10H	2
		3	C3D20RH	C3D10H	2
5	14	1	C3D10H	C3D10H	2
		2	C3D20H	C3D10H	2
		3	C3D20RH	C3D10H	2
6	15	1	C3D10H	C3D10H	2
		2	C3D20H	C3D10H	2
		3	C3D20RH	C3D10H	2
7	9	1	C3D10H	C3D10H	2
		2	C3D20H	C3D10H	2
		3	C3D20RH	C3D10H	2
8	13	1	C3D10H	C3D10H	2
		2	C3D20H	C3D10H	2
		3	C3D20RH	C3D10H	2
9	12	1	C3D10H	C3D10H	2
		2	C3D20H	C3D10H	2
		3	C3D20RH	C3D10H	2
10	13	1	C3D10H	C3D10H	2
		2	C3D20H	C3D10H	2
		3	C3D20RH	C3D10H	2

Table A2. Tetrahedral (tetra) and hexahedral (hexa) AAA wall mesh with two layers of elements through the wall thickness generated in this study: number of elements, percentage of poor-quality elements with respect to the total number of elements for each part, and values of mesh quality measures.

Patient no.	Wall mesh topology	No. of elements	Poor-quality elements (%)	Min. interior element angle	Max. interior element angle	Max. volumetric skew	Min. Jacobin ratio
1	Tetra	718,738	0%	3.32	173.37	1.00 (1 element)	–
	Hexa	21,090	0%	49.73	131.16	–	0.86
2	Tetra	849,479	0%	1.47	177.05	1.00 (16 elements)	–
	Hexa	35,616	0%	45.31	140.64	–	0.73
3	Tetra	662,876	0%	15.29	144.15	0.97	–
	Hexa	27,072	0.25%	25.92	162.42	–	0.57 (5 elements)
4	Tetra	855,299	0.01%	0.50	179.00	1.00 (48 elements)	–
	Hexa	46,472	0.11%	37.21	146.44	–	0.72
5	Tetra	527,213	0%	2.15	175.70	1.00 (8 elements)	–
	Hexa	25,894	0.06%	43.65	143.47	–	0.75
6	Tetra	624,894	0%	16.42	139.88	0.96	–
	Hexa	26,578	0%	58.70	126.16	–	0.84
7	Tetra	656,219	0%	6.15	142.29	0.97	–
	Hexa	32,200	0%	52.94	130.73	–	0.87
8	Tetra	863,937	0%	0.63	178.75	1.00 (25 elements)	–
	Hexa	44,968	0.02%	43.19	140.19	–	0.74
9	Tetra	745,591	0%	17.16	139.80	0.95	–
	Hexa	49,720	0.19%	38.91	146.40	–	0.67
10	Tetra	864,558	0%	16.35	138.73	0.96	–
	Hexa	52,884	0%	47.97	136.25	–	0.86

Table A3. Solution convergence and computation time for displacement-based (elements C3D10, C3D8, C3D8R, C3D20, and C3D20R) and hybrid (elements C3D10H, C3D8H, C3D8RH, C3D20H, and C3D20RH) finite element (FE) formulations. MPS = maximum principal stress, GCI=grid convergence index (Equation 2), p=order of convergence (Equation 3), and f_estimated= estimated asymptotic solution. Refer to Table A1 for more information about the models. For models 4, 6, 8, 12, 20, and 28, p, GCI and f_estimated could not be obtained due to the solution inconsistency.

Model no.	FE formulation	No. of elements across AAA wall thickness	No. of elements	No. of nodes	Peak MPS (MPa)	Computation time (sec)	p	GCI	festimated
4	C3D10	1	72,507	137,862	1.194662	13	–	–	–
6		2	718,738	1,117,412	1.923035	259
8		4	2,164,776	3,462,038	1.916670	2512
5	C3D10H	1	72,507	137,862	0.555902	17	2.227	0.0433	0.6923
7		2	718,738	1,117,412	0.575515	287
9		4	2,164,776	3,462,038	0.667334	2825
14	C3D8	2	27,972	42,291	0.267113	4	0.748	0.0613	0.3192
22		4	225,676	283,210	0.288184	36
30		8	1,801,840	2,031,084	0.300730	991
15	C3D8H	2	27,972	42,291	0.267113	4	0.748	0.0613	0.3192
23		4	225,676	283,210	0.288184	45
31		8	1,801,840	2,031,084	0.300730	4084
16	C3D8R	2	27,972	42,291	0.264955	3	0.608	0.0831	0.3224
24		4	225,676	283,210	0.284721	30
32		8	1,801,840	2,031,084	0.297690	893
17	C3D8RH	2	27,972	42,291	0.264955	4	0.608	0.0831	0.3224
25		4	225,676	283,210	0.284721	41
33		8	1,801,840	2,031,084	0.297690	4277
10	C3D20	1	3,465	45,095	0.349787	3	0.226	0.2580	0.23529
18		2	27,972	154,734	0.337121	20
26		4	225,676	1,075,083	0.322297	401
11	C3D20H	1	3,465	45,095	0.307018	3	1.950	0.0012	0.3129
19		2	27,972	154,734	0.311382	26
27		4	225,676	1,075,083	0.312511	520
12	C3D20R	1	3,465	45,095	0.312899	3	–	–	–
20		2	27,972	154,734	0.311790	17
28		4	225,676	1,075,083	0.312540	308
13	C3D20RH	1	3,465	45,095	0.311511	3	3.574	0.0003	0.3127
21		2	27,972	154,734	0.311593	23
29		4	225,676	1,075,083	0.312570	568

Table A4. Computation time, peak maximum principal stress (MPS), relative absolute difference in the peak value of MPS between tetrahedral and hexahedral meshes for the studied AAAs. C3D10H is fully integrated quadratic tetrahedron using hybrid formulation. C3D20RH is quadratic hexahedron using hybrid formulation and reduced integration. Relative absolute difference was calculated as the absolute difference between the peak MPS obtained using C3D10H (MPS_C3D10H) and C3D20RH (MPS_C3D20RH) elements divided by the peak MPS obtained using C3D20RH elements: $\frac{\left|{\mathrm{\text{MPS}}}_{\mathrm{C}3\mathrm{D}10\mathrm{H}}-{\mathrm{\text{MPS}}}_{\mathrm{C}3\mathrm{D}20\mathrm{H}\mathrm{R}}\right|}{{\mathrm{\text{MPS}}}_{\mathrm{C}3\mathrm{D}20\mathrm{H}\mathrm{R}}}·100\mathrm{\%}.$

Patient ID	Model no.	Element type of AAA wall	Computation time (sec)	Peak MPS (MPa)	Peak MPS relative absolute difference (%)
1	1	C3D10H	880	0.1574	2.5%
	3	C3D20RH	157	0.1535
2	1	C3D10H	1641	0.2237	2.4%
	3	C3D20RH	296	0.2184
3	1	C3D10H	1196	0.1888	0.1%
	3	C3D20RH	133	0.1886
4	1	C3D10H	3048	0.2141	0.5%
	3	C3D20RH	1355	0.2130
5	1	C3D10H	503	0.2265	3.4%
	3	C3D20RH	152	0.2188
6	1	C3D10H	753	0.2230	1.1%
	3	C3D20RH	181	0.2206
7	1	C3D10H	823	0.1154	1.1%
	3	C3D20RH	221	0.1141
8	1	C3D10H	1795	0.2145	2.4%
	3	C3D20RH	519	0.2094
9	1	C3D10H	1415	0.2087	0.1%
	3	C3D20RH	373	0.2090
10	1	C3D10H	2139	0.1803	4.0%
	3	C3D20RH	501	0.1731