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

correct stress in checking coelascence with Hill

parent 78e2e9c9
No related branches found
No related tags found
No related merge requests found
...@@ -2370,3 +2370,293 @@ void mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::I1J2J3_getYieldNormal(STen ...@@ -2370,3 +2370,293 @@ void mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::I1J2J3_getYieldNormal(STen
} }
} }
}; };
void mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const
{
// check condition if check coalescence is carried out
if (getNumOfYieldSurfaces() > 1)
{
Msg::Error("check coalescenece has not implemented for num yield surface >1");
return;
}
// Get ipvs
const IPCoalescence* q0Thom = &q0->getConstRefToIPCoalescence();
IPCoalescence* q1Thom = &q1->getRefToIPCoalescence();
// check active
if (q0Thom->getCoalescenceOnsetFlag() and q1->dissipationIsActive())
{
q1Thom->getRefToCoalescenceActiveFlag() = true;
}
else
{
q1Thom->getRefToCoalescenceActiveFlag() = false;
}
bool willUseNormalToCheck = false;
if (_checkCoalescenceWithNormal)
{
if (q1->getLocation() == IPVariable::INTERFACE_MINUS or q1->getLocation() == IPVariable::INTERFACE_PLUS)
{
willUseNormalToCheck = true;
const SVector3& normal = q1->getConstRefToCurrentOutwardNormal();
if (normal.norm() <=0.)
{
Msg::Error("coalescence cannot check with zero normal!!!");
return;
};
}
else
{
willUseNormalToCheck = false;
if (!_withBothYieldSurfacesInBulk)
{
return;
}
}
}
if (q0Thom->getCoalescenceOnsetFlag())
{
// if onset already occurs
q1Thom->getRefToCoalescenceOnsetFlag() = q0Thom->getCoalescenceOnsetFlag();
q1Thom->getRefToCoalescenceMode() = q0Thom->getCoalescenceMode();
// copy data
q1Thom->setIPvAtCoalescenceOnset(*q0Thom->getIPvAtCoalescenceOnset());
q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft();
q1Thom->getRefToLodeParameterOffset() = q0Thom->getLodeParameterOffset();
q1Thom->getRefToYieldOffset() = q0Thom->getYieldOffset();
}
else if (q1->dissipationIsActive())
{
const STensor3& kCor = q1->getConstRefToCorotationalKirchhoffStress();
const STensor3& F = q1->getConstRefToDeformationGradient();
const double R = q1->getCurrentViscoplasticYieldStress();
double J = STensorOperation::determinantSTensor3(F);
static STensor3 corSig;
corSig = kCor;
if (_stressFormulation == CORO_CAUCHY)
{
corSig *= (1./J);
}
// check at interface
bool interfaceNeckActive = true;
bool interfaceShearActive = true;
if (willUseNormalToCheck)
{
double CftThomason;
_mlawCoales->getCoalescenceLaw()->computeConcentrationFactor(q0,q1,CftThomason);
double CftShear;
_mlawShear->getCoalescenceLaw()->computeConcentrationFactor(q0,q1,CftShear);
SVector3 normal = q1->getConstRefToCurrentOutwardNormal();
normal.normalize();
// Compute Sig from Kcor following: Sig = invFeT * sigCor * FeT = Fe*sigCor*invFe
static STensor3 invFep, FeSigCor , Fe, Sig;
STensorOperation::inverseSTensor3(q1->getConstRefToFp(), invFep);
STensorOperation::multSTensor3(F, invFep, Fe);
STensorOperation::inverseSTensor3(Fe, invFep);
STensorOperation::multSTensor3(Fe, corSig, FeSigCor);
STensorOperation::multSTensor3SecondTranspose(FeSigCor, invFep, Sig);
// traction vector
static SVector3 tractionForce;
STensorOperation::multSTensor3SVector3(Sig,normal,tractionForce);
double tauN = STensorOperation::dot(tractionForce,normal);
// Deduce tangent part
double tauTanSq = 0.;
for (int i=0; i<3; i++)
{
double temp = tractionForce(i)-tauN*normal(i);
tauTanSq += temp*temp;
}
// Compute effective stress
double tauEff = 0.;
if (tauN > 0.)
{
tauEff = sqrt(tauN*tauN + _beta*tauTanSq);
}
else
{
tauEff = sqrt(_beta*tauTanSq);
}
if (tauEff - CftThomason*R >0)
{
interfaceNeckActive = true;
}
else
{
interfaceNeckActive = false;
}
if (sqrt(3.*tauTanSq) - CftShear*R >0)
{
interfaceShearActive = true;
}
else
{
interfaceShearActive = false;
}
}
static STensor3 sigEff;
STensorOperation::multSTensor43STensor3(_aniM, corSig, sigEff);
double Gurson = _mlawGrowth->I1J2J3_yieldFunction(sigEff,R,q0,q1,T);
double Thomason = _mlawCoales->I1J2J3_yieldFunction(sigEff,R,q0,q1,T);
double Shear = _mlawShear->I1J2J3_yieldFunction(sigEff,R,q0,q1,T);
// active plasticity
// check coalescence is active
bool& CoalesFlag = q1Thom->getRefToCoalescenceOnsetFlag();
int& CoalesMode = q1Thom->getRefToCoalescenceMode();
CoalesFlag = false;
CoalesMode = 0;
if (_useTwoYieldRegularization)
{
double FG = Gurson +1.;
double FT = Thomason +1.;
double FS = Shear +1.;
if (FG < 0. or FT < 0. or FS < 0)
{
Msg::Error("FT,FS FG must be all positive but FT=%e, FS = %e, FG=%e",FT,FS,FG);
}
// check necking
if (std::max(FT,FS)> FG)
{
double ratio = FG/FT;
if (pow(ratio,_yieldRegExponent) < _coalesTol && interfaceNeckActive)
{
CoalesFlag = true;
CoalesMode = 1;
printf("Necking coalescence at Gurson = %e, Thomason =%e, Shear = %e\n",
Gurson,Thomason,Shear);
}
else
{
ratio = FG/FS;
if (pow(ratio,_yieldRegExponent) < _coalesTol && interfaceShearActive)
{
CoalesFlag = true;
CoalesMode = 2;
printf("Shear coalescence at Gurson = %e, Thomason =%e, Shear = %e\n",
Gurson,Thomason,Shear);
}
else
{
CoalesFlag = false;
CoalesMode = 0;
}
}
}
}
else
{
if (Thomason > std::max(Gurson,Shear) && interfaceNeckActive)
{
CoalesFlag = true;
CoalesMode = 1;
printf("Necking coalescence at Gurson = %e, Thomason =%e, Shear = %e\n",
Gurson,Thomason,Shear);
}
else if (Shear > std::max(Gurson,Thomason) && interfaceShearActive)
{
CoalesFlag = true;
CoalesMode = 2;
printf("Shear coalescence at Gurson = %e, Thomason =%e, Shear = %e\n",
Gurson,Thomason,Shear);
}
else
{
CoalesFlag = false;
CoalesMode = 0;
}
}
// if coalesnce occurs
if (CoalesFlag)
{
q1Thom->setIPvAtCoalescenceOnset(*q1);
// if offset is used to satisfy Thomason yield surface at the onset of coalescence
if (_useTwoYieldRegularization)
{
// never offset with yield regulation
q1Thom->getRefToCrackOffsetOnCft() = 1.;
q1Thom->getRefToLodeParameterOffset() = 1.;
q1Thom->getRefToYieldOffset() = 0.;
}
else
{
if (_withCftOffset)
{
if (CoalesMode==1)
{
q1Thom->getRefToCrackOffsetOnCft() = 1. + Thomason;
}
else if (CoalesMode==2)
{
q1Thom->getRefToCrackOffsetOnCft() = 1. + Shear;
}
q1Thom->getRefToLodeParameterOffset() = 1.;
q1Thom->getRefToYieldOffset() = 0.;
#ifdef _DEBUG
printf("coalescence ocurrs using Cft offset = %e\n",q1Thom->getCrackOffsetOnCft());
#endif //_DEBUG
}
else if (_withLodeOffset)
{
// Lode offset is not used for material law, because law already accounts for Lode variable
q1Thom->getRefToCrackOffsetOnCft() = 1.;
q1Thom->getRefToLodeParameterOffset()= 1.;
q1Thom->getRefToYieldOffset() = 0.;
}
else if (_withYieldOffset)
{
q1Thom->getRefToCrackOffsetOnCft() = 1.;
q1Thom->getRefToLodeParameterOffset() = 1.;
if (CoalesMode==1)
{
q1Thom->getRefToYieldOffset() = Thomason;
}
else if (CoalesMode==2)
{
q1Thom->getRefToYieldOffset() = Shear;
}
#ifdef _DEBUG
printf("coalescence ocurrs using yield offset = %e\n",q1Thom->getYieldOffset());
#endif //_DEBUG
}
else
{
q1Thom->getRefToCrackOffsetOnCft() = 1.;
q1Thom->getRefToLodeParameterOffset() = 1.;
q1Thom->getRefToYieldOffset() = 0.;
#ifdef _DEBUG
printf("coalescence ocurrs, not offset is used\n");
#endif //_DEBUG
}
}
}
else
{
q1Thom->getRefToCrackOffsetOnCft() = 1.;
q1Thom->getRefToLodeParameterOffset() = 1.;
q1Thom->getRefToYieldOffset() = 0.;
}
}
};
void mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::forceCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const
{
Msg::Error("mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::forceCoalescence is not implemented!!!");
};
bool mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::checkCrackCriterion(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const STensor3& Fcur, const SVector3& normal, double delayFactor, const double* T) const
{
Msg::Error("mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::checkCrackCriterion is not implemented!!!");
}
\ No newline at end of file
...@@ -311,6 +311,9 @@ class mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS : public mlawNonLocalPorou ...@@ -311,6 +311,9 @@ class mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS : public mlawNonLocalPorou
bool diff=false, STensor43* DNpDkcor=NULL, STensor3* DNpDR=NULL, std::vector<STensor3>* DNpDY=NULL, bool diff=false, STensor43* DNpDkcor=NULL, STensor3* DNpDR=NULL, std::vector<STensor3>* DNpDY=NULL,
bool withThermic=false, STensor3* DNpDT=NULL bool withThermic=false, STensor3* DNpDT=NULL
) const; ) const;
virtual void checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const;
virtual void forceCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const;
virtual bool checkCrackCriterion(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const STensor3& Fcur, const SVector3& normal, double delayFactor, const double* T) const;
}; };
#endif // MLAWNONLOCALPOROUSCOALESCENCE_H_ #endif // MLAWNONLOCALPOROUSCOALESCENCE_H_
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment