由于C++没有封装矩阵类,所以还是用到了《计算机常用数值算法与程序》(C++)一书中的头文件“Matrix.h”。对《有限单元法》书中的例题进行了编程验算,编程水平太菜。程序冗杂得不行了...
题目:给出了一个勾三股四的单元三角形,弹性模量E和泊松比v已知,单元厚度t=1。求单元的刚度矩阵;
思路:按照书本中的步骤,基本是单刚矩阵Ke=B'DBtA;其中,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;
原文链接: https://www.cnblogs.com/mt-luo/p/4153546.html
欢迎关注
微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍
原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/157520
非原创文章文中已经注明原地址,如有侵权,联系删除
关注公众号【高性能架构探索】,第一时间获取最新文章
转载文章受原作者版权保护。转载请注明原作者出处!