diff --git a/NonLinearSolver/internalPoints/CMakeLists.txt b/NonLinearSolver/internalPoints/CMakeLists.txt
index 174256314067b55dc46aedc208c330c14e29b92f..719159682cf4933391e62982bc429875a2b40f37 100644
--- a/NonLinearSolver/internalPoints/CMakeLists.txt
+++ b/NonLinearSolver/internalPoints/CMakeLists.txt
@@ -60,6 +60,7 @@ set(SRC
ipNonLinearTVM.cpp
ipNonLinearTVE.cpp
ipNonLinearTVP.cpp
+ ipMullinsEffect.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
diff --git a/NonLinearSolver/internalPoints/ipMullinsEffect.cpp b/NonLinearSolver/internalPoints/ipMullinsEffect.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5ae3728b9a0566549831a3e3747000b1f8b4dc1b
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipMullinsEffect.cpp
@@ -0,0 +1,45 @@
+//
+// Description: storing class for Mullins Effect
+//
+//
+// Author: <Ujwal Kishore J.>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "ipMullinsEffect.h"
+
+IPMullinsEffect::IPMullinsEffect(): eta(0), DetaDpsi(0), DDetaDpsipsi(0), DetaDT(0.), DDetaDTT(0.)
+{
+
+}
+
+IPMullinsEffect::IPMullinsEffect(const IPMullinsEffect &source)
+{
+ eta = source.eta;
+ DetaDpsi = source.DetaDpsi;
+ DDetaDpsipsi = source.DDetaDpsipsi;
+ DetaDT = source.DetaDT;
+ DDetaDTT = source.DDetaDTT;
+}
+
+IPMullinsEffect &IPMullinsEffect::operator=(const IPMullinsEffect &source)
+{
+ eta = source.eta;
+ DetaDpsi = source.DetaDpsi;
+ DDetaDpsipsi = source.DDetaDpsipsi;
+ DetaDT = source.DetaDT;
+ DDetaDTT = source.DDetaDTT;
+ return *this;
+}
+
+void IPMullinsEffect::restart()
+{
+ restartManager::restart(eta);
+ restartManager::restart(DetaDpsi);
+ restartManager::restart(DDetaDpsipsi);
+ restartManager::restart(DetaDT);
+ restartManager::restart(DDetaDTT);
+ return;
+}
diff --git a/NonLinearSolver/internalPoints/ipMullinsEffect.h b/NonLinearSolver/internalPoints/ipMullinsEffect.h
new file mode 100644
index 0000000000000000000000000000000000000000..153ce88e3afeeee2b7948d67ed85b568340a4e81
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipMullinsEffect.h
@@ -0,0 +1,50 @@
+//
+// Description: storing class for Mullins Effect
+//
+//
+// Author: <Ujwal Kishore J.>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef IPMULLINSEFFECT_H_
+#define IPMULLINSEFFECT_H_
+
+#include "ipHardening.h"
+
+class IPMullinsEffect{
+
+ protected:
+ double eta; // scaling variable
+ double DetaDpsi;
+ double DDetaDpsipsi;
+ double DetaDT;
+ double DDetaDTT;
+
+ public:
+ IPMullinsEffect();
+ IPMullinsEffect(const IPMullinsEffect &source);
+ virtual IPMullinsEffect &operator=(const IPMullinsEffect &source);
+ virtual ~IPMullinsEffect(){}
+ virtual void restart();
+
+ virtual IPMullinsEffect * clone() const {return new IPMullinsEffect(*this);};
+
+ virtual double getEta() const {return eta;}
+ virtual double getDetaDpsi() const {return DetaDpsi;}
+ virtual double getDDetaDpsipsi() const {return DDetaDpsipsi;};
+ virtual double getDetaDT() const{ return DetaDT;};
+ virtual double getDDetaDTT() const{ return DDetaDTT;};
+ virtual void set(const double &_eta, const double &_DetaDpsi, const double& _DDetaDpsipsi, const double& _DetaDT, const double& _DDetaDTT)
+ {
+ eta=_eta;
+ DetaDpsi=_DetaDpsi;
+ DDetaDpsipsi = _DDetaDpsipsi;
+ DetaDT = _DetaDT;
+ DDetaDTT = _DDetaDTT;
+ }
+};
+
+
+#endif // IPMULLINSEFFECT_H_
diff --git a/NonLinearSolver/materialLaw/CMakeLists.txt b/NonLinearSolver/materialLaw/CMakeLists.txt
index 9da69791fd17edb91f649091ca268300272dc4f0..3e2cda868ea7b93d773217f4527863a57b99ee3d 100644
--- a/NonLinearSolver/materialLaw/CMakeLists.txt
+++ b/NonLinearSolver/materialLaw/CMakeLists.txt
@@ -99,6 +99,7 @@ set(SRC
mlawNonLinearTVM.cpp
mlawNonLinearTVE.cpp
mlawNonLinearTVP.cpp
+ mullinsEffect.cpp
# src/op_eshelby.cpp
# Headers without cpp change this ??
ID.h
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
index 379f0dbfb6575aaacc7ac58c8d276ed763be27eb..ca9ba845f2927e5a77f3762ad2110c97756ac202 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
@@ -525,10 +525,10 @@ double mlawNonLinearTVE::setTemp(const double T0, const double T1, const int par
return T_set;
}
-double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const double Tc) const{
+double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const double T) const{
// Update the Properties to the current temperature (see below which ones are being used)
- double KT, GT, AlphaT; getK(KT,Tc); getG(GT,Tc); getAlpha(AlphaT,Tc);
+ double KT, GT, AlphaT; getK(&q,KT,T); getG(&q,GT,T); getAlpha(AlphaT,T);
double Psy_mech = 0.;
@@ -538,7 +538,7 @@ double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const dou
STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
// Get thermal strain
- double Eth = 3.*AlphaT*(Tc-_Tinitial);
+ double Eth = 3.*AlphaT*(T-_Tinitial);
// Get Equilibrium freeEnergy_mech
Psy_mech = KT*0.5*(trEe-Eth)*(trEe-Eth) + GT*STensorOperation::doubledot(devEe,devEe);
@@ -546,7 +546,7 @@ double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const dou
// Add Branch freeEnergy_mech
if ((_Ki.size() > 0) or (_Gi.size() > 0)){
- std::vector<double> KiT,GiT,AlphaiT; KiT.resize(_N,0.); GiT.resize(_N,0.); AlphaiT.resize(_N,0.); getKi(KiT,Tc); getGi(GiT,Tc); getAlphai(AlphaiT,Tc);
+ std::vector<double> KiT,GiT,AlphaiT; KiT.resize(_N,0.); GiT.resize(_N,0.); AlphaiT.resize(_N,0.); getKi(&q,KiT,T); getGi(&q,GiT,T); getAlphai(AlphaiT,T);
for (int i=0; i<_N; i++){
static STensor3 devEebranch;
@@ -560,7 +560,7 @@ double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const dou
devEebranch -= devGamma;
if (_modelFlag == 2){ // Additional volumetric terms
- trGamma += 3*(AlphaT-AlphaiT[i])*(Tc-_Tinitial);
+ trGamma += 3*(AlphaT-AlphaiT[i])*(T-_Tinitial);
}
Psy_mech += KiT[i]*0.5*(trEe-trGamma)*(trEe-trGamma) + GiT[i]*STensorOperation::doubledot(devEebranch,devEebranch);
@@ -570,9 +570,9 @@ double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const dou
return Psy_mech;
}
-void mlawNonLinearTVE::corKirInfinity( const STensor3& devEe, const double& trEe, const double T1, STensor3& CorKirDevInf, double& CorKirTrInf) const{
+void mlawNonLinearTVE::corKirInfinity(const IPNonLinearTVE *q1, const STensor3& devEe, const double& trEe, const double T1, STensor3& CorKirDevInf, double& CorKirTrInf) const{
double KT, GT, AlphaT;
- getK(KT,T1); getG(GT,T1); getAlpha(AlphaT,T1);
+ getK(q1,KT,T1); getG(q1,GT,T1); getAlpha(AlphaT,T1);
CorKirDevInf = 2*GT*devEe;
CorKirTrInf = KT*(trEe - 3*AlphaT*(T1-_Tinitial));
}
@@ -587,7 +587,7 @@ void mlawNonLinearTVE::getTVEdCorKirDT(const IPNonLinearTVE *q0, IPNonLinearTVE
// Update the Properties to the current temperature (see below which ones are being used)
double KT, GT, AlphaT, DKDT, DGDT, DAlphaDT, DDKDT, DDGDT, DDAlphaDT;
- getK(KT,T,&DKDT,&DDKDT); getG(GT,T,&DGDT,&DDGDT); getAlpha(AlphaT,T,&DAlphaDT,&DDAlphaDT);
+ getK(q1,KT,T,&DKDT,&DDKDT); getG(q1,GT,T,&DGDT,&DDGDT); getAlpha(AlphaT,T,&DAlphaDT,&DDAlphaDT);
// Add corKirinf terms to DcorKirDT
STensor3 corKirDevInf, devEe;
@@ -604,7 +604,7 @@ void mlawNonLinearTVE::getTVEdCorKirDT(const IPNonLinearTVE *q0, IPNonLinearTVE
STensorOperation::decomposeDevTr(DEe,devDE,trDE);
// Add Infinity terms
- corKirInfinity( devEe, trEe, T, corKirDevInf, corKirTrInf );
+ corKirInfinity(q1, devEe, trEe, T, corKirDevInf, corKirTrInf );
for (int i=0; i<3; i++){
DcorKirDT(i,i) += corKirTrInf*DKDT/KT - 3*KT*(DAlphaDT*(T-_Tinitial) + AlphaT);
DDcorKirDT(i,i) += corKirTrInf*DDKDT/KT - 6*DKDT*(DAlphaDT*(T-_Tinitial) + AlphaT)- 3*KT*(DDAlphaDT*(T-_Tinitial) + 2*DAlphaDT);
@@ -620,6 +620,7 @@ void mlawNonLinearTVE::getTVEdCorKirDT(const IPNonLinearTVE *q0, IPNonLinearTVE
// Update the Branch Moduli to the current temperature
std::vector<double> KiT, GiT, DKiDT, DGiDT, DDKiDT, DDGiDT;
KiT.resize(_N,0.);GiT.resize(_N,0.);DKiDT.resize(_N,0.);DGiDT.resize(_N,0.);DDKiDT.resize(_N,0.);DDGiDT.resize(_N,0.);
+ getKi(q1,KiT,T); getGi(q1,GiT,T);
static STensor3 DOiDT;
static STensor3 DDOiDT;
@@ -672,10 +673,10 @@ void mlawNonLinearTVE::getTVEdCorKirDT(const IPNonLinearTVE *q0, IPNonLinearTVE
STensorOperation::zero(DOGammai); STensorOperation::zero(DOGammaiDT); STensorOperation::zero(DDOGammaiDT);
if (_TemFuncOpt == 0){ // Use Shift Factor
- DKiDT[i] = -_Ki[i]/(_ki[i]) * Ddt_shiftDT*expdtkby2; // Removed _ki/2 and _gi/2
- DDKiDT[i] = -_Ki[i]/(_ki[i]) * (DDdt_shiftDT - 1/(_ki[i])*pow(Ddt_shiftDT,2))*expdtkby2;
- DGiDT[i] = -_Gi[i]/(_gi[i]) * Ddt_shiftDT*expdtgby2;
- DDGiDT[i] = -_Gi[i]/(_gi[i]) * (DDdt_shiftDT - 1/(_gi[i])*pow(Ddt_shiftDT,2))*expdtgby2;
+ DKiDT[i] = - KiT[i]/(_ki[i]) * Ddt_shiftDT*expdtkby2; // Removed _ki/2 and _gi/2
+ DDKiDT[i] = - KiT[i]/(_ki[i]) * (DDdt_shiftDT - 1/(_ki[i])*pow(Ddt_shiftDT,2))*expdtkby2;
+ DGiDT[i] = - GiT[i]/(_gi[i]) * Ddt_shiftDT*expdtgby2;
+ DDGiDT[i] = - GiT[i]/(_gi[i]) * (DDdt_shiftDT - 1/(_gi[i])*pow(Ddt_shiftDT,2))*expdtgby2;
}
for (int j=0; j<3; j++){
@@ -739,12 +740,11 @@ void mlawNonLinearTVE::getMechSourceTVE(const IPNonLinearTVE *q0, IPNonLinearTVE
const double& DKDTsum, const double& DGDTsum, const STensor43& DEeDFe,
double *Wm, STensor3 *dWmdF, double *dWmdT) const{
- // I am stupid
+ // Initialize
double mechSrcTVE, dmechSrcTVEdT = 0.;
STensor3 dmechSrcTVEdE(0.);
STensor3 dmechSrcTVEdF(0.);
-
// Make Identity Tensors
static const STensor43 I4(1.,1.);
@@ -772,7 +772,8 @@ void mlawNonLinearTVE::getMechSourceTVE(const IPNonLinearTVE *q0, IPNonLinearTVE
// Update the Branch Moduli to the current temperature
std::vector<double> KiT, GiT, DKiDT, DGiDT, DDKiDT, DDGiDT;
KiT.resize(_N,0.);GiT.resize(_N,0.);DKiDT.resize(_N,0.);DGiDT.resize(_N,0.);DDKiDT.resize(_N,0.);DDGiDT.resize(_N,0.);
-
+ getKi(q1,KiT,T); getGi(q1,GiT,T);
+
static STensor3 DOiDT;
static STensor3 DDOiDT;
static STensor3 DOGammaiDT;
@@ -824,11 +825,10 @@ void mlawNonLinearTVE::getMechSourceTVE(const IPNonLinearTVE *q0, IPNonLinearTVE
STensorOperation::zero(DOGammai); STensorOperation::zero(DOGammaiDT); STensorOperation::zero(DDOGammaiDT);
if (_TemFuncOpt == 0){ // Use Shift Factor
- getKi(KiT,T); getGi(GiT,T);
- DKiDT[i] = -_Ki[i]/(_ki[i]) * Ddt_shiftDT*expdtkby2; // Removed _ki/2 and _gi/2
- DDKiDT[i] = -_Ki[i]/(_ki[i]) * (DDdt_shiftDT - 1/(_ki[i])*pow(Ddt_shiftDT,2))*expdtkby2;
- DGiDT[i] = -_Gi[i]/(_gi[i]) * Ddt_shiftDT*expdtgby2;
- DDGiDT[i] = -_Gi[i]/(_gi[i]) * (DDdt_shiftDT - 1/(_gi[i])*pow(Ddt_shiftDT,2))*expdtgby2;
+ DKiDT[i] = - KiT[i]/(_ki[i]) * Ddt_shiftDT*expdtkby2; // Removed _ki/2 and _gi/2
+ DDKiDT[i] = - KiT[i]/(_ki[i]) * (DDdt_shiftDT - 1/(_ki[i])*pow(Ddt_shiftDT,2))*expdtkby2;
+ DGiDT[i] = - GiT[i]/(_gi[i]) * Ddt_shiftDT*expdtgby2;
+ DDGiDT[i] = - GiT[i]/(_gi[i]) * (DDdt_shiftDT - 1/(_gi[i])*pow(Ddt_shiftDT,2))*expdtgby2;
}
for (int j=0; j<3; j++){
@@ -878,9 +878,9 @@ void mlawNonLinearTVE::getMechSourceTVE(const IPNonLinearTVE *q0, IPNonLinearTVE
STensorOperation::zero(DGammai); STensorOperation::zero(DDGammaiDT);
DQiDT -= DOiDT; DDQiDT -= DDOiDT; DDGammaiDT -= DDOGammaiDT;
for (int j=0; j<3; j++){
- Qi(j,j) += (_Ki[i]*trEe - trOi[i]);
+ Qi(j,j) += (KiT[i]*trEe - trOi[i]);
for (int k=0; k<3; k++){
- Qi(j,k) += (2*_Gi[i]*devEe(j,k) - devOi[i](j,k));
+ Qi(j,k) += (2*GiT[i]*devEe(j,k) - devOi[i](j,k));
}
}
DGammai += DEe;
@@ -893,8 +893,8 @@ void mlawNonLinearTVE::getMechSourceTVE(const IPNonLinearTVE *q0, IPNonLinearTVE
static STensor43 temp_DDQiDTDE;
static STensor43 temp_DDgammiDtDE;
- double Ki_term = (_Ki[i]*(1-expdtkby2)+T*DKiDT[i]);
- double Gi_term = (_Gi[i]*(1-expdtgby2)+T*DGiDT[i]);
+ double Ki_term = (KiT[i]*(1-expdtkby2)+T*DKiDT[i]);
+ double Gi_term = (GiT[i]*(1-expdtgby2)+T*DGiDT[i]);
isotropicHookTensor( Ki_term, Gi_term , temp_DDQiDTDE );
double Ki_term_2 = 1/3*expdtkby2;
@@ -1000,7 +1000,7 @@ void mlawNonLinearTVE::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
// Update the Properties to the current temperature (see below which ones are being used)
double KT, GT, AlphaT_0, AlphaT, DKDT, DGDT, DAlphaDT;
- getK(KT,T_set,&DKDT); getG(GT,T_set,&DGDT); getAlpha(AlphaT_0,T0); getAlpha(AlphaT,T_set,&DAlphaDT);
+ getK(q1,KT,T_set,&DKDT); getG(q1,GT,T_set,&DGDT); getAlpha(AlphaT_0,T0); getAlpha(AlphaT,T_set,&DAlphaDT);
// Get LogStrain
static STensor3 DE, devDE, devEe;
@@ -1046,7 +1046,7 @@ void mlawNonLinearTVE::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
DKiDT.resize(_N,0.); DGiDT.resize(_N,0.); DAlphaiDT.resize(_N,0.); DDKiDT.resize(_N,0.); DDGiDT.resize(_N,0.);
if (_TemFuncOpt == 1){ // Use Scalar Functions
- getKi(KiT,T_set,&DKiDT,&DDKiDT); getGi(GiT,T_set,&DGiDT,&DDGiDT); getAlphai(AlphaiT_0,T0); getAlphai(AlphaiT,T_set,&DAlphaiDT);
+ getKi(q1,KiT,T_set,&DKiDT,&DDKiDT); getGi(q1,GiT,T_set,&DGiDT,&DDGiDT); getAlphai(AlphaiT_0,T0); getAlphai(AlphaiT,T_set,&DAlphaiDT);
}
@@ -1064,9 +1064,9 @@ void mlawNonLinearTVE::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
// Maxwell
if (_TemFuncOpt == 0){ // Use Shift Factor
- getGi(GiT,T_set); // Removed 2*_gi[i]
- DGiDT[i] = -_Gi[i]/(_gi[i]) * Ddt_shiftDT*expdtgby2;
- DDGiDT[i] = -_Gi[i]/(_gi[i]) * (DDdt_shiftDT - 1/(_gi[i])*pow(Ddt_shiftDT,2))*expdtgby2;
+ getGi(q1,GiT,T_set); // Removed 2*_gi[i]
+ DGiDT[i] = -GiT[i]/(_gi[i]) * Ddt_shiftDT*expdtgby2;
+ DDGiDT[i] = -GiT[i]/(_gi[i]) * (DDdt_shiftDT - 1/(_gi[i])*pow(Ddt_shiftDT,2))*expdtgby2;
}
// For Single Convolution
@@ -1151,9 +1151,9 @@ void mlawNonLinearTVE::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
// Maxwell
if (_TemFuncOpt == 0){ // Use Shift Factor
- getKi(KiT,T_set); getAlphai(AlphaiT,T_set); // Removed 2*_ki[i]
- DKiDT[i] = -_Ki[i]/(_ki[i]) * Ddt_shiftDT *expdtkby2;
- DDKiDT[i] = -_Ki[i]/(_ki[i]) * (DDdt_shiftDT - 1/(_ki[i])*pow(Ddt_shiftDT,2)) *expdtkby2;
+ getKi(q1,KiT,T_set); getAlphai(AlphaiT,T_set); // Removed 2*_ki[i]
+ DKiDT[i] = -KiT[i]/(_ki[i]) * Ddt_shiftDT *expdtkby2;
+ DDKiDT[i] = -KiT[i]/(_ki[i]) * (DDdt_shiftDT - 1/(_ki[i])*pow(Ddt_shiftDT,2)) *expdtkby2;
}
// Single Convolution
@@ -1262,14 +1262,13 @@ void mlawNonLinearTVE::constitutive(const STensor3& F0,
const IPNonLinearTVE *q0=dynamic_cast<const IPNonLinearTVE *>(q0i);
IPNonLinearTVE *q1 = dynamic_cast<IPNonLinearTVE *>(q1i);
- double Tc = _Tref;
+ double T = _Tref;
static SVector3 gradT, temp2;
static STensor3 temp3;
static STensor33 temp33;
static double tempVal;
- static STensor43 dFpdF, dFedF;
- static STensor3 dFpdT, dFedT, dGammadF;
- double dGammadT;
+ static STensor43 dFedF;
+ static STensor3 dFedT;
if(elasticTangent!=NULL) Msg::Error("mlawNonLinearTVE elasticTangent not defined");
predictorCorrector_ThermoViscoElastic(F0,F,P,q0,q1,Tangent,dFedF,dFedT,_Tref,_Tref,gradT,gradT,temp2,temp3,temp3,temp2,temp33,tempVal,
@@ -1431,7 +1430,7 @@ void mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(
// Update the Properties to the current temperature (see below which ones are being used)
double KT, GT, AlphaT, KThConT, CpT, DKDT, DGDT, DAlphaDT, DKThConDT, DCpDT, DDKDT, DDGDT, DDAlphaDT;
- getK(KT,T_set,&DKDT,&DDKDT); getG(GT,T_set,&DGDT,&DDGDT); getAlpha(AlphaT,T_set,&DAlphaDT,&DDAlphaDT); getKTHCon(KThConT,T_set,&DKThConDT); getCp(CpT,T_set,&DCpDT);
+ getK(q1,KT,T_set,&DKDT,&DDKDT); getG(q1,GT,T_set,&DGDT,&DDGDT); getAlpha(AlphaT,T_set,&DAlphaDT,&DDAlphaDT); getKTHCon(KThConT,T_set,&DKThConDT); getCp(CpT,T_set,&DCpDT);
// Calculate/update Log Strain
static STensor3 C;
@@ -1520,7 +1519,7 @@ void mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(
double& trDE = q1->getRefToTrDE();
STensorOperation::decomposeDevTr(DEe,devDE,trDE);
- corKirInfinity( devEe, trEe, T, corKirDevInf, corKirTrInf );
+ corKirInfinity( q1, devEe, trEe, T, corKirDevInf, corKirTrInf );
for (int i=0; i<3; i++){
// DcorKirDT(i,i) += corKirTr/3/(Ke+Kde*(T-T_temp))*(DKe + DKde*(T-T_temp)) - 3*KT*(DAlphaDT*(T-_Tinitial) + AlphaT); // 1st formulation
@@ -1541,7 +1540,7 @@ void mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(
KiT.resize(_N,0.);GiT.resize(_N,0.);AlphaiT.resize(_N,0.);DKiDT.resize(_N,0.);DGiDT.resize(_N,0.);DAlphaiDT.resize(_N,0.);DDKiDT.resize(_N,0.);DDGiDT.resize(_N,0.);DDAlphaiDT.resize(_N,0.);
if (_TemFuncOpt == 1){
- getKi(KiT,T_set,&DKiDT,&DDKiDT); getGi(GiT,T_set,&DGiDT,&DDGiDT); getAlphai(AlphaiT,T_set,&DAlphaiDT,&DDAlphaiDT);
+ getKi(q1,KiT,T_set,&DKiDT,&DDKiDT); getGi(q1,GiT,T_set,&DGiDT,&DDGiDT); getAlphai(AlphaiT,T_set,&DAlphaiDT,&DDAlphaiDT);
}
static STensor3 DOiDT;
@@ -1600,11 +1599,11 @@ void mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(
STensorOperation::zero(DOGammai); STensorOperation::zero(DOGammaiDT); STensorOperation::zero(DDOGammaiDT);
if (_TemFuncOpt == 0){ // Use Shift Factor
- getKi(KiT,T_set); getGi(GiT,T_set); getAlphai(AlphaiT,T_set);
- DKiDT[i] = -_Ki[i]/(_ki[i]) * Ddt_shiftDT*expdtkby2; // Removed _ki/2 and _gi/2
- DDKiDT[i] = -_Ki[i]/(_ki[i]) * (DDdt_shiftDT - 1/(_ki[i])*pow(Ddt_shiftDT,2))*expdtkby2;
- DGiDT[i] = -_Gi[i]/(_gi[i]) * Ddt_shiftDT*expdtgby2;
- DDGiDT[i] = -_Gi[i]/(_gi[i]) * (DDdt_shiftDT - 1/(_gi[i])*pow(Ddt_shiftDT,2))*expdtgby2;
+ getKi(q1,KiT,T_set); getGi(q1,GiT,T_set); getAlphai(AlphaiT,T_set);
+ DKiDT[i] = - KiT[i]/(_ki[i]) * Ddt_shiftDT*expdtkby2; // Removed _ki/2 and _gi/2
+ DDKiDT[i] = - KiT[i]/(_ki[i]) * (DDdt_shiftDT - 1/(_ki[i])*pow(Ddt_shiftDT,2))*expdtkby2;
+ DGiDT[i] = - GiT[i]/(_gi[i]) * Ddt_shiftDT*expdtgby2;
+ DDGiDT[i] = - GiT[i]/(_gi[i]) * (DDdt_shiftDT - 1/(_gi[i])*pow(Ddt_shiftDT,2))*expdtgby2;
}
for (int j=0; j<3; j++){
@@ -1660,9 +1659,9 @@ void mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(
STensorOperation::zero(DGammai); STensorOperation::zero(DDGammaiDT);
DQiDT -= DOiDT; DDQiDT -= DDOiDT; DDGammaiDT -= DDOGammaiDT;
for (int j=0; j<3; j++){
- Qi(j,j) += (_Ki[i]*trEe - trOi[i]);
+ Qi(j,j) += (KiT[i]*trEe - trOi[i]);
for (int k=0; k<3; k++){
- Qi(j,k) += (2*_Gi[i]*devEe(j,k) - devOi[i](j,k));
+ Qi(j,k) += (2*GiT[i]*devEe(j,k) - devOi[i](j,k));
}
}
DGammai += DEe;
@@ -1678,8 +1677,8 @@ void mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(
static STensor43 temp_DDQiDTDE;
static STensor43 temp_DDgammiDtDE;
- double Ki_term = (_Ki[i]*(1-expdtkby2)+T*DKiDT[i]);
- double Gi_term = (_Gi[i]*(1-expdtgby2)+T*DGiDT[i]);
+ double Ki_term = (KiT[i]*(1-expdtkby2)+T*DKiDT[i]);
+ double Gi_term = (GiT[i]*(1-expdtgby2)+T*DGiDT[i]);
isotropicHookTensor( Ki_term, Gi_term , temp_DDQiDTDE );
double Ki_term_2 = 1/3*expdtkby2;
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.h b/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
index 9afdddbf5a9a02afeeeb2245871a79f064fa7457..c8eed1b2bd409aed251beddcdc1018920e736bc3 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
@@ -22,49 +22,6 @@
#include "scalarFunction.h"
#include "mlawHyperelastic.h"
-/* Patch Notes
-
-1) 23/12/21 - NonLinearTVE adopted from mlawHyperelastic; Pure ViscoElasticity with Temp Functions for Moduli and Thermal Props (latter Still Unused)
- [mLAW FORMED -- .h , .cpp]
- Changed mu to G
-.......1.0.1) 24/12/21 - Created new IPVariable, updated dG3DMaterialLaw, update dG3DIPVariable, updated mlaw.h,
-.......1.0.2) 28/12/21 - Added the mlaw.cpp and ip.cpp to the CMakeFiles.txt in materiallaw/ and internalpoints/ respectively.
-.......1.0.3) 28/12/21 - (Tested against previous results of mlawHyperelastic in Octave const Temp results --- Works fine!!).
-
-...1.1) 29/12/21 - Updated with ThermoMechanics - (UNTESTED so far).
-.......1.1.1) 29/12/21 - Updated dG3DIPVariable with inheritance from ThermoMechanicsDG3DIPVariableBase Pure virtual class (and LinearThermoMechanicsDG3DIPVariable).
- - Noted that the Pure virtual base class may cause problems when introducing nonLocal Variables.
-.......1.1.2) 29/12/21 - Updated dG3DMaterialLaw with adaptations from LinearThermoMechanicsDG3DMaterial.
-.......1.1.3) 29/12/21 - Updated ipNonLinearTVE with adaptations from LinearThermoMechanicsDG3DMaterial.
-.......1.1.4) 29/12/21 - Updated mlawNonLinearTVE with adaptations from LinearThermoMechanicsDG3DMaterial.
-
-...1.2) 30/12/21 - Add Temperature BCs, add temperature extradof from dG3DIpvariable (ThermoMechanicsDG3DIPVariableBase),
- change the Stress calculation Function in dG3DMaterialLaw;
-.......1.2.1) 05/01/22 - linearK and stiff_alphadilatation added to ipvcur pointer in dG3Dmateriallaw.cpp through direct member functions in mlaw....TVE.h
- - IC Working, ScalarFluxBC Not Working!!!!!!!!!
-.......1.2.2) 06/01/22 - CleanUp!
- - Removed thermal props unused member functions from everywhere.
- - Remember to add tempFunc multiplication where needed! Note that my Tg scalarFunctions will likely not be normalised functions.
-.......1.2.3) 06/01/22 - Attempted addition of coupled ThermoMech features from ThermoElastic part of FullJ2ThermoMech...dG3DMaterialLaw
-
-...1.3) 15/02/22 - Addition of TTSP
-.......1.3.1) 22/02/22 - Addition of WLF shift factor and shifted time member functions
-.......1.3.2) 04/03/22 - Addition of convolution single/double integrals in the predictor functions
- - CLEANUP 1.0 -> Making mlaw and ipvariable dependent classes of mlawHyperElastic
-.......1.3.3) 14/03/22 - Thermal and mechanical source code
-.......1.3.4) 24/03/22 - Tangent and mechanical source derivatives
-.......1.3.5) 13/04/22 - CLEANUP 2.0
-.......1.3.6) 25/04/22 - Fully Corrected TVE derivatives.
-...1.4) 03/06/22 - CLEANUP 3.0 Renamed class to mlawNonLinearTVE
-.......1.4.1) 04/06/22 - Removed all Unused Code
-.......1.4.2) 07/06/22 - Added calculateStress boolean in ThermoViscoElasticPredictor -> to just update TVE flow in TVP predictorCorrector functions
- *
-.......1.?.?) ??/??/22 - Addition of new class for Psi_ve
- *
-
-2) 04/06/22 - NonLinearTVP adopted from mlawPowerYieldHyper; Temp Functions for Hardening Moduli and Yield Strengths
- [mLAW FORMED -- .h , .cpp, in the mlawNonLinearTVE]
-*/
class mlawNonLinearTVE : public mlawPowerYieldHyper{
protected:
@@ -133,10 +90,10 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
double *t_shift = NULL, double *Ddt_shiftDT = NULL, double *DDdt_shiftDT = NULL,
const int opt_time = 1, const int opt_order = 2) const;
virtual double linearInterpTemp(const double t_interp, const double t0, const double t1, const double T0, const double T1) const;
- virtual double freeEnergyMechanical(const IPNonLinearTVE& q, const double Tc) const ;
+ virtual double freeEnergyMechanical(const IPNonLinearTVE& q, const double T) const ;
virtual double setTemp(const double T0, const double T1, const int para) const ;
- void corKirInfinity( const STensor3& devEe, const double& trEe, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const;
+ void corKirInfinity(const IPNonLinearTVE *q1, const STensor3& devEe, const double& trEe, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const;
void ThermoViscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
const IPNonLinearTVE *q0, IPNonLinearTVE *q1,
double& Ke, double& Ge, double& Kde, double& Gde,
@@ -274,112 +231,128 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
return Beta_i;
}
- virtual void getK(double& KT, const double Tc, double *dK = NULL, double *ddK = NULL) const{
- KT = _K*_temFunc_K->getVal(Tc);
+ virtual void getK(const IPNonLinearTVE *q1, double& KT, const double T, double *dK = NULL, double *ddK = NULL) const{
+ double K = _K*q1->getConstRefToElasticBulkPropertyScaleFactor();
+ KT = _K*_temFunc_K->getVal(T);
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(T); // 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(T);
}
}
- virtual void getG(double& GT, const double Tc, double *dG = NULL, double *ddG = NULL) const{
- GT = _G*_temFunc_G->getVal(Tc);
+ virtual void getG(const IPNonLinearTVE *q1, double& GT, const double T, double *dG = NULL, double *ddG = NULL) const{
+ double G = _G*q1->getConstRefToElasticShearPropertyScaleFactor();
+ GT = _G*_temFunc_G->getVal(T);
if(dG!=NULL){
- *dG = _G*_temFunc_G->getDiff(Tc);
+ *dG = _G*_temFunc_G->getDiff(T);
}
if(ddG!=NULL){
- *ddG = _G*_temFunc_G->getDoubleDiff(Tc);
+ *ddG = _G*_temFunc_G->getDoubleDiff(T);
}
}
- virtual void getAlpha(double& alphaT, const double Tc, double *dAlpha = NULL, double *ddAlpha = NULL) const{
- alphaT = _scalarAlpha*_temFunc_Alpha->getVal(Tc);
+ virtual void getAlpha(double& alphaT, const double T, double *dAlpha = NULL, double *ddAlpha = NULL) const{
+ alphaT = _scalarAlpha*_temFunc_Alpha->getVal(T);
if(dAlpha!=NULL){
- *dAlpha = _scalarAlpha*_temFunc_Alpha->getDiff(Tc);
+ *dAlpha = _scalarAlpha*_temFunc_Alpha->getDiff(T);
}
if(ddAlpha!=NULL){
- *ddAlpha = _scalarAlpha*_temFunc_Alpha->getDoubleDiff(Tc);
+ *ddAlpha = _scalarAlpha*_temFunc_Alpha->getDoubleDiff(T);
}
}
- virtual void getKTHCon(double& KThConT,const double Tc, double *dKThCon = NULL) const{
- KThConT = _scalarK*_temFunc_KThCon->getVal(Tc);
+ virtual void getKTHCon(double& KThConT,const double T, double *dKThCon = NULL) const{
+ KThConT = _scalarK*_temFunc_KThCon->getVal(T);
if(dKThCon!=NULL){
- *dKThCon = _scalarK*_temFunc_KThCon->getDiff(Tc);
+ *dKThCon = _scalarK*_temFunc_KThCon->getDiff(T);
}
}
- virtual void getCp(double& CpT, const double Tc, double *dCp = NULL, double *ddCp = NULL) const{
- CpT = _Cp*_temFunc_Cp->getVal(Tc);
+ virtual void getCp(double& CpT, const double T, double *dCp = NULL, double *ddCp = NULL) const{
+ CpT = _Cp*_temFunc_Cp->getVal(T);
if(dCp!=NULL){
- *dCp = _Cp*_temFunc_Cp->getDiff(Tc);
+ *dCp = _Cp*_temFunc_Cp->getDiff(T);
}
if(ddCp!=NULL){
- *ddCp = _Cp*_temFunc_Cp->getDoubleDiff(Tc);
+ *ddCp = _Cp*_temFunc_Cp->getDoubleDiff(T);
}
}
// The 3rd argument is a pointer to a vector.
- virtual void getKi(std::vector<double>& Ki_T, const double Tc, std::vector<double>* dKi=NULL, std::vector<double>* ddKi=NULL) const{
+ virtual void getKi(const IPNonLinearTVE *q1, std::vector<double>& Ki_T, const double T, std::vector<double>* dKi=NULL, std::vector<double>* ddKi=NULL) const{
+
+ std::vector<double> Ki;
+ Ki.clear(); Ki.resize(_N,0.);
+ for (int i=0; i<_Ki.size();i++){
+ Ki[i] = _Ki[i]*q1->getConstRefToElasticBulkPropertyScaleFactor();
+ }
+
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(T);
if(dKi!=NULL){
- dKi->at(i) = _Ki[i]*_temFunc_Ki[i]->getDiff(Tc);
+ dKi->at(i) = _Ki[i]*_temFunc_Ki[i]->getDiff(T);
}
if(ddKi!=NULL){
- ddKi->at(i) = _Ki[i]*_temFunc_Ki[i]->getDoubleDiff(Tc);
+ ddKi->at(i) = _Ki[i]*_temFunc_Ki[i]->getDoubleDiff(T);
}
}
}
// The 3rd argument is a pointer to a vector.
- virtual void getGi(std::vector<double>& Gi_T, const double Tc, std::vector<double>* dGi=NULL, std::vector<double>* ddGi=NULL) const{
+ virtual void getGi(const IPNonLinearTVE *q1, std::vector<double>& Gi_T, const double T, std::vector<double>* dGi=NULL, std::vector<double>* ddGi=NULL) const{
+
+ 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(T);
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(T); // *(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(T);
}
}
}
// The 3rd argument is a pointer to a vector.
- virtual void getAlphai(std::vector<double>& Alphai_T, const double Tc, std::vector<double>* dAlphai=NULL, std::vector<double>* ddAlphai=NULL) const{
+ virtual void getAlphai(std::vector<double>& Alphai_T, const double T, std::vector<double>* dAlphai=NULL, std::vector<double>* ddAlphai=NULL) const{
Alphai_T.resize(_N,0.);
for (int i=0; i<_Alphai.size();i++){
- Alphai_T[i] = _Alphai[i]*_temFunc_Alphai[i]->getVal(Tc);
+ Alphai_T[i] = _Alphai[i]*_temFunc_Alphai[i]->getVal(T);
if(dAlphai!=NULL){
- dAlphai->at(i) = _Alphai[i]*_temFunc_Alphai[i]->getDiff(Tc); // *(dAlphai+i) is equivalent to dAlphai[i]
+ dAlphai->at(i) = _Alphai[i]*_temFunc_Alphai[i]->getDiff(T); // *(dAlphai+i) is equivalent to dAlphai[i]
}
if(ddAlphai!=NULL){
- ddAlphai->at(i) = _Alphai[i]*_temFunc_Alphai[i]->getDoubleDiff(Tc);
+ ddAlphai->at(i) = _Alphai[i]*_temFunc_Alphai[i]->getDoubleDiff(T);
}
}
}
// Deprecated member function - Currently implementing scalar functions directly in the .py file.
- virtual void temFuncGlassTransition(const double Tc, const double opt, const double Prop0, const double Prop1, double& Prop, double *dProp = NULL) const {
+ virtual void temFuncGlassTransition(const double T, const double opt, const double Prop0, const double Prop1, double& Prop, double *dProp = NULL) const {
// Prop0 = At room temp or glassy state
// Prop1 = In rubbery state - or the lowest possible value after glass transition
if (opt == 0) { // LinearFunction - This is only for test, Prop1 = 0;
- Prop = Prop0*(1 - 0.1/Prop0*(Tc - 298));
+ Prop = Prop0*(1 - 0.1/Prop0*(T - 298));
if (dProp!= NULL) {
*dProp = -0.1;
}
} else if (opt == 1) { // Gibson et al
double Tg = 373; double k = 0.0975; double m = 0;
- Prop = 0.5*(Prop0+Prop1) -0.5*(Prop0-Prop1) * tanh(k*(Tc-Tg)) - m*(Tc-Tg);
+ Prop = 0.5*(Prop0+Prop1) -0.5*(Prop0-Prop1) * tanh(k*(T-Tg)) - m*(T-Tg);
if (dProp!= NULL) {
- *dProp = -0.5*k*(Prop0-Prop1) * (1-pow(tanh(k*(Tc-Tg)),2)) - m;
+ *dProp = -0.5*k*(Prop0-Prop1) * (1-pow(tanh(k*(T-Tg)),2)) - m;
}
} else if (opt == 2) { // Huang et al
double Tg = 373; double k = 0.0994; double m = 2.968;
- Prop = Prop0 - (Prop0-Prop1) / (pow((m*exp(-m*k*(Tc-Tg))+1),(1/m)));
+ Prop = Prop0 - (Prop0-Prop1) / (pow((m*exp(-m*k*(T-Tg))+1),(1/m)));
if (dProp!= NULL) {
- *dProp = -(Prop0-Prop1)*(pow(m,2.0)*k*exp(-m*k*(Tc-Tg)))/( m* ( pow( (m*exp(-m*k*(Tc-Tg))+1),(1+1/m) ))); // Correct this derivative - CHANGE!!
+ *dProp = -(Prop0-Prop1)*(pow(m,2.0)*k*exp(-m*k*(T-Tg)))/( m* ( pow( (m*exp(-m*k*(T-Tg))+1),(1+1/m) ))); // Correct this derivative - CHANGE!!
}
}
}
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index 47241a4449d7010345a9dd9b1b4129a3e3d634dc..8622e8b2efb3cfe1c7238296aeb94a453564a617 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -342,13 +342,13 @@ void mlawNonLinearTVP::getFreeEnergyTVM(const IPNonLinearTVP *q0, IPNonLinearTVP
//
// psiVE
- double KT, GT, cteT, dKdT, ddKdTT, dGdT, ddGdTT, DcteDT, DDcteDTT; getK(KT,T,&dKdT,&ddKdTT); getG(GT,T,&dGdT,&ddGdTT); getAlpha(cteT,T,&DcteDT,&DDcteDTT);
+ double KT, GT, cteT, dKdT, ddKdTT, dGdT, ddGdTT, DcteDT, DDcteDTT; getK(q1,KT,T,&dKdT,&ddKdTT); getG(q1,GT,T,&dGdT,&ddGdTT); getAlpha(cteT,T,&DcteDT,&DDcteDTT);
const STensor3& E = q1->getConstRefToElasticStrain();
static STensor3 devKeinf, DdevKeinfDT, DDdevKeinfDTT, devEe;
static double trKeinf, DtrKeinfDT, DDtrKeinfDTT, trEe;
STensorOperation::decomposeDevTr(E,devEe,trEe);
- mlawNonLinearTVE::corKirInfinity(devEe,trEe,T,devKeinf,trKeinf);
+ mlawNonLinearTVE::corKirInfinity(q1,devEe,trEe,T,devKeinf,trKeinf);
DtrKeinfDT = trKeinf*dKdT/KT - 3*KT*(DcteDT*(T-_Tinitial) + cteT);
DDtrKeinfDTT = trKeinf*ddKdTT/KT - 6*dKdT*(DcteDT*(T-_Tinitial) + cteT)- 3*KT*(DDcteDTT*(T-_Tinitial) + 2*DcteDT);
@@ -632,8 +632,9 @@ 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 *= 0.8;
+ // mechSourceTVP -= STensorOperation::doubledot(tempX,alpha1);
+ // mechSourceTVP -= (R*Dgamma)/this->getTimeStep();
}
*Wm = mechSourceTVP;
}
@@ -686,11 +687,11 @@ 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) += 0.8 * ( 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();
}
}
}
@@ -700,11 +701,11 @@ void mlawNonLinearTVP::getMechSourceTVP(const STensor3& F0, const STensor3& F,
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 += 0.8*( 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();
}
}
*dWmdT = dmechSourceTVPdT;
diff --git a/NonLinearSolver/materialLaw/mullinsEffect.cpp b/NonLinearSolver/materialLaw/mullinsEffect.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..022b622a808ce228e9b79149ce1691940ad1238d
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mullinsEffect.cpp
@@ -0,0 +1,132 @@
+//
+// Description: Define mullins effect for softening in polymers in (pseudo-)hyperelastic regime
+//
+//
+// Author: <Ujwal Kishore J.>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#include "mullinsEffect.h"
+
+mullinsEffect::mullinsEffect(const int num,const bool init): _num(num), _initialized(init){}
+
+mullinsEffect::mullinsEffect(const mullinsEffect &source)
+{
+ _num = source._num;
+ _initialized = source._initialized;
+}
+
+mullinsEffect& mullinsEffect::operator=(const mullinsEffect &source)
+{
+ _num = source._num;
+ _initialized = source._initialized;
+ return *this;
+}
+
+
+negativeExponentialScaler::negativeExponentialScaler(const int num, double r, double m, bool init):
+ mullinsEffect(num,init),_r(r),_m(m){
+
+ _temFunc_r= new constantScalarFunction(1.);
+ _temFunc_m= new constantScalarFunction(1.);
+};
+
+negativeExponentialScaler::negativeExponentialScaler(const negativeExponentialScaler& src):
+ mullinsEffect(src), _r(src._r), _m(src._m){
+
+ _temFunc_r = NULL;
+ if (src._temFunc_r != NULL){
+ _temFunc_r = src._temFunc_r->clone();
+ }
+
+ _temFunc_m = NULL;
+ if (src._temFunc_m != NULL){
+ _temFunc_m = src._temFunc_m->clone();
+ }
+};
+
+negativeExponentialScaler& negativeExponentialScaler::operator =(const mullinsEffect& src){
+ mullinsEffect::operator=(src);
+ const negativeExponentialScaler* psrc = dynamic_cast<const negativeExponentialScaler*>(&src);
+ if(psrc != NULL)
+ {
+ _r = psrc->_r;
+ _m = psrc->_m;
+
+ if (psrc->_temFunc_r != NULL){
+ if (_temFunc_r != NULL){
+ _temFunc_r->operator=(*psrc->_temFunc_r);
+ }
+ else{
+ _temFunc_r = psrc->_temFunc_r->clone();
+ }
+ }
+
+ if (psrc->_temFunc_m != NULL){
+ if (_temFunc_m != NULL){
+ _temFunc_m->operator=(*psrc->_temFunc_m);
+ }
+ else{
+ _temFunc_m = psrc->_temFunc_m->clone();
+ }
+ }
+ }
+ return *this;
+};
+
+negativeExponentialScaler::~negativeExponentialScaler(){
+ if (_temFunc_r != NULL) {delete _temFunc_r; _temFunc_r = NULL;};
+ if (_temFunc_m != NULL) {delete _temFunc_m; _temFunc_m = NULL;};
+};
+
+void negativeExponentialScaler::setTemperatureFunction_r(const scalarFunction& Tfunc){
+ if (_temFunc_r != NULL) delete _temFunc_r;
+ _temFunc_r = Tfunc.clone();
+}
+
+void negativeExponentialScaler::createIPVariable(IPMullinsEffect* &ipv) const
+{
+ if(ipv != NULL) delete ipv;
+ ipv = new IPMullinsEffect();
+}
+
+void negativeExponentialScaler::mullinsEffectScaling(double psi0, double _psi, IPMullinsEffect &ipv) const{
+
+ double psi = _psi;
+
+ double eta = 0., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.;
+
+ if (psi<psi0){
+ eta = 1. - _r*( 1. - exp( - sqrt(_m*(psi0-psi)) ) );
+ }
+ else{
+ eta = 1.;
+ }
+
+ ipv.set(eta,DetaDpsi,DDetaDpsipsi,DetaDT,DDetaDTT);
+}
+
+void negativeExponentialScaler::mullinsEffectScaling(double psi0, double _psi, double T, IPMullinsEffect &ipv) const{
+
+ double psi = _psi;
+ double r = _r*_temFunc_r->getVal(T);
+ double m = _m*_temFunc_m->getVal(T);
+
+ double eta = 0., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.;
+
+ if (psi<psi0){
+ eta = 1. - r*( 1. - exp( - sqrt(m*(psi0-psi)) ) );
+ }
+ else{
+ eta = 1.;
+ }
+
+ ipv.set(eta,DetaDpsi,DDetaDpsipsi,DetaDT,DDetaDTT);
+}
+
+mullinsEffect *negativeExponentialScaler::clone() const
+{
+ return new negativeExponentialScaler(*this);
+}
+
diff --git a/NonLinearSolver/materialLaw/mullinsEffect.h b/NonLinearSolver/materialLaw/mullinsEffect.h
new file mode 100644
index 0000000000000000000000000000000000000000..bbca308bc5dda7e811531773d623c4fae7ddffe3
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mullinsEffect.h
@@ -0,0 +1,71 @@
+//
+// Description: Define mullins effect for softening in polymers in (pseudo-)hyperelastic regime
+//
+//
+// Author: <Ujwal Kishore J.>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#ifndef MULLINSEFFECT_H_
+#define MULLINSEFFECT_H_
+
+#include "ipMullinsEffect.h"
+#include "fullMatrix.h"
+#include "scalarFunction.h"
+
+class mullinsEffect {
+ public:
+ enum scalingFunction{negativeExponentialScaler};
+
+ protected:
+ int _num; // number of law (must be unique !)
+ bool _initialized; // to initialize law
+
+ public:
+ #ifndef SWIG
+ mullinsEffect(const int num, const bool init=true);
+ virtual ~mullinsEffect(){}
+ mullinsEffect(const mullinsEffect &source);
+ virtual mullinsEffect& operator=(const mullinsEffect &source);
+ virtual int getNum() const{return _num;}
+ virtual scalingFunction getType() const=0;
+ virtual void createIPVariable(IPMullinsEffect* &ipv) const=0;
+ virtual void initLaws(const std::map<int,mullinsEffect*> &maplaw)=0;
+ virtual const bool isInitialized() {return _initialized;}
+ virtual void mullinsEffectScaling(double psi0, double _psi, IPMullinsEffect &ipv) const=0;
+ virtual void mullinsEffectScaling(double psi0, double _psi, double T, IPMullinsEffect &ipv) const=0;
+ virtual mullinsEffect * clone() const=0;
+ #endif
+
+};
+
+// eta = 1 - r* ( 1 - exp(-sqrt(m*(psi0_max-psi0))) ) -> Ricker et al, 2021 (modified Zuniga, 2005)
+class negativeExponentialScaler : public mullinsEffect{
+ #ifndef SWIG
+ protected:
+ double _r;
+ double _m;
+ scalarFunction* _temFunc_r;
+ scalarFunction* _temFunc_m;
+ #endif // SWIG
+
+ public:
+ negativeExponentialScaler(const int num, double r, double m, bool init = true);
+ void setTemperatureFunction_r(const scalarFunction& Tfunc);
+ void setTemperatureFunction_m(const scalarFunction& Tfunc);
+ #ifndef SWIG
+ virtual ~negativeExponentialScaler();
+ negativeExponentialScaler(const negativeExponentialScaler& src);
+ virtual negativeExponentialScaler& operator =(const mullinsEffect& src);
+ virtual scalingFunction getType() const{return mullinsEffect::negativeExponentialScaler; };
+ virtual void createIPVariable(IPMullinsEffect* &ipv) const;
+ virtual void initLaws(const std::map<int,mullinsEffect*> &maplaw) {}; //nothing now, we will see if we use the mapping
+ virtual void mullinsEffectScaling(double psi0, double _psi, IPMullinsEffect &ipv) const;
+ virtual void mullinsEffectScaling(double psi0, double _psi, double T, IPMullinsEffect &ipv) const;
+ virtual mullinsEffect * clone() const;
+ #endif // SWIG
+};
+
+
+#endif // MULLINSEFFECT_H_
diff --git a/NonLinearSolver/nlsolver/scalarFunction.h b/NonLinearSolver/nlsolver/scalarFunction.h
index 94ee94541c0d4f5afee402805e59385f04056510..5e8b88af214b72647957642800e2525a7fb14812 100644
--- a/NonLinearSolver/nlsolver/scalarFunction.h
+++ b/NonLinearSolver/nlsolver/scalarFunction.h
@@ -574,6 +574,39 @@ class negExpWLFshiftFactor : public scalarFunction
#endif // SWIG
};
+// function y = exp(- C(x-Tref) )
+class negativeExponentialFunction : public scalarFunction
+{
+ protected:
+ double _C, _Tref;
+
+ public:
+ negativeExponentialFunction(double C, double Tref) : scalarFunction(),_C(C),_Tref(Tref){}
+ #ifndef SWIG
+ negativeExponentialFunction(const negativeExponentialFunction& src) : _C(src._C),_Tref(src._Tref){}
+ virtual negativeExponentialFunction& operator = (const scalarFunction& src){
+ scalarFunction::operator=(src);
+ const negativeExponentialFunction* psrc = dynamic_cast<const negativeExponentialFunction*>(&src);
+ if (psrc!= NULL){
+ _C = psrc->_C;
+ _Tref = psrc->_Tref;
+ }
+ return *this;
+ }
+ virtual ~negativeExponentialFunction(){};
+ virtual double getVal(const double x) const {
+ return exp(-_C*(x-_Tref));
+ };
+ virtual double getDiff(const double x) const {
+ return (-_C*exp(-_C*(x-_Tref)));
+ };
+ virtual double getDoubleDiff(const double x) const {
+ return (pow(_C,2)*exp(-_C*(x-_Tref)));
+ };
+ virtual scalarFunction* clone() const {return new negativeExponentialFunction(*this);};
+ #endif // SWIG
+};
+
#endif // SCALARFUNCTION_H_