The flow of bile in the human cystic duct

ARTICLE IN PRESS

Journal of Biomechanics 37 (2004) 1913–1922

The flow of bile in the human cystic duct R.C. Ooia, X.Y. Luoa,*, S.B. China, A.G. Johnsonb, N.C. Birdb a

Department of Mechanical Engineering, The University of Sheffield, Mappin Street, Sheffield, S1 3JD, UK b Academic Surgical Unit, Royal Hallamshire Hospital, Sheffield, S10 2JF, UK Accepted 23 February 2004

Abstract Clinical studies suggest that the flow of bile in the biliary system may be a contributing factor in the pathogenesis of cholelithiasis, but little is known about its transport mechanism. This paper reports a numerical study of steady flow in human cystic duct models. In order to obtain parametric data on the effects of various anatomical features in the cystic duct, idealised models were constructed, first with staggered baffles in a channel to represent the valves of Heister and lumen. The qualitative consistency of these findings are validated by modelling two of the real cystic ducts obtained from operative cholangiograms. Three-dimensional (3D) models were also constructed to further verify the two-dimensional (2D) results. It was found that the most significant geometric parameter affecting resistance is the baffle clearance (lumen size), followed by the number of baffles (number of folds in the valves of Heister), whilst the least significant ones are the curvature of the cystic duct and the angle between the neck and the gallbladder. The study presented here forms part of a larger project to understand the functions of the human cystic duct, especially the influence of its various anatomical structures on the resistance to bile flow, so that it may aid the assessment of the risk of stone formation in the gallbladder. r 2004 Elsevier Ltd. All rights reserved. Keywords: Biliary system; Cystic duct; Bile flow; Gallstones; Computational Fluid Dynamics

1. Introduction The human biliary system consists of an organ and ductal system that creates, transports, stores, and releases bile into the duodenum to aid digestion of fats. The anatomy comprises the liver, gallbladder, and biliary tract (cystic, hepatic and common bile ducts), Fig. 1. The most common biliary diseases are cholelithiasis and cholecystitis. Cholecystectomy is the most commonly performed abdominal operation in the West, with some 60,000 operations estimated in the UK each year (Calvert et al., 2000) costing the National Health Service approximately d40 million per annum. The three essential factors in the pathophysiological genesis of cholelithiasis are believed to be: 1. super-saturation of bile with cholesterol, 2. presence of calculi nucleating agents, 3. reduction in gallbladder motility. *Corresponding author. Tel.: +44-114-2227752; fax: +44-1142227890. E-mail address: x.y.luo@sheffield.ac.uk (X.Y. Luo). 0021-9290/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2004.02.029

Holzbach et al. (1973) showed that supersaturated bile (Factors 1 and 2) frequently exist in healthy people and so research began to focus on the efficiency of the bile transport mechanisms (Factor 3). The gallbladder differs from other hollow organs in only having a single conduit (the cystic duct) for filling and emptying. The pathophysiological functions of cystic duct and its contribution to the genesis of cholelithiasis have been discussed extensively over the years (Lichtenstein and Ivy, 1937; Scott and Otto, 1979; Otto et al., 1979; Pitt et al., 1981,1982; Doty et al., 1983; and Sharp et al., 1990). Recent investigations have pointed out that cholelithiasis is related to abnormal cystic duct configurations (Bornman et al., 1988; Shaw et al., 1993; Kubota et al., 1993; Castelain et al., 1993; Caroli-Bosc et al., 1997; Deenitchin et al., 1998). In particular, Deenitchin et al. (1998) showed that, statistically, subjects with gallstones had longer or narrower cystic ducts than those without. Patients with Cystic Duct Syndrome (non-calculus partial obstruction) have also been found to have a low gallbladder ejection fraction

ARTICLE IN PRESS R.C. Ooi et al. / Journal of Biomechanics 37 (2004) 1913–1922

1914

Nomenclature a b c D Dn Cp L n P Q

baffle thickness (m) distance between consecutive baffles (m) baffle clearance (m) duct diameter or (m)  channel height  Repffiffiffiffiffiffiffiffiffiffiffi Dean number D=2r 2 pressure coefficient (DP/12rV2) (Euler Number) cystic duct length between the first and the last baffles (m) number of baffles static pressure (Pa) volumetric flow rate (m3 s1)

r R Rn Rd Re V DP g m y r co

radius of curvature (m) =DP/Q, cystic duct resistance (kg m4 s1) resistance for duct with n number of baffles (kg m4 s1) =(Rn/R0) dimensionless cystic duct resistance =(rVD/m) Reynolds number average fluid velocity at inlet (m s1) pressure drop (Pa) curvature of cystic duct (degree) dynamic viscosity (Pa s) angle of bend, cystic duct inlet (degree) fluid density (kg m3) the ratio of the flow rate in the major recirculating eddy to that of through-flow

acting as a sphincter. Bile flow through two- and threedimensional (2D and 3D) models of idealised and actual cystic ducts has been simulated to characterise the influence of anatomical structure on flow resistance. Further development of such numerical simulation models may also provide insight into the pathogenesis of gallstone formation and the origins of biliary pain.

2. The geometry of the cystic duct models

Fig 1. Operative cholangiograms: (i) Patient A—with gallstones; (ii) Patient B—without gallstone. (Image kindly provided by Department of Radiology, Royal Hallamshire Hospital, Sheffield.) (Magnification factor x1).

(Fink-Bennett et al., 1985). These studies suggest a link between complex cystic duct anatomy and prolonged retention of bile in gallbladder. It is generally accepted that prolonged stasis of bile in the gallbladder is a significant contributing factor to gallstone formation (Holzbach et al., 1973, Catnach et al., 1993). This suggests that the fluid mechanics, in particular the relationship between cystic duct geometries and resistance to bile flow, of the biliary system may play an important role in gallstone formation. As far as the authors are aware, there has been no reported numerical study on the relationship between flow, structure, the anatomy of the cystic duct and the corresponding resistance to flow from the gallbladder. The overall aim of the present research is to understand the role of fluid mechanics in the human biliary system. The present paper focuses on the role of the cystic duct

The anatomy of the human cystic duct is complex and varies between individuals. Apart from duct length and diameter, the presence of the ‘‘valves of Heister’’ in the lumen complicates the cystic duct geometry. The socalled valves consist of several semi-lunar or crescentshaped folds, Fig. 1. These folds were reported by Mentzer (1926) as ‘‘leaflets’’ and some of them appeared to be arranged in a spiral manner. The role of the valve, as to whether it acts as an active or passive impedance device, has been debated in the literature (Lichtenstein and Ivy, 1937; Scott and Otto, 1979; Otto et al., 1979). The cystic duct diameter ranges from 2 to 5 mm and its length from 10 to 60 mm. The number of folds in the cystic duct varies from between 2 and 14. To understand the importance of all the geometrical factors, the present study proceeded from simple, idealised 2D models first. 2.1. 2D idealised cystic duct models The idealised model employed here approximates the cystic duct to a 2D channel (the duct), with alternating baffles representing the folds (valves), and clearance, c, between a baffle and the wall representing the lumen, Fig. 2(i). Baffle clearance ratio, c/D, represents the lumen size relative to the channel height. Operative cholangiograms suggest that these to vary between 0.2 and 0.6. The geometry of the model follows the relationship:L ¼ bðn  1Þ þ an:

ARTICLE IN PRESS R.C. Ooi et al. / Journal of Biomechanics 37 (2004) 1913–1922

1915

2.2. 2D models of patients’ biliary system More realistic 2D models of the human biliary system were then constructed using intra-operative cholangiograms. Fig. 1(i) shows the cholangiogram of Patient A (with gallstones) whereas Fig. 1(ii) shows that for the ‘normal’ Patient B (without stones). The image was scanned to construct the geometry and the computational mesh of the flow domain, Figs. 2(iv) and (v). Bile flows from the gallbladder to the cystic duct and also from the liver to the hepatic duct. Both streams flow to the common bile duct. Table 2 lists the essential geometrical dimensions of the models for the two patients. 2.3. 3D cystic duct models

Fig 2. Computational domain of idealised cystic duct: (i) straight cystic duct model n ¼ 10 and c=D ¼ 0:5; (ii) curved duct n ¼ 10; c=D ¼ 0:3 and g ¼ 90 ; (iii) Duct with bent inlet y ¼ 180 and c=D ¼ 0:7; (iv) computational domain and mesh for Patient A; (v) computational domain and mesh for Patient B; (vi) 3D baffle model with semi-circular baffles; (vii) The 3-D spiral model n ¼ 10; c=D ¼ 0:5:

Since the duct length, L, is fixed in all models, increasing the number of baffles, n, is equivalent to reducing the distance between the baffles (the width of the fold in the valve), b. The effect of cystic duct curvature on flow structures was modelled on the same idealised ducts but with different curvature, g. In order to simulate the effects of the angle between the neck and body of the gallbladder a 180 bend at the inlet (to represent the neck) was incorporated into the model. Examples of straight, curved and bent inlet duct geometries are shown in Fig. 2, and Table 1 lists the geometric configurations of the 21 idealised cystic duct models simulated here. All the simulations were carried out with the duct length, L ¼ 50 mm; the duct diameter, D ¼ 5 mm; and the baffle thickness, a ¼ 1 mm: These data were estimated from a full-scale cholangiogram of Patient A in Fig. 1.

Resin models cast from freshly excised cystic ducts have been made to provide realistic 3D geometries. However, the model is subject dependent and extracting the relatively small dimensions with acceptable resolution proved extremely time consuming. Instead, qualitative understanding of 3D effects have been obtained from dimensional models representing two extreme conditions. The first comprises a straight tube with baffles that have been projected from the 2D idealised model, Fig. 2(vi). This results in a flow domain where fluid flows through a series of semi-circular baffles, with a clearance ratio of 0.5. This is the ‘3D baffle model’. Another 3D model was constructed with a spiral representing the valves of Heister, Fig. 2(vii). This ‘3D spiral model’, when projected onto a 2D plane is similar to the 2D model with a clearance ratio of 0.5, but with the baffle slightly inclined. Autopsy observations by Mentzer (1926) showed that the structure of the valves of Heister may be a combination of these two geometries.

3. The mathematical model and numerical methods 3.1. Assumptions Ultrasonographic imaging of the rate of change of the gallbladder volume (Dodds et al., 1985) suggests that human gallbladder empties at an average rate of 1 ml/ min. The average flow rate in the duct is about 0.5– 1.0 ml/min for fasting gallbladders and 2.0–3.0 ml/min after meal (Howard et al., 1991). The Reynolds number, Re, based on mean velocity and mean duct diameter, thus varies from 1 to 40. Therefore, laminar flow may be assumed. In addition, as the average time to empty a gallbladder with a mean volume of 35 ml of bile is about 30 min (Dodds et al., 1989), the flow is sufficiently slow to be considered as steady state.

ARTICLE IN PRESS R.C. Ooi et al. / Journal of Biomechanics 37 (2004) 1913–1922

1916

Table 1 Geometric parameters of the 2D idealised models Set

Description

I

Straight duct

II III IV V All D = 5 mm, a = 1 mm

n

c/D

L/D

g

y

No. of models

0.5

10

0

0

8

Varying c/D Curved duct

0, 2, 4, 6, 8, 10, 12, 14 10 10

0.3, 0.7, 0.9 0.5

10 10

0 0

3 4

Curved duct (varying c/D) Bent gallbladder neck (varying c/D)

10 10

0.7, 1.0 0.5, 0.7, 0.9, 1.0

10 10

0 30 , 45 , 60 , 90 90 0

0 180 Total

2 4 21

Table 2 Geometrical details of the two patients Patient

D (mm)

L/D

n

c/D

y

Boundary condition

A

Cystic duct Hepatic duct Common duct Gallbladder fundus

5.0 4.9 (inlet) 4.0 (outlet) 10.0

8.0 — — —

12 — — —

0.4–0.9 — — —

45 — — —

— Inlet velocity Pressure outlet Inlet velocity

B

Cystic duct (neck inlet) Hepatic duct

9.0 4.8 (inlet)

3.1 —

6 —

0.7–0.8 —

45 —

— Inlet velocity

Common duct Gallbladder fundus

2.5 (outlet) 23.6

— —

— —

— —

— —

Pressure outlet Inlet velocity

Ultrasonographic and cholangiographic images show that the cystic duct may move with respiration. However, in order to simplify the simulation at this initial stage of investigation, we assumed a rigid, no-slip wall for all the models. There have been few studies on the rheological properties of human bile. Gottschalk and Lochner (1990) reported, from 33 samples, that post-operative T-tube human bile is a Maxwell fluid. Coene et al. (1994) found that 11 of 36 hepatic bile samples displayed nonNewtonian behaviour. There appear to be no reports on the viscosity of gallbladder bile. Our preliminary measurements of fresh human bile collected at 37 C suggest that for a healthy person without gallstones, the gallbladder bile is a Newtonian fluid with a viscosity of 1 mPa s and a density of 1000 kg m3 (Luo et al., 2004). The governing equations for the model thus reduce to the familiar Navier–Stokes and the Continuity Equations for steady-state laminar flow with constant fluid property.

3.2. Boundary conditions Boundary conditions for the 2D and 3D models were assumed as follows. Duct lengths of five times the duct diameter were attached to the first and the last baffles. A uniform inlet velocity, obtained from volumetric flow

rate of the bile, was imposed on the flow inlet, and a zero static pressure condition was employed at the flow exit. The 2D models of patients’ biliary system included the gallbladder, which in vivo expels bile to the cystic duct by contraction. Since a rigid wall was assumed in the present simulations, a uniform velocity profile was imposed at the gallbladder fundus to simulate the effect of gallbladder contraction. A uniform velocity profile was also assumed at the inlet to the common hepatic duct. The mean inlet velocities at the gallbladder fundus and hepatic duct were obtained from clinically measured volumetric flow rates of 1 and 0.7 ml/min, respectively (Winkelstein and Aschner, 1926; Dodds et al., 1985). Details of the boundary conditions are given in Table 2. 3.3. Numerical methods Computational grids were generated using GAMBIT and the governing equations were solved with a proprietary computational fluid dynamics software, FLUENT5 (Fluent Inc., 1998). Numerical diffusion error was minimised by employing a second-order central differencing based upwind scheme. Convergence of the solution was judged by monitoring the scaled residuals and flow parameters at critical points as well as successively reducing the value of the convergence criterion. Convergence was achieved when the ratio of the total mass imbalance in all the cells to the mass flow

ARTICLE IN PRESS R.C. Ooi et al. / Journal of Biomechanics 37 (2004) 1913–1922

rate drops below 4–5 orders of magnitude from the initial estimation.

12

n = 14 Patient A

10

Patient B

3.4. Grid independence and validations

4. Results The characteristic parameter employed by previous workers in this field (Otto et al., 1979; Pitt et al., 1981,1982; Doty et al., 1983; Sharp et al., 1990) is the cystic duct resistance, R, defined as a ratio of pressure drop (DP) across the duct to its flow rate (Q), to describe the impedance exerted by the geometry on the flow. However, R is a dimensional quantity and interpretation of its relationship with the non-dimensional Reynolds number may prove somewhat difficult. Hence a nondimensional resistance, Rd=Rn/R0, where Rn is the resistance of a cystic duct model with ‘n’ baffles, and R0 is the corresponding resistance without baffles (n ¼ 0), is used in the present parametric study to reflect the influence of the cystic duct anatomy. 4.1. Effects of anatomical structure A parametric study involving 21 2D simulations (see Table 1) was carried out to investigate the effects of various anatomical features of the cystic duct on the flow conditions. Fig. 3(i) shows the typical effects of

12

Rd

8 6

10

4

8 6 4 2

2 0 0

10

20

30

40

Re

(i)

20

c/D = 0.3 Patient A Patient B

15

Rd

The resolution of the uniform-sided triangular (for 2D models), hexahedral (for the 3D baffle model) and tetrahedron (for the 3D spiral model) grids were determined by refining the cell size until grid independence was achieved. The optimum dimension of the side of the triangular grid was 0.25 mm, resulting in a minimum of 9000 control volumes for 2D models, 72,000 for 3D baffle model and 400,000 for 3D spiral model for grid-independent solutions. For similar flow conditions, the flow structures for the idealised cystic duct simulations were found to be qualitatively similar to those obtained by Berner et al. (1984) from flow visualisation in a rectangular duct with staggered baffles. Kelkar and Patankar (1987) obtained numerical solutions for a channel flow with staggered fins of negligible fin thickness and with stream-wise periodic boundaries at the inlet and outlet. For the flow conditions in the 2D idealised model described by c=D ¼ 0:5; b=D ¼ 1:09 and Re ¼ 50; the present work found a value of co (the ratio of the flow rate in the major recirculating eddy to that of through-flow) of 0.227, whereas Kelkar and Patankar (1987) reported for a similar conditions a value of 0.23. The difference may be attributed mainly to the finite width of the baffle employed here. Further details of these validations tests are given in the Appendix A.

1917

10

0.5 5 0.7 0.9

0

0

(ii)

10

20

30

40

Re

Fig 3. Dimensionless resistance vs. Reynolds number for the 2D ideal and patient models: (i) effects of n, c/D = 0.5; (ii) effects of c/D, n ¼ 10:

increasing the number of baffles or folds in a straight cystic duct on Rd. For a straight duct of a fixed length and baffle width, increasing the number of baffles is equivalent to reducing the baffle gap width. When the number of baffles is less than 10, Rd increases only slightly and linearly with increasing Re. At larger number of baffles, Rd increases rapidly and non-linearly with increasing Re. This effect is more pronounced at higher Re and larger number of baffles. At low Re, the resistance depends primarily on the number of baffles. The flow is able to negotiate the baffles with no or insignificant flow separation. At a given Re, the flow follows an increasingly more corrugated path with increasing number of baffles. Increasing Re causes flow separation behind each baffle, leading to higher losses. Thus duct resistance is expected to be higher in a duct with more baffles or small gaps between the baffles. Fig. 3(ii) suggests that clearance ratio (for a fixed number of baffles) has similar effect as the number of baffles on Rd in the straight cystic duct model. At a large clearance ratio of about c=DX0:7; Rd increases only slightly with increasing Re. The flow passage is relatively

ARTICLE IN PRESS 1918

R.C. Ooi et al. / Journal of Biomechanics 37 (2004) 1913–1922

Fig 5. Effects of the Reynolds number in straight cystic duct: (i) pressure along the duct centreline; (ii) stream function contours. Fig 4. Effects of the geometry on stream function: (i) number of baffles, c=D ¼ 0:5 and Re ¼ 40; (ii) baffle clearance, n ¼ 10 and Re ¼ 40:

unobstructed and the flow negotiates the baffles with minimal increase in resistance. At smaller clearance, Rd increases rapidly and non-linearly with increasing Re. At low clearance and high Re, the flow has to accelerate through the narrow gaps, thus creating large recirculations downstream of the baffles. Sections of streamlines for different clearance and number of baffles are shown in Fig. 4. The angle between the gallbladder and the cystic duct, referred to here as the neck angle, varies between individuals. Simulations were carried out with the most extreme case of a 180 neck angle with n ¼ 10; Re ¼ 40; and 1:0pc=Dp0:5: It was found that the neck exerts negligible additional resistance as compared to the straight duct, and for duct without baffles (n ¼ 0), the neck only adds 8% resistance to the duct. Fig. 5 shows the pressure distributions for c=D ¼ 0:5; n ¼ 10; and various values of Re: The pressure decreases linearly from the inlet until the flow encounters the first baffle. The baffle produces a narrowing of the flow passage, and continuity requires the flow to accelerate through the clearance. Thus, there is a corresponding sharp decrease in the pressure. After the baffle, a slight expansion follows where pressure recovers a little as the flow decelerates. This pattern is repeated over each baffle. The curvature of the human cystic duct also varies between individuals. An anatomical range of curvatures, from 30 to 90 , was simulated here, with the conditions defined in Sets III and IV in Table 1. The Dean number, Dn ¼ 0:5ReOðD=2rÞ; for these conditions ranges between 0.06 and 3.96 which are within the regime where secondary flow may not be important (Ito, 1969), so 2D simulation may be justified. It was found that the ratio of the resistance of the curved duct over the

straight one decreases with the increase of the baffle height from 27% (duct with no baffles) to 0.6% (where c=D ¼ 0:5; the most serious obstruction investigated), for Re ¼ 40; and n ¼ 10: This suggests that the energy dissipation due to baffles overrides that due to the curvature effect. 4.2. Models of patient’s biliary system Due to the rigid wall and steady flow assumptions, the present models do not simulate the contractile activity of gallbladder. The flow simulated here is similar to that of cholangiogram contrast medium which is usually injected through the fundus. Figs. 1, 2(iv) and (v) show that Patient A had a longer but more uniform cystic duct, with more baffles, and a smaller c/D than those in Patient B. The latter duct also narrows considerably between the gallbladder neck and the hepatic duct, making it more like a nozzle, a structure known to have small flow losses and hence resistance. The dimensionless resistance of these two patients’ biliary systems were computed over a range of Reynolds numbers and shown in Fig. 3 with those of the 2D idealised models. It is gratifying to note that the relationship between Rd and Re for these two patients follows the trend of the 2D idealised models, thus lending credence to the value of employing a simple model to establish parametric data. It is also seen that the resistance of Patient A is higher that that of Patient B for all values of Re computed. Therefore in order to maintain the same flow, Patient A would have required a greater pressure difference between the gallbladder and the common bile duct than that for Patient B. Since Patient A had gallstones but not Patient B, these results suggest a tentative relationship between bile flow resistance, anatomical structure of the cystic duct and the presence of gallstones.

ARTICLE IN PRESS R.C. Ooi et al. / Journal of Biomechanics 37 (2004) 1913–1922 12

2-D 3-D n=2 n = 10 n = 14

10

Rd

8 6 4 2 0 0

10

20

30

40

Re Fig 6. Comparison of the variation of Rd vs. Re between the 2D models and 3D for different numbers of baffles, c=D ¼ 0:5 for all cases.

6

1919

spiral folds (3D spiral model). The 2D and 3D-baffle models show that Rd increases more rapidly with increasing Re than that for the 3D spiral model and that the two-dimensional model has a higher resistance than both the 3D models. This is due to significant energy dissipation caused by flow separations behind the rectangular baffles, especially at high Re. At low Re, the resistance in the 3D spiral model is higher than that of the 3D baffle model. This is because although R0 is kept the same for both models, flow in the former has to negotiate a longer flow path through the spirals, thus its non-dimensional resistance is higher. However, as Re is increased, the resistance in the 3D baffle model increases more rapidly due to flow separations, which are absent in the 3D spiral model. It is worth noting that in dimensional terms, the maximum pressure difference required to deliver normal bile at a flow of 9 ml/min (Re ¼ 40) for 2D and 3D are about 7 and 5 Pa (both 3D baffle and 3D spiral models), respectively.

2-D 3-D Spiral

5

5. Discussion

3-D Baffle

Rd

All n = 10 4

3

2 0

10

20

30

40

Re Fig 7. Effects of the assumed geometry of the valve of Heister on Rd vs. Re.

4.3. 3D models The geometry of the semi-circular baffles in the 3D baffle model dictates that the flow is predominantly in the axial and radial directions, whilst that for the 3D spiral model would have the flow following the spiral direction and so the circumferential velocity component becomes important. Fig. 6 compares the dimensionless resistance of the 3D baffle model with that of the 2D models over the same range of Re and number of baffles. For a given number of baffles, the 3D model shows a smaller resistance than the 2D ones over the range of Re. The difference increases significantly with a larger number of baffles because the 2D baffles produce larger flow separation than 3D ones, due to their more squared shape. The effects of the assumptions made in the anatomy of the valves of Heister and the number of flow dimensions are examined in Fig. 7. It shows the resistances through a cystic duct where the valves of Heister were modelled as (i) 2D alternating baffles, (ii) 3D semi-circular baffles (3D baffle model) and (iii) 3D

The resistances of the biliary models from the two patients and also those of the 3D models fit within the parametric range of the idealised duct models, thus lending credence to employing a 2D model to study bile flow in cystic duct. The viscosity of the bile considered here is believed to be representative of an average person without gallstones. Hence the trends in the change in resistance to geometric factors may resemble those in the early stage of gallbladder diseases. The maximum pressure difference computed for our models is only about 7.7 Pa. This value appears underestimated, considering that pressure readings across the canine cystic duct have been estimated at 100–200 Pa (Winkelstein and Aschner, 1926). There may be several reasons for this: firstly, our experimental study (Luo et al., 2004) has revealed that, although the bile of a healthy person is about the same as that for water, the viscosity of bile from patients with chronic gallbladder disease can be hundred times higher than this. Therefore the resistance, or the pressure required to deliver the bile flow, in these patients will be significantly higher, since pressure increases linearly with viscosity. Secondly, the lumen of the cystic duct in some patients can be much smaller than that employed here. (We have used the upper limit, 5 mm.). Therefore, the measured pressure of a patient can easily be two orders of magnitude higher, since a 1% reduction in the duct diameter can increase the pressure drop by 4%. The canine bile viscosity (and cystic duct geometries) were not reported by Winkelstein and Aschner and so an estimate of their effects on pressure drop cannot be made. Thirdly, and most importantly, our rigid duct models are basically in the ‘open’ state and so do not require ‘opening’ or ‘inflating’

ARTICLE IN PRESS R.C. Ooi et al. / Journal of Biomechanics 37 (2004) 1913–1922

to allow the passage of bile. Similarly, no pressure is required to maintain the ‘open’ state. By contrast, in the whole gallbladder preparation, an ‘opening’ pressure is required to initiate bile flow but, once established, no significant increase in pressure is required to maintain flow (Bird, 1998). Therefore, in the absence of matched geometrical parameters between the patient anatomy and the CFD models, only qualitative data can be obtained from these rigid wall computations. Thus, subsequent studies will be quantitative pressure estimates using 3D models obtained from scans of resin casts of the cystic duct following operative removal. In addition, we will examine the fluid–structure interactions between the elastic wall of the duct and bile flow (Ooi, 2003). This parametric study on an idealised cystic duct model shows that the most significant geometric factor affecting the non-dimensional resistance is the duct clearance, followed by the number of baffles, whilst the least significant ones are curvature of the cystic duct and the angle of the neck. At large clearance (c=D > 0:7) or small number of baffles (no6) the resistance to flow increases almost linearly with increasing Re. Conversely, resistance increases rapidly and non-linearly with increasing Re. Increasing resistance implies increasing difficulty in evacuating bile from the gallbladder. The non-dimensional resistance reported in the parametric study is always larger than unity. The baffles in the cystic duct models create obstacles to the flow, hence a larger pressure gradient is required for a given flow than in that of a duct without baffles. This suggests that the folds projecting into the lumen of the cystic duct may function as a flow resistor (passively) or regulator (if any contractile activity takes place in the muscle layer of the duct). Dodds et al. (1989) suggested that it minimises the luminal pressure in the common bile duct thus avoiding excessive distension during the expulsion of gallbladder bile. The simulations here indicate that patients with a high resistance cystic duct (e.g. with small baffle clearances and large number of baffles) may face a greater risk of forming gallstones, by facilitating stasis of gallbladder bile. The resistance, R, used by previous workers and the non-dimensional one used here, Rd, do not reflect the influence of viscosity on the resistance of the cystic duct. In fact, Fig. 3 may give the misleading impression that increasing viscosity (decreasing Re) may lead to reduced resistance. Most engineering applications employ pressure coefficient or the Euler number, (Cp =DP/12rV2, the ratio of pressure difference to dynamic pressure) to show flow resistance. The higher the pressure coefficient, the greater the pressure difference for a given flow and hence the greater resistance. The data in Fig. 3(i) are recalculated in terms of Cp in Fig. 8 and show that increasing the viscosity increases the pressure coefficient and hence resistance.

104

103 n=14

Cp

1920

2

10

n=0

101

1 0

10

20

30

40

Re Fig 8. Pressure coefficient vs. Reynolds number; c=D ¼ 0:5:

It should be noted that the good agreement between the 2D and 3D models does not consider the effect of secondary flow induced by curvature, although this may be insignificant due to the low Dean Number employed. Further work is required to determine a more realistic 3D geometry of the biliary system, the effects of the elastic properties the cystic duct wall and rheological properties of gallbladder and hepatic bile, so that 3D unsteady flow in the human biliary system could be studied and the role of fluid mechanics further quantified. The simulations performed here were part of an effort to establish a non-invasive method for investigating the role of fluid mechanics in the pathogenesis of gallstone formation. They have demonstrated the value of CFD techniques as an additional tool in the study of the pathogenesis of gallstone formation.

6. Conclusion Numerical simulations of bile flow in cystic duct models have been conducted to investigate the functions of the cystic duct in gallstone formation. The pattern of cystic duct resistance with geometric configurations of the duct was analysed. The reliability of these simple models were confirmed by the simulations of two biliary systems with their geometry mapped from operative cholangiograms, as well as of the 3D models. It was found that the most significant geometric factor responsible for increasing cystic duct resistance is the baffle clearance, followed by the number of baffles, whilst the least important are the angles between the duct and gallbladder. Thus it may be hypothesised that patients with a small lumen and large number of valves of Heister are more likely to develop gallstones, since the higher resistance through the duct will promote stasis of bile in the gallbladder. Further work is required to measure the effects of elasticity of the cystic duct and possibly non-Newtonian behaviour of bile when stones

ARTICLE IN PRESS R.C. Ooi et al. / Journal of Biomechanics 37 (2004) 1913–1922

1921

0.01 A

B

0 0.047

0.049

0.051

0.053

-0.01

y = 1.5mm -0.02

mesh size =0.33 mesh size =0.25

-0.03

A

B

mesh size =0.20

-0.04

(a)

(b) 2.4

(Pa)

Static Pressure (Pa)

3.7

3.2

2.7 0.048

0.053

2.0 1.8 1.6

A

1.4 1.2 B

1.0 0.054

0.058

(c)

2.2

0.056

0.058

0.06

0.062

(d)

Fig 9. (a) Detailed flow pattern between two baffles; (b) y-velocity profile (m/s) at y ¼ 1:5 mm; (c) static pressure plot across A–B; (d) Axial position, x(m).

are formed before a full account of the role of fluid mechanics in gallstone formation can be made.

3.5E-03

Y-Velocity (m/s)

2.5E-03

Acknowledgements The authors would like to thank Swann–Morton Foundation for financially supporting RCO and the UK CVCP for providing funding through an ORS award for RCO.

1.5E-03

5.0E-04

-5.0E-04

0.08

1st Order Upwind 1st Order Power Law 2nd Order Upwind 2nd Order QUICK

0.09

axial position (m)

-1.5E-03

(a)

Appendix A. Grid independence, error estimation and validation check A.1. Grid independence, numerical diffusion, error estimation and validation check A.1.1. Grid independence check was carried for ducts with c/D=0.3, n=14 and Re=40 Fig. 9(a) shows the detailed flow pattern between two baffles. Effects of the grid size is shown by comparing the velocity and pressure along a line at y ¼ 1:5 mm from the duct wall. Fig. 9(b) shows the velocity and Fig 9(c) the static pressure along that line. The grid sizes (the side of the triangular cell) for the three grids were 0.33, 0.25 and 0.20 mm.

3.0E-02

Y-Velocity (m/s)

2.0E-02

1st Order Upwind 1st Order Power Law 2nd Order Upwind 2nd Order QUICK

1.0E-02

0.0E+00 0.067

0.068

axial position (m) 0.069 0.070

0.071

(b) Fig 10. (a) Y-velocity plot at duct centre line after the last baffle; (b) Yvelocity plot across A–B as shown in Fig. 9a (Re ¼ 40; n ¼ 10;c=D ¼ 0:5).

ARTICLE IN PRESS 1922

R.C. Ooi et al. / Journal of Biomechanics 37 (2004) 1913–1922

For 2D case both Fig 9(b) and 9(c) suggest that a grid size of 0.25 mm is sufficient to attain grid independent simulation. For three dimensional case, the pressure is more sensitive to grid resolution (compared to velocity). Fig. 9(d) also suggests that a grid size of 0.25 mm provides grid independence. A.1.2. Effects of discretisation schemes on numerical diffusion Comparison of the effects of different numerical schemes for 3-D baffle model is provided in Fig. 10. The comparison shows that the more sophisticated QUICK scheme has no significance discrepancies with second order Upwind scheme. Therefore second order Upwind scheme is used for the modelling.

References Berner, C., Durst, F., McEligot, D.M., 1984. Flow around baffles. Transactions of the ASME—Journal of Heat Transfer 106, 743–749. Bird, N.C., Higham, S.E., Ahmed, R., Majeed, A.W., 1998. A perfused, whole, human gallbladder model for measuring motility in vitro. Neurogastroenterology and Motility 10 (5), 454. Bornman, P.C., Kottler, R.E., Terblanche, J., Kingsnorth, A.N., Kringe, J.E.J., Marks, I.N., 1988. Does low entry of cystic duct predispose stones in the common bile duct? British Medical Journal 297, 31–32. Calvert, N.W., Troy, G.P., Johnson, A.G., 2000. Laparoscopic cholecystectomy: A good buy? A cost comparison with smallincision (mini) cholecystectomy. European Journal of Surgery 166, 782–786. Caroli-Bosc, F., Demarquay, J., Conio, M., Deveau, C., Hastier, P., Harris, A., Dumas, R., Delmont, J., 1997. Is biliary lithogenesis affect by length and implantation of cystic duct? Digestive Diseases and Sciences 42 (10), 2045–2051. Castelain, M., Grimaldi, C., Harris, A.G., Caroli-Bosc, F., Hastier, P., Dumas, R., Delmont, J., 1993. Relationship between cystic duct diameter and the presence of cholelithiasis. Digestive Diseases and Sciences 38 (12), 2220–2224. Catnach, S.M., Anderson, J.V., Trembath, R.C., et al., 1993. Effect of octreotide on gallstone prevalence and gallbladder motility in acromegaly. Gut 34, 270–273. Coene, P.P.L.O., Groen, A.K., Davids, P.H.P., Hardeman, M., Tytgat, G.N.J., Huibregtse, K., 1994. Bile viscosity in patients with biliary drainage. Scandinavian Journal of Gastroenterology 29, 757–763. Deenitchin, G.P., Yoshida, J., Chijiiwa, K., Tanaka, M., 1998. Complex cystic duct is associated with cholelithiasis. HPB Surgery 11, 33–37. Dodds, W.J., Groh, W.J., Darweesh, R., Lawson, T.L., Kishk, S., Kern, M.K., 1985. Sonographic measurement of gallbladder volume. American Journal of Roentgenology 145, 1009–1011. Dodds, W.J., Hogan, W.J., Geenen, J.E., 1989. Motility of the biliary system. In: Schultz, S.G. (Ed.), Handbook of Physiology: the

Gastrointestinal system, Vol 1, Section 6, Part 2 (28). American Physiological Society, Bethesda, MD, pp. 1055–1101. Doty, J.E., Pitt, H.A., Kuchenbecker, S.L., DenBesten, L., 1983. Impaired gallbladder emptying before gallstone formation in the prairie dog. Gastroenterology 85, 168–174. Fink-Bennett, D., DeRidder, P., Kolozsi, W., Gordon, R., Rapp, J., 1985. Cholecystokinin cholescintigraphic findings in the cystic duct syndrome. Journal of Nuclear Medicine 26, 1123–1128. Fluent Incorporated, 1998. Fluent5 User’s Guide, 1999. GAMBIT Modelling Guide, Lebanon NH. Gottschalk, M., Lochner, A., 1990. Behaviour of postoperative viscosity of bile fluid from T-drainage. A contribution to cholelithogenesis. Gastoenterologisches Journal 50 (2), 65–67. Holzbach, R.T., Marsh, M., Olszewski, M., Holan, K., 1973. Cholesterol solubility in bile. Evidence that supersaturated bile is frequent in healthy man. Journal of Clinical Investigation 52, 1467–1479. Howard, P.J., Murphy, G.M., Dowling, R.H., 1991. Gallbladder emptying patterns in response to a normal meal in healthy subjects and patients with gallstones: ultrasound study. Gut 32, 1406–1411. Ito, H., 1969. Laminar flow in curved pipes. Zeitschrift fur . Angewandte Mathematik und Mechanik 49, 653–663. Kelkar, K.M., Patankar, S.V., 1987. Numerical prediction of flow and heat transfer in a parallel plate channel with staggered fins. Transactions of the ASME—Journal of Heat Transfer 109, 25–30. Kubota, Y., Yamaguchi, T., Tani, K., Takaoka, K., Fujimura, K., Ogura, M., Yamamoto, S., Mizuno, T., Inoue, K., 1993. Anatomical variation of pancreatobiliary ducts in biliary stone diseases. Abdominal Imaging 18, 145–149. Lichtenstein, M.E., Ivy, A.C., 1937. The function of the ‘‘valves’’ of heister. Sugery 1, 38–52. Luo, X.Y., Chin, S.B., Ooi, R.C., Clubb, M., Johnson, A.G., Bird, N., 2004. The rheological properties of human bile, Fourth International Conference on Fluid Mechanics. Dalian, China. Mentzer, S.H., 1926. The valves of heister. Archives of Surgery 13, 511–522. Ooi, R.C., Luo, X.Y., Chin, S.B., Johnson, A.G., Bird, Ni.C., 2003. Fluid–structure interaction simulation of the human cystic duct. ASME, Summer Bioengineering Conference, FL. Otto, W.J., Scott, G.W., Rodkiewicz, C.M., 1979. A comparison of resistances to flow through the cystic duct and the sphincter of oddi. Journal of Surgical Research 27, 68–72. Pitt, H.A., Doty, J.E., DenBesten, L., Kuchenbecker, S.L., 1982. Stasis before gallstone formation: altered gallbladder compliance or cystic resistance? American Journal of Surgery 143, 144–149. Pitt, H.A., Roslyn, J.J., Kuchenbecker, S.L., Doty, J.E., DenBesten, L., 1981. The role of cystic duct resistance in the pathogenesis of cholesterol gallstones. Journal of Surgical Research 30, 508–514. Scott, G.W., Otto, W.J., 1979. Resistance and sphincter-like properties of the cystic duct. Surgery Gynecology and Obstetrics 149, 177–182. Sharp, K.W., Ross, C.B., Tillman, V.N., Williams, L.F., 1990. Changes in gallbladder volume do not affect cystic duct resistance. Archives of Surgery 125, 460–462. Shaw, M.J., Dorsher, P.J., Vennes, J.A., 1993. Cystic duct anatomy: an endoscopic perspective. The American Journal of Gastroenterology 88 (12), 2102–2106. Winkelstein, A., Aschner, P.W., 1926. The mechanism of the flow of bile from the liver into the intestines. American Journal of the Medical Sciences 171, 104–111.