Skip to content
Snippets Groups Projects
Commit da67628d authored by Markus Holzer's avatar Markus Holzer Committed by Helen Schottenhamml
Browse files

[BugFix] Wrong factor in Casson model

parent 21b1d786
No related branches found
No related tags found
1 merge request!119[BugFix] Wrong factor in Casson model
Show whitespace changes
Inline Side-by-side
doc/notebooks/11_tutorial_Non_Newtonian_Flow.ipynb
+ 7
− 19
View file @ da67628d
%% Cell type:markdown id: tags:
# Casson model for simulating Non Newtonian blood flow
%% Cell type:code id: tags:
``` python
import math
from lbmpy.session import *
from lbmpy.non_newtonian_models import CassonsParameters
from lbmpy.relaxationrates import lattice_viscosity_from_relaxation_rate
```
%% 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)
%% Cell type:code id: tags:
``` python
stencil = LBStencil(Stencil.D3Q27)
W = 41
L = 1 * W
domain_size = (L, W, W)
omega = 1.0
nu = lattice_viscosity_from_relaxation_rate(omega)
driving_force = 0.0001
```
%% 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
%% Cell type:code id: tags:
``` python
# yield stress
sigma_y = 2*1e-6
parameters = CassonsParameters(yield_stress=sigma_y, omega_min=0.2, omega_max=1.98)
```
%% Cell type:code id: tags:
``` python
dh = ps.create_data_handling(domain_size=domain_size, periodicity=(True, False, False))
src = dh.add_array('src', values_per_cell=len(stencil))
dh.fill('src', 0.0, ghost_layers=True)
dst = dh.add_array('dst', values_per_cell=len(stencil))
dh.fill('dst', 0.0, ghost_layers=True)
velField = dh.add_array('velField', values_per_cell=dh.dim)
dh.fill('velField', 0.0, ghost_layers=True)
```
%% 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
%% Cell type:code id: tags:
``` python
lbm_config = LBMConfig(stencil=Stencil.D3Q27, method=Method.SRT,
relaxation_rate=1, compressible=False, force=(driving_force / L, 0,0),
cassons=parameters,
output={'velocity': velField}, kernel_type='stream_pull_collide')
method = create_lb_method(lbm_config=lbm_config)
```
%% Cell type:code id: tags:
``` python
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)
kernel_init = ast_init.compile()
dh.run_kernel(kernel_init)
```
%% Cell type:markdown id: tags:
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}
$$
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:
$$
\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}$
%% Cell type:code id: tags:
``` python
lbm_optimisation = LBMOptimisation(symbolic_field=src, symbolic_temporary_field=dst)
update = create_lb_update_rule(lbm_config=lbm_config,
lbm_optimisation=lbm_optimisation)
ast_kernel = ps.create_kernel(update, target=dh.default_target, cpu_openmp=False)
kernel = ast_kernel.compile()
```
%% Cell type:code id: tags:
``` python
update
```
%% Output
# you can look at the complete equation set by removing the comment below
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:
``` python
def pipe_geometry_callback(x, y, z):
global W
radius = W / 2
y_mid = W / 2
z_mid = W / 2
return (y - y_mid) ** 2 + (z - z_mid) ** 2 > radius ** 2
```
%% Cell type:code id: tags:
``` python
bh = LatticeBoltzmannBoundaryHandling(method, dh, src.name, name="bh")
peridocity = LBMPeriodicityHandling(stencil=stencil, data_handling=dh, pdf_field_name=src.name)
wall = NoSlip("wall")
bh.set_boundary(wall, mask_callback=pipe_geometry_callback)
plt.figure(dpi=200)
plt.boundary_handling(bh, make_slice[0.5, :, :])
```
%% Output
%% Cell type:code id: tags:
``` python
def timeloop(timeSteps):
for i in range(timeSteps):
bh()
peridocity()
dh.run_kernel(kernel)
dh.swap("src", "dst")
```
%% Cell type:code id: tags:
``` python
timeloop(2000)
plt.figure(dpi=200)
plt.scalar_field(dh.gather_array('velField')[domain_size[0] // 2, :, :, 0]);
plt.colorbar()
```
%% Output
<matplotlib.colorbar.Colorbar at 0x130dec1f0>
<matplotlib.colorbar.Colorbar at 0x12dd99400>
%% Cell type:code id: tags:
``` python
def poiseuille_profil(r, nu, a, dpdx, dim):
dim_scale = 2 if dim == 2 else 1
return -dim_scale/(4 * nu) * (a**2 - r**2) * dpdx
def poiseuille_profil_nN(r, rc, nu, a, dpdx, dim):
dim_scale = 2 if dim == 2 else 1
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))
else:
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
dpdx = -driving_force / L
if stencil.D == 2:
r_c = abs(1/ math.sqrt(2)*sigma_y/dpdx)
elif stencil.D == 3:
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_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:
``` python
plt.figure(dpi=200)
plt.scalar_field(dh.gather_array('velField')[:, domain_size[1] // 2, :, 0]);
plt.colorbar()
```
%% Output
<matplotlib.colorbar.Colorbar at 0x163128160>
<matplotlib.colorbar.Colorbar at 0x12f07a610>
%% Cell type:code id: tags:
``` python
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)
plt.legend(["Non-Newtonian flow simulation", "Analytic non-Newtonian flow", "Analytic Newtonian flow"])
```
%% Output
<matplotlib.legend.Legend at 0x1631c7fd0>
<matplotlib.legend.Legend at 0x12f026a30>
%% Cell type:code id: tags:
``` python
a = dh.gather_array(velField.name)[domain_size[0]//2, domain_size[1]//2, :, 0]
b = analytical_solution_nN
assert np.max((a - b) / b) < 0.07
```
......
lbmpy/non_newtonian_models.py
+ 1
− 1
View file @ da67628d
......@@ -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
# 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
omega_cond = sp.Piecewise((omega_min, new_omega < omega_min), (omega_max, new_omega > omega_max), (new_omega, True))
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment