Phase Separation

To simulate multiphase fluids, we need long range interactions between fluid 'particles'.

Shan-Chen Model

For our purposes interactions with nearest neighbor particle densities f will be sufficient to simulate the basic phenomena of multiphase fluid interactions. For single component multiphase fluids (e.g., water/water vapor) an attractive (cohesive) force F between nearest neighbor fluid particles is needed and, for the D2Q9 model, is induced as follows: $$ \boldsymbol{F}(\boldsymbol{x},t)=-G\Psi(\boldsymbol{x},t)\sum_{a=1}^{8}w_a\psi(\boldsymbol{x+e_a}\Delta t,t)\boldsymbol{e_a} $$

where G is the interaction strength, w_a is 1/9 for a={1, 2, 3, 4}, is 1/36 for a={5, 6, 7, 8}, and psi (Greek letter) is the interaction potential:

$$ \Psi(\rho)=\Psi_0exp(\frac{-\rho_0}{\rho}) $$

psi₀ and rho₀ are arbitrary constants. This interaction potential is special in that its "...behavior is consistent with that of an isothermal process..." (Shan and Chen, 1994)

Phase (Liquid-Vapor) Separation Case

In this example, we simulate the phase (Liquid-Vapor) separation. The simulation initialized with an average density of 200 mu lu^-2 with a random variation incorporated via the 'static' initial condition. At the initial condition, we set psi0=4 and rho0=200. This initial density falls on the negatively-sloped, non-physical part of the G=-120 Equations of State (EOS). Therefore, it is unstable and leads to phase separation. G<0 for attraction between particles and the force is stronger when the density is higher. Thus, dense regions (liquid) experience a stronger cohesive force than vapor, which leads to surface tension phenomena.

In this case, the phase separation ultimately leads to a single droplet in a vapor atmosphere. It should be noted that, whether liquid drops or vapor bubbles are formed depends on the total mass in the domain and consequently on the initial density selected.

When phase separation occurs, there is a strong tendency for the interfaces formed to minimized their total area (or length in 2D). This is a straightforward consequence of free energy minimization and occurs in part by geometric rearrangement into the minimum surface area volume (a sphere in 3D or a circle in 2D).

The results at different time steps (0, 100, 200, ..., 100000) are displayed in:

Phase Separation

Code Explanation

Setup a random variation incorporated via the 'static' initial condition.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37


template<typename T, template<typename U> class Descriptor>
class RandomInitializer : public BoxProcessingFunctional2D_L<T,Descriptor> 
{
public :
    RandomInitializer(T rho0_, T maxRho_) : rho0(rho0_), maxRho(maxRho_)
    { };
    virtual void process(Box2D domain,BlockLattice2D<T,Descriptor>& lattice)
    {
        for (plint iX = domain.x0; iX <= domain.x1; ++iX) {
            for (plint iY = domain.y0; iY <= domain.y1; ++iY)
            {
                T rho = rho0 + ((T)rand()/(T)RAND_MAX)*maxRho;
                Array<T,2> zeroVelocity (0.,0.);
                
                iniCellAtEquilibrium(lattice.get(iX,iY), rho, zeroVelocity);
                
                lattice.get(iX,iY).setExternalField (
                        Descriptor<T>::ExternalField::densityBeginsAt,
                        Descriptor<T>::ExternalField::sizeOfDensity, &rho );
                lattice.get(iX,iY).setExternalField (
                        Descriptor<T>::ExternalField::momentumBeginsAt,
                        Descriptor<T>::ExternalField::sizeOfMomentum, &zeroVelocity[0] );
            }
        }
    };
    virtual RandomInitializer<T,Descriptor>* clone() const
    {
        return new RandomInitializer<T,Descriptor>(*this);
    };
    virtual void getTypeOfModification(std::vector<modif::ModifT>& modified) const
    {
        modified[0] = modif::staticVariables;
    };
    
private :
    T rho0, maxRho;
};

1
2
3


    // Use a random initial condition, to activate the phase separation.
    applyProcessingFunctional(new RandomInitializer<T,DESCRIPTOR>(rho0,deltaRho), 
                              lattice.getBoundingBox(),lattice);

The selection of the parameters G, rho0, and psi0:

1
2
3


    const T rho0 = 200.0;
    const T deltaRho = 1.0;
    const T psi0 = 4.0;

The implementation of Shan/Chen interaction potential:

1
2
3
4
5
6


    // Add the data processor which implements the Shan/Chen interaction potential.
    plint processorLevel = 1;
    integrateProcessingFunctional (
            new ShanChenSingleComponentProcessor2D<T,DESCRIPTOR> (
                G, new interparticlePotential::PsiShanChen94<T>(psi0,rho0) ),
            lattice.getBoundingBox(), lattice, processorLevel );

Full Codes

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142


/* This file is part of the Palabos library.
 *
 * The Palabos softare is developed since 2011 by FlowKit-Numeca Group Sarl
 * (Switzerland) and the University of Geneva (Switzerland), which jointly
 * own the IP rights for most of the code base. Since October 2019, the
 * Palabos project is maintained by the University of Geneva and accepts
 * source code contributions from the community.
 * 
 * Contact:
 * Jonas Latt
 * Computer Science Department
 * University of Geneva
 * 7 Route de Drize
 * 1227 Carouge, Switzerland
 * jonas.latt@unige.ch
 *
 * The most recent release of Palabos can be downloaded at 
 * <https://palabos.unige.ch/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "palabos2D.h"
#include "palabos2D.hh"
#include <cstdlib>
#include <iostream>

using namespace plb;
using namespace std;

typedef double T;

#define DESCRIPTOR descriptors::ShanChenD2Q9Descriptor

template<typename T, template<typename U> class Descriptor>
class RandomInitializer : public BoxProcessingFunctional2D_L<T,Descriptor> 
{
public :
    RandomInitializer(T rho0_, T maxRho_) : rho0(rho0_), maxRho(maxRho_)
    { };
    virtual void process(Box2D domain,BlockLattice2D<T,Descriptor>& lattice)
    {
        for (plint iX = domain.x0; iX <= domain.x1; ++iX) {
            for (plint iY = domain.y0; iY <= domain.y1; ++iY)
            {
                T rho = rho0 + ((T)rand()/(T)RAND_MAX)*maxRho;
                Array<T,2> zeroVelocity (0.,0.);
                
                iniCellAtEquilibrium(lattice.get(iX,iY), rho, zeroVelocity);
                
                lattice.get(iX,iY).setExternalField (
                        Descriptor<T>::ExternalField::densityBeginsAt,
                        Descriptor<T>::ExternalField::sizeOfDensity, &rho );
                lattice.get(iX,iY).setExternalField (
                        Descriptor<T>::ExternalField::momentumBeginsAt,
                        Descriptor<T>::ExternalField::sizeOfMomentum, &zeroVelocity[0] );
            }
        }
    };
    virtual RandomInitializer<T,Descriptor>* clone() const
    {
        return new RandomInitializer<T,Descriptor>(*this);
    };
    virtual void getTypeOfModification(std::vector<modif::ModifT>& modified) const
    {
        modified[0] = modif::staticVariables;
    };
    
private :
    T rho0, maxRho;
};


int main(int argc, char *argv[])
{
    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/");

    srand(global::mpi().getRank() + 3);

    // For the choice of the parameters G, rho0, and psi0, we refer to the book
    //   Michael C. Sukop and Daniel T. Thorne (2006), 
    //   Lattice Boltzmann Modeling; an Introduction for Geoscientists and Engineers.
    //   Springer-Verlag Berlin/Heidelberg.
    
    const T omega = 1.0;
    const int nx   = 400;
    const int ny   = 400;
    const T G      = -120.0;
    const int maxIter  = 100001;
    const int saveIter = 100;
    const int statIter = 100;
    
    const T rho0 = 200.0;
    const T deltaRho = 1.0;
    const T psi0 = 4.0;

    MultiBlockLattice2D<T, DESCRIPTOR> lattice (
            nx,ny, new ExternalMomentBGKdynamics<T, DESCRIPTOR>(omega) );
            
    lattice.periodicity().toggleAll(true);

    // Use a random initial condition, to activate the phase separation.
    applyProcessingFunctional(new RandomInitializer<T,DESCRIPTOR>(rho0,deltaRho), 
                              lattice.getBoundingBox(),lattice);

    // Add the data processor which implements the Shan/Chen interaction potential.
    plint processorLevel = 1;
    integrateProcessingFunctional (
            new ShanChenSingleComponentProcessor2D<T,DESCRIPTOR> (
                G, new interparticlePotential::PsiShanChen94<T>(psi0,rho0) ),
            lattice.getBoundingBox(), lattice, processorLevel );

    lattice.initialize();
    
    pcout << "Starting simulation" << endl;
    for (int iT=0; iT<maxIter; ++iT) {
        if (iT%statIter==0) {
            std::unique_ptr<MultiScalarField2D<T> > rho( computeDensity(lattice) );
            pcout << iT << ": Average rho fluid one = " << computeAverage(*rho) << endl;
            pcout << "Minimum density: " << computeMin(*rho) << endl;
            pcout << "Maximum density: " << computeMax(*rho) << endl;
        }
        if (iT%saveIter == 0) {
            ImageWriter<T>("leeloo").writeScaledGif (
                    createFileName("rho", iT, 6), *computeDensity(lattice) );
        }

        lattice.collideAndStream();
    }
}

Reference

[1] \palabos\examples\codesByTopic\shanChenMultiPhase\segregation2D.cpp

[2] Michael C. Sukop and Daniel T. Thorne (2006), "Chapter 6.3 Phase (Liquid-Vapor) Separation and Interface Minimization", Lattice Boltzmann Modeling: an Introduction for Geoscientists and Engineers. Springer-Verlag Berlin/Heidelberg.

Phase Separation

Shan-Chen Model

Phase (Liquid-Vapor) Separation Case

Code Explanation

Full Codes

Reference

相关文章：