Constitutive (physical) equation for interface. More...

#include <users_modules/basic_finite_elements/cohesive_interface/src/CohesiveInterfaceElement.hpp>

Collaboration diagram for CohesiveElement::CohesiveInterfaceElement::PhysicalEquation:

Public Member Functions
	PhysicalEquation (MoFEM::Interface &m_field)

virtual	~PhysicalEquation ()

MoFEMErrorCode	iNitailise (const FEMethod *fe_method)
	Initialize history variable data. More...

double	calcG (int gg, MatrixDouble gap_loc)
	Calculate gap opening. More...

MoFEMErrorCode	getKappa (int nb_gauss_pts, const FEMethod *fe_method)
	Get pointer from the mesh to histoy variables \(\kappa\). More...

MoFEMErrorCode	calcDglob (const double omega, MatrixDouble &R)
	Calculate stiffness material matrix. More...

MoFEMErrorCode	calcOmega (const double kappa, double &omega)
	Calculate damage. More...

MoFEMErrorCode	calcTangetDglob (const double omega, double g, const VectorDouble &gap_loc, MatrixDouble &R)
	Calculate tangent material stiffness. More...

virtual MoFEMErrorCode	calculateTraction (VectorDouble &traction, int gg, CommonData &common_data, const FEMethod *fe_method)
	Calculate tractions. More...

virtual MoFEMErrorCode	calculateTangentStiffeness (MatrixDouble &tangent_matrix, int gg, CommonData &common_data, const FEMethod *fe_method)
	Calculate tangent stiffness. More...

virtual MoFEMErrorCode	updateHistory (CommonData &common_data, const FEMethod *fe_method)
	Update history variables when converged. More...

Public Attributes
MoFEM::Interface &	mField

bool	isInitialised

double	h

double	youngModulus

double	beta

double	ft

double	Gf

Range	pRisms

Tag	thKappa

Tag	thDamagedPrism

double	E0

double	g0

double	kappa1

double *	kappaPtr

int	kappaSize

MatrixDouble	Dglob

MatrixDouble	Dloc

Detailed Description

Constitutive (physical) equation for interface.

This is linear degradation model. Material parameters are: strength \(f_t\), interface fracture energy \(G_f\), elastic material stiffness \(E\). Parameter \(\beta\) controls how interface opening is calculated.

Model parameter is interface penalty thickness \(h\).

Definition at line 51 of file CohesiveInterfaceElement.hpp.

Constructor & Destructor Documentation

◆ PhysicalEquation()

CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::PhysicalEquation ( MoFEM::Interface & m_field )

inline

Definition at line 55 of file CohesiveInterfaceElement.hpp.

55 :

56 mField(m_field),isInitialised(false) {};

◆ ~PhysicalEquation()

virtual CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::~PhysicalEquation ( )

inlinevirtual

Definition at line 58 of file CohesiveInterfaceElement.hpp.

58 {

59 }

Member Function Documentation

◆ calcDglob()

MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcDglob	(	const double	omega,
		MatrixDouble &	R
	)

inline

Calculate stiffness material matrix.

\[ \mathbf{D}_\textrm{loc} = (1-\Omega) \mathbf{I} E_0 \]

where \(E_0\) is initial interface penalty stiffness

\[ \mathbf{D}_\textrm{glob} = \mathbf{R}^\textrm{T} \mathbf{D}_\textrm{loc}\mathbf{R} \]

Definition at line 136 of file CohesiveInterfaceElement.hpp.

                                                                  {
       MoFEMFunctionBeginHot;
       Dglob.resize(3,3);
       Dloc.resize(3,3);
       Dloc.clear();
       double E = (1-omega)*E0;
       Dloc(0,0) = E;
       Dloc(1,1) = E;
       Dloc(2,2) = E;
       Dglob = prod( Dloc, R );
       Dglob = prod( trans(R), Dglob );
       MoFEMFunctionReturnHot(0);
     }

◆ calcG()

double CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcG	(	int	gg,
		MatrixDouble	gap_loc
	)

inline

Calculate gap opening.

\[ g = \sqrt{ g_n^2 + \beta(g_{s1}^2 + g_{s2}^2)} \]

Definition at line 95 of file CohesiveInterfaceElement.hpp.

                                               {
       return sqrt(pow(gap_loc(gg,0),2)+beta*(pow(gap_loc(gg,1),2)+pow(gap_loc(gg,2),2)));
     }

◆ calcOmega()

MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcOmega	(	const double	kappa,
		double &	omega
	)

inline

Calculate damage.

\[ \Omega = \frac{1}{2} \frac{(2 G_f E_0+f_t^2)\kappa}{(ft+E_0 \kappa)G_f} \]

Definition at line 157 of file CohesiveInterfaceElement.hpp.

                                                                {
       MoFEMFunctionBeginHot;
       omega = 0;
       if(kappa>=kappa1) {
         omega = 1;
         MoFEMFunctionReturnHot(0);
       } else if(kappa>0) {
         double a = (2.0*Gf*E0+ft*ft)*kappa;
         double b = (ft+E0*kappa)*Gf;
         omega = 0.5*a/b;
       }
       MoFEMFunctionReturnHot(0);
     }

◆ calcTangetDglob()

MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calcTangetDglob	(	const double	omega,
		double	g,
		const VectorDouble &	gap_loc,
		MatrixDouble &	R
	)

inline

Calculate tangent material stiffness.

Definition at line 173 of file CohesiveInterfaceElement.hpp.

                                                                                                             {
       MoFEMFunctionBeginHot;
       Dglob.resize(3,3);
       Dloc.resize(3,3);
       double domega = 0.5*(2*Gf*E0+ft*ft)/((ft+(g-ft/E0)*E0)*Gf) - 0.5*((g-ft/E0)*(2*Gf*E0+ft*ft)*E0)/(pow(ft+(g-ft/E0)*E0,2)*Gf);
       Dloc.resize(3,3);
       //r0
       Dloc(0,0) = (1-omega)*E0 - domega*E0*gap_loc[0]*gap_loc[0]/g;
       Dloc(0,1) = -domega*E0*gap_loc[0]*beta*gap_loc[1]/g;
       Dloc(0,2) = -domega*E0*gap_loc[0]*beta*gap_loc[2]/g;
       //r1
       Dloc(1,0) = -domega*E0*gap_loc[1]*gap_loc[0]/g;
       Dloc(1,1) = (1-omega)*E0 - domega*E0*gap_loc[1]*beta*gap_loc[1]/g;
       Dloc(1,2) = -domega*E0*gap_loc[1]*beta*gap_loc[2]/g;
       //r2
       Dloc(2,0) = -domega*E0*gap_loc[2]*gap_loc[0]/g;
       Dloc(2,1) = -domega*E0*gap_loc[2]*beta*gap_loc[1]/g;
       Dloc(2,2) = (1-omega)*E0 - domega*E0*gap_loc[2]*beta*gap_loc[2]/g;
       Dglob = prod(Dloc,R);
       Dglob = prod(trans(R),Dglob);
       MoFEMFunctionReturnHot(0);
     }

◆ calculateTangentStiffeness()

virtual MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calculateTangentStiffeness	(	MatrixDouble &	tangent_matrix,
		int	gg,
		CommonData &	common_data,
		const FEMethod *	fe_method
	)

inlinevirtual

Calculate tangent stiffness.

Definition at line 231 of file CohesiveInterfaceElement.hpp.

       {
       MoFEMFunctionBeginHot;
  
       try {
         if(!isInitialised) {
           ierr = iNitailise(fe_method); CHKERRG(ierr);
           isInitialised = true;
         }
         if(gg==0) {
           ierr = getKappa(common_data.gapGlob.size1(),fe_method); CHKERRG(ierr);
         }
         double g = calcG(gg,common_data.gapLoc);
         double kappa = fmax(g-g0,kappaPtr[gg]);
         double omega = 0;
         ierr = calcOmega(kappa,omega); CHKERRG(ierr);
         //std::cerr << gg << " " << omega << std::endl;
         int iter;
         ierr = SNESGetIterationNumber(fe_method->snes,&iter); CHKERRG(ierr);
         if((kappa <= kappaPtr[gg])||(kappa>=kappa1)||(iter <= 1)) {
           ierr = calcDglob(omega,common_data.R[gg]); CHKERRG(ierr);
         } else {
           ublas::matrix_row<MatrixDouble > g_loc(common_data.gapLoc,gg);
           ierr = calcTangetDglob(omega,g,g_loc,common_data.R[gg]);
         }
         tangent_matrix.resize(3,3);
         noalias(tangent_matrix) = Dglob;
         //std::cerr << "t " << tangent_matrix << std::endl;
       } catch (const std::exception& ex) {
         std::ostringstream ss;
         ss << "throw in method: " << ex.what() << std::endl;
         SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
       }
       MoFEMFunctionReturnHot(0);
     }

◆ calculateTraction()

virtual MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::calculateTraction	(	VectorDouble &	traction,
		int	gg,
		CommonData &	common_data,
		const FEMethod *	fe_method
	)

inlinevirtual

Calculate tractions.

\[ \mathbf{t} = \mathbf{D}_\textrm{glob}\mathbf{g} \]

Definition at line 203 of file CohesiveInterfaceElement.hpp.

       {
       MoFEMFunctionBeginHot;
  
       if(!isInitialised) {
         ierr = iNitailise(fe_method); CHKERRG(ierr);
         isInitialised = true;
       }
       if(gg==0) {
         ierr = getKappa(common_data.gapGlob.size1(),fe_method); CHKERRG(ierr);
       }
       double g = calcG(gg,common_data.gapLoc);
       double kappa = fmax(g-g0,kappaPtr[gg]);
       double omega = 0;
       ierr = calcOmega(kappa,omega); CHKERRG(ierr);
       //std::cerr << gg << " " << omega << std::endl;
       ierr = calcDglob(omega,common_data.R[gg]); CHKERRG(ierr);
       traction.resize(3);
       ublas::matrix_row<MatrixDouble > gap_glob(common_data.gapGlob,gg);
       noalias(traction) = prod(Dglob,gap_glob);
       MoFEMFunctionReturnHot(0);
     }

◆ getKappa()

MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::getKappa	(	int	nb_gauss_pts,
		const FEMethod *	fe_method
	)

inline

Get pointer from the mesh to histoy variables \(\kappa\).

Definition at line 104 of file CohesiveInterfaceElement.hpp.

                                                                         {
       MoFEMFunctionBeginHot;
       EntityHandle ent = fe_method->numeredEntFiniteElementPtr->getEnt();
  
       rval = mField.get_moab().tag_get_by_ptr(thKappa,&ent,1,(const void **)&kappaPtr,&kappaSize);
       if(rval != MB_SUCCESS || kappaSize != nb_gauss_pts) {
         VectorDouble kappa;
         kappa.resize(nb_gauss_pts);
         kappa.clear();
         int tag_size[1];
         tag_size[0] = nb_gauss_pts;
         void const* tag_data[] = { &kappa[0] };
         rval = mField.get_moab().tag_set_by_ptr(thKappa,&ent,1,tag_data,tag_size);  CHKERRG(rval);
         rval = mField.get_moab().tag_get_by_ptr(thKappa,&ent,1,(const void **)&kappaPtr,&kappaSize);  CHKERRG(rval);
       }
       MoFEMFunctionReturnHot(0);
     }

◆ iNitailise()

MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::iNitailise ( const FEMethod * fe_method )

inline

Initialize history variable data.

Create tag on the prism/interface to store damage history variable

Definition at line 72 of file CohesiveInterfaceElement.hpp.

                                                          {
       MoFEMFunctionBeginHot;
  
       double def_damaged = 0;
       rval = mField.get_moab().tag_get_handle(
         "DAMAGED_PRISM",1,MB_TYPE_INTEGER,thDamagedPrism,MB_TAG_CREAT|MB_TAG_SPARSE,&def_damaged
       ); MOAB_THROW(rval);
       const int def_len = 0;
       rval = mField.get_moab().tag_get_handle("_KAPPA",def_len,MB_TYPE_DOUBLE,
       thKappa,MB_TAG_CREAT|MB_TAG_SPARSE|MB_TAG_VARLEN,NULL); CHKERRG(rval);
       E0 = youngModulus/h;
       g0 = ft/E0;
       kappa1 = 2*Gf/ft;
       MoFEMFunctionReturnHot(0);
     }

◆ updateHistory()

virtual MoFEMErrorCode CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::updateHistory	(	CommonData &	common_data,
		const FEMethod *	fe_method
	)

inlinevirtual

Update history variables when converged.

Definition at line 272 of file CohesiveInterfaceElement.hpp.

       {
       MoFEMFunctionBeginHot;
  
  
       if(!isInitialised) {
         ierr = iNitailise(fe_method); CHKERRG(ierr);
         isInitialised = true;
       }
       ierr = getKappa(common_data.gapGlob.size1(),fe_method); CHKERRG(ierr);
       bool all_gauss_pts_damaged = true;
       for(unsigned int gg = 0;gg<common_data.gapGlob.size1();gg++) {
         double omega = 0;
         double g = calcG(gg,common_data.gapLoc);
         double kappa = fmax(g-g0,kappaPtr[gg]);
         kappaPtr[gg] = kappa;
         ierr = calcOmega(kappa,omega); CHKERRG(ierr);
         //if(omega < 1.) {
         all_gauss_pts_damaged = false;
         //}
       }
       if(all_gauss_pts_damaged) {
         EntityHandle ent = fe_method->numeredEntFiniteElementPtr->getEnt();
         int set_prism_as_demaged = 1;
         rval = mField.get_moab().tag_set_data(thDamagedPrism,&ent,1,&set_prism_as_demaged); CHKERRG(rval);
       }
       MoFEMFunctionReturnHot(0);
     }