## Abstract

The right atrium (RA) combines flows from the inferior (IVC) and superior vena cava (SVC). Here RA mixing is simulated using computational fluid dynamics, comparing four modeling approaches. A patient-averaged model (11 M cells) was created from four volunteers. We compared: (1) unsteady k–*ω* Reynolds-averaged Navier–Stokes (URANS) (2) implicit large eddy simulation with second-order upwind convection scheme (iLES-SOU) (3) iLES with bounded-central difference convection scheme (iLES-BCD) and (4) LES with wall-adapting local eddy-viscosity (LES-WALE). A constant inlet flow rate of 6 L/min was applied with both IVC/SVC contributions ranging from 30–70%. A higher density mesh (37 M cells) was also simulated for models 2 and 4 (equal IVC/SVC flow) to assess the accuracy of models 1–4. Results from the 11 M cell LES-WALE model showed good agreement with the 37 M cell meshes. All four 11 M cell models captured the same large-scale flow structures. There were local differences in velocity, time-averaged wall shear stress, and IVC/SVC mixing when compared to LES-WALE, particularly at high SVC flow. Energy spectra and velocity animations from the LES-WALE model suggest the presence of transitional flow. For the general flow structures, all four methods provide similar results, though local quantities can vary greatly. On coarse meshes, the convection scheme and subgrid-scale (SGS) model have a significant impact on results. For RA flows, URANS should be avoided and iLES models are sensitive to convection scheme unless used on a highly resolved grid.

## Introduction

The right atrium (RA) combines flows from the coronary sinus, inferior vena cava (IVC), and superior vena cava (SVC). These flows then pass through the tricuspid valve to the right ventricle. RA hemodynamics has been well characterized in several 4D-magnetic resonance imaging (4D-MRI) imaging studies [1–3] and relatively few computational fluid dynamics (CFD) studies [4–6]. A particular focus of the CFD studies to date has been to assess the performance of cannulas for veno-venous extracorporeal membrane oxygenation. In this context, CFD is a powerful design tool allowing for the testing and iterative improvement of cannulas in the complete range of potential operating conditions.

Historically, computational studies of the hemodynamics of healthy vessels have assumed laminar flow [7,8], neglecting the effect of unresolved eddies. However, turbulence models are commonplace in CFD studies of diseases like atherosclerosis [10], aortic valve stenosis [11], or defects like aortic coarctation [12–14]. Of the few CFD studies of RA hemodynamics, some have assumed fully laminar flow throughout the domain, solving the Navier–Stokes equations on relatively coarse meshes (<1.5 M cells) with no additional modeling for subgrid-scale eddies [5,6,9]. It is common to quote relatively low Reynolds numbers (Re < 2300) at the model boundaries as a justification for the omission of any turbulence modeling. The assumptions of the Reynolds criterion (Newtonian viscosity, steady flow, and straight circular pipe geometry), however, do not hold well for blood flows and hence turbulence can be expected at much lower Reynolds numbers [8,10,11]. Other studies of the RA have applied some form of turbulence modeling [4,12]. Observations from in vivo imaging data often suggest the presence of turbulent or disturbed flow in the RA though these qualitative descriptions fail to characterize whether this is transient or fully developed turbulence [2,13]. Detailed CFD analysis of flow unsteadiness can provide quantitative information about the nature of the flow, including assessing potential turbulent structures.

In this study we apply unsteady Reynolds-averaged Navier–Stokes (URANS), implicit large eddy simulation (LES) (with both second-order upwind (SOU) and bounded-central difference (BCD) convection schemes) and LES with subgrid-scale (SGS) modeling and assess how well each approach capture the flow scales appearing in the RA using the same grid resolution.

## Methodology

### Patients and Imaging.

Four healthy volunteers (three female, one male, and aged 58.3±3.5 years) were imaged at the Department of Radiology at Karolinska University Hospital. Two computed tomography acquisitions were made on a dual source multidetector Somatom definition flash (Siemens Healthcare, Forchheim, Germany) with the injection of iodinated contrast media and saline to capture the brachial and iliac veins, SVC, IVC, and RA. These images had a slice thickness of 0.75 mm. All patients gave informed consent and ethical approval was obtained from the Swedish Ethical Review Authority (Ethical permit 2018/438-31). This imaging data has been used in a previous study of RA physiology [15]. This study is unique in its focus on modeling and flow characterization.

### Three-Dimensional Reconstruction.

A three-dimensional reconstruction of the veins and RA was created for each subject in the open-source software 3dslicer (version 4.10.2) by a single analyst. Manual editing to remove imaging artifacts and minor veins was performed in star-ccm+ (version 15.06.011, Siemens, Munich, Germany). From these reconstructions, a patient-average model was created, with the maximum and minimum vein diameters reflecting the mean values for the group. Similarly, the dimensions of the RA and the distance from the iliac and brachial veins to the RA were scaled to reflect the mean.

### Meshing.

where *ɛ* = relative error, *r* = mesh refinement ratio, *p* = order of convergence and $f1,\u2009f2,\u2009f3$ are the hemodynamic quantities calculated by the three meshes of decreasing cell count.

Three meshes of increasing density (2.5 M, 8.8 M, and 20 M cells) were simulated by applying an LES with a wall-adapting local eddy-viscosity (WALE) SGS model and a BCD convection scheme. GCI values for the time-averaged velocity at the line probes entering the RA are acceptably low for both, the medium and fine meshes giving mean GCI values of 2.29% and 0.25%, respectively (Fig. 1(c)). Time-averaged velocity convergence was better for flows entering the RA from the SVC (probe 1) compared to the IVC (probe 2). Convergence of wall shear stress (WSS) was also assessed to be sufficient and is presented in the Supplemental Materials (Fig. S1) on the ASME Digital Collection. When applied to the full geometry, these parameters resulted in an 11.2 M cell mesh. The average grid size in the RA, and so the filter size applied to the LES simulations, was 0.22 mm. To assess mesh resolution in the LES-WALE model, the proportion of resolved turbulent kinetic energy to total turbulent kinetic energy (resolved + modeled) was calculated over 1 s of physical time. The mean of this ratio was 97.4% in the RA for the 50/50 flow scenario, indicating a high-resolution mesh capable of resolving the majority of the energy containing eddies.

A further refined mesh of 37 M cells was simulated at the 50/50% flow scenario for the implicit LES with second-order upwind convection scheme (iLES-SOU) and LES-WALE models. With an increasing mesh density, the solution will approach that of a direct numerical simulation (DNS) (here 37 M cells resulted in an average grid size of 0.08 mm in the RA). Consequently, the choice of discretization scheme and SGS model is expected to be of less importance. The 37 M cell mesh results were compared with the 11 M cell results (also at 50/50% flow) to establish which model should serve as the point of comparison for the investigation of all IVC/SVC flow rates.

### Boundary Conditions.

Inlets at the coronary sinus, brachial, iliac, and renal veins were cut perpendicular to the vessel centerline and extruded to at least 10 times the arterial diameter (Fig. 1(a)). The tricuspid valve was assumed to be circular; of average area (2.5 cm^{2}) [14] and was also extruded to 10 times its diameter. Constant inlet flow rates were applied at all inlets and caval flow rates were distributed by inlet area. Total inlet flow was 6 L/min. The natural balance of flows into the RA in healthy adults is approximately 65% from the IVC and 35% from the SVC (18) though this can vary greatly, especially in the critically ill. The proportion of flow entering the RA from the SVC and IVC was therefore varied from 30-70% each, in increments of 5%. The coronary sinus was assigned an additional 2% of total cardiac inflow [16]. A zero-pressure condition was applied at the tricuspid valve. A non-Newtonian Quemada viscosity model for blood was applied with a density of 1050 kg/m^{3} and a constant hematocrit of 35%.

### Flow Modeling.

The complete range of inlet conditions was simulated using:

SST

*k*–ω URANS.iLES-SOU—commonly used when the laminar flow is expected.

iLES-BCD.

LES model with wall-adapting local eddy-viscosity (LES-WALE) subgrid model and BCD convection scheme.

**I**is the identity tensor, $T\xaf$ is the mean viscous stress tensor, and $fb$ are the resultant body forces. $TRANS$ can be represented by Boussinesq's approximation (Eq. (5))

**S**is the mean strain rate tensor. Turbulent eddy viscosity can then be calculated by several methods, most commonly the

*k*–

*ω*model first proposed by Wilcox et al. [18] (Eq. (6)) and the

*k*–

*ɛ*model [4] (Eq. (7))

where *k* is the turbulent kinetic energy, *T* is the turbulent time scale (a function of the specific dissipation rate, *ω* in Eq. (6) or the turbulent dissipation rate, *ɛ* in Eq. (7)), $C\mu $ is a model coefficient and $f\mu $ is a damping function [17]. To harness the relative strengths of the two methods, Menter's SST *k*–*ω* model utilizes the *k*–*ω* form at the inner boundary layer and *k*–*ɛ* elsewhere [19].

where $\Delta $ is the grid filter width and *S* is the modulus of the mean strain rate tensor (of $v\u0303$). LES requires greater computational resources than URANS, needing to resolve the time and length scales of nonfiltered eddies, but allows for the observation of transient structures. Instead of using LES coupled with an SGS, the computational grid can act as a filter, i.e., iLES, in which the grid resolution directly determines the length scales captured by the simulation. Assuming laminar conditions, several studies in the literature have adopted an approach similar to iLES but using a coarse grid resolution, truncating the simulation and resolving only part of the flow spectra omitting all turbulent flow scales.

All models were implemented in the commercial CFD solver star-ccm+ with the segregated flow solver which applies the semi-implicit method for pressure-linked equations algorithm for pressure and velocity coupling [17]. The flow field was initialized with steady SST *k*–*ω* RANS for (1.) and (4.) whilst the steady laminar flow model in star-ccm+ was used to initialize (2.) and (3.). Simulations were run for 1.5 s of physical time with a constant time-step of $10\u22125$s. Momentum and continuity residuals converged to <$10\u22125$ and *y*+ values were <1 for all simulations. These were run on the Mahti Supercomputer (CSC, Espoo, Finland). Simulation times for the URANS, iLES-SOU, iLES-BCD, and LES-WALE models were approximately 17.2, 12.3, 15.8, and 31.9 thousand core hours, respectively.

#### Hemodynamic Metrics.

Results for time-averaged wall shear stress (TAWSS) and velocity were averaged over the final 1 s of simulation. To assess the impact of turbulence modeling on the mixing of IVC and SVC flows within the RA, a passive scalar was applied to the inlet flows. This passive scalar was designated a value of 1 at the SVC and 0 for the IVC and can be thought of as a virtual contrast agent applied to each flow. This was also averaged over the final 1 s of the simulation. A plane through the center of the RA was defined to compare the time-averaged velocity and passive scalar quantities. The surface of the RA was divided into the auricle and central atrium segments to compare TAWSS results at these locations. The known thrombogenicity of the auricle region in the left atrium [21,22], and reports of thrombosis in the right auricle made the local shear conditions here of particular interest [23]. The energy spectra were calculated for several points of interest in the geometry for the 30/70, 50/50, and 70/30 flow scenarios. Velocity components were collected at these locations every 0.001 s for an additional 2 s of simulation (1000 Hz sampling rate).

## Results

On the 37 M cell mesh, the LES-WALE and iLES SOU models showed good agreement. Comparing the 11 M and 37 M grids, iLES-SOU showed an appreciable difference in the direction of the SVC jet entering the RA (Fig. 2). On the coarse mesh, the SVC jet was less diffuse and directed toward the center of the RA. The fine mesh results showed better agreement with the 11 M cell iLES-BCD and LES-WALE models. The coarse LES-WALE model was therefore used as the point of comparison for all subsequent coarse mesh results run for all IVC/SVC flow rates.

### Right Atrium Structures.

Figure 3(a) shows the time-averaged velocity streamlines for the IVC/SVC% = 65/35 flow scenario, the expected ratio for healthy adults reported by Kuzo et al. (ten subjects, mean age = 30 years) [24]. All four models showed a large vortical structure in the RA and a smaller counter-rotating structure in the auricle, with the orientation and position of these structures differing slightly. Time-averaged velocity values at the midplane differ significantly from LES-WALE for both, iLES and the URANS, particularly when SVC flow dominated (25–45%) (Figs. 3(b) and 3(c)). Except for the 55/45% scenario, iLES-BCD had the lowest relative difference compared to LES-WALE. Both, URANS and the iLES showed lower velocities (compared to LES-WALE) at a region of flow separation (point 1, Fig. 3(b)) at SVC flow rates >55% and higher velocities in the central atrium (Point 2, Fig. 3(b)). Time-averaged velocity results for all simulations are presented in full in the Supplemental Materials (Fig. S2) on the ASME Digital Collection.

### Wall Shear Stress.

TAWSS values computed with the URANS and the iLES models differed greatly from LES-WALE (Figs. 4(a) and 4(b)). This relative difference was generally reduced as IVC inflow increased except for the 60/40 flow scenario. Here the auricle TAWSS values obtained from URANS diverged sharply from LES-WALE values (117%) (Fig. 4(b)). As expected, TAWSS results from the iLES-BCD model were closest to the LES-WALE results (*P* < 0.005). TAWSS contour plots for all simulations are presented in full in the Supplemental Materials (Fig. S3) on the ASME Digital Collection.

### Blood-in-Blood Diffusion.

Relative differences in the calculation of the passive scalar representing IVC/SVC mixing were generally high when the SVC flow rate was high for both, the URANS and iLES (SOU and BCD) simulations as compared to LES-WALE (Figs. 5(a) and 5(b)). As the proportion of flow entering from the IVC increased, URANS and iLES results began to reflect those obtained from LES-WALE more closely (Fig. 5(b)). Both iLES models showed particularly good agreement (<5% surface average difference) when IVC flow accounted for >55% of total inflow. Underestimation of SVC blood at point 3 in Fig. 5(c) was very high for both, the URANS and the iLES models at most flow rates compared to LES-WALE, however, this also improved in the iLES models when IVC flow is >55%. The URANS and iLES passive scalar values were much closer to LES-WALE values at points 1 and 2, particularly at high proportions of IVC flow. IVC/SVC mixing results for all simulations are presented in full in the Supplemental Materials (Fig. S4) on the ASME Digital Collection.

### Flow Length Scales and Turbulence.

Energy spectra displayed in Fig. 6 differ between the iLES-SOU and LES-WALE models. The slope of the energy cascade is steeper at points 3 and 5 indicating significant dissipation at these locations. Energy spectra for the 30/70 and 70/30 flow conditions are included in the Supplemental Materials (Fig. S5) on the ASME Digital Collection. Time and length scales for the WALE model are presented in Fig. 7(a) and a comparison is made between the effect of the turbulent viscosity (left) and dynamic viscosity (from the Quemada model) (right) in the LES-WALE model, under 50/50 flow conditions (Fig. 7(b)). The impact of turbulent viscosity is clearly concentrated in the near-wall region, particularly at the region of flow separation.

## Discussion

In this study, we simulated a four-patient-averaged model of the RA and vena cavae under a range of SVC and IVC flow rates using URANS, implicit LES with two convection schemes, and LES with a WALE SGS model. Whilst the entire range of IVC/SVC flow rates were run on a more computationally efficient 11 M cell mesh, we compared results at the 50/50 flow scenario with results from a 37 M cell mesh. On the denser mesh, iLES and LES-WALE models converged on a common solution and showed good agreement with our coarse mesh LES-WALE model, which serves as our point of comparison for all subsequent simulations. For the four CFD models, we compared results for time-averaged velocity, TAWSS, and a passive scalar representing the mixing of SVC and IVC flows. All four models produced the same large-scale structures though local velocities differed greatly. Agreement in these velocity fields improved as IVC flow increased. The relative differences (compared to LES-WALE) in TAWSS calculations in the RA were high for both, the URANS and iLES models. Mean TAWSS in the RA showed better agreement across the four models as IVC flow increased. Relative differences in the calculation of the passive scalar representing SVC and IVC mixing were significant for both, the URANS and the iLES models as compared to LES-WALE, though the iLES-BCD results were similar to LES-WALE when IVC flow was ≥50%. The region of flow separation where the SVC flow diverts away from the atrial wall was particularly susceptible to high relative differences. The energy spectra for both LES models at several points in the RA suggest the presence of transitional flow.

### Large Vortex Structures.

Computational fluid dynamics studies have great potential in the analysis of RA hemodynamics because of their relatively high spatial resolution compared to 4D-MRI and the fact that boundary conditions can be manipulated to reflect a range of vena cava flow rates. In this study, all modeling approaches captured the same large-scale structures (Fig. 3(a)). Though orientations differed slightly between simulations, a large rotating structure comprising both, SVC and IVC flow was visible in the central atrium, consistent with previous imaging studies [1,2]. In addition to this, the four models all show a smaller counter-rotating structure in the auricle. The small variances in the location of the structures are the likely reason for high relative differences when comparing velocity fields, though these were not so high when the IVC flow rate was high. As the proportion of flow entering from the SVC was increased, a region of flow separation (Fig. 3(b), point 1) became more prominent. Though this region has not previously been studied in the RA, the separation of flow occurring downstream from a stenosed carotid artery serves as a useful comparison. Johari et al. compared the results from SST *k*–*ω* RANS with LES (Smagorinsky SGS model) and experimental particle image velocimetry for a case of carotid artery stenosis noting reasonable agreement in velocity results between all three methods [25]. The velocity profiles obtained from LES however more closely reflected the experimental data throughout the flow separation region. iLES-BCD appears to better characterize the flow separation zone compared to iLES-SOU on the 11 M cell mesh (except for the 45/55% flow scenario). BCD is known to be more accurate for complex turbulent flows [17,26]. The impact of the convection scheme is clear in Fig. 2 where the solution from the iLES-SOU model was much closer to the LES-WALE result when the cell count was increased. The results from this study suggest that URANS and both iLES models provide a reasonable indication of large vortex structures in the RA but velocities at specific locations, on coarser grids, are likely to be impacted by the SGS turbulence model and choice of convection scheme, especially at high SVC flow rates.

### Wall Shear Stress.

Wall shear stress is a common metric to report in CFD studies of blood flow. When exposed to high shear stress red blood cells can become damaged (hemolysis) impacting their ability to deliver oxygen [27,28]. Furthermore, both, abnormally high and low wall shear stresses impact the layer of endothelial cells which line our vessels [29]. Extremes of wall shear stress play a role in the development of thrombosis and emboli [22,30]. The results from this study show that the modeling approach can have a vast impact on WSS calculations. The relative difference when comparing the URANS, iLES-SOU, and iLES-BCD with LES-WALE showed consistently high differences in both, the auricle and atrium (absolute difference >20%). Given the high relative differences in the velocity fields, high spatial differences in WSS are expected. The surface average results, presented in Fig. 4(c), showed that at high SVC flow rates TAWSS is consistently underestimated by URANS and iLES-SOU. As can be seen in Fig. 2 and Fig. S2 available in the Supplemental Materials on the ASME Digital Collection these models failed to accurately describe the incoming jet of SVC blood which according to the 37 M result should remain closer to the atrium wall, producing higher shear rates. The mischaracterization of this SVC jet is also the likely reason for the very high auricle TAWSS difference reported for the URANS model at the 60/40 flow condition. At this flow rate, the SVC jet remained more focused compared to all other models. Further to this, high relative differences between LES and URANS results may be related to high instability associated with a transition to turbulence, a known weakness of this type of modeling (see Turbulence in the Right Atrium section). These results suggest that the convection scheme applied clearly impacts the TAWSS calculation and to a lesser degree so too does the inclusion of a WALE model.

### Mixing of Flows.

The data from this study provides some insight into how well each model characterizes the mixing of SVC and IVC flows. The relative differences (compared to the LES-WALE) in the IVC/SVC scalar calculations were high when SVC flow was high and decreased as the proportion of IVC flow increased. For IVC flow rates >55%, what might be considered physiological conditions for a healthy adult [24], the inclusion of the WALE SGS model or the convection scheme used had little impact (<5%). Perhaps the most noticeable difference was in the flow separation region, where the SVC inflow deviates from the atrial wall (indicated in red in Fig. 5(a)). Here the calculation of the IVC/SVC scalar with the URANS and iLES-SOU models showed very high relative differences compared to LES-WALE (>50%) with an underestimation of the proportion of SVC flow. The energy spectra in this region (Fig. 6, point 5) suggests dissipation through turbulent eddies in the LES-WALE model whilst the SOU model remains in the inertial subrange, causing a lack of mixing. Whilst the impact of the WALE model in this region may account for some of these differences (Fig. 7), results from the refined meshes suggest that the SOU convection scheme does not capture the recirculation zone accurately. The limitations of URANS in instances of flow separation are well-documented and are apparent here where the proportion of SVC blood is similarly underestimated [31]. These results underline the importance of the modeling approach in accurately characterizing the mixing of flows in the RA.

### Turbulence in the Right Atrium.

A key difference between URANS, implicit LES, and LES with WALE is how transitional flows are represented. Transitional flow has been reported in previous CFD simulations of vascular disease [10,32,33] and could be reasonably expected for the range of flow rates in this study. At sufficiently high resolution, as in DNS (but also LES), transitional flows are well represented without additional turbulence modeling. The 11 M cell iLES flow models in this study however are not sufficiently refined to capture all length scales but perform well in capturing the unsteady dynamics. URANS models have been developed for a wider range of engineering applications, where fully turbulent flow is more common. For this reason, they are not considered well suited to modeling cardiovascular flows [34]. Conversely, LES-WALE models are well suited to transitional flows. The lack of clear peaks in the energy spectra in Fig. 6 suggests a high degree of instability. This is particularly noticeable at points 2 and 5 in the shear layer where some periodic eddies might be expected in a fully turbulent flow. These instabilities become clear in the animation (Supplemental Video available in the Supplemental Materials on the ASME Digital Collection) of the 50/50 flow scenario. This data suggests the existence of transitional flow in the RA when SVC and IVC inlet flows are equal.

### Recommendations for Future Studies.

The higher density meshes in Fig. 2 provide a worthwhile point of comparison for the four models analyzed. Being far more refined, these models are no longer highly sensitive to the convection scheme used or the application of an SGS model and are approaching a DNS solution where all significant turbulence is resolved. These results showed that the URANS (Fig. S2 available in the Supplemental Materials on the ASME Digital Collection) and iLES-SOU failed to accurately characterize the incoming SVC flow. Using a BCD convection scheme better captured this flow on the coarser mesh and the addition of the WALE model provided some further improvement. Overall, these results suggest that for efficient and accurate simulation, both, the BCD convection scheme and an SGS model should be applied. Many commercial CFD solvers offer laminar flow models, essentially implicit LES with a first or second-order upwind convection scheme. Whilst these may provide accurate results on highly resolved meshes, a significant error is introduced when using a less dense mesh. URANS does not appear to be well suited to the characterization of flows in the RA.

### Limitations.

All simulations in this study implement a rigid-wall boundary as in previous CFD studies [4,5]. Though the superior and inferior vena cava are highly compliant and there is a contraction in the RA wall during the cardiac cycle, a comparison of rigid-wall CFD with 4D-MRI results has previously shown reasonable agreement [35]. IVC and SVC waveforms have some pulsatility, though not as much as in arteries. In this investigation, constant inlet flow rates were used. A comparison of constant and pulsatile inlet conditions shows little impact on time-averaged velocity results in the RA (Fig. S6 available in the Supplemental Materials on the ASME Digital Collection). Just as a BCD convection scheme produced more accurate results on the 11 M cell mesh compared to SOU, third-order convection schemes may be more suitable to capture complex flows developing in the RA.

## Conclusions

In this study, we apply four approaches to simulating the flow in the RA. The results show that general flow structures are appreciable using either URANS or iLES with SOU or BCD convection scheme. For high spatial accuracy in velocity, TAWSS, or IVC/SVC mixing, the WALE subgrid model and choice of convection scheme have a significant impact. Overall, an implicit LES approach with a BCD convection scheme and WALE SGS model provided the most accurate results. Velocity animations and the energy spectra suggest the existence of transitional flow in the RA.

## Funding Data

Partnership for Advanced Computing in Europe (PRACE) for awarding us access to the Supercomputer Mahti at the CSC-IT Center for Science, Finland (Grant No. 17DECI0055; Funder ID: 10.13039/501100001943).