Citation
Numerical modeling of non-Newtonian flows using central moment lattice Boltzmann method

Material Information

Title:
Numerical modeling of non-Newtonian flows using central moment lattice Boltzmann method
Creator:
Adam, Saad
Place of Publication:
Denver, CO
Publisher:
University of Colorado Denver
Publication Date:
Language:
English

Thesis/Dissertation Information

Degree:
Doctorate ( Doctor of philosophy)
Degree Grantor:
University of Colorado Denver
Degree Divisions:
College of Engineering and Applied Sciences, CU Denver
Degree Disciplines:
Engineering and applied science
Committee Chair:
Welch, Samuel
Committee Members:
Premnath, Kannan
Calvisi, Michael
Jenkins, Peter
Carey, Varis

Notes

Abstract:
Fluid mechanics of non-Newtonian fluids, which arise in numerous settings, are characterized by non-linear constitutive models that pose challenges for computational methods. Lattice Boltzmann methods (LBMs) o↵ers some computational advantages due to its kinetic basis and its simpler and local stream-and-collide procedure enabling efficient simulations. In this thesis, we will develop and apply advanced formulations of the LBM with improved accuracy and numerical stability for computations of non-Newtonian fluids involving broader parameter ranges. In this regard, we present cascaded LBMs based on central moments and using multiple relaxation times to simulate non-Newtonian power-law fluids in two- (2D) and three- (3D) dimensions and including forcing terms. The relaxation times of the second order moments are varied locally based on the local shear rate that is based on the second invariant of the strain rate tensor and parameterized by the consistency coefficient and the power law index of the constitutive relation of the power law fluid. The 2D and 3D cascaded LBM formulations are based on two-dimensional, nine velocity (D2Q9) lattice, and three-dimensional, nineteen velocity (D3Q19) lattice, respectively, where the e↵ects of collisions as well as the impressed forces are represented in terms of appropriate changes in central moments of di↵erent orders. Numerical validation study of these new LBMs for various benchmark rheological flows that encompass both shear thinning and shear thickening behavior and including comparisons against recent numerical solutions for non-Newtonian cubic cavity flow, are presented. Furthermore, numerical stability comparisons of the proposed advanced LBM against the LBM based on a single relaxation time model are made. Numerical results demonstrate the accuracy, second order grid convergence and improvements in stability for simulations of non-Newtonian fluid flows with the use of the cascaded LBM. Finally, as an illustrative application, the e↵ects of various characteristic parameters on the wake interactions and drag forces for non-Newtonian flow of power law fluids over a pair of cylinders in tandem arrangement are investigated.

Record Information

Source Institution:
University of Colorado Denver
Holding Location:
Auraria Library
Rights Management:
Copyright Adam Saad. Permission granted to University of Colorado Denver to digitize and display this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.

Downloads

This item is only available as the following downloads:


Full Text

PAGE 1

NUMERICALMODELINGOFNON-NEWTONIANFLOWS USINGCENTRALMOMENTLATTICEBOLTZMANN METHOD by SAADADAM BS,UniversityofTobruk,2001 MS,UniversityofColoradoatBoulder,2010 Athesissubmittedtothe FacultyoftheGraduateSchoolofthe UniversityofColoradoinpartialfulÞllment oftherequirementsforthedegreeof DoctorofPhilosophy EngineeringandAppliedScienceProgram 2018

PAGE 2

ii ThisthesisfortheDoctorofPhilosophydegreeby SaadAdam hasbeenapprovedforthe EngineeringandAppliedScienceProgram SamuelWelch,Chair KannanPremnath,Advisor MichaelCalvisi PeterJenkins VarisCarey Date:December15,2018

PAGE 3

iii Adam,Saad(PhD,EngineeringandAppliedSciences) NUMERICALMODELINGOFNON-NEWTONIANFLOWSUSINGCENTRALMOMENT LATTICEBOLTZMANNMETHOD ThesisdirectedbyKannanPremnath ABSTRACT Fluidmechanicsofnon-Newtonianßuids,whichariseinnumeroussettings,arecharacterized bynon-linearconstitutivemodelsthatposechallengesforcomputationalmethods.Lattice Boltzmannmethods(LBMs)o!erssomecomputationaladvantagesduetoitskineticbasisand itssimplerandlocalstream-and-collideprocedureenablinge"cientsimulations.Inthisthesis, wewilldevelopandapplyadvancedformulationsoftheLBMwithimprovedaccuracyand numericalstabilityforcomputationsofnon-Newtonianßuidsinvolvingbroaderparameterranges. Inthisregard,wepresentcascadedLBMsbasedoncentralmomentsandusingmultiple relaxationtimestosimulatenon-Newtonianpower-lawßuidsintwo-(2D)andthree-(3D) dimensionsandincludingforcingterms.Therelaxationtimesofthesecondordermomentsare variedlocallybasedonthelocalshearratethatisbasedonthesecondinvariantofthestrain ratetensorandparameterizedbytheconsistencycoe"cientandthepowerlawindexofthe constitutiverelationofthepowerlawßuid.The2Dand3DcascadedLBMformulationsare basedontwo-dimensional,ninevelocity(D2Q9)lattice,andthree-dimensional,nineteenvelocity (D3Q19)lattice,respectively,wherethee!ectsofcollisionsaswellastheimpressedforcesare representedintermsofappropriatechangesincentralmomentsofdi!erentorders.Numerical validationstudyofthesenewLBMsforvariousbenchmarkrheologicalßowsthatencompassboth shearthinningandshearthickeningbehaviorandincludingcomparisonsagainstrecentnumerical solutionsfornon-Newtoniancubiccavityßow,arepresented.Furthermore,numericalstability comparisonsoftheproposedadvancedLBMagainsttheLBMbasedonasinglerelaxationtime modelaremade.Numericalresultsdemonstratetheaccuracy,secondordergridconvergenceand

PAGE 4

iv improvementsinstabilityforsimulationsofnon-Newtonianßuidßowswiththeuseofthe cascadedLBM.Finally,asanillustrativeapplication,thee!ectsofvariouscharacteristic parametersonthewakeinteractionsanddragforcesfornon-Newtonianßowofpowerlawßuids overapairofcylindersintandemarrangementareinvestigated. Theformandcontentofthisabstractareapproved.Irecommenditspublication. Approved:KannanPremnath

PAGE 5

v DEDICATION Thisthesisisdedicatedtomyfamilyfortheirunderstanding,motivationandpatienceduringthe challengesofgraduateschoolandlife.Iamtrulythankfulforhavingallofyouinmylife.This workisalsodedicatedtomysupervisor,Dr.KannanPremnath,forthetimethathehasspent withmeasIcarriedoutthisproject.

PAGE 6

vi ACKNOWLEDGEMENTS Thisworkcouldnothavebeencompletedbutforthepatronageandassistanceofmanypeople. Firstandtheforemost,IwouldliketoexpressmyprofoundgratitudetowardsmyadvisorKannanPremnathforhisablestguidanceandvaluablesuggestionsthroughoutthedurationofthis project.Hisadviceandencouragementinspiredmetoimproveandimplementmyworksuccessfully.Hiskindnessandhelpareappreciated.Iamverythankfultothemembersofmythesisdefensecommittee,Dr.SamuelWelch,Dr.PeterJenkins,Dr.VarisCareyandDr.MichaelCalvisi, forspendingmuchtimereadingthisthesisandprovidingmanyconstructivecomments.Special thanksarealsoduetoallthefriendsfortheirincessantencouragement.Lastbutnottheleast,I appreciatetheundyingsupportofmyfamilytowhomthisworkissincerelydedicated.

PAGE 7

vii TABLEOFCONTENTS CHAPTER I INTRODUCTION1 1.1Motivation......................................... 1 1.2ParallelComputing.................................... 4 1.2.1ClassiÞcationofParallelPlatforms........................ 5 1.2.2ParallelizationofLBMalgorithms........................ 6 1.3Focusofthiswork..................................... 7 1.4OutlineofThesis...................................... 8 II TWO-DIMENSIONALCASCADEDLBMFORNON-NEWTONUANFLUIDFLOWS10 2.1Introduction......................................... 10 2.2Macroscopicgoverningequationsforthedynamicsofnon-Newtonianßuidmotion.. 14 2.3CascadedLBMbasedonCentralmomentsfornon-Newtonianßuidßow....... 16 2.4ResultsandDiscussion................................... 25 2.4.1Channelßowofpower-lawandBinghamßuids................. 25 2.4.22Dlid-drivencavityßowofnon-Newtonianßuids................ 29 2.4.3Non-Newtonianßowofpowerlawßuidsoverasinglecircularcylinder.... 39 2.4.4Non-Newtonianßowofpowerlawßuidsoverapairofcylindersintandem arrangement.................................... 46 2.5Comparisonofnumericalstabilityofdi!erentcollisionmodels............. 55 2.6SummaryandConclusions................................. 55 III THREE-DIMENSIONALCASCADEDLBMFORNON-NEWTONUANFLUIDFLOWS58 3.1Introduction......................................... 58 3.2Governingequationsforthree-dimensionalgeneralizednon-Newtonianßuidßows.. 61 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

PAGE 8

viii 3.3Three-dimensionalcascadedLBMfornon-Newtonianßuidßows........... 63 3.3.1SelectionofBasisVectorsforMoments..................... 63 3.3.2ContinuousCentralMoments:DistributionFunction,itsLocalAttractor andForcing..................................... 66 3.3.3CascadedLattice-BoltzmannEquationwithForcingTerms.......... 68 3.3.4VariousDiscreteCentralMomentsandGalileanInvarianceMatchingPrinciple......................................... 69 3.3.5VariousDiscreteRawMomentsandSourceTermsinParticleVelocitySpace 71 3.3.6StructureoftheCascadedCollisionOperator.................. 73 3.3.7OperationalStepsoftheCentralMomentLBM................. 78 3.4Resultsanddiscussion................................... 79 3.4.1Validationofthe3DcascadedLBM....................... 79 3.5Comparisonofnumericalstabilityofdi!erentcollisionmodels............. 83 3.6SummaryandConclusions................................. 84 IV SUMMARY,CONCLUSIONSANDOUTLOOK95 BIBLIOGRAPHY 97 APPENDIX A 103 1.1Chapman-EnskoganalysisofthecentralmomentLBM:Cubicvelocityerrors, equilibrium-momentscorrectionsandstrainratetensorformulas........... 103 1.1.1IdentiÞcationofcubicvelocityerrors....................... 103 1.1.2Eliminationofcubicvelocityerrorsviacorrectionstomomentequilibria... 108 1.1.3Strainratetensorcomponentsvianon-equilibriummoments.......... 113 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

PAGE 9

ix B 115 2.1OrthogonalMatrixoftheMomentBasis K fortheD3Q19Model........... 115 2.2RawMomentsofSourcesfortheD3Q19Model..................... 116 2.3ProjectionsoftheRawSourceMomentstotheOrthogonalBasisVectorsforthe D3Q19Model........................................ 117 2.4SourceTermsinParticleVelocitySpacefortheD3Q19Model............. 115 2.5Thechangeofdi!erentmomentsundercollisionfortheD3Q19lattice........ 116 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

PAGE 10

x LISTOFTABLES TABLE 2.1Dependenceofthetotaldragcoe"cient CD ontheReynoldsnumberReandthe powerlawindex n duetotheßowofpowerlawßuidspastanunconÞnedcirculatorresultsobtainedusingthecascadedLBmethodfordi!erentgridresolutions withthebenchmarkdataofBhartietal(2006)...................... 45 2.2Dependenceofthetotaldragcoe"cients C d 1 and C d 2 onthegapratio G andthe powerlawindex n for Re =5duetotheßowofpowerlawßuidspastanunconÞnedpairofcircularcylindersinatandemarrangement.Comparisonbetween thecomputedresultsobtainedusingthecascadedLBmethodwiththebenchmark dataofPatiletal[1].................................... 53 2.3Dependenceofthetotaldragcoe"cients C d 1 and C d 2 onthegapratio G andthe powerlawindex n for Re =10duetoßowofpowerlawßuidspastanunconÞned pairofcircularcylindersinatandemarrangement.ComparisonbetweenthecomputedresultsobtainedusingthecascadedLBmethodwiththebenchmarkdataof Patiletal[1]......................................... 53 2.4Dependenceofthetotaldragcoe"cients C d 1 and C d 2 onthegapratioGandthe powerlawindex n for Re =20duetotheßowofpowerlawßuidspastaunconÞnedpairofcircularcylindersinatandemarrangement.Comparisonbetween thecomputedresultsobtainedusingthecascadedLBmethodwiththebenchmark dataofPatiletal[1].................................... 54

PAGE 11

xi LISTOFFIGURES FIGURE 2.1Schematicsofthenon-Newtonianbenchmarkßowproblems.(a)Pressure-driven ßowbetweentwoparallelplates,(b)lid-drivencavityßow,(c)cross-ßowovera circularcylinder,and(d)ßowoverapairofcircularcylindersintandemarrangement............................................. 26 2.2ComparisonofthenormalizedvelocityproÞlesofpowerlawßuidsinachannel computedusingthecascadedLBmethodandtheanalyticalsolutionfor n =0 . 5 , 1 . 0 and1 . 5at Re =100..................................... 28 2.3ComparisonofthenormalizedvelocityproÞlesofBinghamßuidsinachannel computedusingthecascadedLBmethodandtheanalyticalsolutionfor ! B = 0 , 4 ! 10 ! 4 and8 ! 10 ! 4 atRe=100........................... 29 2.4Relativeglobalerroragainstthevariationinthegridspacing,plottedinlog-log scales,forthesimulationofpowerlawßuidßowinachannelat n =0 . 5and1 . 0 obtainedusingthecascadedLBmethod.......................... 30 2.5Computedvelocitycomponent u/U p alongtheverticalcenterlineandcomponent v/U p alongthehorizontalcenterlineobtainedusingthecascadedLBmethod(lines) andcomparedwiththebenchmarksolutionofNeofytou[2](symbols)forlid-driven cavityßowofpowerlawßuids.Reynoldsnumber Re =100andthepowerlaw index n =0 . 5 , 1 . 0 , 1 . 5.................................... 32 2.6Computedvelocitycomponent u/U p alongtheverticalcenterlineandcomponent v/U p alongthehorizontalcenterlineobtainedusingthecascadedLBmethod(lines) andcomparedwiththebenchmarksolutionsofNeofytou[2](symbols)forliddrivencavityßowofBinghamßuids.Reynoldsnumber Re =100andBingham number Bn =0 . 1 , 1 . 0 . ................................... 33

PAGE 12

xii 2.7StreamlinecontoursforliddrivencavityßowatÞxed n =0 . 5and(a)Re=100, (b)Re=400,(c)Re=1000computedusingthecascadedLBmethod......... 35 2.8StreamlinecontoursforliddrivencavityßowatÞxed n =1 . 0and(a)Re=100, (b)Re=400,(c)Re=1000computedusingthecascadedLBmethod......... 36 2.9StreamlinecontoursforliddrivencavityßowatÞxed n =1 . 5and(a)Re=100, (b)Re=400,(c)Re=1000computedusingthecascadedLBmethod......... 37 2.10StreamlinecontoursforliddrivencavityßowofpowerlawßuidsatÞxedRe=400 (a) n =0 . 75,(b) n =1 . 0,(c) n =1 . 25computedusingthecascadedLBmethod.... 37 2.11StreamlinecontoursforliddrivencavityßowofBinghamßuidsatÞxed Re =400 and(a) Bn =0,(b) Bn =0 . 01,(c) Bn =0 . 1computedusingthecascadedLB method............................................ 38 2.12Locationoftheprimaryvortexforlid-drivencavityßowat n =1.Re=100 , 400 and1000computedusingthecascadedLBmethod................... 38 2.13Locationoftheprimaryvortexforlid-drivencavityßowat n =1,Re=100 , 400 and1000computedusingthecascadedLBmethod................... 39 2.14locationoftheprimaryvortexforlid-drivencavityßowatRe=400for(a)power lawßuidsfor n =0 . 75 , 1 . 0 , 1 . 25 , (b)Binghamßuidsfor Bn =0 , 0 . 1 , 1 . 0computed usingthecascadedLBmethod............................... 40 2.15Streamlinecontoursoftheßowofashearthinningßuid( n =0 . 6)atdi!erentvaluesReynoldsnumberReof5 , 10 , 20 , 20and30overanunconÞnedcircularcylinder computedusingthecascadedLBmethod........................ 42 2.16Instantaneousstreamlineandvelocitycontoursoftheßowofapowerlawßuidat Re=100anddi!erentvalueofthepowerlawindex n of0 . 6 , 1 . 0and1 . 4overan unconÞnedcircularcylindercomputedusingthecascadedLBmethod......... 43 2.17Streamlinecontoursatdi!erenttimeinstantsoftheßowofapowerlawßuidat Re =100, n =1 . 4overanunconÞnedcircularcylindercomputedusingthecascadedLBmethod...................................... 44

PAGE 13

xiii 2.18ComparisonbetweenthenumericalresultsofthecascadedLBmethodandbenchmarkdata[3]forthetotaldragcoe"cient C D asafunctionofReynoldsnumber Re atdi!erentvaluesofthepowerlawindex n ...................... 46 2.19StreamlinecontoursfortheßowofpowerlawßuidsatRe=5atdi!erentvaluesof thegapratio G =2and4andofthepowerlawindex n =0 . 4 , 1 . 0and1 . 8ineach caseoveranunconÞnedpairofcylindersintandemarrangementcomputedusing thecascadedLBmethod.................................. 48 2.20Streamlinecontoursfortheßowofpowerlawßuidsat Re =10atdi!erentvalues ofthegapratio G =2and4andofthepowerlawindex n =0 . 4 , 1 . 0and1 . 8in eachcaseoveranunconÞnedpairofcylindersintandemarrangementcomputed usingthecascadedLBmethod.............................. 49 2.21StreamlinecontoursfortheßowofpowerlawßuidsatRe=20atdi!erentvaluesof thegapratio G =2and4andofthepowerlawindex n =0 . 4 , 1 . 0and1 . 8ineach caseoveranunconÞnedpairofcylindersintandemarrangementcomputedusing thecascadedLBmethod.................................. 51 2.22Instantaneousstreamlinecontoursfortheßowofpowerlawßuidsat Re =150 atdi!erentvaluesofthegapratio G =2and4andofthepowerlawindex n = 0 . 4 , 1 . 0and1 . 8ineachcaseoveranunconÞnedpairofcylindersintandemarrangementcomputedusingthecascadedLBmethod................... 52 2.23ComparisonofthemaximumReynoldsnumberfornumericalstabilityofSRTLB scheme,cascadedLBschemewithoutGIcorrections(centralmomentI),andcascadedLBschemewithGIcorrections(centralmomentsII)forsimulationofnonNewtonianliddrivencavityßowwithn=0.5,1.0,1.5.................. 56 3.1Three-dimensional,nineteenparticlevelocity(D3Q19)lattice.............. 85

PAGE 14

xiv 3.2Relativeglobalerroragainstthevariationinthegridspacing,plottedinlog-log scales,forthesimulationofpowerlawßuidßowinachannelat n =0 . 8, n =1 . 0 and1 . 5obtainedusingthecascadedLBmethod..................... 85 3.3Schematicsofthenon-Newtonianbenchmarkßowproblems.(a)ßuidßowsina lid-drivencavity,(b)3Dfullydevelopedßowinasquareduct............. 86 3.4comparisonofthenormalizedvelocityproÞlesofpowerlawßuidsinachannel computedusing3DCascaded-LBMandtheanalyticalsolutionfor n =0 . 8 , 1 . 0 , and1 . 5.atRe=100.................................... 87 3.5Flowthroughasquareductwithsidelength2asubjectedtoaconstantbody force:ComparisonofsurfacecontoursofthevelocityÞeldforReynoldsnumberRe =20(a)computedbytheD3Q15formulationofthecentralmomentLBMwith forcingtermwith(b)analyticalsolution(seeEq.3.4.1.2)................ 88 3.6Flowthroughasquareductwithsidelength2asubjectedtoaconstantbodyforce: ComparisonofvelocityproÞlescomputedbytheD3Q15formulationofthecentralmomentLBMwithforcingterm(symbols)withanalyticalsolution(lines)(see Eq.(3.4.1.2))atdi!erentlocationsintheductcross-sectionforReynoldsnumber Re=20........................................... 89 3.7StreamlinecontoursoftheductvelocityÞeldforReynoldsnumberRe=20.... 89 3.8Flowthroughasquareductwithsidelength2asubjectedtoaconstantbodyforce: ComparisonofvelocityproÞlescomputedbytheD3Q19formulationofthecentralmomentLBMwithforcingterm(symbols)withanalyticalsolution(lines)(see Eq.(3.4.1.2))atdi!erentlocationsintheductcross-sectionforReynoldsnumber Re=20........................................... 90 3.9StreamlinecontoursoftheductvelocityÞeldforReynoldsnumberRe=80.... 90

PAGE 15

xv 3.10Computedvelocitycomponent u/U p alongtheverticalcenterlineandcomponent v/U p alongthehorizontalcenterlineobtainedusingthecascadedLBmethod( lines)andcomparedwiththebenchmarksolutionofJin(2017)(symbols)forliddrivencavityßowofpowerlawßuids.Reynoldsnumber Re =100 , 400 , 1000and thepowerlawindex n =1 . 5................................. 91 3.11Computedvelocitycomponent u/U p alongtheverticalcenterlineandcomponent v/U p alongthehorizontalcenterlineobtainedusingthecascadedLBmethod( lines)andcomparedwiththebenchmarksolutionofJin(2017)(symbols)forliddrivencavityßowofpowerlawßuids.Reynoldsnumber Re =100 , 400 , 1000and thepowerlawindex n =1 . 0................................. 92 3.12Computedvelocitycomponent u/U p alongtheverticalcenterlineandcomponent v/U p alongthehorizontalcenterlineobtainedusingthecascadedLBmethod( lines)andcomparedwiththebenchmarksolutionofJin(2017)(symbols)forliddrivencavityßowofpowerlawßuids.Reynoldsnumber Re =100 , 400 , 1000and thepowerlawindex n =0 . 8................................. 93 3.13ComparisonofthemaximumReynoldsnumberfornumericalstabilityofSRTLB scheme,andcascadedLBschemeforsimulationofnon-NewtonianLiddrivencavity ßowwith n =0 . 8 , 1 . 0;1 . 5................................. 94

PAGE 16

1 CHAPTERI INTRODUCTION 1.1 Motivation Non-Newtonianßuidßowswithnonlinearrheologicalbehavioroccurinanumberofapplications suchasinchemical,biologicalandmaterialsprocessingcontexts.ItsbehaviorisusuallydeÞned asmaterialsforwhichtherateofshearisdependentonlyonthecurrentmagnitudeoftheshear stress.Suchmaterialscanbecharacterizedaspurelyviscousandtime-independentorgeneralized Newtonianßuids(GNF)[4].Therearemanydi!erentcharacteristicsthata!ectthebehaviorofa ßuidandthatmayalsobeusedtoprovideacategorizationofalltheseßuids.Ofprimaryimportanceindeterminingtheßowbehaviorofaßuidistheviscosity.Onthebasisofviscosity,ßuids canberoughlysubdividedintoNewtonianandnon-Newtonian.ForNewtonianßuids,theviscosityistakentobeconstantoversomerangeofshearrate(thoughitmayvaryslightlywithtemperatureandpressure).Fornon-Newtonianßuids,theviscositydependsontheshearrateitself. Therearemanyviscositymodelsorconstitutiverelationstochoosefromandtheyvaryintheir complexity.Someexamplesofconstitutiverelationsaregivenbythepowerlawmodel,Binghamßuidmodel,Herschel-Bulkleymodel,CrossmodelandCarreaumodel.Agoodoverviewof someofthemorewidelyusedviscositymodelscanbefoundinBarnesetal.[5]orBirdetal.[6]. Studyofrheologicalßowsusinganalyticaltechniquesisextremelychallengingduetotheirinherentnonlinearityandsolutionsfromsuchapproachesaregenerallylimitedtohighlyidealizedconditions.Thus,computationalmethodsplayanessentialroleforinvestigationsofnon-Newtonian ßows.Inthepast,Þnite-di!erence,Þnite-elementorÞnite-volumemethodsweredevelopedinthis regard.However,obtainingstableandaccuratesolutionsofnon-Newtonianßowsespeciallyon complexgeometriesandinthree-dimensionsisstillverychallengingandcomputationallyvery intensive.

PAGE 17

2 ThisthesisfocusesonfurtherdevelopingandapplyingadvancedformulationsofthelatticeBoltzmannmethod(LBM)fornon-Newtonianßows.TheLBMisanewclassofcomputationalßuid dynamics(CFD)techniquewhichisemployedtosimulateavarietyofßuidmechanicsproblems andhasaclearandrigorousfoundationbasedonkinetictheory.ItinvolvesthelatticeBoltzmannequation(LBE),whichcanbederivedasamuchsimpliÞedformoftheBoltzmannequation.Thisequationdescribesthestatisticaldistributionsofparticlesinßuidsatamesoscopic scale.Thisscaleisintermediatebetweenthemicroscopic(ormolecular)andthemacroscopic scales.ThedistributionofapopulationofparticlesaremodiÞedbythemovementorthestreamingandcollisionbehaviorofparticlesinaßowthatoccuronalattice.Atmacroscopicscales, thecollectivebehavioroftheparticlepopulationscorrespondstothemacroscopicdynamicsofa ßuidthatcanbewelldescribedbytheNavier-Stokesequations.TheLBMhasseveralcomputationalandphysicaladvantages.ThestreamingstepoftheLBMislinearandallthenonlinearity isrepresentedlocallyinthecollisionstepviaanequilibriumdistributionfunction;ontheother hand,thenonlinearitypresentintheconvectiontermoftheNavier-Stokesequationsisnonlocal,whichpresentschallengesinthedevelopmentoftheirnumericalsolutionschemes.Assuch, intraditionalCFDmethods,thepressureÞeldisobtainedbythesolutionofatimeconsuming elliptic(Poisson-type)equation,whichisavoidedintheLBM,whereitisobtainedlocally.The streamingstepoftheLBMisbasedonperfect-shiftadvection,which,whencombinedwitharelaxationmodelforthecollisionstep,canleadtosolutionsthataresecondorderaccuratewith relativelylownumericaldissipation,whencomparedtoconventionalschemesbasedonequivalent orderofaccuracy.ThecollisionstepoftheLBMisbasedonakineticrelaxationmodel,which canbedesignedtorepresentcomplexßuidßowphysics;inaddition,theycanbetailoredtoimprovenumericalstabilityoftheLBMforsimulationsofßowsencompassingwiderrangeofgoverningparameters.Furthermore,boundaryconditionsforcomplexgeometriescanberepresented atrelativeeaseusingtheLBMbyemployingsimplekineticrulesfortheparticlepopulations.Finally,sincetheLBMisbasedonlocalcomputations,itisnaturallyamenableforimplementation onparallelcomputersthatenablee"cientsolutionsforlargesizeßowproblems.

PAGE 18

3 However,theLBMalsohassomelimitations.Somemajordisadvantagesincludethefollowing. Itischallengingtosimulatehighspeedcompressibleßowsonstandardlattices,astheydonot haveenoughparticlespeedstorepresenttheenergyequationusingasingledistributionfunction. Asaresult,somemajormodiÞcationstotheLBM,e.g.,usinglargerparticlevelocitysets,are neededtosimulatecompressibleßows,whichmayunderminethesimplicityande"ciencyofthe approach.Also,themethodencountersdi"cultiesandissubjecttonumericalinstabilitieswhen thereisstrongcouplingbetweenßuiddensityandtemperaturevariations.Inaddition,whenthe goalistorepresenttrulyincompressibleßow,theLBMbecomesnumericallysti!toresolvethe propagationofacousticwavesatverylowMachnumberswithrelativelysmallsteps.TheLBM is,byconstruction,fullyexplicitintimeandmarchesthroughthetransientstorepresentsteadystateßows,whichcanbetimeconsuming.TheschemeisdesignedonuniformCartesianmeshes, anditischallengingtoextenditforresolutionofproblemsrequiringstretchedgrids.Moreover, theLBMcomputationsarememory-intensivewiththestreamingofparticlepopulationsincurringlargeamountofmemoryfetchoperations.Overcomingtheseandvariousotherlimitations arepartofcurrentresearchinthisÞeld. WewillnowbrießyreviewthecomputationalelementsoftheLBM.Thelatticeiscomprisedof discreteparticledirections " e ! ,andalongeachofthesedirections,weassociateadiscreteparticle distributionfunction f ! .Inthesimplestoftheformulation,referredtoasthesinglerelaxation timemodel,thecollisionstepoftheLBMisrepresentedbymeansofarelaxationofthisdiscrete distributionfunctiontoitslocalequilibrium f eq ! overacharacteristictime # ! ( f ! ( " x,t ))=( f eq ! " f ! ) / !,f eq ! = f eq ! ( # ,"u )(1.1) Thischaracteristictimecanberelatedtotherateofmomentumdi!usion,i.e.,themomentum di!usivityoftheßuid.Thepost-collisiondistributionfunction ! f ! canthenberepresentedas ! f ! ( " x,t )= f ! ( " x,t )+ # ! ( " x,t )(1.2) Then,theparticlepopulationsadvectalongthecorrespondingparticledirectionsinperfectshift

PAGE 19

4 manneroveratimestep,whichisreferredtoasthestreamingstep.Thiscanberepresentedas f ! ( " x + " e ! $ t ,t + $ t )= ! f ! ( " x,t )(1.3) Atthecompletionofthecollide-and-streamsteps,themacroscopicvariablesorthehydrodynamicÞeldssuchastheßuiddensityandthemomentumareobtainedaszerothandÞrstmomentsoftheparticledistributionfunction,whilethepressureisrelatedtothedensityviaan equationofstateasfollows: # = " ! f ! , #"u = " ! f ! " e ! ,p = c 2 s # (1.4) Neartheboundaries,theboundaryconditionsonthemacroscopicvariablesaree!ectivelyrepresentedbymeansoftheabovepicturebasedonastatisticalformulationoftheparticledynamics,suchastheso-calledbouncebackschemetorepresenttheno-slipboundarycondition.A detailedreviewofvariouscollisionmodelsandboundaryconditionsarepresentedintherecent monographbyRef.[7]. 1.2 ParallelComputing Thecomputationalcostofsolvinglargesizeproblemsisveryhighanditdemandshigh-speed processorsandlargememoryinthesystem.Themodernmulti-coreprocessorshavetheabilityto performseveralcomputationsinparallel,knownasparallelprocessing,whichcanbeusedtoreducethecomputationtimedramatically.Aserialprocessorsysteminonehandsolvesaproblem sequentiallywhileasystemwithparallelprocessorscansolvemultiplepartsofalargeproblem simultaneouslyandcombinetheresultbycommunicatingandcooperatingwitheachotherto producetheÞnalresult.ParallelprocessingisanimportantdevelopmentincomputingtechnologyandcurrentlybeingusedliterallyineveryarenaofscientiÞcapplication.Withtheadventof parallelprocessing,wenotonlycansolveacomputationallyintensivecomplexthree-dimensional problemsfasterbutalsogainedthepowertoincreasethegridresolutionofthespatialandtemporaldimensions.ParallelcomputerscanbeclassiÞedinmanywaysbasedonthedataandin-

PAGE 20

5 structionstreams,memoryallocation,e.g.,theprocessorshavingseparatememoryoroneshared memory.CommonlytheclassiÞcationofparallelcomputersisknownasFlynnsclassiÞcation. 1.2.1 ClassiÞcationofParallelPlatforms 1.2.1.1 BasedonMemoryDistribution • SharedMemoryMachines:Acomputersystemwithsharedmemoryhavelargeglobal memoryandasmallnumberofprocessorsareconnectedtothememorymodulesby meansofsomekindofinterconnectionnetwork.Alltheprocessorsaccesstheglobalmemoryinthesamewayandthatswhyitisalsoknownasauniformmemoryaccess(UMA) machine.Theconvenienceinprogrammingtheprocessinsuchamachineisoneofthe mostimportantadvantagesofsharedmemorymachines.ThemostimportantdisadvantageofthismachineisthatitcanthandlealargenumberofprocessorsbecauseitisdifÞculttoscalethecentralizedmemoryandtheinter-connectionnetworktogether.Most machineusingsharedmemoryarchitecturehaveatmost64processors.SGIPowerOnyx, SUNmultiprocessors,andIBMSP2highnodearesomeexampleofsharedmemorymachines. • DistributedMemoryMachines:Inthedistributedmemorysettingeachprocessorhasits ownlocalmemory.Duringaparallelcomputation,whenoneprocessorneedsinformation whichisnotinitsownlocalmemoryitneedstocommunicatewithotherprocessorcontainingthedata.Thesemachinescanhavemanymoreprocessorthanasharedmemory machinedoes.Adistributedmemorymachinecanhaveseveralhundredoreventhousandsofprocessors.IBMSP2,IntelParagonandNetworkofWorkstations(NOWs)or Clustersaresomeexamplesofsuchmachines.

PAGE 21

6 1.2.1.2 BasedonInstructionandDataStreams • SingleInstructionandSingleDataStream(SISD):Theconventionalserialcomputers belongstothiscategory.Theyhavethesingleprocessorandsinglememoryandoperate ondatafrommemoryoneatatime. • SingleInstructionandMultipleDataStream(SIMD):Insuchsystems,multipleprocessingelementsworkunderthecontrolofasinglecontrolunit.Asingleinstructionisexecutedbymultipleprocessorsondi!erentdatasets.SIMDcanbeeitherasharedmemory machineoradistributedmemorymachine.Everyprocessormustbeallowedtocomplete itsinstructionbeforethenextinstructionistakenforexecution.ExamplesofSIMDincludesILLIAC-IV,PEPE,BSPetc. • MultipleInstructionandMultipleDataStream(MIMD):ThisorganizationismoreßexiblebecausetheMISDmachinescanexecutemultipleinstructionsondi!erentprocessors workingondi!erentdatasets.Intherealsense,MIMDistheorganizationthatissaidto betheparallelcomputer.AllmultiprocessorssystemsfallundertheclassiÞcation.ExampleofMIMDorganizationincludesC.mmp,BurroughsD825,Cray-2,S1,CrayX-MP etc. 1.2.2 ParallelizationofLBMalgorithms Asdiscussedearlier,thestandardformulationoftheLBMconsistsoftwosteps,i.e.,thecollisionstepandthestreamingstep.Thecollisionstepiscompletelylocalasalltherequiredinformationforperforminganycollisionmodelareavailableatagivenlatticenode.Thestreaming stepinvolvesperfectshiftadvectiontotheneighboringnodesalongthelatticelinks,requiring onlylocalizedcommunicationinthecontextofparallelcomputing.TheLBMalgorithmscanbe readilyparallelizedusingthestandarddomaindecompositiontechnique.Thatis,partsofeach

PAGE 22

7 ofthecomputationaldomainisallocatedtoadi!erentprocessor,eachofwhichperformsthe collisionstepsimultaneouslyandthencommunicatesthedataoftheparticledistributionfunctioninthevicinityoftheneighboringprocessorbeforeperformingthestreamingstep.While thefocusofthisdissertationonfurtherdevelopingandapplyingadvancedLBMformulations fornon-Newtonianßowsandnotonparallelization,thereisanextensiveliteraturethatdemonstratenear-linearparallelscalingoftheLBM[7].ParallelizationofLBMcanbeimplemented usingprogramminglanguagelibrariessuchastheMPI,OpenMPorCUDA.Assuch,thenew methodspresentedinthisresearchcanbeusedtosolvelargesizenon-Newtonianßowproblems moree"cientlyusingparallelcomputersinfuturestudies. 1.3 Focusofthiswork Asreviewedanddiscussedinthebodyofthisthesis,theexistingliteratureontheLBMfornonNewtonianßowsistypicallybasedonlessrobustcollisionmodelsthataregenerallyfocusedon two-dimensions(2D),andtheirapplicationsto3Drheologicalßowsarerelativelyrare.Wewill exploittheadditionaldegreesoffreedomavailableindesigningthecollisionsteptoconstruct morerobustLBMfornon-Newtonianßowsthanthoseavailableintheliterature.ThespeciÞc objectivesofthisthesisandnewcontributionsarisingfromthisworkareasfollows. (i) WewilldevelopanadvancedcascadedLBMbasedoncentralmomentsandmultiplerelaxationtimesfor2Dnon-Newtonianßows;inthisregard,wewillconsideratwo-dimensional, ninevelocity(D2Q9)latticeforrepresentingpowerlaw(modelingshearthinningorshear thickeningbehavior)andBinghamßuid(representingyieldbehavior)modelsastherepresentativeconstitutiverelations. (ii) WewillconstructanadvancedcascadedLBMbasedoncentralmomentsandmultiple relaxationtimesfor3Dnon-Newtonianßows;inthisregard,wewillemployathreedimensional,nineteenvelocity(D3Q19)latticeforsimulationofpowerlawßuidßows.

PAGE 23

8 (iii) Wewillvalidatethenew2Dand3DLBMsmentionedabovetodemonstratetheiraccuracyforavarietyofcomplexnon-Newtonianßowsbycomparisonsagainstanalytical/numericalbenchmarksolutions,theirorderofgridconvergenceandimprovements innumericalstabilityinaccessingtheßowregimeswithbroaderrangeofcharacteristic parameterswhencomparedtootherformulationsofthecollisionsmodelsfortheLBM. (iv) Finally,wewillstudythephysicsofcomplexnon-Newtonianßowofpowerlawßuidsover apairofcylindersintandemarrangement.Inthisregard,wewillinvestigatethewake interactionsandtheire!ectsontheresultingforces,orthedragcoe"cients,oneach cylinder,andelucidatetheinßuenceofthevariouscharacteristicparameters,suchasthe powerlawindexandthegapratiobetweenthecylinders.Suchprototypicalßowsareof interestinvariousengineeringapplications. 1.4 OutlineofThesis Theorganizationofthisthesis,whichisdividedintofourchapters,isasfollows.Followingthis introduction(Chapter1),inChapter2,wewillpresentthe2DcascadedLBMbasedonthecentralmomentsfornon-Newtonianpowerlawßuids.Inthischapter,wewillpresentaliterature reviewoftheexistingLBMfornon-Newtonianßows,discussthederivationofthe2Dcascaded LBMtorepresentthepowerlawandBinghamßuids,andaccompaniedbyaChapman-Enskog analysistorecovertheemergentNavier-Stokesequationswithnonlinearconstitutiverelations, andperformnumericalvalidationandcomparativestudies.Furthermore,wewillstudythephysics ofthenon-Newtonianßowoverasingleandapairofcylindersintandemarrangementunder variousßowconditionsandcharacteristicparameters.Chapter3willpresentadetailedderivationofthenew3DcascadedLBMfornon-NewtonianßuidsusingaD3Q19latticeandanumericalvalidationstudyforvariousbenchmarkproblems,includingcomparisonsagainstarecent numericalsolutionforsheardrivennon-Newtonianßowinacubiccavity.Furthermore,wewill

PAGE 24

9 demonstratetheimprovementsinnumericalstabilityachievedwiththeuseofour3Dformulationwhencomparedtootherexistingcollisionmodels.Finally,asummaryandconclusionsof thisworkarepresentedinChapter4,alongwithindicationsofthepossibledirectionoffuture researchinthisarea.

PAGE 25

10 CHAPTERII TWO-DIMENSIONALCASCADEDLBMFORNON-NEWTONUANFLUID FLOWS 2.1 Introduction Non-Newtonianßowsofcomplexßuidsariseinabroadrangeofsituationsinnaturalsettings andengineeringapplicationsincludingmaterials,chemical,petroleum,food,hydrologyandgeologicalprocesses.Theirrheologicalbehaviorisgenerallycharacterizedbynonlinearconstitutiverelations[8,9].Analyticalsolutionsareextremelychallengingtoobtainandareoftenlimitedtohighlyidealizedconditionsowingtotheirnonlinearity.Inthisregard,withadvancesin computingcapabilities,numericalmethodsthatprovidestableandaccuratesolutionsofnonNewtonianßuidßowsespeciallyincomplexgeometriesplayanimportantrole.Inthepast,numericalschemesbasedonÞnitedi!erenceorÞniteelementtechniqueshavebeenappliedand studiedforvariousrheologicalßuidßows[10,11].Morerecently,apowerfultechniquearising fromadi!erentperspective,viz.,thelatticeBoltzmann(LB)method,hasshownremarkablepotentialtotacklesuchcomplexßuidmechanicsproblems. ThelatticeBoltzmann(LB)methodisanovelemergingapproachincomputationalßuiddynamicsbasedonkinetictheory.Inparticular,itmaybeobtainedasaspecialdiscretizationofthe Boltzmannequationandcanbecharacterizedasamesoscopicmethodinvolvingthedynamics ofthedistributionofparticlepopulationsthatrespectlocalconservationanddiscretesymmetryprinciples.Themotionoftheparticledistributionfunctionsisrepresentedasaperfectshift streamingprocessalongthelatticedirections,whichisfollowedbyalocalcollisionstepmodeled asarelaxationprocesstocertainequilibriainLBformulations.Thebehaviorofthemacroscopic ßuidmotionrepresentedbythehydrodynamicvariablescanthenbeobtainedviatakingtheki-

PAGE 26

11 neticmomentsoftheparticledistributionfunctions.TheconsistencyoftheLBschemestothe Navier-StokesequationsmaybeanalyticallyshownbyaChapman-EnskogexpansionoraTaylor seriesexpansionwithappropriatelyprescribedscalingrelationshipbetweenthediscretizedspace stepandthetimestep.Algorithmically,duetoitsinherentlyparallelizationcapabilities,appealingphysicalfeatureasanaturalframeworktoutilizekineticmodelsforcomplexßows,andits easeforrepresentationofboundaryconditionsinheritedfromitsparticlebaseddescriptionhas madeitauniqueandpowerfultechniqueforßowsimulations.Indeed,theadvancesintheLB methodsandtheirapplicationstoavarietyofcomplexßowshavebeendiscussedinvariousreviews[12,13,14,15]andmonographs[16,17,7]. ThecollisionstepintheLBmethod,whoserelaxationmodelcanbetunedtorepresentavariety ofphysicsassociatedwiththemotionofcomplexßuids,alsohasapredominantroleinitsnumericalstability.Oneoftheearliestcollisionmodelsisthesinglerelaxationtime(SRT)model[18], whichispopularowingtoitssimplicityandcaseofimplementation.However,itisknowntobe susceptibletonumericalinstabilities,whenrequiredtoencompasswiderrangeofvariationsin ßuidproperties,suchastheviscosity,orforsimulationsofhighReynoldsnumberßows.Onthe otherhand,byconsideringthemultiplerelaxationtime(MRT)model[19],inwhichdi!erentraw momentsofvariousordersrelaxatdi!erentratesinthecollisionstep,asigniÞcantimprovement inthenumericalstabilityoftheLBmethodcanbeachieved.Furtherenhancementsintheßow modelingcapabilitiesandimprovementsinnumericalstabilizationaremadebypossiblebythe introductionofthecascadedLBmethod[20]. ThecascadedLBmethodisamulti-parametricscheme,whichisbasedonconsideringrelaxation intermsofcentralmoments,whichareobtainedbyshiftingtheparticlevelocitybythelocal ßuidvelocity,intheconstructionofthecollisionoperator.Sincetheyareconstructedforamovingframeofreference,therelaxationprocessatsuccessivelyhigherordersrepresentacascaded structure.Itcanberecastasalinearrelaxationtoageneralizedequilibriumintherestorlattice frameofreference[21].Inordertoextenditscapabilitiesforpracticalapplications,forcingterms

PAGE 27

12 basedoncentralmomentstoincorporatetheinßuenceofimpressedforcesontheßuidmotion wereintroducedintwo-dimensions(2D)[22]andthreedimensions(3D)[23].Detailedstudies onthenumericalpropertiesofsuchadvancedcollisionmodelsandimprovementsachievedwere demonstratedin[24,25].Avariationofthisapproach,inwhichadiscreteratherthanacontinuousequilibriumisadoptedforcentralmomentsbasedrelaxationispresentedin[26].Recently, anapproachtoacceleratetheconvergenceofthecascadedLBmethodhasbeenconstructedand studiedin[27]andfurtherextendedtotacklethermalßowsin2D[28]and3D[29].Itmaybe notedsincethecentralmomentsareobtainedbyshiftingtheparticlevelocitywiththelocalßuid velocity,itprovidesanaturalsettingtomaintainGalileaninvarianceofthechosensetofmoments,andprovidesadditionalfreedomtotimethehigherorderrelaxationtimes,bothofwhich canenhancethequalityofsolutionsandthenumericalstabilityinsimulations.Indeed,recently animprovedpreconditionedcascadedLBformulationwithoutcubicvelocityerrorsonastandard latticewasconstructedandanalyzed[30].Todate,thecascadedLBmethodhasonlybeenappliedandstudiedforNewtonianßuidßowsintheliterature. OneoftheearliestapplicationsoftheLBmethodsfornon-Newtonianßowswaspresentedin[31], whichusedaSRTmodeltorepresentpowerlawßuids.Simulationsofviscoplasticßowswithfree surfaceutilizingaregularizedBinghammodelinaLBformulationwasdiscussedin[32].Subsequently,Boeketal.[33]investigatedpower-lawbasedßuidßowthroughaporousmediabyconsideringaLBschemeadoptingasinglerelaxationtimeapproximation.TheuseofotherrheologicalmodelssuchastheCrossmodelofviscosityintheLBframeworkwasconsideredbyKehrwald etal.[34]andthetruncatedpower-lawmodelforshearthinningandshearthickeningßuidswith improvedaccuracyandstabilitybyGabbanellietal.[35].Asecond-orderspatialaccuracyofthe LBmethodtosimulatenon-Newtonianßuidßowswasdemonstratedin[36].Alatticekinetic formulationwithavariableparameterintroducedtodeterminetheshear-dependentviscosity oftheßuidwasproposedin[37],whichÞxedtherelaxationtimetounityandemployedaÞnite di!erenceapproximationtoevaluatethestrainratetensortoimprovestability.Animplicitfor-

PAGE 28

13 mulationoftheLBcollisionmodeltorepresentyield-stressßuidssuchastheBinghamßuidwas presentedin[38].Thee!ectivenessoftheLBapproachtorepresentgeneralizedNewtonianßuidswasstudiedin[39].Adetailedinvestigationofthenon-Newtonianßowsinporousmediawas presentedin2Din[40]andin3Din[41].FurtheradvancesintheLBsimulationsofBingham plasticßuidßowswerereportedin[42,43,44]andbloodßows,forexample,in[43].Itcanbe seenthatmostoftheLBstudiesonnon-NewtonianßowsemployedtheSRTmodelforthecollisionstep,which,asdiscussedearlier,hassigniÞcantlimitationsrelatedtonumericalstability issues.Morerecently,Chaietal.[45]presentedaMRT-LBformulationbasedonrawmoments tostudygeneralizedNewtonianßuidswithimprovedstability.Also,noteworthyaretherecent accuracystudiesthatdeterminetheconditionsontheLBschemessoastomaintainsecond-order accuracyinsimulationsofnon-Newtonianßowsnumerically[46]andanalytically[47]. Inthiswork,wepresent,fortheÞrsttime,aninvestigationoftheadvancedcentralmoments basedcascadedLBmethodandusingmultiplerelaxationtimesforthesimulationofnon-Newtonian ßuids.Inthisregard,weconsidergeneralizedNewtonianßuidsrepresentedbyapower-lawbased constitutiverelationencompassingbothshearthinningandshearthickeningbehavior.Inaddition,wealsoconsideryield-stressßuidsbasedonaregularizedBinghammodel.Thecomputation ofthestrainratetensor,whichareneededinthesemodels,areperformedlocallyinthecascaded LBformulationviausingsecondordernon-equilibriummoments.Inthisapproach,weincorporateadditionalGalileaninvariance(GI)correctiontermstoeliminatecubicvelocityerrorsfor simulationofnon-Newtonianßowsonastandardtwo-dimensional,ninevelocity(D2Q9)lattice withhigherÞdelityandimprovedstability.WethenperformthevalidationofthecascadedLB methodforawiderangeofbenchmarkrheologicalßowproblems,includingthenon-Newtonian poiseuilleßow,lid-drivencavityßow,andßowoverasinglecylinderandoverapairofcylinders intandemarrangement,inordertodemonstrateitsaccuracy,second-ordergridconvergenceand improvementinnumericalstability.Inaddition,inthecaseofthenon-Newtonianßowovera pairofcylinders,westudythewakeinteractionsandtheire!ectsontheresultingforcesineach

PAGE 29

14 cylinderandelucidatetheinßuenceofvariouscharacteristicmodelßowparameters,includingthe powerlawindex.Thisthesisisorganizedasfollows.Inthenextsection(Sec.II),wepresentthe macroscopicgoverningequationsforthedynamicsofthenon-Newtonianßuidmotion,including theattendantconstitutiverelations.SectionIIIpresentsthecascadedLBmethodusingcentral momentsandmultiplerelaxationtimes,andincludingtheGIcorrectiontermsforthesimulation ofnon-Newtonianßows.Numericalresultsandtheirdiscussionforvariousbenchmarkßowproblemandaninvestigationofthenon-Newtonianßuidmechanicoftheßowoverapairofcylinders aregiveninSec.IV.SectionIVpresentsacomparativestudyofdi!erentcollisionmodelsinthe LBmethodtodemonstratestabilityimprovementwiththeuseofthecentralmomentsbased formulationfornon-Newtonianßow.Finally,inSecIV,summaryandconclusionsarisingfrom thisworkaregiven.Inaddition,aChapman-enskoganalysisofthecubicvelocityerrorsforthe D2Q9latticeandaderivationoftheresultingGIcorrectiontermsarepresentedintheappendix. 2.2 Macroscopicgoverningequationsforthedynamicsofnon-Newtonianßuidmotion InordertorepresentthemotionofthegeneralizedNewtonianßuids(GNF)in2D,wewritethe followingmathematicalmodelgivenintermsofthecontinuityandmomentumequations,respectively,as %# % t + % % x j ( # u j )=0 , (2.1a) % % t ( # u i )+ % % x j ( # u i u j )= " % P % x i + %& ij % x j + Fi, (2.1b) where # isthefuiddensity, u i isthevelocityÞeldwith i # ( x,y ), P isthepressure F i isthe cartesiancomponentoftheimposedbodyforce,and & ij istheviscousstresstensorinducedwithin theGNF,whichcanberepresentedas & ij = µ ( | ú ' ij | )ú ' ij . (2.2)

PAGE 30

15 Here,ú ' ij istheshearratetensorrelatedtothestrainratetensor S ij via ú ' ij =2 S ij ,S ij = 1 2 ( % j u i + % i u j ) , (2.3) and µ ( | ú ' ij | )istheshearratedependente!ectiveorapparentdynamicviscosityoftheGNF,where themagnitudeoftheshearrate | ú ' ij | isrelatedtothesecondinvariantofthesymmetricstrain ratetensor S ij as | ú ' ij | = # 2 S ij S ij . (2.4) Here,theusualconventionofthesummationoftherepeatedindicesisassumed.Inthisstudy, weconsidertheconstitutiverelationsforthepowerlawßuidandtheBinghamßuidasrepresentativeofaclassofGNFfortheirnumericalsolutionusingthecascadedLBmethoddiscussedin thenextsection.Theconstitutiverelationforthepowerlawßuidcanbewrittenas & ij = µ p | ú ' ij | n ! 1 ú ' ij , (2.5) wherethemodelparameter µ p istermedastheconsistencycoe"cientofthepowerlawßuidand theexponent n isdesignatedasthepowerlawindex.Basedonthevaluesofthislattermodelparameter,di!erenttypesofnon-Newtonianßuidscanberepresented: n< 1modelsshearthinningorpseudo-plasticßuid,while,bycontrast n> 1correspondstoshearthickeningßuids; when n =1,theabovedegeneratestotheNewtonianconstitutiverelation.ComparingEqs.(2.2) and(2.5),theshearratedependentapparentviscosityofpowerlawßuidscanberepresentedas µ power ( | ú ' ij | )= µ p | ú ' ij | n ! 1 (2.6) Ontheotherhand,inordertomodeltheviscoplasticbehaviorofyieldstressßuids,weconsider thefollowingconstitutiverelationfortheBinghamßuid & ij = $ & B | ú ' ij | + µ p % ú ' ij , | & ij | > & B , (2.7a) ú ' ij =0 , | & ij | < & B . (2.7b)

PAGE 31

16 Here, & B representstheyieldstressand µ B thepseudo-plasticviscosityoftheBinghamßuid.In general,theaboveconstitutiverelationexhibitsdiscontinuityat & B ,whichcanresultinnumericalstabilityissuesinsimulations.Inordertoovercomethischallenge,acertainformofsmoothingorregularizationcanbeapplied.Inthiswork,weconsiderapopularexponentialdecaybased regularizationduetoPapanastasiou,whichcanbewrittenas[48] & ij = & µ B + & B | ú ' ij | (1 " e ! m | ú " ij | ) ' ' ij . (2.8) Here, m isthestressgrowthexponentandEq.(2.8)isespeciallyagoodapproximationtothe constitutiverelationfortheBinghamßuid(Eq.(2.8))forlarge m .Clearly,forthespecialcase, where & B =0Eq.(2.8)modelsaNewtonianßuid.Finally,itfollowsfromEqs.(2.2)and(2.8) thattheapparentviscosityoftheregularizedBinghamßuidcanberepresentedas µ Bigham ( | ú ' ij | )= µ B + & B | ú ' ij | (1 " e ! m | ú " ij | ) . (2.9) 2.3 CascadedLBMbasedonCentralmomentsfornon-Newtonianßuidßow Inthissection,wewilldescribethecascadedLBformulationbasedoncentralmomentsandmultiplerelaxationtimestosolvethenon-Newtonianßuidßowmodelspresentedearlier.In2D,we considerthetwo-dimensional,ninevelocity(D2Q9)lattice,andstartwiththefollowingmoment basis: T = ( | # $ , | e ! x $ , | e ! y $ , | e 2 ! x + e 2 ! y $ , | e 2 ! x " e 2 ! y $ , | e ! x e ! y $ , | e 2 ! x e ! y $ , | e ! x e 2 ! y $ , | e 2 ! x e 2 ! y $ , (2.10)where | # $ = || " e ! | 0 $ =(1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1)   , | e ! x $ =(0 , 1 , 0 , " 1 , 0 , 1 , " 1 , " 1 , 1)   , | e ! y $ =(0 , 0 , 1 , 0 , " 1 , 1 , 1 , " 1 , " 1)   , | e 2 ! x + e 2 ! y $ =(0 , 1 , 1 , 1 , 1 , 2 , 2 , 2 , 2)   ,

PAGE 32

17 Theseninenon-orthogonalbasisvectorsareobtainedfromthemonomials e m ! x e n ! y atsuccessively increasingordersoftheexponents m and n toacompletesetassupportedbytheD2Q9lattice. Here,wehaveconsideredtheusualbra % . | and | . $ notationstorepresentthe9-dimensionalrow andcolumnvectors,respectively.Wethentransformtheabovebasisintoanequivalentsetoforthogonalbasis(Eq.(2.10))vectorsviatheGram-Schmidtprocedureandarrangedintheincreasingorderofmomentsasfollows: | K 0 $ = | # $ , | K 1 $ = | e ! x $ , | K 2 $ = | e ! y $ , | K 3 $ =3 | e 2 ! x + e 2 ! y $" 4 | # $ , | K 4 $ = | e 2 ! x " e 2 ! y $ , | K 5 $ = | e ! x e ! y $ , | K 6 $ = " 3 | e 2 ! x e ! y $ +2 | e ! y $ , | K 7 $ = " 3 | e ! x e 2 ! y $ +2 | e ! x $ , | K 8 $ =9 | e 2 ! x e 2 ! y $" 6 | e 2 ! x + e 2 ! y $ +4 | # $ . (2.11) Thesearegroupedtoobtainthefollowingorthogonaltransformationmatrix K ,whichisthen usedtoconstructthecascadedLBmethod: K = ( | K 0 $ , | K 1 $ , | K 2 $ , | K 3 $ , | K 4 $ , | K 5 $ , | K 6 $ , | K 7 $ , | K 8 $ ) ) (2.12) Itcanbeenumeratedas K = * + + + + + + + + + + + + + + + + + + + + + + + + + + + , 100 " 400004 110 " 11002 " 2 101 " 1 " 1020 " 2 1 " 10 " 1100 " 2 " 2 10 " 1 " 1 " 10 " 20 " 2 111201 " 1 " 11 1 " 1120 " 1 " 111 1 " 1 " 1201111 11 " 120 " 11 " 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . / (2.13) Thecollisionoperatorinthisapproachisbasedonrelaxationofcentralmomentstotheirequilibria.Thisisconstructedbasedonmatchingtheresultsobtainedforthecontinuouscentralmomentsfromkinetictheorytothoseforthediscretecentralmomentssupportedbythelattice.

PAGE 33

18 Hence,deÞningthecontinuousequilibriumcentralmomentsoforder( m + n )as 0 $ M x m y n = 1 " !" 1 " !" f M ( ( x " u x ) m ( ( y " u y ) n d ( x d ( y , whichyields | 0 $ M x m y n $ =( 0 $ M 0 , 0 $ M x , 0 $ M y , 0 $ M xx , 0 $ M yy , 0 $ M xy , 0 $ M xxy , 0 $ M xyy , 0 $ M xxyy )   =( # , 0 , 0 ,c 2 s # ,c 2 s # , 0 , 0 , 0 ,c 4 s # )   . (2.14) Here, f M = f M ( # ,"u, " ( )denotestheMaxwell-distributionfunctioninthecontinuousvelocity space " ( =( ( x , ( y ): f M ( # ,"u, " ( )= # 2 $ c 2 s exp 23 " ( % & ! % u ) 2 2 c 2 s 45 ,where c s isthespeedofsound,whichisa freeparameter,andwetypicallyset c 2 s =1 / 3. Inaddition,toimposethee!ectofanexternalforceÞeld " F =( F x ,F y ),wedeÞnethecontinuous centralmomentsofthesourcesoforder( m + n )as 0 % F x m y n = 1 " !" 1 " !" & f F ( ( x " u x ) m ( ( y " u y ) n d ( x d ( y , where& f F isthechangeinthedistributionfunctionduetotheforceÞeld.FromRef.[22],it followsthat | 0 % M x m y n $ =( 0 % M 0 , 0 % M x , 0 % M y , 0 % M xx , 0 % M yy , 0 % M xy , 0 % M xxy , 0 % M xyy , 0 % M xxyy )   =(0 ,F x ,F y , 0 , 0 , 0 , 0 , 0 , 0)   . (2.15) Basedontheabove,wenowpresentthecollisionoperatorandthesourcetermforthecascaded LBequation(LBE).Inthisregard,inordertomaintainasecondorderaccuracy,thetrapezoidal ruleisusedtorepresentthesourceterminthecascadedLBE,whichcanbewrittenas f ! ( " x + " e ! $ t ,t + $ t )= f ! ( " x,t )+ # ! ( % x,t ) + 1 2 ( S ! ( % x,t ) + S ! ( % x + % e ! ' t ,t + ' t ) ) . (2.16) Here, f ! isthediscretedistributionfunctionalongtheparticledirection ) : | f ! $ =( f 0 ,f 1 , ááá ,f 8 )   and S ! isthediscretesourcetermrepresentedas | S ! $ =( S 0 ,S 1 ,S 2 , ááá ,S 8 )   .thecollisionterm # ! canbewrittenas # ! = # ! ( " x,t )=( Ká 0 g ) ! , (2.17) where 0 g = | 0 g ! $ =( 0 g 0 , 0 g 1 , ááá , 0 g 8 )   isthevectorrepresentingthechangeinmomentsofdi!erent ordersundercollision,whichwillbeprescribedlatter.SinceEq.(2.16)issemi-implicit,byapply-

PAGE 34

19 ingthestandardvariabletransformation f = f ! " 1 2 S ! [49],Eq.(2.16)canbee!ectivelyrewritten asane!ectivelytimeexplicitmethodas f ! ( " x + " e ! $ t ,t + $ t )= f ! ( " x,t )+ # ! ( " x,t )+ S ! ( " x,t ) , (2.18) wherewehaveconsideredtheusuallatticeunitsbytakingthetimestep $ t equaltounity.Next, toobtainthestructureofthechangeofdi!erentmomentsundercascadedcollision 0 g andthe sourceterm S ,wedeÞnethefollowingdiscretecentralmoments: 6 7 7 7 7 7 7 7 7 8 ö * x m y n ö * eq x m y n ö & x m y n ö ø * x m y n 9 : : : : : : : : ; = " ! 6 7 7 7 7 7 7 7 7 8 f ! f eq ! S ! ø f ! 9 : : : : : : : : ; ( e ! x " u x ) m ( e ! y " u y ) n . (2.19) Bymatchingthediscreteandcontinuouscentralmomentsofequilibriaandsourcesatdi!erent orders,i.e., 0 * eq x m y n = 0 $ M x m y n and 0 & x m y n = 0 % F x m y n ,weget | 0 * eq x m y n $ =( 0 * eq 0 , 0 * eq x , 0 * eq y , 0 * eq xx , 0 * eq yy , 0 * eq xy , 0 * eq xxy , 0 * eq xyy , 0 * eq xxyy )   =( # , 0 , 0 ,c 2 s # ,c 2 s # , 0 , 0 , 0 ,c 4 s # )   , (2.20) | 0 & x m y n $ =( 0 & 0 , 0 & x , 0 & y , 0 & xx , 0 & yy , 0 & xy , 0 & xxy , 0 & xyy , 0 & xxyy )   =(0 ,F x ,F y , 0 , 0 , 0 , 0 , 0 , 0)   , (2.21) Inpractice,thecomputationsinthecascadedLBEareperformedbymeansofwhichdiscreteraw moments,whicharedeÞnedas 6 7 7 7 7 7 7 7 7 8 ö * ! x m y n ö * eq ! x m y n ö & ! x m y n ö ø * ! x m y n 9 : : : : : : : : ; = " ! 6 7 7 7 7 7 7 7 7 8 f ! f eq ! S ! ø f ! 9 : : : : : : : : ; e m ! x e n ! y (2.22) Thevariousnon-conservedmoments 0 * ! x m y n ,whichareneededintheformulationofthecascaded

PAGE 35

20 collisiontermmaybewrittenas 0 * ! xx = " ! f ! e 2 ! x = f 1 + f 3 + f 5 + f 6 + f 7 + f 8 , (2.23a) 0 * ! yy = " ! f ! e 2 ! y = f 2 + f 4 + f 5 + f 6 + f 7 + f 8 , (2.23b) 0 * ! xy = " ! f ! e ! x e ! y = f 5 " f 6 + f 7 " f 8 , (2.23c) 0 * ! xxy = " ! f ! e 2 ! x e ! y = f 5 + f 6 " f 7 " f 8 , (2.23d) 0 * ! xyy = " ! f ! e ! x e 2 ! y = f 5 " f 6 " f 7 + f 8 , (2.23e) 0 * ! xxyy = " ! f ! e 2 ! x e 2 ! y = f 5 + f 6 + f 7 + f 8 . (2.23f) where 0 * x m y n = 0 * x m y n " 1 2 0 & x m y n .Byapplyingthebinomialtheoremtothecentralmomentsof thesourcetermatdi!erentorders,thecorrespondingrawmomentscanbederived,andcanbe writtenas 0 & ! 0 =0 , 0 & ! x = F x , 0 & ! y = F y , 0 & ! xx =2 F x u x , 0 & ! yy =2 F y u y , 0 & ! xy = F x u y + F y u x , 0 & ! xxy = F y u 2 x +2 F x u x u y , 0 & ! xyy = F x u 2 y +2 F y u y u x , 0 & ! xxyy =2 F x u x u 2 y +2 F y u y u 2 x . Inordertoobtaintheexpressionsforthesourcetermsinthevelocityspace S ! ,weÞrstevaluate thesourcemomentsprojectedtotheorthogonalmomentspace,i.e. 0 m s ( = % K ( | S ! $ ,where + = 0 , 1 , 2 ,..., 8: 0 m s 0 =0 , 0 m s 1 = F x , 0 m s 2 = F y , 0 m s 3 =6( F x u x + F y u y ) , 0 m s 4 =2( F x u x " F y u y ) , 0 m s 5 =( F x u y + F y u x ) , 0 m s 6 =(2 " 3 u 2 x ) F y " 6 F x u x u y , 0 m s 7 =(2 " 3 u 2 y ) F x " 6 F y u y u x , 0 m s 8 =6 ( (3 u 2 y " 2) F x u x +(3 u 2 x " 2) F y u y ) . Then,using( Ká S ) ! =( 0 m s 0 , 0 m s 1 , 0 m s 2 ,..., 0 m s 8 )   & | 0 m s ! $ )andinvertingitbyexploitingtheorthogonalpropertyof K ,wecanexplicitlywritetheÞnalexpressionsforthesourcetermsinthe

PAGE 36

21 velocityspace S as S 0 = 1 9 < " m s 3 + m s 8 = , (2.24a) S 1 = 1 36 < 6 m s 1 " m s 3 +9 m s 4 +6 m s 7 " 2 m s 8 = , (2.24b) S 2 = 1 36 < 6 m s 2 " m s 3 " 9 m s 4 +6 m s 6 " 2 m s 8 = , (2.24c) S 3 = 1 36 < " 6 m s 1 " m s 3 +9 m s 4 " 6 m s 7 " 2 m s 8 = , (2.24d) S 4 = 1 36 < " 6 m s 2 " m s 3 " 9 m s 4 " 6 m s 6 " 2 m s 8 = , (2.24e) S 5 = 1 36 < 6 m s 1 +6 m s 2 +2 m s 3 +9 m s 5 " 3 m s 6 " 3 m s 7 + m s 8 = , (2.24f) S 6 = 1 36 < " 6 m s 1 +6 m s 2 +2 m s 3 " 9 m s 5 " 3 m s 6 +3 m s 7 + m s 8 = , (2.24g) S 7 = 1 36 < " 6 m s 1 " 6 m s 2 +2 m s 3 +9 m s 5 +3 m s 6 +3 m s 7 + m s 8 = , (2.24h) S 8 = 1 36 < 6 m s 1 " 6 m s 2 +2 m s 3 " 9 m s 5 +3 m s 6 " 3 m s 7 + m s 8 = . (2.24i) Inordertoconstructthecascadedoperator,wealsoneedtherawmomentsofthechangeofmomentsundercollision 0 g atdi!erentorders,i.e., > ! ( Ká 0 g ) ! e m ! x e n ! y = > ( % K ( | e m ! x e n ! y $ 0 g ( .Since themass(zerothordermoment)andmomentum(Þrstordermoment)arecollisioninvariants,we have 0 g 0 = 0 g 1 = 0 g 2 =0.Then,byexploitingtheorthogonalpropertyofthebasisvectors K ( ,it followsthat > ! ( Ká 0 g ) ! e 2 ! x =6 0 g 3 +2 0 g 4 , > ! ( Ká 0 g ) ! e 2 ! y =6 0 g 3 " 2 0 g 4 , > ! ( Ká 0 g ) ! e ! x e ! y =4 0 g 5 , > ! ( Ká 0 g ) ! e 2 ! x e ! y =4 0 g 2 " 4 0 g 6 , > ! ( Ká 0 g ) ! e ! x e 2 ! y =4 0 g 1 " 4 0 g 7 , > ! ( Ká 0 g ) ! e 2 ! x e 2 ! y =8 0 g 3 +4 0 g 8 . Basedontheabove,theoverallproceduretoobtainthechangeofmomentsundercollision 0 g ( as follows:Startingfromthelowestorder,thenon-conserved,post-collisioncentralmoments(i.e., thesecondordermomentcomponents),wetentativelyequatethemtothecorrespondingcentral momentsequilibria.Oncetheexpressionsforthechangeofmoments 0 g ( ( + ' 3)areobtained, wethendiscardtheequilibriumassumptionandallowforarelaxationofcentralmomentsduringcollision.Thisisachievedbymultiplyingthoseresultingtermsin 0 g ( thatarenotyetinthe

PAGE 37

22 post-collisionstatesbyacorrespondingrelaxationparameter , ( .Thedetailsofsuchaprocedure arepresentedin[50,22].Inaddition,duetothediscretenessofthe D 2 Q 9lattice,itintroduces aliasinge!ectsordegeneraciesintherepresentationofthethirdorderdiagonalmoments,i.e., 0 * ! xxx = 0 * ! x and 0 * ! yyy = 0 * ! y .Thisleadstonon-Galileaninvariant(GI)cubicvelocityerrors,resultinginananisotropicstresstensor,whichcouldhavee!ectonthequalityofthenumerical solutionthatcanbeobtained,includingnegativeimplicationsonthestabilityofthenumerical scheme,regardlessofthecollisionmodelusedintheLBmethod.Inordertoeliminatesuchcubicvelocityerrors,correctionseithertotherelaxationtimes[51]ortotheequilibria[52],[53]of thesecondordermomentsneedtobeconsidered.Inthiswork,fortheÞrsttime,weextendthis approachofintroducingtheGIcorrectionstothesecondorderequilibriummomentsforthesolutionofnon-NewtonianßowsinacascadedLBformulation.TheÞnalexpressionsforthechange ofmomentsundercollision 0 g ( canthenbesummarizedasfollows: 0 g 0 = 0 g 1 = 0 g 2 =0 , (2.25a) 0 g 3 = , 3 12 ? 2 3 # + # ( u 2 x + u 2 y ) " ( 0 * ! xx + 0 * ! yy ) " 1 2 ( 0 & ! xx + 0 & ! yy )+ < 3 x % x u x + 3 y % y u y = @ , (2.25b) 0 g 4 = , 4 4 ? # ( u 2 x " u 2 y ) " ( 0 * ! xx " 0 * ! yy ) " 1 2 ( 0 & ! xx " 0 & ! yy )+ < 4 x % x u x + 4 y % y u y = @ , (2.25c) 0 g 5 = , 5 4 ? # u x u y " 0 * ! xy " 1 2 0 & ! xy @ , (2.25d) 0 g 6 = , 6 4 ? 2 # u 2 x u y + 0 * ! xxy " 2 u x 0 * ! xy " u y 0 * ! xx " 1 2 0 & xxy @ " 1 2 u y (3 0 g 3 + 0 g 4 ) " 2 u x 0 g 5 , (2.25e) 0 g 7 = , 7 4 ? 2 # u x u 2 y + 0 * ! xyy " 2 u y 0 * ! xy " u x 0 * ! yy " 1 2 0 & xyy @ " 1 2 u x (3 0 g 3 " 0 g 4 ) " 2 u y 0 g 5 , (2.25f) 0 g 8 = , 8 4 ? 1 9 # +3 # u 2 x u 2 y " 3 0 * ! xxyy " 2 u x 0 * ! xyy " 2 u y 0 * ! xxy + u 2 x 0 * ! yy + u 2 y 0 * ! xx +4 u x u y 0 * ! xy 4 " 1 2 0 & ! xxyy @ " 2 0 g 3 " 1 2 u 2 y (3 0 g 3 + 0 g 4 ) " 1 2 u 2 x (3 0 g 3 " 0 g 4 ) " 4 u x u y 0 g 5 " 2 u y 0 g 6 " 2 u x 0 g 7 . (2.25g) Here,thepresenceofthediagonalvelocitygradientterms,i.e., % x u x and % y u y andtheirasso-

PAGE 38

23 ciatedcoe"cients 3 x , 3 y , 4 x and 4 y inthechangeofmomentsundercollision 0 g 3 and 0 g 4 forthe relaxationofthesecondordercentralmomentscorrespondtothecorrectionsintroducedtoeliminatenon-GIcubicvelocityerrors.Thesecoe"cientsarethenwrittenasfollows 3 x = " 3 $ 1 , 3 " 1 2 % # u 2 x , 3 y = " 3 $ 1 , 3 " 1 2 % # u 2 y , (2.26a) 4 x = " 3 $ 1 , 4 " 1 2 % # u 2 x , 4 y = " 3 $ 1 , 4 " 1 2 % # u 2 y . (2.26b) AChapman-Enskoganalysisofthecentralmomentsbasedcascadedthatderivesthenon-GI cubicvelocitytermsandtheireliminationviacorrectionstermsinthesecondordermoment equilibriaarepresentedinAppendixA.Asaresult,theconsistencyofthecascadedLBEtothe Navier-Stokesequationsisshown,wheretherelaxationtimes , 3 , , 4 and , 5 arerelatedtothe bulkandshearkinematicviscosities,respectively,as ( = c 2 s A 1 ) 3 " 1 2 B , . = c 2 s A 1 ) j " 1 2 B j =4 , 5 . Then,forthesimulationofnon-NewtonianßowsusingthecascadedLBmethodbasedonthe constitutiverelationsinvolvingtheshearrate | ú ' ij | dependentapparentviscosities µ ( | ú ' ij | )= #. ( | ú ' ij | )ofeitherthepowerlaworBinghamßuids(Sec.3.2),therelaxationtimes , 4 and , 5 can bemadetoadjustlocallyasfollows: , 4 = , 4 ( " x,t )= , 5 = , 5 ( " x,t )= & µ ( | ú ' ij | ) # c 2 s + 1 2 ' ! 1 , (2.27) where µ ( | ú ' ij | )canbeobtainedfromeitherEq.(2.6)or(2.9).Themagnitudeoftheshearrate | ú ' ij | isrelatedtothesecondinvariantofthestrainratetensor S ij via | ú ' ij | = # 2 S ij S ij = C 2( S 2 xx +2 S 2 xy + S 2 yy ) . (2.28) Inthepresentwork,thecomponentsofthestrainratetensor S xx , S yy and S yy areobtainedlocallybasedonthesecondordernon-equilibriummoments,ratherthanusingÞnitedi!erenceapproximations.AderivationoftheseexpressionsispresentedinAppendixA,whichtakesintoaccountthee!ectofthecorrectionsfortheeliminationofthenon-GIterms.TheÞnalexpressions forthestrainratetensoraresummarizedhereasfollows: S xx = % x u x = [÷ c 2 ö m neq 3 + c 2 ö m neq 4 ] [ " c 1 ÷ c 2 " ÷ c 1 c 2 ] , (2.29a)

PAGE 39

24 S yy = % y u y = [ " c 1 ö m neq 4 +÷ c 1 ö m neq 3 ] [ " c 1 ÷ c 2 " ÷ c 1 c 2 ] , (2.29b) S xy = 1 2 ( % x u y + % y u x )= " 1 2 c 3 ö m neq 5 , (2.29c) where ö m neq 3 = 0 * neq ! xx + 0 * neq ! yy =( 0 * ! xx + 0 * ! yy ) " $ 2 3 # + # ( u 2 x + u 2 y ) % , (2.30a) ö m neq 4 = 0 * neq ! xx + 0 * neq ! yy =( 0 * ! xx + 0 * ! yy ) " # ( u 2 x " u 2 y ) , (2.30b) ö m neq 5 = 0 * neq ! xx = 0 * ! xx " # u x u y . (2.30c) and c 1 = 2 # 3 , 3 & 1 " 9 4 u 2 x , 3 ' ,c 2 = 2 # 3 , 3 & 1 " 9 4 u 2 y , 3 ' (2.31a) ÷ c 1 = 2 # 3 , 4 & 1 " 9 4 u 2 x , 4 ' , ÷ c 2 = 2 # 3 , 4 & 1 " 9 4 u 2 y , 4 ' (2.31b) c 3 = # 3 , 5 . (2.31c) TheremainingrelaxationparametersinEqs.(2.25a)Ð(2.25g)i.e., , 6 , , 7 and , 8 inßuencenumericalstability.Ananalysisofthee!ectoftheseparametersispresentedin[54].Inthispresent work,weset , 3 = , 6 = , 7 = , 8 =1forsimplicity,andlocallyvary , 4 and , 5 usingEqs.(2.27)Ð (2.31c)tosimulateincompressiblenon-Newtonianßuidßows. Then,Eq.( ?? )forthecascadedLBEcanberecastasthefollowingcollisionandstreamingsteps: ! f ! ( " x,t )= f ! ( " x,t )+( Ká 0 g ) ! ( " x,t )+ S ! ( " x,t ) , (2.32a) f ! ( " x + " e ! ,t + $ t )= ! f ! ( " x,t ) , (2.32b) Expandingthecollisionterm( Ká 0 g ) ! ,thecomponentsofthepost-collisiondistributionfunction

PAGE 40

25 ! f ! canbeexplicitlyexpressedas ! f 0 = f 0 +[ 0 g 0 " 4( 0 g 3 " 0 g 8 )]+ S 0 , (2.33a) ! f 1 = f 1 +[ 0 g 0 + 0 g 1 " 0 g 3 + 0 g 4 +2( 0 g 7 " 0 g 8 )]+ S 1 , (2.33b) ! f 2 = f 2 +[ 0 g 0 + 0 g 2 " 0 g 3 " 0 g 4 +2( 0 g 6 " 0 g 8 )]+ S 2 ,, (2.33c) ! f 3 = f 3 +[ 0 g 0 " 0 g 1 " 0 g 3 + 0 g 4 " 2( 0 g 7 + 0 g 8 )]+ S 3 ,, (2.33d) ! f 4 = f 4 +[ 0 g 0 " 0 g 2 " 0 g 3 " 0 g 4 " 2( 0 g 6 + 0 g 8 )]+ S 4 , (2.33e) ! f 5 = f 5 +[ 0 g 0 + 0 g 1 + 0 g 2 +2 0 g 3 + 0 g 5 " 0 g 6 " 0 g 7 + 0 g 8 ]+ S 5 , (2.33f) ! f 6 = f 6 +[ 0 g 0 " 0 g 1 + 0 g 2 +2 0 g 3 " 0 g 5 " 0 g 6 + 0 g 7 + 0 g 8 ]+ S 6 , (2.33g) ! f 7 = f 7 +[ 0 g 0 " 0 g 1 " 0 g 2 +2 0 g 3 + 0 g 5 + 0 g 6 + 0 g 7 + 0 g 8 ]+ S 7 , (2.33h) ! f 8 = f 8 +[ 0 g 0 + 0 g 1 " 0 g 2 +2 0 g 3 " 0 g 5 + 0 g 6 " 0 g 7 + 0 g 8 ]+ S 8 . (2.33i) Theßuidßowvariables,suchasthedensity # ,velocityÞeld " u andpressure P arethenobtained as # = 8 " ! =0 f ! ,# u i = 8 " ! =0 f ! e ! i + 1 2 F i ,p = c 2 s # (2.34) 2.4 ResultsandDiscussion WewillnowpresentanumericalvalidationofthecascadedLBmethodbasedoncentralmomentsforavarietyofnon-Newtonian2Dbenchmarkßowproblems.Theseinclude(a)thepressuredrivenßowbetweentwoparallelplates,(b)lid-drivencavityßow,(c)crossßowoveracircular cylinder,and(d)theßowoverapairofcircularcylindersintandemconÞguration(seeFig.2.1). 2.4.1 Channelßowofpower-lawandBinghamßuids First,afully-developedpressuredrivenßowofpowerlawßuidsbetweentwoparallelplateswith aspacingheight H isconsidered.TheanalyticalsolutionforthevelocityproÞle u ( y )forthis

PAGE 41

26 Y X (a) Y X U (b) 10 D 30 D 20 D (c) Lu 10D G Ld U Y X (d) FIGURE2.1:Schematicsofthenon-Newtonianbenchmarkßowproblems.(a)Pressure-driven ßowbetweentwoparallelplates,(b)lid-drivencavityßow,(c)cross-ßowoveracircularcylinder, and(d)ßowoverapairofcircularcylindersintandemarrangement.

PAGE 42

27 problemcanbewrittenas u y = u max D 1 " $ | y | H/ 2 % ( n +1) /n E , (2.35) wherethemaximumvelocity u max atthecenterofthechannelisgivenby u max = n n +1 $ " 1 µ p % P % x % 1 /n $ H 2 % ( n +1) /n . (2.36) Here,thepressuregradient " %p/%x isrepresentedbyaconstantbodyforceinthex-direction, i.e.( F x ,F y )=( " %p/%x, 0)inthesimulations.Aperiodicboundaryconditionisappliedalong thestreamwiseßowdirection,andtheno-slipboundaryconditionsatthechannelwallsarerepresentedusingthestandardhalf-waybouncebackscheme.TheReynoldsnumberforthisßowmay bewrittenas Re = # H n u 2 ! n max /µ p .Weconsideredagridresolutionof3 ! 101forresolvingthe domainofthepowerlawßuidswith Re =100andthevaluesofthepowerlawindex n equalto 0 . 5 , 1 . 0and2 . 0,therebyencompassingtheshearthinning,Newtonianandshearthickeningßuids.Figure(2.2)presentsacomparisonoftheconvergedsteadystatenormalizedvelocityproÞles computedusingthecascadedLBmethod(solidlines)withanalyticalsolution(symbols)fordifferentvaluesof n .anexcellentagreementbetweenthesetwosolutionsareseenforthesimulation ofpowerlawßuids.Next,weconsiderthesimulationofaBinghamßuidßowbetweentwoparallelplatesofheight H .TheanalyticalsolutionforthefullydevelopedvelocityproÞleforthisßow canberepresentedas u ( y )= F G G G H G G G I 1 2 µ B < " * P * x = 3 < H 2 = 2 " y 2 + 4 " + B µ B < H 2 " y + = , 0 ( | y | ( y + , 1 2 µ B < " * P * x = 3 < H 2 = 2 " y 2 4 " + B µ B < H 2 " | y | = ,y + < | y | ( H 2 , (2.37) Here, y + = ! B / ( " %p/%x )isthelocation,wheretheßuidbeginstoyieldarisingfromabalance betweentheyieldforceandappliedpressureforce.Itfollowsthatthemaximumßowvelocity u max canbewrittenas u max = 1 2 µ B $ " % P % x % D $ H 2 % 2 " y 2 + E " ! B µ B $ H 2 " y + % (2.38)

PAGE 43

28 ! 0.5 0 0.5 0 0.2 0.4 0.6 0.8 1 1.2 y/L u/U p n=0.5, Analytical n=0.5, Numerical n=1.0, Analytical n=1.0, Numerical n=1.5,Analytical n=1.5,Numerical FIGURE2.2:ComparisonofthenormalizedvelocityproÞlesofpowerlawßuidsinachannel computedusingthecascadedLBmethodandtheanalyticalsolutionfor n =0 . 5 , 1 . 0and1 . 5 at Re =100. TheReynoldsnumberforthisßowconÞguration,then,becomes Re = # u max H/µ B .Asbefore, wetreatthepresenceofapressuregradientasabodyforceinaperiodicdomainofsize3 ! 101 inourcomputations.Figure2.3showsthenormalizedvelocityproÞlesforthefollowingvaluesof theyieldstress ! B =0,4 ! 10 ! 4 and8 ! 10 ! 4 at Re =100.Itcanbeseenthatthenumerical solutionobtainedusingthecascadedLBmethodisinverygoodagreementwiththeanalytical solutionforBinghamplasticßowthroughachannel. InordertoassesstherateofspatialgridconvergenceofthecascadedLBmethodforthesimulationofnon-Newtonianßuid,wenowstudythevariationintherelativeglobalerrorinthevelocityÞeld E r measuredunderadiscrete L 2 normgivenasfollows: E r = J K K L > < u c ( i,j ) " u a ( i,j ) = 2 > u 2 a ( i,j ) , (2.39) where u c and u a representthecomputedandanalyticalsolutions,respectively,andthesummationiscarriedoverthewholeßowdomain.Inordertostudygridconvergenceofthenumerical solutiontotheincompressiblenon-Newtonianßow,weadoptthedi!usivescalingasthegridres-

PAGE 44

29 ! 0.5 0 0.5 0 0.2 0.4 0.6 0.8 1 1.2 y/L u/U p ! b =0, Analytical ! b =0, Numerical ! b =0.00004, Analytical ! b =0.00004, Numerical ! b =0.00008,Analytical ! b =0.00008,Numerical FIGURE2.3:ComparisonofthenormalizedvelocityproÞlesofBinghamßuidsinachannelcomputedusingthecascadedLBmethodandtheanalyticalsolutionfor ! B =0 , 4 ! 10 ! 4 and8 ! 10 ! 4 atRe=100. olution $ x =1 /N ,where N isthenumberofgridnodes,isvaried.Thatis,themaximumßow velocityreducesproportionaltothegridspacingasthegridresolutionisincreased.Figure2.4 presentstherelationglobalerrorasafunctionofthegridspacinginthelog-logscale,forthe simulationoftheßowofpowerlawßuidsinachannelfor n =0 . 5and1 . 0.Itcanbeseenthat theslopeoftheresultingcurvesare " 2 . 0indicatingthecascadedLBmethodissecondorderaccurateinspaceforthecomputationofpowerlawßuidßows. 2.4.2 2Dlid-drivencavityßowofnon-Newtonianßuids Lid-drivencavityßowisanimportantcomplexßowbenchmarkproblemwithavarietyofßow regimescharacterizedbythepresenceofsingularities,vortexstructuresofdi!erentscalesand shapes,andrecirculatingpatternsthatevolveintoßowinstabilitiesastheReynoldsnumberis increased.Inthecaseofnon-Newtonianßuids,additionalcharacteristicparameters,suchasthe powerlawindexortheBinghamnumber(seebelow)representingtheyieldingbehaviorsigniÞ-

PAGE 45

30 ! 4 ! 3.8 ! 3.6 ! 3.4 ! 3.2 ! 3 ! 2.8 ! 2.6 ! 2.4 ! 8 ! 7.5 ! 7 ! 6.5 ! 6 ! 5.5 ! 5 ! 4.5 ! 4 ! 3.5 Ln (Er) Ln ( ! x) n=0.5 n=1 slope= ! 2 FIGURE2.4:Relativeglobalerroragainstthevariationinthegridspacing,plottedinlog-log scales,forthesimulationofpowerlawßuidßowinachannelat n =0 . 5and1 . 0obtainedusing thecascadedLBmethod. cantlycontributetothedynamicsofsuchsheardrivenßows.Somelimitednumericalstudiesfor thisnon-NewtonianßowconÞgurationareasfollows.Reference[55]performedasimulationof lid-drivencavityßowsofpowerlawßuidsusingap-versionoftheÞniteelementmethod.Benchmarknumericalsolutionsofliddrivencavityßowsofavarietyofßuids,includingpowerlawand Binghamplasticßuids,obtainedusingathirdorderupwardÞnitevolumemethodwerereported in[2].SomeotherversionsofÞniteelementsbasednumericalcomputationswereperformedand studiedforthisßowproblembyRefs.[56,57].Inaddition,morerecently,Ref.[45]presented simulationsofsuchnon-Newtonianlid-drivencavityßowsusingavariantoftheLBmethod.The geometryoftheßowestablishedbythemotionofthelidatavelocity U p inasquarecavityof sidelength L isshowninFig.2.1b.Asbefore,theno-slipboundaryconditionsarerepresented usingthestandardhalf-waybouncebackscheme,wherethewallvelocityforthetopboundary isincorporatedbyamomentumcorrectiontothebouncebackscheme.Inthecaseofpowerlaw ßuids,thecharacteristicReynoldsnumberisdeÞnedas Re = # L n U 2 ! n p /µ p ,whilefortheBinghamßuids,theReynoldsnumberandtheBinghamnumberarerepresentedby Re = # U p L/µ B

PAGE 46

31 and Bn = & B L/ ( µ B U p ),respectively.Figure2.5presentsthecomputedsteadystatevelocity components u and v alongtheverticalandhorizontalcenterlines,respectively,forpowerlawßuidsat Re =100anddi!erentvaluesofthepowerlawindex n .Agridresolutionof201 ! 201was employedandshearthinning( n =0 . 5),Newtonian( n =1)andshearthickening( n =1 . 5)ßuids wereconsidered.Also,showninthisÞgurearethebenchmarknumericalresultsofNeofytou[2]. ItcanbeseenthatthecascadedLBmethodisinverygoodagreementwiththebenchmarkdata forallthecasesofthepowerlawßuidsconsidered.Also,itseenthatasthepowerlawindex n increases,themagnitudeofthepeakverticalvelocityalongthehorizontalincreases,anditslocationisshiftedclosetotheverticalwalltotheleftsideofthecavity.Furthermore,thehorizontal velocitycomponentisfoundtohavedramaticallydi!erentbehaviorfortheshearthickeningßuid, whencomparedtotheshearthinningßuid,withthelatterexhibitingsharpervariationsacross theverticalcenterlineofthecavity.Alltheseobservationsareinqualitativeandquantitative agreementwiththebenchmarkresults[2]. Furthermore,weperformedsimulationoftheßowofBinghamßuidsinaliddrivencavityusinga gridresolutionof401 ! 401at Re =100and Bn =0 . 1and1 . 0,wheretheBinghamnumber Bn representstherelativeroleoftheyieldstressfortheviscoplasticßuidsconsidered.Thischooseof gridresolutionleadstogridindependentresults.ThecomputedvelocityproÞlesalongthecenterlinesofthecavityusingthecascadedLBschemeshowninFig.(2.6)areagainingoodagreement withthebenchmarknumericalsolutionsofNeofytou[2].Astheyieldingbehaviorvariesasthe Bn increasesfrom0 . 1to1 . 0,thehorizontalvelocitycomponentexhibitsaninßexioninitsproÞle,whilethepeakvelocityoftheverticalcomponentdecreases.Thesefeatures,whichwereÞrst reportedin[2],iswellreproducedbythepresentcascadedLBscheme. TheinßuenceoftheReynoldsnumberonthestreamlinepattemforÞxed n =0 . 5(shearthickening)byconsidering Re =100 , 400and1000isshowninFig.(2.7).Itcanbeseenthatthe centeroftheprimaryvortexgetsshiftedtowardsthetoprightcorneras Re isdecreased.Furthermore,thesizeofthesecondvorticalstructuresatthebottomcornersbecomelargerwithan

PAGE 47

32 FIGURE2.5:Computedvelocitycomponent u/U p alongtheverticalcenterlineandcomponent v/U p alongthehorizontalcenterlineobtainedusingthecascadedLBmethod(lines)andcomparedwiththebenchmarksolutionofNeofytou[2](symbols)forlid-drivencavityßowofpower lawßuids.Reynoldsnumber Re =100andthepowerlawindex n =0 . 5 , 1 . 0 , 1 . 5.

PAGE 48

33 (a)Re=100 (b)Re=1000 FIGURE2.6:Computedvelocitycomponent u/U p alongtheverticalcenterlineandcomponent v/U p alongthehorizontalcenterlineobtainedusingthecascadedLBmethod(lines)andcomparedwiththebenchmarksolutionsofNeofytou[2](symbols)forlid-drivencavityßowofBinghamßuids.Reynoldsnumber Re =100andBinghamnumber Bn =0 . 1 , 1 . 0 .

PAGE 49

34 increasein Re .Thesetrendsareconsistentwiththeclassicallid-drivencavityßowresultspresentedin[58].Furthermore,Figs.2.8and2.9presentthestreamlinecontoursat Re =100 , 400 and1000for n =1(Newtonian)and n =1 . 5(shearthickening)ßuids,respectively.Theoverall ßowpatternsforvariousReynoldsnumbersandpowerlawindicesareinqualitativeagreement withpriorbenchmarkstudies[2,58].Thee!ectofthepowerlawindex n forÞxedRe=400 onthestreamlinecontoursforthepowerlawßuidßowinacavityisstudiedinFig.(2.10).The centeroftheprimaryvortexisseentomovetowardstothetoplidat n isincreased.Inaddition,theshearthickeningßuidexhibitsrelativelylargersecondaryvorticesatthebottomleftand rightcornerwhencomparedtotheNewtonianortheshearthinningßuid. TheinßuenceoftheyieldingbehaviorortheBinghamnumber Bn ( Bn =0,0.01and0.1)for Þxed Re =400onthestreamlinepatternsfortheidealviscoplastic(Bingham)ßuidisshown inFig.2.11.Ingeneral,as Bn isincreased,theprimaryvortexisseentomovefromthecentralzonetowardsthetoprightcornerofthelid.Also,therelativemagnitudeoftheyieldstress, i.e., Bn isfoundtohaveamarkede!ectonthesecondaryrecirculatingeddiesinthecavity.In particular,as Bn isprogressivelyincreased,thesizeofthesecondaryverticalstructuresisobservedtobecomesmaller.Thesevariationsoftheßowpatternswith Bn isinagreementwiththe benchmarksolutions[2].Aquantitativecomparisonofthelocationofthecenteroftheprimary vortexismadefor n =1atdi!erentReynoldsnumbers( Re =100 , 400 , 100)usingagridresolutionof401 ! 401inFig.2.13.FromthisÞgure,itisevidentthatastheReincreased,theposition ofthecenterofthemainvortexmovestowardsthegeometriccenterofthecavity.Thecomputed valuesofthelocationsoftheprimaryvortexatdi!erentReareinverygoodagreementwiththe benchmarkdata[2]. Inaddition,Fig.2.14apresentsthelocationoftheprimaryvortexforpowerlawßuidsataÞxed Re =400andfordi!erentvalueofthepowerlawindex n ( n =0.75,1.0,1.25),whileFig.(2.14b) showsitfortheBinghamßuidsatthesame-Þxed Re andfordi!erentvaluesoftheBingham number Bn ( Bn =0 , 0 . 1 , 1 . 0)usingagridresolutionof401 ! 401.Itisevidentifeitherthe

PAGE 50

35 X Y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (a)Re=100 X Y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (b)Re=400 X Y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (c)Re=1000 FIGURE2.7:StreamlinecontoursforliddrivencavityßowatÞxed n =0 . 5and(a)Re=100, (b)Re=400,(c)Re=1000computedusingthecascadedLBmethod.

PAGE 51

36 X Y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (a)Re=100 X Y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (b)Re=400 X Y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (c)Re=1000 FIGURE2.8:StreamlinecontoursforliddrivencavityßowatÞxed n =1 . 0and(a)Re=100, (b)Re=400,(c)Re=1000computedusingthecascadedLBmethod.

PAGE 52

37 X Y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (a)Re=100 X Y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (b)Re=400 X Y 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (c)Re=1000 FIGURE2.9:StreamlinecontoursforliddrivencavityßowatÞxed n =1 . 5and(a)Re=100, (b)Re=400,(c)Re=1000computedusingthecascadedLBmethod. Re400n7_5-eps-converted-to.pdf (a) n = 0 . 75 Re400n1-eps-converted-to.pdf (b) n = 1 . 0 Re400n1_25-eps-converted-to.pdf (c) n = 1 . 25 FIGURE2.10:StreamlinecontoursforliddrivencavityßowofpowerlawßuidsatÞxedRe= 400(a) n =0 . 75,(b) n =1 . 0,(c) n =1 . 25computedusingthecascadedLBmethod.

PAGE 53

38 figs/Re400bn0-eps-converted-to.pdf (a) Bn = 0 figs/Re400bn0_01-eps-converted-to.pdf (b) Bn = 0 . 01 figs/Re400bn0_1-eps-converted-to.pdf (c) Bn = 0 . 1 FIGURE2.11:StreamlinecontoursforliddrivencavityßowofBinghamßuidsatÞxed Re =400 and(a) Bn =0,(b) Bn =0 . 01,(c) Bn =0 . 1computedusingthecascadedLBmethod. 0.52 0.54 0.56 0.58 0.6 0.62 0.55 0.6 0.65 0.7 0.75 X/L Y/L Re=100( Neofytou (2005)) Cascaded MRT Re=400( Neofytou (2005)) Cascaded MRT Re=1000( Neofytou (2005)) Cascaded MRT FIGURE2.12:Locationoftheprimaryvortexforlid-drivencavityßowat n =1.Re=100 , 400 and1000computedusingthecascadedLBmethod.

PAGE 54

39 powerindex n increases,i.e.,iftheßuidisshearthickeningratherthanshearthinningorthe Binghamnumber Bn increases,thelocationoftheprimaryvortexprogressivelymovestoward thedirectionofthetoprightcornerofthecavity,whichareconsistentwiththeÞndingsofthe Neofytou[2]. 0.52 0.54 0.56 0.58 0.6 0.62 0.55 0.6 0.65 0.7 0.75 X/L Y/L Re=100( Neofytou (2005)) Cascaded MRT Re=400( Neofytou (2005)) Cascaded MRT Re=1000( Neofytou (2005)) Cascaded MRT FIGURE2.13:Locationoftheprimaryvortexforlid-drivencavityßowat n =1,Re=100 , 400 and1000computedusingthecascadedLBmethod. 2.4.3 Non-Newtonianßowofpowerlawßuidsoverasinglecircularcylinder Anothercomplexnon-Newtonianßuidmechanicsbenchmarkproblemisthecrossßowoveracircularcylinder.Weconsiderthe2Dßowofincompressiblepowerlawßuidwithauniformvelocity U " overanunconÞnedcircularcylinderofdiameter D .Thecircularcylinderiskeptanupstream distanceof L u =10 D fromtheinlettoitscenterandtheoutletislocatedatadownstreamdistanceof L d =20 D fromthecenterofthecylinder.Thecylinderissymmetricallyplacedalong theverticaldirectionofthecomputationaldomain,whichisspanningdistanceof20 D .Sucha conÞgurationwasnumericallystudiedusingthesolutionoftheNavier-Stokesequationswitha nonlinearconstitutiverelationforthepowerlawßuidsusingaÞnitevolumemethodinSoareset

PAGE 55

40 0.55 0.551 0.552 0.553 0.554 0.555 0.556 0.557 0.57 0.58 0.59 0.6 0.61 0.62 0.63 X/L Y/L n=1.25 n=1.0 n=0.75 (a) 0.59 0.6 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.66 0.68 0.7 0.72 0.74 0.76 X/L Y/L Bn=0 Bn=0.1 Bn=1 (b) FIGURE2.14:locationoftheprimaryvortexforlid-drivencavityßowatRe=400for(a) powerlawßuidsfor n =0 . 75 , 1 . 0 , 1 . 25 , (b)Binghamßuidsfor Bn =0 , 0 . 1 , 1 . 0computedusing thecascadedLBmethod.

PAGE 56

41 al[59]andBhartietal[3].Inthiswork,weusethebenchmarkdataofBhartietal[3]forcomparison. ForthesimulationusingthecascadedLBmethod,inthiswork,weimposeuniformßowvelocity U " attheinlet,convectiveboundaryconditionsattheoutletandfree-slipboundaryconditionsareappliedatthetwosideboundariesofthecomputationaldomain.Theno-slipboundaryconditiononthecurvedboundariesofthecircularcylinderisimposedbyapplyingtheinterpolatedbouncebackscheme[60],whereasecondorderinterpolationisusedinthisstudy. Thisrheologicalßowproblemischaracterizedbythepowerlawindex n andtheReynoldsnumber Re givenby Re = # D n U 2 ! n " /µ p .Weconsideragridresolutionof2001 ! 501.Figure2.15 presentsthestreamlinecontoursforthesteadyßowofashearthinningßuid( n =0 . 6)atvarious Re of5 , 10 , 20,and30.Itcanbeseenthatatalow Re of5,thereisasymmetricßowaround thecircularcylinderwithoutanynoticeablerecirculationzonebehindit.However,astheReincreases,thesizeoftherecirculatingzoneanditslengthisfoundtoprogressivelyincrease,with thepointofseparationmovingforwardonthesurfaceofthecylinder.Theseobservationsare consistentwiththebenchmarkresultsreportedinBhartietal[3].Theinßuenceofthetypeof non-NewtonianßuidorthepowerlawindexnonthestreamlinepatternandtheßowÞeldata higher Re of100isdisplayedinFig.2.16.Figure2.17showssnapshotsofstreamlinecontoursat di!erenttimeinstantsforthecasewith Re =150and n =1 . 4.Thisdemonstratesthattheßow becomesunsteadyforthiscase. ItisevidentthattheßowandthewakepatternissigniÞcantlya!ectedbythepowerlawindex n .Ingeneral,thewakeregionisfoundtogrowwith n ,withtheshearthickeningßuidexhibitingaconsiderablylongerrecirculatingzonewhencomparedtotheshearthinningßuid.Forsuch moderate Re ,thenonlinearityintheßowstemsfrombothinertialandviscoustermsinthemomentumequation.Inthecaseofshearthinningßuid,viscouse!ectsdominateevenfaraway formthecylindercontributingtoasmallerrecirculatingzone,while,bycontrast,forshearthickeningßuid,inertiale!ectsmainlycontributetotheoverallßowbehavior.

PAGE 57

42 (a)Re=5,n=0.6 (b)Re=10,n=0.6 (c)Re=20,n=0.6 (d)Re=30,n=0.6 FIGURE2.15:Streamlinecontoursoftheßowofashearthinningßuid( n =0 . 6)atdi!erent valuesReynoldsnumberReof5 , 10 , 20 , 20and30overanunconÞnedcircularcylindercomputed usingthecascadedLBmethod.

PAGE 58

43 (a)Re=100,n=0.6 (b)Re=100,n=0.6 (c)Re=100,n=1 (d)Re=100,n=1 (e)Re=100,n=1.4 (f)Re=100,n=1.4 FIGURE2.16:Instantaneousstreamlineandvelocitycontoursoftheßowofapowerlawßuid atRe=100anddi!erentvalueofthepowerlawindex n of0 . 6 , 1 . 0and1 . 4overanunconÞned circularcylindercomputedusingthecascadedLBmethod.

PAGE 59

44 X Y 400 600 800 100 200 300 400 (a)t=20000 X Y 400 600 800 100 200 300 400 (b)t=100000 X Y 400 600 800 100 200 300 400 (c)t=200000 X Y 400 500 600 700 800 900 0 200 400 (d)t=300000 FIGURE2.17:Streamlinecontoursatdi!erenttimeinstantsoftheßowofapowerlawßuidat Re =100, n =1 . 4overanunconÞnedcircularcylindercomputedusingthecascadedLBmethod.

PAGE 60

45 TABLE2.1:Dependenceofthetotaldragcoe"cient CD ontheReynoldsnumberReandthe powerlawindex n duetotheßowofpowerlawßuidspastanunconÞnedcirculatorresultsobtainedusingthecascadedLBmethodfordi!erentgridresolutionswiththebenchmarkdataof Bhartietal(2006). Source Re =5 Re =10 Re =20 Re =30 Re =40 Mesh( N # M ) n C D C D C D C D C D Bharti etal [3] 0 . 6 4 . 6505 2 . 909 1 . 95555 1 . 5409 1 . 3368 GridA1 1201 # 301 0 . 6 4 . 83232 4 . 386302 2 . 81125 1 . 864631 1 . 562475 GridA2 2001 # 501 0 . 6 4 . 819913 3 . 232486 2 . 170854 1 . 713212 1 . 419786 GridA3 4001 # 1001 0 . 6 4 . 807917 2 . 991372 1 . 952909 1 . 576768 1 . 368863 Bharti etal [3] 1 . 0 4 . 0374 2 . 7652 2 . 0455 1 . 6869 1 . 5249 GridA1 1201 # 301 1 . 0 4 . 811624 3 . 32588 2 . 385816 1 . 748944 1 . 55013 GridA2 2001 # 501 1 . 0 4 . 476988 2 . 958855 2 . 116828 1 . 722815 1 . 528644 GridA3 4001 # 1001 1 . 0 4 . 025268 2 . 79901 2 . 048933 1 . 710154 1 . 526487 Bharti etal [3] 1 . 4 3 . 5212 2 . 6568 2 . 0489 1 . 7741 1 . 6061 GridA1 1201 # 301 1 . 4 3 . 861923 2 . 882866 2 . 199001 1 . 856652 1 . 654297 GridA2 2001 # 501 1 . 4 3 . 747664 2 . 724362 2 . 089265 1 . 846678 1 . 64216 GridA3 4001 # 1001 1 . 4 3 . 60933 2 . 704615 2 . 079232 1 . 795661 1 . 620551 Akeyparameterofquantitativeinterestisthehydrodynamicforceduetosuchcomplexrheologicalßowsonthecircularcylinder.Thetotaldragforce F D perunitcylinderlengthisdue totheviscousandthepressureforcesactingonthecylinderduetonon-Newtonianßuidmotion,andcanberewrittenintermsofthedimensionlesstotaldragcoe"cient C D asfollows: C D = F D / ( 1 2 # U 2 " D ).Adimensionalanalysisrevealsthatitis,ingeneral,afunctionof Re and n ,i.e., C D = C D ( Re,n ).Inthiswork,weemploythemomentumexchangemethod[61] toobtainthedragforce F D (andhence C D )forvarioussetsof Re and n .Itisfoundthatthe computationsof C D issigniÞcantlyinßuencedbythechoiceofthegridresolution.Hence,simulationswereperformedforthreedi!erentchoicesofgridresolutions:lowresolution(GridA1) 1201 ! 1301,mediumresolution(GridA2)2001 ! 501andhighresolution(GridA3)4001 ! 1001. Thecoe"cientofdrag C D iscomputedusingthecascadedLBschemefortheabovegridresolutionsfor Re =5 , 10 , 20 , 30and40,whereforeachcasethepowerlawindex n isvariedacross 0 . 6 , 1 . 0and1 . 4.TheconvergedresultsalongwiththebenchmarkdataofBharti ?? .Ingeneral, thecomputedresultsobtainedwiththeÞnestresolution(GridA3)arefoundtobeinagood agreementwiththebenchmarkdata[3].As Re increases,foraÞxed n ,thecoe"cientoftotaldrag C D decreases.Interestingly,for Re ( 15, C D isfoundtodecreasewithanincreasein thepowerlawindex n ;ontheotherhand,theoppositetrendisobservedfor15
PAGE 61

46 whereanincreasein n isaccompaniedbyanincreaseinthetotaldragcoe"cient C D .Sucha non-monotonicdependenceof C D on n dependingonthe Re rangeasreportedinbenchmarkresults[3]iswellreproducedbythecascadedLBmethodbasedoncentralmoments.Theseobservationsaresummarizedasaplotbetween C D and Re ,parameterizedfor n ,encompassingshear thinning,Newtonianandshearthickeningßuids,andpresentedinFig.2.18. FIGURE2.18:ComparisonbetweenthenumericalresultsofthecascadedLBmethodand benchmarkdata[3]forthetotaldragcoe"cient C D asafunctionofReynoldsnumber Re atdifferentvaluesofthepowerlawindex n . 2.4.4 Non-Newtonianßowofpowerlawßuidsoverapairofcylindersintandemarrangement Asalastbenchmarkprobleminthisstudy,weconsiderbothsteadyandunsteadyßowofpower lawßuidsaroundtwocircularcylindersinatandemconÞguration.Aschematicofthegeometric arrangementofthisproblemisillustratedinFig.2.1d.Eachcylinderisconsideredtobeofthe samediameter D andthespacingbetweenthem L inthetandemarrangementisrepresentedby meansofthegapratio G = L/D .ThisßowproblemischaracterizedbymeansoftheReynolds number Re = # D n U 2 ! n " /µ p ,where U " istheimposedfreestreamvelocity,powerlawindex n andthegapratio G .Asintheprevioussection,weconsideranunconÞnedßowconÞguration, andtheinßowandoutßowconditions,free-slipconditionsonthesides,aswellastheno-slip boundaryconditionsonthecurvedcylinderwallsasimposedasdiscussedearlier.Benchmark

PAGE 62

47 resultsbasedonthenumericalsolutionusingasÞnitevolumemethodarepresentedinPatilet al[1],whichwillbeusedasareferenceformakingcomparisonsinthisstudy. Weconsideragridresolutionof2141 ! 2141toperformnon-NewtonianßowoverapairofcylindersunderavarietyofconditionsusingthecascadedLBmethod.Inthisregard,thechoices madeforthemagnitudesofthecharacteristicparametersareasfollows:Reisvariedbetween 5to150,thegapratio G istakenaseither2or4,thepowerlawindex n =0 . 4 , 1 . 0and1 . 8, spanningacrossbothshearthinningandshearthickeningßuids.Figure2.19presentsthesteady streamlinecontoursforalow Re of5atgapratios G of2and4,whereineachcase,thepower lawindex n isconsideredequalto0 . 4 , 1 . 0and1 . 8. Itcanbeseenthatatthisrelativelylow Re ,sincetheinertiale!ectsarenegligible,thegaprate G playsaminorroleforNewtonianandshearthickeningßuids;ontheotherhand,forshear thinningßuids( n =0 . 4),sinceviscouse!ectsdominate,itisfoundthatevenforsuchlow Re , thereisaninteractionbetweenthetwocylinderswiththehighergapratio( G =4)casealready exhibitingarecirculationzonefortheupstreamcylinder,whilesuchapatternisabsentforthe G =2case(seeFigs.2.19aand2.19b).Inaddition,itisinterestingtoseethatwhiletheßow issymmetricforNewtonianßuids(seeFigs.2.19cand2.19d),thesymmetryisbrokenforshear thickeningßuidswiththepresenceofalargerandasmallerrecirculatingstandingeddiesbehind theupstreamanddownstreamcylinder,respectively(seeFigs.2.19eand2.19f). Thestreamlinecontours Re =10forsimilarvariationsin G and n asinthepreviouscaseare showninFig.2.20.Forthishigher Re case,whenthegapratio G is2,itcanbeseenthatthere isasigniÞcantinterferencee!ectswithinthegapofthetwocylindersfortheNewtoniancase ( n =1)aswellastheshearthickeningcase( n =1 . 8).itisalsofoundthattherecirculation zonegrowswithanincreasein n .Infact,for G =2,therecirculatingeddiesoftheupstream cylinderiscompletelyconÞnedwithinthegap(seeFig.2.20e).Figure2.21displaysthestreamlinepatternsfor Re =20fortheinteractingßowbetweentwocylindersinatandemarrangement

PAGE 63

48 (a)Re=5,G=2,n=0.4 (b)Re=5,G=4,n=0.4 X Y 950 1000 1050 1100 1150 1200 950 1000 1050 1100 1150 1200 (c)Re=5,G=2,n=1.0 (d)Re=5,G=4,n=1.0 X Y 950 1000 1050 1100 1150 1200 950 1000 1050 1100 1150 1200 (e)Re=5,G=2,n=1.8 (f)Re=5,G=4,n=1.8 FIGURE2.19:StreamlinecontoursfortheßowofpowerlawßuidsatRe=5atdi!erentvalues ofthegapratio G =2and4andofthepowerlawindex n =0 . 4 , 1 . 0and1 . 8ineachcaseoveran unconÞnedpairofcylindersintandemarrangementcomputedusingthecascadedLBmethod.

PAGE 64

49 (a)Re=10,G=2,n=0.4 (b)Re=10,G=4,n=0.4 (c)Re=10,G=2,n=1.0 (d)Re=10,G=4,n=1.0 (e)Re=10,G=2,n=1.8 (f)Re=10,G=4,n=1.8 FIGURE2.20:Streamlinecontoursfortheßowofpowerlawßuidsat Re =10atdi!erentvalues ofthegapratio G =2and4andofthepowerlawindex n =0 . 4 , 1 . 0and1 . 8ineachcaseoveran unconÞnedpairofcylindersintandemarrangementcomputedusingthecascadedLBmethod.

PAGE 65

50 for G =2and4and n =0 . 6 , 1 . 0and1 . 4.Withhigher Re thaninthepreviouscases,relatively largerinertiale!ectscanbeexpected,especiallyfortheshearthickeningßuidcase( n =1 . 8).In particularitcanbenoticedthattheupstreamcylinder'srecirculatingzonescharacterizedbythe presenceoftwostandingvorticesgrowlargeenoughtobecompletelyconÞnedwithinthespacing betweenthecylinderforboth G =2and4(seeFigs.2.21eand2.21f,whereforthelattercase, asigniÞcantlyelongatedrecirculatingeddiesappear.Bycontrast,inthecaseofshearthinning ßuids( n =0 . 4),arelativelyweakandnarrowrecirculatingzoneappearswhen G =2behindthe upstreamcylinderduetotheinterferencee!ectfromthedownstreamcylinder.This,however, becomesevenlesssigniÞcantasthegapbetweenthecylindersisincreased(seeFig.2.21b). Finally,theinstantaneousstreamlinepatternsforamuchhigher Re of150arepresentedinFig.2.22. Forsucharelativelyhigh Re ,theßowbecomesunsteadyandvortexsheddingoccursbehindeach cylinderandthestrongwakeinterferencee!ectsineachcasewiththebreakingoftheup-down symmetryintheßowcanbeexpectedandareobserved.Forexample,fortheNewtonianßuid case( n =1),whileapairofrecirculatingeddiesarestillpresentfor G =2,whenthespacingis furtherincreased,vonKarmanvorticesshedfromtheupstreamcylinderbecomesmixedwiththe downstreamwakeviainertiale!ects.Ontheotherhand,inthecaseofshearthinningßuidscase ( n =0 . 6),therecirculatingpairofvorticaleddiesstillpersistevenforthehigherspacingratio of4,albeitwiththeirgrossdistortionduetosigniÞcantlongrangeviscouse!ectsatlowershear rates.Clearly,thereisacomplexinterplaybetweenthenonlinearitiesduetoinertialandviscous forcesinnon-Newtonianßuidßows,stronggeometricalandwakeinterferenceforhigh Re power lawßuidßowovertwocylindersinatandemarrangement. Wenowmakeaquantitativecomparisonbetweenthetotaldragcoe"cientsforvariousrangesof Re,G and n obtainedusingthecascadedLBmethodandreferencenumericalsolution[1].The totaldragcoe"cientfortheupstreamcylinder(1)isdeÞnedas C d 1 = F D 1 / < 1 2 # U 2 " D = ,where F D 1 ,isthetotalforceduetoviscousstressesandpressureactingonitssurface.Similarly,the totaldragcoe"cientinthedownstreamcylinder(2)isobtainedas C d 2 = F D 2 / < 1 2 # U 2 " D = .These

PAGE 66

51 (a)Re=20,G=2,n=0.4 (b)Re=20,G=4,n=0.4 (c)Re=20,G=2,n=1.0 (d)Re=20,G=4,n=1.0 (e)Re=20,G=2,n=1.8 (f)Re=20,G=4,n=1.8 FIGURE2.21:StreamlinecontoursfortheßowofpowerlawßuidsatRe=20atdi!erentvalues ofthegapratio G =2and4andofthepowerlawindex n =0 . 4 , 1 . 0and1 . 8ineachcaseoveran unconÞnedpairofcylindersintandemarrangementcomputedusingthecascadedLBmethod.

PAGE 67

52 (a)Re=150,G=2,n=0.4 (b)Re=150,G=4,n=0.4 (c)Re=150,G=2,n=1.0 (d)Re=150,G=4,n=1.0 (e)Re=150,G=2,n=1.8 (f)Re=150,G=4,n=1.8 FIGURE2.22:Instantaneousstreamlinecontoursfortheßowofpowerlawßuidsat Re =150 atdi!erentvaluesofthegapratio G =2and4andofthepowerlawindex n =0 . 4 , 1 . 0and 1 . 8ineachcaseoveranunconÞnedpairofcylindersintandemarrangementcomputedusingthe cascadedLBmethod.

PAGE 68

53 twoquantitieswhicharerelatedtothehydrodynamicforcesarecomputedusingthemomentum exchangemethodindicatedearlier. TABLE2.2:Dependenceofthetotaldragcoe"cients C d 1 and C d 2 onthegapratio G andthe powerlawindex n for Re =5duetotheßowofpowerlawßuidspastanunconÞnedpairofcircularcylindersinatandemarrangement.Comparisonbetweenthecomputedresultsobtained usingthecascadedLBmethodwiththebenchmarkdataofPatiletal[1]. Reynoldsnumber Re =5 Gapratio n 0 . 4 0 . 6 1 . 0 1 . 4 1 . 8 G =2 C d 1 4 . 188244 3 . 88316 3 . 407905 3 . 240358 2 . 966944 C d 1 4 . 1589 3 . 7886 3 . 3859 3 . 135 2 . 954 C d 2 2 . 578822 2 . 047688 1 . 248763 0 . 848029 0 . 563664 C d 2 2 . 5945 2 . 0389 1 . 2796 0 . 845 0 . 5892 G =4 C d 1 4 . 608887 4 . 137481 3 . 590397 3 . 185736 2 . 976830 C d 1 4 . 5518 4 . 0817 3 . 5153 3 . 1783 2 . 9596 C d 2 3 . 233944 2 . 500625 1 . 585037 1 . 033909 0 . 723546 C d 2 3 . 266 2 . 5263 1 . 6006 1 . 0819 0 . 9027 G =10 C d 1 4 . 980241 4 . 416323 3 . 834372 3 . 369549 3 . 023957 C d 1 4 . 9278 4 . 4037 3 . 7094 3 . 2971 3 . 0393 C d 2 4 . 119225 3 . 192998 2 . 141775 1 . 412249 0 . 912218 C d 2 4 . 2745 3 . 3171 2 . 1354 1 . 4508 0 . 904 TABLE2.3:Dependenceofthetotaldragcoe"cients C d 1 and C d 2 onthegapratio G andthe powerlawindex n for Re =10duetoßowofpowerlawßuidspastanunconÞnedpairofcircular cylindersinatandemarrangement.Comparisonbetweenthecomputedresultsobtainedusing thecascadedLBmethodwiththebenchmarkdataofPatiletal[1]. Reynoldsnumber Re =10 Gapratio n 0 . 4 0 . 6 1 . 0 1 . 4 1 . 8 G =2 C d 1 2 . 596081 2 . 519294 2 . 465175 2 . 424163 2 . 391814 Cd 1 2 . 5522 2 . 5104 2 . 4711 2 . 4334 2 . 3931 C d 2 1 . 263788 1 . 037479 0 . 704358 0 . 490213 0 . 350302 C d 2 1 . 2705 1 . 0622 0 . 7392 0 . 5228 0 . 3771 G =4 C d 1 2 . 767161 2 . 675699 2 . 611985 2 . 492631 2 . 389944 Cd 1 2 . 7185 2 . 6322 2 . 5192 2 . 4451 2 . 3931 C d 2 1 . 7012819 1 . 402391 0 . 999164 0 . 686184 0 . 466835 C d 2 1 . 7159 1 . 4205 0 . 9953 0 . 7066 0 . 6307 G =10 C d 1 2 . 945272 2 . 867031 2 . 704029 2 . 528601 2 . 451890 C d 1 2 . 9014 2 . 7974 2 . 6345 2 . 5225 2 . 4406 C d 2 2 . 4235551 1 . 999466 1 . 410124 0 . 946226 0 . 664660 Cd 2 2 . 4034 2 . 0057 1 . 423 1 . 0169 0 . 6746 Table2.2presentsthecomputedtotaldragcoe"cients C d 1 and C d 2 for G =2 , 4and10,when n =0 . 4 , 0 . 6 , 1 . 0 , 1 . 4and1 . 8obtainedusingthecascadedLBschemefor Re =5.Inaddition, thedragcoe"cients C d 1 and C d 2 forsimilarrangesofvariationofparameters G and n asabove for Re =10and20arepresentedinTables2.3and2.4,respectively.Ingeneral,itisfoundthe totaldragcoe"cientsforboththeupstreamanddownstreamcylinderscomputedusingthecas-

PAGE 69

54 TABLE2.4:Dependenceofthetotaldragcoe"cients C d 1 and C d 2 onthegapratioGandthe powerlawindex n for Re =20duetotheßowofpowerlawßuidspastaunconÞnedpairofcircularcylindersinatandemarrangement.Comparisonbetweenthecomputedresultsobtained usingthecascadedLBmethodwiththebenchmarkdataofPatiletal[1]. Reynoldsnumber Re =20 Gapratio n 0 . 4 0 . 6 1 . 0 1 . 4 1 . 8 G =2 C d 1 1 . 684060 1 . 754817 1 . 8531 1 . 9261 2 . 027656 C d 1 1 . 656 1 . 733 1 . 8531 1 . 9261 1 . 9673 C d 2 0 . 587634 0 . 512071 0 . 375855 0 . 264022 0 . 191637 Cd 2 0 . 6054 0 . 5297 0 . 3938 0 . 2871 0 . 2032 G =4 C d 1 1 . 710787 1 . 786806 1 . 887693 1 . 936106 2 . 048884 C d 1 1 . 704 1 . 7613 1 . 8555 1 . 9561 1 . 9517 C d 2 0 . 892138 0 . 781178 0 . 56704 0 . 394438 0 . 293014 C d 2 0 . 9253 0 . 8025 0 . 5923 0 . 4256 0 . 2956 G =10 C d 1 1 . 821105 1 . 883853 1 . 928191 1 . 960231 1 . 986042 C d 1 1 . 8069 1 . 8547 1 . 9211 1 . 9561 1 . 9699 C d 2 1 . 390043 1 . 231276 0 . 915376 0 . 634812 0 . 446591 C d 2 1 . 4324 1 . 2578 0 . 9502 0 . 698 0 . 4995 cadedLBmethodareingoodquantitativeagreementwiththebenchmarknumericalresultsof Patiletal[1]forvariousrangesof Re,G and n considered.ItcanbeseenthatforanyÞxed Re and n considered,anincreaseinthegapratio G correspondstoincreaseinthedragcoe"cients ofboththecylinders( C d 1 and C d 2 ).Interestingly,inanalogywiththesinglecylindercasepresentedearlier,increaseinthepowerlaw n leadstodecreaseintotaldragcoe"cient C d 1 ofthe upstreamcylinderforthelower Re conditions(i.e., Re =5,10);bycontrast,forthehigher Re case( Re =20),thetotaldragcoe"cient C d 1 oftheupstreamcylinderincreasesas n isincreased. Ontheotherhand,forthedownstreamcylinder,thetotaldragcoe"cient C d 2 isalwaysfound todecreasewithanincreasein n foranyÞxed G ,regardlessoftheRechosen.Furthermore,for anyÞxed G and n ,boththetotaldragcoe"cients C d 1 and C d 2 decreaseasthe Re isincreased. Finally,whencomparedtosinglecylindercase,itsdragcoe"cient C D (seeTable ?? )isalways foundtobehigherthanthedragcoe"cientfortheupstreamcylinder C d 2 ,which,inturn,is greaterthanthatforthedownstreamcylinderforthetwincylindercaseinatandemarrangementunderotherwisesimilarconditions(seeTable2.2,2.3,and2.4).Thatis, C D >C d 1 >C d 2 , whicharisefromtheconÞnementandwakeinterferencee!ectsinthelattercase.AlltheseobservationsareinqualitativeandquantitativeagreementwiththereferencenumericalresultsofPatil etal[1].

PAGE 70

55 2.5 Comparisonofnumericalstabilityofdi!erentcollisionmodels Wenowmakeacomparisonofthenumericalstabilityofdi!erentLBformulationforthesimulationsofthe2Dliddrivencavityßowatdi!erentgridresolutions.Inthisregard,weconsider thefollowingcollisionmodels:thepopularsinglerelaxationtime(SRT)model,thecascaded MRTmodelwithoutincorporatingcubicvelocityGalileaninvariance(GI)correctionsandthe cascadedMRTmodelwithcubicvelocityGIcorrections.ThecascadedLBschemeswithand withoutGIcorrectionsarereferredtoascentralmomentIandcentralmomentII,respectively, inwhatfollows.IntheSRTmodel,therelaxationtimethatcontrolsthelocalßuidviscosityis referredtoas ! ,whileforthecascadedLBformulations,wetake , 4 = , 5 = 1 + ,whiletherelaxationtimesforthebulkviscosityandhigherordermomentsaresettounityforsimplicity ( i.e.,, 3 = , 6 = , 7 = , 8 =1).Witheachofthecollisionmodels,forachosengridresolution andaÞxedlidvelocity,therelaxationparameter ! isdecreaseduntilthecomputationbecomes unstable.Figure2.23presentstheresultsintermsofthemaximumReynoldsnumberachieved forthepowerlawindex n =0 . 5 , 1 . 0and1 . 5(encompassingshearthinning,Newtonianandshear thickeningßuids)withtheuseofdi!erentcollisionmodelsforgridresolutionsof101 2 , 201 2 and 401 2 ,respectively.Ingeneral,itcanbeseenthatthecascadedLBformulationsissigniÞcantly morestablewhencomparedtotheSRTLBscheme.Itisfoundthat,ingeneral,thesimulation ofshearthinningßuidsismorechallengingthanthatofothertypeofßuids.Inaddition,itisevidentthatwiththeuseofGIcorrections(centralmomentII),themaximumachievableReynolds numberwiththecascadedLBmethodisconsiderablygreaterthenthatwithoutsuchGIcorrections(centralmomentsI)forsimulationofnon-Newtonianßows. 2.6 SummaryandConclusions Inthiswork,acascadedlatticeBoltzmann(LB)formulationforthecomputationofnon-Newtonian ßuidmotionisdiscussed.Thecollisionstepinthisapproachisbasedontherelaxationofcentral

PAGE 71

56 401x401 201x201 101x101 0 200 400 600 800 1000 1200 1400 1600 1800 2000 Maximum Re Grid Resolution SRT Central Moment I Central Moment II (a) n =0 . 5 401x401 201x201 101x101 0 2000 4000 6000 8000 10000 12000 Maximum Re Grid Resolution SRT Central Moment I Central Moment II (b) n =1 . 0 401x401 201x201 101x101 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 x 10 4 Maximum Re Grid Resolution SRT Central Moment I Central Moment II (c) n =1 . 5 FIGURE2.23:ComparisonofthemaximumReynoldsnumberfornumericalstabilityofSRT LBscheme,cascadedLBschemewithoutGIcorrections(centralmomentI),andcascadedLB schemewithGIcorrections(centralmomentsII)forsimulationofnon-Newtonianliddrivencavityßowwithn=0.5,1.0,1.5.

PAGE 72

57 momentsofdi!erentordersatdi!erentrates.Inordertorepresentthenon-Newtonianßuidßow behavior,weconsideredtheirrepresentationvianonlinearconstitutiverelationsofpowerlawand Binghamßuids.Theformercanmodelshearthinningaswellasshearthickeningbehavior,while thelatterrepresentsviscoplasticßuidsendowedwithayieldstress.Thee!ectofthesemodels areincorporatedintothecascadedLBmethodvialocalvariationoftherelaxationtimesofthe secondordermomentsrelatedtothemagnitudeoftheshearrate.Thelatterisexpressedasan invariantofthestrainratetensor,whichisobtainedlocallyfromnon-equilibriumsecondorder moments,whileaccountingfortheadditionalcorrectionsnecessarytoeliminatethecubicvelocity errorsonthestandardD2Q9lattice. ThecascadedLBmethodisvalidatedagainstavailableanalyticaland/orpriornumericalsolutionsfornon-Newtonianchannelßow,lid-drivencavityßow,ßowpastasinglecylinderaswellas overapairofcylindersintandemarrangement.Thecomputedresultsarefoundtobeinvery goodagreementwithbenchmarksolutionsforalltheproblemsconsidered.Inparticular,the complexbehaviorofthedragcoe"cientonthepowerlawindexdependentontherangeofthe Reynoldsnumbersfortheßowoverasinglecylinderandwakeinteractionandinterferencee!ects fortheßowoverapairofcylindersarewellreproducedbythepresentLBapproach.Inaddition,asecondorderspatialconvergenceofthemethodisdemonstratednumerically.Finally,an improvementinthenumericalstabilityduetotheuseofcascadedLBscheme,especially,byincludingcorrectionstoeliminatecubicvelocityerrors,whencomparedtoothercollisionmodelfor simulationsofnon-Newtonianßuidßowsisshown.

PAGE 73

58 CHAPTERIII THREE-DIMENSIONALCASCADEDLBMFORNON-NEWTONUAN FLUIDFLOWS 3.1 Introduction Inthelastfewdecades,thelatticeBoltzmannMethod(LBM)hasbecomeapreferredmethod forsimulatingcomplicatedphysical,chemical,andßuidmechanicsproblems[12],[16],[62].Itis akinetic-basedapproachforßuidßowcomputations.Thisiswhyitisespeciallyusefulforcomputingßuidßowswithmultiplecomponentsinvolvinginterfacialdynamicsnonlinearconstitutive modelsandcomplexboundaries[63].OtherCFDmethodshavemanyissueswhenitcomesto solvingmacroscopichydrodynamicequations,particularlyforßuidinterfacesandphasetransitions.BecausetheLBMliesonthescaleofmicroscopicormicroscopiclevel,challengesthathave beenencounteredwithusingconventionalCFDmethodsarenotfoundwiththeLBMmethod [ ? ,7].ThelatticeBoltzmannEquation(LBE)canbeobtainedthroughseveralways.Onecan derivetheLBEthroughdiscretevelocitymodelsortheBoltzmannkineticequation.Inaddition, thereareseveralapproachestoderiveNavier-StokesequationsfromtheLBE.WithChapmanEnskogexpansionbeingmorepopularamongtherest,otherapproachesincludeasymptoticexpansion,extendedTaylorseriesexpansionandorderofmagnitudeanalysis. ThelatticeBoltzmannmethodsarecomprisedoftwofundamentalsteps,whicharethestreamingstepandcollisionstep.ThestreamingstepisthsameinvariousmodelsofLBE.However, becausethecollisionstepismorecomplicated,andresearcheshavebeendevotingconsiderable e!ortsintoÞndingthemostsuitablecollisionmodelforthelatticeBoltzmannmethod.Among thosedi!erentcollisionmodels,thesimplestandcommonlyusedistheso-calledsingle-relaxationtime(SRT)model[63],[18],whichisbasedontheBhatnagar-Gross-Krook(BGK)approxima-

PAGE 74

59 tion[64].Ontheotherhand,amultiple-relaxation-time(MRT)modelhasbeendevelopedin ordertoimprovenumericalstability[19].MRThasbeenconÞrmedasamorestablecollision modelinvariousproblems.IntheMRTcollisionmodel,variousmomentsarerelaxedtotheir equilibriumstatesatdi!erentratesduringthecollisionstep.Recently,anewclassofcollision model,theso-calledcascadedLBM,hasbeenintroducedbyGeieretal.[65].Lateron,Asinari[32] reinterpretedthisapproachbasedongeneralizedequilibrium.PremnathandBanerjee[20]incorporatedforcingtermsincascadedlatticeBoltzmannapproachofthemethodofcentralmoments andsystematicallyderivedthecascadedLBMapproach.Insuchapproach,Galileaninvarianceis enforcedtothelatticeBoltzmannequation(LBE)basedontherelaxationofvariouscentralmoments.Thisinvolvescomputingmomentswhichareshiftedbythemacroscopicßuidvelocity.In anotherwords,themomentsareprescribedinamovingframeofreference.Comparatively,the momentsinthepriorapproachesarecomputedinarestframeofreference,whicharetermedas therawmoments. Invariousrecentstudies[24,25],thecascadedLBMbasedoncentralmomentsandmultiplerelaxationtimeshasbeenshowntobesigniÞcantlymorestablewhencomparedtotheSRTcollisionmodelbasedLBMforthesimulationofNewtonianßuidsßows.Threedimensionalcentral momentLBMsincludingforcingtermsforD3Q15andD3Q27latticesweresystematicallyderivedin[23].Avariantofsuchaformulationinvolvingdiscretemomentequilibriaratherthan thecontinuousMaxwellianequilibriaispresentedin[66].ApreconditionedcascadedLBMto acceleratesteadyßowsimulationswithimprovedaccuracywasrecentlydevelopedby[29].Furthermore,thecascadedLBMhasrecentlybeenextendedtosimulatethermalconvectiveßowsin two-dimensions(2D)[67]andthree-dimensions(3D)[29].Inaddition,amodiÞedforcingscheme formulationandimplementationofthecascadedLBMhasbeenpresentedin[67]. Complexßuidßowswithnon-Newtonianconstitutiverelationsrepresentanimportantclassof ßowswithnumerousapplicationsinvariousscenariosincludinginengineering,materialsand foodprocessing,andgeophysicalprocesses.TheextensionofLBMtonon-Newtonianßuidßows

PAGE 75

60 hasreceivedsigniÞcantbutlimitedattentionsofar,whichhasbeenreviewedin[15],andsome examplesofsuchstudiesincludetheworksof[35,45,46].Whilemostofthepriorstudieshave demonstratedtheapplicationsofdi!erentLBMformulationsto2Dnon-Newtonianßuidßows, verylittlefocushasbeengiventothedevelopmentandvalidationofLBM,particularlywith usingadvancedcollisionmodels,for3Dnon-NewtonianßuidsinvolvingwelldeÞnedbenchmark problems.Indeed,duetocomplexityofhandlingnonlinearrheologicalbehaviorviamoregeneral constitutiverelationswithattendantnumericalstabilityissuesandtheneedforhighgridresolutionswithmuchhighercomputationaldemands,thereareonlyveryfewstudiesin3Devenin thecontextofcomputationalßuiddynamics(CFD)usingconventionalnumericalschemes.For thewell-deÞnedproblemof3Dnon-Newtoniancubiccavityßow,Refs.[68,69]presentedsomeresultsatcoarsegridresolutionsand,recently,Ref.[70]reportedbenchmarkqualityresultsusinga fractionalstepbasedÞnitevolumemethod. Inthepresentwork,wepresentnewcascadedLBMwithforcingtermsin3Dforthesimulation ofnon-Newtonianßowsrepresentedbypowerlawßuidbasedrheology.Inthisregard,weconstructacollisionmodelbasedoncentralmomentsandmultiplerelaxationtimes,onathreedimensionalnineteenvelocity(D3Q19)lattice.Here,therelaxationtimesforthesecondorder momentsthatemulatetheviscousbehavioroftheßuidareadjustedbythelocalshearrateand parametrizedbytheconsistencycoe"cientandthepowerlawindexofthegeneralizedNewtonian(i.e.,powerlaw)ßuid.Theforcingtermsthatmodelthelocalimpressedforcesarerepresentedassourcetermstothe3DcascadedLBMthatareconstructedusingtheire!ectonappropriatecentralmomentsandthendiscretizedusingatrapezoidalruletoachievesecondorderaccuracy.Wethenpresentavalidationstudyofour3DcascadedLBschemeforthenon-Newtonian ßowinachannel,ductandacubiccavity.Inparticular,forthelattercaseinvolvingcomplex nonlinearßuidmotionwithpower-lawrheologyinacubiccavity,wewillmakeadirectcomparisonofthevelocityproÞlesagainsttherecentbenchmarksolution[70]forvariousReynoldsnumbersandpowerlawindexmagnitudesencompassingbothshearthinningandshearthickening

PAGE 76

61 ßuidbehavior.Finally,asanadditionalobjective,wepresentacomparativestudyofournew3D cascadedLBMformulationagainsttheSRTformulationforsimulationofnon-Newtonianßows anddemonstrateanimprovementinnumericalstabilityachievedbytheformer,whencompared tothelatterapproach. Thispaperisorganizedasfollows.Inthenextsection(Sec.3.2),wepresentanoverviewofthe macroscopicgoverningequationsforgeneralizednon-Newtonianßuidßowsandlisttheattendantnon-linearconstitutiverelation.AdetailedexpositiontheforcingtermsfortheD3Q19latticeisgiveninSec.3.3.Thisincludesthechangesofdi!erentmomentsundercascadedcollision aswellasduetoimpressedforces,andthelocalrelaxationtimesforthesecondordermoments parameterizedforthepowerlawßuid.Additionaldetailsrelatedtothisderivationaregivenin theappendices(AppendixA-2.5).Validationstudiesinvolvingvariouscanonicalnon-Newtonian ßuidßowproblems,includingthoseinacubiccavityarepresentedanddiscussedinSec.3.4.Section3.5makesanumericalstabilitycomparisonofour3DcentralmomentLBMagainstaSRT LBMforsimulationofpowerlawßuids.TheconclusionsareÞnallysummarizedinSec.3.6. 3.2 Governingequationsforthree-dimensionalgeneralizednon-Newtonianßuid ßows Inthepresentwork,weconsiderthesimulationofthree-dimensional(3D)generalizedNewtonian ßows(GNF).Themathematicalmodel,whichisgivenintermsofthecontinuityandmomentum equations,canbewrittenas %# % t + % % x j ( # u j )=0 , (3.1a) % % t ( # u i )+ % % x j ( # u i u j )= " % P % x i + %& ij % x j + Fi, (3.1b) where # and u i arethelocalßuiddensityandthevelocityÞeld,respectively,with i # ( x,y,z ). Here, P , F i and & ij representthepressure,Cartesiancomponentoftheimposedbodyforce,and

PAGE 77

62 viscousstresstensorinducedwithintheGNF,respectively.Theviscousstresstensorcanberepresentedas & ij = µ ( | ú ' ij | )ú ' ij . (3.2) whereú ' ij istheshearratetensor.Itisrelatedtothestrainratetensor S ij by ú ' ij =2 S ij ,S ij = 1 2 ( % x j u i + % x i u j ) , (3.3) and µ ( | ú ' ij | )isassumedtobeane!ectiveorapparentdynamicviscosityoftheGNF,wherethe magnitudeoftheshearrate | ú ' ij | isrelatedtothesecondinvariantofthesymmetricstrainrate tensor S ij as | ú ' ij | = # 2 S ij S ij . (3.4) Here,theusualconventionofthesummationoftherepeatedindicesisassumed.Thus, | ú ' ij | = C 2 ( S 2 xx + S 2 yy + S 2 zz +2( S 2 xy + S 2 yz + S 2 xz ) ) .Inthisstudy,weconsidertheconstitutiverelationsforthepowerlawßuidasarepresentative ofaclassofGNFfortheirnumericalsolutionusingthe3DcascadedLBmethodpresentedinthe nextsection.Theconstitutiverelationforthepowerlawßuidcanbedescribedmathematically asfollows: & ij = µ p | ú ' ij | n ! 1 ú ' ij , (3.5) wherethemodelparameter µ p andtheexponent n areknownastheconsistencycoe"cientof thepowerlawßuidandpowerlawindexofßuid,respectively.Basedonthevaluesofpower-law index n ,thefollowingthreedi!erenttypesofnon-Newtonianßuidscanberepresented: n< 1 correspondstoshearthinningorpseudo-plasticßuid,whereas n> 1modelsshearthickeningßuids,and n =1reducestotheNewtonianconstitutiverelation.ComparingEqs.(3.2)and(3.5), thee!ectiveviscosityforthepower-lawmodelisgivenbythefollowingexpression µ power | ú ' ij | = µ p | ú ' ij | n ! 1 . (3.6)

PAGE 78

63 3.3 Three-dimensionalcascadedLBMfornon-Newtonianßuidßows 3.3.1 SelectionofBasisVectorsforMoments Forsimplicity,withoutlosinggenerality,weconsiderthethree-dimensionalnineteenvelocity (D3Q19)latticemodel,whichisshowninFig.3.1.Theparticlevelocityforthislatticemodel ") e ! maybewrittenas ") e ! = F G G G G G H G G G G G I (0 , 0 , 0) , ) =0 ( ± 1 , 0 , 0) , (0 , ± 1 , 0) , (0 , 0 , ± 1) , ) =1 , ááá , 6 ( ± 1 , ± 1 , 0) , ( ± 1 , 0 , ± 1) , (0 , ± 1 , ± 1) , ) =7 , ááá , 18 (3.7) Here,weuseGreekandLatinsubscriptsfortheparticlevelocitydirectionsandCartesiancoordinatedirections,respectively.BydeÞnition,therawmomentsintheLBMaredeÞnedinterms ofthedistributionfunction f ! as > 18 ! =0 e m ! x e n ! y e p ! z f ! ,where m , n and p areintegers.Forconvenience,weconsidertheDirac'sbra-ketnotationinthispapertorepresentthebasisvectors. Now,wedeÞneasetofthefollowingnineteennon-orthogonalbasisvectorsbasedonthecombinationofthemonomials e m ! x e n ! y e p ! z inanascendingorder. | T 0 $ = | # $&|| ") e ! | 0 $ , | T 1 $ = | e ! x $ , | T 2 $ = | e ! y $ , | T 3 $ = | e ! z $ , | T 4 $ = | e ! x e ! y $ , | T 5 $ = | e ! x e ! z $ , | T 6 $ = | e ! y e ! z $ , | T 7 $ = | e 2 ! x " e 2 ! y $ , | T 8 $ = | e 2 ! x " e 2 ! z $ ,

PAGE 79

64 | T 9 $ = | e 2 ! x + e 2 ! y + e 2 ! z $ , | T 10 $ = | e ! x < e 2 ! x + e 2 ! y + e 2 ! z = $ , | T 11 $ = | e ! y < e 2 ! x + e 2 ! y + e 2 ! z = $ , | T 12 $ = | e ! z < e 2 ! x + e 2 ! y + e 2 ! z = $ , | T 13 $ = | e ! x < e 2 ! y " e 2 ! z = $ , | T 14 $ = | e ! y < e 2 ! z " e 2 ! x = $ , | T 15 $ = | e ! z < e 2 ! x " e 2 ! y = $ , | T 16 $ = | e 2 ! x e 2 ! y + e 2 ! x e 2 ! z + e 2 ! y e 2 ! z $ , | T 17 $ = | e 2 ! x e 2 ! y + e 2 ! x e 2 ! z " e 2 ! y e 2 ! z $ , | T 18 $ = | e 2 ! x e 2 ! y " e 2 ! x e 2 ! z $ . Here,thecompenentsofthebasisvectors | # $ , | e ! x $ , | e ! y $ and | e ! z $ ,whichareusedtobuildthe momentbasisaregivenby | # $& || ") e ! | 0 $ =(1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1)   , | e ! x $ =(0 , 1 , " 1 , 0 , 0 , 0 , 0 , 1 , " 1 , 1 , " 1 , 1 , " 1 , 1 , " 1 , 0 , 0 , 0 , 0)   , | e ! y $ =(0 , 0 , 0 , 1 , " 1 , 0 , 0 , 1 , 1 , " 1 , " 1 , 0 , 0 , 0 , 0 , 1 , " 1 , 1 , " 1)   , | e ! z $ =(0 , 0 , 0 , 0 , 0 , 1 , " 1 , 0 , 0 , 0 , 0 , 1 , 1 , " 1 , " 1 , 1 , 1 , " 1 , " 1)   . Theabove | T j $ , j =0 , 1 , 2 ,.., 18 , canbetransformedintothefollowingequivalentsetoforthogonalbasisvectorsusingthestandardGram-Schmidtprocedure,whichcanbewrittenas | K 0 $ = | # $& || ") e ! | 0 $ , | K 1 $ = | e ! x $ , | K 2 $ = | e ! y $ , | K 3 $ = | e ! z $ , | K 4 $ = | e ! x e ! y $ ,

PAGE 80

65 | K 5 $ = | e ! x e ! z $ , | K 6 $ = | e ! y e ! z $ , | K 7 $ = | e 2 ! x " e 2 ! y $ , | K 8 $ = | e 2 ! x + e 2 ! y + e 2 ! z $" 3 | e 2 ! z $ , | K 9 $ =19 | e 2 ! x + e 2 ! y + e 2 ! z $" 30 | # $ , | K 10 $ =5 | e ! x < e 2 ! x + e 2 ! y + e 2 ! z = $" 9 | e ! x $ , | K 11 $ =5 | e ! y < e 2 ! x + e 2 ! y + e 2 ! z = $" 9 | e ! y $ , | K 12 $ =5 | e ! z < e 2 ! x + e 2 ! y + e 2 ! z = $" 9 | e ! z $ , | K 13 $ = | e ! x e 2 ! y " e ! x e 2 ! z $ , | K 14 $ = | e 2 ! z e ! y " e ! y e 2 ! x $ , | K 15 $ = | e 2 ! x e ! z " e 2 ! y e ! z $ , | K 16 $ =21 | e 2 ! x e 2 ! y + e 2 ! x e 2 ! z + e 2 ! y e 2 ! z $" 16 | e 2 ! x + e 2 ! y + e 2 ! z $ +12 | # $ , | K 17 $ =3 | e 2 ! x e 2 ! y + e 2 ! x e 2 ! z " 2 e 2 ! y e 2 ! z $" 2 | 2 e 2 ! x " e 2 ! y " e 2 ! z $ , | K 18 $ =3 | e 2 ! x e 2 ! y " e 2 ! x e 2 ! z $" 2 | e 2 ! y " e 2 ! z $ . Thus,setofbasisvectorcanbeusedtoformanorthogonalmatrix | K $ as K =[ | K 0 $ , | K 1 $ , | K 2 $ , | K 3 $ , | K 4 $ , | K 5 $ , | K 6 $ , | K 7 $ , | K 8 $ | K 9 $ , | K 10 $ , | K 11 $ , | K 12 $ , | K 13 $ , | K 14 $ , | K 15 $ , | K 16 $ | K 17 $ , | K 18 $ ] (3.8) whosecomponentsaresummarizedinAppendixA.

PAGE 81

66 3.3.2 ContinuousCentralMoments:DistributionFunction,itsLocalAttractorandForcing First,letusnowdeÞnethecontinuouscentralmomentoftheparticledistributionfunction f shiftedbythelocalßuidvelocityoforder( m + n + p )as 0 $ x m y n z p = 1 " !" 1 " !" 1 " !" f ( ( x " u x ) m ( ( y " u y ) n ( ( z " u z ) p d ( x d ( y d ( z . (3.9) Here,andintherestofthispaper,weemploy"hat"overasymboltorepresentthevaluesinthe spaceofmoments.Thecorrespondingcentralmomentlocalattractorcanbewrittenas 0 $ at x m y n z p = 1 " !" 1 " !" 1 " !" f at ( ( x " u x ) m ( ( y " u y ) n ( ( z " u z ) p d ( x d ( y d ( z . (3.10) where f at istheattractororthelocalequilibriumofthedistributionfunction f .Similarly,the changesinthecontinuouscentralmomentsduetoimpressedforcescanbeformallywrittenas 0 % F x m y n z p = 1 " !" 1 " !" 1 " !" & f F ( ( x " u x ) m ( ( y " u y ) n ( ( z " u z ) p d ( x d ( y d ( z , (3.11) Here,& f F isthechangeinthedistributionfunctionduetoforcing,whichwillbeobtainedlater. Now,weconsiderthelocalMaxwell-Boltzmanndistributionfunctioninthecontinuousparticle velocityspace( ( x , ( y , ( z )astheattractor. f M & f M ( # , ") u,( x , ( y , ( z )= # 2 / c 3 2 s exp * + , " A ") ( " ") u B 2 2 c 2 s . / , (3.12) wherewechoose c 2 s =1 / 3.LetusnowdeÞnecontinuousMaxwelliancentralmomentsas 0 $ M x m y n z p = 1 " !" 1 " !" 1 " !" f M ( ( x " u x ) m ( ( y " u y ) n ( ( z " u z ) p d ( x d ( y d ( z . (3.13) Since f M isanevenfunction, 0 $ M x m y n z p * =0isnotzerowhen m , n and p areeven,andotherwise iszero.Usingthefactthat,wecanwritethisquantityintheincreasingorderofmomentsas: 0 $ M 0 = # , 0 $ M i =0 , 0 $ M ii = c 2 s # ,

PAGE 82

67 0 $ M ij =0 ,i * = j, 0 $ M ijj =0 ,i * = j, 0 $ M ijk =0 ,i * = j * = k, 0 $ M iijj = c 4 s # ,i * = j. (3.14) ItmaybenotedthatRef.[71]using 0 $ at x m y n z p = 0 $ M x m y n z p forallorderscanleadtosomeinconsistencyinrecoveringthemacroscopicßuidequations.Toproceedfurtherandobtain 0 $ at x m y n z p ,we nowdeÞnethefollowingpost-collisioncontinuouscentralmomentoforder( m + n + p ): ! 0 $ x m y n z p = 1 " !" 1 " !" 1 " !" ! f ( ( x " u x ) m ( ( y " u y ) n ( ( z " u z ) p d ( x d ( y d ( z . (3.15) Then,letusconsiderthefactorizedformforattractorsas[71] 0 $ at i = ! 0 $ i =0 , 0 $ at ij = ! 0 $ i ! 0 $ j =0 , 0 $ at iij = ! 0 $ ii ! 0 $ j =0 , 0 $ at ijk = ! 0 $ i ! 0 $ j ! 0 $ k =0 , 0 $ at iijj = ! 0 $ ii ! 0 $ jj . (3.16) Inordertorecoverthemomentumßuxandpressuretensorintheßuiddynamicalequations,the diagonalcomponentsofthesecond-ordercentralmomentsshouldbeobtainedfromtheMaxwellian centralmoments.Thatmeans,wehavetoset 0 $ at ii = c 2 s # .Then,the19independentcomponents ofthelocalfactorizedcentralmomentattractorscanbeformallywrittenas 0 $ at 0 =0 , 0 $ at x = 0 $ at y = 0 $ at z =0 , 0 $ at xx = 0 $ at yy = 0 $ at zz = c 2 s # , 0 $ at xy = 0 $ at xz = 0 $ at yz =0 , 0 $ at xyy = 0 $ at xzz = 0 $ at xxy = 0 $ at yzz = 0 $ at xxz = 0 $ at yyz = 0 $ at xyz =0 , 0 $ at xxyy = ! 0 $ xx ! 0 $ yy ,

PAGE 83

68 0 $ at xxzz = ! 0 $ xx ! 0 $ zz , 0 $ at yyzz = ! 0 $ yy ! 0 $ zz . (3.17) SimilarlyforthecontinuoussourcecentralmomentsduetoforceÞelds,onepossiblechoiceisobtainedbychoosingthatbasedonthelocalMaxwellian,i.e.& f F = !$ F # á ( !$ & ! !$ u ) c 2 s f M ,which,however,leadstoaliasinge!ectsforhigherordermoments[22].Tocircumventthisissue,asimple choiceinvolvesde-aliasinghigherordermomentswhilepreservingitsnecessarye!ectontheÞrstordercentralmoments[22]whichisextendedto3Dinthiswork. Inthiscase,thecontinuoussourcecentralmomentscanbespeciÞedas 0 % F x m y n z p = F G G G G G G G G H G G G G G G G G I F x ,m =1 ,n =0 ,p =0 F y ,m =0 ,n =1 ,p =0 F z ,m =0 ,n =0 ,p =1 0 , Otherwise. (3.18) 3.3.3 CascadedLattice-BoltzmannEquationwithForcingTerms First,letuswritethe3DcascadedlatticeBoltzmannequation(LBE)withforcingtermsbydeÞningadiscretedistributionfunctionsupportedbythediscreteparticlevelocityset ") e ! as f = | f ! $ =( f 0 ,f 1 ,f 2 ,...,f 18 )   ,asourceterm S = | S ! $ =( S 0 ,S 1 ,S 2 ,...,S 18 )   aswellasacollisionoperatoras ! c = | # c ! $ =(# c 0 , # c 1 , # c 2 ,..., # c 18 )basedoncentralmoments.ThecascadedLBE canbewrittenasadiscreteversionofthecontinuousBoltzmannequationbyintegratingalong theparticlecharacteristicsas[22] f ! ( ") x + ") e ! ,t +1)= f ! ( ") x,t )+ # c ! ( !$ x,t ) + 1 t +1 t S ! ( !$ x + !$ e ! , ,t + , ) d -. (3.19) Here,thecollisionterm# c ! canbewrittenintermsofboththeunknownchangeofmomentunder collision 0 g andtheorthogonalmatrixofthemomentbasisas[20] # c ! & # c ! ( f , 0 g )=( Ká 0 g ) ! , (3.20)

PAGE 84

69 where 0 g = | 0 g ! $ =( 0 g 0 , 0 g 1 , 0 g 2 ,..., 0 g 18 )   isthevectoroftheunknownchangeofmomentunder collisions,whichwillbederivedlater.Theßuiddynamicalvariables,i.e.,thelocaldensityand momentum,aredeterminedbyusingthedistributionfunctionas # = 18 " ! =0 f ! = % f ! | # $ , (3.21) # u i = 18 " ! =0 f ! e ! i = % f ! | e ! i $ ,i # x,y,z. (3.22) Weconsiderasemi-implicitdiscretizationoftheintegralofthesourcetermviathetrapezoidalin Eq.(3.19).Inordertoeliminatethee!ectsofimplicitness,atransformation f ! = f ! " 1 2 S ! is usedtoreduceitto[22] f ! ( ") x + ") e ! ,t +1)= f ! ( ") x,t )+ # c ! ( !$ x,t ) + S ! ( !$ x,t ) . (3.23) Thederivationofthechangeofmomentundercollision 0 g andthesourceterm S ! involvesmatchingthediscretecentralmomentsandthecontinuouscentralmomentsoftheattractors(equilibria)andthesourceterm.Thesewillbecarriedoutinwhatfollows. 3.3.4 VariousDiscreteCentralMomentsandGalileanInvarianceMatchingPrinciple Tofacilitatethedeterminationofthestructureofthecascadedcollisionoperator 0 g andthesource terms S ! ,wenowdeÞnethefollowing discrete centralmomentsofthedistributionfunction,its attractor,andthesourceterm,respectively: 0 * x m y n z p = % ( e ! x " u x ) m ( e ! y " u y ) n ( e ! z " u z ) p | f ! $ , 0 * at x m y n z p = % ( e ! x " u x ) m ( e ! y " u y ) n ( e ! z " u z ) p | f at ! $ , 0 & x m y n z p = % ( e ! x " u x ) m ( e ! y " u y ) n ( e ! z " u z ) p | S ! $ . (3.24) WealsodeÞnediscretecentralmomentsbasedonthepost-collision( ! f ! )andtransformed( f ! ) distributionfunctions,anditscombination ! f ! tofacilitatesubsequentcalculations: ! 0 * x m y n z p = % ( e ! x " u x ) m ( e ! y " u y ) n ( e ! z " u z ) p | ! f ! $ ,

PAGE 85

70 0 * x m y n z p = % ( e ! x " u x ) m ( e ! y " u y ) n ( e ! z " u z ) p | f ! $ , ! 0 * x m y n z p = % ( e ! x " u x ) m ( e ! y " u y ) n ( e ! z " u z ) p | ! f ! $ . (3.25) Owingtothetransformationofthedistributionfunctionasgiveninthelastsection,itfollows that 0 * x m y n z p = 0 * x m y n z p " 1 2 0 & x m y n z p . (3.26) Now,weinvokeanimportantstepbyequatingthediscretecentralmomentsoftheattractorsof thedistributionfunctionandthesourcetermswiththeircorrespondingcontinuouscentralmomentsderivedinSec.3.3.2.Thus,wehave 0 * at x m y n z p = 0 $ at x m y n z p , (3.27) 0 & x m y n z p = 0 % F x m y n z p . (3.28) Inotherwords,thediscretelocalcentralmomentsattractorsandthediscretecentralmomentsof sourcetermsbecome,respectively, 0 * at 0 =0 , 0 * at x = 0 * at y = 0 * at z =0 , 0 * at xx = 0 * at yy = 0 * at zz = c 2 s # , 0 * at xy = 0 * at xz = 0 * at yz =0 , 0 * at xyy = 0 * at xzz = 0 * at xxy = 0 * at yzz = 0 * at xxz = 0 * at yyz = 0 * at xyz =0 , 0 * at xxyy = ! 0 * xx ! 0 * yy , 0 * at xxzz = ! 0 * xx ! 0 * zz , 0 * at yyzz = ! 0 * yy ! 0 * zz . (3.29) 0 & 0 =0 , 0 & x = F x , 0 & y = F y , 0 & z = F z , 0 & x m y n z p =0 ,m,n,p> 1 . (3.30) BasedonEq.(3.26),wegetthefollowingattractorsintermsofthetransformedcentralmoments

PAGE 86

71 as 0 * at 0 =0 , 0 * at x = " 1 2 F x , 0 * at y = " 1 2 F y , 0 * at z = " 1 2 F z , 0 * at xx = 0 * at yy = 0 * at zz = c 2 s # , 0 * at xy = 0 * at xz = 0 * at yz =0 , 0 * at xyy = 0 * at xzz = 0 * at xxy = 0 * at yzz = 0 * at xxz = 0 * at yyz = 0 * at xyz =0 , 0 * at xxyy = ! 0 * xx ! 0 * yy , 0 * at xxzz = ! 0 * xx ! 0 * zz , 0 * at yyzz = ! 0 * yy ! 0 * zz . (3.31) 3.3.5 VariousDiscreteRawMomentsandSourceTermsinParticleVelocitySpace Inadditiontocentralmoments,thecorrespondingmomentsintherestorlatticeframeofreference,i.e.,rawmomentsareneededforobtaininganexecutablecascadedLBM.Thetoolthatwe useforthispurposeisthebinomialtheorem.Inthisregard,wewillusethefollowingnotations tospecifyvarious discreteraw moments: 0 * ! x m y n z p = " ! f ! e m ! x e n ! y e p ! z = % e m ! x e n ! y e p ! z | f ! $ , (3.32) 0 * ! x m y n z p = " ! f ! e m ! x e n ! y e p ! z = % e m ! x e n ! y e p ! z | f ! $ , (3.33) 0 & ! x m y n z p = " ! S ! e m ! x e n ! y e p ! z = % e m ! x e n ! y e p ! z | S ! $ . (3.34) wherethesuperscript"prime"( % )isusedforthevariousrawmomentsinordertodistinguishthe rawmomentsfromthecentralmomentsthataredesignatedwithouttheprimes.Furthermore, similartoEq.(3.26),therelation 0 * ! x m y n = 0 * ! x m y n z p " 1 2 0 & ! x m y n z p issatisÞed.Basedontheabove, Þrst,wewritetherawmomentsofthedistributionfunctionofdi!erentorderssupportedbythe particlevelocityset 0 * ! x m y n z p = % f ! | e m ! x e n ! y e p ! z $ intermsoftheknownquantities.

PAGE 87

72 First,theconservedtransformedrawmomentsfollowsdirectlyfromthedeÞnitionas 0 * ! 0 = % f ! | # $ = # , 0 * ! x = % f ! | e ! x $ = # u x " 1 2 F x , 0 * ! y = % f ! | e ! y $ = # u y " 1 2 F y , 0 * ! z = % f ! | e ! z $ = # u z " 1 2 F z . (3.35) Thenextstepistotransformthecentralmomentsofthesourceterms(Eq.(3.30))intermsof rawmomentsbyusingthedeÞnitions,i.e.Eq.(3.24)and(3.34),whichbythebinomialtheorem, readilyyieldstheexpressionsthatareenumeratedinAppendix2.2.Thesemomentsshouldbe relatedtothediscretesourcetermsinparticlevelocityspacesothatanoperationalGalileaninvariantapproachcanbederived.Toaccomplishthis,weÞrstobtainasetofintermediatequantities 0 m s ( ,whicharetheprojectionsofthesourcetermstotheorthogonalmatrixofthemoment basis K ,i.e. 0 m s ( = % K ( | S ! $ , + =0 , 1 , 2 ,..., 18,whichcanbeobtainedfromtheaboveusing Eqs.(3.8)and(2.-17)giveninAppendix2.2.Thedetailsof 0 m s ( areprovidedinAppendix2.3. Itisnotedthat 0 m s ( isequivalenttothefollowingmatrixformulation K   S =( Ká S ) ! =( % K 0 | S ! $ , % K 1 | S ! $ , % K 2 | S ! $ ,..., % K 18 | S ! $ ) =( 0 m s 0 , 0 m s 1 , 0 m s 2 ,..., 0 m s 18 )   & | 0 m s ! $ , (3.36) whichcanbeexactlyinvertedbyusingthefollowingorthogonalpropertyof K ,i.e. K ! 1 = K   á D ! 1 ,where D isthediagonalmatrixgivenby[22] D =diag( % K 0 | K 0 $ , % K 1 | K 1 $ , % K 2 | K 2 $ ,..., % K 18 | K 18 $ ) . (3.37) Exploitingthisfact,thelinearsystem(Eq.(3.36))canbesolvedexactlytoyieldtheexpressions fortheGalileaninvariantsourcetermsinthevelocityspace S ! intermsof 0 m s ( ,orequivalently theforce ") F andvelocityÞelds ") u .TheÞnalresultsof S ! ,where ) =0 , 1 , 2 ,..., 18aresummarizedinAppendix2.4. Finally,toobtainthechangeofdi!erentmomentsundercollision 0 g ( inthenextsection,weneed toevaluatetheexpressionsforitsrawmomentsofvariousordersprojectedtothemomentbasis

PAGE 88

73 matrix K ,i.e., " ! ( Ká 0 g ) ! e m ! x e n ! y e p ! z = " ( % K ( | e m ! x e n ! y e p ! z $ 0 g ( . (3.38) SincethezerothandÞrstordermomentsareconservedundercollision(collisioninvariants),it followsthat 0 g 0 = 0 g 1 = 0 g 2 = 0 g 3 =0.Again,exploitingtheorthogonalpropertyofthematrix K , thenon-conservedchangeofmomentsundercollision( 0 g ( )athigherorders( + =5 , 6 ,... 18)canbe obtainedwhicharedetailedinAppendix2.5. ThecascadedLBEgiveninEq.(3.23)cannowberewrittenintheformoftheusualcollisionand streamingsteps,respectively,as ! f ! ( ") x,t )= f ! ( ") x,t )+ # c ! ( !$ x,t ) + S ! ( !$ x,t ) , (3.39) f ! ( ") x + ") e ! ,t +1)= ! f ! ( ") x,t ) , (3.40) wherethesymbol"tilde"( + )intheÞrstequationreferstothepost-collisionstateofthedistributionfunction.Furthermore,weÞnallygettheßuiddensityandtheßuidvelocitybytakingthe zerothandÞrstmomentsofthedistributionfunctionsas # = 26 " ! =0 f ! = % f ! | # $ , (3.41) # u i = 26 " ! =0 f ! e ! i + 1 2 F i = % f ! | e ! i $ + 1 2 F i ,i # x,y,z. (3.42) 3.3.6 StructureoftheCascadedCollisionOperator Theexpressionsforthechangeofdi!erentmomentsundercollision 0 g ofthe3Dcentralmoment collisionoperator# c ! & # c ! ( f , 0 g )canbeobtainedasfollows.ForthezerothandÞrstordermomentsconservedundercollision, 0 g 0 = 0 g 1 = 0 g 2 = 0 g 3 =0.Inessence,theprocedurestartsfrom thelowest-ordercentralmomentsthatarenon-collisionalinvariants(i.e., 0 * xy andhigher),and theirpost-collisioncentralmoments(i.e., ! 0 * xy , ! 0 * xz and ! 0 * yz )aresuccessivelysetequaltothecorrespondingattractorsgiveninEq.(3.31)(i.e., 0 * at xy , 0 * at xz and 0 * at yz ,respectively).

PAGE 89

74 Thisintermediatestepcanprovidetentativeexpressionsfor 0 g basedonanequilibriumassumption,whichismodiÞedtoallowforcollisionasarelaxationprocess.Theyaremultipliedwith acorrespondingrelaxationparameterthatresultsintheÞnalexpressionsforthechangeofmomentsundercollision 0 g ! foragivenorder[20].Inthisstep,therelaxationparameterneedsto multiplyonlywiththosetermsthatarenotyetinpost-collisionstates(i.e.,thosenotinvolving 0 g ( ,where + =0 , 1 , 2 ,...,) " 1)foragiven 0 g ! .Thentheresultscanbetransformedinterms ofrawmomentsbyusingthebinomialtheoremtoyieldexpressionsusefulforcomputations.In ordertosummarizetheÞnalexpressionsofthenon-conservedchangeofdi!erentmomentsunder collision,wedeÞne 0 0 ! x m y n z p = 0 * ! x m y n z p + 0 & ! x m y n z p , (3.43) wheretheexpressionsfor 0 * ! x m y n z p and 0 & ! x m y n z p aredeÞnedinSec.3.3.5.Inthefollowing,for brevity,wepresentonlytheÞnalresults.Fortheabovethreeo!-diagonalcentralmoments,we get 0 g 4 = , 4 4 & " 0 0 ! xy + # u x u y + 1 2 ( 0 & ! x u y + 0 & ! y u x ) ' , (3.44) 0 g 5 = , 5 4 & " 0 0 ! xz + # u x u z + 1 2 ( 0 & ! x u z + 0 & ! z u x ) ' , (3.45) 0 g 6 = , 6 12 & " 0 0 ! yz + # u y u z + 1 2 ( 0 & ! y u z + 0 & ! z u y ) ' . (3.46) 0 g 7 = , 7 12 3 " ( 0 0 ! xx " 0 0 ! yy )+ # ( u 2 x " u 2 y )+( 0 & ! x u x " 0 & ! y u y ) 4 , (3.47) 0 g 8 = , 8 36 3 " ( 0 0 ! xx + 0 0 ! yy " 2 0 0 ! zz )+ # ( u 2 x + u 2 y " 2 u 2 z ) +( 0 & ! x u x + 0 & ! y u y " 2 0 & ! z u z ) 4 , (3.48) 0 g 9 = , 9 126 3 " ( 0 0 ! xx + 0 0 ! yy + 0 0 ! zz )+ # ( u 2 x + u 2 y + u 2 z ) +( 0 & ! x u x + 0 & ! y u y + 0 & ! z u z )+ # 4 . (3.49) 0 g 10 = , 10 8 3 " ( 0 0 ! xyy + 0 0 ! xzz )+2( u y 0 0 ! xy + u z 0 0 ! xz )+ u x ( 0 0 ! yy + 0 0 ! zz ) " 2 # u x ( u 2 y + u 2 z ) " 1 2 0 & ! x ( u 2 y + u 2 z ) " u x ( 0 & ! y u y + 0 & ! z u z ) '

PAGE 90

75 +( u y 0 g 4 + u z 0 g 5 )+ 1 4 u x ( " 3( 0 g 7 + 0 g 8 )+42 0 g 9 ) , (3.50) 0 g 11 = , 11 24 3 " ( 0 0 ! xxy + 0 0 ! yzz )+2( u x 0 0 ! xy + u z 0 0 ! yz )+ u y ( 0 0 ! xx + 0 0 ! zz ) " 2 # u y ( u 2 x + u 2 z ) " 1 2 0 & ! y ( u 2 x + u 2 z ) " u y ( 0 & ! x u x + 0 & ! z u z ) ' +( u x 0 g 4 + u z 0 g 6 )+ 3 4 u y ( 0 g 7 " 0 g 8 +14 0 g 9 ) , (3.51) 0 g 12 = , 12 8 3 " ( 0 0 ! xxz + 0 0 ! yyz )+2( u x 0 0 ! xz + u y 0 0 ! yz )+ u z ( 0 0 ! xx + 0 0 ! yy ) " 2 # u z ( u 2 x + u 2 y ) " 1 2 0 & ! z ( u 2 x + u 2 y ) " u z ( 0 & ! x u x + 0 & ! y u y ) ' +( u x 0 g 5 + u y 0 g 6 )+ 1 2 u z (3 0 g 8 +21 0 g 9 ) , (3.52) 0 g 13 = , 13 8 3 " ( 0 0 ! xyy " 0 0 ! xzz )+2( u y 0 0 ! xy " u z 0 0 ! xz )+ u x ( 0 0 ! yy " 0 0 ! zz ) " 2 # u x ( u 2 y " u 2 z ) " 1 2 0 & ! x ( u 2 y " u 2 z ) " u x ( 0 & ! y u y " 0 & ! z u z ) ' +( u y 0 g 4 " u z 0 g 5 )+ 3 4 u x ( " 0 g 7 +3 0 g 8 ) , (3.53) 0 g 14 = , 14 8 3 " ( 0 0 ! xxy " 0 0 ! yzz )+2( u x 0 0 ! xy " u z 0 0 ! yz )+ u y ( 0 0 ! xx " 0 0 ! zz ) " 2 # u y ( u 2 x " u 2 z ) " 1 2 0 & ! y ( u 2 x " u 2 z ) " u y ( 0 & ! x u x " 0 & ! z u z ) ' +3( u x 0 g 4 " u z 0 g 6 )+ 3 4 u y ( 0 g 7 +3 0 g 8 ) , (3.54) 0 g 15 = , 15 8 3 " ( 0 0 ! xxz " 0 0 ! yyz )+2( u x 0 0 ! xz " u y 0 0 ! yz )+ u z ( 0 0 ! xx " 0 0 ! yy ) " 2 # u z ( u 2 x " u 2 y ) " 1 2 0 & ! z ( u 2 x " u 2 y ) " u z ( 0 & ! x u x " 0 & ! y u y ) ' +( u x 0 g 5 + u y 0 g 6 )+ u z $ 3 2 0 g 8 + 21 2 0 g 9 % , (3.55) 0 g 16 = , 16 12 3 " ( 0 0 ! xxyy + 0 0 ! xxzz + 0 0 ! yyzz )+2 u y ( 0 0 ! xxy + 0 0 ! yzz )+2 u z ( 0 0 ! xxz + 0 0 ! yyz ) " u 2 x ( 0 0 ! yy + 0 0 ! zz )) " u 2 y ( 0 0 ! zz + 0 0 ! xx ) " u 2 z ( 0 0 ! xx + 0 0 ! yy ) " 4( u x u y 0 0 ! xy + u x u z 0 0 ! xz + u y u z 0 0 ! yz ) +3 # ( u 2 x u 2 y + u 2 x u 2 z + u 2 y u 2 z )+( ! 0 * xx ! 0 * yy + ! 0 * xx ! 0 * zz + ! 0 * yy ! 0 * zz )+ u 2 x ( u y 0 & ! y + u z 0 & ! z )+ u 2 y ( u x 0 & ! x + u z 0 & ! z )

PAGE 91

76 + u 2 z ( u x 0 & ! x + u y 0 & ! y ) 4 " 4 3 u x u y 0 g 4 " 4 3 u x u z 0 g 5 " 4 3 u y u z 0 g 6 + 1 2 ( u 2 x " u 2 y ) 0 g 7 + 1 2 ( u 2 x + u 2 y " 2 u 2 z ) 0 g 8 +( " 7( u 2 x + u 2 y + u 2 z ) " 8) 0 g 9 + 4 3 u x 0 g 10 + 4 3 u y 0 g 11 + 4 3 u z 0 g 12 , (3.56) 0 g 17 = , 17 12 3 " ( 0 0 ! xxyy + 0 0 ! xxzz + 0 0 ! yyzz )+2 A u x ( 0 0 ! xyy + 0 0 ! xzz )+2 u y ( 0 0 ! yzz + 0 0 ! xxy )+ 2 u z ( 0 0 ! xxz + 0 0 ! yyz ) B " u 2 x ( 0 0 ! yy + 0 0 ! zz ) " u 2 y ( 0 0 ! xx + 0 0 ! zz ) " u 2 z ( 0 0 ! xx + 0 0 ! yy ) " 4( u x u y 0 0 ! xy + u x u z 0 0 ! xz + u y u z 0 0 ! yz )+( ! 0 * xx ! 0 * yy + ! 0 * xx ! 0 * zz + ! 0 * yy ! 0 * zz ) +3 # ( u 2 x u 2 y + u 2 x u 2 z + u 2 y u 2 z )+ u 2 x ( u y 0 & ! y + u z 0 & ! z )+ u 2 y ( u x 0 & ! x + u z 0 & ! z ) + u 2 z ( u x 0 & ! x + 4 3 u y 0 & ! y ) ' " 4 3 u x u y 0 g 4 " 4 3 u x u z 0 g 5 " 4 3 u y u z 0 g 6 + 1 2 ( u 2 x " u 2 y ) 0 g 7 + 1 2 ( u 2 x + u 2 y " 2 u 2 z ) 0 g 8 +( " 7( u 2 x + u 2 y + u 2 z ) " 8) 0 g 9 + 4 3 u x 0 g 10 + 4 3 u y 0 g 11 + 4 3 u z 0 g 12 , (3.57) Here, 0 * at xxyy = ! 0 * xx ! 0 * yy , 0 * at xxzz = ! 0 * xx ! 0 * zz ,and 0 * at yyzz = ! 0 * yy ! 0 * zz (seeEq.(3.31)).Thisyieldsthe correspondingchangeofmomentsundercollisionas 0 g 18 = , 18 24 3 " ( 0 0 ! xxyy + 0 0 ! xxzz " 2 0 0 ! yyzz )+2 A u x 0 0 ! xyy + u x 0 0 ! xzz + u y 0 0 ! xxy + u z 0 0 ! xxz " 2( u y 0 0 ! yzz + u z 0 0 ! yyz ) B " u 2 x 0 0 ! yy " u 2 x 0 0 ! zz " u 2 y 0 0 ! xx " u 2 z 0 0 ! xx +2 u 2 y 0 0 ! zz + +2 u 2 z 0 0 ! yy " 4( u x u y 0 0 ! xy + u x u z 0 0 ! xz " 2 u y u z 0 0 ! yz )+( ! 0 * xx ! 0 * yy + ! 0 * xx ! 0 * zz " 2 ! 0 * yy ! 0 * zz ) +3 # ( u 2 x u 2 y + u 2 x u 2 z " 2 u 2 y u 2 z )+ u 2 x u y 0 & ! y + u 2 x u z 0 & ! z + u 2 y u x 0 & ! x + u 2 z u x 0 & ! x " 2 u 2 y u z 0 & ! z " 2 u 2 z u y 0 & ! y ) 4 " 2 3 u x u y 0 g 4 " 2 3 u x u z 0 g 5 + 4 3 u y u z 0 g 6 +( 1 4 ( u 2 x " u 2 y " 3 u 2 z " 2) 0 g 7 + 1 4 ( u 2 x " 5 u 2 y + u 2 z " 2) 0 g 8 + 7 4 ( " 2 u 2 x + u 2 y + u 2 z ) 0 g 9 + 2 3 u x 0 g 10 " 1 3 u y 0 g 11 " 1 3 u z 0 g 12 " u y 0 g 14 + u z 0 g 15 , (3.58) Forcalculating 0 g 16 through 0 g 18 intheaboveequations,weneedthepostcollisionstates ! 0 * xx , ! 0 * yy and ! 0 * zz .ThesecanbeobtainedfromEq.(3.26)asfollows. ! 0 * xx = ! 0 * xx + 1 2 0 & xx ,

PAGE 92

77 ! 0 * yy = ! 0 * yy + 1 2 0 & yy , ! 0 * zz = ! 0 * zz + 1 2 0 & zz , wherethesecond-ordertransformedcentralmoments,inturn,canberelatedtothecorrespondingrawmoments,whichareknown,as ! 0 * xx = ! 0 * ! xx " # u 2 x " F x u x , ! 0 * yy = ! 0 * ! yy " # u 2 y " F y u y , ! 0 * zz = ! 0 * ! zz " # u 2 z " F z u z . Notethatthesecanalsobewrittenintermsof 0 0 ! x m y n z p as ! 0 * xx = 3 0 0 ! xx +6 0 g 7 +6 0 g 8 +42 0 g 9 4 " # u 2 x " F x u x , ! 0 * yy = 3 0 0 ! yy " 6 0 g 7 +6 0 g 8 +42 0 g 9 4 " # u 2 y " F y u y , ! 0 * zz = 3 0 0 ! xx " 12 0 g 8 +42 0 g 9 4 " # u 2 z " F z u z . where , 4 , , 5 and , 1 6arerelaxationparameters.Similartothe2DcentralmomentLBMwith sourceterms[22],wecanapplytheChapman-Enskogexpansiontotheabove3Dformulationto showthatitsemergentdynamicscorrespondstotheNavier-Stokesequationsrepresentingßuid motioninthepresenceofgeneralforceÞelds.Someoftherelaxationparametersinthecollisionmodelcanberelatedtothetransportcoe"cients.Forexample,thosecorrespondingtothe second-ordermomentscontrolshearviscosity . oftheßuid.Thatis, . = c 2 s < 1 ) " " 1 2 = where , = , j where j =4 , 5 , 6 , 7 , 8.Therestoftheparameterscanbeseteitherto1(i.e.,equilibration)oradjustedindependentlytocarefullycontrolandimprovenumericalstabilitybymeans ofalinearstabilityanalysis,whileallsatisfyingtheusualbounds0 < , ( < 2.Inthiswork, forsimlicity,wesettherelaxationparametersforhigherordermomentstounity,i.e., , j =1, j =9 , 10 , 11 ,..., 18. Then,forthesimulationofnon-Newtonianßowsofpowerlawßuidsusingtheabove3Dcascaded LBM,wherethenonlinearconstitutiverelationprescribetheapparentviscositydependentonthe

PAGE 93

78 shearrate | ú ' ij | ,i.e., µ ( | ú ' ij | )= #. ( | ú ' ij | )(seeSec.3.2),therelaxationtimesofthesecondorder moments , 4 , , 5 , , 6 , , 7 and , 8 arelocallyadjustedasfollows: , 4 ( " x,t )= , 5 ( " x,t )= , 6 ( " x,t )= , 7 ( " x,t )= , 8 ( " x,t )= ( µ ( | ú ' ij | ) # c 2 s + 1 2 ) ! 1 , (3.59) Where µ ( | ú ' ij | )canbeobtainedfromEq.(3.6),whichisparametrizedbytheconsistencycoe"cient µ p andthepowerlawindex n .Themagnitudeoftheshearrate | ú ' ij | isrelatedtothesecondinvariantofthestrainratetensor S ij via | ú ' ij | = # 2 S ij S ij = C 2( S 2 xx + S 2 yy + S 2 zz +2 S 2 xy +2 S 2 xz +2 S 2 yz ) .Thestrainratetensorcomponents S xx ,S yy ,S zz ,S xy ,S yz and S xz canbeobtainedusingeither Þrstordernon-equilibriummomentsorÞnitedi!erenceapproximations. 3.3.7 OperationalStepsoftheCentralMomentLBM Finally,byexpanding( Ká 0 g ) ! inthecollisionstepofthe3DcascadedLBM(seeEqs.20and39) andusingthesourcetermsinthevelocityspace S ( (SeeAppendix2.4),thepostcollisiondistrbutionfunction f ( ,where + =0 , 1 , 2 ,..., 18,canbeobtainedasfollows: ! f 0 = f 0 +[ 0 g 0 " 30 0 g 9 +12 0 g 16 ]+ S 0 , ! f 1 = f 1 +[ 0 g 0 + 0 g 1 + 0 g 7 + 0 g 8 " 11 0 g 9 " 4 0 g 10 " 4 0 g 16 " 4 0 g 17 ]+ S 1 , ! f 2 = f 2 +[ 0 g 0 " 0 g 1 + 0 g 7 + 0 g 8 " 11 0 g 9 +4 0 g 10 " 4 0 g 16 +2 0 g 17 " 2 0 g 18 ]+ S 2 , ! f 3 = f 3 +[ 0 g 0 + 0 g 2 " 0 g 7 + 0 g 8 " 11 0 g 9 +4 0 g 11 " 4 0 g 16 +2 0 g 17 " 2 0 g 18 ]+ S 3 , ! f 4 = f 4 +[ 0 g 0 " 0 g 2 " 0 g 7 + 0 g 8 " 11 0 g 9 +4 0 g 11 " 4 0 g 16 +2 0 g 17 +2 0 g 18 ]+ S 4 , ! f 5 = f 5 +[ 0 g 0 + 0 g 3 " 2 0 g 8 " 11 0 g 9 " 4 0 g 12 " 4 0 g 16 +2 0 g 17 +2 0 g 18 ]+ S 5 , ! f 6 = f 6 +[ 0 g 0 " 0 g 3 " 2 0 g 8 " 11 0 g 9 +4 0 g 12 " 4 0 g 16 +2 0 g 17 +2 0 g 18 ]+ S 6 , ! f 7 = f 7 +[ 0 g 0 + 0 g 1 + 0 g 2 + 0 g 4 +2 0 g 8 +8 0 g 9 + 0 g 10 + 0 g 11 + 0 g 13 " 0 g 14 + 0 g 16 + 0 g 17 + 0 g 18 ]+ S 7 , ! f 8 = f 8 +[ 0 g 0 " 0 g 1 + 0 g 2 " 0 g 4 +2 0 g 8 +8 0 g 9 " 0 g 10 + 0 g 11 " 0 g 13 " 0 g 14 + 0 g 16 + 0 g 17 + 0 g 18 ]+ S 8 ,

PAGE 94

79 ! f 9 = f 9 +[ 0 g 0 + 0 g 1 " 0 g 2 " 0 g 4 +2 0 g 8 +8 0 g 9 + 0 g 10 " 0 g 11 + 0 g 13 + 0 g 14 + 0 g 16 + 0 g 17 + 0 g 18 ]+ S 9 , ! f 10 = f 10 +[ 0 g 0 " 0 g 1 " 0 g 2 + 0 g 4 +2 0 g 8 +8 0 g 9 " 0 g 10 " 0 g 11 " 0 g 13 + 0 g 14 + 0 g 16 + 0 g 17 + 0 g 18 ]+ S 10 , ! f 11 = f 11 +[ 0 g 0 + 0 g 1 + 0 g 3 + 0 g 5 + 0 g 7 " 0 g 8 +8 0 g 9 + 0 g 10 + 0 g 12 " 0 g 13 + 0 g 15 + 0 g 16 + 0 g 17 " 0 g 18 ]+ S 11 , ! f 12 = f 12 +[ 0 g 0 " 0 g 1 + 0 g 3 " 0 g 5 + 0 g 7 " 0 g 8 +8 0 g 9 " 0 g 10 + 0 g 12 + 0 g 13 + 0 g 15 + 0 g 16 + 0 g 17 " 0 g 18 ]+ S 12 , ! f 13 = f 13 +[ 0 g 0 + 0 g 1 " 0 g 3 " 0 g 5 + 0 g 7 " 0 g 8 +8 0 g 9 + 0 g 10 " 0 g 12 " 0 g 13 " 0 g 15 + 0 g 16 + 0 g 17 " 0 g 18 ]+ S 13 , ! f 14 = f 14 +[ 0 g 0 " 0 g 1 + 0 g 3 + 0 g 6 " 0 g 7 " 0 g 8 +8 0 g 9 " 0 g 10 " 0 g 12 + 0 g 13 " 0 g 15 + 0 g 16 + 0 g 17 " 0 g 18 ]+ S 14 , ! f 15 = f 15 +[ 0 g 0 + 0 g 2 + 0 g 3 + 0 g 6 " 0 g 7 " 0 g 8 +8 0 g 9 + 0 g 11 + 0 g 12 + 0 g 14 " 0 g 15 + 0 g 16 " 2 0 g 17 ]+ S 15 , ! f 16 = f 16 +[ 0 g 0 + 0 g 2 + 0 g 3 + 0 g 6 " 0 g 7 " 0 g 8 +8 0 g 9 + 0 g 11 + 0 g 12 + 0 g 14 " 0 g 15 + 0 g 16 " 2 0 g 17 ]+ S 16 , ! f 17 = f 17 +[ 0 g 0 + 0 g 2 + 0 g 3 + 0 g 6 " 0 g 7 " 0 g 8 +8 0 g 9 + 0 g 11 + 0 g 12 + 0 g 14 " 0 g 15 + 0 g 16 " 2 0 g 17 ]+ S 17 , ! f 18 = f 18 +[ 0 g 0 + 0 g 2 + 0 g 3 + 0 g 6 " 0 g 7 " 0 g 8 +8 0 g 9 + 0 g 11 + 0 g 12 + 0 g 14 " 0 g 15 + 0 g 16 " 2 0 g 17 ]+ S 18 , (3.41) Theabovepost-collisionstateallowsthecompletionofthestreamingstepviaEq.(3.40),followingwhichthehydrodynamicÞeldsofthenon-Newtonian3Dßuidmotioncanbeobtainedfrom Eqs.(3.41)and(3.42). 3.4 Resultsanddiscussion 3.4.1 Validationofthe3DcascadedLBM Tovalidatethepresented3DcascadedLBM,weconsiderthefollowingbenchmarkproblems forwhicheithertheanalyticalsolutionsand/ornumericalresultsavailable:(i)non-Newtonian Poiseuilleßowofapowerlawßuid,(ii)three-dimensionalfullydevelopedßowinasquareduct, and(iii)power-lawßuidßowinacubiclid-drivencavity(seeFig.3.3).FortheÞrsttwoproblems, thereexistanalyticalsolutions,whileforthethirdproblem,weusetherecentnumericalresults of[70]forcomparison.

PAGE 95

80 3.4.1.1 poiseuilleßowofpower-lawßuid First,weconsidertheßowofpower-lawßuidsbetweentwoparallelplatesseparatedbyaheight H drivenaconstantbodyforce.Periodicboundaryconditionsareusedattheinletandoutlets inthestreamwisedirectionandalongthespanwisedirectionandahalf-waybounce-backboundaryconditionisappliedontheÞxedwalls.Theßowisdrivenbyabodyforceof F x =1 . 0 ! 10 ! 6 , whichrepresentsthepressuregradient " * P * x .ThecharacteristicdimensionlessReynoldsnumber ReforthisproblemmaybewrittenasRe= # ( H 2 ) n u 2 ! n max /µ p ,where u max isthemaximumvelocityoftheßowoccurringmidwaybetweenthetwoplates(seebelow).Weconsideragridresolutionof3 ! 3 ! 156forresolvingthedomainoftheßowofpowerlawßuidswithRe=100and n =0 . 8 , 1 . 0and1 . 5,therebyencompassingbothshearthinningandshearthickeningßuids.The analyticalsolutionforthisnon-Newtonianßowproblemisgivenby u ( z )= n n +1 $ " 1 µ p % P % x % 1 /n D $ H 2 % ( n +1) /n " | z | ( n +1) /n E , (3.41) wherethemaximumßuidvelocityisobtainedusing u max = n n +1 $ " 1 µ p % P % x % 1 /n $ H 2 % ( n +1) /n . (3.41) Figure3.4presentsacomparisonbetweenthecomputedvelocityproÞlesfordi!erentvaluesof thepowerlawindexobtainedusingthe3DcentralmomentLBMimplementedonaD3Q19latticewiththeanalyticalsolution.Itisclearthatthenumericalsimulationfordi!erentpowerlaw ßuidsareinverygoodagreementwiththeanalyticalsolution. Inaddition,Fig.3.2demonstratesthatour3DcentralmomentLBMissecondorderaccurate forthesimulationofpowerlawßuids.Inthisgridconvergencetest,wehaveuseddi!usivescalingtoestablishasymptoticconvergenceofourschemetotheincompressiblemacroscopicnonNewtonianßuidßowequationsofshearthinning,Newtonianandshearthickeningßuids.

PAGE 96

81 3.4.1.2 Three-dimensionalfullydevelopedßowinasquareduct Theßuidßowthrougha3Dsquareduct(seeFig.3.3b)isconsideredasthenextexample.While analyticalsolutionforsuchaproblemdoesnotexistforpowerlawßuidatallvaluesofthepower lawindex,anexactsolutionfortheNewtonianßuidßowcase( n =1)exists(seebelow),which couldbeusedtotestour3DcascadedLBformulationonaD3Q19lattice.Thecrosssectionof thesquareductisdeÞnedby " a ( y ( a and " a ( z ( a ,where a isthehalfwidthoftheduct, and x beingthestreamwiseßowdirection.Periodicboundaryconditionsareappliedattheinlet andoutlets,whereasano-slipboundaryconditionisadoptedonthewallboundariesbyusingthe standardhalf-waybouncebackapproach.TheReynoldsnumberReforthisproblemisdeÞned basedonthemaximumßuidvelocityandtheducthalfwidth.Theßowisdrivenbyapplyinga bodyforceof F x =1 . 0 ! 10 ! 6 ,andagridresolutionof3 ! 45 ! 45isemployedtorepresentthe ßowdomain.TheanalyticalsolutionforthevelocityÞeldthefullydevelopedsquareductßow canbeexpressedintermsofthefollowingorthogonalFourierseries:[72] u ( y,z )= 16 a 2 F x #./ 3 " " n =1 ( " 1) ( n ! 1) * , 1 " cosh A (2 n ! 1) $ z 2 a B cosh A (2 n ! 1) $ 2 B / cos A (2 n ! 1) $ y 2 a B (2 n " 1) 3 , (3.41) Figures3.5a3.5bshowthecomparisonofthesurfacecontoursofthevelocityÞeldcomputedusingthe3DcentralmomentLBMusingaD3Q19latticeforaReynoldsnumberof80andthe analyticalsolution(seeEq.(3.4.1.2)).Itisclearthatthecomputedresultsareingoodagreementwiththeexactsolution,whichshowsaparaboloiddistributionforthevelocityproÞle.In ordertogetabetterperspectiveforcomparison,Figs.3.6And3.8showthecomparisonsofthe computedvelocityproÞlesatdi!erentlocationsintheductforRe=20and80,respectively, whichshownexcellentagreementofthe3DcentralmomentLBMwiththeanalyticalsolution.In addition,theiso-speedcontoursofthevelocityÞeldatthesetwoRearepresentedinFigs.3.7a And3.9a,respectively,whichareconsistentwithexpectations.

PAGE 97

82 3.4.1.3 Power-lawßuidßowinalid-drivencubiccavity The3Dlid-drivenßowinacubiccavityisaclassicproblempopulartovalidatenumericalmethods.ThegeometricconÞgurationofthisproblemisshowninFig.3.3a,whereapower-lawßuid isenclosedwithinacubiccavityofsidelength H .Theuppersurfaceofthecavityorthelidis subjectedtoauniformvelocity U p ,whilealltheremainingwallsofthecavityaremaintainedstationary.Asnotedintheintroduction,onlyveryfewstudiesexistonthe3Dnon-Newtonianßows ofpowerlawßuidsincubiccavitiesdrivenbyitstoplid.Arecentworkreportedby[70]that presentshighqualityresultsbasedonÞneenoughgridresolutionsusingaÞnitevolumemethod willbeusedasabenchmarksolutionforcomparison.Weemployagridresolutionof128 ! 128 ! inoursimulations,whichwasfoundtoleadtogridindependentresultsfortheReynoldsnumber casesconsideredinthisregard.Theno-slipboundaryconditionsonthewallsareimposedusing thestandardhalf-waybouncebackscheme,withthemotionofthelidrepresentedbyapplyinga momentumcorrectiontotherespectivebouncebackcondition.ThecharacteristicReynoldsnumberforthisproblemisdeÞnedbyRe= # H n U 2 ! n p /µ p .Weconsiderthesimulationofcubiccavity ßowatthreedi!erentReynoldsnumbersofRe=100 , 400and1000.Ineachcase,weconsider threetypesofnon-Newtonianßuidswithpowerlawindex n =1 . 5 , 1 . 0and0 . 8,therebyencompassingbothshearthickeningandshearthinningcases.Inallthesecasestheßowwasfoundto bestationaryandhenceallthesimulationswererununtilsteadystatewasreached.ThenumericalresultsforthecomputedproÞlesofthevelocitycomponentsalongthehorizontalandvertical centerlinesintheplane z =0 . 5 H attheabovethreeReforthepowerlaw n =1 . 5 , 1 . 0and0 . 8 areshowninFigs.3.10,3.11,and3.12,respectively.WhilefortheNewtonianßuid,theviscosityis aconstant,thelocalviscosityvariesappreciablydependingonthestrainratetensorthatchanges duetovariationsinthevorticalßowstructureinsidethecubiccavity.Asaresult,inthecaseof shearthickeningßuid,astheßowencountershighshearzonesnearthemovinglid,itse!ective viscosityincreasesintheseregions,di!usingitsmomentum(andconverselyfortheshearthinningcase)andalteringthevelocityproÞlesintherespectivecases.Hence,theboundarylayer

PAGE 98

83 becomesthickerforßuidwith n> 1whencomparedtothecasewhere n< 1.Inaddition,is foundthatasthepowerlawindex n increases,thepeakmagnitudeofthevelocitycomponent v alongthehorizontalcenterlineincreases.TheseobservationsareconsistentwiththerecentÞndingsof[70].Overall,thecomputedresultsobtainedusingthepresent3DcentralmomentLBM areingoodquantitativeagreementwiththebenchmarkresultsof[70]. 3.5 Comparisonofnumericalstabilityofdi!erentcollisionmodels Wenowmakeadirectcomparisonofthenumericalstabilitybetweenthepopularsinglerelaxationtime(SRT)modelandthecascadedMRTmodelpresentedinthiswork,bothusingthe D3Q19lattice,forthesimulationof3Dlid-drivencavityßowsimulations.Duetoitbeingashear drivenßowandaccompaniedbygeometricsingularity,thisßowconÞgurationservesasagood benchmarkproblemtotestthestabilityofthenumericalschemes.Forboththeseapproaches, foragivenresolutionandaÞxedlidvelocity,therelaxationparameter ! forthesecondorder momentsthatcontrolsviscositywasdecreasedgraduallyuntilthecomputationbecomesunstable.Hence,theinstabilitymeansthattheglobalerrorofthevelocityÞeldbecomesexponentially large.Figure3.13showsthemaximumReynoldsnumberthatcouldbeattainedbeforethecomputationsbecomeunstableforboththecollisionmodels.Resultsarepresentedforthreedi!erent valuesofthepowerlawindex( n =0 . 8 , 1 . 0and1 . 5)encompassingbothshearthinningandshear thickeningßuidsatvarioussetofgridresolutions(48 ! 48 ! 48,64 ! 64 ! 64,96 ! 96 ! 96and 128 ! 128 ! 128).ItisclearthatthecascadedMRTcomputationsachievesigniÞcantlyhigher ReynoldsnumbersthantheSRTmodelforthesamegridresolutionandforagiventypeofßuid orapowerlawindex.Hence,the3DcascadedformulationismorestablethantheSRT-LBMto simulate3Dnon-Newtonianßows.

PAGE 99

84 3.6 SummaryandConclusions Complexßuidßowssatisfyingnonlinearconstitutiverelationshipsariseinnumerousengineering,chemicalandmaterialsprocessingapplicationsandgeophysicalsituations.Simulationsof suchßowsposenumericalchallengesespeciallyin3D,whoseadequateresolutionareassociated withhighcomputationalcosts.LatticeBoltzmannmethods(LBM)areinherentlyparallelapproachesforßowsimulations.Priore!ortsusingsuchkineticschemeshavemainlyfocusedon computingnon-Newtonianßowsin2D,andutilizinglessrobustunderlyingcollisionmodels.In thiswork,wepresentanew3DcascadedLBMforthesimulationofnon-Newtonianpowerlaw ßuidsonaD3Q19lattice.Thee!ectofcollisionsareprescribedintermsofchangesindi!erent centralmomentsusingmultiplerelaxationtimes.Inparticular,therelaxationtimesofthesecondordermomentsarerelatedtothelocalshearrateandparametrizedbythemodelcoe"cient andexponentofthepowerlawßuids.Inadditions,thee!ectoftheimpressedforcesonthenonNewtonianßuidmotionisintroducedintermsofsourcetermstothe3DcascadedLBEviaappropriatechangesinthecentralmoments.Weperformedvalidationofour3Dcentralmoment LBMforvariousbenchmarkproblemsincludingthecomplexnon-NewtonianßowatReynolds numbersof100,400and1000,andwiththepowerlawindexsetequalto0.8,1.0and1.5.Comparisonsagainstrecentnumericalbenchmarksolutionsdemonstrateverygoodaccuracy.Inaddition,our3DcascadedLBMisshowntohavesuperiornumericalstabilitywhencomparedtothe SRT-LBMforsimulationof3Dsheardrivenßowsofbothshearthinningandshearthickening ßuids.

PAGE 100

85 FIGURE3.1:Three-dimensional,nineteenparticlevelocity(D3Q19)lattice. 10 1 10 2 10 3 10 ! 3 10 ! 2 10 ! 1 Grid Resolution Relative Error n=1.0 n=1.5 n=0.8 FIGURE3.2:Relativeglobalerroragainstthevariationinthegridspacing,plottedinlog-log scales,forthesimulationofpowerlawßuidßowinachannelat n =0 . 8, n =1 . 0and1 . 5obtained usingthecascadedLBmethod.

PAGE 101

86 (a) (b) FIGURE3.3:Schematicsofthenon-Newtonianbenchmarkßowproblems.(a)ßuidßowsina lid-drivencavity,(b)3Dfullydevelopedßowinasquareduct.

PAGE 102

87 ! 50 0 50 0 0.2 0.4 0.6 0.8 1 z u(z) n=0.8 n=1 n=1.5 FIGURE3.4:comparisonofthenormalizedvelocityproÞlesofpowerlawßuidsinachannel computedusing3DCascaded-LBMandtheanalyticalsolutionfor n =0 . 8 , 1 . 0 , and1 . 5.at Re=100.

PAGE 103

88 ! 20 ! 10 0 10 20 ! 20 ! 10 0 10 20 0 0.05 0.1 y Computed z u(y,z) (a) ! 20 ! 10 0 10 20 ! 20 ! 10 0 10 20 0 0.05 0.1 y Analytical z u(y,z) (b) FIGURE3.5:Flowthroughasquareductwithsidelength2asubjectedtoaconstantbody force:ComparisonofsurfacecontoursofthevelocityÞeldforReynoldsnumberRe=20(a)computedbytheD3Q15formulationofthecentralmomentLBMwithforcingtermwith(b)analyticalsolution(seeEq.3.4.1.2).

PAGE 104

89 ! 20 ! 15 ! 10 ! 5 0 5 10 15 20 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 y= ! a/2 y= ! 3a/8 y= ! a/4 y= ! a/8 z u(y,z) FIGURE3.6:Flowthroughasquareductwithsidelength2asubjectedtoaconstantbody force:ComparisonofvelocityproÞlescomputedbytheD3Q15formulationofthecentralmoment LBMwithforcingterm(symbols)withanalyticalsolution(lines)(seeEq.(3.4.1.2))atdi!erent locationsintheductcross-sectionforReynoldsnumberRe=20 y z ! 20 ! 15 ! 10 ! 5 0 5 10 15 20 ! 20 ! 15 ! 10 ! 5 0 5 10 15 20 (a) FIGURE3.7:StreamlinecontoursoftheductvelocityÞeldforReynoldsnumberRe=20

PAGE 105

90 ! 20 ! 15 ! 10 ! 5 0 5 10 15 20 0 0.01 0.02 0.03 0.04 0.05 0.06 y= ! a/2 y= ! 3a/8 y= ! a/4 y= ! a/8 z u(y,z) FIGURE3.8:Flowthroughasquareductwithsidelength2asubjectedtoaconstantbody force:ComparisonofvelocityproÞlescomputedbytheD3Q19formulationofthecentralmoment LBMwithforcingterm(symbols)withanalyticalsolution(lines)(seeEq.(3.4.1.2))atdi!erent locationsintheductcross-sectionforReynoldsnumberRe=20 y z ! 20 ! 15 ! 10 ! 5 0 5 10 15 20 ! 20 ! 15 ! 10 ! 5 0 5 10 15 20 (a) FIGURE3.9:StreamlinecontoursoftheductvelocityÞeldforReynoldsnumberRe=80

PAGE 106

91 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.2 0 0.2 0.4 0.6 0.8 1 ! U ! profile U/U p Y/H Present (n=1.5) Jin et al (2017) (n=1.5) 0 0.2 0.4 0.6 0.8 1 ! 0.3 ! 0.2 ! 0.1 0 0.1 0.2 V ! profile " X/H V/U (a)Re=100 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.2 0 0.2 0.4 0.6 0.8 1 ! U ! profile U/U p Y/H Present (n=1.5) Jin et al (2017) (n=1.5) 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.3 ! 0.2 ! 0.1 0 0.1 0.2 0.3 V ! profile " X/H V/U (b)Re=400 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.2 0 0.2 0.4 0.6 0.8 1 ! U ! profile U/U p Y/H Present (n=1.5) Jin et al (2017) (n=1.5) 0 0.2 0.4 0.6 0.8 1 ! 0.5 ! 0.4 ! 0.3 ! 0.2 ! 0.1 0 0.1 0.2 0.3 V ! profile " X/H V/U (c)Re=1000 FIGURE3.10:Computedvelocitycomponent u/U p alongtheverticalcenterlineandcomponent v/U p alongthehorizontalcenterlineobtainedusingthecascadedLBmethod(lines)andcomparedwiththebenchmarksolutionofJin(2017)(symbols)forlid-drivencavityßowofpowerlaw ßuids.Reynoldsnumber Re =100 , 400 , 1000andthepowerlawindex n =1 . 5.

PAGE 107

92 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.2 0 0.2 0.4 0.6 0.8 1 ! U ! profile U/U p Y/H Present (n=1) Jin et al (2017) (n=1) 0 0.2 0.4 0.6 0.8 1 ! 0.25 ! 0.2 ! 0.15 ! 0.1 ! 0.05 0 0.05 0.1 0.15 V ! profile " X/H V/U (a)Re=100 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.2 0 0.2 0.4 0.6 0.8 1 ! U ! profile U/U p Y/H Present (n=1) Jin et al (2017) (n=1) 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.3 ! 0.2 ! 0.1 0 0.1 0.2 0.3 V ! profile " X/H V/U (b)Re=400 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.2 0 0.2 0.4 0.6 0.8 1 ! U ! profile U/U p Y/H Present (n=1) Jin et al (2017) (n=1) 0 0.2 0.4 0.6 0.8 1 ! 0.5 ! 0.4 ! 0.3 ! 0.2 ! 0.1 0 0.1 0.2 0.3 V ! profile " X/H V/U (c)Re=1000 FIGURE3.11:Computedvelocitycomponent u/U p alongtheverticalcenterlineandcomponent v/U p alongthehorizontalcenterlineobtainedusingthecascadedLBmethod(lines)andcomparedwiththebenchmarksolutionofJin(2017)(symbols)forlid-drivencavityßowofpowerlaw ßuids.Reynoldsnumber Re =100 , 400 , 1000andthepowerlawindex n =1 . 0.

PAGE 108

93 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.2 0 0.2 0.4 0.6 0.8 1 1.2 ! U ! profile U/U p Y/H 0 0.2 0.4 0.6 0.8 1 ! 0.25 ! 0.2 ! 0.15 ! 0.1 ! 0.05 0 0.05 0.1 0.15 V ! profile " X/H V/U (a)Re=100 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.2 0 0.2 0.4 0.6 0.8 1 ! U ! profile U/U p Y/H 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.3 ! 0.2 ! 0.1 0 0.1 0.2 0.3 V ! profile " X/H V/U (b)Re=400 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.2 0 0.2 0.4 0.6 0.8 1 ! U ! profile U/U p Y/H 0 0.2 0.4 0.6 0.8 1 ! 0.4 ! 0.3 ! 0.2 ! 0.1 0 0.1 0.2 0.3 V ! profile " X/H V/U (c)Re=1000 FIGURE3.12:Computedvelocitycomponent u/U p alongtheverticalcenterlineandcomponent v/U p alongthehorizontalcenterlineobtainedusingthecascadedLBmethod(lines)andcomparedwiththebenchmarksolutionofJin(2017)(symbols)forlid-drivencavityßowofpowerlaw ßuids.Reynoldsnumber Re =100 , 400 , 1000andthepowerlawindex n =0 . 8.

PAGE 109

94 48x48x48 64x64x64 96x96x96 128x128x128 0 100 200 300 400 500 600 700 800 Maximum Re Power ! law index (n=0.8) Grid Resolution SRT Central Moment MRT (a)Re=100 48x48x48 64x64x64 96x96x96 128x128x128 0 500 1000 1500 2000 2500 Maximum Re Power ! law index (n=1) Grid Resolution SRT Central Moment MRT (b)Re=400 48x48x48 64x64x64 96x96x96 128x128x128 0 100 200 300 400 500 600 700 800 Maximum Re Power ! law index (n=1.5) Grid Resolution SRT Central Moment MRT (c)Re=1000 FIGURE3.13:ComparisonofthemaximumReynoldsnumberfornumericalstabilityofSRT LBscheme,andcascadedLBschemeforsimulationofnon-NewtonianLiddrivencavityßowwith n =0 . 8 , 1 . 0;1 . 5

PAGE 110

95 CHAPTERIV SUMMARY,CONCLUSIONSANDOUTLOOK Inthischapter,webrießyÞrstsummarizethemainaccomplishmentsofthisthesisandsomerecommendationsforfuturework.Fluidßowsexhibitingnonlinearrheologicalbehavioroccurinnumerousengineering,chemicalandmaterialsprocessingapplications,andgeophysicalcontexts. Numericalsimulationsofsuchßowsarechallengingandcomputationallyintensive,especially in3Dwhereveryfewbenchmarkqualitynumericalsolutionsexist.LatticeBoltzmannmethods (LBMs)arebasedonkineticmodelsandinvolvelocalstream-and-collidesteps.Priore!ortson LBMfornon-Newtonianßowshavelargelyfocusedonproblemsin2Dusinglessrobustcollision models,andtheirapplicationsto3Darerelativelyrare.Inthisresearch,wehaveformulatedadvancedcascadedLBMsin2Dand3Dandincludingsourcetermsforimpressedforcesforsimulationsofnon-Newtonianßowsofpowerlawßuids.Theyarebasedoncentralmomentsandusing multiplerelaxationtimes.Therelaxationtimesofthesecondordermomentsarerelatedtothe localshearrateandparametrizedbythemodelcoe"cientandexponentofthepowerlawßuids, andthosefortheothermomentscanbetunedtoimprovenumericalstability.Inaddition,theeffectoftheimpressedforcesonnon-Newtonianßuidmotionisprescribedintermsofsourceterms tothesecascadedLBMsbymeansofappropriatechangesinthecentralmoments,whicharethen discretizedusingatrapezoidalrule.Numericalvalidationstudiesofthenew2Dand3Dcascaded LBMswereperformedforvariousbenchmarkproblemsencompassingbothshearthinningand shearthickeningbehaviorandtheircomparisonsagainstanalytical/numericalsolutions.These includethenon-Newtonianßowofpowerlawßuidsinacubiccavitydriventhemotionofthetop lidReynoldsnumbersof100,400and1000,andwiththepowerlawindexsetequalto0.8,1.0 and1.5.Comparisonsagainstvariousbenchmarksolutionsdemonstratetheexcellentaccuracy ofthecascadedLBMsforsimulationofnon-Newtonianßows.Itwasalsodemonstratedthatthe cascadedLBMapproachexhibitssecondorderaccuracyingridconvergence.Furthermore,by

PAGE 111

96 meansofsimulationsof2Dsquarecavityand3Dcubiccavityßowsofbothshearthinningand shearthickeningßuids,thesuperiornumericalstabilitycharacteristicsofthe2Dand3Dcascaded LBMs,respectively,whencomparedtothesinglerelaxationtimeLBMwereshown.Investigationsofthepowerlawßuidßowoverasinglecylinderrevealedanon-monotonicdependenceof thedragcoe"cientontheReynoldsnumber.Finally,thee!ectsofvariouscharacteristicparametersincludingthepowerlawindexandthegapratioonthewakeinteractionsanddragforces arisingfromthepowerlawßuidßowoverapairofcylindersinatandemarrangement. Inthelightofthisresearch,therearenumberofdirectionsinwhichtheproposedmethodscan befurtherextendedandapplied.Somerecommendationsforfutureworkincludethefollowing: • ApplicationsofthenewcascadedLBMsdevelopedinthisworkforfundamentalinvestigationofvariouscomplexnon-Newtonianßowsindi!erentgeometricconÞgurations andßowconditions.Examplesincludetheexternalßowofpowerlawßuidsoverasetof cylindersorspheresandinternalßowthroughchannels,ductsandpipesathighReynolds numbers,encompassingtheturbulentßowregime. • Extensionsofthe2Dand3DcascadedLBMsforviscoelasticßowsendowedwithmemory e!ects.Thiswouldinvolvesolutionsofadditionaltransportequationsfortheviscoelastic stresstensor(e.g.,representedbytheOldroyd-Bconstitutiveequation),andtheircouplingtotheßuidmotionviacorrectionstothesecondordermomentsofthecascaded LBMsolverfortheßuidmotion. • FurtherdevelopmentofthecascadedLBMindi!erentdimensionsfornon-Newtonian ßowswithheattransferaswellastwo-phasesystems.

PAGE 112

BIBLIOGRAPHY [1] R.C.Patil,R.P.Bharti,andR.P.Chhabra,"Steadyßowofpowerlawßuidsoverapair ofcylindersintandemarrangement," Industrial & Engineering Chemistry Research,vol.47, no.5,pp.1660Ð1683,2008. [2] P.Neofytoul,"A3rdorderupwindÞnitevolumemethodforgeneralisedNewtonianßuid ßows," Int. J. Numer Meth Fluids,vol.36,pp.644Ð680,2005. [3] R.P.Bharti,R.Chhabra,andV.Eswaran,"Steadyßowofpowerlawßuidsacrossacircularcylinder," The Canadian journal of chemical engineering,vol.84,no.4,pp.406Ð421, 2006. [4] R.P.Chhabra, Bubbles, drops, and particles in non-Newtonian ßuids.CRCpress,2006. [5] H.A.Barnes,J.F.Hutton,andK.Walters, An introduction to rheology.Elsevier,1989. [6] R.B.Bird,"Transportphenomena," Applied Mechanics Reviews,vol.55,no.1,pp.R1ÐR4, 2002. [7] T.Kr¬uger et al.,"ThelatticeBoltzmannmethod:Principlesandpractice.sl:Springerinternationalpublishing,"2017. [8] J.Harris, Rheology and non-Newtonian ßow.LongmanPublishingGroup,1977. [9] M.DevilleandT.B.Gatski, Mathematical modeling for complex ßuids and ßows.Springer Science&BusinessMedia,2012. [10] R.G.OwensandT.N.Phillips, Computational rheology,vol.14.WorldScientiÞc,2002. [11] R.Chhabra, Bubbles, Drops, and Particles in Non-Newtonian Fluids.CRCpress,BocaRaton,2007. [12] S.ChenandG.D.Doolen,"LatticeBoltzmannmethodforßuidßows," Annual review of ßuid mechanics,vol.30,no.1,pp.329Ð364,1998. [13] C.AidunandJ.Clausen,"Lattice-Boltzmannmethodforcomplexßows," Annu. Rev. Fluid Mech.,vol.42,pp.439Ð472,2010. [14] L.Luo,M.Krafczyk,andW.Shyy,"latticeBoltzmannmethodsforcomputationalßuiddymamics,"vol.56,pp.651Ð660,2010. [15] T.N.PhillipsandG.W.Roberts,"LatticeBoltzmannmodelsfornon-newtonianßows," IMA journal of applied mathematics,vol.76,no.5,pp.790Ð816,2011.

PAGE 113

98 [16] S.SucciandS.Succi, The lattice Boltzmann equation: for ßuid dynamics and beyond.Oxforduniversitypress,2001. [17] Z.GuoiandC.Shu, Lattice Boltzmann Method and Its Applications in Engineering,vol.3. WorldScientiÞc,Singapore,2013. [18] Y.Qian,D.d'Humi`eres,andP.Lallemand,"LatticeBGKmodelsforNavier-Stokesequation," EPL (Europhysics Letters),vol.17,no.6,p.479,1992. [19] D.d'Humi`eres,"MultipleÐrelaxationÐtimelatticeBoltzmannmodelsinthreedimensions," Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences,vol.360,no.1792,pp.437Ð451,2002. [20] M.Geier,A.Greiner,andJ.G.Korvink,"CascadeddigitallatticeBoltzmannautomatafor highreynoldsnumberßow," Physical Review E,vol.73,no.6,p.066705,2006. [21] P.Asinari,"GeneralizedlocalequilibriuminthecascadedlatticeBoltzmannmethod," Phys. Rev. E,vol.78,p.016701,2008. [22] K.N.PremnathandS.Banerjee,"IncorporatingforcingtermsincascadedlatticeBoltzmannapproachbymethodofcentralmoments," Physical Review E,vol.80,no.3, p.036702,2009. [23] K.N.PremnathandS.Banerjee,"Onthethree-dimensionalcentralmomentlatticeBoltzmannmethod," Journal of Statistical Physics,vol.143,no.4,pp.747Ð794,2011. [24] M.Geier,M.Schonherr,A.Pasquali,andM.Krafczyk,"ThecumulantlatticeBoltzmannequationinthreedimensions:Theoryandvalidation," Comp. Math. Appl.,vol.704, pp.507Ð547.,2015. [25] Y.Ning,K.Premnath,andD.Patil,"NumericalstudyofthepropertiesofthecentralmomentlatticeBoltzmannmethod," Int. J. Numer Meth Fluids,vol.82,pp.59Ð90,2015. [26] A.DeRosis,"Non-orthogonalcentralmomentsrelaxingtoadiscreteequilibrium:Ad2q9 latticeBoltzmannmodel," EPL (Europhysics Letters),vol.116,no.4,p.44003,2017. [27] F.HajabdollahiandK.N.Premnath,"ImprovingthelowmachnumbersteadystateconvergenceofthecascadedlatticeBoltzmannmethodbypreconditioning," Computers & Mathematics with Applications,2017. [28] F.M.Elseid,S.W.Welch,andK.N.Premnath,"AcascadedlatticeBoltzmannmodelfor thermalconvectiveßowswithlocalheatsources," International Journal of Heat and Fluid Flow,vol.70,pp.279Ð298,2018. [29] F.HajabdollahiandK.N.Premnath,"Centralmoments-basedcascadedlatticeBoltzmann methodforthermalconvectiveßowsinthree-dimensions," International Journal of Heat and Mass Transfer,vol.120,pp.838Ð850,2018.

PAGE 114

99 [30] F.HajabdollahiandK.N.Premnath,"Centralmoments-basedcascadedlatticeBoltzmannmethodforthermalconvectiveßowsinthree-dimensions," Int. J. Heat Mass Transfer, vol.120,pp.838Ð850,2018. [31] E.AharonovandD.H.Rothman,"Non-newtonianßow(throughporousmedia):AlatticeBoltzmannmethod," Geophysical Research Letters,vol.20,no.8,pp.679Ð682,1993. [32] I.Ginzburg,"Equilibrium-typeandlink-typelatticeBoltzmannmodelsforgenericadvection andanisotropic-dispersionequation," Advances in water resources,vol.28,no.11,pp.1171Ð 1195,2005. [33] E.S.Boek,J.Chin,andP.V.Coveney,"LatticeBoltzmannsimulationoftheßowofnonnewtonianßuidsinporousmedia," International Journal of Modern Physics B,vol.17, no.01n02,pp.99Ð102,2003. [34] D.Kehrwald,"LatticeBoltzmannsimulationofshear-thinningßuids," Journal of statistical physics,vol.121,no.1-2,pp.223Ð237,2005. [35] S.Gabbanelli,G.Drazer,andJ.Koplik,"LatticeBoltzmannmethodfornon-newtonian (power-law)ßuids," Physical review E,vol.72,no.4,p.046312,2005. [36] J.Boyd,J.Buick,andS.Green,"Asecond-orderaccuratelatticeBoltzmannnon-newtonian ßowmodel," Journal of physics A: Mathematical and General,vol.39,no.46,p.14241, 2006. [37] M.Yoshino,Y.-h.Hotta,T.Hirozane,andM.Endo,"Anumericalmethodforincompressiblenon-newtonianßuidßowsbasedonthelatticeboltzmannmethod," Journal of non-newtonian ßuid mechanics,vol.147,no.1-2,pp.69Ð78,2007. [38] A.Vikhansky,"Lattice-Boltzmannmethodforyield-stressliquids," Journal of Non-Newtonian Fluid Mechanics,vol.155,no.3,pp.95Ð100,2008. [39] O.Malaspinas,G.Courbebaisse,andM.Deville,"Simulationofgeneralizednewtonianßuids withthelatticeBoltzmannmethod," International Journal of Modern Physics C,vol.18, no.12,pp.1939Ð1949,2007. [40] S.Sullivan,L.Gladden,andM.Johns,"Simulationofpower-lawßuidßowthroughporous mediausinglatticeboltzmanntechniques," Journal of Non-Newtonian Fluid Mechanics, vol.133,no.2-3,pp.91Ð98,2006. [41] S.Sullivan,A.Sederman,M.Johns,andL.Gladden,"VeriÞcationofshear-thinninglbsimulationsincomplexgeometries," Journal of non-newtonian ßuid mechanics,vol.143,no.2-3, pp.59Ð63,2007. [42] C.-H.WangandJ.-R.Ho,"LatticeBoltzmannmodelingofBinghamplastics," Physica A: Statistical Mechanics and its Applications,vol.387,no.19-20,pp.4740Ð4748,2008.

PAGE 115

100 [43] G.Tang,S.Wang,P.Ye,andW.Tao,"Binghamßuidsimulationwiththeincompressible latticeBoltzmannmodel," Journal of Non-Newtonian Fluid Mechanics,vol.166,no.1-2, pp.145Ð151,2011. [44] M.Ohta,T.Nakamura,Y.Yoshida,andY.Matsukuma,"LatticeBoltzmannsimulations ofviscoplasticßuidßowsthroughcomplexßowchannels," Journal of Non-Newtonian Fluid Mechanics,vol.166,no.7-8,pp.404Ð412,2011. [45] Z.Chai,B.Shi,Z.Guo,andF.Rong,"Multiple-relaxation-timelatticeBoltzmannmodel forgeneralizednewtonianßuidßows," Journal of Non-Newtonian Fluid Mechanics,vol.166, no.5-6,pp.332Ð342,2011. [46] D.Conrad,A.Schneider,andM.B¬ohle,"Accuracyofnon-newtonianlatticeBoltzmannsimulations," Journal of Computational Physics,vol.301,pp.218Ð229,2015. [47] Z.YangandW.-A.Yong,"AsymptoticanalysisofthelatticeBoltzmannmethodforgeneralizednewtonianßuidßows," Multiscale Modeling & Simulation,vol.12,no.3,pp.1028Ð1045, 2014. [48] T.C.Papanastasiou,"Flowsofmaterialswithyield," Journal of Rheology,vol.31,no.5, pp.385Ð404,1987. [49] X.He,S.Chen,andR.Zhang,"AlatticeBoltzmannschemeforincompressiblemultiphaseßowanditsapplicationinsimulationofrayleighÐtaylorinstability," Journal of Computational Physics,vol.152,no.2,pp.642Ð663,1999. [50] M.Geier,J.Greiner,andF.Korvink,"CascadeddigitallatticeBoltzmannautomatafor highReynoldsnumberßow," Phys. Rev. E,vol.73,p.066705,2006. [51] P.J.Dellar,"LatticeBoltzmannalgorithmswithoutcubicdefectsingalileaninvarianceon standardlattices," Journal of Computational Physics,vol.259,pp.270Ð283,2014. [52] M.Geier,M.Sch¬onherr,A.Pasquali,andM.Krafczyk,"ThecumulantlatticeBoltzmann equationinthreedimensions:Theoryandvalidation," Computers & Mathematics with Applications,vol.70,no.4,pp.507Ð547,2015. [53] F.HajabdollahiandK.N.Premnath,"Galilean-invariantpreconditionedcentral-moment latticeBoltzmannmethodwithoutcubicvelocityerrorsfore"cientsteadyßowsimulations," Physical Review E,vol.97,no.5,p.053303,2018. [54] Y.Ning,K.N.Premnath,andD.V.Patil,"Numericalstudyofthepropertiesofthecentral momentlatticeBoltzmannmethod," International Journal for Numerical Methods in Fluids, vol.82,no.2,pp.59Ð90,2016. [55] B.C.BellandK.S.Surana,"p-versionleastsquaresÞniteelementformulationfortwodimensional,incompressible,non-newtonianisothermalandnon-isothermalßuidßow," International journal for numerical methods in ßuids,vol.18,no.2,pp.127Ð162,1994.

PAGE 116

101 [56] D.Vola,L.Boscardin,andJ.Latch«e,"LaminarunsteadyßowsofBinghamßuids:anumericalstrategyandsomebenchmarkresults," Journal of Computational Physics,vol.187, no.2,pp.441Ð456,2003. [57] F.ZinaniandS.Frey,"Galerkinleast-squaresmultiÞeldapproximationsforßowsofinelastic non-newtonianßuids," Journal of Fluids Engineering,vol.130,no.8,p.081507,2008. [58] U.Ghia,K.N.Ghia,andC.Shin,"High-resolutionsforincompressibleßowusingthe Navier-Stokesequationsandamultigridmethod," Journal of computational physics,vol.48, no.3,pp.387Ð411,1982. [59] A.Soares,J.Ferreira,andR.Chhabra,"Flowandforcedconvectionheattransferincrossßowofnon-newtonianßuidsoveracircularcylinder," Industrial & Engineering Chemistry Research,vol.44,no.15,pp.5815Ð5827,2005. [60] M.Bouzidi,M.Firdaouss,andP.Lallemand,"MomentumtransferofaBoltzmann-lattice ßuidwithboundaries," Physics of ßuids,vol.13,no.11,pp.3452Ð3459,2001. [61] R.Mei,D.Yu,W.Shyy,andL.-S.Luo,"ForceevaluationinthelatticeBoltzmannmethod involvingcurvedgeometry," Physical Review E,vol.65,no.4,p.041203,2002. [62] S.Chen,S.Dawson,G.Doolen,D.Janecky,andA.Lawniczak,"Latticemethodsandtheir applicationstoreactingsystems," Computers & chemical engineering,vol.19,no.6-7, pp.617Ð646,1995. [63] H.Chen,S.Chen,andW.H.Matthaeus,"RecoveryoftheNavier-Stokesequationsusinga lattice-gasBoltzmannmethod," Physical Review A,vol.45,no.8,p.R5339,1992. [64] P.L.Bhatnagar,E.P.Gross,andM.Krook,"Amodelforcollisionprocessesingases.i. smallamplitudeprocessesinchargedandneutralone-componentsystems," Physical review, vol.94,no.3,p.511,1954. [65] I.V.Karlin,A.Ferrante,andH.C. ¬ Ottinger,"PerfectentropyfunctionsofthelatticeBoltzmannmethod," EPL (Europhysics Letters),vol.47,no.2,p.182,1999. [66] A.DeRosis,"AlternativeformulationtoincorporateforcingtermsinalatticeBoltzmann schemewithcentralmoments," Physical Review E,vol.95,no.2,p.023311,2017. [67] L.Fei,K.H.Luo,andQ.Li,"Three-dimensionalcascadedlatticeBoltzmannmethod:simpliÞedimplementationandconsistentforcingscheme," arXiv preprint arXiv:1801.04586, 2018. [68] M.ReddyandJ.Reddy,"Finite-elementanalysisofßowsofnon-newtonianßuidsinthreedimensionalenclosures," International journal of non-linear mechanics,vol.27,no.1,pp.9Ð 26,1992.

PAGE 117

102 [69] R.N.Elias,M.A.Martins,andA.L.Coutinho,"Paralleledge-basedsolutionofviscoplastic ßowswiththesupg/pspgformulation," Computational Mechanics,vol.38,no.4-5,pp.365Ð 381,2006. [70] K.Jin,S.Vanka,R.Agarwal,andB.Thomas,"Gpuacceleratedsimulationsofthreedimensionalßowofpower-lawßuidsinadrivencube," International Journal of Computational Fluid Dynamics,vol.31,no.1,pp.36Ð56,2017. [71] M.Geier,A.Greiner,andJ.Korvink,"AfactorizedcentralmomentlatticeBoltzmann method," The European Physical Journal Special Topics,vol.171,no.1,pp.55Ð61,2009. [72] F.M.WhiteandI.CorÞeld, Viscous ßuid ßow,vol.3.McGraw-HillNewYork,2006.

PAGE 118

103 APPENDIXA 1.1 Chapman-EnskoganalysisofthecentralmomentLBM:Cubicvelocityerrors, equilibrium-momentscorrectionsandstrainratetensorformulas DuetothealiasingissuesarisingfromtheÞniteparticlevelocitysetassociatedwiththestandard D2Q9lattice,forexample,thethirdorderdiagonalmomentsdegeneratestotheÞrstordermoments:e.g., 0 * ! xxx = > ! f ! e 3 ! x = > ! f ! e ! x = 0 * ! x .Thisleadstonon-Galileaninvariantcubic velocityerrorsintheemergentmacroscopicequationsforanycollisionmodelusedwiththeLBM, andappropriatecorrectionsneedtobeprescribedtoeliminatethem.Inthisappendix,wewill performaChapman-Enskog(C-E)analysisinordertoidentifythecubicvelocityerrorsforthe 2DcentralmomentbasedLBM,andthenprovidecorrectionsviaextendedmomentequilibria, and,Þnally,deriveexpressionsforthestrainratetensorbasedonnon-equiliriummoments.They arebasedontheapproachpresentedinRef.[53]. 1.1.1 IdentiÞcationofcubicvelocityerrors. Westartwiththefollowingnon-orthogonalsetofeigenvectors T ,whichformsasabasisforthe orthogonaltransformationmatrix K usedinthecascadedcascadedcollisionoperatorgiveninthe maintext: T = ( | e ! | 0 , | e ! x $ , | e ! y $ , | e 2 ! x + e 2 ! y $ , | e 2 ! x " e 2 ! y $ , | e ! x e ! y $ , | e 2 ! x e ! y $ , | e ! x e 2 ! y $ , | e 2 ! x e 2 ! y $ ) , (1.0) Basedonthis,wethendeÞnetherelationsbetweentherawmomentsofthedistributionfunctions,theirequilibriaandthesourcetermsas 0 m = T f , 0 m eq = T f eq , 0 S = T S (1.0)

PAGE 119

104 where 0 f = A 0 f 0 , 0 f 1 , 0 f 2 ,..., 0 f 8 B   , 0 f eq = A 0 f eq 0 , 0 f eq 1 , 0 f eq 2 ,..., 0 f eq 8 B   and 0 S = A 0 S 0 , 0 S 1 , 0 S 2 ,..., 0 S 8 B   .The componentsofthevariousrawmomentscanthenbewrittenas 0 m = A 0 m 0 , 0 m 1 , 0 m 2 ,..., 0 m 8 B   = A 0 * ! 0 , 0 * ! x , 0 * ! y , 0 * ! xx + 0 * ! yy , 0 * ! xx " 0 * ! yy , 0 * ! xy , 0 * ! xxy , 0 * ! xyy , 0 * ! xxyy B   , 0 m eq =( 0 m eq 0 , 0 m eq 1 , 0 m eq 2 ,..., 0 m eq 8 )   = A 0 * eq ! 0 , 0 * eq ! x , 0 * eq ! y , 0 * eq ! xx + 0 * eq ! yy , 0 * eq ! xx " 0 * eq ! yy , 0 * eq ! xy , 0 * eq ! xxy , 0 * eq ! xyy , 0 * eq ! xxyy B   , 0 S = A 0 S 0 , 0 S 1 , 0 S 2 ,..., 0 S 8 B   = A 0 & ! 0 , 0 & ! x , 0 & ! y , 0 & ! xx + 0 & ! yy , 0 & ! xx " 0 & ! yy , 0 & ! xy , 0 & ! xxy , 0 & ! xyy , 0 & ! xxyy B   . ThecentralmomentbasedcascadedLBMcanberewrittenintermsoftherelaxationofvariousrawmomentstotheirgeneralizedequilibria,andinordertoestablishconsistencywiththe Navier-Stokesequations(NSE),itissu"cienttoconsidertermsonlyuptosecondorderinthe Machnumber[21,22].InordertofacilitatetheC-Eanalysis,weadoptthisstrategyhere.Hence, themomentrelaxationformulationofsuchaLBMcanbewrittenas ! f ( " x + " e ! $ t )= f ( " x,t )+ T ! 1 & " 0 ' ( 0 m " 0 m eq )+ $ I " 1 2 0 ' % 0 S ' , (1.-3) Wheretherawmomentequilibria 0 m eq andsourcemoments 0 S canbewrittenas 0 * eq ! 0 = # , 0 * eq ! x = # u x , 0 * eq ! y = # u y , 0 * eq ! xx = 1 3 # + # u 2 x , 0 * eq ! yy = 1 3 # + # u 2 y , 0 * eq ! xy = # u x u y , 0 * eq ! xxy = 1 3 # u y + # u 2 x u y , 0 * eq ! xyy = 1 3 # u x + # u x u 2 y , 0 * eq ! xxyy = 1 9 # + 1 3 # ( u 2 x + u 2 y )+ # u 2 x u 2 y , (1.-10) and 0 & ! 0 =0 ,

PAGE 120

105 0 & ! x = F x , 0 & ! y = F y , 0 & ! xx =2 F x u x , 0 & ! yy =2 F y u y , 0 & ! xy = F x u y + F y u x , 0 & ! xxy = F y u 2 x +2 F x u x u y , 0 & ! xyy = F x u 2 y +2 F y u y u x , 0 & ! xxyy =2 F x u x u 2 y +2 F y u y u 2 x , (1.-17) and 0 , isthediagonalrelaxationtimematrixgivenby 0 , = diag (0 , 0 , 0 , , 3 , , 4 , , 5 , , 6 , , 7 , , 8 ) Then,weperformthestandardC-Eanalysisbyexpandingthemomentsandthetimederivative atdi!erentordersofasmallbookkeepingperturbationparameter 1 as 0 m = " " n =0 1 n 0 m ( n ) , % t = " " n =0 1 n % t n , whereusually 1 = $ t .Inaddition,thestreamingoperatorinthelefthandsideofEq.(1.1.1)can beexpandedviaaTaylorexpansionas f ( " x + " e ! 1 ,t + 1 )= n " n =0 1 n n ! ( % t + " e ! á " ) n f ( " x,t ) . (1.-17) AftersubstitutingtheaboveexpansionsinEq.(1.1.1)andconvertingallthetermsintothemomentspaceviaEq.(1.1.1),weobtainthefollowingsequenceofmomentequationsatsuccessively higherordersin 1 as O ( 1 0 ): 0 m (0) = 0 m eq , (1.-16) O ( 1 1 ):( % t 0 + 0 E i % i ) 0 m (0) = " 0 ' 0 m (1) + 0 S , (1.-16) O ( 1 2 ): % t 1 0 m (0) +( % t 0 + 0 E i % i ) & I " 1 2 0 ' ' 0 m (1) = " 0 ' 0 m (2) , (1.-16)

PAGE 121

106 where 0 E i = T ( e ! i I ) T ! 1 ,i # x,y .Then,the O ( 1 )equationsuptothosefortheevolutionofthe secondordermomentsneededforthederivationoftheNSEcanbewrittenas % t 0 # + % x ( # u x )+ % y ( # u y )=0 , (1.-15) % t 0 ( # u x )+ % x $ 1 3 # + # u 2 x % + % y ( # u x u y )= F x , (1.-15) % t 0 ( # u y )+ % x ( # u x u y )+ % y $ 1 3 # + # u 2 y % = F y , (1.-15) % t 0 $ 2 3 # + # < u 2 x + u 2 y = % + % x $ 4 3 # u x + # u x u 2 y % + % y $ 4 3 # u y + # u 2 x u y % = " , 3 0 m (1) 3 +2 F x u x +2 F y u y , (1.-15) % t 0 < # < u 2 x " u 2 y == + % x $ 2 3 # u x " # u x u 2 y % + % y $ " 2 3 # u y + # u 2 x u y % = " , 4 0 m (1) 4 +2 F x u x " 2 F y u y , (1.-15) % t 0 ( # u x u y )+ % x $ 1 3 # u y + # u 2 x u y % + % y $ 1 3 # u x + # u x u 2 y % = " , 5 0 m (1) 5 + F x u y + F y u x . (1.-15) Analogously,the O ( 1 2 )equationsuptothosefortheevolutionoftheÞrstordermomentsthat willbeusedtorecovertorecovertheNSEreadas % t 1 # =0 , (1.-14) % t 1 ( # u x )+ % x & 1 2 $ 1 " 1 2 , 3 % 0 m (1) 3 + 1 2 $ 1 " 1 2 , 4 % 0 m (1) 4 ' + % y &$ 1 " 1 2 , 5 % 0 m (1) 5 ' =0 , (1.-14) % t 1 ( # u y )+ % x &$ 1 " 1 2 , 5 % 0 m (1) 5 ' + % y & 1 2 $ 1 " 1 2 , 3 % 0 m (1) 3 " 1 2 $ 1 " 1 2 , 4 % 0 m (1) 4 ' =0 , (1.-14) IntheaboveEqs.(1.-14),(1.-14),thesecondordernon-equilibriummoments 0 m (1) 3 = 0 * ! (1) xx + 0 * ! (1) yy , 0 m (1) 4 = 0 * ! (1) xx " 0 * ! (1) yy ,and 0 m (1) 5 = 0 * ! (1) xy areunknowns,whichwillnowbedeterminedbysimplifying Eqs.(1.-15)and(1.-15),respectively,inthefollowing. Thus,startingfromEq.(1.-15), 0 m (1) 3 canberepresentedas 0 m (1) 3 = 1 , 3 & " % t 0 $ 2 3 # + # ( u 2 x + u 2 y ) % " % x $ 4 3 # u x + # u x u 2 y %

PAGE 122

107 " % y ( # u x ) " % y $ 4 3 # u y + # u 2 x u y % +2( F x u x + F y u y ) ' . (1.-14) Forsimplifying,theaboveequation(Eq.(1.-14)),weneedexpressionsofthehigherorderterms oftheform % t 0 < # u 2 i = and % i A # u i u 2 j B ,where i,j # x,y .FromEq.(1.-15), % t 0 ( # u x )= " 1 3 % x # " % x < # u 2 x = " % y ( # u x u y )+ F x ,andusingthisalongwithEq.(1.-15),i.e., % t 0 # = " % x ( # u x ) " % y ( # u y ) in % t 0 < # u 2 x = =2 u x % t 0 ( # u x )+ u 2 x % t 0 # toreplacethetimederivativesinfavorofspatialderivatives, weget % t 0 < # u 2 x = =2 u x & " 1 3 % x # " % x < # u 2 x = " % y ( # u x u y )+ F x ' + u 2 x [ % x ( # u x )+ % y ( # u y )] . Similarly, % t 0 < # u 2 y = =2 u y & " 1 3 % y # " % y < # u 2 y = " % x ( # u x u y )+ F y ' + u 2 y [ % x ( # u x )+ % y ( # u y )] . Simplifyingtheabove,wehave " % t 0 < # u 2 x = = 2 3 u x % x # +2 u x % x < # u 2 x = +2 u x % y ( # u x u y ) -2F x u x " u 2 x % x ( # u x ) " u 2 y % y ( # u y ) , (1.-13) " % t 0 < # u 2 y = = 2 3 u y % y # +2 u y % y < # u 2 y = +2 u y % x ( # u x u y ) -2F y u y " u 2 y % x ( # u x ) " u 2 y % y ( # u y )(1.-13)Also,inviewofEq.(1.-14),itfollowsthat " % x < # u x u 2 y = = " u 2 y % x ( # u x ) " 2 # u x u y % x u y . (1.-12) " % y < # u 2 x u y = = " u 2 x % y ( # u y ) " 2 # u y u x % y u x . (1.-12) Now,collectingallthehigherordertermsin 0 m (1) 3 giveninEq.(1.-14)togetherusingEqs.(1.13),(1.-13),(1.-12)and(1.-12),andsimplifyingbyretainingallcubicvelocityterms,whileneglectingtermsofnegligiblysmallhigherorders,weget " % t 0 < # u 2 x = " % t 0 < # u 2 y = " % x < # u x u 2 y = " % y < # u 2 x u y = = 2 3 ( u x % x # + u y % y # ) " 2( F x u x + F y u y )+3 # u 2 x % x u x +3 # u 2 y % y u y . (1.-12)

PAGE 123

108 Then,usingtheaboveequationEq.(1.-12)alongwith % t 0 # = " % x ( # u x ) " % y ( # u y )inEq.(A10), thesecond-orderdiagonalnon-equilibriummoment 0 m (1) 3 becomes 0 m (1) 3 = " 2 , 4 # ( % x u x " % y u y )+ E 4 g , (1.-12) where E 3 g . 3 # , 3 u 2 x % x u x . (1.-12) Followingasimilarprocedureforthenon-equilibriummomentcomponent 0 m (1) 4 startingfrom Eq.(1.-15)andafterconsiderablesimpliÞcation,weget 0 m (1) 4 = " 2 , 4 # ( % x u x " % y u y )+ E 4 g , (1.-12) where E 4 g . 3 # , 4 u 2 y % y u y . (1.-12) Inaddition,theo!-diagonalsecond-ordernon-equilibriummomentcomponent 0 m (1) 4 followsfrom Eq.(1.-15),as 0 m (1) 4 = " # 3 , 5 ( % x u y " % y u x ) . (1.-12) Intheabove,theterms E 3 g and E 4 g giveninEqs.(1.1.1)and(1.1.1),respectively,representthe spuriousnon-GalileaninvarianttermsthatarisefromthealiasingofthethirdorderdiagonalmomentsforthestandardD2Q9lattice.Inthenextsection,wewillpresentastrategytoe!ectively eliminatetheseerrors. 1.1.2 Eliminationofcubicvelocityerrorsviacorrectionstomomentequilibria Inordertoremovethenon-Galileaninvarianterrors E 3 g and E 4 g identiÞedintheprevioussection, weproposetoaddcorrectionsviaextendedmomentequilibriatothesecondorderdiagonalcom-

PAGE 124

109 ponentsas 0 m eq = * + + + + + + + + + + + + + + + + + + + + + + + + + + + , 0 m eq (0) 0 0 m eq (1) 1 0 m eq (1) 2 0 m eq (1) 3 0 m eq (1) 4 0 m eq (1) 5 0 m eq (1) 6 0 m eq (1) 7 0 m eq (1) 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . / + $ t * + + + + + + + + + + + + + + + + + + + + + + + + + + + , 0 0 0 0 m eq (1) 3 0 m eq (1) 4 0 0 0 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . / = * + + + + + + + + + + + + + + + + + + + + + + + + + + + , 0 * eq ! 0 0 * eq ! x 0 * eq ! y 0 * eq ! xx + 0 * eq ! yy 0 * eq ! xx " 0 * eq ! yy 0 * eq ! xy 0 * eq ! xxy 0 * eq ! xyy 0 * eq ! xxyy . . . . . . . . . . . . . . . . . . . . . . . . . . . . / + $ t * + + + + + + + + + + + + + + + + + + + + + + + + + + + , 0 0 0 3 x % x u x + 3 y % y u y 4 x % x u x " 4 y % y u y 0 0 0 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . / wherethecorrectiontothemomentequilibria 0 m eq (1) 3 and 0 m eq (1) 4 aregivenby 0 m eq (1) 3 = 3 x % x u x + 3 y % y u y (1.-11) 0 m eq (1) 4 = 4 x % x u x " 4 y % y u y (1.-11) Here,thecoe"cients 3 x , 4 y , 4 x and 4 y appearingasperfectorstothevelocitygradienttermswill bedeterminedviaamodiÞedC-Eanalysisinthefollowing.Thus,weconsidertheC-Eexpansionsusingthemomentequilibriacorrections,i.e., 0 m eq = 0 m eq (0) + 0 m eq (1) $ t ,as 0 m = 0 m eq (0) + 1 0 m eq (1) + 1 2 0 m eq (2) + ...,% t = % t 0 + 1% t 1 + 1 2 % t 2 + .... (1.-11) Then,followingthesameprocedureasintheprevioussection,wegetasequenceofmodiÞedmomentequationsatdi!erentorderin 1 as O ( 1 0 ): 0 m (0) = 0 m eq (0) , (1.-10a) O ( 1 1 ):( % t 0 + 0 E i % i ) 0 m (0) = " 0 ' 3 0 m (1) " 0 m eq (1) 4 + 0 S, (1.-10b) O ( 1 2 ): % t 1 0 m (0) +( % t 0 + 0 E i % i ) & I " 1 2 0 ' ' 0 m (1)

PAGE 125

110 + A % t 0 + 0 E i % i B & 1 2 0 ' 0 m eq (1) ' = " 0 ' 0 m (2) . (1.-10b) Whilethe O ( 1 )equationfortheÞrstordermomentsarethesameasthosegiveninEqs.(1.-15)Ð (1.-15),thepresenceofthecorrection 0 m eq (1) modiÞestheevolutionequationsforthecomponents ofthesecondordermoments,andhencetheEqs.(1.-15)Ð(1.-15)aretobereplacedbythefollowing: % t 0 $ 2 3 # + # ( u 2 x + u 2 y ) % + % x $ 4 3 # u x + # u x u 2 y % + % y $ 4 3 # u y + # u 2 x u y % = " , 3 0 m (1) 3 + , 3 0 m eq (1) 3 +2 F x u x +2 F y u y , (1.-9) % t 0 < # ( u 2 x " u 2 y ) = + % x $ 2 3 # u x " # u x u 2 y % + % y $ " 2 3 # u y + # u 2 x u y % = " , 4 0 m (1) 4 + , 4 0 m eq (1) 4 +2 F x u x " 2 F y u y , (1.-9) % t 0 ( # u x u y )+ % x $ 1 3 # u y + # u 2 x u y % + % y $ 1 3 # u x + # u x u 2 y % = " , 5 0 m (1) 5 + , 5 0 m eq (1) 5 + F x u y + F y u x , (1.-9) Then,theleading O ( 1 2 )momentequationsgiveninEqs.(1.-14)Ð(1.-14),owingtoEq.(1.-10b) arereplacedwith % t 1 # =0 , (1.-8a) % t 1 ( # u x )+ % x & 1 2 $ 1 " 1 2 , 3 % 0 m (1) 3 + 1 2 $ 1 " 1 2 , 4 % 0 m (1) 4 ' + % y &$ 1 " 1 2 , 5 % 0 m (1) 5 ' + % x & 1 4 , 3 0 m eq (1) 3 + 1 4 , 4 0 m eq (1) 4 ' + % y & 1 2 , 5 0 m eq (1) 5 ' =0 , (1.-8a) % t 1 ( # u y )+ % x &$ 1 " 1 2 , 5 % 0 m (1) 5 ' + % y & 1 2 $ 1 " 1 2 , 3 % 0 f (1) 3 " 1 2 $ 1 " 1 2 , 4 % 0 m (1) 4 ' + % x & 1 2 , 5 0 m eq (1) 5 ' + % y & 1 4 , 3 0 m eq (1) 3 " 1 4 , 4 0 m eq (1) 4 ' =0 , (1.-8a) FromEq.(1.-9),itfollowsthatthesecondordernon-equilibriummomentcomponent 0 m (0) 3 isthen modiÞedbytheuseofthecorrectionterm 0 m eq (0) 3 to ö m (1) 3 = 1 , 3 3 " 2 3 % t 0 # " % t 0 ( # u 2 x ) " % t 0 ( # u 2 y ) " 4 3 % x ( # u x ) " % x ( # u x u 2 y ) " 4 3 % y ( # u y )

PAGE 126

111 " % y ( # u 2 x u y )+2 F x u x +2 F y u y 4 +ö m eq (1) 3 , ThetermswithinthesquarebracketsintheaboveequationareidenticaltothosegiveninEqs.(1.1.1) and(1.1.1).Hence,wehave ö m (1) 3 = " 2 3 , 3 # ( % x u x + % y u y )+ E 3 g +ö m eq (1) 3 (1.-9) Similarly,itfollowsfromEq.(1.-9)and(1.1.1)and(1.1.1)thatthenon-equilibriummoment componentö m (1) 4 reducesto ö m (1) 4 = " 2 3 , 4 # ( % x u x " % y u y )+ E 4 g +ö m eq (1) 4 (1.-8) Theformofthenon-Galileaninvariantterms, E 3 g and E 4 g ,intheabovetwoequationsaregiven inEqs.(1.1.1)and(1.1.1),respectively,whilethemomentcorrectiontermsö m eq (1) 3 andö m eq (1) 4 are prescribedwithasyetunknowncoe"cientsinEqs.(1.-11)and(1.-11),respectively.Thencombiningleading O ( 1 )equationsgiveninEqs.(1.-15)Ð(1.-15)with 1 timesthecorresponding O ( 1 2 ) equations(Eqs.(1.-8a)Ð(1.-8a)),weget % t # + % x ( # u x )+ % y ( # u u )=0(1.-7) % t ( # u x )+ % x ( 1 3 # + # u 2 x )+ % y ( # u x u y )= F x " 1% x & 1 2 A 1 " , 3 2 B ö m (1) 3 + 1 2 A 1 " , 4 2 B ö m (1) 4 ' " 1% y 3A 1 " , 5 2 B ö m (1) 5 4 " 1% x 3 , 3 4 ö m eq (1) 3 + , 4 4 ö m eq (1) 4 4 " 1% y 3 , 5 2 ö m eq (1) 5 4 (1.-7) % t ( # u y )+ % y ( 1 3 # + # u 2 y )+ % x ( # u x u y )= F y " 1% x 3A 1 " , 5 2 B ö m eq (1) 5 4 " 1% y & 1 2 A 1 " , 3 2 B ö m (1) 3 " 1 2 A 1 " , 4 2 B ö m (1) 4 ' " 1% x 3 , 5 2 ö m eq (1) 5 4 " 1% y 3 , 3 4 ö f eq (1) 3 " , 4 4 ö m eq (1) 4 4 (1.-7) Inordertodeterminethemomentequilibriacorrections 0 m eq (1) 3 and 0 m eq (1) 4 weÞrstsubstitutefor thenon-equilibriummoments 0 m (1) 3 , 0 m (1) 4 and 0 m (1) 5 giveninEqs.(1.-9),(1.-8)and(1.1.1),re-

PAGE 127

112 spectively,intherighthandsideofthe x -momentumeqn(Eq.(1.-7)),whichthenbecomes = F x + 1% x D 1 3 $ 1 , 3 " 1 2 % # ( % x u x + % y u y )+ 1 3 $ 1 , 3 " 1 2 % # ( % x u x " % y u y ) E + 1% y D 1 3 $ 1 , 5 " 1 2 % # ( % x u y + % y u x ) E " 1% x D 1 2 A 1 " , 3 2 B E 3 g + 1 2 A 1 " , 4 2 B E 4 g E " 1% x D 1 2 0 m eq (1) 3 + 1 2 0 m eq (1) 4 E (1.-7) ItisclearthattheÞrsttwolinesintheaboveequation(Eq.(1.1.2))arerelatedtothephysicsof ßuidßowscorrespondingtothebodyforceandthestrainratecomponents,whiletheremaining ternsshouldsumtozero.Hencetheconditionstoeliminatethenon-Galileaninvariantterms E 3 g and E 4 g are A 1 " , 3 2 B E 3 g + 0 m eq (1) 3 =0(1.-6) A 1 " , 4 2 B E 4 g + 0 m eq (1) 4 =0(1.-6) whichconstraintherequiredmomentequilibriacorrections 0 m eq (1) 3 and 0 m eq (1) 4 . Then,substitutingEqs.(1.1.1),(1.1.1),(1.-11)and(1.-11)intheaboveequations,Eqs.(1.-6) and(1.-6),wecandeterminethecoe"cients 3 x , 3 y , 4 x and 4 y ,whicharesummarizedasfollows: 3 x = " 3 $ 1 , 3 " 1 2 % # u 2 x , 3 y = " 3 $ 1 , 3 " 1 2 % # u 2 y (1.-6) and 4 x = " 3 $ 1 , 4 " 1 2 % # u 2 x , 4 y = " 3 $ 1 , 4 " 1 2 % # u 2 y (1.-6) Itmaybenotedthattheaboveconstraintrelations(Eqs.(1.-6)and(1.-6))areidenticallysatisÞedtocorrectlyrecoverthe y -momentumaswell.Hence,withtheuseofthemomentequilibria corrections,i.e.,Eqs.(1.-11),(1.-11),(1.1.2)and(1.1.2),therecoveredmacroscopicßuidßow equationswithoutcubicvelocityerrorsarisingassimpliÞcationsofEqs.(1.-7)Ð(1.-7)aregiven by % t # + ! á j =0 , (1.-5)

PAGE 128

113 % x j x + ! á ( j u x )= " % x p + % x [ . 4 (2 % x j x " ! á j )+ . 3 ! á j ] + % y [ . 5 ( % x j y + % y j x )]+ F x , (1.-5) % y j y + ! á ( j u y )= " % y p + % x [ . 5 ( % x j y + % y j x )] + % y [ . 4 (2 % y% y " ! á j )+ . 3 ! á j ]+ F y , (1.-5) whichcorrespondtotheNSE,where p = 1 3 # isthepressureÞeld, j = # u ,andthebulkviscosity ( . 3 )andshearviscosity( . 4 = . 5 )aregivenby . 3 = 1 3 $ 1 , 3 " 1 2 % , . 4 = 1 3 $ 1 , 4 " 1 2 % , . 5 = 1 3 $ 1 , 5 " 1 2 % . (1.-5) 1.1.3 Strainratetensorcomponentsvianon-equilibriummoments Thecomponentsofthestrainratetensor % x u x , % y u y and( % x u x + % y u y )canberelatedtothelocalsecondordernon-equilibriummoments,whileaccountingforGalileaninvariantcorrections,as follows.First,fromEqs.(1.1.1),(1.-9),and(1.1.2),weget " c 1 % x u x " c 2 % y u y = 0 m (1) 3 . (1.-5) where c 1 = 2 # 3 , 3 & 1 " 9 4 u 2 x , 3 ' ,c 2 = 2 # 3 , 3 & 1 " 9 4 u 2 y , 3 ' . (1.-5) Analogously,from(1.1.1),(1.-8),and(1.1.2),weobtain " ÷ c 1 % x u x +÷ c 2 % y u y = 0 m (1) 4 . (1.-5) where ÷ c 1 = 2 # 3 , 4 & 1 " 9 4 u 2 x , 4 ' , ÷ c 2 = 2 # 3 , 4 & 1 " 9 4 u 2 y , 4 ' (1.-5) SolvingEqs.(1.1.3)and(1.1.3)forthediagonalcomponentsofthestrainratetensor,weget % x u x = 3 ÷ c 2 0 m (1) 3 + c 2 0 m (1) 4 4 [ " c 1 ÷ c 2 " ÷ c 1 c 2 ] , (1.-5) % y u y = 3 " c 1 0 m (1) 4 +÷ c 1 0 m (1) 3 4 [ " c 1 ÷ c 2 " ÷ c 1 c 2 ] . (1.-5)

PAGE 129

114 Finally,fromEq.(1.1.1),theo!-diagonalcomponentofthestrainratetensorcanbewrittenas % x u y + % y u x = " 1 c 3 0 m (1) 5 , (1.-5) where c 3 = # 3 , 5 . (1.-5) Thenon-equilibriummoments 0 m (1) 3 , 0 m (1) 4 and 0 m (1) 5 neededinEqs.(1.1.3),(1.1.3)and(1.1.3)can becomputedvia 0 m (1) 3 = " ! < e 2 ! x + e 2 ! y = f ! " & 2 3 # + # ( u 2 x + u 2 y ) ' , (1.-4) 0 m (1) 4 = " ! < e 2 ! x " e 2 ! y = f ! " # ( u 2 x " u 2 y ) , (1.-4) 0 m (1) 5 = " ! e ! x e ! y f ! " # u x u y . (1.-4)

PAGE 130

115 APPENDIXB 2.1 OrthogonalMatrixoftheMomentBasis K fortheD3Q19Model Thecomponentsoftheorthogonalmatrix K ofthemomentbasiswhichisderivedinSec.3.3.1 (seeEq.(3.8))isthuswrittenas K = ! " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # 100000000 ! 300000001200 110000011 ! 11 ! 400000 ! 4 ! 40 1 ! 10000011 ! 11400000 ! 4 ! 40 1010000 ! 11 ! 110 ! 40000 ! 42 ! 2 10 ! 10000 ! 11 ! 11040000 ! 42 ! 2 10010000 ! 2 ! 1100 ! 4000 ! 422 100 ! 10000 ! 2 ! 11004000 ! 422 11101000281101 ! 10111 1 ! 110 ! 100028 ! 110 ! 1 ! 10111 11 ! 10 ! 1000281 ! 10110111 1 ! 1 ! 10100028 ! 1 ! 10 ! 110111 11010101 ! 18101 ! 10111 ! 1 1 ! 1010 ! 101 ! 18 ! 10110111 ! 1 110 ! 10 ! 101 ! 1810 ! 1 ! 10 ! 111 ! 1 1 ! 10 ! 10101 ! 18 ! 10 ! 110 ! 111 ! 1 1011001 ! 1 ! 180110 ! 1 ! 11 ! 20 10 ! 1100 ! 1 ! 1 ! 180 ! 110 ! 1 ! 11 ! 20 101 ! 100 ! 1 ! 1 ! 1801 ! 10111 ! 20 10 ! 1 ! 1001 ! 1 ! 180 ! 1 ! 10 ! 111 ! 20 $ % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % &

PAGE 131

116 2.2 RawMomentsofSourcesfortheD3Q19Model ThediscreterawmomentsofsourcesofvariousorderscanbewrittenintermsoftheCartesian componentsoftheforceÞeldas 0 & ! 0 = % S ! | # $ =0 , 0 & ! x = % S ! | e ! x $ = F x , 0 & ! y = % S ! | e ! y $ = F y , 0 & ! z = % S ! | e ! z $ = F z , 0 & ! xx = % S ! | e 2 ! x $ =2 F x u x , 0 & ! yy = % S ! | e 2 ! y $ =2 F y u y , 0 & ! zz = % S ! | e 2 ! z $ =2 F z u z , 0 & ! xy = % S ! | e ! x e ! y $ = F x u y + F y u x , 0 & ! xz = % S ! | e ! x e ! z $ = F x u z + F z u x , 0 & ! yz = % S ! | e ! y e ! z $ = F y u z + F z u y , 0 & ! xyy = % S ! | e ! x e 2 ! y $ = F x u 2 y +2 F y u y u x , 0 & ! xzz = % S ! | e ! x e 2 ! z $ = F x u 2 z +2 F z u z u x , 0 & ! xxy = % S ! | e 2 ! x e ! y $ = F y u 2 x +2 F x u x u y , 0 & ! yzz = % S ! | e ! y e 2 ! z $ = F y u 2 z +2 F z u z u y , 0 & ! xxz = % S ! | e 2 ! x e ! y $ = F z u 2 x +2 F x u x u z , 0 & ! yyz = % S ! | e 2 ! y e ! z $ = F z u 2 y +2 F y u y u z , 0 & ! xxyy = % S ! | e 2 ! x e 2 ! y $ =2 F x u x u 2 y +2 F y u y u 2 x , 0 & ! xxzz = % S ! | e 2 ! x e 2 ! z $ =2 F x u x u 2 z +2 F z u z u 2 x , 0 & ! yyzz = % S ! | e 2 ! y e 2 ! z $ =2 F y u y u 2 z +2 F z u z u 2 y . (2.-17)

PAGE 132

117 2.3 ProjectionsoftheRawSourceMomentstotheOrthogonalBasisVectorsfor theD3Q19Model Theorthogonalprojectionsoftherawsourcemomentsinthevelocityspacecanbewrittenas 0 m s 0 = % K 0 | S ! $ =0 , 0 m s 1 = % K 1 | S ! $ = F x , 0 m s 2 = % K 2 | S ! $ = F y , 0 m s 3 = % K 3 | S ! $ = F z , 0 m s 4 = % K 4 | S ! $ =( F x u y + F y u x ) , 0 m s 5 = % K 5 | S ! $ =( F x u z + F z u x ) , 0 m s 6 = % K 6 | S ! $ =( F y u z + F z u y ) , 0 m s 7 = % K 7 | S ! $ =2( F x u x " F y u y ) , 0 m s 8 = % K 8 | S ! $ =2( F x u x + F y u y " 2 F z u z ) , 0 m s 9 = % K 9 | S ! $ =38( F x u x + F y u y + F z u z ) , 0 m s 10 = % K 10 | S ! $ =5 ( F x (3 u 2 x + u 2 y + u 2 z )+2 u x ( F y u y + F z u z ) ) " 9 F x , 0 m s 11 = % K 11 | S ! $ =5 ( F y ( u 2 x +3 u 2 y + u 2 z )+2 u y ( F x u x + F z u z ) ) " 9 F y , 0 m s 12 = % K 12 | S ! $ =5 ( F z ( u 2 x + u 2 y +3 u 2 z )+2 u z ( F x u x + F y u y ) ) " 9 F z , 0 m s 13 = % K 13 | S ! $ = ( F x ( u 2 y " u 2 z )+2 u x ( F y u y " F z u z ) ) , 0 m s 14 = % K 14 | S ! $ = ( F y ( u 2 x " u 2 z )+2 u y ( F z u z " F x u x ) ) , 0 m s 15 = % K 15 | S ! $ = ( F z ( u 2 x " u 2 y )+2 u z ( F x u x " F y u y ) ) , 0 m s 16 = % K 16 | S ! $ =42 ( F x u x ( u 2 y + u 2 z )+ F y u y ( u 2 x + u 2 z )+ F z u z ( u 2 x + u 2 y ) ) " 32[ F x u x + F y u y + F z u z ] , 0 m s 17 = % K 17 | S ! $ =6 ( F x u x ( u 2 y + u 2 z )+ F y u y ( u 2 x " 2 u 2 z )+ F z u z ( u 2 x " 2 u 2 y ) ) " 4[2 F x u x " ( F y u y + F z u z )] ,

PAGE 133

115 0 m s 18 = % K 18 | S ! $ =6 ( F x u x ( u 2 y + u 2 z )+ F y u y ( u 2 x " 2 u 2 z )+ F z u z ( u 2 x " 2 u 2 y ) ) " 4(2 F x u x " F y u y " F z u z ) . (2.-37) 2.4 SourceTermsinParticleVelocitySpacefortheD3Q19Model ThesourcetermsintheparticlevelocityspacecanbeobtainedbysolvingEq.(3.36)givenas follows(thevaluesof l 1 , l 2 , l 3 , l 4 , l 5 ,and l 6 are110,160,190,665,1197,and1260,respectively) S 0 = 1 399 [21 0 m s 0 " 5 0 m s 9 +19 0 m s 16 ] , S 1 = 1 20 l 5 [ l 6 0 m s 0 +2 l 5 0 m s 1 +3 l 4 0 m s 7 + l 4 0 m s 8 " l 1 0 m s 9 " 2 l 5 0 m s 10 " 2 l 3 0 m s 16 " 2 l 4 0 m s 17 ] , S 2 = 1 20 l 5 [ l 6 0 m s 0 " 2 l 5 0 m s 1 +3 l 4 0 m s 7 + l 4 0 m s 8 " l 1 0 m s 9 +2 l 5 0 m s 10 " 2 l 3 0 m s 16 " 2 l 4 0 m s 17 ] , S 3 = 1 20 l 5 [ l 6 0 m s 0 +2 l 5 0 m s 2 " 3 l 4 0 m s 7 + l 4 0 m s 8 " l 1 0 m s 9 " 2 l 5 0 m s 11 " 2 l 3 0 m s 16 + l 4 0 m s 17 " 3 l 4 0 m s 18 ] , S 4 = 1 20 l 5 [ l 6 0 m s 0 " 2 l 5 0 m s 2 " 3 l 4 0 m s 7 + l 4 0 m s 8 " l 1 0 m s 9 +2 l 5 0 m s 11 " 2 l 3 0 m s 16 + l 4 0 m s 17 " 3 l 4 0 m s 18 ] , S 5 = 1 20 l 5 [ l 6 0 m s 0 +2 l 5 0 m s 3 " 2 l 4 0 m s 8 " l 1 0 m s 9 " 2 l 5 0 m s 12 " 2 l 3 0 m s 16 + l 4 0 m s 17 +3 l 4 0 m s 18 ] , S 6 = 1 20 l 5 [ l 6 0 m s 0 " 2 l 5 0 m s 3 2 " 2 l 4 0 m s 8 " l 1 0 m s 9 +2 l 5 0 m s 12 " 2 l 3 0 m s 16 + l 4 0 m s 17 +3 l 4 0 m s 18 ] , S 7 = 1 40 l 5 [2 l 6 0 m s 0 +4 l 5 0 m s 1 +4 l 5 0 m s 2 +10 l 5 0 m s 4 +4 l 4 0 m s 8 + l 2 0 m s 9 + l 5 0 m s 10 + l 5 0 m s 11 +9 l 4 0 m s 13 " 9 l 4 0 m s 14 + l 3 0 m s 16 + l 4 0 m s 17 +3 l 4 0 m s 18 ] , S 8 = 1 40 l 5 [2 l 6 0 m s 0 " 4 l 5 0 m s 1 +4 l 5 0 m s 2 " 10 l 5 0 m s 4 +4 l 4 0 m s 8 + l 2 0 m s 9 " l 5 0 m s 10 + l 5 0 m s 11 " 9 l 4 0 m s 13 " 9 l 4 0 m s 14 + l 3 0 m s 16 + l 4 0 m s 17 +3 l 4 0 m s 18 ] , S 9 = 1 40 l 5 [2 l 6 0 m s 0 " 4 l 5 0 m s 1 " 4 l 5 0 m s 2 " 10 l 5 0 m s 4 +4 l 4 0 m s 8 + l 2 0 m s 9 + l 5 0 m s 10 " l 5 0 m s 11 +9 l 4 0 m s 13 +9 l 4 0 m s 14 + l 3 0 m s 16 + l 4 0 m s 17 +3 l 4 0 m s 18 ] , S 10 = 1 40 l 5 [2 l 6 0 m s 0 " 4 l 5 0 m s 1 " 4 l 5 0 m s 2 +10 l 5 0 m s 4 +4 l 4 0 m s 8 + l 2 0 m s 9 " l 5 0 m s 10 " l 5 0 m s 11 " 9 l 4 0 m s 13 +9 l 4 0 m s 14 + l 3 0 m s 16 + l 4 0 m s 17 +3 l 4 0 m s 18 ] , S 11 = 1 40 l 5 [2 l 6 0 m s 0 " 4 l 5 0 m s 1 +4 l 5 0 m s 3 +10 l 5 0 m s 5 +6 l 4 0 m s 7 " 2 l 4 0 m s 8 + l 2 0 m s 9 + l 5 0 m s 10 + l 5 0 m s 12 " 9 l 4 0 m s 13 +9 l 4 0 m s 14 + l 3 0 m s 16 + l 4 0 m s 17 " 3 l 4 0 m s 18 ] ,

PAGE 134

116 S 12 = 1 40 l 5 [2 l 6 0 m s 0 " 4 l 5 0 m s 1 +4 l 5 0 m s 3 " 10 l 5 0 m s 5 +6 l 4 0 m s 7 " 2 l 4 0 m s 8 + l 2 0 m s 9 " l 5 0 m s 10 + l 5 0 m s 12 +9 l 4 0 m s 13 +9 l 4 0 m s 14 + l 3 0 m s 16 + l 4 0 m s 17 " 3 l 4 0 m s 18 ] , S 13 = 1 40 l 5 [2 l 6 0 m s 0 +4 l 5 0 m s 1 " 4 l 5 0 m s 3 " 10 l 5 0 m s 5 +6 l 4 0 m s 7 " 2 l 4 0 m s 8 + l 2 0 m s 9 + l 5 0 m s 10 " l 5 0 m s 12 " 9 l 4 0 m s 13 " 9 l 4 0 m s 14 + l 3 0 m s 16 + l 4 0 m s 17 " 3 l 4 0 m s 18 ] , (2.-57) S 14 = 1 40 l 5 [2 l 6 0 m s 0 " 4 l 5 0 m s 1 " 4 l 5 0 m s 3 +10 l 5 0 m s 5 +6 l 4 0 m s 7 " 2 l 4 0 m s 8 + l 2 0 m s 9 " l 5 0 m s 10 " l 5 0 m s 12 +9 l 4 0 m s 13 " 9 l 4 0 m s 15 + l 3 0 m s 16 + l 4 0 m s 17 " 3 l 4 0 m s 18 ] , S 15 = 1 40 l 5 [2 l 6 0 m s 0 +4 l 5 0 m s 2 +4 l 5 0 m s 3 +10 l 5 0 m s 6 " 6 l 4 0 m s 7 " 2 l 4 0 m s 8 + l 2 0 m s 9 + l 5 0 m s 11 + l 5 0 m s 12 +9 l 4 0 m s 14 " 9 l 4 0 m s 15 + l 3 0 m s 16 " 2 l 4 0 m s 17 ] , S 16 = 1 40 l 5 [2 l 6 0 m s 0 " 4 l 5 0 m s 2 +4 l 5 0 m s 3 " 10 l 5 0 m s 6 " 6 l 4 0 m s 7 " 2 l 4 0 m s 8 + l 2 0 m s 9 " l 5 0 m s 11 + l 5 0 m s 12 " 9 l 4 0 m s 14 " 9 l 4 0 m s 15 + l 3 0 m s 16 " 2 l 4 0 m s 17 ] , S 17 = 1 40 l 5 [2 l 6 0 m s 0 +4 l 5 0 m s 2 " 4 l 5 0 m s 3 " 10 l 5 0 m s 6 " 6 l 4 0 m s 7 " 2 l 4 0 m s 8 + l 2 0 m s 9 + l 5 0 m s 11 " l 5 0 m s 12 +9 l 4 0 m s 14 +9 l 4 0 m s 15 + l 3 0 m s 16 " 2 l 4 0 m s 17 ] , S 18 = 1 40 l 5 [2 l 6 0 m s 0 " 4 l 5 0 m s 2 " 4 l 5 0 m s 3 +10 l 5 0 m s 6 " 6 l 4 0 m s 7 " 2 l 4 0 m s 8 + l 2 0 m s 9 " l 5 0 m s 11 " l 5 0 m s 12 " 9 l 4 0 m s 14 +9 l 4 0 m s 15 + l 3 0 m s 16 " 2 l 4 0 m s 17 ] . (2.-65) 2.5 Thechangeofdi!erentmomentsundercollisionfortheD3Q19lattice Thechangeofdi!erentmomentsundercollisionfollowfromtheorthogonalpropertyofthemomentbasismatrix K ,whicharegivenby " ! ( Ká 0 g ) ! e ! x e ! y = " ( % K ( | e ! x e ! y $ 0 g ( =4 0 g 4 , " ! ( Ká 0 g ) ! e ! x e ! z = " ( % K ( | e ! x e ! z $ 0 g ( =4 0 g 5 , " ! ( Ká 0 g ) ! e ! y e ! z = " ( % K ( | e ! y e ! z $ 0 g ( =4 0 g 6 ,

PAGE 135

117 " ! ( Ká 0 g ) ! e 2 ! x = " ( % K ( | e 2 ! x $ 0 g ( =6 0 g 7 +6 0 g 8 +42 0 g 9 , " ! ( Ká 0 g ) ! e 2 ! y = " ( % K ( | e 2 ! y $ 0 g ( = " 6 0 g 7 +6 0 g 8 +42 0 g 9 , " ! ( Ká 0 g ) ! e 2 ! z = " ( % K ( | e 2 ! z $ 0 g ( = " 12 0 g 8 +42 0 g 9 , " ! ( Ká 0 g ) ! e ! x e 2 ! y == " ( % K ( | e ! x e 2 ! y $ 0 g ( =4 0 g 10 +4 0 g 13 , " ! ( Ká 0 g ) ! e ! x e 2 ! z = " ( % K ( | e ! x e 2 ! z $ 0 g ( =4 0 g 10 " 4 0 g 13 , " ! ( Ká 0 g ) ! e 2 ! x e ! y = " ( % K ( | e 2 ! x e ! y $ 0 g ( =4 0 g 11 " 4 0 g 14 , " ! ( Ká 0 g ) ! e ! y e 2 ! z = " ( % K ( | e ! y e 2 ! z $ 0 g ( =4 0 g 11 +4 0 g 14 , " ! ( Ká 0 g ) ! e 2 ! x e ! z = " ( % K ( | e 2 ! x e ! z $ 0 g ( =4 0 g 12 +4 0 g 15 , " ! ( Ká 0 g ) ! e 2 ! y e ! z == " ( % K ( | e 2 ! y e ! z $ 0 g ( =4 0 g 12 " 4 0 g 15 , " ! ( Ká 0 g ) ! e 2 ! x e 2 ! y = " ( % K ( | e 2 ! x e 2 ! y $ 0 g ( =8 0 g 8 +32 0 g 9 +4 0 g 16 +4 0 g 17 +4 0 g 18 , " ! ( Ká 0 g ) ! e 2 ! x e 2 ! z = " ( % K ( | e 2 ! x e 2 ! z $ 0 g ( =4 0 g 7 " 4 0 g 8 +32 0 g 9 +4 0 g 16 +4 0 g 17 " 4 0 g 18 , " ! ( Ká 0 g ) ! e 2 ! y e 2 ! z = " ( % K ( | e 2 ! y e 2 ! z $ 0 g ( = " 4 0 g 7 " 4 0 g 8 +32 0 g 9 +4 0 g 16 " 8 0 g 17 . (2.-78)