基于刚体准则的平面框架结构静力非线性分析方法与流程

本发明涉及用于结构数值分析领域，具体涉及一种基于刚体准则的平面框架结构静力非线性分析方法。

背景技术：

钢框架结构极限状态与承载力分析涉及由于构件截面屈服导致的材料非线性与几何非线性的耦合问题，其关键在于非线性单元的构建以及迭代分析方法的确定。常用非线性单元包括以w.f.chen等提出的考虑二阶效应的梁柱单元与纤维元，迭代分析方法则多为基于更新拉格朗日格式(updatedlagrange，简称ul)的增量迭代法。

w.f.chen等提出的考虑二阶效应的梁柱理论与塑性铰理论是构建常用非线性单元的基础。由于轴力与弯矩的耦合作用，通常难以得到杆件的精确屈曲函数，而多采用高阶多项式建立近似梁柱非线性单元；塑性铰模型可以方便有效地描述材料屈服导致的构件截面刚度退化，chan和zhou等使用高阶多项式单元结合塑性铰模型，考虑了屈服截面的转动效应；iu和bradford提出了改进塑性铰模型，引入了轴力屈服弹簧，考虑了屈服截面转动与拉伸效应；vogel采用塑性区域模型分析了截面的屈服过程以及屈服导致的转动与拉伸效应。纤维模型将构件视为若干纤维单元的组合，可以完整考虑构件截面的应力分布及其变化过程，是分析框架结构二阶效应与非弹性效应的精度较高的方法，但纤维模型对于大型复杂结构而言，单元划分数较多，计算成本过大。综上，梁柱理论采用了数学意义上的近似屈曲函数，需要依靠插值函数的阶次或细分单元以提高精度，而纤维模型计算效率难以满足复杂工程需求。因此，两类模型的计算效率都难以满足工程需求。

结构非线性分析通常采用基于ul格式的增量迭代法，即：将荷载分为一系列增量荷载的叠加，在每个增量步内将非线性方程线性化，采用迭代方法将线性化带来的误差减小到可接受的程度，再进行下一增量步的计算。循环反复从而得到结构的受力变形全过程。但对于强非线性问题，若单元结点力计算阶段用到几何刚度矩阵，则几何刚度矩阵是否合理就决定了计算是否准确，故需要推导准确的几何刚度矩阵。前述方法所用结构非线性刚度矩阵尤其是几何刚度矩阵的近似性，所产生的误差往往无法通过迭代得到修正，会导致计算结果的严重失真甚至发散，从而在分析实际框架结构中几何与材料非线性耦合的复杂非线性问题时，无法得到精确的极限承载力与结构变形过程。

技术实现要素：

为解决现有技术中存在的技术问题，本发明提供一种基于刚体准则的平面框架结构静力非线性分析方法，能够精确地得到结构加载过程的荷载位移曲线和极限荷载，所用单元形式简单，物理意义明确，毋需划分大量单元，计算精度与效率高。

具体方案如下：

一种基于刚体准则的平面框架结构静力非线性分析方法，其关键在于包括

以下步骤：

步骤1、建立集中塑性铰弹簧平面梁单元模型；

步骤2、建立满足刚体准则的上述模型对应的单元刚度矩阵；

步骤3、基于更新拉格朗日格式与刚体准则，建立求解方程的增量迭代过程；

步骤4、迭代求解增量方程，得到结构加载过程的荷载位移曲线与极限荷载。

更进一步的，建立集中塑性铰弹簧平面梁单元模型，设框架梁单元的屈服只发生于零长度的塑性铰内，即不考虑材料塑性沿构件长度方向的发展，屈服集中在一个截面内，而塑性铰之间杆件仍为完全弹性。

更进一步的，利用轴力与弯矩耦合的截面精细初始屈服与全屈服曲线，建立框架梁单元截面逐步屈服塑性铰弹簧刚度模型，塑性铰弹簧的刚度系数为

其中，ei/l为梁的弹性弯曲刚度，初始屈服函数为完全屈服函数为fp为截面的轴向承载力，mp为塑性极限弯矩，fx和m分别为单元的轴力和弯矩。

更进一步的，采用以下单元刚度矩阵：

[k]＝[kg]+[kep]

其中，[kg]为单元几何刚度矩阵，[kep]为单元弹塑性刚度矩阵。

更进一步的，利用刚体准则，即：对处于初始平衡状态的刚体单元，当发生刚体移动和转动时，由于单元没有发生变形，单元在原结点力作用下应保持平衡，原单元结点力的大小不变，力的方向随单元转动而转动，和虚功原理建立的单元几何刚度矩阵如下：

当单元发生刚体位移{u}r时，单元虚应变能为零，几何刚度矩阵不会产生结点力增量，即：[kg]{u}r＝0，其中{u}r^t＝<00θr0lθrθr>。

更进一步的，利用截面逐步屈服塑性铰框架梁单元结点弯矩的增量平衡方程，建立单元弹塑性刚度矩阵，形式如下

其中，

k36＝2eizsasbl/(12e²iz²+4eizl(sa+sb)+sasbl²)；

k66＝eizsb(12eiz+4sal)/(12e²iz²+4eizl(sa+sb)+sasbl²)；

其中，e为弹性模量，a为单元截面面积，l为单元长度，iz为截面惯性矩，sa为单元a端的弹簧刚度，sb为单元b端的弹簧刚度。

更进一步的，增量迭代过程分为预测阶段、单元结点力计算阶段和误差检查阶段，预测阶段是从结构增量平衡方程中解得结构位移增量即求解其中，表示第i增量步内第j-1迭代步所形成的结构总体刚度矩阵，由单元刚度矩阵一一对应集成；为第j迭代步待计算的结构位移增量；为第j迭代步的荷载增量，为节点不平衡力，数值为第j-1迭代步的荷载向量与结构内力向量的差值；单元结点力计算阶段是由单元结点位移增量计算出单元结点力误差检查阶段首先将各单元结点力组装得到整体坐标下的单元内力并与荷载比较得到结点不平衡力若结点不平衡力满足收敛要求，则进行下一增量步计算，否则重复本增量步内的预测阶段和结点力计算阶段直至满足收敛要求。

更进一步的，在单元结点力计算阶段，单元结点位移增量由单元的刚体转动和单元的自然变形两部分组成，单元结点力由下式计算：

其中，为初始状态下的单元结点力列阵，将其乘以几何刚度矩阵，确定单元大转动的作用，{f}是自然变形产生的单元结点力增量。

与现有结构非线性分析技术相比，本发明方法由于将刚体准则与考虑截面逐步屈服的塑性铰弹簧模型相结合，具有以下显著优点：

1、单元几何刚度矩阵形式简洁；

2、所采用的刚度矩阵均为线性刚度矩阵，形式简单；

3、增量方程推导简洁，物理意义明确；

4、适于柔性平面钢框架结构的几何非线性和材料非线性耦合问题，划分单元数少，能极大提高实际工程结构极限承载力分析的计算精度与效率。

附图说明

本发明的附图说明如下：

图1集中塑性铰弹簧平面梁单元模型示意图；

图2弯矩和轴力共同作用下的修正初始和全屈服曲线图；

图3平面梁单元发生刚体转动示意图；

图4增量-迭代过程的三阶段示意图；

图5平面框架梁单元的变形示意图；

图6为有初始倾角的单层单跨平面钢框架结构示意图；

图7为包括本发明在内的5种分析方法得出的荷载位移曲线图；

图8为多层多跨平面钢框架结构示意图；

图9为包括本发明在内的5种分析方法得出的荷载位移曲线图；

图10为本发明的方法流程图。

具体实施方式

下面结合实施例和附图对本发明作进一步说明。

本发明包括以下步骤：

步骤1、建立集中塑性铰弹簧平面梁单元模型；

如图1所示的集中塑性铰弹簧平面梁单元，在平面弹性梁单元两端引入塑性铰弹簧，设框架梁单元的屈服只发生于零长度的塑性铰内，即不考虑材料塑性沿构件长度方向的发展，屈服集中在一个截面内，而塑性铰之间杆件仍为完全弹性。塑性铰弹簧刚度表示为：

其中ei/l为梁的弹性弯曲刚度，φy(m,fx)为初始屈服函数，φp(m,fx)为完全屈服函数，fx和m分别为单元的轴力和弯矩。

如图2，采用轴力与弯矩耦合的截面精细初始屈服与全屈服曲线，初始屈服函数和完全屈服函数分别为：

式中，fp为截面的轴向承载力，mp为塑性极限弯矩，fx和m分别为单元的轴力和弯矩。

步骤2、建立满足刚体准则的上述模型对应的单元刚度矩阵；

步骤1)、单元刚度矩阵表达式的确定

图1所示的弹塑性梁单元的单元刚度矩阵[k]可以表示为：

[k]＝[kg]+[kep](4)

其中，[kg]为单元几何刚度矩阵，[kep]为单元弹塑性刚度矩阵。

根据虚功原理与刚体准则，推导得到图3所示单元做刚体运动对应的单元几何刚度[kg]为

其中，a为单元截面面积，l为单元长度，iz为单元截面惯性矩，fxb为单元的轴力，mza和mzb分别为单元a端和b端的弯矩。

利用截面逐步屈服塑性铰框架梁单元结点弯矩的增量平衡方程，建立单元弹塑性刚度矩阵：

其中，

k36＝2eizsasbl/(12e²iz²+4eizl(sa+sb)+sasbl²)；

k66＝eizsb(12eiz+4sal)/(12e²iz²+4eizl(sa+sb)+sasbl²)；

其中，e为弹性模量，a为单元截面面积，l为单元长度，iz为截面惯性矩，sa为单元a端的弹簧刚度,sb为单元b端的弹簧刚度。

在弹塑性刚度矩阵[kep]表达式中，塑性铰弹簧刚度sa、sb的确定方法为：

首先通过(2)-(3)式计算单元截面的初始屈服函数值φy和完全屈服函数值φp，判断截面的屈服状态，确定塑性铰弹簧的刚度值：

若初始屈服函数值小于1，单元截面处于完全弹性状态，弹簧刚度保持初始值不变；

若初始屈服函数值大于1，完全屈服函数值小于1时，通过(1)式重新计算塑性铰弹簧刚度sa、sb。

若完全屈服函数值大于1，单元截面达到全屈服状态，弹簧刚度sa、sb退化为0。

步骤3、基于更新拉格朗日格式与刚体准则，建立求解方程的增量迭代过程；

步骤1)将步骤2中建立的各单元刚度矩阵[k]一一对应组集得到结构总体刚度矩阵[k]，基于更新的拉格朗日列式建立结构的增量平衡方程，从中解得结构位移增量即求解

其中，表示第i增量步内第j-1迭代步所形成的结构总体刚度矩阵；为第j迭代步待计算的结构位移增量；为第j迭代步的荷载增量；为节点不平衡力，为第j-1迭代步的荷载向量与结构内力向量的差值。

步骤2)计算各单元在当前单元坐标系下的单元位移增量运用刚体准则的思想，由单元结点位移增量计算出单元结点力即计算下式

其中，为初始状态下的单元结点力列阵，将其乘以几何刚度矩阵，考虑几何大变形的作用，仅改变方向而不改变其数值大小，再叠加自然变形产生的单元结点力{f}，得到变形后的单元结点力

步骤3)将各单元结点力组装得到整体坐标下的结构内力并与荷载比较得到结点不平衡力若结点不平衡力满足收敛要求，则进行下一增量步计算，否则重复本增量步内的预测和结点力计算直至满足收敛要求。

步骤4、迭代求解增量方程，得到结构加载过程的荷载位移曲线以及极限荷载

重复步骤2、3的增量迭代过程，直至结构位移达到预设值，提取结构加载过程中每一增量步得到的荷载与位移值，得到结构的荷载位移曲线以及极限荷载。

为了验证本发明的精确度，现采用四种现有的非线性分析方法和本发明的分析方法对如图5所示的单层单跨平面钢框架结构进行分析。

其中，图5所示的单层单跨平面钢框架结构的高度为5米，跨度为4米的单层单跨平面钢框架结构，初始倾角为ψ＝1/400。已知，钢材的弹性模量e＝205000mpa,屈服强度σy＝235mpa，框架柱与梁的型钢规格分别为heb300与hea340。两柱顶同时承受大小为2800kn的竖向荷载，左侧柱顶承受大小为35kn的水平荷载。实现步骤为：

步骤1、划分每根杆件为一个单元，运用计算机语言matlab编写程序建立相应结构模型；

步骤2、根据本发明，计算单元弹塑性刚度矩阵[kep]和几何刚度矩阵[kg]；

步骤3、基于更新拉格朗日列式，建立结构的增量刚度方程；

步骤4、运用刚体准则迭代求解增量方程，得到结构加载过程的荷载位移曲线以及极限荷载。

图7中，绘出了本发明得到的荷载位移曲线以及常用算法的结果。

表1中，列出了本发明计算得到的极限荷载因子及常用算法的结果。

表1单层单跨框架结构极限荷载因子

[1]vogel塑性区域方法。参见vogelu.calibratingframes[j].stahlbau.1985,10(10):1-7.(vogelu，校准框架，钢结构，1985,10(10):1-7.)

[2]vogel塑性铰方法。参见vogelu.calibratingframes[j].stahlbau.1985,10(10):1-7.(vogelu，校准框架，钢结构，1985,10(10):1-7.)

[3]bradford修正塑性铰方法。参见iuck,bradfordma.higher-ordernon-linearanalysisofsteelstructures.partii:refinedplastichingeformulation[j].advancedsteelconstruction.2012,8(2):183-198.(iuck,bradfordma，钢结构高阶非线性分析.第二部分：修正塑性铰方法，高等钢结构，2012,8(2):183-198.)

[4]商业软件ansys二阶弹塑性方法

采用现有的三种非线性分析方法和本发明的分析方法对如图8所示的四层两度平面钢框架结构进行非线性分析,结构的几何与材料参数如图8所示。

图9中，绘出了本发明得到的荷载位移曲线以及常用算法的结果。

表2中，列出了本发明计算得到的极限荷载因子及常用算法的结果。

表2两跨四层框架结构极限荷载因子

[1]kukreti和zhou改进塑性铰法。参见kukretiar,zhouf.eight-boltendplateconnectionanditsinfluenceonframebehavior[j].engineeringstructures.2006,28(11):1483-1493.(kukretiar,zhouf.八螺栓端板连接及其对框架性能的影响，工程结构，2006,28(11):1483-1493.)

[2]tndoan-ngoc梁柱理论。参见doan-ngoct,dangx,chuq,etal.second-orderplastic-hingeanalysisofplanarsteelframesusingcorotationalbeam-columnelement[j].journalofconstructionalsteelresearch.2016,121:413-426.(doan-ngoct,dangx,chuq,etal.采用协同转动梁柱单元的平面钢框架二阶塑性铰分析，建筑钢材研究杂志，2016,121:413-426.)。

[3]商业软件ansys二阶弹塑性方法(一根杆件划分3个单元)

[4]商业软件ansys二阶弹塑性方法(一根杆件划分10个单元)

通过对比现有方法及商业软件的计算结果可见，对于平面框架结构几何非线性与材料非线性耦合的问题，本发明的计算精度高，所用单元形式简单，物理意义明确，毋需划分大量单元，计算效率高。

最后需要说明的是，上述描述仅仅为本发明的优选实施例，本领域的普通技术人员在本发明的启示下，在不违背本发明宗旨及权利要求的前提下，可以做出多种类似的表示，这样的变换均落入本发明的保护范围之内。