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

조회 수: 16(최근 30일)
Kostas 24 Nov 2019
Commented: Kostas 24 Nov 2019
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?
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

로그인 to comment.

### 채택된 답변

David Goodmanson 24 Nov 2019
David Goodmanson 님이 편집함. 24 Nov 2019
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).

#### 댓글 수: 1

Kostas 24 Nov 2019
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!

로그인 to comment.