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

2 views (last 30 days)
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
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)];
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);

Answers (2)

William Rose
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.
Mathew Smith
Mathew Smith on 26 Feb 2023
Edited: Mathew Smith on 26 Feb 2023
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.

Sign in to comment.

William Rose
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:
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]];
That gets you started.
William Rose
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

Sign in to comment.


Find more on Programming 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!

Translated by