diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index 3206b2ac6331291a40ab66f8d8ca75bcd8328af7..44ae251b26e11b5160339404eace358e6a74b83e 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -152,7 +152,7 @@ void mlawNonLinearTVP::getChabocheCoeffs(fullVector<double>& coeffs, const doubl
coeffs.resize(2);
- coeffs(0) = 0; // 1; //DEBUGGING
+ coeffs(0) = 1; //DEBUGGING
coeffs(1) = 0;
// add dependency on p later
@@ -1125,7 +1125,6 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static fullVector<double> m(2);
double Px(0.);
getChabocheCoeffs(m,0,q1);
- getPlasticPotential(ptildepr, PhiEqpr, Px); // CHANGE -> this is wrong, change to backstress instead of phi
double DfDGamma = 0.;
double dfdDgamma = 0.;
@@ -1133,28 +1132,19 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
double v = 1.;
// Initialise Cx, Gt, Kt
- double Cx = m(0)*Dgamma + m(1)*Px - 1/Hb * dHbdT * (T-T0) + 1; // CHECK
- double Gt = Ge;
- double Kt = Ke;
+ double Cx = m(0)*Dgamma + m(1)*Px + 1;
+ if (Hb>0.){Cx -= 1/Hb * dHbdT * (T-T0);}
+ 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
// Initialise ptilde and devPhi (phi)
- ptilde = ptildepr + 1/3*trXn; // current effective pressure -> first invariant
- ptilde *= 1/v;
- devPhi = devPhipr + devXn;
- devPhi *= 1/u;
-
- if (abs(Cx) > 0.0){
- Gt += pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
- Kt += pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
-
- ptilde -= 1/3*trXn*1/Cx;
- devPhi -= devXn;
- devPhi *= 1/Cx;
- }
+ ptilde = ptildepr; // current effective pressure -> first invariant
+ devPhi = devPhipr;
double PhiEqpr2 = 1.5*devPhi.dotprod();
double PhiEqpr = sqrt(PhiEqpr2); // -> second invariant
double PhiEq = PhiEqpr; // current effective deviator
+ getPlasticPotential(ptildepr, PhiEqpr, Px); // CHANGE -> this is wrong, change to backstress instead of phi
// Get f and A
double f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
@@ -1168,12 +1158,11 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
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){
@@ -1199,18 +1188,17 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static STensor3 DdevPhidDgamma, term2;
STensorOperation::zero(DdevPhidDgamma);
- if (abs(Cx) > 0.0){
- // 1st term of He
- DdevPhidDgamma = 3*Gamma*pow(kk,1)*Hb*devPhipr; // pow(kk,2) CHANGE
- DdevPhidDgamma -= devXn*(1 + 6*Gamma*(Ge + pow(kk,1)*Hb/2)); // pow(kk,2) CHANGE
- DdevPhidDgamma *= 1/(Cx*Cx*u*u) * dCxdDgamma;
+ // 1st term of He
+ DdevPhidDgamma = 3*Gamma*pow(kk,1)*Hb*devPhipr; // pow(kk,2) CHANGE
+ DdevPhidDgamma -= devXn*(1 + 6*Gamma*(Ge + pow(kk,1)*Hb/2)); // pow(kk,2) CHANGE
+ DdevPhidDgamma *= 1/(Cx*Cx*u*u) * dCxdDgamma;
- // 2nd term of He
- term2 = 3*Gamma*pow(kk,1)*Cx*devPhi; // pow(kk,2) CHANGE
- term2 *= 1/(Cx*Cx*u) * dHb;
+ // 2nd term of He
+ term2 = 3*Gamma*pow(kk,1)*Cx*devPhi; // pow(kk,2) CHANGE
+ term2 *= 1/(Cx*Cx*u) * dHb;
- DdevPhidDgamma -= term2;
- }
+ DdevPhidDgamma -= term2;
+
double Stemp = STensorOperation::doubledot(devPhi,DdevPhidDgamma);
double He = 1.5*Stemp/PhiEq;
@@ -1218,59 +1206,79 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
Hp -= (2*_b*Gamma*pow(kk,1)*Cx*ptilde)*1/(3*Cx*Cx*v) * dHb;
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.);
+ DfDGamma -= pow(etaOverDt,_p)*_p*pow(Gamma,(_p-1.));
- double dGamma = -f/DfDGamma; // NR
+ // double dGamma = -f/DfDGamma; // NR
+
+ double dfdGamma = - (_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 dAdDgamma = (12*PhiEq*He + 8/3*_b*_b*ptilde*Hp)/(2*A);
+ double dgdDgamma = 1 - kk*Gamma*dAdDgamma;
+ double dgdGamma = - dDgammaDGamma;
+ double g = Dgamma - kk*Gamma*A;
+
+ double dGamma = (-f + dfdDgamma*g/dgdDgamma)/(-dfdDgamma*dgdGamma/dgdDgamma + dfdGamma); // NR
+ double dDgamma = -(g+dgdGamma*dGamma)/dgdDgamma; // NR
+
if (Gamma + dGamma <=0.){
Gamma /= 2.; // NR
}
else
Gamma += dGamma;
-
+
+ if (Dgamma + dDgamma <=0.){
+ Dgamma /= 2.;
+ }
+ else
+ Dgamma += dDgamma;
+
//Msg::Error("Gamma = %e",Gamma);
+ //Msg::Error("Dgamma = %e",Dgamma);
+
+ // Update gamma and all Hardening moduli
+ updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
+ hardening(q0,q1,T);
+ getYieldCoefficients(q1,a);
+ if (q1->_ipKinematic != NULL){
+ Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
+ dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter derivative (dHb/dgamma ??)
+ dHbdT = Hb * q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
+ ddHbddT = Hb * q1->_ipKinematic->getDDRDDT();
+ }
+ //a.print("a update");
- // Correct Phi, A, delgamma
- Cx = m(0)*Dgamma + m(1)*Px - 1/Hb * dHbdT * (T-T0) + 1;
+ // Update Cx, u, v
+ Cx = m(0)*Dgamma + m(1)*Px + 1;
+ if (Hb>0.){Cx -= 1/Hb * dHbdT * (T-T0);}
+ Gt = Ge + pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
+ Kt = Ke + pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
u = 1.+6.*Gt*Gamma;
v = 1.+2.*_b*Kt*Gamma;
- ptilde = ptildepr + 1/3*trXn; // current effective pressure -> first invariant
- ptilde *= 1/v;
- devPhi = devPhipr + devXn;
- devPhi *= 1/u;
-
- devPhi = (devPhipr + devXn * (1-1/Cx));
+ // Update Phi
+ devPhi = devPhipr + devXn*(1-1/Cx);
+ ptilde = ptildepr + 1/3*trXn*(1-1/Cx);
devPhi *= 1./u;
PhiEq = sqrt(1.5*devPhi.dotprod());
- ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v;
-
- if (abs(Cx) > 0.0){
- Gt += pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
- Kt += pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
-
- ptilde -= 1/3*trXn*1/Cx;
- devPhi -= devXn;
- devPhi *= 1/Cx;
- }
+ ptilde *= 1./v;
+ // Update A
A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
- Dgamma = kk*Gamma*A;
-
+
+ // 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");
-
+
+ // Update f
f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
double viscoTerm = etaOverDt*Gamma; // I THINK temperature dependency of viscosity should be present in eta, but this relation remains the same
if (Gamma>0 and etaOverDt>0) f-= pow(viscoTerm,_p);
@@ -1290,14 +1298,14 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
q1->_DgammaDt = Dgamma/this->getTimeStep();
- // update effective stress tensor
- // NEW
- devPhi = (devPhipr + devXn * (1-1/Cx));
- devPhi *= 1./u; //NEW
- ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v; //NEW
- // NEW
+ // Correct Phi (effective stress tensor)
+ devPhi = devPhipr + devXn*(1-1/Cx);
+ ptilde = ptildepr + 1/3*trXn*(1-1/Cx);
+ devPhi *= 1./u;
+ PhiEq = sqrt(1.5*devPhi.dotprod());
+ ptilde *= 1./v;
- // update normal
+ // Correct normal
devN = devPhi;
devN *= 3.;
trN = 2.*_b*ptilde;
@@ -1313,45 +1321,33 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
GammaN *= Gamma;
STensorOperation::expSTensor3(GammaN,_order,expGN,&dexpAdA);
- // update plastic deformation tensor
+ // Correct plastic deformation tensor
STensorOperation::multSTensor3(expGN,Fp0,Fp1);
- // update IP
- updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma); // Why again???? WHAT??
+ // Correct IP gamma
+ updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
- // update elastic deformation tensor, corotational stress
+ // Correct elastic deformation tensor (Ee), 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
+ // Correct 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
+ // Correct backstress
static STensor3 DB; // increment
- DB = (N); // increment
- DB *= (kk*Hb*Gamma); // pow(kk,2) CHANGE
- // NEW
+ DB = (N);
+ DB *= (pow(kk,1)*Hb*Gamma); // pow(kk,2) CHANGE
DB += q0->_backsig;
DB *= 1/Cx;
DB -= q0->_backsig;
- // NEW
- q1->_backsig += DB; // update
-
- /*
- // corotationaal Kirchhoff stress
- q1->_kirchhoff = devPhi;
- q1->_kirchhoff += q1->_backsig;
+ q1->_backsig += DB;
- q1->_kirchhoff(0,0) += (ptilde);
- q1->_kirchhoff(1,1) += (ptilde);
- q1->_kirchhoff(2,2) += (ptilde);
-
- */
- // NEW
+ // Correct Mandel
q1->_ModMandel = devPhi;
q1->_ModMandel += q1->_backsig;
@@ -1359,15 +1355,12 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
q1->_ModMandel(1,1) += (ptilde);
q1->_ModMandel(2,2) += (ptilde);
- // NEW
-
}
else{
q1->getRefToDissipationActive() = false;
}
}
-
const STensor3& KS = q1->_kirchhoff;
// second Piola Kirchhoff stress
// static STensor3 S;
@@ -1619,7 +1612,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
temp8 *= (ptildepr + 1/3*trXn*(1-1/Cx));
temp8 *= 1/pow(v,2);
- DphiPDCepr -= temp7; // CHECK sign
+ DphiPDCepr += temp7; // CHECK sign - corrected
DphiPDCepr -= temp8;
// DdevphiDCepr
@@ -1631,8 +1624,8 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
STensorOperation::zero(temp13);
STensorOperation::zero(temp14);
- temp13 = devXn; // This must be temp12?
- temp13 *= 1/pow(Cx,2); // This must be temp12?
+ temp12 = devXn; // This must be temp12?
+ temp12 *= 1/pow(Cx,2); // This must be temp12?
STensorOperation::prod(temp12, dCxdCepr, 1., temp10);
temp10 *= 1/u;
@@ -1647,7 +1640,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
STensorOperation::prod(temp13, temp14, 1., temp11);
temp11 *= 1/pow(u,2);
- DdevphiDCepr -= temp10; // CHECK sign
+ DdevphiDCepr += temp10; // CHECK sign - corrected
DdevphiDCepr -= temp11;
//9. get DtrNDCepr DdevNdCepr
@@ -1992,7 +1985,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static double temp18;
static STensor3 temp19;
- dPhiPdT = dPhipprDT/v - 1/3*trXn/pow(Cx,2)*dCxdT/v; // no correction to this // CHECK the sign of the 2nd term
+ dPhiPdT = dPhipprDT/v + 1/3*trXn/pow(Cx,2)*dCxdT/v; // no correction to this // CHECK the sign of the 2nd term -corrected
temp18 = (ptildepr+1/3*trXn*(1-1/Cx))*(2*_b*Gamma*dKxdT)/pow(v,2); // Gamma derivatives will change this
dPhiPdT -= temp18;
@@ -2000,7 +1993,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
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); // CHECK the sign of the 2nd term
+ dDevPhiDT(i,j) = (dDevPhiprDT(i,j) + devXn(i,j)/pow(Cx,2)*dCxdT)/u - temp19(i,j); // CHECK the sign of the 2nd term -corrected
}
}
@@ -2607,848 +2600,3 @@ void mlawNonLinearTVP::predictorCorrector_TVP_AssociatedFlow(const STensor3& F,
// FILL THIS
}
};
-
-/*
-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& 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;
- 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;
- STensorOperation::multSTensor3(Dp,Fp_dot,Fpinv);
-
- // 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
- getIsotropicHardeningForce(q0,q1,T0,T,Gamma,dFpdT,dFpdF,dPhiPdF,ddRdTT,ddRdFdT);
-
- double eta(0.),Deta(0.),DetaDT(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);
- double etaOverDt = eta/this->getTimeStep();
- double viscoTerm = etaOverDt*Gamma;
-
- const double& R0 = 0;
-
- // double& R = q1->getRefToIsotropicHardeningForce(); //change
- const double& dRdT_0 = 0;
- // double& dRdT = q1->_dIsoHardForcedT;
- double R(0.);
- double dRdT(0.);
- R = -a(1)/a(2) * (trMe-trX1);
- dRdT = (-DaDT(1)/a(2) + a(1)/pow(a(2),2) *DaDT(2))*(trMe-trX1) + (-a(1)/a(2))*(dTrMeDT-dTrXdT);
- R -= a(0)/a(2);
- dRdT += (-DaDT(0)/a(2) + a(0)/pow(a(2),2)*DaDT(2));
- if (Gamma>0 and etaOverDt>0){
- R -= 1/a(2)*pow(eta*Gamma/this->getTimeStep(),_p);
- dRdT += (1/pow(a(2),2)*DaDT(2)*pow(eta*Gamma/this->getTimeStep(),_p)
- - 1/a(2)*_p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(DetaDT*Gamma/this->getTimeStep()
- + eta*dGammaDT/this->getTimeStep()) );
- }
-
- R -= T*dRdT;
-
- // mechSourceTVP
- if(Wm!=NULL){
- double 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;
-
- 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,j,k,s)*Fpinv(s,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++){
- dDpdFp(i,j,k,l) = 0.;
- for (int m=0; m<3; m++)
- dDpdFp(i,j,k,l) += 1/this->getTimeStep()*_I4(i,j,k,m)*Fpinv(m,l) + Fp_dot(i,m)*dFpinvdFp(m,j,k,l);
- }
-
- // 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 k=0; k<3; k++)
- for (int l=0; l<3; l++)
- dDpdT(i,j) += dDpdFp(i,j,k,l)*dFpdT(k,l);
- }
-
- // 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) += 2*dXdT(i,j)/(T-T0) - 2*(X1(i,j)-X0(i,j))/pow((T-T0),2);
- 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))/(T-T0);
- }
- }
- }
-
- // 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;
- }
-
-
- // dRdF
- // static STensor3 dRdF;
- const STensor3& dRdF_0 = 0.; //q0->_dIsoHardForcedF;
- // STensor3& dRdF = q1->_dIsoHardForcedF;
- static STensor3 dRdF;
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- dRdF(i,j) = (-Da(1)/a(2) + a(1)/pow(a(2),2) *Da(2))*(trMe-trX1)*DgammaDF(i,j) + (-a(1)/a(2))*dPhiPdF(i,j) + (-Da(0)/a(2) + a(0)/pow(a(2),2)*Da(2))*DgammaDF(i,j);
- }
-
- if (Gamma>0 and etaOverDt>0){
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- dRdF += (1/pow(a(2),2)*Da(2)*pow(eta*Gamma/this->getTimeStep(),_p)*DgammaDF(i,j)
- - 1/a(2)*_p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(eta*DGammaDF(i,j)/this->getTimeStep()));
- }
- }
-
- // ddRdTT, ddRdTF -> make all these -- CHANGE!!
- static STensor3 ddRdTF;
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- ddRdTF(i,j) = (dRdF(i,j)-dRdF_0(i,j))/(T-T0);
- }
- // STensorOperation::zero(ddRdTF); // ddGammaDTF is very difficult to estimate
-
- static double ddRdTT;
- ddRdTT = 2*dRdT/(T-T0) -2*(R-R0)/pow((T-T0),2); // 0.; // ddGammaDTT is very difficult to estimate
-
-
- // dMedF is available
- // const STensor43 dMedF = q1->_DModMandelDF;
-
- if(dWmdF!=NULL){
- STensor3 dmechSourceTVPdF;
-
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- dmechSourceTVPdF(i,j) = - (dRdF(i,j) - T*ddRdTF(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(k,l) + Me(i,j)*dDpdF(i,j,k,l)
- - (dXdF(i,j,k,l) - T*ddXdTF(i,j,k,l))*alpha1(k,l) - tempX(i,j)*dAlphadF(i,j,k,l)/this->getTimeStep();
- }
- }
-
- *dWmdF = dmechSourceTVPdF;
- }
-
- if(dWmdT!=NULL){
- double dmechSourceTVPdT(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;
- }
-}; */
-
-/*void mlawNonLinearTVP::getFreeEnergyTVM(const IPNonLinearTVP *q0, const 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;
- const STensor3& ddMedTT(0.); // CHANGE!!
- const STensor3& ddXdTT(0.); // CHANGE!!
-
- STensor3 devMe, devX, dDevXdT, dDevMeDT, ddDevXdTT, ddDevMeDTT;
- double trMe, trX, dTrXdT, dTrMeDT, ddTrXdTT, ddTrMedTT;
- STensorOperation::decomposeDevTr(q1->_ModMandel,devMe,trMe);
- STensorOperation::decomposeDevTr(q1->_backsig,devX,trX);
- STensorOperation::decomposeDevTr(q1->_DbackSigDT,dDevXdT,dTrXdT);
- STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
- STensorOperation::decomposeDevTr(ddXdTT,ddDevXdTT,ddTrXdTT);
- STensorOperation::decomposeDevTr(ddMedTT,ddDevMeDTT,ddTrMedTT);
-
- // 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 -> make IP for this CHANGE!!!
-
- // get Isotropic Hardening stuff
- double eta(0.),Deta(0.),DetaDT(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);
- double etaOverDt = eta/this->getTimeStep();
- double viscoTerm = etaOverDt*Gamma;
-
- double R1(0.), R2(0.), R3(0.), dR1dT(0.), dR2dT(0.), dR3dT(0.), ddR1dTT(0.), ddR2dTT(0.), ddR3dTT(0.);
- R1 = -a(1)/a(2) * (trMe-trX);
- dR1dT = (-DaDT(1)/a(2) + a(1)/pow(a(2),2) *DaDT(2))*(trMe-trX) + (-a(1)/a(2))*(dTrMeDT-dTrXdT);
- R2 = -a(0)/a(2);
- dR2dT = (-DaDT(0)/a(2) + a(0)/pow(a(2),2)*DaDT(2));
- if (Gamma>0 and etaOverDt>0){
- R3 -= 1/a(2)*pow(eta*Gamma/this->getTimeStep(),_p);
- dR3dT += (1/pow(a(2),2)*DaDT(2)*pow(eta*Gamma/this->getTimeStep(),_p)
- - 1/a(2)*_p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(DetaDT*Gamma/this->getTimeStep()
- + eta*dGammaDT/this->getTimeStep()) );
- }
-
- // 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);
-
- 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) + X(i,j)*DgammaDT*1/(pow(kk*Hb,2))*dHbdT;
- }
-
-
- // finally
- STensorOperation::doubleContractionSTensor3(X,alpha,psiVP);
- psiVP *= 0.5;
- psiVP += 0.5*( 1/_HR[0]*R1*a(1)/a(2) + 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*(0.);
- }
- DpsiVPdT += 0.5*( 1/_HR[0]*(dR1dT*a(1)/a(2) + R1*(DaDT(1)/a(2) - a(1)/pow(a(2),2) *DaDT(2))) + 2/_HR[1]*R2*dR2dT + 2/_HR[2]*R3*dR3dT );
-
- DDpsiVPdTT += 0.; // 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;
- }
-
-};*/
-
-/* OLD PLasticity Correction step
-
- * // 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 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);
- getTVEdCorKirDT(q0,q1,T0,T);
-
- // NEW
- static STensor3 Me, Mepr, Kepr; // Ke is corKir
-
- getModifiedMandel(Ce,q0,q1);
- Me = q1->_ModMandel;
- Mepr = Me;
- Kepr = q1-> _kirchhoff;
- // const STensor3& dKeprDT = q1->getConstRefToDcorKirDT();
- // NEW
-
- // NEW
- static STensor3 devXn;
- static double trXn;
- STensorOperation::decomposeDevTr(q0->_backsig,devXn,trXn); // needed for chaboche model
- // NEW
-
- static STensor3 PhiPr;
- PhiPr = q1->_ModMandel; // _kirchhoff; //NEW
- 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.;
-
- // NEW
- // Change u and v inside iterations
- double u = 1.;
- double v = 1.;
- // NEW
-
- double A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
-
-
- double dDgammaDGamma = 0.;
- double Dgamma = 0.; // eqplastic strain
-
-// NEW
- double Gt = 0.; // Ge + kk*Hb/2.;
- double Kt = 0.; // Ke + kk*Hb/3.;
- 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;
- Gt = Ge + pow(kk,2) * Hb/(2*Cx);
- Kt = Ke + pow(kk,2) * Hb/(3*Cx);
-// NEW
-
- 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();
-
- // NEW
- // DEFINE THESE INSIDE PLASTIC ITERATIONS -> Because dgamma is available here
- // need to calculate m(0) and m(1), preferably like yield coeffs (a(i)), they are functions of gamma
- // make sure that Px is already multiplied to m(1) or something
- // put these inside the iterations? Because dgamma isnt available here -> may have to put it in IP?
- // add temp funcs for dHbdT
-
- // Cx = m(0)*(q0->_epspbarre + Dgamma) + m(1)*Px - 1/Hb * dHbdT + 1;
- // Gt = Ge + pow(kk,2) * Hb/(2*Cx);
- // Kt = Ke + pow(kk,2) * Hb/(3*Cx);
- // NEW
-
- // NEW
- double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A);
- // double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A); // -> the same
- // NEW
-
- 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);
- // NEW
- 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; // CHANGE THIS DERIVATIVE -> tensor subtracting scalar
- 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;
- // NEW
- 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;
- // NEW
- // PhiEq = PhiEqpr/u; // -> original
- // ptilde = ptildepr/v; // -> original
- devPhi = (devPhipr + devXn * (1-1/Cx));
- devPhi *= 1./u;
- PhiEq = sqrt(1.5*devPhi.dotprod());
- ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v;
- // NEW
-
- 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; // I THINK temperature dependency of viscosity should be present in eta, but this relation remains the same
- 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
- // NEW
- devPhi = (devPhipr + devXn * (1-1/Cx));
- devPhi *= 1./u; //NEW
- ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v; //NEW
- // NEW
-
- // 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 *= (kk*Hb*Gamma);
- // NEW
- DB += q0->_backsig;
- DB *= 1/Cx;
- DB -= q0->_backsig;
- // NEW
- q1->_backsig += DB; // update
-
- /*
- // corotationaal Kirchhoff stress
- q1->_kirchhoff = devPhi;
- q1->_kirchhoff += q1->_backsig;
-
- q1->_kirchhoff(0,0) += (ptilde);
- q1->_kirchhoff(1,1) += (ptilde);
- q1->_kirchhoff(2,2) += (ptilde);
-
- */
- /*
- // NEW
- q1->_ModMandel = devPhi;
- q1->_ModMandel += q1->_backsig;
-
- q1->_ModMandel(0,0) += (ptilde);
- q1->_ModMandel(1,1) += (ptilde);
- q1->_ModMandel(2,2) += (ptilde);
-
- // NEW
-
- }
- else{
- q1->getRefToDissipationActive() = false;
- }
- }
-
-
- const STensor3& KS = q1->_kirchhoff;
- // second Piola Kirchhoff stress
- // static STensor3 S;
- // STensorOperation::multSTensor3STensor43(KS,DlnDCe,S);
-
- // NEW
- getModifiedMandel(Ce, q0, q1); // update Mandel // this step is probably wrong -> Mandel is already being updated above in the NR loop. Should corKir come from Mandel then?
- // WHAT ???
- const STensor3& MS = q1->_ModMandel;
- // second Piola Kirchhoff stress
- static STensor3 S, Ceinv;
- STensorOperation::inverseSTensor3(Ce,Ceinv);
- STensorOperation::multSTensor3(Ceinv,MS,S);
- // NEW
-
- 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);
- }
-
-
- // defo energy
- // q1->getRefToElasticEnergy()=deformationEnergy(*q1);
- q1->getRefToViscousEnergyPart()=viscousEnergy(*q0,*q1)+q0->getConstRefToViscousEnergyPart();
- q1->getRefToPlasticEnergy() = q0->plasticEnergy();
- if (Gamma > 0){
- // NEW
- // double dotKSN = dot(KS,N);
- //q1->getRefToPlasticEnergy() += Gamma*dotKSN;
- double dotMSN = dot(MS,N);
- q1->getRefToPlasticEnergy() += Gamma*dotMSN; // = Dp:Me
- // NEW
- }
-
-
- // 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);
-
- // thermalEnergy -> assuming stiff = false
- double CpT;
- if (stiff){
- double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
- getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
- CpT = -T*DDpsiTVMdTT;
- }
- else{
- getCp(CpT,T);
- }
-
- q1->_thermalEnergy = CpT*T;
-
- // thermalSource
- if (this->getTimeStep() > 0.){
- thermalSource = -CpT*(T-T0)/this->getTimeStep();
- }
- else
- thermalSource = 0.;
-
- // NEW - mechSourceTVP
- double& Wm = q1->getRefToMechSrcTVP();
- double& Wm_TVE = q1->getRefToMechSrcTVE();
-
- // Have to do this to update DcorKirDT in IP, and mechSourceTVE and derivatives
- mechanicalSource = 0.;
- mlawNonLinearTVE::getMechSourceTVE(q0,q1,T0,T,DKDTsum,DGDTsum,_I4,&Wm_TVE);
- // STensor3 temp101 = q1->_DcorKirDT;
- mechanicalSource += Wm_TVE;
-
- STensor3 P_TVE(0.);
- SVector3 fluxT_TVE(0.);
- double thermalSource_TVE(0.);
- double mechanicalSource_TVE(0.);
- /* mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(q0->getConstRefToFe(), q1->getConstRefToFe(), P_TVE, q0, q1, Tangent, dFedF, dFedT,
- T0, T, gradT0, gradT, fluxT_TVE, dPdT, dfluxTdgradT, dfluxTdT, dfluxTdF,
- thermalSource_TVE, dthermalSourcedT, dthermalSourcedF,
- mechanicalSource, dmechanicalSourcedT, dmechanicalSourceF,
- stiff);
-
-
- getMechSourceTVP(q0,q1,T0,T,_I,_I4,_I,&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);
-
- // 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.;
- }
-
- */
\ No newline at end of file
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 074c235685f1464342f4e05f1a75d5052b929888..7ddfc5fff55a96e7349214256723bce6e42889d0 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -245,7 +245,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
int numCurlDof = ipvcur->getNumConstitutiveCurlVariable();
_stressLaw->stress(ipvcur,ipvprev,false,checkfrac,dTangent);
-#if 0
+#if 1
_stressLaw->stress(ipvcur,ipvprev,true,checkfrac,dTangent);
ipvcur->getConstRefToTangentModuli().print("dPdE Analytique");
for (int k=0; k< numNonLocalVars; k++)
@@ -607,7 +607,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
}
}
-#if 1
+#if 0
ipvcur->getConstRefToTangentModuli().print("dPdE Numerique");
for (int k=0; k< numNonLocalVars; k++)
{