SDOF with frequency sweep and variable natural frequency

Hi, I'm trying to understand if the syntax of my code is correct and the results I'm seeing are right; the problem I'm studying doesn't get a lot of attention and the literature about it is scarce, so my options are limited and I tried everything but I just cannot find where the error is IF there is.

Basically I'm modelling a simple SDOF system with the equations expressed in the State Space, so instead of a single second order ODE, I have a system of two first order ODEs. The frequency is swept from a starting frequency f0 to an end frequency f1 linearly, instead of a constant frequency (like is done usually). This causes beats in the response after crossing a resonance.

I have implemented this case into my code and checked the results with a paper that considers the LTI system.

Problems arise the moment I try to implement a varying stiffness k(t) linearly with time (the natural frequency of the system gets higher or lower with time based on the case considered, and the stiffness used is linked to it by a simple formula k=(Omega(t)^2)*mass). The system becomes an LTV which I solved using the example in ode45's page with time dependant coefficients.

To know if my code works I tried to confront the results with a paper I found from the 90s, this article uses a dimensionless formulation. I start from this dimensionless fromulation and its parameters and I retrieve the dimensional parameters that I need, then to confront the results I bring the calculated solution back again in the dimensionless formulation. I'm going to attach the code and the article to understand better the problem.

The fact is that, I cannot reproduce the graphs in figure 3, the amplitude is right but not the x values at which the peak shows up, so probably I'm doing something wrong with the solution of the system.

Code:

m=0.1;
fn0=1000;                                                     % Natural frequency 0 rpms
OmegaN_0= 2*pi*fn0; 

f0=0;
f1=5000; 

F0=1;
Zeta_f= 0.025;                                                % Damping
alpha=0.004;                                                  % Adimensional Sweep Rate: OmegaF_dot/OmegaR^2
beta=[-0.002,0.0,0.002];                                      % Adimensional Stiffening Rate: Omegan_dot/OmegaR^2
OmegaR=(alpha.*OmegaN_0)./(alpha-beta);                       % Sweep-Stiffening crossing
OmegaF_dot=(OmegaR.^2).*alpha;                                % Acceleration Rate
OmegaN_dot=(OmegaR.^2).*beta;                                 % Stiffening Rate

c=2*m*Zeta_f.*OmegaR;
q_char=F0./(m.*OmegaR.^2);

t_end=(f1-f0)./(OmegaF_dot/(2*pi));

for i=1:length(beta)

    c_i=c(i);

% Initial conditions
    y0 = [0;0];
    dt=0.000001;
    
    opts = odeset(RelTol=1e-9,AbsTol=1e-9);
    [t_sol, y_sol] = ode89(@(t,y) state(t,y,m,c_i,f0,F0,OmegaF_dot,OmegaN_dot,OmegaN_0,i),0:dt:t_end(i),y0,opts);

% Output
    x = y_sol(:,1);
    v = y_sol(:,2);
    
    f=f0+(OmegaF_dot(i)/(2*pi))*t_sol;
    OmegaF=2*pi*f;
    eta_tau{i}=OmegaF/OmegaR(i); 
    
    Q{i}=x/q_char(i);
    Q_env{i}=abs(hilbert(Q{i}));

end


function dydt=state(t,y,m,c_i,f0,F0,OmegaF_dot,OmegaN_dot,OmegaN_0,i)
 
    F_t=F0*cos((0.5*OmegaF_dot(i)*(t^2))+(2*pi*f0)*t);
% stiffening   
   
    OmegaN=OmegaN_0+OmegaN_dot(i)*((2*pi*f0)/(OmegaF_dot(i)))+OmegaN_dot(i)*t;        % natural omega      [rad/s]
    k=(OmegaN^2)*m;                                                                   % stiffness          [N/m]
   
    dydt=zeros(2,1);
    dydt(1)=y(2);
    dydt(2)=(1/m)*(F_t-c_i*y(2)-k*y(1));
end

I used the same parameters used in the article, also to give context, everything is normalized with respect to OmegaR (crossing between the excitation frequency and natural frequency straight lines) into said paper (page 2 and 5 are sufficient for the scope of the post, everything else is useless for this matlab code).

I just need a check on the syntax for the LTV system solution, the article is present to better understand everything.

Thanks in advance.

[Torsten]

It would help if you included the plot of the figure in the MATLAB code that you want to compare with the figure in the article.

Did you try to reproduce the dimensionless results of the article directly with a suitable MATLAB code ?

[Lucrezia]

Hi Torsten, yes I tried to reproduce the graph converting my result into a dimensionless plot, the code is:

% Plot

subplot(3,1,1)

plot(eta_tau{1}, Q_env{1}, 'Color', colors(1,:), 'LineWidth', 0.5);

grid minor

xlabel('Frequency [Hz]');

ylabel('Displacement x [m]');

title('SDOF sweep through transient');

xlim([0 2])

I repeat this for each value of beta, the dimensionless axes are eta_tau (the excitation frequency divide by OmegaR) and the amplitude (divide by a characteristic amplitude q_char). The amplitude is very close to that of the paper, the values on eta_tau not, for beta higher than zero infact the max peak should be on values smaller of eta_tau than for beta<0.

I get the opposite.

[Torsten]

From the article, I get the impression that alpha, zeta_0 and beta are the only parameters that are used to define and solve the dimensionless model.

Equation (1) of the article reads after the transformation

Q''(tau) + c / (m*omega_r)*Q'(tau) + k(tau/omega_r) / (m*omega_r^2) * Q(tau) = f(tau/omega_r) / (m*q^hat *omega_r^2)

which should be equation (8).

So it should be possible to write this differential equation only in terms of alpha, zeta_0 and beta. Is that correct ?

But in your code, you use values for many other quantities besides alpha, zeta_0 and beta to define the differential equation, e.g. m, q_char, OmegaN_0, F0 etc. So for me it's dubious how you want to compare figure 3 with your results if you don't know what the authors set for the quantities that are not mentionned in the article and that you implicitly assumed.

But I might be wrong with my objection because I lack the physical background.

[Lucrezia]

I know that the physical part of the problem might be confusing, but let's say that I'm "sure" about those values, It is the right approach to solve an LTV system? The core of the problem is to solve dydt where F_t and k change with time. The rest is up to me to figure it out.

Ps: If you are wondering why I use all those term, basically I'm converting their equation to a "standard" formulation by doing the inverse formula, nothing more.

[Torsten]

If I transform equation (8) back to the standard form (1), I get

m * qdotdot + c * qdot + m * (omega0 + omegadot*t)^2 * q = fhat * cos(Theta0 + Omega0*t + 0.5*Omegadot*t^2)

Thus

dydt(1) = y(2)
dydt(2) = (-c * y(2) - m * (omega0 + omegadot*t)^2 * y(1) + fhat * cos(Theta0 + Omega0*t + 0.5*Omegadot*t^2)) / m

And my question was if it really suffices to prescribe alpha, zeta_0 and beta to determine all those constants of the standard equation in order to compare your results with figure (3). E.g. you specify OmegaN_0, m, omega0, F0. You cannot be sure that the authors chose the same values.

By the way: I get omega_R = beta*OmegaN_0./(beta-alpha) instead of OmegaR=(alpha.*OmegaN_0)./(alpha-beta). I assume you set the omega0 in omega(t) = omega0 + omegadot*t equal to 0 ?

[Lucrezia]

F0 is the force amplitude, it doesn't really matters on the response, especially in the dimensionless plot they use, omega0 is the starting natural frequency I put as 1000 Hz (in my code N stands for natural and F for Forcing). OmegaR is the crossing between the two straight lines as showed in the paper and I double checked the formula from another paper, should be right as it is. Tomorrow I'm going to try to solve it with the parameters used by the authors without converting them in dimensional form...I still don't understand if I'm doing it right or not.

[Torsten]

Tomorrow I'm going to try to solve it with the parameters used by the authors without converting them in dimensional form.

In my opinion, you can't reproduce the results in figure (3) with only knowing what the authors set for alpha, beta and zeta0. But maybe you are right and the missing parameters don't matter for the result about the maximum amplitude.

I still don't understand if I'm doing it right or not.

In order to say whether you do it right or wrong I must understand what you are doing.

The authors say that omega_r is where the linear functions

Omega(t) = Omega0 + Omegadot*t
omega(t) = omega0 + omegadot*t

cross.

This happens for

t* = ( Omega0-omega0)/(omegadot-Omegadot)

Thus

omega_r = omega0 + omegadot * t* = omega0 + omegadot* (Omega0-omega0)/(omegadot-Omegadot)

I cannot follow why this equals

OmegaR = (alpha.*OmegaN_0)./(alpha-beta); 

[Lucrezia]

it's the same, try it with: omega=10+0.5*t and Omega=5*t, thus Omega0=0.

[Torsten]

Yes, it works for Omega0 = 0. But who knows whether the authors used Omega0 = 0 for their figure ?

Ok, I looked through your code, and for Omega0 = 0 and f0 = 0, I agree that you correctly implemented the dimensional equation.

But the equation

m * qdotdot + c * qdot + m * (omega0 + omegadot*t)^2 * q = fhat * cos(Theta0 + Omega0*t + 0.5*Omegadot*t^2)

shows that you had to assume values for m, omega0, fhat, Theta0 and Omega0 (=0). At least for me it's not obvious that the results should be qualitatively similar to figure (3) of the article.