elmt_frame2d.c

有限元分析源代码

/*
 *  ============================================================================= 
 *  ALADDIN Version 1.0 :
 *       elmt_frame2d.c : Two-dimensional Beam-Column 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.
 *                                                                    
 *  -------------------------------------------------------------------
 *  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 "fe_database.h"
#include "symbol.h"
#include "fe_functions.h"
#include "elmt.h"

/*
#define DEBUG 
*/

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

#ifdef __STDC__
ARRAY *elmt_frame_2d(ARRAY *p, int isw)
#else
ARRAY *elmt_frame_2d(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;

double cs, sn, xl,xx,yy,  eps,chi,gam, mzc,mass;
double vl1,vl2,tl1,tl2, fx1,fx2,fy1,fy2,mz1,mz2,e6,e2;
double  eplas , alp, sig;
int NS_Sub_Incr;
DIMENSIONS *dp_length, *dp_force, *dp_moment;
DIMENSIONS *dp_stress, *dp_degree, *dp_temperature;

double sum, temp;
int    i, j, k, ii;
int    UNITS_SWITCH;

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

    UNITS_SWITCH = CheckUnits();

    switch(isw) {
        case PROPTY:
             /* beam element :   material properties  */

             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 2d 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 */
            
            G.value = p->work_material[1].value = E.value/(1.0 - 2.0*nu) ;
            if( UNITS_SWITCH == ON )  G.dimen = E.dimen;

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

             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;

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

             break;
        case CHERROR:
             break;
        case STRESS_LOAD:
             break;
        case STRESS_UPDATE:
             break;
	case PRESSLD:
             break;
	case LOAD_MATRIX:
        case STRESS:

             cs = p->coord[0][1].value - p->coord[0][0].value;
             sn = p->coord[1][1].value - p->coord[1][0].value;

             xl = sqrt(cs * cs + sn * sn);
             cs = cs/xl; 
             sn = sn/xl;
             p->length.value = xl;

             xx  = 0.5*(p->coord[0][0].value + p->coord[0][1].value);            /* xx = 0.5(x1+x2) */
             yy  = 0.5*(p->coord[1][0].value + p->coord[1][1].value);            /* yy = 0.5(y1+y2) */
             eps = (cs*(p->displ->uMatrix.daa[0][1] - p->displ->uMatrix.daa[0][0])
                    +sn*(p->displ->uMatrix.daa[1][1] - p->displ->uMatrix.daa[1][0])) / xl;

                    /* eps = (u2-u1)cos(phi)+(v2-v1)sin(phi) */

             chi = (p->displ->uMatrix.daa[2][1] - p->displ->uMatrix.daa[2][0]) / xl;
             gam = -12 *(-sn*(p->displ->uMatrix.daa[0][0] - p->displ->uMatrix.daa[0][1])
                   + cs*(p->displ->uMatrix.daa[1][0]- p->displ->uMatrix.daa[1][1]))
                   / (xl*xl*xl)-6.*(p->displ->uMatrix.daa[2][0] + p->displ->uMatrix.daa[2][1])/xl/xl;

             /* local deflections */

             vl1 = -sn * p->displ->uMatrix.daa[0][0] + cs * p->displ->uMatrix.daa[1][0];
             vl2 = -sn * p->displ->uMatrix.daa[0][1] + cs * p->displ->uMatrix.daa[1][1];
             tl1 = p->displ->uMatrix.daa[2][0];
             tl2 = p->displ->uMatrix.daa[2][1];

             /* computer axial forces fx, transverse forces fy, and moment mz */

             e6 = 6* EIzz /xl/xl;
             e2 = 2* EIzz /xl;

             /* Elasto_Plastic load case                            */
             /* Considering uniaxial element with x_displ, Fx force */

             fx1 = - EA * eps;
             fy1 = 2*e6/xl*(vl1-vl2) + e6*(tl1 + tl2);
             fx2 = - fx1;
             fy2 = - fy1;
             mz1 = e6 * ( vl1 - vl2) + e2 * ( 2 * tl1 + tl2 );
             mz2 = e6 * ( vl1 - vl2) + e2 * (     tl1 + 2 * tl2 );
             mzc = (mz1 - mz2)/2;

             /* Add FEF if elmt loaded */
 
             if( p->elmt_load_ptr != NULL) { 
                 p = sld07(p, STRESS);

                 /* Add FEF to joint forces */

                 fx1 = fx1  - p->nodal_loads[0].value;
                 fy1 = fy1  - p->nodal_loads[1].value;
                 mz1 = mz1  - p->nodal_loads[2].value;
                 fx2 = fx2  - p->nodal_loads[3].value;
                 fy2 = fy2  - p->nodal_loads[4].value;
                 mz2 = mz2  - p->nodal_loads[5].value;
             } 

             /* Assign force values */

             if( UNITS_SWITCH == ON ) {
                if(UNITS_TYPE == SI) {
                    dp_length = DefaultUnits("m");
                    dp_force  = DefaultUnits("N");
                }
                else if(UNITS_TYPE == US) {
                    dp_length = DefaultUnits("in");
                    dp_force  = DefaultUnits("lbf");
                }
            
                /* node no 1 */
                UnitsCopy( p->nodal_loads[0].dimen, dp_force );
                UnitsCopy( p->nodal_loads[1].dimen, dp_force );
                UnitsMultRep( p->nodal_loads[2].dimen, dp_force, dp_length );

                /* node no = 2 */
                UnitsCopy( p->nodal_loads[3].dimen, p->nodal_loads[0].dimen );
                UnitsCopy( p->nodal_loads[4].dimen, p->nodal_loads[1].dimen );
                UnitsCopy( p->nodal_loads[5].dimen, p->nodal_loads[2].dimen );
             }

             p->nodal_loads[0].value = fx1; 
             p->nodal_loads[1].value = fy1; 
             p->nodal_loads[2].value = mz1;
             p->nodal_loads[3].value = fx2; 
             p->nodal_loads[4].value = fy2; 
             p->nodal_loads[5].value = mz2; 

             if(isw == LOAD_MATRIX ) {
                p->nodal_loads[0].value = fx1*cs - fy1*sn; 
                p->nodal_loads[1].value = fx1*sn + fy1*cs; 
                p->nodal_loads[2].value = mz1;
                p->nodal_loads[3].value = fx2*cs - fy2*sn; 
                p->nodal_loads[4].value = fx2*sn + fy2*cs; 
                p->nodal_loads[5].value = mz2; 
             }

             if(isw == STRESS && PRINT_STRESS == ON) {
                printf("\n");
                printf("Elmt No %3d : \n", p->elmt_no);
                if( UNITS_SWITCH == ON ) {
                    printf("Coords (X,Y) = (%10.3f %s, %10.3f %s)\n", 
                        xx/dp_length->scale_factor,dp_length->units_name,
                        yy/dp_length->scale_factor,dp_length->units_name);

                    printf("exx = %13.5e , curva = %13.5g , gamma = %13.5e\n", eps,chi, gam);
                    printf("\n");

                    /* node i */
                    printf(" Fx1 = %13.5e %s  Fy1 = %13.5e %s  Mz1 = %13.5e %s\n",
                        p->nodal_loads[0].value/p->nodal_loads[0].dimen->scale_factor,
                        p->nodal_loads[0].dimen->units_name,
                        p->nodal_loads[1].value/p->nodal_loads[1].dimen->scale_factor,
                        p->nodal_loads[1].dimen->units_name,
                        p->nodal_loads[2].value/p->nodal_loads[2].dimen->scale_factor,
                        p->nodal_loads[2].dimen->units_name);

                    /* node j */
                    printf(" Fx2 = %13.5e %s  Fy2 = %13.5e %s  Mz2 = %13.5e %s\n",
                        p->nodal_loads[3].value/p->nodal_loads[3].dimen->scale_factor,
                        p->nodal_loads[3].dimen->units_name,
                        p->nodal_loads[4].value/p->nodal_loads[4].dimen->scale_factor,
                        p->nodal_loads[4].dimen->units_name,
                        p->nodal_loads[5].value/p->nodal_loads[5].dimen->scale_factor,
                        p->nodal_loads[5].dimen->units_name);
                        printf("\n");
                    /* Member Forces */
                    printf(" Axial Force : x-direction = %13.5e %s\n",
                        -p->nodal_loads[0].value/p->nodal_loads[0].dimen->scale_factor,
                         p->nodal_loads[0].dimen->units_name);
                    printf(" Shear Force : y-direction = %13.5e %s\n",
                         p->nodal_loads[1].value/p->nodal_loads[1].dimen->scale_factor,
                         p->nodal_loads[1].dimen->units_name);
                    printf("\n");

                    free((char *) dp_length->units_name);
                    free((char *) dp_length);