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

Added tension-compression asymmetry to TVP

parent 693192a2
No related branches found
No related tags found
1 merge request!406Correct NLbranchOption form 2 to 3; Added the memberfunction setNLBranchTYPE...
Show whitespace changes
Inline Side-by-side
NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
+ 89
− 25
View file @ 1c1e02ce
......@@ -381,10 +381,18 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
// update A, B -> ExtraBranch calculations are within
double Ke(0.), Ge(0.), Ke_pr(0.), Ge_pr(0.);
static STensor43 Ge_Tensor_pr, Ge_Tensor, Bd_stiffnessTerm;
STensorOperation::zero(Ge_Tensor_pr);
STensorOperation::zero(Ge_Tensor);
static STensor43 Ge_Tensor_pr, Ge_Tensor, Bd_stiffnessTerm, Ge_TrEeTensor_pr, Ge_TrEeTensor;
STensorOperation::zero(Ge_Tensor_pr); STensorOperation::zero(Ge_Tensor);
STensorOperation::zero(Ge_TrEeTensor_pr); STensorOperation::zero(Ge_TrEeTensor);
if (_extraBranchNLType != TensionCompressionRegularisedType && _ExtraBranch_TVE_option != 3){
mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm);
}
else{ // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm,&Ge_TrEeTensor);
Ge_TrEeTensor_pr = Ge_TrEeTensor;
}
Ge_Tensor = _I4;
Ge_Tensor *= Ge*2; // Because the function does not do it
Ge_Tensor += Bd_stiffnessTerm;
......@@ -465,7 +473,12 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
devPhi = devPhipr;
// Initialise the rest
if (_extraBranchNLType != TensionCompressionRegularisedType && _ExtraBranch_TVE_option != 3){
getDho(Gamma,Cepr,Ceinvpr,Kepr,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv);
}
else{ // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
getDho(Gamma,Cepr,Ceinvpr,Kepr,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,&Ge_TrEeTensor);
}
double PhiEqpr2 = 1.5*devPhi.dotprod();
double PhiEqpr = sqrt(PhiEqpr2); // -> second invariant
......@@ -524,9 +537,9 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
DdevPhidGamma_RHS(i,j) += ( -0.5*2*_b*Ke*(-ptilde-Gamma*dPhiPdGamma)*_I(i,j) - 3*pow(kk,1)*Hb/Cxdev*devPhi(i,j) ); // pow(kk,2) DEBUG
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DdevPhidDgamma_RHS(i,j) += 0.5* 2.*_b/3.*Gamma*Dho(i,j,k,l)*_I(k,l)*Hp;
DdevPhidGamma_RHS(i,j) += 0.5* ( -3.*Ge_Tensor(i,j,k,l)*devPhi(k,l) + 2.*_b/3.*Gamma*Dho(i,j,k,l)*_I(k,l)*dPhiPdGamma );
DdevPhidGamma_RHS(i,j) += 0.5*Dho(i,j,k,l)*N(k,l); // DEBUG
DdevPhidDgamma_RHS(i,j) += 0.5* 2.*_b/3.*Gamma*(-Ge_TrEeTensor(i,j,k,l) + Dho(i,j,k,l))*_I(k,l)*Hp;
DdevPhidGamma_RHS(i,j) += 0.5* ( -3.*Ge_Tensor(i,j,k,l)*devPhi(k,l) + 2.*_b/3.*Gamma*(-Ge_TrEeTensor(i,j,k,l) + Dho(i,j,k,l))*_I(k,l)*dPhiPdGamma );
DdevPhidGamma_RHS(i,j) += 0.5*(-Ge_TrEeTensor(i,j,k,l) + Dho(i,j,k,l))*N(k,l); // DEBUG
}
}
for (int i=0; i<3; i++)
......@@ -610,8 +623,15 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
}
// Update Phi
if (_extraBranchNLType != TensionCompressionRegularisedType && _ExtraBranch_TVE_option != 3){
getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA,
Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv);
}
else{ // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA,
Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,&Ge_TrEeTensor);
}
PhiEq = sqrt(1.5*devPhi.dotprod());
// Update A
......@@ -655,8 +675,14 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
}
// Correct Phi
if (_extraBranchNLType != TensionCompressionRegularisedType && _ExtraBranch_TVE_option != 3){
getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA,
Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv);
}
else{ // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA,
Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,&Ge_TrEeTensor);
}
PhiEq = sqrt(1.5*devPhi.dotprod());
// Correct Normal, H = expGN
......@@ -827,6 +853,8 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
DcorKirDEe(i,j,k,l) += Ge_Tensor(i,j,p,q)*_Idev(p,q,k,l);
}
}
DcorKirDEe_pr += Ge_TrEeTensor_pr; // TC Assymmetry
DcorKirDEe += Ge_TrEeTensor; // TC Assymmetry
static STensor43 DcorKir_pr_DCepr, DcorKirDCepr;
STensorOperation::multSTensor43(DcorKirDEe_pr,DlnDCepr,DcorKir_pr_DCepr);
......@@ -940,7 +968,8 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
for (int p=0; p<3; p++)
for (int q=0; q<3; q++){
DdevphiDCepr_RHS(i,j,k,l) += 0.5* ( Dho(i,j,p,q)* ( 2.*_b*Gamma/3.*_I(p,q)*DphiPDCepr(k,l) ) );
DdevphiDCepr_RHS(i,j,k,l) += 0.5* ( Ge_TrEeTensor(i,j,p,q)*0.5*DlnDCepr(p,q,k,l) +
(-Ge_TrEeTensor(i,j,p,q) + Dho(i,j,p,q))* ( 2.*_b*Gamma/3.*_I(p,q)*DphiPDCepr(k,l) ) );
for (int r=0; r<3; r++)
for (int s=0; s<3; s++){
......@@ -1015,8 +1044,9 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
for (int p=0; p<3; p++)
for (int q=0; q<3; q++){
DdevphiDCepr_RHS(i,j,k,l) += (0.5* ( - 3.*Ge_Tensor(i,j,p,q)*devPhi(p,q)*DGDCepr(k,l) +
Dho(i,j,p,q)* ( N(p,q)*DGDCepr(k,l) + 2.*_b*Gamma/3.*_I(p,q)*DphiPDCepr(k,l) ) ));
DdevphiDCepr_RHS(i,j,k,l) += (0.5* ( - 3.*Ge_Tensor(i,j,p,q)*devPhi(p,q)*DGDCepr(k,l)
+ (-Ge_TrEeTensor(i,j,p,q) + Dho(i,j,p,q)) * ( N(p,q)*DGDCepr(k,l) + 2.*_b*Gamma/3.*_I(p,q)*DphiPDCepr(k,l) )
+ Ge_TrEeTensor(i,j,p,q) *0.5*DlnDCepr(p,q,k,l) ) );
for (int r=0; r<3; r++)
for (int s=0; s<3; s++){
......@@ -1332,7 +1362,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DdevphiDT_RHS(i,j) += 0.5*( -3*dGammaDT*Ge_Tensor(i,j,k,l)*devPhi(k,l)
+ Dho(i,j,k,l)*( dGammaDT*N(k,l) + 2.*_b/3.*dPhiPdT*_I(k,l) )
+ (-Ge_TrEeTensor(i,j,k,l) + Dho(i,j,k,l))*( dGammaDT*N(k,l) + 2.*_b/3.*dPhiPdT*_I(k,l) )
+ Ce(i,k)*DcorKirDT(k,l)*Ceinv(l,j) );
}
}
......@@ -1388,7 +1418,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DdevphiDT_RHS(i,j) += 0.5*( -3*dGammaDT*Ge_Tensor(i,j,k,l)*devPhi(k,l)
+ Dho(i,j,k,l)*( dGammaDT*N(k,l) + 2.*_b/3.*Gamma*dPhiPdT*_I(k,l) )
+ (-Ge_TrEeTensor(i,j,k,l) + Dho(i,j,k,l))*( dGammaDT*N(k,l) + 2.*_b/3.*Gamma*dPhiPdT*_I(k,l) )
+ Ce(i,k)*DcorKirDT(k,l)*Ceinv(l,j) );
}
}
......@@ -1768,7 +1798,10 @@ void mlawNonLinearTVENonLinearTVP::get_G1_Tensor(const STensor3& Cepr, const STe
void mlawNonLinearTVENonLinearTVP::getDho(const double& Gamma, const STensor3& Cepr, const STensor3& Ceinvpr, const STensor3& KS,
const double& Ke, const STensor43& Ge_Tensor,
const double& kk, const double& Hb, const double& Cxtr, const double& Cxdev,
const STensor3& expGN, const STensor43& dexpAdA, STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv) const{
const STensor3& expGN, const STensor43& dexpAdA,
STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv, STensor43* Ge_TrEeTensor) const{
// Ge_TrEeTensor -> Dependence of devCorKir on TrEe. DO NOT CHANGE ITS VALUE!
// Get H = expGN is a symmetric tensor, HT, Hinv
static STensor3 HT, Hinv, HTinv;
......@@ -1861,6 +1894,9 @@ void mlawNonLinearTVENonLinearTVP::getDho(const double& Gamma, const STensor3& C
for (int q=0; q<3; q++)
DcorKirDEe(i,j,k,l) += Ge_Tensor(i,j,p,q)*_Idev(p,q,k,l);
}
if (Ge_TrEeTensor!= NULL){
DcorKirDEe += (*Ge_TrEeTensor);
}
// Dho -> 2nd term
for (int i=0; i<3; i++)
......@@ -1924,12 +1960,19 @@ void mlawNonLinearTVENonLinearTVP::getDho(const double& Gamma, const STensor3& C
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
Dho2_u(i,j,k,l) += 3*Gamma*pow(kk,1)*Hb/Cxdev * _I4(i,j,k,l); // DEBUG
Dho2_u(i,j,k,l) += 3*Gamma*pow(kk,1)*Hb/Cxdev * _I4(i,j,k,l); // DEBUG pow(kk,2)
for (int p=0; p<3; p++)
for (int q=0; q<3; q++){
Dho2_u(i,j,k,l) += 1.5*Gamma*(Ge_Tensor(i,j,p,q)*_I4(p,q,k,l) - Dho(i,j,p,q)*_I4(p,q,k,l));
}
}
if (Ge_TrEeTensor!= NULL){
static STensor43 Ge_TrEeTensor_temp;
STensorOperation::zero(Ge_TrEeTensor_temp);
Ge_TrEeTensor_temp = (*Ge_TrEeTensor);
Ge_TrEeTensor_temp *= (1.5*Gamma);
Dho2_u += Ge_TrEeTensor_temp;
}
STensorOperation::inverseSTensor43(Dho2_u,Dho2_u_inv);
// Get Dho2_v, Dho2_v_inv
......@@ -1944,7 +1987,8 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
double& Ke, double& Ge, STensor43& Ge_Tensor,
double& ptilde, STensor3& devPhi,
STensor3& Phi, STensor3& N, STensor3& expGN, STensor43& dexpAdA,
STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv) const{
STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv,
STensor43* Ge_TrEeTensor) const{
// This function calculates Phi iteratively.
// At every iteration it will calculate Ee, corKir and N
......@@ -1978,7 +2022,12 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
// Initialise corKir
static STensor3 KS;
static STensor43 Bd_stiffnessTerm;
if (Ge_TrEeTensor == NULL){
mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm);
}
else{ // TC asymmetry
mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm,Ge_TrEeTensor);
}
Ge_Tensor = _I4;
Ge_Tensor *= Ge*2; // *2 because the function doesnt do it
Ge_Tensor += Bd_stiffnessTerm;
......@@ -2023,7 +2072,12 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
}
// Initialise Dho, Dho2, Dho2inv
if (Ge_TrEeTensor == NULL){
getDho(Gamma,Cepr,Ceinvpr,KS,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv);
}
else{ // TC asymmetry
getDho(Gamma,Cepr,Ceinvpr,KS,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,Ge_TrEeTensor);
}
// Initialise DPhi
static STensor3 DPhi;
......@@ -2091,7 +2145,12 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
Ee(i,j) = Eepr(i,j) - Gamma*N(i,j);
// update corKir
mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,false,Bd_stiffnessTerm);
if (Ge_TrEeTensor == NULL){
mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm);
}
else{ // TC asymmetry
mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm,Ge_TrEeTensor);
}
Ge_Tensor = _I4;
Ge_Tensor *= Ge*2; // *2 because the function doesnt do it
Ge_Tensor += Bd_stiffnessTerm;
......@@ -2118,7 +2177,12 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
J(i,j) += ( Phi(i,j) - MS(i,j) + q1->_backsig(i,j) );
// update Dho2, Dho2inv
if (Ge_TrEeTensor == NULL){
getDho(Gamma,Cepr,Ceinvpr,KS,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv);
}
else{ // TC asymmetry
getDho(Gamma,Cepr,Ceinvpr,KS,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,Ge_TrEeTensor);
}
// update J_tol
/*
......
NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.h
+ 4
− 2
View file @ 1c1e02ce
......@@ -27,12 +27,14 @@ class mlawNonLinearTVENonLinearTVP : public mlawNonLinearTVP{
double& Ke, double& Ge, STensor43& Ge_Tensor,
double& ptilde, STensor3& devPhi,
STensor3& Phi, STensor3& N, STensor3& expGN, STensor43& dexpAdA,
STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv) const;
STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv,
STensor43* Ge_TrEeTensor = NULL) const;
virtual void getDho(const double& Gamma, const STensor3& Cepr, const STensor3& Ceinvpr, const STensor3& KS,
const double& Ke, const STensor43& Ge_Tensor,
const double& kk, const double& Hb, const double& Cxtr, const double& Cxdev,
const STensor3& expGN, const STensor43& dexpAdA, STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv) const;
const STensor3& expGN, const STensor43& dexpAdA,
STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv, STensor43* Ge_TrEeTensor = NULL) const;
virtual void get_G1_Tensor(const STensor3& Cepr, const STensor3& expGN, const STensor43& DCeinvprDCepr, const STensor3& KS, STensor43& G1) const;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment