有限元教材PFF程序使用说明(C++版).pdf

《有限元教材PFF程序使用说明(C++版).pdf》由会员分享，可在线阅读，更多相关《有限元教材PFF程序使用说明(C++版).pdf（20页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、1.8 平面刚架程序设计平面刚架程序设计 1.8.1 程序说明及总流程图程序说明及总流程图 依据前述平面刚架矩阵位移法原理，本节使用 Visual C+6.0 编译器编制了平面刚架先处理法计算程序 PFF，该程序未使用 C+面向对象的编程方法，而仍沿用兼容 C 的结构化程序设计方法，以便和本科生先修课程 C 语言相容。本程序不建议使用 TC+3.0 及更低版本的编译器编译。一一PFF 的主要功能和特点的主要功能和特点(1)输入单元编号、结点编号和结点位移分量统一编码，单元材料信息；(2)先处理法形成结构刚度矩阵和综合结点荷载列阵；(3)Gauss 消元法解线性代数方程组；(4)计算并输出结点位

2、移和单元杆端力；二二PFF 的使用方法的使用方法(1)输入输出数据文件 PFF 从原始数据文件中读取结构的离散化信息，经计算后，将结果输出到结果数据文件中。运行编译好的 PFF 可执行程序，会提示“Please input primary data file name!”，这时输入事先准备好的原始数据文件名，回车后，提示“Please input calculation result file name!”，输入一个文件名，PFF 会自动以该文件名生成结果数据文件。原始数据文件最好同 PFF 可执行程序放置在同一文件夹中，结果数据文件不应是已有的文件，而必须由 PFF 自动生成。两文件建议使用

3、英文和数字命名。在输入文件名时，应包含后缀。C 语言以 0 为数组各维的开始下标，称为 0 基数组；而一般习惯是以 1 为数组各维的开始下标，称为 1 基数组。PFF 源程序编制时，仍以 C 语言为标准，使用 0 基数组，即 PFF程序内部处理机制是 0 基的。但是为方便读者，原始数据文件和输出数据文件则按一般习惯，使用 1 基数组对结点、单元、结点位移分量、结点荷载和非结点荷载进行编号，即 PFF 程序外部表现是 1 基的。(2)原始数据文件的格式 下面依行序说明原始数据文件输入数据的含义：总体信息 ne,nj,n,np,nf,e ne单元总数；nj结点总数；n结点位移未知量总数；np结点荷

4、载总数；nf非结点荷载总数；e弹性模量 结点坐标(xi,yi)(i=0;inj;i+)xi结点 i+1 的 x 坐标；yi结点 i+1 的 y 坐标。依结点编号顺序，先 x 坐标，后 y 坐标，存储所有结点坐标。每个结点坐标占一行。单元信息(iji0,iji1,ai,zii)(i=0;ine;i+)iji0单元 i+1 的始结点码；iji1单元 i+1 的末结点码；ai单元 i+1 的截面积；zii单元 i+1 的截面惯性矩。依单元编号顺序，存储各单元的始末结点及截面积 A 和惯性矩 I。结点位移分量统一编码(jnij(j=0;j3;j+)(i=0;inj;i+)jnij结点 i+1 的结点位

5、移分量统一编码。依结点编号的顺序，存储各结点的结点位移分量统一编码。每结点位移分量均按 u、v、的顺序依次输入，被约束的位移分量输 0，未被约束的位移分量则顺序编号输入。结点荷载(pjij(j=0;j3;j+)(i=0;inp;i+)pjij第 i+1 号结点荷载；其中 pji0第 i+1 号结点荷载作用的结点编号。pji1第 i+1 号结点荷载作用的方向代码，整体坐标系 x、y 和转角方向分别用 1.0，2.0 和 3.0 表示。pji2第 i+1 号结点荷载的大小，带符号，与整体坐标系一致（力偶顺时针）为正，反之为负。依上面的顺序，存储全部结点荷载。非结点荷载(pfij(j=0;j4;j+

6、)(i=0;i0 i=0,np-1 l!=0 pl=pji2 nf0 i=0,nf-1 m=int(pfi0)调用 lsc 调用 eff 调用 elv 按式（1-24）求EeF j=0,5 l=lvj l!=0 pl=pl+pej 否 是 第一步：处理结点荷载 否（无非结点荷载）是 否（无结点荷载）是 否 是 11 返回对结点荷载循环 荷载作用结点号 对非结点荷载循环 荷载作用单元号 求单元 m 的常数 荷载方向代码 荷载定位编码 结点荷载送入 F 求PeF 形成单元定位向量 EeF存入数组 pe6 对eEF的分量循环 EeF的第 j 个分量在F中的地址 排除定位向量的零分量 将eEF叠加到F

7、 第二步：处理非结点荷载 图 1.19 子程序 jlp 的框图 m=0,ne-1 调用 lsc 调用 esm 调用 elv i=0,5 l=lvi d(i)=0.0 l!=0 是 否 di=pl eee FKeeFTFnf0 i=0,nf-1 l=int(pfi0)m=l 11调用 eff 得 fo6 j=0,5 fj=fj+foj 2否（无非结点荷载）是 否 是 形成单元定位向量 对e 六个分量循环 若l0，则l是e 中第 i 个位移分量 取单元杆端位移 存在数组 fd6中 存在数组 f6中 对单元循环 叠加单元固端力求单元固端力 若 m=l，则荷载作用在单元m上荷载作用单元号eF求Ke 2

8、对单元循环 求单元常数 eF第一步：计算由结点位移引起的杆端力 第二步：有非结点荷载的单元，还需叠加单元固端力 图 1.20 子程序 mvn 的框图（计算杆端力部分）1.8.3 算例及源程序算例及源程序 一算例一算例 下面以一个例子，说明原始数据文件的输入方法，以及结果数据文件的含义。例例 1.4 已知如图 1.21a 所示的刚架，各杆，。为该刚架准备 PFF 程序使用的原始数据文件。72=3 10 kN/mE20.16mA 40.002mI 3m3m6m6m20kN m15kN10kN5kN/m1xy31(1,2,3)2(0,0,0)3(4,5,6)4(4,5,7)5(8,9,10)6(0,

9、0,11)7(12,0,13)245(a)(b)结构离散化 图 1.21 解解 按本章 1.2.1 的方法进行离散化，对单元、结点和结点位移未知量进行编号，如图1.21b 所示。接下来建立一个文本文件，按上述原始数据文件的格式，输入原始数据如下（右侧为注释，不输入）ne,nj,n,np,nf,e 5,7,13,2,2,3.0e7(xi,yi)(i=0;inj;i+)0,0 6,0 0,6 0,6 6,6 0,12 6,12(iji0,iji1,ai,zii)(i=0;ine;i+)1,2,0.16,0.002 1,3,0.16,0.002 4,5,0.16,0.002 3,6,0.16,0.0

10、02 5,7,0.16,0.002(jnij(j=0;j3;j+)(i=0;inj;i+)1,2,3 0,0,0 4,5,6 4,5,7 8,9,10 0,0,11 12,0,13(pjij(j=0;j3;j+)(i=0;inp;i+)1,3.0,-20.0 3,1.0,10.0(pfij(j=0;j4;j+)(i=0;inf;i+)1,2,15.0,3.0 5,1,5.0,6.0 结果数据文件的含义结果数据文件的含义 仍以例 1.4 为例，PFF 输出的结果数据文件如下（括号中的汉字部分为说明）：PLANE FRAME STRUCTURE ANALYSIS（总体信息）NE=5,NJ=7,N=

11、13,NP=2,NF=2,E=3e+007 COORDINATES OF JOINT（结点坐标信息）Joint x y （结点号）（x坐标）（y坐标）1 0 0 2 6 0 3 0 6 4 0 6 5 6 6 6 0 12 7 6 12 INFORMATION OF ELEMENTS（单元信息）Element Joint-i Joint-j A ZI （单元号）（始结点号）（末结点号）（横截面积）（截面惯性矩）1 1 2 0.16 0.002 2 1 3 0.16 0.002 3 4 5 0.16 0.002 4 3 6 0.16 0.002 5 5 7 0.16 0.002 INFORMAT

12、ION OF JOINT DEGREES OF FREEDOM（结点位移分量统一编码）Joint u v ceta （结点号）(u方向)（v方向）（方向）1 1 2 3 2 0 0 0 3 4 5 6 4 4 5 7 5 8 9 10 6 0 0 11 7 12 0 13 JOINT LOAD（结点荷载）Joint XYM Load （结点号）（方向代码）（荷载大小）1 3 -20 3 1 10 NON-JOINT LOAD（非结点荷载）Element Type Load c （单元号）（荷载类型）（荷载大小）（荷载位置参数c）1 2 15 3 5 1 5 6（以上六部分直接来自原始数据文件，

13、可作校核用）（以下两部分为计算结果）JOINT DISPLACEMENTS（结点位移）Joint u v ceta （节点号）(u方向)（v方向）（方向）1 -1.6151e-005 -1.6406e-005 6.6171e-004 2 0.0000e+000 0.0000e+000 0.0000e+000 3 -6.7499e-003 -1.7578e-005 2.9077e-004 4 -6.7499e-003 -1.7578e-005 -1.4939e-003 5 -6.7874e-003 1.8750e-005 3.0061e-003 6 0.0000e+000 0.0000e+000

14、 -1.8329e-003 7 -3.8324e-002 0.0000e+000 6.0061e-003 MEMBER-END FORCES OF ELEMENTS（单元杆端内力）（注意正符号规定不同于传统位移法，轴力剪力与单元坐标系一致时为正）（i 代表单元始端，j 代表单元末端）Element FN FQ M（单元号）（杆端轴力）（杆端剪力）（杆端弯矩）1(i)-12.9212 -0.937619 15.0542 (j)12.9212 -14.0624 24.3201 2(i)0.937619 -12.9212 -35.0542 (j)-0.937619 12.9212 -42.4729

15、3(i)30 15 0 (j)-30 -15 90 4(i)-14.0624 7.07882 42.4729 (j)14.0624 -7.07882 0 5(i)15 -30 -90 (j)-15 0 0 根据结果数据文件，可绘出结构的内力图，如图 1.22 所示。15.0524.3235.0542.4790.0090.0022.5022.500.9414.0612.9215.007.0830.0012.920.9430.0014.0615.00图(kN)图(kNm)MFQ(a)(b)图(kN)(c)NF 图 1.22 例 1.4 所示结构内力图 二二PFF 的源程序的源程序 PFF 源程序中

16、，凡是“/*”和“*/”之间，以及同一行内“/”之后的内容，均为注释内容，不会被编译、执行。/*STATIC ANALYSIS FOR PLANE FRAME*/*(METHOD OF FORMER TREATMENT)*/#include stdio.h#include math.h#include stdlib.h FILE*fpin,*fpout;const int ARRSIZE=50;/Array size const int JDOF=3;/Joint degree of freedom const int JPE=2;/Joints per element /*Main prog

17、ram reads the control parameters and organizes*/*the whole calculation by calling subroutines.*/main()char indat30,outdat30;int ne,nj,n,np,nf;double e;double xARRSIZE,yARRSIZE,aARRSIZE,ziARRSIZE,pjARRSIZE3,pfARRSIZE4;double pJDOF*ARRSIZE,tkJDOF*ARRSIZEJDOF*ARRSIZE;int ijARRSIZEJPE,jnARRSIZEJDOF;int

18、input(int ne,int nj,int np,int nf,double*x,double*y,int(*ij)JPE,double*a,double*zi,int(*jn)JDOF,double(*pj)3,double(*pf)4);int tsm(int ne,int nj,int n,double e,double*x,double*y,int(*ij)JPE,double*a,double*zi,int(*jn)JDOF,double(*tk)JDOF*ARRSIZE);int jlp(int ne,int nj,int n,int np,int nf,double*x,do

19、uble*y,int(*ij)JPE,int(*jn)JDOF,double(*pj)3,double(*pf)4,double*p);int gauss(double(*a)JDOF*ARRSIZE,double*b,int n);int mvn(int ne,int nj,int n,int nf,double e,double*x,double*y,int(*ij)JPE,double*a,double*zi,int(*jn)JDOF,double(*pf)4,double*p);printf(Please input primary data file name!n);scanf(%s

20、,indat);if(fpin=fopen(indat,r)=NULL)printf(Cannot open input file.n);return(1);printf(Please input calculation result file name!n);scanf(%s,outdat);if(fpout=fopen(outdat,w)=NULL)printf(Cannot open result file.n);return(1);fscanf(fpin,%d,%d,%d,%d,%d,%lg,&ne,&nj,&n,&np,&nf,&e);fprintf(fpout,PLANE FRAM

21、E STRUCTURE ANALYSISnn);fprintf(fpout,NE=%d NJ=%d N=%d NP=%d NF=%d E=%lgn,ne,nj,n,np,nf,e);input(ne,nj,np,nf,x,y,ij,a,zi,jn,pj,pf);tsm(ne,nj,n,e,x,y,ij,a,zi,jn,tk);jlp(ne,nj,n,np,nf,x,y,ij,jn,pj,pf,p);gauss(tk,p,n);mvn(ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p);fclose(fpin);fclose(fpout);printf(nCalculation

22、finished!n);printf(Please check result file%s for more information.n,outdat);return(0);/*Read and print primary data.*/int input(int ne,int nj,int np,int nf,double*x,double*y,int(*ij)JPE,double*a,double*zi,int(*jn)JDOF,double(*pj)3,double(*pf)4)int i,j;for(i=0;inj;i+)fscanf(fpin,%lg,%lg,&xi,&yi);fpr

23、intf(fpout,nCOORDINATES OF JOINTSn);fprintf(fpout,Joint x yn);for(i=0;inj;i+)fprintf(fpout,%-4d%-7lg%-7lgn,i+1,xi,yi);/Base-0 adjustment.for(i=0;ine;i+)fscanf(fpin,%d,%d,%lg,%lg,&iji0,&iji1,&ai,&zii);fprintf(fpout,nINFORMATION OF ELEMENTSn);fprintf(fpout,Element Joint-i Joint-j A ZIn);for(i=0;ine;i+

24、)fprintf(fpout,%-4d%-4d%-4d%-7lg%-7lgn,i+1,iji0,iji1,ai,zii);/Base-0 adjustment.for(i=0;inj;i+)for(j=0;jJDOF;j+)fscanf(fpin,%d,&jnij);fprintf(fpout,nINFORMATION OF JOINT DEGREES OF FREEDOMn);fprintf(fpout,Joint u v cetan);for(i=0;inj;i+)fprintf(fpout,%-4d,i+1);/Base-0 adjustment.for(j=0;j0)for(i=0;i

25、np;i+)for(j=0;j3;j+)fscanf(fpin,%lg,&pjij);fprintf(fpout,nJOINT LOADn);fprintf(fpout,Joint XYM Loadn);for(i=0;inp;i+)for(j=0;j0)for(i=0;inf;i+)for(j=0;j4;j+)fscanf(fpin,%lg,&pfij);fprintf(fpout,nNON-JOINT LOADn);fprintf(fpout,Element Type Load cn);for(i=0;inf;i+)for(j=0;j4;j+)fprintf(fpout,%-7lg,pfi

26、j);fprintf(fpout,n);return(0);/*Calculate length,sine and cosine of a member.*/int lsc(int m,int ne,int nj,double*x,double*y,int(*ij)JPE,double*bl,double*si,double*co)int i,j;double dx,dy;i=ijm0-1,j=ijm1-1;/Base-0 adjustment.dx=xj-xi;dy=yj-yi;*bl=sqrt(dx*dx+dy*dy);*si=dy/*bl;*co=dx/*bl;return(0);/*S

27、et up element location vector.*/int elv(int m,int ne,int nj,int(*ij)JPE,int(*jn)JDOF,int*lv)int i,j,k;i=ijm0-1,j=ijm1-1;/Base-0 adjustment.for(k=0;kJDOF;k+)lvk=jnik;lvk+JDOF=jnjk;return(0);/*Calculate element stiffness matrix referred*/*to global coordinate system.*/int esm(int m,int ne,double e,dou

28、ble*a,double*zi,double bl,double si,double co,double ekJPE*JDOFJPE*JDOF)int i,j;double c1,c2,c3,c4,s1,s2,s3,s4,s5,s6;c1=e*am/bl;c2=2.0*e*zim/bl;c3=3.0*c2/bl;c4=2.0*c3/bl;s1=c1*co*co+c4*si*si;s2=(c1-c4)*si*co;s3=c3*si;s4=c1*si*si+c4*co*co;s5=c3*co;s6=c2;ek00=s1;ek01=s2;ek02=-s3;ek03=-s1;ek04=-s2;ek05

29、=-s3;ek11=s4;ek12=s5;ek13=-s2;ek14=-s4;ek15=s5;ek22=2.0*s6;ek23=s3;ek24=-s5;ek25=s6;ek33=s1;ek34=s2;ek35=s3;ek44=s4;ek45=-s5;ek55=2.0*s6;for(i=0;iJPE*JDOF-1;i+)for(j=i+1;jJPE*JDOF;j+)ekji=ekij;return(0);/*Assemble total stiffness matrix.*/int tsm(int ne,int nj,int n,double e,double*x,double*y,int(*i

30、j)JPE,double*a,double*zi,int(*jn)JDOF,double(*tk)JDOF*ARRSIZE)double ekJPE*JDOFJPE*JDOF;int lvJPE*JDOF;int i,j,l,k,m;double bl,si,co;int lsc(int m,int ne,int nj,double*x,double*y,int(*ij)JPE,double*bl,double*si,double*co);int esm(int m,int ne,double e,double*a,double*zi,double bl,double si,double co

31、,double ekJPE*JDOFJPE*JDOF);int elv(int m,int ne,int nj,int(*ij)JPE,int(*jn)JDOF,int*lv);for(i=0;in;i+)for(j=0;jn;j+)tkij=0.0;for(m=0;mne;m+)lsc(m,ne,nj,x,y,ij,&bl,&si,&co);esm(m,ne,e,a,zi,bl,si,co,ek);elv(m,ne,nj,ij,jn,lv);for(l=0;lJPE*JDOF;l+)i=lvl;if(i!=0)for(k=0;kJPE*JDOF;k+)j=lvk;if(j!=0)tki-1j

32、-1=tki-1j-1+eklk;/Base-0 adjustment.return(0);/*Calculate element fixed-end forces.*/int eff(int i,double(*pf)4,int nf,double bl,double*fo)int no;double q,c,b,c1,c2,c3;int j;no=int(pfi1);q=pfi2;c=pfi3;b=bl-c;c1=c/bl;c2=c1*c1;c3=c1*c2;for(j=0;jJPE*3;j+)foj=0.0;switch(no)case 1:fo1=-q*c*(1.0-c2+c3/2.0

33、);fo2=-q*c*c*(0.5-2.0*c1/3.0+0.25*c2);fo4=-q*c*c2*(1.0-0.5*c1);fo5=q*c*c*c1*(1.0/3.0-0.25*c1);return(0);case 2:fo1=-q*b*b*(1.0+2.0*c1)/bl/bl;fo2=-q*c*b*b/bl/bl;fo4=-q*c2*(1.0+2.0*b/bl);fo5=q*c2*b;return(0);case 3:fo1=6.0*q*c1*b/bl/bl;fo2=q*b*(2.0-3.0*b/bl)/bl;fo4=-6.0*q*c1*b/bl/bl;fo5=q*c1*(2.0-3.0*

34、c1);return(0);case 4:fo1=-q*c*(0.5-0.75*c2+0.4*c3);fo2=-q*c*c*(1.0/3.0-0.5*c1+0.2*c2);fo4=-q*c*c2*(0.75-0.4*c1);fo5=q*c*c*c1*(0.25-0.2*c1);return(0);case 5:fo0=-q*c*(1.0-0.5*c1);fo3=-0.5*q*c*c1;return(0);case 6:fo0=-q*b/bl;fo3=-q*c1;return(0);return(0);/*Form total joint load vector.*/int jlp(int ne

35、,int nj,int n,int np,int nf,double*x,double*y,int(*ij)JPE,int(*jn)JDOF,double(*pj)3,double(*pf)4,double*p)int i,j,k,l,m;int ii;double bl,si,co;int lvJPE*JDOF;double fo6,peJPE*JDOF;int lsc(int m,int ne,int nj,double*x,double*y,int(*ij)JPE,double*bl,double*si,double*co);int elv(int m,int ne,int nj,int

36、(*ij)JPE,int(*jn)JDOF,int*lv);int eff(int i,double(*pf)4,int nf,double bl,double*fo);for(i=0;i0)for(i=0;inp;i+)j=int(pji0);k=int(pji1);l=jnj-1k-1;/Base-0 adjustment.if(l!=0)pl-1=pji2;/Base-0 adjustment.fprintf(fpout,nPJ=n);for(ii=0;ii0)for(i=0;inf;i+)m=int(pfi0)-1;/Base-0 adjustment.lsc(m,ne,nj,x,y,

37、ij,&bl,&si,&co);eff(i,pf,nf,bl,fo);elv(m,ne,nj,ij,jn,lv);pe0=-fo0*co+fo1*si;pe1=-fo0*si-fo1*co;pe2=-fo2;pe3=-fo3*co+fo4*si;pe4=-fo3*si-fo4*co;pe5=-fo5;for(j=0;jJPE*JDOF;j+)l=lvj;if(l!=0)pl-1=pl-1+pej;/Base-0 adjustment.fprintf(fpout,nP=n);for(ii=0;iin;ii+)fprintf(fpout,%12.4lg,pii);fprintf(fpout,n);

38、return(0);/*Solution of simultaneous equations by the GAUSS*/*elimination method.*/int gauss(double(*a)JDOF*ARRSIZE,double*b,int n)int i,j,k;double a1;for(k=0;kn-1;k+)for(i=k+1;in;i+)a1=aki/akk;for(j=k+1;j-1;i-)/Base-0 adjustment.for(j=i+1;jn;j+)bi=bi-aij*bj;bi=bi/aii;return(0);/*Print joint displac

39、ements.Calculate and print*/*member-end forced of elements.*/int mvn(int ne,int nj,int n,int nf,double e,double*x,double*y,int(*ij)JPE,double*a,double*zi,int(*jn)JDOF,double(*pf)4,double*p)int i,j,l,m,ii;int lvJPE*JDOF;double bl,si,co;double fo6,dJPE*JDOF,fdJPE*JDOF,f6;double ekJPE*JDOFJPE*JDOF;int

40、lsc(int m,int ne,int nj,double*x,double*y,int(*ij)JPE,double*bl,double*si,double*co);int esm(int m,int ne,double e,double*a,double*zi,double bl,double si,double co,double ekJPE*JDOFJPE*JDOF);int elv(int m,int ne,int nj,int(*ij)JPE,int(*jn)JDOF,int*lv);int eff(int i,double(*pf)4,int nf,double bl,doub

41、le*fo);fprintf(fpout,nJOINT DISPLACEMENTSn);fprintf(fpout,Joint u v cetan);for(j=0;jnj;j+)for(i=0;iJDOF;i+)di=0.0;l=jnji;if(l!=0)di=pl-1;/Base-0 adjustment.fprintf(fpout,%-4d,j+1);/Base-0 adjustment.for(ii=0;iiJDOF;ii+)fprintf(fpout,%12.4le,dii);fprintf(fpout,n);fprintf(fpout,nMEMBER-END FORCES OF E

42、LEMENTSn);fprintf(fpout,Element FN FQ Mn);for(m=0;mne;m+)lsc(m,ne,nj,x,y,ij,&bl,&si,&co);esm(m,ne,e,a,zi,bl,si,co,ek);elv(m,ne,nj,ij,jn,lv);for(i=0;iJPE*JDOF;i+)l=lvi;di=0.0;if(l!=0)di=pl-1;/Base-0 adjustment.for(i=0;iJPE*JDOF;i+)fdi=0.0;for(j=0;j0)for(i=0;inf;i+)l=int(pfi0);if(m=l-1)/Base-0 adjustment.eff(i,pf,nf,bl,fo);for(j=0;jJPE*3;j+)fj=fj+foj;fprintf(fpout,%4d(i),m+1);/Base-0 adjustment.for(i=0;i3;i+)fprintf(fpout,%12.6lg,(fabs(fi)1.0e-12)?0.0:fi);fprintf(fpout,n);fprintf(fpout,(j);for(i=0;iJDOF;i+)fprintf(fpout,%12.6lg,(fabs(fi+3)1.0e-12)?0.0:fi+3);fprintf(fpout,n);return(0);