From 7746f65dc030adb6edabd331507efb238dcf8f4a Mon Sep 17 00:00:00 2001 From: FLE_Knight <ujwalkishore.jinaga@mail.polimi.it> Date: Wed, 24 May 2023 18:27:38 +0200 Subject: [PATCH] made some changes to mullins Effect, gives unacceptable results for now --- NonLinearSolver/internalPoints/ipField.h | 2 + .../internalPoints/ipMullinsEffect.cpp | 2 +- .../internalPoints/ipNonLinearTVE.cpp | 6 ++- .../internalPoints/ipNonLinearTVE.h | 1 + .../materialLaw/mlawNonLinearTVE.cpp | 53 +++++++++++++------ .../materialLaw/mlawNonLinearTVE.h | 6 +-- .../materialLaw/mlawNonLinearTVP.cpp | 19 ++++--- NonLinearSolver/materialLaw/mullinsEffect.cpp | 24 +++++---- NonLinearSolver/materialLaw/mullinsEffect.h | 10 ++-- dG3D/src/dG3DIPVariable.cpp | 25 ++++++++- dG3D/src/dG3DMaterialLaw.cpp | 2 +- 11 files changed, 103 insertions(+), 47 deletions(-) diff --git a/NonLinearSolver/internalPoints/ipField.h b/NonLinearSolver/internalPoints/ipField.h index 7fa320ea0..0948b41e1 100644 --- a/NonLinearSolver/internalPoints/ipField.h +++ b/NonLinearSolver/internalPoints/ipField.h @@ -333,6 +333,8 @@ class IPField : public elementsField { EQUIVALENT_EIGENSTRESS,EQUIVALENT_EIGENSTRAIN, EIGENSTRESS_XX, EIGENSTRESS_YY, EIGENSTRESS_ZZ, EIGENSTRESS_XY, EIGENSTRESS_XZ, EIGENSTRESS_YZ, EIGENSTRAIN_XX, EIGENSTRAIN_YY, EIGENSTRAIN_ZZ, EIGENSTRAIN_XY, EIGENSTRAIN_XZ, EIGENSTRAIN_YZ, + // moduli scalers e.g. mullins effect + MAX_DEFO_ENERGY, BULK_SCALE_FACTOR, SHEAR_SCALE_FACTOR, UNDEFINED}; enum Operator { MEAN_VALUE=1, MIN_VALUE, MAX_VALUE, CRUDE_VALUE}; static std::string ToString(const int i); diff --git a/NonLinearSolver/internalPoints/ipMullinsEffect.cpp b/NonLinearSolver/internalPoints/ipMullinsEffect.cpp index 5ae3728b9..f53d5384b 100644 --- a/NonLinearSolver/internalPoints/ipMullinsEffect.cpp +++ b/NonLinearSolver/internalPoints/ipMullinsEffect.cpp @@ -10,7 +10,7 @@ #include "ipMullinsEffect.h" -IPMullinsEffect::IPMullinsEffect(): eta(0), DetaDpsi(0), DDetaDpsipsi(0), DetaDT(0.), DDetaDTT(0.) +IPMullinsEffect::IPMullinsEffect(): eta(1.0), DetaDpsi(0.), DDetaDpsipsi(0.), DetaDT(0.), DDetaDTT(0.) { } diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp b/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp index 9a047d211..fce0b2dc4 100644 --- a/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp +++ b/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp @@ -18,7 +18,7 @@ IPNonLinearTVE::IPNonLinearTVE(const J2IsotropicHardening* comp, const kinematicHardening* kin, const int N, const mullinsEffect* mullins): IPHyperViscoElastoPlastic(comp,trac,kin,N), // IPHyperViscoElastic(N), _T(298.15),_devDE(0.),_trDE(0.),_DE(0.),_DcorKirDT(0.),_DDcorKirDTT(0.), _mechSrcTVE(0.),_DmechSrcTVEdT(0.),_DmechSrcTVEdF(0.), - _thermalEnergy(0.),_dtShift(0.),_DdtShiftDT(0.),_DDdtShiftDDT(0.){ + _thermalEnergy(0.),_dtShift(0.),_DdtShiftDT(0.),_DDdtShiftDDT(0.), _psiMax(0.){ _devOGammai.clear(); _trOGammai.clear(); _devDOGammaiDT.clear(); @@ -66,7 +66,7 @@ IPNonLinearTVE::IPNonLinearTVE(const IPNonLinearTVE& src): IPHyperViscoElastoPla _mechSrcTVE(src._mechSrcTVE),_DmechSrcTVEdT(src._DmechSrcTVEdT),_DmechSrcTVEdF(src._DmechSrcTVEdF), _C(src._C),_D(src._D), _E(src._E), _F(src._F), _G(src._G), _thermalEnergy(src._thermalEnergy), - _dtShift(src._dtShift),_DdtShiftDT(src._DdtShiftDT),_DDdtShiftDDT(src._DDdtShiftDDT){ + _dtShift(src._dtShift),_DdtShiftDT(src._DdtShiftDT),_DDdtShiftDDT(src._DDdtShiftDDT), _psiMax(src._psiMax){ if (src._ipMullinsEffect != NULL) _ipMullinsEffect = dynamic_cast<IPMullinsEffect*>(src._ipMullinsEffect->clone()); @@ -110,6 +110,7 @@ IPNonLinearTVE& IPNonLinearTVE::operator=(const IPVariable& src) _dtShift = psrc->_dtShift; _DdtShiftDT = psrc->_DdtShiftDT; _DDdtShiftDDT = psrc->_DDdtShiftDDT; + _psiMax = psrc->_psiMax; if ( psrc->_ipMullinsEffect != NULL) { if (_ipMullinsEffect == NULL){ @@ -173,6 +174,7 @@ void IPNonLinearTVE::restart(){ restartManager::restart(_dtShift); restartManager::restart(_DdtShiftDT); restartManager::restart(_DDdtShiftDDT); + restartManager::restart(_psiMax); if (_ipMullinsEffect != NULL) restartManager::restart(_ipMullinsEffect); diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVE.h b/NonLinearSolver/internalPoints/ipNonLinearTVE.h index 4b55144a0..82c653b90 100644 --- a/NonLinearSolver/internalPoints/ipNonLinearTVE.h +++ b/NonLinearSolver/internalPoints/ipNonLinearTVE.h @@ -62,6 +62,7 @@ class IPNonLinearTVE : public IPHyperViscoElastoPlastic{ // Mullins Effect IP IPMullinsEffect* _ipMullinsEffect; + double _psiMax; // maximum strain energy public: IPNonLinearTVE(const J2IsotropicHardening* comp, diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp index a814172f6..706f2ff61 100644 --- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp @@ -534,11 +534,11 @@ double mlawNonLinearTVE::setTemp(const double T0, const double T1, const int par return T_set; } -double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const double Tc) const{ +double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const double T0, const double T) const{ // Update the Properties to the current temperature (see below which ones are being used) double KT, GT, AlphaT; - getK(&q,KT,Tc); getG(&q,GT,Tc); getAlpha(AlphaT,Tc); + getK(&q,KT,T); getG(&q,GT,T); getAlpha(AlphaT,T); double Psy_mech = 0.; @@ -548,7 +548,7 @@ double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const dou STensorOperation::decomposeDevTr(q._Ee,devEe,trEe); // Get thermal strain - double Eth = 3.*AlphaT*(Tc-_Tinitial); + double Eth = 3.*AlphaT*(T-_Tinitial); // Get Equilibrium freeEnergy_mech Psy_mech = KT*0.5*(trEe-Eth)*(trEe-Eth) + GT*STensorOperation::doubledot(devEe,devEe); @@ -558,45 +558,57 @@ double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const dou std::vector<double> KiT,GiT,AlphaiT; KiT.resize(_N,0.); GiT.resize(_N,0.); AlphaiT.resize(_N,0.); - getKi(&q,KiT,Tc); getGi(&q,GiT,Tc); getAlphai(AlphaiT,Tc); + getKi(&q,KiT,T); getGi(&q,GiT,T); getAlphai(AlphaiT,T); for (int i=0; i<_N; i++){ static STensor3 devEebranch; devEebranch = devEe; + + // viscous overstrain = _OGammai + // viscous strain (_Gammai) = viscous overstrain - Ee + double trGamma = q._trOGammai[i]; + trGamma -= trEe; + static STensor3 devGamma; + devGamma = q._devOGammai[i]; + devGamma -= devEe; + + devEebranch -= devGamma; // temporary tensor - // devEebranch -= q._A[i]; - // Psy_mech += KiT[i]*0.5*(trEe-q._B[i])*(trEe-q._B[i])+GiT[i]*STensorOperation::doubledot(devEebranch,devEebranch); - - double trGamma = q._trOGammai[i]; - static STensor3 devGamma = q._devOGammai[i]; - devEebranch -= devGamma; if (_modelFlag == 2){ // Additional volumetric terms - trGamma += 3*(AlphaT-AlphaiT[i])*(Tc-_Tinitial); + trGamma += 3*(AlphaT-AlphaiT[i])*(T-_Tinitial); } Psy_mech += KiT[i]*0.5*(trEe-trGamma)*(trEe-trGamma) + GiT[i]*STensorOperation::doubledot(devEebranch,devEebranch); } } - + + // Get thermal Energy + // Psy_mech += _Cp*( (T-T0) - T*log(T/T0) ); // this probably does nothing since we are interested only in the elastic energy + return Psy_mech; } void mlawNonLinearTVE::calculateMullinsEffectScaler(const IPNonLinearTVE* q0, IPNonLinearTVE *q1, const double T) const{ + // Initialise + double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.; + + // Msg::Error(" Inside calculateMullinsEffectScaler, before updating eta = %e!!",eta); + // update eta - Mullins Effect Scalar if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){ - _mullinsEffect->mullinsEffectScaling(q0->_elasticEnergy, *q0->_ipMullinsEffect, q1->_elasticEnergy,*q1->_ipMullinsEffect); + _mullinsEffect->mullinsEffectScaling(q1->_psiMax, *q0->_ipMullinsEffect, q0->_elasticEnergy, q1->_elasticEnergy, *q1->_ipMullinsEffect); } - + // get eta - Mullins Effect Scalar (at unloading intitiation-> eta = 1.; unloading end-> eta = eta_min; reloading -> eta = eta_min;) - double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.; if (q1->_ipMullinsEffect != NULL){ eta = q1->_ipMullinsEffect->getEta(); DetaDpsi = q1->_ipMullinsEffect->getDetaDpsi(); DetaDT = q1->_ipMullinsEffect->getDetaDT(); DDetaDTT = q1->_ipMullinsEffect->getDDetaDTT(); } + // Msg::Error(" Inside calculateMullinsEffectScaler, after getting and updating eta = %e!!",eta); // update Bulk/Shear PropertyScaleFactor q1->getRefToElasticBulkPropertyScaleFactor() = eta; @@ -1029,7 +1041,7 @@ void mlawNonLinearTVE::ThermoViscoElasticPredictor(const STensor3& Ee, const STe double& DKe, double& DGe, double& DKde, double& DGde, double& KTsum, double& GTsum, double& DKDTsum, double& DGDTsum, const double T0, const double T1, const bool calculateStress) const{ - + // Mullins Effect Scaler -> should be before any calls to get functions if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){ mlawNonLinearTVE::calculateMullinsEffectScaler(q0, q1, T1); @@ -1502,7 +1514,14 @@ void mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic( STensorOperation::multSTensor3(F,secondPK,P); // 1st PK // elastic energy - q1->_elasticEnergy = freeEnergyMechanical(*q1,T); // deformationEnergy(*q1,T); + q1->_elasticEnergy = freeEnergyMechanical(*q1,T0,T); + if (q1->_elasticEnergy > q0->_psiMax){ + q1->_psiMax = q1->_elasticEnergy; + } + else{ + q1->_psiMax = q0->_psiMax; + } + // Msg::Error(" Inside TVE predCorr, after updating q1->_elasticEnergy = %e!!",q1->_elasticEnergy); // thermal energy q1->_thermalEnergy = CpT*T; diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.h b/NonLinearSolver/materialLaw/mlawNonLinearTVE.h index 90ae75b93..26073b30d 100644 --- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.h +++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.h @@ -77,7 +77,7 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{ double *t_shift = NULL, double *Ddt_shiftDT = NULL, double *DDdt_shiftDT = NULL, const int opt_time = 1, const int opt_order = 2) const; virtual double linearInterpTemp(const double t_interp, const double t0, const double t1, const double T0, const double T1) const; - virtual double freeEnergyMechanical(const IPNonLinearTVE& q, const double Tc) const ; + virtual double freeEnergyMechanical(const IPNonLinearTVE& q, const double T0, const double T) const ; virtual double setTemp(const double T0, const double T1, const int para) const ; virtual void calculateMullinsEffectScaler(const IPNonLinearTVE* q0, IPNonLinearTVE *q1, const double T) const; @@ -229,6 +229,7 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{ virtual void getK(const IPNonLinearTVE *q1, double& KT, const double Tc, double *dK = NULL, double *ddK = NULL) const{ double K = _K*q1->getConstRefToElasticBulkPropertyScaleFactor(); + // Msg::Error(" Inside getK, after getting and updating eta; K = %e!!",K); KT = K*_temFunc_K->getVal(Tc); if(dK!=NULL){ *dK = K*_temFunc_K->getDiff(Tc); // Here, *dK can be imagined as *(dK+0) or dK[0] if dK were a vector. See dKi below. @@ -274,13 +275,12 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{ // The 3rd argument is a pointer to a vector. virtual void getKi(const IPNonLinearTVE *q1, std::vector<double>& Ki_T, const double Tc, std::vector<double>* dKi=NULL, std::vector<double>* ddKi=NULL) const{ - std::vector<double> Ki; Ki.clear(); Ki.resize(_N,0.); for (int i=0; i<_Ki.size();i++){ Ki[i] = _Ki[i]*q1->getConstRefToElasticBulkPropertyScaleFactor(); } - + // Msg::Error(" Inside getKi, after getting and updating eta; Ki[0] = %e!!",Ki[0]); Ki_T.resize(_N,0.); for (int i=0; i<_Ki.size();i++){ Ki_T[i] = Ki[i]*_temFunc_Ki[i]->getVal(Tc); diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp index 873300d96..a67d67427 100644 --- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp @@ -364,7 +364,7 @@ void mlawNonLinearTVP::getFreeEnergyTVM(const IPNonLinearTVP *q0, IPNonLinearTVP const std::vector<STensor3>& devDDOiDTT = q1->getConstRefToDevDDOiDTT(); const std::vector<double>& trDDOiDTT = q1->getConstRefToTrDDOiDTT(); - psiVE = mlawNonLinearTVE::freeEnergyMechanical(*q1,T); + psiVE = mlawNonLinearTVE::freeEnergyMechanical(*q1,T0,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 @@ -3249,7 +3249,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3& double CpT; if (stiff){ double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT(); - getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT); + // getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT); // CpT = -T*DDpsiTVMdTT; getCp(CpT,T); } @@ -3281,14 +3281,21 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3& double& psiTVM = q1->getRefToFreeEnergyTVM(); // getFreeEnergyTVM(q0,q1,T0,T,&psiTVM,NULL); - q1->_elasticEnergy = mlawNonLinearTVE::freeEnergyMechanical(*q1,T); // deformationEnergy(*q1,T); + q1->_elasticEnergy = mlawNonLinearTVE::freeEnergyMechanical(*q1,T0,T); + if (q1->_elasticEnergy > q0->_psiMax){ + q1->_psiMax = q1->_elasticEnergy; + } + else{ + q1->_psiMax = q0->_psiMax; + } + // Msg::Error(" Inside TVP after correction, after updating q1->_elasticEnergy = %e!!",q1->_elasticEnergy); if (stiff){ - const double& DDpsiTVMdTT_0 = q0->getConstRefToDDFreeEnergyTVMdTT(); - double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT(); - getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT); + // 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; diff --git a/NonLinearSolver/materialLaw/mullinsEffect.cpp b/NonLinearSolver/materialLaw/mullinsEffect.cpp index 47f32a8b9..14bda0551 100644 --- a/NonLinearSolver/materialLaw/mullinsEffect.cpp +++ b/NonLinearSolver/materialLaw/mullinsEffect.cpp @@ -85,44 +85,46 @@ void negativeExponentialScaler::setTemperatureFunction_r(const scalarFunction& T _temFunc_r = Tfunc.clone(); } +void negativeExponentialScaler::setTemperatureFunction_m(const scalarFunction& Tfunc){ + if (_temFunc_m != NULL) delete _temFunc_m; + _temFunc_m = Tfunc.clone(); +} + void negativeExponentialScaler::createIPVariable(IPMullinsEffect* &ipv) const { if(ipv != NULL) delete ipv; ipv = new IPMullinsEffect(); } -void negativeExponentialScaler::mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, IPMullinsEffect &ipv) const{ +void negativeExponentialScaler::mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, IPMullinsEffect &ipv) const{ // at unloading intitiation-> eta = 1.; unloading end-> eta = eta_min; reloading -> eta = eta_min; double eta0 = ipvprev.getEta(); - double psi = _psi; - double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.; - if (psi<psi0){ - eta = 1. - _r*( 1. - exp( - sqrt(_m*(psi0-psi)) ) ); + if (_psi<psi_max && _psi<_psi0){ + eta = 1. - _r*( 1. - exp( - sqrt(_m*(psi_max-_psi)) ) ); } else{ - eta = eta0; + eta = 1.; } ipv.set(eta,DetaDpsi,DDetaDpsipsi,DetaDT,DDetaDTT); } -void negativeExponentialScaler::mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, double T, IPMullinsEffect &ipv) const{ +void negativeExponentialScaler::mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, double T, IPMullinsEffect &ipv) const{ double eta0 = ipvprev.getEta(); - double psi = _psi; double r = _r*_temFunc_r->getVal(T); double m = _m*_temFunc_m->getVal(T); - double eta = 0., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.; + double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.; - if (psi<psi0){ - eta = 1. - r*( 1. - exp( - sqrt(m*(psi0-psi)) ) ); + if (_psi<psi_max && _psi<_psi0){ + eta = 1. - r*( 1. - exp( - sqrt(m*(psi_max-_psi)) ) ); } else{ eta = eta0; diff --git a/NonLinearSolver/materialLaw/mullinsEffect.h b/NonLinearSolver/materialLaw/mullinsEffect.h index 9e06ee651..dd10cdd9f 100644 --- a/NonLinearSolver/materialLaw/mullinsEffect.h +++ b/NonLinearSolver/materialLaw/mullinsEffect.h @@ -33,14 +33,14 @@ class mullinsEffect { virtual void createIPVariable(IPMullinsEffect* &ipv) const=0; virtual void initLaws(const std::map<int,mullinsEffect*> &maplaw)=0; virtual const bool isInitialized() {return _initialized;} - virtual void mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, IPMullinsEffect &ipv) const=0; - virtual void mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, double T, IPMullinsEffect &ipv) const=0; + virtual void mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, IPMullinsEffect &ipv) const=0; + virtual void mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, double T, IPMullinsEffect &ipv) const=0; virtual mullinsEffect * clone() const=0; #endif }; -// eta = 1 - r* ( 1 - exp(-sqrt(m*(psi0_max-psi0))) ) -> Ricker et al, 2021 (modified Zuniga, 2005) +// eta = 1 - r* ( 1 - exp(-sqrt(m*(psi_max-psi))) ) -> Ricker et al, 2021 (modified Zuniga, 2005) class negativeExponentialScaler : public mullinsEffect{ #ifndef SWIG protected: @@ -61,8 +61,8 @@ class negativeExponentialScaler : public mullinsEffect{ virtual scalingFunction getType() const{return mullinsEffect::negativeExponentialScaler; }; virtual void createIPVariable(IPMullinsEffect* &ipv) const; virtual void initLaws(const std::map<int,mullinsEffect*> &maplaw) {}; //nothing now, we will see if we use the mapping - virtual void mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, IPMullinsEffect &ipv) const; - virtual void mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, double T, IPMullinsEffect &ipv) const; + virtual void mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, IPMullinsEffect &ipv) const; + virtual void mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, double T, IPMullinsEffect &ipv) const; virtual mullinsEffect * clone() const; #endif // SWIG }; diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp index 88095add6..e48cf6989 100644 --- a/dG3D/src/dG3DIPVariable.cpp +++ b/dG3D/src/dG3DIPVariable.cpp @@ -5251,6 +5251,18 @@ double NonLinearTVEDG3DIPVariable::get(const int i) const{ else if (i == IPField::THERMALFLUX_Z){return getConstRefToFlux()[0](2);} else if (i == IPField::FIELD_SOURCE){return getConstRefToFieldSource()(0);} else if (i == IPField::MECHANICAL_SOURCE){return getConstRefToMechanicalSource()(0);} + else if (i == IPField::DEFO_ENERGY) + { + return _ipViscoElastic->_elasticEnergy; + } + else if (i == IPField::MAX_DEFO_ENERGY) + { + return _ipViscoElastic->_psiMax; + } + else if (i == IPField::BULK_SCALE_FACTOR) + { + return _ipViscoElastic->getConstRefToElasticBulkPropertyScaleFactor(); + } else { @@ -5453,7 +5465,18 @@ double NonLinearTVPDG3DIPVariable::get(const int i) const{ { return _ipViscoElastoPlastic->getConstRefToKinematicHardening().getR(); } - + else if (i == IPField::DEFO_ENERGY) + { + return _ipViscoElastoPlastic->_elasticEnergy; + } + else if (i == IPField::MAX_DEFO_ENERGY) + { + return _ipViscoElastoPlastic->_psiMax; + } + else if (i == IPField::BULK_SCALE_FACTOR) + { + return _ipViscoElastoPlastic->getConstRefToElasticBulkPropertyScaleFactor(); + } else { return dG3DIPVariable::get(i); diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp index 1c42e2af3..f7f61561b 100644 --- a/dG3D/src/dG3DMaterialLaw.cpp +++ b/dG3D/src/dG3DMaterialLaw.cpp @@ -616,7 +616,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV } } -#if 1 +#if 0 ipvcur->getConstRefToDeformationGradient().print("F Numerique"); // FLE ipvcur->getConstRefToFirstPiolaKirchhoffStress().print("P Numerique"); // FLE ipvcur->getConstRefToTangentModuli().print("dPdE Numerique"); -- GitLab