含分數階導數元件非線性蠕變模型的二次開發

張 濤，何利軍

（南昌航空大學土木與建筑學院，江西 南昌330063）

含分數階導數元件非線性蠕變模型的二次開發

張 濤，何利軍

（南昌航空大學土木與建筑學院，江西 南昌330063）

通過GDS應力三軸儀對湛江黏土在恒定圍壓不同偏應力下的三軸蠕變試驗，得到幾組不同偏應力下軸向應變與時間關系數據，鑒于傳統線性元件組合模型無法精確描述黏土蠕變的非線性問題，提出采用分數階導數元件分別替換傳統Burger模型中兩個牛頓粘壺線性元件，建立起含兩個分數階導數元件的非線性蠕變模型，采用C++語言將該模型在FLAC3D中進行二次開發，通過建立一個三軸數值算例，將不同偏應力下數值計算結果與試驗數據進行對比，結果表明提出的非線性蠕變模型相比Burger模型更適合描述湛江黏土的蠕變過程。

分數階導數；二次開發；蠕變

巖土的蠕變問題一直是工程基礎處理、邊坡防護和圍巖開挖與支護過程中需要重視的棘手難題，而軟黏土在工程實踐中屬于地質不良的土體，其在我國沿海和內陸江河沖積流域分布較為廣泛，給工程建設帶來諸多潛在危害。

軟土的蠕變是巖土流變學研究的重點內容之一，許多研究人員依據大量的試驗和工程實踐得到了一些可觀的成果：如楊愛武[1]對天津濱海吹填軟土進行三軸不固結不排水和單軸壓縮蠕變試驗，結果表明兩類試驗歷時曲線都具有非線性特征；雷華陽[2]對濱海軟土進行不同加荷方式下的蠕變試驗，得到應變與應力和時間關系，結果表明濱海軟土具有明顯的非線性特征且其蠕變變形受其結構制約和影響。這些研究成果使得人們對土體蠕變的認識不斷加深，在本構模型建立方面也逐漸由單純的線性理論模型試圖描述土體蠕變特性，逐漸發展為半經驗半理論模型和非線性理論模型。

傳統蠕變模型都是通過粘、彈、塑三類線性元件以串并聯方式組合得到，用以描述對應土體的粘彈、粘塑、彈塑等特性，但這些模型具有一定的局限性，只適用于描述線性蠕變過程，而實際軟土蠕變具有很強的非線性特性。本文結合前人建立非線性蠕變本構模型研究方法和在FLAC3D中實現二次開發經驗，提出采用分數階導數元件（阿貝爾元件）替換傳統Burger模型中的牛頓粘壺，建立起變參數流變本構模型，并在FLAC3D中完成該模型的二次開發工作，通過建立模擬三軸蠕變數值算例，比較不同偏應力下數值結果與試驗數據的吻合程度。

湛江黏土是一種具有一定典型性的強結構性黏土[3]，前期通過GDS應力三軸儀對湛江黏土在恒定圍壓不同偏應力下的三軸蠕變試驗，得到幾組不同偏應力下軸向應變與時間關系數據[4]，采用各種模型進行了描述[5]，由于傳統模型都是線性模型，基于分數階軟體元件構建非線性模型[6]，在理論上具有更好的描述效果。本文在FLAC3D中基于C++進行二次開發，相對于基于FISH語言二次開發而言[7]，其二次開發效果具有更好的穩定性及方便工程技術人員采用的特點。

1 分數階導數

1.1 分數階導數簡介

分數階導數理論相比于整數階導數理論起步較晚。且應用范圍并不廣泛，在人們長期的科學實踐應用過程中，研究人員發現有些問題采用整數階并不能很好的描述其發展過程，而分數階能夠給出合理的解釋且參數簡潔明確。直到上世紀80年代，分數階導數才得到廣泛應用，主要集中在天氣和氣候研究、醫學圖像處理、工程中非線性問題處理等，很好地解釋了自然科學以及工程領域一些非經典現象。由于分數階導數計算的復雜性及無法向整數階那樣得到精確的解析解，現今對分數階導數還沒有明確的定義，但應用廣泛的主要有三種，分別是Riemann-Liouville定義、Grunwald-Letnikov定義和Caputo定義[8-9]，三類定義之間既有聯系又有區別。

分數階相比于傳統整數階理論，優點主要在于：① 能夠很好的體現系統函數發展過程；② 與非線性模型相比，物理意義更明確，參數更簡潔；③ 克服了經典整數階理論與試驗數據不吻合的嚴重缺點，使用較少參數即可獲得很好的效果[10]。

1.2 分數階導數元件

分數階導數元件最主要的特點就是非線性[11-12]，相比傳統線性牛頓粘壺元件，更適合描述土體、橡膠、瀝青等具有粘性的復雜材料體，該元件最早由挪威數學家阿貝爾提出，又名Abel（阿貝爾）元件，如圖1。

圖1 阿貝爾元件Fig.1 Abel element

對上述元件，其模型本構關系為

1.3 分數階導數模型

本文建立的兩元件分數階導數模型，是將Abel元件分別代替經典Burger模型中Maxwell粘壺和Kelvin粘壺而得到（見圖2，圖3）。

圖2 Burger模型Fig.2 Burgermodel

圖3 Sburger模型Fig.3 Sburgermodel

對于Burger模型其由Maxwell模型和Kelvin模型串聯而成，分別包含一個虎克彈簧體和一個牛頓粘壺，可以用來描述土體的粘彈性，其一維本構方程如下

上式及后式中均統一約定，E表示彈性模量，η表示粘滯系數，G表示剪切模量，下標M和K分別代表Mawell部分參數和Kelvin部分參數

相應三維本構方程如

上式中σ1=σ1-σ3，σ11=σ1+2σ3，其它參數同上。

對于兩元件分數階導數模型本構方程的推導，分兩部分，第一部分導出分數階Maxwell部分本構關系；第二部分導出分數階Kelvin部分本構關系；上述兩部分串聯便得到了兩元件分數階導數Burger模型，所以這兩部分本構關系的和即為兩元件分數階導數Burger模型本構關系，具體推導如下：

分數階Maxwell部分由一個虎克彈簧體和一個阿貝爾(Abel)元件串聯而成，基于前述Abel元件本構關系，很容易給出該部分一維本構關系如下

式中：ηa為阿貝爾粘壺基礎粘滯系數，r為分數階階數。

分數階Kelvin部分由一個虎克彈簧體和一個Abel元件串聯而成，其本構關系推導如下[13]

式中：Q＝EK／ηK，D為微分算子，r為分數階階數，若r＝1即為整數階。

將Heaviside單位階躍函數θ（t）代入式（5）并令θ（t）＝σ（t），表示在t＞0時在模型兩端施加恒定的單位應力，通過Laplace變換再反演便得到了式（6）的蠕變柔量表達式

上式中Y1＝1／EK，由展開級數再逐項反演可得，因此蠕變柔量可進一步寫為

由于上式含有無窮級數和Gamma函數，形式復雜，不便于直接應用，需對上式進一步簡化。在上式無窮級數內添加n=-1項，同時在無窮級數外減去此項，并進一步令r=1-β，0≤β＜1，實際表明軟粘土蠕變過程與時間密切相關，為了保留蠕變時間的指數關系，需對無窮級數內的Gamma函數進行簡化，令β=0，則上式變為

式（9）即為分數階Kelvin模型的蠕變柔量，若β＝0即成為線性粘彈性模型，由式（9）和式（4）相加便得到兩元件分數階Burger模型一維本構方程，如下式

相應三維本構方程形式如下

若對式（11）作個簡單變換，令ηK＝ηKtβ，ηM=ηaΓ（1＋r）t1-r，便可進一步得到簡化后的下式

觀察式（12）與Burger三維本構關系式式（3）進行對比發現，經過替換后的Sburger模型本構關系同Burger模型三維本構關系完全一樣，兩者的區別之一在于分數階導數模型中粘滯系數ηK＝ηKtβ和ηM=ηaΓ（1＋r）t1-r是隨時間成指數形式變化，而Burger模型中這兩類粘滯系數始終為常數，基于這一特點，采用C++語言將分數階導數模型在FLAC3D中實現二次開發成為可能。

2 分數階導數模型二次開發

2.1 二次開發要點

得出兩類模型差異后，將Sburger模型在FLAC3D中進行二次開發，主要工作是通過C++語言建立循環結構[14]，實現參數ηK和ηM隨蠕變時間呈指數關系變化。借鑒FLAC3D中內置的Burger模型源程序(包含. cpp源文件和.h頭文件)編寫規則和各函數功能，以該源程序為藍本通過if..else條件選擇結構和for循環語句實現兩類粘滯系數隨時間的變化關系。并將更改后的分數階模型程序文件以項目形式添加入visual studio 2008中生成Dll(動態鏈接庫)文件[15]，并命名為userSburger.dll。然后將該鏈接復制到FLAC3D內置模型程序文件目錄下，通過在FLAC3D窗口輸入相關命令，完成所開發模型同FLAC3D的接口、載入，初步實現該模型二次開發。

模型開發過程中對程序的更改主要集中在基本參數的定義、常量與過程量的定義、基類的描述、模型注冊ID號和Run（）函數中本構模型計算關系的修改等。

2.2 二次開發驗證

在開發模型同FLAC3D進行接口、載入后，需要對模型開發的正確性及用于數值計算的可靠性進行驗證。驗證方法分為兩部分，一是判斷所開發模型程序編寫是否正確，即該模型能否用于數值計算；其次判斷所開發模型是否實現兩類粘滯系數隨時間成指數形式變化，這是開發成功的關鍵。首先通過建立一個簡單算例，在Sburger模型中對兩個分數階階數γ和β取值0退化成Burger模型形式，得ηK=ηK和ηM=ηa；保持其它參數不變，計算結果如圖4。

圖4 Sburger模型退化后計算結果Fig.4 Calculation of deteriorated Sburgermodel

從圖4可以看出，對Sburger模型的分數階階數取值為0后，保持其他參數相同，其數值計算結果與內置Burger模型完全吻合。說明兩類分數階模型按照上述方式通過取特殊值可退化為經典Burger模型，也進一步表明開發的兩類分數階模型是可以在FLAC3D中作為本構模型來選取并參與數值計算。

圖5 Sburger模型中不同階數數值計算曲線Fig.5 Numerical calculation curves of different orders in Sburgermodel

從圖5可得：對Sburger模型階數取不同值，數值計算曲線存在明顯差異，其基本規律為當β取小，γ取大時，試樣在等速蠕變階段的變形速率較大；而β取大，γ取小時，試樣在等速蠕變階段的變形速率較小；出現這一規律的原因在于β取值越小使得ηK隨時間增加的變化率下降，即單位時間內的粘滯系數增長率降低，導致在恒定偏應力作用下，其對應蠕變變形量增加；反之β取值越大，單位時間內粘滯系數增長率較高，使得對應蠕變變形量降低。同理對于另一階數γ，其取相同值所得ηM的結果正好與β相反。綜上，兩個分數階階數的不同取值，將影響兩類粘滯系數隨時間變化增長速率的快慢，從而在恒定偏應力下，使得等速蠕變階段蠕變變形存在顯著差異。

3 三軸數值算例

3.1 三軸蠕變試驗數據

兩類模型的基本參數均依據三軸蠕變試驗數據通過各模型三維本構關系擬合得到，表1為湛江軟粘土試樣基本物理力學性質[16]，圖6為該試樣在50 kPa恒定圍壓不同偏應力下，柱體軸向應變量隨時間變化關系。

表1 湛江軟粘土試樣基本物理力學性質Tab.1 Basic physical andmechanical properties of Zhanjiang clay samp le

圖6 三軸蠕變試驗曲線Fig.6 Curves of triaxial creep test

圖6是對試樣進行分級加載條件下得到，而只有分別加載才能得到更加真實的土體蠕變曲線，且考慮到巖土材料流變的非線性特性。需將分級加載下的蠕變數據轉化為分別加載下的蠕變數據，這里采用陳宗基先生提出并由其學生發展的“陳氏法”[17]進行處理，整理后的蠕變數據如圖7。

圖7 經“陳氏法”整理后三軸蠕變試驗曲線Fig.7 Curves of triaxial creep test under Chen’smethod

3.2 模型建立

為了對比分數階導數模型與Burger模型描述湛江軟粘土蠕變過程的可靠性，在FLAC3D中建立一個圓柱體模型(如圖8)算例模擬三軸蠕變過程。圓柱高為2個單位，長和寬均為0．5個單位，半徑為0．25個單位，先對模型施加50 kPa的圍壓和重力加速度，然后模型各方向位移清零，再在上頂面施加豎直方向壓力，蠕變時步取為1e-4，分別采用Burger模型和本文開發的分數階導數Sburger模型進行計算。模型參數取值依據相應本構關系基于三軸蠕變試驗數據擬合得到。

圖8 柱體模型示意圖Fig.8 Schematic diagram of cylindermodel

3.3 參數選取

采用1Stopt數學軟件及孫均[18]提出的擬合方法擬合得到兩類模型的基本參數，大致可分為兩部分，第一部分當t→o時，只有Maxwell彈簧存在變形，可擬合出參數K和GM，第二部分當t→∞時，蠕變達到穩定，將此時本構關系式減去原本構關系式便得到包含粘性和粘彈性本構關系式，與對應的三軸試驗數據進行回歸擬合。得到Burger模型在不同偏應力下參數如表2，同理可得分數階Sburger模型參數如表3。并由擬合相關系數R可以看出，各參數擬合情況較好。

表2 恒定圍壓不同偏應力下burger模型相關參數Tab.2 Parameters of Burgermodel under different partial stress of constant confining pressure

表3 恒定圍壓不同偏應力下Sburger模型相關參數Tab.3 Parameters of Sburgermodel under different partial stress of constant confining pressure

由表2和表3中兩個模型的相關擬合系數R可知，Sburger模型的擬合效果要優于Burger模型。其原因在于湛江黏土的蠕變包含有非線性變形部分，而傳統線性模型無法準確描述非線性問題，所以采用非線性本構關系擬合參數的效果要比線性本構關系好；其次，Burger模型在30 kPa偏應力下，擬合系數較高，這是在低偏應力下，黏土蠕變主要產生線性彈性、粘性和粘彈性變形；而當偏應力增加，隨著蠕變變形的進一步增長，非線性變形也在不斷上升成為后期變形的主要增長部分，是導致65～135 kPa偏應力下Burger模型擬合效果較低的原因之一。另外，從參數擬合效果來看，在數學角度上能一定程度說明Sburger模型本構關系要比經典Burger模型更能描述試驗曲線。

將上述兩模型采用對應參數分別在不同偏應力下進行三軸數值計算，并將兩者數值計算結果同三軸蠕變試驗數據進行對比，如圖9～圖12。

由圖9～圖12不難看出，開發兩元件分數階導數模型比傳統經典Burger模型，更適合描述湛江軟粘土的初始蠕變和穩定蠕變階段。從而驗證了本文提出采用分數階元件替換經典整數階模型中牛頓粘壺元件，建立變參數流變模型，描述軟土非線性蠕變思路的正確性；也進一步說明，本文對分數階導數模型在

FLAC3D中實現二次開發的正確性和可行性。

圖9 偏應力30 kPa下數值計算曲線Fig.9 Numerical calculation curve under the partial stress 30 kPa

圖10 偏應力65 kPa下數值計算曲線Fig.10 Numerical calculation curve under the partial stress 65 kPa

圖11 偏應力100 kPa下數值計算曲線Fig.11 Numerical calculation curve under the partial stress 100 kPa

圖12 偏應力135 kPa下數值計算曲線Fig.12 Numerical calculation curve under the partial stress 135 kPa

4 結論

從本文提出分數階元件替換傳統線性粘壺元件，建立變參數蠕變模型，并在FLAC3D中完成二次開發和驗證的結果來看，分數階元件相對于傳統元件，對蠕變曲線的刻畫能力更強，預示著對工程的數值模擬結果更加準確。由于土體蠕變變形具有廣泛的非線性特性，分數階元件相對于整數階元件，將更具備普遍意義。伴隨著工程建設的需要和蠕變模型自身的發展，本文研究成果將在未來的巖土工程設計計算中被越來越多地接納和采用。

[1]楊愛武，張兆杰，孔令偉.不同應力狀態下軟黏土蠕變特性試驗研究[J].巖土力學，2014，35（S2）：53-60.

[2]雷華陽，賈亞芳，李肖.濱海軟土非線性蠕變特性的試驗研究[J].巖石力學與工程學報，2013，31（1）：2806-2816.

[3]拓勇飛，孔令偉，郭愛國，等.湛江地區結構性軟土的賦存規律及其工程特性[J].巖土力學，2004（12）：1879-1884.

[4]蔡羽，孔令偉，郭愛國，等.剪應變率對湛江強結構性黏土力學性狀的影響[J].巖土力學，2006（8）：1235-1240.

[5]何利軍，蔡羽，吳應紅.湛江軟黏土的流變本構模型研究[J].合肥工業大學學報：自然科學版，2011（4）：565-570.

[6]何利軍，孔令偉，吳文軍，張先偉，蔡羽.采用分數階導數描述軟黏土蠕變的模型[J].巖土力學，2011（S2）：239-243.

[7]何利軍.分數階導數蠕變模型及在FLAC3D中的數值實現[J].南昌航空大學學報：自然科學版，2015（4）：35-39.

[8]李秀.幾類粘彈性材料的分數階積分模型及其數值解法[D].寧夏：寧夏大學，2013：8-24.

[9]王玉嬌.分數階導數及其應用[D].太原：太原理工大學，2014：10-44.

[10]陳文，孫洪廣，李西成.力學與工程問題的分數階導數建模[M].北京：科學出版社，2010：124-140.

[11]郭佳奇，喬春生，徐沖，等.基于分數階微積分的Kelvin-Voigt流變模型[J].中國鐵道科學，2009，30（4）：1-6.

[12]何志磊，朱珍德，朱明禮，等.基于分數階導數的非定常蠕變本構模型研究[J].巖土力學，2016，37（3）：737-744.

[13]王志方，張國忠，劉剛.膠凝原油的黏彈性流變模型[J].高校化學工程學報，2008，22（2）：351-355.

[14]古萬榮著.Visual C++完全自學手冊[M].北京：機械工業出版社，2009：29-111.

[15]陳育民，徐鼎平.FLAC/FLAC3D基礎與工程實例[M].2版.北京：中國水利水電出版社，2013：50-165.

[16]孔令偉，何利軍，張先偉.湛江黏土的蠕變模型與變參數塑性元件[J].巖土力學，2012，33（8）：2241-2246.

[17]陳宗基，康文法.巖石的封閉應力、蠕變和擴容及本構方程[J].巖石力學與工程學報，1991，10（4）：299-212.

[18]孫鈞著.巖土材料流變及其工程應用[M].北京：中國建筑工業出版社，1999：146-158.

Secondary Development of Non-linear Creep Modelw ith Fractional Order Derivative Elements

Zhang Tao，He Lijun
（School of Civil Engineering and Architecture，Nanchang Hangkong University，Nanchang 330063，China）

Triaxial creep test of Zhanjiang clay was conducted at constant confining pressure and under different deviatoric stress through GDS,which obtained several groups of data concerning the relationship between axial strain and time under different deviatoric stress.In view of the fact that traditional linear elements of nonlinear creepmodel does not accurately describe the clay problem,a fractional order derivative elementwas proposed to replace the two linear elements of the Newton’s pot in the Burgermodel respectively.Then a non-linear creepmodel with two fractional derivatives was established.By using C++language,the secondary development of themodel in the FLAC3D was completed.Finally,comparison wasmade between numerical results under different stress and experimental data by establishing a triaxial numerical example.The results show that compared with the Burgermodel,the non-linear creepmodel ismore favorable for the description of creep process for Zhanjiang clay.

fractional derivative；secondary development；creep

1005-0523（2017）05-0021-08

T U433

A

2017－05－09

南昌航空大學校級基金項目（EA201111317）

張濤（1991一），男，碩士研究生，主要研究方向為蠕變模型的二次開發。

何利軍（1977一），男，講師，博士，研究方向為結構性土體本構模型研究。

（責任編輯 王建華）