Skip to content
Snippets Groups Projects
Commit b79dca56 authored by Julien Leclerc's avatar Julien Leclerc
Browse files

modify internal variabes functions chi and cft

parent 0214b099
No related branches found
No related tags found
2 merge requests!50Coal jl2,!49Transfer Coal jl2 into master
Show whitespace changes
Inline Side-by-side
NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.cpp
+ 51
− 26
View file @ b79dca56
......@@ -20,7 +20,7 @@ mlawNonLocalPorousThomasonLaw::mlawNonLocalPorousThomasonLaw(const int num,const
const double tol, const bool matrixbyPerturbation, const double pert) :
mlawNonLocalPorosity(num,E,nu,rho,fVinitial,j2IH,cLLaw,tol,matrixbyPerturbation,pert),
_lambda0(lambda0), _n(10.),
_chi_failure(0.99), _chi_saturated(0.9995), _chi_decay(_chi_saturated-_chi_failure), _fV_Failure(0.95),
_chi_switch(0.99), _chi_saturated(0.9995), _chi_decay(_chi_saturated-_chi_switch), _fV_Failure(0.95),
_initialGuessMethod(0), _maxite(20), _theta(1.), _kappa(1.)
{
// modify _coalescenceLaw from default type created by mlawNonLocalPorosity::mlawNonLocalPorosity to adapted one
......@@ -35,7 +35,7 @@ mlawNonLocalPorousThomasonLaw::mlawNonLocalPorousThomasonLaw(const int num,const
mlawNonLocalPorousThomasonLaw::mlawNonLocalPorousThomasonLaw(const mlawNonLocalPorousThomasonLaw &source) :
mlawNonLocalPorosity(source),
_lambda0(source._lambda0), _n(source._n),
_chi_failure(source._chi_failure), _chi_saturated(source._chi_saturated), _chi_decay(source._chi_decay), _fV_Failure(source._fV_Failure),
_chi_switch(source._chi_switch), _chi_saturated(source._chi_saturated), _chi_decay(source._chi_decay), _fV_Failure(source._fV_Failure),
_initialGuessMethod(source._initialGuessMethod), _maxite(source._maxite),
_theta(source._theta), _kappa(source._kappa)
{
......@@ -82,20 +82,29 @@ void mlawNonLocalPorousThomasonLaw::createIPState(IPStateBase* &ips,const bool*
};
double mlawNonLocalPorousThomasonLaw::computeLigamentRatio(const double hatP, const double tildefV) const
void mlawNonLocalPorousThomasonLaw::computeLigamentRatio(const double hatP, const double tildefV, double& chi, double& dChidhatP, double& dChidtildefV) const
{
double Chi_prime = pow(1.5*tildefV*_lambda0 , 1./3.)*exp( hatP/3.*_kappa);
if(Chi_prime < _chi_failure)
dChidhatP = Chi_prime*_kappa/3.;
dChidtildefV = Chi_prime/tildefV/_kappa/3.;
if(Chi_prime <= _chi_switch)
{
return Chi_prime;
chi = Chi_prime;
}
else
{
return (_chi_saturated + (_chi_failure-_chi_saturated)*exp(-(Chi_prime-_chi_failure)/_chi_decay));
double expfactor = exp((_chi_switch-Chi_prime)/_chi_decay);
chi = _chi_saturated - _chi_decay*expfactor;
dChidhatP *= expfactor;
dChidtildefV *= expfactor;
}
};
void mlawNonLocalPorousThomasonLaw::computeConcentrationFactor(double& Cft, double& dCftDChi, const double Chi) const
void mlawNonLocalPorousThomasonLaw::computeConcentrationFactor(const double Chi, double& Cft, double& dCftDChi) const
{
// Compute ConcentrationFactor and its derivatives
double alpha = 0.1;
......@@ -109,6 +118,8 @@ void mlawNonLocalPorousThomasonLaw::computeConcentrationFactor(double& Cft, doub
dCftDChi = (-2.)*ChiBis*Cft_2 + Cft_1*(2.*alpha*(1./ChiBis - 1.)*(-1./ChiBis/ChiBis) - 0.5*beta*pow(ChiBis,-1.5));
dCftDChi *= 1.5;
}
else
{
......@@ -121,18 +132,23 @@ void mlawNonLocalPorousThomasonLaw::computeConcentrationFactor(double& Cft, doub
}
};
bool mlawNonLocalPorousThomasonLaw::yieldFunction(const double tau_eq, const double p, const IPNonLocalPorosity *ipv, const double fVstar) const
bool mlawNonLocalPorousThomasonLaw::yieldFunction(const double tau_eq, const double p, const IPNonLocalPorosity *ipv, const double fVtilde) const
{
double dChidhatP;
double dChidtildefV;
// Check yield criterion
// Chi
const double hatP = ipv->getMatrixPlasticStrain();
//const double chi = q0->getConstRefToIPCoalescence().getLigamentRatio();
const double chi = computeLigamentRatio(hatP, fVstar);
double chi = ipv->getConstRefToIPCoalescence().getLigamentRatio();
computeLigamentRatio(hatP, fVtilde, chi, dChidhatP, dChidtildefV);
// Cft
double Cft = 0.;
double dc = 0.;
computeConcentrationFactor(Cft, dc, chi);
computeConcentrationFactor(chi, Cft, dc);
// Yield
double yield = ipv->getConstRefToIPJ2IsotropicHardening().getR();
......@@ -151,7 +167,7 @@ void mlawNonLocalPorousThomasonLaw::computeResidual(fullVector<double> &res, dou
const double tpr_eq, const double p_pr,
const double Chin, const double fVn,
const double tildefV, const double Hatpn,
const double yield0, double& dChidChiPrime,
const double yield0, double& dChidtildefV,
double& dfVdDeltaD, double& dfVdDeltaQ, double& dfVdDeltaHatP) const
{
double G = _mu;
......@@ -181,6 +197,12 @@ void mlawNonLocalPorousThomasonLaw::computeResidual(fullVector<double> &res, dou
// Compute DeltaChi and Cft
double Chi1, dChidDeltaHatP;
computeLigamentRatio(Hatpn+DeltaHatP, tildefV, Chi1, dChidDeltaHatP, dChidtildefV);
double Cft, dCftdChi;
computeConcentrationFactor(Chi1, Cft, dCftdChi);
DeltaChi = Chi1-Chin;
/*
double Chi_prime = pow(1.5*tildefV*_lambda0 , 1./3.)*exp( (Hatpn+DeltaHatP)/3.*_kappa);
double dChidDeltaHatP = pow(1.5*tildefV*_lambda0 , 1./3.) * exp( (Hatpn+DeltaHatP)/3.*_kappa) /3.*_kappa;
double Chi1 = Chi_prime;
......@@ -195,10 +217,10 @@ void mlawNonLocalPorousThomasonLaw::computeResidual(fullVector<double> &res, dou
dChidChiPrime = (_chi_saturated-_chi_failure)*exp(-(Chi_prime-_chi_failure)/_chi_decay)/_chi_decay;
}
dChidDeltaHatP *= dChidChiPrime;
DeltaChi = Chi1-Chin;
double Cft, dCftdChi;
computeConcentrationFactor(Cft, dCftdChi, Chi1);
*/
// Compute DeltafV
......@@ -296,7 +318,7 @@ bool mlawNonLocalPorousThomasonLaw::plasticCorrector(const STensor3& F1, const d
double DeltaChi = 0.;
// derivatives
double dChidChiPrime = 1.;
double dChidtildefV = 1.;
double dfVdDeltaD = 0.;
double dfVdDeltaQ = 0.;
double dfVdDeltaHatP = 0.;
......@@ -316,7 +338,9 @@ bool mlawNonLocalPorousThomasonLaw::plasticCorrector(const STensor3& F1, const d
const double fV0 = q0->getLocalPorosity();
fV1 = fV0;
// chi
const double Chi0 = computeLigamentRatio(hatP0,fVtilde); // should be moved somewhere else before in yield testing
double Chi0;
double dChidHatP;
computeLigamentRatio(hatP0,fVtilde,Chi0,dChidHatP,dChidtildefV); // should be moved somewhere else before in yield testing
double& Chi1 = q1->getRefToIPCoalescence().getRefToLigamentRatio();
Chi1 = Chi0;
......@@ -358,7 +382,7 @@ bool mlawNonLocalPorousThomasonLaw::plasticCorrector(const STensor3& F1, const d
static fullMatrix<double> J(3,3); J.setAll(0.);
// Compute residues
computeResidual(res, res_norm, J, stiff, q0, q1, Delta_d, Delta_q, DeltaHatP, DeltafV, DeltaChi,
tau_preq, p_pr, Chi0, fV0, fVtilde, hatP0, yield0, dChidChiPrime,
tau_preq, p_pr, Chi0, fV0, fVtilde, hatP0, yield0, dChidtildefV,
dfVdDeltaD, dfVdDeltaQ, dfVdDeltaHatP);
// Iterative NR procedure
......@@ -471,7 +495,7 @@ bool mlawNonLocalPorousThomasonLaw::plasticCorrector(const STensor3& F1, const d
// Recompute residues before checking for convergence
computeResidual(res, res_norm, J, stiff, q0, q1, Delta_d, Delta_q, DeltaHatP, DeltafV, DeltaChi,
tau_preq, p_pr, Chi0, fV0, fVtilde, hatP0, yield0, dChidChiPrime,
tau_preq, p_pr, Chi0, fV0, fVtilde, hatP0, yield0, dChidtildefV,
dfVdDeltaD, dfVdDeltaQ, dfVdDeltaHatP);
// Control of iteration number
......@@ -482,7 +506,7 @@ bool mlawNonLocalPorousThomasonLaw::plasticCorrector(const STensor3& F1, const d
//Msg::Error(" Initial residu = %e, %e, %e", resinit(0), resinit(1), resinit(2));
Msg::Error(" Deltad = %e ,Deltaq = %e ,DeltaHatP = %e ,DeltafV = %e ,DeltaChi = %e", Delta_d, Delta_q, DeltaHatP, DeltafV, DeltaChi);
double Cft, dCftdChi;
computeConcentrationFactor(Cft, dCftdChi, Chi1);
computeConcentrationFactor(Chi1, Cft, dCftdChi);
Msg::Error(" tau_preq = %e ,p_pr = %e ,chi0 = %e ,fV0 = %e, Cft = %e", tau_preq/yield0, p_pr/yield0, Chi0, fV0, Cft);
Msg::Error("No convergence for plastic correction in coalescence !!");
return false;
......@@ -586,7 +610,7 @@ bool mlawNonLocalPorousThomasonLaw::plasticCorrector(const STensor3& F1, const d
if (stiff)
{
computeResidual(res, res_norm, J, stiff, q0, q1, Delta_d, Delta_q, DeltaHatP, DeltafV, DeltaChi,
tau_preq, p_pr, Chi0, fV0, fVtilde, hatP0, yield0, dChidChiPrime,
tau_preq, p_pr, Chi0, fV0, fVtilde, hatP0, yield0, dChidtildefV,
dfVdDeltaD, dfVdDeltaQ, dfVdDeltaHatP);
// Compute internal variables derivatives from residual derivatives
......@@ -624,8 +648,8 @@ bool mlawNonLocalPorousThomasonLaw::plasticCorrector(const STensor3& F1, const d
}
double Cft, dCftdChi;
computeConcentrationFactor(Cft, dCftdChi, Chi1);
double dres0dfVtilde = -yield/yield0*dCftdChi*pow(1.5*fVtilde*_lambda0 , 1./3.)*exp(hatP1/3.*_kappa)*_kappa/3./fVtilde*dChidChiPrime;
computeConcentrationFactor(Chi1, Cft, dCftdChi);
double dres0dfVtilde = -yield/yield0*dCftdChi*dChidtildefV;
double dres1dfVtilde = 0.;
double dres2dfVtilde = 0.;
......@@ -898,7 +922,8 @@ void mlawNonLocalPorousThomasonLaw::constitutive(
// Coalesence
double& Chi = ipvcur->getRefToLigamentRatio();
Chi = computeLigamentRatio(ipvprev->getMatrixPlasticStrain(), tildefV);
double dChidHatP, dChidtildefV;
computeLigamentRatio(ipvprev->getMatrixPlasticStrain(), tildefV, Chi, dChidHatP, dChidtildefV);
ipvcur->getRefToIPCoalescence().getRefToCoalescenceActiveFlag() = false;
......
NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.h
+ 4
− 4
View file @ b79dca56
......@@ -23,7 +23,7 @@ protected:
// Numerical parameters
double _n; // interpolation exponent
double _chi_failure; // maximal value of chi before onset of saturated regime
double _chi_switch; // maximal value of chi before onset of saturated regime
double _chi_saturated; // maximal value of chi allowed asymptotically
double _chi_decay; // rate of saturation (= characteristic lenght of exp. decay)
double _fV_Failure; // maximal value of fV allowed before total failure
......@@ -125,10 +125,10 @@ public:
};
// Specific functions to this material law
virtual double computeLigamentRatio(const double hatP, const double tildefV) const; // need to be public
virtual void computeLigamentRatio(const double hatP, const double tildefV, double& chi, double& dChidhatP, double& dChidtildefV) const; // need to be public
protected:
inline virtual void computeConcentrationFactor(double& Cft, double& dCftDChi, const double Chi) const;
inline virtual void computeConcentrationFactor(const double Chi, double& Cft, double& dCftDChi) const;
inline virtual void computeResidual(fullVector<double> &res, double& res_norm, fullMatrix<double> &J, const bool stiff,
const IPNonLocalPorosity *q0, IPNonLocalPorosity *q1,
......@@ -137,7 +137,7 @@ protected:
const double tpr_eq, const double p_pr,
const double Chin, const double fVn,
const double tildefV, const double Hatpn,
const double yield0, double& dChidChiPrime,
const double yield0, double& dChidtildefV,
double& dfVdDeltaD, double& dfVdDeltaQ, double& dfVdDeltaHatP) const;
#endif //SWIG
......
NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp
+ 2
− 1
View file @ b79dca56
......@@ -194,7 +194,8 @@ void mlawNonLocalPorousCoupledLaw::constitutive(
// Coalesence
double& Chi = ipvcur->getRefToLigamentRatio();
Chi = _mlawCoales->computeLigamentRatio(ipvprev->getMatrixPlasticStrain(), tildefV);
double dChidHatP, dChidfvtilde;
_mlawCoales->computeLigamentRatio(ipvprev->getMatrixPlasticStrain(), tildefV, Chi, dChidHatP, dChidfvtilde);
//ipvcur->getRefToIPCoalescence().getRefToCoalescenceActiveFlag() = ipvprev->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment