In openseespy, i am trying to get mass, stiffness and damping matrices for mdof lumped model which have isolator in midstory. GimmeMCK does not give me right damping matrix but it gives right mass and stiffness matrices. How can i get right damping matrix?
Openseespy code given below
Code: Select all
# -*- coding: utf-8 -*-
"""
Created on Mon Jan 30 14:49:55 2023
@author: Berkan
"""
import openseespy.opensees as ops
import numpy as np
ops.wipe()
x1 = 0
x2 = 1
x3 = 2
x4 = 3
x5 = 4
x6 = 5
x7 = 6
x8 = 7
x9 = 8
x10 = 9
x11 = 10
x12 = 11
x13 = 12
x14 = 13
x15 = 14
x16 = 15
x17 = 16
x18 = 17
x19 = 18
x20 = 19
x21 = 20
k = 100000 #stiffness
m = 100000/980.665 #mass
ro = 0.05#damping ratio
A = 1 #Area of section
# =============================================================================
# OpenSees Analysis
# =============================================================================
#create Model
ops.model('basic','-ndm',1,'-ndf',1)
# create nodes
ops.node(1,x1)
ops.node(2,x2,'-mass',m)
ops.node(3,x3,'-mass',m)
ops.node(4,x4,'-mass',m)
# ops.node(5,x5,'-mass',m)
# ops.node(6,x6,'-mass',m)
# ops.node(7,x7,'-mass',m)
# ops.node(8,x8,'-mass',m)
# ops.node(9,x9,'-mass',m)
# ops.node(10,x10,'-mass',m)
# ops.node(11,x11,'-mass',m)
# ops.node(12,x12,'-mass',m)
# ops.node(13,x13,'-mass',m)
# ops.node(14,x14,'-mass',m)
# ops.node(15,x15,'-mass',m)
# ops.node(16,x16,'-mass',m)
# ops.node(17,x17,'-mass',m)
# ops.node(18,x18,'-mass',m)
# ops.node(19,x19,'-mass',m)
# ops.node(20,x20,'-mass',m)
# ops.node(21,x21,'-mass',m)
# set boundary condition
ops.fix(1,1)
# =============================================================================
# SubStructure
# =============================================================================
# define materials
ops.uniaxialMaterial('Elastic',1,k,0)
# define elements
ops.element("Truss",1,1,2,A,1,'-doRayleigh',1)
ops.element("Truss",2,2,3,A,1,'-doRayleigh',1)
ops.element("Truss",3,3,4,A,1,'-doRayleigh',1)
# =============================================================================
# FINDING DAMPING MATRIX
# =============================================================================
# eigenvalue = ops. eigen('-fullGenLapack', 2)
eigenvalue = ops. eigen('-genBandArpack', 2)
naturalfrequency = np.sqrt(eigenvalue)
a0 = ro * (2*naturalfrequency[0]*naturalfrequency[1])/(naturalfrequency[0]+naturalfrequency[1])
a1 = ro * (2)/(naturalfrequency[0]+naturalfrequency[1])
ops.rayleigh(a0, a1, 0.0, 0.0)
ops.wipeAnalysis()
ops.system('FullGeneral')
ops.analysis('Transient')
# Mass
ops.integrator('GimmeMCK',1.0,0.0,0.0)
ops.analyze(1,0.0)
# Number of equations in the model
N = ops.systemSize() # Has to be done after analyze
M = ops.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
ops.integrator('GimmeMCK',0.0,0.0,1.0)
ops.analyze(1,0.0)
K = ops.printA('-ret')
K = np.array(K)
K.shape = (N,N)
#damping
ops.integrator('GimmeMCK', 0.0, 1.0, 0.0)
ops.analyze(1, 0.0)
C = ops.printA('-ret')
C = np.array(C)
C.shape = (N,N)
print('\n', 'M=', '\n', M)
print('\n', 'K=', '\n', K)
print('\n', 'C=', '\n', C)
print('\n', 'N=', '\n', N)
C1 = a0*M + a1*K
print('\n', 'C1=', '\n', C1)
# =============================================================================
# Isolation layer
# =============================================================================
# ops.node(4,x4,'-mass',m)
ops.node(5,x5,'-mass',m)
kiso = 56.0564
ops.uniaxialMaterial('Elastic',3,kiso,0)
ro_iso = 0.5
ops.element("Truss",4,4,5,A,3,'-doRayleigh',1)
eigenvalue = ops. eigen('-fullGenLapack', 2)
naturalfrequency = np.sqrt(eigenvalue)
a0 = ro_iso * (2*naturalfrequency[0]*naturalfrequency[1])/(naturalfrequency[0]+naturalfrequency[1])
a1 = ro_iso * (2)/(naturalfrequency[0]+naturalfrequency[1])
ops.rayleigh(a0, a1, 0.0, 0.0)
# ops.rayleigh(0, 0, 0.0, 0.0)
ops.wipeAnalysis()
ops.system('FullGeneral')
ops.analysis('Transient')
# Mass
ops.integrator('GimmeMCK',1.0,0.0,0.0)
ops.analyze(1,0.0)
# Number of equations in the model
N = ops.systemSize() # Has to be done after analyze
M = ops.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
ops.integrator('GimmeMCK',0.0,0.0,1.0)
ops.analyze(1,0.0)
K = ops.printA('-ret')
K = np.array(K)
K.shape = (N,N)
#damping
ops.integrator('GimmeMCK', 0.0, 1.0, 0.0)
ops.analyze(1, 0.0)
C = ops.printA('-ret')
C = np.array(C)
C.shape = (N,N)
print('\n', 'M=', '\n', M)
print('\n', 'K=', '\n', K)
print('\n', 'C=', '\n', C)
print('\n', 'N=', '\n', N)
C1 = a0*M + a1*K
print('\n', 'C1=', '\n', C1)
# =============================================================================
# SupStructure
# =============================================================================
# ops.node(5,x5,'-mass',m)
ops.node(6,x6,'-mass',m)
ops.node(7,x7,'-mass',m)
ops.node(8,x8,'-mass',m)
ops.node(9,x9,'-mass',m)
ops.node(10,x10,'-mass',m)
ops.node(11,x11,'-mass',m)
ops.node(12,x12,'-mass',m)
ops.node(13,x13,'-mass',m)
ops.node(14,x14,'-mass',m)
ops.node(15,x15,'-mass',m)
ops.node(16,x16,'-mass',m)
ops.node(17,x17,'-mass',m)
ops.node(18,x18,'-mass',m)
ops.node(19,x19,'-mass',m)
ops.node(20,x20,'-mass',m)
ops.node(21,x21,'-mass',m)
ops.uniaxialMaterial('Elastic',2,1*k,0)
ops.element("Truss",5,5,6,A,2,'-doRayleigh',1)
ops.element("Truss",6,6,7,A,2,'-doRayleigh',1)
ops.element("Truss",7,7,8,A,2,'-doRayleigh',1)
ops.element("Truss",8,8,9,A,2,'-doRayleigh',1)
ops.element("Truss",9,9,10,A,2,'-doRayleigh',1)
ops.element("Truss",10,10,11,A,2,'-doRayleigh',1)
ops.element("Truss",11,11,12,A,2,'-doRayleigh',1)
ops.element("Truss",12,12,13,A,2,'-doRayleigh',1)
ops.element("Truss",13,13,14,A,2,'-doRayleigh',1)
ops.element("Truss",14,14,15,A,2,'-doRayleigh',1)
ops.element("Truss",15,15,16,A,2,'-doRayleigh',1)
ops.element("Truss",16,16,17,A,2,'-doRayleigh',1)
ops.element("Truss",17,17,18,A,2,'-doRayleigh',1)
ops.element("Truss",18,18,19,A,2,'-doRayleigh',1)
ops.element("Truss",19,19,20,A,2,'-doRayleigh',1)
ops.element("Truss",20,20,21,A,2,'-doRayleigh',1)
eigenvalue = ops. eigen('-genBandArpack', 2)
naturalfrequency = np.sqrt(eigenvalue)
a0 = ro * (2*naturalfrequency[0]*naturalfrequency[1])/(naturalfrequency[0]+naturalfrequency[1])
a1 = ro * (2)/(naturalfrequency[0]+naturalfrequency[1])
ops.rayleigh(a0, a1, 0.0, 0.0)
ops.wipeAnalysis()
ops.system('FullGeneral')
ops.analysis('Transient')
# Mass
ops.integrator('GimmeMCK',1.0,0.0,0.0)
ops.analyze(1,0.0)
# Number of equations in the model
N = ops.systemSize() # Has to be done after analyze
M = ops.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
ops.integrator('GimmeMCK',0.0,0.0,1.0)
ops.analyze(1,0.0)
K = ops.printA('-ret')
K = np.array(K)
K.shape = (N,N)
#damping
ops.integrator('GimmeMCK', 0.0, 1.0, 0.0)
ops.analyze(1, 0.0)
C = ops.printA('-ret')
C = np.array(C)
C.shape = (N,N)
print('\n', 'M=', '\n', M)
print('\n', 'K=', '\n', K)
print('\n', 'C=', '\n', C)
print('\n', 'N=', '\n', N)
C1 = a0*M + a1*K
print('\n', 'C1=', '\n', C1)