elmt_frame3d.c

有限元分析源代码

    cx = p->coord[0][1].value - p->coord[0][0].value;
    cy = p->coord[1][1].value - p->coord[1][0].value;

    switch(type) {
        case 1: /* plane truss */
	case 2:
	case 3:
	case 4: /* plane frame, with axial force */
             el = sqrt(cx * cx + cy * cy );
             cx = cx/el; 
             cy = cy/el;
             if(type == 4) {
                t[0][0] = cx;
                t[0][1] = cy;
                t[1][0] = -cy;
                t[1][1] = cx;
                if(type == 1) break;
                   t[2][2] = 1;
             }
             else {
                t[0][0] = cx;
                t[0][2] = cy;
                t[1][1] = 1;
                t[2][0] = -cy;
                t[2][2] = cx;
             } 
             break;
        case 5: /* space truss */
        case 6: /* space frame */
             cz = p->coord[2][1].value - p->coord[2][0].value;
             el = sqrt(cx * cx + cy * cy + cz * cz);
             cx = cx/el; 
             cy = cy/el;
             cz = cz/el;
             if(type == 5){
                t[0][0] = cx*cx;
                t[0][1] = cx*cy;
                t[0][2] = cx*cz;
                t[1][1] = cy*cy;
                t[1][2] = cy*cz;
                t[2][2] = cz*cz;
                for(i = 2; i <= 3; i++){
                    ii = i -1;
                    for(j = 1; j <= ii; j++) 
                        t[i-1][j-1] = t[j-1][i-1];
                }
             }
             else {
                csa = cos(alpha);
                sna = sin(alpha);
                t[0][1] = cy;
                if(cx == 0.0 && cz == 0.0){  /* vertical space frame */
                   t[1][0] = -cy*csa;
                   t[1][2] = sna;
                   t[2][0] = cy * sna;
                   t[2][2] = csa;
                }
		else {
                   elone = sqrt(cx * cx + cz * cz);
                   t[0][0] = cx;
                   t[0][2] = cz;
                   t[1][0] = (-cx*cy*csa - cz*sna)/elone;
                   t[1][2] = (-cy*cz*csa + cx*sna)/elone;
                   t[1][1] = csa * elone;
                   t[2][0] = ( cx*cy*sna - cz*csa)/elone;
                   t[2][1] = -sna * elone;
                   t[2][2] = ( cy*cz*sna + cx*csa)/elone;
               }
            }
            break;
        default:
            break;
    }

    if(type == 6){
       nf = 4;
       ndff = 3;
    }
    else{
       nf   = 2;
       ndff = p->dof_per_node;
    }

    for(i=1;i<=nf;i++) {
        jj = (i-1) * ndff; 
        kk = (i-1) * ndff;
        for(j=1;j<=ndff;j++){
            jj = jj+ 1;
            for(n=1;n<=ndff;n++){
                kk = kk  + 1;
                rot[jj-1][kk-1] = t[j-1][n-1];
            }
            kk = (i-1) * ndff;
        }
    }

    return(rot);
}


/* ===================================================================== */
/* MATRIX *rotate3d(s, r ,ndf)                                           */
/* Rotation matrix  for 3d frame element                                 */
/* Multiplies  the stiffness matrix "s" by the transformation matrix "r" */
/* (ie. Kg = Tt*Kl*T)  Return Kg                                         */
/* ===================================================================== */

#ifdef __STDC__
double **rotate3d(double **s, double **r , int no_dof)
#else
double **rotate3d(s, r ,no_dof)
double **s, **r;
int    no_dof;
#endif
{
int    i,j,n;
double t;
double **rt, **asj, **sj, **sij, **kj;

     if(r[0][0] == 1.0) return(s);

  /* transpose of rotation matrix  r */

     rt = MatrixAllocIndirectDouble(no_dof, no_dof);
     for(i = 1; i <= no_dof; i++)
       for(n =1; n<= no_dof; n++)
           rt[n-1][i-1] = r[i-1][n-1];

  /* sii terms of transformed stiffness matrix s */

     kj = MatrixAllocIndirectDouble(no_dof, no_dof);
     for(i = 1; i <= no_dof; i++)
         for(n =1; n<= no_dof; n++)
             kj[i-1][n-1] = s[i-1][n-1];

     asj = (double **) dMatrixMult( kj, no_dof, no_dof, r, no_dof, no_dof);
     sj  = (double **) dMatrixMult( rt, no_dof, no_dof, asj, no_dof, no_dof);

     for(i=1;i<= no_dof; i++)
         for(n=1;n<=no_dof; n++)
             s[i-1][n-1] = sj[i-1][n-1];

     MatrixFreeIndirectDouble(asj, no_dof);
     MatrixFreeIndirectDouble( sj, no_dof);

  /* sjj terms of transformed stiffness matrix s */

     for(i = 1; i <= no_dof; i++)
         for(n =1; n<=no_dof;n++)
             kj[i-1][n-1] = s[i+no_dof-1][n+no_dof-1];

     asj = (double **) dMatrixMult( kj, no_dof, no_dof, r, no_dof, no_dof);
     sj  = (double **) dMatrixMult( rt, no_dof, no_dof, asj, no_dof, no_dof);

     for(i=1;i<= no_dof; i++)
         for(n=1;n<=no_dof; n++)
             s[i+no_dof-1][n+no_dof-1] = sj[i-1][n-1];

     MatrixFreeIndirectDouble(asj, no_dof);
     MatrixFreeIndirectDouble( sj, no_dof);

  /* sij terms of transformed stiffness matrix s */

     for(i = 1; i <= no_dof; i++)
         for(n =1; n<= no_dof; n++)
             kj[i-1][n-1] = s[i-1][n+no_dof-1];

     asj = (double **) dMatrixMult( kj, no_dof, no_dof, r, no_dof, no_dof);
     sj  = (double **) dMatrixMult( rt, no_dof, no_dof, asj, no_dof, no_dof);

     for(i=1;i<= no_dof; i++)
         for(n=1;n<= no_dof ;n++)
             s[i-1][n+no_dof-1] = s[n+no_dof-1][i-1] = sj[i-1][n-1];

     MatrixFreeIndirectDouble(asj, no_dof);
     MatrixFreeIndirectDouble( sj, no_dof);
     MatrixFreeIndirectDouble( kj, no_dof);
     MatrixFreeIndirectDouble( rt, no_dof);

     return(s);
}


/* ================================================== */
/* Equivalent Loading Procedure for 3 D frame element */
/* ================================================== */

#ifdef __STDC__
ARRAY *sld05(ARRAY *p, int task)
#else
ARRAY *sld05(p,task)
ARRAY *p;
int   task;
#endif
{
ELEMENT_LOADS   *elsptr;
ELOAD_LIB          *elp;
double  P, a ,b ;
double  L, load[15];
double  px,py,pz, mx,my,mz, bx,by,bz, ze,ze2,ze3;
double  X1,Y1,Z1,MX1,MY1,MZ1,
        X2,Y2,Z2,MX2,MY2,MZ2, 
        MCZ, MCY, MCZT, MCYT;
double  f1,f2,f3,f4,f5,f6,
        df1x, df3x, df5x,df2x, df6x;
double  **rmat, **rt;
int     *da;
int     inc,i,j;

/* printf(">>>>In sld05 : 3d elmt    \n");  */

       /* Initialize total load */

       for(inc=1; inc <= 12; inc++)
          load[inc-1] = 0.0;
       MCZT = 0.0;
       MCYT = 0.0;

    switch(task){
        case PRESSLD: case STRESS:
             L = p->length.value;  elsptr  =  p->elmt_load_ptr;

             for(j=1; j<= elsptr->no_loads_faces; j++) {
                 elp  =  &elsptr->elib_ptr[j-1];

                 P =  elp->P.value;
                 a =  elp->a.value;
                 b =  elp->b.value;

                 /* Pt loads */

                 px =  elp->px.value;
                 py =  elp->py.value;
                 pz =  elp->pz.value;

                 /* moments */

                 mx =  elp->mx.value;
                 my =  elp->my.value;
                 mz =  elp->mz.value;

                 /* distributed loading */

                 bx =  elp->bx.value;
                 by =  elp->by.value;
                 bz =  elp->bz.value;

                 if(a > L) /* error message */
                    printf(">>ERROR in sld; Elmt Load dist. 'a' > Elmt Length; El_no= %d\n",p->elmt_no);  

                 if (elp->type == -1) { /* Distributed loading Condition */
                     inc = 0;

                 /* set default values */
                 if(b == 0.0) b = L; /* dist loading acts on entire length */

                 /* first calc f(b) */
                    ze = b/L;    
SHP_START:
                    ze2 = ze * ze; ze3 = ze2 * ze;
                    f1 =    1   -  ze2/2;
                    f2 =    ze3*ze/2  - ze3 + ze ;
                    f3 =    (ze3*ze/4 - 2*ze3/3 + ze2/2) * L;
                    f4 =    ze2/2 ;
                    f5 =    -ze *ze3/2 +  ze3;
                    f6 =    (ze3*ze/4  - ze3/3) * L;
                    inc++;
                 
                    if(inc == 1) {
                       /* temp hold of values f(b) */
                       X1 = f1; Y1 = f2; Z1 = f3; 
                       X2 = f4; Y2 = f5; Z2 = f6; 
                       ze = a/L;
                       goto SHP_START;
                    }
                    else{
                       /* f() = f(b) - f(a)  */
                       f1 = X1 - f1; f2 =  Y1 - f2; f3 = Z1 - f3; 
                       f4 = X2 - f4; f5 =  Y2 - f5; f6 = Z2 - f6; 
                    }

                    X1 = bx * f1 * L;
                    Y1 = by * f2  * L;
                    Z1 = bz * f2 * L;
                    MX1 = 0.0 ;
                    MY1 =-bz * f3 * L;
                    MZ1 = by * f3 * L;

                    X2 = bx * f4 * L;
                    Y2 = by * f5 * L;
                    Z2 = bz * f5 * L;
                    MX2 = 0.0 ;
                    MY2 =-bz * f6 * L;
                    MZ2 = by * f6 * L;

                    if (task == STRESS){
                        /* +ve  simply support moment at center */
                        if(b==L && a== 0.0) {/* udl acting on entire length */
                           MCZ = -by * (L * L)/8;   
                           MCY = -bz * (L * L)/8;   
                        }
                    else { /* approximate mom at center */
                        MCZ = -by *(b-a)* (L - (a+b)/2)/2;   
                        MCY = -bz *(b-a)* (L - (a+b)/2)/2;   
                    }
                 }
               }      /* end of dist loading */

               else {                 /* Concentrated Loading Condition */
                   /* shape functions */
                  ze = a/L;      ze2 = ze * ze;        ze3 = ze2 * ze;

                  f1 =     1  -  ze;
                  f2 =     2 * ze3  -  3 * ze2 + 1;
                  f3 =    (ze3 - 2 * ze2 + ze) * L;
                  f4 =     ze ;
                  f5 =    -2 * ze3  +  3 * ze2;
                  f6 =    (ze3 - ze2 ) * L;

                  if(my != 0.0 || mz != 0.0){
                     /* derivatives of shape function */  
                     df2x =  6 *( ze2  -  ze) / L;
                     df3x =  3 *ze2  - 4*ze  + 1;
                     df5x = - df2x; 
                     df6x =  3 *ze2  - 2*ze;
                  }

                  X1 = px * f1 + 0;
                  Y1 = py * f2 + mz *  df2x;
                  Z1 = pz * f2 - my *  df2x;
                  MX1 =       mx * f1 ;
                  MY1 =-pz * f3 + my *  df3x;
                  MZ1 = py * f3 + mz *  df3x;

                  X2 = px * f4 + 0;
                  Y2 = py * f5 + mz *  df5x;
                  Z2 = pz * f5 - my *  df5x;
                  MX2 =       mx * f4 ;
                  MY2 =-pz * f6 + my *  df6x;
                  MZ2 = py * f6 + mz *  df6x;

                  if(task == STRESS) { /* +ve  simply support moment at center */
                     MCY = -pz * (L -a)/2 ;   
                     MCZ = -py * (L -a)/2 ;   
                  }
              }

              /* Add Contributation to Total Equivalent Load */

                 load[0] = load[0] + X1;
                 load[1] = load[1] + Y1;
                 load[2] = load[2] + Z1;
                 load[3] = load[3] + MX1;
                 load[4] = load[4] + MY1;
                 load[5] = load[5] + MZ1;
                 load[6] = load[6] + X2;
                 load[7] = load[7] + Y2;
                 load[8] = load[8] + Z2;
                 load[9] = load[9] + MX2;
                 load[10] = load[10] + MY2;
                 load[11] = load[11] + MZ2;

                 if(task == STRESS) { /* mid pt moment */
                    MCYT = MCYT+ MCY;
                    MCZT = MCZT+ MCZ;
                 }
              }
              break;
         case STRESS_LOAD:
              for(i = 0; i <= 11; i++) 
              load[i] = p->nodal_loads[i].value;
              break;
         default:
              break;
    }
  
    /* Rotate to Global */
 
    rmat = (double **) MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
    rmat = (double **) tmat(rmat,6,p);
    rt   = (double **)MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
    for( i=1 ; i<=p->size_of_stiff ; i++ )
       for( j=1 ; j<=p->size_of_stiff ; j++ )
          rt[j-1][i-1] = rmat[i-1][j-1];

    for (inc=1; inc <= p->size_of_stiff; inc++){
         p->nodal_loads[inc-1].value = 0.0;
         for (j=1; j <= p->size_of_stiff; j++)
              p->nodal_loads[inc-1].value = p->nodal_loads[inc-1].value + rt[inc-1][j-1]* (double) load[j-1];
    }

    /* mid pt moment */
/*
    if(task == STRESS){
       p->nodal_loads[13].value = MCYT; 
       p->nodal_loads[14].value = MCZT; 
    }
*/
 
    MatrixFreeIndirectDouble(rmat, p->size_of_stiff);
    MatrixFreeIndirectDouble(rt, p->size_of_stiff);

    return(p);
}