AAA geometry
Ten virtual aneurysm models were generated with the CAD software ProEngineer Wildfire (Parametric Technology Corporation, Needham, MA). The models differ in degree of asymmetry and wall heterogeneity, and are comprised by a fluid domain, ΩF, representing the aortic lumen and a solid domain, ΩS, representing the AAA wall. The fluid domain is characterized by circular cross sections parallel to the x-y plane, with an undilated diameter, d = 2 cm, at the inlet and outlet sections and a maximum diameter, D = 3 d, at the midsection of the AAA sac. The asymmetry of the model is governed by the β parameter given by Eq. (1).
where r and R are the radii measured at the midsection of the AAA sac from the longitudinal z-axis to the posterior and anterior walls, respectively, as shown in the inset of Figure 1(a). Thus, β = 1.0 yields an axisymmetric aneurysm. The geometry of the fluid domain is given by Eq. (2), which defines the diameter of each cross section, φ(z), and the deviation of its centroid from the z-axis, δ(z):
The geometry of the solid domain is given by an AAA wall with either a (i) uniform thickness (UW) or a (ii) variable thickness (VW). Both types of wall designs model the thickness of the aneurysm as material extruding normally from the surface enclosing the lumen. The UW model has a thickness given by λ = 1.5 mm, while the VW model is given by λ(z) = 1.5d/φ(z) (in mm) at each cross-section. It follows that the local wall thickness varies between 0.5 mm and 1.5 mm (with a mean of 1.0 mm for the AAA sac), inversely proportional to the local diameter of the cross-section. For a ruptured AAA, wall thickness can be as low as 0.23 mm at the rupture site with a surface-wide average of 1.45 mm [16]. The asymmetry of a virtual AAA is given by βε {1.0, 0.8, 0.6, 0.4, 0.2}, thus yielding a total of ten geometries, five of which have a uniform wall thickness (UW) and the other five variable wall thickness (VW). β = 1.0 corresponds to azymuthal symmetry, and β = 0.2 is an AAA for which only the anterior wall is dilated while the posterior wall is nearly flat. Figure 1(a) shows β = 0.6 with fluid domain ΩF entities in a uniform thickness model. Figure 1(b) illustrates qualitatively the variable thickness domain ΩS at three different wall locations for the same β = 0.6 model.
Governing equations and boundary conditions
The governing equations for the fluid domain are the continuity and Navier-Stokes equations with the assumptions of homogenous, incompressible, and Newtonian flow. Since the fluid domain is deformable in an FSI problem, an Arbitrary Lagrangian-Eulerian (ALE) formulation has been adopted. The ALE formulation [17] introduces a moving coordinate system to model the deformation of the fluid domain. The momentum and mass conservation equations governing the flow are given by Eq. (3) in ALE form:
ρfΔ·v = 0 (3b)
where ρf is the fluid density, τf is the fluid stress tensor,
are the body forces per unit volume, v is the fluid velocity vector, and
is the moving coordinate velocity, respectively. In the ALE formulation,
is the relative velocity of the fluid with respect to the moving coordinate velocity. Blood is modeled to have a density ρf = 1.05 g/cm3 and a dynamic viscosity μ = 3.85 cP. The governing equation for the solid domain is the momentum conservation equation given by Eq. (4). In contrast to the ALE formulation of the fluid equations, a Lagrangian coordinate system is adopted:
where ρs is the AAA wall density, τs is the solid stress tensor,
are the body forces per unit volume, and
is the local acceleration of the solid. The AAA wall is assumed to be an isotropic, linear, elastic solid with a density ρs = 2.0 g/cm3, a Young's Modulus E = 2.7 MPa and a Poisson's ratio υ = 0.45. The wall material implemented in this work represents a tissue of average characteristics for the aneurysmal abdominal aorta, i.e. a linearization of the stress-strain curve as reported in [13]. Previous studies have shown that aneurysm tissue is a non-linear, isotropic, hyperelastic material [18]. Hence, the constitutive linearity of the AAA wall is a simplification of the FSI and stress analyses in this study. In the FSI formulation, the AAA wall is assumed to undergo large displacements and small strains.
The boundary of the fluid domain is divided into the following regions for the assignment of boundary conditions: inlet (
), outlet (
), and the fluid-structure interaction interface (
), as shown in Figure 1(a). The applied boundary conditions on the non-FSI regions are (i) a time dependent fully developed velocity profile on
and (ii) a time dependent normal traction (due to luminal pressure) on
. These are presented by Eq. (5) as follows:
where τnn is the normal traction, u(t) and p(t) are the time dependent velocity and pressure waveforms shown in Figure 2,
designates the normal of the respective boundary, and I is the standard identity matrix. Time dependency, as introduced by u(t) and p(t), is given by Fourier series representations of the waveforms generalized in Eq. (6):
N is the number of harmonics used to reproduce the in vivo measurements of luminal velocity (N = 18), u(t), and pressure (N = 7), p(t), respectively. These waveforms are triphasic pulses appropriate for normal hemodynamics conditions in the infrarenal segment of the human abdominal aorta first reported by Mills et al [19]. The use of an input transient velocity based on normal physiology is justified by the fact that the inlet boundary condition is applied above the proximal neck of the aneurysm, an undilated segment of the abdominal aorta. For average resting conditions, blood flow in the abdominal aorta is generally laminar [20, 21]; flow deceleration achieved after peak systole induces laminar disturbed flow conditions and vortex formation even under simulated exercise conditions [22–24]. Inlet peak systolic flow occurs at t = 0.304 seconds and outlet peak pressure at t = 0.4 seconds. The time-average Reynolds number is Rem = 410, which is characteristic of a patient in resting conditions [25]. Rem is calculated as
, where
is the time-averaged, mean inlet velocity and d is the inlet diameter. The Womersley number,
, characterizes the flow frequency ω (ω = 2π/T and T = 1.0 seconds), the geometry and the fluid viscous properties, and is α = 13.1, a typical value for the human abdominal aorta under resting conditions [26]. The amplitude coefficient of the velocity waveform is defined as γ = Repeak/Rem = 5.25.
Figure 1(b) shows the boundary of the solid domain divided into inlet (
), outlet (
) and the fluid-structure interface (
) regions. The FSI interfaces
and
are identical, coupling the fluid and solid domains. The boundary conditions on the non-FSI regions of the solid domain impose zero translation on the ends
and
as given by Eq. (7). This corresponds to completely fixing the ends of the domain, simulating the tethering of the aorta by the surrounding tissue and organs. Numerical experimentation with the objective of minimizing the stresses at the proximal and distal necks yielded placement of the inlet and outlet sections at a distance 1.5d apart from the aneurysm sac.
The boundary condition at the outer wall surface corresponds to a reference zero normal traction, as the peritoneum and surrounding tissues do not exert any significant pressure on the arterial wall. There are no published data on normal forces exerted by internal organs and tissue on the wall abdominal aorta. The final set of boundary conditions is applied to the FSI interfaces
and
as follows: (i) displacements of the fluid and solid domain FSI boundaries must be compatible, (ii) tractions at these boundaries must be at equilibrium and (iii) fluid obeys the no-slip condition. These conditions are given by Eq. (8):
where d, τ, and
are displacement vectors, stress tensors, and boundary normals with the subscripts s and f indicating a property of the fluid and solid, respectively.
In addition to the FSI methodology, we investigated alternative computational solid stress (CSS) techniques for approximating the AAA wall stresses. In these analyses we disregard the blood flow and attempt to obtain comparatively accurate results by applying a spatially-uniform pressure function onto the inner wall. The CSST method (transient CSS) applies the transient function p(t) from Eq. 4(b) to simulate the effect of luminal pressure acting on the inner wall. Similarly, the CSSS method (static CSS) applies p(t = 0.4) in a quasi-static formulation to obtain the stresses at peak systolic pressure. In these two approaches we utilize only the solid domain ΩS as shown in Figure 1(b), with prescribed zero translation at the proximal and distal ends as given by Eq. (7), and with a pressure boundary condition as given by Eq.(9):
Numerical discretization
The software Adina (v8.0, ADINA R&D, Inc., Cambridge, MA) was utilized for the numerical simulation of fluid-structure interaction (FSI) between the wall and the lumen, as well as the alternative CSST and CSSS analyses involving only the aneurysmal wall, as described in [15]. The Finite Element Method (FEM) is used to solve the governing equations, which discretizes the computational domain into finite elements that are interconnected by nodal points. In this work we make use of linear hexahedral, eight-node elements to discretize the fluid and solid domains. The mesh generator Gridgen (Pointwise Inc., Fort Worth, TX) is used to develop the finite element grids. The ten aneurysm models are all composed of 17,280 hexahedral elements (19,093 nodes) for the fluid and 5,760 hexahedral elements (8,784 nodes) for the solid domains. Figure 3 illustrates the fluid and solid meshes for β = 0.6. Mesh sensitivity analyses were conducted with four additional mesh sizes ranging from 12,480 to 44,928 fluid elements and 3,840 to 13,824 solid elements. Independence in mesh size was obtained for the primary variables (velocity components, fluid pressure and structural displacements) within 5% relative error for the 4th mesh (32,256 fluid and 10,752 solid elements). However, the mesh used in the present study was chosen due to its adequate compromise between acceptable CPU simulation times (71 CPU-hours on average) and moderate relative errors of the primary variables at randomly selected nodal points (15% on average). In this regard, it is important to understand that this work is a baseline study for the development of a computational, pre-operative planning tool for physicians, and as such treatment decisions must be made within a reasonable turnaround time. Additionally, mass conservation in the fluid domain was met for a relative error (comparison between volume flow rate at the inlet and outlet sections, and rate of change of volume within the geometry) of 1%.
Equations (3–4) are reduced to weak form by following the standard Galerkin procedure [27]. In short, the fluid domain employs special Flow-Condition-Based-Interpolation (FCBI) hexahedral elements which use constant functions to interpolate velocity and bi-linear functions to interpolate pressure and displacements on each fluid element. The solid domain employs Mixed-Interpolation hexahedral elements, preferred in modeling nearly incompressible media, which use constant functions to interpolate pressure and bilinear functions to interpolate displacements on each solid element. The discretized equations for the fluid and solid elements are assembled into one system of equations, coupling the fluid and solid meshes. A sparse matrix solver based on Gaussian elimination is used for solving this system.
The FSI methodology utilizes an implicit time integration scheme, first applied on the fluid-structure interface of the fluid domain (
) where the coordinate system of the fluid and solid domains is Lagrangian. The results are then utilized to solve each domain entirely. Pulsatile flow is simulated over five to eight cycles with a time step size ΔtFSI = 5·10-3seconds until periodic convergence is achieved. Figure 4 shows convergence for the fluid domain (4a) and the solid domain (4b) for β = 0.2 in terms of five nodal point values of axial velocity and displacement magnitude, respectively. For the purpose of comparison, Figure 4(c) illustrates the convergence of displacement magnitude for the CSST analysis at the same five selected nodal points. The simulations were performed on a Tru64 Unix operating system using up to eight 1.15 GHz EV7 processors and in-memory computing. The CSS approaches only utilize the solid domain; hence, the final matrix assembly consists of only solid element equations and the computational times are significantly less when compared with the FSI computational times. For a consistent comparison, the CSST system of equations is solved with the same solution methods as in FSI but with a time step size
seconds. The CSSS system of equations is solved with the sparse matrix solver for a steady state solution with a constant and uniform pressure p(t = 0.4) applied on the inner AAA wall. The computational times for the CSSS simulations are negligible in comparison with the CSST (3 CPU-hours on average) and FSI (71 CPU-hours on average) simulations. A fully-coupled FSI method is computationally more expensive due to the memory required to adapt the matrices for a moving mesh algorithm and, thus, must be weighed against the clinical benefit of obtaining a complete characterization of the flow dynamics and wall stresses of the aneurysm sac.