diff --git a/doc/sphinx/lbmpy.bib b/doc/sphinx/lbmpy.bib
index c9d012c1062817ae97e0652c1f391894c532f826..80148f03e295f2e04214438b0a6491e3a6fcf091 100644
--- a/doc/sphinx/lbmpy.bib
+++ b/doc/sphinx/lbmpy.bib
@@ -124,10 +124,22 @@ url = {https://doi.org/10.1016/j.jcp.2017.05.040},
doi = {10.1016/j.jcp.2017.05.040},
journal = {J. Comput. Phys.},
month = {nov},
-pages = {862–888},
+pages = {862-888},
numpages = {27}
}
+@article{Alihussein2024,
+author = {Alihussein, Hussein and Prasannakumar, Adarsh and Geier, Martin and Krafczyk, Manfred},
+title = {{Direct cumulant lattice Boltzmann simulations of transitional flow in gyroidal structures including experimental validation}},
+journal = {Computers & Mathematics with Applications},
+volume = {157},
+pages = {159-172},
+year = {2024},
+issn = {0898-1221},
+doi = {10.1016/j.camwa.2023.12.029},
+url = {https://www.sciencedirect.com/science/article/pii/S0898122123005904},
+}
+
@Article{Coreixas2019,
title = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
author = {Coreixas, Christophe and Chopard, Bastien and Latt, Jonas},
diff --git a/src/lbmpy/boundaries/boundaryconditions.py b/src/lbmpy/boundaries/boundaryconditions.py
index d1c67eb20f9f360d4fbaed5b2200b20a4756cf81..f92e5c45193e9055210ee0f6a24e245cf77a5952 100644
--- a/src/lbmpy/boundaries/boundaryconditions.py
+++ b/src/lbmpy/boundaries/boundaryconditions.py
@@ -900,9 +900,9 @@ class ExtrapolationOutflow(LbBoundary):
array is used which stores them.
.. math ::
- f_{\overline{1}jkxyzt} = f_{\overline{1}jk(x - \Delta x)yz(t - \Delta t)} c \theta^{\frac{1}{2}}
+ f_{\overline{1}jkxyzt} = f_{\overline{1}jk(x - \Delta x)yz(t - \Delta t)} c \theta^{\frac{1}{2}}
\frac{\Delta t}{\Delta x} + \left(1 - c \theta^{\frac{1}{2}} \frac{\Delta t}{\Delta x} \right)
- f_{\overline{1}jk(x - \Delta x)yzt}
+ f_{\overline{1}jk(x - \Delta x)yzt} + \gamma_\rho w_{ijk} \Delta \rho
Args:
@@ -917,11 +917,16 @@ class ExtrapolationOutflow(LbBoundary):
positional arguments, specifying the initial density on boundary nodes
initial_velocity: tuple of floating point constants or callback taking spatial coordinates (x, y [,z]) as
positional arguments, specifying the initial velocity on boundary nodes
+ density: The ExtrapolationOutflow does not impose any coniditon on the density (or pressure). A simple
+ correction term to impose a density condition is proposed by :cite:`Alihussein2024`
+ density_correction_factor: If density condition is added to the boundary (see last term of equation) then this
+ value determines how fast the boundary pushes the density to the imposed value.
"""
def __init__(self, normal_direction, lb_method, dt=1, dx=1, name=None,
streaming_pattern='pull', zeroth_timestep=Timestep.BOTH,
- initial_density=None, initial_velocity=None, data_type='double'):
+ initial_density=None, initial_velocity=None, data_type='double',
+ density=None, density_correction_factor=0.01):
self.lb_method = lb_method
self.stencil = lb_method.stencil
@@ -947,6 +952,8 @@ class ExtrapolationOutflow(LbBoundary):
self.equilibrium_calculation = None
self.data_type = data_type
+ self.density = density
+ self.density_correction_factor = density_correction_factor
if initial_density and initial_velocity:
equilibrium = lb_method.get_equilibrium(conserved_quantity_equations=AssignmentCollection([]))
@@ -989,12 +996,17 @@ class ExtrapolationOutflow(LbBoundary):
# Initial fluid cell PDF values
entry['pdf'] = pdf_acc.read_pdf(boundary_data.pdf_array, domain_cell, inv_dir)
entry['pdf_nd'] = get_boundary_cell_pdfs(domain_cell, outflow_cell, inv_dir)
+ if self.density:
+ entry['rho'] = self.density - 1.0 if self.lb_method.conserved_quantity_computation.zero_centered_pdfs else self.density
+
@property
def additional_data(self):
"""Used internally only. For the ExtrapolationOutflow information of the previous PDF values is needed. This
information is stored in the index vector."""
data = [('pdf', create_type(self.data_type)), ('pdf_nd', create_type(self.data_type))]
+ if self.density:
+ data.append(('rho', create_type(self.data_type)))
return data
@property
@@ -1011,7 +1023,11 @@ class ExtrapolationOutflow(LbBoundary):
list containing NeighbourOffsetArrays
"""
- return [NeighbourOffsetArrays(lb_method.stencil)]
+ code_nodes = [NeighbourOffsetArrays(lb_method.stencil)]
+ if self.density:
+ code_nodes.append(LbmWeightInfo(lb_method, self.data_type))
+ return code_nodes
+
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
subexpressions = []
@@ -1021,9 +1037,23 @@ class ExtrapolationOutflow(LbBoundary):
neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
tangential_offset = tuple(offset - normal for offset, normal in zip(neighbor_offset, self.normal_direction))
+ correction_term = 0
+ if self.density:
+ cqc = lb_method.conserved_quantity_computation
+ density_symbol = sp.Symbol("rho")
+ pdf_field_accesses = [f_out[tangential_offset](i) for i in range(len(lb_method.stencil))]
+ density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol})
+ subexpressions.append(density_equations.all_assignments)
+
+ weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
+ weight_of_direction = weight_info.weight_of_direction
+ weight = weight_of_direction(dir_symbol, lb_method)
+ correction_term = self.density_correction_factor * weight * index_field[0]('rho')
+
+
interpolated_pdf_sym = sp.Symbol('pdf_inter')
interpolated_pdf_asm = Assignment(interpolated_pdf_sym, (index_field[0]('pdf') * (self.c * dtdx))
- + ((sp.Number(1) - self.c * dtdx) * index_field[0]('pdf_nd')))
+ + ((sp.Number(1) - self.c * dtdx) * index_field[0]('pdf_nd')) + correction_term)
subexpressions.append(interpolated_pdf_asm)
asm = Assignment(f_in.center(inv_dir[dir_symbol]), interpolated_pdf_sym)
@@ -1035,6 +1065,10 @@ class ExtrapolationOutflow(LbBoundary):
asm = Assignment(index_field[0]('pdf_nd'), interpolated_pdf_sym)
boundary_assignments.append(asm)
+ if self.density:
+ asm = Assignment(index_field[0]('rho'), density_symbol)
+ boundary_assignments.append(asm)
+
return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)