CN112214817A - 考虑层间条件以及横观各向同性的多层位移响应确定方法 - Google Patents

考虑层间条件以及横观各向同性的多层位移响应确定方法 Download PDF

Info

Publication number
CN112214817A
CN112214817A CN202011046882.7A CN202011046882A CN112214817A CN 112214817 A CN112214817 A CN 112214817A CN 202011046882 A CN202011046882 A CN 202011046882A CN 112214817 A CN112214817 A CN 112214817A
Authority
CN
China
Prior art keywords
layer
pavement
mechanical model
hankel
laplace
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011046882.7A
Other languages
English (en)
Other versions
CN112214817B (zh
Inventor
张军辉
范海山
刘杰
张石平
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Changsha University of Science and Technology
Original Assignee
Changsha University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Changsha University of Science and Technology filed Critical Changsha University of Science and Technology
Priority to CN202011046882.7A priority Critical patent/CN112214817B/zh
Publication of CN112214817A publication Critical patent/CN112214817A/zh
Application granted granted Critical
Publication of CN112214817B publication Critical patent/CN112214817B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Structural Engineering (AREA)
  • Computational Mathematics (AREA)
  • Civil Engineering (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Architecture (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种考虑层间条件以及横观各向同性的多层位移响应确定方法，先输入一个N(N≥1)层层状路面体系的不同结构层的材料参数，所述材料参数包括：水平模量、垂直模量、水平泊松比、垂直泊松比、水平粘性系数、垂直粘性系数、相邻两层的层间滑移系数、密度以及厚度；然后对N层层状路面体系表面施加动力荷载，根据路面结构层的层数以及是否存在刚性基岩判断当前N层层状路面体系所属的路面力学模型，并依据当前N层层状路面体系的不同结构层的材料参数计算其所属的路面力学模型在Laplace‑Hankel域下的位移响应解，进而计算得到当前N层层状路面体系在时域范围内的表面位移响应。计算速度快、计算稳定性好，计算结果准确。

Description

考虑层间条件以及横观各向同性的多层位移响应确定方法
技术领域
本发明属于道路工程技术领域，涉及一种考虑层间接触条件及横观各向同性的多层动力体系位移响应确定方法。
背景技术
多层层状动力体系的力学响应计算一直是路面力学的研究热点，也是实现精准模量反算需要解决的难题。目前，常用FWD(Falling Weight Deflectometer，落锤式弯沉仪)对路面结构性能进行检测，而模量反算也正是在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域下的位移响应解
Figure BDA0002708265480000021
步骤S3、基于当前N层层状路面体系所属的路面力学模型在Laplace-Hankel域下的位移响应解
Figure BDA0002708265480000022
计算当前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域下的位移响应解
Figure BDA0002708265480000031
路面力学模型依据当前N层层状路面体系的层数以及是否存在刚性基岩分为单层半无限空间力学模型、单层带基岩力学模型、多层半无限空间力学模型、多层带基岩力学模型：如果N层层状路面体系的路面结构层的层数为1且不存在基岩，则其所属的路面力学模型为单层半无限空间力学模型；如果N层层状路面体系的路面结构层的层数为多层且不存在基岩，则其所属的路面力学模型为多层半无限空间力学模型；如果N层层状路面体系的路面结构层的层数为1且存在基岩，则其所属的路面力学模型为单层带基岩力学模型；如果N层层状路面体系的路面结构层的层数为多层且存在基岩，则其所属的路面力学模型为多层带基岩力学模型。
先通过式(8)和式(9)计算得到C11、C13、C33、C44，再按照式(7)计算中间参数
Figure BDA0002708265480000032
再按式(6)计算中间参数α和δ：
Figure BDA0002708265480000033
Figure BDA0002708265480000041
Figure BDA0002708265480000042
Figure BDA0002708265480000043
其中，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域下的位移响应解
Figure BDA0002708265480000044
a.计算单层半无限空间力学模型在Laplace-Hankel域下的位移响应解
Figure BDA0002708265480000045
先按式(5)计算中间参数β和ψ，最后按照式(4)计算在特定ξ和s下的
Figure BDA0002708265480000046
Figure BDA0002708265480000047
Figure BDA0002708265480000051
其中，
Figure BDA0002708265480000052
为依次对表面荷载p(r,t)的r和t进行Hankel变换以及Laplace变换后的荷载表达式，
Figure BDA0002708265480000053
r表示所求解的位置距离荷载中心的水平距离，t表示时间；J0(ξr)表示自变量为ξr时零阶贝塞尔函数对应的函数结果。由于当前N层层状路面体系所属的力学模型为单层半无限空间力学模型，因此式(6)～(9)中n＝1，即将当前不带基岩的单层路面体系的材料参数带入式(6)～(9)中计算。
b.计算单层带基岩力学模型在Laplace-Hankel域下的位移响应解
Figure BDA0002708265480000054
先按式(13)计算中间参数A11、A13、A21、A23、γ1、γ2，按式(12)计算中间参数Δ，再按式(11)计算与边界条件相关的待定系数A、B、C、D，最后按照式(10)计算
Figure BDA0002708265480000055
Figure BDA0002708265480000056
其中，待定系数A、B、C、D按照式(11)计算：
Figure BDA0002708265480000057
中间参数Δ按照式(12)计算：
Figure BDA0002708265480000058
式(12)中，A11、A12、A13、A14、A21、A23、γ1、γ2均为中间参数，其按照式(13)计算：
Figure BDA0002708265480000061
由于当前N层层状路面体系所属的力学模型为单层带基岩力学模型，因此式(6)～(9)以及式(11)～(12)中n＝1，即将当前带基岩的单层路面体系的厚度带入式(6)～(9)以及式(11)～(12)中计算待定系数A、B、C、D以及中间参数Δ。
c.计算多层半无限空间力学模型在Laplace-Hankel域下的位移响应解
Figure BDA0002708265480000062
先按式(20)计算图1所示的多层路面力学体系的第1至N-1层的中间参数K1～K6后，然后将其带入传递矩阵
Figure BDA0002708265480000063
的计算公式即式(17)～(19)中分别计算其第1至N-1层的考虑层间条件的传递矩阵，然后将该N层层状路面体系中第N层的材料参数带入式(16)求取中间参数M1～M7，得到中间矩阵M；再按式(15)计算综合传递矩阵[Lij]，最后按式(14)计算
Figure BDA0002708265480000064
Figure BDA0002708265480000065
Figure BDA0002708265480000066
Figure BDA0002708265480000071
Figure BDA0002708265480000072
Figure BDA0002708265480000081
Figure BDA0002708265480000082
Figure BDA0002708265480000091
d.计算多层带基岩力学模型在Laplace-Hankel域下的位移响应解
Figure BDA0002708265480000092
先按式(20)计算带基岩的多层层状路面体系的第1至N层的中间参数K1～K6后，按式(17)～(19)计算其第1至N-1层的考虑层间条件的传递矩阵
Figure BDA0002708265480000093
然后将其第N层的材料参数(水平模量EhN、垂直模量EvN、水平泊松比μhN、垂直泊松比μvN、层间滑移系数αN、密度ρN以及厚度hN，)带入Tn中计算该带基岩的多层层状路面体系的第N层的传递函数TN，再按式(22)计算考虑基岩影响的多层综合传递矩阵[Rij]，最后按式(21)计算
Figure BDA0002708265480000094
Figure BDA0002708265480000095
Figure BDA0002708265480000096
步骤S3、基于当前N层层状路面体系所属的路面力学模型在Laplace-Hankel域下的位移响应解
Figure BDA0002708265480000097
计算当前N层层状路面体系在时域范围内的表面位移响应，具体实现过程如下：
先采用式(3)对当前层状路面体系对应的路面力学模型在Laplace-Hankel域下的位移响应解
Figure BDA0002708265480000098
进行Hankel逆变换，然后采用式(2)计算中间结果，最后采用式(1)对Hankel逆变换的结果进行Laplace逆变换，求得当前N层层状路面体系在时域范围内的表面位移响应：
Figure BDA0002708265480000101
Figure BDA0002708265480000102
Figure BDA0002708265480000103
其中，其中，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是单位虚数，j2＝-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个积分节点值；
Figure BDA0002708265480000104
表示ξ＝xik、h＝0时在Hankel-Laplace域中的竖向位移解；
Figure BDA0002708265480000105
表示在Laplace域内的竖向位移解，m为序列数；Re为取复数实部运算；J0(xikr)表示自变量为xikr时零阶贝塞尔函数对应的函数结果。式(3)为采用高斯插值积分进行Hankel逆变换，以上计算过程均采用编程方式实现，这里仅介绍数值计算方法。所述数值计算方法均为编程实现，不限定编程语言，本发明实施例计算的数据均通过MATLAB 2016a计算得到。
关于步骤S2中全部公式的推导方法如下所示：
1)单层横观各向同性通解公式推导
由于本发明实施例的理论模型为轴对称模型，从轴对称动力平衡控制方程出发对单层通解公式进行推导：
Figure BDA0002708265480000111
其中，σ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是时间。
并且结合几何方程：
Figure BDA0002708265480000112
以及物理方程：
Figure BDA0002708265480000113
其中，ε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)所示。
联立(23)、(24)、(25)三式，可以得到：
Figure BDA0002708265480000114
其中，
Figure BDA0002708265480000115
表示在Hankel-Laplace域内的竖向位移。
对式(26)进行Laplace-Hankel变换后，可以得到：
Figure BDA0002708265480000121
其中，
Figure BDA0002708265480000122
是对ur(r,z,t)进行Hankel积分变换而得到的，
Figure BDA0002708265480000123
是对uz(r,z,t)进行Hankel积分变换而得到的，
Figure BDA0002708265480000124
Figure BDA0002708265480000125
是对
Figure BDA0002708265480000126
进行Laplace变换而得到的，
Figure BDA0002708265480000127
是对
Figure BDA0002708265480000128
进行Laplace变换而得到的，
Figure BDA0002708265480000129
ξ是Hankel积分变换的常数，s是Laplace积分变换的常数。
求解式(27)，并结合式(7)进行化简，按照现代常微分方程组求解理论可以解得：
Figure BDA00027082654800001210
其中，相关参数表达式如式(6)、(7)所示。
通过(24)、(25)两式，结合Hankel-Laplace积分变换，可以得到下列关系式：
Figure BDA00027082654800001211
其中，
Figure BDA00027082654800001212
是对σz(r,z,t)进行Hankel-Laplace积分变换得到的，
Figure BDA00027082654800001213
是对τzr(r,z,t)进行Hankel-Laplace积分变换得到，
Figure BDA00027082654800001214
和
Figure BDA00027082654800001215
的表达式类似于
Figure BDA00027082654800001216
和
Figure BDA00027082654800001217
结合式(28)、(29)，可得到下列表达式：
Figure BDA0002708265480000131
因此，单层横观各向同性的通解公式可以由下列矩阵表达式表示：
Figure BDA0002708265480000132
其中，参数K1～K6的形式如式(20)所示。
2)半无限空间即单层体系位移响应计算公式推导
带入半无限空间的边界条件：
Figure BDA0002708265480000133
联立式(31)和(32)，可以得到下式：
Figure BDA0002708265480000134
求解方程组(33)，可得
Figure BDA0002708265480000135
表达式如下：
Figure BDA0002708265480000136
Figure BDA0002708265480000137
3)带基岩的单层体系公式推导
带入该情况下的边界条件：
Figure BDA0002708265480000141
其中，h0为刚性下卧层上覆土层厚度。
联合式(31)和(34)，可得：
Figure BDA0002708265480000142
求解矩阵方程组(35)，可得
Figure BDA0002708265480000143
表达式如下：
Figure BDA0002708265480000144
Figure BDA0002708265480000145
Figure BDA0002708265480000146
Figure BDA0002708265480000151
4)传递矩阵公式推导
传递矩阵按照下式进行推导：
Figure BDA0002708265480000152
结合式(31)和(36)，可得传递矩阵元素如下所示：
Figure BDA0002708265480000161
Figure BDA0002708265480000162
5)层间接触条件计算：
按式(37)来表征层间滑移状态：
Figure BDA0002708265480000171
式中，
Figure BDA0002708265480000172
表示在Hankel-Laplace域内第n层底面的位移解，
Figure BDA0002708265480000173
表示在Hankel-Laplace域内第n+1层顶面的位移解，
Figure BDA0002708265480000174
表示在Hankel-Laplace域内第n层底面的剪应力解；其中，当αn＝0时，为完全连续状态，当αn＝0.99时，表示完全滑动状态；βn为与相邻两层材料的有关系数，与模量与泊松比有关，按下式进行计算：
Figure BDA0002708265480000175
由式(37)可得相邻结构层上下表面力学响应存在如下关系：
Figure BDA0002708265480000176
式中，Cn为考虑层间接触状态的矩阵，
Figure BDA0002708265480000177
将式(39)带入式(36)，可得相邻两层间的应力应变传递关系：
Figure BDA0002708265480000178
令
Figure BDA0002708265480000179
Tn为第n层的传递矩阵，
Figure BDA00027082654800001710
称
Figure BDA00027082654800001711
为考虑接触条件的传递矩阵，此外，对于第N层而言，有
Figure BDA00027082654800001712
6)多层横观各向同性体系计算公式推导：
多层系统的边界条件如下：
Figure BDA0002708265480000181
式中，
Figure BDA0002708265480000182
为Hankel-Laplace域内的第一层的竖向应力，
Figure BDA0002708265480000183
为Hankel-Laplace域内的第一层的水平切向应力，
Figure BDA0002708265480000184
为Hankel-Laplace域内的第N层的水平位移，
Figure BDA0002708265480000185
为Hankel-Laplace域内的第N层的垂直位移。
将式(36)变为下面所示的表达式：
Figure BDA0002708265480000186
其中，参数M1～M7表达式如式(16)所示。
结合式(40)、(41)、(42)，第N层的控制参数必须满足AN＝CN＝0：
Figure BDA0002708265480000187
A,B,C,D为通解公式中与力学模型以及边界条件相关的参数，根据实际情况的不同而不同，只要求解出来A,B,C,D待定参数，就可以得到特定条件下的解。半无限空间求解思想是无穷远处待定系数A和C为零(z为无穷大时，位移为零，因此解的正向次幂前的系数必定为零，否则不满足边界条件)，因此需要将应力响应向量通过左乘[M]矩阵变为系数向量；而存在刚性基岩时，边界条件为第N层一定深度处的位移为零，因此不需要进行转换，直接令应力响应向量中位移为零构建矩阵方程，两者边界条件的引入不一致，因此求解存在区别。
其中，矩阵[Lij]的元素可通过式(15)、(16)计算得到，BN，CN分别为第N层的控制参数，通过求解式(43)，可得
Figure BDA0002708265480000188
表达式如下：
Figure BDA0002708265480000189
7)多层横观各向同性带刚性下卧层体系计算公式推导：
层间连续条件与多层横观各向同性体系相同，但系统的边界条件如下：
Figure BDA0002708265480000191
结合传递矩阵式(36)以及边界条件式(44)，可得下列方程组：
Figure BDA0002708265480000192
其中，
Figure BDA0002708265480000193
表示在Hankel-Laplace域内的第N层的水平切应力，
Figure BDA0002708265480000194
表示在Hankel-Laplace域内的第N层的竖向应力；矩阵[Rij]可由式(22)计算得到，求解式(45)可得
Figure BDA0002708265480000195
表达式如下：
Figure BDA0002708265480000196
在计算过程中，很多情况下将出现数据溢出的问题导致计算无法进行，按照下述方法对步骤S2中的计算进行修正：
所述多层体系
Figure BDA0002708265480000197
的计算中，由于存在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的计算结果：
Figure BDA0002708265480000198
Figure BDA0002708265480000199
其中，c代表为防止数据溢出的控制参数。
所述
Figure BDA0002708265480000201
的计算中，由于需要传递矩阵
Figure BDA0002708265480000202
的反复相乘，且传递矩阵往往数据过大，造成计算结果的不稳定性。这里对式(19)中的
Figure BDA0002708265480000203
采用式(48)进行修正计算，计算结束后采用式(48)计算的
Figure BDA0002708265480000204
替换式(19)的计算结果：
Figure BDA0002708265480000205
其中，
Figure BDA0002708265480000206
为
Figure BDA0002708265480000207
的最小矩阵元素，
Figure BDA0002708265480000208
为
Figure BDA0002708265480000209
中第i行第j列的矩阵元素。通过除以自身最小的元素对矩阵进行整体缩放，可以有效防止由于数值过大而数据溢出无法计算，在保留传递矩阵法计算速度的同时大大降低出现数据溢出的可能性。
在进行Hankel-Laplace逆变换时，需要确定高斯积分积分上限χ以及积分区间长度ΔL，按照如下经验公式进行计算，所述经验公式为本发明实施例大量试算的结果：
Figure BDA00027082654800002010
Figure BDA00027082654800002011
在求解Ψk时，由于计算机存储精度的限制，可能会出现一些不可计算的情况(如不可避免的会出现Ψk分母为零导致最终计算结果无效的情况)，因此在采用式(1)对Hankel逆变换的结果进行Laplace逆变换前，需要先对Ψk进行有效性检验并标记出其中的出现问题的数据，然后对出现问题的Ψk进行修正计算，最后再利用进行修正计算后的Ψk结合式(1)对Hankel逆变换的结果进行Laplace逆变换。采用式(51)对少量离散的问题数据进行修正，采用式(52)对连续出现的问题数据进行修正：
Figure BDA00027082654800002012
Figure BDA00027082654800002013
修正完毕后，使用计算所得
Figure BDA00027082654800002014
替换之前的问题数据Ψk。
采用有限元软件ABAQUS进行对比计算。所设路面结构为四层，其中路面-基层、底基层-土基的层间接触状态分别假设为完全连续以及完全滑动，ABAQUS中采用摩擦系数加以约束；路面、基层、底基层的材料假设为横观各向同性，土基假设为粘弹性；路面结构参数如表1所示，进行ABAQUS建模得到ABAQUS有限元模型。
表1路面结构计算参数表
Figure BDA0002708265480000211
注：T＝60ms；Na＝240；Δt＝0.5ms，荷载为加载时间为30ms，峰值为707kPa的半正弦波，并且荷载半径r0＝0.15m，荷载表达式为：p(r,t)＝707sin(2πt/0.06)，(0s≤t≤0.03s)；求解的表面位移响应位置分别为r＝0,300,600,1600mm。
整个ABAQUS有限元模型采用轴对称方式进行建模，为了模拟层状体系水平方向的无限性，对于边界部分采用无限元CIN3D8单元，有限尺寸部分采用C3D8单元，计算结果如图3所示，图3中MATLAB为采用本发明实施例方法所得计算结果，ABAQUS为采用有限元模型的计算结果，ABAQUS有限元可以作为一种广泛的检验计算正确性的方法，可以看到，两者对于同一个力学模型的计算结果较为接近，因此本发明实施例提供的方法具有较好的计算结果，能够满足考虑层间接触条件以及横观各向同性动力体系表面位移相应计算的相关要求。在计算中，采用Intel i7 8750H/8G笔记本测试得到ABAQUS完成以此计算需要2个小时，而采用本发明的方法仅需2min，大大降低了计算成本。
传统传递矩阵法存在数据溢出问题，即计算结果超过计算机存储极限，因此只能计算多层半空间体系的结果，但是经本发明实施例改进后，对于多层带基岩的结果也可以进行计算，改进后的传递矩阵法大大增加了原有方法的适用范围，提高了传递矩阵法的适用性。
以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。

Claims (10)

1.考虑层间条件以及横观各向同性的多层位移响应确定方法，其特征在于，按照以下步骤进行：
步骤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域下的位移响应解
Figure FDA0002708265470000011
步骤S3、基于当前N层层状路面体系所属的路面力学模型在Laplace-Hankel域下的位移响应解
Figure FDA0002708265470000012
计算当前N层层状路面体系在时域范围内的表面位移响应。
2.根据权利要求1所述的考虑层间条件以及横观各向同性的多层位移响应确定方法，其特征在于，所述步骤S2中的路面力学模型包括单层半无限空间力学模型、单层带基岩力学模型、多层半无限空间力学模型、多层带基岩力学模型；如果N层层状路面体系的路面结构层的层数为1且不存在基岩，则其所属的路面力学模型为单层半无限空间力学模型；如果N层层状路面体系的路面结构层的层数为多层且不存在基岩，则其所属的路面力学模型为多层半无限空间力学模型；如果N层层状路面体系的路面结构层的层数为1且存在基岩，则其所属的路面力学模型为单层带基岩力学模型；如果N层层状路面体系的路面结构层的层数为多层且存在基岩，则其所属的路面力学模型为多层带基岩力学模型。
3.根据权利要求2所述的考虑层间条件以及横观各向同性的多层位移响应确定方法，其特征在于，所述步骤S2中N层层状路面体系所属的力学模型为单层半无限空间力学模型时，其在Laplace-Hankel域下的位移响应解
Figure FDA0002708265470000013
的计算值为：
Figure FDA0002708265470000014
Figure FDA0002708265470000015
其中，
Figure FDA0002708265470000021
为对动力荷载p(r,t)的r和t进行Hankel变换以及Laplace变换后的荷载，
Figure FDA0002708265470000022
r表示所求解的位置距离荷载中心的水平距离，t表示时间；J0(ξr)表示自变量为ξr时零阶贝塞尔函数对应的函数结果；C11、C13、C33、C44均为描述轴对称横观各向同性物理方程中描述应力应变关系的参数；β、ψ、α、δ、
Figure FDA0002708265470000023
Figure FDA0002708265470000024
和
Figure FDA0002708265470000025
均为中间参数；ξ为Hankel积分变换参数，s为Hankel积分变换参数；
所述步骤S2中N层层状路面体系所属的力学模型为单层带基岩力学模型时，其在Laplace-Hankel域下的位移响应解
Figure FDA0002708265470000026
的计算为：
Figure FDA0002708265470000027
其中，待定系数A、B、C、D按照式(11)计算，中间参数Δ按照式(12)计算，式(12)中，A11、A12、A13、A14、A21、A23、γ1、γ2均为中间参数，其按照式(13)计算：
Figure FDA0002708265470000028
Figure FDA0002708265470000029
Figure FDA0002708265470000031
其中，C13、C33、C44均为描述轴对称横观各向同性物理方程中描述应力应变关系的参数；α、δ、
Figure FDA0002708265470000032
和
Figure FDA0002708265470000033
均为中间参数；ξ为Hankel积分变换参数；由于当前N层层状路面体系所属的力学模型为单层带基岩力学模型，因此式(11)～(12)中n＝1，即将当前路面体系的厚度带入式(11)～(12)中计算待定系数A、B、C、D以及中间参数Δ。
4.根据权利要求2所述的考虑层间条件以及横观各向同性的多层位移响应确定方法，其特征在于，所述步骤S2中N层层状路面体系所属的力学模型为多层半无限空间力学模型时，其在Laplace-Hankel域下的位移响应解
Figure FDA0002708265470000034
的计算为：
先按式(20)计算N层层状路面体系的第1至N-1层的中间参数K1～K6后，并按式(17)～(19)分别计算其第1至N-1层的考虑层间条件的传递矩阵
Figure FDA0002708265470000035
然后将该N层层状路面体系中第N层的材料参数带入式(16)求取中间参数M1～M7，得到中间矩阵M；再按式(15)计算综合传递矩阵[Lij]，最后按式(14)计算
Figure FDA0002708265470000036
Figure FDA0002708265470000037
Figure FDA0002708265470000038
Figure FDA0002708265470000041
Figure FDA0002708265470000042
Figure FDA0002708265470000051
Figure FDA0002708265470000052
Figure FDA0002708265470000061
其中，
Figure FDA0002708265470000062
为对动力荷载p(r,t)的r和t进行Hankel变换以及Laplace变换后的荷载，
Figure FDA0002708265470000063
r表示所求解的位置距离荷载中心的水平距离，t表示时间；J0(ξr)表示自变量为ξr时零阶贝塞尔函数对应的函数结果；Tn为N层层状路面体系中第n层的传递矩阵，r0表示上附动力荷载p(r,t)的作用半径；C13、C33、C44均为描述轴对称横观各向同性物理方程中描述应力应变关系的参数；α、δ、
Figure FDA0002708265470000064
和
Figure FDA0002708265470000065
均为中间参数；ξ为Hankel积分变换参数。
5.根据权利要求4所述的考虑层间条件以及横观各向同性的多层位移响应确定方法，其特征在于，所述步骤S2中N层层状路面体系所属的力学模型为多层带基岩力学模型时，其在Laplace-Hankel域下的位移响应解
Figure FDA0002708265470000066
的计算为：
先按式(20)计算带基岩的多层层状路面体系的第1至N层的中间参数K1～K6后，按式(17)～(19)计算其第1至N-1层的考虑层间条件的传递矩阵
Figure FDA0002708265470000067
然后将其第N层的材料参数带入Tn中计算该带基岩的多层层状路面体系的第N层的传递函数TN，再按式(22)计算考虑基岩影响的多层综合传递矩阵[Rij]，1≤i,j≤4最后按式(21)计算
Figure FDA0002708265470000068
Figure FDA0002708265470000069
Figure FDA0002708265470000071
6.根据权利要求3～5任一项所述的考虑层间条件以及横观各向同性的多层位移响应确定方法，其特征在于，所述中间参数α、δ、
Figure FDA0002708265470000072
和
Figure FDA0002708265470000073
通过式(6)～(7)计算，C11、C13、C33、C44通过式(8)～(9)计算：
Figure FDA0002708265470000074
Figure FDA0002708265470000075
Figure FDA0002708265470000076
Figure FDA0002708265470000077
其中，Ehn(s)表示KELVIN粘弹性模型在Laplace积分变换下在水平方向上的结果即水平方向的粘弹性算子，Evn(s)表示KELVIN粘弹性模型在Laplace积分变换下在竖直方向上的结果即竖直方向的粘弹性算子；
如果当前N层层状路面体系所属的力学模型为单层半无限空间力学模型或单层带基岩力学模型，则式(6)～(9)中n＝1，即将当前路面体系的材料参数带入式(6)～(9)中计算α、δ、
Figure FDA0002708265470000081
C11、C13、C33、C44。
7.根据权利要求3～5任一项所述的考虑层间条件以及横观各向同性的多层位移响应确定方法，其特征在于，所述步骤S2中需要对
Figure FDA0002708265470000082
采用式(46)进行修正计算，计算结束后采用式(46)计算的
Figure FDA0002708265470000083
一一对应替换
Figure FDA0002708265470000084
Figure FDA0002708265470000085
的计算结果：
Figure FDA0002708265470000086
Figure FDA0002708265470000087
其中，c代表为防止数据溢出的控制参数。
8.根据权利要求4或5所述的考虑层间条件以及横观各向同性的多层位移响应确定方法，其特征在于，所述步骤S2中对式(19)中的
Figure FDA0002708265470000088
需要采用式(48)进行修正计算，计算结束后采用式(48)计算的
Figure FDA0002708265470000089
替换式(19)的计算结果：
Figure FDA00027082654700000810
其中，
Figure FDA00027082654700000811
为
Figure FDA00027082654700000812
的最小矩阵元素，
Figure FDA00027082654700000813
为
Figure FDA00027082654700000814
中第i行第j列的矩阵元素。
9.根据权利要求1～5任一项所述的考虑层间条件以及横观各向同性的多层位移响应确定方法，其特征在于，所述步骤S3的实现过程如下：
先采用式(3)对当前层状路面体系对应的路面力学模型在Laplace-Hankel域下的位移响应解
Figure FDA00027082654700000815
进行Hankel逆变换，然后采用式(2)计算中间结果，最后采用式(1)对Hankel逆变换的结果进行Laplace逆变换，求得当前N层层状路面体系在时域范围内的表面位移响应：
Figure FDA00027082654700000816
Figure FDA0002708265470000091
Figure FDA0002708265470000092
其中，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是单位虚数，j2＝-1；N0是高斯积分的分段总数，N0＝[χ/ΔL]，[·]是向上取整函数，即[χ/ΔL]表示对χ/ΔL向上取整，χ是数值积分的上限，ΔL是高斯积分区段划分长度；Ak是20节点高斯插值的第k个积分节点的权重；xik是高斯积分对应的i个积分区段的第k个积分节点值，xik＝(i-1)ΔL+xk，xk为20节点高斯积分的第k个积分节点值；
Figure FDA0002708265470000093
表示ξ＝xik、h＝0时在Hankel-Laplace域中的竖向位移解；
Figure FDA0002708265470000094
表示在Laplace域内的竖向位移解，m为序列数；Re为取复数实部运算；J0(xikr)表示自变量为xikr时零阶贝塞尔函数对应的函数结果。
10.根据权利要求9所述的考虑层间条件以及横观各向同性的多层位移响应确定方法，其特征在于，所述高斯积分积分上限χ以及积分区间长度ΔL，按照式进(49)～(50)进行计算：
Figure FDA0002708265470000095
Figure FDA0002708265470000096
在采用式(1)对Hankel逆变换的结果进行Laplace逆变换前，需要先对Ψk进行有效性检验并标记出其中的出现问题的数据，然后对出现问题的Ψk进行修正计算，最后再利用进行修正计算后的Ψk结合式(1)对Hankel逆变换的结果进行Laplace逆变换；
对所述出现问题的Ψk进行修正计算时，采用式(51)对少量离散的问题数据进行修正，采用式(52)对连续出现的问题数据进行修正：
Figure FDA0002708265470000101
Figure FDA0002708265470000102
修正完毕后，使用计算所得
Figure FDA0002708265470000103
替换之前的问题数据Ψk。
CN202011046882.7A 2020-09-29 2020-09-29 考虑层间条件以及横观各向同性的多层位移响应确定方法 Active CN112214817B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011046882.7A CN112214817B (zh) 2020-09-29 2020-09-29 考虑层间条件以及横观各向同性的多层位移响应确定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011046882.7A CN112214817B (zh) 2020-09-29 2020-09-29 考虑层间条件以及横观各向同性的多层位移响应确定方法

Publications (2)

Publication Number Publication Date
CN112214817A true CN112214817A (zh) 2021-01-12
CN112214817B CN112214817B (zh) 2022-07-08

Family

ID=74051414

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011046882.7A Active CN112214817B (zh) 2020-09-29 2020-09-29 考虑层间条件以及横观各向同性的多层位移响应确定方法

Country Status (1)

Country Link
CN (1) CN112214817B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112765825A (zh) * 2021-01-27 2021-05-07 长沙理工大学 考虑横观各向同性及模量非线性分布的路基表面位移响应确定方法及应用、设备及存储介质
CN113609447A (zh) * 2021-08-20 2021-11-05 浙江大学建筑设计研究院有限公司 椭圆荷载沥青路面结构表面处力学响应快速计算方法
CN114896714A (zh) * 2022-03-31 2022-08-12 暨南大学 一种基于锥凹形压头接触的涂层结构界面破坏预测方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070143020A1 (en) * 2005-12-05 2007-06-21 Schlumberger Technology Corporation Transversely isotropic model for wellbore stability analysis in laminated formations
FR2992066A1 (fr) * 2012-06-15 2013-12-20 Commissariat Energie Atomique Capteur de champ magnetique a force de laplace
CN107764644A (zh) * 2017-09-30 2018-03-06 交通运输部公路科学研究所 基于路面材料模量应力和应变依赖模型的沥青路面结构分析当量方法
CN110749442A (zh) * 2019-01-29 2020-02-04 石家庄铁道大学 Laplace小波自适应稀疏表示的滚动轴承故障特征提取方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070143020A1 (en) * 2005-12-05 2007-06-21 Schlumberger Technology Corporation Transversely isotropic model for wellbore stability analysis in laminated formations
FR2992066A1 (fr) * 2012-06-15 2013-12-20 Commissariat Energie Atomique Capteur de champ magnetique a force de laplace
CN107764644A (zh) * 2017-09-30 2018-03-06 交通运输部公路科学研究所 基于路面材料模量应力和应变依赖模型的沥青路面结构分析当量方法
CN110749442A (zh) * 2019-01-29 2020-02-04 石家庄铁道大学 Laplace小波自适应稀疏表示的滚动轴承故障特征提取方法

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112765825A (zh) * 2021-01-27 2021-05-07 长沙理工大学 考虑横观各向同性及模量非线性分布的路基表面位移响应确定方法及应用、设备及存储介质
CN112765825B (zh) * 2021-01-27 2022-07-08 长沙理工大学 考虑模量非线性分布的路基表面位移响应确定方法
CN113609447A (zh) * 2021-08-20 2021-11-05 浙江大学建筑设计研究院有限公司 椭圆荷载沥青路面结构表面处力学响应快速计算方法
CN113609447B (zh) * 2021-08-20 2023-09-29 浙江大学建筑设计研究院有限公司 椭圆荷载沥青路面结构表面处力学响应快速计算方法
CN114896714A (zh) * 2022-03-31 2022-08-12 暨南大学 一种基于锥凹形压头接触的涂层结构界面破坏预测方法
CN114896714B (zh) * 2022-03-31 2024-05-17 暨南大学 一种基于锥凹形压头接触的涂层结构界面破坏预测方法

Also Published As

Publication number Publication date
CN112214817B (zh) 2022-07-08

Similar Documents

Publication Publication Date Title
CN112214817B (zh) 考虑层间条件以及横观各向同性的多层位移响应确定方法
Le et al. Surface corrections for peridynamic models in elasticity and fracture
Liu et al. Isogeometric analysis of laminated composite and functionally graded sandwich plates based on a layerwise displacement theory
Kakouris et al. Phase‐field material point method for brittle fracture
Shlyannikov et al. Characterization of crack tip stress fields in test specimens using mode mixity parameters
Pierron et al. The virtual fields method: extracting constitutive mechanical parameters from full-field deformation measurements
Gao et al. Ductile tearing in part-through cracks: experiments and cell-model predictions
Daxini et al. A review on recent contribution of meshfree methods to structure and fracture mechanics applications
CN103955604A (zh) 一种含裂纹金属梯度材料剩余强度预测方法
CN103838913B (zh) 曲线箱梁弯桥的有限单元法
CN106991219A (zh) 一种考虑三维分形的法向界面刚度预测方法
CN109117504B (zh) 一种双向功能梯度曲壳振动分析方法
Grover et al. An efficient C0 finite element modeling of an inverse hyperbolic shear deformation theory for the flexural and stability analysis of laminated composite and sandwich plates
Kumari et al. Three-dimensional static analysis of Levy-type functionally graded plate with in-plane stiffness variation
CN114861519B (zh) 复杂地质条件下初始地应力场加速优化反演方法
Majidi et al. On the use of the extended finite element and incremental methods in brittle fracture assessment of key-hole notched polystyrene specimens under mixed mode I/II loading with negative mode I contributions
CN114330067A (zh) 一种软基水闸有限元模型修正方法
Rakočević et al. A computational method for laminated composite plates based on layerwise theory
Peters Yield functions taking into account anisotropic hardening effects for an improved virtual representation of deep drawing processes
Ali et al. Analytical and numerical buckling analysis of rectangular functionally-graded plates under uniaxial compression
Hosseini et al. A meshless collocation method on nonlinear analysis of functionally graded hyperelastic plates using radial basis function
Rao et al. Fuzzy logic-based expert system to predict the results of finite element analysis
Helfen et al. Numerical multiscale modelling of sandwich plates
CN112651153A (zh) 一种确定晶体塑性有限元模型材料参数的方法
Kinner-Becker et al. A simulation-based analysis of internal material loads and material modifications in multi-step deep rolling

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant