How to set up Runge Kutta Simulation for Numerical Vehicle Dynamics?

Hey all, I'm struggling to set up my RK4 equations to simulate some simple vehicle dynamics. Essentially, what I want to do is to simulate the time passed, distance travelled, velocity, velocity^2 and acceleration numerically for a vehicle. The vehicle will start from stationary and travel a straight track with varying inclines until it reaches an end distance of about 1000m.

Drive Force: Engine: The engine is either on and off (depending on distance along the track, program will check 'status' between each step calculation), and if the engine is on, the force provided will be constant without delay.

The resistance forces I will have are:

Air Drag: = 0.5 * Density * Drag Area * Velocity^2 (Only Velocity^2 will change, all other variables constant)

Rolling Resistance = Mass * Gravity * Coefficient of rolling resistance (All values will be constant)

Incline = Mass * Gravity * Sine(incline in rads) (The incline will be obtained from an array which holds data for distance along track and incline at that point. Data is for every metre, therefore incline only changes every 1m)

The equation of motion on the vehicle in the direction of travel will be as follows :

Ma = FEngine - FAirDrag - FRolling - FIncline

However I struggle to realise how to turn this into Runge Kutta 4th order equations that will calculate distance travelled, velocity, velocity^2 and acceleration at each step calculation, say h = 0.1s.

I've seen examples of how it's done when there is just one force, where they use separation of variables. However I can't really get my head around doing it when there is more forces. I would think it should not be too hard given that all forces apart from AirDrag are either constant, or calculated to a numerical value from a separate data array.

Any help would be appreciated! Many thanks in advance

Answer 1

[Jan]

It depends how your RK4 code is written. Remember that you have to convert the higher order equation to a system of equations of order 1 at first. In your case these are 2 equations, one for the velocity and one for the acceleration. All you "simulate" is the time, the position and the velocity. To get velocity^2 and acceleration you can simply evaluate the function with the formerly calculated position and velocity.

To do this, your RK4 must accept a vector as input. Because you did not post any code, it is hard to suggest any real code. But simply try to implement this advice and show us the code.

[dada5555]

Hi Jan, thank you for your answer. By watching a video online, I managed to set up RK4 code for velocity. But I'm struggling to add position in to it. I did the velocity equation of motion by replacing acceleration as dv/dt. However, I' m not sure how to do position within that RK4. I assume it has something to do with doing it in vector form like you suggested, but I'm not sure how to do that?

Here is my code:

function myrk4()

clear all

dt = 1; %Step Time

%%%INITIAL CONDTIONS
t(1) = 0; 
x(1) =0;
v0 = 0;

%%%%RK4 SOLVER
%Set Initial Conditions
v4(1) = v0;
t4(1) = 0;
n= 1;

% Solve until 50 seconds passed
while t4 < 200
t4(n+1) = t4(n) +dt; % Time 

%RK4 Weighted Averages
k1 = f(v4(n),t4(n));
k2 = f(v4(n) +k1*dt/2, t4(n) + dt/2);
k3 = f(v4(n) +k2*dt/2, t4(n) + dt/2);
k4 = f(v4(n) +k3*dt, t4(n) + dt);

phi = (1/6) * (k1 +2*k2 +2*k3 +k4); %Sum of RK4 Averages

%Calculates next velocity
v4(n+1) = v4(n) + phi*dt;

n= n +1; 

end

%%%PLOT velocity versus time
plot(t4,v4, 'b-', 'LineWidth',2)
hold on
end

%%%Equation of Motion in respect to dv/dt AKA vdot
  function vdot = f(v, t)
  %Set Constants 
  M= 90;              % Mass in Kg
  rho = 1.225;        % Air Density in kg/m^3
  gradient = 0;       % Incline kept as zero for simplicity atm
  Rolling = 0.003;    % Rolling Resistance
  g = 9.81;           % Gravitational acceleration
  DragArea = 0.05;    % Drag Area in m^2
  FEngine = 80;       % Engine Force, kept  constant for simplicity atm

    %%Equation of motion in respect to dv/dt
    vdot = (1/M)*(FEngine - 0.5*rho*DragArea*v^2) - (g*(Rolling + sin(gradient)));

end

%

[Torsten]

If you have the velocity, you only need to add a second ODE

d/dt (position) = velocity

to get the position.

[dada5555]

Thanks all, adding Torsten's ODE and implementing it in the right order did the trick!

[Ponmalar M]

input t in the function f(v,t) is not used? how it will be when we write for position?

[Jan]

@Ponmalar M: Usually the function to be integrated is written as f(t,y), where y is a vector containing the position and the velocity as y(1) and y(2).