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

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

1 引 言

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

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

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

本文中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。∶表示兩個(gè)張量的雙點(diǎn)積，表示張量積，張量s的模定義為s＝

2 本構(gòu)方程

2．1 黏彈性

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

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

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

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

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

式中s和p分別為偏應(yīng)力張量和平均應(yīng)力，Jg（0）和Jk（0）分別為t＝0時(shí)所對(duì)應(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)摩擦角計(jì)算得到。q＝，l0可由d的初值d0等參數(shù)計(jì)算，d為與硬化相關(guān)的參量。

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

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

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

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

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

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

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

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

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

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

3．1 黏彈性

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

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

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

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

式（16）為考慮黏彈性應(yīng)變歷史和在每個(gè)時(shí)間步長(zhǎng)內(nèi)逐段線性化后用增量形式表示的黏彈性本構(gòu)關(guān)系。對(duì)于平均應(yīng)力p和黏彈性體積應(yīng)變?chǔ)舦ve，可根據(jù)相應(yīng)的黏彈性模型，采用上述類似方法，得到n＋1時(shí)刻平均應(yīng)力的表達(dá)式和與彈性問(wèn)題相對(duì)應(yīng)的黏彈性體積模量為

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

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

3．2 黏塑性

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

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

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

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

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

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

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

把試應(yīng)變和實(shí)際應(yīng)變?cè)隽看胧剑?5，17）后，并考慮式（10，23，24），得到實(shí)際應(yīng)力和試應(yīng)力的關(guān)系分別為

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

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

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

式中

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

3．4 一致切向矩陣

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

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

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

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

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

需要說(shuō)明的是，在一致切向矩陣的公式中，采用了四階張量Is，而不是四階單位張量I，是因?yàn)橛邢拊浖胁捎玫氖枪こ虘?yīng)變。兩者在張量運(yùn)算過(guò)程中無(wú)差別，其差別僅體現(xiàn)在最后的計(jì)算上［3，4］。

4 算例分析和討論

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

4．1 算例分析

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

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

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

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

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

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

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

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

4．2 討 論

時(shí)間步長(zhǎng)的大小是影響計(jì)算精度的一個(gè)重要因素［3，4，15，16］。上述算例的時(shí)間跨度為910d，時(shí)間步長(zhǎng)Δt≈4d，加載過(guò)程中每時(shí)步的荷載增量約為11．11kPa。表1給出了普通個(gè)人電腦上本文算法、文獻(xiàn)［4］的黏彈－塑性算法和 ABAQUS黏彈性［14］計(jì)算的CPU時(shí)間比較?？梢钥闯觯疚牡酿棧に苄杂?jì)算較ABQUS黏彈性計(jì)算約慢14%，盡管采用了隱式時(shí)間差分算法，其計(jì)算效率仍較高。

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

圖4 土中累計(jì)等效黏塑性應(yīng)變?cè)茍DFig．4 Contours of equivalent viscoplastic strain in soil

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

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

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

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

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

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