有限元源程序

1、/*程序开始*#include#include #define NE 6#define NJ 6#define NZ 6#define NPJ 0#define NPF 1#define NJ3 18#define DD 12#define E0 3.0000E7#define A0 0.5#define I0 4.16667E-2#define PI 3.141592654FILE *fp;这是输入参数的初始化和定义全局变量*intjmNE+13=0,0,0,0,1,3,032,035,0,6,5,0,5,4,024;double gcNE+1=0.0,5.0,5.0,5.0,5.0,5.0

2、,5.0; /* 杆长数组*/double gjNE+1=0.0,90.0,90.0,0.0,90.0,90.0,0.0; /* 杆角数组*/double mjNE+1=0.0,0.5,0.5,0.5,0.5,0.5,0.5; /* 面积数组*/double gxNE+1=0.0,I0,I0,I0,I0,I0,I0; /* 各杆惯性矩数组*/int zcNZ+1=0,1,2,3,16,17,18;double pfNPF+15=0,0,0,0,0,0,-45,5.0,6.0,5.0,pjNPJ+13=0.0,0.0,0.0;double kzNJ3+1NJ3+1,pNJ3+1;double p

3、e7,f7,f07,t77;double ke77,kd77;/*kz-整体刚度矩阵/*ke-整体坐标下的单元刚度矩阵/*kd-局部坐标下的单元刚度矩阵/* t-坐标变换/*这是函数的声明void jdugd(int);void zb(int);void gdnl(int);void dugd(int);/*主程序开始void main(void)int i,j,k,e,d,dh,h,ii,jj,hz,a1,b1,m,l,dl,zl,z,j0;double cl,wy7;int IM,IN,jn;文件定义if(fp=fopen(c:youxianyuan.txt,w)=NULL) printf

4、(non);elseprintf(yesn);fp=fopen(c:youxianyuan.txt,w);/功能：形成矩阵pif(NPJ0)for(i=1;i0)for(i=1;i=NPF;i+)hz=i;gdnl(hz);e=(int)pfhz3;zb(e);for(j=1;j=6;j+)pej=0.0;for(k=1;k*for(e=l;e=NE;e+)dugd(e);for(i=l;i=2;i+)for(ii=l;ii=3;ii+)h=3*(i-l)+ii;dh=3*(jmei-l)+ii;for(j=l;j=2;j+)for(jj=l;jjO)kzdhdl=kzdhdl+kehl;*引

5、入边界条件*for(i=l;i=NZ;i+)z=zci;kzzl=1.0;for(j=2;jDD)jO=DD;else if(z=DD)jO=z;for(j=2;j=j0;j+)kzz-j+lj=O.O;pz=O.O;/*高斯消元法解方程组*/*消元*/for(k=1;kk+DD-1)IM=k+DD-1;else if(NJ3=k+DD-1)IM=NJ3;IN=k+1;for(i=IN;i=IM;i+)l=i-k+1;cl=kzkl/kzk1;jn=DD-l+1;for(j=1;j=1;i-)if(DDNJ3-i+1)j0=NJ3-i+1;elsej0=DD;for(j=2;j=j0;j+)h

6、=j+i-1;pi=pi-kzij*ph;pi=pi/kzi1;printf(n);printf(n) ;printf(NJ=U=V=CETA= n);fprintf(fp,n);fprintf(fp,n);fprintf(fp,NJ=U=V=CETA= n);for(i=1;i=NJ;i+)printf(%-9d %-12.11f %-12.11f %-12.11fn”,i,p3*i-2,p3*i-1,p3*i);fprintf(fp,”%-9d %-12.11f %-12.11f %-12.11fn”,i,p3*i-2,p3*i-1,p3*i);printf(n) ;fprintf(fp,

7、n);/*根据E的值输出相应E单元的N,Q,M(A,B)的结果*printf(E=N=Q=M= n);fprintf(fp,E=N=Q=M= n);计算轴力N,剪力Q,弯矩M*/for(e=1;e=NE;e+)jdugd(e);zb(e);for(i=1;i=2;i+)for(ii=1;ii=3;ii+)h=3*(i-1)+ii;dh=3*(jmei-1)+ii;wyh=pdh;for(i=1;i=6;i+)fi=0.0;for(j=1;j=6;j+)for(k=1;k0)for(i=1;i=NPF;i+) if(pfi3=e)hz=i;gdnl(hz);for(j=1;j=6;j+)fj=f

8、j+f0j;printf(%-4d(A) fprintf(fp,%-4d(A) printf( (B) fprintf(fp, (B)%-9.5f%-9.5f%-9.5f%-9.5f%-9.5f%-9.5f%-9.5f%-9.5f%-9.5fn,e,f1,f2,f3);%-9.5fn,e,f1,f2,f3);%-9.5fn,f4,f5,f6);%-9.5fn,f4,f5,f6);return;/*主程序结束*/gdnl()函数： void gdnl(int hz)int ind,e;double g,c,l0,d;g=pfhz1;c=pfhz2;e=(int)pfhz3;ind=(int)pf

9、hz 4;l0=gce;d=l0-c;f01=0.0;f02=7/20*g;f03=-g*l0/20;f04=0.0;f05=-3*g/20;f06=g*l0/30;/zb()函数： 功能：构成坐标变换矩阵void zb(int e)double ceta,co,si;int i,j;ceta=(gje*PI)/180;co=cos(ceta);si=sin(ceta);t11=co;t12=si;t21=-si;t22=co;t33=1.0;for(i=1;i=3;i+)for(j=1;j=3;j+)ti+3j+3=tij;/jduge()函数:void jdugd(int e)double a0,l0,j0;int i,j;a0=mje;l0=gce;j0=gxe;for(i=0;i=6;i+)for(j=0;j=6;j+)kdij=0.0;kd11=E0*a0/l0;kd22=12*E0*j0/pow(l0,3);kd32=6*E0*j0/pow(l0,2);kd33=4*E0*j0/l0;kd41=-kd11;kd44=kd11;kd52=-kd22;kd53=-kd32;kd55=kd22;kd62=kd32;kd63=2*E0*j0/l0;kd65=-kd32;kd66=kd33;for(i=1;i=6;i+)for(j=1;j=i;j+)kdji=