Problem with a 3-span continuous beam bridge

[gaoyuanw]

Hi guys, I am encountered with a problem.
I am modeling a 3-span bridge and tried to applied lateral loads & ground motion.
The thing is that the X-direction(along the longitudinal length) loading is fine for both lateral & ground motion.
But when it comes to the Y-direction loading, Transient and Static analysis can't be completed.
I tried to apply lateral loads to pier 1 along the Y-direction, it succeeded. However, when I tried to apply Y-direction lateral loads to pier 2-4, it failed.

WARNING BandGenLinLapackSolver::solve() -factorization failed, matrix singular U(i,i) = 0, i= 24
WARNING Linear::solveCurrentStep() -the LinearSOE failed in solve()
StaticAnalysis::analyze() - the Algorithm failed at iteration: 1 with domain at load factor 0.2
OpenSees > analyze failed, returned: -3 error flag

I did the Eigen analysis, the results showed that the MODES are normal. I don't know where the problem is.

Hope u guys are safe and well. Looking forward to your reply!

[selimgunay]

Are you modeling your piers with fixed supports or simply supported at the bottom?

[gaoyuanw]

selimgunay wrote:
> Are you modeling your piers with fixed supports or simply supported at the
> bottom?

Dear selimgunay:

Thanks for your reply !
I was modeling with fixed supports.... I somehow solved the problem by using only 2 nodes (1 element) for the pier(dispBeamColumn), initially there are 3 elements.... But I don't know why.

I am now encountered with a new problem. I am trying to simulate soil-pile interaction for this bridge by using P-Y,Q-Z and T-z springs. I created pile nodes, string nodes, zero length elements to simulate strings, and equal dof to connect pile and strings.
I can now apply static loads, but I can't apply ground motion. I have tried several algorithms....
I am wondering which kind of pattern should I use (Uniform or MultipleSupport)? And where to apply the ground motion.....(Because I didn't simulate the soil, only the interaction........)

My Python codes are as follows: (Only One-Pier-model was created for simplification)


%matplotlib notebook
from openseespy.postprocessing.Get_Rendering import *
import openseespy.opensees as op
from openseespy.opensees import *
import numpy as np
import matplotlib
import math
import ReadRecord
wipe()
########################################## Environment
model('basic','-ndm',3,'-ndf', 6)

#单位：KN,M,T
#混凝土重度25KN/m**3
Pier1X=0
BeamNodeMass=100
beamMass = [BeamNodeMass,BeamNodeMass,BeamNodeMass,0,0,0]
PilecapMass = [96,96,96,0,0,0]
Pier13Mass = [34/2,34/2,34/2,0,0,0]
z1=10.000 #主梁质心坐标
z2=9.100 #墩顶坐标
z3=1.100 #墩底坐标
z4=0.000 #承台质心坐标
pcwidth= 2.000#桩间距
pclength=2.000 #桩间距（mm)
#Bridge Parameter(unknown),以地面（承台中心）为z向0点

op.node(1+10000,Pier1X,0,z1) #主梁质点（10000、20000代表一号桩……）

#node(2,0,0,z1) #支座节点
op.node(3+10000,Pier1X,0,z2) #墩顶
op.node(4+10000,Pier1X,0,z3) #墩底
op.node(5+10000,Pier1X,0,z4) #承台中心
#主梁

op.mass(1+10000,*beamMass)
op.mass(3+10000,*Pier13Mass)
op.mass(4+10000,*Pier13Mass)
op.mass(5+10000,*PilecapMass)

#主梁质点，承台质点

#-------------定义支座(固定),
#uniaxialMaterial('Elastic',1,5000000000)
#supmatTags=[1,1,1,1,1,1]
#supdir=[1,2,3,4,5,6]
#element('zeroLength',1,1,2,'-mat',*supmatTags,'-dir',*supdir)
#固定支座定义结束--------------



#-------------定义刚臂(梁中心和墩顶)
vecxz0=[0,1,1]
op.geomTransf('Linear',1,*vecxz0)
op.rigidLink('beam',10001,10003)

#刚臂单元定义结束--------------

print("主梁质点+梁墩刚臂定义完成")

#-------------Beginning of Pier Element setup
#Define materials for nonlinear columns
#Pier Cross section:Cover Concrete, Core Concrete, Steel Bars

#Core Concrete Para

#Kpa
fpc1=-31600
ec01=-0.002
fpcu1=-26860
ecu1=-0.0038

#Cover Concrete Para
fpc2=-26333
ec02=-0.001
fpcu2=-0.0
ecu2=-0.002

#Steel Bars Para
fy1=400*1000
E1=2*10**8
b=0.0

#Material for Pier
op.uniaxialMaterial('Concrete01',2,fpc1,ec01,fpcu1,ecu1)
op.uniaxialMaterial('Concrete01',3,fpc2,ec02,fpcu2,ecu2)
op.uniaxialMaterial('Steel01',4,fy1,E1,b)

#Cross-section parameter for pier
PierY=2.2 #横桥向墩宽
PierX=1.5 #纵桥向墩宽
Barspos1=0.300
Cover01=0.250
As1=0.000490625 #(φ25)
Barsnum=58
#Cross-section Fiber Assembly

op.section('Fiber',1,'-torsion',2)
op.patch('rect',2,10,1,-0.5,-0.85,0.5,0.85)
op.patch('rect',3,10,1,0.5,-1.1,0.75,1.1)
op.patch('rect',3,10,1,-0.75,-1.1,-0.5,1.1)
op.patch('rect',3,2,1,-0.5,-1.1,0.5,-0.85)
op.patch('rect',3,2,1,-0.5,0.85,0.5,1.1)


#核心混凝土+混凝土保护层
op.layer('straight',4,Barsnum/5*1.5,As1,0.45,-0.8,0.45,0.
op.layer('straight',4,Barsnum/5*1.5,As1,-0.45,-0.8,-0.45,0.
op.layer('straight',4,Barsnum/5*1,As1,-0.45,-0.8,0.45,-0.
op.layer('straight',4,Barsnum/5*1,As1,-0.45,0.8,0.45,0.

print("墩柱截面设置完毕")
#局部坐标-全局坐标转换
transfType1 = 'Linear'
transfTag1 = 2
vecxz1=[0.0,1.0,1.0] #与X的叉乘是局部的y

op.geomTransf(transfType1,transfTag1,*vecxz1)
op.beamIntegration('Lobatto',1,1,5) #积分：Lobatto，1,1号截面，5个积分点
op.element('dispBeamColumn',10002,10003,10004,transfTag1,1)
#op.element('dispBeamColumn',10003,10056,10057,transfTag1,1)
#op.element('dispBeamColumn',10004,10057,10004,transfTag1,1)
#op.element('dispBeamColumn',20002,20003,20004,transfTag1,1)
#op.element('dispBeamColumn',20003,20056,20057,transfTag1,1)
#op.element('dispBeamColumn',20004,20057,20004,transfTag1,1)
#op.element('dispBeamColumn',30002,30003,30004,transfTag1,1)
#op.element('dispBeamColumn',30003,30056,30057,transfTag1,1)
#op.element('dispBeamColumn',30004,30057,30004,transfTag1,1)
#op.element('dispBeamColumn',40002,40003,40004,transfTag1,1)
#element('dispBeamColumn',40003,40056,40057,transfTag1,1)
#element('dispBeamColumn',40004,40057,40004,transfTag1,1)


#Pier setup finished------------------------
print("墩柱单元设置完毕")

#--------------承台与墩底连接刚臂
geomTransf('Linear',3,*vecxz0)
element('elasticBeamColumn',10005,10004,10005,1400,2*10**8,76903069,5,1429,1867,3)

#承台质点与墩底刚臂单元定义结束--------------


print("上部结构建立完成")



plot_model()



#-------------桩、土弹簧(py,qz,tz)建模，其中利用华盛顿大学的既有代码计算侧向土的抵抗力。
L1= 1.100 #高于地面的桩长
L2= 55.000 #插入土体的桩长
D1= 1.0 #桩径
PileEleNum= 70 #桩单元个数
EleSize=(L1+L2)/PileEleNum #桩单元尺寸
PileNodeNum= PileEleNum+1 #桩节点个数

PileNum=4 #桩根数
XLength=2.000 #X向桩间距
YLength=2.000 #Y向桩间距
######################土弹簧节点建模

#计数器
count =0

#桩节点土弹簧:node 从桩顶开始编号(1101x向,1201y向,1301z向)
#土弹簧节点的自由度为3,均为平动
model('basic','-ndm',3,'-ndf',3)
pile1_nodes = dict()
for j in range(1):
for i in range(PileNodeNum):
if (i*EleSize>L1 and i*EleSize<=(L1+L2)):
node((j+1)*10000+i+1101,j*25+0.5*XLength,0.5*YLength,-i*EleSize)
#node((j+1)*10000+i+1201,j*25+0.5*XLength,0.5*YLength,-i*EleSize)
node((j+1)*10000+i+1301,25*j+0.5*XLength,0.5*YLength,-i*EleSize)
node((j+1)*10000+i+2101,25*j+0.5*XLength,-0.5*YLength,-i*EleSize)
#node((j+1)*10000+i+2201,25*j+0.5*XLength,-0.5*YLength,-i*EleSize)
node((j+1)*10000+i+2301,25*j+0.5*XLength,-0.5*YLength,-i*EleSize)
node((j+1)*10000+i+3101,25*j-0.5*XLength,-0.5*YLength,-i*EleSize)
#node((j+1)*10000+i+3201,25*j-0.5*XLength,-0.5*YLength,-i*EleSize)
node((j+1)*10000+i+3301,25*j-0.5*XLength,-0.5*YLength,-i*EleSize)
node((j+1)*10000+i+4101,25*j-0.5*XLength,0.5*YLength,-i*EleSize)
#node((j+1)*10000+i+4201,25*j-0.5*XLength,0.5*YLength,-i*EleSize)
node((j+1)*10000+i+4301,25*j-0.5*XLength,0.5*YLength,-i*EleSize)
pile1_nodes[i+11101]=(0.5*XLength,0.5*YLength,-i*EleSize)
pile1_nodes[i+11201]=(0.5*XLength,0.5*YLength,-i*EleSize)
pile1_nodes[i+11301]=(0.5*XLength,0.5*YLength,-i*EleSize)
count = count +1

print("Finished creating all spring nodes...")
print(count)
# 在地面下的土弹簧节点个数
NodeEmbedNum = count

#土弹簧固定性:tz和qz、y向的yz全部固定，x方向的yz约束2个，之后用主从约束固定到桩身上
for j in range(1):
for i in range(PileNodeNum):
if (i*EleSize>L1 and i*EleSize<=(L1+L2)):
fix((j+1)*10000+i+1101,0,0,1)
#fix((j+1)*10000+i+1201,1,0,1)
fix((j+1)*10000+i+1301,1,1,1)
fix((j+1)*10000+i+2101,0,0,1)
#fix((j+1)*10000+i+2201,1,0,1)
fix((j+1)*10000+i+2301,1,1,1)
fix((j+1)*10000+i+3101,0,0,1)
#fix((j+1)*10000+i+3201,1,0,1)
fix((j+1)*10000+i+3301,1,1,1)
fix((j+1)*10000+i+4101,0,0,1)
#fix((j+1)*10000+i+4201,1,0,1)
fix((j+1)*10000+i+4301,1,1,1)

print("Finished creating all spring node fixities...")

#-----土体性质
#soil unit weight (kN/m^3)
gamma = 17.0
# soil internal friction angle (degrees)
phi = 36.0
# soil shear modulus at pile tip (kPa)
Gsoil = 150000.0

# select pult definition method for p-y curves
# API (default) --> 1
# Brinch Hansen --> 2
puSwitch = 1

# variation in coefficent of subgrade reaction with depth for p-y curves
# API linear variation (default) --> 1
# modified API parabolic variation --> 2
kSwitch = 1

# effect of ground water on subgrade reaction modulus for p-y curves
# above gwt --> 1
# below gwt --> 2
gwtSwitch = 1



##########################################################
# 华盛顿大学获得py，qz，tz弹簧参数的源代码 #
##########################################################
# #
# Procedure to compute ultimate lateral resistance, p_u, #
# and displacement at 50% of lateral capacity, y50, for #
# p-y springs representing cohesionless soil. #
# Converted to openseespy by: Pavan Chigullapally #
# University of Auckland #
# #
# Created by: Hyung-suk Shin #
# University of Washington #
# Modified by: Chris McGann #
# Pedro Arduino #
# Peter Mackenzie-Helnwein #
# University of Washington #
# #
###########################################################

# references
# American Petroleum Institute (API) (1987). Recommended Practice for Planning, Designing and
# Constructing Fixed Offshore Platforms. API Recommended Practice 2A(RP-2A), Washington D.C,
# 17th edition.
#
# Brinch Hansen, J. (1961). "The ultimate resistance of rigid piles against transversal forces."
# Bulletin No. 12, Geoteknisk Institute, Copenhagen, 59.
#
# Boulanger, R. W., Kutter, B. L., Brandenberg, S. J., Singh, P., and Chang, D. (2003). Pile
# Foundations in liquefied and laterally spreading ground during earthquakes: Centrifuge experiments
# and analyses. Center for Geotechnical Modeling, University of California at Davis, Davis, CA.
# Rep. UCD/CGM-03/01.
#
# Reese, L.C. and Van Impe, W.F. (2001), Single Piles and Pile Groups Under Lateral Loading.
# A.A. Balkema, Rotterdam, Netherlands.



def get_pyParam ( pyDepth, gamma, phiDegree, b, pEleLength, puSwitch, kSwitch, gwtSwitch):

#----------------------------------------------------------
# define ultimate lateral resistance, pult
#----------------------------------------------------------

# pult is defined per API recommendations (Reese and Van Impe, 2001 or API, 1987) for puSwitch = 1
# OR per the method of Brinch Hansen (1961) for puSwitch = 2

pi = 3.14159265358979
phi = phiDegree * (pi/180)
zbRatio = pyDepth / b

#-------API recommended method-------

if puSwitch == 1:

# obtain loading-type coefficient A for given depth-to-diameter ratio zb
# ---> values are obtained from a figure and are therefore approximate
zb = []
dataNum = 41
for i in range(dataNum):
b1 = i * 0.125
zb.append(b1)
As = [2.8460, 2.7105, 2.6242, 2.5257, 2.4271, 2.3409, 2.2546, 2.1437, 2.0575, 1.9589, 1.8973, 1.8111, 1.7372, 1.6632, 1.5893, 1.5277, 1.4415, 1.3799, 1.3368, 1.2690, 1.2074, 1.1581,
1.1211, 1.0780, 1.0349, 1.0164, 0.9979, 0.9733, 0.9610, 0.9487, 0.9363, 0.9117, 0.8994, 0.8994, 0.8871, 0.8871, 0.8809, 0.8809, 0.8809, 0.8809, 0.8809]

# linear interpolation to define A for intermediate values of depth:diameter ratio
for i in range(dataNum):
if zbRatio >= 5.0:
A = 0.88
elif zb[i] <= zbRatio and zbRatio <= zb[i+1]:
A = (As[i+1] - As[i])/(zb[i+1] - zb[i]) * (zbRatio-zb[i]) + As[i]

# define common terms
alpha = phi / 2
beta = pi / 4 + phi / 2
K0 = 0.4

tan_1 = math.tan(pi / 4 - phi / 2)
Ka = math.pow(tan_1 , 2)

# terms for Equation (3.44), Reese and Van Impe (2001)
tan_2 = math.tan(phi)
tan_3 = math.tan(beta - phi)
sin_1 = math.sin(beta)
cos_1 = math.cos(alpha)
c1 = K0 * tan_2 * sin_1 / (tan_3*cos_1)

tan_4 = math.tan(beta)
tan_5 = math.tan(alpha)
c2 = (tan_4/tan_3)*tan_4 * tan_5

c3 = K0 * tan_4 * (tan_2 * sin_1 - tan_5)

c4 = tan_4 / tan_3 - Ka

# terms for Equation (3.45), Reese and Van Impe (2001)
pow_1 = math.pow(tan_4,8)
pow_2 = math.pow(tan_4,4)
c5 = Ka * (pow_1-1)
c6 = K0 * tan_2 * pow_2

# Equation (3.44), Reese and Van Impe (2001)
pst = gamma * pyDepth * (pyDepth * (c1 + c2 + c3) + b * c4)

# Equation (3.45), Reese and Van Impe (2001)
psd = b * gamma * pyDepth * (c5 + c6)

# pult is the lesser of pst and psd. At surface, an arbitrary value is defined
if pst <=psd:
if pyDepth == 0:
pu = 0.01

else:
pu = A * pst

else:
pu = A * psd

# PySimple1 material formulated with pult as a force, not force/length, multiply by trib. length
pult = pu * pEleLength

#-------Brinch Hansen method-------
elif puSwitch == 2:
# pressure at ground surface
cos_2 = math.cos(phi)

tan_6 = math.tan(pi/4+phi/2)

sin_2 = math.sin(phi)
sin_3 = math.sin(pi/4 + phi/2)

exp_1 = math.exp((pi/2+phi)*tan_2)
exp_2 = math.exp(-(pi/2-phi) * tan_2)

Kqo = exp_1 * cos_2 * tan_6 - exp_2 * cos_2 * tan_1
Kco = (1/tan_2) * (exp_1 * cos_2 * tan_6 - 1)

# pressure at great depth
exp_3 = math.exp(pi * tan_2)
pow_3 = math.pow(tan_2,4)
pow_4 = math.pow(tan_6,2)
dcinf = 1.58 + 4.09 * (pow_3)
Nc = (1/tan_2)*(exp_3)*(pow_4 - 1)
Ko = 1 - sin_2
Kcinf = Nc * dcinf
Kqinf = Kcinf * Ko * tan_2

# pressure at an arbitrary depth
aq = (Kqo/(Kqinf - Kqo))*(Ko*sin_2/sin_3)
KqD = (Kqo + Kqinf * aq * zbRatio)/(1 + aq * zbRatio)

# ultimate lateral resistance
if pyDepth == 0:
pu = 0.01
else:
pu = gamma * pyDepth * KqD * b

# PySimple1 material formulated with pult as a force, not force/length, multiply by trib. length
pult = pu * pEleLength

#----------------------------------------------------------
# define displacement at 50% lateral capacity, y50
#----------------------------------------------------------

# values of y50 depend of the coefficent of subgrade reaction, k, which can be defined in several ways.
# for gwtSwitch = 1, k reflects soil above the groundwater table
# for gwtSwitch = 2, k reflects soil below the groundwater table
# a linear variation of k with depth is defined for kSwitch = 1 after API (1987)
# a parabolic variation of k with depth is defined for kSwitch = 2 after Boulanger et al. (2003)

# API (1987) recommended subgrade modulus for given friction angle, values obtained from figure (approximate)

ph = [28.8, 29.5, 30.0, 31.0, 32.0, 33.0, 34.0, 35.0, 36.0, 37.0, 38.0, 39.0, 40.0]

# subgrade modulus above the water table
if gwtSwitch == 1:
k = [10, 23, 45, 61, 80, 100, 120, 140, 160, 182, 215, 250, 275]

else:
k = [10, 20, 33, 42, 50, 60, 70, 85, 95, 107, 122, 141, 155]

dataNum = 13
for i in range(dataNum):
if ph[i] <= phiDegree and phiDegree <= ph[i+1]:
khat = (k[i+1]-k[i])/(ph[i+1]-ph[i])*(phiDegree - ph[i]) + k[i]

# change units from (lb/in^3) to (kN/m^3)
k_SIunits = khat * 271.45

# define parabolic distribution of k with depth if desired (i.e. lin_par switch == 2)
sigV = pyDepth * gamma

if sigV == 0:
sigV = 0.01

if kSwitch == 2:
# Equation (5-16), Boulanger et al. (2003)
cSigma = math.pow(50 / sigV , 0.5)
# Equation (5-15), Boulanger et al. (2003)
k_SIunits = cSigma * k_SIunits

# define y50 based on pult and subgrade modulus k

# based on API (1987) recommendations, p-y curves are described using tanh functions.
# tcl does not have the atanh function, so must define this specifically

# i.e. atanh(x) = 1/2*ln((1+x)/(1-x)), |x| < 1

# when half of full resistance has been mobilized, p(y50)/pult = 0.5
x = 0.5
log_1 = math.log((1+x)/(1-x))
atanh_value = 0.5 * log_1

# need to be careful at ground surface (don't want to divide by zero)
if pyDepth == 0.0:
pyDepth = 0.01

y50 = 0.5 * (pu/ A)/(k_SIunits * pyDepth) * atanh_value
# return pult and y50 parameters
outResult = []
outResult.append(pult)
outResult.append(y50)

return outResult

#########################################################################################################################################################################

#########################################################################################################################################################################

###########################################################
# #
# Procedure to compute ultimate tip resistance, qult, and #
# displacement at 50% mobilization of qult, z50, for #
# use in q-z curves for cohesionless soil. #
# Converted to openseespy by: Pavan Chigullapally #
# University of Auckland #
# Created by: Chris McGann #
# Pedro Arduino #
# University of Washington #
# #
###########################################################

# references
# Meyerhof G.G. (1976). "Bearing capacity and settlement of pile foundations."
# J. Geotech. Eng. Div., ASCE, 102(3), 195-228.
#
# Vijayvergiya, V.N. (1977). "Load-movement characteristics of piles."
# Proc., Ports 77 Conf., ASCE, New York.
#
# Kulhawy, F.H. ad Mayne, P.W. (1990). Manual on Estimating Soil Properties for
# Foundation Design. Electrical Power Research Institute. EPRI EL-6800,
# Project 1493-6 Final Report.

def get_qzParam (phiDegree, b, sigV, G):

# define required constants; pi, atmospheric pressure (kPa), pa, and coeff. of lat earth pressure, Ko
pi = 3.14159265358979
pa = 101
sin_4 = math.sin(phiDegree * (pi/180))
Ko = 1 - sin_4

# ultimate tip pressure can be computed by qult = Nq*sigV after Meyerhof (1976)
# where Nq is a bearing capacity factor, phi is friction angle, and sigV is eff. overburden
# stress at the pile tip.
phi = phiDegree * (pi/180)

# rigidity index
tan_7 = math.tan(phi)
Ir = G/(sigV * tan_7)
# bearing capacity factor
tan_8 = math.tan(pi/4+phi/2)
sin_5 = math.sin(phi)
pow_4 = math.pow(tan_8,2)
pow_5 = math.pow(Ir,(4*sin_5)/(3*(1+sin_5)))
exp_4 = math.exp(pi/2-phi)

Nq = (1+2*Ko)*(1/(3-sin_5))*exp_4*(pow_4)*(pow_5)
# tip resistance
qu = Nq * sigV
# QzSimple1 material formulated with qult as force, not stress, multiply by area of pile tip
pow_6 = math.pow(b, 2)
qult = qu * pi*pow_6/4

# the q-z curve of Vijayvergiya (1977) has the form, q(z) = qult*(z/zc)^(1/3)
# where zc is critical tip deflection given as ranging from 3-9% of the
# pile diameter at the tip.

# assume zc is 5% of pile diameter
zc = 0.05 * b

# based on Vijayvergiya (1977) curve, z50 = 0.125*zc
z50 = 0.125 * zc

# return values of qult and z50 for use in q-z material
outResult = []
outResult.append(qult)
outResult.append(z50)

return outResult

#########################################################################################################################################################################

#########################################################################################################################################################################
##########################################################
# #
# Procedure to compute ultimate resistance, tult, and #
# displacement at 50% mobilization of tult, z50, for #
# use in t-z curves for cohesionless soil. #
# Converted to openseespy by: Pavan Chigullapally #
# University of Auckland #
# Created by: Chris McGann #
# University of Washington #
# #
###########################################################

def get_tzParam ( phi, b, sigV, pEleLength):

# references
# Mosher, R.L. (1984). "Load transfer criteria for numerical analysis of
# axial loaded piles in sand." U.S. Army Engineering and Waterways
# Experimental Station, Automatic Data Processing Center, Vicksburg, Miss.
#
# Kulhawy, F.H. (1991). "Drilled shaft foundations." Foundation engineering
# handbook, 2nd Ed., Chap 14, H.-Y. Fang ed., Van Nostrand Reinhold, New York

pi = 3.14159265358979

# Compute tult based on tult = Ko*sigV*pi*dia*tan(delta), where
# Ko is coeff. of lateral earth pressure at rest,
# taken as Ko = 0.4
# delta is interface friction between soil and pile,
# taken as delta = 0.8*phi to be representative of a
# smooth precast concrete pile after Kulhawy (1991)

delta = 0.8 * phi * pi/180

# if z = 0 (ground surface) need to specify a small non-zero value of sigV

if sigV == 0.0:
sigV = 0.01

tan_9 = math.tan(delta)
tu = 0.4 * sigV * pi * b * tan_9

# TzSimple1 material formulated with tult as force, not stress, multiply by tributary length of pile
tult = tu * pEleLength

# Mosher (1984) provides recommended initial tangents based on friction angle
# values are in units of psf/in
kf = [6000, 10000, 10000, 14000, 14000, 18000]
fric = [28, 31, 32, 34, 35, 38]

dataNum = len(fric)


# determine kf for input value of phi, linear interpolation for intermediate values
if phi < fric[0]:
k = kf[0]
elif phi > fric[5]:
k = kf[5]
else:
for i in range(dataNum):
if fric[i] <= phi and phi <= fric[i+1]:
k = ((kf[i+1] - kf[i])/(fric[i+1] - fric[i])) * (phi - fric[i]) + kf[i]


# need to convert kf to units of kN/m^3
kSIunits = k * 1.885

# based on a t-z curve of the shape recommended by Mosher (1984), z50 = tult/kf
z50 = tult / kSIunits

# return values of tult and z50 for use in t-z material
outResult = []
outResult.append(tult)
outResult.append(z50)

return outResult

###################################################################################################################
# 获得土弹簧参数代码结束，在下方调用 #
###################################################################################################################

print("UW OKKK")
for i in range(PileNodeNum):
if(i*EleSize>L1):
pyDepth = i*EleSize-L1
pyParam = get_pyParam(pyDepth,gamma,phi,D1,EleSize,puSwitch,kSwitch,gwtSwitch)
pult = pyParam[0]
y50 = pyParam[1]
uniaxialMaterial('PySimple1',i+1101,2,pult,y50,0.0) #x方向的py
uniaxialMaterial('PySimple1',i+2101,2,pult,y50,0.0)
uniaxialMaterial('PySimple1',i+3101,2,pult,y50,0.0)
uniaxialMaterial('PySimple1',i+4101,2,pult,y50,0.0)



#t-z土弹簧,q-z土弹簧
for i in range(PileNodeNum):
if(i*EleSize>L1 and i*EleSize<L1+L2):
pyDepth = i*EleSize-L1
sigV = gamma * pyDepth
tzParam = get_tzParam(phi,D1,sigV,EleSize)
tult = pyParam[0]
z50 = pyParam[1]
uniaxialMaterial('TzSimple1',i+1301,2,tult,z50,0.0)
uniaxialMaterial('TzSimple1',i+2301,2,tult,z50,0.0)
uniaxialMaterial('TzSimple1',i+3301,2,tult,z50,0.0)
uniaxialMaterial('TzSimple1',i+4301,2,tult,z50,0.0)
if(i*EleSize==L1+L2):
sigVq = gamma * L2
qzParam = get_qzParam(phi,D1,sigVq,Gsoil)
qult = qzParam [0]
z50q = qzParam [1]
uniaxialMaterial('QzSimple1',i+1301,2,qult,z50q,0.0)
uniaxialMaterial('QzSimple1',i+2301,2,qult,z50q,0.0)
uniaxialMaterial('QzSimple1',i+3301,2,qult,z50q,0.0)
uniaxialMaterial('QzSimple1',i+4301,2,qult,z50q,0.0)


print("Finished creating all p-y, t-z, and q-z spring material objects...")

######################土弹簧material建立成功

#--------------------建立土弹簧Zero-Length Elements
#???问题：x向和y向py弹簧，通过主从关系连接到桩上，z向弹簧，通过zero length 和x，y向弹簧连接？（感觉可能不收敛？？,目前选的这个
#???还是，其中一个方向，例如x向与z向连，通过，然后equal DOF与桩；然后pile上的点和y向zerolength连。
for j in range(1):
for i in range(PileNodeNum):
if(i*EleSize>L1):
element('zeroLength',(j+1)*10000+i+1101,(j+1)*10000+i+1101,(j+1)*10000+i+1301,'-mat',i+1101,i+1101,i+1301,'-dir',1,2,3)
#element('zeroLength',(j+1)*10000+i+1201,(j+1)*10000+i+1201,(j+1)*10000+i+1301,'-mat',i+1101,i+1301,'-dir',2,3)
element('zeroLength',(j+1)*10000+i+2101,(j+1)*10000+i+2101,(j+1)*10000+i+2301,'-mat',i+2101,i+2101,i+2301,'-dir',1,2,3)
#element('zeroLength',(j+1)*10000+i+2201,(j+1)*10000+i+2201,(j+1)*10000+i+2301,'-mat',i+2101,i+2301,'-dir',2,3)
element('zeroLength',(j+1)*10000+i+3101,(j+1)*10000+i+3101,(j+1)*10000+i+3301,'-mat',i+3101,i+3101,i+3301,'-dir',1,2,3)
#element('zeroLength',(j+1)*10000+i+3201,(j+1)*10000+i+3201,(j+1)*10000+i+3301,'-mat',i+3101,i+3301,'-dir',2,3)
element('zeroLength',(j+1)*10000+i+4101,(j+1)*10000+i+4101,(j+1)*10000+i+4301,'-mat',i+4101,i+4101,i+4301,'-dir',1,2,3)
#element('zeroLength',(j+1)*10000+i+4201,(j+1)*10000+i+4201,(j+1)*10000+i+4301,'-mat',i+4101,i+4301,'-dir',2,3)


print("Finished creating all zero-Length elements for springs...")

#建立pile nodes
model('basic', '-ndm', 3, '-ndf', 6)
#for i in range(4):
# node((i+1)*10000+7,i*25+0.5*pcwidth,0.5*pclength,z4)
# node((i+1)*10000+8,i*25+0.5*pcwidth,-0.5*pclength,z4)
# node((i+1)*10000+9,i*25-0.5*pcwidth,-0.5*pclength,z4)
# node((i+1)*10000+10,i*25-0.5*pcwidth,0.5*pclength,z4)

PileMass = 276.32 #t
PileNodeMass=[PileMass/PileNodeNum,PileMass/PileNodeNum,PileMass/PileNodeNum,0,0,0]
for j in range(1):
for i in range(PileNodeNum):
zCoord = -i*EleSize
node((j+1)*10000+i+101,25*j+0.5*XLength,0.5*YLength,zCoord)
node((j+1)*10000+i+201,25*j+0.5*XLength,-0.5*YLength,zCoord)
node((j+1)*10000+i+301,25*j-0.5*XLength,-0.5*YLength,zCoord)
node((j+1)*10000+i+401,25*j-0.5*XLength,0.5*YLength,zCoord)
fix((j+1)*10000+i+101,0,0,0,0,0,1)
fix((j+1)*10000+i+201,0,0,0,0,0,1)
fix((j+1)*10000+i+301,0,0,0,0,0,1)
fix((j+1)*10000+i+401,0,0,0,0,0,1)
mass((j+1)*10000+i+101,*PileNodeMass)
mass((j+1)*10000+i+201,*PileNodeMass)
mass((j+1)*10000+i+301,*PileNodeMass)
mass((j+1)*10000+i+401,*PileNodeMass)


print("Finished creating all pile nodes & fixities...")

#建立土弹簧和Pile 之间的Equal Dof 主从约束
for j in range(1):
for i in range (PileNodeNum):
if(i*EleSize>L1):
equalDOF((j+1)*10000+i+101,(j+1)*10000+i+1101,1,2,3)
#equalDOF((j+1)*10000+i+101,(j+1)*10000+i+1201,1,2,3)
equalDOF((j+1)*10000+i+201,(j+1)*10000+i+2101,1,2,3)
#equalDOF((j+1)*10000+i+201,(j+1)*10000+i+2201,1,2,3)
equalDOF((j+1)*10000+i+301,(j+1)*10000+i+3101,1,2,3)
#equalDOF((j+1)*10000+i+301,(j+1)*10000+i+3201,1,2,3)
equalDOF((j+1)*10000+i+401,(j+1)*10000+i+4101,1,2,3)
#equalDOF((j+1)*10000+i+401,(j+1)*10000+i+4201,1,2,3)


print("Finished creating all equal degrees of freedom...")
#-------------Beginning of Pile Element setup
#Define materials for nonlinear columns
#Pile Cross section:Cover Concrete, Core Concrete, Steel Bars


secTag = 2
E = 25000000.0
A = 0.785
Iz = 0.049
Iy = 0.049
G = 9615385.0
J = 0.098

matTag = 10
section('Elastic', 2, E, A, Iz, Iy, G, J)

# elastic torsional material for combined 3D section
uniaxialMaterial('Elastic', 10, 1e10)

# create combined 3D section
secTag3D = 3
op.section('Aggregator', secTag3D, 10, 'T', '-section', 2)

#局部坐标-全局坐标转换
transfType2 = 'Linear'
transfTag2 = 5
vecxz2=[0.0,1.0,1.0] #与全局X的叉乘是局部的y

geomTransf(transfType2,transfTag2,*vecxz2)
beamIntegration('Lobatto',2,2,2) #积分：Lobatto，1,1号截面，5个积分点
for j in range(1):
for i in range(PileEleNum):
element('dispBeamColumn',(j+1)*10000+i+101,(j+1)*10000+i+101,(j+1)*10000+i+102,transfTag2,2,'-mass',PileMass/PileEleNum)
element('dispBeamColumn',(j+1)*10000+i+201,(j+1)*10000+i+201,(j+1)*10000+i+202,transfTag2,2,'-mass',PileMass/PileEleNum)
element('dispBeamColumn',(j+1)*10000+i+301,(j+1)*10000+i+301,(j+1)*10000+i+302,transfTag2,2,'-mass',PileMass/PileEleNum)
element('dispBeamColumn',(j+1)*10000+i+401,(j+1)*10000+i+401,(j+1)*10000+i+402,transfTag2,2,'-mass',PileMass/PileEleNum)


print('Pile setup finished')

#-------------承台与桩顶连接，可提离的模拟
#可考虑：水平约束，竖向约束（只压不拉），摩擦，碰撞
#目前考虑，水平约束，竖向约束，和摩擦

#--------------承台中心与四周连接刚臂
vecxz3=[0.5*pcwidth,0.5*pclength,1.0]
vecxz4=[0.5*pcwidth,-0.5*pclength,1.0]
vecxz5=[-0.5*pcwidth,-0.5*pclength,1.0]
vecxz6=[-0.5*pcwidth,0.5*pclength,1.0]
geomTransf('Linear',20,*vecxz3)
geomTransf('Linear',21,*vecxz4)
geomTransf('Linear',22,*vecxz5)
geomTransf('Linear',23,*vecxz6)
for i in range(1):
element('elasticBeamColumn',(i+1)*10000+6,(i+1)*10000+5,(i+1)*10000+101,1400,32500000,13541667,2733,1429,1867,20)
element('elasticBeamColumn',(i+1)*10000+7,(i+1)*10000+5,(i+1)*10000+201,1400,32500000,13541667,2733,1429,1867,21)
element('elasticBeamColumn',(i+1)*10000+8,(i+1)*10000+5,(i+1)*10000+301,1400,32500000,13541667,2733,1429,1867,22)
element('elasticBeamColumn',(i+1)*10000+9,(i+1)*10000+5,(i+1)*10000+401,1400,32500000,13541667,2733,1429,1867,23)

print('承台中心与四周刚臂单元定义结束')


#GAP单元，摩擦单元，水平约束单元建立
#Pile Head 节点号：101,201,301,401；承台底节点，5,6,7,8




print('桩顶与承台连接定义完毕')
plot_model()


###定义recorder
recorder('Node','-file','shBnode101Disp.out','-time','-node',10001,'-dof',1,'disp')
recorder('Node','-file','shBnode103Disp.out','-time','-node',10003,'-dof',1,'disp')
recorder('Node','-file','shBnode101DispY.out','-time','-node',10001,'-dof',2,'disp')
recorder('Node','-file','shBnode103DispY.out','-time','-node',10003,'-dof',2,'disp')
recorder('Node','-file','SoilReaction.out','-time','-nodeRange',11103,PileNodeNum+11103-2,'-dof',1,2,3,'reaction')
#定义自重荷载：主梁，墩柱，承台，桩(KN)
gravitBeam1 = 3960
gravitBeam2 = 1440
gravityPier = 324
gravityPileCap = 960
gravityPile =PileMass*10
timeSeries('Linear',1)
pattern('Plain',1,1)
for i in range(4):
load(10000+i*100+101,0,0,(-gravitBeam2-gravityPier-gravityPileCap)*0.25,0,0,0)


for j in range(1):
for i in range (PileNodeNum):
load((j+1)*10000+i+101,0,0,-gravityPile/PileNodeNum,0,0,0)
load((j+1)*10000+i+201,0,0,-gravityPile/PileNodeNum,0,0,0)
load((j+1)*10000+i+301,0,0,-gravityPile/PileNodeNum,0,0,0)
load((j+1)*10000+i+401,0,0,-gravityPile/PileNodeNum,0,0,0)

constraints('Transformation') # how it handles boundary conditions
numberer('Plain') # renumber dof's to minimize band-width (optimization), if you want to
system('BandGeneral') # how to store and solve the system of equations in the analysis
algorithm('Linear') # use Linear algorithm for linear analysis
integrator('LoadControl', 0.1) # determine the next time step for an analysis, # apply gravity in 10 steps
analysis('Static') # define type of analysis static or transient
analyze(10) # perform gravity analysis
loadConst('-time', 0.0) # hold gravity constant and restart time

#设置地震荷载
G =9.8
dt,nPts=ReadRecord.ReadRecord('elCentro.at2','elCentro.dat')

timeSeries('Path',2,'-dt',dt,'-filePath','elCentro.dat','-factor',G)
pattern('MultipleSupport',2)
#pattern('UniformExcitation',2,1,'-accel',2)
groundMotion(1, 'Plain','-accel',1, '-int','Trapezoidal')
imposedMotion(10101,1,1)
imposedMotion(10201,1,1)
imposedMotion(10301,1,1)
imposedMotion(10401,1,1)
##

Tol =1*10**-8
maxNumIter = 10
wipeAnalysis()
constraints('Transformation')
numberer('RCM')
system('BandGeneral')
test('EnergyIncr', Tol, maxNumIter)
algorithm('ModifiedNewton')

NewmarkGamma = 0.5
NewmarkBeta = 0.25
integrator('Newmark', NewmarkGamma, NewmarkBeta)
analysis('Transient')


TmaxAnalysis = dt*nPts # maximum duration of ground-motion analysis

Nsteps = int(TmaxAnalysis/ dt)

ok = analyze(Nsteps, dt)
tCurrent = getTime()

test01 = {1:'NormDispIncr', 2: 'RelativeEnergyIncr', 4: 'RelativeNormUnbalance',5: 'RelativeNormDispIncr', 6: 'NormUnbalance'}
algorithm01 = {1:'KrylovNewton', 2: 'SecantNewton' , 4: 'RaphsonNewton',5: 'PeriodicNewton', 6: 'BFGS', 7: 'Broyden', 8: 'NewtonLineSearch'}

for i in test01:
for j in algorithm01:

if ok != 0:
if j < 4:
algorithm(algorithm01[j], '-initial')

else:
algorithm(algorithm01[j])

test(test01[i], Tol, 10)
ok = analyze(Nsteps, dt)
print(test01[i], algorithm01[j], ok)
if ok == 0:
break
else:
continue

reactions()
Nodereactions1 = dict()
Nodereactions2 = dict()
for i in range(PileEleNum):
Nodereactions1[i] = eleForce(i+10101,5)
Nodereactions2[i] = eleForce(i+10201,5)
print(Nodereactions1[i])

wipe()

[jennyanams]

I am very impressed with this trial software and do want to purchase it at the right price.
https://www.couponstechie.com/promo-codes/kaspersky