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?