Skip to content
Snippets Groups Projects
Commit efbfb40b authored by Van Dung NGUYEN's avatar Van Dung NGUYEN
Browse files

new failure cri

parent 109d4f0c
No related branches found
No related tags found
1 merge request!39Master
......@@ -140,8 +140,134 @@ bool CritialPlasticDeformationCriteration::isFailed(const IPVariable* ipvprev, c
else{
double Phif = this->getFailureCriterion(ipvprev,ipv);
double r = viscoipv->getFailureOnset();
double Ff = Phif - r;
if (Ff >=0) return true;
if (Phif >= r) return true;
else return false;
}
};
void ThreeCriticalStressCriterion::getFactor(const double sigt, const double sigc, const double sigs, double& aa, double& bb, double& cc) const{
double pc = -sigc/3.;
double pt = sigt/3.;
double mat[2][2];
mat[0][0] = pc;
mat[0][1] = sigc-sigs;
mat[1][0] = pt;
mat[1][1] = sigt-sigs;
double det = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
double inv[2][2];
if(fabs(det) > 0){
double ud = 1. / det;
inv[0][0] = mat[1][1] * ud;
inv[1][0] = -mat[1][0] * ud;
inv[0][1] = -mat[0][1] * ud;
inv[1][1] = mat[0][0] * ud;
}
else{
Msg::Fatal("Singular matrix 3x3 ThreeCriticalStressCriterion::ThreeCriticalStressCriterion");
}
double vec[2];
vec[0] = (sigs-sigc)*pc;
vec[1] = (sigs-sigt)*pt;
double sol[2];
sol[0] = inv[0][0]*vec[0]+inv[0][1]*vec[1];
sol[1] = inv[1][0]*vec[0]+inv[1][1]*vec[1];
// compute factor a, b, c
aa = sol[0]*sol[1];
bb = sol[1];
cc = sigs - aa/bb;
//printf(" aa = %e bb = %e cc = %e\n",aa,bb,cc);
};
ThreeCriticalStressCriterion::ThreeCriticalStressCriterion(const int num, const RateFailureLaw& trac, const RateFailureLaw& comp, const RateFailureLaw& shear): FailureCriterionBase(num){
_compFailure = comp.clone();
_tracFailure = trac.clone();
_shearFailure = shear.clone();
_isConstant = (comp.isConstant() && trac.isConstant() && shear.isConstant());
// esitmated a,b, c
if (_isConstant){
double sigc,sigt,sigs,dR;
_compFailure->get(0,sigc,dR);
_tracFailure->get(0,sigt,dR);
_shearFailure->get(0,sigs,dR);
getFactor(sigt,sigc,sigs,a,b,c);
}
};
ThreeCriticalStressCriterion::ThreeCriticalStressCriterion(const ThreeCriticalStressCriterion& src): FailureCriterionBase(src), a(src.a),b(src.b),c(src.c), _isConstant(src._isConstant){
_compFailure = NULL;
if (src._compFailure){
_compFailure= src._compFailure->clone();
}
_tracFailure = NULL;
if (src._tracFailure){
_tracFailure = src._tracFailure->clone();
}
_shearFailure = NULL;
if (src._shearFailure){
_shearFailure = src._shearFailure->clone();
}
};
ThreeCriticalStressCriterion::~ThreeCriticalStressCriterion(){
if (_compFailure){
delete _compFailure;
_compFailure = NULL;
}
if (_tracFailure){
delete _tracFailure;
_tracFailure = NULL;
}
if (_shearFailure){
delete _shearFailure;
_shearFailure = NULL;
}
}
double ThreeCriticalStressCriterion::getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const{
const IPHyperViscoElastoPlastic* viscoipv = dynamic_cast<const IPHyperViscoElastoPlastic*>(ipv);
if (viscoipv == NULL){
Msg::Error("PowerFailureCriterion::isFailed is implemented for IPHyperViscoElastoPlastic only");
return 1e10;
}
else{
if (!_isConstant){
double sigc,sigt,sigs,dR;
const IPHyperViscoElastoPlastic* viscoipvprev = static_cast<const IPHyperViscoElastoPlastic*>(ipvprev);
_compFailure->get(viscoipvprev->_DgammaDt,sigc,dR);
_tracFailure->get(viscoipvprev->_DgammaDt,sigt,dR);
_shearFailure->get(viscoipvprev->_DgammaDt,sigs,dR);
getFactor(sigt,sigc,sigs,a,b,c);
}
// compute failure r and DrDF
const STensor3& K = viscoipv->_kirchhoff;
double pres = K.trace()/3.;
STensor3 devK = K.dev();
double Keq = sqrt(1.5*devK.dotprod());
// get coeficinet
//K.print("Kir stess");
//Msg::Info(" p = %e KEq = %e, cr = %e",pres,Keq,Keq - a/(pres+b) -c);
return (Keq - a/(pres+b) -c);
}
};
bool ThreeCriticalStressCriterion::isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const{
const IPHyperViscoElastoPlastic* viscoipv = dynamic_cast<const IPHyperViscoElastoPlastic*>(ipv);
if (viscoipv == NULL){
Msg::Error("PowerFailureCriterion::isFailed is implemented for IPHyperViscoElastoPlastic only");
return false;
}
else{
double Phif = this->getFailureCriterion(ipvprev,ipv);
double r = viscoipv->getFailureOnset();
if (Phif >= r) return true;
else return false;
}
};
\ No newline at end of file
......@@ -97,4 +97,29 @@ class CritialPlasticDeformationCriteration : public FailureCriterionBase{
#endif //SWIG
};
class ThreeCriticalStressCriterion : public FailureCriterionBase{
protected:
#ifndef SWIG
RateFailureLaw* _compFailure;
RateFailureLaw* _tracFailure;
RateFailureLaw* _shearFailure;
bool _isConstant;
mutable double a, b, c; // cache for surface sigVM = c + a/(p+b) , a,b,c estimated from sigc, sigt, sigs
private:
void getFactor(const double sigt, const double sigc, const double sigs, double& aa, double& bb, double& cc) const;
#endif //SWIG
public:
ThreeCriticalStressCriterion(const int num, const RateFailureLaw& trac, const RateFailureLaw& comp, const RateFailureLaw& shear);
#ifndef SWIG
ThreeCriticalStressCriterion(const ThreeCriticalStressCriterion& src);
virtual ~ThreeCriticalStressCriterion();
virtual double getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const;
virtual bool isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const;
virtual FailureCriterionBase* clone() const {return new ThreeCriticalStressCriterion(*this);};
#endif //SWIG
};
#endif // FAILURECRITERION_H_
......@@ -26,6 +26,7 @@ class RateFailureLaw {
virtual double getRefStrength() const = 0;
virtual void get(const double rate, double& R, double& dR) const = 0;
virtual RateFailureLaw* clone() const = 0;
virtual bool isConstant() const = 0;
#endif // SWIG
};
......@@ -51,7 +52,7 @@ class ConstantFailureLaw : public RateFailureLaw{
dR = 0.;
};
virtual RateFailureLaw* clone() const {return new ConstantFailureLaw(*this);};
virtual bool isConstant() const {return true;};
#endif // SWIG
};
......@@ -69,6 +70,7 @@ class PowerFailureLaw : public RateFailureLaw{
virtual double getRefStrength() const{return _X0;}
virtual void get(const double rate, double& R, double& dR) const ;
virtual RateFailureLaw* clone() const {return new PowerFailureLaw(*this);};
virtual bool isConstant() const {return false;};
#endif // SWIG
};
......@@ -88,6 +90,7 @@ class LograrithmicFailureLaw : public RateFailureLaw{
virtual double getRefStrength() const{return _Xr;}
virtual void get(const double rate, double& R, double& dR) const ;
virtual RateFailureLaw* clone() const {return new LograrithmicFailureLaw(*this);};
virtual bool isConstant() const {return false;};
#endif // SWIG
};
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment