Clear Filters
Clear Filters

Generating codim-1 bifurcation diagram of sliding friction oscillator

57 views (last 30 days)
I am performing continuation analysis in MATCONT and I am finding it difficult to generate co-dim bifurcation diagram.
Here is my linearized system of equations which represents a 3-dof clutch dynamics with sliding type friction:
x1' = x2
x2' = -(de+Fn*mukin)*x2/je + (Fn*mukin*x4)/je - ((x2-x4)*mukin + mus)/je
x3' = x4
x4' = (Fn*mukin*x2)/jd - (ds+Fn*mukin)*x4/jd +ds*x6/jd + ks*x7/jd
x5' = x6
x6' = ds*x4/jv -ds*x6/jv - ks*x7/jv
x7' = x6 - x4
This system dynamics represents clutch dynamics with the sliding friction between inertias je and jd. jv is the load inertia. ks is stiffness between jd and jv. de,ds is damping coefficient. Fn is normal force acting on jd and je during engagement of clutch. mukin is slope of friction curve usually between [-0.01 to 0.01] and mus is static coefficient of friction 0.136 for this analysis.
je = 1.57
jd = 0.0066
jv = 4
ks = 56665.5
ds,de = 0.1
Fn = 8000
Please guide me how to perform the continuation of this system in MATCONT.

Answers (1)

Umar
Umar on 1 Sep 2024

Hi @Junaid Ali ,

You asked, “Please guide me how to perform the continuation of this system in MATCONT. “

Please see my response to your comments below.

Here’s how you can set up and perform the continuation analysis using MATCONT. So, I can assume that you have MATCONT already installed and added to your MATLAB path. If not, please click the following links for more information.

https://www.mathworks.com/matlabcentral/answers/1883987-how-to-load-matcont-data-file-and-plot

https://www.mathworks.com/matlabcentral/answers/2148254-matcont-load-state-space-a-b-c-d-matrices-for-continuation-analysis

https://www.mathworks.com/matlabcentral/answers/99328-how-do-i-create-a-bifurcation-diagram-in-matlab

https://www.mathworks.com/matlabcentral/answers/770418-how-to-plot-bifurcation-diagram-in-matlab

https://dercole.faculty.polimi.it/tds/matcont.html

You have provided the following system of equations as shown in your comments.

function dx = clutchDynamics(t, x, params)
  % Unpack parameters
  de = params.de;
  ds = params.ds;
  Fn = params.Fn;
  mukin = params.mukin;
  mus = params.mus;
  je = params.je;
  jd = params.jd;
  jv = params.jv;
  ks = params.ks;
    % State variables
    x1 = x(1);
    x2 = x(2);
    x3 = x(3);
    x4 = x(4);
    x5 = x(5);
    x6 = x(6);
    x7 = x(7);
    % Define the system of equations
    dx = zeros(7, 1);
    dx(1) = x2;
    dx(2) = -(de + Fn * mukin) * x2 / je + (Fn * mukin * x4) / je - ((x2 - x4) * 
    mukin + mus) / je;
    dx(3) = x4;
    dx(4) = (Fn * mukin * x2) / jd - (ds + Fn * mukin) * x4 / jd + ds * x6 / jd +
    ks * x7 / jd;
    dx(5) = x6;
    dx(6) = ds * x4 / jv - ds * x6 / jv - ks * x7 / jv;
    dx(7) = x6 - x4;
  end

You also defined the parameters to be used in your system. I would suggest to implement them in a separate MATLAB script or function:

params.de = 0.1;
params.ds = 0.1;
params.Fn = 8000;
params.mukin = linspace(-0.01, 0.01, 100); % Varying mukin for continuation
params.mus = 0.136;
params.je = 1.57;
params.jd = 0.0066;
params.jv = 4;
params.ks = 56665.5;

To use MATCONT, you need to set up the initial conditions and the continuation parameters. Here’s how you can do it:You need to find an equilibrium point of your system which can be done by solving the system of equations for a specific value of mukin and then initialize MATCONT with the equilibrium point and specify the parameters for continuation. If you are not familiar with MATCONT commands below, I will suggest to review them by clicking the link below.

https://webspace.science.uu.nl/~kouzn101/NBA/ManualMatcontMNov2017.pdf

% Initial conditions (example)
x0 = [0; 0; 0; 0; 0; 0; 0]; % Initial guess for the state variables
mukin0 = 0; % Initial value of mukin
% Define the continuation parameters
options = contset; % Create a default options structure
options = contset(options, 'MaxNumPoints', 100); % Maximum number of 
points
options = contset(options, 'Singularities', 1); % Track singularities
% Start continuation
[x, v, s] = cont(@clutchDynamics, x0, mukin0, options, params);

Once you have performed the continuation, now you can plot the bifurcation diagram which can be done using MATLAB,

figure;
plot(x(1, :), x(2, :)); % Plotting x1 vs x2 as an example
xlabel('x1');
ylabel('x2');
title('Bifurcation Diagram of 3-DOF Clutch Dynamics');
grid on;

Now, to analyze co-dimension bifurcations, you may need to vary more than one parameter. My suggestion would be setting up a second parameter in a similar manner to mukin and use MATCONT to track bifurcations in a multi-parameter space.

  • Additional Suggestions*

In my opinion, it may be beneficial to perform a sensitivity analysis on the parameters such as mukin, mus to understand their influence on system behavior. Also, depending on your results, consider investigating different types of bifurcations such as saddle-node or Hopf bifurcations, which can provide deeper insights into system stability and a linear stability analysis around the bifurcation points to deduce the nature of the bifurcations.

Hope this should help you get started with your project now.

  25 Comments
Junaid Ali
Junaid Ali on 25 Sep 2024 at 14:04
Edited: Junaid Ali on 25 Sep 2024 at 14:51
Thank you @Umar
I get it what you did. I think my next step shoud be defining a jacobian function and calculate eigen values at each equilibrium point using that jacobian. And see if I can identify any Hopf point.
% Define system parameters
params.je = 1.57;
params.jd = 0.0066;
params.jv = 4;
params.ks = 56665.5;
params.de = 0.01;
params.ds = 0.01;
params.Te = 100;
params.Tv = 80;
params.mu0 = 0.11;
params.muk = 0.056; % Initial value for muk
params.np = 14;
params.Rm = 0.0506;
params.Fn = 15000;
params.Kg = 1000;
params.R = 50; % Initial value for R
% Initial conditions
x0 = [0.01; 0.01; 0.01; 0.01; 0.01]; % Small perturbation
ts = 0:0.001:1;
% Solve ODE
[t,z] = ode45(@(t,z)clutch3dof(t,z,params), ts, x0);
figure;
plot(ts, z);
title('System Response Over Time');
xlabel('Time (s)');
ylabel('State Variables');
grid on;
% Finding equilibrium
options = optimoptions('fsolve', 'Display', 'off', 'TolFun', 1e-6, 'TolX', 1e-6, 'MaxFunctionEvaluations', 20000, 'MaxIterations', 20000);
[x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x, params), x0, options);
if exitflag > 0
disp('Equilibrium found:');
disp(x_eq);
else
warning('Failed to find equilibrium. Check initial conditions or parameter values.');
end
Equilibrium found:
1.0e+03 * 2.0000 1.9500 1.9500 -0.0000 -0.0150
% Bifurcation analysis
params.muk_range = linspace(-0.1, 0.1, 10); % Range for muk
params.R_range = linspace(0, 150, 10); % Range for R
bifurcation_points_muk = []; % Initialize bifurcation points for muk
bifurcation_points_R = []; % Initialize bifurcation points for R
% Loop through muk range
for muk = params.muk_range
params.muk = muk;
[x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x, params), x0, options);
if exitflag > 0
bifurcation_points_muk = [bifurcation_points_muk; muk, x_eq(2)];
end
end
% Loop through R range
for R = params.R_range
params.R = R;
[x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x, params), x0, options);
if exitflag > 0
bifurcation_points_R = [bifurcation_points_R; R, x_eq(2)];
end
end
% Plot bifurcation points for muk
figure;
plot(bifurcation_points_muk(:,1), bifurcation_points_muk(:,2), 'o-');
title('Bifurcation Diagram: y2 vs muk');
xlabel('muk');
ylabel('y2 (Equilibrium State)');
grid on;
% Plot bifurcation points for R
figure;
plot(bifurcation_points_R(:,1), bifurcation_points_R(:,2), 'o-');
title('Bifurcation Diagram: y2 vs R');
xlabel('R');
ylabel('y2 (Equilibrium State)');
grid on;
figure
% Initialize previous eigenvalue storage
E_vals_prev = [];
% Loop through values of muk
for muk = params.muk_range
params.muk = muk;
% Calculate eigenvalues for the Jacobian matrix
E_vals = eig(Jacobian2DOF(x_eq, params));
% Normalize the eigenvalues to lie within the unit circle
E_vals_normalized = E_vals ./ abs(E_vals);
% Plot the current normalized eigenvalues
hold on
plot(real(E_vals_normalized), imag(E_vals_normalized), '*b', 'MarkerSize', 8);
% If this is not the first iteration, plot arrows to show movement
if ~isempty(E_vals_prev)
for i = 1:length(E_vals)
% Plot arrows indicating movement from previous to current eigenvalue
quiver(real(E_vals_prev(i)), imag(E_vals_prev(i)), ...
real(E_vals_normalized(i)) - real(E_vals_prev(i)), ...
imag(E_vals_normalized(i)) - imag(E_vals_prev(i)), ...
0, 'MaxHeadSize', 0.5, 'Color', 'k'); % Arrows in red
end
end
% Store the current eigenvalues for the next iteration
E_vals_prev = E_vals_normalized;
% Plot the unit circle for reference
theta = linspace(0, 2*pi, 100);
plot(cos(theta), sin(theta), '--k'); % Unit circle
% Set up axis labels, grid, and aspect ratio
xlabel('Real Part');
ylabel('Imaginary Part');
axis equal;
grid on;
end
% Equilibrium function
function eq = equilibrium_function(x, params)
dxdt = clutch3dof(0, x, params);
eq = dxdt;
end
function jacob = Jacobian2DOF(y_eq,params)
je = params.je;
jd = params.jd;
jv = params.jv;
ks = params.ks;
de = params.de;
ds = params.ds;
mu0 = params.mu0;
muk = params.muk;
np = params.np;
Rm = params.Rm;
Fn = params.Fn;
Kg = params.Kg;
y1 = y_eq(1);
y2 = y_eq(2);
y3 = y_eq(3);
y4 = y_eq(4);
y5 = y_eq(5);
% First row of Jacobian (A1)
A1_1 = (-muk * np * Rm * (Kg * y5 + Fn) * tanh(2 * y1 - 2 * y2) + ...
(muk * (y2 - y1) - mu0) * np * Rm * (Kg * y5 + Fn) * (2 - 2 * tanh(2 * y1 - 2 * y2)^2) - de) / je;
A1_2 = (muk * np * Rm * (Kg * y5 + Fn) * tanh(2 * y1 - 2 * y2) + ...
(muk * (y2 - y1) - mu0) * np * Rm * (Kg * y5 + Fn) * (-2 + 2 * tanh(2 * y1 - 2 * y2)^2)) / je;
A1_3 = 0;
A1_4 = 0;
A1_5 = ((muk * (y2 - y1) - mu0) * np * Rm * Kg * tanh(2 * y1 - 2 * y2)) / je;
% Second row of Jacobian (A2)
A2_1 = (muk * np * Rm * (Kg * y5 + Fn) * tanh(2 * y1 - 2 * y2) - ...
(muk * (y2 - y1) - mu0) * np * Rm * (Kg * y5 + Fn) * (2 - 2 * tanh(2 * y1 - 2 * y2)^2)) / jd;
A2_2 = (-muk * np * Rm * (Kg * y5 + Fn) * tanh(2 * y1 - 2 * y2) - ...
(muk * (y2 - y1) - mu0) * np * Rm * (Kg * y5 + Fn) * (-2 + 2 * tanh(2 * y1 - 2 * y2)^2) - ds) / jd;
A2_3 = ds / jd;
A2_4 = ks / jd;
A2_5 = -(muk * (y2 - y1) - mu0) * np * Rm * Kg * tanh(2 * y1 - 2 * y2) / jd;
% Third row of Jacobian (A3)
A3_1 = 0;
A3_2 = ds / jv;
A3_3 = -ds / jv;
A3_4 = -ks / jv;
A3_5 = 0;
% Fourth row of Jacobian (A4)
A4_1 = 0;
A4_2 = -1;
A4_3 = 1;
A4_4 = 0;
A4_5 = 0;
% Fifth row of Jacobian (A5)
A5_1 = -1;
A5_2 = 1;
A5_3 = 0;
A5_4 = 0;
A5_5 = 0;
jacob = [
A1_1, A1_2, A1_3, A1_4, A1_5;
A2_1, A2_2, A2_3, A2_4, A2_5;
A3_1, A3_2, A3_3, A3_4, A3_5;
A4_1, A4_2, A4_3, A4_4, A4_5;
A5_1, A5_2, A5_3, A5_4, A5_5];
end
function dydt = clutch3dof(t,y,params)
je = params.je;
jd = params.jd;
jv = params.jv;
ks = params.ks;
de = params.de;
ds = params.ds;
Te = params.Te;
Tv = params.Tv;
mu0 = params.mu0;
muk = params.muk;
np = params.np;
Rm = params.Rm;
Fn = params.Fn;
R = params.R;
Kg = params.Kg;
dydt = zeros(5,1);
dydt(1) = (Te - de*y(1) - np*Rm*(Fn+Kg*y(5))*(mu0 + muk*(y(1) - y(2)))*tanh(2*(y(1) - y(2))))/je;
dydt(2) = (ks*y(4) + ds*(y(3) - y(2)) + np*Rm*(Fn+Kg*y(5))*(mu0 + muk*(y(1) - y(2)))*tanh(2*(y(1) - y(2))))/jd;
dydt(3) = (-ks*y(4) - ds*(y(3) - y(2))- Tv)/jv ;
dydt(4) = y(3) - y(2);
dydt(5) = R - (-y(2) + y(1));
end
Umar
Umar on 26 Sep 2024 at 1:49
Hi @Junaid Ali,
Your methodology appears sound, and with careful execution of these steps, you should be able to derive meaningful insights about the dynamics of your system. Pay close attention to the stability conditions indicated by the eigenvalues, especially when identifying potential Hopf points, as they can indicate transitions between different dynamic behaviors in your system. Good luck with your analysis!

Sign in to comment.

Categories

Find more on General Applications in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!