0.   引言

机场道面结构对于飞机轮载作用的力学响应是进行道面厚度设计的基本依据。刚性道面是由板-接缝-基础组成的3D结构体系, 然而现有刚性道面设计方法, 包括中国民用机场刚性道面设计方法, 美国联邦航空局(FAA)刚性道面传统设计方法^[1]和针对新一代大型飞机的设计新方法^[2]等, 采用的结构分析模型都是基于弹性地基板理论或弹性多层体系理论, 难以考虑接缝传荷等因素, 存在不少局限性。有限元法能够考虑有限尺寸板、接缝传荷、任意荷载形式、板底的接触或支承条件, 已成为分析道面结构响应的主流方法。国内外许多学者都对有限元法在刚性道面(或路面)结构响应分析中的应用进行了大量研究, 如Tabatabaie、Huang、Guo等基于Winkler地基上的弹性薄板理论, 分别采用梁单元、弹簧单元、梁及弹簧单元模拟接缝的传荷作用, 开发了ILLI-SLAB、KENSLABS、JSLAB等刚性路面2D有限元分析程序^[3-5]。然而这些模型大多是基于Winkler地基上弹性薄板理论的2D模型, 存在诸多局限性, 如不能考虑应力沿板厚的分布、混凝土板下的多层基础体系等。另一方面, 机场刚性道面结构层厚度大, 接缝传荷机理复杂, 作用的飞机起落架构型复杂、机轮多, 需要建立更加符合实际的道面结构有限元分析模型。

为此, 本文应用ABAQUS大型通用有限元软件, 首先建立Winkler地基上四边自由的单块板, 通过采用不同的单元类型, 划分不同的网格密度, 选取不同的模型平面尺寸, 进行计算结果(应力、挠度)的收敛性分析, 以合理选取这些模型参数; 然后按照结点的贡献面积, 将接缝刚度分配到连接接缝两侧板上对应结点的弹簧单元上, 建立模拟集料嵌锁型接缝和传力杆型接缝的统一模型, 并通过对比所建模型计算的接缝传荷系数与已有回归模型的预估结果, 验证了模型的合理性。

1.   单元类型、网格密度和平面尺寸

1.1   单元类型和网格密度

道面各结构层都采用六面体实体单元进行模拟。ABAQUS单元类型库中包含各种类型的三维实体单元, 其基本性质见表 1。其中, C3D27和C3D27R为变结点单元, 即单元结点数可以通过移除或添加体中心及六个面中心的结点, 使单元结点数在21~27的范围内变化。

表  1  三维六面体单元的基本性质

Table  1.  Properties of 3D hexahedral elements

下载: 导出CSV 

| 显示表格

通过选取表 1中的6种三维实体单元, 每种单元类型划分4种网格密度, 建立总计24个有限元模型进行计算结果的收敛性分析, 以确定最合适的单元类型和网格密度, 提高求解精度和缩短计算时间。收敛性分析采用的道面结构和材料参数见表 2, 采用的荷载为B767-200主起落架上的一个轮载, 荷载参数见表 3, 荷载作用位置为道面板中, 轮印面积简化为正方形。

表  2  道面结构和材料参数

Table  2.  Structural and material parameters of pavement

下载: 导出CSV 

| 显示表格

表  3  飞机荷载参数

Table  3.  Loading data of aircrafts

下载: 导出CSV 

| 显示表格

混凝土板的网格划分为粗、较粗、较细和细4种密度等级(图 1), 对应的h/e(h为板厚, e为单元边长)分别为2、4、6、8。基础不建实体模型, 直接采用ABAQUS接触功能模块中的Elastic Foundation进行模拟。

图  1  网格密度

Figure  1.  Densities of element meshes

下载: 全尺寸图片 幻灯片

图  2  最大应力收敛性

Figure  2.  Convergences of maximum stresses

下载: 全尺寸图片 幻灯片

图  3  最大挠度收敛性

Figure  3.  Convergences of maximum deflections

下载: 全尺寸图片 幻灯片

以无量纲表征的最大应力σ_ih²/P_s(σ_i为板中最大应力)和最大挠度ω_ikl²/P_s(ω_i为板中最大挠度)的计算结果分别见图 2、3。从图中可以看出: 采用线性单元计算得到的σ_ih²/P_s和ω_ikl²/P_s随网格密度的增大迅速收敛, 在网格划分很细(h/e=8)的情况下, 挠度计算结果与二次单元基本相同, 但应力计算结果与二次单元仍有较大差距, 最大误差约为15.4%。采用二次单元在网格划分最粗与最细的情况下, 最大应力σ_ih²/P_s的误差不超过2.5%, 而最大挠度ω_ikl²/P_s完全相同, 表明即使划分很粗的网格, 采用二次单元计算得到的最大应力σ_ih²/P_s仍具有较高的精度。考虑到有限元模型的单元和结点数量随着网格密度的增大迅速增多, 需要的计算时间相应增加, 网格划分并非越细越好。ABAQUS中Serendipity单元不能用于接触分析, 因此, 在考虑各结构层的层间接触时, 应选择Lagrangian单元, 进一步考虑到减缩积分单元C3D27R的计算时间要小于完全积分单元C3D27, 因此, 本文确定采用二次减缩积分单元C3D27R模拟道面各结构层材料, 单元密度满足单元尺寸约为板厚的一半即可, 不要求划分十分细密的网格。

1.2   平面尺寸

一些多轮起落架上的机轮分布范围较广, 有限元分析结果会受到混凝土板平面尺寸大小的影响。为此, 通过调整道面板边长L在3l~7l(l为1.235 m)内变化, 分析表 3中3种机型的起落架(分别代表单轴双轮D、双轴双轮DT和三轴双轮TDT起落架)作用下道面结构力学响应的变化规律, 确定混凝土板的最小尺寸, 使道面板结构响应在混凝土板平面尺寸继续增大时不再有明显变化, 对于采用固定尺寸的混凝土板, 相当于确定有限元模型应包含的板块数量。3种构型起落架的轮印分布和作用位置见图 4, 其中尺寸单位为m。

图  4  轮印分布

Figure  4.  Wheel configurations

下载: 全尺寸图片 幻灯片

3种构型起落架作用下, 不同平面尺寸的混凝土板板底弯拉应力σ_y随横坐标x的变化规律分别见图 5~7(应力提取路径通过最大弯拉应力点)。

图  5  单轴双轮作用下应力σ_y随x的变化

Figure  5.  Changes of σ_ywith x under single axle-biwheel

下载: 全尺寸图片 幻灯片

图  6  双轴双轮作用下应力σ_y随x的变化

Figure  6.  Changes of σ_ywith x under biaxle-biwheel

下载: 全尺寸图片 幻灯片

3种构型起落架作用下, 除混凝土板边长L为3l外, 其他几种平面尺寸计算得到的板底水平拉应力的分布没有明显差异。混凝土板边长分别为4l与7l时, 计算得到的板底最大水平拉应力分别为3.74、3.88 MPa(单轴双轮)、3.73、3.87 MPa(双轴双轮)、4.61、4.77 MPa(三轴双轮), 两者非常接近, L为4l相对于L为7l的计算结果误差仅为3.61%、3.62%和3.35%, 因此, 当混凝土板边长L≥4l时, 模型平面尺寸的影响就可以忽略不计。本文建立的有限元模型为具有接缝的两块5 m×5 m的混凝土板系统, 即使是在三轴双轮起落架的作用下, 该模型的平面尺寸大小仍然是适宜的。

图  7  三轴双轮作用下应力σ_y随x的变化

Figure  7.  Changes of σ_ywith x under triaxle-biwheel

下载: 全尺寸图片 幻灯片

2.   接缝传荷作用的模拟

2.1   接缝传荷刚度的确定

道面接缝通过传力杆系统、集料嵌锁系统或传力杆及集料嵌锁系统传递荷载, 传递荷载的类型都以剪力为主, 因此, 无论是集料嵌锁型接缝还是传力杆型接缝, 均可采取在接缝两侧对应结点设置弹簧单元的方法模拟接缝的剪力传递作用, 从而采用统一的接缝模型模拟两类接缝的传荷作用。接缝传荷能力以接缝刚度q表征, 它表示单位长度接缝内产生单位竖向位移差的剪力。

对于传力杆型接缝, 接缝刚度q可采用Friberg理论方法计算确定^[6]。传力杆的传荷刚度包括混凝土对传力杆的支承和传力杆自身两部分, 反映混凝土对传力杆支承的弹簧刚度D_CI为

$\begin{array}{l}D{}_{\text{C}\text{Ι}}=\frac{4\beta {}^3}{2+\beta b}E{}_{\text{d}}{\rm I}{}_{\text{d}}(1)\\ \beta =\left(\frac{{\rm K}d}{4E{}_{\text{d}}{\rm I}{}_{\text{d}}}\right){}^{1/4}(2)\end{array}$

式中: β为传力杆-混凝土的相对刚度; K为混凝土对传力杆的支承模量, 取4.07×10⁵ MN·m^-3[7]; d为传力杆直径; E_d为传力杆的弹性模量; I_d为传力杆截面的惯性矩, 圆筋时I_d为πd⁴/64;b为接缝缝隙宽。

反映传力杆自身的剪切弹簧刚度C为

$\begin{array}{l}C=\frac{E{}_{\text{d}}{\rm I}{}_{\text{d}}}{b{}^3(1+\phi )}(3)\\ \phi =\frac{12E{}_{\text{d}}{\rm I}{}_{\text{d}}}{G{}_{\text{d}}A{}_{\text{z}}b{}^2}(4)\end{array}$

式中: G_d为传力杆的剪切模量; A_z为传力杆的有效截面面积, 圆筋时A_z为0.225πd²。

包含以上两部分的组合剪切刚度D为

$D=1/\left(\frac{1}{D{}_{C{\rm I}}}+\frac{1}{12C}\right)(5)$

通过式(5)可计算得到一根传力杆的组合刚度, 则接缝单位长度的刚度q为

$q=D/s(6)$

式中: s为接缝中传力杆的间距。

对于集料嵌锁型接缝, 接缝刚度q可根据接缝两侧的弯沉测试得到挠度传荷系数, 由文献[8, -9]中建立的接缝刚度系数与挠度传荷系数的关系式估算, 该方法也同样适用于传力杆型接缝。Crovetti建立的无因次的接缝刚度系数q/kl与挠度传荷系数L_TE的相关关系模型^[8]为

$L{}_{\text{Τ}\text{E}}=\frac{1}{1+1.2(q/kl){}^{-0.849}}(7)$

Zollinger等建立的相关关系模型^[9]为

$L{}_{\text{Τ}\text{E}}=1/[1+10{}^{\frac{0.214-0.183\frac{a}{l}-\lg \left(\frac{q}{kl}\right)}{1.180}}](8)$

式中: a为正方形轮印边长的1/2或圆形轮印的半径。

2.2   接缝传荷刚度的分配

在有限元分析中, 道面实体模型经过离散化生成单元后, 混凝土板侧面的结点可分为板角结点、板边结点和板中结点(图 8), 通过在对应结点设置弹簧单元SPRING2, 并使弹簧刚度方向与接缝剪力方向一致, 即可实现接缝传荷作用的模拟。按照3种不同位置结点对接缝刚度的贡献面积, 将接缝的总刚度分配到每个弹簧单元上。

图  8  板侧面结点布置

Figure  8.  Location of nodes in slab side face

下载: 全尺寸图片 幻灯片

3种位置结点的刚度贡献面积分别为

板角结点

$A{}_{\text{C}}=a{}_{\text{n}}b{}_{\text{n}}$

板边结点

$A{}_{\text{E}}=2a{}_{\text{n}}b{}_{\text{n}}$

板中结点

$A{}_{\text{Ι}}=4a{}_{\text{n}}b{}_{\text{n}}$

板角、板边和板中结点的贡献面积比为1∶2∶4, 因而板角、板边和板中结点的分配刚度可分别设为k_j、2k_j和4k_j。3种位置结点的数量分别为

板角结点

${\rm N}{}_{\text{C}}=4$

板边结点

${\rm N}{}_{\text{E}}=2(n{}_{\text{R}}+n{}_{\text{C}}-4)$

板中结点

${\rm N}{}_{\text{Ι}}=(n{}_{\text{R}}-2)(n{}_{\text{C}}-2)$

3种位置结点的分配刚度与数量的乘积之和应等于接缝的总刚度, 即

$\begin{array}{l}k{}_{\text{j}}{\rm N}{}_{\text{C}}+2k{}_{\text{j}}{\rm N}{}_{\text{E}}+4k{}_{\text{j}}{\rm N}{}_{\text{Ι}}=q\lambda \\ k{}_{\text{j}}=\frac{q\lambda }{4(n{}_{\text{R}}-1)(n{}_{\text{C}}-1)}\end{array}$

式中: λ为接缝长度; n_R、n_C为板侧面结点的行数及列数; a_n、b_n为相邻结点之间的水平与竖向距离。

2.3   接缝模型的验证

图  9  接缝有限元模型

Figure  9.  Finite element model for joint

下载: 全尺寸图片 幻灯片

图  10  应力σ_y随接缝刚度的变化

Figure  10.  Changes of bending stresses with joint stiffnesses

下载: 全尺寸图片 幻灯片

图  11  弯沉ω随接缝刚度的变化

Figure  11.  Changes of deflections with joint stiffnesses

下载: 全尺寸图片 幻灯片

图  12  ABAQUS模型验证

Figure  12.  Verification for ABAQUS model

下载: 全尺寸图片 幻灯片

为了验证本文建立的ABAQUS接缝传荷模型是否合理, 建立了Winkler地基上接缝连接的两块5 m×5 m的道面板的有限元模型, 分析不同接缝刚度下, 接缝一侧道面板板边中部受荷时接缝附近道面结构的力学响应, 比较ABAQUS模型计算的挠度传荷系数L_TE与Zollinger等关系式^[9]是否一致。道面结构和材料参数同表 1(板的尺寸改为5.00 m×5.00 m×0.36 m)。每个单元的平面尺寸为17.86 cm×17.86 cm, 在竖向沿混凝土板厚划分为2层, 每个单元高度为18 cm。接缝有限元模型见图 9。

调整q/kl在10^-4~10⁵的范围内变化, 模拟接缝从几乎不具有传荷能力到挠度传荷系数接近100%。不同接缝刚度对应的应力σ_y随离接缝位置x的变化规律见图 10, 板表面弯沉ω随离接缝位置x的变化规律见图 11, 接缝传荷系数L_TE与接缝刚度系数的对比分析结果见图 12。由图 12可见, ABAQUS计算得到的各点L_TE都基本位于关系曲线上, 表明本文建立的接缝模型合理。

3.   多层基础体系的模拟

ABAQUS可以十分方便地建立基层、底基层、地基等基础结构层的实体模型以及进行层间接触分析。接触是一种典型的与状态相关的非线性行为, 接触状态可以是完全粘结、完全光滑或部分粘结, 当为完全光滑或部分粘结时, 层间界面处往往会发生滑移或张开等位移不连续现象。

在ABAQUS中可以通过定义接触面或接触单元来模拟接触问题, 其中定义接触面的方法相对容易。采用定义接触面的方法进行接触分析的主要步骤包括: 定义接触面; 定义接触属性和接触; 定义边界条件, 消除模型的刚体位移和避免过约束。其中接触属性包括接触面之间的法向作用和切向作用。对于法向作用, ABAQUS中接触压力和间隙的默认关系是“硬接触”, 即接触面之间能够传递的接触压力大小不受限制; 当接触面变为零或负值时, 两个接触面分离, 并且去掉相应结点上的接触约束。对于切向作用, ABAQUS中常用的摩擦模型为库伦摩擦, 即使用摩擦系数来表示接触面之间的摩擦特性。库伦摩擦的计算公式为

τ_crit=μp_n

式中: τ_crit为临界切向力; μ为摩擦系数; p_n为法向接触压力。另外, 在ABAQUS中, Winkler地基模型可直接采用Interaction功能模块中的Elastic Foundation, 而不用建立实体模型, 十分方便。

不同构型起落架作用下混凝土板底拉应力σ_y和板表面弯沉ω分别见图 13、14。可见, 多轮荷载作用引起的σ_y分布受各自机轮分布形式的影响, 呈双峰或三峰形状, 且靠近接缝处的峰值要大于离接缝较远处的峰值, 主要原因是距离接缝位置越近, 道面结构相对越薄弱, 轮载作用下产生的σ_y越大。混凝土板表面弯沉ω受多轮荷载叠加效应的影响, 在起落架覆盖范围内, 弯沉ω较大且变化相对平缓, 而在起落架覆盖范围外, 弯沉ω迅速减小。

4.   实例分析

采用本文建立的有限元模型, 分析表 3中3种机型的一个主起落架作用在道面板边中部时的结构响应。接缝刚度根据机场道面传力杆的布置要求, 按Friberg理论计算为1 006.6 MN·mm^-2。实例中道面的结构和材料参数见表 4。

表  4  实例中道面结构和材料参数

Table  4.  Structural and material parameters in example

道面结构层	参数	数值
水泥混凝土面层	弯拉强度/MPa	5.0
	弯拉弹性模量/MPa	36 000
	泊松比	0.15
	厚度/cm	36
水泥稳定碎石基层	回弹模量/MPa	1 500
	泊松比	0.25
	厚度/cm	30
天然砂粒垫层及土基	反应模量/(MN·m-3)	30

下载: 导出CSV 

| 显示表格

图  13  不同构型起落架的σ_y

Figure  13.  σ_yunder various landing assemblies

下载: 全尺寸图片 幻灯片

图  14  不同构型起落架的ω

Figure  14.  ω under various landing assemblies

下载: 全尺寸图片 幻灯片

5.   结语

(1) 应用ABAQUS建立刚性道面结构有限元模型时, 各结构层宜采用二次积分单元C3D27或C3D27R进行模拟, 相应网格密度满足单元尺寸约为板厚的一半即可, 模型的平面尺寸需大于4倍相对刚度半径。一般情况下, 建立具有接缝的两块5 m×5 m的道面板系统是合适的。

(2) 接缝刚度可根据Friberg理论或挠度传荷系数与接缝刚度系数之间的相关关系确定。接缝刚度分配可按照结点的贡献面积, 分配到连接接缝两侧板上对应结点的弹簧单元SPRING2上, 并设置SPRING2的弹簧方向与接缝传递的剪力方向一致, 模拟传力杆型接缝或集料嵌锁型接缝的传荷作用。

(3) 建立的有限元模型可以全面考虑接缝的传荷作用、大型飞机复杂起落架的多轮荷载作用、板下多层基础结构、各结构层层间接触条件等因素, 具有普遍适用性。