Pore-Scale Modeling of Non-Newtonian Shear-Thinning Fluids in Blood Oxygenator Design May 2016Journal of Biomechanical Engineering 138(5):1-14

Pore-scale modelling of non-Newtonianshear-thinning ﬂuids in blood oxygenator designKenny W. Q. Low ∗Advanced Sustainable Manufacturing Technologies (ASTUTE) Project,College of Engineering, Swansea University, Bay Campus,Fabian Way, Swansea SA1 8EN,United KingdomEmail: k.w.q.low@swansea.ac.ukRaoul van LoonCollege of Engineering, Swansea University, Bay Campus,Fabian Way, Swansea SA1 8EN,United KingdomEmail: r.vanloon@swansea.ac.ukSamuel A. RollandAdvanced Sustainable Manufacturing Technologies (ASTUTE) Project,College of Engineering, Swansea University, Bay Campus,Fabian Way, Swansea SA1 8EN,United KingdomEmail: s.rolland@swansea.ac.ukJohann SienzAdvanced Sustainable Manufacturing Technologies (ASTUTE) Project,College of Engineering, Swansea University, Bay Campus,Fabian Way, Swansea SA1 8EN,United KingdomEmail: j.sienz@swansea.ac.ukABSTRACTThis paper reviews and further develops pore-scale computational ﬂow modelling techniques used for creepingﬂow through orthotropic ﬁbre-bundles used in blood oxygenators. Porous model signiﬁcantly reduce geometricalcomplexity by taking a homogenisation approach to model the ﬁbre bundles. This signiﬁcantly simpliﬁes meshingand can avoid large time-consuming simulations. Analytical relationships between permeability and porosity existfor Newtonian ﬂow through regular arrangements of ﬁbres, and are commonly used in macro-scale porous modelsby introducing a Darcy viscous term in the ﬂow momentum equations. To this extent, veriﬁcation of analyticalNewtonian permeability-porosity relationships have been conducted for parallel and transverse ﬂow through squareand staggered arrangements of ﬁbres. Similar procedures are then used to determine the permeability-porosityrelationship for non-Newtonian blood. Results demonstrate that modelling non-Newtonian shear-thinning ﬂuidsin porous media can be performed via a generalised Darcy equation with a porous medium viscosity decomposedinto a constant term and a directional expression through least squares ﬁtting. This concept is then investigatedfor various non-Newtonian blood viscosity models. The proposed methodology is conducted with two differentporous model approaches; homogeneous and heterogeneous, and validated against a high ﬁdelity model. Results of∗Corresponding author.BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 1 the heterogeneous porous model approach yields improved pressure and velocity distribution which highlights theimportance of wall effects.1 IntroductionPatients’ lives are maintained during open heart surgery and other cardiopulmonary procedures by oxygenating theirblood in cardiopulmonary bypass (CPB) devices. In addition, patients with failing lungs are supported by extracorporealmembrane oxygenation (ECMO) machines [1]. The number of patients is growing rapidly every year in both adults [1, 2]and children [3]. Both CPB and ECMO devices employ mass transfer through gas-permeable hollow ﬁbres to provide bloodwith oxygen and to remove carbon dioxide. Several types of oxygenators have been developed and early records of suchdevices show that a ﬁlm oxygenator was ﬁrst used [4]. An increase in efﬁciency was obtained in the next generation of bloodoxygenators, i.e. the bubble oxygenators, by passing bubbles of oxygen directly through the blood. Later developmentsproduced membrane oxygenators using sheets or ﬁbres and presently hollow ﬁbre membrane oxygenators are the mostcommon type of oxygenators in the market due to their larger exchange surface area for a given volume. However, suchoxygenators can only support low levels of patient activity and are typically bulky pieces of equipment. These limitations ofthe current designs motivate research to design better performing solutions to improve patients’ quality of life.To operate successfully, these devices must be designed to give acceptable pressure drops, stress distribution withinthe blood and residence time distribution. Experimental research is a costly and time consuming design route and alsoaccess to human blood for tests is subject to restrictive regulations. Animal blood is more available and considered safer butdifferences exist in their physical and biological properties compared to human blood. This makes computational modellingan appealing approach due to its ability to rapidly evaluate design variations and reducing cost for the iterative phase ofa product. Computational Fluid Dynamics (CFD) is an effective tool for investigating the blood ﬂow behaviour throughoxygenators. It provides qualitative and quantitative predictions of ﬂuid ﬂows and gives an insight into the ﬂow patterns thatmight be difﬁcult or impossible to study experimentally. Investigation of blood ﬂow involves modelling oxygenators directlywith the presence of the individual physical ﬁbres inside the device. Due to stringent mesh requirements and the presence ofhundreds and thousands of ﬁbres, this would be however, labour intensive in setting up the simulation. To overcome this, asimple approach is to describe the ﬁbrous media inside the oxygenator as a porous medium with the assumption that the ﬁbrebundles are arranged in a structured manner. Such approach is known as the porous model approach. This would eliminatethe time-consuming process of meshing such complex geometry whilst reducing the computational time. In addition tothis, blood ﬂow behaviour at a microscopic level needs to be studied in order to characterise the ﬁbrous medium propertiesand scaled up towards the macroscopic level. The porous model approach here has been applied to model medical devicessuch as blood oxygenator [5–12], dialyser [13] and pump-oxygenator [14]. Blood mainly comprises plasma, red blood cells(RBC), white blood cells (WBC) and platelets. The unique interaction and mechanical behaviour of these cells make bloodbehave as a non-Newtonian ﬂuid. However, under high shear rates, blood tends to behave as a Newtonian ﬂuid [15]. Therehas been active research in modelling the deformation and interactions of various cells in blood in order to enable a betterunderstanding of their respective rheological impact [16–20] and biological processes [21, 22]. Considering the discretecomponents of blood in ﬂow modelling at the scale of a blood oxygenator is currently unfeasibly expensive. Therefore,blood is commonly treated as a single-phase homogeneous continuum model where the rheological properties of the wholeblood are described through non-Newtonian blood viscosity models, developed based on experimental ﬁttings of blood,which depends on the shear rate [21]. Furthermore, there have been works in modelling blood as a two-phase [23–25];plasma and RBC, or three-phase ﬂuids [26]; plasma, RBC and WBC, which involves modelling the plasma (and WBC)as a Newtonian ﬂuid and RBCs as a modiﬁed non-Newtonian ﬂuid based on the haematocrit. This method is used toinvestigate the aggregation of different ﬂuid compositions in blood ﬂow which is not the current scope of this work. In bloodoxygenator applications, blood is conventionally treated in the literature as a single-phase Newtonian ﬂuid for simplicity[5, 7, 8, 11, 12, 14]. Hence, this present work aims to demonstrate how a porous medium approach can be incorporated tomodel non-Newtonian blood ﬂow inside oxygenators.Darcy’s law is widely used for describing Newtonian ﬂuids through porous media at low Reynolds number. Nonetheless,several approaches have been developed for describing non-Newtonian ﬂow through porous media. The common consider-ation is to employ the same Darcy equation with the assumption that the non-Newtonian effects are being lumped togetherinto one coefﬁcient, i.e. the porous medium viscosity [27], which contains the bulk viscosity and the porous medium char-acteristic effects [28–32]. Furthermore, Savins [33] and Sochi [34] focussed on classifying the solution methods for allcharacterised rheology types and highlighted various methods in predicting non-Newtonian ﬂows in porous medium. Thecurrent available literature describes a wide range of complex ﬂuid behaviour in a general porous medium which is beyondthe scope of the present work where it is aimed to investigate the micro-structural effects for blood ﬂow (non-Newtonianshear-thinning ﬂuid) in a uniformly arranged circular ﬁbrous medium. Various studies investigate ﬂuid ﬂow and the correla-tion between porosity and permeability in uniformly arranged circular ﬁbres. These studies have been performed analytically,numerically as well as experimentally throughout the last sixty years. It mostly highlighted the restricted porosity range ofearly models [35–43] but also gathered useful information to evaluate historical agreement between authors. More recently,BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 2 Yazdchi et al. [44] looked at micro-structural effects for ﬂuid ﬂow across ﬁbre beds but no investigations were done forﬂow along ﬁbres. Tamayol et al. [45–49] published a series of works building a description of permeability against porosityfor Newtonian ﬂuid ﬂows along and across ﬁbre beds. These established descriptions in their work [45, 47] for ﬂow alongand across ﬁbres are the focus to validate CFD models for Newtonian ﬂuids. Nevertheless, these work focused solely onNewtonian ﬂuids. Therefore, in this present context, we looked to comprehend the approach from the non-Newtonian ﬂowin general porous media into the ﬂow study of uniformly arranged ﬁbres to provide blood ﬂow macroscopic properties intothe porous model approach.The objective of the current work is to characterise the non-Newtonian blood ﬂow through a network of uniformlyarranged ﬁbres at the microscopic level and use the outcome to build a macro-scale model based on a porous media formu-lation, thus producing an efﬁcient modelling route for ﬁbre bundles. First, an outline of the method is given followed by thegoverning equations of ﬂuid ﬂow in the unit cell structure. Next, the governing equations for Newtonian and non-Newtonianﬂuid ﬂow in orthotropic porous media are described in order to demonstrate the proposed up-scaling methodology from mi-croscopic level to macroscopic level. The numerical setup for the ﬂow simulation of the unit cell structure is then described.The ﬁndings are laid out in Section 2. The ﬁndings are then presented in Section 3, which includes: 1) the validation workfor Newtonian ﬂuids to ensure that the numerical simulation could reproduce the results from the literature, 2) identiﬁca-tion of porous medium viscosities in the porous models for axial/transverse ﬂows through staggered/square arrangementsand 3) validation of the porous methodology by studying the three-dimensional (3D) ﬂow through a small oxygenator andcomparing it against simulations that include individual ﬁbres.2 Materials and MethodsModelling the ﬁbre based blood oxygenator using a porous model approach implies that the ﬁbres are described asa porous medium and assumed to be uniformly arranged. From there, due to periodicity, a representative unit cell can beemployed to investigate the blood ﬂow behaviour. Two different types of ﬁbre packing are presented, i.e a square arrangement(Sq.) and a staggered arrangement (St.). The unit cell deﬁnition for both conﬁgurations are depicted in Fig. 1. The porosity,or void fraction, of the unit cell,εcan be determined fromε=1−πd24S2Sq.1−πd22√3S2St.(1)where the variable term Sis the distance between the centerlines of adjacent ﬁbres (also known as the spacing of the ﬁbres)and dis the ﬁbre diameter. In this work, the porosity considered for square arrangement ranges from 0.22 then 0.3 to 0.9 inincrements of 0.1. For staggered arrangements, the porosity investigated ranges from 0.1 to 0.9 in similar increments.The modelling approach relies on the assumption that non-Newtonian ﬂow in a porous medium can be characterisedbased on the ﬂow rate. In order to achieve this characterisation, ﬂow is modelled at microscopic level to deﬁne the pressuredrop as a function of the ﬂow rate and porosity. Postulating that the permeability is a property of the medium and not theﬂow, any difference in pressure drop across a given domain (unit cell) is due to the ﬂuid property, as it is the only propertyof the ﬂuid involved in governing creeping ﬂows through porous media. This ﬂuid property is known as the equivalentviscosity,µeq, and will be a target to establish a suitable ﬂow dependent shear rate related formulation, called the porousmedium viscosity,µpm, based on different ﬂow rates. The method can be veriﬁed by modelling the ﬂow of non-Newtonianﬂuids through porous media, using the characteristic shear rates to compute the porous medium ﬂuid viscosity matrix to beused in the porous model governing equations.2.1 Governing equationsAccording to the ﬁndings from Zierenberg et al. [50], the Reynolds number for several commercial cross-ﬂow bloodoxygenators in typical operating conditions are found to be ranging from 0.6 to 7. Therefore, ﬂuid ﬂow in the unit cell isconsidered to be incompressible and in low Reynolds number regime (Re ≈0.1−10) [51]. The governing equations of massand momentum can then be expressed in the Cartesian coordinate system as∇∇∇·uuu=0 (2)∇∇∇·(ρuuu⊗uuu) = −∇∇∇P+∇∇∇·τττ(3)BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 3 (a)(b)(c)Fig. 1. (a) A graphical description of the ﬁbre bundles in the oxygenators which can be represented by a unit cell of ordered (b) square and(c) staggered arrays of cylinderswhere ∇∇∇≡ {∂/∂x,∂/∂y,∂/∂z}T,uuu≡ {ux,uy,uz}T,ρ,Pandτττdenote the gradient operator, velocity vector, ﬂuid density,pressure and viscous stress tensor. The termτττcan be expressed asτττ=2µDDD;DDD=12∇∇∇uuu+ (∇∇∇uuu)T(4)whereµis the dynamic viscosity of the ﬂuid and DDDis the rate-of-strain tensor. The superscript ‘T’ denotes the transposeoperator. For non-Newtonian ﬂuid, the viscosity is related to the second invariant of the rate-of-strain tensor ( ˙γ=√2DDD:DDD).Blood is known as a non-Newtonian shear-thinning ﬂuid in nature. There are a number of viscosity models being formulatedand reviewed by Yilmaz and Gundogdu [52]. The Quemada model was chosen in this work as this model has been usedsuccessfully in the literature [25, 53, 54]. It shows sufﬁcient ﬂexibility through its explicit parameters that could accountfor various important elements affecting the blood viscosity [53]. For this reason, it is adopted in the present work and itsformulation is shown below where the shear rate dependent viscosity is given by [55]µ=µp1−12k0+k∞√˙γr1+√˙γrHct−2(5)whereµpis the blood plasma viscosity, k0is the intrinsic viscosity at zero shear, k∞is the intrinsic viscosity at inﬁnite shear,Hct is the haematocrit index and ˙γris expressed as˙γr=˙γ˙γcrit(6)where ˙γcrit denotes the critical shear rate. The parameters in the Quemada model are extracted from the work of Marcinkowska-Gapi´nska et al. [53] whereµp=0.00128 Pa·s, k0=4.01, k∞=1.77, Hct =0.45 and ˙γcrit =4.2 s-1.2.2 Newtonian ﬂow in porous mediumThe volume averaged velocity, known as the superﬁcial velocity, is deﬁned as⟨uuu⟩=1VVfuuudV; (7)where ⟨uuu⟩≡ {⟨ux⟩,uy,⟨uz⟩}T,V,Vfare the velocity vector, total volume and volume of ﬂuid. According to the Darcyequation, the relationship between ⟨uuu⟩through a porous media and ∇∇∇Pis linear in low Re ﬂow regime. The general expres-BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 4 sion in the Cartesian coordinate system is given as⟨uuu⟩=−KKKµ∇∇∇P(8)whereµ,KKKand ∇∇∇Pdenote the dynamic viscosity, permeability tensor and pressure gradient vector respectively. The perme-ability tensor is a symmetric 3 ×3 tensor (Ki j =Kji). There are six unknowns in KKKwhich needs to be determined. Next, thecorresponding components in KKKcan be rearranged in the form of Voigt notation asFFF3×6˜KKK6×1=UUU3×1(9)where FFF=∇xP0 0 ∇yP∇zP00∇yP0∇xP0∇zP0 0 ∇zP0∇xP∇yP;˜KKK=KxxKyyKzzKxyKxzKyzand UUU=⟨ux⟩uy⟨uz⟩The matrix FFFwill contain pressure gradient terms with respect to their Cartesian coordinate directions (x-, y- and z- axes)labelled as a subscript. In order to solve the full ˜KKK, the number of equations needs to be increased which requires constructingFFFand UUUas a family of cases, labelled ‘a’, making Eq. 9 asFFF3a×6˜KKK6×1=UUU3a×1(10)The cases are performed through numerical experiments where mass ﬂow rate, ˙m, is prescribed in one direction (axial ordiagonal). From these experiments, ⟨ui⟩and ∇iPare extracted and a least squares approach is then adopted to obtain sixaccurate permeability values. It should be noted that the accuracy of ˜KKKis dependent of the number of cases constructed,hence, minimising the error in the least square process. Alternatively, KKKcan be diagonalised intoKKKd=K10 00K200 0 K3(11)where the diagonal terms denote the principal permeabilities (1, 2 and 3) accompanied with three respective eigenvectors.2.3 Non-Newtonian ﬂow in porous mediumModelling non-Newtonian ﬂuids in porous media is not straightforward and appropriate modiﬁcations need to be per-formed to the Darcy equation. A common approach is to treat the non-Newtonian ﬂuid as a generalised Newtonian ﬂuidmodel, based on the hypotheses that the ﬂow is steady, and that the same equations used for Newtonian ﬂuids can be appliedto non-Newtonian ﬂuids [27]. Consequently, the constantµterm in Eq. 8 is being replaced by a so-called porous mediumviscosity [28–32]. Furthermore, similarly to non-Newtonian viscosity models, the porous medium viscosity is also relatedto the porous medium shear rate. The general expression of the Darcy equation for non-Newtonian ﬂuids yields the form⟨uuu⟩=−µµµ−1pm ˙γγγpmKKK∇∇∇P(12)whereµµµpm and ˙γγγpm are the porous medium viscosity and porous medium shear rate matrices respectively taking into accountthe ﬁbre direction and arrangement. Both matrices are in diagonal form [28,30–32], which is expressed asµµµpm =diag(µ1(˙γ1),µ2(˙γ2),µ3(˙γ3)). The porous medium shear rate matrix is expressed, in component-wise, as˙γγγpm =1√εαααuuudKKK−1/2(13)BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 5 whereαααis a set of diagonalised geometry calibration constants to allow the up-scaling from micro-scale to macro-scaleand uuudis the ﬂuid velocity matrix arranged in diagonalised form. They contain information on the effect of local geometryinside the porous medium on the local ﬂuid velocity, velocity gradients and the ﬂuid viscosity [32]. The denominator termis known as the characteristic length of the porous medium. The computation of each diagonal component of the geometrycalibration is performed separately. Firstly, ⟨ui⟩and ∇iPare obtained through numerical simulation where ﬂow is introducedin one direction. Next, by utilising the appropriate component of the diagonalised permeability tensor from the Newtoniansimulation or analytical expression, an equivalent viscosity (µeq) can be calculated. This is possible because porosity is aproperty of the media only and not ﬂow dependent. Finally, a suitable geometry calibration constant can be computed toproduce a porous medium viscosity that matches the calculatedµeq. This process here is repeated for the two remainingdirections. Eventually, the requirement to perform any microscopic ﬂuid simulations can be eliminated using this approach.Lopez et al. [30] reported the range ofαto be in the range of 1 to 15 in their experimental work for shear-thinning ﬂuids.However, correlations were difﬁcult to establish and deﬁne especially for low porosity geometries. In this work, an alternativeapproach is presented where the porous medium viscosity matrix (µµµpm) is decomposed into two parameters; one is a constantcharacteristic term and the other contains directional information. This can be expressed asµµµpm =µcˆµµµ(14)whereµcdenotes the scalar characteristic viscosity and similarly is related to a characteristic shear rate, ˙γc, which is writtenbased on dimensional grounds as˙γc=¯u(Sx−d)(Sy−d)ε(15)The terms ¯u,Sxand Sydenote the mean velocity, spacing along the x−and y−axes respectively. The matrix ˆµµµcon-tains the directional effects fromµcwhich depends on the nature of the pore geometry and yields a diagonalised form(ˆµµµ=diag(ˆµ1,ˆµ2,ˆµ3)). The goal here is to obtain a correlation for each diagonal term in ˆµµµfor each ﬁbre conﬁguration withthe aim to subsequently eliminate the requirement of performing micro-scale numerical simulations. Upon obtaining ⟨ui⟩and ∇iPfrom the numerical simulations,µeq can be calculated for eachεby assuming that the permeability is obtained fromexisting analytical (Newtonian ﬂow) expressions. Next, ˙γcis computed from Eq. 15 andµcis calculated subsequently basedon the chosen non-Newtonian viscosity model. Once this is completed, we could potentially obtain a value from ˆµµµwhereµpm,ii matchesµeq,ii. This process is repeated for eachεand correlations could then be determined via least squares ﬁtting.Finally, utilising Eq. 14,µµµpm can be computed which provides the up-scaled property for micro-structural non-Newtonianblood ﬂow effects in the entire porous medium.2.4 Numerical setupThe model geometries are constructed and discretised using DesignModeler and Meshing modules of ANSYS Work-bench v14.5 [56]. A structured mesh is used throughout this work and a mesh arrangement for square and staggered conﬁg-uration for 0.22/0.1, 0.5 and 0.9 porosities are illustrated in Fig. 2. Since the unit cell is repeatable throughout the porousmedia, a periodic boundary condition is seemed appropriate where the ﬂow entering the domain is treated as identical tothe ﬂow exiting the computational domain resulting in a periodic fully-developed ﬂow. A mass ﬂow, ˙m, is prescribed at theinlet and will vary with eachεin order to consistently maintain Re. No-slip boundary conditions are imposed along the ﬁbrewalls. A commercial cell-centered ﬁnite volume solver, ANSYS Fluent®v14.5 [56], is used to solve the ﬂuid conservationequations. The continuity and momentum equations are coupled through an iterative pressure-correction scheme (PISO)where the maximum value of the convergence criterion was set to a value of 1×10−10 on scaled residuals for continuity andmomentum. For blood simulations, the Quemada model is written as a user-deﬁned function (UDF). ANSYS CFD-Post [56]is used for post-processing purposes.A mesh sensitivity study is performed to ensure that the computational grid has minimal impact on the results. Sincedifferent porosities would alter the size of the domain, mesh sensitivity is conducted for eachε. Water is used in this studywhereρandµvalues are taken as 1,000 kg/m3and 0.001 Pa·s respectively. Two types of ﬂow are considered which aretransverse (perpendicular to the ﬁbres) and parallel ﬂows (along the ﬁbres). Transverse ﬂow simulations, in this case alongthe x-axis, are performed in two-dimensional (2D) whereas parallel ﬂow simulations along the z-axis are conducted in 3Dwith a one element thick mesh. For transverse ﬂow, the inlet velocity can be calculated asux=⟨ux⟩S(S−d)Sq.⟨ux⟩√3/2S(√3/2S−d/2)St.(16)BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 6 Fig. 2. Computational grid for ordered (a) square (porosity 0.22, 0.50 and 0.90) and (b) staggered unit cell (porosity 0.10, 0.50 and 0.90).For parallel ﬂow, the inlet velocity in both arrangements are identical and can be computed asuz=⟨uz⟩ε(17)The ¯uare computed from Re with das the characteristic length. The variables ⟨uuu⟩and ∇∇∇Pof both ﬂow directions are thevariables evaluated for mesh convergence. The relative error, efor each variable is deﬁned aseψMi=|ψMi−ψN|ψN;i∈ {1,2,...N−1}(18)whereψis the variable, subscripts iand Ndenote the mesh number and ﬁnest mesh respectively. The convergence criterionis set to be approximately 0.5 % or less to avoid repeating this study for different ﬂuids and ﬂow regimes. Table 1 shows themesh which ensures good precision at a reasonable computational cost for square and staggered arrangements.3 Results and Discussions3.1 Effect of off-diagonal terms in permeability tensorThis section demonstrates the computation of KKK. There are altogether 18 cases (numerical simulations) performedwhere the ﬁrst three ﬂow directions are purely applied in x-, y- and z- axes respectively. The remaining cases are basedon each ﬁxed plane (xy-, xz- and yz- planes) where an additional ﬁve ﬂow directions are applied at 15°, 30°, 45°, 60° and75° respectively. The computed KKKvalues are compared against those considering ﬂow purely in x-, y- and z- directionand they are tabulated in Table 2 alongside their respective principal permeability values. From the results in Table 2, itis demonstrated that the off-diagonal terms in KKK, obtained by the procedure outlined in section 2.2, have very minimalcontributions (≈10−12 to 10−13) leading to minor alterations in the principal permeability values. As expected, we cantherefore treat the principal permeabilities, K1=Kxx,K2=Kyy and K3=Kzz, to be orthogonal and coinciding with thedeﬁned coordinate axes.3.2 Newtonian ﬂuid validationThe validation studies on Newtonian ﬂuids are performed with two objectives. The ﬁrst one is to ensure that the resultsreported in the literature could be reproduced. The second one is to ensure that a change in ﬂuid properties does not affectthe permeability as theoretically predicted from the Darcy equation, which already accounts forµ. To achieve these goals,the simulations are conducted with two ﬂuids used in the development of medical devices. The ﬁrst one is water, used forBIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 7 Table 1. Total number of elements required to achieve mesh independent results for square and staggered arrangements.Square arrangementPorosity 0.22 0.30 0.50 0.70 0.90Div. Narrow 9 12 18 24 30Div. Fibre 33 30 24 18 12Elements 5,400 1,728 2,376 2,880 3,240e∇xP2.83 ×10−33.82 ×10−32.31 ×10−31.72 ×10−31.81 ×10−3e⟨ux⟩1.39 ×10−41.04 ×10−47.08 ×10−59.66 ×10−59.39 ×10−5e∇zP1.54 ×10−31.03 ×10−37.03 ×10−46.07 ×10−46.04 ×10−4e⟨uz⟩1.42 ×10−41.16 ×10−47.79 ×10−55.85 ×10−53.76 ×10−5Staggered arrangementPorosity 0.10 0.30 0.50 0.70 0.90Div. Narrow 12 16 18 24 30Div. Fibre 48 40 24 18 12Elements 15,552 12,000 9,216 9,216 7,680e∇xP2.44 ×10−31.90 ×10−31.73 ×10−31.15 ×10−31.67 ×10−3e⟨ux⟩1.99 ×10−56.05 ×10−67.70 ×10−51.17 ×10−52.15 ×10−4e∇zP2.28 ×10−31.50 ×10−31.73 ×10−31.15 ×10−31.47 ×10−3e⟨uz⟩4.81 ×10−51.60 ×10−52.28 ×10−51.86 ×10−51.03 ×10−5Table 2. Permeability tensor components (left) and principal permeability components (right) for porosity 0.5. Comb. 1 refers to ﬂow purelyin x−,y−and z−axes and Comb. 2 refers to the 18 cases.Comp. Comb. 1 Comb. 2 Comp. Comb. 1 Comb. 2Kxx 4.27 ×10−10 4.26 ×10−10 K14.27 ×10−10 4.24 ×10−10Kyy 4.27 ×10−10 4.27 ×10−10 K24.27 ×10−10 4.28 ×10−10Kzz 1.61 ×10−91.61 ×10−9K31.61 ×10−91.61 ×10−9Kxy – 2.07 ×10−12Kxz – 2.01 ×10−13Kyz – 2.61 ×10−13ﬁrst approach experiments due to its wide availability, and the second ﬂuid is a glycerol and water solution with a viscositymimicking that of blood. The glycerol-mix solution hasρof 1091.6 kg/m3andµof 0.03 Pa·s. The analytical formulationsfrom Tamayol and Bahrami [47] and Gebart [43] are chosen as a baseline to validate the transverse ﬂows through orderlysquare arranged ﬁbres. The non-dimensional transverse ﬂow permeability (K∗1=K1/d2) results of square and staggeredconﬁgurations for both Newtonian ﬂuids are plotted in Fig. 3(a). It can be seen that the predicted K∗1values are in excellentcorrelation with the analytical results of Gebart. The comparison between numerical results and the analytical expressionsof Tamayol and Bahrami, however, are only in good agreement for the square arrangement only. The error is computed tobe approximately 40% for the staggered conﬁguration. As for parallel ﬂows, the analytical expressions of Drummond andTahir [39] and Tamayol and Bahrami [45] are selected. Figure 3(b) illustrates the non-dimensional parallel ﬂow permeability(K∗3=K3/d2) results in square and staggered conﬁgurations for Newtonian ﬂuids. According to both ﬁgures, it can beBIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 8 (a) Transverse Flow (b) Parallel FlowFig. 3. Comparison of non-dimensional permeability for square and staggered conﬁgurations due to (a) transverse and (b) parallel ﬂowsagainst analytical results based on Re =0.1. Numerical results are indicated as (#) and (2) for square and staggered conﬁgurationsrespectively. Analytical expression abbreviations is indicated as (Sq) and (St) for square and staggered conﬁgurations respectively.observed that the numerical results provide exceptional agreement with the analytical models of Tamayol and Bahrami. Theanalytical models of Drummond and Tahir gives excellent agreement as well for high and mid-rangeεbut will start to deviateat lowε.3.3 Newtonian vs. Non-Newtonian velocity ﬁeldsFigures 4(a) to 4(c) and 4(d) to 4(f) illustrate the normalised velocity magnitude contours for transverse ﬂow past squareand staggered arrangements, respectively, at selectedεranges (lowestε, 0.5 and 0.9). For both conﬁgurations, it can beobserved that the results are evident where velocity increases as the ﬂuid approaches the narrowing region and reduces asit draws near the widening region. These incremental and decremental effects of the velocities through different regionsvary signiﬁcantly asεreduces. Furthermore, regions of recirculation are visible in both arrangements for low to moderateporosities, which is illustrated by the streamlines superimposed onto the contours. Corner eddies will occur in the ﬁrstinstance, as seen in Fig. 4(e), and will increase in size as the ﬁbre gap decreases up to a point where they merge to formone central eddy which can be visualised in Fig. 4(b). The size of the central eddy will continue to change its shape (fromﬂat to tall) as the gap decreases. Two or more counter rotating eddies will then form when the gap is reduced until a certainvalue which can be seen in Figs. 4(a) and 4(d) respectively. Flow separation in the staggered arrangement occurs slightlylater compare to the square conﬁguration. This is because the staggered arrangement forces a ﬂow path whereby the ﬂowfollows the surface of the ﬁbres along a greater arc than the purely tangential ﬂow inherent to the square arrangement.Separation is consequently delayed. In commercial blood oxygenators, where ﬂow occurs mainly in transverse direction(across the ﬁbres), the recirculation phenomena (stagnant regions) will promote the formation of clots [57] thus affectingthe overall performance of mass transfer. Furthermore, blood experiences high shear stress when being forced throughsmall gap in between ﬁbres which could potentially cause haemolysis [58]. Figures 4(g) to 4(i) and 4(j) to 4(l) show thenormalised velocity magnitude contours for parallel ﬂow past square and staggered orderly arranged ﬁbres at lowest, 0.5and 0.9 porosities. Again, the velocity distribution for both cases are as expected where maximum velocity will occur in theregion of lowest resistance. If the gap between the ﬁbres is sufﬁciently small, stagnation or extremely low velocity regions(∥uuu∥ ≈ 0) appearing in those gaps could potentially initiate blood clotting [57]. Focussing on transverse ﬂow, a comparisonof velocity proﬁles for Newtonian and non-Newtonian ﬂuids, in the narrowing region, is illustrated in Fig. 5(a) for squareand staggered arrangements respectively. It is shown that the velocity proﬁle in the square arrangement appears more ﬂat(uniform) than that of the staggered arrangement. Their actual differences relative to the non-Newtonian velocity magnitudenormalised by its maximum value are plotted in Fig. 5(b). Based on the results, it is evident that for low porosity medium, thevelocity proﬁle still resembles a parabolic proﬁle due to high shear-thinning effects and the absence of a region unaffectedby the walls, leading to minimal difference between Newtonian and non-Newtonian ﬂuids. The difference, however, startsincreasing where the proﬁle of the non-Newtonian ﬂuid ﬂow begins to ﬂatten, tending towards a plug-ﬂow proﬁle withsteeper gradients near the wall. These steep gradients can be seen at both ends of the plot asεincreases.3.4 Identiﬁcation of the non-Newtonian porous-equivalent viscosity,µµµpmIn this section, the porous medium viscosities (in different directions),µµµpm, are determined as homogenised coefﬁcientsfor the entire porous medium representing the non-Newtonian ﬂow effects found in the micro-structure. To achieve this, theBIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 9 (a)ε=0.22 – Transv. Flow (b)ε=0.50 – Transv. Flow (c)ε=0.90 – Transv. Flow(d)ε=0.10 – Transv. Flow (e)ε=0.50 – Transv. Flow (f)ε=0.90 – Transv. Flow(g)ε=0.22 – Parall. Flow (h)ε=0.50 – Parall. Flow (i)ε=0.90 – Parall. Flow(j)ε=0.10 – Parall. Flow (k)ε=0.50 – Parall. Flow (l)ε=0.90 – Parall. FlowFig. 4. Normalised velocity magnitude contours of variousεfor transverse and parallel ﬂows past square and staggered conﬁgurations.BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 10 0 0.2 0.4 0.6 0.8 100.511.52x 10−3 (a)ε=0.90 – Square & Staggered0 0.2 0.4 0.6 0.8 1−0.08−0.06−0.04−0.0200.020.040.060.08 (b) Normalised actual differenceFig. 5. (a & b) Comparison of normalised velocity proﬁle at the narrowing for Newtonian and non-Newtonian ﬂuids of variousεfor transverseﬂow past square and staggered conﬁgurations. (c) Normalised actual difference relative to non-Newtonian results of velocity proﬁles for thelowest porosity, 0.5 and 0.9 of square and staggered conﬁgurations.procedure laid out in Section 2.3 is followed where ⟨ui⟩and ∇iPare extracted initially from the numerical simulations ofthe micro-structure. In these calculations, viscosity is known through the Quemada model and analytical expressions for thenormalised permeability, KKK, are adopted according to Gebart [43] (transverse ﬂow) and Tamayol and Bahrami [45] (parallelﬂow) as they agreed well with the numerical results for Newtonian ﬂuids. The analytical expressions areK∗1=49π√21−εcrit(1−ε)−12.5;εcrit =1−π/4 Sq.49π√61−εcrit(1−ε)−12.5;εcrit =1−π/2√3St.(19)K∗3=116ϕ−1.479 −lnϕ+2ϕ−ϕ2/2−0.0186ϕ4; Sq.116ϕ−1.498 −lnϕ+2ϕ−ϕ2/2−0.0018ϕ4; St.(20)whereεcrit andϕ(≡1−ε)are known as the critical porosity and solid fraction respectively. The matrixµµµeq can now becalculated through Eq. 12. The resulting components ofµµµeq for each conﬁguration and ﬂow direction are illustrated in Fig.6(a). The range varies approximately from 0.0035 Pa·s up to 0.01 Pa·s in a quadratic manner asεincreases. In transverseﬂow, blood is required to pass through a combination of narrow and wide gaps in between ﬁbres resulting in a higher overallshear than in parallel ﬂow. Hence,µeq,ii is observed to be lower for transverse ﬂow in both ﬁbre conﬁgurations. For lowε, higher velocities are observed at the narrowings in which subsequently leads to high shear exposure causing the bloodto behave like a Newtonian ﬂuid. Lower shear is experienced by the blood in domains of highεwhere the permeability ishigh and the velocity ratio between the narrowing and widening regions is relatively similar, leading to highε. The absolutedifference between both ﬂow directions normalised with their respective parallel ﬂow equivalent viscosity are plotted inFig. 6(b). It can be seen that the staggered arrangement has smaller variation in theirµeq,ii (≈3.6% −9.6%) whereas forthe square arrangement, there are larger differences (≈4.5% −21%). Both normalised differences increase asεdecreases.Next, equatingµµµeq withµµµpm, we can obtain a set of ˆµµµbased on Eq. 14. We then sought to generalise these ﬁndings asmuch as possible so it can be used in a wider set of applications and good correlations were found between ˆµµµandε. Hence,four separate linear expressions are ﬁtted via least squares with a R2correlation of more than 0.95. These correlations areillustrated in Fig. 7 for the Quemada model.3.5 Porous medium viscosities for other non-Newtonian modelsThe previous section (Section 3.4) focused its analysis on the Quemada model. There is a variety of blood viscositymodels being developed and at present, there is no commonly agreed model due to the fact that none of these viscosityBIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 11 0 0.2 0.4 0.6 0.8 100.0020.0040.0060.0080.010.012 (a) Equivalent Viscosity0 0.2 0.4 0.6 0.8 100.0020.0040.0060.0080.010.012 (b) Normalised ErrorFig. 6. (a) Equivalent viscosity computation for square and staggered conﬁgurations with transverse and parallel ﬂows. (b) Normalisedabsolute difference based on their respective parallel ﬂow equivalent viscosities.0 0.2 0.4 0.6 0.8 100.511.52 Fig. 7. The components of ˆµµµaccording to their respective ﬂow direction (transverse and parallel ﬂows) for square, (#) & (2), and staggered,(△)&(3) conﬁgurations. The data are ﬁtted via least squares and expressed in solid lines with the equation and R2correlation.models fully express the behaviour of blood rheology [52]. In this current section, popular blood viscosity models arecompared, namely the Carreau-Yasuda, Cross and Modiﬁed Casson, and are adopted from published work [15, 59]. Theshear rate dependent viscosity expressions for the three non-Newtonian models with their respective coefﬁcients for bloodare tabulated in Table 3 and their shear thinning behaviour in comparison with the Quemada model is illustrated in Fig.8(a). A similar procedure is applied to derive the correlations for these non-Newtonian models, Figs. 8(b), 8(c) and 8(d)are obtained for Carreau-Yasuda, Cross and Modiﬁed Casson models respectively. It can be observed that all correlationsare ﬁtted as linear functions in each ﬂow direction and ﬁbre conﬁguration. When all computedµµµpm are plotted with theequivalent shear rate, computed from the viscosity model based on the equivalent viscosity, the range consistently fallsbetween the shear rate value of 1 s−1to 1000 s−1(see Fig. 9) where there is a slight variation for the lowest porosity of eachconﬁguration. It should be noted that these blood viscosity models are developed based on experimental correlations of bloodand contain designated coefﬁcients. Therefore, the correlations will need to be reassessed if a different set of coefﬁcients areused. Nonetheless, the Quemada model will continue to be used in the subsequent section (Section 4) as this non-Newtonianblood model contains ﬂexible independent coefﬁcients that accounts for important factors affecting the viscosity of blood,such as plasma viscosity, haematocrit, red cell aggregation and deformation [53], compared to the other models.BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 12 Table 3. Non-Newtonian viscosity models for blood (µ), given in Pa·s, as a function of shear rate (˙γ), given in s−1Blood model Dynamic Viscosity,µCarreau-Yasuda model [15, 52]µ=µ∞+ (µ0+µ∞)[1+ (λ˙γ)p](n−1)/pµ∞=0.00345 Pa·s;µ0=0.056 Pa·s;λ=1.902 s; p=0.22Cross model [15, 52]µ=µ∞+(µ0−µ∞)1+(λ˙γ)pµ∞=0.00345 Pa·s;µ0=0.056 Pa·s;λ=1.007 s; p=1.028Modiﬁed Casson model [52, 59]µ=√µp+√τ0√λ+√˙γ2µp=0.002654 Pa·s;τ0=0.0436 Pa·s;λ=0.02181 s4 Methodology ValidationTo demonstrate the effectiveness of this up-scaling, or homogenisation, methodology for non-Newtonian shear-thinningﬂuids, ﬂow inside a small experimental oxygenator without explicitly modelling the ﬁbres, i.e. the porous model, is simulatedand compared against a simulation where every ﬁbre is explicitly modelled, i.e. the high ﬁdelity model. The region of interestis illustrated in Fig. 10(a) with the high ﬁdelity model and porous model shown in Figs. 10(b) and 10(c), respectively. Thecross-sectional dimensions are 17 mm by 8.5 mm. The high ﬁdelity model consists of 380µm diameter ﬁbres arranged in asquare manner, in line with the main ﬂow direction, with a gap (distance between two ﬁbres) of 256µm. There are altogether338 ﬁbres in the device, with 26 ﬁbres in the horizontal direction and 13 layers along the vertical direction. The distancebetween the left and right walls onto the centre of the ﬁrst ﬁbre is approximately 550µm and 434µm for the distance forthe top and bottom walls.The governing equations of ﬂuid ﬂow for the high ﬁdelity model are the incompressible Navier-Stokes equations. Incontrast, the governing equations for the porous model are similar but with an additional Darcy viscous loss term introducedin the momentum equation to describe the ﬁbrous region. The governing equations can be expressed as∇∇∇·uuu=0 (21)∇∇∇·(ρuuu⊗uuu) = −∇∇∇P+∇∇∇·τττ−µµµpmKKK−1εuuu(22)whereµµµpm,KKKare the porous medium viscosity matrix described in Eq. 14 and the analytical Newtonian permeability tensor,respectively. In this section, the porous medium is modelled using two approaches: the fully homogeneous approach, knownas porous model 1, and the heterogeneous approach, denoted as porous model 2 from hereon. The homogeneous porousmodel implies that the porosity mimics the remaining cross-sectional ﬂuid domain deducted from the total number of ﬁbres.On the other hand, a heterogeneous porous model contains 4 different regions of varying porosity, permeabilities and porousmedium viscosity. These are the left and right wall regions (from wall to ﬁrst ﬁbre centre), the bottom and top wall regions,the corner regions and the interior region of the oxygenator. The porous medium viscosity for the interior region is basedon the correlation obtained in Section 3.4. The parameters required for the remaining regions near the wall are based onadditional microscopic numerical ﬂow simulations as well as determination of the porous medium viscosities which areperformed using the existing deﬁned methodology. This would give a more realistic representation of the permeability andporous medium viscosity near the wall.For this validation, the Quemada model is used to describe the blood in the high ﬁdelity model whilst the porous mediumviscosity is used for the porous model. All ﬂuid densities are 1050 kg/m3. A mass ﬂow rate is prescribed at the inlet. Fibreand device walls are treated as no-slip. However, for porous model 2, the walls situated in the ﬁbrous region are treated asfree slip as the porous model used in the adjacent region already accounts for the wall friction losses. The numerical setup isperformed in similar fashion as explained in Section 2.4. The value of convergence criterion is set to a value of 1 ×10−6onscaled residuals for continuity and momentum. The simulations are executed in parallel mode and a total of 8 processors areused on a local 64-bit Windows Intel Xeon machine.The mesh sensitivity study is conducted for the porous model only as computer resources are limited, the full geometrywas therefore modelled using the best computational effort available. Water is used in this study with a ﬂow rate of 150mL/min prescribed at the inlet. Pressure at the inlet (Pin), pressure difference across the inlet and outlet of the porousmedium (∆P) and maximum velocity magnitude (∥uuu∥max) at the outlet are the variables evaluated for mesh convergence.BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 13 10−4 10−2 10010210410−310−210−1100 (a) Non-Newtonian Models (see Table 3)0 0.2 0.4 0.6 0.8 100.511.522.5 (b) Carreau-Yasuda Model0 0.2 0.4 0.6 0.8 100.511.522.5 (c) Cross Model0 0.2 0.4 0.6 0.8 100.511.52 (d) Modiﬁed Casson ModelFig. 8. (a) Comparison of dynamic viscosity (µ) as a function of shear rate ( ˙γ) for various non-Newtonian shear-thinning blood viscositymodels. The components in ˆµµµsubjected to their respective ﬂow direction (transverse and parallel ﬂows) for square, (#)&(2), and staggered,(△)&(3) conﬁgurations of (b) Carreau-Yasuda model, (c) Cross model and (d) Modiﬁed Casson model. The data are ﬁtted via least squaresand expressed in solid lines with the equation and R2correlation.Meshes M1 M2 M3 M4 M5Elements 292,283 1,602,869 2,925,113 4,652,417 10,259,842ePin 0.0045 0.0030 0.0016 0.0010 –e∆P0.0067 0.0032 0.0020 0.0011 –e∥uuu∥max 0.1084 0.0209 0.0116 0.0044 –Table 4. Total number of elements required to achieve mesh independent results for porous model approach.The relative errors based on the ﬁnest mesh for each variable are being evaluated. Table 4 shows the mesh sensitivity studyfor porous model. It can be seen that all three variables decrease when the total number of elements increase. Mesh M3 ischosen (normalised errors are approximately less than or around 1%) to ensure good precision at a reasonable computationalcost.Figures 11(a), 11(b) and 11(c) illustrate the pressure contour normalised by the maximum pressure of the high ﬁdelitymodel. It is observed that pressure decreases across the device and their cross-sectional distribution gradually achieveuniformity in the ﬁbrous region. This illustrates that the ﬂow is well distributed throughout the cross section of the device. Itcan be seen that the pressure contour trend of porous model 1 does not agree well with the high ﬁdelity model at the upstreamregion but gradually achieve proper trends downstream midway through the device. Conversely, the pressure contour trendsin porous model 2 is in excellent agreement with the high ﬁdelity model especially in the porous region. The ∆Pbetween theBIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 14 10010110210310410510−310−210−1 Fig. 9. Porous medium viscosity (µpm) data points for square and staggered conﬁgurations plotted against equivalent shear rate ( ˙γeq) ofvarious non-Newtonian shear-thinning blood viscosity models. Solid line represents non-Newtonian viscosity models from Table 3.(a) Fluid Domain(b) High Fidelity Model(c) Porous ModelFig. 10. (a) Oxygenator domain of interest being modelled in two ways; (b) one with the inclusion of physical ﬁbres (high ﬁdelity model) and(c) without the presence of ﬁbres (porous model).upstream and downstream of the ﬁbrous medium are computed to be 427.13 Pa, 394.53 Pa and 431.70 Pa for high ﬁdelitymodel, porous model 1 and porous model 2 respectively which shows the accuracy of porous model 2, error of 1% relativeto high ﬁdelity model compared to porous model 1 (7.6%).Flow patterns are observed in three positions of the ﬁbrous region namely ‘after inlet’ (Fig. 12(a)), ‘midway’ (Fig.12(b)) and ‘before outlet’ (Fig. 12(c)). The velocity magnitudes are normalised based on their respective maximum velocitymagnitude at various positions and depending on the modelling approach (Fig. 12(d) to Fig. 12(l)). It is shown that bothporous modelling approaches display the anticipated contour trends with the blood entering the ﬁbrous regions at the topleft corner where the ﬂow picks the shortest possible path. Further down the domain the blood is distributed more uniformlybefore it migrates towards the bottom right corner near the exit. After careful inspection of the ﬂow pattern between the twoporous modelling approaches, porous model 2 appears to show better ﬂow trends than porous model 1. This can be visualisedclearly when comparing Figs. 12(e), 12(h) and 12(k), with further help of the overlaying ﬁbres in both porous modellingapproaches. Porous model 2 shows lower ﬂow rate in both top and bottom walls and coincides well with the trends in thehigh ﬁdelity model. This is because the ﬁbre distance is small (244µm) and as a result, the velocity magnitude are not ashigh compared to the interior (256µm), left and right walls (360µm). Nevertheless, the visualisation of the ﬂow distributionfor both porous modelling approaches comparing against the high ﬁdelity model might not be straightforward or intuitive tovisualise and further post-processing of the results are therefore pursued in the following.In order to further verify both porous modelling approaches with the high ﬁdelity model in more detail, ˙mpast theBIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 15 (a) High Fidelity Model (b) Porous Model 1 (c) Porous Model 2Fig. 11. Pressure contours of (a) High ﬁdelity model, (b) Porous model 1 and (c) Porous model 2 normalised by the maximum pressure ofthe high ﬁdelity model.aforementioned three monitoring positons are being evaluated. Each monitoring position consists of 36 grouped cells, dueto computational resource constraints, and they are based on a set of unit cells and wall unit cells. The interior groupedcells consist of 3 ×2 unit cells. The top and bottom grouped cells comprise of the same amount of unit cells as the interiorgrouped cells with 3 additional wall unit cells whereas the left and right grouped cells contain a 2 ×2 unit cells and 2 wallunit cells. Finally, the corner grouped cells contain a 2 ×2 unit cells with 4 wall unit cells and one corner unit cell. Therelative error based on the high ﬁdelity model normalised by the area-weighted ˙mtotal is expressed ase˙m,i=∑Mj=1˙mpm,j−˙mhf,j˙mtotal ∑Mj=1AjAtotal ;i∈ {1,2,...L}(23)where Mand Lrefer to the total amount of unit cells in a grouped cell and wall unit cells and 36 grouped cells, respectively.These errors are plotted in Figs. 13(a) to 13(c) and Figs. 13(d) to 13(f) for porous model 1 and 2 respectively. From theresults, it is shown that the inclusion of the wall effects is important and this can be observed clearly in the results of porousmodel 1. For porous model 1, it is evident that ﬂow, ˙m, at left wall near the inlet is under-predicted and e˙mappears to beextremely high (25%) especially in top left corner (see Fig. 13(a)). This clearly shows that the permeability values prescribedat the left and right wall are too low. The same analysis applies to the area just before outlet as illustrated in Fig. 13(c). Aneven more clear illustration is observed in the centre region of the domain, Figure 13(b), in which the under-prediction of ˙moccurs at the left and right wall (e˙m≈11%). As an overall result for the oxygenator, the mass ﬂows in the top and bottomwalls as well as the interior are over-predicted by approximately 5%. In porous model 2, microscopic numerical simulationsare performed for the wall regions from which more realistic values can be found for the permeability and porous mediumviscosity near the wall. This leads to signiﬁcant improvements in ﬂow distribution. At the region near the inlet, the error inmass ﬂow has reduced substantially, from 25% to 7.6% (see Fig. 13(d)). As a result, the ﬂow is better distributed across thedomain and also reduces the errors in the centre region of the domain. Focussing on the other end of the wall regions, thereappears to be an increase or less signiﬁcant reduction in e˙m. However, it should be noted that their contribution towards thetotal cross-sectional mass ﬂow is insigniﬁcant (approximately one order of magnitude lower). Again, similar interpretationcan be applied to the region of ‘before outlet’ (Fig. 13(d)). Drastic improvement can also be seen in the ‘midway’ region(Fig. 13(e)) and e˙merrors are reduced to approximately 1% whereas the rest of the regions are reduced to 2% and 2.8% forthe interior including the corners and top and bottom wall regions respectively.The porous model approach using the present proposed up-scaling methodology for non-Newtonian blood ﬂow provedto be a good approximation in predicting the overall pressure drop across the ﬁbrous medium and contour trends for pressureand velocity. The porous model approach however suffers from some shortcomings. The devices in real life contain ﬁbresBIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 16 (a) After Inlet (b) Midway (c) Before Outlet(d) High Fidelity Model (e) High Fidelity Model (f) High Fidelity Model(g) Porous Model 1 (h) Porous Model 1 (i) Porous Model 1(j) Porous Model 2 (k) Porous Model 2 (l) Porous Model 2Fig. 12. (a – c) Flow monitoring positions in the device. (d – e) High ﬁdelity model, (g – i) Porous model 1 and (j – l) Porous model 2 velocitymagnitude contours, normalised by their respective maximum magnitude and modelling approaches, at these ﬂow monitoring positions. Notethat the ﬁbres in the porous model are for illustration purposes.arranged in a non-uniform manner that might affect the performance. The porous model, in contrast, assumes that the ﬁbresare arranged in a uniform manner. Furthermore, it has also been demonstrated the importance of wall effect which requirecareful consideration as the wall gaps in actual devices are little-known. Also, local shear stresses for haemolysis evaluationproved to be another limitation of this approach. Nevertheless, this approach are considered to be an efﬁcient method forhigh turnover of design iterations in a short period of time where the porous medium characteristics are required withoutchanging the mesh of the device resulting in the ease of setting up the simulation for evaluating the oxygenator performance.BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 17 (a) Porous Model 1 (b) Porous Model 1 (c) Porous Model 1(d) Porous Model 2 (e) Porous Model 2 (f) Porous Model 2Fig. 13. Relative error normalised by the area weighted mass ﬂow of the high ﬁdelity model at three monitoring positions for (a – c) PorousModel 1 and (d – f) Porous Model 2. Note that the ﬁbres in the porous model are for illustration purposes.5 ConclusionsAn approach is presented for characterising a non-Newtonian shear-thinning ﬂuid in a (ﬁbre bed) porous mediumwhereby a porous medium viscosity replaces the constant viscosity of Newtonian ﬂuids. The porous medium viscosityis decomposed into two distinctive parts: one is a constant characteristic term and the other contains directional effectsdepending on the pore geometry. Two unit cell conﬁgurations are considered describing the transverse and parallel ﬂowthrough a ﬁbre bundle. The full range of porosities is then explored for each conﬁguration. Microscopic scale simulationsare initially performed using Newtonian ﬂuids where the permeability tensor obtained are orthotropic in nature and showexcellent agreement against analytical expressions. In non-Newtonian ﬂuids, a set of directional viscosity values is ﬁttedusing the method validated on the Newtonian ﬂuds.The effectiveness of this proposed methodology is then tested on an experimental oxygenator using Darcy’s law andcompared against the high ﬁdelity model. Two different porous model approaches are adopted with one being a homogeneousporous medium and the other being heterogeneous to account for larger gaps between the device’s wall and the ﬁbre bundle.The pressure and velocity data throughout the device show a good agreement between high-ﬁdelity and porous models but asigniﬁcant improvement is made when accounting for the wall-ﬁbre gaps. Further mass ﬂow rate validations on the inclusionof near-wall effects show that the homogeneous porous model approach tends to under-predict the ﬂow in near-wall regionscontaining larger gaps compared to the interior regions and vice versa for near-wall regions with a smaller gaps. In contrast,the heterogeneous porous model approach yields a much improved ﬂow distribution due to the more detailed description ofthe volume’s porosity variations. This highlights the importance of including the wall effects to capture the ﬂow features andrequires further investigation.Nevertheless, the proposed up-scaling or homogenisation methodology for non-Newtonian blood ﬂow proves to be aneffective and accurate way to describe the ﬂow through the ﬁbre bundles of a blood oxygenator and to predict the pressuredrop across the device. It should be noted that careful consideration of the near-wall regions greatly improves the predictions.Therefore, future works include a ﬁner investigation of the wall effects and expanding to higher Reynolds number.AcknowledgementsThe authors would like to acknowledge the ﬁnancial support of the Welsh Government, the European Regional Devel-opment Fund (ERDF) and the funding provided via the ASTUTE 2020 operation. They would also like to thank HaemairLtd. for their cooperation and assistance throughout the project.BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 18 NomenclatureA=Area.Atotal =Total area.d=Fibre diameter.DDD=Rate-of-strain tensor.e=Relative error.Hct =Haematocrit.KKK=Permeability tensor.KKKd=Diagonalised permeability tensor.k0=Rate of rouleaux zero shear.k∞=Rate of rouleaux at inﬁnite shear.˙m=Mass ﬂow rate.P=Pressure.Re =Reynolds number.S=Fibre spacing.Sx=Fibre spacing in x-axis.Sy=Fibre spacing in y-axis.S=Fibre spacing.uuud=Diagonalised velocity matrix.ux=Velocity in x-axis.uz=Velocity in z-axis.¯u=Mean velocity.uuu=Velocity vector.⟨uuu⟩=Volume averaged velocity vector.V=Total volume.Vf=Volume of ﬂuid.ααα=Shift factor matrix.ε=Porosity.εcrit =Critical porosity.˙γ=Shear rate.˙γc=Characteristic shear rate.˙γr=Shear rate ratio.˙γcrit =Critical shear rate.˙γpm =Porous medium shear rate matrix.∇∇∇=Gradient operator.µ=Dynamic viscosity.µc=Characteristic viscosity.µµµeq =Equivalent viscosity matrix.µp=Blood plasma viscosity.µµµpm =Porous medium viscosity matrix.ˆµµµ=Directional viscosity matrix.ϕ=Solid fraction.ψ=Variable.τττ=Viscous stress tensor.BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 19 References[1] Turner, D. A., and Cheifetz, I. M., 2013. “Extracorporeal membrane oxygenation for adult respiratory failure”. Respi-ratory Care, 58(6), pp. 1038–1052.[2] Schmidt, M., Hodgson, C., and Combes, A., 2015. “Extracorporeal gas exchange for acute respiratory failure in adultpatients: a systematic review”. Critical Care, 19(99), pp. 1–14.[3] Brogan, T. V., Zabrocki, L., Thiagarajan, R. R., Rycus, R. H., and Bratton, S. L., 2012. “Prolonged extracorporealmembrane oxygenation for children with respiratory failure”. Pediatric Critical Care Medicine, 13(4), pp. e249–e254.[4] Gibbon Jr., J. H., 1937. “Artiﬁcial maintenance of circulation during experimental occlusion of pulmonary artery”.Archives of Surgery, 34(6), pp. 1105–1131.[5] Gartner, M. J., Wilhelm, C. R., Gage, K. L., Fabrizio, M. C., and Wagner, W. R., 2000. “Modeling ﬂow effects onthrombotic deposition in a membraneoxygenator”. Artiﬁcial Organs, 24(1), pp. 29–36.[6] Gage, K. L., Gartner, M. J., Burgreen, G. W., and Wagner, W. R., 2002. “Predicting membrane oxygenator pressuredrop using Computational Fluid Dynamics”. Artiﬁcial Organs, 26(7), pp. 600–607.[7] Mazaheri, A. R., and Ahmadi, G., 2006. “Uniformity of the ﬂuid ﬂow velocities within hollow ﬁber membranes ofblood oxygenation devices”. Artiﬁcial Organs, 30(1), pp. 10–15.[8] Zhang, J., Nolan, T., Zhang, T., Grifﬁth, B. P., and Wu, Z. J. J., 2007. “Characterization of membrane blood oxygenationdevices using computational ﬂuid dynamics”. Journal of Membrane Science, 288(12), pp. 268 – 279.[9] Bhavsar, S. S., Schmitz-Rode, T., and Steinseifer, U., 2011. “Numerical modeling of anisotropic ﬁber bundle behaviorin oxygenators”. Artiﬁcial Organs, 35(11), pp. 1095–1102.[10] Sano, Y., Adachi, J., and Nakayama, A., 2013. “A porous media theory for characterization of membrane bloodoxygenation devices”. Heat and Mass Transfer, 49(7), pp. 973–984.[11] Pelosi, A., Sheriff, J., Stevanella, M., Fiore, G. B., Bluestein, D., and Redaelli, A., 2014. “Computational evalua-tion of the thrombogenic potential of a hollow-ﬁber oxygenator with integrated heat exchanger during extracorporealcirculation”. Biomechanics and Modeling in Mechanobiology, 13(2), pp. 349–361.[12] Consolo, F., Fiore, G. B., Pelosi, A., Reggiani, S., and Redaelli, A., 2015. “A numerical performance assessment of acommercial cardiopulmonary by-pass blood heat exchanger”. Medical Engineering & Physics, 37(6), pp. 584 – 592.[13] Sano, Y., and Nakayama, A., 2012. “A porous media approach for analyzing a countercurrent dialyzer system”. Journalof Heat Transfer, 134(7), pp. 1–11.[14] Zhang, J., Taskin, M. E., Koert, A., Zhang, T., Gellman, B., Dasse, K. A., Gilbert, R. J., Grifﬁth, B. P., and Wu,Z. J., 2009. “Computational design and in vitro characterization of an integrated maglev pump-oxygenator”. ArtiﬁcialOrgans, 33(10), pp. 805–817.[15] Cho, Y. I., and Kensey, K. R., 1991. “Effects of the non-newtonian viscosity of blood on ﬂows in a diseased arterialvessel. part 1: Steady ﬂows”. Biorheology, 28(3-4), pp. 241 – 262.[16] Boryczko, K., Dzwinel, W., and Yuen, D. A., 2003. “Dynamical clustering of red blood cells in capillary vessels”.Journal of Molecular Modeling, 9(1), pp. 16–33.[17] Fedosov, D. A., Caswell, B., and Karniadakis, G. E., 2010. “A multiscale red blood cell model with accurate mechanics,rheology, and dynamics”. Biophysical Journal, 98(10), pp. 2215 – 2225.[18] Bessonov, N., Babushkina, E., Golovashchenko, S. F., Tosenberger, A., Ataullakhanov, F., Panteleev, M., Tokarev, A.,and Volpert, V., 2014. “Numerical modelling of cell distribution in blood ﬂow”. Mathematical Modelling of NaturalPhenomena, 9, pp. 69–84.[19] Fedosov, D. A., Noguchi, H., and Gompper, G., 2013. “Multiscale modeling of blood ﬂow: from single cells to bloodrheology”. Biomechanics and Modeling in Mechanobiology, 13(2), pp. 239–258.[20] Gambaruto, A. M., 2015. “Computational haemodynamics of small vessels using the Moving Particle Semi-implicit(MPS) method”. Journal of Computational Physics, 302, pp. 68 – 96.[21] Popel, A. S., and Johnson, P. C., 2005. “Microcirculation and hemorheology”. Annual Review of Fluid Mechanics,37(1), pp. 43–69.[22] Lipowsky, H. H., 2005. “Microvascular rheology and hemodynamics”. Microcirculation, 12, pp. 5 – 15.[23] Das, B., Enden, G., and Popel, A. S., 1997. “Stratiﬁed multiphase model for blood ﬂow in a venular bifurcation”.Annals of Biomedical Engineering, 25(1), pp. 135–153.[24] Jung, J., Lyczkowski, R. W., Panchal, C. B., and Hassanein, A., 2006. “Multiphase hemodynamic simulation ofpulsatile ﬂow in a coronary artery”. Journal of Biomechanics, 39(11), pp. 2064 – 2073.[25] Kim, Y. H., VandeVord, P. J., and Lee, J. S., 2008. “Multiphase non-Newtonian effects on pulsatile hemodynamics ina coronary artery”. International Journal for Numerical Methods in Fluids, 58(7), pp. 803–825.[26] Jung, J., and Hassanein, A., 2006. “Three-phase CFD analytical modeling of blood ﬂow”. Medical Engineering &Physics, 30(1), pp. 91–103.[27] Bird, R. B., Stewart, W. E., and Lightfoot, E. N., 2006. Transport Phenomena, 2nd ed. John Wiley & Sons, Inc., NewYork.[28] Sorbie, K. S., Clifford, P. J., and Jones, E. R. W., 1989. “The rheology of pseudoplastic ﬂuids in porous media usingBIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 20 network modeling”. Journal of Colloid and Interface Science, 130(2), pp. 508 – 534.[29] Pearson, J. R. A., and Tardy, P. M. J., 2002. “Models for ﬂow of non-Newtonian and complex ﬂuids through porousmedia”. Journal of Non-Newtonian Fluid Mechanics, 102(2), pp. 447 – 473.[30] Lopez, X., Valvatne, P. H., and Blunt, M. J., 2003. “Predictive network modeling of single-phase non-Newtonian ﬂowin porous media”. Journal of Colloid and Interface Science, 264(1), pp. 256 – 265.[31] Perrin, C. L., Tardy, P. M. J., Sorbie, K. S., and Crawshaw, J. C., 2006. “Experimental and modeling study of Newtonianand non-Newtonian ﬂuid ﬂow in pore network micromodels”. Journal of Colloid and Interface Science, 295(2), pp. 542– 550.[32] Tosco, T., Marchisio, D. L., Lince, F., and Sethi, R., 2013. “Extension of the darcy-forchheimer law for shear–thinningﬂuids and validation via pore-scale ﬂow simulations”. Transport in Porous Media, 96(1), pp. 1–20.[33] Savins, J. G., 1969. “Non-Newtonian ﬂow through porous media”. Industrial & Engineering Chemistry, 61(10),pp. 18–47.[34] Sochi, T., 2010. “Non-Newtonian ﬂow in porous media”. Polymer, 51(22), pp. 5007–5023.[35] Happel, J., 1959. “Viscous ﬂow relative to arrays of cylinders”. AIChE Journal, 5(2), pp. 174–177.[36] Hasimoto, H., 1959. “On the periodic fundamental solutions of the stokes equations and their application to viscousﬂow past a cubic array of spheres”. Journal of Fluid Mechanics, 5, 2, pp. 317–328.[37] Sangani, A., and Acrivos, A., 1982. “Slow ﬂow through a periodic array of spheres”. International Journal of Multi-phase Flow, 8(4), pp. 343–360.[38] Sangani, A., and Acrivos, A., 1982. “Slow ﬂow past periodic arrays of cylinders with application to heat transfer”.International Journal of Multiphase Flow, 8(3), pp. 193–206.[39] Drummond, J. E., and Tahir, M. I., 1984. “Laminar viscous ﬂow through regular arrays of parallel solid cylinders”.International Journal of Multiphase Flow, 10(5), pp. 515 – 540.[40] Sparrow, E. M., and Loefﬂer, A. L., 1959. “Longitudinal laminar ﬂow between cylinders arranged in regular array”.AIChE Journal, 5(3), pp. 325–330.[41] Bruschke, M. V., and Advani, S. G., 1993. “Flow of generalized Newtonian ﬂuids across a periodic array of cylinders”.Journal of Rheology (1978-present), 37(3), pp. 479–498.[42] der Westhuizen, J. V., and Plessis, J. P. D., 1996. “An attempt to quantify ﬁbre bed permeability utilizing the phaseaverage Navier-Stokes equation”. Composites Part A: Applied Science and Manufacturing, 27(4), pp. 263 – 269.[43] Gebart, B. R., 1992. “Permeability of unidirectional reinforcements for RTM”. Journal of Composite Materials, 26(8),pp. 1100–1133.[44] Yazdchi, K., Srivastava, S., and Luding, S., 2011. “Microstructural effects on the permeability of periodic ﬁbrousporous media”. International Journal of Multiphase Flow, 37(8), pp. 956–966.[45] Tamayol, A., and Bahrami, M., 2010. “Parallel ﬂow through ordered ﬁbers: an analytical approach”. Journal of FluidsEngineering, 132(11), pp. 1–7.[46] Tamayol, A., and Bahrami, M., 2008. “Numerical investigation of ﬂow in ﬁbrous porous media”. In ECI InternationalConference on Heat Transfer and Fluid Flow in Microscale, Whistler, Canada, Sep, pp. 21–26.[47] Tamayol, A., and Bahrami, M., 2011. “Transverse permeability of ﬁbrous porous media”. Phys. Rev. E, 83, Apr,p. 046314.[48] Tamayol, A., and Bahrami, M., 2009. “Analytical determination of viscous permeability of ﬁbrous porous media”.International Journal of Heat and Mass Transfer, 52(9-10), pp. 2407 – 2414.[49] Tamayol, A., and Bahrami, M., 2011. “In-plane gas permeability of proton exchange membrane fuel cell gas diffusionlayers”. Journal of Power Sources, 196, pp. 3559 – 3564.[50] Zierenberg, J. R., Fujioka, H., Hirschl, R. B., Bartlett, R. H., and Grotberg, J. B., 2007. “Pulsatile blood ﬂow andoxygen transport past a circular cylinder”. Journal of Biomechanical Engineering, 129(2), pp. 202–215.[51] Edwards, D. A., Shapiro, M., BarYoseph, P., and Shapira, M., 1990. “The inﬂuence of Reynolds number upon theapparent permeability of spatially periodic arrays of cylinders”. Physics of Fluids A, 2(1), pp. 45–55.[52] Yilmaz, F., and Gundogdu, M. Y., 2008. “A critical review on blood ﬂow in large arteries; relevance to blood rheology,viscosity models, and physiologic conditions”. Korea-Australia Rheology Journal, 20(4), pp. 197–211.[53] Marcinkowska-Gapi´nska, A., Gapinski, J., Elikowski, W., Jaroszyk, F., and Kubisz, L., 2007. “Comparison of threerheological models of shear ﬂow behavior studied on blood samples from post-infarction patients”. Medical & Biolog-ical Engineering & Computing, 45(9), pp. 837–844.[54] Johnston, B. M., Johnston, P. R., Corney, S., and Kilpatrick, D., 2004. “Non-Newtonian blood ﬂow in human rightcoronary arteries: steady state simulations”. Journal of Biomechanics, 37(5), pp. 709 – 720.[55] Quemada, D., 1978. “Rheology of concentrated disperse systems II. a model for non-Newtonian shear viscosity insteady ﬂows”. Rheologica Acta, 17(6), pp. 632–642.[56] ANSYS IN C., 2011. ANSYS Documentations., release 14.5 ed. USA.[57] Fasano, A., Santos, R., and Sequeira, A., 2012. “Blood coagulation: A puzzle for biologists, a maze for mathemati-cians”. In Modeling of Physiological Flows, D. Ambrosi, A. Quarteroni, and G. Rozza, eds., Vol. 5 of MS&A Modeling,BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 21 Simulation and Applications. Springer Milan, Milano, pp. 41–75.[58] Leverett, L., Hellums, J., Alfrey, C., and Lynch, E., 1972. “Red blood cell damage by shear stress”. BiophysicalJournal, 12(3), pp. 257 – 273.[59] Buchanan Jr., J. R., Kleinstreuer, C., and Comer, J. K., 2000. “Rheological effects on pulsatile hemodynamics in astenosed tube”. Computers & Fluids, 29(6), pp. 695 – 724.BIO-15-1501 Low et al. [DOI: 10.1115/1.4032801] 22

Public Last updated: 2024-09-03 07:40:27 AM

Comments