Skip to content

Commit

Permalink
Modified FiberExpPowSBM and FiberPowLinearSBM to use FEParamDouble pr…
Browse files Browse the repository at this point in the history 
…operties. Introduced ReactionRateRuberti to model strain-dependent fiber degradation
  • Loading branch information
@gateshian
gateshian committed Apr 22, 2024
1 parent 26a430e commit 654e815
Show file tree
Hide file tree
Showing 7 changed files with 265 additions and 44 deletions.
2 changes: 2 additions & 0 deletions FEBioMix/FEBioMix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ SOFTWARE.*/
#include "FECarterHayes.h"
#include "FEReactionRateConst.h"
#include "FEReactionRateHuiskes.h"
#include "FEReactionRateRuberti.h"
#include "FEReactionRateNims.h"
#include "FEReactionRateExpSED.h"
#include "FEReactionRateSoluteAsSBM.h"
Expand Down Expand Up @@ -452,6 +453,7 @@ void FEBioMix::InitModule()
REGISTER_FECORE_CLASS(FESFDSBM , "spherical fiber distribution sbm");
REGISTER_FECORE_CLASS(FEReactionRateConst , "constant reaction rate" );
REGISTER_FECORE_CLASS(FEReactionRateHuiskes , "Huiskes reaction rate" );
REGISTER_FECORE_CLASS(FEReactionRateRuberti , "Ruberti reaction rate" );
REGISTER_FECORE_CLASS(FEReactionRateNims , "Nims reaction rate" );
REGISTER_FECORE_CLASS(FEReactionRateExpSED , "exp-sed reaction rate" );
REGISTER_FECORE_CLASS(FEReactionRateSoluteAsSBM , "solute-as-sbm reaction rate");
Expand Down
31 changes: 20 additions & 11 deletions FEBioMix/FEFiberExpPowSBM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,7 @@ mat3ds FEFiberExpPowSBM::Stress(FEMaterialPoint& mp)

// initialize material constants
double rhor = spt.m_sbmr[m_lsbm];
double ksi = FiberModulus(rhor);
double ksi = FiberModulus(mp, rhor);

// deformation gradient
mat3d &F = pt.m_F;
Expand Down Expand Up @@ -134,14 +134,17 @@ mat3ds FEFiberExpPowSBM::Stress(FEMaterialPoint& mp)
// only take fibers in tension into consideration
if (In_1 >= eps)
{
double alpha = m_alpha(mp);
double beta = m_beta(mp);

// get the global spatial fiber direction in current configuration
nt = F*n0;

// calculate the outer product of nt
mat3ds N = dyad(nt);

// calculate strain energy derivative
Wl = ksi*pow(In_1, m_beta-1.0)*exp(m_alpha*pow(In_1, m_beta));
Wl = ksi*pow(In_1, beta-1.0)*exp(alpha*pow(In_1, beta));

// calculate the fiber stress
s = N*(2.0*Wl/J);
Expand All @@ -162,7 +165,7 @@ tens4ds FEFiberExpPowSBM::Tangent(FEMaterialPoint& mp)

// initialize material constants
double rhor = spt.m_sbmr[m_lsbm];
double ksi = FiberModulus(rhor);
double ksi = FiberModulus(mp, rhor);

// deformation gradient
mat3d &F = pt.m_F;
Expand Down Expand Up @@ -190,6 +193,9 @@ tens4ds FEFiberExpPowSBM::Tangent(FEMaterialPoint& mp)
// only take fibers in tension into consideration
if (In_1 >= eps)
{
double alpha = m_alpha(mp);
double beta = m_beta(mp);

// get the global spatial fiber direction in current configuration
nt = F*n0;

Expand All @@ -198,8 +204,8 @@ tens4ds FEFiberExpPowSBM::Tangent(FEMaterialPoint& mp)
tens4ds NxN = dyad1s(N);

// calculate strain energy 2nd derivative
double tmp = m_alpha*pow(In_1, m_beta);
Wll = ksi*pow(In_1, m_beta-2.0)*((tmp+1)*m_beta-1.0)*exp(tmp);
double tmp = alpha*pow(In_1, beta);
Wll = ksi*pow(In_1, beta-2.0)*((tmp+1)*beta-1.0)*exp(tmp);

// calculate the fiber tangent
c = NxN*(4.0*Wll/J);
Expand All @@ -222,7 +228,7 @@ double FEFiberExpPowSBM::StrainEnergyDensity(FEMaterialPoint& mp)

// initialize material constants
double rhor = spt.m_sbmr[m_lsbm];
double ksi = FiberModulus(rhor);
double ksi = FiberModulus(mp, rhor);

// loop over all integration points
double In_1;
Expand All @@ -244,12 +250,15 @@ double FEFiberExpPowSBM::StrainEnergyDensity(FEMaterialPoint& mp)
// only take fibers in tension into consideration
if (In_1 >= eps)
{
double alpha = m_alpha(mp);
double beta = m_beta(mp);

// calculate strain energy derivative
if (m_alpha > 0) {
sed = ksi/(m_alpha*m_beta)*(exp(m_alpha*pow(In_1, m_beta))-1);
if (alpha > 0) {
sed = ksi/(alpha*beta)*(exp(alpha*pow(In_1, beta))-1);
}
else
sed = ksi/m_beta*pow(In_1, m_beta);
sed = ksi/beta*pow(In_1, beta);
}

return sed;
Expand Down Expand Up @@ -284,7 +293,7 @@ double FEFiberExpPowSBM::Tangent_SE_Density(FEMaterialPoint& mp)
{
FESolutesMaterialPoint& spt = *mp.ExtractData<FESolutesMaterialPoint>();
double rhor = spt.m_sbmr[m_lsbm];
return StrainEnergy(mp)*m_g/rhor;
return StrainEnergy(mp)*m_g(mp)/rhor;
}

//-----------------------------------------------------------------------------
Expand All @@ -293,6 +302,6 @@ mat3ds FEFiberExpPowSBM::Tangent_Stress_Density(FEMaterialPoint& mp)
{
FESolutesMaterialPoint& spt = *mp.ExtractData<FESolutesMaterialPoint>();
double rhor = spt.m_sbmr[m_lsbm];
return Stress(mp)*m_g/rhor;
return Stress(mp)*m_g(mp)/rhor;
}

16 changes: 8 additions & 8 deletions FEBioMix/FEFiberExpPowSBM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ SOFTWARE.*/
#pragma once
#include <FEBioMech/FEElasticMaterial.h>
#include <FEBioMech/FERemodelingElasticMaterial.h>
#include <FECore/FEModelParam.h>
#include "febiomix_api.h"

//-----------------------------------------------------------------------------
Expand Down Expand Up @@ -57,7 +58,7 @@ class FEBIOMIX_API FEFiberExpPowSBM : public FEElasticMaterial, public FERemodel
double Density(FEMaterialPoint& pt) override;

//! return fiber modulus
double FiberModulus(double rhor) { return m_ksi0*pow(rhor/m_rho0, m_g);}
double FiberModulus(FEMaterialPoint& pt, double rhor) { return m_ksi0(pt)*pow(rhor/m_rho0(pt), m_g(pt));}

//! Create material point data
FEMaterialPointData* CreateMaterialPointData() override;
Expand All @@ -77,13 +78,12 @@ class FEBIOMIX_API FEFiberExpPowSBM : public FEElasticMaterial, public FERemodel
mat3ds Tangent_Stress_Density(FEMaterialPoint& pt) override;

public:
double m_alpha; // coefficient of (In-1) in exponential
double m_beta; // power of (In-1) in exponential
double m_ksi0; // fiber modulus ksi = ksi0*(rhor/rho0)^gamma
double m_rho0; // rho0
double m_g; // gamma
// int m_sbm; //!< global id of solid-bound molecule
int m_lsbm; //!< local id of solid-bound molecule
FEParamDouble m_alpha; // coefficient of (In-1) in exponential
FEParamDouble m_beta; // power of (In-1) in exponential
FEParamDouble m_ksi0; // fiber modulus ksi = ksi0*(rhor/rho0)^gamma
FEParamDouble m_rho0; // rho0
FEParamDouble m_g; // gamma
int m_lsbm; //!< local id of solid-bound molecule

public:
FEVec3dValuator* m_fiber; //!< fiber orientation
Expand Down
43 changes: 26 additions & 17 deletions FEBioMix/FEFiberPowLinearSBM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,12 +105,15 @@ mat3ds FEFiberPowLinearSBM::Stress(FEMaterialPoint& mp)
FEElasticMaterialPoint& pt = *mp.ExtractData<FEElasticMaterialPoint>();
FESolutesMaterialPoint& spt = *mp.ExtractData<FESolutesMaterialPoint>();

double lam0 = m_lam0(mp);
double beta = m_beta(mp);

// initialize material constants
double rhor = spt.m_sbmr[m_lsbm];
double E = FiberModulus(rhor);
double I0 = m_lam0*m_lam0;
double ksi = E/4/(m_beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-m_beta);
double b = ksi*pow(I0-1, m_beta-1) + E/2/sqrt(I0);
double E = FiberModulus(mp, rhor);
double I0 = lam0*lam0;
double ksi = E/4/(beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-beta);
double b = ksi*pow(I0-1, beta-1) + E/2/sqrt(I0);

// deformation gradient
mat3d &F = pt.m_F;
Expand Down Expand Up @@ -146,7 +149,7 @@ mat3ds FEFiberPowLinearSBM::Stress(FEMaterialPoint& mp)

// calculate the fiber stress magnitude
sn = (In < I0) ?
2*In*ksi*pow(In-1, m_beta-1) :
2*In*ksi*pow(In-1, beta-1) :
2*b*In - E*sqrt(In);

// calculate the fiber stress
Expand All @@ -166,11 +169,14 @@ tens4ds FEFiberPowLinearSBM::Tangent(FEMaterialPoint& mp)
FEElasticMaterialPoint& pt = *mp.ExtractData<FEElasticMaterialPoint>();
FESolutesMaterialPoint& spt = *mp.ExtractData<FESolutesMaterialPoint>();

double lam0 = m_lam0(mp);
double beta = m_beta(mp);

// initialize material constants
double rhor = spt.m_sbmr[m_lsbm];
double E = FiberModulus(rhor);
double I0 = m_lam0*m_lam0;
double ksi = E/4/(m_beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-m_beta);
double E = FiberModulus(mp, rhor);
double I0 = lam0*lam0;
double ksi = E/4/(beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-beta);

// deformation gradient
mat3d &F = pt.m_F;
Expand Down Expand Up @@ -207,7 +213,7 @@ tens4ds FEFiberPowLinearSBM::Tangent(FEMaterialPoint& mp)

// calculate modulus
cn = (In < I0) ?
4*In*In*ksi*(m_beta-1)*pow(In-1, m_beta-2) :
4*In*In*ksi*(beta-1)*pow(In-1, beta-2) :
E*sqrt(In);

// calculate the fiber tangent
Expand All @@ -229,12 +235,15 @@ double FEFiberPowLinearSBM::StrainEnergyDensity(FEMaterialPoint& mp)
FEElasticMaterialPoint& pt = *mp.ExtractData<FEElasticMaterialPoint>();
FESolutesMaterialPoint& spt = *mp.ExtractData<FESolutesMaterialPoint>();

double lam0 = m_lam0(mp);
double beta = m_beta(mp);

// initialize material constants
double rhor = spt.m_sbmr[m_lsbm];
double E = FiberModulus(rhor);
double I0 = m_lam0*m_lam0;
double ksi = E/4/(m_beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-m_beta);
double b = ksi*pow(I0-1, m_beta-1) + E/2/sqrt(I0);
double E = FiberModulus(mp, rhor);
double I0 = lam0*lam0;
double ksi = E/4/(beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-beta);
double b = ksi*pow(I0-1, beta-1) + E/2/sqrt(I0);

// loop over all integration points
double In;
Expand All @@ -258,8 +267,8 @@ double FEFiberPowLinearSBM::StrainEnergyDensity(FEMaterialPoint& mp)
{
// calculate strain energy density
sed = (In < I0) ?
ksi/m_beta*pow(In-1, m_beta) :
b*(In-I0) - E*(sqrt(In)-sqrt(I0)) + ksi/m_beta*pow(I0-1, m_beta);
ksi/beta*pow(In-1, beta) :
b*(In-I0) - E*(sqrt(In)-sqrt(I0)) + ksi/beta*pow(I0-1, beta);
}

return sed;
Expand Down Expand Up @@ -294,7 +303,7 @@ double FEFiberPowLinearSBM::Tangent_SE_Density(FEMaterialPoint& mp)
{
FESolutesMaterialPoint& spt = *mp.ExtractData<FESolutesMaterialPoint>();
double rhor = spt.m_sbmr[m_lsbm];
return StrainEnergy(mp)*m_g/rhor;
return StrainEnergy(mp)*m_g(mp)/rhor;
}

//-----------------------------------------------------------------------------
Expand All @@ -303,6 +312,6 @@ mat3ds FEFiberPowLinearSBM::Tangent_Stress_Density(FEMaterialPoint& mp)
{
FESolutesMaterialPoint& spt = *mp.ExtractData<FESolutesMaterialPoint>();
double rhor = spt.m_sbmr[m_lsbm];
return Stress(mp)*m_g/rhor;
return Stress(mp)*m_g(mp)/rhor;
}

16 changes: 8 additions & 8 deletions FEBioMix/FEFiberPowLinearSBM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ SOFTWARE.*/
#pragma once
#include <FEBioMech/FEElasticMaterial.h>
#include <FEBioMech/FERemodelingElasticMaterial.h>
#include <FECore/FEModelParam.h>
#include "febiomix_api.h"

//-----------------------------------------------------------------------------
Expand Down Expand Up @@ -57,7 +58,7 @@ class FEBIOMIX_API FEFiberPowLinearSBM : public FEElasticMaterial, public FERemo
double Density(FEMaterialPoint& pt) override;

//! return fiber modulus
double FiberModulus(double rhor) { return m_E0*pow(rhor/m_rho0, m_g);}
double FiberModulus(FEMaterialPoint& pt, double rhor) { return m_E0(pt)*pow(rhor/m_rho0(pt), m_g(pt));}

//! Create material point data
FEMaterialPointData* CreateMaterialPointData() override;
Expand All @@ -77,13 +78,12 @@ class FEBIOMIX_API FEFiberPowLinearSBM : public FEElasticMaterial, public FERemo
mat3ds Tangent_Stress_Density(FEMaterialPoint& pt) override;

public:
double m_E0; // fiber modulus E = E0*(rhor/rho0)^gamma
double m_lam0; // stretch ratio at end of toe region
double m_beta; // power law exponent in toe region
double m_rho0; // rho0
double m_g; // gamma
// int m_sbm; //!< global id of solid-bound molecule
int m_lsbm; //!< local id of solid-bound molecule
FEParamDouble m_E0; // fiber modulus E = E0*(rhor/rho0)^gamma
FEParamDouble m_lam0; // stretch ratio at end of toe region
FEParamDouble m_beta; // power law exponent in toe region
FEParamDouble m_rho0; // rho0
FEParamDouble m_g; // gamma
int m_lsbm; //!< local id of solid-bound molecule

public:
FEVec3dValuator* m_fiber; //!< fiber orientation
Expand Down
Loading

0 comments on commit 654e815

Please sign in to comment.