基于模擬退火分數(shù)階導數(shù)伯格斯模型的人工凍土蠕變研究

周仁戰(zhàn),姚兆明

1.安徽理工大學土木建筑學院，安徽 淮南 232001；2.淮南市建筑管理處，安徽 淮南 232008)

基于模擬退火分數(shù)階導數(shù)伯格斯模型的人工凍土蠕變研究

周仁戰(zhàn)1,2,姚兆明1

1.安徽理工大學土木建筑學院，安徽 淮南 232001；2.淮南市建筑管理處，安徽 淮南 232008)

準確掌握人工凍土的蠕變特性對控制凍結法施工的安全具有重要意義。目前常用的凍土蠕變本構模型是經(jīng)驗模型和整數(shù)階元件模型，模型中的參數(shù)大多缺乏明確的物理意義或者參數(shù)繁多。以伯格斯蠕變模型為基礎，建立了參數(shù)簡單且物理意義明確的人工凍土分數(shù)階導數(shù)蠕變方程，并給出了模擬退火算法優(yōu)化模型參數(shù)的過程?；谑覂热渥冊囼灲Y果，分析了溫度、加載應力和土質對人工凍土單軸蠕變特性的影響。并將試驗結果與本次建立的蠕變模型計算結果進行比較發(fā)現(xiàn)，兩者吻合良好。研究成果對目前國內外煤礦立井表土層施工中凍土爆破方面的研究、凍結管參數(shù)的選取、凍結壁壁厚及凍結溫度的設計具有重要的基礎作用。

人工凍土蠕變；模擬退火；分數(shù)階導數(shù)；伯格斯模型

凍土的蠕變一直是凍土力學研究和凍土工程設計施工所關注的問題，隨著寒區(qū)工農(nóng)業(yè)的發(fā)展以及凍結施工技術在工程上的推廣應用,對凍土強度及蠕變的研究力度越來越大，尤其是人工凍土蠕變。因此準確建立人工凍土蠕變模型，明確其本構方程中相關參數(shù)的物理意義，對目前國內外煤礦立井表土層施工中凍土爆破方面的研究、凍結管參數(shù)的選取、凍結壁壁厚及凍結溫度的設計具有重要的作用，也對其它工程領域中人工凍土蠕變的預測、工程施工方案的設計和施工過程中的安全隱患排除具有極強的現(xiàn)實指導意義。

自 20 世紀 60 年代以來，凍土的流變本構理論研究取得了較大的發(fā)展。其中經(jīng)驗模型和整數(shù)階元件模型是目前常用的凍土蠕變本構模型。經(jīng)驗模型通過對人工凍土蠕變試驗數(shù)據(jù)進行最小二乘法擬合,得到應變與時間的函數(shù)關系，該模型公式簡單，針對性強,可以很好地對蠕變曲線進行擬合,但模型參數(shù)缺乏明確的物理意義[1]。整數(shù)階元件模型是由彈簧元件、摩擦元件和阻尼器元件按各種組合而得。整數(shù)階元件蠕變模型的構成具有明確的物理力學意義，且直觀易懂[3]。整數(shù)階微積分能描述Euclid空間里的巖土體蠕變本構模型。整數(shù)階伯格斯元件模型、賓厄姆元件模型及西原元件模型等是常用的人工凍土經(jīng)典蠕變模型[2]。然而，模擬一個完整的人工凍土蠕變模型，往往需要很多元件進行組合，進而導致所建模型參數(shù)繁多。

本構模型參數(shù)辨識關系到人工凍土蠕變本構模型計算的精度。本構模型參數(shù)識別的過程是首先利用實驗數(shù)據(jù)建立合適的本構模型；然后用優(yōu)化算法回歸分析得到本構模型中待確定的參數(shù)。由此可見，參數(shù)的回歸結果取決于優(yōu)化算法的選擇。隨著人工智能的發(fā)展,很多學者將模擬退火算法應用到巖土本構模型和參數(shù)辨識中,并取得了豐富的成果[4-5]。模擬退火法是根據(jù)物理學中固體物質的退火過程與一般組合優(yōu)化問題的相似性，模擬高溫金屬降溫的熱力學過程。模擬退火方法在初始溫度很高、溫度下降緩慢的條件下，能以概率1收斂到全局最優(yōu)值，從而使其具有全局優(yōu)化的能力[6]。

為此，本文將巖土蠕變模型建立在伯格斯模型基礎上，用分數(shù)單元代替伯格斯模型中串聯(lián)的黏壺，建立分數(shù)階伯格斯模型，采用模擬退火算法優(yōu)化模型參數(shù)來模擬人工凍土單軸蠕變特性，得到的模型不僅具有經(jīng)典元件蠕變模型的優(yōu)點,即模型參數(shù)有明確物理力學意義，而且使用較少的串聯(lián)、并聯(lián)元件就能達到很好的模擬效果[7-10]。

1 人工凍土分數(shù)階導數(shù)伯格斯蠕變模型

1.1 伯格斯模型

伯格斯模型是目前應用較廣泛的整數(shù)階元件模型，模型能反映材料的彈性、黏彈性和彈塑性等力學行為，模型理論推導較成熟。馬克斯威爾體和開爾文體串聯(lián)而成的伯格斯模型如圖1所示。

圖1 伯格斯模型

伯格斯模型蠕變方程為[3]

(1)

式中,EB表示伯格斯模型中串聯(lián)彈簧的彈性模量；EK表示伯格斯模型中并聯(lián)彈簧的彈性模量；ηB表示伯格斯模型中串聯(lián)黏壺的黏滯系數(shù)；ηK表示伯格斯模型中并聯(lián)黏壺的黏滯系數(shù)；t表示時間，t=0時，只有彈簧起作用，e=sc/EB。t→∞時，e∞→∞但其變化速率為常值。

1.2 分數(shù)階導數(shù)伯格斯模型

1) 分數(shù)階導數(shù)

分數(shù)階導數(shù)是一個研究任意階次導數(shù)算子特性及應用的數(shù)學問題。分數(shù)階導數(shù)就是把導數(shù)的階次推廣到分數(shù)甚至復數(shù)領域，但因其計算復雜且物理意義不明確，長期以來只是進行數(shù)學上的理論研究。B.B.Mandelbrot在1982年指出分數(shù)維存在于自然界和許多技術科學中，自此，分數(shù)階導數(shù)才獲得了快速的發(fā)展[11-12]。

函數(shù)f(t)的r階分數(shù)階導數(shù)定義為[13-14]

(2)

Ir(t)=0,t≤0

(3)

(4)

G(1+z)=zG(z),Re(z)>0

(5)

式中，Dr為整數(shù)階導數(shù)算子的推廣，定義為分數(shù)階導數(shù)；0≤r<1，G為Gamma 函數(shù)。Abel核Ir(t)與函數(shù)f(t)的廣義Stieltjes卷積就是函數(shù)f(t)的分數(shù)階導數(shù)。將應力與應變的分數(shù)階導數(shù)成正比定義為Abel黏壺的本構關系，則Abel黏壺本構關系如下式所示

s=ηDre(t)=ηIr(t)*de(t)=

(6)

蠕變柔量

(7)

(8)

Drf(t)的Laplace變換為

(9)

(m-1

若m=1，并設初始條件為零，則

L[Drf(t)=srF(s)-D-(1-r)f(0)=srF(s)

(10)

2) 分數(shù)階導數(shù)蠕變元件

用Abel體代替在經(jīng)典模型理論中的Newton黏性流體就形成了分數(shù)階導數(shù)蠕變模型[15]，Abel體元件如圖2所示。

圖2 Abel體元件

Abel黏壺的本構方程為

(11)

由式(11)可見，當r=1，Abel體就成為牛頓體元件，能反映理想流體的力學性質；當r=0，Abel體就成為是彈簧元件，能反映理想固體的力學性質?？梢?，當r∈(0,1)時，Abel體能反映處于理想流體和理想固體之間的材料力學性質，人工凍土蠕變特性介于理想流體和理想固體間的某種特性。

當s=sc常數(shù)時，材料的蠕變特性就能用Abel體元件來描述。對式(11)進行分數(shù)階積分，依據(jù)Riemann-Liouville型分數(shù)階微積分算子理論，可以得出[15]

(12)

3) 分數(shù)階導數(shù)伯格斯模型蠕變方程

圖3元件組合為分數(shù)階導數(shù)伯格斯模型。

圖3 分數(shù)階導數(shù)伯格斯模型

分數(shù)階導數(shù)伯格斯模型的本構方程為

(13)

開爾文模型本構關系為

(14)

Abel黏壺本構關系為

(15)

彈簧的本構關系為

(16)

根據(jù)模型元件的串聯(lián)關系，分數(shù)階伯格斯蠕變模型的總應變?yōu)?/p>

e(t)=e1(t)+e2(t)+e3(t)

(17)

分數(shù)階導數(shù)伯格斯模型蠕變方程為

(18)

式中,模型中串聯(lián)彈簧的彈性模量用EB表示；分數(shù)階Newton體黏滯系數(shù)用ξ表示；分數(shù)導數(shù)階數(shù)用r表示；并聯(lián)彈簧的彈性模量用EK表示；并聯(lián)黏壺的黏滯系數(shù)用ηK表示；t表示時間。

1.3 模擬退火算法優(yōu)化模型參數(shù)

模擬退火算法，是一種直接的非線性反演，其獨特的優(yōu)點是能夠避免反演陷入目標函數(shù)的局部極小。算法的思想來源于模擬固體退火降溫的過程。在金屬熱加工工藝中，退火是指金屬熱處理過程，是將金屬材料加熱到某一高溫狀態(tài)，然后讓其慢慢冷卻下來的過程。從統(tǒng)計熱力學的觀點來說，隨著溫度的降低，物質的能量將逐漸趨近于一個較低的狀態(tài)，并最終達到某種平衡。模擬退火算法就是在金屬退火機理基礎上建立起來的一種全局最優(yōu)化方法，能夠以隨機搜索技術從概率的意義上找出目標函數(shù)的全局最小點。

模擬退火算法優(yōu)化模型參數(shù)過程如下[10]

1)初始化退火溫度Tk(令k=0)，產(chǎn)生隨機初始解集{x0}。在優(yōu)化Burgers模型參數(shù)時

{x0}=[EK0,ηK0,EB0,ηB0] (19)

在優(yōu)化分數(shù)階伯格斯模型參數(shù)時

{x0}=[EK0,ηK0,EB0,r0,η0] (20)

2)在溫度Tk反復操作，一直達到溫度Tk的平衡值，在解x的范圍內產(chǎn)生新的可行解x′；計算x′的目標函數(shù)f(x′)和x的目標函數(shù)f(x)的ηf差值；按概率min{1,exp(-ηf/Tk)}>random[0,1]接收x′，其中random[0,1]是[0,1]區(qū)間內的隨機數(shù)。

3)退火操作。Tk+1=CKk,k→k+1，其中C(0,1)。當達到假定的收斂數(shù)據(jù)，則完成退火過程；相反，轉到步驟2。最優(yōu)值的優(yōu)化方向由退火溫度控制，同時劣質解由概率exp(-ηf/Tk)來接收，因此局部極值點可以通過此種算法得以避免。初始溫度只要設置足夠高，設置足夠慢的退火過程，全局最優(yōu)解就能得到。

2 試驗研究

2.1 人工凍土單軸抗壓、蠕變試驗

人工凍土單軸抗壓、蠕變試驗土樣為黏土、砂質黏土。砂質黏土和黏土取樣深度分別在180～200m和210～230m。凍土單軸試驗在WDT-100型微機控制電液伺服試驗機上進行。黏土、砂質黏土-5℃、-10℃、-15℃單軸強度試驗值見表1。在三個溫度水平(-5℃、-10℃、-15℃)進行蠕變試驗，試驗取3級荷載sc，分別為sc=0.3ss、sc=0.5ss、sc=0.7ss,ss為人工凍土單軸抗壓強度，蠕變試驗結果見表2。

表1 凍土單軸抗壓強度

表2 人工凍土單軸蠕變試驗

2.2 人工凍土蠕變影響因素分析

1) 溫度對人工凍土蠕變特性的影響

從圖4中可以看出，凍結砂質黏土在加載系數(shù)為0.5時，在不同溫度條件下均在較短時間內蠕變達到穩(wěn)定，但隨著養(yǎng)護溫度降低，凍土的蠕變穩(wěn)定值明顯下降，-5℃、-10℃、-15℃蠕變穩(wěn)定時應變值分別為7.48%、6.92%、5.70%；黏土在加載系數(shù)為0.7時，不同溫度條件下均出現(xiàn)非穩(wěn)定破壞，當溫度低被破壞時應變值較小。表明溫度越低凍土蠕變的應變越小。

圖4 人工凍土蠕變特性受溫度的影響

2) 加載應力對人工凍土蠕變特性的影響

圖5 人工凍土蠕變特性受加載系數(shù)的影響

從圖5中可以看出兩種土質在小加載系數(shù)時短時間內就達到穩(wěn)定，但隨著加載系數(shù)加大，達到蠕變穩(wěn)定時間相應出現(xiàn)滯后，同時蠕變穩(wěn)定時應變值加大。如圖5(a)所示砂質黏土在-5℃養(yǎng)護溫度時，在加載系數(shù)為0.3條件下，在1h左右即可達到穩(wěn)定應變3.10%；而在加載系數(shù)為0.5條件下，需在5h左右才達到穩(wěn)定應變7.48%。在大加載系數(shù)條件下，表現(xiàn)出蠕變非穩(wěn)定加速。

3) 土質對人工凍土蠕變特性的影響

圖6 人工凍土蠕變特性受土質的影響

從圖6中可以看出不同溫度條件下凍結砂質黏土與凍結黏土蠕變變化規(guī)律一致，但黏土整體蠕變應變值要大于砂質黏土蠕變應變值，-10℃養(yǎng)護溫度加載系數(shù)為0.5條件下，黏土蠕變應變值為8.81%，而砂質黏土蠕變應變值僅為6.92%，差值達到1.83%。

2.3 人工凍土蠕變模型模擬結果分析

伯格斯模型的參數(shù)和分數(shù)階導數(shù)伯格斯模型參數(shù)均采用模擬退火算法來優(yōu)化。模擬退火最大內部循環(huán)數(shù)取90，冷卻系數(shù)取0.99，能量轉化公式為Botzmann函數(shù)，能量變換為單步變換。黏土、砂質黏土在-5℃、-10℃、-15℃下，加載值分別為0.3ss、0.5ss、0.7ss。伯格斯蠕變模型參數(shù)見表3。

表3 伯格斯模型參數(shù)

伯格斯模型計算值和單軸蠕變試驗值對比見圖7，文中只畫出部分伯格斯模型計算值和試驗值的關系。

圖7 蠕變試驗數(shù)值和Burgers模型計算值的關系

從圖7可以看出，在擬合曲線的拐點附近，Burgers模型與試驗數(shù)據(jù)偏差較大。隨著加載時間的延長，這種偏差越來越大?？梢?，在足夠長的時間里，Burgers模型將與試驗數(shù)據(jù)嚴重偏離。

黏土、砂質黏土在溫度-5℃、-10℃、-15℃和加載系數(shù)0.3ss、0.5ss、0.7ss下的分數(shù)階導數(shù)伯格斯蠕變模型參數(shù)見表4。

表4 分數(shù)階導數(shù)伯格斯模型參數(shù)

分數(shù)階導數(shù)伯格斯模型計算值和單軸蠕變試驗值對比如圖8所示。

由圖8可知，含分數(shù)算子的伯格斯蠕變模型更能精確地擬合試驗數(shù)據(jù)，相關系數(shù)均在0.99左右。這是由分數(shù)算子的核函數(shù)為冪律函數(shù)，且其指數(shù)可為任意實數(shù)時有更大適應范圍所決定的。

圖8 蠕變試驗值與分數(shù)階Burgers模型計算值關系

3 結論

以伯格斯蠕變模型為基礎，建立了人工凍土分數(shù)階導數(shù)蠕變方程，并給出了模擬退火算法優(yōu)化模型參數(shù)的過程；結合不同溫度、不同加載水平的凍土蠕變試驗實例得到如下結論：

(1)溫度越低凍土蠕變的應變越??；不同溫度條件下凍結砂質黏土與凍結黏土蠕變變化規(guī)律一致，但黏土整體蠕變應變值要大于砂質黏土蠕變應變值；隨著加載系數(shù)增大，凍土蠕變越難出現(xiàn)穩(wěn)定加速過程，且在大加載系數(shù)條件下，只表現(xiàn)出蠕變非穩(wěn)定加速。

(2)在分析人工凍土單軸蠕變受溫度、土質及加載應力對土質影響規(guī)律的基礎上，利用分數(shù)階導數(shù)和整數(shù)階元件模型建立了分數(shù)階伯格斯蠕變模型模擬人工凍土單軸蠕變特性，該模型具有物理意義明確及參數(shù)少的優(yōu)點，是模擬人工凍土單軸蠕變特性計算的一種新方法。

(3)分數(shù)階伯格斯模型由分數(shù)階導數(shù)元件和整數(shù)階導數(shù)元件結合而成，利用模擬模擬退火法來優(yōu)化蠕變方程參數(shù)。將人工凍土蠕變模型計算結果和試驗結果進行比較發(fā)現(xiàn)，兩者吻合良好，相關系數(shù)均在0.99左右，能對人工凍土的蠕變規(guī)律進行很好的模擬。

[1] 尹曉文,傅 強,馬昆林.凍土三軸蠕變非線性數(shù)學模型研究[J].冰川凍土,2013,35(1):171-176.

[2] 李強,王奎華,謝康和.人工凍土流變模型的識別與參數(shù)反演[J]. 巖石力學與工程學報, 2004, 23(11): 1 895-1 899.

[3] 康永剛,張秀娥.巖石蠕變的非定常分數(shù)伯格斯模型[J].巖土力學, 2011,32(11): 3 237-3 241.

[4] 陳炳瑞,馮夏庭,丁秀麗,等.基于模式搜索的巖石流變模型參數(shù)識別[J].巖石力學與工程學報,2005,24(2):207-211.

[5] 朱珍德,徐衛(wèi)亞. 巖體黏彈性本構模型辨識及其工程應用[J].巖石力學與工程學報,2002,21(11):1 605-1 609.

[6] 高鷹,謝勝利.基于模擬退火的粒子群優(yōu)化算法[J].計算機工程與應用,2004,40(1)47.

[7] SHI MIU NOBUYUKI, ZHANG WEI. Fractional calculus approach to dynamic problem of viscoelastic material[J].JSME, 1999,42(4):827-830.

[8] KOELLER R. C. Applications of fractional calculus to the theory of viscoelasticity[J]. Journal of Applied Mechanics, 1984,51(2):299-307.

[9] OLDHAM K. B., SPANIER J.. The fractional calculus[M]. New York: Academic Press, 1974.

[11] SCHIESSEL H, METZKER R, BLUMEN A, et al. Generalized viscoelastic models: their fractional equations with solutions[J]. Journal of Physics A: Mathematical and General, 1995, 28(23): 6 567-6 584.

[11] SCHIESSEL H, BLUMEN A. Mesoscopic pictures of the sol-gel transition: ladder models and fractal networks[J]. Macromolecule, 1995, 28(11): 4 013- 4 019.

[12] HEYMANS N, BAUWENS J C. Fractal rheological models and fractional differential equations for viscoelastic behavior[J]. Rheologica Acta, 1994, 33(3):210-219.

[13] 孫海忠,張衛(wèi).一種分析軟土黏彈性的分數(shù)階導數(shù)開爾文模型[J].巖土力學,2007, 28(9):1 983- 1 986.

[14] MANDELBROT B B. The fractal geometry of nature[M]. San Francisco:Freeman,1982.

[15] 殷德順,任俊娟,和成亮等.一種新的巖土流變模型元件[J].巖石力學與工程學報,2007,26(9):1 899-1 903.

(責任編輯：李 麗，編輯：丁 寒)

Study of Artificial Frozen Soil Creep Based on Simulated Annealing Fractional Order Derivative Burgers Creep Model

ZHOU Ren-zhan1,2, YAO Zhao-ming1

(1. School of Civil Engineering and Architecture, Anhui University of Science and Technology, Huainan Anhui 232001, China; 2. Department of Architecture Administration of Huainan City, Huainan Anhui 232008, China)

A proper estimation of creep characteristic of artificial frozen soil is of great significance in controlling the safety of freezing method construction. The most commonly used model for frozen soil creep constitutive model is the empirical and integer order component model, but these model parameters are mostly lack of clear physical meaning or there are too many parameters. This article was based on the Burgers creep model in simulating annealing optimization, and an artificial frozen soil fractional derivative creep equation was established. The equation was of simple parameter and clear physical meaning. Furthermore, the optimizing process of the model parameters by using simulated annealing algorithm was given. Based on the laboratory creep test results, the influence of temperature, loading stress and soil property on uniaxial creep property of artificial frozen soil was analyzed. Compared the simulation results of artificial frozen soil creep model with test results, it can be found that the two results are in good agreement. The research results play a significant role in the study of frozen soil blast, selection of freezing pipe parameters, design of frozen wall thickness and freezing temperature?in vertical shaft topsoil construction of coal mine shaft in the current domestic and foreign research.

artificial frozen soil creep; simulated annealing;fractional order derivative; burgers model

2016-05-23

周仁戰(zhàn)(1982-)，男，安徽碭山人，高級工程師，博士，研究方向：巖土工程、地下工程。

TU443

A

1672-1098(2016)06-0030-09