偶應力對層狀巖體結構面邊界層效應的影響

張敦福，王相玉,，朱家明,，李術才，朱維申

（1.山東大學 土建與水利學院，濟南 250061；2.清華大學 工程力學系，北京 100084；3.中國科學院力學研究所，北京 100190）

1 引 言

巖體是自然地質體，內部存在斷層、層理、節理、裂隙等。巖石中兩個以上的塊體相互連接，荷載通過其連接界面傳遞而產生變形。由于接觸問題屬于力-位移混合邊界問題，采用混合有限元（雜交單元）方法具有優勢[1]。求解界面接觸問題有多種方法，如接觸單元法、Lagrange乘子法、罰方法等。完整巖塊在有限元計算中可用常規的實體單元代表。節理是層狀巖體中的不連續面，其厚度接近于0。對于重要節理，通常采用節理單元計算；對于不重要的節理，可采用等效單元計算。

經典連續介質理論認為，材料某點處應力僅是該點的應變以及該點的變形歷史上的函數，而與該點以外的其他點的應力無關[2]。但由于連續性假設不能嚴格滿足，因此，將連續介質力學應用于巖土介質時，應力和應變等分量代表的僅是相當小但不是無窮小體積上的統計平均值。

在應變梯度不大的情況下，使用統計平均值來替代連續介質力學理論解可以較為恰當地描述介質的力學反映。但當材料出現較高的應變梯度時，在相當小體積上，應變呈高階非線性，經典理論所代表的統計平均值就不能真實地反映出材料在相當小的體積上的強度和變形行為。實質上，梯度項的出現反映這樣一個事實：在某種尺度下的微結構相互作用使得變形是非局部的，應變梯度及內部長度描述的是不均質材料微結構之間的相互影響及作用。

偶應力理論在20世紀90年代初被引入巖土工程領域。Mühlhaus[3]的工作涉及層狀巖體、節理巖體、層狀半無限土體、巖煤和砂礫等。陳勝宏[4]、余成學[5]等率先開展研究，在這方面取得了開創性的成果。潘一山等[6]把微結構應變梯度理論應用于巖體力學，利用試驗研究了巖石破壞的尺度效應，又提出了巖石失穩破壞的應變梯度模型，并給出了節理巖體內部特征長度的數值。其他學者也對這一理論進行了不同方面的研究。偶應力理論的有限單元法要求位移解滿足C1連續，而一般單元只能滿足C0連續，陳萬吉[7]等很多學者在構造新的滿足 C1連續的單元方面做了許多工作。Taiji Adachi等[8]采用罰函數強制宏觀轉角與微觀轉角相同提出了一種有限單元，但由于計算結果不理想，只好對將組成一個四邊形的兩個三角形單元中的一個實施懲罰，以提高計算精度。Fleck等[9]采用罰函數法放松了對旋轉自由度的約束，取得了較好的效果。肖其林等[10]、李雷等[11]構造了非協調單元和雜交元，并用它分析了偶應力問題。楊樂等[12]、劉俊等[13]、冀賓等[14]將Cosserat介質元引入均勻化技術，建立層狀巖體的復合等效本構關系，進行有限元數值模擬，驗證了新模型的有效性和可行性。唐欣薇等[15]利用位移漸近展開技術和均勻化理論，建立了多尺度力學框架下的有限元平衡方程，重點考察了單胞尺寸對宏觀力學性能的影響。

本文基于偶應力理論，采用有限單元法，研究層狀巖體中界面邊界層效應問題，分析了偶應力的尺度效應、材料特征長度、彈性模量等參數對界面邊界層效應的影響。為了提高計算結果的精度，采用高次矩形單元中的Serendipity單元（S單元）。

2 有限元插值函數及位移插值

2.1 經典彈性理論

根據插值函數求法的不同，矩形單元可以分為Lagrange單元（L單元）和Serendipity單元（S單元）。由于這兩族單元的插值函數求法不同，所以單元節點的分布也不相同，即L單元的單元內部包含節點，而S單元的內部沒有節點，這也是Serendipity族單元的特點。由于Lagrange族高次單元的內部節點不影響單元間的連續性，它們能在單元范圍上取消，以便減少單元矩陣的階數，所以本文采用了圖1的S單元，以避免出現在L單元內部的中間節點。

圖1 S單元插值函數計算示意圖Fig.1 An illustration of interpolation function for S-element

S單元的節點插值函數具有下列性質：

式中：Ni為節點插值函數； i,j=1,2,...,8；δij為克羅內克爾符號。

對于節點1的插值函數N1，在節點2，3，…，8上應取0值，在節點1上取單位值。相應的N1在x/a-1=0，y/b-1=0和x/a+y/b+1=0邊界上為0。因此，假設N1為

類似地，可以得到N1，N2，…，N8的表達式。現將插值函數一并列出

由插值函數的定義可知，單元內任一點的位移可以這樣給出：

式中：ui、vi分別為節點i的位移分量；[N]=[N1?I N2?I…N8?I]（I為2階單位矩陣）為插值函數矩陣；{δ}e為單元e上8個節點的位移。

2.2 彈性偶應力理論

彈性偶應力理論與經典彈性理論的區別主要在于引入了偶應力mx、my，即單位面積上的彎矩。平面微元體有3個自由度：2個平動分量和1個轉動(旋轉)分量。

由于偶應力的存在，剪應力互等定理不再成立，即τxy≠τyx。為了便于分析，將剪應力τxy和τyx寫成對稱分量τs和反對稱分量τα兩部分。即

這樣就可以很清楚地描述τs和τα的作用。即，τs使微元體產生剪切變形，τα使微元體產生宏觀轉動。如圖2所示。

老機八從懷里掏出兩塊銀元遞給大勇：老子死了就幫俺捎給俺娘，跟俺娘說一聲：娘，俺不能盡孝了。對了，老子萬一活著呢，你得還把錢還給俺。

圖2 微元體的變形Fig.2 Deformation of micro element

由于偶應力理論中單元有3個自由度，所以位移插值模式與以前有所不同，寫成矩陣形式為

式中：u、v分別為節點i沿x軸和y軸方向的位移；ω為單元轉角；[N]=[N1?I N2?I…N8?I ]（I為3階單位矩陣）稱為插值函數矩陣，是3×24階矩陣；{δ}e為單元e上8個節點的位移，是24×1階矩陣。所以，應變矩陣 [B ]=[B1B2… B8]。其中

3 層狀巖體結構面邊界層效應

3.1 偶應力理論和經典彈性理論結果比較

由于偶應力理論中剪應變γxy≠γyx，剪應力τxy≠τyx，所以偶應力理論下剪應變用(γxy+γyx)/2與經典彈性理論下剪應變γxy相比較；剪應力用(τxy+τyx)/2與經典彈性理論下剪應力τxy相比較。

如圖3所示，層狀巖體的兩端受均布壓力q=10 MPa和剪力τ=1 MPa。材料1的彈性模量與泊松比分別為E1=20 GPa、μ1=0.25，材料2的彈性模量與泊松比分別為E2= 40 GPa，μ2=0.25。材料1、2的特征長度為 l1=l2= 1 mm，第二剪切模量為 GC=0.5G，G為剪切彈性模量。

圖3 層狀巖體Fig.3 Layered rock mass

兩種材料在節理面上滿足以下假設：兩種材料的界面黏合性好，無滑動、張開或嵌入。

定義圖4所示路徑，x坐標為0.025 m，y坐標為0.052～0.048 m。將該路徑各應力、應變值的經典彈性理論和偶應力理論曲線繪于圖5～10。

從圖5～10可知，偶應力使得正應變、剪應變、主應變的絕對值減小。剪應變尺度效應明顯，主應變和正應變尺度效應不明顯。剪應變在邊界層兩側出現了一個過渡區域，邊界層兩側變化趨勢一致。但剪應力在邊界層附近不再連續，邊界層兩側變化趨勢一致。

3.2 彈性模量E2/E1的影響

取E2分別為40、60 GPa，其他參數同3.1節，考察 E2/E1對剪應力和剪應變的影響。圖11給出了E2/ E1對剪應力、剪應變的影響曲線。

圖4 計算路線Fig.4 Calculating route

圖5 剪應力Fig.5 Shear stress

圖6 剪應變Fig.6 Shear strain

圖7 正應力σx及σyFig.7 Normal stresses σxand σy

圖8 正應變εx和εyFig.8 Normal strains εxand εy

圖9 主應力σ1和σ3Fig.9 Principal stresses σ1and σ3

圖10 主應變ε1和ε3Fig.10 Principal strains ε1and ε3

圖11 E2/E1對剪應力、剪應變的影響Fig.11 Influences of E2/E1on shear stress and shear strain

從圖11可以看出，彈性模量只對剪應力、剪應變的數值有影響，而不影響其變化趨勢，即不影響邊界層的范圍。剪應力不再連續。

3.3 第二剪切彈性模量GC的影響

3.4 泊松比的影響

兩種材料的泊松比μ依次取為0.18、0.25、0.32，其他參數與3.1節相同，考察μ對剪應力和剪應變的影響。圖13給出了泊松比對剪應力、剪應變的影響曲線。從圖可看出，泊松比不影響剪應力、剪應變過渡區域大小，不改變變化趨勢，只影響剪應力、剪應變的大小。剪應力不再連續。

3.5 特征長度的影響

兩種材料的特征長度l依次取為 0.5、1.0、2 mm，其他參數同3.1節，圖14給出了特征長度對剪應力、剪應變的影響曲線。

圖12 GC對剪應力及剪應變的影響Fig.12 Influences of GCon shear stress and shear strain

圖13 μ 對剪應力及剪應變的影響Fig.13 Influences of μ on shear stress and shear strain

圖14 l對剪應力及剪應變的影響Fig.14 Influences of l on shear stress and shear strain

從圖 14可見，考慮偶應力后，剪應力、剪應變在邊界附近出現了一個過渡區域。隨著特征長度減小，剪應力、剪應變的過渡區域減小，在過渡區內斜率變大。剪應力不再連續。

4 偶應力對錨桿界面剪應力的影響

巖體中的拉拔錨桿，受力為400 kN，尺寸如圖15所示。材料1為錨桿：直徑為0.12 m，彈性模量E1=210 GPa，泊松比μ1=0.25；材料2為巖體：彈性模量E2=40 GPa，μ2=0.14。材料1、2的特征長度分別為l1=0.1 mm、l2=1 mm；兩種材料的第二剪切模量取GC/G=0.5。

圖15 拉拔錨桿Fig.15 Drawing bolt

兩種材料在節理面上滿足以下假設：兩種材料的界面粘合性好、無嵌入。

圖16給出了錨桿界面剪應力分布。由圖可知，拉拔錨桿錨固段界面剪應力的分布集中在錨固段前部，呈先增大后衰減的趨勢，最大值點到端頭的距離約等于錨桿直徑的一倍，這些結論與文獻[16－19]一致。圖17為沿x軸的剪應變及剪應力分布。由圖可知，考慮偶應力后，界面上剪應力分布沒有明顯變化，偶應力對剪應力峰值有明顯影響。剪應力、剪應變在界面附近出現了一個過渡區域，并且剪應力不再連續。

圖16 錨桿界面剪應力Fig.16 Shear stress in the interface between bolt and rock

圖17 沿x軸的剪應變及剪應力Fig.17 Shear stress and shear strain along x-axis

5 結 論

（1）在邊界層內，存在一定的轉動梯度效應，偶應力不為0，界面處剪應力不再連續。

（2）偶應力對層狀巖體結構面剪應變有顯著影響，尺度效應明顯。在偶應力理論下，應變量減小，并在邊界層附近出現過渡區，剪應變突變現象得以改善，但剪應力不再連續。

（3）材料特征長度對剪應變和剪應力的過渡區域大小有影響。隨著特征長度的減小，過渡區域減小，剪應變和剪應力在過渡區域內斜率增大，即應變梯度增大，尺度效應更加明顯。

（4）第二剪切模量 GC對剪應變和剪應力的過渡區域大小沒有影響。隨著 GC增大，剪應變和剪應力在過渡區域內斜率增大，即應變梯度增大，尺度效應更加明顯。

（5）泊松比、彈性模量對剪應變和剪應力的過渡區域大小沒有影響，只影響界面處剪應變和剪應力的大小，不改變剪應變和剪應力的變化趨勢。

（6）錨桿界面剪應力的分布集中在錨固段前部，呈先增大后衰減的趨勢，最大值點到端頭的距離約等于錨桿直徑的一倍。偶應力不改變錨桿界面上的剪應力分布，但對峰值有明顯的影響，峰值變小，說明錨桿的抵抗外載的能力變大。

[1]張楚漢. 論巖石、混凝土離散-接觸-斷裂分析[J]. 巖石力學與工程學報,2008,27(2): 217－235.ZHANG Chu-han. Discrete-contact-fracture analysis of rock and concrete[J]. Chinese Journal of Rock Mechanics and Engineering,2008,27(2): 217－235.

[2]趙冰,李寧,盛國剛,等. 應變梯度理論在巖土力學中的進展述評[J]. 長沙交通學院學報,2005,21(1): 47－52.ZHAO Bing,LI Ning,SHENG Guo-gang,et al. A survey of strain gradient applied in rock and soil mechanics[J].Journal of Changsha Communications University,2005,21(1): 47－52.

[3]MüHLHAUS H. B. Stress and couple stress in a layered half plane with surface loading[J]. International Journal for Numerical and Analytical Methods in Geomechanics,1990,14(8): 545－563.

[4]陳勝宏,王鴻儒,熊文林. 節理巖體的偶應力的影響研究[J]. 水利學報,1990,(1): 44－48.CHEN Sheng-hong,WANG Hong-ru,XIONG Wen-lin.The effect of couple stress on jointed rock mass[J].Journal of Hydraulic Engineering,1990,(1): 44－48.

[5]余成學,熊文林,陳勝宏. 層狀巖體的彈黏塑性Cosserat介質理論及其工程應用[J]. 水利學報,1996,(4): 10－17.SHE Cheng-xue,XIONG Wen-lin,CHEN Sheng-hong.Elastovisco-plastic Cosserat theory of layered rockmass and its application to engineering[J]. Journal of Hydraulic Engineering,1996,(4): 10－17.

[6]潘一山,袁旭東,章夢濤. 巖石失穩破壞的應變梯度模型[J]. 巖石力學與工程學報,1998,17(增刊): 760－765.PAN Yi-shan,YUAN Xu-dong,ZHANG Meng-tao.Model of gradient strain about failure rockmass with instability[J]. Chinese Journal of Rock Mechanics and Engineering,1998,17(Supp.): 760－765.

[7]陳萬吉. 應變梯度理論有限元 C0-1分片檢驗及其變分基礎有限元增強型分片檢驗[J]. 大連理工大學學報,2004,44(4): 474－477.CHEN Wan-ji. Finite element methods in strain gradient theory: C0-1patch test and its variational basic[J]. Journal of Dalian University of Technology,2004,44(4): 474－477.

[8]ADACHI T,TOMTA Y,TANAKA M. Computational simulation of deformation behavior of 2D lattice continuum[J]. International Journal of Mechanical Sciences,1998,40(9): 857－866.

[9]FLECK N A,SHU J. Microbuckle initiation in fibre composites: A finite element study[J]. Journal of the Mechanics and Physics of Solids,1995,43(12): 1887－1918.

[10]肖其林,凌中,吳永禮. 偶應力問題的雜交-混合元分析[J]. 計算力學學報,2003,20(4): 427－433.XIAO Qi-lin,LING Zhong,WU Yong-li. Hybrid-mixed finite element analysis of couple stress problems[J].Chinese Journal of Computational Mechanics,2003,20(4): 427－433.

[11]李雷,吳長春. 基于Cosserat理論的應變梯度非協調數值研究[J]. 工程力學,2004,21(5): 166－171.LI Lei,WU Chang-chun. Incompatible finite element for materials with strain gradient effects[J]. Engineering Mechanics,2004,21(5): 166－171.

[12]楊樂,吳德倫,李德萬,等. 層狀巖體的 Cosserat介質均勻化[J]. 應用力學學報,2010,27(1): 96－100.YANG Le,WU De-lun,LI De-wan,et al. Homogenization method for layered rock mass as Cosserat medium[J]. Chinese Journal of Applied Mechanics,2010,27(1): 96－100.

[13]劉俊,黃銘. 考慮偶應力影響的雙組節理巖體等效模型[J]. 長江科學院院報,2010,27(9): 43－46.LIU Jun,HUANG Ming. Equivalent model of rock-masses with two sets of joints considering influence of couple stresses[J]. Journal of Yangtze River Scientific Research Institute,2010,27(9): 43－46.

[14]冀賓,陳萬吉,趙杰. 基于偶應力理論剪切帶問題的彈塑性有限元分析[J]. 力學學報,2009,41(2): 192－199.JI Bin,CHEN Wan-ji,ZHAO Jie. Elastoplastic finite element analysis of shear band with couple stress theory[J]. Chinese Journal of Theoretical and Applied Mechanics,2009,41(2): 192－199.

[15]唐欣薇,張楚漢. 基于均勻化理論的混凝土宏細觀力學特性研究[J]. 計算力學學報,2009,26(6): 876－881.TANG Xin-wei,ZHANG Chu-han. Study of concrete in macro and meso scale mechanical properties based on homogenization theory[J]. Chinese Journal of Computational Mechanics,2009,26(6): 876－881.

[16]趙同彬,尹延春,譚云亮,等. 錨桿界面力學試驗及剪應力傳遞規律細觀模擬分析[J]. 采礦與安全工程學報,2011,28(2): 220－224.ZHAO Tong-bin,YIN Yan-chun,TAN Yun-liang,et al.Mechanical test of bolt interface and microscopic simulation of transfer law for shear stress[J]. Journal of Mining & Safety Engineering,2011,28(2): 220－224.

[17]賀若蘭,張平,李寧,等. 拉拔工況下全長粘結錨桿工作機理[J]. 中南大學學報,2006,37(2): 401－407.HE Ruo-lan,ZHANG Ping,LI Ning,et al. Working mechanism of fully grouted bolt in pull-out working state[J]. Journal of Central South University,2006,37(2): 401－407.

[18]陳浩,任偉中,李丹,等. 深埋隧道錨桿支護作用的數值模擬與模型試驗研究[J]. 巖土力學,2011,32(增刊1):719－786.CHEN Hao,REN Wei-zhong,LI Dan,et al. Numerical simulation and model test study of mechanism of bolt in deep tunnel[J]. Rock and Soil Mechanics,2011,32(Supp.1): 719－786.

[19]張幼振,石智軍,張晶. 巖石錨桿錨固段荷載分布試驗研究[J]. 巖土力學,2010,31(增刊2): 184－188.ZHANG You-zhen,SHI Zhi-jun,ZHANG Jing.Experimental study of load distribution of anchoring section for rock bolts[J]. Rock and Soil Mechanics,2010,31(Supp. 2): 184－188.