From 45d1ff31dd24faf920fc4c4c5239285083daa0c5 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Fri, 4 Nov 2016 06:35:18 +0100 Subject: [PATCH] squash! Added: Tangent matrix for monolithic coupling of elasticity and phasefield. Added: Option to return derivatives of the stress degradation function. Added: Parsing and printing of the stabilization parameter alpha, this is the same quantity as stabk in the CahnHilliard integrand. The monolithic coupling matrix was incorrect --- FractureElasticityVoigt.C | 58 ++++++++++++++++++----------------- FractureElasticityVoigt.h | 2 +- Test/TestFractureElasticity.C | 2 +- 3 files changed, 32 insertions(+), 30 deletions(-) diff --git a/FractureElasticityVoigt.C b/FractureElasticityVoigt.C index 22b0066..983fb3e 100644 --- a/FractureElasticityVoigt.C +++ b/FractureElasticityVoigt.C @@ -28,8 +28,9 @@ bool FractureElasticityVoigt::evalStress (double lambda, double mu, double Gc, const SymmTensor& epsil, double* Phi, - SymmTensor* sigma, Matrix* dSdE, - bool postProc, bool printElm) const + SymmTensor* sigma, SymmTensor* sigP, + Matrix* dSdE, bool postProc, + bool printElm) const { PROFILE3("FractureEl::evalStress"); @@ -90,6 +91,12 @@ bool FractureElasticityVoigt::evalStress (double lambda, double mu, double Gc, *sigma = C0*trEps; *sigma += 2.0*mu*(Gc*ePos + eNeg); } + if (sigP) + { + // Evaluate the non-degraded tensile stress tensor + *sigP = lambda*trEps; + *sigP += 2.0*mu*ePos; + } // Evaluate the tensile energy Phi[0] = mu*(ePos*ePos).trace(); @@ -113,6 +120,7 @@ bool FractureElasticityVoigt::evalStress (double lambda, double mu, double Gc, std::cout <<"M("<< 1+a <<") =\n"<< M[a]; std::cout <<"ePos =\n"<< ePos <<"eNeg =\n"<< eNeg; if (sigma) std::cout <<"sigma =\n"<< *sigma; + if (sigP) std::cout <<"sigP =\n"<< *sigP; std::cout <<"Phi = "<< Phi[0]; if (postProc) std::cout <<" "<< Phi[1] <<" "<< Phi[2] <<" "<< Phi[3]; std::cout << std::endl; @@ -242,7 +250,7 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt, size_t nstrc = (nsd+1)*nsd/2; Matrix Bmat, dSdE(nstrc,nstrc); - SymmTensor eps(nsd), sigma(nsd); + SymmTensor eps(nsd), sigma(nsd), sigP(nsd); bool lHaveStrains = false; if (eKm || eKg || iS || m_mode == SIM::RECOVERY) @@ -269,12 +277,12 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt, if (noSplit) { // Evaluate the constitutive matrix and the stress tensor at this point - if (!material->evaluate(dSdE,sigma,myPhi[fe.iGP],fe,X,eps,eps)) + if (!material->evaluate(dSdE,sigP,myPhi[fe.iGP],fe,X,eps,eps)) return false; // Degrade the stresses and strain energy isotropically myPhi[fe.iGP] *= Gc; - sigma *= Gc; + sigma = Gc*sigP; } else { @@ -290,7 +298,7 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt, // Evaluate the stress state at this point, with degraded tensile part if (!this->evalStress(lambda,mu,Gc,eps,&myPhi[fe.iGP],&sigma, - eKm ? &dSdE : nullptr)) + eKc ? &sigP : nullptr, eKm ? &dSdE : nullptr)) return false; } } @@ -312,6 +320,14 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt, if (eM) // Integrate the mass matrix this->formMassMatrix(elMat.A[eM-1],fe.N,X,fe.detJxW); + if (iS && lHaveStrains) + { + // Integrate the internal forces + sigma *= -fe.detJxW; + if (!Bmat.multiply(sigma,elMat.b[iS-1],true,true)) // ES -= B^T*sigma + return false; + } + if (eS) { // Integrate the load vector due to gravitation and other body forces @@ -320,28 +336,14 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt, this->formCrackForce(elMat.b[eS-1],elMat.vec,fe,X); } - if (!lHaveStrains) - return true; - - Vector BtSig; - sigma *= -fe.detJxW; - if (eKc && !Bmat.multiply(sigma,BtSig,true)) // BtSig = -B^T*sigma - return false; - - if (iS) - { - // Integrate the internal forces - if (eKc) - elMat.b[iS-1].add(BtSig); // ES += BtSig - else if (!Bmat.multiply(sigma,elMat.b[iS-1],true,true)) // ES -= B^T*sigma - return false; - } - - if (eKc) + if (eKc && lHaveStrains) { // Integrate the coupling matrix between displacement and phase-field DOFs - BtSig *= -this->getStressDegradation(fe.N,elmInt.vec,1); // g'(c) - if (!elMat.A[eKc-1].outer_product(BtSig,fe.N,true)) // Kc += B^T*sigma*N*dGc + Vector BtSig; + sigP *= this->getStressDegradation(fe.N,elmInt.vec,1)*fe.detJxW; // *= g'(c) + if (!Bmat.multiply(sigP,BtSig,true)) // BtSig = B^T*g'(c)*sigP + return false; + if (!elMat.A[eKc-1].outer_product(BtSig,fe.N,true)) // Kc += BtSig*N return false; } @@ -353,7 +355,7 @@ bool FractureElasticityVoigt::evalStress (double lambda, double mu, double Gc, const SymmTensor& epsilon, double* Phi, SymmTensor& sigma) const { - return this->evalStress(lambda,mu,Gc,epsilon,Phi,&sigma,nullptr,true); + return this->evalStress(lambda,mu,Gc,epsilon,Phi,&sigma,nullptr,nullptr,true); } @@ -445,7 +447,7 @@ bool FractureElasticNorm::evalInt (LocalIntegral& elmInt, // Evaluate the tensile-degraded strain energy if (!p.evalStress(lambda,mu,Gc,eps,Phi, dirDer ? &sigma : nullptr, - nullptr,true,printElm)) + nullptr,nullptr,true,printElm)) return false; } diff --git a/FractureElasticityVoigt.h b/FractureElasticityVoigt.h index fd9460f..aa4c239 100644 --- a/FractureElasticityVoigt.h +++ b/FractureElasticityVoigt.h @@ -62,7 +62,7 @@ class FractureElasticityVoigt : public FractureElasticity //! \brief Evaluates the stress tensor and its derivative w.r.t. the strains. bool evalStress(double lambda, double mu, double Gc, const SymmTensor& epsilon, double* Phi, - SymmTensor* sigma, Matrix* dSdE, + SymmTensor* sigma, SymmTensor* sigP, Matrix* dSdE, bool postProc = false, bool printElm = false) const; unsigned short int eKc; //!< Zero-based index to element coupling matrix diff --git a/Test/TestFractureElasticity.C b/Test/TestFractureElasticity.C index a55c65a..8358e9c 100644 --- a/Test/TestFractureElasticity.C +++ b/Test/TestFractureElasticity.C @@ -29,7 +29,7 @@ public: virtual ~FracEl() {} bool calcStress(double lambda, double mu, double Gc, const SymmTensor& eps, double& Phi, SymmTensor& sigma, Matrix& dSdE) const - { return evalStress(lambda,mu,Gc,eps,&Phi,&sigma,&dSdE); } + { return evalStress(lambda,mu,Gc,eps,&Phi,&sigma,nullptr,&dSdE); } bool calcStress(double lambda, double mu, double Gc, const SymmTensor& eps, double& Phi, SymmTensor& sigma, Tensor4& dSdE) const { return FractureElasticity::evalStress(lambda,mu,Gc,eps,&Phi,sigma,&dSdE); }