四节点单元刚度矩阵fortran编程

平面四节点单元的刚度矩阵的fortran编程parameter(E0=2.0e06,u=0.3,h=2)!参数的选定real,dimension(2,2)::kii,kij,kji,kjj,kim,kin,kjm,kjn,kmi,kmj,kmm,kmn,kni,knj,knm,knnreal,dimension(3,8)::B!单元总应变矩阵real,dimension(2,4)::z！节点的坐标real,dimension(3,2)::Bi,Bj,Bm,Bn！单元应变矩阵real,dimension(8,8)::Kee!单元刚度矩阵real,dimension(3,3)::D!矩阵D为3×3矩阵的定义real,x,ydataz/xi,yi,xj,xm,ym,xn,yn/a=z(1,2)-z(1,1)=xj-xib=z(2,3)-z(2,2)=ym-yjNi=0.25*(1+xi*x/a)*(1+yi*y/b)！形函数的定义Nj=0.25*(1+xj*x/a)*(1+yj*y/b)Nm=0.25*(1+xm*x/a)*(1+ym*y/b)Nn=0.25*(1+xn*x/a)*(1+yn*y/b)Bi(1,1)=xi*(1+yi*y/b)/(4*a)！Bi矩阵的定义Bi(1,2)=0Bi(2,1)=0Bi(2,2)=yi*(1+xi*x/a)/(4*b)Bi(3,1)=yi*(1+xi*x/a)/(4*b)Bi(3,2)=xi*(1+yi*y/b)/(4*b)Bj(1,1)=xj*(1+yj*y/b)/(4*a)Bj(1,2)=0Bj(2,1)=0Bj(2,2)=yj*(1+xj*x/a)/(4*b)Bj(3,1)=yj*(1+xj*x/a)/(4*b)Bj(3,2)=xj*(1+yj*y/b)/(4*b)Bm(1,1)=xm*(1+ym*y/b)/(4*a)Bm(1,2)=0Bm(2,1)=0Bm(2,2)=ym*(1+xm*x/a)/(4*b)Bm(3,1)=ym*(1+xm*x/a)/(4*b)Bm(3,2)=xm*(1+ym*y/b)/(4*b)!Bm矩阵的定义Bn(1,1)=xn*(1+yn*y/b)/(4*a)Bn(1,2)=0Bn(2,1)=0Bn(2,2)=yn*(1+xn*x/a)/(4*b)Bn(3,1)=yn*(1+xn*x/a)/(4*b)Bn(3,2)=xn*(1+yn*y/b)/(4*b)!Bn矩阵的定义D(1,1)=(E0/(1-u*u))*1!D矩阵的定义D(2,2)=(E0/(1-u*u))*1D(1,2)=(E0/(1-u*u))*uD(2,1)=(E0/(1-u*u))*uD(3,1)=0D(3,2)=0D(1,3)=0D(2,3)=0D(3,3)=(E0/(1-u*u))*((1-u)/2)K0ii=matmul(transpose(Bi),D,Bi)！矩阵相乘K0ij=matmul(transpose(Bi),D,Bj)！transpose（Bi）为矩阵的转置K0im=matmul(transpose(Bi),D,Bm)K0in=matmul(transpose(Bi),D,Bn)K0ji=matmul(transpose(Bj),D,Bi)K0jj=matmul(transpose(Bj),D,Bj)K0jm=matmul(transpose(Bj),D,Bm)K0jn=matmul(transpose(Bj),D,Bn)K0mi=matmul(transpose(Bm),D,Bi)K0mj=matmul(transpose(Bm),D,Bj)K0mm=matmul(transpose(Bm),D,Bm)K0mn=matmul(transpose(Bm),D,Bn)K0ni=matmul(transpose(Bn),D,Bi)K0nj=matmul(transpose(Bn),D,Bj)K0nm=matmul(transpose(Bn),D,Bm)K0nn=matmul(transpose(Bn),D,Bn)!weijifenqianKeejuzhenzijuzhensubroutinejfen(f0,a1,b1)！一重积分辛普森的子程序parameter(n=100)reala1,b1,h,x,s,fo,f1,siintegern,ih=(b1-a1)/nx=a;s=0doi=1,nx=x+hf1=1.0/(1+x)si=(f0+f1)*h/2s=s+sif0=f1enddoendsubroutinedoi=1,2doj=1,2calljfen(fs,a,b,f,eps,s)！调用fortran的二重积分子程序Kii(i,j)=jfen(K0ii(i,j),a,b,a,0.001,b)！对Kii的每个元素进行二重积分的运算Kij(i,j)=jfen(K0ij(i,j),a,b,a,0.001,b)Kim(i,i)=jfen(K0im(i,j),a,b,a,0.001,b)Kin(i,j)=jfen(K0in(i,j),a,b,a,0.001,b)Kji(i,j)=jfen(K0ji(i,j),a,b,a,0.001,b)Kjj(i,j)=jfen(K0jj(i,j),a,b,a,0.001,b)Kjm(i,j)=jfen(K0jm(i,j),a,b,a,0.001,b)Kjn(i,j)=jfen(K0jn(i,j),a,b,a,0.001,b)Kmi(i,j)=jfen(K0mi(i,j),a,b,a,0.001,b)Kmj(i,j)=jfen(K0mj(i,j),a,b,a,0.001,b)Kmm(i,j)=jfen(K0mm(i,j),a,b,a,0.001,b)Kmn(i,j)=jfen(K0mn(i,j),a,b,a,0.001,b)Kni(i,j)=jfen(K0ni(i,j),a,b,a,0.001,b)Knj(i,j)=jfen(K0nj(i,j),a,b,a,0.001,b)Knm(i,j)=jfen(K0nm(i,j),a,b,a,0.001,b)Knn(i,j)=jfen(K0nn(i,j),a,b,a,0.001,b)enddoenddoKee=merge(kii,kij,kji,kjj,kim,kin,kjm,kjn,Kmi,Kmj,Kmm,Kmn,Kni,Knj,Knm,Knn,mask)!对各个doi=1,8doj=1,8print*,Kee(i,j)！Kee的输出enddoenddoendsubroutinejfen(fs,a,b,f,eps,s)！在网上找出的二重积分计算的子程序externalfs,f！子程序参考徐士良老师编写的子程序集doubleprecisiona,b,f,s,h,s1,s2,ts1,ts2,x,g,s0,cn=1h=0.5*(b-a)c=(b-a)*1.0e-6callsimp1(a,fs,f,eps,s1)callsimp1(b,fs,f,eps,s2)ts1=h*(s1+s2)10x=a-hTs2=0.5*ts1Do20j=1,nX=x+2*hCallsimp1(x,fs,f,eps,g)Ts2=ts2+h*g20continueS=(4*ts2-ts1)/0.3n=n+nif(n.ge.16)thenif(abs(s-s0).le.eps*abs(s)+1.0))returnendifs0=sts1=ts2h=0.5*hif(abs(h).ge.abs(c))goto10returnendsubroutinesimp1(x,fs,f,eps,g)doubleprecisionx,f,g,y1,y2,h,c,ts1,ts2,y,g0n=1callfs(x,y1,y2)h=0.5*(y2-y1)c=(y2-y1)*1.0e-06ts1=h*(f(x,y1)+f(x,y2))10y=y1-hts2=0.5*ts1do20i=1,ny=y+2*hts2=ts2+h*f(x,y)20continueG=(4*ts2-ts1)/3.0n=n+nif(n.ge.16)thenif(abs(g-g0).le.eps*(abs(g)+1.0))returnendifg0=gts1=ts2h=0.5*hif(abs(h).ge.abs(c))goto10returnend