Elastic/plastic simple harmonic oscillator

elastic-plastic.py

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Function that returns dz/dt = [dx/dt, dv/dt]
def model(z, t):
    x = z[0]
    v = z[1]
    xp = z[2]
    
    stress = k*(x-xp) 

    def yield_curve(tau):
        lambd = (tau - mu * p) / eta
        if lambd < 0: lambd = 0
        return lambd
    
    plastic_strain_rate = yield_curve(stress)

    # For the simplified model, lambda is exactly the plastic strain
    dxp = plastic_strain_rate

    dxdt = v
    dvdt = -k * (x-xp) + loading(t)
    return [dxdt, dvdt, dxp]

def loading(t):
    # give the system a unit pulse between t0 and t1
    t0 = 1
    t1 = 2
    if (t>t0) & (t<t1): return 1
    else: return 0
    
    

# Initial conditions
x0  = 0.0  # initial displacement
v0  = 0.0  # initial velocity
xp0 = 0.0  # initial plastic velocity
z0 = [x0, v0, xp0]  # initial state vector

# Time points
t = np.linspace(0, 12, 1000)  # time array from 0 to 10 with 1000 points

# Parameters
k = 1  # spring constant, units stress per length
p = 1    # pressure, assumed constant
mu = 0.4 # friction coffficient
eta = 1  # viscoplastic viscosity

plt.figure(figsize=(10, 6))
for mu in (0.4,100):

    # Solve the ODE
    z = odeint(model, z0, t)

    # Extracting position and velocity
    x = z[:, 0]
    v = z[:, 1]

    plt.plot(t, x, label=f'mu={mu}')

load = []
for TT in t: load.append(loading(TT))
plt.plot(t,load,'-k',label='loading')
plt.xlabel('Time (s)')
plt.ylabel('Displacement')
plt.legend()
plt.grid(True)
plt.show()