Molecular dynamic simulation of wet water vapour transport in porous medium
E. G. Nikonov M. Popovičová



 

Molecular dynamic simulation of wet water vapour transport in porous medium
topic JINR LIT No. 05-6-1118-2014/2019, protocol No. 4596-6-17/19. HybriLIT resources.
    E.G. Nikonov^1,  M. Popovičová^2

10cm




^1Joint Institute for Nuclear Research,
 141980 Dubna, Moscow Region, Russia
 email: e.nikonov@jinr.ru 



^2University of Prešov,
 str. Konštantinova 16, 080 01 Prešov,  Slovakia
 email: maria.popovicova@unipo.sk




    
========================================================================================================================================================================================================================================================================




Studies of the transport of wet water vapour are relevant for various areas of human activity, including the construction and production of building materials, mining, agriculture, environmental safety of technological processes, scientific research. In particular, one of the methods for extracting highly viscous bitumen grades of oil from the subterranean depths is based on the dilution of the contents of the porous medium by means of pumping coolants. The simplest and environmentally friendly coolant is wet steam containing small drops of water. In addition, with the help of heated water vapour, the filtering elements of collectors are cleaned from sediments of the solid phase (for example, paraffins, gas hydrates and ice floes) on the walls of the porous medium.
For realistic modeling of the processes of filtration and heat and mass transfer during the injection of wet steam into a porous medium, it is necessary to investigate the characteristics of the interaction of saturated water vapour with individual through-type pores. In this paper, a study was carried out through mathematical modelling of the dependence of the diffusion rate of wet water vapour on the pressure difference outside the pore that occurs when water vapour is injected into a porous medium. The dependencies of the diffusion rate on the pore cross section, the magnitude of vapour adsorption on the pore walls, as well as the effect of water vapour temperature on all these processes were also investigated. Practical interest is the study of the influence of the rate of cooling of water vapour on the diffusion rate and the adsorption of water vapour on the wall of the pores. The calculations were carried out using a hybrid type model that combines molecular dynamics and macro-diffusion approaches to describe the interaction of water vapour with individual pores.

Keywords: porous media, molecular dynamics, wet vapour transport




§ INTRODUCTION


Investigations of wet water vapour transport phenomena are very actual for various areas of human activity, including the construction and production of building materials <cit.>, mining, agriculture, environmental safety of technological processes, scientific research. In particular, one of the methods for extracting highly viscous bitumen grades of oil from the subterranean depths is based on the dilution of the contents of the porous medium by means of pumping coolants. The simplest and environmentally friendly coolant is wet steam containing small drops of water <cit.>. In addition, with the help of heated water vapour, the filtering elements of collectors are cleaned from sediments of the solid phase (for example, paraffins, gas hydrates and ice floes) on  walls of the porous medium.

One of the most interesting phenomenon for scientific research in various fields of engineering as, e.g., building physics and earth sciences is  an absorption of moisture in unsaturated porous media <cit.>. It has been shown <cit.> that, under certain assumptions, the moisture transfer obeys a diffusion law in which the local moisture content is the driving potential, while the diffusivity is a function of the moisture content itself. This diffusion
law is suitable to predict the transient moisture absorption in homogeneous, capillary active, construction materials and soils if gravity is negligible when compared to the capillary forces. Based on the foregoing, for a correct description of porous materials behavior during wetting and drying, it is necessary to know with a sufficient degree of accuracy the dependence of the diffusivity on moisture content in the porous material at given time moment.

For realistic modelling of wet steam penetration in and out porous medium, it is necessary to investigate the characteristics of the interaction of saturated water vapour with individual through-type pores. It could be useful for understanding of the same processes in the bulk porous media with close to uniform distribution of pores or capillaries. In this paper, studies of the dependence of the diffusivity on the moisture content of a porous material are carried out by using molecular dynamic simulations. All calculations have been made for the case of two-dimensionl model of individual through-type pores.  There was also investigated a dependence of the diffusion rate of wet water vapour on the pressure difference outside the pore that occurs when water vapour is injected into porous medium or when the porous material simply  dries or gets wet. Dependences of the diffusion rate on the pore cross section, the magnitude of water vapour adsorption on the pore walls, as well as the effect of water vapour temperature on all these processes were also researched. Practical interest is the study of the influence of the rate of cooling of water vapour on the diffusion rate and the adsorption of water vapour on the wall of the pores. Calculations were carried out using software for molecular dynamic simulations developed in <cit.>.



§ GENERAL PROBLEM

In accordance with the general approach to the description of diffusion processes water vapour passage through a porous media can be done by diffusion equation in the following form.

    ∂ w(𝐫,t)/∂ t=∇·[ D(w,𝐫) ∇ w(𝐫,t)],

where w(𝐫,t) is the density of the water vapour (or humidity) at location 𝐫 and time t and D(w,𝐫) is the generalized diffusion coefficient (or diffusivity) for density w at location 𝐫. If the diffusion coefficient depends on the water vapour density then the equation is nonlinear, otherwise it is linear.

More generally, when D is a symmetric positive definite matrix, the equation describes anisotropic diffusion, which is written (for three dimensional diffusion) as:

    ∂ w(𝐫,t)/∂ t=∑_i=1^3∑_j=1^3∂/∂ x_i[D_ij(w,𝐫)∂ w(𝐫,t)/∂ x_j].


If D is constant, then the equation reduces to the following linear differential equation:

    ∂ w(𝐫,t)/∂ t= D ∇^2 w(𝐫,t).


In the case of isothermal absorption in a material which can be assumed homogeneous at the
macroscopic scale and where gravity is negligible when compared to the capillary forces the diffusion of water vapour through the porous media can be described by one-dimensional moisture transfer model with the following diffusion equation <cit.>.

    ∂ w(x,t)/∂ t= ∂/∂ x(D(w)∂ w(x,t)/∂ x).

Here w(x,t) denotes the local volumetric water content (m^3/m^3) and D(w) the diffusivity
(m^2/s), while the initial and boundary conditions are given by the following equations:

    w(0,x)   =   w_0<1
    
    w(t,0)   =   1
    lim_x→∞ w(t,x)    =    w_0

To solve the problem given by the equation  (<ref>), the diffusivity D(w), which characterizes the
material behavior, has to be known. Experimental data show that this function presents generally a nearly
exponential growth. In fact, approximated solutions for an exponential diffusivity curve have been already published in various studies <cit.>. In all approaches currently existing in the case of the dependence of the diffusivity  on the local volumetric water content in a porous material, direct methods are used to solve the diffusion equation. These can be numerical methods based on the finite difference or finite element methods. This may be a method of approximating the original nonlinear problem and solving it either analytically or numerically under the assumption that the solution of this approximating problem in some sense tends to the solution of the original problem. One of the way to solve the problem  (<ref>) analytically is the change of the original diffusivity function by multiple step function which is equal a constant in each step interval <cit.>. This approach proved to be very effective for solving problem for the one-dimensional diffusion model (<ref>). 

In this paper it is proposed a method for solving problem (<ref>) based on the use of molecular dynamic simulation to obtain the dependence of the diffusion coefficient on the local volumetric water content in a porous material. As assumed above it is considered isothermal absorption in a material which can be assumed homogeneous at the macroscopic scale. It can be assumed that the behavior of water vapour in a homogeneous porous medium should be highly similar to the behavior of water vapour in an individual through-type pore (Fig. <ref>) with specially selected initial and boundary conditions for equation (<ref>). For further consideration the equation (<ref>) is used for description water vapour evolution in the pore with the following initial and boundary conditions.

    w(t,0)=   w_1   
    
    w(t,l_x)=   w_2   
    
    w(0,x)=   w_2-w_1/l_xx+w_1




 


§ MOLECULAR DYNAMIC MODEL

Molecular dynamic simulation of microscopic system evolution is based on representation of object modeling as a system of interacting particles (atoms or molecules). The evolution of the system is a result of a motion of the particles mentioned above. Particle coordinates in each subsequent time step are calculated by integration of equations of motion. These equations contain potentials of particle interactions with each other and with an external environment.

In classical molecular dynamics, the behavior of an individual particle is described by the Newton equations of motion <cit.>, which can be written in the following form

    m_id^2 r⃗_⃗i⃗/dt^2=f⃗_⃗i⃗,

where i  - a particle number, (1≤ i ≤ N), N  - the total number of particles, m_i  - particle mass, r⃗_⃗i⃗ - coordinates of position,  f⃗_⃗i⃗ - the resultant of all forces acting on the particle. This resultant force has the following representation

    f⃗_⃗i⃗ = -∂ U(r⃗_⃗1⃗,…,r⃗_⃗N⃗)/∂r⃗_⃗i⃗ + f⃗_⃗i⃗^ex,

where U  - the potential of particle interaction, f⃗_⃗i⃗^ex - a force caused by external fields. 

It is used the following geometrical configuration (Fig. <ref>) for molecular dynamic model of water vapour diffusion through the pore.


It takes the following boundary conditions for further molecular dynamic simulations. The periodic boundary conditions have been fixed for boundaries Γ_11=[-l_x≤ x ≤ 0; -l_y], Γ_12=[-l_x; -l_y≤ y≤ ly], Γ_13=[-l_x≤ x ≤ 0; l_y], Γ_31=[l_x≤ x ≤ 2l_x; -l_y], Γ_32=[2l_x; -l_y≤ y≤ ly], Γ_33=[l_x≤ x ≤ 2l_x; l_y]. There are no boundaries between Ω_1, Ω_2 and Ω_3 domains. Γ_21=[0≤ x ≤ l_x; -l_y], Γ_22=[0≤ x ≤ l_x; l_y] are absolutely non transparent for water molecules boundaries (”walls”) of the pore.

For a simulation of particle interaction in Ω_1, Ω_2 and Ω_3 domains is used the Lennard-Jones potential <cit.>

    U(r)=4ε[(σ/r)^12-(σ/r)^6]

with σ = 3.17 and ε = 6.74· 10^-3 eV. 

Here r – the distance between the centers of the particles, ε – the depth of the potential well, σ – the distance, where the energy of interaction becomes equal to zero. Parameters ε and σ are characteristic for each sort of atom. The minimum of the potential is reached when r_min = σ√(2).
It is the most used to describe the evolution of water in liquid and saturated vapour form.

The interaction between particles and the ”walls” are described by one of the well-known potential (<ref>) so called Lennard-Jones 9-3 wall potential <cit.>.

    U_9-3(ξ)=2πρ_wσ_wf^3ε_wf[2/45(σ_wf/ξ)^9-1/3(σ_wf/ξ)^3]

Here σ_wf and ε_wf are the Lennard-Jones length and energy parameters for the interaction of a water molecule with a single constituent particle of the solid, and ξ is the perpendicular distance between the fluid molecule and the material surface. In the case of two-dimensional model of a through-type pore ξ is Δ y, the distance from the water molecule to the pore wall along the y coordinate. ρ_w is the volume density of the solid.









Equations of motion (<ref>) were integrated by Velocity Verlet method <cit.>.
Within the framework of the Velocity Verlet method integrating the equations of motion is performed as follows:


  * At the beginning of each step values r(t), v(t), f(t) at the time t are defined or calculated in the previous step.


  * First, the coordinates of the particle's new location are calculated at time t+Δ t , then the velocities of particles  are calculated at time t+Δ t/2


    r(t+Δ t)     =     r(t)+Δ t  v(t)+ (Δ t)^2/2  a(t), 
    
    v(t+Δ t/2)     =     v(t)+Δ t/2 f(t)/m.



  * Next, the forces f(t) acting on the particles are recalculated at time t+Δ t.


  * Finally, the values of velocities are calculated at time t+Δ t


    v(t+Δ t) = v(t+Δ t/2)+Δ t/2 f(t+Δ t)/m.



One of the important thing for the molecular dynamic simulations to account for the effects of energy interchanging with the ambient is a special algorithm which is called a thermostat. In this work Berendsen thermostat <cit.> is used for temperature calibration and control in all domains. 
This thermostat uses alternating nonlinear friction in the equations of motion and is realized by the following equations.

    dr_i(t)/dt   =   v_i(t), 
    
     dv_i(t)/dt   =   f_i(t)/m_i-λ(t)v_i(t).

The coefficient of the velocity recalculation λ(t) at every time step t 

    λ(t)=[1+ Δ t/τ_B(T_0/T(t)-1) ]^1/2.
 
depends on the so called ”rise time” of the thermostat τ_B which belongs to the interval [0.1,2].

 
τ_B describes strength of the coupling of the system to a hypothetical heat bath. For increasing τ_B, the coupling weakens, i.e. it takes longer to achieve given temperature T_0 from current temperature T(t).

The Berendsen algorithm is simple to implement and it is very efficient for reaching the desired temperature from far-from-equilibrium configurations. 

The Andersen thermostat <cit.> is used for simulation of interaction of water molecules with walls of the pore to account for the sticking effect of vapour particles to pore walls. In addition the Andersen thermostat is the simplest one which does correctly sample the NVT ensemble. This thermostat acts on particles entered into a thin layer of width Δ l_y near the pore walls. And at each step, some prescribed number of particles is selected, and their momenta (actually, their velocities) are drawn from a Gaussian distribution at the prescribed temperature.

This is intended to imitate collisions with bath particles at a specified T. The strength of the coupling to the heat bath is specified by a collision frequency, ν. For each particle, a random variate is selected between 0 and 1. If this variate is less than νΔ t, where Δ t is a time step, then that particle's momenta are reset. In this case, the Velocity Verlet integration scheme with the use of the Andersen thermostat can be represented as follows.


  * Calculation of particle position for the time moment t+Δ t.

    𝐫(t+Δ t)=𝐫(t)+𝐯(t)Δ t +𝐟(t)/2m(Δ t)^2


  * Calculation of particle velocity for the time moment t+Δ t/2.

    𝐯(t+Δ t/2)=𝐯(t)+𝐟(t)/2mΔ t


  * And the final calculation of particle velocity for the time moment t+Δ t.

    𝐯(t+Δ t)=𝐯(t+Δ t/2)+𝐟(t+Δ t)/2mΔ t


  * Particle velocity reset.

If ξ < νΔ t, where ξ∈ U[0;1] and U[0;1] numerical uniform distribution, then the value of the velocity 𝐯(t+Δ t) is set to the random number ζ∈ N[1;T], where N[1;T] is Gaussian distribution with variance equal T.



The pressure in the pore was controlled using the formula based on virial equation <cit.>.

    P = 1/3V(⟨ 2K⟩-⟨∑_i<j r_ij· f(r_ij)⟩).

Here V is the pore volume, ⟨ 2K⟩ is the doubled kinetic energy averaged over the ensemble, f(r_ij) is the force between  particles i and j at a distance r_ij.
 


§ COMPUTER SIMULATIONS

In this paper it is used a flat pore model with dimensions of 1μ m × 1μ m × d, where d is the thickness of one layer of water molecules equal to the diameter of the molecule. Since it is assumed that the motion of molecules along the z axis does not occur, this model is actually two-dimensional. The concentration of molecules, thus, is calculated as three-dimensional one, considering that water molecules are located in the volume V = 1μ m × 1μ m × d. Everywhere in this article it is used reduced units which is marked by ^*. 

If we consider that water molecules are located in the volume V^* = 3160 × 3160 × 1 then initial concentrations were obtained from the density of water vapour at the appropriate pressure and density at a given temperature using known tabulated data. Specifically for temperature  T_0=35^oC (T_0^*=3.94) density of saturated water vapour is 0.03962 kg/m^3. In this case  the  pore  contains 420 water molecules which are evenly distributed in the pore at the initial time moment. Initial velocities are randomly selected from [0,1],  then the velocities are set so that the total velocity (momentum) is equal to zero, and then they are changed to achieve the desired initial temperature T_0.

We consider that on the left side of the pore (Fig. <ref>) we have 100 % saturated water vapour for temperature 35^oC. It means that the pressure p_L = 5.622  kPa and density is equal to 0.03962  kg/m^3. On the right side of the pore  temperature is the same, but pressure is lower p_R = 2.337  kPa that means  that there is less number of water molecules (there is not 100 % saturated water vapour) in the right part of the pore. In reduced units we have p_L^* = 1.65183· 10^-4, p_R^* = 0.686647· 10^-4. We assume that pressure changes his value linear from the highest to the lowest value inside the pore	


    p = (1.65183 - 0.965183/3160x)· 10^-4, 0≤ x ≤ 3160.
			

The difference between the pressure values on the left and right ends of the pore and on the right leads to the appearance of external force acting on each particle in the pore. This force is also a linear function of the variable x. It  can be considered as the second part of the sum (<ref>). In general we can compute this force as  F = p· S, where p is pressure and S is the surface to which the pressure is applied. In our example, the surface is a cross-section of the water molecule S=π(d/2)^2 or S^*=π/4 and external force is


    F = π/4(1.65183 - 0.965183/3160x)· 10^-4 .


The conclusion is that every particle in the pore is subject to the force from interaction  with all other particles in the pore by Lennard-Jones potential  and the force interaction with walls of pore by Lennard-Jones 9-3 potential  and external force caused by difference of pressures. 


On the left and the right side of pore  (Fig. <ref>) are implement periodic boundary conditions and up and down side of pore is from solid material from silicon with σ=0.420  nm and ε=32.5818· 10^-23 J, thus σ_w^*=1.327 and ε_w^*=0.3018. We used Lorentz-Berthelot rules  <cit.> to compute parameters for Lennard-Jones 9-3 potential for interaction particles and wall,σ_wf^*=1.1635 and ε_wf^*=0.5494.  Density of pore wall is ρ_w=2330  kg/m^3,   ρ_w^*=21.00541.

Our run has 2000000 time steps with length dt^* = 0.005, we reach result time 16663  ps. The Berendsen thermostat works with parameter τ_B=0.5  ps. The Andersen thermostat uses ν=1 and it works in layer Δ l_y=2σ.

Time evolution of averaged over pore volume diffusivity is presented in Figure <ref>.


Time evolution of the ratio of the number of stuck and moving molecules to the total number of molecules in the pore is presented in Figure <ref>.


Dependence of the average diffusivity on the ratio of the number of stuck and non stuck  molecules to the total number of water molecules are presented in Figures <ref> and <ref> correspondingly.




The two-dimensional density profile D(x,t) is presented in Figure <ref>.




§ CONCLUSIONS

As follows from the simulation results behavior of wet water vapour in the through-type pore is very similar to moisture behavior in homogeneous porous material. Diffusivity is almost exponential time function (Fig. <ref>) and also non linearly, up to model discretization, depends on the δ_adh=N_adh/N_tot, the ratio of stuck water molecules number to total number of water molecules in the pore (Fig. <ref>). It is very similar to dependency of the diffusivity on local volumetric water content inside the pore. Further improvement of the model should lead to more realistic description of wet water vapour behavior inside the through-type pore and will be able to more accurately describe the experimental data for isothermal absorption in a material which can be assumed homogeneous at the macroscopic scale and where gravity is negligible when compared to the capillary forces.

99

[Andersen, 1980]
Andersen1980
Andersen, H. C., Molecular dynamics simulations at constant pressure and/or temperature. // J. Chem. Phys. — 1980, — vol. 72(4), — P. 2384–2393.

[Bianchi Janetti, 2018]
Janetti2018
Bianchi Janetti M., Moisture Absorption in Capillary Active Materials: Analytical Solution for a Multiple Step Diffusivity Function.  // Transp. Porous Med. — 2018, — Vol. 125, — Iss. 3, — pp. 633–645 and references therein.

[Berendsen, 1984]
Berendsen1984
H.J.C.Berendsen, J.P.M.Postma, W.F.van Gunsteren, A.DiNola, J.R.Haak, Molecular dynamics with coupling to an external bath.  // J. Chem. Phys. — 1984, — vol. 81, — P. 3684–3690.
Michele Bianchi Janetti

[Gould, 2005]
Gould Gould H., Tobochnik J., Christian W., An Introduction to Computer Simulation Methods, Chapter 8. Third edition, 2005, pp. 267-268.

[Frenkel and Smith, 2002]
Frenkel2002 Frenkel D., Smith B., Understanding molecular simulation : from algorithms to applications. Second edition, Academic Press, 2006, 658 pp. 

[Krus, 1996]
Krus1996 Krus M., Moisture transport and storage coefficients of porous mineral building materials: Theoretical principals and new test method., Fraunhofer IRB Verlag, Stuttgart, 1996.–175 pp. 

[Lennard-Jones, 1924]
LJ
J. E. Lennard-Jones, On the Determination of Molecular Fields. // Proc. Roy. Soc. — 1924, — vol. A 106, — P. 463–477.

[Nezbeda, 2003]
Nezbeda2003
Nezbeda, Ivo, Jiří Kolafa a Miroslav Kotrla, Úvod do počítačových simulací: metody Monte Carlo a molekulární dynamiky. 2. upr. vyd. Praha: Karolinum, 2003, 204 s.

[NPP, arXiv:1709]

NPP1709
Nikonov E.G., Pavluš M., Popovičová M., 2D microscopic and macroscopic simulation of water and porous material interaction. — 2017, — arXiv:1709.05878 [physics.flu-dyn]

[NPP, arXiv:1708]

NPP1708
Nikonov E.G., Pavluš M., Popovičová M., Molecular dynamic simulation of water vapor interaction with blind pore of dead-end and saccate type. — 2017, — arXiv:1708.06216 [physics.flu-dyn]

[Shagapov, 2004]
Shagapov2004
Shagapov V.Sh., Rakhmatillin I.R., Nasyrova L.A., On the theory of injection of wet steam into a porous medium. // High Temperature —    2004, — vol. 42:6, — p. 951–960.

[Siderius D.W. and Gelb L.D., 2011]
Siderius2011
Siderius D.W.  and Gelb L.D., Extension of the Steele 10-4-3 potential for adsorption calculations
in cylindrical, spherical, and other pore geometries. // J. Chem. Phys —2011, — vol. 135, 084703, — p. 1–7.

[Verlet, 1967]
Verlet1967
Verlet L., Computer 'experiments' on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. // Phys. Rev. —    1967, — vol. 159, — p. 98–103.



