Problem
svMultiPhysics currently supports three fluid viscosity models: Newtonian, Carreau-Yasuda, and Casson (get_viscosity() in Code/Source/solver/fluid.cpp). It does not support the Quemada model, which accounts for both hematocrit and shear-rate dependence of blood viscosity. Quemada is used in cardiovascular CFD where blood rheology may affect flow patterns — e.g., hepatic flow distribution in Fontan (TCPC) simulations — and we need it to reproduce/compare against Fontan simulations run in other solvers.
In addition, the existing non-Newtonian test cases (tests/cases/fluid/carreau_yasuda, tests/cases/fluid/casson) are regression tests: they compare against a stored result_001.vtu from a pulsatile pipe run, but do not verify the models against a known solution. So a model that is coded incorrectly but stable would still pass.
Proposed solution
1. Add the Quemada model
Quemada viscosity [1, 2]:
$$\mu(\dot\gamma) = \mu_F \left(1 - \frac{k \phi}{2}\right)^{-2}$$
where the intrinsic viscosity $k$ depends on shear rate:
$$k(\dot\gamma) = \frac{k_0 + k_\infty \sqrt{\dot\gamma/\dot\gamma_c}}{1 + \sqrt{\dot\gamma/\dot\gamma_c}}$$
with plasma viscosity $\mu_F$, hematocrit $\phi$, and parameters $k_0$ (low-shear limit of $k$), $k_\infty$ (high-shear limit of $k$), and critical shear rate $\dot\gamma_c$. The last three typically depend on hematocrit.
[1] Quemada D. "Rheology of concentrated disperse systems and minimum energy dissipation principle. I. Viscosity-concentration relationship." Rheol Acta 1977;16(1):82–94. doi:10.1007/BF01516932
[2] Quemada D. "Rheology of concentrated disperse systems II. A model for non-newtonian shear viscosity in steady flows." Rheol Acta 1978;17(6):632–642. doi:10.1007/BF01522036
- Add
viscType_Quemada to FluidViscosityModelType and the name map in consts.cpp
- Add parameter parsing in
read_files.cpp / Parameters
- Implement
mu, mu_s, and mu_x (derivative w.r.t. shear rate, for the tangent) in get_viscosity()
- Document parameters and add a README, following the existing Casson/CY test cases
2. Verification test cases for all non-Newtonian models
Add steady, fully developed flow in a straight tube driven by a constant pressure gradient, compared against analytical/reference solutions for velocity profile and flow rate:
-
Quemada: Popel AS, Enden G. "An analytical solution for steady flow of a Quemada fluid in a circular tube." Rheol Acta 1993;32(4):422–426. doi:10.1007/BF00435088 (PMC). Closed-form velocity profile and flow rate in terms of $k_0$, $k_\infty$, $\dot\gamma_c$. The solution reduces to Newtonian (parabolic) and Casson in limiting cases, so one setup gives three checks.
-
Casson: closed-form plug-flow solution for a Casson fluid in a tube.
-
Carreau-Yasuda: no closed form; compare against a high-resolution 1D numerical solution of the radial momentum equation.
-
Newtonian: Poiseuille.
Suggested metrics: L2 / max error of the axial velocity profile along a radial line, and relative error in flow rate, with mesh refinement to show convergence.
This verifies that the constitutive models are coded correctly. It does not validate which blood model is physiologically appropriate; that is out of scope here.
3. Review of existing implementations (to be confirmed by the tests above)
While reading get_viscosity(), a few things in the Casson branch look worth checking:
mu_x is computed as 2*mu_o*(mu_o+mu_i)/gamma, which is positive, while Casson viscosity decreases with shear rate (d mu/d gamma < 0). Compare with the CY branch, where mu_x has the sign of (n-1).- Below the low-shear-rate threshold, viscosity is held constant but
mu_x is non-zero.
gamma is passed by reference and overwritten with lam below the threshold, which changes the caller's value.
These may affect only nonlinear convergence (tangent) rather than the converged solution, but the verification cases should make that clear.
Additional context
Related: #633 (templated viscous stress models) touches the same code area.
Code of Conduct
Problem
svMultiPhysics currently supports three fluid viscosity models: Newtonian, Carreau-Yasuda, and Casson (
get_viscosity()inCode/Source/solver/fluid.cpp). It does not support the Quemada model, which accounts for both hematocrit and shear-rate dependence of blood viscosity. Quemada is used in cardiovascular CFD where blood rheology may affect flow patterns — e.g., hepatic flow distribution in Fontan (TCPC) simulations — and we need it to reproduce/compare against Fontan simulations run in other solvers.In addition, the existing non-Newtonian test cases (
tests/cases/fluid/carreau_yasuda,tests/cases/fluid/casson) are regression tests: they compare against a storedresult_001.vtufrom a pulsatile pipe run, but do not verify the models against a known solution. So a model that is coded incorrectly but stable would still pass.Proposed solution
1. Add the Quemada model
Quemada viscosity [1, 2]:
where the intrinsic viscosity$k$ depends on shear rate:
with plasma viscosity$\mu_F$ , hematocrit $\phi$ , and parameters $k_0$ (low-shear limit of $k$ ), $k_\infty$ (high-shear limit of $k$ ), and critical shear rate $\dot\gamma_c$ . The last three typically depend on hematocrit.
[1] Quemada D. "Rheology of concentrated disperse systems and minimum energy dissipation principle. I. Viscosity-concentration relationship." Rheol Acta 1977;16(1):82–94. doi:10.1007/BF01516932
[2] Quemada D. "Rheology of concentrated disperse systems II. A model for non-newtonian shear viscosity in steady flows." Rheol Acta 1978;17(6):632–642. doi:10.1007/BF01522036
viscType_QuemadatoFluidViscosityModelTypeand the name map inconsts.cppread_files.cpp/Parametersmu,mu_s, andmu_x(derivative w.r.t. shear rate, for the tangent) inget_viscosity()2. Verification test cases for all non-Newtonian models
Add steady, fully developed flow in a straight tube driven by a constant pressure gradient, compared against analytical/reference solutions for velocity profile and flow rate:
Suggested metrics: L2 / max error of the axial velocity profile along a radial line, and relative error in flow rate, with mesh refinement to show convergence.
This verifies that the constitutive models are coded correctly. It does not validate which blood model is physiologically appropriate; that is out of scope here.
3. Review of existing implementations (to be confirmed by the tests above)
While reading
get_viscosity(), a few things in the Casson branch look worth checking:mu_xis computed as2*mu_o*(mu_o+mu_i)/gamma, which is positive, while Casson viscosity decreases with shear rate (d mu/d gamma < 0). Compare with the CY branch, wheremu_xhas the sign of(n-1).mu_xis non-zero.gammais passed by reference and overwritten withlambelow the threshold, which changes the caller's value.These may affect only nonlinear convergence (tangent) rather than the converged solution, but the verification cases should make that clear.
Additional context
Related: #633 (templated viscous stress models) touches the same code area.
Code of Conduct