【正文】 這 個 剛度矩陣的表達式為： ? ?31 32 34 3752 54 5764 6787k k k kk k kKkkk?????????對 稱 () 上述剛度矩陣表達的意思是，支座位移在非支座結(jié)點上產(chǎn)生的 力?？梢钥闯?，求解的 關(guān)鍵在于如何確定方程中 結(jié)點荷載向量 ??P 的表達式。若要在 結(jié)點 4上 施加一個 水平位移 ，就需要 將結(jié)點 4的 x 方向 自由度約束住 ，如圖 (b)所 示 ， 然后 對結(jié)點 4的施加 一個水平向的 支座位移 ?，F(xiàn)以一個單元為例來說明整個支座位移加載的實現(xiàn)過程。通過以前的一些分析計算結(jié)果 可以看出，位移加載在一定程度上可以解決上述問題，所以本文在割線剛度法的基礎(chǔ)上，對剪切板進行了位移加載方式的研究。對于剪切板的模擬，過去通常采用的是力加載模式。tmy? h(10) 混凝土最大拉應(yīng)變 39。cmxy? h(8) 混凝土 X 方向最大拉應(yīng)變 39。cmx? h(6) 混凝土 Y 方向最大壓應(yīng)變 39。pcy? h(4) 混凝土塑性應(yīng)變 39。 混凝土單元在某個迭代收斂 步處所設(shè)置的歷史變量如下表所示： 表 混凝土單元材料參數(shù) h(*) Table Material Parameters of the Concret Element Parameter Description h(1) 混凝土主應(yīng)力角 ? h(2) 混凝土 X 方向塑性應(yīng)變 39。比如結(jié)構(gòu)在經(jīng)過了加卸載之后，會產(chǎn)生不可恢復(fù)的塑性應(yīng)變，而下一荷載步通過公式計算得到的塑性應(yīng)變值很有可能小于上一荷載步計算所得值，若程序沒有記憶功能，那么小值將會把大值覆蓋，得出不合理的結(jié)果。 單元剛度矩陣的形成 切線剛度法采用的是 FEAPpv 程序默認(rèn)的高斯點數(shù)，即用 2 2 的高斯點進行積分，其具體的單元剛度矩陣形成過程見圖 混凝土應(yīng)力： ζ cx , ζ cy, η cxy 混凝土的材料矩陣： [D c ] 鋼筋混凝土的總應(yīng)力 ζ x= ζ cx+ ρ x ζ sx ζ y= ζ cy+ ρ y ζ sy η xy= η cxy鋼筋混凝土的總材料剛度矩陣 [D] = [D c ] + ρ x [D s ] x + ρ y [D s ] y 鋼筋應(yīng)力： ζ sx , ζ sy 鋼筋的材料矩陣： [D s ] x , [D s ] y 混凝土應(yīng)變： ε cx = ε x ε cy = ε y γ cxy = γ xy 調(diào)用混凝土模型 鋼筋應(yīng)變： ε sx = ε x ε sy = ε y 調(diào)用鋼筋模型 ε x, ε y, γ xy 對于割線剛度法，在每個單元上僅取其形心點 處的應(yīng)變作為整個單元的平均應(yīng)變，一定程度上簡化了單元剛度形成的計算過程。圖 所示的是單元應(yīng)力和材料矩陣形成的過程圖。c? D(5) 混凝土受拉開裂強度 crf D(7) 鋼筋彈性模量 sE D(8) 鋼筋的應(yīng)變強化模量 b D(10) x向鋼筋的屈 服強度 yxf D(11) y 向鋼筋的屈服強度 yyf D(12) x向鋼筋的配筋率 x? D(13) y 向鋼筋的配筋率 y? D(14） 單元的厚度 t D(181) X 向鋼筋直徑 D(182) Y 向鋼筋直徑 D(183) 混凝土 骨料尺寸 單 元應(yīng)力和材料矩陣的形成過程 鋼筋混凝土單元的應(yīng)力包含兩個部分：鋼筋應(yīng)力和混凝土應(yīng)力。 D(*)數(shù)組中各個元素的具體代表值見表 表 鋼筋混凝土單元材料參數(shù) Table Material Parameters of the RC Element Parameter Description D(1) 混凝土的初始彈性模量 cE D(2） 混凝土的初始泊 松比 ? D(3) 混凝土抗壓峰值強度 39。其中 ISW是 FEAPpv程序中起控制作用的任務(wù)參數(shù)，相應(yīng)于不同取值的 ISW的具體任務(wù)見表 表 FEAPpv單元子程序任務(wù)選項 [32] Table Task Options for FEAPpv Element Subprogram[32] isw task Description Access Command 1 Input d(*) parameters Mesh:MATE,n 2 Check elements Soln:CHECk 3 Compute tangent/residual Soln:TANG Store in S/r Soln:TANG 4 Output element variables Soln:FORM,REAC Plot:REAC 6 Compute residual Soln:FORM,REAC Plot:REAC 8 Nodal projections Soln:STRE,NODE Plot:STRE,PSTR 12 History update Soln:TIME 14 Initialize history BATCh,INTEr f a l s et r u ef a l s et r u ef a l s et r u et r u ef a l s et r u ef a l s e返回輸出單元結(jié)點的應(yīng)力i s w = 8形成單元結(jié)點的恢復(fù)力列陣i s w = 6輸出單元的應(yīng)力、應(yīng)變i s w = 4i s w = 3i s w = 1形成單元剛度矩陣輸入材料信息d ( * ) f al set r uef al set r uet r uef al set r uef al se返回輸出單元結(jié)點的應(yīng)力i sw =8輸出單元的應(yīng)力、應(yīng)變i sw =4i sw =3i sw =1形成單元剛度矩陣輸入材料信息d(*) 通過對比兩個結(jié)構(gòu)圖，也可以看出，割線剛度法省去了計算恢復(fù)力的步驟，更為簡便。 表 FEAPpv單元子程序的啞元 [32] Table Arguments of FEAPpv Element Subprogram[32] Parameter Description d(*) Element data parameters ul(ndf,nen,j) Element nodal solution parameters nen is number of nodes on an element (max) xl(ndm,nen) Element nodal reference coordinates ix(nen) Element global node numbers tl(nen) Element nodal temperature values s(nst,nst) Element matrix (., stiffness, mass) r(ndf,nen) Element vector (., residual, mass) may also be used as r(nst) ndf Number unknowns (max) per node ndm Space dimension of mesh nst Size of element arrays S and R . Normally nst = ndf*nen isw Task parameter to control putation integer numnp, numel, nummat, nen, neq,ipr mon /cdata/ numnp, numel, nummat, nen, neq,ipr real*8 dm integer n, ma, mct, iel, nel mon /eldata/ dm, n, ma, mct, iel, nel integer nh1,nh2,nh3 mon /hdata/ nh1,nh2,nh3 integer ior,iow mon /iofile/ ior,iow integer ndf,ndm,nen1,nst,nneq mon /sdata/ ndf,ndm,nen1,nst,nneq real*8 hr integer mr mon / / hr(1),mr(1000) 圖 FEAPpv 單元公共塊 [32] Figure FEAPpv Element Common Blocks[32] FEAPpv中單元子程序的簡單結(jié)構(gòu)圖如圖 所示。 單元子程序中的各種信息，是通過啞元和公共變量進行傳遞的。 鋼筋混凝土基于切線剛度法的 二維 非線性分析 具體 流程圖 如下所示： 開始讀入模型信息讀入求解控制信息形成結(jié)構(gòu)結(jié)點位移列陣和結(jié)構(gòu)的等效結(jié)點荷載列陣已知單元結(jié)點位移和位移增量調(diào)用鋼筋混凝土單元形成單元剛度矩陣單元剛度集成結(jié)構(gòu)總體剛度 矩陣、單元結(jié)點恢復(fù)力集成結(jié)構(gòu)結(jié)點的恢復(fù)力由結(jié)構(gòu)的結(jié)點荷載和恢 解方程，求得結(jié)構(gòu)的結(jié)點位移增量收斂判斷復(fù)力求得殘差更新結(jié)構(gòu)的結(jié)點位移列陣 輸出計算結(jié)果結(jié)束計算是否達到最大迭代步數(shù)不收斂收斂是否達到目標(biāo)加載步是否是進入下一迭代步否進入下一荷載步更新結(jié)構(gòu)的結(jié)點位移列陣用結(jié)構(gòu)總剛和外荷載形成的結(jié)點力解方程，求得結(jié)構(gòu)的結(jié)點全量 位移用單元剛度集成結(jié)構(gòu)總體剛 度矩陣用已知的單元結(jié)點位移求解單元應(yīng)變并調(diào)用鋼筋混凝土本構(gòu)形成單元剛度矩陣形成結(jié)構(gòu)結(jié)點位移列陣和結(jié)構(gòu)的等效結(jié)點荷載列陣進入下一荷載步結(jié)束計算是否達到目標(biāo)加載步是 輸出計算結(jié)果收斂判斷收斂否讀入求解控制信息讀入模型信息開始不收斂是否達到最大迭代步數(shù)進入下一迭代步否是 鋼筋混凝土基于割線剛度法的非線性分析流程圖 本文所研究的鋼筋混凝土基于割線剛度的二維非線性有限元分析流程圖如下所示 可以看出 ，割線剛度法和切線剛度法算法流程的差別在于： 1. 割線剛度法 不需要計算單元恢復(fù)力以及結(jié)構(gòu)的殘余力 ，從而不必分配大型的數(shù)組用以存取恢復(fù)力和殘余力，節(jié)省了計算時間和儲存空間，提高了效 率； 2. 割線剛度法直接利用外荷載在結(jié)點上產(chǎn)生的力和集成的結(jié)構(gòu)總體剛度矩陣來求解結(jié)構(gòu)的全量位移，而非切線剛度法中通過殘余力求得的增量位移； 3. 割線剛度法以單元割線剛度元素作為收斂條件，相比于切線剛度法的收斂判斷準(zhǔn)則更為簡潔。圖 所示為整個宏命令的結(jié)構(gòu)圖： 鋼筋混凝土基于切線剛度法的非線性分析流程圖 FEAPpv 的源程序采用的是典型的基于切線剛度的有限元計算方法。 FEAPpv 的宏（ MACR）命令流包含與 6 個子程序，分別是 PMACR、PMACIO、 PMACR PMACR PMACR3 和 PMACR4，其中 PMACR4 是一個用戶自定義宏命令子程序，不使用時為空，且不會對整個程序運行有影響。下圖為控制程序的簡化流程圖。表 列出了一些主要數(shù)組的名稱、指針編號和它們各個維度的大小。 FEAPpv 程序可分為三個基本模塊： 1. 數(shù)據(jù)輸入模塊和先處理器 2. 求解模塊 3. 輸出模塊 以下是簡化的 FEAPpv 程序流程圖： 上圖中的程序輸入模塊將讀入輸入文件中的材料（ MATE）、幾何（ COOR， ELEM）、邊界條件（ BOUN）、荷載（ FORC）以及位移（ DISP）信息，用以建立所有的單元數(shù)組。所提供的材料模型適用于彈性、粘滯彈性、塑性以及熱傳導(dǎo)本構(gòu)方程，用戶也可以根據(jù)自己的需要添加另外的模型。 FEAPpv 程序包含一個有限的單元庫。 FEAPpv 程序自身包含了足夠多的命令來應(yīng)用于結(jié)構(gòu)力學(xué)、流體力學(xué)、熱傳導(dǎo)以及許多其他需要用微分方程來模擬問題的領(lǐng)域，靜態(tài)和動態(tài)問題都可以用該程序求解。一個問題的求解，是通過使用命令語句的形式來達到的，其中的求解算法則由用戶定義并編寫。在研究和應(yīng)用領(lǐng)域中，需要不斷對其進行改進來解決新的問題和滿足新的分析要求。 本文程 序選用的是第二種形式的收斂準(zhǔn)則。其中一種是以各個單元的割線剛度誤差的最大值作為判斷標(biāo)準(zhǔn)，即 ? ?1 1 2 21 2 1 2 1 2 1 2, , , , . . . . . . , , . . . . . . , ,i i n nc c c c c c c c EE E E E E E E E ? ??? ? ? ? ? ? ? ? ? ?, 其中，?是向量的最大范數(shù)， 12,iiccEE??分別是第 i 個單元的混凝土兩個主