ODE45 calculating the total energy in the system and checking the solver
조회 수: 3 (최근 30일)
이전 댓글 표시
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
댓글 수: 0
채택된 답변
David Goodmanson
2019년 11월 24일
편집: David Goodmanson
2019년 11월 24일
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).
댓글 수: 3
Salah Djerouni
2020년 6월 17일
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 :
djerouni100@gmail.com
thanks in advance .
jasem kamel
2023년 5월 29일
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 :
jasem.m.kamel@gmail.com
thanks a lot in advance .
추가 답변 (0개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!