轴力杆件有限元分析
#include
#include
#include
#include
/////////////////////单元、节点与荷载类的定义/////////////////////////////////////////////
// 节点类
class Node
{
public:
int nd; // 表示节点号，从0开始编号
int id1; // 表明节点的x方向的位移的约束性质，为1则被约束，为0则自由
int id2; // 表明节点的y方向的位移的约束性质，为1则被约束，为0则自由
double x; // 表明节点的x坐标
double y; // 表明节点的y坐标
void Boundary(int nNode); // 边界条件处理
}
// 杆件单元类
class Truss
{
public:
int ne; // 单元号
int nd1; // 单元节点1
int nd2; // 单元节点2
double E; // 弹性模量
double A; // 截面面积
void Stiff(); // 计算单元刚度矩阵
void Assemble(); // 将单元刚度矩阵叠加到结构刚度阵
double Force(); // 计算杆件的内力
private:
static double k[4][4]; // 保存单元刚度矩阵
}
// 荷载类
class Load
{
public：
int nl; // 荷载编号
int nd; // 荷载所作用节点
int dir; // 荷载方向 弯=1,；X方向=2；Y方向=3
double value; // 荷载大小
void Assemble(); // 组装荷载矩阵
}
/////////////////////有限元数据的储存空间/////////////////////////////////////////////
Node nodes[250]; // 输入节点数据
Truss elements[1000]; // 输入单元数据
Load loads[500]; // 输入荷载数据
double K[500][500]; // 结构刚度矩阵
double load[500]; // 结构等效节点力
double disp[500]; // 计算出来的节点位移
/////////////////////单元类的实现/////////////////////////////////////////////
// 计算杆件单元的长度
/* double Truss::Length()
{
double dx=nodes[nd2].x-nodes[nd1].x;
double dy=nodes[nd2].y-nodes[nd1].y;
return sqrt(dx*dx+dy*dy);
}*/
// 计算单元刚度矩阵
void Truss::Stiff()
{
double dx=nodes[nd2].x-nodes[nd1].x;
double dy=nodes[nd2].y-nodes[nd1].y;
double len=(dx*dx+dy*dy);
double cos=dx/len;
double sin=dy/len;
double EAL=E*A/len;
k[0][0]=cos*cos*EAL;
k[1][1]=sin*sin*EAL;
k[0][1]=cos*sin*EAL;
k[1][0]=cos*sin*EAL;
for(int i=0;i<2;i++) // 刚度阵矩阵转置
{
for(int j=0;j<2;j++)
{
k[i][j+2]=-k[i][j];
k[i+2][j]=-k[i][j];
k[i+2][j+2]=k[i][j];
}
}
}
// 将单元刚度阵叠加到结构刚度阵
void Truss::Assemble()
{
Stiff(); // 计算单元刚度矩阵
// 获得方程编号
int id[4];
id[0]=nodes[nd1].nd*2;
id[1]=id[0]+1;
id[2]=nodes[nd2].nd*2;
id[3]=id[2]+1;
for(int i=0;i<4;i++)
{
for(int j=0;j<4;j++)
{
int row=id[i];
int col=id[j];
K[row][col]+=k[i][j];
}
}
}
// 计算杆件内力
double Truss::Force()
{
double dx=nodes[nd2].x-nodes[nd1].x;
double dy=nodes[nd2].y-nodes[nd1].y;
double len= sqrt(dx*dx+dy*dy);
d
ouble cos=dx/len;
double sin=dy/len;
double EAL=E*A/len;
double u1=cos*disp[2*nd1]+sin*disp[2*nd1+1];
double u2=cos*disp[2*nd2]+sin*disp[2*nd2+1];
return EAL*(u2-u1);
}
//////////////////////荷载类的实现//////////////////////////////////////////
// 叠加形成的荷载矩阵
void Load::Assemble()
{
int nfun=nd*2+dir;
load[nfun]+=value;
}
/////////////////////高斯消去法求解方程//////////////////////////////////////
//高斯消去法求解方程
void Guass(int neq)
{
for(int s=0;s
disp[s]=load[s];
}
// 消元
for(s=0;s
for(int j=s+1;j
for(int i=s+1;i
for(int j=s+1;j
}
}
// 回代
for(int i=neq-1;i>=0;i--)
{
for(int j=i+1;j
}
///////////////////////节点类的实现///////////////////////////////////////////////
void Node::Boundary(int nNode)
{
int neq=2*nNode;
if(1==id1)
{
for(int i=0;i
load[2*nd]=0.0;
}
if(1==id2)
{
for(int i=1;i
load[2*nd]=0.0;
}
}
/////////////////////////////////桁架有限元主程序////////////////////////////////////////
void main(int argc,char** argv)
{
int nNode; // 节点数目
int nEle; // 单元数目
int nLoad; // 节点荷载数目
int nd; // 节点号
int ne; // 单元号
int nl; // 节点荷载号
// 输入节点数据
cin>>nNode>>nEle>>nLoad;
cout<<"节点数目="<
for(int i=0;i
cin>>nd;
nodes[nd].nd=nd;
cin>>nodes[nd].id1>>nodes[nd].id2>>nodes[nd].x>>nodes[nd].y;
cout<< nodes[nd].nd<<"\t"<
// 输入单元数据
cout<<"\n单元数据\n 单元号\t 弹性模量\t 面积\t 左节点\t 右节点\n";
for(i=0;i
cin>>ne;
elements[ne].ne=ne;
cin>>elements[ne].E>>elements[ne].A>>elements[ne].nd1>>elements[ne].nd2;
cout<< elements[ne].ne<<"\t"<
// 输入节点荷载
cout<<"\n节点荷载\n 序号\t 节点\t 方向\t 大小\n";
for(i=0;i
cin>>nl;
loads[nl].nl=nl;
cin>>loads[nl].nd>>loads[nl].dir>>loads[nl].value;
cout<< loads[nl].nl<<"\t"<
// 叠加荷载矩阵
for(i=0;i
loads[i].Assemble();
}
// 叠加结构刚度矩阵
for(i=0;i
{
elements[i].Assemble();
}
// 边界条件处理，结构刚度举证充0置1，荷载阵列置0
for(i=0;i
nodes[i].Boundary(nNode);
}
// 解方程
Gauss(2*nNode);
// 输出位移
cout<<"\n 节点位移\n 节点\tx 位移\ty 位移\n";
for(i=0;i
cout<
// 计算杆件内力
cout<<"\n 杆件内力\n 单元\t 轴力\n";
for(i=0;i
cout<
}