14.4.2. Effective Stress Site Response Analysis of a Layered Soil ColumnΒΆ

  1. The original model can be found here.

  2. The Python code is converted by Harsh Mistry from The University of Manchester, UK (harsh.mistry@manchester.ac.uk), and shown below, which can be downloaded here.

  3. The Input motion (velocity-time) can be downloaded here.

  1# -*- coding: utf-8 -*-
  2"""
  3Created on Fri Jan 29 16:39:41 2021
  4
  5@author: harsh
  6"""
  7import numpy as np
  8import math as mm
  9import opensees as op
 10import time as tt
 11##################################################################
 12#                                                                #
 13# Effective stress site response analysis for a layered          #
 14# soil profile located on a 2% slope and underlain by an         #
 15# elastic half-space.  9-node quadUP elements are used.          #
 16# The finite rigidity of the elastic half space is               #
 17# considered through the use of a viscous damper at the          #
 18# base.                                                          #
 19#                                                                #
 20#    Converted to openseespy by: Harsh Mistry                    #
 21#                                The University of Manchester    #
 22#                                                                #
 23#   Created by:  Chris McGann                                    #
 24#                HyungSuk Shin                                   #
 25#                Pedro Arduino                                   #
 26#                Peter Mackenzie-Helnwein                        #
 27#              --University of Washington--                      #
 28#                                                                #
 29# ---> Basic units are kN and m   (unless specified)             #
 30#                                                                #
 31##################################################################
 32#-----------------------------------------------------------------------------------------
 33#  1. DEFINE SOIL AND MESH GEOMETRY
 34#-----------------------------------------------------------------------------------------
 35
 36op.wipe()
 37nodes_dict = dict()
 38
 39#---SOIL GEOMETRY
 40# thicknesses of soil profile (m)
 41soilThick = 30.0
 42# number of soil layers
 43numLayers = 3
 44# layer thicknesses
 45layerThick=[20.0,8.0,2.0]
 46
 47# depth of water table
 48waterTable = 2.0
 49
 50# define layer boundaries
 51layerBound=np.zeros((numLayers,1))
 52layerBound[0]=layerThick[0];
 53for i in range(1,numLayers):
 54    layerBound[i]=layerBound[i-1]+layerThick[i]
 55    
 56#---MESH GEOMETRY
 57# number of elements in horizontal direction
 58nElemX = 1
 59# number of nodes in horizontal direction
 60nNodeX =2 * nElemX+1
 61# horizontal element size (m)
 62sElemX = 2.0
 63
 64# number of elements in vertical direction for each layer
 65nElemY = [40,16,4]
 66
 67# total number of elements in vertical direction
 68nElemT = 60
 69
 70sElemY = np.zeros((numLayers,1))
 71# vertical element size in each layer
 72for i in range(numLayers):
 73    sElemY[i] = [layerThick[i-1]/nElemY[i-1]]
 74    print('size:',sElemY[i])
 75    
 76# number of nodes in vertical direction
 77nNodeY = 2 * nElemT+1
 78
 79# total number of nodes
 80nNodeT = nNodeX * nNodeY
 81
 82#-----------------------------------------------------------------------------------------
 83#  2. CREATE PORE PRESSURE NODES AND FIXITIES
 84#-----------------------------------------------------------------------------------------
 85op.model('basic', '-ndm', 2, '-ndf', 3) 
 86
 87count = 1
 88layerNodeCount = 0
 89dry_Node=np.zeros((500,1))
 90node_save=np.zeros((500,1))
 91# loop over soil layers
 92for k in range(1,numLayers+1):
 93    # loop in horizontal direction
 94    for i in range(1,nNodeX+1,2):
 95       if k==1:
 96           bump = 1
 97       else:
 98           bump = 0
 99       j_end=2 * nElemY[k-1] + bump    
100       for j in range(1,j_end+1,2):
101            xCoord  = (i-1) * (sElemX/2)
102            yctr    = j + layerNodeCount
103            yCoord  = (yctr-1) * (np.float(sElemY[k-1]))/2
104            nodeNum = i + ((yctr-1) * nNodeX)
105            op.node(nodeNum, xCoord, yCoord)
106
107          # output nodal information to data file
108            nodes_dict[nodeNum] = (nodeNum, xCoord, yCoord)
109            node_save[nodeNum] = np.int(nodeNum)
110          # designate nodes above water table
111            waterHeight = soilThick - waterTable
112            if yCoord >= waterHeight:
113                dry_Node[count] = nodeNum
114                count = count+1
115    layerNodeCount = yctr + 1
116
117dryNode=np.trim_zeros(dry_Node)
118Node_d=np.unique(node_save)        
119Node_d=np.trim_zeros(Node_d)
120np.savetxt('Node_record.txt',Node_d)    
121print('Finished creating all -ndf 3 nodes')
122print('Number of Dry Nodes:',len(dryNode))    
123
124# define fixities for pore pressure nodes above water table
125for i in range(count-1):
126    n_dryNode=np.int(dryNode[i])
127    op.fix(n_dryNode, 0, 0, 1)
128    
129op.fix(1, 0, 1, 0)
130op.fix(3, 0, 1, 0)    
131print('Finished creating all -ndf 3 boundary conditions...')
132
133# define equal degrees of freedom for pore pressure nodes
134for i in range(1,((3*nNodeY)-2),6):
135    op.equalDOF(i, i+2, 1, 2)
136    
137print("Finished creating equalDOF for pore pressure nodes...")
138
139#-----------------------------------------------------------------------------------------
140#  3. CREATE INTERIOR NODES AND FIXITIES
141#-----------------------------------------------------------------------------------------
142op.model('basic', '-ndm', 2, '-ndf', 2) 
143
144xCoord = np.float(sElemX/2)
145
146# loop over soil layers
147layerNodeCount = 0
148
149for k in range(1,numLayers+1):
150    if k==1:
151        bump = 1
152    else:
153        bump = 0
154    j_end=2 * nElemY[k-1] + bump    
155    for j in range(1,j_end+1,1):
156        yctr = j + layerNodeCount
157        yCoord  = (yctr-1) * (np.float(sElemY[k-1]))/2
158        nodeNum = (3*yctr) - 1 
159        op.node(nodeNum, xCoord, yCoord)
160        # output nodal information to data file
161        nodes_dict[nodeNum] = (nodeNum, xCoord, yCoord)
162    
163    layerNodeCount = yctr   
164
165# interior nodes on the element edges
166# loop over layers
167layerNodeCount = 0
168
169for k in range(1,numLayers+1):
170    # loop in vertical direction
171    for j in range(1,((nElemY[k-1])+1)):
172        yctr     = j + layerNodeCount;
173        yCoord   = np.float(sElemY[k-1])*(yctr-0.5)
174        nodeNumL = (6*yctr) - 2
175        nodeNumR = nodeNumL + 2
176        
177        op.node(nodeNumL ,0.0, yCoord)
178        op.node(nodeNumR , sElemX, yCoord)
179        
180        # output nodal information to data file
181        nodes_dict[nodeNumL] = (nodeNumL ,0.0, yCoord)
182        nodes_dict[nodeNumR] = (nodeNumR , sElemX, yCoord)
183    layerNodeCount = yctr
184
185print("Finished creating all -ndf 2 nodes...")        
186
187# define fixities for interior nodes at base of soil column
188op.fix(2, 0, 1)
189print('Finished creating all -ndf 2 boundary conditions...')
190
191# define equal degrees of freedom which have not yet been defined
192for i in range(1,((3*nNodeY)-6),6):
193    op.equalDOF(i  , i+1, 1, 2)
194    op.equalDOF(i+3, i+4, 1, 2)
195    op.equalDOF(i+3, i+5, 1, 2)
196
197op.equalDOF(nNodeT-2, nNodeT-1, 1, 2)
198print('Finished creating equalDOF constraints...')
199
200#-----------------------------------------------------------------------------------------
201#  4. CREATE SOIL MATERIALS
202#-----------------------------------------------------------------------------------------
203
204# define grade of slope (%)
205grade = 2.0
206slope = mm.atan(grade/100.0)
207g     = -9.81
208
209xwgt_var = g * (mm.sin(slope))
210ywgt_var = g * (mm.cos(slope))
211thick = [1.0,1.0,1.0]
212xWgt  = [xwgt_var, xwgt_var, xwgt_var] 
213yWgt  = [ywgt_var, ywgt_var, ywgt_var] 
214uBulk = [6.88E6,  5.06E6, 5.0E-6]
215hPerm = [1.0E-4, 1.0E-4, 1.0E-4]
216vPerm = [1.0E-4, 1.0E-4, 1.0E-4]
217
218
219# nDMaterial PressureDependMultiYield02
220# nDMaterial('PressureDependMultiYield02', matTag, nd, rho, refShearModul, refBulkModul,\
221#    frictionAng, peakShearStra, refPress, pressDependCoe, PTAng,\
222#        contrac[0], contrac[2], dilat[0], dilat[2], noYieldSurf=20.0,\
223#            *yieldSurf=[], contrac[1]=5.0, dilat[1]=3.0, *liquefac=[1.0,0.0],e=0.6, \
224#                *params=[0.9, 0.02, 0.7, 101.0], c=0.1)
225
226op.nDMaterial('PressureDependMultiYield02',3, 2, 1.8, 9.0e4, 2.2e5, 32, 0.1, \
227                                      101.0, 0.5, 26, 0.067, 0.23, 0.06, \
228                                      0.27, 20, 5.0, 3.0, 1.0, \
229                                      0.0, 0.77, 0.9, 0.02, 0.7, 101.0)
230
231op.nDMaterial('PressureDependMultiYield02', 2, 2, 2.24, 9.0e4, 2.2e5, 32, 0.1, \
232                                      101.0, 0.5, 26, 0.067, 0.23, 0.06, \
233                                      0.27, 20, 5.0, 3.0, 1.0, \
234                                      0.0, 0.77, 0.9, 0.02, 0.7, 101.0)
235    
236op.nDMaterial('PressureDependMultiYield02',1, 2, 2.45, 1.3e5, 2.6e5, 39, 0.1, \
237                                      101.0, 0.5, 26, 0.010, 0.0, 0.35, \
238                                      0.0, 20, 5.0, 3.0, 1.0, \
239                                      0.0, 0.47, 0.9, 0.02, 0.7, 101.0)    
240
241print("Finished creating all soil materials...")
242
243#-----------------------------------------------------------------------------------------
244#  5. CREATE SOIL ELEMENTS
245#-----------------------------------------------------------------------------------------
246
247for j in range(1,nElemT+1):
248    nI = ( 6*j) - 5
249    nJ = nI + 2
250    nK = nI + 8
251    nL = nI + 6
252    nM = nI + 1
253    nN = nI + 5
254    nP = nI + 7
255    nQ = nI + 3
256    nR = nI + 4
257    
258    lowerBound = 0.0
259    for i in range(1,numLayers+1):
260        if j * sElemY[i-1] <= layerBound[i-1] and j * sElemY[i-1] > lowerBound:
261            # permeabilities are initially set at 1.0 m/s for gravity analysis,
262            op.element('9_4_QuadUP', j, nI, nJ, nK, nL, nM, nN, nP, nQ, nR, \
263                           thick[i-1], i, uBulk[i-1], 1.0, 1.0, 1.0, xWgt[i-1], yWgt[i-1])
264                
265        lowerBound = layerBound[i-1]
266            
267print("Finished creating all soil elements...")
268#-----------------------------------------------------------------------------------------
269#  6. LYSMER DASHPOT
270#-----------------------------------------------------------------------------------------
271
272# define dashpot nodes
273dashF =  nNodeT+1
274dashS =  nNodeT+2
275
276op.node(dashF,  0.0, 0.0)
277op.node(dashS,  0.0, 0.0)
278
279# define fixities for dashpot nodes
280op.fix(dashF, 1, 1)
281op.fix(dashS, 0, 1)
282
283# define equal DOF for dashpot and base soil node
284op.equalDOF(1, dashS,  1)
285print('Finished creating dashpot nodes and boundary conditions...')
286
287# define dashpot material
288colArea      = sElemX * thick[0]
289rockVS       = 700.0
290rockDen      = 2.5
291dashpotCoeff = rockVS * rockDen
292
293#uniaxialMaterial('Viscous', matTag, C, alpha)
294op.uniaxialMaterial('Viscous', numLayers+1, dashpotCoeff * colArea, 1)
295
296# define dashpot element
297op.element('zeroLength', nElemT+1, dashF, dashS, '-mat', numLayers+1, '-dir', 1)
298
299print("Finished creating dashpot material and element...")
300
301#-----------------------------------------------------------------------------------------
302#  7. CREATE GRAVITY RECORDERS
303#-----------------------------------------------------------------------------------------
304
305# create list for pore pressure nodes
306load_nodeList3=np.loadtxt('Node_record.txt')
307nodeList3=[]
308
309for i in range(len(load_nodeList3)):
310    nodeList3.append(np.int(load_nodeList3[i]))
311# record nodal displacment, acceleration, and porepressure
312op.recorder('Node','-file','Gdisplacement.txt','-time','-node',*nodeList3,'-dof', 1, 2, 'disp')
313op.recorder('Node','-file','Gacceleration.txt','-time','-node',*nodeList3,'-dof', 1, 2, 'accel')
314op.recorder('Node','-file','GporePressure.txt','-time','-node',*nodeList3,'-dof', 3, 'vel')
315
316# record elemental stress and strain (files are names to reflect GiD gp numbering)
317op.recorder('Element','-file','Gstress1.txt','-time','-eleRange', 1,nElemT,'material','1','stress')
318op.recorder('Element','-file','Gstress2.txt','-time','-eleRange', 1,nElemT,'material','2','stress')
319op.recorder('Element','-file','Gstress3.txt','-time','-eleRange', 1,nElemT,'material','3','stress')
320op.recorder('Element','-file','Gstress4.txt','-time','-eleRange', 1,nElemT,'material','4','stress')
321op.recorder('Element','-file','Gstress9.txt','-time','-eleRange', 1,nElemT,'material','9','stress')
322op.recorder('Element','-file','Gstrain1.txt','-time','-eleRange', 1,nElemT,'material','1','strain')
323op.recorder('Element','-file','Gstrain2.txt','-time','-eleRange', 1,nElemT,'material','2','strain')
324op.recorder('Element','-file','Gstrain3.txt','-time','-eleRange', 1,nElemT,'material','3','strain')
325op.recorder('Element','-file','Gstrain4.txt','-time','-eleRange', 1,nElemT,'material','4','strain')
326op.recorder('Element','-file','Gstrain9.txt','-time','-eleRange', 1,nElemT,'material','9','strain')
327
328print("Finished creating gravity recorders...")
329
330#-----------------------------------------------------------------------------------------
331#  8. DEFINE ANALYSIS PARAMETERS
332#-----------------------------------------------------------------------------------------
333
334#---GROUND MOTION PARAMETERS
335# time step in ground motion record
336motionDT = 0.005
337# number of steps in ground motion record
338motionSteps = 7990
339
340#---RAYLEIGH DAMPING PARAMETERS
341# damping ratio
342damp = 0.02
343# lower frequency
344omega1 = 2 * np.pi * 0.2
345# upper frequency
346omega2 = 2 * np.pi * 20
347# damping coefficients
348a0 = 2*damp*omega1*omega2/(omega1 + omega2)
349a1 = 2*damp/(omega1 + omega2)
350print("Damping Coefficients: a_0 = $a0;  a_1 = $a1")
351
352#---DETERMINE STABLE ANALYSIS TIME STEP USING CFL CONDITION
353# maximum shear wave velocity (m/s)
354vsMax = 250.0
355# duration of ground motion (s)
356duration = motionDT*motionSteps
357# minimum element size
358minSize = sElemY[0]
359
360for i in range(2,numLayers+1):
361    if sElemY[i-1] <= minSize:
362        minSize = sElemY[i-1]
363
364# trial analysis time step
365kTrial = minSize/(vsMax**0.5)
366# define time step and number of steps for analysis
367if motionDT <= kTrial:
368    nSteps = motionSteps
369    dT     = motionDT
370else:
371    nSteps = np.int(mm.floor(duration/kTrial)+1)
372    dT     = duration/nSteps
373    
374
375print("Number of steps in analysis: $nSteps")
376print("Analysis time step: $dT")
377
378#---ANALYSIS PARAMETERS
379# Newmark parameters
380gamma = 0.5
381beta  = 0.25
382
383#-----------------------------------------------------------------------------------------
384#  9. GRAVITY ANALYSIS
385#-----------------------------------------------------------------------------------------
386# update materials to ensure elastic behavior
387op.updateMaterialStage('-material', 1, '-stage', 0)
388op.updateMaterialStage('-material', 2, '-stage', 0)
389op.updateMaterialStage('-material', 3, '-stage', 0)
390
391op.constraints('Penalty', 1.0E14, 1.0E14)
392op.test('NormDispIncr', 1e-4, 35, 1)
393op.algorithm('KrylovNewton')
394op.numberer('RCM')
395op.system('ProfileSPD')
396op.integrator('Newmark', gamma, beta)
397op.analysis('Transient')
398
399startT = tt.time()
400op.analyze(10, 5.0E2)
401print('Finished with elastic gravity analysis...')
402
403# update material to consider elastoplastic behavior
404op.updateMaterialStage('-material', 1, '-stage', 1)
405op.updateMaterialStage('-material', 2, '-stage', 1)
406op.updateMaterialStage('-material', 3, '-stage', 1)
407
408# plastic gravity loading
409op.analyze(40, 5.0e2)
410
411print('Finished with plastic gravity analysis...')
412
413#-----------------------------------------------------------------------------------------
414#  10. UPDATE ELEMENT PERMEABILITY VALUES FOR POST-GRAVITY ANALYSIS
415#-----------------------------------------------------------------------------------------
416
417# choose base number for parameter IDs which is higer than other tags used in analysis
418ctr = 10000.0
419# loop over elements to define parameter IDs 
420for i in range(1,nElemT+1):
421    op.parameter(np.int(ctr+1.0), 'element', i, 'vPerm')
422    op.parameter(np.int(ctr+2.0), 'element', i, 'hPerm')
423    ctr   = ctr+2.0
424
425# update permeability parameters for each element using parameter IDs
426ctr = 10000.0
427for j in range(1,nElemT+1):
428    lowerBound = 0.0
429    for i in range(1,numLayers+1):
430        if j * sElemY[i-1] <= layerBound[i-1] and j*sElemY[i-1] > lowerBound:
431            op.updateParameter(np.int(ctr+1.0), vPerm[i-1])
432            op.updateParameter(np.int(ctr+2.0), hPerm[i-1])
433        lowerBound = layerBound[i-1]
434    ctr = ctr+2.0
435
436print("Finished updating permeabilities for dynamic analysis...")
437
438#-----------------------------------------------------------------------------------------
439#  11. CREATE POST-GRAVITY RECORDERS
440#-----------------------------------------------------------------------------------------
441
442# reset time and analysis
443op.setTime(0.0)
444op.wipeAnalysis()
445op.remove('recorders')
446
447# recorder time step
448recDT = 10*motionDT
449
450# record nodal displacment, acceleration, and porepressure
451op.recorder('Node','-file','displacement.txt','-time', '-dT',recDT,'-node',*nodeList3,'-dof', 1, 2, 'disp')
452op.recorder('Node','-file','acceleration.txt','-time', '-dT',recDT,'-node',*nodeList3,'-dof', 1, 2, 'accel')
453op.recorder('Node','-file','porePressure.txt','-time', '-dT',recDT,'-node',*nodeList3,'-dof', 3, 'vel')
454
455# record elemental stress and strain (files are names to reflect GiD gp numbering)
456op.recorder('Element','-file','stress1.txt','-time', '-dT',recDT,'-eleRange', 1,nElemT,'material','1','stress')
457op.recorder('Element','-file','stress2.txt','-time', '-dT',recDT,'-eleRange', 1,nElemT,'material','2','stress')
458op.recorder('Element','-file','stress3.txt','-time', '-dT',recDT,'-eleRange', 1,nElemT,'material','3','stress')
459op.recorder('Element','-file','stress4.txt','-time', '-dT',recDT,'-eleRange', 1,nElemT,'material','4','stress')
460op.recorder('Element','-file','stress9.txt','-time', '-dT',recDT,'-eleRange', 1,nElemT,'material','9','stress')
461op.recorder('Element','-file','strain1.txt','-time', '-dT',recDT,'-eleRange', 1,nElemT,'material','1','strain')
462op.recorder('Element','-file','strain2.txt','-time', '-dT',recDT,'-eleRange', 1,nElemT,'material','2','strain')
463op.recorder('Element','-file','strain3.txt','-time', '-dT',recDT,'-eleRange', 1,nElemT,'material','3','strain')
464op.recorder('Element','-file','strain4.txt','-time', '-dT',recDT,'-eleRange', 1,nElemT,'material','4','strain')
465op.recorder('Element','-file','strain9.txt','-time', '-dT',recDT,'-eleRange', 1,nElemT,'material','9','strain')
466
467print("Finished creating all recorders...")
468
469#-----------------------------------------------------------------------------------------
470#  12. DYNAMIC ANALYSIS
471#-----------------------------------------------------------------------------------------
472op.model('basic', '-ndm', 2, '-ndf', 3)
473
474# define constant scaling factor for applied velocity
475cFactor = colArea * dashpotCoeff
476
477# define velocity time history file
478velocityFile='velocityHistory';
479data_gm=np.loadtxt('velocityHistory.txt')
480#motionSteps=len(data_gm)
481#print('Number of point for GM:',motionSteps)
482
483# timeseries object for force history
484op.timeSeries('Path', 2, '-dt', motionDT, '-filePath', velocityFile+'.txt', '-factor', cFactor)
485op.pattern('Plain', 10, 2)
486op.load(1, 1.0, 0.0, 0.0)
487
488print( "Dynamic loading created...")
489
490op.constraints('Penalty', 1.0E16, 1.0E16)
491op.test('NormDispIncr', 1e-3, 35, 1)
492op.algorithm('KrylovNewton')
493op.numberer('RCM')
494op.system('ProfileSPD')
495op.integrator('Newmark', gamma, beta)
496op.rayleigh(a0, a1, 0.0, 0.0)
497op.analysis('Transient')
498
499# perform analysis with timestep reduction loop
500ok = op.analyze(nSteps,dT)
501
502# if analysis fails, reduce timestep and continue with analysis
503if ok !=0:
504    print("did not converge, reducing time step")
505    curTime = op.getTime()
506    mTime = curTime
507    print("curTime: ", curTime)
508    curStep = curTime/dT
509    print("curStep: ", curStep)
510    rStep   = (nSteps-curStep)*2.0
511    remStep = np.int((nSteps-curStep)*2.0)
512    print("remStep: ", remStep)
513    dT = dT/2.0
514    print("dT: ", dT)
515
516    ok = op.analyze(remStep, dT)
517    # if analysis fails again, reduce timestep and continue with analysis    
518    if ok !=0:
519        print("did not converge, reducing time step")
520        curTime = op.getTime()
521        print("curTime: ", curTime)
522        curStep = (curTime-mTime)/dT
523        print("curStep: ", curStep)
524        remStep = np.int((rStep-curStep)*2.0)
525        print("remStep: ", remStep)
526        dT = dT/2.0
527        print("dT: ", dT)
528    
529        ok = op.analyze(remStep, dT)    
530
531endT = tt.time()
532print("Finished with dynamic analysis...")
533print("Analysis execution time: ",(endT-startT))
534op.wipe()