BC example 1

Description

This boundary condition is used to impose on a the displacement field given by the linear solution at a crack tip. The crack is loaded under mode I or II. This boundary condition is limited to 2D plane problems. It will work on axisymmetric problems , however its physical meaning will be uncertain.

The displacement field is given by:

Mode I:

\[\begin{split}\begin{aligned} u_1 = {K_I\over2\mu}\sqrt{{r\over2\pi}}\cos{\theta\over2} \left(\kappa-1+2\sin^2{\theta\over2}\right) \\ u_2 = {K_I\over2\mu}\sqrt{{r\over2\pi}}\sin{\theta\over2} \left(\kappa+1-2\cos^2{\theta\over2}\right) \\ \end{aligned}\end{split}\]

Mode II:

\[\begin{split}\begin{aligned} u_1 = {K_{II}\over2\mu}\sqrt{{r\over2\pi}}\sin{\theta\over2} \left(\kappa+1+2\cos^2{\theta\over2}\right) \\ u_2 = -{K_{II}\over2\mu}\sqrt{{r\over2\pi}}\cos{\theta\over2} \left(\kappa-1-2\sin^2{\theta\over2}\right) \\ \end{aligned}\end{split}\]

with

\[\begin{split}\kappa = \begin{cases} 3-4\nu & plane strain \\ (3-\nu)/(1+\nu) & plane stress \\ \end{cases}\end{split}\]

For elastic problems:

\[\begin{split}J = \begin{cases} (1-\nu^2)\left(K_I^2+K^2_{II}\right)/E \hskip2cm plane strain \\ \left(K_I^2+K^2_{II}\right)/E \hskip2cm plane stress \\ \end{cases}\end{split}\]

Code Listing

class BC_K_FIELD : public BC {
     NNSET* set;
     VECTOR center;
     double poisson, young, kappa, two_mu;
     int mode;
  public :
     BC_K_FIELD();
     virtual void initialize(ASCII_FILE&,MESH&,int);
     virtual ~BC_K_FIELD() {}
     virtual void _update(MESH&);
     RTTI_INFO;
};

IMPL_RTTI_INFO(BC,BC_K_FIELD);

DECLARE_OBJECT(BC,BC_K_FIELD,K_field)

BC_K_FIELD::BC_K_FIELD() {}

void BC_K_FIELD::initialize(ASCII_FILE& file,MESH& mesh,int)
{ set = (NNSET*)mesh.find_nset(file.getSTRING());
  center  = file.getVECTOR(2); if(!file.ok) VEC_REQ("center",file);
  young   = file.getdouble();  if(!file.ok) DBL_REQ("young",file);
  poisson = file.getdouble();  if(!file.ok) DBL_REQ("poisson",file);
  STRING str=file.getSTRING();
  if     (str=="plane_stress") kappa = (3.-poisson)/(1.+poisson);
  else if(str=="plane_strain") kappa = 3.-4.*poisson;
  else INPUT_ERROR("expecting plane_stress or plane_strain");
  two_mu = young/(1.+poisson);
  str=file.getSTRING();
  if     (str=="I") mode = 1;
  else if(str=="II") mode=2;
  else INPUT_ERROR("expecting I or II to specify loading mode");
  base_read(file);
}

void BC_K_FIELD::_update(MESH&)
{ static DOF_TYPE Displ_name[2]={DOF::translate("U1",DOF_NODE_ID),
                                 DOF::translate("U2",DOF_NODE_ID)};
  double Komu=compute_increment()[0]/two_mu, r, thetas2, du1, du2, fact;
  VECTOR axx(2), vr; axx[0]=1.; axx[1]=0.;
  for(int i=0;i<!(*set);i++) {
     vr = (*set)[i]->coord-center;
     r = norm(vr);
     thetas2 = get_angle(axx,vr)/2.;
     fact = Komu*sqrt(r/2./M_PI);
     if(mode==1) {
         du1 = fact*cos(thetas2)*
               (kappa-1.+2*sin(thetas2)*sin(thetas2));
         du2 = fact*sin(thetas2)*
               (kappa+1.-2.*cos(thetas2)*cos(thetas2));
     } else {
         du1 = fact*sin(thetas2)*
               (kappa+1.+2.*cos(thetas2)*cos(thetas2));
         du2 = -fact*cos(thetas2)*
               (kappa-1.-2*sin(thetas2)*sin(thetas2));
     }
     (*set)[i]->add_fixed_dof(Displ_name[0],du1);
     (*set)[i]->add_fixed_dof(Displ_name[1],du2);
  }
}