SpyderのOpenseesを使用してプッシュオーバー分析を行っています。 しかし、直線的なパフォーマンス曲線しか得られません。 問題がどこにあるのかを見つけるのを手伝っていただけますか?
私が試したこと:
print("==========================") from openseespy.opensees import * import os, sys import numpy as np import math as math import matplotlib.pyplot as plt print("Starting RCFrameGravity example") wipe() # Create ModelBuilder (with two-dimensions and 3 DOF/node) model('basic', '-ndm', 2, '-ndf', 3) # Create nodes # ------------ # Set parameters for overall model geometry height = 1300.0 # Create nodes # tag, X, Y node(1, 0.0, 0.0) node(2, 0.0, 0.0) node(3, 0.0, height) # Fix supports at base of columns # tag, DX, DY, RZ fix(1, 1, 1, 1) fix(2, 1, 1, 0) #%% # Define materials uniaxialMaterial('Hysteretic', 1, 0.006, 75514606.25, 0.03, 75514606.25*1.1, -0.006, -75514606.25, -0.03, -75514606.25*1.1, 1.0, 1.0, 0.0, 0.0, 0.7) #%% # Define elements # ---------------------- # Geometry of column elements # tag geomTransf('PDelta', 1) # Number of integration points along length of element np = 10 # Lobatto integratoin beamIntegration('Lobatto', 1, 1, np) # Create the coulumns element('zeroLength', 1, 1, 2, '-mat',1, '-dir',6) bc = 250.0 dc = 250.0 A = bc*dc Iy = bc*dc**3/12 Ec = 26692.696 element('elasticBeamColumn', 2, 2, 3, A, Ec, Iy, 1) # import vfo rendering module import vfo.vfo as vfo # render the model after defining all the nodes and elements # vfo.plot_model(show_eletags="yes") # plot mode shape # vfo.plot_modeshape(modenumber=3) #%% # Define gravity loads # -------------------- # a parameter for the axial load P = 200000.0 # Create a Plain load pattern with a Linear TimeSeries timeSeries('Linear', 1) pattern('Plain', 1, 1) # Create nodal loads at nodes 3 # nd FX, FY, MZ # load(3, -P, 0.0, 0.0) load(3, 0.0, -P, 0.0) # ------------------------------ # End of model generation # ------------------------------ # ------------------------------ # Start of analysis generation # ------------------------------ # Create the system of equation, a sparse solver with partial pivoting system('ProfileSPD') # Create the constraint handler, the transformation method constraints('Transformation') # Create the DOF numberer, the reverse Cuthill-McKee algorithm numberer('RCM') # Create the convergence test, the norm of the residual with a tolerance of # 1e-12 and a max number of iterations of 50 test('NormDispIncr', 1.0e-12, 50, 3) # Create the solution algorithm, a Newton-Raphson algorithm algorithm('Newton') # Create the integration scheme, the LoadControl scheme using steps of 0.1 integrator('LoadControl', 0.1) # Create the analysis object analysis('Static') # ------------------------------ # End of analysis generation # ------------------------------ # ------------------------------ # Finally perform the analysis # ------------------------------ # perform the gravity load analysis, requires 50 steps to reach the load level analyze(50) # Print out the state of nodes 3 # print node 3 # Print out the state of element 1 and 2 # print ele 1 2 u3 = nodeDisp(3, 1) results = open('results.out', 'a+') results.close() print("==========================") #%% # ---------------------------------------------------- # Start of Model Generation & Initial Gravity Analysis # ---------------------------------------------------- # Set the gravity loads to be constant & reset the time in the domain loadConst('-time', 0.0) # ---------------------------------------------------- # End of Model Generation & Initial Gravity Analysis # ---------------------------------------------------- # ---------------------------------------------------- # Start of additional modelling for lateral loads # ---------------------------------------------------- # Define lateral loads # -------------------- # Set some parameters H = 200000.0 # Reference lateral load # Set lateral load pattern with a Linear TimeSeries pattern('Plain', 2, 1) # Create nodal loads at nodes 3 # nd FX FY MZ load(3, H, 0.0, 0.0) # ---------------------------------------------------- # End of additional modelling for lateral loads # ---------------------------------------------------- # ---------------------------------------------------- # Start of modifications to analysis for push over # ---------------------------------------------------- # Set some parameters dU = 1.0 # Displacement increment # Change the integration scheme to be displacement control # node dof init Jd min max integrator('DisplacementControl', 3, 1, dU, 1, dU, dU) # ---------------------------------------------------- # End of modifications to analysis for push over # ---------------------------------------------------- # ------------------------------ # Start of recorder generation # ------------------------------ # Stop the old recorders by destroying them # remove recorders # Create a recorder to monitor nodal displacements # recorder Node -file node3.out -time -node 3 -dof 1 2 3 disp # Create a recorder to monitor element forces in columns # recorder EnvelopeElement -file ele12.out -time -ele 1 2 forces # -------------------------------- # End of recorder generation # --------------------------------- # ------------------------------ # Finally perform the analysis # ------------------------------ # Set some parameters maxU = 250.0 # Max displacement currentDisp = 0.0 ok = 0 test('NormDispIncr', 1.0e-12, 1000) algorithm('ModifiedNewton', '-initial') topdispl = [] #empty list to record top node displacement ndispl = {1: [], 2: [], 3: []} eforce = {(1,1): [], (1,2):[], (1,3):[], (1,4):[], (1,5):[], (1,6):[], (2,1): [], (2,2):[], (2,3):[], (2,4):[], (2,5):[], (2,6):[]} while ok == 0 and currentDisp < maxU: ok = analyze(1) print('currentDisp =',currentDisp) # if the analysis fails try initial tangent iteration if ok != 0: print("modified newton failed") break # print "regular newton failed .. lets try an initail stiffness for this step" # test('NormDispIncr', 1.0e-12, 1000) # # algorithm('ModifiedNewton', '-initial') # ok = analyze(1) # if ok == 0: # print "that worked .. back to regular newton" # test('NormDispIncr', 1.0e-12, 10) # algorithm('Newton') currentDisp = nodeDisp(3, 1) #record each step topdispl.append(nodeDisp(3,1)) for ntag in [1,2,3]: ndispl[ntag].append(nodeDisp(ntag,1)) for etag in [1,2]: eforce[etag,1].append(eleForce(etag,1)) eforce[etag,2].append(eleForce(etag,2)) eforce[etag,3].append(eleForce(etag,3)) eforce[etag,4].append(eleForce(etag,4)) eforce[etag,5].append(eleForce(etag,5)) eforce[etag,6].append(eleForce(etag,6)) xpoints=ndispl[3] ypoints=eforce[2,1] plt.plot(xpoints, ypoints) results = open('results.out', 'a+') results.close() # Print the state at node 3 # print node 3 print("==========================")
コメント