Fixing boundary condition for zero length section

[Prafullamalla]

I am confused to fix boundary for zerolengthsection. Should I fix x and y direction and release the rotation?
ops.node(3, 0.0,0.0)
ops.node(1, 0.0,0.0);
ops.node(2, 0.0, H)
ops.fix (3, 1, 1, 1);
# Fix for zerolengthsection
ops.fix (1, 1, 1, 0);
ops.element ('zeroLengthSection' ,2, 3, 1,.....

[mhscott]

You should make the local x-axis of the zeroLengthSection [0.1.0], then you only have to fix the global X DOF of node 1.

I also don't recommend using zeroLengthSection in this way, but that's another story.
https://portwooddigital.com/2022/04/24/ ... gth-logic/

[Prafullamalla]

Thankyou. I understood the local axis x will orient along the global Y-axis with the zeroLengthSection [0,1,0]. But what is the reason, why we are restraining in X direction?

[Prafullamalla]

There is problem if I dont restraing Node 3 in global Y direction

Code: Select all

# -*- coding: utf-8 -*-
"""
Created on Fri Jan 20 09:46:58 2023

@author: praf_
"""
import openseespy.opensees as ops
import opsvis as opsv
import numpy as np
import matplotlib.pyplot as plt
from math import asin, sqrt
import os

# Units in kip/in
# column.tcl
# This model is a 2D model

# SET UP ----------------------------------------------------------------------------
ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3); 

# -----------------------------------------------------
# define section geometry
Dcol =24.0; # Column Diameter (in)
LCol =96.0; # column length (in) 
P =147.0 ; # weight (kips)
Weight= P
g =386.22; #[inch/s2]
# calculated parameters
Mass = P/g; # nodal mass

# nodal coordinates:
ops.node(3, 0.0,0.0)
ops.node(1, 0.0,0.0); # node#, X, Y
ops.node(2, 0.0, LCol)


# Single point constraints -- Boundary Conditions
ops.fix (3, 1, 1, 1);
# Fix all degrees of freedom except axial and bending

ops.fix (1,1,0, 0);



# we need to set up parameters that are particular to the model.
IDctrlNode =2; # node where displacement is read for displacement control
IDctrlDOF =1; # degree of freedom of displacement read for displacement control
iSupportNode ="3"; # define support node, if needed.

# Define Column geometry
#-------------------------------------------------------------------
coverSec =0.874; # Column cover to reinforcing steel (in)
numBarsSec =22; # number of longitudinal-reinforcement bars in the column section
dlSec =0.626; # [in] Steel diameter T15.9
PI =3.14;
barAreaSec = PI*dlSec**2/4; #area of longitudinal-reinforcement bars # (in^2)


# MATERIAL parameters
#-------------------------------------------------------------------
IDconcC =1; # material ID tag -- confined core concrete
IDconcU =2; # material ID tag -- unconfined cover concrete
IDreinf =3; # material ID tag -- reinforcement

# Define uniaxial materials
#-------------------------------------------------------------------
# Mander concrete model for confined and unconfined concrete
# characteristics of materials test of Lehman et al. (PEER 98-01)

# confined concrete (core)
#-------------------------
fc =4.496; # [Ksi] CONCRETE Compressive Strength
ec =0.002294; # concrete strain at maximum strength ??????
ecu= 0.008101; # [Ksi] concrete strain at crushing strength ????
Ec =3822.0; # [Ksi]Concrete Elastic Modulus
fct =0.5046; #[Ksi] maximum tensile strength of concrete
et =0.000132; # ultimate tensile strain of concrete
#set beta 0.1; #value defining the exponential curve parameter to define the residual stress

#uniaxialMaterial Concrete04 $matTag $fc $ec $ecu $Ec <$fct $et> <$beta>;
#uniaxialMaterial Concrete04 $IDconcC $fc $ec $ecu $Ec $fct $et;
#uniaxialMaterial Concrete04 $IDconcC [expr -$fc] [expr -$ec] [expr -$ecu] $Ec $fct $et;
ops.uniaxialMaterial ('Concrete04', IDconcC, -fc ,-ec ,-ecu ,Ec ,fct, et);



# unconfined concrete (cover)
#----------------------------
fcU =4.496; # [Ksi] CONCRETE Compressive Strength
ecU =0.002; # concrete strain at maximum strength (#I did not find the experimentation values)
ecuU= 0.0035; # [Ksi] concrete strain at crushing strength (#I did not find the experimentation values)
EcU =3822.0; # [Ksi]Concrete Elastic Modulus
fctU= 0.5046; #[Ksi] maximum tensile strength of concrete
etU =0.000132; # ultimate tensile strain of concrete
#set betaU 0.1; #value defining the exponential curve parameter to define the residual stress

# uniaxialMaterial Concrete04 $matTag $fc $ec $ecu $Ec <$fct $et> <$beta>;
# niaxialMaterial Concrete04 $IDconcU [expr -$fcU] [expr -$ecU] [expr -$ecuU] $EcU $fctU $etU;
ops.uniaxialMaterial ('Concrete04', IDconcU, -fcU, -ecU, -ecuU, EcU, fctU, etU);

IDconcU2 =11
ops.uniaxialMaterial ('Concrete01', IDconcU2, -fcU, -ecU, -fcU*0.99, -ecuU)

# ----------- Reinforcing Steel --------
Es =29000.0 ; # modulus of steel (Ksi)??????
Fy =67.0; # [Ksi] yield stress Lehman et al.2000 (specimen 415.)
Fu= 91.37; #[Ksi] ultimate stress Lehman et al.2000
B =0.01; # strain-hardening ratio (Lehman et al. 2000)
R0 =18.5; # control the transition from elastic to plastic branches
cR1 =0.925; # control the transition from elastic to plastic branches
cR2 =0.15; # control the transition from elastic to plastic branches
#set a1 0.04;
#set a2 1.0;
#set a3 0.04;
#set a4 1.0;

#uniaxialMaterial Steel02 $IDreinf $Fy $Es $B $R0 $cR1 $cR2 $a1 $a2 $a3 $a4;
ops.uniaxialMaterial ('Steel02', IDreinf ,Fy, Es, B, R0, cR1, cR2);

print("materials defined")

#===============================================================Bond-Slip=================================================================
alphaBS= 0.4 ; #alphaBS is a parameter used in the local bond-slip relation and can be taken as 0.4
Pow = (1.0/alphaBS)

#Eq for a model with units (ksi, kips, inches)
#set Sy [expr (((((($dlSec*$Fy*1000.0)/(4000.0*(sqrt($fc*-1000.0))))*((2.0*$alphaBS)+1.0))*pow($Pow,1))*0.1)+0.013)]
Sy =0.08337;
Su = Sy*35.0
b =0.5
R =0.7
bondslipMat= 5

#uniaxialMaterial Bond_SP01 $matTag $Fy $Sy $Fu $Su $b $R
ops.uniaxialMaterial ('Bond_SP01' ,bondslipMat, Fy, Sy, Fu ,Su, b, R)

print ("materials defined ")

#----------------------------------------------------------------
# Generate a circular reinforced concrete section
# with one layer of steel evenly distributed around the perimeter and a confined core.
# confined core.
# by: Michael H. Scott, 2003
# Notes
# The center of the reinforcing bars are placed at the inner radius
# The core concrete ends at the inner radius (same as reinforcing bars)
# The reinforcing bars are all the same size
# The center of the section is at (0,0) in the local axis system
# Zero degrees is along section y-axis
#Following recommendations using uniformly distibutions ratdial disretization shema

ri= 0.0; # inner radius of the section, only for hollow sections
ro = Dcol/2; # overall (outer) radius of the section
nfCoreR =10; # number of radial divisions in the core (number of "rings")
nfCoreT= 20; # number of theta divisions in the core (number of "wedges")
nfCoverR= 1; # number of radial divisions in the cover
nfCoverT= 20; # number of theta divisions in the cover

# Define ELEMENTS & SECTIONS
#----------------------------
ColSecTag =1; # assign a tag number to the column section

# Define the fiber section
#-------------------------
rc = ro-coverSec; # Core radius
bcent= 6.39 ; # expr[$coversec+$Dsec/2]
rb = ro-bcent;

fib_sec_3 = [['section', 'Fiber', 1, '-GJ', 1.0e6],
             ['patch', 'circ', IDconcC, nfCoreT, nfCoreR, 0., 0., ri, rc, 0, 360],
             ['patch', 'circ', IDconcC,  nfCoreT, nfCoverR, 0., 0., rc, ro, 0, 360],
             ['layer', 'circ', IDreinf, numBarsSec, barAreaSec, 0, 0, rc],
             ]


matcolor = ['r', 'lightgrey', 'gold', 'w', 'w', 'w']
opsv.plot_fiber_section(fib_sec_3, matcolor=matcolor)
plt.axis('equal')

ops.section('Fiber', ColSecTag)

# Create the concrete core fibers
ops.patch('circ', IDconcC, nfCoreT, nfCoreR, 0., 0., ri, rc, 0, 360)
ops.patch('circ', IDconcU2,  nfCoreT, nfCoverR, 0., 0., rc, ro, 0, 360)
ops.layer('circ', IDreinf, numBarsSec, barAreaSec, 0, 0, rc)

print ("section geometry defined of Column")
Area=PI/4*Dcol**2
Iy=Iz=PI/64*Dcol**4
Jxx=PI/32*Dcol**4


# Create the concrete core fibers
SecTagBS= 3
ops.section('Fiber', SecTagBS)
ops.patch('circ', IDconcC, nfCoreT, nfCoreR, 0., 0., ri, rc, 0, 360)
ops.patch('circ', IDconcC,  nfCoreT, nfCoverR, 0., 0., rc, ro, 0, 360)
ops.layer('circ', bondslipMat, numBarsSec, barAreaSec, 0, 0, rc)
print ( "section geometry defined of Bond Slip")

#Aggregator
#------------------------
shearsec =4;
TagsecAgg =6;
AreaSec = PI*Dcol**2/4; #gross cross-sectional area (in^2)
G = Ec/2.4;
GA = G*AreaSec;
E_mod=Ec
G_mod=G
#Assign shear to the model of flexure
#------------------------------------

ops.uniaxialMaterial ('Elastic', shearsec, 1e8)
print ("shear defined")

#section Aggregator $TagsecAgg $shearsec Vy -section $ColSecTag
ops.section ('Aggregator', TagsecAgg, shearsec, 'Vy', '-section' ,SecTagBS)

print ("section aggregator defined")


# define geometric transformation
#--------------------------------
#performs a linear geometric transformation of beam stiffness and resisting force
#from the basic system to the global-coordinate system

# Define column elements
# ----------------------
geomTag    = 1
colIntTag  = 1
# Geometry of column elements
#                tag Corotational
ops.geomTransf('Linear', geomTag)

nIP = 3                             # Number of integration points along length of element
# Lobatto integration
ops.beamIntegration('Lobatto', colIntTag, ColSecTag, nIP)
# element connectivity:
eleType = 'forceBeamColumn'                 # accuracy improved by increasing ip or elements
#eleType = 'elasticBeamColumn' 
#eleType = 'nonlinearBeamColumn' 

ops.element(eleType, 1, 1, 2, colIntTag, geomTag)




print ("element nonlinear defined")

######ZeroElement of Bond-Slip
#element zeroLengthSection $EleTag $ibNode $jbNode $SecTagBS
ops.element ('zeroLengthSection' ,2, 3, 1, TagsecAgg ,'-orient' ,0, 1, 0)

# Define gravity loads
#  a parameter for the axial load

# Create a Plain load pattern with a Linear TimeSeries
ops.timeSeries ('Linear', 1)
ops.pattern    ('Plain', 1, 1)

# Create nodal loads at nodes 3 & 4
#    nd  FX,  FY, MZ
ops.load(2, 0.0, -P , 0.0)
# ------------------------------
# End of model generation
# ------------------------------

# ------------------------------
# Start of analysis generation
# ------------------------------
# Gravity Analysis
# Create the system of equation, a sparse solver with partial pivoting
Tol = 1.0e-6
ops.system      ('BandGeneral')
ops.constraints ('Plain')
ops.numberer    ('RCM')
ops.test        ('NormDispIncr', Tol, 10, 3)
ops.algorithm   ('Newton')
ops.integrator  ('LoadControl', 0.1)
ops.analysis    ('Static')
ops.analyze     (10)

#-----maintain constant gravity loads and reset time to zero------
# Set the gravity loads to be constant & reset the time in the domain
ops.loadConst  ('-time', 0.0)
ops.printModel('-node',2)

print ("Model Built")


# Display Model
opsv.plot_model()


print("Pushover Analysis Started!")

# Set some parameters
LateralLoad = 1.0  # Reference lateral load of 1kN
# Set lateral load pattern with a Linear TimeSeries
ops.pattern('Plain', 2, 1)

# Create nodal loads at nodes 3 & 4
#    nd    FX  FY  MZ
ops.load(IDctrlNode, LateralLoad, 0.0, 0.0)

# ----------------------------------------------------
# Start of push over analysis
# ----------------------------------------------------
print("Starting RC column analysis")
maxU = 0.3  # Max displacement
print("Maximum top node displacement", maxU)
# Set some parameters
dU    = 0.002       # Displacement increment
dUMin = 0.002
dUMax = 0.0020

(Tol, MaxIter) = (1.0e-6, 6)      # Tolerance and Maximum iteration (1.0e-12, 10) 
(controlNode, controlDOF) =(2, 1)   # Control node and control DOF
ops.system       ('BandGeneral')
ops.constraints  ('Plain')
ops.numberer     ('RCM')
ops.test         ('NormDispIncr', Tol, MaxIter)
ops.analysis     ('Static');

# Change the integration scheme to be displacement control
#                             node dof init Jd min max
ops.integrator   ('DisplacementControl', controlNode, controlDOF, dU, 1, dUMin,dUMax)
ops.algorithm    ('Newton')

# Set some parameters
currentDisp = ops.nodeDisp(controlNode, controlDOF)
ok          = 0
# Arrays for plotting
Uplot = [0]
LpPlot = [0]
# perform the analysis
Nsteps=int(maxU/dU)

                         # Counter to store data
while (ok == 0 and currentDisp < maxU):
    ops.analyze(1)
    # if the analysis fails try initial tangent iteration
    if ok != 0:
        print("regular newton failed .. lets try an initail stiffness for this step")
        ops.test('NormDispIncr', 1.0e-6, 10,0)
        ops.algorithm('ModifiedNewton', '-initial')
        ok = ops.analyze(1)
        ops.algorithm('Newton')
    currentDisp = ops.nodeDisp(controlNode, controlDOF)
    Uplot.append(currentDisp)
    factorLoad  =ops.getTime()
    LpPlot.append(factorLoad)
    print( "Horizontal load", factorLoad, 'kips,','displacement at control node',currentDisp,'inch')

#P PLOTTING CURVE    
# x axis values
x = Uplot
# corresponding y axis values
y = LpPlot

import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots(figsize=(4, 3))
#plt.ylim(-25,25)
#plt.xlim(-45,45)
ax.plot(x,y)
ax.plot(x,y, color='r', linewidth = 1.0)

ax.set_title('Load Vs Displacement curve', fontsize=10)
ax.set_xlabel('Disp [inch]', fontsize=10)
ax.set_ylabel('Load [kips]', fontsize=10)
#leg = ax.legend()
plt.tick_params(labelsize=6)
plt.tight_layout()
plt.show()
ops.wipe