在許多工程應用中，梁都是由彈性構件進行支承的，比如一些軸通常由滾珠、滾柱或軸頸軸承進行支承，一些較大的梁由彈性墻進行支承，還有一大類梁支承在地基上，一般將這類問題稱為Winkler地基。

單排滾珠軸承可以看做是：在每個軸承處都有一個節點，并將軸承剛度加到單元剛度矩陣中對應垂直自由度的對角位置上（見圖1a）：而對于滾柱或軸頸軸承，還要考慮相應的轉動（力矩）剛度。 對于較寬的軸頸軸承和Winkler地基，我們使用支承介質的單位長度上的剛度s來描述（見圖1b）。在支承介質所作用的長度范圍內，總勢能將會多出下面一項。

（1）

圖1 彈性支承 當我們將代人離散化模型后，上式將變為

（2）

其中H為Hermite插值函數

（3）

從公式（2）中的求和公式里，可以看出其中的單元剛度矩陣項，即

（4）

因為常規的梁單元剛度矩陣也是通過最小勢能原理就行推導，對應的單元應變能和剛度矩陣的表達式為

（5）

（6）

對公式（4）進行積分之后，可得

（6）

對于具有彈性地基支撐的單元，這一剛度矩陣需要加入到公式（6）梁單元對應的傳統的單元剛度矩陣中，即。而矩陣就是彈性地基的一致剛度矩陣。 如果考慮到軸向變形，則為

（6）

對應的只需要用零元素擴充軸向自由度，因為梁的彎曲和軸向變形不相互耦合。 這次課提供的彈性地基梁的有限元案例，是這種雙排樁支護結構的matlab有限元編程計算分析。模型的計算簡圖如圖2所示。

圖2 模型計算簡圖 從結構上分析，雙排樁支護結構如同嵌入土中的門式剛架，與單排樁的懸臂結構、多支點的混合支護結構、重力式擋土結構等支護形式的受力機理有明顯的差異。開挖面下的樁體受到側向的地基抗力在模型中簡化為土彈簧反力，樁底的約束由具體的土質條件和計算模型來確定，綜合考慮上述幾個方面的問題后，包括樁頂結點的處理，彈性地基梁單元剛度矩陣的確定以及荷載列陣的確定等，就可以建立基本的平衡方程求解樁身各結點位移和樁身內力。 對于坑底下部的彈性地基梁單元的剛度矩陣常規方法是按集中剛度的原則先確定在單元結點處的等效彈簧剛度然后疊加到總剛。如果將上述計算的彈簧剛度按結點位置相應的疊加到總剛那么只能加到總體剛度矩陣對角線元素上。 事實上這樣處理彈性地基梁單元的方法不是十分合理。上述方法只是在結點設置彈簧沒有考慮到被動區土壓力也是沿樁體連續分布的并不等效于各個集中彈簧力的作用并且等效剛度不是簡單的將某個區域的剛度求和。為此本文提出一種思路就是考慮單元兩結點間布滿彈簧再將這些分布的彈簧反力轉化為等效結點力建立彈性地基梁單元結點力與結點位移的關系最終可得到基于一般梁單元修正后的彈性地基梁單元剛度矩陣，如下式所示：

（7）

剛度矩陣具體推導過程請看論文《雙排樁開挖過程的改進有限元分析方法》。

本算例用到的基本參數

（1）土性指標：土體粘聚力 c=12k Pa；內摩擦角 φ=25°；土體重度 γ=19.2k N/m3。土體平均壓縮模量 Es=5MPa，不考慮地下水位的影響。采用單一的土層計算。

（2）基坑開挖深度為 9.0m，前后排樁呈矩形布置，樁直徑為0.8m，樁彈性模量E1=3.0×107kN/m^2，樁間距為2m，前排樁入土深度為 11m，樁長為 20m。

（3）連梁截面尺寸 b×h=800mm×600mm，連梁彈性模量 E2=3.0×107k N/ m^2，連梁之間的距離等于兩樁間距，兩排樁的排距為 2.0m，樁頂與連梁按剛接考慮。

（4）彈簧的反力系數，m=4000kN/m^3，樁底采用單鏈桿支承約束，以此替代樁土之間摩擦力的作用，水平向不約束。

（5）土壓力采用朗肯主動土壓力計算，并考慮 10k Pa 的地面施工超載，坑底以上為三角形的分布，基坑底面以下為矩形分布。計算彈性地基梁的剛度矩陣相應的代碼如下：

%計算彈性地基梁單元剛度矩陣函數
function [Ke] = FrameElementKe2(A,E,I,R,BarLength)
   global m D
   L=BarLength;
   k=m*D*L;
   k1=k;
   k2=k;
   ke=[E*A/L       0             0          -E*A/L       0                0
       0        12*E*I/L^3+4*k1*L/15+k2*L/12     6*E*I/L^2     0       -12*E*I/L^3+k1*L/12+k2*L/15   6*E*I/L^2;
       0        6*E*I/L^2+k1*L^2/30+k2

由于不同基坑開挖深度為9m，所以九米以上的樁可以認為傳統梁，九米以下的樁為彈性地基梁。因此在進行單元剛度矩陣計算和組裝時要根據單元所在位置分別進行計算和組裝，具體的matlab代碼如下：

%遍歷所有單元，將各單元剛度陣分塊組裝到總體剛度陣
for iEle =1:EleCount
%該單元的兩個節點的編號
n1=ele(iEle,2);n2=ele(iEle,3);
%計算坐標變換矩陣
R=CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],BarLength(iEle));
%計算單元剛度矩陣 Ke=R'*ke*R;局部坐標系下的單元剛度陣轉換為全局坐標下的單元剛度陣
if y(n1)<=deep     
ke= FrameElementKe2(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
else     
ke= FrameElementKe1(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
end
%將各單元剛度分塊組裝到總剛相應位置
eleDof=[n1*3-2:n1*3,n2*3-2:n2*3];
K(eleDof,eleDof)=K(eleDof,eleDof)+ke; 
end

最終計算得到的前后排樁隨深度的撓度和彎矩分別如圖3所示，與文獻中的結果基本一致，另外本案例還對雙排樁支撐結構的彎矩圖進行繪制，如圖4所示。

圖3 前后排樁不同深度處水平位移

圖4 前后排樁不同深度處彎矩

圖5 彎矩圖

審核編輯：郭婷