Select Git revision
-
Giannis Nikiteas authoredGiannis Nikiteas authored
mlawNonLocalPorousCoupled.cpp 35.93 KiB
//
// C++ Interface: material law
//
// Description: Porous material law with Gurson void growth and coalescence
//
// Author: <J. Leclerc>, (C) 2017
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "mlawNonLocalPorousCoupled.h"
#include "nonLinearMechSolver.h"
mlawNonLocalPorousCoupledLaw::mlawNonLocalPorousCoupledLaw(const int num,const double E,const double nu, const double rho,
const double q1, const double q2, const double q3,
const double fVinitial, const double lambda0, const J2IsotropicHardening &j2IH, const CLengthLaw &cLLaw,
const double tol, const bool matrixbyPerturbation, const double pert) :
mlawNonLocalPorosity(num, E, nu, rho, fVinitial, j2IH, cLLaw, tol, matrixbyPerturbation, pert),
_withCrackTransition(false)
{
_mlawGrowth = new mlawNonLocalDamageGurson(num, E, nu, rho, q1, q2, q3, fVinitial, j2IH, cLLaw, tol, matrixbyPerturbation, pert);
_mlawCoales = new mlawNonLocalPorousThomasonLaw(num, E, nu, rho, fVinitial, lambda0, j2IH, cLLaw, tol, matrixbyPerturbation, pert);
// modify current coalescence law
mlawNonLocalPorosity::setCoalescenceLaw(*_mlawCoales->getCoalescenceLaw());
_failedTol = 0.99;
};
mlawNonLocalPorousCoupledLaw::mlawNonLocalPorousCoupledLaw(const int num,const double E,const double nu, const double rho,
const double q1, const double q2, const double q3,
const double fVinitial, const double lambda0, const double kappa, const J2IsotropicHardening &j2IH, const CLengthLaw &cLLaw,
const double tol, const bool matrixbyPerturbation, const double pert) :
mlawNonLocalPorosity(num, E, nu, rho, fVinitial, j2IH, cLLaw, tol, matrixbyPerturbation, pert),
_withCrackTransition(false)
{
_mlawGrowth = new mlawNonLocalDamageGurson(num, E, nu, rho, q1, q2, q3, fVinitial, j2IH, cLLaw, tol, matrixbyPerturbation, pert);
_mlawCoales = new mlawNonLocalPorousThomasonLaw(num, E, nu, rho, fVinitial, lambda0, kappa, j2IH, cLLaw, tol, matrixbyPerturbation, pert);
// modify current coalescence law
mlawNonLocalPorosity::setCoalescenceLaw(*_mlawCoales->getCoalescenceLaw());
_failedTol = 0.99;
};
mlawNonLocalPorousCoupledLaw::mlawNonLocalPorousCoupledLaw(const mlawNonLocalPorousCoupledLaw &source) :
mlawNonLocalPorosity(source),_withCrackTransition(source._withCrackTransition)
{
_mlawGrowth = NULL;
if(source._mlawGrowth != NULL) _mlawGrowth = dynamic_cast<mlawNonLocalDamageGurson*>(source._mlawGrowth->clone());
_mlawCoales = NULL;
if(source._mlawCoales != NULL) _mlawCoales = dynamic_cast<mlawNonLocalPorousThomasonLaw*>(source._mlawCoales->clone());
};
void mlawNonLocalPorousCoupledLaw::setNonLocalMethod(const int i){
mlawNonLocalPorosity::setNonLocalMethod(i);
_mlawGrowth->setNonLocalMethod(i);
_mlawCoales->setNonLocalMethod(i);
};
// Options settings
void mlawNonLocalPorousCoupledLaw::setOrderForLogExp(const int newOrder)
{
mlawNonLocalPorosity::setOrderForLogExp(newOrder); _mlawGrowth->setOrderForLogExp(newOrder);
_mlawCoales->setOrderForLogExp(newOrder);
Msg::Info("mlawNonLocalPorousCoupledLaw :: order = %d is used to approximate log and exp",newOrder);
};
void mlawNonLocalPorousCoupledLaw::setViscosityLaw(const viscosityLaw& vico)
{
mlawNonLocalPorosity::setViscosityLaw(vico);
_mlawGrowth->setViscosityLaw(vico);
_mlawCoales->setViscosityLaw(vico);
};
void mlawNonLocalPorousCoupledLaw::setNucleationLaw(const NucleationLaw& added_gdnLaw)
{
mlawNonLocalPorosity::setNucleationLaw(added_gdnLaw);
_mlawGrowth->setNucleationLaw(added_gdnLaw);
_mlawCoales->setNucleationLaw(added_gdnLaw);
};
void mlawNonLocalPorousCoupledLaw::setScatterredInitialPorosity(double f0min, double f0max)
{
mlawNonLocalPorosity::setScatterredInitialPorosity(f0min,f0max);
_mlawGrowth->setScatterredInitialPorosity(f0min, f0max);
_mlawCoales->setScatterredInitialPorosity(f0min, f0max);
Msg::Info("mlawNonLocalPorousCoupledLaw :: Add scatter in initial porosity");
_randMinbound = f0min;
if(f0min < 0.){Msg::Fatal("mlawNonLocalPorousCoupledLaw :: Wrong scatter parameter fV0_min");}
_randMaxbound = f0max;
if(f0max < f0min){Msg::Fatal("mlawNonLocalPorousCoupledLaw :: Wrong scatter parameter fV0_max");}
};
void mlawNonLocalPorousCoupledLaw::setCoalescenceLaw(const CoalescenceLaw& added_coalsLaw)
{
mlawNonLocalPorosity::setCoalescenceLaw(added_coalsLaw);
_mlawCoales->setCoalescenceLaw(added_coalsLaw);
};
void mlawNonLocalPorousCoupledLaw::setYieldSurfaceExponent(const double newN)
{
_mlawCoales->setYieldSurfaceExponent(newN);
};
void mlawNonLocalPorousCoupledLaw::setRoundedYieldSurfaceMethod(const int method){
_mlawCoales->setRoundedYieldSurfaceMethod(method);
};
void mlawNonLocalPorousCoupledLaw::setOnsetTriaxialityForRoundedYieldSurface(const double newTc){
_mlawCoales->setOnsetTriaxialityForRoundedYieldSurface(newTc);
};
void mlawNonLocalPorousCoupledLaw::setCfTOffsetMethod(const int MethodNumber){
_CfTOffsetMethod = MethodNumber;
_mlawCoales->setCfTOffsetMethod(MethodNumber);
if (_CfTOffsetMethod == 0 or _CfTOffsetMethod == 1 or _CfTOffsetMethod == 2 or _CfTOffsetMethod == 3){
Msg::Info("mlawNonLocalPorousCoupledLaw::useCfTOffset: Offset method %d is used",_CfTOffsetMethod);
}
else{
Msg::Fatal("mlawNonLocalPorousCoupledLaw::useCfTOffset: Offset method %d is not existing",_CfTOffsetMethod);
}
};
void mlawNonLocalPorousCoupledLaw::setCrackTransition(const bool fl){
_withCrackTransition = fl;
};
void mlawNonLocalPorousCoupledLaw::setBalanceEquationType(const int type){
_mlawCoales->setBalanceEquationType(type);
};
void mlawNonLocalPorousCoupledLaw::setShearPorosityGrowthFactor(const double k) {
mlawNonLocalPorosity::setShearPorosityGrowthFactor(k);
_mlawGrowth->setShearPorosityGrowthFactor(k);
_mlawCoales->setShearPorosityGrowthFactor(k);
};
void mlawNonLocalPorousCoupledLaw::setFailureTolerance(const double NewFailureTol){
if ((NewFailureTol > 1.) or (NewFailureTol < 0.)){
Msg::Fatal("mlawNonLocalPorousCoupledLaw::setFailureTolerance: Uncorrect value of tolerance failure");
}
_failedTol = NewFailureTol;
_mlawGrowth->setFailureTolerance(NewFailureTol);
_mlawCoales->setFailureTolerance(NewFailureTol);
};
void mlawNonLocalPorousCoupledLaw::setLocalRegularizedFunction(const scalarFunction& fct){
mlawNonLocalPorosity::setLocalRegularizedFunction(fct);
_mlawGrowth->setLocalRegularizedFunction(fct);
_mlawCoales->setLocalRegularizedFunction(fct);
};
void mlawNonLocalPorousCoupledLaw::setCorrectedRegularizedFunction(const scalarFunction& fct){
mlawNonLocalPorosity::setCorrectedRegularizedFunction(fct);
_mlawGrowth->setCorrectedRegularizedFunction(fct);
_mlawCoales->setCorrectedRegularizedFunction(fct);
};
void mlawNonLocalPorousCoupledLaw::setPostBlockedDissipationBehavior(const int method){
mlawNonLocalPorosity::setPostBlockedDissipationBehavior(method);
_mlawGrowth->setPostBlockedDissipationBehavior(method);
_mlawCoales->setPostBlockedDissipationBehavior(method);
};
void mlawNonLocalPorousCoupledLaw::clearCLengthLaw(){
mlawNonLocalPorosity::clearCLengthLaw();
_mlawGrowth->clearCLengthLaw();
_mlawCoales->clearCLengthLaw();
};
void mlawNonLocalPorousCoupledLaw::setCLengthLaw(const CLengthLaw& clength){
mlawNonLocalPorosity::setCLengthLaw(clength);
_mlawGrowth->setCLengthLaw(clength);
_mlawCoales->setCLengthLaw(clength);
};
//-------------------added
void mlawNonLocalPorousCoupledLaw::setReferenceT(const double referenceT){
mlawNonLocalPorosity::setReferenceT(referenceT);
_mlawGrowth->setReferenceT(referenceT);
_mlawCoales->setReferenceT(referenceT);
};
void mlawNonLocalPorousCoupledLaw::setThermalExpansionCoefficient(const double alp){
mlawNonLocalPorosity::setThermalExpansionCoefficient(alp);
_mlawGrowth->setThermalExpansionCoefficient(alp);
_mlawCoales->setThermalExpansionCoefficient(alp);
};
void mlawNonLocalPorousCoupledLaw::setSubStepping(const bool fl, const int maxNumStep){
mlawNonLocalPorosity::setSubStepping(fl,maxNumStep);
_mlawGrowth->setSubStepping(fl,maxNumStep);
_mlawCoales->setSubStepping(fl,maxNumStep);
};
void mlawNonLocalPorousCoupledLaw::setTemperatureFunction_rho(const scalarFunction& Tfunc){
mlawNonLocalPorosity::setTemperatureFunction_rho(Tfunc);
_mlawGrowth->setTemperatureFunction_rho(Tfunc);
_mlawCoales->setTemperatureFunction_rho(Tfunc);
};
void mlawNonLocalPorousCoupledLaw::setTemperatureFunction_E(const scalarFunction& Tfunc){
mlawNonLocalPorosity::setTemperatureFunction_E(Tfunc);
_mlawGrowth->setTemperatureFunction_E(Tfunc);
_mlawCoales->setTemperatureFunction_E(Tfunc);
};
void mlawNonLocalPorousCoupledLaw::setTemperatureFunction_nu(const scalarFunction& Tfunc){ mlawNonLocalPorosity::setTemperatureFunction_nu(Tfunc);
_mlawGrowth->setTemperatureFunction_nu(Tfunc);
_mlawCoales->setTemperatureFunction_nu(Tfunc);
};
void mlawNonLocalPorousCoupledLaw::setTemperatureFunction_alp(const scalarFunction& Tfunc){
mlawNonLocalPorosity::setTemperatureFunction_alp(Tfunc);
_mlawGrowth->setTemperatureFunction_alp(Tfunc);
_mlawCoales->setTemperatureFunction_alp(Tfunc);
};
void mlawNonLocalPorousCoupledLaw::setStressTriaxialityFunction_kw(const scalarFunction& fct){
mlawNonLocalPorosity::setStressTriaxialityFunction_kw(fct);
_mlawGrowth->setStressTriaxialityFunction_kw(fct);
_mlawCoales->setStressTriaxialityFunction_kw(fct);
};
void mlawNonLocalPorousCoupledLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
//bool inter=true;
//const MInterfaceElement *iele = static_cast<const MInterfaceElement*>(ele);
//if(iele==NULL) inter=false;
double fvInit = getInitialPorosity();
IPVariable* ipvi = new IPNonLocalPorosity(fvInit,_j2IH,&_cLLaw,&_gdnLawContainer,_coalescenceLaw);
IPVariable* ipv1 = new IPNonLocalPorosity(fvInit,_j2IH,&_cLLaw,&_gdnLawContainer,_coalescenceLaw);
IPVariable* ipv2 = new IPNonLocalPorosity(fvInit,_j2IH,&_cLLaw,&_gdnLawContainer,_coalescenceLaw);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
};
double mlawNonLocalPorousCoupledLaw::yieldFunction(const double kcorEq, const double pcor, const double R, const double yieldfV,
const IPNonLocalPorosity* q0, const IPNonLocalPorosity* q1, const double* T) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
return _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
}
else return _mlawGrowth->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
};
void mlawNonLocalPorousCoupledLaw::computeResidual(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
const double yieldfV, const double DyieldfVDfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
const double DeltafV, const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
const double* T
) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
_mlawCoales->computeResidual(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,
DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
}
else {
_mlawGrowth->computeResidual(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,
DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
}
}
void mlawNonLocalPorousCoupledLaw::computeDResidual(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
const double yieldfV, const double DyieldfVDfV, const double DyieldfVDtildefV,
const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr, const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T, const double* DRDT,
const STensor3* DdevKcorprDT , const double* DpcorprDT, const double* DDeltafVDT,
double *dres0dT, double *dres1dT, double *dres2dT
) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
_mlawCoales->computeDResidual(dres0dEpr,dres1dEpr,dres2dEpr,
dres0dtildefV,dres1dtildefV,dres2dtildefV,
devKcorpr,kcorEqpr,pcorpr,R,
yieldfV,DyieldfVDfV,DyieldfVDtildefV,
DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
q0, q1,
DeltaHatD,DeltaHatQ,DeltaHatP,
T,DRDT,
DdevKcorprDT ,DpcorprDT,DDeltafVDT,
dres0dT,dres1dT,dres2dT);
}
else {
_mlawGrowth->computeDResidual(dres0dEpr,dres1dEpr,dres2dEpr,
dres0dtildefV,dres1dtildefV,dres2dtildefV,
devKcorpr,kcorEqpr,pcorpr,R,
yieldfV,DyieldfVDfV,DyieldfVDtildefV,
DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
q0, q1,
DeltaHatD,DeltaHatQ,DeltaHatP,
T,DRDT,
DdevKcorprDT ,DpcorprDT,DDeltafVDT,
dres0dT,dres1dT,dres2dT);
}
};
void mlawNonLocalPorousCoupledLaw::computeDResidualLocal(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
const double yieldfV, const double DyieldfVDfV,
const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T, const double* DRDT,
const STensor3* DdevKcorprDT , const double* DpcorprDT, const double* DDeltafVDT,
double *dres0dT, double *dres1dT, double *dres2dT
) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
_mlawCoales->computeDResidualLocal(dres0dEpr, dres1dEpr, dres2dEpr,
devKcorpr, kcorEqpr, pcorpr,R,
yieldfV, DyieldfVDfV,
DDeltafVDEpr, DdevKcorprDEpr, DpcorprDEpr,
q0, q1,
DeltaHatD, DeltaHatQ,DeltaHatP,
T, DRDT,
DdevKcorprDT , DpcorprDT,DDeltafVDT,
dres0dT,dres1dT,dres2dT);
}
else{
_mlawGrowth->computeDResidualLocal(dres0dEpr, dres1dEpr, dres2dEpr,
devKcorpr, kcorEqpr, pcorpr,R,
yieldfV, DyieldfVDfV,
DDeltafVDEpr, DdevKcorprDEpr, DpcorprDEpr,
q0, q1,
DeltaHatD, DeltaHatQ,DeltaHatP,
T, DRDT,
DdevKcorprDT , DpcorprDT,DDeltafVDT,
dres0dT,dres1dT,dres2dT); }
}
double mlawNonLocalPorousCoupledLaw::getOnsetCriterion(const IPNonLocalPorosity* q1) const{
return q1->getConstRefToIPCoalescence().getCrackOffsetOnCft();
}
void mlawNonLocalPorousCoupledLaw::nucleation(const double p1, const double p0, const double f0, IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
_mlawCoales->nucleation(p1,p0,f0,q1,q0);
}
else{
_mlawGrowth->nucleation(p1,p0,f0,q1,q0);
}
};
void mlawNonLocalPorousCoupledLaw::checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const
{
// Get ipvs
const IPThomasonCoalescence* q0Thom = dynamic_cast<const IPThomasonCoalescence*>(&q0->getConstRefToIPCoalescence());
IPThomasonCoalescence* q1Thom = dynamic_cast<IPThomasonCoalescence*>(&q1->getRefToIPCoalescence());
if (q1Thom == NULL){
Msg::Fatal("mlawNonLocalPorousThomasonLaw::computeResidual: IPThomasonCoalescence must be used!");
}
// Get stress state
const STensor3& kCor = q1->getConstRefToCorotationalKirchhoffStress();
static STensor3 devKcor;
static double kcorEq, pcor;
STensorOperation::decomposeDevTr(kCor,devKcor,pcor);
pcor /= 3.;
kcorEq = sqrt(1.5*devKcor.dotprod());
// Get needed data for testing criterion
const double R = q1->getCurrentViscoplasticYieldStress();
const double yieldfV = q1->getYieldPorosity();
const double R0 = _j2IH->getYield0();
const double Cft = q1Thom->getConcentrationFactor();
// Coupling management
if (_withCrackTransition){
// If crack transition is used,
// At the interface: coalescence is used, namely as crack insertion criterion
// In the bulk, nothing is done (== always Gurson)
q1Thom->getRefToYieldOffset() = 0.;
if (q1->getLocation() == IPVariable::INTERFACE_MINUS or q1->getLocation() == IPVariable::INTERFACE_PLUS){
// At the interface:
// Check if transition has already occurs or not
if (!q1Thom->getCoalescenceOnsetFlag()){
// If Coalescence has not yet occured at least once
// Check if the normal traction force is sufficient to obtain coalescence in normal direction of the interface
// Get normal traction on the interface from tau stress
double sigmax = 0.;
static SVector3 Norm;
Norm = q1->getConstRefToReferenceOutwardNormal();
if (Norm.norm() > 0){ // is it usefulll ???
Norm.normalize();
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
sigmax += Norm(i)*kCor(i,j)*Norm(j);
}
}
}
else{
Msg::Fatal("mlawNonLocalPorousCoupledLaw::updatePlasticState: zero normal");
}
// Check if the normal traction force is sufficient to obtain coalescence in normal direction of the interface
// and manage offset
// Why is the offset needed?
// We want to ensure a continuity with the Thomason yield surface at coalescence onset
// The offset correspond to a factor to obtain fT = 0 with the actual stress state
// -actual stress state: yieldThomason = sthingOfStress - Cft*R/R0
// -correction to satidfy fT=0 with actual stress: yieldThomason + Cft*R/R0 - CftOffset*Cft*R/R0 = 0
// -solving for offset gives: CftOffset = (yieldThomason + Cft*R/R0)/(Cft*R/R0) = 1+yieldThomason*R0/(Cft*R)
// NB: if _CfTOffsetMethod > 1: fT>0,
// Compute Thomason yield surface
double yieldThomason = _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
// Offset management
if(_CfTOffsetMethod == 2 or _CfTOffsetMethod == 3){
// Method 2: Offset is always computed to ensure contuinity if coalescence is forced (offset can be <1)
q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
}
if ((sigmax - Cft*R > _tol) and (yieldThomason > _tol)){
// If coalescence occurs, update Cft offset value and others variables
// Update onset coalescence variables
// Activate coalescence onset
q1Thom->getRefToCoalescenceOnsetFlag() = true;
// Set onset variables to the actual values (yield porosity and geometrical parameters)
q1Thom->getRefToPorosityAtCoalescenceOnset() = q1->getYieldPorosity();
q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q1Thom->getLigamentRatio();
q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q1Thom->getAspectRatio();
q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q1Thom->getShapeFactor();
q1Thom->getRefToCoalescenceActiveFlag() = true;
q1Thom->getRefToAccelerateRate() = 1.;
// Offset management
if (_CfTOffsetMethod == 1){
// Method 1: offset computed at coalescence onset only
q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
Msg::Info("Coalescence occurs for the first time: Offset =%e",q1Thom->getRefToCrackOffsetOnCft());
}
else if (_CfTOffsetMethod == 0){
// Method 0: Offset is not used: always equal to 1.
q1Thom->getRefToCrackOffsetOnCft() = 1.;
}
}
else{
// Coalescence never occurs and is not occuring at the next step
q1Thom->getRefToCoalescenceOnsetFlag() = false;
// Reset onset variables to their default value (yield porosity and geometrical parameters)
q1Thom->getRefToPorosityAtCoalescenceOnset() = 0.;
q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = 0.;
q1Thom->getRefToAspectRatioAtCoalescenceOnset() = 0.;
q1Thom->getRefToShapeFactorAtCoalescenceOnset() = 0.;
q1Thom->getRefToCoalescenceActiveFlag() = false;
q1Thom->getRefToAccelerateRate() = 1.;
// Offset management
if (_CfTOffsetMethod == 1 or _CfTOffsetMethod == 0){
q1Thom->getRefToCrackOffsetOnCft() = 1.;
}
}
}
else{
// Coalescence has occured at least once
q1Thom->getRefToCoalescenceOnsetFlag() = true; // already done....
// Update onset value (yield porosity and geometrical parameters) q1Thom->getRefToPorosityAtCoalescenceOnset() = q0Thom->getPorosityAtCoalescenceOnset();
q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q0Thom->getLigamentRatioAtCoalescenceOnset();
q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q0Thom->getAspectRatioAtCoalescenceOnset();
q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q0Thom->getShapeFactorAtCoalescenceOnset();
// Cft onset - Method 0-1-2: keep it always constant once the coalescence occurs
q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft();
// Determine which mode should be used for the next time step.
if(!q0Thom->getCoalescenceOnsetFlag()){
// First time that coalescence is detected: Thomason == ON
// This case is only appearing with with selective update
q1Thom->getRefToCoalescenceActiveFlag() = true;
q1Thom->getRefToAccelerateRate() = 1.0;
q1Thom->getRefToYieldOffset() = _mlawGrowth->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
}
else{
// Not the first time that coalescence appears
if (q0Thom->getCoalescenceActiveFlag()){
// If Thomason is/was used for this time step : check if Gurson is more restrictive
double yieldGurson = _mlawGrowth->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
q1Thom->getRefToYieldOffset() = yieldGurson;
if(yieldGurson > _tol){
Msg::Info("Change yield surface to Gurson fG =%e",yieldGurson);
q1Thom->getRefToCoalescenceActiveFlag() = true; // Gurson will be used at the next time step
q1Thom->getRefToAccelerateRate() = 1.0;
}
else{
q1Thom->getRefToCoalescenceActiveFlag() = true; // Thomason will be used at the next time step
q1Thom->getRefToAccelerateRate() = 1.0;
}
}
else{
// If Gurson was/is used for this time step : check if Thomason is more restrictive
double yieldThomason = _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
q1Thom->getRefToYieldOffset() = yieldThomason;
if(yieldThomason > _tol){
Msg::Info("Change yield surface to Thomason fT = %e",yieldThomason);
q1Thom->getRefToCoalescenceActiveFlag() = true; // Thomason will be used at the next time step
q1Thom->getRefToAccelerateRate() = 1.0;
if (_CfTOffsetMethod == 3){
q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft() + q1Thom->getYieldOffset()*R0/(Cft*R);
if (q1Thom->getRefToCrackOffsetOnCft() > 2.){Msg::Error("High predicted value of Cft Offset = %e",q1Thom->getRefToCrackOffsetOnCft());};
}
}
else{
q1Thom->getRefToCoalescenceActiveFlag() = false; // Gurson will be used at the next time step
q1Thom->getRefToAccelerateRate() = 1.0;
}
}
}
}
}
else
{
// In the bulk:
// Coalescence never occurs
q1Thom->getRefToCoalescenceOnsetFlag() = false;
// Reset onset variables to their default value
// Porosity
q1Thom->getRefToPorosityAtCoalescenceOnset() = 0.;
// Geometrical parameters
q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = 0.;
q1Thom->getRefToAspectRatioAtCoalescenceOnset() = 0.;
q1Thom->getRefToShapeFactorAtCoalescenceOnset() = 0.;
q1Thom->getRefToCoalescenceActiveFlag() = false; q1Thom->getRefToCrackOffsetOnCft() = 1.;
q1Thom->getRefToAccelerateRate() = 1.;
q1Thom->getRefToYieldOffset() = 0.;
}
}
else{
// If crack transition is not used: no distinction between bulk and interface ipvs
// if Thomason yield surface is first used, it always is used
// onset and active are the same
if (!q0Thom->getCoalescenceOnsetFlag()){
double yieldThomason = _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
if (yieldThomason > _tol){
q1Thom->getRefToCoalescenceCriterion() = 0.;
// update offset
if (_CfTOffsetMethod == 1){
q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
Msg::Info("coalescence occurs, use Thomason yield surface in next step with offset = %e",q1Thom->getCrackOffsetOnCft());
}
else if(_CfTOffsetMethod == 0){
q1Thom->getRefToCrackOffsetOnCft() = 1.;
}
q1Thom->getRefToCoalescenceActiveFlag() = true;
q1Thom->getRefToCoalescenceOnsetFlag() = true;
//Msg::Info("coalescene yield surface is used");
// initial coalesence onset value
q1Thom->getRefToPorosityAtCoalescenceOnset() = yieldfV;
q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q1Thom->getLigamentRatio();
q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q1Thom->getAspectRatio();
q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q1Thom->getShapeFactor();
q1Thom->getRefToAccelerateRate() = 1.;
}
else{
q1Thom->getRefToCoalescenceCriterion() = yieldThomason; // save data
q1Thom->getRefToCoalescenceActiveFlag() = false;
q1Thom->getRefToCoalescenceOnsetFlag() = false;
// reset to defaut value
q1Thom->getRefToPorosityAtCoalescenceOnset() = 0.;
q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = 0.;
q1Thom->getRefToAspectRatioAtCoalescenceOnset() = 0.;
q1Thom->getRefToShapeFactorAtCoalescenceOnset() = 0.;
q1Thom->getRefToCrackOffsetOnCft() = 1.;
q1Thom->getRefToAccelerateRate() = 1.;
}
}
else{
q1Thom->getRefToCoalescenceCriterion() = 0.;
q1Thom->getRefToCoalescenceActiveFlag() = true;
q1Thom->getRefToCoalescenceOnsetFlag() = true;
q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q0Thom->getLigamentRatioAtCoalescenceOnset();
q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q0Thom->getAspectRatioAtCoalescenceOnset();
q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q0Thom->getShapeFactorAtCoalescenceOnset();
q1Thom->getRefToPorosityAtCoalescenceOnset() = q0Thom->getPorosityAtCoalescenceOnset();
q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft();
q1Thom->getRefToAccelerateRate() = q0Thom->getAccelerateRate();
}
}
};
void mlawNonLocalPorousCoupledLaw::checkFailed(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const{
// Get ipvs
const IPThomasonCoalescence* q0Thom = dynamic_cast<const IPThomasonCoalescence*>(&q0->getConstRefToIPCoalescence());
IPThomasonCoalescence* q1Thom = dynamic_cast<IPThomasonCoalescence*>(&q1->getRefToIPCoalescence());
if (q1Thom == NULL){ Msg::Fatal("mlawNonLocalPorousCoupledLaw::checkFailed: IPThomasonCoalescence must be used!");
}
// Manage failure flag:
// If it's previously failed, it always be failed. Otherwise, the criterion will depends on the active plastic mode
if (q0->isFailed()){
q1->setFailed(true);
}
else
{
q1->setFailed(false);
// If coalescence plastic mode is currently used: criterion based on Cft
double yieldfV = q1->getYieldPorosity();
double Chi = q1Thom->getLigamentRatio();
double Cft = q1Thom->getConcentrationFactor();
double CftOffest = q1Thom->getCrackOffsetOnCft();
if (_CfTOffsetMethod==2 or _CfTOffsetMethod==3){
if (!q0->getConstRefToIPCoalescence().getCoalescenceOnsetFlag()){CftOffest = 1.;};
}
if (Cft*CftOffest < (1.-_failedTol)*5.){
q1->setFailed(true);
Msg::Info("mlawNonLocalPorousCoupledLaw::checkFailed: an ip failed with Cft at Cft= %e, chi= %e, yieldfV= %e", Cft, Chi, yieldfV);
}
else if (Chi > _failedTol){
q1->setFailed(true);
Msg::Info("mlawNonLocalPorousCoupledLaw::checkFailed: an ip failed with Chi at Cft= %e, chi= %e, yieldfV= %e", Cft, Chi, yieldfV);
}
else if (yieldfV > _failedTol*_mlawGrowth->getFailureYieldPorosityValue() ){
q1->setFailed(true);
Msg::Info("mlawNonLocalPorousCoupledLaw::checkFailed: an ip failed with fVy at Cft= %e, chi= %e, yieldfV= %e", Cft, Chi, yieldfV);
}
}
};
void mlawNonLocalPorousCoupledLaw::localPorosityGrowth(double& DeltafV,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity *q1) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
_mlawCoales->localPorosityGrowth(DeltafV,DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,q0,q1);
}
else{
_mlawGrowth->localPorosityGrowth(DeltafV,DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,q0,q1);
}
};
// Local porosity growth (with der.)
void mlawNonLocalPorousCoupledLaw::localPorosityGrowth(double& DeltafV, double& DDeltafVDDeltaHatD,double& DDeltafVDDeltaHatQ, double& DDeltafVDDeltaHatP, STensor3& DDeltafVDEpr,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const STensor43& DdevKcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity *q1,
const double* T, const STensor3* DdevKcorprDT, double* DDeltafVDT) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
_mlawCoales->localPorosityGrowth(DeltafV,DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,DDeltafVDEpr,
DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,DdevKcorprDEpr,q0,q1,
T,DdevKcorprDT,DDeltafVDT);
}
else{
_mlawGrowth->localPorosityGrowth(DeltafV,DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,DDeltafVDEpr,
DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,DdevKcorprDEpr,q0,q1,
T,DdevKcorprDT,DDeltafVDT);
}
};
void mlawNonLocalPorousCoupledLaw::computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J, const double kcorEqpr, const double pcorpr, const double R, const double H,
const double yieldfV, const double DyieldfVDfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
const double DeltafV, const double* T
) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
_mlawCoales->computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
}
else {
_mlawGrowth->computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
}
};
void mlawNonLocalPorousCoupledLaw::computeDResidual_multipleNonLocalVariables(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
std::vector<double> &dres0dtildefV, std::vector<double> &dres1dtildefV, std::vector<double> &dres2dtildefV,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
const double yieldfV, const double DyieldfVDfV, const std::vector<double>& DyieldfVDtildefV,
const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T, const double* DRDT,
const STensor3* DdevKcorprDT, const double* DpcorprDT, const double* DDeltafVDT,
double *dres0dT, double *dres1dT, double *dres2dT
) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
_mlawCoales->computeDResidual_multipleNonLocalVariables(dres0dEpr,dres1dEpr,dres2dEpr,
dres0dtildefV,dres1dtildefV,dres2dtildefV,
devKcorpr,kcorEqpr,pcorpr,R,
yieldfV,DyieldfVDfV,DyieldfVDtildefV,
DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
q0, q1,
DeltaHatD,DeltaHatQ,DeltaHatP,
T,DRDT,
DdevKcorprDT ,DpcorprDT,DDeltafVDT,
dres0dT,dres1dT,dres2dT);
}
else {
_mlawGrowth->computeDResidual_multipleNonLocalVariables(dres0dEpr,dres1dEpr,dres2dEpr,
dres0dtildefV,dres1dtildefV,dres2dtildefV,
devKcorpr,kcorEqpr,pcorpr,R,
yieldfV,DyieldfVDfV,DyieldfVDtildefV,
DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
q0, q1,
DeltaHatD,DeltaHatQ,DeltaHatP,
T,DRDT,
DdevKcorprDT ,DpcorprDT,DDeltafVDT,
dres0dT,dres1dT,dres2dT);
}
};