CN112765825B - Determination of Subgrade Surface Displacement Response Considering Nonlinear Modulus Distribution - Google Patents

Abstract

The invention discloses a roadbed surface displacement response determining method considering modulus nonlinear distribution, which comprises the following steps: obtaining roadbed material parameters, applying dynamic load on the roadbed surface, and constructing a roadbed modulus function which is distributed in a vertical nonlinear manner; calculating a surface displacement response solution of the roadbed in the Laplace-Hankel domain, wherein the surface displacement response solution simultaneously considers the transverse isotropy and the nonlinear distribution of the roadbed modulus along the vertical direction; and calculating the surface displacement response of the current roadbed model in a time domain range based on a displacement response solution of the current roadbed mechanics model in a Laplace-Hankel domain. The invention fully reflects the non-uniform distribution of the modulus of the actual roadbed along the depth direction, simultaneously considers the transverse isotropy of the roadbed, improves the accuracy of the roadbed surface displacement response, has high precision and high speed, and can obtain the distribution condition of the modulus of the roadbed along the depth direction by combining with the roadbed nondestructive detection technology.

Description

Determination of Subgrade Surface Displacement Response Considering Nonlinear Modulus Distribution

technical field

The invention belongs to the technical field of road engineering, relates to a method for determining the displacement response of a subgrade surface considering the nonlinear distribution of modulus, and in particular relates to a method for determining the displacement response of a subgrade surface considering the transverse isotropy and the nonlinear distribution of the modulus and its application , equipment and storage media.

Background technique

With the development of non-destructive testing technology, some soil-based rapid non-destructive testing instruments have appeared one after another, among which the portable drop weight deflectometer (PFWD) is the most widely used. PFWD performs the inversion calculation of subgrade parameters by applying impact load on the top surface of the subgrade, combined with a certain mechanical model and parameter adjustment algorithm. However, the current mainstream subgrade parameter inversion calculation is still based on the isotropic linear elastic system, which cannot reflect the actual characteristics of the subgrade. At present, a large number of research results show that the subgrade has strong stress-dependent characteristics, and as the depth increases, the subgrade modulus generally presents an upward trend with the depth. In addition, the transverse isotropic characteristics of the subgrade caused by natural settlement cannot be ignored in the actual calculation.

Some existing technologies are aimed at the pavement structure, the material modulus in each layer is considered to be uniform and does not change with space, but the modulus of the roadbed varies in the vertical direction, and the modulus has nonlinear distribution characteristics in the vertical direction. Existing methods are difficult to accurately determine the displacement response of the subgrade surface. Part of the existing technology is based on the consideration of the homogeneous modulus of the subgrade, iterative calculation is performed through the nonlinearity of the strain and modulus of the structural layer, which is essentially a mathematical iterative method, and the mechanical model does not consider this kind of modulus. Due to the nonlinear distribution characteristics, the iterative method can only approach the real results, the calculation accuracy is low, and repeated iterations are required, and the calculation is slow.

SUMMARY OF THE INVENTION

In order to solve the above problems, the present invention provides a method for determining the displacement response of the subgrade surface considering the nonlinear distribution of modulus, which fully reflects the non-uniform distribution of the modulus of the actual subgrade along the depth direction, and at the same time considers the transverse isotropy of the subgrade, which improves the performance of the subgrade. To determine the accuracy of the displacement response of the subgrade surface, this method has the characteristics of high precision and high speed. Combined with the subgrade non-destructive testing technology, the distribution of the subgrade modulus along the depth direction can be obtained, which solves the problems existing in the prior art.

The technical scheme adopted in the present invention is a method for determining the displacement response of the subgrade surface considering the nonlinear distribution of modulus, which is specifically carried out according to the following steps:

Step S1, obtaining subgrade material parameters, the subgrade material parameters include: the vertical modulus E _v0 of the subgrade top surface, the subgrade modulus E _v∞ at the infinite depth of the subgrade, the subgrade modulus non-uniformity coefficient α, the subgrade modulus ratio _ne , horizontal Poisson’s ratio μ _h , vertical Poisson’s ratio μ _v , subgrade density ρ and subgrade thickness H;

ξ和s分别表示Hankel变换以及Laplace积分变换系数；Step S2, apply a dynamic load p(r, t) to the subgrade surface, and establish a cylindrical coordinate system with the center of the dynamic load as the coordinate center, where r represents the radial coordinate, z represents the vertical coordinate, and φ represents the circumferential coordinate. The nonlinearly distributed subgrade modulus function E _v (z)=E _v0 +(E _v∞ -E _v0 )(1-e ^-αz ); the calculation also considers the transverse isotropy and the vertical nonlinearity of the subgrade modulus Solution of Surface Displacement Response of Distributed Subgrade in Laplace-Hankel Domain

ξ and s represent Hankel transform and Laplace integral transform coefficients, respectively;

计算当前路基模型在时域范围内的表面位移响应。Step S3, the displacement response solution based on the current subgrade mechanical model in the Laplace-Hankel domain

Calculate the surface displacement response of the current subgrade model in the time domain.

The above method is used to obtain the vertical modulus distribution state of the subgrade on site.

An electronic device comprising:

memory for storing instructions executable by the processor; and

A processor for executing the instructions to implement the above method.

A computer-readable storage medium is used for storing a computer program, and when the computer program is executed, the above-mentioned method can be implemented.

The beneficial effects of the present invention are:

The embodiment of the present invention describes the modulus distribution of the subgrade along the depth direction through the constructed E _v (z) function, and at the same time introduces a transverse isotropic physical equation to derive a variable coefficient partial differential equation system considering the characteristics above the subgrade. And through Hankel-Laplace integral transformation and Frobenius method, an analytical solution and a corresponding determination method for the dynamic response of the subgrade considering the transverse isotropy of the subgrade and the nonlinear distribution of the modulus along the vertical direction are proposed. In the calculation of the dynamic response of the subgrade, the transverse isotropy of the subgrade and the nonlinear distribution characteristics of the modulus along the depth direction are fully considered, which can more truly reflect the actual state of the subgrade than the traditional elastic system. Compared with traditional methods such as finite element method, it has the characteristics of high precision and fast speed.

The determination method according to the embodiment of the present invention can obtain the modulus distribution describing the subgrade, so as to evaluate the subgrade stiffness, which is more accurate and faster than the traditional method; combined with the nondestructive testing of the subgrade, the vertical modulus curve of the subgrade can be obtained, The vertical modulus change curve of long-term monitoring can evaluate subgrade humidification, and provide a more reliable theoretical basis for subgrade design, testing, evaluation, and research.

Description of drawings

In order to explain the embodiments of the present invention or the technical solutions in the prior art more clearly, the following briefly introduces the accompanying drawings that need to be used in the description of the embodiments or the prior art. Obviously, the accompanying drawings in the following description are only These are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can also be obtained according to these drawings without creative efforts.

FIG. 1 is a schematic diagram of a subgrade mechanics system with isotropy in transverse view and nonlinear distribution of modulus in the vertical direction under the action of a rigid bearing plate according to an embodiment of the present invention.

FIG. 2 is a flowchart of a method for determining a dynamic response according to an embodiment of the present invention.

FIG. 3 is a diagram of an axisymmetric model of ABAQUS in an embodiment of the present invention.

4a-4b are comparison diagrams between the method for determining the dynamic response according to the embodiment of the present invention and the method using the finite element model ABAQUS.

5 is a comparison diagram of the numerical simulation calculation result obtained by the humidity-stress dependent finite element method and the actual value of the modulus distribution along the depth direction of the subgrade structure obtained by fitting according to the embodiment of the present invention.

Detailed ways

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

An embodiment of the present invention provides a method for determining the displacement response of a subgrade surface considering the nonlinear distribution of modulus, as shown in FIG. 2 , which is specifically performed according to the following steps:

Step S1, obtaining subgrade material parameters, the subgrade material parameters include: the vertical modulus E _v0 of the subgrade top surface, the subgrade modulus E _v∞ at the infinite depth of the subgrade, the subgrade modulus non-uniformity coefficient α, the subgrade modulus ratio _ne , horizontal Poisson’s ratio μ _h , vertical Poisson’s ratio μ _v , subgrade density ρ and subgrade thickness H; the subgrade modulus non-uniformity coefficient α characterizes the change speed of subgrade modulus with depth direction; the subgrade modulus ratio n _e is uniformly distributed in the entire subgrade space; the thickness H can be a finite value or an infinite value.

Step S2, as shown in Figure 1, apply a dynamic load p(r, t) to the subgrade surface, and establish a cylindrical coordinate system with the center of the dynamic load as the coordinate center, where r represents the radial coordinate, z represents the vertical coordinate, and φ represents the circumferential direction. Coordinate, construct the subgrade modulus function E _v (z)=E _v0 +(E _v∞ -E _v0 )(1-e ^-αz ) distributed along the vertical nonlinearity, R in Figure 1 represents the radius of the rigid bearing plate; Calculate the solution of the surface displacement response of the subgrade in the Laplace-Hankel domain considering the transverse isotropy and the vertical nonlinear distribution of the modulus

的计算为：Subgrade Surface Displacement Response Solution in Laplace-Hankel Domain of Subgrade Mechanical Model Considering Transverse Isotropy and Vertical Nonlinear Modulus Distribution

is calculated as:

In Equation (1), A ₁ , A ₂ , A ₃ , and A ₄ are boundary condition parameters, which are obtained by solving the following matrix equation [Equation (2)]. Among them, if the subgrade thickness H→+∞, there is A ₃ =A ₄ =0, then the formula (2) can be degenerated into a two-dimensional matrix equation;

为动荷载p(r,t)进行Hankel-Laplace变换后的结果，

ξ以及s分别为Hankel积分变换以及Laplace积分变换系数；E_v0表示路基顶面模量；In formula (2)

The result of Hankel-Laplace transformation for the dynamic load p(r,t),

ξ and s are the Hankel integral transformation and Laplace integral transformation coefficients respectively; E _v0 represents the subgrade top surface modulus;

In the formula, α is the non-uniform coefficient of the subgrade modulus; e is a constant, Ξ is obtained by calculating the subgrade top surface modulus E _v0 and the subgrade infinite depth modulus E _v∞ , calculated according to formula (7); Ξ represents the subgrade top The dimensionless coefficient of the relative size of the surface and bottom modulus, the larger the value, the greater the difference between the upper and lower modulus of the subgrade;

C₁₃＝kn_eμ_v(1+μ_h)，

C₄₄＝0.5(1+μ_v)^-1；

k表示中间变量，用于计算横观各向同性的各项参数，简化表达式；n_e为路基模量比；μ_v，μ_h为路基垂直和水平方向上的泊松比，n是级数序列的序号。C ₁₁ , C ₁₂ , C ₁₃ , C ₃₃ , and C ₄₄ are parameters that characterize the transverse isotropic characteristics of subgrade materials in axisymmetric problems, C ₁₁ = _kne (1- _ne μ _v ),

C ₁₃ = _kne μ _v (1+μ _h ),

C ₄₄ =0.5(1+μ _v ) ⁻¹ ;

k represents the intermediate variable, which is used to calculate the parameters of transverse isotropy, and the expression is simplified; n _e is the subgrade modulus ratio; μ _v , μ _h are the Poisson’s ratio in the vertical and horizontal directions of the subgrade, and n is the grade Number sequence number.

On the basis of the subgrade mechanics model shown in Fig. 1, combined with the dynamic balance equations of the axisymmetric problem [Equations (26)-(27)], the physical equations that characterize the axisymmetric transverse isotropy, and the geometric equations [Equations (28)] ~(31)], the partial differential equation system shown in equations (32)-(33) can be obtained.

Among them, ur = ur ( _r , _z , t), uz = uz ( _r , _z , t) are the displacement components along the r direction and the z direction respectively; S2 obtains the displacement in the Hankel-Laplace domain Analytical solution, S3 obtains the displacement response results in the time domain through inverse integral transformation; ρ is the density of the roadbed; t is the time; σ _r =σ _r (r,z,t),σ _z =σ _z (r,z ,t), σ _φ =σ _φ (r,z,t) are the stress components in the r direction, z direction and φ direction of the cylindrical coordinate system respectively; τ _zt =τ _zt (r,z,t) is the zr plane shear stress on.

定义关于x的零阶Hankel变换：

定义关于x的一阶Hankel变换：

ξ和s分别表示Hankel变换以及Laplace积分变换系数。In order to solve the system of partial differential equations shown in equations (32) to (33), the Laplace transform of t in equations (32) to (33) is firstly performed; The first-order Hankel transform; the zero-order Hankel transform about r is performed on the formula (33) after Laplace transform; the variable coefficient ordinary differential equation system is obtained as formulas (34)-(35), and the independent variable z is omitted in the formula. where, define the Laplace transform with respect to x:

Define the zeroth-order Hankel transform with respect to x:

Define the first-order Hankel transform with respect to x:

ξ and s represent Hankel transform and Laplace integral transform coefficients, respectively.

表示对u_r进行Hankel-Laplace一阶变换后的结果，

表示对u_z进行Hankel-Laplace零阶变换后的结果。

represents the result of Hankel- _Laplace first-order transformation of ur,

Indicates the result of Hankel-Laplace zero-order transformation of u _z .

To solve the system of ordinary differential equations with variable coefficients shown in equations (34) to (35), the following variables are introduced:

ψ=Ξe ^-αz (36)

Ξ is calculated according to formula (7), and in the subsequent derivation, the variable ψ is used to replace the independent variable z in the formula, and the expression of E _v (ψ) can be obtained as shown in formula (38). Using Ψ to perform variable substitution for z, E _v (ψ) has the same meaning as E _v (z), except that the independent variable is replaced, and its value remains unchanged. The independent variable is z in (34)~(35) and ψ in (40)~(43); the purpose is to eliminate the exponential term that exists in E _v (z) and change it to be as shown in (39) Algebraic form; only variable substitution is done, the expression itself is not changed.

At the same time, applying the derivative theory, the relationship of formula (39) can be obtained; where f(z) represents any function with z as an independent variable:

Substituting equations (38)-(39) into equations (34)-(35), the system of variable-coefficient ordinary differential equations with ψ independent variable as shown in (40)-(43) can be obtained.

Among them, equations (40) to (41) correspond to the situation of E _v0 <E _v∞ , and equations (42) to (43) correspond to the situation of E _v0 >E _v∞ . The independent variable ψ is omitted in the formula; among them, δ is an intermediate calculation parameter, δ=E _v0 /E _v∞ .

表示

的一阶导，

表示

的二阶导；

表示

的一阶导，

表示

的二阶导；Based on the infinite series theory, the infinite series is used for approximation. It is assumed that the ordinary differential equations shown in equations (40) to (43) have displacement solutions as shown in equations (44) to (45); where

express

the first derivative of ,

express

the second derivative of ;

express

the first derivative of ,

express

the second derivative of ;

表示z方向位移分量的无穷级数中对应m_i的第n项级数的系数，

表示r方向位移分量的无穷级数中对应m_i的第n项级数的系数；a_n为

的简写，b_n为

的简写；将如式(44)～(45)所示的级数形式位移解

带入到式(40)～(43)中，可获得如式(46)～(49)所示的等式关系。式(46)～(47)对应E_v0<E_v∞的情况，式(48)～(49)对应E_v0>E_v∞的情况。式中，θ为中间计算参数，θ＝ρs²/E_v∞；In the formula, m is the parameter, and n is the serial number of the series;

represents the coefficient of the n-th series corresponding to m _i in the infinite series of displacement components in the z direction,

Represents the coefficient of the nth series corresponding to m _i in the infinite series of displacement components in the r direction; a _n is

shorthand for , b _n is

The abbreviation of

Substituting into the equations (40) to (43), the equational relationships shown in the equations (46) to (49) can be obtained. Equations (46) to (47) correspond to the case of E _v0 <E _v∞ , and equations (48) to (49) correspond to the case of E _v0 >E _v∞ . In the formula, θ is an intermediate calculation parameter, θ=ρs ² /E _v∞ ;

Let n=0, and compare the two sides of equations (46) to (49), the following relationship can be obtained:

[α ² C ₄₄ m ² -C ₁₁ ξ ² -θ]a ₀ +αξ(C ₁₃ +C ₄₄ )mb ₀ =0 (50)

αξ(C ₁₃ +C ₄₄ )ma ₀ -[α ² C ₃₃ m ² -C ₄₄ ξ ² -θ]b ₀ =0 (51)

If equations (50) to (51) are established, m must satisfy equation (52).

In the complex number domain, the equation for m [Eq. (52)] must have four solutions, in order:

Θ₃＝(C₄₄ξ²+θ)(C₁₁ξ²+θ)；其中θ为中间计算参数，θ＝ρs²/E_v∞；ρ为路基密度，θ表达式中的s为Hankel积分变换以及Laplace积分变换系数，积分变换后，相当于时间t就变成了s，可当成一个变量看待。In formula (8)～(11), m _i is an intermediate calculation parameter, Θ ₁ , Θ ₂ , Θ ₃ are respectively intermediate calculation parameters; Θ ₁ =α ² C ₃₃ C ₄₄ ;

Θ ₃ =(C ₄₄ ξ ² +θ)(C ₁₁ ξ ² +θ); where θ is the intermediate calculation parameter, θ=ρs ² /E _v∞ ; ρ is the subgrade density, and s in the θ expression is the Hankel integral Transformation and the Laplace integral transformation coefficient, after the integral transformation, the equivalent time t becomes s, which can be regarded as a variable.

表达式。

分别表示无穷级数中对应m_i的第n项级数的系数，为特定的系数序列；m has four solutions. According to formula (44), it can be seen that the final solution should be in the form of a combination of four infinite series, in order to distinguish the coefficients an in the above four _infinite series, so add the subscript to distinguish, formula The superscripts (1), (2), (3), and (4) in (12) correspond to the four solutions of m respectively. Substituting equation (12) into equation (50), the equation (13) can be obtained

expression.

respectively represent the coefficients of the nth series corresponding to m _i in the infinite series, which are specific coefficient sequences;

的递推公式，见式(14)、(15)。According to (46)-(49), the corresponding series coefficients are equal, that is, the series of the same power should be correspondingly equal on both sides of the equation, thus we can obtain

The recurrence formula of , see equations (14) and (15).

Calculated according to the recurrence relation shown in equations (12) to (15)

For the expressions after parameter expansion in equations (14) and (15), see (16) to (21):

为中间计算参数，按式(16)～(21)计算。其中，

的计算随着E_v0与E_v∞的相对关系计算有所区别，这里分为两种情况依次进行考虑，当E_v0<E_v∞时，按式(20)计算，当E_v0>E_v∞时，按式(21)计算，式(21)中，δ为中间计算参数，δ＝E_v0/E_v∞。In formulas (12) to (15),

is an intermediate calculation parameter, calculated according to formulas (16) to (21). in,

The calculation is different with the relative relationship between E _v0 and E _v∞ . There are two cases to consider in turn. When E _v0 <E _v∞ , calculate according to formula (20), when E _v0 >E _v When _∞ , calculate according to formula (21), in formula (21), δ is an intermediate calculation parameter, δ=E _v0 /E _v∞ .

引入中间参数

根据(46)～(49)中对应级数系数相等就可以获得递推公式，因为n＝0的结果已经知道了，通过这个可以获得n＝1的结果，依次递推计算，n＝2，n＝3，…都可以计算出来。In formula (20)

Introduce intermediate parameters

According to the corresponding series coefficients in (46) to (49), the recursive formula can be obtained, because the result of n=0 is already known, and the result of n=1 can be obtained through this, and the recursive calculation in turn, n=2, n=3, ... can be calculated.

According to the theory of ordinary differential solution, let the solution of the equation be in series form, and m has four corresponding results. Therefore, there are four infinite series solutions to the equation. The final general solution formula can be expressed by the linear combination of the above four infinite series solutions. , that is, the final general displacement solution is shown in Equation (53) and Equation (1):

There are four undetermined coefficients A ₁ , A ₂ , A ₃ , and A ₄ in the general solution. To solve for its value, it is necessary to bring in certain boundary conditions. According to Figure 1, the boundary conditions shown in equations (54) to (57) are introduced; where p(r, t) is the expression of the load on the top surface of the subgrade.

u _r | _{z = H} = 0 (56)

u _z | _{z = H} = 0 (57)

The Laplace transform of t and the 0-order Hankel transform of r are performed on equations (54) and (57); the Laplace transform of t and the first-order Hankel transform of r are performed on equations (55) and (56). On this basis, the transformation relations shown in equations (7) and (38) are introduced, and the boundary conditions shown in (58) to (61) can be obtained.

The general solution expressions of equations (53) and (1) are put into the boundary conditions shown in equations (58) to (61), and the fourth-order matrix equation shown in equation (2) can be obtained. 2), the undetermined coefficients A ₁ , A ₂ , A ₃ , and A ₄ can be obtained from the matrix equation shown in There must be A ₃ =A ₄ =0, then the matrix equation shown in equation (2) can be further degenerated into a second-order matrix equation. The determined A ₁ , A ₂ , A ₃ , and A ₄ are then substituted into equations (53) and (1), and the analytical solution results in the Hankel-Laplace domain can be obtained. Further, perform inverse Hankel transform and inverse Laplace transform on the result to obtain its time domain solution. The calculation process is shown in Figure 2.

进一步根据式(2)～(6)建立四阶矩阵方程，在此基础上求解系数A₁,A₂,A₃,A₄。最终按(1)式计算在Hankel-Laplace域内的路基表面动力响应解

In the actual calculation, firstly, according to the vertical modulus E _v0 of the top surface of the subgrade, the modulus E _v∞ at the infinite depth of the subgrade, the non-uniformity coefficient α of the subgrade modulus, the subgrade modulus ratio _ne , the horizontal Poisson's ratio μ _h , Parameters such as vertical Poisson’s ratio μ _v , subgrade density ρ and subgrade thickness H are calculated according to formulas (8) to (11), m ₁ , m ₂ , m ₃ , and m ₄ ; further coefficients are calculated according to formulas (12) to (21). series

Further, a fourth-order matrix equation is established according to equations (2) to (6), and the coefficients A ₁ , A ₂ , A ₃ , and A ₄ are solved on this basis. Finally, the dynamic response solution of the subgrade surface in the Hankel-Laplace domain is calculated according to the formula (1).

计算当前路基模型在时域范围内的表面位移响应。Step S3, the displacement response solution based on the current subgrade mechanical model in the Laplace-Hankel domain

Calculate the surface displacement response of the current subgrade model in the time domain.

进行Hankel逆变换，然后通过式(22)对Hankel逆变换的结果进行Laplace逆变换，求得当前路基力学体系在时域范围内的表面位移响应：The displacement response of the current subgrade mechanical model in the Laplace-Hankel domain is solved by formula (23).

Perform the inverse Hankel transform, and then perform the inverse Laplace transform on the result of the Hankel inverse transformation by formula (22), to obtain the surface displacement response of the current subgrade mechanical system in the time domain:

为

则f(r,0,t_l)为u_z(r,0,t_l)；u_z(r,0,t_l)为当前路基力学模型在时域范围内的表面位移响应，其对应的时间步为t_l，t_l＝l·T/N_a，T为求解总时长，N_a是求解的总的时间增量步数，l是增量步序列，l≤N；r是计算点位距离荷载中心的水平距离；Δt是时间增量步长；式(22)表示采用Durbin法进行Lpalace逆变换，L和γ是Durbin法计算参数，L·N_a＝50～5000，γ·T＝5～10；j是单位虚数，j²＝-1；对于Durbin法，最后1/4的时间求解会产生较大的系统误差，通常采用延长求解总时长T来解决，本发明通过测试，取2倍时间长度可以满足精度要求(对Hankel-Laplace域内的解进行逆变换时，需要给定计算的时长，比如需要计算60ms动力响应，由于Durbin法存在一定的缺陷，因此在实际计算时，计算时长参数一般取2～3倍，本发明实施例取2倍时间长度，即计算时长参数T取120ms，但最后的结果仅取前60ms)。χ是数值积分的上限，A_k是20节点高斯插值的第k个积分节点的权重；x_ik是高斯积分对应的i个积分区段的第k个积分节点值，x_ik＝(i-1)ΔL+x_k，ΔL是高斯积分区段划分长度，x_k为20节点高斯积分的第k个积分节点值；N₀是高斯积分的分段总数，N₀＝ceil(χ/ΔL)，通过对计算(χ/ΔL)并向上取整获得；

表示任意函数f(·)在Hankel-Laplace变换后的结果；

表示任意函数f(·)在Laplace变换后的结果；Re(·)为取复数实部运算；J₀(x_ikr)表示自变量为(x_ikr)时零阶贝塞尔函数对应的函数结果。式(1)、式(22)为Durbin法(一阶自相关检验方法)进行Laplace逆变换，式(23)为采用高斯插值积分进行Hankel逆变换，以上计算过程均采用编程方式实现，这里仅介绍数值确定方法。所述数值确定方法均为编程实现，不限定编程语言，本发明实施例计算的数据均通过MATLAB 2016a计算得到。in,

for

Then f(r, 0, t _l ) is u _z (r, 0, t _l ); u _z (r, 0, t _l ) is the surface displacement response of the current subgrade mechanical model in the time domain, and its corresponding The time step is t _l , t _l =l·T/N _a , T is the total solution time, N _a is the total number of time increment steps for the solution, l is the sequence of incremental steps, l≤N; r is the calculation point The horizontal distance of the bit from the load center; Δt is the time increment step size; Equation (22) indicates that the Durbin method is used to perform the Lpalace inverse transformation, L and γ are the parameters calculated by the Durbin method, L·N _a =50～5000, γ·T =5～10; j is a unit imaginary number, j ² =-1; for the Durbin method, the last 1/4 time solution will produce a large systematic error, which is usually solved by extending the total solution time T. The present invention passes the test, Taking twice the time length can meet the accuracy requirements (when performing inverse transformation of the solution in the Hankel-Laplace domain, the calculation time needs to be given, for example, 60ms dynamic response needs to be calculated. Due to the certain defects of the Durbin method, in actual calculation, The calculation duration parameter generally takes 2 to 3 times, and the embodiment of the present invention takes 2 times the time length, that is, the calculation duration parameter T takes 120ms, but the final result only takes the first 60ms). χ is the upper limit of the numerical integration, A _k is the weight of the k-th integration node of the 20-node Gaussian interpolation; x _ik is the k-th integration node value of the i integration section corresponding to the Gaussian integration, x _ik = (i-1 )ΔL+x _k , ΔL is the division length of the Gaussian integral segment, x _k is the value of the kth integral node of the 20-node Gaussian integral; N ₀ is the total number of segments of the Gaussian integral, N ₀ =ceil(χ/ΔL), Obtained by calculating (χ/ΔL) and rounding up;

Represents the result of the Hankel-Laplace transformation of an arbitrary function f( );

Represents the result of any function f(·) after Laplace transformation; Re(·) is the operation of taking the real part of complex numbers; J ₀ (x _ik r) represents the zero-order Bessel function corresponding to the independent variable (x _ik r). function result. Equation (1) and Equation (22) are the inverse Laplace transform using the Durbin method (first-order autocorrelation test method), and Equation (23) is the inverse Hankel transform using Gaussian interpolation. Introduce the numerical determination method. The numerical determination methods are all implemented by programming, and the programming language is not limited. The data calculated in the embodiments of the present invention are all calculated by MATLAB 2016a.

It has been verified that when the upper limit n of the infinite series sequence in equation (1) is 100, the upper limit of the Gaussian integral integral χ in equation (23), and the integral interval length ΔL are calculated according to formulas (24) to (25), the result can be converged to exact solution.

Effect verification:

The finite element software ABAQUS was used for comparative calculation. The parameters of the subgrade are shown in Table 1. In ABAQUS, the self-defined subprogram UMAT is used to realize the nonlinear distribution of the subgrade along the depth direction. p(r,t) is the load corresponding to a rigid bearing plate with a diameter of 30cm. The load is in the form of a half-sine wave, the peak value is 100kPa, and the action time is 20ms; the calculation time is set to 40ms, and the time step increment is 1ms. . When applying the present invention for calculation, according to the parameter recommended values of formulas (24) to (25), the calculation parameters are set as: n=100, χ=300, ΔL=10.

Table 1 Subgrade structure calculation parameter table

name unit type Modulus/MPa n_e μ_v(=μ_h) ρ/kg·m^-3 H/m carrier plate CAX4 2.1×10⁵ 1.0 0.25 7500 - roadbed CAX4/CINAX4 100+200(1-e^-2z) 1.0 0.35 1816 1.0/2.0

As shown in Figure 3, the entire ABAQUS finite element model is modeled in an axisymmetric manner. In order to simulate the infinity of the layered system in the horizontal direction, the infinite element CINAX4 element is used for the boundary part, and the CAX4 element is used for the finite size part. The calculation results are shown in Figure 4. In Figure 4, MATLAB is the calculation result obtained by using the method of the embodiment of the present invention, and ABAQUS is the calculation result using the finite element model. It can be seen that the calculation results of the two for the same mechanical model are basically the same, so the method provided by the embodiment of the present invention has better calculation accuracy, and can satisfy the subgrade considering the transverse isotropy and the nonlinear distribution of the modulus along the vertical direction. Requirements for dynamic response calculations. In the calculation, using the Intel i7 8750H/8G notebook test, it takes 5 minutes for ABAQUS to complete the calculation, while the method of the present invention only takes 2 minutes, which greatly reduces the calculation cost.

When calculating the uneven distribution of roadbed according to the previous method, it is necessary to obtain local meteorological data. The meteorological data mainly includes data such as wind speed, temperature, and precipitation. The calculation software is used to calculate the uneven distribution of the humidity field of the roadbed. On this basis, the non-uniform modulus field of the subgrade is calculated considering the stress state of the subgrade, the data of the non-uniform humidity field, and the coupling equation of the subgrade soil dynamic rebound modulus and wet force. The calculation consumption is large and it is not suitable for promotion.

The present invention adopts E _v (z)=E _v0 +(E _v∞ -E _v0 )(1-e ^-αz ) to approximately describe the vertical modulus distribution of the roadbed. In the formula, E _v0 , E _v∞ and α are approximate matching parameters. In order to verify the rationality of the formula, the meteorological data of Nanchang City in the past 16 years are used, combined with the soil-modulus wet-mechanical coupling equation to analyze the depth direction of the center of the roadway in the initial state, 1-year operation, 2-year operation and equilibrium humidity state. The modulus of the subgrade was calculated, and E _v (z) combined with the least squares method was used for fitting. It can be seen from Figure 5 that the fitting effect is good, indicating that E _v (z) can ideally characterize the operation process of the subgrade. The modulus distribution of the subgrade in the depth direction.

Regardless of the operating time of the subgrade, the E _v (z) formula can be used to fit this, and the fitting effect is good, indicating that E _v (z) can ideally characterize the subgrade along the depth direction of the subgrade during operation. Modulus distribution law.

On the basis of verifying the rationality of E _v (z), the embodiment of the present invention further proposes a method for determining the displacement response of the subgrade surface based on the nonlinear distribution of the subgrade, and the determination method can be used to inversely calculate the modulus distribution of the subgrade. ,Specifically:

Measure the displacement response on site, record the displacement response and load data on site; randomly assume a series of E _v (z) parameters E _v0 , E _v∞ , α and the subgrade modulus ratio _ne within a certain range, where E _v0 The range of E v∞ is 20MPa～100MPa, the range of E _v∞ is 100MPa～300MPa, the range of the subgrade modulus non-uniformity coefficient α is 0.05～5.0, and the range of _ne is 1～4; and the rest parameters usually take the following values: Poisson The ratio is 0.35 and the density is 1.8～2.0g/cm ³ . On the basis of a series of assumed E _v (z) parameters, the obtained surface load data is taken as p(r,t) in S2, combined with the Using the calculation method proposed by the present invention, a series of theoretical displacement response calculation results can be obtained, which are compared and matched with the actual displacement curve. According to the matching results, the parameters E _v0 , E _v∞ , α, _ne are adjusted by parameter adjustment algorithms (such as genetic algorithm, swarm intelligence algorithm, etc.), so that the theoretical displacement response is as close as possible to the field measured displacement response. The parameter value corresponding to the most matching theoretical displacement response calculation result is used as the inversion result of the actual subgrade, so as to obtain the modulus distribution law of the subgrade along the depth direction, and realize the subgrade stiffness evaluation.

In addition, a large number of theoretical displacement responses can be obtained in advance to construct a database through the calculation method proposed in the present invention, and the field measured displacement curve can be quickly matched by combining with the currently common BP neural network and database search method, so that the vertical displacement curve of the field roadbed can be quickly obtained. Vector modulus distribution state. The E _v (z) parameters and _ne corresponding to the theoretical displacement response closest to the actual displacement response are obtained, so as to obtain the modulus distribution function along the depth direction of the actual subgrade structure. Different from the traditional method, which can only obtain a single structural modulus, the calculation by the present invention can more accurately reflect the vertical modulus distribution state of the roadbed.

The results of the vertical modulus distribution of the subgrade can be obtained through climatic data and a series of iterative calculations of finite elements, and then the E _v (z) of the present invention can be used to fit the vertical modulus distribution of the subgrade in theory. A series of loads p(r, t) are used to obtain a series of theoretical displacement response calculation results, which are then used for theoretical analysis and calculation. Such as the impact of overloading, calculating the subgrade work area and so on.

The invention is combined with the non-destructive testing of the roadbed, so that the vertical modulus curve of the roadbed can be obtained, and the long-term monitoring of the vertical modulus change curve can evaluate the wetness of the roadbed. The method for evaluating the subgrade humidification through the present invention: subgrade humidification causes the overall reduction of the vertical modulus curve of the subgrade; by regularly detecting the subgrade, and comparing the vertical modulus distribution curves of the subgrade in different periods, it can be determined The humidity state is qualitatively evaluated; if the modulus of the top of the subgrade decreases, it means that the overall humidity of the subgrade increases; if the modulus of the bottom is basically unchanged, it means that the bottom has reached the equilibrium humidity, while the upper part has not yet reached the equilibrium humidity; The non-destructive testing of the subgrade with the mechanical model of the mass-elastic system cannot be realized. The traditional mechanical model can only obtain one modulus, which approximately represents the modulus of the entire subgrade, so the depth distribution law cannot be obtained.

Advantages of the present invention:

1. At present, the non-uniform modulus distribution of subgrade has been recognized by scholars in the field, but in practical application, the traditional methods mostly need to use finite element method to solve the dynamic response calculation considering the non-uniform distribution of subgrade modulus along the vertical direction. Although this method can obtain more accurate results, the calculation speed is slow, and it needs to be re-modeled every time. It is convenient to perform calculations in the field. The present invention introduces the nonlinear distribution state of the modulus along the depth direction, and combines the current ordinary differential and partial differential solution theory to carry out the relevant derivation of the analytical solution. It is realized by programming, the calculation is relatively simple, the complicated process of finite element is avoided, and the rapid calculation on site can be realized to meet the relevant requirements of actual calculation.

2. The prior art 1 (the multi-layer displacement response determination method considering interlayer contact conditions and transverse isotropy) divides the pavement structure into N layers, but the inner modulus of each layer is uniform. The present invention only considers the subgrade, and the subgrade modulus distribution is not uniform, and the modulus in the mechanical model is substantially different. The invention provides a method for determining the road surface displacement under the impact load of non-destructive testing equipment such as PFWD for the subgrade with transverse isotropy and nonlinear modulus distribution. The quantity change curve can be evaluated, and then the internal humidification of the roadbed can be obtained, which overcomes the defect that only the structural modulus can be measured in the traditional non-destructive testing.

If the method for determining the displacement response of the subgrade surface considering the nonlinear distribution of modulus according to the embodiment of the present invention is implemented in the form of a software function module and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on such understanding, the technical solution of the present invention can be embodied in the form of a software product in essence, or the part that contributes to the prior art or the part of the technical solution. The computer software product is stored in a storage medium, including Several instructions are used to cause a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the method for determining the displacement response of the subgrade surface considering the nonlinear distribution of the modulus according to the embodiment of the present invention. The aforementioned storage medium includes: a U disk, a removable hard disk, a ROM, a RAM, a magnetic disk, or an optical disk and other mediums that can store program codes.

The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the protection scope of the present invention. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention are included in the protection scope of the present invention.

Claims (10)

1. a method for determining the subgrade surface displacement response considering the nonlinear distribution of modulus, it is characterized in that, specifically carry out according to the following steps: Step S1, obtaining subgrade material parameters, the subgrade material parameters include: the vertical modulus E _v0 of the subgrade top surface, the subgrade modulus E _v∞ at the infinite depth of the subgrade, the subgrade modulus non-uniformity coefficient α, the subgrade modulus ratio _ne , horizontal Poisson’s ratio μ _h , vertical Poisson’s ratio μ _v , subgrade density ρ and subgrade thickness H; Step S2, apply a dynamic load p(r, t) to the subgrade surface, and establish a cylindrical coordinate system with the center of the dynamic load as the coordinate center, where r represents the radial coordinate, z represents the vertical coordinate, and φ represents the circumferential coordinate. The nonlinear distribution of the subgrade modulus function E _v (z) = E _v0 + (E _{v∞ -} E _v0 )(1-e ^-αz ); the calculation also considers the transverse isotropy and the vertical nonlinearity of the subgrade modulus Solution of Surface Displacement Response of Distributed Subgrade in Laplace-Hankel Domain

ξ and s are Hankel integral transform and Laplace integral transform coefficients respectively; Step S3, the displacement response solution based on the current subgrade mechanical model in the Laplace-Hankel domain

Calculate the surface displacement response of the current subgrade model in the time domain. 2. a kind of subgrade surface displacement response determination method considering nonlinear distribution of modulus according to claim 1, is characterized in that, in described step S2,

In formula (1-1), A ₁ , A ₂ , A ₃ , and A ₄ are boundary condition parameters,

Represents the coefficient of the n-th series corresponding to m _i in the infinite series of displacement components in the r direction, i=1, 2, 3, 4; n is the serial number of the series, and m ₁ ～m ₄ are intermediate parameters; Ξ A dimensionless coefficient representing the relative magnitudes of the subgrade top and bottom moduli;

In formula (1-3)～(1-6), Θ ₁ , Θ ₂ , Θ ₃ are respectively intermediate calculation parameters; Θ ₁ =α ² C ₃₃ C ₄₄ ;

Θ ₃ =(C ₄₄ ξ ² +θ)(C ₁₁ ξ ² +θ); Among them, θ is an intermediate calculation parameter, θ=ρs ² /E _v∞ ; ξ and s are the Hankel integral transformation and Laplace integral transformation coefficients, respectively, after the integral transformation, the time t becomes s; C ₁₁ , C ₁₂ , C ₁₃ , C ₃₃ , C ₄₄ are parameters that characterize the transverse isotropic characteristics of subgrade materials in axisymmetric problems; C ₁₁ = _kne (1- _ne μ _v ),

C ₁₃ = _kne μ _v (1+μ _h ),

C ₄₄ =0.5(1+μ _v ) ⁻¹ ;

k represents the intermediate variable. 3. A method for determining the displacement response of a subgrade surface considering the nonlinear distribution of modulus according to claim 2, wherein in the step S2, the boundary condition parameters A ₁ , A ₂ , A ₃ , and A ₄ are specifically according to Determined by:

is the result of Hankel-Laplace transformation of the load p(r,t), Ψ _1i , Ψ _2i , Ψ _3i , Ψ _4i are the intermediate calculation parameters,

Represents the coefficient of the n-th series corresponding to m _i in the infinite series of displacement components in the z-direction. 4. a kind of subgrade surface displacement response determination method considering nonlinear distribution of modulus according to claim 3, is characterized in that, in described step S2,

Specifically determined as follows:

In formulas (1-12) to (1-15),

is the intermediate calculation parameter, calculated according to formula (1-16)～(1-21); when E _v0 <E _v∞ , according to formula (1-20), when E _v0 >E _v∞ , according to formula ( 1-21) calculation, in formula (1-21), δ is an intermediate calculation parameter, δ=E _v0 /E _v∞ ;

5 . The method for determining the displacement response of a subgrade surface considering the nonlinear distribution of modulus according to claim 3 , wherein the method for determining the intermediate parameters m ₁ to m ₄ is: 6 . On the basis of the subgrade mechanics model, combined with the dynamic balance equation of the axisymmetric problem, the physical equation representing the axisymmetric transverse isotropy, and the geometric equation, the partial differentials shown in equations (1-22) to (1-23) can be obtained. equation set:

Among them, ur = ur ( _r , z, t), u _z = _uz ( _r , z, t) are the displacement components along the radius r direction of the cylindrical coordinate system and the z direction of the subgrade depth; t is the time; σ _r =σ _r (r,z,t),σ _z =σ _z (r,z,t),σ _φ =σ _φ (r,z,t) are respectively in the r direction, the z direction and the φ direction. stress component; τ _zr is the shear stress in the zr plane; First, perform Laplace transform on t in equations (1-22) to (1-23); perform first-order Hankel transform on r on the Laplace transformed equation (1-22); 23) The zero-order Hankel transform about r is carried out; the system of variable coefficient ordinary differential equations such as equations (1-24) ~ (1-25) is obtained, and the independent variable z is omitted in the formula;

in,

represents the result of Hankel- _Laplace first-order transformation of ur,

Represents the result of Hankel-Laplace zero-order transformation of u _z ; In order to solve the system of ordinary differential equations with variable coefficients shown in equations (1-24) to (1-25), the variable ψ is introduced: ψ=Ξe ^-αz (1-26) Replacing the independent variable z in the subgrade modulus function E _v (z)=E _v0 +(E _v∞ -E _v0 )(1-e ^-αz ) with the variable ψ, the expression of E _v (ψ) is as follows (1 -27) shown:

At the same time, through the derivative theory, formula (1-28) can be obtained; where f(z) represents any function with z as an independent variable:

Substituting equations (1-27) to (1-28) into equations (1-24) to (1-25), we can obtain the equations shown in (1-29) to (1-32) with the ψ argument Ordinary differential equations with variable coefficients; where equations (1-29)～(1-30) correspond to the case of E _v0 <E _v∞ , and equations (1-31)～(1-32) correspond to the situation where E _v0 >E _v∞ situation; the independent variable ψ is omitted in the formula; wherein, δ is an intermediate calculation parameter, δ=E _v0 /E _v∞ ;

In order to solve the system of ordinary differential equations shown in equations (1-29) to (1-32), based on the infinite series theory, an infinite series is used for approximation, assuming that equations (1-29) to (1-32) are shown in The ordinary differential equations of , there are displacement solutions as shown in equations (1-33) ~ (1-34);

In the formula, m is the parameter, and n is the serial number of the series;

Bringing it into equations (1-29) to (1-32), the equations shown in equations (1-35) to (1-38) can be obtained; equations (1-35) to (1-36) ) corresponds to the case of E _v0 <E _v∞ , and equations (1-37) to (1-38) correspond to the case of E _v0 >E _v∞ ;

Let n=0, and compare the two sides of equations (1-35) to (1-38), the following relationship can be obtained: [α ² C ₄₄ m ² -C ₁₁ ξ ² -θ]a ₀ +αξ(C ₁₃ +C ₄₄ )mb ₀ =0 (1-39) αξ(C ₁₃ +C ₄₄ )ma ₀ -[α ² C ₃₃ m ² -C ₄₄ ξ ² -θ]b ₀ =0 (1-40) If equations (1-39) to (1-40) are established, then m must satisfy equation (1-41);

In the complex number domain, there must be four solutions to equation (1-41) about m, see equations (1-3) to (1-6). 6 . The method for determining the displacement response of a subgrade surface considering the nonlinear distribution of modulus according to claim 5 , wherein the boundary conditions are shown in formulas (1-42) to (1-45):

u _r | _{z = H} = 0 (1-44) u _z | _{z = H} = 0 (1-45) Perform Laplace transform on t and 0-order Hankel transform on r for equations (1-42) and (1-45); perform Laplace transform on t for equations (1-43) and (1-44) and The first-order Hankel transform of r; on this basis, the transformation relations shown in equations (1-2) and (1-27) are introduced, and equations (1-46) to (1-49) can be obtained. boundary conditions;

for

The first derivative of ,

for

The first-order derivative of , the general solution expression is brought into the boundary conditions shown in equations (1-46) to (1-49), and the fourth-order matrix equation shown in equation (1-7) can be obtained. From the matrix equation shown in (1-7), the undetermined coefficients A ₁ , A ₂ , A ₃ , and A ₄ can be obtained. 7. A method for determining the displacement response of a subgrade surface considering the nonlinear distribution of modulus according to claim 1, wherein the step S3 is specifically: by formula (1-51), the current subgrade mechanics model is paired in Laplace -Displacement response solution in Hankel domain

That is, formula (1-1), perform inverse Hankel transform, and then perform inverse Laplace transform on the result of Hankel inverse transformation through formula (1-50), and obtain the surface displacement response of the current mechanical model in the time domain u _z (r ,0,t _l );

in,

express

f(r,0,t _l ) represents u _z (r,0,t _l ); the time step is t _l , t _l =l·T/N _a , T is the total solution time, and _Na is the total solution time The number of time increment steps, l is the incremental step sequence, l≤N; r is the horizontal distance between the calculation point and the load center; Δt is the time increment step; Transformation, L and γ are the parameters calculated by the Durbin method, j is a unit imaginary number, j ² =-1; χ is the upper limit of the numerical integration, A _k is the weight of the kth integration node of the 20-node Gaussian interpolation; x _ik is the Gaussian integral The k-th integral node value of the corresponding i integral sections, x _ik =(i-1)ΔL+x _k , ΔL is the division length of the Gaussian integral section, and x _k is the k-th integral node of the 20-node Gaussian integral value; N ₀ is the total number of segments of the Gaussian integral, N ₀ =ceil(χ/ΔL), obtained by calculating (χ/ΔL) and rounding up; Re(·) is the operation of taking the complex real part; J ₀ ( x _ik r) represents the function result corresponding to the zero-order Bessel function when the independent variable is (x _ik r). 8. A method according to any one of claims 1-7 for obtaining the vertical modulus distribution state of the subgrade on site. 9. An electronic device comprising: memory for storing instructions executable by the processor; and A processor for executing the instructions to implement the method of any of claims 1-7. 10. A computer-readable storage medium for storing a computer program, which can implement the method of any one of claims 1-7 when the computer program is executed.

Images