Performance of Lunar Regolith Shield under Meteoroid Impact: Uncertainties of a Numerical Prediction

Abstract

The protection of the future lunar base from radiation, thermal stresses and impacts of meteoroids can be achieved by several options, generally consisting in a shielding system, made of either regolith or combined materials. In the present paper, the incidence of two sources of uncertainty on the FEM calculation of stress propagation through the covering regolith layer has been assessed. First, the investigation has pointed out, for a given impact and a given constitutive model, the uncertainty in the stress prediction related to the strain interval adopted for the soil stiffness measurement. Thereafter, calculation has been performed, for a given stiffness value, changing the assumed collision duration of one order of magnitude for equal impact momentum, that is, changing the maximum impact force too, correspondingly. The simulation has been performed based on physical and mechanical parameters of DNA-1A lunar simulant. The results provide indication of the relative importance of the calculation assumptions, which could address the design of a regolith shield.

1. Introduction

The protection from extreme temperatures, solar/cosmic radiation and impacts of meteoroids will greatly condition the design of the future permanent moonbase [1,2,3]. This goal can be achieved by several options, generally consisting in a shielding system, for example using an inflatable habitat [3,4] as a support on which a dome can be constructed. By exploiting local material, an efficient radiation shielding is achievable by manufacturing structures with a sufficient wall thickness [5]. Under this perspective, using the D-Shape printing technology [6] with the lunar simulant DNA-1, Cesaretti et al. [5] demonstrated that the granular material may be processed to create bricks useful to build an outpost on the Moon. Fateri et al. [7], using JSC-1A and JSC-2A simulants, investigated their suitability for a sintering process. Nowadays, the researchers’ interest towards these manufacturing processes is growing and further technological solutions are proposed [8].

Although the main habitat of a lunar base, i.e., the pressurized habitat for crew, will benefit from design solutions based on achievements of the additive manufacturing for processing regolith, the shield of temporary or secondary structures, such as underground silos and systems for power storage, could be realized adopting more simple solutions. For example, layers of regolith deposited and compacted by rovers remain a suitable option for covering and protecting temporary structures. Therefore, theoretical, experimental, and numerical studies should still be addressed at modelling the stress propagation, related to impacts, through regolith earthworks and the penetration of solar radiation and heat within them.

Due to the extreme environmental conditions, on site adjustment of the construction procedure will not be an option. It means that the design must be based on highly reliable prediction of the static and dynamic response of the geotechnical works, to preserve the astronauts’ safety. For the numerical simulation of the construction steps, physical and mechanical parameters experimentally deduced on lunar regolith simulants may be implemented in the calculation codes, provided that all uncertainties are clearly stated.

When dealing with the numerical simulation of an impulsive load on a regolith bed, such as the load resulting from high kinetic energy impact of meteoroid on the regolith shield of a lunar habitat, three main sources of uncertainty can be identified:

• The characteristics of the impact, in terms of force-time history.
• The geotechnical properties of the soil, controlling the load propagation through the layer.
• Uncertainties related to the calculation model.

In a rough and preliminary way, the maximum impact force related to a collision can be deduced from the momentum under the simple assumption of perfectly inelastic impact, provided that the duration of collision and the shape of the force-time function are known. A theoretical approach for modelling the impact load acting on a structure led to analytical solutions [9,10] for the impact force and duration by considering a projectile striking a simply supported beam. Gionet [11] compared these analytical solutions with results from nonlinear finite element analysis, to investigate the structural response of a frame-membrane lunar habitat. However, a meteoroid striking a regolith layer is quite different problem and requires, to be fully described, size and morphology of the resulting crater. Furthermore, the dynamics of a projectile penetrating sand [12] should not be ignored.

When assessing the mechanical performance of the regolith shield, a “design impact” must be selected, based on probabilistic considerations concerning the meteoroid flux [13,14,15,16], i.e., concerning mass and velocity of the impactor and thus the impact momentum, and on physical considerations concerning the collision duration and the shape of the force-time function. Several experimental works have been carried out, mostly under terrestrial gravity and for low to medium kinetic energy, to measure the collision time, defined as the interval between the initial contact and the moment the impactor comes to rest. This information has been obtained, for example, by direct measurements of the acceleration of objects impacting granular media [17] (Figure 1). Efforts have been made by researchers to model the drag force related to impacts [18,19]. Scaling laws for penetration depth and for collision duration have been obtained by Goldman and Umbanhowar [17]. According to the authors, collision time t_c is independent of collision velocity v_c for sufficiently high v_c. Ricchetti at al. [20] logged in the laboratory, by means of piezoresistive sensor, the stress variation with time at the bottom of a lunar simulant bed during low kinetic energy impact due to free fall of sphere, an information indirectly related to the collision time. In addition, strain measurements can give information on the duration of an impact: Rickman et al. [21], looking at the micrometeoroid/orbital debris impact detection for spacecrafts, used fiber Bragg gratings to measure the in-plane strain fields of an aluminium plate, and the related time history, caused by projectile traveling at hypervelocity speeds.

Deducing a collision time for meteoroid impact on the Moon is however still a challenge, also due to the difficulty of in situ recording of the impulse history during real impact events. The seismic excitation caused by meteoroid impact may indirectly provide information on the collision: Gudkova et al. [22], based on Apollo-recorded impacts, found relationships between the momentum transfer of the impactor and the duration of the related seismic excitation process. Indirect information on the impact duration and magnitude may also derive from telescope observation, considering duration and brightness of the associated flash, although the light curve duration is not strictly a collision time.

Another source of uncertainty is the set of geotechnical parameters of the employed simulant. Although similar to the real lunar soils in terms of composition, the reliability of simulants is necessarily limited. First, the range for geotechnical parameters of the real lunar regolith is not so narrow, whereas a specific simulant possesses a unique grain size distribution. Secondly, even assuming perfect correspondence between real and simulated regolith in terms of granulometry and mineralogy, a terrestrial laboratory generally does not reproduce the space environment. According to Zhou et al. [23], lunar submicroscopic iron (SMFe), as one of the main products of space weathering, is a good indicator of maturity. It is however difficult to reproduce the space weathering on Earth and, on the other hand, it is questionable that this weathering reflects on the mechanical properties of lunar soils in a relevant way. More easily, a laboratory assessment of simulant reliability could concern the response of the material to thermal cycles characterized by temperature variations such as those experienced by the regolith on the Moon surface. In this respect, permanent bulk mechanical changes could derive from fatigue of the single particles, if they are prone to crushing or fracturing. Even dealing with a “thermally aged” simulant, under different temperatures the material can exhibit different elastic moduli, the reason why assuming a stiffness value measured at 20 °C implies some uncertainties in the numerical prediction for real lunar soils. Moreover, for a given temperature, it is well known that an elastic modulus exhibits a decay (Figure 2) with increasing the strain interval over which the stiffness measurement is made [24,25,26]. From an engineering point of view, this problem is relevant because the designer needs to assume the expected service strain regime of the geotechnical structure before calculation.

Numerical simulation of the effects of impacts on granular bed have been made using Discrete Element Method (DEM), for example by Torres-Cisneros et al. [27]. Ho and Masuya [28] approached the calculation of impact load by rockfall on energy absorbing sand-cushioning layer using Finite Element Method (FEM). DEM-FEM combined approached for simulating an impact problem has been employed too, such as in the study from Marzulli et al. [29]. In the present paper, the incidence of two sources of uncertainty on the FEM calculation of stress propagation through the regolith layer has been assessed. First, the investigation has pointed out, for a given impact and a given constitutive model, the uncertainty in the stress prediction related to the strain interval adopted for the soil stiffness measurement. Thereafter, calculation has been performed, for a given stiffness value, changing the assumed collision duration of one order of magnitude for equal impact momentum, that is, changing the maximum impact force too, correspondingly. The simulation has been performed based on physical and mechanical parameters of DNA-1A lunar simulant [30,31], having same origin of the simulant DNA created by Monolite U.K. [5], which is suitable to mimic low-Ti mare soil [32].

2. Material

A detailed description of the compositional features of DNA-1A can be found in Marzulli and Cafaro [30], who also studied its mechanical behaviour by element testing. At the particle scale, a micromechanical study on this simulant was carried out by Sandeep et al. [33]. The repose angle of the material, which is about 40°, reflects the grain size, shape, and roughness. The shear strength parameters measured during direct shear tests seem to vary depending on the normal stresses range employed for their determination: at low confining pressure, i.e., less than 10 kPa, Marzulli and Cafaro measured on self-weight compacted specimens a friction angle of 56° and a cohesion of 2–3 kPa, whereas at pressures higher than 10 kPa they found null cohesion and 44–47° friction angle. It should be underlined, however, that the mechanical properties deduced for simulants under terrestrial conditions could differ from those characterizing the material in the lunar environment [34]. For lunar grains, van der Waals forces need to be considered [35]. Experimental data obtained by Bromwell [36] and Nelson [37] under ultrahigh vacuum (10–7 Pa) and high temperature (394 K) showed an increase friction angle of 13° and an increase cohesion of 1.1 kPa testing simulated lunar soil. Increased cohesion of lunar simulant in high vacuum is also mentioned in Johnson et al. [38,39]. More experimental effort is desirable to detect the effect of temperature and high vacuum on the elastic moduli of lunar soils, that is at present an important source of uncertainty.

The dilatancy angles of lunar simulants, when not measured during laboratory tests, can be deduced as a function of the difference between peak and critical state friction angles and based on empirical equations linking Relative Density, D_r, and mean stress at failure to dilatancy [40], as made by Alshibli and Hasan [41], who proposed a statistical model for JSC-1A based on triaxial compression tests. Huang and Zheng [42], for example, discuss this calculation approach for a simulant at low confining pressure.

Marzulli and Cafaro [30] measured the Young modulus of DNA-1A at different confining pressures and over different axial strain interval (

Table 1. Young modulus values of DNA-1A for different confining pressures and strain intervals of measurement (Replotted with permission from Marzulli and Cafaro ©2019 American Society of Civil Engineers. Permission conveyed through Copyright Clearance Center, Inc.).

E (kPa)	p’ (kPa)	Axial Strain Interval (%)
42,065	158	0.17
19,814	102	0.25
9537	102	0.56
7828	67	0.53
7183	60	0.27
6057	43	0.20
34,129 a	41	6.2094 × 10−3
33,077 a	21	6.01343 × 10−3

^a With local strain gauge.

). The modulus greatly depends on both factors. Values of stiffness of other lunar simulants reported in the literature may be very different, probably depending on the measurement technique and adopted strain interval. For example, Colwell et al. [43] mention Young modulus of about 8 MPa at D_r = 66% for MLS-1 [44] and of 65–110 MPa at D_r = 60% for JSC-1 [45]. Since the elastic moduli of the material play a role in the load propagation through the layer, it is crucial to measure the material stiffness over proper strain interval when modelling impulsive phenomena.

Under low kinetic energy impact in the laboratory, realized by free fall of sphere, the lunar simulant DNA-1A has revealed a crater morphology with ray system (Figure 3a,b). This phenomenon has been already studied, for example by Sabuwala et al. [46]. Numerical simulations for granular material performed by Kadono et al. [47], who studied the crater-ray formation, showed that clear mesh pattern appeared at lower coefficients of restitution between particles but was less clear at larger one, suggesting that the inelastic collisions between particles cause the clear mesh-pattern formation of impact ejecta.

3. Methods

The problem here analysed is sketched in Figure 4. A regolith covering is supposed to absorb the impact of meteoroids and prevent the buried structure from damages. In this work, a thickness of 5.4 m is assumed for the regolith layer [14], which would provide protection from radiation and meteoroids of maximum mass of 37 kg, or diameter of 52 cm or smaller [48], reducing the flux of a penetrating impact on the structure to between 10⁻⁸ and 10⁻⁷ impacts/km²/yr [13,15]. Lower thickness could be employed in case of composite shields, adopting regolith in combination with other materials. The Relative Density of the regolith shield is assumed to be D_r = 79%. A second 50 cm thick layer immediately below the regolith cover identifies the vault of the habitat to be shielded. Assuming that the latter is manufactured by 3D-printing [5] using in-situ regolith, the parameters assigned to it in the numerical simulation (“lunar” unit weight, γ = 3.6 kN/m³; Young modulus, E = 2,350,000 kN/m²; Poisson ratio, ν = 0.2) reflect a level of stiffness that could plausibly be imparted to the material for satisfying structural requirements. This layer has been modelled as Linear Elastic Non-Porous in the FEM code. Furthermore, the conventional bedrock condition under the entire model was adopted.

The geotechnical characteristics of the regolith shield layer employed in the FEM calculation (

Table 2. Geotechnical parameters of the regolith layer employed in the FEM calculation.

Friction angle	54°
Cohesion	0
Dilatancy angle	20°
Earth pressure coefficient at rest	0.3
Poisson’s ratio	0.2 [49]
Void ratio	0.94
Density	1.39 g/cm3
Young modulus	33 MPa; 6 MPa
Rayleigh damping parameter α	0.1232 [20]
Rayleigh damping parameter β	0.6410 × 10−3 [20]

) have been assumed to be constant with depth. They have been chosen based on the mechanical behaviour exhibited by DNA-1A lunar simulant during triaxial tests [30]. In particular, the shear behaviour of a dry specimen initially confined under isotropic pressure of about 20 kPa has been of reference. The Young modulus has been assumed to be equal to two different values, i.e., 33 MPa and 6 MPa, taken at two levels of axial strain during the triaxial compression test, i.e., 0.006% and 0.2%, respectively. These stiffness values have been here calculated by further interpretation of the experimental data from Marzulli and Cafaro. It is worth noting that the first value can be measured on a triaxial specimen using local transducers for the axial strain measurement (Figure 2), directly placed on the specimen membrane, that is, using not conventional devices for the triaxial test.

The numerical simulation has been performed by means of the FEM code Plaxis 2D CE V22 (©2022 Bentley Systems, Inc. 2022). In the geometrical model, the layer has length of 48 m Laterally, viscous boundaries [50] have been adopted, to absorb the outgoing wave energy. The constitutive model adopted for the soil is the Mohr-Coulomb model. The unit weight assigned to the regolith layer in Plaxis has been calculated based on the lunar gravity, starting from density value. To assign initial stress state to the soil, the K₀-procedure has been used. The estimate of K₀ could be made by several relations proposed in the literature. In this case, the value of the at-rest coefficient (Table 2) has been deduced using a relationship by Guo [51], based on previous works by Bauer [52,53] and von Wolffersdorff [54], which expresses K₀ as a function of the Critical State friction angle, φ’_cv, and the Relative Density, D_r, of the soil. The Poisson’s ratio, ν = 0.2, has been assumed based on the value indicated by Slyuta [49] for lunar regolith (as mean value to the depth of 200 cm), consistent with the value obtained by the elastic relationship between K₀ and ν.

The numerical simulation has been carried out both without and with damping. Damping has been modelled using the Rayleigh parameters α and β (Table 2) adopted by Ricchetti et al. [20] for modelling an impact problem involving DNA-1A. They deduced these parameters using data from Senetakis et al. [55] concerning a different volcanic soil. Direct determination of these parameters on DNA-1A is desirable for future research. The impact of the meteoroid is modelled as a punctual force variable with time according to a triangular function. The impact force acts vertically.

The 2D FE domain is discretized by 2108 elements under plane strain conditions. The plane-strain condition corresponds to negligible displacement in the z-direction, therefore, along the longitudinal axis, the strain is assumed to be zero (ε_z =  0). To ensure that the waves propagating through the layer are not filtered at high frequencies, Kuhlemeyer and Lysmer [56] suggested a maximum size of the element equal to one-eighth of the wavelength associated with the highest frequency of the input signal. Furthermore, the time step should be smaller than the ratio of the distance between two nodes and the wave propagation velocity of the soil [57]. In the present study these suggestions have been considered: the average element area in the middle zone of the layer, that is the zone most influenced by the load diffusion, was about 0.04 m², and the time step was 1.5 × 10⁻⁴ s. The tolerated error for stress was 0.01 kPa.

The impact momentum adopted for this calculation comes from an impactor velocity of 18 Km/s, consistent with data from literature [13,14,58], and a meteoroid mass of 1 Kg. For simplicity, in this work a perfectly inelastic collision has been assumed. In the dynamic FEM calculation, a triangular shape of the force-time function has been assumed and, consequently, the required input parameters are the maximum value of the impact force and the time interval taken by the impact. As mentioned in the Introduction, it has been investigated, for equal impact momentum, the effect of two different hypotheses on the collision time, i.e., 0.01 s and 0.1 s, leading to maximum forces of 3600 kN and 360 kN, respectively. Finally, soil melting and vaporization [59], although expected in case of high energy impacts, are phenomena ignored in this work.

4. Results

4.1. Incidence of the Strain Interval Adopted for Stiffness Determination

The vertical stress with time at a depth of about 4.5 m along the impact vertical is shown in Figure 5, for the two values of Young modulus implemented in the FEM calculation, without (α = 0; β = 0; Figure 5a) and with (Figure 5b) Rayleigh damping. The comparison has been made for equal collision time (0.1 s) and, thus, for equal maximum impact force (360 kN). The transient effect of the impact on vertical stress is evident in all cases, but differences in terms of maximum value of the vertical stress have been found (

Table 3. The results of the FEM calculation for the three combinations of Young modulus and collision duration (or maximum force).

Investigated Hypotheses of Young Modulus and Collision Duration (or Max Force)	Max Vertical Stress [kN/m²]
	without Damping		with Damping
	Point A	Point B	Point A	Point B
33 MPa (0.1 s—360 kN)	73.0	47.3	65.8	44.8
6 MPa (0.1 s—360 kN)	123.8	98.6	102.3	74.5
33 MPa (0.01 s—3600 kN)	1043.2	788.6	632.6	435.7

). For no damping, the maximum stresses calculated in the two stiffness options, i.e., 73 kPa and 124 kPa for the higher and the lower stiffness, respectively, give evidence of an important discrepancy related to the definition of the strain interval adopted for the soil stiffness measurement. For the assumed damping (see Rayleigh parameters in Table 2), the calculated maximum stresses are 66 kPa and 102 kPa for the higher and the lower stiffness, respectively. The discrepancy calculated from the average value is about 52% in case of no damping and 43% with damping: although relevant, discrepancies arising from this kind of uncertainty could be managed when designing the regolith shield, for example by devoted safety factors, provided that the elastic moduli decay with strains is experimentally deduced with accuracy and suitable measurement devices.

In Figure 6a–d the stress evolution through the layer related to impact for the case of stiffer regolith is shown by means of stress contours, with reference to four different times. The plots give evidence of the propagation of the stress perturbation. A little stress perturbation with respect to the initial geostatic condition is detected even after the end of the impulsive force—time history (plot d).

It has been also calculated (Table 3), at the same depth of about 4.5 m, the maximum impact-induced stress along a 2 m far vertical (point B in Figure 4), which is clearly expected to be lower even assuming an elastic continuum half-space. Based on the FEM results, looking just as an example at the case of 33 MPa Young modulus and 0.1 s collision time, the calculated vertical stress at the point B decreases to 65% of the maximum stress calculated along the loading vertical (point A) in case of no damping, and to 68% when implementing damping.

4.2. Incidence of the Assumed Collision Duration

The vertical stress with time resulting from the impact, at a depth of about 4.5 m, is shown in Figure 7, for the two hypotheses of collision time and thus, as previously discussed, for the two corresponding values of maximum impact force at the ground, i.e., 360 kN and 3600 kN. The comparison is made for the case of Young modulus equal to 33 MPa, without (Figure 7a) and with (Figure 7b) Rayleigh damping. The transient effect of the impact on vertical stress is evident in all cases, but an important difference in terms of maximum value of the vertical stress has been found (Table 3). In this case, the difference is not only in terms of maximum stress but also in terms of duration of the stress perturbation (Figure 7), which is shorter for the case of 0.01 s collision duration, as expected. Varying the collision duration of one order of magnitude, from 0.1 s to 0.01 s, has caused a difference in terms of maximum vertical stress induced at a depth of 4.5 m of one order of magnitude too. At the point A, calculating the stress without damping, the maximum vertical stress is about 73 kPa for the case of 0.1 s collision duration and becomes about 1043 kPa when 0.01 s collision duration is assumed. The variation in case of calculation with damping is from about 66 kPa to about 633 kPa. Looking at the stress perturbation at the point B, the same variation in the assumed collision duration leads to corresponding stresses of approximately 47 kPa and 789 kPa without damping, and of 45 kPa and 436 kPa in case of damping (Table 3). Thus, it seems that uncertainties concerning the collision time should not be tolerated when simulating the propagation with depth of the load caused by meteoroid impact on the regolith shield hypothesized in this work. Indeed, the resulting stress uncertainty would be very important and much higher than any reasonable safety factor for the design.

5. Conclusions

The relative incidence of two sources of uncertainties on the calculation of stresses induced by meteoroid impact on a lunar regolith shield has been investigated. In particular, the effects of assumptions concerning the soil stiffness and the duration of collision have been pointed out. From the designer’s point of view, changing the assumed collision time in the calculation of meteoroid impacts should necessarily imply a change in the maximum impact force. The reason for this lies in the fact that the “design impact” should be defined on a probabilistic base, by selecting both mass and velocity of the impactor, that is, by selecting a given momentum. The uncertainty related to the stiffness of a regolith layer has been considered to arise from the different strain intervals the designer can adopt for the measurement of the elastic moduli.

Based on the FEM results, the following conclusions can be drawn:

-

For the assumed impact, the discrepancy in terms of maximum stress at the layer bottom related to the implemented regolith stiffness is 43–52% of the average value, depending on the damping assumption.

-

This discrepancy is comparable with the percentage differences in terms of maximum stress (at the bottom) between the loading vertical and a vertical 2m far: for the case of stiffness 33 MPa and force 360 kN, for example, they are 38–43%, depending on the damping assumption.

-

Changing the collision duration of one order of magnitude, for a given stiffness, leads to changes in maximum stress at the bottom of one order of magnitude.

Concluding, the results suggest that the uncertainty of the impulsive stress prediction related to the implemented soil stiffness, although not negligible, could be managed by suitable safety factors to avoid damage to the buried structure, whereas the collision duration to be implemented in the calculation should be defined with very great accuracy. In this respect, the combined use of numerical simulation and physical modelling of impact processes is probably the strategy to be pursued.