We want to compute the response, using the constant acceleration algorithm plus MNR, of an Elasto Plastic (EP) system... but how we can confirm or reject our results?
It turns out that computing the exact response of an EP system with a single degree of freedom is relatively simple.
Here we discuss a program that computes the analytical solution of our problem.
The main building blocks of the program will be two functions that compute, for the elastic phase and for the plastic phase, the analytical functions that give the displacement and the velocity as functions of time.
We are definining a function that, for a linear dynamic system, returns not the displacement or the velocity at a given time, but rather a couple of functions of time that we can use afterwards to compute displacements and velecities at any time of interest.
The response depends on the parameters of the dynamic system $m,c,k,$ on the initial conditions $x_0, v_0,$ and on the characteristics of the external load.
Here the external load is limited to a linear combination of a cosine modulated, a sine modulated (both with the same frequency $\omega$) and a constant force,
but that's all that is needed for the present problem.
The particular integral being
substituting in the equation of motion and equating all the corresponding terms gives the undetermined coefficients in $\xi(t)$, then evaluation of the general integral and its time derivative for $t=0$ permits to find the constants in the homogeneous part of the integral.
The final step is to define the displacement and the velocity function, according to the constants we have determined, and to return these two function to the caller
In [267]:
def resp_elas(m,c,k, cC,cS,w, F, x0,v0):
wn2 = k/m ; wn = sqrt(wn2) ; beta = w/wn
z = c/(2*m*wn)
wd = wn*sqrt(1-z*z)
# xi(t) = R sin(w t) + S cos(w t) + D
det = (1.-beta**2)**2+(2*beta*z)**2
R = ((1-beta**2)*cS + (2*beta*z)*cC)/det/k
S = ((1-beta**2)*cC - (2*beta*z)*cS)/det/k
D = F/k
A = x0-S-D
B = (v0+z*wn*A-w*R)/wd
def x(t):
return exp(-z*wn*t)*(A*cos(wd*t)+B*sin(wd*t))+R*sin(w*t)+S*cos(w*t)+D
def v(t):
return (-z*wn*exp(-z*wn*t)*(A*cos(wd*t)+B*sin(wd*t))
+wd*exp(-z*wn*t)*(B*cos(wd*t)-A*sin(wd*t))
+w*(R*cos(w*t)-S*sin(w*t)))
return x,v
In this case the equation of motion is
the homogeneous response is
and the particular integral, for a load described as in the previous case, is again
Having computed $R, S,$ and $D$ from substituting $\xi$ in the equation of motion, $A$ and $B$ by imposing the initial conditions,we can define the displacement and velocity functions and, finally, return these two functions to the caller.
In [268]:
def resp_yield(m,c, cC,cS,w, F, x0,v0):
# csi(t) = R sin(w t) + S cos(w t) + Q t
Q = F/c
det = w**2*(c**2+w**2*m**2)
R = (+w*c*cC-w*w*m*cS)/det
S = (-w*c*cS-w*w*m*cC)/det
# x(t) = A exp(-c t/m) + B + R sin(w t) + S cos(w t) + Q t
# v(t) = - c A/m exp(-c t/m) + w R cos(w t) - w S sin(w t) + Q
#
# v(0) = -c A / m + w R + Q = v0
A = m*(w*R + Q - v0)/c
# x(0) = A + B + S = x0
B = x0 - A - S
def x(t):
return A*exp(-c*t/m)+B+R*sin(w*t)+S*cos(w*t)+Q*t
def v(t):
return -c*A*exp(-c*t/m)/m+w*R*cos(w*t)-w*S*sin(w*t)+Q
return x,v
We need to find when
to individuate the three ranges of different behaviour
We can use the simple and robust algorithm of bisection to find the roots for
In [269]:
def bisect(f,val,x0,x1):
h = (x0+x1)/2.0
fh = f(h)-val
if abs(fh)<1e-8 : return h
f0 = f(x0)-val
if f0*fh > 0 :
return bisect(f, val, h, x1)
else:
return bisect(f, val, x0, h)
In [270]:
mass = 1000. # kg
k = 40000. # N/m
zeta = 0.03 # damping ratio
fy = 2500. # N
In [271]:
damp = 2*zeta*sqrt(k*mass)
xy = fy/k # m
In [272]:
t1 = 0.3 # s
w = pi/t1 # rad/s
Po = 6000. # N
In [273]:
x0=0.0 # m
v0=0.0 # m/s
x_next, v_next = resp_elas(mass,damp,k, 0.0,Po,w, 0.0, x0,v0)
In [274]:
t_yield = bisect(x_next, xy, 0.0, t1)
print t_yield, x_next(t_yield)*k
In [275]:
t_el = linspace( 0.0, t_yield, 201)
x_el = vectorize(x_next)(t_el)
v_el = vectorize(v_next)(t_el)
# ------------------------------
figure(0);grid()
plot(t_el,x_el,
(0,0.25),(xy,xy),'--b',
(t_yield,t_yield),(0,0.0699),'--b')
title("$x_{el}(t)$")
xlabel("Time, s")
ylabel("Displacement, m")
# ------------------------------
figure(1);grid()
plot(t_el,v_el)
title("$\dot x_{el}(t)$")
xlabel("Time, s")
ylabel("Velocity, m/s")
Out[275]:
In [276]:
x0=x_next(t_yield)
v0=v_next(t_yield)
print x0, v0
now, the load must be expressed in function of a restarted time,
In [277]:
cS = Po*cos(w*t_yield)
cC = Po*sin(w*t_yield)
print Po*sin(w*0.55), cS*sin(w*(0.55-t_yield))+cC*cos(w*(0.55-t_yield))
Now we generate the displacement and velocity functions for the yielded phase, please note that the yielded spring still exerts a constant force $f_y$ on the mass, and that this fact must be (and it is) taken into account.
In [278]:
x_next, v_next = resp_yield(mass, damp, cC,cS,w, -fy, x0,v0)
At this point I must confess that I have already peeked the numerical solution, hence I know that the velocity at $t=t_1$ is grater than 0 and I know that the current solution is valid in the interval $t_y\leq t\leq t_1$.
In [279]:
t_y1 = linspace(t_yield, t1, 101)
x_y1 = vectorize(x_next)(t_y1-t_yield)
v_y1 = vectorize(v_next)(t_y1-t_yield)
In [280]:
figure(3) ; grid()
plot(t_el,x_el, t_y1,x_y1,
(0,0.25),(xy,xy),'--b',
(t_yield,t_yield),(0,0.0699),'--b')
xlabel("Time, s")
ylabel("Displacement, m")
# ------------------------------
figure(4) ; grid()
plot(t_el, v_el, t_y1, v_y1)
xlabel("Time, s")
ylabel("Velocity, m/s")
Out[280]:
In the next phase, still it is $\dot x> 0$ so that the spring is still yielding, but now $p(t)=0$, so we must compute two new state functions, starting as usual from the initial conditions (note that the yielding force is still applied)
In [281]:
x0 = x_next(t1-t_yield)
v0 = v_next(t1-t_yield)
print x_0, v_0
x_next, v_next = resp_yield(mass, damp, 0, 0, w, -fy, x0, v0)
t2 = t1 + bisect( v_next, 0.0, 0, 0.3)
print t2
t_y2 = linspace( t1, t2, 101)
x_y2 = vectorize(x_next)(t_y2-t1)
v_y2 = vectorize(v_next)(t_y2-t1)
print x_next(t2-t1)
# ------------------------------
figure(5) ; grid()
plot(t_el,x_el, t_y1,x_y1, t_y2, x_y2,
(0,0.25),(xy,xy),'--b',
(t_yield,t_yield),(0,0.0699),'--b')
xlabel("Time, s")
ylabel("Displacement, m")
# ------------------------------
figure(6) ; grid()
plot(t_el, v_el, t_y1, v_y1, t_y2, v_y2)
xlabel("Time, s")
ylabel("Velocity, m/s")
Out[281]:
The only point worth commenting is the constant force that we apply to our system.
The force-displacement relationship for an EP spring is
taking the negative, constant part of the last expression into the right member of the equation of equilibrium we have a constant term, as follows
In [282]:
x0 = x_next(t2-t1) ; v0 = 0.0
print x0, x_y2[-1]
x_next, v_next = resp_elas(mass,damp,k, 0.0,0.0,w, k*x0-fy, x0,v0)
t_e2 = linspace(t2,4.0,201)
x_e2 = vectorize(x_next)(t_e2-t2)
v_e2 = vectorize(v_next)(t_e2-t2)
In [283]:
# ------------------------------
figure(7) ; grid()
plot(t_el, x_el, '-b',
t_y1, x_y1, '-r',
t_y2, x_y2, '-r',
t_e2, x_e2, '-b',
(0.6, 4.0), (x0-xy, x0-xy), '--y')
title("In blue: elastic phases.\n"+
"In red: yielding phases.\n"+
"Dashed: permanent plastic deformation.")
xlabel("Time, s")
ylabel("Displacement, m")
Out[283]:
In [284]:
def p(t):
if t<0.3:
return 6000.0*sin(pi*t/0.3)
return 0.0
In [284]: