Reputation: 1
i am testing a simple model using OpenSeespy to check whether i am moving in correct direction but when i make any progress, my kernel keeps dying. i am providing the code here Please, if you all can help, it will be great!
#model setup
import openseespy.opensees as ops
import matplotlib.pyplot as plt
# Wipe any existing model data
ops.wipe()
# Define a 2D model with 3 degrees of freedom per node
ops.model('basic', '-ndm', 2, '-ndf', 3)
# Define Nodes
# Soil nodes
# Base soil nodes (fully fixed)
ops.node(1, 0.0, 0.0) # Node 1
ops.node(2, 5.0, 0.0) # Node 2
ops.node(3, 5.0, 5.0) # Node 3
ops.node(4, 0.0, 5.0) # Node 4
# Intermediate soil nodes (free)
ops.node(5, 2.5, 0.0) # Node 5
ops.node(6, 2.5, 5.0) # Node 6 (Contact Point and Beam Base)
# Additional soil nodes for refined soil model
ops.node(7, 2.5, 10.0) # Node 7 (Beam Top)
ops.node(8, 2.5, 15.0) # Node 8
# Define Materials
# Soil material (Elastic Isotropic)
E_soil = 50e6 # Young's Modulus in Pascals (50 MPa) - Adjusted for softer soil
nu_soil = 0.3 # Poisson's Ratio
ops.nDMaterial('ElasticIsotropic', 1, E_soil, nu_soil)
# Beam section (Elastic)
E_beam = 210e9 # Young's Modulus in Pascals (210 GPa)
A_beam = 0.01 # Cross-sectional Area in m²
I_beam = 1e-5 # Moment of Inertia in m⁴
ops.section('Elastic', 1, E_beam, A_beam, I_beam)
# Define Elements
# Soil Elements (SSPquad) - Corrected Definitions
ops.element('SSPquad', 1, 1, 5, 6, 4, 1, 2.5, 5.0) # Element 1
ops.element('SSPquad', 2, 5, 2, 3, 6, 1, 2.5, 5.0) # Element 2
ops.element('SSPquad', 3, 6, 7, 8, 4, 1, 0.1, 10.0) # Element 3
# Beam Element (DispBeamColumn) - Correct Definition
ops.geomTransf('Linear', 1)
ops.beamIntegration('Legendre', 1, 1, 5) # beamIntTag=1, geomTransfTag=1, 5 Gauss points
ops.element('dispBeamColumn', 10, 6, 7, 1, 1) # Element 10 connects Node 6 to Node 7
# Apply Boundary Conditions
# Fix soil nodes at the base (all DOFs)
ops.fix(1, 1, 1, 1) # Node 1
ops.fix(2, 1, 1, 1) # Node 2
ops.fix(5, 1, 1, 1) # Node 5
# Fix top soil nodes (all DOFs to prevent rigid body motions)
ops.fix(3, 1, 1, 1) # Node 3
ops.fix(4, 1, 1, 1) # Node 4
# Fix rotational DOF for Node 6 (Contact Point and Beam Base)
ops.fix(6, 0, 0, 1) # Node 6: fix DOF3, free DOFs1 and 2
# Beam Top Node (7): Fix rotation, allow horizontal and vertical movement
ops.fix(7, 0, 0, 1) # Node 7: fix DOF3, free DOFs1 and 2
# Define Analysis Parameters
# Rayleigh Damping
alphaM = 0.0 # Mass proportional damping
betaK = 1e-2 # Increased stiffness proportional damping for better numerical stability
ops.rayleigh(alphaM, betaK, 0.0, 0.0)
# Define analysis settings
ops.constraints('Transformation') # Constraint handler
ops.numberer('RCM') # Node numbering
ops.system('UmfPack') # Alternative solver for better stability
ops.test('NormDispIncr', 1e-5, 200) # Increased number of iterations
ops.algorithm('ModifiedNewton') # Alternative algorithm for improved convergence
ops.integrator('LoadControl', 0.0005) # Smaller load step size for finer increments
ops.analysis('Static') # Analysis type
# Apply Loads
# Define load pattern
ops.timeSeries('Linear', 1)
ops.pattern('Plain', 1, 1)
ops.load(7, 100.0, 0.0, 0.0) # Apply horizontal load at Beam Top Node (7)
# Run Analysis
num_steps = 2000 # Increased number of steps for smaller load increments
analysis_result = ops.analyze(num_steps)
# Post-Processing and Visualization
if analysis_result == 0:
# Retrieve and print displacements
disp_beam_top = ops.nodeDisp(7)
disp_beam_base = ops.nodeDisp(6)
print(f'Beam Top Node Displacement (7): {disp_beam_top}')
print(f'Beam Base Node Displacement (6): {disp_beam_base}')
# Define a scale factor for visualization
disp_scale = 1000 # Adjust as needed for better visibility
# --- Plot Soil Deformation ---
soil_elements = [
[1, 5, 6, 4], # Element 1
[5, 2, 3, 6], # Element 2
[6, 7, 8, 4], # Element 3
]
plt.figure(figsize=(10, 8))
plt.title('Soil Elements Deformation')
for elem_nodes in soil_elements:
# Undeformed coordinates
x_undeformed = [ops.nodeCoord(node)[0] for node in elem_nodes + [elem_nodes[0]]]
y_undeformed = [ops.nodeCoord(node)[1] for node in elem_nodes + [elem_nodes[0]]]
# Deformed coordinates
x_deformed = [ops.nodeCoord(node)[0] + disp_scale * ops.nodeDisp(node)[0] for node in elem_nodes + [elem_nodes[0]]]
y_deformed = [ops.nodeCoord(node)[1] + disp_scale * ops.nodeDisp(node)[1] for node in elem_nodes + [elem_nodes[0]]]
# Plot undeformed and deformed shapes
plt.plot(x_undeformed, y_undeformed, 'b--', linewidth=1)
plt.plot(x_deformed, y_deformed, 'r-', linewidth=2)
plt.xlabel('X Coordinate (m)')
plt.ylabel('Y Coordinate (m)')
plt.legend(['Undeformed', 'Deformed'])
plt.grid(True)
plt.axis('equal')
plt.show()
# --- Plot Beam Deformation ---
beam_node_tags = [6, 7]
beam_coords_x = [ops.nodeCoord(node)[0] for node in beam_node_tags]
beam_coords_y = [ops.nodeCoord(node)[1] for node in beam_node_tags]
beam_disp_x = [ops.nodeCoord(node)[0] + disp_scale * ops.nodeDisp(node)[0] for node in beam_node_tags]
beam_disp_y = [ops.nodeCoord(node)[1] + disp_scale * ops.nodeDisp(node)[1] for node in beam_node_tags]
plt.figure(figsize=(10, 8))
plt.title('Beam Deformation')
plt.plot(beam_coords_x, beam_coords_y, 'b--', linewidth=1, label='Undeformed')
plt.plot(beam_disp_x, beam_disp_y, 'r-', linewidth=2, label='Deformed')
plt.xlabel('X Coordinate (m)')
plt.ylabel('Y Coordinate (m)')
plt.legend()
plt.grid(True)
plt.axis('equal')
plt.show()
else:
print('Analysis failed.')
ops.printModel()
i was expecting the displacement This is a very simple model. i want to make more user-interactive model but i am stuck on this only
Upvotes: 0
Views: 50