📄 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,
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -