Add Carreau–Yasuda non-Newtonian viscosity model

docs/src/refs.bib

@@ -195,6 +195,18 @@ @TechReport{Carpenter1994
   number      = {NASA-TM-109112},
   month       = jun,
 }
+@Article{Carreau1972,
+  author    = {Carreau, Pierre J.},
+  title     = {Rheological Equations from Molecular Network Theories},
+  journal   = {Transactions of the Society of Rheology},
+  year      = {1972},
+  volume    = {16},
+  number    = {1},
+  pages     = {99--127},
+  issn      = {0038-0032},
+  doi       = {10.1122/1.549276},
+  publisher = {Wiley},
+}
 
 @Article{Clausen2013,
   author    = {Clausen, Jonathan R.},
@@ -810,6 +822,18 @@ @Article{Tafuni2018
   publisher = {Elsevier {BV}},
 }
 
+@Article{VahabiSadeghy2014,
+  author    = {Vahabi, M. and Sadeghy, K.},
+  title     = {Simulating Bubble Shape during its Rise in Carreau--Yasuda Fluids Using {WC-SPH} Method},
+  journal   = {Nihon Reoroji Gakkaishi},
+  year      = {2014},
+  volume    = {41},
+  number    = {5},
+  pages     = {319--329},
+  month     = jan,
+  doi       = {10.1678/rheology.41.319},
+}
+
 @Article{Valizadeh2015,
   author    = {Valizadeh, Alireza and Monaghan, Joseph J.},
   title     = {A study of solid wall models for weakly compressible SPH},
@@ -862,6 +886,33 @@ @Article{Westerhof2008
   publisher = {Springer Science and Business Media LLC},
 }
 
+@Article{Yasuda1981,
+  author    = {Yasuda, K. and Armstrong, R. C. and Cohen, R. E.},
+  title     = {Shear flow properties of concentrated solutions of linear and star branched polystyrenes},
+  journal   = {Rheologica Acta},
+  year      = {1981},
+  volume    = {20},
+  number    = {2},
+  pages     = {163--178},
+  issn      = {0035-4511},
+  doi       = {10.1007/BF01513059},
+  publisher = {Springer},
+}
+
+@Article{Zhang2017,
+  author    = {Zhang, Li and Ban, Xiaojuan and Wang, Xiaokun and Liu, Xing},
+  title     = {A Symmetry Particle Method towards Implicit Non-Newtonian Fluids},
+  journal   = {Symmetry},
+  year      = {2017},
+  volume    = {9},
+  number    = {2},
+  pages     = {26},
+  month     = feb,
+  issn      = {2073-8994},
+  doi       = {10.3390/sym9020026},
+  publisher = {MDPI AG},
+}
+
 @Article{Zhang2025,
   author    = {Zhang, Shuoguo and Fan, Yu and Wu, Dong and Zhang, Chi and Hu, Xiangyu},
   journal   = {Physics of Fluids},

docs/src/systems/fluid.md

@@ -155,6 +155,47 @@ The inter-particle-averaged shear stress is defined as:
 
 with the dynamic viscosity of each particle given by `` \eta_a = \rho_a \nu_a ``, where `` \nu_a `` is the kinematic viscosity.
 
+#### ViscosityCarreauYasuda
+
+`ViscosityCarreauYasuda` implements the Carreau–Yasuda non-Newtonian viscosity model,
+originally proposed by [Carreau (1972)](@cite Carreau1972) and extended by
+[Yasuda et al. (1981)](@cite Yasuda1981). In this model, the kinematic viscosity
+depends on the local shear rate. This makes it suitable for shear-thinning and
+shear-thickening fluids, such as polymer solutions or blood-like fluids.
+Instead of prescribing a single constant viscosity, the apparent viscosity
+smoothly transitions between a low-shear plateau and a high-shear plateau.
+
+In SPH, this can be incorporated by evaluating a shear-rate-dependent
+viscosity locally and using it in the standard viscous discretization. A Newtonian
+fluid is recovered as a special case when the parameters are chosen such that the
+viscosity becomes independent of the shear rate. ([Zhang et al. (2017)](@cite Zhang2017);
+[Vahabi & Sadeghy (2014)](@cite VahabiSadeghy2014)).
+
+
+##### Mathematical Formulation
+
+In the Carreau–Yasuda model, the kinematic viscosity ``\nu`` depends on the shear-rate magnitude ``\dot\gamma`` as
+```math
+\nu(\dot\gamma) = \nu_\infty + (\nu_0 - \nu_\infty)
+\left[ 1 + (\lambda \dot\gamma)^a \right]^{\frac{n-1}{a}}.
+```
+where
+
+- ``\nu_0``: zero-shear kinematic viscosity,
+- ``\nu_\infty``: infinite-shear kinematic viscosity,
+- ``\lambda``: time constant,
+- ``a``: Yasuda parameter,
+- ``n``: power-law index (``n < 1`` for shear-thinning, ``n > 1`` for shear-thickening),
+- ``\dot\gamma``: shear-rate magnitude.
+
+In this implementation the shear-rate magnitude is approximated per particle pair as
+``\dot\gamma \approx \frac{\lVert \mathbf{v}_{ab} \rVert}{\lVert \mathbf{r}_{ab} \rVert + \epsilon}``,
+with ``\mathbf{v}_{ab}`` the relative velocity, ``\mathbf{r}_{ab}`` the position difference,
+and ``\epsilon`` a small regularization parameter.
+
+All viscosities here are kinematic viscosities (m²/s); dynamic viscosity is obtained internally
+via ``\eta = \rho \nu``. A Newtonian fluid is recovered for ``n = 1`` and
+``\nu_0 = \nu_\infty``
 
 ```@autodocs
 Modules = [TrixiParticles]

src/TrixiParticles.jl

@@ -81,7 +81,7 @@ export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        WendlandC6Kernel, SpikyKernel, Poly6Kernel, LaguerreGaussKernel
 export StateEquationCole, StateEquationIdealGas, StateEquationAdaptiveCole
 export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS,
-       ViscosityMorrisSGS
+       ViscosityMorrisSGS, ViscosityCarreauYasuda
 export DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari, DensityDiffusionAntuono
 export tensile_instability_control
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,

src/io/io.jl

@@ -341,3 +341,14 @@ function add_system_data!(system_data, shifting_technique::ParticleShiftingTechn
     system_data["shifting_technique"] = Dict{String, Any}()
     system_data["shifting_technique"]["model"] = type2string(shifting_technique)
 end
+
+function add_system_data!(system_data, viscosity::ViscosityCarreauYasuda)
+    system_data["viscosity_model"] = Dict{String, Any}()
+    system_data["viscosity_model"]["model"]   = type2string(viscosity)
+    system_data["viscosity_model"]["nu0"]     = viscosity.nu0
+    system_data["viscosity_model"]["nu_inf"]  = viscosity.nu_inf
+    system_data["viscosity_model"]["lambda"]  = viscosity.lambda
+    system_data["viscosity_model"]["a"]       = viscosity.a
+    system_data["viscosity_model"]["n"]       = viscosity.n
+    system_data["viscosity_model"]["epsilon"] = viscosity.epsilon
+end

src/schemes/fluid/viscosity.jl

@@ -447,3 +447,68 @@ function kinematic_viscosity(system, viscosity::ViscosityMorrisSGS, smoothing_le
                              sound_speed)
     return viscosity.nu
 end
+
+@doc raw"""
+    ViscosityCarreauYasuda(; nu0, nu_inf, lambda, a, n, epsilon=0.01)
+
+Non-Newtonian viscosity model based on the Carreau–Yasuda law [Carreau (1972)](@cite Carreau1972), [Yasuda et al. (1981)](@cite Yasuda1981).
+
+See [the docs on viscosity](@ref viscosity_sph) for an overview and comparison of implemented viscosity models.
+
+# Keywords
+- `nu0`:     Zero-shear kinematic viscosity.
+- `nu_inf`:  Infinite-shear kinematic viscosity.
+- `lambda`:  Time constant of the Carreau–Yasuda law.
+- `a`:       Yasuda parameter controlling the transition shape.
+- `n`:       Power-law index (shear-thinning/thickening behavior).
+- `epsilon`: Parameter to prevent singularities in the shear-rate approximation.
+"""
+struct ViscosityCarreauYasuda{ELTYPE}
+    nu0     :: ELTYPE  # zero-shear kinematic viscosity
+    nu_inf  :: ELTYPE  # infinite-shear kinematic viscosity
+    lambda  :: ELTYPE  # time constant
+    a       :: ELTYPE  # Yasuda parameter
+    n       :: ELTYPE  # power-law index
+    epsilon :: ELTYPE  # regularization
+end
+
+function ViscosityCarreauYasuda(; nu0, nu_inf, lambda, a, n, epsilon=0.01)
+    ViscosityCarreauYasuda(nu0, nu_inf, lambda, a, n, epsilon)
+end
+
+@propagate_inbounds function (viscosity::ViscosityCarreauYasuda)(particle_system,
+                                                                 neighbor_system,
+                                                                 v_particle_system,
+                                                                 v_neighbor_system,
+                                                                 particle, neighbor,
+                                                                 pos_diff, distance,
+                                                                 sound_speed,
+                                                                 m_a, m_b,
+                                                                 rho_a, rho_b,
+                                                                 grad_kernel)
+    epsilon = viscosity.epsilon
+
+    smoothing_length_particle = smoothing_length(particle_system, particle)
+    smoothing_length_neighbor = smoothing_length(particle_system, neighbor)
+    smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
+
+    v_a = viscous_velocity(v_particle_system, particle_system, particle)
+    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
+    v_diff = v_a - v_b
+
+    gamma_dot = norm(v_diff) / (distance + epsilon)
+
+    # Compute Carreau-Yasuda effective viscosity
+    (; nu0, nu_inf, lambda, a, n) = viscosity
+    nu_eff = nu_inf + (nu0 - nu_inf) * (1 + (lambda * gamma_dot)^a)^((n - 1) / a)
+    nu_a = nu_eff
+    nu_b = nu_eff
+
+    return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel,
+                                 m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon)
+end
+
+function kinematic_viscosity(system, viscosity::ViscosityCarreauYasuda,
+                             smoothing_length, sound_speed)
+    return viscosity.nu0
+end

test/schemes/fluid/viscosity.jl

@@ -135,4 +135,59 @@
         @test isapprox(dv[1], -2.0328356978041036e-5, atol=6e-15)
         @test isapprox(dv[2], 6.776118992680346e-5, atol=6e-15)
     end
+
+    @testset verbose=true "`ViscosityCarreauYasuda`" begin
+        v = fluid.velocity
+        v .= 0.0
+        v[1, 1] = v_diff[1]
+        v[2, 1] = v_diff[2]
+
+        m_a = 0.01
+        m_b = 0.01
+
+        # ------------------------------------------------------------------
+        # 1) Newtonian limit: should match the constant-viscosity behavior
+        # ------------------------------------------------------------------
+        nu = 7e-3
+        viscosity = ViscosityCarreauYasuda(nu0     = nu, 
+                                           nu_inf  = nu,
+                                           lambda  = 1.0,
+                                           a       = 2.0,
+                                           n       = 1.0,
+                                           epsilon = 0.01)
+        system_wcsph = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(),
+                                                   state_equation, smoothing_kernel,
+                                                   smoothing_length; viscosity=viscosity)
+
+        grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff, 
+                                                           distance, 1)
+        dv = viscosity(system_wcsph, system_wcsph,
+                       v, v, 1, 2, pos_diff, distance,
+                       sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+
+        @test isapprox(dv[1], -1.0895602048035410e-5, atol=6e-15)
+        @test isapprox(dv[2],  3.6318673493451364e-5, atol=6e-15)
+
+        # ------------------------------------------------------------------
+        # 2) Shear-thinning case: fixed (precomputed) values
+        # ------------------------------------------------------------------
+        viscosity = ViscosityCarreauYasuda(nu0     = 3.5e-6,
+                                           nu_inf  = 1.0e-6,
+                                           lambda  = 3.313e-2,
+                                           a       = 2.0,
+                                           n       = 0.3,
+                                           epsilon = 0.01)
+        system_wcsph = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(),
+                                                   state_equation, smoothing_kernel,
+                                                   smoothing_length; viscosity=viscosity)
+
+        grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff, 
+                                                           distance, 1)
+        dv = viscosity(system_wcsph, system_wcsph,
+                       v, v, 1, 2, pos_diff, distance,
+                       sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+
+        @test isapprox(dv[1], -5.33743497379846e-9, atol=6e-15)
+        @test isapprox(dv[2],  1.7791449912661534e-8, atol=6e-15)
+    end
 end