Casson model for simulating Non Newtonian blood flow

import math
from lbmpy.session import *
import pystencils as ps

from lbmpy.non_newtonian_models import CassonsParameters
from lbmpy.relaxationrates import lattice_viscosity_from_relaxation_rate

The primary source for the implementation of the Casson model can be found here

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

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

# yield stress
sigma_y = 2*1e-6

parameters = CassonsParameters(yield_stress=sigma_y, omega_min=0.2, omega_max=1.98)

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)

To activate the Casson model in the derivation of the LBM equations, simply pass the parameter class to the LBMConfig

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)

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)

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}\)

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)
kernel = ast_kernel.compile()

# you can look at the complete equation set by removing the comment below

# update

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

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, ps.make_slice[0.5, :, :])

def timeloop(timeSteps):
    for i in range(timeSteps):
        bh()
        peridocity()
        dh.run_kernel(kernel)
        dh.swap("src", "dst")

timeloop(2000)

plt.figure(dpi=200)
plt.scalar_field(dh.gather_array('velField')[domain_size[0] // 2, :, :, 0]);
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7f68c817ca70>

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)]

plt.figure(dpi=200)
plt.scalar_field(dh.gather_array('velField')[:, domain_size[1] // 2, :, 0]);
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7f68c5ab4c80>

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"])

<matplotlib.legend.Legend at 0x7f68c55c1f70>

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