Select Git revision
mlawNonLocalPorousCoupled.cpp
-
Julien Leclerc authoredJulien Leclerc authored
mlawNonLinearTVP.cpp 96.91 KiB
//
// C++ Interface: Material Law
//
// Description: Non-Linear Thermo-Visco-Mechanics (Thermo-ViscoElasto-ViscoPlasto Law) with Non-Local Damage Interface (soon.....)
//
// Author: <Ujwal Kishore J - FLE_Knight>, (C) 2022
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "STensorOperations.h"
#include "nonLinearMechSolver.h"
#include "mlawNonLinearTVP.h"
mlawNonLinearTVP::mlawNonLinearTVP(const int num,const double E,const double nu, const double rho, const double tol,
const double Tinitial, const double Alpha, const double KThCon, const double Cp,
const bool matrixbyPerturbation, const double pert, const bool thermalEstimationPreviousConfig):
mlawNonLinearTVE(num, E, nu, rho, tol, Tinitial, Alpha, KThCon, Cp, matrixbyPerturbation, pert, thermalEstimationPreviousConfig),
_TaylorQuineyFactor(0.9), _HR(1.){
// by default, no temperature dependence
_temFunc_Sy0 = new constantScalarFunction(1.);
_temFunc_H = new constantScalarFunction(1.);
_temFunc_Hb = new constantScalarFunction(1.);
_temFunc_visc = new constantScalarFunction(1.);
};
// , mlawPowerYieldHyper(src)
mlawNonLinearTVP::mlawNonLinearTVP(const mlawNonLinearTVP& src): mlawNonLinearTVE(src) {
_TaylorQuineyFactor = src._TaylorQuineyFactor;
_HR = src._HR;
_temFunc_Sy0 = NULL; // temperature dependence of initial yield stress
if (src._temFunc_Sy0!= NULL)
_temFunc_Sy0 = src._temFunc_Sy0->clone();
_temFunc_H = NULL; // temperature dependence of hardening stress
if (src._temFunc_H != NULL)
_temFunc_H = src._temFunc_H->clone();
_temFunc_Hb = NULL; // temperature dependence of kinematic hardening
if (src._temFunc_Hb != NULL)
_temFunc_Hb = src._temFunc_Hb->clone();
_temFunc_visc = NULL; // temperature dependence of viscosity
if (src._temFunc_visc != NULL)
_temFunc_visc = src._temFunc_visc->clone();
};
mlawNonLinearTVP& mlawNonLinearTVP::operator=(const materialLaw& source){
mlawNonLinearTVE::operator=(source);
// mlawPowerYieldHyper::operator=(source);
const mlawNonLinearTVP* src =dynamic_cast<const mlawNonLinearTVP*>(&source);
if(src != NULL){
_TaylorQuineyFactor = src->_TaylorQuineyFactor;
_HR = src->_HR;
if(_temFunc_Sy0 != NULL) delete _temFunc_Sy0; // temperature dependence of initial yield stress
if (src->_temFunc_Sy0!= NULL)
_temFunc_Sy0 = src->_temFunc_Sy0->clone();
if(_temFunc_H != NULL) delete _temFunc_H; // temperature dependence of hardening stress
if (src->_temFunc_H != NULL)
_temFunc_H = src->_temFunc_H->clone();
if(_temFunc_Hb != NULL) delete _temFunc_Hb; // temperature dependence of kinematic hardening
if (src->_temFunc_Hb != NULL) _temFunc_Hb = src->_temFunc_Hb->clone();
if(_temFunc_visc != NULL) delete _temFunc_visc; // temperature dependence of viscosity
if (src->_temFunc_visc != NULL)
_temFunc_visc = src->_temFunc_visc->clone();
}
return *this;
};
mlawNonLinearTVP::~mlawNonLinearTVP(){
if (_temFunc_Sy0 != NULL) delete _temFunc_Sy0; _temFunc_Sy0 = NULL;
if (_temFunc_H != NULL) delete _temFunc_H; _temFunc_H= NULL;
if (_temFunc_Hb != NULL) delete _temFunc_Hb; _temFunc_Hb= NULL;
if (_temFunc_visc != NULL) delete _temFunc_visc; _temFunc_visc= NULL;
};
// CHECK IF IPSTATE NEEDS the same definition as in mlawPowerYieldHyper --- WHAT??
// NEW
void mlawNonLinearTVP::setIsotropicHardeningCoefficients(const double HR1, const double HR2, const double HR3){
_HR.clear();
_HR.resize(3,0.);
_HR[0] = HR1;
_HR[1] = HR2;
_HR[2] = HR3;
};
void mlawNonLinearTVP::hardening(const IPNonLinearTVP* q0, IPNonLinearTVP* q1, const double& T) const{
//Msg::Error("epspCompression = %e, epspTRaction = %e, epspShear = %e",q->_epspCompression,q->_epspTraction,q->_epspShear);
if (_compression != NULL && q1->_ipCompression != NULL){
_compression->hardening(q0->_epspCompression, *q0->_ipCompression, q1->_epspCompression,T,*q1->_ipCompression);
}
if (_traction!= NULL && q1->_ipTraction != NULL){
_traction->hardening(q0->_epspTraction,*q0->_ipTraction, q1->_epspTraction,T,*q1->_ipTraction);
}
if (_kinematic!= NULL && q1->_ipKinematic != NULL)
_kinematic->hardening(q0->_epspbarre,*q0->_ipKinematic,q1->_epspbarre,T,*q1->_ipKinematic);
};
void mlawNonLinearTVP::getModifiedMandel(const STensor3& C, const IPNonLinearTVP *q0, IPNonLinearTVP *q1) const{
STensor3& M = q1->getRefToModifiedMandel();
const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
STensor3 Cinv;
STensorOperation::inverseSTensor3(C, Cinv);
STensor3 temp1, temp2;
STensorOperation::multSTensor3(C, corKir, temp1);
STensorOperation::multSTensor3(temp1, Cinv, temp2);
M = 0.5*(corKir + temp2);
};
void mlawNonLinearTVP::checkCommutavity(STensor3& commuteCheck,const STensor3& Ce, const STensor3& S, const IPNonLinearTVP *q1) const{
const STensor3& M = q1->getConstRefToModifiedMandel();
const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
STensor3 temp1, temp2, temp3, temp4, temp6;
STensorOperation::multSTensor3(Ce, S, temp1);
STensorOperation::multSTensor3(S, Ce, temp2);
STensorOperation::multSTensor3(corKir, Ce, temp3);
STensorOperation::multSTensor3(Ce, corKir, temp4);
// Check S commuteCheck = temp1;
commuteCheck -= temp2;
// Check corKir (just for fun)
temp6 = temp3;
temp6 -= temp4;
// RANK OF A ZERO TENSOR is ZERO - maybe use this?
};
void mlawNonLinearTVP::getChabocheCoeffs(fullVector<double>& coeffs, const double& opt, const IPNonLinearTVP *q1) const{
double p = q1->_epspbarre; // eq.plastic strain
coeffs.resize(2);
coeffs(0) = 1; //DEBUGGING
coeffs(1) = 0;
// add dependency on p later
};
void mlawNonLinearTVP::getPlasticPotential(const double& phiP, const double& phiE, double& Px) const{
Px = pow(phiE,2) + _b*pow(phiP,2);
};
void mlawNonLinearTVP::getYieldCoefficientTempDers(const IPNonLinearTVP *q1, const double& T, fullVector<double>& DcoeffsDT, fullVector<double>& DDcoeffsDTT) const{
double sigc(0.), Hc(0.), dsigcdT(0.), ddsigcdTT(0.);
sigc = q1->_ipCompression->getR();
Hc = q1->_ipCompression->getDR();
dsigcdT = q1->_ipCompression->getDRDT();
ddsigcdTT = q1->_ipCompression->getDDRDTT();
double sigt(0.), Ht(0.), dsigtdT(0.), ddsigtdTT(0.);
sigt = q1->_ipTraction->getR();
Ht = q1->_ipTraction->getDR();
dsigtdT = q1->_ipTraction->getDRDT();
ddsigtdTT = q1->_ipTraction->getDDRDTT();
DcoeffsDT.resize(3);
DDcoeffsDTT.resize(3);
double m = sigt/sigc;
double DmDT = (dsigtdT*sigc- sigt*dsigcdT)/(sigc*sigc);
double DDmDTT = (ddsigtdTT*sigc - sigt*ddsigcdTT)/(sigc*sigc) - 2*(dsigtdT*sigc- sigt*dsigcdT)/(sigc*sigc*sigc)*dsigcdT;
double Da1Dm = 3./sigc*(_n*pow(m,_n-1.)/(m+1.) - (pow(m,_n)-1.)/(m+1.)/(m+1.));
double DDa1Dmm = 3./sigc*(_n*(_n-1)*pow(m,_n-2.)/(m+1.) - _n*(pow(m,_n)-1.)/(m+1.)/(m+1.) - _n*pow(m,_n-1)/(m+1.)/(m+1.) + 2.*(pow(m,_n)-1.)/(m+1.)/(m+1.)/(m+1.));
double term = _n*pow(m,_n) + _n*pow(m,_n-1.) - pow(m,_n);
double dterm = pow(_n,2)*pow(m,_n-1.) + _n*(_n-1)*pow(m,_n-2.) - _n*pow(m,_n-1);
DcoeffsDT(2) = -_n*pow(sigc,-_n-1.)*dsigcdT;
DcoeffsDT(1) = Da1Dm*DmDT - 3.*(pow(m,_n)-1.)/(m+1.)/(sigc*sigc)*dsigcdT;
DcoeffsDT(0) = ((_n*pow(m,_n-1)+1.)/(m+1) - (pow(m,_n)+m)/(m+1.)/(m+1.))*DmDT;
DDcoeffsDTT(2) = -_n*pow(sigc,-_n-1.)*ddsigcdTT -_n*(-_n-1)*pow(sigc,-_n-2.)*dsigcdT*dsigcdT;
DDcoeffsDTT(1) = -1/sigc*dsigcdT*DcoeffsDT(1) + 3/sigc*( (dterm/(m+1.)/(m+1.) -2*term/(m+1.)/(m+1.)/(m+1.))*DmDT*DmDT + term/(m+1.)/(m+1.)*DDmDTT
+ 1/sigc*( (1/sigc*pow(dsigtdT,2) - ddsigcdTT) * (pow(m,_n)-1.)/(m+1.) + dsigtdT * term/(m+1.)/(m+1.) ) );
DDcoeffsDTT(0) = ((_n*pow(m,_n-1)+1.)/(m+1) - (pow(m,_n)+m)/(m+1.)/(m+1.))*DDmDTT
+ ((_n*(_n-1)*pow(m,_n-2))/(m+1) + 2*(pow(m,_n)+m)/(m+1.)/(m+1.)/(m+1.))*DmDT*DmDT;
};
void mlawNonLinearTVP::getFreeEnergyTVM(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double& T0, const double& T,
double *psiTVM, double *DpsiTVMdT, double *DDpsiTVMdTT) const{
// get some Things
const double& Gamma = q1->_Gamma; const double& DgammaDT = q1->_DgammaDT;
const double& dGammaDT = q1->_DGammaDT;
const STensor3& Me = q1->_ModMandel;
const STensor3& X = q1->_backsig;
const STensor3& dXdT = q1->_DbackSigDT;
STensor3 devMe, devX, dDevXdT, dDevMeDT;
double trMe, trX, dTrXdT, dTrMeDT;
STensorOperation::decomposeDevTr(q1->_ModMandel,devMe,trMe);
STensorOperation::decomposeDevTr(q1->_backsig,devX,trX);
STensorOperation::decomposeDevTr(q1->_DbackSigDT,dDevXdT,dTrXdT);
STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
// Start
double psiVE, psiVP, psiT = 0.;
double DpsiVEdT, DpsiVPdT, DpsiTdT = 0.;
double DDpsiVEdTT, DDpsiVPdTT, DpsiTdTT = 0.;
double CpT, DCpDT, DDCpDTT;
mlawNonLinearTVE::getCp(CpT,T,&DCpDT,&DDCpDTT);
// psiT
psiT = CpT*((T-_Tinitial)-T*log(T/_Tinitial));
DpsiTdT = -CpT*log(T/_Tinitial) - psiT/CpT*DCpDT;
DpsiTdTT = -CpT/T - psiT/CpT*DDCpDTT;
//
// psiVP
// get R, dRdT
fullVector<double> a(3), Da(3), DaDT(3), DDaDTT(3);
getYieldCoefficients(q1,a);
getYieldCoefficientDerivatives(q1,q1->_nup,Da);
getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
double sigc(0.), Hc(0.), dsigcdT(0.), ddsigcdTT(0.);
sigc = q1->_ipCompression->getR();
Hc = q1->_ipCompression->getDR();
dsigcdT = q1->_ipCompression->getDRDT();
ddsigcdTT = q1->_ipCompression->getDDRDTT();
double R1(0.), R2(0.), R3(0.), dR1dT(0.), dR2dT(0.), dR3dT(0.), ddR1dTT(0.), ddR2dTT(0.), ddR3dTT(0.);
const std::vector<double>& IsoHardForce = q1->_IsoHardForce;
const std::vector<double>& dIsoHardForcedT = q1->_dIsoHardForcedT;
const std::vector<double>& ddIsoHardForcedTT = q1->_ddIsoHardForcedTT;
R1 = IsoHardForce[0]; R2 = IsoHardForce[1]; R3 = IsoHardForce[2];
dR1dT = dIsoHardForcedT[0]; dR2dT = dIsoHardForcedT[1]; dR3dT = dIsoHardForcedT[2];
ddR1dTT = ddIsoHardForcedTT[0]; ddR2dTT = ddIsoHardForcedTT[1]; ddR3dTT = ddIsoHardForcedTT[2];
// get Kinematic Hardening stuff
double Hb(0.), dHbdT(0.), ddHbddT(0.);
if (q1->_ipKinematic != NULL){
Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
dHbdT = q1->_ipKinematic->getDRDT();
ddHbddT = q1->_ipKinematic->getDDRDDT();
}
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
static STensor3 alpha;
alpha = X;
alpha *= 1/(pow(kk,2)*Hb);
// ddXdTT
static STensor3 ddXdTT; // ddGammaDTT and ddQDTT are very difficult to estimate
STensorOperation::zero(ddXdTT);
double delT = T-T0;
if (delT>1e-10){
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){ ddXdTT(i,j) += 2*dXdT(i,j)/(T-T0) - 2*(X(i,j)-q0->_backsig(i,j))/pow((T-T0),2);
}
}
static STensor3 dAlphadT, ddAlphadTT;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dAlphadT(i,j) = dXdT(i,j)*1/(pow(kk,2)*Hb) - X(i,j)*1/(pow(kk*Hb,2))*dHbdT;
ddAlphadTT(i,j) = ddXdTT(i,j)*1/(pow(kk,2)*Hb) - 2*dXdT(i,j)*1/(pow(kk*Hb,2))*dHbdT
- X(i,j)*( -2/(pow(kk*Hb,3))*dHbdT*dHbdT + 1/(pow(kk*Hb,2))*ddHbddT);
}
// finally
STensorOperation::doubleContractionSTensor3(X,alpha,psiVP);
psiVP *= 0.5;
psiVP += 0.5*( 1/_HR[0]*R1*a(1)*sigc + 1/_HR[1]*pow(R2,2) + 1/_HR[2]*pow(R3,2) );
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
DpsiVPdT += 0.5*(dXdT(i,j)*alpha(i,j)+X(i,j)*dAlphadT(i,j));
DDpsiVPdTT += 0.5*(ddXdTT(i,j)*alpha(i,j) + 2*dXdT(i,j)*dAlphadT(i,j) + X(i,j)*ddAlphadTT(i,j));
}
DpsiVPdT += 0.5*( 1/_HR[0]*(dR1dT*a(1)*sigc + R1*(DaDT(1)*sigc + a(1)*dsigcdT)) + 2/_HR[1]*R2*dR2dT + 2/_HR[2]*R3*dR3dT );
DDpsiVPdTT += 0.5*( 1/_HR[0]*(ddR1dTT*a(1)*sigc + 2*dR1dT*(DaDT(1)*sigc + a(1)*dsigcdT) + R1*(DDaDTT(1)*sigc + 2*DaDT(1)*dsigcdT + a(1)*ddsigcdTT))
+ 2/_HR[1]*R2*ddR2dTT + 2/_HR[2]*R3*ddR3dTT + 2/_HR[1]*dR2dT*dR2dT + 2/_HR[2]*dR3dT*dR3dT); // CHANGE!! 2nd derivative of R and X are difficult
//
// 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);
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);
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);
DdevKeinfDT = devKeinf;
DdevKeinfDT *= dGdT/GT;
DDdevKeinfDTT = devKeinf;
DDdevKeinfDTT *= ddGdTT/GT;
const std::vector<STensor3>& devOi = q1->getConstRefToDevViscoElasticOverStress();
const std::vector<double>& trOi = q1->getConstRefToTrViscoElasticOverStress();
const std::vector<STensor3>& devDOiDT = q1->getConstRefToDevDOiDT();
const std::vector<double>& trDOiDT = q1->getConstRefToTrDOiDT();
const std::vector<STensor3>& devDDOiDTT = q1->getConstRefToDevDDOiDTT();
const std::vector<double>& trDDOiDTT = q1->getConstRefToTrDDOiDTT();
psiVE = mlawNonLinearTVE::freeEnergyMechanical(*q1,T);
DpsiVEdT = 1/9 * (1/KT * trKeinf * DtrKeinfDT - 1/(2*pow(KT,2)) * dKdT * pow(trKeinf,2));
DDpsiVEdTT = 1/9 * (1/KT * (pow(DtrKeinfDT,2)+trKeinf*DDtrKeinfDTT) - 1/pow(KT,2)*dKdT*trKeinf*DtrKeinfDT
- 1/(2*pow(KT,2)) *(ddKdTT*pow(trKeinf,2) + 2*dKdT*trKeinf*DtrKeinfDT)
+ 1/pow(KT,3)*pow(dKdT,2)*pow(trKeinf,2));
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
DpsiVEdT += 2/GT * devKeinf(i,j)*DdevKeinfDT(i,j) - 1/pow(GT,2) * dGdT * devKeinf(i,j)*devKeinf(i,j);
DDpsiVEdTT += 2/GT*(DdevKeinfDT(i,j)*DdevKeinfDT(i,j) + devKeinf(i,j)*DDdevKeinfDTT(i,j) - dGdT*devKeinf(i,j)*DdevKeinfDT(i,j))
- 1/pow(GT,2) * (ddGdTT*devKeinf(i,j)*devKeinf(i,j) + 2*dGdT*devKeinf(i,j)*DdevKeinfDT(i,j));
}
if ((_Ki.size() > 0) or (_Gi.size() > 0)){
for (int k=0; k<_Gi.size(); k++){
DpsiVEdT += 1/9 * 1/_Ki[k] * trOi[k] * trDOiDT[k];
DDpsiVEdTT += 1/9 * 1/_Ki[k] * (trOi[k]*trDDOiDTT[k] + pow(trDOiDT[k],2));
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){ DpsiVEdT += 2/_Gi[k] * devOi[k](i,j) * devDOiDT[k](i,j);
DDpsiVEdTT += 2/_Gi[k] * (devOi[k](i,j) * devDDOiDTT[k](i,j) + pow(devDOiDT[k](i,j),2));
}
}
}
// FINALLY
if(psiTVM!=NULL){
*psiTVM = psiT + psiVE + psiVP;
}
if(DpsiTVMdT!=NULL){
*DpsiTVMdT = DpsiTdT + DpsiVEdT + DpsiVPdT;
}
if(DDpsiTVMdTT!=NULL){
*DDpsiTVMdTT = DpsiTdTT + DDpsiVEdTT + DDpsiVPdTT;
}
};
void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double T0, const double& T,
const STensor3& DphiPDF, std::vector<double>* ddRdTT, std::vector<STensor3>* ddRdFdT) const{
const double& Gamma_0 = q0->_Gamma;
const double& Gamma = q1->_Gamma;
const double& DgammaDT = q1->_DgammaDT;
const double& dGammaDT_0 = q0->_DGammaDT;
const double& dGammaDT = q1->_DGammaDT; // = 0. for debugging //DEBUGGING
const STensor3& DgammaDF = q1->_DgammaDF;
const STensor3& DGammaDF = q1->_DGammaDF;
const double& gamma0 = q0->_epspbarre;
const double& gamma1 = q1->_epspbarre;
static STensor3 devMe, devX1, dDevXdT, dDevMeDT;
static double trMe, trX1, dTrXdT, dTrMeDT;
static double trMe_0, trX1_0, dTrXdT_0, dTrMeDT_0;
STensorOperation::decomposeDevTr(q1->_ModMandel,devMe,trMe);
STensorOperation::decomposeDevTr(q0->_ModMandel,devMe,trMe_0);
STensorOperation::decomposeDevTr(q1->_backsig,devX1,trX1);
STensorOperation::decomposeDevTr(q0->_backsig,devX1,trX1_0);
STensorOperation::decomposeDevTr(q1->_DbackSigDT,dDevXdT,dTrXdT);
STensorOperation::decomposeDevTr(q0->_DbackSigDT,dDevXdT,dTrXdT_0);
STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
STensorOperation::decomposeDevTr(q0->_DModMandelDT,dDevMeDT,dTrMeDT_0);
trMe /= 3; trX1 /= 3; dTrMeDT /= 3; dTrXdT /= 3;
trMe_0 /= 3; trX1_0 /= 3; dTrMeDT_0 /= 3; dTrXdT_0 /= 3;
// dTrMeDT = 0.; dTrXdT = 0.; // for debugging only //DEBUGGING
// get Dgamma
double Dgamma = gamma1 - gamma0;
// get R
double eta(0.),Deta(0.),DetaDT(0.),DDetaDTT(0.);
fullVector<double> a(3), Da(3), DaDT(3), DDaDTT(3);
getYieldCoefficients(q1,a);
getYieldCoefficientDerivatives(q1,q1->_nup,Da);
getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
if (_viscosity != NULL)
_viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT,DDetaDTT);
double etaOverDt = eta/this->getTimeStep();
double viscoTerm = etaOverDt*Gamma;
const std::vector<double>& R0 = q0->_IsoHardForce;
std::vector<double>& R = q1->_IsoHardForce;
const std::vector<double>& dRdT_0 = q0->_dIsoHardForcedT;
std::vector<double>& dRdT = q1->_dIsoHardForcedT;
R.resize(3,0.); dRdT.resize(3,0.);
double sigc(0.), Hc(0.), dsigcdT(0.), ddsigcdTT(0.);
sigc = q1->_ipCompression->getR();
Hc = q1->_ipCompression->getDR();
dsigcdT = q1->_ipCompression->getDRDT();
ddsigcdTT = q1->_ipCompression->getDDRDTT();
R[0] = - a(1) * (trMe-trX1) * sigc;
dRdT[0] = - DaDT(1)*(trMe-trX1)*sigc - a(1)*( (dTrMeDT-dTrXdT)*sigc + (trMe-trX1) * dsigcdT);
R[1] = - a(0) * sigc;
dRdT[1] = - (DaDT(0) * sigc + a(0) * dsigcdT);
if (Gamma>0 and etaOverDt>0){
R[2] = - pow(eta*Gamma/this->getTimeStep(),_p) * sigc;
dRdT[2] = - _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(DetaDT*Gamma + eta*dGammaDT)/this->getTimeStep() * sigc - pow(eta*Gamma/this->getTimeStep(),_p) * dsigcdT;
}
// dRdF
// static STensor3 dRdF;
const std::vector<STensor3>& dRdF_0 = q0->_dIsoHardForcedF;
std::vector<STensor3>& dRdF = q1->_dIsoHardForcedF;
dRdF.resize(3,0.);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dRdF[0](i,j) = -Da(1)*(trMe-trX1)*DgammaDF(i,j)*sigc - a(1)*(DphiPDF(i,j)*sigc + (trMe-trX1)*Hc*DgammaDF(i,j));
dRdF[1](i,j) = -Da(0)*DgammaDF(i,j)*sigc - a(0)*Hc*DgammaDF(i,j);
}
if (Gamma>0 and etaOverDt>0){
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dRdF[2](i,j) = - _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(eta*DGammaDF(i,j)/this->getTimeStep())*sigc - pow(eta*Gamma/this->getTimeStep(),_p) * Hc*DgammaDF(i,j);
}
}
// ddRdTT
double delT = T-T0;
if(ddRdTT!=NULL){
ddRdTT->at(0) = - DDaDTT(1)*(trMe-trX1)*sigc - a(1)*( 2*(dTrMeDT-dTrXdT)*dsigcdT + (trMe-trX1) * ddsigcdTT );
ddRdTT->at(1) = -( DDaDTT(0)*sigc + 2*DaDT(0)*dsigcdT + a(0)*ddsigcdTT);
if (Gamma>0 and etaOverDt>0){
ddRdTT->at(2) = - _p*(_p-1)*pow(eta*Gamma/this->getTimeStep(),_p-2)*pow(((DetaDT*Gamma + eta*dGammaDT)/this->getTimeStep()),2) * sigc
- _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*((DDetaDTT*Gamma + 2*DetaDT*dGammaDT)/this->getTimeStep() * sigc +
2*(DetaDT*Gamma + eta*dGammaDT)/this->getTimeStep()* dsigcdT)
- pow(eta*Gamma/this->getTimeStep(),_p) * ddsigcdTT;
}
if (delT>1e-10){
double DphiPDT = dTrMeDT- dTrXdT; double DphiPDT_0 = dTrMeDT_0 - dTrXdT_0;
double DDphiPDTT = (DphiPDT-DphiPDT_0)/delT; // (DphiPDT - (trMe-trMe_0+trX1_0-trX1)/delT)/delT; // 2*DphiPDT/delT - 2*(trMe-trMe_0+trX1_0-trX1)/pow(delT,2); // DEBUGGING
ddRdTT->at(0) -= a(1) * DDphiPDTT * sigc;
if (Gamma>0 and etaOverDt>0){
double DDGammaDTT = (dGammaDT-dGammaDT_0)/delT; // 2*dGammaDT/delT - 2*(Gamma-Gamma_0)/pow(delT,2);
ddRdTT->at(2) -= _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(eta*DDGammaDTT)/this->getTimeStep() * sigc;
}
}
q1->_ddIsoHardForcedTT = *ddRdTT;
}
// ddRdTT = 2*dRdT/(T-T0) - 2*(R-R0)/pow((T-T0),2); // ddGammaDTT is very difficult to estimate
// ddRdTF
if(ddRdFdT!=NULL){
for (int k=0; k<3; k++){
STensorOperation::zero(ddRdFdT->at(k));
if(delT>1e-10){ ddRdFdT->at(k) = dRdF[k];
ddRdFdT->at(k) -= dRdF_0[k];
ddRdFdT->at(k) *= 1/delT;
}
}
}
};
void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double T0, const double& T,
const STensor3& dFpdT, const STensor43& dFpdF, const STensor3& DphiPDF,
double *Wm, STensor3 *dWmdF, double *dWmdT) const{
const double& Gamma = q1->_Gamma;
const double& DgammaDT = q1->_DgammaDT;
const double& dGammaDT = q1->_DGammaDT;
const STensor3& DgammaDF = q1->_DgammaDF;
const STensor3& DGammaDF = q1->_DGammaDF;
const double& gamma0 = q0->_epspbarre;
const double& gamma1 = q1->_epspbarre;
const STensor3& Fp0 = q0->_Fp;
const STensor3& Fp1 = q1->_Fp;
const STensor3& Me = q1->_ModMandel;
const STensor3& X1 = q1->_backsig;
const STensor3& X0 = q0->_backsig;
const STensor3& dXdT = q1->_DbackSigDT; // Just for debugging = 0. //DEBUGGING
const STensor3& dXdT_0 = q0->_DbackSigDT;
const STensor43& dXdF_0 = q0->_DbackSigDF;
const STensor43& dXdF = q1->_DbackSigDF;
const STensor3& dMeDT = q1->_DModMandelDT;
const STensor3& dMeDT_0 = q0->_DModMandelDT;
const STensor43& dMedF = q1->_DModMandelDF;
const STensor43& dMedF_0 = q0->_DModMandelDF;
static STensor3 devMe, devX1, dDevXdT, dDevMeDT;
static double trMe, trX1, dTrXdT, dTrMeDT;
STensorOperation::decomposeDevTr(q1->_ModMandel,devMe,trMe);
STensorOperation::decomposeDevTr(q1->_backsig,devX1,trX1);
STensorOperation::decomposeDevTr(q1->_DbackSigDT,dDevXdT,dTrXdT);
STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
// get Dgamma
double Dgamma = gamma1 - gamma0;
// get Hb
double Hb(0.), dHb(0), dHbdT(0.), ddHbddT(0.);
if (q1->_ipKinematic != NULL){
Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
dHb = q1->_ipKinematic->getDDR();
dHbdT = q1->_ipKinematic->getDRDT();
ddHbddT = q1->_ipKinematic->getDDRDDT();
}
// get Dp
static STensor3 Fpinv, Fp_dot, Dp;
STensorOperation::inverseSTensor3(Fp1,Fpinv);
Fp_dot = Fp1;
Fp_dot -= Fp0;
Fp_dot *= 1/this->getTimeStep();
STensorOperation::multSTensor3(Fp_dot,Fpinv,Dp);
// get alpha
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
static STensor3 alpha0, alpha1;
alpha0 = X0; alpha1 = X1;
alpha0 *= 1/(pow(kk,2)*Hb); alpha1 *= 1/(pow(kk,2)*Hb); // how to get alpha0-current or previous Hb??? WHAT?
alpha1 -= alpha0;
alpha1 *= 1/this->getTimeStep();
// get TempX - convenient backStress tensor STensor3 tempX;
tempX = X1;
tempX -= T*dXdT;
// get R, dRdT, dRdF
std::vector<double> ddIsoHardForcedTT; ddIsoHardForcedTT.resize(3,0.);
std::vector<STensor3> ddIsoHardForcedFdT; ddIsoHardForcedFdT.resize(3,0.);
getIsotropicHardeningForce(q0,q1,T0,T,DphiPDF,&ddIsoHardForcedTT,&ddIsoHardForcedFdT);
const std::vector<double>& IsoHardForce = q1->_IsoHardForce;
const std::vector<double>& dIsoHardForcedT = q1->_dIsoHardForcedT;
const std::vector<STensor3>& dIsoHardForcedF = q1->_dIsoHardForcedF;
double R(0.), dRdT(0.), ddRdTT(0.);
STensor3 dRdF, ddRdFdT; STensorOperation::zero(dRdF); STensorOperation::zero(ddRdFdT);
for (int i=0; i<3; i++){
R += IsoHardForce[i];
dRdT += dIsoHardForcedT[i];
ddRdTT += ddIsoHardForcedTT[i];
dRdF += dIsoHardForcedF[i];
ddRdFdT += ddIsoHardForcedFdT[i]; //DEBUGGING
}
R -= T*dRdT;
// mechSourceTVP
if(Wm!=NULL){
double mechSourceTVP = 0.;
if (this->getTimeStep() > 0){
mechSourceTVP += STensorOperation::doubledot(Me,Dp);
mechSourceTVP -= STensorOperation::doubledot(tempX,alpha1);
mechSourceTVP -= (R*Dgamma)/this->getTimeStep();
}
*Wm = mechSourceTVP;
// Second Term TBD - need dXdT, dRdT -> added above
}
// Some Generic Conversion Tensors
static STensor43 dDpdFp, dFpinvdFp;
static STensor3 dFpinvdT;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dFpinvdFp(i,j,k,l) = 0.;
for (int m=0; m<3; m++)
for (int s=0; s<3; s++)
dFpinvdFp(i,j,k,l) -= Fpinv(i,m)*_I4(m,s,k,l)*Fpinv(s,j);
}
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dDpdFp(i,j,k,l) = 0.;
for (int m=0; m<3; m++)
dDpdFp(i,j,k,l) += 1/this->getTimeStep()*_I4(i,m,k,l)*Fpinv(m,j) + Fp_dot(i,m)*dFpinvdFp(m,j,k,l);
}
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dFpinvdT(i,j) = 0.;
for (int m=0; m<3; m++)
for (int s=0; s<3; s++)
dFpinvdT(i,j) -= Fpinv(i,m)*dFpdT(m,s)*Fpinv(s,j);
}
// dDpdF and dDpdT static STensor43 dDpdF;
// STensorOperation::zero(dDpdF);
STensorOperation::multSTensor43(dDpdFp, dFpdF, dDpdF);
static STensor3 dDpdT;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dDpdT(i,j) = 0.;
for (int m=0; m<3; m++){
dDpdT(i,j) += dFpdT(i,m)*Fpinv(m,j)*1/this->getTimeStep() + Fp_dot(i,m)*dFpinvdT(m,j); // dDpdFp(i,j,k,l)*dFpdT(k,l); // maybe use Dp = Fp_dot Fpinv
}}
// dXdF is available
// const STensor43 dXdF = q1->_DbackSigDF;
// ddXdTF -> make it -- CHANGE!!
static STensor43 ddXdTF; // ddGammaDTF and ddQDTF are very difficult to estimate
static STensor3 ddXdTT; // ddGammaDTT and ddQDTT are very difficult to estimate
STensorOperation::zero(ddXdTF);
STensorOperation::zero(ddXdTT);
double delT = T-T0;
if (delT>1e-10){
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
ddXdTT(i,j) += (dXdT(i,j)-dXdT_0(i,j))/delT; // 2*dXdT(i,j)/(T-T0) - 2*(X1(i,j)-X0(i,j))/pow((T-T0),2); // = 0. in DEBUGGING
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
ddXdTF(i,j,k,l) += (dXdF(i,j,k,l)-dXdF_0(i,j,k,l))/delT; // = 0. in DEBUGGING
}
}
}
// dAlphadF - from dXdF
static STensor43 dAlphadF;
STensorOperation::zero(dAlphadF);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dAlphadF(i,j,k,l) += dXdF(i,j,k,l)*1/(pow(kk,2)*Hb) + X1(i,j)*DgammaDF(k,l)*1/(pow(kk*Hb,2))*dHb;
}
// dAlphadT - from dXdT
static STensor3 dAlphadT;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dAlphadT(i,j) = dXdT(i,j)*1/(pow(kk,2)*Hb) + X1(i,j)*DgammaDT*1/(pow(kk*Hb,2))*dHbdT;
}
// dMedF is available
// const STensor43 dMedF = q1->_DModMandelDF;
if(dWmdF!=NULL){
STensor3 dmechSourceTVPdF(0.);
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();
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)*dDpdF(i,j,k,l)
- (dXdF(i,j,k,l) - T*ddXdTF(i,j,k,l))*alpha1(i,j) - tempX(i,j)*dAlphadF(i,j,k,l)/this->getTimeStep();
}
}
}
*dWmdF = dmechSourceTVPdF;
}
if(dWmdT!=NULL){ double dmechSourceTVPdT(0.);
if (this->getTimeStep() > 0){
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dmechSourceTVPdT += dMeDT(i,j)*Dp(i,j) + Me(i,j)*dDpdT(i,j)
+ T*ddXdTT(i,j)*alpha1(i,j) - tempX(i,j)*dAlphadT(i,j)/this->getTimeStep()
+ T*ddRdTT*Dgamma/this->getTimeStep() - R*DgammaDT/this->getTimeStep() ;
}
}
*dWmdT = dmechSourceTVPdT;
}
};
// NEW
void mlawNonLinearTVP::constitutive( // This function is just a placeholder, defined in the pure virtual class -> does nothing, its never called.
const STensor3& F0, // previous deformation gradient (input @ time n)
const STensor3& F1, // current deformation gradient (input @ time n+1)
STensor3 &P1, // current 1st Piola-Kirchhoff stress tensor (output)
const IPVariable *q0i, // array of previous internal variables
IPVariable *q1i, // current array of internal variable (in ipvcur on output),
STensor43 &Tangent, // tangents (output)
const bool stiff, // if true compute the tangents
STensor43* elasticTangent,
const bool dTangent,
STensor63* dCalgdeps) const{
static SVector3 gradT, temp2;
static STensor3 temp3;
static STensor33 temp33;
static double tempVal;
static STensor43 dFpdF, dFedF;
static STensor3 dFpdT, dFedT;
const IPNonLinearTVP *q0 = dynamic_cast<const IPNonLinearTVP *> (q0i);
IPNonLinearTVP *q1 = dynamic_cast<IPNonLinearTVP *>(q1i);
if(elasticTangent != NULL) Msg::Error("mlawNonLinearTVP elasticTangent not defined");
predictorCorrector_ThermoViscoPlastic(F0,F1,P1,q0,q1,Tangent,dFpdF, dFpdT, dFedF, dFedT,
_Tref,_Tref,gradT,gradT,temp2,temp3,temp3,temp2,temp33,
tempVal,tempVal,temp3,
tempVal,tempVal,temp3,
stiff, NULL);
};
void mlawNonLinearTVP::constitutive(
const STensor3& F0, // initial deformation gradient (input @ time n)
const STensor3& F1, // updated deformation gradient (input @ time n+1)
STensor3& P, // updated 1st Piola-Kirchhoff stress tensor (output)
const IPVariable *q0i, // array of initial internal variable
IPVariable *q1i, // updated array of internal variable (in ipvcur on output),
STensor43 &Tangent, // constitutive mechanical tangents (output)
const double& T0, // previous temperature
const double& T, // temperature
const SVector3 &gradT0, // previous temperature gradient
const SVector3 &gradT, // temperature gradient
SVector3 &fluxT, // temperature flux
STensor3 &dPdT, // mechanical-thermal coupling
STensor3 &dfluxTdgradT, // thermal tengent
SVector3 &dfluxTdT,
STensor33 &dfluxTdF, // thermal-mechanical coupling
double &thermalSource, // - Cp*dTdt
double &dthermalSourcedT, // thermal source
STensor3 &dthermalSourcedF,
double &mechanicalSource, // mechanical source--> convert to heat
double &dmechanicalSourcedT,
STensor3 &dmechanicalSourceF,
const bool stiff,
STensor43* elasticTangent) const{
const IPNonLinearTVP *q0 = dynamic_cast<const IPNonLinearTVP *>(q0i); IPNonLinearTVP *q1 = dynamic_cast<IPNonLinearTVP *>(q1i);
// Declaring here, coz I dont know better.
static STensor43 dFedF, dFpdF;
static STensor3 dFpdT, dFedT;
if (!_tangentByPerturbation){
this->predictorCorrector_ThermoViscoPlastic(F0, F1, P, q0, q1, Tangent, dFpdF, dFpdT, dFedF, dFedT,
T0, T, gradT0, gradT, fluxT, dPdT, dfluxTdgradT, dfluxTdT, dfluxTdF,
thermalSource, dthermalSourcedT, dthermalSourcedF,
mechanicalSource, dmechanicalSourcedT, dmechanicalSourceF,
stiff, elasticTangent);
}
else{
this->predictorCorrector_ThermoViscoPlastic(F0, F1, P, q0, q1, T0, T, gradT0, gradT, fluxT, thermalSource, mechanicalSource, elasticTangent);
if (stiff)
this->tangent_full_perturbation(F0,F1,P,q0,q1,T0, T, gradT0, gradT, fluxT, thermalSource, mechanicalSource,
Tangent, dFpdF, dFpdT, dFedF, dFedT,
dPdT, dfluxTdgradT, dfluxTdT, dfluxTdF,
dthermalSourcedT, dthermalSourcedF,
dmechanicalSourcedT, dmechanicalSourceF);
}
};
void mlawNonLinearTVP::predictorCorrector_ThermoViscoPlastic(
const STensor3& F0, // initial deformation gradient (input @ time n)
const STensor3& F1, // updated deformation gradient (input @ time n+1)
STensor3 &P, // updated 1st Piola-Kirchhoff stress tensor (output)
const IPNonLinearTVP *q0, // array of initial internal variable
IPNonLinearTVP *q1, // updated array of internal variable (in ipvcur on output),
const double& T0, // previous temperature
const double& T, // temperature
const SVector3 &gradT0, // previous temeprature gradient
const SVector3 &gradT, // temeprature gradient
SVector3 &fluxT, // temperature flux)
double &thermalSource,
double& mechanicalSource,
STensor43* elasticTangent) const{
// temp variables
static STensor43 Tangent, dFpdF, dFedF;
static STensor3 dPdT, dfluxTdgradT, dthermalSourcedF, dmechanicalSourceF, dFpdT, dFedT;
static STensor33 dfluxTdF;
static SVector3 dfluxTdT;
static double dthermalSourcedT, dmechanicalSourcedT;
predictorCorrector_ThermoViscoPlastic(F0,F1,P,q0,q1,Tangent,dFpdF,dFpdT,dFedF,dFedT,
T0,T,gradT0,gradT,fluxT,dPdT,dfluxTdgradT,dfluxTdT,dfluxTdF,
thermalSource,dthermalSourcedT,dthermalSourcedF,
mechanicalSource,dmechanicalSourcedT,dmechanicalSourceF,false,elasticTangent);
};
void mlawNonLinearTVP::tangent_full_perturbation(
const STensor3& F0,
const STensor3& F1, // updated deformation gradient (input @ time n+1)
STensor3& P, // updated 1st Piola-Kirchhoff stress tensor (output)
const IPNonLinearTVP *q0, // array of initial internal variable
IPNonLinearTVP *q1, // updated array of internal variable (in ipvcur on output),
const double& T0, // previous temperature
const double& T, // temperature
const SVector3 &gradT0, // previous temeprature gradient
const SVector3 &gradT, // temeprature gradient
SVector3 &fluxT, // temperature flux)
double &thermalSource,
double &mechanicalSource,
STensor43 &Tangent, // mechanical tangents (output)
STensor43 &dFpdF, // plastic tangent STensor3 &dFpdT, // plastic tangent
STensor43 &dFedF, // elastic tangent
STensor3 &dFedT, // elastic tangent
STensor3 &dPdT, // mechanical-thermal coupling
STensor3 &dfluxTdgradT, // thermal tengent
SVector3 &dfluxTdT,
STensor33 &dfluxTdF, // thermal-mechanical coupling
double &dthermalSourcedT, // thermal source
STensor3 &dthermalSourcedF,
double &dmechanicalSourcedT,
STensor3 &dmechanicalSourceF) const{
static STensor3 Fplus, Pplus;
static SVector3 fluxTPlus, gradTplus;
static double thermalSourcePlus;
static double mechanicalSourcePlus;
static IPNonLinearTVP qPlus(*q0);
// perturb F
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
Fplus = F1;
Fplus(i,j) += _perturbationfactor;
predictorCorrector_ThermoViscoPlastic(F0,Fplus,Pplus,q0,&qPlus,T0,T,gradT0,gradT,fluxTPlus,thermalSourcePlus,mechanicalSourcePlus, NULL);
q1->_DgammaDF(i,j) = (qPlus._epspbarre - q1->_epspbarre)/_perturbationfactor;
q1->_DirreversibleEnergyDF(i,j) = (qPlus._irreversibleEnergy - q1->_irreversibleEnergy)/_perturbationfactor;
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
Tangent(k,l,i,j) = (Pplus(k,l) - P(k,l))/(_perturbationfactor);
dFpdF(k,l,i,j) = (qPlus._Fp(k,l)-q1->_Fp(k,l))/_perturbationfactor;
dFedF(k,l,i,j) = (qPlus._Fe(k,l)-q1->_Fe(k,l))/_perturbationfactor;
}
dfluxTdF(k,i,j) = (fluxTPlus(k) - fluxT(k))/_perturbationfactor;
}
q1->getRefToDElasticEnergyPartdF()(i,j)=(qPlus.getConstRefToElasticEnergyPart()-q1->getConstRefToElasticEnergyPart())/_perturbationfactor;
q1->getRefToDViscousEnergyPartdF()(i,j)=(qPlus.getConstRefToViscousEnergyPart()-q1->getConstRefToViscousEnergyPart())/_perturbationfactor;
q1->getRefToDPlasticEnergyPartdF()(i,j)=(qPlus.getConstRefToPlasticEnergyPart()-q1->getConstRefToPlasticEnergyPart())/_perturbationfactor;
dthermalSourcedF(i,j) = (thermalSourcePlus - thermalSource)/_perturbationfactor;
dmechanicalSourceF(i,j) = (mechanicalSourcePlus - mechanicalSource)/_perturbationfactor;
}
}
// perturb gradT
double gradTpert = _perturbationfactor*T0/1e-3;
for (int i=0; i<3; i++){
gradTplus = gradT;
gradTplus(i) += gradTpert;
predictorCorrector_ThermoViscoPlastic(F0,F1,Pplus,q0,&qPlus,T0,T,gradT0,gradTplus,fluxTPlus,thermalSourcePlus,mechanicalSourcePlus, NULL);
for (int k=0; k<3; k++){
dfluxTdgradT(k,i) = (fluxTPlus(k) - fluxT(k))/gradTpert;
}
}
// perturb T
double Tplus = T+ _perturbationfactor*T0;
predictorCorrector_ThermoViscoPlastic(F0,F1,Pplus,q0,&qPlus,T0,Tplus,gradT0,gradT,fluxTPlus,thermalSourcePlus,mechanicalSourcePlus, NULL);
q1->_DgammaDT = (qPlus._epspbarre - q1->_epspbarre)/(_perturbationfactor*T0);
q1->_DirreversibleEnergyDT = (qPlus._irreversibleEnergy - q1->_irreversibleEnergy)/(_perturbationfactor*T0);
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
dFpdT(k,l) = (qPlus._Fp(k,l) - q1->_Fp(k,l))/(_perturbationfactor*T0);
dFedT(k,l) = (qPlus._Fe(k,l) - q1->_Fe(k,l))/(_perturbationfactor*T0);
dPdT(k,l) = (Pplus(k,l) - P(k,l))/(_perturbationfactor*T0); }
dfluxTdT(k) = (fluxTPlus(k) - fluxT(k))/(_perturbationfactor*T0);
}
q1->getRefToDElasticEnergyPartdT() = (qPlus.getConstRefToElasticEnergyPart()-q1->getConstRefToElasticEnergyPart())/(_perturbationfactor*T0);
q1->getRefToDViscousEnergyPartdT() = (qPlus.getConstRefToViscousEnergyPart()-q1->getConstRefToViscousEnergyPart())/(_perturbationfactor*T0);
q1->getRefToDPlasticEnergyPartdT() = (qPlus.getConstRefToPlasticEnergyPart()-q1->getConstRefToPlasticEnergyPart())/(_perturbationfactor*T0);
dthermalSourcedT = (thermalSourcePlus - thermalSource)/(_perturbationfactor*T0);
dmechanicalSourcedT = (mechanicalSourcePlus - mechanicalSource)/(_perturbationfactor*T0);
};
void mlawNonLinearTVP::predictorCorrector_ThermoViscoPlastic(
const STensor3& F0, // initial deformation gradient (input @ time n)
const STensor3& F1, // updated deformation gradient (input @ time n+1)
STensor3& P, // updated 1st Piola-Kirchhoff stress tensor (output)
const IPNonLinearTVP *q0, // array of initial internal variable
IPNonLinearTVP *q1, // updated array of internal variable (in ipvcur on output),
STensor43 &Tangent, // mechanical tangents (output)
STensor43 &dFpdF, // plastic tangent
STensor3 &dFpdT, // plastic tangent
STensor43 &dFedF, // elastic tangent
STensor3 &dFedT, // elastic tangent
const double& T0, // previous temperature
const double& T, // temperature
const SVector3 &gradT0, // previoustemeprature gradient
const SVector3 &gradT, // temeprature gradient
SVector3 &fluxT, // temperature flux
STensor3 &dPdT, // mechanical-thermal coupling
STensor3 &dfluxTdgradT, // thermal tengent
SVector3 &dfluxTdT,
STensor33 &dfluxTdF, // thermal-mechanical coupling
double &thermalSource, // - cp*dTdt
double &dthermalSourcedT, // thermal source
STensor3 &dthermalSourcedF,
double &mechanicalSource, // mechanical source--> convert to heat
double &dmechanicalSourcedT,
STensor3 &dmechanicalSourceF,
const bool stiff,
STensor43* elasticTangent) const{
if (_tangentByPerturbation){
if (_nonAssociatedFlow){
this->predictorCorrector_TVP_nonAssociatedFlow(F0,F1,q0,q1,P,false,Tangent,dFedF,dFpdF,dFedT,dFpdT,
T0,T,gradT0,gradT,fluxT,dPdT,dfluxTdgradT,dfluxTdT,dfluxTdF,
thermalSource,dthermalSourcedT,dthermalSourcedF,
mechanicalSource,dmechanicalSourcedT,dmechanicalSourceF);
}
else{
this->predictorCorrector_TVP_AssociatedFlow(F1,q0,q1,P,false,Tangent,dFedF,dFpdF,T0,T);
}
// if (stiff){ // WHAT????? WHY TWICE? (The same line happens above as STIFF is always TRUE)
// tangent_full_perturbation(Tangent,dFedF,dFpdF,P,F,q0,q1);
// }
}
else{
if (_nonAssociatedFlow){
this->predictorCorrector_TVP_nonAssociatedFlow(F0,F1,q0,q1,P,stiff,Tangent,dFedF,dFpdF,dFedT,dFpdT,
T0,T,gradT0,gradT,fluxT,dPdT,dfluxTdgradT,dfluxTdT,dfluxTdF,
thermalSource,dthermalSourcedT,dthermalSourcedF,
mechanicalSource,dmechanicalSourcedT,dmechanicalSourceF);
}
else{
this->predictorCorrector_TVP_AssociatedFlow(F1,q0,q1,P,stiff,Tangent,dFedF,dFpdF,T0,T);
}
}
/*
// compute mechanical tengent
if (elasticTangent!= NULL)
{
STensor3 Ptmp(0.);
static IPNonLinearTVP q1tmp(*q1);
q1tmp.operator=(*(dynamic_cast<const IPVariable*> (q1)));
mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(q0->getConstRefToFe(), q1->getConstRefToFe(), Ptmp, q0, &q1tmp, *elasticTangent, dFedF, dFedT,
T0, T, gradT0, gradT, fluxT, dPdT, dfluxTdgradT, dfluxTdT, dfluxTdF,
thermalSource, dthermalSourcedT, dthermalSourcedF,
mechanicalSource, dmechanicalSourcedT, dmechanicalSourceF,
true);
// mechSource requires taylorQuinney Term - CHANGE
// UPDATE FOR TVP
// thermalSource derivatives -DEFINE LATER
double CpT;
if (!stiff){ // if stiff, we dont need this technically
getCp(CpT,T);
}
else{
// thermalSource TVP
CpT = -T*q1->getConstRefToDDFreeEnergyTVMdTT();
}
// thermalSource
if (this->getTimeStep() > 0.){
thermalSource = -CpT*(T-T0)/this->getTimeStep();
}
else{thermalSource = 0.;}
// mechSource - TVP
mechanicalSource += q1->getConstRefToMechSrcTVP();
dmechanicalSourcedT += q1->getConstRefTodMechSrcTVPdT();
dmechanicalSourceF += q1->getConstRefTodMechSrcTVPdF();
}*/
};
void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3& F0, const STensor3& F, const IPNonLinearTVP *q0, IPNonLinearTVP *q1,
STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, STensor3& dFedT, STensor3& dFpdT,
const double T0,
const double T,
const SVector3 &gradT0, // previous temperature gradient
const SVector3 &gradT, // temperature gradient
SVector3 &fluxT, // temperature flux
STensor3 &dPdT, // mechanical-thermal coupling
STensor3 &dfluxTdgradT, // thermal tengent
SVector3 &dfluxTdT,
STensor33 &dfluxTdF, // thermal-mechanical coupling
double &thermalSource, // - Cp*dTdt
double &dthermalSourcedT, // thermal source
STensor3 &dthermalSourcedF,
double &mechanicalSource, // mechanical source--> convert to heat
double &dmechanicalSourcedT,
STensor3 &dmechanicalSourceF) const{
// compute elastic predictor
STensor3& Fp1 = q1->_Fp;
const STensor3& Fp0 = q0->_Fp;
Fp1 = Fp0; // plastic deformation tensor
q1->_epspbarre = q0->_epspbarre; // plastic equivalent strain
q1->_epspCompression = q0->_epspCompression;
q1->_epspTraction = q0->_epspTraction;
q1->_epspShear = q0->_epspShear;
q1->_backsig = q0->_backsig; // backstress tensor
q1->_DbackSigDT = q0->_DbackSigDT; // DbackSigDT
q1->_DgammaDt = 0.;
static STensor3 Fpinv, Ce, Fepr;
STensorOperation::inverseSTensor3(Fp1,Fpinv);
STensorOperation::multSTensor3(F,Fpinv,Fepr);
STensorOperation::multSTensor3FirstTranspose(Fepr,Fepr,Ce);
// NEW
static STensor3 Cepr, Ceinvpr;
Cepr = Ce;
STensorOperation::inverseSTensor3(Cepr,Ceinvpr);
// NEW
static STensor3 invFp0; // plastic predictor
invFp0= Fpinv;
STensor3& Fe = q1->_Fe;
Fe = Fepr;
static STensor43 DlnDCepr, DlnDCe;
static STensor63 DDlnDDCe;
static STensor3 expGN;
static STensor43 dexpAdA; // estimation of dexpA/dA
STensor3& Ee = q1->_Ee;
STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCepr,&DDlnDDCe);
Ee *= 0.5;
DlnDCe = DlnDCepr;
// update A, B
double Ke, Ge = 0.; double Kde, Gde = 0.; double KTsum, GTsum = 0.;
double DKe, DGe = 0.; double DKde, DGde = 0.; double DKDTsum, DGDTsum = 0.;
ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T);
getTVEdCorKirDT(q0,q1,T0,T);
static STensor3 Me, Mepr, Kepr; // Ke is corKir
getModifiedMandel(Ce,q0,q1);
Me = q1->_ModMandel;
Mepr = Me;
Kepr = q1-> _kirchhoff;
static STensor3 devXn;
static double trXn;
STensorOperation::decomposeDevTr(q0->_backsig,devXn,trXn); // needed for chaboche model
static STensor3 PhiPr;
PhiPr = q1->_ModMandel;
PhiPr -= q1->_backsig;
static STensor3 devPhipr,devPhi; // effective dev stress predictor
static double ptildepr,ptilde;
STensorOperation::decomposeDevTr(PhiPr,devPhipr,ptildepr);
ptildepr/= 3.;
ptilde = ptildepr; // current effective pressure -> first invariant
devPhi = devPhipr;
double PhiEqpr2 = 1.5*devPhipr.dotprod();
double PhiEqpr = sqrt(PhiEqpr2); // -> second invariant
// plastic poisson ratio
q1->_nup = (9.-2.*_b)/(18.+2.*_b);
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
// double Gamma = 0.; // flow rule parameter
double& Gamma = q1->_Gamma; // flow rule parameter
double PhiEq = PhiEqpr; // current effective deviator
// hardening this->hardening(q0,q1,T);
static fullVector<double> a(3), Da(3); // yield coefficients and derivatives in respect to plastic deformation
this->getYieldCoefficients(q1,a);
double Hb(0.), dHbdT(0.), ddHbddT(0.);
if (q1->_ipKinematic != NULL){
Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
dHbdT = Hb * q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
ddHbddT = Hb * q1->_ipKinematic->getDDRDDT();
}
//a.print("a init");
static STensor3 devN; // dev part of yield normal
double trN = 0.; // trace part of yield normal
static STensor3 N; // yield normal
STensorOperation::zero(devN);
STensorOperation::zero(N);
double f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
double DfDGamma = 0.;
double dfdDgamma = 0.;
double u = 1.;
double v = 1.;
double A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
double dDgammaDGamma = 0.;
double Dgamma = 0.; // eqplastic strain
static fullVector<double> m(2);
getChabocheCoeffs(m,0,q1);
double Cx(0.), Px(0.);
getPlasticPotential(ptilde, PhiEq, Px);
Cx = m(0)*(q0->_epspbarre + Dgamma) + m(1)*Px - 1/Hb * dHbdT + 1;
double Gt = Ge + pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
double Kt = Ke + pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
if (q1->dissipationIsBlocked()){
q1->getRefToDissipationActive() = false;
}
else{
if (f>_tol){
q1->getRefToDissipationActive() = true;
// plasticity
int ite = 0;
int maxite = 100; // maximal number of iters
//Msg::Error("plasticity occurs f = %e",f);
//double f0 = fabs(f);
while (fabs(f) >_tol or ite <1){
double eta(0.),Deta(0.),DetaDT(0.);
if (_viscosity != NULL)
_viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT);
double etaOverDt = eta/this->getTimeStep();
double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A);
dDgammaDGamma = kk*(A+Gamma*dAdGamma); // mistake in the paper (VD 2016)
this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*ptilde -Da(0);
double dCxdDgamma = m(0);
double Stemp = STensorOperation::doubledot(devPhi,devXn);
double He = (PhiEq*6*Gamma*kk*kk*Hb - Stemp*3.*1./PhiEq)/(2*Cx*Cx*u) * dCxdDgamma;
double Hp = ptilde*(2*_b*Gamma*kk*kk*Hb - trXn)/(3*Cx*Cx*v) * dCxdDgamma; dfdDgamma += a(2)*_n*pow(PhiEq,(_n-1))*He - a(1)*Hp;
if (Gamma>0 and etaOverDt>0)
dfdDgamma -= _p*pow(etaOverDt,_p-1.)*Deta/this->getTimeStep()*pow(Gamma,_p); // THIS term is absent in the paper!! WHAT??
DfDGamma = dfdDgamma*dDgammaDGamma - (_n*a(2)*6.*Gt)*pow(PhiEq,_n)/u + a(1)*ptilde*2.*_b*Kt/v;
if (Gamma>0 and etaOverDt>0)
DfDGamma -=pow(etaOverDt,_p)*_p*pow(Gamma,_p-1.);
double dGamma = -f/DfDGamma; // WHAT
if (Gamma + dGamma <=0.){
Gamma /= 2.; // WHAT
}
else
Gamma += dGamma;
//Msg::Error("Gamma = %e",Gamma);
u = 1.+6.*Gt*Gamma;
v = 1.+2.*_b*Kt*Gamma;
devPhi = (devPhipr + devXn * (1-1/Cx));
devPhi *= 1./u;
PhiEq = sqrt(1.5*devPhi.dotprod());
ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v;
A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
Dgamma = kk*Gamma*A;
//Msg::Error("it = %d, u=%e, v=%e, Dgamma=%e",ite,u,v,Dgamma);
updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
hardening(q0,q1,T);
getYieldCoefficients(q1,a);
//a.print("a update");
f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
double viscoTerm = etaOverDt*Gamma;
if (Gamma>0 and etaOverDt>0) f-= pow(viscoTerm,_p);
ite++;
//if (ite> maxite-5)
//Msg::Error("it = %d, DfDGamma = %e error = %e dGamma = %e, Gamma = %e",ite,DfDGamma,f,dGamma,Gamma);
if (fabs(f) <_tol) break;
if(ite > maxite){
Msg::Error("No convergence for plastic correction in mlawNonLinearTVP nonAssociatedFlow Maxwell iter = %d, f = %e!!",ite,f);
P(0,0) = P(1,1) = P(2,2) = sqrt(-1.);
return;
}
};
q1->_DgammaDt = Dgamma/this->getTimeStep();
// update effective stress tensor
devPhi = (devPhipr + devXn * (1-1/Cx));
devPhi *= 1./u;
ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v;
// update normal
devN = devPhi;
devN *= 3.;
trN = 2.*_b*ptilde;
N = devN;
N(0,0) += trN/3.;
N(1,1) += trN/3.;
N(2,2) += trN/3.;
// estimate exp(GammaN)
// static STensor3 expGN;
static STensor3 GammaN;
GammaN = N;
GammaN *= Gamma;
STensorOperation::expSTensor3(GammaN,_order,expGN,&dexpAdA);
// update plastic deformation tensor
STensorOperation::multSTensor3(expGN,Fp0,Fp1);
// update IP
updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma); // Why again???? WHAT??
// update elastic deformation tensor, corotational stress
STensorOperation::inverseSTensor3(Fp1,Fpinv);
STensorOperation::multSTensor3(F,Fpinv,Fe);
STensorOperation::multSTensor3FirstTranspose(Fe,Fe,Ce);
STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCe,&DDlnDDCe);
Ee *= 0.5; // Here We have assumed that, Ce commutes with Q (i.e. with M)
// update A, B
ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T); //,false);
getTVEdCorKirDT(q0,q1,T0,T); // update dCorKirdT
// DKDTsum and DGDTsum = sum of bulk/shear moduli derivatives wrt T
// backstress
static STensor3 DB; // increment
DB = (N); // increment
DB *= (pow(kk,1)*Hb*Gamma); // pow(kk,2) CHANGE
DB += q0->_backsig;
DB *= 1/Cx;
DB -= q0->_backsig;
q1->_backsig += DB; // update
q1->_ModMandel = devPhi;
q1->_ModMandel += q1->_backsig;
q1->_ModMandel(0,0) += (ptilde);
q1->_ModMandel(1,1) += (ptilde);
q1->_ModMandel(2,2) += (ptilde);
}
else{
q1->getRefToDissipationActive() = false;
}
}
// update Stress
const STensor3& KS = q1->_kirchhoff;
getModifiedMandel(Ce, q0, q1); // update Mandel
const STensor3& MS = q1->_ModMandel;
static STensor3 S, Ceinv; // get 2nd PK
STensorOperation::inverseSTensor3(Ce,Ceinv);
STensorOperation::multSTensor3(Ceinv,q1->_ModMandel,S);
// first PK
for(int i=0; i<3; i++)
for(int j=0; j<3; j++){
P(i,j) = 0.;
for(int k=0; k<3; k++)
for(int l=0; l<3; l++)
P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
}
// check Commutavity
static STensor3 Check;
checkCommutavity(Check,Ce,S,q1);
// defo energy
// q1->getRefToElasticEnergy()=deformationEnergy(*q1);
q1->getRefToViscousEnergyPart()=viscousEnergy(*q0,*q1)+q0->getConstRefToViscousEnergyPart(); q1->getRefToPlasticEnergy() = q0->plasticEnergy();
if (Gamma > 0){
double dotMSN = dot(MS,N);
q1->getRefToPlasticEnergy() += Gamma*dotMSN; // = Dp:Me
}
// fluxT
double KThConT, DKThConDT; mlawNonLinearTVE::getKTHCon(KThConT,T,&DKThConDT);
double J = 1.;
STensor3 Finv(0.);
if (_thermalEstimationPreviousConfig){ // ADD _thermalEstimationPreviousConfig
STensorOperation::inverseSTensor3(F0,Finv);
J = STensorOperation::determinantSTensor3(F0);
}
else{
STensorOperation::inverseSTensor3(F,Finv);
J = STensorOperation::determinantSTensor3(F);
}
static STensor3 Cinv;
STensorOperation::multSTensor3SecondTranspose(Finv,Finv,Cinv);
STensorOperation::multSTensor3SVector3(Cinv,gradT,fluxT);
fluxT *= (-KThConT*J);
// ThermSrc and MechSrc after dPdF and dPdT - But need the following tensors for the mechSrcTVP function call
static STensor3 DphiPDF;
STensorOperation::zero(dFpdF);
STensorOperation::zero(dFpdT);
// I didnt make this
if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
q1->getRefToIrreversibleEnergy() = q1->defoEnergy();
}
else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
(this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
q1->getRefToIrreversibleEnergy() = q1->plasticEnergy();
}
else{
q1->getRefToIrreversibleEnergy() = 0.;
}
if (stiff){
// dP/dF
//1. get DtrKeprDCepr and DdevKeprDCepr, also DKeprDCepr
static STensor3 DpDCepr;
static STensor43 DdevKDCepr;
static STensor43 DdevKeprDCepr, DdevMeprDCepr, DKeprDCepr;
static STensor3 DpKeprDCepr, DpMeprDCepr; // K = corKir; pK -> pressure of corKir
STensorOperation::multSTensor3STensor43(_I,DlnDCepr,DpKeprDCepr);
DpKeprDCepr*= (0.5*Ke);
STensorOperation::multSTensor43(_Idev,DlnDCepr,DdevKeprDCepr);
DdevKeprDCepr*= Ge;
// DKeprDCepr
DKeprDCepr = DdevKeprDCepr;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DKeprDCepr(i,j,k,l) += _I(i,j)*DpKeprDCepr(k,l);
}
// initiate here - corrected corKir derivatives wrt Cepr
DpDCepr = DpKeprDCepr; DdevKDCepr = DdevKeprDCepr;
//2. get DCeprDCepr and DCeinvprDCepr
static STensor43 DCeprDCepr, DCeinvprDCepr;
DCeprDCepr = _I4;
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DCeinvprDCepr(i,s,k,l) = 0.;
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
DCeinvprDCepr(i,s,k,l) -= Ceinvpr(i,m)*_I4(m,j,k,l)*Ceinvpr(j,s);
}
//3. get DtrMeprDCepr and DdevMeprDCepr
static STensor43 DMeprDCepr;
static STensor3 temp1;
static STensor43 temp2, temp3;
STensorOperation::multSTensor3(Kepr,Ceinvpr,temp1);
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
temp2(i,s,k,l) = 0.;
temp3(i,s,k,l) = 0.;
for (int m=0; m<3; m++){
temp2(i,s,k,l) += _I4(i,m,k,l)*temp1(m,s);
temp3(i,s,k,l) += DKeprDCepr(i,m,k,l)*Ceinvpr(m,s) + Kepr(i,m)*DCeinvprDCepr(m,s,k,l);
}
}
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DMeprDCepr(i,s,k,l) = 0.;
for (int m=0; m<3; m++)
DMeprDCepr(i,s,k,l) += Cepr(i,m)*temp3(m,s,k,l);
}
DMeprDCepr += temp2;
DMeprDCepr += DKeprDCepr;
DMeprDCepr *= 0.5;
// Because, TrMepr = TrKepr
DpMeprDCepr = DpKeprDCepr;
//
STensorOperation::zero(DdevMeprDCepr);
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DdevMeprDCepr(i,s,k,l) += DMeprDCepr(i,s,k,l) - _I(i,s)*DpMeprDCepr(k,l);
}
//4. get DdevphiDCepr and DphiPprDCepr
static STensor3 DphiPprDCepr, DphiPDCepr;
static STensor43 DdevphiprDCepr, DdevphiDCepr;
STensorOperation::zero(DphiPprDCepr);
STensorOperation::zero(DphiPDCepr);
STensorOperation::zero(DdevphiprDCepr);
STensorOperation::zero(DdevphiDCepr);
DphiPprDCepr = DpMeprDCepr; DdevphiprDCepr = DdevMeprDCepr;
DdevphiDCepr = DdevphiprDCepr;
DdevphiDCepr *= 1./u;
DphiPDCepr = DphiPprDCepr;
DphiPDCepr *= 1./v;
//5. get DphiEprDCepr from DdevphiDCepr
static STensor3 DphiEDCepr, DphiEDdevPhi;
STensorOperation::zero(DphiEDdevPhi);
DphiEDdevPhi = devPhi;
DphiEDdevPhi *= 1.5/(PhiEq);
STensorOperation::multSTensor3STensor43(DphiEDdevPhi,DdevphiDCepr,DphiEDCepr);
// 6. to 11. (inside if loop-> Gamma > 0.)
static STensor3 dAdCepr, dfDCepr;
static STensor3 DGDCepr;
static STensor3 DgammaDCepr;
static STensor3 DtrNDCepr;
static STensor43 DdevNDCepr;
static STensor43 dFpDCepr;
static STensor43 dCedCepr, dCeinvdCepr;
static STensor43 dXdCepr;
static STensor3 dTrXdCepr;
if (Gamma >0){
// plastic
//6. get dAdCepr from DdevphiDCepr
double fact = a(2)*_n*pow(PhiEq,_n-1.);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
dAdCepr(i,j) = (4.*_b*_b*ptilde/(A*3.))*DphiPDCepr(i,j);
dfDCepr(i,j) = -a(1)*DphiPDCepr(i,j);
dAdCepr(i,j) += (6./A)*PhiEq*DphiEDCepr(i,j);
dfDCepr(i,j) += fact*DphiEDCepr(i,j);
}
}
//7. get DGDCepr, DgammaDCepr
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
DGDCepr(i,j) = (-dfDCepr(i,j)-dfdDgamma*kk*Gamma*dAdCepr(i,j))/DfDGamma;
}
}
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
DgammaDCepr(i,j) = kk*Gamma*dAdCepr(i,j)+ kk*dDgammaDGamma*DGDCepr(i,j);
}
}
//8. update DdevphiDCepr and DphiEDCepr
static STensor3 dCxdCepr, dPxdCepr;
STensorOperation::zero(dPxdCepr);
STensorOperation::zero(dCxdCepr);
dPxdCepr *= m(1);
dCxdCepr = DgammaDCepr;
dCxdCepr *= m(0);
dCxdCepr += dPxdCepr;
// DphiEDCepr
static STensor3 temp7, temp8, temp9;
STensorOperation::zero(temp7);
STensorOperation::zero(temp8);
STensorOperation::zero(temp9);
temp7 = dCxdCepr;
temp7 *= (1/3*1/(pow(Cx,2))*trXn/v);
temp9 = dCxdCepr;
temp9 *= (2/3*_b*Gamma*pow(kk,1)*Hb/pow(Cx,2)); // pow(kk,2) CHANGE
temp8 = DGDCepr;
temp8 *= (2*_b*(Ke+pow(kk,1)*Hb/(3*Cx))); // pow(kk,2) CHANGE
temp8 -= temp9;
temp8 *= (ptildepr + 1/3*trXn*(1-1/Cx));
temp8 *= 1/pow(v,2);
DphiPDCepr += temp7;
DphiPDCepr -= temp8;
// DdevphiDCepr
static STensor43 temp10, temp11;
static STensor3 temp12, temp13, temp14;
// STensorOperation::zero(temp10);
// STensorOperation::zero(temp11);
STensorOperation::zero(temp12);
STensorOperation::zero(temp13);
STensorOperation::zero(temp14);
temp12 = devXn;
temp12 *= 1/pow(Cx,2);
STensorOperation::prod(temp12, dCxdCepr, 1., temp10);
temp10 *= 1/u;
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
temp13(i,j) = devPhipr(i,j) + (1-1/Cx)*devXn(i,j);
temp14(i,j) = 6*DGDCepr(i,j)*( Ge + pow(kk,1)*Hb/(2*Cx) ) - 6*Gamma*(pow(kk,1)*Hb/(2*pow(Cx,2))*dCxdCepr(i,j)); // pow(kk,2) CHANGE
}
}
STensorOperation::prod(temp13, temp14, 1., temp11);
temp11 *= 1/pow(u,2);
DdevphiDCepr += temp10;
DdevphiDCepr -= temp11;
//9. get DtrNDCepr DdevNdCepr
DtrNDCepr = DphiPDCepr;
DtrNDCepr *= (2*_b);
DdevNDCepr = DdevphiDCepr;
DdevNDCepr *= 3;
// 10. get dFpDCepr and dCeinvdCepr
// dFpDCepr
static STensor43 temp15;
static STensor3 temp15_tr;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
temp15_tr(i,j) = Gamma/3.*DtrNDCepr(i,j);
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
temp15(i,j,k,l) = N(i,j)*DGDCepr(k,l)+ Gamma*DdevNDCepr(i,j,k,l)+ Gamma/3.*_I(i,j)*DtrNDCepr(k,l);
}
static STensor43 EprFp0;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
EprFp0(i,j,k,l) = 0.; for (int s=0; s<3; s++){
EprFp0(i,j,k,l) += dexpAdA(i,s,k,l)*Fp0(s,j);
}
}
STensorOperation::multSTensor43(EprFp0,temp15,dFpDCepr);
// dCeinvdCepr; Here, A = GN
static STensor43 dexpAdCepr, dexpAinvdCepr;
STensorOperation::multSTensor43(dexpAdA,temp15,dexpAdCepr);
static STensor3 expGNinv;
STensorOperation::inverseSTensor3(expGN,expGNinv);
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dexpAinvdCepr(i,s,k,l) = 0.;
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
dexpAinvdCepr(i,s,k,l) -= expGNinv(i,m)*dexpAdCepr(m,j,k,l)*expGNinv(j,s);
}
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dCedCepr(i,s,k,l) = 0.;
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
dCedCepr(i,s,k,l) += expGNinv(m,i)*_I4(m,j,k,l)*expGNinv(j,s) - expGNinv(m,i)*Cepr(m,j)*dexpAinvdCepr(j,s,k,l)
- dexpAinvdCepr(m,i,k,l)*Cepr(m,j)*expGNinv(j,s); // assuming GN is symmetric
} // This derivative is wrong somehow. Use d()/dF below instead.
// 10.1 get dXdCepr
const STensor3 backSig = q1->_backsig;
static double trX;
static STensor3 devX;
STensorOperation::decomposeDevTr(backSig,devX,trX);
STensorOperation::zero(dXdCepr);
STensorOperation::zero(dTrXdCepr);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dTrXdCepr(i,j) = pow(kk,1)*Hb*temp15_tr(i,j)/Cx - trX*dCxdCepr(i,j)/Cx; // pow(kk,2) CHANGE
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
dXdCepr(i,j,k,l) = pow(kk,1)*Hb*temp15(i,j,k,l)/Cx - q1->_backsig(i,j)*dCxdCepr(k,l)/Cx; // CHECK, MIGHT HAVE TO DIFFERENTIATE Hb wrt gamma
// pow(kk,2) CHANGE
}
// 11. update DpDCepr, DdevKDCepr
// (DpKeprDCepr DdevKeprDCepr)
// update
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
DpDCepr(i,j) -= Ke*(DGDCepr(i,j)*trN+Gamma*DtrNDCepr(i,j));
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
DdevKDCepr(i,j,k,l) -= 2.*Ge*(DGDCepr(k,l)*devN(i,j)+Gamma*DdevNDCepr(i,j,k,l));
}
}
}
}
}
else{ // elastic
STensorOperation::zero(DgammaDCepr);
STensorOperation::zero(dFpDCepr);
STensorOperation::zero(DtrNDCepr);
STensorOperation::zero(DdevNDCepr);
STensorOperation::zero(DGDCepr);
STensorOperation::zero(dXdCepr);
STensorOperation::zero(dTrXdCepr);
// STensorOperation::zero(dCedCepr);
dCedCepr = _I4; // DOUBT _ CHECK
}
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dCeinvdCepr(i,s,k,l) = 0.;
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
dCeinvdCepr(i,s,k,l) -= Ceinv(i,m)*dCedCepr(m,j,k,l)*Ceinv(j,s);
}
// 12. get dKcorDcepr
static STensor43 dKcorDcepr;
dKcorDcepr = DdevKDCepr;
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
dKcorDcepr(i,j,k,l) += _I(i,j)*DpDCepr(k,l);
}
}
}
}
//13. get CeprToF conversion tensor
static STensor43 CeprToF;
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
CeprToF(i,j,k,l) = 2.*Fepr(k,i)*invFp0(j,l);
}
}
}
}
static STensor43 DKcorDF;
STensorOperation::multSTensor43(dKcorDcepr,CeprToF,DKcorDF);
// 14. get dSdCepr and dSdF; also get dMedCepr and dMedF
// Done below -> need dCedF
STensor43& dXdF = q1->getRefToDbackStressdF();
static STensor3 dTrXdF;
STensorOperation::multSTensor43(dXdCepr,CeprToF,dXdF);
// DphiPDF - use this for the pressure term in R
// for the pressure term in R
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
DphiPDF(i,j) = 0.;
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DphiPDF(i,j) += DphiPDCepr(k,l)*CeprToF(k,l,i,j);
}
}
// 15. get DgammaDF
STensor3& DgammaDF = q1->_DgammaDF;
STensor3& DGammaDF = q1->_DGammaDF;
if (Gamma > 0){
STensorOperation::multSTensor3STensor43(DgammaDCepr,CeprToF,DgammaDF);
STensorOperation::multSTensor3STensor43(DGDCepr,CeprToF,DGammaDF);
STensorOperation::multSTensor43(dFpDCepr,CeprToF,dFpdF);
}
else{
STensorOperation::zero(DgammaDF);
STensorOperation::zero(DGammaDF);
// STensorOperation::zero(dFpdF);
}
// 16. Everything else and Tangent
static STensor43 DinvFpDF, dCedF, dCeinvDF; //
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DinvFpDF(i,s,k,l) = 0.;
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
DinvFpDF(i,s,k,l) -= Fpinv(i,m)*dFpdF(m,j,k,l)*Fpinv(j,s);
}
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dFedF(m,j,k,l) = _I(m,k)*Fpinv(l,j);
for (int s=0; s<3; s++)
dFedF(m,j,k,l) += F(m,s)*DinvFpDF(s,j,k,l);
}
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dCedF(i,j,k,l) = 0.;
for (int p=0; p<3; p++){
dCedF(i,j,k,l) += Fe(p,i)*dFedF(p,j,k,l) + dFedF(p,i,k,l)*Fe(p,j); // This Works!
}
}
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dCeinvDF(i,s,k,l) = 0.;
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
dCeinvDF(i,s,k,l) -= Ceinv(i,m)*dCedF(m,j,k,l)*Ceinv(j,s);
}
// 17. Tangent - dPdF
static STensor43 dSdCepr, dSdF, dMedCepr;
STensor43& dMedF = q1->getRefToDModMandelDF();
/* Something wrong with dCedCepr - good luck debugging
* for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dMedCepr(i,j,k,l) = dKcorDcepr(i,j,k,l); // CHANGED
for (int m=0; m<3; m++)
for (int s=0; s<3; s++)
dMedCepr(i,j,k,l) += dCedCepr(i,m,k,l)*KS(m,s)*Ceinv(s,j) + Ce(i,m)*(dKcorDcepr(m,s,k,l)*Ceinv(s,j) + KS(m,s)*dCeinvdCepr(s,j,k,l));
// CHECK!! the last terms
} dMedCepr *= 0.5;
STensorOperation::multSTensor43(dMedCepr,CeprToF,dMedF);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dSdCepr(i,j,k,l) = 0.;
for (int m=0; m<3; m++){
// dSdCepr(i,j,k,l) += dKcorDcepr(i,m,k,l)*Ceinv(m,j) + KS(i,m)*dCeinvdCepr(m,j,k,l);
// dSdCepr(i,j,k,l) += dCeinvdCepr(i,m,k,l)*KS(m,j) + Ceinv(i,m)*dKcorDcepr(m,j,k,l);
dSdCepr(i,j,k,l) += dCeinvdCepr(i,m,k,l)*MS(m,j) + Ceinv(i,m)*dMedCepr(m,j,k,l);
// KS is corKir
}}
// dSdCepr *= 0.5;
STensorOperation::multSTensor43(dSdCepr,CeprToF,dSdF);*/
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dMedF(i,j,k,l) = DKcorDF(i,j,k,l); // CHANGED
for (int m=0; m<3; m++)
for (int s=0; s<3; s++)
dMedF(i,j,k,l) += dCedF(i,m,k,l)*KS(m,s)*Ceinv(s,j) + Ce(i,m)*(DKcorDF(m,s,k,l)*Ceinv(s,j) + KS(m,s)*dCeinvDF(s,j,k,l));
// CHECK!! the last terms
}
dMedF *= 0.5;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dSdF(i,j,k,l) = 0.;
for (int m=0; m<3; m++){
// dSdF(i,j,k,l) += DKcorDF(i,m,k,l)*Ceinv(m,j) + KS(i,m)*dCeinvDF(m,j,k,l);
// dSdF(i,j,k,l) += dCeinvDF(i,m,k,l)*KS(m,j) + Ceinv(i,m)*DKcorDF(m,j,k,l);
dSdF(i,j,k,l) += dCeinvDF(i,m,k,l)*MS(m,j) + Ceinv(i,m)*dMedF(m,j,k,l);
// KS is corKir
}}
// dSdF *= 0.5;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
Tangent(i,j,k,l) = 0.;
for (int m=0; m<3; m++){
for (int n=0; n<3; n++){
Tangent(i,j,k,l) += dFedF(i,m,k,l)*S(m,n)*Fpinv(j,n);
Tangent(i,j,k,l) += Fe(i,m)*dSdF(m,n,k,l)*Fpinv(j,n);
Tangent(i,j,k,l) += Fe(i,m)*S(m,n)*DinvFpDF(j,n,k,l);
}
}
}
// dP/dT
// 1. get dMeprDT from dKeprDT -> DcorKirprDT from TVE
static STensor3 dMeprDT, dDevMeprDT, DcorKirDT;
double dpMeprDT(0.);
const STensor3& dKeprDT = q1->getConstRefToDcorKirDT();
STensorOperation::zero(dMeprDT);
static STensor3 temp16, temp17;
STensorOperation::multSTensor3(dKeprDT, Ceinvpr, temp16);
STensorOperation::multSTensor3(Cepr, temp16, temp17);
dMeprDT = 0.5*(dKeprDT + temp17);
STensorOperation::decomposeDevTr(dMeprDT,dDevMeprDT,dpMeprDT);
// dpMeprDT *= 1/3; WHAT IS THE PROBLEM HERE? WHY DOES IT GIVE ZERO?
dpMeprDT = dMeprDT.trace()/3;
DcorKirDT = dKeprDT; // update later
// 2. get dPhipprDT and dDevPhiprDT
static double dPhipprDT;
static STensor3 dDevPhiprDT;
dPhipprDT = dpMeprDT;
dDevPhiprDT = dDevMeprDT;
// 3. get dCxdT, dGxdT, dKxdT, dXdT
static double dCxdT, dGxdT, dKxdT;
dCxdT = 1/(pow(Hb,2))*pow(dHbdT,2)*(T-T0) - 1/Hb*ddHbddT*(T-T0) - 1/Hb*dHbdT;
dGxdT = DGDTsum + pow(kk,1)/(2*Cx)*dHbdT - (pow(kk,1)*Hb)/(2*pow(Cx,2))*dCxdT; // pow(kk,2) CHANGE
dKxdT = DKDTsum + pow(kk,1)/(3*Cx)*dHbdT - (pow(kk,1)*Hb)/(3*pow(Cx,2))*dCxdT; // pow(kk,2) CHANGE
// need this for DmechsourceDT -> initialise here (X is backStress)
STensor3& dXdT = q1->getRefToDbackStressdT();
STensorOperation::zero(dXdT);
dXdT -= (q1->_backsig)*(dCxdT);
dXdT *= 1/Cx;
// 4. get dPhiPdT and dDevPhiDT
static double dPhiPdT;
static STensor3 dDevPhiDT;
STensorOperation::zero(dDevPhiDT);
static double temp18;
static STensor3 temp19;
dPhiPdT = dPhipprDT/v + 1/3*trXn/pow(Cx,2)*dCxdT/v; // no correction to this
temp18 = (ptildepr+1/3*trXn*(1-1/Cx))*(2*_b*Gamma*dKxdT)/pow(v,2); // Gamma derivatives will change this
dPhiPdT -= temp18;
temp19 = (devPhipr + devXn*(1-1/Cx))*(6*Gamma*dGxdT);
temp19 *= 1/pow(u,2);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
dDevPhiDT(i,j) = (dDevPhiprDT(i,j) + devXn(i,j)/pow(Cx,2)*dCxdT)/u - temp19(i,j);
}
}
// 5. get dPhiEdT
static double dPhiEdT;
STensorOperation::doubleContractionSTensor3(DphiEDdevPhi,dDevPhiDT,dPhiEdT);
// 6. to 11. (inside if loop-> Gamma > 0.)
static double dAdT, dfdT;
static double dGammaDT;
static double DgammaDT;
static double DtrNdT;
static STensor3 DdevNdT;
double eta(0.),Deta(0.),DetaDT(0.);
if (_viscosity != NULL)
_viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT);
double etaOverDt = eta/this->getTimeStep();
if (Gamma >0){
// plastic
// 6. get DaDT, dAdT and dfdT
fullVector<double> DaDT(3), DDaDTT(3);
getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
double a1, a2, a3;
double a4, a5, a6;
a1 = a(0); a2 = a(1); a3 = a(2);
a4 = DaDT(0); a5 = DaDT(1); a6 = DaDT(2);
dAdT = (4.*_b*_b*ptilde/(A*3.))*dPhiPdT + (6./A)*PhiEq*dPhiEdT;
dfdT = DaDT(2)*pow(PhiEq,_n) + a(2)*_n*pow(PhiEq,_n-1.)*dPhiEdT - DaDT(1)*ptilde - a(1)*dPhiPdT - DaDT(0);
if (this->getTimeStep()>0){
dfdT -= (DetaDT*Gamma/this->getTimeStep())*_p*pow((eta*Gamma/this->getTimeStep()),(_p-1));
}
// 7. get dGammaDT, DgammaDT
dGammaDT = (-dfdT-dfdDgamma*kk*Gamma*dAdT)/DfDGamma;
DgammaDT = kk*Gamma*dAdT + kk*dDgammaDGamma*dGammaDT;
// 8. update DdevphiDT and DphiPDT
static double dPxdT, temp18_2;
static STensor3 temp19_2;
dPxdT = 0.*m(1);
dCxdT += m(0)*DgammaDT;
dCxdT += dPxdT;
dPhiPdT += temp18;
temp18_2 = (ptildepr+1/3*trXn*(1-1/Cx))*(2*_b*(dGammaDT*(Ke+pow(kk,1)*Hb/(3*Cx)) + Gamma*dKxdT))/pow(v,2); // pow(kk,2) CHANGE
dPhiPdT -= temp18_2;
dDevPhiDT += temp19;
temp19_2 = (devPhipr+devXn*(1-1/Cx))*(6*(dGammaDT*(Ge+pow(kk,1)*Hb/(2*Cx)) + Gamma*dGxdT)); // pow(kk,2) CHANGE
temp19_2 *= 1/pow(u,2);
dDevPhiDT -= temp19_2;
// 9. get DtrNdT DdevNdT
DtrNdT = dPhiPdT;
DtrNdT *= (2*_b);
DdevNdT = dDevPhiDT;
DdevNdT *= 3;
// 10. get dFpdT
static STensor3 temp20;
STensorOperation::zero(temp20);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
temp20(i,j) += N(i,j)*dGammaDT + Gamma*DdevNdT(i,j) + 1/3*Gamma*_I(i,j)*DtrNdT;
static STensor43 EprFp0;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
EprFp0(i,j,k,l) = 0.;
for (int s=0; s<3; s++){
EprFp0(i,j,k,l) += dexpAdA(i,s,k,l)*Fp0(s,j);
}
}
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
dFpdT(i,j) = 0.;
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
dFpdT(i,j) += EprFp0(i,j,k,l)*temp20(k,l);
}
}
}
}
// 10.1 dXdT - backstress temperature derivative -> correction
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
dXdT(i,j) += 1/Cx*((pow(kk,1)*(dHbdT*Gamma*N(i,j) + Hb*temp20(i,j)))); // pow(kk,2) CHANGE
}
}
// 11. DcorKirDT
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
DcorKirDT(i,j) -= ( DKDTsum*Gamma*trN + Ke*(dGammaDT*trN+Gamma*DtrNdT))*_I(i,j);
DcorKirDT(i,j) -= 2.*( DGDTsum*Gamma*devN(i,j) + Ge*(dGammaDT*devN(i,j)+Gamma*DdevNdT(i,j)));
}
}
}
else{
// elastic
STensorOperation::zero(DgammaDT);
// STensorOperation::zero(dFpdT);
STensorOperation::zero(DtrNdT);
STensorOperation::zero(DdevNdT);
STensorOperation::zero(dGammaDT);
}
// 11.1
q1->_DGammaDT = dGammaDT;
q1->_DgammaDT = DgammaDT;
// 12. get dKcorDT
// done above
// 13. get DinvFpdT, dFedT, dCedT, dCeinvDT
static STensor3 DinvFpdT, dFedT, dCedT, dCeinvdT;
STensorOperation::zero(DinvFpdT);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int p=0; p<3; p++)
for (int q=0; q<3; q++){
DinvFpdT(i,j) -= Fpinv(i,p)*dFpdT(p,q)*Fpinv(q,j);
}
STensorOperation::zero(dFedT);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++){
dFedT(i,j) += F(i,k)*DinvFpdT(k,j);
}
STensorOperation::zero(dCedT);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
for (int p=0; p<3; p++){
dCedT(i,j) += Fe(p,i)*dFedT(p,j) + dFedT(p,i)*Fe(p,j); // This Works!
}
}
STensorOperation::zero(dCeinvdT);
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dCeinvdT(i,s) -= Ceinv(i,k)*dCedT(k,l)*Ceinv(l,s);
}
// 13.1 get dMeDT -> needed for DmechSourceDT
STensor3& dMeDT = q1->getRefToDModMandelDT();
// STensorOperation::zero(dMeDT); for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dMeDT(i,j) = DcorKirDT(i,j);
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dMeDT(i,j) += dCedT(i,k)*KS(k,l)*Ceinv(l,j) + Ce(i,k)*(DcorKirDT(k,l)*Ceinv(l,j) + KS(k,l)*dCeinvdT(l,j));
}
}
dMeDT *= 0.5;
// 14. get dSdT
static STensor3 dSdT;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dSdT(i,j) = 0.;
for (int m=0; m<3; m++){
// dSdT(i,j) += DcorKirDT(i,m)*Ceinv(m,j) + KS(i,m)*dCeinvdT(m,j) + dCeinvdT(i,m)*KS(m,j) + Ceinv(i,m)*DcorKirDT(m,j);
dSdT(i,j) += Ceinv(i,m)*dMeDT(m,j) + dCeinvdT(i,m)*MS(m,j);
}
}
// dSdT *= 0.5;
// 15. Finally, get dPdT
// static STensor3 dPdT;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dPdT(i,j) = 0.;
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dPdT(i,j) += (dFedT(i,k)*S(k,l)*Fpinv(j,l) + Fe(i,k)*dSdT(k,l)*Fpinv(j,l) + Fe(i,k)*S(k,l)*DinvFpdT(j,l));
}
}
}
if(stiff){
// TVP - flux, thermalEnergy, thermalSource
if (!_thermalEstimationPreviousConfig){
static STensor3 DJDF;
static STensor43 DCinvDF;
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
DJDF(i,j) = J*Finv(j,i);
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
DCinvDF(i,j,k,l) = 0.;
for (int p=0; p<3; p++){
for (int a=0; a<3; a++){
for (int b=0; b<3; b++){
DCinvDF(i,j,k,l) -= 2.*F(k,p)*Cinv(i,a)*Cinv(j,b)*_I4(a,b,p,l);
}
}
}
}
}
}
}
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
dfluxTdF(i,j,k) += (DJDF(j,k)*fluxT(i)/J);
for (int l=0; l<3; l++){
dfluxTdF(i,j,k) -= (J*DCinvDF(i,l,j,k)*gradT(l)*KThConT);
}
}
}
}
} }
// thermSrc and MechSrc
// thermalEnergy
double CpT;
if (stiff){
double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
// CpT = -T*DDpsiTVMdTT;
getCp(CpT,T);
}
else{
getCp(CpT,T);
}
q1->_thermalEnergy = CpT*T;
// thermalSource
if (this->getTimeStep() > 0.){
thermalSource = -CpT*(T-T0)/this->getTimeStep();
}
else
thermalSource = 0.;
// mechSourceTVP
double& Wm = q1->getRefToMechSrcTVP();
double& Wm_TVE = q1->getRefToMechSrcTVE();
mechanicalSource = 0.;
mlawNonLinearTVE::getMechSourceTVE(q0,q1,T0,T,DKDTsum,DGDTsum,_I4,&Wm_TVE); // mechSourceTVE
mechanicalSource += Wm_TVE;
getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,&Wm);
mechanicalSource += Wm;
// freeEnergy and elastic energy
double& psiTVM = q1->getRefToFreeEnergyTVM();
getFreeEnergyTVM(q0,q1,T0,T,&psiTVM,NULL);
q1->_elasticEnergy = mlawNonLinearTVE::freeEnergyMechanical(*q1,T); // deformationEnergy(*q1,T);
if (stiff){
const double& DDpsiTVMdTT_0 = q0->getConstRefToDDFreeEnergyTVMdTT();
double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
// double CpT = -T*DDpsiTVMdTT;
// double CpT_0 = -T0*DDpsiTVMdTT_0;
// thermal source derivatives
double DCpDT(0.);
mlawNonLinearTVE::getCp(CpT,T,&DCpDT);
static STensor3 DCpDF;
STensorOperation::zero(DCpDF); // CHANGE for DCpDF
if (this->getTimeStep() > 0){
dthermalSourcedT = -DCpDT*(T-T0)/this->getTimeStep() - CpT/this->getTimeStep();
// dthermalSourcedT = -CpT/this->getTimeStep() - (CpT-CpT_0)/this->getTimeStep(); // thermalSource = -CpT*(T-T0)/this->getTimeStep();
for(int i = 0; i< 3; i++){
for(int j = 0; j< 3; j++){
dthermalSourcedF(i,j) = -DCpDF(i,j)*(T-T0)/this->getTimeStep();
}
}
}
else{
dthermalSourcedT = 0.;
STensorOperation::zero(dthermalSourcedF);
}
// mechSourceTVP derivatives
// conversion tensor for mechSourceTVE Terms
static STensor43 DCeDFe, DEeDFe;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DCeDFe(i,j,k,l) = 0.;
for (int p=0; p<3; p++){
DCeDFe(i,j,k,l) += 2*F(p,i)*dFedF(p,j,k,l);
}
}
STensorOperation::multSTensor43(DlnDCe,DCeDFe,DEeDFe);
DEeDFe *= 0.5;
double& dWmdT_TVE = q1->getRefTodMechSrcTVEdT();
STensor3& dWmdF_TVE = q1->getRefTodMechSrcTVEdF();
mlawNonLinearTVE::getMechSourceTVE(q0,q1,T0,T,DKDTsum,DGDTsum,DEeDFe,NULL,&dWmdF_TVE,&dWmdT_TVE);
dmechanicalSourcedT = 0.;
STensorOperation::zero(dmechanicalSourceF);
dmechanicalSourceF += dWmdF_TVE;
dmechanicalSourcedT += dWmdT_TVE;
double& dWmdT = q1->getRefTodMechSrcTVPdT();
STensor3& dWmdF = q1->getRefTodMechSrcTVPdF();
getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,NULL,&dWmdF,&dWmdT);
dmechanicalSourceF += dWmdF;
dmechanicalSourcedT += dWmdT;
}
};
void mlawNonLinearTVP::predictorCorrector_TVP_AssociatedFlow(const STensor3& F, const IPNonLinearTVP *q0, IPNonLinearTVP *q1,
STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, const double T0, const double T) const{
// compute elastic predictor
STensor3& Fp1 = q1->_Fp;
const STensor3& Fp0 = q0->_Fp;
Fp1 = Fp0; // plastic deformation tensor
q1->_epspbarre = q0->_epspbarre; // plastic equivalent strain
q1->_epspCompression = q0->_epspCompression;
q1->_epspTraction = q0->_epspTraction;
q1->_epspShear = q0->_epspShear;
q1->_backsig = q0->_backsig; // backstress tensor
q1->_DgammaDt = 0.; // plastic rate --> failure
static STensor3 Fpinv, Ce, Fepr;
STensorOperation::inverseSTensor3(Fp1,Fpinv);
STensorOperation::multSTensor3(F,Fpinv,Fepr);
STensorOperation::multSTensor3FirstTranspose(Fepr,Fepr,Ce);
static STensor3 invFp0; // plastic predictor
invFp0= Fpinv;
STensor3& Fe = q1->_Fe;
Fe = Fepr;
static STensor43 DlnDCepr, DlnDCe;
static STensor63 DDlnDDCe;
static STensor43 dexpAdA; // estimation of dexpA/dA
STensor3& Ee = q1->_Ee;
STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCepr,&DDlnDDCe);
Ee *= 0.5;
DlnDCe = DlnDCepr;
// update A, B
// double Ge, Ke;
double Ke, Ge = 0.; double Kde, Gde = 0.; double KTsum, GTsum = 0.;
double DKe, DGe = 0.; double DKde, DGde = 0.; double DKDTsum, DGDTsum = 0.;
ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T);
static STensor3 devKpr; // dev corotational kirchoff stress predictor
static double ppr; // pressure predictor
STensorOperation::decomposeDevTr(q1->_kirchhoff,devKpr,ppr);
ppr /= 3.;
double keqpr = sqrt(1.5*devKpr.dotprod());
static STensor3 devK;
devK= devKpr; // dev corotational kirchoff stress
double p = ppr; // pressure
// hardening
this->hardening(q0,q1,T);
fullVector<double> a(3), Da(3); // yield coefficients and derivatives in respect to plastic deformation
this->getYieldCoefficients(q1,a);
static STensor3 devN; // dev part of yield normal
double trN = 0.; // trace part of yield normal
static STensor3 N; // yield normal
STensorOperation::zero(devN);
STensorOperation::zero(N);
double sigVM = keqpr;
double Dgamma = 0.;
double Gamma = 0.;
double g0 = 0.;
double z = a(1)/a(2)*pow(sigVM,1.-_n);
double A = sqrt(1.5+1.*z*z/(3.*_n*_n));
static fullMatrix<double> invJ(2,2); // inversed jacobian
invJ.setAll(0.);
double dAdz =0.;
double dzdDgamma = 0.;
double dzdsigVM = 0.;
double dg0dsigVM = -1./(3.*Ge);
if (q1->dissipationIsBlocked()){
q1->getRefToDissipationActive() = false;
}
else{
double f = a(2)*pow(sigVM,_n) - (a(1)*ppr+a(0));
if (f>_tol){
q1->getRefToDissipationActive() = true;
// plasticity
int ite = 0;
int maxite = 1000; // maximal number of iters
double val = 1.5*a(2)*_n*pow(3.*fabs(p),_n-2.);
q1->_nup = (val*p+ a(1)/3.)/(val*2.*p - a(1)/3.);
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
A*=kk;
fullVector<double> R(2); // residual of eqs to estimate Dgamma, Dpression and Gamma
R(0) = Dgamma - g0*A;
R(1) = a(2)*pow(sigVM,_n) - a(1)*(ppr+Ke*z*g0/_n) - a(0);
double res = R.norm();
while(res>_tol or ite<1){
this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
static fullMatrix<double> J(2,2);
J.setAll(0.);
dAdz = kk*z/(3.*_n*_n*A);
dzdDgamma = (Da(1)*a(2) - a(1)*Da(2))/(a(2)*a(2))*pow(sigVM,1.-_n);
dzdsigVM = a(1)/a(2)*(1.-_n)*pow(sigVM,-_n);
J(0,0) = 1. - g0*dAdz*dzdDgamma;
J(0,1) = -dg0dsigVM*A-g0*dAdz*dzdsigVM;
J(1,0) = Da(2)*pow(sigVM,_n) - Da(1)*(ppr+Ke*z*g0/_n) - Da(0) - a(1)*Ke*dzdDgamma*g0/_n;
J(1,1) = a(2)*_n*pow(sigVM,_n-1) - a(1)*Ke/_n*(dzdsigVM*g0+z*dg0dsigVM);
double detJ = J(0,0)*J(1,1) - J(1,0)*J(0,1);
if (detJ == 0.) Msg::Error("the corrected system can not be solved: mlawNonLinearTVP::predictorCorrector_associatedFlow");
invJ(0,0) = J(1,1)/detJ;
invJ(0,1) = -J(0,1)/detJ;
invJ(1,0) = -J(1,0)/detJ;
invJ(1,1) = J(0,0)/detJ;
double solDgamma = -invJ(0,0)*R(0)-invJ(0,1)*R(1);
double solsigVM = -invJ(1,0)*R(0)-invJ(1,1)*R(1);
// bool isSolved = J.luSolve(residual,sol);
// if (!isSolved) Msg::Error("the corrected system can not be solved: mlawPowerYieldHyper::predictorCorrector");
//sol.print("sol");
if (solDgamma >0.1)
Dgamma+= 0.1;
else if (Dgamma+ solDgamma< 0.)
Dgamma /= 2.;
else
Dgamma += solDgamma;
if (sigVM + solsigVM< 0.)
sigVM /= 2.;
else
sigVM += solsigVM;
updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
this->hardening(q0,q1,T);
this->getYieldCoefficients(q1,a);
g0 = (keqpr- sigVM)/(3.*Ge);
z = a(1)/a(2)*pow(sigVM,1.-_n);
A = kk*sqrt(1.5+1.*z*z/(3.*_n*_n));
R(0) = Dgamma - g0*A;
R(1) = a(2)*pow(sigVM,_n) - a(1)*(ppr+Ke*z*g0/_n) - a(0);
res = R.norm();
ite++;
//Msg::Info("iter = %d, res = %e, Dgamma = %e, sigVM/keqpr = %e",ite,res, Dgamma, sigVM/keqpr);
if(ite > maxite){
Msg::Error("No convergence for plastic correction in mlawPowerYieldHyper iter = %d, res = %e!!",ite,res);
P(0,0) = P(1,1) = P(2,2) = sqrt(-1.);
break;
}
}
Gamma = z*g0/(_n*a(1));
// update
p += Ke*Gamma*a(1);
double ff =a(2)*_n*pow(sigVM,_n-2.);
devK*= 1./(1.+3.*Ge*Gamma*ff);
// estimate yield normal
devN = devK;
devN *= 1.5*ff;
trN = -a(1);
N = devN;
N(0,0) += trN/3.;
N(1,1) += trN/3.;
N(2,2) += trN/3.;
// estimate exp(GammaN)
static STensor3 expGN;
static STensor3 GammaN;
GammaN = N;
GammaN *= Gamma;
STensorOperation::expSTensor3(GammaN,_order,expGN,&dexpAdA);
// update plastic deformation tensor
STensorOperation::multSTensor3(expGN,Fp0,Fp1);
// update IP
updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
// update elastic deformation tensor, corotational stress
STensorOperation::inverseSTensor3(Fp1,Fpinv);
STensorOperation::multSTensor3(F,Fpinv,Fe);
STensorOperation::multSTensor3FirstTranspose(Fe,Fe,Ce);
STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCe,&DDlnDDCe);
Ee *= 0.5;
//
ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T,false);
}
else{
q1->getRefToDissipationActive() = false;
}
}
// corotational Kirchhoff stress tenor
STensor3& KS = q1->_kirchhoff;
KS = devK;
KS(0,0) += p;
KS(1,1) += p;
KS(2,2) += p;
// first Piola Kirchhoff stress
static STensor3 S;
S*=(0.);
for(int i=0; i<3; i++)
for(int j=0; j<3; j++)
for(int k=0; k<3; k++)
for(int l=0; l<3; l++)
S(i,j)+= KS(k,l)*DlnDCe(k,l,i,j);
STensorOperation::zero(P);
for(int i=0; i<3; i++)
for(int j=0; j<3; j++)
for(int k=0; k<3; k++)
for(int l=0; l<3; l++)
P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
// defo energy
q1->_elasticEnergy=deformationEnergy(*q1);
q1->getRefToPlasticEnergy() = q0->plasticEnergy();
if (Gamma > 0){
double dotKSN = dot(KS,N);
q1->getRefToPlasticEnergy() += Gamma*dotKSN;
}
if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
q1->getRefToIrreversibleEnergy() = q1->defoEnergy();
}
else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
(this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
q1->getRefToIrreversibleEnergy() = q1->plasticEnergy();
}
else{
q1->getRefToIrreversibleEnergy() = 0.;
}
if (stiff){
// FILL THIS
}
};