40 views (last 30 days)

Hi all,

I have written a code of a damped mass-spring system, that takes as inputs the relevant mass, damping, stiffness and excitation force matrices/vectors, and outputs the displacement and the velocity of each mass.

What I would like to do is to be able to extract the acceleration as well (i.e. the value of dxdt(m:n) as written in the odefcn function at the bottom of the code) in addition to the displacement and velocity outputs.

Now I know that I can find the gradient of the velocity, or reconstruct my equation of motion using my solutions for the displacement and velocity and solve for the acceleration. However, having tried, I would like to avoid these approaches, as the actual problem that I am being faced with (which has significantly different equations of motion) has noisy velocity vectors, and the equation is quite complicated to reconstruct and solve for the acceleration (doable but I'd rather avoid it).

So her is the code:

clear

clc

% Inputs

% Degrees of Freedom

N = 3 ;

% Masses (kg)

m = repmat(100, 1, N) ;

% Spring coefficients (N/m)

k = repmat(15, 1, N+1) ;

% Damping coefficients (Ns/m)

c = repmat(5, 1, N+1) ;

% Maximum force (N)

f0 = repmat(100, 1, N) ;

% Phase angle between force (rad)

phi = linspace(0, pi, N) ;

% Frequency of force (Hz)

freq = 0.2 ;

% Time of simulation (s)

tmax = 100 ;

% Mass initial positions (m)

xi = zeros(1, N) ;

% ODE Inputs

% Degrees of Freedom

N = length(m) ;

% Mass Matrix

M = diag(m) ;

% Stiffness Matrix

k1 = diag(k(1:end-1) + k(2:end)) ;

k2 = diag( -k(2:end-1), 1) ;

k3 = diag( -k(2:end-1), -1) ;

K = k1 + k2 + k3 ;

% Damping Matrix

c1 = diag(c(1:end-1) + c(2:end)) ;

c2 = diag( -c(2:end-1), 1) ;

c3 = diag( -c(2:end-1), -1) ;

C = c1 + c2 + c3 ;

% ODE Initial Conditions

% Time

tspan = [0 tmax] ;

% Displacement and velocity

x0 = [ xi zeros(1, N) ] ;

% Solution

[t, xSol] = ode45(@(t, x) odefcn(t, x, M, K, C, f0, freq, phi, N), tspan, x0) ;

% Results

x = xSol(:, 1:N) ; % Displacement (m)

xdot = xSol(:, N+1:end) ; % Velocity (m/s)

% ODE Function

function dxdt = odefcn(t, x, M, K, C, f0, freq, phi, N)

% Indices

m = N + 1 ;

n = N * 2 ;

% Preallocate

dxdt = zeros(n, 1) ;

% ODE Equation

dxdt(1:N) = x(m:n) ;

dxdt(m:n) = - M \ K * x(1:N) - M \ C * x(m:n) + M \ f0' .* sin(2*pi*freq * t + phi') ;

end

Thanks for your help in advance,

KMT.

James Tursa
on 13 Sep 2019

Edited: James Tursa
on 13 Sep 2019

John D'Errico
on 14 Sep 2019

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.