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:

Answer 1

[William Rose]

@Joerg Fricke-Schmidt,

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.

[Joerg Fricke-Schmidt]

Hello William,

thanks a lot for your answer.

1. Yes, my system sketch does not fit to the first order equation. That's the point. I'm stuck in adapting it to physical problem I want to solve.
2. The plate displacement is sinusoidal over time with a constant frequency. I need to calculate if at max (constant) frequency and the weakest spring mass condition the mass will start loosening the contact to the surface. And yes the mass will oscillate with the same driving frequency as long as the contact force between mass an driving plate is >0 under dynamic condition. That is what I need to calculate.
3. Yes the sinusoidal plate discplacement is independent of the load on it! That's the point where I stuck to put it into equations on a spring-mass-oscillator.
4. Sorry for beeing not specific enough with my equation. I have the damping term in it but the damping coefficient is "0". The system I want to calculate is undamped.
5. My equation is not yet considering pre-loaded condition. I do not know yet how to put that into a equation of motion. The equation was my start to describe the problem but now I'm stuck. The "A" term is a force driving the mass. But this will not work for my problem since the plate will drive the mass independend from the load on it. And this is exactly my question. How is a displacement driven mass described. Maybe I cannot use the equation of motion at all and need another approach.

I hope I could myself clear a little bit better. I will read your linked example. I would appreciate if you have any other idea to describe that physical problem.

Thanks a lot. Regards, Joerg

[William Rose]

@Joerg Fricke-Schmidt,

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]

Hello William,

seems that's what I was looking for. I will become comfortable with the Euler approach again. I can barely remember...

Any idea why the mass is less following the moving down displacement when k becomes larger?

Regrads, Joerg

[William Rose]

@Joerg Fricke-Schmidt,

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.

[William Rose]

@Joerg Fricke-Schmidt,

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]

Hey William,

thanks a lot for the great update. I was already starting to calculate this in Excel based on your earlier advices. Your new details will help me a lot to go forward. I will update the calculation with the pre-loaded condition. This should impact the motion of the mass below "0" position.

Thanks you very much!

Regards,

Joerg

[William Rose]

@Joerg Fricke-Schmidt, you're welcome.