Convergence Issues with SFI_MVLEM Element in OpenSeesPy Model

Forum for asking and answering questions related to use of the OpenSeesPy module

Moderators: silvia, selimgunay, Moderators

Post Reply
DJERRAD
Posts: 2
Joined: Thu Aug 05, 2021 1:11 am

Convergence Issues with SFI_MVLEM Element in OpenSeesPy Model

Post by DJERRAD »

Hello OpenSees users,
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

Thank you in advance for your assistance, and I look forward to your feedback.
Best regards,
Post Reply