首页 > 代码库 > 单元刚度矩阵 C++
单元刚度矩阵 C++
由于C++没有封装矩阵类，所以还是用到了《计算机常用数值算法与程序》（C++）一书中的头文件“Matrix.h”。对《有限单元法》书中的例题进行了编程验算，编程水平太菜。程序冗杂得不行了...
题目：给出了一个勾三股四的单元三角形，弹性模量E和泊松比v已知，单元厚度t=1。求单元的刚度矩阵；
思路：按照书本中的步骤，基本是单刚矩阵Ke=B‘*D*B*t*A;其中，B为单元应变矩阵，DB=S为单元应力矩阵。一步一步进行求解，程序如下：
#include<iostream>#include<cmath>#include"Matrix.h"//#include"Comm.h"using namespace std;void main(void){ double aijm(double x2,double x3,double y2,double y3); double bijm(double y2,double y3); double cijm(double x2,double x3); //求解参数ai,bj,am...的函数声明 double xi,xj,xm,yi,yj,ym; //i，j，m点逆时针排列 double ai,aj,am,bi,bj,bm,ci,cj,cm; //读取三角单元的三个顶点坐标 cout<<"xi,yi="; cin>>xi>>yi; cout<<"xj,yj="; cin>>xj>>yj; cout<<"xm,ym="; cin>>xm>>ym; //求解参数ai,aj,am,... ai=aijm(xj,xm,yj,ym); aj=aijm(xi,xm,yi,ym); am=aijm(xi,xj,yi,yj); bi=bijm(yj,ym); bj=bijm(ym,yi); bm=bijm(yi,yj); ci=cijm(xj,xm); cj=cijm(xi,xm); cm=cijm(xi,xj); //输出参数 cout<<"ai="<<ai<<"\taj="<<aj<<"\tam="<<am<<endl; cout<<"bi="<<bi<<"\tbj="<<bj<<"\tbm="<<bm<<endl; cout<<"ci="<<ci<<"\tcj="<<cj<<"\tcm="<<cm<<endl; /***********************************************/ double E=2*pow(10,5); //弹性模量MPa double v=0.2; //泊松比 double t=1; //单元厚度设为1 double A=3*4/2; //三角形单元的面积 /***********************************************/ /*求三角形单元的应变矩阵B=[Bi,Bj,Bm]*/ const double b[3][6]= { {bi,0,bj,0,bm,0}, {0,ci,0,cj,0,cm}, {ci,bi,cj,bj,cm,bm} }; matrix<double> B(&b[0][0],3,6); //前面没有加const，否则后面不能进行tanspose转置 B=B/(2*A); //得到应变矩阵B； cout<<endl<<"B="<<endl; MatrixLinePrint(B); //输出应变矩阵B； /*求应力矩阵S=DB=[si,sj,sm];*/ double E0,v0,CONSTANT; //平面应力问题 E0=E; v0=v; CONSTANT=E0/(2*A*(1-v0*v0)); //参量CONSTANT /************************************************* const double ssi[3][2]= { {bi,v0*ci}, {v0*bi,ci}, {(1-v0)*ci/2,(1-v0)*bi/2} }; const matrix<double> Si(&ssi[0][0],3,2); cout<<"Si="<<endl; MatrixLinePrint(Si); const double ssj[3][2]= { {bj,v0*cj}, {v0*bj,cj}, {(1-v0)*cj/2,(1-v0)*bj/2} }; const matrix<double> Sj(&ssj[0][0],3,2); cout<<"Sj="<<endl; MatrixLinePrint(Sj); const double ssm[3][2]= { {bm,v0*cm}, {v0*bm,cm}, {(1-v0)*cm/2,(1-v0)*bm/2} }; const matrix<double> Sm(&ssm[0][0],3,2); cout<<"Sm="<<endl; MatrixLinePrint(Sm); ****************************************/ //合成应力矩阵S； const double s[3][6]= { {bi,v0*ci,bj,v0*cj,bm,v0*cm}, {v0*bi,ci,v0*bj,cj,v0*bm,cm}, {(1-v0)*ci/2,(1-v0)*bi/2,(1-v0)*cj/2,(1-v0)*bj/2,(1-v0)*cm/2,(1-v0)*bm/2} }; matrix<double> S(&s[0][0],3,6); //定义应力矩阵S； S=S*CONSTANT; //得到应力矩阵S； cout<<"S="<<endl; MatrixLinePrint(S); //输出应力矩阵 /**************求B应变矩阵的转置BT********************/ matrix<double> BT(6,3); MatrixTranspose(B,BT); //将应变矩阵B转置，得到BT cout<<"BT="<<endl; MatrixLinePrint(BT); //输出转置应变矩阵 /**************求单元刚度矩阵Ke***********************/ matrix<double> Ke(6,6); //定义单元刚度矩阵 Ke=(BT*S)*t*A; //得到单元的刚度矩阵 cout<<"Ke="<<endl; MatrixLinePrint(Ke); //输出单刚矩阵}/*********************************************************/double aijm(double x2,double x3,double y2,double y3){ return x2*y3-x3*y2;}/*********************************************************/double bijm(double y2,double y3){ return y2-y3;}/*********************************************************/double cijm(double x2,double x3){ return -x2+x3;}
输入：
输出结果如下：
对比书中的结果一致，可以验证单元刚度矩阵的三个性质：对称性，奇异性，主元素大于0;
单元刚度矩阵 C++
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。