Half-car model with CoupledZeroLength Elements

[bennuDJ]

Hello,

I’m working on building a half-car model using CoupledZeroLength elements. The model has 4 degrees of freedom (DoF), including the vertical displacements of the individual masses and the lateral rotation about the vehicle’s center of mass.

https://imgur.com/a/hn5qoax

I’ve attached the code I’ve developed so far. Unfortunately, it doesn’t produce the correct solution. Could you please provide guidance on how to properly model this half-car system?

The initial goal is to test the model's performance under a standard static load. Once this is validated, I plan to conduct a transient analysis by applying multi-support ground-imposed displacements to simulate road input.

Thank you!

Code: Select all

import openseespy.opensees as ops
import matplotlib.pyplot as plt
import opsvis as opsv

# Parameters
kt = 175.e4  # Tire stiffness (N/m)
ks = 6.e6    # Suspension stiffness (N/m)
cs = 1.e4    # Suspension damping (N*s/m)

mu = 900     # Unsprung mass (kg)
mb = 10500   # Sprung mass (kg)
ib = 50000   # Sprung mass moment of inertia (kg*m^2)
g = 9.81     # Gravity (m/s^2)

# Total load due to gravity
N = (mb + 2 * mu) * g  # Static load

# Define the OpenSees model
ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3)  # 2D model with 3 DOFs per node

# Define nodes (X, Y positions)
# Base nodes (unsprung masses)
ops.node(0, 0, 0)       # Front axle base
ops.fix(0, 1, 1, 1)     # Fixed

ops.node(100, 5, 0)     # Rear axle base
ops.fix(100, 1, 1, 1)   # Fixed

# Unsprung masses (tires)
ops.node(1, 0, 0)       # Front unsprung mass
ops.node(101, 5, 0)     # Rear unsprung mass

# Sprung mass (rigid body)
ops.node(2, 0, 0)       # Front connection to body
ops.node(102, 5, 0)     # Rear connection to body
ops.node(10, 2.5, 0)    # Center of mass of the body
ops.fix(2, 1, 0, 0)    # Front connection to body free in Y and rotation
ops.fix(102, 1, 0, 0)    # Rear connection to body free in Y and rotation
ops.fix(10, 1, 0, 0)    # Center mass free in Y and rotation

# Define mass
ops.mass(10, 0, mb, ib)

ops.mass(1, 0, mu, 0)
ops.mass(101, 0, mu, 0)


# Define materials
ops.uniaxialMaterial('Elastic', 1, kt)       # Tire stiffness
ops.uniaxialMaterial('Elastic', 2, ks, cs)  # Suspension stiffness

# Define elements
# Front axle
ops.element('CoupledZeroLength', 1, 0, 1, 1, 2, 1)  # Tire and suspension
ops.element('CoupledZeroLength', 2, 1, 2, 1, 2, 2)  # Suspension to sprung mass

# Rear axle
ops.element('CoupledZeroLength', 3, 100, 101, 1, 2, 1)  # Tire and suspension
ops.element('CoupledZeroLength', 4, 101, 102, 1, 2, 2)  # Suspension to sprung mass

# Rigid diaphragm to connect front and rear sprung masses
ops.rigidDiaphragm(1, 10, *[2, 102])

# Visualization of the model
opsv.plot_model()
plt.show()

# Static analysis setup
# Create TimeSeries
ops.timeSeries("Linear", 1)

# Create a plain load pattern
ops.pattern("Plain", 1, 1)

# Apply static vertical load at the center of mass
ops.load(10, 0, -N, 0)  # Apply vertical load (-N)

# Static analysis configuration
ops.test('NormDispIncr', 1e-6, 10, 0)
ops.algorithm("KrylovNewton")       # Linear algorithm
ops.system("UmfPack")         # System of equations
ops.constraints("Transformation")      # Constraints
ops.integrator("LoadControl", 1.0)  # Load control with step size 1.0
ops.analysis("Static")        # Static analysis

ops.analyze(1)

[mhscott]

I don't think the 'rigidDiaphragm' is doing what you want. Try using 'rigidLink -beam' instead.

https://portwooddigital.com/2023/04/02/ ... nto-a-bar/

[bennuDJ]

Dear Prof. Scott,

Thank you for your response. Using a rigid beam to connect the center of mass to the oscillators provides an appropriate solution. Attached is the working code.

Code: Select all

import openseespy.opensees as ops
import opsvis as opsv
import matplotlib.pyplot as plt
import numpy as np

# =============================================================================
# Material Definition
# =============================================================================
ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3)

# Material 1 - Tire stiffness
k1 = 175.e4  # Stiffness for material 1
ops.uniaxialMaterial('Elastic', 1, k1)

# =============================================================================
# Model Generation
# =============================================================================
# Geometry parameters
axle_space = 5  # Axle spacing (inches)

# Node and Element Definitions
# Row 1 (Left axle)
ops.node(0, 0, 0)  # Fixed node for row 1
ops.fix(0, 1, 1, 1)
ops.node(1, 0, 0)  # Free node for row 1
ops.element('zeroLength', 1, 0, 1, '-mat', 1, 1, '-dir', 1, 2)  # Tire element between node 0 and node 1

# Row 2 (Right axle)
ops.node(100, -axle_space, 0)  # Fixed node for row 2
ops.fix(100, 1, 1, 1)
ops.node(101, -axle_space, 0)  # Free node for row 2
ops.element('zeroLength', 2, 100, 101, '-mat', 1, 1, '-dir', 1, 2)  # Tire element between node 100 and node 101

# Center Node
ops.node(10, -axle_space * 0.5, 0)
ops.mass(10, 0, 10500, 50000)

# Rigid Links
ops.rigidLink('-beam', 10, 1)
ops.rigidLink('-beam', 10, 101)

# Visualize Model
opsv.plot_model()
plt.show()

# =============================================================================
# Transient Analysis Setup
# =============================================================================
# Time Parameters
total_time = 2  # seconds
time_steps = 1000
time = np.linspace(0, total_time, time_steps)  # Time vector
dt = time[1] - time[0]

# Road Disturbance Profiles
def disturbance_function(t, freq, start, end):
    v = 0.025 * (1 - np.cos(2 * np.pi * freq * t))
    v[(t < start) | (t > end)] = 0
    return v.tolist()

freq = 4  # Frequency in Hz
v1 = disturbance_function(time, freq, 2 / freq, 3 / freq)
v2 = disturbance_function(time, freq, 3 / freq, 4 / freq)

# Time Series for Ground Motion
ops.timeSeries('Path', 1, '-dt', dt, '-values', *v2)
ops.timeSeries('Path', 2, '-dt', dt, '-values', *v1)

# Apply Ground Motion
ops.pattern('MultipleSupport', 1)
ops.groundMotion(1, 'Plain', '-disp', 1)
ops.imposedMotion(0, 2, 1)  # Apply to node 0 (DOF 2 - y-direction)
ops.groundMotion(2, 'Plain', '-disp', 2)
ops.imposedMotion(100, 2, 2)  # Apply to node 100 (DOF 2 - y-direction)

# =============================================================================
# Analysis Settings
# =============================================================================
ops.constraints('Transformation')
ops.numberer('RCM')
ops.system('UmfPack')
ops.test('NormDispIncr', 1e-6, 10)
ops.algorithm('Newton')
ops.integrator('Newmark', 0.5, 0.25)
ops.analysis('Transient')

# Perform Transient Analysis
node_1_disp, node_101_disp = [], []
node_0_react, node_100_react = [], []

for step in range(len(time)):
    ok = ops.analyze(1, dt)
    if ok != 0:
        print(f"Analysis failed at step {step}, trying ModifiedNewton algorithm.")
        ops.algorithm('ModifiedNewton')
        ok = ops.analyze(1, dt)
        if ok != 0:
            print(f"Analysis failed at step {step} after retrying.")
            break
    ops.reactions()
    node_0_react.append(ops.nodeReaction(0, 2))
    node_100_react.append(ops.nodeReaction(100, 2))
    node_1_disp.append(ops.nodeDisp(1, 2))
    node_101_disp.append(ops.nodeDisp(101, 2))

# =============================================================================
# Post-Processing
# =============================================================================
# Plot Road Disturbance
plt.figure(figsize=(10, 6))
plt.plot(time, v1, label='Road Disturbance 1', color='b')
plt.plot(time, v2, label='Road Disturbance 2', color='r')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.title('Road Disturbance Profile')
plt.legend()
plt.grid(True)
plt.show()

# Plot Node Displacements
plt.figure(figsize=(10, 6))
plt.plot(time, node_1_disp, label='Node 1 Displacement', color='g')
plt.plot(time, node_101_disp, label='Node 101 Displacement', color='r')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.title('Node Displacements')
plt.legend()
plt.grid(True)
plt.show()

# Plot Node Reactions
plt.figure(figsize=(10, 6))
plt.plot(time, node_0_react, label='Node 0 Reaction', color='g')
plt.plot(time, node_100_react, label='Node 100 Reaction', color='r')
plt.xlabel('Time (s)')
plt.ylabel('Reaction (N)')
plt.title('Node Reactions')
plt.legend()
plt.grid(True)
plt.show()

However, I am encountering another issue. In this case, I can only use one "zeroLength" element to form a single axis and connect the axis components free nodes via rigid link beam. When I attempt to connect two "zeroLength" elements for a single axis (to simulate sprung and unsprung masses) and link them together, the analysis fails, and I receive the following error message:

Analysis failed at step 0, trying ModifiedNewton algorithm.
Analysis failed at step 0 after retrying.
WARNING: numeric analysis returns 1 -- Umfpackgenlinsolver::solve
WARNING AcceleratedNewton::solveCurrentStep() - the LinearSysOfEqn failed in solve()
DirectIntegrationAnalysis::analyze() - the Algorithm failed at time 0.002002
OpenSees > analyze failed, returned: -3 error flag
WARNING: numeric analysis returns 1 -- Umfpackgenlinsolver::solve
WARNING ModifiedNewton::solveCurrentStep() - the LinearSysOfEqn failed in solve()
DirectIntegrationAnalysis::analyze() - the Algorithm failed at time 0.002002
OpenSees > analyze failed, returned: -3 error flag

I followed this example (https://portwooddigital.com/2023/11/11/ ... er-models/) for linking the elements (which works for a 1-dimensional element), but it does not function in this case. Could you please suggest how I might resolve this issue?

Here is the non-working code that results with an error:

Code: Select all

import openseespy.opensees as ops
import opsvis as opsv
import matplotlib.pyplot as plt
import numpy as np

# =============================================================================
# Material Definition
# =============================================================================
ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3)

# Material 1 - Tire stiffness
k1 = 175.e4  # Stiffness for material 1
ops.uniaxialMaterial('Elastic', 1, k1)
k2 = 6.e6  # Stiffness for material 1
ops.uniaxialMaterial('Elastic', 2, k2)

# =============================================================================
# Model Generation
# =============================================================================
# Geometry parameters
axle_space = 5  # Axle spacing (inches)

# Node and Element Definitions
# Row 1 (Left axle)
ops.node(0, 0, 0)  # Fixed node for row 1
ops.fix(0, 1, 1, 1)
ops.node(1, 0, 0)  # Free node for row 1 (unsprung)
ops.mass(1, 0, 900, 0)
ops.node(2, 0, 0) # Free node for row 1 (sprung)

ops.element('zeroLength', 1, 0, 1, '-mat', 1, 1, '-dir', 1, 2)  # Tire element between node 0 and node 1
ops.element('zeroLength', 2, 1, 2, '-mat', 2, 2, '-dir', 1, 2)  # Suspension element between node 1 and node 0

# Row 2 (Right axle)
ops.node(100, -axle_space, 0)  # Fixed node for row 2
ops.fix(100, 1, 1, 1)
ops.node(101, -axle_space, 0)  # Free node for row 2 (unsprung)
ops.mass(101, 0, 900, 0)
ops.node(102, -axle_space, 0) # Free node for row 1 (sprung)

ops.element('zeroLength', 3, 100, 101, '-mat', 1, 1, '-dir', 1, 2)  # Tire element between node 100 and node 101
ops.element('zeroLength', 4, 101, 102, '-mat', 2, 2, '-dir', 1, 2)  # Tire element between node 100 and node 101


# Center Node
ops.node(10, -axle_space * 0.5, 0)
ops.mass(10, 0, 10500, 50000)

# Rigid Links
ops.rigidLink('-beam', 10, 2)
ops.rigidLink('-beam', 10, 102)

# Visualize Model
opsv.plot_model()
plt.show()

# =============================================================================
# Transient Analysis Setup
# =============================================================================
# Time Parameters
total_time = 2  # seconds
time_steps = 1000
time = np.linspace(0, total_time, time_steps)  # Time vector
dt = time[1] - time[0]

# Road Disturbance Profiles
def disturbance_function(t, freq, start, end):
    v = 0.025 * (1 - np.cos(2 * np.pi * freq * t))
    v[(t < start) | (t > end)] = 0
    return v.tolist()

freq = 4  # Frequency in Hz
v1 = disturbance_function(time, freq, 2 / freq, 3 / freq)
v2 = disturbance_function(time, freq, 3 / freq, 4 / freq)

# Time Series for Ground Motion
ops.timeSeries('Path', 1, '-dt', dt, '-values', *v2)
ops.timeSeries('Path', 2, '-dt', dt, '-values', *v1)

# Apply Ground Motion
ops.pattern('MultipleSupport', 1)
ops.groundMotion(1, 'Plain', '-disp', 1)
ops.imposedMotion(0, 2, 1)  # Apply to node 0 (DOF 2 - y-direction)
ops.groundMotion(2, 'Plain', '-disp', 2)
ops.imposedMotion(100, 2, 2)  # Apply to node 100 (DOF 2 - y-direction)

# =============================================================================
# Analysis Settings
# =============================================================================
ops.constraints('Transformation')
ops.numberer('RCM')
ops.system('UmfPack')
ops.test('NormDispIncr', 1e-6, 10)
ops.algorithm('Newton')
ops.integrator('Newmark', 0.5, 0.25)
ops.analysis('Transient')

# Perform Transient Analysis
node_1_disp, node_101_disp = [], []
node_2_disp, node_102_disp = [], []
node_0_react, node_100_react = [], []

for step in range(len(time)):
    ok = ops.analyze(1, dt)
    if ok != 0:
        print(f"Analysis failed at step {step}, trying ModifiedNewton algorithm.")
        ops.algorithm('ModifiedNewton')
        ok = ops.analyze(1, dt)
        if ok != 0:
            print(f"Analysis failed at step {step} after retrying.")
            break
    ops.reactions()
    node_0_react.append(ops.nodeReaction(0, 2))
    node_100_react.append(ops.nodeReaction(100, 2))
    node_1_disp.append(ops.nodeDisp(1, 2))
    node_101_disp.append(ops.nodeDisp(101, 2))
    node_2_disp.append(ops.nodeDisp(2, 2))
    node_102_disp.append(ops.nodeDisp(102, 2))

# =============================================================================
# Post-Processing
# =============================================================================
# Plot Road Disturbance
plt.figure(figsize=(10, 6))
plt.plot(time, v1, label='Road Disturbance 1', color='b')
plt.plot(time, v2, label='Road Disturbance 2', color='r')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.title('Road Disturbance Profile')
plt.legend()
plt.grid(True)
plt.show()

# Plot Node Displacements
plt.figure(figsize=(10, 6))
plt.plot(time, node_1_disp, label='Node 1 Displacement', color='g')
plt.plot(time, node_101_disp, label='Node 101 Displacement', color='r')
plt.plot(time, node_2_disp, label='Node 2 Displacement', color='g', linestyle=':')
plt.plot(time, node_102_disp, label='Node 102 Displacement', color='r', linestyle=':')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.title('Node Displacements')
plt.legend()
plt.grid(True)
plt.show()

# Plot Node Reactions
plt.figure(figsize=(10, 6))
plt.plot(time, node_0_react, label='Node 0 Reaction', color='g')
plt.plot(time, node_100_react, label='Node 100 Reaction', color='r')
plt.xlabel('Time (s)')
plt.ylabel('Reaction (N)')
plt.title('Node Reactions')
plt.legend()
plt.grid(True)
plt.show()

Thank you!

[bennuDJ]

I’ve found a solution by generating an additional point at the center of mass (Point 11) and connecting it rigidly to Points 1 and 101, where the unsprung mass of the wheel is located.