Time-dependent system ODE with initial half-sine pulse shows incorrect results
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(2 个)
William Rose
2023-2-25
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 个评论
William Rose
2023-2-26
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
2023-2-26
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 个评论
William Rose
2023-2-27
编辑:William Rose
2023-2-27
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
2023-2-27
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
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!