elmt_shell_8n.c

有限元分析源代码

         break;
      }

#ifdef DEBUG
       printf("*** In elmt_shell() : print equiv_nodal_load\n\n");
       MatrixPrintIndirectDouble(p->equiv_nodal_load);
       printf("*** In elmt_shell() : end case EQUIV_NODAL_LOAD\n");
#endif

             break;

        case STRESS_UPDATE:

        /* update stress for given displacement */
        /* and displacement incremental         */
             break;
        case STRESS:
        case STRESS_LOAD:
        case LOAD_MATRIX:
        /* ====================================== */
        /* Form internal nodal load vector due to */
        /* stress at previous load step           */
        /* ====================================== */

#ifdef DEBUG
       printf("****** In elmt_shell_8nodes_implicit() : \n");
       printf("       enter cases: STRESS_UPDATE, LOAD_MATRIX, STRESS_LOAD \n");     
       printf(" elemt_state = %d \n", p->elmt_state);
#endif
          nodal_load      = MatrixAllocIndirectDouble(p->size_of_stiff, 1); 
          nodal_load_temp = MatrixAllocIndirectDouble(p->size_of_stiff, 1); 

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

          EE  = E.value;

        /* =================================================*/
        /* nodal load due to the inital stress or stress at */
        /* previous step                                    */
        /* =================================================*/

      /* Integration loop over thickness */

          integ_coord = dVectorAlloc(no_integ_pts+1);
          weight      = dVectorAlloc(no_integ_pts+1);

          gauss(integ_coord,weight,no_integ_pts);

          for(ii = 1; ii <= no_integ_pts; ii++) {
              nodal_load_temp 
              = Shell_Nodal_Load_8Node(nodal_load_temp, p, co_coord,
                             integ_coord[ii], ii, E.value, nu, isw);
              for(i = 1; i<= p->size_of_stiff; i++) {
                  nodal_load[i-1][0] += nodal_load_temp[i-1][0]*weight[ii];  
              }

          } /* gaussian integration ends */

          free((char *) integ_coord);
          free((char *) weight);

      size = p->size_of_stiff;	

      for(i= 1; i<= p->size_of_stiff; i++) {
          p->nodal_loads[i-1].value = nodal_load[i-1][0];  
      }

       /***************************************************/
       /* Transform nodal_load vector from local lamina   */
       /* coordinate system to global coordinate system   */
       /***************************************************/

      /* =====================================*/
      /* Calculate the inverse of T_matrix    */
      /*     T_matrix inverse [T]^-1 =        */
      /*        = T_matrix trnspose : [T]^t   */
      /* =====================================*/

      /* store the direction matrix : T_matrix */

   if( isw == LOAD_MATRIX ) {
             e1_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);
             e2_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);
             e3_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);

             Lamina_Sys_8node(p, e1_ptr, e2_ptr, e3_ptr);

    Direction_Matrix = MatrixAllocIndirectDouble(6,6);
    for ( i = 1; i <= 6; i++) {
        if(i <= 3) {
          Direction_Matrix[0][i-1] = e1_ptr[i-1][0];
          Direction_Matrix[1][i-1] = e2_ptr[i-1][0];
          Direction_Matrix[2][i-1] = e3_ptr[i-1][0];
          Direction_Matrix[3][i-1] = 0.0;
          Direction_Matrix[4][i-1] = 0.0;
          Direction_Matrix[5][i-1] = 0.0;
        }
        else {
          j = i - 3;
          Direction_Matrix[0][i-1] = 0.0;
          Direction_Matrix[1][i-1] = 0.0;
          Direction_Matrix[2][i-1] = 0.0;
          Direction_Matrix[3][i-1] = e1_ptr[j-1][0];
          Direction_Matrix[4][i-1] = e2_ptr[j-1][0];
          Direction_Matrix[5][i-1] = e3_ptr[j-1][0];
        }
    }

    T_matrix     = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
    T_Transpose  = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);

      for (i = 1; i <= p->size_of_stiff; i++)  {
           for(j = 1; j <= p->nodes_per_elmt; j++) {
               for( k = 1; k <= p->dof_per_node; k++) {
                    n1  = p->dof_per_node*(j-1);
                    n2  = p->dof_per_node*j;
                    n   = n1 + k;
                    if(i > n1 && i <= n2) {
                       ii = i - n1;
                       T_matrix[i-1][n-1] = Direction_Matrix[ii-1][k-1];
                    }
                    else 
                       T_matrix[i-1][n-1] = 0.0;
                }
           }
       }

      for( i = 1; i <= p->size_of_stiff; i++) {
         for( j = 1; j <= p->size_of_stiff; j++) 
           T_Transpose[i-1][j-1] = T_matrix[j-1][i-1];
      }
     
      size = p->size_of_stiff;	
      nodal_load_temp = dMatrixMultRep(nodal_load_temp,
                        T_Transpose, size, size, nodal_load, size,1);

      for(i= 1; i<= p->size_of_stiff; i++) {
          p->nodal_loads[i-1].value = nodal_load_temp[i-1][0];  
      }

      MatrixFreeIndirectDouble(e1_ptr,3);
      MatrixFreeIndirectDouble(e2_ptr,3);
      MatrixFreeIndirectDouble(e3_ptr,3);
      MatrixFreeIndirectDouble(Direction_Matrix,6);
      MatrixFreeIndirectDouble(T_matrix,p->size_of_stiff);
      MatrixFreeIndirectDouble(T_Transpose,p->size_of_stiff);
    }

      MatrixFreeIndirectDouble(nodal_load,p->size_of_stiff);
      MatrixFreeIndirectDouble(nodal_load_temp,p->size_of_stiff);


   /* ------------NODAL LOAD UNITS ------------------------*/

      switch( UNITS_SWITCH ) {
        case ON:
           if(UNITS_TYPE == SI || UNITS_TYPE == SI_US ) {
               dp_force  = DefaultUnits("N");
               dp_length = DefaultUnits("m");
           }
           else {
               dp_force  = DefaultUnits("lbf");
               dp_length = DefaultUnits("in");
           }

           /* node no 1 */
           UnitsCopy( p->nodal_loads[0].dimen, dp_force );
           UnitsCopy( p->nodal_loads[1].dimen, dp_force );
           UnitsCopy( p->nodal_loads[2].dimen, dp_force );
           UnitsMultRep( p->nodal_loads[3].dimen, dp_force, dp_length );
           UnitsCopy( p->nodal_loads[4].dimen, p->nodal_loads[3].dimen );

           /* node no > 1 */
           for(i = 2; i <= p->nodes_per_elmt; i++) {    
               for(j = 1; j <= p->dof_per_node; j++) {
                   k = p->dof_per_node*(i-1)+j;
                   if(j <= 3) 
                      UnitsCopy( p->nodal_loads[k-1].dimen, p->nodal_loads[0].dimen );
                   else 
                      UnitsCopy( p->nodal_loads[k-1].dimen, p->nodal_loads[3].dimen );
               }
           }
           free((char *) dp_force->units_name);
           free((char *) dp_force);
           free((char *) dp_length->units_name);
           free((char *) dp_length);

         break;
      case OFF:
         break;
      default:
         break;
    }

#ifdef DEBUG
printf("*** In elmt_shell_8nodes_implicit() : end case LOAD_MATRIX, STRESS_UPDATE, STRESS_LOAD \n");
#endif
	     break;

        case MASS_MATRIX:  /* form mass matrix */

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

      Shell_8Node_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);

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

    return(p);
}


/* ============================================*/
/* Equivalent Load due to distributed pressure */
/* for 8 node shell element                    */
/* ============================================*/

ARRAY *sld108(p,task)
ARRAY *p;
int task;
{
#ifdef DEBUG
      printf("**** enter sld108(): \n");
#endif
       /* add details */

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


void Stress_Update_8Node(p, co_coord, B1_matrix, integ_pt)
ARRAY                   *p;
double          **co_coord;
double          **B1_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_8Node(): \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); /* stress incremental            */

#ifdef DEBUG
      printf(" In Stress_Update_8Node(): 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, 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);
     dMatrixPrint("strain_incr", strain_incr, dof, 1);
     dMatrixPrint("stress", stress, dof, 1);
#endif

     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;