Added Tangent_ShearViscosity_StrainRate member to all visous fluid materials.

SteveMaas1978 · SteveMaas1978 · commit 34df599c24e5 · 2026-03-06T07:00:42.000-07:00
diff --git a/FEBioFluid/FEBinghamFluid.cpp b/FEBioFluid/FEBinghamFluid.cpp
@@ -99,6 +99,18 @@ double FEBinghamFluid::ShearViscosity(FEMaterialPoint& pt)
     return mu;
 }
 
+double FEBinghamFluid::Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp)
+{
+	FEFluidMaterialPoint& vt = *mp.ExtractData<FEFluidMaterialPoint>();
+	mat3ds D = vt.RateOfDeformation();
+	double dmu = 0;
+	double gdot = sqrt(2 * (D.sqr()).tr());
+	if (gdot > 0) {
+		dmu = m_tauy / pow(gdot, 2) * ((1 + m_n * gdot) * exp(-m_n * gdot) - 1);
+	}
+	return dmu;
+}
+
 //-----------------------------------------------------------------------------
 //! bulk viscosity
 double FEBinghamFluid::BulkViscosity(FEMaterialPoint& pt)
diff --git a/FEBioFluid/FEBinghamFluid.h b/FEBioFluid/FEBinghamFluid.h
@@ -49,6 +49,9 @@ class FEBIOFLUID_API FEBinghamFluid : public FEViscousFluid
     
     //! dynamic viscosity
     double ShearViscosity(FEMaterialPoint& mp) override;
+
+	//! derivative of shear viscosity w.r.t. strain rate
+	double Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp) override;
     
     //! bulk viscosity
     double BulkViscosity(FEMaterialPoint& mp) override;
diff --git a/FEBioFluid/FECarreauFluid.cpp b/FEBioFluid/FECarreauFluid.cpp
@@ -95,6 +95,19 @@ double FECarreauFluid::ShearViscosity(FEMaterialPoint& pt)
     return mu;
 }
 
+//! derivative of shear viscosity w.r.t. strain rate
+double FECarreauFluid::Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp)
+{
+	FEFluidMaterialPoint& vt = *mp.ExtractData<FEFluidMaterialPoint>();
+	mat3ds D = vt.RateOfDeformation();
+	double gdot = sqrt(2 * (D.sqr()).tr());
+	double lamg2 = m_lam * m_lam * gdot * gdot;
+
+	double dmu = 2 * (m_mu0 - m_mui) * (m_n - 1) * m_lam * m_lam * pow(1 + lamg2, (m_n - 3) * 0.5);
+
+	return dmu;
+}
+
 //-----------------------------------------------------------------------------
 //! bulk viscosity
 double FECarreauFluid::BulkViscosity(FEMaterialPoint& pt)
diff --git a/FEBioFluid/FECarreauFluid.h b/FEBioFluid/FECarreauFluid.h
@@ -50,6 +50,9 @@ class FEBIOFLUID_API FECarreauFluid :	public FEViscousFluid
     //! dynamic viscosity
     double ShearViscosity(FEMaterialPoint& mp) override;
     
+	//! derivative of shear viscosity w.r.t. strain rate
+	double Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp) override;
+
     //! bulk viscosity
     double BulkViscosity(FEMaterialPoint& mp) override;
     
diff --git a/FEBioFluid/FECarreauYasudaFluid.cpp b/FEBioFluid/FECarreauYasudaFluid.cpp
@@ -99,6 +99,20 @@ double FECarreauYasudaFluid::ShearViscosity(FEMaterialPoint& pt)
     return mu;
 }
 
+//! derivative of shear viscosity w.r.t. strain rate
+double FECarreauYasudaFluid::Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp)
+{
+	FEFluidMaterialPoint& vt = *mp.ExtractData<FEFluidMaterialPoint>();
+	mat3ds D = vt.RateOfDeformation();
+	double gdot = sqrt(2 * (D.sqr()).tr());
+	double lamga = pow(m_lam * gdot, m_a);
+
+	double dmu = (m_a >= 2) ? 2 * (m_mu0 - m_mui) * (m_n - 1) * pow(m_lam, m_a) * pow(gdot, m_a - 2)
+		* pow(1 + lamga, (m_n - m_a - 1) / m_a) : 0;
+
+	return dmu;
+}
+
 //-----------------------------------------------------------------------------
 //! bulk viscosity
 double FECarreauYasudaFluid::BulkViscosity(FEMaterialPoint& pt)
diff --git a/FEBioFluid/FECarreauYasudaFluid.h b/FEBioFluid/FECarreauYasudaFluid.h
@@ -50,6 +50,9 @@ class FEBIOFLUID_API FECarreauYasudaFluid :	public FEViscousFluid
     //! dynamic viscosity
     double ShearViscosity(FEMaterialPoint& mp) override;
     
+	//! derivative of shear viscosity w.r.t. strain rate
+	double Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp) override;
+
     //! bulk viscosity
     double BulkViscosity(FEMaterialPoint& mp) override;
     
diff --git a/FEBioFluid/FECrossFluid.cpp b/FEBioFluid/FECrossFluid.cpp
@@ -96,6 +96,19 @@ double FECrossFluid::ShearViscosity(FEMaterialPoint& pt)
     return mu;
 }
 
+//! derivative of shear viscosity w.r.t. strain rate
+double FECrossFluid::Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp)
+{
+	FEFluidMaterialPoint& vt = *mp.ExtractData<FEFluidMaterialPoint>();
+	mat3ds D = vt.RateOfDeformation();
+	double gdot = sqrt(2 * (D.sqr()).tr());
+	double lamg = m_lam * gdot;
+
+	double dmu = -2 * (m_mu0 - m_mui) * m_m * pow(m_lam, m_m) * pow(gdot, m_m - 2) / pow(1 + pow(lamg, m_m), 2);
+
+	return dmu;
+}
+
 //-----------------------------------------------------------------------------
 //! bulk viscosity
 double FECrossFluid::BulkViscosity(FEMaterialPoint& pt)
diff --git a/FEBioFluid/FECrossFluid.h b/FEBioFluid/FECrossFluid.h
@@ -50,6 +50,9 @@ class FEBIOFLUID_API FECrossFluid :	public FEViscousFluid
     //! dynamic viscosity
     double ShearViscosity(FEMaterialPoint& mp) override;
     
+	//! derivative of shear viscosity w.r.t. strain rate
+	double Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp) override;
+
     //! bulk viscosity
     double BulkViscosity(FEMaterialPoint& mp) override;
     
diff --git a/FEBioFluid/FENewtonianFluid.cpp b/FEBioFluid/FENewtonianFluid.cpp
@@ -106,6 +106,12 @@ double FENewtonianFluid::ShearViscosity(FEMaterialPoint& mp)
     return m_mu;
 }
 
+//! derivative of shear viscosity w.r.t. strain rate
+double FENewtonianFluid::Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp)
+{
+	return 0.0;
+}
+
 //-----------------------------------------------------------------------------
 //! bulk viscosity
 double FENewtonianFluid::BulkViscosity(FEMaterialPoint& mp)
diff --git a/FEBioFluid/FENewtonianFluid.h b/FEBioFluid/FENewtonianFluid.h
@@ -57,6 +57,9 @@ class FEBIOFLUID_API FENewtonianFluid :	public FEViscousFluid
     //! dynamic viscosity
     double ShearViscosity(FEMaterialPoint& mp) override;
 
+	//! derivative of shear viscosity w.r.t. strain rate
+	double Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp) override;
+
     //! bulk viscosity
     double BulkViscosity(FEMaterialPoint& mp) override;
 
diff --git a/FEBioFluid/FEPowellEyringFluid.cpp b/FEBioFluid/FEPowellEyringFluid.cpp
@@ -96,6 +96,20 @@ double FEPowellEyringFluid::ShearViscosity(FEMaterialPoint& pt)
     return mu;
 }
 
+//! derivative of shear viscosity w.r.t. strain rate
+double FEPowellEyringFluid::Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp)
+{
+	FEFluidMaterialPoint& vt = *mp.ExtractData<FEFluidMaterialPoint>();
+	mat3ds D = vt.RateOfDeformation();
+	double gdot = sqrt(2 * (D.sqr()).tr());
+	double lamg = m_lam * gdot;
+
+	double dmu = (lamg < 1e-3) ? -2 * (m_mu0 - m_mui) * m_lam * m_lam / 3. :
+		(2 * (m_mu0 - m_mui) * (gdot / sqrt(1 + pow(lamg, 2)) - asinh(lamg) / m_lam)) / pow(gdot, 3);
+
+	return dmu;
+}
+
 //-----------------------------------------------------------------------------
 //! bulk viscosity
 double FEPowellEyringFluid::BulkViscosity(FEMaterialPoint& pt)
diff --git a/FEBioFluid/FEPowellEyringFluid.h b/FEBioFluid/FEPowellEyringFluid.h
@@ -50,6 +50,9 @@ class FEBIOFLUID_API FEPowellEyringFluid :	public FEViscousFluid
     //! dynamic viscosity
     double ShearViscosity(FEMaterialPoint& mp) override;
     
+	//! derivative of shear viscosity w.r.t. strain rate
+	double Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp) override;
+
     //! bulk viscosity
     double BulkViscosity(FEMaterialPoint& mp) override;
     
diff --git a/FEBioFluid/FEQuemadaFluid.cpp b/FEBioFluid/FEQuemadaFluid.cpp
@@ -99,6 +99,17 @@ double FEQuemadaFluid::ShearViscosity(FEMaterialPoint& pt)
     return mu;
 }
 
+//! derivative of shear viscosity w.r.t. strain rate
+double FEQuemadaFluid::Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp)
+{
+	FEFluidMaterialPoint& vt = *mp.ExtractData<FEFluidMaterialPoint>();
+	mat3ds D = vt.RateOfDeformation();
+	double gdot = sqrt(2 * (D.sqr()).tr());
+	double grsqrt = sqrt(gdot / m_gc);
+	double dmu = (gdot > 0) ? 4 * m_mu0 * m_H / m_gc * (m_k0 - m_ki) * (1 + grsqrt) / grsqrt / pow(-2 * (1 + grsqrt) + m_H * (m_k0 + m_ki * grsqrt), 3) : 0.0;
+	return dmu;
+}
+
 //-----------------------------------------------------------------------------
 //! bulk viscosity
 double FEQuemadaFluid::BulkViscosity(FEMaterialPoint& pt)
diff --git a/FEBioFluid/FEQuemadaFluid.h b/FEBioFluid/FEQuemadaFluid.h
@@ -51,6 +51,9 @@ class FEBIOFLUID_API FEQuemadaFluid :	public FEViscousFluid
     //! dynamic viscosity
     double ShearViscosity(FEMaterialPoint& mp) override;
     
+	//! derivative of shear viscosity w.r.t. strain rate
+	double Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp) override;
+
     //! bulk viscosity
     double BulkViscosity(FEMaterialPoint& mp) override;
     
diff --git a/FEBioFluid/FEViscousFluid.h b/FEBioFluid/FEViscousFluid.h
@@ -49,9 +49,12 @@ class FEBIOFLUID_API FEViscousFluid : public FEMaterialProperty
     
     //! tangent of stress with respect to rate of deformation tensor D
     virtual tens4ds Tangent_RateOfDeformation(FEMaterialPoint& mp) = 0;
-    
+
     //! dynamic viscosity
     virtual double ShearViscosity(FEMaterialPoint& mp) = 0;
+
+	//! derivative of shear viscosity w.r.t. strain rate
+	virtual double Tangent_ShearViscosity_StrainRate(FEMaterialPoint& mp) { return 0.0; }
     
     //! bulk viscosity
     virtual double BulkViscosity(FEMaterialPoint& mp) = 0;
