Time-dependent system ODE with initial half-sine pulse shows incorrect results

6 ビュー (過去 30 日間)
Mathew Smith
Mathew Smith 2023 年 2 月 25 日
コメント済み: William Rose 2023 年 2 月 27 日
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);
Error using matlab.internal.math.getdimarg
Dimension argument must be a positive integer scalar within indexing range.

Error in cumtrapz>getDimArg (line 91)
dim = matlab.internal.math.getdimarg(dim);

Error in cumtrapz (line 49)
dim = min(ndims(y)+1, getDimArg(dim));

Error in solution (line 7)
-D/m*y(2)-K/m*y(1)+ cumtrapz(cumtrapz(t,h(t)))*K/m + cumtrapz(t,h(t))*D/m + h(t)];

Error in solution (line 14)
[t,F]=ode45(@(t,y) ODExvf(t,y,m,D,K,h),tspan,F0);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
plot(t,F);

回答 (2 件)

William Rose
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
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.
Mathew Smith
Mathew Smith 2023 年 2 月 26 日
編集済み: Mathew Smith 2023 年 2 月 26 日
Hello William,
the terms on the right side have units of force: k*x [N/m*m] and c*v [Ns/m*m/s]. I have clarified it above. You understood correctly, that there is base excitation. However, it is not excited by half-sine-wave displacement, but half-sine-wave acceleration. I can do it with displacement pulse (although my solution is more complicated than your suggestion). I am not able to include there half-sine-wave acceleration base excitation.
I am sorry for not clear description of the problem and I thank you for your effort.

サインインしてコメントする。


William Rose
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
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
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 ExchangeProgramming についてさらに検索

製品


リリース

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by