From 092304abaef1236ca19c5e98e3d2af49861710bf Mon Sep 17 00:00:00 2001
From: ludovic noels <l.noels@ulg.ac.be>
Date: Sun, 4 Jun 2023 18:55:56 +0200
Subject: [PATCH] add energy integral
---
.../materialLaw/mlawHyperelastic.cpp | 51 +++++++++++++------
.../materialLaw/mlawHyperelastic.h | 4 +-
2 files changed, 37 insertions(+), 18 deletions(-)
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
index 50a6486a3..ec66346a5 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
@@ -98,7 +98,7 @@ void mlawHyperViscoElastic::setViscoElasticData(const std::string filename){
fclose(file);
};
-void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &B_d, STensor3 &dB_vddev) const
+void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &intA, double &B_d, STensor3 &dB_vddev, double &intB) const
{
int method =1;
if(method ==0)
@@ -114,18 +114,37 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
dB_vddev=dev;
dB_vddev*=2.*getPiDevCorrection()*getDevCorrection()*pow(exp(getThetaDevCorrection()*dev.dotprod())-1.,getPiDevCorrection()-1.)* exp(getThetaDevCorrection()*dev.dotprod())*getThetaDevCorrection();
+ Msg::Error("mlawHyperViscoElastic::evaluatePhiPCorrection: need to evaluate the volume energy");
}
else
{
A_v = getVolumeCorrection()*(tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection())+tanh(getZetaVolumeCorrection()));
dA_vdE = getVolumeCorrection()*getXiVolumeCorrection()*2./3.*(1.-tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection())*tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection()));
+ intA=0.;
+ if (getXiVolumeCorrection()!=0.)
+ {
+ double u=tr*tr;
+ intA=0.5*getVolumeCorrection()*(3./getXiVolumeCorrection()*(log(cosh(getXiVolumeCorrection()/3.*u-getZetaVolumeCorrection())) -log(cosh(-getZetaVolumeCorrection())) ) +tanh(getZetaVolumeCorrection())*u ) ;
+
+ }
B_d = getDevCorrection()*(tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection())+tanh(getPiDevCorrection()));
STensorOperation::zero(dB_vddev);
dB_vddev=dev;
dB_vddev*=2.*getPiDevCorrection()*getDevCorrection()*(1.-tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection())*tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection()) );
+ intB=0.;
+ if(getThetaDevCorrection()!=0)
+ {
+ double u=dev.dotprod();
+ intB=0.5*getDevCorrection()*(1./getThetaDevCorrection()*(log(cosh(getThetaDevCorrection()*u-getPiDevCorrection())) -log(cosh(-getPiDevCorrection())) ) +tanh(getPiDevCorrection())*u ) ;
+
+ }
+
+
+
+
}
}
@@ -140,10 +159,10 @@ double mlawHyperViscoElastic::deformationEnergy(const IPHyperViscoElastic& q) co
static STensor3 devEe;
STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
- double Av, dAv,B_d;
+ double Av, dAv,B_d, intA,intB;
static STensor3 dB_d;
- evaluatePhiPCorrection(trEe, devEe, Av, dAv, B_d, dB_d);
- Psy = getUpdatedBulkModulus(&q)*(1.+Av)*(0.5*trEe*trEe)+getUpdatedShearModulus(&q)*(1.+B_d)*STensorOperation::doubledot(devEe,devEe); //this is not correct we should do an integral
+ evaluatePhiPCorrection(trEe, devEe, Av, dAv, intA, B_d, dB_d, intB);
+ Psy = getUpdatedBulkModulus(&q)*(0.5*trEe*trEe+intA)+2.*getUpdatedShearModulus(&q)*(0.5*STensorOperation::doubledot(devEe,devEe)+intB); //this is not correct we should do an integral
//
//
for (int i=0; i<_N; i++){
@@ -417,7 +436,7 @@ void mlawHyperViscoElastic::updateViscoElasticFlow(const IPHyperViscoElastic *q0
};
void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1,
- double& Ke, double& Ge, double& A_v, double& B_d, double& dApr, STensor3& dB_vddev) const{
+ double& Ke, double& Ge, double& A_v, double& B_d, double& dApr, STensor3& dB_vddev, double &intA, double &intB) const{
if ((_Ki.size() > 0) or (_Gi.size() > 0)){
static STensor3 DE, devDE;
static double trDE;
@@ -455,7 +474,7 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
double p, trEe;
STensorOperation::decomposeDevTr(Ee,devE,trEe);
- evaluatePhiPCorrection(trEe, devE, A_v, dApr, B_d, dB_vddev);
+ evaluatePhiPCorrection(trEe, devE, A_v, dApr, intA, B_d, dB_vddev, intB);
devK = devE;
devK *= (2.*getUpdatedShearModulus(q1)*(1.+B_d)); // deviatoric part
@@ -547,7 +566,7 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
STensorOperation::decomposeDevTr(Ee,devE,trEe);
STensorOperation::zero(dB_vddev);
- evaluatePhiPCorrection(trEe, devE, A_v, dApr, B_d, dB_vddev);
+ evaluatePhiPCorrection(trEe, devE, A_v, dApr, intA, B_d, dB_vddev, intB);
devK = devE;
devK *= (2.*getUpdatedShearModulus(q1)*(1.+B_d)); // deviatoric part
@@ -602,10 +621,10 @@ void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0,
STensorOperation::decomposeDevTr(E,devE,trE);
double Ke, Ge;
- double A_v_pr=0., B_d_pr = 0., dApr=0.;
+ double A_v_pr=0., B_d_pr = 0., dApr=0., intApr=0., intBpr=0.;
static STensor3 dB_devpr;
STensorOperation::zero(dB_devpr);
- viscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
+ viscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr,intApr,intBpr);
const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
static STensor3 secondPK;
@@ -894,10 +913,10 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
- double A_v_pr, B_d_pr, dApr;
+ double A_v_pr=0., B_d_pr=0., dApr=0.,intApr=0.,intBpr=0.;
static STensor3 dB_devpr;
STensorOperation::zero(dB_devpr);
- viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
+ viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr, intApr, intBpr);
STensorOperation::decomposeDevTr(Ee,devEepr,trEepr);
//std::cout << "trEepr = "<<trEepr <<std::endl;
@@ -980,7 +999,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
double A_v=0.;
double B_d = 0.;
- double dAdtrEe=0.;
+ double dAdtrEe=0., intA=0., intB=0.;
static STensor3 dBddevEe;
STensorOperation::zero(dBddevEe);
@@ -1078,7 +1097,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
while (1){
- evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, B_d, dBddevEe);
+ evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, intA, B_d, dBddevEe, intB);
J=getUpdatedBulkModulus(q1)*A_v*trEe-ptilde_Vol;
dJdptilde_Vol = getUpdatedBulkModulus(q1)*( dAdtrEe*dtrEedptilde*trEe+A_v*dtrEedptilde)-1.;
@@ -1121,7 +1140,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
while (1){
- evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, B_d, dBddevEe);
+ evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, intA, B_d, dBddevEe, intB);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
@@ -1735,10 +1754,10 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c
- double A_v_pr=0., B_d_pr = 0., dApr=0.;
+ double A_v_pr=0., B_d_pr = 0., dApr=0., intApr=0., intBpr=0.;
static STensor3 dB_devpr;
STensorOperation::zero(dB_devpr);
- viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
+ viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr, intApr, intBpr);
static STensor3 devKpr; // dev corotational kirchoff stress predictor
static double ppr; // pressure predictor
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.h b/NonLinearSolver/materialLaw/mlawHyperelastic.h
index 4118280b8..08d5d5531 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.h
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.h
@@ -59,14 +59,14 @@ class mlawHyperViscoElastic : public materialLaw{
void updateViscoElasticFlow(const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1, double& Ke, double& Ge) const;
void viscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1,
- double& Ke, double& Ge, double& A_v, double& B_d, double& dAprd, STensor3& dB_vddev) const;
+ double& Ke, double& Ge, double& A_v, double& B_d, double& dAprd, STensor3& dB_vddev, double &intA, double &intB) const;
void isotropicHookTensor(const double K, const double G, STensor43& L) const;
void predictorCorrector_ViscoElastic(const STensor3& F0, const STensor3& F, STensor3&P, const IPHyperViscoElastic *q0_, IPHyperViscoElastic *q1,
const bool stiff, STensor43& Tangent) const;
- void evaluatePhiPCorrection(double trEe, const STensor3 &devEpr, double &A_v, double &dAprdEepr, double &B_d, STensor3 &dB_vddevEe) const;
+ void evaluatePhiPCorrection(double trEe, const STensor3 &devEpr, double &A_v, double &dAprdEepr, double &intA, double &B_d, STensor3 &dB_vddevEe, double &intB) const;
//void evaluateA_dA(double trEe, double &A_v, double &dAprdEepr) const;
#endif // SWIG
--
GitLab