From 74a2219da204982216630d931093b8424fcbe2e4 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Mon, 3 Dec 2018 09:56:35 +0100
Subject: [PATCH] lbmpy phase-field: high density entropic model works :)
---
phasefield/analytical.py | 12 +++--
phasefield/eos.py | 63 ++++++++++++++------------
phasefield/high_density_ratio_model.py | 51 +++++++++++++++++++++
stencils.py | 4 ++
4 files changed, 97 insertions(+), 33 deletions(-)
create mode 100644 phasefield/high_density_ratio_model.py
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 5328a73c..e067c923 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -26,6 +26,7 @@ def symbolic_order_parameters(num_symbols):
def free_energy_functional_3_phases(order_parameters=None, interface_width=interface_width_symbol, transformed=True,
include_bulk=True, include_interface=True, expand_derivatives=True,
kappa=sp.symbols("kappa_:3")):
+ """Free Energy of ternary multi-component model :cite:`Semprebon2016`. """
kappa_prime = tuple(interface_width ** 2 * k for k in kappa)
c = sp.symbols("C_:3")
@@ -287,12 +288,13 @@ def extract_gamma(free_energy, order_parameters):
continue
if len(diff_factors) != 2:
- raise ValueError("Could not new_filtered Λ because of term " + str(product))
+ raise ValueError("Could not determine Λ because of term " + str(product))
indices = sorted([order_parameters.index(d.args[0]) for d in diff_factors])
- result[tuple(indices)] += prod(e for e in product if e.func != Diff)
+ increment = prod(e for e in product if e.func != Diff)
if diff_factors[0] == diff_factors[1]:
- result[tuple(indices)] *= 2
+ increment *= 2
+ result[tuple(indices)] += increment
return result
@@ -352,7 +354,7 @@ def pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transfo
op = order_parameters
def get_entry(i, j):
- p_if = pressure_tensor_interface_component_new(free_energy, op, dim, i, j) if include_interface else 0
+ p_if = pressure_tensor_interface_component(free_energy, op, dim, i, j) if include_interface else 0
if include_bulk:
p_b = pressure_tensor_bulk_component(free_energy, op) if i == j else 0
else:
@@ -385,5 +387,5 @@ def force_from_pressure_tensor(pressure_tensor, functions=None, pbs=None):
def pressure_tensor_bulk_sqrt_term(free_energy, order_parameters, density, c_s_sq=sp.Rational(1, 3)):
- pbs = sp.sqrt(density * c_s_sq - pressure_tensor_bulk_component(free_energy, order_parameters))
+ pbs = sp.sqrt(sp.Abs(density * c_s_sq - pressure_tensor_bulk_component(free_energy, order_parameters)))
return pbs
diff --git a/phasefield/eos.py b/phasefield/eos.py
index fcc8a408..132f818b 100644
--- a/phasefield/eos.py
+++ b/phasefield/eos.py
@@ -1,5 +1,37 @@
import sympy as sp
+from pystencils.cache import disk_cache
+
+
+# ---------------------------------- Equations of state ----------------------------------------------------------------
+
+def carnahan_starling_eos(density, gas_constant, temperature, a, b):
+ """Carnahan Starling equation of state.
+
+ a, b are parameters specific to this equation of state
+ for details see: Equations of state in a lattice Boltzmann model, by Yuan and Schaefer, 2006
+ """
+ e = b * density / 4
+ fraction = (1 + e + e ** 2 - e ** 3) / (1 - e) ** 3
+ return density * gas_constant * temperature * fraction - a * density ** 2
+
+
+def carnahan_starling_critical_temperature(a, b, gas_constant):
+ return 0.3773 * a / b / gas_constant
+
+
+def van_der_walls_eos(density, gas_constant, temperature, a, b):
+ pressure = sp.Symbol("P")
+ vdw = sp.Equality((pressure + a * density ** 2) * (1 / density - b), gas_constant * temperature)
+ return sp.solve(vdw, pressure)[0]
+
+
+def van_der_walls_critical_temperature(a, b, gas_constant):
+ return 8 * a / 27 / b / gas_constant
+
+
+# ----------------------------- Functions operating on equation of states ----------------------------------------------
+
def eos_from_free_energy(free_energy, density):
"""Compute equation of state from free energy"""
@@ -14,12 +46,10 @@ def free_energy_from_eos(eos, density, integration_constant):
density: symbolic! density parameter
integration_constant:
"""
- return (sp.integrate(eos / (density**2), density) + integration_constant) * density
-
-
-# ---------------------------------- Equations of state ----------------------------------------------------------------
+ return (sp.integrate(eos / (density ** 2), density) + integration_constant) * density
+@disk_cache
def maxwell_construction(eos, tolerance=1e-4):
"""Numerical Maxwell construction to find Ï_gas and Ï_liquid for a given equation of state.
@@ -51,11 +81,10 @@ def maxwell_construction(eos, tolerance=1e-4):
max_p, min_p = eos.subs(rho, max_rho), eos.subs(rho, min_rho)
shift_max = max_p * 0.999
shift_min = max(0, min_p)
-
+
c = (shift_max + shift_min) / 2
deviation = tolerance * 2
while abs(deviation) > tolerance:
- print("Deviation", deviation, "Shift", c)
zeros = sp.solve(eos - c)
integral_bounds = (min(zeros), max(zeros))
deviation = get_deviation(float(integral_bounds[0]), float(integral_bounds[1]), float(c))
@@ -66,25 +95,3 @@ def maxwell_construction(eos, tolerance=1e-4):
c = (shift_max + shift_min) / 2
return integral_bounds
-
-
-# To get final free energy:
-# - from maxwell construciton $\rho_{min}$ and $\rho_{max}$
-# - remove slope from free energy function: C determined by $C = - \frac{d}{dÏ} F(C=0) $
-# - energy shift = $F(Ï_{liquid})$ or $F(Ï_{gas})$ (should be equal)
-# - final free energy := $F - F(Ï_{liquid})$
-
-
-def carnahan_starling_eos(density, gas_constant, temperature, a, b):
- """Carnahan Starling equation of state.
-
- a, b are parameters specific to this equation of state
- for details see: Equations of state in a lattice Boltzmann model, by Yuan and Schaefer, 2006
- """
- e = b * density / 4
- fraction = (1 + e + e**2 - e**3) / (1 - e)**3
- return density * gas_constant * temperature * fraction - a * density ** 2
-
-
-def carnahan_starling_critical_temperature(a, b, gas_constant):
- return 0.3773 * a / b / gas_constant
diff --git a/phasefield/high_density_ratio_model.py b/phasefield/high_density_ratio_model.py
new file mode 100644
index 00000000..5d11f4dc
--- /dev/null
+++ b/phasefield/high_density_ratio_model.py
@@ -0,0 +1,51 @@
+import sympy as sp
+from pystencils.fd import Diff
+from lbmpy.phasefield.eos import free_energy_from_eos
+from pystencils.sympyextensions import remove_small_floats
+
+
+def free_energy_high_density_ratio(eos, density, density_gas, density_liquid, c_liquid_1, c_liquid_2, lambdas, kappas):
+ """ Free energy for a ternary system with high density ratio :cite:`Wohrwag2018`
+
+ Args:
+ eos: equation of state, has to depend on exactly on symbol, the density
+ density: symbol for density
+ density_gas: numeric value for gas density (can be obtained by `maxwell_construction`)
+ density_liquid: numeric value for liquid density (can be obtained by `maxwell_construction`)
+ c_liquid_1: symbol for concentration of first liquid phase
+ c_liquid_2: symbol for concentration of second liquid phase
+ lambdas: pre-factors of bulk terms, lambdas[0] multiplies the density term, lambdas[1] the first liquid and
+ lambdas[2] the second liquid phase
+ kappas: pre-factors of interfacial terms, order same as for lambdas
+
+ Returns:
+ free energy expression
+ """
+ assert eos.atoms(sp.Symbol) == {density}
+ # ---- Part 1: contribution of equation of state, ψ_eos
+ symbolic_integration_constant = sp.Dummy(real=True)
+ psi_eos = free_energy_from_eos(eos, density, symbolic_integration_constant)
+ # accuracy problems in free_energy_from_eos can lead to complex solutions for integration constant
+ psi_eos = remove_small_floats(psi_eos, 1e-14)
+
+ # integration constant is determined from the condition ψ(Ï_gas) == ψ(Ï_liquid)
+ equal_psi_condition = psi_eos.subs(density, density_gas) - psi_eos.subs(density, density_liquid)
+ solve_res = sp.solve(equal_psi_condition, symbolic_integration_constant)
+ assert len(solve_res) == 1
+ integration_constant = solve_res[0]
+ psi_eos = psi_eos.subs(symbolic_integration_constant, integration_constant)
+
+ # energy is shifted by ψ_0 = ψ(Ï_gas) which is also ψ(Ï_liquid) by construction
+ psi_0 = psi_eos.subs(density, density_gas)
+
+ # ---- Part 2: standard double well potential as bulk term, and gradient squares as interface term
+ def f(c):
+ return c ** 2 * (1 - c) ** 2
+
+ f_bulk = (lambdas[0] / 2 * (psi_eos - psi_0) +
+ lambdas[1] / 2 * f(c_liquid_1) +
+ lambdas[2] / 2 * f(c_liquid_2))
+ f_interface = (kappas[0] / 2 * Diff(density) ** 2 +
+ kappas[1] / 2 * Diff(c_liquid_1) ** 2 +
+ kappas[2] / 2 * Diff(c_liquid_2) ** 2)
+ return f_bulk + f_interface
diff --git a/stencils.py b/stencils.py
index 17a6b066..98bd3b1e 100644
--- a/stencils.py
+++ b/stencils.py
@@ -31,6 +31,10 @@ get_stencil.data = {
'braunschweig': ((0, 0),
(-1, 1), (-1, 0), (-1, -1), (0, -1),
(1, -1), (1, 0), (1, 1), (0, 1)),
+ 'uk': ((0, 0),
+ (1, 0), (-1, 0), (0, 1), (0, -1),
+ (1, 1), (-1, -1), (-1, 1), (1, -1),
+ )
},
'D3Q15': {
'walberla':
--
GitLab