LBE Exercise 2: Phase separation and equations of state


As well as modelling systems with multiple fluids, LBE is able to model multiple phases (e.g. vapour and liquid at equilibrium) for a single fluid. One method of doing this is to apply the Shan/Chen pseudopotential method [Shan1993], which calculates a function of density (pseudopotential, \(\psi\)) for each lattice site. The force acting on the fluid at a given lattice site is related to the gradients of the pseudopotential (calculated using a stencil) and is applied during the collision step.

One effect of the interaction forces is it modifies the equation of state for the lattice fluid:

\[p = \rho c_s^2 + \frac{1}{2} g c_s^2 \psi^2\]

where \(g\) is the interaction strength used as a multiplier for the forces. It is possible to apply a known equation of state by rearranging the above equation to find the required pseudopotential [Yuan2006]:

\[\psi = \sqrt{\frac{2\left(p \left(\rho, T\right) - \rho c_s^2\right)}{g c_s^2}}\]

where the interaction strength \(g\) is set to \(\pm 1\) to ensure the square root in the above expression is evaluated with a positive value.

An alternative method of applying a different equation of state is to use a free-energy LBE method [Swift1996], which modifies the local equilibrium distribution function to incorporate terms for the bulk pressure and interfacial tension (the latter related to density gradients). While this method is more complicated to apply, it is designed to be thermodynamically consistent and it is possible to specify the interfacial tension between phases (which is generally an emergent property for Shan/Chen).


We are going to take a look at how well the Shan/Chen pseudopotential and free-energy LBE approaches can model vapour/liquid equilibrium for a fluid behaving according to the Peng-Robinson cubic equation of state [Peng1976]. The LBE system we are starting with sets the universal gas constant \(R\) to 1, the attraction parameter \(a = \tfrac{2}{49}\), the finite volume parameter \(b = \tfrac{2}{21}\), the acentric (non-sphericity) factor \(\omega = 0.344\) and the system temperature to 0.055 (below the critical value). To encourage the fluid to separate into phases, we can set a random ‘noise’ in the initial density to produce starting density gradients.


Copy the lbin.sys and input files from the directory DEMO/LBE/2D_PhaseSeparation in your copy of DL_MESO to your working directory.

These simulation input files represent a periodic box of a single fluid with a density (or rather, a small range of densities) that lie between the expected densities of the liquid and vapour phases. These phases should form spontaneously, with the liquid phase collecting either into a circular drop or a rectangular layer across the box (depending on how much fluid exists in the box).

Follow the instructions in Setting up DL_MESO (the section ‘Running DL_MESO’) to run the simulation and produce lbout*.vts files.


  1. Run DL_MESO_LBE in your working directory with the supplied input files for 20,000 timesteps and visualise with Paraview. You should see the fluid separate into two static phases: the high-density fluid can be considered liquid, while the low-density fluid is vapour.
    • Try loading the lbin.sys file in the Java GUI (or a text editor) and varying the noise intensity, initial density and/or fluid relaxation time: does this change the obtained liquid and vapour densities?
    • Apply the Plot Over Line filter to a region with a boundary between the phases: how does the fluid density vary as it goes from phase to phase?
  2. Take a look at the fluid velocities in the system at the end of the simulation. Can you identify where the highest velocities are situated (relative to the two phases)? Given the speed of sound in lattice units is the reciprocal of the square root of three (approximately 0.577), are the velocity magnitudes reasonable?
  3. What is the critical temperature? (Refer to Table 5.1 in chapter 5 of the DL_MESO User Manual to find the expression in terms of the equation of state parameters.) Try varying the system temperature towards this value:
    • What happens to the liquid and vapour densities as the temperature changes?
    • What happens to the maximum fluid velocities in the system when you change the system temperature?
  4. One possible way to reduce fluid velocities is to reduce the forces acting on the fluid by rescaling the pressure used to calculate the pseudopotentials [Liu2010]. For the Peng-Robinson equation of state, reducing the universal gas constant and attraction parameter by the same factor should preserve the density ratio of the two phases (if not the actual values). Try reducing these two coefficients to half or quarter of their original values.
    • What effect does this have on fluid velocities at equilibrium?
    • How does this strategy change the width of the interface between the phases?
  5. Now try modifying the mesoscale interaction method from Shan/Chen to Swift free energy. In addition to the equation of state parameters, you will also need to specify an interfacial tension between the phases: I would suggest a value of about 0.01.
    • Run the simulation and visualise the results with Paraview. Compare the results with an equivalent Shan/Chen simulation run: do you still get the same liquid and vapour densities? What has happened to the velocity magnitudes?
    • Try modifying the interfacial tension: what does this do to the shape of the liquid drop? Can you run the simulation with zero interfacial tension?
    • Does rescaling the pressure have a similar effect for free-energy LBE calculations as it does for Shan/Chen?
[Shan1993]X Shan and H Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Physical Review E, 47 1815-1819, 1993, doi: 10.1103/PhysRevE.47.1815.
[Yuan2006]P Yuan and L Schaefer, Equations of state in a lattice Boltzmann model, Physics of Fluids, 18 042101, 2006, doi: 10.1063/1.2187070.
[Swift1996]MR Swift, E Orlandini, WR Osborn and JM Yeomans, Lattice Boltzmann simulations of liquid-gas and binary fluid systems, Physical Review E, 54, 5041-5052, 1996, doi: 10.1103/PhysRevE.54.5041.
[Peng1976]DY Peng and DB Robinson, A new two-constant equation of state, Industrial & Engineering Chemistry Fundamentals, 15, 59-64, 1976, doi: 10.1021/i160057a011.
[Liu2010]M Liu, Z Yu, T Wang, J Wang and LS Fan, A modified pseudopotential for a lattice Boltzmann simulation of bubbly flow, Chemical Engineering Science, 65, 5615-5623, 2010, doi: 10.1016/j.ces.2010.08.014.