基于有限单元柔度法的材料与几何双重非线性空间梁柱单元.pdf

http:/ 基于有限单元柔度法的材料与几何双重非线性 空间梁柱单元1 1陈滔1，黄宗明11重庆大学土木工程学院(400045)Email: 摘 要：摘 要：本文从有限单元柔度法的基本思想出发，基于完全拉格朗日格式(TL 格式)，建立了能考虑材料与几何双重非线性的一般化空间梁柱单元，适用于满足 EulerBernoulli 梁柱二阶分析理论假定的杆系结构非线性分析。另外，在梁柱单元截面分析中引入纤维模型，使其适用于解决钢筋混凝土这类复合材料结构的非线性分析问题。对钢筋混凝土双向偏心受压柱和钢筋混凝土框架试验结果的模拟分析表明，本文所提方法是正确、可靠的，能有效地分析框架柱的材料与几何双重非线性问题。关键词：关键词：钢筋混凝土柱 有限单元柔度法 梁柱单元 几何非线性 二阶效应 1引言 1引言 钢筋混凝土柱在变轴力及双向弯曲共同作用下的拟静力试验较为接近其在地震作用下的实际受力状况，是钢筋混凝土柱抗震性能研究中较为理想、接近实际的试验方案，但也是试验技术实现上最为困难的方案1,因而绝大多数情况下研究者们只能采用定轴力与双向弯曲或变轴力与单向弯曲组合的试验加载方案。因此，建立一种能有效考虑任意加载路径以及变轴力与双向弯曲共同作用影响的梁柱单元分析模型，模拟和替代部分结构试验，对于结构抗震性能的深入研究具有重要意义。建立这样的力学分析模型困难主要来自两个方面：其一是如何反映材料非线性所引起的变轴力与双向弯曲之间的强烈耦合作用；其二是如何反映框架柱在变形较大时的几何非线性效应，即工程应用中所惯称的二阶效应。目前对前一问题已有一些研究成果2-4，问题基本得到解决。这当中基于有限单元柔度法与纤维模型的梁柱单元优势突出，但在国内并未引起足够重视。对于后一问题，研究者们在钢筋混凝土双向偏心受压长柱的稳定性分析中也已有所探讨5-6，目前普遍采用基于有限单元刚度法的非线性空间梁柱单元。该模型以假定的单元位移插值函数为出发点，在分析中单元内部的位移场分布总是满足该假定的位移分布形式，而对单元内部截面力场的分布并未要求满足平衡条件。当单元实际位移场分布形式与假定位移场分布形式有较大差异时，计算效果便会受到影响，甚至不能反映问题的真实本质，因此往往通过加密计算网格（一根杆件离散为多个单元）的方法来逼近真实位移场。这种处理方式明显加大了求解问题的规模，牺牲了计算的效率，甚至会引起数值分析上的不稳定性及收敛性问题的发生。1本课题得到高等学校博士学科点专项科研基金（项目编号：20020611014）和国家自然科学基金（项目编号：59978055）资助 -1-http:/ 现有的基于有限单元柔度法的梁柱单元尚不能考虑框架柱在变形较大时的二阶效应。由于在材料非线性问题中，基于有限单元柔度法的梁柱单元较之基于有限单元刚度法的梁柱单元优势明显7，发展基于有限单元柔度法并能考虑材料与几何双重非线性的梁柱单元就成为非线性梁柱单元研究的重要方向。本文在已有研究工作的基础上7，从EulerBernoulli梁柱二阶分析理论基本假定出发，结合有限单元柔度法的基本思想，首次提出了能考虑材料与几何双重非线性的一般化空间梁柱单元。通过与文6钢筋混凝土双向偏心受压柱和文9钢筋混凝土门式框架试验结果的计算机模拟分析对比，对本文所提方法的正确与可靠性进行了验证。2梁柱单元分析 2梁柱单元分析 2.1 单元力与变形定义 2.1 单元力与变形定义 本文以下讨论基于 EulerBernoulli 梁柱二阶分析理论，即满足小应变、小转动以及平截面假定。应力应变度量方式采用工程应力与工程应变。有限单元描述格式采用完全拉格朗日格式(TL 格式)。分别用和eF ed表示梁柱单元的杆端力和杆端位移矢量。单元位移场矢量()()()()Texwxvxuxu=(1)(),(),(xwxvxu分别为单元杆轴上任一点处的zyx,方向的位移。若忽略杆件的剪切及扭转变形，则截面力矢量sxF)(和截面变形矢量可分别表示为 sxd)(TyzxsxMxMxNxF)()()()(=(2)(xNx为截面轴力，分别为绕截面轴的弯矩；)()(xMxMzy、zy、()()()()()()()()()TTyzxsdxxwddxxvddxxdwdxxdvdxxduxxxxd+=2222222121 (3)x为截面轴向应变，()()xxzy、分别为绕截面轴的曲率。zy、2.2 基本理论 2.2 基本理论 在柔度法中需首先给定单元截面力形函数矩阵，用以描述单元截面力场。对于图 1(a)所示的无刚体位移模式下的梁柱单元，可以证明根据单元变形后位形下的平衡关系及单元力学边界条件可得到式（4）以描述单元截面力场矢量。注意到该平衡方程与单元处于何种受力状态无关，即使单元处于材料与几何双重非线性阶段仍能成立（二阶理论假定范围内）。为表达清楚简化起见，在坐标系下讨论的有关各量均带上划线sxF)(zyxttt，而在坐标系（TL坐标系）下讨论的有关各量则不带上划线。zyx000-2-http:/()()efsFxNxF1=（4）图 1 空间梁柱单元力与变形(a)无刚体位移模式下的梁柱单元 (b)TL 坐标系下的梁柱单元(a)Beam-column element without rigid-body modes(b)Beam-column element in TL local coordinates Fig.1 Forces and associated deformations for spatial beam-column element 式中()1xNf为单元截面力形函数矩阵()()()=LxLxxwLxLxxvxNf100001000011 （5）F为无刚体位移模式下单元杆端力矢量 TyjyizjzixeMMMMNF=（6）与 eF对应的无刚体位移模式下单元杆端变形矢量 ed为 TyjyizjzieLd=（7）需说明的是，由于不考虑扭转变形与其余变形间的耦合作用，因而式（6）、（7）以及以下讨论均不计扭转变形及相应的扭矩,只在最终的刚度方程中直接按线弹性方式加入相应项。注意到式（3）为单元变形协调的微分形式，可以证明在考虑单元几何边界条件后，由式（3）可以得到单元变形协调的积分形式如下 ()()=LsTdexdxdxNd （8）式中 -3-http:/()()()=LxLxxwLxLxxvxNd100210012100001 （9）注意到在小应变和小转动假定下上述积分区间长度可不必区分变形前后的差异。L若采用柔度表述的一般化截面本构关系形式：()()()sssxFdxfxdd=（10）式中()sxf为截面切线柔度矩阵。它既可以由直接给出的截面力与变形关系确定，即由截面层次直接确定；也可以由更为细化的纤维模型确定，即由纤维材料应力应变关系集成确定7。联合式（4）、（8）和（10），经推导整理后可得到 eeeddFdf=(11)式中单元切线柔度矩阵 ef ()()()()()()()xdxAxdxxAxxdxNxfxNfzzLyyLfsTLde+=2121 (12)()()()()()()()()()()()()+=zjyjziyixyjyizjzixxffMxwMxwMxwMxwNxwMxvMxvMxvMxvNxvNxNxN000001 （13）()()()()()()=00000000000000000000yjyizjzixyMxwMxwMxwMxwNxwxA （14）()()()()()()=00000000000000000000yjyizjzixzMxvMxvMxvMxvNxvxA （15）对单元切线柔度矩阵求逆便可得到单元切线刚度矩阵 1=eefk （16）-4-http:/ 单元刚度矩阵 ek在上述推导中除采用梁柱二阶分析理论基本假定外，并未采用其他任何假定，因而在梁柱二阶分析范围内单元刚度矩阵是准确的，与单元所处的非线性状态无关，包括梁柱存在变截面情形。不过，在实际运用中，单元刚度矩阵的准确程度还依赖于截面模型化的方式，即截面切线刚度矩阵()sxk（()()1=ssxfxk）。2.3 单元刚度矩阵2.3 单元刚度矩阵 ek的计算 的计算 在确定单元刚度矩阵 ek时需要计算单元位移对单元杆端力的导数()()()()()()TyjyizjzixeMxvMxvMxvMxvNxvFxv=(17)()()()()()()TyjyizjzixeMxwMxwMxwMxwNxwFxw=(18)由于单元切线刚度矩阵计算时一般采用数值积分方式，如经典 Gauss 数值积分法，因而问题转化为计算()eFxv、()eFxw在若干积分控制点处的导数。已知个积分控制点处的截面曲率n()、iyx()izx()，可采用 Lagrange 多项式插值方法求得单元内任一点ni,1=x处曲率()、xy()xz为()()()()TnyyynyGLL211121=（19）()()()()TnzzznzGLL211121=（20）式中为范德蒙矩阵 G=121211211111nnnnniiinGLMMMMMLMMMMML （21）()LxniLxii=,1 应用单元位移与曲率的关系式（3）及单元几何边界条件可得()()()izGHv1=（22）()()()iyGHw1=（23）式中()()()()(+=+1322116121nnnLHL)（24）-5-http:/()()()()TnyyyiyL21=（25）()()()()TnzzzizL21=（26）由式（10）知()()()()()()()()()()()()()()()()()()=yzxssyzxsdMdMdNfffffffffFdfddddd333231232221131211 (27)将式（4）代入上式，整理后可得()()eTeyyFdFd=(28)()()eTezzFdFd=(29)式中()()()()()()TeyaaaaaF54321=()()()()()()TezbbbbbF54321=()()()()()()()()+=xxxxNwNwfNvNvffa3332311()()()()()()zixzixMwNfMvNfa+=333212()()()()()zjxzjxMwNfMvNfa+=33323()()()()()()+=yixyixMwNfMvNfa143332()()()()()+=yjxyjxMwNfMvNfa33325()()()()()()()()+=xxxxNwNwfNvNvffb2322211()()()()()()zixzixMwNfMvNfb+=232212()()()()()zjxzjixMwNfMvNfb+=23223 -6-http:/()()()()()()+=yixyixMwNfMvNfb142322()()()()()+=yjxyjxMwNfMvNfb23225 联合式(22)、(23)、(28)和(29)可得()()eTeFdFvvd=（30）()()eTeFdFwwd=（31）式中()()()()()TnxbbbGHNv111211L=()()()()()TnzibbbGHMv222211L=()()()()()TnzjbbbGHMv333211L=()()()()()TnyibbbGHMv444211L=()()()()()TnyjbbbGHMv555211L=()()()()()TnxaaaGHNw111211L=()()()()()TnziaaaGHMw222211L=()()()()()TnzjaaaGHMw333211L=()()()()()TnyiaaaGHMw444211L=()()()()()TnyjaaaGHMw555211L=利用上述导数公式计算各积分控制点的相应导数值，便可得到一系列线性方程组。此时已可采用迭代的方式求解方程组，由于积分控制点数一般取为 36 个因而方程组阶数较低，采用直接求解的方式精度与效率均较高，推荐采用。2.4 刚体位移的添加 2.4 刚体位移的添加 -7-http:/ 由于基于柔度法的梁柱单元推导时参考位形选取变形后的待求位形并在的无刚体位移模式下进行，即在 Euler 坐标系下建立单元力与变形的关系，因而有限单元方程建立采用TL 格式时，还需在两个坐标系之间建立相应变换关系，见图 1。在考虑小转动假定下，单元杆端位移矢量的变换为 eRBMeddTdd=（32）式中变换矩阵 RBMT=01010000010000010001010010001000001000001010001000010001LLLLLLLLLwwLvvLwwLvvTijijijijRBM 单元杆端力矢量的变换关系为 eFeTRBMeddTFdTFd+=（33）式中变换矩阵 FT=000000000000000000000000000000000000000100000100000010000010000000000000000000000000000000000000000000000000000100000100000010000010000000000000LNTxF 因此，在 TL 坐标系下建立单元杆端力与杆端位移之间的关系如下 eeeFdddk=（34）式中 FRBMeTRBMeTTkTk+=（35）3分析实例 3分析实例 算例 1：J.K.Kim 钢筋混凝土双向偏心受压柱6-8-http:/ J.K.Kim 完成了如图 2 所示一系列矩形或方形截面钢筋混凝土偏心受压柱试验工作，与以往早期研究者不同，除把问题定位于双向偏心外，而且量测了包括后峰值反应在内的完整非线性受力过程。为了保证所测量数据的可靠性，除采用先进的电液伺服加载系统外（INSTRON 8506），每组试件均采用完全相同的两根试件共同测定。共完成 8 组试件，分别为方形截面 3 组（荷载偏心角分别为），矩形截面 5 组（荷载偏心角ooo45,30,0分别为）。荷载偏心角示意见图 2。所有柱子两端加载边界条件都为铰结，且偏心距均为 40mm。ooooo90,60,45,30,0 Fig.2 Details of specimens and loading Fig.3 Axial forcelateral deflection curves for RS30 图 2 试件尺寸与加载 图 3 RS30 轴力挠度曲线 试验根据混凝土圆柱体测定了混凝土的抗压强度为MPa、初始弹性模量MPa 及劈拉强度 MPa。上述各值均为平均值，由 48 个圆柱体试验结果平均得到。纵筋采用3 变形钢筋（直径27=cf24300=cE4.3=spf525.9=dmm,屈服强度436=yfMPa）,箍筋采用光面钢筋（直径mm,屈服强度7.4=d165=yvf MPa），纵筋配筋率分别为矩形截面2.14%，方形截面 2.85%，箍筋间距为 100mm。混凝土保护层厚度为 18mm。本文给出了其中一组试件 RS30（矩形截面，荷载偏心角）的分析结果对比。其余试件组的模拟分析结果有类似分析精度，限于篇幅此处不一一列出，具体见文8。图 3中所示分析结果均只采用 1 个单元模拟。从图 3 中可看出，柔度法能很好地模拟柱子受力全过程，即使对于柱子处于不稳定受力阶段的模拟仍有足够高的分析精度。经分析发现若仅考虑材料非线性，柱子分析承载力值将高出实测值达 30，表明该柱已发生非弹性屈曲，注意到该柱长细比并不十分突出。o30算例 2：钢筋混凝土门式框架9 考虑到目前尚缺钢筋混凝土空间框架的有效试验数据，对钢筋混凝土框架的整体结构层次的验证选取了由原重庆建筑大学魏巍、刘毅等人所完成的两榀门式钢筋混凝土平面框架的试验结果。两榀试验框架均先行施加竖向静力荷载，后施加水平静力荷载至破坏。试验采用 -9-http:/ Fig 4 Details of test frames Fig 5 Computing schematic diagram 图 4 试验框架尺寸及配筋 图 5 计算简图 Fig 8 correlations on computing methods(KJ-1)Fig 9 correlations on computing methods(KJ-2)图 8 计算方法对比(KJ-1)图 9 计算方法对比(KJ-2)Fig 6 KJ-1 Lateral load-displacement Fig 7 KJ-2 Lateral load-displacement 图 6 KJ-1 荷载位移曲线 图 7 KJ-2 荷载位移曲线 MTS 电液伺服系统进行加载，获得了包括结构反应下降段在内的完整的非线性全过程受力-10-http:/ 信息，为检验本文方法在结构强非线性阶段的精度与适用性提供了可能。试验框架的配筋及尺寸见图 4。两榀试验框架编号分别为 KJ-1 和 KJ-2，KJ-1 最终破坏形式被设计为大偏压破坏（右柱），而 KJ-2 最终破坏形式被设计为小偏压破坏（右柱），通过加大右柱所受轴力实现。更为详尽的试验介绍见文献9。模拟分析中均采用 1 柱 1 个单元，横梁离散化为 3个单元，如图 5 所示两个集中荷载之间为 1 个单元。VRF从图 6 和图 7 模拟分析结果来看，柔度法分析结果与试验结果取得了相当一致的分析结果。表明所建议方法能很好反映梁柱构件包括二阶效应在内的非弹性受力行为，并且计算效率与精度都明显优于传统的基于刚度的两节点梁柱单元分析方法，这必将为分析大型复杂结构非弹性地震动力反应创造有利的计算条件。从刚度法所计算的结构承载力结果高出试验承载力近 50来看，主要是由于对材料非线性分析的不准确所造成，这一点可通过图 8 和图 9明显看出，这表明在钢筋混凝土框架非线性分析中材料非线性是问题的主要方面，这点在结构的正常使用阶段更是如此（大致在层间位移角达到 1200 之前）；但对于结构的破坏阶段重力二阶效应（几何非线性）所占比重逐渐加大，直至通常认为的倒塌界限时（层间位移角达到 150），仅考虑材料非线性可能会带来 50(KJ-1)和 100%(KJ-2)的分析误差。当然这不是一般性结论，仅是该具体试验结果的数值分析结论。但至少可以定性得出的结论是，在结构的强非线性阶段仅考虑结构的材料非线性的影响将可能会明显高估结构的承载力，同时也无法真正反映出结构已处于不稳定受力阶段的真实力学行为，从而掩盖了结构倒塌的真正原因。4.结论 4.结论 本文从有限单元柔度法的基本思想出发，基于完全拉格朗日格式(TL 格式)，建立了能考虑材料与几何双重非线性的一般化空间梁柱单元。所提出的梁柱单元可适用于满足 EulerBernoulli 梁柱二阶分析理论假定的杆系结构非线性分析，可采用任意的截面本构关系。算例表明，本文所提方法是正确、可靠的，能有效地分析框架柱的材料与几何双重非线性问题，为分析大型复杂结构非线性地震动力反应创造了有利的计算条件。参考文献 参考文献 1 Bousias S.N.,Verzeletti G.,Fardis M.N.,and Gutierrez E.Load-path effects in column biaxial bending with axial force J.Journal of Engineering Mechanics,ASCE,1995,121(5):596-605.2 Zeris C.A.,and Mahin S.A.Behavior of reinforced concrete structures subjected to biaxial excitation J.Journal of Structural Engineering,ASCE,1991,117(9):2657-2673.3 Spacone E.,Filippou F.C.,and Taucer F.F.Fiber beam-column model for non-linear analysis of R/C frames:part I.Formulation J.Earthquake Engineering and Structural Dynamics,1996,25:711-725.4 Spacone E.,Filippou F.C.,and Taucer F.F.Fiber beam-column model for non-linear analysis of R/C frames:part II.Applications J.Earthquake Engineering and Structural Dynamics,1996,25:727-742.5 叶英华等.双向偏心受压钢筋混凝土柱稳定性分析 J.建筑结构学报,1998,19(2):44-49.6 Kim J.K.,and Yang J.K.The behaviour of reinforced concrete columns subjected to axial force and biaxial bending J.Engineering Structures,2000,23:1518-1528.7 黄宗明，陈滔.基于有限单元柔度法和刚度法的非线性梁柱单元比较研究 J.工程力学,2003,20(5):24-31.8 陈滔.基于有限单元柔度法的钢筋混凝土框架三维非弹性地震反应分析 D.重庆大学.博士学位论-11-http:/ 文.2003.9 刘毅.钢筋混凝土框架柱二阶效应与稳定中若干问题的讨论 D.重庆建筑大学.硕士学位论文.1999.MATERIAL AND GEOMETRICALLY NONLINEAR SPATIAL BEAM-COLUMN ELEMENT BASED ON THE FINITE ELEMENT FLEXIBILITY METHOD CHEN Tao HUANG ZongMing College of Civil Engineering,Chongqing University,Chongqing,PRC,400045 Abstract:Based on the basic idea of the finite element flexibility method and total Lagrangian formulation(TL formulation),the general spatial beam-column element for dual nonlinearity of material and geometry is established in this paper.The proposed beam-column element is applicable to the general nonlinear analysis of framed structures that satisfy the basic assumptions of the second-order theory for Euler-Bernoulli beam-columns.Moreover,Fiber-model is introduced into the sectional analysis of the proposed beam-column element so that it can be applicable to the nonlinear analysis of composite material structures such as reinforced concrete.Numerical simulations of test examples including reinforced concrete two-way eccentrically loaded columns and reinforced concrete frames shows that the proposed formulation is exact and reliable，and can effectively solve the dual nonlinear problem of material and geometry for frame columns.Keywords:reinforced concrete;finite element flexibility method;beam-column element;geometrically nonlinear;second-order effect 陈滔：男。1973 年生。博士。主要研究方向是建筑结构抗震与钢筋混凝土结构基本性能。-12-