From 962f1134dd7eeb87236d362e1a6c634c3897c790 Mon Sep 17 00:00:00 2001
From: brendan-waters <brendan.waters@sydney.edu.au>
Date: Tue, 12 Nov 2024 17:14:59 +1100
Subject: [PATCH] iMEM progress
---
src/lbmpy/boundaries/__init__.py | 8 +-
src/lbmpy/boundaries/boundaryconditions.py | 351 +++++++++++++++++--
src/lbmpy/boundaries/wall_function_models.py | 59 ++++
3 files changed, 393 insertions(+), 25 deletions(-)
diff --git a/src/lbmpy/boundaries/__init__.py b/src/lbmpy/boundaries/__init__.py
index e522740..8814eed 100644
--- a/src/lbmpy/boundaries/__init__.py
+++ b/src/lbmpy/boundaries/__init__.py
@@ -1,12 +1,12 @@
from lbmpy.boundaries.boundaryconditions import (
- UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow, WallFunctionBounce,
+ UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow, WallFunctionBounce, iMEMBounceBack,
ExtrapolationOutflow, NeumannByCopy, NoSlip, NoSlipLinearBouzidi, QuadraticBounceBack, StreamInConstant, FreeSlip)
from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
-from lbmpy.boundaries.wall_function_models import MoninObukhovSimilarityTheory, LogLaw, MuskerLaw, SpaldingsLaw
+from lbmpy.boundaries.wall_function_models import MoninObukhovSimilarityTheory, LogLaw, MuskerLaw, ExplicitMuskerLaw, SpaldingsLaw
-__all__ = ['NoSlip', 'NoSlipLinearBouzidi', 'QuadraticBounceBack', 'FreeSlip', 'WallFunctionBounce',
+__all__ = ['NoSlip', 'NoSlipLinearBouzidi', 'QuadraticBounceBack', 'FreeSlip', 'WallFunctionBounce','iMEMBounceBack',
'UBB', 'FixedDensity',
'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
'DiffusionDirichlet', 'NeumannByCopy', 'StreamInConstant',
'LatticeBoltzmannBoundaryHandling',
- 'MoninObukhovSimilarityTheory', 'LogLaw', 'MuskerLaw', 'SpaldingsLaw']
+ 'MoninObukhovSimilarityTheory', 'LogLaw', 'MuskerLaw', 'ExplicitMuskerLaw', 'SpaldingsLaw']
diff --git a/src/lbmpy/boundaries/boundaryconditions.py b/src/lbmpy/boundaries/boundaryconditions.py
index d1c67eb..457aa51 100644
--- a/src/lbmpy/boundaries/boundaryconditions.py
+++ b/src/lbmpy/boundaries/boundaryconditions.py
@@ -240,7 +240,7 @@ class QuadraticBounceBack(LbBoundary):
self.init_wall_distance = init_wall_distance
self.equilibrium_values_name = "f_eq"
self.inv_dir_symbol = TypedSymbol("inv_dir", create_type("int32"))
-
+ self.unique_id = "_out" # Needed for iMEM has no effect on QBB
super(QuadraticBounceBack, self).__init__(name, calculate_force_on_boundary)
@property
@@ -299,31 +299,24 @@ class QuadraticBounceBack(LbBoundary):
pdf_field_accesses = [f_out(i) for i in range(len(lb_method.stencil))]
pdf_symbols = [sp.Symbol(f"pdf_{i}") for i in range(len(lb_method.stencil))]
- f_xf = sp.Symbol("f_xf")
- f_xf_inv = sp.Symbol("f_xf_inv")
- q = sp.Symbol("q")
- feq = sp.Symbol("f_eq")
- weight = sp.Symbol("w")
- weight_inv = sp.Symbol("w_inv")
- v = [TypedSymbol(f"c_{i}", self.data_type) for i in range(lb_method.stencil.D)]
- v_inv = [TypedSymbol(f"c_inv_{i}", self.data_type) for i in range(lb_method.stencil.D)]
+
+ f_xf = f_out(dir_symbol)
+ f_xf_inv = f_out(inv_dir[dir_symbol])
+
+ q = sp.Symbol(f"q")
+ feq = sp.Symbol(f"f_eq{self.unique_id}")
+
+ weight = weight_of_direction(dir_symbol, lb_method)
+ weight_inv = weight_of_direction(inv, lb_method)
+
one = sp.Float(1.0)
half = sp.Rational(1, 2)
subexpressions = [Assignment(pdf_symbols[i], pdf) for i, pdf in enumerate(pdf_field_accesses)]
- subexpressions.append(Assignment(f_xf, f_out(dir_symbol)))
- subexpressions.append(Assignment(f_xf_inv, f_out(inv_dir[dir_symbol])))
subexpressions.append(Assignment(q, index_field[0]('q')))
- subexpressions.append(Assignment(weight, weight_of_direction(dir_symbol, lb_method)))
- subexpressions.append(Assignment(weight_inv, weight_of_direction(inv, lb_method)))
- for i in range(lb_method.stencil.D):
- offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
- subexpressions.append(Assignment(v[i], offset[i]))
-
- for i in range(lb_method.stencil.D):
- offset = NeighbourOffsetArrays.neighbour_offset(inv, lb_method.stencil)
- subexpressions.append(Assignment(v_inv[i], offset[i]))
+ v = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
+ v_inv = NeighbourOffsetArrays.neighbour_offset(inv, lb_method.stencil)
cqc = lb_method.conserved_quantity_computation
rho = cqc.density_symbol
@@ -683,7 +676,7 @@ class WallFunctionBounce(LbBoundary):
# using wall function model
wall_law_assignments = self.wall_function_model.shear_stress_assignments(
density_symbol=rho_for_tau_wall, velocity_symbol=u_m, shear_stress_symbol=tau_w,
- wall_distance=(wall_distance - sp.Rational(1, 2) * self.dy),
+ wall_distance=(wall_distance + sp.Rational(1, 2) * self.dy),
u_tau_target=self.target_friction_velocity)
result.append(wall_law_assignments)
@@ -721,6 +714,322 @@ class WallFunctionBounce(LbBoundary):
# end class WallFunctionBounce
+class iMEMBounceBack(LbBoundary):
+ """
+ Wall function based on the bounce back idea, cf. :cite:`Asmuth2021`.
+
+ Args:
+ lb_method: LB method which is used for the simulation
+ pdfs: Symbolic representation of the particle distribution functions.
+ normal_direction: Normal direction of the wall. Currently, only straight and axis-aligned walls are supported.
+ wall_function_model: Wall function that is used to retrieve the wall stress :math:`tau_w` during the simulation.
+ See :class:`lbmpy.boundaries.wall_treatment.WallFunctionModel` for more details
+ mean_velocity: Optional field or field access for the mean velocity. As wall functions are typically defined
+ in terms of the mean velocity, it is recommended to provide this variable. Per default, the
+ instantaneous velocity obtained from pdfs is used for the wall function.
+ sampling_shift: Optional sampling shift for the velocity sampling. Can be provided as symbolic variable or
+ integer. In both cases, the user must assure that the sampling shift is at least 1, as sampling
+ in boundary cells is not physical. Per default, a sampling shift of 1 is employed which
+ corresponds to a sampling in the first fluid cell normal to the wall. For lower friction
+ Reynolds numbers, choosing a sampling shift >1 has shown to improve the results for higher
+ resolutions.
+ dt: time discretisation. Usually one in LB units
+ dx: space discretisation. Usually one in LB units (assumes constant in x,y,z i.e. dx=dy=dz)
+
+ target_friction_velocity: A target friction velocity can be given if an estimate is known a priori. This target
+ friction velocity will be used as initial guess for implicit wall functions to ensure
+ convergence of the Newton algorithm.
+ name: Optional name of the boundary.
+ data_type: Floating-point precision. Per default, double.
+ """
+
+ def __init__(self, lb_method, pdfs, normal_direction, wall_function_model,
+ inst_velocity=None, clipped_flag=None, mean_velocity=None, sampling_shift=1, dt=1, dx=1,
+ target_friction_velocity=None, BBObject=None,
+ name=None, data_type='double'):
+
+ self.stencil = lb_method.stencil
+
+ self.pdfs = pdfs
+
+ self.wall_function_model = wall_function_model
+
+ if mean_velocity:
+ if isinstance(mean_velocity, Field):
+ self.mean_velocity = mean_velocity.center_vector
+ elif isinstance(mean_velocity, Field.Access):
+ self.mean_velocity = mean_velocity.field.neighbor_vector(mean_velocity.offsets)
+ else:
+ raise ValueError("Mean velocity field has to be a pystencils Field or Field.Access")
+ else:
+ self.mean_velocity = None
+
+ if isinstance(inst_velocity, Field):
+ self.inst_velocity = inst_velocity.center_vector
+ elif isinstance(inst_velocity, Field.Access):
+ self.inst_velocity = inst_velocity.field.neighbor_vector(inst_velocity.offsets)
+
+ if isinstance(clipped_flag, Field):
+ self.clipped_flag = clipped_flag.center_vector
+ elif isinstance(clipped_flag, Field.Access):
+ self.clipped_flag = clipped_flag.field.neighbor_vector(clipped_flag.offsets)
+
+ if not isinstance(sampling_shift, int):
+ self.sampling_shift = TypedSymbol(sampling_shift.name, np.uint32)
+ else:
+ assert sampling_shift >= 1, "The sampling shift must be greater than 1."
+ self.sampling_shift = sampling_shift
+
+ self.dt = dt
+ self.dx = dx
+
+ self.data_type = data_type
+
+ self.target_friction_velocity = target_friction_velocity
+
+ if len(normal_direction) - normal_direction.count(0) != 1:
+ raise ValueError("Only normal directions for straight walls are supported for example (0, 1, 0) for "
+ "a WallFunctionBounce applied to the southern boundary of the domain")
+
+ self.normal_direction = normal_direction
+
+ assert all([n in [-1, 0, 1] for n in self.normal_direction]), \
+ "Only -1, 0 and 1 allowed for defining the normal direction"
+
+ tangential_component = [int(not n) for n in self.normal_direction]
+
+ self.normal_axis = tangential_component.index(0)
+ self.tangential_axis = list(range(len(self.normal_direction)))
+
+ self.tangential_axis.remove(self.normal_axis)
+
+ self.dim = self.stencil.D
+
+ self.inv_dir_symbol = TypedSymbol("inv_dir", create_type("int32"))
+
+ self.BounceBack = BBObject if BBObject else NoSlip()
+
+ if name is None:
+ name = f"iMEM : {offset_to_direction_string([-x for x in normal_direction])}"
+
+ super(iMEMBounceBack, self).__init__(name, calculate_force_on_boundary=False)
+
+ @property
+ def additional_data(self):
+ """Used internally only. For the NoSlipLinearBouzidi boundary the distance to the obstacle of every
+ direction is needed. This information is stored in the index vector."""
+ return [('q', create_type(self.data_type))]
+
+ def get_additional_code_nodes(self, lb_method):
+ inv_directions = [str(self.stencil.index(inverse_direction(direction))) for direction in self.stencil]
+ dtype = self.inv_dir_symbol.dtype
+ name = self.inv_dir_symbol.name
+
+ return [NeighbourOffsetArrays(self.stencil),
+ LbmWeightInfo(lb_method, self.data_type),
+ TranslationArraysNode([(dtype, name, inv_directions), ], {self.inv_dir_symbol})]
+
+ def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
+ ## Step 0: Needed symbols for offsets and indices
+ # neighbour offset symbols are the stencil directions defined in stencils.py:L130ff.
+ neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, self.stencil)
+
+ weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
+ weight_of_direction = weight_info.weight_of_direction
+
+ u_m = sp.Symbol("u_m")
+ u_m_w = sp.Symbol("u_m_w")
+ tau_w = sp.Symbol("tau_w")
+ wall_stress = sp.symbols("tau_w_x tau_w_y tau_w_z")
+
+ cqc = lb_method.conserved_quantity_computation
+
+ result = []
+
+
+ # Step 2: Compute Bounce-Back pdf's
+ c_f_ijk = self.stencil.stencil_entries
+
+ pdfs_in_wall_dir = []
+ for f in range(self.stencil.Q):
+ if c_f_ijk[f][self.normal_axis] == - self.normal_direction[self.normal_axis]:
+ pdfs_in_wall_dir.append(f)
+ print(f, c_f_ijk[f], c_f_ijk[f][self.normal_axis])
+
+ pdf_symbols = [sp.Symbol(f"pdf_{i}") for i in range(len(lb_method.stencil))]
+ pdf_field_accesses = [f_out(i) for i in range(len(lb_method.stencil))]
+ result += [Assignment(pdf_symbols[i], pdf) for i, pdf in enumerate(pdf_field_accesses)]
+
+ f_ijk_field_accesses = [pdf_symbols[i] for i in pdfs_in_wall_dir]
+
+ f_ijk_bb_field_accesses = [sp.Symbol(f"f_bb_{i}") for i in pdfs_in_wall_dir]
+
+ f_ijk_bb_mainexpressions = []
+ f_ijk_bb_subexpressions = []
+
+ unique_lhs = set( pdf_symbols )
+
+ for f_bb, i in zip(f_ijk_bb_field_accesses, pdfs_in_wall_dir):
+ self.BounceBack.unique_id = i
+ bb_rule = self.BounceBack(f_out, f_in, i, inv_dir, lb_method, index_field, None)
+ f_ijk_bb_mainexpressions.append( Assignment(f_bb, bb_rule.main_assignments[0].rhs ))
+
+ for subexpressions in bb_rule.subexpressions:
+ if subexpressions.lhs not in unique_lhs:
+ unique_lhs.add( subexpressions.lhs )
+ f_ijk_bb_subexpressions.append(subexpressions)
+
+ result += f_ijk_bb_subexpressions
+ result += f_ijk_bb_mainexpressions
+
+ print(f_ijk_field_accesses)
+ print(f_ijk_bb_field_accesses)
+
+ # Step 1: Compute u_tau using reference point
+ # Adjusted normal direction for 2D case
+ normal_direction = (list(self.normal_direction) + [0])[:3]
+
+ if (not self.mean_velocity):
+ pdf_center_vector = sp.Matrix([0] * self.stencil.Q)
+
+ for i in range(self.stencil.Q):
+ pdf_center_vector[i] = self.pdfs[neighbor_offset[0] + normal_direction[0],
+ neighbor_offset[1] + normal_direction[1],
+ neighbor_offset[2] + normal_direction[2]](i)
+
+ eq_equations = cqc.equilibrium_input_equations_from_pdfs(pdf_center_vector)
+ result.append(eq_equations.all_assignments)
+
+ # reference velocity for wall stress calculation
+ if self.mean_velocity:
+ u_for_tau_wall = tuple(u_mean_i.get_shifted(
+ self.sampling_shift * normal_direction[0],
+ self.sampling_shift * normal_direction[1],
+ self.sampling_shift * normal_direction[2]
+ ) for u_mean_i in self.mean_velocity)
+
+ rho_for_tau_wall = sp.Float(1)
+ else:
+ rho_for_tau_wall = cqc.density_symbol
+ u_for_tau_wall = cqc.velocity_symbols
+
+ u_mag = Assignment(u_m, sp.sqrt(sum([u_for_tau_wall[i]**2 for i in self.tangential_axis])))
+ result.append(u_mag)
+
+ # Assumes wall is at distace 0.5 from the fluid cell center
+ wall_distance = (self.sampling_shift + sp.Rational(1, 2) * self.dx)
+
+ # Call wall model instance to get shear stress
+ wall_law_assignments = self.wall_function_model.shear_stress_assignments(
+ density_symbol=rho_for_tau_wall, velocity_symbol=u_m,
+ shear_stress_symbol=tau_w, wall_distance=wall_distance,
+ u_tau_target=self.target_friction_velocity)
+ result += wall_law_assignments
+
+ # for i in self.tangential_axis:
+ # result.append(Assignment(wall_stress[i], u_for_tau_wall[i] / u_m * tau_w ))
+
+ # Could also be:
+ u_norm = sp.Symbol("u_norm")
+ u_mag = Assignment(u_norm, sp.sqrt(sum([self.mean_velocity[i]**2 for i in self.tangential_axis])))
+ result.append(u_mag)
+ for i in self.tangential_axis:
+ result.append(Assignment(wall_stress[i], cqc.velocity_symbols[i] / u_norm * tau_w ))
+
+ # Step 3: Compute F_f and F_wm eqs. 24 and 27
+ # F_fluid
+ F_f = sp.symbols("F_fx F_fy F_fz")
+
+ # Δx^2 for 2D or Δx^3 for 3D, see Kruger, 2016, 10.1007/978-3-319-44649-3 Eqs. 5.79-5.80
+ coeff = (self.dx ** self.stencil.D) / self.dt
+
+ # Eq. 24 in Asmuth2021:
+ result += [ Assignment( F_f[j], coeff * sum( c_f_ijk[i][j] *
+ (f_ijk_field_accesses[idx] + f_ijk_bb_field_accesses[idx])
+ for idx, i in enumerate(pdfs_in_wall_dir) )
+ ) for j in self.tangential_axis
+ ]
+
+ # F_wm
+ # Eq. 27 in Asmuth2021:
+ F_wm = sp.symbols("F_wm_x F_wm_y F_wm_z")
+ result += [Assignment( F_wm[i], wall_stress[i] * self.dx ** 2 ) for i in self.tangential_axis]
+
+ # F_uw
+ # Eq. 26 in Asmuth2021:
+ F_uw = sp.symbols("F_uvw_x F_uvw_y F_uvw_z")
+ result += [Assignment( F_uw[i], F_wm[i] - F_f[i] ) for i in self.tangential_axis]
+
+ # Step 4: Compute u_w
+ wi = lb_method.weights
+ wall_velocity = sp.symbols("u_wx u_wy u_wz")
+
+ c_s_sq = sp.Rational(1, 3)
+
+ dx4_on_dt2 = coeff * ( self.dx / self.dt )
+
+ rho = cqc.density_symbol if cqc.compressible else 1
+
+ # rho = 1
+ velocity_expression = [ dx4_on_dt2 * sum( c_f_ijk[i][k] *
+ (- 2 * wi[i] * rho * sum( wall_velocity[j] * c_f_ijk[i][j] for j in range(self.stencil.D)
+ ) / c_s_sq) for i in pdfs_in_wall_dir
+ ) for k in range(self.stencil.D) ]
+
+ result += [ Assignment( wall_velocity[i], sp.solve(velocity_expression[i] - F_uw[i], wall_velocity[i])[0] ) for i in self.tangential_axis]
+ result.append(Assignment( wall_velocity[self.normal_axis], 0 ))
+
+ clipped_velocity = sp.symbols("clip_u_wx clip_u_wy clip_u_wz")
+
+ u_for_clipping = tuple( 0.5 * u_mean_i.get_shifted( normal_direction[0], normal_direction[1], normal_direction[2]
+ ) for u_mean_i in self.inst_velocity)
+
+ result += [ Assignment( clipped_velocity[i],
+ sp.Piecewise((sp.Min(u_for_clipping[i], sp.Max(-u_for_clipping[i], wall_velocity[i])), sp.Gt(u_for_clipping[i], -u_for_clipping[i])),
+ (sp.Max(u_for_clipping[i], sp.Min(-u_for_clipping[i], wall_velocity[i])), True))
+ ) for i in self.tangential_axis]
+ result.append(Assignment( clipped_velocity[self.normal_axis], wall_velocity[self.normal_axis] ))
+
+ #! NOTE: For debugging only ... so we can see if clipping is being used or not
+ result += [ Assignment( self.clipped_flag[i], clipped_velocity[i] - wall_velocity[i]) for i in range(self.stencil.D)]
+
+ # # Step 5: Add f_ijk_uw and f_ijk_bb
+ # Standard NoSlip Bounce-back
+ self.BounceBack.unique_id = "_out"
+ bb_rule = self.BounceBack(f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, None)
+
+ for subexpressions in bb_rule.subexpressions:
+ if subexpressions.lhs not in unique_lhs:
+ unique_lhs.add( subexpressions.lhs )
+ result.append(subexpressions)
+
+ f_ijk_bb = bb_rule.main_assignments[0].rhs
+
+ # velocity = [v_i for v_i in wall_velocity]
+ velocity = [v_i for v_i in clipped_velocity]
+ shifted_vel_eqs = cqc.equilibrium_input_equations_from_init_values(velocity=velocity)
+ shifted_vel_eqs = shifted_vel_eqs.new_without_subexpressions()
+ velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.velocity_symbols).main_assignments]
+
+ weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
+ weight_of_direction = weight_info.weight_of_direction
+ f_ijk_uw = 2 / c_s_sq * rho * sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) * weight_of_direction(
+ dir_symbol, lb_method)
+
+ # Step 6: Return Assignment
+ # Eq. 20 in Asmuth2021: Note, eq. 20 should not have the 2x on LHS and should be - not +
+ result.append(Assignment(f_in(inv_dir[dir_symbol]), f_ijk_bb - f_ijk_uw) )
+
+ # for ass in result:
+ # print("")
+ # print(ass)
+
+ return result
+
+# end class iMEMBounceBack
+
+
+
class UBB(LbBoundary):
r"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle. Furthermore, a density
at the wall can be implied. The boundary condition is implemented with the following formula:
diff --git a/src/lbmpy/boundaries/wall_function_models.py b/src/lbmpy/boundaries/wall_function_models.py
index 4f62778..44ef439 100644
--- a/src/lbmpy/boundaries/wall_function_models.py
+++ b/src/lbmpy/boundaries/wall_function_models.py
@@ -52,6 +52,65 @@ class MoninObukhovSimilarityTheory(ExplicitWallFunctionModel):
return [Assignment(shear_stress_symbol, u_tau ** 2 * density_symbol)]
+class ExplicitMuskerLaw(ExplicitWallFunctionModel):
+ def __init__(self, viscosity, kappa=0.41, name="Explicit Musker"):
+ self.viscosity = viscosity
+ self.kappa = kappa
+
+ super(ExplicitMuskerLaw, self).__init__(name=name)
+
+ def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
+ velocity_symbol: sp.Symbol, wall_distance, u_tau_target=None):
+
+ # Constants
+ c = 0.23
+ s = 8.347e-4
+
+ one_third = sp.Float(1/3)
+
+ l, a, alpha, beta, beta_sq, u = sp.symbols('l a alpha beta beta_sq u')
+ u_tau = sp.Symbol("eval_utau")
+
+ u_tau_init = u_tau_target if u_tau_target else velocity_symbol / sp.Float(100)
+
+ u_plus = velocity_symbol / u_tau_init
+ y_plus = (wall_distance * u_tau_init) / self.viscosity
+
+ r_plus = u_plus * y_plus
+
+ sqrt_r_plus = sp.sqrt( r_plus )
+
+ log_arg = sp.Symbol("log_arg")
+ log_law = (1/self.kappa) * sp.ln( ( sqrt_r_plus + a ) / a )
+
+ coeff = sp.Symbol("coeff")
+ coeff_exp = alpha / ( a + 4 * alpha )
+
+ first_arg = sp.Symbol("first_arg")
+ first_arg_exp = ( a - 4 * alpha ) * sp.ln( ( a * ((sqrt_r_plus - alpha) ** 2 + beta_sq )) / ( 2 * alpha * ( sqrt_r_plus + a ) **2 ))
+
+ second_arg = sp.Symbol("second_arg")
+ second_arg_exp = ( 2 * alpha * (5*a - 4*alpha) / beta ) * ( sp.atan( (sqrt_r_plus - alpha)/beta ) + sp.atan( alpha/beta ) )
+
+ u_p = log_arg + coeff * ( first_arg + second_arg_exp )
+
+ assignments = [Assignment(l, (( 0.5 * sp.sqrt( ( 4*s + 27*c**3) / (27*c) ) / (c*s) ) + \
+ (( 2*s + 27*c**3 ) / ( 54*s*c**3)))**one_third ),
+ Assignment(a, l + (1 / (9 * l * c**2)) + (1 / (3*c)) ),
+ Assignment(alpha, 0.5*(a - 1/c )),
+ Assignment(beta_sq, 2*a*alpha - alpha**2 ),
+ Assignment(beta, sp.sqrt( beta_sq ) ),
+ Assignment(log_arg, log_law ),
+ Assignment(coeff, coeff_exp ),
+ Assignment(first_arg, first_arg_exp ),
+ Assignment(second_arg, second_arg_exp ),
+ Assignment(u, u_p ),
+ Assignment(u_tau, velocity_symbol / u ),
+ Assignment(shear_stress_symbol, u_tau ** 2 * density_symbol)]
+
+ return assignments
+
+
class ImplicitWallFunctionModel(WallFunctionModel, ABC):
"""
Abstract base class for implicit wall functions that require a Newton procedure to solve for the wall shear stress.
--
GitLab