Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine http://pih.sagepub.com/
A novel periodic boundary condition for computational hemodynamics studies Fereshteh Bahramian and Hadi Mohammadi Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine 2014 228: 643 originally published online 11 July 2014 DOI: 10.1177/0954411914542170 The online version of this article can be found at: http://pih.sagepub.com/content/228/7/643
Published by: http://www.sagepublications.com
On behalf of:
Institution of Mechanical Engineers
Additional services and information for Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine can be found at: Email Alerts: http://pih.sagepub.com/cgi/alerts Subscriptions: http://pih.sagepub.com/subscriptions Reprints: http://www.sagepub.com/journalsReprints.nav Permissions: http://www.sagepub.com/journalsPermissions.nav Citations: http://pih.sagepub.com/content/228/7/643.refs.html
>> Version of Record - Aug 8, 2014 OnlineFirst Version of Record - Jul 11, 2014 What is This?
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
Original Article
A novel periodic boundary condition for computational hemodynamics studies
Proc IMechE Part H: J Engineering in Medicine 2014, Vol. 228(7) 643–651 Ó IMechE 2014 Reprints and permissions: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/0954411914542170 pih.sagepub.com
Fereshteh Bahramian1 and Hadi Mohammadi1,2
Abstract In computational fluid dynamics models for hemodynamics applications, boundary conditions remain one of the major issues in obtaining accurate fluid flow predictions. For major cardiovascular models, the realistic boundary conditions are not available. In order to address this issue, the whole computational domain needs to be modeled, which is practically impossible. For simulating fully developed turbulent flows using the large eddy simulation and dynamic numerical solution methods, which are very popular in hemodynamics studies, periodic boundary conditions are suitable. This is mainly because the computational domain can be reduced considerably. In this study, a novel periodic boundary condition is proposed, which is based on mass flow condition. The proposed boundary condition is applied on a square duct for the sake of validation. The mass-based condition was shown to obtain the solution in 15% less time. As such, the mass-based condition has two decisive advantages: first, the solution for a given Reynolds number can be obtained in a single simulation because of the direct specification of the mass flow, and second, simulations can be made more quickly.
Keywords Finite volume method, computational fluid dynamics, vascular hemodynamics, turbulence, direct numerical solution, large eddy simulation, heart valve prostheses
Date received: 16 December 2013; accepted: 29 May 2014
Introduction Boundary conditions remain one of the major issues in obtaining accurate blood flow predictions using computational fluid dynamics (CFD) since it is impossible and most often unnecessary to simulate the whole domain. Instead, simulations are usually conducted on a subregion of interest having a certain boundary with the surrounding environment.1–4 Particularly in hemodynamics studies for cardiovascular applications, assigning appropriate boundary conditions is challenging in unsteady conditions.5–8 For simulating fully developed turbulent flows using large eddy simulation (LES) and direct numerical solution (DNS), periodic boundary conditions are particularly useful, in that the computational domain can be reduced considerably since even a generally unsteady turbulent flow can be approximated as repeating itself over a certain distance. The use of periodic conditions for such cases also eliminates the need for precise inlet conditions, which are often not known to the detail required for numerical simulation. Traditionally, periodic boundary conditions have been based on prescribing the pressure gradient;
however, several iterations at different imposed pressure gradients may be required before the desired Reynolds number is achieved. Also, considering the mean pressure gradient as a constant volumetric source results in variations of the mass flux with time. This is because of an imbalance between the total wall shear stress and the imposed pressure gradient at each time step.9 Also, the mass flow rate can instead be specified, allowing the pressure field to evolve as part of the solution field. This is convenient since only one simulation is required to obtain the solution for a particular Reynolds number. 1
Division of Mechanical Engineering, School of Engineering, Faculty of Applied Science, The University of British Columbia, Kelowna, BC, Canada 2 Biomedical Engineering Graduate Program, Faculty of Applied Science, The University of British Columbia, Vancouver, BC, Canada Corresponding author: Hadi Mohammadi, Division of Mechanical Engineering, School of Engineering, Faculty of Applied Science, The University of British Columbia, Okanagan, Kelowna, BC, V1V 1V7, Canada. Email: [email protected]
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
644
Proc IMechE Part H: J Engineering in Medicine 228(7)
Figure 1. The concept of periodic boundary condition for hemodynamics applications. (a) using conventional CFD approach and (b) using periodic boundary conditions.
Periodic boundary condition is significant for hemodynamics applications. Let us discuss a standard model, which is considered for the analysis of blood flow through an idealized stenosis. The geometry of interest is an idealized blood vessel with an axisymmetric stenosis. The degree of stenosis is 75% (area reduction). Let us assume that D is the diameter of the non-stenosed pipe. In a standard framework, the length of the computational domain is set to 24D, which is 5D for the inlet length, 2D for the length of the stenosed section and 17D for the outlet length. ‘‘Model A’’ in Figure 1(a) precisely demonstrates ‘‘the standard CFD model.’’ ‘‘Model B’’ in Figure 1(b) uses the same dimensions except for the outlet length, which is significantly reduced. This is because a small fraction of the outlet length is considered, and using the boundary conditions developed in this study, the same concepts are being applied in order to reach a particular condition in the outlet, for example, fully developed turbulent profile of velocity, but in a fraction of the central processing unit (CPU) time. In this work, we compare the performance of pressure- versus mass-based periodic boundary conditions for solving a fully developed turbulent duct flow using LES. The proposed boundary condition is applied on a square duct for the sake of validation (Figure 1). The two conditions are compared in terms of the accuracy of the solutions and the time required to obtain the solutions. Fully developed turbulent flow in a square duct was chosen as a test case because of the simplicity of the geometry and the availability of experimental and numerical data with which to compare. The Smagorinsky10 sub-grid scale (SGS) model with a damping function11 was used in the simulations. The remaining sections of this article outline the mathematical model, the numerical scheme and the results of the present comparisons.
Mathematical model The governing equations solved for LES are12 ∂ ui ∂ 1 ∂P ∂t ij ∂2 ui + +n ui uj = ∂xj r ∂xj ∂xj ∂t ∂xj ∂xj
ð1Þ
and ∂ ui =0 ∂xi
ð2Þ
where an over-bar indicates spatial filtering. In this study, the top hat filter has been used in three directions. This type of filter is usually employed in physical space. ui and P are the filtered velocity and pressure, respectively, and r and v are the density and kinematic viscosity, respectively. The effect of the sub-grid structures on the resolved scales is described by t ij = ui uj ui uj
ð3Þ
To close equations (1) and (2), this term must be modeled using resolved quantities. Since the goal of this study is a comparison between two different periodic boundary conditions, the Smagorinsky model has been chosen because of its ease of implementation. This model is based on the assumption of a balance between the energy production and dissipation in the kinetic energy equation for SGS. The SGS stresses are modeled using the Smagorinsky eddy viscosity model as given below12,13 t ij
dij t kk = 2nT Sij 3
where nT can be written as 2 þ3 1 ed263 S nT = Cs D S = 2Sij S 1=2 ij
ð4Þ
ð5Þ ð6Þ
here, Cs is the Smagorinsky constant, = (Dx Dy Dz)1=3 is the filter width and D dþ = 1=(xþ2 + yþ2 ) is the distance in wall units. As found by Piomelli et al.,14 the optimum value of Cs is around 0.1. Also, the above wall-damping function improves the near-wall asymptotic behavior of the SGS eddy viscosity. Several previous numerical simulations of turbulent flow in a square duct have been reported in the literature, demonstrating that the Smagorinsky model is suitable for this flow12,13,15 and in good agreement with DNS.16–18 There are two reasons for the success of the
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
Bahramian and Mohammadi
645
Smagorinsky model. First, it results in adequate diffusion and dissipation to produce stable numerical computations. Second, the low-moment statistics of the large eddies is not very sensitive to the SGS motion.19
Numerical scheme The two different conditions were implemented into a three-dimensional (3D), unstructured and parallel CFD code developed in this study. The code solves the incompressible Navier–Stokes equations using secondorder spatial discretization and a three-time-level method for time discretization. The equations are solved in a collocated, finite volume framework, and coupling between the mass and momentum equations was maintained using a variant of the Rhie and Chow method of interpolation. The diffusion term is discretized using a semi-implicit discrete formulation for the gradient operation. The convection term is implemented implicitly with an additional explicit term, which includes a flux limiter function. The fully implicit, nonlinear system of discretized equations was solved using biconjugate gradient stabilized (BICGSTB) method. The code was parallelized using the Portable, Extensible Toolkit for Scientific computation (PETSc) framework with a single-program multiple-data (SPMD) message-passing model and was run on the SHARCNET distributed memory parallel computer network.
Boundary condition The duct considered was D 3 D in cross section and 6.4D long and was discretized as by Ham et al.,12 using 48 volumes along the duct length and 64 3 64 finite volumes in cross section. Wall conditions were imposed on the four sides of the duct, while the inlet and outlet
were assumed to be periodic with one another. Simulations were run using the mass-based periodic condition and the more traditional pressure gradient– based condition. Using the mass-based condition, the velocity profile is obtained from the outlet of the domain, corrected to provide the desired mass flow rate and imposed at the duct inlet. Information for the pressure profile at the outlet is obtained from the mid-section, corrected to provide the desired outlet pressure force and then imposed at the outlet. The pressure field inside the domain and the pressure level at the inlet are obtained as part of the solution. Using the pressurebased condition, the total pressure is separated into mean and fluctuating parts. The mean part is considered as a constant source of momentum, and the fluctuating part is periodically applied at the streamwise periodic boundaries.
Results and discussion Because of the nature of the computational tests presented herein, the results set out to show two important facts that the results obtained using the two different boundary conditions are both accurate and similar. Both aspects are verified by comparing the present results to results obtained in previous studies using the same LES model. Before results of the simulations are presented, we briefly describe the measures taken to ensure that the present results are time-step and grid independent and numerically converged. Simulations of square duct flow were run for a Reynolds number (based on the bulk velocity and the side length) of 5810, the same as that used by Ham et al.12 as shown in Figure 2. To capture the effect of turbulent fluctuations in the duct, the non-dimensional time-step of 6.2 3 1023 was initially applied to the simulation based on the Re number and the length of
Figure 2. The square duct considered in this study with dimensions of 1.0 3 1.0 3 3.5 m3.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
646
Proc IMechE Part H: J Engineering in Medicine 228(7)
the geometry under study. As a next step, the smaller time-steps of Dt = 3.1 3 1023 and 1.55 3 1023 were considered to assess the time-step independence of the simulations. Table 1 summarizes the results of the time-step study and shows that the results for Dt = 3.1 3 1023 are slightly different from Dt = 6.2 3 1023 and in good agreement with those of Dt = 1.55 3 1023. Table 1 also shows that Dt = 3.1 3 1023 can provide results with high accuracy. As such, Dt = 3.1 3 1023 was chosen for all subsequent calculations, giving a maximum Courant–Friedrichs–Lewy (CFL) number of 3. Although this number is larger than the limit of most semi-implicit methods, it is deemed suitable for the fully implicit time discretization used herein. The sensitivity of the results to grid size was tested by doubling the grid resolution in both transverse and spanwise directions. Despite the increased grid resolution, the computed results shown in Figure 3 were very close to each other. It is worth mentioning that complete grid independence is difficult to achieve with LES since the range of the resolved scales is dependent on grid size. To remove this dependency, the mesh would have to be very fine, which is not worthwhile, since then LES becomes DNS. The accepted practice is to show that the average flow features like mean velocity and fluctuations intensity do not change significantly with the grid.20 Convergence of the CFD code was assessed on the basis of residuals normalized by the maximum values Table 1. Mean flow data at three different time-steps. Time-step
Ucl
Vmax
Wmax
1.55 3 1023 3.1 3 1023 6.2 3 1023
20.809 20.817 21.417
0.3810 0.3792 0.6460
0.4124 0.4186 0.9343
of the independent variable. To minimize the error in the results and to ensure that suitable convergence was obtained, the simulation was carried out at three different convergence criteria of 1024, 1025 and 1026. Table 2 compares several parameters of the flow such as maximum and average velocities in three directions and shows that there is no significant difference between the results obtained with maximum residuals of 1025 versus 1026, while the maximum residuals of 1024 cannot produce enough accurate results. This is also evident in Figure 4, which compares mean streamwise velocity and streamwise normal stress \ u2 . profiles for the different convergence criteria. Thus, the maximum normalized residuals of 1025 were chosen to all subsequent reported results. Unless otherwise stated, all the flow variables shown in the following plots are normalized with respect to the friction velocity. Initial conditions for the simulations were based on an analytical 1/7th power law velocity profile with a random 10% perturbation in all three directions of the velocity field. Simulations were run for several thousand time-steps until the initial transients diminished and a statistically stationary turbulent flow was obtained. The evolution of the flow was monitored by observing the instantaneous pressures and velocities at several points in the square duct domain; points were monitored near the walls and in the core of the flow. Figures 5 and 6 show the instantaneous velocity at a point close to the wall and demonstrate that a statistically stationary flow is obtained with the mass-based periodic condition after approximately 3000 time-steps, while use of the pressure gradient–based condition required nearly 5000 time-steps. Similar trends were observed at all other points monitored in the flow field. Once the statistically stationary result was obtained, simulations were continued for at least 5000 additional time-steps from which the time-average and turbulence statistics were obtained.
Figure 3. Comparison of (a) mean streamwise velocity and (b) streamwise normal stress \ u2 . for different convergence criteria.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
Bahramian and Mohammadi
647
Table 2. Mean flow data at four different flow field errors. Maximum iteration error
1024
1025
1026
Time (s) Ub Ucl Ucl/Ub Vmax Wmax
78,911 16.247 20.780 1.279 0.4183 0.4105
150,003 16.136 20.817 1.290 0.3792 0.4186
297,996 16.138 20.850 1.292 0.3810 0.4152
Figure 4. Comparison of (a) mean streamwise velocity and (b) streamwise normal stress \ u2 . at different grid resolutions.
Figure 5. Evolution of the instantaneous velocity field close to the wall using the pressure-based periodic condition.
Figure 6. Evolution of the instantaneous velocity field close to the wall using the mass-based periodic condition.
As Table 3 shows, in terms of the total number of solver iterations, the mass-based condition required roughly 20% fewer iterations to produce solutions with 5000 time-steps of useful data for averaging. (Note that the number of iterations does not scale precisely with the number of time-steps because the internal solver iterations required for an iteration of the
numerical algorithm can be quite different.) Using eight processors on SHARCNET, a solution using the mass-based condition was obtained approximately 15% faster than the same solution obtained using the pressure-based condition. Thus, a significant timesaving was realized using the mass-based periodic condition.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
648
Proc IMechE Part H: J Engineering in Medicine 228(7)
Table 3. Comparison of timing parameters of two different types of boundary conditions. Boundary condition type
Iteration per time-step
Total iterations
CPU time
Mass-based Pressure-based
5 7
28,000 34,000
42 h 48 h
CPU: central processing unit.
Figure 7. (a) Instantaneous secondary velocity vectors and (b) contours of instantaneous velocity on a transverse plane.
Figure 8. (a) Mean secondary velocity vectors and (b) contours of mean streamwise velocity in a transverse plane.
Instantaneous velocity vectors and contours are shown for one instant in time in Figure 7 for a transverse plane located in the middle of the duct. The fluctuating component of velocity was obtained by subtracting the computed mean velocity from the
resolved instantaneous velocity at each control volume. The ensemble-averaged mean velocity contours and the secondary velocity vectors obtained using the massbased condition and extracted from a cross section near the middle of the square duct are shown in Figure 8.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
Bahramian and Mohammadi
649
This plot clearly shows the characteristic corner vortices. The level of symmetry is as good or better than most results reported elsewhere,15,16 which indicates that sufficient time-steps have been used to obtain the statistical average. Figure 9 shows that axial velocities along the lower wall bisector (Y = 0), obtained using the two different methods, are effectively identical away from the wall and match identically with those obtained by Madabhushi and Vanka.15 While the result of massbased method is identical to Vanka’s axial velocity through the whole spanwise domain, the deviation between the pressure-based method and Vanka’s result
Figure 9. Profiles of mean streamwise velocity normalized with bulk velocity.
is about 1% close to the wall. It is worth noting that the accuracy of Madabhushi and Vanka’s15 results has been verified by DNS study of Gavrilakis et al.16 The distribution of secondary velocities (V) at three different locations of Y along the wall bisector is shown in Figure 10. Again, the plot indicates that solutions obtained using the two periodic boundary conditions are essentially the same. The quantitative agreement between two periodic boundary conditions used in this study falls within the maximum error of 1% away from the wall, while close to the wall, both method predict almost the same results. As seen in Figure 9, the computed results of two methods are slightly different from the results obtained by Madabhushi and Vanka15 everywhere. The maximum error (;8%) happens away from the wall at Y = 20.25 and 20.2 \ Z \ 20.02. Figures 11 and 12 show that the root-mean-square (RMS) intensities of the large-scale fluctuations in streamwise, spanwise and normal directions are effectively independent of the type of periodic boundary conditions and agree well with data obtained by Madabhushi and Vanka,15 with no loss of accuracy in the near-wall region. In Figure 11(a), the streamwise RMS velocities of two compared methods are in a very good agreement with each other and also Vanka’s data everywhere except at Y \ 20.25, where the results of pressure-based and mass-based periodic boundary condition methods are overestimated and underestimated, respectively, compared to Vanka’s data but still remain within an acceptable error of 1%. The spanwise and normal RMS velocities shown in Figure 11(b) and (c) appear to be a little higher away from the wall with maximum error of 5%. However, around the
Figure 10. Profiles of V at different Y locations.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
650
Proc IMechE Part H: J Engineering in Medicine 228(7)
Figure 11. Profiles of (a) turbulence normal stress \ u2 . , (b) normal stress in the Y direction \ v2 . , (c) normal stress in the Z direction \ w2 . and (d) turbulence kinetic energy along the lower wall bisector (Y = 0.5). TKE: turbulence kinetic energy.
maximum point, results obtained from both methods are identical to Vanka’s data. Figure 11(d) shows turbulence kinetic energies (TKEs) derived from both simulations. Here, a similar behavior to streamwise RMS velocity (Figure 11(a)) is found for TKE. The distributions of viscous, SGS and resolved turbulent shear stresses are shown in Figure 12. Comparing the magnitude of the three components of shear stresses on the mass-based method with the corresponding results from the pressure-based method and those of Madabhushi and Vanka,15 it is obvious that despite differences within 7%, all terms of shear stress are in a good agreement.
Conclusion Two different periodic boundary conditions have been used to solve for fully developed turbulent flow in a square cross-sectional duct using LES. The mass-based
Figure 12. Profiles of resolved, viscous and sub-grid scale (SGS) shear stress t zx along the lower wall bisector (Y = 0.5).
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
Bahramian and Mohammadi
651
condition is compared to the more traditional approach based on the specification of a constant streamwise pressure gradient. The results obtained from the two conditions are shown to be the same in terms of both mean flow and turbulence statistics. The mass-based condition was shown to obtain the solution in 15% less time. As such, the mass-based condition has two decisive advantages: first, the solution for a given Reynolds number can be obtained in a single simulation because of the direct specification of the mass flow, and second, simulations can be made more quickly.
Acknowledgements This study received the computing support from SHARCNET. Declaration of conflicting interests The authors declare that there is no conflict of interest. Funding This study was funded by the University of British Columbia (Start-up Grant) and NSERC (Discovery Grant). References 1. Park S, Lee SW, Lim OK, et al. Computational modeling with fluid-structure interaction of the severe m1 stenosis before and after stenting. Neurointervention 2013; 8: 23–28. 2. Cebral JR, Mut F, Weir J, et al. Association of hemodynamic characteristics and cerebral aneurysm rupture. AJNR Am J Neuroradiol 2011; 32: 264–270. 3. Cebral JR, Mut F, Raschi M, et al. Aneurysm rupture following treatment with flow-diverting stents: computational hemodynamics analysis of treatment. AJNR Am J Neuroradiol 2011; 32: 27–33. 4. Younis HF, Kaazempur-Mofrad MR, Chan RC, et al. Hemodynamics and wall mechanics in human carotid bifurcation and its consequences for atherogenesis: investigation of inter-individual variation. Biomech Model Mechanobiol 2004; 3: 17–32. 5. Mohammadi H and Mequanint K. Prosthetic aortic heart valves: modeling and design. Med Eng Phys 2011; 33: 131–147.
6. Mohammadi H and Bahramian F. Boundary conditions in simulation of stenosed coronary arteries. Cardiovasc Eng 2009; 9(3): 83–91. 7. Mohammadi H, Bahramian F and Wan WK. Advanced modeling strategy for the analysis of heart valve leaflet tissue mechanics using high-order finite element method. Med Eng Phys 2009; 31(9): 1110–1117. 8. Mohammadi H, Ahmadian MT and Wan WK. Timedependent analysis of leaflets in mechanical aortic bileaflet heart valves in closing phase using the finite strip method. Med Eng Phys 2006; 28: 122–133. 9. Moser RD, Kim J and Mansour NN. Direct numerical simulation of turbulent channel flow up to Ret = 590. Phys Fluids 1999; 11(4): 943–945. 10. Smagorinsky J. General circulation experiments with the primitive equations. Mon Weather Rev 1963; 91: 99–164. 11. Piomelli U, Ferziger JH and Moin P. Models for large eddy simulations of turbulent channel flows including transpiration. Report TF32, 1987. Stanford, CA: Department of Mechanical Engineering, Stanford University. 12. Ham FE, Lien FS and Strong AB. Large-eddy simulation of the square duct flow using a parallel PC cluster. In: Proceedings of the 8th annual conference of the CFD Society of Canada, Montre´al, Canada, 11–13 June 2000. 13. Xu HY and Pollard A. Simulation of square and annular duct turbulent flows. In: Proceedings of the 5th annual conference of the CFD Society of Canada, Victoria, Canada, 25–27 May 1997. 14. Piomelli U, Moin P and Ferziger JH. Model consistency in large eddy simulation of turbulent channel flows. Phys Fluids 1988; 31(7): 1884–1891. 15. Madabhushi RK and Vanka SP. Large eddy simulation of turbulence-driven secondary flow in a square. Phys Fluids A: Fluid 1991; 3(11): 2734–2745. 16. Gavrilakis S. Numerical simulations of low-Reynoldsnumber turbulent flow through a straight square duct. J Fluid Mech 1992; 244: 101–129. 17. Huser A and Biringen S. Direct numerical simulations of turbulent flow in a square duct. J Fluid Mech 1993; 57: 65–95. 18. Breuer M and Rodi W. Large-eddy simulation of turbulent flow through straight and curved ducts. ERCOFTAC Bull 1994; 22: 26–29. 19. Wilcox DC. Turbulence modeling for CFD. La Can˜ada Flintridge, CA: DCW Industries, Inc., 1998. 20. Mittal R, Simmons SP and Najjar F. Numerical study of pulsatile flow in a constricted channel. J Fluid Mech 1991; 3(11): 2734–2745.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
A novel periodic boundary condition for computational hemodynamics studies Fereshteh Bahramian and Hadi Mohammadi Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine 2014 228: 643 originally published online 11 July 2014 DOI: 10.1177/0954411914542170 The online version of this article can be found at: http://pih.sagepub.com/content/228/7/643
Published by: http://www.sagepublications.com
On behalf of:
Institution of Mechanical Engineers
Additional services and information for Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine can be found at: Email Alerts: http://pih.sagepub.com/cgi/alerts Subscriptions: http://pih.sagepub.com/subscriptions Reprints: http://www.sagepub.com/journalsReprints.nav Permissions: http://www.sagepub.com/journalsPermissions.nav Citations: http://pih.sagepub.com/content/228/7/643.refs.html
>> Version of Record - Aug 8, 2014 OnlineFirst Version of Record - Jul 11, 2014 What is This?
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
Original Article
A novel periodic boundary condition for computational hemodynamics studies
Proc IMechE Part H: J Engineering in Medicine 2014, Vol. 228(7) 643–651 Ó IMechE 2014 Reprints and permissions: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/0954411914542170 pih.sagepub.com
Fereshteh Bahramian1 and Hadi Mohammadi1,2
Abstract In computational fluid dynamics models for hemodynamics applications, boundary conditions remain one of the major issues in obtaining accurate fluid flow predictions. For major cardiovascular models, the realistic boundary conditions are not available. In order to address this issue, the whole computational domain needs to be modeled, which is practically impossible. For simulating fully developed turbulent flows using the large eddy simulation and dynamic numerical solution methods, which are very popular in hemodynamics studies, periodic boundary conditions are suitable. This is mainly because the computational domain can be reduced considerably. In this study, a novel periodic boundary condition is proposed, which is based on mass flow condition. The proposed boundary condition is applied on a square duct for the sake of validation. The mass-based condition was shown to obtain the solution in 15% less time. As such, the mass-based condition has two decisive advantages: first, the solution for a given Reynolds number can be obtained in a single simulation because of the direct specification of the mass flow, and second, simulations can be made more quickly.
Keywords Finite volume method, computational fluid dynamics, vascular hemodynamics, turbulence, direct numerical solution, large eddy simulation, heart valve prostheses
Date received: 16 December 2013; accepted: 29 May 2014
Introduction Boundary conditions remain one of the major issues in obtaining accurate blood flow predictions using computational fluid dynamics (CFD) since it is impossible and most often unnecessary to simulate the whole domain. Instead, simulations are usually conducted on a subregion of interest having a certain boundary with the surrounding environment.1–4 Particularly in hemodynamics studies for cardiovascular applications, assigning appropriate boundary conditions is challenging in unsteady conditions.5–8 For simulating fully developed turbulent flows using large eddy simulation (LES) and direct numerical solution (DNS), periodic boundary conditions are particularly useful, in that the computational domain can be reduced considerably since even a generally unsteady turbulent flow can be approximated as repeating itself over a certain distance. The use of periodic conditions for such cases also eliminates the need for precise inlet conditions, which are often not known to the detail required for numerical simulation. Traditionally, periodic boundary conditions have been based on prescribing the pressure gradient;
however, several iterations at different imposed pressure gradients may be required before the desired Reynolds number is achieved. Also, considering the mean pressure gradient as a constant volumetric source results in variations of the mass flux with time. This is because of an imbalance between the total wall shear stress and the imposed pressure gradient at each time step.9 Also, the mass flow rate can instead be specified, allowing the pressure field to evolve as part of the solution field. This is convenient since only one simulation is required to obtain the solution for a particular Reynolds number. 1
Division of Mechanical Engineering, School of Engineering, Faculty of Applied Science, The University of British Columbia, Kelowna, BC, Canada 2 Biomedical Engineering Graduate Program, Faculty of Applied Science, The University of British Columbia, Vancouver, BC, Canada Corresponding author: Hadi Mohammadi, Division of Mechanical Engineering, School of Engineering, Faculty of Applied Science, The University of British Columbia, Okanagan, Kelowna, BC, V1V 1V7, Canada. Email: [email protected]
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
644
Proc IMechE Part H: J Engineering in Medicine 228(7)
Figure 1. The concept of periodic boundary condition for hemodynamics applications. (a) using conventional CFD approach and (b) using periodic boundary conditions.
Periodic boundary condition is significant for hemodynamics applications. Let us discuss a standard model, which is considered for the analysis of blood flow through an idealized stenosis. The geometry of interest is an idealized blood vessel with an axisymmetric stenosis. The degree of stenosis is 75% (area reduction). Let us assume that D is the diameter of the non-stenosed pipe. In a standard framework, the length of the computational domain is set to 24D, which is 5D for the inlet length, 2D for the length of the stenosed section and 17D for the outlet length. ‘‘Model A’’ in Figure 1(a) precisely demonstrates ‘‘the standard CFD model.’’ ‘‘Model B’’ in Figure 1(b) uses the same dimensions except for the outlet length, which is significantly reduced. This is because a small fraction of the outlet length is considered, and using the boundary conditions developed in this study, the same concepts are being applied in order to reach a particular condition in the outlet, for example, fully developed turbulent profile of velocity, but in a fraction of the central processing unit (CPU) time. In this work, we compare the performance of pressure- versus mass-based periodic boundary conditions for solving a fully developed turbulent duct flow using LES. The proposed boundary condition is applied on a square duct for the sake of validation (Figure 1). The two conditions are compared in terms of the accuracy of the solutions and the time required to obtain the solutions. Fully developed turbulent flow in a square duct was chosen as a test case because of the simplicity of the geometry and the availability of experimental and numerical data with which to compare. The Smagorinsky10 sub-grid scale (SGS) model with a damping function11 was used in the simulations. The remaining sections of this article outline the mathematical model, the numerical scheme and the results of the present comparisons.
Mathematical model The governing equations solved for LES are12 ∂ ui ∂ 1 ∂P ∂t ij ∂2 ui + +n ui uj = ∂xj r ∂xj ∂xj ∂t ∂xj ∂xj
ð1Þ
and ∂ ui =0 ∂xi
ð2Þ
where an over-bar indicates spatial filtering. In this study, the top hat filter has been used in three directions. This type of filter is usually employed in physical space. ui and P are the filtered velocity and pressure, respectively, and r and v are the density and kinematic viscosity, respectively. The effect of the sub-grid structures on the resolved scales is described by t ij = ui uj ui uj
ð3Þ
To close equations (1) and (2), this term must be modeled using resolved quantities. Since the goal of this study is a comparison between two different periodic boundary conditions, the Smagorinsky model has been chosen because of its ease of implementation. This model is based on the assumption of a balance between the energy production and dissipation in the kinetic energy equation for SGS. The SGS stresses are modeled using the Smagorinsky eddy viscosity model as given below12,13 t ij
dij t kk = 2nT Sij 3
where nT can be written as 2 þ3 1 ed263 S nT = Cs D S = 2Sij S 1=2 ij
ð4Þ
ð5Þ ð6Þ
here, Cs is the Smagorinsky constant, = (Dx Dy Dz)1=3 is the filter width and D dþ = 1=(xþ2 + yþ2 ) is the distance in wall units. As found by Piomelli et al.,14 the optimum value of Cs is around 0.1. Also, the above wall-damping function improves the near-wall asymptotic behavior of the SGS eddy viscosity. Several previous numerical simulations of turbulent flow in a square duct have been reported in the literature, demonstrating that the Smagorinsky model is suitable for this flow12,13,15 and in good agreement with DNS.16–18 There are two reasons for the success of the
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
Bahramian and Mohammadi
645
Smagorinsky model. First, it results in adequate diffusion and dissipation to produce stable numerical computations. Second, the low-moment statistics of the large eddies is not very sensitive to the SGS motion.19
Numerical scheme The two different conditions were implemented into a three-dimensional (3D), unstructured and parallel CFD code developed in this study. The code solves the incompressible Navier–Stokes equations using secondorder spatial discretization and a three-time-level method for time discretization. The equations are solved in a collocated, finite volume framework, and coupling between the mass and momentum equations was maintained using a variant of the Rhie and Chow method of interpolation. The diffusion term is discretized using a semi-implicit discrete formulation for the gradient operation. The convection term is implemented implicitly with an additional explicit term, which includes a flux limiter function. The fully implicit, nonlinear system of discretized equations was solved using biconjugate gradient stabilized (BICGSTB) method. The code was parallelized using the Portable, Extensible Toolkit for Scientific computation (PETSc) framework with a single-program multiple-data (SPMD) message-passing model and was run on the SHARCNET distributed memory parallel computer network.
Boundary condition The duct considered was D 3 D in cross section and 6.4D long and was discretized as by Ham et al.,12 using 48 volumes along the duct length and 64 3 64 finite volumes in cross section. Wall conditions were imposed on the four sides of the duct, while the inlet and outlet
were assumed to be periodic with one another. Simulations were run using the mass-based periodic condition and the more traditional pressure gradient– based condition. Using the mass-based condition, the velocity profile is obtained from the outlet of the domain, corrected to provide the desired mass flow rate and imposed at the duct inlet. Information for the pressure profile at the outlet is obtained from the mid-section, corrected to provide the desired outlet pressure force and then imposed at the outlet. The pressure field inside the domain and the pressure level at the inlet are obtained as part of the solution. Using the pressurebased condition, the total pressure is separated into mean and fluctuating parts. The mean part is considered as a constant source of momentum, and the fluctuating part is periodically applied at the streamwise periodic boundaries.
Results and discussion Because of the nature of the computational tests presented herein, the results set out to show two important facts that the results obtained using the two different boundary conditions are both accurate and similar. Both aspects are verified by comparing the present results to results obtained in previous studies using the same LES model. Before results of the simulations are presented, we briefly describe the measures taken to ensure that the present results are time-step and grid independent and numerically converged. Simulations of square duct flow were run for a Reynolds number (based on the bulk velocity and the side length) of 5810, the same as that used by Ham et al.12 as shown in Figure 2. To capture the effect of turbulent fluctuations in the duct, the non-dimensional time-step of 6.2 3 1023 was initially applied to the simulation based on the Re number and the length of
Figure 2. The square duct considered in this study with dimensions of 1.0 3 1.0 3 3.5 m3.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
646
Proc IMechE Part H: J Engineering in Medicine 228(7)
the geometry under study. As a next step, the smaller time-steps of Dt = 3.1 3 1023 and 1.55 3 1023 were considered to assess the time-step independence of the simulations. Table 1 summarizes the results of the time-step study and shows that the results for Dt = 3.1 3 1023 are slightly different from Dt = 6.2 3 1023 and in good agreement with those of Dt = 1.55 3 1023. Table 1 also shows that Dt = 3.1 3 1023 can provide results with high accuracy. As such, Dt = 3.1 3 1023 was chosen for all subsequent calculations, giving a maximum Courant–Friedrichs–Lewy (CFL) number of 3. Although this number is larger than the limit of most semi-implicit methods, it is deemed suitable for the fully implicit time discretization used herein. The sensitivity of the results to grid size was tested by doubling the grid resolution in both transverse and spanwise directions. Despite the increased grid resolution, the computed results shown in Figure 3 were very close to each other. It is worth mentioning that complete grid independence is difficult to achieve with LES since the range of the resolved scales is dependent on grid size. To remove this dependency, the mesh would have to be very fine, which is not worthwhile, since then LES becomes DNS. The accepted practice is to show that the average flow features like mean velocity and fluctuations intensity do not change significantly with the grid.20 Convergence of the CFD code was assessed on the basis of residuals normalized by the maximum values Table 1. Mean flow data at three different time-steps. Time-step
Ucl
Vmax
Wmax
1.55 3 1023 3.1 3 1023 6.2 3 1023
20.809 20.817 21.417
0.3810 0.3792 0.6460
0.4124 0.4186 0.9343
of the independent variable. To minimize the error in the results and to ensure that suitable convergence was obtained, the simulation was carried out at three different convergence criteria of 1024, 1025 and 1026. Table 2 compares several parameters of the flow such as maximum and average velocities in three directions and shows that there is no significant difference between the results obtained with maximum residuals of 1025 versus 1026, while the maximum residuals of 1024 cannot produce enough accurate results. This is also evident in Figure 4, which compares mean streamwise velocity and streamwise normal stress \ u2 . profiles for the different convergence criteria. Thus, the maximum normalized residuals of 1025 were chosen to all subsequent reported results. Unless otherwise stated, all the flow variables shown in the following plots are normalized with respect to the friction velocity. Initial conditions for the simulations were based on an analytical 1/7th power law velocity profile with a random 10% perturbation in all three directions of the velocity field. Simulations were run for several thousand time-steps until the initial transients diminished and a statistically stationary turbulent flow was obtained. The evolution of the flow was monitored by observing the instantaneous pressures and velocities at several points in the square duct domain; points were monitored near the walls and in the core of the flow. Figures 5 and 6 show the instantaneous velocity at a point close to the wall and demonstrate that a statistically stationary flow is obtained with the mass-based periodic condition after approximately 3000 time-steps, while use of the pressure gradient–based condition required nearly 5000 time-steps. Similar trends were observed at all other points monitored in the flow field. Once the statistically stationary result was obtained, simulations were continued for at least 5000 additional time-steps from which the time-average and turbulence statistics were obtained.
Figure 3. Comparison of (a) mean streamwise velocity and (b) streamwise normal stress \ u2 . for different convergence criteria.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
Bahramian and Mohammadi
647
Table 2. Mean flow data at four different flow field errors. Maximum iteration error
1024
1025
1026
Time (s) Ub Ucl Ucl/Ub Vmax Wmax
78,911 16.247 20.780 1.279 0.4183 0.4105
150,003 16.136 20.817 1.290 0.3792 0.4186
297,996 16.138 20.850 1.292 0.3810 0.4152
Figure 4. Comparison of (a) mean streamwise velocity and (b) streamwise normal stress \ u2 . at different grid resolutions.
Figure 5. Evolution of the instantaneous velocity field close to the wall using the pressure-based periodic condition.
Figure 6. Evolution of the instantaneous velocity field close to the wall using the mass-based periodic condition.
As Table 3 shows, in terms of the total number of solver iterations, the mass-based condition required roughly 20% fewer iterations to produce solutions with 5000 time-steps of useful data for averaging. (Note that the number of iterations does not scale precisely with the number of time-steps because the internal solver iterations required for an iteration of the
numerical algorithm can be quite different.) Using eight processors on SHARCNET, a solution using the mass-based condition was obtained approximately 15% faster than the same solution obtained using the pressure-based condition. Thus, a significant timesaving was realized using the mass-based periodic condition.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
648
Proc IMechE Part H: J Engineering in Medicine 228(7)
Table 3. Comparison of timing parameters of two different types of boundary conditions. Boundary condition type
Iteration per time-step
Total iterations
CPU time
Mass-based Pressure-based
5 7
28,000 34,000
42 h 48 h
CPU: central processing unit.
Figure 7. (a) Instantaneous secondary velocity vectors and (b) contours of instantaneous velocity on a transverse plane.
Figure 8. (a) Mean secondary velocity vectors and (b) contours of mean streamwise velocity in a transverse plane.
Instantaneous velocity vectors and contours are shown for one instant in time in Figure 7 for a transverse plane located in the middle of the duct. The fluctuating component of velocity was obtained by subtracting the computed mean velocity from the
resolved instantaneous velocity at each control volume. The ensemble-averaged mean velocity contours and the secondary velocity vectors obtained using the massbased condition and extracted from a cross section near the middle of the square duct are shown in Figure 8.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
Bahramian and Mohammadi
649
This plot clearly shows the characteristic corner vortices. The level of symmetry is as good or better than most results reported elsewhere,15,16 which indicates that sufficient time-steps have been used to obtain the statistical average. Figure 9 shows that axial velocities along the lower wall bisector (Y = 0), obtained using the two different methods, are effectively identical away from the wall and match identically with those obtained by Madabhushi and Vanka.15 While the result of massbased method is identical to Vanka’s axial velocity through the whole spanwise domain, the deviation between the pressure-based method and Vanka’s result
Figure 9. Profiles of mean streamwise velocity normalized with bulk velocity.
is about 1% close to the wall. It is worth noting that the accuracy of Madabhushi and Vanka’s15 results has been verified by DNS study of Gavrilakis et al.16 The distribution of secondary velocities (V) at three different locations of Y along the wall bisector is shown in Figure 10. Again, the plot indicates that solutions obtained using the two periodic boundary conditions are essentially the same. The quantitative agreement between two periodic boundary conditions used in this study falls within the maximum error of 1% away from the wall, while close to the wall, both method predict almost the same results. As seen in Figure 9, the computed results of two methods are slightly different from the results obtained by Madabhushi and Vanka15 everywhere. The maximum error (;8%) happens away from the wall at Y = 20.25 and 20.2 \ Z \ 20.02. Figures 11 and 12 show that the root-mean-square (RMS) intensities of the large-scale fluctuations in streamwise, spanwise and normal directions are effectively independent of the type of periodic boundary conditions and agree well with data obtained by Madabhushi and Vanka,15 with no loss of accuracy in the near-wall region. In Figure 11(a), the streamwise RMS velocities of two compared methods are in a very good agreement with each other and also Vanka’s data everywhere except at Y \ 20.25, where the results of pressure-based and mass-based periodic boundary condition methods are overestimated and underestimated, respectively, compared to Vanka’s data but still remain within an acceptable error of 1%. The spanwise and normal RMS velocities shown in Figure 11(b) and (c) appear to be a little higher away from the wall with maximum error of 5%. However, around the
Figure 10. Profiles of V at different Y locations.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
650
Proc IMechE Part H: J Engineering in Medicine 228(7)
Figure 11. Profiles of (a) turbulence normal stress \ u2 . , (b) normal stress in the Y direction \ v2 . , (c) normal stress in the Z direction \ w2 . and (d) turbulence kinetic energy along the lower wall bisector (Y = 0.5). TKE: turbulence kinetic energy.
maximum point, results obtained from both methods are identical to Vanka’s data. Figure 11(d) shows turbulence kinetic energies (TKEs) derived from both simulations. Here, a similar behavior to streamwise RMS velocity (Figure 11(a)) is found for TKE. The distributions of viscous, SGS and resolved turbulent shear stresses are shown in Figure 12. Comparing the magnitude of the three components of shear stresses on the mass-based method with the corresponding results from the pressure-based method and those of Madabhushi and Vanka,15 it is obvious that despite differences within 7%, all terms of shear stress are in a good agreement.
Conclusion Two different periodic boundary conditions have been used to solve for fully developed turbulent flow in a square cross-sectional duct using LES. The mass-based
Figure 12. Profiles of resolved, viscous and sub-grid scale (SGS) shear stress t zx along the lower wall bisector (Y = 0.5).
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
Bahramian and Mohammadi
651
condition is compared to the more traditional approach based on the specification of a constant streamwise pressure gradient. The results obtained from the two conditions are shown to be the same in terms of both mean flow and turbulence statistics. The mass-based condition was shown to obtain the solution in 15% less time. As such, the mass-based condition has two decisive advantages: first, the solution for a given Reynolds number can be obtained in a single simulation because of the direct specification of the mass flow, and second, simulations can be made more quickly.
Acknowledgements This study received the computing support from SHARCNET. Declaration of conflicting interests The authors declare that there is no conflict of interest. Funding This study was funded by the University of British Columbia (Start-up Grant) and NSERC (Discovery Grant). References 1. Park S, Lee SW, Lim OK, et al. Computational modeling with fluid-structure interaction of the severe m1 stenosis before and after stenting. Neurointervention 2013; 8: 23–28. 2. Cebral JR, Mut F, Weir J, et al. Association of hemodynamic characteristics and cerebral aneurysm rupture. AJNR Am J Neuroradiol 2011; 32: 264–270. 3. Cebral JR, Mut F, Raschi M, et al. Aneurysm rupture following treatment with flow-diverting stents: computational hemodynamics analysis of treatment. AJNR Am J Neuroradiol 2011; 32: 27–33. 4. Younis HF, Kaazempur-Mofrad MR, Chan RC, et al. Hemodynamics and wall mechanics in human carotid bifurcation and its consequences for atherogenesis: investigation of inter-individual variation. Biomech Model Mechanobiol 2004; 3: 17–32. 5. Mohammadi H and Mequanint K. Prosthetic aortic heart valves: modeling and design. Med Eng Phys 2011; 33: 131–147.
6. Mohammadi H and Bahramian F. Boundary conditions in simulation of stenosed coronary arteries. Cardiovasc Eng 2009; 9(3): 83–91. 7. Mohammadi H, Bahramian F and Wan WK. Advanced modeling strategy for the analysis of heart valve leaflet tissue mechanics using high-order finite element method. Med Eng Phys 2009; 31(9): 1110–1117. 8. Mohammadi H, Ahmadian MT and Wan WK. Timedependent analysis of leaflets in mechanical aortic bileaflet heart valves in closing phase using the finite strip method. Med Eng Phys 2006; 28: 122–133. 9. Moser RD, Kim J and Mansour NN. Direct numerical simulation of turbulent channel flow up to Ret = 590. Phys Fluids 1999; 11(4): 943–945. 10. Smagorinsky J. General circulation experiments with the primitive equations. Mon Weather Rev 1963; 91: 99–164. 11. Piomelli U, Ferziger JH and Moin P. Models for large eddy simulations of turbulent channel flows including transpiration. Report TF32, 1987. Stanford, CA: Department of Mechanical Engineering, Stanford University. 12. Ham FE, Lien FS and Strong AB. Large-eddy simulation of the square duct flow using a parallel PC cluster. In: Proceedings of the 8th annual conference of the CFD Society of Canada, Montre´al, Canada, 11–13 June 2000. 13. Xu HY and Pollard A. Simulation of square and annular duct turbulent flows. In: Proceedings of the 5th annual conference of the CFD Society of Canada, Victoria, Canada, 25–27 May 1997. 14. Piomelli U, Moin P and Ferziger JH. Model consistency in large eddy simulation of turbulent channel flows. Phys Fluids 1988; 31(7): 1884–1891. 15. Madabhushi RK and Vanka SP. Large eddy simulation of turbulence-driven secondary flow in a square. Phys Fluids A: Fluid 1991; 3(11): 2734–2745. 16. Gavrilakis S. Numerical simulations of low-Reynoldsnumber turbulent flow through a straight square duct. J Fluid Mech 1992; 244: 101–129. 17. Huser A and Biringen S. Direct numerical simulations of turbulent flow in a square duct. J Fluid Mech 1993; 57: 65–95. 18. Breuer M and Rodi W. Large-eddy simulation of turbulent flow through straight and curved ducts. ERCOFTAC Bull 1994; 22: 26–29. 19. Wilcox DC. Turbulence modeling for CFD. La Can˜ada Flintridge, CA: DCW Industries, Inc., 1998. 20. Mittal R, Simmons SP and Najjar F. Numerical study of pulsatile flow in a constricted channel. J Fluid Mech 1991; 3(11): 2734–2745.
Downloaded from pih.sagepub.com at UNIV WASHINGTON LIBRARIES on October 28, 2014
