CE 521 - Structural Dynamics

Instructor: Prof. Wang, Ph.D, P.E

Homework 6
(Numerical Integration)

1. Given: A water tower (SDF system) has the following dynamic properties as shown in Figure
(a). It is subjected to the half-sine-wave loading shown in Figure (b). Neglect damping and
assume that the tower is initiated at rest.

Required:
1) Using the Trapezoidal Method, determine the displacement response history of the system
to the given half-cycle sine pulse force for the first 3.0 second by using the time interval
∆t = 0.1 second. Include your MATLAB code in your submission.
2) Repeat (1) by using the time interval ∆t = 0.05 second.
3) Evaluate the theoretical solution. Compare with the numerical solution in 1) and 2) above.
Display the difference between your numerical answer and analytical answer.
4) Using the Simpson’s Method, determine the displacement response history of the system to
the given half-cycle sine pulse force for the first 2.0 second by using the time interval
∆t = 0.1 second.
5) Submit a figure using MATLAB containing all four plots as required above.
>> % MATLAB Code for Water Tower Displacement Response to Half-Sine
Wave Force

% Given Parameters
w = 978.8; % Weight of the tower in kips
k = 100; % Stiffness in kips/in
g = 386.4; % Acceleration due to gravity in in/s^2
m = w / g; % Mass of the tower in kips-s^2/in
wn = sqrt(k / m); % Natural frequency of the system in rad/s

% Define the force function P(t)

P_max = 100; % Peak force in kips
T_half = 0.6; % Duration of half sine wave in seconds
force = @(t) (t <= T_half) .* (P_max * sin(pi * t / T_half)); % Piecewise function

% Time setup
t_end = 3.0; % Total time for response calculation
dt1 = 0.1; % First time step interval in seconds
dt2 = 0.05; % Second time step interval in seconds
time1 = 0:dt1:t_end; % Time vector for dt = 0.1 sec
time2 = 0:dt2:t_end; % Time vector for dt = 0.05 sec

%% 1) Displacement response using Trapezoidal Method with dt = 0.1

u1 = zeros(size(time1)); % Initialize displacement vector for dt = 0.1
v1 = zeros(size(time1)); % Initialize velocity vector for dt = 0.1

for i = 1:length(time1)-1
t = time1(i);
t_next = time1(i+1);
% External force at current and next time step
P_current = force(t);
P_next = force(t_next);

% Trapezoidal method for acceleration

a_current = (P_current - k * u1(i)) / m;
u1(i+1) = u1(i) + dt1 * v1(i) + (dt1^2 / (2 * m)) * (P_current - k * u1(i));
v1(i+1) = v1(i) + (dt1 / 2) * (a_current + (P_next - k * u1(i+1)) / m);
end

%% 2) Displacement response using Trapezoidal Method with dt = 0.05

u2 = zeros(size(time2)); % Initialize displacement vector for dt = 0.05
v2 = zeros(size(time2)); % Initialize velocity vector for dt = 0.05

for i = 1:length(time2)-1
t = time2(i);
t_next = time2(i+1);
% External force at current and next time step
P_current = force(t);
P_next = force(t_next);

% Trapezoidal method for acceleration

a_current = (P_current - k * u2(i)) / m;
u2(i+1) = u2(i) + dt2 * v2(i) + (dt2^2 / (2 * m)) * (P_current - k * u2(i));
v2(i+1) = v2(i) + (dt2 / 2) * (a_current + (P_next - k * u2(i+1)) / m);
end
%% 3) Theoretical solution
u_theoretical = zeros(size(time1)); % Placeholder for theoretical
displacement solution
for i = 1:length(time1)
t = time1(i);
if t <= T_half
u_theoretical(i) = (P_max / k) * (1 - cos(wn * t));
else
u_theoretical(i) = (P_max / k) * (1 - cos(wn * T_half)) * cos(wn * (t
- T_half));
end
end

%% 4) Plot the results

figure;
hold on;
plot(time1, u1, 'b', 'DisplayName', 'Numerical Solution (dt = 0.1 s)');
plot(time2, u2, 'r--', 'DisplayName', 'Numerical Solution (dt = 0.05 s)');
plot(time1, u_theoretical, 'g', 'DisplayName', 'Theoretical Solution');
xlabel('Time (s)');
ylabel('Displacement (in)');
legend;
title('Displacement Response of Water Tower');

%% Plot the difference between numerical and theoretical solutions

diff_dt1 = u1 - u_theoretical; % Difference for dt = 0.1
diff_dt2 = interp1(time2, u2, time1) - u_theoretical; % Interpolated
difference for dt = 0.05

figure;
plot(time1, abs(diff_dt1), 'b', 'DisplayName', '|Numerical - Theoretical|
(dt = 0.1 s)');
hold on;
plot(time1, abs(diff_dt2), 'r--', 'DisplayName', '|Numerical - Theoretical|
(dt = 0.05 s)');
xlabel('Time (s)');
ylabel('Absolute Difference (in)');
legend;
title('Difference Between Numerical and Theoretical Solutions');

>>
 Trapezoidal Numerical Integration MATLAB
For Displacement Response Histo
Eora mit in ku
m
t

I tming It é
A sinwalt e de
If
figit f fn sina.coswpt coswot.sinwofd
n

fiigltse coswntdt fiftontfigure sina.tt

f.fr
A B
A F T dT E FLI Tg T COSWDT

B gleldt E
g T
ig t Sin WDT

Trapezoidal Method
A Flo 2 f 20 7 24130 f 24 Cn Dots f n ot
20 t 3413017 29 In Dot not
B 910 2g g
St 0.1 sec
% MATLAB Code for Water Tower Displacement Response
using Trapezoidal Method

% Given Parameters
w = 978.8; % Weight of the tower in kips
k = 100; % Stiffness in kips/in
g = 386.4; % Acceleration due to gravity in in/s^2
m = w / g; % Mass of the tower in kips-s^2/in
wn = sqrt(k / m); % Natural frequency of the system in
rad/s P_max: Sets the maximum force of the
half-sine wave.
% Define the force function P(t)
P_max = 100; % Peak force in kips T_half: Defines the duration (0.6 seconds)
T_half = 0.6; % Duration of half sine wave in of the half-sine wave.
seconds
force = @(t) (t <= T_half) .* (P_max * sin(pi * t / T_half)); % force: Uses an anonymous function to
Piecewise force function define P(t), which applies the half-sine force
until T_half and is zero afterward. This
models an impulse-type load that acts over
a short time.
% Time setup
t_end = 3.0; % Total time for response calculation t_end: Sets the simulation duration to 3
dt1 = 0.1; % First time step interval in seconds seconds.
dt2 = 0.05; % Second time step interval in seconds
time1 = 0:dt1:t_end; % Time vector for dt = 0.1 sec dt1, dt2: Defines the two time intervals,
time2 = 0:dt2:t_end; % Time vector for dt = 0.05 sec Δt=0.1 sec and Δt=0.05 sec, for numerical
analysis.

time1, time2: Creates time vectors for each

interval, which will be used to step through
each time point in the calculations.

%% 1) Displacement response using Trapezoidal Method with

u1, v1: Initializes zero vectors for
dt = 0.1
displacement and velocity,
u1 = zeros(size(time1)); % Initialize displacement vector for
respectively, over time1 (for
dt = 0.1
Δt=0.1 sec). These will store the
v1 = zeros(size(time1)); % Initialize velocity vector for dt =
calculated values for each time step.
0.1

for i = 1:length(time1)-1
% Current time and next time step Loop (for i = 1:length(time1)-1): Iterates
t = time1(i); through each time point in time1 to
t_next = time1(i+1); calculate displacement and velocity.

% External force at current and next time step t, t_next: Sets the current time (t) and the
P_current = force(t); next time (t_next) in time1.
P_next = force(t_next);
P_current, P_next: Calculates the external
force at the current (t) and next (t_next)
time steps using the force function.
 % Trapezoidal method for displacement and velocity a_current: Calculates the current
a_current = (P_current - k * u1(i)) / m; % Current acceleration acceleration using
u_pred = u1(i) + dt1 * v1(i) + (dt1^2 / (2 * m)) * (P_current - k * u1(i));
v_pred = v1(i) + (dt1 / 2) * (a_current + (P_next - k * u_pred) / m);

u_pred: Predicts the displacement at

the next time step (t_next) using the
trapezoidal method:

v_pred: Predicts the velocity at the

next time step using trapezoidal
integration for velocity as well.

% Update displacement and velocity

u1(i+1) = u_pred; u1(i+1), v1(i+1): Stores the predicted
v1(i+1) = v_pred; displacement and velocity values at the
end next time step in the u1 and v1 arrays.
These updated values are used in the
next iteration.

%% 2) Displacement response using Trapezoidal Method

with dt = 0.05 u2, v2: Initializes zero vectors for
u2 = zeros(size(time2)); % Initialize displacement vector displacement and velocity for time steps
for dt = 0.05 defined in time2 (for Δt=0.05 sec).
v2 = zeros(size(time2)); % Initialize velocity vector for dt
= 0.05
for i = 1:length(time2)-1 Loop for dt = 0.05: Repeats the same
% Current time and next time step process as in the previous loop but for
t = time2(i); the finer time step (time2 vector).
t_next = time2(i+1);
t, t_next, P_current, and P_next: Sets
% External force at current and next time step up the current and next time values
P_current = force(t); and forces as before, but with the
P_next = force(t_next); time2 vector.

% Trapezoidal method for displacement and velocity Similar calculations for acceleration,
a_current = (P_current - k * u2(i)) / m; % Current displacement prediction, velocity
acceleration prediction, and updating u2 and v2
u_pred = u2(i) + dt2 * v2(i) + (dt2^2 / (2 * m)) * (P_current - k * arrays for each time step in time2.
u2(i));
v_pred = v2(i) + (dt2 / 2) * (a_current + (P_next - k * u_pred) /
m);

% Update displacement and velocity

u2(i+1) = u_pred;
v2(i+1) = v_pred;
end
time 1- time step @ 0.1 sec time 2- time step @ 0.05 sec
u1- @ time step u2- @ time step
0.1 sec 0.05 sec
v1- @ time step
v2- @ time step
0.1 sec
0.05 sec