elmt_shell_4n.c

有限元程序

                 stress_dev[i-1][0] = p->stress->uMatrix.daa[i-1][ii-1] - p->LC_ptr->back_stress[i-1][ii-1];

               A_2 += stress_dev[i-1][0]*stress_dev[i-1][0]; 
            }

       /* Step 2: comparision                        */

            R = p->LC_ptr->R[ii-1];

            if(A_2 <= R*R) { /* ELASTIC DEFORMATION */
               p->mater_matrix->uMatrix.daa 
               = MATER_MAT_SHELL(p->mater_matrix->uMatrix.daa,E, nu);
            }
            else {           /* PLASTIC DEFORMATION */
        
              m1  = MATER_MAT_SHELL(m1, E, nu);

              /* Step 3: normal direction on yield surface */ 
              /*         stored temporarily on stress_dev  */
              /* Norm_Direc */

              for(i = 1; i <= dof; i++) {
                 stress_dev[i-1][0] /= R; 
                 Norm_Transpose[0][i-1] = stress_dev[i-1][0];
              }
                 
              /* Step 4: Calculate [Cep]                   */ 
              /*         [Cep] = [C] -[C]*[n][n]^t         */
              
               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);

               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] - mater_temp[i-1][j-1];
                  }
               }
         
            } /* 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);
  }

  /* Elastic Plastic Materials      */
  /* with combined strain hardening */

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

     m1              = MatrixAllocIndirectDouble(p->dof_per_node, p->dof_per_node);
     stress_dev      = MatrixAllocIndirectDouble(dof, 1);
     Norm_Transpose  = MatrixAllocIndirectDouble(1,dof);
     mater_temp      = MatrixAllocIndirectDouble(dof, dof);
    
     switch(p->elmt_state) {
       /* Elastic deformation */
       case 0: 
         p->mater_matrix->uMatrix.daa 
         = MATER_MAT_SHELL(p->mater_matrix->uMatrix.daa, E, nu);
       break;
       case 1:
       /* plastic deformation */

       /* Step 1: Calculate the spherical stress and */
       /*         deviatric stress of trial stress   */
       /*         also the radiaus square A_2        */

         mean_stress = (p->stress->uMatrix.daa[0][ii-1]- p->LC_ptr->back_stress[0][ii-1]
                       +p->stress->uMatrix.daa[1][ii-1]- p->LC_ptr->back_stress[1][ii-1])/3.0;
         
         A_2 = 0.0;
         for(i = 1; i <= dof; i++) {
             if(i <= 2)           
               stress_dev[i-1][0] = p->stress->uMatrix.daa[i-1][ii-1]
                                  - mean_stress - p->LC_ptr->back_stress[i-1][ii-1];
             A_2 += stress_dev[i-1][0]*stress_dev[i-1][0]; 
         }
         A = sqrt(A_2);

       /* Step 2: comparision   */

         R = p->LC_ptr->R[ii-1];

         if(A_2 <= R*R) { /* ELASTIC DEFORMATION */
            p->mater_matrix->uMatrix.daa 
            = dMatrixCopyRep(p->mater_matrix->uMatrix.daa,
              m1, p->dof_per_node, p->dof_per_node);
         }
         else {  /* PLASTIC DEFORMATION */
        
           m1  = MATER_MAT_SHELL(m1, E, nu);

           /* Step 3: normal direction on yield surface */ 
           /*         stored temporarily on stress_dev  */
           for(i = 1; i <= dof; i++) {
               stress_dev[i-1][0] /= R;  /* Norm_Direc */
               Norm_Transpose[0][i-1] = stress_dev[i-1][0];
           }

       /* 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];