Finite Element Dam

Problem Statement

Note

This example demonstrates one technique that may be used to create finite element (FE) blocks in 3DEC. In this case, the finite elements compose a dam with a base and abutments.

The project file for this example may be viewed/run in 3DEC.[1] The data files used are shown at the end of this example.

3DEC Model Geometry

The model building starts with the construction of the dam segments as five ordinary 3DEC blocks. This is done in a FISH function (geodam) using the block create prism command. The FISH function also sets up FISH parameters that are used for dimensions in later commands. This geometry is shown in Figure 1.

../../../../_images/dampanels.png

Figure 1: Dam panels

Each of the 3DEC blocks contains eight vertices. After the geometry has been created as 3DEC blocks, the blocks are then converted to finite element blocks using the command feblock generate.

Each of the top central blocks have been converted into three finite element bricks with 20 nodes each. The bottom central blocks have each been converted into one degenerate finite element brick. The end blocks each contain two degenerate finite element bricks.

The next step is to change the face group names on the dam blocks. This is done using a FISH function (markdamfaces) to set the face groups for the block faces. The block face groups are then used by the block generate from-topo command to extrude the existing blocks down to the base of the model. The extruded blocks for this step are shown in Figure 2.

../../../../_images/damextrude1.png

Figure 2: Extrusion of panels

The base blocks are then extruded in the y-direction to form the upstream and downstream base blocks as shown in Figure 3.

../../../../_images/damextrude2.png

Figure 3: Blocks extruded in the upstream and downstream directions

The upstream and downstream abutments are then created using block create prism commands. The result of these commands is shown in Figure 4.

../../../../_images/damabut.png

Figure 4: Dam with abutments

A fault is cut into the right abutment, and the blocks are zoned as shown in Figure 5.

../../../../_images/damfebz.png

Figure 5: Zoned model

Dam Loading

The assignment of properties and the loading of the dam is shown in the file dam-loading.dat.

The next step is to assign properties to the rock, dam and joints. Initially, the joints are kept elastic to prevent movement during the cycling to obtain gravitational equilibrium. The dam is also given an artificially low modulus.

Gravity is defined, and in-situ stresses are applied to the rock and the dam. Then the dam is hidden, and the block insitu topograph command is used to initialize the valley stresses. The stresses are then removed from the dam.

At this point we want to run the model to equilibrium. We can use the feblock gravity delete command to remove the gravitational load from the dam for the finite elements. We do not want the dam loading the valley yet. We also apply the boundary conditions to the model and assign some monitoring histories.

Now we want to remove any stress from the dam that may have been generated by the valley floor and walls, and put the gravity back in the finite element model. We also reset the monitoring histories.

In this stage we use a boundary condition to add water pressures to the dam (block face apply stress).

After cycling 5000 steps we piut the joints back to the proper strength values and use a FISH function (ppjoint) to add pore pressures to the joints.

Finally we reduce the friction on the fault to let it slip. We also remove the cohesion between the rock and the concrete to simulate a cracked contact. The final result showing the slipping fault is shown in Figure 6.

../../../../_images/damdisp.png

Figure 6: Displacements in rock and dam

Data Files

Note

For clarity, all commands and keywords are written in their entirey and do not use command abbreviation when presented in example data files in documentation.

dam-geometry.dat

model new
; ---------------------------------------------------------------------
; --- FE block example - Arch dam analysis ---
; ---------------------------------------------------------------------
model configure feblock
model configure dynamic
;
block tolerance .001
;
; --- GEOMDAM ---
; FISH function to define dam geometry
; (simplified 5 block cylindrical arch shape)
;
fish define geomdam

  rad = 100.0
  xcent = 0.0
  ycent = -rad
  thick = 10.0
  angtot = 100.0 * math.degrad
  nb = 5
  angb = angtot / nb
  zz2 = 80.0
  zz1 = 0.0

  loop i (1,nb)

     ang1 = -0.5*angtot + (i-1)*angb
     ang2 = ang1 + angb
     xx1 = xcent + (rad-thick)*math.sin(ang1)
     yy1 = ycent + (rad-thick)*math.cos(ang1)
     xx2 = xcent + (rad-thick)*math.sin(ang2)
     yy2 = Ycent + (rad-thick)*math.cos(ang2)
     xx3 = xcent + (rad)*math.sin(ang2)
     yy3 = ycent + (rad)*math.cos(ang2)
     xx4 = xcent + (rad)*math.sin(ang1)
     yy4 = ycent + (rad)*math.cos(ang1)

     command
       block create prism face-1 [xx1] [yy1] [zz1] [xx2] [yy2] [zz1] ...
                                 [xx3] [yy3] [zz1] [xx4] [yy4] [zz1] ...
                          face-2 [xx1] [yy1] [zz2] [xx2] [yy2] [zz2] ...
                                 [xx3] [yy3] [zz2] [xx4] [yy4] [zz2] ...
                          group ['block='+string(i)]
    endcommand

     if i = 1
       xa1 = xx2
       ya1 = yy2
       xb1 = xx3
       yb1 = yy3
       xd1 = xx1
       yd1 = yy1
     endif
     if i = 5
       xa2 = xx1
       ya2 = yy1
       xb2 = xx4
       yb2 = yy4
       xd2 = xx2
       yd2 = yy2
     endif
     if i = 2
       xe1 = xx2
       ye1 = yy2
       xf1 = xx3
       yf1 = yy3
     endif
     if i = 4
       xe2 = xx1
       ye2 = yy1
       xf2 = xx4
       yf2 = yy4
     endif

; abutment points
     if i = 1
       xxright = xx4
       yyright = yy4
       xxrightdown = xx1
       yyrightdown = yy1
       xxrightd = xx1
       yyrightd = yy1
     endif
     if i = 5
       xxleft = xx3
       yyleft = yy3
       xxleftdown = xx2
       yyleftdown = yy2
       xxleftd = xx2
       yyleftd = yy2
     endif

; central cantilever
     if i = 3
       xxrcen = xx1
       yyrcen = yy1
       xxlcen = xx2
       yylcen = yy2
     endif

  endloop

  za = 20.0
  xc1 = -20.0
  yc1 = -thick
  zd = 60.0
  zd07 = 0.8*zd
  xc2 = 20.0
  yc2 = -thick
  command
; horizontal cut
    block hide range group 'block=1' 
    block hide range group 'block=5' 
    block cut joint-set points [xa1] [ya1] [za] ...
                               [xb1] [yb1] [za] ...
                               [xc1] [yc1] [za]
    block join range group 'block=3'
;   

; right abutment
    block hide range group 'block=3'
    block hide range group 'block=4'
    block delete range position-z [zz1] [za]
    block create polyhedron face [xe1] [ye1] [zz1] ...
                                 [xa1] [ya1] [za] ...
                                 [xe1] [ye1] [za] ...
                            face [xf1] [yf1] [za] ...
                                 [xb1] [yb1] [za] ...
                                 [xf1] [yf1] [zz1] ...
                            face [xb1] [yb1] [za] ...
                                 [xa1] [ya1] [za] ...
                                 [xe1] [ye1] [zz1] ...
                                 [xf1] [yf1] [zz1] ...
                            face [xa1] [ya1] [za] ...
                                 [xb1] [yb1] [za] ...
                                 [xf1] [yf1] [za] ...
                                 [xe1] [ye1] [za] ...
                            face [xe1] [ye1] [za] ...
                                 [xf1] [yf1] [za] ...
                                 [xf1] [yf1] [zz1] ...
                                 [xe1] [ye1] [zz1] ...
                            group 'block=2'   
    block join range group 'block=2'
    block hide off
                                  

    block hide range group 'block=1'  not
    block cut joint-set points [xa1] [ya1] [za] ...
                               [xb1] [yb1] [za] ...
                               [xd1] [yd1] [zd]
    block delete range position-z [zz1] [zd07]
    block hide off
;

; left abutment
    block hide range group 'block=4' not
    block delete range position-z [zz1] [za]
    block create polyhedron face [xa2] [ya2] [za] ...
                                 [xe2] [ye2] [zz1] ...
                                 [xe2] [ye2] [za] ...
                            face [xb2] [yb2] [za] ...
                                 [xf2] [yf2] [za] ...
                                 [xf2] [yf2] [zz1] ...
                            face [xa2] [ya2] [za] ...
                                 [xb2] [yb2] [za] ...
                                 [xf2] [yf2] [zz1] ...
                                 [xe2] [ye2] [zz1] ...
                            face [xb2] [yb2] [za] ...
                                 [xa2] [ya2] [za] ...
                                 [xe2] [ye2] [za] ...
                                 [xf2] [yf2] [za] ...
                            face [xf2] [yf2] [za] ...
                                 [xe2] [ye2] [za] ...
                                 [xe2] [ye2] [zz1] ...
                                 [xf2] [yf2] [zz1] ...
                            group 'block=4'
    block join range  group 'block=4'
    block hide off;
    
    block hide range group 'block=5'  not    
    block cut joint-set points [xa2] [ya2] [za] ...
                               [xb2] [yb2] [za] ...
                               [xd2] [yd2] [zd]
    block delete range position-z [zz1] [zd07]
    block hide off
 endcommand
  xx45 = xa2
  yy45 = ya2
  zz45 = za
  xx45d = xb2
  yy45d = yb2
  zz45d = za

  izzright = block.gp.near(xxright,yyright,64.0)
  izzleft = block.gp.near(xxleft,yyleft,64.0)
  zzright = block.gp.pos.z(izzright)
  zzleft = block.gp.pos.z(izzleft)
  ii=io.out(' zzright '+string(zzright))
  ii=io.out(' zzleft '+string(zzleft))

; model boundaries
  xxrbou = -150.0
  xxlbou = 150.0
  yyubou = 150.0
  yydbou = -180.0
  zztopbou = 80.0
  zzbotbou = -100.0

end
[geomdam]
block group 'dam'

model save 'dampanels'
;
; create FE mesh in dam blocks
;
feblock generate 20 20 20  range group 'block=2' 
feblock generate 20 20 20  range group 'block=3' 
feblock generate 20 20 20  range group 'block=4' 
feblock generate 20 20 13  range group 'block=1' 
feblock generate 20 20 13  range group 'block=5' 
;
; --- MARKDAMFACES ---
; mark dam faces (fe faces and zone faces)
fish define markdamfaces
; foundation joint
  ii=io.out(' found. joint')
  loop foreach ib block.list
    section
    if block.group(ib,'block') # "1" & block.group(ib,'block') # "5"
        if block.pos.z(ib) > za
          exit section
        endif
    endif
; fe faces
    loop foreach ifef block.feb.facelist(ib) 
      if feblock.face.normal.z(ifef) < -0.1
        feblock.face.group(ifef) = 'face_'+block.group(ib,'block')
        feblock.face.region(ifef) = int(block.group(ib,'block'))
      endif
    endloop
; zone faces
    loop foreach ifa block.facelist(ib,false) 
      if block.face.normal.z(ifa) < -0.1
        block.face.group(ifa) = 'face_'+block.group(ib,'block')
      endif
    endloop
    endsection
  endloop
; upstream face
  ii=io.out(' upstream face')
  loop foreach ib block.list
    section
; fe faces
    loop foreach ifef block.feb.facelist(ib)
      vnx = feblock.face.normal.x(ifef)
      vny = feblock.face.normal.y(ifef)
      vnz = feblock.face.normal.z(ifef)
      if math.abs(vnz) < 0.1
        angxy = math.atan2(vnx,vny) / math.degrad
        if math.abs(angxy) < 41.0
            feblock.face.group(ifef) = "upstreamface"
        endif
      endif
    endloop
; zone faces
    loop foreach ifa block.facelist(ib)
      vnx = block.face.normal.x(ifa)
      vny = block.face.normal.y(ifa)
      vnz = block.face.normal.z(ifa)
      if math.abs(vnz) < 0.1
        angxy = math.atan2(vnx,vny) / math.degrad
        if math.abs(angxy) < 41.0
            block.face.group(ifa) = "upstreamface"
        endif
      endif
    endloop
    endsection
  endloop
end
[markdamfaces]

;
; --- create rock blocks ---
; below dam
block generate from-topography direction 0 0 -1 plane origin 0 0 [zzbotbou] ...
  normal 0 0 1 group 'downstream' range group 'face_1'
block generate from-topography direction 0 0 -1 plane origin 0 0 [zzbotbou] ...
  normal 0 0 1 group 'downstream' range group 'face_2'
block generate from-topography direction 0 0 -1 plane origin 0 0 [zzbotbou] ...
  normal 0 0 1 group 'downstream' range group 'face_3'
block generate from-topography direction 0 0 -1 plane origin 0 0 [zzbotbou] ...
  normal 0 0 1 group 'downstream' range group 'face_4'
block generate from-topography direction 0 0 -1 plane origin 0 0 [zzbotbou] ...
  normal 0 0 1 group 'downstream' range group 'face_5'

model save 'damextrude1'

;
block hide range group "dam"
block generate from-topography direction 0 1 0  plane origin 0 [yyubou] 0 ...
      normal 0 1 0 group 'upstream' range cylinder end-1 [xcent] [ycent]  ...
      -100 end-2 [xcent] [ycent] 100 radius [rad-0.1*thick] [rad-1.00*thick]
block generate from-topography direction 0 -1 0 plane origin 0 [yydbou] 0  ...
      normal 0 1 0 group 'downstream' range cylinder end-1 [xcent] [ycent] ...
      -100 end-2 [xcent] [ycent] 100 radius [rad-1.25*thick] [rad-1.00*thick]
block hide off

model save 'damextrude2'
;
; abutments downstream
block create prism  face-1 [xxright] [yyright] [zzright] ...
                           [xxrightd] [yyrightd] [zd] ...
                           [xxrightd] [yyrightd] [zztopbou] ...
                           [xxright] [yyright] [zztopbou] ...
                    face-2 [xxright] [yydbou] [zzright] ...
                           [xxrightd] [yydbou] [zd] ...
                           [xxrightd] [yydbou] [zztopbou] ...
                           [xxright] [yydbou] [zztopbou] ...
                    group 'downstream'
block create prism  face-1 [xxleft] [yyleft] [zzleft] ...
                           [xxleftd] [yyleftd] [zd] ...
                           [xxleftd] [yyleftd] [zztopbou] ...
                           [xxleft] [yyleft] [zztopbou] ...
                    face-2 [xxleft] [yydbou] [zzleft] ...
                           [xxleftd] [yydbou] [zd] ...
                           [xxleftd] [yydbou] [zztopbou] ...
                           [xxleft] [yydbou] [zztopbou] ...
                    group 'downstream'
;
; abutments downstream, below foundation joint
block create prism  face-1 [xxright] [yyright] [zzright] ...
                           [xxrightd] [yyrightd] [zd] ...
                           [xxrightd] [yyrightd] [zzbotbou] ...
                           [xxright] [yyright] [zzbotbou] ...
                    face-2 [xxright] [yydbou] [zzright] ...
                           [xxrightd] [yydbou] [zd] ...
                           [xxrightd] [yydbou] [zzbotbou] ...
                           [xxright] [yydbou] [zzbotbou] ...
                    group 'downstream'
block create prism  face-1 [xxleft] [yyleft] [zzleft] ...
                           [xxleftd] [yyleftd] [zd] ...
                           [xxleftd] [yyleftd] [zzbotbou] ...
                           [xxleft] [yyleft] [zzbotbou] ...
                    face-2 [xxleft] [yydbou] [zzleft] ...
                           [xxleftd] [yydbou] [zd] ...
                           [xxleftd] [yydbou] [zzbotbou] ...
                           [xxleft] [yydbou] [zzbotbou] ...
                    group 'downstream'
                    
;
; left and right boundary blocks
; upstream
block create brick [xxrbou] [xxright] [yyright] ...
                   [yyubou] [zzbotbou] [zztopbou] ...
                   group 'upstream'
block create brick [xxleft] [xxlbou] [yyright] ...
                   [yyubou] [zzbotbou] [zztopbou] ...
                   group 'upstream'
; downstream
block create brick [xxrbou] [xxright] [yydbou] ...
                   [yyright] [zzbotbou] [zztopbou] ...
                   group 'downstream'
block create brick [xxleft] [xxlbou] [yydbou] ...
                   [yyright] [zzbotbou] [zztopbou] ...
                   group 'downstream'
;
model save 'damabut'
;
; join rock blocks
block join range group 'upstream'
block join range group 'downstream'
;
; create wedgelength in left abutment
block hide range group 'downstream' not
block cut joint-set dip 90 dip-direction -80 origin [xxleft] [yyleft] 80
block cut joint-set dip 90 dip-direction -65 origin 70 -60 80
block cut joint-set dip 0 origin [xx45] [yy45] [zz45]
block hide off
;
block join range group 'block_2'
block join range group 'block_4'
block join range group 'block_5'
;
; mesh rock blocks
block zone generate edgelength 50
;
block contact group "dam"  range group-intersection 'dam' 'dam'
block contact group "foundation" range group-intersection 'dam' 'downstream'
block contact group "fault_1" range plane dip 90 dip-direction -80 ...
                                    origin [xxleft] [yyleft] 80 distance 0.001
block contact group "fault_2" range plane dip 90 dip-direction -65 ...
                                    origin 70 -60 80 distance 0.001
block contact group "fault_3" range plane dip 0 ...
                                    origin [xx45] [yy45] [zz45]  distance 0.001
;
model save 'damfebz'
;

dam-loading.dat

model restore 'damfebz'
;
; --- assign materials, properties ---
;
block zone cmodel assign elastic
;
; rock
; E=20000 MPa, v=0.2, K=11111, G=8333
block zone property density 0.0027 bulk 11111 shear 8333
;
; concrete
; use E/100 forigin insitu stage
; E=200 MPa, v=0.2, K=111.11, G=83.33
block zone property density 0.0024 bulk 111.11 shear 83.33  range group 'dam'
;
; rock joints
block contact jmodel assign mohr
;
block contact property stiffness-normal 10000 stiffness-shear 5000 ...
  friction 38 range group 'Default=None'
;
; dam joints (elastic forigin now)
block contact property stiffness-normal 10000 stiffness-shear 5000 ...
  cohesion 1e20 tension 1e20 range group 'dam'
;
; foundation joint (elastic forigin now)
block contact property stiffness-normal 10000 stiffness-shear 5000 ...
  cohesion 1e20 tension 1e20 range group 'foundation'
;
; --- stage 1 : insitu stresses ---
; (dam blocks soft and elastic)
;
; gravity
model gravity 0 0 -10
;
block insitu stress -1.0e-4 -1.0e-4 -1.0e-4 0 0 0
;
; hide dam
block hide range group 'dam'
;
; topographic stresses
block insitu topography ratio-x 0.3 ratio-y 0.3
;
; show dam
block hide off
;
model save 'damfebs0'
;
model large-strain off
; remove dam stresses
block zone initialize stress 0 0 0 0 0 0 range group 'dam'
block contact reset stress range group 'dam'
block contact reset stress range group 'foundation'
;
; remove dam weight
feblock gravity delete
;
; boundary
block gridpoint apply velocity-y 0. range position-y -180
block gridpoint apply velocity-y 0. range position-y 150
block gridpoint apply velocity-x 0. range position-x -150
block gridpoint apply velocity-x 0. range position-x 150
block gridpoint apply velocity-z 0. range position-z -100
;
; hist
model history mechanical unbalanced-maximum
block history displacement-x position 0 0 80
block history displacement-y position 0 0 80
block history displacement-z position 0 0 80
block hide range group 'dam'
block history displacement-x position 0 0 0
block history displacement-y position 0 0 0
block history displacement-z position 0 0 0
block hide off
;
model dynamic active off
block mechanical damping local
model solve
model save 'damfebs1'
;

; --- stage 2 : dam self weight ---
; (real dam properties, elastic joints)
model restore 'damfebs1'
;
block gridpoint initialize displacement 0 0 0
block gridpoint initialize velocity 0 0 0
block contact reset displacement
;
; remove dam stresses
block zone initialize stress 0 0 0 0 0 0 range group 'dam'
block contact reset stress range group 'dam'
block contact reset stress range group 'foundation'
;
; concrete properties (real E)
; E=20000 MPa, v=0.2, K=11111, G=8333
block zone property bulk 11111 shear 8333 range group 'dam'
;
; hist
history delete
model history mechanical unbalanced-maximum
block history displacement-x position 0 0 80
block history displacement-y position 0 0 80
block history displacement-z position 0 0 80
block hide range group 'dam'
block history displacement-x position 0 0 0
block history displacement-y position 0 0 0
block history displacement-z position 0 0 0
block hide off
block history displacement-x position 70 -40 66
block history displacement-y position 70 -40 66
block history displacement-z position 70 -40 66
block history displacement-x position 57 -51 40
block history displacement-y position 57 -51 40
block history displacement-z position 57 -51 40
;
; dam self weight
feblock gravity
model solve
;
model save 'damfebs2'
;
; --- stage 3 - water load ---
; (elastic joints)
model restore 'damfebs2'
;
; boundary, fix base
block gridpoint apply velocity-x 0 velocity-y 0 velocity-z 0 ...
  range position-z -100
;
; water load
block face apply stress -0.8 -0.8 -0.8 0 0 0 ...
  gradient-z 0.01 0.01 0.01 0 0 0 range group 'upstreamface'
;
model solve
;
model save 'damfebs3'
; --- stage 4 - water load ---
; (nonlinear dam & foundation joints)
;
block contact property stiffness-normal 10000 stiffness-shear 5000 ...
  cohesion 0 tension 0 friction 45 range group 'dam'
block contact property stiffness-normal 10000 stiffness-shear 5000 ...
  cohesion 3 tension 2 friction 40 range group 'foundation'
;
model solve
model save 'damfebs4'
;
; --- stage 5 : pp in rock joints ---
model restore 'damfebs4'
;
fish define ppjoint
  loop foreach ic block.contact.list
  ; skip dam and foundation joints
    b1 = block.contact.b1(ic)
    b2 = block.contact.b2(ic)
    if block.group(b1) = 'dam'| block.group(b2) = 'dam'
      continue
    endif
    loop foreach icx block.contact.subcontactlist(ic)
      xx = block.subcontact.pos.x(icx)
      yy = block.subcontact.pos.y(icx)
      zz = block.subcontact.pos.z(icx)
      ar = block.subcontact.area(icx)
      iver = block.subcontact.gp(icx)
      ppold = block.subcontact.pp(icx)
    ;pp only in joints in vicinity of dam
      if xx < -82.1 | xx > 82.1
        continue
      endif
      if yy < -60 |yy > 0.1 
        continue
      endif
      ppnew = (80-zz)*0.01
      fpnew = math.max(ppnew*ar,0.0)
      ; set pp force in contact
      block.subcontact.force.pp(icx) = fpnew
    endloop
  endloop
end
[ppjoint]
;
;
block update subcontact 0
;
model solve
model save 'damfebs5'
;
; --- stage 6 - reduce friction on rock joints ---
model restore 'damfebs5'
;
block contact property friction 20 range group 'None'
;
block contact property friction 20 range group 'fault_1'
block contact property friction 20 range group 'fault_2'
; assume foundation cracked
block contact property cohesion 0 tension 0 range group 'foundation'
;
model cycle 5000
model save 'damfebs6'
;
program return

Endnote