下载文档原格式
八节点等参单元平面有限元分析程序土木工程学院2011.2目录1.概述 (1)2.编程思想 (2)2.1.八节点矩形单元介绍 (2)2.2.有限元分析的模块组织 (5)3.程序变量及函数说明 (6)3.1.主要变量说明： (6)3.2.主要函数说明 (7)4.程序流程图 (8)5.程序应用与ANSYS分析的比较 (9)5.1.问题说明 (9)5.2.ANSYS分析结果 (10)5.3.自编程序分析结果 (12)5.4.结果比较分析 (12)参考文献 (14)附录程序源代码 (15)《计算力学》课程大作业1.概述通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。
具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。
随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。
有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。
因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。
因此，必须寻找新的编程语言。
随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。
现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。
C语言最初是为操作系统、编译器以及文字处理等编程而发明的。
5．等参单元本章包括以下内容： 5.1等参单元的基本概念 5.2四边形八节点等参单元 5.3等参单元的单元分析 5.4六面体等参单元5.1等参单元的基本概念在进行有限元分析时，单元离散化会带来计算误差，主要采用两种方法来降低单元离散化产生的误差：1）提高单元划分的密度，被称为h 方法（h-method ）；2）提高单元位移函数多项式的阶次，被称为p 方法（p-method ）。
在平面问题的有限单元中，我们可以选择四结点的矩形单元，如图5-1所示，该矩形单元在x 及y 方向的边长分别为2a 和2b 。
图5-1 四结点矩形单元同第三章的方法类似，将单元的位移模式选为，xy a y a x a a u 4321+++= xy a y a x a a v 8765+++=（5-1）可得到，p p m m j j i i u N u N u N u N u +++=p p m m j j i i v N v N v N v N v +++=（5-2）形态函数为， )1)(1(41b y a x N i --=)1)(1(41b y a x N j -+=)1)(1(41b y a x N m ++= )1)(1(41by ax N p +-=（5-3）上述单元位移模式满足位移模式选择的基本要求： 1）反映了单元的刚体位移和常应变， 2）单元在公共边界上位移连续。
在矩形单元的边界上，坐标x 和y 的其中一个取常量，因此在边界上位移是线性分布的，由两个结点上的位移确定。
与三结点三角形单元相比，四结点矩形单元的位移模式是坐标的二次函数，能够提高计算精度，但也有显著的缺点，两种单元的比较如下。
表5-1 三结点三角形单元与四结点矩形单元比较如果任意形状的四边形四结点单元采用矩形单元的位移模式，则在公共边界上不满足位移连续性条件。
为了既能得到较高的计算精度，又能适应复杂的边界形状，可以采用坐标变换。
图5-2任意四结点四边形单元图5-3四结点正方形单元在图5-2所示的任意四边形单元上，用等分四条边的两族直线分割四边形，以两族直线的中心为原点，建立局部坐标系),(ηξ，沿ξ及η增大的方向作为ξ轴和η轴，并令四条边上的ξ及η值分别为1±。
第二节 八节点等参数单元一、等参数单元图2.1八节点等参单元现在介绍一种常用的八结点等参数平面单元。
八结点等参单元的构成，为处理结构的曲边边界提供了优良的条件．在单元划分时内部单元可取为八结点直四边形元，边界单元的边界边可取为曲边，这相当于用三点构成的抛物线去逼近原结构曲线边界，这要比用三角形元和任意四边形单元逼近曲线边界减少离散化过程带来的误差．因此八结点等参元的引入不仅可提高单元内部插值精度，还能较好地处理曲线边界．首先研究一个边长等于2的正方形单元(图2.1)。
在其形心处按置一个局部坐标系ξη，单元各结点的坐标(i i ,ξη)分别为士1或0,因此单元四条边界的方程式都能用简单的公式写出。
选取位移模式:222212345678u a a a a a a a a ξηξηξηξηξη=+++++++222212345678v b b b b b b b b ξηξηξηξηξη=+++++++ (2.1)根据8个节点的位移信息，可确定形函数：()()()()000011111,2,3,44i N i ξηξη=+++-= ()()()201115,72i N i ξη=-+=()()()201116,82i N i ηξ=-+=上面的形函数可统一写成()()()22000011114i i i N ξηξηξη=+++- ()()()222011112i i ξξηη+--+()()()222011112i i ηηξξ+--+ (2.2)同前一样，用上式的形状函数构成真实单元的位移函数()1,n i i i u N u ξη==∑, ()1,ni i i v N v ξη==∑ (2.3)和坐标变换式()1,n i i i x N x ξη==∑, ()1,ni i i y N y ξη==∑ (2.4)通过上式的坐标变换，使得图2.1 (b)所示,ξη平面上的八个结点分别映射成图2.1 (a)所示xy 平面上的八个结点，它们的坐标是,,i i x y (i=1，2，…，8)。
program stress!character*64 fname1,fname2dimension ake(16,16)dimension ld(1000),is(16),a(30000)common/sol1/nn,ne,nm,nccommon/sol2/nf,nd,n,nfdcommon/shu1/mea(200,9)common/shu2/x(1000),y(1000)common/shu3/aeu(10,3)common/shu4/p(1000)common/shu5/jz(1000,2)!write (*,'(a22)') 'input data file:'!read(*,'(a)') fname1!write (*,'(a22)') 'output data file:'!read(*,'(a)') fname2open(1,file='fname1_re1.txt')open(2,file='fname2_re2.txt')read(1,*) nn,ne,nm,ncnf=2nd=8nfd=nf*ndn=nn*nfdo 50 i=1,nn50 read(1,*) k,x(i),y(i),p(2*i-1),p(2*i)read(1,*) ((jz(i,j),j=1,2),i=1,nc)do 100 i=1,neread(1,*) ie,(mea(i,j),j=1,9)100 continueread(1,*) ((aeu(i,j),j=1,3),i=1,nm)!400 format(4i5)!405 format(i5,4f15.6)!410 format(5x,i5,5x,i5)!415 format(5x,i5,5x,9i5)!420 format(3f15.6)write(2,460) nn,ne,nm,ncwrite(2,465) (i,x(i),y(i),p(2*i-1),p(2*i),i=1,nn)write(2,470) ((jz(i,j),j=1,2),i=1,nc)write(2,475) (i,(mea(i,j),j=1,9),i=1,ne)write(2,480) (i,(aeu(i,j),j=1,3),i=1,nm)460 format(/5x,'总的节点数，单元数，材料数，约束数'//(10x,4i5))465 format(/5x,'节点号',5x,'横坐标',5x,'纵坐标',5x,'节点载荷'//(5x,i5,4f15.6))470 format(/5x,'约束节点号及节点类型'//(5x,i5,5x,i5))475 format(/'节点编号矩阵：单元号',15x,'八个节点的整体编号',12x,'材料号'//(10x,i5,8i5,5x,i5))480 format(/3x,'材料的材料号',5x,'杨氏模量',10x,'泊松比',10x,'厚度'//(5x,i5,5x,3e15.6))call fld(nt,ld)do 500 i=1,nt500 a(i)=0do 600 ie=1,necall ke(ie,ake)call fis(ie,is)do 560 i=1,nfddo 560 j=1,nfdif(is(i)-is(j)) 560,520,520520 ni=is(i)ij=ld(ni)-ni+is(j)a(ij)=a(ij)+ake(i,j)560 continue600 continuecall fcc(nt,ld,a)call decom(nt,a,ld)!write(*,*) 'kakka'write(2,850) (i,p(2*i-1),p(2*i),i=1,nn)!write(*,*) 'kakka'call str!write(*,*) 'kakka'850 format(//25x,'节点位移'//5x,'node no.',12x,'u',12x,'v',/(7x,i5,2e15.7))end!c 应力输出子程序subroutine strcommon /sol1/nn,ne,nm,nccommon/shu3/aeu(10,3)common/shu4/p(1000)common/bmdmp/bmatx(3,16),dmatx(3,3)common/shu1/mea(200,9)dimension s(3,16),st(3),ve(16),is(16)dimension ss(8),tt(8)data ss/-1,1,1,-1,0,1,0,-1/data tt/-1,-1,1,1.,-1,0,1,0/call zero(3,3,dmatx)call zero(3,16,bmatx)do 800 ie=1,newrite(2,850) ienm1=mea(ie,9)yong=aeu(nm1,1)poiss=aeu(nm1,2)thick=aeu(nm1,3)call modps(yong,poiss)call fis(ie,is)do 750 i=1,16ni=is(i)ve(i)=p(ni)750 continuedo 300 i=1,8exisp=ss(i)etasp=tt(i)k=mea(ie,i)call sfr2(exisp,etasp)call jacob2(ie,djacb)call bmatpscall dot(3,3,16,dmatx,bmatx,s)call dot(3,16,1,s,ve,st)write(2,880) i,k,st(1),st(2),st(3)c1=(st(1)+st(2))/2.0c2=sqrt((st(1)-st(2))**2/4.0+st(3)**2)psi=c1+c2psj=c1-c2s1=amax1(psi,psj,0.0)s3=amin1(psi,psj,0.0)s2=st(1)+st(2)-s1-s3!write(*,*) s1,s2,s3c3=(s1-s2)**2+(s2-s3)**2+(s1-s3)**2seq4=sqrt(c3/2.0)!write(2,*) seq4300 continue800 continue880 format(i5,5x,i5,5x,f15.5,5x,f15.5,5x,f15.5)850 format(//10x,'第',i3,'个单元的节点应力'/2x,'局部编号',2x,'整体编号',5x,'u方向应力',12x,'v方向应力',12x,'剪应力')returnend! c 这是一个形成局部编号与整体编号对应关系is数组的子程序SUBROUTINE FIS(IE,IS)common/shu1/mea(200,9)dimension is(nfd)DO 150 ID=1,NDDO 100 JF=1,NFI1=(ID-1)*NF+JFIS(I1)=(MEA(IE,ID)-1)*NF+JF100 CONTINUE150 CONTINUERETURNEND!c 这是一个形成变带宽存储指示数组ld的子程序SUBROUTINE FLD(NT,LD)common /sol1/nn,ne,nm,nccommon/sol2/nf,nd,n,nfdcommon/shu1/mea(200,9)dimension ld(n)LD(1)=1DO 300 K=1,NNNMIN=1000DO 200 I=1,NEDO 150 J=1,NDIF(MEA(I,J).NE.K) GOTO 150DO 100 L=1,NDIF (MEA(I,L).LT.NMIN) NMIN=MEA(I,L) 100 CONTINUE150 CONTINUE200 CONTINUEDO 250 I=1,NFJ=(K-1)*NF+IIF(J.NE.1) LD(J)=LD(J-1)+(K-NMIN)*NF+I 250 CONTINUE300 CONTINUENT=LD(N)RETURNEND! c 这是一个处理约束的子程序SUBROUTINE FCC(NT,LD,A)common /sol1/nn,ne,nm,nccommon/shu5/jz(1000,2)dimension ld(n),a(nt)DO 350 K=1,NCI=JZ(K,1)J=JZ(K,2)NI=(I-1)*NF+JNJ=LD(NI)A(NJ)=1.0E25350 CONTINUERETURNEND! c 这是一个解线性方程组的子程序SUBROUTINE DECOM(NT,A,LD)common/shu4/p(1000)common/sol2/nf,nd,n,nfddimension a(nt),ld(n)DO 20 I=1,NLDI=LD(I)IF (I.EQ.1) GOTO 10IO=LDI-IIM1=I-1MI=1-IO+LD(IM1)IF(MI.GT.IM1) GOTO 10DO 8 J=MI,IM1JO=LD(J)-JIJ=IO+JIF(J.EQ.1) GOTO 6JM1=J-1MJ=1-JO+LD(JM1)MIJ=MAX0(MI,MJ)IF(MIJ.GT.JM1) GOTO 6S=0.0DO 2 K=MIJ,JM1IK=IO+KJK=JO+KS=S+A(IK)*A(JK)2 CONTINUEA(IJ)=A(IJ)-S6 p(I)=p(I)-A(IJ)*p(J)8 CONTINUE10 ALDI=A(LDI)IF(I.EQ.1.OR.MI.GT.IM1) GOTO 13DO 12 J=MI,IM1IJ=IO+JLDJ=LD(J)S=A(IJ)A(IJ)=S/A(LDJ)ALDI=ALDI-S*A(IJ)12 CONTINUE13 A(LDI)=ALDIp(I)=p(I)/ALDI20 CONTINUEDO 40 LDI=2,NI=N-LDI+2IO=LD(I)-IMI=1-IO+LD(I-1)IM1=I-1IF(MI.GT.IM1) GOTO 40DO 30 J=MI,IM1IJ=IO+Jp(J)=p(J)-A(IJ)*p(I)30 CONTINUE40 CONTINUERETURNEND! c 这是一个形成单元刚度矩阵的子程序subroutine ke(ie,ake)dimension posgp(2)common/sol1/nn,ne,nm,nccommon/bmdmp/bmatx(3,16),dmatx(3,3) common/shu3/ aeu(10,3)common/shu1/mea(200,9)common/shu2/x(1000),y(1000) dimension bt(16,3),bd(16,3)dimension ake(16,16),ake1(16,16) posgp(1)=-0.5773502692posgp(2)=0.5773502692nm1=mea(ie,9)young=aeu(nm1,1)poiss=aeu(nm1,2)thick=aeu(nm1,3)call modps(young,poiss)call zero(16,16,ake)do 80 iagus=1,2do 80 jagus=1,2exisp=posgp(iagus)etasp=posgp(jagus)call sfr2(exisp,etasp)call jacob2(ie,djacb)call bmatpscall tran(3,16,bmatx,bt)call dot(16,3,3,bt,dmatx,bd)call dot(16,3,16,bd,bmatx,ake1)do 70 i=1,16do 70 j=1,1670 ake(i,j)=ake(i,j)+djacb*thick*ake1(i,j)80 continue!write(2,*) ie!do 90 i=1,16!write(*,*) ie!write(2,10) (ake(i,j),j=1,16)! 90 continue! 10 format(16e15.3)returnend! c 这是一个形成弹性矩阵的子程序subroutine modps(y,p)common/bmdmp/bmatx(3,16),dmatx(3,3) call zero(3,3,dmatx)!write(*,*) y,pconst=y/(1.0-p*p)dmatx(1,1)=constdmatx(2,2)=constdmatx(1,2)=const*pdmatx(2,1)=dmatx(1,2)dmatx(3,3)=const*(1.0-p)/2.0! do 100 i=1,3!write(*,*) (dmatx(i,j),j=1,3)!100 continuereturnend! c 这是一个计算形函数当前积分点及形函数对局部坐标的导数值的子程序subroutine sfr2(s,t)common/sd/sh(8),deriv(2,8)s2=s*2.0t2=t*2.0ss=s*stt=t*tst=s*tstt=s*t*tsst=s*s*tst2=s*t*2.0sh(1)=(-1.0+st+ss+tt-sst-stt)/4.0sh(2)=(-1.0-st+ss+tt-sst+stt)/4.0sh(3)=(-1.0+st+ss+tt+sst+stt)/4.0sh(4)=(-1.0-st+ss+tt+sst-stt)/4.0sh(5)=(1.0-t-ss+sst)/2.0sh(6)=(1.0+s-tt-stt)/2.0sh(7)=(1.0+t-ss-sst)/2.0sh(8)=(1.0-s-tt+stt)/2.0deriv(1,1)=(t+s2-st2-tt)/4.0deriv(1,2)=(-t+s2-st2+tt)/4.0deriv(1,3)=(t+s2+st2+tt)/4.0deriv(1,4)=(-t+s2+st2-tt)/4.0deriv(1,5)=-s+stderiv(1,6)=(1.0-tt)/2.0deriv(1,7)=-s-stderiv(1,8)=(-1.0+tt)/2.0deriv(2,1)=(s+t2-ss-st2)/4.0deriv(2,2)=(-s+t2-ss+st2)/4.0deriv(2,3)=(s+t2+ss+st2)/4.0deriv(2,4)=(-s+t2+ss-st2)/4.0deriv(2,5)=(-1.0+ss)/2.0deriv(2,6)=-t-stderiv(2,7)=(1.0-ss)/2.0deriv(2,8)=-t+streturnend! c 这是一个形成jacobi矩阵的子程序subroutine jacob2(ie,djacb)dimension xjacm(2,2),xjaci(2,2)common/sd/sh(8),deriv(2,8)common/car/cartd(2,8)common/shu1/mea(200,9)common/shu2/x(1000),y(1000)do 200 j=1,2xjacm(1,j)=0.0xjacm(2,j)=0.0do 200 k=1,8if(j.eq.1) thenxjacm(1,j)=xjacm(1,j)+deriv(1,k)*x(mea(ie,k))xjacm(2,j)=xjacm(2,j)+deriv(2,k)*x(mea(ie,k))elsexjacm(1,j)=xjacm(1,j)+deriv(1,k)*y(mea(ie,k))xjacm(2,j)=xjacm(2,j)+deriv(2,k)*y(mea(ie,k))endif200 continuedjacb=xjacm(1,1)*xjacm(2,2)-xjacm(1,2)*xjacm(2,1)!write(*,*) xjacm(1,1) ,xjacm(1,2),xjacm(2,1),xjacm(2,2)!write(*,*) ie!write(*,*) djacbif(djacb.lt.1.e-6) thenwrite(6,6100) ie6100 format('单元编号为:',i4,'的jcobi行列式的值小于或等于零!') stop '程序运行被中断于子程序jacob2'endifxjaci(1,1)=xjacm(2,2)/djacbxjaci(2,2)=xjacm(1,1)/djacbxjaci(1,2)=-xjacm(1,2)/djacbxjaci(2,1)=-xjacm(2,1)/djacbdo 300 i=1,2do 300 k=1,8cartd(i,k)=0.0do 300 j=1,2cartd(i,k)=cartd(i,k)+xjaci(i,j)*deriv(j,k)300 continuereturnend! c 这是一个形成应变矩阵的子程序subroutine bmatpscommon/car/cartd(2,8)common/bmdmp/bmatx(3,16),dmatx(3,3)call zero(3,16,bmatx)ngash=0do 200 i=1,8mgash=ngash+1ngash=mgash+1bmatx(1,mgash)=cartd(1,i)bmatx(1,ngash)=0.0bmatx(2,mgash)=0.0bmatx(2,ngash)=cartd(2,i)bmatx(3,mgash)=cartd(2,i)bmatx(3,ngash)=cartd(1,i)200 continuereturnendSUBROUTINE ZERO(M,N,A)DIMENSION A(M,N)DO 10 I=1,MDO 10 J=1,NA(I,J)=0.010 CONTINUERETURNENDSUBROUTINE TRAN(M,N,A,B)DIMENSION A(M,N),B(N,M)DO 20 I=1,MDO 20 J=1,NB(J,I)=A(I,J)20 CONTINUERETURNENDSUBROUTINE DOT(M,N,L,A,B,C)DIMENSION A(M,N),B(N,L),C(M,L)DO 50 I=1,MDO 50 J=1,LC(I,J)=0.0DO 40 K=1,NC(I,J)=C(I,J)+A(I,K)*B(K,J)40 CONTINUE50 CONTINUERETURNEND。
cpe8单元等参元形状函数-回复什么是cpe8单元等参元形状函数？在计算力学领域中，等参元是指一种在有限元分析中使用的数学函数。
这些函数用于近似描述结构物体的形状或者表示变量的分布。
CPE8单元是指一个八节点的六面体元素，其中CPE代表六面体元素的类型，8代表该元素具有八个节点。
等参元形状函数则是用于表示这些节点在六面体元素中的几何属性和作用的数学函数。
等参元形状函数的主要作用是描绘和定义六面体元素的形状和几何特征。
这些函数由节点坐标和插值方法决定，可以通过非常重要的节点坐标来表示六面体元素的形状。
等参元形状函数通常是多项式形式，方便于在有限元分析中进行数值计算。
对于CPE8单元而言，它的八个节点可以通过等参元形状函数来定义。
每个节点的坐标用局部坐标系表示，用参数ξ，η和ζ表示。
在CPE8单元中，等参元形状函数通常用Lagrange插值多项式表示，这些多项式是基于节点坐标的线性和二次函数。
具体来说，CPE8单元的等参元形状函数通常包括一个常数项、一个线性项和一个二次项。
这些形状函数的具体表达式可以通过如下方式计算：N1(ξ,η,ζ) = a + bξ+ cη+ dζ+ eξη+ fηζ+ gξζ+ hξηζN2(ξ,η,ζ) = a' + b'ξ+ c'η+ d'ζ+ e'ξη+ f'ηζ+ g'ξζ+ h'ξηζN3(ξ,η,ζ) = a'' + b''ξ+ c''η+ d''ζ+ e''ξη+ f''ηζ+ g''ξζ+ h''ξηζ...N8(ξ,η,ζ) = a''' + b'''ξ+ c'''η+ d'''ζ+ e'''ξη+ f'''ηζ+ g'''ξζ+ h'''ξηζ在这些表达式中，N1至N8代表八个节点的等参元形状函数，而a至h、a'至h'''则是与对应节点相关的系数。
