Analytical and Numerical Study ofPlug Flow Inside Round/ConcentricMicrochannelsbyYadi CaoB. Eng., University of Science and Technology of China, 2016A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)September 2018c© Yadi Cao, 2018The following individuals certify that they have read, and recommend to the College ofGraduate Studies for acceptance, a thesis entitled:Analytical and Numerical Study of Plug Flow Inside Round/Concentric Microchannelssubmitted by Yadi Cao in partial fulfillment of the requirements forthe degree of Master of Applied Sciencein Mechanical EngineeringDr. Sunny Ri Li, School of EngineeringSupervisorDr. Joshua Brinkerhoff, School of EngineeringSupervisory Committee MemberDr. Clarie Yu Yan, School of EngineeringSupervisory Committee MemberDr. Kenneth Chau, School of EngineeringUniversity ExamineriiAbstractPlug flow can generate recirculating flow between interfaces of immiscible fluids. De-pending on the phase used to segment flow, it can be gas-liquid plug flow or liquid-liquidplug flow. The recirculating flow can enhance heat transfer as compared to the continuouspipe flow, especially in the laminar regime such as the microchannel flow. The present workfocuses on the flow and heat transfer of liquid plugs with low Reynolds numbers.The flow is modeled by applying Stokes simplification, and the solution is obtained bysolving fourth-order partial differential equation sets. Solutions of two types of plug floware obtained: 1) the gas-liquid plug flow in the concentric microchannel; 2) the liquid-liquidplug flow in the circular microchannel. For the gas-liquid plug flow study, the flow patternsinside the liquid phase including the volume ratio of the inner and outer vortexes, the ratioof maximum-to-minimum stream functions, the averaged recirculation flux as well as theskin friction coefficient are investigated in details. Correlations for predicting the maximumand minimum of the stream function are developed. For the liquid-liquid plug flow study,the influences of plug lengths and the viscosity ratio upon the cap vortexes and the overallskin friction coefficient are studied in details.The heat transfer of the gas-liquid plug flow in the concentric microchannel is simulatednumerically in MATLAB. Three types of thermal boundary conditions are investigated. Thedeveloping process of the thermal field can be explained using a simple thermal network foreach boundary condition. The influences of parameters including the plug aspect ratio, thechannel inner-outer radius ratio and the Peclect number upon the thermal conductance andheat transfer enhancement to the single-phase flow are investigated systematically. Then asimplified model for the fully developed thermal field is extracted for the quick calculationneed in the design work. The results obtained from about 12,000 cases form a databasethat can be used in the future design work of heat exchanger based upon this kind of flow.iiiLay SummaryPlug flows can enhance the heat transfer inside the microchannels by generating recir-culating flow between interfaces of immiscible fluids. In this thesis, the analytical solutionsof the two kinds of plug flow are obtained, and the physics behind it is revealed. The flowfield results are used in heat transfer simulation under a wide range of working conditions.The obtained database of both flow field and heat transfer results can help the design workof heat exchanger based on this configuration in the future.ivPrefaceThe work outlined in this thesis was conducted by Yadi Cao under the supervision of Dr.Sunny Ri Li. It was supported by Natural Sciences and Engineering Research Council ofCanada (NSERC). All the presented research work was finished in the Thermal Managementand Multiphase Flows Laboratory in the School of Engineering at the University of BritishColumbia (Okanagan campus). Portions of this thesis have been published in journals aswell as in conference proceedings.• A part of chapter 2 has been published as a proceeding at a conference. Cao, Y.D., Li,R. Liquid Plug in Gas Flow in Annular Channel. 3rd American Society of Thermaland Fluids Engineers, Mar. 2018. (TFEC-2018-24337).• A part of chapter 2 has been accepted by the Physics of fluids. Cao, Y.D., Li, R.A Liquid Plug Moving in an Annular Pipe – Flow Analysis. Physics of fluids. Jul.2018. (POF18-AR-01493R).vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xviiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xixChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Challenges faced by cooling in microchannel . . . . . . . . . . . . . . . . . . 21.3 Numerical/experimental studies of plug flow in microchannel . . . . . . . . 61.4 Analytical models of plug flows . . . . . . . . . . . . . . . . . . . . . . . . . 91.5 Motivation and Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.6 Organization of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11Chapter 2: Analytical study of gas-liquid plug flow in concentric microchan-nel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13viTABLE OF CONTENTS2.1 Mathematic modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.1.1 Governing equations and boundary conditions . . . . . . . . . . . . . 132.1.2 Analytical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.2 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.2.1 Verify solution accuracy using CFD Results . . . . . . . . . . . . . . 192.2.2 Validate using continuous flow(1-D model) . . . . . . . . . . . . . . . 212.2.3 Two asymmetric vortexes . . . . . . . . . . . . . . . . . . . . . . . . 222.2.4 Quantify radial transport using the averaged recirculation flux . . . 252.2.5 Recirculation period . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.2.6 The skin friction coefficient . . . . . . . . . . . . . . . . . . . . . . . 292.3 Summary of this chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30Chapter 3: Numerical study of heat transfer inside plug flow in chapter 2at asymmetric boundary conditions . . . . . . . . . . . . . . . . . 323.1 Modeling for transient heat transfer process in plug flow . . . . . . . . . . . 323.1.1 Governing equation and initial/boundary conditions . . . . . . . . . 323.1.2 Nondimensionalization . . . . . . . . . . . . . . . . . . . . . . . . . . 353.2 Modeling for the asymptotic case - continuous Stokes flow heat transfer . . 363.3 Heat transfer performance at inner iso-flux boundary (IFB) . . . . . . . . . 383.3.1 Processes of heat transfer and its simplified thermal network . . . . 393.3.2 Influences of the inner-outer radius ratio η . . . . . . . . . . . . . . . 423.3.3 Influences of the aspect ratio β . . . . . . . . . . . . . . . . . . . . . 443.3.4 Influences of the Peclet number Pe . . . . . . . . . . . . . . . . . . . 443.4 Heat transfer performance at the inner iso-thermal boundaries (ITB), com-parison to IFB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.5 Heat transfer performance at the outer iso-flux boundaries (OFB), compari-son to IFB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.6 Results for steady heat transfer under fully developed stage . . . . . . . . . 523.6.1 Simplification based upon control volume method . . . . . . . . . . . 523.6.2 Some results and potential future design work . . . . . . . . . . . . . 543.7 Summary of this chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62viiTABLE OF CONTENTSChapter 4: Analytical study of liquid-liquid plug flow in circular microchan-nel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 644.1 Mathematic modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 644.1.1 Governing equations and boundary conditions . . . . . . . . . . . . . 644.1.2 Analytical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.2 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.2.1 Validate using previous numerical study [1] . . . . . . . . . . . . . . 714.2.2 Validate using the existing gas-liquid model [2] . . . . . . . . . . . . 724.2.3 Influence of viscosity ratio µ2/µ1 . . . . . . . . . . . . . . . . . . . . 744.2.4 Influence of plug lengths β1, β2 . . . . . . . . . . . . . . . . . . . . . 774.2.5 Skin friction coefficients Cf . . . . . . . . . . . . . . . . . . . . . . . 794.3 Summary of this chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81Chapter 5: Conclusions, limitations and potential future work . . . . . . . 84Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95Appendix A: Correlations for long plugs (β > 2) in chapter 2 . . . . . . . . . . . 96viiiList of TablesTable 3.1 Set-up for numerical solution of heat transfer . . . . . . . . . . . . . 36Table 3.2 Set-up for relaxation ratio (SOR) for steady state heat transfer . . . 54ixList of FiguresFigure 1.1 Schematic show for the recirculation in plug flow. The contour plotis from the results in chapter 2 in this thesis. . . . . . . . . . . . . . 5Figure 1.2 Schematic show for stages of the heat transfer of gas-liquid plug flowin a slit microchannel in [3]. For plug flow: (I) the thermal entranceregion. (II) the transition region. (III) the fully developed region.For single-phase flow there is no transition region. Reprinted withpermission. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6Figure 2.1 Schematic show for gas-liquid plug flow in concentric microchannel 13Figure 2.2 The sensitivity study for the analytical solution. Variation of Cf withthe growing series number. β = 2, η = 0.5. . . . . . . . . . . . . . . 20Figure 2.3 The velocity vector plot of (a). the analytical model and (b). theCFD results, β = 1, η = 0.50 . . . . . . . . . . . . . . . . . . . . . . 21Figure 2.4 The comparison between plugs with different aspect ratio β and the1-D model of (a) the axial velocity uˆz. (b) The stream function ψˆ.The sample line is at zˆ = 0.25 lˆ, and the inner-outer radius ratio isη = 0.90. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22Figure 2.5 The contour plots for stream function in plugs with different η, thecross shape marker are the locations for vortex centers. . . . . . . . 23Figure 2.6 The volume ratio between the inner and the outer vortex Vˆin,v/Vˆout,v. 24Figure 2.7 The index for relative location in radial direction of vortex centers. 25Figure 2.8 The stream function ratio between the inner and the outer vortexcenters ψˆmin/ψˆmax. . . . . . . . . . . . . . . . . . . . . . . . . . . . 26xLIST OF FIGURESFigure 2.9 The recirculation flux rˆuˆr of the inner vortex and the outer vortex.Main plot is the averaged value, the small subplots show where thesample line is (right-top), and a detailed distribution on the sampleline for β = 2, η = 0.5 (right-bottom), respectively. . . . . . . . . . 28Figure 2.10 The distribution of recirculation time tˆ on the central of the plug(zˆ = lˆ/2). The double-arrows mark the boundaries between theinner and the outer vortexes. β = 1, 2, 3, 4, η = 0.5. . . . . . . . . 29Figure 2.11 The skin friction coefficient (a). versus β for selected η and (b).versus η for selected β . . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 3.1 The schematic show of boundary conditions (a). The inner iso-flux boundaries (IFB). (b). The iso-thermal boundaries (ITB). (c).The outer iso-flux boundaries (OFB). The contours are taken attˆ = (1, 6, 25, 40) (eq. (2.43)) under these 3 boundary types corre-spondingly. (β, η,Pe) = (2, 0.5, 100). . . . . . . . . . . . . . . . . . . 34Figure 3.2 Heat transfer process at IFB, (β, η,Pe) = (2, 0.5, 100). (a). (I) ∼(IV) Sequence of dimensionless temperature Tˆ distribution at tˆI∼IV =(1, 6, 25, 40). (b). Variation of dimensionless temperature at theinner wall Tˆin,w, at the outer wall Tˆout,w and the mean valueˆ¯T ,and the mixing index γ. The black solid points correspond to theexamples shown in (a). (c). Simplified thermal network for IFB. . . 39Figure 3.3 (a). Development of the thermal conductance σˆplug, (b). influenceof the inner-outer radius ratio η upon σˆplug and the enhancement byplug flow comparing to the single-phase flow heat transfer σˆplug/σˆsp.(β,Pe) = (2, 100) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Figure 3.4 Influence of plug dimensionless length β upon (a). the thermal con-ductance σˆplug, (b). the enhancement by plug flow comparing to thesingle-phase flow heat transfer σˆplug/σˆsp. The arrows in figures markthe direction of increasing β. Pe = 100. . . . . . . . . . . . . . . . . 45xiLIST OF FIGURESFigure 3.5 Influence of Peclet number Pe upon (a). the thermal conductanceσˆplug, (b). the enhancement by plug flow comparing to the single-phase flow heat transfer σˆplug/σˆsp. The arrows in figures mark thedirection of increasing Pe. β = 2. . . . . . . . . . . . . . . . . . . . 47Figure 3.6 Heat transfer process at ITB, (β, η,Pe) = (2, 0.5, 100). (a). Variationof the mean temperature ˆ¯T , and the dimensionless heat flow at theinner wall Qˆin, at the outer wall Qˆout and at the boundary betweentwo vortexes Qˆb. The black solid points correspond to the examplesshown in (b)., (b). (I) ∼ (IV) Sequence of dimensionless temperatureTˆ distribution at tˆI∼IV = (0.5, 3, 9.5, 25). . . . . . . . . . . . . . . . 48Figure 3.7 Comparison between IFB and ITB (a) (b). the thermal networks,(c). the thermal conductance σˆplug, (d). the enhancement by plugflow comparing to the single-phase flow heat transfer σˆplug/σˆsp. Theworking condition is the same as in fig. 3.3. . . . . . . . . . . . . . . 49Figure 3.8 Comparison between IFB and OFB (a) ∼ (b). the thermal networks,(c). the thermal conductance σˆplug, (d). the enhancement by plugflow comparing to the single-phase flow heat transfer σˆplug/σˆsp. Theworking condition is the same as in fig. 3.3. . . . . . . . . . . . . . . 51Figure 3.9 Temporal derivative of the thermal field ∂Tˆ∂tˆat tˆ = 50 ,at IFB, allother parameters are same as in fig. 3.2 . . . . . . . . . . . . . . . . 52Figure 3.10 Comparison between (a). thermal field under simplified model forfully developed stage, and (b). that at tˆ = 60 using transient method,at OFB, all other parameters are same as in fig. 3.2 . . . . . . . . . 53Figure 3.11 I of II: The sequence of plots for Nu under Pe = 1, 2, 5, 10, 20, 50 atIFB condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56Figure 3.12 II of II: The sequence of plots for Nu under Pe = 100, 200, 500, 1000, 2000at IFB condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Figure 3.13 I of II: The sequence of plots for Nu under Pe = 1, 2, 5, 10, 20, 50 atITB condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58xiiLIST OF FIGURESFigure 3.14 II of II: The sequence of plots for Nu under Pe = 100, 200, 500, 1000, 2000at ITB condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59Figure 3.15 I of II: The sequence of plots for Nu under Pe = 1, 2, 5, 10, 20, 50 atOFB condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60Figure 3.16 II of II: The sequence of plots for Nu under Pe = 100, 200, 500, 1000, 2000at OFB condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61Figure 4.1 Schematic show for liquid-liquid plug flow in microchannel with roundcross-section. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Figure 4.2 Schematic show for Stokes flow in a cavity with moving walls, ademonstration of the problem in [4]. F ,G are generic boundaryconditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66Figure 4.3 The sensitivity study for the analytical solution. Variation of aver-aged Cf with the growing series number. β1 = 1, β2 = 1, µ2/µ1 = 16. 71Figure 4.4 Comparison between plots of flow field (a) from simulation in [1] (re-printed with permission) and (b) from the analytical solution in thiswork. β1 = 6.0, β2 = 4.8, µ2/µ1 = 1.2076. . . . . . . . . . . . . . . . 72Figure 4.5 Comparison between contours of stream function (a) from the an-alytical solution of the gas-liquid model in [2] (here only the equa-tions in this citation are used to plot the contour) and (b) from theanalytical solution of liquid-liquid model in this work. (c) The dis-tribution of axial velocity uˆz on 3 sample lines zˆ = 1.00, 0.40, 0.20.β1 = 1, β2 = 2, µ2/µ1 = 64 (β1, µ2/µ1 are not needed for (a)). . . . 73Figure 4.6 Schematic show of the reason for cap vortexes when µ2/µ1 is suffi-ciently large. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74Figure 4.7 Sequence of stream lines with increasing µ2/µ1 = (1, 2, 4, 16), whereβ1 = β2 = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 4.8 The influence of increasing µ2/µ1 upon (a) ψˆmax the intensity of mainvortex centers , (b) ψˆmin the intensity of cap vortex centers in plug1 and (c) Cˆcap the capacity (volume) of the cap vortex in plug 1. . . 76xiiiLIST OF FIGURESFigure 4.9 Sequence of stream lines with increasing β1 = (1, 2, 4), where β2 = 1and µ2/µ1 = 16. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77Figure 4.10 The influence of varying β1 upon (a) ψˆmax the intensity of mainvortex centers , (b) ψˆmin the intensity of cap vortex centers in plug1 and (c) Cˆcap the capacity (volume) of the cap vortex in plug 1. . . 78Figure 4.11 The influence of varying β2 upon (a) ψˆmax the intensity of mainvortex centers , (b) ψˆmin the intensity of cap vortex centers in plug1 and (c) Cˆcap the capacity (volume) of the cap vortex in plug 1. . . 79Figure 4.12 The skin friction coefficients for plug 1, plug 2 and their mean valuewith varying µ2/µ1, β1 = β2 = 1. . . . . . . . . . . . . . . . . . . . . 80Figure 4.13 Variations of mean skin friction Cf coefficients with increasing µ2/µ1,β1 = β2 = 1 under (a)β2 = 1, β1 = 1 ∼ 10, (b)β1 = 1, β2 = 1 ∼ 10 . . 82Figure A.1 The volume ratio of two vortexes Vˆin,v/Vˆout,v, and the ratio of max/minstream functions |ψˆmin/ψˆmax| when β > 2. The solid lines are thecorrelations for the points obtained from analytical results in chapter 2. 96Figure A.2 The total of volumetric recirculation rate of two vortexes ψˆmax−ψˆminbecomes independent of β when β > 2. . . . . . . . . . . . . . . . . 97Figure A.3 The refresh frequency defined in eq. (A.3). The solid lines are thecorrelations for the points obtained from analytical results in chapter 2. 98Figure A.4 Comparison between the correlation and the analytical results. Thesolid lines are the correlations for the points obtained from analyticalresults in chapter 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100xivList of SymbolsSymbols DefinitionsParametersα thermal diffusivityβ aspect ratio ratio (eq. (2.11))C capacity for extensive propertiesCf skin friction coefficient (eq. (2.32))cp specific heatDh hydraulic diameterη inner-outer radius ratio (eq. (2.11))Ff friction force (eq. (2.29))γ mixing index (eq. (3.24))k thermal conductivityl length, characteristic length (eq. (1.3))λ surface tensionµ dynamic viscosityP static pressureψ stream functionQ heat flowq′′ heat fluxRc radius of curvatureρ densityσ thermal conductance (eqs. (3.25) and (3.26))t timexvList of Symbolsu, u, U velocity scalar, velocity vector, plug moving/mean velocityV volumez, r axial index, radial indexPopular dimensionless parametersCa Capillary number (eq. (2.1))De Dean number (eq. (1.2))Ec Eckert number (eq. (3.2))f Fanning friction factor (eq. (2.31))Nu Nusselt number (eq. (3.32))Pe Peclet number (eq. (3.10))Pr Prandtl number (eq. (1.4))Re Reynolds number (eq. (1.1))We Weber number (eq. (1.3))Special values and coefficientsA,B,C,D,E, F constant coefficientsω, χ eigenvaluesSpecial functions and operators∆ difference operatorH finite Hankel transformation (eq. (4.13))I modified Bessel function of the firstJ Bessel function of the first kindK modified Bessel function of the second kindL−1 a linear operator (eq. (2.6))∇ Nabla operator∇2 Laplace operatorS finite Fourier transformation (eq. (2.16))xviList of SymbolsSuperscriptsˆ dimensionless value¯ averaged value (eq. (3.23))Subscriptscap (for) cap vortexes (inside plug)f fluid (eqs. (3.25) and (3.26)) or frictional (eq. (2.29))i, j makers for phases in chapter 4 (i, j = 1, 2 ∧ i 6= j)in, out (at) inner side, (at) outer sidel, m, n series numberplug (for) plugrec (for) recirculation (inside plug)sp, 1D (for) single phased flow, (for) 1 dimensional flowstd a standard value (for nondimensionalization)w, v (of) wall, (of) vortexxviiAcknowledgmentsI gratefully appreciate the patient guidance from Dr. Ri Li.The project is fully funded by the Natural Sciences and Engineering Research Council(NSERC).I would also like to especially thank Dr. Xuan Gao for providing the initial motivationfor the subject of this thesis.It is an honor to have Dr. Joshua Brinkerhoff and Dr. Yu Yan as the committeemembers. Their advice during the annual report helped optimize this work.xviiiDedicationTo the great love and the unconditional support frommy parents.xixChapter 1Introduction1.1 BackgroundDriven by the fast development of applications including microelectronics,micro-reactors,micro-electro-mechanical systems(MEMS), the study of microfluidics has drawn much atten-tion in recent years [5]. The applications of microfluidics includes environmental detection[6], chemical reaction [7–9], and drug delivery [10]. Despite the efficiency brought by MEMSand microfluidics, one of the biggest challenges faced by engineers is the heat generationin these devices. Traditional cooling methods have faced challenges owing to increasingheat flux and thermal resistance caused by the complicated internal structure of micro-sizesystems. Thus, the large heat flux removal methods for these devices have gained muchattention in recent years [11]. Many modifications of traditional heat transfer methods havebeen done to suit the requirement of micro-sized systems. These modified methods includespray cooling [12–14], jet impingement cooling [15], falling film cooling [16] and heat pipes[17]. However, it has been pointed out that these methods are neither easily controllablenor cost-friendly [11, 18]. Other alternatives are needed for low-cost cooling methods.The heat exchanger using microchannels has the advantage of small volume and largesurface-volume ratio because of the small scale [19]. It also can be parallel configured andpackaged to cooling capacity. Tuckerman and Pease [20] introduced the concept of themicrochannel heat exchanger and they found out the heat removal limit for single phaseflow was 0.79 kW/cm2. However, this limit is relatively lower comparing to the heat fluxof current super-powered computing chips [21–23]. New efficient methods are needed forthermal design need in the future.The heat transfer limit in the work of Tuckerman and Pease [20] can be explained usinglow Reynolds number (Re) (defined below) in microchannels. The Reynolds number (Re)11.2. Challenges faced by cooling in microchannelis often used to determine the flow patterns in different situations, and it is defined usingthe ratio of inertial force and viscous force. In the case of internal flow, it can be writtenas belowRe =ρDhUµ, (1.1)where ρ is the density of the fluid, Dh is the hydraulic diameter of the channel, and Uis the mean velocity of the fluid, and µ is the dynamic viscosity of the fluid.In microchannel, it requires extremely high pumping power to make Re high enough(∼ 2, 000) to reach the turbulent regime owing to the small scale of Dh. Thus the laminarflow is often preferred, and heat transfer will be limited by the thermal diffusion processif there is no vortex inside to disturb the flow field. Since the extensive quantities likeenergy can be transported with flowing fluid, the mixing process including heat transfercan be enhanced by vortexes which disturb the flow field. Lots of previous studies havebeen focusing on generating vortexes inside single-phase flow. Another alternative, whichhas also gained attention, is to involve flow boiling inside the microchannel to make full useof the latent heat. In the next section, previous works related to both the options mentionedabove will be reviewed.1.2 Challenges faced by cooling in microchannelMost studies have been aiming to generate vortexes in the single-phase flow by eitherputting vortex promoters or varying geometries of the channel. However, their studies stillrequire quite high Re (100 ∼ 5, 600) in order to generate strong vortexes and have anapparent heat transfer enhancement. Hence, these methods will face a similar challenge ofextremely high pumping power as triggering turbulent flow in microchannel does. Thus theapplication of these methods can be limited. The reviews of some selected works related tothe above methods are presented below.Vortex generators are obstacles with different geometry in the channel, and they canproduce transverse flow due to flow separation on the interface between the obstacles andthe fluid. In the numerical study of Meis, Varas, Velazquez, and Vega[24], obstacles with21.2. Challenges faced by cooling in microchanneldifferent shapes and orientations were installed inside a microchannel as vortex promoters.The Reynolds number required in their work above vary within 600 ∼ 1, 200, and thisdesign has a slightly higher heat transfer enhancement than pumping power penalty. In theexperimental study of Wang, Houshmand, Elcock and Peles [25], Wang, Nayebzadeh, Yu,Shin and Peles[26], micropillar of different shapes/pin fins were put into the microchannel,respectively. The micro particle image velocimetry (µPIV) technique was used to observeflow structure downstream of a pillar, while resistance temperature detector (RTD) wasused to measure the spatially averaged temperature. The triangular pillar has the bestenhancement with Nusselt number vary within 17.7 ∼ 88.9, while the corresponding rangefor Reynolds number is 100 ∼ 5, 600 in order to generate vortexes in the downstream of thepillar. Nevertheless, they do not include information about pressure loss.Another way of introducing vortexes is by making the channel curved to obtain Deanvortexes. Dean vortexes are a kind of secondary flows (which always appear in a pair ofcounter-rotating vortexes) generated by the imbalance between centrifugal forces and radialpressure [27]. The dimensionless Dean number (De) can be used to describe the strengthof Dean vortexes. It is defined as belowDe = Re√Dh2Rc, (1.2)where Re is the Reynolds number, Dh is the hydraulic diameter, and Rc is the radiusof curvature of the path of the channel.In the numerical studies of Sui, Teo, Lee, Chew, and Shu[28], wavy channel shapes wereapplied to generate Dean vortexes inside the microchannel. The Reynolds number requiredin these two works is within 100 ∼ 400 and 300 ∼ 650 respectively. The enhancement inNusselt number is about 1.67 ∼ 2.02 times depending on the tunnel shape, while the pressureloss penalty is slightly lower than the heat transfer enhancement. In the numerical studyof Xia, Jiang, Wang, Zhai, and Ma[29], besides applying wavy channel, reentrant cavitiesat the wall were added to provide more heat transfer area. Nevertheless that the drag areaincreased simultaneously, and the authors did not mention the efficiency of add-on cavities.Chai, Xia, Wang, Zhou, and Cui[30] investigated both numerically and experimentallythe flow heat transfer inside microchannels with periodically expansion-constriction cross-31.2. Challenges faced by cooling in microchannelsections. The interval for Re in their work is within 300 ∼ 750, and the Nusselt numbercan be enhanced by up to 180%, they did not provide details about pressure loss as well.Experimental studies about flow in microchannel involving boiling have also been carriedout in recent years due to the considerable latent heat of fluid. Despite many applicationsof cooling with phase change in electronics, the physics behind phase change is not clear,which leads to high complexity to flow control [18]. The reviews of selected works relatedto above are presented below.The heat removal in a mini-channel can be as high as 10 kW/cm2 [31], where DI (de-ionized) water was used and the mass flow rate is 5, 000 ∼ 134, 000 kg/m2s−1. However,in the microchannel, the mass flow rate is much lower, and the critical heat flux (CHF)usually cannot reach this ultra high value. In the work of Deng, Wan, Qin, Zhang andChu[32], Krishnamurthy and Peles[33], structured microchannel with micro pin fins (SM-MPF) were fabricated using a laser micro-milling method. Lots of tiny reentrant cavitiesincrease nucleation site density (the number of sites where nucleation happens in a givenarea) for boiling. The mass flow rate is 200 ∼ 300 kg/m2s−1 and CHF can be 0.01 ∼0.11 kW/cm2, which results in 10% ∼ 175% enhancement comparing to flow boiling instraight microchannel depending on working liquids. Liu, Li, Liu and Gau[34] and Choi,Shin, Yu and Kim[35] investigated the influence of wettability of the wall upon boiling inthe microchannel and found out that the heat transfer in hydrophobic microchannels ishigher than that in hydrophilic ones. In hydrophobic channels, the nucleation site densityis observed to be higher, and more departure bubbles can disturb the flow field inside theliquid film, making the liquid film unstable and enhance the heat transfer. Jaikumar andKandlikar[36] aligned nucleating regions with non-nucleating ones in the feeder microchannelto separate pathways for returning liquid and upward floating vapor. The separate pathwaysled the fluid to refill the boiling nucleation to hinder the dry-out and thus can enhance theboiling efficiency. They found out an optimal alignment for enhancement, and the criticalheat flux can be 0.39 kW/cm2 under this alignment.Overall, introducing vortexes in the single-phase flow has the requirement of high Rewhile the mechanism of flow boiling is still not well understood, which makes flow boilinghard to predict and control. Thus, stable multiphase flow without phase change is preferred41.2. Challenges faced by cooling in microchannelfor prediction and easier control [18]. Among all the multiphase flow types without phasechange, the plug flow has drawn attention in recent years [37]. The simple introduction toplug flow and the vortexes inside is made below.Figure 1.1: Schematic show for the recirculation in plug flow. The contour plot is from theresults in chapter 2 in this thesis.Plug flow is also called slug flow, Taylor flow or segmented flow, it is a flow structurewhere another kind of immiscible fluid segments the liquid into separate plugs. The liquidplug has a shape of plug and nearly takes all the cross-section of the channel [38]. In thereference sticking to a moving plug, the non-slip wall becomes a moving wall with a velocityin the transverse direction. The moving wall drives the fluid by viscous force, and then theflowing fluid bounces back at the interface at the rear end of the plug, then travels alongthe middle axial of the plug and bounces back again at the frontal end and finishes a circleof recirculation(fig. 1.1). Comparing to using vortex promoter or generating Dean vortexes,this method prefers the lower speed of the plug[39], and thus extremely high pumping poweris not required. To explain this, a new dimensionless parameter, Weber number is writtenas below, where l is the characteristic length and λ is the surface tension,We =ρ u2 lλ. (1.3)Weber number (We) represents the ratio between inertia and surface tension. When thevelocity of the plug is small, We is small, and the interface tends to be more solid comparingto the momentum brought by the circulation, and thus the flow pattern inside tends to bemore stable. Besides this advantage of not requiring high pumping power, plug flow also hasthe advantage of vortex size. Comparing to vortex promoters where vortexes are generated51.3. Numerical/experimental studies of plug flow in microchannelonly after the obstacles, the vortexes in a liquid plug nearly takes up the whole cross-section[37] which indicates the possibility of better mixing ability.1.3 Numerical/experimental studies of plug flow inmicrochannelThe heat transfer inside the plug flow has been investigated experimentally in the recentdecades. In the experimental work of [40], Nitrogen gas was injected coaxially into DI (de-ionized) water for stable two-phase flow. The high-speed camera and thermocouples alongthe channel were used to measure the flow field and the temperature field respectively. Theenhancement of the Nusselt number is 176% while the pressure loss penalty is around 22%.Though the intrusive thermocouples may disrupt the flow field and furthermore influencethe heat transfer results. What is more, thermocouples can only collect the temperature atlimited points instead of alongside the whole channel. Thus, the non-intrusive measuringtechnique that can collect data for an area is preferred.Figure 1.2: Schematic show for stages of the heat transfer of gas-liquid plug flow in a slitmicrochannel in [3]. For plug flow: (I) the thermal entrance region. (II) the transitionregion. (III) the fully developed region. For single-phase flow there is no transition region.Reprinted with permission.61.3. Numerical/experimental studies of plug flow in microchannelA high-resolution infrared(IR) camera was used in [41] to collect the continuous tem-perature distribution of a stainless steel tube which contained plug flow inside. They foundout the developing stages of the plug flow heat transfer with constant input heat flux canbe summarized into following (also see the schematic show in fig. 1.2):• In the entry region (about one plug length long), the Nusselt number for short or longplugs reaches the plug fully developed asymptotic limit, or the continuous flow limitrespectively;• In the transition region (about one period of internal circulation of the plug), Nusseltnumber oscillates for short plugs (validated in [3, 42]);• In the fully developed stage, Nusselt number gradually becomes stable.Despite their explicit reveal of the stages of plug flow heat transfer, the IR camera canonly measure the surface temperature of the outer channel wall, the information inside forboth flow and thermal field are missing.The most precise optical measuring technique of the flow field to date is the microparticle image velocimetry (µPIV), and one can find literature reviews of visualizationin microfluidics in the work of Wereley and Meinhart[43], Sinton[44]. The laser-inducedfluorescence (LIF) method is ideal to measure the temperature field inside the flow fieldnon-intrusively, which was adopted in the work of Ross, Gaitan and Locascio[45], Ghaini,Mescher, and Agar[1]. A third order polynomial correlation was set up to show the tem-perature distribution by the dye particle density in the microchannel, and the precision isabout 2.5 ∼ 3.4◦C [45]. However, these optical methods meet challenges when the light tiltsthrough the interfaces with different refractive indexes, which happens for curved channelwalls [46]. The materials have to be chosen carefully to match the refraction coefficientsin order to avoid the failure. To our best knowledge, there are very few investigations intoplug flow heat transfer using both µPIV and LIF, and it could be a research point in thefuture.A review for numerical work for the plug flow heat transfer are included in [47]. Moststudies had focused on the plug flow heat transfer inside of circular [48, 49] or slit channels71.3. Numerical/experimental studies of plug flow in microchannel[50–52], which can be a result of dramatically increasing the computational time when thefield is unable to be simplified using symmetricity. Though, some works conducted in the3-D domain also use techniques (such as the volume of fluid method, lattice Boltzmannmethod, level set method and moving grid method) to capture the interface between fluids[53–56]. Though these methods can capture details with high precision, the calculation timefor them is too long. In the real design/optimization work where additional simulationsbesides the flow field will be conducted (for multiple times before finding out an optimumdesign), it is preferred using existing flow field results [3, 42, 57]or at least by pre-definingthe interface shape [58] to save calculation time.Despite many investigations of both gas-liquid and liquid-liquid plug flow heat transfer,it has been pointed out in the review works recently that significant gaps exist concerningboth the measured values and the correlations for heat transfer coefficient and pressuredrop between different pieces of literature, and there are little agreements among them[18, 47, 59]. The only scenario where some agreement is observed is the gas-liquid plug flowin the circular microchannel [59]. Most of the correlations did not reveal a clear picture forthe physics of heat transfer, and were merely a collection of obtained results, and thus cannot be applied when it exceeds their original conditions.Since there still lacks reliable correlations for the plug flow heat transfer under variousconditions, fast simulation and dense database can be an alternative. Unlike flow in macrosize and under high speed, the Stokes simplification for slow viscous flow is very reliablebecause Reynolds number for plug flow is small (Re 1) in microchannel [38]. Thus, theonly challenge of building up a dense database remains as calculation time. The normalprocedure for simulating fluid flow heat transfer is by solving the flow field and thermal fieldsimultaneously, which will result in a considerable amount of calculation time. Fortunatelyin the microchannel flow scenarios, the temperature difference does not need to be toomuch to reach a considerable temperature gradient owing to the small scale (i.e., the meantemperature difference between the inlet and the outlet is less than 9 K, the differencebetween the wall and the mean value is less than 2 K in [41]), the thermal properties canusually be assumed to be constant. If the Prandtl number (defined below, where k isthermal conductivity of the fluid) is large (Pr 1), the momentum diffusivity (ν) is much81.4. Analytical models of plug flowslarger than thermal diffusivity (α), then the flow and thermal field can be decoupled to savecalculation time [3, 42, 60].Pr =να=cpµk. (1.4)Moreover, the analytical method can be applied in the microchannel with certain kindsof simple geometry to calculate the flow field with some simplifications, which has nearlyno calculation time and therefore are preferred for building up the database. In the nextsection, a review of all existing analytical methods of plug flow in microchannels is prepared.1.4 Analytical models of plug flowsThe Graetz problem describes the steady-state heat transfer in fully developed internalflow, which implies the zero or the constant temperature derivative in the flow directiondepending on boundary conditions. It applies the assumption of the constant thermalproperties such as ρ, µ, k and cp and incompressible Newtonian flow. This model has beenapplied in studying plug flow heat transfer, and the literature related will be presentedbelow.In the work of Shojaeian and Kosar [61, 62], a simple analytical model was adopted basedon the Graetz problem while the velocity in the channel is assumed to be constant and uni-form. Muzychka, Walsh, and Walsh [63] also made full of the Graetz problem by assumingthe velocity is constant and uniform. They managed to conclude the influence of tunnelgeometries by extracting a few key parameters of the cross-section shape. These methodscould help estimate the heat transfer performance roughly. However, the assumption ofuniform and constant velocity field drops the recirculation inside the plug and can highlyunderestimate the mixing ability of plug flow. A better way is to consider the vortexesformed inside plugs by modeling and solving the flow field in the 2-D cross-section.In order to analytically obtain the vortex inside plugs, the Stokes flow assumption,which implies Reynolds number is very small (Re 1) and drops the convection of themomentum, was firstly used in gas-liquid plug flows by Sivasamy, Che, Wong, Nguyen, andYobas [38]. In their work, the 2-D flow field was built and analytically solved using a fourth91.5. Motivation and Objectiveorder partial differential equation set (PDEs). This model has been proved to be effectiveand convenient in microchannels with different shapes, such as in a slit microchannel [38],in a curved microchannel [64], and in a microchannel with circular cross-section [2]. Somenumerical studies of heat transfer based upon the analytical flow field were carried out[3, 42], and the flow field results helped to understand the mechanism behind the heattransfer enhancement. The analytical model for the liquid plug train in a 2-D channel wasproposed in [57]. The basic idea is the same as the gas-liquid plug flow in [38], thoughthe friction force between two immiscible liquids increased the complexity on the interface.One of the simplified results of the liquid plug train, the liquid-liquid plug flow, was appliedfor the study of heat transfer performance for liquid-liquid plugs at asymmetric boundaryconditions in [60].1.5 Motivation and ObjectiveThis thesis analytically and numerically studies the plug flow and heat transfer of thisflow type.With the fast developing of applications such as MEMS, the size of these devices shrinksdramatically, and the heat flux also increases considerably. Traditional flow cooling meth-ods, which usually generate vortexes in turbulent flow regime, have faced the challengebrought by small characteristic length and small Reynolds number under micro scales. Theplug flow shows promising potential for enhancing heat transfer because it can generate well-defined vortexes even in the Stokes regime (Re 1) resulting from the interface betweenanother kind of immiscible fluid. Hence, plug flow based heat exchanger can be a solutionto the future thermal management in MEMS. However, according to the previous section,there still lacks the analytical studies of plug flow in microchannels with many other kindsof geometry. For instance, one of the most commonly seen geometry in the traditional heatexchanger, the concentric tube, is not studied before. Moreover, in the previous studies,heat transfer simulations were mostly conducted using the transient method which leads toa quite long calculation time even provided with pre-stored flow field results and makes ithard and costly to build/extend the database.101.6. Organization of the ThesisThus, the most primary objective of this thesis is to find the analytical solutions of plugflow in different geometries that people have not studied. The analytical flow field resultscan help understand the mechanism of heat transfer enhancement. The second objectiveis to systematically investigate the heat transfer within plug flows in these geometries, aswell as the influences of geometrical parameters or other parameters upon the heat transferperformance. At last, optimally a simplified heat transfer model can be extracted to savecalculation time so that I can build up a database/an empirical correlation to cover amassive amount of working conditions in an acceptable amount of time. The pre-storeddatabase/pre-defined correlation can help the optimization without the need of runningsimulations for multiple times by the designers. Only a few times of accurate simulationare needed for calibration after they have found the optimal designs.1.6 Organization of the ThesisIn chapter 2, I find the analytical solution of gas-liquid plug flow inside the concentricmicrochannel. The flow patterns like locations of vortex center, the stream functions ofvortex centers, the radial transport velocity and recirculation period are investigated. Thepumping power in the form of skin friction coefficient is studied.In chapter 3, I carry out the heat transfer simulations in plug flow in chapter 2 atthree different kinds of commonly seen asymmetric boundary conditions: The inner iso-flux boundaries (shorten as IFB), the iso-thermal boundaries (shorten as ITB), and theouter iso-flux boundaries (shorten as OFB). Simplified thermal networks for these heattransfer processes are extracted. The influences of plug aspect ratio, the inner-outer radiusratio as well as Peclect number upon heat transfer enhancement are analyzed. Finally, inthis chapter, a simplified model for heat transfer at the fully developed stage is found tosave calculation time. An extensive database containing results from about 12,000 workingconditions are built up for future design work.In chapter 4, I find the analytical solution of liquid-liquid plug flow inside microchannelwith a round cross-section by modeling and solving the flow field in two immiscible fluidssimultaneously. The focus is on flow patterns, among which stands out the secondary vortex111.6. Organization of the Thesisinside the plug with low viscosity when the viscosity ratio is far from 1. The secondary vortexat the cap of the plug is resulting from the momentum transfer due to interface continuity.The influence of plug lengths and viscosity ratio upon the skin friction coefficients areanalyzed in preparation for calculating the pressure loss in real design work.In chapter 5, a summary of this thesis is made. The limitations of this work as wellas those of the analytical modeling are also presented. Some suggestions for the potentialwork are provided in the end.In appendix A, my supervisor Dr. Ri Li and I find out a correlation for calculatingmaximum and minimum stream functions in the plug flow in chapter 2. However, thecorrelations for heat transfer performance or enhancement is not established due to limitedtime. Hence, this part is recorded in an appendix as a secondary outcome of the thesis.12Chapter 2Analytical study of gas-liquid plugflow in concentric microchannel2.1 Mathematic modeling2.1.1 Governing equations and boundary conditionsFigure 2.1: Schematic show for gas-liquid plug flow in concentric microchannelThe plug flow in the micro concentric tube is modeled under a cylindrical coordinatesystem, which is moving at the same speed as the liquid plug, using a fourth order PDE. Thegoverning equations have the following assumptions: 1. the flow field is fully developed, 2.the liquid plug takes up the whole cross-section, 3. the interface between the liquid and thegas is flat, 4. the whole flow field is rotationally symmetric, 5. the fluid is Newtonian withuniform and isotropic physical properties, 6. the fluid obeys Stokes assumption (Re 1).Before starting the modeling procedures, I would like to talk firstly about the gap132.1. Mathematic modelingbetween these assumptions and real scenario. To my best knowledge, the Stokes assumptionis the most appropriate method to date for analytically calculating the flow field in plugflows. The request of low Reynolds number (Re 1) can be easily reached under micro-scales. Meanwhile, it can capture the vortexes formed because the flow field is solved in the2-D domain.However, some assumptions above such as that the whole cross-section is taken up bythe plug (no liquid film) can hardly be reached in real scenarios. It has been studied in[41, 51, 52, 54, 65] that the liquid film’s thickness is very small and negligible only whenCa 1 (capillary number Ca is defined below, the ratio of viscous force and the surfacetension). Otherwise, the liquid film will be thicker and the plug shape will deform owing tothe viscous force. The existence of liquid film can also lead to the slippery boundary betweenthe plug boundary and the wall, which makes the modeling and deriving the solution morecomplicated.Ca =µUλ. (2.1)The interface is also hard to be flat in the real cases. The influence of curvature uponthe flow field was investigated in [2], where the author set the contact angels at two ends tobe the same. It was found out that the flow field was not much affected when the contactangel is between 45◦ ∼ 135◦. Under these cases the flat end assumption is valid and thecontribution of surface tension upon, for instance, the pressure loss, can be de-coupled fromthe contributions from the frictional force. However, in real cases, the contact angle can bedifferent at two ends, and it is also different at two walls if the channel is not symmetric(i.e., the meandering channel or the concentric channel in this work) which leads to toomany combinations to validate.As mentioned earlier in section 1.3, the capture of these detailed features can only bedone by accurate multi-phase flow simulation with surface capturing methods if they arenot pre-defined, which are also very time-consuming. Negotiation has to be made betweenthe accuracy and efficiency. Hence, in this work, I will still adopt the analytical methodwith the Stokes flow model and the above assumptions because I aim to build either aquick method or a dense database for the design need. Though it is recommended that in142.1. Mathematic modelingthe real design work, the analytical results in this chapter, the results from the simplifiedmodel of heat transfer in chapter 3 and the analytical results in chapter 4 can be used forprimary searching to narrow down the range of optimal designs quickly. Then a few timesof accurate simulation can be conducted to calibrate the influences brought by liquid filmor curved interfaces.Following the assumptions above, the continuity equation is,∇ • u = 0. (2.2)In the momentum equation, the convection term is neglected (Re 1 ),µ∇2u−∇P = 0. (2.3)By applying the Stokes stream function, the continuity equation eq. (2.2) is satisfiedautomatically, and the momentum equation eq. (2.3) is simplified as below by taking curveat both sidesL4−1ψ = 0, (2.4)where the Stokes stream function ψ is defined asur = −1r∂ψ∂z, uz =1r∂ψ∂r, (2.5)and the operator L4−1 isL4−1 = (∂2∂z2+∂2∂r2− 1r∂r)(∂2∂z2+∂2∂r2− 1r∂r). (2.6)Three assumptions are made here to obtain the boundary conditions: First, because ofthe continuity of the streamline, the stream function can be set to a constant value at allboundaries. For simplicity, it is set to zero. Second, there is no slip between the walls andthe liquid plug. Thus, the axial velocity at walls is −U in the moving reference. At last,there is no shear stress caused by the gas because usually, gas has low viscosity. As a result,the boundary condition can be written as,152.1. Mathematic modelingψ(0, r) = 0, ψ(l, r) = 0, ψ(z, r1) = 0, ψ(z, r2) = 0, (2.7)1r∂2ψ∂z2(0, r) = 0,1r∂2ψ∂z2(l, r) = 0, (2.8)1r∂ψ∂r(z, r1) = −U, 1r∂ψ∂r(z, r2) = −U. (2.9)The nondimensionalization is done before solving the PDEs, and the aspect ratio β, theinner-outer radius ratio η and the eigenvalues ωn are also defined byzˆ ≡ zr1, rˆ ≡ rr1, uˆz ≡ uzU, uˆr ≡ urU, ψˆ ≡ ψUr21, (2.10)η =r2r1= rˆ2, β =lˆ1− η , ωn =npilˆ, (2.11)Substitute eqs. (2.10) and (2.11) into eqs. (2.4) and (2.7) to (2.9) to obtain the dimen-sionless governing equation and boundary conditionsLˆ4−1ψˆ = 0, (2.12)ψˆ(0, rˆ) = 0, ψˆ(lˆ, rˆ) = 0, ψˆ(zˆ, 1) = 0, ψˆ(zˆ, η) = 0, (2.13)1rˆ∂2ψˆ∂zˆ2(0, rˆ) = 0,1rˆ∂2ψˆ∂zˆ2(lˆ, rˆ) = 0, (2.14)1rˆ∂ψˆ∂rˆ(zˆ, 1) = −1, 1rˆ∂ψˆ∂rˆ(zˆ, η) = −1. (2.15)2.1.2 Analytical solutionStrong periodicity can be observed focusing the boundary conditions at the two plugends in eqs. (2.13) and (2.14) Thus, the finite Fourier transformation can be applied here.162.1. Mathematic modelingDue to the zero value of the function at these ends, the Sine transformation is chosenSn[ψˆ] = 2lˆ∫ lˆ0ψˆ sin (npilˆzˆ)dzˆ = gn(rˆ). (2.16)The anti-transformation can be applied to obtain the original function, which is theuniversal solution to the PDEs.ψˆ =∞∑n=1ψˆn =∞∑n=1sin (npilˆzˆ)gn(rˆ). (2.17)The properties of the Sine transformation are listed below.Sn[∂2ψˆ∂zˆ2] = −ω2ngn +2ωnlˆ[ψˆ(0, rˆ)− (−1)nψˆ(lˆ, rˆ)] = −ω2ngn, (2.18)Sn[∂4ψˆ∂zˆ4] = −ω2nSn[∂2ψˆ∂zˆ2] +2ωnlˆ[∂2ψˆ∂zˆ2(0, rˆ)− (−1)n∂2ψˆ∂zˆ2(lˆ, rˆ)] = ω4ngn, (2.19)Sn[∂ψˆ∂rˆ] =dgndr. (2.20)Substitue eqs. (2.16) and (2.18) to (2.20) into eqs. (2.12) to (2.15), the system thenbecomes a 4th order, linear and homogeneous ordinary differential equation one (ODEs).(d2dr− 1rddr− ω2n)(d2dr− 1rddr− ω2n)gn = 0, (2.21)gn(η) = 0, gn(1) = 0, (2.22)g′n(η) = −2η[1− (−1)n]lˆωn, g′n(1) = −2[1− (−1)n]lˆωn, (2.23)The universal solution of the ODEs eq. (2.21) isgn(rˆ) = Anrˆ2I2(ωnrˆ) +BnrˆI1(ωnrˆ) + Cnrˆ2K2(ωnrˆ) +DnrˆK1(ωnrˆ), (2.24)where Iν ,Kν are the νth order of first and second modified Bessel functions respectively.172.1. Mathematic modelingAlso, obviously, for series n being even numbers, only zero solution will be obtained. In thischapter, the series n is by default set to be odds without further notice. For odd series, theconstant coefficients can be mounted using the boundary condition eqs. (2.22) and (2.23).The system is formed with n linear and homogeneous subsystems each with a dimension of4× 4, which can all be solved easily by Cramer’s rule [66].I2(ωn) I1(ωn) K2(ωn) K1(ωn)ηI2(ωnη) I1(ωnη) ηK2(ωnη) K1(ωnη)I1(ωn) I0(ωn) −K1(ωn) −K0(ωn)ηI1(ωnη) I0(ωnη) −ηK1(ωnη) −K0(ωnη)AnBnCnDn =00− 4lˆω2n− 4lˆω2n. (2.25)A series solution of the eqs. (2.12) to (2.15) can be traced back finally by puttingeq. (2.24) into eq. (2.17). Then, the velocity field is obtained using the definition of thestream function in eq. (2.5).ψˆ =∑n=1,3,5...sin(ωnzˆ)[Anrˆ2I2(ωnrˆ) +BnrˆI1(ωnrˆ) + Cnrˆ2K2(ωnrˆ) +DnrˆK1(ωnrˆ)], (2.26)uˆz =∑n=1,3,5...ωn sin(ωnzˆ)[AnrˆI1(ωnrˆ) +BnI0(ωnrˆ)− CnrˆK1(ωnrˆ)−DnK0(ωnrˆ)], (2.27)uˆr = −∑n=1,3,5...ωn cos(ωnzˆ)[AnrˆI2(ωnrˆ) +BnI1(ωnrˆ) + CnrˆK2(ωnrˆ) +DnK1(ωnrˆ)]. (2.28)As mentioned earlier, the pumping power is a key factor in evaluating the performanceand efficiency of a heat exchanger. Thus, the pressure difference ∆P and the skin frictioncoefficient Cf are derived in the following for later study.To be noted that there are too many factors affecting the skin friction coefficient Cf to beall included (like the capillary number Ca, Weber number We and the interface curvature).182.2. Results and discussionsHere, only the friction force caused by the radial gradient of the internal recirculation isconsidered. The influences of these interface geometry parameters are neglected due to theassumption of the two flat ends.Ff = 2piµ(r∫ l0∂uz∂rdz)∣∣∣∣r1r2, (2.29)∆P =Ffpi(r21 − r22). (2.30)The Fanning friction factor and the skin friction coefficient are defined as below, whereDh =4pi(r21−r22)2pi(r1+r2)= 2(r1 − r2),f =14Dh∆Pl12ρV2=Dh∆P2ρlV 2, (2.31)Cf = fRe. (2.32)eqs. (1.1) and (2.29) to (2.32) yield the result for the Fanning skin friction coefficientCf =8β(1 + η)∑n=1,3,5...ωn[Anrˆ2I0(ωnrˆ) +BnrˆI1(ωnrˆ) + Cnrˆ2K0(ωnrˆ) +DnrˆK1(ωnrˆ)]∣∣∣∣η1.(2.33)The sensitivity study is carried out, where the skin friction coefficient Cf is used as theobjective. It is clearly shown in fig. 2.2 that Cf becomes stable after the series number ishigher than 101. For security, for all of the calculations, I set the series number to 101 toensure accuracy.2.2 Results and discussions2.2.1 Verify solution accuracy using CFD ResultsA comparison to the numerical simulation conducted in the Fluent 16.2 (ANSYS, inc.USA) is made to verify the accuracy of the analytical results. The inner-outer radius ratio192.2. Results and discussionsFigure 2.2: The sensitivity study for the analytical solution. Variation of Cf with thegrowing series number. β = 2, η = 0.5.is η =, 0.5, 0.75 , and the aspect ratio is β = 1. The frontal and rear ends of the plug areset to be planar to fit the assumption of the analytical model. The laminar flow model isapplied, the Reynolds number is set to be 0.5. Water is chosen as the operating fluid. Thenumber of meshes for all three cases are about 5×105, which passes the mesh independencecheck. The outer and the inner walls are set to moving wall condition with −10mm/s. Forthe two ends between the liquid plug and the gas, no shear wall condition is applied. Theabsolute criteria for convergence are set to be 10−3 for continuity and velocity in all threedirections. Since the model is axisymmetric, I randomly choose a slice which passes theaxial. Then the data for the flow field on these slices are exported for comparison.From fig. 2.3.(a) and .(b), the velocity field of the analytical model and that of the CFDsimulation corresponds well with each other. No obvious difference could be found. Thecomparison here can be treated as a verification to the analytical solution’s accuracy.202.2. Results and discussionsFigure 2.3: The velocity vector plot of (a). the analytical model and (b). the CFD results,β = 1, η = 0.502.2.2 Validate using continuous flow(1-D model)To validate the analytical model using continuous flow, the flow field of fully developedStokes flows in the concentric annulus was calculated in [67],uˆ1D(rˆ) =2(rˆ2 − 1) ln η − 2(η2 − 1) ln rˆη2 − 1− (η2 + 1) ln η − 1. (2.34)The stream function can be integrated from eq. (2.34) and the constant term could bebounded using the zero value boundary condition.ψˆ1D(1) = ψˆ1D(η) = 0, (2.35)ψˆ1D(rˆ) =[(rˆ2 + η2 − 1) ln η − (2η2 − 1) ln rˆ]rˆ2 + rˆ2 ln rˆ − η2 ln η2(η2 − 1− (η2 + 1) ln η) . (2.36)The skin friction coefficient for 1-D model is also be derived in a similar fashion (eqs. (1.1),(2.29) to (2.32) and (2.34))Cf,1D =16(1− η)2 ln ηη2 − 1− (η2 + 1) ln η . (2.37)Especially, the skin friction coefficient can be calculated under 2 asymptotic situationswhere the inner radius is zero (in a circular tube, η → 0, Cf,1D → 16) and is infinitely close212.2. Results and discussionsto the outer radius (in a slit channel, η → 1, Cf,1D → 24), respectively. The results for Cf,1Dwill be used as a validation between continuous flow and very long plugs in section 2.2.6.Figure 2.4: The comparison between plugs with different aspect ratio β and the 1-D modelof (a) the axial velocity uˆz. (b) The stream function ψˆ. The sample line is at zˆ = 0.25 lˆ,and the inner-outer radius ratio is η = 0.90.The comparison between the axial velocity and the stream function on the sample line(zˆ = 0.25 lˆ) of the plug flow model and these of the continuous 1-D model is plotted infig. 2.4. When the plug is short (i.e., β = 1.0), a maximum difference of 0.2 and 0.0025 canbe observed for uˆz and ψˆ, respectively. When β = 2.0, the difference becomes very small butstill visible. When the plug is long enough (i.e., β = 10.0), the axial velocity distribution,as well as that of stream function, will be infinitely close to those of the 1-D model. Thus,the result here can be considered as a validation under the asymptotic scenarios.2.2.3 Two asymmetric vortexesThe vortexes in plug flow are forced vortexes resulting from the driving force of themoving wall (in the plug reference) and the interface between fluids. Hence, the flow patterninside the plug can be influenced by the geometry parameters including the inner-outerradius ratio η of the channel and the aspect ratio β of the plug.A group of contours of the stream function is presented in fig. 2.5. Two asymmetricvortexes can be observed. The one closer to the inner wall is referred to as the inner vortex,while the other is called the outer vortex. Obvious variation such as the size of the two222.2. Results and discussionsvortexes, the location of their centers and the maximum/minimum stream function valuescan be observed with growing η. Thus, discussions upon these phenomena will be conductedin this subsection.Figure 2.5: The contour plots for stream function in plugs with different η, the cross shapemarker are the locations for vortex centers.Since the inner and the outer vortex have different orientations, they can be classifiedusing the signal of their stream functions. Due to the rotationally symmetric system, I canset the polar angel interval to be 0 ∼ 1 for simplification. Then, the size or the volumetaken up by the inner vortex can be numerically integrated using the formula below, wheresgn is the signal function,Vˆin,v =∫ lˆ0∫ 1η1− sgn(ψˆ)2rˆdrˆdzˆ. (2.38)Similarly, the volume of the outer vortex can be written as,Vˆout,v =∫ lˆ0∫ 1η1 + sgn(ψˆ)2rˆdrˆdzˆ. (2.39)The total volume of the plug is,Vˆplug = Vˆin,v + Vˆout,v =∫ lˆ0∫ 1ηrˆdrˆdzˆ =(1− η2)lˆ2. (2.40)The volume ratio between the inner and the outer vortex Vˆin,v/Vˆout,v is plotted in fig. 2.6.Vˆin,v/Vˆout,v is always smaller than 1 and it is a increasing function of η. When η → 1 twovortexes become nearly identical to each other, the volume ratio Vˆin,v/Vˆout,v becomes closeto 1, which is like the plug flow in a slit channel. Vˆin,v/Vˆout,v is also a decreasing function232.2. Results and discussionsof the aspect ratio β, especially when it is short (β < 2). This indicates the inner vortexcan take up more space when β < 2. When β > 2, volume ratio becomes nearly stable andthe influence of β is not obvious.Figure 2.6: The volume ratio between the inner and the outer vortex Vˆin,v/Vˆout,v.The volume change of vortexes can also lead to the swift locations of their centers. Thecenters always locate at the middle (axial direction) of the plug, and the radial index arewhere stream function reaches maximum/minimum value (rˆψmax , rˆψmin). Define the indexfor the relative locations in the radial direction of vortex centers as below,Iout,v =rˆψmax − η1− η , Iin,v =rˆψmin − η1− η . (2.41)The indexes Iout,v and Iin,v are plotted in fig. 2.7. When η increases, both Iout,v andIin,v increases, which means both centers move towards the outer wall. As shown in fig. 2.6,the inner vortex takes more volume of the whole and expands when η increases, it is naturalfor the inner vortex center to shift away from the inner wall. Meanwhile, the outer vortexshrinks and the center should move towards the outer wall for the similar reason. Moreover,the indexes also become independent of the aspect ratio when β > 2.242.2. Results and discussionsFigure 2.7: The index for relative location in radial direction of vortex centers.The strength difference between the two vortexes is also of interest. Thus, the absolutevalue of the ratios of minimum and maximum stream functions |ψˆmin/ψˆmax| are plotted infig. 2.8. |ψˆmin/ψˆmax| has nearly the same variations to those of the volume ratio Vˆin,v/Vˆout,v(in fig. 2.6). Thus there is no need to describe these variations again. However, the ratiosof both volumes and the stream functions only reveals the difference between the innerand the outer vortex in the same plug, and it can not be used to describe the magnitudesof circulation cross multiple cases with different geometry parameters. The magnitudes ofcirculation inside the plug thus will be discussed in the next subsection.2.2.4 Quantify radial transport using the averaged recirculation fluxChe, Wong, and Nguyen [3] have pointed out in their research that the heat transferenhancement inside the plug flow (in a slit channel) is due to the internal recirculationand the transverse flow near the front and the rear ends. They plotted out the transversevelocity distribution on the sample lines which cross the vortex centers for investigation.Here I adopt the similar strategy by putting sample lines through the two vortex centers.However, instead of using the radial velocity, the product of velocity and radius rˆuˆr (referredto as the recirculation flux) is studied owing to the changing cross-section area at different rˆ252.2. Results and discussionsFigure 2.8: The stream function ratio between the inner and the outer vortex centersψˆmin/ψˆmax.in the cylindrical coordinator. Moreover, the averaged value of rˆuˆr at inner/outer vortexesare studied instead of their distributions for the convenient comparison between multiplecases. The averaged recirculation flux can be calculated by concerning the volumetric flowrate (the change of stream function defined by eq. (2.5)) on the right half of sample lines,ψˆmax/min = ψˆmax/min − 0 =∫ lˆ/2lˆ(−rˆuˆr)dlˆ∣∣∣∣rˆψˆmax/min=∫ lˆlˆ/2rˆuˆrdlˆ∣∣∣∣rˆψˆmax/min.Reorganize the equation above and define the averaged recirculation flux (rˆ ˆ¯ur)out,v/in,vof two vortexes, where the signal of it represents the direction of rotation only.(rˆ ˆ¯ur)out,v/in,v =ψˆmax/minlˆ − lˆ/2 =2ψˆmax/minlˆ. (2.42)As shown in fig. 2.9 recirculation flux at both vortexes are decreasing functions of aspectratio β, which is referred to as the diminishing effect when the plug length increases. In thesmall subplot in fig. 2.9, the recirculation at the outer always vortex is stronger than that of262.2. Results and discussionsthe inner vortex. When the inner-outer radius ratio η increases, the averaged recirculationflux of the inner vortexes becomes stronger, and that of the outer vortex becomes weaker.When η → 1, they become nearly the same in absolute value because the channel tends tobe a slit one and the flow field becomes symmetric.Overall, the averaged recirculation flux (rˆ ˆ¯ur)out,v/in,v can describe and compare themagnitudes of vortexes under multiple cases. This parameter can help understand the heattransfer enhancement.2.2.5 Recirculation periodAs mentioned earlier in section 1.3, the recirculation period is an important factor forevaluate procedures of heat transfer in plug flow such as estimating the entrance region.However, the definition of the recirculation period is not consistent. For example, it wasroughly defined as one plug length plus twice of the hydraulic diameter of the tube (lplug +2Dh) in [41] because this definition can suit his explanation for experiment results, whilemore accurate trajectory simulations were conducted in [3, 42] to record the distribution ofrecirculation period. Here I applied the latter way to investigate the recirculation period.The recirculation periods are recorded by numerically tracing some passive points whichare initially put on the central plane until they finish one period of recirculation. Thenondimensionalization is conducted as below,tˆ =tU(r1 − r2) . (2.43)The distributions of recirculation period tˆ on the central plane(in axial direction, zˆ = lˆ/2)are plotted in fig. 2.10. When β increases, the amplitude of tˆrec increases exponentiallyowing to the diminishing effect mentioned earlier, which can lead to weaker heat transferenhancement. The double-arrow marks are the boundaries between the inner and the outervortexes, and it also moves slightly towards the inner wall, which indicates a smaller volumeratio Vˆin,v/Vˆout,v. This matches well with the conclusion made in fig. 2.6.272.2. Results and discussionsFigure 2.9: The recirculation flux rˆuˆr of the inner vortex and the outer vortex. Main plotis the averaged value, the small subplots show where the sample line is (right-top), and adetailed distribution on the sample line for β = 2, η = 0.5 (right-bottom), respectively.282.2. Results and discussionsFigure 2.10: The distribution of recirculation time tˆ on the central of the plug (zˆ = lˆ/2).The double-arrows mark the boundaries between the inner and the outer vortexes. β =1, 2, 3, 4, η = 0.5.2.2.6 The skin friction coefficientThe influences of β are presented in fig. 2.11.(a). under selected η. With the increasingβ, Cf drops down dramatically owing to the diminishing effect of the recirculation. The dashlines in fig. 2.11.(a). are the results for the 1-D model calculated through eq. (2.37). The1-D results are infinitely close to those of long plugs and it can be treated as a validation.When β > 46 the gap between the Cf of the plug flow and the 1-D model is lower than1%. Thus, Cf can be simply evaluated using the results from 1-D model (eq. (2.37)) whenthe plug is long enough (β > 46). An asymptotic case with η = 10−5 is also conducted tocompare with the Cf in plug flow with circular cross-sections [2]. And, the results show aperfect agreement.The influences of η are presented in fig. 2.11.(b). under selected β. Overall, Cf increasesas η grows. The trend is more obvious for relatively longer plugs. For instance, Cf increasesabout 12.5%, 25% and 50% for β = 2, 10, 100 when η increases from 0.05 to 0.95, respec-tively. For long plugs the internal recirculation is quite weak, thus the influence of different292.3. Summary of this chaptercurvature of the 2 walls plays an important role. While for short plugs (β = 1) the flow isalmost characterized by the recirculation, thus Cf is less dependent on η.Figure 2.11: The skin friction coefficient (a). versus β for selected η and (b). versus η forselected β2.3 Summary of this chapterThe plug flow has two well-defined vortexes due to the immersible interfaces betweenthe liquid and the gas, which can potentially enhance the mixing ability and heat transfer.In this work, the plug flow in the concentric microchannel is modeled using a 4th orderPDE set. The series solution for the set is found, and the flow field is investigated indetail. Focuses are made upon the influence of the geometry parameters (the aspect ratioβ, and the inner-outer radius ratio η) upon the flow pattern such as vortex center locationI, intensity ratio |ψˆmin/ψˆmax|, the volume ratio Vˆin,v/Vˆout,v, the averaged recirculation fluxrˆ ˆ¯ur of two vortexes and the recirculation period tˆ. The influences upon the skin frictioncoefficient Cf are also studied. These findings are summarized as below:• The flow pattern is greatly influenced by the inner-outer radius ratio η. The higherη makes the flow tend to be like plug flow in a 2-D microchannel with symmetricpatterns (Vˆin,v/Vˆout,v → 1, |ψˆmin/ψˆmax| → 1). At lower η, however, two differentvortexes were generated. The outer vortex is similar to that in the microchannel withround cross-section which takes up most of the volume in the plug, while the inner302.3. Summary of this chaptervortex is weaker and takes up the smaller volume (Vˆin,v/Vˆout,v < 1, |ψˆmin/ψˆmax| < 1).• The radial transport phenomena can be quantified using the averaged recirculationflux. At the inner vortex, the averaged recirculation flux grows stronger with increas-ing η, while that at outer vortex becomes weaker. The averaged recirculation fluxesof both vortexes decrease dramatically when β increases, which leads to exponentiallyincreasing recirculation period with increasing β. This is referred to as the diminishingeffect of recirculation with increasing β.• The skin friction coefficient Cf can be influenced by both β and η. Overall Cf dropsnearly exponentially with growing β. It increases with growing η, though this trend ismore obvious for relatively longer plugs. For instance, Cf increases about 12.5%, 25%and 50% for β = 2, 10, 100 respectively when η increases from 0.05 to 0.95.31Chapter 3Numerical study of heat transferinside plug flow in chapter 2 atasymmetric boundary conditions3.1 Modeling for transient heat transfer process in plugflow3.1.1 Governing equation and initial/boundary conditionsThe governing equation for transient heat transfer for a single plug can be written asbelow, where α = k/ρcp represents the thermal diffusivity, k is the conductivity and cp isthe specific heat of the fluid.∂T∂t+1r∂(urTr)∂r+∂(uzT )∂z= α[1r∂∂r(r∂T∂r) +∂2T∂z2]. (3.1)It should be noted that the influences of viscous dissipation, body forces, and pressuredifference are all neglected by assuming the Eckert number is small (Ec 1). The physicalmeaning of Ec is the ratio of the kinetic energy to the enthalpy driving force (cp∆T )for heat transfer (∆T is the temperature difference between the wall and the mean valueof fluid). This assumption is reasonable and can be validated by calculating Ec usingparameters from previous studies. For example, in the experimental work of Asthana,Zinovik, Weinmueller and Poulikakos[22] the velocity is U = 1 m/s, the working liquid iswater with cp = 4186J/(K · kg), ∆T ∼= 2.55K. Thus the Eckert number is Ec = 0.9E−5 1and the effects of above can be neglected.323.1. Modeling for transient heat transfer process in plug flowEc =U2cp∆T. (3.2)In the real cases, the initial condition should be a steady-state temperature distribution,which is determined by the boundaries before the heating part of the channel. Though usingthe real initial condition can help improve the accuracy of the developing process of thefield, it requires much more input pieces of information and adds more complexity to ourinvestigation. Hence for the initial condition in this work, the whole plug is set to have auniform temperature distribution.T = Ti. (3.3)There are many kinds of boundary conditions in reality, and it is costly to evaluatethe heat transfer performance at all of these boundary conditions. Hence, three typicalboundary types are chosen to study in this thesis and they are listed as below.• The inner iso-flux boundaries (shorten as IFB). It represents the plug enters a heatexchanger with constant heat flux at the inner wall, while the outer wall is isolated.• The iso-thermal boundaries (shorten as ITB). It represents the plug enters a heatexchanger with constant inner and outer wall temperatures.• The outer iso-flux boundaries (shorten as OFB). This boundary type is basically thesame as IFB, the only difference is the orientation of the heat transfer that constantheat flux is loaded from the outer wall, while the inner wall is isolated.The schematic show for transient heat transfer is plotted in fig. 3.1. As the liquid plugmoves towards the downstream of the channel, it gradually accumulates heat. The thermalfield inside also develops. The development of the thermal field, as well as heat transferability, is of interest.The boundary condition for IFB− k∂T∂r(z, r2) = q′′, −k∂T∂r(z, r1) = 0. (3.4)333.1. Modeling for transient heat transfer process in plug flowFigure 3.1: The schematic show of boundary conditions (a). The inner iso-flux boundaries(IFB). (b). The iso-thermal boundaries (ITB). (c). The outer iso-flux boundaries (OFB).The contours are taken at tˆ = (1, 6, 25, 40) (eq. (2.43)) under these 3 boundary typescorrespondingly. (β, η,Pe) = (2, 0.5, 100).343.1. Modeling for transient heat transfer process in plug flowFor ITBT (z, r2) = Th, T (z, ri) = Ti. (3.5)For OFB− k∂T∂r(z, r2) = 0, −k∂T∂r(z, r1) = −q′′. (3.6)The heat transfer at two ends of the plug is neglected because its segmented by gas.Thus the adiabatic condition at the two ends of the plug is as below,∂T∂z(0, r) = 0,∂T∂z(l, r) = 0. (3.7)3.1.2 NondimensionalizationThe nondimensionalization for governing equations and boundary conditions could beconducted using the parameters belowTˆ ≡ T − Ti∆Tstd, (3.8)where the standard temperature difference ∆Tstd =q′′(r1−r2)k is for the IFB and OFB,∆Tstd = Th − Ti for ITB.Substitute eqs. (2.10), (2.11), (2.43) and (3.8) into eqs. (3.1) and (3.3) to (3.6) to obtainthe the dimensionless governing equations∂Tˆ∂tˆ+ (1− η)( uˆrrˆ∂(Tˆ rˆ)∂rˆ+ uˆz∂Tˆ∂zˆ) =(1− η)2Pe[1rˆ∂∂rˆ(rˆ∂Tˆ∂rˆ) +∂2Tˆ∂zˆ2], (3.9)where Pe is the Peclet number. Pe represents the ratio between advection ability anddiffusion ability,Pe =U(r1 − r2)α. (3.10)The initial conditionTˆ = 0. (3.11)353.2. Modeling for the asymptotic case - continuous Stokes flow heat transferGridsTˆ : vertex centered, uniform gridsuˆ: staggered grids, pre-stored analytical resultsNr, Nz 200, β ∗ 200∆tˆ 18 min(1/Nr,Pe/N2r )Advection terms Linear upwind differential scheme (LUDS)Diffusion terms Central differential scheme (CS)Temporal terms First order forward in time (1st FT)Table 3.1: Set-up for numerical solution of heat transferThe boundary condition for IFB∂Tˆ∂rˆ(zˆ, η) = − 11− η ,∂Tˆ∂rˆ(zˆ, 1) = 0. (3.12)For ITBT (zˆ, η) = 1, T (zˆ, 1) = 0. (3.13)For OFB∂Tˆ∂rˆ(zˆ, η) = 0,∂Tˆ∂rˆ(zˆ, 1) =11− η . (3.14)And, the adiabatic condition at the two ends of the plug∂Tˆ∂zˆ(0, rˆ) = 0,∂Tˆ∂zˆ(lˆ, rˆ) = 0. (3.15)The eqs. (3.9) and (3.11) to (3.15) are discretized using the finite volume method, thencoded to be simulated in MATLAB 2016A. As illustrated in section 1.3, the flow field andthe thermal field is decoupled under the scenarios where Pr 1. The simulation uses thepre-stored flow field results obtained by the analytical method in chapter 2. More detailsabout the setup of the numerical solving can be examined in table 3.1.3.2 Modeling for the asymptotic case - continuous Stokesflow heat transferWhen discussing the heat transfer performance of the plug flow configuration, usuallythat of the single-phase flow is also derived for comparison. For instance, the ratio of the363.2. Modeling for the asymptotic case - continuous Stokes flow heat transferNusselt number Nu of the two kinds of configurations can be used to show the enhancementby plug flow.To obtain the fully developed thermal field of the continuous Stokes flow heat transferinside the concentric tube, firstly a Euler reference sticking to the ground is used for sim-plification. Under this reference, the fully developed thermal field is a function of spatialindexes (zˆ, rˆ) only, and the axial derivative ∂Tˆ∂zˆ becomes constant and independent of theradial index rˆ [68]. The derivation of the solution for the fully developed thermal field atIFB is shown as below, while those at other two boundary types can also be obtained in asimilar fashion.With the help of the control volume method, it is found out that∂Tˆ∂zˆ=2ηPe(1− η2) . (3.16)A simplified version of eq. (3.9) is used to describe the thermal field by substituteeq. (3.16) into it, plus by setting the radial velocity uˆr = 0, the axial diffusion term∂2Tˆ∂zˆ2= 0,and the temporal term ∂Tˆ∂tˆ= 0.2η(1 + η)(1− η)2 (uˆz rˆ) =1rˆ∂∂rˆ(rˆ∂Tˆ∂rˆ). (3.17)To be noted here the uˆz is the dimensionless velocity to the ground, thus it equals touˆ1D(rˆ) + 1 in the Euler reference sticking to the plug in eq. (2.34). For simplification, theform of the velocity field is reorganized for later derivation.uˆz = Erˆ2 + F ln rˆ − E, (3.18)where the constant coefficients E,F areE =2 ln ηη2 − η2 ln η − ln η − 1 , F =2(1− η2)η2 − η2 ln η − ln η − 1 . (3.19)Substitute eqs. (3.18) and (3.19) into eq. (3.17), then integrate twice to obtain the fullydeveloped thermal field in the form of difference between local temperature and that ofthe inner wall Tˆ − Tˆin,w. Since in the internal flow, heat transfer performance is usually373.3. Heat transfer performance at inner iso-flux boundary (IFB)described using dimensionless parameters (like Nusselt number Nu or thermal resistance σ)based upon temperature difference and spatial derivative, the distribution of Tˆ − Tˆin,w isenough for calculation.Tˆ − Tˆin,w = 2η(1 + η)(1− η)2 (E16rˆ4 +F4rˆ2 ln rˆ − E + F4rˆ2 +E + F4ln rˆ)∣∣∣∣rˆη. (3.20)Using a similar fashion, the results under the inner iso-thermal and the outer iso-fluxboundaries can also be obtained.For ITB,Tˆ = logη rˆ. (3.21)For OFB,Tˆ − Tˆout,w = 2(1 + η)(1− η)2 (E16rˆ4 +F4rˆ2 ln rˆ− E + F4rˆ2 +E + F + 2− 2η4ln rˆ)∣∣∣∣rˆ1. (3.22)The calculation was conducted in MATLAB 2016A, and results were stored in thecompatible style as these from plug flows, so they can be quoted for the same post-processing(such as calculating Nu), and then compared with each other. In the following sections,subscript sp represents that the results are taken from single-phase/continuous Stokes flowheat transfer.3.3 Heat transfer performance at inner iso-flux boundary(IFB)In this section, the heat transfer process at IFB is presented versus the dimensionlesstime tˆ. Then a simplified thermal network is introduced, where several dimensionless pa-rameters such as the dimensionless thermal conductance σplug are defined to evaluate theheat transfer performance. The enhancement by applying plug flow configuration is calcu-383.3. Heat transfer performance at inner iso-flux boundary (IFB)lated by comparing to the single-phase flow in the form of conductance ratio σplug/σsp, andinfluences of dimensionless inputs η, β and Pe are presented and discussed.3.3.1 Processes of heat transfer and its simplified thermal networkFigure 3.2: Heat transfer process at IFB, (β, η,Pe) = (2, 0.5, 100). (a). (I) ∼ (IV) Sequenceof dimensionless temperature Tˆ distribution at tˆI∼IV = (1, 6, 25, 40). (b). Variation ofdimensionless temperature at the inner wall Tˆin,w, at the outer wall Tˆout,w and the meanvalue ˆ¯T , and the mixing index γ. The black solid points correspond to the examples shownin (a). (c). Simplified thermal network for IFB.A typical sequence of the developing thermal field of a plug (β, η,Pe) = (2, 0.5, 100)at IFB is plotted in fig. 3.2. After entering the heating channel with a uniform initialtemperature, the thermal boundary layer first develops within the inner vortex with the393.3. Heat transfer performance at inner iso-flux boundary (IFB)help of both advection and diffusion (fig. 3.2 (a). (I).). The mean temperature ˆ¯T (definedbelow) and the inner wall average temperature ˆ¯Tin,w both increase (fig. 3.2 (b).). At thisbeginning stage, there is nearly no heat input into the outer vortex yet, thus no obviousincrease is observed for the outer wall average temperatureˆˆTout,w.ˆ¯T =∫ r1r2∫ l0 ruzTdzdr∫ r1r2∫ l0 ruzdzdr=∫ 1η∫ lˆ0 rˆuˆzTˆdzˆdrˆ∫ 1η∫ lˆ0 rˆuˆzdzˆdrˆ. (3.23)The expansion of the thermal layer continues and reaches the boundary between theinner and the outer vortexes (fig. 3.2 (a). (II).). Then there occurs an obvious heat transferbetween the inner and the outer vortex. The heat transferred through the vortex boundaryis considered as the input into the outer vortex, and it initiates the increase of the outerwall average temperature ˆ¯Tout,w (fig. 3.2 (b).) with the help of the advection and diffusionof the outer vortex.After developing with the help of the convection in the inner vortex, the heat transferthrough the vortex boundary and the convection in the outer vortex, the relative shapeof the temperature contour gradually becomes steady (fig. 3.2 (a). (III) and (IV).). Thisyields the steady and identical growth rate of ˆ¯T , ˆ¯Tin,w andˆ¯Tout,w. The mixing index γ isdefined to describe the relative location of ˆ¯T between the interval of ˆ¯Tin,w andˆ¯Tout,w. Atthis stage, the mixing index γ becomes steady due to the identical growth rate of ˆ¯T , ˆ¯Tin,wand ˆ¯Tout,w. Thus, the fully thermal developed field is considered to be reached here.γ =ˆ¯T − ˆ¯Tout,wˆ¯Tin,w − ˆ¯Tout,w. (3.24)The result in fig. 3.2 reveals the process of heat transfer at IFB: (1). The convectioninside the inner vortex. (2). The development of the heat transfer through the boundarybetween the inner and the outer vortexes. (3). The convection in both vortexes and (4).The fully thermal developed field. Under the fully thermal developed field, a simplifiedthermal network is proposed in fig. 3.2 (c). to describe the heat transfer performance insidethe plug flow. The thermal resistance 1/σˆin,w and 1/σˆout,w are defined as below to evaluatethe overall heat transfer ability between the inner wall and the fluid, the outer wall and403.3. Heat transfer performance at inner iso-flux boundary (IFB)the fluid respectively. Since the outer wall is adiabatic under this boundary condition, theroute for 1/σˆout,w is cut, and all the heat flows into the capacitor (plug) to increase themean temperature ˆ¯T .σˆin,w−f = − ηβ∫ lˆ0(∂Tˆ∂r )rˆ=ηdzˆˆ¯Tin,w − ˆ¯T. (3.25)σˆout,w−f =1β∫ lˆ0(∂Tˆ∂r )rˆ=1dzˆˆ¯Tout,w − ˆ¯T. (3.26)Substitute eq. (3.12) into eq. (3.25), the thermal resistance under this boundary con-dition (IFB) is obtained below. The simple relationship in eq. (3.27) indicates that thehigher the conductance is, the relatively lower inner wall temperature (comparing to themean temperature) is achieved, and it is less possible to cause device failure due to localhigh temperature.σˆplug = σˆin,w−f =η1− η1ˆ¯Tin,w − ˆ¯T. (3.27)Since both the heat input and the total plug capacity are proportional to the plug lengthlˆ. It is better to define Cˆplug as the thermal capacity per length for comparison betweencases with different lengths. Then Cˆplug can be calculated using control volume method(with the results in eq. (3.16)),Cˆplug =Pe(1− η2)2=PeVˆpluglˆ. (3.28)Since the properties has been assumed to be constant and uniform, the thermal capac-ity for the volume taken up by the inner and the outer vortex can also be calculated bycomparing eqs. (2.38) to (2.40) and (3.28),Cˆin,v =PeVˆinlˆ. (3.29)Cˆout,v =PeVˆoutlˆ. (3.30)413.3. Heat transfer performance at inner iso-flux boundary (IFB)It is easily validated that Cˆin,v/Cˆout,v = Vˆin,v/Vˆout,v. This actually applies to anyextensive physical properties as long as they are distributed uniformly. Thus in the laterdiscussions, no specific difference needs to be addressed between one kind of capacity ratioand another such as the thermal capacity ratio and the mass capacity (volume) ratio.In the following subsections, the influences of the plug aspect ratio β, the inner-outerradius ratio η, and the Peclect number Pe are taken into account. Considering that in theheat exchanger design work, the thermal field is fully developed for the most part of thedevice, only influences at this stage are investigated. Thus, the mixing index γ is recordedto assure that the thermal fully developed field is reached before the investigation.3.3.2 Influences of the inner-outer radius ratio ηIn fig. 3.3 (a). there presents the developing of the thermal conductance σˆplug undervarying inner-outer radius ratios η, while the other two parameters are fixed at (β,Pe) =(2, 200) during the simulations. σˆplug experiences a decreasing process after the starting ofthe heat transfer and reaches a constant value after a period of time. The dimensionless timefor sampling here is tˆ = 58, and the conductance σˆplug for all thermal fields are independentof tˆ afterwards, thus the fully developed stage is reached.The influence of the inner-outer radius ratio η is plotted in fig. 3.3 (b). The thermalconductance σˆplug increases with increasing η, because (1). the increasing averaged recircu-lation flux (rˆ ˆ¯ur) at the inner vortex enhances the advection near the inner wall (fig. 2.9),(2). the capacity ratio Cˆin,v/Cˆout,v of the inner vortex increases (fig. 2.6), and the plug meantemperature ˆ¯T is therefore closer to the average inner wall temperature ˆ¯Tin,w,ˆ¯Tin,w − ˆ¯T isrelatively lower, σˆplug is therefore higher (eq. (3.27)).The conductance for single-phase flow heat transfer σˆsp and the enhancement denotedby σˆplug/σˆsp are plotted in fig. 3.3 (b). It is found that the enhancement to single-phaseflow σˆplug/σˆsp first increases when η increases, it reaches a peak value of 2.78 at η = 0.55,and then decreases when η increases furthermore. All of the enhancements at different ηare larger than 1.46 (obtained at η=0.90), which means there exist enhancement to thesingle-phase flow for all the inner-outer radius ratio η.423.3. Heat transfer performance at inner iso-flux boundary (IFB)Figure 3.3: (a). Development of the thermal conductance σˆplug, (b). influence of theinner-outer radius ratio η upon σˆplug and the enhancement by plug flow comparing to thesingle-phase flow heat transfer σˆplug/σˆsp. (β,Pe) = (2, 100)433.3. Heat transfer performance at inner iso-flux boundary (IFB)3.3.3 Influences of the aspect ratio βIn fig. 3.4 it plots the influence of plug dimensionless length β upon the thermal con-ductance σˆplug, while the Peclet number Pe is fixed at Pe=200 during the simulations. Infig. 3.3 (a). σˆplug decreases when β increases, while its increasing trend with increasing ηremains still.This decreasing effect with growing β is already predicted earlier in both section 2.2.4 andsection 3.3.1. Since the recirculation velocity, which determines the strength of advection,decreases with growing β. Moreover, the increasing length also increases the distance foreach recirculation trajectory, which highly increases the recirculation period tˆrec (fig. 2.10),reduces the recirculation frequency, and reduces the heat transfer enhancement.In fig. 3.4 (b). the enhancement comparing to the single-phase flow heat transferσˆplug/σˆsp decreases when β increases. When η = 0.9, σˆplug/σˆsp drops slightly and re-mains about 1.42 when β increases from 1 to 10. While at η=0.1, σˆplug/σˆsp drops from 2.02to 1.10 when β experiences the same increase. This indicates that for the lower inner-outerradius ratio η, the enhancement tends to be more affected by the plug length β. One of thereason can be that, at lower η the capacity ratio σˆplug/σˆsp decreases more obviously withthe growing β (fig. 2.6), and the plug mean temperature ˆ¯T is thus away from inner walltemperature ˆ¯Tin,w, the differenceˆ¯Tin,w− ˆ¯T increases more, and σˆplug drops more (eq. (3.27)).Under each β, the enhancement still increases first, reaches the peak and then decreasesagain with the growing η. The peak value of each σˆplug/σˆsp curve drops due to the dimin-ishing of the internal recirculation, while the η to reach the peak value also increases whenβ increases. This interaction of β and η makes the enhancement to single-phase flow morecomplex, and thus the findings here can be used to search for the optimum design.3.3.4 Influences of the Peclet number PeThe influence of Pe is plotted in fig. 3.5 while the dimensionless plug length is fixed atβ=2 during the simulations. In fig. 3.5 (a). the thermal conductance σˆplug increases whenPe grows. Meanwhile, it is observed in fig. 3.5 (b). that the enhancement to the single-phaseflow σˆplug/σˆsp also increases greatly when Pe grows. The mathematic explanation behindthis is obvious that the total capacity of the plug Cˆplug = Pe(1 − η2)/2 is proportional to443.3. Heat transfer performance at inner iso-flux boundary (IFB)Figure 3.4: Influence of plug dimensionless length β upon (a). the thermal conductanceσˆplug, (b). the enhancement by plug flow comparing to the single-phase flow heat transferσˆplug/σˆsp. The arrows in figures mark the direction of increasing β. Pe = 100.453.4. Heat transfer performance at the inner iso-thermal boundaries (ITB), comparison to IFBthe increasing Pe (eq. (3.27)), the amplitude of temperature difference ˆ¯Tin,w − ˆ¯T then shalldecreases with comparable scale, and the conductance therefore becomes higher (eq. (3.27)).The physics behind this explanation is that higher Pe means relatively higher heat capacityor convection ability, more portion of heat is carried by the hot fluid directly into thecolder area away from the wall. Instead of diffusion in a single direction (inner wall-outerwall), this portion of heat diffuses into the vortex center/the outer vortex when advectedby the fluid on the closed stream line, the diffusion area is larger, and thus the overall heattransfer performance is stronger. In fig. 3.5 (b). η to reach the peak value increases whenPe increases. The result is similar to that in fig. 3.3 (b). and is therefore useful for thedesign work.Despite the findings in this section at IFB. In the real heat exchangers, there existmore types of boundary conditions and the whole process can be more complex. Thus, thediscussions of the heat transfer at 2 other typical boundary conditions are carried out inthe following 2 sections. Lots of the conclusion at this 3 boundary conditions are similar,and more focus will be made upon the differences between them.3.4 Heat transfer performance at the inner iso-thermalboundaries (ITB), comparison to IFBUnlike the inner iso-flux boundary condition, the outer wall here is not adiabatic andthere exists the heat output through the outer wall. The developments of the heat inputQˆin, the heat through the vortex boundary Qˆb and the heat output Qˆout for a typical plugare plotted in fig. 3.6 (a)., while the other working conditions are the same as in fig. 3.2The heat input Qˆin reaches a high value as soon as the heat transfer starts and graduallydecreases. Heat through the vortex boundary Qˆb initiates at about tˆ = 0.5 and graduallyincreases. At last, heat output through the outer wall Qˆout initiates at tˆ = 3.0 and graduallyincreases. Three heat transfer rates obey the following relationship Qˆout ≤ Qˆb ≤ Qˆin, whichyield positive heat accumulation and temperature growing in both vortexes (Qˆin − Qˆb ≥ 0and Qˆb − Qˆout ≥ 0). The equitation is obtained when the thermal field develops longenough (tˆ > 26.0), and there exists zero heat accumulation in the plug (Qˆout = Qˆb = Qˆin),463.4. Heat transfer performance at the inner iso-thermal boundaries (ITB), comparison to IFBFigure 3.5: Influence of Peclet number Pe upon (a). the thermal conductance σˆplug, (b).the enhancement by plug flow comparing to the single-phase flow heat transfer σˆplug/σˆsp.The arrows in figures mark the direction of increasing Pe. β = 2.473.4. Heat transfer performance at the inner iso-thermal boundaries (ITB), comparison to IFBthe mean temperature ˆ¯T will no longer increase as well. Since the mixing index can besimplified into the plug mean temperature at ITB γ =ˆ¯T−Tˆout,wTˆin,w−Tˆout,w =ˆ¯T−01−0 =ˆ¯T , the meantemperature ˆ¯T can be used as the index to assure the thermal fields are fully developedduring the simulations.Figure 3.6: Heat transfer process at ITB, (β, η,Pe) = (2, 0.5, 100). (a). Variation of themean temperature ˆ¯T , and the dimensionless heat flow at the inner wall Qˆin, at the outerwall Qˆout and at the boundary between two vortexes Qˆb. The black solid points correspondto the examples shown in (b)., (b). (I) ∼ (IV) Sequence of dimensionless temperature Tˆdistribution at tˆI∼IV = (0.5, 3, 9.5, 25).The thermal network for the iso-thermal boundary condition is then extracted andplotted in fig. 3.7 (b). Comparing to that under the inner iso-flux boundary conditionin fig. 3.7 (a). the power source is changed from constant flux into constant temperature483.4. Heat transfer performance at the inner iso-thermal boundaries (ITB), comparison to IFBˆ¯Tout,w = 1. There exists both the resistance between the inner wall and the fluid, andthat between the fluid and the outer wall since the outer wall is not adiabatic. The totalresistance is thus 1/σˆin,w−f + 1/σˆout,w−f = (σˆin,w−f + σˆout,w−f )/(σˆin,w−f σˆout,w−f ) , whichyields the overall plug conductance σˆplug = (σˆin,w−f σˆout,w−f )/(σˆin,w−f + σˆout,w−f ) underthe iso-thermal boundary condition. At the stage of the fully thermal developed field, thereis zero flow into the capacitor Cˆplug, which exactly corresponds to zero heat accumulationin the whole plug.Figure 3.7: Comparison between IFB and ITB (a) (b). the thermal networks, (c). thethermal conductance σˆplug, (d). the enhancement by plug flow comparing to the single-phaseflow heat transfer σˆplug/σˆsp. The working condition is the same as in fig. 3.3.The plug conductance, as well as the enhancement comparing to the single-phase flow,are then plotted in fig. 3.7 (c). and (d). respectively. It can be concluded that changingthe boundary condition to the iso-thermal type only lowers the amplitude of the thermalconductance σˆplug and the enhancement σˆplug/σˆsp, while their trends versus varying workingconditions like η are not much affected. The reason can be the existence of an additional493.5. Heat transfer performance at the outer iso-flux boundaries (OFB), comparison to IFBresistor between the outer wall and the fluid 1/σˆout,w−f , which has been discussed above.3.5 Heat transfer performance at the outer iso-fluxboundaries (OFB), comparison to IFBIf the heat is loaded at the outer wall, the overall heat transfer performance shouldvary as well, because the heat path has reversed, and because the two vortexes are notsymmetric to each other. In the thermal network (fig. 3.8 (a). and (b).), the overall plugconductance is not the same one at two boundary types. For the inner iso-flux boundarycondition σˆplug = σˆin,w−f is the conductance between the inner wall and the fluid, whileσˆplug = σˆout,w−f is that between the outer wall and the fluid for the outer iso-flux boundarycondition.The plug conductance for both boundary types are plotted in fig. 3.8 (c). that σˆplug isalways higher under the outer iso-flux boundary condition. This is due to the stronger outervortex with higher radial transport velocity (fig. 2.9) and higher capacity (fig. 2.6). Thedifference of σˆplug between two boundary types gradually decreases when η grows, becausethe difference between the two vortexes gradually diminishes.In fig. 3.8 (d). it is found that the enhancement comparing to the single-phase flowσˆplug/σˆsp is a decreasing function of the inner-outer radius ratio η under the outer iso-fluxboundary condition. This is a result of both decreasing (1). the recirculation flux (rˆ ˆ¯ur)(fig. 2.9), and (2). The capacity of the outer vortex Cˆout,v (fig. 2.6). When η → 1 the twovortexes become nearly identical to each other, thus σˆplug/σˆsp under two boundary typeshave a trend to converge with each other, both of which can be replaced using the resultsin a slit channel.The results in this section enhance the fact that different thermal conductance and en-hancement to the single-phase flow vary with the different boundary types [18], it also varieswhen the heat is loaded from the outside or the inside due to the asymmetric vortexes. Sincethe correlations for plug flow heat transfer is far away from being sufficient for designingthe plug flow-based heat exchanger [59], if any correlations are to be proposed in the futurework, these effects should be taken into consideration.503.5. Heat transfer performance at the outer iso-flux boundaries (OFB), comparison to IFBFigure 3.8: Comparison between IFB and OFB (a) ∼ (b). the thermal networks, (c). thethermal conductance σˆplug, (d). the enhancement by plug flow comparing to the single-phaseflow heat transfer σˆplug/σˆsp. The working condition is the same as in fig. 3.3.513.6. Results for steady heat transfer under fully developed stage3.6 Results for steady heat transfer under fully developedstage3.6.1 Simplification based upon control volume methodAn interesting observation is in section 3.3.1, that the growth rate of ˆ¯T , ˆ¯Tin,w andˆ¯Tout,wis identical at fully developed stage. It indicates a possibility that the conclusion for singlephased flow mentioned in [68], that ∂Tˆ∂Zˆis uniform inside the channel, also suits for the2-D scenario in plug flow heat transfer. To examine that, the distribution of ∂Tˆ∂tˆat thefully developed stage as below under the same condition as in fig. 3.2 is plotted out. Thestandard deviation of ∂Tˆ∂tˆis also calculated.Figure 3.9: Temporal derivative of the thermal field ∂Tˆ∂tˆat tˆ = 50 ,at IFB, all other param-eters are same as in fig. 3.2As shown in fig. 3.9, the distribution of ∂Tˆ∂tˆis quite uniform, the value is about 6.7 ×10−3 = 2×0.5100×(1+0.5) =2ηPe(1+η) , and the standard deviation of∂Tˆ∂tˆis as low as 1.4431× 10−5,which is 0.2202% of the mean value of ∂Tˆ∂tˆ. This can be a validation to the hypothesisthat the uniform distribution of ∂Tˆ∂tˆis obtained in the plug flow heat transfer under thefully developed stage. With the help of above conclusion, a demonstration of deriving thesolution of fully developed thermal field at IFB will be presented below, while the other twoboundary conditions can be processed using a similar fashion.523.6. Results for steady heat transfer under fully developed stageAt IFB, the growth rate can be calculated by the control volume method, which is2ηPe(1+η) . The governing equation under fully developed stage is2ηPe(1− η2) + (1− η)(uˆrrˆ∂(Tˆ rˆ)∂rˆ+ uˆz∂Tˆ∂zˆ) =(1− η)2Pe[1rˆ∂∂rˆ(rˆ∂Tˆ∂rˆ) +∂2Tˆ∂zˆ2]. (3.31)Figure 3.10: Comparison between (a). thermal field under simplified model for fully devel-oped stage, and (b). that at tˆ = 60 using transient method, at OFB, all other parametersare same as in fig. 3.2Where the definition for the the dimensionless temperature Tˆ = T∆Tstd −2ηPe(1+η) tˆ−const,which is independent of the developing time tˆ. As mentioned earlier, when evaluatingthe heat transfer performance the relative shape, or say the spatial distribution instead533.6. Results for steady heat transfer under fully developed stageRange of Pe SOR0 ∼ 500 1− (3/2500)Pe500 ∼ 1000 200/Pe− Pe/4000 + 1/8> 1000 75/PeTable 3.2: Set-up for relaxation ratio (SOR) for steady state heat transferof the absolute value is important. Thus, taking away the term 2ηPe(1+η) tˆ will not effectthe calculation of, for example, σˆin,w−f in eq. (3.25). Other boundary conditions are thesame as in eq. (3.12), though it should be noted that there lack the boundary conditionfor temperature values. To bound the system, the temperature at a certain location isartificially set to zero (Tˆ (η, 0) = 0). Then eq. (3.31) is discretized using the finite volumemethod (details provided in table 3.1), though the time step is not involved in this problem.The equation is iterated using the Jacobi method with a relaxation ratio (SOR) based uponPe. Details about SOR are provided in table 3.2.A comparison between the result of this method at OFB and that from the transientmethod at tˆ = 60 is plotted as below. There exist nearly no difference in the relative shapebetween two contours, the results based on spatial distribution such as σˆ will not be affectedat all. Thus this can be treated as a validation.With the help of this simplified model, calculating the thermal field under 3× 19× 19×11 = 11, 913 kinds of commonly seen working conditions (coefficients are for numbers ofboundary conditions, β, η, Pe respectively) was conducted in an acceptable time, and theresults are stored in a database for later study.3.6.2 Some results and potential future design workA popular dimensionless parameter Nusselt number (Nu) is chosen to describe the resultsunder all 11, 913 working conditions. The definition of Nu isNu =hDhk. (3.32)Where h is the convective heat transfer coefficient, and Dh =4pi(r21−r22)2pi(r1+r2)= 2(r1−r2). ForIFB and ITB,543.6. Results for steady heat transfer under fully developed stageh =∫ l0(−k ∂T∂r )r=r1dz∫ l0(Tin,w − T¯ )r=r1dz. (3.33)For OFB,h =∫ l0(−k ∂T∂r )r=r2dz∫ l0(Tin,w − T¯ )r=r2dz. (3.34)Substitute eqs. (3.4) to (3.6) and (3.24) into eqs. (3.33) and (3.34)For IFB,Nu =2T¯in,w − ˆ¯T. (3.35)For ITB,Nu = − 2β∫ lˆ0(∂Tˆ∂rˆ )rˆ=ηdzˆ1− ˆ¯T. (3.36)For OFB,Nu =2T¯out,w − ˆ¯T. (3.37)It is easily validated that Nu/Nusp = σ/σsp, because one can derive from, for instance,eqs. (3.25) and (3.35) that Nu = f(η)σ and Nusp = f(η)σsp, where f(η) only depends on theboundary for sampling. Thus, mathematically speaking its equivalent to use either of themto describe the heat transfer enhancement. Here Nu is chosen only due to its popularityamong engineers and there is no need to analyze the thermal network in this subsection.The distribution of Nu versus β and η is plotted using the mesh plot function in MATLAB2016A, where each subplot is under a certain Pe.These plots along with the datasets can be quoted in the future design work, which ispreferred especially for optimization/inverse problem. For example, when Pe for a particularworking condition is near anyone in our database, one can efficiently use the interpolation toobtain the required outputs. Even if Pe exceeds any of those in the database, calculating thethermal field is quite fast with the help of analytical flow field solution and our simplified553.6. Results for steady heat transfer under fully developed stageFigure 3.11: I of II: The sequence of plots for Nu under Pe = 1, 2, 5, 10, 20, 50 at IFBcondition.563.6. Results for steady heat transfer under fully developed stageFigure 3.12: II of II: The sequence of plots for Nu under Pe = 100, 200, 500, 1000, 2000 atIFB condition.573.6. Results for steady heat transfer under fully developed stageFigure 3.13: I of II: The sequence of plots for Nu under Pe = 1, 2, 5, 10, 20, 50 at ITBcondition.583.6. Results for steady heat transfer under fully developed stageFigure 3.14: II of II: The sequence of plots for Nu under Pe = 100, 200, 500, 1000, 2000 atITB condition.593.6. Results for steady heat transfer under fully developed stageFigure 3.15: I of II: The sequence of plots for Nu under Pe = 1, 2, 5, 10, 20, 50 at OFBcondition.603.6. Results for steady heat transfer under fully developed stageFigure 3.16: II of II: The sequence of plots for Nu under Pe = 100, 200, 500, 1000, 2000 atOFB condition.613.7. Summary of this chaptermodel. The advantage would be more evident with the increasing amounts of workingconditions needed in order to search for an optimal design.3.7 Summary of this chapterThe numerical simulation for the heat transfer process at three kinds of boundary con-ditions in gas-liquid plug flow in tube-in-tube microchannels is carried out in MATLAB2016A. The heat transfer process is analyzed and extracted into simplified thermal net-works. The mixing index γ is used to assure the simulation is continuing until the thermalfield is fully developed. At the fully thermal developed field, the influences of the plug lengthβ, the inner-outer radius ratio η and the Peclet number Pe upon the plug thermal resis-tance σˆplug and upon the enhancement to the single-phase flow σˆplug/σˆsp are investigated.The difference between heat transfer performance between the inner iso-flux boundaries andother two boundaries are presented and explained using their thermal networks. Summaryof findings are listed as follows:• The process of heat transfer inside the plug flow at a certain boundary type can besimplified using a thermal network. Different boundary conditions determine whethersome routes containing resistors should be cut, they also determine the location andthe type of the heat source, as well as the orientation of the heat transfer path.• At the inner iso-flux condition (IFB), growing η leads to the higher thermal capacityratio of the inner vortex Cˆin,v/Cˆout,v as well as higher recirculation flux (rˆ ˆ¯ur) at theinner vortex, which lowers the temperature gap between the inner wall and the meanvalue of the plug ˆ¯Tin,w− ˆ¯T and increases the thermal conductance σˆplug. The variationof the enhancement to single-phase flow heat transfer with varying η is not singular,and the peak enhancement is reached at about η = 0.55 when (β,Pe) = (2, 100).• Longer plugs have lower σˆplug and σˆplug/σˆsp, which is caused by lower (rˆ ˆ¯ur) anddramatically longer recirculation period tˆrec. With growing β, the η corresponding tothe peak of σˆplug/σˆsp also increases.• Higher Pe increases the thermal capacity of the whole plug Cˆplug, which leads to more623.7. Summary of this chapterportion of heat being directly advected into the plug and enhances the combiningeffect of advection and diffusion, σˆplug and σˆplug/σˆsp are enhanced, therefore.• At the iso-thermal boundary condition (ITB), lower amplitudes of both σˆplug andσˆplug/σˆsp are reached comparing to those at the inner iso-flux boundary condition.However, their trends versus different working conditions like η are not affected. Thisvariance in the amplitudes is caused by an additional resistor between the outer walland the fluid 1/σˆout,w−f .• At the outer iso-flux boundaries (OFB), both σˆplug and σˆplug/σˆsp are always higherthan those under the inner iso-flux boundary condition, because both the capacity ofthe outer vortex Cˆout,v and the recirculation flux (rˆ ˆ¯ur) of the outer vortex are largerthan those of the inner vortex. When η → 1, all these parameters tend to convergeat both IFB and OFB because now the plug is like the plug inside of a slit channel.The difference between these three boundary types enhances the fact that differentcorrelations should be developed when the boundary condition type changes, or theheat transfer path reverses. It can be a focus of the future work.• A simplified model for heat transfer at the fully thermal developed stage was intro-duced to save simulation time. About 12, 000 cases were simulated in an acceptabletime, Nusselt number of which were extracted and mesh plotted for design work needin the future.63Chapter 4Analytical study of liquid-liquidplug flow in circular microchannelThe liquid-liquid plug flow in microchannel has different boundary conditions comparingto the gas-liquid plug flow. In gas-liquid flow, the viscous force at the two ends of the liquidplug is neglected since the low viscosity of the gas plug. In liquid-liquid plug flow, thisinterfacial interaction between the two plugs can no longer be ignored. Also, it makes theboundary condition, along with the analytical solutions, more complicated.In this chapter, firstly the focus is on the mathematical modification towards the bound-ary conditions, as well as towards the process of obtaining the analytical solutions in theliquid-liquid plug flow in the microchannel with a round cross-section. Then, some discus-sions are made upon phenomena caused by the interface interaction. Finally, the influenceof the viscosity ratio of the two plugs upon the skin friction coefficient will be presented.4.1 Mathematic modeling4.1.1 Governing equations and boundary conditionsApply the same basic assumptions in chapter 2, and use two moving coordinators stickingto each plug respectively, the liquid-liquid plug flow can be modeled using two PDEs with theplug index i = 1, 2; j = 2, 1 in the substrate respectively. Here the nondimensionalizationprocess has been skipped since it is basically the same as in chapter 2.{Lˆ4−1ψˆ}i = 0, (4.1)Two of the assumptions about boundary conditions can be still adapted from chapter 2:644.1. Mathematic modelingFigure 4.1: Schematic show for liquid-liquid plug flow in microchannel with round cross-section.1. The stream function is zero at all boundaries, 2. there is no slip between the walls andliquid plug.{ψˆ(0, rˆ)}i = 0, {ψˆ(β, rˆ)}i = 0, {ψˆ(zˆ, 1)}i = 0 (4.2){ 1rˆ∂ψˆ∂rˆ(zˆ, 1)}i = −1, (4.3)However, at the two ends of each plug, there exists viscous force due to the velocitygradients. They can be described using: 1. continuity for velocities at two sides of theinterface, 2. continuity for the viscous force at two sides of the interface.{ 1rˆ∂ψˆ∂zˆ(0, rˆ)}i = { 1rˆ∂ψˆ∂zˆ(β, rˆ)}j , (4.4)µiµj{ 1rˆ∂2ψˆ∂zˆ2(0, rˆ)}i = { 1rˆ∂2ψˆ∂zˆ2(β, rˆ)}j . (4.5)From eqs. (4.4) and (4.5) it can be seen that PDEs for plug 1 and 2 are coupled together,thus two systems should be solved simultaneously.654.1. Mathematic modeling4.1.2 Analytical solutionTo obtain the series solution of eqs. (4.1) to (4.5), similar fashion can be found from amore general problem from [4]: the Stokes flow in a 2-D cavity given boundary conditionson both top/bottom walls and front/rear walls. In this work, the authors used 2 infiniteseries under 2 orthogonal basis functions of x, y respectively to describe the flow field.Any distribution in the x direction, such as boundary condition on the top or on thebottom, has an unique position in the functional space containing the orthogonal basisM(x) = 〈M1(x),M2(x) · · · 〉 = 〈Mm(x)〉. Where M(x) is found through the eigenvalueproblem. Then the index function for each basis can be found by substituting the basisinto original PDEs, denoted as g(y) = 〈g1(y), g2(y) · · · 〉 = 〈gm(y)〉. Coefficients for g(y)can be bounded by getting the inner product of the boundary condition and each basisMm(x) respectively, because of the orthogonality, only the terms of series m will be kept.Similarly, any distribution in the y direction has an unique position in the functional spacecontaining the orthogonal basis L(y) = 〈Ll(y)〉, which comes with an index function in thisspace f(x) = 〈fl(x)〉. The sum of 2 infinite series M(x) • g(y) + L(y) • f(x) thus can fitany distributions on both direction x, y, thus can fit any given boundary conditions on bothtop/bottom and front/rear walls.Figure 4.2: Schematic show for Stokes flow in a cavity with moving walls, a demonstrationof the problem in [4]. F ,G are generic boundary conditions.664.1. Mathematic modelingThus here the sum of 2 infinite series is used to describe the flow field,ψˆ = M(rˆ) • g(zˆ) +L(zˆ) • f(rˆ).The orthogonal basis L(zˆ) must obey∂2∂zˆ2L = −Ω2L, (4.6)where Ω is the eigenvalue matrix,Ω =ω1ω2. . . .The boundary condition for this set of basis isL(0) = L(β) = 0. (4.7)The basis from eqs. (4.6) and (4.7) is L(zˆ) = 〈sin(ω1zˆ), sin(ω2zˆ) · · · 〉 = 〈sin(ωlzˆ)〉, whereωl = lpi/β. This has been obtained already in chapter 2. Then the index vector should obey(∂2∂rˆ2− 1rˆ∂rˆ−Ω2)( ∂2∂rˆ2− 1rˆ∂rˆ−Ω2)f = 0. (4.8)From eq. (4.8), the index function can be obtainedf(rˆ) = 〈fl(rˆ)〉 = 〈 rˆI1(ωlrˆ)rˆ2I2(ωlrˆ) •ElFl〉,which is part of that in chapter 2. The terms of rˆK1(ωlrˆ), rˆ2K2(ωlrˆ) are dropped becausehere the flow is in the circular tube in stead of concentric tube, and these terms approachthe infinity when rˆ = 0.The orthogonal basis in rˆ direction can also be found in a similar fashion, which is(∂2∂rˆ2− 1rˆ∂rˆ)M = −X2M , (4.9)674.1. Mathematic modelingwhere the eigenvalue matrix isX =χ1χ2. . . .The boundary condition isM(0) = M(1) = 0. (4.10)The solution of eqs. (4.9) and (4.10) is M(zˆ) = 〈rˆJ1(χ1rˆ), rˆJ1(χ2rˆ) · · · 〉 = 〈rˆJ1(χmrˆ)〉,where χm is the mth zero root of J1(χ) and Jν is the νth order of Bessel function of thefirst kind. The index vector here obeys(∂2∂zˆ2−X2)( ∂2∂zˆ2−X2)g = 0. (4.11)From eq. (4.11), the index function can be obtainedg(zˆ) = 〈gm(zˆ)〉 = 〈sinh(χmzˆ)cosh(χmzˆ)zˆ sinh(χmzˆ)zˆ cosh(χmzˆ) •AmBmCmDm〉.Thus, the universal solution is the sum of 2 series,ψˆ = M(rˆ) • g(zˆ) +L(zˆ) • f(rˆ) =∞∑l=1sin(ωlzˆ)[ElrˆI1(ωlrˆ) + Flrˆ2I2(ωlrˆ)]+∞∑m=1rˆJ1(χmrˆ)[Am sinh(χmzˆ) +Bm cosh(χmzˆ) + Cmzˆ sinh(χmzˆ) +Dmzˆ cosh(χmzˆ)],(4.12)where A ∼ F are constant coefficients. Substitute eq. (4.12) into boundary conditionseqs. (4.2) to (4.5), then apply the finite Fourier transformation eq. (2.16) when dealing684.1. Mathematic modelingwith distribution in zˆ direction, and apply the Hankel transformation (defined below) whendealing with distribution in rˆ direction.Hm[ψˆ] =∫ 10ψˆJ1(χmrˆ)drˆ. (4.13)The zero value for stream function on the wall,{ElI1(ωl) + FlI2(ωl)}i = 0. (4.14)The zero value for stream function on the front/rear ends of the plug,{Bm}i = 0, {Am sinh(χmβ) + Cmβ sinh(χmβ) +Dmβ cosh(χmβ)}i = 0. (4.15)Continuity for velocity at two sides of the interface,{ J2(χm)22χm cosh(χmβm)sinh(χmβm) + χmβm cosh(χmβm)cosh(χmβm) + χmβm sinh(χmβm) •AmCmDm+∞∑l=1ωl(−1)l ∫ 10 rˆJ1(χmrˆ)I1(ωlrˆ)drˆ∫ 10 rˆ2J1(χmrˆ)I2(ωlrˆ)drˆ •ElFl}i= { J2(χm)22χm01 •AmCmDm+∞∑l=1ωl ∫ 10 rˆJ1(χmrˆ)I1(ωlrˆ)drˆ∫ 10 rˆ2J1(χmrˆ)I2(ωlrˆ)drˆ •ElFl}j . (4.16)The non-slip boundary on the wall{βωl2I0(ωl)I1(ωl) •ElFl+ ∞∑m=1χmJ0(χm)∫ β0 sinh(χmzˆ) sin(ωlzˆ)dzˆ∫ β0 zˆ sinh(χmzˆ) sin(ωlzˆ)dzˆ∫ β0 zˆ cosh(χmzˆ) sin(ωlzˆ)dzˆ •AmCmDm}i694.1. Mathematic modeling= {−1− (−1)lωl}i (4.17)Continuity for viscous force at two sides of the interface,µiµj{χ2m sinh(χmβ)2χm cosh(χmβ) + βχ2m sinh(χmβ)2χm sinh(χmβ) + βχ2m cosh(χmβ) •AmCmDm}i = {2χmCm}j . (4.18)By truncating both series, the infinity large system formed by eqs. (4.14) to (4.18) canbe truncated into a system of (N × 5× 2)× (N × 5× 2), where N is the maximum order ofthe truncated series, and all (N × 5× 2) coefficients can be obtained by solving this linearsystem. To reduce the rank of the matrix, it is recommended to reorganize the equationsstarting from eqs. (4.14), (4.15) and (4.18), where the coefficients for one series are nottangled with these from the other one. Choose an independent variable from each series atan order, thus all the other 3 coefficients at this order can be expressed by them. Then,substitute into eqs. (4.16) and (4.18), and the system is reduced into (N×2×2)×(N×2×2).In the practice of this thesis, Cm, Fl are chosen as the two sets of independent variables, andthen calculations for other three sets of coefficients are conducted after solving the system.This can save a considerate amount of calculation time.After bounding the coefficients in eq. (4.12), the flow field as well as the skin frictioncoefficients for each plug can be calculated. Refer to chapter 2 for the procedures since theyare basically identical.uˆz =∞∑l=1ωlsin(ωlzˆ)[ElI0(ωlrˆ) + FlrˆI1(ωlrˆ)]+∞∑m=1χmJ0(χmrˆ)[Am sinh(χmzˆ) + Cmzˆ sinh(χmzˆ) +Dmzˆ cosh(χmzˆ)], (4.19)uˆr = −∞∑l=1ωlcos(ωlzˆ)[ElI1(ωlrˆ) + FlrˆI2(ωlrˆ)]704.2. Results and discussions−∞∑m=1J1(χmrˆ){[Amχm + Cmχmzˆ +Dm] cosh(χmzˆ) + [Cm +Dmχmzˆ] sinh(χmzˆ)}, (4.20)Cf =4β∞∑l=1ωl[1− (−1)l][ElI1(ωl) + FlI0(ωl)]. (4.21)The sensitivity study is also carried out here. The objective is the averaged skin frictioncoefficient (eq. (4.23)). It is shown in fig. 4.3 that the averaged skin friction coefficient isstable after series number is larger than 61. Thus, 61 is chosen as the maximum seriesnumber.Figure 4.3: The sensitivity study for the analytical solution. Variation of averaged Cf withthe growing series number. β1 = 1, β2 = 1, µ2/µ1 = 16.4.2 Results and discussions4.2.1 Validate using previous numerical study [1]To validate the analytical solution, a case from [1] is chosen for comparison. The workingliquid for plug 1 is n-butyl acetate (l1 = 3.0 mm, µ1 = 7.370 × 10−4 Pa · s) while that for714.2. Results and discussionsplug 2 is water (l2 = 2.4 mm, µ2 = 8.090× 10−4 Pa · s). The inner diameter for the channelis 1.0 mm (r = 0.5 mm). Thus, the dimensionless parameters here is β1 = 6.0, β2 =4.8, µ2/µ1 = 1.2076. The comparison between the results from [1] and from the analyticalsolution is plotted in the fig. 4.4. As it shows, no obvious difference can be seen in twosubplots, most of the main characters such as well defined internal circulation and the turnaround of streamlines near the ends of plugs are well captured by our analytical solution.Thus, it can be treated as an validation.Figure 4.4: Comparison between plots of flow field (a) from simulation in [1] (re-printed withpermission) and (b) from the analytical solution in this work. β1 = 6.0, β2 = 4.8, µ2/µ1 =1.2076.4.2.2 Validate using the existing gas-liquid model [2]The gas-liquid plug flow model is a simplification of the liquid-liquid model, where theviscosity is neglected in the gas phase. Hence, two models should converge when µ2/µ1 islarge enough and phase 1 becomes the gas phase in the liquid-liquid model.The gas-liquid plug flow can either be obtained by setting the inner-outer radius ratioinfinity close to 1 (η → 1) in chapter 2 (though some numerical error will occur owing tothe large value of Bessel function K near the inner wall). Thus, it is preferred to obtainthe results from the previous work of Che, Wong, Neng and Nguyen [2] directly. The case724.2. Results and discussionsfor validate is β = 2 for the gas-liquid model and β1 = 1, β2 = 2, µ2/µ1 = 64 for theliquid-liquid model. The viscosity ratio µ2/µ1 = 64 is chosen because throughout my work,I find lots of variations with increasing viscosity ratio becomes stable when µ2/µ1 > 64(These variations will be discussed in the next a few subsections).Figure 4.5: Comparison between contours of stream function (a) from the analytical solutionof the gas-liquid model in [2] (here only the equations in this citation are used to plot thecontour) and (b) from the analytical solution of liquid-liquid model in this work. (c) Thedistribution of axial velocity uˆz on 3 sample lines zˆ = 1.00, 0.40, 0.20. β1 = 1, β2 =2, µ2/µ1 = 64 (β1, µ2/µ1 are not needed for (a)).The comparison between two results is plotted in fig. 4.5. No obvious difference can beseen from two contour plots. The axial velocity distribution is also nearly identical. Thus,this can be treated as a validation under the condition when the viscosity ratio µ2/µ1 is farfrom 1.734.2. Results and discussions4.2.3 Influence of viscosity ratio µ2/µ1As mentioned in [57, 59], the gradient of uˆr near the interface is the key for momentumtransfer as well as for generating cap vortexes. When µ2/µ1 is larger than 1, the gradient ofuˆr is different at two sides of the interface owing to the continuity of shear force. The plug1 with smaller viscosity will have larger ∂uˆr/∂zˆ, this might lead to the transverse velocitynear the interface because of the continuity of velocity, and in another word, might lead tosecondary vortex near this interface. A demonstration is shown in fig. 4.6.Figure 4.6: Schematic show of the reason for cap vortexes when µ2/µ1 is sufficiently large.744.2. Results and discussionsThe most natural way of strengthening or weakening the shear force is to vary theviscosity ratio of two plugs µ2/µ1. In fig. 4.7, the influence of viscosity ratio µ2/µ1 of twoplugs upon the flow pattern as well as secondary vortexes are presented, where the lengthsof two plugs are mounted at β1 = 1, β2 = 1. In fig. 4.7 (a) when µ2/µ1 is 1, an well definedvortex (main vortexes) can be seen in each plug with a shape of bullet head (head is towardsthe center of channel rˆ = 0). The ’cache’ zone can be seen near the interface of two fluids,where the velocity gradually becomes nearly zero in both plugs. In fig. 4.7 (b) when µ2/µ1is 2, the shape of the main vortex in the plug 1 becomes like triangle with curved head andthe ’cache’ zone takes more space. The main vortex of plug 2 however expands, the head(at rˆ = 0) of the curved triangle becomes wider. The cap vortexes, which develop insidethe ’cache’ zone and firstly appear when µ2/µ1 approaches 4 as in fig. 4.7 (c), graduallybecome stronger when µ2/µ1 increases to 16 and squeeze the space of original main vortexas in fig. 4.7 (d).Figure 4.7: Sequence of stream lines with increasing µ2/µ1 = (1, 2, 4, 16), where β1 = β2 = 1.In fig. 4.8 (b), the intensity of secondary vortex center ψˆmin is plotted against varyingµ2/µ1. Since the negative symbol here means the rotation of the vortex, thus the decreasingof ψˆmin in fig. 4.8 (b) yields the strengthening cap vortexes however in the inversed rotation754.2. Results and discussionsas that of the main vortex. The decreasing trend of ψˆmin is obvious when µ2/µ1 firstlyincreases because of stronger momentum transfer by shear force at the interface, and thenit becomes stable when ψˆmin reaches 64, where the plug 1 can be treated as a gas plug.Figure 4.8: The influence of increasing µ2/µ1 upon (a) ψˆmax the intensity of main vortexcenters , (b) ψˆmin the intensity of cap vortex centers in plug 1 and (c) Cˆcap the capacity(volume) of the cap vortex in plug 1.Similarly, in fig. 4.8 (c) an increasing trend of the volume taken up by the cap vortexesCˆcap (defined below) can be firstly observed, then it becomes stable when µ2/µ1 becomeslarger than 64. The variation of the intensity of main vortexes ψˆmax then can be easilyunderstood as in fig. 4.8 (a). When µ2/µ1 increases, the ’cache’ zone in plug 1 gradually de-velops and becomes cap vortexes which completes with its main vortex, thus ψˆmax decreasesuntil µ2/µ1 becomes 64. The inversed process and trending happen in plug 2 naturally.764.2. Results and discussionsCˆcap =∫ β0∫ 101− sgn(ψˆ)2rˆdrˆdzˆ. (4.22)4.2.4 Influence of plug lengths β1, β2As discussed in chapter 2 the circulation inside the gas-liquid plug flow can be highlyaffected by the lengths or say aspect ratios of plugs. Thus, the influences of the length ofplug 1 and plug 2 are investigated, respectively.Figure 4.9: Sequence of stream lines with increasing β1 = (1, 2, 4), where β2 = 1 andµ2/µ1 = 16.The sequence of streamlines of varying β1 is plotted inside fig. 4.9, where β2 = 1 andµ2/µ1 = 16. Two obvious findings are listed: 1. The cap vortexes shrink slightly in sizeand, thus 2. the shape of the main vortex in plug 1 varies from a triangle to a like trapezoid.However, no visible change can be observed in plug 2.The detailed influences of β1 are then plotted in fig. 4.10. From fig. 4.10 (a) it can beseen that the influence of β1 upon the main vortex in plug 1 is like that in gas-liquid plugflow. ψˆmax in plug 1 firstly increases, reaches the 1-d asymptotic limit and then becomesstable when β1 increases. The influence upon ψˆmax in plug 2 is not visible, which suitwell to the observation made earlier in fig. 4.9. In fig. 4.10 (b) and (c) it can be seen thatcap vortexes ψˆmin become slightly weaker while its capacity Cˆcap firstly decreases when β1774.2. Results and discussionsFigure 4.10: The influence of varying β1 upon (a) ψˆmax the intensity of main vortex centers, (b) ψˆmin the intensity of cap vortex centers in plug 1 and (c) Cˆcap the capacity (volume)of the cap vortex in plug 1.784.2. Results and discussionsincreases from 1 to 2, then increases again when β1 increases from 2 to 10.Figure 4.11: The influence of varying β2 upon (a) ψˆmax the intensity of main vortex centers, (b) ψˆmin the intensity of cap vortex centers in plug 1 and (c) Cˆcap the capacity (volume)of the cap vortex in plug 1.Similarly, the influences of varying β2 are plotted out in fig. 4.10. However, no obviousinfluence of any is observed besides that upon the main vortex in plug 2 itself.4.2.5 Skin friction coefficients CfThe skin friction coefficients for two plugs can be quoted from eq. (4.21), respectively.The total pressure loss can be calculated by repeating the calculation for that introducedby each plug using their skin friction coefficients. However, an overall coefficient to evaluate794.2. Results and discussionsthe pumping power without calculating for multiple times is preferred. It is easily foundout from the definition of skin friction coefficient that,∆P ∝ µβ.And since the total pressure loss is∆P = ∆P1 + ∆P2.Figure 4.12: The skin friction coefficients for plug 1, plug 2 and their mean value withvarying µ2/µ1, β1 = β2 = 1.Thus, the definition of the overall friction coefficient is as below. In fig. 4.12, the frictioncoefficients for plug 1, plug 2 and their mean value are presented, where the plug lengths aremounted at β1 = β2 = 1. Identical value is obtained when µ1 = µ2 naturally. When µ2/µ1804.3. Summary of this chapterincreases from 1 to 128, Cf,1 increases from 76.81 to 99.22, while Cf,2 decreases from 76.81to 54.41. The mean value for a two-plug period Cf is closer to the variation of Cf,2, whichdecreases from 76.81 to 54.76. The value of Cf at extremely high µ2/µ1 is very close tothat of gas-liquid flow in [2], since here the flow field for plug 2 is already close to gas-liquidflow field with the same length β2.Cf =µ1β1Cf,1 + µ2β2Cf,2µ1β1 + µ2β2. (4.23)The detailed influences of plug lengths β1 and β2 are presented in fig. 4.13. In fig. 4.13(a) It can be concluded that: 1. The mean skin friction coefficient Cf is a decreasingfunction of β1. 2. The influence of β1 is more obvious when the viscosity of two plugs arestill comparable (µ2/µ1 is small). When µ2/µ1 is large, as have seen in fig. 4.12 the meanCf can be close to Cf,2, and the influence caused by plug 1 is limited. 3. Cf is a decreasingfunction of µ2/µ1 when β1 is small, while it becomes an increasing function of µ2/µ1 whenβ1 > 2. The influence of β2 is simple as in fig. 4.13 (b) that Cf is a decreasing function ofboth β2 and µ2/µ1.4.3 Summary of this chapterThe interaction between the interface of liquid-liquid plug flow caused by shear force cancomplicate the flow field. Secondary vortexes near the caps of the fluid with significantlysmaller viscosity will be generated. In this chapter, 2 PDEs was setup simultaneously tomodel the liquid-liquid plug flow in microchannels with round cross-section. Detailed studyof the relationship between main and cap vortexes was carried out, the influence of pluglengths and viscosity ratio upon flow pattern and skin friction coefficient were discussed.Findings are summarized as below:• The cap vortexes would gradually appear when µ2/µ1 increases. It is formed owingto the velocity continuity on the interface, as well as the high gradient of velocityin the plug 1 with small viscosity caused by the continuity of shear force. The capvortexes become stronger as µ2/µ1 increases, the capacity (volume) of which alsoexpands simultaneously.814.3. Summary of this chapterFigure 4.13: Variations of mean skin friction Cf coefficients with increasing µ2/µ1, β1 =β2 = 1 under (a)β2 = 1, β1 = 1 ∼ 10, (b)β1 = 1, β2 = 1 ∼ 10824.3. Summary of this chapter• The intensity of main vortexes in two plugs show inversed trends to each other withincreasing µ2/µ1. ψˆmax in plug 1 drops since it competes with the cap vortexes inside,while that in plug 2 grows with increasing µ2/µ1 naturally.• The lengths of plugs will mainly influence their main vortex. The length of the plugwith small viscosity β1 has an impact upon the cap vortexes in itself, while the lengthof another plug β2 barely has an impact upon the cap vortexes.• The skin friction coefficient Cf is a decreasing function of all β1, β2 and µ2/µ1 but forone case, where β1 is large and where Cf becomes an increasing function of µ2/µ1.• All of the trends mentioned here become stable when µ2/µ1 becomes large sufficientlyand converges to that under the gas-liquid plug flow scenario with the same lengths,which is owing to infinitely small velocity gradient near the ends in plug 2, and henceits flow field is infinitely close to that in gas-liquid plug flow.83Chapter 5Conclusions, limitations andpotential future workPlug flow creates well defined and stable vortexes, which are owing to the circulatingfluid inside caused by interfaces between two immiscible fluids. It is very efficient in heattransfer enhancement comparing to single-phase flow. Meanwhile, it is also much easier tocontrol compared to flow with phase change. Inside microchannel, the dominance of viscousforce makes the momentum equation linear and solvable (the Stokes assumption), whichcan help save a considerable amount of time for the design work with a massive number ofcases to search within.In this thesis, I concentrated on finding analytical solutions to plug flow in microchannelswith different geometry and find out two unrevealed series solutions in the end: 1. gas-liquid plug flow in the concentric microchannel and 2. liquid-liquid plug flow in the circularmicrochannel. Then I systematically investigated the influences of multiple inputs uponboth flow pattern and skin friction coefficient in details.An application of the gas-liquid plug flow field in the concentric tube, the heat transferat three kinds of typical boundary conditions were simulated and investigated in detail.Enhancements by the plug flow heat transfer to single-phase flow heat transfer were revealedfor all types of boundary types. Three simplified thermal networks of heat transfer at eachboundary type were extracted and used to explain heat transfer performance under differentworking conditions. At the end of this numerical study, a simplified model used for the fullydeveloped thermal field was extracted with the help of the control volume method. Around12,000 cases were simulated in an acceptable amount of time, whose Nusselt numbers werestored for direct quoting to help the design work in the future.84Chapter 5. Conclusions, limitations and potential future workHowever, there still stand some limitations of this work due to the limited time andfacilities. I will talk about these limitations and give some suggestions to enhance theresults of this thesis and to create new findings in the future.First of all, though some validation and verification from either asymptotic cases orprevious researches have been carried out, no experiment in the author’s lab has been con-ducted relevant to this topic. As have illustrated in the introduction, there stand challengesto both advanced non-intrusive, micro-scale temperature measuring facilities and perfectlymatched materials to avoid failure caused by optical refraction. Unfortunately, there werelimited experienced hands relative to either at the beginning stage of this study, which ledto the hard decision on not doing experiments. However, if one can accurately and quanti-tatively capture the plug flow field in complicated, non-straight microchannels, merely theirmethodology will be a significant contribution.Secondly, by the very end of finishing this work, my supervisor and I found a correlationwhich can directly calculate the maximum/minimum stream functions. The correlation isnot redundant to the analytical solution since it saves the time of calculating and searchingthe whole field and only focuses on the key of enhancing mixing or heat transfer. Undernumerous working conditions, it was also found out the maximum/minimum stream func-tions (in fig. A.4) has a similar trend as the heat transfer enhancement at IFB and OFB (infig. 3.8). Though the very dense database and quick simulation based upon simplified modelis already great for engineering design need, this can still be a chance to reveal a quantita-tive connection between the plug flow field and thermal performance. Unfortunately owingto limited time (see details in appendix A), the correlation is only developed for the flowfield. To build the connection between circulation flux and the heat transfer enhancementcan be a focus in the future.Lastly, the limitation of the analytical model has been mentioned earlier in chapter 2.The influence of surface tension cannot be included in the analytical modeling. However, Ihave come up with an idea (thanks to Dr. Brinkerhoff’s hint) of treating the liquid film as aliquid-liquid interface, when the plug is not deformed too much. The boundary conditionshere may not be interpreted as the non-slip moving wall anymore, however can still bedescribed using a known distribution when the information about the liquid film is given85Chapter 5. Conclusions, limitations and potential future worksuch as the film thickness (which can set the calculation domain) and the slippery velocitycompared to the wall (which can help to determine the first order and second order boundaryconditions). Then by applying similar fashion in chapter 4 the analytical solution can stillbe obtained. The modification to heat transfer model can be even more straightforwardsince the flow is in laminar region, the heat transfer in the liquid film can be seen as aGreatz problem or even as pure conduction if the slippery is not strong. Overall, to includethe liquid film in the current analytical modeling can be a potential point for the futurework.86Bibliography[1] A. Ghaini, A. Mescher, and D. W. Agar, “Hydrodynamic studies of liquid–liquid slugflows in circular microchannels,” Chemical engineering science, vol. 66, no. 6, pp. 1168–1178, 2011.[2] Z. Che, T. N. Wong, and N.-T. Nguyen, “An analytical model for plug flow in micro-capillaries with circular cross section,” International Journal of Heat and Fluid Flow,vol. 32, no. 5, pp. 1005–1013, 2011.[3] ——, “Heat transfer enhancement by recirculating flow within liquid plugs in mi-crochannels,” International Journal of Heat and Mass Transfer, vol. 55, no. 7-8, pp.1947–1956, 2012.[4] V. Meleshko, “Steady stokes flow in a rectangular cavity,” Proc. R. Soc. Lond. A, vol.452, no. 1952, pp. 1999–2022, 1996.[5] H. A. Stone, A. D. Stroock, and A. Ajdari, “Engineering flows in small devices: mi-crofluidics toward a lab-on-a-chip,” Annu. Rev. Fluid Mech., vol. 36, pp. 381–411,2004.[6] S. Lansing, J. Vı´quez, H. Mart´ınez, R. Botero, and J. Martin, “Quantifying electric-ity generation and waste transformations in a low-cost, plug-flow anaerobic digestionsystem,” ecological engineering, vol. 34, no. 4, pp. 332–348, 2008.[7] M. T. Kreutzer, F. Kapteijn, J. A. Moulijn, and J. J. Heiszwolf, “Multiphase monolithreactors: chemical reaction engineering of segmented flow in microchannels,” ChemicalEngineering Science, vol. 60, no. 22, pp. 5895–5916, 2005.[8] S. Lansing, J. F. Martin, R. B. Botero, T. N. Da Silva, and E. D. Da Silva, “Methane87Bibliographyproduction in low-cost, unheated, plug-flow digesters treating swine manure and usedcooking grease,” Bioresource technology, vol. 101, no. 12, pp. 4362–4370, 2010.[9] B. K. Yen, A. Gu¨nther, M. A. Schmidt, K. F. Jensen, and M. G. Bawendi, “A micro-fabricated gas–liquid segmented flow reactor for high-temperature synthesis: the caseof cdse quantum dots,” Angewandte Chemie International Edition, vol. 44, no. 34, pp.5447–5451, 2005.[10] R. Sista, Z. Hua, P. Thwar, A. Sudarsan, V. Srinivasan, A. Eckhardt, M. Pollack, andV. Pamula, “Development of a digital microfluidic platform for point of care testing,”Lab on a Chip, vol. 8, no. 12, pp. 2091–2104, 2008.[11] T. Dixit and I. Ghosh, “Review of micro-and mini-channel heat sinks and heat ex-changers for single phase fluids,” Renewable and Sustainable Energy Reviews, vol. 41,pp. 1298–1311, 2015.[12] X. Gao and R. Li, “Effects of nozzle positioning on single-phase spray cooling,” Inter-national Journal of Heat and Mass Transfer, vol. 115, pp. 1247–1257, 2017.[13] J. Kim, “Spray cooling heat transfer: the state of the art,” International Journal ofHeat and Fluid Flow, vol. 28, no. 4, pp. 753–767, 2007.[14] L. Lin and R. Ponnappan, “Heat transfer characteristics of spray cooling in a closedloop,” International Journal of Heat and Mass Transfer, vol. 46, no. 20, pp. 3737–3746,2003.[15] K. Jambunathan, E. Lai, M. Moss, and B. Button, “A review of heat transfer data forsingle circular jet impingement,” International journal of heat and fluid flow, vol. 13,no. 2, pp. 106–115, 1992.[16] K. Daou, R. Wang, and Z. Xia, “Desiccant cooling air conditioning: a review,” Renew-able and Sustainable Energy Reviews, vol. 10, no. 2, pp. 55–77, 2006.[17] S. Launay, V. Sartre, and J. Bonjour, “Parametric analysis of loop heat pipe operation:a literature review,” International Journal of Thermal Sciences, vol. 46, no. 7, pp. 621–636, 2007.88Bibliography[18] T. Bandara, N.-T. Nguyen, and G. Rosengarten, “Slug flow heat transfer withoutphase change in microchannels: A review,” Chemical Engineering Science, vol. 126,pp. 283–295, 2015.[19] M. Asadi, G. Xie, and B. Sunden, “A review of heat transfer and pressure drop char-acteristics of single and two-phase microchannels,” International Journal of Heat andMass Transfer, vol. 79, pp. 34–53, 2014.[20] D. B. Tuckerman and R. F. W. Pease, “High-performance heat sinking for vlsi,” IEEEElectron device letters, vol. 2, no. 5, pp. 126–129, 1981.[21] A. Asthana, I. Zinovik, C. Weinmueller, and D. Poulikakos, “Significant heat transferenhancement in microchannels with a segmented flow of two immiscible liquids,” in2010 14th International Heat Transfer Conference. American Society of MechanicalEngineers, 2010, pp. 141–149.[22] ——, “Significant nusselt number increase in microchannels with a segmented flow oftwo immiscible liquids: An experimental study,” International Journal of Heat andMass Transfer, vol. 54, no. 7-8, pp. 1456–1464, 2011.[23] A. Bar-Cohen, M. Arik, and M. Ohadi, “Direct liquid cooling of high flux micro andnano electronic components,” Proceedings of the IEEE, vol. 94, no. 8, pp. 1549–1570,2006.[24] M. Meis, F. Varas, A. Vela´zquez, and J. Vega, “Heat transfer enhancement in micro-channels caused by vortex promoters,” International Journal of Heat and Mass Trans-fer, vol. 53, no. 1-3, pp. 29–40, 2010.[25] Y. Wang, F. Houshmand, D. Elcock, and Y. Peles, “Convective heat transfer andmixing enhancement in a microchannel with a pillar,” International Journal of Heatand Mass Transfer, vol. 62, pp. 553–561, 2013.[26] Y. Wang, A. Nayebzadeh, X. Yu, J.-H. Shin, and Y. Peles, “Local heat transfer in a mi-crochannel with a pin fin—experimental issues and methods to mitigate,” InternationalJournal of Heat and Mass Transfer, vol. 106, pp. 1191–1204, 2017.89Bibliography[27] W. Dean and J. Hurst, “Note on the motion of fluid in a curved pipe,” Mathematika,vol. 6, no. 1, pp. 77–85, 1959.[28] Y. Sui, C. Teo, P. S. Lee, Y. Chew, and C. Shu, “Fluid flow and heat transfer in wavymicrochannels,” International Journal of Heat and Mass Transfer, vol. 53, no. 13-14,pp. 2760–2772, 2010.[29] G. Xia, J. Jiang, J. Wang, Y. Zhai, and D. Ma, “Effects of different geometric structureson fluid flow and heat transfer performance in microchannel heat sinks,” InternationalJournal of Heat and Mass Transfer, vol. 80, pp. 439–447, 2015.[30] L. Chai, G. Xia, L. Wang, M. Zhou, and Z. Cui, “Heat transfer enhancement in mi-crochannel heat sinks with periodic expansion–constriction cross-sections,” Interna-tional Journal of Heat and Mass Transfer, vol. 62, pp. 741–751, 2013.[31] I. Mudawar and M. B. Bowers, “Ultra-high critical heat flux (chf) for subcooled waterflow boiling—i: Chf data and parametric effects for small diameter tubes,” Interna-tional Journal of Heat and Mass Transfer, vol. 42, no. 8, pp. 1405–1428, 1999.[32] D. Deng, W. Wan, Y. Qin, J. Zhang, and X. Chu, “Flow boiling enhancement ofstructured microchannels with micro pin fins,” International Journal of Heat and MassTransfer, vol. 105, pp. 338–349, 2017.[33] S. Krishnamurthy and Y. Peles, “Flow boiling heat transfer on micro pin fins entrenchedin a microchannel,” Journal of Heat Transfer, vol. 132, no. 4, p. 041007, 2010.[34] T. Y. Liu, P. Li, C. Liu, and C. Gau, “Boiling flow characteristics in microchannelswith very hydrophobic surface to super-hydrophilic surface,” International Journal ofHeat and Mass Transfer, vol. 54, no. 1-3, pp. 126–134, 2011.[35] C. Choi, J. S. Shin, D. I. Yu, and M. H. Kim, “Flow boiling behaviors in hydrophilicand hydrophobic microchannels,” Experimental Thermal and Fluid Science, vol. 35,no. 5, pp. 816–824, 2011.[36] A. Jaikumar and S. Kandlikar, “Pool boiling enhancement through bubble induced90Bibliographyconvective liquid flow in feeder microchannels,” Applied Physics Letters, vol. 108, no. 4,p. 041604, 2016.[37] A. R. Betz and D. Attinger, “Can segmented flow enhance heat transfer in microchannelheat sinks?” International Journal of Heat and Mass Transfer, vol. 53, no. 19-20, pp.3683–3691, 2010.[38] J. Sivasamy, Z. Che, T. N. Wong, N.-T. Nguyen, and L. Yobas, “A simple method forevaluating and predicting chaotic advection in microfluidic slugs,” Chemical engineer-ing science, vol. 65, no. 19, pp. 5382–5391, 2010.[39] P. Angeli and A. Gavriilidis, “Hydrodynamics of taylor flow in small channels: a re-view,” Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Me-chanical Engineering Science, vol. 222, no. 5, pp. 737–751, 2008.[40] Y. Lim, S. Yu, and N. Nguyen, “Flow visualization and heat transfer characteristics ofgas–liquid two-phase flow in microtube under constant heat flux at wall,” InternationalJournal of Heat and Mass Transfer, vol. 56, no. 1-2, pp. 350–359, 2013.[41] P. A. Walsh, E. J. Walsh, and Y. S. Muzychka, “Heat transfer model for gas–liquid slugflows under constant flux,” International Journal of Heat and Mass Transfer, vol. 53,no. 15-16, pp. 3193–3201, 2010.[42] Z. Che, T. N. Wong, and N.-T. Nguyen, “Heat transfer in plug flow in cylindricalmicrocapillaries with constant surface heat flux,” International Journal of ThermalSciences, vol. 64, pp. 204–212, 2013.[43] S. T. Wereley and C. D. Meinhart, “Recent advances in micro-particle image velocime-try,” Annual Review of Fluid Mechanics, vol. 42, pp. 557–576, 2010.[44] D. Sinton, “Microscale flow visualization,” Microfluidics and Nanofluidics, vol. 1, no. 1,pp. 2–21, 2004.[45] D. Ross, M. Gaitan, and L. E. Locascio, “Temperature measurement in microfluidicsystems using a temperature-dependent fluorescent dye,” Analytical chemistry, vol. 73,no. 17, pp. 4117–4123, 2001.91Bibliography[46] B.-J. Kim, Y. Z. Liu, and H.-J. Sung, “Micro piv measurement of two-fluid flow withdifferent refractive indices,” Measurement science and technology, vol. 15, no. 6, p.1097, 2004.[47] V. Talimi, Y. Muzychka, and S. Kocabiyik, “A review on numerical studies of slugflow hydrodynamics and heat transfer in microtubes and microchannels,” InternationalJournal of Multiphase Flow, vol. 39, pp. 88–104, 2012.[48] R. Gupta, D. F. Fletcher, and B. S. Haynes, “On the cfd modelling of taylor flow inmicrochannels,” Chemical Engineering Science, vol. 64, no. 12, pp. 2941–2950, 2009.[49] K. Fukagata, N. Kasagi, P. Ua-arayaporn, and T. Himeno, “Numerical simulation ofgas–liquid two-phase flow and convective heat transfer in a micro tube,” InternationalJournal of Heat and Fluid Flow, vol. 28, no. 1, pp. 72–82, 2007.[50] Z. Dai, Z. Guo, D. F. Fletcher, and B. S. Haynes, “Taylor flow heat transfer in mi-crochannels—unification of liquid–liquid and gas–liquid results,” Chemical EngineeringScience, vol. 138, pp. 140–152, 2015.[51] V. Talimi, Y. Muzychka, and S. Kocabiyik, “Slug flow heat transfer in square mi-crochannels,” International Journal of Heat and Mass Transfer, vol. 62, pp. 752–760,2013.[52] J. Zhang, D. F. Fletcher, and W. Li, “Heat transfer and pressure drop characteris-tics of gas–liquid taylor flow in mini ducts of square and rectangular cross-sections,”International Journal of Heat and Mass Transfer, vol. 103, pp. 45–56, 2016.[53] A. Gu¨nther, M. Jhunjhunwala, M. Thalmann, M. A. Schmidt, and K. F. Jensen, “Mi-cromixing of miscible liquids in segmented gas- liquid flow,” Langmuir, vol. 21, no. 4,pp. 1547–1555, 2005.[54] Z. Che, T. N. Wong, N.-T. Nguyen, and C. Yang, “Three dimensional features of con-vective heat transfer in droplet-based microchannel heat sinks,” International Journalof Heat and Mass Transfer, vol. 86, pp. 455–464, 2015.92Bibliography[55] S. K. R. Cherlo, S. Kariveti, and S. Pushpavanam, “Experimental and numerical in-vestigations of two-phase (liquid- liquid) flow behavior in rectangular microchannels,”Industrial & Engineering Chemistry Research, vol. 49, no. 2, pp. 893–899, 2009.[56] Q. Li, K. H. Luo, Q. Kang, Y. He, Q. Chen, and Q. Liu, “Lattice boltzmann methods formultiphase flow and phase-change heat transfer,” Progress in Energy and CombustionScience, vol. 52, pp. 62–105, 2016.[57] Z. Che, N.-T. Nguyen, T. N. Wong et al., “Analysis of chaotic mixing in plugs movingin meandering microchannels,” Physical Review E, vol. 84, no. 6, p. 066309, 2011.[58] M. Kashid, I. Gerlach, S. Goetz, J. Franzke, J. Acker, F. Platte, D. Agar, and S. Turek,“Internal circulation within the liquid slugs of a liquid- liquid slug-flow capillary mi-croreactor,” Industrial & engineering chemistry research, vol. 44, no. 14, pp. 5003–5010,2005.[59] A. Abdollahi, R. N. Sharma, and A. Vatani, “Fluid flow and heat transfer of liquid-liquid two phase flow in microchannels: A review,” International Communications inHeat and Mass Transfer, vol. 84, pp. 66–74, 2017.[60] Z. Che, T. N. Wong, N.-T. Nguyen, and C. Yang, “Asymmetric heat transfer in liquid–liquid segmented flow in microchannels,” International Journal of Heat and MassTransfer, vol. 77, pp. 385–394, 2014.[61] M. Shojaeian and A. Kos¸ar, “Convective heat transfer of non-newtonian power-lawslip flows and plug flows with variable thermophysical properties in parallel-plate andcircular microchannels,” International Journal of Thermal Sciences, vol. 100, pp. 155–168, 2016.[62] M. Shojaeian, M. Yildiz, and A. Kos¸ar, “Heat transfer characteristics of plug flowswith temperature-jump boundary conditions in parallel-plate channels and concentricannuli,” International Journal of Thermal Sciences, vol. 84, pp. 252–259, 2014.[63] Y. Muzychka, E. Walsh, and P. Walsh, “Simple models for laminar thermally develop-93Bibliographying slug flow in noncircular ducts and channels,” Journal of Heat Transfer, vol. 132,no. 11, p. 111702, 2010.[64] Z. Che, T. N. Wong, and N.-T. Nguyen, “An analytical model for a liquid plug movingin curved microchannels,” International Journal of Heat and Mass Transfer, vol. 53,no. 9-10, pp. 1977–1985, 2010.[65] C. Meyer, M. Hoffmann, and M. Schlu¨ter, “Micro-piv analysis of gas–liquid taylorflow in a vertical oriented square shaped fluidic channel,” International Journal ofMultiphase Flow, vol. 67, pp. 140–148, 2014.[66] Cramer and Gabriel, Introduction a l’analyse des lignes courbes algebriques par GabrielCramer... chez les freres Cramer & Cl. Philibert, 1750.[67] W. E. Langlois and M. O. Deville, Slow viscous flow. Springer, 1964.[68] T. L. Bergman, F. P. Incropera, D. P. DeWitt, and A. S. Lavine, Fundamentals of heatand mass transfer. John Wiley & Sons, 2011.94Appendix95Appendix ACorrelations for long plugs (β > 2)in chapter 2The analysis in this section will focus on long plugs, which have relatively large aspectratios, i.e. β > 2. Because as shown in chapter 2, the influence of β upon the volume ratioof the two vortexes Vˆin,v/Vˆout,v and the ratio of maximum and minimum stream functions|ψˆmin/ψˆmax| are negligible when β > 2. Thus these two parameters are dependent on ηonly. We will derive the correlation for them when β > 2 with the help of the followingfigure,Figure A.1: The volume ratio of two vortexes Vˆin,v/Vˆout,v, and the ratio of max/min streamfunctions |ψˆmin/ψˆmax| when β > 2. The solid lines are the correlations for the pointsobtained from analytical results in chapter 2.96Appendix A. Correlations for long plugs (β > 2) in chapter 2The solid lines in fig. A.1 are the correlation made. The volume ratio Vˆin,v/Vˆout,v showsa relationship with the radius ratio of the channel η, which isVˆin,v/Vˆout,v = η. (A.1)Similarly, fig. A.1 shows the ratio of max/min stream functions |ψˆmin/ψˆmax| varies withη as below,|ψˆmin/ψˆmax| = η1.16. (A.2)To further derive the correlations, define the total volumetric recirculation rate of theplug using the summation of those of two vortexes, which is |ψˆmax|+ |ψˆmin| = ψˆmax− ψˆmin.Figure A.2: The total of volumetric recirculation rate of two vortexes ψˆmax− ψˆmin becomesindependent of β when β > 2.In the small subplot in fig. A.2 the total recirculation rate grows and then becomesstable when plug length lˆ increases. This stable value is marked with a subscript max. The97Appendix A. Correlations for long plugs (β > 2) in chapter 2ratio (ψˆmax − ψˆmin)/(ψˆmax − ψˆmin)max becomes infinity close to 1 as well when β > 2.Hence, ψˆmax − ψˆmin can also be correlated using η only when β > 2.Figure A.3: The refresh frequency defined in eq. (A.3). The solid lines are the correlationsfor the points obtained from analytical results in chapter 2.Define the refresh period τˆ as below, which is the ratio between the volume of the wholeplug Vˆplug and the total recirculation rate ψˆmax− ψˆmin, here the result in eq. (2.40) is used,τˆ =Vˆplugψˆmax − ψˆmin=(1− η2)lˆ2(ψˆmax − ψˆmin). (A.3)If τˆ−1 is plotted against the plug length lˆ = β(1 − η) for multiple cases obtained fromchapter 2 with β > 2 as in fig. A.3, it can be clearly concluded that τˆ−1 is dependent on lˆonly. A simple inverse model is adopted here,τˆ−1 = 0.1973lˆ−1. (A.4)Then the total recirculation rate ψˆmax − ψˆmin is concluded from eqs. (A.3) and (A.4),98Appendix A. Correlations for long plugs (β > 2) in chapter 2ψˆmax − ψˆmin = 0.0986(1− η2). (A.5)The correlations for ψˆmax/min are concluded by combining eqs. (A.2) and (A.5),ψˆmax =0.0986(1− η2)1 + η1.16, (A.6)ψˆmin = −0.0986(1− η2)η1.161 + η1.16. (A.7)An important set of parameters for quantifying vortexes, the averaged recirculation fluxby the two vortexes (rˆ ˆ¯ur)out,v/in,v are defined in chapter 2. Thus they can be calculated bycombining eqs. (2.42) and (A.6) or eqs. (2.42) and (A.7), respectively.(rˆ ˆ¯ur)out,v =0.1973(1− η2)(1 + η1.16)lˆ, (A.8)(rˆ ˆ¯ur)in,v = −0.1973(1− η2)η1.16(1 + η1.16)lˆ. (A.9)A comparison between the analytical results in chapter 2 and these from correlationseqs. (A.5) to (A.7) for long plugs (β > 2) is plotted as below. For the inner vortex, thecorrelation in eq. (A.7) shows good agreement with all analytical results. For the outervortex, the correlation in eq. (A.7) agrees well with the data points except for plugs withsmall inner radius (η <∼ 0.2). As a result, eq. (A.5) predicts the total plug circulation ratewith good accuracy except for small inner radius (η <∼ 0.2). Overall the correlations canwell predict parameters such as maximum or minimum stream functions and the averagedrecirculation flux, the appropriate range for applying these correlations is β > 2 ∧ η >∼ 0.2.99Appendix A. Correlations for long plugs (β > 2) in chapter 2Figure A.4: Comparison between the correlation and the analytical results. The solid linesare the correlations for the points obtained from analytical results in chapter 2.100
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Analytical and numerical study of plug flow inside...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Analytical and numerical study of plug flow inside round/concentric microchannels Cao, Yadi 2018
pdf
Page Metadata
Item Metadata
Title | Analytical and numerical study of plug flow inside round/concentric microchannels |
Creator |
Cao, Yadi |
Publisher | University of British Columbia |
Date Issued | 2018 |
Description | Plug flow can generate recirculating flow between interfaces of immiscible fluids. Depending on the phase used to segment flow, it can be gas-liquid plug flow or liquid-liquid plug flow. The recirculating flow can enhance heat transfer as compared to the continuous pipe flow, especially in the laminar regime such as the microchannel flow. The present work focuses on the flow and heat transfer of liquid plugs with low Reynolds numbers. The flow is modeled by applying Stokes simplification, and the solution is obtained by solving fourth-order partial differential equation sets. Solutions of two types of plug flow are obtained: 1) the gas-liquid plug flow in the concentric microchannel; 2) the liquid-liquid plug flow in the circular microchannel. For the gas-liquid plug flow study, the flow patterns inside the liquid phase including the volume ratio of the inner and outer vortexes, the ratio of maximum-to-minimum stream functions, the averaged recirculation flux as well as the skin friction coefficient are investigated in details. Correlations for predicting the maximum and minimum of the stream function are developed. For the liquid-liquid plug flow study, the influences of plug lengths and the viscosity ratio upon the cap vortexes and the overall skin friction coefficient are studied in details. The heat transfer of the gas-liquid plug flow in the concentric microchannel is simulated numerically in MATLAB. Three types of thermal boundary conditions are investigated. The developing process of the thermal field can be explained using a simple thermal network for each boundary condition. The influences of parameters including the plug aspect ratio, the channel inner-outer radius ratio and the Peclect number upon the thermal conductance and heat transfer enhancement to the single-phase flow are investigated systematically. Then a simplified model for the fully developed thermal field is extracted for the quick calculation need in the design work. The results obtained from about 12,000 cases form a database that can be used in the future design work of heat exchanger based upon this kind of flow. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2018-09-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0372129 |
URI | http://hdl.handle.net/2429/67225 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Engineering, School of (Okanagan) |
Degree Grantor | University of British Columbia |
GraduationDate | 2018-11 |
Campus |
UBCO |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2018_november_cao_yadi.pdf [ 32MB ]
- Metadata
- JSON: 24-1.0372129.json
- JSON-LD: 24-1.0372129-ld.json
- RDF/XML (Pretty): 24-1.0372129-rdf.xml
- RDF/JSON: 24-1.0372129-rdf.json
- Turtle: 24-1.0372129-turtle.txt
- N-Triples: 24-1.0372129-rdf-ntriples.txt
- Original Record: 24-1.0372129-source.json
- Full Text
- 24-1.0372129-fulltext.txt
- Citation
- 24-1.0372129.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0372129/manifest