Skip to content
Snippets Groups Projects
Commit 5b5c7c89 authored by Helen Schottenhamml's avatar Helen Schottenhamml
Browse files

Merge branch 'FixNonNewton' into 'master'

[BugFix] Wrong factor in Casson model

See merge request !119
parents 21b1d786 da67628d
Branches
Tags
1 merge request!119[BugFix] Wrong factor in Casson model
Pipeline #39390 passed
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Casson model for simulating Non Newtonian blood flow # Casson model for simulating Non Newtonian blood flow
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import math import math
from lbmpy.session import * from lbmpy.session import *
from lbmpy.non_newtonian_models import CassonsParameters from lbmpy.non_newtonian_models import CassonsParameters
from lbmpy.relaxationrates import lattice_viscosity_from_relaxation_rate from lbmpy.relaxationrates import lattice_viscosity_from_relaxation_rate
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
The primary source for the implementation of the Casson model can be found [here](https://doi.org/10.1007/s10955-005-8415-x) The primary source for the implementation of the Casson model can be found [here](https://doi.org/10.1007/s10955-005-8415-x)
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
stencil = LBStencil(Stencil.D3Q27) stencil = LBStencil(Stencil.D3Q27)
W = 41 W = 41
L = 1 * W L = 1 * W
domain_size = (L, W, W) domain_size = (L, W, W)
omega = 1.0 omega = 1.0
nu = lattice_viscosity_from_relaxation_rate(omega) nu = lattice_viscosity_from_relaxation_rate(omega)
driving_force = 0.0001 driving_force = 0.0001
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
The only parameter of the model is the so-called yield_stress. The main idea is that no strain rate is observed below some stress. However, this leads to the problem that the modified relaxation rate might no longer lead to stable LBM simulations. Thus an upper and lower limit for the shear relaxation rate must be given. All the parameters are combined in the `CassonsParameters` dataclass The only parameter of the model is the so-called yield_stress. The main idea is that no strain rate is observed below some stress. However, this leads to the problem that the modified relaxation rate might no longer lead to stable LBM simulations. Thus an upper and lower limit for the shear relaxation rate must be given. All the parameters are combined in the `CassonsParameters` dataclass
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# yield stress # yield stress
sigma_y = 2*1e-6 sigma_y = 2*1e-6
parameters = CassonsParameters(yield_stress=sigma_y, omega_min=0.2, omega_max=1.98) parameters = CassonsParameters(yield_stress=sigma_y, omega_min=0.2, omega_max=1.98)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
dh = ps.create_data_handling(domain_size=domain_size, periodicity=(True, False, False)) dh = ps.create_data_handling(domain_size=domain_size, periodicity=(True, False, False))
src = dh.add_array('src', values_per_cell=len(stencil)) src = dh.add_array('src', values_per_cell=len(stencil))
dh.fill('src', 0.0, ghost_layers=True) dh.fill('src', 0.0, ghost_layers=True)
dst = dh.add_array('dst', values_per_cell=len(stencil)) dst = dh.add_array('dst', values_per_cell=len(stencil))
dh.fill('dst', 0.0, ghost_layers=True) dh.fill('dst', 0.0, ghost_layers=True)
velField = dh.add_array('velField', values_per_cell=dh.dim) velField = dh.add_array('velField', values_per_cell=dh.dim)
dh.fill('velField', 0.0, ghost_layers=True) dh.fill('velField', 0.0, ghost_layers=True)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
To activate the Casson model in the derivation of the LBM equations, simply pass the parameter class to the LBMConfig To activate the Casson model in the derivation of the LBM equations, simply pass the parameter class to the LBMConfig
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
lbm_config = LBMConfig(stencil=Stencil.D3Q27, method=Method.SRT, lbm_config = LBMConfig(stencil=Stencil.D3Q27, method=Method.SRT,
relaxation_rate=1, compressible=False, force=(driving_force / L, 0,0), relaxation_rate=1, compressible=False, force=(driving_force / L, 0,0),
cassons=parameters, cassons=parameters,
output={'velocity': velField}, kernel_type='stream_pull_collide') output={'velocity': velField}, kernel_type='stream_pull_collide')
method = create_lb_method(lbm_config=lbm_config) method = create_lb_method(lbm_config=lbm_config)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
init = pdf_initialization_assignments(method, 1.0, velocity=velField.center_vector, pdfs=src.center_vector) init = pdf_initialization_assignments(method, 1.0, velocity=velField.center_vector, pdfs=src.center_vector)
ast_init = ps.create_kernel(init, target=dh.default_target) ast_init = ps.create_kernel(init, target=dh.default_target)
kernel_init = ast_init.compile() kernel_init = ast_init.compile()
dh.run_kernel(kernel_init) dh.run_kernel(kernel_init)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
The main expression for the Casson model is: The main expression for the Casson model is:
$$ $$
\lvert \dot{\gamma} \rvert = 2 \sqrt{S_{\alpha \beta} S_{\alpha \beta}} = \frac{3 \sigma}{\tau_\mu \nu^2 \rho \Delta t} \lvert \dot{\gamma} \rvert = 2 \sqrt{S_{\alpha \beta} S_{\alpha \beta}} = \frac{3 \sigma}{\tau_\mu \nu^2 \rho \Delta t}
$$ $$
where $\lvert \dot{\gamma} \rvert$ is the shear strain rate, $\sigma$ is the shear stress, $\nu$ is a constant viscosity $\rho$ is the density and $\Delta t$ is the lattice timestep size. The strain rate tensor is computed in the same way as in `Tutorial 06: Modifying a LBM method: Smagorinsky model`. where $\lvert \dot{\gamma} \rvert$ is the shear strain rate, $\sigma$ is the shear stress, $\nu$ is a constant viscosity $\rho$ is the density and $\Delta t$ is the lattice timestep size. The strain rate tensor is computed in the same way as in `Tutorial 06: Modifying a LBM method: Smagorinsky model`.
In order to calculate the adapted viscosity the following equation is used: In order to calculate the adapted viscosity the following equation is used:
$$ $$
\sqrt{\frac{\mu}{\eta}} = \frac{1}{1 - \theta} \left[1 + \sqrt{\theta \left[ 1 + \frac{\rho \Delta t \nu^2}{\eta} \frac{3}{2} \left( 1 - \theta \right) \right]} \right] \sqrt{\frac{\mu}{\eta}} = \frac{1}{1 - \theta} \left[1 + \sqrt{\theta \left[ 1 + \frac{\rho \Delta t \nu^2}{\eta} \frac{3}{2} \left( 1 - \theta \right) \right]} \right]
$$ $$
where $\theta = \frac{\sigma_y}{\sigma}$ where $\theta = \frac{\sigma_y}{\sigma}$
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
lbm_optimisation = LBMOptimisation(symbolic_field=src, symbolic_temporary_field=dst) lbm_optimisation = LBMOptimisation(symbolic_field=src, symbolic_temporary_field=dst)
update = create_lb_update_rule(lbm_config=lbm_config, update = create_lb_update_rule(lbm_config=lbm_config,
lbm_optimisation=lbm_optimisation) lbm_optimisation=lbm_optimisation)
ast_kernel = ps.create_kernel(update, target=dh.default_target, cpu_openmp=False) ast_kernel = ps.create_kernel(update, target=dh.default_target, cpu_openmp=False)
kernel = ast_kernel.compile() kernel = ast_kernel.compile()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
update # you can look at the complete equation set by removing the comment below
```
%% Output
AssignmentCollection: dst_C^10, dst_C^24, dst_C^4, dst_C^26, dst_C^7, velField_C^1, dst_C^17, dst_C^6, dst_C^25, velField_C^2, dst_C^8, dst_C^19, dst_C^12, dst_C^21, dst_C^14, dst_C^18, velField_C^0, dst_C^9, dst_C^16, dst_C^3, dst_C^5, dst_C^2, dst_C^1, dst_C^11, dst_C^20, dst_C^22, dst_C^13, dst_C^15, dst_C^0, dst_C^23 <- f(src_BE^13, src_TE^17, src_T^6, src_BSE^20, src_TW^18, src_E^3, src_TNE^26, src_TSW^23, src_BW^14, src_BNE^22, src_BSW^19, src_C^0, src_N^2, src_S^1, src_TS^15, src_NE^9, src_B^5, src_SE^7, src_TN^16, src_TNW^25, src_W^4, src_SW^8, src_NW^10, src_TSE^24, src_BS^11, src_BN^12, src_BNW^21) # update
```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def pipe_geometry_callback(x, y, z): def pipe_geometry_callback(x, y, z):
global W global W
radius = W / 2 radius = W / 2
y_mid = W / 2 y_mid = W / 2
z_mid = W / 2 z_mid = W / 2
return (y - y_mid) ** 2 + (z - z_mid) ** 2 > radius ** 2 return (y - y_mid) ** 2 + (z - z_mid) ** 2 > radius ** 2
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
bh = LatticeBoltzmannBoundaryHandling(method, dh, src.name, name="bh") bh = LatticeBoltzmannBoundaryHandling(method, dh, src.name, name="bh")
peridocity = LBMPeriodicityHandling(stencil=stencil, data_handling=dh, pdf_field_name=src.name) peridocity = LBMPeriodicityHandling(stencil=stencil, data_handling=dh, pdf_field_name=src.name)
wall = NoSlip("wall") wall = NoSlip("wall")
bh.set_boundary(wall, mask_callback=pipe_geometry_callback) bh.set_boundary(wall, mask_callback=pipe_geometry_callback)
plt.figure(dpi=200) plt.figure(dpi=200)
plt.boundary_handling(bh, make_slice[0.5, :, :]) plt.boundary_handling(bh, make_slice[0.5, :, :])
``` ```
%% Output %% Output
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def timeloop(timeSteps): def timeloop(timeSteps):
for i in range(timeSteps): for i in range(timeSteps):
bh() bh()
peridocity() peridocity()
dh.run_kernel(kernel) dh.run_kernel(kernel)
dh.swap("src", "dst") dh.swap("src", "dst")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
timeloop(2000) timeloop(2000)
plt.figure(dpi=200) plt.figure(dpi=200)
plt.scalar_field(dh.gather_array('velField')[domain_size[0] // 2, :, :, 0]); plt.scalar_field(dh.gather_array('velField')[domain_size[0] // 2, :, :, 0]);
plt.colorbar() plt.colorbar()
``` ```
%% Output %% Output
<matplotlib.colorbar.Colorbar at 0x130dec1f0> <matplotlib.colorbar.Colorbar at 0x12dd99400>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def poiseuille_profil(r, nu, a, dpdx, dim): def poiseuille_profil(r, nu, a, dpdx, dim):
dim_scale = 2 if dim == 2 else 1 dim_scale = 2 if dim == 2 else 1
return -dim_scale/(4 * nu) * (a**2 - r**2) * dpdx return -dim_scale/(4 * nu) * (a**2 - r**2) * dpdx
def poiseuille_profil_nN(r, rc, nu, a, dpdx, dim): def poiseuille_profil_nN(r, rc, nu, a, dpdx, dim):
dim_scale = 2 if dim == 2 else 1 dim_scale = 2 if dim == 2 else 1
if abs(r) >= rc: if abs(r) >= rc:
p = a**2-r**2-8/3*rc**(1/2)*(a**(3/2)-abs(r)**(3/2))+2*rc*(a-abs(r)) p = a**2-r**2-8/3*rc**(1/2)*(a**(3/2)-abs(r)**(3/2))+2*rc*(a-abs(r))
else: else:
p = a**2-8/3*rc**(1/2)*a**(3/2)+2*rc*a-1/3*rc**2 p = a**2-8/3*rc**(1/2)*a**(3/2)+2*rc*a-1/3*rc**2
return -dim_scale/(4*nu)*dpdx*p return -dim_scale/(4*nu)*dpdx*p
dpdx = -driving_force / L dpdx = -driving_force / L
if stencil.D == 2: if stencil.D == 2:
r_c = abs(1/ math.sqrt(2)*sigma_y/dpdx) r_c = abs(1/ math.sqrt(2)*sigma_y/dpdx)
elif stencil.D == 3: elif stencil.D == 3:
r_c = abs(1/ math.sqrt(2)*2*sigma_y/dpdx) r_c = abs(1/ math.sqrt(2)*2*sigma_y/dpdx)
analytical_solution = [poiseuille_profil(r, nu, W/2, dpdx, stencil.D) for r in range(math.ceil(-W/2), math.floor(W/2)+1)] analytical_solution = [poiseuille_profil(r, nu, W/2, dpdx, stencil.D) for r in range(math.ceil(-W/2), math.floor(W/2)+1)]
analytical_solution_nN = [poiseuille_profil_nN(r, r_c, nu, W/2, dpdx, stencil.D) for r in range(math.ceil(-W/2), math.floor(W/2)+1)] analytical_solution_nN = [poiseuille_profil_nN(r, r_c, nu, W/2, dpdx, stencil.D) for r in range(math.ceil(-W/2), math.floor(W/2)+1)]
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
plt.figure(dpi=200) plt.figure(dpi=200)
plt.scalar_field(dh.gather_array('velField')[:, domain_size[1] // 2, :, 0]); plt.scalar_field(dh.gather_array('velField')[:, domain_size[1] // 2, :, 0]);
plt.colorbar() plt.colorbar()
``` ```
%% Output %% Output
<matplotlib.colorbar.Colorbar at 0x163128160> <matplotlib.colorbar.Colorbar at 0x12f07a610>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
plt.plot(dh.gather_array(velField.name)[domain_size[0]//2, domain_size[1]//2, :, 0]) plt.plot(dh.gather_array(velField.name)[domain_size[0]//2, domain_size[1]//2, :, 0])
plt.plot(analytical_solution_nN) plt.plot(analytical_solution_nN)
plt.plot(analytical_solution) plt.plot(analytical_solution)
plt.legend(["Non-Newtonian flow simulation", "Analytic non-Newtonian flow", "Analytic Newtonian flow"]) plt.legend(["Non-Newtonian flow simulation", "Analytic non-Newtonian flow", "Analytic Newtonian flow"])
``` ```
%% Output %% Output
<matplotlib.legend.Legend at 0x1631c7fd0> <matplotlib.legend.Legend at 0x12f026a30>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
a = dh.gather_array(velField.name)[domain_size[0]//2, domain_size[1]//2, :, 0] a = dh.gather_array(velField.name)[domain_size[0]//2, domain_size[1]//2, :, 0]
b = analytical_solution_nN b = analytical_solution_nN
assert np.max((a - b) / b) < 0.07 assert np.max((a - b) / b) < 0.07
``` ```
......
...@@ -58,7 +58,7 @@ def add_cassons_model(collision_rule, parameter: CassonsParameters, omega_output ...@@ -58,7 +58,7 @@ def add_cassons_model(collision_rule, parameter: CassonsParameters, omega_output
# rhs of equation 14 in https://doi.org/10.1007/s10955-005-8415-x # rhs of equation 14 in https://doi.org/10.1007/s10955-005-8415-x
# Note that C_2 / C_4 = 3 for all configurations thus we directly insert it here # Note that C_2 / C_4 = 3 for all configurations thus we directly insert it here
eq14 = one / (one - theta) * (one + sp.sqrt(theta * (one + rho / eta * sp.Rational(3, 2) * (one - theta)))) eq14 = one / (one - theta) * (one + sp.sqrt(theta * (one + rho / eta * sp.Rational(1, 6) * (one - theta))))
new_omega = one / tau new_omega = one / tau
omega_cond = sp.Piecewise((omega_min, new_omega < omega_min), (omega_max, new_omega > omega_max), (new_omega, True)) omega_cond = sp.Piecewise((omega_min, new_omega < omega_min), (omega_max, new_omega > omega_max), (new_omega, True))
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment