My statics Model Runs, but when I try to apply a thermal it isn't working, anyone knows how to solve it?
import openseespy.opensees as op
import vfo.vfo as vfo
import math
import numpy as np
import ArcCircle as Arc
import Soil as sl
import time
import opsvis as opsv
import matplotlib.pyplot as plt
# ---------------
# Units
# ---------------
m = 1
cm = 10**-2*m
mm = 10**-3*m
N = 1
KN = 10e2*N
Pa = 1 *N/m**2
MPa = 10**6*Pa
GPa = 10**9*Pa
inches = 0.0254*m
ft = 12*inches
kip = 4.45*10**3*N
ksi = 6.89*10**6*Pa
# ---------------
# Input Variables
# ---------------
G = 22*KN/m**3
Gs = 28*KN/m**3
Gw = 9.81*KN/m**3
# Look at the soil profile above ARC
MaxElev = 461.0
# It's the x location from when the arc starts and afterwards you have no increase in soil elevation anymore
flatzone = 5.4
# Elevation of Arc: note not center that is the surface
ArcMaxElev= 451.65
# Start Elavation of concrete Arc
SElevArc = 448
# Water depth from the surface
waterd = 3.65
# Horizontal Distance between nodes for the arc in meters
ARCNodedis = 0.9 * m
# Wall Mesh
Wallmesh = 3
# ---------------
# Opensees Anaylsis
# ---------------
def GetSections():
# Remove existing model
op.wipe()
# set modelbuilder
op.model('basic', '-ndm', 2, '-ndf', 3)
# Mat Tags for steel elements
# Steel material
SteelmatTag = 1
E = 210 * GPa
b = 0.01
# Define materials
op.uniaxialMaterial("Steel01Thermal", SteelmatTag, E, E, b)
# Mat Tags for reinforced concrete elements
CoreConcmatTag = 2
CoverConcmatTag = 3
SteelReinmatTag = 4
# Concrete parameters
fpc = -30 * MPa
fpcu = 0.0 # Recommended to be set as 0
epsc = -0.002
ecu = -0.005
# Core concrete (confined)
# uniaxialMaterial('Concrete01', matTag, fpc, epsc0, fpcu, epsU)
op.uniaxialMaterial('Concrete01', CoreConcmatTag, fpc, epsc, fpcu, ecu)
# Cover concrete (unconfined)
# uniaxialMaterial('Concrete01', matTag, fpc, epsc0, fpcu, epsU)
op.uniaxialMaterial('Concrete01', CoverConcmatTag, fpc, epsc, fpcu, ecu)
# Steel parameters for reinforcing steel
fy = 420 * MPa # Yield stress
Er = 200 * GPa # Young's modulus
# uniaxialMaterial('Steel01', matTag, fy, E0, b)
op.uniaxialMaterial('Steel01Thermal', SteelReinmatTag, fy, Er, 0.01)
# Sections
# HEB 300 section
HEBsecTag = 1
HEBd = 0.3 * m
HEBtw = 0.011 * m
HEBbf = 0.3 * m
HEBtf = 0.019 * m
# section('WFSection2d', secTag, matTag, d, tw, bf, tf, Nfw, Nff)
op.section('WFSection2d', HEBsecTag, SteelmatTag, HEBd, HEBtw, HEBbf, HEBtf, 15, 15)
# IPE 180 section
IPEsecTag = 2
IPEd = 180 * mm
IPEtw = 5.3 * mm
IPEbf = 91 * mm
IPEtf = 8 * mm
# section('WFSection2d', secTag, matTag, d, tw, bf, tf, Nfw, Nff)
op.section('WFSection2d', IPEsecTag, SteelmatTag, IPEd, IPEtw, IPEbf, IPEtf, 15, 15)
# Concrete Wall sections
LWallsecTag = 3
SWallsecTag = 4
ArcsecTag = 5
# Dimensions of wall that will be made as a column
LcolWidth = 1.1 * m # large width
ScolWidth = 0.9 * m # smaller width
ArcWidth = 0.65 * m # Arc Width
colDepth = 2.0 * m
cover = 4.0 * cm
# Fibers
ver_conc_fib = 10 # Core Vertical Conc Fibers
hor_conc_fib = 20# Core Horizontal Conc Fibers
# Area of steel for bar phi 18 mm
As = math.pi * (18 * mm) ** 2 / 4
# Large Concrete Wall section
y1 = LcolWidth / 2.0
z1 = colDepth / 2.0
# Number of bars per layer
num_layer_top = 14
num_layer_mid = 2
num_layer_bot = 14
op.section('FiberThermal', LWallsecTag)
# Create the concrete core fibers
# patch('rect', matTag, numSubdivIJ, numSubdivJK, *crdsI, *crdsJ, *crdsK, *crdsL)
op.patch('rect', CoreConcmatTag, ver_conc_fib, hor_conc_fib, cover - y1, cover - z1, y1 - cover, z1 - cover)
# Create the concrete cover fibers (top, bottom, left, right)
op.patch('rect', CoverConcmatTag, ver_conc_fib, 1, -y1, z1 - cover, y1, z1)
op.patch('rect', CoverConcmatTag, ver_conc_fib, 1, -y1, -z1, y1, cover - z1)
op.patch('rect', CoverConcmatTag, 1, hor_conc_fib, -y1, cover - z1, cover - y1, z1 - cover)
op.patch('rect', CoverConcmatTag, 1, hor_conc_fib, y1 - cover, cover - z1, y1, z1 - cover)
# Create the reinforcing fibers (left, middle, right)
# layer('straight', matTag, numFiber, areaFiber, *start, *end)
op.layer('straight', SteelReinmatTag, num_layer_top, As, y1 - cover, z1 - cover, y1 - cover, cover - z1)
op.layer('straight', SteelReinmatTag, num_layer_mid, As, 0.0, z1 - cover, 0.0, cover - z1)
op.layer('straight', SteelReinmatTag, num_layer_bot, As, cover - y1, z1 - cover, cover - y1, cover - z1)
# Smaller Concrete Wall section
y1 = ScolWidth / 2.0
z1 = colDepth / 2.0
# Number of bars per layer
num_layer_top = 14
num_layer_mid = 2
num_layer_bot = 14
op.section('FiberThermal', SWallsecTag)
# Create the concrete core fibers
# patch('rect', matTag, numSubdivIJ, numSubdivJK, *crdsI, *crdsJ, *crdsK, *crdsL)
op.patch('rect', CoreConcmatTag, ver_conc_fib, hor_conc_fib, cover - y1, cover - z1, y1 - cover, z1 - cover)
# Create the concrete cover fibers (top, bottom, left, right)
op.patch('rect', CoverConcmatTag, ver_conc_fib, 1, -y1, z1 - cover, y1, z1)
op.patch('rect', CoverConcmatTag, ver_conc_fib, 1, -y1, -z1, y1, cover - z1)
op.patch('rect', CoverConcmatTag, 1, hor_conc_fib, -y1, cover - z1, cover - y1, z1 - cover)
op.patch('rect', CoverConcmatTag, 1, hor_conc_fib, y1 - cover, cover - z1, y1, z1 - cover)
# Create the reinforcing fibers (left, middle, right)
# layer('straight', matTag, numFiber, areaFiber, *start, *end)
op.layer('straight', SteelReinmatTag, num_layer_top, As, y1 - cover, z1 - cover, y1 - cover, cover - z1)
op.layer('straight', SteelReinmatTag, num_layer_mid, As, 0.0, z1 - cover, 0.0, cover - z1)
op.layer('straight', SteelReinmatTag, num_layer_bot, As, cover - y1, z1 - cover, cover - y1, cover - z1)
# Arc Concrete Wall section
y1 = ArcWidth / 2.0
z1 = colDepth / 2.0
# Number of bars per layer
num_layer_top = 14
num_layer_mid = 2
num_layer_bot = 14
op.section('FiberThermal', ArcsecTag)
# Create the concrete core fibers
# patch('rect', matTag, numSubdivIJ, numSubdivJK, *crdsI, *crdsJ, *crdsK, *crdsL)
op.patch('rect', CoreConcmatTag, ver_conc_fib, hor_conc_fib, cover - y1, cover - z1, y1 - cover, z1 - cover)
# Create the concrete cover fibers (top, bottom, left, right)
op.patch('rect', CoverConcmatTag, ver_conc_fib, 1, -y1, z1 - cover, y1, z1)
op.patch('rect', CoverConcmatTag, ver_conc_fib, 1, -y1, -z1, y1, cover - z1)
op.patch('rect', CoverConcmatTag, 1, hor_conc_fib, -y1, cover - z1, cover - y1, z1 - cover)
op.patch('rect', CoverConcmatTag, 1, hor_conc_fib, y1 - cover, cover - z1, y1, z1 - cover)
# Create the reinforcing fibers (left, middle, right)
# layer('straight', matTag, numFiber, areaFiber, *start, *end)
op.layer('straight', SteelReinmatTag, num_layer_top, As, y1 - cover, z1 - cover, y1 - cover, cover - z1)
op.layer('straight', SteelReinmatTag, num_layer_mid, As, 0.0, z1 - cover, 0.0, cover - z1)
op.layer('straight', SteelReinmatTag, num_layer_bot, As, cover - y1, z1 - cover, cover - y1, cover - z1)
# op.patch('quad') is needed for arc at 1 h 07 in vid
return HEBsecTag, IPEsecTag, LWallsecTag, SWallsecTag, ArcsecTag
def GetModel(HEBsecTag, IPEsecTag, LWallsecTag, SWallsecTag, ArcsecTag, Wallmesh):
# Node Locations
# Left wall th = 1.1
x1 = 0.
y1 = 0.
xx1 = 0.
yy1 = 4.58 * m
# Right Wall th = 1.1
x2 = 17.1 * m
y2 = 0.
xx2 = x2 * m
yy2 = yy1 * m
# Left Well th = 0.9
x3 = 0.
y3 = 7. * m
pt1 = (x3, y3)
# Right Wall th = 0.9
x4 = x2 * m
y4 = y3 * m
pt3 = (x4, y4)
# Center of Arc
x5 = x2/2 * m
y5 = 10.275 * m
pt2 = (x5, y5)
center, radius = Arc.circle_center_and_radius(pt1, pt2, pt3)
# create nodes
# LSnode = Left start node
LSnode = 10
# LLnode = Left Largewall node
LLnode = 20
# Lswnode = Left smallwall node
Lswnode = 30
# RSnode = Right start node
RSnode = 40
# RLnode = Right Largewall node
RLnode = 50
# Rswnode = Right smallwall node
Rswnode = 60
# Cnode = Connection node of tie rods
Cnode = 70
# Start Node for Arc
SArcNV = 71
# Create a list to store the values, now named 'node_array'
node_array = [LSnode, LLnode, Lswnode, RSnode, RLnode, Rswnode, Cnode, SArcNV]
# Left Wall nodes
op.node(LSnode, x1, y1)
op.node(LLnode,xx1,yy1)
op.node(Lswnode, x3, y3)
# Right wall nodes
op.node(RSnode, x2, y2)
op.node(RLnode,xx2,yy2)
op.node(Rswnode, x4, y4)
# Center node of Steel Sections of IPE and HEB
op.node(Cnode, x5, y3)
# Number of elements for walls
Wallmesh = Wallmesh
i = 1
while i < Wallmesh:
yn1 = (yy1/Wallmesh) * i
yn2 = yy1 + ((y3-yy1)/Wallmesh) * i
op.node(LSnode + i,0., yn1)
op.node(LLnode + i,0., yn2)
op.node(RSnode + i,x2, yn1)
op.node(RLnode + i,x2, yn2)
i += 1
# nodes for arc
Nodevalue = SArcNV
# Middle Node for Arc
MArcNV = 0
# Max value multiplied by 10 to better keep increments
ARCNodeit = int(x2*10)
fARCNodedis = int(ARCNodedis*10)
flag = False
for i in range(fARCNodedis, ARCNodeit, fARCNodedis):
inc = i/10
ArcY = Arc.find_y_on_circle(inc, center, radius)
op.node(Nodevalue, inc, ArcY)
Nodevalue += 1
if abs(x5 - inc) < ARCNodedis and flag == False:
op.node(Nodevalue, x5, y5)
MArcNV = Nodevalue
Nodevalue += 1
flag = True
# set boundary condition
op.fix(LSnode, 1, 1, 1)
op.fix(RSnode, 1, 1, 1)
# Beam integrator
HEBIntTag = 1
IPEIntTag = 2
LColIntTag = 3
SColIntTag = 4
ArcIntTag = 5
op.beamIntegration('Lobatto', HEBIntTag, HEBsecTag, 5)
op.beamIntegration('Lobatto', IPEIntTag, IPEsecTag, 5)
op.beamIntegration('Lobatto', LColIntTag, LWallsecTag, 10)
op.beamIntegration('Lobatto', SColIntTag, SWallsecTag, 10)
op.beamIntegration('Lobatto', ArcIntTag, ArcsecTag, 10)
# Geometric Transform
HEBTranfTag = 1
IPETranfTag = 2
LColTranfTag = 3
SColTranfTag = 4
ArcTranfTag = 5
op.geomTransf('Linear', HEBTranfTag)
op.geomTransf('Linear', IPETranfTag)
op.geomTransf('Linear', LColTranfTag)
op.geomTransf('Linear', SColTranfTag)
op.geomTransf('Linear', ArcTranfTag)
# Large Col.
# Assign starting element value
i = 1
ii = Wallmesh
for k in range(len(node_array)):
j = node_array[k]
if j == Lswnode or j >= Rswnode:
continue
if j < LLnode or j >= RSnode and j < RLnode:
for i in range (i ,ii, 1):
# op.element('Truss', eleTag, *eleNodes, A, matTag, <'-rho', rho>, <'-cMass', cFlag>, <'-doRayleigh', rFlag>)
op.element('forceBeamColumn', i, j, j + 1, LColTranfTag, LColIntTag, '-mass', 0.0)
i += 1
j += 1
op.element('forceBeamColumn', i, j, node_array[k + 1], LColTranfTag, LColIntTag, '-mass', 0.0)
i += 1
ii += Wallmesh
else:
for i in range (i ,ii, 1):
# op.element('Truss', eleTag, *eleNodes, A, matTag, <'-rho', rho>, <'-cMass', cFlag>, <'-doRayleigh', rFlag>)
op.element('forceBeamColumn', i, j, j + 1, SColTranfTag, SColIntTag, '-mass', 0.0)
i += 1
j += 1
op.element('forceBeamColumn', i, j, node_array[k + 1], SColTranfTag, SColIntTag, '-mass', 0.0)
i += 1
ii += Wallmesh
# Last wall element number
Lwallele = i - 1
# op.element('forceBeamColumn', 1, LSnode, LLnode, LColTranfTag, LColIntTag, '-mass', 0.0)
# op.element('forceBeamColumn', 2, LLnode, Lswnode, SColTranfTag, SColIntTag, '-mass', 0.0)
# op.element('forceBeamColumn', 3, RSnode, RLnode, LColTranfTag, LColIntTag, '-mass', 0.0)
# op.element('forceBeamColumn', 4, RLnode, Rswnode, SColTranfTag, SColIntTag, '-mass', 0.0)
# Tie Rod
# op.element('Truss', eleTag, *eleNodes, A, matTag, <'-rho', rho>, <'-cMass', cFlag>, <'-doRayleigh', rFlag>)
op.element('dispBeamColumnThermal', i, Lswnode, Cnode, HEBTranfTag, HEBIntTag, '-mass', 0.0)
op.element('dispBeamColumnThermal', i + 1, Cnode, Rswnode, HEBTranfTag, HEBIntTag, '-mass', 0.0)
op.element('dispBeamColumnThermal', i + 2, Cnode, MArcNV, IPETranfTag, IPEIntTag, '-mass', 0.0)
# Arc
# Starting element value of Arc
SEleArc = i + 3
# op.element('Truss', eleTag, *eleNodes, A, matTag, <'-rho', rho>, <'-cMass', cFlag>, <'-doRayleigh', rFlag>)
op.element('forceBeamColumn', SEleArc, Lswnode, SArcNV, ArcIntTag, ArcTranfTag, '-mass', 0.0)
# Arc Nodes total value
if round(x2 % ARCNodedis, 10) == 0:
Antv = int((x2//ARCNodedis))
else:
Antv = int((x2//ARCNodedis) + 1)
elements = i + 4
# x,y = op.nodeCoord(72)
# print(x,y)
for i in range(1,Antv,1):
op.element('forceBeamColumn', elements, Cnode + i, SArcNV + i, ArcIntTag, ArcTranfTag, '-mass', 0.0)
elements += 1
op.element('forceBeamColumn', elements, Rswnode, Nodevalue - 1, ArcIntTag, ArcTranfTag, '-mass', 0.0)
# number of total elements
totelements = elements
return totelements, SEleArc, Wallmesh, ARCNodedis, y3, x5
def GetRecorders():
# Record Results
# op.recorder('Node', '-file', filename, '-xml', filename, '-binary', filename, '-tcp', inetAddress, port, '-precision', nSD=6, '-timeSeries', tsTag, '-time', '-dT', deltaT=0.0, '-closeOnWrite', '-node', *nodeTags=[], '-nodeRange', startNode, endNode, '-region', regionTag, '-dof', *dofs=[], respType)
op.recorder('Node', '-file', "NodeDispHEB.out",'-time','-node', 7, '-dof', 1, 2, 3, 'disp')
op.recorder('Node', '-file', "NodeDisp5.out",'-time','-node', 70, '-dof', 1, 2, 3, 'disp')
op.recorder('Node', '-file', "ReactionHEB.out",'-time','-node', 1, 3, '-dof', 1, 2, 3, 'reaction')
op.recorder('Element', '-file', "ElementsHEB.out",'-time','-ele', 13, 14, 15, 'forces')
op.recorder('Element', '-file', "Arc.out",'-time','-ele', 17, 'forces')
op.recorder('Node', '-file', "wall.out",'-time','-node', 30, '-dof', 1, 2, 3, 'disp')
#def load(totelements, SEleArc, Wallmesh, G, Gs):
def RunAnalysisStatic(totelements, SEleArc, Wallmesh, ARCNodedis, y3, x5, G, Gs):
# ------------------------------
# Start Plot
# ------------------------------
ModelName = '2D_Model'
LoadCaseName = 'Static'
vfo.createODB(ModelName, LoadCaseName)
# ------------------------------
# Apply Loads
# ------------------------------
# timer start : op.start()
# create TimeSeries
op.timeSeries("Linear", 1)
# create a plain load pattern
op.pattern("Plain", 1, 1)
# Create the nodal load - command: load nodeID xForce yForce
# op.load(nodeTag, *loadValues)
# ops.eleLoad('-ele',eleTag,'-type','beamUniform',wya,wxa,aOverL,bOverL,wyb,wxb)
# Apply Load on Arc only surchage no lateral load
for i in range(SEleArc, totelements, 1):
x = (i - SEleArc)*ARCNodedis
h = sl.Height_Soil_Conc(x, x5, MaxElev, flatzone, SElevArc) * m
S = sl.surcharge_soil(G, Gs, h, MaxElev, ArcMaxElev)*2 # N/m
op.eleLoad('-ele', i, '-type', '-beamUniform', -S, 0.0, 0.0, 1.0, -S, 0.0)
h = sl.Height_Soil_Conc(0, x5, MaxElev, flatzone, SElevArc) * m
S = sl.surcharge_soil(G, Gs, h, MaxElev, ArcMaxElev)*2 # N/m
op.eleLoad('-ele', totelements, '-type', '-beamUniform', S, 0.0, 0.0, 1.0, S, 0.0)
# op.eleLoad('-ele', 8, '-type', '-beamUniform', -S, 0.0, 0.0, 1.0, -S, 0.0)
# op.eleLoad('-ele', 9, '-type', '-beamUniform', S, 0.0, 0.0, 1.0, S, 0.0)
# S in KN/m^2
S = S/2
ka = 0.3
# At h = 0 m & h = 7.05 m
h0 = 0*m
h7 = y3
# LatPa in KN/m so multiply by 2 meter
LatPah0= sl.lat_soil_pres(G,Gs,ka,h0,S,Gw,waterd)*2
LatPah7= sl.lat_soil_pres(G,Gs,ka,h7,S,Gw,waterd)*2
# print(LatPah0,LatPah7)
# Negative Side
for i in range(1,Wallmesh + 1, 1):
# Negative Side
# Choosing Element
k,j = op.eleNodes(i)
# Taking the lower node Coord
x,y = op.nodeCoord(k)
# Taking the upper node Coord
xx, yy = op.nodeCoord(j)
# Inverting the depth for the interpolation
y = h7-y
yy = h7 - yy
# print(i,k,j,x,y,xx,yy,sl.lin_int(yy, h0, LatPah0, h7, LatPah7))
op.eleLoad('-ele', i, '-type', '-beamUniform', -sl.lin_int(y, h0, LatPah0, h7, LatPah7), 0.0, 0.0, 1.0, -sl.lin_int(yy, h0, LatPah0, h7, LatPah7), 0.0)
k,j = op.eleNodes(i + Wallmesh)
x,y = op.nodeCoord(k)
xx, yy = op.nodeCoord(j)
y = h7-y
yy = h7 - yy
op.eleLoad('-ele', i + Wallmesh, '-type', '-beamUniform', -sl.lin_int(y, h0, LatPah0, h7, LatPah7), 0.0, 0.0, 1.0, -sl.lin_int(yy, h0, LatPah0, h7, LatPah7), 0.0)
# Positive Side
k,j = op.eleNodes(i + Wallmesh * 2)
x,y = op.nodeCoord(k)
xx, yy = op.nodeCoord(j)
y = h7-y
yy = h7 - yy
op.eleLoad('-ele', i + Wallmesh * 2, '-type', '-beamUniform', sl.lin_int(y, h0, LatPah0, h7, LatPah7), 0.0, 0.0, 1.0, sl.lin_int(yy, h0, LatPah0, h7, LatPah7), 0.0)
k,j = op.eleNodes(i + Wallmesh * 3)
x,y = op.nodeCoord(k)
xx, yy = op.nodeCoord(j)
y = h7-y
yy = h7 - yy
op.eleLoad('-ele', i + Wallmesh * 3, '-type', '-beamUniform', sl.lin_int(y, h0, LatPah0, h7, LatPah7), 0.0, 0.0, 1.0, sl.lin_int(yy, h0, LatPah0, h7, LatPah7), 0.0)
# ------------------------------
# Start of analysis generation
# ------------------------------
# create SOE
op.system("BandSPD")
# create DOF number
op.numberer("RCM")
# create constraint handler
op.constraints("Plain")
# create integrator
op.integrator("LoadControl", 0.05)
# create algorithm
op.algorithm("Newton")
# create analysis object
op.analysis("Static")
# perform the analysis
op.initialize()
ok = op.analyze(34)
# timer stop: op.stop()
vfo.plot_model(ModelName,show_nodetags="yes",show_eletags="yes")
vfo.plot_deformedshape(ModelName, LoadCaseName, scale=100, tstep=1, overlap='yes')
op.wipe()
def RunThermalFunction(totelements, SEleArc, Wallmesh, ARCNodedis, y3, x5, G, Gs):
# ------------------------------
# Start Plot
# ------------------------------
ModelName = '2D_Model'
LoadCaseName = 'Static'
vfo.createODB(ModelName, LoadCaseName)
# timer start : op.start()
#Timeseries Tag 1
tsTag = 1
# create TimeSeries
op.timeSeries("Linear", tsTag)
#Timeseries Tag 2
tsTag = 2
# create TimeSeries 2
op.timeSeries("Linear", tsTag)
# ------------------------------
# Apply Loads
# ------------------------------
# Pattern Tag
patternTag = 1
# create a plain load pattern
op.pattern("Plain", patternTag, tsTag)
# Create the nodal load - command: load nodeID xForce yForce
# op.load(nodeTag, *loadValues)
# ops.eleLoad('-ele',eleTag,'-type','beamUniform',wya,wxa,aOverL,bOverL,wyb,wxb)
# Apply Load on Arc only surchage no lateral load
for i in range(SEleArc, totelements, 1):
x = (i - SEleArc)*ARCNodedis
h = sl.Height_Soil_Conc(x, x5, MaxElev, flatzone, SElevArc) * m
S = sl.surcharge_soil(G, Gs, h, MaxElev, ArcMaxElev)*2 # N/m
op.eleLoad('-ele', i, '-type', '-beamUniform', -S, 0.0, 0.0, 1.0, -S, 0.0)
h = sl.Height_Soil_Conc(0, x5, MaxElev, flatzone, SElevArc) * m
S = sl.surcharge_soil(G, Gs, h, MaxElev, ArcMaxElev)*2 # N/m
op.eleLoad('-ele', totelements, '-type', '-beamUniform', S, 0.0, 0.0, 1.0, S, 0.0)
# op.eleLoad('-ele', 8, '-type', '-beamUniform', -S, 0.0, 0.0, 1.0, -S, 0.0)
# op.eleLoad('-ele', 9, '-type', '-beamUniform', S, 0.0, 0.0, 1.0, S, 0.0)
# S in KN/m^2
S = S/2
ka = 0.3
# At h = 0 m & h = 7.05 m
h0 = 0*m
h7 = y3
# LatPa in KN/m so multiply by 2 meter
LatPah0= sl.lat_soil_pres(G,Gs,ka,h0,S,Gw,waterd)*2
LatPah7= sl.lat_soil_pres(G,Gs,ka,h7,S,Gw,waterd)*2
# print(LatPah0,LatPah7)
# Negative Side
for i in range(1,Wallmesh + 1, 1):
# Negative Side
# Choosing Element
k,j = op.eleNodes(i)
# Taking the lower node Coord
x,y = op.nodeCoord(k)
# Taking the upper node Coord
xx, yy = op.nodeCoord(j)
# Inverting the depth for the interpolation
y = h7-y
yy = h7 - yy
# print(i,k,j,x,y,xx,yy,sl.lin_int(yy, h0, LatPah0, h7, LatPah7))
op.eleLoad('-ele', i, '-type', '-beamUniform', -sl.lin_int(y, h0, LatPah0, h7, LatPah7), 0.0, 0.0, 1.0, -sl.lin_int(yy, h0, LatPah0, h7, LatPah7), 0.0)
k,j = op.eleNodes(i + Wallmesh)
x,y = op.nodeCoord(k)
xx, yy = op.nodeCoord(j)
y = h7-y
yy = h7 - yy
op.eleLoad('-ele', i + Wallmesh, '-type', '-beamUniform', -sl.lin_int(y, h0, LatPah0, h7, LatPah7), 0.0, 0.0, 1.0, -sl.lin_int(yy, h0, LatPah0, h7, LatPah7), 0.0)
# Positive Side
k,j = op.eleNodes(i + Wallmesh * 2)
x,y = op.nodeCoord(k)
xx, yy = op.nodeCoord(j)
y = h7-y
yy = h7 - yy
op.eleLoad('-ele', i + Wallmesh * 2, '-type', '-beamUniform', sl.lin_int(y, h0, LatPah0, h7, LatPah7), 0.0, 0.0, 1.0, sl.lin_int(yy, h0, LatPah0, h7, LatPah7), 0.0)
k,j = op.eleNodes(i + Wallmesh * 3)
x,y = op.nodeCoord(k)
xx, yy = op.nodeCoord(j)
y = h7-y
yy = h7 - yy
op.eleLoad('-ele', i + Wallmesh * 3, '-type', '-beamUniform', sl.lin_int(y, h0, LatPah0, h7, LatPah7), 0.0, 0.0, 1.0, sl.lin_int(yy, h0, LatPah0, h7, LatPah7), 0.0)
# Pattern Tag
patternTag = 2
# MaxTemp
maxtemp = 1000.0
# create a plain load pattern
op.pattern("Plain", patternTag, tsTag)
# Apply Thermal Load
op.eleLoad('-ele', 14, '-type', '-beamThermal', 1000.0, -0.05, 1000.0, 0.05)
# ------------------------------
# Start of analysis generation
# ------------------------------
# Pattern 1
# create SOE
op.system("BandSPD")
# create DOF number
op.numberer("RCM")
# create constraint handler
op.constraints("Plain")
# create integrator
op.integrator("LoadControl", 0.05)
# create algorithm
op.algorithm("Newton")
# create analysis object
op.analysis("Static")
# perform the analysis
op.initialize()
if op.analyze(20)!= 0:
print("First analysis failed.")
# Remove LoadPattern 1
# You can either set the factor of TimeSeries 1 to 0 or leave it unchanged
op.loadConst('-pattern', 1, '-timeSeries', 1, '-fact', 1.0) # Continue with same factor or reduce if needed
# Set TimeSeries 2 to start applying load
op.loadConst('-pattern', 2, '-timeSeries', 2, '-fact', 1.0) # Now TimeSeries 2 is active
# Pattern 2
# Time Increment
incrtemp = 0.01
# create SOE
op.system("BandGeneral")
# create DOF number
op.numberer("Plain")
# Test
op.test('NormDispIncr', 1.0e-3, 100, 1)
# create constraint handler
op.constraints("Plain")
# create integrator
op.integrator("LoadControl", incrtemp)
# create algorithm
op.algorithm("Newton")
# create analysis object
op.analysis("Static")
nstep = 100
temp = [0.0]
disp = [0.0]
for i in range(nstep):
if op.analyze(1) < 0:
print("Second analysis failed.")
break
temp.append(op.getLoadFactor(patternTag)*maxtemp)
disp.append(op.nodeDisp(70,1))
plt.plot(temp,disp,'-o')
plt.xlabel('Temperature')
plt.ylabel('Nodal displacement')
plt.grid()
plt.show()
# # perform the analysis
# op.initialize()
# ok = op.analyze(1)
# timer stop: op.stop()
vfo.plot_model(ModelName,show_nodetags="yes",show_eletags="yes")
vfo.plot_deformedshape(ModelName, LoadCaseName, scale=100, tstep=1, overlap='yes')
op.wipe()
HEBsecTag, IPEsecTag, WallsecTag, SWallsecTag, ArcsecTag = GetSections()
totelements, SEleArc, Wallmesh, ARCNodedis, y3, x5 = GetModel(HEBsecTag, IPEsecTag, WallsecTag, SWallsecTag, ArcsecTag, Wallmesh)
GetRecorders()
#RunAnalysisStatic(totelements, SEleArc, Wallmesh, ARCNodedis, y3, x5, G, Gs)
RunThermalFunction(totelements, SEleArc, Wallmesh, ARCNodedis, y3, x5, G, Gs)
SectionForceDeformation::getTemperatureStress(double *dataMixed) - should not be called
Moderators: silvia, selimgunay, Moderators