In [1]:
%config InlineBackend.figure_format = 'svg'
A single degree of freedom system, with mass m = 500 kg, stiffness k = 72 kN m−1 and damping ratio ζ = 0.05 is at rest when it is subjected to an external force
where $p_0$ = 5.4 kN and $\omega_n ={}$ 12 rad s${}^\text{−1}$ is the undamped natural frequency of the dynamical system.
In [2]:
m = 500.0 ; k = 72000.0 ; z = 0.05 ; p0 = 5400.0 ; wn = 12.0 ;
c = 2*m*wn*z
wd = wn*sqrt(1.0-z*z)
a = wn/3
def p(t,p0=5400.0): return p0*exp(-a*t)
print 'Damping coefficient: ', c, 'N s/m'
We can find the exact response for t ≥ 0 using superposition of the general integral and a particular solution.
The particular integral for our loading must be $\xi(t)=G\exp(-at)$; substituting in the equation of motion and differentiating we have
solving for $G$ we have
In [3]:
G = p0/(a*a*m-a*c+k)
print 'Particular integral coefficient: ', G, 'm'
The general integral is given by the integral of the homogenous equation plus the particular integral we have just determined:
The system starts from initial rest conditions, i.e.,
Solving the two equations above for B and A gives
In [4]:
B = -G ; A = (a*G+z*wn*B)/wd
print 'Coefficients of the general integral'
print ' A: ', A, 'm'
print ' B: ', B, 'm'
In [5]:
def x_ex(t): return exp(-z*wn*t)*(A*sin(wd*t)+B*cos(wd*t))+G*exp(-a*t)
time = linspace(0.0, 3.0, 301)
plot(time,x_ex(time),label=r'$x(t)$',)
plot(time,G*exp(-a*time), label=r'$\xi(t)$')
xlabel('Time, s')
ylabel('Deflections, m')
legend(); grid()
In [6]:
h = 0.02
t_stop = 3.0
The linear acceleration algorithm solves the incremental equation of equilibrium where the velocity and acceleration increments, $\Delta v$ and $\Delta a$, are written in terms of the displacement increment $\Delta x$ and the initial conditions $v_0$, $a_0$.
This results in an equation of incremental equilibrium in $\Delta x$, with a modified load increment $\Delta p$ (depending on $v_0$ and $a_0$) and modified stiffness. These are the relevant coefficients:
In [7]:
a0_coef = 3*m + c*h/2
v0_coef = 3*c + 6*m/h
k_m_c = 3*c/h + 6*m/h/h
k_mod = k + 3*c/h + 6*m/h/h
We initialize the solution, and also a container for storing the iterates of the solution
In [8]:
container = []
t0 = 0.0 ; p0 = p(t0)
x0 = 0.0 ; v0 = 0.0
Now we propagate the solution
In [9]:
while 1:
a0 = (p0-c*v0-k*x0)/m
container.append([t0,p0,x0,v0,a0])
if t0 > t_stop: break
t1 = t0 + h
p1 = p(t1)
dp = (p1-p0) + v0_coef*v0 + a0_coef*a0
dx = dp/k_mod
dv = 3*dx/h - 3*v0 - a0*h/2
x0 = x0 + dx
v0 = v0 + dv
t0 = t1
p0 = p1
For simplicity in plotting we put the data collected in our container in a 2D data sructure,
In [10]:
tpxva = array(container)
Now it's time to plot our results
1. The numerical response, 0 < t < 3.0 s
2. A comparison of exact and numerical for 0 < t < 0.3 s
3. Another comparison, at the end of the interval, 2.7 s < t < 3.0 s
In [11]:
plot(tpxva[:,0], tpxva[:,2], label=r'$x(t)$, num.'); grid()
xlabel('Time, s') ; ylabel('Deflections, m')
dummy = axis(xmax=3.0)
In [12]:
plot(time, x_ex(time), label=r'$x(t)$, exact',)
plot(tpxva[:,0], tpxva[:,2], '.b', label=r'$x(t)$, numeric')
xlabel('Time, s') ; ylabel('Deflections, m')
axis(ymin=0.0,xmax=0.30001) ; grid()
In [13]:
plot(time, x_ex(time), label=r'$x(t)$, exact',)
plot(tpxva[:,0], tpxva[:,2], '.b', label=r'$x(t)$, numeric')
xlabel('Time, s') ; ylabel('Deflections, m')
axis(ymin=-0.00, ymax=+0.02, xmin=2.7, xmax=3.0) ; grid()
From the last figure, it is apparent that the numerical procedure, while really accurate in maxima estimation, introduces a slight period elongation in the SDOF system response.
We want to define a function that returns the tangent stiffness and the cumulated plastic deformation as a function of the current displacement and velocity and the cumulated plastic deformation itself.
Such a function depends on the elastic stiffness and the yield force, so here we have a function
that defines such a function in terms of k and fy.
In [14]:
def make_kt(k,fy):
"make_kt(k,fy) returns a function kt(u,v,up) returning kt, up"
u_y = fy/k
def kt(u,v,up, uy=u_y):
f=k*(u-up)
if (-fy)<f<fy: return k,up
if fy<=f and v>0: up=u-uy;return 0,up
if fy<=f and v<=0: up=u-uy;return k,up
if f<=(-fy) and v<0: up=u+uy;return 0,up
else: up=u+uy;return k,up
return kt
kt = make_kt(k, 6000.0)
In [15]:
def step(t0, p0, u0, v0, a0, k0, up, fy=6000.0, maxiter=30):
t1 = t0 + h
p1 = p(t1)
Dp = p1 - p0
Dp_= Dp + v0_coef*v0 + a0_coef*a0
k_ = k0 + k_m_c
# we prepare the machinery for the modified Newton-Raphson
# algorithm. if we have no state change in the time step, then the
# N-R algorithm is equivalent to the standard procedure
u_init=u0; v_init=v0 # initial state
f_spring=k*(u0-up) # the force in the spring
DR=Dp_ # the unbalanced force, initially equal to the
# external load increment
j = 0
while j < maxiter:
Du=DR/k_ # the disp increment according to the initial stiffness
u_next = u_init + Du
v_next = v_init + 3*Du/h - 3*v_init - a0*h/2
# we are interested in the total plastic elongation
oops,up=kt(u_next,v_next,up)
# because we need the spring force at the end
# of the time step
f_spring_next=k*(u_next-up)
# so that we can compute the fraction of the
# incremental force that's equilibrated at the
# end of the time step
df=f_spring_next-f_spring+(k_-k0)*Du
# and finally the incremental forces unbalanced
# at the end of the time step
DR=DR-df
# finish updating the system state
u_init=u_next; v_init=v_next; f_spring=f_spring_next
# if the unbalanced load is small enough (the
# criteria used in practical programs are
# energy based) exit the loop - note that we
# have no plasticization/unloading DR==0 at the
# end of the first iteration
if abs(DR)<fy*1E-6 :
break
else:
print "#", t0, j, u_next, DR/fy, up
j = j+1
# now the load increment is balanced by the spring force and
# increments in inertial and damping forces, we need to compute the
# full state at the end of the time step, and to change all
# denominations to reflect the fact that we are starting a new time step.
Du=u_init-u0
Dv=3*Du/h-3*v0-0.5*a0*h
u1=u0+Du ; v1=v0+Dv
#print t1, v1, v_init
return t1, p1, u1, v1
In [16]:
h = 0.005
a0_coef = 3*m + c*h/2
v0_coef = 3*c + 6*m/h
k_m_c = 3*c/h + 6*m/h/h
k_mod = k + 3*c/h + 6*m/h/h
We initialize the state of the system, note that the quantities that can be computed in terms of the current state are placed inside
the execution loop. We initialize also a container where we will store the current state and additional before advancing
the solution with a call to the step function
In [17]:
# initialization of the system state and a bit more
t0=0.0 ; p0 = p(t0)
u0=0.0 ; v0 = 0.0; up = 0.0
container = []
# solution loop
while True:
k0, up = kt(u0, v0, up)
a0=(p0-c*v0-k*(u0-up))/m
# print t0, p0, u0, v0, a0, k0, up
container.append((t0, p0, u0, v0, a0, k0, up))
if t0 > t_stop: break
t0, p0, u0, v0 = step( t0, p0, u0, v0, a0, k0, up)
We put all our data in an array and finally we plot our results. Note that the EP response is quite different from the elastic response, even if the maximum responses are almost identical, because in the EP case we had to put a relatively large permanent deformation in the spring to reach the maximum deflection and this permanent deformation has an evident effect on the (almost, the load is now sooo small...) free response.
In [18]:
tpxvaupkt=array(container)
uy=6000.0/72000.0 ; umax=uy+up
plot( tpxva[:,0], tpxva[:,2], label='EL (num)')
plot( tpxvaupkt[:,0], tpxvaupkt[:,2], label='EP (num)')
plot( (0.1,0.35), (uy,uy), '-k', label=r'$u_y$')
plot( (0.1,0.35), (umax,umax), ':r')
plot( (0.5,3.0), (up,up), ':r', label='$u_p$')
plot( (0.0,3.0), (0,0), '-k')
legend()
axis(xmax=3)
xlabel('Time, s') ; ylabel('Deflections, m')
grid()
In [19]:
plot(tpxvaupkt[:,0],tpxvaupkt[:,6],'-k')
plot(tpxvaupkt[:,0],tpxvaupkt[:,6],'.r')
axis(xmin=0.16, xmax=0.30)
xlabel('Time, s') ; ylabel('Cumulated plastic deformation, m')
grid()
In [20]:
plot(tpxvaupkt[:,0],tpxvaupkt[:,2], '-k', label='EP (num)')
plot(tpxva[:,0],tpxva[:,2], '-r', label='EL (num)')
axis(xmin=0.16,xmax=0.30,ymin=0.00) ; xlabel('Time, s') ; ylabel('Deflections, m'); legend(loc=5) ; grid()