Direct Numerical Simulation Method Research Articles

Numerical Study on Bubble Rising in Complex Channels Saturated with Liquid Using a Phase-Field Lattice-Boltzmann Method

Packed bed reactors have been widely applied in industrial production, such as for catalytic hydrogenation. Numerical simulations are essential for the design and scale-up of packed beds, especially direct numerical simulation (DNS) methods, such as the lattice-Boltzmann method (LBM), which are the focus of future researches. However, the large density difference between gas and liquid in packed beds often leads to numerical instability near phase interface when using LBM. In this paper, a lattice-Boltzmann (LB) model based on diffuse-interface phase-field is employed to simulate bubble rising in complex channels saturated with liquid, while the numerical problems caused by large liquid-to-gas density ratio are solved. Among them, the channel boundaries are constructed with regularly arranged circles and semicircles, and the bubbles pass through the channels accompanied by deformation, breakup, and coalescence behaviors. The phase-field LB model is found to exhibit good numerical stability and accuracy in handing the problem of the bubbles rising through the high-density liquid. The effects of channel structures, gas-liquid physical properties, and operating conditions on bubble deformation, motion velocity, and drag coefficient are simulated in detail. Moreover, different flow patterns are distinguished according to bubble behavior and are found to be associated with channel structure parameters, gravity Reynolds number (ReGr), and Eötvös number (Eo).

Pore-scale investigation of immiscible gas-assisted gravity drainage

Gas-assisted gravity drainage (GAGD) is an effective method of oil recovery that is influenced by the properties of the fluids and formations involved. In this paper, a direct numerical simulation method is employed to investigate immiscible GAGD in an oil-wet porous medium. The interface between oil and gas is tracked via the phase-field method. A series of numerical simulations are performed over a large range of values of various factors (gravity force, capillary force, viscous force, viscosity ratio, and porous medium properties) to investigate the gas flooding process in a porous medium. The results show the oil–gas interface as a transition zone that migrates during the GAGD process. Gravity improves oil–gas interfacial stability because the continuous oil film gravity-assisted hydraulic connection effect can overcome the capillary force in a small pore. The oil displacement process is dominated by gravity instead of the capillary and viscous forces when the gravitational number exceeds 500. Finally, the pore-scale dimensionless number (Npore) enables a quantitative analysis of the effects of various factors on GAGD. Npore helps predict GAGD oil recovery.

Numerical solution of the conjugate heat transfer problem for turbulent liquid flow in a tube using the large eddy simulation method

This article presents the results of turbulent flow and temperature fluctuations simulation in the liquid and the tube wall in a conjugate formulation using the hybrid Reynolds-averaged Navier-Stokes / large eddy simulation (RANS/LES) method. Calculations are performed for water, air and mercury. The results show that the “conjugate” LES model implemented in in-house CFD code ANES provides a good reproduction of the available experimental results on temperature fluctuations near the wall for cases of different wall thicknesses and thermal activity of the wall material and the fluid. An approximate model for the wall is proposed and can be used as a boundary condition for calculations by direct numerical simulation (DNS) and LES methods in the non-conjugate formulations.

Wall shear-stress extraction by an optical flow algorithm with a sub-grid formulation

In this study, we developed a novel optical-flow algorithm for determining the wall shear-stress on the surface of objects. The algorithm solves the thin-oil-film equation using a numerical scheme that recovers local features neglected by smoothing filters. A variational formulation with a smoothness constraint was applied to extract the global shear-stress fields. The algorithm was then applied to scalar images generated using direct numerical simulation (DNS) method, which revealed that the errors were smaller than those of conventional methods. The application of the proposed algorithm to recover the wall shear-stress on a low-aspect-ratio wing and on an axisymmetric boattail model taken as examples in this study showed a strong potential for analysing shear-stress fields. Compared to the methods used in previous studies, proposed method reveals more local features of separation line and singular points on object surface.

Pore-Network Stitching Method: A Pore-to-Core Upscaling Approach for Multiphase Flow

Pore-network modeling is a widely used predictive tool for pore-scale studies in various applications that deal with multiphase flow in heterogeneous natural rocks. Despite recent improvements to enable pore-network modeling on simplified pore geometry extracted from rock core images and and its computational efficiency compared to direct numerical simulation methods, there are still limitations to modeling a large representative pore-network for heterogeneous cores. These are due to the technical limits on sample size to discern void space during X-ray scanning and computational limits on pore-network extraction algorithms. Thus, there is a need for pore-scale modeling approaches that have the natural advantages of pore-network modeling and can overcome these limitations, thereby enabling better representation of heterogeneity of the 3D complex pore structure and enhancing the accuracy of prediction of macroscopic properties. This paper addresses these issues with a workflow that includes a novel pore-network stitching method to provide large-enough representative pore-network for a core. This workflow uses micro-CT images of heterogeneous reservoir rock cores at different resolutions to characterize the pore structure in order to select few signature parts of the core and extract their equivalent pore-network models. The space between these signature pore-networks is filled by using their statistics to generate realizations of pore-networks which are then connected together using a deterministic layered stitching method. The output of this workflow is a large pore-network that can be used in any flow and transport solver. We validate all steps of this method on different types of natural rocks based on single-phase and two-phase flow properties such as drainage relative permeability curves of carbon dioxide and brine flow. Then, we apply the stochastic workflow on two large domain problems, connecting distant pore-networks and modeling a heterogeneous core. We generate multiple realizations and compare the average results with properties from a defined reference pore-network for each problem. We demonstrate that signature parts of a heterogeneous core, which are a small portion of its entire volume, are sufficient inputs for the developed pore-network stitching method to construct a representative pore-network and predict flow properties.

DeePore: A deep learning workflow for rapid and comprehensive characterization of porous materials

DeePore22GitHub Repository: https://github.com/ArashRabbani/DeePore is a deep learning workflow for rapid estimation of a wide range of porous material properties based on the binarized micro–tomography images. By combining naturally occurring porous textures we generated 17,700 semi–real 3–D micro–structures of porous geo–materials with size of 2563 voxels and 30 physical properties of each sample are calculated using physical simulations on the corresponding pore network models. Next, a designed feed–forward convolutional neural network (CNN) is trained based on the dataset to estimate several morphological, hydraulic, electrical, and mechanical characteristics of the porous material in a fraction of a second. In order to fine–tune the CNN design, we tested 9 different training scenarios and selected the one with the highest average coefficient of determination (R2) equal to 0.885 for 1418 testing samples. Additionally, 3 independent synthetic images as well as 3 realistic tomography images have been tested using the proposed method and results are compared with pore network modelling and experimental data, respectively. Tested absolute permeabilities had around 13% relative error compared to the experimental data which is noticeable considering the accuracy of the direct numerical simulation methods such as Lattice Boltzmann and Finite Volume. The workflow is compatible with any physical size of the images due to its dimensionless approach and can be used to characterize large–scale 3–D images by averaging the model outputs for a sliding window that scans the whole geometry.

Energy transport characteristics of converging Richtmyer–Meshkov instability

In this paper, the Richtmyer–Meshkov (RM) instability in spherical and cylindrical converging geometries with a Mach number of about 1.5 is investigated by using the direct numerical simulation method. The heavy fluid is sulfur hexafluoride, and the light fluid is nitrogen. The shock wave converges from the heavy fluid into the light fluid. The main focus is on the energy transport characteristics in the mixing layer during the entire development process from early instability to late-time turbulent mixing. First, the turbulence kinetic energy transport equation is analyzed, and it is found that the production and dissipation mechanisms of the turbulence induced by the spherical and cylindrical converging RM instabilities in the mixing layer are the same. The turbulent diffusion terms are crucial in the whole development processes of the mixing layers. Before the reflected shock waves transit the interfaces, the dissipation terms can be ignored relative to other terms, and after that, the dissipation terms are close to the production terms and play an important role. The compressibility terms are approximate to the production terms and promote the production of turbulence kinetic energy in the later stage. The viscous diffusion terms can be ignored throughout the process. Then, the enstrophy transport equation is researched, and it is found that, in the mixing layers, the baroclinicity terms play a leading role in the early stage, while the vortex stretching terms play a leading role in the later stage, and the vortex stretching term of the spherical converging geometry develops faster than that of the cylindrical converging geometry. The compressibility terms are positive in the early stage, which promote the production of enstrophy. After the reflected shock waves transit the interfaces, the compressibility terms become negative, which inhibit the production of enstrophy. In addition, the results of the present direct numerical simulation also show that the density fluctuation spectra in the centers of the mixing layers of the spherical and cylindrical converging RM instabilities present the obvious −5/3 scaling law.

Direct Numerical Simulation of Interaction between Transient Flow and Blade Structure in a Modern Low-Pressure Turbine

A feature of a modern aeronautical Low-Pressure Turbine (LPT) is the high blade loadings with complex, transient and separated flow regimes. Most existing research have focused only on analysing the transient flow and flow separation in such turbines. The aerodynamics of a modern LPT, however, can be significantly influenced by the interaction between the unsteady flow field and the blade structure motion in a complex non-linear fashion which could lead to aeroelastic instabilities such as flutter. Therefore, the understanding of the mechanism of the interaction between the flow field unsteadiness and the blade structure in a modern LPT is essential to examine the vibration stress levels to ensure the blade mechanical integrity. The novelty of this paper, first and foremost, is using a high-fidelity direct numerical simulation method to explore the mechanism of flutter and forced response in a modern LPT, T106A turbine, and to study the effects of various sources of unsteadiness on the aeroelastic instabilities of the blade. Secondly, this paper investigates and identifies the adequate working ranges of the harmonic balance method, which has been widely used for the aeromechanical analysis of turbomachines, on predicting the behaviour of the highly unsteady flow due to fluid-structure interaction in an LPT. Another emphasis of this paper is the determination of the capability of the Unsteady Reynolds Averaged Navier–Stokes (URANS) model for the aeroelasticity analysis of an LPT involving the highly unsteady flow. This paper will bridge a key gap in the knowledge of aeroelasticity modelling and analysis of modern LPTs.

DNS of secondary flows over oscillating low-pressure turbine using spectral/hp element method

This paper investigates the secondary vortex flows over an oscillating low-pressure turbine blade using a direct numerical simulation (DNS) method. The unsteady flow governing equations over the oscillating blade are discretized and solved using a spectral/hp element method. Themethodemploys high-degree piecewise polynomial basis functions which results in a very high-orderfinite element approach. The results show that the blade oscillation can significantly influence the transitional flow structure and the wake profile. It was observed that the separation point over vibrating T106A blades was delayed 4.71% compared to the stationary one at Re = 51,800. Moreover, in the oscillating case, the separated shear layers roll up, break down and shed from the trailing edge. However, the blade vibration imposes additional flow disturbances on the suction surface of the blade before leaving from the trailing edge. Momentum thickness calculations revealed that after flow separation point, the momentum thickness grows rapidly which is due to the inverse flow gradients which generate vortex flows in this area. It was concluded that the additional vortex generations due to the blade vibrations cause higher momentum thickness increment compared to the conventional stationary LPT blade.

Particle Velocity Fluctuations in Viscous Gas with Random Velocity as the Sum of Two Correlated Color Noises

The article focuses on methods for studying the phenomenon of two-phase turbulent flows. The turbulence effect on the movement of solid particles in a viscous gas is under study. Dynamics of particles movement in a gas is written in the Stokes approximation, which allows us to suppose the dynamic relaxation time to be a constant value.The random gas velocity is modeled by the sum of two correlated random noises. It is shown that this approach makes it possible to model noise of any structural complexity. The paper describes two research methods based on fundamentally different Euler and Lagrange approaches to the description of a continuous medium. The first approach uses a well-known generalization of the spectral analysis technique for random processes, a popular method for studying turbulence. The second approach implementation is based on the modern generalizations of the theory of numerical algorithms for solving stochastic ordinary differential equations. The spectral method is used to obtain analytical expressions of correlation functions and variance of random processes describing the velocity of gas and solid particles. The qualitative difference between the correlation of fluctuations of modulated random velocities and the behavior of correlations in the case of a single-component gas velocity composition is analyzed. A method of direct numerical simulation for studied processes based on the numerical solution of a stochastic ordinary differential equations system is proposed and analyzed in detail. An array of statistical data obtained as a result of direct numerical modeling is collected and processed. Analytical results are compared qualitatively with numerical results. The influence of input parameters on the character of turbulent flow is studied. The dynamic relaxation time has a significant effect on the complexity of the autocorrelation function of the particle velocity and the response function of particles to gas velocity fluctuations. It is shown that the obtained functions tend to the known results of the standard theory. The considered methods for describing two-phase turbulent flows hold promise for further research.

Simulation of thermal explosion on the base of solution of the equation for the probability density function of particle temperature fluctuations

The influence of temperature fluctuations of gas on the thermal explosion of particles with internal exothermic reactions is considered. The dependence of the rate of a chemical reaction on the concentration of a gas reagent is taken into account. Based on a rigorous system of heat and mass transfer equations, equations for the average temperature and concentration of the reagent inside the particle are obtained. A study of the dynamics of changes in the concentration of reagents and temperature inside the particle in a random field of temperature is carried out on the basis of the probability density function. Based on the solution of the systems of stochastic ordinary differential equations, a method of direct numerical simulation of random temperatures and reagent concentration in a particle is developed. A method for solving the equation for the probability density function is created. The results of comparing the empirical probability density function obtained by direct numerical simulation and the results of solving the equation for the probability density function are presented.

A consistent method for direct numerical simulation of droplet evaporation

A consistent method for direct numerical simulations of evaporation processes is presented. It is implemented into the CFD code FS3D which solves the incompressible Navier-Stokes equations and is based on a Volume-of-Fluid method. In the new framework convection of the fluid and the gaseous phase is conducted first, whereas the physics of the phase change is captured subsequently in a fully consistent loop.We present several validation cases and show, that the method is able to simulate evaporation processes for a large range of ambient conditions, ranging from supercooled droplets to extreme environments with high evaporation rates. Comparisons to analytical models from literature and other experimental investigations are in very good agreement.Then, the method is applied for several three-dimensional direct numerical simulations of evaporating droplets. We present simulations of supercooled droplets in an airflow and compare them to experimental data. Furthermore, an oscillating droplet with evaporation is simulated and evaluated regarding the influence of the phase change on the droplet oscillation behavior. The potential of the method is also demonstrated for a free falling droplet in a high pressure and temperature atmosphere. Vapor fields and temperature regimes are evaluated in this case and show very good agreement with experiments.

Direct Numerical Simulation of Laminar Natural Convection in a Square Cavity at Different Inclination Angle

The effect of angle of rotation on laminar natural convection inside the square cavity have been observed in this research. It was assumed that left and right walls heated isothermally, whereas other two walls act as adiabatic. This problem was solved by assuming 2-D and by Direct Numerical Simulation (DNS) method using ANSYS Fluent 16.0. A series of DNS simulation were carried out for different inclination angle (𝜃 = 0°~90°) of the cavity at Ra = 103 & 104. It was observed that at average Nusselt number increase up to some value of angle of inclination after that it decrease though this variation is not significant.

Secondary instability of stationary Görtler vortices originating from first/second Mack mode

This work investigates the origination of the secondary instability in Görtler vortices using the linear stability theory, BiGlobal analysis, three-dimensional linear parabolized stability equations (3DLPSEs), and direct numerical simulation (DNS). The flow over a concave wall suffering from the Görtler instability and first/second Mack mode instability is selected. Furthermore, this work simulates the evolution of infinitesimal Mack mode disturbance in a flow perturbed by finite-amplitude Görtler vortices by using DNS and 3DLPSE methods. The 3DLPSE approach accurately predicts the process of Mack mode disturbance evolving into the secondary instability of Görtler vortices, and a perfect agreement with results by DNS is obtained. The results indicate that the secondary instability of stationary Görtler vortices can originate from the first/second Mack mode. The evolutions of first/second Mack mode with different spanwise wavenumbers are performed based on 3DLPSE and compared against the BiGlobal method. The results show that the shape functions and growth rates of disturbances always tend to the results of dominant modes obtained by the BiGlobal method. Because the dominant mode might shift from one to another, the overall evolution cannot be predicted only by the BiGlobal method based on a fixed mode. According to our computations, it is deduced that the Mack modes with the same frequency and symmetric characteristics would finally develop into the secondary instability with similar shapes.

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ГИБКОГО РОТОРА С ДВУМЯ ШАРОВЫМИ АВТОБАЛАНСИРАМИ

Ball auto-balancing devices can to compensate changes of unbalance "on the move" only for rotors operating at supercritical speeds. For automatic balancing of such rotors, classified as flexible rotors, several auto-balancers located in different cross sections of the shaft are necessary. This makes it necessary to account bending fluctuations on studies of dynamics of the rotor with auto-balancers, that is especially important in the design of the real rotors. In view of the complexity of experimental studies of such rotors in the article the method of direct numerical simulation of the dynamics of the flexible rotor system – supports – auto-balances is considered. The methodological basis of this method is the use of a discrete multi-mass rotor model, which is equivalent in dynamic characteristics to a real rotor, and also the equations of dynamics of the system discrete rotor – supports – auto-balancers, obtained in the direct form of recording. For definition of discrete masses and a matrix of coefficients of influence of stiffness of rotor cross-sections it is supposed to use calculations for finite-element model of a real rotor by existing software complexes of the engineering analysis. The mathematical model of the system dynamics obtained by the Lagrange method takes into account the non-stationarity of the rotor rotation speed, the influence of gravity and the rolling friction of the balls in the auto-balancer cages. Verification of the mathematical model was performed by reproducing the published data using a computational model for a two-support single-disk three-mass rotor with a two-ball auto-balancer. For a four-mass rotor with two two-ball auto-balancers, the results of numerical simulation of dynamics for the modes of acceleration, steady-state rotation and deceleration are presented. It is shown that for the system under consideration, only partial auto-balancing takes place in the steady rotation mode, including after a stepwise increase of the imbalance.

Numerical simulation of supersonic flow past a plate with surface material sublimation

A method of direct numerical simulation and a method of solving the boundary-layer equations were applied to parameters of a supersonic boundary layer for a flow past a flat plate (Mach number M = 2) for the case of a plate coated with a sublimation material. The sublimating material is naphthalene (C10 H8). Comparison of results from these two approaches — numerical simulation and solution of a boundary layer under the assumption on the local self-similarity — demonstrated a satisfactory agreement between them. Calculations demonstrated that a higher surface temperature produces a higher mass rate of evaporation. Meanwhile, the total heat flux to the solid wall decreases and the wall temperature is lower than for the case of zero sublimation. Since the molecular mass of naphthalene is by several times higher than the molecular mass of air and due to evaporation-induced wall cooling, we observe a higher density of the mixture of air with the sublimating substance vapor near the wall. This may facilitate a higher stability of the supersonic boundary layer and delays the flow transition to the turbulent state.

Multi-Scale Analysis of Fretting Fatigue in Heterogeneous Materials Using Computational Homogenization

This paper deals with modeling of the phenomenon of fretting fatigue in heterogeneous materials using the multi-scale computational homogenization technique and finite element analysis (FEA). The heterogeneous material for the specimens consists of a single hole model (25% void/cell, 16% void/cell and 10% void/cell) and a four-hole model (25% void/cell). Using a representative volume element (RVE), we try to produce the equivalent homogenized properties and work on a homogeneous specimen for the study of fretting fatigue. Next, the fretting fatigue contact problem is performed for 3 new cases of models that consist of a homogeneous and a heterogeneous part (single hole cell) in the contact area. The aim is to analyze the normal and shear stresses of these models and compare them with the results of the corresponding heterogeneous models based on the Direct Numerical Simulation (DNS) method. Finally, by comparing the computational time and % deviations, we draw conclusions about the reliability and effectiveness of the proposed method.

Characteristics of heat flow as well as process of heat conduction and transport in partitioned thermal convection

How to improve the heat transfer efficiency of the system is always a hot issue in thermal convection research field. It is found that when the partitions are added at equal distances to the classical physical model -- Rayleigh-Bénard convection system and gaps are left between the horizontal plates and partitions, the heat transfer efficiency of the system increases significantly with the number of partitions increasing. The heat transfer efficiency can reach up to 3.1 times that of the non-partition device with the specific geometric parameters. In this paper, the Direct Numerical Simulation (DNS) method is used to simulate the partitioned convection system. The mechanism of the heat transfer enhancement of the system is analyzed by studying the characteristics of the heat flow as well as the heat conduction and transport in the system. After the flow in partitioned convection system is fully developed, the fluid in each channel moves vertically in alternating direction and the upward(downward) channel has a higher(lower) temperature than the average temperature of the cell. Due to the symmetry of the system, only the bottom region of the low temperature channel, the bottom region of the high temperature channel and the gap region connecting these two channels are selected for research. By discussing the lateral and longitudinal transport processes of heat flow in the above three regions, the heat flux in the channels and gap areas of the system are studied by quantitative analysis. The results show that the low-temperature fluid in the channel impacts on the bottom plate and then moves into the gap area; the fluid is continuously heated by the bottom plate and flows out of gap area with high temperature. Finally, the fluid converges in the heat transfer channel and forms a longitudinal jet. The external input heat flux of the area that has only horizontal heat transport is large, accounting for 92% of the total heat flux which is obtained from the bottom plate. The heat flux of gap area accounts for 64% of the total heat flux, but the external input heat flux of the area that has only longitudinal heat transport is smaller. The convection system with different geometric parameters has the effect of enhancing heat transfer efficiency, and enhancement capability of the system is stronger when the height of gaps is smaller. When the number of partitions <i>n</i> = 11 and the height of gap <i>d</i> = 0.01, the value of the global <i>Nu</i> number is largest, <i>Nu</i> = 82, which is more than three times that of the system without partitions.

Two-phase flow of CO2-brine in a heterogeneous sandstone: Characterization of the rock and comparison of the lattice-Boltzmann, pore-network, and direct numerical simulation methods

Understanding the physics of two-phase flow of CO2 and brine in porous geological formations is essential to sequestration of carbon dioxide in deep saline reservoirs, as well as the older problem of enhanced oil recovery from hydrocarbon reservoirs by CO2 injection. A pilot CO2 injection in Decatur, Illinois, was undertaken, with the injection zone being the highly saline and heterogeneous Mt. Simon sandstone, in order to better understand the feasibility of full-scale sequestration process. This paper reports the results of an extensive study of the morphology of the sandstone and its heterogeneity, and simulation of single-phase and two-phase flow of CO2 and brine in the formation's three-dimensional images. As we demonstrate by extensive analysis, the formation is much more heterogeneous than the typical sandstone, such as Berea sandstone. In addition to characterizing the morphology of the sandstone and computing its important flow characteristics, an important goal of the study is to compare the accuracy and computational efficiency of three distinct simulation approaches, namely, the lattice-Boltzmann (LB) approach, direct numerical simulation (DNS) of the governing equations of fluid flow that uses the finite-volume method coupled with the OpenFOAM simulator, and pore-network (PN) simulation. After validating the simulators by comparing the computed relative permeabilities that they produce for Berea sandstone, we simulate displacement of brine by CO2 at low and relatively high capillary numbers, and compute the relative permeabilities and other quantities of interest. We demonstrate that all the three methods provide consistent relative permeability-saturation functions that are in close agreement with one another. However, although the LB and DNS both produce similar relative permeabilities, the DNS approach is computationally more efficient because it simulates drainage by only a single set of computations over the entire saturation range, whereas the LB simulation requires separate simulation for each set saturation. Thus, the question of what method to use for simulating such flow processes at the scale of core plugs should mainly be addressed based on the computational time that one can afford and the computational resources that one has access to. Another important question addressed is the effect of the resolution of the computational grids or lattices used, particularly when one uses the LB method with voxelized images of porous media. We show that, unlike many claims in the past, one may need many lattice units per voxel in order to obtain reliable, lattice-independent results.

Characterization of nucleation of methane hydrate crystals: Interfacial theory and molecular simulation

HypothesisSolutions of water and methane gas at favorable thermodynamic conditions lead to the formation of crystalline methane hydrates. In natural and industrial environments, the nucleation process might occur in the solution’s bulk or at the solid-liquid and liquid-gas interfaces, which evolve into distinct morphologies. A complete molecular level understanding and material characterization of preferred nucleation sites and morphologies is required to inhibit or promote crystallization, as required. MethodologyComputational simulations are utilized in this work in combination with analytical theory to calculate the supersaturation and interfacial tension as the driving force and suppressor, respectively, in the hydrate crystal formation process. We employ accurate molecular dynamics (MD) techniques to obtain critical thermodynamic and mechanical properties, and subsequently, analyze the formation using the classical nucleation theory (CNT). FindingsWe report the interfacial tension at all possible interfaces in water-methane gas solutions. We apply both our direct numerical simulation method and Antonow’s rule to find the tension at the methane hydrate and gas interface, and importantly conclude that Antonow’s rule overestimates the values. We calculate the work of formation and nucleation rate of the methane hydrate with and without additives. The nucleation probabilistically forms in the ranked order of film-shaped, cap-shaped, lens-shaped, and homogeneous. We postulate that the premelting of hydrate crystals at the hydrate-gas interface creates an intermediate quasi-liquid layer, which works in favor of the lens-shaped formation compared to homogeneous cases. However, the subtle difference in surface energy indicates high concentration of water and gas molecules at the interface is the main reason behind lens-shaped clustering. We lastly show that ice properties cannot be used to approximate the hydrate formation work.