Discrete Modeling of Projectile Penetration into Calcareous Sand Subgrade Runway Pavement Structure and Genetic Algorithm Optimization for Critical Striking Velocity

1. Introduction

The airport has been the hot-target of military attack owing to its large footprint, obvious targets and prominent military status, which can effectively neutralize the enemy’s capability to conduct air operations [1]. Mainly composed of a surface layer, a base layer and a compacted soil layer, the airport runway is a typical layered structure [2]. Although some studies proved that asphalt concrete surfaces could absorb the impact energy [3], the concrete pavements is widely used in the airports due to the fact that the high temperature and aircraft load generated by takeoffs of jet airplanes can damage and deteriorate asphalt surfaces [4].

Calcareous sand is a sediment primarily composed of the skeletal remains of marine organisms, with calcium carbonate accounting for over 97% of its weight [5,6]. It is widely found in the coral reefs of the Nansha Islands in the South China Sea and is characterized by high intragranular porosity, significant particle angularity, irregular shape, and a tendency to break easily [7,8,9]. Moreover, calcareous sand exhibits higher apparent cohesion than quartz sand and demonstrates considerable load-bearing capacity after compaction [10], making it an ideal backfilling material for the subgrade of airport runways in marine engineering. Additionally, due to the remote location of the South China Sea islands, using calcareous sand as a foundation material for island construction projects effectively reduces transportation costs and minimizes the reliance on expensive imported materials [11]. It is believed that the island airport with calcareous sand subgrade runway pavement structure is vulnerable to cluster munitions, whereby the projectile penetrate to a certain depth with detonation as shown in Figure 1. The firt stage damage composed of penetration is critical for the second stage internal explosion, which needs to be comprehensively studied for damage evalutation and munition optimization.

The destructive effect and impact resistance performance of airport runways structured with multi-layer pavements under impacting loads have been a significant area of research. Han et al. conducted the field penetration and blast tests on airport pavements subjected to blast loading from a cluster bomb unit and assessed the condition of damage by numerical simulation, suggesting the rapid hardening high-strength concrete for the damaged concrete pavements repair [12]. Wu et al. developed a new muti-layer pavement system consisting of asphalt concrete layer reinforced with geogrid, high strength concrete layer and engineered cementitious composites layer, and the system was compared with conventional concrete pavements through field blast tests and numerical modeling, revealing that the new multi-layer pavement had better blast resistance [3]. Zha et al. used the finite element simulation to compare the damage contours and dynamic responses of the anti-blasting asphalt pavement with a reinforced concrete layer and the conventional asphalt pavement, which showed that the additional reinforced concrete layer improved the blast resistance performance by about 20% [13]. To investigate the damage of scaled-down target of the airport runway, Wei et al. carried out experiment and numerical simulation to analyze the destructive effect of static explosion of charge after penetration and explosion of prefabricated hole, which showed that the crater size produced by the former explosion mode was obviously larger than that produced by the latter under the same explosion depth and explosive charge [14]. A three-dimensional graphics processing unit (GPU)-accelerated smoothed particle hydrodynamics (SPH) method was developed by Chen et al. [15,16], and was applied to the simulation of soil fragmentation and fracture propagation of the concrete-soil multilayered medium subjected to underground explosion involving millions of particles. Comparing with the experimental data, the SPH model was able to reproduce the damage pattern of concrete slabs with different depths of explosives burial. Although the dynamic impact responses of airport runway with multi-layer pavements has been extensively studied through experiments and numerical simulation, there is still a knowledge gap in the research on pavement structures with calcareous sand as the subgrade, which hinders the development of island infrastructure engineering.

As the pavement structure surface layer, the concrete panel or plate under the impactive loading conditions has been extensively studied [17,18,19,20]. The experiments and numerical simulations of the structural behavior of full-scale reinforced concrete slabs under blast loading were conducted by Castedo et al. [21], finding that the addition of steel fibers or polypropylene fibers to concrete slabs could effectively improve the tensile strength and blast resistance compared to simple reinforced concrete slabs. Feng et al. introduced the impact resistance theory to establish a semi-empirical analytical model for projectile perforation on steel-concrete-steel sandwich panels and an analytical model for perforation on steel-concrete panels, respectively, to predict the residual velocities and the target damage patterns, analyzing the dynamic response of the projectile on concrete-steel panels of different thicknesses [22,23]. The fiber reinforced concrete runway pavement proposed by Ali et al. [24] was proved much better than the normal concrete pavement, and the researchers came up with finite element model to analyze the performance of improved concrete runway pavement under impact loading based on the results of drop weight impact test. Hardened steel ovoid nose projectiles with 19 mm diameter and 3 caliber-radius-head of 0.5 and 1 kg mass and 11.8 and 23.7 length-to-diameter ratios were used to impact plain and reinforced concrete targets in conjunction with numerical simulations in order to find out the behavior of the concrete against the impact loading as well as to explore the length-to-diameter ratio of the projectile effect on the ballistic performance [25].

Discrete Element Method (DEM) has been widely adopted to quantitatively analyze the deformation and damage processes occurring in concrete or geotechnical materials under both macro and micro mechanical regimes [26,27]. The unique advantage of DEM is that it is independent of numerical mesh and allows large deformation and fracture of particles from the microscopic layer without deformation constraints, which is suitable for dealing with discontinuous problems [28,29]. Geng et al. applied DEM to analyze the infiltration distribution changes of tailings slopes and tailings landslides in different rainy seasons by constructing discrete units with specific elastic modulus and strength properties with the combination of the conversion formulas for the discrete units and numerical simulation experiments, realizing the automatic modeling of the discrete units [30]. Mechtcherine and Shyshko presented a numerical approach based on DEM to establish a link between the yield stress of the simulated concrete and the model parameters as a parametric reference for modeling the behavior of fresh concrete in different working processes, which showed that the numerical analysis agreed well with the final shape of the concrete in the slump-flow test [31]. DEM can also be used to study the physical process of impact-induced rock fragmentation in rockfall analysis [32]. A normal impact fragmentation model of synthetic spherical rock block under different impact loading rates was carried out by Shen et al. [33], which could reproduce the whole process of rock fragmentation under the impact, being in good agreement with the experimental observations. Geng et al. used matDEM software to develop an asteroid impact sampling model under microgravity conditions to explore the effects of the initial impact velocity and projectile shape on the internal characteristics of regolith and ejecta after impact, verifying the feasibility of discrete element simulation of the impact process [34]. Lattice Discrete Particle Model (LDPM) was proposed to simulate concrete at the aggregate level and characterize the formation and evolution of cracks in concrete under loading, investigating the structural size and geometry effect on the strength and fracture process of concrete, which was proved to have a good agreement with experimental data and performed well in the prediction of splitting test results [35]. Therefore, it is advisable to apply DEM to model and investigate the dynamic responses of the airport runways with concrete and calcareous sand layers with matDEM as an effective tool.

Over the past decades, extensive research has been conducted on the damage caused by various weapons to airport runways [36]. However, nowadays, the use of cluster munitions to create multiple craters on the runways, achieving the blockade of airports, has become the primary attacking strategy [37]. Based on the analysis of the crater characterization and the mechanism of implosion in the runways obtained from the experiments and simulations, the damage modes of implosion in airfield runways can be classified into three categories: open crater mode, bulge crater mode and hidden crater mode [38], as shown in Figure 1. In particular, for anti-runway munitions, the bulge crater, cuased by explosion near interface, is the severest damage of the above three destruction modes, creating a relatively long blockade of the runway and being the most difficult and time-consuming to repair. Hence, the critical striking velocity corresponding to the interface location DOP is deemed to be the optimal penetration scenario.

In order to assess the damage effects of island airport runway pavements under projectile penetration, this paper adopted DEM simulation to numerically model the penetration of rigid projectiles into runway pavements with calcareous sand subgrade. The material parameters were calibarated with uniaxial compression and triaxial compression tests for sequent penetration model validation. Furthermore, the penetration depth of airport runway structures with concrete pavement and calcareous sand base were numerically evaluated. Combining with back-propagation (BP) neural network and genetic algorithm, the critical velocity of airport runway penetration was predicted via intelligent optimization of pavement penetration. The results may shed some light on damage assessment as well as anti-airfield runway cluster munitions design.

2. Discrete Modeling of Granular Material

2.1. Discrete Element Method

In this study, a DEM model of granular material was established using MatDEM, achieving by stacking and cementing a large number of particles that conformed to Newtonian equations of motion, while the contact and interaction between the particles were simulated by a breakable elastic spring that acted only on the contact points between neighboring particles [26,39]. The forces between particles can be categorized into normal and tangential forces.

Figure 2b depicts the connection of two particles in the normal direction, including the normal force and normal deformation, which can be simulated by a normal spring. When the neighboring particles are connected to each other, the tensile or compressive spring force acting on the particle can be expressed as follows:

$$F_{\mathrm{n}}=K_{\mathrm{n}}·X_{\mathrm{n}},X_{\mathrm{n}}<X_{\mathrm{b}}$$

(1)

where $K_{\mathrm{n}}$ is the normal stifness; $X_{\mathrm{n}}$ is the normal displacement relative to the equilibrium position.

When the normal relative displacement $X_{\mathrm{n}}$ between the two particles is greater than the breaking displacement $X_{\mathrm{b}}$, the normal spring between the particles breaks and the tension force becomes non-existent, at which time the normal spring force between the particles is zero. Under intact bonding conditions, the maximum normal force between the particles $F_{\mathrm{nmax}}$ can be calculated as Equation (2).

$$F_{\mathrm{nmax}}=K_{\mathrm{n}}·X_{\mathrm{b}}$$

(2)

where $X_{\mathrm{b}}$ is the fracutre displacement.

When the state of the two particles is in compressive contact, repulsive force is generated between them. The normal spring force can be expressed by Equation (3) and is negative.

$$F_{\mathrm{n}}=K_{\mathrm{n}}·X_{\mathrm{n}},X_{\mathrm{n}}<0$$

(3)

Similarly, the shear force and shear deformation between particles can be simulated by a tangential spring [40], as shown in Figure 2c, and the expression for the shear force ($F_{\mathrm{s}}$) is as follows:

$$F_{\mathrm{s}}=K_{\mathrm{s}}·X_{\mathrm{s}}$$

(4)

where $K_{\mathrm{s}}$ is the shear stifness; $X_{\mathrm{s}}$ is the tangential relative displacement.

The destruction of tangential springs follows the Moho-Coulomb criterion, which is able to determine the maximum shear force of the intact bond, as expressed by Equation (5).

$$F_{\mathrm{smax}}=F_{\mathrm{s}0}-{\mu }_{\mathrm{p}}·F_{\mathrm{n}}$$

(5)

where $F_{\mathrm{smax}}$ is the maximum shear force; $F_{\mathit{s}\mathit{0}}$ is the inter-particle shear force; ${\mu }_{\mathrm{p}}$ is the inter-particle coeffcient of friction.

When $F_{\mathrm{n}}$ is zero, $F_{\mathrm{s}0}$ reaches the maximum, and the shear force increases with the increasing of $F_{\mathrm{n}}$. When the external force exceeds the maximum shear, the tangential connection breaks and only sliding friction exists between the particles as follows Equation (6).

$$F_{\mathrm{smax}’}=-{\mu }_{\mathrm{p}}·F_{\mathrm{n}}$$

(6)

2.2. Calibration for Concrete and Calcareous Sand

Mesoscopic mechanical parameters of DEM, including normal stiffness ($K_n$), shear stifness ($K_s$), breaking displacement ($X_b$), initial shear resistance ($F_{s0}$) and friction coefcient (${\mu }_p$), can be roughly estimated from five basic macroscopic mechanical properties theoretically, including young’s modulus (E), Poisson’s ratio ($\nu$), compressive strength ($f_c$), tensile strength ($f_t$) and internal friction coefcient ($\mu$). The mechanical properties of elements are determined by the conversion formulas. In the regular packing model, there are analytical solutions between the micro mechanical parameters of the elements and the macro mechanical properties of the model, i.e. the conversion formulas:

$$K_{\mathrm{n}}={\displaystyle \frac{\sqrt{2}Ed}{4(1-2\nu )}}$$

(7)

$$K_{\mathrm{s}}=K_{\mathrm{n}}{\displaystyle \frac{(1-5\nu )}{(1+\nu )}}$$

(8)

$$X_{\mathrm{b}}={\displaystyle \frac{3K_n+K_s}{6\sqrt{2}{K}_n(K_n+K_s)}}{f}_t{d}^2$$

(9)

$$F_{\mathrm{s}0}={\displaystyle \frac{1-\sqrt{2}\mu }{6}}{f}_c{d}^2$$

(10)

$${\mu }_{\mathrm{p}}={\displaystyle \frac{-2\sqrt{2}+\sqrt{2}I}{2+2I}}$$

(11)

where I is the frictial parater defined as $I={[{(1+{\mu }^2)}^{1/2}+\mu ]}^2$.

In this section, uniaxial and triaxial compression tests were conducted on concrete and calcareous sand to calibrate the discrete element method parameters, thereby enhancing the credibility of the research on penetration resistance. The uniaxial compressive strength is a crucial indicator, reflecting the strength limit and stress-strain relationship under loading conditions, and is essential for analyzing the bearing capacity and deformation characteristics of concrete structures.

According to Table 1, the density, Young’s modulus, Poisson’s ratio, tensile strength, and compressive strength of concrete and calcareous sand were used as inputs for generating DEM particles. The concrete primarily consisted of coarse aggregate with a particle size of approximately 10 mm and fine aggregate with a particle size of about 1 mm. In contrast, calcareous sand had an irregular particle size, with an average ranging from 0.5 mm to 0.9 mm. The DEM parameters for the concrete shown in Table 2, under the condition of 23 MPa for a cubic specimen, are as follows: normal stiffness $K_n=5.42\times {10}^9$ N/m, shear stiffness $K_s=5.38\times {10}^8$ N/m, fracture displacement $X_b=3.76\times {10}^{-7}$ m, initial shear resistance $F_{s0}=6.31\times {10}^4$ N and coefficient of friction ${\mu }_p$ = 0.19. For the calcareous sand sample, the required DEM parameters are: normal stiffness $K_n=1.01\times {10}^9$ N/m, shear stiffness $K_s=1.14\times {10}^8$ N/m, fracture displacement $X_b=3.02\times {10}^{-5}$ m, initial shear resistance $F_{s0}=2.87\times {10}^4$ N and coefficient of friction ${\mu }_p$ = 0.14.

Based on the experimental data obtained by Warren et al. [41], numerical simulation of uniaxial compression was carried out for 15 cm × 15 cm × 15 cm cubic concrete specimens. Fixed constraint was conducted at the bottom of the specimen and the load was applied at the top at a rate of 2.4 N/ms [42]. The loading condition and the final vertical fracture damage of the model are illustrated in Figure 3. The simulated stress-strain curve behaved consistently with the curve obtained from the uniaxial compression test, as shown in Figure 4a. Initially, the concrete worked in the elastic stage and the stress increased proportionally with the strain. With the increase of compressive stress, micro-cracks inside the concrete gradually developed and plastic strain began to occur, the stress-strain curve no longer maintained a linear relationship and the slope gradually decreased. Upon reaching the peak stress, visible cracks began to appear on the surface of the specimen, and under continued loading, these cracks rapidly penetrated from the surface through the interior of the material, leading to complete destruction of the specimen.

To assess the mechanical properties and strength characteristics of concrete and calcareous sand under more complex stress states, triaxial compressive responses were conducted using axial displacement control in numerical simulations while maintaining constant lateral confining pressure. Figure 4b displays the axial stress versus axial strain curves for concrete under limiting pressures of 50 MPa, 100 MPa, and 200 MPa. Figure 5 shows the corresponding curves for calcareous sand at limiting pressures of 200 kPa, 400 kPa, 800 kPa, and 1200 kPa. The stress-strain curves of calcareous sands exhibited nonlinearities, which were mainly caused by particle crushing and the tansformation of the relative positions of sand particles during the loading process [43]. The curves of both calcareous sand and concrete had the similar variabilities, which showed that the longitudinal strength and deformability of the specimens of the two materials improved with the increase of transverse pressure. The peak values of the stress-strain curves rose, and the decline trend after the maximum stress value gradually slowed down with the increase of the ultimate pressure value, while the peak stress, peak strain and elastic modulus also increased significantly. Additionally, the increase in lateral pressure restricted the transverse deformation of the specimens, thereby limiting both lateral expansion and the propagation of internal microcracks. Strain hardening arose from the heterogeneity of the granular skeleton of concrete and calcareous sand materials [44], and the axial response exhibited a transition from strain softening to strain hardening as the microstructure within the materials changed with the increasing confining pressure [45].

2.3. Validation of Penetration Model

The validation of the discrete element method in simulating uniaxial and triaxial compression tests for concrete and calcareous sand has been validated in the previous section, demonstrating its potential for more complex loading simulations of these materials. In this section, the numerical performance of DEM for the penetration modeling of concrete and calcareous sand were validated by both the depth of penetration and projectile deceleration.

As dimensioned in Figure 6, the ogival nose projectile for concrete penetration with 76.2 mm caliber, 3.0 CRH and a total length of 530.73 mm was machined from 4340 RC45 steel, and the nominal mass of the projectile was 13 kg. The conical nose projectile for calcareous sand was internally hollow, with a diameter of 14.5 mm and a length of 87 mm, which was made of high-strength alloy steel with a density of 7850 $kg/m^3$, as mentioned by Miao et al. [46]. The body of projectile was set to be a rigid material due to the negligible wear and tear during the penetration process.

The concrete penetration target was cylindrical with a diameter of 1.83 m, which was large enough to avoid rear and lateral boundary effect, and meanwhile, the target was laterally constrained. The calcareous sand penetration target had a cross-sectional size of 0.6 × 0.6 m and a thickness of 1.7 m, with fixed lateral boundries. Penalty contact was adopted to compute the normal contact force between the granular material particles and the rigid projectile [47]. The experimental data for penetration into concrete and calcareous sand, as modeled in this paper, were derived from Forrestal [48] and Miao [46], respectively. The numerical results of DOP values for both concrete and calcareous sand under impact tests are presented in Table 3 and Table 4, showing that penetration into concrete or calcareous sand was validated by examining the DOP, with a maximum relative error not greater than 5.6 %.

The deceleration of simulations and tests for projectiles penetrating concrete targets at three different striking velocities of 238 m/s, 276 m/s, and 370 m/s were comparatively presented in Figure 7 with G denoting the gravitational acceleration. In terms of general trends, the DEM simulation results matched well with the experimental deceleration curves. The deceleration histories for all striking velocities exhibited a rapid increase during the initial projectile nose penetration phase followed by a nearly flat plateau corresponding to the tunneling phase, where the deceleration values remained relatively steady but fluctuated slightly due to the penalty function contact model between the concrete particles and the projectile. Additionally, higher impact velocities extended the duration of the tunneling phase and resulted in higher deceleration peak values. At last, the projectile suffered with increasing resistance inside the target, causing a sharp drop in deceleration to zero after the plateau phase. The cross sections views of the simulation of projectile penetration into a concrete target at different moments are illustrated in Figure 8, where the depth of penetration gradually increased with time. The concrete particles were pushed away by the projectile nose, which were squeezed to cause concrete scattering fragments around the crater.

Obtained from the numerical simulations, the projectile deceleration histories during calcareous sand penetration are plotted in Figure 9. In the early stage of the projectile penetration, the deceleration increased rapidly in a short period of time until the maximum value, which was affected by the velocity. The higher striking velocities resulted in higher peak deceleration, which was about 32% higher at 710 m/s compared with that at 411 m/s. After the peak point, the deceleration of the projectile declined rapidly and subsequently entered a platform phase. At this time, the target static resistance played a dominant role in the penetration resistance, decreasing the deceleration slowly [49].

3. Pavement Structure Penetration Simulation

Given the specific characteristics of the airport runway, the penetration process is divided into two distinct phases: the concrete penetration phase and the calcareous sand penetration phase. In this section, the previously validated DEM model was employed to conduct numerical simulations of runway pavement penetration with various striking velocities and concrete thicknesses, with the aim of providing a comprehensive assessment of the airport runway structure under impact.

The model of projectiel penetration into airport runway pavement structure is depicted as Figure 10. In this paper, the Russian-made BETAR-25 anti-runway cluster munition was selected for numerical simulation, with ogival nose projectile body with a total weight of 9-13 kg, a diameter of 76 mm, 3.0 CRH, and an overall length of 480 mm [50]. With stricking velocity of 170 m/s-280 m/s, the projectile is supposed to impact the target with normal penetration.

The island airport runway model was constructed as a two-layer composite structure with a concrete layer $h_c$ for the surface and a calcareous sand layer $h_s$ for the base, measuring 1.8 m × 1.8 m area, with lateral constraints. The model particles were randomly generated using Fuller’s coefficient 0.5 in the range of aggregate particle size of 5-15 mm.

Table 5 shows the simulation results of DOP of various velocities of projectile and thickness of concrete layer. Among all the results, three representitive striking velocities (200 m/s, 230 m/s and 260 m/s) as well as three different concrete thickness (20 cm, 30 cm, 40 cm) have been selected and plotted in Figure 11 with 3 different colors. During the surface penetration phase, the velocity of the projectile into the concrete decreased while the deceleration initially increased. As the resistance reached its maximum value, the deceleration also peaked and was maintained in a short plateau period with slight fluctuations. However, as the projectile penetrated deeper, it was subjected to the resistance, leading to the decline of velocities and decelerations. After entering the calcareous sand layer, which offered significantly lower resistance than concrete, the reduction in velocities became slow down. In addition, when the projectile navigated through the pit-opening stage in the calcareous sand layer, the deceleration rose slightly, and after the pit-opening stage, the resistance acting on the projectile head became smaller due to the degradation of the tensile stress between the surface of the projectile body and the calcareous sand particles [45]. On the other hand, the lower striking velocities further reduced the resistant force. These two main factors both contribute to the decrease of the deceleration at a low rate and the stabilization with slight fluctuations.

Since the impact velocity and mass of the projectile were constant, the peak of the deceleration curves remained almost unchanged regardless of the target thickness, although thicker targets led to longer penetration durations. However, the projectile was unable to penetrate the concrete layer at a thickness of 40 cm. Variations in striking velocities resulted in changes to peak deceleration, with higher velocities producing greater peak deceleration values and shorter penetration durations.

4. Intelligent Optimization of Pavement Penetration

4.1. GABP Neural Network Proxy Model

The back propagation neural network (BPNN) is a multi-layer feedforward network that employs the error backpropagation algorithm [51], which is considered one of the most classic and extensively utilized models in neural networks [52,53]. BPNN excels in handling non-linear relationships and in the training process of continuous iteration, the primary concept involves computing the error between the output data and the expected data, followed by adjusting the weights and thresholds in a retrograde fashion to minimize the error [54,55]. The conventional BP neural network algorithm relies on the gradient descent technique for optimization. However, the random initialization of weights and thresholds, coupled with the issue of wide variation in initial value selection, may lead the network to get stuck in local optima [56]. To improve the optimization efficiency of pavement penetration, this work adopted a novel proxy model approach, combining genetic algorithms with (GA) neural network models, improving the forecasting capabilities of the overall model. The GABP proxy model flowchart is shown in Figure 12, implying the following detailed steps:

Step 1 Data preprocessing and network structure determination. Divide the original data into training, testing, and validation sets, adopt the single hidden layer BPNN model, and configuring genetic algorithm parameters.

Step 2 Population initialization and interations. Encode the initial value according to the fitness function and generate new individuals by carrying out selection, crossover, and mutation operations.

Step3 New population updation. Recalculate the fitness value for the updated group, update the global optimal individual meeting the fitness value conditions until reaching the maximum number of iterations.

Step 4 BP neural network training. Assign the optimal individual for BPNN, obtain the optimal weights and thresholds, perform training simulations on error evaluation calculations and backpropagation until convergence.

4.1.1. Data Processing

The airport runway penetration scenarios were obtained from a series of simulations, encompassing projectile weight, pavement concrete thickness, striking velocity, and penetration depth as the features. A total of 52 samples were gathered which shown in Table 5, with 80% of the data randomly assigned to the training set and the rest allocated to the testing set. To convert input and output data to values within the range of [0,1], the transformation formula is utilized as follows:

$${\overline{R}}_i={\displaystyle \frac{R_i-R_{min}}{R_{max}-R_{min}}}$$

(12)

where $R_{\mathrm{i}}$, $R_{\min }$, $R_{\max }$ represents the input data, the minimum value and the maximum of the data, respectively.

4.1.2. Development of GABP Proxy Model

The BP neural network structure is shown in Figure 13. The input layer comprised three nodes corresponding to three eigenvectors: projectile weight, pavement concrete thickness and striking velocity. The output layer consisted of one node representing penetration depth. The number of nodes in the hidden layer was determined by the neurons in the input and output layers, with five neurons chosen for the hidden layer in this study.

The initial parameters were configured as 1000 times training iterations with a target of 0.005 and a learning rate set to 0.1 using the TRAINLM training function. In the genetic algorithm, 500 evolutions were performed with a population size of 200. The selection factor was set to 0.09, while the crossover factor was randomly chosen within the range of 0 to 1. A non-uniform mutation approach was applied, utilizing a non-uniform probability distribution to adjust parent parameters. This implies that the mutation probability evolved over time, potentially decreasing as iterations increased or as the population included a greater number of fitter individuals. The initial parameters for the GABP proxy model are presented in Table 6.

4.1.3. Assessment of GABP Proxy Model

In order to evaluate the accuracy of the training and prediction performance, two statistical parameters, which are the root mean square error ($RMSE$) and determination coefficient ($R^2$) are served as the error discrimination criteria[57]. The root mean square error ($RMSE$) indicates the predictive capability of the model while determination coefficientassesses the fit of the model to the data. The definitions for $RMSE$ and $R^2$ are shown as follows:

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left(P_i-R_i\right)}^2}$$

(13)

$$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(P_i-R_i)}^2}{{\sum }_{i=1}^n{(\overline{R}-R_i)}^2}}$$

(14)

where n is the number of samples, $P_{\mathrm{i}}$ is the predicted value, $R_{\mathrm{i}}$ is the real experimental value and $\overline{R}$ is the average of the real values. The smaller $RMSE$ value indicates better predictive ability of the model;similarly, the value of $R^2$ closer to 1 indicates a better fitting effect of the model.

Figure 14 displays the scatter plot of predicted versus actual values for both networks, which are closely related to each other. The GABP neural network proxy model performance peaked at epoch 10, with a regression correlation coefficient exceeding 0.99 across all sets. The average error of the training and testing sets shown in Table 7 are 0.0886 m and 0.0971 m. The subsequent prediction results demonstrated that the BP neural network optimized by the GA algorithm not only significantly validated the accuracy of the neural network model but also exhibited superior performance in predicting the airport runway penetration depth.

4.2. Genetic Algorithm Search for Critical Striking Velocity

Based on crater morphology and crack extension, the damage patterns of airport runways can be classified into three categories: open crater, bulge, and hidden crater patterns [58], as illustrated in Figure 1. From the perspective of post-damage repair, the open crater and hidden crater damage modes are generally considered easier to repair. In wartime emergencies, these types of damage can often be addressed by simply clearing and filling the damaged area, followed by the placement of a steel plate to allow temporary aircraft operations. However, when rumble damage occurs in the surface layer concrete, the presence of ring cracks penetrating the surface layer significantly increases the volume of material to be removed, making repairs much more challenging. Consequently, genetic algorithms have been used to predict the critical velocity at which concrete layer is just perforated, in order to maximize explosion damage to the pavement structure.

Genetic algorithm is a computational method that randomly searches for optimal solutions by simulating the natural evolutionary process [59,60], which draws on natural selection and genetic mechanisms in biology. The genetic algorithm mainly includes population initialization, fitness function evaluation and genetic operator determination. The fitness evaluation function reflects the ability of individuals to adapt to the environment. Genetic operators mainly include selection, crossover and mutation [61]. Through the iterative repetition of these steps, the algorithm gradually searches for the optimal solution.

The initial population size was set to 100, with a maximum of 100 iterations and the crossover probability and mutation probability were set at 0.8 and 0.2, respectively. The optimal penetration velocity at the point of maximum damage was determined using the condition of just penetrating the concrete surface as the constraint boundary. Additionally, the larger adaptation degree represented the optimal penetration speed corresponding to it.

As can be seen from Figure 15, we used the fitness function as a criterion for the loss value. Observing the Loss value, the concrete layer of 20 cm had converged in the 22nd iteration, the concrete layer of 30 cm had converged in the 8th iteration and the concrete layer of 40 cm had converged in the 3rd iteration. The fitness function proved that the GA neural network proxy model provided a good fit and prediction of the actual values and the model had a strong explanatory power. When the thickness of the concrete surface layer was 20 cm, 25 cm, 30 cm, 35 cm and 40 cm, the optimal speed was 282.19 m/s, 308.79 m/s, 328.68 m/s, 345.34 m/s, 378.39 m/s, respectively which shown in Table 8.

5. Concluding Remarks

It is of great interest to assess the damage of the island airport pavement with multi-layer structures under impact loading. In this study, DEM was employed to model rigid projectile penetration through runway pavements composed of concrete and calcareous sand, with the model calibrated using uniaxial and triaxial compression tests. The critical impact velocity causing maximum explosion damage was determined using GA optimization via GABP proxy model. The main conclusions are drawn as follows:

(1) Based on analysis of normal and shear interactions between adjacent particles by DEM, the numerical model for hard projectile penetration on concrete and calcareous sand targets could be validated against the test data.

(2) Deceleration history of projectile penetration in concrete pavement and calcareous sand subgrade multi-layer structure suggested a pulse due to concrete layer peforation resistance which is followed by a long tail resulting from low strength calcareous sand medium static resistance.

(3) The BP neural network optimized by the GA algorithm was capable of predicting the depth of penetration of airport runways pavement structures with high accuracy.

(4) Calling the GABP proxy model, the GA optimization was performed for most destructive effect of concrete surface layer with thickness of 20 cm, 25 cm, 30 cm, 35 cm and 40 cm whereby the critical velocities were predicted as 282.19 m/s, 308.79 m/s, 328.68 m/s, 345.34 m/s, 378.39 m/s, respectively.

[custom]