PC-SAFT vapor-liquid DFT solver

[DiegoTMelfi]

I believe this is related to #87 but this issue is seen in PC-SAFT and vapor-liquid interfaces so I posting it as its own thing.

It seems that the DFT solver fails for components (or mixtures involving components) with very low/negligible vapor pressure (e.g. polymers).

For example, for polyethylene glycol:

MW_PEG = 3e2

identifier = Identifier(name='PEG')
psr = PcSaftRecord(m=0.0289*MW_PEG, sigma= 3.5497, epsilon_k=241.48, kappa_ab=0.001, 
                    epsilon_k_ab=1425.12, na=1, nb=1, mu = 0)
PEG = PureRecord(identifier, molarweight=MW_PEG, model_record=psr)
parameters = PcSaftParameters.new_pure(PEG)

functional = HelmholtzEnergyFunctional.pcsaft(parameters)
EoS = EquationOfState.pcsaft(parameters)
cp = State.critical_point(EoS)
vle = PhaseEquilibrium.pure(functional, 298.15*KELVIN)
interface = PlanarInterface.from_tanh(
    vle=vle, 
    n_grid=1024, 
    l_grid=300*ANGSTROM, 
    critical_temperature=cp.temperature)
initial_density = interface.density
surface_tension = interface.solve().surface_tension

When MW_PEG > 300 (where the PEG vapor phase density is in the order of µmol/m³ or lower), errors such as:

RuntimeError: DFT did not converge within the maximum number of iterations.
RuntimeError: PureFMTAssocFunctional encountered illegal values during the iteration.
RuntimeError: The matrix appears to be singular.

Are retrieved.

I tried playing with the different solvers, solver options, and the l_grid/n_grid parameters, but couldn't find anything that would work consistently.

Are there any workarounds for cases where the component (or one of the components in a mixture) is virtually absent in one of the phases?

[prehner]

Hey, the DFT calculations become numerically harder if you investigate components with very low vapor pressure, because the density in the gas phase becomes close to 0.

You can try to push the limits to some extent with different solver options. When you do that, you might want to check out the content of profile.solver_log for cases that still converge.

When you reach actual polymers (I guess the transition is somewhat fuzzy) the question arises whether the assumptions in the PC-SAFT DFT model are still reasonable for these kinds of molecules.

[g-bauer]

Hey @DiegoTMelfi.

Adding to what Philipp mentioned - using a Newton solver after a Picard iteration with lower damping coefficient and more grid points worked for me (feos v0.5.1). If speed is a concern, you might be able to find an optimum of smaller grid size and smaller damping coefficient.

from feos.si import *
from feos.dft import *
from feos.pcsaft import *

MW_PEG = 3e2

identifier = Identifier(name='PEG')
psr = PcSaftRecord(m=0.0289*MW_PEG, sigma= 3.5497, epsilon_k=241.48, kappa_ab=0.001, 
                    epsilon_k_ab=1425.12, na=1, nb=1, mu = 0)
PEG = PureRecord(identifier, molarweight=MW_PEG, model_record=psr)
parameters = PcSaftParameters.new_pure(PEG)

functional = HelmholtzEnergyFunctional.pcsaft(parameters)
cp = State.critical_point(functional)
vle = PhaseEquilibrium.pure(functional, 298.15*KELVIN)
interface = PlanarInterface.from_tanh(
    vle=vle, 
    n_grid=2048, 
    l_grid=300*ANGSTROM, 
    critical_temperature=cp.temperature)
initial_density = interface.density
solver = DFTSolver().picard_iteration(damping_coefficient=0.01).newton()
surface_tension = interface.solve(solver).surface_tension

Note that you can use a HelmholtzEnergyFunctional as equation of state, so there is no need to create an EquationOfState object to calculate the VLE.

[g-bauer]

Are there any workarounds for cases where the component (or one of the components in a mixture) is virtually absent in one of the phases?

For mixtures, you can use PhaseEquilibrium.tp_flash and specify the index of the non-volatile component in the system. In DFT, non-volatile components can currently not be specifically considered.

[DiegoTMelfi]

Thank you for the quick responses, @prehner and @g-bauer!
Speed is not a concern, but I was hoping for some alternative that could handle really high molecular weights (in the order of thousands) where the non-volatile component is virtually absent. Thank you again!