戚春香，李 瑤

（中國民航大學(xué)機(jī)場學(xué)院，天津 300300）

軸對稱圓形均布荷載作用下無限大薄板的彎矩問題作為彈性力學(xué)的經(jīng)典問題之一，受到眾多學(xué)者關(guān)注[1-4]，彈性小撓度基本理論應(yīng)用較多，借助復(fù)雜函數(shù)處理方法，求得圓形荷載作用下無限大板的撓度或彎矩解析解，但實(shí)際工程中不存在“無限大”的情況，建立可代替“無限大”情況的模型具有一定的現(xiàn)實(shí)意義。

應(yīng)用ABAQUS軟件進(jìn)行板的彎矩研究時，無法建立理論中的無限大板，只能建立具有一定尺寸的道面板模型，因此存在邊界效應(yīng)問題，即板的尺寸影響板的彎矩大小，若尺寸選用不當(dāng)，結(jié)果將產(chǎn)生較大誤差，無法進(jìn)行有效模擬分析。國內(nèi)外許多學(xué)者都對有限元法應(yīng)用于剛性道面結(jié)構(gòu)響應(yīng)進(jìn)行了探討，Tabatabaie[5]、Guo等[6]基于Winkler地基上的彈性薄板理論，分別采用梁單元、梁和彈簧單元來模擬接縫的傳荷作用，但由于采用2D模型，不能考慮應(yīng)力沿板厚的分布，因此存在一定的局限性；周正峰等[7]應(yīng)用ABAQUS對Winkler地基上四邊自由的單塊板在飛機(jī)輪載作用下的應(yīng)力收斂性進(jìn)行分析，確定了有限元模型應(yīng)采用的平面尺寸，但沒有考慮荷載作用范圍的影響，結(jié)果具有一定的局限性；譚林等[8]應(yīng)用ABAQUS對單軸均布拉伸載荷作用下有限寬的中心圓孔板進(jìn)行分析，將其與板寬無限大情形下的解析解進(jìn)行比較，給出解析解的適用范圍，但其采用2D單元，具有一定局限性，且沒有得出定量關(guān)系，無法為后續(xù)的模擬提供理論支持；戚春香等[9]基于Winkler地基模型理論，利用有限元軟件模擬了三維道面板模型的撓度傳荷系數(shù)，但并未深入考慮邊界效應(yīng)影響。

通過ABAQUS模擬軸對稱圓形均布荷載作用下板的彎矩，將所得數(shù)值解與理論公式計算軸對稱圓形均布荷載作用下無限大薄板的彎矩所得解析解進(jìn)行比較，研究圓形均布荷載作用下邊界效應(yīng)對道面板彎矩的影響，確定了ABAQUS模擬軸對稱圓形均布荷載作用下板的尺寸與相對剛度半徑l、加載直徑D的關(guān)系，為應(yīng)用ABAQUS模擬無限大板時確定板的尺寸提供數(shù)值參考，同時為進(jìn)一步研究基于邊界效應(yīng)的道面板相關(guān)問題，建立有限元模型提供理論參考。

1 Winkler地基模型理論

根據(jù)Winkler地基模型理論[10-11]，在無限大薄板上施加軸對稱圓形均布荷載并求解徑向單位長度彎矩的計算過程如下。

1.1 計算相對剛度半徑

混凝土板相對剛度半徑l的表達(dá)式為

其中：Ec為混凝土板彈性模量（MN/m2）；k為地基反應(yīng)模量（MN/m3）；v為混凝土板泊松比；h為混凝土板板厚（m）。

1.2 計算徑向單位長度彎矩

徑向單位長度彎矩Mt的表達(dá)式為

其中：R為圓形均布荷載作用半徑（m）；Q為圓形均布荷載合力（N），Q=qπR2，q 為圓形均布荷載；l為相對剛度半徑（m）；r為計算點(diǎn)到荷載中心的距離（m）；t為積分參數(shù)；J0（x）為第1類0階貝塞爾函數(shù)；J1（x）為第1類1階貝塞爾函數(shù)。

2 基本參數(shù)和建模過程

2.1 參數(shù)選取

土基參數(shù)取值：地基反應(yīng)模量k=50 MN/m3。

混凝土板參數(shù)取值：選用正方形板，彈性模量Ec=40 000 MN/m2；泊松比 v=0.15；板厚 h共取 13組，根據(jù)式（1）計算可得相對應(yīng)的13組相對剛度半徑，如表1所示。

表1 板厚及相對剛度半徑Tab.1 Slab thickness and relative rigidity radius m

軸對稱圓形均布荷載參數(shù)取值：均布荷載q=1 MN/m2；加載直徑 D 取 0.2、0.4、1.0、2.0、3.0 m 共5組數(shù)據(jù)，對應(yīng)邊長L取值如表2所示（D為0.2、0.4、1.0 m時，對應(yīng)的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的圓形加載盤，厚度設(shè)置為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個。模型數(shù)量總計1 378個。

2）材料屬性

正方形混凝土板與圓形加載盤材料參數(shù)均為Ec=40 000 MN/m2，v=0.15。

3）相互作用

正方形混凝土板上表面與圓形加載盤下表面采用tie連接。正方形混凝土板下表面設(shè)置彈性地基，k=50 MN/m3。

4）施加荷載

在圓形加載盤上表面施加q=1MN/m2的均布荷載。

5）網(wǎng)格劃分

采用S8R單元，為保持良好鏈接，正方形混凝土板以圓形加載盤的直徑D為基準(zhǔn)，按大小劃分網(wǎng)格，選用雙向向內(nèi)遞增的方式，最小密度為D/20，最大密度為D/2，其中D為加載直徑，此時分析結(jié)果已無明顯變化。

采用S8R單元，圓形加載盤按個數(shù)劃分網(wǎng)格，如圖1所示，圓形內(nèi)部直角邊OA、OB、OC、OD均劃分為20個網(wǎng)格，外部圓弧AB、BC、CD、DA均劃分為40個網(wǎng)格，此時分析結(jié)果已無明顯變化。

6）提取數(shù)據(jù)

分別提取模型中心點(diǎn)即O處的徑向單位長度彎矩Mt。

3 有限元結(jié)果與解析解比較

3.1 數(shù)值解中心點(diǎn)彎矩變化及其與解析解的比較

對于不同加載直徑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時，中心點(diǎn)彎矩Mt隨邊長L的變化及其與解析解的比較。同理，圖 2（b）、（c）、（d）、（e）分別反映D為0.4 m、1.0 m、2.0 m、3.0 m時各相對剛度半徑下Mt隨L的變化及其與解析解的比較。

圖2中不同相對剛度半徑曲線中心點(diǎn)彎矩Mt的數(shù)值解隨正方形混凝土板邊長L的變化均呈現(xiàn)相同的規(guī)律；不同加載直徑下，各數(shù)值解曲線與解析解之間的差值也呈現(xiàn)出相同的變化趨勢。具體如下：

圖2 數(shù)值解與解析解的比較分析Fig.2 Comparison of numerical and analytic solutions

1）其他條件一定時，隨著混凝土板邊長L的增加，數(shù)值解先逐漸上升到一定峰值，此時數(shù)值解與解析解兩曲線之間的差值最大，此時L為不合理邊長；而后數(shù)值解逐漸降低，最終逐漸收斂于解析解，但仍有微小波動。

2）其他條件一定時，相對剛度半徑l越大，荷載中心徑向單位長度彎矩Mt越大。

3）圖2中各曲線前面數(shù)值解明顯小于解析解的部分，誤差均很大，此時混凝土板模型所取邊長L小于板厚h的10倍，不能稱之為薄板，而理論上式（2）不適用于非薄板的條件，因此利用ABAQUS模擬Winkler地基模型無限大板施加軸對稱圓形均布荷載時，必須滿足模型為薄板的前提條件，即L>10 h。

3.2 彎矩峰值邊長與收斂邊長定量分析

提取圖2中各曲線峰值點(diǎn)和初始收斂點(diǎn)的邊長L，其與加載直徑D的關(guān)系如圖3所示，與相對剛度半徑l的關(guān)系如圖4所示。

圖3 邊長與加載直徑關(guān)系圖Fig.3 Relation between side length and loaded circle diameter

3.2.1 峰值點(diǎn)邊長分析

由圖3峰值點(diǎn)L-D關(guān)系趨勢可得，加載直徑D與峰值點(diǎn)邊長L成冪函數(shù)關(guān)系，即

圖4 邊長與相對剛度半徑關(guān)系圖Fig.4 Relation between side length and relative rigidity radius

由圖4峰值點(diǎn)L-l關(guān)系趨勢可知相對剛度半徑l與峰值點(diǎn)邊長L成線性關(guān)系，再根據(jù)式（3）所得指數(shù)關(guān)系，可得L/D0.25與l成線性關(guān)系，各個加載直徑D下的L/D0.25與相對剛度半徑l的關(guān)系如圖5所示，為統(tǒng)一量綱添加a=1 m。

圖 5 a0.25·L/D0.25與 l關(guān)系圖Fig.5 Relation between a0.25·L/D0.25and l

由圖5可得

其中：c、d為系數(shù)，不同加載直徑D對應(yīng)c、d取值如表3所示。

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

由表3可得，系數(shù)c與直徑D存在冪函數(shù)關(guān)系，即

系數(shù)d與直徑D存在線性關(guān)系，即

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

其中：a=1 m，b為待定系數(shù)。由于多次擬合計算產(chǎn)生誤差較大，因此將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關(guān)系圖示Fig.6 Relation between L·D0.05/a0.05-0.5D1.3/a0.3and l

由圖6可得，當(dāng)荷載中心徑向單位長度彎矩Mt數(shù)值解達(dá)到最大峰值時，正方形混凝土板邊長L、相對剛度半徑l、加載直徑D存在關(guān)系（相關(guān)系數(shù)R2為0.987）為

即

其中，a=1 m。

3.2.2 收斂點(diǎn)邊長分析

由圖3、圖4可得，收斂點(diǎn)與峰值點(diǎn)邊長L與相對剛度半徑l、加載直徑D存在相似關(guān)系，同理可得：當(dāng)荷載中心徑向單位長度彎矩Mt數(shù)值解開始收斂于解析解時，正方形混凝土板邊長L、相對剛度半徑l、加載直徑D存在關(guān)系（相關(guān)系數(shù)R2為0.972）為

即

其中，a=1 m。

3.2.3 最合理邊長表達(dá)式驗(yàn)證

收斂點(diǎn)邊長即為有限元建模最適宜邊長，表達(dá)式如式（11）所示，為保證其合理性，驗(yàn)證如下。

在其他條件相同時，選取不同的相對剛度半徑以及加載直徑，代入式（11）可得該情況下最合理邊長，如表4所示。

分別利用ABAQUS計算表4所示邊長模型所得彎矩數(shù)值解，并利用式（2）計算其對應(yīng)解析解，將兩者進(jìn)行誤差分析，結(jié)果如表5所示。

表4 最合理邊長Tab.4 Most appropriate side length m

表5 誤差分析Tab.5 Error analysis

由表5可得，模型邊長根據(jù)式（11）計算取值所得彎矩數(shù)值解與解析解誤差不超過0.3%，該最適宜邊長表達(dá)式具有合理性和準(zhǔn)確性。

4 結(jié)語

應(yīng)用ABAQUS有限元軟件模擬軸對稱圓形均布荷載作用下板的邊界效應(yīng)對彎矩的影響，同時將數(shù)值解與解析解相比較，得到以下結(jié)論。

1）圓形均布荷載作用下邊界效應(yīng)對道面板彎矩影響明顯，在其他條件一定的情況下，隨著混凝土板邊長的增加，模擬所得彎矩數(shù)值解呈先上升后下降的趨勢，并最終逐漸收斂于圓形均布荷載作用下無限大薄板彎矩解析解。

2）圓形均布荷載作用下考慮邊界效應(yīng)的影響，混凝土道面板模型最不合理邊長與板的相對剛度半徑呈線性關(guān)系，與加載直徑呈冪函數(shù)關(guān)系，進(jìn)一步得到基于相對剛度半徑和加載直徑的有限元模型最不適邊長表達(dá)式，即峰值點(diǎn)邊長表達(dá)式，此時ABAQUS模擬所得彎矩數(shù)值解與理論公式所得解析解誤差最大，所得結(jié)果不具有可信度。

3）圓形均布荷載作用下考慮邊界效應(yīng)的影響，混凝土道面板模型最適宜邊長與板的相對剛度半徑呈線性關(guān)系，與加載直徑呈冪函數(shù)關(guān)系，進(jìn)一步得到基于相對剛度半徑和加載直徑的有限元模型最適宜邊長表達(dá)式，經(jīng)驗(yàn)證，該表達(dá)式所得邊長對應(yīng)的彎矩數(shù)值解與解析解誤差很小，可應(yīng)用于有限元建模。