Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

null matrix exported for C for modal damping, GimmeMCK #105

Open
jad418 opened this issue Feb 24, 2023 · 1 comment
Open

null matrix exported for C for modal damping, GimmeMCK #105

jad418 opened this issue Feb 24, 2023 · 1 comment

Comments

@jad418
Copy link

jad418 commented Feb 24, 2023

Hello,

I am attempting to do a damping study where I need the damping matrices for different damping models. Currently when I run GimmeMCK for modal damping, I get a null matrix as the result. Per Michael Scott's suggestion, I updated OpenseesPy to the most recent version, but I am still having the same issue with my code.

from openseespy.opensees import *
import numpy as np
import matplotlib.pyplot as plt
import os
import vfo.vfo as vfo

import scipy.linalg as slin

wipe()

ndm = 2;
model('basic', '-ndm', 2, '-ndf', 3);

#establish units
inch = 1;
kip = 1;
sec = 1;

Dependent units

sq_in = inchinch;
ksi = kip/sq_in;
ft = 12
inch;

Constants

g = 386.2inch/(secsec);
#pi = math.acos(-1);

#establish the geometry of the 8-story SMF. Use centerline dimension model initially
NStories = 3.0;
NBays = 1.0;
WidthBay = 20.0ft;#in
HtStoryTyp = 20.0
ft;#in

ColLine1 = 0.0;
ZeroLengthColLine1 = ColLine1+0.0;
ZeroLengthColLine2 = ZeroLengthColLine1 + WidthBay;
ColLine2 = ZeroLengthColLine2+0.0;

Base = 0.0;
Floor1 = HtStoryTyp;
Floor2 = 2HtStoryTyp;
Floor3 = 3
HtStoryTyp;

#Establish the node geometry for the centerline model
node(10, ColLine1,Base);
node(20, ColLine2,Base);
node(11, ColLine1,Floor1);
node(101, ZeroLengthColLine1,Floor1);
node(201, ZeroLengthColLine2,Floor1);
node(21, ColLine2,Floor1);
node(12, ColLine1,Floor2);
node(102, ZeroLengthColLine1,Floor2);
node(202, ZeroLengthColLine2,Floor2);
node(22, ColLine2,Floor2);
node(13, ColLine1,Floor3);
node(103, ZeroLengthColLine1,Floor3);
node(203, ZeroLengthColLine2,Floor3);
node(23, ColLine2,Floor3);

#Define the fixities, setup so that the base of columns have fixed-bases
fix(10,1,1,1);
fix(20,1,1,1);

#apply the masses to the column nodes
FloorColMass1 = 50kip/g;
FloorColMass2 = 50
kip/g;
FloorColMass3 = 50*kip/g;

mass(10,0.0,0.0,0.0);
mass(20,0.0,0.0,0.0);
mass(11,FloorColMass1,0.0,0.0);
mass(21,FloorColMass1,0.0,0.0);
mass(12,FloorColMass2,0.0,0.0);
mass(22,FloorColMass2,0.0,0.0);
mass(13,FloorColMass3,0.0,0.0);
mass(23,FloorColMass3,0.0,0.0);

#Add EQ DOf Contraint
#Floor 1
equalDOF(11, 101, 1);
equalDOF(101, 201, 1);
equalDOF(201, 21, 1);
#Floor 2
equalDOF(12, 102, 1);
equalDOF(102, 202, 1);
equalDOF(202, 22, 1);
#Floor 3
equalDOF(13, 103, 1);
equalDOF(103, 203, 1);
equalDOF(203, 23, 1);

#yielding bi-linear material
matType = 'Steel01'
matTag1 = 1
Fy = 110**9ksi;
E0 = 110**9ksi;
b = 0.00000001;# null to be EPP?
matArgs1 = [Fy, E0, b]
uniaxialMaterial(matType, matTag1, *matArgs1)

#elastic material
matType = 'Steel01'
matTagE1 = 2
Fy = 29000ksi;
E0 = 29000
ksi;
b = 0.02;
matArgs2 = [Fy, E0, b]
uniaxialMaterial(matType, matTagE1, *matArgs2)

#establish the beam and column elements

Assign geometric transformation

ColTransfTag=3
BeamTranfTag=4

geomTransf('PDelta', ColTransfTag)
geomTransf('Linear', BeamTranfTag)

#define iteration limits
maxIter=10;
tol=1e-12;

Area = 100inch;
E_mod = 29000
ksi;
Iz = 4000inchinchinchinch;

#create columns
element('elasticBeamColumn', 31011, 10,11, Area, E_mod, Iz, ColTransfTag)
element('elasticBeamColumn', 32021, 20,21, Area, E_mod, Iz, ColTransfTag)
element('elasticBeamColumn', 31112, 11,12, Area, E_mod, Iz, ColTransfTag)
element('elasticBeamColumn', 32122, 21,22, Area, E_mod, Iz, ColTransfTag)
element('elasticBeamColumn', 31213, 12,13, Area, E_mod, Iz, ColTransfTag)
element('elasticBeamColumn', 32223, 22,23, Area, E_mod, Iz, ColTransfTag)

#create beams
element('elasticBeamColumn', 4101201, 101,201, Area, E_mod, Iz, ColTransfTag)
element('elasticBeamColumn', 4102202, 102,202, Area, E_mod, Iz, ColTransfTag)
element('elasticBeamColumn', 4103203, 103,203, Area, E_mod, Iz, ColTransfTag)

#Create zero-length rotational springs
element('zeroLength', 511101, 11,101, '-mat', matTag1, '-dir', 3)
element('zeroLength', 521201, 21,201, '-mat', matTag1, '-dir', 3)
element('zeroLength', 512102, 12,102, '-mat', matTag1, '-dir', 3)
element('zeroLength', 522202, 22,202, '-mat', matTag1, '-dir', 3)
element('zeroLength', 513103, 13,103, '-mat', matTag1, '-dir', 3)
element('zeroLength', 523203, 23,203, '-mat', matTag1, '-dir', 3)

Nmodes = 3
w2 = eigen('-fullGenLapack',Nmodes)

w = np.sqrt(w2)
print(w)
zeta = 0.03

damping = [zeta]*Nmodes

modalDampingQ(*damping)

alphaM = 2zeta(w[0]*w[2])/(w[0]+w[2]);

betaK = 2*zeta/(w[0]+w[2]);

rayleigh(alphaM, betaK, 0.0, 0.0)

zeta = 0.03
damping = [zeta]*Nmodes
modalDampingQ(*damping)

Tol = 1e-8
maxNumIter = 10
GMdirection = 1
GMfact = 1.0
GMfatt = gGMfact
DtAnalysis = 0.005
sec # time-step Dt for lateral analysis
TmaxAnalysis = 30.0*sec # maximum duration of ground-motion analysis
print("1")

Set some parameters

record = 'A-TMZ000'

import ReadRecord

Permform the conversion from SMD record to OpenSees record

dt, nPts = ReadRecord.ReadRecord(record+'.at2', record+'.dat')
#print(dt, nPts)

Uniform EXCITATION: acceleration input

IDloadTag = 400 # load tag
TStag = 2222222;
timeSeries('Path', TStag, '-dt', dt, '-filePath', record+'.dat', '-factor', GMfatt)
pattern('UniformExcitation', IDloadTag, GMdirection, '-accel', TStag)

vfo.plot_model(model='3storyModalOnlyElastic', show_nodes='yes', show_nodetags='no', show_eletags='no', font_size=10, setview='3D', elementgroups=None, line_width=1, filename=None)

recorder('Node', '-file', "NodeDisp", '-xml', "NodeDisp", '-binary', "NodeDisp", '-tcp', '-precision', nSD=6, '-timeSeries', tsTag, '-time', '-dT', deltaT=0.0, '-closeOnWrite', '-node', *nodeTags=[], '-nodeRange', startNode, endNode, '-region', regionTag, '-dof', *dofs=[], respType)

allNodes = getNodeTags()
print(allNodes)
ColLine1Nodes = [10,11,12,13];
ColLine2Nodes = [20,21,22,23];
recorder('Node','-file','NodalDispColYielding1','-time','-dT',dt,'-node',*ColLine1Nodes,'-dof',1,'disp')
recorder('Node','-file','NodalDispColYielding2','-time','-dT',dt,'-node',*ColLine2Nodes,'-dof',1,'disp')

wipeAnalysis()
constraints('Transformation')
numberer('RCM')
system('FullGeneral')

test('EnergyIncr', Tol, maxNumIter)

test('NormUnbalance', 0.001, 100, 0)
algorithm('KrylovNewton')

algorithm('ModifiedNewton')

NewmarkGamma = 0.5
NewmarkBeta = 0.25
integrator('Newmark', NewmarkGamma, NewmarkBeta)
print("check")
analysis('Transient')
analyze(nPts,dt)
print("2")

Mass

integrator('GimmeMCK',1.0,0.0,0.0)
analyze(1,0.0)

Number of equations in the model

N = systemSize() # Has to be done after analyze

M = printA('-ret') # Or use ops.printA('-file','M.out')
M = np.array(M) # Convert the list to an array
M.shape = (N,N) # Make the array an NxN matrix

Stiffness

integrator('GimmeMCK',0.0,0.0,1.0)
analyze(1,0.0)
K = printA('-ret')
K = np.array(K)
K.shape = (N,N)

Damping

integrator('GimmeMCK',0.0,1.0,0.0)
analyze(1,0.0)
C = printA('-ret')
C = np.array(C)
C.shape = (N,N)

modalDamping = damping ;# Or whatever
Nmodes = len(modalDamping)

print('Analysis M C and K Run')

@mhscott
Copy link
Collaborator

mhscott commented Feb 24, 2023

Please post a minimal working example (MWE)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants