3
$\begingroup$

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:

enter image description here

The plot I am supposed to get:

enter image description here

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?

$\endgroup$
  • 1
    $\begingroup$ It would help if you include pictures of what you expect the responses to be and what you're actually getting from your code. $\endgroup$ – Paul Feb 21 '16 at 14:38
  • $\begingroup$ I am not allowed to post more than two pictures, so I included the two of the phase angle. $\endgroup$ – Charcuterie Charizard Feb 21 '16 at 21:56
  • 1
    $\begingroup$ 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. $\endgroup$ – Biswajit Banerjee Feb 22 '16 at 22:24
  • $\begingroup$ I added the missing code for plotting the phase angle. $\endgroup$ – Charcuterie Charizard Feb 23 '16 at 20:19
4
$\begingroup$

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

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ I have some questions about your explanation: - In the critical damping coeffcient, you have an extra term m for the mass which is unknown for me. - Is there a reason you are increasing the damping to 0.3? - Due to the initial condition of u_0 = 0, the curve has to start at a displacement that equals 0? $\endgroup$ – Charcuterie Charizard Feb 23 '16 at 20:29
  • $\begingroup$ And where is the term phi_d coming from? $\endgroup$ – Charcuterie Charizard Feb 23 '16 at 20:43
  • 1
    $\begingroup$ 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. $\endgroup$ – Biswajit Banerjee Feb 23 '16 at 20:52
  • $\begingroup$ It is obvious to search A and B, using the total solution instead of only the homogeneous one. I changed it in my code. Unfortunately, I am still not getting the results I expect. $\endgroup$ – Charcuterie Charizard Feb 24 '16 at 18:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.