elmt_frame2d.c

有限元分析源代码

                    free((char *) dp_force->units_name);
                    free((char *) dp_force);
                } else {
                    printf("Coords (X,Y) = (%10.3f , %10.3f )\n", xx, yy);
                    printf("exx = %13.5e , curva = %13.5g , gamma = %13.5e\n", eps,chi, gam);
                    printf("\n");

                    /* node i */
                    printf(" Fx1 = %13.5e   Fy1 = %13.5e   Mz1 = %13.5e \n",
                        p->nodal_loads[0].value, 
                        p->nodal_loads[1].value, 
                        p->nodal_loads[2].value);
                    /* node j */
                    printf(" Fx2 = %13.5e   Fy2 = %13.5e   Mz2 = %13.5e \n",
                        p->nodal_loads[3].value,
                        p->nodal_loads[4].value,
                        p->nodal_loads[5].value);
                    printf("\n");
                    /* Member Forces */
                    printf(" Axial Force : x-direction = %13.5e \n", -p->nodal_loads[0].value);
                    printf(" Shear Force : y-direction = %13.5e \n",  p->nodal_loads[1].value);
                    printf("\n");
                }
             }
             break;
        case STIFF: /* form element stiffness */

#ifdef DEBUG
       printf("*** In elmt_frame_2d() : start case STIFF\n", isw);
#endif

             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;

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

             if(p->nodes_per_elmt == 2) /* elastic elment use 2-node elmt */
                p->stiff = beamst(p, p->stiff, EA, EIzz, xl, cs, sn,
                                  p->size_of_stiff, p->dof_per_node);

#ifdef DEBUG
      for(i = 1; i <= p->size_of_stiff; i++) {
         sum = 0.0;
         temp = 0.0;
         for(j = 1; j <= p->dof_per_node; j++) {
            for(k = 1; k <= p->nodes_per_elmt; k++) {
                ii = p->dof_per_node*(k-1)+j;
                if(j <= 2)
                  sum += p->stiff->uMatrix.daa[i-1][ii-1];
                if(j > 1)
                  temp += p->stiff->uMatrix.daa[i-1][ii-1];
            }
         }
         printf(" Force[%d] = %lf\n", i, sum);
         printf(" Moment[%d] = %lf\n", i, temp);
      }
#endif

             

#ifdef DEBUG
       MatrixPrintIndirectDouble(p->stiff);
       printf("*** In elmt_frame_2d() : end case STIFF\n", isw);
#endif

             break;
        case MASS_MATRIX:

#ifdef DEBUG
       printf("*** In elmt_frame_2d() : start case MASS\n", isw);
       printf("                : EA      = %8.2f\n", EA);
       printf("                : EI      = %8.2f\n", EIzz);
       printf("                : Density = %8.2f\n", p->work_material[5].value);
#endif


             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;

        /* Calculate mass = m_bar = mass/length                   */
        /* in units of (kg/m) or ((lbf-sec^2/in)/in)=(lb/in)      */
        /* if no units, then assume gravity g = 9.80665 m/sec^2   */

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

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

             p->stiff = beamms(p, p->stiff, p->type, mass, xl, cs, sn,
                               p->size_of_stiff, p->dof_per_node);

#ifdef DEBUG
       printf("mass\n");
       MatrixPrintIndirectDouble( p->stiff);
       printf("*** In elmt_frame_2d() : end case MASS\n", isw);
#endif
             break;
        default:
             break;
    }
#ifdef DEBUG
       printf("*** leaving elmt_frame_2d() \n");
#endif

    return(p);
}


/* ==================== */
/* Beam Stiffness       */
/* ==================== */

#ifdef __STDC__
MATRIX *beamst(ARRAY *p, MATRIX *s, double EA, double EI, double length, double cs, double sn, int size_of_stiff, int no_dof)
#else
MATRIX *beamst(p, s, EA, EI, length, cs, sn, size_of_stiff, no_dof)
ARRAY  *p;
MATRIX *s;
double  EA, EI, length, cs, sn;
int     size_of_stiff, no_dof ;
#endif
{
int             i, j, k;
int    length1, length2;
double                t;
DIMENSIONS     *d1, *d2;

/***********************************************/
/* Elastic Stiffness Matrix with two node at   */
/* each element                                */
/***********************************************/
    i = no_dof + 1;
    j = no_dof + 2; 
    k = no_dof + 3;

#ifdef DEBUG
 printf(" in beamst(): EA = %lf length = %lf \n", EA, length);
 printf(" in beamst(): EI = %lf cs= %lf sn = %lf \n", EI, cs, sn);
#endif


    t = EA/length; 
 
    s->uMatrix.daa[0][0]     =  t;
    s->uMatrix.daa[i-1][i-1] =  t;
    s->uMatrix.daa[0][i-1]   = -t;
    s->uMatrix.daa[i-1][0]   = -t;

    t = 12 * EI /(length*length*length);

    s->uMatrix.daa[1][1]     =  t;
    s->uMatrix.daa[j-1][j-1] =  t;
    s->uMatrix.daa[1][j-1]   = -t;
    s->uMatrix.daa[j-1][1]   = -t;

    t = (EI+ EI) / length;

    s->uMatrix.daa[2][2]   = s->uMatrix.daa[k-1][k-1] = t+t;
    s->uMatrix.daa[2][k-1] = s->uMatrix.daa[k-1][2]   = t;

    t = 6 * EI/length/length;
    s->uMatrix.daa[1][2]   = s->uMatrix.daa[2][1]       = t;
    s->uMatrix.daa[1][k-1] = s->uMatrix.daa[k-1][1]     = t;
    s->uMatrix.daa[2][j-1] = s->uMatrix.daa[j-1][2]     = -t;
    s->uMatrix.daa[j-1][k-1] = s->uMatrix.daa[k-1][j-1] = -t;

    /* ==================== */
    /* units buffer         */
    /* ==================== */

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

       UnitsMultRep( &(s->spColUnits[0]), d1, d2 );
       UnitsCopy( &(s->spColUnits[1]), &(s->spColUnits[0]) );
       UnitsMultRep( &(s->spColUnits[2]), &(s->spColUnits[0]), d2 );

       ZeroUnits( &(s->spRowUnits[0]) );
       ZeroUnits( &(s->spRowUnits[1]) );
       UnitsCopy( &(s->spRowUnits[2]), d2 ); 

       UnitsCopy( &(s->spColUnits[3]), &(s->spColUnits[0]) ); 
       UnitsCopy( &(s->spColUnits[4]), &(s->spColUnits[1]) ); 
       UnitsCopy( &(s->spColUnits[5]), &(s->spColUnits[2]) ); 

       UnitsCopy( &(s->spRowUnits[3]), &(s->spRowUnits[0]) ); 
       UnitsCopy( &(s->spRowUnits[4]), &(s->spRowUnits[1]) ); 
       UnitsCopy( &(s->spRowUnits[5]), &(s->spRowUnits[2]) ); 

       free((char *) d1->units_name);
       free((char *) d1);
       free((char *) d2->units_name);
       free((char *) d2);
    }

    s->uMatrix.daa = (double **) rotate(s->uMatrix.daa, cs, sn, size_of_stiff, no_dof);

#ifdef DEBUG
  printf("flag: in beamst() : STIFF after rotation \n");
  MatrixPrintIndirectDouble( s );
#endif
    return(s);

}

/* =================== */
/* Beam Mass Matrix    */
/* =================== */

#ifdef __STDC__
MATRIX *beamms(ARRAY *p, MATRIX *s, int mtype, double mass, double xl, double cs, double sn, int nst, int ndf)
#else
MATRIX *beamms(p, s, mtype, mass, xl, cs, sn, nst, ndf)
ARRAY  *p;
MATRIX *s;
double mass, xl, cs, sn;  /* mass, length, cosine, ans sine */
int nst, ndf,   mtype;
#endif
{
int    i, j, k, n, i1, j1;
int      length1, length2;
double  t, s1, s2, s3, dv;
double sg[5], ag[5], ba[3], rba[3], bb[5], rbb[5];
DIMENSIONS  *d1, *d2, *d3;

#ifdef DEBUG
       printf("*** Enter beamms() : matrix type = %10d  \n", mtype);
       printf("                   : t           = %10.3e\n", (mass * xl/2));
#endif

    t = mass * xl/2;
    switch(mtype) {
	case LUMPED:
             s->uMatrix.daa[0][0] = s->uMatrix.daa[1][1]  =  t;
             s->uMatrix.daa[2][2] =  t*xl*xl/12.0;   /* Using 16.0 in original version */

             s->uMatrix.daa[ndf][ndf] = t;
             s->uMatrix.daa[ndf+1][ndf+1] = t;
             s->uMatrix.daa[ndf+2][ndf+2] = t*xl*xl/12.0;   /* Using 16.0 in original version */

	     break;
	case CONSISTENT:

             gauss(sg, ag, 4);

             for(n =1; n<=4; n++) {

                 dv = 0.5 * xl * mass * ag[n];

                 s1 = (1 + sg[n])/2;
                 s2 = s1 * s1;
                 s3 = s1 * s2;

                 ba[1] = 1 - s1;
                 ba[2] = s1;
                 bb[1] = 1 - 3*s2 + s3 + s3;
                 bb[2] = xl* (s1- s2 - s2 + s3);
                 bb[3] = 3*s2 - s3 - s3;
                 bb[4] = xl  *  (-s2 + s3 );

                 rba[1] = ba[1] * dv;
                 rba[2] = ba[2] * dv;
                 for(i=1; i<= 4; i++)
                     rbb[i] = bb[i] * dv;

                 i = ndf;
                 s->uMatrix.daa[0][0] = s->uMatrix.daa[0][0] + ba[1] * rba[1];
                 s->uMatrix.daa[0][i] = s->uMatrix.daa[0][i] + ba[1] * rba[2];
                 s->uMatrix.daa[i][i] = s->uMatrix.daa[i][i] + ba[2] * rba[2];

                 for(i=1; i<= 4; i++){
                     i1 = i + 1;
                     if(i1 > 3 ) i1 = i1 + 1;
                     for(j=1; j<= 4; j++){
                         j1 = j  + 1;
                         if(j1 > 3) j1 = j1 +1;
                         s->uMatrix.daa[i1-1][j1-1] = s->uMatrix.daa[i1-1][j1-1] + bb[i]* rbb[j];
                     }
                 }
             }

             for(i = 0; i < 4; i++)
                 for(j = 0; j < i; j++)
                     s->uMatrix.daa[i][j] = s->uMatrix.daa[j][i];