A study on launch site ground of vehicle-mounted missile based on elastic layer theory

Abstract

With the widespread use of the vehicle-mounted missile launching system, the research of launch site ground becomes very important. The loading of missile launch system under the supporting disk is similar to inverted bell, and its radius and size of the missile launching are much more than a car. So the ready-made equivalent modulus formula under the action of circular uniformly distributed loading is not suitable for the calculation of launch site ground of vehicle-mounted missile. Therefore, based on the elastic layered system and the Gauss–Legendre numerical calculation method, a new calculation method of equivalent modulus is proposed, as well as the regression formulation for calculation of equivalent modulus of foundation is also shown. Based on the above analysis, a launch site ground mechanical model is presented. The thin plate model on elastic half-space foundation was formulated taking into account normal direction and the double-layer elastic system model was proposed for in tangential direction. Validity of the model was done by coring the programming calculation with the finite element numerical simulation for selected cases.

Introduction

With the development of detection technology and the increasing of precision strike capability, weapon system’s survival ability is becoming very important. Fast-moving launch system will undoubtedly improve the weapon system survivability. For trajectory missile, fast-moving missile launch system is more flexible than underground Wells launch system, which can move quickly in a short period of time from one position to another position. Fast-moving missile launch system has two ways, which are railway launch system and vehicle-mounted launch system. With the size of missile becoming smaller and smaller and the improvement of transportation equipment performance, vehicle-mounted launch system received more attention.

In the process of vehicle-mounted launch, constraint force fluctuations from ground to launch platform will cause strong vibration of components, which will affect the precision of a missile launch, seriously threatens the safety of launch. But, in the field of weapon launch presently, launch site ground model mainly focused on the soil model.^1–3 Proposed methods are not suitable for rigid concrete ground having layered structure. Therefore, based on elastic half-space foundation model and elastic layered system theory, this article reports the launch site ground mechanical model.

Equivalent modulus of the launch site ground

Model description

In the calculation of the deflection of the field vehicle-mounted launcher system, the elastic foundation plate theory^4–7 could be used. The floor structure was regarded as the elastic concrete thin plate surface layer and the foundation. The supporting layers which were under the foundation were characterized by only one modulus. In order to obtain this modulus, the multilayer structure is equivalent to a homogeneous body (elastic half-space) applying the elastic layer theory. Therefore, the parameter is called the equivalent modulus.

The literature^8–10 shows that the calculations of the equivalent modulus are based on the assumption of circular uniformly distributed loads, and the derived formula is only applicable under these conditions. The structure of the launcher supporting disk stiffness is much larger than the field, the loading under the rigid supporting disk is in a shape of inverted bell type and the radius of loading distribution and size of the missile launching are bigger than dimensions of vehicle. Therefore, the existing equivalent formulas for modulus are not suitable for the calculation of the launch site ground. In this section, based on the theory of elastic layered system, the formula for equivalent modulus under the action of rigid supporting disk is considered, which is using equivalent principle deflection. The equivalent modulus regression formula is established by means of a lot of number of calculations, which makes equivalent modulus to reach more accurate and has better extension.

Equivalent modulus of double-layer system

In Figure 1, q(r) represents the loading; h is the upper thickness of the double-layer elastic system; $E_1,{\mu }_1$ , $E_0,{\mu }_0$ express Young’s modulus and Poisson’s ratio of upper and lower layer, respectively; $E_t,{\mu }_t$ are Young’s modulus and Poisson’s ratio of the elastic half-space system; $w_f,w_t$ follow deflection of the double-layer elastic system and the half-layer elastic system, respectively.

Figure 1.

Double-layer elastic system and elastic half-space body under the rigid supporting disk.

As shown in Figure 1, double-layer elastic system is deflected by the rigid supporting disk. If this layer is considered as an elastic half-space layer,^11–13 a deflection value is also obtained under the same load. Equivalent modulus of the foundation could be captured, applying deflection equivalent principle.

Vertical displacement of the elastic half-space system under loading in a form of inverted bell

The formula for the loading represented by inverted bell is expressed as follows ¹⁴

$$q\left(r\right)=\{\begin{array}{c}\frac{1}{2}{q}_0{(1-\frac{r^2}{{\delta }^2})}^{-1/2}\hfill \\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}r⩽\delta \hfill \\ \multicolumn{1}{c}{r>\delta }\end{array}$$

(1)

where $q_0=Q/(\pi {\delta }^2)$ is the stress due to concentration of uniform load Q on the circular plane, $r,\delta$ represent the integral radius and the effect radius of the loading.

The Hankel integral transformation was used for the formula (1) and it was integrated into the integral expression of the surface displacement of the elastic half-space

$$w_t=\frac{2(1-{\mu }_t^2)}{E_t}\underset{0}{\overset{\infty }{\int }}{J}_0\left(\xi r\right)\overline{q}\left(\xi \right)d\xi$$

(2)

where $J_0(\xi r)$ expresses Bessel function; $\overline{q}(\xi )$ represents the Hankel integral transformation of the loading; $\xi$ is the intermediate variables.

The values of displacements within the zone followed by their radius were equal because of the rigid supporting disk, r = 0, the formula (2). According to the equation $\underset{r\to 0}{\lim }{J}_0(\xi r)=1$ , the deflection for the maximum value of loading vector is of

$$w_t=\frac{\pi (1-{\mu }_t^2)}{2E_t}{q}_0\delta$$

(3)

Vertical displacement of the double-layer system under loading in a form of inverted bell

According to the calculations for the elastic layer system,^12,15 the deflection value of the double-layer elastic system under the loading represented by equation (1) was determined. The associated items are not presented because of $A_2,C_2$ are equal to zero

$$C_1-(1-2{\mu }_1)A_1-D_1+(1-2{\mu }_1)B_1=1$$

(4)

$$C_1+2{\mu }_1{A}_1+D_1-2{\mu }_1{B}_1=0$$

(5)

$$\begin{array}{c}{e}^{\xi h}[C_1-(1-2{\mu }_1-\xi h)A_1]\hfill \\ -e^{-\xi h}[D_1+(1-2{\mu }_1+\xi h)B_1]\hfill \\ +e^{-\xi h}[D_2+(1-2{\mu }_0+\xi h)B_2]=0\hfill \end{array}$$

(6)

$$\begin{array}{c}{e}^{\xi h}[C_1+(2{\mu }_1+\xi h)A_1]\hfill \\ -e^{-\xi h}[D_1-(2{\mu }_1-\xi h)B_1]\hfill \\ -e^{-\xi h}[D_2-(2{\mu }_0-\xi h)B_2]=0\hfill \end{array}$$

(7)

$$\begin{array}{c}\frac{E_0(1+{\mu }_1)}{E_1(1+{\mu }_0)}\{e^{\xi h}[C_1+(2{\mu }_1+\xi h)A_1]\hfill \\ -e^{-\xi h}[D_1-(2{\mu }_1-\xi h)B_1]\}\hfill \\ +e^{-\xi h}[D_2-(1-\xi h)B_2]=0\hfill \end{array}$$

(8)

$$\begin{array}{c}\frac{E_0(1+{\mu }_1)}{E_1(1+{\mu }_0)}\{e^{\xi h}[C_1-(2-4{\mu }_1-\xi h)A_1]\hfill \\ +e^{-\xi h}[D_1+(2-4{\mu }_1+\xi h)B_1]\}\hfill \\ -e^{-\xi h}[D_2+(2-4{\mu }_0+\xi h)B_2]=0\hfill \end{array}$$

(9)

where $A_1,B_1,C_1,D_1,B_2,D_2$ are the integration constants obtained, based on boundary conditions.

It could be seen that the system of linear equations was quite complex. In this article, the integral coefficients were solved by the symbolic computation in MATLAB. After calculating the coefficients $A_1,B_1,C_1,D_1,B_2,D_2$ . And the coefficients $A_1,B_1,C_1,D_1$ were multiplied by the Hankel integral $\overline{q}(\xi )$ , respectively, and then substituted into formulas (4)–(9). The amount of deflection of the top layer was obtained by the integral operation eliminating $\xi$ . Employing z = 0, r = 0, the deflection under the loading for its maximum is of the $w_f$ formula

$$\begin{array}{c}{w}_f=-\frac{1+{\mu }_1}{E_1}{q}_0\delta \underset{0}{\overset{\infty }{\int }}\frac{\sin \left(\xi \delta \right)}{\xi }\hfill \\ [C_1-(2-4{\mu }_1)A_1+D_1+(2-4{\mu }_1)B_1]d\xi \hfill \end{array}$$

(10)

The establishment of equivalent relation

According to the equivalent principle deflection, the deflection of the elastic half-space system is equal to the deflection of the double-layer elastic system, $w_t=w_f$ , the formula can be simplified as follows

$$\begin{array}{c}{E}_t=-\pi (1-{\mu }^2)E_1/\hfill \\ \left\{(1+{\mu }_1)\underset{0}{\overset{\infty }{\int }}\frac{\sin \left(\xi \delta \right)}{\xi }[C_1-(2-4{\mu }_1)A_1+D_1+(2-4{\mu }_1)B_1]d\xi \right\}\hfill \end{array}$$

(11)

where $\pi =\mathrm{PI}$ .

The relationship between $E_t$ and $E_1$ , and $E_0$ and h can be established by formula (11). Because of the infinite integral of the Bessel function in the expression, the integral result was difficult to obtain. In this article, the Gauss–Legendre numerical integration method is used, and the discussion of integral method and integral upper limit is presented.

Examples and methods of validation

We have assumed that the mechanical parameters of two structural layers of the double-layer system were identical, which could be theoretically considered as an elastic half-space. In this way, it was possible to compare the deflection of the elastic half-space and the double layer. Therefore, the correctness of the equivalent parameters could be verified. Their values are listed in Table 1.

Table 1.

Mechanical parameters.

Parameters	h (cm)	$E_1=E_0(\mathrm{MPa})$	${\mu }_1={\mu }_0$	${\mu }_t$	$\delta$	$q_0(\mathrm{MPa})$
Values	30	2800	0.3	0.3	50	0.35

Taking the mechanical parameters in formula (10), the deflection of the double-layer elastic system was obtained

$$w_f=0.087472\left(\mathrm{mm}\right)$$

Employing formula (3) at the same values of the parameters considered, the deflection of the elastic half-layer system was also calculated

$$w_f=0.089339\left(\mathrm{mm}\right)$$

It could be seen that the two values determined were very close, the error was of −2.09%, expressing that the double-layer elastic equivalent modulus formula is correct.

The regression formula of comprehensive equivalent modulus

In order to make the calculation to be generality, a large number of rigid concrete fields with different structures are selected for calculation. Since Poisson’s ratio of the material changes little, the effect on the deflection is small. Assume that Poisson’s ratio of the upper and lower layer is consistent. The specific values are shown in Table 2.

Table 2.

The coefficient of the double-layer elastic foundation.

	Foundation	Soil
Thickness (cm)	10–40	–
Young’s modulus (MPa)	100–30,000	20–50
Poisson’s ratio	0.3	0.3

According to the supporting disk’s size of vehicle-mounted missile system and the load characteristics, assume that a radius of 50 cm loading was applied to the launch site ground surface, which form was formula (1) and the stress due to the load and plane was of $q_0=0.35\mathrm{MPa}$ , The equivalent modulus was calculated by formula (11), and then the results are processed by multivariate nonlinear regression. The regression curve was assumed to have the following form

$$\begin{array}{c}y=a+\mathrm{blnx}+c{\left(\mathrm{lnx}\right)}^2\hfill \\ +\mathrm{dhlnx}+\mathrm{eh}{\left(\mathrm{lnx}\right)}^2+f{\left(\mathrm{lnx}\right)}^3\hfill \end{array}$$

(12)

where $y=E_t/E_0,x=E_1/E_0$ , and a, b, c, d, e, and f fitting coefficients.

The fitting coefficient is shown in Table 3.

Table 3.

Fitting coefficients.

Fitting coefficients	a	b	c	d	e	f
Values	1	0.26	−0.1376	1.6	0.6	0.012

Figure 2 shows the results for the equivalent elastic modulus determined at h = 35 cm. It could be seen from the figure that the method proposed in this article makes up for the deficiency of the equivalent modulus calculated by the standard formula ¹⁶ when $E_1/E_0$ is less than 2.5. When the $E_1/E_0$ is greater than 8, course of the equivalent modulus was consistent with that of the standard formula, and the solution of the regression formula is located lower than that of the standard formula.

Figure 2.

Course of proportion of the integrated equivalent modulus for double-layer elastic system.

Mechanics model of launch site ground

There are two types of movement of the launching equipment relative to the launch site ground during the launching process: perpendicular and lateral movements. Therefore, during the launching of the equipment, the surface of the ground is mainly subjected to compression and friction. Under compression, the ground deflects and provides a vertical upward support force; the shear force produces a lateral displacement and provides resistance to movement in the opposite direction.

In the above, the foundation is equivalent to the elastic half-space applying the foundation’s equivalent modulus, then the launch site ground can be seen as a double-layer system. The compression responses of the field under compression and shear force were analyzed, respectively.

Normal motion control equation

Concrete road stiffness under compression is much larger than the foundation layer and the soil stiffness, which had a good ability for spreading the forces. According to the conventional test, it was found that the deflection of the surface layer was small under the loading emission, and the thickness of the surface layer usually reaches the following values 200–300 mm. The deflection of the surface layer was significantly lower than the thickness, so the problem of rigid concrete surface layer can be effectively captured based on the hypothesis of small deflection under typical exploitation conditions. Therefore, the elastic half-space model was adopted for the loading in a form of inverted bell, which can be seen in Figure 3.

Figure 3.

The elastic half-space model.

Assuming h, $E_1,{\mu }_1$ express the layer’s thickness of Young’s modulus and Poisson’s ratio; $E_0,{\mu }_0$ represent the foundation’s Young’s modulus and Poisson’s ratio; the compression loading q(r) was applied to the surface layer; the top surface of the foundation was selected to be the zone for the pressure p(r).

According to the Kirchhoff layer theory, the following assumptions were made:

1. The normal line of the layers before and after deformation remains a straight line and perpendicular to the surface.

2. The normal line of the layers neither elongates nor shortens, and the points are not parallel to the surface displacement.

3. In the deformation process, the layer and the foundation were always in close contact, the vertical displacement of the top of the foundation and that of the surface layer is equal.

Control equation and integral solution

The differential equation for plate bending is presented in a lot of papers. ¹⁷ Using the cylindrical coordinates, the elastic surface equation of the plate is of

$$\mathrm{D}{\nabla }^2{\nabla }^{2}w\left(r\right)=q\left(r\right)-p\left(r\right)$$

(13)

where ${\nabla }^2={\partial }^2/{\partial }_r^2+\partial /(r\times {\partial }_r)+{\partial }^2/{\partial }_z^2$ Laplace operator; $\mathrm{D}=E{h}^3/[121-{\mu }^2]$ bending rigidity of the plate; $w$ represents deflection of the plate.

w(r) and p(r) are unknown quantities. According to the assumptions, the vertical displacement of the foundation top surface was equal to the vertical displacement of the surface layer. From the literature,^11–13 the elastic half-space displacement expression under the axisymmetric vertical load is taken as follows

$$\mathrm{w}=\frac{2(1-{\mu }^2)}{E}\underset{0}{\overset{\infty }{\int }}{J}_0\left(\xi r\right)\overline{p}\left(\xi \right)d\xi$$

(14)

Substituting equation (14) into equation (13), and performing Hankel integral transformation to the final equation, which would get $\overline{p}(\xi )$ expression, and applying $\overline{p}(\xi )$ into equation (14), we can get the plate deflection

$$\mathrm{w}=\frac{2(1-{\mu }^2)}{E}\underset{0}{\overset{\infty }{\int }}\frac{J_0\left(\xi r\right)}{1+l^3{\xi }^3}\overline{q}\left(\xi \right)d\xi$$

(15)

where $l=[2D1-{\mu }^2/E{]}^{1/3}$ .

Deflection of the launch site ground under the action of uniform load and loading in a form of inverted bell

Conditions for uniform loading

$$\mathrm{q}\left(\mathrm{r}\right)=\{\begin{array}{c}{q}_0\hfill \\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}r⩽\delta \hfill \\ \multicolumn{1}{c}{r>\delta }\end{array}$$

(16)

The formulas (1), (16) Hankel integral transformation, and formula (15) enabled to obtain:

Deflection of the launch site ground under the action of uniform loading

$$\mathrm{w}=\frac{2(1-{\mu }^2)}{E}{q}_0\delta \underset{0}{\overset{\infty }{\int }}\frac{J_0\left(\xi r\right)J_1\left(\delta \xi \right)}{\xi (1+l^3{\xi }^3)}d\xi$$

(17)

Deflection of the launch site ground under the action of the loading represented by inverted bell shape

$$\mathrm{w}=\frac{2(1-{\mu }^2)}{E}{q}_0\delta \underset{0}{\overset{\infty }{\int }}\frac{J_0\left(\xi r\right)\sin \left(\delta \xi \right)}{\xi (1+l^3{\xi }^3)}d\xi$$

(18)

Tangential motion control equation

Control equation and integral solution

There was an axial symmetry of the role of the unidirectional horizontal loading acted on the ground surface and having radius of δ. There was no loading outside the radius, which can be seen in Figure 4.

Figure 4.

The double-layer elastic system with unidirectional horizontal loading.

When m = 1, the expression of the stress and displacement state components of the elastic half-space under axisymmetric circular distribution unidirectional horizontal loading could be obtained

$$\{\begin{array}{c}u=\frac{1+\mu }{2E}(U_2-V_0)\cos \theta \hfill \\ \multicolumn{1}{c}{\begin{array}{c}v=\frac{1+\mu }{2E}(U_2+V_0)\sin \theta \hfill \\ \multicolumn{1}{c}{w=-[\frac{1+\mu }{E}\underset{0}{\overset{\infty }{\int }}{J}_1\left(\xi r\right)\{[C-(2-4\mu -\xi z)A]e^{\xi z}}\end{array}}\\ \multicolumn{1}{c}{+[D+(2-4\mu +\xi z)B]e^{-\xi z}\}]\cos \theta }\\ \multicolumn{1}{c}{\begin{array}{c}{\sigma }_z=-[\underset{0}{\overset{\infty }{\int }}\xi J_1\left(\xi r\right)\{[C-(1-2\mu -\xi z)A]e^{\xi z}\hfill \\ \multicolumn{1}{c}{+[D+(1-2\mu +\xi z)B]e^{-\xi z}\}]\cos \theta }\\ \multicolumn{1}{c}{\begin{array}{c}{\tau }_{z\theta }=\frac{1}{2}(T_2+S_0)\sin \theta \hfill \\ \multicolumn{1}{c}{{\tau }_{\mathrm{zr}}=\frac{1}{2}(T_2-S_0)\cos \theta }\end{array}}\end{array}}\end{array}$$

(19)

where

$$\begin{array}{c}{U}_2=\underset{0}{\overset{\infty }{\int }}{J}_2\left(\xi r\right)\hfill \\ \{[C+(1+\xi z)A+2F]e^{\xi z}-[D-(1-\xi z)B-2I]e^{-\xi z}\}d\xi \hfill \\ V_0=\underset{0}{\overset{\infty }{\int }}{J}_0\left(\xi r\right)\hfill \\ \{[C+(1+\xi z)A-2F]e^{\xi z}-[D-(1-\xi z)B+2I]e^{-\xi z}\}d\xi \hfill \\ T_2=\underset{0}{\overset{\infty }{\int }}\xi J_2\left(\xi r\right)\hfill \\ \{[C+(2\mu +\xi z)A+F]e^{\xi z}-[D-(2\mu -\xi z)B-I]e^{-\xi z}\}d\xi \hfill \\ S_0=\underset{0}{\overset{\infty }{\int }}\xi J_2\left(\xi r\right)\hfill \\ \{[C+(2\mu +\xi z)A-F]e^{\xi z}-[D-(2\mu -\xi z)B+I]e^{-\xi z}\}d\xi \hfill \end{array}$$

• Surface stress boundary conditions (z = 0)

In the cylindrical coordinate system, the stress boundary conditions were expressed of

$${\tau }_{z\theta }=\{\begin{array}{c}q\left(r\right)\sin \theta \hfill \\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}(0⩽r⩽\delta )\hfill \\ \multicolumn{1}{c}{(\delta ⩽r⩽\infty )}\end{array}$$

(20)

$${\tau }_{\mathrm{zr}}=\{\begin{array}{c}-q\left(r\right)\cos \theta \hfill \\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}(0⩽r⩽\delta )\hfill \\ \multicolumn{1}{c}{(\delta ⩽r⩽\infty )}\end{array}$$

(21)

$${\sigma }_z=0(0⩽r⩽\delta )$$

(22)

In order to facilitate the calculation, the above expression was rewritten

$$\{\begin{array}{c}\frac{{\tau }_{z\theta }}{\sin \theta }+\frac{{\tau }_{\mathrm{zr}}}{\cos \theta }=0\hfill \\ \multicolumn{1}{c}{\frac{{\tau }_{z\theta }}{\sin \theta }-\frac{{\tau }_{\mathrm{zr}}}{\cos \theta }=2q\left(r\right)}\\ \multicolumn{1}{c}{{\sigma }_z=0(0⩽r⩽\infty ,\mathrm{z}=0)}\end{array}$$

(23)

Using equation (19) in formula (23), and Hankel integral conversion as well as $F_i=\overline{F}\overline{q}(\mathrm{\epsilon })(F=A,B,C,D,i=1,2,\dots ,N)$ , and employing ${\overline{A}}_i,{\overline{B}}_i,{\overline{C}}_i,{\overline{D}}_i$ to $A_i,B_i,C_i,D_i(i=1,2,\dots ,N)$ the following system of equation

$$\{\begin{array}{c}{C}_1+2{\mu }_1{A}_1+F_1+D_1-2{\mu }_1{B}_1-I_1=0\hfill \\ \multicolumn{1}{c}{C_1+2{\mu }_1{A}_1-F_1+D_1-2{\mu }_1{B}_1+I_1=2}\\ \multicolumn{1}{c}{C_1-(1-2{\mu }_1)A_1-D_1-(1-2{\mu }_1)B_1=0}\end{array}$$

(24)

where $\overline{q}(\xi )={\int }_0^{\infty }\mathrm{rq}(r)J_0(\xi r)\mathrm{dr}$ Hankel integral conversion of the loading concentration.

• Interlayer bonding conditions (z = h)

According to the continuous interface assumption, the upper and lower interfaces should have the same stress state, the stress and displacement at the junction are continuous, with the expression described as

$$\{\begin{array}{cc}{w}_{z1}=w_{z2}\hfill & \begin{array}{cc}{\mu }_{z1}={\mu }_{z2}\hfill & v_{z1}=v_{z2}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{{\sigma }_{z1}={\sigma }_{z2}}& \begin{array}{cc}{\tau }_{\mathrm{zr}1}={\tau }_{\mathrm{zr}2}\hfill & {\tau }_{z\theta 1}={\tau }_{z\theta 2}\hfill \end{array}\hfill \end{array}$$

(25)

The stress and displacement components in formula (19) were introduced in formula (25), and Hankel integral changes were performed on both sides of each equation. By adding or subtracting the equations, six linear equations were obtained

$$\begin{array}{c}\multicolumn{1}{c}{\begin{array}{c}{e}^{\xi h}[C_1-(1-2{\mu }_1-\xi h)A_1]\hfill \\ \multicolumn{1}{c}{-e^{-\xi h}[D_1+(1-2{\mu }_1+\xi h)B_1]}\end{array}}\\ \multicolumn{1}{c}{\begin{array}{c}=e^{\xi h}[C_2-(1-2{\mu }_0-\xi h)A_2]\hfill \\ \multicolumn{1}{c}{-e^{-\xi h}[D_2+(1-2{\mu }_0+\xi h)B_2]}\end{array}}\end{array}$$

(26)

$$\begin{array}{c}\begin{array}{c}{e}^{\xi h}[C_1+(2{\mu }_1+\xi h)A_1]\hfill \\ +e^{-\xi h}[D_1-(2{\mu }_1-\xi h)B_1]\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{c}=e^{\xi h}[C_2+(2{\mu }_0+\xi h)A_2]\hfill \\ \multicolumn{1}{c}{+e^{-\xi h}[D_2-(2{\mu }_0-\xi h)B_2]}\end{array}}\end{array}$$

(27)

$$\begin{array}{c}\begin{array}{c}\frac{E_0(1+{\mu }_1)}{E_1(1+{\mu }_0)}\{e^{\xi h}[C_1+(1+\xi h)A_1]\hfill \\ -e^{-\xi h}[D_1-(1-\xi h)B_1]\}\hfill \end{array}\hfill \\ =e^{\xi h}[C_2+(1+\xi h)A_2]-e^{-\xi h}[D_2-(1-\xi h)B_2]\hfill \end{array}$$

(28)

$$\begin{array}{c}\begin{array}{c}\frac{E_0(1+{\mu }_1)}{E_1(1+{\mu }_0)}\{e^{\xi h}[C_1-(2-4{\mu }_1-\xi h)A_1]\hfill \\ +e^{-\xi h}[D_1+(2-4{\mu }_1+\xi h)B_1]\}\hfill \end{array}\hfill \\ =e^{\xi h}[C_2-(2-4{\mu }_1-\xi h)A_2]-e^{-\xi h}[D_2+(2-4{\mu }_0+\xi h)B_2]\hfill \end{array}$$

(29)

$$F_1{e}^{\xi h}-I_1{e}^{-\xi h}=F_2{e}^{\xi h}-I_2{e}^{-\xi h}$$

(30)

$$\frac{E_0(1+{\mu }_1)}{E_1(1+{\mu }_0)}(F_1{e}^{\xi h}+I_1{e}^{-\xi h})=F_2{e}^{\xi h}-I_2{e}^{-\xi h}$$

(31)

• Additional condition of elastic half-space $(\mathrm{z}\to \infty )$

$$A_2=0,C_2=0,F_2=0$$

(32)

Equations (23), (24), (26)–(32) formed an algebraic system of 12 equations. Through the computation applying symbols, the linear algebraic equations were solved and the integral coefficients of each layer were obtained. As can be seen from the matrix coefficient, the integral coefficient is only related to the intermediate variable ξ, the remaining variables are known. After obtaining the integral coefficient, multiply the Hankel integral of the loading and take it into the formula (19), then the theoretical solution of double-layer system under the uniaxial horizontal load with symmetric circular distribution of the load concentration was obtained.

Tangential displacement of the launch site ground surface under the action of uniform and not uniform loading

(1), (16) Hankel integral transformation, and taken into the tangential control equation of the theoretical solution, z = 0

$$\mathrm{u}=\frac{1+{\mu }_1}{2E}(U_2-V_0)$$

(33)

Deflection of the launch site ground under the action of uniform loading was expressed as follows

$$\{\begin{array}{c}{U}_2=q_0\delta \underset{0}{\overset{\infty }{\int }}\frac{J_1\left(\xi \delta \right)J_2\left(\xi \delta \right)}{\xi }(C+A+2F-D+B+2I)d\xi \hfill \\ \multicolumn{1}{c}{V_0=q_0\delta \underset{0}{\overset{\infty }{\int }}\frac{J_1\left(\xi \delta \right)J_0\left(\xi \delta \right)}{\xi }(C+A-2F-D+B-2I)d\xi }\end{array}$$

(34)

Deflection of the launch site ground under the action of loading in a form of inverted bell was presented

$$\{\begin{array}{c}{U}_2=q_0\delta \underset{0}{\overset{\infty }{\int }}\frac{J_1\left(\xi \delta \right)\sin \left(\xi \delta \right)}{2\xi }(C+A+2F-D+B+2I)d\xi \hfill \\ \multicolumn{1}{c}{V_0=q_0\delta \underset{0}{\overset{\infty }{\int }}\frac{J_1\left(\xi \delta \right)\sin \left(\xi \delta \right)}{2\xi }(C+A-2F-D+B-2I)d\xi }\end{array}$$

(35)

Numerical integral method with infinite integral of Bessel function

In the above sections, the expressions of stress and displacement had infinite integral of Bessel oscillation function and exponential function integrals in the solution of elastic half-space and double-layer elastic system, which can only be calculated by numerical integration. There are two types of numerical methods for the theory of elastic layer. One is the Simpson numerical integration method and the other is the Gaussian numerical integration method. In order to verify the accuracy and validity of the numerical method, the following expression should be considered

$$\underset{0}{\overset{\infty }{\int }}\frac{J_0(r\xi )J_1(\delta \xi )}{1+{\xi }^3}d\xi$$

(36)

Taking r = 0.9 and δ = 0.4, the plot of the integrand can be reached as shown in Figure 4.

It could be seen from Figure 5 that the integrand with two Bessel functions was oscillatory and rapidly reached zero. At ξ = 6, the value of the integrand is diminished up to −9.879e–5, and while the upper limit $x_s=10\delta /h_1\approx 17.39$ , it could be seen that the value of the upper limit of the integral is correct.

Figure 5.

The plot of the integrand function.

In order to compare the computational efficiency of the Gauss–Legendre integral and the Simpson integral for solving the launch site ground problem, numerical integration of the expression (36) was performed using the two methods. The results were shown in Table 4.

Table 4.

Integration results of Simpson and Gauss–Legendre method.

Simpson method	Number of nodes	10	30	40	50	70
	Approximation	0.1117	0.1066	0.1066	0.1065	0.1065
	Difference from the previous one	–	−4.78%	0%	−0.094%	0%
Gauss–Legendre method	Number of nodes	10	15	20	25	30
	Approximation	0.1097	0.1063	0.1063	0.1065	0.1065
	Difference from the previous one	–	−1.5%	0.28%	−0.094%	0%

As can be seen from Table 4, when the number of nodes was 30, the approximate value of Simpson’s integral method began to be stable, and the number of calculations was close to 90 times. When the number of nodes reached 15, the approximation of the Gauss–Legendre method began to stabilize, and the number of calculations was 15. It could be seen that the Gauss–Legendre method was superior to the Simpson method with respect to time needed to achieve final results and accuracy. This article adopted Gauss–Legendre numerical integration method.

Model validation

Validity of the model was proved by comparing the programming calculation with the finite element numerical simulation for selected cases. The type of the load was formula (16), while the acting radius of the loading was of δ = 440 mm, and resultant pressure was of $q_0=0.4\mathrm{MPa}$ . Parameters of a launch site ground are listed in Table 5. The equivalent modulus of the foundation was calculated by the equivalent modulus method proposed in this article. The calculated modulus was equivalent to that of the concrete surface in Table 6.

Table 5.

Mechanical parameters of the launch site ground.

The layer number	The name of structural material	Young’s modulus (MPa)	Poisson’s ratio	Thickness (cm)
1	PCC	30,000	0.15	22
2	Water stable gravel	4200	0.35	20
3	Lime soil	1800	0.35	20
4	Compacted soil foundation	40	0.35	–

PCC: plain cement concrete.

Table 6.

Mechanical parameters of the equivalent launch site ground.

The layer number	The name of structural material	Young’s modulus (MPa)	Poisson’s ratio	Thickness (cm)
1	PCC	30,000	0.15	22
2	Elastic half-layer foundation	291.52	0.35	–

PCC: plain cement concrete.

The normal control equation was determined by MATLAB programming. The finite element calculation was carried out in ABAQUS software. Figure 6 presents a surface deflection curve of the cross section of the field under uniform loading, where the deflection value was negative for down deflection and its value was the theoretical approach for the numerical solution. From the figure, it could be seen that the deflection was dependent on tangential distance, reaching the minimum for the tangential distance of 0 mm. The deflection increased with tangential distance increasing, this is consistent with the results of Johnson. ¹⁴ The deflection value of the finite element solution at the load center and the loading edge were 0.385 and 0.32 mm, the deflection value of the theoretical solution at the load center and the loading edge were 0.397 and 0.344 mm. If the finite element numerical solution was used as the reference, the calculation errors of the center and the edge were 0.462% and 7.5%, respectively, which indicates on small differences.

Figure 6.

Deflection of the launch site ground under uniform loading.

Figure 7 shows a three-dimensional (3D) shell view of the ground deflection calculated by MATLAB. Figure 8 illustrates 3D deformation diagram from the finite element numerical calculation. From the surface of the floor, deformation map could be more intuitive to see that both were axisymmetric deformation.

Figure 7.

A shell of the ground deflection in 3D coordinate system.

Figure 8.

The ground deflection in 3D coordinate system.

Conclusion

The main work of this article is to build the mechanics model of launch site ground.

Based on the features of the contact structure between the supporting disk and the launch site ground, the equivalent modulus calculation method of the elastic layer system was elaborated and verified. The calculation program of the equivalent modulus of the foundation of the launch site ground was compiled. The regression formula for the modulus was obtained employing a large number of different structures.

Based on the above analysis and according to the characteristics of the launching equipment relative and the launch site ground, the mechanical model of the launch site ground was proposed. The elastic half-space elastic model and the elastic double-layer model were established in the normal direction and tangential direction, respectively. The deflection of the ground under a symmetrical uniform loading was calculated and compared with the results obtained by the use of ABAQUS software. The results showed that the calculation errors for the theoretical solution were 0.462% and 7.5%, respectively.