mus_nonNewtonian_module Module

This module keeps all information about the nonNewtonian models. Contains routines which calculates non-Newtonian kinematic viscosity according to non-Newtonian model.

It supports three non-Newtonian models: Power law, Carrear Yasuda and Casson. All these models are described in Ashrafizaadeh, M., & Bakhshaei, H. (2009). A comparison of non-Newtonian models for lattice Boltzmann blood flow simulations. Computers and Mathematics with Applications, 58(5), 1045–1054.

For further information about the theory visit the non-newtonian theory page.


Uses


Variables

Type Visibility Attributes Name Initial
integer, public, parameter :: powerLaw = 1

Identifier for Power-Law model
integer, public, parameter :: casson = 2

Identifier for Casson model
integer, public, parameter :: carreauYasuda = 3

Identifier for Carreau-Yasuda model

Abstract Interfaces

abstract interface

Interface to calculate kinematic viscosity for non-Newtonian model. Viscosity is computed from shear rate which is computed from strain rate which is computed from nonEquilibrium PDF which in turn is computed from pre-collision PDF

  • private subroutine proc_calc_nNwtn_visc_fromPreColPDF(nNwtn, viscKine, omega, state, neigh, auxField, densPos, velPos, nSize, nSolve, nScalars, nAuxScalars, layout, convFac)

    Arguments

    Type IntentOptional Attributes Name
    class(mus_nNwtn_type), intent(in) :: nNwtn

    contains non-Newtonian model parameters loaded from config file
    real(kind=rk), intent(inout) :: viscKine(:)

    output is kinematic viscosity from non-Netonian model
    real(kind=rk), intent(in) :: omega(:)

    Kinematic viscosity omega from last timestep
    real(kind=rk), intent(in) :: state(:)

    state array
    integer, intent(in) :: neigh(:)

    neigh array to obtain precollision pdf
    real(kind=rk), intent(in) :: auxField(:)

    Auxiliary field variable array
    integer, intent(in) :: densPos

    position of density in auxField
    integer, intent(in) :: velPos(3)

    position of velocity components in auxField
    integer, intent(in) :: nSize

    number of elements in state array
    integer, intent(in) :: nSolve

    Number of element to solve in this level
    integer, intent(in) :: nScalars

    number of scalars in state array
    integer, intent(in) :: nAuxScalars

    number of scalars in auxField array
    type(mus_scheme_layout_type), intent(in) :: layout

    scheme layout
    type(mus_convertFac_type), intent(in) :: convFac

    conversion factor to convert lattice to physical units

Derived Types

type, public ::  mus_nNwtn_type

The nonNewtonian fluid feature description type

Read more…

Components

Type Visibility Attributes Name Initial
logical, public :: active = .false.

Indicator whether nonNewtonian feature is active maybe not useful. schemeHeader%kind can used to check if nNwtn is active
character(len=labellen), public :: label

nonNewtonian fluid model label
integer, public :: model

nonNewtonian fluid model identifier
type(mus_nNwtn_PL_type), public :: PL

Power Law (PL) model type
type(mus_nNwtn_CY_type), public :: CY

Carreau-Yasuda (CY) model type
type(mus_nNwtn_CS_type), public :: CS

Casson model type
procedure(proc_calc_nNwtn_visc_fromPreColPDF), public, pointer :: calcVisc => null()

this procedure compute kinematic viscosity in lattice unit on current level from preCollision PDF based on non-Newtonian model. It uses shear-rate to compute viscosity. Non-newtonian model is given in dynamic viscosity in physical unit so it is dimensionalized using viscDyna in physics%fac and lattice kinematic viscosity = lattice dynamic viscosity / rho (local density) for compressible model and lattice kinematic viscosity = lattice dynamic viscosity / rho0 for incompressible model.

type, private ::  mus_nNwtn_PL_type

The nonNewtonian power law model parameter

Read more…

Components

Type Visibility Attributes Name Initial
real(kind=rk), public :: n = 0.5_rk

exponentiation parameter
real(kind=rk), public :: visc0 = 0.0035_rk

Unit consistency index. Dynamic viscosity parameter when shear rate equals to 1. i.e. backgroud viscosity for power-law model
real(kind=rk), public :: nMinus1

parameter for computation

type, private ::  mus_nNwtn_CY_type

The nonNewtonian power law model parameter

Read more…

Components

Type Visibility Attributes Name Initial
real(kind=rk), public :: n = 0.2128_rk

model parameter
real(kind=rk), public :: a = 0.64_rk

model parameter
real(kind=rk), public :: lambda = 8.2_rk

model parameter
real(kind=rk), public :: visc0 = 0.16_rk

model parameter, dynamic viscosity at zero shear-rate
real(kind=rk), public :: viscInf = 0.0035_rk

model parameter, dynamic viscosity in infinity shear-rate
real(kind=rk), public :: nMinus1Div_a = (0.2128_rk-1._rk)/0.64_rk

calculated parameter for later usage, nMinus1Div_a = (n-1)/a

type, private ::  mus_nNwtn_CS_type

The nonNewtonian power law model parameter

Read more…

Components

Type Visibility Attributes Name Initial
real(kind=rk), public :: k0 = 0.1937_rk

model parameter
real(kind=rk), public :: k1 = 0.055_rk

model parameter

Functions

private function viscPhy_PL(me, shearRate)

nonNewtonian power-law model

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_type), intent(in) :: me

nonNewtonian parameters
real(kind=rk), intent(in) :: shearRate

Return Value real(kind=rk)

private function viscPhy_CS(me, shearRate)

nonNewtonian Casson model

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_type), intent(in) :: me

nonNewtonian parameters
real(kind=rk), intent(in) :: shearRate

Return Value real(kind=rk)

private function viscPhy_CY(me, shearRate)

nonNewtonian Carreau-Yasuda model

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_type), intent(in) :: me

nonNewtonian parameters
real(kind=rk), intent(in) :: shearRate

Return Value real(kind=rk)


Subroutines

public subroutine mus_nNwtn_load(me, conf, parent)

Read in the nonNewtonian table from Lua file and dump parameters to logUnit Specificly, this routine calls each model parameter loader.

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_type), intent(out) :: me

nonNewtonian type
type(flu_State), intent(inout) :: conf

lua state
integer, intent(in), optional :: parent

parent handle

public subroutine mus_nNwtn_save2lua(me, conf)

write nonNewtonian fluid parameters into a lua file

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_type), intent(in) :: me

nonNewtonian parameters
type(aot_out_type) :: conf

public subroutine mus_nNwtn_dump2outUnit(me, outUnit)

Dump nonNewtonian fluid parameters to outUnit

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_type), intent(in) :: me

nonNewtonian parameters
integer, intent(in) :: outUnit

public subroutine mus_assign_nNwtnVisc_ptr(nNwtn, schemeHeader)

This routine assigns function pointer to compute non-Newtonian viscosity

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_type), intent(inout) :: nNwtn

non-Newtonian type
type(mus_scheme_header_type), intent(in) :: schemeHeader

scheme header

public subroutine calcVisc_CY(nNwtn, viscKine, omega, state, neigh, auxField, densPos, velPos, nSize, nSolve, nScalars, nAuxScalars, layout, convFac)

Calculate kinematic viscosity from nonNewtonian Carreau-Yasuda model. $\mu = \mu_\inf + (\mu_0-\mu_\inf)(1+(\lambdashearRate)a)^((n-1)/a)$

Read more…

Arguments

Type IntentOptional Attributes Name
class(mus_nNwtn_type), intent(in) :: nNwtn

contains non-Newtonian model parameters loaded from config file
real(kind=rk), intent(inout) :: viscKine(:)

output is physical kinematic viscosity will be overwritten by non-Netonian model
real(kind=rk), intent(in) :: omega(:)

Kinematic viscosity omega from last timestep
real(kind=rk), intent(in) :: state(:)

state array
integer, intent(in) :: neigh(:)

neigh array to obtain precollision pdf
real(kind=rk), intent(in) :: auxField(:)

Auxiliary field variable array
integer, intent(in) :: densPos

position of density in auxField
integer, intent(in) :: velPos(3)

position of velocity components in auxField
integer, intent(in) :: nSize

number of elements in state array
integer, intent(in) :: nSolve

Number of element to solve in this level
integer, intent(in) :: nScalars

number of scalars in state array
integer, intent(in) :: nAuxScalars

number of scalars in auxField array
type(mus_scheme_layout_type), intent(in) :: layout

scheme layout
type(mus_convertFac_type), intent(in) :: convFac

conversion factor to convert lattice to physical units

private subroutine mus_nNwtn_PL_load(me, conf, tHandle)

Read in the nonNewtonian Power Law (PL) model parameters from Lua file

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_PL_type), intent(out) :: me

nonNewtonian type
type(flu_State), intent(inout) :: conf

lua state
integer, intent(inout) :: tHandle

nonNewtonian table handle

private subroutine mus_nNwtn_CY_load(me, conf, tHandle)

Read in the nonNewtonian Carreau-Yasuda (CY) model parameters from Lua file

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_CY_type), intent(out) :: me

nonNewtonian type
type(flu_State), intent(inout) :: conf

lua state
integer, intent(inout) :: tHandle

nonNewtonian table handle

private subroutine mus_nNwtn_CS_load(me, conf, tHandle)

Read in the nonNewtonian Casson model parameters from Lua file

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_CS_type), intent(out) :: me

nonNewtonian type
type(flu_State), intent(inout) :: conf

lua state
integer, intent(inout) :: tHandle

nonNewtonian table handle

private subroutine mus_nNwtn_PL_save(me, conf)

write nonNewtonian Power Law (PL) parameters into a lua file

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_PL_type), intent(in) :: me

nonNewtonian parameters
type(aot_out_type) :: conf

private subroutine mus_nNwtn_CY_save(me, conf)

write nonNewtonian (CY) parameters into a lua file

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_CY_type), intent(in) :: me

nonNewtonian parameters
type(aot_out_type) :: conf

private subroutine mus_nNwtn_CS_save(me, conf)

write nonNewtonian Casson parameters into a lua file

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_CS_type), intent(in) :: me

nonNewtonian parameters
type(aot_out_type) :: conf

private subroutine mus_nNwtn_PL_dump(me, outUnit)

Dump nonNewtonian Power Law (PL) parameters to outUnit

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_PL_type), intent(in) :: me

nonNewtonian parameters
integer, intent(in) :: outUnit

private subroutine mus_nNwtn_CY_dump(me, outUnit)

Dump nonNewtonian (CY) parameters to outUnit

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_CY_type), intent(in) :: me

nonNewtonian parameters
integer, intent(in) :: outUnit

private subroutine mus_nNwtn_CS_dump(me, outUnit)

Dump nonNewtonian (CY) parameters to outUnit

Arguments

Type IntentOptional Attributes Name
type(mus_nNwtn_CS_type), intent(in) :: me

nonNewtonian parameters
integer, intent(in) :: outUnit

private subroutine calcVisc_PL(nNwtn, viscKine, omega, state, neigh, auxField, densPos, velPos, nSize, nSolve, nScalars, nAuxScalars, layout, convFac)

Calculate kinematic viscosity from nonNewtonian power-law model. $\mu = K shearRate^(n-1)$. Shear rate is computed from strain rate which is computed from nonEquilibrium PDF which in turn computed from pre-collision PDF

Arguments

Type IntentOptional Attributes Name
class(mus_nNwtn_type), intent(in) :: nNwtn

contains non-Newtonian model parameters loaded from config file
real(kind=rk), intent(inout) :: viscKine(:)

output is physical kinematic viscosity will be overwritten by non-Netonian model
real(kind=rk), intent(in) :: omega(:)

Kinematic viscosity omega from last timestep
real(kind=rk), intent(in) :: state(:)

state array
integer, intent(in) :: neigh(:)

neigh array to obtain precollision pdf
real(kind=rk), intent(in) :: auxField(:)

Auxiliary field variable array
integer, intent(in) :: densPos

position of density in auxField
integer, intent(in) :: velPos(3)

position of velocity components in auxField
integer, intent(in) :: nSize

number of elements in state array
integer, intent(in) :: nSolve

Number of element to solve in this level
integer, intent(in) :: nScalars

number of scalars in state array
integer, intent(in) :: nAuxScalars

number of scalars in auxField array
type(mus_scheme_layout_type), intent(in) :: layout

scheme layout
type(mus_convertFac_type), intent(in) :: convFac

conversion factor to convert lattice to physical units

private subroutine calcVisc_CS(nNwtn, viscKine, omega, state, neigh, auxField, densPos, velPos, nSize, nSolve, nScalars, nAuxScalars, layout, convFac)

Calculate kinematic viscosity from nonNewtonian Casson model. $\mu = (k0 + k1 * sqrt(shearRate))^2/shearRate$

Read more…

Arguments

Type IntentOptional Attributes Name
class(mus_nNwtn_type), intent(in) :: nNwtn

contains non-Newtonian model parameters loaded from config file
real(kind=rk), intent(inout) :: viscKine(:)

output is physical kinematic viscosity will be overwritten by non-Netonian model
real(kind=rk), intent(in) :: omega(:)

Kinematic viscosity omega from last timestep
real(kind=rk), intent(in) :: state(:)

state array
integer, intent(in) :: neigh(:)

neigh array to obtain precollision pdf
real(kind=rk), intent(in) :: auxField(:)

Auxiliary field variable array
integer, intent(in) :: densPos

position of density in auxField
integer, intent(in) :: velPos(3)

position of velocity components in auxField
integer, intent(in) :: nSize

number of elements in state array
integer, intent(in) :: nSolve

Number of element to solve in this level
integer, intent(in) :: nScalars

number of scalars in state array
integer, intent(in) :: nAuxScalars

number of scalars in auxField array
type(mus_scheme_layout_type), intent(in) :: layout

scheme layout
type(mus_convertFac_type), intent(in) :: convFac

conversion factor to convert lattice to physical units

private subroutine calcVisc_incomp_PL(nNwtn, viscKine, omega, state, neigh, auxField, densPos, velPos, nSize, nSolve, nScalars, nAuxScalars, layout, convFac)

Calculate kinematic viscosity from nonNewtonian power-law model for incompressible model $\mu = K shearRate^(n-1)$

Read more…

Arguments

Type IntentOptional Attributes Name
class(mus_nNwtn_type), intent(in) :: nNwtn

contains non-Newtonian model parameters loaded from config file
real(kind=rk), intent(inout) :: viscKine(:)

output is physical kinematic viscosity will be overwritten by non-Netonian model
real(kind=rk), intent(in) :: omega(:)

Kinematic viscosity omega from last timestep
real(kind=rk), intent(in) :: state(:)

state array
integer, intent(in) :: neigh(:)

neigh array to obtain precollision pdf
real(kind=rk), intent(in) :: auxField(:)

Auxiliary field variable array
integer, intent(in) :: densPos

position of density in auxField
integer, intent(in) :: velPos(3)

position of velocity components in auxField
integer, intent(in) :: nSize

number of elements in state array
integer, intent(in) :: nSolve

Number of element to solve in this level
integer, intent(in) :: nScalars

number of scalars in state array
integer, intent(in) :: nAuxScalars

number of scalars in auxField array
type(mus_scheme_layout_type), intent(in) :: layout

scheme layout
type(mus_convertFac_type), intent(in) :: convFac

conversion factor to convert lattice to physical units

private subroutine calcVisc_incomp_CS(nNwtn, viscKine, omega, state, neigh, auxField, densPos, velPos, nSize, nSolve, nScalars, nAuxScalars, layout, convFac)

Calculate kinematic viscosity from nonNewtonian Casson model for incompressible model. $\mu = (k0 + k1 * sqrt(shearRate))^2/shearRate$

Read more…

Arguments

Type IntentOptional Attributes Name
class(mus_nNwtn_type), intent(in) :: nNwtn

contains non-Newtonian model parameters loaded from config file
real(kind=rk), intent(inout) :: viscKine(:)

output is physical kinematic viscosity will be overwritten by non-Netonian model
real(kind=rk), intent(in) :: omega(:)

Kinematic viscosity omega from last timestep
real(kind=rk), intent(in) :: state(:)

state array
integer, intent(in) :: neigh(:)

neigh array to obtain precollision pdf
real(kind=rk), intent(in) :: auxField(:)

Auxiliary field variable array
integer, intent(in) :: densPos

position of density in auxField
integer, intent(in) :: velPos(3)

position of velocity components in auxField
integer, intent(in) :: nSize

number of elements in state array
integer, intent(in) :: nSolve

Number of element to solve in this level
integer, intent(in) :: nScalars

number of scalars in state array
integer, intent(in) :: nAuxScalars

number of scalars in auxField array
type(mus_scheme_layout_type), intent(in) :: layout

scheme layout
type(mus_convertFac_type), intent(in) :: convFac

conversion factor to convert lattice to physical units

private subroutine calcVisc_incomp_CY(nNwtn, viscKine, omega, state, neigh, auxField, densPos, velPos, nSize, nSolve, nScalars, nAuxScalars, layout, convFac)

Calculate kinematic viscosity from nonNewtonian Carreau-Yasuda model for incompressible model. $\mu = \mu_\inf + (\mu_0-\mu_\inf)(1+(\lambdashearRate)a)^((n-1)/a)$

Read more…

Arguments

Type IntentOptional Attributes Name
class(mus_nNwtn_type), intent(in) :: nNwtn

contains non-Newtonian model parameters loaded from config file
real(kind=rk), intent(inout) :: viscKine(:)

output is physical kinematic viscosity will be overwritten by non-Netonian model
real(kind=rk), intent(in) :: omega(:)

Kinematic viscosity omega from last timestep
real(kind=rk), intent(in) :: state(:)

state array
integer, intent(in) :: neigh(:)

neigh array to obtain precollision pdf
real(kind=rk), intent(in) :: auxField(:)

Auxiliary field variable array
integer, intent(in) :: densPos

position of density in auxField
integer, intent(in) :: velPos(3)

position of velocity components in auxField
integer, intent(in) :: nSize

number of elements in state array
integer, intent(in) :: nSolve

Number of element to solve in this level
integer, intent(in) :: nScalars

number of scalars in state array
integer, intent(in) :: nAuxScalars

number of scalars in auxField array
type(mus_scheme_layout_type), intent(in) :: layout

scheme layout
type(mus_convertFac_type), intent(in) :: convFac

conversion factor to convert lattice to physical units