Time-dependent system ODE with initial half-sine pulse shows incorrect results
2 views (last 30 days)
Show older comments
Hello,
I have damped spring-mass system with initial half-sine pulse (acceleration) h loaded from a base.
Motion of equation should be the following:
m*d2(y)/d2t + c*d(y)/dt + k*y = k*(h)dtdt + c*(h)dt
k*(h)dtdt = stiffness*double integrated acceleration h = stiffness*displacement [N/m*m]
c*(h)dt = damping_coeff*integrated acceleration h = damping_coeff*velocity [Ns/m*m/s]
*Previous part has been edited for clarification.
Here are my concerns:
1. I am not sure how to introduce here the acceleration pulse in the equations,
2. for integration of k*(h)dtdt + c*d(h)dt I am using cumtrapz function. Is it OK?
Best regards
Mathew
t1=0.5;
a=0; h0=1;
h = @(t) a*(t>t1)+ h0*sin(pi.*t./t1).*(t<=t1);
ODExvf = @(t,y,m,D,K,h) [y(2);
-D/m*y(2)-K/m*y(1)+ cumtrapz(cumtrapz(t,h(t)))*K/m + cumtrapz(t,h(t))*D/m + h(t)];
m=1;
D=10;
K=100;
tspan=0.01:0.05:5;
F0=[0 0];
[t,F]=ode45(@(t,y) ODExvf(t,y,m,D,K,h),tspan,F0);
plot(t,F);
0 Comments
Answers (2)
William Rose
on 25 Feb 2023
I would not use cumtrapz() or (worse) nested calls to cumtrapz(cumtrapz()).
I don;t understand the right hand side of the equation
m*d(y)/dt + c*d(y)/dt + k*y = k*(h)dtdt + c*(h)dt
It seems you use dt*dt after k*h to get the units of force, and likewise you use dt after c*h to get to units of force. It is always nice when the units are correct, but I don't know physically what this means. Why is the acceleration multiplied by the spring constant k and by the dashpot constant c? Maybe your system supplies as much force as necessary to acheive a specific acceleration. If so, then during that period, you don;t need to solve a complicated differential equaiton - you can just differentiate once and twice to get velocity and position. When the pulse disconnects, then you solve the unforced system, with a zero on the right hand side.
Also, I assume you meant to write
m*d^2(y)/dt^2 on the left hand side.
3 Comments
William Rose
on 26 Feb 2023
THank you for clarifying the equation.
Now that it is clear, a problem is also clear: each term on the left side has units of force. The terms on the right side have units of displacement and velocity. Or maybe displacement and displacement*. Either way, all terms on both sides must have the same units.
*You wrote, in your clarification, "c*(h)dt = once integrated displacement = velocity", but I suspect you meant to write "c*(h)dt = once integrated velocity = displacement".
I still am not sure I understand what you are trying to simulate. I think you have a mass attached to a base by a spring and a dashpot in parallel. The base undergoes a half-sine-wave displacement pulse. This brief movement of the base will set the mass into motion (if it was previously stationary), and it will disrupt the motion of the mass (if it was already moving).
However, you said in your comment "There is no acceleration of the base...." This statement is inconsistent with the problem as I described it. If the base goes through a half-sine-wave of motion then it DOES accelerate. Therefore my understanding of the system is wrong, or your statement that there is no acceleration of the base is wrong.
William Rose
on 26 Feb 2023
If I am correct that the base moves according to a prescribed function of time, then the differential equation can be written
where y is the mass position and y_b is the base position. Your code will include the half-pulses yb=A*sin(wt) and vb=dyb/dt=A*w*cos(wt). To make these functions zero after an initial half-pulse, do this:
yb=A*sin(w*t).*(t<pi/w);
vb=A*w*cos(w*t).*(t<pi/w);
function dydt = myode(t,y)
m=1; c=1; k=1; % physical system properties
A=1; w=1; % amplitude, frequency of base displacement
yb=A*sin(w*t).*(t<pi/w); % displacement pulse
vb=A*w*cos(w*t).*(t<pi/w); % velocity pulse
dydt=[y(2);[figure out this part]];
end
That gets you started.
12 Comments
William Rose
on 27 Feb 2023
Edited: William Rose
on 27 Feb 2023
I like doing it with yb(t) and vb(t) defined as funtions. But I understand that you want to do it by numerical integration of the acceleration pulse. You asked why the most recent code you posted does not work. Here things I would try:
ic() must have 4 elements, not 2.
dydt=zeros(4,1) not zeros(4:1)
Eliminate yb from odefcn() and from the main script.
The y vector is:
y=[yb;vb;y;v];
which looks recursive, but you know what I mean: y and v on the right side are the mass position and velocity. Rewrite the formulas in odefcn(t,y) accordingly.
Define m, c, k inside odefcn().
William Rose
on 27 Feb 2023
For example:
function dydt=odefcn(t,y)
% Next: Unpack y, makes the code below easier to understand
yb=y(1); vb=y(2); ym=y(3); vm=y(4);
m=1; c=1; k=1; % physical constants (kg, N-s/m, N/m)
G=1; t1=0.5; % base acceleration amplitude, duration (m/s^2, s)
w=pi/t1; % base acceleration frequency (rad/s)
accBase=G*sin(w*t).*(t<t1); % base acceleration
dydt = [0;0;0;0];
dydt(1)=vb; % d(yb)/dt
dydt(2)=accBase; % d(vb)/dt
dydt(3)=vm; % d(ym)/dt
dydt(4)=(-c*(vm-vb) - k*(ym-yb))/m; % d(vm)/dt
end
See Also
Categories
Find more on Ordinary Differential Equations 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!