Applied Science, Faculty of
Engineering, School of (Okanagan)
DSpace
UBCO
Jalil Khabbazi, Arash
Master of Applied Science - MASc
University of British Columbia
Hydrogen, as an alternative source of energy, has the potential to minimize greenhouse gas (GHG) emissions. Through different sources and industrial processes, it is feasible to produce hydrogen. Therefore, blending certain percentages of hydrogen into the existing natural gas pipeline infrastructure can greatly decrease the carbon intensity of current gas grids while also promoting a more sustainable energy system. Nevertheless, the significant density contrast of almost eight to nine times between hydrogen and natural gas introduces a challenge as buoyant hydrogen tends to stratify and become trapped near the top wall once injected via a Tee junction. In this study, the hydrogen blending into the natural gas distribution pipelines using Tee junctions is numerically investigated, and the mixing homogeneity (uniformity) length is quantified through the coefficient of variation (CoV) of the hydrogen mole fraction. The effect of secondary (side) flow jet intensity on turbulent mixing in a distribution-pressure (DP) system is first looked into by altering the secondary pipe diameter while maintaining the same flow rate. Consequently, the differences between vertical top-side, horizontal, and vertical bottom-side injection configurations in an intermediate-pressure (IP) case are studied. The ideal and real gas equations of state (EoS) scenarios are next investigated in the IP case. Overall, a greater secondary flow jet intensity results in a shorter mixing homogeneity length, as investigated in the DP case. The mixing homogeneity length is shortened even further if the secondary flow penetrates deep into the primary (main) flow of natural gas without getting trapped near the top wall from the outset. Vertical bottom-side injection also results in a mixing homogeneity length that is roughly four times shorter than horizontal injection and five times shorter than vertical top-side injection in the IP case. When industry standards are followed, injecting from the bottom definitely outperforms other injection setups. Furthermore, using real gas EoS, such as Soave-Redlich-Kwong (SRK), becomes crucial when the operating pressure rises, resulting in a longer mixing homogeneity length.
Mixing Gaseous Hydrogen intoNatural Gas DistributionPipelinesbyArash Jalil KhabbaziA 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)November 2023© Arash Jalil Khabbazi, 2023The following individuals certify that they have read, and recommend tothe College of Graduate Studies for acceptance, a thesis entitled:Mixing Gaseous Hydrogen into Natural Gas DistributionPipelinessubmitted by Arash Jalil Khabbazi in partial fulfilment of the require-ments of the degree of Master of Applied ScienceRi Li, School of EngineeringSupervisorJoshua Brinkerhoff, School of EngineeringSupervisory Committee MemberSina Kheirkhah, School of EngineeringSupervisory Committee MemberKenneth Chau, School of EngineeringUniversity ExamineriiAbstractHydrogen, as an alternative source of energy, has the potential to min-imize greenhouse gas (GHG) emissions. Through different sources and in-dustrial processes, it is feasible to produce hydrogen. Therefore, blendingcertain percentages of hydrogen into the existing natural gas pipeline infras-tructure can greatly decrease the carbon intensity of current gas grids whilealso promoting a more sustainable energy system. Nevertheless, the signif-icant density contrast of almost eight to nine times between hydrogen andnatural gas introduces a challenge as buoyant hydrogen tends to stratify andbecome trapped near the top wall once injected via a Tee junction.In this study, the hydrogen blending into the natural gas distributionpipelines using Tee junctions is numerically investigated, and the mixinghomogeneity (uniformity) length is quantified through the coefficient of vari-ation (CoV) of the hydrogen mole fraction. The effect of secondary (side)flow jet intensity on turbulent mixing in a distribution-pressure (DP) sys-tem is first looked into by altering the secondary pipe diameter while main-taining the same flow rate. Consequently, the differences between verticaltop-side, horizontal, and vertical bottom-side injection configurations in anintermediate-pressure (IP) case are studied. The ideal and real gas equationsof state (EoS) scenarios are next investigated in the IP case.Overall, a greater secondary flow jet intensity results in a shorter mixinghomogeneity length, as investigated in the DP case. The mixing homogeneitylength is shortened even further if the secondary flow penetrates deep intothe primary (main) flow of natural gas without getting trapped near thetop wall from the outset. Vertical bottom-side injection also results in aiiimixing homogeneity length that is roughly four times shorter than horizontalinjection and five times shorter than vertical top-side injection in the IP case.When industry standards are followed, injecting from the bottom definitelyoutperforms other injection setups. Furthermore, using real gas EoS, such asSoave-Redlich-Kwong (SRK), becomes crucial when the operating pressurerises, resulting in a longer mixing homogeneity length.ivLay SummaryHydrogen, the most abundant chemical element in the universe, is analternative to other energy sources. Blending hydrogen into the existing nat-ural gas pipelines has the potential to reduce greenhouse gas (GHG) emis-sions while improving energy storage efficiency. Nonetheless, challenges suchas acceptable hydrogen concentrations in natural gas grid infrastructure andthe effects of hydrogen on pipeline safety, durability, and integrity should beconsidered. Due to the substantial density difference between the two gases,hydrogen tends to become trapped towards the top wall when injected via aTee junction. Thus, this thesis provides a numerical investigation of hydro-gen and natural gas mixing evolution in natural gas distribution pipelines.Given the significant thermophysical differences between the two gases, it iscritical to understand how hydrogen is injected into the natural gas pipeline.Furthermore, the operating pressure of the pipeline is significant in analyzingthis mixing behavior.vPrefaceThe numerical investigations provided in this thesis are the result of re-search undertaken by the author, Arash Jalil Khabbazi, as part of the Masterof Applied Science program in Mechanical Engineering at the University ofBritish Columbia (UBC) under the direct supervision of Dr. Sunny Ri Li.The project is supported by the Natural Sciences and Engineering ResearchCouncil of Canada (NSERC). Model design, code development, numericalsimulations, and manuscript preparation are performed by Arash Jalil Khab-bazi with insightful guidelines from Sunny Ri Li. The user-defined equationof state (EoS) codes for the binary mixture are developed in C by extendingthe standard codes in the ANSYS Fluent’s manual for pure substances. Thefollowing components have been published, presented, or submitted to reflectthe research findings dissemination:− Parts of Chapter 1 have been presented at the 13th International GreenEnergy Conference, 2021, titled: “The blending and transmission ofhydrogen and natural gas in transmission and distribution pipelines.”− Parts of Chapters 2 and 3 have been presented at the Canadian Soci-ety for Mechanical Engineering (CSME) International Congress, 2022,titled: “Green hydrogen supply to urban infrastructure and buildingsthrough blending into the existing grid.”- The Best Presentation of the Advanced Energy Systems Symposium.− Parts of Chapter 4 have been published as a proceeding and presentedat the Canadian Society for Mechanical Engineering (CSME) Interna-tional Congress, 2023, titled: “Blending of hydrogen into a natural gasvidistribution pipeline in British Columbia through a Tee junction forreducing GHG emissions.”- The Best Paper of the Congress.− Parts of Chapters 1 to 5 are submitted to the International Journalof Hydrogen Energy, 2023, titled: “Mixing hydrogen into natural gasdistribution pipeline system through Tee junctions.”Also, the user-defined EoS codes generated and used in this work aremade available for public use on GitHub1. For more information, pleasecontact the author through the email provided on GitHub.1GitHub: https://github.com/arashjkh?tab=repositoriesviiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . .xviiList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . xixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . .xxivDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xxviChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Previous Studies . . . . . . . . . . . . . . . . . . . . . . . . . 41.2.1 Thermophysical Properties . . . . . . . . . . . . . . . . 41.2.2 Steady-state Flow Assumption . . . . . . . . . . . . . . 81.2.3 Transient Flow Assumption . . . . . . . . . . . . . . . 10viiiTABLE OF CONTENTS1.2.4 Mixing and Injection . . . . . . . . . . . . . . . . . . . 131.3 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171.4 Objective and Organization of Thesis . . . . . . . . . . . . . . 17Chapter 2: Mathematical and Computational Details . . . . . 192.1 Flow Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.2 Quantifying the Mixing Homogeneity . . . . . . . . . . . . . . 222.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 232.4 Numerical Discretization and Validation . . . . . . . . . . . . 242.5 Computational Domain . . . . . . . . . . . . . . . . . . . . . . 28Chapter 3: Thermodynamic Relations . . . . . . . . . . . . . . 333.1 SRK Equation of State (EoS) Model . . . . . . . . . . . . . . 333.2 Thermodynamic Properties . . . . . . . . . . . . . . . . . . . 343.3 Mixture Thermodynamic Relations . . . . . . . . . . . . . . . 37Chapter 4: Numerical Study of the Distribution-pressure (DP)Pipeline . . . . . . . . . . . . . . . . . . . . . . . . . 394.1 Case Description . . . . . . . . . . . . . . . . . . . . . . . . . 394.2 Numerical Results and Discussion . . . . . . . . . . . . . . . . 42Chapter 5: Numerical Study of the Intermediate-pressure (IP)Pipeline . . . . . . . . . . . . . . . . . . . . . . . . . 475.1 Case Description . . . . . . . . . . . . . . . . . . . . . . . . . 475.2 Numerical Results and Discussion . . . . . . . . . . . . . . . . 495.2.1 Quantifying the Mixing Homogeneity . . . . . . . . . . 495.2.2 Comparing Ideal Gas and Real Gas EoS . . . . . . . . 53Chapter 6: Concluding Remarks . . . . . . . . . . . . . . . . . . 576.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 576.2 Key Outcomes . . . . . . . . . . . . . . . . . . . . . . . . . . . 576.3 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . 58ixTABLE OF CONTENTSBibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68Appendix A: On the Analysis of Vertical Low-pressure DistributionPipes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68A.1 Case Description . . . . . . . . . . . . . . . . . . . . . . . 68A.2 Analytical Methodology . . . . . . . . . . . . . . . . . . . 69A.3 Maximum Flow Rate Condition . . . . . . . . . . . . . . 73A.4 Minimum Flow Rate Condition . . . . . . . . . . . . . . 77A.5 Concluding Remarks . . . . . . . . . . . . . . . . . . . . 80Appendix B: Supplementary Data . . . . . . . . . . . . . . . . . . . 81Appendix C: Ideal Gas EoS Code . . . . . . . . . . . . . . . . . . . 89xList of TablesTable 2.1 Summary of the grid study. . . . . . . . . . . . . . . . 29Table 4.1 Mixing homogeneity length in terms of z/D for all cases. 44Table 5.1 Mixing homogeneity length in terms of z/D for all in-jection configurations using both ideal gas and SRKreal gas EoS. . . . . . . . . . . . . . . . . . . . . . . . . 54xiList of FiguresFigure 1.1 Hydrogen delivery pathways. . . . . . . . . . . . . . . 3Figure 1.2 Schematic of the Tee junction with the main flow ofnatural gas and the side flow of hydrogen. . . . . . . . 16Figure 2.1 (a) Schematic of the entire pipe; (b) A view of theTee junction. The notation χH2 represents the hydro-gen mole fraction. The boundary notations and gridstructure are shown. The origin of the coordinate sys-tem is located at the center of the branch. . . . . . . . 20Figure 2.2 A non-scaled view of the computational domain andthe grid near the outlet. . . . . . . . . . . . . . . . . . 25Figure 2.3 Pressure histories recorded in the simulation of Test 2using the current SRK EoS. . . . . . . . . . . . . . . . 26Figure 2.4 The validation exercise using different studies for pres-sure, p (MPa), versus decompression wave speed, W(m/s). The dashed lines correspond to the numeri-cal simulations conducted by Liu et. [1], the blacksolid line shows the experimental results conductedby Botros et al. [2], and the red solid line representsthe present work. . . . . . . . . . . . . . . . . . . . . . 27Figure 2.5 (a) 3-D view of the computational grid; (b) cross-sectional (z=0) view of the computational domain. . . 28xiiLIST OF FIGURESFigure 2.6 Top wall friction factor λ up to 50D from the injectionpoint along the primary pipe for five cases. For clarity,the length z is indicated in terms of the primary pipediameter D. . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 2.7 Top wall hydrogen mole fraction χH2 up to 300D fromthe branch center along the primary pipe for all cases. 31Figure 2.8 The centerline velocity magnitude, u (m/s), up to300D from the branch center along the primary pipefor all cases. . . . . . . . . . . . . . . . . . . . . . . . 32Figure 4.1 A view of the entire pipe with a length of 302D (a)and its Tee junction (b). . . . . . . . . . . . . . . . . . 40Figure 4.2 A view of hydrogen distribution at the branch center(z/D = 0) for secondary pipes of varying sizes (a) NPS3.0; (b) NPS 2.0; (c) NPS 1.25; (d) NPS 1.0; (e) NPS0.75; and (f) NPS 0.5. For clarity, the width x andheight y are indicated in terms of the primary pipediameter D. . . . . . . . . . . . . . . . . . . . . . . . . 41Figure 4.3 Cross-sectional view of the hydrogen mole fraction af-ter 0.5D from the injection point for secondary pipesof varying sizes (a) NPS 3.0; (b) NPS 2.0; (c) NPS1.25; (d) NPS 1.0; (e) NPS 0.75; and (f) NPS 0.5. . . . 42Figure 4.4 Cross-sectional view of hydrogen mole fraction after2D from the branch center for secondary pipes of vary-ing sizes (a) NPS 3.0; (b) NPS 2.0; (c) NPS 1.25; (d)NPS 1.0; (e) NPS 0.75; and (f) NPS 0.5. . . . . . . . . 43Figure 4.5 Cross-sectional view of hydrogen distribution after 10Dfrom the injection location for secondary pipes of vary-ing sizes (a) NPS 3.0; (b) NPS 2.0; (c) NPS 1.25; (d)NPS 1.0; (e) NPS 0.75; and (f) NPS 0.5. . . . . . . . . 44xiiiLIST OF FIGURESFigure 4.6 Cross-sectional view of hydrogen distribution after 159Dfrom the injection point for secondary pipes of varyingsizes (a) NPS 3.0; (b) NPS 2.0; (c) NPS 1.25; (d) NPS1.0; (e) NPS 0.75; and (f) NPS 0.5. . . . . . . . . . . . 45Figure 4.7 Variation in CoV along the cross-sections of the pri-mary pipe for all six cases. . . . . . . . . . . . . . . . 46Figure 5.1 Inclined view of hydrogen and methane flow in thepipeline for (a) vertical top-side, (b) horizontal, and(c) vertical bottom-side injection configurations. Themole fraction of hydrogen χH2 on cross-sectional slicesof the pipes is shown on the contour map. The slicesare evenly spaced at 0.5D from the branch center upto 10D. . . . . . . . . . . . . . . . . . . . . . . . . . . 48Figure 5.2 Cross-sectional view of hydrogen distribution after 0.5ddownstream of the injection point for (a) vertical top-side, (b) horizontal, and (c) vertical bottom-side injec-tion configurations. The contour map shows the molefraction of hydrogen χH2 . The width x and height yare expressed in terms of the primary pipe diameterD for clarity. . . . . . . . . . . . . . . . . . . . . . . . 50Figure 5.3 The variation of hydrogen mole fraction χH2 along thetop, center, bottom, left, and right lines within the pri-mary pipe, and CoV of hydrogen mole fraction on eachz-slice along the primary pipe for (a) vertical top-side,(b) horizontal, and (c) vertical bottom-side injectionscenarios. . . . . . . . . . . . . . . . . . . . . . . . . . 51Figure 5.4 Cross-sectional view of hydrogen distribution after 0.5Ddownstream of the injection point for (a) vertical top-side, (b) horizontal, and (c) vertical bottom-side in-jection configurations. . . . . . . . . . . . . . . . . . . 52xivLIST OF FIGURESFigure 5.5 Cross-sectional representation of hydrogen mole frac-tion after 2D downstream of injection location for (a)vertical top-side, (b) horizontal, and (c) vertical bottom-side injection designs. . . . . . . . . . . . . . . . . . . 52Figure 5.6 Cross-sectional illustration of hydrogen mole fractionafter 10D downstream of injection point for (a) verticaltop-side, (b) horizontal, and (c) vertical bottom-sideinjection cases. . . . . . . . . . . . . . . . . . . . . . . 53Figure 5.7 The CoV of hydrogen mole fraction on each cross-section along the primary pipe for (a) vertical top-side,(b) horizontal, and (c) vertical bottom-side injectioncases. The red lines indicate the SRK EoS scenario,whereas the blue lines reflect the ideal gas EoS sce-nario. The black dashed line represents the homoge-neous state, which has a Coefficient of Variation (CoV)of 1.5%. . . . . . . . . . . . . . . . . . . . . . . . . . . 54Figure 5.8 Cross-sectional display of hydrogen mole fraction af-ter 54D downstream of injection place for (a) verticaltop-side, (b) horizontal, and (c) vertical bottom-sideinjection instances. . . . . . . . . . . . . . . . . . . . . 55Figure 5.9 Cross-sectional display of hydrogen mole fraction after229D downstream of injection location for (a) verticaltop-side, (b) horizontal, and (c) vertical bottom-sideinjection scenarios. . . . . . . . . . . . . . . . . . . . . 55Figure 5.10 Cross-sectional illustration of hydrogen distributionafter 280D downstream of the branch center for (a)vertical top-side, (b) horizontal, and (c) vertical bottom-side injection cases. . . . . . . . . . . . . . . . . . . . 56Figure 6.1 (a) 2-D view of the whole computational grid; (b)Zoomed view of the computational domain. . . . . . . 73xvLIST OF FIGURESFigure 6.2 Hydrogen mole fraction along the inlet and outlet forboth cases of 5 mole% and 20 mole%. . . . . . . . . . 74Figure 6.3 (a) Hydrogen mole fraction contour for (a) 5 mole%case and (b) 20 mole% case. . . . . . . . . . . . . . . . 75Figure 6.4 (a) 2-D view of the whole computational grid; (b)Zoomed view of the computational domain. . . . . . . 77Figure 6.5 Hydrogen mole fraction along the inlet and outlet forboth cases of 5 mole% and 20 mole%. . . . . . . . . . 78Figure 6.6 (a) Hydrogen mole fraction contour for (a) 5 mole%case and (b) 20 mole% case. . . . . . . . . . . . . . . . 79Figure 6.7 A summary of the parameters used for the DP casewith a side pipe of NPS 3.0. . . . . . . . . . . . . . . . 82Figure 6.8 A summary of the parameters used for the DP casewith a side pipe of NPS 2.0. . . . . . . . . . . . . . . . 83Figure 6.9 A summary of the parameters used for the DP casewith a side pipe of NPS 1.25. . . . . . . . . . . . . . . 84Figure 6.10 A summary of the parameters used for the DP casewith a side pipe of NPS 1.0. . . . . . . . . . . . . . . . 85Figure 6.11 A summary of the parameters used for the DP casewith a side pipe of NPS 0.75. . . . . . . . . . . . . . . 86Figure 6.12 A summary of the parameters used for the DP casewith a side pipe of NPS 0.5. . . . . . . . . . . . . . . . 87Figure 6.13 A summary of the parameters used for the IP case. . . 88xviList of AbbreviationsAA area-averagedACV automatic control valveAGA8 American Gas Association Report No. 8B.C. British ColumbiaBL Berthelot and LorenzCFD computational fluid dynamicsCH4 methaneCO2 carbon dioxideCoV coefficient of variationDNN deep neural networkDP distribution-pressureEoS equation of stateFEM finite element methodGERG-2008 Groupe-Europeen-de-Recherches-Gazie`res 2008GH2 gaseous hydrogenGHG greenhouse gasH2 hydrogenH2NG hydrogen and natural gas mixtureH2S hydrogen sulfideHENG hydrogen-enriched natural gasHHV higher heating valueID inner diameterIP intermediate-pressureLH2 liquid hydrogenxviiLHV lower heating valueMAPE mean absolute percentage errorMOP maximum operating pressureN2 nitrogenNIST National Institute of Standards and TechnologyNPS nominal pipe sizeOD outer diameterOTEC ocean thermal energy conversionPR76 Peng-Robinson 1976PT pressure transducerPtG power-to-gasRCV rapid connect valveRK Redlich-KwongSch. scheduleSRK Soave-Redlich-KwongSTD standard wallTI turbulence intensityTP transmission-pressureUDRGM user-defined real gas modelUK United KingdomxviiiList of SymbolsParametersα(T ) Soave’s alpha function (Eq. 3.4)χ mole fractionχH2 hydrogen mole fractionχi species mole fractionm˙ mass flow raten˙ molar flow rateV˙ volumetric flow rateR pipe radius (Eq. 2.14)µ temperature-dependant dynamic viscosityµt turbulent viscosityω acentric factor (Eq. 3.5)χH2 average hydrogen mole fraction (Eq. 2.8)jm average mass flux (Eq. 2.14)R molar gas constantΦ viscous dissipation function (Eq. 2.3)ψi species property (Eq. 3.34)ψm mixture property (Eq. 3.34)ρ densityρideal ideal gas density (Eq. 3.11)ρreal real gas density (Eq. 3.11)τadv advection timexixτdiff diffusion timeτstf stratification timeε turbulence dissipation rate (Eq. 2.6)a energy parameterA pipeline cross-sectional areaai, aj species energy parameteraij binary energy parameteram mixture energy parametera0 Soave-Redlich-Kwong EoS energy constanta1, a2, a3 cubic real gas coefficientsb co-volume parameterbi, bj species co-volume parameterbij binary co-volume parameterbm mixture co-volume parameterc sound speed (Eq. 3.25)cp specific heat at constant pressurecdepp specific heat at constant pressure departure termcidealp ideal gas specific heat at constant pressureC1ε, C2ε, C3ε,C1, C2, σε, σkk − ε realizable model-specific constants (Eqs. 2.6 and2.7)cp,1, cp,2, cp,3,cp,4, cp,5ideal gas specific heat at constant pressure coefficients(Eq. 3.21)D primary (main) pipe diameterd secondary (side) pipe diameterDeff effective diffusivityDi,m mass diffusivity of the ith component in the mixture (Eq.2.5)Dt turbulent diffusivityF real gas property function (Eq. 3.12)F dep property departure function (Eq. 3.12)xxF ideal ideal gas property function (Eq. 3.12)fi body force per unit volume (the gravitational force act-ing on the flow) (Eq. 2.2)g gravitational forceGb the generation of turbulence kinetic energy due to buoy-ancy (Eqs. 2.7 and 2.6)Gk the production of turbulence kinetic energy (Eq. 2.6)h enthalpy (Eq. 3.13)H pipe heighthdep enthalpy departure termhideal ideal gas enthalpyIAA area-averaged turbulence intensity (Eq. 2.15)jm mass flux (Eq. 2.14)K thermal conductivityk turbulence kinetic energy (Eq. 2.6)kij, lij binary interaction parametersL characteristic lengthM molecular weight (molar mass)m, n fitting coefficients for fully-developed mass flux profile(Eq. 2.14)Mi species molecular weightMm mixture molecular weightn empirical substance characteristic constant (Eq. 3.5)patm standard pressure (Eq. 3.26)pcr critical pressureR mass-specific gas constantr radial distance from the inlet face center (Eq. 2.14)ReD pipe flow Reynolds number (Eq. 2.15)s entropy (Eq. 3.14)sdep entropy departure termxxisideal ideal gas entropyT temperatureTcr critical temperatureTr reduced temperature (Eq. 3.26)U flow velocityui, uj velocity componentv specific volumeW decompression wave speedxi, xj Cartesian spacial indexYi, Yj species mass fractionY1, Y2 mass fraction of species 1, mass fraction of species 2YM the contribution of the fluctuating dilatation in com-pressible turbulence to the overall dissipation rate (Eq.2.6)Z compressibility factorPopular Dimensionless NumbersSct turbulent Schmidt numberFr Froude numberPe Peclet numberRe Reynolds numberSc Schmidt numberSuperscripts¯ 1) averaged value, 2) molar gas propertydep departure termideal ideal gas termSubscriptsCH4 methane’sxxii1, 2 species 1, species 2 (for binary mixtures)AA area-averagedadv advectionatm atmosphericcr criticaldiff diffusioneff effectiveH2 hydrogen’si,m (of) the ith component (in) the mixturei, j 1) Einstein summation indices, 2) species (component)indicesij binary indexm mixture indexNG natural gas’sout outlet indexp 1) (at) constant pressure, 2) primary (main) inlet indexr reduceds secondary (side) inlet indexstf stratificationt turbulentxxiiiAcknowledgementsFirst and foremost, I would like to express my gratitude and appreciationto my wonderful supervisor, Dr. Sunny Ri Li, for all their support andguidance during the course of this project. Sunny’s unwavering belief in mehas been an endless source of motivation for me throughout this journey,regardless of the challenges I have encountered.I acknowledge the support of the Natural Sciences and Engineering Re-search Council of Canada (NSERC), project number ALLRP 557088 – 20(Cette recherche a ´et´e financ´ee par le Conseil de recherches en sciencesna- turelles et en g´enie du Canada (CRSNG), num´ero de projet ALLRP557088 – 20). I would also like to acknowledge the Digital Research Allianceof Canada for granting us access to computational resources for carrying outnumerical modelings.I am also honored to have had Dr. Joshua Brinkerhoff and Dr. SinaKheirkhah as my committee members. Their useful advice during the annualreport and informative talks about my future aspirations have proven to bereally beneficial.I want to give special recognition to our esteemed industrial partner,FortisBC, for their amazing work in advancing renewable gas, especially hy-drogen, infrastructure in British Columbia, Canada. Gratitude particularlygoes to the Renewable Gas Supply team, John Quinn, Vincent Chou, andMatthew Hill, for all their support and guidance during the course of thisresearch project.Furthermore, I am grateful for all the support I had from my friend,Mahdi Azar Barenji, during difficult times. I am also thankful for my dearestxxivcolleague and friend, Mojtaba Zabihi, who provided me with support andencouragement during all the challenging moments.Last but not least, my heartfelt thanks go to my devoted mother, Elmira,my father, Atabak, and my lovely sister, Aylin, for their endless love andencouragement. They are the strongest, nicest, and loveliest individuals Ihave ever known, and I hope to make them proud one day by fulfilling mygoals.xxvDedicationAli, my beloved late grandfather, whose loving glance from above inspiresme with unbounded pride and motivation.Shalaleh, my wonderful cousin, a true embodiment of resilience and achieve-ment and a nature fanatic. On loving memory of Shalaleh and her dear par-ents, who left us far too young, their loss has left an everlasting void in ourhearts.xxviChapter 1Introduction1.1 BackgroundMoving toward a prosperous economy and cleaner climate demands build-ing on existing measures to significantly lower greenhouse gas (GHG) emis-sions in the long term. Canada, alongside many other countries, is committedto achieving net-zero GHG emissions by 2050 [3]. Several climate plans havebeen legislated to date. For instance, the 2030 emissions reduction plan out-lines that the government will meet its enhanced Paris Agreement target toreduce emissions by 40-45% from 2005 levels by 2030 [4]. The importance ofrenewable and low-carbon gases such as hydrogen and biomethane is unde-niable in this regard. According to a recent joint study in British Columbia(B.C.) [5], replacing conventional natural gas with made-in-B.C. renewableand low-carbon gases plays a crucial role in achieving CleanBC targets forefficient and affordable decarbonization of the existing infrastructures [6].An ideal hydrogen economy in the future entails more development of itstechnical aspects, including production, transmission, distribution, and stor-age [7]. Hydrogen, besides its great potential for reducing GHG emissions,also plays an important role in the electricity sector. The intermittent na-ture of renewable energy sources such as wind and solar in power generationcould endanger the stability and reliability of power systems [8]. A high-capacity storage and distribution system for electrical energy is, therefore,necessary [9]. Given this, the surplus electricity generated from renewableenergy sources is converted to green hydrogen through water electrolysis,called Power-to-Gas (PtG) technology [10]. The advantage of the green hy-11.1. BACKGROUNDdrogen alternative is that its synthesis from water creates just oxygen as abyproduct, and its recombination with oxygen from the air produces wateronce more when electricity is required, like in fuel cells [11]. The ”renewable”or ”green” terms for hydrogen stand for its production method, which is fromrenewable energy sources. For instance, an ocean thermal energy conversion(OTEC) process can be used to harness the solar energy stored in oceansand produce hydrogen consequently [12]. Extensive research pertaining tohydrogen production routes has been carried out in [13, 14].Hydrogen generated from PtG or other methods has to be transportedto end users as gaseous hydrogen (GH2) or liquid hydrogen (LH2). GH2 andLH2 trailer delivery, GH2 and LH2 rail delivery, LH2 cargo ships, and GH2Pipeline delivery are the most well-established pathways (see Figure 1.1).Yang and Ogden [15] demonstrated that pipeline delivery surpasses GH2 andLH2 trucks, where demand is relatively high for hydrogen. Additionally,pipeline transmission is considered to be the only feasible option for large-scale GH2 transport between distant locations [16]. Therefore, it is plausibleto imagine dedicated hydrogen infrastructure in the long term, supported byrising production and consumption rates. While creating dedicated pipelinesfor hydrogen parallel to the gas pipelines seems feasible, lack of space is chal-lenging [17]. Moreover, the availability of current natural gas infrastructureprovides an economical and secure opportunity to integrate hydrogen [18].As such, the most cost-effective, efficient solution during the current transi-tion period is using the existing natural gas infrastructure, in which the twogaseous substances are transported as one mixture called H2NG or HENG(hydrogen-enriched natural gas). Hydrogen and external renewable carbonsources could also be merged to generate renewable methane through themethanation process [19], which is an alternative for blending hydrogen intothe existing grids [20].Due to the different thermophysical properties of hydrogen and naturalgas, the H2NG mixture changes significantly in its transmission conditions.21.1. BACKGROUNDFig. 1.1 – Hydrogen delivery pathways.These changes also vary from low pressures to high pressures. An in-depthunderstanding of the process feasibility, associated hazards, risks, and safetymeasures is, thus, crucial [21]. That is to say, the high contents of hydrogencould give rise to multiple issues, including but not limited to grid integrity,safety, leakage, and the necessity of modifying or altering pipelines, storagetanks, compressors, or even end-user devices [22]. Assuming high H2 con-centrations, for instance, the distance between compressor stations shouldincrease to prevent the over-pressure [23]. Besides, the scarcity of empiri-cal studies on the impact of methane on the susceptibility of pipeline steelto hydrogen embrittlement entails a better understanding of material dura-bility [24]. In this regard, preventing major issues calls for a technicallyreasonable amount of injected hydrogen. This hydrogen percentage varies31.2. PREVIOUS STUDIESacross different regions, such as the United States, United Kingdom, Europe,or Canada, owing to different natural gas compositions as well as govern-mental acts and provincial laws. A case-by-case assessment, therefore, isinevitable [25]. In B.C., Canada, a range of 5-20 mol% H2 will be consid-ered in the distribution pipeline system. It is worth noting that the naturalgas pipeline infrastructure consists of transmission and distribution systems.The distribution system falls into two sub-categories: distribution-pressure(DP) and intermediate-pressure (IP) pipelines. The DP system operateswithin a pressure range of 2.5 kPag to 700 kPag, while the IP system runswithin a pressure range of 701 kPag to 2070 kPag. On the other hand, atransmission-pressure (TP) pipeline within the transmission system operateswithin a pressure range of 2071 kPag to 14893 kPag.1.2 Previous StudiesNumerical and analytical techniques play an important role in modelingpipeline flow phenomena owing to the high expenses and hazards of conduct-ing corroborative experimental tests. Throughout the literature, a variety oftopics has encouraged many since the outset of 21st century, including butnot limited to thermophysical properties, steady and unsteady flow modeling,mixing and injection, to name a few.1.2.1 Thermophysical PropertiesHydrogen and natural gas differ significantly in terms of their thermo-physical characteristics. This dissimilarity between the two gases affects theoverall thermophysical properties and energetic aspects of the H2NG mix-ture. Since natural gas compositions in the pipelines may differ, the poten-tial hydrogen concentration in the natural gas grid must be evaluated andtested precisely for different gas pipelines and natural gas cases before theintroduction of hydrogen.41.2. PREVIOUS STUDIESThe density of hydrogen is approximately nine times lower than that ofnatural gas, and its higher heating value (HHV) is one-third that of naturalgas. Accordingly, Haeseldonckx and D’haeseleer [17] demonstrated that byconsidering an H2NG mixture flow rate of three times higher than that ofnatural gas, it is quite possible to transport the approximately same amountof energy via H2NG mixture compared to natural gas. On the other hand,they investigated the different amounts of injected hydrogen with respect tochanges in the Wobbe index indicator. The Wobbe index, which determinesthe interchangeability of the fuel gases, is a crucial factor in combustionprocesses [26]. Depending on the type of natural gas, the Wobbe indexhas the lowest values in the range of 65-85 vol% hydrogen. As a result,they reported that hydrogen concentrations should not exceed 17 vol% sincegreater concentrations can cause a variety of technical and chemical issues inpipes.In a Dutch case study, Schouten et al. [27] showed that to transport thesame amount of energy in the presence of 9 mol% hydrogen in the mixture,a higher flow rate of 9% is inevitable. As such, a higher flow rate entailsan increase in pipeline pressure. They also discovered that injecting 10%hydrogen into rich gas, which has a higher Wobbe index than hydrogen,lowers the Wobbe index of the H2NG mixture. Surprisingly, this is also truefor lean gas, despite the fact that hydrogen has a higher Wobbe index thanlean gas. The data also indicated that adding 75-80 mol% hydrogen reducesthe Wobbe index of the H2NG mixture. It is also worth noting that theterms lean gas and rich gas are used for natural gases with different amountsof hydrocarbon richness, and therefore, different HHVs [26].Lehner et al. [10] claimed that hydrogen injection into the natural gaspipeline infrastructure gives rise to several questions. They presented a briefreview of recent studies about natural gas grid transport capacity and theeffects of injection on gas characteristics (e.g., HHV, LHV, and Wobbe in-dex). The findings suggested that the admissible concentration of hydrogen is51.2. PREVIOUS STUDIESclosely correlated with natural gas quality in the pipeline. As a result, a 5-15vol% hydrogen blend is achievable. Furthermore, when hydrogen is injected,both the HHV and the Wobbe index decrease. It was also demonstratedthat injecting up to 10% hydrogen reduces energy flow by 5-6%. To compen-sate for this reduction, the flow rate of H2NG must be raised, resulting in apressure drop and more energy consumption in compressors.In a Polish case study, Blacharski et al. [28] examined the variations ofpseudo-compressibility factor (Z), density, viscosity, lower and higher heat-ing values, the Wobbe indices, and transmitted energy amount per normalcubic meter. Since the maximum operating pressure in the Poland gas gridis 8.4 MPa, the mentioned parameters were studied for the pressure from101.325 kPa to 10 MPa and hydrogen concentration range of 0-30 mol%.According to their findings, increased hydrogen molar fraction and decreas-ing pipeline pressure both lead to an increase in the pseudo-compressibilityfactor. They established that increasing hydrogen content reduces gas densityby using both the Clapeyron equation and the derived results for the pseudo-compressibility factor. Additionally, hydrogen injection of up to 100% in apressure of 10 MPa results in the reduction of both heating values (HHV andLHV). They also reported that the reduction in density compared to purenatural gas is negligible at low pressures (0.4%), while it is possible to havea 2.4% reduction at high pressures (10 MPa). According to earlier research,Wobbe indices decline within the hydrogen concentration range of 0-79 mol%but increase in the region of 79-100 mol% hydrogen, even if heating valuesdecrease in that range.The effects of blending 1-10 vol% hydrogen into existing natural gaspipelines on the thermophysical properties of the H2NG mixture were pre-sented by Deymi-Dashtebayaz et al. [29]. They discovered that the concen-tration of hydrogen in the H2NG mixture is inversely related to the higherand lower heating values (HHV and LHV) as well as the relative density ofthe H2NG mixture. Furthermore, the Wobbe index decreases with low hy-61.2. PREVIOUS STUDIESdrogen injections (1-10 vol%). They also determined that in the case of theH2NG mixture, turbo-compressors must spend more energy to compensatefor the effects of pressure drop in order to convey the same amount of en-ergy to end-users. The use of the American Gas Association Report No. 8(AGA8) equation of state (EoS) instead of the previously utilized empiricalEoS, consideration of other components of natural gas than methane in cal-culations, and evaluation of the effects of hydrogen injection on the criticalthermo-physical properties of H2NG were all innovative aspects of this work.A comprehensive description of equations of state will be provided later.Abd et al. [30] studied the impact of up to 10% hydrogen concentra-tion on pipeline flow and natural gas thermophysical characteristics. As aresult, adding 10% hydrogen to the pipeline reduces density, increases pres-sure losses, and decreases temperature drop (11.78%, 5.39%, and 6.14%,respectively). With the addition of hydrogen to the natural gas pipeline,the critical pressure of the mixture rises, the critical temperature falls, andthe bubble curve of the phase envelope approaches infinity. A major part ofthis work is its discovery of density change, in which injection of hydrogenup to 2 vol% enhances mixture density; nevertheless, injection of 2-10 vol%hydrogen decreases mixture density.Considering the aforementioned findings, the injection of hydrogen intoexisting natural gas pipelines is not self-evident. In reality, it requires athorough examination of the characteristics of natural gas pipelines as wellas the qualities of H2NG mixtures with various hydrogen levels. Thus, futurework on local natural gas pipelines must take into account all of these factorsin order to produce rational conclusions. Furthermore, because the amount ofpermissible hydrogen concentration varies by scenario, this must be assessedin future research to avoid potential risks or issues.71.2. PREVIOUS STUDIES1.2.2 Steady-state Flow AssumptionIn general, establishing an analytical solution for the differential equationof a steady-state flow is possible since the term of time does not exist in asteady-state flow equation. Several studies have studied steady-state H2NGflows in natural gas pipelines in recent years. Their primary goals havebeen to calculate pressure and temperature changes throughout the pipeline,determine mass flow rate and energy transport rate, and investigate com-pressors and decompressors. In most situations, a constant temperature wasassumed to calculate the pressure drop along the pipeline.In addition to their previous findings, Schouten et al. [27] investigatedthe effects of 10 and 25 mol% hydrogen addition on temperature changeat the injection site and pressure drop along a high-pressure natural gaspipeline. The results demonstrated that the pressure drop decreases whenthe same flow volume is transported. In contrast, if the same amount ofenergy is assumed to be delivered, the volume of the H2NG mixture willincrease, raising the pressure drop. An intriguing discovery is that as thetemperature rises, so will the pressure drop. This supports the idea that high-temperature flows require greater volumes to transfer comparable amounts ofenergy. They also demonstrated that when hydrogen is injected into naturalgas, the temperature reduces. This temperature drop will be increased asthe pressure in the pipeline rises. Using rich gas instead of lean gas alsocontributes to the temperature drop caused by injection. This outstandingresearch looked into a variety of issues of hydrogen injection into existingnatural gas pipelines. As such, it is regarded as a pioneering contribution tothis field.Tabkhi et al. [31] provided an analytical model for steady-state H2NGmixture flow in existing natural gas pipelines, assuming a pseudo-constantcompressibility factor between two points in the pipeline. As previouslystated, transferring the same amount of energy in the presence of hydrogennecessitates a larger H2NG mixture flow rate, resulting in increased pressure81.2. PREVIOUS STUDIESloss. The pressure gradient is an important metric in determining compres-sor energy usage. As a result, they used their steady-state mathematicalmodel to optimize energy use in a complex case involving multiple compres-sor stations and pipelines. According to the results of their optimization,the maximum acceptable hydrogen fraction in an existing gas pipeline is 6mass%.Wodek et al. [18] analyzed pressure drop, density decrease, and changesin compressibility factor and flow velocity over a 200-kilometer-long high-pressure pipeline for hydrogen molar fractions of 0% (pure methane), 10%,20%, and 100% (pure hydrogen). They used the Renouard equation, a com-monly used equation for high-pressure gas pipelines, to compute pressuredrops assuming a constant-average temperature for the H2NG mixture. Theycame to the conclusion that raising the concentration of hydrogen, increasingpipe diameter, and decreasing flow rate reduce pressure loss. According totheir observations, the compressibility factor increases when the pipeline andhydrogen are added. The density of a gas can be calculated using its tem-perature, pressure, and compressibility factor. As a result, the density dropcalculations revealed that adding hydrogen increases density drop along thepipeline. It is worth mentioning that the density drop in pure hydrogen is mi-nor due to small pressure losses. It was also demonstrated that the increasein flow velocity down the pipeline is inversely proportional to the concentra-tion of hydrogen. Furthermore, the flow increase in the case of pure hydrogenwas demonstrated to be minimal. The assumption of constant-average tem-perature in this work is notable since it may minimize the mistakes in theresults.Kuczynski et al. [32] calculated the thermophysical changes in the flow ofa steady-state H2NG mixture through high-pressure natural gas pipelines ofvarying diameters. They assumed constant outlet pressure and inlet temper-ature with a 15% hydrogen injection. They discovered that the presence ofhydrogen reduces pressure drop and increases temperature rise, which is con-91.2. PREVIOUS STUDIESsistent with previous studies. As a result, despite the negative contributionof hydrogen to the energy characteristics of the H2NG mixture, the mixturemay be carried over longer distances than pure natural gas. Because hydro-gen has a lower density than natural gas, transporting the same flow rates ofthe H2NG mixture necessitates larger pipeline dimensions. In addition, thepresence of hydrogen raises the compressibility factor, raises the isobaric heatcapacity, and lowers the Joule-Thomson effect. The Joule-Thomson effect isthe isenthalpic pressure drop coinciding with a temperature change. Unlikemany other real gases, hydrogen has a negative Joule-Thomson coefficient,in which the pressure reduction results in a temperature rise. This paper isnotable for providing realistic findings for both hydrogen and H2NG mixturescenarios and comparing the results to one another.In a nutshell, applying a steady-state assumption for flow appears rea-sonable in a number of situations where end-user demand is constant, ortransience is not significant enough. This is because transient flow modelingincurs extra costs in terms of time, energy, and expenses.1.2.3 Transient Flow AssumptionHydraulic transient occurs due to periodic alterations in end-user demandor sudden changes in the pipeline infrastructure, including crack propagationin pipelines or the immediate closure of valves. Because time is a componentin transient analysis, thorough research of numerical modeling is requiredto produce practical solutions for complex flow equations. Some researchersassumed isothermal flow along the pipeline, which necessitated the solutionof two mass and momentum equations. A non-isothermal assumption, on theother hand, necessitates the use of a third equation, the energy conservationequation, in order to derive reasonable responses for the temperature of theH2NG mixture as a function of time and pipeline length.Elaoud and Hadj-Ta¨ıeb [33] investigated the fluctuating pressure changesof H2NG flow due to a sudden closure in a downstream rapid connect valve101.2. PREVIOUS STUDIES(RCV) using two different numerical approaches. To prevent compressorsfrom failure due to sudden pressure elevations, they also implemented anupstream automatic control valve (ACV), which starts to close by sensingperturbations resulting from the propagation of compression waves by rapidclosure of RCV. They showed that the pressure increases significantly dueto a surplus mass flow into the pipeline during the closure period of ACV.It was also demonstrated that blending hydrogen into natural gas signifi-cantly increases this maximum pressure resulting from the sudden closure ofthe valve. Because such a transient behavior can result in high pressures,these findings can benefit researchers in considering the dynamic behavior ofpipeline pressure in addition to steady-state pressure in analyzing pipelinesafety, integrity, leakage, and material durability. The findings of this studycan also be used to determine the amount of hydrogen that can be injectedinto present gas pipelines while keeping the maximum allowable pressurein mind. Another noteworthy component of this paper is that the authorsused two distinct numerical methods to acquire the same results and confirmtheir conclusions. It is also worth mentioning that applying an expressionto calculate the average density of a mixture is unrealistic and can lead tomistakes.In a Polish case study, Uilhoorn [34] studied non-isothermal unsteadyH2NG flow in a high-pressure gas pipeline in the case of varying hydrogenfractions. They investigated the impacts of hydrogen injection on pressurechanges, temperature changes, and energetic aspects of the pipeline and com-pared their results in both steady-state and transient-state conditions withone another. They concluded that hydrogen addition reduces average pres-sure drop in the case of the same gas demand and increases it consideringthe same energy demand. On the other hand, increasing hydrogen concentra-tion contributes to elevating the average temperature drop in both scenarios.Furthermore, it was shown that hydrogen addition reduces linepack in termsof both gas and energy demand. As an interesting finding, it was shown111.2. PREVIOUS STUDIESthat despite pressure profiles, which are almost the same in both steady andtransient states, temperature profiles are quite different from one anotherbecause it takes time for the temperature to reach an equilibrium phase ina transient state. This work successfully demonstrated that an isothermalsteady-state flow assumption does not suffice for analyzing real-life end-userapplications since end-user demand generally fluctuates. Therefore, it is vitalto assume non-isothermal unsteady-state flows in cases where exact temper-ature predictions are required.Under certain conditions, an appropriate EoS can produce accurate phys-ical properties, particularly in high-pressure transient cases. Hafsi et al. [35]studied the unsteady flow of the H2NG mixture in a looped gas network usingvan der Waals EoS and linear mixing rules via finite element method (FEM).Based on their findings, increasing hydrogen concentration leads to greateroverpressure, resulting from pressure oscillations due to transient flow. Be-cause the steady-state pressure of the H2NG mixture (boundary conditionfor transient flow) is lower than that of natural gas, it can compensate forthe overpressure caused by hydrogen addition. Overpressure may be prob-lematic for pipelines in two scenarios: 1) when transient produces on theupstream, and 2) when steady-state pressure drop is reduced in a way thatit equals with that of natural gas (compressors usually do this). The authorsinvestigated the latter and demonstrated that the same steady-state pressuredrop for the H2NG mixture as for natural gas might endanger the pipeline.In other words, high internal pressure and pressure oscillations can respec-tively lead to pipeline failure and pipeline fatigue by exceeding its allowablecircumferential stress. They also showed that 30% hydrogen is allowable toensure a safe pipeline. Implementing a finite element method and the find-ings regarding pipeline durability immediately after calculations about fluidicaspects make this work unique.Using multiple EoS such as Groupe-Europeen-de-Recherches-Gazie`res 2008(GERG-2008) [36], Peng-Robinson 1976 (PR76) [37], and AGA8 [38], Liu et121.2. PREVIOUS STUDIESal. [39] proposed a computation fluid dynamics (CFD) model for predictingthe decompression wave speed of high-pressure H2NG mixtures and success-fully validated their findings with experimental tests conducted by Botros etal. [2]. In a similar study, the impacts of using different combining rules forthe H2NG mixture through PR76 EoS were explored [1].In summary, when demand in usage sites changes or a sudden changein the pipeline takes place, it is self-evident that a transient model must beused. An unstable model can be improved further by making isothermal ornon-isothermal assumptions. In general, using a non-isothermal assumptionof pipeline flow complicates mathematical computations far more than usingan isothermal assumption. The case study and available resources (time,energy, and money) for the research determine if it is worthwhile.1.2.4 Mixing and InjectionThe previous research studies mostly assume a homogeneous H2NG mix-ture with certain hydrogen percentages. However, considering a homoge-neous mixture right after the injection process is unrealistic. Depending onthe operating pressure of the pipeline, from low to high, a density differenceof eight to nine times between hydrogen and methane is evident. Therefore,due to the buoyant nature of hydrogen, it tends to rise and get trapped inthe vicinity of the top wall upon the injection [40], which calls into ques-tion the functionality, durability, and safety of the grid. That is to say, thecurrent steel pipelines are specifically designed for natural gas transmission,and hydrogen can embrittle carbon steel and permeate through it, which isnot safe [41]. Moreover, in most industries including paper pulp, pharmaceu-tical, polymers, chemical, and biochemical, efficient mixing is critical. Theefficacy of the mixing operations determines the quality of the final products[42]. Similarly, the current end-user devices require a uniform mixture ofgases with certain additions of hydrogen to achieve optimal efficiency. Forinstance, with a focus on the performance of domestic gas appliances, Sun131.2. PREVIOUS STUDIESet al. [43] found that the amount of added hydrogen should be below 23%.Given this, further research is required for hydrogen and natural gas mixingwithin the pipelines, known as in-line mixing.Among many in-line mixing options, Tee junctions and static mixers haveproven to be efficient for single-phase turbulent flows. Having said this, whenthe side jet has a high enough momentum, Tee junctions prove to be an ef-fective and economical option [44]. Tee junction is the subject of this studyfor various reasons, including high momentum in the side flow, its fundamen-tal nature, simplicity, cost-efficiency, and broad application. However, staticmixers are also frequently used in industry. passive static mixers are moredesirable to many industrial processes than mechanically agitated systemsbecause they consume less energy, as the mixer does not move while the fluidmoves.Numerous scholars have investigated the usage of different types of staticmixers in a variety of industrial applications in the literature. A comprehen-sive study on various static mixers can be found in [45]. Through the non-intrusive electrical resistance tomography (ERT) approach, Jegatheeswaranet al. [46] evaluated the mixing of non-Newtonian fluids in a chaotic SMXstatic mixer and showed that the static mixer aided in enhancing the rate ofdiffusion by increasing fluid stretching. A similar CFD work was carried outin [47]. In another study, Liu et al. [48] investigated the impacts of four dif-ferent types of mixers on the mixing process of natural gas with various gascomponents in the pipeline, with the goal of enhancing mixing uniformity.As such, the spiral vane static mixer had the greatest influence on mixing.Moreover, the effect of adding a helical static mixer in the gas network wasstudied by Kong et al. [49] and the effect of mixing unit number, torsionangle, and mixing unit distance on the was explored. Accordingly, the use ofthe gas mixer will enhance gas mixing while minimizing pressure loss in thepipe. Furthermore, Liu et al. [50] investigated the hydrogen and natural gasmixing in an SMX static mixer. They looked at how mixing elements, gas141.2. PREVIOUS STUDIESflow rate, intake pressure, and hydrogen volume fraction affected the resultsand concluded that four elements are a balanced choice for optimal pressureloss and homogeneous mixing.In the context of Tee junctions, Forney and Lee [51] published one of themost fundamental studies on the optimal pipeline mixing at a Tee junction.Eames et al. [40] performed a numerical modeling of pure hydrogen injec-tion into a horizontal natural gas pipeline at a Tee junction and showed thatthe injection rate, side pipe diameter, and injection location significantly af-fect the mixing process within transmission-pressure pipelines in the UnitedKingdom (UK). Lindner et al. [52] aimed to optimize gas mixing using CFDsimulations by altering the configuration of a main pipe and two inclinedside inlets. They discovered that increasing side pipe angles, especially ob-tuse angles, lowered the required pipe length for establishing a homogeneousgas mixing. Villuendas et al. [53] proposed several CFD studies to sim-ulate the H2NG mixture through a pipe setup, which includes a blendingstation, the pumping and injection procedure, and various pipeline geome-tries, in order to predict the flow characteristics that may emerge under newgas circumstances through high-pressure pipes. Su et al. [54] employed CFDsimulations to analyze the mixing dynamics of hydrogen and methane withina Tee junction in a natural gas pipeline. Their study focused on the influenceof multiple factors on the overall homogeneity of the mixing process. A deepneural network (DNN) model based on CFD simulation data was developedto quickly predict the coefficient of variation (CoV) of a hydrogen concen-tration under different operating circumstances. The schematic of the Teejunction for the current problem, i.e., the main flow of natural gas and theside flow of hydrogen, is depicted in Figure 1.2.151.2.PREVIOUSSTUDIESFig. 1.2 – Schematic of the Tee junction with the main flow of natural gas and the side flow ofhydrogen.161.3. MOTIVATION1.3 MotivationAs a general consensus, most of the works in the literature mainly fo-cus on studying transmission system pipelines with high operating pressure.However, investigations of distribution pipeline infrastructure are rare, if notnon-existent.Additionally, the previous works assume a homogeneous H2NG mixturewith certain percentages of hydrogen. However, a homogeneous mixtureassumption right after the injection process is unrealistic.In terms of predicting the mixing homogeneity length inside the pipeline,the ideal gas and real gas EoS could vary greatly depending on the operat-ing pressure. No research has yet explored or compared these instances inthe hydrogen mixing process within an intermediate-pressure pipeline in anatural gas distribution system.While earlier investigators looked into the injection setting from differentlocations, the precise technique to quantify this phenomenon using the CoVcontinues to be unknown.1.4 Objective and Organization of ThesisThe primary purpose of this Thesis is to numerically evaluate the hydro-gen injection into existing natural gas distribution pipeline systems, wherethe mixing homogeneity length is quantified by defining the coefficient of vari-ation (CoV). Furthermore, the research gaps identified in the earlier sectionare tackled successfully in this Thesis.In Chapter 2, the details for the turbulent multi-species steady-state flowmodel are provided. Furthermore, in order to quantify the mixing homogene-ity and provide a quantity for the intensity of segregation, the CoV is defined.Next, the governing equations that link the inlet and outlet boundaries areexpressed. The computational details and domain, the grid study, as well asthe validation exercise, are represented at the end.171.4. OBJECTIVE AND ORGANIZATION OF THESISIn Chapter 3, the Soave-Redlich-Kwong (SRK) EoS and associated ther-modynamic properties are defined. The mixture relations are introduced toextend the pure gas expressions to those of a mixture.In Chapter 4, ideal gas EoS for a binary mixture of hydrogen and methaneis employed, and the impact of changing the secondary pipe diameter on themixing homogeneity in DP pipelines is quantified via CoV initially.In Chapter 5, an IP pipeline case is proposed, and a comparison is madebetween ideal gas and real gas scenarios, where the real gas SRK EoS for themixture is developed in conjunction with the van der Waals linear mixingrules. The effect of vertical top-side, horizontal, and bottom-side injectionsis also investigated.In Chapter 6, this thesis is summarized and concluded. Several recom-mendations for the prospective work are provided at the conclusion.In Appendix A, a preliminary study for a vertical subdivision of a distribu-tion pipeline has been carried out. The analysis investigates the likelihood ofnatural gas and hydrogen stratification in low-pressure vertical pipes carryinggas to residential buildings at maximum and minimum flow rate conditions.In Appendix B, the supplementary data and a summary of the flow pa-rameters are embedded.Appendix C includes the user-defined ideal gas EoS codes for the binarymixture of hydrogen and methane. The section presents the C codes for theideal gas in conjunction with van der Waals linear mixing rules for the binarymixtures of hydrogen and methane. For more information regarding the realgas codes, it is recommended to reach out to the author directly throughemail, which is provided in the GitHub1.1GitHub: https://github.com/arashjkh?tab=repositories18Chapter 2Mathematical andComputational Details2.1 Flow ModelFigure 2.1 shows a physical problem with a primary (main) pipeline anda secondary (side) one through which hydrogen is injected. The steady anddynamic flow behavior within a pipeline is governed by conservation equa-tions. Thus, transport equations are required further to evaluate the mixingquality of non-reacting multi-species gases. The steady-state model withEinstein summation notation is shown below, where the continuity equationis according to [55], such that∂(ρui)∂xi= 0 (2.1)where ρ is density, ui is the velocity component, and xi is the Cartesianspacial index. The momentum equation is∂(ρuiuj)∂xj= − ∂p∂xi+∂∂xj[µ(∂ui∂xj+∂uj∂xi)]+ ρfi (2.2)where p is pressure, µ dynamic viscosity, and fi the body force per unitvolume (the gravitational force acting on the flow). The energy equation,192.1. FLOW MODELFig. 2.1 – (a) Schematic of the entire pipe; (b) A view of the Tee junc-tion. The notation χH2 represents the hydrogen mole fraction. Theboundary notations and grid structure are shown. The origin of thecoordinate system is located at the center of the branch.which is derived from the first law of thermodynamics, is∂(ρcpujT )∂xj=∂∂xj(K∂T∂xj)+ uj∂p∂xj+ µΦ (2.3)where T is temperature, cp the specific heat at constant pressure, K thermalconductivity, and Φ the viscous dissipation function.An equation of state is necessary to link density and other gas variables.The choice between ideal gas and real gas EoS is determined by the operatingpressure and temperature of the system. Therefore, the gas equation of stateisp = ρZRT (2.4)202.1. FLOW MODELwhere R is the mass-specific gas constant, and Z denotes the compressibilityfactor that measures the deviation of a real gas from ideal gas behavior,with Z = 1 representing the ideal gas condition. A detailed analysis of thecompressibility factor and real gas EoS is discussed later in Chapter 3.Gas mixing is a multi-species problem, hence an additional species trans-port equation is needed to solve for species mass fraction Yi at different places,expressed as [56], such that∂(ρujYi)∂xj=∂∂xj[(ρDi,m +µtSct)∂Yi∂xj](2.5)whereDi,m is the mass diffusivity of the ith component in the mixture (laminardiffusivity coefficient), which is a fluid property. The turbulent Schmidtnumber is Sct = µt/ρDt taken as 0.7, where µt is the turbulent viscosity,and Dt is the turbulent diffusivity, which is a turbulent flow characteristic.In turbulent flow conditions, Di,m is generally not important, as the effectivediffusivity Deff is dominated by Dt. It should be noted that for a multi-component mixture,∑ni=1 Yi is always unity. Therefore, solving Eq. 2.5 forone of the species is sufficient in binary mixtures, as Y1+Y2 = 1 relates bothto one another.Due to the high Reynolds numbers of the primary and secondary flows,turbulent mixing becomes a crucial flow characteristic as the pipe pressurerises. Various turbulence models have been implemented in the literature, de-pending on the application and computational resources. The k−εmodel [57]is particularly useful when the turbulent length scale is significantly largerthan the molecular length scale, as in pipe flows. Thus, two additional trans-port equations for the turbulence dissipation rate ε and its kinetic energy kmust be defined. In this paper, the realizable k−ε model [58] is implementedbecause it has demonstrated superior performance across a broad range offlows. The transport equations for the turbulence kinetic energy k and its212.2. QUANTIFYING THE MIXING HOMOGENEITYdissipation rate ε in k − ε realizable model are∂(ρk)∂t+∂(ρkuj)∂xj=∂∂xj[(µ+µtσk)∂k∂xj]+Gk +Gb − pε− YM(2.6)∂(ρε)∂t+∂(ρεuj)∂xj=∂∂xj[(µ+µtσε)∂ε∂xj]+ ρC1Sε− ρC2 ε2k +√vε+ C1εεkC3εGb(2.7)where Gk is the production of turbulence kinetic energy; Gb is the generationof turbulence kinetic energy due to buoyancy; YM is the contribution of thefluctuating dilatation in compressible turbulence to the overall dissipationrate; and the C1ε, C2ε, C3ε, C1, C2, σε, and σk are model-specific constants[58].2.2 Quantifying the Mixing HomogeneityThe coefficient of variation (CoV) is a quantity to measure the degree ofsegregation and is widely employed to evaluate the mixing homogeneity [59].It is calculated by dividing the standard deviation of the concentration databy the mean concentration measurements. The mole fraction is used insteadof mole concentration in this study as it gives a more reasonable measurefor relative proportions of the components regardless of the thermodynamicproperties of the mixture. Accordingly, the CoV is specified on every cross-section along the primary pipe such thatCoV =1χH2√1A∫A(χH2 − χH2)2 dA (2.8)where χH2 is the hydrogen mole fraction on various points of the cross-sections, χH2 the average hydrogen mole fraction on each cross-section, andthe cross-section area is given by A = πD2/4. Generally, a CoV between222.3. BOUNDARY CONDITIONS0.5% and 5% indicates a homogeneous mixture in industrial processes [44].To ensure a homogeneous mixture, a more stringent value of 1.5% is selectedfor the cases in this study.2.3 Boundary ConditionsThere are two inlets and one outlet in the system. The mass flow ratem˙ correlates to the molar flow rate n˙ through n˙ = m˙/M , where M is themolecular weight (molar mass). The total mole conservation and hydrogenspecies mole conservation within the system isn˙out = n˙p + n˙s (2.9)χH2,outn˙out = n˙s (2.10)where ’p’, ’s’, and ’out’ subscripts denote the primary inlet, secondary inlet,and outlet conditions, respectively; and χH2,out represents the average hydro-gen mole fraction at the outlet. Therefore, combining Eqs. 2.9 and 2.10 leadtoχH2,out =n˙sn˙p + n˙s(2.11)Rewriting Eq. 2.4 on a molar basis leads ton˙ =pV˙ZRT(2.12)where R is the molar gas constant and relates to the mass-specific gas con-stant R through R = R/M . The volumetric flow rate V˙ can be defined asV˙ = m˙/ρ. Combining 2.11 and 2.12 givesχH2,out =11 +pppsZsZpTsTpV˙pV˙s(2.13)232.4. NUMERICAL DISCRETIZATION AND VALIDATIONwhich relates inlet properties and the average hydrogen mole fraction at theoutlet. The inlet pressure ratio pp/ps is negligible for short pipe entry lengths.Additionally, the compressibility factor ratio Zs/Zp can be neglected if idealgas conditions prevail.A power-law profile for the fully-developed mass flux is defined by con-sidering a radially constant density at both inlets such thatjmjm=(n + 1n)[1−( rR)m] 1n(2.14)where jm stands for the mass flux, jm the average mass flux, R the piperadius, r the radial distance from the inlet face center, and the fitting coeffi-cients m and n are selected based on the experimental profiles, set to 2 and7, respectively [60].Determining a suitable turbulence intensity (TI) is essential in highlyturbulent flows. But it is difficult to precisely set a TI profile, so an averagevalue is typically used instead. Basse [61] suggested that an area-averaged(AA) power-law fit has the best correlation with the experimental results in[62]. Therefore,IAA = 0.3173×ReD−0.1095 (2.15)where IAA is the area-averaged turbulence intensity and ReD the pipe flowReynolds number.2.4 Numerical Discretization and ValidationThe numerical simulation is performed using ANSYS Fluent 2022R1, acommercial CFD software package. The coupled method, which simultane-ously solves for both velocity and pressure fields, outperforms a segregatedmethod in robustness and accuracy, particularly in compressible flow cases.Thus, it is utilized here in conjunction with a second-order discretization forgradients, pressure, and diffusion terms, as well as a second-order upwind242.4. NUMERICAL DISCRETIZATION AND VALIDATIONscheme for the convection terms to ensure accuracy. The realizable k − εturbulence model with full buoyancy effects is used. The species transportmodel with no volumetric reactions is employed for tracking the species ofmethane and hydrogen.ANSYS Fluent includes real gas EoS models such as Redlich-Kwong(RK), PR76, and the National Institute of Standards and Technology (NIST)real gas models that employ REFPROP. However, the built-in EoS has lim-ited temperature ranges and poor processing efficiency when modeling gasmixtures [1]. The user-defined real gas model (UDRGM) in ANSYS Fluentallows users to customize real gas models to meet their specific modelingneeds. Accordingly, two UDRGM codes for ideal gas and SRK real gas bi-Fig. 2.2 – A non-scaled view of the computational domain and the gridnear the outlet.252.4. NUMERICAL DISCRETIZATION AND VALIDATIONnary mixture of hydrogen and methane were developed in C programminglanguage and compiled into ANSYS Fluent consequently.To validate the accuracy of SRK UDRGM code, the second shock tuberupture test (8.28 mole% of H2) for H2NG decompression at high pressuresperformed by Botros et al. [2] at TransCanada Pipeline’s Gas DynamicsTest Facility is examined. The shock tube has an inner diameter (ID) of 15inches and a length of 42 meters. Two CFD simulations using PR76 andGERG-2008 UDRGM carried out by Liu et al. [1] are also compared withthe current research and the experimental study. The computational domainis shown in Figure 2.2, and details and steps are chosen according to [1].Botros et al. [2] and Liu et al. [1] used pressure curves measured atFig. 2.3 – Pressure histories recorded in the simulation of Test 2 usingthe current SRK EoS.262.4. NUMERICAL DISCRETIZATION AND VALIDATIONthe first five pressure transducer (PT) locations (PT1, PT1A, PT1B, PT2,and PT3, with a distance of 29.5 mm, 92.4 mm, 102.8 mm, 200 mm, and350 mm, respectively) to calculate the decompression wave speed W . In theexperiment, pressure transducers were used to monitor the dynamic pressurefluctuations at these locations. Pressure changes are tracked for up to 3 msin this study, as shown in Figure 2.3. Given this, Figure 2.4 indicates thatusing the current SRK code, the overall trend for decompression wave speedversus pressure is consistent with previous findings. However, the curve isslightly shifted towards the right and bottom compared to the other results.Fig. 2.4 – The validation exercise using different studies for pressure, p(MPa), versus decompression wave speed, W (m/s). The dashed linescorrespond to the numerical simulations conducted by Liu et. [1], theblack solid line shows the experimental results conducted by Botros etal. [2], and the red solid line represents the present work.272.5. COMPUTATIONAL DOMAINThe reason for this discrepancy lies in the fact that natural gas is composedof methane and other heavier hydrocarbons. The previous studies have con-sidered natural gas and hydrogen (8.28 mole%), but the current SRK codeonly takes into account the binary mixture of methane and hydrogen (also,8.28 mole%). This difference in mixture composition leads to slightly greaterdecompression wave speeds in the SRK code results, which is reasonable.2.5 Computational DomainFluent meshing software is used to construct a polyhedral grid throughoutthe domain, as shown in Figure 2.5. Polyhedral meshes can adapt to chang-ing flow conditions and adhere to complex geometry. They offer reasonableaccuracy while also drastically reducing computing costs. As a result, a gridstudy is undertaken for the base case with a length of 300D, a main pipeof NPS 6.0, and a secondary pipe of NPS 3.0 to test the sensitivity of theFig. 2.5 – (a) 3-D view of the computational grid; (b) cross-sectional(z=0) view of the computational domain.282.5. COMPUTATIONAL DOMAINTable 2.1 – Summary of the grid study.16.20M 8.63M 4.59M 2.36M 1.17MTop wall friction factor MAPE 0.981 1.369 1.9537 2.566 -Top wall H2 mole fraction MAPE 0.63 0.658 2.902 2.255 -Centerline velocity magnitude MAPE 0.279 0.381 0.476 0.611 -simulation results to the mesh resolution. ”NPS” refers to nominal pipe sizeand provides an approximation of inner diameter for small pipe sizes [63].All pipes in this work have been allocated a standard wall (STD) thicknessof 40s, where ”Sch.” represents Schedule, a pipe wall thickness standard.Five cases with various cell numbers are investigated in this manner. Thetop wall friction factor, the hydrogen mole fraction on the top surface, andthe velocity magnitude on the centerline are considered for the grid study.The mean absolute percentage error (MAPE) is used to measure the relativedifference between consecutive cases, and the findings are summarized inTable 2.1.As illustrated in Figure 2.6, the top wall friction factor, specifically for thedotted curves (1.17M and 2.36M elements), is inconsistent at first. However,the curves approach an asymptotic value after some distance from the injec-tion point. Thus, a denser mesh is utilized 1D upstream and 2D downstreamof the branch center to ensure the results are accurate (see Figure 2.5(a)).Furthermore, it can be observed from Figure 2.7 that the MAPE decreasesfor the hydrogen mole fraction on the top wall as the number of meshesincreases from 1.17M up to 16.20M. The results for the centerline velocitymagnitude can be observed from Figure 2.8. Clearly, the cases with 1.17Mand 2.36M cells are not performing well. Overall, the values between 8.63Mand 4.59M elements are within an acceptable range compared to previouscases. Therefore, the 4.59M cell case is picked for efficiency.292.5. COMPUTATIONAL DOMAINλFig. 2.6 – Top wall friction factor λ up to 50D from the injection pointalong the primary pipe for five cases. For clarity, the length z is indi-cated in terms of the primary pipe diameter D.302.5. COMPUTATIONAL DOMAINFig. 2.7 – Top wall hydrogen mole fraction χH2 up to 300D from thebranch center along the primary pipe for all cases.312.5. COMPUTATIONAL DOMAINFig. 2.8 – The centerline velocity magnitude, u (m/s), up to 300D fromthe branch center along the primary pipe for all cases.32Chapter 3Thermodynamic RelationsUnlike the ideal gas equation of state, a real gas EoS takes into accountthe size of the molecules and the inter-molecular forces. There are two maintypes of real gas EoS: 1) cubic EoS, including PR76 [37] and Soave-Redlich-Kwong (SRK) [64]; 2) EoS with complicated structures, such as GERG-2008[36]. Equations with complex structures offer more accurate results comparedto simple cubic EoS. However, applying a complex EoS is challenging, and thecalculation time is significantly higher [1]. When high-pressure models arerequired in the gas-processing industries, both SRK EoS and PR76 EoS arethe models of choice [65]. The SRK EoS is employed in this study accordingly.3.1 SRK Equation of State (EoS) ModelThe SRK EoS for a pure gas is given by [64], such thatp =RTv − b +av(v + b)(3.1)where v is the specific volume, a the energy parameter, and b the co-volumeparameter. For a pure gas, a and b area = 0.42747R2T 2crpcrα(T ) (3.2)b = 0.08664RTcrpcr(3.3)333.2. THERMODYNAMIC PROPERTIESwhere Tcr is the critical temperature, pcr is the critical pressure, and α(T ) isthe Soave’s alpha function represented asα(T ) =[1 + n(1−√TTcr)]2(3.4)where the empirical constant n is a substance characteristic expressed asn = 0.48 + 1.574ω − 0.176ω2 (3.5)where ω is the acentric factor, which is -0.217 and 0.011 for hydrogen andmethane, respectively.3.2 Thermodynamic PropertiesTo solve for the density, Eq. 3.1 can be expressed in the cubic form suchthatv3 + a1v2 + a2v + a3 = 0 (3.6)where a1, a2, and a3 area1 = −RTp(3.7)a2 = −(b2 +RTbp− ap)(3.8)a3 = −abp(3.9)Thus, the gas density isρ =Mv(3.10)Denoting the real gas density obtained through Eq. 3.1 as ρreal and ideal343.2. THERMODYNAMIC PROPERTIESgas density as ρideal, the compressibility factor Z isZ =ρidealρreal(3.11)Enthalpy, entropy, and specific heat are all calculated using the appropri-ate properties of an ideal gas and the departure functions [66]. Representingthe property departure function as F dep and the property function calculatedthrough ideal gas relations as F ideal, the real gas property function F isF = F ideal + F dep (3.12)The enthalpy h and entropy s are thereforeh = hideal + hdep (3.13)s = sideal + sdep (3.14)where the departure terms arehdep = pv −RT +Tdadt− abln(v + bv)(3.15)sdep = Rln(v − bv0)+dadTbln(v + bv)(3.16)wheredadT= − 1T[(1 + n)√aa0 − a] (3.17)with a0 = 0.42747R2T 2cr/pcr. The ideal gas terms arehideal =∫ TT0cidealp (T ) dT + h0(T0) (3.18)353.2. THERMODYNAMIC PROPERTIESsideal =∫ TT0cidealp (T )TdT +R lnp0p+ s0(T0, p0) (3.19)where the variable with 0 subscript represents the reference state conditions.To calculate the specific heat at constant pressure cp, Eq. 3.12 can beextended to cp such thatcp = cidealp + cdepp (3.20)The ideal gas specific heat at constant pressure cidealp iscidealp = cp,1 + cp,2T + cp,3T2 + cp,4T3 + cp,5T4 (3.21)where the constants from cp,1 to cp,5 are given in [67]. The specific heatdeparture term cdepp iscdepp =(∂hdep∂T)p(3.22)wherecdepp = p(∂v∂T)p−R +Td2adT 2bln(v + bv)−(TdadT− a) dvdTv(v + b)(3.23)whered2adT 2=a0n2TTcr(n+√aTcra0T)(3.24)The sound speed c isc = v√√√√√− cpcp(∂v∂p)T+ T(∂v∂T)p2 (3.25)The temperature-dependant dynamic viscosity µ is defined in [68], suchthatµ = 6.3× 10−7M1/2(pcr/patm)2/3Tcr1/6(T3/2rTr + 0.8)(3.26)363.3. MIXTURE THERMODYNAMIC RELATIONSwhere patm is the standard pressure 101.325 kPa, and the reduced tempera-ture is Tr = T/Tcr.The thermal conductivity K is expressed according to [69], such thatK = µ(cp +10.4M)(3.27)3.3 Mixture Thermodynamic RelationsMixing and combining rules are required to extend the EoS to gas mix-tures. The van der Waals one-fluid mixing principles are a popular methodfor a mixture containing p components. Thus, a and b in Eq. 3.1 are extendedto am and bm througham =p∑i=1p∑j=1YiYjaij (3.28)bm =p∑i=1p∑j=1YiYjbij (3.29)where Yi and Yj are mass fraction of a component in the mixture; am themixture energy parameter; and bm the mixture co-volume parameter. Thebinary energy parameter aij and co-volume parameter bij can be obtained viaclassical Berthelot [70] and Lorenz [71] (BL) combining rules, such thataij =√aiaj(1− kij) (3.30)bij =12(bi + bj)(1− lij) (3.31)where kij and lij are called binary interaction parameters. The former is setto 0.0156 for a hydrogen-methane binary substance (i ̸= j) and zero for purehydrogen or methane (i = j). The latter is non-zero only when studyingspecial or complex polar cases, thus it is set to zero (lij = 0) in this context.373.3. MIXTURE THERMODYNAMIC RELATIONSAs such, Eqs. 3.28 and 3.29 reduce toam =p∑i=1p∑j=1YiYj√aiaj(1− kij) (3.32)bm =p∑i=1Yibi (3.33)By applying Eqs. 3.32 and 3.33 to Eqs. 3.15 and 3.16, the departurefunctions can be extended to mixtures. Therefore, the mixture enthalpyhm and entropy sm can be obtained by evaluating the ideal gas terms forthe mixture. Using a linear mixing rule, the mixture and single componentproperties relate throughψm =p∑i=1Yiψi (3.34)where ψm can be replaced by Rm, cidealp,m , hidealm , sidealm , µm, and km and standsfor a mixture property, and ψi denotes that of a single component. Oncehm is calculated, the mixture specific heat cp,m can be determined using Eq.3.20.Additionally, the mole fraction of the ith component χi relates to themass fraction of the ith component Yi through χi = YiMm/Mi. Given this,The molecular weight M relates to the mass-specific gas constant R in Eq.2.4 as outlined in [72] through R = R/M , where R is the universal gasconstant.38Chapter 4Numerical Study of theDistribution-pressure (DP)Pipeline4.1 Case DescriptionLocated in the interior region of B.C., a distribution-pressure NPS 6.0pipeline exiting a gate station with an inner diameter D of 6.065 inches isevaluated (see Figure 4.1). The length L of the primary pipeline is 302D,which is roughly 46 meters long and connects to the first branch of thedownstream distribution network. As a result, it is a good representation ofthe maximal length L required for homogeneous blending.Six secondary pipes for the Tee junction ranging from NPS 0.5 to NPS3.0 are chosen to reflect actual pipe sizes routinely used in the industry. Theoperating distribution pressure in this study is 420 kPag, which is decreasedfrom a maximum transmission pressure of 6.62 MPag within the station. Theinlet pressures are nearly equivalent, with no pressure drop due to the shortentry lengths. The primary and secondary flows have temperatures of 276.5K and 288.15 K, respectively.The flow rate and other gas qualities vary throughout the year, fromlow summer to high winter consumption. As a result, an average mass flowrate of approximately 0.5 kg s-1 passes through the primary pipe. Thus, theReynolds number of the main flow is approximately 397,722, which is highlyturbulent. Consequently, an average hydrogen mole fraction of 20 mole% is394.1. CASE DESCRIPTIONFig. 4.1 – A view of the entire pipe with a length of 302D (a) and itsTee junction (b).selected to ensure the maximized hydrogen injection. Through Eq. 2.13, thesecondary pipe flow rate can be calculated using the average hydrogen molefraction in the system and the flow rate at the primary pipe. Given this,the Reynolds number of the side flow ranges from approximately 34750 to171405, which is again turbulent.Since the operating pressure of 420 kPag is relatively low, especially whencompared to the IP and TP pipelines, the Zs/Zp term can be neglected inEq. 2.13. In other words, the compressibility factor equals one. Therefore,the density contrast of almost eight times in ideal gas conditions originatesfrom the different molecular weights of methaneMCH4 = 16.04 kg kmol-1 andhydrogen MH2 = 2.02 kg kmol-1. A summary of the selected parameters anddata is embedded in Appendix B. That is to say, the parameters for the side404.1. CASE DESCRIPTIONpipes of NPS 3.0, NPS 2.0, NPS 1.25, NPS 1.0, NPS 0.75, and NPS 0.5 canbe observed in Figure 6.7, Figure 6.8, Figure 6.9, Figure 6.10, Figure 6.11,and Figure 6.12, respectively.Fig. 4.2 – A view of hydrogen distribution at the branch center (z/D = 0)for secondary pipes of varying sizes (a) NPS 3.0; (b) NPS 2.0; (c) NPS1.25; (d) NPS 1.0; (e) NPS 0.75; and (f) NPS 0.5. For clarity, the widthx and height y are indicated in terms of the primary pipe diameter D.414.2. NUMERICAL RESULTS AND DISCUSSION4.2 Numerical Results and DiscussionFigure 4.2 illustrates a cross-sectional view (z-slice) of hydrogen distribu-tion at the branch center (z/D = 0) for all cases. Considering a constant massflow rate, the jet penetration intensity increases as the secondary pipe diam-eter decreases. Moving half the primary pipe diameter (0.5D) downstream ofthe flow, as indicated in Figure 4.3, a low secondary flow momentum causesthe incoming flow to deflect when it comes to contact with the primary flow.When the hydrogen is blended, its velocity is greatly accelerated in theprimary flow direction, ultimately achieving the primary flow velocity. Be-cause of the significant density difference between hydrogen and methane,hydrogen is buoyant and becomes trapped towards the top wall. This isFig. 4.3 – Cross-sectional view of the hydrogen mole fraction after 0.5Dfrom the injection point for secondary pipes of varying sizes (a) NPS3.0; (b) NPS 2.0; (c) NPS 1.25; (d) NPS 1.0; (e) NPS 0.75; and (f) NPS0.5.424.2. NUMERICAL RESULTS AND DISCUSSIONevident from the first three cases in Figure 4.3(a-c), where the hydrogen istrapped near the top wall of the pipe. This additionally confirms the poormixing quality, which is comparable to an empty or coaxial tube, where sim-ple mixing occurs due to turbulence eddies that promote axial and radialmass interchange. In contrast, the mixing is expected to be enhanced signif-icantly in the last three cases (see Figure 4.3(d-f)), where the secondary flowpenetrates deeply into the primary flow.The progression of mixing after two times the primary pipe diameter (2D)downstream of the injection point is depicted in Figure 4.4, where it can beobserved that the hydrogen is diffusing only from the bottom side in the firstthree cases (a-c), while it diffuses from all sides in the last three cases (d-f).Furthermore, as illustrated in Figure 4.4(f), the side stream collides with thelower inner surface of the pipe, where hydrogen is diverted and flows up alongFig. 4.4 – Cross-sectional view of hydrogen mole fraction after 2D fromthe branch center for secondary pipes of varying sizes (a) NPS 3.0; (b)NPS 2.0; (c) NPS 1.25; (d) NPS 1.0; (e) NPS 0.75; and (f) NPS 0.5.434.2. NUMERICAL RESULTS AND DISCUSSIONthe inner wall of the pipe.Figure 4.5 shows the mixing progression ten times the primary pipe di-ameter (10D) downstream of the injection center. The last three cases (Fig-ure 4.5(d-f)) have nearly achieved homogeneity. However, the first threecases (Figure 4.5(a-c)) continue to display the dark blue region for hydrogenmole fraction, or zero mole% hydrogen, near the bottom wall. Table 4.1summarizes the mixing homogeneity length for all scenarios.Fig. 4.5 – Cross-sectional view of hydrogen distribution after 10D fromthe injection location for secondary pipes of varying sizes (a) NPS 3.0;(b) NPS 2.0; (c) NPS 1.25; (d) NPS 1.0; (e) NPS 0.75; and (f) NPS0.5.Table 4.1 – Mixing homogeneity length in terms of z/Dfor all cases.NPS 3.0 NPS 2.0 NPS 1.25 NPS 1.0 NPS 0.75 NPS 0.5z/D 159 154 144 39 52 34444.2. NUMERICAL RESULTS AND DISCUSSIONIn Figure 4.6, the mixing progression is illustrated after 159D from theinjection point, which is the fully mixed state for the widest secondary pipecase (NPS 3.0). The uniform blue tint in all of them indicates that all of theinstances are entirely mixed.Figure 4.7 compares the CoV versus the dimensionless length z/D for ev-ery case and verifies the prior findings. That is, the curves grow considerablysteeper as the secondary pipe diameter decreases, showing that the mixinghomogeneity length diminishes. In Figure 4.7, there is a large disparity be-tween the NPS 1.0 and NPS 1.25 cases, indicating that the three pipes withthe smallest diameters achieve homogeneity in much shorter lengths. Thisobservation conforms with earlier interpretations. Due to buoyancy forces,hydrogen remains along the top wall in the first three scenarios, resulting ina stratified flow of methane and hydrogen for some distance after the injec-Fig. 4.6 – Cross-sectional view of hydrogen distribution after 159D fromthe injection point for secondary pipes of varying sizes (a) NPS 3.0; (b)NPS 2.0; (c) NPS 1.25; (d) NPS 1.0; (e) NPS 0.75; and (f) NPS 0.5.454.2. NUMERICAL RESULTS AND DISCUSSIONtion. As such, turbulent diffusion is vital in driving the mixing process inthe latter stages. However, the secondary flow momentum flux is substan-tially increased in the last three cases, where buoyancy forces tend to elevatethe hydrogen flow simultaneously, which significantly increases the mixingprocess.Based on the results, the secondary stream momentum must be suffi-ciently high to achieve distributive mixing. However, it is critical to keep itfrom becoming too high, as this would lead to the incoming jet hitting the farwall, which would result in backmixing. Furthermore, increasing the velocityof the hydrogen can result in increased noise levels, which necessitates the useof practical velocity limits for injections in above-grade pipelines. Accord-ingly, a maximum velocity of 30 m/s is allowed, which FortisBC recommendsfor modeling hydrogen injection in distribution lines. As a result, it is criticalto balance the demand for efficient mixing with practical considerations.Fig. 4.7 – Variation in CoV along the cross-sections of the primary pipefor all six cases.46Chapter 5Numerical Study of theIntermediate-pressure (IP)Pipeline5.1 Case DescriptionThe intermediate-pressure pipeline within the distribution system is lo-cated in the coastal region of B.C., which is studied to illustrate the realgas effects and quantify the mixing homogeneity length for different injec-tion configurations (see Figure 5.1). The infrastructure includes an NPS 30pipeline with an inner diameter D of 29.25 inches, which is operating at 1200kPag with a potential future maximum operating pressure (MOP) of 2070kPag. Here, MOP is considered to more accurately capture the distinctionbetween ideal gas and real gas scenarios. The line is supplied through a singlefeed gate station from an upstream high-pressure transmission pipeline. TheNPS 30 IP pipeline extends 20 kilometers, where its pressure is reduced atmultiple points along its route and fed into the lower pressure 420 kPag (60psig) DP network that supplies gas to the downstream distribution networkserving residential, commercial, and industrial customers. In order to ensurea homogeneous mixing is achieved, a pipeline with a length L of 402D isconsidered here. Because this system will operate with a maximum of fivemole% hydrogen, a secondary pipe of NPS 3 for the Tee junction is proposed.Similar to the previous case, the pressure drop at the inlets is negligible, asthe entry lengths are considered short. The temperature of the primary and475.1. CASE DESCRIPTIONsecondary flows is 288.15 K and 313.15 K, respectively.The large seasonal range of the mass flow through the line highlights theneed to study a specific average mass flow rate for the purpose of settingup the CFD model. Accordingly, an average mass flow rate of 14.78 kg s-1flows within the primary pipeline, which is measured based on the availabledaily data of FortisBC for 2022-2023. Based on this flow rate, the Reynoldsnumber is about 2,353,911, which is very turbulent. The fluid within themain pipe is considered to be pure methane, whereas the one injected is acombination of methane and hydrogen with an 85 mole% hydrogen content.Having the primary mass flow rate and the average hydrogen mole fraction,which is five mole% for this case, the flow rate in the secondary pipe canFig. 5.1 – Inclined view of hydrogen and methane flow in the pipelinefor (a) vertical top-side, (b) horizontal, and (c) vertical bottom-sideinjection configurations. The mole fraction of hydrogen χH2 on cross-sectional slices of the pipes is shown on the contour map. The slices areevenly spaced at 0.5D from the branch center up to 10D.485.2. NUMERICAL RESULTS AND DISCUSSIONbe found by Eq. 2.13. Similarly, the Reynolds number of the side flow isapproximately 389,046, which is turbulent.As opposed to the prior scenario, when the pressure was low, the pres-sure of 2.07 MPag in this case is considerably high. As a result, employinga real gas EoS rather than an ideal one is inevitable. In other words, theZs/Zp term in Eq. 2.13 is no longer negligible. A significant density contrastof up to nine times can be observed depending on the specific EoS consid-ered. More information about this occurrence will be provided later. Also, amore detailed table of the parameters and data is provided in Appendix B,Figure 6.13.5.2 Numerical Results and DiscussionAccording to Figure 5.1, three Tee junction cases with vertical top-side,vertical bottom-side, and horizontal configurations are studied. Figure 5.2displays a cross-sectional (z-slice) view of the three injection configurationsat the edge of the Tee junction, i.e., half the secondary pipe diameter (0.5d)downstream of the branch center (z/D = 0). Previously, Eames et al. [40]demonstrated that the mixing is highest when injected from the bottom,allowing hydrogen to ascend through the pipe flow. Building upon it, in thisstudy, the mixing homogeneity length for various injection configurations isquantified using CoV.5.2.1 Quantifying the Mixing HomogeneityThe variation of hydrogen mole fraction along different lines (top, center,bottom, right, and left) within the primary pipe is shown in Figure 5.3. Foreach slice along the pipe, the coefficient of variation (CoV) of the hydrogenmole fraction is exhibited. A fully-mixed criterion of 1.5% is also definedto determine the level of mixing achieved. Figure 5.3(c) verifies the previ-ous findings, as the mixing is significantly improved in vertical bottom-side495.2. NUMERICAL RESULTS AND DISCUSSIONFig. 5.2 – Cross-sectional view of hydrogen distribution after 0.5d down-stream of the injection point for (a) vertical top-side, (b) horizontal,and (c) vertical bottom-side injection configurations. The contour mapshows the mole fraction of hydrogen χH2. The width x and height y areexpressed in terms of the primary pipe diameter D for clarity.injection. In fact, the mixing length for bottom-side injection is nearly fivetimes shorter than that for top-side injection (Figure 5.3(a)) and four timesshorter than that for horizontal injection (Figure 5.3(b)). Additionally, aninteresting fluctuating phenomenon for the hydrogen mole fraction along thecenter, left, and right lines within the pipe is observed in Figure 5.3(a andb). When hydrogen is injected, it tends to penetrate the primary flow but isdeflected and lifted toward the top wall. The raised hydrogen collides withthe top wall and is deflected once more. This incidence occurs for some dis-tances following the injection location and is eventually dampened when afully mixed state is attained.505.2.NUMERICALRESULTSANDDISCUSSION≥≤≥≤≥≤Fig. 5.3 – The variation of hydrogen mole fraction χH2 along the top, center, bottom, left, and rightlines within the primary pipe, and CoV of hydrogen mole fraction on each z-slice along the primarypipe for (a) vertical top-side, (b) horizontal, and (c) vertical bottom-side injection scenarios.515.2. NUMERICAL RESULTS AND DISCUSSIONThe mixing evolution can be better understood by observing the cross-sectional contours of the hydrogen mole fraction at different regions. Accord-ing to Figure 5.4, moving 0.5D downstream of the injection point, hydrogentends to rise towards the top wall right after injection takes place. Whenmoving 2D downstream of the injection center, this phenomenon can be out-lined more clearly, as demonstrated in Figure 5.5.Moving away from the branch center by 10D, as shown in Figure 5.6, itis clear that the hydrogen remains trapped near the top wall for the verticalFig. 5.4 – Cross-sectional view of hydrogen distribution after 0.5D down-stream of the injection point for (a) vertical top-side, (b) horizontal, and(c) vertical bottom-side injection configurations.Fig. 5.5 – Cross-sectional representation of hydrogen mole fraction af-ter 2D downstream of injection location for (a) vertical top-side, (b)horizontal, and (c) vertical bottom-side injection designs.525.2. NUMERICAL RESULTS AND DISCUSSIONFig. 5.6 – Cross-sectional illustration of hydrogen mole fraction after10D downstream of injection point for (a) vertical top-side, (b) hori-zontal, and (c) vertical bottom-side injection cases.top-side design, as opposed to the vertical bottom-side injection.5.2.2 Comparing Ideal Gas and Real Gas EoSTo illustrate the difference between ideal gas and real gas EoS cases, allinjection configurations are considered. Table 5.1 summarizes the mixinghomogeneity length for the three injection configurations using ideal gas andSRK real gas EoS. According to Figure 5.7(a-c), it is evident that real gasEoS predicts a longer homogeneity length as opposed to ideal gas EoS. Aspreviously stated, depending on the pressure, the density contrast is eighttimes when using ideal gas EoS, but it is more when using a real gas EoS, suchas SRK EoS. As a result, employing the SRK EoS increases the buoyancyof hydrogen, resulting in a rapid upward movement and trapping toward thetop wall. As such, when the pressure within the pipeline system reachessubstantial levels, it is critical to use a real gas EoS, such as SRK, ratherthan an ideal gas EoS. This is required in order to model the behavior ofgases under these conditions accurately.535.2.NUMERICALRESULTSANDDISCUSSIONFig. 5.7 – The CoV of hydrogen mole fraction on each cross-section along the primary pipe for (a)vertical top-side, (b) horizontal, and (c) vertical bottom-side injection cases. The red lines indicatethe SRK EoS scenario, whereas the blue lines reflect the ideal gas EoS scenario. The black dashedline represents the homogeneous state, which has a Coefficient of Variation (CoV) of 1.5%.Table 5.1 – Mixing homogeneity length in terms of z/D for all injection configurations usingboth ideal gas and SRK real gas EoS.Top (ideal) Top (SRK) Horizontal (ideal) Horizontal (SRK) Bottom (ideal) Bottom (SRK)z/D 266 280 219 229 50 54545.2. NUMERICAL RESULTS AND DISCUSSIONFig. 5.8 – Cross-sectional display of hydrogen mole fraction after 54Ddownstream of injection place for (a) vertical top-side, (b) horizontal,and (c) vertical bottom-side injection instances.Fig. 5.9 – Cross-sectional display of hydrogen mole fraction after 229Ddownstream of injection location for (a) vertical top-side, (b) horizontal,and (c) vertical bottom-side injection scenarios.After calculating the mixing homogeneity lengths for various injectionscenarios, examining the cross-sectional contours of hydrogen distributionwould be useful. Given the relatively high pressure within the pipeline, usingreal gas EoS would make sense. As a result, the fully-mixed condition forthe vertical bottom-side injection scenario is illustrated by Figure 5.8. Thevertical top-side (a) and horizontal (b) injection configurations continue todisplay non-uniform hydrogen distribution.The hydrogen distribution contour is shown in Figure 5.9 at z/D = 229,555.2. NUMERICAL RESULTS AND DISCUSSIONFig. 5.10 – Cross-sectional illustration of hydrogen distribution after280D downstream of the branch center for (a) vertical top-side, (b)horizontal, and (c) vertical bottom-side injection cases.where the horizontal injection scenario reaches the uniformly mixed state.Finally, 280D downstream of the injection point, the vertical top-sideinjection case reaches homogeneity, with the hydrogen mole fraction contourbecoming fully uniform (see Figure 5.10).56Chapter 6Concluding Remarks6.1 SummaryIn this Thesis, we have undertaken a CFD analysis to study the blendingof hydrogen into existing natural gas distribution pipeline systems throughTee junctions. Two pipeline cases were considered, a distribution-pressureone (420 kPag) in the interior region of B.C. in Canada and an intermediate-pressure case (2070 kPag) within the coastal region of B.C. in Canada. Thesimulations were performed on ANSYS Fluent by incorporating user-definedideal gas and SRK real gas models for the binary mixture of hydrogen andmethane using van der Waals linear mixing rules, and the real gas code wasvalidated by comparing the previous experimental and numerical cases forhydrogen and natural gas decompression in high-pressure steel tubes.6.2 Key OutcomesIn the DP case, where the ideal gas behavior prevails, the effect of chang-ing the secondary pipe diameter on the mixing homogeneity length was quan-tified using CoV for hydrogen mole fraction. The study revealed that whenjet penetration is low, indicated by a great secondary pipe diameter anda constant flow rate, hydrogen stratification is significant due to buoyancyforces and gets trapped near the top wall. Thus, the mixing homogeneitylength and the secondary jet intensity are directly related. It was also shownthat the mixing homogeneity length gets shortened significantly if the sec-ondary flow penetrates deep into the primary flow without getting trapped576.3. FUTURE DIRECTIONSnear the top wall from the outset. This phenomenon occurs because thebuoyancy forces tend to elevate the hydrogen flow at the same time whilethe turbulent mixing is taking place, enhancing the overall mixing process.Moreover, the difference between using an ideal gas EoS and a real gasone in predicting the mixing homogeneity length for an IP case was observed.When the operating pressure increased, using a real gas EoS was inevitable.Accordingly, the SRK real gas model tended to predict a greater mixinglength as opposed to an ideal gas one because the density contrast betweenhydrogen and methane increases in this case. To illustrate better the previousfindings in the literature, the difference between injection configurations onthe mixing length is quantified through CoV. It was shown that verticalinjection from the bottom leads to a mixing homogeneity length of almostfour times shorter than that for horizontal injection and five times shorterthan that for vertical top-side injection. This clearly shows that injectingfrom the bottom surpasses other injection configurations, provided that theindustry standards are followed.6.3 Future DirectionsOur future study will primarily focus on the buoyancy impact of hydro-gen for other injection designs and pipe configurations, such as the verticaldistribution gas line versus the horizontal gas pipeline described here. Also,the possibility of using static mixers can be studied and compared with thecurrent scenarios and findings.58Bibliography[1] Xiong Liu, Guillaume Michal, Ajit Godbole, and Cheng Lu. Decom-pression modelling of natural gas-hydrogen mixtures using the Peng-Robinson equation of state. International Journal of Hydrogen Energy,46(29):15793–15806, apr 2021.[2] KK Botros, S Igi, and J Kondo. Measurements of decompression wavespeed in natural gas containing 2-8%(mole) hydrogen by a specializedshock tube. In International Pipeline Conference, volume 50275, pageV003T05A001. American Society of Mechanical Engineers, 2016.[3] Government of Canada. Canadian Net-Zero Emissions AccountabilityAct, jun 2021.[4] Environment Canada and Climate Change Canada. 2030 Emissions Re-duction Plan: Canada’s Next Steps for Clean Air and a Strong Economy.Technical report, 2022.[5] BC Bioenergy Network, FortisBC, and Province of British Columbia.B.C. Renewable and Low-carbon Gas Supply Potential Study. Technicalreport, jan 2022.[6] The Government of British Columbia. CleanBC: Roadmap to 2030.Technical report, 2021.[7] Joan M. Ogden. Prospects for building a hydrogen energy infrastructure.Annual Review of Energy and the Environment, 24:227–279, nov 1999.59[8] Meltem Peker, Ayse Selin Kocaman, and Bahar Y. Kara. Benefits oftransmission switching and energy storage in power systems with highrenewable energy penetration. Applied Energy, 228:1182–1197, oct 2018.[9] U. Bu¨nger, J. Michalski, F. Crotogino, and O. Kruck. Large-scale under-ground storage of hydrogen for the grid integration of renewable energyand other applications. In Compendium of Hydrogen Energy, pages 133–163. Elsevier, jan 2016.[10] Markus Lehner, Robert Tichler, Horst Steinmu¨ller, and Markus Koppe.Power-to-Gas: Technology and Business Models. Springer InternationalPublishing, 2014.[11] Xianguo Li. Principles of Fuel Cells. Taylor Francis, 2005.[12] Pouria Ahmadi, Ibrahim Dincer, and Marc A. Rosen. Multi-objectiveoptimization of an ocean thermal energy conversion system for hydrogenproduction. International Journal of Hydrogen Energy, 40(24):7601–7608, jun 2015.[13] Ibrahim Dincer and Canan Acar. Review and evaluation of hydrogenproduction methods for better sustainability. International Journal ofHydrogen Energy, 40(34):11094–11111, sep 2014.[14] J. D. Holladay, J. Hu, D. L. King, and Y. Wang. An overview of hydrogenproduction technologies. Catalysis Today, 139(4):244–260, jan 2009.[15] Christopher Yang and Joan Ogden. Determining the lowest-cost hy-drogen delivery mode. International Journal of Hydrogen Energy,32(2):268–286, feb 2007.[16] Tory Borsboom-Hanson, Shashank Reddy Patlolla, Omar E. Herrera,and Walter Me´rida. Point-to-point transportation: The economics of hy-drogen export. International Journal of Hydrogen Energy, 47(74):31541–31550, aug 2022.60[17] Dries Haeseldonckx and William D’haeseleer. The use of the natural-gas pipeline infrastructure for hydrogen transport in a changing marketstructure. International Journal of Hydrogen Energy, 32(10-11):1381–1386, jul 2007.[18] Tomasz W lodek, Mariusz Laciak, Katarzyna Kurowska, and LukaszWe¸grzyn. Thermodynamic analysis of hydrogen pipeline transportation– selected aspects. AGH Drilling, Oil, Gas, 33(2):379, 2016.[19] Manuel Go¨tz, Jonathan Lefebvre, Friedemann Mo¨rs, Amy McDanielKoch, Frank Graf, Siegfried Bajohr, Rainer Reimert, and Thomas Kolb.Renewable Power-to-Gas: A technological and economic review. Renew-able Energy, 85:1371–1390, jan 2016.[20] Luis M. Romeo, Marco Cavana, Manuel Bailera, Pierluigi Leone, Begon˜aPen˜a, and Pilar Lisbona. Non-stoichiometric methanation as strategyto overcome the limitations of green hydrogen injection into the naturalgas grid. Applied Energy, 309:118462, mar 2022.[21] Zine labidine Messaoudani, Fotis Rigas, Mahar Diana Binti Hamid, andChe Rosmani Che Hassan. Hazards, safety and knowledge gaps on hy-drogen transmission via natural gas grid: A critical review. InternationalJournal of Hydrogen Energy, 41(39):17511–17525, oct 2016.[22] M. Penev M.W. Melaina, O. Antonia. Blending hydrogen into naturalgas pipelines networks: A review of key issues. Technical Report March,National Renewable Energy Lab. (NREL), Golden, CO (United States),2013.[23] Jingxuan Liu, Lin Teng, Bin Liu, Peng Han, and Weidong Li. Analysisof hydrogen gas injection at various compositions in an existing naturalgas pipeline. Frontiers in Energy Research, 9:685079, 2021.61[24] Cuiwei Liu, Hongchao Yang, Cailin Wang, Huimin Zhang, Rui Ding,Lina Ai, Xin Fan, Rui Zhang, Xiusai Xu, Yuanxing Ning, Y. FrankCheng, and Yuxing Li. Effects of CH4 and CO on hydrogen embrittle-ment susceptibility of X80 pipeline steel in hydrogen blended naturalgas. International Journal of Hydrogen Energy, apr 2023.[25] Joan Ogden, Amy Myers Jaffe, Daniel Scheitrum, Zane McDonald, andMarshall Miller. Natural gas as a bridge to hydrogen transportationfuel: Insights from the literature. Energy Policy, 115:317–329, apr 2018.[26] James G. Speight. Handbook of Natural Gas Analysis. John Wiley Sons,Inc., Hoboken, NJ, USA, aug 2018.[27] J. A. Schouten, J. P.J. Michels, and R. Janssen-Van Rosmalen. Effectof H2-injection on the thermodynamic and transportation properties ofnatural gas. International Journal of Hydrogen Energy, 29(11):1173–1180, sep 2004.[28] Tomasz Blacharski, Piotr Janusz, Maciej Kaliski, and Lukasz Zabrzeski.The effect of hydrogen transported through gas pipelines on the perfor-mance of natural gas grid. AGH Drilling, Oil, Gas, 33(2):515, 2016.[29] Mahdi Deymi-Dashtebayaz, Amir Ebrahimi-Moghadam, Seyyed ImanPishbin, and Mahdi Pourramezan. Investigating the effect of hydrogeninjection on natural gas thermo-physical properties with various com-positions. Energy, 167:235–245, jan 2019.[30] Ammar Ali Abd, Samah Zaki Naji, Tye Ching Thian, and Mohd RosleeOthman. Evaluation of hydrogen concentration effect on the naturalgas properties and flow performance. International Journal of HydrogenEnergy, 46(1):974–983, jan 2021.[31] F. Tabkhi, C. Azzaro-Pantel, L. Pibouleau, and S. Domenech. A math-ematical framework for modelling and evaluating natural gas pipeline62networks under hydrogen injection. International Journal of HydrogenEnergy, 33(21):6222–6231, nov 2008.[32] Szymon Kuczynski, Mariusz Laciak, Andrzej Olijnyk, Adam Szurlej,and Tomasz W lodek. Thermodynamic and technical issues of hydro-gen and methane-hydrogen mixtures pipeline transmission. Energies,12(3):569, feb 2019.[33] Sami Elaoud and Ezzeddine Hadj-Ta¨ıeb. Transient flow in pipelines ofhigh-pressure hydrogen-natural gas mixtures. International Journal ofHydrogen Energy, 33(18):4824–4832, sep 2008.[34] F. E. Uilhoorn. Dynamic behaviour of non-isothermal compressible nat-ural gases mixed with hydrogen in pipelines. International Journal ofHydrogen Energy, 34(16):6722–6729, aug 2009.[35] Zahreddine Hafsi, Sami Elaoud, and Manoranjan Mishra. A computa-tional modelling of natural gas flow in looped network: Effect of up-stream hydrogen injection on the structural integrity of gas pipelines.Journal of Natural Gas Science and Engineering, 64:107–117, apr 2019.[36] O. Kunz and W. Wagner. The GERG-2008 wide-range equation ofstate for natural gases and other mixtures: An expansion of GERG-2004. Journal of Chemical and Engineering Data, 57(11):3032–3091,nov 2012.[37] Ding-Yu Peng and Donald B Robinson. A new two-constant equation ofstate. Industrial & Engineering Chemistry Fundamentals, 15(1):59–64,1976.[38] Kenneth E Starling and Jeffrey L Savidge. Compressibility factors ofnatural gas and other related hydrocarbon gases. American Gas Asso-ciation, Operating Section, 1992.63[39] Bin Liu, Xiong Liu, Cheng Lu, Ajit Godbole, Guillaume Michal, andLin Teng. Decompression of hydrogen—natural gas mixtures in high-pressure pipelines: CFD modelling using different equations of state.International Journal of Hydrogen Energy, 44(14):7428–7437, mar 2019.[40] I. Eames, M. Austin, and A. Wojcik. Injection of gaseous hydrogeninto a natural gas pipeline. International Journal of Hydrogen Energy,47(61):25745–25754, jun 2022.[41] L. Briottet, I. Moro, and P. Lemoine. Quantifying the hydrogen embrit-tlement of pipeline steels for safety considerations. International Journalof Hydrogen Energy, 37(22):17616–17623, nov 2012.[42] Paulo E. Arratia and Fernando J. Muzzio. Planar laser-induced fluo-rescence method for analysis of mixing in laminar flows. Industrial andEngineering Chemistry Research, 43(20):6557–6568, sep 2004.[43] Mengxiao Sun, Xiaomei Huang, Yelong Hu, and Shan Lyu. Effects onthe performance of domestic gas appliances operated on natural gasmixed with hydrogen. Energy, 244:122557, apr 2022.[44] Arthur W. Etchells and Chris F Meyer. Mixing in Pipelines. In Handbookof Industrial Mixing, pages 391–477. John Wiley Sons, Ltd, jan 2004.[45] R. K. Thakur, Ch Vial, K. D.P. Nigam, E. B. Nauman, and G. Djelveh.Static mixers in the process industries - a review, aug 2003.[46] Sinthuran Jegatheeswaran, Farhad Ein-Mozaffari, and Jiangning Wu.Efficient mixing of yield-pseudoplastic fluids at low Reynolds numbers inthe chaotic SMX static mixer. Chemical Engineering Journal, 317:215–231, jun 2017.[47] Sinthuran Jegatheeswaran, Farhad Ein-Mozaffari, and Jiangning Wu.Process intensification in a chaotic SMX static mixer to achieve an64energy-efficient mixing operation of non-newtonian fluids. Chemical En-gineering and Processing: Process Intensification, 124:1–10, feb 2018.[48] Qi Liu, Yang Liu, Shanji Li, Henan Xu, Guohua Liang, Chong Sun,and Jingyun Sun. Analysis of the static mixer effect on natural gasmixing process in a pipeline. Flow Measurement and Instrumentation,85:102146, jun 2022.[49] Mingmin Kong, Shuaiming Feng, Qi Xia, Chen Chen, Zhouxin Pan, andZengliang Gao. Investigation of mixing behavior of hydrogen blended tonatural gas in gas network. Sustainability, 13(8):4255, 2021.[50] Yongzhen Liu, Anas Rao, Fanhua Ma, Xigui Li, Jinhua Wang, andQiuhong Xiao. Investigation on mixing characteristics of hydrogen andnatural gas fuel based on SMX static mixer. Chemical EngineeringResearch and Design, 197:738–749, sep 2023.[51] L. J. Forney and H. C. Lee. Optimum dimensions for pipeline mixingat a T-junction. AIChE Journal, 28(6):980–987, nov 1982.[52] Gert Lindner, Sonja Schmelter, Regine Model, Andreas Nowak, VolkerEbert, and Markus Ba¨r. A Computational Fluid Dynamics Study on theGas Mixing Capabilities of a Multiple Inlet System. Journal of FluidsEngineering, Transactions of the ASME, 138(3), mar 2016.[53] Teresa Villuendas, Carlos Montan˜e´s, Antonio Go´mez, Alberto CerezoAlarco´n, M Dolores Storch de Gracia, and Javier Sa´nchez-La´ınez. Ad-vancing in the decarbonized future of natural gas transmission net-works through a cfd study. International Journal of Hydrogen Energy,47(35):15832–15844, 2022.[54] Yue Su, Jingfa Li, Wangyi Guo, Yanlin Zhao, Jianli Li, Jie Zhao, andYusheng Wang. Prediction of mixing uniformity of hydrogen injec-65tion in natural gas pipeline based on a deep learning model. Energies,15(22):8694, nov 2022.[55] Frank White. Viscous Fluid Flow. McGraw-Hill Mechanical Engineer-ing, 3rd edition, 2005.[56] Adrian Bejan. Convection Heat Transfer. John wiley & sons, 2013.[57] Brian Edward Launder and Dudley Brian Spalding. Lectures in Mathe-matical Models of Turbulence. Academic Press, 1972.[58] Tsan-Hsing Shih, William W Liou, Aamir Shabbir, Zhigang Yang, andJiang Zhu. A new k− ε eddy viscosity model for high reynolds numberturbulent flows. Computers & fluids, 24(3):227–238, 1995.[59] Alena Kukukova´, Benjamin Noe¨l, Suzanne M Kresta, and Joelle Aubin.Impact of sampling method and scale on the measurement of mixingand the coefficient of variance. AIChE journal, 54(12):3068–3083, 2008.[60] Amgad Salama. Velocity profile representation for fully developed tur-bulent flows in pipes: a modified power law. Fluids, 6(10):369, 2021.[61] Nils T Basse. Turbulence intensity scaling: A fugue. Fluids, 4(4):180,2019.[62] Marcus Hultmark, Margit Vallikivi, Sean Christoper Collison Bailey,and AJ Smits. Turbulent pipe flow at extreme reynolds numbers. Phys-ical review letters, 108(9):094501, 2012.[63] ASME. ASME B36.19M-2018: Stainless Steel Pipe. ASME, 2018.[64] Giorgio Soave. Equilibrium constants from a modified Redlich-Kwongequation of state. Chemical Engineering Science, 27(6):1197–1203, jun1972.66[65] Saeid Mokhatab, William A Poe, and John Y Mak. Handbook of naturalgas transmission and processing: principles and practices. Gulf profes-sional publishing, 2018.[66] Adrian Bejan. Advanced Engineering Thermodynamics. John Wiley &Sons, 2016.[67] Michael J Moran, Howard N Shapiro, Daisie D Boettner, and Mar-garet B Bailey. Fundamentals of Engineering Thermodynamics. JohnWiley & Sons, 2010.[68] Nicholas P Cheremisinoff. Fluid Flow Pocket Handbook. Gulf PublishingCompany, 1984.[69] ERG Eckert and RM Drake Jr. Analysis of Heat and Mass Transfer.McGraw-Hill Company, 1972.[70] D Berthelot and P Sacerdote. Sur le me´lange des gaz et la compressibilite´des me´langes gaseux. Compt. Rendus, 128:820–822, 1899.[71] Hendrik Antoon Lorentz. Ueber die anwendung des satzes vom virial inder kinetischen theorie der gase. Annalen der physik, 248(1):127–136,1881.[72] Richard Edwin Sonntag, Claus Borgnakke, Gordon John Van Wylen,and Steve Van Wyk. Fundamentals of Thermodynamics, volume 6. Wi-ley New York, 1998.[73] ANSYS Fluent 2020R2. 8.2.6. User-Defined Real Gas Models(UDRGM). 2020.67AppendicesAppendix A: On the Analysis of VerticalLow-pressure Distribution PipesThis study aims to understand the H2NGmixture behavior in low-pressurevertical distribution pipelines. Due to the significant density difference be-tween hydrogen and methane, it is critical to figure out whether hydrogenand methane get stratified (separated) in vertical flow conditions.A.1 Case DescriptionThe pipeline system in this analysis is classified as a Vertical Subdivisionwithin the engineering design group. A 45-story residential condominiumtower with a range of end-use appliances in each unit is considered. Theparticular tower is served by a 60 mm (NPS 2) carbon steel distributionmain at a distribution pressure (DP) of 420 kPag. The DP supply is thenmetered and reduced in pressure to 1.75 kPag at a meter assembly at thebase of the tower. A single main header at 1.75 kPag serves the entire tower,with each unit branching off to a meter assembly on each floor. Followingthe pressure reduction, the supply line is increased to 168mm (NPS 6) toaccount for the pressure reductions. On the 11th floor, the header decreasesto 114 mm (NPS 4). On the 23rd floor, the header reduces in size to 60 mm(NPS 2). It is important to note that the gauge pressure of the gas will beimpacted by the increase in elevation of each tower floor, with the top floorappliances 135 meters above the lower appliances.68In this study, only 11 floors of the tower are investigated. Therefore, anNPS 6 pipe with an inner diameter D of 0.154056 meters, with an approxi-mate vertical pipe length L of 33.6 meters is considered. With this said, thevertical pipe length is approximately 237D. The pipe thickness is accordingto the previous chapters, i.e., a standard wall thickness of STD or 40s, where”Sch.” represents Schedule, a pipe wall thickness standard. The operatingpressure of the system is 1.75kPag, which is classified as low pressure.Moreover, two dynamic flow cases are considered: 1) the maximum flowrate of 155 Sm3/hr (standard cubic meters per hour), and 2) a low amountof 1 Sm3/hr. Both hydrogen mole fractions of 5 % and a maximum of 20 %are studied better to understand the impact of hydrogen presence within thislow-pressure natural gas system. It is worth noting that since the pressureis very low (Z ≃ 1), the ideal gas EoS is used, and there is no need for realgas EoS such as the SRK one used in the previous chapters.A.2 Analytical MethodologyTo identify the impacts of different factors such as advection, diffusion,and buoyancy (the main driver of stratification), three factors should be com-pared and studied: 1) advection vs. diffusion, 2) advection vs. stratification,and 3) stratification vs. diffusion. With this said, it is critical to introduceseveral dimensionless numbers.The first one is the Peclet number Pe, the ratio of a physical quantityadvection with the flow over the diffusion rate of that quantity by a gradient.For a mass transfer or species problem, it is defined as the multiplicationof the Reynolds number Re by the Schmidt number Sc (Pe = Re × Sc).ThereforePe =ULDeff(6.1)where U is the flow velocity, L is the characteristic length. The Schmidt69number Sc is provided asSc =µρDeff(6.2)In turbulent flow conditions, the turbulent Schmidt number Sct is usedand is defined asSct =µtρDt(6.3)As stated previously, turbulent diffusion is the main driver of diffusion inturbulent conditions, and molecular diffusion is relatively negligible. Thus,Deff equals to Dt.The Peclet number Pr can be defined in temporal form considering theadvection and diffusion occurrence times as followsPe =L2/DeffL/U=τdiffτadv(6.4)where τdiff and τadv stand for diffusion time and advection time, respec-tively.In the current analysis, the Peclet number Pr can be extended toPe =UDDeff=D2/DeffD/U=τdiffτadv(6.5)where the characteristic length L is the pipe inner diameter D.Now, if Pe≫ 1, the diffusion time is significantly longer than the advec-tion time, meaning that the advection is the main driver of mass transfer, andvice versa. Since in the current analysis, the hydrogen and natural gas arealready pre-mixed, comparing the advection versus diffusion is not sufficient.Hence, further parameters should be introduced.The second important expression is the Froude number Fr which is givenas the ratio of the flow inertia to the external field (such as gravitationalforce). The Froude number Fr is defined as70Fr =U√gL(6.6)where g is the gravitational force. Similar to the Peclet number Pe, theFroud number Fr can also be written in temporal form, such thatFr =√L/gL/U=τstfτadv(6.7)where τstf is the required temporal interval for stratification which is adirect result of the significant density contrast of hydrogen and natural gas,i.e., the buoyancy effect.Unlike the previous formula (Eq. 6.5), where the characteristic length Lwas the pipe diameter D, here it is the pipe height H, as the flow is movingupwards along the vertical pipe. Therefore, the formula can be extended tothe current scenario byFr =U√2(ρCH4 − ρH2)gHρCH4 + ρH2=√(ρCH4 + ρH2)H2(ρCH4 − ρH2)gH/U=τstfτadv(6.8)where ρCH4 is the methane density, ρH2 is the hydrogen density, and H isthe pipe height, as previously mentioned.Eq. 6.8 provides a good measure to compare the advection versus thestratification. Accordingly, if Fr ≫ 1, the stratification is negligible sincethe advection is dominant. This can be verified through the ratio of thetime intervals, as τstf is significantly higher than τadv. On the other hand, ifFr ≪ 1, stratification is stronger than advection. Thus, further investiga-tion, especially for the diffusion versus stratification should be carried out.In light of this analysis, it will become evident whether diffusion is dominantor buoyancy, which is a direct result of the gravitational force.In order to compare diffusion versus stratification, the characteristic tem-71poral intervals of these two should be introduced. Therefore, the requireddiffusion time over the stratification time is defined as followsτdiffτstf=D2/Deff√(ρCH4 + ρH2)H2(ρCH4 − ρH2)g(6.9)Similar to previous time ratios, if τdiff/τstf ≫ 1, the diffusion time isgreater than that of stratification, resulting in a significant buoyancy forceeffect compared to diffusion. On the other hand, if τdiff/τstf ≪ 1, the diffusionis stronger than the stratification. As a final note, if τdiff/τstf ∼ 1, thediffusion and stratification times are almost equal, which renders the case toa stationary vertical gas column at equilibrium.Given these expressions, two long pipes of 237D in length with both highand low flow rates are studied in the following sections.72A.3 Maximum Flow Rate ConditionThis system applies a maximum flow rate of 155 Sm3/hr. This conditiontakes place when most of the residents are using the gas during the winterseason. Based on the flow rate, the Reynolds number Re is about 22,000,which necessitates using an appropriate turbulence model, due to the tur-bulent flow conditions. Owing to the low pressure within the system (1.75kPg), a k−ω SST model is implemented to better model the turbulence andthe velocity profile just next to the pipe wall. The flow rate of 155 Sm3/hrapproximates a mass flow rate of 0.0282875 kg/s within the system. By us-ing Eq. 3.10, the mass flow rate is applied as a turbulent mass flux at theinlet boundary (the bottom of the pipe), and m and n are chosen 2 and 7,respectively. Since the gravitational force is acting along the pipe and in theopposite direction, a 2-D axisymmetric solver is chosen within the ANSYSFluent software to reduce the computational time and costs significantly.For the grid generation, the 36.6-meter pipe is divided into 7500 parts withFig. 6.1 – (a) 2-D view of the whole computational grid; (b) Zoomedview of the computational domain.73no bias factor, and both the inlet and outlet are divided into 50 parts with abias factor of 25 to create an inflation layer (to capture the boundary layerbehavior better). Therefore, an optimal structured quadrilateral grid with375k elements is generated, which results in a Y+ of less than 1. Figure 6.1illustrates the selected grid.For a mixture of 5 mole% hydrogen, the flow velocity U is almost 2.3 m/s,and for 20 mole% hydrogen it is about 2.66 m/s. The mixture density ρ forthe former case is 0.66 kg/s, and for the latter is approximately 0.57 kg/s.The dynamic viscosity µ is calculated using Eq. 3.26, where the first case is0.0000107396 Pa.s and the second case is 0.0000106589 Pa.s.Fig. 6.2 – Hydrogen mole fraction along the inlet and outlet for bothcases of 5 mole% and 20 mole%.74Given these conditions, the Pe number, as defined in Eq. 6.5, for thefirst case is 15238, and for the second case is 15353. Since Pe≫ 1, hydrogenand methane will be transported upstream. This implies that advectiondominates diffusion in this case. This can also be verified by the characteristictemporal interval ratio in Eq. 6.5, where the diffusion time (characteristictemporal interval) is significantly longer than that of the advection, furthershowing that the advection prevails.Consequently, the Fr number for the former case is 0.31, and for the latterone is 0.36. Since the value is close to 1, we can conclude that advection andstratification have the same effect on the system. Therefore, stratification isunlikely in this case. In order to verify this, the effects of diffusion versusstratification should be analyzed.Therefore, by using Eq. 6.9, the relative strength of diffusion and stratifi-cation will be known. For the first scenario with 5 mole% hydrogen injection,the result is 65.94, and for 20 mole% of hydrogen case, the result is 57.35.According to these findings, it is evident that stratification is a little bit domi-nant, however, the value is not that significant to separate the two gases sinceFig. 6.3 – (a) Hydrogen mole fraction contour for (a) 5 mole% case and(b) 20 mole% case.75the gases are uniformly distributed and the advection is significant.According to Figure 6.2, it is clear that in a maximum flow rate of 155Sm3/hr, the hydrogen and methane flow are in a fully mixed condition,showing that both diffusion and advection dominate the buoyancy force.Figure 6.3 further shows that the color of the hydrogen mole fraction contouris unchanged near the outlet (11th floor).76A.4 Minimum Flow Rate ConditionIn this case, the minimum flow rate of 1 Sm3/hr (compared to 155Sm3/hr) is considered. This flow rate is relatively low and is evident whensome of the end-users are using gas, especially in summer times. The flowrate approximates a mass flow rate of 0.0001825 kg/s within the system.With this low flow rate, the advection is relatively small. Since the Reynoldsnumber Re is low (about 140), the flow is laminar, which entails implement-ing a laminar flow condition within the ANSYS Fluent software. Similar tothe previous study, the inlet mass flux (the bottom of the pipe) can be definedthrough Eq. 3.10. However, instead of n = 7, n = 1 is used, which standsfor a fully-developed laminar flow. The variable m is again 2. Similarly,the gravitational force is acting along the opposite of the pipe. The axisym-metric 2-D solver is selected on the ANSYS Fluent software to decrease thecomputational time and costs.For the grid generation, the 36.6-meter pipe is divided into 7500 partsFig. 6.4 – (a) 2-D view of the whole computational grid; (b) Zoomedview of the computational domain.77with no bias factor, and both the inlet and outlet are divided into 25 partswithout a bias factor. Thus, an optimal structured quadrilateral grid with187.5k elements is created. Figure 6.4 shows the chosen grid for this case.For a mixture of 5 mole% hydrogen, the flow velocity U is almost 0.0148m/s, and for 20 mole% hydrogen it is almost 0.017191253 m/s. The mixturedensity ρ for the former case is 0.66 kg/s, and for the latter is approximately0.57 kg/s. The dynamic viscosity µ is calculated using Eq. 3.26, where thefirst case is 0.0000107396 Pa.s and the second case is 0.0000106589 Pa.s.Given that the flow velocity is low, the advection is small. Consequently,the Fr number is also small for both cases. This implies that Fr ≪ 1 andFig. 6.5 – Hydrogen mole fraction along the inlet and outlet for bothcases of 5 mole% and 20 mole%.78advection is not very strong, but it is still there. Therefore, buoyancy force(stratification) versus diffusion should be compared and studied. As men-tioned, in a laminar flow case, molecular diffusion is dominant, and turbulentdiffusion is negligible. Subsequently, by applying Eq. 6.9, the ratio for bothcases is 22.81. As tackled previously, this ratio is not significant enough toresult in the separation of two gases, given that both gases are fully mixedprior to potential separations. In order to verify this analytical finding, anumerical case is modeled with the given conditions.According to Figure 6.5, it is obvious that the hydrogen mole fraction atthe outlet is uniform, which outlines that the flow is fully mixed. Havingno turbulent diffusion, the main driver of diffusion is molecular diffusion inlaminar flow cases. It can be observed that when two gases are in a fullymixed condition, the diffusion tends to maintain them together. To betternote the findings, Figure 6.6 demonstrates that the color of the hydrogenmole fraction contour is unchanged close to the outlet.Fig. 6.6 – (a) Hydrogen mole fraction contour for (a) 5 mole% case and(b) 20 mole% case.79A.5 Concluding RemarksThe two-part study in this chapter shows that hydrogen and methanewill not get stratified as long as the above-mentioned conditions are met.Following the results of this study, it can be concluded that when two gasesare already pre-mixed, there must be either a very long pipe or a very strongexternal force to separate the two gases. In addition, both the advection andthe molecular diffusion play an important role in keeping the gas moleculestogether and act against the buoyancy force which is the primary driver ofstratification.It is critical to understand the different factors such as diffusion, advec-tion, and external force (gravity) in a mixture of hydrogen and methanewith significant density contrast, especially in vertical scenarios where thegravitational force g tends to separate the two flows. It is important to un-derstand this in pipes with longer heights. That is to say, our analysis coversa pipeline with a length of 237D for pipes with dynamic flow conditions,however, it is likely to have stratification in pipes with even longer heights(for instance, tall towers) in a completely stagnant situation (even though itis not very realistic). The stratification might also happen when the gasesare not pre-mixed, i.e., hydrogen is injected through a Tee junction similarto the main chapters of this thesis. Therefore, this entails a comprehensiveanalysis, including both analytical and numerical investigations.80Appendix B: Supplementary DataA summary of the parameters used for Chapters 4 and 5 is provided inthis chapter. For the DP scenario, six cases with different side pipe diameterswere investigated. Accordingly, the parameters for the side pipes of NPS 3.0,NPS 2.0, NPS 1.25, NPS 1.0, NPS 0.75, and NPS 0.5 can be observed inFigure 6.7, Figure 6.8, Figure 6.9, Figure 6.10, Figure 6.11, and Figure 6.12,respectively.For the IP case, since only the side pipe configuration changes, all theflow parameters are the same, and therefore, only one figure suffices. Thecorresponding data can be observed in Figure 6.13.81Fig. 6.7 – A summary of the parameters used for the DP case with aside pipe of NPS 3.0.82Fig. 6.8 – A summary of the parameters used for the DP case with aside pipe of NPS 2.0.83Fig. 6.9 – A summary of the parameters used for the DP case with aside pipe of NPS 1.25.84Fig. 6.10 – A summary of the parameters used for the DP case with aside pipe of NPS 1.0.85Fig. 6.11 – A summary of the parameters used for the DP case with aside pipe of NPS 0.75.86Fig. 6.12 – A summary of the parameters used for the DP case with aside pipe of NPS 0.5.87Fig. 6.13 – A summary of the parameters used for the IP case.88Appendix C: Ideal Gas EoS CodeIn this section, I present a user-defined-real-gas-model (UDRGM) code,which includes the ideal gas equation of state in conjunction with linear vander Waals mixing rules for a binary mixture of hydrogen and methane. Thiscode was initially developed by ANSYS Inc. for a pure substance [73], and Ihave extended it to a binary mixture of hydrogen and methane.Additionally, I have developed a real gas SRK EoS code for binary mix-tures of methane and hydrogen. This code can predict the thermodynamicbehavior of these two gases with high precision. For more information aboutthe code, please send me an email through my GitHub1 page.The user-defined ideal gas EoS code generated and used in this work ismade available for public use on GitHub1 and is as follows:1 /*2 *sccs id: @(#) real_ideal.c 1.10 Copyright1900/11/09 ANSYS , Inc.3 */45 /*6 * Copyright 1988 -1998 ANSYS , Inc.7 * All Rights Reserved8 *9 * This is unpublished proprietary source code ofANSYS , Inc.10 * It is protected by U.S. copyright law as anunpublished work11 * and is furnished pursuant to a written licenseagreement. It12 * is considered by ANSYS , Inc. to be confidential and1GitHub: https://github.com/arashjkh?tab=repositories89may not be13 * used , copied , or disclosed to others except inaccordance with14 * the terms and conditions of the license agreement.15 */161718 /*19 * Warning !!! Including udf.h is for getting definitionsfor20 * ANSYS FLUENT constructs such as Domain. Youmust21 * NOT reference any ANSYS FLUENT globalsdirectly from22 * within this module nor link this against anyANSYS23 * FLUENT libs , doing so will cause dependencieson a24 * specific ANSYS FLUENT binary such as fl551.exe and25 * thus won’t be version independent.26 */272829 /* ALL VARIABLES ARE IN SI UNITS.30 *31 * pressure [Pa]32 * temperature [K]33 * molar_weight [kg/kmol]34 * specific_volume [m3/kg]35 * molar_volume [m3/kmol]36 * universal_gas_constant (molar_gas_constant) [J/kmol -90K]37 * mass_specific_gas_constant (J/kg-K)38 * enthalpy [J/kg]39 * entropy [J/kg-K]40 * specific_heat [J/kg-K]41 *42 */4344 #include "udf.h"45 #include "stdio.h"46 #include "ctype.h"47 #include "stdarg.h"4849 #if RP_DOUBLE50 #define SMLL 1.e-2051 #else52 #define SMLL 1.e-1053 #endif5455 #define NSPECIE_NAME 8056 #define RGASU UNIVERSAL_GAS_CONSTANT /* 8314.34 SI units: J/kmol/K */57 #define PI 3.1415926545859 /* Here input the number of species in the mixture */60 /* THIS IS A USER INPUT */61 #define n_specs 26263 /* OPTIONAL REFERENCE (OFFSET) VALUES FOR ENTHALPY ANDENTROPY */64 #define h_ref 0.065 #define s_ref 0.0916667 static int (* usersMessage)(const char *,...);68 static void (* usersError)(const char *,...);6970 /* REFERENCE STATE */71 #define P_ref 101325.072 #define T_ref 298.1573 // static double P_ref , T_ref;74 static char gas[n_specs ][ NSPECIE_NAME ];7576 /* static property parameters */77 static double cp[5][ n_specs ]; /* specific heatpolynomial coefficients */78 static double mw[n_specs ]; /* molecular weights */79 static double tcrit[n_specs ]; /* critical temperature */80 static double pcrit[n_specs ]; /* critical pressure */81 static double vcrit[n_specs ]; /* critical specificvolume */8283 /* Static variables associated with ideal gas model */84 static double rgas[n_specs], cp_inT_ref[n_specs],h_ideal_ref[n_specs ];8586 void Mw();87 void Cp_Parameters ();88 void Tcrit ();89 void Pcrit ();90 void Vcrit ();9192 double cp_ideal_gas(double temp , int i);93 double viscosity(double temp , int i);94 double thermal_conductivity(double temp , int i);929596 DEFINE_ON_DEMAND(I_do_nothing)97 {98 /*99 This is a dummy function100 must be included to allow for the use of the101 ANSYS FLUENT UDF compilation utility102 */103 }104105 /*******************************************************************/106 /* Mixture Functions*/107 /* These are the only functions called from ANSYS FLUENTCode */108 /*******************************************************************/109110 void MIXTURE_Setup(Domain *domain , cxboolean vapor_phase, char *specielist ,111 int (* messagefunc)(const char *format ,...),112 void (* errorfunc)(const char *format,...))113 {114 /* This function will be called from ANSYS FLUENTafter the115 UDF library has been loaded.116 User must enter the number of species in the mixture93117 and the name of the individual species.118 */119120 int i, j;121 usersMessage = messagefunc;122 usersError = errorfunc;123124 /*125 P_ref = ABS_P(RP_Get_Real ("reference -pressure "),op_pres);126 T_ref = 298.15;127128 Message0 ("\n MIXTURE_Setup: ideal gas equation ofState \n");129 Message0 ("\n MIXTURE_Setup: reference -temperature is%f \n", T_ref);130131 if (P_ref == 0.0){132 Message0 ("\n MIXTURE_Setup: reference -pressurewas not set by user \n");133 Message0 ("\n MIXTURE_Setup: setting reference -pressure to 101325 Pa \n");134 P_ref = 101325.0;135 }136 */137138 /*====================================================*/139 /* ========= User Input Section====================== */140 /*94====================================================*/141 /*142 Define Species name.143 DO NOT use space for naming species144 */145 (void)strcpy(gas[0], "H2");146 (void)strcpy(gas[1], "CH4");147148 /*====================================================*/149 /* ========= End Of User Input Section=============== */150 /*====================================================*/151152 Message0("\n MIXTURE_Setup: RealGas mixtureinitialization \n");153 Message0("\n MIXTURE_Setup: Number of Species = %d \n", n_specs);154 for (i=0; i<n_specs; ++i)155 {156 Message0("\n MIXTURE_Setup: Species [%d] = %s \n", i, gas[i]);157 }158159 /*160 concatenate species name into one string and sendback to fluent161 */95162 strcat(specielist , gas [0]);163 for (i=1; i<n_specs; ++i)164 {165 strcat(specielist , " ");166 strcat(specielist , gas[i]);167 }168169 /* initialize */170 Mw();171 Cp_Parameters ();172 Tcrit();173 Pcrit();174 Vcrit();175176 for (i=0; i<n_specs; ++i)177 {178 rgas[i] = RGASU/mw[i];179 h_ideal_ref[i] = T_ref *(cp[0][i] + T_ref *(0.5* cp[1][i] + T_ref *(1./3.* cp[2][i] + T_ref *(0.25*cp[3][i] + T_ref *0.2*cp[4][i]))));180 cp_inT_ref[i] = cp[0][i]*log(T_ref)+T_ref *(cp[1][i]+T_ref *(0.5* cp[2][i]+ T_ref *(1./3.* cp[3][i]+0.25* cp[4][i]* T_ref)));181 }182183 }184185 /*------------------------------------------------------------*/186 /* FUNCTION: MIXTURE_mw*/96187 /* Returns mixture molar mass*/188 /*------------------------------------------------------------*/189190 double MIXTURE_mw(double yi[])191 {192 double sum = 0.0;193 int i;194195 for (i=0; i<n_specs; ++i)196 sum += yi[i]/mw[i];197198 return (1.0/ sum); /* (kg/kmol) */199 }200201 /*------------------------------------------------------------*/202 /* FUNCTION: MIXTURE_rgas*/203 /* Returns mass -specific mixture gas constant*/204 /*------------------------------------------------------------*/205206 double MIXTURE_rgas(double yi[])207 {208 double rgas_m = 0.0;209 int i;97210211 for (i=0; i<n_specs; ++i)212 rgas_m += yi[i]*rgas[i];213214 return rgas_m; /* (J/kg/K) */215 }216217 /*------------------------------------------------------------*/218 /* FUNCTION: MIXTURE_density*/219 /* Returns density given T and P*/220 /*------------------------------------------------------------*/221222 double MIXTURE_density(cell_t cell , Thread *thread ,cxboolean vapor_phase , double temp , double P, doubleyi[])223 {224 return (P/( MIXTURE_rgas(yi)*temp)); /* (kg/m^3) */225 }226227 /*------------------------------------------------------------*/228 /* FUNCTION: MIXTURE_specific_heat*/229 /* Returns specific heat given T and rho*/98230 /*------------------------------------------------------------*/231232 double MIXTURE_specific_heat(cell_t cell , Thread *thread, double temp , double density , double P, double yi[])233 {234 double cp_ideal_gas = 0.0;235 int i;236237 for (i=0; i<n_specs; ++i)238 cp_ideal_gas += yi[i]*(cp[0][i] + temp*(cp[1][i]+ temp*(cp[2][i] + temp*(cp[3][i] + temp*cp[4][i]))));239240 return cp_ideal_gas; /* (J/kg/K) */241 }242243 /*------------------------------------------------------------*/244 /* FUNCTION: MIXTURE_enthalpy*/245 /* Returns specific enthalpy given T and rho*/246 /*------------------------------------------------------------*/247248 double MIXTURE_enthalpy(cell_t cell , Thread *thread ,double temp , double density , double P, double yi[])249 {99250 double h = 0.0;251 int i;252253 for (i=0; i<n_specs; ++i)254 h += yi[i]*( temp*(cp[0][i] + temp *(0.5* cp[1][i]+ temp *(1./3.* cp[2][i] + temp *(0.25* cp[3][i]+ temp *0.2* cp[4][i]))))255 - h_ideal_ref[i]);256257 return (h + h_ref); /* (J/kg) */258 }259260 /*------------------------------------------------------------*/261 /* FUNCTION: MIXTURE_cp_integral*/262 /* Returns entropy ’s Cp_integral part given T andrho */263 /*------------------------------------------------------------*/264265 double MIXTURE_cp_integral(double temp , double yi[])266 {267 double cp_int = 0.0;268 int i;269270 for (i=0; i<n_specs; ++i)271 cp_int += yi[i]*(cp[0][i]*log(temp)+temp*(cp[1][i]+temp *(0.5* cp[2][i]+temp *(1./3.* cp[3][i]+0.25* cp[4][i]*temp)))100272 - cp_inT_ref[i]);273274 return cp_int; /* (J/kg -K) */275 }276277 /*------------------------------------------------------------*/278 /* FUNCTION: MIXTURE_entropy*/279 /* Returns entropy given T and rho*/280 /*------------------------------------------------------------*/281282 double MIXTURE_entropy(cell_t cell , Thread *thread ,double temp , double density , double P, double yi[])283 {284 return (s_ref + MIXTURE_cp_integral(temp , yi) +MIXTURE_rgas(yi)*log(P_ref/P)); /* (J/kg/K) */285 }286287 /*------------------------------------------------------------*/288 /* FUNCTION: MIXTURE_speed_of_sound*/289 /* Returns s.o.s given T and rho*/290 /*------------------------------------------------------------101*/291292 double MIXTURE_speed_of_sound(cell_t cell , Thread *thread , double temp , double density , double P, doubleyi[])293 {294 double cp = MIXTURE_specific_heat(cell , thread , temp, density , P, yi);295 double rgas = MIXTURE_rgas(yi);296297 return (sqrt(rgas*temp*cp/(cp -rgas))); /* m/s */298 }299300 /*------------------------------------------------------------------------*/301 /* FUNCTION: MIXTURE_rho_t*/302 /* derivative of denisty w.r.t. temperature at constantpressure (drdT|p) */303 /*------------------------------------------------------------------------*/304305 double MIXTURE_rho_t(cell_t cell , Thread *thread , doubletemp , double density , double P, double yi[])306 {307 return (-density/temp); /* (kg/m^3/K) */308 }309310 /*------------------------------------------------------------------------102*/311 /* FUNCTION: MIXTURE_rho_p*/312 /* derivative of denisty w.r.t. pressure at constanttemperature (drdp|T) */313 /*------------------------------------------------------------------------*/314315 double MIXTURE_rho_p(cell_t cell , Thread *thread , doubletemp , double density , double P, double yi[])316 {317 double rgas = MIXTURE_rgas(yi);318 return (1./( rgas*temp)); /* (kg/m^3/Pa) */319 }320321 /*-------------------------------------------------------------------------*/322 /* FUNCTION: MIXTURE_enthalpy_t*/323 /* derivative of enthalpy w.r.t. temperature at constantpressure (dhdT|p) */324 /*-------------------------------------------------------------------------*/325326 double MIXTURE_enthalpy_t(cell_t cell , Thread *thread ,double Temp , double density , double P, double yi[])327 {328 return MIXTURE_specific_heat(cell , thread , Temp ,density , P, yi); /* J/(kg.K) */103329 }330331 /*-------------------------------------------------------------------------*/332 /* FUNCTION: MIXTURE_enthalpy_p*/333 /* derivative of enthalpy w.r.t. pressure at constanttemperature (dhdp|T) */334 /*-------------------------------------------------------------------------*/335336 double MIXTURE_enthalpy_p(cell_t cell , Thread *thread ,double Temp , double density , double P, double yi[])337 {338 /* general form dh/dp|T = (1/ rho)*[ 1 + (T/rho)*drho/dT|p] */339 /* but for ideal gas dh/dp = 0 */340 return 0.0; /* J/(kg.Pa) */341 }342343 /*------------------------------------------------------------*/344 /* FUNCTION: MIXTURE_viscosity*/345 /*------------------------------------------------------------*/346347 double MIXTURE_viscosity(cell_t cell , Thread *thread ,104double temp , double density , double P, double yi[])348 {349 double mu = 0.;350 int i;351352 for (i=0; i<n_specs; ++i)353 mu += yi[i]* viscosity(temp , i);354355 return mu; /* kg/(m.s) */356 }357358 /*------------------------------------------------------------*/359 /* FUNCTION: MIXTURE_thermal_conductivity*/360 /*------------------------------------------------------------*/361362 double MIXTURE_thermal_conductivity(cell_t cell , Thread*thread , double temp , double density , double P,double yi[])363 {364 double kt = 0.;365 int i;366367 for (i=0; i<n_specs; ++i)368 kt += yi[i]* thermal_conductivity(temp , i);369370 return kt; /* W/(m.K) */371 }105372373 /*******************************************************************/374 /* Species Property Definitions*/375 /*******************************************************************/376377 void Mw() /* molecular weight */378 { /* kg/kmol */379 mw[0] = 2.01594; /*H2*/380 mw[1] = 16.04303; /*CH4*/381 }382383 void Pcrit () /* critical pressure */384 { /* Pa */385 pcrit [0] = 1.293 e6; /*H2*/386 pcrit [1] = 4.599 e6; /*CH4*/387 }388389 void Tcrit () /* critical temperature */390 { /* K */391 tcrit [0] = 32.98; /*H2*/392 tcrit [1] = 190.56; /*CH4*/393 }394395 void Vcrit () /* critical specific volume */396 { /* m3/kg */397 vcrit [0] = 0.031846; /*H2*/398 vcrit [1] = 0.006146; /*CH4*/106399 }400401 void Cp_Parameters () /* coefficients of specific heatpolynomials */402 { /* J/kg/K */403 cp [0][0] = 13602.45; /*H2*/404 cp [1][0] = 3.402317;405 cp [2][0] = -0.003358423;406 cp [3][0] = -3.907953e-07;407 cp [4][0] = 1.705345e-09;408409 cp [0][1] = 403.5847; /*CH4*/410 cp [1][1] = 9.057335;411 cp [2][1] = -0.01442509;412 cp [3][1] = 1.580519e-05;413 cp [4][1] = -6.343051e-09;414 }415416 /*------------------------------------------------------------*/417 /* FUNCTION: cp_ideal_gas*/418 /* Returns ideal gas specific heat given T*/419 /*------------------------------------------------------------*/420421 double cp_ideal_gas(double temp , int i)422 {423 double cpi = (cp[0][i]+temp*(cp[1][i]+temp*(cp107[2][i]+temp*(cp[3][i]+temp*cp[4][i]))));424425 return cpi;426 }427428 /*******************************************************************/429 /* FUNCTION: viscosity*/430 /*******************************************************************/431432 double viscosity(double temp , int i)433 {434 double mu , tr , tc , pcatm;435436 tr = temp/tcrit[i];437 tc = tcrit[i];438 pcatm = pcrit[i]/101325.;439 mu = 6.3e-7* sqrt(mw[i])*pow(pcatm ,0.6666)/pow(tc,0.16666) *(pow(tr ,1.5) /(tr+0.8));440441 return mu; /* (kg/m/s) */442 }443444 /*******************************************************************/445 /* FUNCTION: thermal_conductivity*/108446 /*******************************************************************/447448 double thermal_conductivity(double temp , int i)449 {450 double cp , mu;451452 cp = cp_ideal_gas(temp , i);453 mu = viscosity(temp , i);454455 return (cp + 1.25* rgas[i])*mu; /* W/m/K */456 }457458 /*******************************************************************/459 /* Mixture Functions Structure*/460 /*******************************************************************/461462 UDF_EXPORT RGAS_Functions RealGasFunctionList =463 {464 MIXTURE_Setup , /* initialize*/465 MIXTURE_density , /* density*/466 MIXTURE_enthalpy , /* sensibleenthalpy */467 MIXTURE_entropy , /* entropy109*/468 MIXTURE_specific_heat , /* specific_heat*/469 MIXTURE_mw , /*molecular_weight */470 MIXTURE_speed_of_sound , /* speed_of_sound*/471 MIXTURE_viscosity , /* viscosity*/472 MIXTURE_thermal_conductivity , /*thermal_conductivity */473 MIXTURE_rho_t , /* drho/dT |constp */474 MIXTURE_rho_p , /* drho/dp |constT */475 MIXTURE_enthalpy_t , /* dh/dT |const p*/476 MIXTURE_enthalpy_p , /* dh/dp |const T*/477 };110
Thesis/Dissertation
eng
Mechanical Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Li, Ri
Mixing gaseous hydrogen into natural gas distribution pipelines
Text