Research on the identification of modulus parameters for rigid pavement panel foundations

<!-- metadata commented in wiki content


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<big>Research on the identification of modulus parameters for rigid pavement panel foundations</big></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Hui-min Cao<sup>1,2</sup>;Han-jun Zhao<sup>2,3</sup>;Xin-e Yan<sup>1*</sup></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
(1.Xi’an Traffic Engineering Institute(China);2. Northwestern Polytechnical University(China);3.Shanxi Transportation Research Institute Group CO.Ltd(China))</div>
Fund Project:1.Young and middle-aged fund project at Xi’an Traffic Engineering Institute(2023KY-39); 2.The Youth Innovation Team of Shaanxi Universities(Key Technology Innovation Team for Urban Rail Transit Track Bed); 3.The Scientific Research Program Funded by Education Department of Shaanxi Provincial Government (Program No.23JK0532).
-->
==Abstract==

Based on the Pasternak two-parameter foundation model, this study establishes a dynamic model for analyzing the deflection of slabs under moving loads. The variational method and reciprocal equal work theorem are employed to obtain the deflection solution. Additionally, the Matlab program is utilized with practical examples to calculate the foundation response modulus under actual conditions when the pavement panel simultaneously transmits shear and bending forces. By fitting measured dynamic deflections with theoretical predictions using the least square method, the response modulus <math display="inline">K</math> of the foundation is determined. This demonstrates that our proposed model effectively captures the realistic behavior of rigid pavements..

'''Keywords''': Pasternak model, least squares method, foundation reaction modulus

==1. Introduction==

With the rapid development of China's civil aviation industry, higher performance and safety requirements for airport pavement are inevitable.Currently, the commonly used foundation models include the E.Winkler foundation model and the elastic half-space foundation model. Although the Winkler foundation model theory is simple and intuitive, requiring only one elastic parameter to express ideal elastic foundation characteristics and relating displacement of any point solely to applied stress at that point,it fails to consider soil-ground continuity. In contrast, the elastic half-space foundation model overemphasizes stress diffusion between soil and ground, making calculations more cumbersome and impractical for engineering applications. Consequently, mechanics experts have developed several more reasonable foundation calculation models through extensive experimentation. Zhang et al. <span id='citeF-1'></span> [[#cite-1|[1]]] and Chen and Liu <span id='citeF-2'></span>[[#cite-1|[2]]] utilized the Kelvin viscoelastic model under moving load conditions to obtain a relatively accurate foundation reaction modulus. The viscoelastic Winkler dynamic model verified by Xing <span id='citeF-3'></span>[[#cite-3|[3]]] better reflects actual airport operations. Patil et al. <span id='citeF-4'></span>[[#cite-4|[4]]] employed Pasternak's two-parameter soil medium as a model to investigate material nonlinearity effects on pavement response in supporting soil medium. Kumar et al. <span id='citeF-5'></span>[[#cite-5|[5]]] applied a novel finite element-based cyclic response model for rigid pavement design while establishing the period range of moving load within rigid pavement parameters. Airport pavements consist of multiple concrete panels working collectively. This paper selects the Pasternak two-parameter foundation model based on independent parameters:foundation reaction modulus <math display="inline">K</math>  and shear modulus <math display="inline">G</math>; thereby establishing a dynamic rigid pavement model considering joint load transfer capacity.

==2. The formulation of dual-panel motion equations==

===2.1 The formulation of the equation of motion===

Based on the two-parameter foundation, a dynamic model of rigid thin plate pavement on the two-parameter foundation is established, considering boundary shear and bending under moving load <span id='citeF-6'></span>[[#cite-6|[6]]], where  <math display="inline">\mathit{\boldsymbol{Q}}</math> is the foundation reaction force. The dynamic equation of the two-parameter foundation incorporates two independent parameters: reaction modulus <math display="inline">K</math> and shear modulus <math display="inline">G</math>. [[#img-1|Figure 1]] illustrates the schematic diagram of the two-parameter foundation, while [[#img-2|Figure 2]] depicts the forces acting on the shear layer <span id='citeF-7'></span>[[#cite-7|[7]]].

<div id='img-1'></div>
{| class="wikitable" style="margin: 0em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
<!--|style="text-align: center;padding:10px;"| [[Image:Draft_Cao_167492020-picture- 1.svg|center|300px]]
|--->
| style="text-align: center;padding:10px;"| [[File:Cao_et_al_2024a_7391_Draft_Cao_167492020-picture-_1.png|400px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 1'''. Schematic of a two-parameter foundation model
|}


<div id='img-2'></div>
{| class="wikitable" style="margin: 0em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Hui-min_340971489-image2.png|400px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 2'''. Forces on the shear layer
|}


According to the Kirchhoff thin plate theory <span id='citeF-8'></span>[[#cite-8|[8]]] and the theory of elastic mechanics, we establish the motion equation for plates subjected to a two-parameter foundation under dynamic loading as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> D{\nabla }^{4}w+M\frac{{\partial }^{2}w}{\partial {t}^{2}}+Kw-G{\nabla }^{2}w=P  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(1)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math>D=\frac{E{h}^{3}}{12\left( 1-{\mu }^{2}\right) }</math>
|}
|
|}

where:

<math display="inline">D</math> is the flexural stiffness of the plate,

<math display="inline">\mathrm{\mu}</math>, the Poisson's ratio

<math display="inline">M=\rho h</math>, the areal density

<math display="inline">{\nabla }^{2}=\frac{{\partial }^{2}}{\partial {x}^{2}}+\frac{{\partial }^{2}}{\partial {y}^{2}}</math>, Laplace operator in two dimensions

<math display="inline">K</math>, foundation reaction modulus

<math display="inline">G</math>, foundation shear modulus, and

<math display="inline">w</math> is the vertical deflection of the pavement panel.

The aforementioned equation represents the differential equation of motion for the track panel system subjected to a general moving load, denoted as <math display="inline">P\left( x,y,t\right)</math>. In the context of the investigated form of moving load in this study, the equality <math display="inline">P\left( x,y,t\right)</math>  can be expressed as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> P\left( x,y,t\right) ={P}_{0}\delta \left( x-{v}_{x}t\right) \delta \left( y- \eta \right)  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(2)
|}


The function <math display="inline">\delta \,</math> is represented by <math display="inline">\delta \left( x\right)</math>, and the uniform gliding speed of the aircraft is denoted by <math display="inline">{v}_{x}</math>.

===2.2 Boundary conditions (shear and bending)===

The plate dimensions are defined as length <math display="inline"> a </math> and width <math display="inline"> b </math>, with simultaneous transfer of shear force and bending moment between the plates in the actual working state. The boundary shear force is generated by the spring support force, where the elastic support rigidity coefficient <math display="inline">{k}_{1}\,</math> represents the boundary condition parameter. Similarly, the bending moment is generated by the spring-supporting moment, with the elastic supporting rigidity coefficient <math display="inline">{k}_{2}</math> representing another boundary condition parameter. Assuming that <math display="inline">Y=</math><math>0,b</math>  corresponds to the side where the seam is located, we can express its general formula for boundary conditions as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  {V}_{y=0,b}=D{\left[ \frac{{\partial }^{3}w}{\partial {y}^{3}}+\frac{\left( 2-\mu \right) {\partial }^{3}w}{\partial {x}^{2}\partial y}\right] }_{y=0,b}=-{k}_{1}{\left[ w\right] }_{y=0,b}</math> 
|-
|style="text-align: center;" |<math>{M}_{y\left| y=0,b\right.} =D\left( \frac{{\partial }^{2}w}{\partial {y}^{2}}+ \mu \frac{{\partial }^{2}w}{\partial {x}^{2}}\right) \Big| y=0,\quad b= -{k}_{2}\left| \frac{\partial w}{\partial y}\right| _{y=0,b} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(3)
|}

==3. The solution to the equation of motion==

===3.1 Solving time factors===

The variational formulation of Eq. (1) for plate motion is as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>\int_{0}^{a}\int_{0}^{b}\left\{ D{\nabla }^{4}w+M\frac{{\partial }^{2}w}{\partial {t}^{2}}+Kw-G{\nabla }^{2}w-P\right\} \delta wdxdy=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(4)
|}

The lateral load concentration stands out among them. <math display="inline">P=</math><math>P\left( x,y,t\right)</math>, <math display="inline">w=w\left( x,y,t\right)</math>.

Select a function in the specified format as the solution for the variational equation mentioned above

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \left\{ \begin{matrix}\ w(x,y,t)=T(t)W(x,y)    
\\ P(x,y,t)=B(t)W(x,y)\end{matrix} \right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(5)
|}


The mode function of the plate is represented by <math display="inline">W\left( x,y\right)</math>, and the coefficients of the deflection function <math display="inline">w\left( x,y,t\right)</math>  and load <math display="inline">P\left( x,y,t\right)</math>  are denoted as <math display="inline">T\left( t\right)</math>  and <math display="inline">B\left( t\right)</math>,  respectively. The variation corresponding to the deflection function <math display="inline">w\left( x,y,t\right)</math>  is as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \delta w=W\left( x,y\right) T\left( t\right)  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(6)
|}

Substituting Eqs. (5) and (6) into (4) gives:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \left( t\right) +{{w}_{mn}}^{2}T\left( t\right) =\frac{1}{M}B\left( t\right)   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(7)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> {{w}_{mn}}^{2}  =\frac{1}{M}\frac{{I}_{1}}{{I}_{2}}</math>
|-
|style="text-align: center;" | <math>{I}_{1}  =\int_{0}^{a}\int_{0}^{b}\left[ D{\nabla }^{4}W\left( x,y\right) +KW\left( x,y\right) -G{\nabla }^{2}W\left( x,y\right) \right] W\left( x,y\right) dxdy</math>
|-
| style="text-align: center;" |<math>{I}_{2}  =\int_{0}^{a}\int_{0}^{b}{W}^{2}\left( x,y\right) dxdy </math>
|}
|
|}

The computable expression for <math display="inline">T(t)</math> is as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  T\left( t\right) =\left( {c}_{1}\mathrm{cos}\,{w}_{mn}t+{c}_{2}\mathrm{sin}\,{w}_{mn}t\right) +\frac{1}{M{w}_{mn}}\int_{0}^{t}B\left( \tau \right) \mathrm{sin}\,{w}_{mn}\left( t- \tau \right) d\tau </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(8)
|}

where, <math display="inline">{w}_{mn}</math> represents the natural frequency of undamped free vibration for a thin plate. The first term in Eq. (8) denotes the thin plate's free vibration, while the second term represents its forced vibration.

Among them, <math display="inline"> c_1 </math>  and <math display="inline"> c_2 </math> are constants that depend on the initial conditions of motion, while <math display="inline">B(t)</math>  is dependent on the load characteristics. In this study, the load is considered as a concentrated load <math display="inline"> P_0 </math>, neglecting its mass. Consequently, <math display="inline">B(t)</math> can be expressed as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> B\left( t\right) =\frac{{P}_{0}W\left( {v}_{x}t,\eta \right) }{\displaystyle\int_{0}^{a}\int_{0}^{b}{W}^{2}\left( x,y\right) dxdy}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(9)
|}

According to Eqs. (8) and (9), the dynamic deflection of the plate under the influence of a moving concentrated load P, with zero initial conditions (<math display="inline"> c_1 =c_2</math>), can be determined as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  w\left( x,y,t\right) =\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{W\left( x,y\right) }{M{w}_{mn}\displaystyle \int_{0}^{a}\int_{0}^{b}{W}^{2}\left( x,y\right) dxdy}</math><math>\int_{0}^{t}{P}_{0}W\left( {v}_{x}\tau ,\eta \right) \mathrm{sin}\,{w}_{mn}\left( t- \tau \right) d\tau </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(10)
|}

===3.2 Solving coordinate factors===

In order to address the analytical solution for the bending behavior of diverse boundary plates subjected to moving loads, a four-sided simply supported rectangular plate with dimensions, geometry, and material properties that precisely match those of the actual system is considered as the fundamental model <span id='citeF-9'></span>[[#cite-9|[9]]], as depicted in [[#img-3|Figure 3]].

<div id='img-3'></div>
{| class="wikitable" style="margin: 0em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Hui-min_340971489-image3.png|492px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 3'''. Basic system of board
|}


The plate is subjected to a unit concentration force at the flow coordinate <math display="inline">\left( \zeta ,\eta \right)</math>, and the displacement is solved using a trigonometric series

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> {w}_{1}=\frac{4}{abD}\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{\mathrm{sin}\,{k}_{m}\zeta \mathrm{sin}\,{k}_{n}\eta }{({{k}_{m}}^{2}+{{k}_{n}}^{2}{)}^{2}}\mathrm{sin}\,{k}_{m}x\mathrm{sin}\,{k}_{n}y  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(11)
|}

where, <math>{k}_{m}=\frac{m\pi }{a},\quad {k}_{n}=\frac{n\pi }{b}</math>.

The solution  <math display="inline">{w}_{1}</math>  is referred to as the fundamental solution corresponding to the rectangular plate system.

The boundary condition can be simplified as a simply supported boundary by incorporating the corresponding support settlement and constraint moment under actual boundary conditions. The additional support settlement is assumed to be:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> {W}_{x0}=\sum _{n=1}^{\infty }{A}_{n}\mathrm{sin}\,{k}_{n}y\qquad ,\qquad {W}_{xa}=\sum _{n=1}^{\infty }{B}_{n}\mathrm{sin}\,{k}_{n}y,  </math>
|-
| style="text-align: center;" | <math>{W}_{y0}=\sum _{m=1}^{\infty }{C}_{m}\mathrm{sin}\,{k}_{m}x \qquad , \qquad {W}_{yb}=\sum _{m=1}^{\infty }{D}_{m}\mathrm{sin}\,{k}_{m}x.</math>                
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(12)
|}

The additional constraint torque shall be determined as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> {M}_{x0}=\sum_{n=1}^{\infty }{E}_{n}\mathrm{sin}\,{k}_{n}y \qquad ,\qquad {M}_{xa}=\sum _{n=1}^{\infty }{F}_{n}\mathrm{sin}\,{k}_{n}y,  </math>
|-
| style="text-align: center;" | <math>{M}_{y0}=\sum_{m=1}^{\infty }{G}_{m}\mathrm{sin}\,{k}_{m}x \qquad , \qquad {M}_{yb}=\sum _{m=1}^{\infty }{H}_{m}\mathrm{sin}\,{k}_{m}x.</math>                
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(13)
|}

The reciprocal theorem of work between the fundamental system and the real system can be utilized to derive:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{align}\displaystyle W\left( \xi ,\eta \right) = & \displaystyle\int_{0}^{a}\int_{0}^{b}\left[ \left( M{{w}_{mn}}^{2}-K\right) W\left( x,y\right) +G{\nabla }^{2}W\left( x,y\right) \right] {w}_{1}dxdy\\
&\displaystyle +\int_{0}^{a}{W}_{y0}{V}_{y0}dx+\int_{0}^{b}{W}_{x0}{V}_{x0}dy-\displaystyle\int_{0}^{a}{W}_{yb}{V}_{yb}dx- \displaystyle\int_{0}^{b}{W}_{xa}{V}_{xa}dy\\
& \displaystyle + \int_{0}^{a}{M}_{y0}\frac{\partial {w}_{1}}{\partial y}\Bigg| y=  0 dx\int_{0}^{b}{M}_{x0}\displaystyle \frac{\partial {w}_{1}}{\partial x}\Bigg| x=0dy\\
& \displaystyle -\int_{0}^{a}{M}_{yb}\frac{\partial {w}_{1}}{\partial y}\Bigg| y= b dx-\displaystyle \int_{0}^{b}{M}_{xa}\frac{\partial {w}_{1}}{\partial x }\Bigg| x=  a dy \end{align}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(14)
|}

Set:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  W\left( x,y\right) =\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{A}_{mn}\mathrm{sin}\,{k}_{m}x\mathrm{sin}\,{k}_{n}y </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(15)
|}

Substituting Eqs. (12), (13) and (15) into (14) yields:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math>{A}_{mn}=\frac{\frac{2}{a}{k}_{m}\left\{ D\left[ {{k}_{m}}^{2}+\left( 2-\mu \right) {{k}_{n}}^{2}\right] \left( {A}_{n}-{B}_{n}\mathrm{cos}\,m\pi \right) +\left( {E}_{n}-{F}_{n}\mathrm{cos}\,m\pi \right) \right\} }{D{\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) }^{2}-\left[ M{{w}_{mn}}^{2}-K-G\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) \right] }</math> 
|-
| style="text-align: center;" |<math>\qquad +\frac{\frac{2}{b}{k}_{n}\left\{ D\left[ {{k}_{n}}^{2}+\left( 2-\mu \right) {{k}_{m}}^{2}\right] \left( {C}_{m}-{D}_{m}\mathrm{cos}\,n\pi \right) +\left( {G}_{m}-{H}_{m}\mathrm{cos}\,n\pi \right) \right\} }{D{\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) }^{2}-\left[ M{{w}_{mn}}^{2}-K-G\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) \right] }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(16)
|}

Substituting <math display="inline">{A}_{mn}</math>  into <math display="inline">W\left( x,y\right)</math>   can obtain the mode function expression of the triangular series type:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>W\left( x,y\right) =\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{{A}_{n}\left[ {{k}_{m}}^{2}+\left( 2-\mu \right) {{k}_{n}}^{2}\right] \frac{2}{a}{k}_{m}D}{D{\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) }^{2}-\left[ M{{w}_{mn}}^{2}-K-G\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) \right] }\mathrm{sin}\,{k}_{m}x\mathrm{sin}\,{k}_{n}y</math> <math>\quad \quad \quad +\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{-{B}_{n}\mathrm{cos}\,m\pi \left[ {{k}_{m}}^{2}+\left( 2-\mu \right) {{k}_{n}}^{2}\right] \frac{2}{a}{k}_{m}D}{D{\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) }^{2}-\left[ M{{w}_{mn}}^{2}-K-G\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) \right] }\mathrm{sin}\,{k}_{m}x\mathrm{sin}\,{k}_{n}y</math> <math>\quad \quad \quad +\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{{C}_{m}\left[ {{k}_{m}}^{2}+\left( 2-\mu \right) {{k}_{n}}^{2}\right] \frac{2}{b}{k}_{n}D}{D{\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) }^{2}-\left[ M{{w}_{mn}}^{2}-K-G\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) \right] }\mathrm{sin}\,{k}_{m}x\mathrm{sin}\,{k}_{n}y</math> <math>\quad \quad \quad +\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{-{D}_{m}\mathrm{cos}\,n\pi \left[ {{k}_{n}}^{2}+\left( 2-\mu \right) {{k}_{m}}^{2}\right] \frac{2}{b}{k}_{n}D}{D{\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) }^{2}-\left[ M{{w}_{mn}}^{2}-K-G\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) \right] }\mathrm{sin}\,{k}_{m}x\mathrm{sin}\,{k}_{n}y</math> <math>\quad \quad \quad +\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{{E}_{n}\frac{2}{a}{k}_{m}}{D{\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) }^{2}-\left[ M{{w}_{mn}}^{2}-K-G\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) \right] }\mathrm{sin}\,{k}_{m}x\mathrm{sin}\,{k}_{n}y</math> <math>\quad \quad \quad +\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{-{F}_{n}\frac{2}{a}{k}_{m}\mathrm{cos}\,m\pi }{D{\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) }^{2}-\left[ M{{w}_{mn}}^{2}-K-G\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) \right] }\mathrm{sin}\,{k}_{m}x\mathrm{sin}\,{k}_{n}y</math> <math>\quad \quad \quad +\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{{G}_{m}\frac{2}{b}{k}_{n}}{D{\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) }^{2}-\left[ M{{w}_{mn}}^{2}-K-G\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) \right] }\mathrm{sin}\,{k}_{m}x\mathrm{sin}\,{k}_{n}y</math> <math>\quad \quad \quad +\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{-{H}_{m}\frac{2}{b}{k}_{n}\mathrm{cos}\,n\pi }{D{\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) }^{2}-\left[ M{{w}_{mn}}^{2}-K-G\left( {{k}_{m}}^{2}+{{k}_{n}}^{2}\right) \right] }\mathrm{sin}\,{k}_{m}x\mathrm{sin}\,{k}_{n}y</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(17)
|}


Set:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| <math>W\left( x,y\right) ={w}_{1}+{w}_{2}+{w}_{3}+{w}_{4}+{w}_{5}+{w}_{6}+{w}_{7}+{w}_{8}</math>
|}

The mode function of the triangular series can only satisfy the condition that the boundary deflection is zero while ensuring homogeneity in the corresponding boundary shear force and bending moment. Therefore, <math display="inline">{w}_{1}</math>, <math display="inline">{w}_{2}</math>, <math display="inline">{w}_{3}</math> and <math display="inline">{w}_{4}</math> undergo hyperbolic transformations in one of the following directions. Similarly, <math display="inline">{w}_{5}</math>, <math display="inline">{w}_{6}</math>, <math display="inline">{w}_{7}</math> and <math display="inline">{w}_{8}</math> are transformed accordingly.  We can get:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \begin{align} {w}_{1}= & \sum _{n=1}^{\infty }\frac{{A}_{n}}{\left( {{\alpha }_{n}}^{2}-{{\beta }_{n}}^{2}\right) }\left\{ \begin{matrix}\quad \left[ \frac{{{\alpha }_{n}}^{2}sh{\alpha }_{n}\left( a-x\right) }{sh{\alpha }_{n}a}-\frac{{{\beta }_{n}}^{2}sh{\beta }_{n}\left( a-x\right) }{sh{\beta }_{n}a}\right] \\\ +\left( 2-\mu \right) {{k}_{n}}^{2}\left[ -\frac{sh{\alpha }_{n}\left( a-x\right) }{sh{\alpha }_{n}a}+\frac{sh{\beta }_{n}\left( a-x\right) }{sh{\beta }_{n}a}\right] \end{matrix}\right\} \mathrm{sin}\,{k}_{n}y.\\
{w}_{2}=& \sum _{n=1}^{\infty }\frac{{B}_{n}}{\left( {{\alpha }_{n}}^{2}-{{\beta }_{n}}^{2}\right) }\left\{ \begin{matrix}\quad\left[ \frac{{{\alpha }_{n}}^{2}sh{\alpha }_{n}x}{sh{\alpha }_{n}a}-\frac{{{\beta }_{n}}^{2}sh{\beta }_{n}x}{sh{\beta }_{n}a}\right] \\\ +\left( 2-\mu \right) {{k}_{n}}^{2}\left[ -\frac{sh{\alpha }_{n}x}{sh{\alpha }_{n}a}+\frac{sh{\beta }_{n}x}{sh{\beta }_{n}a}\right] \end{matrix}\right\} \mathrm{sin}\,{k}_{n}y.\\
{w}_{3}= & \sum _{m=1}^{\infty }\frac{{C}_{m}}{\left( {{\alpha }_{m}}^{2}-{{\beta }_{m}}^{2}\right) }\left\{ \begin{matrix}\quad\left[ \frac{{{\alpha }_{m}}^{2}sh{\alpha }_{m}\left( b-y\right) }{sh{\alpha }_{m}b}-\frac{{{\beta }_{m}}^{2}sh{\beta }_{m}\left( b-y\right) }{sh{\beta }_{m}b}\right] \\\ +\left( 2-\mu \right) {{k}_{m}}^{2}\left[ -\frac{sh{\alpha }_{m}\left( b-y\right) }{sh{\alpha }_{m}b}+\frac{sh{\beta }_{m}\left( b-y\right) }{sh{\beta }_{m}b}\right] \end{matrix}\right\} \mathrm{sin}\,{k}_{m}x.\\
{w}_{4}= & \sum _{m=1}^{\infty }\frac{{D}_{m}}{\left( {{\alpha }_{m}}^{2}-{{\beta }_{m}}^{2}\right) }\left\{ \begin{matrix}\quad\left[ \frac{{{\alpha }_{m}}^{2}sh{\alpha }_{m}y}{sh{\alpha }_{m}b}-\frac{{{\beta }_{m}}^{2}sh{\beta }_{m}y}{sh{\beta }_{m}b}\right] \\\ +\left( 2-\mu \right) {{k}_{m}}^{2}\left[ -\frac{sh{\alpha}_{m}y}{sh{\alpha }_{m}b}+\frac{sh{\beta }_{m}y}{sh{\beta }_{m}b}\right] \end{matrix}\right\} \mathrm{sin}\,{k}_{m}x.\\
{w}_{5}= & \sum _{n=1}^{\infty }\frac{{E}_{n}}{D\left( {{\alpha }_{n}}^{2}-{{\beta }_{n}}^{2}\right) }\left[ - \frac{sh{\alpha }_{n}\left( a-x\right) }{sh{\alpha }_{n}a}+\frac{sh{\beta }_{n}\left( a-x\right) }{sh{\beta }_{n}a}\right] \mathrm{sin}\,{k}_{n}y.\\
{w}_{6}=& \sum _{n=1}^{\infty }\frac{{F}_{n}}{D\left( {{\alpha }_{n}}^{2}-{{\beta }_{n}}^{2}\right) }\left[ - \frac{sh{\alpha }_{n}x}{sh{\alpha }_{n}a}+\frac{sh{\beta }_{n}x}{sh{\beta }_{n}a}\right] \mathrm{sin}\,{k}_{n}y.\\
{w}_{7}= & \sum _{m=1}^{\infty }\frac{{G}_{m}}{D\left( {{\alpha }_{m}}^{2}-{{\beta }_{m}}^{2}\right) }\left[ - \frac{sh{\alpha }_{m}\left( b-y\right) }{sh{\alpha }_{m}b}+\frac{sh{\beta }_{m}\left( b-y\right) }{sh{\beta }_{m}b}\right] \mathrm{sin}\,{k}_{m}x.\\
{w}_{8}=& \sum _{m=1}^{\infty }\frac{{H}_{m}}{D\left( {{\alpha }_{m}}^{2}-{{\beta }_{m}}^{2}\right) }\left[ -\frac{sh{\alpha }_{m}y}{sh{\alpha }_{m}b}+\frac{sh{\beta }_{m}y}{sh{\beta }_{m}b}\right] \mathrm{sin}\,{k}_{m}x.\end{align}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(18)
|}

The coordinate factor after the aforementioned transformation not only satisfies the condition of zero boundary deflection but also fulfills equations (12) and (13). In this scenario, the actual system must adhere to boundary conditions (3). The substitution of <math display="inline">{B}_{n}</math>, <math display="inline">{C}_{m}</math>, <math display="inline">{D}_{m}</math>, <math display="inline">{E}_{n}</math>, <math display="inline">{F}_{n}</math>, <math display="inline">{G}_{m}</math> and <math display="inline">{H}_{m}</math> into Eq. (3) yields eight equations, the solutions of which provide the corresponding values. Substituting these values into Eq. (17) determines the coordinate factor  <math display="inline">W\left( x,y\right)</math>.

By substituting the coordinate factor <math display="inline">W\left( x,y\right)</math>  into Eq. (10), one can derive the expression for  <math display="inline">W\left( x,y,t\right)</math>  and obtain an analytical formula for the deflection of the pavement panel

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> W\left( x,y,t\right) =\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }\frac{W\left( x,y\right) }{M{w}_{mn}\int_{0}^{a}\displaystyle\int_{0}^{b}{W}^{2}\left( x,y\right) dxdy}\displaystyle\int_{0}^{t}{P}_{0}W\left( {v}_{x}\tau ,\eta \right) \mathrm{sin}\,{w}_{mn}\left( t-\tau \right) d\tau  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(19)
|}

==4. The identification of foundation modulus==

According to the analytical expression of deflection, the deflection value of the plate can be determined when the reaction modulus <math display="inline">K</math> of the foundation is fixed. Conversely, if the measured deflection at each measuring point is known, a fitting process can be conducted by establishing an objective function to compare and match the measured dynamic deflection with theoretical dynamic deflection. The program <span id='citeF-10'></span>[[#cite-10|[10]]] is implemented using Matlab for identifying the response modulus of the foundation. During parameter optimization, an objective function is formulated based on the least squares criterion:


<div style="text-align: right; direction: ltr; margin-left: 1em;">
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| <math>\mathrm{min}\,\epsilon (K)=\mathrm{min}\,\sqrt{\frac{1}{n}\sum _{i=1}^{n}({w}_{i}-\overline{{w}_{i}}{)}^{2}}</math>
|}</div>


Among them, <math display="inline">\mathrm{min}\,\epsilon (K)</math>,  for the objective function, <math display="inline">{w}_{i}</math> and  <math display="inline">\overline{{w}_{i}}</math> represent the theoretically calculated and measured deflection values of the measurement points, respectively, and <math display="inline">n</math> represents the number of measurement points on the pavement. Adjust <math display="inline">K</math> one by one to minimize the objective function, and the <math display="inline">K</math> value at this point is the identified foundation response modulus.

==5. Verification of the case==

The model and inversion recognition method presented in this paper were utilized to compute the pavement surface of an airport located in Nanjing, with geometric dimensions of a=6m, b=4m, elastic modulus <math display="inline">E=3.5\times 1{0}^{4\, }</math>MPa, Poisson ratio <math display="inline">\mu =</math><math>0.167</math>, density <math display="inline">\rho =2500\, \frac{{\rm kN}}{{\rm m}^{2}}</math>, thickness <math display="inline">h=</math><math>15.75</math>mm. Additionally, considered were foundation shear modulus <math display="inline">G=35.17 \frac{{\rm N}}{{\rm cm}^{3}}</math>, rigidity coefficients <math display="inline">k1=1.5706\times 1{0}^{4}\frac{{\rm kN}}{{\rm m}}</math>, <math display="inline">k2=1.5830\times 1{0}^{4}\frac{{\rm kN}}{{\rm m}}</math> as well as loading equipment for aircraft (<math display="inline">{P}_{0} =215000</math>N).

When the velocity reaches 2.59 m/s and <math display="inline">y=1.1</math>m, [[#tab-1|Table 1]] presents both the measured and theoretical deflection values, while [[#img-4|Figure 4]] illustrates the fitting curve.

<div class="center" style="font-size: 75%;">'''Table 1'''.  Inversion results of foundation modulus (velocity = <math display="inline">2.59</math> m/s and <math display="inline">y=1.1</math>m)</div>

<div id='tab-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! style="text-align: left;"| Coordinate
|  style="text-align: center;"|0.0 
|  style="text-align: center;"|0.5 
|  style="text-align: center;"|1.0 
|  style="text-align: center;"|1.5 
|  style="text-align: center;"|2.0 
|  style="text-align: center;"|2.5 
|  style="text-align: center;"|3.0 
|-
!  style="text-align: left;"|Measured deflection(mm)
|  style="text-align: center;"|0.395 
|  style="text-align: center;"|0.365 
|  style="text-align: center;"|0.335 
|  style="text-align: center;"|0.270 
|  style="text-align: center;"|0.195 
|  style="text-align: center;"|0.115 
|  style="text-align: center;"|0.050 
|-
!  style="text-align: left;"|Theoretical deflectio(mm)
|  style="text-align: center;"|0.402 
|  style="text-align: center;"|0.388 
|  style="text-align: center;"|0.349 
|  style="text-align: center;"|0.285 
|  style="text-align: center;"|0.201 
|  style="text-align: center;"|0.104 
|  style="text-align: center;"|0.073 
|-
!  style="text-align: left;"|Deviation value(mm)
|  style="text-align: center;"|-0.007
|  style="text-align: center;"|-0.023
|  style="text-align: center;"|-0.014
|  style="text-align: center;"|-0.015
|  style="text-align: center;"|-0.006
|  style="text-align: center;"|0.011
|  style="text-align: center;"|-0.023
|-
!  style="text-align: left;"|K(N/cm3)
|  colspan='7'  style="text-align: center;"|58.33
|}


<div id='img-4'></div>
{| class="wikitable" style="margin: 0em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Hui-min_340971489-image4.png|468px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 4'''. Comparison curve between theoretical deflection and measured deflection
|}


When the velocity reaches 2.79 m/s and <math display="inline">y=0.68</math>m, [[#tab-2|Table 2]] presents both the measured and theoretical deflection values, while [[#img-5|Figure 5]] illustrates the corresponding fitting curve.

<div class="center" style="font-size: 75%;">'''Table 2'''.  Inversion results of foundation modulus (velocity = <math display="inline">2.79</math> m/s and <math display="inline">y=0.68</math>m) </div>

<div id='tab-2'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! style="text-align: left;"| Coordinate
|  style="text-align: center;"|0.0 
|  style="text-align: center;"|0.5 
|  style="text-align: center;"|1.0 
|  style="text-align: center;"|1.5 
|  style="text-align: center;"|2.0 
|  style="text-align: center;"|2.5 
|  style="text-align: center;"|3.0 
|-
!  style="text-align: left;"|Measured deflection(mm)
|  style="text-align: center;"|0.5 
|  style="text-align: center;"|0.47 
|  style="text-align: center;"|0.397 
|  style="text-align: center;"|0.315 
|  style="text-align: center;"|0.225 
|  style="text-align: center;"|0.145 
|  style="text-align: center;"|0.071 
|-
!  style="text-align: left;"|Theoretical deflection(mm)
|  style="text-align: center;"|0.485 
|  style="text-align: center;"|0.469 
|  style="text-align: center;"|0.420 
|  style="text-align: center;"|0.343 
|  style="text-align: center;"|0.242 
|  style="text-align: center;"|0.127 
|  style="text-align: center;"|0.060 
|-
!  style="text-align: left;"|Deviation value(mm)
|  style="text-align: center;"|0.015
|  style="text-align: center;"|0.001
|  style="text-align: center;"|-0.023
|  style="text-align: center;"|-0.028
|  style="text-align: center;"|-0.017
|  style="text-align: center;"|0.018
|  style="text-align: center;"|0.011
|-
|  style="text-align: left;"|K(N/cm3)
|  colspan='7'  style="text-align: center;"|58.40
|}


<div id='img-5'></div>
{| class="wikitable" style="margin: 0em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Hui-min_340971489-image5.png|456px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 5'''. Comparison curve between theoretical deflection and measured deflection
|}


According to [[#tab-1|Tables 1]] and [[#tab-2|2]], as well as [[#img-4|Figures 4]] and [[#img-5|5]], it is evident that the theoretical deflection values align closely with the measured deflection values. This demonstrates that the dynamic model established in this paper accurately reflects real-world conditions. Furthermore, the stability and convergence of the parameter inversion method used have been proven.

==6. Conclusions==

(1) This study established a more accurate dynamic model for airport pavement systems, taking into account reasonable assumptions and simplifications of the foundation, loads, boundary conditions, and pavement panels. The model is a rigid thin plate pavement dynamic model that considers boundary shear and bending dual parameter foundation under moving loads. It accurately reflects the actual working state of rigid pavement and can be applied to the design of asphalt and overlay thickness on aircraft rigid pavement;

(2) The methods of variational method and reciprocal work theorem were used in this study to solve the differential equations of motion of the built dynamic system. An analytical solution for deflection was obtained, and an inversion calculation program based on the principle of least squares was compiled. Structural parameters of the rigid pavement panel were inverted and identified based on measured dynamic deflection. The stability, accuracy, and reliability of the identification method were verified through examples.

==Acknowledgments==

The authors would like to thank The Youth Innovation Team of Shaanxi Universities(Key Technology Innovation Team for Urban Rail Transit Track Bed);Youth fund project at Xi’an Jiaotong Engineering Institute(2023KY-39) and The Scientific Research Program Funded by Education Department of Shaanxi Provincial Government (Program No.23JK0532).Hui-min Cao and Han-jun Zhao contributed equally to this work.

==References==
<div class="auto" style="text-align: left;width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">

<div id="cite-1"></div>
[[#citeF-1|[1]]]  Zhang W., Chen X.,  Zuo X. Research on airport rigid pavement parameter identification on dynamic loading (in Chinese). Computer Simulation, 29(09):73-76+202, 2012.

<div id="cite-2"></div>
[[#citeF-2|[2]]] Chen X., Liu H. Dynamic identification simulation of foundation response modulus of airport (in Chinese). Computer Simulation, 26(07):54-57, 2009. 

<div id="cite-3"></div>
[[#citeF-3|[3]]]  Xing Y., Liu H. Foundation parameter identification research of airport rigid pavement (in Chinese). Central South Highway Engineering, (04):52-54, 2006.

<div id="cite-4"></div>
[[#citeF-4|[4]]]  Patil V., Sawant V., Deb K. Finite element analysis of rigid pavement on a nonlinear two parameter foundation model. International Journal of Geotechnical Engineering, 6(3):275-286, 2012. DOI: 10.3328/IJGE.2012.06.03.274-286

<div id="cite-5"></div>
[[#citeF-5|[5]]] Kumar Y., Trivedi A., Shukla S.K. Deflections governed by the cyclic strength of rigid pavement subjected to structural vibration due to high-velocity moving loads.  Journal of Vibration Engineering & Technologies, 12:3543–3562, 2024. 

<div id="cite-6"></div>
[[#citeF-6|[6]]]  Taheri M.R., Zaman M.M . Effects of a moving aircraft and temperature differential on response of rigid pavements. Computers & Structures, 57(3):503-511, 1995. DOI:10.1016/0045-7949(94)00625-D

<div id="cite-7"></div>
[[#citeF-7|[7]]]  Worku A., Lulseged A. Kinematic pile-soil interaction using a rigorous two-parameter foundation model. Soil Dynamics and Earthquake Engineering, 165, 107701, 2023.

<div id="cite-8"></div>
[[#citeF-8|[8]]]  Li C.J., Zhang Y., Jia Y.M., Chen J. The polygonal scaled boundary thin plate element based on the discrete Kirchhoff theory. Computers & Mathematics with Applications, 97:223-236, 2021.

<div id="cite-9"></div>
[[#citeF-9|[9]]] Alisjahbana S.W., Wangsadinata W. Dynamic response of rigid roadway pavement under moving traffic loads. The 8th International Structural Engineering and Construction Conference, Sydney, Australia, 2015.

<div id="cite-10"></div>
[[#citeF-10|[10]]]  Heitmann S., Aburn M.J., Breakspear M. The brain dynamics toolbox for Matlab. Neurocomputing, 315:82-88, 2018. DOI:10.1016/j.neucom.2018.06.026