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.

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

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

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

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

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

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

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

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

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) {};

CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::mField

MoFEM::Interface & mField

Definition CohesiveInterfaceElement.hpp:53

CohesiveElement::CohesiveInterfaceElement::PhysicalEquation::isInitialised

bool isInitialised

Definition CohesiveInterfaceElement.hpp:54

◆ ~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);
    }