六结点矩形单元的有限元公式推导word资料13页.docx
- 文档编号:5139199
- 上传时间:2022-12-13
- 格式:DOCX
- 页数:10
- 大小:53.16KB
六结点矩形单元的有限元公式推导word资料13页.docx
《六结点矩形单元的有限元公式推导word资料13页.docx》由会员分享,可在线阅读,更多相关《六结点矩形单元的有限元公式推导word资料13页.docx(10页珍藏版)》请在冰豆网上搜索。
六结点矩形单元的有限元公式推导word资料13页
六结点矩形单元的有限元公式推导
一、a型6结点矩形单元
1.1单元位移模式及插值函数的构造
首先引入局部坐标
、
其中
利用教材中Serendipity四边形单元的构造方法,假设开始只有4的个角结点,对应这些结点的插值函数可以用双一次拉格朗日多项式构造,即
当增加结点5、结点6时,与它对应的插值函数可表示为y方向2次,x方向1次拉格朗日多项式的乘积,即
此时
在结点5、结点6处不为0,为了使其继续满足这个要求,将其修正为
可以证明这样构造的插值函数满足
和
这些基本要求,继而可将位移函数表示成结点位移的函数
矩阵形式为
1.2应变矩阵的推导
1.3应力矩阵的推导
1.4单元刚度矩阵的推导
利用最小位能原理建立有限元方程,求得单元刚度矩阵为
其中
其中
形函数的偏导数如下表
二、b型6结点矩形单元
2.1单元位移模式及插值函数的构造
同样引入局部坐标
、
其中
假设开始只有4个角结点,对应这些结点的插值函数可以用双一次拉格朗日多项式构造,即
当增加结点5时,与它对应的插值函数可表示为y方向2次,x方向1次拉格朗日多项式的乘积,当增加结点6时,与它对应的插值函数可表示为x方向2次,y方向1次拉格朗日多项式的乘积,即
此时
在结点5、结点6处不为0,为了使其继续满足这个要求,将其修正为
可以证明这样构造的插值函数满足
和
这些基本要求,继而可将位移函数表示成结点位移的函数
矩阵形式为
2.2应变矩阵的推导
2.3应力矩阵的推导
2.4单元刚度矩阵的推导
利用最小位能原理建立有限元方程,求得单元刚度矩阵为
其中
其中
其中形函数的偏导数如下表
三、a型6结点矩形单元与4结点矩形单元的计算比较
3.1计算示例与参数
为了验证上文中6结点矩形单元有限元公式的正确性,此处对图示结构用6结点矩形单元与4结点矩形单元分别进行受力分析。
用6结点矩形单元进行有限元分析的程序用C语言自行编写,而用4结点矩形单元进行有限元分析则在ANSYS中进行。
图示结构为一正方形薄板,边长为10m,左下角、右下角受到固定约束,左上角受到1N的拉力,右上角收到1N的压力,共划分为100个单元。
弹性模量E、薄板厚度t均取为1,泊松比取为0.25。
具体图示如下。
4结点矩形单元受力分析示意图(利用ANSYS)
a型6结点矩形单元受力分析示意图
3.2计算结果及对比
利用ANSYS软件计算的位移结果如下图所示:
由于使用自己编写的程序计算的位移结果数据量太大,所以在此仅列出若干结点的位移值并与ANSYS的计算结果相比较。
典型节点如下图所示。
在典型结点处,使用6结点单元和4结点单元求得的位移结果的比较见下表。
X方向位移
Y方向位移
结点编号
6结点
4结点(ANSYS)
误差
6结点
4结点
(ANSYS)
误差
1
13.871
12.689
1.182
13.871
12.689
1.182
2
7.4918
6.7298
0.762
4.9145
4.5615
0.353
3
4.2021
3.6963
0.5058
1.4431
1.3240
0.1191
4
1.4279
1.1584
0.2695
-1.3311
-1.2139
-0.1172
5
-0.83084
-0.88378
0.05294
-3.4081
-3.0521
-0.356
6
13.871
12.689
1.182
-13.871
-12.689
-1.182
7
7.4918
6.7298
0.762
-4.9145
-4.5615
-0.353
8
4.2021
3.6963
0.5058
-1.4431
-1.3240
-0.1191
9
1.4279
1.1584
0.2695
1.3311
1.2139
0.1172
10
-0.83084
-0.88378
0.05294
3.4082
3.0521
0.3561
3.3结果分析
从计算得到的位移结果可以发现,使用6结点与4结点单元进行分析得到的位移误差较小,这验证了本文之前所推导得到的有限元公式的正确性。
但是,两者之前仍然具有一定的误差。
这可能是由于划分单元的个数过小造成的。
四、6结点单元C语言计算程序
#include
#include
#include
FILE*fpElist,*fpNlist,*fpDlist,*fpFlist,*fpOutput;
FILE*fpNodeData,*fpElemData;//输出文件,用于ANSYS的输入
#defineNMAX1000//最大单元数减一
doublet=1;//薄板厚度
doubleEy=1;//杨氏模量
doublev=0.25;//泊松比
inti,j,k;//循环变量
doubledbTemp,dbTemp1;//临时变量
intNumOfNode;//节点总数
intE[NMAX+1][7],NumOfElem;//单元信息矩阵,单元总数
doubleK[NMAX+1][NMAX+1],kk[13][13];//总体刚度矩阵,单元刚度矩阵
intD[NMAX*3];//节点位移约束矩阵,1表示位移为0,0表示位移不受约束;排列顺序为:
节点1的X方向,节点1的Y方向,节点2的X方向。
。
。
doubleF[NMAX*3];//节点荷载矩阵,排列顺序为:
节点1的X方向受力,节点1的Y方向受力,节点2的X方向受力。
。
。
intNumOfVari=0,V[NMAX*3];//自由变量总数,自由变量->所有变量
doubleU[NMAX*3];//节点位移变量,结果
doubleaa[6][6]={
{1./5,1./30,-1./30,-1./5,-1./5,1./5},
{1./30,1./5,-1./5,-1./30,-1./5,1./5},
{-1./30,-1./5,1./5,1./30,1./5,-1./5},
{-1./5,-1./30,1./30,1./5,1./5,-1./5},
{-1./5,-1./5,1./5,1./5,8./15,-8./15},
{1./5,1./5,-1./5,-1./5,-8./15,8./15}};
doublebb[6][6]={
{1./4,5./12,5./12,1./4,-1./3,-1./3},
{-5./12,-1./4,-1./4,-5./12,1./3,1./3},
{5./12,1./4,1./4,5./12,-1./3,-1./3},
{-1./4,-5./12,-5./12,-1./4,1./3,1./3},
{1./3,-1./3,-1./3,1./3,0,0},
{-1./3,1./3,1./3,-1./3,0,0}};
doublecc[6][6]={
{1./4,-5./12,5./12,-1./4,1./3,-1./3},
{5./12,-1./4,1./4,-5./12,-1./3,1./3},
{5./12,-1./4,1./4,-5./12,-1./3,1./3},
{1./4,-5./12,5./12,-1./4,1./3,-1./3},
{-1./3,1./3,-1./3,1./3,0,0},
{-1./3,1./3,-1./3,1./3,0,0}};
doubledd[6][6]={
{19./9,13./9,13./18,19./18,-16./9,-8./9},
{13./9,19./9,19./18,13./18,-16./9,-8./9},
{13./18,19./18,19./9,13./9,-8./9,-16./9},
{19./18,13./18,13./9,19./9,-8./9,-16./9},
{-16./9,-16./9,-8./9,-8./9,16./9,8./9},
{-8./9,-8./9,-16./9,-16./9,8./9,16./9}};
doublea_b=1;//矩形单元的长宽比a/b
voidmain()
//建立单元信息矩阵,100个单元,11*21个结点
NumOfElem=100;
NumOfNode=11*21;
for(i=0;i<=9;i++)
for(j=1;j<=10;j++)
E[i*10+j][1]=22*i+j+23;
E[i*10+j][2]=22*i+j+1;
E[i*10+j][3]=22*i+j;
E[i*10+j][4]=22*i+j+22;
E[i*10+j][5]=22*i+j+12;
E[i*10+j][6]=22*i+j+11;
for(i=1;i<=NumOfElem;i++)
printf("%d)%d%d%d%d%d%d\n",i,E[i][1],E[i][2],E[i][3],E[i][4],E[i][5],E[i][6]);
//计算单元刚度矩阵
for(i=1;i<=6;i++)
for(j=1;j<=6;j++)
kk[2*i-1][2*j-1]=aa[i-1][j-1]+(1-v)/2*dd[i-1][j-1];
kk[2*i-1][2*j]=v*bb[i-1][j-1]+(1-v)/2*cc[i-1][j-1];
kk[2*i][2*j-1]=v*cc[i-1][j-1]+(1-v)/2*bb[i-1][j-1];
kk[2*i][2*j]=dd[i-1][j-1]+(1-v)/2*aa[i-1][j-1];
printf("单刚\n");
for(i=1;i<=12;i++)
for(j=1;j<=12;j++)
kk[i][j]=kk[i][j]*Ey/(1-v*v);
printf("%lf",kk[i][j]);
printf("\n");
//输出单元刚度矩阵的结果文件
if((fpOutput=fopen("DanGang.txt","w+"))==NULL)
printf("Cannotopenthisfile.\n");
exit(0);
for(i=1;i<=12;i++)
for(j=1;j<=12;j++)
fprintf(fpOutput,"%7.4lf",kk[i][j]);
fprintf(fpOutput,"\n");
fclose(fpOutput);
//根据单元刚度矩阵装配整体刚度矩阵
for(k=1;k<=NumOfElem;k++)
for(i=1;i<=6;i++)
for(j=1;j<=6;j++)
K[E[k][i]*2-1][E[k][j]*2-1]+=kk[i*2-1][j*2-1];
K[E[k][i]*2-1][E[k][j]*2]+=kk[i*2-1][j*2];
K[E[k][i]*2][E[k][j]*2-1]+=kk[i*2][j*2-1];
K[E[k][i]*2][E[k][j]*2]+=kk[i*2][j*2];
//输出整体刚度矩阵的结果文件
if((fpOutput=fopen("ZhengGang.txt","w+"))==NULL)
printf("Cannotopenthisfile.\n");
exit(0);
for(i=1;i<=2*NumOfNode;i++)
for(j=1;j<=2*NumOfNode;j++)
fprintf(fpOutput,"%7.4lf",K[i][j]);
fprintf(fpOutput,"\n");
fclose(fpOutput);
//设置节点位移约束矩阵
D[1]=1;D[2]=1;
D[221*2-1]=1;D[221*2]=1;
//设置节点荷载矩阵
F[11*2-1]=1;
F[231*2-1]=-1;
//解线性方程组
//将0位移分量的平衡方程划去
NumOfVari=0;//自由变量总数
for(i=1;i<=2*NumOfNode;i++)
if(D[i]==0)
NumOfVari++;
V[NumOfVari]=i;
//利用高斯消元法求解方程组:
分两步
//第一步消元
for(i=1;i<=NumOfVari-1;i++)
if(K[V[i]][V[i]]<0.00001)
printf("主元为零\n");
for(j=i+1;j<=NumOfVari;j++)
dbTemp=K[V[j]][V[i]]/K[V[i]][V[i]];
for(k=i;k<=NumOfVari;k++)
if(k==i)
K[V[j]][V[k]]=0;
else
K[V[j]][V[k]]=K[V[j]][V[k]]-dbTemp*K[V[i]][V[k]];
F[V[j]]=F[V[j]]-dbTemp*F[V[i]];
//第二步回代
for(i=1;i<=2*NumOfNode;i++)
U[i]=0;
for(i=NumOfVari;i>=1;i--)
U[V[i]]=F[V[i]];
for(j=NumOfVari;j>i;j--)
U[V[i]]-=K[V[i]][V[j]]*U[V[j]];
U[V[i]]=U[V[i]]/K[V[i]][V[i]];
//输出节点位移文件NodeData.txt
if((fpNodeData=fopen("NodeData.txt","w+"))==NULL)
printf("Cannotopenthisfile.\n");
exit(0);
fprintf(fpNodeData,"%4d\n",NumOfNode);
for(i=1;i<=NumOfNode;i++)
fprintf(fpNodeData,"%4d,%10.6lf,%10.6lf\n",i,U[i*2-1],U[i*2]);
//节点编号X方向位移Y轴方向位移
希望以上资料对你有所帮助,附励志名言3条:
:
1、世事忙忙如水流,休将名利挂心头。
粗茶淡饭随缘过,富贵荣华莫强求。
2、“我欲”是贫穷的标志。
事能常足,心常惬,人到无求品自高。
3、人生至恶是善谈人过;人生至愚恶闻己过。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 结点 矩形 单元 有限元 公式 推导 word 资料 13