Solve Displacement Driven Mass/Spring System
显示 更早的评论
I try to find a way to calculate the transient behavior of a mass pre-loaded by a spring but driven by a harmonic mass displacement. Single mass oscillators are described with either driving forces at the mass center or displacements at the root point. But my physical problem is that a mass is harmonically driven by a surface touching the mass. The turning point (x = 0) is NOT the equilibrium state. The spring is pre-loaded.
I want to calculate the driving displacement frequency when the contact force between mass and surface is < 0 (mass starts loosing the contact when the surface is moving down again).
All I can manage so far is to solve the standard equation for harmonic force driven mass point oscillation. The equation is not considering a spring pre-loaded condition (I have no idea to manage that either).
Any guidance to go forward?
Thanks.
standard equation of motion to a force driven mass oscillator without pre-loaded condition (this does not fit to my problem!):
equ = @(t,x) [x(2); (-d/m*x(2)-k/m*x(1)+A/m*sin(wd*t))];
% xi are start conditions [x0, v0]
xi = [x0; v0];
% options structure for ODE solvers; % solver statistics
opts = odeset('stats','on','RelTol',1e-8);
delt = .0001;
tspan = 0:delt:0.03;
% solving differential equation using ode78 solver w/ high quality
tic, [t,x_]=ode78(equ, tspan, xi, opts); toc
x = x_(:,1); % x-displacement vector
xdot = x_(:,2); % x-dot or velocity vector
my physical problem to solve:
回答（1 个）
William Rose
2023-11-7
2 个投票
You have written a system of first order differential equations, but I don't think it describes the system you have drawn and described. I say that becaus ethe system that you have described suggests that the mass is sometimes in contact with the driving plate, and sometimes not. If that is true, then your differential equation should have an "if" statement, or equivalent. (More on that below.)
You say "I want to calculate the driving displacement frequency...". Isn't the frequency of the driving plate a constant? If so, then it will constrain the mass to oscillate at the same frequency. If the frequency of the driving plate is not constant, then what determines the driving plate frequency?
Does the driving plate have a sinusoidal displacement pattern which is independent of the load on it? If yes, then the position of the mass equals the position of the plate when they are in contact. This will be challenging to model, because it implies infinite acceleration when the mass "collides" with the plate and suddenly changes velocity from its "free" velocity to the plate velocity.
Your equation for dx(2)/dt (also known as dv/dt) is
Term 1 on the right side looks like a damping term, where d is the damping coefficient. Is that correct? You did not mention damping in your description. Term 2 on the right side is the force from a spring which produces zero force when x=0. You say "The turning point (x = 0) is NOT the equilibrium state. The spring is pre-loaded." Please explain. Perhaps you are referring to gravity acting on the mass, to pull the mass down from the "zero" position of the spring. But there is no gravity term in your equation for dv/dt. Term 3 is an external forcing term. The units only work out right if A hase units of force. Is that correct? If so, then it means the mass is subjected to an external force, not an externally imposed displacement. Please clarify.
Differential equations with "if" statements are tricky. I have dealt with this problem when modeling the cardiovascular system, which includes one-way valves. (See here for a recent example.) In my experience, the most effective way to deal with the "if" statement is to convert it to a steep but smooth (i.e. differentiable) function that goes from 0 to 1.
7 个评论
Joerg Fricke-Schmidt
2023-11-8
William Rose
2023-11-8
As I mentioned above, there is the possibility of infinite acceleration, when the mass collides with the plate. Therefore I would use first order forward Euler method, with some if statements, rather than Matlab's ode45() or ode15s() or similar routine. I tried the first order Euler approach, with if statements, in an Excel spreadsheet. I used mass=1 kg, k=10 or 100 or 1000 N/m, A=0.1 m, and w=2*pi radians/s. I got the results shown in the plots below. The only thing that changes in the three plots below is the value of the spring constant, k.
k=10 N/m below
k=100 N/m below
k=1000 N/m below
If it can be done in a simple Excel spreadsheet, it can be done in Matlab.
Joerg Fricke-Schmidt
2023-11-8
William Rose
2023-11-8
The spring pulls upward when the mass displacement is <0. When k is larger, the spring pulls upward more strongly when the mass displacement is negative, therefore the mass follows the plate less closely when the plate displacement is negative.
Here is the way I think about the problem.
Compute the plate position xp(t), and velocty vp(t), without differential equations.
Mass position update: x(t+dt)=max(x(t)+v(t)*dt, xp(t+dt))
Note that x(t+dt)+v(t)*dt is the Euler method, also known as first order method. I added the "max(...)" to account for possible pushing by the plate.
Mass velocity update:
If the mass was not touching the plate on the previous timestep, then update the velocity by Euler method: v(t+dt)=v(t)+a(t)*dt, where a(t)=acceleration=-(k/m)*x(t)
Otherwise, update the velocity by whichever is greater: the Euler above, or the current plate velocity.
The code below implements the ideas above.
% define constants
A=0.1; % plate displacement amplitude (m)
w=2*pi; % plate frequency (radian/s)
m=1; % mass (kg)
k=100; % spring constant (N/m)
dt=0.001; % time step (s)
tmax=3; % max time (s)
% compute plate position, velocity
t=0:dt:tmax; % time vector
N=length(t);
xp=A*sin(w*t);
vp=A*w*cos(w*t);
% compute mass position (m), mass velocity (m/s)
x=zeros(1,N); v=zeros(1,N);
for i=2:N
x(i)=max(xp(i),x(i-1)+v(i-1)*dt); % position update
if x(i-1)>xp(i-1)
v(i)=v(i-1)-k*x(i-1)*dt/m; % velocity update
else
v(i)=max(v(i-1)-k*x(i-1)*dt/m,vp(i)); % velocity update
end
end
% display results
figure;
plot(t,x,'-r.',t,xp,'-b');
grid on; xlabel('Time (s)'); ylabel('Position (m)')
legend('Mass','Plate');
titlestr=sprintf('Mass-Spring-Plate, A=%.2f m, k=%.0f N/m',A,k);
title(titlestr)
Forward Euler method, used above, uses the first order Taylor series approximation. The error at each time step is proportional to the time step squared, and the cumulative error is (potentially) proportional to the time step size. A fourth order method such as ode45() uses a fourth order Taylor series approximation, so the error at each time step is proportional to the 5th power of the time step. Therefore ode45() can acheive a much higher level of accuracy than first order, with the same time step size. Forward Euler can also be unstable, i.e. it can lead to solutions that oscillate and grow to infinity, if the step size is too large. Backward Euler is stable, but sometimes it is not practical to compute the necessary equations. These drawback of the Euler method are reasons why people use ode45(). But ode45() does not deal well with derivatives that are discontinuous or infinite, which can happen when you have if() statements and max() or min() functions in the derivatives. One can reduce the time step size to a tiny value, to achieve a desired level of accuracy with a first order method. (If you are unsure, divide the time step size by 2, and re-run. If the solution did not change significantly, then you need not reduce the time step size further.) With today's fast computers with large amounts of memory, this approach is often easier and faster than using a more elegant high-order method.
Good luck with your work.
Joerg Fricke-Schmidt
2023-11-9
William Rose
2023-11-9
@Joerg Fricke-Schmidt, you're welcome.
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
2023-11-7
2023-11-9
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
