Antonios Vytiniotis Master Thesis 2009 PDF

Contents
1 Introduction ............................................................................................................................ 4
1.1 Liquefaction ...................................................................................................................... 4
1.2 Centrifuge Experiments.................................................................................................... 7
1.3 Pre-fabricated Vertical Drains .......................................................................................... 8
1.4 Objectives ......................................................................................................................... 9
2 Background ........................................................................................................................... 10
2.1 Centrifuge testing ........................................................................................................... 10
2.1.1 Basic Principles ........................................................................................................ 10
2.1.2 Scaling laws ............................................................................................................. 10
2.1.3 Other similitude problems ...................................................................................... 12
2.2 Analytical and numerical research towards modeling drainage-based soil improved
methods .................................................................................................................................... 14
3 Centrifuge Experiment .......................................................................................................... 21
3.1 Experiment outline ......................................................................................................... 21
3.2 Model Preparation ......................................................................................................... 22
3.3 Scale Factors ................................................................................................................... 23
3.4 Limitations of recorded data and model........................................................................ 24
4 Finite Element model ............................................................................................................ 34
4.1 Governing Equations ...................................................................................................... 34
4.2 Verification through a one-dimensional analytical solution .......................................... 36
4.2.1 Verification of the Analytical Solution .................................................................... 37
4.2.2 Comparison of the full formulation and the u-p formulation ................................ 38
4.2.3 Validation of OPENSEES .......................................................................................... 41
4.3 Constitutive soil model................................................................................................... 42
4.4 Non-linear geometry issues ........................................................................................... 52
4.5 Appendix ........................................................................................................................ 53
4.5.1 One-dimensional dynamic response of a fully saturated soil column ................... 53
1
4.5.2 Verification of Opensees ......................................................................................... 59
5 Drain Elements ...................................................................................................................... 71
5.1 Truss theory .................................................................................................................... 71
5.2 Darcy Weisbach Equation .............................................................................................. 73
5.3 Laminar Flow .................................................................................................................. 74
5.3.1 Flow equations ........................................................................................................ 74
5.3.2 Finite Element approximation ................................................................................ 75
5.3.3 Opensees Implementation ..................................................................................... 78
5.4 Turbulent Flow ............................................................................................................... 79
5.4.1 Flow equations ........................................................................................................ 79
5.4.2 Finite Element Approximation ................................................................................ 80
5.4.3 Opensees Implementation ..................................................................................... 85
5.5 Storage Capacity ............................................................................................................. 85
5.6 Scaling laws for PV-drains .............................................................................................. 86
5.6.1 Application of scaling laws to the SSKN01 centrifuge experiment......................... 93
5.6.2 Plane strain drain transmissivity properties ........................................................... 94
5.7 Appendix ........................................................................................................................ 96
5.7.1 One dimensional Problems ..................................................................................... 96
5.7.2 Plane Strain Consolidation ...................................................................................... 96
5.7.3 Laminar Drains Source Code ................................................................................... 99
5.7.4 Fully Turbulent Flow Drains Source Code ............................................................. 114
6 FEM Simulation ................................................................................................................... 131
6.1 Finite element model description ................................................................................ 131
6.1.1 Model Setup .......................................................................................................... 131
6.1.2 Selection of model parameters ............................................................................ 134
6.2 Base case analysis......................................................................................................... 141
6.2.1 Predicted Excess Pore Pressures........................................................................... 142
6.2.2 Predicted accelerations ........................................................................................ 143
6.2.3 Predicted displacements ....................................................................................... 144
2
6.3 Effect of different approximations in PV-drains simulations ....................................... 163
6.3.1 Drain resistance .................................................................................................... 163
6.3.2 Drain stiffness ....................................................................................................... 165
6.3.3 Drain storage capacity .......................................................................................... 170
6.3.4 Drain turbulence ................................................................................................... 175
7 Summary – Conclusions ...................................................................................................... 178
7.1 Simulating Vertical Drains ............................................................................................ 178
7.2 Similitude Issues ........................................................................................................... 179
7.3 Validation ..................................................................................................................... 179
7.4 Future research ............................................................................................................ 180
8 Bibliography ........................................................................................................................ 181
3
1 Introduction
1.1 Liquefaction
Seismic risk is a very significant design factor in many areas of the planet. Structures should be
designed in such a way so that both the superstructure and the foundation and soil beneath the
structure can withstand the expected cyclic loads. In this section, the problems of the substruc-
ture are discussed.
The most important problem of the response of saturated sandy soils under cyclic loading is
liquefaction. By liquefaction we refer to the situation where the effective stress of the soil is
Figure 1-1 Tilting of the Kawagishi-cho apartment buildings near the Shinano river bank due to the 1964 Niigata earthquake
4
reduced significantly because of the contractive tendency of medium and loose sands under
monotonic or cyclic loading. This situation leads to significant shear displacements, differential
displacements below structures, and eventually to significant damage, if the structure is not
designed against this occurring shear strength loss.
5
A very major incident, the Niigata earthquake (1964), brought liquefaction to the attention of
engineers. During this event many typical problems due to liquefaction have been observed
such as bearing capacity failures, lateral spreading, and sand boils. Very famous is the tilting of
the Kawagishi-cho apartment buildings near the Shinano river bank (Figure 1-1), which suffered
very small structural damage. The Kobe earthquake (1995) also increased awareness and
showed that an earthquake can significantly damage permanently the operations of a port. The
damages in the Kobe port were enormous (Figure 1-2).
Figure 1-2 Significant damage in the Kobe port facilities due to the 1995 Kobe earthquake
6
1.2 Centrifuge Experiments
Class A and Class B predictions (Lambe, 1973) in real-scale problems are very hard to perform
since even the exact input motion in an earthquake event is not known with great precision.
Also, Class C predictions usually also referred to as back analysis are mistrusted since one can
change the model parameters to match the predictions. Testing scale models is very important
in order to validate numerical models for geotechnical earthquake engineering. By means of
these tests engineers can have data to perform Class A and Class B predictions.
Figure 1-3 The UC-Davis geotechnical centrifuge
The shaking table has been used with great success for structural engineering problems. On the
other hand in geotechnical problems gravity is the most important external force and is not
7
modeled correctly in a scaled model, in which densities are kept constant but the linear dimen-
sions reduced(Zienkiewicz, Chan, Pastor, Schrefler, & Shiomi, 1999). The need to test our model
in a stress field similar to the real one led to into performing tests in a scale model in a centri-
fuge, similar to the one shown in Figure 1-3.
This need comes from the fact the response of the soil is path dependent and depends on the
initial stress state. Centrifuge testing facilities gave us the opportunity to develop and test with
great accuracy, new building improvements, numerical models, and eventually to improve
modern building codes and simulation capabilities.
1.3 Pre-fabricated Vertical Drains
Many liquefaction risk mitigation techniques are currently available. One of these remediation
Figure 1-4 Earthquake drain installation
8
techniques using drainage as a means of reducing liquefaction risk is the installation of PV-
drains (earthquake drains). PV drains are perforated, corrugated plastic pipes encased in a geo-
synthetic fabric (geo-textile), ranging from 75 to 200mm in diameter. The use of PV-drains is
very common in facilitating (accelerating) the consolidation of clay layers in earthwork con-
structions (e.g. levees, foundations). During an earthquake (or more generally, during cyclic
loading) the PV-drains installed in higher permeability sands can help increase drainage and
hence prevent the buildup of excess pore pressures that causes liquefaction and leads to large
ground movements, and significant damage to buildings above the liquefied layers.
1.4 Objectives
This project has five primary objectives:
• Derive an uncoupled theory for mechanical deformation and flow inside a prefabricated
vertical drain.
• Use this theory to discuss similitude issues between drains in model and prototype
scale.
• Implement the theory of mechanical and flow behavior of a drain using the FEM soft-
ware framework of Opensees.
• Validate the implemented theory with centrifuge experiments performed at UC-Davis
• Explain different assumptions that can be used for the design of the drains.
9
2 Background
2.1 Centrifuge testing
2.1.1 Basic Principles
Scale models inside a centrifuge are usually built either inside rigid boxes or inside laminar box-
es. The second one has the advantage that its sides can move like shear beams allowing to si-
mulate far-field response of soil with minimum interference. These models are instrumented
with wireless or wired sensors that, for geotechnical earthquake engineering purposes, meas-
ure pore pressure, acceleration, and displacements around the model. Video recording is also
widely used to examine phenomena happening during the tests (e.g. mechanisms of failure).
After the initial construction phase, and instrumentation installation, the scale model is in-
stalled on a shaking table, which in turn is installed on the arm of the centrifuge. Next, scale
model starts to rotate in a fairly high angular speed, which produces an almost uniform field of
usually up to 100g in the model area. Usually, some spin-ups are performed to test the model,
the instrumentation, and calibrate the sensors. Then, during a spin-up one can apply by means
of the shaking table an input motion to the model and record the output field of interest.
2.1.2 Scaling laws
In a centrifuge test, it is important to relate parameters measured in the scale model (M) to the
prototype (P) model, the real-scale model that is modeled. In the following we present the ba-
10
sic scaling laws. The equilibrium equations of the coupled pore-pressure displacement problem,
if we neglect the convective terms, for both the model and the prototype are:
(1.1)
(1.2)
Where σ is the stress matrix, x is the length, ρ is the material total density, b is the body force
vector, α is the acceleration of the soil skeleton, and w is the velocity of water relative to the
soil skeleton.
If we assume that the length is scaled by a factor of N:
(1.3)
then the acceleration is scaled:
(1.4)
Which can only be possible if the time is also scaled by:
(1.5)
The body forces
(1.6)
The average acceleration of water
11
(1.7)
The hydraulic conductivity is scaled by the following factor, in order to match diffusion:
(1.8)
Due to the scaling in permeability, in the centrifuge model we use either a different soil materi-
al, or a fluid with different viscosity but same density as water (e.g. silicon oil). Using a different
soil is problematic since it would be difficult to match it with the prototype stress-strain rela-
tion. Using a fluid with different viscosity is less problematic but still introduces small errors due
to different bulk modulus compared to the prototype fluid.
An important limitation of the centrifuge test is that it is unsuitable to model the response of
dynamic events under partially saturated conditions, since one need to use an alternate fluid in
order to achieve dynamic compatibility of diffusion and inertial behavior, something that
changes surface tension, drying, and wetting characteristics.
2.1.3 Other similitude problems
When modeling dynamic coupled pore pressure displacement problems other similitude prob-
lems that people should have in mind. Undrained shear strength of clayey soil has been ob-
served to increase about 5% to 15% for every log cycle of strain rate (Randolph, Cassidy,
Gourvenec, & Erbrich, 2005). In terms of the grain size effect, it is not well known how it effects
footing settlement even though a rule of thumb of
12
(1.9)
should eliminate these problems both for strip ( B=width) and circular (B=diameter) footings.
Also, for piles of diameter B there should be no significant grain size effects in case:
(1.10)
The development of shear banding seems to be significantly affected by the grain size, as
shown from model tests on a trap door and on a cavity collapse problem(Stone & Muir, 1992).
One should account for the non-uniformity of the centrifuge acceleration on the vertical stress
in soil samples(Schofield, 1980). The silo effect also can strongly reduce the vertical stress in
depth. Modern silo theories showing drastic dependence of the vertical stress field on the ratio
of the height and diameter of the box can be used(Garnier, 2001). The horizontal deflection of
the container wall can diminish the K0 value. In order for this to be avoided, one should try to
limit the horizontal deflections to less than H/200 (where H is the height of the utilized box)
(Ternet, 1999).
The homogeneity of the pluviated sand is also not an issue, since, by using appropriate tech-
niques and equipment, variation on dry unit weight within a container may be less than ±0.5%
(±0.1 kN/m3)(Garnier, 2001). As far as the in-flight in-situ CPT testing is concerned, we should
always make sure that the ratio S/B (S distance of the CPT to the nearest container wall, and B
the diameter of the cone) is greater than 10, and that the ratio D/B, where D is the diameter of
the container, is greater than 30:
13
(1.11)
(1.12)
Various other similitude aspects exist relative to pile driving, tunneling, erosion, sedimentation,
electro-osmosis, non-aqueous phase transport in soils and fractures, and heat transfer that are
very important but are not discussed in this thesis.
2.2 Analytical and numerical research towards modeling drainage-
based soil improved methods
No researcher has up to this point numerically simulated in 3-D the coupled pore-pressure dis-
placement response of soil layers improved with pre-fabricated vertical drains, because of the
complexity and computational cost. On the other hand, work has been performed using various
simplifications, such as reducing the modeling dimensions, simplifying the PDE’s to be solved,
and assuming infinite drain permeability.
14
Figure 2-1 Axisymmetric to plane strain equivalence, in the case the distance between the drains is the same in the 3D and in
the plane strain model
Hird et al, (1992) have investigated the use of two-dimensional, plane strain approximations of
consolidation for arrays of PV drains. The authors propose an equivalence that allows a true-
radial consolidation problem around a drain to be simulated with in plane strain, using an equal
strain and uniform vertical strain assumption. Their results enable selection of equivalent drain
spacing or equivalent soil permeability based on the same average degree of consolidation
within the soil layer. Their findings are also applicable for PV drains in high permeability sand
deposits. Their methodology is presented in Figure 2-1 for the case that the distance between
the drains is the same in the 3D real life scenario and in the 2D plane strain approximation.
Since, their relation keeps the same degree of consolidation in the two modeling spaces, it
should be noted that discrepancies on the exact pore pressure profiles are expected.
15
Figure 2-2 Design charts for groups of gravel drains against liquefaction (Seed & Booker, 1977)
Seed & Booker, (1977) were the first to analyze the role of PV drains for liquefaction risk mitiga-
tion. Their analyses consider the buildup and concurrent dissipation of pore pressures for drains
of infinite permeability. Their methodology presented in the form of design charts is explained
in Figure 2-2.
After the revolutionary work of Seed and Booker experimental results have been performed in
order to verify their method. Gravel drains have been created in a sandy soil to a depth of 11m,
with various spacing ratios of 0.25, 0.333, and 0.417. The experimental results in Figure 2-3 are
printed on top of the Seed and Booker design chart, and the read values from the class A pre-
diction are also shown.
16
An important conclusion is drawn: a drain’s finite permeability cannot be ignored, since even a
slight increase in well resistance can bring considerable increase in the pore pressure
ratio(Onoue, Mori, & Takano, In-situ experiment and analysis on well resistance of gravel
drains, 1987). Onoue (1988) has used the axi-symmetric diffusion equation, and an empirical
pore pre-
Figure 2-3 Read values from diagram presented by (Seed & Booker, 1977) and measured values of excess pore water pres-
sure ratio (Onoue, Mori, & Takano, In-situ experiment and analysis on well resistance of gravel drains, 1987)
ssure generation model, in order to include the effect of well resistance in design charts. It is
found that when the cycle ratio (number of cycles divided by the number of cycles to reach li-
quefaction) is close to one (1) then the greater the time the more significant the difference be-
tween the analytical results considering the vertical flow and those disregarding it. On the con-
trary, when the cycle ratio is equal or greater to two (2), then there is no practical difference
should we choose to consider vertical drainage or not. For this reason, for cycle ratio of 1
17
Figure 2-4 Relationship between coefficient of well resistance and spacing ratio
Onoue produced one chart ignoring vertical flow (Figure 2-5) and one without ignoring it
(Figure 2-4). For cycle ratios of 2,3, and 4 only one chart has been produced (Figure 2-4).
Similar work has also been performed, by using one-dimensional finite elements, with the same
axi-symmetric diffusion equation, and pore pressure generation model with Onoue, in order to
include also the storage effect of the drain(Pestana, Hunt, Goughnour, & Kammerer, Effect of
Storage Capacity on Vertical Drain Performance in Liquefiable Sand Deposits). In many cases,
18
the phreatic level of water inside a drain is not at the top of the drain. So, as the excess pore
pressure develops, the water level will first rise up to the top of the drain (water will be stored),
and then it will get out of the drain; Storage capacity is thus the amount of water that will be
stored in the drain before flow out of the drain starts. Analyses of the effect of storage capacity
have been performed, with perfect drains, drains with finite permeability, drains with variable
initial water level, and drains with presence of reservoir (storage capacity) of varying size. From
these analyses, a combined plot showing the effect of storage capacity drain permeability is
presented (Figure 2-6).
Figure 2-5 Relation between coefficient of well resistance and spacing ratio in the case where vertical direction de-watering
is disregarded
A two-dimensional finite difference code, FLAC, and an advanced elasto-plastic bounding sur-
face plasticity soil model(Andrianopoulos, 2006), have been utilized for the simulation of soil
improved with infinite permeability stone columns(Papadimitriou, Moutsopoulou, Bouckovalas,
& Brennan, 2007). Their results show that the design charts of Seed and Booker for thin layers
with insignificant drain resistance have inherent conservatism due to the fact that these charts
19
are based on undrained tests and ignore the drainage happening suddenly in the initial stage of
cycling. Of course this is counteracted from the unconservatism of not including the finite drain
permeability and the storage capacity.
Figure 2-6 Pore pressure ratio with storage and varying drain resistance, where ks is the permeability of the soil, kd is the
permeability of the drain. The drain spacing is s/d=5 and the water level inside the drain begins at 1m depth.
20
3 Centrifuge Experiment
3.1 Experiment outline
In order to evaluate the effectiveness of PV drains for liquefaction remediation (Earthquake
drains) two centrifuge tests were performed at the centrifuge facility at UC Davis, at March
2007 (SSK01) (Kamai, et al., 2008) and at October 2007(RNK01)(Kamai, et al., 2008) . They com-
pare the response of two similar facing slopes with a central channel. Beneath the left side
slope, lies a layer of loose sand containing an array of PV-drains, while beneath the right side
slope lies a layer of untreated loose sand. The two sides were symmetrically sloped 3° towards
a 200 mm wide central channel and were both comprised of three distinct layers (Figure 3-1):
(1) a bottom layer of dense Nevada Sand, overlain by (2) a liquefiable layer of loose Nevada
Sand, which was overlain by (3) a layer of compacted Yolo Loam. The Yolo Loam was used
to
Figure 3-1 Conceptual diagram of the PV drains centrifuge model
21
impede the vertical dissipation of pore water pressure out of the liquefiable layer. Shaking
events were applied to the model at a centrifugal acceleration of approximately 15 g. All events
were applied transverse to the river channel in the (longitudinal) direction. In the first (SSK01)
test the excitation is harmonic, and in the second test (RNK01) the excitation is an actual rec-
orded acceleration time history.
The tests have two important goals. The first is to be a proof-of-concept, showing that, in-
principle, PV drains can be effective in reducing liquefaction resistance. On the other hand the
most important scope of the test is to validate numerical models, which could be used next to
perform parametric studies or to be used directly by practitioners in design. Additionally the
second test (RNK01) compared to the first (SSK01) deals with also two new questions: investi-
gate the effect of the drain stiffness in the soil (soil pinning), and examine the difference in re-
sponse of the soil to a real earthquake excitation compared to the harmonic excitation in the
first test. In this thesis we concentrate on the first test SSK01 which due to the harmonic excita-
tion applied is better for validation of our numerical model.
3.2 Model Preparation
At first the PV-drains were prepared. They are nylon tubes (7mm inside diameter) in which
holes (1.5 mm diameter) were drilled every 5mm to allow for drainage of the excess pore pres-
sure. The tubes were wrapped twice with a Precision Woven Polypropylene Mesh, to keep sand
grains from clogging or entering the tube during liquefaction. The crust material was con-
structed using natural Yolo Loam, which was sun dried and sieve to pass a #10 sieve. Then wa-
ter was added to reach the optimum water content of this soil (15%).
22
The drains were at first placed on a triangular grid using wire as a support. Their placement po-
sition is shown in Figure 3-2, where one can see that small changes in their installation geome-
try have taken place during pluviation. Then a large box pluviator was used to construct to
create the layer of the dense (Dr=84%) sand, and a barrel pluviator was used to construct the
rest of the model of the loose sand layer (Dr=40%). The crust was placed in three layers.
Finally, the model was flooded with CO2 and placed under a vacuum of approximately 90kPa to
remove air within the soil. Then water was slowly dripped into saturation troughs, slowing satu-
rating the model from bottom to top. Saturation was targeted so that the entire liquefiable
layer would be saturated.
During construction of the scale model three types of sensors were installed. Pore pressure
transducers to monitor excess pore pressure (Figure 3-5,Figure 3-6), displacement transducers
(Figure 3-7,Figure 3-8) to measure horizontal and vertical displacements, and accelerometers
(Figure 3-3,Figure 3-4) to measure amplification or de-amplification of the input motion in vari-
ous locations. Instruments installed were both wired or wireless. From the above figures one
can see that especially on the top right side (the side that liquefied the most) the pore pressure
transducers and the accelerometers moved significantly downwards during shaking.
3.3 Scale Factors
The factors used to convert the data to prototype scale are indicated in Table 3-1.
Table 3-1 Scale Factors
Quantity Prototype/Model Dimension
23
Time 15/1
Displacement, Length 15/1
Acceleration, Gravity 1/15
Pressure, Stress 1/1
Permeability 15/1
3.4 Limitations of recorded data and model
The relative densities of the sand layers are assumed to be the same with the ones the pluvia-
tors were calibrated for. The relative density can also be estimated from the volume and weight
of each layer. Such calculations indicate that the loose sand layer may not have Dr=40% but
Dr=30%. Slight variations in the measured height of the layer (±2mm) will produce through large
variation in relative density (±10%). For this reason, the actual relative density was assumed to
be very close to the target, based on the pluviator calibrations prior to model preparation.
Apart from that, one has to keep in mind that on the top right side (the untreated side), pore
pressure transducers sunk significantly during liquefaction so their results can only be used if
one takes into consideration their vertical movement.
24
Figure 3-2 Location of the drains before and after pluviation
25
Figure 3-3 Accelerometer locations: cross section (as built and after test)
26
Figure 3-4 Accelerometer locations: plan view (as built and after shaking)
27
Figure 3-5 Pore pressure transducers: cross section (before and after shaking)
28
Figure 3-6 Pore pressure transducers: plan view (before and after shaking)
29
Figure 3-7 Displacement transducers: cross section
30
Figure 3-8 Displacement transducers: plan view
31
Figure 3-9 Initial and final deformed shape: cross section
32
Figure 3-10 Model pictures – from top left corner, clockwise: (1) drains placed before pluviation (2) surface markers on un-
treated side (3) sand boil – cross section through crust (4)&(5) the untreated side after the test – cracking and sand boils.
33
4 Finite Element model
This chapter presents description and verification for the numerical procedures used in Open-
sees.
4.1 Governing Equations
This section describes the equilibrium and mass conservation equations that are to be solved.
These equations are coupled and should be solved simultaneously. The equations that predict
the response of a soil fully saturated with fluid are (Biot, 1956, and Zienkiewicz et al, 1999):
Conservation of momentum:
(4.1)
Diffusion equation:
(4.2)
Conservation of mass:
(4.3)
Where S is the divergence matrix, σ is the stress matrix, ρ is the total density of the fluid, w is
the average (Darcy) velocity of the percolating water (thus relative to the soil skeleton), b are
the body forces, R are the viscous drag forces, ρf is the fluid density, n is the porosity, ε is the
strain. Q is defined by:
34
(4.4)
Where Kf is the fluid compressibility and Ks is the soil skeleton compressibility, and α has to do
with the definition of the effective stresses:
(4.5)
These equations turn out to be difficult to solve in their complete form. Simplifications are
possible, by assuming that the flow of water inside the soil skeleton is relatively slow. We can
thus neglect the relative acceleration of water to the soil skeleton. Under this assumption, we
get the following set of equations, which is referred to as the u-p approximation. This approxi-
mation can be used in consolidation problems, and in coupled dynamic pore pressure dis-
placement analyses with very small errors (Zienkiewicz, Chan, Pastor, Schrefler, & Shiomi,
1999).
(4.6)
(4.7)
Another simplification assumes that we can ignore only the convective terms of the fluid acce-
leration; this results in the so-called u-p-U approximation:
(4.8)
(4.9)
where U is the water displacement relative to the soil skeleton.
35
Current practice uses Flac, ABAQUS, and DYNAFLOW, OPENSEES for solution of advanced geo-
technical earthquake engineering problems. A brief summary of their capabilities and limita-
tions for the applications having to do with geotechnical earthquake engineering is provided in
Table 4-1. In this research we use Opensees, “an object-oriented software framework for simu-
lation applications in earthquake engineering using finite element methods”(Mazzoni,
McKenna, & Fenves, 2005), for its capabilities, modularity, and opensource development.
Table 4-1 Comparison between available software for geotechnical earthquake engineering simulations
Software Pros Cons
ABAQUS Implicit; Good pre- and post No dynamic coupled pore pressure-
u-p processing; Modularity; FE solver; displacement elems;
DYNAFLOW Implicit; Fully coupled elements; FE No pre/post-processing; No modularity;
u-p-U solver;
FLAC Un-coupled finite difference solver; Finite difference; Explicit code with less
u-p free field boundaries; Modularity; error controlling than implicit codes;
OPENSEES Fully coupled implicit FE solver; Ad- No pre/post-processing; Opensource;
u-p, u-p-U vanced constitutive models; Modulari-
ty; Opensource;
4.2 Verification through a one-dimensional analytical solution
In order to verify Opensees we compare the steady state response of the response of a satu-
rated soil column, where sinusoidal pressure and free drainage is applied on top. The stress
36
strain response is assumed to be linear. The geometry of the problem is presented in Figure
4-1.
4.2.1 Verification of the Analytical Solution
In Figure 4-2 we can see the results of the analytical solution for four (4) different angular fre-
quencies ω = 0.1, 1, 10, 100rad/s. From these plots we see that the solution follows the boun-
dary conditions. Also, the higher the frequency the smaller the displacement at the top; for
Figure 4-1 Geometry of the verification problem
lower frequencies the solution converges to the static solution. From the maximum pore water
pressure plot we can see that for low frequencies drainage makes the pore water pressure at
the bottom to be less than unity (where unity is the amplitude of the applied pressure on top).
37
Figure 4-2 Comparison of maximum displacement, displacement at t=5s, maximum excess pore water pressure and pore
water pressure for various angular frequencies
4.2.2 Comparison of the full formulation and the u-p formulation
In this section we compare the two analytical formulations, the full formulation accounting for
pore fluid acceleration terms and the convective terms, and the u-p formulation. A sample
analysis has been performed for a loading frequency of ω=10rad/s and permeability
k=0.001m/s, k = 0.1m/s, and k = 0.2m/s thus changing the effect of the pore water velocity from
negligible to important.
38
Figure 4-3 Comparison of the fully coupled solution and the u-p formulation for
k=0.001m/s
Figure 4-5 we see profiles of the absolute maximum displacement and pore water vs depth.
From these figures, we see that for a frequency ω=10rad/s, which is reasonable for geotechnic-
al earthquake engineering problems, the results of the two formulations are very close for a
39
permeability up to k=0.1m/s which is very large even for a loose sand. On the other hand, we
see that for k=0.2m/s differences start to become significant. Of course, one should note that
this is a p-wave propagation problem, and different behavior in 1-D shear wave strain propaga-
tion should be expected. In a shear wave propagation problem smaller discrepancies between
the two formulations should be expected, since the coupling between shear deformation and
volumetric response introduces a smaller driver force for the relative to the soil skeleton
movement of water, than the direct effect a p-wave has on the displacement of water around
the soil mass.
Figure 4-4 Comparison of the fully coupled solution and the u-p formulation for k=0.1m/s
40
4.2.3 Validation of OPENSEES
In this section we are testing OPENSEES by comparing its results to an analytical solution. We
use the reference analysis described on the previous section with ω=10rad/s and k=0.2m/s
since in this case the fully coupled and the u-p formulation differ noticeably. The results are
presented in Figure 4-6. We can see that the pore pressure profile calculated through Opensees
with quad-up elements matches greatly the analytical solution for the u-p formulation. In terms
of the vertical displacement vs depth there are minor discrepancies, but still the Opensees solu-
tion is closer to the analytical u-p solution, and also the displacements are so small (of the order
of 10-4) that a small numerical error could cause this small difference.
41
Figure 4-6 Verification of the coupled pore pressure displacement solver in Opensees, snapshot at maximum displacement
on top
4.3 Constitutive soil model
In order to accurately predict the cyclic response of the sand deposits inside the laminar box
one needs to use an accurate constitutive soil model. In this analysis two soil models have been
used, both of them implemented in the standard version of OPENSEES. The constitutive soil
model for the sand is the PressureDependMultiYield02 (Yang, Elgamal, & Parra, 2003) and the
model for the Yolo Loam is the PressureIndependMultiYield (Mazzoni et al, 2005).
The pressudependmultiyield02 model is a multi-yield surface elasto-plastic model, based on
similar work from Prevost (1985). The model has parameters that do not depend on the stress
state, but depend on the voids ratio. It can simulate shear-induced volume contraction of
42
dilation, cyclic mobility, and liquefaction (e.g. undrained behavior of sands or silts under cyclic
loading). The pressurerindependmultiyield02 model is an elasto-plastic model. The volumetric
stress strain reponse is linear-elastic. Plasticity occurs only in the deviatoric stress-strain
response and is insensitive to confinement (e.g. undrained behavior of clays).
In Figure 4-17 the geometry of a simple verification problem is presented. It consists of one
quad-up element, a plane strain element solving the coupled pore-pressure displacement using
the u-p approximation. The element is a four noded element using bilinear isoparametric
formulation. Each element node has 2 degrees of freedom (DOFs) for the displacement of the
soil skeleton and 1 DOF for pore pressure. The bottom nodes are fixed. The top nodes are tied
to have the same vertical and horizontal displacements. A cyclic shear load is applied at the top
of the element with a natural frequency of ω=1rad/s. In Table 4-2,
Table 4-3 and Table 4-4 we see the parameters used for these models(Mazzoni, McKenna, &
Fenves, 2005). The permeability of the soil is selected as k=3x10-5m/s, the top of the element is
drained and the bottom is impermeable. The test includes also inertia effects according to the
density of each layer. With this choice of parameters there is partial drainage during the test.
The results are presented in the form of timehistories of the horizontal and vertical stresses,
stress paths, and stress-strain curves. By comparing Figure 4-8 to Figure 4-10 with Figure 4-11
to Figure 4-13 we see that with the same parameters the model is able to capture the stress
level dependency of the shear behavior. By comparing Figure 4-11 to Figure 4-13 with Figure
4-14 to Figure 4-16 we see that the model predicts different behavior for different sand densi-
ties by means of changing the model parameters. In Figure 4-17 to Figure 4-19 we see the be-
43
havior of the pressure-independent material, which behaves as almost perfectly plastic after
the level of stress exceeds the cohesion of the layer. We note that at every simulation the con-
stitutive model predicts that the lateral earth pressure ratio becomes very close to 1.
Figure 4-7 Geometry of the Cyclic Simple Shear Test
Table 4-2 Yolo Loam model Properties
PressureIndependMultiYield (Soft Clay)
Ρ 1.3
Gref 13000
Kref 65000
C 18.0
γpeak 0.1
Table 4-3 Dense Sand model properties
PressureDependendMultiyield02 (for a dense sand, Dr=80%)

Ρ 2.07 cpressuredependence 0.5
44
Gref 130000 ψphasetransformation 26.0
Kref 260000 c1 0.013
Φ 36.5 c3 0.0
γpeak 0.1 d1 0.3
pref 80 d3 0.0
Table 4-4 Medium-loose sand model properties
PressureDependendMultiyield02 (for a medium-loose sand, Dr=40%)

Ρ 1.98 cpressure dependence 0.5
Gref 90000 ψphase transformation 26.0
Kref 220000 c1 0.067
Φ 32.0 c3 0.23
γpeak 0.1 d1 0.06
pref 80 d3 0.27
45
Vertical Effective Stress vs Time (Soils Convention)
Vertical Stress (kPa)
60
40
20
0 5 10 15
Time (s)
Horizontal Effective Stress vs Time (Soils Convention)
Horizontal Stress (kPa)
60
40
20
0 5 10 15
Time (s)
f
Figure 4-8 Vertical and Horizontal Effective Stress timehistories for a DSS test of a Medium-Loose Sand for σ'v,initial=50kPa
Stress Path (Soils convention) Stress Path (Soils Convention)

40 40
(σ 'v-σ 'h)/2) (kPa)
(σ 'v-σ 'h)/2) (kPa)
20 20
0 0
-20 -20
-40 -40
0 50 -40 -20 0 20 40
(σ 'v+σ 'h)/2) (kPa) τ xy (kPa)
Figure 4-9 Stress paths for a DSS test of a Medium-Loose Sand for σ'v,initial=50kPa
46
Stress Strain Curve Stress Path (Soils Convention)
40 40
20 20
τ xy (kPa)
τ xy (kPa)
0 0
-20 -20
-40 -40
-0.02 -0.01 0 0.01 0.02 0 50
Strain Mean effective Stress (kPa)
Figure 4-10 Stress path and stress strain curve for a DSS test of a Medium-Loose Sand for σ'v,initial=50kPa

60
40
20
0
0 5 10 15
Time (s)
60
40
20
0
0 5 10 15
Time (s)
Figure 4-11 Vertical and Horizontal Effective Stress timehistories for a DSS test of a Medium-Loose Sand for σ'v,initial=20kPa
47
40
(σ 'v-σ 'h)/2) (kPa)
40
(σ 'v-σ 'h)/2) (kPa)

20 20
0 0
-20 -20
-40 -40
0 20 40 -40 -20 0 20 40
Figure 4-12 Stress paths for a DSS test of a Medium-Loose Sand for σ'v,initial=20kPa

40 40
20 20
τ xy (kPa)
τ xy (kPa)
0 0
-20 -20
-40 -40
-0.1 0 0.1 0.2 0.3 0 20 40
Figure 4-13 Stress path and stress strain curve for a DSS test of a Medium-Loose Sand for σ'v,initial=20kPa
48
Vertical Stress (kPa) Vertical Effective Stress vs Time (Soils Convention)
60
40
20
0 5 10 15
Time (s)
60
40
20
0 5 10 15
Time (s)
Figure 4-14 Vertical and Horizontal Effective Stress timehistories for a DSS test of a Dense Sand for σ'v,initial=50kPa

40 40
(σ 'v-σ 'h)/2) (kPa)
(σ 'v-σ 'h)/2) (kPa)
20 20
0 0
-20 -20
-40 -40
0 50 -40 -20 0 20 40
Figure 4-15 Stress paths for a DSS test of a Dense Sand for σ'v,initial=50kPa
49
40 40
20 20
τ xy (kPa)
τ xy (kPa)
0 0
-20 -20
-40 -40
-4 -2 0 2 4 0 50
Strain -3 Mean effective Stress (kPa)
x 10
Figure 4-16 Stress path and stress strain curve for a DSS test of a Dense Sand for σ'v,initial=50kPa

50
40
30
0 5 10 15
Time (s)
50
40
30
0 5 10 15
Time (s)
Figure 4-17 Vertical and Horizontal Effective Stress timehistories for a DSS test of a Soft Clay for σ'v,initial=50kPa
50
40 40
(σ 'v-σ 'h)/2) (kPa)
(σ 'v-σ 'h)/2) (kPa)

20 20
0 0
-20 -20
-40 -40
0 20 40 -40 -20 0 20 40
(σ 'v+σ 'h)/2) (kPa) τ (kPa)
xy
Figure 4-18 Stress paths for a DSS test of a Soft Clay for σ'v,initial=50kPa

40 40
20 20
τ xy (kPa)
τ xy (kPa)
0 0
-20 -20
-40 -40
-0.5 0 0.5 1 1.5 0 20 40
Figure 4-19 Stress path and stress strain curve for a DSS test of a Soft Clay for σ'v,initial=50kPa
51
4.4 Non-linear geometry issues
A very important feature of Opensees we need to include into our analysis when modeling a
centrifuge test is a compression only connection with the laminar box that contains the soil. For
this reason Opensees has an elastic-no-tension uniaxial material, that can be used with any li-
near (e.g. truss) or zerolength element (element that represents force deformation relationship
between two nodes at the same coordinates). This material gives a convergence error upon
tension since the implicit module of the code is used. For this reason, it should always be
coupled in parallel with an elastic material with very small elastic modulus, when convergence
problems are encountered. It should be noted that the definition of no-tension is easy when we
use it to define force deformation properties in a truss. If we define a zero-length element to
have an elastic-no-tension connection, then the order of the nodes when defining the zero-
length element actually affects which direction is constrained and which is not, because in prin-
ciple in this type of element tension and compression cannot be defined.
Figure 4-20 Elastic-no-tension deformation-force curve
52
4.5 Appendix
4.5.1 One-dimensional dynamic response of a fully saturated soil column
The reference problem chosen is finding the steady state dynamic response of a fully saturated
soil column in which free drainage and a sinusoidal pressure are applied at the top. The materi-
al is linear elastic and the pore water is considered to be incompressible.
4.1
4.2
4.3
4.4
4.5
4.6
where σ is the vertical stress, σ’ is the effective stress, ε is the vertical strain, u is the vertical
displacement, D is the one-dimensional compression modulus, ρ is the density of the total
composite, ρf is the density of the fluid phase, w is the pore water displacement relative to the
soil skeleton, Kf is the compressibility of the fluid phase, p is the pore water pressure and n is
the porosity. From these equations we can get a system of ordinary differential equations of u
and w:
4.7
53
4.8
Next we apply separation of variables to sole the above system. In the case of a sinusoidal exci-
tation:
4.9
the solution of the σ, u, w variables is of the form of
4.10
where X is a space function. Using the above we have:
4.11
4.12
4.13
4.14
Note that the results will be complex numbers. We define:
4.15
4.16
4.17
4.18
54
4.19
4.20
4.21
4.22
The previous equations take the following form:
4.23
4.24
The above equations can be written in a canonical form:
4.25
4.26
And we define:
4.27
4.28
4.29
55
4.30
The solution of the system of PDE’s is:
4.31
4.32
Where λi’s are the solutions of the characteristic equation:
4.33
Next, we apply the boundary conditions. At the surface:
4.34
4.35
4.36
And at the bottom,
4.37
4.38
4.39
The boundary conditions give us the following system of equations that should be solved:
4.40
4.41
56
4.42
4.43
Finally, we need to estimate p, the pore pressure vs depth:
4.44
4.45
4.46
4.47
4.48
In order to find a solution to u-p approximation of this problem, we omit the relative accelera-
tion of the pore fluid to the soil skeleton:
4.49
4.50
The solution is the same as above, except for the fact than A, B, Γ, and Δ are defined:
4.51
4.52
4.53
57
4.54
58
4.5.2 Verification of Opensees
4.5.2.1Comparing full to u-p formulations

In this section we present some extra figures for comparison of the full to the u-p formulation.
59
60
4.5.2.2Matlab Code Solving Analytically the full formulation of the dynamic response of
a soil column
function [u max_u p max_p w H L]=coupl_1d(t)
% This function solves analytically the 1d coupled (pore
% pressure-displacement) dynamic response problem of a soil column.
% t is the time of interest
% It uses the real part of exp(i*omega*t) so the excitation is a cosine
K_f=2200000; % Volumetric compressibility of the fluid

n=0.333; % Porosity
E=30000; % Elastic modulus of the soil skeleton
v=0.3; % Poisson’s ratio
rho_f=1.; % Fluid density
rho=2.; % Average density of multi-phase medium
L=10; % Height of soil column
omega=10.; % Natural Frequency of the Applied Load
q=1; % Amplitude of the Applied Load
g=10; % Acceleration of gravity
kappa=0.2; % Permeability (Hydraulic Conductivity)
rho_dry=rho-n*rho_f
e=n/(1-n)
D_oned=E*(1-v)/((1+v)*(1-2*v))
k=(K_f/n)/(D_oned+K_f/n)
V_c2=(D_oned+K_f/n)/rho;
beta=rho_f/rho
sqrt(V_c2)
T=2*pi/omega
T_star=2*L/sqrt(V_c2);
Pi_1=(2/beta/pi)*kappa*T/g/(T_star^2)
Pi_2=pi^2*(T_star/T)^2
% Now we solve the system of differential equations:

% Look at the theory
% The equivalent equations are:
% d2u/dz2=A*u+B*w
% d2w/dz2=C*u+D*w
% The general solution is:
% u=Ci*b*e^(lambda_i*z) (Einstein summation convention)
A=(beta*Pi_2-Pi_2)/(1-k);
B=(beta/n*Pi_2-i/Pi_1-beta*Pi_2)/(1-k);
C=-beta*Pi_2-k/(1-k)*(beta*Pi_2-Pi_2);
C=C/k;
D=-beta/n*Pi_2+i/Pi_1-k/(1-k)*(beta/n*Pi_2-i/Pi_1-beta*Pi_2);
D=D/k;
% Now we create the characteristic polynomial

% The solution of the characteristic polynomial are the lambda's
61
P_char=[1 0 -(A+D) 0 (A*D-B*C)];
lambda=roots(P_char);
if (A*D-B*C)==0
'Beware A*D-B*C=0'
end
if (A-D)^2+4*B*C==0
'Beware (A-D)^2+4*B*C=0'
end
% We solve the system L_m*X=R

% The solution of this system is the C_i's
L_m=[exp(lambda(1)) exp(lambda(2)) exp(lambda(3)) exp(lambda(4));...
(lambda(1)^2-A)*exp(lambda(1)) (lambda(2)^2-A)*exp(lambda(2)) (lamb-
da(3)^2-A)*exp(lambda(3)) (lambda(4)^2-A)*exp(lambda(4));...
lambda(1) lambda(2) lambda(3) lambda(4);...
(lambda(1)^2-A)*lambda(1) (lambda(2)^2-A)*lambda(2) (lambda(3)^2-
A)*lambda(3) (lambda(4)^2-A)*lambda(4)];
R=[0 ; 0; q*L/D_oned/B; -q*L/D_oned];

X=L_m\R;
% t=1 % The absolute time
n_inc=1000;
L_inc=L/n_inc;
z=0;
i_n=0;
H=0;
% We calculate the displacements

while z<=1
i_n=i_n+1;
z=z+L_inc/L;
u(i_n)=0;
for i_it=1:4
u(i_n)=u(i_n)+X(i_it)*B*exp(lambda(i_it)*z);
end
temp_u=u(i_n);
u(i_n)=abs(u(i_n));
max_u(i_n)=u(i_n);
u(i_n)=u(i_n)*real(exp(i*(omega*t-phase(temp_u))));
H(i_n)=z*L;
end
% We calculate the fluid displacement

z=0;
i_n=0;
while z<=1
i_n=i_n+1;
z=z+L_inc/L;
w(i_n)=0;
for i_it=1:4
w(i_n)=w(i_n)+X(i_it)*(lambda(i_it)^2-A)*exp(lambda(i_it)*z);
62
end
temp_w=w( i_n);
w(i_n)=abs(w(i_n));
w(i_n)=w(i_n)*real(exp(i*(omega*t-phase(temp_w))));
end
% We calculate the pore pressure

z=0;
i_n=0;
while z<=1
i_n=i_n+1;
z=z+L_inc/L;
p(i_n)=0;
for i_it=1:4
p(i_n)=p(i_n)+omega^2*rho_f*X(i_it)*B/lambda(i_it)*(exp(lambda(i_it)*z)-
1)+...
(rho_f*omega^2/n-rho_f*g*i*omega/kappa)*X(i_it)*(lambda(i_it)^2-
A)/lambda(i_it)*...
(exp(lambda(i_it)*z)-1);
end
p(i_n)=p(i_n)*L;
temp_p=p(i_n);
p(i_n)=abs(p(i_n));
max_p(i_n)=p(i_n);
p(i_n)=p(i_n)*real(exp(i*(omega*t-phase(temp_p))));
end
63
4.5.2.3Matlab Code Solving Analytically the u-p formulation of the dynamic response of
a soil column
function [u max_u p max_p w H L]=coupl_1d_up(t)
% This function solves analytically the 1d coupled (pore
% pressure-displacement) dynamic response problem of a soil column.
% This function follows the u-p formulation ignoring the acceleration terms
% for the pore fluid.
% t is the time of interest
% It uses the real part of exp(i*omega*t) so the excitation is a cosine
K_f=2200000; % Volumetric compressibility of the fluid

n=0.333; % Porosity
E=30000; % Elastic modulus of the soil skeleton
v=0.3; % Poisson’s ratio
rho_f=1; % Fluid density
rho=2.; % Average density of multi-phase medium
L=10; % Height of soil column
omega=10.; % Natural Frequency of the Applied Load
q=1; % Amplitude of the Applied Load
g=10; % Acceleration of gravity
kappa=0.2; % Permeability (hydraulic conductivity)
rho_dry=rho-n*rho_f
e=n/(1-n)
D_oned=E*(1-v)/((1+v)*(1-2*v))
k=(K_f/n)/(D_oned+K_f/n)
V_c2=(D_oned+K_f/n)/rho
beta=rho_f/rho
sqrt(V_c2)
T=2*pi/omega
T_star=2*L/sqrt(V_c2);
Pi_1=(2/beta/pi)*kappa*T/g/(T_star^2)
Pi_2=pi^2*(T_star/T)^2
% Now we solve the system of differential equations:

% Look at the theory
% The equivalent equations are:
% d2u/dz2=A*u+B*w
% d2w/dz2=C*u+D*w
% The general solution is:
% u=Ci*b*e^(lambda_i*z) (Einstein summation convention)
A=(beta*Pi_2-Pi_2)/(1-k);
B=(-i/Pi_1)/(1-k);
C=-beta*Pi_2-k/(1-k)*(beta*Pi_2-Pi_2);
C=C/k;
D=i/Pi_1-k/(1-k)*(-i/Pi_1);
D=D/k;
64
% Now we create the characteristic polynomial
% The solution of the characteristic polynomial are the lambda's
P_char=[1 0 -(A+D) 0 (A*D-B*C)];
lambda=roots(P_char);
if (A*D-B*C)==0
'Beware A*D-B*C=0'
end
if (A-D)^2+4*B*C==0
'Beware (A-D)^2+4*B*C=0'
end
% We solve the system L_m*X=R

% The solution of this system is the C_i's
L_m=[exp(lambda(1)) exp(lambda(2)) exp(lambda(3)) exp(lambda(4));...
(lambda(1)^2-A)*exp(lambda(1)) (lambda(2)^2-A)*exp(lambda(2)) (lamb-
da(3)^2-A)*exp(lambda(3)) (lambda(4)^2-A)*exp(lambda(4));...
lambda(1) lambda(2) lambda(3) lambda(4);...
(lambda(1)^2-A)*lambda(1) (lambda(2)^2-A)*lambda(2) (lambda(3)^2-
A)*lambda(3) (lambda(4)^2-A)*lambda(4)];
R=[0 ; 0; q*L/D_oned/B; -q*L/D_oned];

X=L_m\R;
% t=1 % The absolute time
n_inc=1000;
L_inc=L/n_inc;
z=0;
i_n=0;
H=0;
% We calculate the displacements

while z<=1
i_n=i_n+1;
z=z+L_inc/L;
u(i_n)=0;
for i_it=1:4
u(i_n)=u(i_n)+X(i_it)*B*exp(lambda(i_it)*z);
end
temp_u=u(i_n);
u(i_n)=abs(u(i_n));
max_u(i_n)=u(i_n);
u(i_n)=u(i_n)*real(exp(i*(omega*t-phase(temp_u))));
H(i_n)=z*L;
end
% We calculate the fluid displacement

z=0;
i_n=0;
while z<=1
i_n=i_n+1;
z=z+L_inc/L;
w(i_n)=0;
65
for i_it=1:4
w(i_n)=w(i_n)+X(i_it)*(lambda(i_it)^2-A)*exp(lambda(i_it)*z);
end
temp_w=w( i_n);
w(i_n)=abs(w(i_n));
w(i_n)=w(i_n)*real(exp(i*(omega*t-phase(temp_w))));
end
% We calculate the pore pressure

z=0;
i_n=0;
while z<=1
i_n=i_n+1;
z=z+L_inc/L;
p(i_n)=0;
for i_it=1:4
p(i_n)=p(i_n)+omega^2*rho_f*X(i_it)*B/lambda(i_it)*(exp(lambda(i_it)*z)-
1)+...
(-rho_f*g*i*omega/kappa)*X(i_it)*(lambda(i_it)^2-
A)/lambda(i_it)*...
(exp(lambda(i_it)*z)-1);
end
p(i_n)=p(i_n)*L;
temp_p=p(i_n);
p(i_n)=abs(p(i_n));
max_p(i_n)=p(i_n);
p(i_n)=p(i_n)*real(exp(i*(omega*t-phase(temp_p))));
end
66
4.5.2.4Opensees tcl/tk Code to test the u-p approximation
# This analysis solve the dynamic response of a soil column
# upon the application of a constant load on top
# The problem is saturated
wipe
model BasicBuilder -ndm 2 -ndf 3
# Below we define variables for the analysis

set E 30000
set nu 0.3
set numXele 1; # number of elements in x (H) direction
set numYele 40; # number of elements in y (V) direction
set xSize .5; # Element size in x direction
set ySize .25; # Element size in z direction
set numXnode [expr $numXele+1]
set numYnode [expr $numYele+1]
set smass 2.;
set peak_shear_strain 10000.
set c 100.
set G [expr $E/(2*(1+$nu))]
set B [expr $E/(3*(1-2*$nu))]
set i 1
set j 1
set pi 3.141593
# Define material
nDMaterial PressureIndependMultiYield 1 2 $smass $G $B $c $peak_shear_strain\
0. 100. 0. 1
# Define nodes
for {set i 1} {$i <= $numXnode} {incr i 1} {
for {set j 1} {$j <= $numYnode} {incr j 1} {
set xdim [expr ($i-1)*$xSize]
set ydim [expr ($j-1)*$ySize]
set nodeNum [expr $i + ($j-1)*$numXnode]
node $nodeNum $xdim $ydim
}
}
# define elements
set k 0.2
set k [expr $k/10/1.] ;#actual value used in computation
set gravX 0.0
set gravY 0.0
set press 0.0
set bulk_f 2.2e6
set n_por 0.333
set bulk [expr $bulk_f/$n_por]
for {set i 1} {$i <= $numXele} {incr i 1} {
for {set j 1} {$j <= $numYele} {incr j 1} {
set eleNum [expr $i + ($j-1)*$numXele]
set n1 [expr $i + ($j-1)*$numXnode]
set n2 [expr $i + ($j-1)*$numXnode + 1]
set n4 [expr $i + $j*$numXnode + 1]
set n3 [expr $i + $j*$numXnode]
element quadUP $eleNum $n1 $n2 \
$n4 $n3 1.0 1 $bulk 1. $k $k $gravX $gravY $press
}
}
67
#vees
# Fix Base:
fix $i 1 1 0
}
# Fix X Direction: (To all but the top nodes)

for {set i [expr $numXnode+1]} {$i <= [expr $numXnode*($numYnode-1)]} {incr i 1} {
fix $i 1 0 0
}
# Drainage on top: (And x-Fixity)

fix [expr $i + ($numYnode-1)*$numXnode] 1 0 1;
}
set omega 10.
# Calculate the period from the timeseries

set T [expr 2*$pi/$omega]
set Timeseries "Sine 0.0 50. $T -shift [expr $pi/2]"
pattern Plain 1 $Timeseries {
load [expr $i + ($numYnode-1)*$numXnode] 0. -.25 0.;
}
}
#build
recorder Node -file output_disp.txt -time \
-node 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39\
41 43 45 47 49 51 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 -dof 1 2 disp
recorder Node -file porepress.txt -time \
-node 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39\
41 43 45 47 49 51 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 -dof 3 vel
recorder Node -file output_react.txt -time -node 1 2 -dof 1 2 reaction
#recorder Element -file output_stress.txt material 1 stiffness

set gamma 0.5
test NormDispIncr 1.0e-5 10 0;
algorithm Newton
integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4] \
0.00 0.0 0.002 0.0
analysis Transient
set startT [clock seconds]
analyze 500 0.01
set endT [clock seconds]
puts "Execution time: [expr $endT-$startT] seconds."
68
4.5.2.5Opensees tcl/tk code to test the constitutive model
model BasicBuilder -ndm 2 -ndf 3
node 1 0.000 0.000
node 2 1.000 0.000
node 3 0.000 1.000
node 4 1.000 1.000
# Depending on the material one should comment/uncomment the respective material

# Dense
# nDMaterial PressureDependMultiYield02 1 2 2.1 130000 260000 36.5 0.1 101. .5 26.
0.013 0.0 0.3 0.0
# Loose
# nDMaterial PressureDependMultiYield02 1 2 1.7 60000 160000 31. 0.1 101. .5 31. 0.087
0.18 0. 0.0
# Medium-Loose
# nDMaterial PressureDependMultiYield02 1 2 1.8 90000 220000 32. 0.1 101. .5 26. 0.067
0.23 0.27 0.77
# Soft Clay
nDMaterial PressureIndependMultiYield 1 2 2.1 13000 65000 18. 0.1
# Define Elements
set gravX 0.0
set gravY 0.
set press 0.
set bulk_f 2.2e9
set n_por 0.333
set bulk [expr $bulk_f/$n_por]
element quadUP 1 3 1 2 4 1.0 1 $bulk 1. [expr 0.00003/9.81/1.] [expr 0.00003/9.81/1.]
$gravX $gravY fix 1 1 1 0
fix 2 1 1 0
fix 3 0 0 1
fix 4 0 0 1
equalDOF 3 4 1 2
# Depending on the excitation one should change the Timeseriesini factor

set Timeseriesini "Constant -factor 50."
# pattern UniformExcitation 1 1 -accel $Timeseries;

pattern Plain 1 $Timeseriesini {
load 4 0. -.5 0
load 3 0. -.5 0
}
# Set material to elastic for gravity loading

updateMaterialStage -material 1 -stage 0
# GRAVITY APPLICATION (elastic behavior)

set gamma 1.5
# create the SOE, ConstraintHandler, Integrator, Algorithm and Numberer

integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4] \
0.00 0.0 0.002 0.0
constraints Transformation
algorithm Newton
numberer RCM
system ProfileSPD
analysis Transient
analyze 3 5.e5
puts "End of Elastic Phase of Gravity Application"
69
# Set material to elasto-plastic for the rest of the loading
updateMaterialStage -material 1 -stage 1
analyze 5 5.e5
puts "End of Gravity Application. Starting Dynamic Excitation..."
# rezero time
wipeAnalysis
setTime 0.0
set omega 1.
set pi 3.141593
# Calculate the period from the timeseries

set T [expr 2*$pi/$omega]
set Timeseries "Sine 0.0 50. $T -shift 0.0 -factor 20."
pattern Plain 2 $Timeseries {
load 4 1. 0. 0
load 3 1. 0. 0
}
set gamma 0.6

integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4] 0.00 0.0 0.002 0.0
constraints Transformation
algorithm KrylovNewton
numberer RCM
system ProfileSPD
# create the analysis object

analysis Transient
# Recorders
recorder Node -file output_disp.txt -time -dof 1 2 disp
recorder Node -file output_pore.txt -time -dof 3 vel
recorder Element -file stress_1.txt -time -eleRange 1 1 material 1 stress
recorder Element -file strain_1.txt -time -eleRange 1 1 material 1 strain
# perform analysis
set startT [clock seconds]
analyze 1500 0.01
set endT [clock seconds]
puts "Execution time: [expr $endT-$startT] seconds."
70
5 Drain Elements
In this section we present the implementation of drain elements. These elements uncouple the
flow of water inside them and the mechanical deformation of the pipe material. Two element
are discussed, the first predicting laminar flow inside the drain, and the second predicting fully
turbulent flow inside the drains. In terms of the mechanical part the element acts as a truss.
These elements would allow engineers to study the effect of installation of finite-permeability
drains in soil layers both in situations involving laminar (slow) flow inside the drains (e.g. consol-
idation in clays), and in situations involving turbulent (fast) flow inside the drains (e.g. fast flow
during a seismic event in earthquake drains). These finite elements allow for taking into ac-
count the effect of storage capacity. This effect is very important to accurately predict the re-
sponse of sand layers improved with earthquake drains, when they are below non-permeable
clay layers. Finally, the mechanical effects of soil pinning around the PV drains, causing localiza-
tions, can also be studied.
5.1 Truss theory
A truss finite element is being used to predict the mechanical part of the drain element. Simple
truss theory is very well established in a finite element context. The definitions for this element
are presented in Figure 5-1. We define the force vector in global coordinates:
5.1
71
and the displacements vector:
5.2
The stiffness matrix for this truss element is:
5.3
With the above definitions the equilibrium is defined as:
5.4
Figure 5-1 Coordinates, displacements, and forces definitions for truss element
72
5.2 Darcy Weisbach Equation
The Darcy – Weisbach equation is a widely used phenomenological equation used in hydraulics.
It relates pressure loss due to friction along a given length of a pipe to the average velocity of
the fluid flow. The pressure loss can be calculated:
5.5
Where λ is a dimensionless coefficient of laminar or turbulent flow, L is the length of the pipe, D
is the diameter of the pipe, ρ is the density of the water, V is the average velocity of the flow. λ
is equivalent to the Darcy friction factor (f) and can be estimated for both laminar and turbulent
flow from the Moody diagram.
Figure 5-2 Estimation of the Darcy friction factor for laminar and turbulent flow
73
5.3 Laminar Flow
5.3.1 Flow equations
Laminar flow occurs when a fluid flows smoothly or in regular paths. In laminar flow, sometimes
called streamline flow, the velocity, pressure, and other flow properties at each point in the flu-
id remain constant. Laminar flow over a horizontal surface may be thought of as consisting of
thin layers, or laminae, all parallel to each other. Inside a pipe, the fluid in contact with the pipe
is stationary, but all the other layers slide over each other.
In laminar flow, it is possible to simplify the Darcy-Weisbach equation to a linear relation be-
tween the pressure loss and the average velocity of water (or flow). For laminar flow:
5.6
where Re is the Reynolds number.
For a circular tube:
5.7
where ν is the kinematic viscosity of water.
Using the above equations:
5.8
5.9
74
5.10
5.11
5.12
5.13
where μ is the dynamic viscosity of water (or any fluid in general), and
5.14
5.3.2 Finite Element approximation
Using this formulation we can define a one-dimensional element in a two dimensional space
relating water pressure to flow. We first define the vector of flow (equivalent to the vector of
external forces in the truss element):
5.15
And the vector of pressures:
5.16
By using equation 5.10 we can define a transmissivity matrix:
5.17
Then the equilibrium equation in this element is defined as:
75
Figure 5-3 Coordinates, pressure, and flow definitions for one-dimensional pipe element
5.18
Next, we assume uncoupled behavior between the mechanical behavior of the truss and the
flow taking place inside the pipe. We now define the generalized force vector in global coordi-
nates:
5.19
and the generalized displacements vector:
5.20
By combining the previous forms we have the generalized stiffness matrix:
76
5.21
Where:
5.22
5.23
5.24
5.25
The generalized equilibrium has the form:
5.26
Figure 5-4 Coordinates, pressure, flow, force, and displacement definitions for the uncoupled drain element
77
5.3.3 Opensees Implementation
The finite element is implemented in the Opensees Software framework, to be used together
with quad_up elements, in order to predict drainage inside soil layers by means of drains(e.g.
PV-drains or stone columns). A new command is added in the interpreter that takes the argu-
ments:
element Pipelin2 eleid node1 node2 Material Area Cl γw
eleid is the id of the element, a unite integer number assigned to this element, node1 and
node2 are the start and end nodes of the element, Material is a an integer number already de-
fined assigning a specific constitutive material for the mechanical response, Area is the area of
the drain that contributes to the mechanical behavior, Cl is defined:
5.27
Units of the Cl are:
5.28
where L is length, F is the force, T is the time
Finally, γw is the density of the water times the acceleration of gravity (e.g. negative when
pointing downwards).
78
5.4 Turbulent Flow
5.4.1 Flow equations
Turbulent flow is a flow regime characterized by chaotic property changes. For fully turbulent
flow according to the Moody Diagram (Figure 5-2) λ is a constant. From the Darcy-Weisbach
equation:
5.29
5.30
5.31
5.32
5.33
which can be written in a rate form:
5.34
5.35
79
5.4.2 Finite Element Approximation
One can define a one-dimensional element in a two dimensional space relating water pressure
to flow, with the same conventions used in Figure 5-3. In the finite element approximation of
the turbulent flow regime for the drains we are using a consistent Jacobian formulation.
We first define the vector of flow (equivalent to the vector of external forces in the truss ele-
ment):
5.36
5.37
5.38
And the vector of pressures:
5.39
5.40
5.41
During the nth step of the integration:
5.42
And at the n+1th step:
80
5.43
So:
5.44
We also define Δp1(n) and Δp2(n) as the increments of pore pressure in node 1 and 2 of the drain
element respectively. Also:
5.45
Using this we can write the consistent transmissivity matrix:
5.46
5.47
This formulation has a disadvantage. Δp(n) might be very low or zero, and this could cause nu-
merical errors, or no convergence. When Δp(n) is very small then we can instead use the conti-
nuum Jacobian. By using equation 6.35 we can define a continuous transmissivity matrix:
5.48
81
5.49
As we can see a problem still rises in the continuum Jacobian matrix, when ΔP is very close to
zero, or else when i, the hydraylic gradient, is very small. Remembering that we only need the
continuous Jacobian when the Δp(n) is very close to zero then, we derive one more scheme to
be used numerically when both (a) i is very small (b) Δp(n) is very small. Under these circums-
tances we linearize the Q vs i equation, and we assume a linear region of size 2dc.
5.50
Which similar to the laminar flow drain gives:
5.51
Next, we assume uncoupled behavior between the mechanical behavior of the truss and the
flow taking place inside the pipe. We now define the generalized force vector in global coordi-
nates:
5.52
82
5.53
and the generalized displacements vector:
5.54
5.55
The generalized equilibrium has the form:
5.56
By combining the previous forms we have the generalized stiffness matrix:
5.57
Where:
5.58
83
5.59
5.60
And k33 is when Δp(n)>dc:
5.61
When Δp(n)<dc and i>dc:
5.62
And Δp(n)<dc and i<dc:
5.63
In Figure 5-5 a summary of the used approximations and regimes is shown.
Figure 5-5 Various regimes and definitions used to integrate the Darcy-Weisbach equation for fully turbulent flow
84
5.4.3 Opensees Implementation
The finite element is implemented in the Opensees Software framework, to be used together
with quad_up elements, in order to predict drainage inside soil layers by means of drains of any
type. A new command is added in the interpreter that takes the arguments:
element Pipe3 eleid node1 node2 Material Area Ct γw dc
eleid is the id of the element, a unite integer number assigned to this element, node1 and
node2 are the start and end nodes of the element, Material is a an integer number already de-
fined assigning a specific constitutive material for the mechanical response, Area is the area of
the drain that contributes to the mechanical behavior, Ct is:
5.64
The units of Ct are defined:
5.65
And γw is the density of the water times the acceleration of gravity (e.g. negative when pointing
downwards) and dc defines the initial linear gradient regime, or the size of any gradient regime
that can be considered almost zero.
5.5 Storage Capacity
In typical applications the water level inside the drain is not at the surface. Thus, the water level
inside it first has to increase, before allowing dissipation of excess pore pressure. This is a signif-
85
icant effect that could reduce the effectiveness and the applicability of the installation of PV
drains as a liquefaction risk mitigation technique.
Since, the drain elements calculate flow inside the drain, caused by the pore pressure the
quad_up elements around the drain apply, one can keep track the volume of water inside the
pipe. A varying water level inside the pipe means a variable boundary condition at the top of
the pipe. The effect of storage capacity can be modeled by simply keeping track of the volume
of water inside the drain and changing the boundary condition on the top of the drain. This
procedure has the disadvantage of assuming perfect drainage inside the part of the drain that
stores water, but this error is insignificant compared to (1) added accuracy due to the finite
permeability of the rest of the drain, and (2) added accuracy due to precise estimation of the
storage capacity of the drain.
5.6 Scaling laws for PV-drains
The flow coming out a drain in the prototype is related to the flow coming out of the model
drain in the following way:
5.66
So for laminar flow:
5.67
5.68
86
5.69
5.70
And for turbulent flow:
5.71
5.72
5.73
5.74
Eventually, we conclude that the parameter relating the pore pressure gradient to the flow of
water has to be N3 times larger in the prototype scale when we deal with laminar flow drains
and N2.5 times larger when we deal with turbulent flow. This is consistent with well known scal-
ing laws for flow inside the soil mass. For laminar flow the permeability in the prototype should
be N times larger and for turbulent flow the permeability should stay the same. These results
are easily concluded should one use the scaling laws for drains and insert an equivalent per-
meability for the laminar drain (to see that this changes by a factor of N for the prototype) and
the kinematic viscosity for the turbulent drains (to see that this does not changes as we scale up
the model). From the above analysis a problem arises having to do with the Reynolds number:
5.75
87
5.76
5.77
5.78
5.79
So a flow that is laminar in the model scale might be turbulent in the prototype scale. This is
important because the experiments are not representing reality correctly and they are not on
the safe side. Using a different fluid to match diffusion makes the situation much worse since
then:
5.80
When designing a scale model with drains one should design the transmissivity of the scale
drains so that it matches the extra resistance caused by the turbulence in the real-scale drains,
since it is very common that the model scale flow is mostly laminar and the prototype scale
flow is mostly turbulent. A methodology is proposed at this point that would allow selection of
pipes for scale modeling purposes.
In prototype scale we need to find a laminar flow parameter Cl that gives the smallest error
compared with the flow coming out of a drain with parameter Ct, for a specific range of i (pore
pressure gradient). We define a function that is the square of the difference of the laminar flow
calculations to the turbulent flow calculations.
88
5.81
5.82
5.83
Next we define the integral of function f:
5.84
5.85
5.86
In order for the flow calculation for turbulent flow to match the calculation for the laminar flow
we need to minimize F.
5.87
5.88
5.89
5.90
5.91
If we replace the expressions for Cl and Ct we have:
89
5.92
5.93
5.94
5.95
5.96
5.97
The above diameter is the design diameter for the model drains, given the parameters of the
prototype drains. This relation has been produced under the assumption that the flow is lami-
nar in the model scale and turbulent in the model scale. This allows for the model scale laminar
flow to be as similar as it can be to the prototype scale turbulent flow by minimizing the
squares of the distances between the flow calculated by the two theories. This relation can be
inverted, in order to examine what is the prototype drain that the model drain is representing.
5.98
Also, a reasonable range for PV drains that would not make us loose accuracy in the range of
small gradient, but sufficient to capture large gradients would be:
90
5.99
A reasonable value, considering the fact that in a normal setup very few excess pore pressures
will develop along the line of the drain would be:
5.100
5.101
So the aforementioned relationships simplify as:
5.102
5.103
Also, one needs to notice that a fluid of different viscosity might be used to match diffusion be-
tween the scale model and the prototype, according to the following equation:
5.104
So we re-write the equations in the following form:
5.105
91
5.106
5.107
5.108
So, under the assumption that:
5.109
We have:
5.110
5.111
5.112
5.113
92
5.6.1 Application of scaling laws to the SSKN01 centrifuge experiment
In this section we use the aforementioned scaling laws for the centrifuge experiment in order to
examine similitude issues between the model and prototype drains. We know for our case that:
5.114
5.115
5.116
5.117
(from Moody’s diagram for ε/d=2.4*10-5)
5.118
5.119
5.120
Which gives, using eqn. 6.91:
5.121
5.122
Which in this circumstance is very similar to the prototype dimension of the drains
(Dp=15*0.007m=0.105m), so we will not expect large discrepancies for this specific problem for
using either the laminar flow or the turbulent flow equation. It is though a little bit larger, which
means that under this range of i's the assumption of full turbulence in the drains makes them
93
less permeable. On the other hand, the drains used in this experiment have a very smooth sur-
face on the contrary to real drains that might have a non-smooth surface due to bacteria
growth and material deposition. Also, the maximum gradient used here might not be always
that small depending on the permeability of the drains, the drain spacing, and the permeability
of the soil. These factors can increase significantly the prototype dimension that the drains cor-
respond to. Thus, there should be caution when we use results from the centrifuge experiment
to interpret what happens in reality. The turbulent flow solution since this gives more accurate
pore pressures, which are needed to predict correctly the soil response since it is path depen-
dent.
5.6.2 Plane strain drain transmissivity properties
Following Hird’s solution, the axi-symmetric transmissivity of a laminar flow drain is scaled, in
order to model the effect of the drain in a plane strain model, given that we keep the drain
spacing the same between the 3D and 2D models:
5.123
where R is the drain spacing in 3D.
The above equation has been derived from solving analytically the consolidation problem of a
plane strain unit cell with a laminar flow drain. Solving analytically the consolidation problem of
a plane strain unit cell with a fully turbulent flow drain is much harder. Thus, in order to employ
a plane strain model to model the axisymmetric real life scenario we go back to eqn. 5.76. For a
given fully turbulent flow drain we can find an equivalent laminar flow drain. We can convert
the laminar flow drain transmissivity to the plane strain one and then convert back to the fully
94
turbulent flow drain transmissivity. Since Cl and Ct are linearly dependent on each other, then
the same scaling factor as in equation 5.123 applies to fully turbulent flow drains:
5.124
Where w is the width of the plane strain finite element grid.
95
5.7 Appendix
5.7.1 One dimensional Problems
For verification purposes tests on the element level have been performed using pressure con-
trol and flow control, using one (1) and two (2) elements. The simulations performed are sum-
marized on the next table.
Table 5-1 Validation tests for laminar and turbulent flow drains
Pore Pressure Control Flow Control
One element
• •
Two Elements
• •
5.7.2 Plane Strain Consolidation
Also the consolidation of a plane strain unit cell, with elastic material governing the response of
the soil skeleton has been examined. The geometry used is shown in Figure 5-6. The boundaries
are all impermeable except for the top of the drain. In Figure 5-7 we can see results from this
simulation. These results compare the directly calculated vertical settlement at the top of the
drain to the indirectly calculated vertical settlement at the top of the drain, which was esti-
mated by integrating the flow of water coming out of the drain, denoting basically a mass equi-
librium inside the soil mass. We can see for the laminar drain that the directly calculated vertic-
96
al displacement is very close to the indirectly calculated one. The directly calculated vertical
displacement is a bit more than the indirectly calculated one due some flow coming out from
the soil mass rather than from the drain. On the other hand, for the turbulent drains we can see
that the results match very closely each other, this time because the permeability of the drain is
very large, and almost all of the flow comes out of drain (very small flow is coming out of the
soil mass). Under the current implementation, the turbulent flow drains appear to work much
better using the Krylov subspace acceleration technique, since other methods are shown to
leave a small value of flow at the end of the analysis.
Figure 5-6 Geometry of the plane strain consolidation verification problem. Zero pore pressure boundary condition is only
applied at point A.
97
Figure 5-7 Validation for laminar and fully turbulent flow drains
98
5.7.3 Laminar Drains Source Code
5.7.3.1Class Definition
/* ****************************************************************** **
** OpenSees - Open System for Earthquake Engineering Simulation **
** Pacific Earthquake Engineering Research Center **
** **
** **
** (C) Copyright 1999, The Regents of the University of California **
** All Rights Reserved. **
** **
** Commercial use of this program without express permission of the **
** University of California, Berkeley, is strictly prohibited. See **
** file 'COPYRIGHT' in main directory for information on usage and **
** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. **
** **
** Developed by: **
** Frank McKenna ([email protected]) **
** Gregory L. Fenves ([email protected]) **
** Filip C. Filippou ([email protected]) **
** **
** ****************************************************************** */
// $Revision: 1.00 $
// $Date: 2008/07/18 18:05:53 $
// $Source: /usr/local/cvs/OpenSees/SRC/element/pipe/Pipe.h,v $
// Written: Antonios Vytiniotis

// Created: 07/08
// Revision: A
//
// Description: This file contains the definition for the Pipelin2. A
Pipelin2 object
// provides the abstraction of the small deformation bar element plus
predicts the
// uncoupled pore pressure change according to Darcy Weisbach equation for
laminar flow.
// Each pipe object is associated with a material object dealing with the
axial compressibility
// of the drain. This Pipelin2 element will work in 2d problems in a 3DOF
domain.
//
// What: "@(#) Pipelin2.h, revA"
#ifndef Pipelin2_h
#define Pipelin2_h
#include <Element.h>
99
#include <Matrix.h>
class Node;
class Channel;
class UniaxialMaterial;
class Pipelin2:public Element {

public:
//constructors
Pipelin2 (int tag, int Nd1, int Nd2, UniaxialMaterial &theMaterial,
double A, double C_3, double Grav=0.0);
Pipelin2();
//destructor
~Pipelin2();
//public methods to obtain information about dof & connectivity

int getNumExternalNodes(void) const;
const ID &getExternalNodes(void);
int getNumDOF(void);
Node **getNodePtrs(void);
//public methods to set the state of the element

void setDomain(Domain *theDomain);
int commitState(void);
int revertToLastCommit(void);
int revertToStart(void);
int update(void);
//public methods to obtain stiffness, mass, damping, and residual

information
const Matrix &getTangentStiff(void);
const Matrix &getInitialStiff(void);
const Matrix &getDamp(void);
const Matrix &getMass(void);
void zeroLoad(void);
int addLoad(ElementalLoad *theLoad, double loadFactor);
int addInertiaLoadToUnbalance(const Vector &accel);
const Vector &getResistingForce(void);
const Vector &getResistingForceIncInertia(void);
//public methods for output

int sendSelf(int commitTag, Channel &theChannel);
int recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker
&theBroker);
int displaySelf(Renderer &theViewer, int displayMode, float fact);
void Print(OPS_Stream &s, int flag=0);
Response *setResponse(const char **argv, int argc, OPS_Stream &s);
int getResponse(int responseID, Information &eleInformation);
//protected:
private:
//private member function - only availabe to objects of the class
double computeCurrentStrain(void) const;
//private attributes - a copy for each object of the class
100
UniaxialMaterial *theMaterial; //pointer to a material
ID externalNodes; // contains the id's of end nodes
Matrix trans; //hold the transformation matrix
// Vector *theLoad; // pointer to the load vector P
double L; //length of Pipe based on undeformed configuration

double C_3;
double A;
double d_y_class;
int eletag;
double Gamma; //weight per unit volume
Node *end1Ptr, *end2Ptr; //two pointer to the trusses nodes
Node *theNodes[2]; //two pointer to the trusses nodes
in a matrix form (AV)
//private class attribute

static Matrix trussK;
static Matrix trussD;
static Matrix trussM;
static Vector trussR;
};
#endif
5.7.3.2Class Implementation
/* ****************************************************************** **
** **
** **
** **
** **
** Developed by: **
** **
** ****************************************************************** */
// $Date: 2008/07/18 18:05:53 $
// $Source: /usr/local/cvs/OpenSees/SRC/element/Pipe/Pipelin2.cpp,v $

// Created: 07/08
// Revision: A
//
// Description: This file contains the implementation for the Pipelin2 class.
101
//
#include "Pipelin2.h"
#include <Information.h>
#include <Parameter.h>
#include <Domain.h>
#include <Node.h>
#include <Channel.h>
#include <FEM_ObjectBroker.h>
#include <UniaxialMaterial.h>
#include <Renderer.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <ElementResponse.h>
#include <Matrix.h>
#include <Vector.h>
#include <ElasticMaterial.h>
// initial the class wide variables

Matrix Pipelin2::trussK(6,6);
Matrix Pipelin2::trussM(6,6);
Matrix Pipelin2::trussD(6,6);
Vector Pipelin2::trussR(6);
Pipelin2::Pipelin2(int tag,
int Nd1, int Nd2,
UniaxialMaterial &theMat,
double a, double c3, double g)
:Element(tag,ELE_TAG_Pipelin2),
theMaterial(0),
externalNodes(2),
trans(1,4),L(0.0), A(a), C_3(c3), Gamma(g),
end1Ptr(0), end2Ptr(0), eletag(tag),
d_y_class(0.0)
{
//create a copy of the material object
theMaterial=theMat.getCopy();
//fill in the ID containing external node info with node id's

externalNodes(0)=Nd1;
for (int i=0; i<2; i++)
theNodes[i] = 0;
trussR.Zero();
}
//constructor which should be invoked by an FE_ObjectBroker only

Pipelin2::Pipelin2()
:Element(0,ELE_TAG_Pipelin2),
102
theMaterial(0),
externalNodes(2),
trans(1,4), L(0.0), A(0.0), C_3(0.0), Gamma(0.0),
end1Ptr(0), end2Ptr(0),
d_y_class(0.0)
{
for (int i=0; i<2; i++)
theNodes[i] = 0;
}
Pipelin2::~Pipelin2()
{
if (theMaterial !=0)
delete theMaterial;
}
int Pipelin2::getNumExternalNodes(void) const

{
return 2;
}
const ID &
Pipelin2::getExternalNodes(void)
{
return externalNodes;
}
int
Pipelin2::getNumDOF(void){
return 6;
}
Node **
Pipelin2::getNodePtrs(void)
{
return theNodes;
}
void
Pipelin2::setDomain(Domain *theDomain)
{
//first ensure nodes exist in Domain and set the node
pointers
int Nd1 =externalNodes(0);
end1Ptr =theDomain->getNode(Nd1);
theNodes[0] = theDomain->getNode(Nd1);
if (theNodes[0]==0)
return;
if (theNodes[1]==0)
return;
// call the DomainComponent class method
103
this->DomainComponent::setDomain(theDomain);
//ensure connected nodes have corrent number of dof's

int dofNd1=theNodes[0]->getNumberDOF();
if ((dofNd1 !=3) || (dofNd2 !=3))
return; //don't go any further otherwise segmentation
fault
//now determine the length & transformation matrix

const Vector &end1Crd=theNodes[0]->getCrds();
double dx= end2Crd(0)-end1Crd(0);

double dy= end2Crd(1)-end1Crd(1);
d_y_class=dy;
L=sqrt(dx*dx+dy*dy);
if (L==0.0)
return;
double cs=dx/L;
double sn=dy/L;
trans(0,0)=-cs;
trans(0,1)=-sn;
trans(0,2)= cs;
trans(0,3)= sn;
// // determine the nodal mass for lumped mass approach

// M=M*A*L/2; //M was set to rho by the constructor
}
int
Pipelin2::commitState()
{
return theMaterial->commitState();
}
int
Pipelin2::revertToLastCommit()
{
return theMaterial->revertToLastCommit();
}
int
Pipelin2::revertToStart()
{
return theMaterial->revertToStart();
}
int
Pipelin2::update()
{
//determine the current strain given trial displacements at
nodes
104
double strain=this->computeCurrentStrain();
//set the strain in the materials

theMaterial->setTrialStrain(strain);
return 0;
}
const Matrix &

Pipelin2::getTangentStiff(void)
{
if (L==0) {//if length ==zero - we zero and return
trussK.Zero();
return trussK;
}
//get the current E from the material for the strain that
was set
// at the material when the update() method was invoked
double E = theMaterial->getTangent();
//form the tangent stiffness matrix

Matrix K_temp(4,4);
K_temp =trans^trans; //This is a temporary matrix
containing the truss stiffness parameters
K_temp *=A*E/L;
// trussK.Zero();
// Truss stiffness components:
trussK(0,0)=K_temp(0,0);
trussK(2,2)=0.;
trussK(5,5)=0.;
trussK(5,2)=0.;
trussK(2,5)=0.;
return trussK;
}
const Matrix &
Pipelin2::getInitialStiff(void)
{
if (L==0) {
trussK.Zero();
return trussK;
}
//get the current strain from the material

double strain = theMaterial->getStrain();
105
//get the current stress from the material
double stress = theMaterial->getStress();
//compute the tangent

double E=stress/strain;

Matrix K_temp(4,4);
K_temp *=A*E/L;
// trussK.Zero();
trussK(2,2)=0.;
trussK(5,5)=0.;
trussK(5,2)=0.;
trussK(2,5)=0.;
// opserr << "Componenents of the matrix K11" << trussK(3,3)

<<endln;
return trussK;
}
const Matrix &

Pipelin2::getDamp(void)
{
//No damping associated with this type of element
trussD.Zero();
trussD(2,2)=-C_3/L;
trussD(5,5)=-C_3/L;
trussD(5,2)=C_3/L;
trussD(2,5)=C_3/L;
double deleteme=C_3/L;
double deleteme2=C_3/L;
return trussD;
}
const Matrix &

Pipelin2::getMass(void)
{
if (L==0){
trussM.Zero();
return trussM;
}
// At this point we have zero lumped mass

trussM.Zero();
return trussM;
106
}
void
Pipelin2::zeroLoad(void)
{
//does nothing - no element load associated with this
object
}
int
Pipelin2::addLoad(ElementalLoad *theLoad, double
loadFactor)
{
opserr <<"MyTruss::addLoad - load type unknown for truss with
tag: " << this->getTag() << endln;
return -1;
}
int
Pipelin2::addInertiaLoadToUnbalance(const Vector &accel)
{
return 0;
}
const Vector &

Pipelin2::getResistingForce()
{
if (L==0) {//if length ==zero - zero and return
trussR.Zero();
return trussR;
}
// R=Ku-Pext
//force =F*transformation
double force = A* theMaterial->getStress();
trussR(0)= trans(0,0)*force;
const Vector &vel1 = theNodes[0]->getTrialVel();

// This is the linear element with total disp (no need for state
params)
// Domain::update(double a, double b);
// double dt=;
//
trussR(2)=trussD(2,2)*vel1(2)+trussD(2,5)*vel2(2)+trussD(2,2)*d_y_class
*Gamma;
// trussR(5)=trussD(5,2)*vel1(2)+trussD(5,5)*vel2(2)-
trussD(5,5)*d_y_class*Gamma;
trussR(2)=-C_3/L*vel1(2)+C_3/L*vel2(2)-
C_3/L*d_y_class*Gamma;
trussR(5)=C_3/L*vel1(2)-
C_3/L*vel2(2)+C_3/L*d_y_class*Gamma;
return trussR;
107
}
const Vector &

Pipelin2::getResistingForceIncInertia()
{
this->getResistingForce();
//No inertia is included in the in this element formulation

return trussR;
}
int
Pipelin2::sendSelf (int commitTag, Channel &theChannel)
{
int dataTag=this->getDbTag();
// Pipelin2 packs it's data into a Vector and sends this to

theChannel
//along with it's dbTag and the commitTag passed in the
arguments
Vector data(6);
data(0)= this->getTag();
data(1)=A;
data(4)=C_3;
data(5)=Gamma;
data(2)=theMaterial->getClassTag();
int matDbTag=theMaterial->getDbTag();
if (matDbTag==0) {
matDbTag =theChannel.getDbTag();
if (matDbTag !=0)
theMaterial->setDbTag(matDbTag);
}
data(3)=matDbTag;
theChannel.sendVector (dataTag, commitTag, data);
theChannel.sendID(dataTag, commitTag, externalNodes);
theMaterial->sendSelf(commitTag, theChannel);
return 0;
}
int
Pipelin2::recvSelf(int commitTag, Channel &theChannel,
FEM_ObjectBroker &theBroker)
{
int dataTag= this->getDbTag();
Vector data(6);
theChannel.recvVector(dataTag, commitTag, data);
this->setTag((int)data(0));
A=data(1);
C_3=data(4);
Gamma=data(5);
108
theChannel.recvID(dataTag, commitTag, externalNodes);
int matClass=data(2);
int matDb = data(3);
theMaterial= theBroker.getNewUniaxialMaterial(matClass);
theMaterial->setDbTag(matDb);
theMaterial->recvSelf(commitTag, theChannel, theBroker);
return 0;
}
int
Pipelin2::displaySelf(Renderer &theViewer, int displayMode,
float fact)
{
const Vector &end1Crd= end1Ptr->getCrds();
const Vector &end1Disp=end1Ptr->getDisp();
Vector v1(3);
Vector v2(3);
for (int i=0; i<2;i++) {
v1(i)=end1Crd(i)+end1Disp(i)*fact;
}
if (displayMode==3) {
//use the strain as the drawing measure
return theViewer.drawLine(v1, v2, strain, strain);
}
else if (displayMode==2){
//otherwise use the material stress
double stress =A*theMaterial->getStress();
return theViewer.drawLine(v1,v2,stress, stress);
}
else{
//use the axial force
double force = A*theMaterial->getStress();
return theViewer.drawLine(v1,v2,force,force);
}
}
void
Pipelin2::Print(OPS_Stream &s, int flag)
{
//compute the strain and axial force in the member
double strain, force;
if (L==0.0) {
strain=0;
force=0.0;
}
else{
strain = theMaterial->getStrain();
force=A*theMaterial->getStress();
}
109
const Vector &vel1 = theNodes[0]->getVel();

trussR(2)=-C_3/L*vel1(2)+C_3/L*vel2(2)-
C_3/L*d_y_class*Gamma;
trussR(5)=C_3/L*vel1(2)-
C_3/L*vel2(2)+C_3/L*d_y_class*Gamma;
if (flag==0) {//print everythin

s<< "Element: " <<this->getTag();
s<< " type: My Truss iNode: "<< externalNodes(0);
s<< " jNode: "<<externalNodes(1);
s<< " Area: "<< A;
if (Gamma!=0) s << "Gamma: "<<Gamma;
s<< " \n\t strain: " <<strain;

s<< " axial load: " <<force;
s<< " \n\t unbalanced load: " <<trussR;
s<< " \t Material: " << *theMaterial;
s<< endln;
} else if (flag==1) {//just print ele id, strain and force
s<< this->getTag() << " " <<strain << " " << force
<<endln;
}
}
Response *
Pipelin2::setResponse(const char **argv, int argc ,
OPS_Stream &s)
{
// we compare arg(0) for known response types for the Truss
//axial force
if(strcmp(argv[0], "axialForce")==0)
return new ElementResponse(this, 1, 0.0);
//a material quantity

else if (strcmp(argv[0], "material")==0)
return theMaterial->setResponse(&argv[1], argc-1, s);
else
return 0;
}
int
Pipelin2::getResponse(int responseID, Information
&eleInformation)
{
switch (responseID){
case -1:
return -1;
110
case 1:
return eleInformation.setDouble (A*theMaterial-
>getStress());
default:
return 0;
}
}
double
Pipelin2::computeCurrentStrain(void) const
{
//determine the strain
const Vector &disp1=end1Ptr->getTrialDisp();
double dLength=0.0;
for (int i=0;i<2;i++)
dLength -= (disp2(i)-disp1(i))*trans(0,i);
double strain =dLength/L;
return strain;
}
5.7.3.3Tcl/tk command interpreter

/* ****************************************************************** **
** **
** **
** **
** **
** Developed by: **
** **
** ****************************************************************** */
// $Revision: 1. $
// $Date: 2008/07/20 19:20:46 $
// $Source: /usr/local/cvs/OpenSees/SRC/element/pipe/TclPipelin2Command.cpp,v
$
// File: ~/element/TclPipelin2Command.C
//
// Written: avytin
// Created: 09/08
111
// Revision: A
//
// Description: This file contains the implementation of the
TclModelBuilder_Pipelin2()
// command.
//
// What: "@(#) TclModelBuilder.C, revA"
#include <stdlib.h>
#include <string.h>
#include <Domain.h>
#include "Pipelin2.h"
#include <TrussSection.h>
#include <TclModelBuilder.h>
#include <CorotTruss.h>
#include <CorotTrussSection.h>
extern void printCommand(int argc, TCL_Char **argv);
int
TclModelBuilder_Pipelin2(ClientData clientData, Tcl_Interp *interp, int argc,
TCL_Char **argv, Domain*theTclDomain,
TclModelBuilder *theTclBuilder,
int eleArgStart){
//make sure at least one other
argument to contain type of system
if (argc!=8 && argc!=9){
interp->result = "WARNING bad
command - Pipelin2 eleId iNode jNode matID Area c_3 Gamma";
return TCL_ERROR;
}
//get the id, x_loc, y_loc

int trussId, iNode, jNode, matID;
double A, C_3, Gamma=0.0;
if (Tcl_GetInt(interp,argv[2],
&trussId)!= TCL_OK){
interp->result = "WARNING
invalid eleId - Pipelin2 eleId iNode jNode matID Area c_3 Gamma";
return TCL_ERROR;
}
if (Tcl_GetInt(interp, argv[3],
&iNode) != TCL_OK) {
invalid iNode - Pipelin2 eleId iNode jNode matID Area c_3 Gamma";
return TCL_ERROR;
}
&jNode) != TCL_OK) {
invalid jNode - Pipelin2 eleId iNode jNode matID Area c_3 Gamma";
return TCL_ERROR;
}
112
&matID) != TCL_OK) {
invalid matID - Pipelin2 eleId iNode jNode matID Area c_3 Gamma";
return TCL_ERROR;
}
if (Tcl_GetDouble(interp, argv[6],
&A) != TCL_OK) {
invalid Area - Pipelin2 eleId iNode jNode matID Area c_3 Gamma";
return TCL_ERROR;
}
&C_3) != TCL_OK) {
invalid C_3 - Pipelin2 eleId iNode jNode matID Area c_3 Gamma";
return TCL_ERROR;
}
&Gamma) != TCL_OK) {
invalid C_3 - Pipelin2 eleId iNode jNode matID Area c_3 gamma";
return TCL_ERROR;
}
UniaxialMaterial *theMaterial =
theTclBuilder->getUniaxialMaterial(matID);
if (theMaterial ==0) {
opserr << "WARNING
TclPipelin2 - Pipelin2 - no Material found with tag ";
opserr << matID << endln;
return TCL_ERROR;
}
//now create the truss and add it

to the domain
Element *theTruss = 0;
theTruss=new
Pipelin2(trussId,iNode,jNode,*theMaterial,A,C_3,Gamma);
if (theTruss==0) {
opserr << "WARNING
TclPipelin2 - Pipelin2 - ran out of memory for node ";
opserr << trussId << endln;
return TCL_ERROR;
}
if (theTclDomain-
>addElement(theTruss)==false) {
delete theTruss;
opserr << "WARNING
TclPipelin2 - Pipelin2 - could not add Pipelin2 to the domain";
return TCL_ERROR;
}
113
//Everything is OK
return TCL_OK;
}
5.7.4 Fully Turbulent Flow Drains Source Code
5.7.4.1Class Implementation
/* ****************************************************************** **
** **
** **
** **
** **
** Developed by: **
** **
** ****************************************************************** */
// $Date: 2008/07/18 18:05:53 $
// $Source: /usr/local/cvs/OpenSees/SRC/element/pipe/Pipe.h,v $

// Created: 07/08
// Revision: A
//
// Description: This file contains the definition for the Pipe3. A Pipe3
object
// provides the abstraction of the small deformation bar element plus
predicts the
// uncoupled pore pressure change according to Darcy Weisbach equation. Each
pipe
// object is associated with a material object dealing with the axial
compressibility
// of the drain. This Pipe3 element will work in 2d problems in a 3DOF
domain.
//
// What: "@(#) Pipe3.h, revA"
#ifndef Pipe3_h
#define Pipe3_h
114
#include <Element.h>
#include <Matrix.h>
class Node;
class Channel;
class UniaxialMaterial;
//#define ELE_TAG_MyTruss 4002
// This is a trial implementation of a simle 2-d Pipe3 element

class Pipe3:public Element {
public:
//constructors
Pipe3 (int tag, int Nd1, int Nd2, UniaxialMaterial &theMaterial, double
A, double C_3, double Gamma=0.0, double D_C=0.0);
Pipe3();
//destructor
~Pipe3();
//public methods to obtain information about dof & connectivity

int getNumExternalNodes(void) const;
const ID &getExternalNodes(void);
int getNumDOF(void);
Node **getNodePtrs(void);
//public methods to set the state of the element

void setDomain(Domain *theDomain);
int commitState(void);
int revertToLastCommit(void);
int revertToStart(void);
int update(void);
//public methods to obtain stiffness, mass, damping, and residual

information
const Matrix &getTangentStiff(void);
const Matrix &getInitialStiff(void);
const Matrix &getDamp(void);
const Matrix &getMass(void);
void zeroLoad(void);
int addLoad(ElementalLoad *theLoad, double loadFactor);
int addInertiaLoadToUnbalance(const Vector &accel);
const Vector &getResistingForce(void);
const Vector &getResistingForceIncInertia(void);
//public methods for output

int sendSelf(int commitTag, Channel &theChannel);
int recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker
&theBroker);
int displaySelf(Renderer &theViewer, int displayMode, float fact);
void Print(OPS_Stream &s, int flag=0);
Response *setResponse(const char **argv, int argc, OPS_Stream &s);
int getResponse(int responseID, Information &eleInformation);
//protected:
private:
115
//private member function - only availabe to objects of the class
double computeCurrentStrain(void) const;
//private attributes - a copy for each object of the class

UniaxialMaterial *theMaterial; //pointer to a material
ID externalNodes; // contains the id's of end nodes
Matrix trans; //hold the transformation matrix
// Vector *theLoad; // pointer to the load vector P
double L; //length of Pipe3 based on undeformed configuration

double C_3;
double A;
double D_C;
double d_y_class;
int eletag;
double Gamma; //weight per unit volume
Node *end1Ptr, *end2Ptr; //two pointer to the trusses nodes
Node *theNodes[2]; //two pointer to the trusses nodes
in a matrix form (AV)
//private class attribute

static Matrix trussK;
static Matrix trussD;
static Matrix trussM;
static Vector trussR;
};
#endif
5.7.4.2Class Definition
/* ****************************************************************** **
** **
** **
** **
** **
** Developed by: **
** **
** ****************************************************************** */
// $Date: 2008/07/18 18:05:53 $
// $Source: /usr/local/cvs/OpenSees/SRC/element/Pipe/Pipe3.cpp,v $
116
// Created: 07/08
// Revision: A
//
// Description: This file contains the implementation for the Pipe3 class.
// The current implementation does not work very well for Newton and
EnergyIncr
// (Flow does not go to zero at end of analysis)
// It works well for DispIncr and KrylovNewton Analysis eventhough some
leakage
// seem to be happening after the end of the analysis. More verification is
needed
// in order to see the effect of the d_c parameter and incrementation and
integration
// schemes. Also, minor controls might be needed to calculate R at very low
i.
//
#include "Pipe3.h"
#include <Information.h>
#include <Parameter.h>
#include <Domain.h>
#include <Node.h>
#include <Channel.h>
#include <FEM_ObjectBroker.h>
#include <UniaxialMaterial.h>
#include <Renderer.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <ElementResponse.h>
#include <Matrix.h>
#include <Vector.h>
#include <ElasticMaterial.h>
// initial the class wide variables

Matrix Pipe3::trussK(6,6);
Matrix Pipe3::trussM(6,6);
Matrix Pipe3::trussD(6,6);
Vector Pipe3::trussR(6);
Pipe3::Pipe3(int tag,
int Nd1, int Nd2,
UniaxialMaterial &theMat,
double a, double c3, double g, double dc)
:Element(tag,ELE_TAG_Pipe3),
theMaterial(0),
externalNodes(2),/*theLoad(0),*/
trans(1,4),L(0.0), A(a), C_3(c3), Gamma(g), D_C(dc),
end1Ptr(0), end2Ptr(0), eletag(tag),
d_y_class(0.0)
117
{
//create a copy of the material object
theMaterial=theMat.getCopy();
//fill in the ID containing external node info with node id's

for (int i=0; i<2; i++)
theNodes[i] = 0;
}
//constructor which should be invoked by an FE_ObjectBroker only

Pipe3::Pipe3()
:Element(0,ELE_TAG_Pipe3),
theMaterial(0),
externalNodes(2),/*theLoad(0),*/
trans(1,4), L(0.0), A(0.0), C_3(0.0), Gamma(0.0),
D_C(0.0),end1Ptr(0), end2Ptr(0),
d_y_class(0.0)
{
for (int i=0; i<2; i++)
theNodes[i] = 0;
}
Pipe3::~Pipe3()
{
if (theMaterial !=0)
delete theMaterial;
}
int Pipe3::getNumExternalNodes(void) const

{
return 2;
}
const ID &
Pipe3::getExternalNodes(void)
{
return externalNodes;
}
int
Pipe3::getNumDOF(void){
return 6;
}
Node **
Pipe3::getNodePtrs(void)
{
return theNodes;
}
void
Pipe3::setDomain(Domain *theDomain)
{
118
//first ensure nodes exist in Domain and set the node
pointers
if (theNodes[0]==0)
return;
if (theNodes[1]==0)
return;
// call the DomainComponent class method

this->DomainComponent::setDomain(theDomain);
//ensure connected nodes have corrent number of dof's

if ((dofNd1 !=3) || (dofNd2 !=3))
return; //don't go any further otherwise segmentation
fault
//now determine the length & transformation matrix

double dx= end2Crd(0)-end1Crd(0);

double dy= end2Crd(1)-end1Crd(1);
d_y_class=dy;
L=sqrt(dx*dx+dy*dy);
if (L==0.0)
return;
double cs=dx/L;
double sn=dy/L;
trans(0,0)=-cs;
trans(0,1)=-sn;
trans(0,2)= cs;
trans(0,3)= sn;
// // determine the nodal mass for lumped mass approach

// M=M*A*L/2; //M was set to rho by the constructor
}
int
Pipe3::commitState()
{
return theMaterial->commitState();
}
int
Pipe3::revertToLastCommit()
119
{
return theMaterial->revertToLastCommit();
}
int
Pipe3::revertToStart()
{
return theMaterial->revertToStart();
}
int
Pipe3::update()
{
//determine the current strain given trial displacements at
nodes
double strain=this->computeCurrentStrain();
//set the strain in the materials

theMaterial->setTrialStrain(strain);
return 0;
}
const Matrix &

Pipe3::getTangentStiff(void)
{
if (L==0) {//if length ==zero - we zero and return
trussK.Zero();
return trussK;
}
//get the current E from the material for the strain that
was set
// at the material when the update() method was invoked
double E = theMaterial->getTangent();

Matrix K_temp(4,4);
K_temp *=A*E/L;
// trussK.Zero();
trussK(2,2)=0.0;
trussK(5,5)=0.0;
120
trussK(5,2)=0.0;
trussK(2,5)=0.0;
return trussK;
}
const Matrix &
Pipe3::getInitialStiff(void)
{
if (L==0) {
trussK.Zero();
return trussK;
}
//get the current strain from the material

//get the current stress from the material

double stress = theMaterial->getStress();
//compute the tangent

double E=stress/strain;

Matrix K_temp(4,4);
K_temp *=A*E/L;
// trussK.Zero();
trussK(2,2)=0.0;
trussK(5,5)=0.0;
trussK(5,2)=0.0;
trussK(2,5)=0.0;
return trussK;
}
const Matrix &

Pipe3::getDamp(void)
{
//No damping associated with this type of element
trussD.Zero();
// Darcy-Weisbach components
const Vector &vel1 = end1Ptr->getTrialVel();
const Vector &vel2 = end2Ptr->getTrialVel();
// const Vector &disp1 = theNodes[0]->getIncrDisp();
121
// const Vector &disp2 = theNodes[1]->getIncrDisp();
double i_a;
i_a=vel2(2)/L-vel1(2)/L+d_y_class*Gamma/L; //Element
hydraulic gradient
// Find the sign of the gradient:
double sign_i_a=1;
if (i_a<0)
sign_i_a=-1;
//Control for singularities in the incrementaly linearized equation (can be
changed)
if (sign_i_a*i_a<D_C)
i_a=sign_i_a*D_C;
// i_a=D_C;
trussD(2,2)=-C_3/2/sqrt(sign_i_a*i_a)/L;
trussD(5,5)=-C_3/2/sqrt(sign_i_a*i_a)/L;
trussD(5,2)=C_3/2/sqrt(sign_i_a*i_a)/L;
trussD(2,5)=C_3/2/sqrt(sign_i_a*i_a)/L;
return trussD;
}
const Matrix &

Pipe3::getMass(void)
{
if (L==0){
trussM.Zero();
return trussM;
}
// At this point we have zero lumped mass

trussM.Zero();
return trussM;
}
void
Pipe3::zeroLoad(void)
{
//does nothing - no element load associated with this
object
}
int
Pipe3::addLoad(ElementalLoad *theLoad, double loadFactor)
{
opserr <<"MyTruss::addLoad - load type unknown for truss with
tag: " << this->getTag() << endln;
return -1;
}
int
Pipe3::addInertiaLoadToUnbalance(const Vector &accel)
{
return 0;
}
122
const Vector &
Pipe3::getResistingForce()
{
if (L==0) {//if length ==zero - zero and return
trussR.Zero();
return trussR;
}
// R=Ku-Pext
//force =F*transformation
double force = A* theMaterial->getStress();

if ((vel2(2)-vel1(2)+d_y_class*Gamma)>0.0)
{
trussR(2)=C_3*sqrt((vel2(2)-
vel1(2)+d_y_class*Gamma)/L);
trussR(5)= -C_3*sqrt((vel2(2)-
double deleteme4=(vel2(2)-vel1(2)+d_y_class*Gamma);
double deleteme5=vel2(2);
}
else
{
trussR(2)=-C_3*sqrt((-vel2(2)+vel1(2)-
d_y_class*Gamma)/L);
trussR(5)= +C_3*sqrt((-vel2(2)+vel1(2)-
double deleteme4=(-vel2(2)+vel1(2)-d_y_class*Gamma);
}
double deleteme2=trussR(2);
double deleteme3=trussR(5);
return trussR;
}
const Vector &

Pipe3::getResistingForceIncInertia()
{
this->getResistingForce();
//No inertia is included in the in this element formulation

return trussR;
}
int
123
Pipe3::sendSelf (int commitTag, Channel &theChannel)
{
int dataTag=this->getDbTag();
// Pipe3 packs it's data into a Vector and sends this to

theChannel
//along with it's dbTag and the commitTag passed in the
arguments
Vector data(7);
data(0)= this->getTag();
data(1)=A;
data(4)=C_3;
data(5)=Gamma;
data(6)=D_C;
data(2)=theMaterial->getClassTag();
int matDbTag=theMaterial->getDbTag();
if (matDbTag==0) {
matDbTag =theChannel.getDbTag();
if (matDbTag !=0)
theMaterial->setDbTag(matDbTag);
}
data(3)=matDbTag;
theChannel.sendVector (dataTag, commitTag, data);
theChannel.sendID(dataTag, commitTag, externalNodes);
theMaterial->sendSelf(commitTag, theChannel);
return 0;
}
int
Pipe3::recvSelf(int commitTag, Channel &theChannel,
FEM_ObjectBroker &theBroker)
{
int dataTag= this->getDbTag();
Vector data(7);
theChannel.recvVector(dataTag, commitTag, data);
this->setTag((int)data(0));
A=data(1);
C_3=data(4);
Gamma=data(5);
D_C=data(6);
theChannel.recvID(dataTag, commitTag, externalNodes);
int matClass=data(2);
int matDb = data(3);
theMaterial= theBroker.getNewUniaxialMaterial(matClass);
theMaterial->setDbTag(matDb);
theMaterial->recvSelf(commitTag, theChannel, theBroker);
return 0;
}
124
int
Pipe3::displaySelf(Renderer &theViewer, int displayMode,
float fact)
{
Vector v1(3);
Vector v2(3);
for (int i=0; i<2;i++) {
}
if (displayMode==3) {
//use the strain as the drawing measure
return theViewer.drawLine(v1, v2, strain, strain);
}
else if (displayMode==2){
//otherwise use the material stress
double stress =A*theMaterial->getStress();
return theViewer.drawLine(v1,v2,stress, stress);
}
else{
//use the axial force
double force = A*theMaterial->getStress();
return theViewer.drawLine(v1,v2,force,force);
}
}
void
Pipe3::Print(OPS_Stream &s, int flag)
{
//compute the strain and axial force in the member
double strain, force;
if (L==0.0) {
strain=0;
force=0.0;
}
else{
strain = theMaterial->getStrain();
force=A*theMaterial->getStress();
}

if ((vel2(2)-vel1(2)+d_y_class*Gamma)>0.0)
{
trussR(2)=C_3*sqrt((vel2(2)-
125
trussR(5)= -C_3*sqrt((vel2(2)-
double deleteme4=(vel2(2)-vel1(2)+d_y_class*Gamma);
}
else
{
trussR(2)=-C_3*sqrt((-vel2(2)+vel1(2)-
trussR(5)= +C_3*sqrt((-vel2(2)+vel1(2)-
}
if (flag==0) {//print everythin

s<< "Element: " <<this->getTag();
s<< " type: My Truss iNode: "<< externalNodes(0);
s<< " jNode: "<<externalNodes(1);
s<< " Area: "<< A;
if (Gamma!=0) s << "Gamma: "<<Gamma;
s<< " \n\t strain: " <<strain;

s<< " axial load: " <<force;
s<< " \n\t unbalanced load: " <<trussR;
s<< " \t Material: " << *theMaterial;
s<< endln;
} else if (flag==1) {//just print ele id, strain and force
s<< this->getTag() << " " <<strain << " " << force
<<endln;
}
}
Response *
Pipe3::setResponse(const char **argv, int argc , OPS_Stream
&s)
{
// we compare arg(0) for known response types for the Truss
//axial force
if(strcmp(argv[0], "axialForce")==0)
return new ElementResponse(this, 1, 0.0);
//a material quantity

else if (strcmp(argv[0], "material")==0)
return theMaterial->setResponse(&argv[1], argc-1, s);
else
return 0;
}
int
Pipe3::getResponse(int responseID, Information
&eleInformation)
{
switch (responseID){
126
case -1:
return -1;
case 1:
return eleInformation.setDouble (A*theMaterial-
>getStress());
default:
return 0;
}
}
double
Pipe3::computeCurrentStrain(void) const
{
//determine the strain
double dLength=0.0;
for (int i=0;i<2;i++)
dLength -= (disp2(i)-disp1(i))*trans(0,i);
double strain =dLength/L;
return strain;
}
5.7.4.3Tcl/tk command interpreter

/* ****************************************************************** **
** **
** **
** **
** **
** Developed by: **
** **
** ****************************************************************** */
// $Revision: 1. $
// $Date: 2008/07/20 19:20:46 $
// $Source: /usr/local/cvs/OpenSees/SRC/element/pipe/TclPipe3Command.cpp,v $
// File: ~/element/TclPipe3Command.C
127
//
// Written: avytin
// Created: 09/08
// Revision: A
//
// Description: This file contains the implementation of the
TclModelBuilder_Pipe3()
// command.
//
// What: "@(#) TclModelBuilder.C, revA"
#include <stdlib.h>
#include <string.h>
#include <Domain.h>
#include "Pipe3.h"
#include <TrussSection.h>
#include <TclModelBuilder.h>
#include <CorotTruss.h>
#include <CorotTrussSection.h>
extern void printCommand(int argc, TCL_Char **argv);
int
TclModelBuilder_Pipe3(ClientData clientData, Tcl_Interp *interp, int argc,
TCL_Char **argv, Domain*theTclDomain,
TclModelBuilder *theTclBuilder,
int eleArgStart){
//make sure at least one other
argument to contain type of system
if (argc!=10){
interp->result = "WARNING bad
command - Pipe3 eleId iNode jNode matID Area c_3 Gamma d_c";
return TCL_ERROR;
}
//get the id, x_loc, y_loc

int trussId, iNode, jNode, matID;
double A, C_3, Gamma, D_C;
if (Tcl_GetInt(interp,argv[2],
&trussId)!= TCL_OK){
invalid eleId - Pipe3 eleId iNode jNode matID Area c_3 Gamma d_c";
return TCL_ERROR;
}
&iNode) != TCL_OK) {
invalid iNode - Pipe3 eleId iNode jNode matID Area c_3 Gamma d_c";
return TCL_ERROR;
}
&jNode) != TCL_OK) {
invalid jNode - Pipe3 eleId iNode jNode matID Area c_3 Gamma d_c";
128
return TCL_ERROR;
}
&matID) != TCL_OK) {
invalid matID - Pipe3 eleId iNode jNode matID Area c_3 Gamma d_c";
return TCL_ERROR;
}
&A) != TCL_OK) {
invalid Area - Pipe3 eleId iNode jNode matID Area c_3 Gamma d_c";
return TCL_ERROR;
}
&C_3) != TCL_OK) {
invalid C_3 - Pipe3 eleId iNode jNode matID Area c_3 Gamma d_c";
return TCL_ERROR;
}
&Gamma) != TCL_OK) {
invalid Gamma - Pipe3 eleId iNode jNode matID Area c_3 Gamma d_c";
return TCL_ERROR;
}
&D_C) != TCL_OK) {
invalid d_c - Pipe3 eleId iNode jNode matID Area c_3 Gamma d_c";
return TCL_ERROR;
}
UniaxialMaterial *theMaterial =
theTclBuilder->getUniaxialMaterial(matID);
if (theMaterial ==0) {
opserr << "WARNING TclPipe3 -
Pipe3 - no Material found with tag ";
opserr << matID << endln;
return TCL_ERROR;
}
//now create the truss and add it

to the domain
// MyTruss *theTruss = new
MyTruss(trussId,iNode,jNode,*theMaterial,A,M);
Element *theTruss = 0;
theTruss=new
Pipe3(trussId,iNode,jNode,*theMaterial,A,C_3,Gamma, D_C);
if (theTruss==0) {
129
Pipe3 - ran out of memory for node ";
return TCL_ERROR;
}
if (theTclDomain-
>addElement(theTruss)==false) {
delete theTruss;
Pipe3 - could not add Pipe3 to the domain";
return TCL_ERROR;
}
//Everything is OK
return TCL_OK;
}
130
5.7.4.4
6 FEM Simulation
6.1 Finite element model description
6.1.1 Model Setup
At first the modeling space for the numerical analyses is selected. The Yolo Loam layer, due to
its very low permeability is modeled with 2DOF quad elements, and total stress analysis is em-
ployed. The elasto-plastic pressureindependmultiyield model is used in this case. The Sand lay-
ers are modeled with 3DOF quad-up elements, and effective stress analysis is used, since we
care about the coupled pore-pressure displacement phenomena. The pressuredependmul-
tiyield02 constitutive model is used to model the mechanical deformation of sand. These two
types of domains (2DOF and 3DOF) are connected to each other with equalDOF objects. The
acceleration is applied as a uniform excitation in the bottom of the model. Finally, the drains
have been treated either as perfect drains (by means of fixed pore pressure boundary condi-
tions), or as finite transmissivity drains (by means of the drain elements presented in the pre-
vious chapter).
131
Figure 6-1 Different domains employed in the FE analysis
The laminar box in the FE model should be modeled. The laminar box fixes the displacements of
one side of the model to the other side. The connection of the soil elements with the laminar
box is modeled as fixity for the sand layer and as a tension-only connection at the layer of Yolo
Loam. This effect is important on the first soil layer (Yolo loam) because the compression forces
are small compared to the rest of the model. Figure 6-3 shows exactly the setup being used to
model the connection of the elements of the model to the laminar box. On the other hand, the
top layer of Yolo Loam has very small effective stresses, so it is bound to move away from the
laminar box. The grid nodes of this layer are connected with a non-linear no-tension zerolength
element to a set of nodes that are on the laminar box. Next the nodes on the laminar box are
connected to the ones on the other side with equalDOF connections in both the x and y direc-
tion. Finally, nodal mass has been added to the nodes on the sides to include the effect of the
mass of the laminar box.
132
Figure 6-2 FE Model Setup
133
Figure 6-3 Boundary conditions of the FE model
6.1.2 Selection of model parameters
6.1.2.1Constitutive soil parameters

As mentioned in the data limitations section of the experiment presentation, the relative densi-
ty of the dense sand is assumed to be 80% and the density of the loose sand is assumed to be
40%. From the Opensees manual (Mazzoni, McKenna, & Fenves, 2005) we then select the
properties for the pressuredependmultiyield02 model. The parameters are shown in Table 6-1
toTable 6-2.
Table 6-1 Dense Sand model properties
PressureDependendMultiyield02 (for a dense sand, Dr=80%)
134
Ρ 2.07 cpressuredependence 0.5
Gref 130000 ψphasetransformation 26.0
Kref 260000 c1 0.013
Φ 36.5 c3 0.0
γpeak 0.1 d1 0.3
pref 80 d3 0.0
Table 6-2 Medium-loose sand model properties
PressureDependendMultiyield02 (for a medium-loose sand, Dr=40%)

Ρ 1.98 cpressure dependence 0.5
Gref 90000 ψphase transformation 26.0
Kref 220000 c1 0.067
Φ 32.0 c3 0.23
γpeak 0.1 d1 0.06
pref 80 d3 0.27
Density of the Yolo loam was measured initially to be 13kN/m3 whereas from water content
tests after the end of shaking it was calculated as 18.5kN/m3. This is due to consolidation hap-
pening during the test. The width of the Yolo loam layer changes during the test. For our ana-
lyses we keep 13kN/m3 by using the initial 1m layer thickness.
The undrained strength of clays is dependent on the vertical effective stress. Its value is mainly
important in order to estimate the slippage between the Sand and the Yolo Loam Layers. In the
setup of this experiment the vertical effective stresses are from 0kPa to 13kPa.
6.1
135
6.2
So, according to Ladd for normally consolidated clays, for a DSS mode of shearing at the bottom
of the layer:
6.3
6.4
Which gives us an indication of the value of the undrained strength of this layer.
On the other hand there are indirect measurements of the Yolo Loam cohesion. During shaking
one has measurements exactly below the Yolo Loam Sand interface and above the Yolo Loam
Sand interface. From Newton’s theory:
6.5
6.6
6.7
So if we plot the relative acceleration between the two layers (Figure 6-4), we find that
6.8
Which means that:
6.9
We select for our analysis an approximate value of:
6.10
Due to the uncertainties in measured quantities.
136
Figure 6-4 Relative acceleration below and above the Loose Sand - Yolo Loam interface for the SSK01-10 phase of shaking
Table 6-3 Yolo Loam model Properties
PressureIndependMultiYield (Soft Clay)
Ρ 1.3
Gref 13000
Kref 65000
C 6.0
γpeak 0.1
6.1.2.2Soil Permeability
The permeability of the soil is measured to be in the model scale:
6.11
So the permeability is at the prototype scale in order to match diffusion N (15) times larger than
the one in model scale:
6.12
137
At the side where the drains are installed, permeability should be scaled according to Hird, so
that the average degree of dissipation of excess pore pressure in the plane strain model is the
same as the average degree of dissipation of excess pore pressure in the real 3-D scenario.
6.13
6.14
6.1.2.3Laminar box parameters

The approximate weight of the FSB3 laminar box is approximately the same to the FSB2 laminar
box properties shown in Table 6-4. This is evenly assigned to the side nodes of the drain as
0.0105Mgr/side node.
Table 6-4 Design details of container FSB2 for the large centrifuge
6.1.2.4PV-Drain parameters
Ths PV-Drain parameters that should be calculated are the transmissivity for laminar flow, the
transmissivity for fully turbulent flow, the stiffness for truss behavior, and the storage capacity.
The inner diameter of the drain in model scale is 7mm and the outer diameter is 9mm. So the
effective area that influences the mechanical truss behavior is:
138
6.15
6.16
6.17
6.18
As far as the storage capacity is affected the volume that needs to fill up is a cylinder with di-
ameter the inner diameter of the drain. Thus the cross section of the drain needed to estimate
the effect of storage capacity is:
6.19
6.20
6.21
6.22
The material for the drains is nylon with an elastic modulus of:
6.23
Now we need to estimate the transmissivity parameters. We estimate Cl initially for the model
scale:
6.24
6.25
6.26
139
6.27
6.28
This is the transmissivity parameter we should use if we were solving the true 3D radial drai-
nage problem in the prototype scale. In reality we solve a plane strain problem, so we need to
scale the transmissivity according to Hird:
6.29
6.30
6.31
which is the parameter we are using.
For the fully turbulent drains we introduce a similar type of analysis. We estimate initially Ct for
the model scale:
6.32
6.33
6.34
6.35
140
6.36
We now scale the transmissivity according to Hird:
6.37
6.38
6.39
which is the parameter we are using in our analysis.
A summary of results is presented in Table 6-5.

Table 6-5 Drain properties in prototype scale
Mechanical Deformation
Storage Capacity
Drain properties
6.2 Base case analysis
This part shows the results of the best approximation to the centrifuge experiment. This analy-
sis assumes that flow in the drains is laminar, and includes the effect of storage capacity. In this
section we compare analytically the simulated results with the experimental results.
141
At first we examine the applicability of the assumption of laminar flow inside the drains. In Fig-
ure 6-5 we examine on the first plot, the flow coming out of the drain. The green line
represents the limit above which the flow becomes turbulent for the prototype scale and the
yellow line represents the limit above which the flow becomes laminar for the model scale. One
can see that most of the flow is taking place in the laminar regime in the model scale, and in the
turbulent regime for the prototype scale. Also, in this Figure we estimate vertical displacement
on the top of the drain by integrating the flow of water coming out of the pipe. As we can see
this indirect method of predicting the displacement gives larger displacements from the direct
calculation of the displacement, due to the flow taking place from the un-treated side to the
treated one.
6.2.1 Predicted Excess Pore Pressures
In Figure 6-6 to Figure 6-9 the predicted pore pressures are compared with the ones experi-
mentally measured. By looking at the plots we can clearly see that the prediction captures well
the measured values. The effect of the drains in decreasing the accumulation of excess pore
pressure during shaking is evident.
Discrepancies though exist between the simulation and the experiment which are attributed to
two main reasons: use of larger permeability and relative density parameters than the ones
that exist in reality.
From points E and F we get two conclusions:
142
• By looking the initial development of excess pore pressure, the initial accumulation is
not captured very well. This is attributed possibly to the fact that the estimations for the
actual density of the soil vary from 30% to 40%. Since the selected value for this analysis
is 40% it is possible that the actual relative density is much smaller, thus faster accumu-
lation of excess pore pressure is expected.
• By examining the dissipation phase of the excess pore pressure it is apparent that the
selected permeability of the sand layer was too high; excess pore pressures dissipate
much faster in the simulation
From point D:
• We fail to capture exactly the measured pore pressures. This is attributed to the fact
that low depth sensors sank significantly during the analysis
From points A, B, C:
• By examining the initial accumulation phase, we can see that accumulation is happening
much slower than in reality. This is attributed one more time to erroneous selection of
relative density and permeability parameters
• Dissipation is one more time happening faster in the simulation, possibly due to exces-
sively large permeability value selection.
6.2.2 Predicted accelerations
143
In Figure 6-9 Figure 6-11 a comparison of the predicted vs. the simulated accelerations is pre-
sented. In general, the analysis does not predict very well the acceleration on the top of the
layer. On the top of the sand layer we fail to get the de-amplification on the un-treated side and
the amplification on the treated side correctly. We are under predicting accelerations in the
treated zone, because the actual excess pore pressures are in reality larger than the ones in si-
mulation, so the actual response is that of a softer material. In the untreated side we are over-
predicting the accelerations again due to the fact that the material is so soft (due to liquefac-
tion) that shear waves cannot propagate. These phenomena are one more time attributed to
the erroneous relative density and permeability parameters.
6.2.3 Predicted displacements
In Figure 6-12Figure 6-14 we compare the simulated to the experimental horizontal displace-
ments. The simulated displacements are smaller than the experimental ones. Figure 6-18 Figure
6-20 show the horizontal displacement profiles. From these, one can see that the small strain
assumption used in the simulation is valid for the simulation of this test. In Figure 6-21 we see
the final predicted deformed shape.
In Figure 6-15Figure 6-17 a comparison between the experimental and simulated vertical dis-
placements is presented. The predicted vertical displacements are similar to the measured
ones. Should we look though to Figure 6-17 we see that close to the laminar box we under pre-
dict settlements, whereas closer to the middle we over predict settlements. Some authors
(Arulanandan & Subico, 1992) suggest that the value of permeability is non-constant during the
144
shaking, and that this strongly affects measured vertical displacements in the centrifuge expe-
riment.
On the contrary, the author believes this is due to the non-uniformity of the gravitational field,
a hypothesis that needs though to be evaluated in more detail. According to this, in the case
that the surface of the soil was horizontal before the experiment, after significant shaking, the
soil surface will become a curve with radius equal to the rotation radius of the centrifuge. This
explains why the experiment compared with the simulation –which assumers uniform vertical
gravitational field - produces larger displacements at the sides and smaller closer to the middle.
145
-3
x 10 Flow coming out of drain No2 vs Time
4
Flow
Flow (m /s)
2 Laminar flow limit (prototype)
3
Laminar flow limit (model)
0
-2
0 2 4 6 8 10 12 14
Time (s)
-3 Volume of water coming out of drain No2 vs Time
x 10
6
Volume (m )
3
0
0 2 4 6 8 10 12 14
Time (s)
-3 Vertical displacement on top of drain No2
x 10
4
Displacement (m)
Indirect
Direct
2
-2
0 2 4 6 8 10 12 14
Time (s)
Figure 6-5 Flow and volume of water coming out of drain No2, and comparison of directly and indirectly predicted displacements on the top of the drain during the first shak-
2
ing event (amax=0.687m/s )
146
Treated side Untreated side
20
A D
10
0
60
Pore Pressure, p (kPa)
B E
50
40
30
80
C F
70
60
50
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Time, t (s) Time, t (s)
Experiment
Simulation
2
Figure 6-6 Comparison of measured and simulated pore pressures for the first phase of shaking (amax=0.687m/s )
147
20 A D
10
0
60
B E
50
40
30
80 C F
70
60
50
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Time, t (s) Time, t (s) Experiment
Simulation
2
Figure 6-7 Comparison of measured and simulated pore pressures for the second phase of shaking (amax=1.071m/s )
148
20 A D
10
0
60 B E
50
40
30
90 C F
80
70
60
50
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Time, t (s) Time, t (s) Experiment
Simulation
2
Figure 6-8 Comparison of measured and simulated pore pressures for the third phase of shaking (amax=2.943m/s )
149
2
1 A D
-1
Horizontal Acceleration, α (m/s )
2
-2
2
1 B E
-1
-2
1
C F
-1
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Experiment
Simulation
2
Figure 6-9 Comparison of measured and simulated horizontal accelerations for the first phase of shaking (amax=0.687m/s )
150
2
A D
-2
2
2
B E
-2
2
C F
-2
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Experiment
Simulation
2
Figure 6-10 Comparison of measured and simulated horizontal accelerations for the second phase of shaking (amax=1.071m/s )
151
10
A D
5
0
2
-5
5 B E
-5
5
C F
-5
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Experiment
Simulation
2
Figure 6-11 Comparison of measured and simulated horizontal accelerations for the third phase of shaking (amax=2.943m/s )
152
0.2
0.1 C F
-0.1
-0.2
0.2
Horizontal, u (m)
0.1 B E
-0.1
-0.2
0.1
A D
-0.1
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Experiment
Simulation
2
Figure 6-12 Comparison of measured and simulated horizontal displacements for the first phase of shaking (amax=0.687m/s )
153
0.3 C F
0.2
0.1
0
Horizontal Displacement, u x (m)
-0.1
0.3 B E
0.2
0.1
-0.1
0.2 A D
0.1
-0.1
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Experiment
Simulation
2
Figure 6-13 Comparison of measured and simulated horizontal displacements for the second phase of shaking (amax=1.071m/s )
154
0.3 C F
0.2
0.1
0
-0.1
0.3 B E
0.2
0.1
-0.1
0.2 A D
0.1
-0.1
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Experiment
Simulation
2
Figure 6-14 Comparison of measured and simulated horizontal displacements for the third phase of shaking (amax=2.943m/s )
155
0.04
C F
0.02
0
Vertical Settlement, uy (m)
0.08
B E
0.06
0.04
0.02
0
0.1 A D
0.05
0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Experiment
Simulation
2
Figure 6-15 Comparison of measured and simulated settlements for the first phase of shaking (amax=0.687m/s )
156
0.06
C F
0.04
0.02
0
B E
0.1
0.05
0.1 A D
0.05
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Experiment
Simulation
2
Figure 6-16 Comparison of measured and simulated settlements for the second phase of shaking (amax=1.071m/s )
157
0.15
0.1 C F
0.05
0
0.25 B E
0.2
0.15
0.1
0.05
0
0.2 A D
0.15
0.1
0.05
0
-0.05
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Experiment
Simulation
2
Figure 6-17 Comparison of measured and simulated settlements for the third phase of shaking (amax=2.943m/s )
158
X-Acceleration (m/s )
2
Applied X-Acceleration on the bottom Sections

1
-1
0 5 10 15
Time (s)
Section A Section B
6 6
5 5
4 4
z-coordinate (m)
z-coordinate (m)
3 3
2 2
1 1
0 0
-0.1 0 0.1 -0.1 0 0.1
X-Displacement (m) X-Displacement (m)
2
Figure 6-18 Horizontal displacements profiles during the first phase of shaking (amax=0.687m/s )
159
2

2
-2
120 125 130 135
Time (s)
Section A Section B
6 6
5 5
4 4
z-coordinate (m)
z-coordinate (m)
3 3
2 2
1 1
0 0
-0.2 0 0.2 -0.2 0 0.2
2
Figure 6-19 Horizontal displacements profiles during the second phase of shaking (amax=1.071m/s )
160
2

5
-5
235 240 245 250 255
Time (s)
Section A Section B
6 6
5 5
4 4
z-coordinate (m)
z-coordinate (m)
3 3
2 2
1 1
0 0
-0.5 0 0.5 -0.5 0 0.5
2
Figure 6-20 Horizontal displacements profiles during the third phase of shaking (amax=2.943m/s )
161
Figure 6-21 Final predicted deformed shape
162
6.3 Effect of different approximations in PV-drains simulations
6.3.1 Drain resistance
In this section we compare two assumptions: (a) perfect drains, and (b) laminar drains. Both
models do not include the effect of storage capacity. We can see from Figure 6-22 that in the
current model there is no difference between the laminar drain and perfect drain assumption.
This is attributed to the fact that the drains are very permeable and the permeability of the soil
is high (due to scaling).
163
2
60 Third Phase (amax=2.943m/s )
40
20
Pore Pressure, p(kPa)
2
60 Second Phase (amax=1.071m/s )
40
20
50 2
First Phase (amax=0.687m/s )
40
30
20
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Perfect Drains
Laminar Drains
Figure 6-22 Comparison of predicted pore pressures using the perfect drains vs. laminar drains assumption
164
6.3.2 Drain stiffness
In this section we compare an analysis with laminar drains including the drain stiffness to one
ignoring the drain stiffness. In Figure 6-23 a comparison of the predicted pore pressures, it is
shown that for the second phase of shaking the predicted pore pressures are higher when the
drain stifness effect is included, whereas at the third phase they are lower. This is attributed to
two competing effects: the drains do not allow the soil to settle, reducing the excess pore
pressure, whereas they impede the soil from undergoing large shear deformations, after the
dilation angle, incresing the excess pore pressure. These effects do not reflect very clearly on
the predicted accelerations (Figure 6-24), but are better illustrated in the predicted horizontal
displacements shown in Figure 6-25.
165
166
2
40
20
2
40
20
50 2
40
30
20
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Incl. Drain Stiffness
Ign. Drain Stifness
Figure 6-23 Comparison of predicted pore pressures including and ignoring the effect of drain stiffness
167
5
2
Third Phase (amax=2.943m/s )
-5
Horizontal Acceleration, ux (m/s )
2
-10
4
2
Second Phase (amax=1.071m/s )
2
0
-2
-4
-6
2 2
0
-2
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Ign. Drain Stiffness
Figure 6-24 Comparison of predicted horizontal accelerations including and ignoring the effect of drain stiffness
168
0.3
2
0.2
0.1
0
-0.1
0.3
2
0.2
0.1
-0.1
2
0.1
-0.1
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Ign. Drain Stiffness
Figure 6-25 Comparison of predicted horizontal displacements including and ignoring the effect of drain stiffness
169
6.3.3 Drain storage capacity
In this section we compare simulations with drains, with and without storage capacity. A com-
parison using laminar drains is presented in Figure 6-26 to Figure 6-28, where the pore pres-
sures, the horizontal accelerations, and the horizontal displacements are plotted. Taking into
account the effect of storage capacity, increases the excess pore pressures, and the predicted
permanent displacements on top of the treated side.
170
171
2
40
20
2
40
20
50 2
40
30
20
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Ignor. Drain Storage
Incl. Drain Storage
Figure 6-26 Comparison of predicted pore pressures illustrating the effect of drain storage
172
2
-5
Horizontal Acceleration, ax (m/s )
2
4
2
2
0
-2
-4
-6
2 2
0
-2
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Ign. Drain Storage
Incl. Drain Storage
Figure 6-27 Comparison of predicted horizontal accelerations illustrating the effect of drain storage
173
0.3
2
0.2
0.1
0
-0.1
0.3
2
0.2
0.1
-0.1
2
0.1
-0.1
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Ignor. Drain Storage
Incl. Drain Storage
Figure 6-28 Comparison of predicted horizontal displacements illustrating the effect of drain storage
174
6.3.4 Drain turbulence
A comparison of the predicted pore pressures, under the assumption that the flow in the drains
is fully turbulent vs. the assumption that the flow in the drains is laminar, is plotted in Figure
6-29. Due to the significant transmissivity of the drains and the large permeability of the soil
used in the centrifuge model, we cannot see discrepancies between the results using laminar
drains and fully turbulent drains.
175
176
2
40
20
2
40
20
50 2
40
30
20
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Laminar Drains
Fully Turbulent Drains
Figure 6-29 Comparison of predicted pore pressures using laminar vs. fully turbulent drains
177
7 Summary – Conclusions
This thesis focuses on three different issues: (a) it establishes a method to estimate the re-
sponse of soil with vertical drains,(b) it discusses similitude issues for centrifuge modeling,(c) it
performs validation analyses of the implemented numerical methods,(d) and addresses design
recommendations for researchers and engineers.
7.1 Simulating Vertical Drains
This thesis used an uncoupled theory of mechanical deformation and flow inside a drain in or-
der to investigate the effect of earthquake drains during cyclic loading of sandy deposits. The
mechanical part of the drain’s response is idealized as a truss. The pore pressure flow part of
the drain’s response is assumed to be either fully turbulent or laminar.
Classes written in C++ have been implemented in the OpenSees FEM framework in order to si-
mulate the drains. Both formulations have been shown to work effectively, and give reasonable
results, within acceptable convergence levels and speed. For the fully turbulent drains, a consis-
tent Jacobian integration scheme is used due to its superior performance compared to the con-
tinuous Jacobian.
The effect of storage capacity has also been successfully modeled. Care should be taken by re-
searchers to achieve strict convergence when they use this feature, since small spurious per-
turbations in the predicted flow inside the drain could significantly affect the storage effect.
178
7.2 Similitude Issues
This thesis investigated similitude issues having to do with the design of the model scale drains.
It is shown that if both in the model scale and in the prototype scale the flow is laminar then
the transmissivity parameter Cl is N3 times larger in the prototype scale. It is also shown that if
both flow are turbulent, then the transmissivity parameter Ct is N2.5 times larger in the proto-
type scale.
Also, it is shown that the Re number in prototype scale and in model scale in the drains is dif-
ferent. It is actually possible that model scale flow is laminar whereas prototype scale flow is
turbulent. For this situation, a methodology is proposed to experimentalists, in order to choose
a model scale diameter for the drains (where flow is laminar) that fits the bet the prototype
scale diameter (where fully turbulent flow is expected). The prototype scale diameter DP of a
turbulent flow drain that corresponds to the model scale DM of a laminar flow drain is:
7.1
7.3 Validation
The implementation of the PV drains has been validated against centrifuge experiments per-
formed at UC-Davis. Results showed great accordance with the experiment. Discrepancies be-
tween the simulated and experimental results are attributed mostly to uncertainties of the soil
permeability and relative density. Validation illustrates that the storage effect is significant for
the numerical analysis of the cyclic response of soils treated with PV-drains. Also, with the de-
179
sign parameters used in the centrifuge model, it is found that the drains behave almost as per-
fect drains. A different set of tests is needed to evaluate the differences between perfect drain,
laminar drain, or fully turbulent drain assumption.
7.4 Future research
The following issues need to be addressed into future research:
• Implement the mechanical deformation of the drain using beam theory rather than
truss theory. This would allow use of the drain elements for simulation of the response
of layers improved with stone columns.
• Improve the constitutive soil model predictions. Most important is the ability to be able
to simulate the response of layers of the same sand at different stress levels, with dif-
ferent void ratios using the same set of parameters.
• Apart from acquiring expertise in the simulation of a geotechnical earthquake engineer-
ing problem consisting of shear wave propagation in soil layers with pre-installed PV
drains, it is very important issue to create design charts for engineers. These should give
recommendations for engineers for various levels of shaking (1), number of cycles (2),
frequency of loading (3), stress level (4), and soil type (5), in order to:
o evaluate the applicability of the improvement method by means of estimating

the maximum pore pressure ratio (umax/σv0)
o estimate the spacing and the types of the drains
o estimate the size of the zone that needs improvement
180
8 Bibliography
Andrianopoulos, K. I. (2006). Numerical Modeling of Static & Dynamic Behavior of Elastoplastic

Soils. Athens: Geotechnical Division NTUA.
Arulanandan, K., & Subico, J. J. (1992). Post liquefaction settlement of sands. Proceedings of the
Wroth Memorial Symposium. England: Oxford University.
Biot, M. (1956). Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. I.

Low-Frequency range. The Journal of the Acoustical Society of America , 168-178.
Conlee, C., Gallagher, P., Kamai, R., Boulanger, R., & Rix, G. (2008). Evaluation of the
Effectiveness of Colloidal Silica for Liquefaction Remediation: Centrifuge Data Report for CTC01.
Center for Geotechnical Modeling, University of California at Davis.
Corral, G. Cyclic Undrained Behavior of Loose Sand Treated by Colloidal Silica. Cambridge.
Dafalias, Y. F., & Manzani, M. T. (2004). Simple Plasticity Sand Model Accounting for Fabric
Change Effects. Journal of Engineering Mechanics , 622-634.
Dafalias, Y. (1994). Overview of constitutive models used in velacs. Verification of Numerical

Procedures for the Analysis of Soil Liquefaction Problems. Rotterdam: A.A. Balkema Publishers.
Elgamal, A., Yang, Z., & Parra, E. (2002). Computational Modeling of Cyclic Mobility and Post-
Liquefaction Site Response. Soil Dynamics and Earthquake Engineering , 259-271.
Gallagher, P. M., & Mitchell, J. (2002). Influence of Colloidal Silica Grout on Liquefaction
Potential and Cyclic Undrained Behavior of Loose Sand. Soil Dynamics and Earthquake
Engineering .
Gallagher, P. M., Pamuk, A., & Abdoun, T. (2007). Stabilization of Liquefiable Soils Using
Colloidal Silica Grout. Journal of Materials in Civil Engineering , 33-40.
Gallagher, P., & Mitchell, J. (2002). Influence of colloidal silica grout on liquefaction potential
and cyclic undrained behavior of loose sand. Soil Dynamics and Earthquake Engineering , 1017-
1026.
Garnier, J. (2001). Physical models in geotechnics:state of the art and recent advances. First
Columb Lecture .
181
Hird, C, C., Pyrah, I. C., & Russell, D. (1992). Finite Element Modelling of Vertical Drains beneath
Embankments on Soft Ground. Geotechnique , 499-511.
Hughes, T., & Pister, K. (1978). Consistent linearization in mechanics of solids and structures.
Computers and Structures 8 , 391-397.
Jeremic, B. (2008). Lecture Notes on Computational Geomechanics: Inelastic Finite Elements for
Pressure Sensitive Materials. Department of Civil and Environmental Engineering, University of
California, Davis.
Kamai, R., Howell, R., Conlee, C., Boulanger, R., Marinucci, A., Rathje, E., et al. (2008).
Evaluation of the Effectiveness of Prefabricated Vertical Drains for Liquefaction Remediation -
Centrifuge Data Report for RNK01. Department of Civil & Environmental Engineering, College of
Engineering, University of California at Davis.
Kamai, R., Kano, S., Conlee, C., Marinucci, A., Boulanger, R., Rathje, E., et al. (2008). Evaluation
of the Effectiveness of Prefabricated Vertical Drains for Liquefaction Remediation: Centrifuge
Data Report for SSK01. Center for Goetechnical Modeling, University of California at Davis.
Kammerer, A. M., Pestana, J. M., & Seed, R. (2002). Undrained response of monterey 0/30 sand
under multidirectional cyclic simple shear loading conditions. University of California, Berkeley.
Kodaka, T., Oka, F., Ohno, Y., Takyu, T., & Yamasaki, N. (2003). Modeling of Cyclic Deformation
and Strength Characteristics of Colloidal Silica Treated Sand. Geomechanics 2003 (pp. 205-216).
ASCE 2005.
Kolymbas, D. (2000). Introduction to Hypoplasticity. Taylor & Francis.
Lambe, T. W. (1973). Predictions in soil engineering. Geotechnique , 149-202.
Li, X. S., & Dafalias, Y. F. (2002). Constitutive Modeling of Inherently Anisotropic Sand Behavior.
Journal of Geotechnical and Geoenvironmental Engineering , 868-880.
Liao, H., Huang, C., & Chao, B. (2003). Liquefaction Resistance of a Colloidal Silica Grouted Sand.
Grouting 2003, (pp. 1305-1313).
Mazzoni, S., McKenna, F., & Fenves, G. L. (2005). Opensees Command Language Manual.
Nunez, I. L. (1988). Driving and tesion loading of piles in sand on a centrifuge. Centrifuge 1988,
(pp. 353-362). Balkema, Rotterdam.
182
Oka, F., Yashima, A., Tateishi, A., Taguchi, Y., & Yamashita, S. (1999). A cyclic elasto-plastic
model constitutive model for sand considering a plastic-strain dependence of the shear
modulus. Geotechnique , 661-680.
Onoue, A. (1988). Diagrams considering Well Resistance for Designing Spacing Ratio of Gravel
Drains. Soils and Foundations , 160-168.
Onoue, A., Mori, N., & Takano, J. (1987). In-situ experiment and analysis on well resistance of
gravel drains. Soils and Foundations , 42-60.
Papadimitriou, A. G., & Bouckovalas, G. D. (2001). Plasticity Model for Sand under Small and
Large Cyclic Strains: a Multi-axial Formulation . Soil Dynamics and Earthquake Engineering , 191-
204.
Papadimitriou, A. G., Bouckovalas, G. D., & Dafalias, Y. F. (2001). Plasticity Model for Sand under
Small and Large Cyclic Strains. Journal of Geotechnical and Geoenvironmental Engineering , 973-
983.
Papadimitriou, A., Moutsopoulou, M.-E., Bouckovalas, G., & Brennan, A. (2007). Numerical
Investigation of Liquefaction Mitigation using Gravel Drains. 4th International Conference on
Earthquake Geotechnical Engineering. Thessaloniki - GREECE 2007: 4ICEGE.
Pestana, J. M., & Whittle, A. J. (2002). Evaluation of a constitutive model for clays and sands:
Part I - sand behaviour. INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL
METHODS IN GEOMECHANICS , 1097-1121.
Pestana, J. M., & Whittle, A. J. (2002). Evaluation of a constitutive model for clays and sands:
Part II - clay behaviour. INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS
IN GEOMECHANICS , 1123-1146.
Pestana, J. M., & Whittle, A. J. (1999). Formulation of a unified constitutive model for clays and
sands. INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS IN
GEOMECHANICS , 1215-1243.
Pestana, J. M., Hunt, C. E., Goughnour, R. R., & Kammerer, A. M. Effect of Storage Capacity on
Vertical Drain Performance in Liquefiable Sand Deposits.
Pestana, J. M., Hunt, C. E., Goughnour, R. R., & Kammerer, A. M. (1997). Effect of Storage
Capacity on Vertical Drain Performance in Liquefiable Sand Deposits.
Prevost, J. (1985). A Simple Plasticity Theory for Frictional Cohesionless Soils. Soil Dynamics and
Earthquake Engineering , 9-17.
183
Randolph, M. F., Cassidy, M. J., Gourvenec, S. M., & Erbrich, C. (2005). Challenges of offshore
geotechnical engineering. Proceedings of the 16th International Conference on Soil Mechanics
and Geotechnical Engineering, (pp. 123-176). Osaka, Japan.
Remaud, D. (1999). Pieux sous charges latérales: Etude expérimentale de l'effet de groupe.
Thése de Doctorat de l'Université de Nantes.
Schofield, A. (1980). Cambridge geotechnical centrifuge operations. Géotechnique , 227-268.
Seed, H. B., & Booker, J. R. (1977). Stabilization of Potentially Liquefiable Sand Deposits using
Gravel Drains. Journal of the Geotechnical Engineering Division , 757-768.
Simo, J., & Hughes, T. (1998). Computational Inelasticity. New York: Springer-verlang.
Spencer, L., Rix, G., & Gallagher, P. (2007). Dynamic Properties of Colloidal Silica Gel and Sand
Mixtures. Fourth International Conference on Earthquake Geotechnical Engineering.
Stone, K., & Muir, W. (1992). Effects of dilatancy and particle size observed in model tests on
sand. Soils and Foundations , 43-47.
Ternet, O. (1999). Reconstitution et caractérization des massifs de sable: application aux essais
en centrifugeuse et en chambre de calibration. Thése de doctorat . Université de Caen.
Towhata, I. (2007). Developments of Soil Improvement Technologies for Mitigation of

Liquefaction Risk. 4th Internationl Conference on Earthquake Geotechnical Engineering (pp.
355-383). Thessaloniki: Springer.
Whittle, A., & Kavvadas, M. (1994). Formulation of MIT-E3 Constitutive Model for
Overconsolidated Clays. Journal of Geotechnical Engineering , 173-198.
Yang, Z., Elgamal, A., & Parra, E. (2003). Computational Model for Cyclic Mobility and
Associated Shear Deformation. Journal of Geotechnical and Geoenvironmental Engineering ,
1119-1127.
Zienkiewicz, O., Chan, A., Pastor, M., Schrefler, B., & Shiomi, T. (1999). Computational
Geomechanics. Chichester: John Wiley & Sons Ltd.
184

Antonios Vytiniotis Master Thesis 2009 PDF

Uploaded by

Copyright:

Available Formats

Antonios Vytiniotis Master Thesis 2009 PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Antonios Vytiniotis Master Thesis 2009 PDF

Uploaded by

Copyright:

Available Formats

Contents

ture are discussed.

designed against this occurring shear strength loss.

damages in the Kobe port were enormous (Figure 1-2).

in order to validate numerical models for geotechnical earthquake engineering. By means of

Figure 1-3 The UC-Davis geotechnical centrifuge

fuge, similar to the one shown in Figure 1-3.

modern building codes and simulation capabilities.

1.3 Pre-fabricated Vertical Drains

Figure 1-4 Earthquake drain installation

This project has five primary objectives:

ware framework of Opensees.

• Validate the implemented theory with centrifuge experiments performed at UC-Davis

2.1 Centrifuge testing

2.1.1 Basic Principles

2.1.2 Scaling laws

If we assume that the length is scaled by a factor of N:

then the acceleration is scaled:

Which can only be possible if the time is also scaled by:

The body forces

The average acceleration of water

to different bulk modulus compared to the prototype fluid.

changes surface tension, drying, and wetting characteristics.

2.1.3 Other similitude problems

footing settlement even though a rule of thumb of

the container, is greater than 30:

very important but are not discussed in this thesis.

2.2 Analytical and numerical research towards modeling drainage-

based soil improved methods

and assuming infinite drain permeability.

the plane strain model

diction are also shown.

presented (Figure 2-6).

improved with infinite permeability stone columns(Papadimitriou, Moutsopoulou, Bouckovalas,

permeability and the storage capacity.

3.1 Experiment outline

In order to evaluate the effectiveness of PV drains for liquefaction remediation (Earthquake

Figure 3-1 Conceptual diagram of the PV drains centrifuge model

orded acceleration time history.

perform parametric studies or to be used directly by practitioners in design. Additionally the

tion applied is better for validation of our numerical model.

3.2 Model Preparation

layer would be saturated.

transducers and the accelerometers moved significantly downwards during shaking.

3.3 Scale Factors

Table 3-1 Scale Factors

Quantity Prototype/Model Dimension

Displacement, Length 15/1

Acceleration, Gravity 1/15

Pressure, Stress 1/1

3.4 Limitations of recorded data and model

one takes into consideration their vertical movement.

4.1 Governing Equations

strain. Q is defined by:

with the definition of the effective stresses:

leration; this results in the so-called u-p-U approximation:

where U is the water displacement relative to the soil skeleton.

lation applications in earthquake engineering using finite element methods”(Mazzoni,

Software Pros Cons

R=[0 ; 0; qL/D_oned/B; -qL/D_oned];

R=[0 ; 0; qL/D_oned/B; -qL/D_oned];