elmt_frame3d.c

有限元分析源代码

/*
 *  ============================================================================= 
 *  ALADDIN Version 1.0 :
 *       elmt_frame3d.c : Three Dimensional Frame 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.
 *                                                                    
 *  ------------------------------------------------------------------- 
 *  Three Dimensional Frame Element                                    
 *                                                                    
 *  Convention for Nodal Forces                                        
 *             +ve M     -  anticlockwise(RHT Rule)                    
 *             +ve X,Y,Z -  along +ve axis                             
 *  Convention for Member End  Forces                                  
 *             +ve M     -  Sagging Moments                           
 *             +ve SF    -  LHS upwards                               
 *             +ve AF    -  Tension(LHS outwards)                     
 *  ------------------------------------------------------------------- 
 *                                                                     
 *  Written by: Mark Austin, Xiaoguang Chen, and Wane-Jang Lin      December 1995
 *  ============================================================================= 
 */

#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"
/*
#define DEBUG
*/


/* ============================================================== */
/*   Element FRAME_3D                                             */
/*   3D   Frame Element                                           */
/*        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] = plate_thickness;
      p->work_section[12] = J;
      p->work_section[13] = rT;
      p->work_section[14] = width;
      p->work_section[15] = tw;                                   */
/* ============================================================== */

/* macro for J constant */

#define J_Const(x,y) ((1 - 0.63 * x/y)*( x * x* x* y/3))

#ifdef __STDC__
ARRAY *elmt_frame_3d(ARRAY *p, int isw)
#else
ARRAY *elmt_frame_3d(p, isw)
ARRAY *p;
int   isw;
#endif
{
static double  nu;
static QUANTITY  fy, G, E, ET, density;
static double  Ixx, Iyy, Izz, Ixy, Ixz, Iyz, bf, tf, A, depth, weight, EA, EIzz, J, rT;

double d1;
double cs, sn,tn,xl,xx,yy,zz,vv,xn,xm,mbar;
double **rot, **trot, **fr, **dlocal;
int    i,j,k,l;

DIMENSIONS  *dp_length, *dp_force, *dp_moment;
DIMENSIONS  *dp_stress, *dp_degree, *dp_temperature;
int          UNITS_SWITCH;

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

   H_Print = 0;
   UNITS_SWITCH = CheckUnits();

   switch(isw) {
       case PROPTY:  /* MAT PROPS */

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

          /* (1)   check  poi_ratio value */

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

          /* (2)   calculate  G value */
            
/*        if(E.value/((1.0 - 2.0*nu)) != p->work_material[1].value) {
              printf(" elmt_frame_3d(): WARNING: G is not equal to E/(1-2nu), check G for homogeneous material \n");
              printf("                : ignore this message for non-homogeneous materials \n");
          }
*/
          G.value = p->work_material[1].value = E.value/(1.0 - 2.0*nu) ;
          if(UNITS_SWITCH == ON)  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;
          weight = p->work_section[6].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;
          J      = p->work_section[12].value;
          rT     = p->work_section[13].value;

          EA     = E.value*A;
          EIzz   = E.value*Izz;

          /* (3) If J value not input, J calculated based on rectangular */
          /*     section of size (bf x depth)                            */

          if(J == 0.0 ) {
             if(bf == 0.0 || depth == 0.0){
                printf("WARNING >> Must give 'J' or ('width' & 'depth') to calculate stiffness");
                exit(1);
             }

             /* Check bf < depth & cal J */

             if(bf < depth ) 
                 J = p->work_section[12].value = J_Const(bf,depth); 
             else
                 J = p->work_section[12].value = J_Const(depth,bf); 
          }
          break;
       case CHERROR:
            break;
       case STRESS_LOAD:
            break;
       case STRESS_UPDATE:
            break;
       case PRESSLD:
            break;

       case STIFF:

            cs = p->coord[0][1].value - p->coord[0][0].value;           /* Cos Term */
            sn = p->coord[1][1].value - p->coord[1][0].value;           /* Sin Term */
            tn = p->coord[2][1].value - p->coord[2][0].value;           /* Tan Term */
            xl = sqrt(cs * cs + sn * sn + tn * tn); /* Calculate  Length */
            p->length.value = xl;

            /* T matrix is made here: 12*12 size */

            rot = (double **) MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
            rot = (double **) tmat(rot, 6, p); 

            if( UNITS_SWITCH==ON )   p->stiff->spColUnits[0].units_type = UNITS_TYPE;

            p->stiff = beamst3d(p, p->stiff, EA, EIzz, E.value*Iyy, G.value*J ,xl, rot, p->size_of_stiff, p->dof_per_node); 

            MatrixFreeIndirectDouble(rot, p->size_of_stiff);
            break;

       case MASS_MATRIX:

            cs = p->coord[0][1].value - p->coord[0][0].value;           /* Cos Term */
            sn = p->coord[1][1].value - p->coord[1][0].value;           /* Sin Term */
            tn = p->coord[2][1].value - p->coord[2][0].value;           /* Tan Term */
            xl = sqrt(cs * cs + sn * sn + tn * tn); /* Calculate  Length */
            p->length.value = xl;

            /* T matrix is made here: 12*12 size */

            rot = (double **) MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
            rot = (double **) tmat(rot, 6, p); 

            /* Assemble Mass Matrix */
 
            /* Calculate mbar =  mass/length                 */
            /* in units of (kg/m) or (lbf*sec^2/in/in)       */
            /* if no units, assume gravity g=9.80665 m/sec^2 */

             if( weight != 0.0 )
               mbar = weight/9.80665;
             else
               if( density.value > 0 )  mbar = A * density.value ;
             else {
                printf("\nError in input: Need density value to calculate mass matrix\n");
                exit(1);
             }

	    /* Calculate radius of gyration , rT  --  m  ,  in   */
            /* original version      :  rT = p->length.value/ 1.414; */

            if( rT==0 && A!=0 )    rT = sqrt( J / A );

            if( UNITS_SWITCH == ON ) 
                p->stiff->spColUnits[0].units_type = UNITS_TYPE;

            p->stiff = beamms3d(p, p->stiff,p->type, mbar, xl, rT, rot, p->size_of_stiff, p->dof_per_node);

            MatrixFreeIndirectDouble(rot, p->size_of_stiff);
            break;

       case STRESS:
       case LOAD_MATRIX:

            cs = p->coord[0][1].value - p->coord[0][0].value;           /* Cos Term */
            sn = p->coord[1][1].value - p->coord[1][0].value;           /* Sin Term */
            tn = p->coord[2][1].value - p->coord[2][0].value;           /* Tan Term */
            xl = sqrt(cs * cs + sn * sn + tn * tn); /* Calculate  Length */
            p->length.value = xl;

            /* T matrix is made here: 12*12 size */

            rot = (double **) MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
            rot = (double **) tmat(rot, 6, p); 

            /* --------------------------- */
            /* Output Stresses and Strains */
            /* --------------------------- */

            cs = rot[0][0];
            rot[0][0] = 1.0;

            if (UNITS_SWITCH==ON)  p->stiff->spColUnits[0].units_type = UNITS_TYPE;
            p->stiff = beamst3d(p, p->stiff, EA, EIzz, E.value*Iyy,
                       G.value*J, xl, rot, p->size_of_stiff, p->dof_per_node); 

            rot[0][0] = cs;

            fr = MatrixAllocIndirectDouble(p->size_of_stiff,1); /* Here size_of_stiff = 12; i.e 12x1 mat */

            for(l = 1; l<= p->nodes_per_elmt; l++) {
                for(k = 1; k <= p->dof_per_node; k++) {
                    j = (l-1)*p->dof_per_node + k;
                    fr[j-1][0] = p->displ->uMatrix.daa[k-1][l-1];
                }
            }

            dlocal = (double **)  dMatrixMult(rot, p->size_of_stiff, p->size_of_stiff, fr, p->size_of_stiff, 1);

            fr = (double **)  dMatrixMultRep(fr, p->stiff->uMatrix.daa, p->size_of_stiff, p->size_of_stiff, dlocal, p->size_of_stiff, 1);

            xx = 0.5 *(p->coord[0][0].value +p->coord[0][1].value);
            yy = 0.5 *(p->coord[1][0].value +p->coord[1][1].value);
            zz = 0.5 *(p->coord[2][0].value +p->coord[2][1].value);

            if(H_Print == YES && isw == STRESS){
               printf( "\n");
               printf( "3D Frame Element   Stresses\n");
               printf( "---------------------------\n");
               H_Print = NO;
            }

            if(isw == STRESS && PRINT_STRESS == ON ) {
               printf("\n");
               printf("Elmt# %3d : ", p->elmt_no);
               if(UNITS_SWITCH == ON)
                  printf(": Coords (X,Y,Z)= (%8.3f %s,%8.3f %s,%8.3f %s)\n",
                            xx/p->coord[0][0].dimen->scale_factor, p->coord[0][0].dimen->units_name,
                            yy/p->coord[1][0].dimen->scale_factor, p->coord[1][0].dimen->units_name,
                            zz/p->coord[2][0].dimen->scale_factor, p->coord[2][0].dimen->units_name);
	       else
                  printf(": Coords (X,Y,Z)= (%8.3f,%8.3f,%8.3f)\n", xx, yy, zz);

               printf("\n");
	    }

            /*--------------------------------------------------------*/
            /* nodal forces   & member end forces                     */
            /*--------------------------------------------------------*/

            if(p->elmt_load_ptr != NULL ) { /* calculate FEF */
               printf("Fixed End Loads; \n");
               p = sld05(p, STRESS);

               /* Add FEF to joint p->[  ]  orces */
                for(i = 1; i <= 12; i++)
                  fr[i-1][0] = fr[i-1][0] - p->nodal_loads[i-1].value;
            }

            /* -------------------- */
            /* Print Element forces */
            /* -------------------- */

            for(j = 1; j <= p->size_of_stiff; j++) 
                p->nodal_loads[j-1].value = fr[j-1][0];

            /* Assign element forces's units */

            if( UNITS_SWITCH == ON ) {
                if(UNITS_TYPE == SI) {
                  dp_length = DefaultUnits("m");
                  dp_force  = DefaultUnits("N");
                }
                if(UNITS_TYPE == US) {
                  dp_length = DefaultUnits("in");
                  dp_force  = DefaultUnits("lbf");
                }
                dp_moment = UnitsMult( dp_force, dp_length );

                for(j=1;j<=3;j++) {
		    UnitsCopy( p->nodal_loads[j-1].dimen, dp_force );
		    UnitsCopy( p->nodal_loads[j-1+p->dof_per_node].dimen, dp_force );
		    UnitsCopy( p->nodal_loads[j-1+3].dimen, dp_moment );
		    UnitsCopy( p->nodal_loads[j-1+3+p->dof_per_node].dimen, dp_moment );
		}
            }

            if(isw == STRESS && PRINT_STRESS == ON ) {
	       switch(UNITS_SWITCH) {
		 case ON:
                  /* node_i */
                  printf("            Fx1 = %13.5e %s\t Fy1 = %13.5e %s\t Fz1 = %13.5e %s\n",
                         p->nodal_loads[0].value/dp_force->scale_factor, dp_force->units_name,
                         p->nodal_loads[1].value/dp_force->scale_factor, dp_force->units_name,
                         p->nodal_loads[2].value/dp_force->scale_factor, dp_force->units_name);
                  printf("            Mx1 = %13.5e %s\t My1 = %13.5e %s\t Mz1 = %13.5e %s\n",
                         p->nodal_loads[3].value/dp_moment->scale_factor, dp_moment->units_name,
                         p->nodal_loads[4].value/dp_moment->scale_factor, dp_moment->units_name,
                         p->nodal_loads[5].value/dp_moment->scale_factor, dp_moment->units_name);
                  printf("\n");
                  /* node_j */
                  printf("            Fx2 = %13.5e %s\t Fy2 = %13.5e %s\t Fz2 = %13.5e %s\n",
                        p->nodal_loads[6].value/dp_force->scale_factor, dp_force->units_name,
                        p->nodal_loads[7].value/dp_force->scale_factor, dp_force->units_name,
                        p->nodal_loads[8].value/dp_force->scale_factor, dp_force->units_name);
                  printf("            Mx2 = %13.5e %s\t My2 = %13.5e %s\t Mz2 = %13.5e %s\n",
                        p->nodal_loads[9].value/dp_moment->scale_factor, dp_moment->units_name,
                        p->nodal_loads[10].value/dp_moment->scale_factor, dp_moment->units_name,
                        p->nodal_loads[11].value/dp_moment->scale_factor, dp_moment->units_name);
                  printf("\n");

                  printf("            Axial Force : x-direction = %13.5e %s \n",
                            -p->nodal_loads[0].value/dp_force->scale_factor,
                             dp_force->units_name);
                  printf("            Shear Force : y-direction = %13.5e %s \n",
                             p->nodal_loads[1].value/dp_force->scale_factor,
                             dp_force->units_name);
                  printf("                        : z-direction = %13.5e %s \n",
                             p->nodal_loads[2].value/dp_force->scale_factor,
                             dp_force->units_name);
                  printf("\n"); 
                 break;
                 case OFF:
                  /* node_i */
                  printf("            Fx1 = %13.5e\t Fy1 = %13.5e\t Fz1 = %13.5e\n",
                           p->nodal_loads[0].value,
                           p->nodal_loads[1].value,
                           p->nodal_loads[2].value); 
                  printf("            Mx1 = %13.5e\t My1 = %13.5e\t Mz1 = %13.5e\n",
                           p->nodal_loads[3].value,
                           p->nodal_loads[4].value,
                           p->nodal_loads[5].value); 
                  printf("\n"); 
                  /* node_j */
                  printf("            Fx2 = %13.5e\t Fy2 = %13.5e\t Fz2 = %13.5e\n",
                           p->nodal_loads[6].value,
                           p->nodal_loads[7].value,