周長東，郭坤鵬，孟令凱，張曉陽

（北京交通大學 土木建筑工程學院，北京100044）

筒倉作為特種結構，廣泛應用于儲藏各類顆粒材料.筒倉與顆粒散料之間的靜力相互作用十分復雜；同時，確定筒倉與散料之間的靜力相互作用既是筒倉荷載與內力設計的前提，又是研究在地震作用下筒倉與散料顆粒相互作用的基礎.關于散料顆粒對倉壁的壓力計算，我國鋼筋混凝土筒倉設計規(guī)范[1]與歐洲[2]（ISO）和美國[3]規(guī)范存在很大的不同；而且各國規(guī)范對于倉壁壓力的計算都有不足[4].對于淺倉而言，我國采用Rankine理論進行倉內散料對倉壁的壓力設計[1]；而歐洲等規(guī)范則采用Janssen理論進行設計[2－3].對于深倉，盡管我國與歐洲等國都采用Janssen理論進行倉內散料對倉壁的壓力設計，但在設計時，各國規(guī)范所采用的側壓比系數(shù)計算方法卻有很大不同[1－3].因此，為了研究各國規(guī)范筒倉倉壁靜壓力設計值是否具有足夠的安全儲備，有必要對倉壁與散料之間的靜態(tài)壓力作用進行研究.

目前，應用于模擬倉壁與散料之間相互作用的數(shù)值方法主要有有限元法（FEM）和離散元法（DEM）兩種.有限元法發(fā)展較為成熟，其計算結果與精度為廣大工程師所接受；離散元法作為一種相對較新的數(shù)值方法，基于對單個顆粒體力學行為進行研究，也可應用于模擬倉壁與散料散粒體之間的相互作用[5].但 Holst等[6－7]認為，離散元法對于模擬倉壁與散料散粒體之間相互作用有很大的離散性，有限元法則能收獲較好的精度.因此，本文采用有限元方法對倉壁與散料散粒體之間靜力相互作用進行數(shù)值模擬.

在采用有限元方法分析筒倉與散料之間的靜力相互作用時，國內外多數(shù)學者[8－9]都假定散粒體為各向同性，并利用Drucker－Prager理論中黏性系數(shù)、內摩擦角和膨脹角三個參數(shù)來考慮散料的塑性變形階段.然而，散料作為散粒體的集合，其應力狀態(tài)往往在未達到Drucker－Prager理論中定義的屈服應力之前就出現(xiàn)散料內部散粒體重排列的現(xiàn)象[10]；同時，Drucker－Prager理論未能考慮散料初始孔隙比、顆粒硬度及顆粒間應變相對滑移幅值對倉壁壓力分布的影響.因此，本文基于考慮散料初始孔隙比、顆粒硬度及顆粒間應變的亞塑性模型理論，利用大型通用有限元軟件ABAQUS對鋼筋混凝土筒倉倉壁與散料散粒體之間的靜態(tài)壓力作用進行三維有限元模擬，并將模擬結果與我國規(guī)范（GB 50077－2003）[1]、歐洲規(guī)范（EN 1991－4）[2]、美國筒倉設計規(guī)范（ACI 313－97）[3]以及既有試驗數(shù)據(jù)[4]和經(jīng)典筒倉壓力理論進行比較.同時，文中還對倉內顆粒材料內摩擦角、摩擦系數(shù)、初始孔隙比、顆粒硬度及顆粒間應變進行參數(shù)分析，研究其對筒倉倉壁與散料散粒體之間的靜態(tài)壓力作用的影響.

1 數(shù)值模型

1.1 亞塑性本構理論

亞塑性本構理論以連續(xù)介質力學張量函數(shù)為工具，直接建立了應力率與應變率之間的關系，能夠更好地描述顆粒材料初始塑性變形及完全發(fā)展的塑性變形，是一種較新的本構方程表征方法.

目前較為常見的亞塑性本構理論有 Wu－Bauer亞塑性模型、Gudehus－Bauer亞塑性模型、von Wolffersdorff考慮臨界狀態(tài)面的亞塑性本構模型[11]以及由Niemunis和Herle基于Von Wolffersdorff修正的亞塑性本構模型[12]等.

本文所采用的是由Niemunis和Herle基于von Wolffersdorff修正的亞塑性本構模型[12].該本構模型在保留von Wolffersdorff本構模型[11]能夠較好模擬顆粒材料力學行為的基礎上，解決了其在小應力循環(huán)下散粒體材料變形過度積累的棘輪效應問題.下面分別介紹von Wolffersdorff亞塑性本構模型[11]和 由 Niemunis和 Herle基 于 von Wolffersdorff修正的亞塑性本構模型[12].

1.1.1 von Wolffersdorff亞塑性本構模型

von Wolffersdorff亞塑性本構模型[11]的表達式為

標量因子a與臨界狀態(tài)極限面形狀有關；φc為材料內摩擦角.

向密性函數(shù)fb的具體表達式為

式中：hs為顆粒硬度；n為按比例壓縮能量法則時所對應的指數(shù)項；ei為在等向壓縮時，顆粒材料在最小密度時對應的上限孔隙比；ei0為零壓狀態(tài)下顆粒材料可能達到的最大孔隙比；ec0為顆粒材料最大孔隙比；ed0為顆粒材料最小孔隙比.

正壓性函數(shù)fe和fd的具體表達式為

式中：ec為散粒體材料臨界孔隙比；ed為在等向壓縮時，顆粒材料在最大密度時對應的下限孔隙比；α和β為確定三軸試驗中內摩擦角峰值時所需要的指數(shù)項.

上述向密性函數(shù)fb，正壓性函數(shù)fe和fd均考慮了顆粒材料密度和平均應力.其3個特征孔隙比ei，ec，ed都隨著平均應力的增加而減小，并滿足下述表達式：

綜上，von Wolffersdorff亞塑性本構模型需要8個材料參數(shù)：φc，ec0，ed0，ei0，hs，α，β和n.

1.1.2 Niemunis和Herle修正的亞塑性本構模型

Niemunis和Herle修正的亞塑性本構模型[12]表達式為

為了描述某一狀態(tài)下顆粒間應變δ的狀態(tài)，引入?yún)?shù)ρ，其表達式為：

式中：R為顆粒間最大應變；δ代表DΔt積累程度的二階張量；‖‖為張量的歐幾里得范數(shù).

當ρ＝1時，

當0＞ρ＞1時，采用插值法確定顆粒材料剛度的一般表達式為

進而顆粒間應變的表達式為

綜上，為了在Niemunis和Herle修正的亞塑性本構模型[12]中引入顆粒間的應變的概念，需要在上述von Wolffersdorff亞塑性本構模型[11]中8個材料參數(shù)基礎上，補充5個材料參數(shù)：R，mR，mT，βr和χ.材料參數(shù)R為顆粒間應變上限值；mR和mT為在應變空間內，應變路徑分別出現(xiàn)90°和180°負載時對應的剛度強化系數(shù)，其關系為mR＞mT＞1；材料參數(shù)βr和χ用來表示剛度變化的程度，ρβr為考慮當前顆粒間應變的加權系數(shù).

1.2 數(shù)值模型方法及其合理性論證

本文采用ABAQUS中的UMAT子程序對鋼筋混凝土倉壁－散料相互作用進行數(shù)值模擬，并將數(shù)值模擬結果與既有試驗數(shù)據(jù)進行對比，從而驗證本文數(shù)值模型的合理性和數(shù)值穩(wěn)定性.

1.2.1 筒倉數(shù)值模型尺寸和材料本構關系模擬

在利用有限元軟件ABAQUS對鋼筋混凝土筒倉進行建模時，靜力作用下倉頂對倉壁散料壓力分布影響很小[13]，故忽略倉頂?shù)挠绊?筒倉高26.0 m，內徑28.0 m，壁厚0.3 m.散料頂點高度22.3 m，與倉壁相交高度21.3 m.利用ABAQUS中三維殼單元S4R模擬筒倉倉壁，利用三維實體單元C3D8模擬倉內散料.筒倉倉壁選用混凝土材料，其彈性模量E為3×1010Pa，泊松比μ為0.2，密度ρ為2 700 kg·m－3；筒倉內散料為小麥，其本構模型采用Niemunis和Herle修正的亞塑性本構模型[12]，其具體材料參數(shù)見表1和表2.

表1 定義亞塑性本構模型所需材料參數(shù)Tab.1 Material parameters of hypoplastic model

表2 基于顆粒間應變的附加材料參數(shù)Tab.2 Additional parameters of inter－granular strainconcept

1.2.2 筒倉壁－散料接觸關系

本文利用有限元軟件ABAQUS中主從接觸面算法模擬筒倉壁與散料之間的接觸關系.

筒倉壁與散料整體設定為剛柔接觸；筒倉壁與散料表面采用面面接觸；為了防止接觸分析工程中從接觸面侵入主接觸面，本文對主從接觸面之間的摩擦模型選用罰函數(shù)摩擦模型，并將摩擦系數(shù)按我國鋼筋混凝土筒倉設計規(guī)范[1]取為0.4；由于本文重點是研究筒倉壁－散料靜力相互作用，認為鋼筋混凝土筒倉倉壁始終處于彈性狀態(tài)，故本文中的數(shù)值模型不考慮幾何非線性.

1.2.3 數(shù)值模型合理性驗證

為了驗證本文所采用的亞塑性本構理論及所建立數(shù)值模型邊界條件的合理性，分別對筒倉實際結構和筒倉模型結構建立有限元數(shù)值模型，并將有限元數(shù)值模擬結果與陳長冰博士論文[15]中提供的筒倉實際結構試驗數(shù)據(jù)進行對比分析，對比結果見圖1.

圖1 數(shù)值模擬結果與筒倉實際結構試驗對比Fig.1 Com parison of numerical and practical results

由圖1可知，采用本文所建立的有限元數(shù)值模型獲得的結果與陳長冰[5]試驗結果基本一致，從而驗證了本文采用亞塑性本構理論模擬倉內散料的力學行為以及數(shù)值模型的合理性.

2 數(shù)值模擬結果與分析

2.1 數(shù)值模擬結果與各國規(guī)范比較

本文在進行筒倉倉內散料參數(shù)化分析時，為了更加直觀地反映各個材料參數(shù)對倉壁所受側向壓力的影響，將倉內散料頂面處理為平面，高度為26.0 m，其余參數(shù)不變.

本文采用表1和表2的材料參數(shù)研究小麥對鋼筋混凝土筒倉倉壁的側向壓力分布.在側向壓力分布上，楊鴻等[13]分析認為最大側向壓力分布在倉底；Juan等[14]指出底部壓力為負值；楊代恒[15]分析認為最大側壓力分布在距離倉底一段距離處；我國[1]、歐洲（ISO）[2]、美國[3]筒倉設計規(guī)范則表明最大壓力出現(xiàn)在底部.本文研究表明：最大側向壓力出現(xiàn)在筒倉倉底，這與文獻[13]和各國規(guī)范結果保持一致，如圖2所示.文獻[14－15]之所以出現(xiàn)倉壁所受最大側壓力未分布在倉底的結果，是因為其對倉底散料節(jié)點進行位移全約束，這種約束形式與實際情況的邊界條件不符.

圖2 數(shù)值模擬結果與各國規(guī)范對比Fig.2 Com parison of numerical results and codes

同時，從圖2可以看出，倉壁不同高度處各國規(guī)范值較數(shù)值模擬值（FEM）大8.7～23.5 kPa，這樣的安全儲備是各國規(guī)范綜合考慮各類不確定因素的結果；其中最小儲備值8.7 kPa約占倉底最大靜壓力值的11%左右.在筒倉靜壓力設計時，這樣的安全儲備是足夠的.

2.2 數(shù)值模擬結果與經(jīng)典筒倉理論比較

目前，比較經(jīng)典的筒倉壓力分布理論有Janssen理論、Rankine理論、修正的Coulomb理論、Airy理論等；現(xiàn)行的各國筒倉設計規(guī)范都以這些理論為基礎進行筒倉壓力設計.因此，本文將倉內散粒為小麥的數(shù)值模擬結果與這些經(jīng)典的筒倉壓力理論比較，對比結果見圖3.

由圖3可知，Janssen理論、Rankine理論及修正的Coulomb理論都將倉內散粒為小麥的數(shù)值模擬結果包絡在內，這說明上述3種經(jīng)典筒倉壓力分布理論具有較高的安全儲備；其中，利用Janssen理論所得到的筒倉壓力分布值安全儲備最高，這也是各國規(guī)范廣泛采用Janssen理論的一個重要原因.

2.3 不同初始孔隙比對倉壁側向壓力影響

在對倉內散料進行材料參數(shù)化分析前，本文假定倉內散料初始孔隙比保持不變.之所以假定筒倉內散料初始孔隙比保持不變，是因為散料通過傳送裝置運至倉頂后自由落體至倉底或落入倉內既有散料的頂面，直至整個筒倉裝料完畢；在整個裝料過程中，大量散料顆粒碰撞，使得倉內散料在裝料結束后密實程度很高，很難再進一步發(fā)生孔隙比的變化.

圖3 數(shù)值模擬結果與經(jīng)典筒倉壓力理論對比Fig.3 Comparison of numerical results and classical theories

本文所提到的初始孔隙比指的是裝料結束后，短時間內倉內散料所對應的空隙比。在隨后的分析中，假定該初始孔隙比保持不變，進而討論裝料后不同情況下散料的初始孔隙比對倉壁靜壓力的影響。根據(jù)筒倉裝料工作原理，文中假定初始孔隙比為常數(shù)，倉內散料的高度和密度不會隨時間變化而變化.

基于初始孔隙比不變這一基本假定，本文在基于修正的von Wolffersdorff亞塑性本構模型進行散料－筒倉靜力相互作用分析時，分別考慮顆粒材料初始孔隙比為0.40，0.45，0.50，0.55，0.60，0.65和0.70時對倉壁側壓的影響，分析結果見圖4.

圖4 初始孔隙比的影響Fig.4 Influence of initial void ratio

分析結果表明，顆粒材料初始孔隙比越大，其對筒倉倉壁的影響則越大.當散料初始孔隙比為0.40和0.45時，其對倉壁的壓力分布較為接近；隨著散料孔隙比的增加，其對筒壁側壓力值增幅明顯.之所以出現(xiàn)倉內散料初始孔隙比越大，倉壁所受靜壓力越大，是因為散料初始孔隙比越大，其向倉壁周圍的變形越大，筒倉倉壁對其的約束作用越強，倉壁所受靜壓力也越大.也就是說，同樣的散料，在初始孔隙比不同的情況下，散料對倉壁的側壓力分布形式也有明顯的差異.因此，在規(guī)范設計時，應當考慮倉內散料中初始孔隙比這個材料參數(shù)的影響.值得注意的是，顆粒材料并不是一直保持其初始孔隙比的大??；當顆粒材料在周圍環(huán)境振動的作用下，其孔隙比會隨著時間的變化而變化.圖4同時也反映出倉內散料在堆積過程中受到環(huán)境擾動，孔隙比的改變對倉壁側壓力分布的影響.散料孔隙比在擾動下會有變小的趨勢；正是由于這種趨勢，散料對筒倉壁的壓力也有減小的趨勢.然而，多數(shù)文獻在進行筒倉倉壁靜力分析時，認為在重力作用下，倉內顆粒材料孔隙比基本保持不變.這樣的假定忽略了散料對倉壁的側壓力是一個變化的過程，易高估筒倉倉壁承受散料對其側向壓力的能力.

2.4 散料與倉壁間摩擦系數(shù)對倉壁側向壓力影響

在實際工程中，散料與筒倉倉壁之間的摩擦系數(shù)一般集中在0.30～0.50范圍內.本文在保持散料其他物理特性參數(shù)不變的情況下，取散料與倉壁間摩擦系數(shù)分別為0.30，0.35，0.40，0.45和0.50進行分析，結果見圖5.

圖5 摩擦系數(shù)的影響Fig.5 Influence of friction coefficient

從圖5可以看出，倉壁在15～20 m高度范圍內，對應的5條倉壁側壓力分布曲線基本重合；倉壁在0～5 m高度范圍內，其對應的5條側壓力分布曲線有所差異，但相差約5 kPa.上述分析結果表明，不同摩擦系數(shù)對倉壁側向壓力影響很小.

2.5 散料內摩擦角對倉壁側向壓力的影響

實際工程中，倉內散料的內摩擦角一般在22°～40°范圍內.本文在保持散料其他物理特性參數(shù)不變的情況下，分別對散料內摩擦角為22°，25°，28°，31°，34°，37°和40°進行分析，結果見圖6.

圖6 散料內摩擦角的影響Fig.6 Influence of internal friction angle

從圖6可以看出，倉內散料的內摩擦角對筒倉壁側向壓力分布影響很大.隨著散料內摩擦角的增大，倉壁靜壓力在減小，這是因為倉內散料內摩擦角越大，散料顆粒間相互咬合力也越大，因此會減小散料顆粒對倉壁的法向變形，即倉壁對與其接觸的散料約束作用變小，倉壁所受靜壓力也變??；同時，在距離倉底約0.5 m高度范圍內，不同內摩擦角的散料對倉壁產(chǎn)生的側向應力差異非常明顯，最大時甚至相差20.3 kPa，其原因是倉內散料在重力作用下出現(xiàn)了較大的塑性變形.因此，進行筒倉倉壁內力設計時，應當充分考慮散料內摩擦角及塑性變形對筒倉壁靜力作用的影響.

2.6 不同顆粒粒徑的筒倉倉壁側向壓力分析

散料顆粒硬度是間接反映顆粒粒徑大小的材料參數(shù)，其范圍一般為1 000（細砂）～30 000 MPa（粗砂和砂礫）.本文在保持其他材料參數(shù)不變下，選取顆粒硬度分別為1 000，1 800，4 000，6 000，8 000，10 000和30 000 MPa進行分析，結果見圖7.

圖7 不同顆粒硬度的影響Fig.7 Influence of different granulate hardness

從圖7可以看出，散料顆粒粒徑越小，倉壁所受側向壓力越大；筒倉底部不同顆粒粒徑所對應的最大與最小側向壓力相差約12 k Pa.這說明散料顆粒粒徑的大小對筒倉倉壁所受側向壓力的影響較大.然而，在現(xiàn)行各國規(guī)范中，并未將散料的顆粒粒徑考慮到筒倉倉壁側向壓力的設計之中.

3 結論

本文基于修正的von Wolffersdorff亞塑性本構理論對筒倉倉壁與散料顆粒材料之間的靜態(tài)壓力作用進行有限元模擬.可得出以下結論：

（1）將筒倉－散料靜力相互作用的數(shù)值模擬結果與既有筒倉試驗數(shù)據(jù)、中歐美筒倉設計規(guī)范及經(jīng)典筒倉壓力理論進行對比，結果表明采用亞塑性本構理論模擬散料能夠較為準確地反映倉壁與散料相互作用的實際狀況.

（2）通過對倉內散料初始孔隙比、內摩擦角、摩擦系數(shù)、顆粒硬度和顆粒間應變進行參數(shù)分析，認為散料受初始孔隙比、內摩擦角、顆粒硬度和顆粒間應變影響較大，建議規(guī)范按照材料參數(shù)分類進行筒倉倉壁靜力壓力設計.

（3）通過對反映粒徑大小的顆粒硬度進行參數(shù)分析，分析結果表明顆粒粒徑大小對筒倉倉壁壓力影響較大，且顆粒粒徑越大，倉壁所受壓力越大；建議規(guī)范給出同類型散料顆粒粒徑最大值，并借助有限元軟件取最大粒徑值進行筒倉靜力壓力設計.

[1] 中華人民共和國建設部.GB 50077－2003鋼筋混凝土筒倉設計規(guī)范[S].北京：中國計劃出版社，2004.China Ministry of Construction.GB 50077－2003 Code for design of reinforced concrete silos [S].Beijing：China Planning Press，2004.

[2] European Committee for Standardization.Eurocode 8 Design of structures for earthquake resistance，part 4：silos，tanks and pipelines [S]. Brussels： European Committee for Standardization，2006.

[3] American Concrete Institute（ACI）.Standard practice for design and construction of concrete silos and stacking tubes for storing granular materials[S].Detroit：American Concrete Institute，1997.

[4] Carson W.Limits of silo design codes[J].Practice Periodical on Structural Design and Construction，2014，10（10）：61.

[5] 陳長冰.筒倉內散體側壓力沿倉壁分布研究 [D].合肥：合肥工業(yè)大學，2006.CHEN Changbing.Research on distribution of lateral bulksolid pressures on silos’walls[D].Hefei：Hefei University of Technology，2006.

[6] Holst J，Ooi J，Rotter J，et al.Numerical modeling of silo filling.Ⅰ：continuum analyses[J].Journal of Engineering Mechanics，1999，125（1）：94.

[7] Holst J，Ooi J，Rotter J，et al.Numerical modeling of silo filling. Ⅱ：discrete element analyses [J].Journal of Engineering Mechanics，1999，125（1）：104.

[8] WAN Xuewen，YANG Zhaojian，SHU Xuefeng，et al.The static contact statuses between granular materials and flatbottomed steel silos[J].Powder Technology，2013，235：1053.

[9] Guaital M，Couto A，Ayuga F.Numerical simulation of wall pressure during discharge of granular material from cylindrical silos with eccentric hoppers [J].Biosystems Engineering，2003，85（1）：101.

[10] Herle I，Gudehus G.Determinations of parameters of a hypoplastic constitutive model from properties of grain assemblies[J].Mechanics of Cohesive Materials，1999，4：461.

[11] Von Wolffersdorf f. A hypoplastic relation for granular materials with a predefined limit state surface[J].Mechanics of Cohesive Materials，1996，1（3）：251.

[12] Niemunis A，Herle I.Hypoplastic model for cohesive soils with elastic strain range [J].Mechanics of Cohesive Materials，1997，2（4）：279.

[13] 楊鴻，楊代恒，趙陽.鋼筒倉散料靜態(tài)壓力的三維有限元模擬[J].浙江大學學報：工學版，2011，45（8）：1423.YANG Hong， YANG Daiheng， ZHAO Yang. Threedimensional finite element simulation of static granular material pressure for steel silos[J].Journal of Zhejiang University：Engineering Science，2011，45（8）：1423.

[14] Juan A，Moran J M，Guerra M I，et al.Es
Tablishing stress state of cylindrical metal silos using finite element method：comparison with ENV 1993 [J].Thin Walled Structures，2006，44（11）：1192.

[15] 楊代恒.鋼筒倉設計研究及散料壓力的數(shù)值模擬 [D].杭州：浙江大學，2008.YANG Daiheng.Design research of steel silos and numerical simulation of bulk materials pressure[D].Hangzhou：Zhejiang University，2008.