王兆毅， 王宪杰，2*， 张帆， 赵颖， 奂瑾， 江鑫禹
(1.云南大学建筑与规划学院， 昆明 650000； 2.浙江大学建筑工程学院， 杭州 310058)
大跨度空间网壳结构得益于其覆盖空间大、造型新颖、自重轻等优点，成为工程师们非常乐意采用的一种结构形式。单层网壳由于厚度薄，刚度相对较小，几何非线性敏感。所以稳定性分析是单层网壳设计中的关键问题。以往屈曲分析时，设计师们不得不依赖于连续化理论(“拟壳法”)将网壳转化为连续壳体结构，通过近似的非线性解析方法来求出壳体的稳定性承载力，但这种方法有较大的局限性，缺乏统一的理论模式，事实上仅对少数特定的壳体才适用[1]。随着计算机的发展及应用，非线性有限元分析方法转而成为实现结构稳定性分析的主要方法,各种改进的弧长法[2]是这方面的一个重要成果，其为结构的荷载-位移全过程路径跟踪提供了最为有效的计算方法，并能准确的确定其稳定性极限承载能力。
以往单层网壳稳定性分析往往停留在宏观层级上，无法反映细观受力机理。然而结构的失效破坏往往源自材料的局部失效，进而导致关键节点或杆件丧失承载能力，最终导致整体结构的破坏。可见，结构的失效过程本质上是宏观层面与细观尺度之间相互耦合影响的结果[3]。如何在跨尺度连接处施以精确的约束耦合关系，在不损失各尺度的自由度的前提下，建立多尺度计算模型并有效模拟宏观-细观的非线性演化过程是非常必要的。陆新征等[4-5]通过Uforms子程序实现了跨尺度界面的连接，通过多尺度模型与梁单元模型时程分析结构对比，验证了多尺度模型可以更加准确地模拟节点的实际受力情况。石永久等[6]探讨了多尺度模型在钢框架抗震计算的应用问题，实现了不同尺度模型间的有效连接，通过与国内外典型节点试验进行分析对比，验证了细观单元与宏观模型界面连接方法的有效性。罗尧治等[7]利用平截面假定推算了跨尺度连接界面的位移协调方程，并给出了与焊接球连接的杆件长的参考值。朱南海等[8-9]用宏观单元对网壳节点及杆件的薄弱区进行离散，建立网壳结构的多尺度模型，进行了多尺度建模下的网壳静力稳定性及地震损伤过程分析。结果均表明了多尺度模型具有更高的精确度。刘春城等[10]提出了一种易于有限元软件实现的跨尺度单元连接法，建立了结构行为一致多尺度模型，验证了该模型的合理性。 Yu等[11]验算了多点约束法界面连接控制方程，研究结果显示多点约束法可用于实现跨尺度界面的有效连接。上述研究均表明：多尺度模型兼顾计算效率与时间成本。但鲜有提及多尺度模型是否影响结构基本特性这一问题，并没有针对结构不同的重点响应区域建立不同的精细单元，难以全面考察多尺度计算模型的优越性。
基于并行一致多尺度法，建立一个能有效反映细观力学性能响应多尺度计算模型，对实际网壳进行模态分析及荷载-位移全过程非线性屈曲分析。通过全梁单元模型与多尺度模型的分析对比，为多尺度计算模型的发展与工程设计应用提供理论支撑。
多尺度模型构建思路大致上能分为两种：信息传递多尺度方法与并行一致多尺度方法。前者主要是提取宏观结构的各项响应作为精细构件的边界条件，旨在实现不同尺度之间力学性能参数的交互。后者建模策略是对结构薄弱区的构件采用细观尺度的精细模型，一般用实体单元或是壳单元离散，宏观结构采用梁单元等杆系单元离散。目的是同时建立细观与宏观尺度模型，实现跨尺度单元之间连接界面的耦合，以达到整体结果的同步计算分析。并行一致多尺度法的关键点是如何在不损失宏观单元自由度的同时不增加精细单元的额外约束。而多点约束法是有效解耦连接界面处位移协调方程的方法之一。
多点约束法的基本思想是在跨尺度连接界面上利用位移协调条件推导出约束方程，进而实现不同单元的耦合。以宏观单元节点的自由度作为控制值，建立针对细观单元各节点相应的自由度及位移值的控制方程，便可用于表征不同单元之间的力学响应及自由度的传递。且在小变形问题中，多点约束法不用迭代，节约了计算时间。在大变形问题中，多点约束法并不会引起连接处刚度突变，保证了计算效率和计算精度，在求解非线性迭代问题中具有很强的鲁棒性。
Abaqus中提供了Tie、Coupling、MPC(multi-point constrain)等约束方法，其中MPC允许在模型不同自由度中强加约束，将其耦合，可以实现不相容单元荷载、位移的有效、精确传递，且更易考察连接处的应力及位移响应。
基于多尺度模型，建立连接面处的局部坐标系如图1所示，在连接面建立基准梁单元，宏观节点处共有6个自由度{u1,u2,u3,u4,u5,u6}；精细单元节点有3个自由度{u1,u2,u3}；多点约束法连接界面节点自由度的约束方程可写为
图1 跨尺度界面连接示意图Fig.1 Schematic diagram of cross-scale interface connection
C(uB,uS)=uB-CuS=0
(1)
式(1)中：uB为宏观单元节点自由度，为主自由度；uS为精细单元节点自由度，为从自由度；C为约束方程系数矩阵；经简化后的多点约束连接式是一个一次多项式，具体方程为
uB=C1uS1+C2uS2+…+CnuSn+C0
(2)
对于解耦线性约束方程，常用位移协调一致法求解，引入局部坐标系xyz，以B为原点，并符合右手螺旋法则。基于平截面原则，任何精细单元节点位移满足关系式
(3)
式(2)中：xS、yS、zS分别为精细单元节点在局部坐标系中的x、y、z坐标；yB、zB为宏观单元节点在局部坐标系的y、z坐标；βi为精细化模型节点对宏观单元节点的影响权重系数，可以取βi=1/n。
由式(3)可知，宏观单元节点在微观单元连接界面的沿杆件轴向位移、沿平面内外切向位移，由精细单元界面上的全部节点决定。
屈曲分析的目的是确定结构从稳定的平衡状态变为不稳定的平衡状态时的临界荷载及其屈曲模态。采用一致缺陷模态法定义初始缺陷时，先对结构进行特征值屈曲分析和缺陷结构的非线性全过程分析。
线性屈曲分析可以得到理想线性结构的理论屈曲强度，无须进行非线性分析，就可得到临界荷载及屈曲模态，可为非线性屈曲分析提供参考荷载值。其控制方程为
(KL+λKσ)ψ=0
(4)
式(4)中：λ为特征值，即通常意义上的荷载因子；ψ为特征位移向量；KL为结构的小位移刚度矩阵；Kσ为参考初应力矩阵。
为考虑初始缺陷对理论屈曲强度的影响，需对结构进行基于大挠度理论的几何非线性屈曲分析，其控制方法为
KT×Δu=F-r
(5)
式(5)中：KT为单元切线刚度矩阵；Δu为位移质量矩阵；F为外力矩阵；r为残余力矩阵。
本文中所用实例为唐山勒泰中心项目[12]上方分割频率为6的凯威特K8单层网壳，网壳跨度为76.8 m，矢高为12 m，矢跨比为1/6.4。杆件为Q345B圆钢管，管径为168 mm，壁厚为6 mm，节点形式为Q345B焊接空心球节点，外径为420 mm，壁厚为12 mm。节点数为185，杆件数为465。外圈节点铰接。在ABAQUS中屈曲分析得出的屈曲因子为施加荷载的放大系数，为方便考察稳定性极限荷载，本例中在网壳节点处施加单位集中荷载。有限元模型如图2所示。
对网壳进行模态分析，得出其振型与自振频率、周期等固有特性(图3、表1)，便于验证多尺度模型的适用性与正确性。
通过计算可知式(4)中λ=247.6，由于施加的是单位荷载，故该网壳理想临界荷载为247.6 kN。特征值屈曲分析不考虑双重非线性及初始缺陷，所得出的结果只是屈曲荷载的上限，是一种理论解，得不到结构的实际承载力。这里关注的是前几阶屈曲模态，尤其是第一阶屈曲模态(图4)，结构的一阶屈曲模态是结构屈曲的位移倾向，是潜在的位移趋势，如果结构的初始缺陷分布与结构的第一阶屈曲模态相吻合，无疑对结构的稳定性产生最不利影响，故一阶屈曲模态可用来模拟网壳结构中的缺陷分布。
图2 有限元模型[12]Fig.2 Finite element model[12]
图3 全梁单元前四阶振型Fig.3 The first four modes of the full beam element
采用“一致缺陷模态法”引入初始缺陷，根据《空间网格结构技术规程》[13]的相关规定：初始缺陷最大计算值按照网壳跨度L的1/300，取256 mm。有限元软件中，实现引入初始缺陷的关键是定义Imperfection函数。分析过程考虑几何非线性，采用两段式理想弹塑性本构，屈服强度为345 MPa，满足von Mises准则。为考察结构实际节点临界荷载，同样施以单位荷载。考虑几何、材料双重非线性，对结构进行荷载-位移全过程稳定性分析。
通过观察图5可知，进入屈服阶段的杆件主要是第四圈部分环向杆及小部分第一圈及第二圈内的斜向杆，其中第四圈环向杆与径向杆相交的节点的塑形程度最为严重，应力集中较为显著。底部环向杆与斜向杆仍基本处于弹性阶段，其应力响应不明显。位移形变量最大的区域集中于第三圈环杆与径向杆相交的节点，在节点处出现局部的坍缩。观察图6可知当节点荷载增大到114.16 kN时，结构的平衡状态变为不稳定的平衡状态，随后随着荷载不断降低，位移不断增大，表明荷载必须不断下降才能维持结构内外力之间的平衡，此阶段说明结构已经产生失稳破坏，承载力不断上升，位移不断增大。故该网壳结构失稳时的节点临界荷载为114.16 kN,临界时刻的结构最大位移为204.8 mm。
由以上分析可以看出，该网壳结构分别有位移响应集中的重点区域，主要为于第三圈环杆与径向杆相交的4个节点。应力响应集中的重点区域，主要为第四圈环向杆与径向杆相交的8个节点。
表1 全梁单元自振频率与周期Table 1 The natural frequency and period of the full beam element
图4 全梁单元一阶屈曲模态Fig.4 The first buckling mode of the full beam element
图5 全梁单元失稳时刻力学响应Fig.5 Mechanical response at the moment of instability of the full beam element
图6 全梁单元荷载位移曲线Fig.6 Load-displacement curve of full beam element
进行多尺度双重非线性静力稳定性分析时，应先对薄弱区域进行精细化建模(图7)。取壳单元建模时，应力分布受边界效应影响显著，节点塑性发展程度亦可能受边界效应干扰，本文中讨论杆件长度对计算结果的影响，采用文献[7]中研究结果作为建模参数，取空心球节点及与球节点相接处0.9倍杆件外径的杆件作为需精细化建模区域[7]。精细单元采用S4R壳单元建模，宏观单元采用B31梁单元建模，多尺度模型在结构尺寸。材料本构与全梁单元一致。表2为多尺度模型分类一览，参照一致单元模型，对多尺度有限元模型进行模态分析及考虑初始几何缺陷的双重非线性屈曲分析。
图7 多尺度模型及精细化节点连接大样Fig.7 Multi-scale model and the large-scale connection of refined nodes
表2 多尺度模型分类依据
利用模态分析得出多尺度模型的结构固有特性，如图8及图9所示，多尺度模型与全梁单元模型的低阶振型存在明显差别，全梁单元前四阶振型存在一定的对称性，但该特性并没有明显地在多尺度模型振型中体现。究其原理是因为多尺度模型建立具体详尽的节点参数，且由具体的球节点替代了单一梁单元交汇节点，导致结构整体刚度分布随着不同数目及区域的精细化节点的建立而发生改变。观察表3、表4和图10可知，多尺度模型1、2的一阶自振频率与全梁一致单元模型相比分别降低了0.42%和0.13%，且前十阶自振频率相差无几，由此可见，多尺度模型的建立不会显著改变结构的自振频率及周期等固有特性，更进一步验证了多尺度计算模型的正确性。
图8 多尺度模型1前4阶振型Fig.8 The first 4 modes of multi-scale model 1
图9 多尺度模型2前4阶振型Fig.9 The first 4 modes of multi-scale model 2
表3 多尺度模型1自振频率及周期Table 3 Multi-scale model 1 natural vibration frequency and period
表4 多尺度模型2自振频率及周期Table 4 Multi-scale model 2 natural vibration frequency and period
图10 三种计算模型前10阶自振频率对比Fig.10 Comparison of the first ten natural frequencies of the three calculation models
同一致单元模型的非线性屈曲一致，同样以“一致缺陷模态法” 引入初始缺陷，考虑双重分线性，对多尺度模型1、2分别进行非线性屈曲分析。
观察图11可看出，多尺度模型一阶线性屈曲模态是基本一致的，说明引入的初始几何缺陷也基本相同。但理想临界荷载却有所降低，多尺度模型1、2临界荷载分别为213.59 kN和190.09 kN，较全梁单元临界荷载分别降低了13.7%和23.2%。
由图12及图13可知，多尺度模型1、2与全梁单元模型在失稳时刻应力分布基本一致，进入屈服阶段的杆件主要是第四圈部分环向杆，外圈杆件基本处于弹性阶段，其中第四圈环向杆与径向杆相交的节点塑形最为显著。但不同的是多尺度模型表现出了节点与杆件处出现明显的应力集中，说明杆件与节点连接处为薄弱区域，这一结果与实际工程情况相吻合。且杆件中也出现了局部区域的屈服及应力集中，此部分杆件发生变形，且出现局部凹陷。
图11 多尺度模型最低阶线性屈曲模态Fig.11 The lowest-order linear buckling mode of the multi-scale model
图12 多尺度模型1失稳时刻应力分布及节点大样Fig.12 Multi-scale model 1 stress distribution and node size at the moment of instability
观察图14及图15可知，多尺度模型与一致单元模型在失稳时刻位移云图分布大致相同，尽管位移最大的点发生转移，但位移形变量最大的区域都集中于第二、三圈环杆与径向杆相交的节点附近。
图13 多尺度模型2失稳时刻应力分布及节点大样Fig.13 Multi-scale model 2 stress distribution and node size at the moment of instability
图14 多尺度模型1失稳时刻位移分布及节点大样Fig.14 Multi-scale model 1 displacement distribution and node size at the moment of instability
图15 多尺度模型2失稳时刻位移分布及节点大样Fig.15 Multi-scale model 2 displacement distribution and node size at the moment of instability
图16 三种计算模型下荷载-位移曲线Fig.16 Load-displacement curves under the three calculation models
在精细化区域部分，杆件的位移比球节点的位移更为显著，故杆件与节点连接处易出现局部形变，位移变形并不协调。
上述分析表明：单层网壳的整体失稳，往往伴随着节点与杆件相接处薄弱区的应力集中和杆件的变形。且由于精细节点的引入，导致结构局部区域刚度变化，但一致全梁单元模型并没有考虑节点局部变形及应力集中，采用梁单元交汇节点会忽略空间力系交汇处的复杂受力及变形。采用多尺度模型，不仅可以模拟细观杆件的力学性能变化，也能得出较为准确的宏观结构力学响应分布。
由图16可知，多尺度模型1、2的节点临界荷载与全梁单元相比分别降低了5.48%和15.06%，且结构稳定承载能力随着引入更多的精细化节点而下降。主要是因为引入具体形状及厚度参数的精细化模型，改变了整体结构的刚度分布，又考虑局部杆件的屈服和节点与杆件相交处的应力集中，使得结构杆件传力途径发生变化，特别是随着节点的变形及杆件的局部塌陷，进一步加剧了结构失稳破坏，使得多尺度模型还没有达到一致单元模型的极限承载力就发生失稳破坏。由于更多布置在应力集中处的精细化节点的引入，节点塑性发展程度加深，故多尺度模型2的稳定性承载能力比多尺度模型1更低。
以上分析表明：多尺度模型不仅能模拟宏观结构的力学响应分布，也能准确模拟出细观节点的塑性开展。
建立了针对应力、位移响应集中的重点区域进行精细化建模的多尺度模型，比较多尺度模型与全梁单元模型的振型、固有特性和在考虑初始缺陷的非线性屈曲分析中的各项响应及稳定性，验证了在大跨度空间结构和网壳结构性能研究中多尺度建模方法的可行性、准确性以及多尺度建模方法在局部细观力学响应分析时的优越性。主要结论如下：
(1)相对于全梁单元模型而言，多尺度模型兼顾计算效率与计算精度，采用多尺度模型不会显著改变结构的自振频率及周期等固有特性,但由于细化了节点参数的精细化节点域的引入会导致结构的刚度分布发生变化，进而导致低阶振型的改变。
(2)多尺度模型与全梁单元模型在失稳时刻的应力分布、位移分布情况基本一致，在宏观结构分析上，多尺度模型也有一定的精度。由于引入精细化模型考虑了局部杆件的屈服和节点的应力集中，使得节点发生变形、杆件出现局部塌陷，导致多尺度模型的承载力比全梁单元有所降低。
(3)多尺度模型相较于全梁单元，既能精确地确定重点响应区节点的应力、位移边界条件，也能较为清楚地显示出大跨度空间结构细部节点的应力开展情况及破坏过程，为空间结构的分级优化设计提供了新的途径。