蘇小卒，王 偉

(1.同濟(jì)大學(xué)土木工程學(xué)院，上海 200092；2.紹興文理學(xué)院土木工程學(xué)院，浙江 312000)

常見的工程材料(如混凝土、鋼材)為依賴應(yīng)變歷史材料(SHDM)，此類材料的正應(yīng)變?cè)隽窟^(guò)程(應(yīng)變加載)中的切線模量，不同于負(fù)應(yīng)變?cè)隽窟^(guò)程(應(yīng)變卸載)中的切線模量，即材料點(diǎn)應(yīng)變?cè)诩虞d和卸載(下文簡(jiǎn)稱“應(yīng)變加卸載”)時(shí)具有不同的應(yīng)力變化規(guī)律。當(dāng)不考慮其他非線性因素時(shí)，SHDM結(jié)構(gòu)(SHDM制成的結(jié)構(gòu))的有限元靜力問題是一個(gè)材料非線性問題。此靜力問題完成求解時(shí)，會(huì)獲得滿足靜力平衡的荷載-位移曲線(以下簡(jiǎn)稱平衡路徑)。平衡路徑具有硬化段，但當(dāng)材料兼具應(yīng)變軟化性質(zhì)時(shí)(如混凝土)，它可能還具有極值點(diǎn)、軟化段。此靜力問題進(jìn)行求解時(shí)：解法一般采用增量迭代法，即總體上采用增量法求解，增量步采用迭代法求解[1]；迭代時(shí)，為保證收斂解為真實(shí)靜力解，用以更新單元內(nèi)高斯積分點(diǎn)(以下簡(jiǎn)稱高斯點(diǎn))應(yīng)力的應(yīng)變須僅經(jīng)歷過(guò)加載或卸載。

隱式靜力算法[2](ISA)、動(dòng)力松弛法[3](DRM)是兩種求解靜力問題增量步的迭代算法。對(duì)迭代中的應(yīng)力更新步驟，必須采用增量步內(nèi)的結(jié)點(diǎn)位移累積增量計(jì)算高斯點(diǎn)應(yīng)力增量[1]，否則會(huì)引入虛假應(yīng)變歷史，從而使得滿足靜力平衡的迭代收斂解為虛假靜力解[4]。此外，靜力問題也可采用顯式擬靜力算法(EQSA)近似求解[5]。

采用ISA迭代求靜力增量步時(shí)，需在預(yù)判高斯點(diǎn)應(yīng)變加卸載狀態(tài)后進(jìn)行總體剛度矩陣K的組裝和求逆，但迭代是否收斂依賴于材料特征和K特征。平衡路徑硬化段內(nèi)，K正定；極值點(diǎn)處，K奇異；軟化段內(nèi)，K負(fù)定。荷載控制類ISA在ABAQUS有限元軟件中為一般非線性問題的推薦算法[6]，但它不能求解K非正定的靜力問題[7]。位移控制類ISA、整體弧長(zhǎng)類ISA不適合求解應(yīng)變軟化結(jié)構(gòu)的靜力問題[2－8]，前者還難以合適選定控制位移分量[9]，后者還難以跨越材料非線性引起的極值點(diǎn)[10]。局部弧長(zhǎng)類ISA目前算例主要考慮受拉應(yīng)變軟化情形[11－13]。因此，對(duì)于復(fù)雜靜力分析問題，文獻(xiàn)[2]建議先將其轉(zhuǎn)化為動(dòng)力問題，再用EQSA給出靜力近似解。

EQSA本質(zhì)上是動(dòng)力問題的顯式時(shí)間積分算法；當(dāng)通過(guò)低加載速率使動(dòng)力系統(tǒng)的慣性效應(yīng)滿足擬靜力判據(jù)時(shí)，顯式時(shí)間積分算法可視為EQSA[2]。采用EQSA近似求解靜力問題時(shí)，會(huì)在真實(shí)時(shí)間域向量式地更新結(jié)點(diǎn)位移a。因此，EQSA的優(yōu)點(diǎn)為：無(wú)需預(yù)判高斯點(diǎn)應(yīng)變加卸載狀態(tài)和構(gòu)造K[6]。但是，EQSA的局限性為：① 計(jì)算成本具有隨結(jié)點(diǎn)質(zhì)量的減小呈指數(shù)比例增大的規(guī)律[8]；② 合理確定粘滯阻尼系數(shù)時(shí)的經(jīng)驗(yàn)依賴性；③a依時(shí)間序列產(chǎn)生，由此給出的近似靜力解中存在虛假應(yīng)變歷史[14]。

DRM是一種靜力分析算法，其首先構(gòu)造一個(gè)包含粘滯阻尼[3]或動(dòng)態(tài)阻尼[15]的虛擬動(dòng)力系統(tǒng)，然后通過(guò)此系統(tǒng)的動(dòng)能阻尼耗散計(jì)算靜力解。采用DRM求解靜力增量步時(shí)：在虛擬時(shí)間域采用顯式時(shí)間積分算式進(jìn)行a的向量式更新，從而保留EQSA的優(yōu)點(diǎn)；可采用結(jié)點(diǎn)虛擬質(zhì)量緩解EQSA的第一個(gè)局限性，采用動(dòng)態(tài)阻尼規(guī)避第二個(gè)局限性。此外，DRM的另一優(yōu)點(diǎn)為：若SHDM結(jié)構(gòu)兼具應(yīng)變軟化性質(zhì)(如混凝土結(jié)構(gòu))，它能解決“應(yīng)變軟化結(jié)構(gòu)因靜力定解問題失去適定性”而導(dǎo)致的迭代不收斂，相關(guān)論證可見作者的先前工作[16]。Rezaieepajand等[17]進(jìn)行了桁架、框架結(jié)構(gòu)靜力問題在施加不同虛擬阻尼時(shí)的DRM分析，結(jié)果表明kDRM(動(dòng)態(tài)阻尼DRM)的計(jì)算成本最少。DRM的結(jié)構(gòu)分析應(yīng)用始于20世紀(jì)五六十年代[18－19]。而后應(yīng)用范圍逐漸擴(kuò)大，詳見文獻(xiàn)[20－27]。但是，SHDM結(jié)構(gòu)的DRM分析報(bào)道較少見諸于文獻(xiàn)。

本文采用kDRM求解SHDM結(jié)構(gòu)的靜力平衡路徑。迭代求解增量步時(shí)，為了避免虛假應(yīng)變歷史效應(yīng)，構(gòu)造了高斯點(diǎn)應(yīng)力更新的非線性彈性增量算法(NEIA)——盡管材料一般不是非線性彈性材料，但在一個(gè)增量步內(nèi)，將材料視為非線性彈性材料。全文首先給出問題的描述，然后詳述問題的求解，接著描述三個(gè)算例及其計(jì)算結(jié)果，其中兩個(gè)SHDM結(jié)構(gòu)模型兼具應(yīng)變軟化特征，最后給出討論、結(jié)論。

1 問題的描述

本文討論的依賴應(yīng)變歷史材料(SHDM)是指材料點(diǎn)應(yīng)變加卸載時(shí)應(yīng)力響應(yīng)不同的材料。在一維情況下，SHDM滿足[1]：

式中：σ為應(yīng)力標(biāo)量；ε為應(yīng)變；ε0為材料點(diǎn)在應(yīng)變加卸載過(guò)程中任意時(shí)刻的應(yīng)變。在多軸應(yīng)力情況下，只要應(yīng)力、應(yīng)變的某一分量(或等效應(yīng)力、等效應(yīng)變)滿足式(1)時(shí)，即為本文所指的SHDM。土木工程材料中常用的塑性本構(gòu)模型、損傷本構(gòu)模型為SHDM的本構(gòu)模型。

結(jié)構(gòu)有限元分析時(shí)，自由度分為兩類：其一，施加外部結(jié)點(diǎn)荷載的加載自由度；其二，未加載自由度。本文在進(jìn)行結(jié)構(gòu)非線性有限元靜力分析時(shí)：結(jié)構(gòu)為邊界受到充分約束而不具有剛體位移的幾何不變體系；主要采用比例加荷載方式。由此，基于全部自由度建立的待解非線性平衡方程組為：

式中：a為結(jié)點(diǎn)位移向量；p為內(nèi)部結(jié)點(diǎn)荷載向量；λ為外部結(jié)點(diǎn)荷載系數(shù)；fref為外部結(jié)點(diǎn)荷載參照向量；ψ為結(jié)點(diǎn)殘值荷載向量；a、p、fref均為r維向量(有限元模型的自由度數(shù))。p按式計(jì)算[1]，其中B為nσ×r應(yīng)變矩陣，σ為nσ×1應(yīng)力向量，nσ為σ內(nèi)的元素?cái)?shù)，Ω 為固體所占空間。

本文采用加荷載方式或加位移方式求解式(2)所描述的問題(以下簡(jiǎn)稱問題2)。求解方式為加荷載時(shí)：已知量為λ；待求量為a。求解方式為加位移時(shí)：已知量為被選取的加載位移aI(下標(biāo)表示第I個(gè)自由度)；待求量為aII——a中除aI外的其他元素，以及λ；記a=[aI,aII]T。

2 問題的求解

本文中：① 問題2總體上用荷載或位移增量法求解；② 增量步用kDRM求解——結(jié)點(diǎn)位移靜態(tài)解借助對(duì)虛擬動(dòng)力系統(tǒng)施加足夠次的動(dòng)態(tài)阻尼而獲得；③ 通過(guò)高斯點(diǎn)應(yīng)力更新的NEIA使得靜態(tài)解符合靜力解特征——高斯點(diǎn)應(yīng)力由真實(shí)應(yīng)變歷史更新。本文中，施加動(dòng)態(tài)阻尼的含義為：置零“有限元離散系統(tǒng)結(jié)點(diǎn)總動(dòng)能極大值時(shí)刻的”結(jié)點(diǎn)速度。

2.1 靜力問題的增量解法

總體上，問題2用荷載或位移增量法求解。求解時(shí)，整個(gè)平衡路徑連續(xù)地分解成序列{Rn}，其中Rn(n=1,2,…)為第n個(gè)增量步。Rn的解為an、狀態(tài)量Cn、λn。an、Cn中的εn、σn，以及λn可寫為：

式中：Δan為結(jié)點(diǎn)位移增量；Δεn為應(yīng)變?cè)隽浚沪う襫為應(yīng)力增量；Δλn為外部結(jié)點(diǎn)荷載系數(shù)增量。

加荷載增量或者加位移增量求解Rn。默認(rèn)先加荷載增量求解R1；R1施加前的初值λ0、a0、C0為已知量，λ0和C0內(nèi)的σ0滿足式(2)。Rn求解要點(diǎn)如下：

1) 已知量：①Rn-1的解λn-1、an-1、Cn-1；② 加荷載時(shí)施加的固定荷載增量系數(shù)Δλ，或加位移時(shí)對(duì)第I個(gè)自由度施加的固定位移增量ΔaI，即：

2) 未知量：① 加荷載時(shí)為an、Cn，其中an為基本未知量；② 加位移時(shí)為an,II、Δλn、Cn，其中an,II為基本未知量。

3) 未知量用kDRM迭代求解，詳見2.2節(jié)。

4) 記i為kDRM求解的迭代步序號(hào)，對(duì)于求得的

① 若滿足收斂判據(jù)

式中，c1為限值，則得問題的解為加位移時(shí)，若此解符合增量變換規(guī)則：

則表明此解處在平衡路徑硬化段，則后面使用加荷載方式求解Rn+1。文獻(xiàn)[22]在采用DRM求解張拉結(jié)構(gòu)平衡路徑時(shí)，即采用式(3)所示的收斂判據(jù)，并取c1=0.001，此判據(jù)和限值符合文獻(xiàn)[1]的建議。

② 加荷載時(shí)，若滿足發(fā)散判據(jù)

式中，c2為限值，則說(shuō)明λnfref超越或臨近平衡路徑的峰值點(diǎn)，此時(shí)更換為加位移方式后重新求解Rn。

2.2 增量步的kDRM解法

2.2.1 kDRM的控制方程

kDRM求解靜力問題的策略為[28]：在虛擬時(shí)間域監(jiān)測(cè)無(wú)阻尼質(zhì)點(diǎn)系統(tǒng)的動(dòng)能，在各結(jié)點(diǎn)靜力平衡前，循環(huán)執(zhí)行“在動(dòng)能極大值時(shí)刻先置零結(jié)點(diǎn)速度，再重啟無(wú)阻尼運(yùn)動(dòng)”。kDRM求解Rn時(shí)，記第k個(gè)施加動(dòng)態(tài)阻尼循環(huán)為Sk，則其內(nèi)結(jié)點(diǎn)運(yùn)動(dòng)問題寫為：

1)tk為在t＞tk-1上的首個(gè)系統(tǒng)動(dòng)能的極大值時(shí)間；tk也為循環(huán)Sk+1的起始時(shí)間。

2)tk-1為循環(huán)Sk-1的結(jié)束時(shí)間；當(dāng)k=1時(shí)，tk-1=0，此時(shí)即為增量步Rn的加荷載(位移)時(shí)間。

3)Mn為元素是已知常量的矩陣，2.2.2小節(jié)給出了Mn的具體構(gòu)造方法。

4)Cn(t)為元素是未知變量εn(t)、σn(t)、sn(t)的集合。εn(t)=Ban(t)；sn(t)的自變量為εn(t)；σn(t)的自變量為sn(t)及εn(t)。所以，Cn(t)是an(t)的函數(shù)。

5)pn(t)為元素是未知變量的向量，算式為[8]：

式中：e表示單元編號(hào)；ne表示系統(tǒng)內(nèi)的單元數(shù)量；p{e}、B{e}、σ{e}、Ω{e}分別為單元e的內(nèi)部結(jié)點(diǎn)荷載、應(yīng)變矩陣、應(yīng)力向量、所占空間；p{e}、B{e}、σ{e}分別為的矩陣(為單元e的結(jié)點(diǎn)自由度數(shù))；T{e}為單元e的轉(zhuǎn)換矩陣，尺寸為(r為系統(tǒng)全部自由度數(shù))。因?yàn)镃(t)是nan(t)的函數(shù)，所以pn(t)也為an(t)的函數(shù)。

6) 對(duì)于λn(t)：① 加荷載時(shí)為已知的常量

② 加位移時(shí)為未知的變量

式中：fref,I、pn,I為fref、pn(t)中與第I個(gè)自由度(加位移自由度)對(duì)應(yīng)的元素值；λn(t)是an(t)的函數(shù)。

由本小節(jié)分析可知：kDRM求解由問題9所描述的Rn的關(guān)鍵是，給出基本未知量an(t)的數(shù)學(xué)算法。因?yàn)閱栴}9是一個(gè)強(qiáng)非線性動(dòng)力問題，所以求解an(t)時(shí)，宜選用顯式時(shí)間積分算法[29]。

2.2.2 控制方程的中心差分算法求解

本文選用中心差分算法(CFDM)求解問題9；求解時(shí)，ih時(shí)刻的迭代式為[30]：

式中：ik-1為Sk的起始時(shí)刻序號(hào)；h為離散時(shí)間域的固定時(shí)間增量；的計(jì)算式為：

其中，按式(11)或式(12)確定。用式(13)迭代求解問題9中的基本未知量an(t)時(shí)，需解決三個(gè)問題：

1) 確定Mn及h。一般采用CFDM的穩(wěn)定條件構(gòu)造Mn、h。Barnes給出的CFDM穩(wěn)定條件[22]為：

式中：下標(biāo)表示結(jié)點(diǎn)序號(hào)；mj為結(jié)點(diǎn)質(zhì)量；kjmax為結(jié)點(diǎn)沿各坐標(biāo)正負(fù)方向位移剛度最大值。本文估算kjmax的方式為：首先給結(jié)點(diǎn)施加微小位移，然后計(jì)算相應(yīng)結(jié)點(diǎn)荷載，最后算出kjmax。本文取

式中：c3為系數(shù)，值大于1；mmin是結(jié)點(diǎn)質(zhì)量下限值，其作用是避免結(jié)點(diǎn)負(fù)剛度而構(gòu)造出負(fù)的結(jié)點(diǎn)質(zhì)量。本文按式(16)逐結(jié)點(diǎn)算得mj后構(gòu)造Mn。

3) 確定Sk的結(jié)束時(shí)刻。時(shí)刻(i-1/2)h、(i-3/2)h的動(dòng)能分別為若滿足：

則認(rèn)為(i-1)h時(shí)離散質(zhì)點(diǎn)系統(tǒng)出現(xiàn)了動(dòng)能極大值，Sk的結(jié)束時(shí)間以及Sk+1的起始時(shí)間取為(i-1)h，以為初始條件啟動(dòng)Sk+1。

至此，基于式(13)能解出問題9中的基本未知量an(t)，其中：① 加荷載時(shí)，由式(11)知，由迭代式(13)直接解出；② 加位移時(shí)，由式(12)、式(13)知，保持不變，由迭代式(13)直接解出，由式(12)解出。

2.2.3 完整的kDRM數(shù)值計(jì)算流程

依據(jù)2.1節(jié)、2.2節(jié)的討論，Rn的完整kDRM求解流程，見表1。

2.3 高斯點(diǎn)應(yīng)力更新：非線性彈性增量算法

依賴應(yīng)變歷史材料(SHDM)在應(yīng)變加載、卸載時(shí)具有不同材料響應(yīng)，故對(duì)結(jié)構(gòu)有限元模型進(jìn)行受力分析時(shí)，s需記錄高斯點(diǎn)加載、卸載時(shí)應(yīng)力-應(yīng)變關(guān)系的狀態(tài)。

結(jié)構(gòu)靜力分析時(shí)，結(jié)果須符合靜力解特征[1]：高斯點(diǎn)應(yīng)力的更新必須基于增量步起始狀態(tài)，即基于而不能基于增量步中間狀態(tài)，即不能基于為此，采用DRM求解增量步的靜力解時(shí)，本文構(gòu)造了高斯點(diǎn)應(yīng)力更新的非線性彈性增量算法——在當(dāng)前增量步內(nèi)，材料視作非線性彈性，從而ih時(shí)刻的由ih時(shí)刻的及增量步起點(diǎn)處的εn-1、sn-1更新，即：

表1 Rn的完整kDRM求解流程Table 1 kDRM algorithm for Rn

式中，fσ為材料點(diǎn)應(yīng)力計(jì)算函數(shù)。按此算法，僅當(dāng)滿足收斂判據(jù)時(shí)才更新C。NEIA不同于結(jié)構(gòu)動(dòng)力分析時(shí)(包括EQSA)的NIA——在當(dāng)前增量步內(nèi)，ih時(shí)刻的由ih時(shí)刻的及(i-1)h時(shí)刻的更新(注：更新)，即：

按此算法，按時(shí)間序列{ih}更新C。

按上述，圖1和圖2示意了依賴應(yīng)變歷史的、一維非彈性本構(gòu)模型的四個(gè)高斯點(diǎn)(1、2、3、4)，在增量步的DRM求解過(guò)程中，采用NEIA、NIA時(shí)的本構(gòu)狀態(tài)遞變歷程。顯然，有限元離散系統(tǒng)靜態(tài)時(shí)，前者為具備真實(shí)應(yīng)變歷史的真實(shí)靜力平衡狀態(tài)，后者可能為具備虛假應(yīng)變歷史的虛假靜力平衡狀態(tài)。

圖1 非線性彈性增量算法下的本構(gòu)遞變Fig.1 A illustration of constitutive-state-changing under nonlinear elastic increment-algorithm

圖2 非線性增量算法下的本構(gòu)遞變Fig.2 A illustration of constitutive-state-changing under nonlinear increment-algorithm

3 算例

下面采用結(jié)點(diǎn)位移更新的kDRM方法，以及分別采用高斯點(diǎn)應(yīng)力更新的NEIA、NIA算法對(duì)兩個(gè)有限元模型進(jìn)行了靜力平衡路徑計(jì)算，計(jì)算結(jié)果證實(shí)了NEIA的優(yōu)越性。

3.1 算例I

模型I為一個(gè)拉壓桿組成的平面屋架(圖3)，編號(hào)1、2、3桿的截面面積A均為1500 mm2，編號(hào)4、5、6、7桿的A均為16000 mm2，其余桿的A均為900 mm2。分析時(shí)，每根桿子視為一個(gè)二結(jié)點(diǎn)線單元，其形函數(shù)采用一次多項(xiàng)式；單元軸向應(yīng)變?chǔ)?(l- lini)/lini，式中l(wèi)、lini分別為已變形的、未變形的軸向長(zhǎng)度；由ε算得單元軸向應(yīng)力σ后，可通過(guò)(結(jié)點(diǎn)的)直接平衡法算得p{e}。

模型I使用《混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范》(GB 50010―2010)C.1節(jié)定義的依賴應(yīng)變歷史的、應(yīng)變硬/軟化的、一維彈塑性本構(gòu)模型σ=σ(ε)，其為奇函數(shù)；當(dāng)ε ≥0時(shí)，σ表達(dá)式為：

式中：E為彈性模量；εy為屈服應(yīng)變；fy為屈服應(yīng)力；k為硬/軟化段斜率(當(dāng)k＞0時(shí)，為應(yīng)變硬化；當(dāng)k=0時(shí)，為理想彈塑性，當(dāng)k＜0時(shí)，為應(yīng)變軟化)；εuy為硬/軟化段起點(diǎn)應(yīng)變；εu為峰值應(yīng)變。各拉壓桿的本構(gòu)模型參數(shù)取值見表2。應(yīng)力-應(yīng)變關(guān)系的示意圖見圖3。

圖3 模型I的幾何信息、離散和本構(gòu)關(guān)系示意Fig.3 The model I: geometry，meshing and constitutive relation illustration

表2 模型I 桿的本構(gòu)關(guān)系參數(shù)及截面面積Table 2 The model I: parameters of constitutive model and cross-section area of bar

模型I的加載結(jié)點(diǎn)編號(hào)、加載自由度見圖3。分析時(shí)：mmin=0.3 kg，c1=2×10-4，c3=4.0，h=0.1 s，加荷載時(shí)Δλ=1；外部結(jié)點(diǎn)荷載參照向量fref、c2、加位移控制時(shí)加位移自由度上施加的固定位移增量ΔaI的配置見表3。計(jì)算中，考慮了幾何非線性，但是未考慮歐拉失穩(wěn)；當(dāng)任一單元內(nèi)出現(xiàn)應(yīng)變大于εu情形時(shí)，終止增量步分析。

3.2 算例II

模型II為一個(gè)懸臂梁(圖4)，厚度為10 mm，自由端承受集中力。此懸臂梁采用形函數(shù)為二次多項(xiàng)式的四結(jié)點(diǎn)面單元進(jìn)行離散，采用4個(gè)高斯點(diǎn)近似計(jì)算p{e}，詳見圖4。

表3 模型I、II、III的分析參數(shù)Table 3 Analysis parameters for model I, II and III

模型II采用描述理想塑性材料的Prandtl-Reuss本構(gòu)模型[31](圖4中三個(gè)灰色單元除外，原因見本節(jié)下文)，相應(yīng)的增量本構(gòu)關(guān)系為：

式中：dσ為應(yīng)力增量；dε為應(yīng)變?cè)隽?；D為材料剛度矩陣。D為：

式中：De、Dp分別為材料彈性、塑性剛度矩陣，具體由彈性模量E、泊松比ν、增量前應(yīng)力σ0構(gòu)造，詳見文獻(xiàn)[32]；為由試探應(yīng)力σtry=σ0+Dedε算得的偏應(yīng)力張量第二不變量，平面應(yīng)力問題中，β為各向同性參數(shù)。本算例中，材料參數(shù)取值為：ν=0.3，E=200000 N·mm-2，β=300 N·mm-2；按此材料參數(shù)計(jì)算，軸向拉壓情形下，材料屈服強(qiáng)度集中力(位置見圖4)作用附近可能存在復(fù)雜的高應(yīng)力，此高應(yīng)力可能超越給定材料參數(shù)的Prandtl-Reuss本構(gòu)模型描述范圍，從而出現(xiàn)分析無(wú)法進(jìn)行的情況；為了避免此情況，本算例中，集中力作用附近的三個(gè)單元(圖4中的灰色單元)采用各向同性彈性本構(gòu)模型，構(gòu)造此模型的ν=0.3，E=200000 N·mm-2。

圖4 模型II的幾何信息和離散Fig.4 The model II: geometry and meshing

模型II的加載結(jié)點(diǎn)編號(hào)、加載自由度見圖4。分析時(shí)：c1=2×10-4，c3=4.0，h=0.1s，mmin=0.3 kg，加荷載時(shí)Δλ=1；fref、c2、ΔaI的配置見表3。當(dāng)1號(hào)結(jié)點(diǎn)的y向位移a1y＞80 mm時(shí)，終止增量步分析。

3.3 算例III

模型III是一個(gè)素混凝土構(gòu)件(圖5)，厚度為150 mm，承受軸心壓力。有限元分析時(shí)，采用與模型Ⅱ相同的單元類型、單元形函數(shù)、p{e}的計(jì)算方式。

圖5 模型III的幾何信息、離散和本構(gòu)關(guān)系示意Fig.5 The model III: geometry，meshing and constitutive relation illustration

模型III使用Mazars提出的單標(biāo)量彈性損傷本構(gòu)模型[32－33]：

式中：Dini是由彈性模量E、泊松比ν構(gòu)造的材料初始剛度矩陣(各向同性)；標(biāo)量ω是損傷變量。ω的計(jì)算公式為：

式中：ωmax是材料點(diǎn)所達(dá)到過(guò)的最大損傷值；α+ω++α-ω-是此材料點(diǎn)的即時(shí)損傷值，其中α+、α-是依賴于ε的系數(shù)，ω-、ω+是受壓、受拉損傷變量。α+、α-的算式見文獻(xiàn)[33]。ω+、ω-的算式為：

式中：Kini、A+、B+、A-、B-均為材料常數(shù)，其值通過(guò)實(shí)驗(yàn)測(cè)定；是值依賴于ε的等效應(yīng)變，計(jì)算式參見文獻(xiàn)[33]。本算例中，材料參數(shù)為：ν=0.2，E=25500 N·mm-2； A+=0.95，A-=1.38，B+=11500，B-=2000，Kini= 9.8×10-5；利用上述參數(shù)可給出材料的單軸受壓應(yīng)力應(yīng)變關(guān)系(圖5為示意圖)，其中峰值強(qiáng)度f(wàn)c=24.4 N·mm-2，與fc對(duì)應(yīng)的應(yīng)變?chǔ)與=1.75×10-3。

模型III的加載結(jié)點(diǎn)編號(hào)、加載自由度見圖5。模型III是一個(gè)受軸心壓力的混凝土構(gòu)件，使用兩種方法分析：第一種方法(解法1)，使用比例加荷載增量分析完整平衡路徑；第二種方法(解法2)，使用比例加荷載增量分析峰值點(diǎn)前平衡路徑，而使用比例加等位移增量分析峰值點(diǎn)后平衡路徑。比例加等位移增量時(shí)，即將相等位移增量施加于加載自由度，此時(shí)平衡方程組為ψ(fn,an)=fn-p(an)=0，式中fn是外部結(jié)點(diǎn)荷載。求解此方程組時(shí)：非加載自由度上的位移增量Δaα為未知量，而荷載增量Δfα=0；加載自由度上的荷載增量位移增量Δaβ為已知量，而Δfβ為未知量；由kDRM求解Δaα、Δfβ的計(jì)算流程可見文獻(xiàn)[34]。無(wú)論是解法1，還是解法2，本文都采用kDRM求解增量步，且kDRM分析參數(shù)中的c1、c3、h、mmin、Δλ取值相同，即c1=2×10-4，c3=4.0，h=0.1 s，mmin=0.3 kg，Δλ=1。采用解法1時(shí)，kDRM分析時(shí)的其余參數(shù)配置見表3。采用解法2時(shí)，僅對(duì)一組參數(shù)配置下的模型IIID進(jìn)行分析，分析時(shí)：c2及求解平衡路徑硬化段時(shí)的fref見表3；求解平衡路徑軟化段時(shí)，對(duì)加載自由度施加的結(jié)點(diǎn)位移增量滿足Δa1y=Δa2y=Δa3y=Δa4y=Δa5y=0.015 mm。當(dāng)?shù)鈺r(shí)，終止計(jì)算。由計(jì)算后的數(shù)據(jù)分析發(fā)現(xiàn)：解法2獲得的解答同樣滿足表3給出的比例加荷載要求。

3.4 計(jì)算結(jié)果及分析

1) 模型I：圖6為1號(hào)結(jié)點(diǎn)y向位移-荷載關(guān)系曲線(f1y-a1y曲線)，其標(biāo)識(shí)了模型I加載的硬化階段、屈服階段(限模型號(hào)IA、IB和IC)、軟化階段；圖7為IA情況下4號(hào)、5號(hào)桿的ε-a1y的關(guān)系曲線。模型II：圖8為f1y-a1y關(guān)系曲線，其標(biāo)識(shí)了模型II加載的硬化階段、屈服階段；圖9為IIA情況下的1號(hào)單元、2號(hào)單元的1號(hào)高斯點(diǎn)的關(guān)系曲線和關(guān)系曲線。模型III：圖10為f3y-a3y關(guān)系曲線，其標(biāo)識(shí)了模型III加載的硬化階段、軟化階段；圖11為IIIB、IIID情況下1號(hào)、2號(hào)單元的εy-a3y的關(guān)系曲線。對(duì)于圖10：圖形★標(biāo)識(shí)IIIA的最終一個(gè)收斂解(圖中未標(biāo)識(shí)臨近★的IIIB、IIIC的最終一個(gè)收斂解)；圖形◆標(biāo)識(shí)IIIa、IIIb、IIIc的最終一個(gè)收斂解；圖形▼標(biāo)識(shí)“增量步迭代終止條件(式(26))滿足、但收斂條件(式(6))不滿足”的迭代解；圖例IIID標(biāo)示采用解法2和NEIA算得的平衡路徑。對(duì)于圖11：圖形●標(biāo)識(shí)最終一個(gè)收斂解；圖形▲或▼標(biāo)識(shí)“增量步迭代終止條件(式(26))滿足、但收斂條件(式(6))不滿足”的迭代解。

2) 從圖6、圖8、圖10可知：① 采用NEIA時(shí)，不同的增量幅度(IA、IB、IC，IIA、IIB、IIC，IIIA、IIIB、IIIC)給出相同的平衡路徑；② 采用NIA時(shí)，不同的增量幅度(Ia、Ib、Ic，IIa、IIb、IIc，IIIa、IIIb、IIIc)給出不同的平衡路徑，增量幅度越小，算得的平衡路徑越接近NEIA算得的平衡路徑。

圖6 模型I的f1y-a1y曲線Fig.6 The model I: f1y-a1ycurve

圖7 模型I的ε-a1y曲線 (IA)Fig.7 The model I: ε-a1ycurve (IA)

圖8 模型II的f1y-a1y曲線Fig.8 The model II f1y-a1ycurve

3) 對(duì)于模型I，由圖7可知(結(jié)合表2)：4號(hào)桿、5號(hào)桿的材料性能階段與相應(yīng)的平衡路徑階段(圖6)相關(guān)，如結(jié)構(gòu)的軟化階段始于5號(hào)單元的應(yīng)變軟化，且同時(shí)伴隨著4號(hào)桿的應(yīng)變卸載，此現(xiàn)象即為形變局部化。

4) 對(duì)于模型II，對(duì)于x=121 mm處的梁橫截面S，假定S應(yīng)變保持平面，則由材料力學(xué)公式得：① 當(dāng)S僅有最外邊緣纖維應(yīng)力等于fy時(shí)(狀態(tài)e)，f1y=59.7 kN；② 當(dāng)S的所有纖維應(yīng)力等于fy時(shí)(狀態(tài)p)，f1y=89.6 kN。依據(jù)圖8、圖9，并結(jié)合圖4給出的單元1的1號(hào)高斯點(diǎn)、單元2的2號(hào)高斯點(diǎn)的y坐標(biāo)，得：① 當(dāng)a1y≈24 mm時(shí)，S剛越過(guò)狀態(tài)e，此時(shí)結(jié)構(gòu)剛度開始降低，f1y≈68.0 kN，此與理論值接近；② 當(dāng)a1y=49 mm時(shí)，S臨近狀態(tài)p，此時(shí)結(jié)構(gòu)剛度臨近零，f1y≈92.0 kN，此與理論值接近。

圖10 模型III的f3y-a3y曲線Fig.10 The model III: f3y-a3ycurve

圖11 模型III的εy-a3y曲線(IIIB、IIID)Fig.11 The model III: εy-a3ycurve(IIIB、IIID)

4 討論

1) 當(dāng)加載結(jié)點(diǎn)的增量幅度超過(guò)一定范圍時(shí)，NEIA可能會(huì)給出不真實(shí)的平衡路徑，如算例I的ID的平衡路徑不同于IA、IB、IC的平衡路徑。從計(jì)算所得的應(yīng)變數(shù)據(jù)可知：ID的5號(hào)桿跨越了應(yīng)變屈服階段。此導(dǎo)致f1y-a1y曲線不存在屈服階段、算得的平衡路徑偏離真實(shí)平衡路徑。由此可知，即使采用NEIA求解平衡路徑，也應(yīng)避免采用過(guò)大的增量幅度。

2) 盡管算例Ⅰ采用了虛構(gòu)的本構(gòu)關(guān)系，算例Ⅱ采用了較為符合實(shí)際的鋼材本構(gòu)關(guān)系，算例III采用的單標(biāo)量彈性損傷本構(gòu)模型不能很好反映混凝土多軸應(yīng)力狀態(tài)下的本構(gòu)關(guān)系[35]，但因?yàn)镈RM所構(gòu)造的虛擬動(dòng)力過(guò)程均符合熱力學(xué)第一定律，動(dòng)能的耗散必將使得動(dòng)力解答逼近具有真實(shí)應(yīng)變歷史的靜力解答。所以，可以認(rèn)為在應(yīng)用上，本文方法不受有限單元類型和材料本構(gòu)模型特征的限制，可進(jìn)一步研究此方法在由實(shí)驗(yàn)標(biāo)定的、更能真實(shí)反映材料復(fù)雜受力情況的本構(gòu)模型的常見結(jié)構(gòu)靜力分析中的應(yīng)用。

3) 在時(shí)間域進(jìn)行結(jié)構(gòu)受力問題分析時(shí)，問題具有適定性[36]?？蓢L試應(yīng)用DRM、NEIA對(duì)依賴應(yīng)變歷史的、應(yīng)變軟化的、變形局部化的鋼筋混凝土結(jié)構(gòu)進(jìn)行靜力分析。

4) 若僅從求解靜力增量解的迭代次數(shù)衡量，DRM算法遠(yuǎn)大于ISA算法。但DRM的優(yōu)越性在于：其一，它規(guī)避了ISA求解靜力解時(shí)的總體剛度矩陣K的組裝和求逆，此運(yùn)算涉及較多計(jì)算成本，且此計(jì)算成本與計(jì)算模型尺寸成指數(shù)比例關(guān)系；其二，它僅涉及向量運(yùn)算，計(jì)算成本與計(jì)算模型尺寸成線性比例關(guān)系。此外，也可構(gòu)造一些降低DRM迭代次數(shù)的技巧，如線性外推策略。

5 結(jié)論

本文提出采用非線性彈性增量算法(NEIA)，更新SHDM結(jié)構(gòu)(依賴應(yīng)變歷史材料制成的結(jié)構(gòu))非線性有限元靜力增量步動(dòng)力松弛法(DRM)迭代求解時(shí)的高斯點(diǎn)應(yīng)力，以此來(lái)規(guī)避迭代解中的虛假應(yīng)變歷史。NEIA不同于非線性增量算法(NIA)：使用NEIA時(shí)，一個(gè)增量步內(nèi)高斯點(diǎn)的應(yīng)變加卸載路徑被施加該增量步增量前的高斯點(diǎn)本構(gòu)狀態(tài)所固定；而使用NIA時(shí)，高斯點(diǎn)的應(yīng)變加卸載路徑隨增量步內(nèi)虛擬時(shí)間而變化。本文從總體上的增量步解、增量步中的迭代解、迭代步中的應(yīng)力更新步驟三個(gè)層次闡述了結(jié)構(gòu)非線性有限元靜力問題的DRM解法。最后，分別采用NEIA、NIA進(jìn)行了三個(gè)算例的DRM求解。據(jù)此，本文的結(jié)論如下：

(1) 采用基于虛擬動(dòng)力過(guò)程構(gòu)造的DRM求解增量步時(shí)，原靜力問題適定可解，且此問題的解等于上述虛擬動(dòng)力過(guò)程的穩(wěn)態(tài)解。采用NIA時(shí)，此穩(wěn)態(tài)解為包含虛假應(yīng)變歷史效應(yīng)的近似解；而采用NEIA時(shí)，此穩(wěn)態(tài)解為消除虛假應(yīng)變歷史效應(yīng)意義上的正確解，從而NEIA能給出具有真實(shí)應(yīng)變歷史的結(jié)構(gòu)靜力平衡路徑。

(2) 在保持本構(gòu)關(guān)系正確的前提下，NEIA給出的平衡路徑不依賴于增量大??；而NIA給出的平衡路徑依賴于增量大小，當(dāng)增量趨小時(shí)，平衡路徑趨于NEIA給出的平衡路徑。

(3) 增量步的DRM解法通過(guò)動(dòng)態(tài)阻尼、收斂發(fā)散判據(jù)、荷載位移增量轉(zhuǎn)換規(guī)則進(jìn)行靜力平衡路徑硬化段、極值點(diǎn)、軟化段的全過(guò)程分析。

(4) 增量步的DRM解法無(wú)需隱式靜力算法的整體剛度矩陣K的組裝、求逆。

本文方法解決了SHDM結(jié)構(gòu)非線性有限元全過(guò)程靜力計(jì)算中的虛假應(yīng)變歷史問題。