Tunnel lining in OpenSeespy

[mehrshadgh]

Hi, thanks for your excellent work, I have learned quite a lot from that. I would like to ask Professor Scott a question about modeling tunnel lining.
I read the thesis of Stephanie Lange "Numerical Response of Concrete-Lined Tunnels Crossing Active Faults". How can be modeled soil spring with a zero-length section? It is not clear. I tried to model with zero-length element and Elastic no tension material, but this did not work. You will be of great help to me if you give me guidance.
The code is here that used zero length element and elastic no tension:

import openseespy.opensees as ops
import openseespy.postprocessing.ops_vis as opsv
import matplotlib.pyplot as plt
from numpy import pi as pi
import numpy as np


from numpy import sinh as sinh
from numpy import cosh as cosh
from numpy import sin as sin
from numpy import cos as cos

def analytical_def(D_E, t, i, L, P):
E_g = 100e6
i = i
D_E= 6.6
t = t
D_I = D_E-2*t
k = (12*E_g/i) * (D_E/2)
E = 26700e6
I = (pi/64)*(D_E**4-D_I**4)
lan = (k/(4*E*I))**0.25
L = L
P = P
d = (
( P*lan/(2*k) ) *
(
( cosh(lan*L) + cos(lan*L) ) /
( sinh(lan*L) + sin(lan*L) )
)
)
return d
D_E = 6.6
t = 0.3
i = 10
L = 300
P = 100e6
d = analytical_def(D_E, t, i, L, P)
print('d = ', d*1000 , ' mm')


#########################################################
def My_model():
N_Node = 61
num_spring = N_Node - 2
#####################################################
# Build a mmodel
ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3) # frame 2D
# Define nodes
ops.node(1, 0.0, 0.0)
ops.node(2, 5.0, 0.0)
ops.node(3, 10.0, 0.0)
ops.node(4, 15.0, 0.0)
ops.node(5, 20.0, 0.0)
ops.node(6, 25.0, 0.0)
ops.node(7, 30.0, 0.0)
ops.node(8, 35.0, 0.0)
ops.node(9, 40.0, 0.0)
ops.node(10, 45.0, 0.0)
ops.node(11, 50.0, 0.0)
ops.node(12, 55.0, 0.0)
ops.node(13, 60.0, 0.0)
ops.node(14, 65.0, 0.0)
ops.node(15, 70.0, 0.0)
ops.node(16, 75.0, 0.0)
ops.node(17, 80.0, 0.0)
ops.node(18, 85.0, 0.0)
ops.node(19, 90.0, 0.0)
ops.node(20, 95.0, 0.0)
ops.node(21, 100.0, 0.0)
ops.node(22, 105.0, 0.0)
ops.node(23, 110.0, 0.0)
ops.node(24, 115.0, 0.0)
ops.node(25, 120.0, 0.0)
ops.node(26, 125.0, 0.0)
ops.node(27, 130.0, 0.0)
ops.node(28, 135.0, 0.0)
ops.node(29, 140.0, 0.0)
ops.node(30, 145.0, 0.0)
ops.node(31, 150.0, 0.0)
ops.node(32, 155.0, 0.0)
ops.node(33, 160.0, 0.0)
ops.node(34, 165.0, 0.0)
ops.node(35, 170.0, 0.0)
ops.node(36, 175.0, 0.0)
ops.node(37, 180.0, 0.0)
ops.node(38, 185.0, 0.0)
ops.node(39, 190.0, 0.0)
ops.node(40, 195.0, 0.0)
ops.node(41, 200.0, 0.0)
ops.node(42, 205.0, 0.0)
ops.node(43, 210.0, 0.0)
ops.node(44, 215.0, 0.0)
ops.node(45, 220.0, 0.0)
ops.node(46, 225.0, 0.0)
ops.node(47, 230.0, 0.0)
ops.node(48, 235.0, 0.0)
ops.node(49, 240.0, 0.0)
ops.node(50, 245.0, 0.0)
ops.node(51, 250.0, 0.0)
ops.node(52, 255.0, 0.0)
ops.node(53, 260.0, 0.0)
ops.node(54, 265.0, 0.0)
ops.node(55, 270.0, 0.0)
ops.node(56, 275.0, 0.0)
ops.node(57, 280.0, 0.0)
ops.node(58, 285.0, 0.0)
ops.node(59, 290.0, 0.0)
ops.node(60, 295.0, 0.0)
ops.node(61, 300.0, 0.0)
#
ops.node(121, 5.0, 0.0)
ops.node(122, 10.0, 0.0)
ops.node(123, 15.0, 0.0)
ops.node(124, 20.0, 0.0)
ops.node(125, 25.0, 0.0)
ops.node(126, 30.0, 0.0)
ops.node(127, 35.0, 0.0)
ops.node(128, 40.0, 0.0)
ops.node(129, 45.0, 0.0)
ops.node(130, 50.0, 0.0)
ops.node(131, 55.0, 0.0)
ops.node(132, 60.0, 0.0)
ops.node(133, 65.0, 0.0)
ops.node(134, 70.0, 0.0)
ops.node(135, 75.0, 0.0)
ops.node(136, 80.0, 0.0)
ops.node(137, 85.0, 0.0)
ops.node(138, 90.0, 0.0)
ops.node(139, 95.0, 0.0)
ops.node(140, 100.0, 0.0)
ops.node(141, 105.0, 0.0)
ops.node(142, 110.0, 0.0)
ops.node(143, 115.0, 0.0)
ops.node(144, 120.0, 0.0)
ops.node(145, 125.0, 0.0)
ops.node(146, 130.0, 0.0)
ops.node(147, 135.0, 0.0)
ops.node(148, 140.0, 0.0)
ops.node(149, 145.0, 0.0)
ops.node(150, 150.0, 0.0)
ops.node(151, 155.0, 0.0)
ops.node(152, 160.0, 0.0)
ops.node(153, 165.0, 0.0)
ops.node(154, 170.0, 0.0)
ops.node(155, 175.0, 0.0)
ops.node(156, 180.0, 0.0)
ops.node(157, 185.0, 0.0)
ops.node(158, 190.0, 0.0)
ops.node(159, 195.0, 0.0)
ops.node(160, 200.0, 0.0)
ops.node(161, 205.0, 0.0)
ops.node(162, 210.0, 0.0)
ops.node(163, 215.0, 0.0)
ops.node(164, 220.0, 0.0)
ops.node(165, 225.0, 0.0)
ops.node(166, 230.0, 0.0)
ops.node(167, 235.0, 0.0)
ops.node(168, 240.0, 0.0)
ops.node(169, 245.0, 0.0)
ops.node(170, 250.0, 0.0)
ops.node(171, 255.0, 0.0)
ops.node(172, 260.0, 0.0)
ops.node(173, 265.0, 0.0)
ops.node(174, 270.0, 0.0)
ops.node(175, 275.0, 0.0)
ops.node(176, 280.0, 0.0)
ops.node(177, 285.0, 0.0)
ops.node(178, 290.0, 0.0)
ops.node(179, 295.0, 0.0)
# Define start and end Constraint
ops.fix( 1, 1, 1, 1)
ops.fix( 61, 1, 1, 1)
#
ops.fix( 121, 1, 1, 1)
ops.fix( 122, 1, 1, 1)
ops.fix( 123, 1, 1, 1)
ops.fix( 124, 1, 1, 1)
ops.fix( 125, 1, 1, 1)
ops.fix( 126, 1, 1, 1)
ops.fix( 127, 1, 1, 1)
ops.fix( 128, 1, 1, 1)
ops.fix( 129, 1, 1, 1)
ops.fix( 130, 1, 1, 1)
ops.fix( 131, 1, 1, 1)
ops.fix( 132, 1, 1, 1)
ops.fix( 133, 1, 1, 1)
ops.fix( 134, 1, 1, 1)
ops.fix( 135, 1, 1, 1)
ops.fix( 136, 1, 1, 1)
ops.fix( 137, 1, 1, 1)
ops.fix( 138, 1, 1, 1)
ops.fix( 139, 1, 1, 1)
ops.fix( 140, 1, 1, 1)
ops.fix( 141, 1, 1, 1)
ops.fix( 142, 1, 1, 1)
ops.fix( 143, 1, 1, 1)
ops.fix( 144, 1, 1, 1)
ops.fix( 145, 1, 1, 1)
ops.fix( 146, 1, 1, 1)
ops.fix( 147, 1, 1, 1)
ops.fix( 148, 1, 1, 1)
ops.fix( 149, 1, 1, 1)
ops.fix( 150, 1, 1, 1)
ops.fix( 151, 1, 1, 1)
ops.fix( 152, 1, 1, 1)
ops.fix( 153, 1, 1, 1)
ops.fix( 154, 1, 1, 1)
ops.fix( 155, 1, 1, 1)
ops.fix( 156, 1, 1, 1)
ops.fix( 157, 1, 1, 1)
ops.fix( 158, 1, 1, 1)
ops.fix( 159, 1, 1, 1)
ops.fix( 160, 1, 1, 1)
ops.fix( 161, 1, 1, 1)
ops.fix( 162, 1, 1, 1)
ops.fix( 163, 1, 1, 1)
ops.fix( 164, 1, 1, 1)
ops.fix( 165, 1, 1, 1)
ops.fix( 166, 1, 1, 1)
ops.fix( 167, 1, 1, 1)
ops.fix( 168, 1, 1, 1)
ops.fix( 169, 1, 1, 1)
ops.fix( 170, 1, 1, 1)
ops.fix( 171, 1, 1, 1)
ops.fix( 172, 1, 1, 1)
ops.fix( 173, 1, 1, 1)
ops.fix( 174, 1, 1, 1)
ops.fix( 175, 1, 1, 1)
ops.fix( 176, 1, 1, 1)
ops.fix( 177, 1, 1, 1)
ops.fix( 178, 1, 1, 1)
ops.fix( 179, 1, 1, 1)
###
matTag_1 = 1
fpc = -24.132e6
epsc0 = -0.003
fpcu = -12.815e6
epsU = -0.007
Lambda = 0.15
ft = 2.855e6
Ets = 1427.5e6
#
matTag_2 = 2
Fy, E0 = 235.4e6, 200.1e9
bs = 0.01
Ro, CR1, CR2 = 18, 0.925, 0.15
a2, a4 = 1.0, 1.0
a1 = a2*Fy/E0
a3 = a4*Fy/E0
sigInit = 0.0
#ops.uniaxialMaterial('Concrete02', matTag_1, fpc, epsc0, fpcu, epsU, Lambda, ft, Ets)
#ops.uniaxialMaterial('Steel02' , matTag_2, Fy, E0, bs, *[Ro, CR1, CR2], a1, a2, a3, a4, 0.0)
#ops.uniaxialMaterial('Concrete01', matTag_1, fpc, epsc0, fpcu, epsU)
#ops.uniaxialMaterial('Steel01' , matTag_2, Fy, E0, bs)


###
Fiber_Tag_1 = 1
G = E0/(2*(1+0.2))
d_E = 6.60
d_I = 6.00
Ixx = pi * d_E**4/64 - pi * d_I**4/64
Iyy = pi * d_E**4/64 - pi * d_I**4/64
J = Ixx + Iyy
GJ = G * J
#nc1, nr1, nc2, nr2, nb1, a_beg2, a_end2 = 16, 3, 16, 2, 60, 0., 360.
#As30 = (pi * 30.0**2 /4 ) * 1e-6
#r1, r2, r3, r4, rb_1, rb_2 = 3.00, 3.10, 3.20, 3.30, 2.90, 3.10

E = 26700e6
matTag_1 = 1
ops.uniaxialMaterial('Elastic', matTag_1, E)
nc2, nr2 = 16, 4
r1 , r2 = 3.00, 3.30
#
matTag_2 = 2
Fy, E0 = 235.4e6, 200.1e9
bs = 0.01
ops.uniaxialMaterial('Steel01' , matTag_2, Fy, E0, bs)
#
nc2, nr2 = 16, 4
r1 , r2 = 3.00, 3.30
nb1, a_beg2, a_end2 = 60, 0., 360.
As30 = (pi * 30.0**2 /4 ) * 1e-6
rb_1, rb_2 = 3.10, 3.20
ops.section('Fiber', Fiber_Tag_1, '-GJ', GJ)
ops.patch('circ', matTag_1, nc2, nr2, 0.0, 0.0, r1, r2, 0.0, 360.0) #
ops.layer('circ', matTag_2, nb1, As30, 0.0, 0.0, rb_1, 0.0, 360.0) # Rebar Inside
ops.layer('circ', matTag_2, nb1, As30, 0.0, 0.0, rb_2, 0.0, 360.0) # Rebar External
###
N = 2
BI_tag = 1
ops.beamIntegration('Lobatto', BI_tag, Fiber_Tag_1, N)
transfTag = 1
ops.geomTransf('Linear', transfTag)
ops.element('forceBeamColumn', 1, *[1,2], 1, 1)
ops.element('forceBeamColumn', 2, *[2,3], 1, 1)
ops.element('forceBeamColumn', 3, *[3,4], 1, 1)
ops.element('forceBeamColumn', 4, *[4,5], 1, 1)
ops.element('forceBeamColumn', 5, *[5,6], 1, 1)
ops.element('forceBeamColumn', 6, *[6,7], 1, 1)
ops.element('forceBeamColumn', 7, *[7,8], 1, 1)
ops.element('forceBeamColumn', 8, *[8,9], 1, 1)
ops.element('forceBeamColumn', 9, *[9,10], 1, 1)
ops.element('forceBeamColumn', 10, *[10,11], 1, 1)
ops.element('forceBeamColumn', 11, *[11,12], 1, 1)
ops.element('forceBeamColumn', 12, *[12,13], 1, 1)
ops.element('forceBeamColumn', 13, *[13,14], 1, 1)
ops.element('forceBeamColumn', 14, *[14,15], 1, 1)
ops.element('forceBeamColumn', 15, *[15,16], 1, 1)
ops.element('forceBeamColumn', 16, *[16,17], 1, 1)
ops.element('forceBeamColumn', 17, *[17,18], 1, 1)
ops.element('forceBeamColumn', 18, *[18,19], 1, 1)
ops.element('forceBeamColumn', 19, *[19,20], 1, 1)
ops.element('forceBeamColumn', 20, *[20,21], 1, 1)
ops.element('forceBeamColumn', 21, *[21,22], 1, 1)
ops.element('forceBeamColumn', 22, *[22,23], 1, 1)
ops.element('forceBeamColumn', 23, *[23,24], 1, 1)
ops.element('forceBeamColumn', 24, *[24,25], 1, 1)
ops.element('forceBeamColumn', 25, *[25,26], 1, 1)
ops.element('forceBeamColumn', 26, *[26,27], 1, 1)
ops.element('forceBeamColumn', 27, *[27,28], 1, 1)
ops.element('forceBeamColumn', 28, *[28,29], 1, 1)
ops.element('forceBeamColumn', 29, *[29,30], 1, 1)
ops.element('forceBeamColumn', 30, *[30,31], 1, 1)
ops.element('forceBeamColumn', 31, *[31,32], 1, 1)
ops.element('forceBeamColumn', 32, *[32,33], 1, 1)
ops.element('forceBeamColumn', 33, *[33,34], 1, 1)
ops.element('forceBeamColumn', 34, *[34,35], 1, 1)
ops.element('forceBeamColumn', 35, *[35,36], 1, 1)
ops.element('forceBeamColumn', 36, *[36,37], 1, 1)
ops.element('forceBeamColumn', 37, *[37,38], 1, 1)
ops.element('forceBeamColumn', 38, *[38,39], 1, 1)
ops.element('forceBeamColumn', 39, *[39,40], 1, 1)
ops.element('forceBeamColumn', 40, *[40,41], 1, 1)
ops.element('forceBeamColumn', 41, *[41,42], 1, 1)
ops.element('forceBeamColumn', 42, *[42,43], 1, 1)
ops.element('forceBeamColumn', 43, *[43,44], 1, 1)
ops.element('forceBeamColumn', 44, *[44,45], 1, 1)
ops.element('forceBeamColumn', 45, *[45,46], 1, 1)
ops.element('forceBeamColumn', 46, *[46,47], 1, 1)
ops.element('forceBeamColumn', 47, *[47,48], 1, 1)
ops.element('forceBeamColumn', 48, *[48,49], 1, 1)
ops.element('forceBeamColumn', 49, *[49,50], 1, 1)
ops.element('forceBeamColumn', 50, *[50,51], 1, 1)
ops.element('forceBeamColumn', 51, *[51,52], 1, 1)
ops.element('forceBeamColumn', 52, *[52,53], 1, 1)
ops.element('forceBeamColumn', 53, *[53,54], 1, 1)
ops.element('forceBeamColumn', 54, *[54,55], 1, 1)
ops.element('forceBeamColumn', 55, *[55,56], 1, 1)
ops.element('forceBeamColumn', 56, *[56,57], 1, 1)
ops.element('forceBeamColumn', 57, *[57,58], 1, 1)
ops.element('forceBeamColumn', 58, *[58,59], 1, 1)
ops.element('forceBeamColumn', 59, *[59,60], 1, 1)
ops.element('forceBeamColumn', 60, *[60,61], 1, 1)
###
matTag_3 = 3
E = 100e6
ops.uniaxialMaterial('ENT', matTag_3, E)

ops.element('zeroLength', 121, *[2, 121], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 122, *[3, 122], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 123, *[4, 123], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 124, *[5, 124], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 125, *[6, 125], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 126, *[7, 126], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 127, *[8, 127], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 128, *[9, 128], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 129, *[10, 129], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 130, *[11, 130], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 131, *[12, 131], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 132, *[13, 132], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 133, *[14, 133], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 134, *[15, 134], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 135, *[16, 135], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 136, *[17, 136], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 137, *[18, 137], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 138, *[19, 138], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 139, *[20, 139], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 140, *[21, 140], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 141, *[22, 141], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 142, *[23, 142], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 143, *[24, 143], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 144, *[25, 144], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 145, *[26, 145], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 146, *[27, 146], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 147, *[28, 147], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 148, *[29, 148], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 149, *[30, 149], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 150, *[31, 150], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 151, *[32, 151], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 152, *[33, 152], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 153, *[34, 153], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 154, *[35, 154], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 155, *[36, 155], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 156, *[37, 156], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 157, *[38, 157], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 158, *[39, 158], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 159, *[40, 159], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 160, *[41, 160], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 161, *[42, 161], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 162, *[43, 162], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 163, *[44, 163], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 164, *[45, 164], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 165, *[46, 165], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 166, *[47, 166], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 167, *[48, 167], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 168, *[49, 168], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 169, *[50, 169], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 170, *[51, 170], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 171, *[52, 171], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 172, *[53, 172], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 173, *[54, 173], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 174, *[55, 174], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 175, *[56, 175], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 176, *[57, 176], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 177, *[58, 177], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 178, *[59, 178], '-mat', *[3,3], '-dir', *[1,2])
ops.element('zeroLength', 179, *[60, 179], '-mat', *[3,3], '-dir', *[1,2])
# Define constant axial load
axialLoad = -100.0e6 # KN
ops.timeSeries('Constant', 1)
ops.pattern('Plain', 1, 1)
ops.load(int( (N_Node+1)/2 ), 0.0, axialLoad, 0.0)
# Static Analysis
ops.wipeAnalysis()
ops.numberer("Plain")
ops.constraints("Plain")
ops.algorithm("Linear")
ops.system("ProfileSPD")
ops.integrator("LoadControl", 100)
ops.analysis("Static")
ops.analyze(100)
###

disp = ops.nodeDisp(int( (N_Node+1)/2 ),2)
return num_spring, disp


out_1, out_2 = My_model()

print('d = ', out_2*1000 , ' mm')

[mhscott]

You should not use the 'Linear' algorithm for nonlinear problems. Use 'Newton'.