From ae2762e51b8a3ff8b129ce1a2bdc568ea23bb6fd Mon Sep 17 00:00:00 2001
From: FLE_Knight <ujwalkishore.jinaga@mail.polimi.it>
Date: Wed, 31 May 2023 20:08:54 +0200
Subject: [PATCH] Modified formulation for mullins effect in mullinsEffect.cpp;
added psiMax variable to ipMullins
---
.../internalPoints/ipMullinsEffect.cpp | 5 +-
.../internalPoints/ipMullinsEffect.h | 5 ++
.../materialLaw/mlawNonLinearTVE.cpp | 24 ++++++--
.../materialLaw/mlawNonLinearTVE.h | 42 ++++++--------
.../materialLaw/mlawNonLinearTVP.cpp | 39 ++++++++-----
NonLinearSolver/materialLaw/mullinsEffect.cpp | 58 ++++++++++++++++---
6 files changed, 122 insertions(+), 51 deletions(-)
diff --git a/NonLinearSolver/internalPoints/ipMullinsEffect.cpp b/NonLinearSolver/internalPoints/ipMullinsEffect.cpp
index f53d5384b..b4bba50b4 100644
--- a/NonLinearSolver/internalPoints/ipMullinsEffect.cpp
+++ b/NonLinearSolver/internalPoints/ipMullinsEffect.cpp
@@ -10,7 +10,7 @@
#include "ipMullinsEffect.h"
-IPMullinsEffect::IPMullinsEffect(): eta(1.0), DetaDpsi(0.), DDetaDpsipsi(0.), DetaDT(0.), DDetaDTT(0.)
+IPMullinsEffect::IPMullinsEffect(): eta(1.0), DetaDpsi(0.), DDetaDpsipsi(0.), DetaDT(0.), DDetaDTT(0.), psiMax(0.)
{
}
@@ -22,6 +22,7 @@ IPMullinsEffect::IPMullinsEffect(const IPMullinsEffect &source)
DDetaDpsipsi = source.DDetaDpsipsi;
DetaDT = source.DetaDT;
DDetaDTT = source.DDetaDTT;
+ psiMax = source.psiMax;
}
IPMullinsEffect &IPMullinsEffect::operator=(const IPMullinsEffect &source)
@@ -31,6 +32,7 @@ IPMullinsEffect &IPMullinsEffect::operator=(const IPMullinsEffect &source)
DDetaDpsipsi = source.DDetaDpsipsi;
DetaDT = source.DetaDT;
DDetaDTT = source.DDetaDTT;
+ psiMax = source.psiMax;
return *this;
}
@@ -41,5 +43,6 @@ void IPMullinsEffect::restart()
restartManager::restart(DDetaDpsipsi);
restartManager::restart(DetaDT);
restartManager::restart(DDetaDTT);
+ restartManager::restart(psiMax);
return;
}
diff --git a/NonLinearSolver/internalPoints/ipMullinsEffect.h b/NonLinearSolver/internalPoints/ipMullinsEffect.h
index 153ce88e3..1aa954613 100644
--- a/NonLinearSolver/internalPoints/ipMullinsEffect.h
+++ b/NonLinearSolver/internalPoints/ipMullinsEffect.h
@@ -21,6 +21,7 @@ class IPMullinsEffect{
double DDetaDpsipsi;
double DetaDT;
double DDetaDTT;
+ double psiMax;
public:
IPMullinsEffect();
@@ -36,6 +37,7 @@ class IPMullinsEffect{
virtual double getDDetaDpsipsi() const {return DDetaDpsipsi;};
virtual double getDetaDT() const{ return DetaDT;};
virtual double getDDetaDTT() const{ return DDetaDTT;};
+ virtual double getpsiMax() const{ return psiMax;};
virtual void set(const double &_eta, const double &_DetaDpsi, const double& _DDetaDpsipsi, const double& _DetaDT, const double& _DDetaDTT)
{
eta=_eta;
@@ -44,6 +46,9 @@ class IPMullinsEffect{
DetaDT = _DetaDT;
DDetaDTT = _DDetaDTT;
}
+ virtual void setPsiMax(const double &_psiMax){
+ psiMax = _psiMax;
+ }
};
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
index 706f2ff61..248de9c7c 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
@@ -591,6 +591,17 @@ double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const dou
void mlawNonLinearTVE::calculateMullinsEffectScaler(const IPNonLinearTVE* q0, IPNonLinearTVE *q1, const double T) const{
+ // Set _psiMax
+ /*
+ if (q1->_elasticEnergy > q0->_psiMax){
+ q1->_psiMax = q1->_elasticEnergy;
+ }
+ else{
+ q1->_psiMax = q0->_psiMax;
+ }
+ Msg::Error(" Inside calculateMullinsEffectScaler, before updating eta, q1->_elasticEnergy = %e!!",q1->_elasticEnergy);
+ Msg::Error(" Inside calculateMullinsEffectScaler, before updating eta, q1->_psiMax = %e!!",q1->_psiMax);*/
+
// Initialise
double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.;
@@ -1040,11 +1051,13 @@ void mlawNonLinearTVE::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
double& Ke, double& Ge, double& Kde, double& Gde,
double& DKe, double& DGe, double& DKde, double& DGde,
double& KTsum, double& GTsum, double& DKDTsum, double& DGDTsum,
- const double T0, const double T1, const bool calculateStress) const{
+ const double T0, const double T1, const bool calculateStress, const bool mullins) const{
- // Mullins Effect Scaler -> should be before any calls to get functions
- if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){
- mlawNonLinearTVE::calculateMullinsEffectScaler(q0, q1, T1);
+ // Mullins Effect Scaler -> should be before any calls to "get" functions
+ if(mullins){
+ if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){
+ mlawNonLinearTVE::calculateMullinsEffectScaler(q0, q1, T1);
+ }
}
// Set Temperature to calculate temperature updated properties
@@ -1515,12 +1528,13 @@ void mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(
// elastic energy
q1->_elasticEnergy = freeEnergyMechanical(*q1,T0,T);
+ /*
if (q1->_elasticEnergy > q0->_psiMax){
q1->_psiMax = q1->_elasticEnergy;
}
else{
q1->_psiMax = q0->_psiMax;
- }
+ }*/
// Msg::Error(" Inside TVE predCorr, after updating q1->_elasticEnergy = %e!!",q1->_elasticEnergy);
// thermal energy
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.h b/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
index 26073b30d..ef2ff1b12 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
@@ -88,7 +88,7 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
double& Ke, double& Ge, double& Kde, double& Gde,
double& DKe, double& DGe, double& DKde, double& DGde,
double& KTsum, double& GTsum, double& DKDTsum, double& DGDTsum,
- const double T0, const double T1, const bool calculateStress = true) const;
+ const double T0, const double T1, const bool calculateStress = true, const bool mullins = true) const;
void isotropicHookTensor(const double K, const double G, STensor43& L) const;
@@ -228,24 +228,23 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
}
virtual void getK(const IPNonLinearTVE *q1, double& KT, const double Tc, double *dK = NULL, double *ddK = NULL) const{
- double K = _K*q1->getConstRefToElasticBulkPropertyScaleFactor();
- // Msg::Error(" Inside getK, after getting and updating eta; K = %e!!",K);
- KT = K*_temFunc_K->getVal(Tc);
+ // double K = _K*q1->getConstRefToElasticBulkPropertyScaleFactor();
+ KT = _K*_temFunc_K->getVal(Tc);
if(dK!=NULL){
- *dK = K*_temFunc_K->getDiff(Tc); // Here, *dK can be imagined as *(dK+0) or dK[0] if dK were a vector. See dKi below.
+ *dK = _K*_temFunc_K->getDiff(Tc); // Here, *dK can be imagined as *(dK+0) or dK[0] if dK were a vector. See dKi below.
}
if(ddK!=NULL){
- *ddK = K*_temFunc_K->getDoubleDiff(Tc);
+ *ddK = _K*_temFunc_K->getDoubleDiff(Tc);
}
}
virtual void getG(const IPNonLinearTVE *q1, double& GT, const double Tc, double *dG = NULL, double *ddG = NULL) const{
- double G = _G*q1->getConstRefToElasticShearPropertyScaleFactor();
- GT = G*_temFunc_G->getVal(Tc);
+ // double G = _G*q1->getConstRefToElasticShearPropertyScaleFactor();
+ GT = _G*_temFunc_G->getVal(Tc);
if(dG!=NULL){
- *dG = G*_temFunc_G->getDiff(Tc);
+ *dG = _G*_temFunc_G->getDiff(Tc);
}
if(ddG!=NULL){
- *ddG = G*_temFunc_G->getDoubleDiff(Tc);
+ *ddG = _G*_temFunc_G->getDoubleDiff(Tc);
}
}
virtual void getAlpha(double& alphaT, const double Tc, double *dAlpha = NULL, double *ddAlpha = NULL) const{
@@ -275,41 +274,38 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
// The 3rd argument is a pointer to a vector.
virtual void getKi(const IPNonLinearTVE *q1, std::vector<double>& Ki_T, const double Tc, std::vector<double>* dKi=NULL, std::vector<double>* ddKi=NULL) const{
- std::vector<double> Ki;
+ /*std::vector<double> Ki;
Ki.clear(); Ki.resize(_N,0.);
for (int i=0; i<_Ki.size();i++){
Ki[i] = _Ki[i]*q1->getConstRefToElasticBulkPropertyScaleFactor();
- }
- // Msg::Error(" Inside getKi, after getting and updating eta; Ki[0] = %e!!",Ki[0]);
+ }*/
Ki_T.resize(_N,0.);
for (int i=0; i<_Ki.size();i++){
- Ki_T[i] = Ki[i]*_temFunc_Ki[i]->getVal(Tc);
+ Ki_T[i] = _Ki[i]*_temFunc_Ki[i]->getVal(Tc);
if(dKi!=NULL){
- dKi->at(i) = Ki[i]*_temFunc_Ki[i]->getDiff(Tc);
+ dKi->at(i) = _Ki[i]*_temFunc_Ki[i]->getDiff(Tc);
}
if(ddKi!=NULL){
- ddKi->at(i) = Ki[i]*_temFunc_Ki[i]->getDoubleDiff(Tc);
+ ddKi->at(i) = _Ki[i]*_temFunc_Ki[i]->getDoubleDiff(Tc);
}
}
}
// The 3rd argument is a pointer to a vector.
virtual void getGi(const IPNonLinearTVE *q1, std::vector<double>& Gi_T, const double Tc, std::vector<double>* dGi=NULL, std::vector<double>* ddGi=NULL) const{
-
- std::vector<double> Gi;
+ /*std::vector<double> Gi;
Gi.clear(); Gi.resize(_N,0.);
for (int i=0; i<_Gi.size();i++){
Gi[i] = _Gi[i]*q1->getConstRefToElasticShearPropertyScaleFactor();
- }
-
+ }*/
Gi_T.resize(_N,0.);
for (int i=0; i<_Gi.size();i++){
- Gi_T[i] = Gi[i]*_temFunc_Gi[i]->getVal(Tc);
+ Gi_T[i] = _Gi[i]*_temFunc_Gi[i]->getVal(Tc);
if(dGi!=NULL){
- dGi->at(i) = Gi[i]*_temFunc_Gi[i]->getDiff(Tc); // *(dGi+i) is equivalent to dGi[i]
+ dGi->at(i) = _Gi[i]*_temFunc_Gi[i]->getDiff(Tc); // *(dGi+i) is equivalent to dGi[i]
}
if(ddGi!=NULL){
- ddGi->at(i) = Gi[i]*_temFunc_Gi[i]->getDoubleDiff(Tc);
+ ddGi->at(i) = _Gi[i]*_temFunc_Gi[i]->getDoubleDiff(Tc);
}
}
}
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index a67d67427..ce2a5bd7d 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -611,7 +611,7 @@ void mlawNonLinearTVP::getMechSourceTVP(const STensor3& F0, const STensor3& F,
std::vector<STensor3> ddIsoHardForcedFdT; ddIsoHardForcedFdT.resize(3,0.);
// getIsotropicHardeningForce(q0,q1,T0,T,DphiPDF,&ddIsoHardForcedTT,&ddIsoHardForcedFdT); // For analytical derivatives
- getIsotropicHardeningForce(q0,q1,T0,T,DphiPDF); // For numerical derivatives
+ // getIsotropicHardeningForce(q0,q1,T0,T,DphiPDF); // For numerical derivatives
const std::vector<double>& IsoHardForce = q1->_IsoHardForce;
const std::vector<double>& dIsoHardForcedT = q1->_dIsoHardForcedT;
const std::vector<STensor3>& dIsoHardForcedF = q1->_dIsoHardForcedF;
@@ -619,11 +619,14 @@ void mlawNonLinearTVP::getMechSourceTVP(const STensor3& F0, const STensor3& F,
double R(0.), dRdT(0.), ddRdTT(0.);
STensor3 dRdF, ddRdFdT; STensorOperation::zero(dRdF); STensorOperation::zero(ddRdFdT);
+ R += IsoHardForce[2];
+ dRdT += dIsoHardForcedT[2];
+ dRdF += dIsoHardForcedF[2];
for (int i=0; i<3; i++){
- R += IsoHardForce[i];
- dRdT += dIsoHardForcedT[i];
+ // R +=
+ // dRdT += dIsoHardForcedT[i];
ddRdTT += 0.; // ddIsoHardForcedTT[i];
- dRdF += dIsoHardForcedF[i];
+ // dRdF += dIsoHardForcedF[i];
ddRdFdT += 0.; // ddIsoHardForcedFdT[i];
}
@@ -634,9 +637,10 @@ void mlawNonLinearTVP::getMechSourceTVP(const STensor3& F0, const STensor3& F,
double mechSourceTVP = 0.;
if (this->getTimeStep() > 0){
mechSourceTVP += STensorOperation::doubledot(Me,Dp);
- mechSourceTVP -= STensorOperation::doubledot(tempX,alpha1);
- mechSourceTVP -= (R*Dgamma)/this->getTimeStep();
+ // mechSourceTVP -= STensorOperation::doubledot(tempX,alpha1);
+ // mechSourceTVP -= (R*Dgamma)/this->getTimeStep();
}
+ mechSourceTVP *= 0.8;
*Wm = mechSourceTVP;
}
@@ -688,27 +692,29 @@ void mlawNonLinearTVP::getMechSourceTVP(const STensor3& F0, const STensor3& F,
if (this->getTimeStep() > 0){
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
- dmechSourceTVPdF(i,j) = - (dRdF(i,j) - T*ddRdFdT(i,j))*Dgamma/this->getTimeStep() - R*DgammaDF(i,j)/this->getTimeStep();
+ // dmechSourceTVPdF(i,j) = - (dRdF(i,j) - T*ddRdFdT(i,j))*Dgamma/this->getTimeStep() - R*DgammaDF(i,j)/this->getTimeStep();
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dmechSourceTVPdF(i,j) += dMedF(i,j,k,l)*Dp(i,j) + Me(i,j)*dGammaNdF(i,j,k,l)/this->getTimeStep();
- dmechSourceTVPdF(i,j) += - (dXdF(i,j,k,l) - T*ddXdTdF(i,j,k,l))*alpha1(i,j) - tempX(i,j)*dAlphadF(i,j,k,l)/this->getTimeStep();
+ // dmechSourceTVPdF(i,j) += - (dXdF(i,j,k,l) - T*ddXdTdF(i,j,k,l))*alpha1(i,j) - tempX(i,j)*dAlphadF(i,j,k,l)/this->getTimeStep();
}
}
- }
+ }
+ dmechSourceTVPdF *= 0.8;
*dWmdF = dmechSourceTVPdF;
}
if(dWmdT!=NULL){
double dmechSourceTVPdT(0.);
if (this->getTimeStep() > 0){
- dmechSourceTVPdT += T*ddRdTT*Dgamma/this->getTimeStep() - R*DgammaDT/this->getTimeStep();
+ // dmechSourceTVPdT += T*ddRdTT*Dgamma/this->getTimeStep() - R*DgammaDT/this->getTimeStep();
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dmechSourceTVPdT += dMeDT(i,j)*Dp(i,j) + Me(i,j)*dGammaNdT(i,j)/this->getTimeStep();
- dmechSourceTVPdT += T*ddXdTT(i,j)*alpha1(i,j) - tempX(i,j)*dAlphadT(i,j)/this->getTimeStep();
+ // dmechSourceTVPdT += T*ddXdTT(i,j)*alpha1(i,j) - tempX(i,j)*dAlphadT(i,j)/this->getTimeStep();
}
}
+ dmechSourceTVPdT *= 0.8;
*dWmdT = dmechSourceTVPdT;
}
};
@@ -3283,11 +3289,14 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
q1->_elasticEnergy = mlawNonLinearTVE::freeEnergyMechanical(*q1,T0,T);
if (q1->_elasticEnergy > q0->_psiMax){
- q1->_psiMax = q1->_elasticEnergy;
- }
+ q1->_psiMax = q1->_elasticEnergy;
+ }
+ else if (q1->_elasticEnergy>q0->_elasticEnergy){
+ q1->_psiMax = q1->_elasticEnergy;
+ }
else{
- q1->_psiMax = q0->_psiMax;
- }
+ q1->_psiMax = q0->_psiMax;
+ }
// Msg::Error(" Inside TVP after correction, after updating q1->_elasticEnergy = %e!!",q1->_elasticEnergy);
diff --git a/NonLinearSolver/materialLaw/mullinsEffect.cpp b/NonLinearSolver/materialLaw/mullinsEffect.cpp
index 14bda0551..c2769d487 100644
--- a/NonLinearSolver/materialLaw/mullinsEffect.cpp
+++ b/NonLinearSolver/materialLaw/mullinsEffect.cpp
@@ -101,17 +101,58 @@ void negativeExponentialScaler::mullinsEffectScaling(double psi_max, const IPMu
// at unloading intitiation-> eta = 1.; unloading end-> eta = eta_min; reloading -> eta = eta_min;
double eta0 = ipvprev.getEta();
-
+ double _psi_max0 = ipvprev.getpsiMax();
+ double psi_cap = _psi;
double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.;
- if (_psi<psi_max && _psi<_psi0){
+ double _psi_max = std::max(_psi,_psi_max0);
+
+ double Term = sqrt(_m*(psi_max-eta*psi_cap)); // Term = 0 if _psi_max = _psi
+
+ if (Term!=0){
+
+ int ite = 0;
+ int maxite = 100;
+ double _tol = 1e-8;
+ double DfDeta(0.), Deta(0.);
+ double f = eta + _r*(1-exp(-Term)) - 1.;
+
+ while (fabs(f) >_tol or ite <1){
+
+ DfDeta = 1 + _r*(-_m*_psi*exp(-Term)/(2*Term));
+ Deta = - f/DfDeta;
+ if (eta + Deta <=0.){
+ eta /= 2.; // NR
+ }
+ else
+ eta += Deta;
+
+ Term = sqrt(_m*(psi_max-eta*psi_cap));
+ f = eta + _r*(1-exp(-Term)) - 1.;
+
+ if (fabs(f) <_tol) break;
+
+ if(ite > maxite){
+ Msg::Error("No convergence for eta in mullinsEffectScaling: ite= %d, f = %e!!",ite,f);
+ eta = sqrt(-1.);
+ return;
+ }
+ }
+ }
+
+ /*
+ if (_psi<psi_max){ // unloading
eta = 1. - _r*( 1. - exp( - sqrt(_m*(psi_max-_psi)) ) );
}
- else{
+ else if (_psi<_psi0){ // reinitiate eta after reloading, before unloading
eta = 1.;
}
+ else{ // reloading
+ eta = 1.;
+ }*/
ipv.set(eta,DetaDpsi,DDetaDpsipsi,DetaDT,DDetaDTT);
+ ipv.setPsiMax(psi_max);
}
void negativeExponentialScaler::mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, double T, IPMullinsEffect &ipv) const{
@@ -123,11 +164,14 @@ void negativeExponentialScaler::mullinsEffectScaling(double psi_max, const IPMul
double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.;
- if (_psi<psi_max && _psi<_psi0){
- eta = 1. - r*( 1. - exp( - sqrt(m*(psi_max-_psi)) ) );
+ if (_psi<psi_max){ // unloading
+ eta = 1. - _r*( 1. - exp( - sqrt(_m*(psi_max-_psi)) ) );
+ }
+ else if (_psi<_psi0){ // reinitiate eta after reloading, before unloading
+ eta = 1.;
}
- else{
- eta = eta0;
+ else{ // reloading
+ eta = 1.;
}
ipv.set(eta,DetaDpsi,DDetaDpsipsi,DetaDT,DDetaDTT);
--
GitLab