本发明属于道路工程技术领域，涉及一种考虑层间接触条件及横观各向同性的多层动力体系位移响应确定方法。

多层层状动力体系的力学响应计算一直是路面力学的研究热点，也是实现精准模量反算需要解决的难题。目前，常用fwd(fallingweightdeflectometer，落锤式弯沉仪)对路面结构性能进行检测，而模量反算也正是在fwd现场测试数据的基础上依靠路面力学模型以及参数调整算法来实现的，因此，要实现精准反算，需要一种快速、准确的力学响应计算方法。而现有层状体系的力学计算中，多采用传递矩阵法、传递系数法以及刚度矩阵法。其中，传递矩阵法以其较小的矩阵规模和易于实现的优势备受青睐，但其存在由于数据溢出导致计算不稳定性的问题。而现有力学模型的选择中，多以层状连续的弹性体系为主；但实际上，由于分层施工以及不同材料层之间所存在的差异性，使得路面结构往往具有层间不完全连续以及横观各向同性的特征，忽略路面结构的以上特征往往会对动力响应计算结果产生较大影响。因此，本发明对传统传递矩阵法进行了改进，并结合路面结构的这些特点，提出了一种考虑层间接触条件及横观各向同性的多层动力体系动力响应计算方法。

本发明实施例的目的在于提供一种考虑层间接触条件及横观各向同性的多层动力体系位移响应确定方法，以解决现有方法计算速度较慢、计算不稳定性，以及目前所采用的层状连续弹性体系力学模型与实际路面结构不符而影响动力响应计算结果的问题。

本发明实施例所采用的技术方案是：考虑层间条件以及横观各向同性的多层位移响应确定方法，按照以下步骤进行：

步骤s1、输入一个n层层状路面体系的不同结构层的材料参数，所述材料参数包括：第n层的水平模量ehn、垂直模量evn、水平泊松比μhn、垂直泊松比μvn、水平粘性系数ηhn、垂直粘性系数ηvn、第n层与(n+1)层的层间滑移系数αn、密度ρn以及厚度hn，1≤n≤n；

步骤s2、对n层层状路面体系表面施加动力荷载p(r,t)，根据路面结构层的层数以及是否存在刚性基岩判断当前n层层状路面体系所属的路面力学模型，并依据当前n层层状路面体系的不同结构层的材料参数计算其所属的路面力学模型在laplace-hankel域下的位移响应解

步骤s3、基于当前n层层状路面体系所属的路面力学模型在laplace-hankel域下的位移响应解计算当前n层层状路面体系在时域范围内的表面位移响应。

本发明实施例的有益效果是，提出了一种考虑层间条件以及横观各向同性的多层位移响应确定方法，对当前n层层状路面体系所属的路面力学模型在laplace-hankel域下的位移响应解的形式进行了改写，改写之后的公式采用矩阵元素相除的形式，然后通过对传递矩阵进行修正以及缩放，使得对改写之后的公式中采用的矩阵中元素进行整体缩小，在不影响最后的应力响应解的同时，通过除以自身最小的元素对考虑了层间接触条件的传递矩阵进行整体缩放，可以有效防止由于数值过大而造成数据溢出无法计算的问题，在保留传递矩阵法计算速度的同时大大降低了出现数据溢出的可能性，有效避免了数据溢出的问题，有效保证计算稳定性。本发明实施例的位移响应确定方法既能够进行多层半空间体系的力学模型的位移计算，又能进行多层带基岩体系的力学模型的位移计算，大大增加了传统传递矩阵法的适用范围，提高了传递矩阵法的适用性。在进行多层体系的位移响应计算时，考虑了层间接触条件及横观各向同性等路面结构特征，使得采用的力学模型与实际路面结构更相符，保证了位移响应计算结果的准确度。相对于整体刚度矩阵法、有限元法等方法，针对现有力学模型的不足提出了一种快速而又准确的位移响应计算方法，有效解决了现有方法计算速度较慢、计算不稳定性，以及目前所采用的层状连续弹性体系力学模型与实际路面结构不符而影响动力响应计算结果的问题。

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例的一个均布荷载作用下横观各向同性轴对称n(n≥1)层层状力学体系示意图。

图2是本发明实施例的动力响应确定方法的流程图。

图3是本发明实施例动力响应确定方法与采用有限元模型abaqus的方法的对比图。

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供一种考虑层间接触条件及横观各向同性的多层动力体系动力响应确定方法，如图2所示，具体按照以下步骤进行：

步骤s1、确定一个n(n≥1)层层状路面体系的材料参数，所述材料参数包括：第n(1≤n≤n)层的水平模量ehn、垂直模量evn、水平泊松比μhn、垂直泊松比μvn、水平粘性系数ηhn、垂直粘性系数ηvn、第n层与(n+1)层的层间滑移系数αn、密度ρn以及厚度hn，如图1所示。此外，需要确定该n层层状体系的最下面一层是否存在刚性基岩，若存在刚性基岩，还需给出基岩埋置深度的厚度。

步骤s2、对n层层状路面体系表面施加动力荷载p(r,t)即输入动力数据p(r,t)，根据路面结构层的层数以及是否存在刚性基岩判断当前n层层状路面体系所属的路面力学模型，并依据当前n层层状路面体系的不同结构层的材料参数计算其所属的路面力学模型在laplace-hankel域下的位移响应解路面力学模型依据当前n层层状路面体系的层数以及是否存在刚性基岩分为单层半无限空间力学模型、单层带基岩力学模型、多层半无限空间力学模型、多层带基岩力学模型：如果n层层状路面体系的路面结构层的层数为1且不存在基岩，则其所属的路面力学模型为单层半无限空间力学模型；如果n层层状路面体系的路面结构层的层数为多层且不存在基岩，则其所属的路面力学模型为多层半无限空间力学模型；如果n层层状路面体系的路面结构层的层数为1且存在基岩，则其所属的路面力学模型为单层带基岩力学模型；如果n层层状路面体系的路面结构层的层数为多层且存在基岩，则其所属的路面力学模型为多层带基岩力学模型。

先通过式(8)和式(9)计算得到c11、c13、c33、c44，再按照式(7)计算中间参数再按式(6)计算中间参数α和δ：

其中，c11、c13、c33、c44均为描述轴对称横观各向同性物理方程中描述应力应变关系的参数，如式(25)所示；式(9)为kelvin模型下的粘弹性算子表达式，ehn(s)和evn(s)表示kelvin粘弹性模型在laplace积分变换下的结果即粘弹性算子，ehn(s)为水平方向的粘弹性算子，evn(s)为竖直方向的粘弹性算子；ξ和s分别为hankel积分变换以及laplace积分变换的参数，在进行单层体系位移响应计算时，式(8)～(9)中n＝1，在进行多层体系位移响应计算时，式(8)～(9)中1≤n≤n。

然后计算当前n层层状路面体系所属的路面力学模型在laplace-hankel域下的位移响应解

a.计算单层半无限空间力学模型在laplace-hankel域下的位移响应解

先按式(5)计算中间参数β和ψ，最后按照式(4)计算在特定ξ和s下的

其中，为依次对表面荷载p(r,t)的r和t进行hankel变换以及laplace变换后的荷载表达式，r表示所求解的位置距离荷载中心的水平距离，t表示时间；j0(ξr)表示自变量为ξr时零阶贝塞尔函数对应的函数结果。由于当前n层层状路面体系所属的力学模型为单层半无限空间力学模型，因此式(6)～(9)中n＝1，即将当前不带基岩的单层路面体系的材料参数带入式(6)～(9)中计算。

b.计算单层带基岩力学模型在laplace-hankel域下的位移响应解

先按式(13)计算中间参数a11、a13、a21、a23、γ1、γ2，按式(12)计算中间参数δ，再按式(11)计算与边界条件相关的待定系数a、b、c、d，最后按照式(10)计算

其中，待定系数a、b、c、d按照式(11)计算：

中间参数δ按照式(12)计算：

由于当前n层层状路面体系所属的力学模型为单层带基岩力学模型，因此式(6)～(9)以及式(11)～(12)中n＝1，即将当前带基岩的单层路面体系的厚度带入式(6)～(9)以及式(11)～(12)中计算待定系数a、b、c、d以及中间参数δ。

c.计算多层半无限空间力学模型在laplace-hankel域下的位移响应解

先按式(20)计算图1所示的多层路面力学体系的第1至n-1层的中间参数k1～k6后，然后将其带入传递矩阵的计算公式即式(17)～(19)中分别计算其第1至n-1层的考虑层间条件的传递矩阵，然后将该n层层状路面体系中第n层的材料参数带入式(16)求取中间参数m1～m7，得到中间矩阵m；再按式(15)计算综合传递矩阵[lij]，最后按式(14)计算

d.计算多层带基岩力学模型在laplace-hankel域下的位移响应解

先按式(20)计算带基岩的多层层状路面体系的第1至n层的中间参数k1～k6后，按式(17)～(19)计算其第1至n-1层的考虑层间条件的传递矩阵然后将其第n层的材料参数(水平模量ehn、垂直模量evn、水平泊松比μhn、垂直泊松比μvn、层间滑移系数αn、密度ρn以及厚度hn，)带入tⁿ中计算该带基岩的多层层状路面体系的第n层的传递函数tⁿ，再按式(22)计算考虑基岩影响的多层综合传递矩阵[rij]，最后按式(21)计算

步骤s3、基于当前n层层状路面体系所属的路面力学模型在laplace-hankel域下的位移响应解计算当前n层层状路面体系在时域范围内的表面位移响应，具体实现过程如下：

先采用式(3)对当前层状路面体系对应的路面力学模型在laplace-hankel域下的位移响应解进行hankel逆变换，然后采用式(2)计算中间结果，最后采用式(1)对hankel逆变换的结果进行laplace逆变换，求得当前n层层状路面体系在时域范围内的表面位移响应：

其中，其中，uz(r,0,tl)为当前n层层状体系所属的路面力学模型在时域范围内的表面位移响应，其对应的时间步为tl，tl＝l·t/na，t为求解总时长，na是求解的总的时间增量步数，l是增量步序列，l≤n；r是计算点位距离荷载中心的水平距离；δt是时间增量步长；式(1)～(2)表示采用一阶自相关检验方法进行laplace逆变换，l和a是一阶自相关检验方法的计算参数，l·na＝50～5000，a·t＝5～10；j是单位虚数，j²＝-1；对于durbin法，最后1/4的时间求解会产生较大的系统误差，通常采用延长求解总时长t来解决，本发明通过测试，取3倍时间长度可以满足精度要求(对hankel-laplace域内的解进行逆变换时，需要给定计算的时长，比如我需要计算80ms动力响应，由于durbin法存在一定的缺陷，因此在实际计算时，计算时长参数一般取2～3倍，本发明实施例取3倍时间长度，即计算时长参数t取240ms，但最后的结果仅取前80ms)n0是高斯积分的分段总数，n0＝[χ/δl]，[·]是向上取整函数，即[χ/δl]表示对χ/δl向上取整，χ是数值积分的上限，δl是高斯积分区段划分长度；ak是20节点高斯插值的第k个积分节点的权重；xik是高斯积分对应的i个积分区段的第k个积分节点值，xik＝(i-1)δl+xk，xk为20节点高斯积分的第k个积分节点值；表示ξ＝xik、h＝0时在hankel-laplace域中的竖向位移解；表示在laplace域内的竖向位移解，m为序列数；re为取复数实部运算；j0(xikr)表示自变量为xikr时零阶贝塞尔函数对应的函数结果。式(3)为采用高斯插值积分进行hankel逆变换，以上计算过程均采用编程方式实现，这里仅介绍数值计算方法。所述数值计算方法均为编程实现，不限定编程语言，本发明实施例计算的数据均通过matlab2016a计算得到。

关于步骤s2中全部公式的推导方法如下所示：

1)单层横观各向同性通解公式推导

由于本发明实施例的理论模型为轴对称模型，从轴对称动力平衡控制方程出发对单层通解公式进行推导：

其中，σr(r,z,t)是r方向上的正应力，σθ(r,z,t)是θ方向上的正应力，σz(r,z,t)是z方向上的正应力；τzr(r,z,t)是z-r平面的剪应力，ur(r,z,t)是r方向上的位移，uz(r,z,t)是z方向上的位移；ρ是密度；t是时间。

其中，εr(r,z,t)是r方向上的正应变，εθ(r,z,t)是θ方向上的正应变，εz(r,z,t)是z方向上的正应变；γzr(r,z,t)是z-r平面上的剪应变；cij(0<i<4，0<j<4)为描述轴对称横观各向同性物理方程中描述应力应变关系的参数，具体表达式如式(8)所示。

其中，是对ur(r,z,t)进行hankel积分变换而得到的，是对uz(r,z,t)进行hankel积分变换而得到的，是对进行laplace变换而得到的，是对进行laplace变换而得到的，ξ是hankel积分变换的常数，s是laplace积分变换的常数。

求解式(27)，并结合式(7)进行化简，按照现代常微分方程组求解理论可以解得：

其中，相关参数表达式如式(6)、(7)所示。

通过(24)、(25)两式，结合hankel-laplace积分变换，可以得到下列关系式：

结合式(28)、(29)，可得到下列表达式：

因此，单层横观各向同性的通解公式可以由下列矩阵表达式表示：

其中，参数k1～k6的形式如式(20)所示。

2)半无限空间即单层体系位移响应计算公式推导

带入半无限空间的边界条件：

联立式(31)和(32)，可以得到下式：

求解方程组(33)，可得表达式如下：

3)带基岩的单层体系公式推导

带入该情况下的边界条件：

其中，h0为刚性下卧层上覆土层厚度。

求解矩阵方程组(35)，可得表达式如下：

传递矩阵按照下式进行推导：

结合式(31)和(36)，可得传递矩阵元素如下所示：

5)层间接触条件计算：

按式(37)来表征层间滑移状态：

式中，表示在hankel-laplace域内第n层底面的位移解，表示在hankel-laplace域内第n+1层顶面的位移解，表示在hankel-laplace域内第n层底面的剪应力解；其中，当αn＝0时，为完全连续状态，当αn＝0.99时，表示完全滑动状态；βn为与相邻两层材料的有关系数，与模量与泊松比有关，按下式进行计算：

由式(37)可得相邻结构层上下表面力学响应存在如下关系：

式中，cn为考虑层间接触状态的矩阵，

将式(39)带入式(36)，可得相邻两层间的应力应变传递关系：

令tⁿ为第n层的传递矩阵，称为考虑接触条件的传递矩阵，此外，对于第n层而言，有

6)多层横观各向同性体系计算公式推导：

多层系统的边界条件如下：

将式(36)变为下面所示的表达式：

其中，参数m1～m7表达式如式(16)所示。

结合式(40)、(41)、(42)，第n层的控制参数必须满足aⁿ＝cⁿ＝0：

a,b,c,d为通解公式中与力学模型以及边界条件相关的参数，根据实际情况的不同而不同，只要求解出来a,b,c,d待定参数，就可以得到特定条件下的解。半无限空间求解思想是无穷远处待定系数a和c为零(z为无穷大时，位移为零，因此解的正向次幂前的系数必定为零，否则不满足边界条件)，因此需要将应力响应向量通过左乘[m]矩阵变为系数向量；而存在刚性基岩时，边界条件为第n层一定深度处的位移为零，因此不需要进行转换，直接令应力响应向量中位移为零构建矩阵方程，两者边界条件的引入不一致，因此求解存在区别。

其中，矩阵[lij]的元素可通过式(15)、(16)计算得到，bⁿ，cⁿ分别为第n层的控制参数，通过求解式(43)，可得表达式如下：

7)多层横观各向同性带刚性下卧层体系计算公式推导：

层间连续条件与多层横观各向同性体系相同，但系统的边界条件如下：

结合传递矩阵式(36)以及边界条件式(44)，可得下列方程组：

其中，表示在hankel-laplace域内的第n层的水平切应力，表示在hankel-laplace域内的第n层的竖向应力；矩阵[rij]可由式(22)计算得到，求解式(45)可得表达式如下：

在计算过程中，很多情况下将出现数据溢出的问题导致计算无法进行，按照下述方法对步骤s2中的计算进行修正：

所述多层体系的计算中，由于存在e^αh与e^δh需要计算，而在计算时，由于α和δ往往过大，造成计算结果过大而溢出。这里需要对e^αh、e^-αh、e^δh、e^-δh采用式(46)进行修正计算，计算结束后采用式(46)计算的e^α′h，e^-α′h，e^δ′h，e^-δ′h一一对应替换e^αh、e^-αh、e^δh、e^-δh的计算结果：

其中，c代表为防止数据溢出的控制参数。

所述的计算中，由于需要传递矩阵的反复相乘，且传递矩阵往往数据过大，造成计算结果的不稳定性。这里对式(19)中的采用式(48)进行修正计算，计算结束后采用式(48)计算的替换式(19)的计算结果：

其中，为的最小矩阵元素，为中第i行第j列的矩阵元素。通过除以自身最小的元素对矩阵进行整体缩放，可以有效防止由于数值过大而数据溢出无法计算，在保留传递矩阵法计算速度的同时大大降低出现数据溢出的可能性。

在进行hankel-laplace逆变换时，需要确定高斯积分积分上限χ以及积分区间长度δl，按照如下经验公式进行计算，所述经验公式为本发明实施例大量试算的结果：

在求解ψk时，由于计算机存储精度的限制，可能会出现一些不可计算的情况(如不可避免的会出现ψk分母为零导致最终计算结果无效的情况)，因此在采用式(1)对hankel逆变换的结果进行laplace逆变换前，需要先对ψk进行有效性检验并标记出其中的出现问题的数据，然后对出现问题的ψk进行修正计算，最后再利用进行修正计算后的ψk结合式(1)对hankel逆变换的结果进行laplace逆变换。采用式(51)对少量离散的问题数据进行修正，采用式(52)对连续出现的问题数据进行修正：

修正完毕后，使用计算所得替换之前的问题数据ψk。

采用有限元软件abaqus进行对比计算。所设路面结构为四层，其中路面-基层、底基层-土基的层间接触状态分别假设为完全连续以及完全滑动，abaqus中采用摩擦系数加以约束；路面、基层、底基层的材料假设为横观各向同性，土基假设为粘弹性；路面结构参数如表1所示，进行abaqus建模得到abaqus有限元模型。

表1路面结构计算参数表

整个abaqus有限元模型采用轴对称方式进行建模，为了模拟层状体系水平方向的无限性，对于边界部分采用无限元cin3d8单元，有限尺寸部分采用c3d8单元，计算结果如图3所示，图3中matlab为采用本发明实施例方法所得计算结果，abaqus为采用有限元模型的计算结果，abaqus有限元可以作为一种广泛的检验计算正确性的方法，可以看到，两者对于同一个力学模型的计算结果较为接近，因此本发明实施例提供的方法具有较好的计算结果，能够满足考虑层间接触条件以及横观各向同性动力体系表面位移相应计算的相关要求。在计算中，采用inteli78750h/8g笔记本测试得到abaqus完成以此计算需要2个小时，而采用本发明的方法仅需2min，大大降低了计算成本。

传统传递矩阵法存在数据溢出问题，即计算结果超过计算机存储极限，因此只能计算多层半空间体系的结果，但是经本发明实施例改进后，对于多层带基岩的结果也可以进行计算，改进后的传递矩阵法大大增加了原有方法的适用范围，提高了传递矩阵法的适用性。

以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。

《线性代数》复习提纲第一部分：基本要求（计算方面）

N阶特殊行列式的计算（如有行和、列和相等）；

矩阵的运算（包括加、减、数乘、乘法、转置、逆等的混合运算）；

求矩阵的秩、逆（两种方法）；解矩阵方程；

含参数的线性方程组解的情况的讨论；

齐次、非齐次线性方程组的求解（包括唯一、无穷多解）；

讨论一个向量能否用和向量组线性表示；

讨论或证明向量组的相关性；

求向量组的极大无关组，并将多余向量用极大无关组线性表示；

将无关组正交化、单位化；

求方阵的特征值和特征向量；

讨论方阵能否对角化，如能，要能写出相似变换的矩阵及对角阵；

通过正交相似变换（正交矩阵）将对称矩阵对角化；

写出二次型的矩阵，并将二次型标准化，写出变换矩阵；

判定二次型或对称矩阵的正定性。

用n^2个元素aij组成的记号称为n阶行列式。

（1）它表示所有可能的取自不同行不同列的n个元素乘积的代数和；

（2）展开式共有n!项，其中符号正负各半；

一阶|α|=α行列式，二、三阶行列式有对角线法则；

N阶（n>=3）行列式的计算：降阶法

定理：n阶行列式的值等于它的任意一行（列）的各元素与其对应的代数余子式乘积的和。