Issues Using Dirac Function Trying to add Impulse Force at specific time

I am working on a project to simulate a car moving across a bridge as a coupled Vehicle-Bridge system and I am trying to add a collision with a rock at the midpoint of the bridge. The collision is modeled as an impulse load and I have been trying to use dirac delta to apply the impulse load at the desired time when the front of the car, where the impact occurs, reaches the middle of the bridge. The vehicle speed is constant, so this time is known.

I am modeling the impact load as fimp = Io*dirac(t-timp); and plugging it into my external force vectors where required from my derivation, but when I run the code, there is no change in the output that was modeled previously, so the impulse force is not being included. I know dirac delta only gives and answer when the input is zero, so for my scenario when t = timp, so I imagine that may be the issue where the times output by ode45 do not perfectly line up with my impact time, so there is never a value when the input is zero. Is there a way to get this to work as I am intending, or am I using dirac incorrectly and need to take a different approach?

MATLAB Code:

%System Parameters
a = 4.5; M = 2.5e3; I = 3.2e2;
k = 5e5; c = 3.6e3; v0 = 18; g = 9.81;
rho = 8e2; EI = 7e8; L = 400;

%Impact Parameters
Io = 1.2e3;
timp = ((L/2)-2*a)/v0;

%Parameter Matrices
%Beam
Mm = rho*L*[1/630 1/2772; 1/2772 1/12012];
Km = EI/(L^3)*[4/5 6/35; 6/35 2/35];

%Rigid Body
Mr = [M 0; 0 I];
Dr = c*[1 -a; -a a^2]+c*[1 a; a a^2];
Kr = k*[1 -a; -a a^2]+k*[1 a; a a^2];

%Admissible function inputs
syms t
Phix1 = [((v0*t/L)^2)*(1-(v0*t/L))^2 ((v0*t/L)^3)*(1-(v0*t/L))^3];
Phix1T = [((v0*t/L)^2)*(1-(v0*t/L))^2; ((v0*t/L)^3)*(1-(v0*t/L))^3];
Phix2 = [(((v0*t+2*a)/L)^2)*(1-((v0*t+2*a)/L))^2 (((v0*t+2*a)/L)^3)*(1-((v0*t+2*a)/L))^3];
Phix2T = [(((v0*t+2*a)/L)^2)*(1-((v0*t+2*a)/L))^2; (((v0*t+2*a)/L)^3)*(1-((v0*t+2*a)/L))^3];

%Dirac delta functions
fimp = Io*dirac(t-timp);

%Combined
Dc = c*Phix1T*Phix1+c*Phix2T*Phix2;
Dbr = c*Phix1T*[1 -a]+c*Phix2T*[1 a];
Drb = transpose(Dbr);

Kc = k*Phix1T*Phix1+k*Phix2T*Phix2;
Kbr = k*Phix1T*[1 -a]+k*Phix2T*[1 a];
Krb = transpose(Kbr);

%Assembly Matrices
Ma = [Mm zeros(2,2); zeros(2,2) Mr];
Ca = [Dc -Dbr; -Drb Dr];
Ka = [Km+Kc -Kbr; -Krb Kr];
Fa = [0; 0; -M*g+fimp; a*fimp];

%State Representation Matrices
Minv = Ma^-1;
p = [0; 0; 0; 0; Minv*Fa];
Ident = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
A = [zeros(4,4) Ident; -Minv*Ka -Minv*Ca];
Asimp = simplify(A);
Afunc = matlabFunction(Asimp);
Pfunc = matlabFunction(p);

%ODE Solver
tspan = [0 (L-2*a)/v0];
iniCon = [0; 0; .035; 0.015/a; 0; 0; 0.15; -0.45/a];
[t,z] = ode45(@(t,z) Afunc(t)*z + Pfunc(t), tspan, iniCon);

%Solve for y1 and y2
y1 = z(:,3)-a*z(:,4);
y2 = z(:,3)+a*z(:,4);

%Plots
subplot(2,1,1)
plot(t,y1)
grid on
title('Displacement of y_{1}(t)')
xlabel('Time (s)')
ylabel('Displacement (m)')

subplot(2,1,2)
plot(t,y2)
grid on
title('Displacement of y_{2}(t)')
xlabel('Time (s)')
ylabel('Displacement (m)')

Answer 1

[Paul]

Hi Michael,

The ODE solvers aren't set up to handle Dirac impulse functions this way.

As I read the code (make sure to double check this), the Dirac impulse only gets into the code through the third and fourth terms of Fa. So:

a. delete fimp from the equations altogether (so you won't need Pfunc, because p will be a double)

b. itegrate the equations using ode45 from t = 0 to t = timp.

c. Now, at timp the impulse hits and the states would be updated by

znew = z + [0; 0; 0 ;0 ;Minv(timp)*[0; 0; Io; a*Io] ] (I think these are the only terms that include the impulse, check that)

the idea being to integrate through the Dirac from timp(-) to timp(+)

d. integrate again using ode45 over the interval from timp to tfinal specifying znew as the initial state

e. combine the z and t outputs from steps b and d into a single array for plotting.

[Michael]

This definitely makes more sense as a better approach, for the znew term, would it just be the last row of the first ode45?

My new script has z1 as a 5377x8 matrix, so I set:

znew = transpose(z1(5377,:)) + [0; 0; 0; 0; Minv*[0; 0; Io; a*Io;]];

And then used znew as the initial conditions for my second ode45 from timp to tfinal. It does give a very slight jump in the displacement plot, as expected, I just want to make sure that the znew I used was what you were describing.

[Paul]

Correct. I should have been more precise.

Also, you can generalize so that the code will work regardless of the size of z1 (like if you change timp or the ode45 settings, for example)

znew = transpose(z1(end,:)) + [0; 0; 0; 0; Minv*[0; 0; Io; a*Io;]];

I should also point out that this approach assumes that each element of Minv*[0; 0; Io; a*Io;] is continuous at t = timp.

[Michael]

Great, thank you so much for your help!

[Walter Roberson]

Note also that in theory Dirac Delta is not a function but rather a "distribution". The interval is infinitesimally wide and the integral over the interval is defined as 1, so the "value" over the interval is 1/delta and as delta goes to zero the value goes to infinity.

Accordingly, when you call dirac() as a function the value for any specific input is either 0 or infinity, with int() knowing the special properties of the distribution. Symbolic integral of dirac() works, numeric gets you 0 or inf, which is Not Helpful.

So even if you managed to put a "waypoint" in to get ode*() to integrate at the place where the input was exactly 0, it wouldn't go well.

Conclusion: do not include dirac() in the numeric formulas. And heaviside() is going to have boundary problems too.

[Paul]

I would say that symbolic integral int of dirac works when it's multiplied against a function that is continuous at the point where the Dirac impulse is applied. Strange things can happen when the function is discontinuous (thread).