Three patients (56–68–69-year-old, 1 female and 2 males) with aortic diseases from the Union Hospital, Tongji Medical College, Huazhong University of Science and Technology were included in this study. In addition to abdominal aneurysm or dissection, two patients (Patients 1 and 2) exhibited arterial stenosis at the origins of certain visceral arteries due to atherosclerosis or calcification. The third patient (Patient 3) presented with a focal narrowing of the descending aorta caused by aortic dissection. Patient details are summarized in Table 1, and more detailed anatomical measurements are available in the Supplementary Material.

The study protocol was approved by the institutional review board of Union Hospital, Tongji Medical College, Huazhong University of Science and Technology. Since the data involved in this study is retrospective, the requirement for informed consent was waived and anonymized data was used.

Computed tomography angiography (CTA) sequences of three patients were acquired using a 3rd generation Dual-Source CT scanner (SOMATOM Force; Siemens AG, Erlangen, Germany; 100 kV, 250 mA s, rotation time: 0.25 s, field of view (FOV): 250 mm, slice thickness: 0.5 mm, reconstruction kernel: Bv40, contrast agent: Iomeron 400; Bracco, Milan, Italy) with in-plane resolution of 0.5 mm and inter-slice distance of 0.5 mm for the reconstruction of the aorta, as shown in Supplementary Figure S1. Pre-treatment brachial blood pressure and heart rate were also measured and are summarized in Table 1.

The aorta geometries were reconstructed from CTA images using the image-processing software DetecModeling (Boea Wisdom, Hangzhou, China) (Ding et al., 2024). The reconstructed aorta, starting from the ascending aorta and extending to the internal and external iliac arteries, was cut and patched to form the inlet and outlet boundaries of the computational domain.

To obtain a patient-specific healthy aortic model, the lesional aorta was manually repaired to generate a near-healthy reference geometry. To ensure the validity of the geometric repair, all modifications were independently performed by two experienced technicians under the guidance of a radiologist and a vascular surgeon. In cases of disagreement, a senior technician provided the final judgment.

The repair process of the diseased aorta is illustrated in Figure 1. Since the CTA scans of the diseased aorta were obtained relatively soon after the onset of the disease, its overall structure remained largely unaffected (Li et al., 2023; Wen et al., 2023). Therefore, the original vascular framework was preserved during the repair process. Furthermore, as the arterial disease had minimal impact on the location and geometry of bifurcation points (Ford et al., 2009), their anatomical structures were retained in the reference model. For cases of aortic dissection, the intimal flap was first filled to merge the true and false lumens. Subsequently, the merged model underwent contraction and smoothing, as the presence of the intimal flap and blood perfusion in the false lumen could contribute to aortic expansion. In cases of abnormal stenosis or dilation, the vessel diameter was adjusted based on the diameter of unaffected segments within the same anatomical region. A comparison between the lesional geometry and the repaired near-healthy reference geometry is provided in Supplementary Figure S2.

The spatial discretization of both the lesional and reference geometries was carried out using the mesh tool integrated in DetecFluid (Boea Wisdom, Hangzhou, China). Given the complex morphology of the aorta, tetrahedral elements were employed for meshing. To ensure grid independence, three different grid resolutions were tested (see Supplementary Material), and a suitable configuration was selected for each aortic model. In the final meshes, the number of cell elements ranged from approximately 2.5–4.4 million across all patient-specific models, providing sufficiently fine resolution to capture the hemodynamic features. Detailed mesh statistics are provided in Supplementary Material.

As illustrated in Figure 2, an iterative framework based on steady-state CFD simulations was developed to estimate the flow division in the downstream branches of diseased aortic models. Steady-state simulations were first conducted on the repaired reference geometry to obtain inlet and outlets pressures. The outlet pressures were subsequently converted into distal resistance models and prescribed at the outlets of the diseased geometry. Subsequently, steady-state simulations were performed on the diseased model to predict the inlet pressure, which was compared against the target value. Based on this discrepancy, the outlet resistances were updated, and the simulations were repeated until convergence was achieved. The boundary conditions employed in each simulation are described in detail below.

In the steady-state CFD simulation of the reference model, the mean systolic flow rate, derived from empirical flow waveforms presented in Figure 2, was prescribed at the inlet, as this approach has been shown to effective predict mean blood pressure and vascular resistance (Li and Mao, 2023). The downstream flow division was assigned based on widely used measurements from previous studies (Je and Dn, 1994; Moore and Ku, 1994). Specifically, 30% of the cardiac output was allocated to the supra-aortic branches, with further distribution among them proportional to their cross-sectional areas. The celiac trunk received 15.5% of the cardiac output, while the superior mesenteric artery, the left renal artery, and the right renal artery each accounted for 10.5%. If the inferior mesenteric artery was present, it received 0.5% of the cardiac output. The remaining flow was equally divided between the left and right iliac arteries, with 30% directed into the internal iliac artery and 70% into the external iliac artery. Although these flow division values are based on imprecise historical measurements, they remain widely accepted in hemodynamic studies (Bonfanti et al., 2017; Bonfanti et al., 2019). Importantly, this does not affect the validation of the proposed method. More accurate flow division data from healthy individuals could be incorporated into future studies as they become available.

Based on the inlet pressure, outlet pressures, and flow rates predicted by the steady-state CFD simulation of the reference model, along with the empirical mean blood pressure values, the resistance at each outlet in the resistance model can be calculated using Equation 1: R i = P mean − P in − P out , i Q i (1) Where R i represents the resistance of the ith branch artery, mmHg·s·ml^-1; P mean is the empirical mean blood pressure, which will be adjusted according to the discrepancy between the predicted inlet pressure and the target inlet pressure of the diseased aorta; P in is the predicted inlet pressure of the reference model, mmHg; P out , i is the predicted pressure at the ith branch outlet of the reference model, mmHg; and Q i is the flow rate at the ith outlet of the reference model, ml·s^-1.

Next, a steady-state CFD simulation was conducted on the diseased aorta. The inlet boundary condition was set as a mass flow rate equal to that of the reference model. At the outlets, the resistance model was imposed, allowing the outlet pressure to vary dynamically with the flow rate. This simulation provided a predicted inlet pressure for the diseased aorta, which was then compared to the target value. Given that only brachial blood pressure was available for each patient, an empirical formula was used to estimate the mean pressure at the aortic inlet, as presented in Equations 2–4 (Stokes et al., 2021): P mean , ao = 0.4 P sys , ao + 0.6 P dia , ao (2) P sys , ao = 0.83 P sys , bra + 0.15 P dia , bra (3) P dia , ao = P dia , bra (4) Where P mean , ao is the mean pressure at the aortic inlet, mmHg; P sys , ao is the aortic systolic pressure, mmHg; P dia , ao is the aortic diastolic pressure, mmHg; P sys , bra is the brachial systolic pressure, mmHg; and P dia , bra is the brachial diastolic pressure, mmHg.

The threshold for the discrepancy between the predicted and target inlet pressures was set at 1%. If the difference exceeded this threshold, the empirical mean pressure in Equation 1 was adjusted, and the steady-state CFD simulation of the diseased model was repeated. The iteration continued until discrepancy fell below 1%, at which point convergence was achieved. At this stage, the outlet resistance model was deemed representative of the patient’s physiological condition, and the predicted flow division in the aortic branches was considered a more patient-specific estimation.

Transient CFD simulations were performed for each diseased aorta to analyze its hemodynamic characteristics. The inlet boundary condition was prescribed using a flow rate waveform (Figure 2) from the literature (Gallo et al., 2016), scaled according to patient-specific physiological data, and imposed with a plug flow profile. At the outlets, the three-element Windkessel was implemented, with its parameters estimated using the fast approach proposed by Li and Mao (Li and Mao, 2023). In the estimation process, the previously predicted patient-specific flow division in the aortic branches was used.

Additionally, to compare and assess the performance of the proposed method, a traditional scheme based on empirical flow division from healthy individuals (Bonfanti et al., 2017; 2019), specifically, the outlet flow rate used in step 1 in Figure 2, was also employed for Windkessel parameter estimation. In the following text, the results obtained using the previously predicted flow rate are labeled as the “proposed method”, while the results from the traditional scheme are labeled as the “traditional method”.

All CFD simulations were conducted by the software DetecFluid (Boea Wisdom, Hangzhou, China). It is a highly efficient GPU-based software for solving the equations of continuity and momentum conservation using a finite-volume methodology (Liu et al., 2025). Blood was assumed to be an incompressible Newtonian fluid with a density of 1,060 kg m^-3 and a dynamic viscosity of 0.0035 Pa s. The vessel wall was defined as rigid with a no-slip boundary condition. Given the relatively high blood velocity and potential turbulence in the aorta, the standard k - ε model (Moukalled et al., 2016) was employed due to its numerical stability and its demonstrated applicability in aortic hemodynamics (Long Ko et al., 2017; Etli et al., 2021; Hohri et al., 2021). The standard wall function was used to resolve the viscous sublayer. A second-order scheme was used for spatial discretization and a first-order implicit scheme for temporal integration. The SIMPLE algorithm (Moukalled et al., 2016) was used to couple pressure and velocity equations. The converged residuals for all variables were 10^–4. Additionally, for the transient CFD simulation, the cardiac cycle was divided into 1,000 time-steps. The transient simulation was performed over five cardiac cycles, and the analysis was conducted based on the data from the final cycle.

Several metrics related to blood pressure, flow patterns, and wall shear stress were used to characterize the aortic hemodynamics. For temporal statistical metrics, every tenth timestep from the final cycle was used.

Figure 3 illustrates the flow division in aortic branches predicted by the traditional and proposed methods. The specific flow rate values for each branch are provided in Supplementary Table S3. A strong correlation was observed between the two methods in predicting aortic branch flow division (Patient 1: r = 0.865, Patient 2: r = 0.973, Patient 3: r = 0.991). However, noticeable differences were found in certain branches, particularly those with severe stenosis, such as the celiac trunk in Patients 1 and 2. The proposed method predicted lower flow rates in these arteries, which is considered more physiologically realistic, resembling the “steal” phenomenon observed in patients with artery occlusions that redistribute blood flow to meet metabolic demands (Yamashiro et al., 2002; Hatzidakis et al., 2015). The Bland-Altman plot further indicates that the differences in predicted flow fractions for some branches approach or exceed the 95% limits of agreement.

The total resistance R _WK and compliance C of the Windkessel models estimated by the traditional and proposed methods differed, as shown in Figure 4. These variations in parameters are related to the severity of the lesions and the differences in target flow rates in each branch. In the proposed method, the flow rate in the stenotic branches decreased and was redistributed to other branches. As a result, the total resistance of most Windkessel models in Patients 1 and 2 decreased, while the resistance of severely stenotic branches, such as the celiac trunk, increased. The variation in capacitance values of the Windkessel model was more complex, as it was influenced by both the pressure changes due to flow rate variations and volume changes associated with vascular lesions. Among all branches, the greatest change in capacitance occurred in the right renal artery of Patient 3, where it decreased from 0.577 mL mmHg^-1 in the traditional method to 0.022 mL mmHg^-1 in the proposed method.

The inlet pressures of the diseased aorta in Patient 1, predicted by the traditional and proposed methods, are compared and shown in Figure 5. The predicted inlet pressures differed from the measured aortic systolic and diastolic pressures. The proposed method performed better in reproducing the inlet blood pressure. In contrast, the traditional method assigned flow rates typical of healthy conditions to the stenotic branches, resulting in non-physiological pressure losses that influenced the simulation outcomes. The predicted systolic and diastolic pressures for Patients 2 and 3 are summarized in Supplementary Table S3. For Patient 3, the predicted inlet pressure from both methods significantly differed from the actual values. This discrepancy may be related to the limitations of the tuning strategy used in Windkessel models (Li and Mao, 2023).

The blood pressure distributions on the selected planes in Patient 1 at five time points are shown in Figure 6. The values predicted by the proposed method were higher than those predicted by the traditional method, and the distribution patterns were similar. Additionally, the blood pressure variations over time during the cardiac cycle were comparable. Similar trends were observed in Patients 2 and 3, as shown in Supplementary Figures S3,S4.

The TMP values for Patients 1 and 3, calculated by both methods, are shown in Figure 7. The values predicted by the proposed method differed from those predicted by the traditional method, with the traditional method generally overestimating the absolute TMP values. The difference was more pronounced and significant (p < 0.01) in Patient 3, primarily due to two severe stenoses in the true lumen of the descending aorta, located around planes d and g. These stenoses had a broader impact, leading to the observed discrepancies. Both methods predicted the same balance points where TMP equaled zero: for Patient 1, this balance point was between planes e and f, and for Patient 3, it was near plane h. The maximum difference in TMP values predicted by the two methods exceeded 5 mmHg, occurring near the iliac bifurcation in Patient 1, around plane j.

The flow patterns in the aorta of Patient 1 are presented here. Results for Patients 2 and 3 can also be found in the Supplementary Material. Figure 8 illustrates the flow patterns in Patient 1 at five representative time points, and Figure 9 shows the velocity distribution on the selected planes. Unlike the blood pressure distributions, the velocity distributions predicted by the two methods were similar. Aorta tears cause a significant reduction in the cross-sectional area of the TL, but most of the blood still flows into the TL (86.5% in the proposed method and 88.4% in the traditional method), resulting in a high blood velocity in the TL. For branches supplied by both TL and FL, there was a slight but not significant change in the proportion of source. In the REIA, the proportion of blood flow from the TL increased from 85.4% in the traditional method to 87.4% in the proposed method.

Like velocity, the distributions of LNH predicted by the two methods are similar. To maintain conciseness, these distributions are not shown here and can be found in the Supplementary Material. We calculated the difference in LNH by subtracting the values predicted by the traditional method from those predicted by the proposed method, as shown in Figure 10. The differences in LNH were scattered sporadically throughout the aorta, with more pronounced variations observed in the FL, which has a relatively large volume. From a temporal perspective, the differences were more pronounced during mid-systolic acceleration and the diastole phase compared to other times. This may be linked to the strong transition flow occurring between systole and diastole (Liu et al., 2012; Poelma et al., 2015).

The distributions of WSS-related metrics in Patient 1 are shown in Figure 11, with those for Patients 2 and 3 available in the Supplementary Material. The distributions predicted by the two methods were similar, especially for the metrics based on cycle-averaged data such as TAWSS, OSI, and RRT. This is understandable, as these indicators are derived from wall shear stress, which is determined by the velocity gradient near the vessel wall. As previously mentioned, the velocity distributions along the aorta trunk predicted by the two methods were comparable. In the TL, high blood velocity resulted in high and unidirectional WSS, characterized by higher TAWSS, lower OSI and RRT, and more uniform AFI. In contrast, the FL exhibited low and oscillating WSS, with lower TAWSS, higher OSI and RRT, and more varied AFI, due to the larger volume and lower flow rate.

We calculated the difference in these metrics by subtracting the values predicted by the traditional method from those predicted by the proposed method, as shown in Figure 11. In the celiac trunk, the TAWSS predicted by the proposed method was relatively lower due to the reduced flow rate. Differences in OSI were observed in areas near the visceral branches, where flow patterns were complex, while differences in RRT were scattered sporadically throughout the FL. For AFI, significant differences were widely distributed throughout the aorta, especially in the FL.

More realistic predictions of blood flow distribution enable more precise identification of malperfusion syndrome in patients. This represents a significant advancement over past studies that used simulations with zero-pressure outlets to assess the risk of malperfusion. Cheng et al. (2013) utilized transient simulations with zero-pressure outlets to analyze the hemodynamic characteristics of four patients with acute TBAD, two of whom developed complications of malperfusion. They suggested that higher flow rate into the FL might be related to malperfusion syndrome. However, their models did not include downstream branch arteries. Furthermore, a high flow rate into the FL does not necessarily result in malperfusion, as a patent FL can still supply blood to the distal arteries (Bonfanti et al., 2019). Alimohammadi et al. (2014) developed a patient-specific simulation tool for aortic dissection, in which the downstream flow division was first determined using a steady-state CFD simulation with zero-pressure outlets to tune the Windkessel model parameters for multi-scale simulations. Their results showed that 39.4% of the inlet flow went through the BT, leaving only 45.0% for the descending aorta, which contributed to lower limb and/or vital organ malperfusion. However, simulations using zero-pressure outlets likely overestimate the risk of malperfusion in downstream branches because the flow division was directly determined by the geometric resistance in the reconstructed aorta model. For comparison, we conducted simulations with zero-pressure outlets for Patient 1, following the methodologies of Cheng et al. (2013) and Alimohammadi et al. (2014). Table 2 illustrates the flow rate in branches predicted by these simulations. The “0Pa_Windkessel_target” values represent the flow division predicted by the steady-state simulation with zero-pressure outlets. The “0Pa_Windkessel” values denote the flow rates predicted by the multi-scale simulation with Windkessel models, whose parameters were estimated using the flow division in “0Pa_Windkessel_target”. The “0Pa_Transient” values show the flow division predicted by the transient simulation with zero-pressure outlets. It is evident that the flow distributions predicted by these simulations are similar, with approximately 48% of the inlet flow passing through the BT, similar to Alimohammadi et al.‘s study (Alimohammadi et al., 2014). The flow rates in the three branches on the arch were much higher than those in distal branches, which is clearly inconsistent with the physiological condition.

Furthermore, more accurate predictions of blood pressure and WSS-related metrics allow for a more precise evaluation of aortic diseases progression risks, including the obstruction of the TL (Bonfanti et al., 2019; Xu et al., 2021), the evolution of the FL (Chen et al., 2013; Xu et al., 2018), the dilation and rupture of aneurysms (Teng et al., 2022), and the formation of luminal thrombosis (Naim et al., 2016).

As introduced before, the proposed method involves manual repair of the diseased aorta model, which can introduce additional errors. Given that the level of smoothing in model reconstruction significantly impacts hemodynamic predictions (Paritala et al., 2023), we conducted a preliminary investigation into the effect of the smoothing factor on the results. Specifically, the repaired aorta of Patient 1 was smoothed to varying degrees, and the flow division in branches was recalculated using the proposed method. The results, presented in Table 3, indicate that smoothed models yielded more consistent flow rates with the proposed method and deviated significantly from the flow division predicted by the traditional method. Future work should focus on developing automated model repair techniques to strengthen the robustness of the proposed method. Ford et al. (2009) developed an objective and automated technique for digitally repairing arteries with saccular aneurysmal lesions. However, for more complex aortic conditions, such as fusiform aneurysms or aortic dissection, efficient automated repair algorithms remain unavailable.

There are several limitations in this study. First, while we showed that the proposed method could reproduce inlet blood pressure closer to the physiological data of patients, some persistent differences remained. These discrepancies arose because the Windkessel model parameters were estimated by a steady-state simulation with mean systolic blood flowrate at the inlet. This strategy is more effective for healthy aortas, where geometric resistances of branches vary similarly under different inflow rates. However, in stenotic branches of diseased aortas, the variation in resistance with inflow rate differs significantly from that of vessels with normal diameters, leading to less satisfactory predictions. In Li and Mao’s study (Li and Mao, 2023), the difference between the reproduced inlet pressure and physiological data also increased after the artificial narrowing of a branch. In addition, estimating Windkessel parameters from steady-state simulations has an inherent drawback: arterial compliance is neglected. As a result, the predicted inlet pressure ( P in ) in Equation 1 is overestimated, leading to an underestimated resistance term ( R i ). When these parameters are subsequently applied in transient simulations, the Windkessel model tends to underestimate the mean arterial pressure across the cardiac cycle. Fortunately, the framework proposed in this study is adaptable to any Windkessel model tuning strategy that requires target flow rates, including those based on 0D simulations (Stokes et al., 2021). In the future, efforts should focus on improving strategies to better reproduce patient physiological data.

In addition, several common assumptions and simplifications in CFD simulations, such as the rigid vessel wall, empirical inflow waveform combined with a plug flow profile, and Reynolds-averaged turbulence models, may reduce accuracy. Previous studies have shown that assuming a rigid wall tends to overestimate flow velocity and wall shear stress (Ene et al., 2014; Bonfanti et al., 2018). Compared with healthy aorta, neglecting wall compliance in diseased aortas with local stenoses, as considered in this study, may introduce larger errors. On the one hand, wall deformation can substantially alter the local resistance at the stenotic site and thereby affect flow distribution; on the other hand, in some pathologies, such as aortic dissection, intimal flap motion may cause dynamic obstruction (Kim et al., 2023), leading to significant changes in hemodynamics. Simplifications in the inlet boundary condition, particularly the use of a plug profile, have also been shown to compromise the accuracy of wall hemodynamic predictions in the ascending aorta (Stokes et al., 2023a). Moreover, the Reynolds-averaged turbulence model adopted in this study (the standard k - ε model) may not adequately capture the transition between laminar and turbulent flow. Performing large-eddy simulations with finer meshes would enable the resolution of smaller-scale vortical structures and yield more accurate wall hemodynamic characteristics (Manchester et al., 2022; Cheng et al., 2025). Nevertheless, while these assumptions and simplifications may limit the absolute accuracy of the simulations, they do not undermine the demonstration that the proposed approach outperforms the traditional one.

Finally, the actual flow rate in a patient’s aortic branches is also regulated by the sympathetic nervous and endocrine systems (Leito-Rocha et al., 2012), which were not explicitly considered in this study. The sympathetic nervous system is known to regulate blood pressure by controlling peripheral resistance and cardiac output (Luo, 2003). In this study, different inlet pressures were assigned to the repaired and diseased aorta, implying that the effect of the sympathetic nervous system was, to some extent, implicitly considered. If physiological parameters such as blood pressure and cardiac output in the healthy state were available, the proposed method could be further refined. Additionally, future studies should validate the accuracy of the proposed method using clinically measured flow data, for example, obtained by four-dimensional flow magnetic resonance imaging (4D-Flow MRI), two-dimensional phase-contrast magnetic resonance imaging (2D PC-MRI), or Doppler ultrasound.