無拉力彈性地基上矩形板屈曲/后屈曲問題的辛求解方法1)

熊斯浚 鄭新然 梁 立 周 超 趙 巖 李 銳,2)

* (大連理工大學工程力學系,工業裝備結構分析優化與CAE 軟件全國重點實驗室,遼寧大連 116024)

?(中國工程物理研究院化工材料研究所,四川綿陽 621999)

引言

彈性地基板模型常用于描述道路橋梁剛性路面、建筑物承載基礎和微電子器件薄膜-基底結構等重要工程結構的力學行為.由于彈性地基板在承受面內載荷時容易發生屈曲/后屈曲失效,因此對該類問題的研究能夠為相關結構的分析和設計提供重要依據.

就彈性地基本身的數學表征而言,已有大量的研究工作.Winkler[1]在1867 年將彈性地基的反作用力通過離散的彈簧來模擬,彈簧的剛度是均勻的,反力與結構和地基的接觸面積成正比,這即是至今仍廣泛使用的Winkler 地基模型.例如,Zhang 等[2]基于改進的Ritz 法研究了Winkler 地基上功能梯度碳納米管增強復合材料板的屈曲問題.Kiani 等[3]解析求解了局部Winkler 地基上圓板在熱環境下的屈曲問題.Pasternak 地基模型是在Winkler 地基模型基礎上改進的典型模型,由于引入了剪切剛度,可以更好地描述結構在地基上的響應,因此同樣得到了廣泛應用.Kiani 等[4]利用一階剪切理論研究了力-熱載荷作用下Pasternak 地基上夾芯板的屈曲問題.Parida 等[5]研究了Pasternak 地基上功能梯度板的振動及屈曲問題.Sobhy[6]研究了濕熱環境下功能梯度夾芯板在Pasternak 地基上的振動與屈曲問題.此外,還有大量其他彈性地基模型應用于各類問題當中[7-10].

傳統的彈性地基模型在拉、壓工況下均可產生反力,相較于工程中真實存在的無拉力彈性地基存在顯著差異,但無拉力彈性地基需要考慮接觸非線性,計算較為復雜.Ma 等[11-15]針對無拉力彈性地基板的屈曲問題開展了研究,對不同的結構以及載荷形式推導了理論解,并開展了相關試驗研究.Shen 等[16]基于二次攝動方法求解了無拉力彈性地基上的層合板在力-熱載荷作用下的后屈曲問題.Zhong 等[17]基于微分求積方法、Newmark 方法以及Newton 迭代法研究了無拉力彈性地基梁的非線性瞬態熱響應.Wang 等[18]基于有限樣條法求解了無拉力彈性地基上管道的屈曲問題.

對于屈曲和后屈曲問題的求解,其本質是高階偏微分方程復雜邊值問題的求解,而各類數值方法無疑是解決上述問題的重要手段.Thom 等[19]基于有限元方法分析了彈性地基上功能梯度裂紋板的屈曲問題.Chaabani 等[20]基于改進的有限元方法求解了彈性地基上多孔功能梯度夾芯板的屈曲問題.Yang 等[21]采用微分求積方法求解了彈性地基上板的熱屈曲及后屈曲問題.Kumar 等[22]基于伽遼金法研究了彈性地基以及多孔性對功能梯度板屈曲的影響.Lopatin 等[7]采用Ritz 法求解了彈性地基上正交各向異性復合材料矩形板的屈曲問題.Liang 等[23-25]結合Koiter 方法與弧長法,在整個后屈曲平衡路徑上使用漸進展開式,實現了不同缺陷對屈曲載荷的影響分析.

盡管數值方法應用廣泛,但并不能動搖解析方法的地位,因為解析解可作為檢驗數值解精度的基準,且由于能夠精確描述各參量之間的關系而常用于快速優化設計,因而具有重要的理論和應用價值.Shen 等[26-27]基于二次攝動方法對彈性地基上板殼結構的后屈曲問題開展了研究.Khorasani 等[28]研究了彈性地基上蜂窩夾芯板的屈曲問題,并利用Navier法對控制方程進行求解.Yaghoobi 等[29]基于Navier法研究了不同類型彈性地基上石墨烯增強夾芯板的振動、屈曲及彎曲問題.Ruocco 等[30]采用Lévy 法求解了Pasternak 地基上納米板的屈曲問題.Xing等[31-33]提出了分離變量法用于求解矩形板的屈曲、振動等問題.近年來,基于鐘萬勰院士和姚偉岸等開創的辛彈性力學求解體系[34],李銳等提出了辛疊加方法,獲得了一系列復雜約束下板殼結構的振動[35]、屈曲[36-39]和彎曲[40-41]解析解.

針對無拉力彈性地基上矩形板的屈曲問題,本文基于辛方法,將板劃分為若干包含強制邊界條件的板,形成子問題,并在辛空間下利用分離變量與辛本征展開對子問題進行求解,通過子問題邊界處的連續條件獲得板的屈曲模態,進而確定板與地基的接觸狀態;通過迭代求解上述過程,獲得子問題劃分及臨界屈曲載荷的收斂結果.在此基礎上,本文結合辛方法與Koiter 攝動法,對無拉力彈性地基上矩形板的后屈曲問題進行求解,給出了不同彈性地基剛度下板的后屈曲平衡路徑.所得屈曲與后屈曲解均與有限元結果吻合良好.本文所發展的方法求解過程嚴格,無需假定解的形式,不僅為其他求解方法提供了對比基準,還可望為相關工程結構的分析與設計提供理論指導.

1 無拉力彈性地基上矩形板的屈曲分析

1.1 屈曲問題數學模型

如圖1 所示,本文的研究對象為單向均布載荷作用下無拉力彈性地基上的簡支矩形板(“S”代表簡支邊界).沿x方向長度為a,沿y方向寬度為b,沿z方向厚度為h.坐標系原點位于矩形板的角點.地基僅在與板互相接觸時產生反力.

圖1 無拉力彈性地基板示意圖Fig.1 Schematic diagram of a plate on a tensionless elastic foundation

圖1 對應的無拉力彈性地基上矩形板屈曲問題的控制方程為

其中w為板中面沿著z方向的撓度;Nx為板中面沿x方向的均布內力;k為彈性地基剛度;彎曲剛度D可表達為

其中E為彈性模量,ν 為泊松比.

上述問題涉及接觸非線性,相較于傳統的Winkler地基板,其求解的關鍵難點在于獲取板與地基的接觸關系,區分接觸區域與非接觸區域.由于當前研究對象邊界以及外載荷的對稱性,使得接觸/非接觸區域分界線垂直于x軸,因此沿加載方向將整塊板按分界線劃分為若干矩形板的拼接[12].將撓度向上的部分視為無地基板,撓度向下的部分則視為Winkler地基板,不同子問題間通過強制邊界條件約束.基于此思路,本文擬通過迭代確定接觸關系,期間對于屈曲問題的解析求解通過辛方法實現.具體求解流程如圖2 所示: 首先求解Winkler 地基上矩形板的屈曲問題,獲得初始的接觸狀態(屈曲模態);根據初始接觸狀態將板劃分為Winkler 地基板及無地基板的子問題,向上變形的部分為無地基板,向下變形的則為Winkler 地基板;對各子問題施加強制邊界條件,利用辛方法求得各子問題的撓度表達式,再通過各子問題邊界處的連續條件求得屈曲載荷及屈曲模態;由上述屈曲模態確定接觸狀態,若相鄰兩次迭代對應的臨界屈曲載荷及接觸/非接觸區域長度(即每個子問題板的長度)的變化率小于0.01%,則判定求解收斂,否則根據最新的屈曲模態結果重新劃分子問題,繼續計算,直至迭代收斂.

圖2 無拉力彈性地基板屈曲問題求解流程圖Fig.2 Solution flowchart for the buckling of a rectangular plate on a tensionless elastic foundation

1.2 導入Hamilton 體系

將需要考慮接觸狀態的原問題簡化為若干Winkler地基板及無地基板的子問題,各子問題是否為Winkler地基板取決于迭代過程的具體求解結果,圖3 為子問題示意圖.

圖3 子問題劃分示意圖Fig.3 Schematic diagram of the division of subproblems

對于Winkler 地基板在單向面內載荷作用下的屈曲問題,其控制方程為高階偏微分方程,傳統的半逆法需要預先假定解的形式,這使得其適用范圍受限.近年來發展起來的辛方法基于嚴格的理論推導,無需假定解的形式,可解析求解此類問題.以下利用辛方法推導獲取子問題的屈曲解答.

彈性地基上矩形薄板屈曲問題的平衡方程如下式所示,其中令k=0 即可獲得無地基板的控制方程

這里板內彎矩Mx和My,扭矩Mxy,剪力Qx和Qy以及等效剪力Vx和Vy的表達式為

由文獻[42] 的推導方法可將控制方程導入Hamilton 體系,形成如下的矩陣形式

其中狀態向量

通過以上步驟,完成了Winkler 地基板屈曲問題向Hamilton 體系的導入.

1.3 方程求解

利用Hamilton 體系下的分離變量法求解式(7),令

其中Y(y)=[w(y),θ(y),-Vx(y),Mx(y)]T為H 的本征向量,僅與y有關.將式(10)代入式(7)可得本征方程

及方程

其中 μin為H 的本征值.利用下述簡支邊界條件求解式(11)

可求得本征值

其中R=Nx/(2D),βn=nπ/b,n=1,2,3,··· .對應的本征向量為

其中i=1,2,3,4 .式(12)的解答為

由式(15)與式(16),狀態向量Z可展開為

其中待定系數Cni由另一個方向的邊界條件確定.首先有

此外,如圖3 所示,原問題被劃分為若干包含強制邊界條件的子問題,在第i條邊界上的強制邊界條件可表示為

其中Eni與Fni為級數展開系數.由式(18)和式(19)可將式(17)中的待定系數Cni求出,需注意的是此時Eni與Fni仍為未知數.將所得帶有Eni與Fni的式(17)代入以下撓度解的連續條件

其中wi及分別為第i個子問題的撓度及彎矩表達式,由此得到一組關于Eni與Fni的聯立線性方程.為了獲得非零解,令方程的系數矩陣行列式為零,從而求出相應的屈曲載荷,同時也可獲得整體屈曲模態.進一步,通過模態結果判斷接觸關系,即可獲得新的子問題分區,繼而按照圖2 所示步驟開展迭代求解,直至滿足收斂條件,從而獲得無拉力彈性地基上矩形板的屈曲解答.

1.4 屈曲數值結果與討論

通過試算發現,彈性地基剛度越大,達到收斂所需迭代次數就越多.同時,存在多個半波的情況才需要考慮接觸問題,而對于僅存在單個半波的情況,其臨界屈曲模態趨向于不與地基發生接觸,最終結果與無地基板相同.以彈性地基剛度=2000,長寬比a/b=2、彎曲剛度D=1 的板作為對象,驗證本文方法迭代求解的收斂性.迭代過程的屈曲模態結果如圖4 所示,其中曲線Iter.0 代表初始的Winkler 地基板計算結果,用于確定初始分區情況.由圖4 中局部放大圖可知,曲線Iter.6 與Iter.7 幾乎重合,表明求解在7 次迭代后達到收斂.在迭代過程中,整體屈曲模態由Winkler 地基板的4 個半波逐漸變為無拉力彈性地基板的3 個半波,可見無拉力彈性地基對矩形板的屈曲行為具有顯著影響.

圖4 屈曲模態迭代的收斂性驗證(kw=2000,a/b=2)Fig.4 Iterative convergence verification of buckling mode shapes(kw=2000,a/b=2)

為驗證本文方法求解的正確性,基于有限元(FEM)分析軟件Abaqus 開展屈曲分析,計算模型中采用S4 單元,網格尺寸為寬度的1/100.由于線性屈曲分析中無法實現無拉力彈性地基設置,因此需按照圖2 所示流程進行迭代計算.與本文解析方法的思路一致,通過程序實現Abaqus 的參數化建模: 首先通過Winkler 地基板確定初始的子問題分區,并根據不同分區的接觸狀態設定彈性地基剛度,進而進行線性屈曲Buckle 分析,根據分析結果重新劃分子問題,并計算新的屈曲載荷及模態,反復迭代直至屈曲載荷及子問題劃分結果收斂,收斂準則與解析方法一致.結果表明,有限元與本文解析方法的迭代結果高度一致,如圖4 所示(其中圓圈代表有限元計算得到的迭代結果).

本文方法與有限元所得屈曲模態結果如圖5 所示,兩者高度吻合.

圖5 無拉力彈性地基上矩形板的臨界屈曲模態(kw=2000,a/b=2)Fig.5 Critical buckling mode shape of a rectangular plate on a tensionless elastic foundation (kw=2000,a/b=2)

為了進一步驗證本文方法,同時獲得更為普適性的計算結果,表1給出了不同長寬比及彈性地基剛度下板的無量綱臨界屈曲載荷.所有臨界屈曲載荷解均與有限元計算結果吻合良好,進一步驗證了本文結果的準確性.

表1 不同長寬比及彈性地基剛度下矩形板的臨界屈曲載荷-Ncrb2/DTable 1 Critical buckling load factors,-Ncrb2/D,of rectangular plates with different aspect ratios and elastic foundation stiffnesses

圖6 對比分析了彈性地基剛度對屈曲模態的影響.在圖6(a)所示含有兩個半波的算例中,隨著彈性地基剛度的增大,半波的幅值差異越發明顯,與地基發生接觸的區域板的變形相對較小.在圖6(b)所示a/b=2.0 的算例中,彈性地基剛度較小時板的屈曲模態呈現出兩個半波,而系數為2000 時則為3 個半波,可見,更高彈性地基剛度下的臨界屈曲模態趨向于呈現更多半波數.如圖6(c)所示,在彈性地基剛度較小時,屈曲模態呈現出對稱性特點,但當系數增大后,轉變為非對稱的屈曲模態.該現象與文獻中結論一致: 即使結構與邊界條件均對稱,屈曲模態也可能呈現非對稱性[15].圖6(d)再次表明,在不同的彈性地基剛度下,屈曲模態的主要變化在于接觸區域板的變形逐漸減小,而與圖6(a)對比可見,隨著板的長寬比增大,半波數呈現出增多的趨勢.

圖6 不同長寬比矩形板的臨界屈曲模態Fig.6 Critical buckling mode shape of rectangular plate with different aspect ratio

針對無拉力彈性地基上的矩形板,其后屈曲分析更加困難.以下基于Koiter 后屈曲攝動展開理論與辛方法,研究無拉力彈性地基上含缺陷薄板的后屈曲問題.攝動法常用于求解非線性方程,其主要思想是將非線性問題的解用某個小量的漸近形式來表示.

內力-應變關系式(4)和下式

2 無拉力彈性地基上矩形板的后屈曲分析

2.1 后屈曲問題控制方程與求解

基于非線性幾何方程式

以及平衡方程式(3),并引入應力函數F

可得薄板后屈曲問題的非線性控制方程為

其中為初始幾何缺陷.需指出,這里的中面內力Nx、Ny及Nxy通過嚴格推導確定,不同于屈曲問題中的均布內力分布設定.利用攝動方法,將解展開為

其中 ξ 為攝動小參數.單軸壓外載荷Px可展開為

其中Pc為臨界屈曲載 荷,ae及be為后屈曲展開系數.令初始幾何缺陷為板的一階屈曲模態

其中 μ=λ/ξ,λ 是缺陷的幅值.將式(25)代入式(24),有方程的零階控制方程

利用邊界條件可得其解答為[16]

由式(24)、式(25)及式(29)可得一階控制方程

一階方程表達形式與式(1) 一致,替換式(17)中 (1+μ)=Nx,即可獲得一階控制方程的解

同理,2 階控制方程為

將其導入Hamilton 體系,有

構造狀態向量

Hamilton 算子矩陣H 為

利用分離變量,令

將式(37)代入式(33)可得

類似于上文求解思路,由邊界條件

可得式(38)中的本征值

對應的本征向量及其一階Jordan 型為

當m=0 時,對應的本征向量及其滿足邊界條件的Jordan 型為

f可按照辛本征向量展開為

其中G=[g00,g01,g02,g03,···g1m,g3m,···,g4m,g2m,···]T.由式(33)、式(37)、式(38)及式(44)可得

由式(46)的展開系數,可得式(45)的解為

式中待定系數需通過邊界條件確定,由式(42)、式(43)及式(47)～ 式(54)可得2 階方程的解答為

具體地,對各子問題用于求解2 階解的邊界條件為

將各階方程的解代入式(25),可得撓度及應力函數的最終解答.令線性屈曲的最大撓度為w*(x0,y0)=h.由于分叉屈曲明顯是對稱的,因此ae=0,對矩形板面內區域 Ω 進行積分,可得后屈曲展開系數be[43]

因此,載荷可表示為

改變攝動小參數 ξ 的值即可獲得相應的變形及載荷,從而獲得后屈曲平衡路徑.

2.2 后屈曲數值結果與討論

本節展示無拉力彈性地基上無缺陷及有缺陷矩形板的后屈曲平衡路徑.以彎曲剛度D=1、長寬比a/b=1.5 和厚度h=0.01 的簡支矩形板作為算例,分別展示不同彈性地基剛度下的后屈曲平衡路徑,并與有限元計算結果進行對比.有限元模型基于Riks 算法,采用S4 單元進行模擬,單元尺寸為寬度的1/100,加載邊采用耦合約束;同時改變網格節點位置引入一階線性屈曲模態缺陷,缺陷幅值 λ 為0.05 及0.1 倍板的厚度.無拉力彈性地基則通過施加如圖7 所表征的非線性彈簧單元實現[12],其一端固定、另一端與矩形板連接,在板內施加密集的連接單元從而模擬無拉力彈性地基.

圖7 無拉力彈性地基的力-位移曲線示意圖Fig.7 Schematic diagram of the force-displacement relation of a tensionless elastic foundation

如圖8 所示,簡支矩形板的后屈曲平衡路徑與有限元模擬結果在小變形階段(即w*/h＜1 時)吻合良好.隨著變形逐漸增大,有限元與本文方法之間的差異有所增加,但總體誤差相對較小,主要原因在于小參數攝動法僅適用于弱非線性問題,而隨著變形的逐漸增大,問題的非線性程度增加,從而導致了一定誤差,但兩者之間整體趨勢保持一致,驗證了本文方法在后屈曲分析時的有效性.

圖8 無缺陷和含缺陷簡支板的后屈曲平衡路徑Fig.8 The post-buckling equilibrium paths of perfect and imperfect simply supported plates

3 結論

在無拉力彈性地基的作用下,矩形板的屈曲/后屈曲行為受到顯著影響.本文首先求解了無拉力彈性地基上矩形板的屈曲問題: 針對接觸非線性難題,將原問題按照接觸狀態進行分區,劃分為若干Winkler 地基板與無地基板的子問題,利用連續條件獲得屈曲載荷及屈曲模態,進而對子問題重新劃分,迭代計算屈曲載荷及屈曲模態,直至結果收斂.進一步,基于Koiter 后屈曲攝動展開理論與辛方法,本文求解了無拉力彈性地基板的后屈曲問題.對于屈曲和后屈曲問題,均通過算例對比,證明了本文求解的準確性,并得到以下結論: (1) 隨著彈性地基剛度的增大,臨界屈曲模態趨向于呈現更多半波數,接觸區域板的變形逐漸減小;(2) 在彈性地基剛度較小時,屈曲模態呈現出對稱性特點,但當系數增大后,轉變為非對稱的屈曲模態;(3) 隨著板的長寬比增大,半波數呈現出增多的趨勢.本文方法具有理論推導嚴格、計算效率高及無需進行復雜建模等優點,可望推廣至更多復雜屈曲/后屈曲問題的求解.