[7] M. Vahab, M. Y. Hussaini and M. Sussman, A moment-of-fluid method for multiphase systems with phase change, in prep.
[6] M. Vahab, C. Pei, M. Y. Hussaini, M Sussman and Y. Lian, An adaptive coupled level set and moment-of-fluid method for simulating droplet impact and solidification on solid surfaces with application to aircraft icing, 54th AIAA Aerospace Sciences Meeting, 2016, San Diego, CA, 2016 [Link][Preprint]
[5] M. Vahab, G. Miller, A front-tracking shock-capturing method for two gases, Communication in Applied Mathematics and Computational Science, 2015, [Link][Preprint]
[4] M. Vahab, A front-tracking shock-capturing method for two fluids, PhD Dissertation, UC Davis, 2014, [Link]
[3] M. Vahab, Design of a high-order front tracking method in 2D, Master Thesis, UC Davis, 2010, [Link]
[2] M. Vahab, Relative velocities of particles suspended in stochastic Kolmogorov turbulence, Master Thesis, Chalmers University, 2008, [Link]
[1] M. Vahab, M. Ebrahmi, Design and Implementation of a Control and Monitoring Environment Based on IEC Standard for Industrial Controller Structures, KNTU, 2006

Research topics

Conservative front-tracking and shock-capturing method

Shock capturing methods are very effective in modeling gas dynamics. However, when two or more materials with different properties are involved in such systems, common shock capturing methods show significant error in interfacial regions due to (1) diffusive solution around contact discontinuities (2) non-conservative modification to resolve material discontinuities. We have developed a numerical method for hyperbolic systems of conservation laws that captures shocks and tracks material interfaces, while keeps numerical conservation to machine precision. The method designed for Cartesian grids, and resolve the small cut-cell issue with an algebraic update equation that does not impose any extra limit on CFL condition than the high order finite volume method in the bulk regions. This method shows second order convergence for numerical tests of two-phase gas dynamics using Euler equations.

Figure: Richtmyer–Meshkov instability, pressure profile. The gases on left and right are sulfur hexafluoride (SF6)  and air respectively. See [4,5] for details.
Figure: Shock-bubble interaction, density profile. The gases inside and outside the bubble helium and air respectively. See [4,5] for details.

Moment-of-fluid method

The volume-of-fluid(VOF) methods are very simple and widely applied to multiphase problems. However, the interface definition is locally undefined. The interface reconstruction methods such as the least squares volume-of-fluid interface reconstruction algorithm (LVIRA), the efficient least squares volume-of-fluid interface reconstruction algorithm (ELVIRA) works efficiently for reconstructing a smooth interface for two-phase systems, yet, smooth out sharp edges and are less applied to multiphase systems. The coupled level set and volume-of-fluid (CLSVOF) methods are using the level set information to come up with better interface reconstruction, however, the multimaterial regions would be difficult to resolve using these methods as well. In an attempt to decrease the ambiguity of the interface reconstruction, one can include the moment-of-fluid (material first moment in the computational cell) as another tracked integrated variable to the geometry data. Therefore, material reconstruction at multimaterial regions is more accurate and geometrical properties are better preserved. One specific issue for MOF methods is that the reconstruction is local to the computational cell, and does not consider the geometrical information from the neighbor cells. So, the reconstruction is less suited for problems that inherently have strong effects from interface configuration in the neighbor cells, such as surface tension forces. We are introducing effective modification to the MOF method by considering neighbor cell geometrical information. The preliminary results has shown that it has better performance than unmodified MOF method for problems with strong surface tension forces.

Figure: (a) Physical domain: a part of the problem domain containing a triple point. Surface tension forces are at equilibrium determining the contact angles. (b) Discretized domain: crosses show cell centers, circles represent cell centroids and horizontal and vertical rectangles depict the MAC grid points in each cell.

Phase change: Solidification

The solidification of a pure material in a macroscopic model can generally be described as two different cases of normal freezing, and freezing with supercooling(undercooling). The normal freezing starts with liquid cooling, and nucleation happens at the freezing/melting temperature. The isothermal main freezing continues until the bulk material changes to solid phase, continuing with solid cooling. When the solidification process includes supercooling, the nucleation does not happen at the melting temperature due to the causes such as energy barrier of crystallization. The nucleation occurs at a temperature less than melting temperature, followed by the recalescence stage in which the temperature of the material jumps back to the melting temperature, and the isothermal main solidification takes place. Similarly, the solid cooling follows when all the bulk material is solidified.

Figure: Evolution of local temperature during the cooling and freezing procedures of a pure material. (a) Normal freezing (b) Freezing with supercooling. Stages are (1) Liquid cooling, (2) Liquid supercooling/undercooling, (3) Nucleation, (4) Recalescence, (5) Main freezing and (6) Solid cooling.

The liquid/solid cooling and main freezing stages are studied and understood extensively, and usually modeled macroscopically. The supercooling, nucleation and recalescence are more debated and discussed recently. One reason is the thermodynamically metastable nature of the the supercooled material. A wide variety of local and non-local causes may initiate nucleation and push the material from supercooled stage to recalescence. Common discussed variables are the critical activation energy for homogeneous nucleation, critical nucleus size, degree of supercooling, proximity to other nucleation sites. My goals is to have an approach that can capture the normal freezing physically, and extend the model for supercooled freezing while considering multistage freezing.

Figure: Cross section profile for the impact and freezing of a water droplet on a flat ice surface for three different Stefan number, which is the ratio of sensible heat to latent heat (a-c) Stefan number=0.075 (high value -> slow freezing process): Drop impact to surface and spread and rebound with little amount of solidification observed. (d-f) Stefan number=0.75 (medium value -> faster freezing process): After drop impact, solidification happens during the spread and rebound stage, can be seen as different cone slope in the final frozen droplet. (g-i) Stefan number=0.75 (small value -> fastest freezing process): After drop impact, the freezing process is do fast that the drop is frozen almost completely during the spread and partial rebound stages, making a non-conical shape. See [6] for details.

Kolmogorov turbulence and suspended particles

The direct numerical simulation of turbulence is a computationally expensive task. Also, if someone wants perform a statistical analysis of some aspects of turbulence, the task becomes costly in many degrees. Therefore, it is convenient to have a model to have enough similarity to turbulent to test a hypothesis, before spending time and resources on an actual direct numerical simulation of turbulent flows. A widely accepted description of turbulence is based on the Kolmogorov's theory, that predict the scaling law for relative velocity and turn over time of eddies at inertial range of flow. One can create a stochastic model based on the Kraichnan’s ensemble in space and modifying the harmonics of the stream function using the Ornstein-Uhlenbeck process in time to mimic the Kolmogorov turbulence. Using this approach, one can vary the spatial and temporal scaling of the flow to study its effects on the desired variables.  Also, one can modify this algorithm to evaluate only a handful of geometrically spaced harmonics in the desired wave number range, instead of all of them, to make its evaluation faster and more efficient while keeping the characteristics similar to the Kolmogorov turbulence. I have implemented the above mentioned stochastic flow, and verified its spatial and temporal scalings. Also, I have studied the distribution of zero-velocity and zero-acceleration points in such flow and came up with the corresponding analytical and computed fractal dimensions. This method can be also used for as a fast and easy test to check for geometrical properties of a turbulent flow field.

Figure: Fractal dimension measure for (a)zero-velocity points, and (b) zero-acceleration points of velocity field in one spatial dimension. The red, blue and black points represent zeroth, first, and second fractal dimensions corresponding to the box-counting, information, and correlation fractal dimensions respectively. It is shown that the fractal dimension of the zero-velocity points is 2/3, while they are distributed in one spatial dimensional. This verifies that the set of zero-velocity points has a fine-scale structure. In the case of zero-acceleration points, the simulation results confirm that the spatial distribution of zero-acceleration points is random, since the fractal dimension of their set is 1, equal to their domain’s physical dimension.

Furthermore, I have studied the behavior of particles suspended in turbulent flow. I identified a critical length scale for which the particle behavior can not be solely be characterized by the flow field around it. Also, I have shown that existence of large eddies may decrease the relative velocity of colliding particles. It is a counter intuitive results that may shed some light on the process of rain formation in the clouds.

Figure: Comparison of relative velocity of particles (blue) and relative velocity of fluid flow (red) at different separations. (a) L=10, (b) L=100, (c) L=1000. Extension of inertial range causes the decreased relative velocity of the particles at smaller distance (colliding particle), which may lead to particle conglomeration.