Mayank Ratnakar
Mayank Ratnakar

Reputation: 1

I have been working on a slope stability model and using opensees and when i seem to make any progress in my i code my kernal keeps dying

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

Answers (0)

Related Questions