Skip to content

GitLab

Commit bf2b31c9 authored by Markus Holzer's avatar Markus Holzer
Browse files

update phase-field allen cahn 

updated the functionality and the testcases of the phase-field model conservative allen cahn
parent 4203fa7f
Hide whitespace changes
Inline Side-by-side
lbmpy/phasefield_allen_cahn/derivatives.py 0 → 100644
View file @ bf2b31c9
def laplacian(phi_field):
r"""
Get a symbolic expression for the laplacian.
Args:
phi_field: the phase-field on which the laplacian is applied
"""
lp_phi = 16.0 * ((phi_field[1, 0, 0]) + (phi_field[-1, 0, 0])
+ (phi_field[0, 1, 0]) + (phi_field[0, -1, 0])
+ (phi_field[0, 0, 1]) + (phi_field[0, 0, -1])) \
+ 1.0 * (
(phi_field[1, 1, 1]) + (phi_field[-1, 1, 1])
+ (phi_field[1, -1, 1]) + (phi_field[-1, -1, 1])
+ (phi_field[1, 1, -1]) + (phi_field[-1, 1, -1])
+ (phi_field[1, -1, -1]) + (phi_field[-1, -1, -1])) \
+ 4.0 * (
(phi_field[1, 1, 0]) + (phi_field[-1, 1, 0])
+ (phi_field[1, -1, 0]) + (phi_field[-1, -1, 0])
+ (phi_field[1, 0, 1]) + (phi_field[-1, 0, 1])
+ (phi_field[1, 0, -1]) + (phi_field[-1, 0, -1])
+ (phi_field[0, 1, 1]) + (phi_field[0, -1, 1])
+ (phi_field[0, 1, -1]) + (phi_field[0, -1, -1])) \
- 152.0 * phi_field[0, 0, 0]
lp_phi = lp_phi / 36
return lp_phi
def isotropic_gradient(phi_field):
r"""
Get a symbolic expression for the isotropic gradient of the phase-field
Args:
phi_field: the phase-field on which the isotropic gradient is applied
"""
grad_phi_x = 16.00 * (phi_field[1, 0, 0] - phi_field[-1, 0, 0])\
+ (phi_field[1, 1, 1] - phi_field[-1, 1, 1] + phi_field[1, -1, 1] - phi_field[-1, -1, 1]
+ phi_field[1, 1, -1] - phi_field[-1, 1, -1] + phi_field[1, -1, -1] - phi_field[-1, -1, -1])\
+ 4.00 * (phi_field[1, 1, 0] - phi_field[-1, 1, 0] + phi_field[1, -1, 0] - phi_field[-1, -1, 0]
+ phi_field[1, 0, 1] - phi_field[-1, 0, 1] + phi_field[1, 0, -1] - phi_field[-1, 0, -1])
grad_phi_y = 16.00 * (phi_field[0, 1, 0] - phi_field[0, -1, 0]) \
+ (phi_field[1, 1, 1] + phi_field[-1, 1, 1] - phi_field[1, -1, 1] - phi_field[-1, -1, 1]
+ phi_field[1, 1, -1] + phi_field[-1, 1, -1] - phi_field[1, -1, -1] - phi_field[-1, -1, -1])\
+ 4.00 * (phi_field[1, 1, 0] + phi_field[-1, 1, 0] - phi_field[1, -1, 0] - phi_field[-1, -1, 0]
+ phi_field[0, 1, 1] - phi_field[0, -1, 1] + phi_field[0, 1, -1] - phi_field[0, -1, -1])
grad_phi_z = 16.00 * (phi_field[0, 0, 1] - phi_field[0, 0, -1]) \
+ (phi_field[1, 1, 1] + phi_field[-1, 1, 1] + phi_field[1, -1, 1] + phi_field[-1, -1, 1]
- phi_field[1, 1, -1] - phi_field[-1, 1, -1] - phi_field[1, -1, -1] - phi_field[-1, -1, -1])\
+ 4.00 * (phi_field[1, 0, 1] + phi_field[-1, 0, 1] - phi_field[1, 0, -1] - phi_field[-1, 0, -1]
+ phi_field[0, 1, 1] + phi_field[0, -1, 1] - phi_field[0, 1, -1] - phi_field[0, -1, -1])
grad = [grad_phi_x / 72, grad_phi_y / 72, grad_phi_z / 72]
return grad
lbmpy/phasefield_allen_cahn/force_model.py
View file @ bf2b31c9
......@@ -14,21 +14,24 @@ class MultiphaseForceModel:
self._rho = rho
def __call__(self, lb_method):
stencil = lb_method.stencil
simple = Simple(self._force)
force_symp = sp.symbols("force_:{}".format(lb_method.dim))
simple = Simple(force_symp)
force = [f / self._rho for f in simple(lb_method)]
moment_matrix = lb_method.moment_matrix
relaxation_rates = sp.Matrix(np.diag(lb_method.relaxation_rates))
mrt_collision_op = moment_matrix.inv() * relaxation_rates * moment_matrix
return -0.5 * mrt_collision_op * sp.Matrix(force) + sp.Matrix(force)
result = -0.5 * mrt_collision_op * sp.Matrix(force) + sp.Matrix(force)
def macroscopic_velocity_shift(self, density):
return default_velocity_shift(self._rho, self._force)
for i in range(0, len(stencil)):
result[i] = result[i].simplify()
# -------------------------------- Helper functions ------------------------------------------------------------------
subs_dict = dict(zip(force_symp, self._force))
for i in range(0, len(stencil)):
result[i] = result[i].subs(subs_dict)
def default_velocity_shift(density, force):
return [f_i / (2 * density) for f_i in force]
return result
lbmpy/phasefield_allen_cahn/kernel_equations.py
View file @ bf2b31c9
from pystencils.field import Field
from lbmpy.creationfunctions import update_with_default_parameters
from lbmpy.fieldaccess import StreamPushTwoFieldsAccessor, CollideOnlyInplaceAccessor
from pystencils.fd.derivation import FiniteDifferenceStencilDerivation
from lbmpy.maxwellian_equilibrium import get_weights
from pystencils import Assignment, AssignmentCollection
from pystencils.simp import sympy_cse_on_assignment_list
from lbmpy.phasefield_allen_cahn.derivatives import laplacian, isotropic_gradient
import sympy as sp
import numpy as np
......@@ -30,15 +33,16 @@ def chemical_potential_symbolic(phi_field, stencil, beta, kappa):
deriv.set_weight((1, 1), sp.Rational(1, 6))
deriv.set_weight((0, -1), sp.Rational(2, 3))
deriv.set_weight((0, 0), sp.Rational(-10, 3))
if len(stencil) == 27:
deriv.set_weight((1, 1, 1), sp.Rational(1, 48))
# assume the stencil is isotropic
res = deriv.get_stencil(isotropic=True)
# assume the stencil is isotropic
res = deriv.get_stencil(isotropic=True)
lap = res.apply(phi_field.center)
else:
lap = laplacian(phi_field)
# get the chemical potential
mu = 4.0 * beta * phi_field.center * (phi_field.center - 1.0) * (phi_field.center - 0.5) - \
kappa * res.apply(phi_field.center)
kappa * lap
return mu
......@@ -50,28 +54,27 @@ def isotropic_gradient_symbolic(phi_field, stencil):
stencil: stencil of the lattice Boltzmann step
"""
dimensions = len(stencil[0])
deriv = FiniteDifferenceStencilDerivation((0, ), stencil)
deriv = FiniteDifferenceStencilDerivation((0,), stencil)
deriv.assume_symmetric(0, anti_symmetric=True)
deriv.assume_symmetric(1, anti_symmetric=False)
if dimensions == 3:
deriv.assume_symmetric(2, anti_symmetric=False)
if len(stencil) == 19:
if len(stencil) == 9:
res = deriv.get_stencil(isotropic=True)
grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), 0]
elif len(stencil) == 19:
deriv.set_weight((0, 0, 0), sp.Integer(0))
deriv.set_weight((1, 0, 0), sp.Rational(1, 6))
deriv.set_weight((1, 1, 0), sp.Rational(1, 12))
elif len(stencil) == 27:
deriv.set_weight((0, 0, 0), sp.Integer(0))
deriv.set_weight((1, 1, 1), sp.Rational(1, 3360))
res = deriv.get_stencil(isotropic=True)
if dimensions == 2:
grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), 0]
else:
res = deriv.get_stencil(isotropic=True)
grad = [res.apply(phi_field.center),
res.rotate_weights_and_apply(phi_field.center, (1, 0)),
res.rotate_weights_and_apply(phi_field.center, (2, 1))]
else:
grad = isotropic_gradient(phi_field)
return grad
......@@ -83,10 +86,10 @@ def normalized_isotropic_gradient_symbolic(phi_field, stencil):
phi_field: the phase-field on which the normalized isotropic gradient is applied
stencil: stencil of the lattice Boltzmann step
"""
isotropic_gradient = isotropic_gradient_symbolic(phi_field, stencil)
tmp = (sum(map(lambda x: x * x, isotropic_gradient_symbolic(phi_field, stencil))) + 1.e-12) ** 0.5
iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
tmp = (sum(map(lambda x: x * x, isotropic_gradient_symbolic(phi_field, stencil))) + 1.e-32) ** 0.5
result = [x / tmp for x in isotropic_gradient]
result = [x / tmp for x in iso_grad]
return result
......@@ -99,8 +102,8 @@ def pressure_force(phi_field, stencil, density_heavy, density_light):
density_heavy: density of the heavier fluid
density_light: density of the lighter fluid
"""
isotropic_gradient = isotropic_gradient_symbolic(phi_field, stencil)
result = list(map(lambda x: -sp.symbols("rho") * ((density_heavy - density_light) / 3) * x, isotropic_gradient))
iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
result = list(map(lambda x: sp.Rational(-1, 3) * sp.symbols("rho") * (density_heavy - density_light) * x, iso_grad))
return result
......@@ -118,51 +121,47 @@ def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy,
stencil = mrt_method.stencil
dimensions = len(stencil[0])
isotropic_gradient = isotropic_gradient_symbolic(phi_field, stencil)
weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
op = sp.symbols("rho")
gamma = mrt_method.get_equilibrium_terms()
gamma = gamma.subs({sp.symbols("rho"): 1})
iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
moment_matrix = mrt_method.moment_matrix
relaxation_rates = sp.Matrix(np.diag(mrt_method.relaxation_rates))
mrt_collision_op = moment_matrix.inv() * relaxation_rates * moment_matrix
rel = mrt_method.relaxation_rates
eq = mrt_method.moment_equilibrium_values
eq = np.array(eq)
g_vals = [lb_velocity_field.center(i) for i, _ in enumerate(stencil)]
m0 = np.dot(moment_matrix, g_vals)
# get the non equilibrium
non_equilibrium = [lb_velocity_field.center(i) - (weights[i] * op + gamma[i] - weights[i])
for i, _ in enumerate(stencil)]
non_equilibrium = np.dot(mrt_collision_op, non_equilibrium)
m = m0 - eq
m = m * rel
non_equilibrium = np.dot(moment_matrix.inv(), m)
stress_tensor = [0] * 6
# Calculate Stress Tensor MRT
for i, d in enumerate(stencil):
stress_tensor[0] += non_equilibrium[i] * (d[0] * d[0])
stress_tensor[1] += non_equilibrium[i] * (d[1] * d[1])
stress_tensor[0] = sp.Add(stress_tensor[0], non_equilibrium[i] * (d[0] * d[0]))
stress_tensor[1] = sp.Add(stress_tensor[1], non_equilibrium[i] * (d[1] * d[1]))
if dimensions == 3:
stress_tensor[2] += non_equilibrium[i] * (d[2] * d[2])
stress_tensor[3] += non_equilibrium[i] * (d[1] * d[2])
stress_tensor[4] += non_equilibrium[i] * (d[0] * d[2])
stress_tensor[2] = sp.Add(stress_tensor[2], non_equilibrium[i] * (d[2] * d[2]))
stress_tensor[3] = sp.Add(stress_tensor[3], non_equilibrium[i] * (d[1] * d[2]))
stress_tensor[4] = sp.Add(stress_tensor[4], non_equilibrium[i] * (d[0] * d[2]))
stress_tensor[5] += non_equilibrium[i] * (d[0] * d[1])
stress_tensor[5] = sp.Add(stress_tensor[5], non_equilibrium[i] * (d[0] * d[1]))
density_difference = density_heavy - density_light
# Calculate Viscous Force MRT
fmx = (0.5 - tau) * (stress_tensor[0] * isotropic_gradient[0]
+ stress_tensor[5] * isotropic_gradient[1]
+ stress_tensor[4] * isotropic_gradient[2]) * density_difference
fmx = (0.5 - tau) * (stress_tensor[0] * iso_grad[0]
+ stress_tensor[5] * iso_grad[1]
+ stress_tensor[4] * iso_grad[2]) * density_difference
fmy = (0.5 - tau) * (stress_tensor[5] * isotropic_gradient[0]
+ stress_tensor[1] * isotropic_gradient[1]
+ stress_tensor[3] * isotropic_gradient[2]) * density_difference
fmy = (0.5 - tau) * (stress_tensor[5] * iso_grad[0]
+ stress_tensor[1] * iso_grad[1]
+ stress_tensor[3] * iso_grad[2]) * density_difference
fmz = (0.5 - tau) * (stress_tensor[4] * isotropic_gradient[0]
+ stress_tensor[3] * isotropic_gradient[1]
+ stress_tensor[2] * isotropic_gradient[2]) * density_difference
fmz = (0.5 - tau) * (stress_tensor[4] * iso_grad[0]
+ stress_tensor[3] * iso_grad[1]
+ stress_tensor[2] * iso_grad[2]) * density_difference
return [fmx, fmy, fmz]
......@@ -177,18 +176,18 @@ def surface_tension_force(phi_field, stencil, beta, kappa):
kappa: coefficient related to surface tension and interface thickness
"""
chemical_potential = chemical_potential_symbolic(phi_field, stencil, beta, kappa)
isotropic_gradient = isotropic_gradient_symbolic(phi_field, stencil)
return [chemical_potential * x for x in isotropic_gradient]
iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
return [chemical_potential * x for x in iso_grad]
def hydrodynamic_force(lb_velocity_field, phi_field, mrt_method, tau,
def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, tau,
density_heavy, density_light, kappa, beta, body_force):
r"""
Get a symbolic expression for the hydrodynamic force
Args:
lb_velocity_field: hydrodynamic distribution function
phi_field: phase-field
mrt_method: mrt lattice boltzmann method used for hydrodynamics
lb_method: Lattice boltzmann method used for hydrodynamics
tau: relaxation time of the hydrodynamic lattice boltzmann step
density_heavy: density of the heavier fluid
density_light: density of the lighter fluid
......@@ -196,10 +195,10 @@ def hydrodynamic_force(lb_velocity_field, phi_field, mrt_method, tau,
kappa: coefficient related to surface tension and interface thickness
body_force: force acting on the fluids. Usually the gravity
"""
stencil = mrt_method.stencil
stencil = lb_method.stencil
dimensions = len(stencil[0])
fp = pressure_force(phi_field, stencil, density_heavy, density_light)
fm = viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy, density_light)
fm = viscous_force(lb_velocity_field, phi_field, lb_method, tau, density_heavy, density_light)
fs = surface_tension_force(phi_field, stencil, beta, kappa)
result = []
......@@ -226,53 +225,127 @@ def interface_tracking_force(phi_field, stencil, interface_thickness):
return result
def get_assignment_list_stream_hydro(lb_vel_field, lb_vel_field_tmp, mrt_method, force_g, velocity, rho):
def get_update_rules_velocity(src_field, u_in, lb_method, force, density):
r"""
Returns an assignment list of the streaming step for the hydrodynamic lattice Boltzmann step. In the assignment list
also the update for the velocity is integrated
Args:
lb_vel_field: source field of velocity distribution function
lb_vel_field_tmp: destination field of the velocity distribution function
mrt_method: lattice Boltzmann method of the hydrodynamic lattice Boltzmann step
force_g: hydrodynamic force
velocity: velocity field
rho: interpolated density of the two fluids
"""
stencil = mrt_method.stencil
Get assignments to update the velocity with a force shift
Args:
src_field: the source field of the hydrodynamic distribution function
u_in: velocity field
lb_method: mrt lattice boltzmann method used for hydrodynamics
force: force acting on the hydrodynamic lb step
density: the interpolated density of the simulation
"""
stencil = lb_method.stencil
dimensions = len(stencil[0])
velocity_symbol_list = [lb_vel_field.center(i) for i, _ in enumerate(stencil)]
velocity_tmp_symbol_list = [lb_vel_field_tmp.center(i) for i, _ in enumerate(stencil)]
moment_matrix = lb_method.moment_matrix
eq = lb_method.moment_equilibrium_values
first_eqs = lb_method.first_order_equilibrium_moment_symbols
indices = list()
for i in range(dimensions):
indices.append(eq.index(first_eqs[i]))
src = [src_field.center(i) for i, _ in enumerate(stencil)]
m0 = np.dot(moment_matrix, src)
update_u = list()
update_u.append(Assignment(sp.symbols("rho"), m0[0]))
g_subs_dic = dict(zip(velocity_symbol_list, velocity_tmp_symbol_list))
u_symp = sp.symbols("u_:{}".format(dimensions))
zw = sp.symbols("zw_:{}".format(dimensions))
for i in range(dimensions):
update_u.append(Assignment(zw[i], u_in.center_vector[i]))
subs_dict = dict(zip(u_symp, zw))
for i in range(dimensions):
update_u.append(Assignment(u_in.center_vector[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
return update_u
def get_collision_assignments_hydro(collision_rule=None, density=1, optimization={}, **kwargs):
r"""
Get collision assignments for the hydrodynamic lattice Boltzmann step. Here the force gets applied in the moment
space. Afterwards the transformation back to the pdf space happens.
Args:
collision_rule: arbitrary collision rule, e.g. created with create_lb_collision_rule
density: the interpolated density of the simulation
optimization: for details see createfunctions.py
"""
params, opt_params = update_with_default_parameters(kwargs, optimization)
lb_method = params['lb_method']
stencil = lb_method.stencil
dimensions = len(stencil[0])
field_data_type = 'float64' if opt_params['double_precision'] else 'float32'
q = len(stencil)
u_in = params['velocity_input']
force = params['force']
a = [0] * dimensions
for i, direction in enumerate(stencil):
for j in range(dimensions):
a[j] += velocity_tmp_symbol_list[i] * direction[j]
if opt_params['symbolic_field'] is not None:
src_field = opt_params['symbolic_field']
elif opt_params['field_size']:
field_size = [s + 2 for s in opt_params['field_size']] + [q]
src_field = Field.create_fixed_size(params['field_name'], field_size, index_dimensions=1,
layout=opt_params['field_layout'], dtype=field_data_type)
else:
src_field = Field.create_generic(params['field_name'], spatial_dimensions=collision_rule.method.dim,
index_shape=(q,), layout=opt_params['field_layout'], dtype=field_data_type)
if opt_params['symbolic_temporary_field'] is not None:
dst_field = opt_params['symbolic_temporary_field']
else:
dst_field = src_field.new_field_with_different_name(params['temporary_field_name'])
moment_matrix = lb_method.moment_matrix
rel = lb_method.relaxation_rates
eq = lb_method.moment_equilibrium_values
first_eqs = lb_method.first_order_equilibrium_moment_symbols
indices = list()
for i in range(dimensions):
indices.append(eq.index(first_eqs[i]))
pull_g = list()
inv_dir = 0
for i, direction in enumerate(stencil):
inv_direction = tuple(-e for e in direction)
pull_g.append(Assignment(lb_vel_field_tmp(i), lb_vel_field[inv_direction](i)))
inv_dir += lb_vel_field[inv_direction](i)
eq = np.array(eq)
g_vals = [src_field.center(i) for i, _ in enumerate(stencil)]
m0 = np.dot(moment_matrix, g_vals)
mf = np.zeros(len(stencil), dtype=object)
for i in range(dimensions):
pull_g.append(Assignment(velocity.center_vector[i], a[i] + 0.5 * force_g[i].subs(g_subs_dic) / rho))
mf[indices[i]] = force[i] / density
m = sp.symbols("m_:{}".format(len(stencil)))
subexpression = [Assignment(sp.symbols("rho"), inv_dir)]
update_m = get_update_rules_velocity(src_field, u_in, lb_method, force, density)
u_symp = sp.symbols("u_:{}".format(dimensions))
for i in range(dimensions):
subexpression.append(Assignment(u_symp[i], velocity.center_vector[i]))
update_m.append(Assignment(u_symp[i], u_in.center_vector[i]))
for i in range(0, len(stencil)):
update_m.append(Assignment(m[i], m0[i] - (m0[i] - eq[i] + mf[i] / 2) * rel[i] + mf[i]))
update_g = list()
var = np.dot(moment_matrix.inv(), m)
if params['kernel_type'] == 'collide_stream_push':
push_accessor = StreamPushTwoFieldsAccessor()
post_collision_accesses = push_accessor.write(dst_field, stencil)
else:
collide_accessor = CollideOnlyInplaceAccessor()
post_collision_accesses = collide_accessor.write(src_field, stencil)
ac_g = AssignmentCollection(main_assignments=pull_g, subexpressions=subexpression)
for i in range(0, len(stencil)):
update_g.append(Assignment(post_collision_accesses[i], var[i]))
ac_g.main_assignments = sympy_cse_on_assignment_list(ac_g.main_assignments)
hydro_lb_update_rule = AssignmentCollection(main_assignments=update_g,
subexpressions=update_m)
return ac_g
return hydro_lb_update_rule
def initializer_kernel_phase_field_lb(phi_field_distributions, phi_field, velocity, mrt_method, interface_thickness):
......@@ -295,7 +368,7 @@ def initializer_kernel_phase_field_lb(phi_field_distributions, phi_field, veloci
gamma = mrt_method.get_equilibrium_terms()
gamma = gamma.subs({sp.symbols("rho"): 1})
gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity.center_vector)})
# create the kernels for the initialization of the g and h field
# create the kernels for the initialization of the h field
h_updates = list()
def scalar_product(a, b):
......
lbmpy/phasefield_allen_cahn/parameter_calculation.py
View file @ bf2b31c9
......@@ -18,16 +18,21 @@ def calculate_parameters_rti(reference_length=256,
reference_length: reference length of the RTI
reference_time: chosen reference time
density_heavy: density of the heavier fluid
capillary_number: capillary number of the simulation
reynolds_number: reynolds number of the simulation
atwood_number: atwood number of the simulation
peclet_number: peclet number of the simulation
density_ratio: density ration of the heavier and the lighter fluid
capillary_number: capillary number of the simulation represents the relative effect of viscous drag
forces versus surface tension forces
reynolds_number: reynolds number of the simulation is the ratio between the viscous forces in a fluid
and the inertial forces
atwood_number: atwood number of the simulation is a dimensionless density ratio
peclet_number: peclet number of the simulation is the ratio of the rate of advection
of a physical quantity by the flow to the rate of diffusion of the same quantity
driven by an appropriate gradient
density_ratio: density ratio of the heavier and the lighter fluid
viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
"""
# calculate the gravitational acceleration and the reference velocity
gravitational_acceleration = reference_length / (reference_time ** 2 * atwood_number)
reference_velocity = math.sqrt(gravitational_acceleration * reference_length)
g = reference_length / (reference_time ** 2 * atwood_number)
reference_velocity = math.sqrt(g * reference_length)
dynamic_viscosity_heavy = (density_heavy * reference_velocity * reference_length) / reynolds_number
dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
......@@ -49,8 +54,62 @@ def calculate_parameters_rti(reference_length=256,
"dynamic_viscosity_light": dynamic_viscosity_light,
"relaxation_time_heavy": relaxation_time_heavy,
"relaxation_time_light": relaxation_time_light,
"gravitational_acceleration": -gravitational_acceleration,
"gravitational_acceleration": -g,
"mobility": mobility,
"surface_tension": surface_tension
}
return parameters
def calculate_dimensionless_rising_bubble(reference_time=18000,
density_heavy=1.0,
bubble_radius=16,
bond_number=1,
reynolds_number=40,
density_ratio=1000,
viscosity_ratio=100):
r"""
Calculate the simulation parameters for a rising bubble. The parameter calculation follows the ideas of Weber and
is based on the Reynolds number. This means the rising velocity of the bubble is implicitly stated with the
Reynolds number
Args:
reference_time: chosen reference time
density_heavy: density of the heavier fluid
bubble_radius: initial radius of the rising bubble
bond_number: also called eötvös number is measuring the importance of gravitational forces compared to
surface tension forces
reynolds_number: reynolds number of the simulation is the ratio between the viscous forces in a fluid
and the inertial forces
density_ratio: density ratio of the heavier and the lighter fluid
viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
"""
bubble_diameter = bubble_radius * 2
g = bubble_diameter / (reference_time ** 2)
density_light = density_heavy / density_ratio
dynamic_viscosity_heavy = (density_heavy * math.sqrt(g * bubble_diameter ** 3)) / reynolds_number
dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
kinematic_viscosity_heavy = dynamic_viscosity_heavy / density_heavy
kinematic_viscosity_light = dynamic_viscosity_light / density_light
relaxation_time_heavy = 3 * kinematic_viscosity_heavy
relaxation_time_light = 3 * kinematic_viscosity_light
surface_tension = (density_heavy - density_light) * g * bubble_diameter ** 2 / bond_number
# calculation of the Morton number
# Mo = gravitational_acceleration * dynamic_viscosity_heavy / (density_heavy * surface_tension ** 3)
parameters = {
"density_light": density_light,
"dynamic_viscosity_heavy": dynamic_viscosity_heavy,
"dynamic_viscosity_light": dynamic_viscosity_light,
"relaxation_time_heavy": relaxation_time_heavy,
"relaxation_time_light": relaxation_time_light,
"gravitational_acceleration": -g,
"surface_tension": surface_tension
}
return parameters
lbmpy_tests/phasefield_allen_cahn/test_codegen_3D.py 0 → 100644
View file @ bf2b31c9
from lbmpy.phasefield_allen_cahn.analytical import analytic_rising_speed
from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_dimensionless_rising_bubble, \
calculate_parameters_rti
from pystencils import fields, AssignmentCollection
from pystencils.simp import sympy_cse
from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
from lbmpy.stencils import get_stencil
from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments