elmt_shell_4n.c

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

           }

       /* Step 4: Calculate [Cep]                         */ 
       /*      [Cep] = [C] -[C]*[n][n]^t*[C]/(2G+H*2/3) */
       /*      [Cep] = [C] -[C]*[n][n]^t*2G/(2G+H*2/3)  */
              
             p->mater_matrix->uMatrix.daa = dMatrixMultRep(p->mater_matrix->uMatrix.daa,
                                            stress_dev, dof,1, Norm_Transpose, 1, dof);

             mater_temp = dMatrixMultRep(mater_temp, m1, dof, dof,
                          p->mater_matrix->uMatrix.daa, dof, dof);

/*
             p->mater_matrix->uMatrix.daa = dMatrixMultRep(p->mater_matrix->uMatrix.daa,
                                            mater_temp, dof, dof, m1, dof, dof);
*/
             for(i = 1; i <= dof; i++) {
                for(j = 1; j <= dof; j++) {
                    p->mater_matrix->uMatrix.daa[i-1][j-1] 
                    = m1[i-1][j-1]- p->mater_matrix->uMatrix.daa[i-1][j-1]
                      *2.0*G/(2.0*G+p->LC_ptr->H[ii-1]*2.0/3.0);
                }
             }
         
         } /* end of else : plastic deformation */

       break;
       default:
         printf(" In MATER_SHELL_UPDATE():  elmt_no  = %d\n", p->elmt_no);
         printf(" elmt_state = 0 : Elastic_deformation \n");
         printf(" elmt_state = 1 : Elastic-plastic_deformation \n");
         printf(" elmt_state = %d: p->elmt_state \n");
         FatalError(" Unknown elmt state ",(char *)NULL);
       break;
     }
     MatrixFreeIndirectDouble(stress_dev, p->dof_per_node);
     MatrixFreeIndirectDouble(Norm_Transpose, 1);
     MatrixFreeIndirectDouble(m1, p->dof_per_node);
     MatrixFreeIndirectDouble(mater_temp, p->dof_per_node);
  }

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


void DISPL_UPDATE(p)
ARRAY    *p;
{
int i,j, k;

   for(i = 1; i <= p->dof_per_node; i++) {
      for(j = 1; j <= p->nodes_per_elmt; j++) {
          p->displ->uMatrix.daa[i-1][j-1] += p->displ_incr->uMatrix.daa[i-1][j-1];
      }
   } 

}

void Stress_Update(p, co_coord, B_matrix, integ_pt)
ARRAY                   *p;
double          **co_coord;
double          **B_matrix;
int               integ_pt;
{
double                      H, R, fy, E, Et, nu;
double          **displ_incr,**m1,**strain_incr;
double          **back_stress_incr, **strain_pl;
double    **stress, **stress_dev, **stress_incr;
double     temp, G, A_2, A, Lambda, mean_stress;
double         eff_pl_strainTemp, eff_pl_strain;
double       Beta1 = 0.05, temp1, effect_stress;
int                    iNo_iter_step, dof, size;
int                      i, j, k, ii, jj, kk, n;
DIMENSIONS                               *dimen;

  /* Given i -state displ, stress mater_matrix, [Cep], */
  /* and strain,  calculate i+1 state displ_incr       */
  /* these i+1 state stress can be used for            */
  /* i+1 state [Cep] matrix                             */


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

  kk   = integ_pt;
  dof  = p->dof_per_node;
  size = p->size_of_stiff;
  E    = p->work_material[0].value;
  nu   = p->work_material[4].value;
  fy   = p->work_material[2].value;
  G    = E/(1.0+nu)/2.0;

  m1         = MatrixAllocIndirectDouble(p->dof_per_node, p->dof_per_node);
  m1         = MATER_MAT_SHELL(m1, E, nu);

  /* Step 1: calculate strain_incremental  */
  /*         from displacement incremental */
  /*         at given integration point    */

   displ_incr       = MatrixAllocIndirectDouble(p->size_of_stiff, 1);
   strain_incr      = MatrixAllocIndirectDouble(dof,1);
   strain_pl        = MatrixAllocIndirectDouble(dof,1);
   stress           = MatrixAllocIndirectDouble(dof,1); /* stress or stress incremental   */
   stress_dev       = MatrixAllocIndirectDouble(dof,1); /* deviatric component of stress  */
   stress_incr      = MatrixAllocIndirectDouble(dof,1); /* stress incremental             */
   back_stress_incr = MatrixAllocIndirectDouble(dof,1); /* back_stress stress incremental */

#ifdef DEBUG
      printf(" In Stress_Update(): materials_name = %s \n",p->material_name);
#endif

  /* ================= */
  /* Elastic material  */
  /* ================= */

  if(p->material_name != NULL && !strcmp(p->material_name,"ELASTIC")) {

     /* Transfer p->displ into vector form */
     /* -----------------------------------*/
     /* There is no needed to calculate    */
     /* incrementally for the case of      */
     /* linear elasticity                  */
     /* -----------------------------------*/

     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->uMatrix.daa[i-1][j-1];
         }
     }

     strain_incr = dMatrixMultRep(strain_incr, B_matrix, dof,
                                  size, displ_incr, size, 1); 
     stress      = dMatrixMultRep(stress, m1, dof, dof,
                                  strain_incr, dof, 1);

     for (i = 1; i <= dof; i++){
         p->stress->uMatrix.daa[i-1][kk-1] += stress[i-1][0];
     }
     SaveRespondBuffer(p, kk);
  }

  /* ====================================*/
  /* Elastic Plastic & Elastic Perfectly */
  /* Plastic Materials                   */
  /* ====================================*/

  if(p->material_name != NULL &&
      (!strcmp(p->material_name, "ELASTIC_PLASTIC") ||
       !strcmp(p->material_name, "ELASTIC_PERFECTLY_PLASTIC")) ) {

     /* Transfer p->displ_incr into vector form */

     switch(p->elmt_state) {
       case 0 :       /* elastic state */
         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->uMatrix.daa[i-1][j-1];
             }
         }
       break;
       case 1 :       /* inelastic state */
         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];
             }
         }
       break;
       default:
         printf("**** In Stress_Update(): 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, B_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);
     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
      for (i = 1; i <= dof; i++)
         printf(" total stress[%d] = %lf \n", i, stress[i-1][0]);
#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;
              }
          }

          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,