在集中荷载作用下水泥混凝土路面板的有限元分析
李志勇          410076
现有一水泥混凝土路面板，宽为4m，厚0.2m，在板上作用一行车载荷，如图一所示，其荷载大小为50KN，材料参数为E=3×104，泊松比v=0.3，密度ρ=2.4×104。因为路面非常狭长，所以取路面的一个横截面作为研究对象，因此此问题属于平面应变问题。有考虑到对称性，可取一半来研究。建立有限元分析模型，并划分单元网格。采用八节点等参单元，单元划分如图二所示，单元数为10，结点数为53，边界节点数为44。

                   P=50KN          P
                               1m


                                               0.2m
                           4m
图一


                                  图二
有限元方法是近三十年来发展起来的一种借助计算机程序来解决工程力学问题的数值计算方法。三十多年来，有限元法的应用由弹性力学平面问题扩展到空间问题、板壳问题；由静力平衡问题扩展到稳定问题、动力问题和波动问题。分析的对象从弹性材料扩展到塑性、粘弹性、粘塑性和复合材料等。有限元的基本思想是将连续的求解区域离散为一组有限个且按一定方式联结在一起的单元的组合体。
一、        基本原理
有限元方法先把结构先分解为有限个较小的单元，即进行所谓的离散化。结构离散化的目的是在较小的范围内分析单元的内力与位移的关系运用虚功原理建立单元刚度矩阵，然后把各单元又集合成原来的结构，这就要求各单元满足原结构的几何条件（包括支承条件和结点处的变形边疆条件）和平衡条件，从而建立整个结构的刚度方程，以求解原结构的位移和应力。
针对该问题的特点，采用八节点等参单元是因为它具有良好的单元特性及对曲线边界的适应性。如图二 (a)、(b)所示。
图二  八结点等参单元
１.位移插值函数（形函数）Ni
                                
    (i=1，…８)                 (1)
式中： ——为单元内任一点位移；
——为已知结点位移；
——为插值函数。
对于八结点等参单元，插值函数为：
                          
                          
                          
                          
                          
                          
                          
                                                      （2）
二、        计算步骤
1)有限元网格的划分，输入计算参数；如划分单元数，结点数及坐标，荷载位置和大小，结构层材料参数等；
  2)计算单元刚度矩阵并形成总刚矩阵；
  3)引入结点荷载与自重；
  4)引入支撑条件；
  5)解方程组，求出结点位移；
  6)计算单元应力；
7)进行精度校核；
8)输出计算结果。
三、        有限元源程序及数据文件
有限元源程序
dimension a(700,100),b(700),xy(350,2),kk1(100,9),kk2(60),
        +         c(10,3),ne(8),xi(8),yi(8),t(3),f(40,2),pp(40,2),np(20)
        open (5,file='f1.txt',status='old')
        open (6,file='f2.txt',status='unknown')


        read (5,*)  nma        !nma ：材料种类数
        do 11 i=1,nma
11    read (5,*) (c(i,j),j=1,3) !        输入各种材料的参数：弹性模量、泊松比、密度。
      write (6,500) nma !输出材料种类数
        
500        format (10x,'the total number of materials:' ,i10)
        write (6,512)
512      format (10x,'the material parameters: ' /10x,'E',15x,
     +        'v',15x, 'r')
        do 513 i=1,nma
        write (6,514) (c(i,j), j=1,3) !输出各材料参数
514   format (5x,3(e10.4,5x))
513   continue
      read (5,*) na        !na为节点总数
        do 12 i=1,na
12    read (5,*) (xy(i,j),j=1,2)!循环输入各节点坐标。xy(i,1)为x坐标值        
c                                                       xy(i,2)为y坐标值
        
             write(6,501) na        !输出节点总数
501   format(10x,'the joint number in FEM net:',i10)
             write (6,502)
502   format(5x,'the joint ordinates'/10x,'x',20x,'y')
             do 503  i=1,na        
       write(6,504) (xy(i,j),j=1,2)  !输出各节点坐标值
504         format(5x,f10.4,10x,f10.4)          
503   continue
             read(5,*) n1  !n1为单元总数
             do 13 i=1,n1
13    read(5,*)  (kk1(i,j),j=1,9)!输入单元节点编号
c         　（8节点，逆时针），最后一列为材料编号
                    
      write (6,505) n1 !输出单元总数
505   format (10x,'the element number in FEM net:',i10)
      write (6,506)
506   format (2x,'NO.', 60x,'mat.aera')
      do 507 i=1,n1
        write (6,508) (kk1(i,j),j=1,9) !输出节点编号及材料信息。
508   format (8(2x,i5),5x,i6)
507   continue
      read (5,*)  kb        !输入边界节点数kb
        read (5,*) (kk2(i), i=1,kb) !输入边界节点编号                        
      do 572 i=1,kb
572   read(5,*) (f(i,j),j=1,2) !输入边界节点位移值，第i列表示第i点的位移值  
         write (6,509) kb !输出边界节点数
509   format(10x,'the joint number in boundary:', i10)
       write (6,510) (kk2(i),i=1,kb)!输出边界各节点编号
510   format (/10x,i7)

    
       do 571 i=1,kb
        write (6,570) kk2(i), (f(i,j),j=1,2)!输出边界各节点位移值
570    format(10x,'No.',i2,10x,'u=',e10.4,5x ,'v=',e10.4)
571   continue
         do 10 i=1,700
        b(i)=0.0  !将各位移列向量初始化为0
        do 10 j=1,110
10    a(i,j)=0.0 ! 将刚度矩阵各元素初始化为0
      write (*,516)
516    format (10x,'belt started')
       call belt (kk1,n1,nd) !调用子程序，计算半带宽nd
        n3=2*na        !计算总行数。为节点数的2倍(每个节点有x、y)
        n4=nd+1        !计算总列数(最后一列为节点荷载)
        write (*,517)
517   format (10x,'calculating the element elastic matrix')
      do 100 i=1,n1 !n1为单元数
        do 30 i1=1,8
30     ne(i1)=kk1(i,i1) !将kk1数组中第i个单元的节点编号给ne数组
      nb=kk1(i,9) !将kk1数组中第i个单元的材料信息给nb
        do 40 i1=1,8
        nee=ne(i1)
        xi(i1)=xy(nee,1) !将xy中第i1个节点的横坐标赋给xi
40    yi(i1)=xy(nee,2) !将xy中第i1个节点的纵坐标赋给yi

    e=C(nb,1)
        v=c(nb,2)
        r=c(nb,3)
        E=E/(1.0-v*v)
        v=v/(1.0-v)
        call s0 (a,ne,xi,yi,E,v,n3,n4,r)!调用子程序，计算单元刚度矩阵
c                                      及节点荷载  
       write(6,518)  i         !显示出是第几个单元
518   format (20x,2hi=,i10)
100        continue
       read (5,*) nhzd
      do 50 i=1,nhzd
        read (5,*) np(i)
        read(5,*) (pp(i,j),j=1,2)
        write(6,51)np(i)
51         format(2x,i3)
      write(6,*) (pp(i,j),j=1,2)
50      continue
        do 60 i=1,nhzd
          i3=np(i)
        m=2*i3
      n=m-1
        a(n,n4)=a(n,n4)+pp(i,1)
        a(m,n4)=a(m,n4)+pp(i,2)
60         continue

c        write(6,*) a(m,n4)
      call s3(a,n3,n4,sm)        !计算刚度矩阵中的最大元素值sm
        write (*,519)
519   format (5x,'the max.number in matrix has been calculated')
    
      amax=sm*100000.0 !将最大值乘以大数
        call bound (a,kk2,n3,n4,amax,kb,f)!调用子程序，引入位移边界条件
        write (6,520)  
520   format (10x,'the boundary has been introduced')

    
       call gauss(a,n3,n4,nd)!调用子程序，半带宽gauss消去法解方程组。
        write (*,521)
521   format (10x,'gauss passed')
      do 350 i2=1,n3 !n3总行数
350   b(i2)=a(i2,n4)!将a数组中每行的最后一列赋给数组b
      write (6,522)
522   format(10x,'calculating the stresses of joints')
       write (5,850)
850  format  (/10x,'the stresses of joints in FEM net')
        
      do 410 i=1,na        !na为节点总数
        do 411 k=1,3
411        t(k)=0.0
      k1=0
        do 420 j=1,n1
        mb=0
        do 430 k=1,8
        ma=kk1(j,k)
        if (ma.ne.i) goto 430
        mb=k
        k1=k1+1
430   continue
      write (*,600) i,j,mb,k1
600   format (5x,'i,j,mb,k1',4i8)
      if (mb.eq.0) goto 420
        do 440 k=1,8
440   ne(k)=kk1(j,k)
      nb=kk1(j,9)
        do 450 k=1,8
        nee=ne(k)
        xi(k)=xy(nee,1)
450   yi(k)=xy(nee,2)
      E=c(nb,1)
       v=c(nb,2)
       E=E/(1-v*v)
       v=v/(1.0-v)
      call s4(b,ne,xi,yi,E,v,n3,mb,t)!调用子程序，计算节点应力
420       continue
c        if (k1.eq.0)  goto 410
        write (*,523) i,k1
523   format (10x,'joint No.:',i5,5x,'number calculated for the',
     +  'joint:',i5)
        do 460 k=1,3
460         t(k)=t(k)/k1
        td=t(3)**2+((t(1)-t(2))**2)/4.0
        write(6,470) i,(t(k),k=1,3)
470   format(2x,2hn=,i3,2x,3hsx=,e10.4,2x,3hsy=,e10.4,2x,4hsxy=,
     + e10.4)

410     continue
      stop
        end
      subroutine s0(a,ne,xi,yi,E,v,n3,n4,r)  !单元刚度矩阵子程序
        dimension a(700,110),ne(8),xi(8),yi(8),z1(3),z2(3),fnx(8),
        1  fny(8),fn(8)
        z1(1)=0.0-sqrt(0.6)
        z1(2)=0.0
        z1(3)=0.0-z1(1)
        z2(1)=5.0/9.0
        z2(2)=8.0/9.0
        z2(3)=z2(1)
        do 200 k1=1,3  !在ξ方向高斯积分
      x=z1(k1)
        do 200 k2=1,3 !在η方向高斯积分
        y=z1(k2)
        call s2(x,y,xi,yi,fnx,fny,fn,yacobi)
        write(*,700) yacobi
700   format(15x,'yacobi:',3x,e10.4)
      do 100 i=1,8
        m=ne(i)
        a(2*m,n4)=a(2*m,n4)-fn(i)*r*yacobi*z2(k1)*z2(k2) !自重产生的节点荷载
        do 100 j=1,8
        n=ne(j)
        if (m-n)110,120,100
120        M21=2*M-1
        M2=M21+1
        N1=1
        goto 150
110   m21=2*M-1
      m2=m21+1
        n1=2*N-2*M+1
150   a(m21,n1)=a(m21,n1)+z2(k1)*z2(k2)*yacobi*E*(fnx(i)*fnx(j)
     +  +0.5*(1.0-v)*fny(i)*fny(j))/(1.0-v*v)
        a(m21,n1+1)=a(m21,n1+1)+z2(k1)*z2(k2)*yacobi*E*(v*fnx(i)*fny(j)
     +  +0.5*(1.0-v)*fny(i)*fnx(j))/(1.0-v*v)
      if(m.eq.n)goto 160
       a(m2,n1-1)=a(m2,n1-1)+z2(k1)*z2(k2)*yacobi*E*(v*fny(i)    
     + *fnx(j)+0.5*(1.0-v)*fnx(i)*fny(j))/(1.0-v*v)
160   a(m2,n1)=a(m2,n1)+z2(k1)*yacobi*E*(fny(i)*fny(j)
     +    +0.5*(1.0-v)*fnx(i)*fnx(j))/(1.0-v*v)
100   continue
200   continue

        return
       end
      subroutine s2(x,y,xi,yi,fnx,fny,fn,yacobi)        !计算yacobi行列式
        dimension xi(8),yi(8),fnx(8),fny(8),fn(8),fn1(8),fn2(8)
        xa=0.0
        xb=0.0
        ya=0.0
        yb=0.0
        fn1(1)=(1.0-y)*(2.0*x+y)*0.25
        fn1(2)=(1.0-y)*(2.0*x-y)*0.25
        fn1(3)=(1.0+y)*(2.0*x+y)*0.25
        fn1(4)=(1.0+y)*(2.0*x-y)*0.25
        fn1(5)=-2.0*x*(1.0-y)*0.5
        fn1(6)=(1.0-y*y)*0.5
        fn1(7)=-2.0*x*(1.0+y)*0.5
        fn1(8)=-(1.0-y*y)*0.5        
        fn2(1)=(1.0-x)*(2.0*y+x)*0.25
        fn2(2)=(1.0+x)*(2.0*y-x)*0.25
        fn2(3)=(1.0+x)*(2.0*y+x)*0.25
        fn2(4)=(1.0-x)*(2.0*y-x)*0.25
        fn2(5)=-(1.0-x*x)*0.5
        fn2(6)=-2.0*y*(1.0+x)*0.5        
        fn2(7)=(1.0-x*x)*0.5
        fn2(8)=-2.0*y*(1.0-x)*0.5
        fn(1)=(1.0-x)*(1.0-y)*(-x-y-1.0)*0.25
        fn(2)=(1.0+x)*(1.0-y)*(x-y-1.0)*0.25
        fn(3)=(1.0+x)*(1.0+y)*(x+y-1.0)*0.25
        fn(4)=(1.0-x)*(1.0+y)*(-x+y-1.0)*0.25
        fn(5)=(1-x*x)*(1.0-y)*0.5
        fn(6)=(1.0-y*y)*(1.0+x)*0.5
        fn(7)=(1.0-x*x)*(1.0+y)*0.5
        fn(8)=(1.0-y*y)*(1.0-x)*0.5
        do 100 i=1,8
        xa=xa+fn1(i)*xi(i)
        xb=xb+fn2(i)*xi(i)
        ya=ya+fn1(i)*yi(i)
100        yb=yb+fn2(i)*yi(i)
        yacobi=xa*yb-xb*ya
        do 200 i=1,8
        fnx(i)=(yb*fn1(i)-ya*fn2(i))/yacobi
200   fny(i)=(-xb*fn1(i)+xa*fn2(i))/yacobi
        return
        end
        subroutine s3(a,n,m,sm)
        dimension a(700,110)
        rm=0.0
        sm=0.0
        do 10 i=1,n
        t=0.0
        do 20 j=1,m-1
        t=t+abs(a(i,j))
        if(abs(a(i,j)).le.rm) goto 20
        ii=i
        jj=j
        rm=abs(a(i,j))
20    continue
      if(t.le.sm)goto 10
        sm=t
10    continue
      write(*,30)ii,jj,rm
30    format(10x,6hmax a(,i3,1h,,i3,2h)=,e13.6)
      write(*,40) sm
40    format(10x,9h/a(i,j)/=,e13.6)
      return
        end

      subroutine s4(b,ne,xi,yi,E,v,n3,i,t)
        dimension b(n3),ne(8),b1(8),b2(8),xi(8),yi(8),fnx(8),fny(8),
        1  s(3),xii(8),yii(8),fn(8),t(3)
      data xii/-1.0,1.0,1.0,-1.0,0.0,1.0,0.0,-1.0/,yii/-1.0,-1.0,
        1    1.0,1.0,-1.0,0.0,1.0,0.0/
      do 100 j=1,8
        nn=ne(j)
        b1(j)=b(2*nn-1)
100        b2(j)=b(2*nn)
        nn=ne(i)
        x=xii(i)
        y=yii(i)
        s(1)=0.0
        s(2)=0.0
        s(3)=0.0
        call s2(x,y,xi,yi,fnx,fny,fn,yacobi)
        do 200 j=1,8
        s(1)=s(1)+E*(fnx(j)*b1(j)+v*fny(j)*b2(j))/(1.0-v*v)
        s(2)=s(2)+E*(v*fnx(j)*b1(j)+fny(j)*b2(j))/(1.0-v*v)
        s(3)=s(3)+0.5*E*(1.0-v)*(fny(j)*b1(j)+fnx(j)*b2(j))/(1.0-v*v)
200        continue
        do 210 j=1,3
210        t(j)=t(j)+s(j)
        return
        end
c    
    subroutine gauss(a,n3,n4,nd)  !半带宽高斯消元
        dimension a(700,110)
         n=n3
           n1=n-1
        do 10 m1=1,n1

      nd2=n4-2
        ndL=n-m1
        if(ndL.gt.nd2) ndL=nd2

      do 10 L=1,ndL
        nmL1=L+1

      a(m1+L,n4)=a(m1+L,n4)-a(m1,n4)*a(m1,nmL1)/a(m1,1)
        nd1=nd-L

      Do 20 j=1,nd1
        nmjL=j+L
      nm1Lj=m1+L
      A(nm1Lj,j)=a(nm1Lj,j)-a(m1,nml1)*a(m1,nmjl)/a(m1,1)

20      continue
10    continue

      a(n,n4)=a(n,n4)/a(n3,1)
      do 30 m1=1,n1
      mm1=n-m1
      ndd=n-mm1+1
      if (ndd.gt.nd)  ndd=nd
      do 40 j=2,ndd
      mj1=mm1+j-1
      a(mm1,n4)=a(mm1,n4)-a(mj1,n4)*a(mm1,j)        !!!!
40    continue
      a(mm1,n4)=a(mm1,n4)/a(mm1,1)
30         continue
       write(6,100)
100   format(10x,'joint displacements:')
      write(6,101)
101   format(5x,'joint No.',10x,'u',10x,'v')
      m=n3/2
      do 110 I=1,m
      j=2*(I-1)
      write(6,120)  I, a(j+1,n4),a(j+2,n4)
120   format(7x,i3,5x,e10.4,5x,e10.4)
110    continue
      return
      end
      subroutine belt (nodem,n,nd)
      dimension  nodem(100,9)
      k=1
      write(*,110) n
110   format(10x,'n=', i10)
      do 120 I=1,n
      write(*,130)  (nodem(I,j),j=1,9)
130   format(9i7)
120   continue
      do 10 I=1,n
      m1=1
      do 20 j=1,8
      mj=nodem(I,j)
      if (mj.ge.m1) m1=mj
20    continue
      m2=10000
      do 30 j=1,8
      mj=nodem(i,j)
      if (mj.le.m2) m2=mj
30    continue
      ki=m1-m2
      if (ki.ge.k) k=ki
10    continue
      nd=(k+1)*2
      if (nd.le.110) goto 40
      write(*,50) nd
50    format(10x,'the belt number is larger than 110:',i5)
      stop
40    continue
      write(*,100) nd
100   format(20x,'nd=',i10)
      return
      end
        subroutine bound(a,kk2,n3,n4,amax,kb,f)
        dimension a(700,110),kk2(kb),f(40,2)
      do 10 I=1,kb
                  n=kk2(i)
           n1=2*n-1
                  n2=2*n
                  a(n1,1)=a(n1,1)*amax
            a(n2,1)=a(n2,1)*amax
            a(n1,n4)=a(n1,1)*f(i,1)
            a(n2,n4)=a(n2,1)*f(i,2)
10    continue
      return
      end