We study the role of the capillary number, Ca and of the surface wettability on the dynamics of the interface between an invading and a defending phase in a porous medium by means of numerical simulations. We employ a hybrid phase field-immersed boundary approach to successfully model the contact line dynamics over the solid objects. Using a phase-field method which naturally incorporates dynamic wetting we eliminate the need for empirical contact line models to address contact line singularity. We map the two dominant modes governing the motion of the interface, namely, capillary fingering, and stable penetration, in the (Ca - theta) plane, with theta the static contact angle prescribed at the solid pores. Capillary fingering dominates at lower values of Ca and pores hydrophobic to the invading phase, while a stable penetration is observed on hydrophillic surfaces. We present new measurements and analyses, including curvature probability density functions and average curvature. We also show that the pressure needed for the invading phase to advance at constant flow rate decreases with the capillary number, and increases with the contact angle at the capillary numbers considered. The latter is due to a significant increase in the length of the interface in the case of capillary fingering. Finally, we show that it is possible to identify the different interfacial modes by measuring the penetration length and velocity during the medium filling.

The intricate nature of non-Newtonian fluid rheology has raised notable attention, particularly in gas–liquid systems, where the dispersed bubbles may generate shear forces and change the shear-dependent viscosity of the surrounding liquid. While the effective shear rate, γ̇eff=vb/Db, is commonly used to approximate the shear-thinning viscosity around spherical bubbles, deviations may arise for deformed bubbles present in real systems. This work combines laboratory experiments and numerical simulations to investigate the evolution of a single rising bubble in three different systems: water, glycerol/water solutions characterizing viscous-Newtonian systems, and carboxymethyl cellulose (CMC) aqueous solutions exhibiting shear-thinning. The experiment was performed with bubble sizes of 1–9mm using imaging techniques. The measured fluid rheology is modeled by the Carreau model, and used in 3D direct numerical simulations based on a diffuse interface approach. The shear-thinning behaviors are found to increase the bubble terminal velocity through two distinct mechanisms: reducing the apparent viscosity around the bubble and promoting the bubble deformation. The extent of the shear-thinning effect depends on the three dominating regimes under which different rheology parameters play a significant role. Finally, empirical models for bubble terminal velocity and drag coefficient are evaluated using two shear-thinning viscosity estimations, based on the effective shear rate and the average shear-thinning viscosity near the bubble interface. The good agreement between experimental and simulation results validates the proposed models.

The extensional dynamics of an elasto-viscoplastic (EVP) fluid is studied by means of numerical simulations modelling an experimental configuration. Specifically, we track the interface between the EVP material and the Newtonian medium using an algebraic volume of fluid method (MTHINC-VOF) and employ a fully Eulerian immersed boundary method (IBM) to model the motion of the piston responsible for the extension of the material. We investigate the role of different values of the yield stress, surface tension at the interface between the EVP material and the surrounding fluid, polymer viscosity ratio, and extension rates on the necking thickness of the material, extensional viscosity, and yielding of the material for two sets of parameter with low and high elasticity. The results of the simulations reveal that when the yield stress of the EVP material is much larger than the viscous stresses, the material undergoes an elastic deformation, regardless of the selected values of the extension rate, interfacial forces, and viscosity ratio. Moreover, by increasing the ratio of the polymeric viscosity to the total viscosity of the system, the EVP material produces stronger strain hardening and reaches the minimum resolvable width sooner. Specific and novel to our study, we show that interfacial forces cannot be ignored when the surface tension coefficient is such that a Capillary number based on the extensional rate is of order 1. For large values of the surface tension coefficient, the EVP material fails sooner, with a clear deviation from the exponential reduction in the neck thickness. Moreover, our results suggest that the role of the yield stress value on the dynamics of the material is more pronounced at lower elasticity.

A numerical scheme to simulate three-phase fluid flows with phase change is proposed. By combining the Cahn-Hilliard model for water-air interface, Allen-Cahn equation for ice and fluid and Navier-Stokes equation for momentum, we solve the evolution of the water-air interface and water-ice interface simultaneously, including the volume expansion associated with solidification and due to the density difference between water and ice. Unlike existing schemes assuming a divergence-free flow field, the proposed continuous formulation allows for density changes while ensuring mass conservation. A Poisson equation for the pressure field is derived from mass conservation with constant coefficients, which can efficiently be solved without any pre-conditioning. The results demonstrate that the volume expansion during the ice formation and the subsequent motion of the water-air interface are successfully captured. A parametric study is carried out to examine the dependence of the icing on different physical and numerical parameters. Computations with flow disturbance of different amplitudes demonstrate the robustness of the computational scheme and the uniqueness of the solution over the parameters considered.

Multi-phase fluid flows are ubiquitous in natural phenomena and different industrial applications such as in food industry, the medical sector, heat exchangers, power generation systems, to name a few.  Understanding the underlining  physics of multi-phase flows  proved to be a challenging task due to presence of sophisticated dynamics, including the evolution of the interface between any pair of phases, thermodynamics and possibility of phase change, interactions between the fluid phases and a solid phase, etc.  Together with theoretical studies and experiments performed on a variety of multi-phase flow problems, numerical simulations have been employed by many researchers to scrutinise different aspects of the problem. During the last decades, a great many studies have been conducted aiming to provide more accurate numerical frameworks for investigating multi-phase flow problems.

Among the various complicated aspects of a multi-phase flow, the present thesis is focused on few characteristics of it the understanding of which requires more considerations and demands improvements in the numerical frameworks. First, we elaborate on the different interface tracking approaches suit the study of different multi-phase flows. In particular, a Volume of Fluid method, a compressible formulation of a diffuse interface approach, a Cahn-Hilliard phase field method, and an Immersed Boundary method are employed to study wetting phenomemna and fluxes at the interface. We have initially investigated biological-relevant membranes, extensional dynamics of a Elasto-viscoplastic material, and droplet spreading over rough surfaces.  In the second part of the thesis, we propose novel numerical methods and setups to investigate the phase change problems in both nanoscale and mesoscale. In particular, we developed a novel numerical method for the solidification problem, a pressure control setup for studying boiling at nanoscale, and a pressure based algorithm for modelling the boiling and evaporation.

We present a fully Eulerian hybrid immersed-boundary/phase-field model to simulate wetting and contact line motion over any arbitrary geometry. The solid wall is described with a volume-penalisation ghost-cell immersed boundary whereas the interface between the two fluids by a diffuse-interface method. The contact line motion on the complex wall is prescribed via slip velocity in the momentum equation and static/dynamic contact angle condition for the order parameter of the Cahn-Hilliard model. This combination requires accurate computations of the normal and tangential gradients of the scalar order parameter and of the components of the velocity. However, the present algorithm requires the computation of averaging weights and other geometrical variables as a preprocessing step. Several validation tests are reported in the manuscript, together with 2D simulations of a droplet spreading over a sinusoidal wall with different contact angles and slip length and a spherical droplet spreading over a sphere, showing that the proposed algorithm is capable to deal with the three-phase contact line motion over any complex wall. The Eulerian feature of the algorithm facilitates the implementation and provides a straight-forward and potentially highly scalable parallelisation. The employed parallelisation of the underlying Navier-Stokes solver can be efficiently used for the multiphase part as well. The procedure proposed here can be directly employed to impose any types of boundary conditions (Neumann, Dirichlet and mixed) for any field variable evolving over a complex geometry, modelled with an immersed-boundary approach (for instance, modelling deformable biological membranes, red blood cells, solidification, evaporation and boiling, to name a few). 

We study the role of surface topology, surface chemistry, and wall superheat temperature on the onset of boiling, bubble nucleation and growth, and the possible formation of an insulating vapour film by means of a novel setup for large-scale MD simulations. To minimise the effects of the system size on the bubble growth and the formation of the vapour film, we perform simulations in a box larger than those previously considered. The effect of the system pressure on bubble nucleation and growth is isolated by imposing a constant force on a moving piston and mechanically controlling the pressure. The simulations reveal that the presence of a nanostructure determines the nucleation site and facilitates the energy transfer from the hot substrate to the water. The surface chemistry, on the other hand, governs the shape of the formed bubble. A hydrophilic surface accelerates the bubble nucleation, however, decelerates the bubble expansion, thus postponing the formation of the film of vapour. Hence, a hydrophilic surface provides better energy transfer from the hot wall to the water. By analysing the system energy, we show that irrespective of wall topology and chemistry, there is a wall temperature for which the amount of transferred energy from the wall is maximum.

We study the rheology of suspensions of capsules with a rigid nucleus at negligible and finite flow inertia by means of numerical simulations. The capsule membrane is modeled as a thin Neo-Hookean hyperelastic material and the nucleus as a rigid particle with radius equal to half the radius of the undeformed spherical capsules. The fluid and solid motion are coupled with an immersed boundary method, validated for both the deformable membrane and the rigid nucleus. We examine the effect of the Reynolds number, capillary number, and volume fraction on the macroscopic properties of the suspensions, comparing with the case of capsules without nuclei. To explain the rheological measurables, we examine the mean capsule deformation, the mean orientation with respect to the flow direction, and the stress budget. The results indicate that the relative viscosity decreases with the capillary number, i.e., increasing deformability, and increases with inertia. The presence of a nucleus always reduces the membrane deformation. Capsules align more in the flow direction at higher capillary numbers and at higher volume fractions, where we also see a significant portion of them oriented with their longer deformed axis in the spanwise direction. When increasing inertia, the alignment with the flow decreases while more capsules orient in the spanwise direction. The first normal stress difference increases with the capillary number and it is always less for the nucleated capsules. Finally, the relative viscosity and the first normal stress difference increase with the capsule volume fraction, an effect more pronounced for the first normal stress difference.

We have performed direct numerical simulation of the turbulent flow of a polymer solution in a square duct, with the FENE-P model used to simulate the presence of polymers. First, a simulation at a fixed moderate Reynolds number is performed and its results compared with those of a Newtonian fluid to understand the mechanism of drag reduction and how the secondary motion, typical of the turbulent flow in non-axisymmetric ducts, is affected by polymer additives. Our study shows that the Prandtl's secondary flow is modified by the polymers: the circulation of the streamwise main vortices increases and the location of the maximum vorticity moves towards the centre of the duct. In-plane fluctuations are reduced while the streamwise ones are enhanced in the centre of the duct and dumped in the corners due to a substantial modification of the quasi-streamwise vortices and the associated near-wall low- and high-speed streaks; these grow in size and depart from the walls, their streamwise coherence increasing. Finally, we investigated the effect of the parameters defining the viscoelastic behaviour of the flow and found that the Weissenberg number strongly influences the flow, with the cross-stream vortical structures growing in size and the in-plane velocity fluctuations reducing for increasing flow elasticity.We have performed direct numerical simulation of the turbulent flow of a polymer solution in a square duct, with the FENE-P model used to simulate the presence of polymers. First, a simulation at a fixed moderate Reynolds number is performed and its results compared with those of a Newtonian fluid to understand the mechanism of drag reduction and how the secondary motion, typical of the turbulent flow in non-axisymmetric ducts, is affected by polymer additives. Our study shows that the Prandtl's secondary flow is modified by the polymers: the circulation of the streamwise main vortices increases and the location of the maximum vorticity moves towards the centre of the duct. In-plane fluctuations are reduced while the streamwise ones are enhanced in the centre of the duct and dumped in the corners due to a substantial modification of the quasi-streamwise vortices and the associated near-wall low- and high-speed streaks; these grow in size and depart from the walls, their streamwise coherence increasing. Finally, we investigated the effect of the parameters defining the viscoelastic behaviour of the flow and found that the Weissenberg number strongly influences the flow, with the cross-stream vortical structures growing in size and the in-plane velocity fluctuations reducing for increasing flow elasticity.