elmt_shell_4n.c

利用语言编写的有限元分析软件

                                            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(UNITS_TYPE) {
           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_matrix[4][k+2] = 0.0;

          h  = p->work_section[11].value;    /* thickness of the shell */

      /* ----------------------------------*/
      /*  Calculate Bi" matrix             */
      /*  Bi" = [ ] 5x2                    */
      /* ----------------------------------*/

          B_matrix[0][k+4] =  0.0;
          B_matrix[0][k+5] =  shp[0][j-1]*h*0.5*z_coord; 

          B_matrix[1][k+4] = -shp[1][j-1]*h*0.5*z_coord; 
          B_matrix[1][k+5] =  0.0;

          B_matrix[2][k+4] = -shp[0][j-1]*h*0.5*z_coord; 
          B_matrix[2][k+5] =  shp[1][j-1]*h*0.5*z_coord; 

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

          B_matrix[4][k+4] =  0.0;
          B_matrix[4][k+5] =  shp[2][j-1];
      }

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

      return(B_matrix);
}


void Load_Curve(p, H, stress, E, fy)
ARRAY         *p;
double    stress;
double        *H;
double     E, fy;
{
double temp1, temp2;
    
#ifdef DEBUG
    printf(" enter Loac_Curve() \n");
#endif
    temp1 = stress/fy;
    temp2 = p->LC_ptr->n -1.0;

    *H = E/(p->LC_ptr->n*p->LC_ptr->alpha*pow(temp1, temp2));

#ifdef DEBUG
    printf(" Inside Load_Curve(): H = %le\n", *H);
    printf(" Leaving Load_Curve() \n");
#endif

 /* =================================================*/
 /* for Ramberg-Osgood stress-strain relations       */
 /* : strain/strain0 =                               */
 /*    stress/stress0+alpha*(stress/stress0)^n       */
 /*  pl_strain = strain0a*alpha*(stress/stress0)^n   */
 /* =================================================*/
 /* H = d(stress)/d(pl_strain) =                     */
 /*    = E/(alpha*n)*(stress/stress0)^(1-n)          */
 /* =================================================*/
}


void BB_Vector(co_coord, BB1, BB2) 
double **co_coord;
double *BB1, *BB2;
{
double co_x31, co_y31, co_z31;
double co_x42, co_y42, co_z42;
int i, j, k;

#ifdef DEBUG
      printf(" enter BB_Vector() \n");
      dMatrixPrint("co_coord", co_coord, 3, 4);
#endif

      /* BB_Vector */

      co_x31 =  co_coord[0][2] - co_coord[0][0];
      co_y31 =  co_coord[1][2] - co_coord[1][0];
      co_z31 =  co_coord[2][2] - co_coord[2][0];

      co_x42 =  co_coord[0][3] - co_coord[0][1];
      co_y42 =  co_coord[1][3] - co_coord[1][1];

      BB1[0] = -co_y42/4.0;
      BB1[1] =  co_y31/4.0;
      BB1[2] =  co_y42/4.0;
      BB1[3] = -co_y31/4.0;

      BB2[0] =  co_x42/4.0;
      BB2[1] = -co_x31/4.0;
      BB2[2] = -co_x42/4.0;
      BB2[3] =  co_x31/4.0;

 /* in Belyschenko's paper                    */
 /* B1I = (1/2) [...]  instead of (1/4) [...] */
 /* B2I = (1/2) [...]  instead of (1/4) [...] */
    

#ifdef DEBUG
  for(i = 1; i <= 4; i++)
      printf(" BB1[%d] = %lf, BB2[%d] = %lf \n", i, BB1[i-1], i, BB2[i-1]);
      printf(" leaving BB_Vector() \n");
#endif
  
}


double **Shell_Nodal_Load_Plane(nodal_load, p, co_coord, z_coord, z_integ_pt, EE, nu, task)
double                 **nodal_load;