网壳结构由按照一定规律布置的空间杆件构成，常因其自重较轻，受力性能良好，造型丰富，易形成跨度较大的空间等优点而被国内外工程广泛应用[1]。稳定性是网壳结构研究中颇为关键的问题，影响网壳稳定性的主要因素包括：杆件的初始缺陷，高跨比，构件布置形式，材料强度和截面尺寸等[2]。构件的初始缺陷包括结点的装配误差，杆件的初始弯曲，焊接残余应力引起的结构变形等。
非线性有限元法是目前对于网壳稳定性问题应用最广泛的方法，即在考虑几何非线性和材料非线性的情况下跟踪结构的荷载位移曲线，对结构进行全过程的分析[2]。非线性求解算法在形式上可分为完全拉格朗日方法和更新拉格朗日方法[3]。前者在计算结构内力时以结构的初始位形作为参考，不会由于结构变形而造成误差累积；而后者需要利用上一加载步收敛时的位移结果进行下一加载步内力的计算，存在误差的累积[4]。
几何精确梁是一种在大位移和大转动情况下仍能保持位形描述准确的理论模型，更适用于完全拉格朗日的计算格式。其最早由 Ressiner[5]应用于平面二维空间，后又经 Simo[6]推广至三维空间。Ibrahimbegovic[7]将几何精确梁理论推广至空间曲梁结构中。弱形式求积元法[8](简称求积元)是由钟宏志提出的一种基于变分原理，数值积分和微分求积法的数值算法，与传统有限元方法相比不需具体指定位移形函数，在解决复杂几何形状(如变截面或曲梁)、复杂荷载条件、非均匀材料等问题时具有明显优势，在结构分析中有较为广阔的应用前景。
网壳的初始缺陷主要包括结点的位置偏差与杆件的初始弯曲，沈世钊和陈昕[2]通过随机正态分布的方法分析了结点的位置误差对于网壳稳定性的影响。Kani和Mcconnel[9]最早提出利用结构的屈曲模态，控制其幅值大小来模拟网壳的初始缺陷，这种方法也被称为一致缺陷模态法，是用来近似估算存在杆件初弯曲时网壳极限承载力的一种方法。该方法方便快捷，但耦合了结点位置偏差与杆件初始弯曲的缺陷，且杆件实际弯曲状态与模态并不相同，故无法直观显示杆件初始弯曲对于结构承载力的影响程度。曹正罡等[10－12]利用有限元软件ANSYS和一致缺陷模态法研究了不同类型含初始缺陷网壳的稳定承载力。严家川等[13]通过用多段直线近似模拟带有初弯曲的杆件，以完成对网壳稳定性的分析，这种方法相对较为准确，但在有限元软件中建模复杂，且大大增加了杆件单元的数量和计算时间。
本文基于几何精确梁理论构造了一种求解空间曲梁结构的弱形式求积元模型，通过引入纤维模型模拟材料非线性，结合柱面弧长法，采用完全拉格朗日的计算格式，实现对空间曲梁结构的弹塑性大位移大转动分析，避免了使用有限元软件分析时对空间曲梁结构的复杂建模过程，并计算了几类网壳在杆件初弯曲方向角随机，弯曲幅值给定情况下的极限承载力，同时通过参数分析研究了不同弯曲幅值与材料屈服强度的影响，并与JGJ 7—2010《空间网格结构技术规程》[14]建议的一致缺陷模态法及其附录E公式的计算结果进行了对比，弥补了现有计算方法无法便捷单独全面体现初弯曲影响的不足。计算结果表明，极限承载力降低的比例主要与弯曲幅值相关，而在现有钢结构生产工艺下杆件初弯曲的幅值较小，对于稳定性的影响不起控制作用，而《空间网格结构技术规程》[14]建议公式应根据网壳具体形式和材料屈服强度做适当修正。
几何精确梁模型示意图如图1所示，先考虑无初始弯曲的杆件，Ei(i=1,2,3)为初始坐标系的三个方向的单位向量，E3沿初始梁长度方向。图中反映了一根梁变形后的形状，模型的描述主要由截面中轴线的位置向量r和该点对应截面的局部坐标系构成，ei(i=1,2,3)为几何精确梁中定位中轴线坐标r处对应截面的三个方向的定位向量，描述了r处对应截面转动情况。r对应着初始构形中E3坐标某一个确定的值，即给定一个E3，找到初始构形中相应位置的截面，通过r找到该截面变形后的截面中心点的位置，同时根据ei(i=1,2,3)确定截面方向。显然，几何精确梁模型是符合平截面假定的。
对于无初始弯曲的杆件，将其作为一个基本求积元，列出其虚功的内力表达形式如下：
式中：ε为平动应变向量，包含两个方向的剪切应变与轴向应变；n为截面与之对应的内力，包含两个方向的剪力与轴力；κ为转动应变向量，包含两个方向的弯曲曲率与扭转曲率；m为截面与之对应的弯矩，包含两个方向的弯矩与扭矩；L为杆件长度；dL为长度微元。
图1 几何精确梁模型
Fig.1 Geometrically exact beam model
将式(1)中积分区域映射到标准域上，利用弱形式求积元法中常采用的Lobatto积分[15]，可得：
式中：ck代表第k个Lobatto积分点的积分权系数，其他含脚标物理量意义亦类似；s为初始弧长的参数，当杆件无初始弯曲时取为L/2。
定义δde为单元上的基本未知量，包含了n个Lobatto积分点对应位置的δr与δθ(θ为r对应截面处三个方向转动的角度)，其表达形式如下：
根据几何精确梁中应变的表示方法[16]，式(1)中平动应变与转动应变及其变分可写成如下形式：
式中，Ω为截面转动张量，其表达形式如下：
式(4a)中，
表示中轴线坐标r对初始弧长求导，下文中涉及求导的公式均表示对初始弧长求导。当杆件无初始弯曲时，初始弧长为杆件的直杆长。式(4b)中
均为四元数的符号，四元数是由四个参量构成的描述空间转动的物理量[17]，
为转动张量Ω对应的转动四元数；
为该转动四元数的共轭；·为转动四元数的乘法定义；
表示纯四元数；有关四元数的运算法则和定义请参考文献[4]。式(5b)中^代表向量的反对称张量矩阵，式(6)中⊗为张量并矢的符号。
利用微分求积法[15]对式(5a)中含导数的基本未知量在积分点处进行近似，可得下列表达式：
式中：Fki为微分求积公式中的权系数；δki为Dirichlet函数。将式(7a)代入式(5a)后再代入式(2)，整理可得到单元内力的表达式：
假设外力方向不随结构变形而变化，则对式(8)进行变分即可得到单元刚度矩阵表达式：
引入纤维模型[4]考虑材料非线性问题，纤维模型如图2所示。
图2 纤维模型
Fig.2 Fiber model
由于模型的局限性，仅考虑弯曲和轴向拉伸的塑性，而对于剪切和扭转项仍视为弹性。内力表达形式以及式(9a)中的内力变分项可表示如下：
式中：AGx和AGy为截面x方向和y方向的剪切面积；x轴和y轴为截面对应两根主轴；kx和ky为截面关于x轴和y轴的剪切修正系数；Js为截面扭转模量。
将式(10c)代入式(10b)后再代入式(9a)，可得到切向刚度矩阵的表达式如下：
根据式(8)和式(10(a))可以得到在考虑几何非线性与材料非线性情况下单元内力的表达形式，式(11)则给出了单元切向刚度矩阵的表达式。利用单元定位向量将单元刚度矩阵组装成整体刚度矩阵，将单元内力与外力组合成结构整体的内力外力，通过柱面弧长法迭代求解即可。具体步骤与一般非线性有限元求解的步骤相同。与传统有限元方法相比，在求解刚度矩阵表达式时不需具体指定位移形函数，而仅需指定Lobatto积分点的个数。
对于有初始曲率的杆件，只需在前述讨论的基础上对式(4a)、式(4b)稍做修正，将其应变与曲率减去初始的应变与曲率[16]。具体公式如下：
初始状态下的物理量用角标0表达，初始弧长的参数s也按下列公式变化：
式中，S()ξ表示初始弧长S为标准域坐标ξ的函数。
为了验证以上空间直梁与空间曲梁弱形式求积元模型在考虑几何非线性与材料非线性情况下的大位移大转动求解方法的正确性，考虑如下算例。
算例1:Lee框架的弹塑性分析
Lee框架为非线性问题中的一个经典算例[18]，其相关的材料参数如图3所示。其中材料的屈服应力为10.44，屈服前弹性模量为720，屈服后模量为72。本文对于竖向柱采用1个单元计算，横向梁按竖向力左、右两侧分为2个单元计算，每单元采用8个积分点。采用纤维模型划分时，水平和纵向均划分为 15层，为便于直观表达位移相对大小，对水平与竖向位移进行无量纲化：
图3 Lee框架
Fig.3 Lee’s frame
本文采用柱面弧长法进行迭代求解，得到的荷载-位移曲线如图4所示，图中列出了Park和Lee[18]基于几何精确梁和屈服面模型研究该问题得到的结果。由图中结果可以看出，本文算法与屈服面模型研究的算法结果较为接近，误差主要来源于剪切塑性的影响。而无量纲化竖向位移达到0.8左右与屈服面模型结果仍保持很好的一致性，证明了弱形式求积元模型和算法在结构发生大位移与大转角时的有效性。
图4 Lee框架荷载-位移曲线
Fig.4 Load-deflection curve of Lee’s frame
算例2:空间弧形悬臂梁受弯扭
如图 5所示一根空间弧形悬臂梁OA位于E2OE3平面上，构形为1/4段圆弧，圆心位于E2轴上的O1点，于A点作用一扭矩M，方向为-E2方向，材料的相关参数均列于图中。本文对于这一段圆弧采用1个单元计算，该单元采用15个积分点。对于圆管截面，采用纤维模型划分时，环向划分为20份，径向划分为5份。本文采用柱面弧长法进行迭代求解，得到的A点荷载-位移曲线见图6。
图5 90°弧形悬臂梁
Fig.5 90°curved cantilever beam
图6 弧形悬臂梁荷载-位移曲线
Fig.6 Load-deflection curve of a curved cantilever beam
同时使用有限元软件 ABAQUS进行算法结果的检验，计算结果亦见图6。ABAQUS在计算时采用以多段直线等效圆弧的方法，分别划分为3段直线与5段直线进行计算，每段直线上划分为5个单元。结果显示，当扭矩超过80000后材料屈服，结构出现明显的非线性特征。将圆弧划分为3段时存在一定误差，本文结果与将圆弧划分为5段时的计算结果吻合较好，该算例证明了本文算法自由度数量较少，准确性较高，在解决具有初始弯曲构形的结构时不需额外增加单元个数，简化了复杂的曲线建模过程，且杆件数量越多，越能体现该方法的优越性。故较适用于由大量曲梁和直杆构成的杆系结构，可以较为便捷地处理存在杆件初弯曲的网壳稳定性问题。
对于网壳结构，钢结构构件在生产过程中常不可避免地出现初始曲率的情况，而每根杆件在空间中弯曲的方向也不尽相同。随机模拟杆件弯曲方向的方法如图7所示：OBA为一段带有初弯曲的杆件，不失一般性可认为OBA位于Z轴与OA确定的平面内，OBA绕OA旋转θ角得到OCA，θ在(0,2π)中随机取值，OCA即为杆件经过随机旋转后所处的位置。
杆件初弯曲的形状如图8所示。实际情况下钢结构杆件的弯曲形态为一段正弦半波，而由于挠度f0与初始杆长l0之比通常很小，故可用一段半径很大的圆弧近似模拟杆件初始弯曲的形状，P为圆弧圆心，而弯曲程度可由挠曲比，即f0/l0控制。对每根杆件给定挠曲比以及随机旋转的角度θ，可以得到一个杆件初弯曲方向随机的网壳。
图7 杆件初弯曲的方向角
Fig.7 Angle of initial curvature of members
图8 杆件初弯曲形状
Fig.8 Shape of initial curvature of members
图9 ～图11分别为Kiewitt-6(简称K6)型网壳[2]、Geodesic dome网壳[19]以及 Schwedler dome网壳[11]。三种网壳跨度、高度、球面半径、边界条件、材料与截面性质以及荷载均完全相同。下面以 K6型网壳进行说明，其由 156根杆件构成，跨度为30 m，高度为5 m，所有结点均落在半径为25 m的球面上，边界条件为最外侧结点铰接，所有杆件作用竖直向下的均布荷载。暂且不考虑杆件的初始弯曲，分别在材料为弹性及理想弹塑性(屈服强度为345 MPa)的情况下绘制出结构部分结点的荷载位移曲线，并与 ABAQUS的计算结果进行了对比，如图12与图13所示。采用本文算法分析时，每单元采用8个积分点，弹塑性情况下对圆管截面采用纤维模型划分，环向划分为 20份，径向划分为 4份。采用ABAQUS分析时，每根杆件划分为8个单元。可以看出本文计算结果与有限元软件的结果较为吻合，进一步验证了模型和算法的正确性。由图12可知，3号结点竖向位移不断增加，而1号结点竖向位移先增大后减小，这是边侧杆件塌陷所导致，并非弹性卸载的结果，而理想弹塑性情况下中部结点竖向位移则始终增加。下面考虑杆件中存在初始缺陷的情况，先以材料为弹性的K6型网壳进行说明，假定每根杆件挠曲比均为 1/100，初弯曲的方向在空间中随机指定，对于3号结点，每种几何构形都对应一条荷载-位移曲线。取100种初始几何构形，在图中将形成一段连续的区域，取其两侧的包络线绘制在图14中，同时与ABAQUS使用一致缺陷模态法得到的计算结果进行了对比。使用一致缺陷模态法计算时，分别取其一阶与二阶屈曲模态，按挠曲比 1/100的幅值(本例中约为 0.05 m)作为初始缺陷进行求解。该算例中一阶屈曲模态左右对称，故 3号结点荷载-位移曲线仅有一种可能形式，而二阶屈曲模态为左右反对称，故存在位于塌陷侧及位于非塌陷侧两种形式。
图9 K6型网壳
Fig.9 K6 latticed shell
图10 Geodesic dome网壳
Fig.10 Geodesic dome latticed shell
图11 Schwedler dome网壳
Fig.11 Schwedler dome latticed shell
图12 K6网壳荷载-位移曲线(弹性)
Fig.12 Load-deflection curve of K6 latticed shell(elastic)
图13 K6网壳弹塑性荷载-位移曲线
Fig.13 Elastoplastic load-deflection curve of K6 latticed shell
图14 3号结点荷载-位移曲线
Fig.14 Load-deflection curve of node 3
从图14所示可看出，随机方向产生的荷载-位移曲线的包络线能基本包裹住一致缺陷模态法的计算曲线，证明了该方法能反映对称失稳及非对称失稳的不同失稳形式，是有效的计算方法。
网壳的初始缺陷除杆件初始弯曲之外还包括结点的位置偏差，已有文献[2]对结点的位置误差已经通过随机变量模拟的方法做过详尽的论述。而我国《空间网格技术规程》[14]第 4.3.3条明确规定进行网壳全过程分析时应考虑初始几何缺陷的影响，即将结构的一阶屈曲模态作为初始缺陷，缺陷最大计算值可按网壳跨度的1/300取值。这是一种近似估算的方法，杆件实际初始弯曲状态与模态并不相同，且该方法耦合了结点位置偏差与杆件初始弯曲的缺陷，故无法直观显示初始弯曲对于结构承载力的削弱程度。为了进一步研究杆件初弯曲的影响，应在排除结点位置误差干扰的情况下进行单独分析。
下面探究在不同网壳形式下材料屈服强度以及初始缺陷弯曲幅值对于网壳承载力的影响。考虑前文所述的三种网壳，本文对于这三种网壳在不同的材料屈服强度(fy为 235 MPa，345 MPa，390 MPa，420 MPa)以及不同弯曲幅值(挠曲比为 1/100，1/200，1/300，1/500)的情况下进行了弹塑性分析。对于每组参数所对应的网壳，每根杆件挠曲比与材料屈服强度均相同且为给定值，而初弯曲的方向在空间中随机。每组网壳均随机取100种初始几何构形，积分点个数与纤维模型划分方式均与无初始缺陷时相同，对其进行荷载-位移曲线全过程分析，三种网壳的极限承载力在不同材料屈服强度以及不同挠曲比下的分布如图15所示。
图15中亦给出了ABAQUS分别在无初始缺陷和使用规范建议的一致缺陷模态法得到的极限承载力，初始缺陷幅值取为网壳跨度的 1/300，本例中为 0.1 m。上述所求极限承载力均按空间网格技术规程4.3.4要求除以安全系数2。此外，根据空间网格技术规程4.3.5和附录E在图中绘出了规范建议公式的承载力计算结果，具体计算公式如下：
图15 极限承载力分布
Fig.15 Distribution of ultimate load
式中：[qks]为容许承载力标准值，为单位面积上的荷载，再按Be为网壳的等效薄膜刚度；De为网壳的等效抗弯刚度；r为球面的曲率半径，本例中均为25 m。Be与De计算公式较为复杂，与网壳构形相关，但与材料屈服强度无关，详见《空间网格技术规程》附录C。再利用总荷载等效的原则将[qks]转化成单位长度上杆件的线荷载[qks]，计算公式如下：
式中：Sq为网壳在水平面上的投影面积，本文中各网壳均为 706.85 m2；Lq为各杆件在水平面上的投影长度之和。
各网壳相关参数计算值见表1。
由图 15所示结果可看出，不同的网壳形式极限承载力有明显差别：相同条件下Schwedler dome网壳承载力较低，K6型网壳承载力较高，Geodesic dome网壳则位于二者之间。而对于某种网壳，挠曲比越大，其极限承载力的变化范围增大，材料屈服强度越高，其变化范围亦越大。但根据与无缺陷的网壳计算结果对比可以看出，其变化范围主要是位于承载力下降的区间中，即杆件初弯曲对于网壳极限承载力有不利影响，挠曲比越大对于承载力的削减作用越显著，且伴随着屈服强度的增加，由于耦合了更大的结构几何非线性变形的影响故进一步对结构产生不利影响。
表1 网壳相关参数
Table 1 Relevant parameters of latticed shells
网壳类型Be/(kN/m)De/(kN/m)Lq/m [qks]/(kN/m)K6 1.98×1057.59×1026.47×1025.35 Geodesic dome 1.98×1057.57×1026.43×1025.38 Schwedler dome 2.09×1058.01×1026.66×1025.50
本文算法得到的随机数据中最低的极限承载力与一致缺陷模态法的计算结果相对于无初始缺陷时承载力降低的比例，见表2。
表2 承载力下降比例
Table 2 Reductive percentage of ultimate load
网壳类型屈服强度/MPa不同挠曲比下承载力下降比例/(%)1/1001/200 1/300 1/500 一致缺陷模态235 10.57 4.45 2.44 0.96 33.67 K6345 11.50 4.95 2.78 1.17 33.12 390 11.92 5.32 2.97 1.20 33.12 420 12.12 5.47 3.23 1.54 33.10 Geodesic Dome 235 9.28 3.50 1.68 0.43 22.06 345 10.49 4.04 2.02 0.69 21.18 390 10.77 4.07 2.00 0.70 20.83 420 10.90 4.08 1.95 0.58 20.62 Schwedler dome 235 11.56 5.53 3.38 1.58 31.47 345 12.55 6.37 4.17 2.12 32.04 390 12.82 6.61 4.42 2.60 32.12 420 12.93 6.70 4.51 2.72 32.10
结合图15与表2亦可看出：《空间网格结构技术规程》[18]建议幅值的一致缺陷模态法对极限荷载做了明显程度的折减，本文三种网壳承载力下降比例均在 20%～35%区间内，这是因为结点位置的安装偏差即使为毫米量级也可能会对承载力造成较大影响[2]，且还要综合考虑钢结构杆件中残余应力，荷载偏心等其他因素，故初始缺陷幅值的取值较大，较为保守。对于K6型网壳与Geodesic dome网壳，规范建议公式计算结果与高屈服强度下一致缺陷模态法的计算结果较为接近，而对于Schwedler dome网壳，规范建议公式计算结果明显偏高，甚至高于大部分只考虑杆件初弯曲下的承载力计算结果，安全性不足，故建议规范根据该种形式的网壳对公式进行一定程度的修正。此外，由于公式中不包含材料屈服强度这一物理量，而由图15知屈服强度越高时网壳承载力也越高，故也应在公式中适当体现该参量的影响。
由表2亦可知：三种网壳承载力降低的比例差别不大，其中 Schwedler dome网壳对缺陷最为敏感，K6型网壳敏感程度较低，Geodesic dome网壳则位于二者之间。承载力降低的比例随挠曲比的增大和材料屈服强度的增大而增大，其中挠曲比为最主要因素，当挠曲比为1/100时对承载力的削弱达到了10%左右。而当挠曲比为1/500时极限承载力下降比例小于 3%。目前国内钢结构生产工艺的挠曲比为1/2000至1/500[20]，故其对承载力的影响较小，相比较结点位置的安装偏差不为主要的控制因素，可将施工工艺的关注焦点放在控制结点位置的安装偏差上。
本文基于几何精确梁理论构造了一种求解空间曲梁结构的弱形式求积元模型，引入纤维模型模拟材料非线性，结合柱面弧长法，采用完全拉格朗日的计算格式，实现对空间曲梁结构的弹塑性大位移大转动分析，通过算例验证了模型和算法的正确性。与传统有限元多段直线计算曲梁的方法相比，本文算法自由度数量较少、准确性较高，不需额外增加单元个数，与传统有限元方法相比不需具体指定位移形函数，且避免了在有限元软件中复杂的曲线建模过程，结构越复杂，越能体现该方法的优越性。故较适用于由大量直梁和曲梁构成的杆系结构，可以便捷地处理存在杆件初弯曲的网壳稳定性问题。
本文对杆件初弯曲方向进行随机模拟，能体现包括对称失稳与非对称失稳在内的不同失稳形式，计算了 K64型网壳、Geodesic dome网壳及Schwedler dome网壳网壳在不同的材料屈服强度，以及不同的初始缺陷幅值的情况下的极限承载力，通过大量的随机数值模拟实验基本确定在各种初始缺陷幅值下的承载力分布的范围，并与《空间网格结构技术规程》[14]建议的一致缺陷模态法及附录E公式的计算结果进行了对比，排除了结点位置偏差的干扰，弥补了现有计算方法难于便捷单独体现初弯曲影响的不足。计算结果表明：
(1) 文中提到的三种网壳形式在相同条件下极限承载力有明显差别，而承载力因杆件初弯曲缺陷降低的比例则差别较小。承载力降低的比例随挠曲比的增大和材料屈服强度的增大而增大。
(2) 规范建议幅值的一致缺陷模态法对极限荷载做了明显程度的折减，主要是综合考虑结点位置的安装偏差、钢结构杆件中残余应力、荷载偏心等其他因素的影响。
(3) 对于K6型网壳与Geodesic Dome网壳，《空间网格结构技术规程》[14]建议公式计算结果与高屈服强度下一致缺陷模态法的计算结果较为接近；而对于Schwedler dome网壳，《空间网格结构技术规程》[14]建议公式计算结果明显偏高，安全性不足，建议规范根据具体结构形式的网壳对公式进行一定程度的修正，也应适当考虑材料屈服强度对网壳稳定性的影响。
(4) 当挠曲比为 1/100时对不同网壳承载力的削弱的比例约为10%左右。而当挠曲比为1/500时极限承载力下降的比例已经小于 3%。目前国内钢结构生产工艺的挠曲比为1/2000至1/500，故杆件初弯曲相较于结点的位置偏差而言并不为极限承载力的控制因素，可将施工工艺的关注焦点放在控制结点位置的安装偏差上。
[1]符立勇. 大跨度单层球面网壳的稳定性分析[D]. 成都: 西南交通大学, 2002.Fu Liyong. The stability analysis of large-span singlelayer spherical lattice shell [D]. Chengdu: Southwest Jiaotong University, 2002. (in Chinese)
[2]沈世钊, 陈昕. 网壳结构稳定性[M]. 北京: 科学出版社, 1999: 53－56, 88－102.Shen Shizhao, Chen Xin. Stability of reticulated shells[M]. Beijing：Science Press, 1999: 53－56, 88－102. (in Chinese)
[3]Bathe K, Bolourchi S. Large displacement analysis of three-dimensional beam structures [J]. International Journal for Numerical Methods in Engineering, 1979,14(7): 961－986.
[4]张闰. 基于几何精确模型的结构非线性弱形式求积元分析[D]. 北京: 清华大学, 2015.Zhang Run. Nonlinear weak form quadrature element analysis of structures [D]. Beijing: Tsinghua University,2015. (in Chinese)
[5]Reissner E. On one-dimensional finite-strain beam theory: The plane problem [J]. Zeitschrift Für Angewandte Mathematik Und Physik Zamp, 1972,23(5): 795－804.
[6]Simo J C. A finite strain beam formulation. The three-dimensional dynamic problem. Part I [J]. Computer Methods in Applied Mechanics ＆ Engineering, 1985,49(1): 55－70.
[7]Ibrahimbegović A. On finite element implementation of geometrically nonlinear Reissner’s beam theory: threedimensional curved beam elements. [J]. Computer Methods in Applied Mechanics ＆ Engineering, 1995,122(1): 11－26.
[8]Zhong H, Yu T. Flexural vibration analysis of an eccentric annular Mindlin plate [J]. Archive of Applied Mechanics, 2007, 77(4): 185－195.
[9]Kani I M, Mcconnel R E. Collapse of shallow lattice domes [J]. Journal of Structural Engineering, 1987,113(8): 1806－1819.
[10]曹正罡, 范峰, 沈世钊. K6型单层球面网壳的弹塑性稳定[J]. 空间结构, 2005, 11(3): 22－26.Cao Zhenggang, Fan Feng, Shen Shizhao. Elasticplasstic stablity of K6 single layer latticed domes [J].Spatial Structures, 2005, 11(3): 22－26. (in Chinese)
[11]曹正罡, 孙瑛, 范峰. 肋环斜杆型球面网壳弹塑性稳定特性[J]. 哈尔滨工业大学学报, 2007, 39(6): 849－852.Cao Zhenggang, Sun Ying, Fan Feng. Elasto-plastic stability behavior of schwedler domes [J]. Journal of Harbin Institute of Technology, 2007, 39(6): 849－852.(in Chinese)
[12]曹正罡, 孙瑛, 范峰, 等. 单层双曲椭圆抛物面网壳弹塑性稳定性能[J]. 建筑结构学报, 2009, 30(2): 70－76.Cao Zhenggang, Sun Ying, Fan Feng, et al. Elasto-plastic stability of single-layer elliptic paraboloidal reticulated shells [J]. Journal of Building Structures, 2009, 30(2):70－76. (in Chinese)
[13]严佳川, 范峰, 曹正罡. 杆件初弯曲对网壳结构弹塑性稳定性能影响研究[J]. 建筑结构学报, 2012, 33(12):63－71.Yan Jiachuang, Fan Feng, Cao Zhenggang. Research on influence of initial curvature of members on elastoplastic stability of reticulated shells [J]. Journal of Building Structures, 2012, 33(12): 63－71. (in Chinese)
[14]JGJ 7－2010, 空间网格结构技术规程[S]. 北京: 中国建筑工业出版社, 2010.JGJ 7－2010, Technical specification for space frame structures [S]. Beijing: China Architecture ＆ Building Press, 2010. (in Chinese)
[15]Zhang R, Zhong H. Weak form quadrature element analysis of planar slender beams based on geometrically exact beam theory [J]. Archive of Applied Mechanics,2013, 83(9): 1309－1325.
[16]肖乃佳. 基于几何精确梁理论的框架的弱形式求积元分析[D]. 北京: 清华大学, 2011.Xiao Naijia. Weak-form quadrature element analysis of frames based on geometrically exact beam theory [D].Beijing: Tsinghua University, 2011. (in Chinese)
[17]Mcrobie F A, Lasenby J. Simo-Vu Quoc rods using Clifford algebra [J]. International Journal for Numerical Methods in Engineering, 2015, 45(4): 377－398.
[18]Park M S, Lee B C. Geometrically non-linear and elastoplastic three-dimensional shear flexible beam element of von-Mises-type hardening material [J].International Journal for Numerical Methods in Engineering, 2015, 39(3): 383－408.
[19]叶康生, 吴可伟. 空间杆系结构的弹塑性大位移分析[J]. 工程力学, 2013, 30(11): 1－8.Ye Kangsheng, Wu Kewei. Elasto-plastic large displacement analysis of spatial skeletal structures [J].Engineering Mechanics, 2013, 30(11): 1－8. (in Chinese)
[20]陈绍蕃, 顾强. 钢结构. 上册, 钢结构基础[M]. 第 3版. 北京: 中国建筑工业出版社, 2014: 88.Chen Shaofan, Gu Qiang. Steel structures volume one:Basic knowledge of steel structures [M]. 3rd ed. Beijing:China Architecture ＆ Building Press, 2014: 88. (in Chinese)
ELASTO-PLASTIC STABILITY ANALYSIS OF LATTICED SHEELS WITH INITIALLY CURVED MEMBERS BY WEAK-FORM QUADRATURE ELEMENTS