非饱和土边坡稳定性计算程序说明文件

一、数据文件

注意：以下文件内不要输入注释文字

1. dat文件

!dx  h2v
1.0  2.0
!nx1 nx2 ny1 ny2
20  20  10  5
!树荷载的左列和右列数，树荷载kN(各单元格相同)
27 32 10.0
!内摩擦角，粘聚力，剪胀角，重度，弹性模量，泊松比, a, m, n, theta_s 其中,重度=ρ_d*gam_w,theta_s=1-ρ_d/G_s
20.0  15.0  0.0  20.0  1.0e5  0.3 15 0.2 4 0.907749
!水位位置（坡顶为0，y方向为负值）
-20
!非饱和抗剪强度公式中κ
2
!收敛准则  计算次数
0.0001 500
!折减次数
6
!折减系数
1.2  1.3  1.4  1.5  1.55 1.6

其中,第三行表示荷载的加荷位置和大小

1
2
3

27 ! 表示起点在第27列（从左至右）
32 ! 表示终点在第32列（从左至右）
30 ! 表示荷载为30kN

对应下图中的红色位置

加载位置

2. 孔压文件x_f.txt

1
2
3

按节点列依次输入节点处的孔压
每个单元格有3列共8个节点，呈323分布，如下图所示，阴影部分为一个单元格，黄色为单元格号
按照节点的1234...顺序竖向填写x_f.txt文件，每节点占一行

17143576814081714357680903png

二、源代码

本程序单元格序号是从上至下，从左至右依次排序,如下例

1 6 
2 7  11
3 8  12 15
4 9  13 16 18 20
5 10 14 17 19 21

在叠加荷载时，因为所有单元格的g都是按照3 2 1 4 6 7 8 5 的顺序保存，因此我们仅需找到对应的单元格号以及146号节点的相对位置

146号节点的相对位置为g(6), g(8), g(10),下文有介绍

接下来只需找到单元格号即可。

计算单元号

if (tree_load/=0) then !荷载不为零再进行运算
  allocate(tree_nodes(nx1+nx2+ny1*ngrad))

  tree_nodes=0
   ! tree_nodes(nx1+nx2+ny1*ngrad)是储存表面单元格的单元号

  itree = 2 !用于计数，记录需加荷的单元格数量
  tree_nodes(1)=1 

  ! 顶面和斜面前ngrad行
  do i=2,nx1+ngrad 
     tree_nodes(itree) = (i-1)*nye+1
     itree = itree + 1
  end do

  ! 斜面ngrad行以后
  do i=1,ny1-1
    do j=1,ngrad
       if (j==1) then
          tree_nodes(itree) = tree_nodes(itree-1)+nye-i+1
       else
          tree_nodes(itree) = tree_nodes(itree-1)+nye-i
       end if
       itree = itree + 1
    end do
  end do

  ! 底面第一列
  tree_nodes(itree) = tree_nodes(itree-1)+ny2+1
  itree = itree + 1

  ! 底面以后
  do i=2,nx2
     tree_nodes(itree) = tree_nodes(itree-1)+ny2
     itree = itree + 1
  end do

  tree_nodes_= tree_nodes(nel_left:nel_right) ! 储存实际需要加荷的单元格的序号

分配荷载到节点

  ALLOCATE(weight_load(8),tree(ndof))
  !荷载在节点上的分布权重
  weight_load = (/0.166667,0.0,0.0,0.666667,0.0,0.166667,0.0,0.0/)

  !计算各节点荷载大小
  tree=zero
  tree(6)=tree_load*weight_load(1)
  tree(8)=tree_load*weight_load(4)
  tree(10)=tree_load*weight_load(6)

  ! 此处注意区分(见下图)
  ! g数组按3x 3y 2x 2y 1x 1y 4x 4y 6x 6y 7x 7y 8x 8y 5x 5y的顺序保存
  ! weight_load按1 2 3 4 5 6 7 8来保存
  ! 则1 4 6号节点的y方向对应g(6),g(8),g(10)

end if

单元节点

叠加荷载

itree=1 ! 作tree_nodes_数组的索引，将序号依次读出

elements_2: DO iel=1,nels

   ...

   gravlo(g)=gravlo(g)-eld*prop(4,etype(iel))

   if (tree_load.ne.0.0) then !荷载不为零时计算

      if ((iel==tree_nodes_(itree)).and.(itree.le.size(tree_nodes_))) then !判断是否为需加荷的单元

         gravlo(g) = gravlo(g) - tree !根据自由度编号加到指定位置

         itree=itree+1
      end if
   end if

END DO elements_2

叠加外部水荷载

叠加外部水荷载的方式与上方叠加外部荷载的方法类似，将水的压力 $\gamma_w*h$ 叠加到相应的表面节点上

IF(wl>g_coord(2,1))THEN  ! 水位淹没整个边坡
      p3=(wl-g_coord(2,1))*gam_w
      DO i=1,nx1
         gravlo(nf(2,g_num(3,tree_nodes(i))))=gravlo(nf(2,g_num(3,tree_nodes(i))))-dx*p3/d6
         gravlo(nf(2,g_num(4,tree_nodes(i))))=gravlo(nf(2,g_num(4,tree_nodes(i))))-dx*two*p3/d3
         gravlo(nf(2,g_num(5,tree_nodes(i))))=gravlo(nf(2,g_num(5,tree_nodes(i))))-dx*p3/d6
      END DO
      DO i=nx1+1,nx1+ny1*ngrad
         if (mod(i-nx1,ngrad)/=0)then
            y5=g_coord(1,g_num(5,tree_nodes(i)))
            p5=(wl-y5)*gam_w
            gravlo(nf(2,g_num(3,tree_nodes(i))))=gravlo(nf(2,g_num(3,tree_nodes(i))))-dx*p5/d6
            gravlo(nf(2,g_num(4,tree_nodes(i))))=gravlo(nf(2,g_num(4,tree_nodes(i))))-dx*two*p5/d3
            gravlo(nf(2,g_num(5,tree_nodes(i))))=gravlo(nf(2,g_num(5,tree_nodes(i))))-dx*p5/d6
         else
            y5=g_coord(2,g_num(5,tree_nodes(i)))
            p5=(wl-y5)*gam_w
            y7=g_coord(2,g_num(7,tree_nodes(i)))
            p7=(wl-y7)*gam_w
            gravlo(nf(2,g_num(3,tree_nodes(i))))=gravlo(nf(2,g_num(3,tree_nodes(i))))-dx*p5/d6
            gravlo(nf(2,g_num(4,tree_nodes(i))))=gravlo(nf(2,g_num(4,tree_nodes(i))))-dx*two*p5/d3
            gravlo(nf(2,g_num(5,tree_nodes(i))))=gravlo(nf(2,g_num(5,tree_nodes(i))))-dx*p5/d6
            gravlo(nf(1,g_num(5,tree_nodes(i))))=gravlo(nf(1,g_num(5,tree_nodes(i))))-dx*p5/d6
            gravlo(nf(1,g_num(6,tree_nodes(i))))=gravlo(nf(1,g_num(6,tree_nodes(i))))-dx*(p5+p7)/d3
            gravlo(nf(1,g_num(7,tree_nodes(i))))=gravlo(nf(1,g_num(7,tree_nodes(i))))-dx*p7/d6
         end if
      END DO

      IF(ny2>0)THEN ! 如果有破角台阶
         p3=(wl-g_coord(2,nn-2*ny2))*gam_w !p3==p5
         DO i=nx1+ny1*ngrad+1,nx1+ny1*ngrad+nx2
            gravlo(nf(2,g_num(3,tree_nodes(i))))=gravlo(nf(2,g_num(3,tree_nodes(i))))-dx*p3/d6
            gravlo(nf(2,g_num(4,tree_nodes(i))))=gravlo(nf(2,g_num(4,tree_nodes(i))))-dx*two*p3/d3
            gravlo(nf(2,g_num(5,tree_nodes(i))))=gravlo(nf(2,g_num(5,tree_nodes(i))))-dx*p3/d6
         END DO
      END IF
   ELSE IF(wl>g_coord(2,nn-2*ny2))THEN ! 水面在斜面位置(低于坡顶,高于坡脚)
      nw=int((-wl)/dx)
      do i=nx1+nw*ngrad+ngrad,nx1+ny1*ngrad
         y5=g_coord(2,g_num(5,tree_nodes(i)))
         y7=g_coord(2,g_num(7,tree_nodes(i)))
         p5=(wl-y5)*gam_w
         p7=(wl-y7)*gam_w
         if(i==nx1+nw*ngrad+ngrad)then ! 如果水面恰好穿过某行单元格
            cx=y5-wl
            cy=cx
            gravlo(nf(1,g_num(5,tree_nodes(i))))=gravlo(nf(1,g_num(5,tree_nodes(i))))           &
               +p7*cy*(cy**2-two*cy*dx+dx**2)/(d6*dx**2)
            gravlo(nf(1,g_num(6,tree_nodes(i))))=gravlo(nf(1,g_num(6,tree_nodes(i))))           &
               -p7*(cy**3-cy**2*dx-cy*dx**2+dx**3)/(d3*dx**2)
            gravlo(nf(1,g_num(7,tree_nodes(i))))=gravlo(nf(1,g_num(7,tree_nodes(i))))           &
               -p7*(dx**3-cy**3)/(d6*dx*dx)
         else if(mod(i-nx1,ngrad)/=0)then ! 方格斜面台阶内部(仅有上表面接触水)
            gravlo(nf(2,g_num(3,tree_nodes(i))))=gravlo(nf(2,g_num(3,tree_nodes(i))))-dx*p5/d6
            gravlo(nf(2,g_num(4,tree_nodes(i))))=gravlo(nf(2,g_num(4,tree_nodes(i))))-dx*two*p5/d3
            gravlo(nf(2,g_num(5,tree_nodes(i))))=gravlo(nf(2,g_num(5,tree_nodes(i))))-dx*p5/d6
         else ! 方格斜面台阶最外部的单元格(上表面和右表面接触水)
            gravlo(nf(2,g_num(3,tree_nodes(i))))=gravlo(nf(2,g_num(3,tree_nodes(i))))-dx*p5/d6
            gravlo(nf(2,g_num(4,tree_nodes(i))))=gravlo(nf(2,g_num(4,tree_nodes(i))))-dx*two*p5/d3
            gravlo(nf(2,g_num(5,tree_nodes(i))))=gravlo(nf(2,g_num(5,tree_nodes(i))))-dx*p5/d6
            gravlo(nf(1,g_num(5,tree_nodes(i))))=gravlo(nf(1,g_num(5,tree_nodes(i))))-dx*p5/d6
            gravlo(nf(1,g_num(6,tree_nodes(i))))=gravlo(nf(1,g_num(6,tree_nodes(i))))-dx*(p5+p7)/d3
            gravlo(nf(1,g_num(7,tree_nodes(i))))=gravlo(nf(1,g_num(7,tree_nodes(i))))-dx*p7/d6
         end if
      END DO
      IF(ny2>0)THEN 
         p3=(wl-g_coord(2,nn-2*ny2))*gam_w
         DO i=nx1+ny1*ngrad+1,nx1+ny1*ngrad+nx2
            gravlo(nf(2,g_num(3,tree_nodes(i))))=gravlo(nf(2,g_num(3,tree_nodes(i))))-dx*p3/d6
            gravlo(nf(2,g_num(4,tree_nodes(i))))=gravlo(nf(2,g_num(4,tree_nodes(i))))-dx*two*p3/d3
            gravlo(nf(2,g_num(5,tree_nodes(i))))=gravlo(nf(2,g_num(5,tree_nodes(i))))-dx*p3/d6
         END DO
      END IF
   END IF

叠加孔压

根据Vanapalli公式, 使用体积含水率、孔压等参数计算出基质吸力项, 叠加到正应力上

17143599364091714359936244png

elements_3: DO iel=1,nels
            bload=zero
            phi=prop(1,iel)
            tnph=TAN(phi*pi/d180)
            phif=ATAN(tnph/srf(iy))*d180/pi
            psi=prop(3,iel)
            tnps=TAN(psi*pi/d180)
            psif=ATAN(tnps/srf(iy))*d180/pi
            cf=prop(2,iel)/srf(iy)
            e=prop(5,iel)
            v=prop(6,iel)
            CALL deemat(dee,e,v)
            num=g_num(:,iel)
            coord=TRANSPOSE(g_coord(:,num))
            do i=1,nod
               call Vanapalli(i,iel,nprops,kapa,prop,theta_w,suc_sig,n_inpore) ! 计算基质吸力项
            end do
            g=g_g(:,iel)
            eld=loads(g)
            gauss_pts_2: DO i=1,nip
               call shape_fun(fun,points,i)
               pore(iel,i)=dot_product(fun,suc_sig)
               CALL bee8(bee,coord,points(i,1),points(i,2),det)
               eps=MATMUL(bee,eld)
               eps=eps-evpt(:,i,iel)
               sigma=MATMUL(dee,eps)
               sigma(1)=sigma(1)-pore(iel,i) ! 此处124为正应力项, 3为剪应力项
               sigma(2)=sigma(2)-pore(iel,i)
               sigma(4)=sigma(4)-pore(iel,i)
               CALL invar(sigma,sigm,dsbar,lode_theta)
               !-----------------------check whether yield is violated-------------------
               CALL mocouf(phif,cf,sigm,dsbar,lode_theta,f)
               IF(f>fmax)fmax=f
               IF(converged.OR.iters==limit)THEN
                  sigma(1)=sigma(1)+pore(iel,i)
                  sigma(2)=sigma(2)+pore(iel,i)
                  sigma(4)=sigma(4)+pore(iel,i)
                  devp=sigma
               ELSE
                  IF(f>=zero.OR.(converged.OR.iters==limit))THEN
                     CALL mocouq(psif,dsbar,lode_theta,dq1,dq2,dq3)
                     CALL formm(sigma,m1,m2,m3)
                     flow=f*(m1*dq1+m2*dq2+m3*dq3)
                     erate=MATMUL(flow,sigma)
                     evp=erate*dt
                     evpt(:,i,iel)=evpt(:,i,iel)+evp
                     devp=MATMUL(dee,evp)
                  END IF
               END IF
               IF(f>=zero)THEN
                  eload=MATMUL(devp,bee)
                  bload=bload+eload*det*weights(i)
               END IF
            END DO gauss_pts_2
            !-----------------------compute the total bodyloads vector----------------
            bdylds(g)=bdylds(g)+bload
            bdylds(0)=zero
         END DO elements_3