Multiphase simulation of blood flow within main thoracic arteries of 8-year-old child with coarctation of the aorta

In the research, a numerical Computational Fluid Dynamics (CFD) model of the pulsatile blood flow was created and analysed. A real geometry of aorta and its thoracic branches of an 8-year old patient diagnosed with a congenital heart defect – coarctation of the aorta was used. The inlet boundary condition was implemented as the User Define Function according to measured values of volumetric blood flow. The blood flow was treated as multiphase using Euler-Euler approach. Plasma was set as the primary and dominant fluid phase, with the volume fraction of 0.585. The morphological elements (RBC and WBC) were set as dispersed phases being the remaining volume fraction.


Introduction
The untreated coarctation of the aorta (CoA)a narrowing of the descending aorta at the site of the ductus arteriosus insertion causes premature death of affected patients in young adulthood [1]. The blood flow in the normal aorta is laminar whereas in the case of CoA it becomes turbulent at the narrowed site and a higher pressure gradient occurs. The negative effect on the cardiovascular system comes from expo-sure of the upper part of the human body to hypertension and blood flow disturbances with all consequences. Typically, patients with coarctation of the aorta suffer from arterial hypertension, coronary artery disease, stroke, endocarditis, heart failure and aortic dissection. And the organs in the lower part of the body are hypoperfused [2]. While considering the indications for a therapeutic intervention in CoA patients physicians require an invasive blood pressure gradient measurement prior and below the narrowed site according to the American College of Cardiology and American Heart Association guidelines for treatment patients with congenital heart disease [3]. However, the invasiveness of these procedure presents a potential risk for the patient. Therefore it would be beneficial for the clinical practice to establish a non-invasive modality of assessing the coarctation significance. This can be accomplished using a combined data from imaging modalities, like Computed Tomography (CT) or Magnetic Resonance Imaging (MRI), together with blood flow modelling by means of Computational Fluid Dynamics (CFD).
Modelling of blood flow using CFD has recently become more popular and wider available. Development of such computer models opens new perspectives for predicting and solving of many potential problems. CFD achieved the status of the promising tool in bioengineering, where this technique is used for visualisation of blood flows within vessels and cavities of the circulatory system. A well validated numerical model can suggest solutions how to construct vascular stents with the best deployment, visualise blood flow within the artery prior to performing a percutaneous coronary intervention and estimate the significance of vascular lesions [4]. Despite many advantages, the CFD methods are still not ready for clinical use. This is mainly caused by the long computational time of such simulations (which for unstable flow phenomena requires using the complex unsteady solvers). This is an extended version of the conference contribution presented at 7th European Thermal-Sciences Conference (Eurotherm2016) [25] The most common approach used during modelling of the blood flow is a single phase simplification of the blood rheological character. In this approach properties of blood are represented by average density and viscosity for plasma and morphotic elements, whereas vessels are represented by rigid walls [5,6]. In literature, some numerical models describe flow within a selected section of the human circulatory system and an additionally implement a coupled lumped model for simulations of botha part and a whole cardiovascular system [7][8][9]. For the CFD simulations, it is crucial to establish a precise geometry of the flow field. The geometric data for the aorta is typically obtained from the CT or MRI scans. The preparation of geometric data for numerical simulations is an important task and requires special attention [10].
There are other mathematical models that can be used for modelling blood flow, like multifluid Euler-Euler (EE) or hybrid Euler-Lagrange (EL) techniques. The multiphase approach is used less frequently than single phase models. For the multiphase blood modelling with Euler-Euler technique, plasma and morphological elements are treated as continuous and interpenetrating phases. The EE approach is typically used for modelling narrow vessels where Fahraeus-Lindqvist effect occurs [11][12][13][14]. It helps to estimate the RBC volume concentration in the aortic lumen while single phase modelling does not allow for this The current publication presents an application of a three-fluid Euler-Euler approach for modelling of blood flow within the aorta and its main thoracic branches of an 8-year-old patient suffering from coarctation of the aorta (in the modelling the flow through the coronary arteries was omitted).

Numerical model
The geometry of the aorta was rendered from patients data acquired using Gadolinium-enhanced Magnetic Resonance Angiography (MRA). The narrowing of the aorta in this patient caused approximately 65% aortic area reduction. The numerical case was based on geometry and data defined in CFD challenge problem for development numerical models capable to model blood flow (OSMSC Corp., USA) [15]. The investigated geometries after conversion to solid model using GeomagicDesignX (3D Systems, Inc., USA) software are shown in two versions in Fig. 1. Case (b) covers additionally in comparison with Case (a) bifurcation of the brachiocephalic trunk (BCT) into the right subclavian artery (RSA) and the right common carotid artery (RCCT).
The final generated mesh consists of about 700,000 of hybrid elements. In the inlet part, mostly hexahedral elements were used in order to stabilise flow and accurately prescribe velocity profile. The mesh was verified using two meshes characterized by different densities.
The crucial point in the modelling of blood flow is to accurately define the boundary condition at the inlet to the computational domain. In presented work, the velocity inlet, based on the volumetric flow profile, has been set to mimic the conditions during the human cardiac cycle. The whole cycle has been divided into 20 points (with linear interpolation between them), where volumetric flow values were taken from a phase-contrast (PC) MRI through-plane velocity encoding [16]. The periodic pulsatile boundary condition was implemented by a set of User Define Functions (UDFs) using piecewise polynomial functions in subsequent time intervals.
The volumetric flow profiles for specific blood phases are depicted in Fig. 2. The hematocrit of 0.405 (vol.) was used in numerical simulations. WBC contained 0.01 of the volumetric share, while plasma constituted the rest of mixture in the amount of 0.585 (volumetric). It is a value within the normal range. The patient specific rheology can be represented in the model by modifying the inflow volume fraction. Preliminary analysis shows that the results are sensitive to the dominating solid phase content i.e. RBC. Two vertical red lines represent characteristic moments of the cardiac cycle (systole and diastole). The simulation covered whole cardiac cycle, while results are discussed for these two specific time instants.
The distribution of blood flow at the outflows were implemented in the model as outflow boundary condition (cf. Table 1) for the two analysed cases of geometry and were based on data published in [16]. Description of the shortcuts used in the table is pointed in Fig. 1 caption. The blood flow quantity on the outflows were calculated from mass conservation principle using the inlet information.
The presented research can be compared with the publication [17]. In this work, the single phase model was presented, in the same geometry as presented here geometry in Case (b). The cited research covered also 3-element Windkessel model attached at each outflow instead of the constant mass flows distribution proposed in [16]. In the current research, the blood distribution in the region of the brachiocephalic trunk bifurcation was estimated proportionally to the mass flow shares calculated in [17]. One can notice that the mass flows distribution differs between this two cases (see Table 1). In authors' opinion, the flow distribution is more probable using a model with Windkessel sub model presented in [16]. Nevertheless, in all cases presented in Table 1, the disproportions of the blood flow occurring between LSA and RSA seems not to be physiological. There are based on external data (volumetric flow distributions and resistances for Windekessel model) and authors assumed their correctness.
The characteristic parameters used in the numerical model definition are presented in Table 2. It contains material properties of specific phases and the coefficients describing walls behaviour. In the numerical simulations, the walls of the domain were treated as rigid.
In the research, the governing equations describing multiphase flow are presented in the following paragraphs. The volume fraction of fluid and solid phases sums up to 1 according to Eq. (1).
The principle of mass conservation in the multiphase flow is described by Eq. (2) for the fluid phase, i.e. Plasma with subscript f and by Eq. (3) for solid phases, where subscript s represents RBC and WBC phases respectively. In the current research, neither mechanism of mass exchange between specific phases nor any mass sources is considered, yielding the right side of continuity equations being zero.
The principle of momentum conservation in the multiphase flow can be written for the fluid phase as in Eq. (4). To Subscript s represents RBC and WBC solid phase respectively. Subscript q is RBC or WBC depending on configuration (while q ≠ s), N s expresses the number of solid phases (here N s = 2). P represents the pressure of fluid and P s is a granular pressure. The g ! is gravitational acceleration, τ ! ! is stress tensor and K is the interphase exchange momentum coefficient, F is depicted as additinal forcess which can be caused by modelling of the virtual mass force, lift force etc.
Granular temperature parameter used in the kinetic theory is associated with the random motion of particles or, in other words, particle velocity fluctuations. It is a component needed for the calculation of the random kinetic energy of solid phases according to (6) based on [12,18,19]. This value depicted as θ s , is used to calculate closure terms, for instance, granular pressure (collision, kinetic and friction), granular viscosity (collision, kinetic and friction), and solid stress tensor. RBC pressure was represented by p s symbol. Restitution coefficient was depicted as e ss , g 0 , ss was radial distribution function which in this configuration was calculated from (8). The diffusion coefficient for granular energy k θ s can be expressed by (7). 3 To solve interaction between phases, the influence of RBC and WBC presence on plasma flow a number of closure terms in Eqs. (4), (5), (6) and (7) have to be defined. The detailed description of the closures terms responsible for modelling granular pressure, granular and bulk viscosities, etc. can be found in [20,21]. The momentum exchange coefficients (K) between plasma and RBC treating as a dense granular flow were calculated applying Gidaspow [20] drag model. Whereas, the WBC was treated as the dilute granular phase, so Wen&Yu drag model [22] was applied to predict interphase momentum exchange. The Gidaspow model was used to predict drag between RBC and plasma phases because the RBC volume fraction exceeded 0.2. In that case the combination of Ergun and Wen&Yu model should be used. In the case of WBCs their volumetric fraction is much less than 0.2 so the Wen&Yu model can be used for drag modelling.
The turbulence effect was taken into account by selecting the standard k-ε model, where the kinetic (k) and dissipation (ε) of energy were calculated per each phase [24]. Within the pulsatile cycle, the domain of turbulent flow changes both in time and space. It is not realistic to account for it by switching in and out the turbulence model according to the local Reynolds number. The critical portion of the analysis is the flow within the coarctation where turbulent flow dominates. Thus, although the flow may be laminar in the remaining parts of the domain, specifically close to diastole, the turbulence model has been invoked in the entire domain. As there is practically no recirculation, isotropic model like k-eps are admissible. The advantage of the mentioned model is its numerical stability. The restitution coefficient (e), shown in Table 2, is the parameter of the mathematical model that describes interactions/collision of RBC and WBC with walls [12,13]. This value is close to 1 while kinetic energy between WBC and RBC before and after their collision almost did not change.

Results
Presented relative static pressures are referred to the pressure on the descending thoracic aorta domain outlet. The pressure field on the domain walls for the characteristic point during systole is presented in Fig. 3. The highest pressures of about 5 kPa occurred in the region of the inlet. During systole, the smallest static pressure can be noticed in the aortic coarctation and its relative pressure reached −5.7 kPa. Thus, the maximum pressure difference in the computational domain for systole was on the level of 11 kPa. Results are consistent for Case (a) and Case (b) for inlet and DA outlet. Differences between cases can be noticed in the bifurcation region of the smaller arteries. Case (a) is characterised by higher pressure in the region of BCT, LCCA. According to the pressure field, presented in Fig. 4, for diastole characteristic point the maximum pressure difference reaches 40 Pa for Case (a). In this case, the highest pressure was at the outlet of the descending aorta and was also referred to its whole field. The lowest pressure for diastole was on the main inlet and before coarctation, In Case (b) during diastole, the pressure differences are higher between specific outflows and inlet, and it can be the results of different densities of specific phases. The static pressure level was significantly smaller in comparison with systole for the both cases. During diastole, the volumetric flow was nearly zero. The pressure gradient was reversed in comparison to systole (as shown earlier).
The velocity field for this specific time instants is illustrated by means of the velocity vectors in Fig. 5. The velocity vectors in the branching region of BT, LCCA and LSA from the aorta is shown during systole on the left-hand side of the figure. To clarify the flow field view the velocity was scaled to 2 m/s. The highest velocity values occurred in the CoA area and reached over 4 m/s. Based on the velocity field analysis, it can be deduced that during systole, the flow was stable without reversed flows at the outflows boundary conditions. From the inlet, the blood was flowed to the specific outflows according to set boundary conditions. The simulated flow looked mostly laminar. Some eddies occurred in aorta distally to the point where the left subclavian artery comes off and proximally to the coarctation of the aorta. The smallest velocity values were present in the close vicinity of the blood vessel walls. Chosen model of turbulence was assumed due to high Reynold numbers in the specific cross section of the blood vessels. In the inlet cross section for the systole, Reynolds no. reached more than 6200. Moreover, in the region of coarctation and outlet from descending aorta, Reynolds no. achieved more than 4200. Additional results with turbulence eddy dissipation in the domain was presented in multi cross sections in Fig. 6. It is visible that the highest values of eddy dissipation are in the region of the coarctation and reach more than 50 m 2 s −3 .
During diastole, the reversed flow occurred only for individual cells wall on each outflow as seen from the numerical simulation. The highest velocities occurred prior the  Fig. 8 Volumetric fraction of the RBC phase expressed in cross sections of the ascending aorta (left) and the narrowing of the descending aorta (right) during systole coarctation of the aorta and reached 0.15 m/s. This maximum velocity is 28 times smaller than during systole. The velocity field analysis can indicate that flow was chaotic and irregular during diastole. In the presented area, numerous eddies are visible.
Case (b) was characterised by similar velocity filed as Case (a), therefore only velocity vectors during systole in the region of BC bifurcation is presented in Fig. 7. In this figure, maximum scale of velocity was also limited to 2 m/s to clarify the presented velocity field. One can notice that the velocity field did not change in comparison to geometry presented at the beginning. The velocity of RBC in the domain consisted with the plasma velocity. Nevertheless, the velocity in Case (b) can be estimated in the additional branch of BC. The area weighted average velocity in this region reached about 1 m/s in the right subclavian artery and 1.3 m/s in the right common carotid artery.
Conducted analysis allows for calculation of the pressure and velocity fields in the numerical domain. Simplified Case (a), without brachiocephalic bifurcation, showed similar results as Case (b). Maximum velocity in the coarctation cross-section varied no more than 0.5%. Therefore, it could be deduced that simplified case is sufficient in the coarctation analysis.
The volumetric fraction field of phase being the greatest fraction beyond plasma was presented in Fig. 8. The simulation results in a cross section of ascending and descending aorta were shown. The smallest values of the RBC volume fraction, dominating in cross-section area, starts from walls vicinity and occur in the plane of coarctation region. Migration of RBC from the walls to the centre direction could be caused by rising granular temperature in the walls vicinity according to multiphase kinetic theory [20]. The maximum values of RBC volume fraction occurred at some distance from the walls and reaches 0.4067. The volumetric fraction in the presented cross section was consisted with case (b).

Discussion and conclusions
The scope of the present study was to show an application of computer model for modelling of blood flow within the aorta and its major thoracic branches of an 8-year old patient diagnosed with a congenital heart defect -coarctation of the aorta. Two main types of morphological elements -RBCs and WBCs were represented by two separate granular phases modelled in the Euler -Euler approach, while the primary phase was the plasma. The fields of pressure, velocity and volume fractions in selected region were presented for the two characteristic time instants of the heart cycle (during systole and diastole). Moreover, two simulations, covered basic geometry and its extended version, were conducted during the research. The second geometry contained bifurcation of the brachiocephalic trunk and was an enrichment in comparison with the first one.
The blood flow distribution differences between left and right subclavian arteries could be verified only by obtaining of medical data. Presented paper was focused on the multiphase Euler-Euler approach application in investigated geometry with specified assumptions.
The future research could be conducted using multiphase approach and Windkessel sub-models on the outflow boundary conditions.