Feasibility of kinetic orbital bombardment

To calculate a trajectory, the initial conditions must be determined first. They serve as a starting point for the re-entry simulation. The focus of this study is orbital re-entry, where a satellite in orbit around Earth manoeuvres the projectile so that it falls towards Earth’s surface. This scenario is shown in Figure 1. In order for the projectile to hit a target on the Earth’s surface, it needs to manoeuvre and enter the atmosphere and go through orbital re-entry. The orbital velocity, V_orbit, is given by 1Vorbit=μ⊕R,$${V_{{\rm{orbit}}}} = \sqrt {{{{\mu _ \oplus }} \over R}} ,$$ with μ_⊕ = M_⊕G, where G and M_⊕ are the gravitational constant and the mass of the Earth, respectively, and R is the distance to the centre of the Earth: 2R=R⊕+h,$$R = {R_ \oplus } + h,$$ where R_⊕ represents the radius of the Earth and h the altitude. Note that the orbital velocity is not a function of vehicle mass, only of the altitude. To change altitude, to ultimately strike the surface of the earth, the velocity of the vehicle has to change. That velocity change is called ΔV, see e.g. Curtis (2005) and Wiesel (2010). In the most energy-efficient manoeuvres, so-called Hohmann transfers, the ΔV is parallel to the flight path. In case ΔV is applied to increases the velocity, the altitude on the opposite side increases to a more elevated orbit. To deorbit a projectile, ΔV reduces the velocity, which causes the opposite side of the orbit to lower until the atmospheric drag decelerates the projectile further and it is captured. The amount of velocity change determines the steepness of the re-entry trajectory.

Fig. 1:

In order to induce this change in velocity the projectile needs to be equipped with a retro-rocket that thrusts the projectile in the opposite direction of the flight direction. This manoeuvre burns propellant, which has to be carried by the satellite. To estimate the amount of propellant necessary for the manoeuvre, the Tsiolkovsky rocket equation is used: 3ΔV=Ispg0lnmimf$${\rm{\Delta }}V = {I_{{\rm{sp}}}}{g_0}\ln \left( {{{{m_{\rm{i}}}} \over {{m_{\rm{f}}}}}} \right)$$ where I_sp is the specific impulse of the propellant and g₀ is the gravitational acceleration at the surface of Earth. The mass m_i is the initial mass, prior to the burn, i.e. the mass of the rod and the rocket motor, including the propellant. The mass m_f is the final mass, i.e. the mass of the rod and the rocket motor after the propellant has been burned.

Figure 2 displays a diagram of the re-entering projectile, with orange arrows representing the velocity vector and blue arrows portraying the forces acting on the projectile.

Fig. 2:

The re-entry model used includes aerodynamic drag. Once the initial conditions are acquired, they are put into the re-entry simulation, which is bound by several assumptions:

The scenario is two-dimensional.

The latter point can be understood as follows. To keep the projectile stable during flight, it would have to be fin-stabilised. Tailfins ensure that, if the projectile were to gain an angle of attack, the aerodynamic force moment on the tail counteracts this motion, thus stabilising its flight and limiting the angle of attack. Furthermore, the projectile is assumed to be unguided during the atmospheric flight.

When the projectile travels at a hypersonic velocity, through the atmosphere, the air molecules ionise and a plasma forms around the projectile. This plasma blocks Electromagnetic (EM)-transmission (Bond 1958), which is essential for communication. This means that no external input can be given to navigate the projectile. Another way to guide the projectile to its target is to use a seeker head, but building a seeker capable of withstanding the high temperature may not be feasible. Finally, flight control would be problematic, since, if the angle of attack were to become too large, the force moment on the projectile might cause it to break up in flight.

The equations of motion for the subsequent suborbital trajectory are given by Hankey (1988). 4aR.=Vsin(γ)$$\mathop R\limits^. = V\sin (\gamma )$$ 4bV.=−μ⊕R2sin(γ)−ρV2SCD2m$$\mathop V\limits^. = - {{{\mu _ \oplus }} \over {{R^2}}}\sin (\gamma ) - {{\rho {V^2}S{C_D}} \over {2m}}$$ 4cθ.=VRcos(γ)$$\mathop \theta \limits^. = {V \over R}\cos (\gamma )$$ 4dγ.=VR2cos(γ)−μ⊕VR2cos(γ)$$\mathop \gamma \limits^. = {V \over {{R^2}}}\cos (\gamma ) - {{{\mu _ \oplus }} \over {V{R^2}}}\cos (\gamma )$$

The parameters m and γ are the projectile mass and the flight path angle (see Figure 2), respectively. The parameters R, V, θ and γ are time derivatives, and conform to Newton’s notation. By numerically integrating these equations with respect to time, the altitude, velocity, pitch and distance can be calculated as a function of time. This is done using MATLAB SimuLink^®.

2.1.1

Aerodynamic properties of the projectile

The drag coefficient of the projectile depends on its shape and the Mach number. As no information about actual orbital projectiles is available, we have to make some assumptions about the projectile. Elbasheer (2014) proposed a length of 6.1 m and a diameter of 0.3 m with a length over diameter (L/D) of 20. In this paper, we consider a range of different projectiles, each essentially cylindrical and with the same L/D ratio with lengths between 0.1 m and 6.1 m. For a tungsten alloy projectile, this corresponds to a mass between about 0.03 kg and 7,600 kg.

Assuming the projectile has a flechette-like shape, a close resemblance can be found in the XM110 projectile with an L/D of 23. Braun (1973) established a database for small arms aerodynamics, which includes aerodynamic data for a XM110 projectile for both laminar and turbulent flow (see Figure 3). Owing to the high velocity of the projectile, it is easy to gravitate towards assuming a fully turbulent flow. The Reynolds number helps to predict the type of flow. The Reynolds number is the ratio of internal forces and viscous forces: a high Reynolds number is associated with turbulent flow and a lower Reynolds number with a laminar flow. The Reynolds number is defined as: 5Re=ρVLμ$${\mathop{\rm Re}\nolimits} = {{\rho VL} \over \mu }$$ where ρ represents the density of the fluid, V is the fluid velocity, L is the characteristic length and μ is the dynamic viscosity, see e.g. Hibbeler (2020). It depends on the temperature of the fluid and is defined as (Chapman and Cowling 1970; Zheng et al. 2016): 6μ=T288.153/2288.15+cT+cμ288.15$$\mu = {\left( {{T \over {288.15}}} \right)^{3/2}}{{288.15 + c} \over {T + c}}{\mu _{288.15}}$$ where T is the fluid temperature, c is a fluid specific parameter and μ_288.15 is the viscosity at a temperature of 288.15 K. Eq. (6) is known as Sutherland’s formula and c is the so-called Sutherland constant (Chapman and Cowling 1970). For air c = 110.4 K and μ_288.15 = 1.822 × 10⁵ Pa/s, the critical Reynolds number indicates a fluid transition from a laminar to a turbulent flow. For a viscous flow over an external surface the critical Reynolds number is Re_cr = 5 × 10⁵ (see e.g. Hibbeler 2020). In the simulations, the flow is turbulent above Re_cr and laminar below Re_cr. To model the behaviour, the temperature is calculated using the international standard atmosphere (Anderson 2016). The simulation uses a switch case where, depending on the Reynolds number, it takes data either from the turbulent flow graph or from the laminar flow graph.

Fig. 3:

By simulating two trajectories, with one assuming fully turbulent flow and the other considering both types of flow, the difference can be observed. To best visualise the effect, a trajectory with a long flight path through the atmosphere is chosen. Thus, the projectile experiences the longest time under the influence of aerodynamic drag. Normally the flow transitions at the critical velocity. This transition means that there is a grey area where it could be either turbulent or laminar. Furthermore, an initial altitude of 400 km and a projectile with L = 8 m are chosen.

The resulting trajectory, flight velocity and pitch can be seen in Figures 4a, 4b and 4c, respectively. As all three figures show, a minor difference can be seen for the different methods. The impact velocity deviates 2.72%, the pitch deviates 4.70%, the flight time deviates 0.28% and the range deviates 0.32%. This is included in the model, although the effect is small.

Fig. 4:

To describe the projectile penetration into concrete and steel a semi-hydrodynamic model was chosen, i.e. the 1D-model of Tate and Alekseevksii. Of course, this 1D-model is less advanced than, for instance, models based on cavity expansion theory. However, it is a classical model that ought to be capable of adequately estimating the penetration depth of long rod projectiles in steel and concrete. The characteristic equation of Tate and Alekseevskii (see e.g. Alekseevskii 1966; Tate 1969; Hohler and Stilp 1990; Walker 2021) is a modified Bernoulli equation: 7a12ρp(v−u)2+Yp=12ρtu2+Rt,$${1 \over 2}{\rho _p}{(v - u)^2} + {Y_p} = {1 \over 2}{\rho _t}{u^2} + {R_t}{\rm{,}}$$ where u is the penetrating velocity and v is the tail velocity of the projectile, and R_t and Y_p are, respectively, the dynamic penetration strengths of the target and the projectile. Eq. (7a) determines the penetration velocity u. This velocity can be explicitly determined because all parameters in this equation are known. In their model, Tate and Alekseevksii introduce a finite projectile length. During penetration the projectile is eroding. The equation for the deceleration of the projectile is: 7bdvdt=−YpρpL,$${{{\rm{d}}v} \over {{\rm{d}}t}} = - {{{Y_p}} \over {{\rho _p}L}},$$ where L(t) is the time dependent projectile length. The rate of change of the projectile length is given by: 7cdLdt=−(v−u)$${{{\rm{d}}L} \over {{\rm{d}}t}} = - (v - u)$$

Eqs (7a–c) are representative of the Tate and Alekseevksii model. These equations are integrated to obtain the penetration depth (P) and length (L) of the projectile as a function of time. The different penetration mechanisms of the Tate–Alekseevskii model are illustrated in Figure 5.

Fig. 5:

Evidently, the Tate–Alekseevskii model has limitations, for projectile impacts are 3D-phenomena. However, the penetration depth is for our purpose the most important parameter, and it can be described using the 1D-model. Moreover, the (finite) speed of sound in the material is not considered, see e.g. Walker (2021). The longitudinal speed of sound in metals is typically 5–6 km/s. As long as the projectile velocity is smaller than the speed of sound in the material, the effect can be considered to be negligible. For further limitations of the model, see e.g. Bavdekar et al. (2017, 2019).

For realistic projectile impact simulations, both oblique impact and the possibility of ricochet have to be considered. First, we will treat ricochet. When the impact angle is too high with respect to the normal of the surface, the projectile skips off the surface. This is called ricochet. It is a complex projectile target interaction, but Tate (1979) developed a simplified model to calculate the critical ricochet angle β_crit. The projectile will ricochet if the impact angle is larger than the critical ricochet angle. The equation is as follows: 9tan3(β)>23ρpv02YpL2+D2LD1+ρpρt$${\tan ^3}(\beta ) > {2 \over 3}{{{\rho _p}v_0^2} \over {{Y_p}}}\left( {{{{L^2} + {D^2}} \over {LD}}} \right)\left( {1 + \sqrt {{{{\rho _p}} \over {{\rho _t}}}} } \right)$$ where β is the angle relative to the normal, ρ_p and ρ_t are the mass densities of the projectile and target, respectively, v₀ is the impact velocity, Y_p is the dynamic penetration strength, and L and D are the length and diameter of the projectile, respectively.

Penetration models, such as the TA-model, implicitly assume perpendicular impact. However, usually, the projectile will hit the target with an oblique impact angle. Since the velocity vector is in the same direction as the body axis of the projectile, the penetration depth can be calculated by simply correcting for the impact angle: 10P=P0cosγ=P0sinβ$$P = {P_0}\cos \gamma = {P_0}\sin \beta $$ where P is the penetration depth, P₀ is the penetration depth if the projectile would have impacted normal to the surface, and is the impact angle, as visualised in Figure 6. Note that γ is the angle with respect to the surface, while β = 90° - γ is, as said, the angle relative to the normal.

Fig. 6:

When a rod impacts the ground, a small fraction of the kinetic energy is transferred into seismic energy. This fraction is known as the seismic efficiency k (E_seism = k E_k). If the seismic efficiency is known, an estimation can be made on the seismic effects by converting seismic energy to the Richter scale (Gutenberg and Richter 1955): 11M=23logEseism −4.8$$M = {2 \over 3}\left( {\log \left( {{E_{{\rm{seism\;}}}}} \right) - 4.8} \right)$$ where M is the earthquake magnitude according to the Richter scale. Güldemeister and Wünnemann (2017) conducted a quantitative analysis of impact-induced signals using numerical modelling. High velocity impacts were tested of spherical steel and iron meteorites into various materials. As a result of these tests, numerous seismic efficiencies were obtained for different target materials ranging between 1.37 × 10^-4 (for 100% saturated sandstone) and 3.39 × 10^-3 (for quartzite). These values for sandstone and quartzite will be the upper and lower bounds for the simulations.

First, impact simulations were conducted for steel and concrete targets (without the flight model). The vital difference between these materials is that for concrete, the dynamic penetration resistance is smaller than the dynamic penetration strength of tungsten (R_t < Y_p), which means that the projectile has both a rigid phase and an eroding phase, depending on the impact velocity (see Figure 5). For steel, the dynamic penetration resistance is greater than the penetration strength of tungsten (R_t > Y_p), meaning that the projectile only experiences eroding penetration when the impact velocity is above the critical velocity. If the impact velocity is below the critical velocity, no penetration occurs (see Figure 5).

In Figure 7, the normalised penetration depth, P/L, is calculated for a range of impact velocities v₀. The thick orange and blue lines indicate the minimum and maximum values from Table 2. What stands out in Figure 7 is that the simulated materials show very different behaviour; the graph for concrete shows an optimum while the graph for steel does not.

Fig. 7:

The full simulation considers the influence of the impact angles and impact velocity. Figure 8 shows the flight time of the trajectory. In this figure a significant difference can be distinguished in flight time. The higher flight altitudes result in longer flight times, obviously because the projectile travels along a longer path for a higher altitude than for lower altitude.

Fig. 8:

The resulting penetration depth is shown in Figure 9. The figure shows that the initial orbital altitude does not have a significant influence on the penetration depth. At a first glance, it may seem remarkable that for projectiles with a length of more than roughly 1 m, the penetration depth is less than the length of the projectile. However, this is a consequence of the relatively shallow impact. The orbital altitude has only a small influence on the impact parameters because the initial orbital velocity V₀ and ΔV required for deorbiting the projectile change with altitude (see Figure 10).

Fig. 9:

Fig. 10:

These impacts are shallow due to the ΔV for the projectiles in Figure 8 being calculated using a Hohmann transfer ellipse, i.e. a minimum ΔV trajectory, to a lower orbit with a perigee at an altitude of 15 km at roughly the opposite side of the Earth. Below that altitude, aerodynamic drag slows the projectile down sufficiently for it to impact the ground. In Figure 11, the flight path angle (γ) for a 0.56 m projectile is shown as a function of time. This shows that for different initial orbital altitudes, the projectile enters the atmosphere at a different angle. When the projectile starts to experience aerodynamic drag, the trajectories start to look similar, resulting in roughly the same impact parameters. These trajectories are all shallow, although insufficiently shallow for ricochet to occur. In the operational context of kinetic orbital bombardment and Prompt Global Strike, a lower orbital altitude would be advisable because of the low flight time.

Fig. 11:

Different trajectories for deorbiting the projectiles are possible, resulting in steeper impacts, but at the cost of a higher ΔV.

Earthquakes can cause significant damage to urban structures. To calculate the possible earthquake magnitudes caused by projectile impact, the magnitudes are calculated for different projectile lengths from an altitude of 400 km, using a minimal ΔV trajectory (113 m/s for h₀ = 400 km) for quartzite and sandstone.

Figure 12 shows a maximum magnitude of 2.5 on the Richter scale. According to the US geological survey (USGS 2010) more than 1.3 million earthquakes of this magnitude happen annually, and damage occurs above a magnitude of 4 or 5 (USGS 2022), meaning that the seismic effects caused by high velocity impacts are insignificant compared to the penetration effects.

Fig. 12:

The ΔV required to deorbit the projectiles significantly impacts the cost of a hypothetical orbital bombardment system. The mass-to-orbit per projectile for different values of ΔV (see Figure 10) and projectile lengths have been calculated to develop an estimate of the mass that a satellite or rocket has to carry. The mass-to-orbit is the mass of the projectile and the mass of the fuel necessary to provide a certain ΔV (see Eq. 3). Here, we neglect the mass of the rocket motor. The ΔV ranges from the minimum ΔV (113 m/s for h₀ = 400 km) to 750 m/s and an I_sp = 300 s, which is a typical value for modern solid rocket motors in vacuum. These results are illustrated in Figures 13a-13c.

Fig. 13:

In spaceflight, the mass of the payload is an indication for the price to bring such a payload into an orbit. In Figure 14, the payload mass for a 4″, a 1 m and a 6.1 m projectile is presented. For a minimum ΔV, the payload masses are, respectively, 0.036 kg, 34.64 kg and 7,862 kg. For reference, the SpaceX Falcon 9 can carry a maximum payload of 22,800 kg, and has a launch costs of $56.5 million in 2013 (SpaceX 2012). This means that 1 kg payload costs roughly $2,500 to launch to a Low Earth Orbit (LEO, 400 km). For the aforementioned projectiles, this means that the launch of one projectile to LEO costs about a $100, almost a $100,000 and close to $20 million, respectively, excluding the mass of the satellite and the empty mass of the retro-rocket. The price of a system is very important, and as presented, the price increase per length is significant. Therefore, it is important to make an estimation of the influence of the length and ΔV of a projectile. Since global coverage will require multiple satellites, presumably with multiple rods, these costs soon become prohibitive.

Fig. 14:

There are two alternative methods to deliver a payload to its target, other than the orbital bombardment that we consider in this work. It can also be dropped from a bomber or launched using an ICBM. In a comparison of the kinetic orbital bombardment system with a suborbital trajectory, the range of both trajectories has been matched. For the suborbital trajectory, we chose a burnout velocity of 6,750 m/s and a launch elevation of 30°, which are typical values for an ICBM. These parameters result in a range of approximately 6,700 km, depending on the projectile length. To match this range of the suborbital flight, the re-entry uses 435.49 m/s of ΔV and an altitude of 400 km. The simulated bomber was flying at a cruising altitude of 15 km with a cruising speed of 1,000 km/h.

For these flight trajectories, the impact velocity and impact angle can be seen in Table 3, for projectile lengths between 0.1016 m (i.e. 4″) and 6.1 m. For the impact angle, 0° is horizontal. For the impact velocity, the lower boundary corresponds to the shorter projectile and the higher boundary corresponds to the longer projectile. For the impact angle, the opposite is true.

The penetration depth P for concrete and steel is illustrated in Figures 14a and 14b, for concrete (left) and steel (right). For concrete targets, the bomber shows better performance for projectiles larger than 0.3 m for both the re-entry and suborbital flight. Below 0.3 m, the suborbital flight shows better performance. The better performance for the bomber is caused by the lower impact velocity since an optimum value of P/L was found for a relatively low velocity (see Figure 7) and steeper impact angle. For steel targets the suborbital methods show better performance. This is mainly caused by the higher impact velocity in combination with a high impact angle.

From these results, it is apparent that an orbital bombardment system is, in principle, less effective than using a conventional method. For steel targets, launching the projectiles with an ICBM gives a better performance and for concrete targets, using a bomber shows better performance. This is not unprecedented. The US used so-called ‘Lazy Dog’ kinetic projectiles dropped by aircraft in the Korean and Vietnam Wars, see e.g. Karmes (2014). More recently, in the prelude to the 2003 invasion of Iraq, the US used concrete-filled laser guided bombs, to limit collateral damage (Copp 2003a, 2003b). An obvious downside of using a bomber is that it has to fly to its destination to deliver the payload and is more susceptible to intercepts.