圆形荷载作用下边界效应对道面板弯矩的影响

戚春香，李 瑶

（中国民航大学机场学院，天津 300300）

轴对称圆形均布荷载作用下无限大薄板的弯矩问题作为弹性力学的经典问题之一，受到众多学者关注[1-4]，弹性小挠度基本理论应用较多，借助复杂函数处理方法，求得圆形荷载作用下无限大板的挠度或弯矩解析解，但实际工程中不存在“无限大”的情况，建立可代替“无限大”情况的模型具有一定的现实意义。

应用ABAQUS软件进行板的弯矩研究时，无法建立理论中的无限大板，只能建立具有一定尺寸的道面板模型，因此存在边界效应问题，即板的尺寸影响板的弯矩大小，若尺寸选用不当，结果将产生较大误差，无法进行有效模拟分析。国内外许多学者都对有限元法应用于刚性道面结构响应进行了探讨，Tabatabaie[5]、Guo等[6]基于Winkler地基上的弹性薄板理论，分别采用梁单元、梁和弹簧单元来模拟接缝的传荷作用，但由于采用2D模型，不能考虑应力沿板厚的分布，因此存在一定的局限性；周正峰等[7]应用ABAQUS对Winkler地基上四边自由的单块板在飞机轮载作用下的应力收敛性进行分析，确定了有限元模型应采用的平面尺寸，但没有考虑荷载作用范围的影响，结果具有一定的局限性；谭林等[8]应用ABAQUS对单轴均布拉伸载荷作用下有限宽的中心圆孔板进行分析，将其与板宽无限大情形下的解析解进行比较，给出解析解的适用范围，但其采用2D单元，具有一定局限性，且没有得出定量关系，无法为后续的模拟提供理论支持；戚春香等[9]基于Winkler地基模型理论，利用有限元软件模拟了三维道面板模型的挠度传荷系数，但并未深入考虑边界效应影响。

通过ABAQUS模拟轴对称圆形均布荷载作用下板的弯矩，将所得数值解与理论公式计算轴对称圆形均布荷载作用下无限大薄板的弯矩所得解析解进行比较，研究圆形均布荷载作用下边界效应对道面板弯矩的影响，确定了ABAQUS模拟轴对称圆形均布荷载作用下板的尺寸与相对刚度半径l、加载直径D的关系，为应用ABAQUS模拟无限大板时确定板的尺寸提供数值参考，同时为进一步研究基于边界效应的道面板相关问题，建立有限元模型提供理论参考。

1 Winkler地基模型理论

根据Winkler地基模型理论[10-11]，在无限大薄板上施加轴对称圆形均布荷载并求解径向单位长度弯矩的计算过程如下。

1.1 计算相对刚度半径

混凝土板相对刚度半径l的表达式为

其中：Ec为混凝土板弹性模量（MN/m2）；k为地基反应模量（MN/m3）；v为混凝土板泊松比；h为混凝土板板厚（m）。

1.2 计算径向单位长度弯矩

径向单位长度弯矩Mt的表达式为

其中：R为圆形均布荷载作用半径（m）；Q为圆形均布荷载合力（N），Q=qπR2，q 为圆形均布荷载；l为相对刚度半径（m）；r为计算点到荷载中心的距离（m）；t为积分参数；J0（x）为第1类0阶贝塞尔函数；J1（x）为第1类1阶贝塞尔函数。

2 基本参数和建模过程

2.1 参数选取

土基参数取值：地基反应模量k=50 MN/m3。

混凝土板参数取值：选用正方形板，弹性模量Ec=40 000 MN/m2；泊松比 v=0.15；板厚 h共取 13组，根据式（1）计算可得相对应的13组相对刚度半径，如表1所示。

表1 板厚及相对刚度半径Tab.1 Slab thickness and relative rigidity radius m

轴对称圆形均布荷载参数取值：均布荷载q=1 MN/m2；加载直径 D 取 0.2、0.4、1.0、2.0、3.0 m 共5组数据，对应边长L取值如表2所示（D为0.2、0.4、1.0 m时，对应的L范围一致）。

表2 D、L取值Tab.2 D and L values m

2.2 有限元模型

1）建立模型

模型如图1所示，下部为正方形混凝土板，上部为圆形加载盘，图中尺寸如表2所示，选用3D shell单元。

图1 有限元模型Fig.1 Finite element model

建立直径D=0.2 m的圆形加载盘，厚度设置为1×10-10m，建立不同板厚h的正方形混凝土板，如表1所示，从而得到不同的模型相对刚度半径l，则在加载直径D=0.2 m时共有13组不同厚度（相对刚度半径）的模型，在同一厚度（相对刚度半径）的模型中，再分别建立不同边长L的正方形混凝土板模型，边长大小如表2所示，共20组。则D=0.2m系列模型共计13×20=260个。

按照同样的方法分别建立D=0.4 m系列模型共13×20=260个、D=1.0m系列模型共13×20=260个、D=2.0 m系列模型共13×24=312个、D=3.0 m系列模型共13×22=286个。模型数量总计1 378个。

2）材料属性

正方形混凝土板与圆形加载盘材料参数均为Ec=40 000 MN/m2，v=0.15。

3）相互作用

正方形混凝土板上表面与圆形加载盘下表面采用tie连接。正方形混凝土板下表面设置弹性地基，k=50 MN/m3。

4）施加荷载

在圆形加载盘上表面施加q=1MN/m2的均布荷载。

5）网格划分

采用S8R单元，为保持良好链接，正方形混凝土板以圆形加载盘的直径D为基准，按大小划分网格，选用双向向内递增的方式，最小密度为D/20，最大密度为D/2，其中D为加载直径，此时分析结果已无明显变化。

采用S8R单元，圆形加载盘按个数划分网格，如图1所示，圆形内部直角边OA、OB、OC、OD均划分为20个网格，外部圆弧AB、BC、CD、DA均划分为40个网格，此时分析结果已无明显变化。

6）提取数据

分别提取模型中心点即O处的径向单位长度弯矩Mt。

3 有限元结果与解析解比较

3.1 数值解中心点弯矩变化及其与解析解的比较

对于不同加载直径D，在不同相对刚度半径l下，ABAQUS有限元模型正方形混凝土板边长L对计算径向单位长度弯矩Mt的影响及其与解析解的比较如图2所示。

图 2（a）反映 D=0.2 m 时，l分别取 0.511 0 m、0.692 7 m、0.859 4 m、1.016 0 m、1.164 9 m、1.307 7 m、1.445 4 m、1.578 9 m、1.708 7 m、1.835 4 m、1.959 1 m、2.080 3 m、2.199 2 m时，中心点弯矩Mt随边长L的变化及其与解析解的比较。同理，图 2（b）、（c）、（d）、（e）分别反映D为0.4 m、1.0 m、2.0 m、3.0 m时各相对刚度半径下Mt随L的变化及其与解析解的比较。

图2中不同相对刚度半径曲线中心点弯矩Mt的数值解随正方形混凝土板边长L的变化均呈现相同的规律；不同加载直径下，各数值解曲线与解析解之间的差值也呈现出相同的变化趋势。具体如下：

图2 数值解与解析解的比较分析Fig.2 Comparison of numerical and analytic solutions

1）其他条件一定时，随着混凝土板边长L的增加，数值解先逐渐上升到一定峰值，此时数值解与解析解两曲线之间的差值最大，此时L为不合理边长；而后数值解逐渐降低，最终逐渐收敛于解析解，但仍有微小波动。

2）其他条件一定时，相对刚度半径l越大，荷载中心径向单位长度弯矩Mt越大。

3）图2中各曲线前面数值解明显小于解析解的部分，误差均很大，此时混凝土板模型所取边长L小于板厚h的10倍，不能称之为薄板，而理论上式（2）不适用于非薄板的条件，因此利用ABAQUS模拟Winkler地基模型无限大板施加轴对称圆形均布荷载时，必须满足模型为薄板的前提条件，即L>10 h。

3.2 弯矩峰值边长与收敛边长定量分析

提取图2中各曲线峰值点和初始收敛点的边长L，其与加载直径D的关系如图3所示，与相对刚度半径l的关系如图4所示。

图3 边长与加载直径关系图Fig.3 Relation between side length and loaded circle diameter

3.2.1 峰值点边长分析

由图3峰值点L-D关系趋势可得，加载直径D与峰值点边长L成幂函数关系，即

图4 边长与相对刚度半径关系图Fig.4 Relation between side length and relative rigidity radius

由图4峰值点L-l关系趋势可知相对刚度半径l与峰值点边长L成线性关系，再根据式（3）所得指数关系，可得L/D0.25与l成线性关系，各个加载直径D下的L/D0.25与相对刚度半径l的关系如图5所示，为统一量纲添加a=1 m。

图 5 a0.25·L/D0.25与 l关系图Fig.5 Relation between a0.25·L/D0.25and l

由图5可得

其中：c、d为系数，不同加载直径D对应c、d取值如表3所示。

表 3 c、d 取值Tab.3 c and d values

由表3可得，系数c与直径D存在幂函数关系，即

系数d与直径D存在线性关系，即

将式（5）、式（6）代入式（4），整理后可得

其中：a=1 m，b为待定系数。由于多次拟合计算产生误差较大，因此将L·D0.05/a0.05-0.5D1.3/a0.3与相对刚度半径l重新拟合以确定b，如图6所示。

图6 L·D0.05/a0.05-0.5D1.3/a0.3与l关系图示Fig.6 Relation between L·D0.05/a0.05-0.5D1.3/a0.3and l

由图6可得，当荷载中心径向单位长度弯矩Mt数值解达到最大峰值时，正方形混凝土板边长L、相对刚度半径l、加载直径D存在关系（相关系数R2为0.987）为

即

其中，a=1 m。

3.2.2 收敛点边长分析

由图3、图4可得，收敛点与峰值点边长L与相对刚度半径l、加载直径D存在相似关系，同理可得：当荷载中心径向单位长度弯矩Mt数值解开始收敛于解析解时，正方形混凝土板边长L、相对刚度半径l、加载直径D存在关系（相关系数R2为0.972）为

即

其中，a=1 m。

3.2.3 最合理边长表达式验证

收敛点边长即为有限元建模最适宜边长，表达式如式（11）所示，为保证其合理性，验证如下。

在其他条件相同时，选取不同的相对刚度半径以及加载直径，代入式（11）可得该情况下最合理边长，如表4所示。

分别利用ABAQUS计算表4所示边长模型所得弯矩数值解，并利用式（2）计算其对应解析解，将两者进行误差分析，结果如表5所示。

表4 最合理边长Tab.4 Most appropriate side length m

表5 误差分析Tab.5 Error analysis

由表5可得，模型边长根据式（11）计算取值所得弯矩数值解与解析解误差不超过0.3%，该最适宜边长表达式具有合理性和准确性。

4 结语

应用ABAQUS有限元软件模拟轴对称圆形均布荷载作用下板的边界效应对弯矩的影响，同时将数值解与解析解相比较，得到以下结论。

1）圆形均布荷载作用下边界效应对道面板弯矩影响明显，在其他条件一定的情况下，随着混凝土板边长的增加，模拟所得弯矩数值解呈先上升后下降的趋势，并最终逐渐收敛于圆形均布荷载作用下无限大薄板弯矩解析解。

2）圆形均布荷载作用下考虑边界效应的影响，混凝土道面板模型最不合理边长与板的相对刚度半径呈线性关系，与加载直径呈幂函数关系，进一步得到基于相对刚度半径和加载直径的有限元模型最不适边长表达式，即峰值点边长表达式，此时ABAQUS模拟所得弯矩数值解与理论公式所得解析解误差最大，所得结果不具有可信度。

3）圆形均布荷载作用下考虑边界效应的影响，混凝土道面板模型最适宜边长与板的相对刚度半径呈线性关系，与加载直径呈幂函数关系，进一步得到基于相对刚度半径和加载直径的有限元模型最适宜边长表达式，经验证，该表达式所得边长对应的弯矩数值解与解析解误差很小，可应用于有限元建模。