elmt_shell_8n.c

有限元分析源代码

                 displ_incr[k-1][0] = p->displ_incr->uMatrix.daa[i-1][j-1];
             }
         }
       break;
       default:
         printf("**** In Stress_Update_8Node(): elmt_no = %d\n", p->elmt_no);
         printf(" elmt_state = %d: p->elmt_state \n");
         FatalError("*****Undefine element state ()", (char *) NULL);
       break;
     }

     /* Calculate stress increment */

     strain_incr = dMatrixMultRep(strain_incr, B1_matrix, dof,
                                  size, displ_incr, size, 1); 

     stress      = dMatrixMultRep(stress, m1, dof, dof,
                                  strain_incr, dof, 1);

#ifdef DEBUG
     dMatrixPrint("displ_incr", displ_incr, size, 1);
#endif

#ifdef DEBUG
     dMatrixPrint("stress", stress, dof, 1);
     dMatrixPrint("strain", strain_incr, dof, 1);
     printf(" \n stress before incremented \n");
     for(i = 1; i <= dof; i++)
         printf(" p->stress->uMatrix.daa[%d][%d]= %lf \n",
               i, kk, p->stress->uMatrix.daa[i-1][kk-1]);
#endif

     switch(p->elmt_state) {

       /* Elastic deformation */
       case 0: 
         for(i = 1; i <= dof; i++)
             p->stress->uMatrix.daa[i-1][kk-1] = stress[i-1][0];

         SaveRespondBuffer(p, kk);

       break;
       case 1:
         
       /* plastic deformation */
       /* Step 1: Calculate the spherical stress and */
       /*         deviatric stress of trial stress   */
       /*         also the radiaus square A_2        */

#ifdef DEBUG
     for (i = 1; i <= dof; i++)
         printf("incrmental stress[%d] = %lf \n", i, stress[i-1][0]);
#endif
         for (i = 1; i <= dof; i++) {
           stress_incr[i-1][0] = stress[i-1][0];
           stress[i-1][0] += p->stress->uMatrix.daa[i-1][kk-1];
         }

#ifdef DEBUG
      printf(" element no  = %d \n", p->elmt_no);
      for (i = 1; i <= dof; i++){
         printf(" total stress[%d] = %lf \n", i, stress[i-1][0]);
         printf("incrmental stress[%d] = %lf \n", i, stress_incr[i-1][0]);
         printf("previous stress[%d] = %lf \n", i, p->stress->uMatrix.daa[i-1][kk-1]);
      }
#endif
         /* Note : sigma_z == 0 for shell elmt */
         
         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);
         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;
              }
          } else {
              for(i = 1; i <= dof; i++) 
                  back_stress_incr[i-1][0] = 0.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.0*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
          /* 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,B1_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_8Node(): 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_8Node() leaving ",
                    p->stress->uMatrix.daa, p->dof_per_node, 12);
     printf(" Leaving Stress_Update_8Node() \n");
#endif
}

double **B_MATRIX_8Node(B_matrix, p, shp, z_coord, J_inverse)
double        **B_matrix;
ARRAY                 *p;
double             **shp;
double           z_coord;
double       **J_inverse;
{
int       i,j, k, ii, n;
double                h;
static double     **a = NULL;
static double     **b = NULL;
static double     **c = NULL;
static double     **d = NULL;
static double     **e = NULL;
static double     **g = NULL;
static double     **e1_ptr = NULL;
static double     **e2_ptr = NULL;
static double     **e3_ptr = NULL;

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

       a = MatrixAllocIndirectDouble(p->nodes_per_elmt,1);
       b = MatrixAllocIndirectDouble(p->nodes_per_elmt,1);
       c = MatrixAllocIndirectDouble(p->nodes_per_elmt,1);
       d = MatrixAllocIndirectDouble(p->nodes_per_elmt,1);
       e = MatrixAllocIndirectDouble(p->nodes_per_elmt,1);
       g = MatrixAllocIndirectDouble(p->nodes_per_elmt,1);

       e1_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);
       e2_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);
       e3_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);

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

       for(i = 1; i <= p->nodes_per_elmt; i++) {
           a[i-1][0] = J_inverse[0][0]*shp[0][i-1] + J_inverse[0][1]*shp[1][i-1] ;  
           b[i-1][0] = J_inverse[1][0]*shp[0][i-1] + J_inverse[1][1]*shp[1][i-1] ;  
           c[i-1][0] = J_inverse[2][0]*shp[0][i-1] + J_inverse[2][1]*shp[1][i-1] ;  

           d[i-1][0] = h*0.5*(a[i-1][0]*z_coord + J_inverse[0][2]*shp[2][i-1]);  
           e[i-1][0] = h*0.5*(b[i-1][0]*z_coord + J_inverse[1][2]*shp[2][i-1]);  
           g[i-1][0] = h*0.5*(d[i-1][0]*z_coord + J_inverse[2][2]*shp[2][i-1]);  
       }

#ifdef DEBUG
       for(i = 1; i <= p->nodes_per_elmt; i++) {
           printf(" a[%d] = %lf \n", i, a[i-1][0]);
           printf(" b[%d] = %lf \n", i, b[i-1][0]);
           printf(" c[%d] = %lf \n", i, c[i-1][0]);

           printf(" d[%d] = %lf \n", i, d[i-1][0]);
           printf(" g[%d] = %lf \n", i, g[i-1][0]);
           printf(" e[%d] = %lf \n", i, e[i-1][0]);
       }
#endif
       Lamina_Sys_8node(p, e1_ptr, e2_ptr, e3_ptr);
       
      for(i = 1; i <= p->dof_per_node + 1; i++) 
          for(j = 1; j <= p->size_of_stiff; j++) 
              B_matrix[i-1][j-1] = 0.0;

      /* ------------------------------------------ */
      /* Bi'  = []6x3                               */
      /* ------------------------------------------ */

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

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