Theses supervised by Prof. Dr. Metin Muradoğlu
11 theses · Koç University
Computational modeling of airway closure
Capillary instability of a two-layer liquid film lining a rigid tube is studied computationally as a model for liquid plug formation and closure of human airways. The two-layer liquid consists of a serous layer, also called periciliary liquid layer (PCL), at the inner side and a mucus layer at the outer side. Together, they form the airway surface liquid (ASL) lining the airway wall and surrounding an air core. Liquid plug formation occurs due to the Plateau-Rayleigh instability when the liquid film thickness exceeds a critical value. Numerical simulations are performed for the entire closure process including the pre- and post-coalescence phases. The mechanical stresses and their gradients on the airway wall are investigated for physiologically-relevant ranges of the mucus-to-serous thickness ratio, the viscosity ratio, the air-mucus and serous-mucus surface tensions encompassing healthy and pathological conditions of a typical adult human lung. The growth rate of the two-layer model is found to be higher in comparison with a one-layer equivalent configuration. This leads to a much sooner closure in the two-layer model than in the corresponding one-layer model. Moreover, it is found that the serous layer generally provides an effective protection to the pulmonary epithelium against the high shear stress excursions and their gradients. Furthermore, a linear stability analysis using lubrication approximations is also performed, and the results are found to be in good qualitative agreement with the simulations. Airway mucus is a highly complex fluid, and displays a range of non-Newtonian characteristics. In order to take effects of elastoviscoplasticity and shear-thinning of the pulmonary mucus into account, the airway closure problem is studied by modeling the mucus as a single layer, i.e., PCL and ASL are considered as a homogeneous one-layer elastoviscoplastic (EVP) liquid, using the Saramito-HB model [Saramito (Journal of Non-Newtonian Fluid Mechanics (2009) 158(1-3), 154-161)]. Firstly healthy, asthma, cystic fibrosis (CF) and chronic obstructive pulmonary disease (COPD) mucus are represented by following a parameter fitting algorithm (similar to [Fraggedakis et al. (Soft Matter (2016) 12(24), 5378-5401)]) by using the experimental data of [Patarin et al. (Scientific Reports (2020) 10(1), 1-10)]. The obtained parameters show that CF mucus has the highest elastic modulus and the yield stress is the highest in COPD case. This generally confirms that the viscoelastic/viscoplastic features of the mucus increases in pulmonary diseases. Afterwards, extensive numerical simulations are performed to see the effects of these diseases on the airway closure. To that end, the diseased mucus parameters are used, and the surface tension and the initial non-dimensonal thickness are varied in a physiologically meaningful range. Then, the healthy case is chosen as the baseline to perform a parametric study. Here, the elastic modulus, the yield stress and the power-law index of the Saramito-HB model are varied to assess the individual effects of these parameters. The results indicate that increasing yield stress and elastic modulus may abate the airway closure. However, if the initial thickness of the liquid layer and the surface tension increase, which is a typical case in these obstructive pulmonary diseases, then the closure starts to occur for the high viscoelastic/viscoplastic cases as well.
Fluid-structure interaction modelling of plug propagation in compliant airways
The lung is the largest organ that consists of many generations of liquid-lined compliant tubes. Cross-section of airways decreases at each generation down to microscopic scale at the terminal bronchioles that connect to alveoli. The liquid lining may undergo Plateau-Rayleigh instability creating a liquid plug that obstructs gas exchange for all connected distal parts of the lung. The closure occurs when the liquid film thickness exceeds a critical threshold value. Liquid film thickness is about 2-4% of airway diameter in healthy conditions but may increase up to 20% in pathological conditions. The liquid plug is ruptured violently to reopen the airways during the inspiration. Both airway closure and reopening create large mechanical stresses on the airway wall and may fatally damage the epithelial cells. The plug propagation and rupture have been studied both experimentally and computationally, but the airway wall has been usually assumed to be rigid in the previous computational models. However, especially the small airways where the plug formation is more likely to occur are highly flexible and compliant tubes that can completely collapse during the closure and reopening events. Moreover, the crackles sound that is detected by a stethoscope and routinely used for diagnostics examinations by physicians are believed to be created during the plug rupture and subsequently modified by strong fluid-structure interactions. It is therefore of fundamental importance to take the airway wall elasticity into account and examine fluid-structure interactions. In the present thesis, a sharp interface immersed boundary method is developed and combined with the front-tracking method to investigate the effects of wall elasticity on liquid plug propagation and rupture in lower airways. The immersed boundary method is first applied to simulate a single-phase flow past a solid sphere in a circular pipe, and the results are found to be in good agreement with the results obtained by Ansys-Fluent. It is also demonstrated that the immersed boundary method preserves the second-order spatial accuracy of the finite volume method used to solve the flow equations. Then extensive simulations are performed to investigate the fluid-structure interactions in plug propagation and rupture in a compliant tube representing the 9th generation of an adult lung. The tube is assumed to be infinitely long, thin-walled, and impermeable. The wall elasticity is modeled using linear elastic theory. The results show that there exist strong fluid-structure interactions during the plug propagation, plug rupture, and post-rupture phases. The present study lays the foundation for future work in which a full three-dimensional airway model will be developed, and acoustic signals generated during the airway closure and reopening will be examined.
Effects of kinematic hardening of mucus polymers in an airway closure model
Formation of a liquid plug inside an airway of a human lung, known as airway closure, is numerically studied by considering the elastoviscoplastic (EVP) properties of the pulmonary mucus covering the airway walls. The airway is modeled as a rigid tube with single layer of liquid coating its wall. As the model for the liquid layer, the Saramito-Herschel-Bulkley (Saramito-HB) model is coupled with the Isotropic Kinematic Hardening model (Saramito-HB-IKH). Aim of the study is to examine the kinematic hardening features of mucus and evaluate its effects on the closure process. The rheological model is fitted to the experimental mucus data under healthy, asthma, chronic obstructive pulmonary disease (COPD) and cystic fibrosis (CF) conditions. Comparing the fitting results for both models with respect to the experimental data, the Saramito-HB-IKH model is shown to be superior to the Saramito-HB model, especially in terms of capturing the loss moduli under low strain amplitude region of a constant frequency sweep test. Using the mucus parameters acquired through the fitting process, airway closure is computationally investigated for a range of initial undisturbed liquid layer thickness and the Laplace number. The evolution of the effective yield stress, which is a dynamically evolving yield stress introduced in the kinematic hardening model, is examined throughout the closure process. The parameter is observed to accumulate more rapidly for the CF case compared to the other cases. Yielded/unyielded regions of the liquid layer and the stresses on the airway wall are similarly examined. Yielding occurs immediately in the Saramito-HB-IKH model, whereas it commences near closure for the Saramito-HB model. Comparing the wall stresses between two models, kinematic hardening is seen to have notable effect on the closure time, especially for the CF case, with the effect being more pronounced at low Laplace numbers and low initial liquid layer thicknesses. The Saramito-HB-IKH model also captures closure behavior for high elastic modulus cases of COPD and CF under conditions which Saramito-HB model is unable to. The sensitivity of the results to the initial condition of the effective yield stress is evaluated, and it is determined that closure is delayed as the initial value increases. The study additionally reports the standalone effects of the rheological properties of mucus on the wall stresses considering their physiological values as baseline. The elastic modulus is found to influence the closure time and the magnitudes of the stresses and their gradients. The polymeric viscosity particularly affects the relaxation of stresses after closure. The yield stress is revealed to have an important factor in determining the closure time. However, its effects are understood to be case dependent and contingent on the initial liquid layer thickness and Laplace number.
Konvektif ortamlarda damlacik buharlaşmasinin analitik ve sayisal olarak incelenmesi
The significance of Stefan flow lies in its impact on the behaviour of evaporating droplets in convective environments, particularly in scenarios involving high-temperature evaporation of fuel droplets relevant to spray combustion. We examine this phenomenon from both numerical and analytical perspectives. An axisymmetric multiphase front-tracking method is utilized to simulate the evaporation of a deformable droplet in a convective environment, and the local Sherwood number is computed on the droplet for the Reynolds number up to 200. Simulations are also performed for a solid sphere using a sharp-interface immersed boundary method. A novel analytical model is derived based on the solution of the laminar boundary layer on a sphere subjected to local evaporation velocity. The results are compared with the existing evaporation models widely used in spray combustion simulations as well as the new model. We demonstrate that at high evaporation rates, the boundary layer-based models can accurately predict the local Sherwood number only until the separation point. However, the generalization of these models for the wake region result is a considerable error. It is argued that for this region a new correction factor should be derived.
Parallel direct numerical simulation of 3D interface-resolved droplet evaporation
Multi-phase flows and phase transitions are widely observed phenomena in natu- ral processes and engineering applications such as atmospheric precipitation, bubbles forming during water boiling, and spray combustion. To study the physics of such phenomena and to broaden their application areas, simulations of multiphase flows with phase change are crucial. Therefore, in this thesis, a parallel front-tracking/finite-difference method is developed for interface-resolved fully 3D simulations of droplet evaporation. In this method, the mass, momentum, and energy conservation equations are solved on a 3D-Cartesian, uniform, fixed, staggered Eulerian grid while a Lagrangian grid is employed to track the interface between the liquid and vapor phases. The method is first validated for the classical d2-law. Then the method is applied to simulate the evaporation of a 3D Leidenfrost droplet.
Computational modeling and analysis of viscoelastic multiphase flows
Interface-resolved direct numerical simulations of multiphase viscoelastic fluid flows have been performed by using a high fidelity finite-difference/front-tracking method to explore some of the most complex dynamics of these intricate flow regimes. Surfactant contamination is also considered. The flow equations are solved fully coupled with the viscoelastic model equations and the equations governing interfacial and bulk surfactant concentrations. The latter coupling is accomplished by a nonlinear equation of state that relates the surface tension to the surfactant concentration. Simulations are first performed to study the phenomenon of negative wakes generated behind a pair of bubbles rising in an otherwise quiescent viscoelastic fluid under buoyancy. It is found that, in addition to the formation of long tails that tilt towards each other, the lateral symmetry is also broken in the two-bubble system. The rise velocity of the parallel bubbles remains slightly smaller than that of the single bubble due to bubble-bubble interactions. Symmetry observed in the case of parallel bubble pairs is broken when bubbles are initialized by an initial offset in vertical direction. Next, lateral migration of a deformable bubble in a viscoelastic pressure driven channel flow is studied. It is found that the forces induced by fluid inertia push a bubble towards wall while elasticity and deformability pull it towards channel center. This interplay between fluid inertia, elasticity, and bubble deformability determines the orientation of bubble migration in the channel and its final equilibrium position. It is shown that secondary flow velocity induced by a bubble migrating towards wall is an order of magnitude higher than the one induced by a solid particle under similar flow conditions. This behavior is attributed to the slip condition on the bubble surface. Then, a novel mechanism of achieving an elasto-inertial turbulent (EIT) state by injecting bubbles into a pressure-driven viscoelastic channel flow under the conditions for which a single-phase flow remains laminar is demonstrated even for the Reynolds number as low as Re = 10 and a Weissenberg number as low as Wi = 1. It is found that curved streamlines across the bubble interfaces provide the necessary condition to trigger an elastic instability that leads to a fully chaotic motion even at very low Reynolds and Weissenberg numbers for which the single-phase flow remains completely laminar. The energy spectra shows a scaling of −2 for this multiphase EIT regime. It is also found that bubbles move towards the channel centerline and form a string-shaped alignment pattern in the core region at lower values of Re = 10 and Wi = 1. Unlike solid particles, increasing shear-thinning effect tends to break up alignment of bubbles. Finally, interface-resolved direct numerical simulations are performed to investigate drag reduction by polymer additives in a bubbly turbulent channel flow in the presence of surfactant contamination for a fixed Reynolds number of Re = 5600. Simulations are carried out for the high-drag-reduction (HDR) and maximum-drag-reduction (MDR) regimes. It is shown that, once flow is made viscoelastic by polymer additives, an HDR regime is achieved with a drag reduction of 49%. It is found that Reynolds stress is reduced in this HDR regime but still remains higher than polymer stress. Injection of contaminated bubbles reduces Reynolds stress even further and drag is not increased. Reynolds stress in single-phase flow of maximum MDR regime reduces below the value of polymer stress with a drag reduction of 68% and flow is dominated by elastic effects. Injection of bubbles in this MDR regime produces non-monotonic effect. A small concentration of contaminated bubbles increases Reynolds stress, turbulence production by mean flow and turbulent kinetic energy whereas a higher concentration of bubbles in the MDR regime reduces all these quantities. However, these effects remain restricted to core region of the channel and do not increase drag. It is found that turbulent kinetic energy spectrum shows a scaling of −2 in both of these multiphase HDR and MDR regimes and this scaling is not affected by a change in bubble concentration.
Bir yakıt damlacığının buharlaşması ve yanmasının doğrudan sayısal simülasyonu
A front-tracking method is developed for the Direct Numerical Simulation (DNS) of droplet evaporation and combustion in a liquid-gas multiphase system. One field formulation is used to solve the flow, energy and species equations with suitable jump conditions. Both phases are assumed to be incompressible; however, the divergence-free velocity field condition is modified to account for the phase change at the interface. Both temperature and species gradient driven phase change processes/models are simulated in 2D planar configuration. For the species gradient driven phase change process, the Clausius-Clapeyron equilibrium relation is used to find the vapor mass fraction and subsequently the evaporation mass flux at the interface. Extensive validation studies are performed using the benchmark cases: The Stefan and the sucking interface problems, d2-law and the wet bulb temperature comparison with the psychrometric chart values. The models are then applied to simulate the evaporation of a single and two droplets systems that move, evaporate, interact and undergo significant deformation in a gravitational field. The implementations have been demonstrated to be grid convergent and the global mass conservation is satisfied for all the studied cases. The temperature gradient based phase change model is then incorporated into the axisymmetric multiphase solver and comprehensive validation studies are performed. This is then further extended to model the burning process following the evaporation as a first step towards the development of a computational framework for the direct numerical simulations of spray combustion. We used single step as well as reduced chemical kinetic mechanism, constant thermodynamic properties, unity Lewis number for all the species and an ideal gas behaviour for a n-heptane droplet combustion; the chemical kinetics being handled by the CHEMKIN. An operator-splitting approach is used to advance temperature and species mass fractions in time. The numerical results of the droplet burning rate, flame temperature, flame standoff ratio and ignition delay times show good agreement with the experimental and previous numerical studies.
Yüzey aktif maddesi ve viskoelastisitenin damacık dinamiğine üzerindeki kombine etkileri
Combined effects of surfactant and viscoelasticity on droplet dynamics are studied in an axisymmetric tube using a front-tracking method. Both rising and sedimenting droplet cases are considered. Viscoelasticity is contained in the droplet fluid while surfactant is contained in the continuous phase. The FENE-CR model is used to account for droplet liquid viscoelasticity. The incompressible Navier-Stokes equations are solved fully coupled with the surfactant evolution and the viscoelastic model equations in the entire computational domain. The surface tension is related to the interfacial surfactant concentration by a non-linear equation of state. Extensive simulations are performed to investigate the effects of both viscoelasticity and surfactant on droplet dynamics in nearly spherical, dimpled ellipsoidal cap and skirted regimes. It is found that the surfactant-induced Marangoni stresses generally counteract the viscous shear stresses to reduce mobility of the interface and thus decrease the terminal velocity of the droplet. This effect is most pronounced in the nearly spherical regime in which the interface is almost fully immobilized making the droplet behave like a solid sphere. In the dimpled ellipsoidal cap and skirted regimes, the effects of surfactant diminish significantly since the surfactant adsorbed at the leading edge is quickly convected toward the rear edge rendering the most part of the interface clean. The viscoelastic stresses are found to be nearly uniform in the nearly spherical regime and thus have negligible influence on drop dynamics. However, the viscoelastic stresses concentrate around the rear stagnation point and induce an indentation there especially in the dimpled ellipsoidal cap and skirted regimes. This effect is significantly amplified in the presence of the surfactant leading to creation of a hole at the drop centerline in extreme cases.
Partikül-çözümlenmiş çok-fazlı akışların simülasyonu için geliştirilen adaptif örgü artırım metodu
Increasing computational resources allow the fast development of numerical methods to simulate complex physical systems. However, numerical analysis of multiphase flows is still a challenging task due to wide disparity in length and time scales. Direct numerical simulation (DNS) methods, such as front-tracking method, provide valuable insights into multiphase flows by resolving all scales. But, in most of the DNS applications the quality of the grid resolution plays a crucial role in development of computational models of such flow systems due to presence of sharp interface between the phases. However the interface usually occupies a small portion of the overall domain, hence applying the resolution around the interface in the entire computational region is not an efficient way. Block-structured adaptive mesh refinement (AMR) method, which was first developed by Berger and Oliger [Berger and Oliger, 1984], offers a local and adaptive refinement of the grid. A properly-nested hierarchy of refinement levels increases the resolution around the region where the predetermined error criterion is exceeded. Almgren et al. [Almgren et al., 1998] combined block-structured AMR with a second-order projection method for incompressible Navier-Stokes equations and increased the efficiency by recursive time refinement algorithm. However, time refinement results in cumbersome synchronization operations to match coarse and fine levels after each time step. On the other hand, their projection method approximately satisfies the divergence-free constraint, hence, maintaining global conservation requires additional algorithmic complexity. Vanella et al. [Vanella et al., 2010] applied the structured AMR method to fluid-solid interaction problem by using a staggered grid arrangement and they ignored time refinement. Even though their multilevel multigrid solver simplified the synchronization step and satisfied divergence-free constraint exactly, using the same time step size in the finest level in all other levels reduced the computational efficiency. The present study combines the subcycled block-structured AMR method with the three-dimensional finite-volume/front-tracking method developed by Unverdi and Tryggvason [Unverdi and Tryggvason, 1992]. The purpose is to accurately resolve multiphase flow by satisfying the divergence-free constraint exactly and to gain the advantage of efficient AMR algorithm with time refinement at the same time. The algorithm presented here also avoids the complexity of synchronization step by using a fully-staggered grid arrangement first proposed by Harlow and Welch [Harlow and Welch, 1965]. The validation is performed by solving two benchmark problems: (1) The Hagen-Poiseuille problem with variable density/viscostiy and (2) mutiphase flow with a stationary bubble. Although the results to benchmark problems are very promising, the algorithm should be extended to simulate moving and deforming bubbles in the future studies. Moreover, to asses the efficiency of the algorithm, the performance analysis should be conducted when the algorithm is completed.
Immersed boundary method for fluid-fluid-solid three phase flows
Multiphase flows can be observed in a broad range of industrial and natural processes, and their ubiquity puts them in a significant position in computational fluid dynamics (CFD). This importance also brings a requirement of understanding the underlying physics in detail and one way to achieve this goal is to do high fidelity simulations based on the first principles in a computational environment. These multiphase systems can contain fluid-fluid interactions, fluid-solid interactions and often both of them, therefore various methods have been developed and successfully applied to a wide range range of problems of practical importance in the last few decades. The multiphase flow systems are particularly of interest in micro-robotic applications such as a chemotherapy delivery operation to a cancerous area by a remote-controlled microswimmer and in biological flows such as the flow in human lungs. The ultimate goal of this thesis is to develop a new high-fidelity numerical method for simulations of fluid-solid interactions and combine it with a front-tracking algorithm to simulate fluid-fluid-solid interactions in three phase systems. In the present study, an immersed boundary method (IBM) is developed to treat fluid-solid interfaces that are fully or partially immersed in a host fluid. Using the front-tracking framework, the solid-fluid interfaces are represented by a Lagrangian grid that consists of marker points located at the vertices of a 3D triangular surface mesh. The direct forcing method developed by Uhlmann [152] is used to impose the no-slip and no-penetration boundary conditions on the interface. The Newton-Euler equations are solved for the rigid body dynamics of the fully immersed undeformable solid particles in a Lagrangian frame and tightly coupled with the continuum flow equations solved on a fixed Eulerian grid. For this purpose, quantities calculated on the Lagrangian grid are smoothed onto the Eulerian grid with a distribution function and the same distribution function is also used to transfer quantities from the interface to the Eulerian grid in a conservative manner similar to the distribution/interpolation operations in the front-tracking method. The hybrid immersed-boundary/front-tracking (IBM/FTM) solver is then validated for several benchmark cases. First, it is applied to simulate a single stationary rigid sphere exposed to a uniform flow at different Reynolds numbers and the results are found to be in good agreement with the results in the literature. Then, sedimentation of a single sphere in a resting fluid is simulated for various Reynolds numbers and the results are compared with the experimental results of Mordant and Pinton [102]. The method is demonstrated to be convergent in terms of spatial and temporal errors. The numerical method is also validated for a wide range of Reynolds numbers studied experimentally by Cate et al. [27]. Finally, the method is used to simulate particulate flows involving many particles and preliminary results are obtained. The method is found to be robust and the results look qualitatively accurate but a more rigorous quantitative comparison must be done for the validation, which is a subject of a future study. The method is totally developed in the front-tracking framework and can be readily used to simulate three phase flows involving fluid-fluid and fluid-solid interfaces.
Kabarcıklı Türbülanslı Akışlarda Yüzey Aktif Maddesi ve Viskoelastisitenin Etkileri
Interface-resolved direct numerical simulations are performed to examine the sole and combined effects of soluble surfactant and viscoelasticity on the structure of a bubbly turbulent channel flow. The incompressible flow equations are solved fully coupled with the FENE-P viscoelastic model and the equations governing interfacial and bulk surfactant concentrations. The latter coupling is achieved through a non-linear equation of state which relates the surface tension to the surfactant concentration at the interface. The two-fluid Navier-Stokes equations are solved using a front-tracking method. First, the effects of soluble surfactant on the lateral migration of a bubble in a pressure-driven channel flow are examined. Extensive computations are performed to investigate the bubble dynamics for a wide range of parameters. It is found that surfactant dramatically changes the bubble dynamics. In the clean case, the bubble position depends on its deformability, characterized by the Eotvos () and the capillary () numbers. The spherical bubble moves towards the wall, while the deformable one migrates away from it. On the other hand, in the presence of the surfactant, even the spherical bubble moves away from the wall. It is also found that the contaminated bubble stays away from the wall for and while it migrates towards the wall for . Also, at high , the onset of path instability is observed for both the clean and the contaminated cases. However, adding surfactant to the system triggers the path instability earlier and amplifies the oscillations afterwards. Next, simulations are performed to examine the combined effects of soluble surfactant and viscoelasticity on the structure of a bubbly turbulent channel flow. It is found that, for the surfactant-free case, bubbles move toward the wall due to inertial lift force, resulting in formation of wall layers and a significant decrease in the flow rate. Conversely, a high-enough concentration of surfactant changes the direction of lateral migration of bubbles, i.e., the contaminated bubbles move toward the core region and spread out across the channel. When viscoelasticity is considered, viscoelastic stresses counteract the Marangoni stresses, promoting formation of bubbly wall-layers and consequently strong decrease in the flow rate. The formation of bubble wall-layers for combined case depends on the interplay of the inertial and elastic, and Marangoni forces. Further, effects of different types of surfactant are examined using the physical sorption kinetics. For Newtonian turbulent bubbly flow, effects of Triton X-100 and 1-Pentanol are examined. It is observed that the sorption kinetics highly affect the dynamics of bubbly flow. Minute amount of Triton X-100 is found to be sufficient to prevent the formation of bubble clusters. On the other hand, 100 times more of 1-Pentanol surfactant is not sufficient to prevent the formation of layers. For viscoelastic turbulent flow, it is found that, polymer drag reduction is completely lost for the surfactant-free case, and the addition of small amount of surfactant (Triton X-100) restores the polymer drag reduction for the viscoelastic turbulent bubbly flows. Lastly, the polymer drag reduction of turbulent poly-dispersed bubbly flow is examined in the presence of soluble surfactant. It is found that the drag increases for the clean viscoelastic poly-dispersed bubbly flow but it is lower than the corresponding Newtonian flow due the lateral migration of large bubbles to the center of channel. In the presence of surfactant, the drag reduction for viscoelastic poly-dispersed bubbly flow is revived.