3D simulation of a viscous flow past a compliant model of arteriovenous-graft annastomosis

Hemodialysis is a common treatment for end-stage renal-disease patients to manage their renal failure while awaiting kidney transplant. Arteriovenous graft (AVG) is a major vascular access for hemodialysis but often fails due to the thrombosis near the vein-graft anastomosis. Almost all of the existing computational studies involving AVG assume that the vein and graft are rigid. As a first step to include vein/graft flexibility, we consider an ideal vein-AVG anastomosis model and apply the lattice Boltzmann-immersed boundary (LB-IB) framework for fluid-structure-interaction. The framework is extended to the case of non-uniform Lagrangian mesh for complex structure. After verification and validation of the numerical method and its implementation, many simulations are performed to simulate a viscous incompressible flow past the anastomosis model under pulsatile flow condition using various levels of vein elasticity. Our simulation results indicate that vein compliance may lessen flow disturbance and a more compliant vein experiences less wall shear stress (WSS). © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Arteriovenous graft for vascular access during dialysis has faced a long-lasting severe problem -blockage of the shunt due to thrombosis caused by occlusion of the outflow venous anastomosis and draining vein [1,2] . The intimal hyperplasia (IH), i.e., the thickening of the innermost layer of vessel wall, is the initial pathological event leading to thrombosis [3][4][5] . The most important factors for IH initialization include the significantly disturbed flows resulting from the incoming pulsatile arterial flow from the graft. Vascular endothelial cells lining the innermost vein wall are sensitive to the wall shear stress (WSS), wall normal stress (WNS), and their gradients (WSSG and WNSG) [6][7][8] . The abnormal values of these biomechanical variables caused by flow disturbances in the vein are thought to induce the initialization and development of IH in literature [9][10][11][12][13] .
Computational fluid dynamics has been used to investigate the flow and force fields in blood flow past the AVG anastomosis by many researchers. The influences of various factors, from AVG attaching angle and graft-to-vein diameter ratio, to non-Newtonian properties of blood and turbulent flows on the flow and force fields have been studied. Longest and Kleinstreuer [14] conducted computational comparison studies of two AVG configurations using user-enhanced finite volume code. Loth et al. [15] computationally investigated using spectral element method the transitional flow at the AVG anastomosis. They studies the distribution and importance of velocity fluctuation and WSS at two Reynolds numbers. Lee et al. [16] studied the influence of various flow division conditions between the two ends of the AVG-hosting vein at five different Reynolds numbers using the spectral element method. Ortega et al. [17] studied the reduction of WSS caused by the dialysis needle in the AVG by a polymer adapter. Manos et al. [18] combined computational fluid dynamics (ANSYS), in-vivo hemodynamics, histology data, and wall mechanics to investigate the flow patterns, WSS, and intimal hyperplasia of an arteriovenous shunt model. Kim et al. [19] computationally studied the flow patterns and WSS distributions in both distal and proximal AVG anastomoses using three different attaching angles. Khruasingkeaw et al. [20] studied the effect of five different attaching angles on the WSS distribution in the anastomosis of an arteriovenous graft using commercial software SolidWorks. The blood was assumed to be power-law fluid. Williams [21] performed computational fluid dynamics analysis for the optimization of AVG configurations to minimize shear stress and shear rate. McNally et al. [22] computationally characterized the flow and WSS fields for an arteriovenous graft anastomosis equipped with a novel modular anastomotic valve device (MAVD) using anastomotic angles of 90 °and 30 °.
However, we have noticed that all of the CFD studies of AVGrelated flows except McNally et al. [23] have assumed that blood vessel and AVG are rigid. McNally et al. [23] computationally (via ANSYS) investigated the reduction of flow disturbance and WSS by the MAVD-attached vein-graft anastomosis controlling on-and-off of the blood flow from the graft. They used fixed vein elasticity and did not study the effect of vein elasticity. However, vein and AVG are flexible and the vein flexibility may differ from patient to patient and from time to time for a specific patient. For example, the vein selected for dialysis may have different degree of flexibility for patients at different age. Furthermore, the vein flexibility may change over time due to IH development. Since forces exerted by the blood flow on the vein wall may depend on the deformability of the vessel [24][25][26] , and more and more simulations of blood flows are utilizing various fluid-structure-interaction methods [27][28][29][30] ; it is desired to gauge possible effect of vein-wall compliance on the flow and force field in AVG anastomosis. Decorato et al. [31] has investigated using the commercial software AN-SYS the differences between a rigid model and a deformable one for blood flow past an arteriovenous fistula (AVF) and found that the time-averaged WSS was overestimated by the rigid model. To the best of our knowledge, there exists no such studies for arteriovenous graft (AVG). As a first step to incorporate structure deformability, we introduce a 3D compliant model for the vein-AVG anastomosis and extend the lattice Boltzmann -immersed boundary (LB-IB) method [32][33][34][35] for fluid-structure interaction to study the influence of vein flexibility on flow and force fields [36] .
The LB-IB method is a lattice-Boltzmann based version of the popular immersed boundary (IB) framework [37,38] . The IB framework is well tested, efficient, and robust. It has long been applied to modeling and simulation of blood flows [39][40][41][42] . In the LB-IB method, the Navier-Stokes equations for motion of fluid and immersed structure are replaced by the lattice Boltzmann (LB) equations [43][44][45][46][47][48][49][50] . The LB equations are relatively simpler to use, easier to handle complex rigid flow domain, and more convenient to include additional fluid/flow physics into a new model [49,50] .
After verification and validation of the method and its implementation, a series of simulations are designed and performed to investigate influence of vein wall flexibility on flow and force fields in the model anastomosis of an arteriovenous graft. We consider five different levels of vein elasticity ranging from very flexible to rigid. To the best of our knowledge, such studies are not present in literature.
The remaining is organized as follows. Section 2 introduces a model problem abstracted from the arteriovenous graft for hemodialysis. Section 3 describes mathematical formulation of the model problem. Section 4 outlines numerical methods for the mathematical formulation. Section 5 performs verification and validation of the numerical methods and its implementation. Section 6 reports main computational results on flow and force fields. Section 7 concludes the paper with a summary.

Model problem
In this section, we introduce an AVG anastomosis model that incorporates the vessel and graft deformability. Modeling and simulation of blood flow in arteriovenous-graft anastomosis are challenging because: 1) The vein and AVG are thin-walled deformable structures; 2) the vein is living tissue made from nonlinearly viscoelastic material whose constitutive law is not yet well known; 3) the surrounding tissue of vein/graft is deformable porous media whose material properties not yet known well; 4) the blood is a multicomponent inhomogeneous non-Newtonian complex fluid; 5) the blood flow past AVG anastomosis could be laminar, transitional or turbulence; 6) the vein and AVG have complicated geometry and relative position and may not lie on the same plane.
To have a physiologically realistic yet mathematically tractable computational model is challenging. Therefore we focus on an idealized AVG-anastomosis model that incorporate the vein/graft elasticity. We treat the blood as a homogeneous incompressible viscous Newtonian fluid (the vein and AVG have relatively large diameters) and treat blood flow as laminar flow at intermediate Reynolds numbers. Important vein-graft configuration parameters include the diameter ratio of graft to vein and the graft attaching angle. In dialysis practice, these parameters are variable because the patient's vein selected for dialysis may have different size and shape. Besides, these parameters may change over time due to development of thrombosis and movement of the patient's forearm. Since influence of these variables have already been studied in literature (see Introduction) and our focus here is on the vein deformability, these parameters are fixed in our model (diameter ratio of 1:1 and attaching angle of π /3), to avoid possible complications. The Reynolds number in dialysis is relatively high (hundreds to thousands) [51] . But our LB-IB method is not stable for high Re flows. However, the Re may be much lower when the patient is at rest (e.g. sleep) during non-dialysis period. Therefore we fix the Reynold number to 100 in our simulations. The pulsatile flow velocity waveform at the graft inlet is taken from a recent study by Quanyu et al. [52] with the maximum speed downscaled to reflect dialysis practice. The steady native blood flow at the vein inlet is chosen to be the mean value of the blood flow at the graft inlet. The vein and graft are modelled as flexible structures with zero thickness and different elasticity by the immersed boundary framework. The vein wall is initially a circular straight cylinder and the graft a circular curved cylinder, both of constant diameter. Each structure (vein or graft) consists of two families of crosslinked flexible fibers that can be compressed/stretched or bent. The two sets of fibers are initially orthogonal to each other and they are mechanically homogeneous along the circumferential and longitudinal directions for each structure. The bending modulus and stretching coefficients of the constitutive fibers are related to the Young's modulus and Poisson ratio [53] . They are different for vein and graft and but are assumed to be the same in longitudinal and circumferential directions. The tissue around the vein/graft is modelled as elastic springs immersed in viscous fluid. Similar approach has been used in [27,30] . The mechanical properties of vein may vary from patient to patient and may change during the period of dialysis due to the development of intimal hyperplasia and formation of thrombosis. Therefore the vein elasticity is treated as a variable and five different levels of compliance are considered in our studies.
Although our AVG-anastomosis model is ideal, some of these assumptions are commonly used in existing studies. For example, blood flow is modelled as laminar flow in [14,19,54] and the vein/graft are circular cylinders on the same plane in [20,23] . The major advantage of our anastomosis model is the ability to simulate vein/graft deformability.
Specifically, we consider a periodic rectangular domain full of a viscous incompressible Newtonian fluid. A distal anastomosis consisting of a straight segment of vein and a curved segment of graft is placed in the center of the domain. The domain is affixed with a right-handed orthogonal coordinate system with the x-axis pointing to the opposite direction of the inflow at the vein inlet, the y-axis pointing from front to rear, and the z-axis pointing from bottom to top, as shown in Fig. 1 . The dimension of the domain is L x × W y × H z . The fluid domain is defined by the Cartesian space: and is associated with the Cartesian The graft (red) and vein (blue) are both modelled as flexible structures made of intersecting circumferential and longitudinal fibers and are defined as C and L represent the arc-length of circumferential fibers and longitudinal fibers, with subscript G for graft and V for vein. The surfaces of the structures are labled with curvilinear Lagrangian coor- on the surfaces. The position of the structure is described by X ( α, t ) ∈ .

Mathematical formulation
The LB-IB formulation for the model problem in dimensionless form is presented in this section.

Equations for fluid
We use Eq. (1) , the Bhatnagar-Gross-Krook (BGK) equation [55] , to describe the motion of fluid and the immersed structure (vein and graft): The function q ( x , v , t ) is the distribution function of singleparticle velocity v at time t and space x . Thereby, the expression represents the probability of finding a particle at moving with a velocity between v and v + d v [32] . The function is the famous BGK approximation [55] to the complex collision operator in the Boltzmann equation, where τ is the relaxation time connecting to the fluid kinematic viscosity ν via the relation: ν = 2 τ −1 6 [32] . In this model, ν is determined by the parameter Reynolds number ( Re ) through the relation: ν = u l ·L l Re , where u l is the characteristic velocity and L l is the characteristic length in lattice Boltzmann units.
To derive the fluid velocity ( u ), one first computes the fluid mass density ( ρ) and momentum ( ρu ) from the velocity distribution function q ( x , v , t ) according to ( 2 ) and ( 3 ).
Notice that it can be shown that our mathematical formulation for the fluid motion (i.e. Eqs. (1) through 3) is equivalent to the classic three-dimensional Navier-Stokes equations for a viscous incompressible Newtonian fluid.

Equations for structure
The elastic structure (vein or graft) is modelled by intersecting longitudinal and circumferential fibers with identical mechanical properties. We take the graft as an example and focus on the longitudinal fibers. In this case, α = (α 1 , α 2 ) ∈ G and α 1 is used to denote a longitudinal fiber, while α 2 is used to denote the arc-length along the fiber at its initial configuration. We use superscript l to identify the variables defined on longitudinal fibers.
We start by evaluating the deformation energy density ( E l ), which consist of a stretching/compression part ( E l s ) and a bending part ( E l b ). The constants K G s and K G b are the stretching/compression coefficient and bending rigidity of the graft respectively.
We can compute the stretching/compressing force density ( F l s ) and the bending force density ( F l b ) by taking the negative of the variational derivatives of E l s and E l b according to ( 6 ) and ( 7 ) respectively.
The circumferential fibers of the graft are treated similarly. We use superscript c to identify the variables defined on circumferential fibers. The stretching/compressing force density ( F c s ) and the bending force density ( F c b ) are computed similarly. Then the Lagrangian force density at any point on the graft can be computed by summing up all the four aforementioned components ( F l s , F l b , F c s and F c b ) as well as the contribution from the virtual springs ( F spring ). Namely, ∀ α ∈ G : where, Fibers modelling the vein are handled similarly and the interacting Lagrangian force densities on the fibers are computed the same way.

Equations for fluid-structure interaction
The fluid-structure interaction is treated by the immersed boundary method [56] . The first step is about the effect of structure on fluid. It describes how to impart the structure-fluidinteracting force to the fluid by computing the Eulerian force density f ( x , t ) from its Lagrangian equivalent ( 10 ): X ( α, t )) d α (10) The structure (or Lagrangian) force density F ( α, t ) is obtained via ( 8 ). The δ( x ) is the Dirac δ-function. The term x − X ( α, t ) is the difference of spatial coordinates between a Eulerian point and a Lagrangian point.
The second step is about the effect of fluid on the immersed structure. It describes how to interpolate the structure velocity U ( α, t ) from the fluid velocity u ( x , t ) as given below ( 11 ): The fluid velocity u ( x , t ) is obtained via ( 2 ) and ( 3 ). This is to state that the structure moves following the fluid.
Finally, the motion of the immersed structure is described by a first-order ODE system ( 12 ):

Initial and boundary conditions
In order to complete the mathematical model, initial and boundary conditions are also required. We denote the union of the left and right boundary (along x -axis) of the computational domain H Z ] by S 1 and denote the remaining boundary by S 2 . We assume the fluid will flow freely (friction free) on S 1 and get bounced back on S 2 : where q b is the velocity distribution function of particle on the channel surface after bounce back. Details are given in Section 4 .
We denote the inlets of the graft and vein by D G and D V , respectively. The velocity on D G is assumed to be pulsatile while that on the D V is assumed to be steady: where s a ( t ) is the profile of the pulsatile flow speed, s b is the constant speed of the steady flow, both are given in Section 6 . The inlets of the vein and graft and the outlet of the vein are fixed. The outer surfaces of the vein-graft are subject to forces exerted by surrounding tissue which are modelled by virtual spring immersed in a viscous fluid.
The initial condition for the fluid is that the fluid velocity is zero everywhere (except the inlets of vein and graft) in the computational domain: The vein-graft is initially located in the center of the computational domain.

Numerical method
We begin this chapter with introducing the construction of Lagrangian mesh for the vein-AVG anastomosis. Then we focus on the discretization of the equations.

Mesh generation for vein-AVG anastomosis
The surfaces of graft and vein are discretized by a Lagrangian mesh consisting of intersecting longitudinal and circumferential elastic fibers as shown in Fig. 2 . The mesh for the vein is in the shape of a cylinder surface while the mesh for the graft consists of 3 parts: G 1 is a straight tube model connecting the suture line (connecting the vein and graft) and the first circle of G 2 , which is a section of a torus, and G 3 is in the same shape of the vein. Both the G 3 and the vein are horizontal so that we can define the two inlets as disks ( D V , D G ) perpendicular to the x-axis. Discretization of G 1 , G 2 , and G 3 are similar to that of the vein. Discretization of the suture line (the intersection of the graft and the vein) takes some extra work. The approach we use is as follows: First project the first circle of G 2 on V along the attaching angle θ while maintaining x 2 ; then discretize the projection so that its nodes are much denser than the mesh points of the vein; finally, for each node on the projection, select a mesh point of the vein that is closest to it. Hence the suture line is discretized by all the selected mesh points. Thus mesh for G 1 is constructed. The mesh points of the vein within the suture line are effectively removed by annihilating the Lagrangian forces associated to those points.

Equation discretization
The BGK Eq. (1) is discretized by the lattice Boltzmann D3Q19 model [57,58] . At each node particles can move along 18 different directions and they may also be allowed to stay. Therefore the particle velocity space v can be discretized by a set of 19 velocities. Along each of the velocity direction, the discretizations in time and space are performed by the method of charactersitic lines. Due to cancellation of discretization errors, the seemingly first order scheme turns out to be second order in time and space.
where w k is the weight and c k is the speed of sound of the model. The external forcing term is treated by the scheme introduced by Guo et al. [59] . After advancing the lattice Boltzmann equations one step forward, the macroscopic variables like density ( ρ) and momentum ( ρu ) can be updated by discretizing ( 2 ) and ( 3 ) as follows:

Structure equations discretization
Standard finite difference method on a non-uniform mesh is used for all other differential equations. For example, the second derivative ( D 2 g )( α) is defined on nodes α as follows: here α n and X n represent the Lagrangian coordinate and Eulerian coordinates of the n th node, respectively, and n denote the initial arc-length between the node n and node n + 1 .
The two integrals ( 10 ) and ( 11 ) are discretized by trapezoidal rule. The Dirac-δ function is discretized in the standard way as in the IB method. The cosine function version is used.

Verification and validation
In this section, we address verification and validation of the numerical methods and their implementations. First we compare our results with existing results on simulation of a 3D viscous flow past a deformable sheet fixed at its mid-line. Then we perform mesh refinement studies on a 3D viscous flow through a deformable tube.
Zhu et al. [32] has studied a 3D viscous flow past a flexible sheet tethered at its mid-line. We have performed the simulation of this problem using our new code. We focus on the drag coefficient. The values of C d are exactly the same up to 8 digits in the first 68 steps. After 10 0,0 0 0 steps, the maximum error in C d is under 3%. Table 1 lists drag coefficient C d for two different cases. The   Fig. 3. Velocity profile on AVG inlet [52] . We perform mesh refinement studies on a 3D viscous flow through a deformable tube. The computational results from a series of simulations with different mesh-width ( h ) are used to estimate the order of accuracy r . The trend of convergence in r towards one is seen on a series of meshes within the capacity of the computer at our dispose (from 10 3 to 320 3 ). Theoretically, the order of accuracy in space and time should be one. This indicates the numerical methods are correct and the implementations are bug free.

Main computational results
In this section, we report the major computational results using varying vein elasticity and perform analyses through flow and force field visualization and comparison of the computational data. The focus is on the effect of vein deformability on flow and force fields in the vein-AVG system. In order to gauge the possible effect of vein flexibility, a series of simulations with a wide range of vein elasticity from very flexible (soft case) to very stiff (rigid case) are performed (notice that the elasticity of the graft is fixed). The model parameters adopted in our simulations are summarized in Table 2 . The case with the least stiffness is named the "soft" (or "deformable") case. The other four cases have a stiffness of 2, 4, 8, and 100 times of the soft case, and are named the "less soft", "less stiff", "stiff", and "rigid" case, respectively. Note that the vein and graft are still elastic in the stiff case. The pulsatile flow velocity profile s a ( t ) (see Fig. 3 ) at the graft inlet is taken from [52] . Note that the velocity is non-dimensionalized for lattice Boltzmann simulation. The steady native blood flow at the vein inlet is chosen to be the mean value of the blood flow at the graft inlet. All parameters on Table 2  and T c = 9300 . Some simulation results (e.g., soft and less soft cases, less stiff and stiff cases, respectively) are similar. Therefore, we focus on only the three typical cases (soft, stiff, and rigid).
Note that since our goal is to gauge the possible effect of vein/graft flexibility, rather than investigate the specific flow and force fields of blood flow inside a given patient-specific arteriovenous graft anastomosis, all model parameters and outcome Table 3 Wall shear stress (gradient) and wall normal stress gradients level.  quantities (e.g. WSS) are kept dimensionless throughout the paper. Given a corresponding real world problem, it is possible to convert the lattice Boltzmann results back to physical units via dimensionless ratios by way of similarity principle [47] . It is done for the WSS for the purpose of comparison with existing results. First, let us examine the averaged force field (WSS, WSSG, and WNSG) for all cases. We point out that the WSS values on the vein and graft computed from our simulations are in line with existing data in literature. The range of the averaged dimensional WSS (converted from the lattice Boltzmann units to physical units) is approximately from 0.35 Pa to 40.5 Pa as the vein flexibility decreases from the greatest (soft) to the least (rigid) in the range listed on Table 2 . The WSS of similar vessel-graft systems was found to vary from less than 1 Pascals to hundreds of Pascals [14,15,17,[19][20][21]23,60] . Table 3 shows the spatially averaged values of the WSS, WSSG, and WNSG on the vein and graft walls for five levels of the vein deformability. "V" denotes the vein and "G" denote the graft. One can see that as the vein becomes stiffer, all of these quantities increase in the vein. But they remain approximately the same in the graft except the rigid case where they increase sharply. Note that the graft elasticity is not changed except in the rigid case.
Our results indicate that a more complaint vein experiences less WSS on average. This is in qualitative agreement with the finding in [31] . Now let us report the soft, stiff, and rigid cases in two pairs for comparison. We compare the soft and stiff cases first. Fig. 4 plots contours of magnitude of velocity and vorticity on selected planes normal to the y-axis at a typical instant (instant a in Fig. 3 ). These figures show that the radial expansion of the vein is less and the flow speed in the vein is greater in the stiff case. Interestingly, more radial expansion in the graft is seen in the stiff case than in the soft case. This may be caused by the less radial expansion in the vein. The flow in the graft feels the "traffic" and slows down. Thus the increased pressure in the graft causes more radial expansion. Fig. 5 plots the WSS, WSSG, and WNSG on the vein/graft walls at instant a. These figures show that the distributions of these variables are quite different in the two cases and their values are greater in the stiff case. Now let us compare the two extreme cases: The soft case with the greatest elasticity and the rigid case with the least elasticity. Fig. 6 visualizes the contours of velocity magnitude on the plane normal to the y-axis at the four instants for rigid (left column) and soft (right column) cases. We notice that the velocity magnitude is smaller in the soft case than that in the rigid case. Also the contours in the soft case looks bigger due to the thicker diameters (caused by radial expansion) of the structures (vein and AVG). We notice that the expansion of the vein in the soft case is not uniform, the vein becomes curved and irregular, particularly in region close to the anastomosis. In both cases, the differences in the four time instants are observable. But the differences for the soft case are less distinct than the rigid case. Fig. 7 visualizes the contours of vorticity in magnitude on the plane normal to the yaxis at four instants for both rigid and deformable cases. Note that apparent differences between the two cases are seen at all four time instants. The contours of the soft case are more or less irregular (less smooth) than the rigid case. This is caused by the motions of the walls of both vein and graft. Now let us show the vorticity components. Our simulations show that the x-component of the vorticity is the least (very small in magnitude) compared to the other two components. Fig. 8 visualize the contours of ω y . It clearly shows the occurrence of secondary flows in the structure, specially in the vein. Similar secondary flows are seen in [61,62] .
From the above figures, we can see some differences among the four time instants in both soft and rigid cases. These differences reflect the influence of flow pulsatility. However, the influence of flow pulsatility on flow field is weaker in the soft case than in the rigid case: The flow patern in the vein of the soft case is quite similar at the four instants. Now we investigate force field by visualizing distributions of WSS, WSSG, and WNSG on the walls of vein and graft. Fig. 9 visualizes the WSS on the vein and AVG walls for both rigid (the first two rows) and soft (the third row) cases. From Fig. 9 , one can notice the obvious difference between the two cases. First of all, the vein and graft are deformed by the flows and become non-uniform in shape (e.g., the vein becomes somewhat curved). Secondly, the WSS on the wall of vein and graft is greater in the rigid case than in the soft case. Also while the differences of WSS among the four typical instants a, b, c, and d in the rigid case are discernible by the eye, they look quite similar in the soft case (therefore only one instant is shown here). The latter two differences are caused by the flexibility of the vein and graft. For a structure (vein and graft here) that can deform with flow, flow disturbances may be lessened to some degree because the structure moves with the flow (instead of completely standing against it). In another word, the structure flexibility can "absorb" impact of the flow. Results of WSSG and WNSG are similar. They are shown in the fourth (WSSG) and fifth (WNSG) rows in Fig. 9 for rigid (left) and soft (right) cases at instant a.
The averaged WSS, WSSG, and WNSG on the vein wall (labeled as "V") and graft wall (labeled as "G") are computed and listed on   Table 4 for both cases. We see slight variations (less than 15%) in the averaged values among the four instants.

Summary and conclusion
We have developed a 3D anastomosis model incorporating vein-graft deformability for blood flow past a distal arteriovenousgraft anastomosis based on the lattice Boltzmann-immersed boundary method using a non-uniform Lagrangian structure mesh for the vein-AVG anastomosis by elastic fibers. The method and its implementation are verified and validated by comparison with existing data and mesh refinement studies.
A series of simulations using five different levels of vein elasticity have been designed and performed to investigate the effect of vein elasticity on the flow and force fields. These flow and force fields have been visualized and analyzed. Our computational data are in agreement with relevant data in literature. The main simulation results are summarized as follows: 1) Vein elasticity may lessen the disturbance of flow pulsatility on the flow and force fields in the vein-graft anastomosis. 2) As the vein gets less elastic, the WSS, WSSG, WNSG and their spatially averaged values all increase in the vein.