蘇小卒，王 偉

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

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

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

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

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

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

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

1 問題的描述

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

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

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

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

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

2 問題的求解

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

2.1 靜力問題的增量解法

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

式中：Δan為結點位移增量；Δεn為應變增量；Δσn為應力增量；Δλn為外部結點荷載系數增量。

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

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

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

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

4) 記i為kDRM求解的迭代步序號，對于求得的

① 若滿足收斂判據

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

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

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

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

2.2 增量步的kDRM解法

2.2.1 kDRM的控制方程

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

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

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

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

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)的函數。

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

式中：e表示單元編號；ne表示系統(tǒng)內的單元數量；p{e}、B{e}、σ{e}、Ω{e}分別為單元e的內部結點荷載、應變矩陣、應力向量、所占空間；p{e}、B{e}、σ{e}分別為的矩陣(為單元e的結點自由度數)；T{e}為單元e的轉換矩陣，尺寸為(r為系統(tǒng)全部自由度數)。因為C(t)是nan(t)的函數，所以pn(t)也為an(t)的函數。

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

② 加位移時為未知的變量

式中：fref,I、pn,I為fref、pn(t)中與第I個自由度(加位移自由度)對應的元素值；λn(t)是an(t)的函數。

由本小節(jié)分析可知：kDRM求解由問題9所描述的Rn的關鍵是，給出基本未知量an(t)的數學算法。因為問題9是一個強非線性動力問題，所以求解an(t)時，宜選用顯式時間積分算法[29]。

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

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

式中：ik-1為Sk的起始時刻序號；h為離散時間域的固定時間增量；的計算式為：

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

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

式中：下標表示結點序號；mj為結點質量；kjmax為結點沿各坐標正負方向位移剛度最大值。本文估算kjmax的方式為：首先給結點施加微小位移，然后計算相應結點荷載，最后算出kjmax。本文取

式中：c3為系數，值大于1；mmin是結點質量下限值，其作用是避免結點負剛度而構造出負的結點質量。本文按式(16)逐結點算得mj后構造Mn。

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

則認為(i-1)h時離散質點系統(tǒng)出現了動能極大值，Sk的結束時間以及Sk+1的起始時間取為(i-1)h，以為初始條件啟動Sk+1。

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

2.2.3 完整的kDRM數值計算流程

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

2.3 高斯點應力更新：非線性彈性增量算法

依賴應變歷史材料(SHDM)在應變加載、卸載時具有不同材料響應，故對結構有限元模型進行受力分析時，s需記錄高斯點加載、卸載時應力-應變關系的狀態(tài)。

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

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

式中，fσ為材料點應力計算函數。按此算法，僅當滿足收斂判據時才更新C。NEIA不同于結構動力分析時(包括EQSA)的NIA——在當前增量步內，ih時刻的由ih時刻的及(i-1)h時刻的更新(注：更新)，即：

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

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

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

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

3 算例

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

3.1 算例I

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

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

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

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

表2 模型I 桿的本構關系參數及截面面積Table 2 The model I: parameters of constitutive model and cross-section area of bar

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

3.2 算例II

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

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

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

式中：dσ為應力增量；dε為應變增量；D為材料剛度矩陣。D為：

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

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

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

3.3 算例III

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

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

模型III使用Mazars提出的單標量彈性損傷本構模型[32－33]：

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

式中：ωmax是材料點所達到過的最大損傷值；α+ω++α-ω-是此材料點的即時損傷值，其中α+、α-是依賴于ε的系數，ω-、ω+是受壓、受拉損傷變量。α+、α-的算式見文獻[33]。ω+、ω-的算式為：

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

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

3.4 計算結果及分析

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

2) 從圖6、圖8、圖10可知：① 采用NEIA時，不同的增量幅度(IA、IB、IC，IIA、IIB、IIC，IIIA、IIIB、IIIC)給出相同的平衡路徑；② 采用NIA時，不同的增量幅度(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) 對于模型I，由圖7可知(結合表2)：4號桿、5號桿的材料性能階段與相應的平衡路徑階段(圖6)相關，如結構的軟化階段始于5號單元的應變軟化，且同時伴隨著4號桿的應變卸載，此現象即為形變局部化。

4) 對于模型II，對于x=121 mm處的梁橫截面S，假定S應變保持平面，則由材料力學公式得：① 當S僅有最外邊緣纖維應力等于fy時(狀態(tài)e)，f1y=59.7 kN；② 當S的所有纖維應力等于fy時(狀態(tài)p)，f1y=89.6 kN。依據圖8、圖9，并結合圖4給出的單元1的1號高斯點、單元2的2號高斯點的y坐標，得：① 當a1y≈24 mm時，S剛越過狀態(tài)e，此時結構剛度開始降低，f1y≈68.0 kN，此與理論值接近；② 當a1y=49 mm時，S臨近狀態(tài)p，此時結構剛度臨近零，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) 當加載結點的增量幅度超過一定范圍時，NEIA可能會給出不真實的平衡路徑，如算例I的ID的平衡路徑不同于IA、IB、IC的平衡路徑。從計算所得的應變數據可知：ID的5號桿跨越了應變屈服階段。此導致f1y-a1y曲線不存在屈服階段、算得的平衡路徑偏離真實平衡路徑。由此可知，即使采用NEIA求解平衡路徑，也應避免采用過大的增量幅度。

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

3) 在時間域進行結構受力問題分析時，問題具有適定性[36]?？蓢L試應用DRM、NEIA對依賴應變歷史的、應變軟化的、變形局部化的鋼筋混凝土結構進行靜力分析。

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

5 結論

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

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

(2) 在保持本構關系正確的前提下，NEIA給出的平衡路徑不依賴于增量大?。欢鳱IA給出的平衡路徑依賴于增量大小，當增量趨小時，平衡路徑趨于NEIA給出的平衡路徑。

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

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

本文方法解決了SHDM結構非線性有限元全過程靜力計算中的虛假應變歷史問題。