elmt_shell_4n.c

有限元程序

         }
         A = sqrt(A_2);
         R = p->LC_ptr->R[kk-1];
         eff_pl_strain = p->effect_pl_strain[kk-1];

       /* Step 2: comparision                        */

        if(A <= R) { /* ELASTIC DEFORMATION */
          for(i = 1; i <= dof; i++) 
             p->stress->uMatrix.daa[i-1][kk-1] = stress[i-1][0];

           SaveRespondBuffer(p, kk);

#ifdef DEBUG1
    printf(" +++++++ ELASTIC DEFORMATION A = %lf R = %lf \n", A, R);
    printf(" at elmt_no = %d, integ_pt = %d\n", p->elmt_no, kk);
#endif
        }
        else {       /* PLASTIC DEFORMATION */
        

       /* Step 3 Estimate number of sub-incrementations needed     */

          if( ABS(p->LC_ptr->beta) < 1E-10){ /* Only for beta = 0 kinematic hardening case */

       /* Step 3.1 Estimate the effective plastic strain increment */

              temp = sqrt(3.0/2.0)*(A-R)/(p->LC_ptr->H[kk-1]+3.0*G);

       /* Step 3.2 Estimate the back stress increment              */
           
              /* Estimate H' */
              if(!strcmp(p->material_name, "ELASTIC_PERFECTLY_PLASTIC")) {
                  H = 0.0;
              }else {
                  if(!strcmp(p->LC_ptr->name, "Ramberg-Osgood")) {
                      effect_stress = A*sqrt(3.0/2.0);
                      Load_Curve(p, &H, effect_stress, E,fy);
                  }
                  if(!strcmp(p->LC_ptr->name, "Bi-Linear")) 
                      H = p->LC_ptr->H[kk-1];
             }
          /* calculate the effective plastic strain incremental */

             temp1 = sqrt(2.0/3.0)*(H + 3.0*G);
             temp  = A/temp1 - R/temp1;            /* eff_pl_strain_incr */

       /* Step 3.3: Compute Lambda  and pl_strain_incr   */
       /*         Lambda = sqrt(3/2)*eff_pl_strain_incr  */
       /* plastic strain incr is now stored in strain_incr */

             Lambda = sqrt(3.0/2.0)*temp;
             for(i = 1; i <= dof; i++) 
                 strain_incr[i-1][0] = Lambda*stress_dev[i-1][0]/A;

          /* Step 3.4:  calculate back stress increment  */
          /* beta = 0 for kinematic hardening */
          /* beta = 1 for isotropic hardening */

              for(i = 1; i <= dof; i++) {
                  back_stress_incr[i-1][0] = H*strain_incr[i-1][0]*2.0/3.0;
              }
          }

          mean_stress = (stress_incr[0][0] - back_stress_incr[0][0]
                        +stress_incr[1][0] - back_stress_incr[1][0])/3.0;
          temp   = 0.0;
          for(i = 1; i <= dof; i++) {
              if(i <= 2)
                 stress_dev[i-1][0] = stress_incr[i-1][0] - mean_stress - back_stress_incr[i-1][0];
              else
                 stress_dev[i-1][0] = stress_incr[i-1][0] - back_stress_incr[i-1][0];

              temp += stress_dev[i-1][0]*stress_dev[i-1][0];
          }
          temp = sqrt(temp);
          iNo_iter_step = (int) (2*temp/R/Beta1) + 1 ;

#ifdef DEBUG1
    printf(" ******Plastic DEFORMATION A  = %lf, R = %lf \n", A, R);
    printf(" ******Plastic DEFORMATION dA = %lf, R = %lf \n", temp, R);
    printf(" at elmt_no = %d, integ_pt = %d\n", p->elmt_no, kk);
    printf(" No of sub-Incremental steps = %d \n",iNo_iter_step);
#endif
#undef DEBUG1
          /* Step 4 Start sub-incrementation  */

         /* copy the plastic strain before sub-incrementation */

         eff_pl_strainTemp = eff_pl_strain;
         for(i = 1; i <= dof; i++) 
             strain_pl[i-1][0] = p->strain_pl->uMatrix.daa[i-1][kk-1];

         switch(iNo_iter_step) {
             case 1: 
               ii = 1;
               Plastic_Deform(p, &H, &R, &eff_pl_strain, stress, stress_dev,
                              strain_incr, E, fy, G, A, iNo_iter_step, dof, ii, kk);
               SaveRespondBuffer(p, kk);
#ifdef DEBUG
   printf(" A= %lf R = %lf H = %lf eff_pl_strain = %lf\n", A, R, H, eff_pl_strain);
#endif
             break;
             default: 
               for(i = 1; i <= p->dof_per_node; i++) {
                   for(j = 1; j <= p->nodes_per_elmt; j++) {
                       k = p->dof_per_node*(j-1)+i;
                       displ_incr[k-1][0] 
                       = p->displ_incr->uMatrix.daa[i-1][j-1]/((double) iNo_iter_step);
                   }
               }
               /* Calculate stress increment */

               strain_incr = dMatrixMultRep(strain_incr,B_matrix, dof, 
                                            size, displ_incr, size, 1);
               stress_incr = dMatrixMultRep(stress_incr, m1, dof, dof,
                                            strain_incr, dof, 1);
               for(ii = 1; ii <= iNo_iter_step; ii++) {
                   /* Trial stress */
                   if(ii == 1) {
                      for(i = 1; i <= dof; i++) 
                          stress[i-1][0] = stress_incr[i-1][0] + 
                                           p->stress->uMatrix.daa[i-1][kk-1];
                   }
                   else
                      for(i = 1; i <= dof; i++) 
                          stress[i-1][0] += stress_incr[i-1][0];

                   mean_stress = (stress[0][0] - p->LC_ptr->back_stress[0][kk-1]
                                 +stress[1][0] - p->LC_ptr->back_stress[1][kk-1])/3.0;
                   A_2 = 0.0;
                   for(i = 1; i <= dof; i++) {
                       if(i <= 2)
                          stress_dev[i-1][0] = stress[i-1][0] - mean_stress
                                             - p->LC_ptr->back_stress[i-1][kk-1];
                       else
                          stress_dev[i-1][0] = stress[i-1][0]
                                             - p->LC_ptr->back_stress[i-1][kk-1];
                       A_2 += stress_dev[i-1][0]*stress_dev[i-1][0];
                   }
                   A = sqrt(A_2);
                   if(A <= R) { /* ELASTIC DEFORMATION */
                      if(i == iNo_iter_step){
                         for(i = 1; i <= dof; i++)
                             p->stress->uMatrix.daa[i-1][kk-1] = stress[i-1][0];
                      }
                   /* go to next sub incremental iteration */

                   }else {   /* PLASTIC DEFORMATION */

                      Plastic_Deform(p, &H, &R, &eff_pl_strain, stress,stress_dev,
                                     strain_incr, E, fy, G, A, iNo_iter_step, dof, ii, kk);
                   }
               } /* end of sub-incremental iteration */
             break;
          } /* end of switch for sub incrementation */
          
          p->LC_ptr->R[kk-1] = R;
          p->LC_ptr->H[kk-1] = H;
          p->effect_pl_strain[kk-1]   = eff_pl_strain;
          p->eff_pl_strain_incr[kk-1] = eff_pl_strain - eff_pl_strainTemp;
          for(i = 1; i <= dof; i++) {
              p->strain_pl_incr->uMatrix.daa[i-1][kk-1]
              = p->strain_pl->uMatrix.daa[i-1][kk-1] - strain_pl[i-1][0];
          }
          SaveRespondBuffer(p, kk);
        }
       break;
       default:
         printf(" In Stress_Update(): elmt_no \n", p->elmt_no);
         printf(" elmt_state = 0 : Elastic_deformation \n");
         printf(" elmt_state = 1 : plastic_deformation \n");
         printf(" elmt_state = %d: p->elmt_state \n");
         FatalError(" Unknown elmt state ",(char *)NULL);
       break;
     }
  }

  /* ASSIGN UNITS TO p ARRAY */
  
   if(CheckUnits() == ON) {
         switch(CheckUnitsType()) {
           case SI:
             dimen = DefaultUnits("Pa");
           break;
           case US:
             dimen = DefaultUnits("psi");
           break;
       }
       for(i = 1; i <= dof; i++)
           p->stress->spRowUnits[i-1] = *DefaultUnits("psi");

       free((char *) dimen->units_name);
       free((char *) dimen);
   }

   MatrixFreeIndirectDouble(m1, dof);
   MatrixFreeIndirectDouble(strain_incr, dof);
   MatrixFreeIndirectDouble(stress, dof); 
   MatrixFreeIndirectDouble(stress_dev,dof);
   MatrixFreeIndirectDouble(stress_incr,dof);
   MatrixFreeIndirectDouble(back_stress_incr,dof);

#ifdef DEBUG
      dMatrixPrint("p->stress in Stress_Update() leaving ", p->stress->uMatrix.daa, p->dof_per_node, 12);
     printf(" Leaving Stress_Update() \n");
#endif
}

void
Plastic_Deform(p, H, R, eff_pl_strain, stress, stress_dev,
               strain_incr, E, fy, G, A, iNo_iter_step, dof, ii, kk)
ARRAY                              *p;
double            *H, *R, fy, E, G, A;
double         **stress, **stress_dev;
double                  **strain_incr;
double                 *eff_pl_strain;
int        iNo_iter_step, dof, ii, kk;
{
double                           temp;
double                         Lambda;
double             eff_pl_strain_incr;
double           temp1, effect_stress;
int                           i, j, k;

#ifdef DEBUG
    printf(" enter Plastic_Deform() \n");
#endif

  /* Step 5 Compute effective incremental */
  /*        plastic strain within each    */
  /*        sub-incrementation            */
              
  /* Estimate H' */
  if(!strcmp(p->material_name, "ELASTIC_PERFECTLY_PLASTIC")) {
     *H = 0.0;
  }else {
     if(!strcmp(p->LC_ptr->name, "Ramberg-Osgood")) {
         effect_stress = A*sqrt(3.0/2.0);

#ifdef DEBUG
    printf(" before H = %le \n", *H);
    printf(" fy = %lf \n", fy);
    printf(" E  = %lf \n", E);
    printf(" effect_stress= %le \n", effect_stress);
#endif
         Load_Curve(p, H, effect_stress, E,fy);
     }
     if(!strcmp(p->LC_ptr->name, "Bi-Linear")) 
         *H = p->LC_ptr->H[kk-1];
  }
  /* calculate the effective plastic strain incremental */

  temp1 = sqrt(2.0/3.0)*(*H + 3.0*G);
  eff_pl_strain_incr  = A/temp1 - (*R)/temp1;
 *eff_pl_strain      += eff_pl_strain_incr;

#ifdef DEBUG3
     printf("========== In Plastic_Deform() : H = %le\n", *H);
     printf("========== In Plastic_Deform() : eff_pl_strain_incr = %le\n", eff_pl_strain_incr);
#endif

  /* Step 6: Compute Lambda  and pl_strain_incr       */
  /*         Lambda = sqrt(3/2)*eff_pl_strain_incr/A  */
  /* plastic strain incr is now stored in strain_incr */

  Lambda = sqrt(3.0/2.0)*eff_pl_strain_incr;
  for(i = 1; i <= dof; i++) {
      strain_incr[i-1][0] = Lambda*stress_dev[i-1][0]/A;
      p->strain_pl->uMatrix.daa[i-1][kk-1] += strain_incr[i-1][0];
  }
  /* Step 7:  Update back stress and Stress */

  /* beta = 0 for kinematic hardening */
  /* beta = 1 for isotropic hardening */

  if( ABS(p->LC_ptr->beta) < 1E-10){
      for(i = 1; i <= dof; i++) {
          temp = (*H)*strain_incr[i-1][0]*2.0/3.0;
          stress[i-1][0] = stress[i-1][0]-2.0*G*strain_incr[i-1][0] + temp;
          p->LC_ptr->back_stress[i-1][kk-1] += temp;
      }
  }else {

  /* Step 8:  Update stress */
  for(i = 1; i <= dof; i++) 
      stress[i-1][0] = stress[i-1][0]-2.0*G*strain_incr[i-1][0];
 }

  if(ii == iNo_iter_step) 
     for(i = 1; i <= dof; i++) 
      p->stress->uMatrix.daa[i-1][kk-1] = stress[i-1][0];

  /* Step 9:  Update R */
#ifdef DEBUG
    printf(" before R = %lf \n", *R);
    printf(" H       = %le \n",  *H);
    printf(" beta    = %lf \n", p->LC_ptr->beta);
    printf(" eff_pl_strain_incr = %le \n", eff_pl_strain_incr);
#endif
    if( ABS(p->LC_ptr->beta -1.0) < 1E-10)
       (*R) += sqrt(2.0/3.0)*(*H)*eff_pl_strain_incr;

#ifdef DEBUG
    printf(" after R = %lf \n", *R);
#endif

#ifdef DEBUG
    printf(" Leaving Plastic_Deform() \n");
#endif

}

double **B_MATRIX_4Node(B_matrix, p, shp, z_coord)
double       **B_matrix;
ARRAY                *p;
double            **shp;
double          z_coord;
{
int       i,j, k, ii, n;
double                h;


/* ======================================================= */
/* B_matrix : strain = Transpose(B_matrix)* nodal_displ.   */
/* ======================================================= */

#ifdef DEBUG 
     printf(" enter B_MATRIX_4Node() \n");
#endif

      /* B_matrix = [B', B"] */

      for(j = 1; j <= p->nodes_per_elmt; j++) {
          k = p->dof_per_node*(j-1)-1;

      /* ------------------------------------------ */
      /* Bi' mattrix is independent of z-coordinate */
      /* Bi' mattrix is estimated first here        */
      /* Bi'  = []5x3                               */
      /* ------------------------------------------ */

          B_matrix[0][k+1] = shp[0][j-1]; 
          B_matrix[0][k+2] = B_matrix[0][k+3] = 0.0;

          B_matrix[1][k+2] = shp[1][j-1]; 
          B_matrix[1][k+1] = B_matrix[1][k+3] = 0.0;

          B_matrix[2][k+1] = shp[1][j-1]; 
          B_matrix[2][k+2] = shp[0][j-1]; 
          B_matrix[2][k+3] = 0.0;

          B_matrix[3][k+3] = shp[1][j-1]; 
          B_matrix[3][k+1] = B_matrix[3][k+2] = 0.0;

          B_matrix[4][k+3] = shp[0][j-1]; 
          B_matrix[4][k+1] = B_