基金项目:国家自然科学基金资助项目(51908356); 浙江省自然科学基金资助项目(LQ19E080013); 中国博士后科学基金资助项目(2019M662056); 绍兴文理学院国际科技合作项目(2019LGGH1005).
第一作者:冯晓东(1987-),男,博士,副教授,研究方向为大跨度空间结构及智能可展结构.E-mail:xiaodong.feng@csu.edu.cn
(1.绍兴文理学院 土木工程学院,浙江 绍兴 312000; 2.绍兴市城投建筑工业化制造有限公司,浙江 绍兴 312000; 3.浙江交工集团股份有限公司,浙江 杭州 310000; 4.浙江省建工集团有限责任公司,浙江 杭州 310000)
(1.School of Civil Engineering,Shaoxing University,Shaoxing 312000,China; 2.Shaoxing City Investment & Construction Industrial Manufacturing Co.,Ltd,Shaoxing 312000,China; 3.Zhejiang Communications Construction Group Co.,Ltd,Hangzhou 310000,China; 4.Zhejiang Construction Engineering Group Co.,Ltd,Hangzhou 310000,China)
tensegrity structure; force density form-finding; geometric symmetry; singular value decomposition
DOI: 10.15986/j.1006-7930.2021.04.006
张拉整体结构是一种由连续的拉索与离散的压杆组合而成的空间索杆张力结构,构件内力以轴力的形式存在,可高效利用结构材料的拉压性能.也正因这种简单的受力形式和高效的材料利用率,使该结构的设计思想被应用于索穹顶结构[1].此外,张拉整体结构具有可折叠性,在产生大变形位移时只产生微小应变,因此在机器人[2]、航天航空[3]及生物学领域[4]都得以广泛应用.
初始预应力为张拉整体结构提供刚度,使其处于一种自平衡状态,这一过程通常称之为找形.在现有找形方法中往往以节点坐标、拓扑连接或初始预应力为变量,并根据变量的不同划分为拓扑找形和预应力找形.找形策略上则大多采用优化算法、动力松弛法、数值迭代法等求解手段获取张拉整体结构的初始形态.在拓扑找形研究中,Gan[5]利用遗传算法寻找非对称的张拉整体结构.谢胜达[6]基于小生境遗传算法对环形张拉整体结构进行拓扑优化.Yin[7]采用直接加杆和递归加杆方式来构造简单多边形张拉整体结构.在预应力找形研究中,早期由陈志华[8]提出张拉整体结构的力密度法,可有效解决带有索元的柔性空间结构初始形态确定问题.Estrada[9]基于力密度法以最小弹性势能为目标函数,采用矩阵相互迭代的方式求解结构的初始力密度.Tran[10-12]结合力密度矩阵的特征值分解和平衡矩阵的奇异值分解,提出了一系列求解张拉整体结构初始形态的方法.尚仁杰[13]则通过几何作图的方式求解张拉整体结构的初始预应力大小.Wang[14]提出一种找形计算框架,无需考虑自应力模态就能处理满足秩亏条件的力密度矩阵.K. Koohestani[15]提出一种基于LU分解的数值计算找形方法.
综上所述,张拉整体结构在找形方面的研究已较为成熟,构形简单的张拉整体结构可快速求解其初始形态,但在求解形态复杂的张拉整体结构的问题上,现有的方法在求解能力上略显不足.因此,本文将基于力密度法和矩阵分析理论,根据结构对称性,对构件分组归类,通过二次奇异值分解的方式求解张拉整体结构的初始力密度,为复杂张拉整体结构的初始形态确定提供一条新的途径.同时本文还提出一种新型张拉整体结构构形,采用不同分组计算该结构的初始状态,并制作实际模型以验证该结构的合理性,丰富了现有张拉整体结构形态的分析理论与应用实践.
本文中张拉整体结构找形计算遵循以下基本假定:
(1)各单元连接方式为铰接;
(2)不考虑自重及外部荷载;
(3)不考虑压杆屈曲以及材料非线性;
(4)不考虑结构的局部或整体失稳.
对于任意张拉整体结构已知有n个自由节点、nf个固定节点、b个单元.根据拓扑连接关系可知其关联矩阵为Cb×n,假定连接单元k两端节点为i和j,则关联矩阵可如下表述:
定义Q=diag(q),其中,q为构件力密度,根据单元的受力情况可建立平衡方程:
式中,Q为力密度矩阵; X,Y,Z(Xf,Yf,Zf)为结构中自有节点(固定节点)在x,y,z方向上的节点坐标矩阵; Cf为固定节点的关联矩阵.
由于张拉整体结构为自平衡结构,节点均为自由节点,因此式(2)可转化为
再将上式简化为
式中,A为平衡矩阵.
结构的自应力模态数s和独立位移模态数m可根据平衡矩阵A求得:
s=b-rA (5)
m=dn-rA (6)
对平衡矩阵A进行奇异值分解得
A=UVWT (7)
U=[u1 u2 L udn] (8)
W=[w1 w2 L wb] (9)
其中U,W为正交矩阵,V为由矩阵A的非负奇异值组成的对角矩阵.且此时的自应力模态可分为两种情况[16]:
(1)自应力模态数为1(-b-rA)
在这种情况下,正交矩阵W可表示为
其中n1为力密度系数向量,满足线性齐次方程(4).
(2)自应力模态数大于1(=b-rA)
在多自应力模态状态下,正交矩阵可表示为
此时,可行自应力可由向量基的线性组合表示为:
q=α1n1+α2n2+…+αsns (12)
在平衡矩阵A的计算过程中并不考虑结构材料特性,经奇异值分解得到的多组向量基可通过不同的线性组合,将存在无数种可行自应力结果,但并非必定满足结构拉压构件关系.为解决此问题,则需事先考虑结构的拉压关系,同时结合几何对称性,将构件划分为数个组别.
定义h为分组数量,则力密度矩阵可表示为
其中:ζi中第i个元素为1,矩阵内其他元素为0; ζ~为力密度拉压关系系数矩阵.
将式(12)代入式(14)得
其中q~中所有元素均为非负数.
因此,式(16)可以写成
A~q~=0 (19)
此时,线性齐次方程的所有解则取决于A的零空间矩阵维数[17].
s~=s+h-rA~ (20)
rA~=rank(A~) (21)
此处所得到的s~与自应力模态数s均可表示为矩阵的秩亏数.
为求解式(19),对A~进行奇异值分解得:
其中U~,W~为正交矩阵,V~为由矩阵A~的非负奇异值组成的对角矩阵.此时的s~可分为三种情况:
(1)s~=0
在该情况下,式(19)无解,则说明几何对称性组别划分不合理,即该张拉整体结构在这种几何对称性下并不处于自平衡状态.为了得到可行解,需要增加组别划分数量以此增大秩亏值直到以下两种情况的出现.
(2)s~=1
对于s~=1,存在唯一解.该组结果不仅满足结构自平衡,而且满足结构的几何对称约束.此时在自平衡状态下杆件均满足压杆受压、拉索受拉.
(3)s~>1
对于s~>1,向量基的线性组合具有无数种可能,即多自应力模态状态.值得注意的是,将一个张拉整体结构划分到s~=1的情况需要花费大量时间.而对于s~>1,只需根据几何对称性和工程要求进行划分,这样的划分不仅能提高设计工作的效率,还可以有效地指导实际施工.
定义χ=[χ1,χ2,…,χh],表示h组可行自应力模态向量基系数.分组后的自应力计算如下式表示:
其中,q~i(i=1,…,h)为从式(23)中rA~+1到s+h部分中选取s+1至s+h部分.
采用内点法求解线性规划问题,并以q~h内所有元素大于0为约束条件,力密度标准差最小为目标函数,实现可行预应力的求解.需要注意的是,对于得到的各组力密度值仍需根据式(25)进行验证,其中本文设计误差需满足k≤10-13
为保证张拉整体结构在可行自应力状态下的几何稳定,可根据弹性刚度矩阵和几何刚度矩阵得到切线刚度矩阵[18-19],其中本文拉索和压杆的轴向刚度分别采用eciaci=125 kN,esiasi=1 200 kN.
L=diag(I) (27)
E=CTdiag(q)C (28)
式中,I为单元长度的矢量表达; e为单元杨氏模量的矢量表达; a为单元横截面积的矢量表达; I3为三维单位矩阵; E为力密度矩阵; 为张量积.
若结构稳定,则KT正定.所以对于任意一个无穷小位移D,矩阵KT的二次型是正定的.
由于KE和KE处于非满秩状态,且零空间包含刚体位移,因此需忽略前6个刚体位移情况,通过最小特征值λ1反映结构的刚度强弱情况.
本文提出一种基于立方多面体结构的新型张拉整体结构.该结构共有12个节点,每个节点均连接6根构件,共有12根压杆(红色)和24根拉索(蓝色),结构拓扑连接方式如图1,具体节点坐标如表1.通过式(4)、(5)可得到该结构的自应力模态数为6,属于多自应力模态结构,因此其初始力密度存在多种可能.采用本文所提出的找形算法求解该张拉整体结构的初始力密度,需事先根据结构的几何对称性对构件进行分组.分组情况的不同可能会导致初始力密度的找形结果不同.
根据几何对称性将索单元分为4组,杆单元分为1组,具体分组如图2所示.经1.3节公式计算得s~=1,因此可说明该结构在此对称分组状态下存在唯一初始力密度,并与Zhou[20-21]所提出的几何对称张拉整体结构找形方法求解结果一致,如表2,其中设计误差k=1.868 0e-1满足设计要求.但文献中所阐述的找形方法无法事先确定结构构件的拉压关系,且找形结果单一.而本文可根据不同的对称分组情况求解出多种初始力密度结果,这也正是本文算法的优势所在.
在上述初始力密度作用下,可根据式(26)考虑该结构的稳定性,其中切向刚度矩阵特征值如表3所示,从表中可以看出最小特征值为正值,则说明结构处于稳定状态.
为进一步阐述本文算法优势,展现多种初始力密度状态的求解能力,可将空间12杆张拉整体结构划分为12组,具体分组情况如图4.此时计算得到的s~=2,则说明该分组情况下具有多种初始力密度,因此需要以分组构件力密度标准差作为目标函数进行优化(图3),优化结果如表4,并根据式(25)求得误差值k=7.474 8e14,满足设计要求.
表3 12杆张拉整体结构切线刚度矩阵特征值(h=5)
Tab.3 Eigenvalues of the tangent stiffness matrix of the 12-strut tensegrity structure(h=5)
对比两种分组结果,由于分组的不同,最终得到的力密度值也会随之改变,从表5中可以看出h=5情况下结构的切线刚度矩阵最小特征值仍为正值,因此可以保证在该力密度作用下结构仍处于稳定状态,同时对比表3,在均满足结构稳定的条件下,12组分组(h=12)情况下结构所形成的刚度要略优于5组分组(h=5)的结构.
现将结构分组增至36组,即单根构件各为一组(以不同颜色表示),如图5.通过计算得到s~=6,且与未分组时的自应力模态数s相同,则说明构件分组的方法并未起到减小自应力模态数的作用.但此时分组的目的仅是为了确定结构内部的拉压关系,进而确保计算结果满足结构杆件的拉压关系.
通过内点法优化(图6),可以看出力密度标准差处于波动下降,在经过60次迭代后趋于稳定,同时一阶最优值在前15次迭代时已趋于0.最终得到的初始力密度值计算设计误差为k=8.892 0e-14,具体力密度数值如表6所示.
表5 空间12杆张拉整体结构切线刚度矩阵特征值(h=12)
Tab.5 Eigenvalues of the tangent stiffness matrix of the 12-struts tensegrity structure in space(h=12)
同时需对结构进行稳定性计算,求解得到切线刚度最小特征值为20.4484大于0,则说明结构稳定,具体数据如表7.单根构件各为一组的情况同样符合非对称结构的分组求解方式.因此,本文算法也可适用于非对称结构的初始力密度求解.
由于张拉整体结构仍处于理论计算阶段,因此仍需制作模型对本文所提出的新型张拉整体结构加以验证.在模型制作过程中采用3 mm直径圆柱木棒为压杆,1 mm直径弹力绳为拉索.在模型制作完成后,尝试改变弹力绳拉伸长度,用于模拟拉索上不同初始力密度的施加状态,从而进一步说明在不同初始力密度状态下,该张拉整体结构的稳定状态.
本文所提出的新型张拉整体结构构形,由12根压杆与24根拉索组成.模型制作之前,应事先对结构尺寸进行确定.在理论模型的计算结果中,对结构构件长度进行化一处理,便可获得实际制作中的构件长度比例为l1~l24:l25~l28:l29~l36=1:21/2:31/2,其中,li为第i个单元的长度.根据h=5的计算结果,索单元力密度比例为qh1:qh2:qh3:qh4=1:0.5:0.5:1,其中,qhi为第i组构件的力密度值.同时根据式(31)~(33)求解各个索单元的初始长度:
其中,qi,j为第i,j个单元的力密度; li,j为第i,j个单元最终成形长度; li,j0为第i,j个单元的初始长度; e为弹性模量; a为横截面积.
联立两个方程可得:
由于理论模型中最终索单元长度均一致,所以li=lj,且设置初始第1组索长为10.0 cm,最终索长为15.0 cm,根据式(33),可求解出第2组初始索长为12.5 cm.而后对结构计算模型进行力学模型简化,完成初始模型制作,具体流程如下:
(1)裁剪16段长度为12.0±2 cm的弹力绳和8段10.0±2 cm的弹力绳,其中2 cm为预留长度,用于拉伸后的长度调节或节点绑扎固定,同时需保证拉伸后弹力绳长度均可达到15.0 cm.再根据计算得到的最终构件长度比例,还需裁剪4段长度为21.2 cm的木棒和8段长度为26.0 cm的木棒,以及12段3 cm塑料胶管,作为节点连接套筒,具体如图7所示,其中忽略木棒受压时的长度变化;
(2)可先制作4组由1根21.2 cm的木棒与2根26.0 cm的木棒用塑胶套筒连接形成的三角构件;
(3)根据理论模型的拓扑连接形式,利用弹力绳将4个三角构件连接成形.其中,索杆节点采用绑扎的方式进行固定,同时剪掉多余部分;
(4)当所有单元都已按照理论模型连接后,需测量各个拉索长度是否达到预设长度.若未达到,则需解开绑扎节点,重新调整弹力绳长度,而后重新绑扎,直至所有拉伸后的弹力绳长度达到预设长度,最终模型如图8所示.
当按照不同分组结果制作模型时,可根据式(33),事先假定一组构件初始长度,通过改变弹力绳的初始长度,进而改变各单元的力密度值.
(1)本文的找形方法对几何对称张拉整体结构具有高效的求解能力,同时也适用于指定拉压构件的几何非对称张拉整体结构的初始力密度找形.
(2)本文算法将空间十二杆张拉整体结构构件分别划分为5组,12组以及36组,以力密度标准差为目标函数,均可求解出最优初始力密度,证实了该算法可有效求解几何对称张拉整体结构的初始力密度.
(3)基于找形结果,利用弹力线模拟拉索,木棒模拟压杆,制作空间十二杆张拉整体结构,验证了本文算法的有效性及合理性.