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.
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.
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.
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.
Figure 4: Dam with abutments
A fault is cut into the right abutment, and the blocks are zoned as shown in Figure 5.
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.
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
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Sep 26, 2024 |
