齊 虎，李云貴，呂西林

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

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

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

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

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

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

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

1 率不相關(guān)彈塑性損傷本構(gòu)模型

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

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

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

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

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

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

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

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

加權(quán)損傷變量：

1.2 塑性應(yīng)變計算

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

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

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

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

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

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

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

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

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

式中：r±表示拉、壓損傷能量釋放率閥值dt表示將r±關(guān)于時間求導(dǎo)；μ±為模型參數(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給出了不同應(yīng)變率作用下混凝土材料的單軸受拉、單軸受壓應(yīng)力—應(yīng)變?nèi)€數(shù)值模擬結(jié)果，圖3給出了不同（應(yīng)變率）加載條件下，混凝土材料抗拉極限承載力、抗壓極限承載力數(shù)值模擬結(jié)果與相應(yīng)靜力強(qiáng)度比值（在圖3 中用ρ表示）.作為對比，圖3中還給出了實驗結(jié)果［14］，可以看出計算結(jié)果和實驗結(jié)果符合很好.材料參數(shù)為：E＝31 000MPa，泊松比ν＝0.2、單軸抗拉強(qiáng)度ft＝3.38 MPa、單軸抗壓強(qiáng)度fc＝27.6 MPa.對模型其他參數(shù)，采用單軸加載進(jìn)行標(biāo)定，其取值分別為：μ＋＝2，μ－＝10 000，a－＝a＋＝4.

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

3 阻尼力的計算

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

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

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

Cauchy粘滯阻尼應(yīng)力σvis可表示為［9］：

彈塑性損傷本構(gòu)關(guān)系為：

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

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

4 一致切線模量

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

將式（7）代入式（14）并對ε求導(dǎo)可得：

將式（6）對εe求導(dǎo)可得：

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

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

5 算例

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

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

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

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

圖6 率不相關(guān)模型計算壩頂相對位移時程比較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可知，考慮應(yīng)變率效應(yīng)后結(jié)構(gòu)頂層位移反應(yīng)峰值比不考慮應(yīng)變率效應(yīng)略大，這點和文獻(xiàn)［2］分析結(jié)果一致.從表2可以看出，考慮應(yīng)變率效應(yīng)后計算時間進(jìn)一步縮短，可見考慮率應(yīng)變在降低材料損傷程度的同時也提高了模型的計算效率和穩(wěn)定性.

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

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

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

6 結(jié)論

本文在提出的靜力彈塑性損傷本構(gòu)模型的基礎(chǔ)上，通過將損傷演化粘滯化使得模型能考慮應(yīng)變加載速率的影響.在模型中引入剛度阻尼應(yīng)力，使其能直接在材料層次考慮剛度阻尼耗能作用.對模型進(jìn)行改進(jìn)使其能考慮受拉塑性應(yīng)變.將建立本構(gòu)模型在ABAQUS中二次開發(fā)，通過Koyna重力壩動力隱式分析算例表明：剛度阻尼的能量耗散作用能延緩結(jié)構(gòu)損傷發(fā)展程度，顯著提高隱式動力計算的收斂性；混凝土受拉性能對計算穩(wěn)定性有一定影響，引入受拉塑性應(yīng)變后模型單軸受拉曲線下降段變得緩和，同時減緩受拉損傷變量的演化，在一定程度上增強(qiáng)模型的穩(wěn)定性，減少計算時間；在模型中考慮應(yīng)變率效應(yīng)后對結(jié)構(gòu)位移反應(yīng)有一定的影響，同時能進(jìn)一步提高模型的穩(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］ 吳建營，李杰.混凝土彈塑性損傷本構(gòu)關(guān)系統(tǒng)一模型［J］.建筑科學(xué)與工程，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］ 吳建營，李杰.考慮應(yīng)變率效應(yīng)的混凝土動力彈塑性損傷本構(gòu)模型［J］.同濟(jì)大學(xué)學(xué)報：自然科學(xué)版，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é)構(gòu)動力學(xué)［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］ 吳建營，李杰.反映阻尼影響的混凝土彈塑性損傷本構(gòu)模型［J］.工程力學(xué)，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.