有限元程序设计002.doc

有限元程序设计作业有限元程序设计作业程序功能说明本程序包含四个子程序，其主要功能分别是：void matmat void matvec 调用结点及单元信息判断并循环形成总刚double matvec double gs调用结点位移及结点信息和对应单元信息，回代求出并输出单元杆端力有限元程序流程图输入单元数 输入结点数根据单元数形成动态数组存储单元及结点信息以及单元截面信息输入非结点荷载数 有形成动态数组存储非结点荷载信息否根据非结点荷载信息判断并形成列向量输入结点荷载数否调用结点及单元信息判断并循环形成总刚形成总荷载列向量输入结点荷载信息并形成结点荷载列向量 有形成动态数组存储结构总刚信息引入支承条件，对总刚及荷载列向量进行处理输入铰结点数及其编号调用结点位移及结点信息和对应单元信息，回代求出并输出单元杆端力解线性方程组，求出结点位移铰结点的处理程序变量说明 NE-存储单元号，单元结点号，两端约束情况nelem-单元数nnode-结点数NM-存储单元弹模，面积，惯性矩，长度，与总体X轴夹角nf-单元自由度nf1-结点自由度 K-结构的总刚度矩阵K0-存储各单元单刚矩阵nelem1-非结点荷载单元数nnode1-有结点荷载结点数k-单元单刚矩阵nnode2-有约束的结点数nnode3-铰结点数 T-单元坐标转换矩阵 Po-非结点荷载列向列量F-单元非结点荷载的固端力nbg-起始结点号nend-终止结点号 NL-单元的非结点荷载信息 LD-单元非结点荷载的大小 Pa-所有单元结点荷载的列向量 LD1-存储各单元荷载信息 P-所有荷载之和列向量NR-存储支座约束的结点号F1-存储各单元的杆端力NLD-存储有结点荷载的结点号 JJ-存储铰结的结点编号 V-返回并存储结点位移列向量 V1-调用各单元的结点位移 源程序/ #include stdafx.h#include #include #define PI 3.14159265double mat366,X6;void matmat(int row,int col,double B166,double B266); void matvec(int row, int col,double mat166,double Y16); double matvec(int row, int col,double * mat1,double * Y1,double * V); double gs(int row,int col,double * A, double * B,double * V); int main(int argc, char* argv)double m=1e8;int nelem,nnode,nelem1,nnode1,nnode2,nnode3,nf,nf1;int i1,i2,i3,nbg,nend,nbg1,nend1,nbg2,nend2,nbg3,nend3,N,N1,N2,n,n1;double L,a,E,A,I,q,T66,T166,k66,mat166;double P16,V16,F16,X16;coutnelem;coutnnode;int * *NE;NE=new int * 5;for(i1=0;i1nelem;i1+)NEi1=new int 5;cout输入:单元编号，起节点号，止节点号，起点是否刚结(是0,否1)，终点是否刚结(是0,否1);for(i1=0;i1nelem;i1+)for(i2=0;i2NEi1i2;double * *NM;NM=new double * 5;for(i1=0;i1nelem;i1+)NMi1=new double 5;cout依次输入单元的:弹模，面积，惯性矩，长度，与x轴夹角;for(i1=0;i1nelem;i1+)for(i2=0;i2NMi1i2;nf1=3*nnode; double * K;K=new double *nf1;for(i1=0;i1nf1;i1+)Ki1=new double nf1;for(i1=0;i1nf1;i1+)for(i2=0;i2nf1;i2+)Ki1i2=0.0;nf=6*nelem; double * K0;K0=new double * nf;for(i1=0;i1nf;i1+)K0i1=new double nf;for(i1=0;i1nelem;i1+) N=NEi10,N1=6*(N-1),nbg=NEi11,nend=NEi12;nbg1=NEi13,nend1=NEi14;E=NMi10,A=NMi11,I=NMi12;L=NMi13,a=NMi14*PI/180.0;if(nbg1=0)&(nend1=0)k01=k02=k04=k05=k10=k13=k20=k23=0; k31=k32=k34=k35=k40=k43=k50=k53=0;k00=k33=(E*A)/L,k03=k30=-(E*A)/L;k11=k44=(12*E*I)/pow(L,3);k14=k41=-(12*E*I)/pow(L,3);k12=k21=k15=k51=(6*E*I)/pow(L,2);k24=k42=k45=k54=-(6*E*I)/pow(L,2);k22=k55=4*E*I/L,k25=k52=2*E*I/L;else if(nbg1=0)&(nend1=1)k01=k02=k04=k05=0;k10= k13=k15=k20=k23=k25=0;k31=k32=k34=k35=0;k40=k43=k45=0;k50=k51=k52=k53=k54=k55=0;k00=k33=(E*A)/L,k03=k30=-(E*A)/L;k11=k44=3*E*I/pow(L,3);k14=k41=-(3*E*I)/pow(L,3);k12=k21=3*E*I/pow(L,2),k24=k42=-3*E*I/pow(L,2);k22=3*E*I/L;else if(nbg1=1)&(nend1=1)for(i2=0;i26;i2+)for(i3=0;i36;i3+)ki2i3=0.0;k00=k33=E*A/L,k03=k30=-E*A/L;elsek01=k02=k04=k05=k10=k12=k13=0;k20=k21=k22=k23=k24=k25=0;k31=k32=k34=k35=0;k40=k42=k43=k50=k52=k53=0;k00=k33=(E*A)/L,k03=k30=-(E*A)/L;k11=k44=3*E*I/pow(L,3);k14=k41=-(3*E*I)/pow(L,3); k15=k51=3*E*I/pow(L,2);k45=k54=-3*E*I/pow(L,2),k55=3*E*I/L;for(i2=0;i26;i2+)N2=N1+i2;for(i3=0;i36;i3+)K0N2i3=ki2i3;for(i2=0;i26;i2+)for(i3=0;i36;i3+)Ti2i3=0.0;T22=T55=1,T01=T34=-sin(a),T10=T43=sin(a);T00=T11=T33=T44=cos(a);matmat(6,6,T,k);for(i2=0;i26;i2+)for(i3=0;i36;i3+)T1i2i3=Ti3i2;for(i2=0;i26;i2+)for(i3=0;i36;i3+)Ti2i3=T1i2i3;for(i2=0;i26;i2+) for(i3=0;i36;i3+)mat1i2i3=mat3i2i3;matmat(6,6,mat1,T);nbg3=3*(nbg-1),nend3=3*(nend-1);for(i2=0;i23;i2+) nbg2=nbg3+i2; for(i3=0;i33;i3+) nend2=nbg3+i3; Knbg2nend2+=mat3i2i3; for(i2=0;i23;i2+) nbg2=nbg3+i2;for(i3=3;i36;i3+) nend2=nend3+i3-3; Knbg2nend2+=mat3i2i3; for(i2=3;i26;i2+)nbg2=nend3+i2-3;for(i3=0;i33;i3+) nend2=nbg3+i3; Knbg2nend2+=mat3i2i3; for(i2=3;i26;i2+)nbg2=nend3+i2-3;for(i3=3;i36;i3+)nend2=nend3+i3-3; Knbg2nend2+=mat3i2i3;double * Po; Po=new double nf1;for(i1=0;i1nf1;i1+)Poi1=0.0; double * * F; F=new double *6;for(i1=0;i1nelem;i1+)Fi1=new double 6;for(i1=0;i1nelem;i1+)for(i2=0;i26;i2+)Fi1i2=0.0; coutnelem1;if(nelem1!=0)int * * NL; NL=new int *3; for(i1=0;i1nelem1;i1+) NLi1=new int3;cout依次输入非结点荷载信息:单元编号，非结点荷载类（均布0，三角形荷载1）,荷载在单元上的布设情况(均布0，三角形荷载:起点为零取1,起点为q取2); for(i1=0;i1nelem1;i1+) or(i2=0;i2NLi1i2; double * LD; LD=new double nelem1; cout依次输入各单元的均布荷载值或三角形荷载的最大值:; for(i1=0;i1LDi1; for(i1=0;i1nelem1;i1+)N=NLi10-1,n=NLi11,n1=NLi12;nbg=NEN1,nend=NEN2,nbg1=NEN3;nend1=NEN4;E=NMN0,A=NMN1,I=NMN2;L=NMN3,a=NMN4*PI/180.0,q=LDi1; for(i2=0;i26;i2+) for(i3=0;i36;i3+) Ti2i3=0.0; T22=T55=1,T01=T34=-sin(a),T10=T43=sin(a);T00=T11=T33=T44=cos(a);if(n=0)if(nbg1=0)&(nend1=0)P10=P13=0,P11=P14=q*L/2;P12=q*pow(L,2)/12,P15=-q*pow(L,2)/12;else if(nbg1=0)&(nend1=1) P10=P13=P15=0,P11=5*q*L/8;P14=3*q*L/8,P12=q*pow(L,2)/8;else if (nbg1=1)&(nend1=0)P10=P12=P13=0,P11=3*q*L/8;P14=5*q*L/8,P15=-q*pow(L,2)/8;elseP10=P12=P13=P15=0,P11=q*L/2,P14=q*L/2;else if(n1=1) if(nbg1=0)&(nend1=0) P10=P13=0,P11=3*q*L/20,P14=7*q*L/20; P12=q*pow(L,2)/30,P15=-q*pow(L,2)/20; else if(nbg1=0)&(nend1=1) P10=P13=P15=0,P11=9*q*L/40; P14=11*q*L/40,P12=7*q*pow(L,2)/120; else if (nbg1=1)&(nend1=0) P10=P12=P13=0,P11=q*L/10; P14=2*q*L/5,P15=-q*pow(L,2)/15; else P10=P12=P13=P15=0,P11=q*L/6,P14=q*L/3;elseif(nbg1=0)&(nend1=0) P10=P13=0,P11=7*q*L/20,P14=3*q*L/20; P12=q*pow(L,2)/20,P15=-q*pow(L,2)/30; else if(nbg1=0)&(nend1=1) P10=P13=P15=0,P11=2*q*L/5; P14=q*L/10,P12=q*pow(L,2)/15; else if (nbg1=1)&(nend1=0) P10=P12=P13=0,P11=11*q*L/40; P14=9*q*L/40,P15=-7*q*pow(L,2)/120; else P10=P12=P13=P15=0,P11=q*L/3,P14=q*L/6;for(i2=0;i26;i2+)FNi2=P1i2;matvec(6,6,T,P1);for(i2=0;i26;i2+)P1i2=-Xi2;nbg3=3*(nbg-1),nend3=3*(nend-1);for(i2=0;i23;i2+)nbg2=nbg3+i2;Ponbg2+=P1i2;for(i2=3;i26;i2+)nbg2=nend3+i2-3;Ponbg2+=P1i2;double * Pa; Pa=new double nf1;for(i1=0;i1nf1;i1+)Pai1=0.0;coutnnode1;if(nnode1!=0)int * NLD;NLD=new int nnode1;cout输入有结点荷载的结点序号:;for(i1=0;i1NLDi1;double * LD1; LD1=new double *3;for(i1=0;i1nnode1;i1+)LD1i1=new double 3;cout按以上结点号顺序依次输入各结点的结点荷载的:x轴方向力，y轴方向力，弯矩;for(i1=0;i1nnode1;i1+)for(i2=0;i2LD1i1i2;for(i1=0;i1nnode1;i1+)N=NLDi1-1,nbg3=3*N;for(i2=0;i23;i2+) nbg2=nbg3+i2; Panbg2=LD1i1i2;delete NLD;for(i1=0;i1nnode1;i1+)delete LD1i1;delete LD1;double * P; P=new double nf1;for(i1=0;i1nf1;i1+)Pi1=Poi1+Pai1;delete Po;delete Pa;coutnnode2;int * NR;NR=new int nnode2;cout输入有支座约束的结点号:;for(i1=0;i1NRi1;for(i1=0;i1nnode2;i1+)N=NRi1-1;for(i2=0;i23;i2+)nbg=3*N+i2;Knbgnbg=m*Knbgnbg,Pnbg=0;coutnnode3;if(nnode3!=0)int * JJ;JJ=new int nnode3;cout输入铰结点的结点号:;for(i1=0;i1JJi1;for(i1=0;i1nnode3;i1+)nbg=3*JJi1-1;Knbgnbg=m,Pnbg=0;delete JJ;double * V; V=new double nf1;gs(nf1,nf1,K,P,V);delete NR;cout单元号 IN IQIM JN JQ JMendl;for(i1=0;i1nelem;i1+) N=NEi10,N1=6*(N-1),nbg=NEi11;nend=NEi12,a=NMi14*PI/180.0;nbg1=3*(nbg-1),nend1=3*(nend-1);for(i2=0;i26;i2+)for(i3=0;i36;i3+)Ti2i3=0.0;T22=T55=1,T01=T34=sin(a),T10=T43=-sin(a);T00=T11=T33=T44=cos(a);for(i2=0;i23;i2+) nbg2=nbg1+i2;V1i2=Vnbg2;for(i2=3;i26;i2+)nend2=nend1+i2-3;V1i2=Vnend2;matvec(6,6,T,V1);for(i2=0;i26;i2+)nbg2=N1+i2;for(i3=0;i36;i3+)ki2i3=K0nbg2i3;for(i2=0;i26;i2+)X1i2=Xi2;matvec(6,6,k,X1);for(i2=0;i23;i2+)nbg2=N-1;F1i2=Fnbg2i2;for(i2=3;i26;i2+)nbg2=N-1;F1i2=Fnbg2i2;for(i2=0;i26;i2+)F1i2=Xi2+F1i2;coutN ;for(i2=0;i26;i2+)coutF1i2 ;coutendl;return 0;void matmat(int row,int col,double B166,double B266) int i1,i2,i3;for(i1=0;i1row;i1+)for(i2=0;i2col;i2+)mat3i1i2=0.0;for(i3=0;i3col;i3+)mat3i1i2+=B1i1i3*B2i3i2; void matvec(int row, int col,double mat166,double Y16) int i1,i2;for(i1=0;i1row;i1+)Xi1=0.0;for(i2=0;i2col;i2+)Xi1+=mat1i1i2*Y1i2; double matvec(int row, int col,double * mat1,double * Y1,double * V) int i1,i2;for(i1=0;i1row;i1+)Vi1=0.0;for(i2=0;i2col;i2+)Vi1+=mat1i1i2*Y1i2; return * V;double gs(int row,int col,double * A, double * B,double * V) int i1,i2,i3,irow;double db,db1;double * k1;k1=new double col;for(i1=0;i1row;i1+) for(i2=i1;i2col;i2+)k1i2=Ai1i2;irow=i1;db=fabs(Ai1i1);for(i2=i1+1;i2row;i2+)if(dbfabs(Ai2i1)db=fabs(Ai2i1);irow=i2;if(db=0.0)cout对角线元素为零,方程组无解.;return false;for(i2=i1;i2col;i2+)Ai1i2=Airowi2;Airowi2=k1i2;db=Ai1i1; db1=Bi1;Bi1=Birow;Birow=db1;for(i2=i1;i2col;i2+) Ai1i2/=db;Bi1/=db;for(i2=i1+1;i2row;i2+)db=Ai2i1;for(i3=i1;i3=0;i1-)for(i2=col-1;i2i1;i2-)Bi1-=Bi2*Ai