Added secant stress, secant tangent, and secant permeability tangent …

…calculations to biphasic solid and biphasic shell domains
febiosoftware · Dec 15, 2024 · 3db8697 · 3db8697
1 parent aa04afe
commit 3db8697
Show file tree

Hide file tree

Showing 10 changed files with 209 additions and 30 deletions.
diff --git a/FEBioMix/FEBioMixPlot.cpp b/FEBioMix/FEBioMixPlot.cpp
@@ -1523,7 +1523,7 @@ bool FEPlotEffectiveElasticity::Save(FEDomain &dom, FEDataStream& a)
 		for (int j=0; j<nint; ++j)
 		{
 			FEMaterialPoint& pt = *el.GetMaterialPoint(j);
-			if      (pb ) s += pb->Tangent(pt);
+            if      (pb ) s += (pb->Tangent(pt)).supersymm();
 			else if (pbs) s += pbs->Tangent(pt);
 			else if (ptp) s += ptp->Tangent(pt);
 			else if (pmp) s += pmp->Tangent(pt);

diff --git a/FEBioMix/FEBiphasic.cpp b/FEBioMix/FEBiphasic.cpp
@@ -177,17 +177,32 @@ mat3ds FEBiphasic::Stress(FEMaterialPoint& mp)
 	return s;
 }
 
+mat3ds FEBiphasic::SecantStress(FEMaterialPoint& mp)
+{
+    FEBiphasicMaterialPoint& pt = *mp.ExtractData<FEBiphasicMaterialPoint>();
+
+    // calculate solid material stress
+    mat3ds s = m_pSolid->SecantStress(mp);
+
+    // add fluid pressure
+    s.xx() -= pt.m_p;
+    s.yy() -= pt.m_p;
+    s.zz() -= pt.m_p;
+
+    return s;
+}
+
 //-----------------------------------------------------------------------------
 //! The tangent is the sum of the elastic tangent plus the fluid tangent. Note
 //! that this function is declared in the base class, so you don't have to 
 //! reimplement it unless additional tangent components are required.
 
-tens4ds FEBiphasic::Tangent(FEMaterialPoint& mp)
+tens4dmm FEBiphasic::Tangent(FEMaterialPoint& mp)
 {
 	FEBiphasicMaterialPoint& pt = *mp.ExtractData<FEBiphasicMaterialPoint>();
 
 	// call solid tangent routine
-	tens4ds c = m_pSolid->Tangent(mp);
+	tens4dmm c = m_pSolid->Tangent(mp);
 
 	// fluid pressure
 	double p = pt.m_p;
@@ -208,7 +223,41 @@ tens4ds FEBiphasic::Tangent(FEMaterialPoint& mp)
 	D[4][4] -= -p;
 	D[5][5] -= -p;
 
-	return tens4ds(D);
+	return tens4dmm(D);
+}
+
+//-----------------------------------------------------------------------------
+//! The tangent is the sum of the elastic tangent plus the fluid tangent. Note
+//! that this function is declared in the base class, so you don't have to
+//! reimplement it unless additional tangent components are required.
+
+tens4dmm FEBiphasic::SecantTangent(FEMaterialPoint& mp)
+{
+    FEBiphasicMaterialPoint& pt = *mp.ExtractData<FEBiphasicMaterialPoint>();
+
+    // call solid tangent routine
+    tens4dmm c = m_pSolid->SecantTangent(mp);
+
+    // fluid pressure
+    double p = pt.m_p;
+
+    // adjust tangent for pressures
+    double D[6][6] = {0};
+    c.extract(D);
+
+    D[0][0] -= -p;
+    D[1][1] -= -p;
+    D[2][2] -= -p;
+
+    D[0][1] -= p; D[1][0] -= p;
+    D[1][2] -= p; D[2][1] -= p;
+    D[0][2] -= p; D[2][0] -= p;
+
+    D[3][3] -= -p;
+    D[4][4] -= -p;
+    D[5][5] -= -p;
+
+    return tens4dmm(D);
 }
 
 //-----------------------------------------------------------------------------
@@ -250,3 +299,89 @@ tens4dmm FEBiphasic::Tangent_Permeability_Strain(FEMaterialPoint& mp)
 {
     return m_pPerm->Tangent_Permeability_Strain(mp);
 }
+
+//! return the material permeability property
+mat3ds FEBiphasic::MaterialPermeability(FEMaterialPoint& mp, const mat3ds E)
+{
+    // Evaluate right Cauchy-Green tensor from E
+    mat3ds C = mat3dd(1) + E*2;
+
+    // Evaluate right stretch tensor U from C
+    vec3d v[3];
+    double lam[3];
+    C.eigen2(lam, v);
+    lam[0] = sqrt(lam[0]); lam[1] = sqrt(lam[1]); lam[2] = sqrt(lam[2]);
+    mat3ds U = dyad(v[0])*lam[0] + dyad(v[1])*lam[1] + dyad(v[2])*lam[2];
+    double J = lam[0]*lam[1]*lam[2];
+
+    // temporarily replace F in material point with U
+    FEElasticMaterialPoint& pt = *mp.ExtractData<FEElasticMaterialPoint>();
+    mat3ds Ui = dyad(v[0])/lam[0] + dyad(v[1])/lam[1] + dyad(v[2])/lam[2];
+    mat3d Fsafe = pt.m_F;
+    double Jsafe = pt.m_J;
+    pt.m_F = U;
+    pt.m_J = J;
+
+    // Evaluate hydraulic permeability
+    mat3ds k = Permeability(mp);
+
+    // Restore original F
+    pt.m_F = Fsafe;
+    pt.m_J = Jsafe;
+
+    // Convert spatial permeability to material permeability
+    mat3ds K = (Ui*k*Ui).sym()*J;
+
+    return K;
+}
+
+//-----------------------------------------------------------------------------
+//! calculate spatial tangent stiffness at material point, using secant method
+tens4dmm FEBiphasic::SecantTangent_Permeability_Strain(FEMaterialPoint& mp)
+{
+    // extract the deformation gradient
+    FEElasticMaterialPoint& pt = *mp.ExtractData<FEElasticMaterialPoint>();
+    mat3d F = pt.m_F;
+    double J = pt.m_J;
+    mat3ds E = pt.Strain();
+    mat3dd I(1);
+
+    // calculate the material permeability at the current deformation gradient
+    mat3ds K = MaterialPermeability(mp,E);
+
+    // create deformation gradient increment
+    double eps = 1e-9;
+    vec3d e[3];
+    e[0] = vec3d(1,0,0); e[1] = vec3d(0,1,0); e[2] = vec3d(0,0,1);
+    tens4dmm Kmm;
+    for (int k=0; k<3; ++k) {
+        // evaluate incremental material permeability
+        mat3ds dE = dyads(e[k], e[k])*(eps/2);
+        mat3ds dK = (MaterialPermeability(mp,E+dE) - K)/eps;
+
+        // evaluate the secant modulus
+        Kmm(0,0,k,k) = dK.xx();
+        Kmm(1,1,k,k) = dK.yy();
+        Kmm(2,2,k,k) = dK.zz();
+        Kmm(0,1,k,k) = Kmm(1,0,k,k) = dK.xy();
+        Kmm(1,2,k,k) = Kmm(2,1,k,k) = dK.yz();
+        Kmm(2,0,k,k) = Kmm(0,2,k,k) = dK.xz();
+        for (int l=0; l<3; ++l) {
+            if (l != k) {
+                // evaluate incremental material permeability
+                mat3ds dE = dyads(e[k], e[l])*(eps/2);
+                mat3ds dK = (MaterialPermeability(mp,E+dE) - K)/eps;
+
+                // evaluate the secant modulus
+                Kmm(0,0,k,l) = Kmm(0,0,l,k) = dK.xx();
+                Kmm(1,1,k,l) = Kmm(1,1,l,k) = dK.yy();
+                Kmm(2,2,k,l) = Kmm(2,2,l,k) = dK.zz();
+                Kmm(0,1,k,l) = Kmm(0,1,l,k) = Kmm(1,0,k,l) = Kmm(1,0,l,k) = dK.xy();
+                Kmm(1,2,k,l) = Kmm(1,2,l,k) = Kmm(2,1,k,l) = Kmm(2,1,l,k) = dK.yz();
+                Kmm(2,0,k,l) = Kmm(2,0,l,k) = Kmm(0,2,k,l) = Kmm(0,2,l,k) = dK.xz();
+            }
+        }
+    }
+
+    return Kmm.pp(F)/J;
+}
diff --git a/FEBioMix/FEBiphasic.h b/FEBioMix/FEBiphasic.h
@@ -119,9 +119,15 @@ class FEBIOMIX_API FEBiphasic : public FEMaterial, public FEBiphasicInterface
 	//! calculate stress at material point
 	mat3ds Stress(FEMaterialPoint& pt);
 
+    //! calculate secant stress at material point
+    mat3ds SecantStress(FEMaterialPoint& pt);
+
 	//! calculate tangent stiffness at material point
-	tens4ds Tangent(FEMaterialPoint& pt);
+	tens4dmm Tangent(FEMaterialPoint& pt);
 
+    //! calculate secant tangent stiffness at material point
+    tens4dmm SecantTangent(FEMaterialPoint& pt);
+
 	//! return the permeability tensor as a matrix
 	void Permeability(double k[3][3], FEMaterialPoint& pt);
 
@@ -131,9 +137,15 @@ class FEBIOMIX_API FEBiphasic : public FEMaterial, public FEBiphasicInterface
     //! return tangent of permeability with strain
     tens4dmm Tangent_Permeability_Strain(FEMaterialPoint& mp);
 
+    //! return secant tangent of permeability with strain
+    tens4dmm SecantTangent_Permeability_Strain(FEMaterialPoint& mp);
+
 	//! return the permeability property
 	FEHydraulicPermeability* GetPermeability() { return m_pPerm; }
 
+    //! return the material permeability property
+    mat3ds MaterialPermeability(FEMaterialPoint& mp, const mat3ds E);
+
 	//! calculate actual fluid pressure
 	double Pressure(FEMaterialPoint& pt);
 

diff --git a/FEBioMix/FEBiphasicShellDomain.cpp b/FEBioMix/FEBiphasicShellDomain.cpp
@@ -34,12 +34,21 @@ SOFTWARE.*/
 #include <FECore/FESolidDomain.h>
 #include <FECore/FELinearSystem.h>
 
+BEGIN_FECORE_CLASS(FEBiphasicShellDomain, FESSIShellDomain)
+    ADD_PARAMETER(m_secant_stress, "secant_stress");
+    ADD_PARAMETER(m_secant_tangent, "secant_tangent");
+    ADD_PARAMETER(m_secant_perm_tangent, "secant_permeability_tangent");
+END_FECORE_CLASS();
+
 //-----------------------------------------------------------------------------
 FEBiphasicShellDomain::FEBiphasicShellDomain(FEModel* pfem) : FESSIShellDomain(pfem), FEBiphasicDomain(pfem), m_dof(pfem)
 {
     m_dofSX = pfem->GetDOFIndex("sx");
     m_dofSY = pfem->GetDOFIndex("sy");
     m_dofSZ = pfem->GetDOFIndex("sz");
+    m_secant_stress = false;
+    m_secant_tangent = false;
+    m_secant_perm_tangent = false;
 }
 
 //-----------------------------------------------------------------------------
@@ -548,15 +557,15 @@ bool FEBiphasicShellDomain::ElementBiphasicStiffness(FEShellElement& el, matrix&
         mat3ds s = ept.m_s;
 
         // get elasticity tensor
-        tens4ds c = m_pMat->Tangent(mp);
-        
+        tens4dmm c = (m_secant_tangent) ? m_pMat->SecantTangent(mp) : m_pMat->Tangent(mp);
+
         // get the fluid flux and pressure gradient
         vec3d gradp = pt.m_gradp + (pt.m_gradp - pt.m_gradpp)*(tau/dt);
 
         // evaluate the permeability and its derivatives
         mat3ds K = m_pMat->Permeability(mp);
-        tens4dmm dKdE = m_pMat->Tangent_Permeability_Strain(mp);
-        
+        tens4dmm dKdE = (m_secant_perm_tangent) ? m_pMat->SecantTangent_Permeability_Strain(mp) : m_pMat->Tangent_Permeability_Strain(mp);
+
         // evaluate the solvent supply and its derivatives
         double phiwhat = 0;
         mat3ds Phie; Phie.zero();
@@ -731,15 +740,15 @@ bool FEBiphasicShellDomain::ElementBiphasicStiffnessSS(FEShellElement& el, matri
         mat3ds s = ept.m_s;
 
         // get elasticity tensor
-        tens4ds c = m_pMat->Tangent(mp);
-        
+        tens4dmm c = (m_secant_tangent) ? m_pMat->SecantTangent(mp) : m_pMat->Tangent(mp);
+
         // get the fluid flux and pressure gradient
         vec3d gradp = pt.m_gradp;
 
         // evaluate the permeability and its derivatives
         mat3ds K = m_pMat->Permeability(mp);
-        tens4dmm dKdE = m_pMat->Tangent_Permeability_Strain(mp);
-        
+        tens4dmm dKdE = (m_secant_perm_tangent) ? m_pMat->SecantTangent_Permeability_Strain(mp) : m_pMat->Tangent_Permeability_Strain(mp);
+
         // evaluate the solvent supply and its derivatives
         double phiwhat = 0;
         mat3ds Phie; Phie.zero();
@@ -942,10 +951,10 @@ void FEBiphasicShellDomain::UpdateElementStress(int iel)
         m_pMat->UpdateSpecializedMaterialPoints(mp, GetFEModel()->GetTime());
 
         // calculate the solid stress at this material point
-        ppt.m_ss = m_pMat->GetElasticMaterial()->Stress(mp);
-        
+        ppt.m_ss = (m_secant_stress) ? m_pMat->GetElasticMaterial()->SecantStress(mp) : m_pMat->GetElasticMaterial()->Stress(mp);
+
         // calculate the stress at this material point
-        pt.m_s = m_pMat->Stress(mp);
+        pt.m_s = (m_secant_stress) ? m_pMat->SecantStress(mp) : m_pMat->Stress(mp);
     }
 }
 

diff --git a/FEBioMix/FEBiphasicShellDomain.h b/FEBioMix/FEBiphasicShellDomain.h
@@ -114,9 +114,16 @@ class FEBIOMIX_API FEBiphasicShellDomain : public FESSIShellDomain, public FEBip
     // That is why I've taken this calculation out of the FEBiphasic class and placed it here.
     vec3d FluidFlux(FEMaterialPoint& mp) override;
 
+protected:
+    bool    m_secant_stress;    //!< use secant approximation to stress
+    bool    m_secant_tangent;   //!< flag for using secant tangent
+    bool    m_secant_perm_tangent;   //!< flag for using secant tangent on permeability
+
 protected:
     int					m_dofSX;
     int					m_dofSY;
     int					m_dofSZ;
 	FEDofList			m_dof;
+
+    DECLARE_FECORE_CLASS();
 };
diff --git a/FEBioMix/FEBiphasicSolidDomain.cpp b/FEBioMix/FEBiphasicSolidDomain.cpp
@@ -36,10 +36,19 @@ SOFTWARE.*/
 #include <FECore/FELinearSystem.h>
 #include "FEBioMix.h"
 
+BEGIN_FECORE_CLASS(FEBiphasicSolidDomain, FESolidDomain)
+    ADD_PARAMETER(m_secant_stress, "secant_stress");
+    ADD_PARAMETER(m_secant_tangent, "secant_tangent");
+    ADD_PARAMETER(m_secant_perm_tangent, "secant_permeability_tangent");
+END_FECORE_CLASS();
+
 //-----------------------------------------------------------------------------
 FEBiphasicSolidDomain::FEBiphasicSolidDomain(FEModel* pfem) : FESolidDomain(pfem), FEBiphasicDomain(pfem), m_dofU(pfem), m_dofSU(pfem), m_dofR(pfem), m_dof(pfem)
 {
 	EXPORT_DATA(PLT_FLOAT, FMT_NODE, &m_nodePressure, "NPR fluid pressure");
+    m_secant_stress = false;
+    m_secant_tangent = false;
+    m_secant_perm_tangent = false;
 
 	if (pfem)
 	{
@@ -641,15 +650,15 @@ bool FEBiphasicSolidDomain::ElementBiphasicStiffness(FESolidElement& el, matrix&
         mat3ds s = ept.m_s;
 
         // get elasticity tensor
-        tens4ds c = m_pMat->Tangent(mp);
-        
+        tens4dmm c = (m_secant_tangent) ? m_pMat->SecantTangent(mp) : m_pMat->Tangent(mp);
+
         // get the fluid flux and pressure gradient
         vec3d gradp = pt.m_gradp + (pt.m_gradp - pt.m_gradpp)*(tau/dt);
 
         // evaluate the permeability and its derivatives
         mat3ds K = m_pMat->Permeability(mp);
-        tens4dmm dKdE = m_pMat->Tangent_Permeability_Strain(mp);
-        
+        tens4dmm dKdE = (m_secant_perm_tangent) ? m_pMat->SecantTangent_Permeability_Strain(mp) : m_pMat->Tangent_Permeability_Strain(mp);
+
         // evaluate the solvent supply and its derivatives
         double phiwhat = 0;
         mat3ds Phie; Phie.zero();
@@ -796,15 +805,15 @@ bool FEBiphasicSolidDomain::ElementBiphasicStiffnessSS(FESolidElement& el, matri
         mat3ds s = ept.m_s;
 
         // get elasticity tensor
-        tens4ds c = m_pMat->Tangent(mp);
-        
+        tens4dmm c = (m_secant_tangent) ? m_pMat->SecantTangent(mp) : m_pMat->Tangent(mp);
+
         // get the fluid flux and pressure gradient
         vec3d gradp = pt.m_gradp;
 
         // evaluate the permeability and its derivatives
         mat3ds K = m_pMat->Permeability(mp);
-        tens4dmm dKdE = m_pMat->Tangent_Permeability_Strain(mp);
-        
+        tens4dmm dKdE = (m_secant_perm_tangent) ? m_pMat->SecantTangent_Permeability_Strain(mp) : m_pMat->Tangent_Permeability_Strain(mp);
+
         // evaluate the solvent supply and its derivatives
         double phiwhat = 0;
         mat3ds Phie; Phie.zero();
@@ -981,10 +990,10 @@ void FEBiphasicSolidDomain::UpdateElementStress(int iel)
         m_pMat->UpdateSpecializedMaterialPoints(mp, GetFEModel()->GetTime());
 
         // calculate the solid stress at this material point
-        ppt.m_ss = m_pMat->GetElasticMaterial()->Stress(mp);
-        
+        ppt.m_ss = (m_secant_stress) ? m_pMat->GetElasticMaterial()->SecantStress(mp) : m_pMat->GetElasticMaterial()->Stress(mp);
+
 		// calculate the stress at this material point
-		pt.m_s = m_pMat->Stress(mp);
+		pt.m_s = (m_secant_stress) ? m_pMat->SecantStress(mp) : m_pMat->Stress(mp);
 	}
 }
 

diff --git a/FEBioMix/FEBiphasicSolidDomain.h b/FEBioMix/FEBiphasicSolidDomain.h
@@ -135,11 +135,18 @@ class FEBIOMIX_API FEBiphasicSolidDomain : public FESolidDomain, public FEBiphas
 	// This function updates the m_nodePressure variable
 	void UpdateNodalPressures();
 
+protected:
+    bool    m_secant_stress;    //!< use secant approximation to stress
+    bool    m_secant_tangent;   //!< flag for using secant tangent
+    bool    m_secant_perm_tangent;   //!< flag for using secant tangent on permeability
+
 protected:
 	int			m_varU, m_varP;	// displacement, pressure field indices
 
 	FEDofList	m_dofU;		// displacement dofs
 	FEDofList	m_dofSU;	// shell displacement dofs
 	FEDofList	m_dofR;		// rigid rotation
 	FEDofList	m_dof;
+
+    DECLARE_FECORE_CLASS();
 };