Question 7.1: The water tank shown in Fig.7. 2a is subjected to the blast ......
The water tank shown in Fig.7. 2a is subjected to the blast loading illustrated in Fig 7.2b. Write a computer program in MATLAB to numerically evaluate the dynamic response of the tower by interpolation of the excitation. Plot the displacement u(t) and velocity \dot{u}(t) response in time interval 0 ≤ t < 0.5 s.
Assume W = 445.5 kN, k = 40913 kN/m, ρ = 0.05; F_0 = 445.5 kN, t_d== 0.05 s.
Use the step size as 0.005 s.
Calculate natural period.
u_{i+1}=a u_i+b \dot{u}_i+c F_i+ d F_{i+1} 7.13a
\dot{u}_{i+1}=a^{\prime} u_i+b^{\prime} \dot{u}_i+c^{\prime} F_i+d^{\prime} F_{i+1} 7.13b
Table 7.1 Coefficients in recurrence formula Eq. 7.13 (ρ < 1)
| a= e ^{-\rho \omega_n \Delta t}\left(\rho \sin \omega_d \Delta t+\sqrt{1-\rho^2} \cos \omega_d \Delta t\right) / \sqrt{\left(1-\rho^2\right)} |
| b= e ^{-\rho \omega_{n \Delta t}}\left(\frac{1}{\omega_{ d }} \sin \omega_d \Delta t\right) |
| c=\frac{1}{k \omega_n \Delta t}\left\{2 \rho+ e ^{-\rho \omega_n \Delta t} \omega_n \Delta t\left[\left(\frac{1-2 \rho^2}{\omega_d \Delta t}-\frac{\rho}{\sqrt{1-\rho^2}}\right) \sin \omega_d \Delta t\right.\right. \left.\left.-\left(1+\frac{2 \rho}{\omega_n \Delta t}\right) \cos \omega_d \Delta t\right]\right\} |
| d=\frac{1}{k \omega_n \Delta t}\left\{\omega_n \Delta t-2 \rho+ e ^{-\rho \omega_n \Delta t} \omega_n \Delta t\left[\left(\frac{2 \rho^2-1}{\omega_d \Delta t}\right) \sin \omega_d \Delta t+\left(\frac{2 \rho}{\omega_n \Delta t}\right) \cos \omega_d \Delta t\right]\right\} |
| a^{\prime}=- e ^{-\rho \omega_n \Delta t}\left(\frac{\omega_n}{\sqrt{1-\rho^2}} \sin \omega_d \Delta t\right) |
| b^{\prime}=e^{-\rho \omega_n \Delta t}\left(-\rho \sin \omega_d \Delta t+\sqrt{1-\rho^2} \cos \omega_d \Delta t\right) / \sqrt{1-\rho^2} |
| c^{\prime}=\frac{1}{k \Delta t}\left\{-1+ e ^{-\rho \omega_n \Delta t}\left[\left(\frac{\omega_{ n } \Delta t}{\sqrt{1-\rho^2}}+\frac{\rho}{\sqrt{1-\rho^2}}\right) \sin \omega_d \Delta t+\cos \omega_d \Delta t\right]\right\} |
| d^{\prime}=\frac{1}{k \Delta t}\left\{1- e ^{-\rho \omega_n \Delta t}\left[\frac{\rho}{\sqrt{1-\rho^2}} \sin \omega_d \Delta t+\cos \omega_d \Delta t\right]\right\} |
K = 40913 kN/m = 40931 × 10³ N/m
m = 445.5 × 1000/9.81 = 45412.8 kgm
T = 2π/ω = 2 π/30.02 = 0.209 s
Δt = 0.005< T/10 = 0.0209 s Hence o.k
The constants are calculated as
a = 0.9888 a′ = -4.4542
b = 0.0049 b′ = 0.9740
c = 1.82 e-10 c′ = 5.4196e-08
d = 9.1305e-011 d′ = 5.467e-08
Using the recurrence formulae u and \dot{u} can be calculated as shown in Table 7.2.
The displacement time history and the velocity time history are plotted as shown in Fig. 7.3 and the program is given below.
Program 7.1: MATLAB program for dynamic response of SDOF using recurrence formulae
Table 7.2 Displacement and velocity values at various times for Example 7.1
| \dot{u} | u | t | \dot{u} | u | t |
| 0.0794 | 0.0066 | 0.055 | 0 | 0 | 0 |
| 0.0479 | 0.0069 | 0.060 | 0.0461 | 0.0001 | 0.005 |
| 0.0158 | 0.0071 | 0.065 | 0.0856 | 0.0004 | 0.010 |
| -0.0161 | 0.0071 | 0.070 | 0.1177 | 0.0010 | 0.015 |
| -0.0473 | 0.0069 | 0.075 | 0.1419 | 0.0016 | 0.020 |
| -0.0769 | 0.0066 | 0.080 | 0.1577 | 0.0024 | 0.025 |
| -0.1094 | 0.0062 | 0.085 | 0.1648 | 0.0032 | 0.030 |
| -0.1291 | 0.0056 | 0.090 | 0.1634 | 0.0040 | 0.035 |
| -0.1506 | 0.0049 | 0.095 | 0.1534 | 0.0048 | 0.040 |
| -0.1684 | 0.0041 | 0.100 | 0.1353 | 0.0055 | 0.045 |
| 0.1096 | 0.0061 | 0.050 |
Script File
%***********************************************************
% DYNAMIC RESPONSE DUE TO EXTERNAL LOAD USING WILSON
RECURRENCE FORMULA
% **********************************************************
m=45412.8;
k=40913000;
wn=sqrt(k/m)
r=0.05;
u(1)=0;
v(1)=0;
tt=.50;
n=100;
n1=n+1
dt=tt/n;
td=.05;
jk=td/dt;
for m=1:n1
p(m)=0.0;
end
t=-dt
% **********************************************************
% ANY EXTERNAL LOADING VARIATION MUST BE DEFINED HERE
% **********************************************************
for m=1:jk+1;
t=t+dt;
p(m)=445500*(td-t)/td;
end
wd=wn*sqrt(1-r^2);
a=exp(-r*wn*dt)*(r*sin(wd*dt)/sqrt(1-r^2)+cos(wd*dt));
b=exp(-r*wn*dt)*(sin(wd*dt))/wd;
c2=((1-2*r^2)/(wd*dt)-r/sqrt(1-r^2))*sin(wd*dt)-(1+2*r/(wn*dt))*cos(wd*dt);
c=(1/k)*(2*r/(wn*dt)+exp(-r*wn*dt)*(c2));
d2=exp(-r*wn*dt)*((2.0*r^2-1)/(wd*dt)*sin(wd.*dt)+2.0*r/
(wn*dt)*cos(wd*dt));
d=(1/k)*(1-2.0*r/(wn*dt)+d2);
ad=-exp(-r*wn*dt)*wn*sin(wd*dt)/(sqrt(1-r^2));
bd=exp(-r*wn*dt)*(cos(wd*dt)-r*sin(wd*dt)/sqrt(1-r^2));
c 1 = e x p ( - r * w n * d t ) * ( ( w n / s q r t ( 1 - r ^ 2 ) + r / ( d t * s q r t ( 1 -
r^2)))*sin(wd*dt)+cos(wd*dt)/dt);
cd=(1/k)*(-1/dt+c1);
d1=exp(-r*wn*dt)*(r*sin(wd*dt)/sqrt(1-r^2)+cos(wd*dt));
dd=(1/(k*dt))*(1-d1);
for m=2:n1
u(m)=a*u(m-1)+b*v(m-1)+c*p(m-1)+d*p(m);
v(m)=ad*u(m-1)+bd*v(m-1)+cd*p(m-1)+dd*p(m);
end
for m=1:n1
s(m)=(m-1)*dt
end
figure(1);
plot(s,u,‘k’);
xlabel(‘ time (t) in seconds’)
ylabel(‘ Response displacement u in m’)
title(‘ dynamic response’)
figure(2);
plot(s,v,‘k’);
xlabel(‘ time (t) in seconds’)
ylabel(‘ Response velocity v in m/sec’)
title(‘ dynamic response’)