/* * =========================================================================== * Use Newmark Integration and modal analysis to compute dynamic time-history * response of supported cantilever beam structure. * * The cantilever is modeled with 6 beam-column elements, having 11 d.o.f. * The time-history response is estimated with "NoEigen" modes. * =========================================================================== * Copyright (C) 1996 by Mark Austin and David Mazzoni. * * This software is provided "as is" without express or implied warranty. * Permission is granted to use this software 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. * * Written By: Mark Austin August 1996 * =========================================================================== */ #include #include #include "matrix.h" #include "vector.h" #include "miscellaneous.h" #include "defs.h" #define DEBUG /* Data structures for frame geometry and elements */ typedef struct { float fXcoord; /* x coordinate */ float fYcoord; /* y coordinate */ } NODE; typedef struct { int iEnd1; /* node no at end 1 of beam element */ int iEnd2; /* node no at end 2 of beam element */ float fE; /* Young's modulus of elasticity */ float fInertia; /* Section moment of inertia */ float fMass; /* Mass per unit length (i.e. kg/m) */ } ELMT; /* * =================================================================== * Functions that define supported cantilever beam structure. * * NODE *mrfBuildNodeArray() : build array of nodal coordinates. * ELMT *mrfBuildElmtArray() : build array of structural elements. * MATRIX *mrfDestinationArray() : define destination array. * =================================================================== */ enum { NoNodes = 7, NoElmts = 6, NoDofPerElmt = 4 }; NODE *mrfBuildNodeArray() { NODE *spNode; int ii; spNode = (NODE *) safeCalloc( NoNodes, sizeof(NODE), __FILE__, __LINE__ ); for (ii = 1; ii <= NoNodes; ii = ii + 1) { spNode[ ii-1 ].fXcoord = (float) 4*(ii-1); /* nodes space 4 m apart */ spNode[ ii-1 ].fYcoord = 0.0; } return ( spNode ); } ELMT *mrfBuildElmtArray() { ELMT *spElmt; int ii; spElmt = (ELMT *) safeCalloc( NoElmts, sizeof(ELMT), __FILE__, __LINE__ ); for (ii = 1; ii <= NoElmts; ii = ii + 1) { spElmt[ ii-1 ].iEnd1 = ii; /* node at end 1 */ spElmt[ ii-1 ].iEnd2 = ii+1; /* node at end 2 */ spElmt[ ii-1 ].fE = 2.00e+11; /* E = 200 GPa */ spElmt[ ii-1 ].fInertia = 6.60e-04; /* I = 6.6 e-4 m^4 */ spElmt[ ii-1 ].fMass = 200.0; /* m = 200 kg/m */ } return ( spElmt ); } MATRIX *mrfDestinationArray() { MATRIX *spDest; spDest = matAllocIndirect("Destination Array", IntegerArray, NoElmts , NoDofPerElmt); spDest->uMatrix.ipp[ 0 ][ 0 ] = 0; spDest->uMatrix.ipp[ 0 ][ 2 ] = 1; spDest->uMatrix.ipp[ 0 ][ 1 ] = 0; spDest->uMatrix.ipp[ 0 ][ 3 ] = 2; spDest->uMatrix.ipp[ 1 ][ 0 ] = 1; spDest->uMatrix.ipp[ 1 ][ 2 ] = 3; spDest->uMatrix.ipp[ 1 ][ 1 ] = 2; spDest->uMatrix.ipp[ 1 ][ 3 ] = 4; spDest->uMatrix.ipp[ 2 ][ 0 ] = 3; spDest->uMatrix.ipp[ 2 ][ 2 ] = 5; spDest->uMatrix.ipp[ 2 ][ 1 ] = 4; spDest->uMatrix.ipp[ 2 ][ 3 ] = 6; spDest->uMatrix.ipp[ 3 ][ 0 ] = 5; spDest->uMatrix.ipp[ 3 ][ 2 ] = 7; spDest->uMatrix.ipp[ 3 ][ 1 ] = 6; spDest->uMatrix.ipp[ 3 ][ 3 ] = 8; spDest->uMatrix.ipp[ 4 ][ 0 ] = 7; spDest->uMatrix.ipp[ 4 ][ 2 ] = 9; spDest->uMatrix.ipp[ 4 ][ 1 ] = 8; spDest->uMatrix.ipp[ 4 ][ 3 ] = 10; spDest->uMatrix.ipp[ 5 ][ 0 ] = 9; spDest->uMatrix.ipp[ 5 ][ 2 ] = 0; spDest->uMatrix.ipp[ 5 ][ 1 ] = 10; spDest->uMatrix.ipp[ 5 ][ 3 ] = 11; return ( spDest ); } /* * ==================================================================================== * main() : Use direct stiffness method to compute joint displacements and distribution * of bending moments and axial forces in moment resitatant frame. * ==================================================================================== */ /* Declarations for structural analysis/finite element functions */ MATRIX * feElementStiffnessMatrix( double, double, double ); MATRIX * feElementMassMatrix( double, double ); int feSizeOfGlobalStiffness ( MATRIX * ); VECTOR * feSkylineProfile( MATRIX *, int ); enum { NoWorkingPointers = 23, NoEigen = 2, NoTimeSteps = 200 }; int main( void ) { MATRIX *spStiff, *spGenStiff, *spKe, *spK, *spEload; MATRIX *spMass, *spGenMass, *spMe, *spM, *spPt; MATRIX *spEigenvalue, *spEigenvector; MATRIX *spDest, *spLoad; MATRIX *spDispl, *spVelocity, *spAccel; MATRIX *spDisplNew, *spVelocityNew, *spAccelNew; MATRIX *spResponse, *spR, *spMasshat; MATRIX *spWork [ NoWorkingPointers ]; VECTOR *spPivot, *spProfile; NODE *spNode; ELMT *spElmt; float fE, fInertia, fDx, fElmtLength; int iRowNo, iColumnNo; int iMin, iMax, iMinRow, iSize; int ii, ij, ik, iNodeEnd1, iNodeEnd2; double dTime = 0.02; double dBeta = 0.25; double dGamma = 0.50; /* [a] : Build nodal coordinate, element property, and destination arrays */ spNode = mrfBuildNodeArray(); spElmt = mrfBuildElmtArray(); spDest = mrfDestinationArray(); /* [b] : Compute (and print) skyline profile for Global Stiffness Matrix */ iSize = feSizeOfGlobalStiffness( spDest ); spProfile = feSkylineProfile( spDest, iSize ); /* [c] : Assemble global mass and stiffness matrices */ spStiff = matAllocSkyline("Global Stiffness", DoubleArray, iSize, iSize, spProfile->uVector.ip ); spMass = matAllocSkyline("Global Mass", DoubleArray, iSize, iSize, spProfile->uVector.ip ); for( ii = 1; ii <= NoElmts; ii = ii + 1) { iNodeEnd1 = spElmt[ii-1].iEnd1; iNodeEnd2 = spElmt[ii-1].iEnd2; fElmtLength = spNode[ iNodeEnd2-1 ].fXcoord - spNode[ iNodeEnd1-1 ].fXcoord; spKe = feElementStiffnessMatrix( spElmt[ii-1].fE,spElmt[ii-1].fInertia,fElmtLength ); spMe = feElementMassMatrix( spElmt[ii-1].fMass, fElmtLength ); for( ij = 1; ij <= NoDofPerElmt; ij = ij + 1) { iRowNo = spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]; if (iRowNo != 0) { for( ik = ij; ik <= NoDofPerElmt; ik = ik + 1) { iColumnNo = spDest->uMatrix.ipp[ ii-1 ][ ik-1 ]; if (iColumnNo != 0) { iMin = MIN( iRowNo, iColumnNo); iMax = MAX( iRowNo, iColumnNo); spStiff->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] += spKe->uMatrix.dpp[ ij-1 ][ ik-1 ]; spMass->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] += spMe->uMatrix.dpp[ ij-1 ][ ik-1 ]; } } } } } /* matPrint( spStiff ); matPrint( spMass ); */ /* [d] : Compute "NoEigen" eigenvalues/eigenvectors */ spEigenvalue = matAllocIndirect( "Eigenvalues", DoubleArray, NoEigen, 1 ); spEigenvector = matAllocIndirect("Eigenvectors", DoubleArray, iSize, NoEigen ); matSubspace( spStiff, spMass, spEigenvalue, spEigenvector , NoEigen ); printf("\n*** Mode Shapes and Natural Periods of Vibration \n"); matPrint ( spEigenvector ); printf("\n"); for( ii = 1; ii <= NoEigen; ii = ii + 1) { printf(" *** Mode %2d : T = %10.5f\n", ii, 2*PI/sqrt(spEigenvalue->uMatrix.dpp[ ii-1 ][ 0 ])); } /* [e] : Initialize working displacement, velocity and acceleration vectors */ spDispl = matAllocIndirect( "Displacement", DoubleArray, NoEigen, 1 ); spVelocity = matAllocIndirect( "Velocity", DoubleArray, NoEigen, 1 ); spAccel = matAllocIndirect( "Acceleration", DoubleArray, NoEigen, 1 ); spPt = matAllocIndirect("External Load P(t)", DoubleArray, iSize, 1 ); /* [f] : Allocate matrix to store time-history response * Column 1 = time (sec). * Column 2 = Displacement of roof (m). * Column 3 = Energy (J). */ spResponse = matAllocIndirect( "Response", DoubleArray, (NoTimeSteps + 1), 3 ); /* [g] : Compute generalised mass and stiffness matrices */ spWork[0] = matTranspose ( spEigenvector ); spWork[1] = matMult ( spStiff, spEigenvector ); spWork[2] = matMult ( spMass, spEigenvector ); spGenStiff = matMult ( spWork[0], spWork[1] ); spGenStiff->cpName = saveString("Generalised Stiffness", __FILE__, __LINE__ ); spGenMass = matMult ( spWork[0], spWork[2] ); spGenMass->cpName = saveString("Generalised Mass", __FILE__, __LINE__ ); matPrint (spGenStiff); matPrint (spGenMass); /* [h] : Compute effective mass matrix */ spWork[3] = matScale (spGenStiff, dTime*dTime*dBeta ); spMasshat = matAdd (spGenMass, spWork[3]); /* [i] : Compute LU decompositio of spMasshat[] */ spPivot = SetupPivotVector( spMasshat ); spMasshat = matLUDecompositionIndirect( spMasshat, spPivot); /* [j] : Enter main loop of Newmark algorithm */ for(ii = 1; ii <= NoTimeSteps; ii++) { /* [j.1] : Update external load vector */ spPt->uMatrix.dpp[ 4 ][ 0 ] = 0.0; if (ii < 51) spPt->uMatrix.dpp[ 4 ][ 0 ] = -1e+4; spEload = matMult (spWork[0], spPt ); spWork[4] = matScale (spAccel, (dTime*dTime/2.0)*(1-2*dBeta) ); spWork[5] = matScale (spVelocity, dTime ); spWork[6] = matAdd (spDispl, spWork[5] ); spWork[7] = matAdd (spWork[6], spWork[4] ); spWork[8] = matMult (spGenStiff, spWork[7] ); spR = matSub (spEload, spWork[8] ); /* [j.2] : Compute new acceleration, velocity and displacement */ spAccelNew = matLUSubstitutionIndirect( spPivot, spMasshat, spR ); spWork[9] = matScale (spAccel, dTime*(1-dGamma)); spWork[10] = matScale (spAccelNew, dTime*dGamma); spWork[11] = matAdd (spWork[9], spWork[10]); spVelocityNew = matAdd (spVelocity, spWork[11]); spWork[12] = matScale (spAccel, (1-2*dBeta)); spWork[13] = matScale (spAccelNew, 2*dBeta ); spWork[14] = matAdd (spWork[12], spWork[13]); spWork[15] = matScale (spWork[14], dTime*dTime/2.0 ); spWork[16] = matAdd (spWork[5], spWork[15]); spDisplNew = matAdd (spDispl, spWork[16]); /* [j.3] : Update response */ matFree( spAccel ); matFree( spVelocity ); matFree( spDispl ); spAccel = matCopy( spAccelNew ); spVelocity = matCopy( spVelocityNew ); spDispl = matCopy( spDisplNew ); /* [j.4] : Compute "kinetic + potential" energy */ spWork[17] = matTranspose ( spDispl ); spWork[18] = matMult ( spGenStiff, spDispl ); spWork[19] = matMult ( spWork[17], spWork[18] ); spWork[20] = matTranspose ( spVelocity ); spWork[21] = matMult ( spGenMass, spVelocity ); spWork[22] = matMult ( spWork[20], spWork[21] ); /* [j.5] : Compute displacement vector */ spWork[23] = matMult ( spEigenvector, spDispl ); /* [j.6] : Save components of time-history response */ spResponse->uMatrix.dpp[ ii ][ 0 ] = ii*dTime; spResponse->uMatrix.dpp[ ii ][ 1 ] = spWork[23]->uMatrix.dpp[4][0]; spResponse->uMatrix.dpp[ ii ][ 2 ] = 0.5*(spWork[19]->uMatrix.dpp[0][0] + spWork[22]->uMatrix.dpp[0][0] ); /* [j.7] : Cleanup memory before moving onto next integration step */ matFree (spEload); matFree (spR); for(ij = 4; ij < NoWorkingPointers; ij++) { matFree (spWork [ij]); } } /* [k] : Print matrix of response components */ matPrint (spResponse ); } /* * ==================================================================== * feSizeOfGlobalStiffness() : Compute size of global stiffness matrix. * * Input : MATRIX * -- destination array. * Output : int iSize -- no columns/rows in global stiffness matrix. * ==================================================================== */ int feSizeOfGlobalStiffness( MATRIX *spDest ) { int ii, ij, iSize = 0; for( ii = 1; ii <= spDest->iNoRows; ii = ii + 1) { for( ij = 1; ij <= spDest->iNoColumns; ij = ij + 1) { iSize = MAX(iSize, spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]); } } return iSize; } /* * ================================================================================= * feSkylineProfile() : Compute size and skyline profile of global stiffness matrix. * * Input : MATRIX * -- destination array. * : int iSize -- size of Global Stiffness Matrix. * Output : VECTOR * -- height of each column in skyline profile. * ================================================================================= */ VECTOR * feSkylineProfile( MATRIX *spDest, int iSize ) { VECTOR * spProfile; int ii, ij, iMinRow, iColumnNo; spProfile = vecAlloc( "Skyline Profile", IntegerArray, iSize); for( ii = 1; ii <= spDest->iNoRows; ii = ii + 1) { iMinRow = MAX(1, spDest->uMatrix.ipp[ ii-1 ][ 0 ]); for( ij = 2; ij <= NoDofPerElmt; ij = ij + 1) { if( spDest->uMatrix.ipp[ ii-1 ][ ij-1 ] != 0) iMinRow = MIN(iMinRow, spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]); } for( ij = 1; ij <= NoDofPerElmt; ij = ij + 1) { iColumnNo = spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]; if((1 + iColumnNo - iMinRow) > spProfile->uVector.ip[ iColumnNo-1 ]) spProfile->uVector.ip[ iColumnNo-1 ] = 1 + iColumnNo - iMinRow; } } return spProfile; } /* * ======================================================================== * feElementStiffnessMatrix() : compute (4x4) beam-column stiffness matrix. * * Input : dE -- Young's modulus of elasticity. * : dInertia -- Section moment of inertia. * : dLength -- Element length. * Output : MATRIX * -- (4x4) beam-column stiffness matrix. * ======================================================================== */ MATRIX *feElementStiffnessMatrix( double dE, double dInertia, double dLength ) { MATRIX *spKe; double dTemp; spKe = matAllocIndirect("[K]" , DoubleArray, NoDofPerElmt, NoDofPerElmt); dTemp = 2.0*dE*dInertia/dLength; spKe->uMatrix.dpp[ 0 ][ 0 ] = dTemp*6.0/dLength/dLength; spKe->uMatrix.dpp[ 1 ][ 0 ] = dTemp*3.0/dLength; spKe->uMatrix.dpp[ 2 ][ 0 ] = -dTemp*6.0/dLength/dLength; spKe->uMatrix.dpp[ 3 ][ 0 ] = dTemp*3.0/dLength; spKe->uMatrix.dpp[ 0 ][ 1 ] = dTemp*3.0/dLength; spKe->uMatrix.dpp[ 1 ][ 1 ] = dTemp*2.0; spKe->uMatrix.dpp[ 2 ][ 1 ] = -dTemp*3.0/dLength; spKe->uMatrix.dpp[ 3 ][ 1 ] = dTemp; spKe->uMatrix.dpp[ 0 ][ 2 ] = -dTemp*6.0/dLength/dLength; spKe->uMatrix.dpp[ 1 ][ 2 ] = -dTemp*3.0/dLength; spKe->uMatrix.dpp[ 2 ][ 2 ] = dTemp*6.0/dLength/dLength; spKe->uMatrix.dpp[ 3 ][ 2 ] = -dTemp*3.0/dLength; spKe->uMatrix.dpp[ 0 ][ 3 ] = dTemp*3.0/dLength; spKe->uMatrix.dpp[ 1 ][ 3 ] = dTemp; spKe->uMatrix.dpp[ 2 ][ 3 ] = -dTemp*3.0/dLength; spKe->uMatrix.dpp[ 3 ][ 3 ] = dTemp*2.0; return ( spKe ); } /* * ===================================================================== * feElementMassMatrix() : compute (4x4) beam-column mass matrix. * * Input : dMass -- Mass per unit length along element (e.g. kg/m). * : dLength -- Element length. * Output : MATRIX * -- (4x4) beam-column stiffness matrix. * ===================================================================== */ MATRIX *feElementMassMatrix( double dMass, double dLength ) { MATRIX *spMe; double dTemp; spMe = matAllocIndirect("[Me]" , DoubleArray, NoDofPerElmt, NoDofPerElmt); dTemp = dMass*dLength/420; spMe->uMatrix.dpp[ 0 ][ 0 ] = dTemp*156; spMe->uMatrix.dpp[ 1 ][ 0 ] = dTemp*22*dLength; spMe->uMatrix.dpp[ 2 ][ 0 ] = dTemp*54; spMe->uMatrix.dpp[ 3 ][ 0 ] = -dTemp*13*dLength; spMe->uMatrix.dpp[ 0 ][ 1 ] = dTemp*22*dLength; spMe->uMatrix.dpp[ 1 ][ 1 ] = dTemp* 4*dLength*dLength; spMe->uMatrix.dpp[ 2 ][ 1 ] = dTemp*13*dLength; spMe->uMatrix.dpp[ 3 ][ 1 ] = -dTemp* 3*dLength*dLength; spMe->uMatrix.dpp[ 0 ][ 2 ] = dTemp*54; spMe->uMatrix.dpp[ 1 ][ 2 ] = dTemp*13*dLength; spMe->uMatrix.dpp[ 2 ][ 2 ] = dTemp*156; spMe->uMatrix.dpp[ 3 ][ 2 ] = -dTemp*22*dLength; spMe->uMatrix.dpp[ 0 ][ 3 ] = -dTemp*13*dLength; spMe->uMatrix.dpp[ 1 ][ 3 ] = -dTemp* 3*dLength*dLength; spMe->uMatrix.dpp[ 2 ][ 3 ] = -dTemp*22*dLength; spMe->uMatrix.dpp[ 3 ][ 3 ] = dTemp* 4*dLength*dLength; return ( spMe ); }