How do I create a graph of the phase angle and response of a harmonic excitation in Matlab?

Question

In my course of dynamics of structures, I am struggling with some basic principles. As a practice, I am plotting (or attempting to plot) the curves that are shown in the course book, with the help of Matlab. For example, I want to plot the phase angle theta as a function of beta, the dimensionless frequency, for a damped SDOF. The code I have is the following:

%Stiffness of the spring in [N/m]
k = 1;
%Mass in [kg]
m = 1;
%Damping constant in [Ns/m]
c = 0.1;
%Initial displacement in [m]
u_0 = 0;
%Initial velocity in [m/s]
v_0 = 0;
%% 
%Resonance frequency in [rad/s]
omega_res = sqrt(k/m);
%Resonance frequency in [Hz]
f_res = (omega_res)/(2*pi);
%Eigenperiod in [s]
T_res = 1/f_res;
%% 
%Critical damping in [Ns/m]
c_c = 2*sqrt(k*m);
%Damping ratio in [/]
xi = c/c_c;
%% 
%Damped eigenfrequency in [rad/s]
omega_d = omega_res * sqrt(1-xi^2);
%Damped eigenfrequency for strongly damped system in [rad/s]
omega_hat = omega_res * sqrt(xi^2-1);
%% 
%Timevector
t = 0:0.1:100;
%Betavector
beta_v = 0:0.001:4;
%D as a function of beta
D = 1./sqrt((1-beta_v.^2).^2+(2*beta_v.*xi).^2);

plot(beta_v,D)
axis([0 4 0 10])
xlabel('Dimensionless frequency [-]')
ylabel('Dynamic amplification factor [-]')

%theta as a function of beta
theta = atan((2*beta_v*xi)./(1-beta_v.^2));

figure
plot(beta_v,angle(theta),'m')
axis([0 4 0 pi])
grid on
xlabel('Dimensionless frequency [-]')
ylabel('Phase angle [rad]')

But plotting gives me the curve of an undamped SDOF. The curve I get:

The plot I am supposed to get:

Does anyone have a clue what is wrong in my approach?

A second struggle occurs when you want to plot the the response for a particular beta. The curve I achieved was not what I expected. My code is the following:

%Response
beta_a = 0.2;
D_a = 1./sqrt((1-beta_a.^2).^2+(2*beta_a.*xi).^2);
theta_a = atan((2*beta_a*xi)./(1-beta_a.^2));
p_hat = 1;
omega = beta_a.*omega_res;
u_p = p_hat/k.*D_a.*exp(i.*(omega.*t-theta_a));
B = u_0 - p_hat./k.*D_a.*exp(-i.*theta_a);
A = (v_0 + xi*omega_res*B - p_hat/k*D_a.*i*omega*exp(-i.*theta_a))/omega_d;
u_h = exp(-xi.*omega_res.*t).*(A.*sin(omega_d.*t) + B.*cos(omega_d.*t));
u = u_p + u_h;
figure
plot(t,real(u))
grid on
xlabel('Time [s]')
ylabel('Displacement [m]')

Can anyone help me with this line of thought? What's wrong?

Answer 1

With your code, and with $c = 0.3$, I get the following results:

The amplitude and phase look OK. But the displacement does not show any damping. To see why read on below.

I'm not sure about your notation. So the equations below may differ from those in your textbook.

ODE

The ODE you are trying to solve is $$ \ddot{u} + 2\xi\omega_{\text{res}}\dot{u} + \omega_{\text{res}}^2 u = \frac{\hat{p}}{m}\text{Re}\left[\exp(i\omega t)\right] $$ where $$ \xi = \frac{c}{c_c} ,~~ c_c = 2 m \omega_{\text{res}}, ~~ \omega_{\text{res}} = \sqrt{\frac{k}{m}} $$ and $c$ is the damping factor, $k$ is the stiffness, $m$ is the mass, $\omega$ is the forcing frequency.

Homogeneous solution

The homogeneous solution is $$ u_h = A\,\exp(-\xi\omega_{\text{res}} t)\sin(\omega_d t + \phi_d) $$ where $\phi_d$ is a phase and $$ \omega_d = \omega_{\text{res}}\sqrt{1 - \xi^2} $$

Notice that the exponential decay term occurs only in the homogeneous solution.

Particular solution

The particular solution is $$ u_p = \frac{\hat{p}}{k} \frac{\text{Re}\left[\exp[i(\omega t - \phi)]\right]}{\sqrt{\left[1 - \beta_a)^2\right]^2 + \left[2\xi\beta_a\right]^2}} $$ where $$ \beta_a = \frac{\omega}{ \omega_{\text{res}}}, ~~\tan\phi = \frac{2\xi\beta_a}{1 - \beta_a^2} \,. $$

Your solution

In your solution the homogeneous term is zero and therefore there is no decay term in your solution. Otherwise your implementation looks correct.