The Coppock solution is a classical analytic solution of internal ballistic equations. This model is found in the classic reference "Theory of the Internal Ballistics of Guns", by J. Corner (1950). The equations are based on the following assumptions:

- The powder mass burn rate is proportional to the chamber pressure and the surface area of the powder.
- The surface area of the powder is computed as it burns.
- The powder burns at the breech pressure.
- The increasing volumes behind the bullet with burning time and bullet movement is taken into account.
- The pressure gradient between breech and bullet base is computed and assumed to be constant depending on the charge weight to bullet weight ratio.
- The powder energy content, geometry, gas heat capacity and ideal gas corrections are used in the solution.
- Frictional and recoil losses are estimated by raising the bullet mass by 5%.
- Heat loss to the barrel is estimated by the Thornhill method and used in the solution.

In order to predict the rate of gas production from a burning propellant grain, we need to make some assumption about how the powder grain burns, both geometrically and thermodynamically. If the burning is a surface effect, the burn rate might be proportional to the number of molecules hitting the surface. This would mean it would be proportional to the pressure and temperature. The temperature of the burning is determined by the energy of combustion and the heat capacity of the products and would be constant. The burn rate could then be expressed as a rate constant relating how fast the flame front burns in a single dimension per unit pressure. A more general expression incorporates a pressure exponent parameter that need not be unity. Piobert postulated that the burning of a powder grain occurs in parallel layers from the surface to the interior. This has been shown to be a pretty good approximation by observation of partially burned grains of powder. The rate of gas production is then related to the surface area of the grain, the rate constant and pressure. The surface area may change with the burning depending on the geometry of the grain. If we designate the minimum burning dimension of the powder grain as "D" and "f" as the fraction of "D" that is left unburned, the reduction of the dimension "D" at any time is (1-f)*D. The burn rate is therefore -D*df/dt which is equal to the linear burn rate constant (beta) times the pressure (P).

-D*df/dt = beta*P (1)

Now the effect of the grain shape must be taken into account. The mass of gas produced will be equal to the mass of powder that has burned which would be proportional to the volume of the powder grain that has burned. If we say "z" is the mass or volume fraction of of the powder burned (or gas produced), then

z = (initial volume - vol. at time "t")/(initial volume)

Now "z" can be calculated as a function of "f" for various shapes. A tubular shape with wall thickness "D" and average radius "R", and length "L":

initial volume = 2*pi*R*D*L

at "t", the wall thickness is D-(1-f)*D = f*D

the length is L-(1-f)*D , so

vol at "t" = 2*pi*R*f*d*(L-(1-f)*D)

and

z = (1-f)(1+f*D/L)

This led to defining a general form function

z = (1-f)*(1+theta*f) (2)

where theta is the form factor and can range from -1 to +1. A form factor of 0 implies a constant area burn surface, such as with a tubular powder with a length much greater than the wall thickness. Long solid cylindrical stick powder such as cordite has a form factor of 1. Multiperforated powder actually increases in surface area with burning and has a negative form factor. Ball powder would follow the equation z = (1-f)^3 and is not easily modeled by the quadratic equation 2. Flake or disk powder has a cubic term but can be approximated by equation (2) with a theta of D/(L/2) where L is the diameter of the disk or length of the flake. The quadratic form of equation (2) allows the analytic solution of the model which would be much more difficult or impossible if we included a cubic term. There are more serious difficulties with modeling ball powders than the form function, however, such as the use of burn rate retardants.

A progressive powder is one with a negative theta, that is the burning area increases with burn time. A degressive one is the opposite, like cordite or spherical powder, where the surface area decreases substantially with burn time. Retardants are additives that are diffused into the powder grain from the outer surface that reduce the initial burn rate of the powder. This is has the effect of rendering the powder more progressive, and is used especially with ball powders to make them less degressive. A progressive powder has the advantage of producing a higher average pressure to peak pressure ratio than degressive powders, thus giving higher bullet speeds for a given peak pressure. In order to model the action of the retardant it would be necessary to introduce "beta" as a function of "f" and solve the system numerically. Coppock's model cannot handle the effect of retardants except by adjusting the value of "theta" and "beta".

One could point out that with the computer power available now, there is no reason to constrain our assumptions to those necessary to solve the equations analytically, and simply do them numerically with a step by step process. This is quite true, but the numeric solution has its own set of difficulties to overcome and a analytic solution is very useful in testing the convergence and accuracy of a numeric model. This analytic solution has been found to offer useful solutions and is a good starting point, so it has been presented here. It also has the advantage of being faster to compute so it can be programmed into a graphing pocket calculator, for example, which would be too slow with a numeric model.

The equation of motion of the shot (bullet) can be written immediately:

Ps(x)*A = w*dv/dt (3)

The pressure at the base of the shot (as a function of x, the shot travel distance) times the area of the shot base would give the accelerating force which equals the mass (w) of the shot times the acceleration.

In freshman chemistry we learned the ideal gas equation:

PV=nRT

which relates pressure, temperature, volume and with the number of moles of gas "n" and the gas constant "R". The equation works best for small molecules at low pressures. Some correction is needed at the high pressures experienced in a gun, so we will add the "covolume" term (b), a simplified version of Van Der Waal's semi-empirical ideal gas corrected equation called the Alfred-Nobel equation:

P*(V-b) = nRT

If we divide both sides by the mass of the gas we get:

P*(V-b)/m = R*n/m*T

If V and b are in units of volume per unit mass of gas the equation becomes

P*(V-b) = R/mw*T

where mw is the molecular weight of the gases (mass per mole). The covolume can be thought of as a correction for the finite volume of the gas molecules.

If the powder is burning at constant volume with no heat loss, the gas equation yields

P(V-b) = R/mw*T0

where T0 is the product gas temperature. The right hand side is a constant that depends on the molecular weight and the gas combustion temperature and is designated as the Force constant "F". This is also called Impetus by some powder companies.

F = R/mw*T0

and has units of energy per unit mass.

The term adiabatic implies a boundary which heat cannot cross. In other words all the energy of an adiabatic expansion is provided by or given to the gas, and results in the gas gaining or losing temperature according to the definition of heat capacity:

E = n * Cv * (T0 - T)

where E is the energy absorbed or provided by the gas, Cv is the heat capacity of the gas at constant volume ( the gas does no work at constant volume) times the temperature difference. "n" is the amount of gas so Cv would be heat capacity per unit amount of gas. The above equation written in differential form would be:

dE = n*Cv*dT

Since work is equal to force times distance or pressure times volume:

dE = P*dV

and n*Cv*dT = -P*dV

The minus sign indicates that the temperature must drop when doing the expansion work. Substituting the ideal gas equation gives:

n*Cv*dT = -nRT*dV/V which separating the variables gives

Cv*dT/T = -R*dV/V

It can be shown that for an ideal gas, Cp - Cv = R, where Cp is the heat capacity at constant pressure. If we apply heat to a gas at constant pressure, it must expand to keep the constant pressure, so it is doing pressure-volume work - which explains why Cp is greater than Cv. The ratio of Cp/Cv is commonly referred to as "gamma". Substituting gamma we have:

R/(gamma-1) = Cv

which is really the definition of gamma in the non-ideal gas case. After integrating we have:

ln(T_{2}/T_{1}) = -(gamma-1)*ln(V_{2}/V_{1}) or

T_{2}/T_{1} = (V_{1}/V_{2})^{(gamma-1)} or from the ideal gas
law

P_{1}*V_{1}^{gamma} = P_{2}*V_{2}^{gamma}

The above results are important when solving the internal ballistic problem after the powder is burnt, which is basically an adiabatic expansion case. Actually there is heat loss to the barrel, but that can be accounted for by adjusting gamma.

The integral of P*dV is the energy of an expansion. From above

P = P_{i}*V_{i}^{gamma} * V^{(-gamma)}

where P is the pressure at volume V and Vi, Pi represent the initial volume and pressure. Integrating P*dV from Vi to V gives the energy of expansion to volume V.

E = integral(P_{i}*V_{i}^{gamma}*V^{(-gamma)}*dV
) from V_{i} to V

= P_{i}*V_{i}/(1-gamma)*[(V/V_{i})^{(1-gamma)} -1 ]

The work done by the gas in the internal ballistic problem goes into heating the barrel, and doing pressure-volume expansion work, part of which results in the acceleration of the bullet. The work done by the gas can then be equated in terms of the heat capacity and temperature difference from above to the pressure volume work and heat loss:

C*z/mw*Cv*(T0-T) = A*integral(P*dx) + Eh

C is the powder charge, z is the mass fraction of gas, so C*z = m, the mass of gas.

m/mw = n.

This energy is equal to the pressure volume work done by the gas plus the heat lost to the barrel (Eh). The Pressure times area (A) is force, and the integral of force over the distance is work. Substituting our equation of state P*(V-b)=RT/mw, eliminating T.

C*z/mw*Cv*T0 = C*z*Cv*P*(V-b)/R + A*integral(P*dx) + Eh

substituting for Cv from:

(gamma-1) = R/Cv

C*z*R/mw*T0/(gamma-1) = C*z*P*(V-b)/(gamma-1) + A*integral(P*dx)+Eh

and introducing force constant F=R/mw*T0

F*C*z/(gamma-1) = P*C*z*(V-b)/(gamma-1) + A*integral(P*dx) + Eh

V is the volume per unit mass, so C*z*V is the total volume occupied by the gas, which is equal to the chamber volume "K" minus the volume of the solid propellant plus the volume of the bore behind the bullet A*x. So if "d" is the density of the propellant, then

C*z*V = K + A*x - C*(1-z)/d

So the energy balance equation becomes:

F*C*z/(gamma-1) = P/(gamma-1)*(K+A*x-C/d-C*z*(b-1/d))+A*int(P*dx)+Eh

Since K-C/d is the initial free space in the case, we can write

K-C/d = A*L, where L is the effective length of the free space

in the chamber as if it were the same cross sectional area as the bore. Then we can write:

- F*C*z/(gamma-1) = P/(gamma-1)*(A*(x+L)-C*z*(b-1/d)) +

A*int(P*dx) + Eh

If the Eh is neglected and the pressure P in the expansion integral is assumed to be the pressure at the base of the shot, the resultant equation is known as Resal's equation.

The pressure in Equation 4 should be taken as the average pressure since it was derived from the equation of state of the gases. We cannot equate this pressure with the pressure at the base of the shot to use with the equation of motion of the shot, since there will be a pressure gradient from the breech of the gun to the base of the shot due to the inertia of the powder gases. Lagrange solved this problem by considering a volume element at position "sigma" where sigma is a fraction of the distance from the breech to the shot base. In this manner the sigma (designated here as "s") is independent of time. If V is the volume at any moment between breech and shot base (the volume of the powder is ignored), then s*V will be the volume behind the cross section at "s". A volume element at "s" will have a mass of p*V*ds where p is the density (rho) of the gas. The velocity of the element will be s*v where v is the velocity of the shot. The momentum of the volume element is then s*v*p*V*ds. The time derivative of the momentum gives the inertial force which is equal to the force from the pressure on the section so:

d(p*V*v*s*ds)/dt = -A*dP/ds * ds or

s*d(p*V*v)/dt = -A*dP/ds

where P is the pressure. Integrating this from s at pressure P to the shot base at pressure Ps, s=1, gives:

A*(P-Ps) = integral( d(p*V*v)/dt *s*ds )

Now substitute the equation of motion A*Ps=w*dv/dt for dt

(P-Ps)/Ps = integral( d(p*V*v)/dv *s/w*ds)

Since the pressure gradient is small in proportion to the absolute pressure, the variation of the density "p" with s is neglected so p*V = C*z and integrating gives:

(P-Ps)/Ps = C/(2w)*d(z*v)/dv*(1-s^2)

If the above equation is solved for the breech pressure Pb, where s=0, and the result combined with the original equation:

P-Ps = (1-s^2)*(Pb-Ps)

The average pressure Pm is determined by the integral

Pm = integral( P*ds) from s=0 to s=1.

If this is done using P from the above equation, the result is:

Pm = 1/3*(2*Pb + Ps)

Also from above when solving for Pb at s=0 z=1 and v=0

Pb/Ps = (1-C/(2*w))

thus giving the relationship between average, breech and shot base pressure:

Ps = Pb/(1+C/(2*w)) = Pm/(1+C/(3*w))

This approximation neglects the fluid friction of the gas which may be important at high velocities, and also neglects the gas density variation, the powder volume and movement, and is really derived based on an all burnt charge. The extra friction towards the muzzle cancels the errors near the breach to some extent, so the net result of the approximation is pretty good. The kinetic energy of the gases can be found by integrating the kinetic energy of volume elements over s similar to what was done above so:

KEgas = integral( 1/2*p*V*s^2*v^2*ds ) from 0 to 1

which using p*V=C*z again gives:

KEgas = 1/6*C*z*v^2

Let us consider the expansion energy term in Equation 4. This energy represents the kinetic energy and rotational energy of the projectile, the frictional losses, the kinetic energy of the gas, losses due to recoil of the gun and even the strain energy in expanding the gun barrel. Another factor is the bullet engraving energy which may be appreciable especially in pistols compared to the available powder energy. Coppock's model assumes that of 5% of the energy of the bullet goes into friction and recoil (1 or 2 % effect). The other losses are very small. The rotational energy is about 0.5%, and the strain energy is less than a percent and are neglected. The friction is taken into account by increasing the mass of the bullet by the 5%. This implies that the friction is proportional to the force on the bullet base, which is not too good an assumption since there would be a component that would be fairly independent of the pressure or position of travel down the bore. The assumption is used since it makes the problem solvable analytically. The engrave force is not taken into account in this model and would be an important factor to include in a more advanced model. The kinetic energy of the powder gases was shown to be 1/6*C*v^2, (at all burnt) so combining these terms:

A*int(P*dx) = 1/2*(w*1.05)*v^2 + 1/6*C*v^2

These terms can be combined into a single kinetic energy term:

1/2*w1*v^2 , where w1 = 1.05*w + 1/3*C

Next, adding Eh, the heat loss to the barrel in Resal's equation:

F*C*z/(gamma-1) = P/(gamma-1)*(A*(x+L)-C*z*(b-1/d)) + 1/2*w1*v^2 + Eh

If we play a similar game with Eh, and write it as a fraction "k" of the kinetic energy term

Eh = k*1/2*w1*v^2

then

F*C*z = P*(A*(x+L)-C*z*(b-1/d)) + (gamma-1)*(1+k)*1/2*w1*v^2

Thus, if we define a new gamma, "gammap" that would be in essence the effective gamma, we could retain the original form of the Resal's equation if:

(gammap-1) = (gamma-1)*(1+k)

From above, k = Eh/Es where Es = 1/2*w1*v^2

Then

(4) F*C*z = P*(A*(x+L)-C*z*(b-1/d)) + (gammap-1)*1/2*w1*v^2

This manipulation implies that the heat loss at any instant is equal to a constant fraction (k) of the shot energy. This assumption has been shown to be good to around 20%. This does not imply that the error in calculating the velocities and pressures are this large, however, since the net energy balance might be quite good. The problem is that to compute

Es = 1/2*w1*v^2

we need to know the final bullet velocity "v". The program achieves this by iterating the calculation several times, using the computed final velocity as an estimate each time until convergence is achieved. The value of Eh, the total heat loss to the barrel is calculated using Thornhill's semi-empirical formula which is a function of the powder charge, (T-T0), and the bore dimensions. Thornhill's equation is the only empirical estimate used in this program.

Thornhill derived a formula based on test firings combined with theoretical results for the total heat lost to the barrel (cal) of a gun as:

Eh = 10.13*H*d*Vol/A

where H=.0127*T*sqrt(d)*R

and T is the maximum temperature of the barrel

d is the bore diameter (in)

A is the area of the bore (in^2)

Vol is the total volume of the gun (in^3)

R is a hydrodynamic roughness factor which ranges from 1.25 in big guns to 1.4 or more in small arms

T = (T0-300)/(1.7+0.38*d^.5*(d^2/C)^0.86)

where ambient is 300 deg K

C is the charge in lbs

This reduces to the following in ft*lbs units:

Eh = 0.397*d^(3/2)*Vol/A*T*R

There are 4 basic equations that describe the system in Coppock's model:

1) The powder burn rate equation.

-D*df/dt = beta*Pb (1)

2) A powder form function that relates mass of gas produced to the linear fraction of the powder grain that has burned.

z = (1-f)*(1+theta*f) (2)

3) The equation of motion.

Ps(x)*A = w*dv/dt (3)

4) The equation of state or energy balance.

(4) F*C*z/(gammap-1) = Pm/(gammap-1)*(A*(x+L)-C*z*(b-1/d)) +

1/2*w1*v^2

It is interesting to consider the effect of including these terms. Equation 4 contains the term (b-1/d). If b = 1/d then the whole term drops out and the equation is simplified considerably. Many early treatments of the problem have made this assumption. In reality, the covolume "b" is around .95 cc/g and 1/d is around 0.62 cc/g and its net effect would be to increase the effective force constant of a simplified equation by about 10%. Note that if one were to include the effect of the increase in volume due to the burning powder solid, and neglect covolume, that the result would be worse than neglecting both effects. The effect of the initial volume displaced by the solid powder is an important effect in either case.

Equation 1 and 3 by substituting for Pb and Ps from the Lagrange relation and integrating yields:

(5) v = A*D*(1-f)/(beta*(w+1/2*C))

giving velocity as a function of fraction burned "f". At "all burnt", f=0 and the velocity can be read directly. The final velocity can be computed from the result of the following "adiabatic" expansion, but one needs to know the "x" value at all burnt in order to figure the expansion ratio from that point. The other 3 equations can be combined and simplified to a differential equation with unitless coefficients by defining the following:

b_av is the breech to average pressure ratio

b_av=(1+C/(2.*w*1.05))/(1+C/(3.*w*1.05))

M = A*A*D*D/(beta^2*b_av^2*w1*F*C)

theta1 = theta + 1/2(gammap-1)*M

Z = 1 - (gammap-1)*M/2 + theta1*f

B = (b-1/d)/(A*L)

This is the resultant diff. eq.

(6) Z*dx/df + M*(x+L) = M*B*L*z

Integrating this gives "x" as a function of "f", from which we get "x" at all burnt (f=0) or "xb".

Equation 4 and 5 give:

(7) A*Pm*(x + L*(1-Bz)) = C*F*Z*(1-f)

giving the pressure as a function of "x" and "f". Combining (6) and (7) and taking dP/df = 0 gives the maximum pressure, but this gets a bit involved. The details are not presented here, but the complete formulas are available in the Mathcad(TM) example. The only thing left is to compute the adiabatic expansion from the all burnt point. The adiabatic expansion equation gives:

P = Pb*(V/Vb)^(-gammap)

where V is the total volume behind the bullet and Vb is the volume at all burnt.

V/Vb = (x + L*(1-B))/(xb + L*(1-B))

and the final velocity is:

v^2 = C*F/w1*(M + Zb*Phi)

where Phi = 2/(gammap-1)*(1-(V/Vb)^(1-gammap)),

and Zb is Z at all burnt.