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