建筑结构在其使用寿命中,极可能遭遇极端荷载,如地震、爆炸、冲击,不可避免地造成建筑结构局部损伤或者破坏,从而极可能导致大规模结构破坏或连续倒塌的连锁反应,如Ronan Point公寓和世界贸易中心大楼事件。钢筋混凝土(RC)框架结构作为一种应用广泛的结构形式,大量科研人员进行了研究。已有的试验研究和理论分析[1―10]表明:钢筋混凝土框架结构倒塌破坏过程由梁机制演变至空腹桁架机制,最后发展至悬链线机制直至破坏。悬链线机制[11]是在结构整体中的水平构件失去抗弯能力后,通过轴力和很大的挠度形成的力矩来抵抗外荷载所产生的弯矩。悬链机制会使结构产生有助于抵抗连续倒塌的附加承载能力,对结构抗连续倒塌能力至关重要。悬链机制处于几何大变形和材料非线性下降段的状态下,需要同时考虑材料非线性和几何非线性,因此对数值分析模型提出了更高的要求。
纤维单元(fiber element)通常被采用模拟RC框架结构的非线性分析,其本质是将构件截面分割划分成许多纵向纤维,每根纤维都处于单轴应力状态,它的优点是可以采用成熟的材料单轴本构关系,并且能够考虑轴力和弯矩的耦合作用。纤维单元有基于力插值和基于位移插值的纤维单元之分[12]。在基于位移插值的单元形函数中,曲率假定为线性分布,对于弹性阶段可以得到精确解,而当构件进入非线性阶段后曲率不再是线性分布,所以非线性分析时一般需将单元进行细分,以减小曲率近似线性分布造成的误差。文献[13]的研究结果表明,一个构件至少要划分为五个单元才能达到较高的精度。Yu和Tan[14]利用Engineer’s Studio软件沿着梁单元利用纤维单元进行建模,每个纤维单元对应一至两个高斯积分点。几何大变形情况下,为解决收敛性问题,将梁划分成多个梁单元。采用Maekawa等[15]研究小组开发的混凝土和钢筋本构模型,此本构模型考虑了钢混之间的粘滑关系。数值结果与试验结果相比,提出的节点模型能够准确预测荷载-变形的趋势。梁益和陆新征等[16]采用MSC Marc有限元软件和清华大学开发的钢筋混凝土纤维梁分析模型THUFIBER[17]作为数值分析平台,可满足几何大变形下结构倒塌求解计算,THUFIBER纤维模型可以对RC杆系构件在复杂受力条件下进行准确模拟。于晓辉等[18]采用 OpenSees软件作为分析平台,为考虑压拱效应和悬链效应,采用 OpenSees软件中基于位移的非线性梁柱单元以及Corotational坐标转化方法来考虑梁在大变形情况下可能引起的附加二阶效应。可见,基于位移的纤维单元可以较好地同时考虑材料非线性和几何非线性。基于力插值的纤维单元(force-based fiber element)采用力插值函数,在不考虑几何非线性、刚体位移及无单元荷载作用的前提下总能保证单元内部任意截面的力场分布均匀且严格满足平衡条件,优点在于不受线性曲率假定分布的限制,所以不需要再细分单元。基于力插值纤维单元的优势是使用一个单元就可模拟一个构件,计算效率大幅提高。基于力插值单元最先由 Mahasuverachai[19]引入,随后Brancaleoni等[20]研究了截面的弯矩-曲率关系,在此基础上Spacone等[21]建议了一种基于力插值单元的状态确定过程,之后被延伸推广,得以广泛使用。Kaba和 Mahin[22]第一次将纤维截面和基于力插值单元结合在一起,用于解决单轴弯曲的非线性问题,随后Zeris和 Mahin[23―24]改进了这一模型,并将其扩展到双轴弯曲。以上研究都是只关注于材料非线性的问题,为了解决基于力插值纤维单元求解几何非线性的问题,Neuenhofer等[25]对单元构造的力插值函数进行修改,使其弯矩不再采用线性插值,而是跟竖向位移有关的函数。为了得到位移的显式表达,构造了基于曲率的位移插值形成(curvature-based displacement interpolation)。然而这一单元只能用于线性材料、几何非线性问题,另外不能解决大转动问题。
目前几何非线性求解问题在固体有限元中按照坐标参考系描述不同分为三种:总体拉格朗日法(total Lagrangian)、更新拉格朗日法(updated Lagrangian)和共旋坐标法(corotational)。总体拉格朗日法始终将初始未变形状态作为参考系,更新拉格朗日将上一步计算的变形状态作为参考系。共旋坐标法是对上面两种方法的一种改进,在每个单元上附上一个跟随单元平动和转动的新坐标系统,将刚体的平动和转动精确扣除,更适用于大转动问题[26]。为了解决基于力插值的纤维单元同时处理材料非线性和几何非线性的问题,本文采用基于共旋坐标法,提出了一种基于共旋坐标法的力插值纤维单元,给出了二维单元形成原理及非线性求解过程,对两个平面算例进行模拟验证。分析结果表明,基于共旋坐标法的力插值纤维单元能够较准确地模拟RC框架结构连续倒塌的悬链线机制。
对于二维平面单元的坐标系如图1所示。对于整体坐标系(X, Y)而言,单元有两个节点I和J,每个节点有3个自由度(2个平动、1个转动),整个单元有6个自由度。整体坐标系下的单元节点力和位移向量分别表示为P和D:
共旋坐标法需要将刚体的平动和转动精确扣除。扣除刚体位移的局部坐标系(x, y)中有3个自由度,包括1个平动、2个转动。
图1 二维单元整体和局部坐标系
Fig.1 Global and local coordinate systems of 2D element
基于共旋坐标法的力插值纤维单元在局部坐标系的变形体内采用平截面假定,沿截面方向划分为纤维,考虑纤维的单轴材料非线性本构关系,然后加上刚体位移,从局部坐标系到整体坐标系的转换中采用共旋坐标法以考虑几何非线性大转动。
扣除刚体位移的局部坐标系(x, y)中有3个自由度(1个平动,2个转动),如图2所示。局部坐标系下的单元节点力和位移向量分别表示为Q和S:
式中: JN 代表节点J处的轴力; IM 、 JM 代表节点I、J处的弯矩;Ju代表节点J处的轴向位移;Iθ、Jθ代表节点I、J处的转角。
单元x截面处的截面力向量可表示为:
式中: N (x)代表截面 x处的轴力; M z (x)代表截面x处绕z轴的弯矩。
基于力插值纤维单元构造的是力插值函数,在不考虑刚体位移,以及无单元荷载作用的前提下认为弯矩为线性插值函数,轴力为常量插值函数:
式中, ()x b 代表力插值函数矩阵。
图2 二维单元局部坐标系下纤维化
Fig.2 Fiber formation in local coordinate systems of 2D element
截面变形向量与力向量的关系表示为:
式中: ()x d 代表截面变形向量; ()x f 代表截面的柔度矩阵,可以由截面刚度矩阵k(x)求逆得到。
假设在平衡状态下,单元两端部作用一广义虚力向量 δQ,此时x截面虚功为:
则整个单元在虚力向量 δQ作用下的虚功可通过沿着单元长度积分可得:
得到:
式中:F为单元的柔度矩阵;K为单元的刚度矩阵。
为了得到截面刚度矩阵k(x),构件截面沿其纵向分割为有限根纤维,每种材料有各自的单轴非线性本构关系。确定每根纤维的轴向应变,这样由本构关系确定其应力,积分得到轴力及弯矩。由此可见,局部坐标系下的变形体非线性只考虑纤维的材料非线性。
单元x截面处的截面变形向量可表示如下:
式中:xε为轴向应变;zφ为绕截面z轴的曲率。
假设x截面有n个纤维,记i处纤维的坐标为(yi, zi),切线模量为E t i,由于每根纤维均处于单轴应力-应变状态,由平截面假定(如图2所示),由截面几何关系(二维情况下)可得 i处纤维的轴向应变可表示为:
得到:
在x截面处有如下关系式成立,即:
使用共旋坐标法有一突出优点,局部坐标系中的单元形成与局部到整体坐标系转换完全无关,因此可以在坐标转换中考虑几何非线性,在局部坐标系的变形体中考虑材料非线性,从而实现单元中材料非线性和几何非线性完全分开考虑。
图 3所示为未变形时的局部坐标系 ˆˆ(,)x y与整体坐标系(,)X Y示意图。首先,单位向量1ˆe为:
式中: IJ X 为节点I和J在整体坐标系的有向向量;L为未变形时的单元长度。
单位向量2ˆe为:
图3 未变形时的局部与整体坐标转换
Fig.3 Transformation between local and global coordinate systems in undeformed configuration
图 4所示为变形后的局部坐标系(x,y)、未变形时的局部坐标系 (xˆ,yˆ)以及整体坐标系(X,Y)示意图。整体坐标系下的单元节点位移向量 D可以写为:
图4 变形后的局部与整体坐标转换
Fig.4 Transformation between local and global coordinate systems in deformed configuration
式中:I U、J U为节点I和J水平位移向量;Iγ、Jγ为节点I和J转角位移。
单位向量 1e为:
式中: IJ U 为变形后节点I和J在整体坐标系的有向向量;l为变形后的单元长度。
单位向量 2e为:
根据图 5所示,局部坐标系的轴向位移可表示为:
转动位移可表示为:
式中,转角α可按式(26)计算得到:
当轴向变形在初始比较小的时候,式(30)的代数解形式计算精度会受影响,因此写成向量形式表达:
式(31)和式(32)在计算转动位移时采用了代数解的表达,因此当刚体大转动而变形体转动小的情况下,ˆα与α数值过于接近,计算精度会受影响,式(31)可重新表示为:
式中:
同样,式(31)也可重新表示为:
可得:
同理,可计算Jθ:
式中:
图5 位移向量的局部与整体坐标转换
Fig.5 Transformation of displacements between local and global coordinate systems
增量位移向量可以由局部坐标系位移向量到整体坐标系位移向量求导得到:
式中,δS、δD分别为局部坐标系和整体坐标系下的增量位移向量。
为了得到需要重新计算。
由式(28)和式(26)可得:
其中,根据式(25)和式(27),有:
由式(26)微分得到1δe:
式(26)两边同时微分,并代入式(29)得到:
式(42)两边同时乘以可以得到:
式(41)代入式(43)中,并由 1e、 2e正交性可得到:
其中,根据式(25)和式(27):
由式(30)和式(31)微分,并将式(39)、式(44)代入:
其中:
根据图6可以精确得到力向量由局部坐标系到整体坐标系的转换关系:
其中:
式中,V为单元内剪力。
把式(48)写成矩阵形式表示:
其中:
可见,力向量转换矩阵与增量位移向量转换矩阵相同。
图6 力向量的局部与整体坐标转换
Fig.6 Transformation of forces between local and global coordinate systems
局部坐标系下单元刚度矩阵增量形式可写为:
整体坐标系下单元刚度矩阵ˆK可由式(46)和式(52)代入式(50)求得:
式中, G K 称为几何刚度矩阵:
式中,“:”为缩略词,代表以下关系式:
式中, r t代表转换矩阵T的第r行。
将式(40)、式(42)、式(44)代入式(47)中,可得到:
同理,可以得到和
因此,式(55)可以写为:
最终几何刚度矩阵 G K 可以写为:
式中:
式中,cos cα=,sin sα=。
将力插值纤维模型与共旋坐标法理论结合在一起即能够充分发挥纤维截面模型精细化优势保证计算的精度,继承力插值单元所具有的独特特点,即无论结构构件处于何种状态下(即使下降软化段),单元内部力的分布也是正确的,同时采用共旋坐标法以考虑几何非线性,将刚体的平动和转动精确扣除。基于共旋坐标法的力插值纤维单元有一突出优点,局部坐标系中的单元形成与局部到整体坐标系转换完全无关,因此可以在坐标转换中考虑几何非线性,在局部坐标系的变形体中考虑材料非线性,从而实现单元中材料非线性和几何非线性完全分开考虑。基于共旋坐标法的力插值纤维单元非线性迭代计算过程可划分为三个层次:截面状态的确定、局部坐标系下单元状态确定以及整体坐标系下单元状态确定。具体如下:
1) 首先由每根纤维的材料本构关系曲线,根据平截面假定得到截面位移、力、刚度矩阵,确定截面状态;
2) 然后在扣除刚体位移,只有变形体的局部坐标系单元内通过构造的力插值函数关系 ()b x,根据1.1节式(6)、式(8)、式(13)和式(14),得到局部坐标系下的单元力、位移以及刚度矩阵,确定局部坐标系下的单元状态。
3) 最后加上刚体位移,通过转换矩阵T,根据1.2节~1.5节公式,得到整体坐标系下的单元力、位移以及刚度矩阵,确定整体坐标系下的单元状态。
局部坐标系下单元状态确定以及整体坐标系下单元状态确定都需要迭代,因此整个计算流程需要2次迭代,采用New ton-Raphson迭代法,其详细非线性迭代计算流程见图7。
图7 非线性迭代计算流程图
Fig.7 Schematic statement of nonlinear iteration algorithm
上述给出了考虑悬链机制RC框架结构连续倒塌分析模型:基于共旋坐标法的力插值纤维单元的原理及形成,通过两个算理以验证此单元在材料非线性及几何大变形下是否具有更好的优势。
本文选取 2013年连续倒塌试验分析竞赛的带楼板梁子结构试验进行数值模拟,如图8所示。梁板截面尺寸、截面配筋以及钢筋混凝土材料特性等试验参数见参考文献[27]。分析中使用基于共旋坐标法的力插值纤维单元,分别就以下三种情况进行分析:
1) 只考虑材料非线性;
2) 考虑材料非线性及P-Δ效应;
3) 考虑材料非线性及几何大变形。
图8 带楼板梁子结构连续倒塌试验
Fig.8 The beam-slab substructure progressive collapse test
混凝土采用Concrete02本构模型,考虑混凝土受拉影响。采用Mander模型[28]考虑箍筋约束对核心混凝土抗压强度有利影响。钢筋采用 Hysteretic(三线式)单轴本构模型,并考虑下降段。截面采用纤维模型,在截面上划分若干混凝土纤维和钢筋纤维,梁板子结构划分成2个单元,如图9所示。
图9 连续倒塌试验分析模型
Fig.9 Modeling approach of progressive collapse test
根据分析结果的荷载-位移曲线,如图10所示,仅考虑材料非线性的情况下,试件梁端进入全塑性阶段之后,试件承载力的数值结果显著降低,可以准确模拟梁机制阶段,但是完全不能模拟悬链线机制。考虑材料非线性及P-Δ效应的情况下,梁机制承载力数值结果较大,因为梁在受竖向力作用产生小挠度情况下,由于轴力作用,会产生压拱作用,采用 P-Δ几何坐标变换考虑挠曲二阶效应,但由于梁机制失效之后,梁全段受拉,采用 P-Δ几何坐标变换会高估试件的梁机制承载力,没有出现明显的承载力下降段及悬链线机制。本文提出的基于共旋坐标法的力插值纤维单元模拟的数值结果可以看出,明显优于前两种情况,能比较精确试件的连续倒塌全过程的梁机制和悬链线机制阶段力学行为。同时分析结果也表明,梁机制阶段主要是材料非线性起控制作用,悬链线机制阶段主要是几何非线性大转动起控制作用。
图10 连续倒塌试验分析结果
Fig.10 Analytical results of progressive collapse test
为了进一步验证本文开发的基于共旋坐标法的力插值纤维单元合理性,本文再次选取了易伟建等[1]为了研究 RC框架倒塌性能而进行的平面框架结构试验进行数值模拟,层高、层数、梁柱截面尺寸、截面配筋及钢筋混凝土材料特性等试验参数见参考文献[1]。采用同算例 1一样的建模方法。
分析结果和试验结果对比如图 11和图 12所示,数值模拟结果与试验结果较吻合,再次证明了本文开发的基于共旋坐标法的力插值纤维单元在模拟RC框架连续倒塌全过程的力学行为具有较高的精度。
图11 试验和模拟的最终破坏状态对比
Fig.11 Comparisons of numerical and experimental failure modes at final damage
图12 试验和分析结果对比
Fig.12 Comparisons between experimental data and analytical results
本文提出的考虑悬链机制RC框架结构连续倒塌分析模型:基于共旋坐标法的力插值纤维单元,解决了基于力插值的纤维单元同时处理材料非线性和几何非线性的问题。基于共旋坐标法的力插值纤维单元的优点是局部坐标系中的单元形成与局部到整体坐标系转换完全无关,因此可以在坐标转换中考虑几何非线性,在局部坐标系的变形体中考虑材料非线性,从而实现单元中材料非线性和几何非线性完全分开考虑。
本文给出了基于共旋坐标法的力插值纤维二维单元形成原理及非线性求解过程。通过两个实例分析表明,基于共旋坐标法的力插值纤维单元能够较准确的模拟 RC框架结构连续倒塌的梁机制和悬链线机制,梁机制阶段主要是材料非线性起控制作用,悬链线机制阶段主要是几何非线性起控制作用。
[1] 易伟建, 何庆锋, 肖岩. 钢筋混凝土框架结构抗倒塌性能的试验研究[J]. 建筑结构学报, 2007, 28(5): 104―109.Yi Weijian, He Qingfeng, Xiao Yan. Collapse performance of RC frame structure [J]. Journal of Building Structures, 2007, 28(5): 104―109. (in Chinese)
[2] 苏幼坡, 王兴国, 宋晓胜, 等. 钢筋混凝土框架梁拱效应的试验研究[J]. 西安建筑科技大学学报(自然科学版), 2009, 41(4): 477―484.Su Youpo, Wang Xingguo, Song Xiaosheng, et al.Experimental investigation on the arching action in reinforced concrete frame-beams [J]. Journal of Xi’an University of Architecture & Technology, 2009, 41(4):477―484. (in Chinese)
[3] 何庆锋, 易伟建. 考虑悬索作用钢筋混凝土梁柱子结构抗倒塌性能试验研究[J]. 土木工程学报, 2011,44(4): 52―59.He Qingfeng, Yi Weijian. Experimental study of the collapse-resistant behavior of RC beam-column sub-structures considering catenary action [J]. China Civil Engineering Journal, 2011, 44(4): 52―59. (in Chinese)
[4] 王浩, 李易, 陆新征, 等. 单层钢筋混凝土框架结构水平向连续倒塌试验研究[J]. 建筑结构学报, 2016,37(10): 65―72.Wang Hao, Li Yi, Lu Xinzhen, et al. Experimental investigation on horizontal progressive collapse of one-story reinforced concrete frame [J]. Journal of Building Structures, 2016, 37(10): 65―72. (in Chinese)[5] Sasani M, Werner A, Kazem i A. Bar fracture modeling in progressive collapse analysis of reinforced concrete structures [J]. Engineering Structures, 2011, 33(2): 401―409.
[6] Qian K, Li B, Ma J X. Load-carrying mechanism to resist progressive collapse of RC buildings [J]. Journal of Structural Engineering, 2014, 141(2): 04014107-1―04014107-14.
[7] 孟宝, 钟炜辉, 郝际平. 不同跨度比下栓焊刚性连接梁柱子结构抗倒塌性能试验研究[J]. 工程力学, 2018,35(1): 79―87.Meng Bao, Zhong Weihui, Hao Jiping. Experimental study on anti-collapse performance for beam-column assemblies w ith bolt and weld rigid connection based on different span ratio [J]. Engineering Mechanics, 2018,35(1): 79―87. (in Chinese)
[8] 王俊杰, 王伟, 孙昕. 压型钢板组合梁中柱子结构的抗连续倒塌试验[J]. 工程力学, 2017, 34(增刊 1):149―153, 178.Wang Junjie, Wang Wei, Sun Xin. Experimental behavior of composite beam-column joints w ith steel profiled decking in a middle-column-removal scenario [J].Engineering Mechanics, 2017, 34(Suppl1): 149―153,178. (in Chinese)
[9] Yu J, Luo L, Li Y. Numerical study of progressive collapse resistance of RC beam-slab substructures under perimeter column removal scenarios [J]. Engineering Structures, 2018, 159: 14―27.
[10] Ren P, Li Y, Lu X, et al. Experimental investigation of progressive collapse resistance of one-way reinforced concrete beam-slab substructures under a m iddlecolumn-removal scenario [J]. Engineering Structures,2016, 118: 28―40.
[11] 何政, 黄国辉. 框架结构悬链线效应研究新进展[J].力学进展, 2012, 42(5): 547―561.He Zheng, Huang Guohui. Progress in studies of catenary action in frame structures [J]. Advances in Mechanics, 2012, 42(5): 547―561. (in Chinese)
[12] Scott M H, Fenves G L. Plastic hinge integration methods for force-based beam-column elements [J].Journal of Structural Engineering, ASCE, 2006, 132(2):244―252.
[13] 杜轲, 孙景江, 许卫晓. 纤维模型中单元、截面及纤维划分问题研究[J]. 地震工程与工程振动, 2012, 32(10):39―46.Du Ke, Sun Jingjiang, Xu Weixiao. The division of element, section and fiber in fiber mode [J]. Earthquake Engineering and Engineering Vibration, 2012, 32 (10):39―46. (in Chinese)
[14] Yu J, Tan K H. Experimental and numerical investigation on progressive collapse resistance of reinforced concrete beam column sub-assemblages [J]. Engineering Structures, 2013, 55(4): 90―106.
[15] Maekawa K, Okamura H, Pimanmas A. Nonlinear mechanics of reinforced concrete [M]. London: Spon Press, 2003.
[16] 梁益, 陆新征, 李易, 等. 3层RC框架的抗连续倒塌设计[J]. 解放军理工大学自然科学版, 2007, 8(6): 659―664.Liang Yi, Lu Xinzheng, Li Yi, et al. Design method to resist progressive collapse for a three story RC frame [J].Journal of PLA University of Science and Technology(Natural Science Edition), 2007, 8(6): 659―664. (in Chinese)
[17] 叶列平, 陆新征, 马千里, 等. 混凝土结构抗震非线性分析模型、方法及算例[J]. 工程力学, 2006, 23(s2):131―140.Ye Lieping, Lu Xinzheng, Ma Qianli, et al. Nonlinear analytical models, methods and examples for concrete structures subject to earthquake loading [J].Engineering Mechanics, 2006, 23(s2): 131―140. (in Chinese)
[18] 于晓辉, 钱凯, 吕大刚. 考虑悬链线效应的钢筋混凝土框架结构抗连续倒塌能力分析 [J]. 建筑结构学报,2017, 38(4): 28―34.Yu Xiaohui, Qian Kai, Lü Dagang. Progressive collapse capacity analysis of reinforced concrete frame structures considering catenary action [J]. Journal of Building Structures, 2017, 38(4): 28―34. (in Chinese)
[19] Mahasuverachai M. Inelastic analysis of piping and tubuar structures [R]. Earthquake Engineering Research Center, University of California, Berkeley. 1982.
[20] Brancaleoni F, Ciampi V, Di Antonio R. Rate-type models for non linear hysteretic structural behavior [M].Palermo, Italy: EUROMECH Colloquium, 1983.
[21] Spacone Enrico, Ciampi V, Filippou F C. A beam element for seism ic damage analysis [R]. Earthquake Engineering Research Center, University of California,Berkeley, August 1992.
[22] Kaba S, Mahin S A. Refined modeling of reinforced concrete columns for seism ic analysis [R]. Earthquake Engineering Research Center, University of California,Berkeley, 1984.
[23] Zeris C A, Mahin S A. Analysis of reinforced concrete beam-columns under uniaxial excitation [J]. Journal of Structural Engineering, ASCE, 1988, 114(4): 804―820.[24] Zeris C A, Mahin S A. Behavior of reinforced concrete structures subjected to biaxial excitation [J]. Journal of Structural Engineering, ASCE, 1991, 117(9): 2657―2673.
[25] Neuenhofer A, Filippou F C. Geometrically nonlinear flexibility-based frame finite element [J]. Journal of Structural Engineering, 1998, 124(6): 704―711.
[26] 魏科. 基于 CR列式的斜拉桥几何非线性分析[D]. 长安大学, 2009.Wei Ke. Geometrical nonlinearity analysis of cable-stayed bridges based on CR formulation [D].Chang’an University, 2009. (in Chinese)
[27] 中国建筑学会抗震防灾分会建筑结构抗倒塌专业委员会, 清华大学. 2013年连续倒塌试验分析竞赛[EB].http://www.collapse-prevention.net/show.asp?ID=20&ad ID=3. 2013.Collapse Prevention Committee, Architectural Society of China, Tsinghua University. 2013 Progressive Collapse Test Competition [EB]. http://www.collapse-prevention.net/show. asp?ID=20&adID=3. 2013. (in Chinese)
[28] Mander J B, Priestley M J N, Park R. Theoretical stress-strain model for confined concrete [J]. Journal of Structural Engineering, 1988, 114(8): 1804―1826.
A PROGRESSIVE COLLAPSE ANALYTICAL MODEL OF RC FRAME STRUCTURES BASED ON COROTATIONAL FORMULATION FOR FORCE-BASED FIBER ELEMENTS