In File carreauDynamicsTemplates.h from complexDynamics directory,
the recursive template to compute the effective viscoisty
template<typename T,int N>
struct ImplicitOmega
{
T operator()(T alpha, T nu0_nuInfOverCs2, T nuInfOverCs2, T nMinusOneOverTwo, T omega0)
{
ImplicitOmega<T,N-1> omega;
T omSqr = omega(alpha,nu0_nuInfOverCs2,nuInfOverCs2,nMinusOneOverTwo,omega0);
omSqr *= omSqr;
return (T)2/((T)1+(T)2*nu0_nuInfOverCs2*
pow((T)1+alpha*omSqr,nMinusOneOverTwo)+2 * nuInfOverCs2);
}
};
template<typename T>
struct ImplicitOmega<T,0>
{
T operator()(T alpha, T nu0_nuInfOverCs2, T nuInfOverCs2, T nMinusOneOverTwo, T omega0)
{
T omegaTmp = (T)2/((T)1+(T)2*nu0_nuInfOverCs2*pow((T)1+
alpha*omega0*omega0,nMinusOneOverTwo)+2 * nuInfOverCs2);
}
};
/// This structure forwards the calls to the appropriate helper class
template<typename T, int N>
struct carreauDynamicsTemplates {
static T fromPiAndRhoToOmega(T alpha, T nu0_nuInfOverCs2, T nuInfOverCs2, T nMinusOneOverTwo, T omega0)
{
ImplicitOmega<T,N> omega;
return omega(alpha,nu0_nuInfOverCs2,nuInfOverCs2,nMinusOneOverTwo,omega0);
}
}; // struct carreauDynamicsTemplatese]
omega is base on the Eq. 5.9 from Dr. Orestis Malaspinas’ Thesis (there are a few typos in Eq5.5 5.8 5.9, hopefully I am not wrong.). There should be a factor 2 before nuInfOverCs2. Certainly, If you don’t set nuInf in Carreau model (the default value is 0.) Then it won’t be an issue.
Plus, Even I set the exponent n = 1 in the code carreauPoiseuille.cpp, it doesn’t reproduce the poiseuille parabolic profile. (This example doesn’t come with velocity profile for the corresponding non-Newtonian fluid for setting-up the inlet boundary condition. I therefore use GuoExternalForceBGKdynamics rather than the default BGKdynamics.). Therefore, other issues might exist in the carreauDynamics.h and carreauDynamics.hh as well.