The parameters of the system, m, k, z, and the yielding force are given, all the remaining quantities are computed.



In [1]:

    
mass = 16e3     
stif = 2.5e6
zeta = 0.05

fy   = 5e3
uy   = fy/stif

wn   = sqrt(stif/mass)
wd   = wn*sqrt(1-zeta**2)
damp = 2*wn*mass*zeta
print mass, stif, fy, uy, wn, wd, damp









    



16000.0 2500000.0 5000.0 0.002 12.5 12.4843652221 20000.0

The initial conditions are assigned



In [2]:

    
x0 =  1e-3
v0 = 25e-3

Next, the constants of integration such that the homogeneous solution satisfies the initial conditions and the definition of a displacement function and a velocity function valid in the elastic phase



In [3]:

    
# x(0) = A ; v(0) = - zeta wn A + wd B => v(0) + zeta wn x(0) = wd B
A = x0 ; B = (v0+zeta*wn*A)/wd
def x_el0(t, A=A, B=B): return exp(-zeta*wn*t)*(A*cos(wd*t)+B*sin(wd*t))
def v_el0(t, A=A, B=B): return exp(-zeta*wn*t)*((B*wd-A*wn*zeta)*cos(wd*t)-(B*wn*zeta+A*wd)*sin(wd*t))

Just to be sure, let's verify that the initial conditions are OK



In [4]:

    
x_el0(0), v_el0(0)









    Out[4]:





(0.001, 0.025000000000000001)

Now, let's plot our elastic response and take note that $x_\text{max}>2\\,$mm${}=u_\text{y}$.



In [5]:

    
t = linspace(0,0.5, 1001)
plot(t*1000,1000*x_el0(t))
grid();xlabel('t/ms');ylabel('x/mm');

The solution that we have is valid only for $x \le u_\text{y}$, so we have to find the time $t_1$ such that $x=u_\text{y}$, and then find $x_1$ and $v_1$, the last valid values for our current solution and the initial values for our next solution.



In [6]:

    
from scipy.optimize import newton
t1 = newton(lambda x:x_el0(x)-.002, 0.0)
x1 = uy
v1 = v_el0(t1)

In the plastic range, until unloading (i.e., for $v\ge0$), the equation of motion is

\begin{align}m\ddot x + c\dot x + f_\text{S}=0\end{align}

and being $f_\text{S}=f_\text{y}$ we have

$$m\ddot x + c\dot x = -f_\text{y};$$

the particular integral is $\xi = -f_\text{y} t/c$ and the general integral can be written

$$x(t) = A \exp(-ct/m) -f_\text{y}t/c + B,$$

where, as usual, the constants of integraton are determined by the initial conditions.

Knowing the constants of integration, you can define the functions that give the displacement and the velocity in the interval from yielding to unloading.



In [7]:

    
# x2(t-t1) = x2(t2) = A exp(-c t2 / m) - fy t2 / c + B ; v2(t2) = - c A / m exp(...) - fy/c
# x2(0) = A+B = x1 ; v2(0) = - A c/m -fy/c = v1
A = - ( v1 + fy/damp ) * mass / damp
B = x1 - A
print A, B
def x_pl(t, A=A, B=B): return A*exp(-damp*t/mass)-fy*t/damp+B
def v_pl(t, A=A, B=B): return -damp*A*exp(-damp*t/mass)/mass -fy/damp









    



-0.208312097754 0.210312097754

Using the newton function we find t2, the time (NB measured starting from t1) at which the velocity is equal to zero,



In [8]:

    
t2 = newton(v_pl,t1)

so that we can compute the position and the velocity at the end of the plastic phase or, equivalently, at the beginning of the unloading. Also, in the following we'll need the total plastic deformation, so it is computed as the difference between the position at the end and the position at the beginning of the plastic phase.



In [9]:

    
x2 = x_pl(t2)
v2 = v_pl(t2)
dp = x2-x1

To plot the first two phases of the response we have to be careful with time ranges, but it;s relatively simple. In blue the elastic, loading phase, in green the plastic phase.



In [10]:

    
t_1 = linspace(0,t1,1001)
t_2 = linspace(0,t2,1001)
plot(t_1,x_el0(t_1), t_2+t1, x_pl(t_2)) ; grid() ;

Now for the elastic unloading, the elastic force is now $f_\text{S}=k(x-\Delta_\text{pl})$ so that the equation of motion can be written

$m\ddot x + c\dot x + kx = k\Delta_\text{pl}$

whose general integral is

$x(t) = \exp(-\zeta\omega_n t)\left(A\cos(\omega_\text{D}t)+B\sin(\omega_\text{D}\right)+\Delta_\text{pl}.$

Next, by the initial conditions we have the constants of integration, and finally the definition of the displacements and velocities in the unloading phase.



In [11]:

    
# x(t3) = exp A cos + exp B sin + dp ; v(t3) = ...
# x3(0) = A + dp = x2 ; v3(0) = B wd - A z wn = v2
A = x2-dp
B = (v2+A*zeta*wn)/wd

def x_el1(t, A=A, B=B): return exp(-zeta*wn*t)*(A*cos(wd*t)+B*sin(wd*t))+dp
def v_el1(t, A=A, B=B): return exp(-zeta*wn*t)*((B*wd-A*wn*zeta)*cos(wd*t)-(B*wn*zeta+A*wd)*sin(wd*t))

Eventually, we can plot the free vibrations of our elastic-perfectly plastic system... note that I've stretched the time axis to let the response damp out almost completely, so that you can readily appreciate that the system is not oscillating around the initial equilibrium position, i.e., zero, but is asymptotically reaching a new position of equilibrium, $\Delta_\text{pl}$.



In [12]:

    
figsize(18,6)
t_3 = linspace(0,8-t1-t2,2001)
plot(t_1,x_el0(t_1), t_2+t1, x_pl(t_2), t1+t2+t_3, x_el1(t_3)) ; grid() ;



In [12]: