有限元法大作业.docx

有限元法大作业.docx

《有限元法大作业.docx》由会员分享，可在线阅读，更多相关《有限元法大作业.docx（48页珍藏版）》请在冰豆网上搜索。

有限元法大作业

有限元法大作业

一平面刚架的程序

用VisualC++编制的平面刚架的源程序如下：

///////////////////////////////////////////////////////程序开始//////////////////////////////////////////////////////////////////

#include"iostream.h"

#include"math.h"

#include"stdlib.h"

#include"conio.h"

//*****************

//声明必要变量

//*****************

#definePI3.141592654

intNE;//单元数

intNJ;//节点数

intNZ;//支承数

intNPJ;//有节点载荷作用的节点数

intNPF;//非节点载荷数

intHZ;//载荷码

intE;//单元码

intfangchengshu;

doubleF[303];//各节点等效总载荷数值

intdym_jdm[100][2];//单元码对应的节点码:

dym_jdm[][0],dym_jdm[][1]分为前后节点总码

intzhichengweizhi[300];//记录支持节点作用点的数组

intfjzhzuoyongdanyuan[100];//非节点载荷作用单元

intfjzhleixing[100];//非节点载荷类型:

1-均布,2-垂直集中,3-平行集中,4-力偶,5-角度集中

doublefjzhzhi[100];//非节点载荷的值

doublefjzhzuoyongdian[100];//非节点载荷在各竿的作用点

doublefjzhjiaodu[100];//非节点载荷作用角度

intjdzhzuoyongdian[100];//节点载荷作用的节点数组

doublejiedianzaihe[101][3];//节点载荷值,其jiedianzaihe[][0]--jiedianzaihe[][2]分别为U,V,M

doublezhengtigangdu[303][303];//整体刚度数组

doublechangdu[100];//各单元竿长数组

doublejiaodu[100];//各单元角度数组

doubletanxingmoliang[100];//各单元弹性模量数组

doubleJ_moliang[100];//各单元J模量数组

doublemianji[100];//各单元面积数组

doubleweiyi[303];//记录各个节点位移的数组

doubledy_weiyi[100][6];//各个单元在局部坐标系中的位移数组dy_weiyi[i][0]-dyweiyi[i][6]分别为第i+1单元的u1,v1,@1,u2,v2,@2

doubledy_neili[100][6];//各个单元在局部坐标系中的固端内力dy_weiyi[i][0]-dyweiyi[i][6]分别为第i+1单元的U1,V1,M1,U2,V2,M2

doublegan_neili[100][6];//各个单元的竿端内力数组,gan_neili[i][6]表示第i+1单元的6内力.

//*******************

//一些处理过程的函数

//*******************

voidinput（）;//数据的输入

voidzonggang（）;//计算总刚度,存放于zhengtigangdu[][]数组中

voidzongzaihe（）;//计算等效总节点载荷

voidzhichengyinru（）;//引入支承条件

voidjsweiyi（）;//求各个节点位移

voidjs_dy_weiyi（）;//求局部坐标系中的位移

voidganduanneili（）;//求竿端内力

voiddy_gangdu（inti,doubledg[6][6]）;//求单元在局部坐标系中的单刚

voidjs_T_T1（inti,doubleT[6][6],doubleT1[6][6]）;//求单元的转换矩阵及其逆阵

//************

//主函数

//************

voidmain（）

{

input（）;

cout<<"输出结果"<

cout<<"＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝”<

zonggang（）;

zongzaihe（）;

zhichengyinru（）;

jsweiyi（）;

cout<<"输出位移:

"<

for（inti=0;i

{

cout<<"u["<

"<

cout<<"v["<

"<

cout<<"@["<

"<

}

js_dy_weiyi（）;

ganduanneili（）;

//输出竿端内力

cout<<"＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝”<

cout<<"输出竿内力:

"<

for（i=0;i

{

cout<<"第"<

"<

cout<<"U1:

"<

"<

"<

cout<<"U2:

"<

"<

"<

}

cout<<"==========================================="<

cout<<"输出完毕,按任意键退出！

"<

getchi（）;

}

//***************************

//求单元在局部坐标系中的单刚

//***************************

voiddy_gangdu（inti,doubledg[6][6]）

{

dg[0][0]=tanxingmoliang[i]*mianji[i]/changdu[i];dg[0][1]=0;dg[0][2]=0;

dg[0][3]=-tanxingmoliang[i]*mianji[i]/changdu[i];dg[0][4]=0;dg[0][5]=0;dg[1][0]=0;

dg[1][1]=12*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]*changdu[i]）;

dg[1][2]=-6*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]）;dg[1][3]=0;

dg[1][4]=-12*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]*changdu[i]）;

dg[1][5]=-6*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]）;dg[2][0]=0;

dg[2][1]=-6*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]）;

dg[2][2]=4*tanxingmoliang[i]*J_moliang[i]/changdu[i];dg[2][3]=0;

dg[2][4]=6*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]）;

dg[2][5]=2*tanxingmoliang[i]*J_moliang[i]/changdu[i];

dg[3][0]=-tanxingmoliang[i]*mianji[i]/changdu[i];dg[3][1]=0;dg[3][2]=0;

dg[3][3]=tanxingmoliang[i]*mianji[i]/changdu[i];dg[3][4]=0;dg[3][5]=0;dg[4][0]=0;

dg[4][1]=-12*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]*changdu[i]）;

dg[4][2]=6*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]）;dg[4][3]=0;

dg[4][4]=12*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]*changdu[i]）;

dg[4][5]=6*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]）;dg[5][0]=0;

dg[5][1]=-6*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]）;

dg[5][2]=2*tanxingmoliang[i]*J_moliang[i]/changdu[i];dg[5][3]=0;

dg[5][4]=6*tanxingmoliang[i]*J_moliang[i]/（changdu[i]*changdu[i]）;

dg[5][5]=4*tanxingmoliang[i]*J_moliang[i]/changdu[i];

}

//*************************

//求单元的转换矩阵及其逆阵

//*************************

voidjs_T_T1（inti,doubleT[6][6],doubleT1[6][6]）

{

intTti,Tti2;

for（Tti=0;Tti<6;Tti++）

{

//计算单元的转换矩阵T

for（Tti2=0;Tti2<6;Tti2++）

{

if（Tti==0&&Tti2==0）

T[Tti][Tti2]=cos（jiaodu[i]）;

elseif（Tti==0&&Tti2==1）

T[Tti][Tti2]=sin（jiaodu[i]）;

elseif（Tti==1&&Tti2==0）

T[Tti][Tti2]=-sin（jiaodu[i]）;

elseif（Tti==1&&Tti2==1）

T[Tti][Tti2]=cos（jiaodu[i]）;

elseif（Tti==2&&Tti2==2）

T[Tti][Tti2]=1;

elseif（Tti==3&&Tti2==3）

T[Tti][Tti2]=cos（jiaodu[i]）;

elseif（Tti==3&&Tti2==4）

T[Tti][Tti2]=sin（jiaodu[i]）;

elseif（Tti==4&&Tti2==3）

T[Tti][Tti2]=-sin（jiaodu[i]）;

elseif（Tti==4&&Tti2==4）

T[Tti][Tti2]=cos（jiaodu[i]）;

elseif（Tti==5&&Tti2==5）

T[Tti][Tti2]=1;

else

T[Tti][Tti2]=0;

if（fabs（T[Tti][Tti2]）<0.0000001）

T[Tti][Tti2]=0;

if（T[Tti][Tti2]>0.99999999）

T[Tti][Tti2]=1;

if（T[Tti][Tti2]<-0.99999999）

T[Tti][Tti2]=-1;

}

}

//计算T的转置矩阵T1

for（Tti=0;Tti<6;Tti++）

{

for（Tti2=0;Tti2<6;Tti2++）

T1[Tti2][Tti]=T[Tti][Tti2];

}

}

//*******************

//计算竿端内力的函数

//*******************

voidganduanneili（）

{

intTi,Ti2;

doublechengji=0;

doubledanyuan_gangdu[6][6];//单元刚度数组

for（inti=0;i

{

dy_gangdu（i,danyuan_gangdu）;

for（Ti=0;Ti<6;Ti++）

{

for（Ti2=0;Ti2<6;Ti2++）

chengji=chengji+danyuan_gangdu[Ti][Ti2]*dy_weiyi[i][Ti2];

gan_neili[i][Ti]=chengji+dy_neili[i][Ti];

chengji=0;

}

}

}

//************************************************************

//计算各个单元在局部坐标系中的位移的函数,存放在dy_weiyi[][]中

//************************************************************

voidjs_dy_weiyi（）

{

intTi,Ti3;

doubleT[6][6],T1[6][6];

doublechengji=0;

for（inti=0;i

{

js_T_T1（i,T,T1）;

intqianjiedian,houjiedian;

qianjiedian=dym_jdm[i][0];

houjiedian=dym_jdm[i][1];

for（Ti=0;Ti<6;Ti++）

{

dy_weiyi[i][Ti]=T[Ti][0]*weiyi[（qianjiedian-1）*3]+T[Ti][1]*weiyi[（qianjiedian-1）*3+1]+T[Ti][2]*weiyi[（qianjiedian-1）*3+2]+T[Ti][3]*weiyi[（houjiedian-1）*3]+T[Ti][4]*weiyi[（houjiedian-1）*3+1]+T[Ti][5]*weiyi[（houjiedian-1）*3+2];

}

}

}

//***********************

//计算得位移的函数

//***********************

voidjsweiyi（）

{

intTi,Ti2;

inti,j,k;

inti_max;

doublemax,zhongjian;

doublegd_kz[303][303];

for（Ti=0;Ti

{

for（Ti2=0;Ti2

gd_kz[Ti][Ti2]=zhengtigangdu[Ti][Ti2];

gd_kz[Ti][NJ*3]=F[Ti];

}

//用高斯主列消去法解方程,先计算增广矩阵的变形矩阵

for（Ti=0;Ti

{

max=fabs（gd_kz[Ti][Ti]）;

i_max=Ti;

for（Ti2=Ti;Ti2

{

if（fabs（gd_kz[Ti2][Ti]）>max）

{

i_max=Ti2-1;

max=gd_kz[Ti2][Ti];

}

}

if（i_max!

=Ti）

{

for（i=Ti;i

{

zhongjian=gd_kz[Ti][i];

gd_kz[Ti][i]=gd_kz[i_max][i];

gd_kz[i_max][i]=zhongjian;

}

}

for（j=Ti+1;j

{

doubleshang=gd_kz[j][Ti]/gd_kz[Ti][Ti];

for（k=0;k

gd_kz[j][k]=-shang*gd_kz[Ti][k]+gd_kz[j][k];

}

}

//反代求值

doublechengji=0;

weiyi[NJ*3-1]=gd_kz[NJ*3-1][NJ*3]/gd_kz[NJ*3-1][NJ*3-1];

for（i=NJ*3-2;i>=0;i--）

{

for（j=i+1;j

chengji=chengji+gd_kz[i][j]*weiyi[j];

weiyi[i]=（gd_kz[i][NJ*3]-chengji）/gd_kz[i][i];

chengji=0;

}

}

//********************

//引入支承条件的函数

//********************

voidzhichengyinru（）

{

inti,Ti;

intweizhi;

fangchengshu=NJ*3-NZ;

for（i=0;i

{

weizhi=zhichengweizhi[i]-1;

for（Ti=0;Ti

zhengtigangdu[weizhi][Ti]=0;

for（Ti=0;Ti

zhengtigangdu[Ti][weizhi]=0;

zhengtigangdu[weizhi][weizhi]=1;

F[weizhi]=0;

}

}

//*************************

//计算等效总节点载荷的函数

//*************************

voidzongzaihe（）

{

inti,Ti;

doubleU1=0,V1=0,M1=0,U2=0,V2=0,M2=0;

doublec,G,l,d;

doubleFf[6];

doubleT[6][6],T1[6][6];//T为坐标转换矩阵,T1为T的转置矩阵

for（i=0;i

F[i]=0;

for（i=0;i<6;i++）

Ff[i]=0;

for（i=0;i

{

for（intjl=0;jl<6;jl++）

dy_weiyi[i][jl]=0;

}

for（i=0;i

{

intdym=fjzhzuoyongdanyuan[i];

U1=V1=M1=U2=V2=M2=0;

for（intt=0;t<6;t++）

Ff[t]=0;

c=fjzhzuoyongdian[i];

l=changdu[（fjzhzuoyongdanyuan[i]-1）];

d=l-c;

G=fjzhzhi[i];

if（fjzhleixing[i]==1）

{

V1=-0.5*G*c*（2-2*c*c/（l*l）+c*c*c/（l*l*l））;

V2=-G*c-V1;

M1=G*c*c*（6-8*c/l+3*c*c/（l*l））/12;

M2=-G*c*c*c*（4-3*c/l）/（12*l）;

U1=U2=0;

}

elseif（fjzhleixing[i]==2）

{

V1=-G*（l+2*c）*d*d/（l*l*l）;

V2=-G*（l+2*d）*c*c/（l*l*l）;

M1=G*c*d*d/（l*l）;

M2=-G*c*c*d/（l*l）;

U1=U2=0;

}

elseif（fjzhleixing[i]==3）

{

U1=-G*d/l;

U2=-G*c/l;

M1=M2=V1=V2=0;

}

elseif（fjzhleixing[i]==4）

{

U1=U2=0;

V1=-6*G*c*d/（l*l*l）;

V2=-V1;

M1=G*d*（2-3*d/l）/l;

M2=G*c*（2-3*c/l）/l;

}

elseif（fjzhleixing[i]==5）

{

V1=-G*sin（fjzhjiaodu[i]）*（l+2*c）*d*d/（l*l*l）;

V2=-G*sin（fjzhjiaodu[i]）*（l+2*d）*c*c/（l*l*l）;

M1=G*sin（fjzhjiaodu[i]）*c*d*d/（l*l）;

M2=-G*sin（fjzhjiaodu[i]）*c*c*d/（l*l）;

U1=U2=0;

U1=-G*cos（fjzhjiaodu[i]）*d/l+U1;

U2=-G*cos（fjzhjiaodu[i]）*c/l+U2;

}

//记录竿的固端内力

dy_neili[dym-1][0]=U1;

dy_neili[dym-1][1]=V1;

dy_neili[dym-1][2]=M1;

dy_neili[dym-1][3]=U2;

dy_neili[dym-1][4]=V2;

dy_neili[dym-1][5]=M2;

js_T_T1（fjzhzuoyongdanyuan[i]-1,T,T1）;

for（Ti=0;Ti<6;Ti++）

Ff[Ti]=-（T1[Ti][0]*U1+T1[Ti][1]*V1+T1[Ti][2]*M1+T1[Ti][3]*U2+T1[Ti][4]*V2+T1[Ti][5]*M2）;

//将载荷转换到整体坐标系中

F[（dym_jdm[dym-1][0]-1）*3]=F[（dym_jdm[dym-1][0]-1）*3]+Ff[0];

F[（dym_jdm[dym-1][0]-1）*3+1]=F[（dym_jdm[d