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
figuresubplot(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 onsubplot(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)% Indexesi = 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
Best Answer