elmt_plate.c

有限元分析源代码

/*
 *  ============================================================================= 
 *  ALADDIN Version 1.0 :
 *         elmt_plate.c : Quadrilateral DKT Plate Bending 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.
 *                                                                    
 *  ------------------------------------------------------------------- 
 *  Element 8 - Quadrilateral DKT Plate Bending Element               
 *            - Written by B.K. Voon (1990)                           
 *                                                                    
 *  Note : 2-dimensional plate element from Berkeley.               
 *       : Only works for elements lying in x-z plane. Here we      
 *         assume that y-axis corresponds to vertical direction.   
 *  ------------------------------------------------------------------- 
 *                                                                     
 *  Written by: Mark Austin, Xiaoguang Chen, and Wane-Jang Lin      December 1995
 *  ============================================================================= 
 */

#include "defs.h"
#include "units.h"
#include "matrix.h"
#include "fe_database.h"
#include "symbol.h"
#include "fe_functions.h"

/*
#define DEBUG 
*/

/* ================================================== */
/* elmt_plate() : Quadrilateral DKT Plate Bending Element */
/* ================================================== */
/*    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_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;                            */
/* ---------------------------------------------------------------*/
/*    p->work_section[i] for i>11 is used to fit each elmts needs.*/
/*    Therefore, p->work_section[i] may be different for differnt */
/*    elmts when i > 11. For DKQ(T) plate element :               */ 
/* ---------------------------------------------------------------*/
/*    p->work_section[12] = D = (E/(1-nu^2))*t^3/12;              */
/*    p->work_section[13] = D*nu;                                 */
/*    p->work_section[14] = D*(1-nu)/2;                           */

/*    p->work_section[16] = flag   ?                              */
/* ============================================================== */
ARRAY *elmt_plate(p,isw)
ARRAY *p;
int isw;
{
static double nu;
static QUANTITY fy, G, E, ET, density;
static double  Ixx, Iyy, Izz, Ixy, Ixz, Iyz, bf, tf, A, depth, weight, EA, thickness;

double ia[3][27],eps[4],sig[4],aa[5],bb[5],cc[5],dd[5],ee[5];
double sg[10],tg[10],wg[10];
double d1,d2,d3,d4,d5,d6,thick,xsj,dn1,dn2,dn3;
double        xx,yy;
double        **bmu;
double        **shp;
double         **bm;
int i,j,l,lint,i1,k;
int    length1, length2, unit_length;
int                     UNITS_SWITCH;
int           ii, jj, kk, k1, k2, k3;
DIMENSIONS  *dimen, *dimen1, *dimen2;

   H_Print = 0;
   UNITS_SWITCH = CheckUnits();

#ifdef DEBUG
       printf("In elmt_plate() : starting with isw = %10d\n", isw);
#endif

   if(isw != PROPTY) {
      shp = MatrixAllocIndirectDouble(3,8);
      bm  = MatrixAllocIndirectDouble(3,12);  
   }

   switch(isw) {
       case PROPTY:   /* Input The Material Properties */
            if(p->work_material[16].value == 0.0 ) {
              /* d[10] has been changed to work_material[16].value    */
              /* what is d[10] any way ??? a flag ? array start at 1  */
              /* p->d[1];   Modulus    : p->d[2];   Poisson   */
              /* p->d[3];   Thickness  : p->d[4];   q-loading */

             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;
             thickness     =  p->work_section[11].value;
             p->work_section[12].value
                           = E.value/(1.0-nu*nu)*thickness*thickness*thickness/12.0; 
             p->work_section[13].value 
                           = nu*p->work_section[12].value; 
             p->work_section[14].value 
                           = 0.5*(1-nu)*p->work_section[12].value; 
             p->work_section[16].value = 1.0;

             switch(UNITS_SWITCH) {
                case 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;

                   dimen            = UnitsPower( p->work_section[11].dimen, 3.0, NO );
                   UnitsMultRep( p->work_section[12].dimen, E.dimen, dimen );
                   free((char *) dimen->units_name);
                   free((char *) dimen);
                   UnitsCopy( p->work_section[13].dimen, p->work_section[12].dimen );
                   UnitsCopy( p->work_section[14].dimen, p->work_section[12].dimen );
                break;
                case OFF:
                break;
                default:
                break;
             }
     

/* ---------------------------------------------------------*/
/* In fortran FEAP version                                  */
/* construct rotation parameters: u_x = 1, u_y = 2          */
/*                      ia[1][iel] = 1,  iel = elmt_type no */
/*                      ia[2][iel] = 2                      */
/* construct rotation parameters: theta_x = 4, theta_y = 5  */
/*                      ir[1][iel] = 4,  iel = elmt_type no */
/*                      ir[2][iel] = 5                      */
/* ---------------------------------------------------------*/
/*-----
               ia[0][p->eiel-1] = 2;         
               ia[1][p->eiel-1] = 3;        
-------*/

#ifdef DEBUG
               printf("In elmt_plate() : modulus   = %10.4f\n",E.value);
               printf("            : poisson   = %10.4f\n",nu);
               printf("            : thickness = %10.4f\n",thickness);
               printf("            : D         = %10.4f\n", p->work_section[12].value);
               printf("            : D*v       = %10.4f\n", p->work_section[13].value);
               printf("            : D(1-v)/2  = %10.4f\n", p->work_section[14].value);
               printf("            : p->work_section[16] = %10.4f\n\n", p->work_section[16].value);
#endif

            }
            break;
       case CHERROR:
       case MASS_MATRIX:
            break;
       case STIFF:   /* Compute The Element Tangent Array */

#ifdef DEBUG
       printf("In elmt_plate() : Start to compute Element Tangent Array\n");
#endif

            jacqud(p->coord,aa,bb,cc,dd,ee);
            l = 3;
            pgauss(l,&lint,sg,tg,wg);

            for(l=1; l <= lint; l++) {

#ifdef DEBUG
       printf("\n\n *** STARTING INT POINT : LINT %3d\n\n", l);
#endif

                shp = qushp8(sg[l-1],tg[l-1],shp,p->coord,&xsj);
                xsj = xsj*wg[l-1];
                bm  = dktqbm(shp,bm,aa,bb,cc,dd,ee);

                /* Compute Weighted Jacobian And D-matrix Constants */

                d1 = p->work_section[12].value * xsj;
                d2 = p->work_section[13].value * xsj;
                d3 = p->work_section[14].value * xsj;

                /* Compute The Element Load Vector */

/* nodal loads 
                for(i = 1; i <= p->nodes_per_elmt; i++) {
                    j = (i-1)*3+1; 
                    p->nodal_loads[j-1].value += (p->d[8])*xsj*(1.-sg[l-1])*(1.-tg[l-1]);
                }
*/
               /* Compute Contribution To Element Stiffness For This Point */

                for(i=1;i<=12;i++) {
                    dn1 = d1*bm[0][i-1]+d2*bm[1][i-1];
                    dn2 = d2*bm[0][i-1]+d1*bm[1][i-1];
                    dn3 = d3*bm[2][i-1];
                    for(j=i;j<=12;j++) 
                        p->stiff->uMatrix.daa[i-1][j-1] += dn1*bm[0][j-1]+dn2*bm[1][j-1]+dn3*bm[2][j-1];
                }
            }

            /* Make Stiffness Symmetric */

            for(i= 2; i<=12; i++) {
                i1 = i-1;
                for(j=1;j<=i1;j++)
                    p->stiff->uMatrix.daa[i-1][j-1] = p->stiff->uMatrix.daa[j-1][i-1];
            }

            /* ========================== feature not included in program yet ============
               Modify Load Vector for Non-Zero Displacement Boundary Conditions : p->u[12][1] 
            for(j=1;j<=p->size_of_stiff; j++) {
                if(p->displ->uMatrix.daa[j-1][0] != 0)
                   for(i=1;i<=p->size_of_stiff;i++)
                       p->nodal_loads[i-1].value -= p->stiff->uMatrix.daa[i-1][j-1]*p->displ->uMatrix.daa[j-1][0];

            } 
              ========================== end of feature =================================== */
           
            /**************************************************/
            /* Assign Units to Stiffness Matrix               */
            /**************************************************/

            /* Initiation of Stiffness Units Buffer           */

            if( CheckUnits() == ON ) {
               if(UNITS_TYPE == SI) {
                  dimen1 = DefaultUnits("Pa");
                  dimen2 = DefaultUnits("m");
               }
               else {
                  dimen1 = DefaultUnits("psi");
                  dimen2 = DefaultUnits("in");
               }

               /* node 1*/

               UnitsMultRep( &(p->stiff->spColUnits[0]), dimen1, dimen2 );
               UnitsMultRep( &(p->stiff->spColUnits[1]), &(p->stiff->spColUnits[0]), dimen2 );
               UnitsCopy( &(p->stiff->spColUnits[2]), &(p->stiff->spColUnits[1]) );

               ZeroUnits( &(p->stiff->spRowUnits[0]) ); 
               UnitsCopy( &(p->stiff->spRowUnits[1]), dimen2 ); 
               UnitsCopy( &(p->stiff->spRowUnits[2]), dimen2 ); 

               /* node i  i > 1*/

               for(i = 2; i <= p->nodes_per_elmt; i++) {