混凝土彈塑性損傷本構模型的動力擴展

齊 虎，李云貴，呂西林

（1.同濟大學 結(jié)構工程與防災研究所，上海200092；2.中國建筑股份有限公司技術中心，北京101320）

混凝土結(jié)構都可能承受沖擊、地震和爆炸等動態(tài)載荷，基于混凝土準靜態(tài)力學性能的分析結(jié)果與實際情況會有較大的偏差［1］，因此，考慮應變率相關的本構模型在研究中逐漸占據(jù)重要的地位.在動力荷載作用下混凝土材料表現(xiàn)出應變率效應，即材料的損傷和非線性程度隨應變率的增加而降低，抗拉強度和抗壓強度則隨應變率的增加而提高等現(xiàn)象［2］.混凝土材料的應變率效應在沖擊荷載（如碰撞、爆炸等）以及地震荷載作用下非常明顯.

研究人員很早就注意到上述現(xiàn)象，提出了多種動力作用下的混凝土本構模型理論，然而，由于這些理論模型本身以及數(shù)值實現(xiàn)算法的復雜性，工程實際中應用最為廣泛的仍然是一些經(jīng)驗模型，如根據(jù)經(jīng)驗提高材料的單軸抗拉強度（或單軸抗壓強度）等簡化處理方法［3－4］.考慮應變加載速率的混凝土動力本構模型通常是在率不相關本構模型的基礎上做如下擴展得到［5］：①將率不相關本構模型破壞面擴大［6］，在實際工程中應用較多，但沒有物理背景；②將彈塑性應力—應變關系粘滯化，使得塑性應力—應變關系率相關［7］，在動力加載過程中混凝土微粒塑性應變的發(fā)展存在滯后；③考慮損傷演化的率相關性［8］，考慮了動力加載情況下混凝土材料損傷的滯后.

文獻［9］通過將塑性應變和損傷變量均粘滯化，將其提出的率不相關彈塑性損傷本構模型進行動力推廣，建立了相應的率相關本構模型，但模型較為復雜，使用不方便.文獻［3］認為：高應變率對微裂紋發(fā)展，即損傷演化的遲滯作用導致了混凝土材料非線性程度降低及動力強度提高.文獻［10］通過對損傷變量進行粘滯化，提出了適合大型鋼筋混凝土結(jié)構分析的混凝土率相關彈塑性損傷本構模型，文獻［2］采用和文獻［10］類似的方法將提出的靜力彈塑性損傷本構模型粘滯化，并通過算例表明應變率效應對實際結(jié)構的動力反應有一定的影響.

損傷力學能夠從本質(zhì)上描述混凝土材料內(nèi)部微裂紋的開裂、發(fā)展所引起的材料宏觀非線性行為，其中損傷變量的發(fā)展演化就是為了模擬材料內(nèi)部微裂紋的開裂發(fā)展，從而在損傷力學中可以將損傷變量粘滯化，通過延緩損傷變量的演化來模擬高應變率對微裂紋發(fā)展的遲滯作用.筆者在以前的研究中，建議了一個混凝土率不相關彈塑性損傷本構模型［11］，該模型能較好地描述混凝土材料在靜力或偽靜力加載下的各種非線性行為.基于以上分析，本文將筆者建議的模型進行動力推廣，使其能夠描述動力加載下混凝土材料的非線性行為.

阻尼是結(jié)構振動過程中一種特有的能量耗散機制［12］.目 前 應 用 最 為 廣 泛 的 是Rayleigh 阻 尼，Rayleigh阻尼假定阻尼矩陣C 是質(zhì)量矩陣M 和剛度矩陣K 的線性組合，數(shù)值計算表明：單純采用質(zhì)量比例阻尼也無法減小高頻噪音的影響［12］.考慮到質(zhì)量比例阻尼本身并無物理依據(jù)，很多學者［9，13］都僅僅考慮剛度比例阻尼的影響.本文將剛度阻尼應力引入到建立的彈塑性損傷本構模型中，使其能在材料層次考慮阻尼的影響.利用該模型對Koyna混凝土重力壩進行了地震動作用下的非線性時程分析，分析結(jié)果與文獻［9］符合較好，表明了建立模型的有效性.

本文首先簡要介紹筆者建議的靜力彈塑性損傷本構模型，然后對其進行擴展使之能夠考慮應變率效應和剛度阻尼的影響，最后通過算例從結(jié)構位移反應和數(shù)值穩(wěn)定性兩個方面探討應變率效應和剛度阻尼對結(jié)構動力反應的影響.

1 率不相關彈塑性損傷本構模型

損傷力學能描述材料的各種非線性本構行為，如剛度退化應變軟化等，目前用結(jié)合塑性力學和損傷力學的彈塑性損傷本構模型來模擬混凝土材料得到了較為廣泛的認同.由于塑性力學本身的復雜性，目前用大多彈塑性損傷本構模型數(shù)值處理復雜計算，效率低穩(wěn)定性不好，且涉及參數(shù)較多，一般只給出了特定情況下的材料參數(shù)取值，因此較難在實際工程中應用.筆者在以前的研究中提出了一個率不相關彈塑性實用本構模型［11］，通過對彈性Helmholtz自由能進行修正，使得模型能較為準確地模擬混凝土材料在雙軸、三軸加載下的非線性本構行為，其彈塑性部分使用經(jīng)驗彈塑性模型［10］降低模型的數(shù)值復雜性，下面首先對筆者開發(fā)的靜力彈塑性損傷本構模型做簡要介紹.

1.1 損傷變量的計算及參數(shù)的確定

損傷能量釋放率采用如下公式計算：

式中：標記“∶”為二階縮并積；Y±表示正、負損傷能量釋放率；ψe表示修正后的彈性Helmholtz自由能；表示受拉、受壓損傷變量表示有效應力張量；的正、負分量為有效柔度張量表示有效應力球量；‖X‖＝X∶X；γ 為本文引入的折減系數(shù)；δ表示向量（1，1，1，0，0，0）.

式中，H（x）為Heaviside函數(shù).

如果式（2）計算得γ＜0，則取γ＝0.式（1）中：

式中：標記“?”為張量積，Pii＝Pi?Pi為二階對稱張量；I為四階一致性張量.

求出損傷能量釋放率后，采用式（5）可計算出受拉、受壓損傷變量d＋、d－：

加權損傷變量：

1.2 塑性應變計算

文獻［10］給出如下塑性應變的經(jīng)驗計算公式：

式中：εp為塑性應變，dεp表示對εp取微分；E表示材料彈性模量；E0表示材料初始彈性剛度張量為單位有效應力張量；參數(shù)βp 控制塑性應變的大小，βp∈（0，1），0 表示線彈性，1 表示理想彈塑性.本文給出如下公式計算βp 取值：

以上經(jīng)驗計算公式（7）—（8）大大簡化了塑性應變的計算，而經(jīng)典塑性理論的流動法則、塑性硬化法則等都沒有得到體現(xiàn).對這種非正統(tǒng)做法的解釋是［10］：①本文模型的建立是為了分析大型實際工程，這種方法計算精度雖然不如經(jīng)典塑性力學，但能大大提高計算效率并增強模型的數(shù)值穩(wěn)定性；②式（7）假設彈性應變方向為塑性流動方向可視為塑性流動因子.

式（7）只考慮受壓塑性應變，為了考慮受拉塑性應變，本文采用如下表達式［2］：

從式（9）中可以看出，參數(shù)βt 和βp 分別控制受壓、受拉塑性變量發(fā)展，當受拉損傷發(fā)展即dd＋＞0時，dεp＞0，塑性應變也得到發(fā)展.

為了進一步提高模型的計算效率和穩(wěn)定性，在工程實用中通常將本構曲線下降段變得緩和.這里引入受拉塑性應變在一定程度上就是為了提高模型穩(wěn)定性，圖1 給出了βt 對混凝土材料單軸受拉曲線的影響，可見隨著βt 的增大，材料受拉骨架曲線變得更加緩和，且損傷變量的發(fā)展也變得緩慢.在從后文算例可以看出，這樣處理能在一定程度上提高本構模型的穩(wěn)定性.

圖1 單軸加載下βt 對模型受拉曲線的影響Fig.1 The influence ofβton tension curves of the model under uniaxial loading

2 損傷能量釋放率Perzyna 粘滯規(guī)則化

應變率無關模型損傷能量釋放率閥值直接取為損傷能量釋放率的歷史最大值，對于應變率相關模型，則需要對損傷能量釋放率進行Perzyna粘滯規(guī)則化，文獻［10］建議：

式中：r±表示拉、壓損傷能量釋放率閥值dt表示將r±關于時間求導；μ±為模型參數(shù).

采用無條件穩(wěn)定，且具有二次收斂精度的梯形算法［3］求解式（10），可得：

式 中：Yn＋1／2＝0.5（Yn＋Yn＋1）；rn＋1／2＝0.5（rn＋rn＋1）.

利用Newton－Raphson算法定義f（r）＝－r＋并給出迭代式：

圖2給出了不同應變率作用下混凝土材料的單軸受拉、單軸受壓應力—應變?nèi)€數(shù)值模擬結(jié)果，圖3給出了不同（應變率）加載條件下，混凝土材料抗拉極限承載力、抗壓極限承載力數(shù)值模擬結(jié)果與相應靜力強度比值（在圖3 中用ρ表示）.作為對比，圖3中還給出了實驗結(jié)果［14］，可以看出計算結(jié)果和實驗結(jié)果符合很好.材料參數(shù)為：E＝31 000MPa，泊松比ν＝0.2、單軸抗拉強度ft＝3.38 MPa、單軸抗壓強度fc＝27.6 MPa.對模型其他參數(shù)，采用單軸加載進行標定，其取值分別為：μ＋＝2，μ－＝10 000，a－＝a＋＝4.

模型計算混凝土材料在應變速率分別為0.000 01、0.000 1、0.001、0.01 和0.1 時 的 單 軸 受壓應力—應變滯回曲線如圖4所示.

3 阻尼力的計算

為了在材料本構模型中直接考慮剛度阻尼的影響，本文在建議的彈塑性損傷本構模型中引入剛度阻尼應力.

圖2 單軸加載下的應變率效應Fig.2 Strain－rate effect under uniaxial loads

圖3 模型計算結(jié)果與實驗結(jié)果比較Fig.3 Comparison between calculation and experimental results

Cauchy粘滯阻尼應力σvis可表示為［9］：

彈塑性損傷本構關系為：

式中，εe，εp分別為彈性和塑性應變.則總應力可表示如下：

圖4 不同加載速率下模型應力—應變曲線Fig.4The strain stress curves of the model under different strain rates

4 一致切線模量

在有限元隱式計算中，通常采用Newton－Raphson算法求解平衡方程，此時需要用到一致切線模量，參照文獻［3］，給出本文模型一致切線模量的推導過程.

將式（7）代入式（14）并對ε求導可得：

將式（6）對εe求導可得：

根據(jù)式（11）—（15）給出的算 法可得dr±＝λdY±［3］.其中

將式（17）—（21）代入式（16）即可求解切線剛度模量.

5 算例

應用本文建立的本構模型，對Koyna混凝土重力壩（圖5）進行了數(shù)值模擬，混凝土材料參數(shù)取和文獻［9］一致，分別為：密度ρ0＝2 643kg·m－3，彈性模量E＝31 027MPa，泊松比ν0＝0.2，βp＝0.5，單軸抗拉強度ft＝2.9 MPa，單軸抗壓強度fc＝24.1 MPa.材料阻尼取其第一振型臨界阻尼的3%.在分析中，壩體和基礎之間假定為剛性連接，混凝土壩體部分網(wǎng)格有760個4節(jié)點平面應力單元組成，地震動引起的水對壩體的壓力則通過附加質(zhì)量法模擬.Koyna大壩在自重及水壓力下的自振頻率見表1.

圖5 Koyna混凝土重力壩模型圖Fig.5 Koyna concrete dam

表1 Koyna大壩在自重及水壓力下的自振頻率Tab.1 The nature frequency of Koyna Dam

從圖6中可以看出本文模型計算位移反應和文獻［9］分析結(jié)果符合較好，本文模型分析結(jié)果略大.本文在分析中發(fā)現(xiàn)如果不考慮剛度阻尼，分析到后期容易出現(xiàn)不收斂.加入剛度阻尼后，本構模型的穩(wěn)定性得到顯著提高，可見阻尼對結(jié)構的能量耗散作用能改善隱式動力計算的收斂性.表2給出不同情況下的計算時間，從中可知加入受拉塑性應變也能在一定程度上縮短計算時間，提高模型的穩(wěn)定性，而對計算結(jié)果影響不大，本文認為這是因為考慮受拉塑性應變后材料受拉下降段變得緩和，且受拉損傷演化得到減緩從而提高模型數(shù)值穩(wěn)定性.在實際分析中，混凝土的受拉性能對分析結(jié)果影響很小，但對計算收斂性的影響卻很大，受拉塑性應變的引入能有效提高模型的穩(wěn)定性和收斂性，從而縮短計算時間.

圖6 率不相關模型計算壩頂相對位移時程比較Fig.6 The comparison of top dam displacements with rate independent model

表2 Koyna大壩在不同情況下所需的計算時間Tab.2 The cost time of Koyna Dam under different conditions

從圖7可知，考慮應變率效應后結(jié)構頂層位移反應峰值比不考慮應變率效應略大，這點和文獻［2］分析結(jié)果一致.從表2可以看出，考慮應變率效應后計算時間進一步縮短，可見考慮率應變在降低材料損傷程度的同時也提高了模型的計算效率和穩(wěn)定性.

圖7 應變率效應對壩頂相對位移的影響Fig.7 The effect of strain rates on dam dop displacement

圖8為三種情況下（βt 均取0.5）模型損傷分布圖，在不考慮剛度阻尼和應變率效應情況下加載到后期，由于模型損傷過度發(fā)展導致計算不收斂，本文在這里給出4.727s時的損傷分布圖，其他兩種情況均取加載結(jié)束時（10s）的損傷分布.從圖8中可以看出剛度阻尼、應變率效應對材料損傷發(fā)展有顯著影響，都能延緩材料損傷發(fā)展程度從而提高計算穩(wěn)定性，如果在實際結(jié)構的計算分析中不考慮其影響會導致結(jié)構損傷過快發(fā)展導致計算結(jié)果失真，有時甚至會出現(xiàn)計算不收斂.

圖8 表2所示不同情況下結(jié)構損傷分布圖Fig.8 The structure damage distribution under different conditions as shown in Table 2

6 結(jié)論

本文在提出的靜力彈塑性損傷本構模型的基礎上，通過將損傷演化粘滯化使得模型能考慮應變加載速率的影響.在模型中引入剛度阻尼應力，使其能直接在材料層次考慮剛度阻尼耗能作用.對模型進行改進使其能考慮受拉塑性應變.將建立本構模型在ABAQUS中二次開發(fā)，通過Koyna重力壩動力隱式分析算例表明：剛度阻尼的能量耗散作用能延緩結(jié)構損傷發(fā)展程度，顯著提高隱式動力計算的收斂性；混凝土受拉性能對計算穩(wěn)定性有一定影響，引入受拉塑性應變后模型單軸受拉曲線下降段變得緩和，同時減緩受拉損傷變量的演化，在一定程度上增強模型的穩(wěn)定性，減少計算時間；在模型中考慮應變率效應后對結(jié)構位移反應有一定的影響，同時能進一步提高模型的穩(wěn)定性.

［1］ Kanehi M B，Zienkiewicz O C，Owen D R J.The visco－plastic approach to problems of plasticity and creep involving geometric nonlinear effects［J］.International Journal for Numerical Methods in Engineering，1978，12：169.

［2］ 吳建營，李杰.混凝土彈塑性損傷本構關系統(tǒng)一模型［J］.建筑科學與工程，2005，22（4）：15.WU Jianying，LI Jie.Unified elastoplastic damage model for concrete［J］.Journal of Architecture and Civil Engineering，2005，22（4）：15.

［3］ 吳建營，李杰.考慮應變率效應的混凝土動力彈塑性損傷本構模型［J］.同濟大學學報：自然科學版，2006，34（11）：1427.WU Jianying，LI Jie.Elastoplastic damage constitutive model for concrete considering strain rate effect under dynamic loading［J］.Journal of Tongji University：Natural Science，2006，34（11）：1427.

［4］ Craig R R J.結(jié)構動力學［M］.常嶺，李振邦譯.北京：人民交通出版社，1996.Craig R R J.Structural dynamics［M］.Translated by CHANG Ling，LI Zhenbang.Beijing：China Communications Press，1996.

［5］ H?u?ler－Combe U，Kitzig M.Modeling of concrete behavior under high strain rates with inertially retarded damage［J］.International Journal of Impact Engineering，2009，36：1106.

［6］ Pandey A K，Kumar R，Paul D K，et al.Strain rate model for dynamic analysis of reinforced concrete structures［J］.Journal of Structural Engineering，2006，132（1）：393.

［7］ Barpi F.Impact behaviour of concrete：a computational approach［J］.Engineering Fracture Mechanics，2004，71（2）：197.

［8］ Suffis A，Lubrecht T A A，Combescure A.Damage model with delay effect：analytical and numerical studies of the evolution of the characteristic damage length［J］.International Journal of Solids and Structures，2003，40：346.

［9］ Lee J，F(xiàn)anves G L.A plastic－damage concrete model for earthquake analysis of dams［J］.Earthquake Engineenng and Structural Dynamics，1998，27（9）：937.

［10］ Faria R，Oliver J，Cervera M.A strain－based plastic viscousdamage model for massive concrete structures ［J］.International Journal of Solids and Structures，1998，35（14）：1533.

［11］ QI Hu，LI Yungui，LüXilin.A practical elastic plastic damage model for concrete［J］.Advanced Materials Research，2011，243－249：313.

［12］ 吳建營，李杰.反映阻尼影響的混凝土彈塑性損傷本構模型［J］.工程力學，2006，23（11）：116.WU Jianying，LI Jie.Elastoplastic damage model for concrete considering damping effects［J］.Engineering Mechanics，2006，23（11）：116.

［13］ Faria R，Vila Pouca N，Delgado R.Seismic benchmark of a R／C wall：numerical simulation and experimental validation［J］.Journal of Earthquake Engineering，2002，6（4）：473.

［14］ Suaris W，Shah S.Rate－sensitive damage theory for brittle solids［J］.ASCE，Journal of Engineering Mechanics，1984，110（6）：985.

［15］ Chopra A K，Chakrabarti P.The Koyna earthquake and the damage to Koyna dam［J］.Bulletin of the Seismological Society of America，1973，63（2）：381.