A Semi-Analytical Solution for Shock Wave Pressure and Radius of Soil Plastic Zone Induced by Lightning Strikes

A semi-analytical solution for forecasting the soil behavior induced by lightning strikes is of great engineering significance to calculate the radius of the soil plastic zone. In this paper, a simplified two-stage method is employed to solve the shock wave pressure and the radius of the soil plastic zone. The solution is verified against experimental data. Using the present model, the major factors dominating the shock wave pressure and the radius of the soil plastic zone are investigated. The results show that (1) the radius of the soil plastic zone (rp) induced by lightning decreases monotonically with cohesion (c) and internal friction angle (φ), while c has a better effect on soil properties than φ does; (2) increasing the initial radius of the plasma channel (ri0) can reduce the pressure (P) and increasing ri0 has a nonnegligible effect on rp; with ri0 increasing by 100%, the radius of the soil plastic zone increases by 47.9–59.7%; (3) the plasma channel length (L) has a significant influence on P and rp, especially when L is at a relatively low level; (4) the rp induced by lightning decreases exponentially with attenuation coefficient (a); (5) the wavefront time is a major factor while the half-value time is a minor factor for the shock wave pressure induced by plasma explosives.


Introduction
Lightning is a dangerous, frequent geophysical phenomenon [1,2], and large amounts of heat and electricity are produced by considerable energy and intense electromagnetic field activities [3,4]. The lightning forms a plasma channel that propagates outward from its starting point with positive and negative charges at each end, known as positive and negative leaders, respectively [3]. For example, lightning caused a 50 m long and 1 m wide breach in the ground in Yuncheng, Shanxi Province, China and caused a pothole in Batok Town, Bukit, Singapore (see Figures 1 and 2, respectively) [5,6]. Lightning struck buried pipelines in Alabama and Cideville, France, both causing pipeline leaks and posing a significant threat to pipeline transportation [7,8]. Lightning struck an airport runway on the west coast of India, causing damage to the runway as well as cracks to the roadbed [9]. In the event of a lightning strike, the soil serves as a significant dissipation medium and its dispersion qualities are critical [10][11][12]. Thus, it is necessary to investigate the principle of soil plastic deformation induced by lightning strikes.
Thus, it is necessary to solve the shock wave pressure induced by lightning. Meanwhile, lightning is a special form of electrical impulse. Some studies in the literature provided analytical solutions to solve the process of plasma explosives. Drabkina et al. [13] carried out innovative work in the arc channel expansion of pulsed gas discharges; however, their model cannot calculate the conductivity and temperature change of the arc channel. Braginskii et al. [14] proposed an arc channel energy equation to study the arc channel expansion, but they assumed that the current is proportional to time, so the calculation can only be used in the arc channel expansion initial stage. Based on the work of Drabkina et al. [13] and Braginskii et al. [14], Engel et al. [15] proposed a method for calculating the channel radius with the arc current as the independent variable, and gave a pressure-time relationship. From the arc channel collision expression of Engel et al. [15], which shows that the arc channel radius increases with time, it is not in line with the actual development of the arc channel law. Numerous analytical solutions have been provided to forecast water shock waves induced by arc discharges [16][17][18][19][20][21]. A physical and mathematical model of electro bursts has been provided to forecast the shock wave pressure generated by the expanding plasma channel [22,23]. Burkin et al. [24] proposed a model for electric pulse rock breaking; however, the model assumes that the plasma channel is a cylinder of constant length, which needs to be improved in subsequent research.  The studies mentioned above have provided analytical solutions of the shock wave pressure induced by arc discharges in air, water, and solid media. However, the equation for the pressure in soil has not been given. A physical and mathematical model of rock breaking has been provided to study the mechanism of high-voltage pulsed rock fragmentation [25][26][27]. Li et al. [28] investigated the electro-dynamic mechanisms of electric breakdown by using the time-transient breakdown model and the finite-difference numerical method. Peng et al. [29] used experiments and numerical simulations to investigate the causes of damage to rocks caused by arc blasts and concluded that the surrounding pressure is the main factor affecting rock cracking, but they did not give a quantitative evaluation index. Pulse discharge technology is an innovative technique that has been used to enhance the bearing capacity of piles, resisting capacity of anchors, and rock-breaking drilling [30][31][32]. Walsh et al. [33] simulated the passage of the pulse, electrical breakdown in the rock, and the mechanical response at the grain-scale by building a multiphysics model, but this model does not implement the mutual coupling between the electric-thermal force. Chen et al. [34] investigated the thermal and impact effects of lightning on rock, using Joule heat theory to obtain the surface temperature distribution of granite, but they did not consider the effect of the current waveform on shockwave pressure. Yan et al. [35] investigated the breakdown process and fragmentation characteristics of anthracite through high-voltage electrical pulse treatment, but they did not explain the phenomenon well.
The purpose of this paper is to provide a semi-analytical solution based on the work of Rao et al. [36,37] for solving the shock wave pressure and the radius of the soil plastic zone (RSPZ) caused by lightning strikes. The proposed solution is verified against the experimental data. Considering the equation of state of the soil, the Runge-Kutta method is applied to solve the governing equation for shock wave pressure. Then, the modified Mohr-Coulomb yield criterion is applied to solve the RSPZ. The major factors dominating the shock wave pressure and the RSPZ are investigated. Soil is a kind of mineral material, and this paper can provide a further understanding of the mechanical properties of soil after lightning strikes and can provide some guidance for lightning mitigation. Figure 3 shows the geometric deformation model of soil by lightning impact. The system is composed of four parts, i.e., lightning, plasma channel, shock wave, and soil. The principle of the plasma channel in soil is introduced by the theory of electrical breakdown of solid dielectrics. When the plasma channel is formed, the energy in the pulse power source is released into the plasma channel (up to 10 4 K) and the channel is heated. Then, the plasma channel is rapidly expanded, generates shock waves to act on the surrounding soil, and finally leads to the deformation of soil [25]. Several suppositions are applied to assist in the solving of the analytical solution: (1) a standard double index current of lightning is assumed to be fixed at expression I(t) [38];

Geometric Model and Simplifying Assumption
(2) the property of soil is assumed to isotropic; (3) the path of the plasma channel is assumed to be a line, i.e., not including branches of a channel.

Lightning Energy Diffusion Model
When a lightning pulse discharge occurs, the energy in the pulse power source is released into the plasma channel and converted into other forms energy. The schematic of energy conversion is shown in Figure 4. The total energy of (W c ) is mainly consumed by resistance loss in the electrical circuit (W R ) and the energy injected into the plasma channel W ch . The energy injected into the plasma channel (W ch ) is converted to the plasma channel heat energy (W cu ) and the shock wave energy (W cw ) to work with channel expansion, wherein the shockwave energy W cw consists of internal energy (W ce ), kinetic energy (W cd ), and reflected wave energy (W cl ) [22]. The channel energy balance equation for W ch , W cu , and W cw is a key equation to connect plasma channel energy and the discharge energy converted into plasma channel heat energy and mechanical work. The relationship between W ch , W cu , and W cw has been described by Burkin et al. [22].
where dW cw = P ch ·dV ch describes the channel work by the expanding channel at its volume V ch = πr 2 ch (t)l ch with pressure P ch , where r ch (t) is the channel radius. W cu = P ch V ch (γ − 1) −1 describes the plasma channel heat energy, where γ is the effective ratio of specific heats.
The governing equation of the energy injected into the plasma channel W ch is [23] where i(t) is the lightning pulse current and R ch (t) is the breakdown channel resistance. The breakdown channel resistance adopts the Weizel-Rompe model of impedance [39]. The equation of impedance of the Weizel-Rompe model is [22,23,39].
where K ch is a spark constant and l ch is the length of the plasma channel.
The double exponential current expressions i(t) is where i m is the peak current, k m is the correction coefficient, α is the wavefront coefficient, and β is the half-value coefficient of the double exponential current.
Substituting the double exponential current expressions i(t) into Equation (3): When solving the problems in the circuits comprising only the spark channel as a load, Equation (5) gives R ch = ∞ at t = 0, which makes the initial conditions uncertain. In this paper, the approximation is where ∆ ≈ πK ch r 0 σ 0 , r 0 is the initial radius of the plasma channel, and σ 0 is the electrical conductivity averaged over the cross-section of the plasma channel.

Stress of Shock Wave Model
A relationship between the pressure on the surface of the channel and the velocity of the channel wall is derived from the Rankine-Hugoniot relations for a shock front [40]. For soil solids, the plasma channel generated by a lightning pulse in the solid can be regarded as an expanding cylindrical piston. Now, consider a piston from which a continuous succession of shock fronts originate as its velocity varies in a continuous manner [16], as shown in Figure 5. The conservation of mass and conservation of momentum at the front is as follows: where ρ is the density, u is the velocity of soil, v is the velocity of the shockwave front, dρ is the change in density, du is the change in velocity, P is the pressure, and dP is the change in pressure. To solve the problem of the stress wave, substituting the Equation (7) into Equation (8) and omitting a low-order factor ρdu: The state equation of Murnaghan is deduced from the definition of bulk modulus [41]. Expanding K(P) by the Maclaurin series, and taking the first two terms of the expansion, the state equation is as follows: where K 0 and K 0 * are determined via experiment, and K 0 is the bulk modulus of soil without pressure.
Substituting Equation (10) into Equation (11) followed by integration: where Ψ and n are determined via experiment, K 0 = ρ 0 v 0 2 , n = 4s − 1, v 0 is the material speed of sound, ρ 0 is the density in the undisturbed state, and s is the material constant [42].
The governing equation for the problem stress wave would be obtained by substituting Equation (9) into Equation (12) as For the solid materials, substituting s = 1.5 into Equation (14) [42]: According to a simple mathematical transformation: The plasma channel expansion rate is: According to the geometric relationship of a cylinder: The volume of the plasma channel is differentiated against time: The expression of Equation (18) is transformed as follows: Substitute Equation (20) into (19) as follows: The shock wave pressure with the plasma channel volume would be obtained by substituting Equation (16) into Equation (21) as Equations (1), (2) and (22) can be solved numerically using the Runge-Kutta methods [43]. Then, the values of the shock wave pressure P and the plasma channel radius r ch (t), which vary with time, can be obtained.

Soil Deformation Model
The cylindrical shockwave comes from the lightning pulse. In a solid medium, stress from the lightning pulse is strongly correlated with stress from the blast load [34]. Figure 6 shows the deformation of soil by shockwave, which divides the soil into three parts: the breakdown channel, plastic zone, and elasticity zone. Soil is an elastomeric material, and the deformation of soil in the elasticity zone is recoverable as the shockwave gradually decays. The main compaction zone of soil is the plastic zone. According to the theory of soil plasticity mechanics, soil stress in the elasticity zone is [44] When r e →∞, where σ r is the radial stress, r is the soil radius, σ r2 is the radial stress when r = r 2 , σ θ is the cyclic stress, r e is the radius of a semi-infinite body, as shown in Figure 6, and σ r2 is the radial stress when r = r e . The plastic zone satisfies the Mohr-Coulomb yield criterion: where c is the soil cohesion and ϕ is the soil friction angle. Due to the short duration time of the lightning strike, the wave propagation in the solid medium resulting from the lightning strike can be treated as a point explosion in an elastic half-space (Chen et al., 2017) [34]. There are some experimental data that indicate the relationship between the loading rate effect with the soil strength, as shown in Figure 7 [45].
Explosives-induced liquefaction is deemed to result from residual pore water pressure increases in saturated soils; meanwhile, explosions can improve the soil yield strength compared to the static loads. Under the dynamic load, the deformation hysteresis effect and poor drainage conditions of soil are closely related to the soil shear strength [46,47]. Note: Ra is the ratio of dynamic strength to static strength of soil, t L is the loading time, and w is the water content.
Explosives can induce a change in the nature of the soil. In this paper, the coefficient λ is induced to improve the soil yield strength. Substituting the coefficient λ into Equation (9): The peak stress during the propagation of a stress wave can be given by Pan et al. [48]: where σ m is the initial peak pressure in the hole wall, r = r r 0 , r 0 is the hole radius, and a is the attenuation coefficient.
Substituting Equations (26) and (27) into Equation (29): According to the continuity of the stress boundary conditions, r = r 2 , the radial stress at the outer boundary of the plastic zone is equal to the radial stress at the inner boundary of the elastic zone. The radius of the plastic zone can be acquired by substituting r = r 2 and Equation (30) into Equation (31) as The soil plastic deformation zone by lightning strike can be acquired by substituting P into Equation (32) as The two-stage approach is employed to solve the radius of the soil plastic zone caused by lightning strikes. The first step, the Runge-Kutta method, is an efficient iterative method for solving nonlinear ordinary differential equations, and it is applied to solve the governing equation for shock wave pressure. The second step, in the process for calculating the radius of the soil plastic zone, has no complex numerical calculation method and it does not require iterative calculations. When the parameters are provided, the shockwave overpressure and the radius of the soil plastic zone are quickly obtained from the procedures.
However, lightning striking soil is an interdisciplinary problem about geotechnical and electrical engineering, as thus far, it has many unresolved issues, and the subject of lightning-induced soil mechanics is in the primary stage of exploration. In this paper, only the shock wave effect caused by lightning strikes is considered, which does not include the thermal effect induced by lightning strikes. When faced with real case problems, these methods can give partial results, i.e., problems caused by shock waves, but cannot handle problems caused by heat. The flowchart for solving the governing equations is as shown in Figure 8:

Comparison and Verification
As far as the authors know, neither experimental data nor an analytical model exists at present for the soil plastic zone induced by lightning strikes. Fortunately, Park et al. [30] provided some experimental data about shock wave pressure induced by a pulse discharge. Thus, the model is compared with the experimental data about shock wave pressure. The parametric analysis of the soil plastic zone induced by lightning is carried out in the next section.

Comparisons with Existing Experiment
The ground borehole expansion induced by pulse discharge technology was tested by Park et al. [30]. The corresponding parameters for calculation are given in Table 1 and Figure 9. According to the data of pulse current reported by Park et al. [30], they could obtain the current as a function of time for the period 0-0.5 ms. The shock wave pressure with time calculated by the proposed analytical model has some inaccuracy with the existing experiment results, while the peak pressure and the trends are in good agreement with the existing experimental results (see Figure 10), indicating that the present model can perform well in depicting shock wave pressure induced by a pulse discharge.

Parameter Analysis
The soil parameters are provided by Xiao et al. [49], as shown in Table 2.  Figure 11a presents the current waveform with time from 0 to 400 µs for three typical systems. Comparing the waveform at 8/20 µs with that at 8/80 µs, that at 8/80 µs has a slower decline. As shown in Figure 11b, the pressure increases followed by a decrease with time, and the radius of the plasma increases monotonically with time. For the waveform of 8/20 µs, the peak pressure is 414.3 MPa, the corresponding peak pressure is 414.2 MPa for the waveform of 8/80 µs, and the peak pressure is 342.6 MPa for the waveform of 20/80 µs. In other words, the wavefront time is a major factor to peak pressure, while the half-value time is a minor factor to peak pressure.

Influence of the Initial Radius of Plasma Channel
Lightning as a source of dynamic load acts on the soil, wherein the initial radius of plasma is an important component for channel expansion [50]. To study the influence of the initial radius of plasma on the soil plastic zone induced by lightning strikes, the RSPZ of the corresponding system is calculated with r i0 from 0.5 mm to 1.0 mm [50]. Figure 12a presents the peak of shockwave pressure curves with r i0 increasing from 0.5 mm to 1.0 mm for four typical systems. The theoretical peak pressure induced by lightning decreases monotonically with the initial radius of plasma. Therefore, Figure 12b presents the RSPZ curves increasing monotonically with r i0 . For r i0 = 0.5 mm, the RSPZs of i m = 20 kA, 40 kA, 80 kA, and 120 kA are 14.69 mm, 18.69 mm, 23.80 mm, and 27.25 mm, respectively, and the corresponding RSPZs are 21.73 mm, 28.93 mm, 37.45 mm, and 43.53 mm for r i0 = 1.0 mm. Namely, with r i0 increasing by 100% (from 0.5 mm to 1.0 mm), the RSPZ increases by 47.9-59.7%. In conclusion, the initial radius of plasma has a significant influence on the RSPZ induced by lightning strikes.

Influence of the Plasma Channel Length
The length of plasma is an important component for channel expansion. However, the plasma length is uncertain [16]. To study the influence of the plasma length on the soil plastic zone induced by lightning strikes, the RSPZ of the corresponding system is calculated with L from 0.01 to 1.0 m [50]. Figure 13a Figure 13b shows the RSPZ with plasma length, L, from 0.01 to 1.0 m for four typical systems. The RSPZ induced by lightning increases monotonically with plasma length. For L = 0.01 m, the RSPZs of i m = 20 kA, 40 kA, 80 kA, and 120 kA are 9.58 mm, 12.03 mm, 14.95 mm, and 16.81 mm, respectively, and the corresponding RSPZs are 16.57 mm, 21.34 mm, 27.38 mm, and 31.60 mm for L = 0.2 m. Namely, the slopes of r p with L range from 36.8 to 77.8. The corresponding RSPZs are 21.87 mm, 28.72 mm, 37.56 mm, and 43.96 mm for L = 1.0 m; in other words, the slopes of r p with L range from 6.6 to 15.5. In conclusion, the RSPZ and peak pressure rapidly increase followed by a slight increase with plasma length.

Influence of the Attenuation Coefficient
The wave propagation in solid medium resulting from the lightning strike can be treated as a point explosion in an elastic half-space ) [34]. The filling medium is a key influencing factor that can affect the effect of decoupled charge blasting [51][52][53]. To study the influence of the attenuation coefficient on the soil plastic zone induced by lightning, the RSPZ of the corresponding system is calculated with the at-tenuation coefficient, a, from 1.2 to 3.0 [53]. Figure 14 shows the RSPZ with attenuation coefficient, a, from 1.2 to 3.0 for nine typical systems. The RSPZ induced by lightning decreases exponentially with attenuation coefficient.

Influence of the Cohesion
To study the influence of cohesion on the change in soil properties induced by lightning strikes, the RSPZ of the corresponding system is calculated with c from 10 kPa to 30 kPa. Figure 15 presents the RSPZ curves with c increasing from 10 kPa to 30 kPa for the nine typical systems. The theoretical RSPZ induced by lightning decreases monotonically with the soil cohesion. As shown in Figure 15, for c = 10 kPa, the RSPZs of (a = 1.

Influence of the Internal Friction Angle
The internal friction angle is a key parameter for soil. To study the influence of internal friction angle on the change in soil properties induced by a lightning strike, the RSPZ of the corresponding system is calculated with ϕ from 10 • to 30 • . Figure 16 presents the RSPZ curves with ϕ increasing from 10 • to 30 • for the nine typical systems. The theoretical RSPZ induced by lightning decreases monotonically with the soil cohesion. As shown in Figure 16, for ϕ = 10 • , the RSPZs of (a = 1.

Conclusions
This paper assumes that the path of the plasma channel is a cylinder and does not have branches of channels. The assumptions in this paper are different from the actual situation. Under lightning impact, the development of the arc channel path in soil has an important impact on the improvement of lightning mitigation measures. However, the mechanism of arc channel development is not clear. With soil as a typical porous medium, the pore water in the soil will also affect the development of the arc channel. During the lightning discharge, a large amount of heat is generated. The thermal and mechanical fields are coupled, which affects the deformation of the soil and the soil properties. All the above-listed cases need to be studied one by one in the future.
This paper presents a semi-analytical solution for the peak pressure and the radius of the soil plastic zone induced by lightning strikes, based on a two-stage method. The Runge-Kutta method is applied to solve the governing equations for shockwave pressure based on mass conservation, momentum conservation, energy conservation, and the state equation of soil. The modified Mohr-Coulomb yield criterion is applied to solve the RSPZ. Therefore, if these procedures are to be implemented in real case problems by geotechnical engineers, they can also be quickly handled. The major factors dominating the shock wave pressure and the RSPZ are studied. The main conclusions are given as follows: (1) The r p induced by lightning decreases monotonically with c and ϕ. Under defined conditions, c increases by 200% (from 10 kPa to 30 kPa), the radius of the soil plastic zone decreases from 25.4% to 31.1%, ϕ increases by 200% (from 10 • to 30 • ), and the radius of the soil plastic zone decreases from 6.3-12.8%. With regard to r p caused by lightning strikes, c has a better effect on soil properties than ϕ does. (2) Increasing the r i0 can reduce the P and increasing r i0 has a nonnegligible effect on r p .
For example, with r i0 increasing by 100%, r p increases by 47.9-59.7%. (3) The P and r p increase monotonically with L. L has a significant influence on P and r p , especially when L is at a relatively low level; for example, when L increases from 0.01 m to 0.2 m, the slopes of r p with L range from 36.8 to 77.8, and the slopes of P with L range from 805.3 to 2315.8. (4) The channel radius r p significantly reduces with attenuation coefficient a. With attenuation coefficient a increasing by one third (from 1.2 to 1.6), the radius of the soil plastic zone decreases from 83.5% to 84.7%, and with attenuation coefficient a increasing by 25% from 2.4 to 3.0, the radius of the soil plastic zone decreases from 51.4% to 52.9%. Thus, it is necessary to improve the attenuation coefficient for decreasing soil damage induced by lightning strikes. (5) The wavefront time is a major factor, while the half-value time is a minor factor for the shock wave pressure induced by plasma explosives. For example, when i m = 20 kA, the peak pressure of 8/80 µs is 414.2 MPa and the peak pressure of 20/80 µs is 342.6 MPa.