THEORETICAL AND COMPUTATIONAL STUDIES OF A LIQUID ENTERING A NETWORK OF PORES By Mohamed Alhaddad A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering-Doctor of Philosophy 2025 ABSTRACT Modeling the rate of fluid penetration into capillaries due to surface tension forces is often based on the Poiseuille flow solution. However, this model does not apply to short capillaries due to non-fully developed conditions at the entrance and exit regions. Improved models are needed for small capillary systems, which are crucial in processes such as oil droplet removal from water using thin membranes. Previous research has addressed deviations from Poiseuille flow near the entrance and moving meniscus, including the use of momentum conservation equations and inertia forces in kinetic models for infinite flow entering capillary tubes. Some studies have considered finite reservoir infiltration, assuming parallel flow lines, but neglected local acceleration due to inertia and gravity effects. This study presents a novel analysis focusing on the dynamic behavior of droplets in pores. It models a finite flow reservoir associated with a droplet and includes drag forces at the capillary channel entrance. The mathematical model incorporates pressure losses due to sudden contraction and viscous dissipation at the tube entrance, which can be significant in low Reynolds number flows. Additionally, it considers energy dissipation due to contact angle hysteresis. The model addresses an apparent anomaly posed by Washburn-Rideal and Levin-Szekely, and is applied to various liquids including water, glycerin, blood, oil, and methanol. It is tested with different geometries and cases, including numerical simulations, showing close agreement with experimental data. Deviations are observed when comparing infinite reservoir flow to finite droplet flow. A parametric study evaluates the effects of dimensionless numbers such as capillary, Reynolds, Weber, and Froude numbers. Results suggest the Weber number's importance over the Capillary number in droplet dynamics. The study also examines finite flow and film penetration in single pores versus pore networks. Computational simulations using ANSYS-FLUENT 23 R2 provide 2D results, using User Defined Functions (UDF) to capture liquid-gas interfaces. These simulations corroborate the mathematical model. Contrary to previous findings, this study demonstrates that contact angle effects are significant in the initial stages of capillary penetration. The proposed solution is valid for very short initial times, applicable to printing, lithographic operations, and filtration systems dealing with oil droplet removal from water using membranes. Copyright by MOHAMED ALHADDAD 2025 To my beloved parents, whose unwavering support and love have been my foundation. To my children, Hodaifa, Raneem, Razan, and Suzanne, who bring joy and inspiration to my life every day. And to my wife, whose companionship and encouragement have been my greatest strength. This work is a testament to your dedication and belief in me. Thank you for always being my pillars of strength. To my dear uncle, Abdullah. Though you are no longer with us, your memory continues to inspire and motivate me. Your wisdom, kindness, and the moments we shared will forever be cherished. This work is also dedicated to your enduring legacy. With all my love and gratitude, v ACKNOWLEDGMENTS I am deeply grateful for the support, friendly supervision and guidance of Professor Andre Benard. This milestone could not have been completed without his help, patience, and valuable advice. I also would like to thank my PhD committee for their valuable contribution and guidance. Thank you, Professor Petty, for the most valuable courses he taught me in multiphase flow phenomenon, and his constant guidance and motivation. Thank you, Professor Volodymyr V. Tarabara, for welcoming me within his research group for the discussion meetings that helped me to understand the filtration systems and membrane technologies. I also would like to thank Professor Abraham Engeda, for his guidance, contribution within my PhD committee, and motivation through my meeting with his PhD students. I also would like to thank Professor Wei Liao, for his guidance and motivation. vi TABLE OF CONTENTS CHAPTER 1: PREVIOUS RESEARCH AND PROBLEM STATEMENT………………….…1 CHAPTER 2: GOVERNING EQUATIONS AND MATHEMATICAL FORMULATION…...40 CHAPTER 3: MATHEMATICAL RESULTS OF DYNAMIC BEHAVIOR OF DROPLETS ENTERING CAPILLARY CHANNELS………………………………………..47 CHAPTER 4: GOVERNING EQUATIONS AND NUMERICAL SIMULATIONS FOR SIMULATING A DROPLET ENTERING A PORE……………………………68 CHAPTER 5: SIMULATION RESULTS AND DYNAMIC BEHAVIOR OF DROPLETS ENTERING CAPILLARY CHANNELS…………………….………………….72 CHAPTER 6: SUMMARY AND CONCLUSION……………………………………………...79 BIBLIOGRAPHY…….………………………………………………………………………….84 APPENDIX: USER DEFINED FUNCTION CODE………………………………….…..….....93 vii CHAPTER 1: PREVIOUS RESEARCH AND PROBLEM STATEMENT 1.1 Introduction The rise of liquids in capillary tubes has been investigated due to its importance for oil recovery, textiles, inkjet printing, fuel cells, and membrane separation. The ability of a fluid to rise in small channels or pores depends on the pressure applied at the entry, the droplet or film size, fluid properties, surface properties, and the pore shape. The classic equation of Washburn and Rideal[3] for the rate of penetration of fluid into capillaries due to surface tension forces is based on the Poiseuille flow across the capillary cross-section. However, such a model of fluid flow cannot apply at the entrance and exit of a capillary. Szekely[8] and Levine [11] have attempted a more rigorous theory of capillary penetration in which the entry flow comes from a large reservoir, and irreversible energy dissipation effects due to circulation during formation of a vena contracta at the entrance to the capillary were considered. They were able to remove the anomaly of the initial infinity in the fluid motion but their treatment of dissipative effects arising from circulation refers to high Reynolds number and flow from infinite reservoir. The flow in pores media or through network of pores have been studied using two-phase flow mechanics to study flow in the unsaturated, or vadose zone, of the subsurface, where empty spaces are partly filled by air and partly by liquid. Two-phase flow mechanics has also been used in filtration systems to characterize the flow of water and oil in oil-water separation design. In recent years, most modern applications and technologies of separation systems are applying two- flow mechanics to modeling the transport of contamination / fouling on membrane and in groundwater[15]. Two-phase flow concepts generally used to describe the interaction mechanisms between liquid-solid, liquid gas, and gas solid in pores media, such as contact angle, interfacial tension, and wettability. 1 In the next sections, important physical concepts needed are summarized to understand the two-phase flow phenomenon with emphasis on capillary flow. 1.2 Interfacial Tension and Capillary Pressure Interfacial energy arises at the contacting interface between two or more fluids, which can be a gas or another liquid barely miscible with the first liquid. This interfacial energy arises from the difference between the inward attraction of the molecules in the interior of each phase and those at the surface of contact [13]. This imbalance in forces exerts a tension on the interface of two immiscible liquids and causes contraction of the interface in small area as possible. This phenomenon is somewhat similar to the behavior of a stretched membrane under tension, and in contact with two fluids on either side of the membrane [15]. The interfacial tension of two immiscible liquids is defined as the amount of energy (N/m) that is required to create a unit area of interface between two immiscible liquids [16], and can be defined as the amount of work that must be performed in order to separate a unit of the fluid from the liquid [15]. In practice, the greater the interfacial tension between two immiscible liquids, the less likely an emulsion will be stable, so emulsions can efficiently be separated after mixing. In addition, interfacial tension decreases with increasing temperatures, surfactants, and gases in solution [17]. The values of interfacial tension of fluids that used in this study are water 72 mN/m, gasoline 21.6mN/m, glycerin 62mN/m, blood 55.89mN/m, alcohol 22mN/m. Capillary pressure arises because of interfacial tension between two contacted immiscible fluids. A discontinuity in pressure exists across the interface separating them. The magnitude of the pressure difference depends on the curvature of the interface considered. Capillary pressure is defined as the difference in pressure between a liquid and a fluid, where the pressures are taken in the two phases as the interface is approached from their respective sides. The Young–Laplace 2 equation describes the capillary pressure between two fluids and relates the capillary pressure to the surface of the wall. Figure 1.1 shows the forces acting at the curved liquid/fluid interface Figure 1.1 shows an infinitesimal element of a curved interface between two immiscible fluids labeled 1 and 2[13], where either both or only one of the two fluids is a liquid. Given the curvature, the pressure in fluid 1 is larger than that in fluid 2. The forces balance in components along the normal to the interface, with a constant surface tension , 𝜎(cid:2869)/(cid:2870) 𝑎𝑛𝑑 𝑟(cid:4593) = 𝑟(cid:4593)(cid:4593) = 𝑟, found that [15]: 𝑝(cid:3030) = 𝑝(cid:2869) − 𝑝(cid:2870) = 𝜎 (cid:3436) 1 𝑟(cid:4593) + 1 𝑟(cid:4593)(cid:4593)(cid:3440) = 2𝜎(cid:2869),(cid:2870) 𝑟 (1.1) where 𝑝(cid:2869) and 𝑝(cid:2870) denote to the pressures in fluid 1 and fluid 2, respectively, 𝑟(cid:4593) and 𝑟(cid:4593)(cid:4593) denote the principal radii of curvature, and 𝑟 is the mean radius of curvature. 3 1.3 Wettability and Contact Angle Wettability is the ability of a liquid to maintain contact with a solid. In a two-immiscible-fluid system, it is defined as the tendency of one fluid to spread preferentially over a solid surface in favor of the second fluid. This tendency is controlled by the balance of adhesive forces (liquid to surface) and cohesive forces (liquid to liquid)[18]. For instance, the behavior of a drop of water on an untreated glass surface can be observed: the droplet spreads on the surface and tends to stretch into a thin film. Under these conditions, the water droplet can be said to be wetting with respect to air. On the other hand, a droplet of mercury is observed to roll over the untreated glass surface with very little spreading. Therefore, in this system, mercury is said to be non-wetting with respect to air [15]. Wettability has been measured using three quantitative methods, contact angle, Amott method (imbibition and forced displacement), and US Bureau of Mines (USBM) [15][19]. In this research, at a specific surface, wettability is depicted by the concept of contact angle method, while Amott and USBM methods measure the average wettability of a core). Wettability is inversely proportional to the contact angle value; large wettability is happening at smaller contact angle [20]. The contact angle characterizes the wettability of liquid and substrate interface, and it is the best wettability measurement method when working on pure fluids and with absence of surfactants or other compounds altering the wettability[19]. In case of a droplet of test liquid formed on a solid surface with surrounding second fluid, as shown in Figure 1.2, the contact angle 𝜃(cid:3031) is often measured at the water phase [21],and it is defined as the angle of the interface plane between the solid surface plane and the tangent plane at the point of contact of the interface with the droplet surface. When 𝜃(cid:3031) = 0, the wetting fluid is said to be perfectly wet the solid surface, and when 𝜃(cid:3031) = 180 the surface is perfectly 4 hydrophobic or/and the liquid is perfectly non-wetting. In case of different fluid than water, the contact angle measured at the denser phase [14]. Figure 1.2 defines the wettability based on the contact angle and surface property In experimental studies, showed that if 𝜃(cid:3031) < 90, the fluid contained in the droplet is defined as the wetting fluid, while the surrounding fluid is non-wetting fluid [15]. In practice and for more specific ranges at 0 ≤ 𝜃(cid:3031) < 65 , the liquid is wetting, and at 105 ≤ 𝜃(cid:3031) < 180 the liquid is non-wetting. At 65 ≤ 𝜃(cid:3031) < 105 is a transient or intermediate case in which the surface has no specific definition regarding the preference for either fluid [21]. Capillary pressure can be defined as the difference between the pressure of the wetting fluid and the pressure of the non-wetting fluid, as illustrated in Figure 1.2, with fluid 1 defined as the non-wetting fluid. In practical applications, the concept of wettability has been applied to two- phase flow transport under specific conditions that require a comprehensive physical explanation. Thus, in the case of two-phase flow displacement in porous media, the fluid that penetrates a pore throat is said to be the wetting fluid, while the fluid repelled by capillary forces is said to be the non-wetting fluid. Therefore, capillary pressure is a measure of the tendency of a porous medium to imbibe the wetting phase or to repel the non-wetting phase [15]. 5 1.4 Viscosity and Density Density is the ratio of mass per unit volume of a substance, and it can be described in terms of Specific Gravity (Relative Density). SG is the ratio of the density of a solution to the density of a reference solution. The reference solution is usually water at 4°C[22]. Most studies show that there is no direct relation between viscosity, density, and surface tension. These three properties are independent of each other; however, viscosity, density, and surface tension are all affected by temperature [23]. In general, for any fluid, with an increase in temperature, the fluid’s density decreases, and thus its viscosity and surface tension decrease as well, making the fluid less viscous. In this study, five fluids have been examined as shown in Table 1.1[15]. Table 1.1 Five fluid properties that have been studied in this work[15] Density Viscosity (mPa Surface tension Fluid Water Gasoline Glycerin Blood (g/ml) 0.997 0.726 1.261 0.994 s) 0.894 0.6 950 3-4 Alcohol(ethanol) 0.789 1.095 (mN/m) 72 21.6 62 55.89 22 1.5 Kinetics of flow in capillary tubes; background and previous research work 1.5.1 Introduction The rise of liquids in capillary tubes has been investigated in several studies for many years due to its importance in applications such as; oil & gas recovery [24], textile[25], ink jet printing[26], fuel cell[27], agriculture[28], medicine[29], the pharmaceutical industry[30], nuclear engineering[31][32][33], fiber industry[34], ceramic industry[35][36][37], food 6 engineering[38], and environmental remediation[15]. In general, the rise of fluids in capillary channels is a topic of interest for any discipline dealing with porous media or materials. This research focuses on the dynamic aspects of capillary infiltration. The dynamics of liquid entering a network of pores, also called the kinetics of infiltration, aim to compute the rate of rise of a droplet of liquid into a capillary tube (pore) and thus gain knowledge of the liquid interface displacement as a function of time. Numerous studies and publications related to capillary penetration in channels can be found in the literature. Lucas, 1918 [1], a German scientist, presented the first paper published in a scientific journal that included a complete equation for the rate of infiltration in capillary tubes. However, because Lucas's paper was written in German, it is traditionally believed that Washburn, 1921 [2], was the first person to develop such a model. Therefore, the ordinary differential equation (ODE) that combines the rate of capillary rise with the total rise in a capillary tube often carries Washburn’s name, and in some publications, the ODE is referred to as the Lucas-Washburn equation [15]. The mathematical model describing capillary rise behavior can be traced back to the eighteenth century. Beneficiaries of this science have appreciated G. K. Mikhailov, who pointed to a letter from Georg Wolfgang Krafft to Euler (dated September 24, 1748), in which Krafft had put forward the puzzle: “The mechanical rule 𝑑𝑐 = 𝑝 𝑑𝑡/𝑚 or 𝑐 𝑑𝑐 = 𝑝 𝑑𝑠/𝑚 assumes that the mass 𝑚 is constant during the entire motion; would it not be possible to create such a rule in which the mass, or rather the moving point, could be variable?… I think that it would then become possible to derive the measure of the rise of fluids into capillary tubes, for in this case the mass that rises 7 increases all the time. Please, Sir, let me know your thoughts regarding this matter at some suitable opportunity” [39]. History has lost Euler’s reply to this letter, this could be a great example of the simulation of a case study associated with the dynamics of the system whose mass varies with time[15]. 1.5.2 Lucas-Washburn Equation Lucas & Washburn [1][2] presented the simplest ODE describing the infiltration of liquid into a capillary tube. Although other forms of the equation were obtained by Bell and Cameron (1906, published in The Journal of Physical Chemistry, 2002)[40] and West, 1911 [41], “they did not start using the Hagen–Poiseuille equation, which was derived independently by Jean Léonard Marie Poiseuille in 1838 and Gotthilf Heinrich Ludwig Hagen, and published by Poiseuille in 1840–41 and 1846 [42]. The theoretical justification of the Poiseuille law was given by George Stokes in 1845”[43] &Wikipedia. Lucas and Washburn both had started from the Hagen-Poiseuille’s law, which relates the steady-state flow of a liquid through a capillary tube [44] 𝑄 = (cid:2872) 𝜋𝑟(cid:2868) 8𝜇𝑙 ∆𝑝 (1.2) Where 𝑄 is the volumetric flow of liquid, 𝑟(cid:2868) and 𝑙 are the radius and the length of the capillary tube, respectively, 𝜇 is the dynamic viscosity of the liquid, and ∆𝑝 is the total driving pressure acting to force the liquid along the capillary. Based on the Equation (1.2), Lucas and Washburn obtained [15] 𝑑𝑧(cid:3030) 𝑑𝑡 = (cid:2870) 𝑟(cid:2868) 8𝜇 ∆𝑝(𝑧(cid:3030)) 𝑧(cid:3030) (1.3) Where 𝑧(cid:3030) is the penetrated length of liquid displacing air inside a capillary tube at time 𝑡, as illustrated in Figure 1.3. Assuming the parabolic profile of the velocity distribution of Hagen- 8 Poiseuille’s flow, 𝑑𝑧(cid:3030) 𝑑𝑡⁄ is a mean velocity across a section of capillary tube, and displaced air (gas phase) velocity has been assumed to be negligible with respect to the viscosity of the liquid[15]. Figure 1.3 shows the schematic of infiltration into a capillary tube with 𝑟(cid:2868) radius and inclination angle 𝜑 to the horizontal With the assumption of the liquid is wetting with respect to the air, the total pressures are acting on the system as show in Figure 1.3 ∆𝑝(𝑧(cid:3030)) = 𝜌𝑔ℎ − 𝜌𝑔𝑧(cid:3030) sin 𝜑 + 2𝜎 cos 𝜃 𝑟(cid:2868) (1.4) where 𝑔 is the gravity, 𝜌 density, surface tension 𝜎 , contact angle 𝜃 , ℎ reservoir head (assumed constant with time), and 𝜑 is the inclination angle of the capillary tube to the horizontal, as illustrated in Figure 1.3. Due to reservoir and the capillary tube are open to the air, the atmospheric pressure was ignored in Equation (1.4). As a result, based on experimental measurements, with 𝑟(cid:2868) known, Equation (1.4) can be used to study the capillary infiltration into porous solids such as powder beds in order to characterize their wettability [2][45][10][46][47]. 9 In the case of 𝜑 = 0 where the capillary tube is horizontal, and ∆𝑝(𝑧(cid:3030)) is independent of 𝑧(cid:3030), Washburn had solved Equation (1.4) for these conditions and obtained 𝑧(cid:3030) = (cid:3496) (cid:2870) 𝑟(cid:2868) 4𝜇 (cid:3436)𝜌𝑔ℎ + 2𝜎 cos 𝜃 𝑟(cid:2868) (cid:3440) (1.5) which the above equation shows that the penetration length of a horizontal capillary tube is proportional to the square root of time. This relation was theoretically approved and experimentally verified with water, alcohol, and benzene for 0.5 mm radius capillary tubes by Bell and Cameron [40]. However, Levy noted from Equation (1.5) that the flow is not driven by the interfacial tension difference. Instead, it results from the laminar flow of a Newtonian liquid in a tube under the action of a constant pressure difference, where the penetration length of the liquid changes with time [15]. On the other hand, in case of 𝜑 = 90(cid:2868), the capillary tube is vertical and located above the reservoir, so ∆𝑝(𝑧(cid:3030)) is dependent of 𝑧(cid:3030). Under these conditions, Lucas and Washburn solved Equation (1.4) and obtained 𝑧(cid:3030) + (ℎ + ℎ(cid:3032)) ln (cid:3436)1 − 𝑧(cid:3030) ℎ + ℎ(cid:3032) (cid:3440) = − (cid:2870)𝜌𝑔 𝑟(cid:2868) 8𝜇 𝑡 where ℎ(cid:3032) is the well-known equilibrium height at 𝑡 → ∞ given by[15] ℎ(cid:3032) = 2𝜎 cos 𝜃 𝜌𝑔𝑟(cid:2868) (1.6) (1.7) In the case of hydrostatic forces are negligible with respect to capillary forces, so the Equations (1.6) & (1.7) both reduce to 𝑧(cid:3030) = (cid:3496) 𝑟(cid:2868)𝜎 cos 𝜃 2𝜇 (1.8) 10 Infiltration into capillary tube has been described by the Equations (1.6) & (1.8). The credit of initial study was given to Washburn, 1921, who has obtained excellent agreement, the following Table 1.2 show the validation of the penetration behavior proportional parabolic with time validated by different authors. The authors found that the flow rate is smaller than expected with presence of bubbles within the liquid. In all these studies has been done with the assumption of initially pre-wetted capillary tube. The studies have been shown the penetration into an initially dry tube could only be predicted correctly in case surface tension property of water equal to 0.0385 𝑁/𝑚 [15]. Table 1.2 Validation of the penetration behavior proportional parabolic with time validated by different authors Author Year Fluid Capillary length Radius (m) (mm) Washburn[2] 1921 0.95 m 0.15 , 0.37 water Pre-wetted Rideal[3] 1922 1.2 0.35 solvents Pre-wetted Malik et al[110] 1979 0.2 , 0.3 Water & alcohol Pre-wetted Fisher and 1979 Very small 0.3 0.15-0.2 water and Pre-wetted Lark[111] µm µm cyclohexane Peek and 1934 0.25-0.36 water Pre-wetted McLean[45] Ligenza & 1951 20-40 µm solvents Pre-wetted Bernstein[49] Studies have shown the successful application of Washburn equation. Washburn obtained reasonable agreement both theoretically and experimentally for capillary rise into a 0.15 mm radius capillary tube using Equation (1.8). Ligenza and Bernstein [49] illustrated very good 11 agreement between Equation (2.6) and experimental data on rise of various solvents into 20-40 µm radius fine vertical capillary tube. Despite the reported successes of most studies validating the Lucas-Washburn equation, high- resolution observations in capillary tubes often conclude that the Lucas-Washburn equation fails to describe the initial stage of wetting liquid infiltration (rapid change case). Initial penetration of the fluid has been observed to be linear with time [50]. At the initial moment of penetration, the penetration rate is much larger than those measured experimentally [15][48][51][52][53]. Siebold (2000) showed the rise of pentane into a 0.191 mm radius vertical capillary tube and compared it with the predictions using Equations (1.6)at ℎ = 0 & (1.8) at 𝜃 = 0 as illustrated in Figure 1.4[48] Figure 1.4 illustrates the rise of pentane into a 0.191 mm radius vertical capillary tube: (1) experimental points; (2) Eq. (1.6) taking into account the hydrostatic pressure; (3) Eq. (1.6) neglecting the hydrostatic pressure[48].Used with permission of Elsevier Science & Technology Journals, from Siebold, Alain, Michel Nardin, Jacques Schultz, André Walliser, and Max Oppliger. "Effect of dynamic contact angle on capillary rise phenomena." Colloids and surfaces A: Physicochemical and engineering aspects 161, no. 1 (2000): 81-87; permission conveyed through Copyright Clearance Center, Inc 12 The graph in Figure 1.4 illustrates the rise of pentane in a vertical capillary tube with a radius of 0.191 mm. It compares experimental data with theoretical predictions from Equation (1.6), both considering and neglecting hydrostatic pressure. The experimental points closely follow the dashed line, indicating that hydrostatic pressure plays a significant role in the capillary rise of pentane. The solid line, which neglects hydrostatic pressure, deviates from the experimental data, showing that ignoring this factor leads to less accurate predictions. The study highlights the importance of considering hydrostatic pressure in capillary rise phenomena, especially when dealing with small capillaries and dynamic contact angles [48]. Using Lucas-Washburn’s equation (1.6) gives infinite initial penetration rate (𝑑𝑧(cid:3030) 𝑑𝑡)⁄ at initial condition 𝑧(cid:3030)(𝑡 = 0) = 0. Washburn was aware of the fact of infinite initial penetration rate, but he stated that the “Lucas-Washburn’s equation (2.3)” is valid whenever Hagen- Poiseuille’s law is applicable at 𝑅𝑒 < 2000. This analysis fails to consider the significant effects of inertia—the forces associated with the change in momentum of the liquid. Inertia forces are significant in the laminar region at Reynolds 𝑅𝑒 < 2000. The length scale of the region in which the Lucas-Washburn’s equation (1.6) is not valid includes his implicit assumption of quasi-static flow[15]. Therefore, Lucas-Washburn’s equation (1.6) is inappropriate to describe the initial conditions where rapid changes in penetration occur, such as the initial conditions of the meniscus rise into a capillary tube. 1.6 Penetration Model with consideration of Inertia Forces 1.6.1 Equation of motion using momentum balance Attempts have been made to address the failure of the Lucas-Washburn equation at the initial stage of wetting liquid infiltration. Scholars have used the momentum conservation equation to include inertia forces and derive kinetic models for liquid penetration into capillary 13 tubes[3][4][5][6][11][50]. Bosanquet[4], Quéré[50] and Pickett[5] considered the momentum difference at an infinitesimal change of time 𝑡 and penetration length 𝑧(cid:3030), in the initial stage of an air-filled capillary tube with radius 𝑟(cid:2868). They obtained the following expression 𝜌 (cid:3436) (cid:2870) (cid:3440) 𝑑𝑧(cid:3030) 𝑑𝑡 + 𝜌𝑧(cid:3030) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) = ∆𝑝(𝑧(cid:3030)) − 8𝜇 (cid:2870) 𝑧(cid:3030) 𝑟(cid:2868) 𝑑𝑧(cid:3030) 𝑑𝑡 Rideal, as well as Levine and Neale, used a similar approach but obtained [15] 𝜌𝑧(cid:3030) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) = ∆𝑝(𝑧(cid:3030)) − 8𝜇 (cid:2870) 𝑧(cid:3030) 𝑟(cid:2868) 𝑑𝑧(cid:3030) 𝑑𝑡 (1.9) (1.10) The Bosanquet model does not consider the momentum of the element of length ∆𝑧(cid:3030) at the reservoir feeding the capillary tube at time t. Their derivation implicitly assumed ∆𝑧(cid:3030) has no momentum[15]. Levine et al. (1976) studied in depth the dynamic behavior of the element entering the capillary tube, which will be discussed in more details in next chapter. Equation (1.10) solved with the assumption of the total driving pressure ∆𝑝(𝑧(cid:3030)) is independent of 𝑧(cid:3030), 𝑧(cid:3030)(𝑑𝑧(cid:3030) 𝑑𝑡)(cid:3047)(cid:2880)(cid:2868) = 0 ⁄ , and horizontal capillary tube 𝑧(cid:3030) = (cid:3496) (cid:2870)∆𝑝 𝑟(cid:2868) 4𝜇 (cid:3428)𝑡 − 𝜏 (cid:3436)1 − 𝑒(cid:2879) (cid:3047) (cid:3099)(cid:3440)(cid:3432) where 𝜏 is a time scale at inertia forces are significant and is given by 𝜏 = (cid:2870) 𝜌𝑟(cid:2868) 8𝜇 (1.11) (1.12) In case of 𝑡 ≫ 𝜏, (1.10) is reduced to (1.5). For the horizontal capillary tube and reservoir head ℎ = 0, so the time scale 𝜏 corresponds to a length scale 𝜁can be given by (1.12) 𝜁 = (cid:3496) 𝜌𝑟(cid:2868) (cid:2871)𝜎 cos 𝜃 16𝜇(cid:2870) (1.13) 14 The effects of inertia can be examined by Equations (1.12) & (1.13). For instance, consider the penetration of water in a 1 mm radius horizontal pre-wetted capillary tube. Using the physical properties of water at 20(cid:2868)𝐶 as shown in Table 1.3[44], The significant effects of inertia can be noted over the length of 10 to100 times the capillary tube radius. In case of vertical capillary tube taking 𝑔 = 9.807 𝑚/𝑠(cid:2870), Equation (1.6) can be solved numerically to get ℎ(cid:3032) which is same as a length scale 𝜁, and compare the result to the Equation(1.7) Table 1.3 Effects of inertia using Equations (1.12) & (1.13) Radius Viscosity mm Pa.s Surface tension N/m Density Kg/m3 𝝉 𝒔 𝜻 mm horizontal 1 1.002E-3 7.28E-2 998.2 0.125 67.3 horizontal 0.1 1.002E-3 7.28E-2 998.2 1.25E-3 2.13 Case 1 2 Table1.4 show that for a 𝑟(cid:2868) = 1 𝑚𝑚 capillary tube radius, 𝜁 = 14.9 𝑚𝑚 equal to ℎ(cid:3032) using Equation(1.7), while for a 𝑟(cid:2868) = 0.1 𝑚𝑚 capillary tube radius, 𝜁 = 2.12 𝑚𝑚 which is not equal to ℎ(cid:3032) = 148.7 using Equation(1.7)[44]. Therefore, the effects of inertia are significant for 𝑟(cid:2868) = 1 𝑚𝑚, while for 𝑟(cid:2868) = 0.1 𝑚𝑚, they are only important with the region less than 1.5% of ℎ(cid:3032). Table 1.4 Effects of inertia using Equations (1.12) & (1.13) Case Radius Viscosity mm Pa.s Surface tension N/m Density 𝝉 𝜻 𝒉𝒆 Kg/m3 𝒔 mm mm 3 Vertical 1 1.002E-3 7.28E-2 4 Vertical 0.1 1.002E-3 7.28E-2 998.2 998.2 14.9 14.9 2.12 148.7 15 At the early stage of capillary penetration, if viscosity and hydrostatic forces are neglected, Equation (1.9) reduced to 𝜌 𝑑 𝑑𝑡 (cid:3436)𝑧(cid:3030) 𝑑𝑧(cid:3030) 𝑑𝑡 (cid:3440) = 2𝜎 cos 𝜃 𝑟(cid:2868) (1.14) Equation.(1.14) can be solved at initial condition 𝑧(cid:3030) = 0 and finite penetration rate 𝑑𝑧(cid:3030) 𝑑𝑡⁄ , obtained [15] where c is the penetration rate given by 𝑧(cid:3030) = 𝑐 𝑡 𝑐 = (cid:3496) 2𝜎 cos 𝜃 𝜌𝑟(cid:2868) (1.15) (1.16) Penetration rate 𝑐 is sometimes referred to as the Bosanquet velocity [15] To solve Equation (1.9), it might seem natural to take the initial conditions 𝑧(cid:3030)((cid:3047)(cid:2880)(cid:2868)) = 0, and (𝑑𝑧(cid:3030) 𝑑𝑡)⁄ (cid:3047)(cid:2880)(cid:2868) = 0. However, taking these initial conditions leads to a singularity at 𝑡 = 0 when a finite force is applied to an infinitesimal mass. To eliminate the singularity, a finite initial velocity is assumed, as given by Equation (1.16)[15]. A drawback is an infinite acceleration at 𝑡 = 0 [15][54], which was solved by Szekely[7] and will be discussed in more details in the coming sections. Reynolds number computed corresponding to the Bosanquet velocity as[15] 𝑅𝑒(𝑐) = (cid:3496) 8𝑟(cid:2868)𝜌𝜎 cos 𝜃 𝜇(cid:2870) (1.17) and found that the flow associated with the capillary infiltration is laminar. Experimental observations of capillary penetration have shown initial behavior as predicted by Equation (1.15), where the penetration length increases linearly with time at the very early 16 stages of liquid infiltration [48][50]. However, the observed penetration rate was found to be less than the predicted velocity 𝑐 [15][50][53][55]. Quéré[50] identified the inconsistency between the observed and predicted velocities as due to[15]:  A difference between the dynamic contact angle and the static contact angle. This difference arises from the characteristics of the flow at the contact line between the liquid/air interface and the capillary wall.  Partial energy losses at the boundary between the capillary tube and its reservoir. These losses are caused by the pressure distribution resulting from the flow from the reservoir into the capillary. Levine and Neale[11] have pointed out to another issue in the both Equations (1.9) and (1.10). Their left-hand side terms indicate that the liquid moves through the capillary tube as a solid plug, suggesting a constant velocity distribution through the capillary tube cross-section[15]. The contradict to the suggested constant velocity profile can be found in left-hand side terms that contains a term corresponding to the Hagen-Poiseuille’s viscous dragging force, which is based on a steady flow parabolic velocity distribution across a section of the capillary tube[15]. Levine and Neale[11], Quéré[50], and others have disregarded another issue: an implicit assumption in their equations states that streamlines are parallel to the walls of the capillary tube and the main direction of flow. Most scientists and researchers have written force balances for the column of liquid in motion instead of using the Navier-Stokes equations. For example, Munson et al.[44] have failed to incorporate this term into account[3][4][5][7][8][50][54]. This issue is a core part of this research and will be explored further in coming sections. 17 1.6.2 Equation of motion using Navier-Stokes Equation Letelier et al.[56] has mad attempt to remove the limitations associated with the Hagen- Poiseuille’s friction approximation in unsteady flow. In their derivations, they kept the assumption of streamline enters parallel to the longitudinal direction of capillary tube. In other words, convective inertia and radial velocity assumed to be zero, and longitudinal velocity profile was not assumed to be parabolic. Letelier et al. used Navier-Stokes Equation and wrote series expansion of velocity distribution and obtained 3 4 𝜌𝑧(cid:3030) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) = ∆𝑝(𝑧(cid:3030)) − 8𝜇 (cid:2870) 𝑧(cid:3030) 𝑟(cid:2868) 𝑑𝑧(cid:3030) 𝑑𝑡 + 1 144 (cid:2870) 𝜌(cid:2870)𝑟(cid:2868) 𝜇 𝑧(cid:3030) 𝑑(cid:2871)𝑧(cid:3030) 𝑑𝑡(cid:2871) + 𝑜 (cid:4678) 𝑑(cid:2872)𝑧(cid:3030) 𝑑𝑡(cid:2872) (cid:4679) (1.18) where 𝑜(𝑑(cid:2872)𝑧(cid:3030) 𝑑𝑡(cid:2872)) ⁄ is the rest of series of order derivatives of z_c in the fourth power and above. Lucas-Washburn equation is obtained from Equation (1.18) by eliminate the order ≥ 2, and Rideal equation is obtained by eliminate the order ≥ 3. The importance of higher order is associated with smaller time ranges. By inspection Equation (1.18) illustrates the term order ≥ 2 can be negligible with respect to the first-order term when the time 𝑡 ≫ 𝜏(cid:2869), where 𝜏(cid:2869) = 3 4 (cid:2870) 𝜌𝑟(cid:2868) 8𝜇 (1.19) Which like Equation (1.12). Similarly, the term order ≥ 3 can be negligible with respect to the second-order term when the time 𝑡 ≫ 𝜏(cid:2870), where 𝜏(cid:2869) = 1 144 3 4 (cid:2870) 𝜌𝑟(cid:2868) 8𝜇 (1.20) Letelier et al.[56] have not validated the derivation experimentally, but their model was evaluated by Batten[15][57]. Batten[57] used data obtained by LeGrand and Rense[51]who investigated the rise of water and ethanol in capillary tubes of radius ranging from 0.242 mm to 0.350 mm [15]. Batten did not claim that Equation(1.18)could predict the capillary rise more 18 accurately than was found by the Bosanquet equation (1.9). However, it should be stated out that Batten neglected the terms ≥ 2. 1.7 Penetration Model with consideration of Inertia Forces Between Capillary tube and Reservoir 1.7.1 Hagenbach Correction The first consideration of the effects of drag forces was attempted by Brittin [60], who added a term to the right side of Bosanquet’s equation (1.9)[4]. This term relates to the friction loss due to cross-sectional contraction at the inlet of the capillary tube. This loss can be attributed to eddies forming at the side of the tube entrance[15][58]. By neglecting the cross-sectional area of the capillary tube compared to the cross-sectional area of the reservoir, the contraction drag force 𝑓(cid:3030)(cid:3031) can be found by 𝑓(cid:3030)(cid:3031) = 𝑒(cid:3092) 2 𝜌𝜋𝑟(cid:2868) (cid:2870) (cid:3436) (cid:2870) (cid:3440) 𝑑𝑧(cid:3030) 𝑑𝑡 (1.21) where 𝑒(cid:3092) is a dimensionless coefficient equal to 0.5[15]. The Vena contracta drag force has been examined by Szekely using the same approaches mentioned above, and he determined a coefficient 𝑒(cid:3092) = 0.45 [15][59]. In general, 𝑒(cid:3092) is determined by the type of connection between the capillary tube and the reservoir. This value will be lower in the case of a well-rounded or “trumpet-shaped” connection type [15]. Conversely, the vena contracta drag force coefficient is larger if the capillary tube is reentrant[44]. In general case, Bosanquet’s equation (1.9) can be rewritten as (cid:2870) (cid:3440) 𝜌 (cid:3436) 𝑑𝑧(cid:3030) 𝑑𝑡 + 𝜌𝑧(cid:3030) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) = ∆𝑝(𝑧(cid:3030)) − 8𝜇 (cid:2870) 𝑧(cid:3030) 𝑟(cid:2868) 𝑑𝑧(cid:3030) 𝑑𝑡 − 𝑒(cid:3092) 2 𝜌 (cid:3436) (cid:2870) (cid:3440) 𝑑𝑧(cid:3030) 𝑑𝑡 (1.22) Brittin [60] used Equation (1.22) and obtained good agreement with one set of data presented by Rense[6][15]. Levine et al.[11] criticized the approach used by Brittin [60]and Szekely, 19 pointing out that the drag force given by Equation (1.21) was only valid at high Reynolds numbers for turbulent flow. In this research, the approach derived by Levine et al.[11] was used, which can be valid for laminar flow in porous media at very low Reynolds. In the case of laminar flow into a capillary tube, there are two types of pressure drops to consider: one in the inlet region and the other corresponding to the conversion of pressure energy to kinetic energy at the inlet of the capillary tube. The first approximation of the inlet pressure can be found by assuming inviscid flow in the reservoir and a constant penetration rate 𝑑𝑧(cid:3030) 𝑑𝑡⁄ . Thus, the inlet pressure 𝑝(cid:2869)(𝑧(cid:3030)) has been found using Bernoulli’s equation as [61][15] 𝑝(cid:2869)(𝑧(cid:3030)) = 𝜌𝑔ℎ − 1 2 𝜌 (cid:3436) 𝑑𝑧(cid:3030) 𝑑𝑡 (cid:2870) (cid:3440) (1.23) With viscous flow conditions within the entry development region, a constant velocity distribution over the cross-section of the circular tube is considered. The boundary layer grows to ensure mass conservation until a parabolic velocity distribution is achieved [15]. Bird et al.[59] defined the length of the developing region in the order of 0.07 𝑟(cid:2868)𝑅𝑒. White [62] suggest 𝑙(cid:3032) = 0.16𝑟(cid:2868)𝑅𝑒 + 1.3𝑟(cid:2868), where the second term 1.3𝑟(cid:2868) is independent of the Reynolds number, making the formula valid for creeping flow at very low Reynolds numbers ≪ 1[63][64]. Studies have shown that the magnitude of kinetic energy associated with a constant velocity distribution is half that of the kinetic energy associated with a parabolic velocity distribution. Additionally, there is a pressure loss due to the transition from constant to parabolic velocity distribution[15] [65]. The drag force corresponding to the pressure loss because of transitioning velocities, can be expressed in the form similar to Equation (1.21). By computing the energy differences using constant and parabolic velocity distributions, the vena contracta drag force coefficient 𝑒(cid:3092) = 1 [15][65]. By combining Equation (1.21) and the second term of left hand side of Equation (1.21) to form a pressure drop 20 ∆𝑝(𝑧(cid:3030)) = (1 + 𝑒(cid:3092)) (cid:2870) (cid:3440) 𝜌 2 (cid:3436) 𝑑𝑧(cid:3030) 𝑑𝑡 (1.24) the coefficient (1 + 𝑒(cid:3092)) is replaced by 𝑚 which is called Hagenbach correction [66][15]. In case of 𝑚 = 2, reduced Equation (1.21) to 𝜌(𝑑𝑧(cid:3030) 𝑑𝑡⁄ )(cid:2870) which is the same term that differentiates Rideal Equation (1.21) from Bosanquet Equation (1.21). This concludes that pressure loss term has already included in Bosanquet[4] analysis but does not in Rideal[3]. Table 1.5 Hagenbach correction from different authors[15] Author Average experimental 𝐦 𝐦 Langhaar Boussinesq Schiller 2.28 2.24 2.16 Atkinson-Goldstein 2.41 Riemann Schiller Hagen 2.248 2.32 2.7 More studies have been done to estimate Hagenbach correction 𝑚 using different approaches as illustrated in Table 1.5. Most of these theories assume that the flow from the reservoir to the capillary tube takes place through a well-rounded entrance, exclude Hagen, the tube entrance was not well rounded but squarely cut off[15]. Levy[15] has summarized three combined effects that attribute to the pressure drop across the contraction from the reservoir to the capillary tube as following:  Conversion of pressure energy to kinetic energy in the reservoir near the inlet of the capillary tube 21  Conversion of pressure energy to kinetic energy due to a change in velocity distributions  Eddy dissipation at the contraction Taking the pressure drop terms given by Equations (1.21) & (1.24), the using Equation (1.4) and Equation (1.10)[15] 𝜌𝑧(cid:3030) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) + 8𝜇 (cid:2870) 𝑧(cid:3030) 𝑟(cid:2868) 𝑑𝑧(cid:3030) 𝑑𝑡 = 𝜌𝑔ℎ + 2𝜎 cos 𝜃 𝑟(cid:2868) − 𝜌𝑔𝑧(cid:3030) sin 𝜑 − 𝑚 2 𝜌 (cid:3436) (cid:2870) (cid:3440) 𝑑𝑧(cid:3030) 𝑑𝑡 (1.25) The above equation was first proposed by Siegel [7]. By taking 𝑚 = 2 in Equation (1.25), the resulting equation is not substantially different from Equation 2.9. Siegel had taken 𝑚 = 5 before he observed that existence of some agreement between prediction and measurement on a micro-gravity rise in a capillary tube with radius 0.95 mm[15]. In case of the flow can be turbulent before entering the capillary tube, Zhmud et al[54] has suggested the possibility of existing a turbulence drag could slow down the penetration rate, and added a term to the right-hand side of Equation (1.25) equal to 0 𝑖𝑓 𝑑𝑧(cid:3030) 𝑑𝑡⁄ ≤ 𝑣(cid:3030)(cid:3045) 𝜙(𝑧(cid:3030)) = (cid:4688) − (cid:2870) 𝑞 (cid:2870) 𝑧(cid:3030) (cid:3436) 𝑟(cid:2868) 𝑑𝑧(cid:3030) 𝑑𝑡 (cid:3440) 𝑖𝑓 𝑑𝑧(cid:3030) 𝑑𝑡⁄ > 𝑣(cid:3030)(cid:3045) (1.26) where 𝑣(cid:3030)(cid:3045) is the critical velocity in which the turbulence begins, 𝑞 is a turbulence coefficient taken to be equal to 0.3 𝑘𝑔 𝑚(cid:2871)⁄ . A very good fit has been obtained by Zhmud et al[54]for the capillary rise of dodecane in a 0.1 mm radius of the capillary tube. However, they noted that the critical velocity 𝑣(cid:3030)(cid:3045) corresponded to a Reynolds number of the order of 2. Therefore, according to empirical approach used by the authors, 𝜙(𝑧(cid:3030)) was not exactly a turbulent drag, but rather it is the second order dissipation correction related to the actual flow pattern[15]. 22 1.7.2 Reservoir Inertia The inlet pressure 𝑝(cid:2869)(𝑧(cid:3030)) has been found using Bernoulli’s equation as shown in Equation (1.24) with assuming constant velocity in the capillary tube, which of course is not true case. These changes in velocities are responsible for reservoir inertia and must be examined. Siegel[7] adopted the study of flow through an orifice done by Morse and Feshbach [15][67]. For the flow rate of 𝜋𝑟(cid:2868) (cid:2870) (𝑑𝑧(cid:3030) 𝑑𝑡⁄ ) through an orifice, the authors illustrated that a plug of effective mass 𝜌𝜆𝜋𝑟(cid:2868) (cid:2871) ; 𝜆 = 𝜋 2⁄ , had to be accelerated when initiating flow. For the inertia term of Equation (1.25), Siegel 1961, found out this inertia was equivalent to increasing in liquid height in the capillary from 𝑧(cid:3030) to 𝑧(cid:3030) + 𝜆𝑟(cid:2868), then he obtained 𝜌(𝑧(cid:3030) + 𝜆𝑟(cid:2868)) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) + 8𝜇 (cid:2870) 𝑧(cid:3030) 𝑟(cid:2868) 𝑑𝑧(cid:3030) 𝑑𝑡 = 𝜌𝑔ℎ + 2𝜎 cos 𝜃 𝑟(cid:2868) − 𝜌𝑔𝑧(cid:3030) sin 𝜑 − 𝑚 2 𝜌 (cid:3436) 𝑑𝑧(cid:3030) 𝑑𝑡 (cid:2870) (cid:3440) (1.27) The empirical choice of 𝜆 = 𝜋 4⁄ leads to a good agreement with Siegel’s experimental data. In addition, initial condition of zero velocity can be compatible with Equation (1.27) because the correction term has removed the singularity associated with Equation (1.9) at 𝑡 = 0. Although the Siegel’s approach of reservoir inertia analysis seems to be empirical, it was independently validated by Szekely by using an innovative method of examining a case of liquid rise into a vertical capillary tube. Instead of using the momentum conservation equation, they analyzed a system of liquid rise into capillary tube by using macroscopic energy balance as[15] 𝑑 𝑑𝑡 (𝐾𝐸 + 𝑃𝐸) = −∆ (cid:4680)(cid:4678) 𝑢(cid:2870) 2 + 𝑔𝑧 + 𝑝 𝜌 (cid:4679) 𝑤(cid:4681) − 𝑊 − 𝐸(cid:3092) (1.28) where 𝐾𝐸 & 𝑃𝐸 are total kinetic energy and potential energy within the system, respectively, the longitudinal liquid velocity 𝑢 = 𝑑𝑧(cid:3030) 𝑑𝑡⁄ assumed to be constant throughout the system, at system pressure 𝑝, mass flow rate 𝑤 = 𝜌𝑢𝐴; 𝐴 cross section area, −𝑊 a rate of work done by the system, and the rate of work dissipated irreversibly 𝐸(cid:3092). The three terms inside the parentheses are 23 respectively, (net input of kinetic energy + the net input of potential energy + the net input of pressure energy) at the cross-sectional boundaries of the system. This case of input energy comes through the boundary with the reservoir and no output energy since the upper boundary is rising along with the liquid[15]. The terms 𝑊 and 𝐸(cid:3092) are respectively, correspond to the rate of works due to viscosity forces and contraction drag force. After the calculations, Szekely obtained 𝜌𝑧(cid:3030) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) = Δ𝑝(𝑧(cid:3030)) − 8𝜇 (cid:2870) 𝑧(cid:3030) 𝑟(cid:2868) 𝑑𝑧(cid:3030) 𝑑𝑡 − 𝑒(cid:3092) 2 𝜌 (cid:3436) 𝑑𝑧(cid:3030) 𝑑𝑡 (cid:2870) (cid:3440) (1.29) It can be noticed that Equation (1.29) includes the form of Rideal’s equation (1.10)[3], in additional to energy loss term similar to that has proposed by Brittin[60]. It should be pointed out that Szekely’s analysis is consistent with Rideal’s analysis with the assumption of input kinetic energy to be equal to 0.5𝜌𝜋𝑟(cid:2868) (cid:2870)(𝑑𝑧(cid:3030) 𝑑𝑡⁄ )(cid:2871)for elements of liquid entering the capillary tube[15]. Szekely has proposed different techniques to treat the driving pressure given by Equation (1.4), they included the pressure distribution due to the flow from reservoir into capillary tube, which they obtained ∆𝑝(𝑧(cid:3030)) = −𝜌𝑔𝑧(cid:3030) + 2𝜎 cos 𝜃 𝑟(cid:2868) + 𝑝(cid:2869)(𝑧(cid:3030)) (1.30) where 𝑝(cid:2869)(𝑧(cid:3030)) cross section inlet pressure. For a vertical case at 𝜑 = +90(cid:2868), 𝑝(cid:2869)(𝑧(cid:3030)) that given by Equation(1.23) and Equation (1.29) is identical to Equation (1.29). 24 Figure 1.5 shows the far and near field of the capillary tube reservoir Szekely has assumed that the connection level of the capillary tube to the reservoir at reservoir head ℎ = 0, they applied an energy balance Equation (1.28) and expression for inlet pressure 𝑝(cid:2869)(𝑧(cid:3030)) has been obtained. The element of liquid flow from reservoir can be seen as a hemispherical body extending through the inlet of the capillary to infinity as shown in Figure 1.5. The centripetal velocity 𝑢(cid:3045), in the far field for any hemispherical liquid of radius 𝑟 ≥ 𝑟(cid:2868), can be obtained from mass conservation equation in spherical coordinates[68] as 𝑢(cid:3045) = 1 2 (cid:4672) 𝑟(cid:2868) 𝑟 (cid:4673) (cid:2870) 𝑑𝑧(cid:3030) 𝑑𝑡 (1.31) Therefore, the corresponding kinetic energy can be measured, and it has been covered through the study of flow towards a spherical sink or away from a point source, for both viscous and inviscid flow[15]. In the near field, within the element of hemispherical liquid of radius 𝑟 ≤ 𝑟(cid:2868), as shown in Figure 1.5, the velocity field is not known. For this region, Szekely assumed that, the velocity within the near region, was equal to the main capillary tube velocity 𝑑𝑧(cid:3030) 𝑑𝑡⁄ , which is the actual 25 velocity of the upper-bond. They obtained an expression for 𝑝(cid:2869)(𝑧(cid:3030)) by deriving the reservoir total kinetic energy 𝑝(cid:2869)(𝑧(cid:3030)) = − 7 6 𝜌𝑟(cid:2868) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) (1.32) However, Szekely presented innovative technique in their energy balance analysis, but it seems to be mistakenly neglecting the outflow of kinetic energy, 0.5𝜌𝜋𝑟(cid:2868) (cid:2870)(𝑑𝑧(cid:3030) 𝑑𝑡⁄ )(cid:2871), from capillary tube to reservoir. By including outflow KE, they should have obtained 𝑝(cid:2869)(𝑧(cid:3030)) = − 7 6 𝜌𝑟(cid:2868) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) − 1 2 𝜌 (cid:3436) 𝑑𝑧(cid:3030) 𝑑𝑡 (cid:3440) (1.33) which is identical to Equation (1.23), but with additional tear of reservoir inertia. Substituting Equation (1.33) into Equation (1.30) and Equation (1.29) gets 𝜌 (cid:3436)𝑧(cid:3030) + 7 6 𝑟(cid:2868)(cid:3440) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) + 8𝜇 (cid:2870) 𝑧(cid:3030) 𝑟(cid:2868) 𝑑𝑧(cid:3030) 𝑑𝑡 = 2𝜎 cos 𝜃 𝑟(cid:2868) − 𝜌𝑔𝑧(cid:3030) − 𝑚 2 𝜌 (cid:3436) 𝑑𝑧(cid:3030) 𝑑𝑡 (cid:2870) (cid:3440) (1.34) which analogous to Siegel’s equation[7] at 𝜑 = +90(cid:2868), ℎ = 0, and 𝜆 = 7 6⁄ . Similar approach was used by Huang et al.[69], but 𝜆 = 3 8⁄ . Based on using incomplete Equation (1.32), Szekely should have obtained the coefficient 𝑒(cid:3092) instead of 𝑚 = 1 + 𝑒(cid:3092). In fact, the authors found the coefficient 2 + 𝑒(cid:3092) as might be an indication of an error in their analysis of energy balance. This error was found by Levine [9] and Sorbie[70], while none of the authors above referred to the error associated with Equation (1.32). Levy[15] pointed out that, Batten[57] extended the energy balance equation to capillary penetration through porous media, but he made the same sign error resulting in an equation similar to that of Szekely. This error had not noticed in several publication[15], might be because 26 of Szekely took 𝑒(cid:3092) = 0.45 as mentioned earlier, so, they got (2 + 𝑒(cid:3092) = 2.45) which is the same order of 2.24 ≤ 𝑚 ≤ 2.41. Figure 1.6 illustrates the penetration rate of water in a vertical 0.133 mm radius capillary tube. A. Experimental data; B. Model of Szekely et al. [1971] given by (2.34) (m = 2.45); C. Lucas- Washburn model[53] Used with permission of Elsevier Science & Technology Journals, from Jeje, Ayodeji A. "Rates of spontaneous movement of water in capillary tubes." Journal of Colloid and Interface Science 69, no. 3 (1979): 420-429; permission conveyed through Copyright Clearance Center, Inc A conclusion can be presented to clarify the effects of reservoir inertia on the penetration into a capillary tube. Within the early stages of the infiltration at 𝑧(cid:3030) is of the order of 𝑟(cid:2868), for the inertia term 𝑧(cid:3030) should be replaced by 𝑧(cid:3030) + 𝜆𝑟(cid:2868). Then, it can be noticed that the reservoir inertia has small influence on the overall behavior of the infiltration, as illustrated in Figure 1.6 [53], Equation (1.34) makes a little improvement of Lucas-Washburn’s equation (1.4), but it still overestimates the initial penetration rate[15]. 27 1.7.3 Couette Correction Levine et al. [9] criticized Szekely et al.[8] approach and they suggested to add a term 𝑚(cid:4593) 𝑅𝑒⁄ mentioned by Oka[66] to Hagenbach correction 𝑚, to be 𝑚 + 𝑚(cid:4593) 𝑅𝑒⁄ , where 𝑚(cid:4593) is a constant, and 𝑅𝑒 is a Reynolds number given by 𝑅𝑒 = 2𝜌𝑟(cid:2868) 𝜇 𝑑𝑧(cid:3030) 𝑑𝑡 Plugging the correction term to Equation (1.21), then the total drag force given by 𝑓(cid:3030)(cid:3031) = 1 2 (cid:4678)𝑚 + 𝑚(cid:4593) 𝑅𝑒 (cid:4679) 𝜌𝜋𝑟(cid:2868) (cid:2870) (cid:3436) 𝑑𝑧(cid:3030) 𝑑𝑡 (cid:2870) (cid:3440) (1.35) (1.36) The term 𝑚(cid:4593) 𝑅𝑒⁄ called Couette correction, and it is an important in measuring a viscosity using rheometer[15][66]. In creeping flow at low Reynolds, additional pressure drop can be counted associated with viscous energy dissipation at the end of tube equivalent to an increase in capillary tube effective length. Using Equation (1.35) and replacing Equation (1.36) in Equation (1.27) in case of ℎ = 0 & 𝜑 = 𝜋 2⁄ , the following expression is found [15] 𝜌(𝑧(cid:3030) + 𝜆𝑟(cid:2868)) 𝑑(cid:2870)𝑧(cid:3030) 𝑑𝑡(cid:2870) + 8𝜇 (cid:2870) (cid:4678)𝑧(cid:3030) + 𝑟(cid:2868) 𝑚(cid:4593) 32 𝑟(cid:2868)(cid:4679) 𝑑𝑧(cid:3030) 𝑑𝑡 = 2𝜎 cos 𝜃 𝑟(cid:2868) − 𝜌𝑔𝑧(cid:3030) − 𝑚 2 𝜌 (cid:3436) (cid:2870) (cid:3440) 𝑑𝑧(cid:3030) 𝑑𝑡 (1.37) Equation (1.37) gives the fact that the supplementary viscous drag is due the change in capillary length from 𝑧(cid:3030) to 𝑧(cid:3030) + 𝑚(cid:4593)𝑟(cid:2868) 32⁄ within the range in which viscous forces apply. As indicated by Sylvester and Rosen[71], Equation (1.36) has been found by superposition of the Hagenbach correction obtained for the inertia-dominant flow, and the Couette correction obtained for the viscous flow, and it was not obtained directly by solving laminar flow equation[15]. 28 Experimental studies on viscous flow through an orifice and through short tubes, measured the values of 𝑚(cid:4593) = 36.7 ± 0.6 [72][73]. These results have a good agreement with that found theoretically for viscous flow through orifice 𝑚(cid:4593) = 12𝜋 ≅ 37.7[74]. Based on the results found by Bond[72][73], Weissberg[75] has used infinite reservoir and derived the upper limit value for 𝑚(cid:4593) = 43.6 by suggesting this value is independent of capillary tube length 𝑙. However, an experimental study has done by Astarita and Greco[76] obtained a large value of 𝑚(cid:4593) = 795, and it suggested that 𝑚(cid:4593) is a very sensitive to the contraction geometry[15]. Their regression analysis was later criticized by Sylvester and Rosen[71] who have obtained 𝑚(cid:4593) = 295 ∓ 50, and they suggested 𝑚(cid:4593) to be decreased with decrease the surface area ratio 𝑟(cid:2868) (cid:2870) 𝑟(cid:3045) (cid:2870)⁄ , 𝑟(cid:3045) is the radius of the reservoir [15]. Consider a very small capillary tube radius compared to reservoir dimensions, 𝑚(cid:4593) is on the order of 10(cid:2870) by taking in account the suggestions made by Sylvester and Rosen[71]. Levy[15] pointed out from Equation (1.37), that increase in capillary tube effective length should not exceed a few radii. In addition, for reservoir inertia correction 𝜆𝑟(cid:2868), Couette correction only valid to early stage of capillary penetration, at large Reynold number and the inertia dominant[15]. This lead that this correction has no significant effects of the theoretical predictions of capillary infiltration, but it may be significant at low initial velocity[15]. The above conclusion was supported by Levine et al.[9] who obtained Equation (1.37) using different approach. They build their analysis applying Navier-Stokes equations on capillary tube with parabolic velocity distribution. Excluding the vena contracta energy loss term, Levine et al.[9] derived an equation identical to Szekely, used combination of Equations (1.29) & (1.30). As what Szekely did for finding 𝑝(cid:2869)(𝑧(cid:3030)), Levine et al.[9] required to come up with an expression describes inlet pressure 𝑝(cid:2869)(𝑧(cid:3030)). Rather than using Szekely approach using energy conservation 29 equation, Levine et al.[9] derived their model of the reservoir pressure far field 𝑟 ≥ 𝑟(cid:2868), from Navier-Stokes equations using reservoir radial velocity given by Equation (1.31). Then, Levine et al.[9] needed to combine the far field pressure with inlet pressure, which it has been done using a momentum conservation equation for the near field system of hemisphere at 𝑟 ≤ 𝑟(cid:2868) as illustrated in Figure 1.5. The near field velocity was unknown, they had to use an approximation approach to find the rate of change of total momentum in the system [9]. They consider the viscous forces acting along the surface 𝑟 = 𝑟(cid:2868), then they were able to obtain an expression covering both the reservoir inertia and the Couette correction. The model was derived by Levine et al.[9] is identical to Equation (1.37) but with 𝜆 = 37 36⁄ , 𝑒(cid:3092) = 2.33, and 𝑚(cid:4593) = 8. Levy[15] believed that Levine et al.[9] should not have neglected the convective inertia term in their analysis of the far field pressure and should have instead obtained 𝑒(cid:3092) = 2.58 as been discussed earlier in this chapter. The value of 𝑚(cid:4593) obtained by authors showed that Couette correction for the effective viscous length is in the order of less than 𝑟(cid:2868), but for accurate study should not be neglected especially in capillary penetration flow from finite reservoir, for instance, dynamics of droplet penetration. 1.7.4 Conclusion on Pressure Drops due to Capillary Tube-Reservoir Interactions Flow in capillary tube faces drag forces at the entrance region, that can potentially act to reduce the penetration rate. for this reason, three kinds of pressure losses effects have been studied. First, pressure drop due to sudden contraction is proportional to square of penetration rate, effects of energy loss have been incorporated to the Hagenbach correction. Second, unsteady penetration rate has been applied, so the reservoir inertia promotes pressure drop, in the amount equal to the increase of infiltration length in the order of 𝑟(cid:2868) within the region in which inertia force is applied. Third, pressure drop due to viscous dissipation at the entrance of 30 capillary tube with very low Reynolds number, incorporated to the Couette correction, is equal to the increase in the penetration length in the order of one or several 𝑟(cid:2868) within the region in which inertia force is applied. It can be concluded that all the three corrections have failed to justify discrepancies between measurements and predictions reported in the literature. Couette corrections are only considerable at the initial stages of the capillary penetration at 𝑧(cid:3030) → 𝑟(cid:2868). After initial stages and depending on interface velocity (contact line), Hagenbach correction might still have significant influence to reduce the penetration rate and need to be calculated in the prediction analysis [15]. 1.8 Penetration Model with consideration of Meniscus Effects 1.8.1 Overview The flow in the vicinity of the advancing meniscus is associated with two problems. First; the Hagen-Poiseuille’s flow has been characterized to be a fully developing flow “parabolic velocity profile”, where meniscus displacement appears to contradict the well-known no-slip boundary condition at the wall of the capillary tube; in which the speed of the fluid at the wall is supposed to be zero[15][77]. Second, dependence on the contact angle on the velocity. As will be discussed in this section, these two issues are intimately connected. In general, these studies are part of comprehensive topic of liquid spreading on solid surfaces, which have been extensively studied [79][80][81][83][85], and has attracted the attention of many investigators both from a theoretical and a practical point of view[90][15]. Many issues remain unresolved that are related to sub-microscopic mechanisms by which a liquid displaces another fluid from a solid surface, especially flow behavior from finite reservoir; water drops as example[15]. 31 1.8.2 Static Contact Angle and Contact Angle Hysteresis Contact angle is located at the point in which the three interfaces of liquid, solid and gas come in contact together. This angle is associated with properties of liquid and solid, and interaction and repulsion forces between liquid and solid. Those forces known as cohesion and adhesion forces which are intermolecular forces. Contact angle can be defined in two cases; static contact angle 𝜃(cid:3046) at which the meniscus at rest intersects the solid at the junction of the three distinct phases, liquid, fluid (can be gas or different liquid) and solid, while dynamic contact angle 𝜃(cid:3031) at the movement of the system. The point at three phases that come in contact together is commonly referred to as the contact line or wetting line[15] Figure 1.7 shows the static contact angle and contact line at Solid, Liquid, and Gas interfaces It can be see that Figure 1.7 illustrates the thermodynamic equilibrium conditions of Liquid- Solid-Gas interfaces have been defined using Young’s equation, that relates static contact angle 𝜃(cid:3046), and surface tensions of liquid σ, liquid-solid 𝜎(cid:3013)/(cid:3046), and gas-solid 𝜎(cid:3008)/(cid:3046). The solid surface assumed to be perfectly flat and homogeneous, Young’s equation wrote as 𝜎 cos 𝜃(cid:3046) = 𝜎(cid:3008)/(cid:3046) − 𝜎(cid:3013)/(cid:3046) (1.38) 32 Even though Equation (1.38) deucedly simple, but it eluded to be experimentally verified due to an incapability of measuring 𝜎(cid:3008)/(cid:3046) and 𝜎(cid:3013)/(cid:3046) [79]. In addition, for the system of Liquid-Gas- Solid, the contact angle is a unique property according to Young’s equation, however practically usually possible examine the range of contact angle for Liquid-Gas-Solid and Liquid-Liquid- Solid [79][91]. The term of hysteresis or contact angle hysteresis it related to solid surface impurities, heterogeneity, and surface roughness[80][88][15]]. Equation (1.38) still can be apply at the microscopic scale of surface roughness, that make possibility of multiple equilibrium configurations in existence of impurities and heterogeneities of surfaces, which can be lead to observe the hysteretic behavior at the macroscopic scale[15]. For the reasons that cause hysteresis phenomenon, it has been found that the contact line is stuck on the solid surface, discontinued at a particular location on the solid surface. This phenomenon is not only for 𝜃 = 𝜃(cid:3046), but whenever 𝜃 sits within a finite interval around 𝜃(cid:3046) [15][80] 𝜃(cid:3045) ≤ 𝜃 ≤ 𝜃(cid:3028) (1.39) where 𝜃(cid:3045) is static receding contact angle and 𝜃(cid:3028) is a static advancing contact angle. static advancing contact angle and static receding contact angle. As shown in Figure 1.8.a, the receding contact angle, 𝜃(cid:3045), is the smallest contact angle can be reached before the wetting line begins to move toward the wetting phase. In contrast, as illustrated in Figure 1.8.b, the static advancing contact angle, 𝜃(cid:3028), is the largest contact angle can be reached before the wetting line begins to move toward the non-wetting phase. Contact angle usually measured once contact line move[88]. 33 Figure 1.8 shows the contact angle hysteresis: a. Static receding contact angle; b Static advancing contact angle For any normal surfaces that have not been specifically treated, the range of difference between 𝜃(cid:3028) − 𝜃(cid:3028) ≥ 10(cid:2868) in the case of Liquid-Gas system[80]. While for Liquid- Liquid systems 𝜃(cid:3028) − 𝜃(cid:3028) ≫ 10(cid:2868). In case of Glycerin-Silicon oil interfaces in horizontal capillary tubes of radius 𝑟(cid:2868) = 1𝑚𝑚, Fermigier and Jenffer[91] observed the magnitude of the interval 𝜃(cid:3028) − 𝜃(cid:3028)was of the order of 60(cid:2868). 1.8.3 Dynamic Contact Angle For liquid displacement by any of Liquid-Fluid-Solid system during flow in capillary channels, the contact line at the Liquid-Fluid interface has observed to move. Measuring dynamic contact angle 𝜃(cid:3031) can be done at the macroscopic scale in the order of few microns, the common measurement techniques have been taken place at the apparent interface intersects the solid surface[83][92]. The variables associate with the dynamic contact angle and influence on it, have been studied under different geometrical configurations as illustrated in Figure 1.9, for instance, these methods include the spreading of drops on a solid surface as shown in Figure 1.9.a [93][94], the flow in capillary tubes of circular cross section; Figure1.9.b[15][91][95][96][97], the wetting lines formed by immersing (withdrawal) a plat; 34 Figure 1.9.c [85][87][93], and the rotation of a horizontal cylinder in a pool of liquid; Figure 1.9.d[15]. Figure 1.9 shows the dynamic contact angles in different geometries used to study them: Spreading drops, Liquid-Fluid displacement in capillary tubes, Steady immersion (withdrawal) Rotation of a horizontal cylinder in a pool of liquid The shape of the interface (contact line) of forced air displacement by various non-volatile liquids (i.e. imbibition) have been examined by Hoffman[95] using capillary tube; 𝑟(cid:2868) = 1 𝑚𝑚 [Levy]. Forced meniscus motion has imposed externally to be at constant velocity, while spontaneous displacement where the wetting process is inherently transient[15]. Under condition of perfect wetting at 𝜃(cid:3046) = 0, Hoffman[95] found that the apparent dynamic contact angle significantly associated with the capillary number at the system controlled by the dominant of interfacial and viscous forces, the capillary number defined as 𝐶𝑎 = 𝜇 𝑢(cid:2868) 𝜎 (1.40) where 𝑢(cid:2868) is interface displacement velocity. Hoffman[95] has found a relationship between apparent contact angle 𝜃(cid:3040) and capillary number plus a shift factor 𝜇 𝑢(cid:2868) 𝜎⁄ + 𝐹(𝜃(cid:3046)), 𝐹(𝜃(cid:3046)) is 35 shift factor, it is very small ≈ 0, generally a constant for any Liquid-Solid-Gas system. This relationship has depicted on a curve can be called a comprehensive curve that characterizes the shape of any liquid-air interface in a motion in the case of only viscous and interfacial forces are important[95]. Figure 1.10 illustrates the effects of flow on apparent contact angle of advancing liquid-air interface[95]Used with permission of Elsevier Science & Technology Journals, from Hoffman, Richard L. "A study of the advancing interface. I. Interface shape in liquid—gas systems." Journal of colloid and interface science 50, no. 2 (1975): 228-241; permission conveyed through Copyright Clearance Center, Inc As illustrated in Figure 1.10, Hoffman[95] examined the relation between capillary number plus shift function ranging from 10(cid:2879)(cid:2873) to around 36, corresponding to apparent contact angle ranging from 0 to 180(cid:2868). Other contribution forced displacement data has obtained by both of Fermigier and Jenffer[91], and Hansen and Toong[98] using similar Hoffman’s technique[95], and other researchers have used alternative experimental techniques that found close fit to 36 comprehensive curve, validate that dynamic contact angles are independent of the flow geometry and measurement method[92][93]. 1.9 Finite Flow in capillary tube Lucas-Washburn’s model is a good approximation in case of a quasi-steady-state laminar flow of a Newtonian flow form infinite reservoir, or a large droplet is deposited on a capillary tubes/pore media, they neglected the inertia effect and viscous drag [99]. However, for flow from finite reservoir/ small drops, the volume and contact angle effects have to be considered especially when it is comparable with the capillary size, and the surface curvature and surface tension may generate a capillary pressure that promotes droplet penetration[100]. The penetration dynamics of a small droplet into capillary tube was first reported by Marmur[99][101][102]. Marmur had studied the effects of surface curvature and surface tension and incorporated these effects into analysis, when a small droplet meets the surface edge of capillary tube, it forms two curvatures at the air-water-solid interfaces: droplet surface outside the pore and the meniscus inside. Marmur found that for non-wetting droplet may penetrate into capillary tube if the initial droplet radius 𝑅(cid:2868) satisfy: 𝑅(cid:2868) < 𝑟(cid:2868) − cos 𝜃(cid:3046) (1.41) where 𝑟(cid:2868) is the capillary tube radius, and 𝜃(cid:3046) is a static contact angle which is 𝜃(cid:3046) > 90(cid:2868) in non- wetting case, while in a wetting case of 𝜃(cid:3046) > 90(cid:2868) the penetration rate is higher in compare with the flow from infinite reservoir. Droplet penetration is driven by the capillary pressure differences between droplet pressure and meniscus in the pore. Laplace pressure of the droplet varies inversely with the droplet radius, the process of penetration enhanced with increasing the pressure differences. In liquid penetration, both the pore radius and droplet volume determine the 37 dynamics behavior of penetration, so to understand and ability to control this behavior needed to increase[103]. Although other forms of the equation were obtained by Bell and Cameron[40] and West[41], but they did not started using Hagen–Poiseuille equation that have been derived independently by Jean Léonard Marie Poiseuille in 1838 and Gotthilf Heinrich Ludwig Hagen and published by Poiseuille in 1840–41 and 1846 [42]. The theoretical justification of the Poiseuille law was given by George Stokes in 1845[43]. Lucas and Washburn both had started from the Hagen-Poiseuille’s law, which relates the steady-state flow of a liquid through a capillary tube[44] 𝑄 = (cid:2872) 𝜋𝑟(cid:2868) 8𝜇𝑙 ∆𝑝 (1.42) Where 𝑄 is the volumetric flow of liquid, 𝑟(cid:2868) and 𝑙 are the radius and the length of the capillary tube, respectively, 𝜇 is the dynamic viscosity of the liquid, and ∆𝑝 is the total driving pressure acting to force the liquid along the capillary. Based on the Equation (1.42), Lucas and Washburn obtained[15] 𝑑𝑧(cid:3030) 𝑑𝑡 = (cid:2870) 𝑟(cid:2868) 8𝜇 ∆𝑝(𝑧(cid:3030)) 𝑧(cid:3030) (1.43) where 𝑧(cid:3030) is the penetrated length of liquid displacing air inside a capillary tube at time 𝑡. Assuming the parabolic profile of the velocity distribution of Hagen-Poiseuille’s flow, 𝑑𝑧(cid:3030) 𝑑𝑡⁄ is a mean velocity across a section of capillary tube, and displaced air (gas phase) velocity has been assumed to be negligible with respect to the viscosity of the liquid[15]. In this study, an analysis of finite fluid flow at low Reynolds for entry flow is presented in which the penetration towards the capillary entrance from a finite reservoir (water droplet or film). Time-dependent solutions of the Navier-Stokes equations are crucial for understanding the 38 rate of liquid penetration in various fluid dynamics scenarios. The Navier-Stokes equations, which describe the motion of viscous fluid substances, are a set of nonlinear partial differential equations that account for the conservation of momentum and mass. When solving these equations for time-dependent problems, we often use numerical methods to handle the complexity of the equations and the boundary conditions. In the context of liquid penetration, these solutions help model how a liquid moves through porous media or capillary tubes over time. The rate of penetration is influenced by factors such as viscosity, surface tension, and the dynamic contact angle. By solving the Navier-Stokes equations, we can predict how quickly a liquid will penetrate a given material, which is essential for applications in fields like material science, biology, and engineering. The next chapter will delve deeper into these solutions, exploring specific mathematical techniques and their applications. It will also discuss how initial and boundary conditions are set up to ensure accurate and stable solutions, and how these solutions can be validated against experimental data 39 CHAPTER 2: GOVERNING EQUATIONS AND MATHEMATICAL FORMULATION 2.1 Overview In this study, an analysis of finite fluid flow at low Reynolds for entry flow is presented in which the penetration towards the capillary entrance from a finite reservoir (water droplet) is considered. Boundary layer theory and time-dependent solutions of Navier-Stokes equations are describing the rate of droplet penetration. A theoretical analysis of various cases is performed, resulting in ordinary differential equations that can be solved relatively rapidly. The cases studied include:  A single droplet entering a pore  A single droplet attached to the surface entering a pore with stick-slip motion of the contact line  A single droplet covering several pores with a receding contact line  Infinitely large and finite film of liquid both drained by a single pore  Infinitely large and finite film drained by several pores Governing equations are derived for the penetration length and the changes in pore or film geometry as the fluid enters the pore(s). Parametric studies are performed to understand the effects of various properties on the solution. The current study presents a novel analysis focusing on the dynamic aspects of droplet behavior into a pore and provides a theoretical study of droplet penetration into a pore using Levin-Szekely theories. The departure from Poiseuille flow in the capillary near the entrance and in the vicinity of the moving meniscus appears to be amenable to mathematical treatment. However, it is noted that the asymptotic solutions provided by Washburn-Rideal and Levin- Szekely are not valid for very short contact times and for flow from a finite reservoir (water 40 droplet), which may typically be encountered in certain printing and lithographic operations. It is a real possibility that a rigorous formulation of the problem would be required for the description of these processes. 2.2 Mathematical Formulation An attempt has been made in this study using N-SEs in combined with boundary layer model during the penetration of the droplet into the pore to obtain a theoretical solution of the dynamic of the droplet. Figure 2.1 shows the system of penetration of droplet with constant initial radius 𝑅(cid:3036) into pore with radius 𝑎 2.2.1 Governed equation of region 1 Continuity and Navier-Stokes equations for incompressible fluid: 𝜕 𝜕𝑧 (𝑟𝑢) + 𝜕 𝜕𝑟 (𝑟𝑣) = 0 (2.1) 41 𝜕𝑢 𝜕𝑡 + 𝑢 𝜕𝑢 𝜕𝑧 + 𝑣 𝜕𝑢 𝜕𝑟 = − 1 𝜌 𝜕𝑃 𝜕𝑧 + 𝜈 (cid:4680) 𝜕(cid:2870)𝑢 𝜕𝑟(cid:2870) + 1 𝑟 𝜕𝑢 𝜕𝑟 + 𝜕(cid:2870)𝑢 𝜕𝑟(cid:2870)(cid:4681) + 𝑔 (2.2) Multiplying by r and using continuity Equation (2.1): 𝜕 𝜕𝑡 (𝑟𝑢) + 𝜕 𝜕𝑧 (𝑟𝑢(cid:2870)) + 𝜕 𝜕𝑟 (𝑟𝑢𝑣) = − 𝑟 𝜌 𝜕𝑃 𝜕𝑧 + 𝜈 (cid:4680) 𝜕 𝜕𝑟 (cid:3436)𝑟 𝜕𝑢 𝜕𝑟 (cid:3440) + 𝜕(cid:2870)𝑢 𝜕𝑧(cid:2870) (cid:4681) + 𝑔𝑟 (2.3) Equation (2.3) can be solved by integration from 𝑟 = 0 𝑡𝑜 𝑟 = 𝑎 , using no-slip boundary condition 𝑢(𝑟 = 𝑎) = 0, and Equation (2.1) with the known of: 1. Volume flux along the tube 𝜋𝑎(cid:2870) (cid:3031)(cid:3053) (cid:3031)(cid:3047) (cid:3028) = 2𝜋 ∫ 𝑟𝑢 𝑑𝑟 (cid:2868) 2. Assuming Poiseuille flow, 𝑢 = 2 (cid:4672)1 − (cid:3045)(cid:3118) (cid:3028)(cid:3118)(cid:4673) (cid:3031)(cid:3053) (cid:3031)(cid:3047) which satisfies the continuity equation 1 3. 𝑢 independent of z, and 𝑑𝑢 𝑑𝑟⁄ ](cid:3045)(cid:2880)(cid:3028) = (− 4 𝑎)⁄ 𝑑𝑧 𝑑𝑡⁄ 4. Poiseuille flow 𝑣 = 0, 𝑝(𝑧, 𝑟, 𝑡) = 𝑝(𝑧, 𝑡) = 𝑝(cid:3028)(cid:3047)(cid:3035) + 2𝛾 𝑅⁄ ; p is independent of 𝑟, R is droplet radius, and atmospheric pressure 𝑝(cid:3028)(cid:3047)(cid:3035) 5. Integrating Equation (2.3) from 𝑧 = 0 to 𝑧 = 𝑧(𝑡) Equation (2.3) now yields 1 2 𝑎(cid:2870)𝑧 𝑑(cid:2870)𝑧 𝑑𝑡(cid:2870) = − 1 𝜌 (cid:4680) 1 2 𝑎(cid:2870)(𝑝(cid:3028)(cid:3047)(cid:3035) + 2𝜎 𝑅⁄ ) − (cid:3505) 𝑟 𝑝(0, 𝑡)𝑑𝑟 (cid:4681) − 4𝒱𝑧 (cid:2868) (cid:3028) 𝑑𝑧 𝑑𝑡 (2.4) − 1 2 𝑔𝑎(cid:2870)𝑧 Figure 2.1 illustrates the system where a droplet with a constant initial radius 𝑅 penetrates a pore with radius a. To solve the Equation (2.4), it is necessary to evaluate the integral related to the pressure at the origin, 𝑧 = 0, denoted as 𝑝(0, 𝑡), as shown in Figure 2.1 2.2.2 Governed equation of region 2 The pressure can be found by introducing region (2) of spherical polar coordinates 𝑅, 𝜃, ∅ as illustrating in Figure 2.2 42 Figure 2.2 shows the system of penetration of droplet with constant initial radius 𝑅(cid:3036) into pore with radius 𝑎 To evaluate 𝑝(0, 𝑡) at the following boundary condition: 1. 𝑝(𝑅, 𝑡) = 𝑝(cid:3028)(cid:3047)(cid:3035) 2. 𝑝(𝑎, 𝑡): will be derived later in Equation2.10, and 3. Velocity field in the region 𝑅 ≤ 𝑎 To evaluate 𝑝(𝑎, 𝑡) by using governed equations in 𝑅 direction as following 𝜕(𝑅(cid:2870)𝜐(cid:3019)) 𝜕𝑅 = 0 (2.5) 𝜌 𝜕𝜐(cid:3019) 𝜕𝑡 + 𝜌𝜐(cid:3019) 𝜕𝜐(cid:3019) 𝜕𝑅 = − 𝜕𝑃 𝜕𝑅 + 𝜇 (cid:4684) 𝜕 𝜕𝑅 (cid:3437) 1 𝑅(cid:2870) 𝜕 𝜕𝑅 (𝑅(cid:2870)𝜐(cid:3019))(cid:3441)(cid:4685) + 𝜌𝑔 (2.6) For region (2), can be solved the momentum equation of droplet by applying the volume conservation between droplet and flow into pore, gets: 4 3 𝜋(𝑅(cid:3036) (cid:2871) − 𝑅(cid:2871)) = 𝜋𝑎(cid:2870)𝑧 Differentiate Equation (2.7) with respect to time 𝑡, where 𝑅(cid:3036) is a constant, gets: −4𝜋𝑅(cid:2870) 𝑑𝑅 𝑑𝑡 = 𝜋𝑎(cid:2870) 𝑑𝑧 𝑑𝑡 ∶ 𝜐(cid:3019) = 𝑑𝑅 𝑑𝑡 and (2.7) (2.8) 43 𝜐(cid:3019) = 𝑑𝑅 𝑑𝑡 = − 𝑎(cid:2870) 4𝑅(cid:2870) 𝑑𝑧 𝑑𝑡 Applying Equation (2.9) in Equation (2.6) and integrate from 𝑎 to 𝑅: 𝑝(𝑎, 𝑡) = 𝑝(𝑅, 𝑡) + 𝜌𝑎(cid:2870) 4 (cid:3428) 1 𝑅 − 1 𝑎 (cid:3432) 𝑑(cid:2870)𝑧 𝑑𝑡(cid:2870) Equation (2.10) satisfy the condition of incompressible flow ∇(cid:2870)𝑃 = 0 (2.9) (2.10) The velocity field in the region 𝑅 ≤ 𝑎 is unknown. Levine [9] introduced an innovative technique to apply for the momentum balance equation to address this unknown for the system in the hemisphere from 0 to 𝑅 ≤ 𝑎 as following 𝑟𝑎𝑡𝑒 𝑜𝑓 𝑐ℎ𝑎𝑛𝑔𝑒 𝑜𝑓 𝑡𝑜𝑡𝑎𝑙 𝑚𝑜𝑚𝑒𝑛𝑡𝑢𝑚 𝑖𝑛 𝑠𝑦𝑠𝑡𝑒𝑚 = 𝑓𝑙𝑢𝑥 𝑜𝑓 𝑚𝑜𝑚𝑒𝑛𝑡𝑢𝑚 𝑒𝑛𝑡𝑟𝑖𝑛𝑔 − 𝑓𝑙𝑢𝑥 𝑜𝑓 𝑚𝑜𝑚𝑒𝑛𝑡𝑢𝑚 𝑙𝑒𝑎𝑣𝑖𝑛𝑔 (2.11) + 𝑠𝑢𝑚 𝑜𝑓 𝑓𝑜𝑟𝑐𝑒𝑠 𝑎𝑐𝑡𝑖𝑛𝑔 𝑜𝑛 𝑡ℎ𝑒 𝑠𝑦𝑠𝑡𝑒𝑚 The force at 𝑅 = 𝑎 in the 𝑧 direction (cid:3095) (cid:2870) 𝑓(cid:2869) = 2𝜋𝑎(cid:2870) (cid:3505) −𝜎(cid:3019)(cid:3019) sin 𝜃 cos 𝜃 𝑑𝜃 (cid:2868) The stress tensor at the hemisphere at 𝑅 = 𝑎 𝜎(cid:3019)(cid:3019) = −𝑝(𝑎, 𝑡) + 2𝜇 (cid:3428) 𝜕𝜐(cid:3019) 𝜕𝑅 (cid:3432) (cid:3019)(cid:2880)(cid:3028) , 𝜎(cid:3019)(cid:3087) = 0 (2.12) (2.13) Plug Equation (2.13) and Equation (2.10) in Equation (1.12) and integration over the hemisphere gets 𝑓(cid:2869) = 𝜋𝑎(cid:2870) (cid:4678)𝑝(𝑅, 𝑡) + 𝜌𝑎(cid:2870) 4 (cid:3428) 1 𝑅 − 1 𝑎 (cid:3432) 𝑑(cid:2870)𝑧 𝑑𝑡(cid:2870) − 𝜇 𝑎 𝑑𝑧 𝑑𝑡 (cid:4679) The corresponding fore term at 𝑧 = 0 in the entry region of the capillary tube is (cid:3028) 𝑓(cid:2870) = −2𝜋 (cid:3505) 𝑟 𝑝(0, 𝑡)𝑑𝑟 (cid:2868) 44 (2.14) (2.15) The flux of momentum in 𝑧 direction entring the hemisphere at 𝑅 = 𝑎 is 𝐹𝑜𝑚 = 2𝜋𝜌𝑎(cid:2870) (cid:3505) (𝜐(cid:3019) (cid:3095) (cid:2870) (cid:2868) (cid:2870))(cid:3019)(cid:2880)(cid:3028) sin 𝜃 cos 𝜃 𝑑𝜃 = 1 4 𝜋𝜌𝑎(cid:2870) (2.16) Levien [9] assumed the momentum leaving at 𝑧 = 0 as following (cid:3028) 2𝜋𝜌 (cid:3505) 𝑟 𝑢(cid:2870) 𝑑𝑟 (cid:2868) = 4 3 𝜋𝜌𝑎(cid:2870) (cid:3436) 𝑑𝑧 𝑑𝑡 (cid:2870) (cid:3440) (2.17) The mean acceleration inside the hemisphere system was approximated by Levine [9] as following 𝑚𝑒𝑎𝑛 𝑎𝑐𝑐𝑒𝑙𝑒𝑟𝑎𝑡𝑖𝑜𝑛 𝑜𝑓 𝑡ℎ𝑒 𝑠𝑦𝑠𝑡𝑒𝑚 = 2 3 𝜋𝜌𝑎(cid:2871) (cid:4678) 19 24 𝑑(cid:2870)𝑧 𝑑𝑡(cid:2870) + 1 8𝑎 (cid:3436) 𝑑𝑧 𝑑𝑡 (cid:2870) (cid:3440) (cid:4679) (2.18) plug Equations (2.14),(2.15),(2.16),(2.17), and (2.18), into Equation (2.11) of momentum balance gets: (cid:3028) (cid:3505) 𝑟 𝑝(0, 𝑡) 𝑑𝑟 = (cid:2868) 𝑎(cid:2870) 2 𝑝(𝑅, 𝑡) + 𝜌𝑎(cid:2871) (cid:3436) 𝑎 8𝑅 − 7 18 (cid:3440) 𝑑(cid:2870)𝑧 𝑑𝑡(cid:2870) − 𝑎𝜇 2 𝑑𝑧 𝑑𝑡 − 7 12 𝜌𝑎(cid:2870) (cid:3436) 𝑑𝑧 𝑑𝑡 (cid:2870) (cid:3440) (2.19) Using Equation (2.19) in Equation (2.4) with mathematical arrangement, the equation of motion of droplet entering the pore obtained as following: (cid:4680)𝑧 − 𝑎(cid:2870) 4𝑅 + 7𝑎 9 (cid:4681) 𝑑(cid:2870)𝑧 𝑑𝑡(cid:2870) − 7 6 (cid:3436) 𝑑𝑧 𝑑𝑡 (cid:2870) (cid:3440) − 𝜇 𝑎(cid:2870)𝜌 [𝑎 + 4𝑧] 𝑑𝑧 𝑑𝑡 − 1 𝜌 [∆𝑃 − 𝜌𝑔𝑧] = 0 (2.20) Where ∆𝑃 is a total pressure on the droplet with equal to 2𝜎 𝑅⁄ . The dimensionless equation of motion is: (cid:3428)𝑧∗ − 𝐶𝑎 𝑅𝑒 8 𝑊𝑒 + 7 (cid:3432) 𝑉∗(cid:4593) − 9 7 6 𝑉∗(cid:2870) − 2𝐹𝑟 𝑅𝑒 [1 + 4𝑧∗] 𝑉∗ − (cid:4680) 2 𝐹𝑟(cid:2870) 𝑊𝑒 − 𝑧∗(cid:4681) = 0 (2.21) 45 Where 𝑉∗(cid:4593) = (cid:3031)(cid:3118)(cid:3053)∗ (cid:3031)(cid:3047)∗(cid:3118) , 𝑉∗ = (cid:3031)(cid:3053)∗ (cid:3031)(cid:3047)∗, Reynolds number, 𝑅𝑒 = (cid:2870)(cid:3028)(cid:3096)(cid:3022) (cid:3091) , Froude number, 𝐹(cid:3019) = (cid:3022) √(cid:3034) (cid:3028) , (cid:3091)(cid:2905) (cid:3097) , Weber number, 𝑊𝑒 = (cid:3096)(cid:2905)(cid:3118)(cid:3019) (cid:3097) , 𝑧 = 𝑎𝑧∗, (cid:3031)(cid:3053)∗ (cid:3031)(cid:3053) = (cid:2869) (cid:3028) , 𝑡 = 𝑡∗(cid:3495) (cid:3028) (cid:3034) , (cid:3031)(cid:3047)∗ (cid:3031)(cid:3047) = Capillary number, 𝐶𝑎 = (cid:3495) (cid:3034) (cid:3028) 46 CHAPTER 3: MATHEMATICAL RESULTS OF DYNAMIC BEHAVIOR OF DROPLETS ENTERING CAPILLARY CHANNELS The behavior of water droplets entering capillary channels has been studied, and presented a developed mathematical model of the finite flow. In this chapter, the penetration rate and penetration distance of the liquid into the pore was captured as a function of time. The penetration results of this mathematical model are validated with previous mathematical models derived by Washburn and Szekely. To validate the results, experimental data has been plotted against the results of mathematical model. The objective was to develop a mathematical model that could predict the droplet at which liquid penetrated the surface pores. 3.1 Mathematical Validation h t g n e l n o i t a r t e n e P 1 0.8 0.6 0.4 0.2 0 0 Present work Szekely Washburn 0.2 0.4 Time 0.6 0.8 1 Figure 3.1 shows the plot of the penetration length vs. time for the results of present work, Szekely, and Washburn, for the droplet radius R=1 mm, pore radius a=0.1 mm, and ΔP=50/a 47 When a water droplet of radius 𝑅 = 1𝑚𝑚 was loaded onto a capillary tube of radius 𝑎 = 0.1𝑚𝑚, and constant physical properties of water. The data results was non-dimensionalized and depicted in the curves Figure 3.1and Figure 3.2. Water droplet was gradually penetrated the pore, by the inspection the curves of droplet penetration length against time, and droplet penetration rate against the time, derived Equation (2.21) differ from Washburn equation in very initial stages of capillary penetration at initial time period of the order of 10-3-10-5 sec. Equation (2.21) plotted against Washburn Equation (1.4) and modified Szekely Equation (1.29). It can be noticed that the higher values of (cid:4672) (cid:3091) (cid:3028)(cid:3118)(cid:3096) [𝑎 + 4𝑧](cid:4673) in Equation (2.21), the steeper of the penetration rate curve, and the shorter is of the time period in which Washburn equation does not applicable. In addition, can be seen an close agreement with modified Szekely Equation (1.29) and follow the same curve behavior but different results due to Szekely solution was derived for flow from infinite reservoir in compare with the current study for flow of finite droplet. e t a r n o i t a r t e n e P 1 0.8 0.6 0.4 0.2 0 0 Droplet Penetration Rate at R=1 mm, a=0.1 mm & ΔP Present work Szekely Washburn 0.2 0.4 Time 0.6 0.8 1 Figure 3.2 shows the plot of the penetration rate vs. time for the results of present work, Szekely, and Washburn, for the droplet radius R=1 mm, pore radius a=0.1 mm, and ΔP=50/a 48 3.2 Experimental Validation To validate the results, experimental data [20] has been plotted against the results of current mathematical model as illustrated on Figure 3.3. The flow dynamics resulted from mathematical model follows the same behavior of flow from experimental results but with significant errors due to the impact of droplet deformation fall from ℎ distance. While the current results evaluated at ℎ = 0 with no deformation, so the differences between curves are due to the droplet deformation due to dynamic pressure. Figure 3.3 shows the plot of the comparison of droplet penetration rate and length of mathematical results vs experimental results Attempt has been made using try and error to estimate the dynamic pressure, Equation (3.1) is a modified mathematical model, but it fits only the current case with the inputs utilized to solve the mathematical model. The comparison results to the experimental results are plotted on Figure 3.4. the close agreement between penetration length curves of mathematical results against the experimental results which validate the mathematical model. 49 1 4 (cid:4680)𝑧 − 𝑎(cid:2870) 4𝑅 + 7𝑎 9 (cid:4681) − 𝑑(cid:2870)𝑧 𝑑𝑡(cid:2870) − 1 𝜌 7 30 1 4 (cid:3428)∆𝑃 + 𝐾𝐸 − 𝜌𝑔𝑧(cid:3432) = 0 (cid:3436) (cid:2870) (cid:3440) 𝑑𝑧 𝑑𝑡 − 2𝜇 𝑎(cid:2870)𝜌 [𝑎 + 4𝑧] 𝑑𝑧 𝑑𝑡 (3.1) where 𝐾𝐸 = 0.5𝜌𝑉 ∗ 𝑣(cid:2870); and 𝑉 & 𝑣 is droplet volume and velocity respectively Figure 3.4 shows the plot of the droplet penetration Rate & Length of Mathematical results (include dynamic pressure) vs, Experimental results Table 3.1 Experimental and mathematical results of droplet penetration length and penetration rate 50 3.3 Transient dynamics of liquids entering a pore (pores) Cases studied Figure 3.5 shows the behavior of droplet entering a pore and dynamics of contact line motion This section explores the behavior of liquids as they enter pores, focusing on the transient dynamics involved. Several cases are studied to understand the various scenarios and their implications as shown in Figure 3.5 and Figure 3.6: 1. Single Droplet Entering a Pore: This case examines how a single droplet of liquid penetrates a pore, considering factors such as the droplet's size, the pore's geometry, and the interaction between the liquid and the pore walls. 2. Single Droplet with Stick-Slip Motion: Here, the focus is on a droplet attached to a surface that enters a pore with stick-slip motion of the contact line. This involves analyzing the intermittent movement of the droplet as it adheres to and slips along the surface. 3. Droplet Covering Multiple Pores: This scenario studies a droplet that covers several pores with a receding contact line. It looks at how the liquid distributes itself among multiple pores and the dynamics of the contact line as it moves. 51 Figure 3.6 shows the behavior of droplet versus film, entering network of pores and dynamics of contact line motion 4. Finite Films Drained by a Single Pore: This case investigates the behavior of finite films of liquid as they are drained by a single pore. It considers the rate of drainage and the changes in the film's geometry. 5. Finite Films Drained by Multiple Pores: Similar to the previous case, but with multiple pores involved. This examines how the liquid film interacts with several pores simultaneously and the resulting dynamics The above-mentioned five cases have been studied and applied to the new mathematical model. This mathematical model does not cover the dynamics of spreading, receding, and stick- slip contact line motion. These dynamics will be studied separately using simulation analysis by ANSYS FLUENT. The mathematical model is valid for cases A, D, and F as shown on Figure 3.6, which will be discussed in more detail in the next section. While cases B, C, and E will be studied numerically. 3.3.1 Case A: Droplet set on an edge Figure 3.7 shows the system of droplet is setting on an edge and entering a pore, the volume conservation between droplet and flow into pore can be determined using volume conservation equation. Equation (2.21) has been solved using various parameters as discussed in the next section. 52 𝜽𝒅 𝑹𝒅 𝟐𝒂 Figure 3.7 shows the droplet set on an edge entering a pore 3.3.1.1 Effect of droplet size entering the pore, 𝒂 = 𝒄𝒐𝒏𝒔𝒕𝒂𝒏𝒕 1 ) m c ( h t g n e L n o i t a r t e n e P 0.8 0.6 0.4 0.2 0 0 Penetration Length at various droplet radius at constant pore radius a=0.01 (cm) Rd=0.01 (cm) Rd=0.02 (cm) Rd=0.05 (cm) Rd=0.1 (cm) Rd=Infinity 4 8 Time (ms) 12 16 20 Figure 3.8 shows the plot of the penetration Length at various droplet radius at constant pore radius a=0.01 (cm) Penetration length and penetration rate have been plotted using Equation (2.21). Figure 3.8 and Figure 3.9 illustrate dynamics behavior of different droplet volumes. High penetration length and high penetration rate happened at lower droplet diameter with is consisted with Laplace pressure that increase with decrease the droplet diameter. 53 The blue curve (𝑅(cid:3031) = 0.01 𝑐𝑚) in Figure 3.8 shows the highest penetration length over time, indicating that smaller droplets penetrate deeper into the pore. As the droplet radius increases, such as in the black curve (𝑅(cid:3031) = 0.1 𝑐𝑚), the penetration length decreases. This decrease in penetration length for larger droplets is due to lower Laplace pressure. Smaller droplets experience higher Laplace pressure, which drives them deeper into the pore. This is consistent with the blue curve showing the highest penetration length. Conversely, larger droplets have lower Laplace pressure, resulting in less penetration, as reflected in the black and red curves. c e s / m c e t a R n o i t a r t e n e P 80 70 60 50 40 30 20 10 0 Penetration Rate at various droplet radius at constant pore radius a=0.01 (cm) Rd=a=0.01 (cm) Rd=0.02 (cm) Rd=0.05 (cm) Rd=0.1 (cm) Rd=Infinity 0 2 4 6 8 10 Time (ms) 12 14 16 18 20 Figure 3.9 shows the plot of the penetration Rate at various droplet radius at constant pore radius a=0.01 (cm) The dynamics behavior of penetration rates shows that higher penetration rates are achieved with smaller droplet diameters due to higher Laplace pressure, as clearly demonstrated by the blue curve in Figure 3.9. This consistency with Laplace pressure, where the pressure increases as the droplet diameter decreases, drives the higher penetration rates observed in the plot. Understanding these dynamics is crucial for applications such as microfluidics, where precise control of droplet behavior in small-scale systems is essential; oil recovery, where enhancing 54 fluid penetration efficiency in porous media is vital; and biomedical engineering, where targeted drug delivery through porous tissues relies on these dynamics. 3.3.1.2 Effect of variant of pore number on droplet entering the pore, 𝒂&𝑹 = 𝒄𝒐𝒏𝒔𝒕𝒂𝒏𝒕 The mathematical model is valid for flow of finite liquid into network of pores, can be the case of track etched membrane. Figure 3.10 and Figure 3.11 show the slight differences in the droplet dynamics behavior. In case of multiple pores, the Laplace pressure decreases, and accordingly, the penetration length and penetration rate will take more time than in case of one pore. c e s / m c e t a R n o i t a r t e n e P 70 60 50 40 30 20 10 0 0 2 4 6 8 10 12 14 16 18 20 Time (ms) One Pore n Pores Figure 3.10 shows the plot of the penetration rate vs. time for the results of present work, for the droplet radius R=1 mm, pore radius a=0.1 mm, and ΔP=50/a, for one pore vs. n pores 55 ) m c ( h t g n e L n o i t a r t e n e P 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 One Pore n Pores 0 2 4 6 8 10 Time (ms) 12 14 16 18 20 Figure 3.11 shows the plot of the penetration rate vs. time for the results of present work, for the droplet radius R=1 mm, pore radius a=0.1 mm, and ΔP=50/a, for one pore vs. n pores 3.3.1.3 Effect of variant of pore number on film entering the pore, 𝒂 =constant, flow from infinite reservoir Filme Penetration Rate at one pore vs. n pores ) c e s / m c ( e t a R n o i t a r t e n e P 70 60 50 40 30 20 10 0 One Pore n Pores 0 20 40 60 80 100 Time (ms) 120 140 160 180 200 Figure 3.12 shows the plot of the film penetration rate versus time for the results of present work, comparing one pore versus multiple pores 56 Filme Penetration Length at one pore vs. n pores ) m c ( h t g n e L n o i t a r t e n e P 4 3.5 3 2.5 2 1.5 1 0.5 0 One Pore n Pores 0 20 40 60 80 100 Time (ms) 120 140 160 180 200 Figure 3.13 shows plot of the film penetration length versus time for the results of the present work, comparing one pore versus multiple pores Figure 3.13 depicts the penetration length of a finite film flow over time within a pore system, the graph contrasts the penetration length between a single pore and multiple pores, underscoring the differences in fluid dynamics between these systems. This comparison is vital for comprehending fluid behavior in porous materials. When finite film flows through a pore and a network of pores, it adheres to the Laplace pressure law, which describes the pressure difference across the curved surface of a droplet due to surface tension. In this scenario, the film is modeled as a spherical cap of a droplet with a large radius. The spherical cap assumption simplifies the complex interactions at the fluid interface, allowing for more accurate predictions of fluid behavior under varying conditions. The penetration rate, as shown in Figure 3.12, is lower during finite film flow through a network of pores compared to a single pore. This is due to the increased resistance and complex pathways within the network, which affect the fluid's movement and distribution. The Brinkman model, an extension of the classical Darcy model, can be applied to account for viscous phenomena in these porous flow systems. This model helps in 57 understanding the multiphase flow applications and the locally mass-conserving methods that are crucial for accurate simulations. 3.3.1.4 Effect of variant of contact angle θ, 𝒂&𝑹 =constant The effects of contact angle through the study of capillary penetration are considerable at the very initial stages of penetration processes which contradict Szekely [1] findings that stated that meaningful results of penetration processes can be taken only about 1-2 sec after the initiation of flow. Figure 3.14 and Figure 3.15 illustrated the significant effects of contact angle on the penetration rate and length. Lower contact angle increases the hydrophilic properties of the pore surface and vice versa. Droplet Penetration Length at R=1 mm, a=0.1 mm 0.3 0.2 0.1 ) m m ( h t g n e L n o i t a r t e n e P 0 0 θ=30 θ=45 θ=60 θ=80 10 20 30 40 50 60 70 80 Time (ms) Figure 3.14 shows the plot of the penetration length versus time at different contact angles for droplet radius R=1 mm, and pore radius a=0.1 mm 58 Droplet Penetration Rate at R=1 mm, a=0.1 mm θ=30 θ=45 θ=60 θ=80 ) m m ( h t g n e L n o i t a r t e n e P 15 10 5 0 0 10 20 30 40 50 60 70 80 Time (ms) Figure 3.15 shows the plot of the penetration rate versus time at different contact angles for a water droplet with radius R=1 mm, and pore radius a=0.1 mm 3.3.1.5 Effect of varying Weber numbers Figure 3.16 presents the droplet dynamics at various Weber numbers (We = 0.07, 0.14, 0.7, 1.4, and 6). The spreading stage started at initial contact times, with the increase of the Weber number, the contact circuit speed of the droplet becomes lower, resulting in a lower Weber number at the same time, and the maximum penetration rate increases. 59 * v 14 12 10 8 6 4 2 0 200 150 100 * z 50 0 We=6 1.4 0.7 0.14 0 20 40 60 80 100 120 t* Figure 3.16 shows the plot of the effect of varying Weber numbers on the penetration length and penetration rate versus time for a water droplet at ρ=1kg/m3, µ=0.01 mPa s, σ=72 dyne/cm, Re=20, Fr=3.2, Ca=1.4E-03, a=0.01cm, U=10 cm/sec 3.3.1.6 Creeping flow Reynolds number Low Reynolds number flow, known as creeping flow or Stokes flow, at the Reynolds number 𝑅𝑒 ≪ 1. Figure 3.17 shows creeping flow at very low Reynolds number, which in this study has been taken as 𝑅𝑒 = 0.03. In Stokes flow, the viscous forces are higher than the inertial forces. However, in the droplet flow into a pore, its initial speed is significant enough that acceleration and inertia cannot be negligible compared to the fluid's viscosity. In infinite flow, when the 60 Reynolds number Re << 1, inertial effects can be ignored, and only viscous resistance is considered. Which is not the case in the finite flow of spherical volumes. 8 6 * v 4 2 0 100 80 60 40 20 0 * z 0 100 200 t* 300 400 Figure 3.17 shows the plot of the creeping flow Reynolds number for a water droplet at R=0.1cm, a=0.015cm, very low Re=0.03, We=14E-08,Ca=139E-08, Fr=2.61E-03 3.3.1.7 Effect of varying liquid properties The properties of five liquids, in addition to water, have been studied with a constant Froude number, as presented in Figure 3.18. It can be noticed that at the initial contact time, the penetration rate and length of the water droplet are higher compared to the ethanol droplet due to the effects of viscous forces and inertial forces. Glycerin, blood, oil, and methanol exhibit similar behavior when compared to each other, with glycerin and blood showing one pattern, and oil and methanol showing another. 61 * v 20 18 16 14 12 10 8 6 4 2 0 Water: Re=2, Fr=10, Ca=1.4E-03, We=0.35 Glycerin 50%water: Re=0.4, Fr=10, Ca=1.08E-02, We=0.51 Blood: Re=0.76, Fr=10, Ca=4.9E-03, We=0.47 Oil: Re=1.14, Fr=10, Ca=5.6E-03, We=0.8 Alcohol Ethanol: Re=13, Fr=10, Ca=5.4E-04, We=0.9 Alcohol methanol: Re=2.7, Fr=10, Ca=2.6E-03, We=0.9 * z 9 8 7 6 5 4 3 2 1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t* Figure 3.18 shows the plot of the effect of varying the penetration length and penetration rate vs time for varying liquid properties at R=0.25cm, a=0.001cm, U=10 cm/sec 3.3.1.8 Effect of variant Froude number Figure 3.19 presents the droplet dynamics at various Froude numbers (Fr = 3.2, 10, and 32). It is noticed that an anomaly region does not follow the properties of the liquid during spreading stage at initial contact time, to the time of maximum penetration rate. With the increase of the Froude number after around one time unite, the penetration length decrease due to the effect of inertia forces. 62 * v 16 12 8 4 0 0 Fr=3.2 10 32 45 30 15 * z 2 4 6 8 10 12 0 14 t* Water at ρ=1kg/m3, µ=0.01 mPa s, σ=72 dyne/cm, U=10 cm/sec, We=14E-2, CA=14E-5 Figure 3.19 shows the plot of the effect of varying Froude number on the penetration length and penetration rate vs time for a water droplet at ρ=1kg/m3, µ=0.01 mPa s, σ=72 dyne/cm, U=10 cm/sec, We=14E-2, CA=14E-5 3.3.1.9 Effect of variant of Reynolds The importance of varying Reynolds numbers on finite flow has been studied and is presented in Figure 3.20. It illustrates the behavior of water droplet penetration rate and penetration length at different Reynolds numbers (5.7, 11.5, 17, 30, and 40). The water droplet flow rate and length increase with increasing Reynolds numbers. 63 40 30 * v 20 10 0 0 600 500 400 * z 300 200 100 0 0 Re=5.7 Re=11.5 Re=17 Re=30 Re=40 Re=5.7 11.5 17 30 40 15 30 45 t* 15 t* 30 Figure 3.20 shows the plot of the effect of varying Reynolds on the penetration length and penetration rate vs time for a water droplet at ρ=1kg/m3, µ=0.01 mPa s, σ=72 dyne/cm Capillary numbers and Weber numbers are the main dimensionless numbers that are used to evaluate the effect of surface tension on the flow. Figure 3.21 shows the droplet penetration rate at varied Capillary number. It noticed that with different calculated capillary numbers, there is no noticeable difference in the droplet dynamics. Since in finite flow at 𝑅𝑒 ≪ 1 that applied in this 64 study, it can be concluded to consider the importance of Weber number rather than capillary number * z 18 16 14 12 10 8 6 4 2 0 0.0014 0.01076 0.0049 0.00565 0.000536 0.00262 0 5 10 15 20 25 30 35 40 45 t* Figure 3.21 shows the plot of the effect of varying the Capillary number on the penetration length and penetration rate vs time for a water droplet at ρ=1kg/m3, µ=0.01 mPa s, σ=72 dyne/cm 3.3.2 Application of the mathematical model on a pore of elliptical cross section A theoretical analysis of finite liquid entering a pore with circular cross section vs entering a pore with elliptical cross section are presented in Figure 3.22 and Figure 3.23. It was shown that in the case of equal cross section areas but different shapes (circular and elliptical) the initial penetration rate resulted through finite flow through a pore with elliptical cross section is faster than in case of circular cross section, but the dynamics behavior changes after around 2ms, so the penetration rate decreases rapidly with elliptical cross section while circular cross section increases. 65 0.8 0.4 ) m c ( h t g n e L n o i t a r t e n e P Circular cross-section Elliptical cross-section circular cross-sections a=0.012 (cm) Capillary pressure using Hydraulic diameter, a=8E-3, b=18E-3 (cm) Capillary pressure using elliptic radiuses, a=8E-3, b=18E-3 (cm) 0 0 5 10 15 Time (ms) 20 25 30 Figure 3.22 shows the plot of the effect of varying cross-section on the penetration length and penetration rate versus time for a water droplet at Circular cross-section versus Elliptical cross- section The initial penetration rate for the circular cross-section (black curve) starts at around 60 𝑐𝑚/𝑠𝑒𝑐 as can be seen in Figure 3.23. The peaks slightly above this value, and then gradually decreases over time to about 10 𝑐𝑚/𝑠𝑒𝑐 at 30 𝑚𝑠, indicating higher initial penetration rates compared to elliptical cross-sections. The elliptical cross-sections (red and blue curves) start lower than the black curve, peak below it, and gradually decrease over time but remain below the black curve throughout the duration. The peak penetration rate for the black curve is higher than the red and blue curves, suggesting that the circular cross-section allows for a higher maximum penetration rate. All curves show a gradual decrease in penetration rate over time, indicating that the driving force, likely Laplace pressure, diminishes as the droplets continue to penetrate the pore. 66 circular cross-sections a=0.012 (cm) Capillary pressure using Hydraulic diameter, a=8E-3, b=18E-3 (cm) Capillary pressure using elliptic radiuses, a=8E-3, b=18E-3 (cm) ) c e s / m c ( e t a R n o i t a r t e n e P 70 60 50 40 30 20 10 0 0 5 10 15 Time (ms) 20 25 30 Figure 3.23 shows the plot of the effect of varying cross-section on the penetration rate and penetration rate versus time for a water droplet at Circular cross-section versus Elliptical cross- section 3.4 Conclusion The calculations presented in this study confirm the appropriateness of the modified Szekely equation and Washburn equation, as an asymptotic solution, for understanding the behavior of capillary penetration. However, this study presented a solution that is valid for very short initial times, which may typically be found in certain printing and lithographic operations. 67 CHAPTER 4: GOVERNING EQUATIONS AND NUMERICAL SIMULATIONS FOR SIMULATING A DROPLET ENTERING A PORE 4.1 Overview In this chapter, the governing equations of Navier-Stokes and the interface tracking method have been studied. Computational simulation methods were used to obtain 2D results from the model using the ANSYS-FLUENT 23 R2. Mesh development for 2D model and its refinement is considered. The computational methods used to simulate droplet dynamics, and the behavior of Stick-slip regimes. User defined function (UDF) and Volume of Fluid (VOF) algorithms are used in the solver to capture the interface between the liquid and the gas. Figure 4.1 shows the governing parameters stick-slip regimes, and receding contact angle Capillary forces, which depend on surface tension and contact angle, drive the movement of droplets in both slip and stick regimes. In the stick regime, pinning occurs when the droplet's contact line temporarily sticks to surface features, affecting its movement and potentially causing hysteresis in the contact angle. The dynamic contact angle changes as the droplet moves, 68 influencing the capillary forces and the overall behavior of the droplet. Understanding these parameters and their interactions is crucial for applications in microfluidics, coating processes, and surface science. Figure 4.1 highlights how different parameters influence the behavior of droplets in stick-slip regimes and their receding contact angles. Understanding these interactions is crucial for applications in surface science and fluid dynamics. 4.2 Mathematical Formulation Laminar, Incompressible, and Newtonian flow has been assumed to solve the continuity equation and conservative of momentum. The two phases’ properties are constant, and the droplet assumed to be spherical setting on the edge. The computational simulation of the droplet setting on a pore is performed using Volume-of-Fluid (VOF) method for Incompressible, and Newtonian two-phase flow. The governing equations of the conservation of mass is: ∇. 𝑉(cid:4652)⃗ = 0 ∂𝑉(cid:4652)⃗ 𝜕𝑡 + 𝑉(cid:4652)⃗. ∇𝑉(cid:4652)⃗ = − 1 𝜌 ∇p + 1 𝜌 ∇. τ + 1 𝜌 𝐹⃗(cid:3020)(cid:3007) + 1 𝜌 𝐹⃗(cid:3003) (4.1) (4.2) Where 𝑉(cid:4652)⃗ is velocity vector, 𝑡 is time, p is the pressure, 𝜌 and τ are the fluid density and share stress tensor, respectively. 𝐹⃗(cid:3020)(cid:3007) is the surface tension force at the Gas-Liquid interface, 𝐹⃗(cid:3003) is all body forces that are acting on the droplet. The governing equations are discretized on a Eulerian grid with structured uniform mesh size. The set-up solutions use Eulerian frame of references to couple the solution with methodology of interface tracking. The Volume-of-Fluid (VOF) method is used in ANSYS-FLUENT 23R2 to track the interface. The phase fraction 𝛼 = 1 is designated for cells filled with full liquid, and 𝛼 = 1 for cells filled with only gas. 0 < 𝛼 < 1 describes the cells where the interface is located. The Volume-of-Fluid (VOF) method is presented by the equation: 69 ∂α 𝜕𝑡 + (𝑉(cid:4652)⃗. ∇)𝛼 = 0 (4.3) The density of the liquids in each computational grid has been calculated based of the fluid fraction in the cell as following: 𝜌(cid:3041) (cid:3030)(cid:3032)(cid:3039)(cid:3039) = 𝛼 × 𝜌(cid:3013)(cid:3036)(cid:3044)(cid:3048)(cid:3036)(cid:3031) + (1 − 𝛼) × 𝜌(cid:3008)(cid:3028)(cid:3046) (4.4) Adding Continuum surface tension force (CSF) to the VOF calculations that were implemented in FLUENT 23R2, and the pressure at the Liquid-Gas interface as following: ∇𝑃 = 𝑃(cid:3013)(cid:3036)(cid:3044)(cid:3048)(cid:3036)(cid:3031) + 𝑃(cid:3008)(cid:3028)(cid:3046) = 𝜎𝑘 (4.4) Where ∇𝑃 is interfacial pressure differences that explain in Equation (1.1). The effect of the surface tension was calculated based on the nondimensional numbers described in Equation (2.21). In this study, for low Rynolds flow, Weber number considered instead of capillary number. The effect of contact angle on Stick-slip regimes has been considered by considering wall adhesion angle in FLUENT 23R2. Advancing static contact angle was applied for the numerical simulation for low Rynolds flow. User Defined Function (UDF) was applied to track the droplet height and the radius of contact circuit. 4.3 Two-Dimensional Solver Methods Solver method in ANSYS-FLUENT using Least Squares Cell-Based Gradient evaluation to discretize the momentum equation. The solution was assumed to be linear. The pressure solution used a spatial discretization scheme PRESTO (Pressure Staggering Option) method. The VOF model for Eulerian multiphase use Geo-Reconstruct scheme. The pressure-based solver is used to solve the two-phase of Liquid-Gas flow modeling using projection method to converge the solution. SIMPLE algorithm is used to solve the relationship 70 between velocity and pressure corrections. The Pressure-Implicit with Splitting of Operators (PISO) is a pressure-velocity coupling scheme as a part of SIMPLE algorithms. 4.4 Geometry and Meshing Numerical simulation was studied on one hole geometry. Numerical code with User Defined Functions (UDF) applications will be validated by modeling liquid droplet entering on capillary hole along the diameter of the droplet while it spreads on the substrate. Initialization method using UDF to enforce the initial spherical phase of liquid into the cubic geometry with structured mesh as shown in Figure 4.2. The Cube dimensions is 4 × 4 × 4mm, contain a droplet with 2mm diameter enforce to enter a 0.1 mm diameter of a pore. Meshing structure used local refinement with resolution ∆𝑥 = 17 × 10(cid:2879)(cid:2874)𝑚 that result 124968 elements and time resolution ∆𝑡 = 10(cid:2879)(cid:2874) 𝑠𝑒𝑐. Two boundary conditions were applied, no-slip boundary condition and pressure-outlet boundary conditions Figure 4.2 illustrates the numerical geometry 71 CHAPTER 5: SIMULATION RESULTS AND DYNAMIC BEHAVIOR OF DROPLETS ENTERING CAPILLARY CHANNELS 5.1 Numerical simulation comparison and validation Numerical simulation of water droplet of a radius 𝑅 = 1 𝑚𝑚 entering onto a capillary tube of radius 𝑎 = 0.1 𝑚𝑚, and constant physical properties of water, was conducted. The data results were depicted in Figure 5.1. The water droplet gradually penetrated the pore. By inspection the curves of droplet penetration length against time, a close agreement with modified Szekely Equation (1.29) can be seen, following the same curve behavior but yielding different results due to the Szekely solution being derived for flow from an infinite reservoir compared with the current study of flow of a finite droplet. 0.4 0.3 ) c e S m / ( 0.2 0.1 Pentration rate Penetration length 0 0 0.01 0.02 0.03 0.04 0.06 0.07 0.08 0.09 0.1 0.05 (Sec) 0.008 0.006 0.004 ) m ( 0.002 0 Figure 5.1 shows the plot of the simulation results for droplet penetration rate and penetration length at zero contact angles, with a water droplet radius of R=1 mm, and pore radius of a=0.1 mm 72 Equation (2.21) was tested against simulation results. When a water droplet with radius 𝑅 = 1𝑚𝑚 was loaded onto a capillary tube with radius 𝑎 = 0.1𝑚𝑚, maintaining constant physical properties of water, the data results are shown in the curves of Figure 5.2 and Figure 5.3. Droplet Penetration Rate at a=0.1mm & R=1mm 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 ) c e s / m ( e t a R n o i t a r t e n e P Mathematical Model Mathematical model Simulation study Simulation Study 0 0 0.005 0.01 Time (Sec) 0.015 0.02 Figure 5.2 shows the plot of the penetration rate versus time for a water droplet with radius R=1 mm, and a pore radius a=0.1 mm, comparing both mathematical results (include dynamic pressure) and numerical simulation results Figure 5.2 shows a water droplet gradually penetrating a pore with a radius of 𝑎 = 0.1𝑚𝑚. By inspecting the curves of droplet penetration rate against time, it is evident that the derived Equation (2.21) differs from the simulation results due to assumptions made during the simulation process. Additionally, there is a close agreement with the simulation results, following the same curve behavior but yielding different results due to errors caused by these assumptions. 73 8% e m u l o v 4% t e l p o r d l a t o t / e m u l o v n o i t a r t e n e P 0% 0 Droplet Penetration at a=0.1mm & R=1mm Simulation Study Simulation study Mathematical model Mathematical Model Time (Sec) 0.02 Figure 5.3 shows the plot of the droplet penetration volume of mathematical results (include dynamic pressure) versus simulation results for a water droplet with radius of R=1mm, and a pore radius of a=0.1mm Figure 5.3 shows the normalized total penetrated volume with respect to the droplet volume. The results indicate that the depth of penetration volume increases with time for both the simulation and the current results, following the same behavior. However, discrepancies between the simulation results and the current results are due to errors caused by assumptions made during the simulation process. On the other hand, the current study demonstrates greater accuracy in comparison to the Szekely equation and the Washburn equation, as shown in Figure5.3 74 5.2 Droplet Stick-slip behavior The effect of contact angle on stick-slip regimes has been considered by incorporating the wall adhesion angle in FLUENT 23R2. An advancing static contact angle was applied for the numerical simulation of low Reynolds flow. A User Defined Function (UDF) was used to track the droplet height and the radius of the contact area. 5.2.1 Oscillation of droplet height and contact circuit stick – slip behavior θ=0 Droplet height Droplet Hight Radius of contact line Radius of contact line 0.002 0.0015 ) m ( 0.001 0.0005 0 0 0.02 0.04 0.06 0.08 0.1 0.12 (Sec) Figure 5.4 shows the plot of the oscillation of droplet height and contact circuit diameter at zero contact angles for a water droplet with a radius R=1 mm, and pore radius a=0.1 mm Figure 5.4 represents the behavior of droplet height oscillation during droplet penetration into the pore. The curve shows the stick-slip motion through the oscillations of droplet height and contact circuit diameter. Changes in oscillation magnitude can be seen in the damping behavior, which is due to surface energy dissipation affecting droplet dynamics. The droplet height and contact circuit diameter oscillate and eventually reach stability at 0.1mm. In fact, 0.1mm is the pore diameter at the last point where the UDF can track the droplet height and contact circuit 75 diameter. The results of these oscillations validate the theory of stick-slip motion as the droplet penetrates into the capillary channel. 5.2.2 Oscillation of droplet height at θ=0 vs θ=900 0.002 0.0015 ) m ( 0.001 0.0005 0 0 Droplet hight at zero contact angle Droplet hight at 90 contact angle 0.02 0.04 (Sec) 0.06 0.08 0.1 Figure 5.5 shows the plot of the oscillation of droplet heights for a water droplet with a radius R=1 mm, and a pore radius of a=0.1 mm at contact angle 00 vs 900 Details of droplet height and contact circuit movement in Figure 5.5 and Figure 5.6, compare the oscillation behavior by changing the contact angle. According to the results of penetration rate in both figures, at θ=00, the penetrated liquid at 1ms is not noticeable compared to θ=900. After t = 1ms, when penetration initiates and it reaches its maximum value, both the decrease in droplet volume and surface energy dissipation affect the oscillation behavior, pushing the droplet into its final damped shape and diminishing in volume. The droplet contact circuit at θ=00 starts at approximately 0.0008 meters and shows fluctuations before stabilizing around 0.0012 meters at 0.04 seconds. This indicates that the droplet's contact area varies initially but eventually reaches a steady state. While the droplet contact circuit at θ=900 starts similarly but rises more sharply, peaking around 0.0014 meters at 76 0.03 seconds, then fluctuates before stabilizing around 0.0012 meters at 0.05 seconds. This suggests that the droplet's contact area increases rapidly before stabilizing. 0.0016 0.0012 ) m ( 0.0008 0.0004 Droplet contact circuit at zero contact angle Droplet Contact Circuit at Zero Contact Angle Droplet contact circuit at 90 contact angle Droplet Contact Circuit at 90 Contact Angle 0 0 0.01 0.02 0.03 0.04 0.06 0.07 0.08 0.09 0.1 0.05 (Sec) Figure 5.6 shows the plot of the oscillation of droplet contact circuits for a water droplet with a radius of R=1 mm, and a pore radius of a=0.1 mm at contact angle of 00 and 900 5.2.3 Droplet penetration Length at θ=0 and θ=90 The effects of contact angle on capillary penetration are considerable during the very initial stages of the penetration process, which contradicts Szekely's findings that meaningful results of penetration processes can only be obtained about 1-2 seconds after the initiation of flow[1]. Figure 5.7 illustrates the significant effects of contact angle on penetration length. A lower contact angle increases the hydrophilic properties of the pore surface, and vice versa. At zero contact angle, the penetration length increases rapidly, reaching approximately 0.006 meters within 0.09 seconds. 77 Figure 5.7 shows the plot of the plot of penetration lengths versus time for a water droplet with a radius of R=1mm and a pore radius of a =0.1mm for the contact angle 00 vs 900 This indicates that the fluid spreads quickly when the contact angle is zero. While at θ=900 contact angle, the penetration length increases more slowly, leveling off around 0.002 meters within the same time frame. This suggests that the fluid spreads less efficiently when the contact angle θ=900. 78 CHAPTER 6: SUMMARY AND CONCLUSION 6.1 Summary and Conclusions This study presents an analysis of finite fluid flow at low Reynolds numbers for entry flow, focusing on the penetration towards the capillary entrance from a finite reservoir (water droplet or film). Time-dependent solutions of the Navier-Stokes equations are solved to determine the rate of liquid penetration. A theoretical analysis of various cases is performed, resulting in ordinary differential equations that can be solved relatively rapidly. The findings from solving the mathematical model are as follows: 1. The mathematical model was plotted against the Washburn and modified Szekely equations. It can be noticed that higher values of (cid:4672) (cid:3091) (cid:3028)(cid:3118)(cid:3096) [𝑎 + 4𝑧](cid:4673) in the mathematical model result in a steeper penetration rate curve and a shorter time period during which the Washburn equation is not applicable. Additionally, there is close agreement with the modified Szekely equation, following the same curve behavior but yielding different results. This discrepancy is due to the Szekely solution being derived for flow from an infinite reservoir, whereas the current study focuses on the flow of a finite droplet. 2. The theoretical analysis of various cases, in addition to numerical simulation, resulted in close agreement with the modified Szekely equation, following the same curve behavior but yielding different results. This discrepancy is due to the Szekely solution being derived for flow from an infinite reservoir, whereas the current study focuses on the flow of a finite droplet. 3. The calculations presented in this study confirm the appropriateness of the modified Szekely equation and the Washburn equation, as asymptotic solutions, for understanding the behavior of capillary penetration. However, this study presents a solution that is valid 79 for very short initial contact times, which may typically be found in certain printing and lithographic operations. 4. A theoretical analysis of finite liquid entering a pore with a circular cross section versus entering a pore with an elliptical cross section was conducted. It was shown that, in the case of equal cross-sectional areas but different shapes (circular and elliptical), the initial penetration rate through a pore with an elliptical cross section is faster than through a circular cross section. However, the dynamic behavior changes after approximately 2 ms, with the penetration rate decreasing rapidly for the elliptical cross section while increasing for the circular cross section. In the case of finite film flowing through a pore and a network of pores, it follows Laplace pressure law. Here, the film is assumed to be a spherical cap of a droplet with a large radius. A lower penetration rate occurs during finite film flow through a network of pores. 5. The results of the current study were validated using experimental data. The experimental data were plotted against the mathematical results, and the flow dynamics from the mathematical model followed the same behavior as the flow from the experimental results. However, significant errors were observed due to the impact of droplet deformation from a fall height ℎ distance. The current results were evaluated at ℎ = 0 with no deformation, so the differences between the curves are attributed to droplet deformation caused by dynamic pressure. 6. The properties of five liquids, in addition to water, have been studied at a constant Froude number. It was noticed that at the initial contact time, the penetration rate and length of the water droplet were higher compared to the ethanol droplet due to the effects of 80 viscous and inertial forces. Glycerin, blood, oil, and methanol follow the same behavior when compared to each other, specifically glycerin and blood, or oil and methanol. 7. The droplet dynamics at various Froude numbers (Fr = 3.2, 10, and 32) have been investigated. An anomaly region was noticed that does not follow the properties of the liquid during the spreading stage at initial contact time, up to the time of maximum penetration rate. With the increase of the Froude number after approximately one time unit, the penetration length decreases due to the effect of inertial forces. 8. Capillary numbers and Weber numbers are the main dimensionless numbers used to evaluate the effect of surface tension on the flow. The study showed with different calculated capillary numbers, there is no noticeable different of the droplet dynamics. Since in finite flow at 𝑅𝑒 ≪ 1 as applied in this study, it can be concluded that the Weber number is more important than the capillary number. 9. The mathematical model was tested at low Reynolds number flow, known as creeping flow or Stokes flow, at 𝑅𝑒 ≪ 1. The present work studied creeping flow at a very low Reynolds number, specifically 𝑅𝑒 = 0.03. In Stokes flow, viscous forces are higher than inertial forces, but in droplet flow into a pore, the initial speed is significant enough that acceleration and inertia cannot be negligible compared to the fluid's viscosity. In infinite flow, when Re << 1, inertial effects can be ignored, and only viscous resistance is considered. This is not the case in the finite flow of spherical volumes. Computational simulation methods were used to obtain 2D results from the model using ANSYS-FLUENT 23 R2. Mesh development for the 2D model and its refinement were considered. Computational methods were used to simulate droplet dynamics and the behavior of 81 stick-slip regimes. User Defined Function (UDF) and Volume of Fluid (VOF) algorithms were used in the solver to capture the interface between the liquid and the gas: 1. Numerical simulation of a water droplet with a radius 𝑅 = 1𝑚𝑚 entering a capillary tube with a radius of 𝑎 = 0.1𝑚𝑚, and constant physical properties of water, showed close agreement with the modified Szekely equation. The simulation followed the same curve behavior but yielded different results because the Szekely solution was derived for flow from an infinite reservoir, whereas the current study focuses on the flow of a finite droplet. The differences from the simulation are due to assumptions made during the simulation process. Additionally, there is close agreement with the simulation results, following the same curve behavior but yielding different results due to errors caused by these assumptions. 2. The results of penetrated volume indicate that the depth of penetration volume increases with time for both the simulation and the current results, following the same behavior. However, discrepancies between the simulation results and the current results are due to errors caused by assumptions made during the simulation process. On the other hand, the current study shows greater accuracy in comparison to the Szekely equation and the Washburn equation 3. The droplet height and contact circuit diameter oscillate and eventually reach stability at 0.1 mm. In fact, 0.1 mm is the pore diameter at the last point where the UDF can track the droplet height and contact circuit diameter. The results of these oscillations validate the theory of stick-slip motion as the droplet penetrates into the capillary channel. 4. The effects of contact angle on capillary penetration are considerable during the very initial stages of the penetration process, which contradicts Szekely's findings that 82 meaningful results of penetration processes can only be obtained about 1-2 seconds after the initiation of flow. 6.2 Future Work Recommendations 1. It is undoubtedly recommended that experimental studies on a droplet setting on a hole substrate be extended to a parallel network of holes, including corner modeling. This includes studying the effect of surface roughness, specifically on corners and edges. 2. Conducting computational thermal studies of finite flow and the network of pores, which change the properties of the fluid and the porous materials. 3. Deriving mathematical model using Navier-Stokes equations including heat flux derivative, to study the effects of temperature on hysteresis in the contact angles through the study of capillary penetration at the very initial stages of penetration processes. 83 BIBLIOGRAPHY [1] Lucas, Richard. "Ueber das Zeitgesetz des kapillaren Aufstiegs von Flüssigkeiten." Kolloid-Zeitschrift 23, no. 1 (1918): 15-22. [2] Washburn, Edward W. "The dynamics of capillary flow." Physical review 17, no. 3 (1921): 273. [3] Rideal, Eric Keightley. "CVIII. On the flow of liquids under capillary pressure." The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 44, no. 264 (1922): 1152-1159. [4] Bosanquet, C. H. "LV. On the flow of liquids into capillary tubes." The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 45, no. 267 (1923): 525-531. [5] Pickett, Gerald. "Rate of Rise of Water in Capillary Tubes." Journal of Applied Physics 15, no. 8 (1944): 623-623. [6] Rense, William A. "Rate of rise of water in capillary tubes." Journal of Applied Physics 15, no. 5 (1944): 436-437. [7] Siegel, Robert. "Transient capillary rise in reduced and zero-gravity fields." (1961): 165- 170. [8] Szekely, Julian, A. W. Neumann, and Y. K. Chuang. "The rate of capillary penetration and the applicability of the Washburn equation." Journal of Colloid and Interface Science 35, no. 2 (1971): 273-278. [9] Levine, S., P. Reed, E. J. Watson, and G. Neale. "A theory of the rate of rise of a liquid in a capillary." Colloid and interface science 3 (1976): 403-419. [10] Levine, Samuel, Paul Reed, Glenn Shutts, and Graham Neale. "Some aspects of wetting/dewetting of a porous medium." Powder Technology 17, no. 2 (1977): 163-181. [11] Levine, Samuel, and Graham H. Neale. "Theory of the rate of wetting of a porous medium." Journal of the Chemical Society, Faraday Transactions 2: Molecular and Chemical Physics 71 (1975): 12-21. [12] Levine, Samuel, James Lowndes, Eric J. Watson, and Graham Neale. "A theory of capillary rise of a liquid in a vertical cylindrical tube and in a parallel-plate channel: Washburn equation modified to account for the meniscus with slippage at the contact line." Journal of Colloid and Interface Science 73, no. 1 (1980): 136-151. [13] Bear, Jacob. Hydraulics of groundwater. Courier Corporation, 2012. [14] Bear, Jacob, and Arnold Verruijt. Modeling groundwater flow and pollution. Vol. 2. Springer Science & Business Media, 1987. 84 [15] Levy, Laurent Claude. "Experimental and theoretical modeling of DNAPL transport in vertical fractured media." PhD diss., Massachusetts Institute of Technology, 2003. [16] Berry, Joseph D., Michael J. Neeson, Raymond R. Dagastine, Derek YC Chan, and Rico F. Tabor. "Measurement of surface and interfacial tension using pendant drop tensiometry." Journal of colloid and interface science 454 (2015): 226-237. [17] Huling, Scott G., and James W. Weaver. Dense nonaqueous phase liquids. Superfund Technology Support Center for Ground Water, Robert S. Kerr Environmental Research Laboratory, 1991. [18] Bedient, Philip B., Hanadi S. Rifai, and Charles J. Newell. Ground water contamination: transport and remediation. Prentice-Hall International, Inc., 1994. [19] Anderson, William. "Wettability literature survey-part 2: Wettability measurement." Journal of petroleum technology 38, no. 11 (1986): 1-246. [20] Hosseini, Saman. "Droplet Impact and Penetration onto Structured Pore Network Geometries." PhD diss., 2015. [21] Sahimi, Muhammad. Flow and transport in porous media and fractured rock: from classical methods to modern approaches. John Wiley & Sons, 2011. [22] Floege, Jurgen, Richard J. Johnson, and John Feehally. Comprehensive Clinical Nephrology E-Book. Elsevier Health Sciences, 2010. [23] Bhasin, Sanjeev Kumar. Pharmaceutical Physical Chemistry: Theory and Practices. Pearson Education India, 2012. [24] Morrow, Norman R., and Geoffrey Mason. "Recovery of oil by spontaneous imbibition." Current Opinion in Colloid & Interface Science 6, no. 4 (2001): 321-337. [25] Pezron, Isabelle, Guillaume Bourgain, and David Quéré. "Imbibition of a fabric." Journal of Colloid and Interface Science 173, no. 2 (1995): 319-327. [26] Tortorich, Ryan P., and Jin-Woo Choi. "Inkjet printing of carbon nanotubes." Nanomaterials 3, no. 3 (2013): 453-468. [27] Holley, Brian M., and Amir Faghri. "Permeability and effective pore radius measurements for heat pipe and fuel cell applications." In ASME International Mechanical Engineering Congress and Exposition, vol. 4711, pp. 353-371. 2004. [28] Matthews, G. P., C. W. Watts, D. S. Powlson, J. C. Price, and W. R. Whalley. "Wetting of agricultural soil measured by a simplified capillary rise technique." European Journal of Soil Science 59, no. 4 (2008): 817-823. 85 [29] Li, Chuxin, Haoyu Dai, Can Gao, Ting Wang, Zhichao Dong, and Lei Jiang. "Bioinspired inner microstructured tube controlled capillary rise." Proceedings of the National Academy of Sciences 116, no. 26 (2019): 12704-12709. [30] Liu, Zhanjie, Yifan Wang, Fernando J. Muzzio, Gerardo Callegari, and German Drazer. "Capillary drop penetration method to characterize the liquid wetting of powders." Langmuir 33, no. 1 (2017): 56-65. [31] Fukano, Tohru, and Akira Kariyasaki. "Characteristics of gas-liquid two-phase flow in a capillary tube." Nuclear Engineering and Design 141, no. 1-2 (1993): 59-68. [32] Bressler, R. G., and P. W. Wyatt. "Surface wetting through capillary grooves." (1970): 126-132. [33] Lee, Shong-Leih, and Hong-Draw Lee. "Evolution of liquid meniscus shape in a capillary tube." (2007): 957-965. [34] Perwuelz, Anne, Pascal Mondon, and Claude Caze. "Experimental study of capillary flow in yarns." Textile Research Journal 70, no. 4 (2000): 333-339. [35] Dang-Vu, Trong, and Jan Hupka. "Characterization of porous materials by capillary rise method." Physicochemical problems of mineral processing 39 (2005): 47-65. [36] Okada, Kiyoshi, Shuhei Uchiyama, Toshihiro Isobe, Yoshikazu Kameshima, Akira Nakajima, and Taisuke Kurata. "Capillary rise properties of porous mullite ceramics prepared by an extrusion method using organic fibers as the pore former." Journal of the European Ceramic Society 29, no. 12 (2009): 2491-2497. [37] Roger, Jérôme, M. Avenel, and L. Lapuyade. "Characterization of SiC ceramics with complex porosity by capillary infiltration: Part B–Filling by molten silicon at 1500° C." Journal of the European Ceramic Society 40, no. 5 (2020): 1869-1876. [38] Saguy, I. Sam, Alejandro Marabi, and Rony Wallach. "New approach to model rehydration of dry food particulates utilizing principles of liquid transport in porous media." Trends in Food Science & Technology 16, no. 11 (2005): 495-506. [39] Kornev, Konstantin G., and Alexander V. Neimark. "Spontaneous penetration of liquids into capillaries and porous membranes revisited." Journal of colloid and interface science 235, no. 1 (2001): 101-113. [40] Bell, James Munsie, and F. K. Cameron. "The flow of liquids through capillary spaces." The Journal of Physical Chemistry 10, no. 8 (2002): 658-674. [41] West, Gilbert D. "On the resistance to the motion of a thread of mercury in a glass tube." Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 86, no. 583 (1911): 20-25. 86 [42] Sutera, Salvatore P., and Richard Skalak. "The history of Poiseuille's law." Annual review of fluid mechanics 25, no. 1 (1993): 1-20. [43] Szabó, István. Geschichte der mechanischen Prinzipien: und ihrer wichtigsten Anwendungen. Vol. 32. Springer-Verlag, 2013. [44] Munson, B. R., D. F. Young, T. H. Okiishi, and W. W. Huebsch. "Fundamentals of Fluid Mechanics, 6th SI ed." (2010). [45] Peek Jr, R. L., and D. A. McLean. "Capillary penetration of fibrous materials." Industrial & Engineering Chemistry Analytical Edition 6, no. 2 (1934): 85- 90. [46] Van Brakel, Jaap, and P. M_ Heertjes. "Capillary rise in porous media. Part I. a problem." Powder Technology 16, no. 1 (1977): 75-81. [47] Siebold, Alain, André Walliser, Michel Nardin, Max Oppliger, and Jacques Schultz. "Capillary rise for thermodynamic characterization of solid particle surface." Journal of colloid and interface science 186, no. 1 (1997): 60-70. [48] Siebold, Alain, Michel Nardin, Jacques Schultz, André Walliser, and Max Oppliger. "Effect of dynamic contact angle on capillary rise phenomena." Colloids and surfaces A: Physicochemical and engineering aspects 161, no. 1 (2000): 81-87. [49] Ligenza, Joseph R., and Richard B. Bernstein. "The rate of rise of liquids in fine vertical capillaries." Journal of The American Chemical Society 73, no. 10 (1951): 4636-4638. [50] Quéré, David. "Inertial capillarity." EPL (Europhysics Letters) 39, no. 5 (1997): 533. [51] LeGrand, Elmo J., and William A. Rense. "Data on rate of capillary rise." Journal of Applied Physics 16, no. 12 (1945): 843-846. [52] Calderwood, G. F. N., and E. W. J. Mardles. "8—The Rate of Flow of Liquids into Capillaries Under the Action of Surface Forces." Journal of the Textile Institute Transactions 46, no. 3 (1955): T161-T170. [53] Jeje, Ayodeji A. "Rates of spontaneous movement of water in capillary tubes." Journal of Colloid and Interface Science 69, no. 3 (1979): 420-429. [54] Zhmud, B. V., F. Tiberg, and K. Hallstensson. "Dynamics of capillary rise." Journal of colloid and interface science 228, no. 2 (2000): 263-269. [55] Kornev, Konstantin G., and Alexander V. Neimark. "Spontaneous penetration of liquids into capillaries and porous membranes revisited." Journal of colloid and interface science 235, no. 1 (2001): 101-113. 87 [56] Letelier, Mario FS, Hans J. Leutheusser, and César Z. Rosas. "Refined mathematical analysis of the capillary penetration problem." Journal of Colloid and Interface Science 72, no. 3 (1979): 465-470. [57] Batten Jr, George L. "Liquid imbibition in capillaries and packed beds." Journal of colloid and interface science 102, no. 2 (1984): 513-518. [58] Brun, Edmond A., André Martinot-Lagarde, and Jean Collab MATHIEU. Mécanique des fluides. No. BOOK. Dunod, 1968. [59] Bird, R. B. (1960). New variational principle for incompressible non‐Newtonian flow. The Physics of Fluids, 3(4), 539-541. [60] Brittin, Wesley E. "Liquid rise in a capillary tube." Journal of Applied Physics 17, no. 1 (1946): 37-44. [61] Tietjens, Oskar Karl Gustav, and Ludwig Prandtl. Applied hydro-and aeromechanics: based on lectures of L. Prandtl. Vol. 2. Courier Corporation, 1957. [62] White, Gilbert F. Natural hazards, local, national, global. Oxford University Press, 1974. [63] Lew, H. S., and Y. C. Fung. "On the low-Reynolds-number entry flow into a circular cylindrical tube." Journal of Biomechanics 2, no. 1 (1969): 105-119. [64] Lew, H. S., and Y. C. Fung. "Entry flow into blood vessels at arbitrary Reynolds number." Journal of biomechanics 3, no. 1 (1970): 23-38. [65] Goldstein, S. (Ed.), Modern Developments in Fluid Dynamics: An Account of Theory and Experiment Relating to Boundary Layers, Turbulent Motion and Wakes, 2 vols., Dover, New York, 1965. [66] Oka, S., The principles of rheometry, in Rheology: Theory and Applications, vol. 3, edited by F. R. Eirich, pp. 17-82, Academic, San Diego, Calif., 1960 [67] Morse, Philip M., and H. Freshbach. "Methods of Theoretical Physics, volume 1 and 2." (1953). [68] Batchelor, Cx K., and G. K. Batchelor. An introduction to fluid dynamics. Cambridge university press, 2000. [69] Huang, Wei, Raghbir S. Bhullar, and Yuan Cheng Fung. "The surface-tension-driven flow of blood from a droplet into a capillary tube." J. Biomech. Eng. 123, no. 5 (2001): 446-454.." [70] Sorbie, K. S., Y. Z. Wu, and S. R. McDougall. "The extended Washburn equation and its application to the oil/water pore doublet problem." Journal of colloid and interface science 174, no. 2 (1995): 289-301. 88 [71] Sylvester, N. D., and S. L. Rosen. "Laminar flow in the entrance region of a cylindrical tube: Part II. Non‐Newtonian fluids." AIChE journal 16, no. 6 (1970): 967-972. [72] Bond, W. N. "The effect of viscosity on orifice flows." Proceedings of the Physical Society of London 33, no. 1 (1920): 225. [73] Bond, W. N. "Viscosity determination by means of orifices and short tubes." Proceedings of the Physical Society of London 34, no. 1 (1921): 139. [74] Roscoe, R. "XXXI. The flow of viscous fluids round plane obstacles." The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 40, no. 302 (1949): 338-351. [75] Weissberg, Harold L. "End correction for slow viscous flow through long tubes." The Physics of Fluids 5, no. 9 (1962): 1033-1036. [76] Astarita, Gianni, and Guido Greco. "Excess pressure drop in laminar flow through sudden contraction. Newtonian liquids." Industrial & Engineering Chemistry Fundamentals 7, no. 1 (1968): 27-31. [77] Rapp, Bastian E. Microfluidics: modeling, mechanics and mathematics. William Andrew, 2016. [78] Dussan V, E. B. "Immiscible liquid displacement in a capillary tube: The moving contact line." AIChE Journal 23, no. 1 (1977): 131-133. [79] Dussan, E. B. "On the spreading of liquids on solid surfaces: static and dynamic contact lines." Annual Review of Fluid Mechanics 11, no. 1 (1979): 371-400. [80] De Gennes, Pierre-Gilles. "Wetting: statics and dynamics." Reviews of modern physics 57, no. 3 (1985): 827. [81] Blake, T. D., Dynamic contact angles and wetting kinetics, in Wettability, edited by J.C. Berg, pp. 251-309, Marcel Dekker, New York, 1993. [82] Blake, Terence D., and Kenneth J. Ruschak. "Wetting: static and dynamic contact lines." In Liquid film coating, pp. 63-97. Springer, Dordrecht, 1997. [83] Kistler, STEPHAN F. "Hydrodynamics of wetting." Wettability 6 (1993): 311-430. [84] Berg, John, ed. Wettability. Vol. 49. CRC Press, 1993. [85] Blake, T. D., and K. J. Ruschak. "A maximum speed of wetting." Nature 282, no. 5738 (1979): 489-491. [86] Blake, T. D., and J. M. Haynes. "Kinetics of liquidliquid displacement." Journal of colloid and interface science 30, no. 3 (1969): 421-423. 89 [87] Blake, T. D., M. Bracke, and Y. D. Shikhmurzaev. "Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle." Physics of fluids 11, no. 8 (1999): 1995-2007. [88] Blake, Terence D., and Kenneth J. Ruschak. "Wetting: static and dynamic contact lines." In Liquid film coating, pp. 63-97. Springer, Dordrecht, 1997. [89] Shikhmurzaev, Yu D. "The moving contact line on a smooth solid surface." International Journal of Multiphase Flow 19, no. 4 (1993): 589-610. [90] Shikhmurzaev, Yulii D. "Moving contact lines in liquid/liquid/solid systems." Journal of Fluid Mechanics 334 (1997): 211-249. [91] Fermigier, Marc, and Patrice Jenffer. "An experimental investigation of the dynamic contact angle in liquid-liquid systems." Journal of colloid and interface science 146, no. 1 (1991): 226-241. [92] Ngan, Clifton G. "On the nature of the dynamic contact angle: an experimental study." Journal of Fluid Mechanics 118 (1982): 27-40. [93] Bracke, M., F. De Voeght, and P. Joos. "The kinetics of wetting: the dynamic contact angle." In Trends in Colloid and Interface Science III, pp. 142-149. Steinkopff, 1989. [94] Foister, Robert T. "The kinetics of displacement wetting in liquid/liquid/solid systems." Journal of colloid and interface science 136, no. 1 (1990): 266-282. [95] Hoffman, Richard L. "A study of the advancing interface. I. Interface shape in liquid— gas systems." Journal of colloid and interface science 50, no. 2 (1975): 228-241. [96] Legait, Benoît, and Pierre Sourieau. "Effect of geometry on advancing contact angles in fine capillaries." Journal of colloid and interface science 107, no. 1 (1985): 14-20. [97] Calvo, A., I. Paterson, R. Chertcoff, M. Rosen, and J. P. Hulin. "Dynamic capillary pressure variations in diphasic flows through glass capillaries." Journal of colloid and interface science 141, no. 2 (1991): 384-394. [98] Hansen, R. J., and T. Y. Toong. "Interface behavior as one fluid completely displaces another from a small-diameter tube." Journal of Colloid and Interface Science 36, no. 3 (1971): 410-413. [99] Marmur, Abraham. "Penetration of a small drop into a capillary." Journal of Colloid and Interface Science 122, no. 1 (1988): 209-219. [100] Yue, Pengtao, and Yuriko Renardy. "Spontaneous penetration of a non-wetting drop into an exposed pore." Physics of fluids 25, no. 5 (2013): 052104. [101] Marmur, Abraham. "The radial capillary." Journal of colloid and interface science 124, no. 1 (1988): 301-308. 90 [102] Marmur, Abraham. "Drop penetration into a thin porous medium." Journal of colloid and interface science 123, no. 1 (1988): 161-169. [103] Choi, Hyunho, Lian Ma, and Hong Liang. "A kinetic study of the spontaneous penetration of a water drop into a hydrophobic pore." Surface Topography: Metrology and Properties 5, no. 1 (2017): 014003. [104] Domenico, Patrick A., and Franklin W. Schwartz. Physical and chemical hydrogeology. Vol. 506. New York: Wiley, 1998. [105] Anderson, William. "Wettability literature survey-part 2: Wettability measurement." Journal of petroleum technology 38, no. 11 (1986): 1-246. [106] Lamanna, Grazia, Simona Tonini, Gianpietro Elvio Cossali, and Bernhard Weigand, eds. Droplet Interactions and Spray Processes. Springer International Publishing, 2020. [107] Tredenick, Eloise C., Troy W. Farrell, W. Alison Forster, and Steven TP Psaltis. "Nonlinear porous diffusion modeling of hydrophilic Ionic agrochemicals in astomatous plant cuticle aqueous pores: a mechanistic approach." Frontiers in plant science 8 (2017): 746. [108] Matyka, Maciej. "Pushing a soft body droplet through porous medium." arXiv preprint arXiv:2004.14656 (2020). [109] Xie, Yanbo, et al. "High-efficiency ballistic electrostatic generator using microdroplets." Nature communications 5.1 (2014): 1-5. [110] Malik, R. S., C. H. Laroussi, and L. W. De Backer. "Experimental investigation of the penetration coefficient in capillary tubes." Soil Science 127, no. 4 (1979): 211-218. [111] Fisher, Leonard R., and Prosper D. Lark. "An experimental study of the Washburn equation for liquid flow in very fine capillaries." Journal of Colloid and Interface Science 69, no. 3 (1979): 486-492. [112] Langhaar, Henry L. "Steady flow in the transition length of a straight tube." J. appl. Mech. 9 (1942). [113] Davis, Stephen H. "On the motion of a fluid-fluid interface along a solid surface." Journal of Fluid Mechanics 65, no. 1 (1974): 71-95. [114] Ramé, Enrique, and Stephen Garoff. "On identifying the appropriate boundary conditions at a moving contact line: an experimental investigation." Journal of Fluid Mechanics 230 (1991): 97-116. [115] Chen, Jing-Den. "Experiments on a spreading drop and its contact angle on a solid." Journal of colloid and interface science 122, no. 1 (1988): 60-72. 91 [116] Siegel, Robert. "Transient capillary rise in reduced and zero-gravity fields." (1961): 165-170. 92 APPENDIX: USER DEFINED FUNCTION CODE User Defined Function code The dynamics of droplet impact and spreading on solid surfaces play a critical role in various engineering and scientific applications, including inkjet printing, spray cooling, and surface coating. Accurate simulation of such multiphase flow phenomena requires detailed tracking of the liquid-gas interface and quantification of key physical parameters such as droplet height, contact line radius, and spreading speed. To achieve this, a User Defined Function (UDF) was developed and implemented in ANSYS Fluent. The UDF consists of two main components: an initialization routine that defines the initial droplet shape and location using Volume of Fluid (VOF) values, and a post-processing routine that computes and records the droplet’s dynamic parameters at the end of each time step. The simulation captures the stick-slip behavior at the contact line and allows for detailed analysis of droplet evolution over time. The outputs from this UDF provide insight into the interface dynamics, enabling validation against experimental data and aiding in the design of optimized surface interactions. The computational methods used to simulate droplet dynamics, and the behavior of Stick-slip regimes. User defined function (UDF) is used in the solver to capture the interface between the liquid and the gas as following: 1. DEFINE_INIT(patching, d)  Purpose: Initializes the droplet by setting the Volume of Fluid (VOF) values in the cells to simulate a droplet positioned at a certain height (y = 0.0016 m) and with a given radius. 93 2. Key logic:  If a cell centroid lies within a sphere of radius 0.0016 m centered at y = 0.0016 m, the liquid phase VOF is set to 1 (droplet region), otherwise it is set to 0 (gas region).  Velocity C_V is set to zero initially, simulating impact from rest. 3. DEFINE_EXECUTE_AT_END(execute_at_end)  Purpose: Post-processing function that runs at the end of each time step to measure and log dynamic quantities:  Maximum height of the droplet interface.  Maximum radius of the contact line at the bottom wall (wall ID 9).  Minimum depth and corresponding velocity at the “hole wall” (wall ID 8). 4. Key features:  Uses C_VOF to track the interface.  PRF_GRHIGH1 and PRF_GRLOW1 ensure global maximum/minimum values across parallel processes.  Outputs data to a file contact.txt in the format:  time, height, radius, depth, speed. The code hase been written as the following: #include "udf.h" /* #define MAX(a,b) (((a)>(b))?(a):(b)) #define MIN(a,b) (((a)>(b))?(b):(a)) */ DEFINE_INIT(patching, d) { cell_t c; Thread* t; Thread* ts; Thread* tp; 94 real xc[ND_ND], dist2; #if !RP_HOST /* loop over all cell threads in the domain */ thread_loop_c(t, d) { ts = THREAD_SUB_THREAD(t, 1); tp = THREAD_SUB_THREAD(t, 0); /* loop over all cells */ begin_c_loop(c, t) { C_CENTROID(xc, c, t); dist2 = xc[0] * xc[0] + (xc[1] - 0.0016) * (xc[1] - 0.0016) + xc[2] * xc[2]; if (dist2 < 0.0016 * 0.0016) { C_VOF(c, ts) = 1.000000; C_VOF(c, tp) = 0.000000; C_V(c, t) = 0.0;// impact speed for droplet from H=5cm } else { } C_VOF(c, ts) = 0.000000; C_VOF(c, tp) = 1.000000; } end_c_loop(c, t) } Message("Patched successfully using UDF!\n"); #endif } DEFINE_EXECUTE_AT_END(execute_at_end) { Domain* d; Thread* t; Thread* tp; Thread* ts; Thread* tf; cell_t c; face_t f; d = Get_Domain(1); /* mixture domain if multiphase */ real Height, Radius, xc[ND_ND], xf[ND_ND]; real Depth, Speed; 95 /* loop over all cell threads in the domain */ Height = -10.; Radius = 0.; Depth = 0.; Speed = 0; #if !RP_HOST /* height */ thread_loop_c(t, d) { tp = THREAD_SUB_THREAD(t, 0); // loop over all cells begin_c_loop(c, t) { C_CENTROID(xc, c, t); if (C_VOF(c, tp) < 0.5) { Height = MAX(Height, xc[1]); } } end_c_loop(c, t) } /* contact line*/ tf = Lookup_Thread(d, 9); // 9 is wall ID of bottom wall tp = THREAD_SUB_THREAD(THREAD_T0(tf), 0); begin_f_loop(f, tf) { if PRINCIPAL_FACE_P(f, tf) { F_CENTROID(xf, f, tf); if (C_VOF(F_C0(f, tf), tp) < 0.5) { Radius = MAX(Radius, xf[0]); } else { } } } end_f_loop(f, tf) 96 /* wall_hole - ID 8*/ tf = Lookup_Thread(d, 8); // 8 is wall ID of hale wall ts = THREAD_SUB_THREAD(THREAD_T0(tf), 1); begin_f_loop(f, tf) { if PRINCIPAL_FACE_P(f, tf) { F_CENTROID(xf, f, tf); if (C_VOF(F_C0(f, tf), tp) < 0.5) { if (Depth > xf[1]) { Depth = xf[1]; Speed = C_V(F_C0(f, tf), ts); } } else { } } } end_f_loop(f, tf) Radius = PRF_GRHIGH1(Radius); Height = PRF_GRHIGH1(Height); Depth = PRF_GRLOW1(Depth); Speed = PRF_GRLOW1(Speed); #endif node_to_host_real(&Radius, 1); node_to_host_real(&Height, 1); node_to_host_real(&Depth, 1); node_to_host_real(&Speed, 1); #if !RP_NODE FILE* fp = NULL; if ((fp = fopen("contact.txt", "a")) != NULL) { fprintf(fp, "%f,\t%f,\t%f,\t%f,\t%f\n", CURRENT_TIME, Height, Radius, fabs(Depth), fabs(Speed)); } fclose(fp); #endif } 97