基于MCFT理論的割線剛度法介紹

1、基于MCFT理論的割線剛度法介紹F. J. Vecchio于1990年在ACI發(fā)表的論文中采用了基于MCFT理論的割線剛度法來(lái)進(jìn)行鋼筋混凝土的有限元分析。通過(guò)比較割線剛度法與基于切線剛度的Newton-Raphson法，可以得出以下幾點(diǎn)結(jié)論：1. 在結(jié)構(gòu)的應(yīng)力-應(yīng)變曲線進(jìn)入到下降段后，對(duì)于切線剛度法，剛度會(huì)出現(xiàn)負(fù)值，從而給計(jì)算處理帶來(lái)難度。而割線剛度法則沒(méi)有這種負(fù)剛度問(wèn)題，但是割線剛度法會(huì)在卸載曲線穿越坐標(biāo)軸時(shí)出現(xiàn)數(shù)值問(wèn)題，這是由于此時(shí)割線剛度極小造成的。2. 以往采用切線剛度法的程序，在更新應(yīng)變時(shí)，對(duì)于一個(gè)單元采用了4個(gè)高斯積分點(diǎn)，每次迭代得到了一個(gè)新的結(jié)構(gòu)位移，就需要在4個(gè)積分點(diǎn)上分別求出

2、它們的應(yīng)變，再以此計(jì)算出4個(gè)新的割線剛度，用以判斷收斂與否。而F. J. Vecchio采用的割線剛度法只涉及了1個(gè)積分點(diǎn)，即中心積分點(diǎn)。根據(jù)附表1可知，單元應(yīng)變和位移的關(guān)系為：其中：令，則可以得到割線剛度法中心積分點(diǎn)的應(yīng)變計(jì)算式為： 每次循環(huán)迭代出的新位移，僅需要換算成中心積分點(diǎn)上的應(yīng)變，通過(guò)一次計(jì)算便可以得出割線剛度用于收斂判斷。相比于切線剛度法，割線剛度法簡(jiǎn)化了計(jì)算過(guò)程。3. 切線剛度法的過(guò)程，需要通過(guò)外荷載節(jié)點(diǎn)力和結(jié)構(gòu)恢復(fù)力求解出一個(gè)殘余力，再以此計(jì)算出增量位移來(lái)判斷是否滿足收斂條件，若不收斂，便更新結(jié)構(gòu)位移并重新計(jì)算結(jié)構(gòu)的切線剛度矩陣。對(duì)于割線剛度法，并不需要求解殘余力，而是直接利用

3、外荷載節(jié)點(diǎn)力和結(jié)構(gòu)割線剛度來(lái)求解出一個(gè)全量位移，之后通過(guò)計(jì)算出的應(yīng)變來(lái)重新計(jì)算結(jié)構(gòu)的割線剛度矩陣。二者相比起來(lái)，割線剛度法不需要啟用一個(gè)新的數(shù)組來(lái)儲(chǔ)存殘余力值，降低了對(duì)程序的存儲(chǔ)空間要求并提高了程序的運(yùn)算速度。從上面三點(diǎn)可以看出，切線剛度法和割線剛度法各有特點(diǎn)和不足，在某些方面，割線剛度法還體現(xiàn)出了一定的優(yōu)勢(shì)：割線剛度法相比于切線剛度法的收斂準(zhǔn)則更為簡(jiǎn)單。割線剛度法相比于切線剛度法的有限元模型更為簡(jiǎn)單。割線剛度法的穩(wěn)定性較好。所以對(duì)割線剛度法進(jìn)行研究佐證是十分有必要的。F. J. Vecchio在其論文中所述的割線剛度法的具體過(guò)程如下：第1步：輸入結(jié)構(gòu)和材料的屬性數(shù)據(jù)第2步：計(jì)算外荷載節(jié)點(diǎn)力第

4、3步：計(jì)算材料的剛度矩陣：1) 確定初始應(yīng)變：，； 2) 計(jì)算主方向應(yīng)變，以及主拉應(yīng)變角；3) 根據(jù)混凝土的本構(gòu)關(guān)系，確定混凝土主拉和主壓應(yīng)力，4) 根據(jù)鋼筋的本構(gòu)關(guān)系，確定鋼筋的主拉和主壓應(yīng)力，；5) 計(jì)算混凝土和鋼筋的割線剛度元素值： (3-1) 第4步：計(jì)算混凝土和鋼筋的材料割線剛度矩陣1) 計(jì)算混凝土主方向上的材料割線剛度矩陣： (3-2)2) 計(jì)算轉(zhuǎn)換矩陣3) 計(jì)算整體坐標(biāo)系（X-Y坐標(biāo)）下的混凝土割線剛度矩陣 (3-3)4) 計(jì)算整體坐標(biāo)系(X-Y坐標(biāo)系)下的鋼筋割線剛度矩陣： (3-4)第5步：計(jì)算組合材料剛度： (3-5)第6步：計(jì)算單元?jiǎng)偠染仃嚕?(3-6)第7步：計(jì)算單元預(yù)

5、應(yīng)變產(chǎn)生的位移，第8步：計(jì)算單元預(yù)應(yīng)變產(chǎn)生的節(jié)點(diǎn)力，第9步：計(jì)算單元節(jié)點(diǎn)力：第10步：計(jì)算結(jié)構(gòu)總剛矩陣： (3-7) 圖 3.5 單元尺寸第11步：計(jì)算結(jié)構(gòu)總剛矩陣的逆矩陣：第12步：計(jì)算節(jié)點(diǎn)位移 Fig3.5 Element dimension (3-8) 第13步：計(jì)算單元應(yīng)變（取單元形心處的應(yīng)變?yōu)槠骄鶓?yīng)變）： (3-9) (3-10) (3-11)注：，為相應(yīng)節(jié)點(diǎn)的水平位移和豎向位移。第14步：檢查單元應(yīng)力：第15步：重新計(jì)算單元割線剛度值：1) 根據(jù)求解出的位移，重新確定主方向上的應(yīng)變值，；3) 根據(jù)混凝土本構(gòu)模型，計(jì)算新的混凝土應(yīng)力，；4) 根據(jù)鋼筋本構(gòu)模型，計(jì)算新的鋼筋應(yīng)力，；5)

6、 重新計(jì)算鋼筋和混凝土的割線剛度元素：，,。第16步 收斂判斷詳細(xì)的算法流程圖如下所示：2.1 修正斜壓場(chǎng)理論（MCFT）介紹修正斜壓場(chǎng)理論（The Modified Compression-Field Theory）是一種用于預(yù)測(cè)鋼筋混凝土膜構(gòu)件在剪力和正應(yīng)力作用下的變形性能的分析模型，它可以確定混凝土和鋼筋的平均應(yīng)力應(yīng)變和局部的應(yīng)力應(yīng)變，以及裂縫的寬度和方向，從而確定構(gòu)件的破壞模式。MCFT是由研究鋼筋混凝土受扭和受剪的斜壓場(chǎng)理論（Compression-Field Theory）發(fā)展而來(lái)。斜壓場(chǎng)理論忽略了開裂混凝土的拉力，而修正斜壓場(chǎng)理論則把裂縫間的拉應(yīng)力考慮了進(jìn)來(lái)。MCFT把開裂混凝土

7、模擬為一種正交各向異性的材料，是一種完全的轉(zhuǎn)角裂縫模型（Smeared Rotating Crack Model）。與固體的不連續(xù)性相悖，開裂混凝土被視作一種裂縫均勻分布的固體連續(xù)體。裂縫可以自由轉(zhuǎn)動(dòng)，并保持和混凝土的主壓應(yīng)力場(chǎng)方向同向。修正斜壓場(chǎng)理論包含了三種關(guān)系：用平均應(yīng)變表示的混凝土和鋼筋的變形協(xié)調(diào)關(guān)系；用平均應(yīng)力表示的混凝土和鋼筋的應(yīng)力平衡關(guān)系；用平均應(yīng)力和平均應(yīng)變表示的混凝土和鋼筋的本構(gòu)關(guān)系。開裂混凝土的本構(gòu)關(guān)系是根據(jù)多倫多大學(xué)所做的鋼筋混凝土剪切板實(shí)驗(yàn)的結(jié)果得到的，所以MCFT采用的是實(shí)驗(yàn)觀測(cè)到的真實(shí)的本構(gòu)模型。同時(shí)，由于裂縫均勻分布以及上述三種關(guān)系用平均應(yīng)力和平均應(yīng)變的概念表述，

8、所以MCFT的一個(gè)重點(diǎn)是要考慮裂縫處的局部應(yīng)力應(yīng)變條件。圖2.1 開裂鋼筋混凝土的修正壓場(chǎng)理論23Fig. 2.1 Modified Compression Field Theory for cracked reinforced concrete23修正斜壓場(chǎng)理論作了以下一些假設(shè)：鋼筋彌散分布于單元之中；裂縫均勻分布并能自由旋轉(zhuǎn)；單元中，剪應(yīng)力和正應(yīng)力均勻分布；應(yīng)力狀態(tài)和應(yīng)變狀態(tài)一一對(duì)應(yīng)，不考慮加載歷史的影響；跨越多條裂縫間距，應(yīng)力應(yīng)變可視作是平均的；主應(yīng)力角和主應(yīng)變角的方向一致；鋼筋和混凝土之間無(wú)粘結(jié)滑移；混凝土和鋼筋有獨(dú)立的本構(gòu)關(guān)系；忽略鋼筋中的剪應(yīng)力。2.1.1 變形協(xié)調(diào)條件在混凝土和鋼

9、筋構(gòu)件中，與平均應(yīng)變相關(guān)的變形協(xié)調(diào)關(guān)系如圖 所示： 圖2.2 混凝土平均應(yīng)變28Fig.2.2 Average concrete strains28因?yàn)榧僭O(shè)鋼筋和混凝土之間沒(méi)有粘結(jié)滑移，所以對(duì)于非預(yù)應(yīng)力鋼筋而言，鋼筋和混凝土中的平均應(yīng)變是相等的。變形協(xié)調(diào)條件可以用一下公式表達(dá)：若剪應(yīng)變值為，則根據(jù)莫爾應(yīng)變圓確定的混凝土平均主拉應(yīng)變和平均主壓應(yīng)變?yōu)椋和ㄟ^(guò)莫爾應(yīng)變圓也可以確定平均主拉應(yīng)變和應(yīng)力與軸的夾角：2.1.2 應(yīng)力平衡關(guān)系從構(gòu)件中取出一個(gè)自由體作為研究對(duì)象，如圖 所示：和方向的力平衡條件要求混凝土平均應(yīng)力和以及鋼筋平均應(yīng)力和的合力要與施加的正應(yīng)力和的合力相平衡。彎矩平衡條件則要求混凝土中的平

10、均剪應(yīng)力要能完全抵抗所施加的剪應(yīng)力（假設(shè)鋼筋沒(méi)有銷栓作用）。平均應(yīng)力的平衡關(guān)系可以表述如下：和分別是和方向的配筋率。因?yàn)殚_裂混凝土是主應(yīng)力方向上的正交各向異性材料，莫爾應(yīng)力圓可以用于建立混凝土平均應(yīng)力和與混凝土平均主拉應(yīng)力之間的關(guān)系式：2.1.3 本構(gòu)關(guān)系本構(gòu)模型需要把變形協(xié)調(diào)條件中的應(yīng)變和平衡條件中的應(yīng)力聯(lián)系起來(lái)。1986年，多倫多大學(xué)的Vecchio和Collins做了30個(gè)板實(shí)驗(yàn)，這些板的尺寸都是，并承受平面內(nèi)的應(yīng)力。通過(guò)對(duì)實(shí)驗(yàn)結(jié)果的分析，他們得出了開裂混凝土受拉和受壓的本構(gòu)模型。對(duì)于受壓混凝土，本構(gòu)關(guān)系是主壓應(yīng)力和主壓應(yīng)變的函數(shù)。上述板的實(shí)驗(yàn)結(jié)果表明，在主拉應(yīng)變同時(shí)存在的情況下，抗壓強(qiáng)

11、度和剛度會(huì)隨著的增大而減小。這一現(xiàn)象被稱為受壓軟化，具體表現(xiàn)在混凝土單軸受壓時(shí)的應(yīng)力-應(yīng)變反應(yīng)的弱化。建議的考慮受壓軟化的公式如下：其中，是單軸受壓混凝土的Hognestad拋物線函數(shù)，常用于一般強(qiáng)度的混凝土。是對(duì)應(yīng)于峰值抗壓應(yīng)力的混凝土圓柱體應(yīng)變（取負(fù)值），由混凝土圓柱體單軸受壓實(shí)驗(yàn)確定。反映的是主拉應(yīng)變的軟化效應(yīng)。對(duì)于受拉混凝土，本構(gòu)關(guān)系是主拉應(yīng)力和主拉應(yīng)變的函數(shù)。首先，必須要確定開裂強(qiáng)度以及相應(yīng)的開裂應(yīng)變。如果信息不足，可以作如下估計(jì)： 其中是混凝土的初始剛度，計(jì)算式為： 在開裂以前，混凝土在受拉狀態(tài)下是線彈性的： 開裂以后，由于混凝土和鋼筋的粘結(jié)作用，拉應(yīng)力將繼續(xù)存在于鋼筋混凝土裂縫間

12、的混凝土中。為了模擬這種受拉硬化現(xiàn)象，混凝土的拉應(yīng)力要隨著混凝土主拉應(yīng)變的增大從抗拉強(qiáng)度逐步減小。MCFT建議的關(guān)系公式如下：對(duì)于處在受拉和受壓狀態(tài)的鋼筋而言，MCFT采用了平均應(yīng)力和平均應(yīng)變的兩折線模型：一個(gè)線彈性的上升段和一個(gè)屈服平臺(tái)。公式如下：其中，是鋼筋的彈性模量，和分別是鋼筋在和方向的屈服應(yīng)力。2.1.4 考慮局部平衡條件有了相容的應(yīng)變條件，根據(jù)本構(gòu)關(guān)系可以確定混凝土和鋼筋中的平均應(yīng)力，以及它們所平衡的施加剪應(yīng)力和正應(yīng)力。然而，如果忽略了構(gòu)件性能由裂縫處的局部鋼筋屈服或沿裂縫的剪切滑移破壞所控制的可能性，那么分析結(jié)果將會(huì)是不保守的。為了把這些可能性考慮進(jìn)來(lái)，MCFT限制了裂縫處的局部

13、應(yīng)力以及混凝土的平均拉應(yīng)力水平。鋼筋混凝土的應(yīng)力場(chǎng)從裂縫間的平均概念到裂縫處的局部概念是不斷變化的。圖 a所示的是裂縫間垂直于主拉應(yīng)力方向的截面上的平均應(yīng)力情況，而圖 b所示的是裂縫處截面上的局部應(yīng)力情況。在裂縫處的截面上，混凝土平均拉應(yīng)力減小到0。為了跨過(guò)裂縫傳遞平均拉應(yīng)力，裂縫處的鋼筋應(yīng)力和應(yīng)變就必須局部提高。在裂縫表面的法線方向根據(jù)力平衡條件，用平均應(yīng)力和局部應(yīng)力建立的數(shù)值等價(jià)關(guān)系式可以表述為：其中，和是裂縫處的局部鋼筋應(yīng)力，和是鋼筋和裂縫截面法線方向的夾角。通過(guò)上面的公式可以明顯看出，混凝土的平均拉應(yīng)力受到了裂縫處鋼筋屈服的限制。如果用鋼筋的屈服強(qiáng)度取代鋼筋局部應(yīng)力，那么公式兩個(gè)括號(hào)中

14、的項(xiàng)就表示了鋼筋的富余強(qiáng)度，也就限制了后開裂混凝土的受拉應(yīng)力不超過(guò)：如圖 a所示，主平面上沒(méi)有剪應(yīng)力。但如果鋼筋與裂縫斜交，那么就會(huì)在裂縫表面產(chǎn)生局部剪應(yīng)力。在裂縫切線方向上根據(jù)力的平衡關(guān)系，可以得到如下公式：不考慮上述方程，局部剪應(yīng)力可以一直增大，直到發(fā)生剪切滑移破壞。剪應(yīng)力的大小會(huì)受到骨料咬合作用的限制，這種骨料咬合作用隨著裂縫寬度的增大和最大骨料尺寸的減小而減小?；赪alraven對(duì)骨料咬合作用的分析，MCFT對(duì)裂縫處剪應(yīng)力大小的限制為：平均裂縫寬度是混凝土主拉應(yīng)變和平均裂縫間距的函數(shù)：和方向的平均裂縫間距和可以根據(jù)粘結(jié)屬性和鋼筋分布情況求出。如果超過(guò)了允許的混凝土最大拉應(yīng)力或者裂縫處

15、的局部剪應(yīng)力，就需要修正應(yīng)變狀態(tài)以求得到較小的混凝土平均拉應(yīng)力值。2.2 混凝土滯回本構(gòu)模型圖 a所示為混凝土受壓的滯回模型，其骨架曲線表示的是結(jié)構(gòu)在單調(diào)加載時(shí)的反應(yīng)。根據(jù)MCFT可知，混凝土單調(diào)本構(gòu)曲線是以Hognestad拋物線為基礎(chǔ)，通過(guò)考慮受壓軟化效應(yīng)，對(duì)其進(jìn)行修正得到的。設(shè)混凝土塑性應(yīng)變?yōu)?，?dāng)再加載時(shí)，混凝土的壓應(yīng)力可以通過(guò)下式計(jì)算：其中，是之前加載得到的最大壓應(yīng)變，是與相應(yīng)的應(yīng)力，是通過(guò)骨架曲線計(jì)算得到的與相應(yīng)的應(yīng)力。如果結(jié)構(gòu)的反應(yīng)落在骨架曲線上，如，那么就將和更新為和。在每一個(gè)荷載步，即時(shí)的塑性應(yīng)變用下式計(jì)算：其中，是與骨架曲線上峰值應(yīng)力相對(duì)應(yīng)的應(yīng)變，如果即時(shí)塑性應(yīng)變大于塑性補(bǔ)償

16、應(yīng)變，那么就將后者更新。任意時(shí)刻的卸載產(chǎn)生的應(yīng)力用如下關(guān)系計(jì)算：其中卸載模量定義為：這一模型隱含的假設(shè)是，當(dāng)混凝土由拉變?yōu)閴簳r(shí)，壓應(yīng)力保持為0，直到裂縫完全閉合，即直到。實(shí)際上，實(shí)驗(yàn)證明這個(gè)假設(shè)是不正確的。由于之后的應(yīng)變總是會(huì)大于0，且受到裂縫剪切滑移的影響，此時(shí)壓應(yīng)力不會(huì)等于0?；炷潦芾臏乇緲?gòu)如圖 b所示。骨架曲線由兩部分組成：開裂前的曲線和開裂后的受拉硬化曲線。設(shè)塑性應(yīng)變?yōu)?，再加載時(shí)的混凝土拉應(yīng)力用下式計(jì)算：其中，是之前加載得到的最大拉應(yīng)變，是與相對(duì)應(yīng)的應(yīng)力，是通過(guò)骨架曲線計(jì)算得到的與相對(duì)應(yīng)的應(yīng)力。對(duì)于第一個(gè)受拉加卸載，塑性補(bǔ)償應(yīng)變?yōu)槌?shù)，其值等于從受壓區(qū)穿軸確定的塑性應(yīng)變大小。接下

17、來(lái)用于骨架曲線計(jì)算的是凈塑性應(yīng)變。在之后的循環(huán)中，塑性補(bǔ)償應(yīng)變用下式計(jì)算：由于目前缺乏合適的模型，從上式可以看出，沒(méi)有正值的補(bǔ)償應(yīng)變。這就使得在加載和卸載的過(guò)程中，受拉滯回曲線都會(huì)經(jīng)過(guò)原點(diǎn)。當(dāng)混凝土反應(yīng)落在骨架曲線上時(shí)，最大拉應(yīng)變和相應(yīng)的應(yīng)力都分別更新到和。卸載時(shí)的應(yīng)力可以通過(guò)下式計(jì)算：其中，卸載模量定義為：2.3 鋼筋的滯回本構(gòu)模型本文的鋼筋的應(yīng)力-應(yīng)變關(guān)系是采用的經(jīng)Filippou等人改進(jìn)后考慮了各向同性應(yīng)變強(qiáng)化的Menegotto-Pinto非線性模型。Menegotto-Pinto模型的表達(dá)式為： （2.64）其中， （2.65） （2.66）方程(2.64)表達(dá)的是斜率分別為和的兩條

18、漸近線的過(guò)渡段，如圖 所示的線（a）和線（b）。點(diǎn)B是兩條漸近線的交點(diǎn)，其應(yīng)力和應(yīng)變分別為和。 A點(diǎn)是荷載的反向點(diǎn)，其應(yīng)力和應(yīng)變分別為和。b為和的比值，叫做應(yīng)變強(qiáng)化率。是影響過(guò)渡曲線形狀的一個(gè)參數(shù)，它反映了鋼筋的包興格（Bauschinger）效應(yīng)。圖2.13 Menegotto-Pinto 鋼筋模型Fig 2.13 Menegotto-Pinto steel model的值取決于當(dāng)前漸進(jìn)線交點(diǎn)B的應(yīng)變和上一次荷載的反向點(diǎn)A的最大或最小應(yīng)變間的應(yīng)變差。Menegotto-Pinto模型中的表達(dá)式為： （2.67）其中， 是第一次加載時(shí)參數(shù)的初始值，和跟一樣，其值由鋼筋試驗(yàn)所確定，將在每次應(yīng)變反

19、向后更新其值。圖2.14 Menegotto-Pinto 鋼筋模型中曲率參數(shù)的定義34Fig 2.14 Definition of curvature parameter in Menegotto-Pinto steel model34上述Menegotto-Pinto模型主要缺點(diǎn)是，它無(wú)法考慮等向應(yīng)變硬化的影響。針對(duì)這個(gè)問(wèn)題，F(xiàn)illipou (1983)建議把線性屈服漸近線上的應(yīng)力變換看作是最大塑性應(yīng)變的函數(shù)： （2.68）其中，和分別是鋼筋屈服時(shí)的應(yīng)變和應(yīng)力，是當(dāng)應(yīng)變反向時(shí)絕對(duì)值最大的應(yīng)變， 與的值由試驗(yàn)確定?？紤]了等向應(yīng)變硬化影響修正后的鋼筋滯回曲線如圖 所示。圖2.15 鋼筋滯回本構(gòu)

20、模型34Fig 2.15 Hysteresis Model for reinforcrment342.4 基于切線剛度的Newton-Raphson迭代法鄧興龍和蔣曉華等編制的非線性有限元程序是采用的基于切線剛度的Newton-Raphson迭代法理論得以實(shí)現(xiàn)的。下面對(duì)其所涉及的高斯積分以及非線性方程組解法中的Newton-Raphson迭代算法、弧長(zhǎng)法和收斂準(zhǔn)則做一些簡(jiǎn)單的介紹。2.4.1 高斯積分基于切線剛度的Newton-Raphson迭代法采用的是積分點(diǎn)的高斯積分，其單元?jiǎng)偠染仃嚭蛦卧刃ЫY(jié)點(diǎn)恢復(fù)力的高斯積分公式介紹如下。至于具體的高斯求積公式的推導(dǎo)過(guò)程，可以參見文獻(xiàn)。用高斯求積公式表

21、述的單元?jiǎng)偠染仃囉?jì)算式為：用高斯求積公式表述的單元等效結(jié)點(diǎn)恢復(fù)力的計(jì)算式為：其中，為應(yīng)變矩陣，為雅克比行列式。和的值與高斯積分點(diǎn)相對(duì)應(yīng)。針對(duì)不同的單元，一般通過(guò)經(jīng)驗(yàn)和試算來(lái)確定最佳的積分點(diǎn)數(shù)。對(duì)于四邊形單元而言，的高斯積分點(diǎn)就可以較好地滿足精度要求。表 中所示的是各個(gè)高斯積分點(diǎn)的坐標(biāo)和加權(quán)系數(shù)，圖 所示的是高斯積分點(diǎn)的分布情況。表2.3 高斯積分中的和 Table 2.3 and in Gauss Integration點(diǎn)數(shù)n積分點(diǎn)坐標(biāo)加權(quán)系數(shù)點(diǎn)數(shù)n積分點(diǎn)坐標(biāo)加權(quán)系數(shù)1024-0.8611363116-0.33998104360.34785484510.65214515492-0.577350

22、26921130.774596669200.5555555560.8888888895-0.9061798459-0.538469310100.23692688510.47862867050.5688888889圖2.9 高斯點(diǎn)的位置Fig 2.9 The location of Gauss Points2.4.2 Newton-Raphson迭代法基于切線剛度的Newton-Rapson迭代法如圖3.6。0000000000P圖3.6 基于切線剛度的Newton-Raphson算法Fig3.6 Newton-Raphson algorithm based on tangent stiffne

23、ssNewton-Raphson是一種變剛度的迭代法。首先，取單元的初始剛度為，根據(jù)所施加的外荷載等效結(jié)點(diǎn)力可以求出位移的第一次近似值為： (3.1)由得到的初始位移可以計(jì)算單元的應(yīng)變值，再結(jié)合本構(gòu)模型求到單元應(yīng)力，進(jìn)而得到相應(yīng)的結(jié)點(diǎn)恢復(fù)力。然后，用外荷載等效結(jié)點(diǎn)力和結(jié)構(gòu)恢復(fù)力求出殘余力，并連同相應(yīng)于初始位移的即時(shí)切線剛度，求出位移增量，即 (3.2) (3.3)得到位移的第二次近似值為： (3.4)以此類推，第k步的計(jì)算過(guò)程為： (3.5) (3.6) (3.7)當(dāng)與的誤差達(dá)到所規(guī)定的限值，或者當(dāng)足夠小的時(shí)候，迭代結(jié)束。 2.4.3 弧長(zhǎng)控制法弧長(zhǎng)控制法的基本原理是，通過(guò)在結(jié)構(gòu)分析的一般方程

24、中添入一個(gè)新的變量，來(lái)實(shí)現(xiàn)對(duì)整個(gè)計(jì)算過(guò)程的自動(dòng)控制。加入變量后，原方程組將變?yōu)椋?（2.99）上述方程組一共有個(gè)未知量：、，但是方程卻只有個(gè)。所以，需要補(bǔ)充如下條件來(lái)完成方程組的求解： （2.100）稱為“弧長(zhǎng)”，是附加的控制變量。這時(shí)，方程組變?yōu)椤；￠L(zhǎng)法的迭代過(guò)程如圖 所示。圖2.14 弧長(zhǎng)法迭代過(guò)程示意圖Fig 2.14 The Arc-length method2.4.4 收斂準(zhǔn)則基于切線剛度的Newton-Raphson迭代法的收斂準(zhǔn)則有三種：力準(zhǔn)則、位移準(zhǔn)則以及能量準(zhǔn)則。選取合適的收斂準(zhǔn)則，有利于提高迭代法的計(jì)算效率，是增量求解策略中一個(gè)重要的組成部分。當(dāng)?shù)Y(jié)束時(shí)，結(jié)構(gòu)應(yīng)當(dāng)處于一個(gè)

25、穩(wěn)定的力平衡狀態(tài)，力準(zhǔn)則就是要求結(jié)構(gòu)的恢復(fù)力和外力之差等于零。但是在數(shù)值計(jì)算中，不可能實(shí)現(xiàn)不平衡力完全為零，這就需要設(shè)定一個(gè)允許值。力準(zhǔn)則的表達(dá)式為：其中，就是給殘余力設(shè)置的一個(gè)允許值。在進(jìn)行位移分析時(shí)，計(jì)算位移應(yīng)該與真實(shí)值接近。位移準(zhǔn)則的表達(dá)式為：其中，是矢量的歐幾里德范數(shù)，是第步迭代的位移增量，為設(shè)置的允許值。能量準(zhǔn)則是用來(lái)度量力和位移同它們的平衡值之間的接近程度的，其表達(dá)式為： （2.115)其中，不等式的左邊表示殘余力在位移增量上所做的功，不等式的右邊是殘余力所做的功的初始值，是為能量設(shè)置的允許值。鄧興龍和蔣曉華等編制的FEAPpv程序采用的是能量收斂準(zhǔn)則。2.5 基于MCFT理論的割

26、線剛度迭代法F. J. Vecchio于1990年在ACI發(fā)表的論文中采用了基于MCFT理論的割線剛度迭代法來(lái)進(jìn)行鋼筋混凝土的有限元分析。通過(guò)比較割線剛度迭代法與基于切線剛度的Newton-Raphson法，可以得出以下幾點(diǎn)結(jié)論：1. 在結(jié)構(gòu)的應(yīng)力-應(yīng)變曲線進(jìn)入到下降段后，對(duì)于切線剛度法，剛度會(huì)出現(xiàn)負(fù)值，從而給計(jì)算處理帶來(lái)難度。而割線剛度迭代法則沒(méi)有這種負(fù)剛度問(wèn)題，但是割線剛度迭代法會(huì)在卸載曲線穿越坐標(biāo)軸時(shí)出現(xiàn)數(shù)值問(wèn)題，這是由于此時(shí)割線剛度極小造成的。2. 以往采用切線剛度法的程序，在更新應(yīng)變時(shí)，對(duì)于一個(gè)單元采用了4個(gè)高斯積分點(diǎn)，每次迭代得到了一個(gè)新的結(jié)構(gòu)位移，就需要在4個(gè)積分點(diǎn)上分別求出它們

27、的應(yīng)變，再以此計(jì)算出4個(gè)新的割線剛度，用以判斷收斂與否。而F. J. Vecchio采用的割線剛度迭代法只涉及了1個(gè)積分點(diǎn)，即中心積分點(diǎn)。根據(jù)附錄1可知，單元應(yīng)變和位移的關(guān)系為：其中：令，則可以得到割線剛度迭代法中心積分點(diǎn)的應(yīng)變計(jì)算式為： 每次循環(huán)迭代出的新位移，僅需要換算成中心積分點(diǎn)上的應(yīng)變，通過(guò)一次計(jì)算便可以得出割線剛度用于收斂判斷。相比于切線剛度法，割線剛度迭代法簡(jiǎn)化了計(jì)算過(guò)程。3. 切線剛度法的過(guò)程，需要通過(guò)外荷載節(jié)點(diǎn)力和結(jié)構(gòu)恢復(fù)力求解出一個(gè)殘余力，再以此計(jì)算出增量位移來(lái)判斷是否滿足收斂條件，若不收斂，便更新結(jié)構(gòu)位移并重新計(jì)算結(jié)構(gòu)的切線剛度矩陣。對(duì)于割線剛度迭代法，并不需要求解殘余力，

28、而是直接利用外荷載節(jié)點(diǎn)力和結(jié)構(gòu)割線剛度來(lái)求解出一個(gè)全量位移，之后通過(guò)計(jì)算出的應(yīng)變來(lái)重新計(jì)算結(jié)構(gòu)的割線剛度矩陣。二者相比起來(lái)，割線剛度迭代法不需要啟用一個(gè)新的數(shù)組來(lái)儲(chǔ)存殘余力值，降低了對(duì)程序的存儲(chǔ)空間要求并提高了程序的運(yùn)算速度。從上面三點(diǎn)可以看出，切線剛度法和割線剛度迭代法各有特點(diǎn)和不足，在某些方面，割線剛度迭代法還體現(xiàn)出了一定的優(yōu)勢(shì)：割線剛度迭代法相比于切線剛度法的收斂準(zhǔn)則更為簡(jiǎn)單。割線剛度迭代法相比于切線剛度法的有限元模型更為簡(jiǎn)單。割線剛度迭代法的穩(wěn)定性較好。所以對(duì)割線剛度迭代法進(jìn)行研究佐證是十分有必要的。2.5.1 算法流程以圖 所示的單元為例，Vecchio在其論文中所述的割線剛度迭代法

29、的具體計(jì)算過(guò)程如下：圖 3.5 單元尺寸Fig3.5 Element dimension第1步：輸入結(jié)構(gòu)和材料的屬性數(shù)據(jù)第2步：計(jì)算外荷載節(jié)點(diǎn)力第3步：計(jì)算材料的剛度矩陣：1) 確定初始應(yīng)變：，； 2) 計(jì)算主方向應(yīng)變，以及主拉應(yīng)變角；3) 根據(jù)混凝土的本構(gòu)關(guān)系，確定混凝土主拉和主壓應(yīng)力，4) 根據(jù)鋼筋的本構(gòu)關(guān)系，確定鋼筋的主拉和主壓應(yīng)力，；5) 計(jì)算混凝土和鋼筋的割線剛度元素值： 第4步：計(jì)算混凝土和鋼筋的材料割線剛度矩陣1) 計(jì)算混凝土主方向上的材料割線剛度矩陣： (3-2)2) 計(jì)算轉(zhuǎn)換矩陣3) 計(jì)算整體坐標(biāo)系（X-Y坐標(biāo)）下的混凝土割線剛度矩陣 (3-3)4) 計(jì)算整體坐標(biāo)系(X-Y坐

30、標(biāo)系)下的鋼筋割線剛度矩陣： (3-4)第5步：計(jì)算組合材料剛度： (3-5)第6步：計(jì)算單元?jiǎng)偠染仃嚕?(3-6)第7步：計(jì)算單元預(yù)應(yīng)變產(chǎn)生的位移第8步：計(jì)算單元預(yù)應(yīng)變產(chǎn)生的節(jié)點(diǎn)力第9步：計(jì)算單元節(jié)點(diǎn)力：第10步：計(jì)算結(jié)構(gòu)總剛矩陣： (3-7) 第11步：計(jì)算結(jié)構(gòu)總剛矩陣的逆矩陣：第12步：計(jì)算節(jié)點(diǎn)位移 (3-8) 第13步：計(jì)算單元應(yīng)變（取單元形心處的應(yīng)變?yōu)槠骄鶓?yīng)變）： (3-9) (3-10) (3-11)注：，為相應(yīng)節(jié)點(diǎn)的水平位移和豎向位移。第14步：檢查單元應(yīng)力：第15步：重新計(jì)算單元割線剛度值：1) 根據(jù)求解出的位移，重新確定主方向上的應(yīng)變值，；3) 根據(jù)混凝土本構(gòu)模型，計(jì)算新的混

31、凝土應(yīng)力，；4) 根據(jù)鋼筋本構(gòu)模型，計(jì)算新的鋼筋應(yīng)力，；5) 重新計(jì)算鋼筋和混凝土的割線剛度元素：，,。第16步：收斂判斷詳細(xì)的算法流程圖如圖 所示。割線剛度迭代法的迭代過(guò)程與基于切線剛度的Newton-Raphson迭代法有一定的可比性，其具體的迭代過(guò)程圖線如圖 所示。每次迭代通過(guò)外荷載與即時(shí)的割線剛度求解出全量位移，直到與真實(shí)位移的誤差小于允許值，迭代計(jì)算結(jié)束。2.5.2 收斂準(zhǔn)則不同于基于切線剛度的Newton-Raphson迭代法，割線剛度迭代法并沒(méi)有采用前面所述及的三種收斂準(zhǔn)則。由于這三種收斂準(zhǔn)則都是貫徹了增量策略的思想，對(duì)于應(yīng)用全量理論的割線剛度迭代法而言，顯然不太合適。因此，割線

32、剛度迭代法采用了一種全新的收斂準(zhǔn)則：以單元的割線剛度作為判斷收斂的條件。如果某前后兩個(gè)迭代步得到的割線剛度誤差小于預(yù)定的限制，那么迭代終止。對(duì)于多個(gè)單元的情況而言，割線剛度收斂準(zhǔn)則可以表現(xiàn)為不同的形式。其中一種是以各個(gè)單元的割線剛度誤差的最大值作為判斷標(biāo)準(zhǔn)，即其中，是向量的最大范數(shù)，分別是第個(gè)單元的混凝土兩個(gè)主方向上的迭代步前后兩次割線剛度誤差，是預(yù)設(shè)的誤差限。另外一種形式是以各個(gè)單元的割線剛度誤差之和作為判斷準(zhǔn)則，數(shù)學(xué)表達(dá)式為：其中，是向量的1-范數(shù)，是預(yù)設(shè)的誤差限。 本文程序選用的是第二種形式的收斂準(zhǔn)則。3.1 FEAPpv簡(jiǎn)介在過(guò)去的幾十年里，有限元法從一種線性結(jié)構(gòu)分析方法演變成了一種求

33、解非線性微分方程的通用方法，F(xiàn)EAPpv便是應(yīng)運(yùn)而生的一個(gè)有限元分析程序（Finite Element Analysis Program）。在研究和應(yīng)用領(lǐng)域中，需要不斷對(duì)其進(jìn)行改進(jìn)來(lái)解決新的問(wèn)題和滿足新的分析要求。FEAPpv程序主要開發(fā)于UNIX平臺(tái)和Window平臺(tái)，其包括一系列模塊用于實(shí)現(xiàn)：1.有限單元模型數(shù)據(jù)的輸入；2. 用于處理大范圍應(yīng)用實(shí)例的求解方法的建立；3. 結(jié)果的圖形和數(shù)值輸出。一個(gè)問(wèn)題的求解，是通過(guò)使用命令語(yǔ)句的形式來(lái)達(dá)到的，其中的求解算法則由用戶定義并編寫。因此，F(xiàn)EAPpv的這一特點(diǎn)使得每一個(gè)用戶可以根據(jù)自己的具體需要來(lái)定義一種特定的求解方式，以此達(dá)到分析研究目的。FE

34、APpv程序自身包含了足夠多的命令來(lái)應(yīng)用于結(jié)構(gòu)力學(xué)、流體力學(xué)、熱傳導(dǎo)以及許多其他需要用微分方程來(lái)模擬問(wèn)題的領(lǐng)域，靜態(tài)和動(dòng)態(tài)問(wèn)題都可以用該程序求解。根據(jù)不同的應(yīng)用需求，用戶也可以自行添加新的模型和命令語(yǔ)句，這些功能可以幫助用戶對(duì)一些特殊結(jié)構(gòu)進(jìn)行單元的劃分，更或者使用戶能從其他程序中導(dǎo)入單元?jiǎng)澐中畔?。FEAPpv程序包含一個(gè)有限的單元庫(kù)。這些單元可以用于一維、二維和三維線性或非線性結(jié)構(gòu)問(wèn)題、固體力學(xué)問(wèn)題以及線性熱傳導(dǎo)問(wèn)題的建模，并調(diào)用一個(gè)材料模型庫(kù)。所提供的材料模型適用于彈性、粘滯彈性、塑性以及熱傳導(dǎo)本構(gòu)方程，用戶也可以根據(jù)自己的需要添加另外的模型。在結(jié)構(gòu)問(wèn)題中，程序能夠通過(guò)單元形成質(zhì)量和幾何剛度

35、矩陣，并計(jì)算單元上的相關(guān)值，如應(yīng)力和應(yīng)變，并能夠把相關(guān)值投影到結(jié)點(diǎn)上，輸出圖形結(jié)果。FEAPpv程序可分為三個(gè)基本模塊：數(shù)據(jù)輸入模塊和先處理器求解模塊輸出模塊以下是簡(jiǎn)化的FEAPpv程序流程圖：上圖中的程序輸入模塊將讀入輸入文件中的材料（MATE）、幾何（COOR， ELEM）、邊界條件（BOUN）、荷載（FORC）以及位移（DISP）信息，用以建立所有的單元數(shù)組。之后，材料參數(shù)、結(jié)點(diǎn)荷載和結(jié)點(diǎn)位移、方程編號(hào)、溫度、單元連接以及坐標(biāo)等分別存儲(chǔ)在程序中的不同數(shù)組里。表 列出了一些主要數(shù)組的名稱、指針編號(hào)和它們各個(gè)維度的大小。表3.1 網(wǎng)格劃分?jǐn)?shù)組名稱、指針編號(hào)和數(shù)組的大小32Table 3.1

36、Mesh Array Names, Numbers and Sizes32NAMENum.dim 1dim 2dim 3DescriptionD25nddnummat-Material parametersF27ndfnumnp2Force and DisplacementID31ndfnumnp2Equation nos.IE32nienummat-Element control, dofs, etc.IX33nen1numel-Element connectionsT38numnp-TemperatureU40ndfnumnp3Solution arrayVEL42ndfnumnpntSo

37、lution rate arrayX43ndmnumnp-Coordinates本文主要涉及和研究FEAPpv的控制程序，至于詳細(xì)的前處理模塊介紹，可以參見文獻(xiàn)。下圖為控制程序的簡(jiǎn)化流程圖?？梢钥闯?，整個(gè)控制程序由網(wǎng)格劃分（PMESH）和宏命令（PMACR）求解兩大部分組成。FEAPpv的宏（MACR）命令流包含與6個(gè)子程序，分別是PMACR、PMACIO、PMACR1、PMACR2、PMACR3和PMACR4，其中PMACR4是一個(gè)用戶自定義宏命令子程序，不使用時(shí)為空，且不會(huì)對(duì)整個(gè)程序運(yùn)行有影響。這幾個(gè)宏命令子程序的具體任務(wù)和作用是：PMACR位于流程上端，負(fù)責(zé)調(diào)用其他幾個(gè)子程序；PMACI

38、O負(fù)責(zé)讀入命令流，并檢查宏命令的循環(huán)和嵌套；PMACR1主要負(fù)責(zé)程序的計(jì)算部分，如計(jì)算剛度和殘余力；PMACR2主要負(fù)責(zé)控制計(jì)算過(guò)程的循環(huán)、荷載步長(zhǎng)以及收斂誤差的大小等；PMACR3則負(fù)責(zé)完成求解和繪圖等其他的功能。圖 所示為整個(gè)宏命令的結(jié)構(gòu)圖：3.2 鋼筋混凝土基于切線剛度法的非線性分析流程圖FEAPpv的源程序采用的是典型的基于切線剛度的有限元計(jì)算方法。切線剛度法的特點(diǎn)是：需要通過(guò)輸入信息來(lái)求得單元的結(jié)點(diǎn)恢復(fù)力，并集成結(jié)構(gòu)的結(jié)點(diǎn)恢復(fù)力，再與外荷載產(chǎn)生的結(jié)點(diǎn)力求殘差，確定出殘余力的大小；根據(jù)殘余力和形成的總剛計(jì)算結(jié)構(gòu)的增量位移，并更新總位移；以力、位移或者能量作為準(zhǔn)則判斷是否收斂。鋼筋混凝土

39、基于切線剛度法的二維非線性分析具體流程圖如下所示：3.3 鋼筋混凝土基于割線剛度法的非線性分析流程圖本文所研究的鋼筋混凝土基于割線剛度的二維非線性有限元分析流程圖如下所示可以看出，割線剛度法和切線剛度法算法流程的差別在于：1. 割線剛度法不需要計(jì)算單元恢復(fù)力以及結(jié)構(gòu)的殘余力，從而不必分配大型的數(shù)組用以存取恢復(fù)力和殘余力，節(jié)省了計(jì)算時(shí)間和儲(chǔ)存空間，提高了效率；2. 割線剛度法直接利用外荷載在結(jié)點(diǎn)上產(chǎn)生的力和集成的結(jié)構(gòu)總體剛度矩陣來(lái)求解結(jié)構(gòu)的全量位移，而非切線剛度法中通過(guò)殘余力求得的增量位移；3. 割線剛度法以單元割線剛度元素作為收斂條件，相比于切線剛度法的收斂判斷準(zhǔn)則更為簡(jiǎn)潔。3.4 單元子程序

40、的數(shù)據(jù)傳遞和結(jié)構(gòu)用戶可以通過(guò)ELMTnn這個(gè)單元子程序來(lái)向FEAPpv的單元庫(kù)添加自己所需要的單元模塊：subroutine ELMTnn (d, ul, xl, ix, tl, s, r, ndf, ndm, nst, isw)其中nn的值從01到05。單元子程序中的各種信息，是通過(guò)啞元和公共變量進(jìn)行傳遞的。表 列出了程序中的一些啞元以及它們所傳遞的數(shù)據(jù)信息，圖 列出了一部分公共塊的名稱及其所含的參數(shù)變量。表3.2 FEAPpv單元子程序的啞元32Table 3.2 Arguments of FEAPpv Element Subprogram32ParameterDescriptiond(*

41、)Element data parametersul(ndf,nen,j)Element nodal solution parameters nen is number of nodes on an element (max)xl(ndm,nen)Element nodal reference coordinatesix(nen)Element global node numberstl(nen)Element nodal temperature valuess(nst,nst)Element matrix (e.g., stiffness, mass)r(ndf,nen)Element ve

42、ctor (e.g., residual, mass)may also be used as r(nst)ndfNumber unknowns (max) per nodendmSpace dimension of meshnstSize of element arrays S and RN.B. Normally nst = ndf*neniswTask parameter to control computationinteger numnp, numel, nummat, nen, neq,iprcommon /cdata/ numnp, numel, nummat, nen, neq,

43、iprreal*8 dminteger n, ma, mct, iel, nelcommon /eldata/ dm, n, ma, mct, iel, nelinteger nh1,nh2,nh3common /hdata/ nh1,nh2,nh3integer ior,iowcommon /iofile/ ior,iowinteger ndf,ndm,nen1,nst,nneqcommon /sdata/ ndf,ndm,nen1,nst,nneqreal*8 hrinteger mrcommon / / hr(1),mr(1000)圖3.5 FEAPpv 單元公共塊32Figure 3.

44、5 FEAPpv Element Common Blocks32FEAPpv中單元子程序的簡(jiǎn)單結(jié)構(gòu)圖如圖 所示。圖 與切線剛度法對(duì)應(yīng)，圖 與割線剛度法對(duì)應(yīng)。其中ISW是FEAPpv程序中起控制作用的任務(wù)參數(shù)，相應(yīng)于不同取值的ISW的具體任務(wù)見表表3.3 FEAPpv單元子程序任務(wù)選項(xiàng)32Table 3.3 Task Options for FEAPpv Element Subprogram32 isw taskDescriptionAccess Command1Input d(*) parametersMesh:MATE,n2Check elementsSoln:CHECk3Compute t

45、angent/residual Soln:TANGStore in S/rSoln:TANG4Output element variablesSoln:FORM,REACPlot:REAC6Compute residualSoln:FORM,REACPlot:REAC8Nodal projectionsSoln:STRE,NODEPlot:STRE,PSTR12History updateSoln:TIME14Initialize historyBATCh,INTEr通過(guò)對(duì)比兩個(gè)結(jié)構(gòu)圖，也可以看出，割線剛度法省去了計(jì)算恢復(fù)力的步驟，更為簡(jiǎn)便。3.4.1 單元子程序中材料參數(shù)的存儲(chǔ)在程序中，單元

46、的各種材料參數(shù)存儲(chǔ)的一維數(shù)組D(*)里。D(*)數(shù)組中各個(gè)元素的具體代表值見表表3.1 鋼筋混凝土單元材料參數(shù)Table 3.1 Material Parameters of the RC ElementParameterDescriptionD(1)混凝土的初始彈性模量 D(2）混凝土的初始泊松比D(3)混凝土抗壓峰值強(qiáng)度D(4)混凝土的受壓峰值應(yīng)變D(5)混凝土受拉開裂強(qiáng)度D(7)鋼筋彈性模量D(8)鋼筋的應(yīng)變強(qiáng)化模量D(10)x向鋼筋的屈服強(qiáng)度D(11)y向鋼筋的屈服強(qiáng)度D(12)x向鋼筋的配筋率D(13)y向鋼筋的配筋率D(14）單元的厚度D(181)X向鋼筋直徑D(182)Y向鋼筋直

47、徑D(183)混凝土骨料尺寸3.4.2 單元應(yīng)力和材料矩陣的形成過(guò)程鋼筋混凝土單元的應(yīng)力包含兩個(gè)部分：鋼筋應(yīng)力和混凝土應(yīng)力。同樣的，鋼筋混凝土單元的材料矩陣也由鋼筋的材料矩陣和混凝土的材料矩陣組合而成。圖 所示的是單元應(yīng)力和材料矩陣形成的過(guò)程圖。對(duì)于切線剛度法而言，圖示的流程是針對(duì)每個(gè)高斯積分點(diǎn)進(jìn)行的，而在割線剛度法中，整個(gè)過(guò)程在單元的中心點(diǎn)得以實(shí)現(xiàn)。3.4.3 單元?jiǎng)偠染仃嚨男纬汕芯€剛度法采用的是FEAPpv程序默認(rèn)的高斯點(diǎn)數(shù)，即用22的高斯點(diǎn)進(jìn)行積分，其具體的單元?jiǎng)偠染仃囆纬蛇^(guò)程見圖對(duì)于割線剛度法，在每個(gè)單元上僅取其形心點(diǎn)處的應(yīng)變作為整個(gè)單元的平均應(yīng)變，一定程度上簡(jiǎn)化了單元?jiǎng)偠刃纬傻挠?jì)算過(guò)

48、程。單元?dú)v史狀態(tài)的存儲(chǔ) 在進(jìn)行低周反復(fù)加載的非線性分析時(shí)，對(duì)單元的歷史狀態(tài)進(jìn)行存儲(chǔ)甚為重要。比如結(jié)構(gòu)在經(jīng)過(guò)了加卸載之后，會(huì)產(chǎn)生不可恢復(fù)的塑性應(yīng)變，而下一荷載步通過(guò)公式計(jì)算得到的塑性應(yīng)變值很有可能小于上一荷載步計(jì)算所得值，若程序沒(méi)有記憶功能，那么小值將會(huì)把大值覆蓋，得出不合理的結(jié)果。所以，在程序中需要設(shè)置一些這種參數(shù)：它們負(fù)責(zé)把單元有用的歷史變量存儲(chǔ)起來(lái)，如上一荷載步的應(yīng)變以及加載過(guò)程中最大的塑性應(yīng)變等?；炷羻卧谀硞€(gè)迭代收斂步處所設(shè)置的歷史變量如下表所示：表3.2 混凝土單元材料參數(shù)h(*)Table 3.2 Material Parameters of the Concret Elemen

49、tParameterDescriptionh(1)混凝土主應(yīng)力角h(2)混凝土X方向塑性應(yīng)變h(3)混凝土Y方向塑性應(yīng)變h(4)混凝土塑性應(yīng)變h(5)混凝土X方向最大壓應(yīng)變h(6)混凝土Y方向最大壓應(yīng)變h(7)混凝土最大壓應(yīng)變h(8)混凝土X方向最大拉應(yīng)變h(9)混凝土Y方向最大拉應(yīng)變h(10)混凝土最大拉應(yīng)變h(11)混凝土X方向上一步應(yīng)變h(12)混凝土Y方向上一步應(yīng)變h(13)混凝土上一步應(yīng)變鋼筋單元在某個(gè)迭代收斂步處所設(shè)置的歷史變量如下表所示：表3.3 鋼筋單元材料參數(shù)h(*)Table 3.3 Material Parameters of the Steel ElementParam

50、eterDescriptionh(18)，h（19）X，Y向鋼筋上一步路徑h(20)，h(21)X,Y向鋼筋上一步應(yīng)變h(22)，h（23）X,Y向鋼筋上一步應(yīng)力h(24)，h（25）X,Y向鋼筋上一步塑性應(yīng)變h(26)，h（27）X,Y向鋼筋漸近線交點(diǎn)應(yīng)變h(28)，h（29）X,Y向鋼筋漸近線交點(diǎn)應(yīng)力h(30)，h（31）X,Y向鋼筋應(yīng)變反向點(diǎn)應(yīng)變h(32)，h（33）X,Y向鋼筋應(yīng)變反向點(diǎn)應(yīng)力h(34)，h（35）X,Y向鋼筋彈性模量h(36)，h（37）X,Y向鋼筋歷史最大應(yīng)變h(38)，h（39）X,Y向鋼筋歷史最小應(yīng)變3.5 位移加載模式30在鋼筋混凝土的非線性有限元分析中，加載控

51、制方法一般可以分為力加載和位移加載兩種。對(duì)于剪切板的模擬，過(guò)去通常采用的是力加載模式。當(dāng)構(gòu)件出現(xiàn)局部軟化、平衡路徑分歧和到達(dá)極值點(diǎn)的情況時(shí)，力加載模式就暴露出了一定的缺陷：不管是采用切線剛度法還是割線剛度法，都有可能得到不收斂的結(jié)果，更不能求得應(yīng)力-應(yīng)變關(guān)系曲線的下降段。通過(guò)以前的一些分析計(jì)算結(jié)果可以看出，位移加載在一定程度上可以解決上述問(wèn)題，所以本文在割線剛度法的基礎(chǔ)上，對(duì)剪切板進(jìn)行了位移加載方式的研究。FEAPpv提供了一種支座位移加載的模式，其具體方法是把結(jié)點(diǎn)需要加位移的自由度先約束住，視為一個(gè)支座，然后在此自由度方向施加支座位移，最后根據(jù)結(jié)構(gòu)剛度矩陣換算成相應(yīng)的結(jié)點(diǎn)力進(jìn)行求解。現(xiàn)以一個(gè)

52、單元為例來(lái)說(shuō)明整個(gè)支座位移加載的實(shí)現(xiàn)過(guò)程。 (a) (b)圖3.6 加支座位移示意圖Fig 3.6 Flow chart of Support Displacement Method如圖3.6(a)所示一個(gè)四邊形單元，結(jié)點(diǎn)1和結(jié)點(diǎn)2為支座。若要在結(jié)點(diǎn)4上施加一個(gè)水平位移，就需要將結(jié)點(diǎn)4的方向自由度約束住，如圖3.6(b)所示，然后對(duì)結(jié)點(diǎn)4的施加一個(gè)水平向的支座位移。設(shè)支座位移向量為，考慮到約束條件，求解方程中的剛度矩陣就會(huì)只剩下第3、5、6、8自由度上的分量， (3.17) (3.18) (3.19)式(3.18)和式(3.19)分別為考慮了約束條件之后得到的剛度矩陣和位移向量的表達(dá)式?？梢钥?/p>

53、出，求解的關(guān)鍵在于如何確定方程中結(jié)點(diǎn)荷載向量的表達(dá)式。在FEAPpv中，結(jié)點(diǎn)荷載是通過(guò)剛度矩陣與支座位移計(jì)算得到的。這個(gè)剛度矩陣的表達(dá)式為： (3.20)上述剛度矩陣表達(dá)的意思是，支座位移在非支座結(jié)點(diǎn)上產(chǎn)生的力。相應(yīng)的結(jié)點(diǎn)荷載向量就可以由下式得到： (3.21)需要注意的是，在采用基于切線剛度的Newton-Raphson迭代法時(shí)，以上求結(jié)點(diǎn)力向量的過(guò)程只需要在每一荷載步的第一次迭代時(shí)進(jìn)行，在此之后的迭代步中，由于求解方程的等號(hào)右端項(xiàng)為不平衡力，因此不需要再重復(fù)計(jì)算支座位移引起的結(jié)點(diǎn)荷載。3.6 二維非線性有限元分析程序VecTor2的介紹 對(duì)于剪切板的位移加載模式，并沒(méi)有現(xiàn)成的實(shí)驗(yàn)數(shù)據(jù)供參考

54、比較，所以需要借助一些其他方法來(lái)驗(yàn)證所編制的采用割線剛度法的FEAPpv程序計(jì)算分析結(jié)果的正確性。VecTor2便是一個(gè)非常好的參照程序。VecTor2是一個(gè)用于鋼筋混凝土結(jié)構(gòu)二維非線性有限元分析的程序，由多倫多大學(xué)在過(guò)去幾十年里致力于研究鋼筋混凝土性能和有限元法應(yīng)用的專家們研究開發(fā)。VecTor2的理論基礎(chǔ)是修正斜壓場(chǎng)理論（MCFT）以及擾動(dòng)應(yīng)力場(chǎng)模型（DSFM），其用于預(yù)測(cè)鋼筋混凝土構(gòu)件反應(yīng)的分析模型承受的是平面內(nèi)的正應(yīng)力和剪應(yīng)力。VecTor2將開裂混凝土模擬為具有彌散和旋轉(zhuǎn)裂縫的正交各向異性材料，并采用割線剛度法來(lái)得到一個(gè)有效可靠的非線性解。VecTor2原本采用的是MCFT通過(guò)鋼筋混

55、凝土平板純剪試驗(yàn)的軟化曲線得到的本構(gòu)關(guān)系，之后針對(duì)各種二階效應(yīng)，如受壓軟化、受拉硬化、受拉軟化以及受拉滑移，加入了一些其他的本構(gòu)模型。此外，VecTor2的分析能力還擴(kuò)展到了對(duì)混凝土的膨脹和約束、粘結(jié)滑移、裂縫剪切滑移變形、鋼筋銷栓作用、鋼筋屈曲等方面的模擬。VecTor2所使用的有限元模型被良好地劃分為一些簡(jiǎn)單單元，這種方法有利于計(jì)算的有效性和數(shù)值的穩(wěn)定性。對(duì)于那些需要通過(guò)相對(duì)仔細(xì)的單元?jiǎng)澐謥?lái)模擬配筋和裂縫形式的鋼筋混凝土結(jié)構(gòu)來(lái)說(shuō)，VecTor2特別適用。VecTor2的單元庫(kù)里包括了連續(xù)應(yīng)變?nèi)切螁卧?、平面?yīng)力矩形單元、用于模擬將鋼筋彌散化的混凝土的四邊形單元、用于模擬分離鋼筋的兩結(jié)點(diǎn)桿單

56、元以及用于模擬粘結(jié)滑移機(jī)理的兩結(jié)點(diǎn)連接單元和四結(jié)點(diǎn)接觸單元。 VecTor2的程序流程圖如下圖所示：由上圖可以看出，VecTor2的程序算法和F. J. Vecchio論文中所述及的割線剛度法過(guò)程是契合的，具有很高的對(duì)比和參考價(jià)值。其中需要注意的一點(diǎn)是，在VecTor2程序里，參數(shù)是由鋼筋和混凝土的塑性應(yīng)變所確定的位移轉(zhuǎn)化而成的結(jié)點(diǎn)力，而F. J. Vecchio采用的割線剛度法則是將單元的預(yù)應(yīng)變轉(zhuǎn)換為作用于結(jié)點(diǎn)的力，然后存儲(chǔ)在這個(gè)向量當(dāng)中。這是二者算法中的不同點(diǎn)。 VecTor2的算法流程為本論文在FEAPpv平臺(tái)上進(jìn)行的割線剛度法的研究和編制指明了方向，也為之后的程序驗(yàn)證工作提供了一個(gè)強(qiáng)有

57、力的輔助工具。3.7 本文所做工作簡(jiǎn)述鄧興龍和蔣曉華等所編制的FEAPpv程序是按照的基于切線剛度的Newton-Raphson迭代法的計(jì)算思路。為了實(shí)現(xiàn)整個(gè)算法的改變，將F. J. Vecchio采用的割線剛度法應(yīng)用于FEAPpv平臺(tái)，本文做了以下的一些主要工作：1. 將FEAPpv程序中的單元切線剛度全部修改為割線剛度。2. 外荷載產(chǎn)生的結(jié)點(diǎn)力的提取和處理。切線剛度迭代法采用殘余力求得的是一個(gè)增量位移，而為了得到結(jié)構(gòu)的全量位移，就需要把外荷載所產(chǎn)生的結(jié)構(gòu)結(jié)點(diǎn)力單獨(dú)提取出來(lái)進(jìn)行位移計(jì)算。由于在過(guò)去的程序中，外荷載結(jié)點(diǎn)力是以一個(gè)完整的形式存儲(chǔ)于數(shù)組里，即不考慮約束條件，若直接代入求解子程序，會(huì)

58、產(chǎn)生錯(cuò)誤。所以此次的工作除了提取外荷載結(jié)點(diǎn)力之外，還根據(jù)結(jié)構(gòu)已知的邊界約束條件，對(duì)外荷載結(jié)點(diǎn)力進(jìn)行處理，最后導(dǎo)入到程序中的求解子程序，計(jì)算出所需要得結(jié)構(gòu)全量位移。3. 位移的更新。切線剛度法用已知位移加上增量位移來(lái)更新結(jié)構(gòu)的位移，對(duì)于割線剛度法，則用求解出來(lái)的全量位移來(lái)直接替代已知位移，作為結(jié)構(gòu)的總位移。4. 收斂準(zhǔn)則的修改。切線剛度迭代法一般選用力、位移或者能量作為收斂準(zhǔn)則。若要應(yīng)用割線剛度法，則要相應(yīng)地對(duì)收斂判斷準(zhǔn)則進(jìn)行修改。本文以單元割線剛度誤差之和作為判斷收斂的條件，并在宏命令求解子程序中加入了與之呼應(yīng)的相關(guān)FORTRAN語(yǔ)句。5. 位移加載方式的程序完善。采用切線剛度迭代法時(shí)，一些參

59、數(shù)是以增量的形式存在于程序中的，如支座位移。在進(jìn)行支座位移加載計(jì)算時(shí)，切線剛度迭代法的做法是：通過(guò)比較某個(gè)荷載步內(nèi)前后兩次迭代步的位移增量的大小，來(lái)判斷是否需要修正不平衡力。若位移增量值大于設(shè)定值，則計(jì)算出一個(gè)增量力來(lái)修正結(jié)構(gòu)殘余力。但對(duì)于割線剛度法來(lái)說(shuō)，增量參數(shù)對(duì)于整個(gè)算法是無(wú)能為力的，需要替換成全量值。本文的工作是把支座位移單獨(dú)提取出來(lái)，存放在一個(gè)新數(shù)組當(dāng)中，然后將其帶入結(jié)點(diǎn)力修正子程序，換算得到一組支座位移產(chǎn)生的全量結(jié)點(diǎn)力，最后把該結(jié)點(diǎn)力向量傳入程序存儲(chǔ)外荷載結(jié)點(diǎn)力的指針數(shù)組，就實(shí)現(xiàn)了對(duì)支座位移加載的處理。如在本文3.7節(jié)位移加載模式中述及的，對(duì)于基于切線剛度的Newton-Raphso

60、n迭代法，僅需要在各個(gè)荷載步的第一個(gè)迭代步計(jì)算支座位移產(chǎn)生的結(jié)點(diǎn)力向量，之后的迭代步由于都是采用增量形式，以不平衡力來(lái)求解方程，所以不需要再重復(fù)計(jì)算支座位移引起的結(jié)點(diǎn)力。而割線剛度法是一種全量方法，求解方程的右端始終都是施加的等效結(jié)點(diǎn)力，所以對(duì)于某個(gè)荷載步而言，其中的每一個(gè)迭代步都應(yīng)該計(jì)算出相應(yīng)的支座位移產(chǎn)生的結(jié)點(diǎn)力向量。此次工作在FEAPpv的位移加載模塊中做了一定的修改，使每一次迭代步都求出即時(shí)的割線剛度矩陣，用以計(jì)算和更新支座位移引起的結(jié)點(diǎn)力，實(shí)現(xiàn)方程的正確求解。4.1 簡(jiǎn)述 利用編制的割線剛度迭代法FEAPpv程序，本文對(duì)幾個(gè)系列的剪切板進(jìn)行了單調(diào)加載和反復(fù)加載下的分析計(jì)算。對(duì)于單調(diào)加