PyLiq1. PDMY, FSP, updateStage problems

[bayram_aygun]

I am a master's student at Rice University/Civil Engineering Dept. I am
doing research on bridge-network reliability under the supervision of Dr.
Leonardo Duenas-Osorio. I have three questions about a simple OpenSees code
which consists of one pile (2 nonlinear beam-column elements) and a soil
column with 1 soil layer (loose sand). I will greatly appreciate if you have
a few minutes to offer some comments about the problems I encounter:


1)I used dummy nodes (2DOF) to make a connection between the pile nodes
(3DOF) and the soil nodes (2 DOF). In order to fix free ends of P-y elements
against rotation and vertical displacement, the soil nodes of corresponding
"zerolength" element must be fixed also against vertical displacements but
that poses the question: Is the simulation less realistic when we take this
approach where the soil cannot deform vertically? Is this step really
necessary?

2)When I write this syntax:

updateMaterials -material 3 bulkModulus 10000.67

Opensees throws an error saying: 'expected integer but got "10000.67" '.
It's not a problem to round the number but is this something that I should
worry about or should I just round it off to 10001. I'm working with
v.1.7.5Beta version -which Frank McKenna provided me with- because the
latest version 1.7.4 apparently has some bugs concerning soil update
properties.

3)I couldn't pull the gravity analysis off with the presence of soil and
pile elements (The penalty constraint causes some trouble and the analysis
doesn't converge. I also tried transformation constraint but it wouldn't
work either.)

First I conducted gravity analysis for the soil column + a "zerolength"
element and then I updated the PyLiq1, PDMY, FSP and defined pile elements
(and nodes of 3DOF). Finally I conducted the transient analysis. Is this the
right process when one's main concern is getting the dynamic response of
soil-pile system? Or the sequence should be different?

Lastly, could you provide me with any examples where PyLiq1, PDMY and FSPM
(FluidSolidPorousMaterial) are used in conjunction with piles? I try to do
the best I can with my code but I must admit that OpenSees soil modelling is
much more mysterious than I thought it would be

Thanks in advance for your time and kindness,

-Bayram AYGUN

[bayram_aygun]

#One Pile-Soil Interaction Model By Bayram Aygun, NOV2007
#Units (m,kN,ton)
#


wipe

set NumSteps 2000
set StepSize [expr 1.0/$NumSteps]

#Build Soil Model

model basic -ndm 2 -ndf 2


#define loose sand
set rho 1.7000; # Saturated soil mass density.
set refShearModul 55000.; # Reference low-strain shear modulus
set refBulkModul 150000.; # Reference bulk modulus
set frictionAng 29; # Friction angle at peak shear strength, in degrees.
set peakShearStra 0.1; # an octahedral shear strain at which the maximum shear strength is reached
set refPress 80. ; # Reference mean effective confining pressure at which Gr, Br, and ?max are defined.
set pressDependCoe 0.5; # A positive constant defining variations of G and B as a function of instantaneous effective confinement p':
set PTAng 29; # Phase transformation angle, in degrees.
set contrac 0.21; # A non-negative constant defining the rate of contraction
set dilat1 0.;
set dilat2 0.; # Non-negative constants defining the rate of
set liquefac1 10.; # Parameters controlling the mechanism of....
set liquefac2 0.02; # ...............
set liquefac3 1.0; #.. cyclic mobility
set e 0.85; # Initial void ratio

# $tag $nd

nDMaterial PressureDependMultiYield 3 2 $rho $refShearModul $refBulkModul $frictionAng $peakShearStra $refPress $pressDependCoe $PTAng 0.17 $dilat1 $dilat2 $liquefac1 $liquefac2 $liquefac3 $e


#define FSP

set combinedBulkModul_1 2.2D6; #Combined undrained bulk modulus for loose soil

#nDMaterial FluidSolidPorousMaterial $tag $nd $soilMatTag $combinedBulkModul

nDMaterial FluidSolidPorous 5 2 3 $combinedBulkModul_1


# Define PyLiq1

set SoilType 2; # Cross-sectional Area in^2; #soilType = 1 Backbone of p-y curve approximates...
#...Matlock (1970) soft clay relation.
#soilType = 2 Backbone of p-y curve approximates...
#.... API (1993) sand relation.

set P_ult_1 65.; #Lateral Capacity(Passive Resistance) of the weak soil layer per unit length
set Y_50_1 0.085; #the displacement at which 50% of the P_ult is mobilized.
set C_d_1 0.07; #Variable that sets the drag resistance within a fully-mobilized gap as Cd*pult.
set C_1 0.0; #Nonzero c values are used to represent radiation damping effects.
set pRes_1 6.5; #Minimum p-y resistance that the material retains

# uniaxialMaterial PyLiq1 $matTag $soilType $pult $Y50 $Cd $c $pRes $solidElem1 $solidElem2

uniaxialMaterial PyLiq1 7 $SoilType $P_ult_1 $Y_50_1 $C_d_1 $C_1 $pRes_1 1 1



updateMaterialStage -material 5 -stage 0 #FSPM
updateMaterialStage -material 3 -stage 0 #PDMY

updateMaterialStage -material 7 -stage 0



#Define Soil Nodes
# x y

node 1 0 0
node 2 10 0
node 3 10 10
node 4 0 10


set inclination 0
set massDen 2.0 ;# solid mass density
set fluidDen 1.0 ;# fluid mass density
set pi 3.1415926535

set unitWeightX [expr ($massDen-$fluidDen)*9.81*sin($inclination/180.0*$pi)] ;# unit weight in X direction
set unitWeightY [expr -($massDen-$fluidDen)*9.81*cos($inclination/180.0*$pi)] ;# unit weight in Y direction


# This command is used to construct a FourNodeQuad element object which uses a bilinear isoparametric formulation.

#element quad $eleTag $iNode $jNode $kNode $lNode $thick $type $matTag <$pressure $rho $b1 $b2>



element quad 1 1 2 3 4 1.0 "PlaneStrain" 5 0.0 0 $unitWeightX $unitWeightY


#dummy node
node 33 10 10

# Begin to define p-y elements

# element zeroLength $eleTag $iNode $jNode -mat $matTag1 $matTag2 ... -dir $dir1 $dir2

element zeroLength 47 33 3 -mat 7 -dir 1

# fix the soil base

fix 1 1 1
fix 2 1 1

#tie nodes 3 and 4

equalDOF 3 4 1 2

# fix nodes 3 and 4 in vertical direction
#fix 3 0 1
#fix 4 0 1

#fix the dummy node in vertical direction

fix 33 0 1

################################
# Gravity Application
###############################
system ProfileSPD
test NormUnbalance 1.0e-5 10 0
constraints Transformation
integrator LoadControl 1 1 1 1

algorithm Newton
numberer RCM

analysis Static

analyze 2


#switch the material to plastic

updateMaterialStage -material 5 -stage 1 #FSPM
updateMaterialStage -material 3 -stage 1 #PDMY

updateMaterials -material 3 bulkModulus 100000

updateMaterialStage -material 7 -stage 1 #PyLiq1


model basic -ndm 2 -ndf 3

node 5 10 0
node 6 10 10
node 7 10 20

#fix the pile tip in every possible DOF
fix 5 1 1 1

# fixing the pile nodes in the vertical direction

fix 6 0 1 0
fix 7 0 1 0

#create geometric transformation for pile elements

geomTransf Linear 1

# define materials for pile (for constitutive nonlinear beam-column elements)

#uniaxialMaterial Concrete02 $matTag $fpc $epsc0 $fpcu $epsU $lambda $ft $Ets
uniaxialMaterial Concrete02 1 -27580. -0.003 -5516. -0.01 0.1 3681. 1930544
#build cover concrete (unconfined)


#uniaxialMaterial Steel02 $matTag $Fy $E $b $R0 $cR1 $cR2 $a1 $a2 $a3 \
#$a4
#use Steel01 next time.!!!!!!!!!!!!!
uniaxialMaterial Steel02 2 460572.64 199949200. 0.01 18 0.925 0.15;# build rein. matr


# Define Fiber Sections

section Fiber 109 {

# patch quad $matTag $numSubdivIJ $numSubdivJK $yI $zI $yJ $zJ $yK $zK $yL $zL
patch quad 1 4 16 -0.30 0.3 -0.3 -0.3 0.3 -0.3 0.3 -0.3

# layer straight $matTag $numBars $areaBar $yStart $zStart $yEnd $zEnd
layer straight 2 16 0.00145 0.25 0.25 0.25 -0.25; # top layer rein
layer straight 2 16 0.00145 -0.25 0.25 -0.25 -0.25; # bottom layer rein

};



#element nonlinearBeamColumn $eleTag $iNode $jNode $numIntgrPts $secTag $transfTag <-mass $massDens>\
element nonlinearBeamColumn 17 5 6 5 109 1 #<-iter $maxIters $tol>
element nonlinearBeamColumn 39 6 7 5 109 1

## fix free ends of P-y elements against rotation and vertical displacements
#fix 3 0 1 1


#Slave the pile node in 3DOF to the dummy node in 2DOF in the x & y direction
#equalDOF $rNodeTag $cNodeTag $dof1 $dof2

equalDOF 33 3 1 2



pattern Plain 1 Linear { ; # define LoadPattern 1. impose load in a linear manner

load 7 0. -50. 0. 0. 0. 0.;

}

###########################
#Now DEFINE TRANSIENT lOADS
#############################

set accelSeries "Series -dt 0.01 -filePath BM68elc.txt -factor 1"; # define acceleration vector from file (dt=0.01 is associated with the input file gm)
pattern UniformExcitation 2 1 -accel $accelSeries; # define where and how (pattern tag, dof) acceleration is applied


#rezero time
setTime 0.0
wipe analysis #This command does not destroy the elements, nodes, materials, etc. It does...
#destroy the solution strategies: the algorithm, analysis, eq'n solver, constraint handler, etc.


constraints Transformation
test NormDispIncr 1.D-12 25 0
algorithm Newton
numberer RCM
system ProfileSPD

set NewmarkGamma 0.3025; # Newmark-integrator gamma parameter (also HHT)
set NewmarkBeta 0.25; # Newmark-integrator beta parameter
integrator Newmark $NewmarkGamma $NewmarkBeta
analysis Transient

recorder Node -file PileTopDisp.out -time -node 7 -dof 1 2 -dT 0.01 disp

[Gang Wang]

I donnot have whole lot of time to look into the code, but it seems to me the code contains some errors and a lot of misunderstanding. Just some answers to your questions:

See the system below:

**************** 7****
**************** |****
**************** |****
**************** |****
4_________ 3 _ (33) 6***
|**********| **** | ****
|**********| **** |****
|**********| **** |****
|**********| **** |****
1 _________| 2 * _5_***


You fixed nodes 1,2, and 5. You created zerLength element between 3 and 33. NOW YOU WANT TO TIE 33 WITH 6 (equalDOF 33 6 1 2 ), instead of 33 and 3 in your code (equalDOF 33 3 1 2), which doesn't make any sense to me.

You also tied 3 with 4, which means the soil can only deform by pure shear mode of deformation, which might be fine. BUT YOU DON'T NEED TO DO THIS. Also fix 33 0 1 shouldn't be there at the first place.

I suggest you to go step by step. Go through static gravity analysis first, before building into a little bit more complicated model, which is now in a state of chaos. I just wonder without a loadPattern, how the gravity can be loaded?

If gravity cannot be properly enforced, it would be a wonder if your model can converge. Since it is pressure-sensitive material. Without enough overburden pressure, the yield surface is just too "small" to converge.

[bayram_aygun]

Hi,

You are right. I should start with a simpler model maybe. As for the gravity analysis, I defined body forces (gravitational) in the element quad ( $b1), so OpenSees calculates the nodal loads consistent with these body forces. I saw this approach in one of Dr. Zhaohui Yang's examples (EX #15).

I will try to put an order into this "chaos" by beginning with a simpler code, I suppose.

Thanks again for your kindness,