I'm currently working on an OpenSees model in Python using the SFI_MVLEM element. However, I'm facing an issue where my model won't converge during the analysis, and I'm seeking some guidance to resolve this problem.
I'd like to add that I have successfully used the MVLEM element with similar settings, and it converged without any issues for many model.
I tried in :
Python 3.8, 3.10 and 3.11
OpenSeesPy 3.5.18 and 3.3.0.1.1
Code: Select all
import openseespy.opensees as ops
import numpy as np
import time
import math
# ------------------------------
# Start of model generation
# -----------------------------
mm = 1.0 # 1 millimeter
N = 1.0 # 1 Newton
sec = 1.0 # 1 second
m = 1000.0 * mm # 1 meter is 1000 millimeters
cm = 10.0 * mm # 1 centimeter is 10 millimeters
kN = 1000.0 * N # 1 kilo-Newton is 1000 Newtons
m2 = m * m # Square meter
cm2 = cm * cm # Square centimeter
mm2 = mm * mm # Square millimeter
MPa = N / mm2 # MegaPascal (Pressure)
kPa = 0.001 * MPa # KiloPascal (Pressure)
GPa = 1000 * MPa # GigaPascal (Pressure)
def generate_cyclic_load(duration=6.0, sampling_rate=50, max_displacement=75):
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
displacement_slope = (max_displacement / 2) / (duration / 2)
cyclic_load = (displacement_slope * t) * np.sin(2 * np.pi * t)
return cyclic_load
def Thomsen_and_Wallace_RW2():
# https://sci-hub.st/10.1061/(asce)0733-9445(2004)130:4(618)
global name
name = 'Thomsen_and_Wallace_RW2'
# Wall Geometry ------------------------------------------------------------------
tw = 102 * mm # Wall thickness
hw = 3.66 * m # Wall height
lw = 1.22 * m # Wall length
lbe = 190 * mm # Boundary element length
lweb = lw - (2 * lbe)
# Material proprieties -----------------------------------------------------------
fc = 41.7 * MPa # Concrete peak compressive stress (+Tension, -Compression)
fy = 414 * MPa # Steel tension yield strength (+Tension, -Compression)
loadcoef = 0.1
DisplacementStep = generate_cyclic_load(duration=10, sampling_rate=50, max_displacement=75)
return tw, hw, lw, lbe, fc, fy, loadcoef, DisplacementStep
ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3)
# ---------------- load geometry and material ------------
validation_model = Thomsen_and_Wallace_RW2()
tw, hw, lw, lbe, fc, fy, loadcoef, DisplacementStep = validation_model
# --------------- Discretization of the model ------------
eleH = 16
eleL = 8
wall_thickness = tw # Wall thickness
wall_height = hw # Wall height
wall_length = lw # Wall width
length_be = lbe # Length of the Boundary Element
length_web = lweb = lw - (2 * lbe) # Length of the Web
m = eleH
n = eleL
eleBE = 2
eleWeb = eleL - eleBE
elelweb = lweb / eleWeb
# --------------- create nodes ---------------
# Loop through the list of node values
for i in range(1, eleH + 2):
ops.node(i, 0, (i - 1) * (hw / eleH))
ops.fix(1, 1, 1, 1) # Fixed condition at node 1
IDctrlNode = eleH + 1 # Control Node (TopNode)
IDctrlDOF = 1 # Control DOF 1 = X-direction
# ---------------------------------------------------------------------------------------
# Define uni-axial materials
# ---------------------------------------------------------------------------------------
sYb = 1 # STEEL Y (boundary element)
sYw = 2 # STEEL Y (web)
sX = 3 # STEEL X (boundary element + web)
# STEEL Y boundary element
fyYbp = fy # fy - tension
fyYbn = fy # fy - compression
bybt = 0.0185 # strain hardening - tension
bybc = 0.02 # strain hardening - compression
# STEEL Y web
fyYwp = fy # fy - tension
fyYwn = fy # fy - compression
bywt = 0.035 # strain hardening - tension
bywc = 0.02 # strain hardening - compression
# STEEL X
fyXt = fy # fy - tension
fyXc = fy # fy - compression
bXt = 0.0185 # strain hardening - tension
bXc = 0.02 # strain hardening - compression
# STEEL misc
Es = 200 * GPa # Young's modulus
R0 = 15.0 # initial value of curvature parameter
a1 = 0.925 # curvature degradation parameter
a2 = 0.15 # curvature degradation parameter
Bs = 0.01 # strain-hardening ratio
cR1 = 0.925 # control the transition from elastic to plastic branches
cR2 = 0.015 # control the transition from elastic to plastic branches
# STEEL Reinforcing steel
ops.uniaxialMaterial('SteelMPF', sYb, fyYbp, fyYbn, Es, bybt, bybc, R0, cR1, cR2) # steel Y boundary
ops.uniaxialMaterial('SteelMPF', sYw, fyYwp, fyYwn, Es, bywt, bywc, R0, cR1, cR2) # steel Y web
ops.uniaxialMaterial('SteelMPF', sX, fyXt, fyXc, Es, bXt, bXc, R0, cR1, cR2) # steel X
# ------------------------------CONCRETE ----------------
concBE = 4
concWeb = 5
# unconfined
fpc = -fc
Ec = 35 * GPa
ec0 = 2 * (fpc / Ec)
ru = 7
xcrnu = 1.039
# confined
Ecc = Ec # Young's modulus
fpcc = fpc * 1.2 # Concrete Compressive Strength, MPa (+Tension, -Compression)
ec0c = 2 * (fpcc / Ec) # strain at maximum compressive stress (-0.0033)
rc = 7.3049 # shape parameter - compression
xcrnc = 1.0125 # cracking strain - compression
ft = 2.05 * MPa # peak tensile stress
et = 0.00008 # strain at peak tensile stress
rt = 1.2 # shape parameter - tension
xcrp = 10000 # cracking strain - tension
ops.uniaxialMaterial('ConcreteCM', concBE, fpcc, ec0c, Ecc, rc, xcrnc, ft, et, rt, xcrp, '-GapClose', 0) # confined concrete -- ec0c
ops.uniaxialMaterial('ConcreteCM', concWeb, fpc, ec0, Ec, ru, xcrnu, ft, et, rt, xcrp, '-GapClose', 0) # unconfined concrete -- ec0
# ----------------------------Shear spring for MVLEM-------------------------------
Ac = lw * tw # Concrete Wall Area
Gc = Ec / (2 * (1 + 0.2)) # Shear Modulus G = E / 2 * (1 + v)
Kshear = Ac * Gc * (5 / 6) # Shear stiffness k * A * G ---> k=5/6
ops.uniaxialMaterial('Elastic', 6, Kshear) # Shear Model for Section Aggregator
Aload = 0.85 * abs(fc) * tw * lw * loadcoef # Aload = 0.07 * abs(fc) * tw * lw * loadcoef
# ---- Steel in Y direction (BE + Web) -------------------------------------------
rouYb = 0.029488496
rouYw = 0.00376992
# ---- Steel in X direction (BE + Web) -------------------------------------------
rouXb = 0.0028274496
rouXw = 0.0028274496
# ----------------------------FSAM Material -------------------------------
matBE = 7
matWeb = 8
ops.nDMaterial('FSAM', matBE, 0.0, sX, sYb, concBE, rouXb, rouYb, 0.2, 0.012) # Boundary (confined concrete only)
ops.nDMaterial('FSAM', matWeb, 0.0, sX, sYw, concWeb, rouXw, rouYw, 0.2, 0.012) # Web (unconfined concrete)
# ------------------------------------------------------------------------
# Define 'MVLEM' or 'SFI-MVLEM' elements
# ------------------------------------------------------------------------
# Set 'MVLEM' parameters thick, width, rho, matConcrete, matSteel
MVLEM_thick = [tw] * n
MVLEM_width = [lbe if i in (0, n - 1) else elelweb for i in range(n)]
MVLEM_rho = [rouYb if i in (0, n - 1) else rouYw for i in range(n)]
MVLEM_matConcrete = [concBE if i in (0, n - 1) else concWeb for i in range(n)]
MVLEM_matSteel = [sYb if i in (0, n - 1) else sYw for i in range(n)]
MVLEM_mat = [matBE if i in (0, n - 1) else matWeb for i in range(n)]
for i in range(eleH):
# ---------------- MVLEM--------------------------------
# ops.element('MVLEM', i + 1, 0.0, *[i + 1, i + 2], eleL, 0.4, '-thick', *MVLEM_thick, '-width', *MVLEM_width, '-rho', *MVLEM_rho, '-matConcrete', *MVLEM_matConcrete, '-matSteel', *MVLEM_matSteel, '-matShear', 6)
# ---------------- SFI_MVLEM--------------------------------
ops.element('SFI_MVLEM', i + 1, *[i + 1, i + 2], eleL, 0.4, '-thick', *MVLEM_thick, '-width', *MVLEM_width, '-mat', *MVLEM_mat)
print('SFI_MVLEM', i + 1, *[i + 1, i + 2], eleL, 0.4, '-thick', *MVLEM_thick, '-width', *MVLEM_width, '-mat', *MVLEM_mat)
# ------------------------------
# Start of analysis generation
# ------------------------------
# -------------GRAVITY-----------------
steps = 10
ops.timeSeries('Linear', 1, '-factor', 1.0) # create TimeSeries for gravity analysis
ops.pattern('Plain', 1, 1)
ops.load(IDctrlNode, *[0.0, -Aload, 0.0]) # apply vertical load
ops.constraints('Transformation')
ops.numberer('RCM')
ops.system('BandGeneral')
ops.test('NormDispIncr', 1.0e-6, 100, 0)
ops.algorithm('Newton')
ops.integrator('LoadControl', 1 / steps)
ops.analysis('Static')
ops.analyze(steps)
ops.loadConst('-time', 0.0)
# -------------CYCLIC-----------------
ops.timeSeries('Linear', 2, '-factor', 1.0)
ops.pattern('Plain', 2, 2)
ops.load(IDctrlNode, *[1.0, 0.0, 0.0]) # apply vertical load
ops.constraints('Transformation') # Transformation 'Penalty', 1e20, 1e20
ops.numberer('RCM')
ops.system("BandGen")
ops.test('NormDispIncr', 1e-8, 100, 0)
ops.algorithm('KrylovNewton')
maxUnconvergedSteps = 20
unconvergeSteps = 0
Nsteps = len(DisplacementStep)
finishedSteps = 0
dispData = np.zeros(Nsteps + 1)
baseShearData = np.zeros(Nsteps + 1)
# Perform cyclic analysis
D0 = 0.0
for j in range(Nsteps):
D1 = DisplacementStep[j]
Dincr = D1 - D0
print(f'Step {j} -------->', f'Dincr = ', Dincr)
if unconvergeSteps > maxUnconvergedSteps:
break
ops.integrator("DisplacementControl", IDctrlNode, 1, Dincr)
ok = ops.analyze(1)
D0 = D1 # move to next step
if ok < 0:
unconvergeSteps = unconvergeSteps + 1
finishedSteps = j + 1
disp = ops.nodeDisp(IDctrlNode, 1)
baseShear = -ops.getLoadFactor(2) / 1000 # Convert to from N to kN
# ------------ Displacement VS Shear ---------------
dispData[j + 1] = disp
baseShearData[j + 1] = baseShear
Best regards,