WANtaroHP (FEM: 2D Seepage Flow Analysis)

Contents

• Outline of a program 'f90_fem_seep.f90'
• Format of input file ('csv' format)
• Format of output file ('csv' format)
• Pressure head analysis around the tunnel
• Zoned embankment dam
• Homogeneous embankment dam
• Programs and sample data
• Expression of result of seepage flow analysis (f90_gmt_seep.f90)

Outline of a program 'f90_fem_seep.f90'

• This is a program for 2D Steady Seepage Flow Analysis. This program can not treat unsteady problem.
• This program can do the saturated-unsaturated seepage flow analysis in the vertical section, and saturated seepage flow analysis in the horizontal and vertical section.
• Von Genuchten model is used as a unsaturated ground model.
• In this program, 3 node element or 4 node element can be used. Mixed model of 3 node elements and 4 node elements can not be solved. User of this program needs to choose either 3 node model or 4 node model.
• Linear triangular element is used as a 3 node element. 1 node has 2 degrees of freedom in horizontal direction and vertical direction.
• 4 node isoparametric element with 4 Gauss points is used as a 4 node element. 1 node has 2 degrees of freedom in horizontal direction and vertical direction.
• Thickness of the element is unit thickness and it is not changeable.
• As boundary conditions, impervious boundary, given total head boundary, given discharge boundary and saturation point boundary can be considered.
• In coordinate system, x-direction is defined as right direction and z-direction is defined as upward direction.
• Only an isotropic material can be treated.
• Nodal head must be inputted as total head and solution of simultaneous equations is described as total head.
• Pressure head is calculated as the difference between total head and position head.
• Simultaneous linear equations are solved using Cholesky method for banded matrix.
• Convergence criterion is that the case change of pressure head becomes less than 1$\times$10$^{-6}$. Upper limit of iteration is set to 100 times.
• Input/Output file name can be defined arbitrarily with the format of 'csv' and those are inputted from command line of a terminal.
• Used language for program is 'Fortran 90' and used compiler is 'GNU gfortran.'

Format of input file ('csv' format)

Comment                              # Comment
nod,NODT,NELT,MATEL,KOH,KOQ,KOU,idan # basic values for analysis
Ak0,alpha,em                         # Material properties
    ....(1 to MATEL)....             #
node-1,node-2,node-3(,node-4),matno  # Element-nodes relationship, material set number
    ....(1 to MATEL)....             #  (counterclockwise order of node numbers)
x,z,hvec                             # Node coordinate, initial value of total head
    ....(1 to NODT)....              # (x: right-direction, z: upward-direction)
nokh,Hinp                            # Node number of given total head boundary, total head of node
    ....(1 to KOH)....               # (Omit data input if KOH=0)
nokq,Qinp                            # Node number of given discharge boundary, discharge of node
    ....(1 to KOQ)....               # (Omit data input if KOQ=0)
noku                                 # Node number of saturation point
    ....(1 to KOU)....               # (Omit data input if KOU=0)

Notice

• Unsaturated properties are $\alpha$ and $m$ on van Genuchten model.
  \begin{align*} &S_e=\left\{\frac{1}{1+\left(\alpha\cdot h_s\right)^n}\right\}^m \qquad\qquad n=\frac{1}{1-m} \quad & (0 < m < 1) \\ &K_r=(S_e)^{0.5}\cdot\left\{1-\left(1-S_e{}^{1/m}\right)^m\right\}^2 & (0 \leqq S_r, K_r \leqq 1) \\ &K=K_r \cdot K_0 \end{align*}
  

  $K$	: Permeability coefficient	$S_e$	: Degree of saturation
  $K_0$	: Saturated permeability coefficient	$h_s$	: Suction head
  $K_r$	: Relative hydraulic conductivity function	$\alpha$	: Scaling parameter
  		$m$	: Non-dimensional parameter

  
  

  Relationship between Degree of saturation and Sunction head

  
  

  Relationship between Degree of saturation and Relative hydraulic conductivityg

  
  
• Regarding the values of $\alpha$ and $m$, Webmaster usually use $\alpha=0.1$ and $m=0.7$.
• Condition of node at which noting is inputted is no discharge. (in-flow discharge + out-flow discharge = 0) If input of in-flow or out-flow discharge is wanted especially, the values of KOQ, nokq and Qinp should be inputted.
• The discharge is known if total head is unknown at a node, and total head is known if discharge is unknown at a node. This condition is maintained except the node with saturation point.
• Nodes which have the possibility becoming saturation point such as downstream slope of embankment dam shall be specified as saturation point boundary with input of KOU and noku.
• It is necessary to input initial total heads for all nodes for iterative process of saturates-unsaturated analysis. It is better to input minimum total head in the domain. If proper initial total heads are not inputted, unsuitable results may be obtained.

Format of output file ('csv' format)

Comment
nod,NODT,NELT,MATEL,KOH,KOQ,KOU,idan
(Each value of above items)
*node characteristics
node,x,z,KOH,KOQ,KOU,Hinp,Qinp
    node : Node number
    x,z  : (x,z) Coordinate of node
    KOH  : Node with given total head (1: yes, 0: no)
    KOQ  : Node with given discharge (1: yes, 0: no)
    KOU  : Node with saturation point (1: yes, 0: no)
    Hinp : Given total head (if KOH=0, value=0)
    Qinp : Given discharge (if KOQ=0, value=0)
  .....(1 to NODT).....
*element characteristics
element,node-1,node-2,node-3(,node-4),Ak0,alpha,em,matno
    element : Element number
    node-1,node-2,node-3(,node-4) : Element-nodes relationship
    Ak0     : Permeability coefficient of element (saturated)
    alpha   : Unsaturated property $\alpha$
    em      : Unsaturated property $m$
    matno   : material set number
  .....(1 to NELT).....
*node value
node,x,z,hvec,qvec,pvec
    node : Node number
    x    : x-coordinate of node
    z    : z-coordinate of node
    hvec : Total head of node
    qvec : Discharge of node
    pvec : Pressure head of node (total head $-$ position head)
  .....(1 to NODT).....
*element value
element,xg,zg,vx,vz,vm,Ak0,Akr,matno
    element : Element number
    xg      : x-coordinate of COG of element
    zg      : z-coordinate of COG of element
    vx      : x-direction mean velocity in global coordinate system
    vz      : z-direction mean velocity in global coordinate system
    vm      : Combined mean velocity in global coordinate system
    Ak0     : Permeability coefficient (saturated)
    Akr     : Relative hydraulic conductivity function of element
    matno   : Material set number
  .....(1 to NELT).....
#,Summary
#,item,QTi,QTo
#,discharge,sum total of in-flow,sum total of out-flow
#,item,ne,vmax,xg,zg
#,vmmax,element number,max.velocity in all elements,x-coordinate of COG,,z-coordinate of COG
#,vimax,element number,max.velocity in in-flow elements,x-coordinate of COG,,z-coordinate of COG
#,vomax,element number,max.velocity in out-flow elements,x-coordinate of COG,,z-coordinate of COG
#,NELT=(Number of elements), NODT=(Number of nodes), nt=(nt), mm=(mm), ib=(ib)
#,nnn=(Number of iteration), icount=(Number of converged degrees of freedom)
#,Calculation time=(Calculation time)
#,Date_time=(date of execution)
    nt : Total degrees of freedom of FE equation
    mm : Dimension of reduced FE equation
    ib : band width of reduced FE equation

Pressure head analysis around the tunnel

• This section shows the example of analysis for tunnel under ground water.
• Domain for analysis is a square area with 200m-side.
• Tunnel is D-shape tunnel with 5m-width and 5m-height.
• Permeability of the ground rock is 1×10−7m/s and overburden of tunnel is 300m.
• Therefore, total head of 300m is acted at perimeters of analysis domain.

Cases for analysis are shown below.

Mesh

Case-0
Case-1
Case-2
Case-3
Case-4

Zoned embankment dam

In this section, examples of analysis for zoned embankment dam are shown. In addition, it is a wellknown fact for zoned embank dam that the free water surface serves as an action which has a permeating point on the core back and does not become such form shown in this analysis.

Mesh diagram
Pressure Head Distribution
Flow Velocity Vector
Pressure Head and Discharge

Homogeneous embankment dam

The seepage flow discharge by seepage flow FEM was compared with it by Casagrande's method.

Casagrande's method

Casagrande's model

First approximation

Second approximation

Since first approximation gives the result of safety side based on trial calculation, first approximation was used on the calculations in this document.

Result of seepage flow FEM

Conditions for analysis

Conditions for analysis and outline of results are shown below.

Sample of Model of Homogeneous embankment dam

Result of FEA

Result of comparison

Following matters are cleared from figures.

• The seepage flow discharge $Q$ by FEM is greater than it by Casagrande's method for each case.
• If the slope becomes steep, the difference of discharge between FEM and Casagrande's method becomes large.
• The distance $a$ from downstream edge of dam to seepage point by FEM is greater than it of Casagrande's method for each.
• If the slope becomes gentle, the difference of distance 'a' between FEM and Casagrande's method becomes large.

Discharge $Q$-Slope relationship on Homogeneous embankment dam	Distance $a$-Slope relationship on Homogeneous embankment dam

Programs and sample data

Fortran programs

Shell scripts

Input data sample for dam analysis

Input data sample for tunnel analysis

Expression of result of seepage flow analysis (f90_gmt_seep.f90)

Outline

This program can show the FEM mesh, pressure head contour map, total head contour map and flow velocity vector map by using GMT.

Some explanations are shown below.

• FEM mesh is drawn by specifying the coordinates of each element and repeating painting out of each element and drawing of each element boundary line.
• The contour map is created in the following procedures.
  • Creation of the mask file which defines the element existence domain by 'grdmask' command.
  • Creation of the grid file for contours by 'surface' command.
  • Creation of the grid file for plot by composition of the mask file and the grid file by 'grdmath' command.
  • Coloring by 'grdimage' command and contour drawing by 'grdcontour' command.
• The vector diagram is created by superimposing the arrow of the flow velocity vector on the colored map of the pressure head or total head generated by contour creation.
• In addition, since a figure will become ugly if the flow velocity vector is drawn on all elements, in this example, the flow velocity vector is drawn once per three elements.

Following files must be prepared before execution of the program.

Output file by FEM analysis is a output by the program 'f90_fem_seep' made by webmaster.

Description of Shell script

Shell script (a_draw_dam.txt) which carries out the process is shown below. Some files which are read from Fortran program are made in this batch file.

001   gfortran -o f90_gmt_seep f90_gmt_seep.f90
002
003   cp 0_cpt_dam.cpt 0_cpt.cpt
004
005   # define of mesh colour (1:MATEL)
006   cat << EOT > _col_mesh.txt
007   255 255 204
008     0 255 255
009   255   0 255
010     0 255 255
011   255 255 204
012   EOT
013   # integer::ind,ibc,iph,ivc,isc
014   # real(4)::baselength,grdspc
015   # real(4)::bav,bfv
016   # real(4)::conint,annint
017   # character::strx*100,strz*100
018   cat << EOT > _inp_base.txt
019   0 1 1 1 1
020   12.0 0.1
021   10.0 5.0
022    1.0 5.0
023   "Distance (m)"
024   "Altitude (m)"
025   EOT
026
027   ./f90_gmt_seep out_fem_zone3m.csv
028   mv -f _fig_gmt_mesh.eps fig_gmt_3mesh_dam.eps
029   mv -f _fig_gmt_conp.eps fig_gmt_3conp_dam.eps
030   mv -f _fig_gmt_vect.eps fig_gmt_3vect_dam.eps
031
032   ./f90_gmt_seep out_fem_zone4m.csv
033   mv -f _fig_gmt_mesh.eps fig_gmt_4mesh_dam.eps
034   mv -f _fig_gmt_conp.eps fig_gmt_4conp_dam.eps
035   mv -f _fig_gmt_vect.eps fig_gmt_4vect_dam.eps
036
037   ./f90_gmt_seep out_test0_3.csv
038   mv -f _fig_gmt_mesh.eps fig_gmt_mesh_earth.eps
039   mv -f _fig_gmt_conp.eps fig_gmt_conp_earth.eps
040   mv -f _fig_gmt_vect.eps fig_gmt_vect_earth.eps
041
042   rm _*

001

To compile the program 'f90_GMTseep.f90'

003

To copy '0_cpt_dam.cpt' which has been prepared beforehand to '0_cpt.cpt.' The contents of file '0_cpt_dam.cpt' is shown below.

-18 255 255 190  0 255 255 125
  0   0 255 255 18   0 125 255
B   255 255 190
F     0 125 255
N   255 255 255  N 255 255 255

006-012

To make the file '_col_mesh.txt.' Mesh colors are defined by this file, 3 figures means 'R G B' in this file. In this example, new file is made using 'cat' hear documents in this script, and since 5 materials are used (parameter MATEL=5), 5 kinds of color is defined.

018-025

To make the file '_inp_base.txt.' This file defines the parameter for GMT.

027

To execute the program 'f90_gmt_seep.' Functions of this program are to make GMT command files, to execute them and to make eps drawings.

Input file 'out_fem_zone4m.csv' which storages the result of FEM analysis is read from command line.

Output drawings are shown below.

028-030

Since output filename is fixed, it is necessary to change the output filename into another name by a 'cp' command etc after output.

042

To delete work files