Skip to content
Snippets Groups Projects
Commit fc427fe8 authored by Ujwal Kishore Jinaga's avatar Ujwal Kishore Jinaga :clown:
Browse files

mullins_framework

parent 356beaad
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
NonLinearSolver/internalPoints/ipTVEMullins.cpp
+ 46
− 0
View file @ fc427fe8
...@@ -17,9 +17,35 @@ IPTVEMullins::IPTVEMullins(const J2IsotropicHardening* comp, ...@@ -17,9 +17,35 @@ IPTVEMullins::IPTVEMullins(const J2IsotropicHardening* comp,
_T(298.15),_devDE(0.),_trDE(0.),_DE(0.),_DcorKirDT(0.),_DDcorKirDTT(0.), _T(298.15),_devDE(0.),_trDE(0.),_DE(0.),_DcorKirDT(0.),_DDcorKirDTT(0.),
_mechSrcTVE(0.),_DmechSrcTVEdT(0.),_DmechSrcTVEdF(0.), _mechSrcTVE(0.),_DmechSrcTVEdT(0.),_DmechSrcTVEdF(0.),
_thermalEnergy(0.),_dtShift(0.),_DdtShiftDT(0.),_DDdtShiftDDT(0.){ _thermalEnergy(0.),_dtShift(0.),_DdtShiftDT(0.),_DDdtShiftDDT(0.){
_devOGammai.clear();
_trOGammai.clear();
_devDOGammaiDT.clear();
_trDOGammaiDT.clear();
_devOi.clear();
_trOi.clear();
_devDOiDT.clear();
_trDOiDT.clear();
_devDDOiDTT.clear();
_trDDOiDTT.clear();
for (int i=0; i < N; i++){
STensor3 el(0.);
_devOGammai.push_back(el);
_trOGammai.push_back(0.);
_devDOGammaiDT.push_back(el);
_trDOGammaiDT.push_back(0.);
_devOi.push_back(el);
_trOi.push_back(0.);
_devDOiDT.push_back(el);
_trDOiDT.push_back(0.);
_devDDOiDTT.push_back(el);
_trDDOiDTT.push_back(0.);
}; };
IPTVEMullins::IPTVEMullins(const IPTVEMullins& src): IPHyperViscoElastoPlastic(src), IPTVEMullins::IPTVEMullins(const IPTVEMullins& src): IPHyperViscoElastoPlastic(src),
_devOGammai(src._devOGammai), _trOGammai(src._trOGammai),
_devDOGammaiDT(src._devDOGammaiDT), _trDOGammaiDT(src._trDOGammaiDT),
_devOi(src._devOi), _trOi(src._trOi), _devDOiDT(src._devDOiDT), _trDOiDT(src._trDOiDT), _devDDOiDTT(src._devDDOiDTT), _trDDOiDTT(src._trDDOiDTT),
_T(src._T),_devDE(src._devDE), _trDE(src._trDE), _DE(src._DE), _DcorKirDT(src._DcorKirDT), _DDcorKirDTT(src._DDcorKirDTT), _T(src._T),_devDE(src._devDE), _trDE(src._trDE), _DE(src._DE), _DcorKirDT(src._DcorKirDT), _DDcorKirDTT(src._DDcorKirDTT),
_mechSrcTVE(src._mechSrcTVE),_DmechSrcTVEdT(src._DmechSrcTVEdT),_DmechSrcTVEdF(src._DmechSrcTVEdF), _mechSrcTVE(src._mechSrcTVE),_DmechSrcTVEdT(src._DmechSrcTVEdT),_DmechSrcTVEdF(src._DmechSrcTVEdF),
_thermalEnergy(src._thermalEnergy), _thermalEnergy(src._thermalEnergy),
...@@ -44,6 +70,16 @@ IPTVEMullins& IPTVEMullins::operator=(const IPVariable& src) ...@@ -44,6 +70,16 @@ IPTVEMullins& IPTVEMullins::operator=(const IPVariable& src)
_dtShift = psrc->_dtShift; _dtShift = psrc->_dtShift;
_DdtShiftDT = psrc->_DdtShiftDT; _DdtShiftDT = psrc->_DdtShiftDT;
_DDdtShiftDDT = psrc->_DDdtShiftDDT; _DDdtShiftDDT = psrc->_DDdtShiftDDT;
_devOGammai = psrc->_devOGammai;
_trOGammai = psrc->_trOGammai;
_devDOGammaiDT = psrc->_devDOGammaiDT;
_trDOGammaiDT = psrc->_trDOGammaiDT;
_devOi = psrc->_devOi;
_trOi = psrc->_trOi;
_devDOiDT = psrc->_devDOiDT;
_trDOiDT = psrc->_trDOiDT;
_devDDOiDTT = psrc->_devDDOiDTT;
_trDDOiDTT = psrc->_trDDOiDTT;
} }
return *this; return *this;
} }
...@@ -68,4 +104,14 @@ void IPTVEMullins::restart(){ ...@@ -68,4 +104,14 @@ void IPTVEMullins::restart(){
restartManager::restart(_dtShift); restartManager::restart(_dtShift);
restartManager::restart(_DdtShiftDT); restartManager::restart(_DdtShiftDT);
restartManager::restart(_DDdtShiftDDT); restartManager::restart(_DDdtShiftDDT);
restartManager::restart(_devOGammai);
restartManager::restart(_trOGammai);
restartManager::restart(_devDOGammaiDT);
restartManager::restart(_trDOGammaiDT);
restartManager::restart(_devOi);
restartManager::restart(_trOi);
restartManager::restart(_devDOiDT);
restartManager::restart(_trDOiDT);
restartManager::restart(_devDDOiDTT);
restartManager::restart(_trDDOiDTT);
} }
NonLinearSolver/internalPoints/ipTVEMullins.h
+ 44
− 0
View file @ fc427fe8
...@@ -30,6 +30,20 @@ class IPTVEMullins : public IPHyperViscoElastoPlastic{ ...@@ -30,6 +30,20 @@ class IPTVEMullins : public IPHyperViscoElastoPlastic{
double _mechSrcTVE; double _mechSrcTVE;
double _DmechSrcTVEdT; double _DmechSrcTVEdT;
STensor3 _DmechSrcTVEdF; STensor3 _DmechSrcTVEdF;
// Viscoelastic Strain
std::vector<STensor3> _devOGammai; // dev viscoelastic strain for each branch
std::vector<double> _trOGammai; // tr viscoelastic strain for each branch
std::vector<STensor3> _devDOGammaiDT;
std::vector<double> _trDOGammaiDT;
// Viscoelastic Stress
std::vector<STensor3> _devOi; // dev viscoelastic stress for each branch
std::vector<double> _trOi; // tr viscoelastic stress for each branch
std::vector<STensor3> _devDOiDT;
std::vector<double> _trDOiDT;
std::vector<STensor3> _devDDOiDTT;
std::vector<double> _trDDOiDTT;
// DcorKirDT // DcorKirDT
STensor3 _DcorKirDT; STensor3 _DcorKirDT;
...@@ -50,6 +64,36 @@ class IPTVEMullins : public IPHyperViscoElastoPlastic{ ...@@ -50,6 +64,36 @@ class IPTVEMullins : public IPHyperViscoElastoPlastic{
virtual IPTVEMullins& operator=(const IPVariable& src); virtual IPTVEMullins& operator=(const IPVariable& src);
virtual ~IPTVEMullins(){}; virtual ~IPTVEMullins(){};
virtual std::vector<STensor3>& getRefToDevViscoElasticOverStrain() {return _devOGammai;};
virtual const std::vector<STensor3>& getConstRefToDevViscoElasticOverStrain() const {return _devOGammai;};
virtual std::vector<double>& getRefToTrViscoElasticOverStrain() {return _trOGammai;};
virtual const std::vector<double>& getConstRefToTrViscoElasticOverStrain() const {return _trOGammai;};
virtual std::vector<STensor3>& getRefToDevDOGammaiDT() {return _devDOGammaiDT;};
virtual const std::vector<STensor3>& getConstRefToDevDOGammaiDT() const {return _devDOGammaiDT;};
virtual std::vector<double>& getRefToTrDOGammaiDT() {return _trDOGammaiDT;};
virtual const std::vector<double>& getConstRefToTrDOGammaiDT() const {return _trDOGammaiDT;};
virtual std::vector<STensor3>& getRefToDevViscoElasticOverStress() {return _devOi;};
virtual const std::vector<STensor3>& getConstRefToDevViscoElasticOverStress() const {return _devOi;};
virtual std::vector<double>& getRefToTrViscoElasticOverStress() {return _trOi;};
virtual const std::vector<double>& getConstRefToTrViscoElasticOverStress() const {return _trOi;};
virtual std::vector<STensor3>& getRefToDevDOiDT() {return _devDOiDT;};
virtual const std::vector<STensor3>& getConstRefToDevDOiDT() const {return _devDOiDT;};
virtual std::vector<double>& getRefToTrDOiDT() {return _trDOiDT;};
virtual const std::vector<double>& getConstRefToTrDOiDT() const {return _trDOiDT;};
virtual std::vector<STensor3>& getRefToDevDDOiDTT() {return _devDDOiDTT;};
virtual const std::vector<STensor3>& getConstRefToDevDDOiDTT() const {return _devDDOiDTT;};
virtual std::vector<double>& getRefToTrDDOiDTT() {return _trDDOiDTT;};
virtual const std::vector<double>& getConstRefToTrDDOiDTT() const {return _trDDOiDTT;};
virtual STensor3& getRefToDcorKirDT() {return _DcorKirDT;}; virtual STensor3& getRefToDcorKirDT() {return _DcorKirDT;};
virtual const STensor3& getConstRefToDcorKirDT() const {return _DcorKirDT;}; virtual const STensor3& getConstRefToDcorKirDT() const {return _DcorKirDT;};
......
NonLinearSolver/materialLaw/mlawTVEMullins.cpp
+ 175
− 10
View file @ fc427fe8
...@@ -121,11 +121,15 @@ void mlawTVEMullins::setViscoElasticNumberOfElement(const int N){ ...@@ -121,11 +121,15 @@ void mlawTVEMullins::setViscoElasticNumberOfElement(const int N){
Msg::Info("Numer of Spring/Dashpot for viscoelastic model: %d",_N); Msg::Info("Numer of Spring/Dashpot for viscoelastic model: %d",_N);
_Ki.clear(); _ki.clear(); _Ki.clear(); _ki.clear();
_Gi.clear(); _gi.clear(); _Gi.clear(); _gi.clear();
_nKi.clear(); _nGi.clear();
_Ki.resize(_N,0.); _Ki.resize(_N,0.);
_ki.resize(_N,0.); _ki.resize(_N,0.);
_Gi.resize(_N,0.); _Gi.resize(_N,0.);
_gi.resize(_N,0.); _gi.resize(_N,0.);
_nKi.resize(_N,0.);
_nGi.resize(_N,0.);
}; };
void mlawTVEMullins::setViscoElasticData(const int i, const double Ei, const double taui){ void mlawTVEMullins::setViscoElasticData(const int i, const double Ei, const double taui){
...@@ -139,6 +143,9 @@ void mlawTVEMullins::setViscoElasticData(const int i, const double Ei, const dou ...@@ -139,6 +143,9 @@ void mlawTVEMullins::setViscoElasticData(const int i, const double Ei, const dou
_ki[i-1] = taui; // /(3.*(1.-2.*_nu)); _ki[i-1] = taui; // /(3.*(1.-2.*_nu));
_Gi[i-1] = GG; _Gi[i-1] = GG;
_gi[i-1] = taui; // /(2.*(1.+_nu)); _gi[i-1] = taui; // /(2.*(1.+_nu));
_nKi[i-1] = KK/_K;
_nGi[i-1] = GG/_G;
Msg::Info("setting: element=%d Ki= %e ki= %e, Gi= %e, gi= %e",i-1,KK,taui,GG,taui); Msg::Info("setting: element=%d Ki= %e ki= %e, Gi= %e, gi= %e",i-1,KK,taui,GG,taui);
} }
...@@ -150,6 +157,7 @@ void mlawTVEMullins::setViscoElasticData_Bulk(const int i, const double Ki, cons ...@@ -150,6 +157,7 @@ void mlawTVEMullins::setViscoElasticData_Bulk(const int i, const double Ki, cons
else{ else{
_Ki[i-1] = Ki; _Ki[i-1] = Ki;
_ki[i-1] = ki; _ki[i-1] = ki;
_nKi[i-1] = KK/_K;
// Msg::Info("setting: element=%d Ki = %e ki = %e, ",i-1,Ki,ki); // Msg::Info("setting: element=%d Ki = %e ki = %e, ",i-1,Ki,ki);
} }
}; };
...@@ -160,6 +168,7 @@ void mlawTVEMullins::setViscoElasticData_Shear(const int i, const double Gi, con ...@@ -160,6 +168,7 @@ void mlawTVEMullins::setViscoElasticData_Shear(const int i, const double Gi, con
else{ else{
_Gi[i-1] = Gi; _Gi[i-1] = Gi;
_gi[i-1] = gi; _gi[i-1] = gi;
_nGi[i-1] = GG/_G;
// Msg::Info("setting: element=%d Gi = %e gi = %e, ",i-1,Gi,gi); // Msg::Info("setting: element=%d Gi = %e gi = %e, ",i-1,Gi,gi);
} }
}; };
...@@ -169,6 +178,7 @@ void mlawTVEMullins::setViscoElasticData(const std::string filename){ ...@@ -169,6 +178,7 @@ void mlawTVEMullins::setViscoElasticData(const std::string filename){
if (file == NULL) Msg::Error("file: %s does not exist",filename.c_str()); if (file == NULL) Msg::Error("file: %s does not exist",filename.c_str());
_Ki.clear(); _ki.clear(); _Ki.clear(); _ki.clear();
_Gi.clear(); _gi.clear(); _Gi.clear(); _gi.clear();
_nKi.clear(); _nGi.clear();
_N = 0; _N = 0;
Msg::Info("reading viscoelastic input data"); Msg::Info("reading viscoelastic input data");
...@@ -176,13 +186,18 @@ void mlawTVEMullins::setViscoElasticData(const std::string filename){ ...@@ -176,13 +186,18 @@ void mlawTVEMullins::setViscoElasticData(const std::string filename){
int ok = fscanf(file,"%d",&_N); int ok = fscanf(file,"%d",&_N);
Msg::Info("Numer of Maxwell elements: %d",_N); Msg::Info("Numer of Maxwell elements: %d",_N);
double KK, k, GG, g; double KK, k, GG, g;
double tempK, tempG = 0.;
for (int i=0; i<_N; i++){ for (int i=0; i<_N; i++){
ok = fscanf(file,"%lf %lf %lf %lf",&KK,&k,&GG,&g); ok = fscanf(file,"%lf %lf %lf %lf",&KK,&k,&GG,&g);
tempK = KK/_K;
tempG = GG/_G;
_Ki.push_back(KK); _Ki.push_back(KK);
_ki.push_back(k); _ki.push_back(k);
_Gi.push_back(GG); _Gi.push_back(GG);
_gi.push_back(g); _gi.push_back(g);
_nKi.push_back(tempK);
_nGi.push_back(tempG);
Msg::Info("Maxwell element %d: K[%d] = %e, k[%d] = %e, G[%d] = %e, g[%d] = %e", Msg::Info("Maxwell element %d: K[%d] = %e, k[%d] = %e, G[%d] = %e, g[%d] = %e",
i,i,KK,i,k,i,GG,i,g); i,i,KK,i,k,i,GG,i,g);
} }
...@@ -203,10 +218,157 @@ double mlawTVEMullins::ShiftFactorWLF(const double T, double *DaDT, double *DDaD ...@@ -203,10 +218,157 @@ double mlawTVEMullins::ShiftFactorWLF(const double T, double *DaDT, double *DDaD
return a_WLF; return a_WLF;
} }
double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double Tc) const{ void mlawTVEMullins::shiftedTime(const double t, const double T0, const double T1, double *t_shift, double *Ddt_shiftDT, double *DDdt_shiftDT,
const int opt_time, const int opt_order) const{
// Shift t -> t is either delta_t or t (final time). This function returns either t^* or delta_t^*.
// opt_time -> 0 or 1 (default), depending on full-step or mid-step integration approximation for shifted time step
// opt_order -> 1 or 2 (default) or 4, order of Guassian quadrature (1 doesnt work, 4 probably works)
double dt = this->getTimeStep();
double t1 = this->getTime();
double t0 = t1 - dt;
// Linearly Interpolate Temperture -> This means dT/dt is constant between tn and tn+1
double ts, t_x0, t_x1, t_x2 = 0.;
double t_x6, t_x7, t_x8, t_x9 = 0.;
if (opt_time == 0){
ts = 2;
t_x0 = (t1+t0)/2;
t_x1 = (t1+t0)/2 - dt/2 * (1/sqrt(3));
t_x2 = (t1+t0)/2 + dt/2 * (1/sqrt(3));
t_x6 = (t1+t0)/2 - dt/2 * sqrt(3/7 - 2/7*sqrt(6/5));
t_x7 = (t1+t0)/2 + dt/2 * sqrt(3/7 - 2/7*sqrt(6/5));
t_x8 = (t1+t0)/2 - dt/2 * sqrt(3/7 + 2/7*sqrt(6/5));
t_x9 = (t1+t0)/2 + dt/2 * sqrt(3/7 + 2/7*sqrt(6/5));
} else {
ts = 4;
t_x0 = (3*t1+t0)/4;
t_x1 = (3*t1+t0)/4 - dt/4 * (1/sqrt(3));
t_x2 = (3*t1+t0)/4 + dt/4 * (1/sqrt(3));
t_x6 = (3*t1+t0)/4 - dt/4 * sqrt(3/7 - 2/7*sqrt(6/5));
t_x7 = (3*t1+t0)/4 + dt/4 * sqrt(3/7 - 2/7*sqrt(6/5));
t_x8 = (3*t1+t0)/4 - dt/4 * sqrt(3/7 + 2/7*sqrt(6/5));
t_x9 = (3*t1+t0)/4 + dt/4 * sqrt(3/7 + 2/7*sqrt(6/5));
}
double T_x0 = linearInterpTemp(t_x0, t0, t1, T0, T1);
double T_x1 = linearInterpTemp(t_x1, t0, t1, T0, T1);
double T_x2 = linearInterpTemp(t_x2, t0, t1, T0, T1);
double T_x6 = linearInterpTemp(t_x6, t0, t1, T0, T1);
double T_x7 = linearInterpTemp(t_x7, t0, t1, T0, T1);
double T_x8 = linearInterpTemp(t_x8, t0, t1, T0, T1);
double T_x9 = linearInterpTemp(t_x9, t0, t1, T0, T1);
// double T_x1 = 1/dt*( (t_1 - t_x1) )*T_0 + 1/dt*( (t_x1 - t_0) )*T_1 ;
// double T_x2 = 1/dt*( (t_1 - t_x2) )*T_0 + 1/dt*( (t_x2 - t_0) )*T_1 ;
if (t_shift!=NULL){
// double a_T1 = ShiftFactorWLF(T1);
if (opt_order == 1){
// 1 point Gauss rule
double a_x0 = ShiftFactorWLF(T_x0);
*t_shift = t/ts * (1/a_x0);
} else if (opt_order == 2){
// 2 point Gauss rule
double a_x1 = ShiftFactorWLF(T_x1);
double a_x2 = ShiftFactorWLF(T_x2);
*t_shift = t/ts * (1/a_x1 + 1/a_x2);
} else {
// 4 point Gauss rule
double a_x6 = ShiftFactorWLF(T_x6);
double a_x7 = ShiftFactorWLF(T_x7);
double a_x8 = ShiftFactorWLF(T_x8);
double a_x9 = ShiftFactorWLF(T_x9);
*t_shift = t/ts * (((18+sqrt(30))/36)*(1/a_x6 + 1/a_x7) + ((18-sqrt(30))/36)*(1/a_x8 + 1/a_x9));
}
}
if (Ddt_shiftDT!=NULL){
// double DaDT_T1 = 0.; double a_T1 = ShiftFactorWLF(T1,&DaDT_T1);
if (dt > 0.){
if (opt_order == 1){
// 1 point Gauss rule
double DaDT_x0 = 0.;
double a_x0 = ShiftFactorWLF(T_x0,&DaDT_x0);
*Ddt_shiftDT = t/ts * (DaDT_x0 * (t_x0-t0)/dt );
} else if (opt_order == 2){
// 2 point Gauss rule
double DaDT_x1, DaDT_x2 = 0.;
double a_x1 = ShiftFactorWLF(T_x1,&DaDT_x1);
double a_x2 = ShiftFactorWLF(T_x2,&DaDT_x2);
*Ddt_shiftDT = t/ts * (DaDT_x1 * (t_x1-t0)/dt + DaDT_x2 * (t_x2-t0)/dt);
} else {
// 4 point Gauss rule
double DaDT_x6, DaDT_x7, DaDT_x8, DaDT_x9 = 0.;
double a_x6 = ShiftFactorWLF(T_x6,&DaDT_x6);
double a_x7 = ShiftFactorWLF(T_x7,&DaDT_x7);
double a_x8 = ShiftFactorWLF(T_x8,&DaDT_x8);
double a_x9 = ShiftFactorWLF(T_x9,&DaDT_x9);
*Ddt_shiftDT = t/ts * ( ((18+sqrt(30))/36)*(DaDT_x6 * (t_x6-t0)/dt + DaDT_x7 * (t_x7-t0)/dt ) +
((18-sqrt(30))/36)*(DaDT_x8 * (t_x8-t0)/dt + DaDT_x9 * (t_x9-t0)/dt ) );
}
} else { *Ddt_shiftDT = 0.;}
}
if (DDdt_shiftDT!=NULL){
// double DDaDT_T1 = 0.; double a_T1 = ShiftFactorWLF(T1,NULL,&DDaDT_T1);
if (dt > 0.){
if (opt_order == 1){
// 1 point Gauss rule
double DDaDT_x0 = 0.;
double a_x0 = ShiftFactorWLF(T_x0,NULL,&DDaDT_x0);
*DDdt_shiftDT = t/ts * ( DDaDT_x0 * pow( ((t_x0-t0)/dt) ,2 ) );
} else if (opt_order == 2){
// 2 point Gauss rule
double DDaDT_x1, DDaDT_x2 = 0.;
double a_x1 = ShiftFactorWLF(T_x1,NULL,&DDaDT_x1);
double a_x2 = ShiftFactorWLF(T_x2,NULL,&DDaDT_x2);
*DDdt_shiftDT = t/ts * ( DDaDT_x1 * pow( ((t_x1-t0)/dt) ,2 ) + DDaDT_x2 * pow( ((t_x2-t0)/dt) ,2 ) );
} else {
// 4 point Gauss rule
double DDaDT_x6, DDaDT_x7, DDaDT_x8, DDaDT_x9 = 0.;
double a_x6 = ShiftFactorWLF(T_x6,NULL,&DDaDT_x6);
double a_x7 = ShiftFactorWLF(T_x7,NULL,&DDaDT_x7);
double a_x8 = ShiftFactorWLF(T_x8,NULL,&DDaDT_x8);
double a_x9 = ShiftFactorWLF(T_x9,NULL,&DDaDT_x9);
*Ddt_shiftDT = t/ts * ( ((18+sqrt(30))/36)*( DDaDT_x6 * pow( (t_x6-t0)/dt, 2 ) + DDaDT_x7 * pow( (t_x7-t0)/dt, 2 )) +
((18-sqrt(30))/36)*( DDaDT_x8 * pow( (t_x8-t0)/dt, 2 ) + DDaDT_x9 * pow( (t_x9-t0)/dt, 2 )) );
}
} else { *DDdt_shiftDT = 0.;}
}
}
double mlawTVEMullins::linearInterpTemp(const double t_interp, const double t0, const double t1, const double T0, const double T1) const{
// Update the Properties to the current temperature (see below which ones are being used) // This function linearly interpolates temperature for the given intermediate time -> t_interp is tn + fraction of step (i.e., midstep = tn + 0.5*dt)
double KT, GT, AlphaT; getK(KT,Tc); getG(GT,Tc); getAlpha(AlphaT,Tc); // Also, linear interp -> constant gradients
double dt = t1 - t0;
double T_interp;
if (dt > 0.){
T_interp = 1/dt*( (t1 - t_interp) )*T0 + 1/dt*( (t_interp - t0) )*T1;
}
else { T_interp = T0;}
return T_interp;
}
double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double T) const{
// Update the Properties to the current temperature
double KT, GT, AlphaT; getK(KT,T); getG(GT,T); getAlpha(AlphaT,T);
double Psy_mech = 0.; double Psy_mech = 0.;
...@@ -216,7 +378,7 @@ double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double ...@@ -216,7 +378,7 @@ double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double
STensorOperation::decomposeDevTr(q._Ee,devEe,trEe); STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
// Get thermal strain // Get thermal strain
double Eth = 3.*AlphaT*(Tc-_Tinitial); double Eth = 3.*AlphaT*(T-_Tinitial);
// Get Equilibrium freeEnergy_mech // Get Equilibrium freeEnergy_mech
Psy_mech = KT*0.5*(trEe-Eth)*(trEe-Eth) + GT*STensorOperation::doubledot(devEe,devEe); Psy_mech = KT*0.5*(trEe-Eth)*(trEe-Eth) + GT*STensorOperation::doubledot(devEe,devEe);
...@@ -247,7 +409,10 @@ double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double ...@@ -247,7 +409,10 @@ double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double
} }
void mlawTVEMullins::corKirInfinity(const STensor3& devEe, const double& trEe, const double T0, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const{ void mlawTVEMullins::corKirInfinity(const STensor3& devEe, const double& trEe, const double T0, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const{
double KT, GT, AlphaT;
getK(KT,T); getG(GT,T); getAlpha(AlphaT,T);
CorKirDevInf = 2*GT*devEe;
CorKirTrInf = KT*(trEe - 3*AlphaT*(T-_Tinitial));
} }
void mlawTVEMullins::getTVEdCorKirDT(const IPTVEMullins *q0, IPTVEMullins *q1, const double& T0, const double& T) const{ void mlawTVEMullins::getTVEdCorKirDT(const IPTVEMullins *q0, IPTVEMullins *q1, const double& T0, const double& T) const{
...@@ -263,7 +428,7 @@ void mlawTVEMullins::getMechSourceTVE(const IPTVEMullins *q0, IPTVEMullins *q1, ...@@ -263,7 +428,7 @@ void mlawTVEMullins::getMechSourceTVE(const IPTVEMullins *q0, IPTVEMullins *q1,
void mlawTVEMullins::ThermoViscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0, void mlawTVEMullins::ThermoViscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
const IPTVEMullins *q0, IPTVEMullins *q1, const IPTVEMullins *q0, IPTVEMullins *q1,
double& Ke, double& Ge, double& Ke, double& Ge,
const double T0, const double T1, const bool calculateStress) const{ const double T0, const double T, const bool calculateStress) const{
// Get LogStrain // Get LogStrain
...@@ -277,7 +442,7 @@ void mlawTVEMullins::ThermoViscoElasticPredictor(const STensor3& Ee, const STens ...@@ -277,7 +442,7 @@ void mlawTVEMullins::ThermoViscoElasticPredictor(const STensor3& Ee, const STens
// Get corKirInfinity // Get corKirInfinity
static double CorKirTrInf; static double CorKirTrInf;
static STensor3 CorKirDevInf; static STensor3 CorKirDevInf;
corKirInfinity(devEe,trEe,T0,T1,CorKirDevInf,CorKirTrInf); corKirInfinity(devEe,trEe,T0,T,CorKirDevInf,CorKirTrInf);
// Initialise CorKirStress - Dev // Initialise CorKirStress - Dev
static STensor3 devK; static STensor3 devK;
...@@ -295,9 +460,9 @@ void mlawTVEMullins::ThermoViscoElasticPredictor(const STensor3& Ee, const STens ...@@ -295,9 +460,9 @@ void mlawTVEMullins::ThermoViscoElasticPredictor(const STensor3& Ee, const STens
double dt_shift_1, dt_shift_2 = 0.; double dt_shift_1, dt_shift_2 = 0.;
double Ddt_shiftDT, DDdt_shiftDT = 0.; double Ddt_shiftDT, DDdt_shiftDT = 0.;
// dt_shift = dt; // dt_shift = dt;
shiftedTime(dt, T0, T1, &dt_shift_1,NULL,NULL,0); shiftedTime(dt, T0, T, &dt_shift_1,NULL,NULL,0);
shiftedTime(dt, T0, T1, &dt_shift_2,&Ddt_shiftDT,&DDdt_shiftDT); shiftedTime(dt, T0, T, &dt_shift_2,&Ddt_shiftDT,&DDdt_shiftDT);
// double aT_1 = ShiftFactorWLF(T1); // double aT_1 = ShiftFactorWLF(T);
std::vector<double> KiT, GiT, DKiDT, DGiDT, DDKiDT, DDGiDT; std::vector<double> KiT, GiT, DKiDT, DGiDT, DDKiDT, DDGiDT;
KiT.resize(_N,0.); GiT.resize(_N,0.); KiT.resize(_N,0.); GiT.resize(_N,0.);
......
NonLinearSolver/materialLaw/mlawTVEMullins.h
+ 7
− 3
View file @ fc427fe8
...@@ -49,8 +49,8 @@ class mlawTVEMullins : public mlawPowerYieldHyper{ ...@@ -49,8 +49,8 @@ class mlawTVEMullins : public mlawPowerYieldHyper{
scalarFunction* _temFunc_Cp; // specific heat temperature function scalarFunction* _temFunc_Cp; // specific heat temperature function
// Scalars - Additional constants // Scalars - Additional constants
std::vector<double> nKi; // branch free energy scalars std::vector<double> _nKi; // normalized branch Bulk Moduli
std::vector<double> nGi; // branch free energy scalars std::vector<double> _nGi; // normalized branch Shear Moduli
// Shift factor - Everything related to it // Shift factor - Everything related to it
double _C1; // WLF Universal Constant 1 - set in .py file double _C1; // WLF Universal Constant 1 - set in .py file
...@@ -61,7 +61,11 @@ class mlawTVEMullins : public mlawPowerYieldHyper{ ...@@ -61,7 +61,11 @@ class mlawTVEMullins : public mlawPowerYieldHyper{
protected: protected:
// Const Functions for constitutive function - To update any variable, pass it as an argument by reference. // Const Functions for constitutive function - To update any variable, pass it as an argument by reference.
virtual double ShiftFactorWLF(const double T, double *DaDT = NULL, double *DDaDDT = NULL) const ; virtual double ShiftFactorWLF(const double T, double *DaDT = NULL, double *DDaDDT = NULL) const ;
virtual double freeEnergyMechanical(const IPTVEMullins& q, const double Tc) const ; virtual void shiftedTime(const double t, const double T0, const double T1,
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 IPTVEMullins& q, const double T) const ;
void corKirInfinity(const STensor3& devEe, const double& trEe, const double T0, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const; void corKirInfinity(const STensor3& devEe, const double& trEe, const double T0, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment