Experimental and Computational Geometry
The idealised aneurysm geometry, as seen in figure 1, consisted of a 25.1 mm diameter inlet and outlet with the aneurysm extending to a maximum diameter of 50.2 mm. The distance from inlet of the aneurysm to its outlet was 75.3 mm. An axisymmetric aneurysm model was chosen for the experimental and computational analysis as it has a well-defined recirculation region within the aneurysmal sac, within which the species concentration varies slowly with time.
Experimental Fluid Flow Boundary Conditions
A constant volumetric flow rate of 3.94 × 10-6 m3/s was applied via a centrifugal pump (MCP-Z, Ismatec Inc.) resulting in a Reynolds number of 800. The resulting Peclet number for this system is 2,560,000 indicating strongly convection-dominated flow. A schematic of the experimental flow system is presented in figure 2.
Species of Interest
The food dye FD&C Blue #1 was used as the transported species. The diffusivity of the dye was calculated using the Stokes-Einstein equation. Assuming a spherical molecule, the characteristic radius of the molecule can be calculated using the molecular weight of the dye. The Stokes-Einstein diffusivity of a molecule of this size in water is 3.125 × 10-10 m2/s, resulting in a Schmidt number of 3200, which is in the range of blood borne species such as oxygen free radicals [13, 16].
Experimental Methodology: The Direct Spectrophometric Concentration Measurement Technique
The technique used in this analysis involves the extraction of small volumes of fluid from the test section using a syringe pump, figure 3, at numerous time-points over the course of the experiment. As the concentration field in the test section is time dependent until equilibrium is reached, the concentration values in each extracted volume would be expected to vary over the time course of the experiment. This technique was developed by Lutostansky (1996) [15] and has been presented in previous publications [17].
Withdrawal of the species sample from the geometry is achieved via a small diameter hole (<1 mm) in the experimental glass model of the idealised aneurysm. This hole was positioned at the location of maximum aneurysm diameter (figure 3). A small-bore syringe needle is placed flush to the inner surface of the glass model. By using a syringe pump (Harvard Apparatus, Holliston, MA, USA) it is possible to accurately extract small volumes of fluid from the test section. During the experiments carried out in this study, a volume of 0.8 ml was withdrawn over the course of one minute. The experiment was 7 minutes in duration and 0.8 ml of fluid was extracted every second minute starting at time equal to zero. Time zero was determined as the time at which the food dye entered the aneurysm section of the experimental model.
Experimental Measurement of Species Concentration
A spectrophotometer (SpectraMax Plus384, Molecular Device) was used to determine the concentration of food dye (FD&C Blue #1) in each sample volume at a wavelength of 630 nm. Both inlet concentration and zero concentration (water) were also determined and the results are presented in terms of a normalised concentration value.
Computational Fluid Dynamics
The computational model was built, subdivided and meshed in Gambit 2.1 (Fluent Inc. Lebanon, NH, USA). A structured quadrahedral grid was created for the idealised aneurysm model. The computational grid employed in this study has an equiangle skew value of less than 0.7 and an aspect ratio magnitude < = 5. Grid independence was established for Sherwood number, and error was specified to within ± 2%. The final grid size was identical for each method, consisting of 280,000 quadrahedral elements in the 2-D axisymmetric model.
The finite volume solver Fluent 6.1 was used to solve the two-dimensional, incompressible, Navier-Stokes equations. Discretisation of the governing momentum equations at each control volume was achieved by the Quadratic Upwind Interpolation for Convective Kinetics (QUICK) scheme, employing a three-point upstream-weighted quadratic interpolation for the cell face values. The Pressure Implicit with Splitting of Operators (PISO) pressure-velocity coupling algorithm was employed for this analysis, due to the time-dependent nature of the concentration profile. A convergence criterion of 1 × 10-3 was defined for the residuals of the continuity, momentum, and species transport equations [1, 18].
Convection-Diffusion Discretisation
The process of discretisation is presented here and is most easily illustrated by considering the steady-state conservation equation for the transport of a scalar quantity, ϕ:
(1)
where ρ equals the fluid density, is the velocity vector, Γϕ is the diffusion co-efficient of ϕ , and S
ϕ
is the source of ϕ per unit volume. This equation is applied to each control volume in the computational domain. Discretisation of equation 1 on any given cell yields:
(2)
and this is the general form of the equations solved by Fluent.
Discrete values of ϕ are stored at cell centres and are used to solve the diffusion term in equation 2. However, the convection term requires face values of the scalar quantity, ϕ
f
, which must be interpolated from the cell centre values. This is achieved by implementing an upwinding scheme. Upwinding involves deriving cell face values from quantities in the cell upstream, or "upwind", relative to the direction of the normal velocity, ν
n
, in equation 2. The upwinding schemes for solution discretisation available in Fluent 6.1 are First-Order Upwind, Second-Order Upwind, Power Law and QUICK. These four upwind schemes will now be introduced but further, more detailed, information is provided by Versteeg and Malalasekera, 1995 [19, 20].
First-Order Upwind Scheme
Quantities derived from this upwind scheme are determined by assuming that the cell-centre values of any field variable represent a cell average value and hold throughout the entire cell, resulting in face values that are identical to cell-centre values. Thus, by selecting this upwind scheme, the face value, ϕ
f
, is equal to the cell-centre value of, ϕ in the upstream cell.
Power Law Scheme
This discretisation scheme involves the interpolation of, ϕ
f
by using the exact solution to the one-dimensional convection-diffusion equation:
(3)
where Γ and ρu are constants across the interval ∂x. Integration of equation 3 yields the following equation that describes how ϕ varies with x:
(4)
where Pe is the Peclet number.
The variation of the scalar quantity between x = 0 and x = L is depicted in figure 4. This figure illustrates that for large Pe values (convection dominated flow), the value of ϕ at = L/2 is approximately equal to the upstream value, an assumption that is identical to the First-Order Upwind scheme described above.
Second-Order Upwind Scheme
When second-order accuracy is desired, a Taylor series expansion is used to calculate ϕ
f
via the following equation:
(5)
where ϕ and ∇ϕ are the upstream cell centred value and its gradient respectively, and is the displacement vector from the upstream cell centroid to the face centroid. The gradient of ϕ needs to be determined in order to solve equation 5. This is achieved using the divergence theorem, which in discrete form is written as
(6)
In equation 6, the face values are computed by averaging ϕ from the two cells adjacent to the face.
QUICK Scheme
The quadratic upstream interpolation for convective kinetics (QUICK) scheme is used to compute a higher-order value of the convected variable ϕ at a face and uses a three-point upstream-weighted quadratic interpolation.
For face e in figure 5 and assuming the flow is from left to right, the face value ϕ
e
can be written as:
(7)
The traditional QUICK scheme is obtained by setting Ψ = 1/8 in equation 7.
Boundary Conditions
Inlet Boundary Conditions
An identical volumetric flow rate to that applied in the experimental analysis was applied at the inlet of the computational model, resulting in a Reynolds number of 800 and a Peclet number of 2,560,000. Due to the time-dependent nature of the concentration profile until equilibrium is reached, an unsteady analysis was modelled implementing a time step size of 1 second. The bulk carrier fluid was assumed to be Newtonian, isothermal, and incompressible [21], with a density of 1000 kg/m3 and a dynamic viscosity of 0.001 Pas. A constant mass fraction of the species of interest was specified at the inlet.
Outlet Boundary Condition
An outflow boundary condition was applied to the outlet of the geometry, allowing the outflow values of pressure and velocity to be extrapolated from the interior of the computational mesh.
Symmetry Boundary Condition
By taking advantage of the symmetry of the experimental model, an axisymmetric boundary condition along the centre axis of the computational model was applied thus reducing computational effort.
Wall Boundary Condition
The model wall was assumed rigid with a no-slip boundary condition. To simulate the glass wall condition of the experimental model, the wall boundary was assumed to be an inert material with zero mass flux normal to the wall surface. Although this wall condition is not the physiologically consistent choice it can be considered physiologically relevant when the permeability of the arterial wall to a species is very low and the event is reaction limited.
Computational Concentration Measurements
The idealisation that the volume of fluid extracted from the experimental geometry represents a hemisphere was assumed in the computational analysis (figure 6(A)). A hemisphere with a volume of 0.8 ml has a radius of 3.37 mm. To obtain the computational results this hemisphere was subdivided into 5 separate volumes of equal height as shown in figure 6(B). The location of the centroid of these five volumes was determined and the concentration value at each of these points was obtained. These values were weighted in accordance to the corresponding volume of fluid in which each centroid was located. A volume weighted average value for concentration was subsequently obtained. This method allowed for a direct comparison to be made between the computationally obtained results to those obtained experimentally.
Four separate computational models were completed with the only variation between the models being the discretisation scheme applied to the solution of the species convection-diffusion equation. The results obtained were compared to determine whether or not the upwinding scheme influenced the computational results. Furthermore, in an attempt to determine the influence of the discretisation scheme selection on the accuracy of the MT solution and to validate the computational analysis of the species convection-diffusion equation, the experimental results were compared to the results obtained from a CFD analysis of the aneurysm geometry.