閆富有， 崔 昊， 張曉婉， 劉忠玉

（鄭州大學 土木工程學院，鄭州450001）

1 引 言

一些巖土或混凝土等材料，具有顯著的流變性，可采用黏彈－黏塑性模型來描述其應力應變關系［1］。黏彈性部分通常由基本黏彈性元件串聯(lián)或并聯(lián)的不同組合來描述，黏塑性部分可采用簡單的Bingham模型或更具一般性的Perzyna等模型［1，2］來描述與應變率相關的塑性流動。在計算過程中，允許應力點偏離屈服面，在應力更新過程中無需讓應力點返回屈服面［2，3］。對于應力點必須返回屈服面，即滿足屈服函數(shù)等于0的條件，稱之為黏彈－塑性模型［4］，以便區(qū)別于本文的黏彈－黏塑性模型。無論是否屈服、加載或卸載，其變形特性均表現(xiàn)為與時間相關的高度非線性，且與應變歷史相關。數(shù)值計算模型通常比較龐大，計算時間長，其計算效率、精度和收斂性以及解的穩(wěn)定性在很大程度上取決于應力更新算法［2－4］。

顯式積分算法是目前普遍應用的黏彈－黏塑性本構(gòu)積分算法［5－8］。計算簡單，但未考慮黏彈性應變歷史，需要嚴格的更小的時間步長限制，應力漂移現(xiàn)象值得商榷［3，4］。為了彌補顯示算法的不足，文獻［4，9，10］基于黏彈性松弛模量Prony級數(shù)的表達式導出應力遞推公式，然后與塑性［4］或黏塑性耦合［9，10］，后者僅限于簡單的 Mises屈服函數(shù)和黏塑性流動模型。對于常用的由Kelvin鏈與彈簧或黏壺串聯(lián)而成的黏彈性模型，其松弛模量并非顯式，復雜模型松弛模量的推演也較難，導致計算程序很難通用化。基于蠕變?nèi)岫葘С鰬Γ瓚冞f推公式的方法源于文獻［11］，后在文獻［12，13］得到應用。該方法避免了松弛模量方法的不足，但把Mises屈服應力導入遞推公式［11］，并假定與應力相關的一些變量在一個時步內(nèi)為常數(shù)，不便應用于其他黏塑性問題。有限元法應力更新算法包括應力的更新和一致切向矩陣的計算，一些文獻僅給出初步的應力遞推算法［12，13］，而無具體計算公式，更沒有表述與算法對應的一致切向矩陣［5，12，13］。研究多采用 Bingham 黏塑性模型［5－13］，而更具一般性的Perzyna模型［3］很少應用于巖土工程中。一些商業(yè)有限元軟件只能進行黏彈性或彈－黏塑性分析［14］，尚不能進行黏彈－黏塑性耦合計算，限制了其應用。

本文針對黏彈性和雙曲型Drucker－Prager屈服函數(shù)［15］、各向同性硬化及 Perzyna黏塑性［2，3］耦合模型，給出完全隱式的應力更新算法和最終公式。首先，基于黏彈性蠕變?nèi)岫?，通過定義與彈性問題相對應的與時間增量相關的黏彈性剪切模量和體積模量，實現(xiàn)遺傳積分方程的線性化，然后與黏塑性耦合。為了保證計算過程的穩(wěn)定性，把Perzyna黏塑性定律轉(zhuǎn)化為與彈塑性問題相對應的一致性條件，建立了黏塑性增量因子在大于0的條件下，由小于其收斂值的一側(cè)逐步逼近其收斂值的單側(cè)逼近算法。最后，通過ABAQUS軟件提供的UMAT材料子程序完成計算，并對其收斂性進行討論。

本文中1表示二階單位張量，（1）ij＝δij，I和Is均為四階張量，（I）ijkl＝δikδjl，（Is）ijkl＝（δikδjl＋δilδjk）／2，δij為 Kronecker delta，i，j，k，l＝1，2，3?！帽硎緝蓚€張量的雙點積，表示張量積，張量s的模定義為s＝

2 本構(gòu)方程

2．1 黏彈性

對于圖1所示的黏彈－黏塑性模型，其黏彈性部分是由一個彈簧、一個黏壺和多個Kelvin體串聯(lián)而成。若僅有一個Kelvin體相串聯(lián)，則為常規(guī)的Burgers體；若無黏壺η0，則為標準線性黏彈性體。ηvp為黏塑性黏性系數(shù)。把應變率張量ε·分解為黏彈應變率和黏塑性應變率［3，4］

式中 上標ve和vp分別代表對應參量的黏彈和黏塑性部分，變量上部的·為相對于時間的變化率。

黏彈性偏應變張量eve和體積應變εvev所對應的黏彈性蠕變?nèi)岫菾g和Jk分別為［11，12］

式中t為時間，下標g和k分別為與剪應力和平均應力所對應的量。λgi＝ηgi／Gi，λki＝ηki／Ki。Gi和Ki，ηgi和ηki，ηg0和ηk0分別為與圖1的Gi，ηi和η0相對應的模型參數(shù)。NG和NK分別為剪應力和平均應力所對應的黏彈性模型中串聯(lián)的Kelvin體個數(shù)，i＝1，2，…，NG／NK。

偏應變和體積應變分別由遺傳積分方程給出

式中s和p分別為偏應力張量和平均應力，Jg（0）和Jk（0）分別為t＝0時所對應的量。

2．2 黏塑性

雙曲線Drucker－Prager塑性屈服函數(shù)為［14］

圖1 黏彈－黏塑性模型Fig．1 A viscoelastic－viscoplastic model

式中β為與內(nèi)摩擦角相關的參量，可由 Mohr－Coulomb與Drucker－Prager模型參數(shù)的轉(zhuǎn)換關系用內(nèi)摩擦角計算得到。q＝，l0可由d的初值d0等參數(shù)計算，d為與硬化相關的參量。

雙曲型Drucker－Prager塑性模型雖然增加了計算過程的復雜性，但保證了屈服函數(shù)的連續(xù)性，使得在屈服面錐頂塑性流動方向具有唯一性，消除了線性Drucker－Prager屈服面錐頂?shù)膽Σ贿B續(xù)性所導致的在其附近應力收斂的困難［15］，提高了計算效率和精度。兩者的接近程度與參數(shù)l0的大小有關，當l0＜0．01時，運用線性和雙曲線屈服函數(shù)所得的結(jié)果相當接近［15］。

把屈服函數(shù)和勢函數(shù)寫成一般形式，分別為

式中η＝tanβ，η－為將內(nèi)摩擦角換為剪脹角Ψ后計算得到的η值。l0和m0均取一較小的正數(shù)，以便更好地接近于線性屈服函數(shù)。

僅考慮各向同性硬化，d為等效黏塑性應變ε－vp的函數(shù)，dd／dε－vp＝H，H 為硬化模量。對于線性硬化，H為常數(shù)，d＝d0＋Hε－vp。等效黏塑性應變率ε－·vp定義為［20］

式中 Δεvp為Δt時間段的黏塑性應變增量，λ·為黏塑性增量因子，σ－為等效應力，本文取為單軸抗壓強度sc。

黏塑性應變率由流動法則給出［3］

式中g(shù)，σ為勢函數(shù)關于應力σ的偏導數(shù)。

采用Perzyna黏塑性模型［2，3］，黏塑性增量因子定義為

式中 模型參數(shù)m為黏塑性敏感性指數(shù)。當m＝1時，即為Bingham黏塑性模型。

3 本構(gòu)積分算法

3．1 黏彈性

把式（2）代入式（4），整理后得

把時間t離散，下標n（n＝0，1，2，…）表示不同時步。t0＝0。第n＋1時步的時間增量記為Δtn＋1＝tn＋1－tn，偏應力增量記為 Δsn＋1＝sn＋1－sn，其他參量類同。在Δtn＋1時步內(nèi)，假設應力線性變化，即ds（τ）／dτ為常數(shù)。對式（13）進行分段積分后，得tn＋1時刻ζgi（tn＋1）的遞推公式為

式中 h（x）＝［exp（x）－1］／x（變量x＝λgiΔtn＋1）。

把式（14）代入式（12），分別取t＝tn＋1和t＝tn，然后兩式相減，并考慮最后一項的積分，得

式（16）為考慮黏彈性應變歷史和在每個時間步長內(nèi)逐段線性化后用增量形式表示的黏彈性本構(gòu)關系。對于平均應力p和黏彈性體積應變εvve，可根據(jù)相應的黏彈性模型，采用上述類似方法，得到n＋1時刻平均應力的表達式和與彈性問題相對應的黏彈性體積模量為

式中 狀態(tài)變量ζki的表達式與式（13，14）相似，僅把其中的s和λgi分別換為p和λki即可。

由式（15，17）可知，下一時刻的應力增量由本時步黏彈性應變增量的影響項和考慮應變歷史后狀態(tài)變量ζgi及ζki的影響項兩部分組成。

3．2 黏塑性

在計算過程中，黏塑性增量因子需滿足式（11），即當屈服函數(shù)f＞0時，產(chǎn)生黏塑性應變。對時間差分后得

一般情況下，m＞1。如果f／d較小或較大，在計算過程中，式（19）往往不穩(wěn)定［3，16］。把式（19）改寫為與一般塑性問題類似的一致性條件為

對于等效黏塑性應變，由式（9）得

屈服函數(shù)和塑性勢函數(shù)對σ的偏導數(shù)分別為

3．3 隱式應力更新算法

在3．1節(jié)中，定義了與時間增量相關的與彈性問題相對應的黏彈性剪切模量和體積模量，其增量形式的本構(gòu)關系在形式上與彈性問題類似，可參考彈－黏塑性本構(gòu)積分算法［3］進行應力更新。

在第n步計算結(jié)束后，對于新的時間段 ［tn，tn＋1］，給定總應變增量Δεn＋1，采用預估－校正兩步法進行計算［3，4，11～13］。首先，把應變增量 Δεn＋1作為黏彈性試應變增量 Δεnve＋，t1ri，即Δεnve＋，1tri＝Δεn＋1。求得試應力增量Δσntri＋1及試應力σntri＋1。用試應力求得屈服函數(shù)f，若f≤0，該試應力即為實際應力，本時步計算結(jié)束；否則，進行黏塑性校正。

把試應變和實際應變增量代入式（15，17）后，并考慮式（10，23，24），得到實際應力和試應力的關系分別為

式中qntri＋1為由試應力求得的qn＋1，qm省略了表示時間的下標n＋1。

把式（25，26）一起代入黏塑性一致性條件式（20），并考慮c2為未知量，建立如下方程組，

把式（28，29）改寫為 N－R迭代格式

式中

求解方程組（30，31）后得到δ（Δλ），更新Δλ，直到r1和r2小于規(guī)定的誤差限（取10－8）后結(jié)束迭代。用收斂后的Δλ，應用式（25，26）即可求得平均應力和偏應力，黏塑性校正過程結(jié)束。然后，計算狀態(tài)變量ζgi，ζki和ε－vp的值，以便下一時刻黏彈性計算。

3．4 一致切向矩陣

對于黏彈性狀態(tài)，其一致切向矩陣Dve可仿照彈性問題給出［3，4］

對于黏塑性狀態(tài)，由一致性條件得

根據(jù)式（27）求得

把式（35）代入式（34），經(jīng)過冗長計算后得

由式（25，26）對平均應力和偏應力微分，并考慮式（35，36）后，得黏彈－黏塑性一致切向矩陣Dvevp的表達式為

需要說明的是，在一致切向矩陣的公式中，采用了四階張量Is，而不是四階單位張量I，是因為有限元軟件中采用的是工程應變。兩者在張量運算過程中無差別，其差別僅體現(xiàn)在最后的計算上［3，4］。

4 算例分析和討論

為了說明上述算法的有效性，采用ABAQUS提供的UMAT子程序［14］編制黏彈－黏塑性應力更新用戶子程序，然后對有限元法計算結(jié)果進行分析和討論。

4．1 算例分析

一地基模型，沿深度z、寬度x和另一方向y分別取10m，20m和1m。沿y方向劃分為等厚度的4個單元，并對該方向的兩側(cè)施加零位移約束，以模擬平面應變問題［3，4］。在地基表面中間部位沿x方向作用一寬度為1m且隨時間線性增加的條形荷載，荷載在180d內(nèi)由0增大到500kPa，此后保持恒定?？紤]對稱性，在x方向取其一半進行計算。模型四周邊界施加水平位移約束，模型底面約束三個方向的位移。分析采用三維有限元法，共劃分2075個結(jié)點和1520個8結(jié)點六面體單元。

黏彈性部分的應力偏量和平均應力均采用一個彈簧、一個黏壺和兩個Kelvin體串聯(lián)而成，為擴展的Burgers模型（較常規(guī)的Burgers模型多為一個串聯(lián)的Kelvin體），NG＝NK＝2。因試驗參數(shù)有限，黏彈性參數(shù)由一維流變試驗結(jié)果擬合得到，分別為 E0＝7．3MPa，η0＝15．6MPa·d，E1＝671．3MPa，η1＝25．3MPa·d，E2＝1246．8MPa，η2＝23570．8MPa·d。取泊松比為0．32，相對應的黏彈性參數(shù)分別為［4，9］

初始黏聚力c＝35kPa，內(nèi)摩擦角φ＝15°，單軸抗壓強度按σc＝2ccosφ／（1．0－sinφ）計算。采用與Mohr－Coulomb屈服面外切園近似的原則確定Drucker－Prager屈服面參數(shù)，即tanβ＝6sinφ／（3－sin φ），d的初值，l0和 m0均為0．001。按線性硬化考慮，硬化模量H＝37，選用關聯(lián)流動法則，Ψ＝φ。有效重度取為9．0kN／m3。黏塑性黏性系數(shù)ηvp＝150MPa·d。

地應力平衡后，施加面力荷載。計算的時間步長Δt≈4d。為了驗證本算法的有效性，取相同的黏彈性參數(shù)，將黏塑性敏感系數(shù)取不同值時，本文的計算結(jié)果與ABAQUS中的黏彈性計算結(jié)果［14］及文獻［4］的黏彈－塑性計算結(jié)果進行比較。

分別取敏感性指數(shù)m＝1，5，10和20進行計算。地基表面荷載中心點處豎向位移和豎向應變εz隨時間變化的結(jié)果比較分別如圖2和圖3所示?？梢钥闯?，在加載過程初期，由于尚未屈服，不同模型的計算結(jié)果基本相同；土體屈服后，不同模型的計算結(jié)果差別很大。黏彈性模型土體豎向位移和應變最小，黏彈－塑性模型的結(jié)果最大，黏彈－黏塑性模型的土體位移介于黏彈性和黏彈－塑性模型的計算結(jié)果之間。不同的黏塑性敏感性指數(shù)m對位移影響較大。當m＝1時（Bingham模型），其黏彈－黏塑性模型的位移最大；隨著m的增大，變形相對減小。這與文獻［3，16］的結(jié)論一致。當m＝20時，黏彈－黏塑性模型的豎向位移和應變幾乎接近于黏彈性模型的結(jié)果。加載結(jié)束后，在恒定的荷載作用下，土體的蠕變行為持續(xù)增加，隨著作用時間的延續(xù)，采用Bingham黏塑性模型計算的位移將逐步趨于黏彈－塑性的結(jié)果。因此，在實際工程中，可采用不同的敏感性指數(shù)試算，通過與實際結(jié)果比較擬合來確定一個合適的m值。

圖2 荷載中心處豎向位移隨時間變化結(jié)果比較Fig．2 Comparison of vertical displacement at the midpoint of loading

圖3 荷載中心處豎向應變隨時間變化結(jié)果比較Fig．3 Comparison of vertical strain component at the midpoint of loading

豎向平面內(nèi)等效塑性應變云圖如圖4所示。由于等效黏塑性應變是根據(jù)應力和黏塑性應變定義的，在荷載對稱面內(nèi)呈現(xiàn)較大的值，塑性區(qū)的深度是有限的，其形狀與彈－黏塑性［2］或黏彈－塑性［4］計算結(jié)果相似。

4．2 討 論

時間步長的大小是影響計算精度的一個重要因素［3，4，15，16］。上述算例的時間跨度為910d，時間步長Δt≈4d，加載過程中每時步的荷載增量約為11．11kPa。表1給出了普通個人電腦上本文算法、文獻［4］的黏彈－塑性算法和 ABAQUS黏彈性［14］計算的CPU時間比較?？梢钥闯?，本文的黏彈－黏塑性計算較ABQUS黏彈性計算約慢14%，盡管采用了隱式時間差分算法，其計算效率仍較高。

為了考察時間步長對計算結(jié)果的影響，取Δt≈2d，所用CPU時間幾乎是Δt≈4d的2倍，兩者的豎向位移之差絕對值的最大值小于1．8×10－3m，一些時間點的相對誤差曲線如圖5所示?？梢钥闯觯诤奢d施加到最大值及其附近時，兩者差別最大。隨著蠕變的增加，誤差很小，并沒有因為時間步長增大一倍而產(chǎn)生計算結(jié)果的漂移。本文采用了完全隱式向后Euler算法，是可以采用較大時間步長而不出現(xiàn)應力漂移的一個重要原因。

圖4 土中累計等效黏塑性應變云圖Fig．4 Contours of equivalent viscoplastic strain in soil

表1 有限元模擬CPU時間比較Tab．1 Comparison of CPU time for FEM simulations

分析迭代過程的收斂性。在每個迭代步，均需滿足Δλ＞0和fn＋1＞0兩個條件，否則，c0和c4就無法計算。在初始屈服或應力點距離屈服面很近時，由于Δλ的值往往很?。ㄈ?0－20量級），Δλ的初始值不能取為0，且需要在Δλ＞0的一側(cè)進行逼近，否則迭代不收斂。該問題是黏塑性隱式算法是否收斂的關鍵，也是不同于一般彈－塑性隱式算法之處［2，3，16］。文獻［16］認為，Δλ取較好的初值才能保證收斂性。由于應力點偏離屈服面的程度不同，該初值往往很難選取。計算表明，如果Δλ的初值取為顯式算法的計算值，應力點偏離屈服面較大時迭代基本都能收斂，但應力點在屈服面附近時往往不收斂。

解決該問題的關鍵是讓Δλ在小于其收斂值的一側(cè)逐漸逼近于其收斂值，且Δλ不能等于0。一方面，Δλ的初值取一個很小但不為0的值，如本文取Δλ＝10－30；另一方面，指定一個較該初值更小的正數(shù)，若迭代過程中Δλ小于該正數(shù)，參數(shù)c0和c4采用上一步的值而不再更新，避免出現(xiàn)雙側(cè)逼近收斂值的情況。反復計算表明，上述處理方法可以保證最終收斂。圖6分別給出了取兩個不同屈服函數(shù)值f時，黏塑性因子的迭代收斂情況。圖6（a）收斂前后的屈服函數(shù)值差別很小，圖6（b）收斂后f＝49．8787，均滿足一致性條件。可以看出，無論應力點靠近屈服面還是偏離屈服面較大，經(jīng)過十多次迭代運算，均能獲得很好的收斂效果。由于已導出其解析式，計算效率和精度均比較高。

圖5 時間步長Δt＝4d與Δt＝2d的相對誤差Fig．5 Relative error betweenΔt＝4dandΔt＝2d

圖6 黏塑性因子收斂曲線（m＝5）Fig．6 Convergence profiles on viscoplastic multiplier from the numerical test（m＝5）

5 結(jié) 論

（1）對于由彈簧、黏壺和Kelvin鏈串聯(lián)而成的黏彈性模型，基于黏彈性蠕變?nèi)岫?，考慮黏彈性應變歷史，假設應力在一個時間步長內(nèi)線性變化，定義了與彈性問題相對應的與時間增量相關的黏彈性剪切模量和體積模量，導出了增量遞推形式的本構(gòu)方程，可方便地與黏塑性模型耦合計算。增加或減少串聯(lián)的黏彈性數(shù)量，計算程序無需修改，實現(xiàn)了計算程序的通用性。

（2）針對黏彈性和考慮各向同性硬化的雙線性型Drucker－Prager屈服函數(shù)和Perzyna黏塑性耦合而成的黏彈－黏塑性模型，采用黏彈性預估和黏塑性校正兩步算法，考慮算法的收斂和穩(wěn)定性，把Perzyna黏塑性流動定律轉(zhuǎn)化為與一般塑性問題相對應的一致性條件，建立了黏塑性增量因子單側(cè)逼近于收斂值的 N－R迭代算法，給出黏彈－Perzyna黏塑性耦合模型的完全隱式算法及最終的計算公式，為黏塑性隱式計算提供了一種行之有效的方法。

（3）算例分析表明，對于相同的模型參數(shù)和加載條件，應用黏彈性模型計算的變形量最小，黏彈－塑性模型的變形量最大，黏彈－Perzyna黏塑性模型的變形介于兩者計算結(jié)果之間。隨著黏塑性敏感性指數(shù)的增大，黏彈－Perzyna黏塑性模型計算的變形量相對減小；當增大到一定值時，其結(jié)果接近于黏彈性模型的計算結(jié)果。