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


```
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


```
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


```
# 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


```
x_el0(0), v_el0(0)
```




    (0.001, 0.025000000000000001)



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


```
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.


```
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.


```
# 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 <code>newton</code> function we find _t2_, the time (NB measured starting from _t1_) at which the velocity is equal to zero,


```
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.


```
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.


```
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

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

whose general integral is

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

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


```
# 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}$.


```
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() ;
```


```

```

3ddb1a6aa01f158de8cbe59dbb87fee4c418d44f
ipynb
Jupyter Notebook
dati_2013/03/Elasto-Plastic_Free_Vibrations.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
dati_2013/03/Elasto-Plastic_Free_Vibrations.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
dati_2013/03/Elasto-Plastic_Free_Vibrations.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
2019-06-23T12:32:39.000Z
2021-08-15T18:33:55.000Z
Qwen/Qwen-72B
1. YES
2. YES