elmt_shell_4n.c

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

#ifdef DEBUG
       printf("*** In elmt_shell_4nodes_implicit() : start case MASS\n");
       printf("                : Density = %14.4e\n", density.value);
       printf("                : shell_thickness = %8.2f\n", h);
       printf("                : Jac             = %8.2f\n", Jac);
#endif
   /*====================================================*/
   /*  CALCULATING MASS MATRIX IN CO_COORDINATE SYSTEM   */ 
   /*====================================================*/

      Shell_4Node_Mass(p, p->stiff, co_coord, density.value, h);

#ifdef DEBUG
       printf("*** In elmt_shell() : end case MASS\n");
#endif
             break;
        default:
             break;
    }

    MatrixFreeIndirectDouble(co_coord, p->no_dimen);
    MatrixFreeIndirectDouble(Direction_Matrix,6);

#ifdef DEBUG
       printf("*** leaving elmt_shell() \n");
#endif

    return(p);
}


/* ============================================*/
/* Equivalent Load due to distributed pressure */
/* for shell Belytschko element                */
/* in local coordinate                         */
/* ============================================*/

ARRAY *sld04(p,task)
ARRAY *p;
int task;
{
ELEMENT_LOADS               *elsptr;
ELOAD_LIB                      *elp;

double                  *nodal_load;
double      Jac, x31, y31, x42, y42; 
double px,py,pz, mx,my,mz, bx,by,bz;

int  i,j,k, ii, jj, kk, n, n1,n2,n3;

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

    /* Initialize total load */

    nodal_load = dVectorAlloc(p->size_of_stiff);

    for(i =1; i <= p->size_of_stiff ;i++)
        nodal_load[i-1] = 0.0;

#ifdef DEBUG
      printf(" **** In sld04(): begin switch()  \n");
#endif

    switch(task){
       case PRESSLD:
#ifdef DEBUG
      printf(" **** In sld04(): case PRESSLD\n");
#endif

       /* ==========================================*/
       /* For distrbuted surface loading:           */
       /* PRESSURE OR TRACTION                      */
       /* nodal equivalent force is calculated by : */
       /* force = integral N^T*traction/pressure dA */
       /* ==========================================*/
       
       /*================================================*/
       /* use one-point integration rule in surface area */
       /* Ni = 1/4.0 [N] = 1/4 [I, I, I, I]              */
       /* [I] = 5x5 unit_matrix                          */
       /* Let traction = [p]5x1                          */
       /* Integ [N]^t*traction dA  = (1/4)*[p]*Jac*4.0   */
       /*                                  [p]           */
       /*                                  [p]           */
       /*                                  [p]           */
       /*                                  [p]           */
       /*================================================*/
         elsptr  =  p->elmt_load_ptr;
         for (j=1; j<= elsptr->no_loads_faces;j++) {
              elp = &elsptr->elib_ptr[j-1];
       
              x31 =  p->coord[0][2].value - p->coord[0][0].value;
              y31 =  p->coord[1][2].value - p->coord[1][0].value;

              x42 =  p->coord[0][3].value - p->coord[0][1].value;
              y42 =  p->coord[1][3].value - p->coord[1][1].value;

              Jac = 0.5*(x31*y42-x42*y31);
 
              for(i = 1; i <= p->size_of_stiff; i++) {
                  n = i/p->dof_per_node;
                  k = i - n*p->dof_per_node;
                  nodal_load[i-1] = 0.25*elp->traction[k-1].value*4.0*Jac; 
              }
         }
       break;
    default:
       break;
    }

    for(i = 1; i <= p->size_of_stiff; i++) {
       p->equiv_nodal_load->uMatrix.daa[i-1][0] = nodal_load[i-1];
    }
    
    free(nodal_load);

#ifdef DEBUG
      printf("**** leaving sld04(): \n");
#endif
    return(p);
}


/* ============================*/
/* Material matrix             */
/* ============================*/

double **MATER_MAT_SHELL(m1, E, nu)
double  **m1;
double E, nu;
{
double           temp;
double       kk = 1.2;       
    
#ifdef DEBUG
     printf(" enter MATER_MAT_SHELL() \n");
     printf(" E = %lf , nu = %lf \n", E, nu);
#endif
    /* Stress : stress_x, stress_y,  stress_xy,  stress_yz,  stress_zx  */
    /* Strain : strain_x, strain_y, 2strain_xy, 2strain_yz, 2strain_zx  */

    
    temp =  E/(1-nu*nu);

    m1[0][0]= m1[1][1] = temp;
    m1[0][1]= m1[1][0] = nu*temp;
    m1[2][2]= E/2.0/(1.0+nu);
    m1[3][3] = m1[4][4] = E/2.0/(1.0+nu)/kk;

    m1[0][2] = m1[2][0] = m1[1][2] = m1[2][1] = 0.0;
    m1[3][0] = m1[3][1] = m1[3][2] = m1[3][4] = 0.0;
    m1[4][0] = m1[4][1] = m1[4][2] = m1[4][3] = 0.0;


#ifdef DEBUG
     dMatrixPrint("m1", m1, 5, 5);
     printf("\n leaving MATER_MAT_SHELL() \n");
#endif
  return (m1);
}

void MATER_SHELL_UPDATE(p, co_coord, E, nu, integ_pt)
ARRAY                              *p;
double                     **co_coord;
double                          E, nu;
int                          integ_pt; /* integration point */
{
double             **m1, **mater_temp;
double **stress_dev, **Norm_Transpose;
double         A_2, mean_stress, R, A;
double                           temp;
int            i, j, k, ii, jj, kk, n;
int                      dof, size, G;

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

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

#ifdef DEBUG
    printf(" In MATER_SHELL_UPDATE(): calculating MATER_MAT_SHELL() \n");
#endif
      dof             = p->dof_per_node;
      size            = p->size_of_stiff;
      ii              = integ_pt;
      G               = E/2.0/(1.0+nu);


#ifdef DEBUG
    printf(" In MATER_SHELL_UPDATE(): update p->material = %s \n", p->material_name);
    printf(" In MATER_SHELL_UPDATE(): integ_pts          = %d \n", integ_pt);
    printf(" In MATER_SHELL_UPDATE(): dof                = %d \n", dof);
    printf(" In MATER_SHELL_UPDATE(): deformation state  = %d \n", p->elmt_state);
#endif

  /* Elastic material */

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

     p->mater_matrix->uMatrix.daa 
     = MATER_MAT_SHELL(p->mater_matrix->uMatrix.daa,E, nu);
   }

  /* Elastic Perfectly Plastic Materials */


  if(p->material_name != NULL &&
     !strcmp(p->material_name,"ELASTIC_PERFECTLY_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);

     /* Elastic Deformation */

     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];
               else
                 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];