elmt_shell_4n.c

有限元分析源代码

/*
 *  ============================================================================= 
 *  ALADDIN Version 1.0 :
 *      elmt_shell_4n.c : STATIC/DYNAMIC SHELL_4N Element
 *                                                                     
 *  Copyright (C) 1995 by Mark Austin, Xiaoguang Chen, and Wane-Jang Lin
 *  Institute for Systems Research,                                           
 *  University of Maryland, College Park, MD 20742                                   
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software for any on any computer system
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  4. This notice is to remain intact.
 *                                                                    
 *  ------------------------------------------------------------------- 
 *
 *  The Algorithms are rederived in order to suit for Implicit Algorithms.
 *  For implicit integration: linear strain and displacement relation is 
 *  employed
 *                                                                    
 *  ------------------------------------------------------------------- 
 *                                                                     
 *  Written by: Xiaoguang Chen                                       January 1995
 *  ============================================================================= 
 */

#include <stdio.h>
#include <string.h>
#include <math.h>

#include "defs.h"
#include "units.h"
#include "matrix.h"
#include "vector.h"
#include "fe_database.h"
#include "symbol.h"
#include "fe_functions.h"
#include "elmt.h"
#include "miscellaneous.h"

/*
#define DEBUG 
#define DEBUG1
#define DEBUG3
*/


/* ============================================================== */
/*   Element SHELL_FOUR_NODES                                     */
/*        Implicit code                                           */
/*        Shell Element:                                          */
/*        material properties array                               */
/*        Input Properties:                                       */
/* ============================================================== */
/*    p->work_material[0] = E;
      p->work_material[1] = G;
      p->work_material[2] = fy;
      p->work_material[3] = ET;
      p->work_material[4] = nu;
      p->work_material[5] = density;
      p->work_material[6] = fu;

      p->work_section[0] = Ixx;
      p->work_section[1] = Iyy;
      p->work_section[2] = Izz;
      p->work_section[3] = Ixy;
      p->work_section[4] = Ixz;
      p->work_section[5] = Iyz;
      p->work_section[6] = weight;
      p->work_section[7] = bf;
      p->work_section[8] = tf;
      p->work_section[9] = depth;                                  
      p->work_section[10] = area;                                 
      p->work_section[11] = thickness;                            */
/* ============================================================== */

ARRAY *elmt_shell(p, isw)
ARRAY    *p;
int     isw;
{

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

   p = elmt_shell_4nodes_implicit(p, isw);

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

   return(p);
}

ARRAY *elmt_shell_4nodes_implicit(p, isw)
ARRAY    *p;
int     isw;
{
static double                                   nu; 
static QUANTITY      fy, G, E, ET, density, lambda; 
static double         Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
static double        bf, tf, A, depth, h, EA, EIzz;

static DIMENSIONS *dp_length, *dp_force, *dp_moment;
static DIMENSIONS *dp_stress, *dp_degree, *dp_temperature;

double                               x21, x31, x42;
double                               y21, y31, y42;
double                               z21, z31, z42;

double                      co_x31, co_y31, co_z31;
double                      co_x42, co_y42, co_z42;

double             EE, Jac,jacobian, k_bar = 0.833;     /* k_bar is shear factor */
double                   elmt_length, aspect_ratio;

double                                  **co_coord;     /*  co_coord is a coordinate rotated with material deformation */
double                   **T_matrix, **T_Transpose;
double                          **Direction_Matrix;

static double                *integ_coord, *weight;

double                           **stress, **displ;
double             **nodal_load, **nodal_load_temp;
double                             diff, sum, temp;
double                      **Cep, **stiff, **mass; 
int         i, j, k, ii, jj, kk, n, n1, n2, k1, k2;
int      surface_pts, surf_integ_pts, no_integ_pts;
int            size, dof, length, length1, length2;
int                                   UNITS_SWITCH;


#ifdef DEBUG
       printf("*** Enter elmt_shell_4nodes_implicit() : isw = %4d\n", isw);
#endif


/* =======================================================*/
/* Calculation of Direction Matrix; Transform coordinate  */
/* into lamina coordiante system                          */
/* =======================================================*/
     
       UNITS_SWITCH = CheckUnits();
       co_coord         = MatrixAllocIndirectDouble(p->no_dimen, p->nodes_per_elmt);
       Direction_Matrix = MatrixAllocIndirectDouble(6,6);

       for( i = 1; i <= 3; i++ ) {
           for ( j = 1; j <= p->nodes_per_elmt; j++) {
               co_coord[i-1][j-1] = p->coord[i-1][j-1].value; /* set initial co_coord */
           }
       }

       Lamina_Sys_Implicit(p, Direction_Matrix, co_coord);

#ifdef DEBUG
        dMatrixPrint("co_coord", co_coord, 3,4);
        dMatrixPrint("Direction_Matrix", Direction_Matrix, 6,6);
#endif

      x21 =  p->coord[0][1].value - p->coord[0][0].value;
      y21 =  p->coord[1][1].value - p->coord[1][0].value;
      z21 =  p->coord[2][1].value - p->coord[2][0].value;

      x31 =  p->coord[0][2].value - p->coord[0][0].value;
      y31 =  p->coord[1][2].value - p->coord[1][0].value;
      z31 =  p->coord[2][2].value - p->coord[2][0].value;

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

      co_x31 =  co_coord[0][2] - co_coord[0][0];
      co_y31 =  co_coord[1][2] - co_coord[1][0];
      co_z31 =  co_coord[2][2] - co_coord[2][0];

      co_x42 =  co_coord[0][3] - co_coord[0][1];
      co_y42 =  co_coord[1][3] - co_coord[1][1];
      co_z42 =  co_coord[2][3] - co_coord[2][1];

      Jac = (x31*y42-x42*y31);
      no_integ_pts = p->integ_ptr->thickness_pts;
      dof          = p->dof_per_node;
      size         = p->size_of_stiff;

/*****************************************************************/
#ifdef DEBUG
    printf(" ***** In elmt_shell_4nodes_implicit(): \n");
    printf("                                      : enter main switch \n");
#endif

    switch(isw) {
        case PROPTY: 

#ifdef DEBUG
    printf(" In elmt_shell_4nodes_implicit(): \n");
    printf("    : enter case of PROPTY\n");
#endif

             /* material properties:  elastic */

             E.value       =  p->work_material[0].value;
             fy.value      =  p->work_material[2].value;
             ET.value      =  p->work_material[3].value;
             nu            =  p->work_material[4].value;
             density.value =  p->work_material[5].value;
             lambda.value  =  nu*E.value/(1.0+nu)/(1.0-2.0*nu);

           /* (1)   check  possion_ratio value */

            if( nu == 0.0 || nu > 0.5 ) {
                printf("WARNING >> ... In shell_Belytschko element() -  nu value = %9.4f,reset to 0.3 !\n",
                       nu);
                nu = 0.3;    /* default poi_ratio value */
            }

            /* (2)   calculate  G value */
            
            G.value = p->work_material[1].value = E.value/(1.0 - 2.0*nu) ;

            if(E.value/((1.0 - 2.0*nu)) != p->work_material[1].value) {
                    printf(" elmt_shell(): WARNING: G is not equal to E/(1-2nu), check G for homogeneous material \n");
                    printf("             : ignore this message for non-homogeneous materials \n");
            }


            if( UNITS_SWITCH == ON ) {
               E.dimen       =  p->work_material[0].dimen;
               fy.dimen      =  p->work_material[2].dimen;
               ET.dimen      =  p->work_material[3].dimen;
               density.dimen =  p->work_material[5].dimen;
               lambda.dimen  =  E.dimen;
               G.dimen = E.dimen;
            }

             Ixx    = p->work_section[0].value;
             Iyy    = p->work_section[1].value;
             Izz    = p->work_section[2].value;
             Ixy    = p->work_section[3].value;
             Ixz    = p->work_section[4].value;
             Iyz    = p->work_section[5].value;
             bf     = p->work_section[7].value;
             tf     = p->work_section[8].value;
             depth  = p->work_section[9].value;
             A      = p->work_section[10].value;
             h      = p->work_section[11].value;    /* thickness of the shell */

             EA     = E.value*A;
             EIzz   = E.value*Izz;
         
#ifdef DEBUG
    printf(" In elmt_shell_4nodes_implicit(): \n");
    printf("    : leaving case of PROPTY\n");
#endif
             break;
             case STIFF: /* form element stiffness */

#ifdef DEBUG
       printf("*** In elmt_shell() : start case STIFF\n");
       printf("                    : Density         = %14.4e\n", density.value);
       printf("                    : shell_thickness = %8.2f\n", h);
       printf("                    : Jac             = %8.2f\n", Jac);
       printf("                    : E               = %lf\n", E.value);
       printf("                    : nu              = %lf\n", nu);
       printf("                    : no_integ_pts in z_direction = %d\n", no_integ_pts);
#endif

                 /* Integration loop over thickness */
      
                 stiff        = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
                 T_matrix     = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
                 T_Transpose  = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
                 integ_coord  = dVectorAlloc(no_integ_pts + 1);
                 weight       = dVectorAlloc(no_integ_pts + 1);

                 size         = p->size_of_stiff;
                 gauss(integ_coord,weight,no_integ_pts);

                 /* Integration loop over surface */

                 for (ii = 1; ii <= no_integ_pts; ii++) {

                      Shell_Stiff_Plane_4node(stiff, p, co_coord, integ_coord[ii], ii, E.value, nu); 	

                    for ( i = 1; i <= size; i++) {
                        for ( j = 1; j <= size; j++) {
                            p->stiff->uMatrix.daa[i-1][j-1] += stiff[i-1][j-1]*h*0.5*weight[ii];     
                        }
                    }
                  
                 }  /* end of gaussian integration through thickness */

#ifdef DEBUG
     printf(" End of integration through the thickness \n");
#endif

#ifdef DEBUG
           dMatrixPrint(" stiff before rotation", p->stiff->uMatrix.daa, p->size_of_stiff, p->size_of_stiff); 
#endif
                 free((char *) integ_coord);
                 free((char *) weight);

                 /***************************************************/
                 /* Transform stiffness matrix from local lamina    */
                 /* coordinate system to global coordinate system   */
                 /***************************************************/

                /* store the direction matrix : T_matrix */

                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;
                        }
                    }
                }

                stiff =  dMatrixMultRep(stiff, p->stiff->uMatrix.daa, size, size, T_matrix, size, size);

                for(i = 1; i <= size; i++){
                    for(j = 1; j <= size; j++) {
                        T_Transpose[i-1][j-1] = T_matrix[j-1][i-1];
                    }
                }
                p->stiff->uMatrix.daa = dMatrixMultRep(p->stiff->uMatrix.daa,
                                 T_Transpose, size, size, stiff, size, size);
    
#ifdef DEBUG
               /* check the symmetry of the stiffness matrix */

                for(i = 1; i <= size; i++){
                    for(j = 1; j <= size; j++) {
                        diff = p->stiff->uMatrix.daa[i-1][j-1] -p->stiff->uMatrix.daa[j-1][i-1];

                        if(diff > 1E-7 && 
                           (ABS(p->stiff->uMatrix.daa[i-1][j-1])) > 1E-7 &&
                           (ABS(diff/p->stiff->uMatrix.daa[i-1][j-1])) > 1E-1) {

                           printf("K[%d][%d] = %lf \n",
                                   i, j, p->stiff->uMatrix.daa[i-1][j-1]);
                           printf("K[%d][%d] = %lf \n",
                                   j, i, p->stiff->uMatrix.daa[j-1][i-1]);
                           printf("diff = %le \n", 
                                  p->stiff->uMatrix.daa[i-1][j-1] -
                                  p->stiff->uMatrix.daa[j-1][i-1]);
                           printf("elmtNo = %d Stiffness matrix IS NOT SYMMETRIC \n",
                                  p->elmt_no);
                           break;
                       }
                       else {
                           if( i == size && j == size)
                              printf("elmtNo = %d Stiffness matrix IS SYMMETRIC \n",
                                      p->elmt_no);
                       }