ODE45 calculating the total energy in the system and checking the solver

Hi all,

I have a simple code written for a 3 degree of freedom damped cart-mass system. My aim is that besides by solving for the instantaneous displacement and velocity, to also check the solver by determining the energy in the system, which should be constant at any point in time t.

In the first case where the damping constant 'c' is set to 0, I get somewhat expected results, as the sum of the potential (PE) and kinetic (KE) energies is pretty much constant along the entire solution time (I assume it is reasonable for the sum of energies not to be exactly constant in this case since the solver isn't perferct, but feel free to enlighten me if that's not the case). However, in the second case, where damping is included, I find that the dissipated energy due to damping is much greater than the sum of the potenial and kinetic energy at any point in time, which does not seem reasonable (I am expecting that: KE+PE at time t + total dissipated energy from time 0 to time t = KE+PE at time 0).

So this brings me to my question: Am I doing something wrong with the solver such that it doesn't take in to account the damping correctly, or am I doing something wrong in the way that I calculate the energy dissipated from the damping?

Thanks for your help in advance,

KMT.

clear ; clc

% User Inputs

    % Coefficients
    m = 3 ;         % Mass, kg
    k = 20 ;        % Spring constant, N/m
    c = 2 ;         % Damper constant, Ns/m

    % Degrees of freedom
    N = 3 ;

    % Solver inputs
    tmax = 30 ;     % Maximum solution time
    vi = 2 ;        % Intiial velocity of 1st DoF

% Matrices
M = diag(repmat(m, 1, N)) ;
K = tridiag(repmat(k, 1, N+1)) ;
C = tridiag(repmat(c, 1, N+1)) ;

% ODE Solution
tspan = [0 tmax] ;
y0 = [zeros(1, N) vi zeros(1, N-1)] ;

[t, y] = ode45(@(t, y) odefcn(t, y, M, K, C, N), tspan, y0) ;

[Ec, En] = energy(y, M, K, C, N) ;

% Plotting
figure

subplot(2, 1, 1)
    yyaxis right
    plot(t, y(:,1:N))
    ylabel('Displacement') ; xlabel('Time')
    grid on

    yyaxis left
    plot(t, y(:,N+1:2*N))
    ylabel('Velocity') ; xlabel('Time') 
    legend(sprintfc('Mass %d', 1:N))
    grid on

subplot(2, 1, 2)
    yyaxis right
    plot(t, Ec)
    ylabel('Kinetic + Potential Energy') ; xlabel('Time')

    yyaxis left
    plot(t, En) 
    ylabel('Dissipated Energy') ; xlabel('Time') 
    grid on

% Energy function
function [Ec, En] = energy(ymat, M, K, C, N) 

% Indexes
i = N ; 
j = N+1 ; 
k = 2*N ; 

% Extract Variables
y  = ymat(:, 1:i)' ;
dy = ymat(:, j:k)' ;

% Energy Equations
KE = 0.5 * dy' * M * dy ;   % Kinetic energy
PE = 0.5 * y' * K * y ;     % Potential energy
FE = abs(y' * C * dy) ;     % Dissipated energy

Ec = diag(KE + PE) ;
En = cumsum(diag(FE)) ; 

end

% ODE function
function dydt = odefcn(~, y, M, K, C, N)

% Indexes
i = N ;
j = N+1 ;
k = 2*N ;

% Preallocate
dydt = zeros(k, 1) ;

% Equation
dydt(1:i) = y(j:k) ;
dydt(j:k) = M \ (-C*y(j:k) -  K*y(1:i)) ;

end

% Tri-diagonal matrix function
function X = tridiag(x_c)

 x(:,:,1) = diag(x_c(1:end-1)) ;
 x(:,:,2) = diag(x_c(2:end)) ;
 x(:,:,3) = diag(-x_c(2:end-1), -1);
 x(:,:,4) = diag(-x_c(2:end-1), 1) ;
 
 X = sum(x,3) ;
 
end

Answer 1

[David Goodmanson]

Hi Kostas,

The problem is in the calculation of the dissipated energy which is not any kind of sum of (y C dy/dt), but rather the sum of

dEn = dy C dy/dt = C (dy/dt)^2 dt.     

I just brought t into the energy function (naturally t gets included in the call to that function) and used the simpleminded cumtrapz function for the integration. The energy function becomes

function [Ec, En] = energy(t, ymat, M, K, C, N)
%                          ^
% % Indexes
i = N ; 
j = N+1 ; 
k = 2*N ; 
% Extract Variables
y  = ymat(:, 1:i)' ;
dy = ymat(:, j:k)' ;
% Energy Equations
KE = 0.5 * dy' * M * dy ;   % Kinetic energy
PE = 0.5 * y' * K * y ;     % Potential energy
FE = abs(dy' * C * dy) ;    % Dissipated energy  <--
Ec = diag(KE + PE) ;
En = cumtrapz(t,diag(FE)) ; %                    <--
end

after which, things work.

p.s. As a matter of personal preference I found the plot of Ec and En to be, not exactly misleading, but with two different vertical scales it's hard to tell by eye what is going on. For a diagnostic I just went to a simple plot of plot(t,Ec,t,En).

[KostasK]

I see, so you found the power and integrated it over time, so my implementation of the theory was wrong.

Works perfectly, thanks for that!

[Salah Djerouni]

Hello KostasK , i'm interested with the code of energy dissipated can you send me the code matlab if you possible sure .

this is my email :

[e-mail]

thanks in advance .

[jasem kamel]

Dear https://www.mathworks.com/matlabcentral/profile/authors/15803667-kostask

i'm interested with the code of energies measurement ( dissipated, Potential and K.E) can you please send me the code matlab if you possible sure .

this is my email :

[e-mail]

thanks a lot in advance .