How to solve a mass-spring-damper system with external time-dependent forces.

45 views (last 30 days)
I have developed a code for this mass-spring-damper system, but right now it only works without any exernal time-dependent forces. Below you can see the code in which I have tried to add a force to the second mass but it doesn't work. I have tried the interp1 method in this code. It would be extremely helpful for me if someone can tell me why it's not working and what can I do to fix it. Thank you.
%%
% MASS MATRIX
% --------------------------------------------------------
n=2; %Nunber of Masses
pointMass = [1 2 ]; % [kg]
massMatrix = zeros(n,n);
for i = 1 : 1 : length(massMatrix)
massMatrix(i, i) = pointMass(i);
end
% ---------------------------------------------------------
%%
%STIFFNESS MATRIX
% --------------------------------------------------------
k=zeros(1,n+1);
k=[1000 2000 3000 ];
stiffnessMatrix=zeros(n,n);
%First row of the stiffness matrix
stiffnessMatrix(1,1)=k(1)+k(2);
stiffnessMatrix(1,2)=-k(2);
%Last row of the stiffness matrix
stiffnessMatrix(n,n)=k(length(k)-1)+k(length(k));
stiffnessMatrix(n,n-1)=-k(length(k)-1);
for i = 2 : n-1
stiffnessMatrix(i,i)= k(i)+k(i+1);
stiffnessMatrix(i,i-1)=-k(i);
stiffnessMatrix(i,i+1)=-k(i+1);
end
% --------------------------------------------------------
%%
%DAMPING MATRIX
% --------------------------------------------------------
d=zeros(1,n+1);
d=[1 2 3 ];
dampingMatrix=zeros(n,n);
%First row of the stiffness matrix
dampingMatrix(1,1)=d(1)+d(2);
dampingMatrix(1,2)=-d(2);
%Last row of the stiffness matrix
dampingMatrix(n,n)=d(length(d)-1)+d(length(d));
dampingMatrix(n,n-1)=-d(length(d)-1);
for i = 2 : n-1
dampingMatrix(i,i)= d(i)+d(i+1);
dampingMatrix(i,i-1)=-d(i);
dampingMatrix(i,i+1)=-d(i+1);
end
% --------------------------------------------------------
%coefficient matrices
O=zeros(length(massMatrix)); I=eye(length(massMatrix));
A=[O I; -inv(massMatrix)*stiffnessMatrix -inv(massMatrix)*dampingMatrix];
B=[O; inv(massMatrix)];
%Force vector
F1=0; F2=10;
w=5; %frequency
ft=linspace(1,4,25);
f1=F1*sin(2*pi*w*ft);
f2=F2*sin(2*pi*w*ft);
f=[f1 f2 ]';
%time step and time vector
dt=0.01;
time = 0:dt:4;
y0= [0.01 0 0 0]; %initial condtiions (x1(0)... xn(0) x1dot(0)... xndot(0))
[tsol,ysol] = ode45(@(t,y) testode_2D(t,y,ft,f,A,B), time, y0);
for p = 1 : n
subplot(n,1,p),
txt = ['x',num2str(p)];
plot(time, ysol(:,p),'DisplayName',txt),
if p == 1
title('Displacement as a function of time','FontSize',30)
end
ylabel('displacement (m)','FontSize',20)
if p == n
xlabel('time (s)' ,'FontSize',20);
end
grid on
legend show
end
%Fucntion-----------------------------------------
function dy=testodeF_2D(t,y,ft,f,A,B)
f = interp1(ft,f,t);
dy= A*y+B*f;
end
%-----------------------------------------

Accepted Answer

William Rose
William Rose on 13 Aug 2021
I am not familiar with the method you used in your code to derive the matrices A and B that describe the system. Therefore I derived formulas for A and B using a different approach. My approach is probably less elegant than your way, but I understand my way. It results in the same values for A and B as what your program gives - which is good! See attached document for my approach.
We had a problem like this in my first semster of college, in 1977. There is a way to diagonalize the system matrix and find eigenvectors. The solution is a linear combination of eigenvectors. If the system is symmetric (in other words, the masses are the same, and the springs and dampers on the sides have the same respective values), then the two eigenvectors are 1. a solution in which the two masses move opposite to one another, like a mirror image, and 2. a solution in which the two masses move in parallel, both left, then both right, etc. In the second solution, the distance between the masses never changes.
  2 Comments
William Rose
William Rose on 13 Aug 2021
There was an error in the PDF document I posted earlier. The matrix B was transposed from what it should have been. The correct version is attached here.
William Rose
William Rose on 13 Aug 2021
Here is a link to an eigenvector analysis of the problem. I googled "coupled masses eigenvector solution". The search returned these lecture notes from Harvard, which is where I was a freshman in 1977, when I first met this problem. It seems they are still teaching the same problem, 44 years later.

Sign in to comment.

More Answers (2)

William Rose
William Rose on 13 Aug 2021
I downloaded your code and saved it as coupledMasses1DextForce0.m. When I ran it, I got an error message about the function definition:
>> coupledMasses1DextForce0
Undefined function or variable 'testode_2D'. ...
This error occurs because the function name in the declaration at the end of the script is "testodeF_2D", which has "F" in it. I removed the "F" from the function name at the end of the script, and I ran the script again. I got a new error message:
>> coupledMasses1DextForce0
Error using interp1>reshapeAndSortXandV (line 424)
X and V must be of the same length. ...
Error in coupledMasses1DextForce0>testode_2D (line 86) ...
This happens because you have defined f so that it is 50x1. It should be 2x25 (given your choice for ft, which, as I explain below, is not a good choice - but we'll get to that later). Thefore I change
f=[f1 f2 ]';
to
f=[f1;f2]';
and re-run. Now I get a new error:
>> coupledMasses1DextForce0
Error using *
Incorrect dimensions for matrix multiplication. ...
Error in coupledMasses1DextForce0>testode_2D (line 87)
dy= A*y+B*f;
This happens because interp1() returns a row vector for interpolated f, but you need a column vector. Therefore we change
dy= A*y+B*f;
to
dy= A*y+B*f';
and re-run. Now the program runs to completion, with no errors. This is good!
However, the figure that is generated shows not data, and the time range is not 0 to 4 seconds as we expect it to be. Here is a screenshot of the figure:
What is going on?
For one thing, the choice of the vector of times, ft, is not good. You have:
ft=linspace(1,4,25);
However, you also specify a time span of 0 to 4 for integration. Therefore your range for ft should be 0 to 4 (or wider). Also, you specify 25 points between t=1 and t=4, i.e. about 8 Hz sampling. That is not enough, because your sine wave frequency is 5 Hz. You shoudl sample at more than twice the sine wave frequency. I suggest
ft=linspace(0,4,201);
which gives 50 samples per second from t=0 to t=4. This is 10 samples per cycle, which is enough for a fairly smooth re-creation of the sinusoid by interpolation.
I am also deleting the 'Fontsize' instructions in the plotting section, since the labels are not displaying well. Then re-run. The plot is not much different. The reason is that the array of results, ysol(), is all NaNs. Why?
To be continued.

William Rose
William Rose on 13 Aug 2021
You do not need to use interpolation to add external forcing. I think a better approach is to compute the forces directly in the function that returns the derivatives. The new function is
function dy=testode_2D(t,y,A,B)
F1=0; %amplitude 1
F2=10; %amplitude 2
w=5; %frequency (cycles/s)
f=[F1;F2]*sin(2*pi*w*t);
dy= A*y+B*f;
end
and the call to ode45() is now
[tsol,ysol] = ode45(@(t,y) testode_2D(t,y,A,B), time, y0);
Also, you define the time span as a vector of evenly spaced times:
time = 0:dt:4;
I recommend defining the beginning and end times only. Let ode45() compute the time steps. It is designed to do that efficiently, in order to mantain a high level of accuracy. Therefore define time as follows:
time = [0,4]; %time span
With these changes, the script gives nice results. This is good! The script produces the figure below.
The script is attached.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!