Time-dependent system ODE with initial half-sine pulse shows incorrect results
2 views (last 30 days)
Show older comments
Mathew Smith on 25 Feb 2023
Commented: William Rose on 27 Feb 2023
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?
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)];
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.
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 on 27 Feb 2023
% 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
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!