3DEC Thermal and Dynamic Manual
3DEC Thermal and Dynamic Manual
©2016
Itasca Consulting Group, Inc. Phone: (1) 612-371-4711
Mill Place Fax: (1) 612-371-4717
111 Third Avenue South, Suite 450 Email: software@itascacg.com
Minneapolis, Minnesota 55401 USA Web: www.itascacg.com
First Edition May 1999
First Revision September 1999
Second Edition January 2003
First Revision August 2005
Third Edition December 2007
Fourth Edition (3DEC Version 5.0) July 2013
First Revision (3DEC Version 5.2) November 2016
Optional Features Contents - 1
TABLE OF CONTENTS
1 THERMAL OPTION
1.1 Numerical Thermal Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-2
1.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-2
1.1.2 Heat Conduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-2
1.1.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-2
1.1.2.2 Mathematical Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . 1-3
1.1.2.3 Numerical Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-6
1.1.2.4 Solving Thermal-Only and Thermal-Mechanical Problems . . . 1 - 13
1.1.3 Heat Convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 16
1.1.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 16
1.1.3.2 Mathematical Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 16
1.1.3.3 Numerical Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 17
1.1.3.4 Stability and Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 18
1.1.3.5 Solving Thermal Convection Problems . . . . . . . . . . . . . . . . . . . . . 1 - 19
1.1.4 Input Instructions for Thermal Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 21
1.1.4.1 3DEC Commands . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 21
1.1.4.2 FISH Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 26
1.1.5 Systems of Units for Thermal Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 27
1.1.6 Verification Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 29
1.1.6.1 Conduction in a Plane Sheet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 29
1.1.6.2 Heating of a Hollow Cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 39
1.1.6.3 Infinite Line Heat Source in an Infinite Medium . . . . . . . . . . . . . 1 - 48
1.1.6.4 One-Dimensional Thermal Transport by Conduction and Convec-
tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 61
1.1.7 Application Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 67
1.1.7.1 Thermal Convection in a Fracture at Depth . . . . . . . . . . . . . . . . . 1 - 67
1.2 Analytical Thermal Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 82
1.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 82
1.2.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 83
1.2.2.1 The Temperature Due to a Distribution of Heat Sources . . . . . . 1 - 83
1.2.2.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 84
1.2.2.3 Thermally Induced Stress Changes . . . . . . . . . . . . . . . . . . . . . . . . 1 - 85
1.2.3 Application in 3DEC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 86
1.2.4 Input Commands . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 86
1.2.4.1 Preparing to Solve a Thermal Problem . . . . . . . . . . . . . . . . . . . . . 1 - 86
1.2.4.2 Defining Heat Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 87
1.2.4.3 Initial and Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 89
1.2.4.4 Thermal Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 90
1.2.4.5 Calculating the Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 90
1.2.4.6 Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 90
1.2.5 Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 91
1.2.5.1 A Single Non-decaying Point Heat Source . . . . . . . . . . . . . . . . . . 1 - 91
1.2.5.2 Superposition of Several Non-decaying Sources . . . . . . . . . . . . . 1 - 93
1.2.5.3 Lines and Grids of Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 96
1.2.5.4 Decaying Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 101
1.2.5.5 Thermomechanical Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 103
1.2.6 An Example Problem – Waste Repository Drift Model . . . . . . . . . . . . . . 1 - 112
1.3 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 117
2 DYNAMIC ANALYSIS
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-1
2.2 Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-2
2.2.1 General Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-2
2.2.2 Rayleigh Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-2
2.2.3 Example of Different Damping Techniques . . . . . . . . . . . . . . . . . . . . . . . . 2-4
2.2.4 Guidelines for Selecting Rayleigh Damping Parameters for Dynamic Anal-
ysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-9
2.3 Natural Modes of Oscillation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 10
2.3.1 Example Problems (from Cundall et al. (1979), pp. 71-73) . . . . . . . . . . 2 - 11
2.4 Wave Transmission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 17
2.5 Partial Density Scaling for Dynamic Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 27
2.5.1 Example of Partial Density Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 27
2.6 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 32
2.6.1 Nonreflecting Boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 32
2.6.2 Free-Field Boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 34
2.6.2.1 Example Using Dynamic Free Field . . . . . . . . . . . . . . . . . . . . . . . 2 - 36
2.7 Application of Dynamic Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 42
2.7.1 Baseline Correction of Input Histories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 43
2.8 Calculation of Natural Frequencies and Modes of Vibration . . . . . . . . . . . . . . . . . 2 - 45
2.9 Verification Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 49
2.9.1 Slip on a Joint Induced by a Propagating Harmonic Shear Wave . . . . . . 2 - 49
2.9.1.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 49
2.9.1.2 Analytic Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 49
2.9.1.3 Numerical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 51
2.9.1.4 Material Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 52
2.9.1.5 Dynamic Loading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 53
2.9.1.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 53
TABLES
Table 1.1 Comparison of features available with numerical vs analytical formulations . . 1-1
Table 1.2 System of SI units for thermal problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 27
Table 1.3 System of Imperial units for thermal problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 27
Table 1.4 Property and parameter values for the example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 62
Table 1.5 Temperature distribution in the vicinity of several sources . . . . . . . . . . . . . . . . . . 1 - 95
Table 2.1 Moduli appropriate to various deformation modes . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 12
Table 2.2 Comparison of theoretical and calculated (3DEC) dynamic period T of oscilla-
tion for three modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 13
Table 2.3 Natural frequencies (Hz) of bending modes for the square pillar . . . . . . . . . . . . 2 - 47
Table 2.4 Medium properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 52
Table 2.5 Discontinuity properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 52
Table 2.6 Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 64
Table 2.7 Discontinuity properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 65
Table 3.1 Summary of structural element liner commands . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-8
Table 3.2 Systems of units – structural elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 - 13
Table 3.3 Nodal coordinates in master element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 - 22
Table 3.4 Coordinates of slave nodes in master element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 - 23
Table 3.5 Definition of the element faces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 - 23
FIGURES
Figure 1.1 3DEC grid for conduction in a plane sheet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 30
Figure 1.2 Comparison of temperatures for the explicit-solution algorithm (Tables 2, 4
and 6 = analytic solution; Tables 1, 3 and 5 = 3DEC solution) . . . . . . . . . . 1 - 34
Figure 1.3 Comparison of temperatures for the implicit-solution algorithm (Tables 2, 4
and 6 = analytic solution; Tables 1, 3 and 5 = 3DEC solution) . . . . . . . . . . 1 - 38
Figure 1.4 3DEC grid for heating of a hollow cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 41
Figure 1.5 Temperature distribution at steady state for heating of a hollow cylinder . . . 1 - 47
Figure 1.6 Radial stress distribution at steady state for heating of a hollow cylinder . . . 1 - 47
Figure 1.7 3DEC model for an infinite line heat source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 49
Figure 1.8 Temperature distribution at 1 year . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 51
Figure 1.9 Radial displacement distribution at 1 year . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 51
Figure 1.10 Radial and tangential stress distributions at 1 year . . . . . . . . . . . . . . . . . . . . . . . 1 - 52
Figure 1.11 Analytical and numerical temperature profiles, compared to conduction solu-
tion at steady state – case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 63
Figure 1.12 Analytical and numerical temperature profiles, compared to conduction solu-
tion at steady state – case 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 63
Figure 1.13 3DEC block geometry for thermal convection in a fracture . . . . . . . . . . . . . . . 1 - 68
Figure 1.14 3DEC block discretization for thermal convection in a fracture . . . . . . . . . . . . 1 - 68
Figure 1.15 Flow velocity in the fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 69
Figure 1.16 Fluid pressure distribution in the fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 70
Figure 1.17 Fluid temperature contours at steady state – case 1 . . . . . . . . . . . . . . . . . . . . . . . 1 - 71
Figure 1.18 History of fluid temperature at 5 monitoring points – case 1 . . . . . . . . . . . . . . . 1 - 72
Figure 1.19 History of rock temperature at 3 monitoring points – case 1 . . . . . . . . . . . . . . . 1 - 72
Figure 1.20 Fluid temperature contours at steady state – case 2 . . . . . . . . . . . . . . . . . . . . . . . 1 - 73
Figure 1.21 History of fluid temperature at 5 monitoring points – case 2 . . . . . . . . . . . . . . . 1 - 74
Figure 1.22 History of rock temperature at 3 monitoring points – case 2 . . . . . . . . . . . . . . . 1 - 74
Figure 1.23 Fluid temperature contours at steady state – case 3 . . . . . . . . . . . . . . . . . . . . . . . 1 - 75
Figure 1.24 History of fluid temperature at 5 monitoring points – case 3 . . . . . . . . . . . . . . . 1 - 76
Figure 1.25 History of rock temperature at 3 monitoring points – case 3 . . . . . . . . . . . . . . . 1 - 76
Figure 1.26 Adiabatic and isothermal boundaries and their 3DEC representations . . . . . . 1 - 85
Figure 1.27 A grid of heat sources, input with the command
GRID (10, 10, 10) (20, 20, 20) (20, 0, 10) &
TYP=1 STR=10 N12=5 N23=6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 89
Figure 1.28 Temperature distribution around a constant heat source . . . . . . . . . . . . . . . . . . . 1 - 93
Figure 1.29 Temperature distribution around a decaying heat source . . . . . . . . . . . . . . . . . . 1 - 102
Figure 1.30 Conceptual representation of 3DEC block for modeling effect of point heat
source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 105
Figure 1.31 3DEC model for instantaneous point source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 107
Figure 1.32 Radial stress (Szz ) on z-axis (x = y = 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 108
Figure 1.33 Tangential stresses Sxx and Syy along z-axis (x = y = 0) after 1 day . . . . . . . . 1 - 109
Figure 1.34 Tangential stresses Sxx and Syy along z-axis (x = y = 0) after 3 days . . . . . . . 1 - 110
Figure 1.35 z-displacement along z-axis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 111
Figure 1.36 Conceptual model of repository drift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 113
Figure 1.37 3DEC model of repository drift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 113
Figure 1.38 Model of drift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 114
Figure 1.39 Temperatures and displacements on drift center line . . . . . . . . . . . . . . . . . . . . . 1 - 114
Figure 2.1 Variation of normalized critical damping ratio with angular frequency . . . . . 2-4
Figure 2.2 Plot of vertical displacement versus time, for a single block contacting on a
rigid base with gravity suddenly applied (no damping) . . . . . . . . . . . . . . . . 2-7
Figure 2.3 Plot of vertical displacement versus time, for a single block contacting on a
rigid base with gravity suddenly applied (mass and stiffness damping) . . 2-7
Figure 2.4 Plot of vertical displacement versus time, for a single block contacting on a
rigid base with gravity suddenly applied (mass damping only) . . . . . . . . . 2-8
Figure 2.5 Plot of vertical displacement versus time, for a single block contacting on a
rigid base with gravity suddenly applied (stiffness damping only) . . . . . . 2-8
Figure 2.6 Plot of velocity spectrum versus frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-9
Figure 2.7 Comparison of fundamental wavelengths for bars with varying end conditions 2 - 11
Figure 2.8 Column of variably sized blocks subjected to triangular-shaped impulse load
at base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 18
Figure 2.9 Input wave (solid) at base and calculated wave (dashed) at top of column of
rigid block model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 19
Figure 2.10 Column of variably sized blocks subdivided into finite difference zones . . . . 2 - 19
Figure 2.11 Input wave (solid) at base and calculated wave (dashed) at top of column of
deformable block model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 20
Figure 2.12 Unfiltered velocity history . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 25
Figure 2.13 Unfiltered power spectral density plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 25
Figure 2.14 Filtered velocity history at 15 Hz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 26
Figure 2.15 Results of filtering at 15 Hz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 26
Figure 2.16 View of model with small, thin blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 28
Figure 2.17 Velocities at the bottom and top of the model, for analysis without any density
scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 29
Figure 2.18 Velocities at the bottom and top of the model, for analysis with partial density
scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 29
Figure 2.19 Model for seismic analysis of surface structures and free-field mesh . . . . . . . 2 - 35
Figure 2.20 Model of dam with free field blocks visible . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 37
Figure 2.21 x-velocity histories at model base and dam crest using viscous boundaries . 2 - 38
Figure 2.22 x-velocities of the base and dam crest using free field boundaries . . . . . . . . . . 2 - 39
Figure 2.23 The baseline correction process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 43
Figure 2.24 Square pillar model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 46
Figure 2.25 Transmission and reflection of incident harmonic wave at a discontinuity . . 2 - 49
Figure 2.26 Geometry for the problem of slip induced by harmonic shear . . . . . . . . . . . . . 2 - 51
EXAMPLES
Example 1.1 Conduction in a plane sheet – explicit solution . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 31
Example 1.2 Conduction in a plane sheet – implicit solution . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 35
Example 1.3 Heating of a hollow cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 41
Example 1.4 Infinite line heat source in an infinite medium . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 52
Example 1.5 Exponential integral function – “exp int.3dfis” . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 59
Example 1.6 FISH function – “fxface.3dfis” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 60
Example 1.7 One-dimensional thermal transport by conduction and convection . . . . . . . 1 - 64
Example 1.8 Thermal convection in a fracture at depth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 77
Example 1.9 Single non-decaying point heat source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 92
Example 1.10 Temperature distribution in the vicinity of several sources . . . . . . . . . . . . . . . 1 - 93
Example 1.11 Data file for point sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 96
Example 1.12 Data file for lines of sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 99
Example 1.13 Data file for grids of sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 100
Example 1.14 Temperatures due to a single decaying heat source . . . . . . . . . . . . . . . . . . . . . 1 - 101
Example 1.15 Thermomechanical verification of 3DEC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 105
Example 1.16 Waste repository drift model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - 115
Example 2.1 Data file for Rayleigh damping example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-5
Example 2.2 Data file for confined compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 13
Example 2.3 Data file for unconfined compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 14
Example 2.4 Data file for shear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 15
Example 2.5 Column of variably sized rigid blocks subjected to impulse load at base . . 2 - 20
Example 2.6 Column of variably sized deformable blocks subjected to impulse load at
base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 22
Example 2.7 No density scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 30
Example 2.8 Partial density scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 31
Example 2.9 Dynamic free field boundary example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 39
Example 2.10 Elastic vibration modes of a square pillar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 47
Example 2.11 MILR3D.3DDAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 56
Example 2.12 MILR3D.3DFIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 59
Example 2.13 DAY3D.3DDAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 67
Example 2.14 DAY3D.3DFIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 69
Example 2.15 VEL INP.3DFIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 72
Example 2.16 ANA SLP.3DFIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 - 74
Example 3.1 Structural liner in tunnel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-9
Example 3.2 Finite element model of a dam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 - 37
1 THERMAL OPTION
Introduction
The thermal option in 3DEC allows the simulation of transient heat conduction. There are two
separate formulations of the thermal logic. The first is a numerical formulation using the explicit
or implicit finite difference method. This method is more accurate for short times, and includes
thermal-mechanical fluid coupling. The second is an analytical formulation that uses superposition
of point heat sources in an infinite medium. This method is suitable for long thermal times, and is
very fast. The comparative features of the two methods are:
1.1.1 Introduction
The thermal option of 3DEC allows simulation of transient heat conduction in materials, and
development of thermally induced displacements and stresses. This option has the following specific
features:
1. Any of the mechanical models may be used with the thermal model.
2. Temperature, flux, convective and adiabatic boundary conditions may be prescribed.
3. Heat sources may be inserted into the material as either point sources or volume sources.
These sources may decay exponentially with time.
4. Both explicit- and implicit-solution algorithms are available.
5. The thermal option provides for one-way coupling to the mechanical stress calculations
through the thermal expansion coefficients.
6. Temperatures can be accessed via FISH.
7. Coupling between fluid flow in fractures and the thermal logic is available.
The thermal option also allows modeling of transient thermal convection in fractures filled with
a moving fluid. In this version of the code, fluid thermal expansion is not considered in the fluid-
thermo-mechanical coupling.
This section contains descriptions of the thermal formulation and the 3DEC input commands for
thermal analysis. Thermal conduction is considered in Section 1.1.2; thermal convection is reviewed
in Section 1.1.3. Recommendations for solving thermal problems are also provided, along with
several verification problems. Refer to these examples as a guide to creating 3DEC models for
thermal analysis and coupled thermal-stress or thermal-fluid flow analysis.
1.1.2.1 Introduction
The logic for heat conduction through a solid is available for deformable blocks. (It is not activated
for rigid blocks.) By default, convection is off and heat is conducted across the blocks using an
interpolation procedure. (The fractures are not seen by the flow of heat in this case.) Fluid-thermal
coupling can also be activated. (This topic is addressed in Section 1.1.3.)
The variables involved in heat conduction in 3DEC are temperature and the three components of
the heat flux. These variables are related through the energy-balance equation and transport laws
derived from Fourier’s law of heat conduction. Substitution of Fourier’s law into the energy-balance
equation yields the differential equation of heat conduction, which may be solved for particular ge-
ometries and properties, given specific boundary and initial conditions. Thermal volumetric-strains
are introduced into the incremental mechanical constitutive laws to account for thermomechanical
coupling.
Conventions and Definitions
As a notation convention, the symbol ai denotes component i of the vector [a] in a Cartesian system
of reference axes; Aij is component (i, j ) of tensor [A]. Also, f,i is used to represent the partial
derivative of f with respect to xi . (f can be a scalar variable or a vector component.)
The Einstein summation convention applies only to indices i, j and k, which take the values 1, 2,
3 for components that involve spatial dimensions. The indices may take any values when used in
matrix equations.
SI units are used to illustrate parameters and dimensions of variables. See Section 1.1.5 for con-
versions to other systems of units.
The following dimensionless numbers are useful in the characterization of transient heat conduction.
Characteristic length:
volume of solid
Lc = (1.1)
surf ace area exchanging heat
Thermal diffusivity:
k
κ= (1.2)
ρ Cv
Characteristic time:
L2c
tc = (1.3)
κ
Energy-Balance Equation
The differential expression of the energy balance has the form
∂ζ
−qi,i + qv = (1.4)
∂t
where qi is the heat-flux vector in [W/m2 ], qv is the volumetric heat-source intensity in [W/m3 ],
and ζ is the heat stored per unit volume in [J /m3 ]. In general, temperature changes may be caused
by changes in both energy storage and volumetric strain, . And the thermal constitutive law relating
those parameters may be expressed as
∂T ∂ζ ∂
= Mth − βth (1.5)
∂t ∂t ∂t
In 3DEC, we consider a particular case of this law for which βth = 0 and Mth = ρC1 v . ρ is the mass
density of the medium in [kg/m3 ], and Cv is the specific heat at constant volume in [J /kg ◦ C].
The hypothesis here is that strain changes play a negligible role in influencing the temperature – a
valid assumption for quasi-static mechanical problems involving solids and liquids. Accordingly,
we may write
∂ζ ∂T
= ρ Cv (1.6)
∂t ∂t
∂T
−qi,i + qv = ρCv (1.7)
∂t
Note that for nearly all solids and liquids, the specific heats at constant pressure and at constant
volume are essentially equal. Consequently, Cv and Cp can be used interchangeably.
Transport Law
The basic law that defines the relation between the heat-flux vector and the temperature gradient is
Fourier’s law. For a stationary, homogeneous, isotropic solid, this constitutive law is given in the
form
qi = −kT,i (1.8)
qn = h(T − Te ) (1.9)
where qn is the component of the flux normal to the boundary in the direction of the exterior normal,
h is the convective heat-transfer coefficient [W/m2 ◦ C], T is the temperature of the boundary surface
and Te is the temperature of the surrounding fluid [◦ C]. Note that in the numerical formulation used
in 3DEC, boundaries are adiabatic by default.
Mechanical Coupling – Thermal Strains
Solving thermal-stress problems requires reformulation of the incremental stress-strain relations,
which is accomplished by subtracting from the total strain increment that portion due to temperature
change. Because free thermal expansion results in no angular distortion in an isotropic material,
the shearing-strain increments are unaffected. The thermal-strain increments associated with the
free expansion corresponding to temperature increment T have the form
where αt [1/◦ C] is the coefficient of linear thermal expansion, and δij is the Kronecker delta.
The differential equation of motion and the rate of strain-velocity relations are based upon mechan-
ical considerations, and are unchanged for thermomechanics.
The energy-balance equation (Eq. (1.7)) and Fourier’s law (Eq. (1.8)) are solved in 3DEC using a
finite-difference approach based on a medium discretization into zones (tetrahedra). The numerical
scheme rests on a nodal formulation of the energy-balance equation. This formulation is equivalent
to the mechanical formulation (presented in Section 1 in Theory and Background) that leads to
the nodal form of Newton’s law. It is obtained by substituting the temperature, heat-flux vector and
temperature gradient for velocity vectors, stress tensors and strain-rate tensors, respectively. The
resulting system of ordinary differential equations is solved using two modes of discretization in
time, corresponding to explicit and implicit formulations. The principal results are summarized
below.
Finite-Difference Approximation to Space Derivatives
By convention, the nodes of a tetrahedron are referred to locally by a number from 1 to 4, face n is
opposite node n, and the superscript (f ) relates to the value of the associated variable on face f .
A linear temperature variation is assumed within a tetrahedron; the temperature gradient, expressed
in terms of nodal values of the temperature by application of the Gauss divergence theorem, may
be written as
1 l (l) (l)
4
T,j = − T nj S (1.11)
3V
l=1
where [n](l) is the exterior unit vector normal to face l, S is the face surface area and V is the
tetrahedron volume.
Nodal Formulation of the Energy-Balance Equation
The energy-balance Eq. (1.7) may be expressed as
qi,i + b∗ = 0 (1.12)
where
∂T
b∗ = ρCv − qv (1.13)
∂t
is the equivalent of the instantaneous “body forces” used in the mechanical node formulation. First
consider a single tetrahedron. Using this analogy, the nodal heat Qne [W ], n = 1,4, in equilibrium
with the tetrahedron heat flux and body forces, may be expressed as
qv V dT n
Qne = Qnt − + mn Cvn (1.14)
4 dt
where
(n)
qi ni S (n)
Qnt = (1.15)
3
and
ρV
mn = (1.16)
4
In principle, the nodal form of the energy-balance equation is established by requiring that, at each
global node, the sum of equivalent nodal heat (−Qne ) from all tetrahedra meeting at the node and
nodal contribution (Qnw ) of applied boundary fluxes and sources be zero.
The components of the tetrahedron heat-flux vector in Eq. (1.15) are related to the temperature
gradient by means of the transport law (see Eq. (1.8)). In turn, the components of the temperature
gradient can be expressed in terms of the tetrahedron’s nodal temperatures, using Eq. (1.11).
In order to save computer time, a zone formulation is adopted in 3DEC. At each node, n, of a
particular zone, the sum, Qnz , of contributions from all tetrahedra belonging to the zone and meeting
at the node is formed and divided by two for two overlays. Local zone matrices, [M], that relate
nodal zone values, Qnz , to nodal zone temperatures, T n (n = 1,4), are assembled. These matrices
being symmetrical, 10 components are computed and saved at the beginning of the computation,
and updated every 10 steps in large-strain mode. By definition of zone matrices, we have
where [C] is the global matrix, and [T ] is, in this context, the global vector of nodal temperatures.
− Qne + Qnw = 0 (1.19)
where, for simplicity of notation, the sign is used to represent summation of the contributions
at global node n of all zones meeting at that node. Using Eq. (1.14) in Eq. (1.19), we obtain, after
some manipulations,
dT n 1
= − Qn
T + Qn
app (1.20)
dt [mCp ]n
where QnT is a function of the nodal temperatures defined in Eq. (1.18), and Qnapp is the known
contribution of applied volume sources, boundary fluxes and point sources, defined as
qv V n
Qnapp =− + Qw (1.21)
4
Eq. (1.20)
is the nodal form of the energy-balance equation at node n; the right-hand side term,
QnT + Qnapp , is referred to as the “out-of-balance heat.” One such expression holds at each global
node involved in the discretization. Together they form a system of ordinary differential equations
that is solved in 3DEC using both explicit and implicit finite-difference schemes. The domain of
application of each scheme is discussed below.
Explicit Finite-Difference Formulation
In the explicit formulation, the temperature at a node is assumed to vary linearly over a time interval,
t: the derivative in the left-hand side of Eq. (1.20) is expressed using forward finite differences;
and the out-of-balance heat is evaluated at time t. Starting with an initial temperature field, nodal
temperatures at incremental time values are updated, provided the temperature is not fixed, using
the expression
n
T<t+t> = T<t>
n
+ T<t>
n
(1.22)
where
n
T<t> = χ n QnT <t> + Qnapp<t> (1.23)
and
t
χn = − (1.24)
[mCp ]n
Numerical stability of the explicit scheme can only be achieved if the timestep remains below a
limiting value.
Stability Criterion
To derive the stability criterion, we consider the situation in which a node, n, in an assembly of
zones is given a temperature perturbation, T0 , from a zero initial state. Using Eq. (1.18), we obtain
Qnapp = Dnn To (1.26)
where Dnn is used to represent the temperature coefficient in the global convective term at node n.
After one thermal timestep, the new temperature at node n is (see Eqs. (1.22) through (1.24))
t
n
T<t> = T0 1 − (Cnn + Dnn ) (1.27)
[mCp ]n
To prevent alternating signs of temperatures as the computation is repeated for successive t, the
coefficient of T0 in the preceding relation must be positive. Such a requirement implies that
[mCp ]n
t < (1.28)
Cnn + Dnn
The value of the timestep used in 3DEC is the minimum nodal value of the right-hand side in
Eq. (1.28), multiplied by the safety factor of 0.8.
To assess the influence of the parameters involved, it is useful to keep in mind the following
representation of the critical timestep. If Lc is the smallest tetrahedron characteristic length, we
may write an expression of the form
−1
1 κ h
tcr = + (1.29)
m L2c ρ Cv Lc
where tcr is the critical timestep, κ, h, ρ and Cv are defined previously in Section 1.1.2.2, and m
is a constant, larger than unity, which depends on the medium geometrical discretization (e.g., m =
6 for a regular discretization in cubes – see Karlekar and Desmond 1982).
The critical timestep in Eq. (1.29) corresponds to a measure of the characteristic time needed for
the diffusion “front” to propagate throughout the tetrahedron. To estimate the time needed for the
front to propagate throughout a particular domain, a similar expression can be used, provided Lc is
interpreted as the characteristic length of the domain under consideration. Consider the case where
no convection occurs. By taking the ratio of characteristic times for the domain and the tetrahedron,
it may be seen that the number of steps needed to model the propagation of the diffusion process
throughout that domain is proportional to the ratio of square characteristic lengths for the domain
and the tetrahedron. That number may be so large that the use of the explicit method alone could
become prohibitive.
On the other hand, the advantage of this first-order method is that the calculated timestep is usually
small enough to follow nodal temperature fluctuations accurately.
Implicit Finite-Difference Formulation
The requirement that t should be restricted in size to ensure stability sometimes results in an ex-
tremely small timestep – of the order of a fraction of a second, especially when transient conduction
in multiple layers is involved. The implicit formulation eliminates this restriction, but it involves
solving simultaneous equations at each timestep.
The implicit formulation in 3DEC uses the Crank-Nicolson method, in which the temperature at
a node is assumed to vary quadratically over the time interval t. In this method, the derivative
dT /dt in Eq. (1.20) is expressed using a central-difference formulation corresponding to a half-
timestep while the out-of-balance heat is evaluated by taking average values at t and t + t. In this
formulation, we have
n
T<t+t> = T<t>
n
+ T<t>
n
(1.30)
where
1
n
T<t> =χ n
QnT <t+t> + QnT <t> + Qnapp (1.31)
2
and
1 n
Qnapp = Qapp<t+t> + Qnapp<t> (1.32)
2
j
QnT <t> = Cnj T<t> (1.33)
and
j
QnT <t+t> = Cnj T<t+t> (1.34)
j
QnT <t+t> = QnT <t> + Cnj T<t> (1.35)
After substitution of Eq. (1.35) into Eq. (1.31), we obtain, using Eq. (1.33),
j 1 j
n
T<t> =χ n
Cnj T<t> + Cnj T<t> + Qnapp (1.36)
2
χn j j
δnj − Cnj T<t> = χ n Cnj T<t> + Qnapp (1.37)
2
For simplicity of notation, we define the known matrix [A] and vector [b<t> ] as
χn
Anj = δnj − Cnj (1.38)
2
and
j
bn<t> = χ n Cnj T<t> + Qnapp (1.39)
j
Anj T<t> = bn<t> (1.40)
Such an equation may be written for each of the nodes involved in the grid. The resulting system
of equations is solved in 3DEC using Jacobi’s iterative method. In this approach, temperature
increments at iteration r + 1 and node n are calculated using the recurrence relation,
1 j <r>
n<r+1>
T<t> = T<t>
n<r>
+ −Anj T<t> + bn<t> (1.41)
Ann
where Einstein’s notation convention applies to j indices only. Using the definition Eq. (1.38) of
[A], Eq. (1.41) takes the form
1 χn j <r>
n<r+1>
T<t> = n<r>
T<t> + χn
Cnj T<t> − T<t>
n<r>
+ bn<t> (1.42)
1− 2 Cnn
2
n<0>
T<t> =0 (1.43)
max T<t>
n<r+1> n<r>
− T<t> < 10−2 max T<t>
n<r>
(1.44)
n n
The magnitude of the timestep must be selected in relation to both convergence and accuracy of the
implicit scheme. Although the Crank-Nicolson method is stable for all positive values of t (for
no convection), the convergence of Jacobi’s method is not unconditionally guaranteed unless the
matrix [A] is strictly diagonally dominant – i.e., when
n
Akj < Akk (1.45)
j =1
j =k
for 1 ≤ i ≤ n (see Dahlquist and Bjorck 1974). According to the definition Eq. (1.38) of Anj and
Eq. (1.24) of χ n , this sufficient condition is always fulfilled for sufficiently small values of t. If
convergence of Jacobi’s method is not achieved, an error message is issued. It is then necessary
to either reduce the magnitude of the implicit timestep, or use the explicit method. (The explicit
timestep may be used as a lower bound.)
Also, although the implicit method is second-order accurate, some insight may be needed to select
the appropriate timestep. Indeed, its value must remain small compared to the wavelength of
any nodal temperature fluctuation. Typically, the explicit method is used earlier in the run or
in its perturbed stages, while the implicit method is preferred for the remainder of the simulation.
(Alternatively, the implicit method could be used with the explicit timestep value for extra accuracy.)
Computation time and computer memory are two factors that must be taken into consideration when
selecting the implicit approach in 3DEC. In the implicit method, a set of equations must be solved
at each timestep requiring a minimum of three iterations. The amount of calculation required for
one iterative step is approximately equal to that needed for one timestep in the explicit scheme.
Also, intermediate values must be stored in the iterative procedure, requiring more memory to be
allocated than in the explicit scheme. Those inconveniences, however, can be more than offset by
the much larger timestep generally permitted by the implicit method, or by the gain in accuracy
allowed.
Thermal-Stress Coupling
The heat transfer may be coupled to thermal-stress calculations at any time during a transient
simulation. The coupling occurs in one direction only (i.e., the temperature may result in stress
changes, but mechanical changes in the body resulting from force application do not result in
temperature change). This restriction is not believed to be of great significance here since the
energy changes for quasi-static mechanical problems are usually negligible. The stress change in
a triangular zone is given by (from Eq. (1.10))
The preceding assumes a constant temperature in each triangular zone, which is interpolated from
the surrounding gridpoints. This stress is added to the zone stress state prior to application of the
constitutive law.
Mechanical properties can also be made a function of temperature change by accessing temperature
and property values via FISH.
1.1.2.4 Solving Thermal-Only and Thermal-Mechanical Problems
3DEC has the ability to perform thermal-only and coupled thermal-mechanical and thermal-fluid
flow analyses. The forms of the coupled interactions are described in Section 1.1.2.4. In all cases,
the first step in the command procedure is to issue a CONFIG thermal command so that extra memory
can be assigned for the thermal calculation. Then, a choice must be made between explicit and
implicit thermal algorithms. By default, the explicit algorithm will be selected, but the implicit
mode of calculation may be activated or deactivated at any stage of the thermal-only or coupled
calculation using the command SET th implicit on or SET th implicit off. In the explicit mode, the
thermal timestep will be calculated automatically, but a smaller timestep can be selected using the
SET thdt command. The magnitude of the timestep must be specified by the user in the implicit
mode. This is done by issuing a SET thdt command. The thermal model and properties must be
specified for any zones that conduct heat. The density must also be specified. Initial and boundary
conditions are assigned to complete the thermal problem setup. Flux boundary conditions, for
instance, are assigned through a range to the boundaries of the thermal domain. The thermal
ρ
tmech = Lc (1.47)
K + (4/3)G
Hence, the ratio of thermal-to-mechanical timestep may be expressed as (see Eq. (1.29), assuming
no convection and m = 1)
tther K + 4/3G Lc
= (1.48)
tmech ρ κ
where κ is the thermal diffusivity. For nonmetals, this property is of the order of 10−6 m2 /s, at
most. For rock and soil, ρ is of the order of 103 kg/m3 , while K + (4/3)G is approximately
1010 N/m2 . Using those orders of magnitude in Eq. (1.48), it may be observed that the ratio of
thermal-to-mechanical time scales is of the order 105 Lc . This ratio remains very large, even if Lc
is of the order of 1 mm. In practice, mechanical effects can be assumed to occur instantaneously
when compared to diffusion effects. This is also the approach adopted in 3DEC, where no time is
associated with any of the mechanical steps taken in association with the thermal steps.
The second issue concerns the thermal-mechanical coupling in 3DEC. This coupling occurs only
in one direction: temperature changes cause thermal strains to occur, which influence the stresses;
while the thermal calculation is unaffected by the mechanical changes taking place.
In most modeling situations, the initial mechanical conditions correspond to a state of equilibrium
that must first be achieved before the coupled analysis is started. If the medium is elastic and the
thermal-mechanical response must be investigated, for instance, at a certain thermal time, a thermal-
only calculation may be performed until the desired time (SET thermal on mech off). The thermal
calculation may then be turned off, and the mechanical calculation on, using the SET thermal off
and SET mech on commands. Following this, the system may be cycled to mechanical equilibrium
before the response is analyzed.
For nonlinear mechanical models (i.e., when plasticity is involved), the thermal changes must be
communicated to the mechanical module at closer time intervals, to respect the path-dependency of
the system. Typically, at small dimensionless thermal time, a certain number of mechanical steps
must be taken for each thermal step to allow the system to adjust according to the different time
scales involved. At large dimensionless thermal time, if the system approaches thermal equilibrium,
several thermal timesteps may be taken without significantly disturbing the mechanical state of the
medium. A corresponding numerical simulation may be controlled manually by alternating between
thermal-only and mechanical-only modes.
In a third approach, the STEP command may be used while both mechanical and thermal modules
are on (SET thermal on mech on). In this option, one mechanical step will be taken for each
thermal step. Here, thermal steps are assumed to be so small that one mechanical step is enough
to re-equilibrate the system mechanically after each thermal step is taken. Section 1.1.4.1 may be
consulted for a complete list of available command options.
1.1.3.1 Introduction
In the convection logic, it is assumed that fluid flow occurs within the fractures, which are satu-
rated with fluid, and that the rock matrix is impermeable to fluid flow (3DEC approach). Heat is
transported by convection by the fluid, and by conduction within it. Heat can also be transported
by conduction in the rock, as described in the previous section. The fluid has its own temperature,
which will generally be different than the rock. Heat transfer between the fluid in the fractures and
the contacting rock (fluid-thermal coupling) may take place according to Newton’s law of cooling.
The logic for heat transfer within the fluid, and coupling to heat transfer within the rock, is presented
in the following section.
1.1.3.2 Mathematical Model Description
Heat convection in the flow planes is described by the following equations. Heat is transported by
conduction in the fracture fluid, according to Fourier’s law:
where qfT is specific heat flux in the fluid, and kfT is fluid thermal conductivity. The energy-balance
equation for the fluid obeys the equation
∂Tf
ρf Cf + ∇ · qfT + ρf Cf q f · ∇Tf + Af h(Tf − Ts ) = 0 (1.50)
∂t
where ρf Cf is fluid density times specific heat; q f is specific fluid discharge; Af is contact area
per unit fluid volume; h is fluid/rock heat transfer coefficient; and Tf , Ts are the temperature of
fluid and solid block, respectively.
For the solid blocks, flow of fluid is neglected; transport of heat obeys Fourier’s law:
q T = −k T ∇T (1.51)
where q T is specific heat flux, and k T is rock thermal conductivity. The energy balance is
∂Ts
ρs Cs + ∇ · qsT − As h(Tf − Ts ) = 0 (1.52)
∂t
where ρs Cs is solid density times specific heat; and As is contact area per unit volume of solid (seen
from the fluid, there is contact on 2 sides: A+ − +
s , As , and Af = As + As ).
−
The coupling between thermal and fluid flow is assumed to be weak for this development (i.e.,
fluid pressure generation by fluid thermal expansion is neglected compared to the other effects).
Also, the dependence of fluid properties (viscosity, density) on temperature is not directly taken
into account in the current formulation. Finally, fluid thermal expansion is not considered in the
thermal-mechanical coupling for the development.
1.1.3.3 Numerical Formulation
The logic for thermal conduction in the solid blocks is reported in Section 1.1.2. We focus here on
the implementation of heat convection in the fractures.
The numerical implementation of the heat equations for the fracture fluid is done using an explicit,
finite volume scheme. The implementation of heat conduction in the fluid is done using a scheme
similar to the one used to solve the flow transport and fluid mass conservation equations in the
fractures (the equations have a similar form). The numerical technique relies on the discretization
of the flow planes into triangular zones, in which temperature is assumed to vary linearly. The
numerical formulation in a flow plane follows along the lines of the FLAC thermal scheme (for
one overlay), and the reader is referred to the Thermal Option section in the Optional Features
volume of the FLAC Manual for a description of the technique. The difference here is that two
additional contributions are added to the out-of-balance flux at the nodes, before new temperatures
are evaluated for the thermal step.
The first contribution accounts for heat convection. It is calculated on a flow-zone basis. A heat
source term of the form ρf Cf Qf · ∇Tf A is computed at each thermal timestep, for a triangular
zone of area, A, using the information of fluid specific discharge, Qf , from the flow calculation,
and temperature gradient for the zone. (The symbol Qf stands for heat flux integrated over fracture
thickness. The temperature gradient is evaluated from nodal temperatures using the Gauss diver-
gence theorem.) The heat source term is divided into three equal parts, and one is added to the
out-of-balance flux at each node of the zone.
The second contribution accounts for heat exchange between the fracture fluid and the rock, accord-
ing to Newton’s law of cooling. It is calculated on a flow-node basis. The principle is as follows.
Let S be the nodal area allocated to the node in the flow plane. The fluid temperature at the flow
node is Tf . The temperature of the rock is Ts+ on one side, and Ts− on the other side. For the fluid,
the surface interaction term h(Tf − Ts+ )S + h(Tf − Ts− )S is added to the out-of-balance flux at the
flow node, at each thermal timestep. The contribution h(Tf − Ts+ )S is distributed, with reversed
sign to the nodes of the contacting block face on one side, and the same is done for h(Tf − Ts+ )S
on the other side. New nodal temperatures are calculated for the step using the standard approach.
The stability of the numerical scheme for solving the energy balance equation Eq. (1.50) is governed
by two numbers: the Peclet number for the grid, Peg , and the grid Courant number, Crg . For uniform
specific heat, Cf , and, in the absence of fluid/rock interaction, the energy balance equation can be
written in the form
∂Tf
= ∇ · (D∇T ) − q f · ∇T (1.53)
∂t
where, by definition,
kfT
D= (1.54)
ρf Cf
The grid Peclet number gives a measure of the relative effects of convection and conduction, and
is defined as
qf
Peg = (1.55)
D/l
q f t
Crg = (1.56)
l
In the one-dimensional case, for instance, the classical constraints are Peg ≤ 2 and Crg ≤ 1 (see, for
example, Perrochet and Berod 1993). In the current implementation, the explicit fluid and thermal
timesteps are identical to those adopted in the conduction logic. The code thus sets no intrinsic
limit on the grid Peclet number. However, the accuracy is likely to be different for different number
values, and to be influenced by grid size and timestep.
Thus, the explicit timestep adopted by the code is not unconditionally stable; more research needs
to be done, and tests conducted, to resolve this issue. One recommendation for the user of this
development is, before embarking on a large run, to exercise and check the stability on the sim-
ple one-dimensional model presented in Section 1.1.6.4, using representative properties, specific
discharge and grid sizes planned for the problem. Although for most practical applications, the
heat conducted in the fluid will be small compared to the amount being advected, it is often ad-
vantageous, from a numerical point of view, to take conduction into account in the simulation. If
instability occurs in the small model (see preceding stability criterion for grid Peclet and Courant
numbers), the grid size should be reduced until stability is achieved. The large model discretization
will then need to be reviewed in light of the stable grid size obtained for the small model.
1.1.3.5 Solving Thermal Convection Problems
To perform a coupled fluid-thermal calculation with 3DEC, fluid and thermal configurations must
be selected, using the command CONFIG thermal gw. The convection mode is selected (using SET
convection on) to allow heat to propagate in the fractures by conduction in the fluid, and by convec-
tion from the flow of fluid. The fluid thermal properties of specific heat and thermal conductivities
are assigned with the commands FLUID fspec heat and FLUID fth cond, respectively. The initial
fluid temperature and thermal boundary conditions for the fluid must be assigned independently
from the initial rock temperature and thermal boundary conditions for the rock. Fluid temperature
is initialized at knots of flow planes using the INSITU fluidtemp command. Three types of fluid
thermal boundary conditions can be applied in the fractures corresponding to fixed temperature
(BOUND fluidtemp), surface heat flux (BOUND flux) and joint heat source (BOUND knflux).
To model heat transfer between fluid and rock, the problem of heat conduction in the rock must be
set up using the approach described in Section 1.1.2. Also, the property of heat transfer coefficient
can be assigned using the FLUID fht coe command (see Section 1.1.7 for an example application).
Two commands (LIST tknot and LIST fluid) are available to print the temperature at knots in the flow
plane and fluid properties, respectively. Finally, a contour plot of fluid temperature can be produced
by issuing the command PLOT fluidtemp.
A list of commands specific to heat convection simulations is included below. (SI units are indicated
in parentheses.)
BOUND fluidtemp <val> <range>
fixes fluid temperature at a given value for knots in range. [C]
BOUND flux <val> <range>
fixes heat flux at a given value for knots in range. Flux is per
unit length of joint trace on boundary surface. [W/m2 ]
BOUND knflux <val> <range>
fixes heat point source at a given value for knots in range.
[W/m3 ]
FLUID fht coe <value> [W/(m2 C)]
fluid heat transfer coefficient (global value)
FLUID fspec heat <value>
fluid specific heat [J/(kg C)]
The following commands are provided to run thermal problems. Note that the thermal options
are invoked by new keywords used with existing commands in the standard mechanical code.
The command CONFIG thermal must be the first thermal command given before any other thermal
commands are invoked.
Zone-Type Keyword
vsource v
A heat-generating source, v, is applied as a volume source of the
specified strength (e.g., in W/m3 ) in each zone in the specified
range. When a new source is applied to a zone with the exist-
ing source, the new source strength replaces the existing source
strength.
Decay of the heat source can be represented by a FISH history
using the history keyword. See the flux keyword for an example.
Face-Type Keywords
convection v1 v2
v1 is the temperature, Te , of the medium to which convection
occurs.
v2 is the convective heat-transfer coefficient, h (e.g., in W/m2 ◦ C).
A convective boundary condition is applied over the range of faces
specified. The history keyword is not active for convection.
flux v
v is the initial flux (e.g., in W/m2 ).
A flux is applied over the range of faces specified. This command is
used to specify a constant flux into (v > 0) or out of (v < 0) a thermal
boundary of the grid. Decay of the flux can be represented by a
FISH history using the optional keyword history. For example,
the following FISH function performs an exponential decay of the
applied flux:
def decay
decay=exp(deconst*(thtime-thini))
end
set thini=0.0 deconst=-1.0
apply flux=10 hist=decay
knflux fixes heat point source at a given value for knots in range.
CONFIG thermal
This command specifies extra memory to be assigned to each zone or gridpoint for
a thermal analysis.
property thermal
thermal properties
tknot prints temperature information at flow knots.
zone thermal
zone temperatures and heat flux
PLOT flowplanecontour
temperature
clock t
t is the computer runtime limit, in minutes (default clock time is
unlimited time)
thtime t
t is the thermal “heating-time” limit for the thermal calculation.
1.1.4.2 FISH Variables
The following scalar variables are available in a FISH function to assist with thermal analysis:
thdt thermal timestep
thtime thermal time
The following 3DEC grid variables can be accessed and modified by a FISH function:
gp temp gridpoint temperature
fk temp flow knot temperature
All thermal quantities must be given in an equivalent set of units; no conversions are performed by
the program. Tables 1.2 and 1.3 present examples of consistent sets of units for thermal parameters.
Decay Constant s
−1 s
−1 s
−1 s
−1
where:1K = 1.8 R;
1J = 0.239 cal = 9.48 × 10−4 Btu;
1J/kg K = 2.39 × 10−4 btu/1b R;
1W = 1 J/s = 0.239 cal/s = 3.412 Btu/hr;
1W/m K = 0.578 Btu/(ft/hr R); and
1W/m2 K = 0.176 Btu/ft2 hr R.
Note that temperatures may be quoted in the more common units of ◦ C (instead of K) or ◦ F (instead
of R), where:
Several verification problems are presented to demonstrate the thermal model in 3DEC. The data
files for these examples are contained in the “\Options\Thermal” directory.
1.1.6.1 Conduction in a Plane Sheet
∞
z 2 −κn2 π 2 t/L2 T2 cos(nπ) − T1 nπz
T (z, t) = T1 + (T2 − T1 ) + e sin (1.57)
L π n L
n=1
∞
T (z, t) z 2 −κn2 π 2 t/L2 sin nπL z
= 1− − e (1.58)
T1 L π n
n=1
The thermal conductivity for this example is 1.6 W/m◦ C, the specific heat is 0.2 J/kg◦ C, the mass
density of the material is 1000 kg/m3 , and the temperature, T , is 100◦ C.
The analytical solution is programmed as a FISH function for direct comparison to the numerical
results at selected thermal times. The analytical and numerical temperature results for these times
are stored in tables.
In the 3DEC model, the sheet is defined as a column of 25 zones. A constant temperature boundary
of 100◦ C is applied at the face located at z = 0, and a constant temperature boundary of 0◦ C is
applied at the face located at z = 1. The model grid is shown in Figure 1.1:
E
pGre00 0 u0 i
0rGGG
pGrelG,lrG03:,G:S,0
:0
r
Example 1.1 contains the 3DEC data file for this problem using the explicit formulation to obtain
the solution. Example 1.2 contains the data file using the implicit formulation. The comparison of
analytical and numerical temperatures at three thermal times for the explicit solution is shown in
Figure 1.2; the comparison for the implicit solution is shown in Figure 1.3. Normalized temperature
(T /T1 ) is plotted versus normalized distance (z/L) in the two figures, where Tables 2, 4 and 6 contain
the analytical solution for temperatures, and Tables 1, 3 and 5 contain the 3DEC solutions. The
three thermal times are 1.455, 7.273 and 72.73 seconds for the explicit solution, and 1.455, 11.45
and 71.45 seconds for the implicit solution. The solution has reached the equilibrium thermal state
by the last time in each case. For both solution formulations, the difference between analytical and
numerical temperatures at steady state is less than 0.1%.
E
pGre00 0 u0 i
0rGGG
pGrelG,lrG03:,G:S,04
r0,0Gi,0
p04 0Gi,0
,0,0riS0
.04 0riS0
S0,0rS00
e04 0rS00
E
pGre00 0 u0 i
0r-e
pGrelG,lrG0-:,G:S.0y
r0,0Gi,0
p0y 0Gi,0
,0,0riS0
.0y 0riS0
S0,0rS00
e0y 0rS00
A hollow cylinder of infinite length is initially at a constant temperature of 0◦ C. The inner radius
of the cylinder is exposed to a constant temperature of 100◦ C, and the outer radius is kept at 0◦ C.
The problem is to determine the temperatures and thermally induced stresses in the cylinder when
the equilibrium thermal state is reached.
Nowacki (1962) provides the solution to this problem in terms of the temperatures and radial,
tangential and axial stresses at the steady-state thermal state:
T (r) ln(b/r)
= (1.59)
Ta ln(b/a)
σr (r) ln(b/r) (b/r)2 − 1
= − − (1.60)
mGTa ln(b/a) (b/a)2 − 1
σt (r) ln(b/r) − 1 (b/r)2 + 1
= − + (1.61)
mGTa ln(b/a) (b/a)2 − 1
σa (r) 2 ln(b/r) − 2(λ+G)
λ
λ 2
= − + (1.62)
mGTa ln(b/a) 2λ + G (b/a)2 − 1
where: T = temperature;
r = radial distance from the cylinder center;
a = inner radius of the cylinder;
b = outer radius of the cylinder;
Ta = temperature at the inner radius;
σr = radial stress;
σt = tangential stress;
σa = axial stress;
3Kαt
m = λ+2G ;
λ = K − 23 G;
K is the bulk modulus;
G is the shear modulus; and
αt is the linear thermal expansion coefficient.
The analytical solutions for temperature and stresses are programmed as FISH functions in the
3DEC data file. The analytical and numerical results can then be compared directly in tables.
The following properties are prescribed for this example:
Geometry
inner radius of cylinder (a) 1.0 m
outer radius of cylinder (b) 2.0 m
Material Properties
density (ρ) 2000 kg/m3
specific heat (Cp ) 880.0 J/kg ◦ C
thermal conductivity (k) 4.2 W/m ◦ C
linear thermal expansion coefficient (αt ) 5.4 × 10−6 / ◦ C
shear modulus (G) 28.0 GPa
bulk modulus (K) 48.0 GPa
A thin quarter-section of the cylinder is modeled with 3DEC. Figure 1.4 shows the 3DEC grid.
A constant-temperature boundary of 100◦ C is specified for the inner radius of the model; the
temperature at the outer radius is specified to be 0◦ C.
The 3DEC model can be run as either an uncoupled or a coupled thermal-mechanical analysis. For a
thermal-elastic analysis, it is more efficient to perform an uncoupled analysis. This is demonstrated
by running both uncoupled and coupled models for this problem. We first run the model in an
uncoupled mode: the thermal calculation is performed first to reach the equilibrium heat-flux state;
then the thermally induced mechanical stresses are calculated. The 3DEC data file is listed in
Example 1.3.
E
pGre00 0 u0 i
07G7r
pGrelG,lrG0-8,G8S70
80
r
rrn = re
x1 = xc + rr * cos(ang)
z1 = zc + rr * sin(ang)
x2 = xc + rrn * cos(ang)
z2 = zc + rrn * sin(ang)
x3 = xc + rrn * cos(angn)
z3 = zc + rrn * sin(angn)
x4 = xc + rr * cos(angn)
z4 = zc + rr * sin(angn)
command
poly prism &
a @x1 @ya @z1 @x2 @ya @z2 @x3 @ya @z3 @x4 @ya @z4 &
b @x1 @yb @z1 @x2 @yb @z2 @x3 @yb @z3 @x4 @yb @z4
endcommand
endloop
end
;
@genb
ns = 1
loop while ns < nn
x = (xtable(tab1,ns) + xtable(tab1,ns+1)) * 0.5
; p_z = z_near(x,0.,0.)
p_z = z_near(x,0.05,0.)
xc = z_x(p_z)
zc = z_z(p_z)
ra2 = xc*xc + zc*zc
ra = sqrt(ra2)
xtable(tab3,ns) = ra
val=(z_sxx(p_z)*xc*xc + z_szz(p_z)*zc*zc + 2.*z_sxz(p_z)*xc*zc)/ra2
ytable(tab3,ns) = val / (c_mmu * t1)
xtable(tab5,ns) = ra
val=(z_sxx(p_z)*zc*zc + z_szz(p_z)*xc*xc - 2.*z_sxz(p_z)*xc*zc)/ra2
ytable(tab5,ns) = val / (c_mmu * t1)
xtable(tab7,ns) = ra
val=z_syy(p_z)
ytable(tab7,ns) = val / (c_mmu * t1)
xtable(9,ns) = ra
val=((z_szz(p_z) - z_sxx(p_z))*xc*zc + z_sxz(p_z)*(xc*xc-zc*zc))/ra2
ytable(9,ns) = val / (c_mmu * t1)
ns = ns + 1
end_loop
end
def ana_solst ; table tab1 must be available
ns = 1
loop while ns < nn
x = (xtable(tab1,ns) + xtable(tab1,ns+1)) * 0.5
p_z = z_near(x,0.,0.)
xc = z_x(p_z)
zc = z_z(p_z)
ra = sqrt(xc*xc + zc*zc)
xtable(tab4,ns) = ra
val = c_b / ra
ytable(tab4,ns) = - (ln(val)*oc1 - (val * val - 1.)*oc2)
xtable(tab6,ns) = ra
ytable(tab6,ns) = - ((ln(val)-1.)*oc1 + (val * val + 1.)*oc2)
xtable(tab8,ns) = ra
ytable(tab8,ns) = - ((2.*ln(val)-oc3)*oc1 + 2.*oc3*oc2)
ns = ns + 1
end_loop
end
;
@num_solt
@ana_solt
table 1 name ’3DEC’
table 2 name ’Analytical’
plot create plot ’Table
plot table 1 style mark 2 xaxis label ’Radius’ &
Yaxis Label ’Temperature (T/Ta)’
@num_solst
@ana_solst
;
table 3 name ’3DEC Radial Stress’
table 4 name ’Analytical Radial Stress’
table 5 name ’3DEC Tangential Stress’
table 6 name ’Analytical Tangential Stress’
table 7 name ’3DEC Axial Stress’
table 8 name ’Analytical Axial Stress’
plot create plot ’Table2’
pl table 3 style mark 4 5 style mark 6 line style dash &
7 style mark 8 line style dot &
xaxis label ’Radius’ yaxis label ’Stress’
;
ret
Numerical and analytical results are compared in Figures 1.5 and 1.6. The figures show plots of
tables for temperature and stress distributions through the cylinder at steady state. The plotted
values are normalized: temperature is normalized by dividing by Ta ; and stress is normalized by
dividing by mGTa . Figure 1.5 shows the temperature distribution at steady state for the numerical
and analytical solutions. The agreement is very good, with an error of less than 0.1%. A comparison
of results for radial, tangential and axial stress distributions at steady state is provided in Figure 1.6.
E
pGre00 0 u0 i
07G7r
pGrelG,lrG0-8,G8S70y
r0,
p0y
dr l t
E
pGre004 0 u0 i
07G7r
pGrelG,lrG0-8,G8S-0R
,0,yT40 0
.0R 0 0
S0,yT40 0
e0R 0 0
70,yT40R 0
-0R 0R 0
Figure 1.6 Radial stress distribution at steady state for heating of a hollow
cylinder
An infinite line heat source with a constant heat-generating rate is located in an infinite elastic
medium with constant thermal properties. Nowacki (1962) provides the solution to this problem
for the transient values of temperature, radial and tangential stress and radial displacement:
T 1
= E1 (ξ ) (1.63)
a 4π
σr 1 1 − e−ξ
= E1 (ξ ) + (1.64)
bG −4π ξ
σt 1 1 − e−ξ
= E1 (ξ ) − (1.65)
bG −4π ξ
ur 1 1 − e−ξ
= r E1 (ξ ) + (1.66)
bL 8π ξ
r2
where: ξ = 4κt ;
r = radial distance to the line source;
k
κ = ρCp ;
a = qk ;
9K
b = αt a 3K+4G ;
L = unit length; and
∞ −u
E1 (ξ ) = ξ e u du is the exponential integral.
The material properties and initial and boundary conditions for this example are defined:
Material Properties
density (ρ) 2000 kg/m3
shear modulus (G) 30 GPa
bulk modulus (K) 50 GPa
specific heat (Cp ) 1000 J/kg ◦ C
thermal conductivity (k) 4 W/m ◦ C
linear thermal expansion coefficient (αt ) 5 × 10−6 /◦ C
Initial/Boundary Conditions
initial uniform temperature 0◦ C
initial stress state no stresses
It is assumed that the material properties are temperature-independent, the thermal output of the
source is constant (no decay), and the line heat source is of infinite length.
1
The 3DEC grid for this problem is a radial sector of 192 of a cylindrical disk. The axis of the line
heat source coincides with the y-axis of the model. The grid is radially graded in the xz-plane; a
single layer of zones is used in the y-direction. The model is shown in Figure 1.7. A FISH function
(fxface.fis) is used to apply the roller boundaries at the inclined boundary.
E
pGre00 0 u0 i
03,GGG
pGrelG,lrG03:,p:p,0
:0
The line heat source is represented in 3DEC by a series of point sources. The intensity of the
point-source strength, qp , is adjusted to produce an equivalent intensity to the line heat source, Q.
For a quarter-symmetry model and one zone thickness, the relation between qp and Q is
QL 1
qp = · (1.67)
2 192
E
pGre00) 0 u0 i
03,GGG
pGrelG,lrG03:,p:p,0d
rr0,y)0 0c MrGGC
rp0 0 0c MrGGC
dlbrl
E
pGre00 0 u0 i
03,GGG
pGrelrrlGp0.4GG4GG0!
80 0 0
30 0 0
d d-te-
d-tT
E
pGre00 0 u0 i
03,GGG
pGrelG,lrG03:,p:p.0y
,0,0 0
.0 0 0
S0,0 0
e0 0 0
dlbel
dlb
end
@gcons
;
def genb
xc = 0.0
zc = 0.0
ya = 0.0
yb = y1
dang = pi / (float(n2) * 2.0)
;
rr = 1.0
dr = 1.0
loop i (2,n1)
dr = dr * 1.1
rr = rr + dr
endloop
dr1 = 500.0 / rr
;
loop i (1, n1)
i1 = i
ang = 0.0
angn = dang
if i = 1
rr = 0.0
dr = dr1
else
rr = rrn
dr = dr * 1.1
endif
rrn = rr + dr
x1 = xc + rr * cos(ang)
z1 = zc + rr * sin(ang)
x2 = xc + rrn * cos(ang)
z2 = zc + rrn * sin(ang)
x3 = xc + rrn * cos(angn)
z3 = zc + rrn * sin(angn)
x4 = xc + rr * cos(angn)
z4 = zc + rr * sin(angn)
if i = 1
command
poly prism &
a @x1 @ya @z1 @x2 @ya @z2 @x3 @ya @z3 &
b @x1 @yb @z1 @x2 @yb @z2 @x3 @yb @z3
endcommand
else
command
save ils_year0.3dsav
rest ils_year0.3dsav
;
; apply roller boundaries on plane at angle DANG
call fxface.3dfis
def fxnnn
ang = dang
fxfacenx = -sin(dang)
fxfaceny = 0.0
fxfacenz = cos(dang)
end
@fxnnn
;
set fishcall 10 @fxface
set mech on
set thermal off
step 20000
save ils_year1.3dsav
;
def cons
c_g = 3.e10 ; shear modulus
c_k = 5.e10 ; bulk modulus
c_al = 5.e-6 ; coefficient of thermal expansion
c_tk = 4. ; conductivity
c_cp = 1e3 ; specific heat
q_q = 1600. ; line source intensity
c_density = 2.e3
kappa = c_tk / (c_density * c_cp)
end_if
gp_ = gp_next(gp_)
end_loop
ib_ = b_next(ib_)
end_loop
nn = 48
end
def ana_soltu
ib_ = block_head
loop while ib_ # 0
gp_ = b_gp(ib_)
loop while gp_ # 0
rad = sqrt(gp_y(gp_)ˆ2 + gp_z(gp_)ˆ2)
if rad < 1.e-4 then
x = gp_x(gp_)
if x = 0.0 then
table(tab2,x) = 0.0
table(tab8,x) = 0.0
else
e_val = x * x * o4c / tv
val = exp_int
table(tab2,x) = val * o4p
table(tab8,x) = (val + (1.-exp(-e_val))/e_val) * x * o8p
if x < 100.0
table(tab12,x) = table(tab2,x)
endif
; nn = nn + 1
end_if
end_if
gp_ = gp_next(gp_)
end_loop
ib_ = b_next(ib_)
end_loop
end
@num_soltu
@ana_soltu
@num_solst
@ana_solst
save ils.3dsav
;
table 1 name ’3DEC temperature’
table 2 name ’analytical temperature’
table 3 name ’3DEC radial stress’
table 4 name ’analytical radial stress’
end_if
end
vny = face_ny(fa_)
vnz = face_nz(fa_)
dd = vnx*fxfacenx + vny*fxfaceny + vnz*fxfacenz
dd1 = 1.0 - dd
if dd1 < fxfacetol
loop j (1,3)
gp_ = face_gp(fa_,j)
;
gpvx = gp_xvel(gp_)
gpvy = gp_yvel(gp_)
gpvz = gp_zvel(gp_)
dd = gpvx*fxfacenx + gpvy*fxfaceny + gpvz*fxfacenz
gp_xvel(gp_) = gp_xvel(gp_) - dd * fxfacenx
gp_yvel(gp_) = gp_yvel(gp_) - dd * fxfaceny
gp_zvel(gp_) = gp_zvel(gp_) - dd * fxfacenz
;
endloop
endif
fa_ = face_next(fa_)
end_loop
end
; =========================================================
This validation problem illustrates the effect of forward and backward convection, at steady-state,
on the conduction temperature profile in a saturated plane sheet. The saturated sheet has thickness
a, and length L, and is modeled in 3DEC by a horizontal fracture of aperture a. The fracture is
located in a block that is kept adiabatic for the simulation. A system of axes is defined, with the
x-axis pointing to the right along the fracture plane, and the origin located on the left boundary of
the model. The temperature boundary conditions correspond to fixed temperatures: T1 at x = 0, and
T2 at x = L. There is a constant fluid flux in the sheet, corresponding to a fixed boundary pressure
P1 to the left, and P2 to the right of the model.
At steady state, the one-dimensional form of the energy balance equation may be expressed as
d 2 T̂ d T̂
2
− Pe =0 (1.68)
d x̂ d x̂
where T̂ = (T − T1 )/T1 , x̂ = x/L, and Pe is the problem Peclet number. The Peclet number is
defined as
Pe = ρf Cf qL/kfT (1.69)
where ρf is fluid density, Cf is fluid specific heat, q is fluid specific discharge and kfT is fluid
thermal conductivity. The specific discharge is given by the formula
a 2 dp
q=− (1.70)
12μ dx
ePe x̂ − 1
T̂ = T̂2 (1.71)
ePe − 1
The values for the model properties and parameters are listed in Table 1.4:
The Peclet number for the simulation has a magnitude of 1.17. The thermal flux is oriented to
the right. Two simulations are carried out: case 1 is for convection and conduction acting in the
same direction (flow to the right); and case 2 for flow to the left (convection and conduction acting
in opposite directions). The 3DEC temperature profile for steady-state conduction-convection is
compared to the analytical solution in Figures 1.11 and 1.12. As may be seen from these figures,
the match is very good for this particular simulation. The steady-state conduction solution (straight
line) is also sketched on the plots. Compared to the conduction solution, the effect of fluid flow is
seen to displace the temperature curves in the direction of the flow.
E
pGre00 0 u0 i
0rGGG
pGrelG,lrG03:,p:pS0d
rG0,
rr0d
rp00
cb2l
E
pGre00 0 u0 i
0rGGG
3gp,gpGre0r:p,:.30!
rG0,
rr0y
rp00
cb2l
endif
_anat = _t1+(_t2-_t1)*term
table(11,_x) = _anat
table(12,_x) = _t1+(_t2-_t1)*_xhat
end_if
pnt = xmem(pnt)
end_loop
end
@_ftemp
table 10 name ’3DEC’
Table 11 name ’Analytic’
Table 12 name ’Steady State’
title
1D conduction/convection
plot table 10 style mark 11 12 xaxis label ’Location’ &
yaxis label ’Temperature’
ret
The stability of the numerical convection algorithm implemented in 3DEC can be appreciated by
exercising the data file for various values of the grid Peclet and Courant numbers. For example,
by varying grid density and specific discharge, the values of grid Peclet and Courant numbers for
which the algorithm is stable may be derived (see Section 1.1.3.4 on numerical stability).
The 3DEC thermal-fluid coupling algorithm was used to analyze thermal conduction and convection
in a single fracture at depth.
The properties were
Rock
Density 2700 kg/m3
Bulk modulus 6 × 109 Pa
Shear modulus 4 × 109 Pa
Thermal conductivity 2.8 W/m/◦ C
Specific heat 800 J/kg/◦ C
Fracture
Normal stiffness 1 × 109 Pa/m
Shear stiffness 1 × 109 Pa/m
Aperture 2 × 10−5 m
Fluid
Density 1000 kg/m3
Bulk modulus 2 × 109 Pa
Viscosity 1 × 10−3 Pa.sec
Thermal conductivity 0.6 W/m/◦ C
Specific heat 4200 J/kg/◦ C
A cubic domain of 0.5 × 0.5 × 0.5 m was modeled, with the fracture placed on the plane at y = 0.
The block geometry is shown in Figure 1.13, and the block discretization in Figure 1.14.
E
pGre00 0 u0 i
0erSGG
pGrelrGlpS0rr3pG3r,0
30
r
E
pGre00 0 u0 i
0erSGG
pGrelrGlpS0rr3pG3r,0
30
r
In the first stage, the in situ conditions were applied without fluid flow or temperature gradients. In
the second stage, fluid flow at the in situ temperature was simulated by applying an inflow into the
fracture of 2 × 10−7 m3 /s/m at the boundary x = −0.25 m, and an equal outflow at the boundary
x = 0.25 m. These boundary conditions produced a uniform flow velocity of 0.01 m/s in the fracture,
as shown in Figure 1.15. The corresponding fluid pressure distribution is presented in Figure 1.16.
E
pGre00 0 u0 i
0erSGG
pGrelG,lrG0-8,,8r509
u
80GiGrGGre3
80rier3G5
E
pGre00 0 u0 i
0erSGG
pGrelrGlpS0rr4p.4r508
i!i
SiG5SGnGe
SiGeGGnGe
SiG.GGnGe
SiGpGGnGe
SiGGGGnGe
.i:-GGnGe
.i:eGGnGe
.i:.GGnGe
.i:pSGnGe
In the third stage, the injection fluid temperature was fixed at 50◦ C at the boundary x = −0.25 m.
At the outflow boundary, the fluid temperature was left free. The rock temperature was fixed at
190◦ C on all of the model outer boundaries.
Three values were considered for the heat transfer coefficient, h:
(i) Case 1: h = 1 W/m2 /◦ C
In this case, the fluid temperature is progressively reduced along the fracture, reaching a value
of about 150◦ C at the outflow boundary. The final distribution of fluid temperatures is shown in
Figure 1.17.
E
pGre00 0 u0 i
0erSGG
pGrelG,lrG0-8,,8r-0+
ri.e-5nGp
ri.SGGnGp
ri.GGGnGp
ri,SGGnGp
ri,GGGnGp
ripSGGnGp
ripGGGnGp
rirSGGnGp
rirGGGnGp
riGSGGnGp
riGGGGnGp
3iSGGGnGr
3iGGGGnGr
-iSGGGnGr
-iGGGGnGr
5iSGGGnGr
5iGGGGnGr
eiSGGGnGr
eiGGGGnGr
SiSGGGnGr
SiGGGGnGr
The history of fluid temperature at 5 points along the fracture is presented in Figure 1.18. The rock
temperature histories at 3 locations along the fracture are shown in Figure 1.19 (on a different time
scale). The maximum drop of rock temperature is about 4◦ C near the center of the joint, as the rock
temperature is kept fixed at 190◦ C on the model boundaries, including the fluid outflow boundary
at x = 0.25 m.
E
pGre00 0 u0 i
0erSGG
pGrelG,lrG0-8,,8r-04
"
rG0" 00#0M0lGirS
rr0" 00#0M0lGirG
rp0" 00#0M00GiGG
r,0" 00#0M0lGirG
r.0" 00#0M0lGipS
$i0
e2ao
E
pGre00 0 u0 i
0erSGG
pGrelG,lrG0-8,,8r-04
"
p0" 00#0M0lGir
,0" 00#0M00GiG
.0" 00#0M00Gir
$i0
e e2uH
E
pGre00 0 u0 i
0erSGG
pGrelG,lrG0-8,,8r307
ri-33pnGp
ri-GGGnGp
ri5GGGnGp
rieGGGnGp
riSGGGnGp
ri.GGGnGp
ri,GGGnGp
ripGGGnGp
rirGGGnGp
riGGGGnGp
3iGGGGnGr
-iGGGGnGr
5iGGGGnGr
eiGGGGnGr
SiGGGGnGr
The history of fluid temperature at the same 5 points along the fracture is presented in Figure 1.21.
The rock temperature histories at 3 locations along the fracture are plotted in Figure 1.22.
E
pGre00v 0 u0 i
0erSGG
pGrelG,lrG0-8,,8r30=
"
rG0" 00#0A0lGirS
rr0" 00#0A0lGirG
rp0" 00#0A00GiGG
r,0" 00#0A0lGirG
r.0" 00#0A0lGipS
$i0
e2ao
E
pGre00v 0 u0 i
0erSGG
pGrelG,lrG0-8,,8r30=
"
p0" 00#0A0lGir
,0" 00#0A00GiG
.0" 00#0A00Gir
$i0
e e2uH
E
pGre00 0 u0 i
0erSGG
pGrelG,lrG0-8,,8r307
ri3GGGnGp
ri-GGGnGp
ri5GGGnGp
rieGGGnGp
riSGGGnGp
ri.GGGnGp
ri,GGGnGp
ripGGGnGp
rirGGGnGp
riGGGGnGp
3iGGGGnGr
-iGGGGnGr
5iGGGGnGr
eiGGGGnGr
SiGGGGnGr
The history of fluid temperature at 5 points along the fracture is presented in Figure 1.24, and
the rock temperature histories at 3 locations are presented in Figure 1.25. (Note that the rock
temperature is kept at 190◦ C on all model boundaries.)
E
pGre00 0 u0 i
0erSGG
pGrelG,lrG0-8,,8pG04
"
rG0" 00#0M0lGirS
rr0" 00#0M0lGirG
rp0" 00#0M00GiGG
r,0" 00#0M0lGirG
r.0" 00#0M0lGipS
$i0
e2ao
E
pGre00 0 u0 i
0erSGG
pGrelG,lrG0-8,,8pG04
"
p0" 00#0M0lGir
,0" 00#0M00GiG
.0" 00#0M00Gir
$i0
e e2uH
;
; --- thermal fluid coupling --- h = 100 ---
rest thf4b.3dsav
title
fluid/rock heat transfer h = 100.0
; (thermal fluid boundary conditions)
; fix temperature at x=-0.250
bou fluidtemp 50 range x -.25
; heat transfer coefficient h
fluid fht_coe 100.0
; run only thermal analysis
; (keep flow and stresses unchanged)
set flow off therm on mech off
; set fluid-rock thermal interaction
set convection on
cycle 60000
save thf4_h100.3dsav
;
plot current plot Temp
plot clear
1.2.1 Introduction
3DEC is a three-dimensional distinct element code for the analysis of problems in solid mechanics.
A simple thermomechanical option has been added to the code. The main purpose of this option is to
model heat sources embedded in the material. Point heat sources are placed either individually, or in
lines and grids, to represent point, line or plane sources. The stress changes caused by temperature
changes (due to thermal expansion) are equilibrated by using the standard mechanical 3DEC code
to take a “snapshot” of the transient stresses at different times.
Note that the thermal option in 3DEC does not provide a general-purpose heat transfer program. Its
purpose is specifically oriented to solving design problems associated with nuclear waste disposal.
There are several advantages to this method of calculating the temperatures:
(1) The infinite thermal boundary is automatically incorporated.
(2) The calculation procedure is extremely quick, compared with finite element
or other numerical methods.
(3) The calculation of the temperature at any time is independent of the temperature
at previous times, so that it takes the same amount of time to calculate the
solution regardless of the time at which the solution is required.
(4) The mechanical boundary conditions can be applied correctly, because the
standard mechanical logic in 3DEC is used for the mechanical part of the
calculation.
(5) Inhomogeneous and anisotropic mechanical properties can be used.
There are two limitations:
(1) The material is assumed to be thermally homogeneous and isotropic with
constant properties.
(2) Only a few restricted thermal boundary conditions can be applied (i.e., adia-
batic and isothermal planes).
In spite of the limitations of this method of solving thermal-mechanical problems, the thermal
coupling in 3DEC is still a powerful aid when solving problems with heat sources, such as nuclear
waste isolation problems.
1.2.2 Theory
The basic theoretical components of thermomechanical computations in 3DEC are the determination
of the temperature distribution due to the heat sources, and the calculation of the stresses induced
by these temperature changes.
1.2.2.1 The Temperature Due to a Distribution of Heat Sources
The simplest possible heat source to consider is the instantaneous point source in an infinite medium.
The temperature distribution resulting from such a source is given by the simple expression (Carslaw
and Jaeger 1959)
−r 2
Q exp 4Kt
T = 3/2
(1.72)
8(πKt) ρ Cp
The temperature at a point due to any spatial and temporal distribution of heat sources can be found
by considering the distribution as a sum of point sources. In the limit, as the number of sources
becomes large, and the strength of each is reduced so that the total heat output is constant, the sum
becomes an integral over time, space (a line, surface or volume) or both.
For example, the temperature due to a time-dependent heat source of strength Q(t) can be written
as
Q(t ) e−r /4K(t − t )
t 2
1
T = dt (1.73)
8(πK)3/2 ρ Cp o (t − t )3/2
For the case that the heat source is constant (i.e., Q(t) = H (t), where H (t) is the Heaviside step
function, Eq. (1.73) becomes (Carslaw and Jaeger 1959))
1 r
T = erfc (1.74)
4π Kr (4Kt)1/2
Qe −r 2 ir
T = exp · Re w(At)1/2 + (1.75)
4πKr 4Kt (4Kt)1/2
2iz
inf e−t
2
where:w(z) = π 0 z2 −t 2
dt, I m{z} > 0 ;
Re{z} = real part of z;
I m{z} = imaginary part of z; and
i = (−1)1/2 .
Line, plane and volume sources can be written as integrals of a point source over a distance, an
area or a volume, respectively. In this way, completely arbitrary distributions of heat sources can
be represented by analytical expressions. In 3DEC, the integration is carried out numerically, so
that line sources are represented as a series of point sources along a line segment, and plane sources
are represented by a grid of point sources. Thus, to determine the temperature at any point at any
instant in time, a sum of terms must be computed, where each term is of the form of the right-hand
side of Eq. (1.74) or Eq. (1.75).
1.2.2.2 Boundary Conditions
Because the point source solution used in 3DEC applies to an infinite domain, only a limited number
of boundary conditions can be specified. Orthogonal isothermal planes and planes of symmetry
can be specified. Symmetric planes are effected by applying an image source of equal magnitude
on the opposite side of the plane of symmetry, while isothermal planes require an image source
of opposite sign (a heat sink). The effect of having two equal sources is that each produces a flux
across the plane, but the fluxes are in opposite directions so that the net flux is zero. Thus, symmetry
planes are equivalent to adiabatic boundaries. If the sources are of opposite sign, on the other hand,
the net flux is from the heat source to the heat sink. In this case, the two sources each contribute to
the temperature at any point on the plane, but because the two contributions are equal in magnitude
but of opposite sign, the effect is a net temperature change of zero (i.e., an isothermal condition).
Figure 1.26 illustrates the application of adiabatic and isothermal conditions on different faces,
and shows the heat source combinations required to achieve these conditions. Each boundary plane
specified doubles the amount of calculation required, so that if three different planes are defined,
eight times as much calculation is needed.
Adiabatic
Isothermal
Adiabatic
=
Isothermal
Isothermal
=
Adiabatic
Figure 1.26 Adiabatic and isothermal boundaries and their 3DEC represen-
tations
V = βT (1.76)
where:β = volumetric coefficient of thermal expansion;
= 3α, where α is the linear expansion coefficient; and
T = temperature change.
If this volume is then recompressed back to its original size, the mean stress change will be
The negative sign reflects the convention that compressive stresses are negative.
The thermally induced stresses need not be in equilibrium since they depend solely on temperature
changes. If they are not in equilibrium, the medium will move to restore equilibrium. For example,
if two thin metal plates are bolted together but separated by a thin insulator and one is heated, the
thermally induced stresses in the heated plate will cause it to expand, producing bending as the
other plate resists the motion. One plate will be in tension, and the other in compression.
Heat source characteristics such as number of components, decay constants and proportion of
each component are input by the user. In addition, the location and initial strengths of arrays of
point sources are input. Symmetric and isothermal planes may be specified, as may the initial
temperature and the time at which the solution is desired. Once all these quantities have been input,
the temperatures are calculated at all gridpoints. Stresses are not calculated at this stage, but when
mechanical cycling is performed, the average temperature change for each zone is used to calculate
a stress change in each zone. The stress state thus produced will not be in equilibrium, but is treated
as an initial condition for the mechanical calculation. The standard 3DEC mechanical calculation
is used to reach equilibrium.
The commands required for solving thermal-mechanical problems with 3DEC are given below. The
commands are presented in groups corresponding to their function. It is assumed that the reader is
familiar with 3DEC and the mechanical calculations it performs, as well as the input commands.
1.2.4.1 Preparing to Solve a Thermal Problem
THERMAL This command must be input before any polygons are specified. It ensures that
room is allocated in 3DEC program memory for the temperatures and other thermal
variables.
The following three commands are used to define source locations, strengths and starting times:
GRID x1 y1 z1 x2 y2 z2 x3 y3 z3 type n str s N12 n12 N23 n23 <tstart ts>
A grid of point sources is set up, each of strength str. The source is of type n and starts
at time ts. The grid need not be rectangular. The points (x1,y1,z1), (x2,y2,z2) and
(x3,y3,z3) are the top-left, top-right and bottom-right points of the grid, as illustrated
in Figure 1.27. The number n12 defines the number of points across the grid from
left to right, and n23 defines the number from top to bottom.
Once again, the source may be outside the 3DEC grid. The start time defaults to
zero if not specified, but the type and strength must be specified. The type, strength
and number of sources may be specified in any order.
LINE xb yb zb xe ye ze type n strength s n12 np <tstart ts>
A “line” source consisting of np point sources, each of strength str is placed along
the line beginning at xb,yb,zb and ending at xe,ye,ze. The source is of type n, and
starts at time ts. As with the point source, the source may be outside the 3DEC grid.
The start time defaults to zero if not specified, but the type and strength must be
specified. The type, strength and number of sources may be specified in any order.
POINT x y z type n strength s <tstart ts>
A point source of initial strength str is placed at point x,y,z, which need not be within
the blocks defined by 3DEC. The source is of type n (with components and decay
constants defined by the SOURCE command) and starts at time ts. The start time
defaults to zero if not specified, but the type and strength must be specified. The
type and strength may be specified in any order.
n !=6
n = 5
Figure 1.27 A grid of heat sources, input with the command
GRID (10, 10, 10) (20, 20, 20) (20, 0, 10) &
TYP=1 STR=10 N12=5 N23=6
ISOTHERM <x><y><z>
This command is used to define isothermal boundary planes. For example, the two
planes defined by x = 0 and y = 0 are isothermal if the command ISOTHERM x y is
issued.
NOTE: A plane cannot simultaneously be an isothermal boundary and a symmetry
plane.
1.2.4.4 Thermal Properties
PROPERTY The thermal diffusivity, conductivity and linear expansion coefficient can be input
using the keywords diffus, cond and thexp. The thermal diffusivity and conductivity
apply to all materials, but the expansion coefficient can be different for different
materials. For this reason, a material number must be given for the latter (see the
PROPERTY command in Section 1 in the Command Reference).
For example, the commands
prop cond=2.1 diffus=2.3e-6
prop mat=1 thexp=3e-6
prop mat=2 thexp=2.7e-6
assign the conductivity and diffusivity of the entire rock mass as 2.1 W/m ◦ C and
2.3e-6 m2 s, but assign different thermal expansion coefficients to different materials.
1.2.4.5 Calculating the Solution
TIME t
The time at which thermal results are to be calculated is input.
THSOLVE The thermal solution is calculated when this command is issued. Temperatures are
calculated at all gridpoints in the model. Stress increments are not calculated at this
point, but are calculated for all zones the next time mechanical cycling is performed.
If the THSOLVE command is not issued but the program detects a change in con-
ditions (such as source parameters or time) when cycling is next performed, it will
automatically update the thermal solution before calculating new stress increments.
1.2.4.6 Output
Temperature contours may be plotted from the vector/contour menu in the graphics mode.
LIST heat
This command provides information about the heat sources, and other information
specific to the thermal option.
1.2.5 Verification
First, the temperature distribution due to a single point source is compared with an analytical
solution. Then, to test the thermal logic thoroughly, different combinations of point, line and grid
sources are examined. These results should be expected to agree extremely well with the analytical
solutions, since 3DEC is essentially just evaluating the analytical solutions. These first examples
merely confirm that the coding is correct, and illustrate the use of the program.
Finally, the correctness of the mechanical solution is demonstrated by comparing the stress state
due to a point source with that of a point source in an infinite medium. This is a true test of the
thermomechanical coupling, since the 3DEC mechanical logic is used to equilibrate the stresses.
This test is particularly rigorous because the analytical solution contains a singularity at the source
itself, producing steep gradients around the source.
1.2.5.1 A Single Non-decaying Point Heat Source
Consider a 30 kW heat source, located at the origin, in an infinite medium with the following
properties:
Figure 1.28 shows the results of the 3DEC calculations compared with the analytical solution, which
is obtained using Eq. (1.74). 3DEC uses an analytical expression for the temperature calculations,
so it would be expected that the temperatures would agree exactly. The slight differences are due to
the methods of calculating the complementary error function, erfc(x). This function is not readily
available in computer languages, but must be approximated by polynomials. The polynomial used
to evaluate the analytical solution differs from that used in 3DEC. The difference is less than one
percent everywhere, except at the last few points for a time of 2e6 seconds, where the temperatures
are so small that the difference is not significant. Note that temperatures are only shown at a distance
of 2 m from the source because, if placed closer than this, the singularity at the source location
would dominate the figure.
500
t = 2e6
400
Temperature (ºC) t = 2e7
300 t = 2e11
200
100
0
0 2 4 6 8 10 12
Distance from Source
This data file illustrates the superposition of four different sources, with different strengths and
starting times. Once again, the solution is compared with the analytic solution (Table 1.5), and,
again, the 3DEC results are in excellent agreement. The largest errors are located where only small
temperature changes occur.
;
point 2 2 10 typ 1 stre 7000 tstart 1e7
point 1 2 10 typ 2 stre 10000 tstart 1e7
point 2 1 10 typ 1 stre 4000 tstart 2e7
point 2 2 5 typ 2 stre 9000 tstart 3e7
;
time 8e6 ; t less than zero for all sources
thso
list temp
set log on
time 1.2e7 ; t=2e6 for two sources
thso
list temp
time 3e7 ; t=2e7 for two,1e7 for one, 0 for other
thso
list temp
set log off
ret
To verify the operation of the commands to generate lines and grids of sources, one problem is
solved with three different sets of data. The first file uses points only, the second uses lines, and the
third uses grids. Note how the use of the grid and line commands drastically reduces the amount
of input required. It can easily be verified that these three data sets generate the same temperature
distribution.
;
point 0.5 1.5 1.0 typ 1 stre 3000 tstart 1e7
point 0.5 1.5 1.0 typ 2 stre 4000 tstart 2e7
point 0.5 1.5 1.0 typ 3 stre 5000 tstart 0
;
point 0.5 0.5 1.5 typ 1 stre 3000 tstart 1e7
point 0.5 0.5 1.5 typ 2 stre 4000 tstart 2e7
point 0.5 0.5 1.5 typ 3 stre 5000 tstart 0
;
point 0.5 1.0 1.5 typ 1 stre 3000 tstart 1e7
point 0.5 1.0 1.5 typ 2 stre 4000 tstart 2e7
point 0.5 1.0 1.5 typ 3 stre 5000 tstart 0
;
point 0.5 1.5 1.5 typ 1 stre 3000 tstart 1e7
point 0.5 1.5 1.5 typ 2 stre 4000 tstart 2e7
point 0.5 1.5 1.5 typ 3 stre 5000 tstart 0
;
point 1.0 0.5 0.5 typ 1 stre 3000 tstart 1e7
point 1.0 0.5 0.5 typ 2 stre 4000 tstart 2e7
point 1.0 0.5 0.5 typ 3 stre 5000 tstart 0
;
point 1.0 1.0 0.5 typ 1 stre 3000 tstart 1e7
point 1.0 1.0 0.5 typ 2 stre 4000 tstart 2e7
point 1.0 1.0 0.5 typ 3 stre 5000 tstart 0
;
point 1.0 1.5 0.5 typ 1 stre 3000 tstart 1e7
point 1.0 1.5 0.5 typ 2 stre 4000 tstart 2e7
point 1.0 1.5 0.5 typ 3 stre 5000 tstart 0
;
point 1.0 0.5 1.0 typ 1 stre 3000 tstart 1e7
point 1.0 0.5 1.0 typ 2 stre 4000 tstart 2e7
point 1.0 0.5 1.0 typ 3 stre 5000 tstart 0
;
point 1.0 1.0 1.0 typ 1 stre 3000 tstart 1e7
point 1.0 1.0 1.0 typ 2 stre 4000 tstart 2e7
point 1.0 1.0 1.0 typ 3 stre 5000 tstart 0
;
point 1.0 1.5 1.0 typ 1 stre 3000 tstart 1e7
point 1.0 1.5 1.0 typ 2 stre 4000 tstart 2e7
point 1.0 1.5 1.0 typ 3 stre 5000 tstart 0
;
point 1.0 0.5 1.5 typ 1 stre 3000 tstart 1e7
point 1.0 0.5 1.5 typ 2 stre 4000 tstart 2e7
point 1.0 0.5 1.5 typ 3 stre 5000 tstart 0
;
point 1.0 1.0 1.5 typ 1 stre 3000 tstart 1e7
set log on
list temp
ret
line 1.5 1.0 0.5 1.5 1.0 1.5 typ 2 stre 4000 tstart 2e7 n12=3
line 1.5 1.0 0.5 1.5 1.0 1.5 typ 3 stre 5000 tstart 0 n12=3
;
line 0.5 1.5 0.5 0.5 1.5 1.5 typ 1 stre 3000 tstart 1e7 n12=3
line 0.5 1.5 0.5 0.5 1.5 1.5 typ 2 stre 4000 tstart 2e7 n12=3
line 0.5 1.5 0.5 0.5 1.5 1.5 typ 3 stre 5000 tstart 0 n12=3
;
line 1.0 1.5 0.5 1.0 1.5 1.5 typ 1 stre 3000 tstart 1e7 n12=3
line 1.0 1.5 0.5 1.0 1.5 1.5 typ 2 stre 4000 tstart 2e7 n12=3
line 1.0 1.5 0.5 1.0 1.5 1.5 typ 3 stre 5000 tstart 0 n12=3
;
line 1.5 1.5 0.5 1.5 1.5 1.5 typ 1 stre 3000 tstart 1e7 n12=3
line 1.5 1.5 0.5 1.5 1.5 1.5 typ 2 stre 4000 tstart 2e7 n12=3
line 1.5 1.5 0.5 1.5 1.5 1.5 typ 3 stre 5000 tstart 0 n12=3
;
;
time 2.2e7
thso
set log on
list temp
ret
The following data file demonstrates the solution for a decaying heat source:
The resulting temperature distribution is shown in Figure 1.29. The solution was checked against
that given by Eq. (1.75), using tabulated values of w(z) found in Abramowitz and Stegun (1972).
3DEC agrees extremely well with the analytical solution.
14
12
10
Temperature (C)
0
0 2 4 6 8 10 12
Distance From Source (m)
The primary purpose of the thermal logic in 3DEC is to examine thermomechanical effects, not
just thermal effects. To test this capability, the response of an infinite medium to an instantaneous
point source is determined.
Nowacki (1962) gives the solution for an instantaneous point source in an infinite medium as
ui = φ,i
(1.78)
σij = 2G φ,ij − δij φ,kk
mQ K
where:φ(R, t) = 4π R erf (β)1/2
;
R = distance from source;
t = time;
α(1+ν)
m = (1−ν) ;
ν = Poisson’s ratio;
G = shear modulus;
β = 4 Kt;
K = thermal diffusivity;
Q = “source strength”
= Qo /ρ Cp ;
Qo = heat released (J );
ρ = density;
Cp = specific heat; and
α = thermal expansion coefficient.
From these equations, we can obtain the stresses and displacements along the z-axis (x = y = 0):
mG 2z 4 z3
σxx = σyy = −V1 + V2 + V3
2πz3 (π)1/2 (π)1/2
mG 4z
σzz = 2V1 − V2
2π z3 (π)1/2
ux = uy = 0
m V2 V1
uz = 1/2
−
2π z (π) 2z
where:V1 (z, t) = Q erf (β)z1/2 ;
V2 (z, t) = (β)Q1/2 exp(−z2 /β) ; and
V3 (z, t) = V2 (z, t)/β.
The following material properties and source parameters are used to solve this problem with 3DEC:
The total heat input is 108 J , which in 3DEC is represented by a 100 kW source acting for 1000
seconds. A cube of side 5 m is used to represent one octant of the medium, as shown in Figure 1.30.
For the times at which the solution is obtained here (1 day and 3 days), the model is large enough
to reduce the effect of the boundaries.
Left Face
Fixed in
x - Direction
y
z Front Face
Fixed in
x z - Direction
Heat Source
Bottom Face
Fixed in
y - Direction
Figure 1.30 Conceptual representation of 3DEC block for modeling effect of
point heat source
cycle 500
sav tmtest3.3dsav
return
The block plot from 3DEC is shown in Figure 1.31. The discretization is finer near the source,
where stress and displacement gradients are greatest.
Figure 1.32 shows the comparison between σzz calculated by 3DEC, and the analytical solution
near the z-axis (x = y = 0). Figures 1.33 and 1.34 show similar plots for σxx and σyy , while
Figure 1.35 shows the displacements. On each of these plots, two different cases are shown: one
where the outer boundaries are free to move; and another where the outer boundaries are fixed
against movement normal to the boundary. The differences between fixed and free boundaries are
not particularly significant, except for the displacements near the outer edge of the block.
-50
Analytic
-100
Free Boundary
Fixed Boundary
-150
Stress (MPa)
-200
-250
-300
-350
-400
0 1 1 2 2 3 3 4 4 5 5
Distance (m)
-10
Analytic
-20 Free Boundary
Fixed Boundary
-30
Stress (MPa)
-40
-50
-60
-70
-80
0 1 1 2 2 3 3 4 4 5 5
Distance (m)
50
-50
Analytic
Free Boundary
-100
Fixed Boundary
Stress (MPa)
-150
-200
-250
-300
-350
-400
0 1 1 2 2 3 3 4 4 5 5
Distance (m)
50
-50
Analytic
-100 Free Boundary
Fixed Boundary
Stress (MPa)
-150
-200
-250
-300
-350
-400
0 1 1 2 2 3 3 4 4 5 5
Distance (m)
Figure 1.33 Tangential stresses Sxx and Syy along z-axis (x = y = 0) after 1
day
10
-10
Analytic
Free Boundary
-20
Fixed Boundary
Stress (MPa)
-30
-40
-50
-60
-70
-80
0 1 2 3 4 5 6
Distance (m)
10
-10
Analytic
-20 Free Boundary
Fixed Boundary
Stress (MPa)
-30
-40
-50
-60
-70
-80
0 1 1 2 2 3 3 4 4 5 5
Distance (m)
Figure 1.34 Tangential stresses Sxx and Syy along z-axis (x = y = 0) after 3
days
0.09
0.08
0.07
Analytic
Free Boundary
0.06
Fixed Boundary
Displacement (m)
0.05
0.04
0.03
0.02
0.01
0.00
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5
Distance (m)
0.030
0.025
Analytic
Free Boundary
Fixed Boundary
0.020
Displacement (m)
0.015
0.010
0.005
0.000
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5
Distance (m)
Figure 1.36 shows a cross section through a drift in a hypothetical nuclear waste repository. A heat
source in the floor of the drift represents a waste canister. The properties of the host rock are:
The initial stress is hydrostatic at 5 MPa, and the heat source is assumed to have only one component
with a decay constant of 4 × 10−10 s−1 and an initial strength of 30 kW.
For this example, it is assumed that there are other heat sources (spaced 15 m apart along the
length of the drift), and that there are other drifts (spaced 20 m apart), also with heat sources in the
floor. The mechanical boundary conditions incorporate this symmetry, but the thermal symmetry
is simulated by grids of regularly spaced heat sources. All but one member of each grid are located
outside the 3DEC grid.
The data file to solve this problem is given below, and Figure 1.37 shows a 3DEC plot of the block.
Figures 1.38 and 1.39 show temperatures and displacements in the grid.
10 m
10 m
3m
1m
3m
3m
10 m
z
y
x
15 m
E
pGre00 0 u0 i
0pGGG
pGrelG,lrG03:,,:pe0
:0
E
pGre004 0 u0 i
0pGGG
pGrelG,lrG03:,,:p-0x
,iGGGGnG,
pi-SGGnG,
piSGGGnG,
pipSGGnG,
piGGGGnG,
ri-SGGnG,
riSGGGnG,
ripSGGnG,
riGGGGnG,
-iSGGGnGp
SiGGGGnGp
piSGGGnGp
GiGGGGnGG
# m$
:0Gir.p8re
:0Sirrr8-
1.3 References
Abramowitz, M., and I. A. Stegun. Handbook of Mathematical Functions. New York: Dover
Publications Inc. (1972).
Carslaw, H. S., and J. C. Jaeger. Conduction of Heat in Solids, 2nd Ed. Oxford: Clarendon Press
(1959).
Crank, J. The Mathematics of Diffusion, 2nd Ed. Oxford: Oxford University Press (1975).
Crouch, S. L., and W. C. McClain. “Interim Report on Development of a Semi-Empirical Numerical
Model for Simulating the Deformational Behavior of a High-Level Radioactive Waste Repository in
Bedded Salt,” University of Minnesota Report to Oak Ridge National Laboratory, ORNL/TM-5462
(1978).
Dahlquist, G., and A. Bjorck. Numerical Methods. Englewood Cliffs, New Jersey: Prentice-Hall
(1974).
Karlekar, B. V., and R. M. Desmond. Heat Transfer, 2nd Ed. St. Paul: West Publishing Co.
(1982).
Nowacki, W. Thermoelasticity. New York: Addison-Wesley (1962).
Perrochet, P., and D. Berod. “Stability of the Standard Crank-Nicolson-Galerkin Scheme Applied
to the Diffusion-Convection Equation: Some New Insights,” in Water Resources Research, 29(9),
3291-3297 (September 1993).
2 DYNAMIC ANALYSIS
2.1 Introduction
2.2 Damping
As described in Section 1 in Theory and Background, 3DEC uses a dynamic algorithm for
problem solution. Natural dynamic systems contain some degree of damping of the vibrational
energy within the system; otherwise, the system would oscillate indefinitely when subjected to
driving forces. Damping is due, in part, to energy loss as a result of slippage along contacts of
blocks within the system, internal friction loss in the intact material, and any resistance caused by
air or fluids surrounding the structure. 3DEC is used to solve two general classes of mechanical
problems: quasi-static and dynamic. Damping is used in the solution of both classes of problems,
but quasi-static problems require more damping. Two types of damping (mass-proportional and
stiffness-proportional) are available in 3DEC. Mass-proportional damping applies a force which is
proportional to absolute velocity and mass, but in the direction opposite to the velocity. Stiffness-
proportional damping applies a force, which is proportional to the incremental stiffness matrix
multiplied by relative velocities or strain rates, to contacts or stresses in zones. In 3DEC, either
form of damping may be used separately or in combination. The use of both forms of damping
in combination is termed Rayleigh damping (Bathe and Wilson 1976). For solution of quasi-static
problems using finite difference schemes, mass-proportional damping is generally used (Otter et al.
1966). 3DEC allows use of an automatic “adaptive” viscous damping scheme developed by Cundall
(1982) for solution of quasi-static problems. These schemes are discussed in Section 1 in Theory
and Background. For dynamic analyses, either mass-proportional or stiffness-proportional, or
both (i.e., Rayleigh), forms of damping may be used, as described in the next section.
In performing dynamic analysis with any code, it is usually necessary to account for energy losses in
the physical system (e.g., heat, hysteresis) which are not accounted for in the numerical algorithm.
In general, very little damping is used for highly elastic systems, and more damping is used for
geomechanical materials, especially soils.
In the continuum analysis of structures, proportional Rayleigh damping is typically used to damp
the natural oscillation modes of the system. In dynamic finite-element analysis, a damping matrix,
C, is formed with components proportional to the mass (M) and stiffness (K) matrices,
C =α M +β K
For a multiple degrees-of-freedom system, the critical damping ratio, ξi , at any angular frequency
of the system, ωi , can be found from (Bathe and Wilson 1976)
α + β ωi2 = 2 ωi ξi (2.1)
or
1α
ξi = + β ωi (2.2)
2 ωi
The critical damping ratio, ξi , is also known as the fraction of critical damping for mode i with
angular frequency, ωi .
Figure 2.1 shows the variation of the normalized critical damping ratio with angular frequency, ωi .
Three curves are given: mass or stiffness components only, and the sum of both components. As
shown, mass-proportional damping is dominant at lower angular frequency ranges, while stiffness-
proportional damping dominates at higher angular frequencies. The curve representing the sum of
both components reaches a minimum at
ξmin = (α β)1/2
(2.3)
ωmin = (α/β)1/2
or
α = ξmin ωmin
(2.4)
β = ξmin / ωmin
β=0
5
α= 0
4
total
ξ i / ξ min
0
0 5 10 15 20 25 30
ωi
Figure 2.1 Variation of normalized critical damping ratio with angular fre-
quency
The required input parameters to specify Rayleigh damping in 3DEC are fmin (input parameter
freq) and ξmin (input parameter fcrit).
For the case shown in Figure 2.1, ωmin = 10 radians per second, and ξmin = 1. Note that the
damping ratio is almost constant over at least a 3:1 frequency range (e.g., from 5 to 15). Since
damping in geologic media is known to be predominately hysteretic (Gemant and Jackson 1937;
Wegel and Walther 1935), and hence independent of frequency, ωmin is usually chosen to lie in
the center of the range of frequencies present in the numerical simulation. In this way, hysteretic
damping is simulated approximately.
From the preceding equations and Figure 2.1, it can be seen that at frequency ωmin (or fmin ), and
only at that frequency, mass damping and stiffness damping each supply half of the total damping.
In order to demonstrate how Rayleigh damping works in 3DEC, the results of the following four
damping cases involving a block sitting on a fixed block with gravity suddenly applied can be
compared:
(a) undamped;
(b) Rayleigh damping (both mass and stiffness damping);
;no damping
title
Plot of vertical displacement vs time (UNDAMPED)
hist label 2 ’Undamped’
hist label 3 ’Time’
plot create plot Hist
plot his 2 vs 3 xaxis label ’Time’ yaxis label ’Vertical displacement’
cy 250
;
;mass and stiffness
rest damp.3dsav
hist n 20
damp 1 16
title
In the first case, with no damping, a natural frequency of oscillation of approximately 16 cycles/sec
can be observed (see Figure 2.2). The theoretical period of oscillation is given by
1 ka 1/2
frequency = = 15.9 cycles/second
2π m
The problem is critically damped if: a fraction of critical damping = 1 is specified; the natural
frequency of oscillation = 16 is specified; and both mass and stiffness damping are used.
The results in Figure 2.3 show that the problem is critically damped using these parameters. If
only mass (Figure 2.4) or stiffness (Figure 2.5) damping is used, then the damping must be doubled
(since each contributes 1/2 to the Rayleigh damping) to obtain critical damping. Figures 2.4 and
2.5 show, again, that the problem is critically damped using mass and stiffness damping with twice
the damping specified (see Example 2.1).
E oo
o oo oc5,-5a
mnplH
m
mnpln
Figure 2.2 Plot of vertical displacement versus time, for a single block con-
tacting on a rigid base with gravity suddenly applied (no damping)
oITSA
SMANmMTmAMoI1S;1TFou
So oro
noTo
mnplo
m
mnpln
Figure 2.3 Plot of vertical displacement versus time, for a single block con-
tacting on a rigid base with gravity suddenly applied (mass and
stiffness damping)
oS;M
SMA2mMDmAMoE©S)©D;o,
So o-
noDo
mnplo
m
mnpln
Figure 2.4 Plot of vertical displacement versus time, for a single block con-
tacting on a rigid base with gravity suddenly applied (mass damp-
ing only)
oMDIF
ISTDmSFmTSoN2I;2F;ou
Io o,
noFo
mnplH
m
mnpln
Figure 2.5 Plot of vertical displacement versus time, for a single block con-
tacting on a rigid base with gravity suddenly applied (stiffness
damping only)
2.2.4 Guidelines for Selecting Rayleigh Damping Parameters for Dynamic Analysis
Range of Predominant
Frequencies
Velocity
Spectrum
Frequency
If the highest predominant frequency is three times greater than the lowest predominant frequency,
then there is a 3:1 span or range that contains most of the dynamic energy in the spectrum. The
idea in dynamic analysis is to adjust ωmin of the Rayleigh damping so that its 3:1 range coincides
with the range of predominant frequencies in the problem. ξmin is adjusted to coincide with the
correct physical damping ratio. The predominant frequencies are neither the input frequencies nor
the natural modes of the system, but a combination of both. The idea is to try to get the right
damping for the important frequencies in the problem.
For some problems involving large movements of blocks, it is improper to use any mass damping
because the block motion might be artificially restricted. Examples of such problems include any
problems involving free flow or fall of blocks under gravity, and impulsive loading of blocks due
to explosions. In such cases it may be appropriate to use only stiffness-proportional damping.
The stiffness-proportional component of Rayleigh damping does affect the critical timestep for the
explicit solution scheme in 3DEC. The controlling timestep, therefore, may need to be reduced
(using the FRACTION command) as the stiffness damping component is increased (see Belytschko
1983). For problems involving free fall and “bounce” of blocks from a fixed base, the coefficient of
restitution is required for accurate modeling. Onishi et al. (1985) provide a method for estimating
stiffness damping parameters based on the coefficient of restitution.
For many problems, the important frequencies are related to the natural mode of oscillation of the
system. Examples of this type of problem include seismic analysis of surface structures such as
dams or dynamic analysis of underground excavations.
For these problems, the fundamental frequency, f , associated with the natural mode of oscillation
is
C
f = (2.6)
γ
For an elastic continuous system, the speed of propagation, Cp , is given by Cp = [(K+4/3 G)/ρ]1/2
for p-waves, and Cs = (G/ρ)1/2 for s-waves, where K = bulk modulus, G = shear modulus and ρ
= density.
The longest wavelength, characteristic length or fundamental wavelength depends on boundary
conditions. Consider a solid bar of length 1 with boundary conditions, as shown in Figure 2.7(a).
The fundamental mode shapes for cases (1), (2) and (3) are as shown in Figure 2.7(b).
If shear motion of the bar gives rise to the lowest natural mode, then Cs is used in the preceding
equation; otherwise, Cp is used if motion parallel to the axis of the bar gives rise to the lowest
natural mode.
In the limit of very high joint stiffness, an assemblage of blocks should resemble a continuum,
both statically and dynamically. Consider the problem of eight square deformable blocks resting
on a rigid base. Three problems can be treated: an unconfined column; a confined column in
compression; and a column in shear.
The column is loaded by applying gravity in either the x- or z-direction. For the dynamic case, the
mass damping is zero, with stiffness-proportional damping as follows:
The case of confined compression is modeled by inhibiting lateral displacement along the vertical
boundaries, which prevents lateral deformation of the blocks. For unconfined compression, lateral
displacement is not inhibited. For the column in shear, vertical motion along all boundaries is
inhibited. Other properties are:
Material Properties
bulk modulus K = 1.5 × 104 | for compression tests
|
shear modulus G = 0.428562 × 104 |
|
Poisson’s ratio 0.4 |
density ρ = 1.0
Table 2.2 compares the theoretical periods and calculated (3DEC) natural periods of oscillation.
The theoretical values for natural period of oscillation are calculated as
natural period, T = 4L (ρ / E ∗ )
def _time
_time = time
end
def z_dis
z_d = gp_zdis(i_vert)
if z_d < v_min then
v_min = z_d
d_t = time
endif
z_dis = z_d
end
def period
ii = out(’Period of oscillation = ’ + string(d_t*2.0))
end
hist n 100 @z_dis @_time
mscale off
damp 0.1 1.0 stiff
hist label 1 ’Vertical Displacement’
hist label 2 ’Time’
solve time 11
@period
plot create plot Hist
plot hist 1 vs 2 xaxis label ’Time’ yaxis label ’Displacement’
ret
def _time
_time = time
end
def z_dis
z_d = gp_zdis(i_vert)
if z_d < v_min then
v_min = z_d
d_t = time
endif
z_dis = z_d
end
def period
ii = out(’Period of oscillation = ’ + string(d_t*2.0))
end
hist n 100 @z_dis @_time
mscale off
damp 0.1 1.0 stiff
hist label 1 ’Vertical Displacement’
hist label 2 ’Time’
solve time 15
@period
plot create plot Hist
plot hist 1 vs 2 xaxis label ’Time’ yaxis label ’Displacement’
ret
The physical stiffness of joints in-situ can have a substantial influence on seismic wave propagation.
Myer et al. (1990) present field and laboratory test results which demonstrate the effect of the
stiffness of dry natural fractures in rock on high frequency attenuation, and changes in travel time
of the seismic wave. It can be important to represent this effect in the discontinuum model if the wave
transmission is be to modeled accurately. However, care must be taken to not introduce a numerical
distortion of the wave which could mask the actual effect of the joints on wave propagation.
Numerical distortion of the propagating wave can occur in a dynamic analysis, whether it is based
on a continuum or discontinuum program, as a function of the modeling conditions. Both the
frequency content of the input wave and the wave speed characteristics of the system will affect the
numerical accuracy of wave transmission. Kuhlemeyer and Lysmer (1973) show that for accurate
representation of wave transmission through a model, the spatial element size must be smaller than
approximately one-tenth to one-eighth of the wavelength associated with the highest frequency
component of the input wave – i.e.,
1
l≤ ·γ (2.7)
8
where γ is the wavelength associated with the highest frequency component for peak velocities
through the medium. For discontinuum codes, this also applies to joint spacing (or block size).
In order to achieve an accurate representation of a stress wave through a distinct element model,
particularly when the joint spacing is variable, the blocks should be made deformable to accom-
modate the element size restriction imposed by frequency and wavelength. This is accomplished
in 3DEC, as discussed in Section 1 in Theory and Background, by subdividing each block into a
mesh of finite-difference elements. These elements are then subject to the Kuhlemeyer and Lysmer
restriction.
The effect of model conditions on numerical distortion of wave transmission is demonstrated by a
simple analysis of a column of blocks subjected to an impulse load applied at the base (Figures 2.8
and 2.10). The block sizes range from 1 m to 5 m (average size of 2 m), the contacts between
blocks have a linearly elastic behavior, and the p-wave speed for the system is 4300 m/sec. A
triangular-shaped impulse load, with a maximum frequency of approximately 200 Hz, is applied
at the base (the solid curve in Figures 2.9 and 2.11). The wavelength associated with the highest
frequency of this system is 21.5 m. Thus, according to Kuhlemeyer and Lysmer, in order to transmit
this wave without distortion, the element size must not exceed approximately 2 m.
A rigid block analysis is done with a constant contact normal stiffness used to produce an average
wave speed of 4300 m/sec (based on the average joint spacing). A highly distorted velocity history is
calculated at the top of the column, as seen by the dashed curve in Figure 2.9. This distortion can be
reduced for this problem by varying the normal stiffness locally to keep the wave speed constant at
contacts between blocks. However, in general, the calculation of effective (local) normal stiffnesses
becomes extremely complex for a multiply jointed system, making this approach impractical.
A deformable block analysis is performed with the maximum size of the finite difference elements
smaller than 2 m (see Figure 2.10). The elastic moduli for the blocks and the contact stiffness are
calculated to produce the given p-wave speed. The distortion in the wave at the top of the column
is now essentially eliminated, as indicated in Figure 2.11. The elastic deformation parameters
represent the physical properties of the blocks and contacts separately in this case, and do not have
to be adjusted locally.
Physically measured values for normal and shear stiffnesses of a geologic structure, such as joints,
faults, bedding planes, etc., are not generally available. It is often necessary to back-calculate
properties based on measured values for the elastic-deformation properties of the intact material
and the wave speed through the jointed system. Formulae relating properties of an equivalent elastic
continuum to properties for intact material and joints are given, for example, by Singh (1973) and
Gerrard (1982). These relations can be used to provide reasonable estimates for joint stiffness
properties in 3DEC, to produce the measured shear and compressional wave speeds of the system.
moHo
moHco
Figure 2.9 Input wave (solid) at base and calculated wave (dashed) at top of
column of rigid block model
Figure 2.10 Column of variably sized blocks subdivided into finite difference
zones
a6©f22
hfu1if©iufa2C©uCfta3
uaTa
ha- a
ca2a
moHo
moHco
Figure 2.11 Input wave (solid) at base and calculated wave (dashed) at top of
column of deformable block model
Example 2.5 Column of variably sized rigid blocks subjected to impulse load at base
new
config dynamic
poly brick 0,10 0,1 0,100
plot create plot Blocks
plot block
plot reset
def varcut
ntot = 36
nc = 1
rat = 1.08
zcut = 5.0
zloc = zcut
loop while nc < ntot
if zloc < 99.0
command
jset dip 0 dd 180 origin 0,0,@zloc
endcommand
endif
if zloc < 50.0 then
zcut = zcut / rat
else
zcut = zcut * rat
endif
zloc = zloc + zcut
nc = nc + 1
endloop
end
@varcut
;
; properties for rigid block model
prop mat=1 dens 1000
prop jmat=1 jkn 10e9 jks 10e9 jcoh 1e10 jten 1e10
prop mat=1 k=10.0e9 g=7.5e9
; impulse load
def i_block
i_block = b_near(5.0,0.5,2.5)
end
@i_block
def pulse
whilestepping
dytime = time
zpulse = vmax / tpeak * dytime
if dytime > tpeak then
zpulse = vmax - (vmax / (tend - tpeak)) * (dytime - tpeak)
endif
if dytime > tend then
zpulse = 0.0
endif
b_zvel(i_block) = zpulse ; velocity assigned to rigid block no. 2
end
set @vmax = 11.0 @tpeak = 0.005 @tend = 0.06
; fix bottom block to apply impulse for rigid block model
fix range x 0,10 y 0,10 z 0,5
; quiet boundary at top for both deformable block model
bound mat 1 zvisc range z 100
bound mat 1
bound xvel 0.0 range x 0.0
bound xvel 0.0 range x 10.0
bound yvel 0.0 range y 0.0
bound yvel 0.0 range y 10.0
; monitor velocities at bottom and top
hist n 10
hist zvel 0,0,0
hist zvel 0,0,95
hist @pulse
hist @dytime
hist label 1 ’Input Wave’
hist label 2 ’Calculated Wave’
hist label 4 ’Time’
; add 5\% stiffness damping
damp 0.05 200 stiff
title
Example problem: dynamic analysis of column shear
plot create plot Hist
plot set jobtitle on
plot hist 1 2 linestyle style dot vs 4 &
xaxis label ’Time’ yaxis Label ’Displacement’
solve time 0.12
save ex2_05.3dsav
ret
Example 2.6 Column of variably sized deformable blocks subjected to impulse load at base
new
config dynamic
poly brick 0,10 0,1 0,100
plot create plot Blocks
plot block
plot reset
def varcut
ntot = 36
nc = 1
rat = 1.08
zcut = 5.0
zloc = zcut
loop while nc < ntot
if zloc < 99.0
command
jset dip 0 dd 180 origin 0,0,@zloc
endcommand
endif
if zloc < 50.0 then
zcut = zcut / rat
else
zcut = zcut * rat
endif
zloc = zloc + zcut
nc = nc + 1
endloop
end
@varcut
;
gen edge 2
; properties for zoned model
prop jmat = 1 dens 1000 jkn 200e9 jks 200e9 jcoh 1e10 jten 1e10
prop mat = 1 k=10.0e9 g=7.5e9
; impulse load
def pulse
whilestepping
dytime = time
zpulse = vmax / tpeak * dytime
if dytime > tpeak then
zpulse = vmax - (vmax / (tend - tpeak)) * (dytime - tpeak)
endif
if dytime > tend then
zpulse = 0.0
endif
pulse = zpulse ; velocity history for zoned model
end
set @vmax = 11.0 @tpeak = 0.005 @tend = 0.06
; velocity boundary for zoned model
bound zvel 1.0 hist @pulse range z -.1 .1
; quiet boundary at top for both deformable block model
bound mat 1 zvisc range z 100
bound mat 1
bound xvel 0.0 range x 0.0
bound xvel 0.0 range x 10.0
bound yvel 0.0 range y 0.0
bound yvel 0.0 range y 10.0
; monitor velocities at bottom and top
hist n 100
hist zvel 0,0,0
hist zvel 0,0,95
hist @pulse
hist @dytime
hist label 1 ’Input Wave’
hist label 2 ’Calculated Wave’
hist label 4 ’Time’
; add 5% stiffness damping
damp 0.05 200 stiff
solve time 0.12
title
Example problem: dynamic analysis of column shear
plot create plot Hist
plot set jobtitle on
For dynamic input with a high peak velocity and short rise-time, the Kuhlemeyer and Lysmer
requirement may necessitate a very fine spatial mesh and a correspondingly small timestep. The
effect is compounded in discontinuum codes because the wave propagation across discontinuities
can produce higher frequency components than are provided in the input wave. The consequence is
that reasonable analyses may be prohibitively time- and memory-consuming, as well as extremely
expensive. In such cases, it may be possible to adjust the input by recognizing that most of the
power for the input history is contained in lower frequency components. By filtering the history and
removing high frequency components, a coarser mesh may be used without significantly affecting
the results.
The filtering procedure can be accomplished with a low-pass filter routine such as the fast Fourier
transform technique. For example, the unfiltered velocity record shown in Figure 2.12 represents
a typical waveform containing a very high frequency spike. The highest frequency of this input
exceeds 50 Hz, but, as shown by the power spectral density plot of Fourier amplitude versus
frequency (Figure 2.13), most of the power (approximately 99%) is made up of components of
frequency 15 Hz or lower. It can be inferred, therefore, that by filtering this velocity history with
a 15 Hz low-pass filter, less than 1% of the power is lost. The input filtered at 15 Hz is shown in
Figure 2.14, and the Fourier amplitudes are plotted in Figure 2.15. The difference in power between
unfiltered and filtered input is less than 1%, while the peak velocity is reduced 38%, and the rise
time is shifted from 0.035 sec to 0.09 sec. Analyses should be performed with input at different
levels of filtering, to evaluate the influence of the filter on model results.
4
Velocity (cm/sec)
3
(Thousands)
-1
0 0.2 0.4
Time (sec)
130
120
110
100
Fourier Amplitude
90
(Times 10E9)
80
70
60
50
40
30
20
10
0
0 2 4 6 8 10 12 14 16 18 20
Frequency
3
2.8
2.6
2.4
2.2
2
Velocity (cm/sec)
1.8
(Thousands)
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
0 0.2 0.4
Time (sec)
Figure 2.14 Filtered velocity history at 15 Hz
130
120
110
100
90
Fourier Amplitude
80
(Times 10E9)
70
60
50
40
30
20
10
0
0 2 4 6 8 10 12 14 16 18 20
Frequency
Density scaling is a technique used in 3DEC in quasi-static calculations that substantially improves
the efficiency in obtaining solutions to large problems. In quasi-static problems, inertial forces are
not important. The gridpoints’ masses can be scaled for optimal numerical convergence without
affecting the solution. In dynamic analyses, however, global scaling cannot be used. Complex
jointed systems often result in very small blocks and/or zones being created during the automatic
meshing procedure. The small blocks/zones require very small timesteps for numerical stability of
the explicit algorithm. This makes some dynamic solutions extremely time-consuming. However,
as these blocks/zones may be very small, with very small masses, it is possible to introduce some
density scaling only for these blocks/zones in such a way that the change of the system inertia is
negligible. This scheme of partial density scaling is implemented in 3DEC in such a way that the
user controls the amount of scaling to be introduced. Given the timestep calculated by the code,
the user specifies the desired timestep with the command MSCALE part dt. This command specifies
that only the amount of density scaling required to achieve the timestep dt is to be applied to the
system. When a CYCLE command is given, a message indicating the number of gridpoint masses
that were scaled, and the amount of additional mass introduced, is printed.
Figure 2.16 shows a simple block system with some low-thickness blocks. The timestep required
for dynamic analysis without any scaling is 1.005e-6 seconds. Using partial density scaling, the
timestep may be increased to 5e-6 seconds, while the total system mass is increased by only 5%.
This information is printed by 3DEC following the use of the MSCALE part command:
no. scaled g.p. masses = 68
min. g.p. scaling factor = 4.038E-02
max. g.p. scaling factor = 1.000E+00
min. g.p. added mass = 0.000E+00
max. g.p. added mass = 3.103E-05
min. block added mass = 0.000E+00
max. block added mass = 1.938E-04
total added mass in model = 8.202E-04
total real mass in model = 1.920E-02
added mass / real mass = 4.272E-02
E
pGre0 0 0! u0 i
0--SS
pGrelG,lrG0.3,r3r90y
30
r
The effect of this amount of partial density scaling was checked by comparing the system response
to a sinusoidal shear applied at the base of the model. A viscous boundary condition was applied
to the top of the model to simulate an extended medium. Figure 2.17 shows the x-velocity applied
at the base, and the velocity obtained at the top of the model, obtained in the run without scaling.
Figure 2.18 shows the same quantities for the run with partial density scaling, with a timestep
about 5 times larger. It can be seen that the wave propagation is not affected by the small amount
of scaling introduced.
m2taH
Figure 2.17 Velocities at the bottom and top of the model, for analysis without
any density scaling
mHtaH
mHtao
Figure 2.18 Velocities at the bottom and top of the model, for analysis with
partial density scaling
config dynamic
poly brick -1 1 -1 1 -1 1
plot create plot Blocks
plot block
plot reset
jset dip 5 dd 45 origin 0 0 0
jset dip 10 dd 40 origin 0 0 0
jset dip 70 dd 95
jset dip 80 dd 10
gen edge 0.5
prop mat 1 dens 0.0024 bulk 33333 shear 20000
prop jmat 1 jkn 500000 jks 500000 jcoh 1e9 jtens 1e9
insitu stress -1e-6 -1e-6 -1e-6 0 0 0
mscale off
damp 0 0 mass
bound yvel 0 zvel 0 range zr -1.0
bound xvel 0.1 hist sin 100 1.0 range zr -1.0
bound xvisc yvisc zvisc range zr 1.0
bound yvel 0 zvel 0 range xr -1.0
bound yvel 0 zvel 0 range xr 1.0
bound yvel 0 zvel 0 range yr -1.0
bound yvel 0 zvel 0 range yr 1.0
bound mat 1
def _time
_time = time
end
hist unbal
hist xvel -1 -1 -1
hist xvel -1 -1 1
hist @_time
hist label 2 ’Velocity at Bottom’
hist label 3 ’Velocity at Top’
hist label 4 ’Time’
plot create plot Hist
The modeling of geomechanics problems involves media which, at the scale of the analysis, are
better represented as unbounded. Deep underground excavations are normally assumed to be
surrounded by an infinite medium, while surface and near-surface structures are assumed to lie on
a half-space. Numerical methods relying on the discretization of a finite region of space require
that appropriate conditions be enforced at the artificial numerical boundaries. In static analyses,
fixed or elastic boundaries (e.g., represented by boundary element techniques) can be realistically
placed at some distance from the region of interest. In dynamic problems, however, such boundary
conditions cause the reflection of outward propagating waves back into the model, and do not allow
the necessary energy radiation. The use of a larger model can minimize the problem, since material
damping will absorb most of the energy in the waves reflected from distant boundaries. However,
this solution leads to large computational costs. The alternative is to use nonreflecting (or absorbing)
boundaries. Several formulations have been proposed. The viscous boundary developed by Lysmer
and Kuhlemeyer (1969) is used in 3DEC. It is based on the use of independent dashpots, and is
nearly totally effective for body waves approaching the boundary at angles of incidence above 30◦ .
For lower angles of incidence, or for surface waves, the energy absorption is only approximate.
However, it has the advantage of being an effective technique which can be used in time-domain
analyses. Its effectiveness has been demonstrated in both finite-element and finite-difference models
(Kunar et al. 1977). A variation of this technique proposed by White et al. (1977) is also widely
used.
More efficient energy absorption (for example, in the case of Rayleigh waves) requires the use of
frequency-dependent dashpots, which can only be used in frequency-domain analyses (e.g., Lysmer
and Waas 1972). These are usually designated as consistent boundaries, and involve the calculation
of dynamic-stiffness matrices coupling all the boundary degrees-of-freedom. Boundary-element
methods may be used to derive these matrices (e.g., Wolf 1985). A comparative study of the
performance of different types of elementary, viscous and consistent boundaries was reported by
Roesset and Ettouney (1977).
A different procedure to obtain efficient absorbing boundaries for use in time-domain studies was
proposed by Cundall et al. (1978). It is based on the superposition of solutions with stress and
velocity boundaries in such a way that reflections are canceled. In practice, it requires adding the
results of two parallel, overlapping grids in a narrow region adjacent to the boundary. This method
has been shown to provide effective energy absorption, but is difficult to implement for a blocky
system with complex geometry, and thus, is not used in 3DEC.
The viscous boundaries proposed by Lysmer and Kuhlemeyer (1969) consist of independent dash-
pots attached to the boundary in the normal and shear directions. They provide viscous normal and
shear tractions given by
tn = −ρ Cp vn
(2.8)
ts = −ρ Cs vs
where:vn and vs are the normal and shear components of the velocity at the boundary;
ρ is the mass density; and
Cp and Cs are the p- and s-wave velocities.
These viscous terms can be introduced directly into the equations of motion of the gridpoints lying
on the boundary. A different approach, however, was implemented in 3DEC, in which the tractions
tn and ts are calculated and applied at every timestep in the same way as the boundary loads. This
alternative scheme allows the viscous boundaries to be used with rigid blocks as well. Tests have
shown that this implementation is equally effective. The only potential problem concerns numerical
stability, because the viscous forces are calculated from velocities lagging by half a timestep. In
practical analyses to date, no reduction of timestep has been required by the use of the nonreflecting
boundaries. Timestep restrictions demanded by high joint stiffnesses or small zones are usually
more important.
Dynamic analysis starts from some in-situ condition. If a velocity is used to provide the static stress
state, this boundary condition can be replaced by nonreflecting boundaries; the boundary reaction
forces should be maintained throughout the dynamic loading phase. If a stress boundary condition
is applied for the static in-situ solution, a stress boundary condition of opposite sign must also be
applied over the same boundary when the nonreflecting boundary is applied for the dynamic phase.
This will allow the correct reaction forces to be in place at the boundary for the dynamic calculation.
Numerical analysis of the seismic response of surface structures such as dams requires the dis-
cretization of a region of the material adjacent to the foundation. The seismic input is normally
represented by plane waves propagating upward through the underlying material. The boundary
conditions at the sides of the model must account for the free-field motion that would exist in the ab-
sence of the structure. In some cases, elementary lateral boundaries may be sufficient. For example,
if only a shear wave were applied on the horizontal boundary AC, shown in Figure 2.19, it would
be possible to fix the boundary along AB and CD in the vertical direction only. These boundaries
should be placed at sufficient distances to minimize wave reflections and achieve free-field condi-
tions. For soils with high material-damping, this condition can be obtained with a relatively small
distance (Seed et al. 1975). However, when the material damping is low, the required distance may
lead to an impractical model. An alternative procedure is to “enforce” the free-field motion in such
a way that boundaries retain their nonreflecting properties (i.e., outward waves originating from
the structure are properly absorbed). This approach was used in the continuum finite-difference
code NESSI (Cundall et al. 1980). A technique of this type, involving the execution of free-field
calculations in parallel with the main-grid analysis, was developed for 3DEC.
The lateral boundaries of the main grid are coupled to the free-field grid by viscous dashpots to
simulate a quiet boundary (see Figure 2.19); the unbalanced forces from the free-field grid are
applied to the main-grid boundary. Both conditions are expressed in Eqs. (2.9), (2.10) and (2.11),
which apply to the free-field boundary along one side boundary plane with its normal in the direction
of the x-axis. Similar expressions may be written for the other sides and corner boundaries:
In this way, plane waves propagating upward suffer no distortion at the boundary because the free-
field grid supplies conditions that are identical to those in an infinite model. If the main grid is
uniform and there is no surface structure, the lateral dashpots are not exercised because the free-field
grid executes the same motion as the main grid. However, if the main-grid motion differs from that
of the free field (due, for example, to a surface structure that radiates secondary waves), then the
dashpots act to absorb energy in a manner similar to quiet boundaries.
B D
free field
free field
A C
seismic wave
Figure 2.19 Model for seismic analysis of surface structures and free-field
mesh
In order to apply the free-field boundary in 3DEC, the model must be oriented such that the base is
horizontal and its normal is in the direction of the y-axis, and the sides are vertical and their normals
are in the direction of either the x- or z-axis. If the direction of propagation of the incident seismic
waves is not vertical, then the coordinate axes can be rotated such that the y-axis coincides with
the direction of propagation. In this case, gravity will act at an angle to the y-axis, and a horizontal
free surface will be inclined with respect to the model boundaries.
The free-field model consists of four plane free-field grids, on the side boundaries of the model and
four column free-field grids at the corners. The plane grids are generated to match the main-grid
zones on the side boundaries, so that there is a one-to-one correspondence between gridpoints in
the free field and the main grid. The four corner free-field columns act as free-field boundaries
for the plane free-field grids. The plane free-field grids are two-dimensional calculations that
assume infinite extension in the direction normal to the plane. The column free-field grids are
one-dimensional calculations that assume infinite extension in both horizontal directions. Both the
plane and column grids consist of standard 3DEC zones, which have gridpoints constrained in such
a way to achieve the infinite extension assumption.
The model should be in static equilibrium before the free-field boundary is applied. The static
equilibrium conditions prior to the dynamic analysis are transferred to the free field automatically
when the command FF apply is invoked. The free-field condition is applied to lateral boundary
gridpoints. All zone data (including model types and current state variables) in the model zones
adjacent to the free field are copied to the free-field region. Free-field stresses are assigned the
average stress of the neighboring grid zone. The dynamic boundary conditions at the base of the
model should be specified before applying the free-field. These base conditions are automatically
transferred to the free field when the free field is applied. Note that the free field is continuous; if
the main grid contains an interface that extends to a model boundary, the interface will not continue
into the free field.
After the FF apply command is issued, the free-field grid will plot automatically whenever blocks
are plotted. Free-field information can be printed with the LIST ff command.
2.6.2.1 Example Using Dynamic Free Field
A simple model of a concrete gravity dam is created. The top boundary of the model is a free
boundary. The base of the model is a viscous (quiet) boundary. Figure 2.21 shows the geometry of
the model with the dam in the center on top. The model is initially run to static equilibrium under
gravity, to equilibrate the body forces and the boundary forces. This must be done prior to applying
the free field boundaries.
E
pGre0 0 0! u0 i
0..e-
pGrelG,lrG0.3,r3p.0#
30
r
p
The next step is to run the model with only viscous boundaries. An impulse shear stress function
is applied to the base of the model. In this model, the boundaries are too close and it can be seen in
Figure 2.21 that there is amplification of the x-velocities at the base, and distortion of the motion
at the dam crest.
E
pGre0 0 0! u0 i
0,e3,
pGrelG,lrG0.4,r4p.0b
p0$0 0 0
,0$0 0 0
i0
crt
cr
Figure 2.21 x-velocity histories at model base and dam crest using viscous
boundaries
Next, the model is rerun using the free field. The free field command creates new blocks at the
boundary of the model and automatically zones them. Again, the impulse shear stress is applied to
the base. The dynamic input consists of a sinusoidal shear stress wave applied at the model base.
Figure 2.22 shows the x-velocity histories for the base and dam crest. In this case, there is no
amplification of the base wave or distortion at the dam crest.
E
pGre0 0 0! u0 i
0..e-
pGrelG,lrG0.3,r3p.0b
p0$0 0 0
,0$0 0 0
i0
crt
cr
Figure 2.22 x-velocities of the base and dam crest using free field boundaries
damp auto
cycle 2000
;
damp 0 0 mass
reset time hist disp
hist unbal
hist xvel 0 0 0
hist xvel 0 -100 0
hist xvel -40 -28 -20
hist label 2 ’X velocity at crest’
hist label 3 ’X velocity at base’
save ff1.3dsav
res ff1.3dsav
;
In 3DEC, the dynamic input can be applied in one of two ways: (a) as a prescribed velocity history;
or (b) as a stress history. Option (a) enforces an exact given velocity history. If only an acceleration
history is available, it must be integrated numerically to produce a velocity history for 3DEC.
The disadvantage of option (a) is that this boundary will not be an absorbing (or nonreflecting)
boundary (i.e., it will reflect back into the model any outgoing stress waves). To avoid this, option
(b) can be used: the velocity record is transformed into a stress record and applied to a nonreflecting
(viscous) boundary. A velocity history may be converted to a stress boundary condition for similar
purposes using the formula
σn = 2 (ρ Cp ) Vn (2.12)
or
σs = 2(ρ Cs ) Vs (2.13)
1/2
Cp = (K + 4/3G) / ρ
and Cs is given by
Cs = (G / ρ)1/2
The factor of two in Eqs. (2.12) and (2.13) accounts for the fact that the applied stress must be
doubled to overcome the effect of the viscous boundary. Note that, in this case, a velocity history
obtained at the boundary may be different than that from the original velocity record because of the
one-dimensional approximations of Eqs. (2.12) and (2.13).
The process of “baseline correction” is performed on time histories so that the final velocity and/or
displacement is zero. For example, the velocity waveform in Figure 2.23(a) might integrate to give
the displacement waveform in Figure 2.23(b). However, it is possible to determine a low frequency
wave (Figure 2.23(c)) which, when added to the original history, produces a final displacement that
is zero (Figure 2.23(d)). The low frequency wave in Figure 2.23(c) can be a polynomial or a sine
function, with free parameters that are adjusted to give the desired results.
The preceding comments really only apply to complex waveforms derived, for example, from field
measurements. When using a synthetic, simple waveform, it is a simple matter to arrange the
synthetic waveform itself such that the final displacement is zero. Normally, in seismic analysis,
the input wave is an acceleration record. The baseline correction procedure can cause both the final
velocity and displacement to be zero. (For information on standard baseline correction procedures,
consult earthquake engineering texts.)
velocity
time
(a) velocity history
displacement
time
(b) displacement history
velocity
time
(c) low frequency velocity wave
displacement
time
(d) resultant displacement history
An alternative to baseline correction of the input record is to apply a displacement shift at the end of
the calculation if there is a residual displacement of the entire model. This can be done by applying
a fixed velocity to the mesh to reduce the residual displacement to zero. This action will not affect
the mechanics of the deformation of the model. Computer codes to perform baseline corrections
are available from several Internet sites. For example, http://nsmp.wr.usgs.gov/processing.html
provides such a code.
3DEC is primarily intended to simulate structural failure. However, for full practical application,
3DEC should also be able to simulate accurately the dynamic response of structures in the elastic
range. The performance of 3DEC block models subjected to dynamic excitation in the elastic
range can be demonstrated by calculating the natural frequencies and modes of vibration for small
amplitude vibrations.
3DEC contains the command DYNAMIC eigenmodes, which can be executed to calculate natural
frequencies for different modes of vibration in a block system. The command is only available
for rigid block models. The calculation of natural frequencies and modes of vibration assumes
that the structural system is elastic. The kinematic variables of the system of rigid blocks are the
6 degrees of freedom of each block, 3 translations and 3 rotations. In a rigid block model, the
deformability is given by the joint stiffness. A global stiffness matrix for the rigid block system
is assembled, which relates the forces and moments applied to the blocks by their neighbors with
the block displacements and rotations. The mass matrix is assumed to be diagonal, for each block
consisting of 3 entries equal to the block mass and the 3 moments of inertia in the 3 coordinate
directions. This assumption involves an approximation as the moments of inertia in the coordinate
directions are not, in general, the principal moments of inertia of the block.
The global stiffness matrix is formed by assembling the elementary stiffness matrix for each sub-
contact. Assuming small displacements, unit displacements and rotations are considered for each of
the 2 blocks in contact. For each of these 12 configurations, the sub-contact normal and shear forces
are calculated based on normal and shear stiffness and sub-contact area. The forces and moments
that result at the centroid of each of the 2 blocks provide the columns of the contact stiffness matrix
for each of the 12 configurations. Adding the elementary sub-contact matrices leads to the stiffness
matrix of the contact between the 2 blocks. The global stiffness matrix is obtained by assembling
all the contact matrices.
A simple vector iteration procedure is used to calculate the eigenvalues. The algorithm gives the first
N eigenvalues requested, although the ordering may not always be exact. For example, when there
are multiple eigenvalues (e.g., for symmetric structures). It is necessary that at least one block be
fixed and that the system have no completely separate blocks. The dynamic (unscaled) masses must
have been calculated. Therefore, the command SET dynamic on must have been given, followed
by a CYCLE 0 command to force the calculation of dynamic masses. The stiffness matrix requires
a large storage. Therefore, it is recommended that the run not be continued after the eigenvalue
calculation (i.e., save the state before it, and continue the dynamic analysis from that saved file).
The response of a 3DEC rigid block model subjected to low level vibration can be verified for the
analysis of bending of pillars and walls. Lemos (2007) presents two verification examples: elastic
vibration modes of a square pillar and of a wall with variable thickness. A version of the square
pillar example is presented here. The verification involves a pillar, 10 m high, with a square section
of 1 × 1 m, composed of 10 blocks, and assumed clamped at the base (Figure 2.24).
The mass density of the pillar material is 2500 kg/m3. The blocks are rigid, and the normal and
shear stiffnesses are 1.0 GPa/m and 0.4 GPa/m, respectively. The data file for the 3DEC model is
listed in Example 2.10. The command DYNAMIC eigenmodes modes 6 is executed to calculate the
first six modes of vibration. The descriptions of these modes calculated in 3DEC are as follows:
mode 1 – 1st bending mode
mode 2 – mode shape orthogonal to mode 1
mode 3 – 2nd bending mode
mode 4 – mode shape orthogonal to mode 2
mode 5 – 1st torsional mode
mode 6 – 3rd bending mode
Note that in order to provide a more accurate representation of bending behavior, additional contacts
are added between rigid blocks in this model. This is achieved with the FACETRIANGLE rad8
command, which adds a center vertex and 4 mid-edge vertices between two rigid blocks for a total
of 9 point contacts. This improves the calculation for bending moments (see Lemos 2007).
β=0
5
α= 0
4
total
ξ i / ξ min
0
0 5 10 15 20 25 30
ωi
The natural frequencies for the first three bending modes of the pillar are compared to values
calculated from Timoshenko beam theory (Chopra 1995, and Ferreira and Fasshauer 2006). The
natural frequencies for the first three bending modes are listed in Table 2.3. The values compare
reasonably closely to the values from Timoshenko beam theory.
Table 2.3 Natural frequencies (Hz) of bending modes for the square pillar
Bending mode 3DEC rigid block model Timoshenko beam theory
facetriangle rad8
; plot mode 1
dyn eigen setm 1
plot create plot Blocks
plot add block
plot add disp arrowhead on
This problem concerns the effects of a planar discontinuity on the propagation of an incident
shear wave. Two homogeneous, isotropic, semi-infinite elastic regions, separated by a planar
discontinuity with a limited shear strength, are shown in Figure 2.25. A normally incident, plane-
harmonic shear wave will cause slip at the discontinuity, resulting in frictional energy dissipation.
Thus, the energy will be reflected, transmitted and absorbed at the discontinuity. The problem is
modeled with 3DEC, and the results are used to determine the coefficients of transmission, reflection
and absorption. These coefficients are compared with ones given by an analytical solution (Miller
1978).
UT
UI UR
A
The coefficients of reflection (R), transmission (T ) and absorption (A) given by Miller (1978) for
the case of uniform material are
ER
R= (2.14)
EI
ET
T = (2.15)
EI
A= 1 − R2 − T 2 (2.16)
where EI , ET and ER represent the energy flux per unit area per cycle of oscillation associated
with the incident, transmitted and reflected waves, respectively. The coefficient A is a measure of
the energy absorbed at the discontinuity. The energy flux EI is given by
t1 +T
EI = σs vs dt (2.17)
t1
σs = ρ c vs (2.18)
Hence,
t1 +T
EI = ρ c vs2 dt (2.19)
t1
The numerical results are determined for four values of the dimensionless frequency ωγ U/τs ,
where γ = (ρG)1/2 , τs = discontinuity cohesion, U = displacement amplitude of incident wave, ρ
= density of the media and G = shear modulus of the media.
The problem geometry modeled by 3DEC is shown in Figure 2.26. The media were modeled with
elastic, fully deformable blocks of height (h/2), width b and length l. The blocks are separated
by a planar discontinuity extending in the xz-plane. The blocks were internally discretized into
tetrahedral zones, as shown in Figure 2.27. Nonreflecting boundary conditions were used on the top
and bottom of the model. Displacements at boundaries along the yz-plane at x = 0 and x = b were
restrained in the y-direction to simulate plane shear wave conditions. Displacements at boundaries
were restrained in the z-direction along the xy-plane at z = 0 and z = l to simulate the plane-strain
condition.
b = 120 m
z y
h = 400 m
x
Planar
discontinuity
l = 30 m A
Harmonic
shear wave
applied
Figure 2.26 Geometry for the problem of slip induced by harmonic shear
The material properties of the elastic media and planar discontinuity are given in Tables 2.4 and
2.5:
The harmonic shear stress applied at the bottom boundary has the following characteristics:
Note that the magnitude of the incident wave must be doubled in the numerical model to account
for the simultaneous presence of the nonreflecting boundary.
2.9.1.6 Results
The analytical solution for wave propagation through a slipping discontinuity (Miller 1978) assumes
a Mohr-Coulomb discontinuity failure criterion with constant cohesion. In 3DEC, however, when
discontinuity shear and/or tensile strength is exceeded, the cohesion and tension are ignored in all
subsequent calculations. Because the analytical solution assumes a constant discontinuity cohesion
regardless of stress history, a FISH function which prevents setting cohesion and tension to zero when
discontinuity shear and/or tensile strength is exceeded was prepared (see the file “MILR3D.FIS” in
Example 2.12).
An initial run assumed that the discontinuity remains elastic by setting the discontinuity cohesion
to 2.5 MPa. A stress wave of amplitude 1 MPa and frequency 1 Hz was applied at the base of the
model. It should be noted that the displacements and velocities are determined at the nodal points
of the tetrahedron. The stresses, however, are determined at the centroid of the tetrahedron. The
time histories of shear stress, velocity and displacement are monitored at Points A (40, 15, −200)
and B (40, 15, 200). The linear history of stresses are monitored close to Points A and B. The shear
stresses at Points A and B are shown in Figure 2.28. From the amplitude of the stress history at
A and B, it is clear that there was perfect transmission of the wave through the discontinuity. It is
also clear from Figure 2.28 that the nonreflecting boundary condition provides perfect absorption
of normally incident waves.
In the next run, the cohesion was lowered to 0.5 MPa to permit slip along the discontinuity. The
recorded shear stresses at Points A and B are shown in Figure 2.29. The peak stress at Point A
is the superposition of the incident wave and the wave reflected from the slipping discontinuity.
Figures 2.30 and 2.31 show the shear stress at Points A and B for a discontinuity cohesion of 0.1
and 0.02 MPa. It can be seen in Figures 2.29 through 2.31 that the shear stress of Point B is limited
by the discontinuity strength at 0.5, 0.1 and 0.02 MPa, respectively.
E
pGre0 0 0! u0 i
0,GGG
pGrelG,lrG0.:,r:,,0T
r0 0T
p0 0m
i030
m^0
m
Figure 2.28 Shear stress vs time at Points A and B for the case of an elastic
discontinuity (cohesion = 2.5 MPa)
E
pGre0 0 0! u0 i
0,GGG
pGrelG,lrG0.:,r:,,0T
r0 0T
p0 0m
i030
m^0
m
Figure 2.29 Shear stress vs time at Points A and B for the case of an elastic
discontinuity (cohesion = 0.5 MPa)
E
pGre0 0 0! u0 i
0,GGG
pGrelG,lrG0.:,r:,.0T
r0 0T
p0 0m
i030
m^0
m
Figure 2.30 Shear stress vs time at Points A and B for the case of an elastic
discontinuity (cohesion = 0.1 MPa)
E
pGre0 0 0! u0 i
0,GGG
pGrelG,lrG0.:,r:,.0T
r0 0T
p0 0m
i030
m^0
m
Figure 2.31 Shear stress vs time at Points A and B for the case of an elastic
discontinuity (cohesion = 0.02 MPa)
The energy flux, EI , is evaluated using Eq. (2.17) at Point A for non-slipping discontinuities (i.e.,
cohesion = 2.5 MPa). ET is evaluated at Point B for the slipping discontinuity (i.e., cohesion = 0.5,
0.1 and 0.02 MPa). ER is determined at Point B by taking the difference in velocity from the runs with
slipping and non-slipping discontinuities. The coefficients of reflection (R), transmission (T ) and
absorption (A) are computed using Eq. (2.15) (see FISH function energy in file “MILR3D.FIS”
in Example 2.12), and are plotted in Figure 2.31, along with the analytical solution (Miller 1978).
The 3DEC results agree very well with the analytical solution (Figure 2.31).
1.0
0.8
A
R
0.6 T
0.4
0.2
0
0.1 1 2 10 50 100 1000
MC7
Js
Figure 2.32 Comparison of transmission, reflection and absorption coeffi-
cients
;
jset dd=90 dip 0 origin 80,0,0
;
gen edge 31
;
prop mat=1 den=2650 k=16.667e9 g=10.0e9
;
; set boundary material property and viscous boundary along xy
; plane at z=-200 and z=200
bound mat=1 range z -200
bound xvisc zvisc range z -200
bound mat=1 range z 200
bound xvisc zvisc range z 200
;
; set roller boundary along xz plane at y=0 and y=30
bound yvel=0 range y 0
bound yvel=0 range y 30
;
; set zvel=0 along yz plane at x=0 and x=80
bound zvel=0 range x 0
bound zvel=0 range x 80
;
; shear stress along xz plane at z=-200
; set sinusoidal wave function for the
; applied stress at the boundary
; freq = 1 Hz
bound hist sin(1.0,5.0) stress 0,0,0,0,2e6,0 range z -200
;
; set histories
hist n=25
; select zone address and shear stress offset
hist sxz 40,15,-200 sxz 40,15,200
hist xvel(40,15,-200) xvel(40,15,200)
hist xdis(40,15,-200) xdis(40,15,200)
hist time
hist label 1 ’Point A’
hist label 2 ’Point B’
;
insitu stress 0,0,-1e-6,0,0,0
call milr3d.3dfis
save milr3d.3dsav
;
tab 1 0,0
tab 2 0,0
tab 3 0,0
prop jmat=1 jkn=10.0e9 jks=10.0e9 jcoh=2.5e6 jten=1e12
@energy
save milr3d_d.3dsav
ret
;
def _speed
speed_ = sqrt(shear_/dens_)
time1_ = height_/speed_
time2_ = time1_ + 1./freq_
agp_ = gp_near(60,-200,30)
bgp_ = gp_near(60, 200,30)
end
set @shear_ 1.e10 @dens_ 2650. @height_ 400. @freq_ 1.
@_speed
def _store1
while_stepping
if time > time1_ then
if time0_ = 0.0 then
time0_ = time
end_if
if time < time2_ then
rtime_ = time-time0_
table(2,rtime_) = gp_xvel(agp_)
table(3,rtime_) = gp_xvel(bgp_)
end_if
end_if
end
;===================================================
;
; calculate coefficients of transmission, reflection
; and absorbtion
;
;===================================================
def energy
dt_ = tdel
items = table_size(1)
factor = dens_ * speed_
Ei = 0.0
Et = 0.0
Er = 0.0
t_n_1 = 0.0
nac = 0
AAA = 0.0
TTT = 0.0
RRR = 0.0
loop i (1,items)
Vsa_e = ytable(1,i)
Vsa_s = ytable(2,i)
Vsb_s = ytable(3,i)
Ei = Ei + factor * Vsa_e * Vsa_e * dt_
Et = Et + factor * Vsb_s * Vsb_s * dt_
Er = Er + factor * (Vsa_s-Vsa_e) * (Vsa_s-Vsa_e) * dt_
end_loop
RRR = sqrt(Er/Ei)
TTT = sqrt(Et/Ei)
AAA = AAA + sqrt(1.0-RRR*RRR-TTT*TTT)
command
set logfile milr3d.log
set log on
list @RRR
list @AAA
list @TTT
set log off
end_command
end
This problem concerns the dynamic behavior of a single discontinuity under explosive loading.
The problem, shown in Figure 2.33, consists of a planar crack of infinite lateral extent in an elastic
medium, and a dynamic load at some distance, h, from the discontinuity. This problem was
modeled using 3DEC to determine the dynamic response of the discontinuity for a line source (in
the z-direction). The slip of the interface at a Point P , obtained by numerical analysis using 3DEC,
is compared with the closed-form solution derived by Day (1985).
Explosive
Line Source
h P
Crack
Plane
X
ε x
Figure 2.33 Problem geometry for an explosive source near a slip-prone dis-
continuity
The closed-form solution for crack slip as a function of time was derived by Day (1985) and is
given by
where:r = (x 2 + h2 )1/2 , distance from the point source to the point on the crack where the slip
is monitored;
H (τ ) = step function;
τ = t − (r/α);
mo = source strength;
α = velocity of pressure wave;
β = velocity of shear wave;
ρ = density;
1/2
ηα = (α −2 − p2 ) , Re ηα ≥ 0;
1/2
ηβ = (β −2 − p2 ) , Re ηβ ≥ 0;
2
R(p) = (1 − 2β 2 p2 ) + 4β 4 ηα ηβ p2 + 2β ηβ γ ; and
2r 1/2 1/2
p = r2 τ + α x + i τ + α
1 r
τ h .
The slip response of the discontinuity for any source history, S(t), can be obtained by convolution
of Figure 2.33 and the source function, S(t):
Figure 2.34 shows the dimensionless analytical results of slip history at a Point P for a smooth step
function, and for the following values of the variables: α 2 = 3β 2 , h = x and γ = 0. The analytic
solution is implemented in FISH function ana slp (listed in Example 2.16).
0.5
Dimensionless Slip
0.4
0.3
0.2
0.1
DIMENSIONLESS Time
Figure 2.35 shows the problem geometry modeled by 3DEC. The source is located at A
(x = 0 and z = 2 h), and the discontinuity is located at z = h. The dynamic input was ap-
plied at the semicircular boundary of radius 0.05 h. The slip movement was monitored at Point P
on the discontinuity.
The continuous medium was modeled with elastic, fully deformable blocks, as shown in Figure 2.36,
and each block was further discretized into tetrahedral zones. In order to generate only one zone
in the y-direction, the thickness of the block in the y-direction and the average edge length of
the tetrahedron were assumed to be the same. The average edge length was 0.065 h. All of the
joints except the discontinuity were “joined” in order to model a continuous elastic medium. The
discontinuity was assigned a high normal stiffness and high tensile strength in order to meet the
assumption implied in the analytical solution.
Nonreflecting boundary conditions were applied along the horizontal boundaries at the top and
bottom of the model and along the vertical boundary at x = 4 h. A symmetric boundary condition
was applied along the vertical boundary at x = 0. In order to simulate plane-strain conditions,
displacement in the y-direction is restrained along xz-planes at y = 0 and y = 0.065 h.
4h
2h
Dynamic Input P
h Discontinuity
h
h
z y
x Non-Reflecting
Boundary
Figure 2.35 Problem geometry and boundary conditions for numerical model
Figure 2.36 3DEC model showing semicircular source and “joined” blocks
used to provide appropriate zoned discretization
The Mohr-Coulomb joint constitutive relation was used in the analysis. The specific 3DEC param-
eters used for the joint relation are listed in Table 2.7:
Radial velocities corresponding to the dynamic solution for a line source in an infinite medium
were enforced at the semicircular boundary. To avoid problems with the singularity at the source,
dynamic input was applied over a surface 0.05 h from the nominal point source.
UDEC analysis with both velocity and pressure input showed that velocity input gives a better
match with the analytical solution than pressure input. The velocity boundary provides a more
accurate representation of the dynamic stress than the pressure boundary, because in pressure input,
the source function is simply scaled by static stress magnitude and neglects the inertial effects of
dynamic stress at the input boundary.
The radial displacement for a line source given by Lemos (1987) is
−1/2
1 t t 2 α2
u=− −1 , t > r/α (2.21)
2 π α r2 r2
−3/2
1 1 t 2 α2
v=− −1 , t > r/α (2.22)
2 π α r2 r2
The actual input velocity record at r = 0.05 h, as shown in Figure 2.37, was obtained by convoluting
Eqs. (2.22) and (2.20).
E
pGre0 0m 0! u0 i
0.GGr
pGrelG,lrG0.:,p:S.09
p00"
i030
m
Figure 2.37 Input radial velocity time history prescribed at r = 0.05 h (di-
mensionless velocity = (h2 ρβ/mo )v, dimensionless time = τβ/ h)
Velocity history at the boundary at r = 0.05 h is calculated in the FISH function vel inp (listed
in Example 2.15).
2.9.2.6 Results
The dimensionless slip at Point P is plotted against the dimensionless time, and is shown in Fig-
ure 2.38. The dimensionless slip is compared with the analytical solution given by Day (1985).
Velocity input was used on the semicircular region at r = 0.05 h for 3DEC. The error at the peak
slip for 3DEC is 1.7%.
The results shown in Figure 2.38 were obtained with a mesh of maximum zone length of 0.065 h.
The slip response on the discontinuity involves higher frequency components because of zero
friction along the discontinuity. This requires a finer mesh for accurate representation.
The dimensionless slip in Figure 2.38 for 3DEC analysis shows a good agreement with the analytical
solution until the dimensionless time of 1.49. The results show, after dimensionless time of 1.49, a
considerable deviation, which can be attributed to boundary reflections.
Nonreflecting boundaries are used along the top, bottom and right-hand side boundaries. Such
viscous boundaries, designed to absorb normally incident p- and s-waves, cannot be fully effective
in this dynamic slip problem because the discontinuity crosses the boundary. Viscous boundaries,
however, are preferable to roller boundaries.
E
pGre0 09 0! u0 i
0.GGr
pGrelG,lrG0.:,p:S.0R
,0R 0
.0,Tm90%
i030
m myr-y
;
call vel_inp.3dfis
call ana_slp.3dfis
cycle 1
call day3d.3dfis
@_load
;
insitu stress -1.0e-9 -1.0e-9 -1.0e-9 0 0 0
;
; set histories
hist n=10
hist @dtime_ @bvel_ @aslip_ @nslip1_ @nslip2_
hist xvel 0.5 20 0.3 yvel 0 20.5 .03 yvel 0 19.5 0.03
hist time
hist label 2 ’Input Velocity’
hist label 3 ’Analytic Solution’
hist label 4 ’3DEC Results’
;
plot create plot InputVelocity
plot hist 2 vs 9 xaxis label ’Time’ yaxis label ’Input Velocity’
cycle 4000
save day3d.3dsav
plot create plot JointSlip
plot hist 3 4 linestyle style dot vs 9 &
xaxis label ’time’ yaxis label ’Joint Slip’
ret
ny_ = (y_-yc_)/d_
command
bo xvel @nx_ yvel @ny_ hist table @_vtab range id @igp_
end_command
end_if
igp_ = gp_next(igp_)
end_loop
ib_ = b_next(ib_)
end_loop
end
;=====================================================
;
; finds contacts closest to the point of interest
;
;=====================================================
def _find
dfmax_ = 1.e30
dbmax_ = 1.e30
icon_ = contact_head
loop while icon_ # 0
icx_ = c_cx(icon_)
loop while icx_ # 0
x_ = cx_x(icx_)
y_ = cx_y(icx_)
z_ = cx_z(icx_)
dxf_= x_-xf_
dyf_= y_-yf_
dzf_= z_-zf_
dxb_= x_-xb_
dyb_= y_-yb_
dzb_= z_-zb_
df_ = sqrt(dxf_*dxf_+dyf_*dyf_+dzf_*dzf_)
db_ = sqrt(dxb_*dxb_+dyb_*dyb_+dzb_*dzb_)
if df_ < dfmax_ then
dfmax_ = df_
icf_ = icx_
end_if
if db_ < dbmax_ then
dbmax_ = db_
icb_ = icx_
end_if
icx_ = cx_next(icx_)
end_loop
icon_ = c_next(icon_)
end_loop
vscale_ = m0_/(h_n_*h_n_*rho_*v_s_)
dscale_ = 0.25*m0_/(h_n_*rho_*v_s_*v_s_)
end
;=====================================================
;
; stores analytical solution in the histories, and
; convertes numerical solution in the dimensionless
; form
;
;=====================================================
def _compare
while_stepping
dtime_ = time*v_s_/h_n_
bvel_ = table(_vtab,time)/vscale_
aslip_ = table(_utab,dtime_)
xslip1_ = cx_xsdis(icf_)
yslip1_ = cx_ysdis(icf_)
zslip1_ = cx_zsdis(icf_)
xslip2_ = cx_xsdis(icb_)
yslip2_ = cx_ysdis(icb_)
zslip2_ = cx_zsdis(icb_)
nslip1_ = sqrt(xslip1_*xslip1_+yslip1_*yslip1_+zslip1_*zslip1_)
nslip2_ = sqrt(xslip2_*xslip2_+yslip2_*yslip2_+zslip2_*zslip2_)
nslip1_ = nslip1_/dscale_
nslip2_ = nslip2_/dscale_
nslip1_max = max(nslip1_max,nslip1_)
aslip_max = max(aslip_max,aslip_)
end
ret
; Input : Za, Zb
; Output: Z = Re(Z) + Im(Z)
Re_z = Re_a + Re_b
Im_z = Im_a + Im_b
end
;
def mult_complex
; multiplication of two complex variables
; Input : Za, Zb
; Output: Z = Re(Z) + Im(Z)
Re_z = Re_a*Re_b - Im_a*Im_b
Im_z = Re_a*Im_b + Im_a*Re_b
end
;
def divi_complex
; division of complex variables Za/Zb
_deno = Re_b * Re_b + Im_b * Im_b
if _deno = 0.0
divi_compex = 1
exit
endif
Re_z = (Re_a*Re_b + Im_a*Im_b)/_deno
Im_z = (Im_a*Re_b - Re_a*Im_b)/_deno
end
;
def sqrt_complex
; squart root of a complex variable
; Input : Zx
; Output: Zr = Re(Zr) + Im(Zr)
; _theta = atan2(Im_x, Re_x) * 0.5
_arg = float(Im_x/Re_x)
_theta = atan(_arg) * 0.5
_sqrtr = sqrt(sqrt(Re_x*Re_x + Im_x*Im_x))
Re_zr = _sqrtr*cos(_theta)
Im_zr = _sqrtr*sin(_theta)
end
;
def ana_slp
;
; Input _nt, _dt, _xd, _hd, gamma, per, rho
;
_dt = float(_dt)
_xd = float(_xd)
_hd = float(_hd)
gamma = float(gamma)
per = float(per)
rho = float(rho)
;
_vs2 = _vs*_vs
_vs4 = _vs2*_vs2
_2vs2 = 2.0 * _vs2
_4vs4 = 4.0 * _vs4
_r = sqrt(_xd*_xd + _hd*_hd)
_r2 = _r * _r
loop _n (1, _nt)
_t = float(_n) * _dt
_tau = _t - _r/_vp
if _tau > 0.0
_t2r2 = sqrt(_t*_t - (_r/_vp)*(_r/_vp))
Re_cp = _t*_xd / _r2
Im_cp = _t2r2*_hd / _r2
Re_a = Re_cp
Im_a = Im_cp
Re_b = Re_cp
Im_b = Im_cp
mult_complex ; Zˆ 2 ---> Re(Z) + Im(Z)
Re_z2 = Re_z
Im_z2 = Im_z
;
Re_x = 1.0/(_vp*_vp) - Re_z2
Im_x = -1.0 * Im_z2
sqrt_complex ; sqrt(Zx)
Re_cetap = Re_zr
Im_cetap = Im_zr
;
Re_x = 1.0/(_vs*_vs) - Re_z2
Im_x = -1.0 * Im_z2
sqrt_complex ; sqrt(Zx)
Re_cetas = Re_zr
Im_cetas = Im_zr
;
Re_a = 1.0 - _2vs2 * Re_z2
Im_a = -1.0 * _2vs2 * Im_z2
Re_b = Re_a
Im_b = Im_a
mult_complex ; (1. - 2.*vsˆ 2*cpˆ 2) ˆ 2
Re_temp1 = Re_z
Im_temp1 = Im_z
;
Re_a = Re_cetap
Im_a = Im_cetap
Re_b = Re_cetas
Im_b = Im_cetas
mult_complex ; cetap * cetas
Re_a = Re_z
Im_a = Im_z
Re_b = Re_z2
Im_b = Im_z2
mult_complex ; cetap * cetas * cpˆ 2
Re_temp2 = _4vs4 * Re_z
Im_temp2 = _4vs4 * Im_z
;
Re_cr = Re_temp1 + Re_temp2
Im_cr = Im_temp1 + Im_temp2
Re_cr = Re_cr + 2.0 * _vs * gamma * Re_cetas
Im_cr = Im_cr + 2.0 * _vs * gamma * Im_cetas
_dut = 2.0*m0*_vs*_vs/(pi*rho*_vp*_vp)
; note Re_a, Im_a store (cetap*cetas)
Re_b = Re_cp
Im_b = Im_cp
mult_complex ; cetap * cetas * cp
;
Re_a = Re_z
Im_a = Im_z
Re_b = Re_cr
Im_b = Im_cr
if divi_complex = 1 ; cetap * cetas * cp / cr
oo = out(’ divided by zero’)
exit
endif
_dut = _dut * Re_z / _t2r2
ytable(_utab, _n) = _dut
else
ytable(_utab, _n) = 0.0
endif
end_loop
;
_nf = int(per/_dt + 0.0001)
_sum = 0.0
loop _n (1, _nf)
_ph = float(_n) * _dt / per
if _ph < 1.0
ytable(_ftab, _n) = sin(pi * _ph)
else
ytable(_ftab, _n) = 0.0
endif
_sum = _sum + ytable(_ftab, _n)
end_loop
;
; du vs. time relation
;
loop _i (1, _nt)
_uf = 0.0
_n = min(_nf, _i)
loop _j (1, _n)
_uf = _uf + ytable(_utab,_i-_j+1)*ytable(_ftab,_j)
end_loop
ytable(_uftab, _i) = _uf / _sum
xtable(_uftab, _i) = float(_i) * _dt
end_loop
;
; Dimensionless relation
;
loop _n (1, _nt)
ytable(_utab, _n) = (4.0*_hd*rho*_vs*_vs/m0)*ytable(_uftab, _n)
xtable(_utab, _n) = float(_n) * _dt * _vs / _hd
end_loop
;
end
@ana_slp
ret
2.10 References
Bathe, K.-J., and E. L. Wilson. Numerical Methods in Finite Element Analysis. Englewood Cliffs,
New Jersey: Prentice-Hall Inc. (1976).
Belytschko, T. “An Overview of Semidiscretization and Time Integration Procedures,” in Compu-
tational Methods for Transient Analysis, Ch. 1, pp. 1-65. New York: Elsevier Science Publishers,
B.V. (1983).
Biggs, J. M. Introduction to Structural Dynamics. New York: McGraw-Hill (1964).
Chopra, A. K. Dynamics of Structures. Prentice Hall (1995).
Cundall, P. A. “Adaptive Density-Scaling for Time-Explicit Calculations,” in Proceedings of the
4th International Conference on Numerical Methods in Geomechanics (Edmonton, Canada,
1982), pp. 23-26 (1982).
Cundall, P. A., et al. “Computer Modeling of Jointed Rock Masses,” U.S. Army, Engineer Water-
ways Experiment Station, Technical Report WES-TR-N-78-4 (August 1978).
Cundall, P. A., et al. “NESSI – Soil Structure Interaction Program for Dynamic and Static Problems,”
Norwegian Geotechnical Institute, Report 51508-9 (December 1980).
Cundall, P. A., et al. “Solution of Infinite Domain Dynamic Problems by Finite Modelling in
the Time Domain,” in Proceedings of the 2nd International Conference on Applied Numerical
Modelling (Madrid, Spain, September 1978), pp. 341-351. London: Pentech Press (1979).
Day, S. M. “Test Problem for Plane Strain Block Motion Codes,” S-Cubed Memorandum (May 1
1985).
Ferreira, A. J. M., and G. E. Fasshauer. “Computation of natural frequencies of shear deformable
beams and plates by an RBF-pseudospectral method,” Comput. Methods Appl. Mech. Engrg., 196,
134-146 (2006).
Gemant, A., and W. Jackson. “The Measurement of Internal Friction in Some Solid Dielectric
Materials,” The London, Edinburgh, and Dublin Philosophical Magazine & Journal of Science,
XXII, 960-983 (1937).
Gerrard, C. M. “Elastic Models of Rock Masses Having One, Two and Three Sets of Joints,” Int.
J. Rock Mech. Min. Sci. & Geomech. Abstr., 19, 15-23 (1982).
Kuhlmeyer, R. L., and J. Lysmer. “Finite Element Method Accuracy for Wave Propagation Prob-
lems,” J. Soil Mech. & Foundations Div., ASCE, 99(SM5), 421-427 (May, 1973).
Kunar, R. R., P. J. Beresford and P. A. Cundall. “A Tested Soil-Structure Model for Surface
Structures,” in Proceedings of the Symposium on Soil-Structure Interaction (Roorkee University,
India, January 1977), Vol. 1, pp. 137-144. Meerut, India: Sarita Prakashan (1977).
Lemos, J. “A Distinct Element Model for Dynamic Analysis of Jointed Rock with Application to
Dam Foundations and Fault Motion.” Ph.D. Thesis, University of Minnesota (June 1987).
Lemos, J. “Numerical Issues in the Representation of Masonry Structural Dynamics with Discrete
Elements,” in Proceedings of the ECCOMAS Thematic Conference on Computational Methods
in Structural Dynamics and Earthquake Engineering (Rethymno, Crete, Greece, 13-16 June
2007). M. Papadrakakis et al., eds. (2007).
Lysmer, J., and R. Kuhlemeyer. “Finite Dynamic Model for Infinite Media,” J. Eng. Mech., Div.
ASCE, 95:EM4, 859-877 (1969).
Lysmer, J., and G. Waas. “Shear Waves in Plane Infinite Structures,” J. Eng. Mech., Div. ASCE,
98, 85-105 (1972).
Miller, R. K. “The Effects of Boundary Friction on the Propagation of Elastic Waves,” Bull. Seis.
Soc. America, 68, 987-998 (1978).
Myer, L. R., L. J. Pyrak-Nolte and N. G. W. Cook. “Effects of Single Fractures on Seismic Wave
Propagation,” in Proceedings of the International Symposium on Rock Joints, pp. 413-422.
Rotterdam: A. A. Balkema (1990).
Ohnishi, Y., et al. “Verification of Input Parameters for Distinct Element Analysis of Jointed
Rock Mass,” in Proceedings of the International Symposium on Fundamentals of Rock Joints
(Björkliden, Sweden, September 1985), pp. 205-214. O. Stephansson, ed. Luleå: CENTEK
Publishers (1985).
Otter, J. R. H., A. C. Cassell and R. E. Hobbs. “Dynamic Relaxation,” Proc. Inst. Civil Eng., 35,
633-665 (1966).
Roesset, J. M., and M. M. Ettouney. “Transmitting Boundaries: A Comparison,” Int. J. Num. &
Analy. Methods Geomech., 1, 151-176 (1977).
Seed, H. B., P. P. Martin and J. Lysmer. “The Generation and Dissipation of Pore Water Pressures
during Soil Liquefaction,” University of California, Berkeley, Earthquake Engineering Research
Center, NSF Report PB-252 648 (August 1975).
Singh, B. “Continuum Characterization of Jointed Rock Masses: Part I – The Constitutive Equa-
tions,” Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 10, 311-335 (1973).
Wegel, R. L., and H. Walther. “Internal Dissipation in Solids for Small Cyclic Strains,” Physics, 6,
141-157 (1935).
White, W., S. Valliappan and I. K. Lee. “Unified Boundary for Finite Dynamic Models,” J. Eng.
Mech., Div. ASCE, 103, 949-964 (1977).
Wolf, J. P. Dynamic Soil-Structure Interaction. New Jersey: Prentice-Hall (1985).
3.1 Introduction
An important aspect of geomechanical analysis and design is the use of structural support to stabilize
a soil or rock mass. The term “support” describes engineered materials used to restrict displacements
in the immediate vicinity of an opening. Interior support consists of linings, steel sets, etc. which
are placed on the interior of an excavation, and in many cases act to truly support, in whole or part,
the weights of individual blocks isolated by discontinuities or zones of loosened rock. Each support
type (liners or FE blocks) is described separately, below.
Tunnel linings are often thin with respect to tunnel diameters, and their characteristic response to
bending deformation may need to be considered. A structural element representation for the tunnel
lining provides a convenient method for including bending effects.
The structural element method is well-documented in structural engineering texts. The use of beam
elements in 2D linear analysis of excavation support is reported by Dixon (1971), Brierley (1975)
and Monsees (1977), among others. Paul et al. (1983) presents analysis using beam elements
which include nonlinear behavior. Analysis of any support structure is initiated by discretization
of the structure into a number of elements whose response to axial, transverse and flexural loads
can be represented in matrix form, such as that shown in Figure 3.1. The rock-structure interface is
represented by springs oriented both radially and tangentially with respect to the support structure.
,S 2 ,T 2
v 2 u 2
L
m2
2,M2
E ,I ,A
,S 1
v 1
,T 1
u 1
m1
1,M1
T1 A u1
12 I SYM.
S1 0 2
v1
L
6I
M1 0 4I
E L
1
=
T2 L -A 0 0 A u2
12 I 6I 12 I
S2 0 - 2 -
L
0 2
v2
L L
6I 6I
M2 0 - 2I 0 - 4I 2
L L
In general, either an implicit or explicit formulation may be used in analyzing the behavior of
a support structure composed of plate-bending elements and interface stiffnesses. In the first
formulation (implicit), a global stiffness matrix is formed for the entire structure. The size of
the stiffness matrix is reduced by deleting free nodes (i.e., those nodes which are not located at
the rock-support interface). This is possible because these nodes are subjected neither to directly
imposed external loads, nor to displacements by the surrounding medium. The resultant efficiency,
however, limits straightforward application of this formulation to quasi-static problems involving
linear elastic behavior. This formulation does not provide information about failure mechanisms or
ultimate capacities of interior supports. However, factors of safety based on lining stresses should
be conservative since they do not take into account the highly indeterminate nature of a lining in
contact with the rock. A detailed description of a two-dimensional formulation, its use with distinct
elements, and numerous other examples are presented by Lorig (1984).
In the second formulation (explicit), local stiffness matrices are used following division of the
structure into triangular plates with the distributed mass of the structure “lumped” at nodal points,
as shown in Figure 3.2. Forces generated in support elements are applied to the lumped masses
which move in response to unbalanced forces and moments in accordance with the equations of
motion. This formulation has the following desirable characteristics: slip between support and
excavation periphery is modeled in a manner identical to block interaction along a discontinuity;
and large displacements with nonlinear material behavior may be readily accommodated. In the
present formulation, the liner is assumed to behave as a linear-elastic material. These capabilities
are illustrated in 2D in Figure 3.3, where a roof block loads and displaces a hypothetical 4-element
interior structural support.
Excavation
Periphery
m
1. Lumped Mass
2. Structural Element
m
3. Interface
Lining Interior
For three-dimensional analysis, the simplest element (which has the generality of being able to
conform to arbitrary boundaries) is the triangle, interconnected with other elements through lumped
masses located at its vertices. The local stiffness matrix is derived by combining a stiffness matrix
for the in-plane (plane-stress) action and a stiffness matrix for the bending action. This is possible
because the displacements prescribed for the in-plane forces do not affect the bending deformation,
and vice versa, as shown in Figure 3.4. The deformations may be treated independently, provided
the local deformations are small (Zienkiewicz 1977, p. 330).
G (M )
zi zi
u (U )
i i E
x
w (W )
i i
G (M )
xi xi
E
The nodal forces are determined from the displacements in the usual way:
F = Ka (3.3)
where the combined stiffness matrix [K] includes a plane-stress stiffness matrix [KP ] and a bending
stiffness matrix [Kb ] – i.e.,
⎡ ⎤
KP 000 0
⎢ 000 0 ⎥
⎢ ⎥
⎢ 00 0 ⎥
[K] = ⎢
⎢ 00
⎥
⎥ (3.4)
⎢ Kb 0 ⎥
⎣ 00 0 ⎦
00 000 0
Note that rotation θzi does not enter into the definition of deformation as a parameter in either mode.
However, if adjacent elements are coplanar, a difficulty, due to lack of restraint, arises if the zero
stiffness is assigned in the θzi direction. As will be discussed later, a fictitious rotational stiffness,
as well as a fictitious couple, Mz , must be introduced.
In the present formulation, the plane-stress stiffness matrix is taken from Desai and Abel (1972,
p. 132); the bending stiffness matrix is taken from Cheung et al. (1968). These stiffness matrices
are derived using a local coordinate system with the centroid of the triangle defined as the origin. In
the present formulations, all global deformations are transformed into the local coordinate system,
and local “forces” then are transformed into global “forces.”
As mentioned previously, a difficulty arises if elements meeting at a node are coplanar. In the
present formulation, a fictitious set of rotational stiffness coefficients is used with all elements,
whether coplanar or not. For a triangular element, these are defined by the matrix
⎡ ⎤ ⎡ ⎤ ⎡ ⎤
Mzi Kmax −Kmax / 2 θzi
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎢ Mzj ⎥ = ⎢ Kmax −Kmax / 2 ⎥ ⎢ θzj ⎥ (3.5)
⎣ ⎦ ⎣ ⎦ ⎣ ⎦
Mzk Sym Kmax θzk
where Kmax is equal to the largest sum of terms in any row of the [Kb ] element matrix.
Element stresses are determined based on total local deformations at the nodes. For the in-plane
stresses, the stresses {σ } are given by
Also, in the present formulation, strains are assumed to be constant within an element, and therefore,
in-plane stresses are also assumed constant.
Bending “stresses,” {σ }T , can also be determined from the “stress” matrix relation for plate-bending
elements in a similar fashion to Eq. (3.6) with
{σ }T = {Mx My Mz } (3.8)
and
The “stress” resultants for the plate-bending element are shown in Figure 3.1. It is assumed that
true stresses vary linearly across the plate thickness – e.g.,
12 Mx
σx = z (3.10)
t3
where:z is measured from the plate mid-plane; and
t is the plate thickness.
The triangular plate bending given by Cheung et al. (1968) tries to represent the average moment
values over its area, rather than follow them by a linear variation (although some terms of such vari-
ation are apparently included). Therefore, in the present formulation, moment values are computed
only at the centroid of each triangular element.
The structural liner elements used in 3DEC require the following input parameters:
(1) Young’s modulus of elasticity of liner material [force/area]
(2) cohesion for contact between liner and host medium [force/area]
(3) friction for contact between liner and host medium [degrees]
(4) normal stiffness for contact between liner and host medium [force/area/disp]
(5) shear stiffness for contact between liner and host medium [force/area/disp]
(6) Poisson’s ratio of liner material
(7) thickness of liner material [length]
(8) tensile limit for contact between liner and host medium [force/area]
All of the commands associated with structural element liners are listed in Table 3.1. See Sec-
tion 1.2.15 in the Command Reference for a detailed explanation of the commands.
This is an example showing the application of a structural liner in the inside of a square tunnel. The
geometry of the tunnel is shown in Figure 3.5. The block in the roof of the tunnel will fall if the
tensile strength in the horizontal joint above it is not sufficient to support its weight. The model is
used to simulate a 20 cm-thick concrete liner. The properties of the liner are
Young’s modulus 15 GPa
Poisson’s ratio .15
thickness .2 m
friction in contact between liner and rock 60 degrees
cohesion in contact between liner and rock .5 MPa
tensile limit of contact between liner and rock .3 MPa
normal stiffness of contact 1.0 GN/m
shear stiffness of contact 1.0 GN/m
The input commands for this example are listed in Example 3.1. Figure 3.6 shows the vertical
displacement history of the roof block with no liner in place. The block is free to fall in this case.
Figure 3.7 shows the vertical displacement history of the same block with the concrete liner in
place.
E
pGre00 0 u0 i
0rpGG
pGrelG,lrG0.:,,:Gr0
:0
r
E
pGre009 0 u0 i
0,pGG
pGrelG,lrG0.:,,:Gr0x
r0l 00,3.4p
i0
e e nm lH
e nm H
E
pGre009 0 u0 i
0,pGG
pGrelG,lrG0.:,,:Gr0x
r0l 00,3.4p
i0
e nm 5
(1) While the interface between the rock surface and the liner can slip or detach,
the liner element itself is linearly elastic and cannot yield or rupture.
(2) The liner placement logic cannot be used to line tunnel intersections.
(3) Block zone discretization must be sufficiently small, so that each node of the
structural liner will fall in a different zone.
(4) Liners have no contact detection logic. Blocks which have no liner nodes
attached may pass through the liner.
Property numbers are assigned to interior support elements with the PROPERTY command. All
quantities must be given in an equivalent set of units (see Table 3.2). The code does not take into
account the weight of the structure when calculating loads. The mass density for structural (beam)
elements may be artificially increased above the actual values to achieve reasonable timesteps for
quasi-static problems. A convenient way to determine how much the mass density may be increased
is to first determine the minimum timestep without any structural elements. Next, structural elements
are introduced and the mass density gradually increased until the timestep approaches the previous
timestep (i.e., timestep without structural elements).
This section describes a finite element formulation which can be used to model tunnel liners and
other structures attached to 3DEC blocks (such as dams). Its implementation within the explicit
solution scheme used by 3DEC, and the logic that was developed to handle the contact between
the special FE blocks and the standard 3DEC blocks, are also presented. The input commands and
keywords added to 3DEC in order to read or generate the FE block mesh and assign material models
and properties are explained.
The new module creates special blocks (FE blocks) that can be added to a 3DEC model composed
of standard polyhedral deformable blocks. The new FE blocks consist of a finite element mesh
composed of three-dimensional brick elements with 20 nodes. The element formulation, presented
in this section, follows the standard methodology found in finite element texts (Hughes 1987).
The geometry of isoparametric finite elements is described by the mapping of a master cubic element
into the actual element shape, as shown in Figure 3.5, which also indicates the node numbering
convention adopted. For the case of the 20-node brick element, the edges are defined by parabolic
segments.
3.3.1.1 Notation Conventions
Uppercase subscripts, such as I , denote element nodes ranging from 1 to 20. Lowercase subscripts,
such as j , k and m, denote coordinate directions and range from 1 to 3. Summation on repeated
indices is implied.
3.3.1.2 Geometry Mapping
The geometry of the master element (a cube with 8 nodes placed at the corners and 12 nodes at
the edge midpoints) is defined in the y-coordinate system, with values ranging from −1 to 1, as
indicated in Figure 3.5.
The coordinates of a generic point of the master element are denoted by
y = {y1 y2 y3 }T (3.11)
x = {x1 x2 x3 }T (3.12)
xI = {xI 1 xI 2 xI 3 }T (3.13)
xk = NI (y) xI k (3.14)
which expresses the coordinates, x, of a generic point as a sum of products of the nodal shape
functions, NI , function of the natural coordinates, y, and the nodal coordinates, xI .
The Jacobian matrix of the transformation is given by
∂ym
[Ymk ] = = ([Jkm ])−1 (3.16)
∂xk
The expressions of the shape functions, NI , and their derivatives with respect to the master element
coordinates, ∂/∂ym (NI ), are given in Section 3.3.7.
3.3.1.3 Displacements and Strains
The displacement field of isoparametric elements is governed by an expression similar to the ge-
ometry mapping (Eq. (3.11)):
uk = NI (y) uI k (3.17)
u = {u1 u2 u3 }T (3.18)
uI = {uI 1 uI 2 uI 3 }T (3.19)
1 ∂uk ∂um
εkm = ( + ) (3.20)
2 ∂xm ∂xk
the derivatives of the displacement field must be calculated from Eq. (3.17),
∂uk ∂NI
= uI k (3.21)
∂xm ∂xm
The derivatives of the shape function in the global coordinates are obtained by application of the
chain rule:
The components of the strain tensor are grouped more conveniently in vector form:
Then, after introducing Eqs. (3.22) and (3.21) in Eq. (3.20), and setting the result in matrix form,
the strains may be expressed as
ε = BI uI (3.24)
⎡ ⎤
NI 1 0 0
⎢ 0 NI 2 0 ⎥
⎢ ⎥
[BI ] = ⎢ 0 0 NI 3 ⎥
⎢ ⎥ (3.25)
⎢ 0 NI 3 NI 2 ⎥
⎣ ⎦
NI 3 0 NI 1
NI 2 NI 1 0
As the force-displacement relations for the element cannot be expressed in closed form, numerical
integration has to be used. The element strains and stresses are evaluated at the Gauss integration
points of the element. Once the strains (and strain increments) are obtained from Eq. (3.24), the
application of the assumed constitutive model produces the new stresses. Formally,
σ = f (ε) (3.26)
FI S = {FI 1 FI 2 FI 3 }T (3.28)
and is calculated as
FI =
S
BI T σ dV (3.29)
v
The integration, extending to the element volume, is performed numerically as a summation of the
values of the integrand at the Gauss points – i.e.,
FI S = BI T σ WG VG (3.30)
G
where VG are the elementary volume terms obtained from the Jacobian matrix of the geometry
transformation, and WG are the standard Gauss point weights (Hughes 1987). In the present
implementation, 2 × 2 × 2 or 3 × 3 × 3 Gauss points may be used. The former corresponds to
reduced integration, and, in some cases of unconstrained systems, may allow spurious mechanisms.
However, since it is faster, and the stresses at the location of the 2 × 2 × 2 Gauss points tend to be
more accurate, it is taken as the default option.
g = {g1 g2 g3 }T (3.31)
are given by
FI G
= NI g dV (3.32)
V
Concentrated forces, such as point contact forces applied on the element faces, originate nodal
forces given by, for node I ,
FI P = NI P (3.33)
where the shape function is evaluated at the point of application of the force
P = {P1 P2 P3 }T (3.34)
The implementation of the 20-node brick element was designed to allow a smooth integration of the
special blocks into the 3DEC logic. In particular, the new blocks were intended to be compatible
with the essential aspects of the contact logic, the solution algorithm and the graphical routines, as
described in the following sections.
The standard deformable blocks in 3DEC have a polyhedral form, with an internal tetrahedral
mesh, leading to an outer surface defined by triangular faces. The contact logic and the graphical
representation are based on triangulated surfaces. In general, the surface of the 20-node brick is
curved, which would require a different, and computationally much more expensive, treatment of
the geometric contact calculations. Therefore, in order to allow the special FE blocks to be handled
by the same routines, it was decided to approximate the block external boundary by means of a
polyhedral surface composed of triangles. As shown in Figure 3.6, eight triangular faces are used
to approximate each face of the brick.
The nodes of the element, either master or slave, perform the same functions as zone gridpoints, and
are similarly treated as vertices for contact purposes. Therefore, the standard 3DEC contact logic
is used, based on sub-contacts located at vertex-to-face and edge-to-edge interaction locations.
The FE block displacements at the sub-contacts are evaluated assuming the approximated triangu-
lated surface, and are therefore obtained from the gridpoint displacements in the same way as for
the polyhedral blocks. The contact forces that are applied to the slave node at the center of the brick
faces must be transferred to the master nodes using Eq. (3.33).
The common-plane logic that is used in 3DEC in the contact detection routines assumes that blocks
are convex. If FE blocks are non-convex, then no contact is detected in the concave portion of their
surface. If a concave face must be in contact with other blocks, then it is necessary to divide the
FE block into convex blocks, which may be joined to behave as a monolithic block, as is done with
regular blocks.
External loads applied with the BOUNDARY command are also applied to FE blocks, based on the
triangulated boundary. The loads applied at the slave nodes are automatically transferred to the
master nodes using Eq. (3.33), as in the case of the contact forces. The only difference with respect
to the regular blocks is in the treatment of gravity loads which must be introduced by means of
Eq. (3.22).
Fixed displacement conditions are applied in 3DEC by prescribing gridpoint velocities using the
BOUNDARY command. This procedure is also used for FE nodes.
The solution procedure in 3DEC is based on an explicit timestepping algorithm that integrates the
nodal equations of motion. The same procedure is used for dynamic and quasi-static analysis. In
the latter case, high damping values are introduced to obtain convergence to static equilibrium or
to a failure mechanism. In the case of deformable blocks, discretized into zones or elements, the
equations of motion of each gridpoint are integrated. The FE blocks are treated in the same way –
i.e., for a given node I , the equations of motion may be expressed as
where m is the nodal mass, α is the viscous damping parameter and FI is the nodal force vector.
The nodal forces are calculated as a sum of four components:
FI = FI S + FI C + FI B + FI G (3.36)
where:FI S are the forces equivalent to the element stresses, calculated by Eq. (3.30);
FI C are obtained from the contact forces, using Eq. (3.33);
FI B are obtained from the applied boundary forces, using Eq. (3.33);
FI G are the gravity forces, given by Eq. (3.32).
The calculation of the timestep required for numerical stability is based on the fictitious discretization
of the element into zones. Since the assemblage of these zones is known to behave in a stiffer manner
than the actual 20-node brick, this approximation provides a safe estimate. For static analysis, the
mass scaling procedure gives the gridpoint masses. For dynamic analysis, the tetrahedral mesh is
still used to give the gridpoint masses, in such a way that the total element mass is preserved. More
accurate schemes exist to provide a diagonal mass matrix for these types of isoparametric elements
(Hughes 1987), as explicit solution methods do not permit the use of the consistent mass matrices
derived from the standard element formulation.
In order to employ FE blocks in a 3DEC model, a configuration option must be selected with the
command CONFIG feblock, so that the needed extra storage is allocated. The FEBLOCK command is
used for most actions involving the special blocks which are invoked by specific keywords. There
are two ways of introducing the special blocks into the model:
(1) reading a previously generated FE mesh from a file; and
(2) automatically generating a mesh inside a brick-shaped block.
In the first case, invoked by the keyword read of the command FEBLOCK, an ASCII file must be
created with an FE mesh described in the typical format of FE codes:
(1) title;
(2) header line with number of nodes (NN) and number of elements (NE);
(3) nodal coordinates – NN lines with x y z coordinates for each node; and
(4) element definition – NE lines with list of nodes for each element.
Further details are presented in Section 3.3.8. Node numbering is irrelevant, given 3DEC ’s solution
procedure. It is possible to supply the full list of 20 nodes for each element, or to read in a mesh of
8-node brick elements (i.e., only with the corner nodes). The code then automatically inserts the
extra mid-edge nodes and transforms these elements into 20-node elements.
The automatic generation of the FE mesh is only possible inside 8-vertex brick blocks. These may
be created with any of the 3DEC block generation and cutting commands, and then the command
FEBLOCK gen transforms them into FE blocks by creating a regular mesh of 20-node elements.
The user only indicates the desired dimensions of the elements in each of the three block axes.
The mechanical behavior of FE blocks may be governed by any of the regular constitutive models
available in 3DEC. The material and constitutive number are assigned with FEBLOCK change.
Different materials may be applied to the various elements belonging to the same block.
For application to arch dams, special routines were developed to allow the adequate geometric
fitting of the concrete structure (represented by one or more FE blocks) and the rock mass in the
foundation (represented by regular polyhedral blocks). Typically, the FE mesh of the concrete shell
will be generated by some preprocessing software for dam engineering applications, and read into
3DEC from a file. Then the new routines are invoked to create a set of regular blocks that fit the
foundation surface. All of these blocks are joined, and afterwards, the rock mass discontinuities
are introduced with the cutting commands.
The first step is to define the foundation surface of the FE mesh. Surface region numbers, assigned to
block faces with FEBLOCK mark sregion, are useful for this purpose. As the faces of the foundation
surface may be curved in general, the keyword linearize is available to set the nodes along the edges
of these faces exactly at the midpoint of the edges, so that they can fit the polyhedral foundation
blocks. Then the keyword base, with suitable geometric arguments, is invoked to create blocks
directly under the foundation faces, as well as blocks extending upstream and downstream.
The isoparametric 20-node brick element is derived from a cubic master element defined in the
y-coordinate system (Figure 3.5). The positions of the master element nodes are listed in Table 3.3.
Nodes 1 to 8 are placed at the corners, and nodes 9 to 20 at the midpoint of the edges.
The coordinates of the slave nodes used in the face triangulation, and numbered from 21 to 27, are
listed in Table 3.4. Nodes 21 to 26 are located at the mid-face points; node 27 is located at the
center of the element.
The definition of the nodes of the 6 element faces is given in Table 3.5. The first 4 nodes correspond
to the corners, the next 4 to the mid-edges, and the last entry to the slave node at mid-face position.
The shape functions for the 20-node brick element are based on 2nd-degree polynomials in the
y-reference system. The expressions and their derivatives with respect to the y-coordinates are
listed below. The constants yI k stand for the coordinates of node I , as given in Table 3.3.
Corner nodes: I = 1 to 8
NI = 18 (1 + yI 1 y1 )(1 + yI 2 y2 )(1 + yI 3 y3 )(yI 1 y1 + yI 2 y2 + yI 3 y3 − 2)
∂NI
∂y1 = 18 yI 1 (1 + yI 2 y2 )(1 + yI 3 y3 )(−2yI 1 y1 + yI 2 y2 + yI 3 y3 − 1)
∂NI
∂y2 = 18 yI 2 (1 + yI 3 y3 )(1 + yI 1 y1 )(−2yI 2 y2 + yI 1 y1 + yI 3 y3 − 1)
∂NI
∂y3 = 18 yI 3 (1 + yI 1 y1 )(1 + yI 2 y2 )(−2yI 3 y3 + yI 1 y1 + yI 2 y2 − 1)
zstress assigns the current Gauss point stresses to the nearest zone of the
parallel discretization, for plotting purposes. 3DEC cannot directly
access the stresses in the FE blocks for plotting purposes. However,
3DEC maintains a fictitious set of zones which correspond in loca-
tion to the FE zones. The zstress keyword transfers the Gauss point
stresses into the fictitious zones so they may be plotted. (This is done
automatically at the end of the CYCLE command.)
In this example we will demonstrate 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 data file used to create this example is listed in Example 3.2. 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 POLYHEDRON prism command. The FISH function also sets up FISH
parameters that are used for dimensions in later commands. This geometry is shown in Figure 3.8:
E
pGre00 0 u0 i
0SSGGG
pGrelG,lrr0.4r,4S.0
40
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 commands
feb gen 20 20 20 range reg 2 3 4
feb gen 20 20 13 range reg 1 5
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 surface region numbers on the dam blocks. This is done using a
FISH function (markdamfaces) to set the sregion (face sreg) for the block faces. The block
face surface regions are then used by the FEBLOCK command to extrude the existing blocks down
to the base of the model. The extruded blocks for this step are shown in Figure 3.9.
feblock base sreg 1 proj z @zzbotbou
feblock base sreg 2 proj z @zzbotbou
feblock base sreg 3 proj z @zzbotbou
feblock base sreg 4 proj z @zzbotbou
feblock base sreg 5 proj z @zzbotbou
E
pGre00 0 u0 i
0G
pGrelG,lrr0.:r,:S.0
:0
The base blocks are then extruded in the y-direction to form the upstream and downstream base
blocks as shown in Figure 3.10.
feblock base sreg 3 axis -1 proj y @yydbou proj z @zzbotbou
feblock base sreg 3 axis 1 proj y @yyubou proj z @zzbotbou reg 21
feblock base sreg 2 axis -3 proj y @yydbou proj z @zzbotbou
feblock base sreg 2 axis 3 proj y @yyubou proj z @zzbotbou reg 21
feblock base sreg 4 axis -3 proj y @yydbou proj z @zzbotbou
feblock base sreg 4 axis 3 proj y @yyubou proj z @zzbotbou reg 21
E
pGre00 0 u0 i
0G
pGrelG,lrr0.:r,:SS0
:0
Figure 3.10 Base blocks extruded in the upstream and downstream directions
The upstream and downstream abutments are then created using POLYHEDRON brick commands.
The result of these commands is shown in Figure 3.11.
; abutments downstream
poly prism a @xxright @yyright @zzright @xxrightd @yyrightd @zd &
@xxrightd @yyrightd @zztopbou @xxright @yyright @zztopbou &
b @xxright @yydbou @zzright @xxrightd @yydbou @zd &
@xxrightd @yydbou @zztopbou @xxright @yydbou @zztopbou &
reg 0
poly prism a @xxleft @yyleft @zzleft @xxleftd @yyleftd @zd &
@xxleftd @yyleftd @zztopbou @xxleft @yyleft @zztopbou &
b @xxleft @yydbou @zzleft @xxleftd @yydbou @zd &
@xxleftd @yydbou @zztopbou @xxleft @yydbou @zztopbou &
reg 0
E
pGre00 0 u0 i
0G
pGrelG,lrr0.:r,:SS0
:0
A fault is cut into the right abutment, and the blocks are zones as shown in Figure 3.12.
find
hide reg 21
hide reg 1 2 3 4 5
jset dip 90 dd -80 or @xxleft @yyleft 80
jset dip 90 dd -65 or 70 -60 80
join reg 2
join reg 3
join reg 4
E
pGre00 0 u0 i
0SSGGG
pGrelG,lrr0.4r,4SS0
40
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.
;
; rock : mat=1
change mat 1
; dam : mat=2
change mat 2 range region 1 2 3 4 5
; rock
; E=20000 MPa, v=0.2, K=11111, G=8333
prop mat 1 dens 0.0027 k 11111 g 8333
; concrete
; use E/100 for insitu stage
; E=200 MPa, v=0.2, K=111.11, G=83.33
prop mat 2 dens 0.0024 k 111.11 g 83.33
; rock joints
prop jmat 1 kn 10000 ks 5000 jfric 38
prop jmat 6 kn 10000 ks 5000 jfric 38
prop jmat 7 kn 10000 ks 5000 jfric 38
Gravity is defined, and in-situ stresses are applied to the rock and the dam. Then the dam is hidden,
and the INSITU topographic command is used to initialize the valley stresses. The stresses are then
removed from the dam.
; --- stage 1 : insitu stresses ---
; (dam blocks soft and elastic)
; gravity
grav 0 0 -10
insitu stress -1.0e-4 -1.0e-4 -1.0e-4 0 0 0
; hide dam
hide mat 2
; topographic stresses
ins topo zup kox 0.3 koy 0.3
; find dam
find
save damfebs0.sav
; remove dam stresses
reset stress range reg 1 2 3 4 5
reset jstr range mint 1 2
reset jstr range mint 2 2
At this point we want to run the model to equilibrium. We can use the FEB 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.
; remove dam weight
feb grav del
; boundary
bound yv 0.0 range y -180
bound yv 0.0 range y 150
bound xv 0 range x -150
bound xv 0 range x 150
bound zv 0 range z -100
; hist
hist unbal
hist xdis 0 0 80
hist ydis 0 0 80
hist zdis 0 0 80
hide mat 2
hist xdis 0 0 0
hist ydis 0 0 0
hist zdis 0 0 0
find
set dyn off
set damp local
cy 15000
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.
; --- stage 2 : dam self weight ---
; (real dam properties, elastic joints)
rest damfebs1.sav
reset disp
reset jdisp
reset damp
; remove dam stresses
reset stress range reg 1 2 3 4 5
reset jstr range mint 2 2
reset jstr range mint 1 2
; concrete properties (real E)
; E=20000 MPa, v=0.2, K=11111, G=8333
prop mat 2 k 11111 g 8333
; hist
reset hist
hist unbal
hist xdis 0 0 80
hist ydis 0 0 80
hist zdis 0 0 80
hide mat 2
hist xdis 0 0 0
hist ydis 0 0 0
hist zdis 0 0 0
find
hist xdis 70 -40 66
hist ydis 70 -40 66
hist zdis 70 -40 66
hist xdis 57 -51 40
hist ydis 57 -51 40
hist zdis 57 -51 40
In this stage we use a boundary condition to add water pressures to the dam:
; --- stage 3 - water load ---
; (elastic joints)
rest damfebs2.sav
reset damp
; boundary, fix base
bound xv 0 yv 0 zv 0 range z -100
; water load
bou str -0.8 -0.8 -0.8 0 0 0 zgrad 0.01 0.01 0.01 0 0 0 range sreg @isregup
pri bou for range sreg 7
step 5000
And we use a FISH function (ppjoint)to add pore pressures to the joints:
; --- stage 5 : pp in rock joints ---
@ppjoint
set nucx 0
step 10000
And 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.
step 5000
step 5000
step 5000
save damfebs6.sav
The final result in which the fault has slipped is shown in Figure 3.13:
E
pGre00m 0 u0 i
0SSGGG
pGrelG,lrr0.4r,4Se0
m 40x
o
40eip..5r
40S
rad = 100.0
xcent = 0.0
ycent = -rad
thick = 10.0
angtot = 100.0 * degrad
nb = 5
angb = angtot / nb
zz2 = 80.0
zz1 = 0.0
loop i (1,nb)
command
poly prism a @xx1 @yy1 @zz1 @xx2 @yy2 @zz1 @xx3 @yy3 @zz1 &
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
hide reg 1 5
jset p3 @xa1 @ya1 @za @xb1 @yb1 @za @xc1 @yc1 @za
join reg 3
;
; right abutment
hide reg 3 4
del zr @zz1 @za
poly 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 &
reg 2
join reg 2
find
;
hide reg 2 3 4 5
jset p3 @xa1 @ya1 @za @xb1 @yb1 @za @xd1 @yd1 @zd
del zr @zz1 @zd07
find
;
; left abutment
hide reg 1 2 3 5
del zr @zz1 @za
poly 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 &
reg 4
join reg 4
find
;
hide reg 1 2 3 4
jset p3 @xa2 @ya2 @za @xb2 @yb2 @za @xd2 @yd2 @zd
del zr @zz1 @zd07
find
endcommand
xx45 = xa2
yy45 = ya2
zz45 = za
xx45d = xb2
yy45d = yb2
zz45d = za
izzright = gp_near(xxright,yyright,64.0)
izzleft = gp_near(xxleft,yyleft,64.0)
zzright = gp_z(izzright)
zzleft = gp_z(izzleft)
ii=out(’ zzright ’+string(zzright))
ii=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
; foundation joint
ii=out(’ found. joint’)
ib = block_head
loop while ib#0
ii=out(’ ib ’+string(ib))
section
if b_region(ib) # 1
if b_region(ib) # 5
if b_z(ib) > za
exit section
endif
endif
endif
; fe faces
ibifef = ib+$kbfef
ii=out(’ ibifef ’+string(ibifef))
ifef = imem(ib+$kbfef)
ii=out(’ ifef ’+string(ifef))
loop while ifef#0
ii=out(’ ifef ’+string(ifef))
ifefvn = index(ifef+$kfefvn)
vnx = fmem(ifefvn)
vny = fmem(ifefvn+1)
vnz = fmem(ifefvn+2)
ii=out(’ vn ’+string(vnx)+’ ’+string(vny)+’ ’+string(vnz))
if vnz < -0.1
nnfefj = nnfefj+1
imem(index(ifef+$kfefreg)) = b_region(ib)
ii=out(’ marked ifef ’+string(ifef)+’ b ’+string(b_region(ib)))
ii=out(’ vn ’+string(vnx)+’ ’+string(vny)+’ ’+string(vnz))
endif
ifefn = imem(index(ifef+$kfefn))
ifef = ifefn
endloop
; zone faces
ifa = b_face(ib)
loop while ifa#0
ii=out(’ ifa ’+string(ifa))
vn = face_n(ifa)
vnx = xcomp(vn)
vny = ycomp(vn)
vnz = zcomp(vn)
ii=out(’ vn ’+string(vnx)+’ ’+string(vny)+’ ’+string(vnz))
if vnz < -0.1
nnzj = nnzj+1
face_sreg(ifa) = b_region(ib)
ii=out(’ marked ifa ’+string(ifa)+’ b ’+string(b_region(ib)))
ii=out(’ vn ’+string(vnx)+’ ’+string(vny)+’ ’+string(vnz))
endif
ifa = face_next(ifa)
endloop
endsection
ib = b_next(ib)
endloop
; upstream face
ii=out(’ upstream face’)
ib = block_head
loop while ib#0
ii=out(’ ib ’+string(ib))
section
; fe faces
ibifef = ib+$kbfef
ii=out(’ ibifef ’+string(ibifef))
ifef = imem(ib+$kbfef)
ii=out(’ ifef ’+string(ifef))
loop while ifef#0
ii=out(’ ifef ’+string(ifef))
ifefvn = index(ifef+$kfefvn)
vnx = fmem(ifefvn)
vny = fmem(ifefvn+1)
vnz = fmem(ifefvn+2)
end
@markdamfaces
save dampanels
save damextrude
; abutments downstream
poly prism a @xxright @yyright @zzright @xxrightd @yyrightd @zd &
@xxrightd @yyrightd @zztopbou @xxright @yyright @zztopbou &
b @xxright @yydbou @zzright @xxrightd @yydbou @zd &
@xxrightd @yydbou @zztopbou @xxright @yydbou @zztopbou &
reg 0
poly prism a @xxleft @yyleft @zzleft @xxleftd @yyleftd @zd &
@xxleftd @yyleftd @zztopbou @xxleft @yyleft @zztopbou &
b @xxleft @yydbou @zzleft @xxleftd @yydbou @zd &
save damabut
join reg 2
join reg 3
join reg 4
save damfebz.3dsav
; rock : mat=1
change mat 1
; dam : mat=2
change mat 2 range region 1 2 3 4 5
; rock
; E=20000 MPa, v=0.2, K=11111, G=8333
prop mat 1 dens 0.0027 k 11111 g 8333
; concrete
; use E/100 for insitu stage
; E=200 MPa, v=0.2, K=111.11, G=83.33
prop mat 2 dens 0.0024 k 111.11 g 83.33
; rock joints
prop jmat 1 kn 10000 ks 5000 jfric 38
prop jmat 6 kn 10000 ks 5000 jfric 38
prop jmat 7 kn 10000 ks 5000 jfric 38
; gravity
grav 0 0 -10
; hide dam
hide mat 2
; topographic stresses
ins topo zup kox 0.3 koy 0.3
; find dam
find
save damfebs0.3dsav
; boundary
bound yr -180.1 -179.9 yv 0.0
bound yr 149.9 150.1 yv 0.0
bound xr -150.1 -149.9 xv 0
bound xr 149.9 150.1 xv 0
bound zr -100.1 -99.9 zv 0
; hist
hist unbal
hist xdis 0 0 80
hist ydis 0 0 80
hist zdis 0 0 80
hide mat 2
hist xdis 0 0 0
hist ydis 0 0 0
hist zdis 0 0 0
find
cy 15000
save damfebs1.3dsav
pri max
pri co sum
reset disp
reset jdisp
reset damp
; hist
reset hist
hist unbal
hist xdis 0 0 80
hist ydis 0 0 80
hist zdis 0 0 80
hide mat 2
hist xdis 0 0 0
hist ydis 0 0 0
hist zdis 0 0 0
find
hist xdis 70 -40 66
hist ydis 70 -40 66
hist zdis 70 -40 66
hist xdis 57 -51 40
hist ydis 57 -51 40
hist zdis 57 -51 40
save damfebs2.3dsav
pri max
pri co sum
rest damfebs2.3dsav
reset damp
; set wr on
; water load
bou str -0.8 -0.8 -0.8 0 0 0 zgrad 0.01 0.01 0.01 0 0 0 range sreg @isregup
step 5000
save damfebs3.3dsav
pri max
pri co sum
reset damp
step 5000
save damfebs4.3dsav
def ppjoint
ic = contact_head
loop while ic#0
icx = c_cx(ic)
jmat = c_mat(ic)
if jmat = 0
jmat = cx_mat(icx)
endif
ib1 = c_b1(ic)
ib2 = c_b2(ic)
ir1 = b_region(ib1)
ir2 = b_region(ib2)
ibm1 = b_mat(ib1)
ibm2 = b_mat(ib2)
vnx = c_nx(ic)
vny = c_ny(ic)
vnz = c_nz(ic)
section
if icx = 0
exit section
endif
xx = cx_x(icx)
yy = cx_y(icx)
zz = cx_z(icx)
ar = cx_area(icx)
iver = cx_vertex(icx)
ppold = cx_pp(icx)
flag = 1
if flag=1
ppnew = (80-zz)*0.01
endif
icx = cx_next(icx)
endloop
endsection
ic = c_next(ic)
endloop
end
@ppjoint
;
set nucx 0
step 5000
step 5000
save damfebs5.3dsav
; xxxxxxxxxxxxxxxxxxxxxxx
prop jmat 1 jfric 20
step 5000
step 5000
step 5000
save damfebs6.3dsav
ret
3.4 References
Brierley, G. “The Performance during Construction of the Liner of a Large, Shallow Underground
Opening in Rock.” Ph.D. Thesis, University of Illinois at Urbana-Champaign (1975).
Cheung, Y. K., I. P. King and O. C. Zienkiewicz. “Slab Bridges with Arbitrary Shape and Support
Conditions – A General Method of Analysis Based on Finite Elements,” Proc. Inst. Civ. Eng., 40,
9-36 (1968).
Desai, C. S., and J. F. Abel. Introduction to the Finite Element Method – A Numerical Method
for Engineering Analysis. New York: Van Nostrand Reinhold Company (1972).
Dixon, J. D. “Analysis of Tunnel Support Structure with Consideration of Support-Rock Interac-
tion,” U.S. Dept. of Interior, Bureau of Mines Investigation, Report RI7526 (June 1971).
Hughes, T. J. R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis.
Prentice Hall (1987).
Lemos, J. V. “Modelling of Arch Dams on Jointed Rock Foundations,” in Prediction and Perfor-
mance in Rock Mechanics & Rock Engineering (Proceedings of ISRM International Symposium
EUROCK ‘96, Turin, Italy, September 1996), Vol. 1, pp. 519-526. G. Barla, ed. Rotterdam: A.
A. Balkema (1996).
Lemos, J. V., et al. “Experimental Study of an Arch Dam on a Jointed Foundation,” in Proceedings
of the Eighth International Congress on Rock Mechanics (Tokyo, Japan, September 1995), Vol.
3, pp. 1263-1266. T. Fujii, ed. Rotterdam: A. A. Balkema (1995).
Lorig, L. J. “A Hybrid Computational Model for Excavation and Support Design in Jointed Media.”
Ph.D. Thesis, University of Minnesota (1984).
Monsees, J. E. “Station Design for the Washington Metro System,” in Proceedings of the Engi-
neering Foundation Conference – Shotcrete Support. ACI Publication SP-54 (1977).
Paul, S. L., et al. “Design Recommendations for Concrete Tunnel Linings,” University of Illinois,
DOT Report No. DOT-TSC-UMTA-83-16 (1983).
Zienkiewicz, O. C. The Finite Element Method. London: McGraw-Hill Book Company (UK)
Limited (1977).