Fast and Accurate Computation of the Displacement Force of Stent Grafts after Endovascular Aneurysm Repair

Purpose: Currently, the displacement force of stent grafts is generally obtained using computational fluid dynamics (CFD), which requires professional CFD knowledge to perform the correct simulation. This study proposes a fast, simple, and clinician-friendly approach to calculating the patient-specific displacement force after endovascular aneurysm repair (EVAR). Methods: Twenty patient-specific post-EVAR computed tomography angiography images were used to reconstruct the patient-specific three-dimensional models, then the displacement forces were calculated using CFD and the proposed approaches, respectively, and their numerical differences were compared and analyzed. Results: Based on the derivation and simplification of the momentum theorem, the patient-specific displacement forces were obtained using the information of the patient-specific pressure, cross-sectional area, and angulation of the two stent graft ends, and the average relative error was no greater than 1.37% when compared to the displacement forces calculated by CFD. In addition, the linear regression analysis also showed good agreement between the displacement force values calculated by the new approach and CFD (R = 0.999). Conclusions: The proposed approach can quickly and accurately calculate the patient-specific displacement force on a stent graft and can therefore help clinicians quickly evaluate the post-EVAR displacement force.


Introduction
An abdominal aortic aneurysm (AAA) is a pathological expansion of the abdominal aortic segment, which is described as the blood vessel diameter expanding by greater than 50% or greater than 3 cm [1]. Currently, endovascular abdominal aortic aneurysm repair (EVAR) is the first preference for the therapy of AAAs [2] due to its higher early remedial impact as compared to normal open surgical treatment [3]. However, compared with normal open surgical repair, EVAR has higher risks of reintervention due to endoleaks and SG migration [4][5][6], especially in the case of complex anatomical morphology, such as short length, large-angle, and conical AAA necks [7][8][9].
The migration of SGs is greatly associated with the magnitude and direction of the displacement force (DF) on the SG [10], and if the magnitude of the DF exceeds the ability of the proximal and distal end of the SG to prevent migration, the SG may migrate, which may lead to proximal seal failure, aneurysm rupture, and limb thrombosis of the SG [11][12][13]. Therefore, the evaluation of post-EVAR DF, especially the maximum DF of the overall graft, is helpful for better stent designs and assessment of EVAR treatment outcomes [14].
Computational fluid dynamics (CFD) has been widely used in the study of AAAs to observe intravascular flow fields and to calculate the DF on the SG following EVAR [15][16][17][18]. However, although the CFD approach can obtain an accurate value of the DF by integrating the SG surface pressure and wall shear stress [19][20][21], its procedure is complicated and tedious, as follows: (1) First, the computational model needs to be reconstructed based on patient-specific CTA images. (2) Second, the 3D models need to be meshed appropriately. (3) Third, time-consuming calculations and simulations with appropriate boundary conditions must be carried out to obtain the pressure field and wall shear stress on the SG. (4) Finally, the pressure and wall shear stress are integrated to obtain the DF. For these procedures, whether 3D model meshing, calculation, or post-processing, professional CFD background knowledge is required, but clinicians generally lack knowledge of CFD techniques. Therefore, a fast, simple, and clinician-friendly approach to obtaining the patient-specific DF on the SG is urgently needed for clinicians to judge whether the SG may migrate after the operation.
Assuming blood is steady and non-viscous, Liffman et al. and Mohan et al. developed a mathematical model for studying the DFs on a planer and symmetrical bifurcated graft using continuity and momentum equations [22,23]. However, due to their simplifying assumptions about the blood and graft geometry, their theoretical model can only give a rough estimation of the axial force acting on a graft; it seems feasible to use the momentum theorem to directly solve the DF. Encouraged by their preliminary studies and aimed at proposing a fast, simple, and clinician-friendly approach obtaining the DF on the patientspecific SG model, the present study first starts from the fundamental theory of fluid mechanisms to derive a reasonable simplified formula based on the momentum theorem, then twenty patient-specific post-EVAR computed tomography angiography (CTA) images were used to reconstruct the patient-specific three-dimensional (3D) models. Finally, the DFs were calculated using CFD and the proposed approach, respectively, and their numerical differences were compared and analyzed.

Materials and Methods
This was a retrospective observational study that was conducted following the principles of the Helsinki Declaration and met the requirements of medical ethics. The Ethical Review Committee of the West China Hospital of Sichuan University (Chengdu, Sichuan, China) approved this research. Since all the data were collected from retrospective anonymized databanks and as our study was purely observational, patient approval and informed consent were waived.

Geometry
Twenty patient-specific SG models, shown in Figure 1, were established using thin-slice CTA images (scanner: SIEMENS/SOMATOM Definition Flash, slice thickness: 1.0 mm, pixel size 0.6133 mm) with the commercially available Mimics software (version 15.0; Materialise, Plymouth, UK).

Governing Equations
The blood was assumed to be incompressible, laminar, homogenous, and Newtonian [24,25], with the corresponding governing equations given as:

Governing Equations
The blood was assumed to be incompressible, laminar, homogenous, and Newtonian [24,25], with the corresponding governing equations given as: where → u and p represent, respectively, the fluid velocity vector and the pressure. ρ and µ are the density 1050 kg/m 3 and dynamic viscosity 3.5 × 10 −3 kg m/s of the blood, respectively [26,27].

Boundary Conditions
The no-slip condition was imposed on the surface of the models according to the rigid wall assumption [28]. It was shown that the maximum DF and peak pressure drop arose at peak pressure throughout the cardiac cycle [29,30]. As a result, a constant velocity of v = 0.1 m/s was applied at the inlets, and the pressure was set as 105 mmHg at the outlets, which characterized the moment of peak pressure throughout the cardiac cycle [31].

Meshing and Computing
Computational meshes were established using ANSYS ICEM 14.5 (ANSYS, Inc., Canonsburg, Pennsylvania, PA, USA). The meshes contained a mixture of tetrahedral and hexahedral volume meshes. As a result of the boundary layer effect of fluid, the mesh near the wall needed to be densified. To decrease the difference in the calculation results caused by the different degrees of mesh refinement, we adopted the Grid Independence Index (GCI) to evaluate the generated mesh [32]. When the GCI of all the test variables was less than 2%, it was considered that there was a small enough spatial discretization error without further mesh refinement [33,34]. The ultimate number of cells for the meshes ranged from 1.5 million to 2.0 million.
The CFD calculations were performed using ANSYS Fluent 14.5 (ANSYS, Canonsburg, Sylvania, Pennsylvania, PA, USA). The calculation model was a laminar flow model. A simple algorithm was used as the calculation method for pressure and velocity coupling; the second-order upwind scheme was used to discretize the momentum equation and energy equation; the finite volume method was used to solve the continuity equation and momentum equation.

Calculation of Displacement Force
The DF was calculated by taking an area integral of the net pressure and wall shear stress (WSS) over the entire wall of the SG: where p is pressure, → τ is the WSS vector, A pdA is the pressure force, and A → τ dA is the WSS force.
The angle between the DF vector and the positive direction of the space rectangular coordinate system can be calculated by the following formula: where F x , F y , and F z are the component sizes of the DF projected on the x-, y-, and z-axes respectively. θ x , θ y , and θ z are the angles between the DF vector and the positive direction of the x-, y-, and z-axes, respectively. The x-direction is the right (positive) and left (negative) direction; the y-direction is the front (positive) and rear (negative) direction; the z-direction is the up (positive) and down (negative) direction.

Governing Equations
The DF can be also calculated by the momentum theorem, which only needs the flow information at the inlet and outlet when the flow is steady. Thus, we may calculate the DF using the momentum theorem as follows: 3) represent the pressure, area, blood flow velocity vector, and unit normal vector of the inlet and outlet sections, respectively, as shown in Figure 2a. The length of the unit normal vector is 1, which mainly represents the spatial angle of the section, so the unit normal vector can also be expressed as where α, β, and γ are the angles between the normal vector and the coordinate axis, as shown in Figure 2b.  In this study, the direction of the unit normal vector was always perpendicular to the inlet and outlet section towards the SG.
If we substitute the boundary condition of CFD into Equation (7), the ratio p i A i /ρA i U 2 i ranges from 76 to 16,396, indicating that the DF caused by the momentum change is negligibly small as compared to the pressure force. Therefore, Equation (7) can be simplified as follows: According to Equation (9), the DF of the SG is decided by the pressure force difference between the two SG ends (Figure 2c). However, it was found that the DF contributed by the pressure drop (ranging from 0.25 to 2.86 mmHg in the current models) only accounted for less than 2% of the total DF. Therefore, Equation (9) may be further simplified to Equation (10).
where p is the patient-specific systolic blood pressure. As a result, as long as the pressure, the cross-sectional area, and the angle of the SG ends are measured, the patient-specific DF can be quickly obtained.

Measurement of the Cross-Sectional Area and Angles of SG Ends
(1) Measurement of the cross-sectional area of the SG ends: The model's surface area was derived directly from the Mimics software. Therefore, after the SG model (model 1) was established, any part (model 2) could be extended from the section to obtain the whole of model 3. Using the Boolean operation between the three models, the cross-sectional area was obtained (Figure 3a).
Bioengineering 2022, 9, x FOR PEER REVIEW 6 of (1) Measurement of the cross-sectional area of the SG ends: The model's surface a ea was derived directly from the Mimics software. Therefore, after the SG model (mod 1) was established, any part (model 2) could be extended from the section to obtain t whole of model 3. Using the Boolean operation between the three models, the cros sectional area was obtained (Figure 3a).  (2) Measurement of the normal vector of the section: We determined the coordinates of three points on the cross-section at both ends of the SG model, and then, the start vector and end vector were determined according to the three points (Figure 3b).

Statistical Analysis
The SPSS software (v21.0, IBM Inc., New York, NY, USA) was used for the statistical analysis. Comparing the relative error (relative error = absolute error/true value × 100%) between the total DF calculated by the momentum theorem and the CFD method, it was considered that the value calculated by the CFD method was the true value in this study. Linear regression analysis was carried out to evaluate the correlation between the simplified formula of the momentum theorem and the CFD method. The coefficient of determination (R-value) was used to measure the fit of the predicted value to the true value. The closer the R-value was to 1, the closer the predicted value was to the true value. A Bland-Altman analysis was carried out to compare the two measurement methods; additionally, the average deviation and limits of agreements (defined as 1.96× standard deviations) are provided [35,36].

Results
The magnitude of the DF calculated using the momentum theorem (Equation (10)) and CFD method (Equation (3)) is shown in Table 1, and the angle of the DF calculated is shown in Table 2. The average relative error of the DF received by the momentum theorem and CFD method was 1.37%. In particular, the maximum relative error in the model of SG6 was 6.91%. On θ x and θ z , the average relative errors of all the models were 1.04% and 1.43% respectively, whilst the maximum relative errors in the model of SG6 were 13.60% and 19.91%, respectively. The relative error of all the models on θ y was small; the average relative error was 0.21%, and the maximum relative error was only 1.04%. F is the displacement force; F x , F y and F z represent the components of F in the x-, y-, and z-directions, respectively. Relative error (%) = absolute value of (CFD value − momentum value)/CFD value × 100%. θx, θy, and θz represent the angles between the force vector F and the positive direction of the rectangular coordinate system. Relative error (%) = absolute value of (CFD value − momentum value)/CFD value × 100%.
A linear regression analysis was carried out on every calculated DF dataset to further evaluate the correlation between the two methods. Figure 4 suggests that there was a precise consistency between the magnitude of the DF calculated by the simplified momentum theorem and the CFD method (R = 0.999), as well as the DF angle (R > 0.993). θx, θy, and θz represent the angles between the force vector F and the positive direction of the rectangular coordinate system. Relative error (%) = absolute value of (CFD value − momentum value)/CFD value *100%.
A linear regression analysis was carried out on every calculated DF dataset to further evaluate the correlation between the two methods. Figure 4 suggests that there was a precise consistency between the magnitude of the DF calculated by the simplified momentum theorem and the CFD method (R = 0.999), as well as the DF angle (R ＞ 0.993).   The Bland-Altman plot ( Figure 5) shows that there was only a small deviation between the magnitude of the DF calculated by the momentum theorem and the CFD method (the upper and lower limits of DF consistency calculated by the two methods were 0.111 N and −0.13 N, respectively); only the DF of SG6 exceeded the limit. The distinction between the DF angles calculated by the two methods was additionally very small (the upper and lower limits of θ x consistency were 5.71 • and −7.19 • , respectively; the upper and lower limits of θ y consistency were 0.68 • and −0.83, respectively; the upper and lower limits of θ z consistency were 5.97 • and −8.05, respectively). The SG6 patient exceeded the limits of θ x and θ z . The SG7 and SG8 patients slightly exceeded the limit on θ y .

Discussions
Endovascular abdominal aortic aneurysm repair (EVAR) is used to insert an interventional device with a stent graft (SG) from the femoral artery and then delivering it to the subrenal position to release the stent. In this way, a new flow channel is constructed to isolate the blood flow, prevent the blood from directly scouring the aneurysm wall, and treat an abdominal aortic aneurysm. Since Volodos et al. first established EVAR technology 1987 [37], it has since become more mature with continuous improvements over more than 30 years to the structure and anchoring method of the SG [38][39][40][41]. However, the use of EVAR is prone to causing endoleaks and stent migration when AAA patients have unfavorable anatomies such as severe aneurysm necks.
At present, CFD is used to obtain the DF on the SG, the acting force produced by the blood to the SG, which is the main cause for the SG's migration. However, the CFDbased DF calculation method requires professional CFD knowledge for accurate simulation. Moreover, the conventional CFD approach is complex and time consuming [42]. In this study, we proposed a simple and fast approach to obtaining the global DF using the momentum equation, which only needs the patient-specific blood pressure, the crosssectional areas, and the angulations of the proximal and distal ends of the SG. The results showed that the DFs calculated by the proposed simplified formula were in good agreement with those obtained by CFD (R > 0.993).
Agreeing with previous CFD studies [23,29,42], the current study showed that the contribution of the change in blood flow velocity to DF is negligibly small. Therefore, the blood momentum inside the SG can be regarded as unchanged, and the DF of an SG is mainly decided by the pressure force difference at its two ends (Equation (9)). This may explain why the DF waveform trend of SGs is always consistent with that of the pressure waveform in previous CFD research [43][44][45]. In addition, although the force is a vector and the change in its magnitude and direction both may contribute to the force difference, this study indicated that the contribution by the pressure drop between the proximal and distal SG accounted for less than 2% of the total DF. That is, the DF on the SG is mainly due to the direction change of pressure forces at the SG ends, and it is reasonable to evaluate the DF using the patient-specific blood pressure (Equation (10)). Furthermore, because the pressure force is the product of the pressure and area, high blood pressure and large areas of SG anchor ends may both cause a large pressure force and, thus, a potentially large DF. As a result, not only is hypertension a risk factor for SG migration after EVAR [22,29,42], but an increase in the area of the proximal and distal SG can also result in an increased risk of adverse events [22,46]. In practice, a proximal oversizing ratio of greater than 30% has been suggested to lead to the risk of stent migration [47,48], and the use of a clock bottom design at the distal end was found to increase the incidence of type Ib endoleaks [49].
This study revealed that the angulation of the proximal and distal SG plays a key role in the total DF on the SG. If the angulations of the proximal and distal SG are not significantly different, the pressure forces at the two ends will cancel each other out, resulting in a small DF on the SG (the SG8, SG9, and SG19 patients; Table 1). On the contrary, if the angulation of the two anchor ends is very different, the DF on the SG will be large (the SG4 and SG18 patients; Table 1). Therefore, minimizing the angle difference between the two SG ends can reduce the risk of displacement following EVAR. Previous studies have indicated that changes in the aorta angulation caused changes in the magnitude and direction of the displacement force [19,50,51]. This may be because the current SG placement generally allows the anchor end to conform to the shape of the blood vessels, so that the DF of the SG will be directly affected by the aortic morphology.
There are still some deficiencies in our research. First of all, there were only 20 patientspecific models reconstructed and tested. Although the results reflect the accuracy of our measurements and calculations to some extent, they need to be verified by a much larger cohort of patients. Secondly, because the CFD calculation showed that the pressure drops of the two SG ends were small, the proposed formula did not account for the contribution of the pressure differences. However, if an implanted SG is excessively distorted, the pressure drop between the proximal and distal ends will increase, and its contribution to the DF will increase. Additionally, the current study showed that if the pressure difference increases by 1mmHg, the calculation error by the simplified equation will increase by 0.8%. Therefore, in the case of the high distortion of an SG, the displacement force caused by a pressure drop should be properly considered. Last but not least, the density and dynamic viscosity used in the CFD simulations represent parameters of cold blood (at a temperature of about 4°C). There may be some errors if the parameters of living humans are not used. However, the simplified formula removes the term regarding density. In addition, the human blood viscosity ranges from 0.0040 ∼ 0.0050 kg·m/s, so such an error should be relatively small and within a controllable range.

Conclusions
In this study, we proposed a simplified formula to calculate the DF of the SG. This approach only requires the section parameters of the 3D-reconstructed SG model, and then, the DF of the SG can be directly obtained by a simple calculation. Compared with the traditional CFD method, this new approach achieved a similar level of accuracy without the need for complex model post-processing or computer simulations.