黏彈－Perzyna黏塑性有限元法應(yīng)力更新隱式算法

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

（鄭州大學(xué) 土木工程學(xué)院，鄭州450001）

1 引 言

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

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

本文針對黏彈性和雙曲型Drucker－Prager屈服函數(shù)［15］、各向同性硬化及 Perzyna黏塑性［2，3］耦合模型，給出完全隱式的應(yīng)力更新算法和最終公式。首先，基于黏彈性蠕變?nèi)岫?，通過定義與彈性問題相對應(yīng)的與時間增量相關(guān)的黏彈性剪切模量和體積模量，實現(xiàn)遺傳積分方程的線性化，然后與黏塑性耦合。為了保證計算過程的穩(wěn)定性，把Perzyna黏塑性定律轉(zhuǎn)化為與彈塑性問題相對應(yīng)的一致性條件，建立了黏塑性增量因子在大于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ù)。把應(yīng)變率張量ε·分解為黏彈應(yīng)變率和黏塑性應(yīng)變率［3，4］

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

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

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

偏應(yīng)變和體積應(yīng)變分別由遺傳積分方程給出

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

2．2 黏塑性

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

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

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

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

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

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

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

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

黏塑性應(yīng)變率由流動法則給出［3］

式中g(shù)，σ為勢函數(shù)關(guān)于應(yīng)力σ的偏導(dǎo)數(shù)。

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

式中 模型參數(shù)m為黏塑性敏感性指數(shù)。當(dāng)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，偏應(yīng)力增量記為 Δsn＋1＝sn＋1－sn，其他參量類同。在Δtn＋1時步內(nèi)，假設(shè)應(yīng)力線性變化，即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）為考慮黏彈性應(yīng)變歷史和在每個時間步長內(nèi)逐段線性化后用增量形式表示的黏彈性本構(gòu)關(guān)系。對于平均應(yīng)力p和黏彈性體積應(yīng)變εvve，可根據(jù)相應(yīng)的黏彈性模型，采用上述類似方法，得到n＋1時刻平均應(yīng)力的表達式和與彈性問題相對應(yīng)的黏彈性體積模量為

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

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

3．2 黏塑性

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

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

對于等效黏塑性應(yīng)變，由式（9）得

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

3．3 隱式應(yīng)力更新算法

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

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

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

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

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

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

式中

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

3．4 一致切向矩陣

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

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

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

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

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

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

4 算例分析和討論

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

4．1 算例分析

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

黏彈性部分的應(yīng)力偏量和平均應(yīng)力均采用一個彈簧、一個黏壺和兩個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，相對應(yīng)的黏彈性參數(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，選用關(guān)聯(lián)流動法則，Ψ＝φ。有效重度取為9．0kN／m3。黏塑性黏性系數(shù)ηvp＝150MPa·d。

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

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

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

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

豎向平面內(nèi)等效塑性應(yīng)變云圖如圖4所示。由于等效黏塑性應(yīng)變是根據(jù)應(yīng)力和黏塑性應(yīng)變定義的，在荷載對稱面內(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所示?？梢钥闯?，在荷載施加到最大值及其附近時，兩者差別最大。隨著蠕變的增加，誤差很小，并沒有因為時間步長增大一倍而產(chǎn)生計算結(jié)果的漂移。本文采用了完全隱式向后Euler算法，是可以采用較大時間步長而不出現(xiàn)應(yīng)力漂移的一個重要原因。

圖4 土中累計等效黏塑性應(yīng)變云圖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就無法計算。在初始屈服或應(yīng)力點距離屈服面很近時，由于Δλ的值往往很?。ㄈ?0－20量級），Δλ的初始值不能取為0，且需要在Δλ＞0的一側(cè)進行逼近，否則迭代不收斂。該問題是黏塑性隱式算法是否收斂的關(guān)鍵，也是不同于一般彈－塑性隱式算法之處［2，3，16］。文獻［16］認為，Δλ取較好的初值才能保證收斂性。由于應(yīng)力點偏離屈服面的程度不同，該初值往往很難選取。計算表明，如果Δλ的初值取為顯式算法的計算值，應(yīng)力點偏離屈服面較大時迭代基本都能收斂，但應(yīng)力點在屈服面附近時往往不收斂。

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

圖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)岫?，考慮黏彈性應(yīng)變歷史，假設(shè)應(yīng)力在一個時間步長內(nèi)線性變化，定義了與彈性問題相對應(yīng)的與時間增量相關(guān)的黏彈性剪切模量和體積模量，導(dǎo)出了增量遞推形式的本構(gòu)方程，可方便地與黏塑性模型耦合計算。增加或減少串聯(lián)的黏彈性數(shù)量，計算程序無需修改，實現(xiàn)了計算程序的通用性。

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

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