beamhinge.cpp

一个计算悬臂梁的有限元vc源码

// BeamHinge.cpp: implementation of the CBeamHinge class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "afxtempl.h"
#include "fstream.h"
#include "OOPFE.h"
#include "Matrix.h"
#include "SparseMatrix.h"
#include "Node.h"
#include "BaseMaterial.h"
#include "BeamHinge.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

#define BEAM_HINGE 2
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CBeamHinge::CBeamHinge()
{

}

CBeamHinge::~CBeamHinge()
{

}

void CBeamHinge::GetStiffness()
{
	double dBuf;
	double dE;
	CTypedPtrArray<CPtrArray,CBaseMaterial*>& apMaterial=*m_papMaterial;

	m_matT=0.0;
	m_matT(2,2)=m_matT(5,5)=1.0;	
	m_matT(0,0)=m_matT(1,1)=m_matT(3,3)=m_matT(4,4)=m_dCos;
	m_matT(0,1)=m_matT(3,4)=m_dSin;
	m_matT(1,0)=m_matT(4,3)=-m_dSin;
	m_matTT.Trans(m_matT);

	dE=apMaterial[m_iMaterialIndex]->GetE();

	m_matKe=0.0;
	dBuf=dE*m_dArea/m_dLength;
	m_matKe(0,0)=m_matKe(3,3)=dBuf;
	m_matKe(0,3)=m_matKe(3,0)=-dBuf;
	dBuf=m_dLength*m_dLength*m_dLength;
	dBuf=3.0*dE*m_dInertia/dBuf;
	m_matKe(1,1)=m_matKe(4,4)=dBuf;
	m_matKe(1,4)=m_matKe(4,1)=-dBuf;
	dBuf=m_dLength*m_dLength;
	dBuf=3.0*dE*m_dInertia/dBuf;
	m_matKe(1,5)=m_matKe(5,1)=dBuf;
	m_matKe(4,5)=m_matKe(5,4)=-dBuf;
	dBuf=3.0*dE*m_dInertia/m_dLength;
	m_matKe(5,5)=dBuf;

	m_matKe=m_matTT*m_matKe;
	m_matKe*=m_matT;
}

int CBeamHinge::GetElementType()
{
	return BEAM_HINGE;
}

void CBeamHinge::CalcuInternalForce(const double* adNodeDisp)
{
	int loop,iBuf;
	double dx,dDx,dMi,dMj,dBuf,dBuf1;
	double dKxi;

	for(loop=0;loop<2;loop++){
		iBuf=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop,0)=adNodeDisp[iBuf];
		iBuf=m_pNode->GetYDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop+1,0)=adNodeDisp[iBuf];
		iBuf=m_pNode->GetRDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop+2,0)=adNodeDisp[iBuf];
	}
	
	GetStiffness();
	m_matNodeDisp=m_matKe*m_matNodeDisp;
	
	GetTransferMatrixGCToLC(m_matT);
	m_matNodeDisp=m_matT*m_matNodeDisp;
	
	m_dN=-m_matNodeDisp(0,0);

	dMi=-m_matNodeDisp(2,0);
	dMj=m_matNodeDisp(5,0);
	m_adM[0]+=dMi;
	m_adM[10]+=dMj;
	m_adQ[0]-=m_matNodeDisp(1,0);
	m_adQ[10]-=m_matNodeDisp(1,0);

	dBuf=(dMi-dMj)/10.0;
	dx=dDx=m_dLength/10.0;
	for(loop=1;loop<10;loop++){
		m_adM[loop]+=dMi-dBuf*(double)loop;
		m_adQ[loop]-=m_matNodeDisp(1,0);
	}

	for(loop=0;loop<2;loop++){
		iBuf=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop,0)=adNodeDisp[iBuf];
		iBuf=m_pNode->GetYDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop+1,0)=adNodeDisp[iBuf];
		iBuf=m_pNode->GetRDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop+2,0)=adNodeDisp[iBuf];
	}
	GetTransferMatrixGCToLC(m_matT);
	m_matNodeDisp=m_matT*m_matNodeDisp;

	m_adFy[0]-=m_matNodeDisp(1,0);
	m_adFy[10]-=m_matNodeDisp(4,0);
	for(loop=1;loop<10;loop++){
		dKxi=loop/10.0;
		dBuf=dKxi*dKxi;
		dBuf1=dKxi*dKxi*dKxi;
		m_adFy[loop]-=(1.0-1.5*dKxi+0.5*dBuf1)*m_matNodeDisp(1,0)
					  +(1.5*dKxi-0.5*dBuf1)*m_matNodeDisp(4,0)
					  +0.5*(dBuf1-dKxi)*m_dLength*m_matNodeDisp(5,0);
	}
	m_dU=m_matNodeDisp(0,0);
	m_dU1=m_matNodeDisp(3,0);	
}

void CBeamHinge::GetMassMatrix()
{
	double dBuf;
	double dDensity;
	CTypedPtrArray<CPtrArray,CBaseMaterial*>& apMaterial=*m_papMaterial;

	m_matT=0.0;
	m_matT(2,2)=m_matT(5,5)=1.0;	
	m_matT(0,0)=m_matT(1,1)=m_matT(3,3)=m_matT(4,4)=m_dCos;
	m_matT(0,1)=m_matT(3,4)=m_dSin;
	m_matT(1,0)=m_matT(4,3)=-m_dSin;
	m_matTT.Trans(m_matT);

	dDensity=apMaterial[m_iMaterialIndex]->GetDensity();
	m_matKe=0.0;
	dBuf=dDensity*m_dArea*m_dLength/820.0;
	m_matKe(0,0)=m_matKe(3,3)=280.0*dBuf;
	m_matKe(0,3)=m_matKe(3,0)=140.0*dBuf;
	m_matKe(1,1)=198.0*dBuf;
	m_matKe(1,4)=m_matKe(4,1)=117.0*dBuf;
	m_matKe(1,5)=m_matKe(5,1)=-33.0*dBuf*m_dLength;
	m_matKe(4,4)=408.0*dBuf;
	m_matKe(4,5)=m_matKe(5,4)=-72.0*dBuf*m_dLength;
	m_matKe(5,5)=16.0*dBuf*m_dLength*m_dLength;

	m_matKe=m_matTT*m_matKe;
	m_matKe*=m_matT;
}