Modeling Techniques and Boundary Conditions in Abdominal Aortic Aneurysm Analysis: Latest Developments in Simulation and Integration of Machine Learning and Data-Driven Approaches

1. Introduction

An abdominal aortic aneurysm (AAA) is formed when the wall of the abdominal aorta weakens, leading to a localized balloon-like structure. By location, AAAs can be classified as supraceliac, juxtarenal, infrarenal and aortoiliac [1]. Infrarenal aneurysms are diagnosed when the diameter of the aorta exceeds 1.5 times of its normal size, approximately about 3 cm [2].The risk of aneurysm rupture increases with the diameter. For the aneurysms larger than 6 cm in diameter, risk of rupture ranges between 10% - 20% annually [3]. However, the diameter cannot solely predict the rupture risk, as clinical studies have documented ruptures in AAAs having diameter less than 5 cm [4]. AAAs develop often without showing any symptoms, resulting in severe cases of AAA being undiagnosed [2].

The aorta wall has thickness around 1.5 mm, composed of three primary layers (intima, media and adventitia), which include smooth muscle cells, collagen and elastin fibers embedded in ground matrix [5]. Their intrinsic mechanical properties enable the aorta to sustain loads at higher pressures. When pathophysiologically examined, an aneurysm can be described as expansion of the vessel wall due to thinning structure and the impact of mechanical stresses. Inflammation, oxidative stress, smooth-muscle-cell death and matrix degradation are the main causes of thin vessel wall structure [6]. Additionally, hemodynamic factors, primarily driven by blood pressure, flow velocity and resistance against flow, mechanically influence the vessel wall. A critical measure of this influence is the wall shear stress (WSS), which represents the frictional force exerted by blood flow along the vessel wall. While WSS acts tangentially to the direction of blood flow, blood pressure exerts perpendicular force to the vessel wall. These combined forces generate peripheral stresses which cause vessels to stretch and deform. The magnitude of these stresses changes throughout the cardiac cycle, making vessel deformation time dependent. Relative wall strain (RWS) quantifies these dynamic effects caused by the peripheral stresses [3].

Blood flow decelerates when it crosses an aneurysm due to the expansion in the cross-sectional area. This deceleration, along with the asymmetric structure and geometric distortions, results in disturbed flow patterns and vortex formation. The thinnest and innermost layer of the vessel wall (intima layer) is a membrane covered with endothelial cells reacting to the disturbed flow and WSS variations by the release of inflammatory markers [5,7]. Consequently, the low velocity regions with vortices and WSS variations promote the accumulation of the platelets in that area, potentially resulting in intraluminal thrombus (ILT) formation that might lead to degradation of the wall. Additionally, the areas of higher pressure indicate where the aneurysm may expand [8].

Numerous modeling and simulation studies have been commonly used in recent years to explore the hemodynamic and biomechanical properties of AAA, and the rupture risk. Rupture tends to occur in areas of low WSS and flow recirculation, rather than in regions of high pressure and WSS which conflict with earlier finite element analysis (FEA) based predictions [9]. AAA wall stress distribution is significantly influenced by aneurysm asymmetry and wall thickness variations. Models with variable wall thickness have shown up to a four-fold increase in peak wall stress compared to those with uniform thickness. Accuracy of wall thickness in computational models to predict AAA rupture risk emphasized by combining FEA and fluid-structure interaction (FSI) [10]. Patient-specific FSI models of AAA have proved that wall stress distribution and WSS are thoroughly affected by the interaction between blood flow and the arterial wall. These models show that high wall stress regions often align with complex flow patterns which carry higher risk in the context of rupture. Considering fluid dynamics via those interactions offers a more realistic simulation of in vivo conditions compared to models that do not include the flow effects [11].

Accurate modeling of AAA is clinically relevant due to its impact on predicting rupture risk, guiding surgical decisions, and improving patient outcomes. The integration of advanced modeling techniques into clinical practice helps avoid unnecessary surgical intervention and improves risk assessment. Rigorous biomechanical analysis is crucial for clinically relevant detection of AAA rupture risk via computational assessment. FSI models that incorporate patient-specific geometries with accurate wall thickness and boundary conditions derived from advanced imaging techniques such as phase-contrast MRI provide higher accuracy in predicting arterial wall stress and strain, which helps to recognize the smaller aneurysms with a high risk of rupture [12,13]. Integration of the machine learning algorithms to the geometric and biomechanical markers significantly enhances the classification accuracy of symptomatic and asymptomatic AAA. It is stated that the generalized additive model classifier achieved the highest accuracy using the combination of seven markers, demonstrating that integrating patient-specific geometric and biomechanical properties provides a more accurate rupture risk assessment than relying solely on the maximum aneurysm diameter [14]. Exceeding the traditional metrics such as maximum diameter [15] lies in bringing the interaction of hemodynamics, wall mechanics, geometric features into play with several modeling and simulation techniques. It is considered that integrating patient-specific data into these models will likely enhance the precision of AAA rupture predictions with the developing computational power and machine learning algorithms, and finally it will improve clinical decision-making and patient outcomes.

In summary, the recent advancements in computational techniques have enabled a more comprehensive understanding of the biomechanical and hemodynamic parameters associated with AAA rupture risk. In computational models, the use of realistic and sophisticated BCs plays a pivotal role in capturing the realistic physiological environment within the AAA. Accurate representation of these BCs at the inlet and outlet of the fluid domain, along with precise modeling of wall mechanics, is essential for enhancing the reliability of simulations. Such developments not only refine predictions of arterial stress and strain but also support the integration of emerging technologies, such as machine learning algorithms, into AAA analysis. These advancements bridge the gap between theoretical models and their practical applications in clinical settings, promoting a more personalized approach to patient care.

This paper provides a comprehensive review of the diverse range of BCs employed in AAA simulations, focusing on their application at the inlet and outlet of the fluid domain as well as wall conditions. Realistic BCs, such as the Womersley profile, Windkessel model, and fluid-structure interaction (FSI) techniques, are critical for accurate AAA simulations but are often absent in standard computational fluid dynamics (CFD) solvers due to their technical complexity and the need for a detailed understanding of the underlying clinical issues. This review aims to address these challenges by exploring the physical and mathematical principles behind these BCs, offering valuable insights to researchers. Additionally, key findings from the literature are summarized, highlighting the limitations and technical barriers of existing models. The paper concludes by discussing recent advancements in predicting hemodynamic and biomechanical parameters and disease progression through machine learning (ML) and data-driven techniques, with a focus on commonly used algorithms, their limitations, and potential future directions.

2. Computational Techniques for AAA Simulations

In the modeling of cardiovascular systems, various techniques have been developed. These can be categorized into low-dimensional (0D and 1D) [16,17,18,19] and high-dimensional (2D and 3D) [20,21,22,23,24] methods, based on the dimensionality of the simulation [18,25,26,27]. Low dimensional techniques are computationally affordable and suitable for capturing general flow features across the entire cardiovascular system [28]. Specifically, 1D models solve the Navier-Stokes equations along the arterial tree to predict the flow profiles and pressure waves as they propagate throughout the system [22,28,29,30,31,32,33]. High dimensional models can resolve detailed hemodynamic and biomechanical features in all spatial dimensions by solving governing equations in the solid and fluid domains. AAA has a unique and complex geometry with a wall thickness around 1.5 mm [34]. Therefore, 3D patient-specific simulations are essential to achieve high-fidelity hemodynamic and biomechanic predictions.

2.1. Patient Specific Modeling Approaches

Patient-specific modeling increases diagnostic and treatment accuracy by providing the examination adapted to individual anatomical differences. Non-invasive imaging methods such as magnetic resonance imaging (MRI) [35,36,37,38], time-resolved three-dimensional ultrasound (3D+t US) [39,40] and computed tomography angiography (CTA) are the typically used to extract patient-specific geometries for AAA [41,42,43,44,45]. When the non-invasive methods are inadequate to determine the patient specific geometry, invasive methods can be applied to provide more detailed information as imaging is performed in a more limited area. For instance, intravascular ultrasound (IVUS) is performed by inserting a catheter containing an ultrasound probe into a vein, allowing the probe to capture images within the vessel. The use of IVUS, especially in cases of AAAs, eliminates the use of contrast agents used in CTA. Patients who underwent CTA during endovascular repair (EVAR) procedures have shown deteriorating renal function over time [46]. Angiography methods are also important tools that can be used to model the geometric structure. For example, 3D rotational angiography (3DRA), which allows visualization of aneurysm structures of the order of 3 mm [47], and other traditional angiography methods can be applied for the aneurysms emerging in larger vessels such as the aorta.

Following patient-specific medical image acquisition in DICOM (Digital Imaging and Communications in Medicine) format, segmentation software such as MIMICS (Materialize, Leuven, Belgium), SimVascular, VMTK and ITK-SNAP reconstruct 3D models. Deep learning tools such as U-Net 3D [48], SegNet [49], 3DResUNet [50], attention-based residual U-Net (ARU-Net) [51] and Context-Aware Cascaded U-Net (CACU-Net) [52] architectures offer significantly faster segmentation process compared to traditional algorithms. The overall workflow of segmenting is represented in Figure 1 Lumen and AAA walls are segmented separately to differentiate between these zones. The range of AAA thickness is within 0.23-4.26 mm, with a median wall thickness of 1.48 mm [34]. The distribution of wall thickness on AAA is heterogeneous and changes from patient to patient. The geometrical structure of the wall has a significant importance to withstand the dynamic loads and careful segmentation of the wall is an important step in computational modeling.

Extracting comprehensive CTA images of the entire vascular network can introduce additional challenges. Moreover, large-scale image-based simulations of the entire arterial network are computationally demanding [16,54]. The common approach is to focus on a specific region of interest by cutting the aorta at certain boundaries. As illustrated in Figure 2, the region of interest is specified and truncated from the rest of the geometry by cutting the domain at the inlet and outlet boundaries during segmentation process [13,55,56,57,58,59,60,61,62,63]. Several segmented 3D AAA models are presented in Figure 3. The supraceliac region (SC) marks the boundary between the descending and abdominal aorta. Downstream from the SC, the aorta branches into the celiac trunk (CT), superior mesenteric (SM), left renal (LR), right renal (RR) and accessory right renal (ARR) arteries. The aorta segment below the renal arteries is called infrarenal region (IR), which bifurcates into the left central iliac (LCI) and right central iliac (RCI) arteries [37,64,65].

The placement of inlet boundary varies depending on the hemodynamic parameters under investigation. Some studies position the inlet by dividing the region of interest from the aorta at the SC [65,66,67], while others cut the aorta from the IR [10,68,69], depending on the type of aneurysm. The number of outlet boundaries varies based on the position of the inlet. When the inlet is positioned in the SC, the model typically includes at least six outlets, corresponding to the CT, SM, LR, RR, LCI and RCI arteries [37,64,65]. In contrast, placing the inlet in the IR reduces the number of outlets to at least two, located at LCI and RCI arteries. However, the number of outlets may increase if additional arterial branches supplying blood to organs are considered. For example, the CT bifurcates into the hepatic (H) and splenic (S) arteries, while central iliac arteries further divide into the external iliac (LEI and REI) and internal iliac (LII and RII) arteries. Therefore, the total number of outlets may vary, depending on the complexity and specific requirements of the simulation. In the following section, numerical analysis is explained by providing details of the governing equations, material properties and boundary conditions used in solid and fluid domains.

2.2. Analysis of the Fluid Domain

Computational fluid dynamics (CFD) is a numerical technique used to determine the main hemodynamic parameters, such as pressure and velocity fields of the overall fluid domain. In this technique, the following governing equations are numerically solved across a computational mesh to characterize the flow properties as given in Eq. (2.1) and Eq. (2.2).

$$\nabla .\mathit{v}=0$$

(2.1)

$${\rho }_f{\displaystyle \frac{\partial \mathit{v}}{\partial t}}+{\rho }_f\left(\left(\mathit{v}-\mathit{w}\right).\nabla \right)\mathit{v}=\nabla .{\mathit{\tau }}_{\mathit{f}}+{\mathit{f}}_{\mathit{f}}$$

(2.2)

Eq. (1) and Eq. (2) are continuity and Navier-Stokes equations, respectively. The term denoted by ${\mathit{f}}_{\mathit{f}}$ is the external force acting per unit volume of fluid. The parameters of $\mathit{v},{\rho }_f$ and${\mathit{\tau }}_{\mathit{f}}$ define the velocity vector, density and stress tensor of the fluid, respectively. The parameter of $\mathit{w}$ denotes the moving coordinate velocity [69]. In ALE formulation, $\mathit{v}-\mathit{w}$ is denoted as the relative velocity of the fluid with respect to the moving coordinate velocity [70,71]. Both pressure and shear stresses are calculated for determining the stress tensors as given in Eqn. (2.3).

$${\mathit{\tau }}_{\mathit{f}}=-p\mathit{I}+2\mu \left(\dot{\gamma }\right)\mathit{\epsilon }$$

(2.3)

$\mathit{I}$defines the identity tensor, $p$ denotes the pressure and $\mu$ denotes the dynamic viscosity. The parameter of $\dot{\gamma }$ defines the shear rate and $\mathit{\epsilon }$ denotes the strain rate tensor given in Eq. (2.4).

$$\mathit{\epsilon }={\displaystyle \frac{1}{2}}\left(\nabla \mathit{v}+{\nabla \mathit{v}}^T\right)$$

(2.4)

Constitutive Models for Blood

As a first approximation, blood can be assumed as a Newtonian fluid with a constant dynamic viscosity across all shear rates ($\mu \left(\dot{\gamma }\right)= 3.5\mathrm{c}\mathrm{P}$), and a constant mass density (${\rho }_f=1.05 g/{cm}^3$) [8,41,68,72]. However, at low shear rates, particularly below 100 s^-1, the non-Newtonian behavior of blood becomes more prominent [73]. At low shear rates, the red blood cells (RBCs) aggregate and form rouleaux, which is rod shaped stacks of individual cells [73]. Rouleaux aggregation disperses as the shear rate increases, reducing the viscosity of blood. The resulting shear-thinning behavior caused by rouleaux disaggregation in blood plasma is the principal cause of the non-Newtonian behavior of blood. In numerical studies, shear-thinning behavior of working fluid, which is blood, can be modeled using several constitutive models. Viscosity functions with bounded and non-zero limiting values of viscosity can be written in the general form given in Eq. (2.5.)

$$\mu \left(\dot{\gamma }\right)={\mu }_{\mathrm{\infty }}+\left({\mu }_0-{\mu }_{\mathrm{\infty }}\right)F\left(\dot{\gamma }\right)$$

(2.5)

In Eq. (2.5), ${\mu }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\mu }_{\infty }$ are the asymptotic viscosity values at zero and infinite shear rates and $F\left(\dot{\gamma }\right)$ is a shear dependent function, satisfying the following natural limit conditions given in Eq. (2.6) and Eq. (2.7).

$$F\left(\dot{\gamma }\right)={\displaystyle \frac{\mu \left(\dot{\gamma }\right)-{\mu }_{\mathrm{\infty }}}{{\mu }_0-{\mu }_{\mathrm{\infty }}}}$$

(2.6)

$$\underset{\dot{\gamma }\to 0}{\mathit{lim}}F\left(\dot{\gamma }\right)=1,\underset{\dot{\gamma }\to \mathrm{\infty }}{\mathit{lim}}F\left(\dot{\gamma }\right)=0$$

(2.7)

Different choices of the function $F\left(\dot{\gamma }\right)$ correspond to different constitutive models for blood, with material constants depending on the factors such as temperature, hematocrit, and plasma. In the literature, eight non-Newtonian models are commonly used to represent the shear-thinning behavior of blood: Carreau, Carreau-Yasuda [74,75,76,77,78], Casson, Quemada [79], Power, Cross, Simplified Cross, and Modified Cross [78,80,81,82,83,84,85,86,87,88], as summarized in Table A in Appendix A.

On the other hand, the behavior of RBCs cannot be fully characterized by viscous phenomenon because they can be regarded as fluid filled elastic cells [89]. Therefore, blood exhibits viscoelastic nature due to RBCs viscous shear thinning and elastic properties [90]. In addition, increasing hematocrit, which is the volumetric ratio of RBCs in blood, causes increasing viscosity of the blood and makes non-Newtonian behavior of blood more significant. However, as the shear rate increases beyond the low shear rate region, the shear-thinning characteristics diminish, and blood exhibits Newtonian behavior [73].

To account for the elastic nature of blood, the viscoelastic Oldroyd-B model is frequently implemented [91,92,93]. In that model, the stress tensor in the linear momentum equation is decomposed into viscous and elastic parts as${\mathit{\tau }}_{\mathit{f}}={\mathit{\tau }}_{\mathit{v}}+{\mathit{\tau }}_{\mathit{e}}$ [73,93]. The viscous part of the Oldroyd-B model ${\mathit{\tau }}_{\mathit{v}}=2{\mu }_v\left(\dot{\gamma }\right)\left(\nabla v+{\nabla v}^T\right)$ can be Newtonian or shear-thinning [90]. The elastic part satisfies the following constitutive equation given in Eq. (2.8).

$${\mathit{\tau }}_{\mathit{e}}+{\lambda }_1\left({\displaystyle \frac{\partial {\mathit{\tau }}_{\mathit{e}}}{\partial t}}+\mathit{v}·\nabla {\mathit{\tau }}_{\mathit{e}}-{\mathit{\tau }}_{\mathit{e}}·\nabla \mathit{v}-{\nabla \mathit{v}}^T·{\mathit{\tau }}_{\mathit{e}}\right)={\mu }_e\left(\nabla \mathit{v}+{\nabla \mathit{v}}^T\right)$$

(2.8)

In Eq. (2.8), ${\mu }_e$ denotes the elastic viscosity coefficient, and ${\lambda }_1$ denotes the relaxation time [87]. For the blood, ${\mu }_e$ is used as $0.0004 \mathrm{P}\mathrm{a}·\mathrm{s}$and ${\lambda }_1$ is used as $0.06$ s [87].

2.3. Analysis of Solid Domain

The finite element analysis (FEA) is a numerical method that allows the calculation of stresses and strains within the solid region to understand mechanical behavior of the solid structure under various conditions. The governing equation provided in Eq. (2.9) is numerically solved across a mesh created through the entire solid domain.

$$\nabla .{\mathit{\tau }}_{\mathit{s}}+{\mathit{f}}_{\mathit{s}}={\rho }_s{\displaystyle \frac{{\partial }^2\mathit{u}}{\partial t^2}}$$

(2.9)

In Eq. (2.9), ${\mathit{\tau }}_{\mathit{s}}$ defines the solid stress tensor, ${\mathit{f}}_{\mathit{s}}$ defines the body forces per unit volume, and $\mathit{u}$ defines the structural displacements. By writing the stress variable in terms of strains using Hooke's law, the differential equation can be written in terms of displacements, as given in Eq. (2.10).

$${\mathit{\tau }}_{\mathit{s}}=\mathit{C}\mathit{\epsilon }$$

(2.10)

In Eq. (2.10), the stiffness tensor is defined as $\mathit{C}$, the solid strains are defined as $\mathit{\epsilon }$, and a linear correlation is defined between the stress and strain. However, in many cardiovascular modeling studies, the solid region to be investigated does not exhibit a linear elastic behavior [94]. Modeling the AAA tissue requires proper defining the mechanical properties of wall and ILT. The constitutive models of AAA wall are developed, calibrated and validated through tissue characterization experiments [95,96,97,98]. To fit stress-strain behaviors, uniaxial tensile tests are conducted using tissue samples gathered during surgery [94]. AAA wall tissue exhibits complex material properties due to its layered, fiber-oriented structure, including hyperelasticity (non-linear stress-strain behavior), viscoelasticity (time-varying response due to relaxation), and anisotropy (direction-dependent characteristics due to distribution of the collagen fibers) [15,99].

2.3.1. Constitutive Models for AAA Wall

As a first approximation, the aneurysm wall can be modeled as a single layer with linear elastic properties and isotropic medium characteristics using elastic modulus of 2.7 MPa, Poisson’s ratio of 0.45, and mass density 2000 $\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$ [57,58,100,101,102]. The elastic modulus ratio for three-layered models of intima: media: adventitia is $~1:3:2$ [103,104] with elastic modulus of 3.2, 4,8 and 1.6 MPa, respectively [105,106].

While linear elastic description can be applied for the healthy aorta or the aneurysm with low curvature [107], effective characterization of non-linear stress-strain properties of AAA wall depends on the tissue's hyperelastic nature [97]. Fitting the hyperelastic constitutive models to the experimental data enables to determine the model coefficients. Although collagen fiber number and orientation influence the mechanical behavior, modeling with isotropic strain energy function is common [39,56,77,89,108]. To improve modeling accuracy, constitutive models including the anisotropy are available, capturing the directional dependency of collagen fibers in the strain energy function [95,96,97,98,107,109,110,111].

• Hyperelastic and Isotropic Models
• Mooney Rivlin Model

In the literature, two- [13,108,112,113] and three- [106] parameter Mooney Rivlin models are commonly used to model non-linear stress-strain nature of AAA tissue and described as Eqs. (2.11), (2.12), and (2.13), respectively. The strain energy per unit volume $\left(W\right)$ is defined as

$$W=C_{10}\left(I_1-3\right)+C_{20}{\left(I_1-3\right)}^2+{\displaystyle \frac{{\left(J-1\right)}^2}{D}}$$

(2.11)

$$I_1=tr\mathit{C}={\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2$$

(2.12)

$$J={\lambda }_1{\lambda }_2{\lambda }_3$$

(2.13)

In Eq. (2.11), $C_{10}$ and $C_{20}$ are the material constants obtained by fitting the experimental data. $I_1$ is the first Green’s strain invariant, $\mathit{C}$ is the right Cauchy-Green strain tensor, and ${\lambda }_1, {\lambda }_2,{\lambda }_3$ are the principal strains that determine the relative changes in length along each principal direction. $J$ is the ratio of deformed elastic volume to the undeformed volume of the material. For incompressible materials, $J=1$due to zero change in material volume, that vanishes the last term of Eq. (2.11). $D$ is the material incompressibility term.

Considering the population averages, [99,114,115] material properties are determined as $C_{10}=17.4 \mathrm{N}/{\mathrm{c}\mathrm{m}}^2$, $C_{20}=188.1 \mathrm{N}/{\mathrm{c}\mathrm{m}}^2$ for a single layered wall structure. However, the material constants show certain variations from patient to patient. The minimum, average, and maximum values of $C_{10}$​ were reported as $15.2, 17.4 \mathrm{a}\mathrm{n}\mathrm{d} 21.9 \mathrm{N}/{\mathrm{c}\mathrm{m}}^2$, respectively, while the corresponding values for $C_{20}$ as 117.6, 188.1, 355.7 $\mathrm{N}/{\mathrm{c}\mathrm{m}}^2$ [56,116]. These studies recommend using population-averaged values in simulations, as variations in the Mooney-Rivlin material constants have been shown to have a negligible impact on peak wall stress and wall displacement [107].

• Yeoh Model

Generalized Yeoh strain energy density function, $W,$ is given in Eq. (2.14).

$$W=\sum _{i=1}^N{C}_{i0}{\left(I_1-3\right)}^i$$

(2.14)

In Eq. (2.14), $I_1$ is the first invariant of the right Cauchy-Green strain, and $C_{i0}$ is stress-like material parameters identified from experimental data [79,117]. For the second order Yeoh model, $W$ is defined as given in Eq. (2.15).

$$W=C_{10}\left(I_1-3\right)+C_{20}{\left(I_1-3\right)}^2$$

(2.15)

In Eq. (2.15), $C_{10}$ is used as $17.4 \mathrm{N}/{\mathrm{c}\mathrm{m}}^2$, and $C_{20}$ is used as $188.1 \mathrm{N}/{\mathrm{c}\mathrm{m}}^2$ [94]. For the fifth order Yeoh model, the material values are calculated from planar biaxial experimental data as $C_{10}=0.5 \mathrm{N}/{\mathrm{c}\mathrm{m}}^2$, $C_{20}=C_{30}=0,$ $C_{40}=370\mathrm{N}/{\mathrm{c}\mathrm{m}}^2$, and $C_{50}=1374 \mathrm{N}/{\mathrm{c}\mathrm{m}}^2$.

• Hyperelastic and Anisotropic Models

To account for the anisotropic nature of the AAA tissue, Fung-type strain energy density function [118] is commonly used in the literature, which can be demonstrated as given in Eqs. (2.16), (2.17), (2.18), and (2.19).

$$W=K\left(e^Q-1\right)$$

(2.16)

$$Q=C_{11}{E}_1^2+C_{22}{E}_2^2+C_{12}{E}_1{E}_2$$

(2.17)

$$E_1={\displaystyle \frac{1}{2}}\left({\lambda }_1^2-1\right)$$

(2.18)

$$E_2={\displaystyle \frac{1}{2}}\left({\lambda }_2^2-1\right)$$

(2.19)

In Eqs. (2.16) and (2.17), $K$ is a material constant related to the stiffness of the tissue, while $C_{11}$, $C_{22}, C_{12}$ are the material constants depending on the specific tissue directions such as circumferential, longitudinal, and radial to capture anisotropy [107,119]. In Eqs. (2.18) and (2.19), $E_1$ and $E_2$ are the Green strains in circumferential and longitudinal directions, while ${\lambda }_1 \mathrm{a}\mathrm{n}\mathrm{d} {\lambda }_2$ are the principal stretches that determine the relative changes in length along two primary axes. By fitting the material parameters to the experimental data, the effect of anisotropic fiber orientations can be modeled. In addition to the Fung-type [120], the models offered by Holzapfel [98], Holzapfel and Gasser [121] and Choi and Vito [122] are commonly applied in literature. Most of the studies combine the isotropic and anisotropic strain energy density functions as given in Eq. (2.20) [94]

$$W= W_{ISO}+ W_{ANISO}$$

(2.20)

In Eq. (2.20), $W_{ISO}$ describes the energy stored withing the extracellular matrix, and $W_{ANISO}$ demonstrates the effect of embedded collagen fibers.

• Constitutive Models for ILT

While ILT consists of two layers [123], with the inner layer forming about one-third and the outer layer about two-thirds of its thickness, modeling it as a single-layer elastic structure is common [56,108,124]. The elastic modulus of ILT is relatively lower than that of the wall. When ILT is initially formed, elastic modulus is nearly equal to 50 kPa, while it may reach up to 200 kPa due to the increasing rigidity of ILT structure with time [56]. Poisson’s ratio and mass density of ILT are generally used as 0.45 and 1050 $\mathrm{k}\mathrm{g}/{\mathrm{c}\mathrm{m}}^3$, respectively. The tensile tests have shown that ILT tissue has an isotropic behavior [107]. To increase accuracy, nonlinear, hyperelastic and isotropic constitutive model for ILT developed by Martino et al. [123] can be used as given in Eqs. (2.21) and (2.22).

$$W=C_{01}\left(I_2-3\right)+C_{02}{\left(I_2-3\right)}^2$$

(2.21)

$$I_2={\displaystyle \frac{1}{2}}\left[{\left(tr\mathit{C}\right)}^2-\left(tr{\mathit{C}}^2\right)\right]={\lambda }_1^2{\lambda }_2^2+{\lambda }_2^2{\lambda }_3^2+{\lambda }_3^2{\lambda }_1^2$$

(2.22)

In Eqs. (2.21) and (2.22), $C_{01}$ and $C_{02}$ are the material parameters, and $I_2$is the second invariant of the left Cauchy-Green stretch tensor. For the ILTs structure, $C_{01}=3.37$ and $2.23 \mathrm{N}/{\mathrm{c}\mathrm{m}}^2$, $C_{02}=3.47$ and $2.24 \mathrm{N}/{\mathrm{c}\mathrm{m}}^2$ for inner and outer layers, respectively. By considering ILT as a single layer, material constants can be taken as $C_{01}=2.80$ and $C_{02}=2.86 \mathrm{N}/{\mathrm{c}\mathrm{m}}^2$. Martino et al. [123] conducted a comparative study using various parameters and concluded that suggested mean parameters are sufficient, eliminating the need for patient-specific material constants.

2.4. Coupling of Solid and Fluid Domains: Fluid-Structure Interaction (FSI)

Modeling the AAA wall as rigid and stationary is common in literature due to its simplicity. In fact, the vascular layers are viscoelastic and subjected to large deformations by the blood flow generated loads. Conversely, the motion of hyperelastic walls alters the shape of fluid domain and hemodynamics within AAA, which necessitates a simultaneous blood flow and wall deformation treatment. Compared to flexible walls, rigid wall assumption obtains elevated WSS values up to 50% [125], underestimates the vortex development [109], and fails to predict pressure wave propagation [126,127]. Therefore, the mutual interaction among wall deformation and blood hemodynamics requires an FSI coupling at the aorta wall.

FSI is a multiphysics technique [71,125,128] that couples fluid domain to solid domain at the wall interface by exchanging the instantaneous fluid velocity and wall deformation values between domains. The domains do not overlap, and the two domains are coupled at the luminal surface by satisfying a set of physical interface conditions to ensure the compatibility of kinematics and tractions. The solid wall domain and ILT structure are modeled by finite-element analysis (FEA) by discretizing the momentum conservation given in Eq. (2.9), while CFD method is used to discretize the Navier-Stokes Equations given in Eq. (2.2) that resolves blood flow patterns. The solid domain moves through space and deforms under hemodynamic loads, altering the shape of the fluid domain and necessitating the use of the Arbitrary Lagrangian-Eulerian (ALE) formulation. The ALE description of fluid domain offers a moving coordinate system to simulate a non-stationary solution domain [129]. In ALE method, the mesh nodes may be displaced or be fixed as Lagrangian and Eulerian descriptions. With this freedom, large distortions within the fluid domain can be managed with more precision [129].

Solving the governing equations within such truncated solid and fluid domains, coupling the solutions through FSI procedure requires proper specification of the conditions at the boundaries of the domain, bringing concerns regarding the accuracy of the simulations. A high-fidelity, physiologically realistic modeling of AAA hemodynamics requires proper treatment to the conditions at the model boundaries. In Figure 4 and Table 1, boundaries are demonstrated for the fluid and solid domains. Appropriate specification of BCs at the inlet, outlets and walls of solid and fluid domains is essential to realistically couple velocity, pressure, stress distribution, kinematics, traction and wave propagation in the upstream and downstream vasculature. The typical modeling techniques commonly employed in the literature for these boundaries are investigated in the following section.

3. Modeling Boundary Conditions in Fluid Domain

The accurate prediction of hemodynamic parameters, such as pressure and flow rate, requires proper coupling of the AAA domain with upstream and downstream hemodynamics at the truncated inlet and outlet boundaries. Proper coupling ensures that the simulated AAA hemodynamics interact with adjacent arterial compartments, closely reflecting physiological conditions. Typically, the time-varying inlet flow rate given in Figure 5 is prescribed for the inlet located at IR [56,68,109]. To enhance model accuracy, patient-specific inlet flow rates can be obtained through Doppler ultrasonography or phase-contrast magnetic resonance imaging (PC-MRI) [130]. At each point, spatial variations in flow velocity can be incorporated through customized velocity profiles based on the time-dependent flow rate. These profiles can be employed as idealized (i.e., flat, parabolic or Womersley) or patient specific. Patient-specific inlet velocity profiles can be measured non-invasively using 4D-Flow MRI [13,130,131], a three-dimensional, time-resolved type of PC-MRI that measures and visualizes the temporal evolution of blood flow within a specific 3D volume [78,81,131]. Alternatively, the inlet boundary can be connected to a lumped description of the heart [132,133,134]. However, using time-dependent pressure at the inlet is not a common approach due to the invasive nature and low precision of intraluminal pressure measurement, which requires catheter placement inside the inlet artery segment.

Applying outlet boundary conditions (BCs) is typically more complex. One challenge is the clinical measurement of intraluminal pressure, which requires catheter placement [13]. However, a significant portion of the cardiovascular system, including capillaries and veins, are located downstream of the aorta and substantially influence the AAA hemodynamics. Acquiring in-vivo hemodynamic parameters at branching arteries using non-invasive techniques like magnetic resonance imaging (MRI) is challenging due to limitations in the spatial and temporal resolutions of the measuring devices and the physiological variations in patients during the imaging [135]. Van’t Veer et al. [136] compared non-invasive brachial cutoff blood pressure measurements with invasive catheter pressure measurements inside the AAA sac, reporting a 5% underestimation of systolic pressure and a 12% overestimation of diastolic pressure. The common approach, especially for geometries with one or two outlets, is to prescribe the time-varying pressure profile illustrated in Figure 5 at the outlets [34]. Various other techniques are available in the literature, including prescribed outlet pressure, flow-split method, lumped parameter models (e.g., 3-element Windkessel), and 1D distributed parameter techniques. The following section critically examines the common methods for defining boundary conditions, their applications, benefits, and drawbacks.

3.1. Inlet BCs

The time-varying flow rate at the inlet of an AAA can be integrated using customized velocity profiles, capturing spatial variations in flow velocity across the inlet section. The pattern of these flow rate waveforms is highly dependent on the inlet location. Depending on the requirements, the inlet can be positioned at the SC or IR locations. Figure 6 demonstrates the measured volumetric flow rates at SC and IR regions [54]. Patient-specific profiles measured from human aorta using 4D flow PC-MRI technique, offer a more accurate representation of individual blood flow patterns, while idealized profiles (e.g., flat, parabolic, or Womersley) can illustrate general flow behaviors and are often used when patient data is unavailable [13,131,137].

The 4D flow PC-MRI images capture blood velocity vector fields in multiple sagittal slices of a healthy human aorta [131,138,139] at various time instants during cardiac cycle. In this method, three-dimensional velocity maps at multiple time instants through the cardiac cycle are generated. The resulting pixel-based time-varying velocity vectors are imposed on each voxel at the inlet section [139]. Chandra et al. [13] proposed a technique for mapping velocity vector field data onto the inlet boundary of a patient-specific AAA model, acquired on the IR plane. The inlet velocity profiles derived from 4D flow PC-MRI have a boundary which changes shape, size, and position due to aortic expansion and contraction during systolic and diastolic phases. However, the cross section of the inlet boundary in the CFD models of AAA geometry remains fixed both temporally and spatially [137]. To address this mismatch, Schwarz-Christoffel mapping aligns the datasets. Directly measuring in vivo inlet velocity profiles is still challenging because of the mismatch generated by such cardiac motion and resolution of measuring devices [140,141]. Moreover, accessing complete high-quality patient-specific geometry and inlet profile data is not always possible due to limited imaging facilities [141,142].

Consequently, patient-specific velocity profiles are not commonly used as inlet BC in AAA simulations. Instead, many studies employ artificial profiles such as flat [58,143], parabolic [8,20,72,89,144,145] and Womersley [146,147,148,149]. Figure 7 shows the general pattern of the idealized velocity profiles supplied at the inlet of AAA models. Several studies have reported no significant difference between flow-MRI derived profiles and artificial ones [150].

The artificial velocity profiles are calculated from time-dependent, patient-specific volumetric flow waveforms represented in Figure 6 at the inlet section of the computational domain. Inlet flow rate waveforms are implemented to calculate the flat and parabolic velocity profiles using Eqs. (3.1) and (3.2)

$$U_{flat}\left(t\right)=U\left(t\right)={\displaystyle \frac{2Q\left(t\right)}{\pi R^2}}$$

(3.1)

$$U_{par}\left(r,t\right)={\displaystyle \frac{2Q\left(t\right)}{\pi R^2}}\left(1-{\left({\displaystyle \frac{r}{R}}\right)}^2\right)$$

(3.2)

where $U_{flat}$ and $U_{par}$ are flat and parabolic velocities, $Q\left(t\right)$ is the time-dependent flow rate, and $R$ is the artery radius. As presented in Figure 7, the flat profile represents a uniform spatial distribution of velocity vector at the inlet interface. The parabolic profile, derived from Poiseuille’s equation, exhibits a parabolic spatial distribution.

These profiles cannot fully capture all the transient effects introduced by the physiological flow rate waveforms, including flow reversal regions 161]. In 1955, Womersley developed an exact solution for incompressible, Newtonian fluid flow through a cylindrical, rigid blood vessel. This solution considers a pressure gradient that is a periodic function of time which enables to capture reverse flow using the equations of motion and continuity. Figure 7 illustrates the velocity profile obtained by Womersley’s solution. To obtain the Womersley profile, it is necessary to write the time-dependent flow rate, $Q\left(t\right)$, in the harmonic form as in Eq. (3.3)

$$Q\left(t\right)=\sum _{n=0}^N{C}_n{e}^{in\omega t}$$

(3.3)

In Eq (3.3), $N$ is the total number of harmonic coefficients [21,150]. $C_n$’s are the Finite Fourier Transform (FFT) coefficients of that flow rate, while the term $n=0$ $\left(C_0\right)$ corresponds to a steady pressure gradient [150]. Fourier series decomposition of the time-dependent flow rate should be performed to obtain the Fourier coefficients of the inlet flow rate waveform. The equation of Womersley velocity profile is given in Eq. (3.4) using inlet flow rate waveform.

$$U_{wom}\left(r,t\right)={\displaystyle \frac{2C_0}{\pi R^2}}\left(1-{\left({\displaystyle \frac{r}{R}}\right)}^2\right)$$

$$+\sum _{n=1}^N{\displaystyle \frac{C_n}{\pi R^2\left(1-\frac{2J_1\left(i^{\frac{3}{2}}{\alpha }_n\right)}{i^{\frac{3}{2}}{\alpha }_n{J}_0\left(i^{\frac{3}{2}}{\alpha }_n\right)}\right)}}\left[1-{\displaystyle \frac{J_0\left({\alpha }_n\frac{r}{R}{i}^{3/2}\right)}{J_0\left({\alpha }_n{i}^{3/2}\right)}}\right]e^{i{\omega }_{n}t}$$

(3.4)

In Eq. (3.4),${\alpha }_n$ and ${\omega }_n$ are the $n^{th}$ term of Womersley number and frequency, respectively. $J_0$and $J_1$ are the Bessel function of the first kind of order zero and first, respectively [152]. In Supplementary Material, derivation of Womersley profile and MATLAB code are provided.

In Eq. (3.4), the first term on the right-hand side equals the steady Poiseuille equation while the second term is the harmonic [152]. Although the Womersley profile presents transient effects of physiological flow, especially for large Womersley number ($\alpha$) values, its application and implementation as an inlet boundary condition can be challenging due to Bessel functions and imaginary numbers that it contains [153,154]. Consequently, most studies in literature utilize flat or parabolic profiles [155]. Current studies [21] show that Womersley profile can be obtained by providing sufficient entrance length, at least $L_{ent}=3D$, and $10D$ is enough for parabolic and flat profiles, respectively.

Helical flow patterns in aortic hemodynamics have been observed in thoracic aorta studies [138,156]. Blood flow forms helical patterns in the ascending and descending aorta, as well as in the upstream sections of AAA. This helical blood flow is a physiological characteristic where blood rotates and advances along the aortic axis [9,156]. These physiological helices potentially help maintain WSS within normal range [102,157], regulate flow, and protect vessels from thrombus deposition [158]. However, most AAA studies neglect the helical features of the incoming flow. To address this, Javadzadegan et al. [102,157] introduced a tangential velocity profile at the inlet to present the helical pattern of the incoming flow as presented in Eq. (3.5)

$$v_{\theta }=\omega r$$

(3.5)

$$\omega ={\displaystyle \frac{U}{R}}C$$

(3.6)

In Eq. (3.5) $U$ is the streamwise velocity, $v_{\theta }$ is the local tangential velocity component, $r$ and $R$ represent axial axis and undilated radius, respectively. $\omega$ is the spiral speed and $C$ is a constant that controls the magnitude of the spiral speed. In the literature, $C$ is typically defined as 1/6 for the aortic flows [159].

3.2. Outlet BCs

The purpose of the outlet BCs is to model the downstream vasculature that includes smaller arteries, arterioles, capillaries, venules, and veins that return blood to the heart [160]. Therefore, the choice of outlet BCs has significant influence on velocity and pressure fields, and wave propagation linked to the wall properties in 3D AAA simulations [161]. Different outlet BCs have been adopted to better produce in vivo hemodynamic conditions in the cardiovascular system including prescribed outlet pressure, flow-split method [162], lumped parameter models such as Windkessel, and resistance [24,25].

3.2.1. Prescribed Outlet Pressure

Prescribing specific pressure at the outlets is common in AAA simulations. Prior studies typically implemented zero-gauge pressure at the outlets by setting $P_{out}=P_{atm}$ [163,164,165], and some current studies continue to use atmospheric pressure at the outlet for simplicity [147,149,166,167,168,169,170,171,172,173,174]. However, using zero-gauge pressure as an outlet boundary condition is insufficient for producing accurate flow and pressure features. This approach assumes that the outlet is open to the atmosphere, neglecting the effects of the posterior vasculature.

To avoid this unrealistic assumption, the common approach is prescribing a time-dependent pressure waveform instead [56,109]. However, invasively measuring pressure for branching arteries is challenging, so a general pressure waveform pattern demonstrated in Figure 5 is typically applied to all outlets. This method is frequently utilized in AAA simulations with an inlet at the IR region, where at most two outlets, LCI and RCI, exist [8,74,76,110,175,176,177,178,179,180,181,182,183]. However, this method is less suitable for AAA simulations with an inlet at the SC. In that case, the increased number of outlets from branching arteries between SC and IR segments brings complexities because the diameters and flow distribution among those arteries differs and the time-varying pressure data cannot be obtained simultaneously for each outlet. Reymond et al. [126,184] have obtained pressure waveforms for all outlet locations at ascending aorta using a 1D model and validated with in-vivo measurements to provide relevant physiological data. Furthermore, this method neglects the effect of arterial compliance, neglecting the downstream wave propagation. For a realistic simulation, matching the inlet flow waveform and outlet pressure distribution is essential and should be taken into consideration.

3.2.2. Flow-Split Method

The flow split method assigns specific flow rates to each outlet. These rates are determined through either formula-based calculations or direct in vivo measurements. Les et al. [37] reported the constant fractions of mean flows to each outlet. As shown in Figure 6a., the flow between SC and IR regions, called upper branch vessel flows (UBVF), is determined by subtracting the mean measured IR flow from the mean measured SC flow, corresponding to 1.31 and 3.51 L/min, respectively [37,64,65,185]. The remaining 2.2 L/min is distributed to upper vessel branches: 33% to the celiac trunk (CT), and 22.3% to superior mesenteric (SMA), the left renal (LR) and right renal (RR) arteries [37]. In the presence of accessory renal arteries, renal flow was divided proportional to the outlet area. CT is then branched into hepatic (H) and splenic (S) arteries, with the flow distributed equally between them [37]. The remaining IR flow is divided equally into the two common iliac arteries, LCI and RCI, and 70% of this flow is diverted to external iliac arteries, while 30% is sent to internal iliac arteries [37,72,133,146]. These percentages are tabulated in Table 2.

To calculate the flow rate split among the outlets, Murray’s law is used [186,187]. The application of Murray’s law to estimate the flow splitting at artery bifurcations has been investigated in the literature. Murray [188] formulated that in branching arteries, the flow in each outlet is proportional to the cross-sectional area of the bifurcated vessel. The general form of the Murrays law is given in Eq. (3.7).

$${\displaystyle \frac{Q^i}{Q_{tot}}}={\displaystyle \frac{r_i^n}{\sum _{i=1}^{N_{outlet}}{r}_i^n}}$$

(3.7)

In Eq. (3.7), $Q^i$is the flow rate at the $i^{th}$ outlet, $Q_{tot}$ is the total flow rate, $r_i$ is the radius of $i^{th}$ outlet, $N_{outlet}$ is the number of outlets. The exponent $n$ varies according to the arterial segment. Generally, it is considered to be 2 for the aortic segment [130,189]. Several studies have used an exponent $n$ of 3 for AAA [45,190,191,192]. However, specifying a fixed flow rate at each outlet is not a realistic boundary condition because the flow division changes during the cardiac cycle, especially in deformable wall simulations where the outlet area changes continuously [193]. Furthermore, this BC neglects wave transmission to downstream vasculature [160].

3.2.3. Lumped Parameter Model

To accurately model downstream vasculature at the outlets of the AAA, coupling the 3D computational domain with reduced order (0D or 1D) is a common practice [135,194]. 0D models, also called lumped models, are governed by a group of ordinary differential equations (ODEs) that assume spatially uniform distribution of blood pressure, $P\left(t\right)$, and flow rate, $Q\left(t\right),$ within cardiovascular compartments at any time instant [195]. In vascular systems, the most widely used 0D models are Windkessel (WK) models [161,196]. This method represents downstream vasculature using lumped-parameter networks similar to electrical circuits consisting of capacitors, resistors, and inductors to relate $Q\left(t\right)$ and $P\left(t\right)$ [189]. The elasticity of arterial walls allows them to expand and store large volumes of blood, then relax and push that blood downstream. This phenomenon, referred to as arterial compliance ($C$), is analogous to the behavior of a capacitor in an electrical circuit. Similarly, as blood flows downstream, the arteries bifurcate and form small-diameter capillaries and veins. The narrowing of arterial diameters generates significant peripheral resistance ($R)$, analogous to electrical resistance [197]. In larger arteries, inertial energy is stored and released due to acceleration or deceleration of fluid, which is expressed as an inertance ($L)$, analogous to electrical inductor. In this framework, blood flow and pressure are analogous to current and voltage, respectively. By applying Kirchoff’s voltage and current laws, ODEs governing the relationships between the $P\left(t\right)$ and $Q\left(t\right)$ are derived, called as WK models [198,199].

WK models are classified based on the number of parameters or circuit elements ($R,$ $C$ and $L)$ [199]. Figure 8 illustrates the most common types: the two-element, three-element, and four-element WK models [198]. In two-element WK model (WK2 or RC), $P\left(t\right)$, is related to $Q\left(t\right)$ and the distal pressure, $P_d\left(t\right)$, through a parallel combination of distal resistance ($R_d)$ and a compliance ($C)$. According to Poiseuille’s law, the resistance is inversely proportional to the fourth power of the artery radius, which makes the smallest arteries and arterioles the primary contributors to resistance in the cardiovascular system [199]. Consequently, the WK2 model accounts only for $R_d$, which represents the resistance of the downstream vasculature caused by the small arteries and the capillary bed [79,200,201]. However, the WK2 model has limitations in accurately describing the pressure-flow rate relationship because it omits the resistance of larger arteries [198,199]. To address this, three- and four-element WK (WK3 and WK4) models were developed [44,198,202]. The WK3 model sufficiently replicates the realistic downstream pressure at the outlet sections that are consistent with experimental data. As a result, the WK3 model is commonly implemented in most hemodynamic simulations [142,155,198,203].

The WK3 model, often called the RCR model, consists of a proximal resistance ($R_p$) in series with a parallel arrangement of a distal resistance ($R_d$) and a compliance ($C$). $R_p$ describes the resistance in the large ascending aorta, proximal to the AAA. Based on the wave transmission theory [198,204] it helps absorb incoming pressure waves and reduce artificial wave reflections [135]. In the literature, the proximal and distal resistances are sometimes denoted as $R_1$ and $R_2$, respectively [133,135,205]. The sum of distal and proximal resistances is known as total arterial resistance, $R_T=R_p+R_d$ [37,197]. The WK3 model describes the relationship between $P\left(t\right)$ and $Q\left(t\right)$ using the ODE in Eq. (3.8).

$$P\left(t\right)=R_{T}Q\left(t\right)-R_{d}C{\displaystyle \frac{dP\left(t\right)}{dt}}+R_p{R}_{d}C{\displaystyle \frac{dQ\left(t\right)}{dt}}+P_d\left(t\right)+R_{d}C{\displaystyle \frac{d{P}_d\left(t\right)}{dt}}$$

(3.8)

In Eq. (3.8), $P_d\left(t\right)$ represents the pressure at which blood flow to smaller arteries and capillaries in the vascular bed ceases [135,206] and is typically assumed to be zero in the literature [130,133,189]. For each outlet $i$, Eq. (3.8) can be written as

$$P^i\left(t\right)=R_T^i{Q}^i\left(t\right)-R_d^i{C}^i{\displaystyle \frac{dP\left(t\right)}{dt}}+R_p^i{R}_d^i{C}^i{\displaystyle \frac{d{Q}^i\left(t\right)}{dt}}i=1,2,3,\dots ,N_{outlet}$$

(3.9)

Figure 9 illustrates the 3D domain of AAA and arterial branches coupled with separate WK3 models at the outlets: S, H, SMA, LR, RR, LEI, LII, REI, and RII. This technique models downstream compartments of the arterial branch outlets in the 3D AAA domain. The parameters $(R_T^i, R_d^i, R_p^i$ and $C^i)$ are constant for each outlet and independent of spatial dimensions [189]. By using $Q^i\left(t\right)$ calculated in the 3D simulation and estimating WK3 parameters for each outlet $i$, $P^i\left(t\right)$ can be determined by solving Eq. (3.9). Various techniques for estimating WK3 parameters exist in the literature. When patient-specific $P^i\left(t\right)$ and $Q^i\left(t\right)$ are available for each outlet, the least-square approach recently proposed by Romarowski et al. [207] is recommended. This method tunes the WK3 parameters to match measured in vivo patient-specific pressure and flow data [130]. However, both $P^i\left(t\right)$ and $Q^i\left(t\right)$ are not always simultaneously available, as most branches are quite narrow to make in vivo measurements.

Some studies have taken the parameters from the literature data to avoid time-consuming parameter estimation process [155,208,209]. However, the parameters utilized should be carefully selected to match accurate patient pressure profiles because the values of WK3 parameters affect the hemodynamic parameters. In an uncertainty estimation study, Boccadifuoco et al. [189] reported that the $P\left(t\right)$waveform at the aorta outlet is mainly affected by capacitance, $C$, while the $P_{max}$is affected by $R_p$. Consequently, iterative solutions or approximations are frequently utilized to tune WK3 parameters [210]. Several studies implemented formulations of 1D modeling for parameter estimation [32,37,72,135,137,197,204,211]. A good agreement has been reported between the simulated pressure and flow rates and the in vivo measurements [37,135,195,196].

In this technique, a global WK3 model is created as the 0D representation of the 3D computational domain, as illustrated in Figure 9. The global WK3 parameters ($R_T^{SC}, R_d^{SC}, R_p^{SC}$ and $C^{SC})$ are defined at the inlet of the computational domain before calculating local WK3 parameters $(R_T^i, R_d^i, R_p^i$ and $C^i)$ for each outlet $i$. An initial estimate of the global compliance ($R_T^{SC})$can be calculated by using the Eqs. (3.10) and (3.11) defined for total resistance in 1D modeling [130,133,135]

$$R_T^{SC}={\displaystyle \frac{P_{mean}^{SC}-P_{out}}{{\overline{Q}}^{SC}}}$$

(3.10)

$$P_{mean}^{SC}={\displaystyle \frac{1}{3}}{P}_{sys}^{SC}+{{\displaystyle \frac{2}{3}}P}_{dia}^{SC}$$

(3.11)

In Eqs. (3.10) and (3.11) $P_{mean}^{SC},P_{sys}^{SC}\mathrm{a}\mathrm{n}\mathrm{d} P_{dia}^{SC}$ are the mean, systolic and diastolic pressures, and ${\overline{Q}}^{SC}$ is the time-averaged mean flow rate at the inlet [37,156,189,205,210]. The outlet pressure ($P_{out})$ is taken as 4.4 kPa in some studies [135], while generally is considered as $P_{out}=0$ [39,45] $.$ In 1D simulations, the global compliance ($C^{SC})$, which is the ratio of a volume change, $∆V,$ and the resulting pressure change, $∆P$, is approximated as in Eq. (3.12).

$$C^{SC}={\displaystyle \frac{∆V}{∆P}}={\displaystyle \frac{Q_{max}^{SC}-Q_{min}^{SC}}{P_{sys}^{SC}-P_{dia}^{SC}}}∆t$$

(3.12)

In Eq. (3.12) $Q_{max}^{SC}\mathrm{a}\mathrm{n}\mathrm{d} Q_{min}^{SC}$ are the maximum and minimum flow rates at the inlet, and $∆t$ is the difference between the maximum and minimum flow rates [135,196,211]. Les et al. [37] determined the system's global resistance and capacitance by iteratively solving for a pressure waveform instead of using Eqs. (3.10-12). This waveform was derived and further iterated as a function of the input SC flow waveform and initial estimation of Windkessel model parameters.

The global parameters ($R_T^{SC}$ and $C^{SC}$) are used to calculate the local WK3 parameters at each outlet segment by utilizing Murray’s law [37,45,156,171,189,191]. Consequently, the total resistance and compliance at $i^{th}$ outlet can be calculated by using Eqs. (3.13) and (3.14), respectively.

$$R_T^i={\displaystyle \frac{\sum _{i=1}^{N_{outlet}}{A}_i}{A_i}}{R}_T={\displaystyle \frac{A_{tot}}{A_i}}{R}_T^{SC}={\displaystyle \frac{{\overline{Q}}^{SC}}{Q^i}}{R}_T^{SC}$$

(3.13)

$$C^i={\displaystyle \frac{A_i}{\sum _{i=1}^{N_{outlet}}{A}_i}}{C}^{SC}={\displaystyle \frac{A_i}{A_{tot}}}{C}^{SC}={\displaystyle \frac{Q_i}{{\overline{Q}}^{SC}}}{C}^{SC}$$

(3.14)

In Eqs. (3.13) and (3.14), $N_{outlet}$ is the number of outlets, $A_i$ and $A_{tot}$ are the area of $i^{th}$ outlet and sum of the area of all the aortic outlets, respectively. Table 2 presents the ratios of time-averaged mean flow rates diverted to sub-branches in the AAA section after the SC region. These values, reported by Les et al. [37], are widely used in WK3 parameter tuning for AAA simulations. The ratio between proximal and total resistance, $R_p^i/R_T^i$ is taken as 5.6% for most of artery segments, with the exception of renal arteries [5,22,37,39,45,133,197,212,213]. Renal arteries, connected to kidneys, have relatively low $R_d^i$ at the vascular beds, resulting in an $R_p^i/R_T^i$ ratio of $28\%$ [133]. Once the total resistance at the $i^{th}$ outlet is determined, the proximal and distal resistances can be calculated as

$$R_T^i=R_p^i+R_d^i;R_p^i=0.056R_T^i;R_d^i=0.944R_T^i$$

(3.15)

The calculation of local from the global WK3 parameters for the downstream compartments of CT, SMA, LR, RR, LCI, RCI, is summarized in Table B in Appendix B. The flowchart in Figure 10 summarizes overall process. Although these 1D formulations are frequently used in AAA simulations, still concerns exist regarding their applicability to detailed geometrical shapes and FSI simulations where the flow rate directed to each branch is not constant due to the arterial deformations [156].

To improve parameter estimation accuracy, various tuning methods have been developed, primarily for ascending aorta and aortic dissections. Jiang et al. [214,215] initially performed a CFD simulation with zero-gauge pressure imposed at all outlets, calculating flow rates and pressures to tune WK3 parameters. In subsequent rounds, they imposed WK3 with estimated values, continuing until the pressure difference was less than 5%. Several studies estimate local parameters through iterative ODE solutions to match physiological pressure distribution within specified systolic and diastolic limits for the AAA [39,212,216,217]. Some of the studies use general values like 120/80 mmHg [38,156,189,201,218,219]. Local WK3 parameters are adjusted until the desired waveform is achieved. Fonken et al. [39,213] measured patient-specific brachial blood pressure in supine position using a brachial cuff. Alimohammadi et al. [205] measured minimum and maximum pressures at all branches of type-B aortic dissection using a transfemoral flush angiographic catheter connected to a pressure transducer. To further improve estimation accuracy, Alimohammadi et al. [205] and Pant et al. [220] applied data assimilation technique and unscented Kalman filters, respectively.

Spilker et al. [221] and Bonfanti et al. [133] used reduced model tuning techniques to enhance initial guess accuracy, and various CFD simulation results have been used for further iterations to optimize local parameter fine-tuning. Li et al. [156] proposed a fast-estimating approach that eliminates the need for CFD simulation iterations. They initially obtain total resistances as input for each interface using 1D formulations, considering systolic and diastolic pressures as 120/80 mmHg. The local WK3 parameters are then optimized using the pattern search algorithm from MATLAB's global optimization toolbox.

After estimating the parameters, coupling of 3D with 0D domain is commonly done explicitly [222,223] or semi-implicitly. The derivative terms in Eq. (3.9) are typically discretized as in Eq. (3.16) by using backward Euler method [205]

$$P_{n+1}^i={\displaystyle \frac{\left(R_p^i+R_d^i+R_p^i\beta \right)Q_{n+1}^i-R_p^i\beta Q_n^i+\beta P_n^i}{1+\beta }}$$

(3.16)

In Eq. (3.16), $\beta =R_d^{i}C/ ∆t$ [224]. To achieve appropriate coupling, the flow rate equality condition for the multiscale model at the interface, $Q_{3D}^i=Q_{0D}^i$, must be satisfied [205]. The instantaneous pressure, $P_{0D}^i\left(t\right)$, is calculated from Eq. (3.16) by using the flow rate for the current solver loop, $Q_{3D}^i\left(t\right)$ and pressure and flow rate from previous time step, $P_{0D}^i(t-∆t)$ and $Q_{3D}^i(t-∆t)$, respectively [205,223]. The resulting pressure, $P_{0D}^i\left(t\right)$ is then supplied back to the 3D domain as the uniform pressure BC at the current time step. For implicit coupling, $P_{0D}^i\left(t\right)$ and $Q_{3D}^i\left(t\right)$ is solved simultaneously in an iteration loop [225]. In explicit coupling, pressure is calculated at the beginning of an iteration by using the flow rate value from the previous time step.

WK3 coupling isn’t available on every solver. SimVascular offers coupled WK2 and WK3 [203], while it is not defined in other software. To make such a multidomain solution in OpenFOAM, heamofoam [226] can be used. To use WK3 BC in ANSYS, a user defined function (UDF) should be created to couple two domains. For CFX, Fortran subroutine can be defined [205,218]. In some studies, 3D-0D coupling is not performed. In several studies [42,227,228], time-dependent pressure waveform has been calculated through WK3 and applied them as prescribed pressure BC at the outlets.

3.2.4. Resistance BC

To relate the flow rate and the pressure at the outlets of the computational domain, some studies prefer to use resistance BC due to its simplicity [24,78,190,229]. In this technique, the resistive property of the downstream vasculature is imposed by a resistance parameter, while the compliance of the arteries is neglected with Eq. (3.17).

$$P_{out}^i\left(t\right)=P_{ref}+R_{out}^i{Q}_{out}^i\left(t\right)i=1,2,3,\dots ,n_{outlet}$$

(3.17)

In Eq. (3.17), $R_{out}^i,$ $P_{out}^i\left(t\right)$ and $Q_{out}^i\left(t\right)$ are the resistance, time-dependent pressure and flow rate at the $i^{th}$ outlet, $P_{ref}$ is the reference pressure related to the venous pressure. $R_{out}^i$ is calculated by distributing the global resistance calculated from Eq. (3.10) into branch arteries with Eq. (3.13). However, resistance BC is not commonly preferred BC in AAA simulations since it severely impacts wave propagation phenomena [160,230].

4. Modeling Boundary Conditions in Solid Domain

4.1. Inlet and Outlet of the Wall

The branch arteries of aorta produce a tethering effect on the AAA. To model this effect, zero rotation and translation conditions are imposed at the inlet and outlet of the wall domain [10]. To further improve the accuracy, Scotti et al. [69] specified 5% axial stretch on these boundaries [231] since the artery walls are physiologically under tension.

4.2. External Wall Boundary

The intra-abdominal pressure imposed by surrounding tissues and organs on the outer surface of the AAA wall requires consideration [232]. Most studies, however, assume a free stress condition with zero pressure on the external AAA wall. This boundary condition can lead to non-physiological wall motion patterns [193], causing several researchers to use a constant intra-abdominal pressure of 12 mmHg [69]. Crosetto et al. [233] introduced a Robin condition which requires appropriate model parameters through curve fitting to the empirical data. Similarly, Moireau et al. [234,235] developed a BC along the outer wall of thoracic aorta (TA) incorporating a viscoelastic term to account for surrounding tissue and organ support.

4.3. FSI Boundary: Coupling Solid and Fluid Domains

FSI modeling can be performed using several numerical approaches: one-way uncoupled [13], one-way coupled [66,236], and two-way coupled, both explicit and implicit [113]. In the one-way uncoupled method, intraluminal pressure load is supplied to the FSI interface, and only solid mechanics equations are solved in the wall domain, without feeding fluid domain with the pressures obtained by solid domain. In two-way coupled methods, solid and fluid domains are simulated simultaneously, interacting through the FSI interface. Comparative studies [236] showed that one-way FSI technique overestimates pressure, WSS, and strain in the fluid domain while underestimating the von Mises stress and displacement. As shown in Eqs. (4.1) -(4.3), three conditions must be met at FSI boundaries: solid and fluid domain displacements must be compatible, boundary tractions must be in equilibrium, and fluid must follow the no-slip condition on the FSI boundary surface.

$${\mathit{d}}_{\mathit{s}}={\mathit{d}}_{\mathit{f}}$$

(4.1)

$${\mathit{\tau }}_{\mathit{s}}·{\widehat{\mathit{n}}}_{\mathit{s}}={\mathit{\tau }}_{\mathit{f}}·{\widehat{\mathit{n}}}_{\mathit{f}}$$

(4.2)

$$\mathit{v}=\dot{{\mathit{d}}_{\mathit{f}}}$$

(4.3)

In Eqs. (4.1)-(4.3), $\mathit{d},\mathit{\tau },\widehat{\mathit{n}}$ and ${\mathit{d}}_{\mathit{f}}$ are displacement vectors, stress tensors and boundary normal vectors at the FSI boundary for solid, s, and fluid, f, domains. In two-way explicit coupling, information exchange among fluid and solid domains occurs explicitly, without iterations within each time step. To obtain converged results, small steps are required. In two-way implicit coupling, the fluid and solid equations are solved simultaneously at each time step, iteratively exchanging information until convergence. While this method requires high computational memory and excessive computational time, it is necessary to reach convergence for improving accuracy.

5. Important Post-Processing Indices

Hemodynamics inside an aneurysm sac can be quantified by several wall shear stress (WSS) descriptors such as time-averaged wall shear stress (TAWSS), oscillatory shear index (OSI), endothelial cell activation potential (ECAP) and relative residence time (RRT). Mathematical definitions of these descriptors are given in Eqs. (5.1)-(5.4) [237,238].

$$\mathrm{T}\mathrm{A}\mathrm{W}\mathrm{S}\mathrm{S}={\displaystyle \frac{1}{T}}{\int }_0^T\left|{\tau }_w\right|dt$$

(5.1)

$$\mathrm{O}\mathrm{S}\mathrm{I}=0.5\left(1-{\displaystyle \frac{\left|\frac{1}{T}{\int }_0^T{\tau }_{w}dt\right|}{\frac{1}{T}{\int }_0^T\left|{\tau }_w\right|dt}}\right)$$

(5.2)

$$\mathrm{E}\mathrm{C}\mathrm{A}\mathrm{P}={\displaystyle \frac{\mathrm{O}\mathrm{S}\mathrm{I}}{\mathrm{T}\mathrm{A}\mathrm{W}\mathrm{S}\mathrm{S}}}$$

(5.3)

$$\mathrm{R}\mathrm{R}\mathrm{T}={\displaystyle \frac{1}{(1-2·\mathrm{O}\mathrm{S}\mathrm{I})·\mathrm{T}\mathrm{A}\mathrm{W}\mathrm{S}\mathrm{S}}}$$

(5.4)

where $T$ and ${\tau }_w$ are the cardiac cycle period and the wall shear stress, respectively.

Several studies also reported a correlation between recirculation zones and WSS descriptors. The utilization of vortex fields may afford a more comprehensive understanding on ILT development and rupture mechanism, compared to WSS descriptors [239,240]. In the literature, vortex structures are commonly quantified by $Q$, $∆$, ${\lambda }_2$, ${\lambda }_{ci}$ (swirling strength) criteria [241,242] and rortex [243]. The equations for vortex identification criteria are given in Appendix C.

For the AAA wall, principal wall stresses, Von Mises stress and wall displacements are critical parameters. Von Mises stress, based on three principal stresses, is a measure of failure prediction as given in Eq. (5.5).

$${\displaystyle \frac{1}{2}}\left[{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2\right]>{\sigma }_y^2$$

(5.5)

In Eq. (5.5), ${\sigma }_1, {\sigma }_2, {\sigma }_3$ are the principal stresses, ${\sigma }_y$ is the uniaxial failure strength of the wall [34]. The term on the left-hand side is the square of the Von Mises stress. Peak wall stress (PWS) is a critical indicator of rupture, occurring when the AAA wall strength is insufficient to withstand the PWS [244].

6. Recent Findings

To accurately simulate the AAA hemodynamics and wall mechanics, researchers have developed and implemented different techniques over a quarter century. Starting from the axisymmetric geometries supplied with simplified boundary conditions and material properties, in silico investigation of the AAA development and rupture characteristics has reached patient-specific geometries with realistic boundary conditions, enhancing the reliability of predictions. Intense research has accomplished to couple different computational domains solved by various methods, such as CFD and FEA at the wall interface by implementing FSI technique, and different dimensions such as 0D and 3D. Increasing use of medical imaging techniques such as CT, MRI and 3D ultrasound enables the extraction of patient-specific geometries with accurate wall thickness distribution with eliminated prestress that ensures precisely capturing the anatomical features of AAA. Utilization of the 4D flow MRI and angiographic catheters enable to provide time-dependent patient-specific inlet velocity profiles and outlet pressure waveforms, respectively. With the help of biaxial and uniaxial tissue tests, accurate constitutive models with patient-specific material constants have been developed.

The utilization of 4D flow MRI inlet velocity profiles is common in the thoracic aorta (TA) studies [81,130,138,142,217], while the implementation of artificial velocity profiles developed from patient-specific inflow waveforms at SC or IR regions is standard practice in AAA simulations. Chandra et al. [13] obtained patient-specific velocity profiles at IR region using 4D flow PC-MRI technique and compared the results with simulations using artificial profiles. Their results showed that WSS, PWS, and strain are influenced by the hemodynamics from the applied inflow boundary condition. A comprehensive biomechanical approach for assessing AAA rupture risk should consider the interaction between the aortic wall and hemodynamics with patient-specific inflow boundary conditions [81,148]. Wei et al. [150] found no significant difference between PC-MRI measured and artificial inlet profiles, though the flat profile showed notable differences in Fontan hemodynamics. Ramazanli et al. [21] compared artificial velocity profiles for AAA hemodynamics and recommended the parabolic velocity profile for its simplicity. However, the flat profile is not recommended by many studies [5].

Prescribing pressure waveform at the outlet of branching arteries remains common in AAA simulations, particularly for 3D domains starting after the IR section with two outlets. Wang et al. [5] compared various inlet and outlet BCs in an AAA domain with two outlets, using a two-way FSI model with hyperelastic and anisotropic properties. They prescribed parabolic and flat profiles at the inlet, along with pressure waveform from the literature and WK3 BCs at the outlet. Their comparison of average TAWSS along AAA, PWS, and maximum displacement showed that WK3 significantly affected all three measurements. The prescribed pressure BC overestimated PWS and displacement while underestimating TAWSS.

Pirola et al. [130] applied a different approach, using patient-specific pressure waveforms measured at each thoracic aorta outlet before MRI geometry extraction. They compared patient-specific outlet pressure BC with WK3 and zero pressure BCs, finding that outlet BC choice significantly impacts aorta hemodynamics. Their results showed that WK3 can accurately reproduce physiological aortic pressure waveforms, suggesting WK3 model when patient-specific pressures are unavailable. Madhavan and Kemmerling's [155] comparison of WK2, WK3, and flow splitting BCs revealed minimal variations between WK2 and WK3, though they found an 18% difference in TAWSS between WK models and flow splitting. Boccadifuoco et al. [189] studied uncertainties in WK3 parameters for both rigid and deformable wall models. Their stochastic analysis showed that compliance $\left(C\right)$ has the greatest impact on hemodynamic predictions compared to distal $\left(R_d\right)$and proximal $\left(R_p\right)$ resistances. However, they concluded that WK3 parameters may not be a major uncertainty source, particularly when considering wall elasticity.

In the literature, the rigid wall assumption is typically applied in a hemodynamic investigation of AAA, resulting in underestimated PWS and overestimated velocity, WSS and OSI values. Integrating FSI at the AAA wall interface increases the simulation accuracy by accounting for wall deformation. To further enhance reliability, researchers have developed constitutive models with patient-specific material constants for wall layers and ILT structure. Comparative studies of the elastic and hyperelastic models show that elastic models underestimate the von Mises stress on the AAA wall [77,106,245] and predict the wall displacement about 20% greater than hyperelastic material [106]. Xenos et al. [124] compared the isotropic Mooney-Rivlin and anisotropic Holzapfel models, finding that the isotropic model underestimates PWS. Similarly, Rissland et al. [246] discovered that isotropic models underestimate the wall deformation. Balzani et al. [111] compared Neo-Hooken isotropic model with anisotropic models, concluding that isotropic models oversimplify arterial mechanical behavior. While anisotropic models are preferred, they significantly increase numerical complexity. To address this, Wang et al. [5] developed a combined strain energy function that models extracellular matrix energy as isotropic and collagen fiber effects as anisotropic.

Though AAA walls vary in thickness, most studies assume constant thickness due to geometry extraction challenges [108,247]. However, several studies show that wall thickness significantly affects wall biomechanical properties [101,248]. Raghavan et al. [249] observed that the minimum wall thickness could reach 0.23 mm around the rupture site. Scotti et al. [10] observed that varying wall thickness leads to an increase of the von Mises stress up to 4 times compared to the uniform thickness. Raut et al. [250] compared maximum principal stress and strain, strain energy density, and displacement magnitude across different thickness models. Their findings strongly support using patient-specific, regionally varying wall thickness from CT scan segmentation, particularly for FEA analysis of AAA.

Addressing pre-stress in CT-reconstructed geometry improves the accuracy of the results. The reconstructed AAA geometry reflects a pressurized aortic state due to the physiological pressure during image acquisition, potentially leading to WSS distribution underestimation. Fonken et al. [39] implemented a backward incremental method (BIM) to estimate pre-stress in measured geometry. Omitting pre-stress leads to increased systolic displacements and decreased systolic wall stresses up to 77.8% and 54.2%, respectively.

In the literature, a vast number of studies have postulated a correlation between the wall shear stress parameters and intraluminal thrombus (ILT) formation [9,37,146,251]. The most well-known approach is that the low and oscillatory WSS may stimulate the endothelial cells and promote the inflammatory process, causing wall-cell adhesion due to platelet accumulation, forming ILT. WSS values lower than 1 Pa are considered as low WSS [41,174]. Most researchers agree that the locations with low WSS, high OSI and high ECAP are prone to thrombus formation and have a higher risk of rupture [238,251]. On the other hand, in some studies reporting contradictory results [1,11,12,13], it is stated that low WSS and high OSI regions do not necessarily coincide with thrombus deposition and atherosclerosis formation [76,146,180].

Several studies in the literature [239,252] remarked the limitations of WSS descriptors as scalar-tensor fields lacking directional information. These descriptors such as OSI, ECAP and RRT can only indicate potential pathology but not possible progress, nor its underlying mechanisms [252]. Saqr et al. [239] highlights the utilization of vector fields, such as vorticity, to visualize hemodynamics. The utilization of vortex fields may afford a more comprehensive understanding on ILT development and rupture mechanism, compared to WSS descriptors [239]. While the effect of wall behavior on WSS parameters has been extensively studied, there is limited research on wall deformation and vortex interactions within the AAA lumen. Kelsey et al. [251] stated that low-velocity recirculation zones are located near the regions where ILT has formed. Varble et al. [253] observed a significant correlation between intracranial aneurysm (IA) rupture and the near-wall vortical patterns exhibiting elevated vorticity levels. Biasetti et al. [76,180] correlated vortical structures with high WSS and attributed bursting of the vortical structures to thrombus deposition at low WSS areas. Several studies have examined the relationship between WSS parameters and vortex structures [240,254]. Zhan et al. [255] utilized the rortex as the vortex identification technique and found a positive correlation between vortices and high WSS regions. However, their study was limited by using zero pressure at the outlet and rigid wall boundary conditions.

7. Potential for Integrating AI and ML in AAA Research

Recent developments in artificial intelligence (AI) have influenced medicine. For that reason, the use of machine learning (ML) and data-driven approaches in AAA diagnosis and treatment is becoming increasingly common. By reducing manual interventions in modeling and analysis processes, ML methods can eliminate the risk of human error. Reducing the time required for analyses that may take weeks or days to order of hours or even real time by ensuring faster execution of risk assessment processes is one of the advantages of data-driven methods.

ML describes the ability of systems to automatically learn from problem-specific data, providing a powerful alternative to manually building analytical models and enabling the automation of complex tasks. In recent years, the use of ML has rapidly expanded beyond computer science including medical simulations [256] and ML has made great strides in sophisticated learning algorithms and efficient pre-processing techniques such as the naissance of deep learning (DL) which is shortly the evolution of artificial neural networks (ANNs) towards deeper neural network architectures and improvement of their learning capabilities [257]. DL methods are based on representation learning, transforming raw data into abstract, and complex functions [257].

DL methods are used to improve computationally intensive methods such as FEA and CFD, making biomechanical simulations faster and more efficient. It is demonstrated [258] that a trained DL model could predict stress distributions with average errors of 0.492% in the Von Mises stress distribution and 0.891% in peak Von Mises stress, significantly reducing the long computing times typically required by FEA. This approach draws attention to DL as an alternative for FEA and CFD in biomechanical and stress analyses, rapidly advancing patient-specific modeling in time-sensitive clinical applications. Similarly, in AAA rupture risk assessment, generalized additive model has shown high accuracy in analyzing geometric and biomechanical markers for the risk prediction [14].

The paradigm of ML and DL is the development of data-driven algorithms. Either structured or unstructured data are utilized to collect and derive the necessary task-related information. DL models have an advantage over traditional ML models due to their increased number of learning layers and higher abstraction levels on large data sets and complex insights [259]. Consequently, data-driven methods, particularly those utilizing deep learning architectures, are increasingly used in AAA diagnosis and treatment studies.

Deep architecture, composed of many non-linear transformations, can compactly represent highly complex functions [260]. Deep neural networks (DNNs) are multi-layered artificial neural networks (ANNs) that can represent complex functions and learn complicated patterns more efficiently with fewer parameters compared to shallow networks. Convolutional neural networks (CNNs) are analogous to traditional ANNs by consisting of neurons that optimize themselves through learning, but they are particularly effective in pattern recognition tasks involving images [261]. CNNs reduce the need for fully connected layers by focusing on learning in convolutional layers, which decreases computational cost and the number of parameters [262].

CNNs have gained significant popularity in computer vision due to their ability to efficiently process spatial information compared to fully connected architectures and have become a preferred method for medical imaging tasks [263,264]. Deep convolutional neural networks (DCNNs) concentrate on the learning in convolutional layers, reducing the need for fully connected layers, which lowers the number of parameters and computational cost [262]. These networks use local connections to identify image features like edges, corners, and textures, enabling better interpretation of spatial information [265].

Segmentation processes for AAA simulations, whether manual or semi-automatic, have benefited from advancements in image processing algorithms, large datasets, and DL techniques. DL based segmentation methods such as ARU-Net [51] and CACU-Net [52] have reduced the time to generate computational models from computed tomography angiography (CTA) scans from approximately 2 hours to 10 minutes, showing reliable performance in creating patient-specific geometries for CFD simulations [266]. Figure 11 illustrates fully automatic segmentation of AAA from CT images by adapting the resnet-based fully convolutional networks (FCN) [267]. The model consists of three steps, extracting the aorta, iliac arteries, and detecting the lumen and other AAA tissues. The automated segmentation results demonstrate a goof agreement with manual segmentation. DCNN achieved a dice similarity of 82% in segmenting intraluminal thrombus (ILT), showing that AI based segmentation has potential to be used for clinical applications, especially in cases with difficult-to-detect boundaries [268].

Applications of AI and ML in AAA Flow Simulations

The application of AI and ML in AAA simulations has gained considerable attention, as several studies have demonstrated their potential in improving both speed and accuracy in predicting hemodynamic parameters and disease progression. For instance, Liang et al. [269,270] showed that DNNs, trained on hemodynamic data derived from CFD simulations, could predict steady-state distributions of pressure, velocity, and velocity magnitude within one second. This achievement highlights the capacity of DNNs to act as rapid surrogates for traditional CFD simulations, significantly reducing computational time. Rościszewski et al. [271] also explored the integration of AI to accelerate CFD simulations, emphasizing the importance of such techniques in clinical settings where fast and accurate predictions are critical. As demonstrated in studies such as Feiger et al. [272] and Hahn et al. [273], at least nine simulations were necessary to train neural networks to predict time-averaged WSS and pressure gradients across stenotic regions, illustrating how ML can be combined with CFD and FEA simulations to develop precise and patient-specific flow prediction models [53]. Similarly, Kim et al. [274] applied neural network architectures to a dataset of 54 patients, testing the model on four key features: radius, ILT thickness, time-averaged WSS, and aneurysm growth rate.

Recurrent neural network (RNN) and CNN based approaches have also been widely used for the prediction of future events in medical applications [275,276]. CFD simulations conducted on patient-specific 3D geometries show that DL frameworks, trained with multi-physical features to overcome small amount of longitudinal AAA data, can outperform traditional methods in predicting AAA growth with improving accuracy and efficiency. These improvements can help clinicians manage AAA progression and improve patient care [274]. Jiang et al. [277] employed a combination of physical vascular adaptation modeling, machine learning (ML) tools, and follow-up scan data to predict the shape evolution of AAAs using a Deep Belief Network (DBN).

Soudah et al. [278] demonstrated a multilayer perceptron (MLP) based approach, creating two neural networks for AAA analysis. The first, a Mesh Neural Network (MNN), generates the aneurysm geometry based on four geometric factors, while the second, the Tension Neural Network, calculates maximum wall stress by combining the MNN output with arterial pressure, achieving an accuracy rate of 95% compared to the finite element method (FEM) based results. However, as MNN results can be noisy due to uneven node distribution, further studies are necessary to account for factors such as ILT effects. Additionally, statistical shape models (SSMs) have been used to describe the shape variability of the aortic arch, with nonlinear regression employed to analyze pressure gradients as a function of flow rate and cross-sectional area [208,273,279].

Joly et al. [208] utilized longitudinal cohort data from CFD simulations and CT scans to study hemodynamic correlations, highlighting ML's potential to handle complex biomechanical phenomena. However, the unpredictable nature of fluid dynamics can negatively affect the AI predictions, which is why combining AI-driven mechanisms with FEM could improve the accuracy of models. For example, Jiang et al. [280] adopted this hybrid approach to estimate aneurysm growth in the aorta, achieving greater accuracy while reducing computation time through optimization techniques and AI.

Recent advancements in AI have also enhanced MRI-based blood flow measurements, particularly with 4D flow MRI, as ML models like CNNs and U-Nets have been applied to automate tasks such as phase-contrast imaging, vessel segmentation, and contour drawing [281]. ML-based super-resolution techniques have been developed to integrate high-resolution CFD simulations with MRI data, resulting in more realistic flow simulations without sacrificing data fidelity [282].

Beside these improvements and advantages, there are various challenges to overcome such as choosing the right implementation, bias and drifting in data, and the mitigation of black-box properties [283]. Further research is necessary to refine AI models for broader clinical applications. For example, CNNs used for segmentation still require significant improvements before they can be widely adopted in clinical settings [266]. Additionally, future research could focus on retraining DNNs using datasets generated by more advanced simulation techniques, such as FSI, to enhance model accuracy and clinical relevance. Exploring larger datasets and more complex architectures could further enhance the predictive capabilities of AI models. Given the scalability and rapid computation time of DNNs, they hold great promise for future clinical applications [269].

8. Conclusions

In this paper, we focused on modeling realistic BCs and application of ML techniques for AAAs, providing a comprehensive overview and summarizing their implications. The main conclusions of the study are provided below.

The numerical models developed to simulate a specific region of the cardiovascular system should accurately reproduce patient-specific hemodynamics and biomechanics in a virtual environment, keeping the computational time within reasonable limits. To accurately introduce the hemodynamic properties at the upstream and downstream areas of the specific region of interest, various BCs have been developed. These conditions, known as inlet and outlet BCs, are artificial boundaries that don't correspond to real boundaries in the original system. They play a crucial role in creating a realistic model of the entire cardiovascular system by accurately connecting the region of interest to the rest of the system while reducing computational demand. On the other hand, the artery wall is a physical boundary in the cardiovascular models. To accurately describe wall deformation under hemodynamic loads, the FSI technique is generally implemented at the solid-fluid interface. This method couples CFD simulation of hemodynamics with FEA simulation of wall mechanics. To precisely characterize the mechanical properties of the artery wall, researchers have developed constitutive models based on tensile tests of the arterial tissue.

Across multiple studies, certain patterns have emerged. Hemodynamic and biomechanical parameters such as TAWSS, OSI, ECAP, RRT, PWS and wall displacement are highly sensitive to the selected boundary conditions at the inlet, outlet, and wall. In the FSI studies of AAA with deformable walls, generally a time-dependent flow rate at the IR section of the aorta is supplied as fully developed, parabolic velocity profile at the inlet. The prescribed time-dependent pressure profile is given at the downstream of the external and internal iliac arteries. More sophisticated BCs at the inlet and outlet such as patient-specific or Womersley profile and 3-element Windkessel model are generally implemented to the pure hemodynamics studies neglecting wall deformation by implementing rigid wall assumption, to decrease the model complexity and computational cost.

Despite the advancements of integration of accurate BCs at the inlet, outlet and wall, integration and coupling of sophisticated BCs having various multiphysics and multidomain aspects is quite challenging. Defining patient-specific velocity profile at the inlet boundary is limited due to the cardiac motion and resolution of measuring devices. Coupling the lumped parameter models to the outlets requires patient-specific parameter estimation that necessitates development of different algorithms or prior CFD simulations. Accurate integration of wall motion is possible with proper discretization and solution of governing equations of mechanical domain using hyperelastic and anisotropic constitutive equations with patient-specific parameters, which requires excessive computational time. As an example, Lan et al. [54] have simulated the complete abdominal region from the SC zone by coupling WK3 BC at the outlet and FSI at the walls. In Figure 12, TAWSS and OSI comparisons of uniaxial, biaxial and rigid wall models are illustrated. Simulation was performed over three cardiac cycles, that take nearly 15 hours using three Intel Xeon Gold 5118 processors interconnected by a 100 GB/s EDR InfiniBand for a total of 72 threads operating at 191 GB RAM and a clock rate of 2.3 GHz. Qin et al.[191] modeled the complete aorta by dividing into four partitions by implementing the Newton-Krylov-Schwartz algorithm, assuming the walls as rigid and resistance BC at the outlet. They reported computational time as 87 hours for 2880 cores for $30.06\times {10}^6$ elements with nearly 70% parallel efficiency, even without extensive coupling of FSI and WK3. Later, they applied a highly parallel framework to study the effect of primary and peripheral branches on the local hemodynamics of abdominal aorta [284]. The results demonstrate that the peripheral branches significantly affect the flow field within the AAA which emphasizes the significance of complete modeling AAA domain, by focusing on the portion starts from SC and accounts for at least the primary branches such as CT, SM, LR, RR, ARR, LCI and, RCI arteries. However, the AAA studies with deformable wall BC typically start from IR region, by eliminating upper arterial branches to prevent excessive computational cost. Therefore, most of them utilize prescribed pressure waveform BC at the outlets, which are located downstream of left and right iliac arteries. Lan et al. [54] utilized a patient-specific computational domain with elastic and isotropic deformable walls, of which inlet section starts from SC region and the WK3 outlet boundary condition has been implemented.

To the best of author’s knowledge, the current literature lacks a model that incorporates FSI at the walls, begins from supraceliac (SC) region using 4D PC-MRI measured patient-specific inlet velocity profiles, and implements WK3 outlet BCs at the aorta branches. Furthermore, the vast majority of the studies using FSI together with WK3 implement hyperelastic and isotropic constitutive models, rather than the models account for wall anisotropy. Despite the significant advancements in computational models of AAA biomechanics, tuning the wall compliance parameter, $C,$ of WK3 model compatible with the hyperelastic and anisotropic nature of arterial walls is missing in the existing literature. Therefore, developing a method for tuning of WK3 parameters to reflect the wall hyperelasticity and anisotropy would be highly valuable for better understanding and predicting the behavior of arterial flow, wave propagation and wall mechanics in this region.

The primary purpose of researchers focusing on AAA hemodynamics and biomechanics is to provide clinicians with a tool that assists them in understanding the growth and potential rupture of the aneurysm during the surgical decision-making process. However, it is still a big challenge to accurately resolve the flow field and wall mechanics within a clinically acceptable computing time. Moreover, the limited resolution of medical imaging and measuring devices prevents accurate description of the local time-dependent velocity profiles and pressure waveforms, patient-specific geometry and wall thickness. These factors limit their practical use in real-time clinical practice, where timely interventions are crucial. Therefore, comparative studies are required to evaluate how simplified and less time-consuming models perform relative to the more accurate, computationally intensive counterparts. To integrate computational models into clinical workflow, development of simplified models which are validated against more accurate and patient-specific models is necessary. At that point, the ML tools and data-driven techniques might be quite helpful to provide clinical predictions due to their efficiency. However, to develop a reliable ML model, highly accurate CFD data sets should be produced by improving the precision of modeling parameters.

Current studies in the literature indicate that ML tools have significant potential for enhancing clinical predictions. This is due to their ability to efficiently process large volumes of data and identify complex patterns that might be overlooked by traditional methods. By leveraging highly accurate CFD data sets and improving the precision of modeling parameters, ML models can achieve a high degree of reliability. Consequently, integrating these advanced ML techniques into clinical workflows can result in more personalized and accurate treatment plans, ultimately improving patient outcomes and healthcare efficiency.