Internal ballistics optimization (2 variables)

Good day guys,

I am from the Czech republic so I would like to apologize for any grammar errors in advance. I am a student of guns and ammunition at the University of Defence and am currently working on my school project, in which I am trying to solve an optimization problem in the field of internal ballistics.

I put all my efforts (the code) bellow as well as a further explenation of what is going on, the overall goal and where I am currently at.

clc,clear,close

% PROBLEM DESCRIPTION
prob = optimproblem("Description",['Optimizing the mass of propellant charge ' ...
    'and thickness of powder grain']);

% VARIABLES
omega = optimvar('omega',1,'LowerBound',2.7e-3,'UpperBound',3.3e-3);
e1 = optimvar('e1',1,'LowerBound',0.21e-3,'UpperBound',0.26e-3);

% DEFINITION OF OBJECTIVE FUNCTION
f = fcn2optimexpr(@internal_ballistics,omega,e1);
p = f(1); 
v0 = f(2);

prob.Objective.first = p;
prob.ObjectiveSense = "min";
prob.Objective.second = v0;
prob.ObjectiveSense = "max";

% CONSTRAINTS
prob.Constraints.vazba1 = p<=300e6;
prob.Constraints.vazba2 = v0>=750;
prob.Constraints.vazba3 = omega<=5.72e-3;

% PROBLEM SOLUTION 
x0.omega = 2.71e-3;
x0.e1 = 0.21e-3;
[sol,optimval,exitflag,output] = solve(prob,x0);

% DISPLAY OF CALCULATION RESULTS
disp(strcat('Optimized mass of propellant charge omega=',num2str(sol.omega(1,1)),'kg'))
disp(strcat('Optimized thickness of powder grain e1=',num2str(sol.e1(1,1)),'m'))
disp(strcat('Resulting chamber pressure p=',num2str(optimval(1,1)*1e-6),'MPa'))
disp(strcat('Resulting velocity v0=',num2str(optimval(2,1)),'m/s'))

function f = internal_ballistics(omega,e1)

% Atmospheric conditions
pa = 1e5; % Ambient pressure [Pa]
Ta = 21+273.15; % Ambient temparature = temparature of propellant [K]

% Weapon parameters
d = 7.62e-3; % Calibre [m]
s = 4.73e-5; % Cross-section area of bore [m^2]
lu = 0.509; % Path of projectile inside the barrel [m]

% Cartridge parameters
% omega = 3.1e-3; % Mass of propellant charge [kg]
mq = 9.6e-3; % Mass of projectile [kg]
c0 = 3.521e-6; % Initial volume of combustion chamber [m^3]
kphi = 1.1; % Coefficient of passive resistance against shell motion [-]
del = omega/c0; % Loading density [kg/m^3]
phi = kphi+1/3*omega/mq; % Fictitious mass of projectile [-]

% Parameters of propellant charge
rho = 1627; % Density of powder mass (bulk density) [kg/m^3]
Tv = 3175; % Explosion temparature [K]
f = 0.73e6; % Specific energy of propellant [J/kg]
Lz = 1.79e-3; % Lenght of powder grain [m]
% e1 = 0.165e-3; % Thickness of powder grain [m]
k1 = 1+e1/Lz; % Coefficient of powder grain [-]
k2 =-e1/Lz; % Coefficient of powder grain [-]
k3 = 0; % Coefficient of powder grain [-]
Ik15 = 0.1702e6; % Pressure impuls of propellant gases [Pa.s]
ikt = 0.0016; % Temparature coefficient [-]
Tnpl = Ta; % Temparature of propellant [K]
Ik = Ik15*(1-ikt*(Tnpl-288.15));

% Thermodynamic parameters of propellant gases
alf = 0.906e-3; % Covolum [m^3/kg]
kappa = 1.2505; % Adiabatic exponent [-]
r = f/Tv; % Gas constant [J/kg/K]
theta = kappa-1; % Heat parameter of powder gas expansion [-]
cv = r/theta; % Specific heat at constant volume [J/kg/K]

% Other parameters
p0 = 4e7; % Ambient pressure [Pa]

%--------------------------------------------------------------------------

% Initial conditions
z0 = 0;                % Burnt thickness of powder grain
V0 = omega/rho;        % Volume of propellant charge 
m0 = pa*(c0-V0)/r/Ta;  % Mass of propellant gases (C0 = c0-V0)
xs0 = 0;               % Path of projectile 
vs0 = 0;               % Velocity of projectile
T0 = Ta;               % Temparature 
y0 = [z0 V0 m0 vs0 xs0 T0];
tspan = 0:2e-4:1;
options = odeset('RelTol',1e-6,'AbsTol',1e-6,'Event',@stopprogram);

[t,y]= ode45(@IntBal,tspan,y0,options);
zz = y(:,1);  % Instaneous burnt thickness of powder grain z = e/e1;
VV = y(:,2);  % [m^3]
mm = y(:,3);  % [kg]
v_s = y(:,4); % [m/s]
x_s = y(:,5); % [m]
TT = y(:,6);  % [K]

CC = c0-VV-alf*mm+x_s*s; % Instaneous change of volume in barrel [m^3]
pp = mm*r.*TT./CC; % Instaneous chamber pressure [Pa]

p_max = max(pp); % Maximum chamber pressure
v0 = max(v_s); % Maximum velocity of projectile
f = [p_max,v0];

save vnibal_7_62x59.mat

function dy = IntBal(~,y)
z = y(1);   
V = y(2);   
m = y(3);   
vs = y(4);  
xs = y(5);  
T = y(6);    

C = c0-V-alf*m+xs*s;
p = m*r*T/C; % (Equation of state)

% Additive constant "q" for determining f_p = q+p
% if p<20e6
%    q=35e6;
% elseif p>=20e6 && p<50e6
%        q=30e6;
% elseif p>=50e6 && p<100e6
%        q=25e6;
% elseif p>=100e6
%        q=20e6;
% end
f_p = p; % Refining pressure function

% Parameter "z" 
if z < 1
   c1 = 1;
else
   c1 = 0;
end

if p < p0 && vs==0
   c2 = 0;
else
   c2 = 1;
end

dzdt = c1*f_p/Ik;                               % Change of burnt thickness of powder grain
dVdt = -c1*V0*(k1 + 2*k2*z + 3*k3*z^2)*dzdt;    % Change of volume of propellant charge 
dmdt = -rho*dVdt;                               % Change of mass of propellant gases
dvsdt = c2*(p-pa)*s/(phi*mq);                   % Change of path of projectile
dxsdt = vs;                                     % Change of velocity of projectile
dCdt = -dVdt - dmdt*alf + s*vs;                 % Change of volume in barrel
dTdt = 1/(m*cv)*(dmdt*(cv*Tv-cv*T)-p*dCdt);     % Change of chamber pressure

dy(1) = dzdt;
dy(2) = dVdt;
dy(3) = dmdt;
dy(4) = dvsdt;
dy(5) = dxsdt;
dy(6) = dTdt;

dy = [dy(1);dy(2);dy(3);dy(4);dy(5);dy(6)];
end

function [value,isterminal,direction]=stopprogram(~,y) % stopping condition
    value = y(5) - lu;
    isterminal=1;
    direction=0;
end
end

Firstly, the code consists of internal ballistics model "function f = internal_ballistics(omega)", onto which is attached the optimization code done via optimization toolbox.

The overall goal is to find optimal values for 2 variables, mass of the propellant charge "omega" and thickness of powder grain "e1", and ofcourse satisfy the enclosed constraints. The ultimate idea is to obtain solution for those 2 variables, which will ensure minimizing the maximum pressure in the barrel "p_max", while at the same time maximizing the maximum velocity of the projectile "v0" but to be honest, I am not sure, if I am able to achieve this by using the optimization toolbox so I would greatly apreciate your opinion on that. What I am focusing on right now is just to find the 2 varibles and minimize "p_max".

So far, the optimization only finds optimal solution for 1 variable, which is "omega" and I do not know, where to go from here (finding those two at the same time).

I am open for discusion, I will appreciate any kind of help and I am open to provide additional information on the model itself if needed.

Thank you very much in advance and have great day

Ludvík Hladký

[Torsten]

If you could define a second optimization variable "e1", where would it enter the equations used for the computation of p and v0 ?

[Ludvík Hladký]

Hello Torsten

Thank you for your question ! I revisited the code and made accurate changes so "e1" enters the equations as it should. I am sorry for that. Now, if you look at the % Parametrs of propellant charge, there you can see "e1" being a part of the Coefficients of powder grains ("k1", "k2").

[Torsten]

clc,clear,close

% PROBLEM DESCRIPTION
prob = optimproblem("Description",'Optimizing the mass of propellant charge','ObjectiveSense',"min");

% VARIABLES
omega = optimvar('omega',1,'LowerBound',2.8e-3,'UpperBound',3.3e-3);
e1 = optimvar('e1',1,'LowerBound',0.165e-3-0.05e-3,'UpperBound',0.165e-3+0.05e-3);

% DEFINITION OF OBJECTIVE FUNCTION
f = fcn2optimexpr(@internal_ballistics,omega,e1);
p = f(1); v0 = f(2);
prob.Objective = p;

% CONSTRAINTS
prob.Constraints.vazba1 = p<=350e6;
prob.Constraints.vazba2 = p>=320e6;
prob.Constraints.vazba3 = v0<=1000;
prob.Constraints.vazba4 = v0>=700;

% PROBLEM SOLUTION 
x0.omega = 2.81e-3;
x0.e1 = 0.165e-3;
[sol,optimval] = solve(prob,x0);

Solving problem using fmincon.

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

sol.omega

ans = 0.0031

sol.e1

ans = 2.1500e-04

disp(strcat('Optimized mass of propellant charge omega=',num2str(sol.omega),'kg'))

Optimized mass of propellant charge omega=0.0030913kg

disp(strcat('Resulting chamber pressure p=',num2str(optimval*1e-6),'MPa'))

Resulting chamber pressure p=334.1865MPa

function f = internal_ballistics(omega,e1)

% Atmospheric conditions
pa = 1e5; % Ambient pressure [Pa]
Ta = 21+273.15; % Ambient temparature = temparature of propellant [K]

% Weapon parameters
d = 7.62e-3; % Calibre [m]
s = 4.73e-5; % Cross-section area of bore [m^2]
lu = 0.509; % Path of projectile inside the barrel [m]

% Cartridge parameters
% omega = 3.1e-3; % Mass of propellant charge [kg]
mq = 9.6e-3; % Mass of projectile [kg]
c0 = 3.521e-6; % Initial volume of combustion chamber [m^3]
kphi = 1.1; % Coefficient of passive resistance against shell motion [-]
del = omega/c0; % Loading density [kg/m^3]
phi = kphi+1/3*omega/mq; % Fictitious mass of projectile [-]

% Parameters of propellant charge
rho = 1627; % Density (specific) of powder mass [kg/m^3]
Tv = 3175; % Explosion temparature [K]
f = 0.73e6; % Specific energy of propellant [J/kg]
Lz = 1.79e-3; % Lenght of powder grain [m]
%e1 = 0.165e-3; % Thickness of powder grain [m]
k1 = 1+e1/Lz; % Coefficient of powder grain [-]
k2 =-e1/Lz; % Coefficient of powder grain [-]
k3 = 0; % Coefficient of powder grain [-]
Ik15 = 0.1702e6; % Pressure impuls of propellant gases [Pa.s]
ikt = 0.0016; % Temparature coefficient [-]
Tnpl = Ta; % Temparature of propellant [K]
Ik = Ik15*(1-ikt*(Tnpl-288.15));

% Thermodynamic parameters of propellant gases
alf = 0.906e-3; % Covolum [m^3/kg]
kappa = 1.2505; % Adiabatic exponent [-]
r = f/Tv; % Gas constant [J/kg/K]
theta = kappa-1; % Heat parameter of powder gas expansion [-]
cv = r/theta; % Specific heat at constant volume [J/kg/K]

% Other parameters
p0 = 4e7; % Ambient pressure [Pa]

%--------------------------------------------------------------------------

% Initial conditions
z0 = 0;                % Burnt thickness of powder grain
V0 = omega/rho;        % Volume of propellant charge 
m0 = pa*(c0-V0)/r/Ta;  % Mass of propellant gases (C0 = c0-V0)
xs0 = 0;               % Path of projectile 
vs0 = 0;               % Velocity of projectile
T0 = Ta;               % Temparature 
y0 = [z0 V0 m0 vs0 xs0 T0];
tspan = 0:1e-6:1;
options = odeset('RelTol',1e-13,'AbsTol',1e-13,'Event',@stopprogram);

[t,y]= ode45(@IntBal,tspan,y0,options);
zz = y(:,1);  % Instaneous burnt thickness of powder grain z = e/e1;
VV = y(:,2);  % [m^3]
mm = y(:,3);  % [kg]
v_s = y(:,4); % [m/s]
x_s = y(:,5); % [m]
TT = y(:,6);  % [K]

CC = c0-VV-alf*mm+x_s*s; % Instaneous change of volume in barrel [m^3]
pp = mm*r.*TT./CC; % Instaneous chamber pressure [Pa]

p_max = max(pp); % Maximum chamber pressure
v0 = max(v_s); % Maximum velocity of projectile
f = [p_max,v0];

%save vnibal_7_62x59.mat

function dy = IntBal(~,y)
z = y(1);   
V = y(2);   
m = y(3);   
vs = y(4);  
xs = y(5);  
T = y(6);    

C = c0-V-alf*m+xs*s;
p = m*r*T/C; % (Equation of state)

% Additive constant "q" for determining f_p = q+p
% if p<20e6
%    q=35e6;
% elseif p>=20e6 && p<50e6
%        q=30e6;
% elseif p>=50e6 && p<100e6
%        q=25e6;
% elseif p>=100e6
%        q=20e6;
% end
f_p = p; % Refining pressure function

% Parameter "z" 
if z < 1
   c1 = 1;
else
   c1 = 0;
end

if p < p0 && vs==0
   c2 = 0;
else
   c2 = 1;
end

dzdt = c1*f_p/Ik;                               % Change of burnt thickness of powder grain
dVdt = -c1*V0*(k1 + 2*k2*z + 3*k3*z^2)*dzdt;    % Change of volume of propellant charge 
dmdt = -rho*dVdt;                               % Change of mass of propellant gases
dvsdt = c2*(p-pa)*s/(phi*mq);                   % Change of path of projectile
dxsdt = vs;                                     % Change of velocity of projectile
dCdt = -dVdt - dmdt*alf + s*vs;                 % Change of volume in barrel
dTdt = 1/(m*cv)*(dmdt*(cv*Tv-cv*T)-p*dCdt);     % Change of chamber pressure

dy(1) = dzdt;
dy(2) = dVdt;
dy(3) = dmdt;
dy(4) = dvsdt;
dy(5) = dxsdt;
dy(6) = dTdt;

dy = [dy(1);dy(2);dy(3);dy(4);dy(5);dy(6)];
end

function [value,isterminal,direction]=stopprogram(~,y) % stopping condition
    value = y(5) - lu;
    isterminal=1;
    direction=0;
end
end

[Geoffrey Kolbe]

1) You need to find out how much powder will fit into the case. This is usually the bulk density of the powder times the case volume. I see that bulk density is not specified, so you may have to make an estimate of the bulk density based on the kernel geometry. That will give you the mass of the powder. Start off with a reasonable value of grain thickness and calculate the bulk density and so the mass of powder that will fill the case,

2) The velocity will be maximised if the case is filled with powder and the powder is all-burnt in the barrel (to maximise the energy from the powder). The pressure will be minimised if the powder is just all burnt when the bullet exits the barrel. Using the mass you get in 1), iterate the grain thickness until the powder is just all-burnt when the bullet exits the barrel.

3) Given the thickness of grain you get from 2), calculate from 1) if the given mass of powder is too little to fill the case or too much to fill the case? If it is too much, reduce the mass slightly and do step 2) again. If it is too little, increase the mass slightly and redo step 2) again.

4) Redo step 3) until the mass does not change and the case is full of powder. You will now have maximise the velocity and minimised the pressure.

[Ludvík Hladký]

Hello Geoffrey,

and thank you so much for your tips ! I followed the steps in 1) and came to realise that the bulk density actually coresponds to the density of powder mass "rho", which is specified in the code. It is rather difficult for me to keep up with all the signs and names for things because in literatures I have seen so far on the topic, every publisher just calls things as they fancy :-D.

Anyway, in the end, the maximum mass of powder, which would fit into the case, should be equal to ≈ 5.73 g (88.4 gr), which is waaay more that is currently needed. But I have used this value as yet another important constraint.

Finally, I think I have also managed to achieve the ultimate goal, which was to min(p_max) and simoultaneously max(v0) for two variables "omega" and "e1". I will update the original code shortly.

[Ludvík Hladký]

Hello again Torsten,

I am not sure, how to respond directly to your comment but I hope you see this. The code you provided seemed a lot less difficult than I have anticipated it to be. I am going to update the original code shortly so you can see where I am at.

The code as it is right now uses genetic algorithm for a solver and gives me the results for optimal values of "omega" and "e1" and also min(p_max) and max(v0).

What I would like to achieve now is displayed on the example picture bellow from SW ANSYS Fluent:

When I used to calculate approximate drag coefficent for bullets, I very much liked that in the "Console", it shows the amount of iterations and their respective values for parameter in each coresponding iteration. Moreover, It also provided graph of the progress of parameters dependent on the iterations.

I would like to ask, if you could help me to transfer this into my example, where in the Command window I would primarly see iterations and for each iteration the change of p_max, v0, omega, e1.

Thank you for your patience and time in advance

[Geoffrey Kolbe]

Hello Ludvik

If you want to maximise the velocity, the case must be full of powder and all the powder must burn in the barrel. So that is two constraints.

You can alter the pressure by altering the grain thickness. The thicker the grain, the longer it takes to burn and so the lower the maximum pressure will be. However it the grain is too thick then the powder will not all be burnt when the bullet exits the barrel. So there is a maximum thickness to the grain so that all the powder burns before the bullet exits the barrel. That is the second contraint again.

However, the amount of powder you can fit in the case will depend on the dimensions of the grains. If rho in the program is the bulk density then what is the propellant density? You will need to know the propellant density so you know how much volume the powder takes up in the case at the start of the computation. You cannot get that from the bulk density. I think you need to check that again.