Laminar-to-turbulence and relaminarization zones detection by simulation of low Reynolds number turbulent blood flow in large stenosed arteries

Abstract

Laminar, turbulent, transitional, or combine areas of all three types of viscous flow can occur downstream of a stenosis depending upon the Reynolds number and constriction shape parameter. Neither laminar flow solver nor turbulent models for instance the k–ω (k–omega), k–ε (k–epsilon), RANS or LES are opportune for this type of flow. In the present study attention has been focused vigorously on the effect of the constriction in the flow field with a unique way. It means that the laminar solver was employed from entry up to the beginning of the turbulent shear flow. The turbulent model (k–ω SST Transitional Flows) was utilized from starting of turbulence to relaminarization zone while the laminar model was applied again with onset of the relaminarization district. Stenotic flows, with 50 and 75% cross-sectional area, were simulated at Reynolds numbers range from 500 to 2000 employing FLUENT (v6.3.17). The flow was considered to be steady, axisymmetric, and incompressible. Achieving results were reported as axial velocity, disturbance velocity, wall shear stress and the outcomes were compared with previously experimental and CFD computations. The analogy of axial velocity profiles shows that they are in acceptable compliance with the empirical data. As well as disturbance velocity and wall shear stresses anticipated by this new approach, part by part simulation, are reasonably valid with the acceptable experimental studies.

Keywords

Laminar-turbulent-relaminarization regimes stenosis part by part simulation numerical analysis new approach

1. Introduction

Atherosclerosis is a disease in which deposits of cholesterol and other fatty substances line the walls of arteries. The earliest microscopic change is the accumulation of lipids in the intima [1]. It contains furring of the arteries by plaque accumulation, may bring on heart attack, stroke, or dysfunction of other organs [2]. The development of cholesterol plaques may partially or completely block the blood’s flow through the arteries. It has been accepted that hemodynamics factors are of importance in the initiation and development of atherosclerotic lesions [3]. Since flow disturbances may occur distal to arterial constrictions caused by atherosclerosis, detection and proper interpretation of disturbed flow patterns is of diagnostic value for moderate to advanced disease states; and it might be inferred from several studies that one hope for identifying mild stenosis may be the recognition of early forms of velocity disturbances [4]. Inside the human body, the normal Reynolds number of blood flow differs from 1 to generally 4000 [5]. Artery stenosis is well-known to be approximately related to the existence of a locally disturbed flow as a result of recirculation zones. It can generate pressure loss which is eminent when the internal diameter is decreased beyond around 50% [6]. The post-stenotic flow major characteristic is the presence of perturbation flow or other forms of disturbances. Beside the stenosis separated flow and the mixing layer created between the jet and the surrounding recirculating fluid turns into turbulent flow [7]. The happening of a turbulent flow in the downstream of a constriction has long been identified [8]. Since then, a large number of studies have been done to scientifically recognize the flow disturbances.

Empirical researches containing stenosis have been performed in animal [9,10], in laboratory case [11–14], and in humans [15] that perturbation have been exhibited distal to the stenosis. Steady flow tests have been done to analyze velocity field in 25, 50 and 75% cross-area reductions [13,14]. As a result of the study the onset of turbulence and relaminarization locations have been indicated. Steady and pulsatile flow characteristics have been investigated in models of arterial stenoses [16]. They found, in the case of steady flow, a depiction of the extent of separated flow regions and a range of pressure drops across the stenoses. Kefayati et al. [17] determined the level of turbulence intensity (TI) with respect to definite geometrical characteristics of the plaque namely stenosis severity, eccentricity, and ulceration. A noticeable difference was displayed in turbulence intensity among cases; increasing degree of constriction severity resulted in increased turbulence intensity.

In vitro, evaluation of local velocity profiles and wall shear stress through stenotic vessels have been faced by special problems, so computational fluid dynamics technique can be applied to measure and scrutinize flow field. Various CFD studies have been done in modeling of stenotic flows. The pioneers in the progress and development of numerical simulation were Lee and Fung [18]. Melaaen [19,20], used k–ε model for solving the Reynolds-averaged Navier–Stokes equations. Pinto et al. [21] have reported a numerical study of the blood flow in vessels with a stenosis using OpenFOAM code that characteristics of pulsatile flow and elasticity of a blood flow were regarded.

Generally the blood flow occurs at low Reynolds numbers in the arteries. Under this condition, it is more frequently begotten the laminar, turbulent, and transitional regimes co-occur in the same stenotic flow field with a large region of separated flow [22]. In the present research the k–ω SST Trans Flows solver results are compared with the Ghalichi et al. [23]. They have used the Wilcox low-Re k–ω turbulence model for a steady flow simulation with 50, 75, and 86% stenoses [24]. Analyses presented that the Wilcox’s k–ω turbulence solver have sufficient accuracy for predicting certain areas of the arteries where laminar, transitional, or turbulent flows coincide. Another CFD study was conducted by Lee et al. [25,26], where steady and pulsatile flows were simulated through series stenoses. Effects of the Reynolds number, Womersley number, constriction ratio and spacing ratio of the stenoses on the turbulent flow field were considered. It was shown that the presence more constrictions could trigger earlier occurrence of turbulence in a low Reynolds number. Wilcox’s standard and transitional k–ω model were applied in 75 and 90% area reduction stenosed vessels with particular emphasis on pulsatile flow. Results showed that there is some promise for simulating complex physiological flows using a relatively simple two-equation turbulence model [27]. Keshavarz-Motamed and Kadem have studied a pulsatile blood flow in a simplified case of the aorta (curved pipe) with coexisting coarctation of the aorta and aortic stenosis. FLUENT 6.3 was implemented for flow field simulations [28].

In spite of the studies to date there is an imperfection of CFD works of arterial flow. The discrepancy between previous studies and this new approach, part by part simulation, has been investigated for steady flow through 50 and 75% stenoses at Reynolds numbers from 500 to 2000 in the present paper although preliminary study of the present work has been published by Tabe et al. [22]. For the first time the authors published the results of part by part simulation including vortex length, axial velocity, disturbance velocity and wall shear stresses.

It could be said that laminar area, flow leaves the stenosis in a jet shape which has a comparatively invariable diameter and velocity before growth the separated shear layer to turbulence [29], were evaluated using laminar viscous solver. Then authors concentrated on the turbulent flow field, the jet starts to distribute due to turbulent mixing [29] up to vanishing perturbations, where turbulent (k–ω SST Trans Flows) model [29] was employed. Moving downstream, laminar viscous model was applied to simulate flow field with starting of relaminarization area (Fig. 1). Figure 2 also shows typical streamlines in 50 and 75% stenoses.

Fig. 1.

Schematic of part by part simulation with switching viscous models to each other.

Fig. 2.

Typical streamlines for 50 and 75% stenosis.

Onset of turbulence and relaminarization for 50 and 75% stenoses at $Re = 1000$ and 2000 are listed in Table 1. Some of these positions were mentioned in the previous studies [7,13,14], and the others were specified with the present study by FFT (Fast Fourier Transform) analysis.

Table 1

Onset of turbulence and relaminarization positions ( $\overline{z} = z / D$ where z is measured from the throat of the stenosis)

	50% stenosis		75% stenosis

	$Re = 1000$	$Re = 2000$	$Re = 1000$	$Re = 2000$
Onset of turbulence	$\overline{z} = 7$	$\overline{z} = 3$	$\overline{z} = 4$	$\overline{z} = 1.5$
Relaminarization	$\overline{z} = 9$	$\overline{z} = 8$	$\overline{z} = 7$	$\overline{z} = 8$

2. Methods

The preceding CFD studies have utilized the different methods e.g. k–epsilon (k–ε), k–omega (k–ω), large eddy simulation (LES), and etc. while flow nature varies in the presence of the constriction. For this reason, in this study the part by part simulation was selected to examine and predict flow field parameters.

2.1. Stenosis geometries

Two stenoses configuration were selected for flow simulation. The geometry of the stenoses, identical to the experimental models [13,14], are: $\begin{matrix} (1) & r (z) / D = 0.5 - A [1 + cos π z / D], \end{matrix}$ where $A = 0.073$ and $A = 0.125$ for 50 and 75% constriction, respectively, $- D ⩽ z ⩽ D$ with $Z = 0$ at the throat of the stenosis (Fig. 3), D is the upstream artery diameter (in the current paper, 1 cm) and $r (z)$ is the local radius that alters with axial position z. Stenosis degree is defined as [ $1 - {(r / R)}^{2}$ ] where r and R are the throat and upstream arterial radii.

2.2. Mesh generation and independence test

The geometry and mesh formation were made in Gambit v2.4.6 (ANSYS) which Quad-Map technique (Quadrilateral elements) was applied for producing surface meshing. The details of this step are fully explained in the author’s preliminary article [22]. Finally 105,000 elements in the axisymmetric case were produced to provide the necessary spatial resolution.

2.3. Governing equations and numerical methods

The Navier–Stokes and continuity equations were applied for laminar and relaminarization regions while the Reynolds Averaged Navier–Stokes equations were based for simulation the steady-state turbulent flow.

The shear-stress transport (SST Trans Flows) k–ω model was developed to combine the robust and precise formulation of the k–ω model in the adjacent wall area with the free-stream independence of the k–ε model in the far field effectively. The governing equations of Menter’s k–ω SST model (k and ω are the turbulence kinetic energy and specific dissipation rate proportional to $ε / k$ , respectively) [30] were defined for the steady-state turbulent flow as follows:

Fig. 3.

The geometry of the stenosis.

Turbulent kinetic energy equation: $\begin{matrix} (2) & {(ρ k u_{i})}_{, i} = {((μ + \frac{μ_{t}}{σ_{k}}) k_{i})}_{, i} + G_{k} - Y_{k} + S_{k} . \end{matrix}$

Specific dissipation rate equation: $\begin{matrix} (3) & {(ρ ω u_{i})}_{, i} = {((μ + \frac{μ_{t}}{σ_{ω}}) ω_{i})}_{, i} + G_{ω} - Y_{ω} + D_{ω} + S_{ω} . \end{matrix}$

In these equations, $G_{k}$ represents the generation of turbulence kinetic energy due to mean velocity gradients and $G_{ω}$ characterizes the generation of ω. $Y_{k}$ and $Y_{ω}$ represent the dissipation of k and ω due to turbulence, $D_{ω}$ describes cross diffusion term, $S_{k}$ and $S_{ω}$ are user-defined source terms, $σ_{k}$ and $σ_{ω}$ have been defined in Menter [30] and imply the turbulent Prandtl numbers for k and ω correspondingly.

2.4. Boundary conditions

The boundary conditions for laminar, turbulent and relaminarization regions considered as follows:

For laminar region. $\begin{array}{l} (4) & u_{z} = U [1 - {(r / R)}^{2}], at the inlet of the vessel, \\ (5) & \partial u_{z} / \partial r = u_{r} = 0, at the axis of the vessel, \\ (6) & u_{z} = u_{r} = 0, at the wall of the vessel . \end{array}$

At the outlet has assumed that the flow is fully developed.

For turbulent region. At the inflow boundary was introduced laminar region’s velocity magnitude at the outflow. k and ω were assumed $k = 1.5 I_{T}^{2} V_{in}^{2}$ and 0.0001 respectively [22]. $I_{T}$ is the turbulence intensity and it specifies from experimental results [13,14]. Furthermore: $\begin{array}{l} (7) & \partial u_{z} / \partial r = u_{r} = 0, at the axis of the vessel, \\ (8) & u_{z} = u_{r} = 0, at the wall of the vessel . \end{array}$

For the k–ω SST Trans Flows scheme the non-slip wall boundary conditions have set as $k = 0$ , and $ω = 6 ν / c_{2} y_{1}$ where $y_{1}$ is the height of the first node above the wall. At the outflow boundary were assumed k and ω equal zero, i.e. $\begin{matrix} (9) & \partial k / \partial z = \partial ω / \partial z = 0 \end{matrix}$

For relaminarization region. Turbulent region’s velocity magnitude was applied at velocity inlet boundary. In addition: $\begin{array}{l} (10) & \partial u_{z} / \partial r = u_{r} = 0, at the axis of the vessel, \\ (11) & u_{z} = u_{r} = 0, at the wall of the vessel . \end{array}$

Again at the outflow was supposed a fully developed flow. The fluid was considered homogeneous, incompressible and Newtonian with a constant kinematic viscosity of 3.5 mPa · s and a density of 1060 kgm⁻³. While the shear stress rates in large arteries are adequately large, Morris et al. showed that the supposition of a Newtonian fluid behavior is reasonable for blood flow [31].

3. Results and discussion

Analyses of results include comparison of velocity profiles, disturbance velocity and wall shear stresses (WSS) were done with experimental data [13,14] and previous CFD studies like [23]. The velocity profiles and wall shear stress were normalized by the mean inlet velocity, $V_{in}$ , and dynamic pressure based on mean velocity inlet, $ρ V_{in}^{2}$ , respectively, while the lengths were normalized by the non-stenotic tube diameter, D.

3.1. Velocity profiles

Figure 4(a)–(d) display normalized axial velocity profiles for 50 and 75% stenosis at Reynolds numbers from 500 to 2000. To access the better simulation results, they were then compared to the results of previous studies i.e. Ahmed and Giddens [13,14] and Ghalichi et al. [23]. CFD results by Tan et al. [7] were also used only for location $z = 7$ in Fig. 4(b). The close scrutiny of the pictured velocity profiles (except in Fig. 4(d) at $z = 4$ and $z = 5$ positions) exhibits that they are reasonably close to the experiments [13,14], compared to the former CFD study [23].

Fig. 4.

Comparison of normalized axial velocity profiles at different axial locations by part by part simulation with Ahmed and Giddens [13,14], Ghalichi et al. [23] and Tan et al. [7]. (a) 75% stenosis, $Re = 500$ ; (b) 75% stenosis, $Re = 1000$ ; (c) 75% stenosis, $Re = 2000$ ; (d) 50% stenosis, $Re = 1000$ .

Fig. 4.

(Continued.)

Ahmed and Giddens observed that onset of turbulence and relaminarization occurred at $z = 4$ and $z = 7$ for the 75% stenosis at $Re = 1000$ , respectively [13]. The same conclusion was also demonstrated in a recent CFD study [7]. It can be seen that the velocity profiles captured by the present paper at $z = 4$ and $z = 7$ are satisfactory.

3.2. Disturbance velocity

For 75% stenosis at $Re = 1000$ and 2000 a diagram of non-dimensional axial disturbance velocity component at different axial positions has been reported in the author’s preliminary article [22]. The axial disturbance velocity, $u^{'}$ , has been interpreted as $\sqrt{(2 / 3) k}$ where k is calculated from $0.5 (\overline{u_{i}^{'} u_{i}^{'}})$ , the turbulent kinetic energy [32].

Fig. 5.

Axial disturbance velocity profiles at different positions for 50% stenosis predicted by k–ω SST Trans Flows. (a) $Re = 1000$ ; (b) $Re = 2000$ . Experimental data of Ahmed and Gidddens [13,14] are shown for comparison.

Figure 5 illustrates axial disturbance velocity for 50% stenosis at $Re = 1000$ and 2000. As can be seen, shortcoming of experimental data and CFD studies for analyzing of prediction results by k–ω SST Trans Flows are clear. Figure 5(a) shows axial disturbance velocity distribution at $Re = 1000$ where k–ω SST Trans Flows has predicted the peak value nearly 35% at $z = 8$ . Although there are not experimental or CFD data for validating presented results in Fig. 5, axial disturbance velocity distributions (see Fig. 4 at [22]) are qualitatively close to the experimental measurements.

3.3. Wall shear stress

Axial variation of wall shear stress is pictured in Fig. 6 for 75% constriction at two Reynolds number that separation flow field characterizes by its negative values [22]. It should be noted that a similar distribution of wall shear stress is concluded by new CFD approach (i.e. part by part simulation) but quantitative differences exist for both Reynolds number.

Fig. 6.

Axial variation of wall shear stress for 75% stenosis (a) $Re = 1000$ ; (b) $Re = 2000$ .

Fig. 7.

Wall shear stress at the stenotic throat as a function of Reynolds number.

The variation of wall shear stress with Reynolds number at the stenosis throat is given in Fig. 7. As expected, the highest values of wall stress occur at the constriction throat. For 75% stenosis, the comparison shows that the predicted wall shear stress is approximately 27% higher than the given experimental data which has predicted 40% by Ghalichi et al. [23]. As well as for 50% constriction, although experimental results aren’t available for comparison, it can be seen that part by part simulation describes wall shear stress very close to the experimental values like 75% stenosis.

4. Conclusions

The topic of interest in the present study is the flow field in the neighborhood of a stenosis. The complicated flow features observed in the stenotic flow, involving flow separation, recirculation, reattachment, strong shear layers, and turbulence in both of steady and unsteady flows [33]. The axisymmetric steady flow structure have been simulated in a cylindrical pipe including 50 and 75% stenosis for a Reynolds number ranging from 500–2000. The main emphasis was being on changing viscous solvers to each other. This means that with starting the throat and increasing axial distance, laminar, turbulent, and again laminar viscous models were applied, respectively. The attained results have been validated with previous experiments and CFD studies in terms of velocity profile, disturbance velocity, and wall shear stress.

In general, the analogy between the anticipated outcomes of the part by part simulation and experiments, as criteria, shows not only there are good matches between them but also there are results improvement comparing with the published numerical data [7]. The velocity profiles have been predicted correctly at different axial positions and in terms of disturbance velocity and wall shear stress. This approach presents results identical to the experiments.

The results of this investigation clearly demonstrate that the part by part simulation is a suitable method for stenosed arterial studies and also can be applied to numerical simulation of any constricted tubes in different applications with low Reynolds number that transition from laminar to turbulence and relaminarization is occurred. In such situations, using high Reynolds turbulence modeling is not applicable for whole region of computation.

Footnotes

Acknowledgements

The authors would like to acknowledge Sahand University of Technology for their support.

Conflict of interest

The authors have no conflict of interest to report.

References

Nazemi

et al., Analysis of particle trajectories in aortic artery bifurcations with stenosis, J. Biomech. Eng. 111 (1989), 311–315.

Mitchell

J.R.A.

and Schwartz

C.J.

, Arterial Disease, F.A. Davis Co, Philadelphia, PA, 1965.

Perktold

et al., Wall shear stress distribution in the human carotid siphon during pulsatile flow, J. Biomech. 21 (1988), 663–671. doi:10.1016/0021-9290(88)90203-5.

Ahmed Saad

et al., Pulsatile poststenotic flow studies with laser Doppler anemometry, J. Biomech. 17 (1984), 695–705. doi:10.1016/0021-9290(84)90123-4.

D.N.

, Blood flow in arteries, Annual Review of Fluid Mechanics 29 (1997), 399–434. doi:10.1146/annurev.fluid.29.1.399.

Young

D.F.

, Fluid mechanics of arterial stenoses, J. Biomech. Eng. 101 (1979), 157–175. doi:10.1115/1.3426241.

Tan

F.P.P.

Wood

N.B.

Tabor

and Xu

X.Y.

, Comparison of LES of steady transitional flow in an idealized stenosed axisymmetric artery model with a rans transitional model, J. Biomech. Eng. 133 (2011), 1–12.

Roach

M.R.

, Poststenotic dilatation in arteries, Cardiovascular Fluid Dynamics 2 (1972), 111–139.

Giddens

D.P.

Mabon

R.F.

and Cassanova

R.A.

, Measurements of disordered flows distal to subtotal vastubescular stenoses in the thoracic aortas of dogs, Circulution Rrs. 39 (1976), 112–119. doi:10.1161/01.RES.39.1.112.

10.

Khalifa

A.M.A.

and Giddens

D.P.

, Analysis of disorder in pulsatile flows with application to poststenotic blood velocity measurement in dogs, J. Biomech. 11 (1978), 129–141. doi:10.1016/0021-9290(78)90005-2.

11.

Clark

, Turbulent velocity measurements in a model of aortic stenosis, J. Biomech. 9 (1976), 677–687. doi:10.1016/0021-9290(76)90169-X.

12.

Cassanova

R.A.

and Giddens

D.P.

, Disorder distal to modelled stenoses in steady and pulsatile flow, J. Biomech. 11 (1978), 441–453. doi:10.1016/0021-9290(78)90056-8.

13.

Saad

A.A.

and Giddens

D.P.

, Velocity measurements in steady flow through axisymmetric stenoses at moderate Reynolds numbers, J. Biomech. 16 (1983), 505–516. doi:10.1016/0021-9290(83)90065-9.

14.

Saad

A.A.

and Giddens

D.P.

, Flow disturbance measurements through a constricted tube at moderate Reynolds numbers, J. Biomech. 16 (1983), 955–963. doi:10.1016/0021-9290(83)90096-9.

15.

Strandness

D.E.J.R.

Didisheim

Clowes

A.W.

and Watson

J.T.

(eds), Vascular Disease – Current Research and Clinical Applications, Grune and Stratton, Orlando, USA, 1987.

16.

Young

D.F.

and Tsai

F.Y.

, Flow characteristics in models of arterial stenoses – I. Steady flow, J. Biomech. 6 (1973), 395–410. doi:10.1016/0021-9290(73)90099-7.

17.

Kefayati

Holdsworth

D.W.

and Poepping

, Turbulence intensity measurements using particle image velocimetry in diseased carotid artery models: Effect of stenosis severity, plaque eccentricity, and ulceration, J. Biomech. 47 (2014), 253–263. doi:10.1016/j.jbiomech.2013.09.007.

18.

Lee

J.S.

and Fung

Y.C.

, Flow in locally constricted tubes at low Reynolds numbers, J. App. Mech. 37 (1970), 9–16. doi:10.1115/1.3408496.

19.

Melaaen

M.C.

, Analysis of fluid flow in constricted tubes and ducts using body fitted non-staggered grids, Int. J. Numer. Methods Fluids 15 (1991), 895–923. doi:10.1002/fld.1650150805.

20.

Zijlema

Segal

and Wesselinh

, Finite volume computation of 2D incompressible turbulent flows in general coordinates on staggered grids, Int. J. Numer. Methods in Fluids 20 (1995), 621–640. doi:10.1002/fld.1650200703.

21.

Pinto

S.I.S.

Doutel

Campos

J.B.L.M.

and Miranda

J.M.

, Bload analog fluid flow in vessels with stenosis: Development of an openfoam code to simulate pulsatile flow and elasticity of the fluid, APCBEE Procedia 7 (2013), 73–79. doi:10.1016/j.apcbee.2013.08.015.

22.

Tabe

Ghalichi

Hossainpour

and Ghasemzadeh

, Numerical simulation of transitional blood flow in large arteries, in: 18th Iranian Conference on Biomedical Eng., 2011, pp. 68–71.

23.

Ghalichi

Champlain

A.D.

Douville

King

and Guidoin

, Low Reynolds number turbulence modeling of blood flow in arterial stenoses, Biorheology 35 (1998), 281–294. doi:10.1016/S0006-355X(99)80011-0.

24.

Wilcox

D.C.

, Turbulence Modeling in CFD, DCW Industries, La Canada, California, 1993.

25.

Lee

T.S.

Liao

and Low

H.T.

, Numerical simulation of turbulent flow through series stenoses, Int. J. Numer. Meth. Fluids 42 (2003), 717–740. doi:10.1002/fld.550.

26.

Lee

T.S.

Liao

and Low

H.T.

, Numerical study of physiological turbulent flows through series arterial stenoses, Int. J. Numer. Meth. Fluids 46 (2004), 315–344. doi:10.1002/fld.755.

27.

Ryval

Straatmana

A.G.

and Steinman

D.A.

, Two-equation turbulence modeling of pulsatile flow in a stenosed tube, J. Biomech. Eng. 126 (2004), 625–635. doi:10.1115/1.1798055.

28.

Keshavarz-Motamed

and Kadem

, 3D pulsatile flow in a curved tube with coexisting model of aortic stenosis and coarctation of the aorta, Medical Engineering & Physics 33 (2011), 315–324. doi:10.1016/j.medengphy.2010.10.017.

29.

Jone

S.A.

, A study of turbulent flow downstream of a constriction in a cylindrical tube at low Reynolds numbers with emphasis on frequency correlations, Ph.D. Thesis, Department of Bioengineering, University of California, 1985.

30.

Menter

F.R.

, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA Journal 32 (1994), 1598–1605. doi:10.2514/3.12149.

31.

Morris

Delassus

Callanan

Walsh

Wallis

Grace

et al., 3-D numerical simulation of blood flow through models of the human aorta, J. Biomech. Eng. 127 (2005), 767–775. doi:10.1115/1.1992521.

32.

White

F.M.

, Viscous Fluid Flow, 2nd edn, University of Rhode Island, 1991.

33.

Paul

M.C.

and Molla

M.M.

, Investigation of physiological pulsatile flow in a model arterial stenosis using large-eddy and direct numerical simulations, Applied Mathematical Modelling 36 (2012), 4393–4413. doi:10.1016/j.apm.2011.11.065.