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

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;
%%
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;
%%
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 [-]')


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?

• It would help if you include pictures of what you expect the responses to be and what you're actually getting from your code. – Paul Feb 21 '16 at 14:38
• I am not allowed to post more than two pictures, so I included the two of the phase angle. – Charcuterie Charizard Feb 21 '16 at 21:56
• Your code for plotting the phase angle seems to be wrong. Could you add that to your question? Otherwise the rest of the code is OK. See my answer below. – Biswajit Banerjee Feb 22 '16 at 22:24
• I added the missing code for plotting the phase angle. – Charcuterie Charizard Feb 23 '16 at 20:19

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 critical damping coeff $c_c = 2\sqrt(k m) = 2 m\sqrt{k/m} = 2 m \omega_{\text{res}}$. The damping coeff is higher just to make the damping more obvious in the figures. When using the initial condition to find $A$ and $B$ make sure that you use $u = u_h + u_p$. The $\phi_d$ term comes if you write your equation in the compact form in the updated answer. – Biswajit Banerjee Feb 23 '16 at 20:52