Computing velocity with Forward Euler method

Hi,

I'm trying to compute the velocity evolution in time of an object being pushed by expanding gas. I'm using a forward euler with a time step of 1e-7 s. While the position evolution that I obtain seem smooth and with physical sense the velocity has abrupt changes and it is 'uncorrect in value'.

The initial conditions are nil position and velocity. p1 = 12 bar, p2 = 0 bar

This is my code so far:

clc
clear
close all

% Parameters  
n = 1; %polytopic exponent
V = 50e-3; % [m^3] tank volume
gamma = 1.4; 
b = (2/(gamma+1))^ (gamma/(gamma-1));           % critical ratio
R = 287;           % [J/kgK] Gas constant
T = 273.15+20;           %[K] Temperature
d = 48e-3;        % [m] bullet diameter
A = pi*d^2/4;      % [m^2] Area
cd = 0.7;
K = 0.868 * cd * A * sqrt(293.15/ (T^2*R));  
Ct = n*R*T/V;    
m = 1;           % [kg] Mass
dt = 1e-7;          % [s]  Time step
N = 1e6;           % Number of steps
% Preallocate
p1 = zeros(1, N);
p2 = zeros(1, N);
x = zeros(1, N+1);   % x(i+1) needed
G = zeros(1, N-1);
v = zeros(1, N);
% Initial conditions
p1(1) = 12*1.01325e5; % [Pa] initial tank pressure
p2(1) = 0; % [Pa] initial barrel pressure
x(1) = 1e-3; % [m] initial bullet position (not exactly zero otherwise impossible division)
v(1) = 0.001;
x(2) = x(1) + v(1)*dt;
for j = 2:N
    ratio = p2(j-1) / p1(j-1);
    if ratio <= b
        G(j-1) = K * p1(j-1);
    else
        G(j-1) = K * p1(j-1) * sqrt(1 - ((ratio - b)/(1 - b))^2);
    end
    % Update p1(t_i)
    p1(j) = p1(j-1) - (G(j-1) / Ct) * dt;
    % Update p2(t_i)
    v(j) = (x(j) - x(j-1)) / dt;
    p2(j) = p2(j-1) + ((G(j-1) * R * T - p2(j-1) * A * v(j)) / (A * x(j-1))) * dt;
    % Update x(t_{i+1})
    x(j+1) = 2*x(j) - x(j-1) + (p2(j) * A / m) * dt^2;
end
% Plot the results
figure;
subplot(2,1,1)
plot(x, 'LineWidth', 1.5)
xlabel('Step #')
ylabel('Bullet Position [m]')
title('Bullet Position')
grid on
subplot(2,1,2)
plot(v, 'LineWidth', 1.5)
xlabel('Step #')
ylabel('Bullet Velocity [m/s]')
title('Bullet Velocity')
grid on

I think the problem is in how I set the position x(2) but I'm struggling to find a way to fix it. If anybody has any suggestions it would be very helpful.

Thanks in advance,

Riccardo

Answer 1

[Torsten]

@Riccardo

Are you sure you can leave temperature fixed with a pressure difference of 12 bar ?

And did you try plotting your results against time:

plot(dt*(0:N),x,'LineWidth', 1.5)
plot(dt*(0:N-1),v,'LineWidth', 1.5)

? You'll see that your results become complex-valued very soon. And plotting complex numbers with plot(x) or plot(v) as you do gives weird plots.

From the documentation for "plot":

plot(Y) plots Y against an implicit set of x-coordinates.

• If Y is a vector, the x-coordinates range from 1 to length(Y).
• If Y is a matrix, the plot contains one line for each column in Y. The x-coordinates range from 1 to the number of rows in Y.

If Y contains complex numbers, MATLAB® plots the imaginary part of Y versus the real part of Y. If you specify both X and Y, the imaginary part is ignored.

clc
clear
close all
% Parameters  
n = 1; %polytopic exponent
V = 50e-3; % [m^3] tank volume
gamma = 1.4; 
b = (2/(gamma+1))^ (gamma/(gamma-1));           % critical ratio
R = 287;           % [J/kgK] Gas constant
T = 273.15+20;           %[K] Temperature
d = 48e-3;        % [m] bullet diameter
A = pi*d^2/4;      % [m^2] Area
cd = 0.7;
K = 0.868 * cd * A * sqrt(293.15/ (T^2*R));  
Ct = n*R*T/V;    
m = 1;           % [kg] Mass
dt = 1e-7;          % [s]  Time step
N = 1e6;           % Number of steps
% Preallocate
p1 = zeros(1, N);
p2 = zeros(1, N);
x = zeros(1, N+1);   % x(i+1) needed
G = zeros(1, N-1);
v = zeros(1, N);
% Initial conditions
p1(1) = 12*1.01325e5; % [Pa] initial tank pressure
p2(1) = 0; % [Pa] initial barrel pressure
x(1) = 1e-3; % [m] initial bullet position (not exactly zero otherwise impossible division)
v(1) = 0.001;
x(2) = x(1) + v(1)*dt;
for j = 2:N
    ratio = p2(j-1) / p1(j-1);
    if ratio <= b
        G(j-1) = K * p1(j-1);
    else
        G(j-1) = K * p1(j-1) * sqrt(1 - ((ratio - b)/(1 - b))^2);
    end
    % Update p1(t_i)
    p1(j) = p1(j-1) - (G(j-1) / Ct) * dt;
    % Update p2(t_i)
    v(j) = (x(j) - x(j-1)) / dt;
    p2(j) = p2(j-1) + ((G(j-1) * R * T - p2(j-1) * A * v(j)) / (A * x(j-1))) * dt;
    % Update x(t_{i+1})
    x(j+1) = 2*x(j) - x(j-1) + (p2(j) * A / m) * dt^2;
end
max(imag(x))

ans = 1.1956e-04

max(imag(v))

ans = 0.0014

% Plot the results
figure;
subplot(2,1,1)
plot(dt*(0:N),[real(x);imag(x)], 'LineWidth', 1.5)
xlabel('Time [s]')
ylabel('Bullet Position [m]')
title('Bullet Position')
grid on
subplot(2,1,2)
plot(dt*(0:N-1),[real(v);imag(v)], 'LineWidth', 1.5)
xlabel('Time [s]')
ylabel('Bullet Velocity [m/s]')
title('Bullet Velocity')
grid on