具有约束支撑的张拉整体结构初始找形方法

冯晓东 赵容舟 黄世荣 刘圣威

(1.绍兴文理学院 土木工程学院，浙江 绍兴312000；2.浙江精工钢结构集团有限公司 技术中心，浙江 绍兴312000)

0 引言

张拉整体结构(Tensegrity structures)作为新型大跨度空间结构凭借其构形独特、质量轻盈、可折叠、形态可控等特点正吸引着世界各国学者的广泛关注，成为当今空间结构领域中最前沿的课题.通常而言，由于杆件中自应力的存在，使得其在任何外力(包括重力)作用之前，整个结构体系能够保持自平衡状态，因此张拉整体结构是典型的自平衡结构.此外，这类结构所具有的一系列潜在特点，刚度可协调性、主动阻尼共振、形态可控性，使得“张拉整体”的概念迅速贯穿于包括但不仅限于建筑土木、智能机器人、航空航天、新材料、生物力学等领域[1-5].

对于这类自平衡结构体系而言，结构的初始几何构形是满足自平衡的必要条件，因此找形问题是对张拉整体结构分析过程中的关键问题，也是后续进行结构模态分析，静、动力分析，能量分析等问题研究的基础.张拉整体结构的初始刚度由预应力提供，其找形的过程同时也是一个找力的过程.设计人员在此阶段即可同时实现两个目标：几何形状和自应力.找形方法可以优先考虑几何方面的要求(以形状参数为主变量)，也可以优先考虑力学方面的要求(以应力参数为主变量)，但截至目前，尚无有效的方法可以仅凭其中一方面的要求实现找形的目的.通常而言，为了确定张拉整体结构初始形态，必须同时用到形状参数和应力参数.根据找形过程中所采用的主变量的不同，找形方法大致分为如下两大类：形状控制法和内力控制法.在张拉整体结构发展初期，大部分学者主要致力于研究前者，以Snelson为代表的大多雕塑家，能很好地阐述“形状控制法”的基本原理.但由于其目标仅仅是制作基于张拉整体思想的模型，因此并不在乎构件是否具有相同的规格尺寸，同时也不在乎制作完成的结构本身的力学特性，而其稳定性则是通过基于实验的试错法及实际的试错过程保证的，因此形状控制法在某些有时候的确可以得到较为理想和精确的结果.力密度法[6]、动力松弛法[7-8]、能量法[9]一系列较为熟知的典型算法则是在此基础上演化而来的.第二类找形方法则是建立在理论模型基础上的，其得到的结果通常可以满足体系本身的受力需求.但采用这类方法时，必须同时确定结构体系的拓扑形态方式及预应力水平，因此，往往可以获得精确的结果，但有时候也存在找形失败的情况.此类找形代表算法主要以搜索可行预应力或预应力优化算法的形式出现[10-15].需要强调的是，应用此类找形方法得到的结果通常是高度规则的几何构型，而不会像试错法那样出现千变万化的结果.此外，绝大部分现有的方法都必须事先假定一系列人为设定的条件，比如为简化找形步骤事先假定结构的几何对称性，或者将力密度系数作为一个符号变量，或者在事先假定若干个单元的初始长度.然而，这些信息在找形之前其实并不都是那么容易确定的.此外，大部分找形方法只适用于寻找自由形态的张拉整体结构的初始构型，而对于有边界支撑的张拉整体结构的初始自应力的确定尚缺乏有效手段.

针对上述问题，本章中提出了一种新颖的数值分析方法，目的在于解决有边界支撑的非自由形态的张拉整体结构初始自应力问题.相比于目前其他现有的方法，文本所呈现的方法，不需要事先假定或定义结构的节点坐标、材料特性、单元长度、对称性或者力密度矩阵的半正定性等一系列附加条件，只需要已知研究对象的空间维数、拓扑关系(即节点和单元的关联性)和单元类型即可，这也正是本文所提出的方法的创新之处.

1 张拉整体结构找形理论基础

1.1 基本假定

针对张拉整体结构的固有特点，找形之前有如下基本假定：

·节点坐标是确定结构几何形态的唯一途径；

·各单元之间的连接方式均为铰接；

·忽略结构自重且不考虑外荷载；

·不考虑结构的局部或整体失稳.

1.2 自平衡方程的推导

考虑一个由b个单元，n个自由节点，nf个固定节点组成的d维(d=2或3)张拉整体结构，其拓扑关系可用关联矩阵CS∈Rb×(n+nf)来表示[16].假定i和j(i

(1)

图1 (a) 二维预应力张拉整体结构 (b) 利用虚拟单元移除支撑后的自平衡张拉整体结构

将关联矩阵CS如下分解为两部分：

(2)

式中，C∈Rb×n为单元与自由节点之间的关联矩阵；Cf∈Rb×nf为单元与固定节点之间的关联矩阵.

如图1(a)为一个简单的二维预应力张拉整体结构，包括5个单元(b=5，4根拉索和1根压杆)和4个节点(2个自由节点，n=2；2个固定节点，nf=2).其关联矩阵如表1所示.

表1 二维张拉整体结构的关联矩阵

Element/NodeCsCCf1234(1)10-10(2)100-1(3)010-1(4)01-10(5)1-100

令x，y，z(∈Rn)和xf，yf，zf(∈Rnf)分别代表自由节点和固定节点在x-，y-和z-方向的节点坐标向量.同时定义q={q1,q2,…,qb}T∈Rb，矩阵中各元素之值为该单元的力fk和单元长度lk的比值，qk=fk/lk(k=1,2…b).因此可将结构的力密度矩阵Q∈Rb×b如下表述：

Q=diag(q)

(3)

于是有铰接结构的平衡方程如下[16]：

CTQCx+CTQCfxf=px

(4)

CTQCy+CTQCfyf=py

(5)

CTQCz+CTQCfzf=pz

(6)

式中，px，py，pz(∈Rn)分别代表x-，y-和z-方向作用在自由节点上的外力向量.

定义E∈Rn×n和Ef∈Rn×nf：

E=CTQC

(7)

Ef=CTQCf

(8)

由式(7)和式(8)可知，对于任意一个张拉整体结构，矩阵E的任意一列或一行的元素值均为0，因此，矩阵E是一个对称的奇异方阵[16].

忽略外力和自重，式(4)-式(8)可以写成如下形式：

Ex=-Efxf

(9)

Ey=-Efyf

(10)

Ez=-Efzf

(11)

通常情况下，固定节点的坐标xf，yf，zf是已知的，这表明可以通过式(9)-式(11)求解未知节点坐标x，y，z.换言之，矩阵E和Ef是恒定的.

根据力密度矩阵E是否为奇异矩阵，可分为如下两种情况.(1)E为非奇异矩阵：索网结构，结构中所有单元均受拉，qk>0(k=1,2…b).(2)E为奇异矩阵：预应力张拉整体结构，由于结构中有压杆(qk<0)存在，导致矩阵E必然秩亏.

在未给定固定节点坐标和力密度矢量的前提下，为了能够执行本方法，同时获得节点坐标和力密度系数这两组参数，这里引入虚拟杆单元来释放固定节点.由此，未知的固定节点可被视为自由节点.但值得注意的是，根据Zhang等[17]的建议，虚拟杆单元和固定单元还是有着本质上的不同，连接固定单元的固定节点坐标是事先指定好的，而连接虚拟杆单元的节点坐标在本方法中是未知的.利用虚拟杆单元这个概念，可以将有支撑预应力张拉整体结构转换为无外界支撑的自平衡结构.

如图1(b)所示，其中粗、细和虚线分别代表压杆、拉索和虚拟杆单元.利用虚拟杆单元6释放固定节点3和4.得到二维自平衡张拉整体结构的关联矩阵C(6×4)，如表2所示.

表2 利用虚拟单元求得的二维自平衡张拉整体结构的关联矩阵

考虑到自适应结构体系的特点，在忽略自重及外荷载的情况下，整个结构体系可视为一个空间自由体系，其几何形状可由相应的节点坐标体现[19].于是可将式(2)-式(6)简写为：

Ex=0

(12)

Ey=0

(13)

Ez=0

(14)

为简化起见,可将式(12)-(14)合并整合为：

(15)

将式(7)和式(8)带入式(12)-式(14)，合并可得张拉整体结构的自平衡方程：

(16)

式中，A∈Rdn×b为张拉整体结构的平衡矩阵.

1.3 秩亏条件

综上所述，矩阵E和A分别代表结构的力密度矩阵和平衡矩阵.于是对于一个d维张拉整体结构，有如下两个必要但不充分的秩亏条件[16].

一是关于半正定矩阵E

nE=n-rE=n-rank(E)≥d+1

(17)

式中nE为矩阵E的秩亏.

二关于平衡矩阵A，同时也是式(16)无解的必要条件

rA=rank(A)

(18)

该秩亏条件能够保证结构的自应力模态数s=b-rA≥1.独立机构位移模态数[18]为m=dn-rA.

需要指出的是，考虑到张拉整体结构是典型的自平衡机构体系，其独立机构位移模态数可分解为刚体位移模态数和无穷小机构位移模态数两部分.其中无穷小机构位移模态数mim可按下式得到：

mim=m-rb

(19)

(20)

式中，m为独立机构位移模态数；rb为刚体位移模态数；d为张拉整体结构的空间维数.

2 找形方法

2.1 力密度矩阵特征值分解

将力密度矩阵E进行特征值分解，如下：

E=ΗΛΗT

(21)

HHT=In

(22)

式中，

H∈Rn×n为正交矩阵，并且其第i列hi∈Rn为矩阵E的特征向量基；

In∈Rn×n为单位矩阵；

Λ∈Rn×n为对角矩阵，其对角线上的元素即为对应的特征值，亦即Λii=λi；

矩阵H的特征向量hi则对应于矩阵Λ的特征值λi.

力密度矩阵E特征值为零的个数等于其零空间的维数[19].此时定义k为矩阵E特征值小于或等于零的个数，则有如下两种情况：

(1)平衡态(k

(23)

(24)

Chi=0

(25)

或

det|ld(ld)T|=0

(26)

式中，ld为b个单元的长度向量，

(27)

(28)

由各相连节点(i,j)组成的单元延n方向的投影长度L∈Rb×n为：

(29)

类似的，需剔除满足式(25)或式(26)的矢量hi.

综上所述，在空间张拉整体结构的初始找形过程中，最有效的方案即为在前四组特征向量基(分别对应于前四个最小特征值)中选取三组可能的特征向量，然后将这些选取的特征值通过本问题出的方法逐步调整直至其缩小为零，由此求出对应的节点坐标矩阵：

(30)

2.2 平衡矩阵的奇异值分解

一旦相应的节点坐标被确定，即可将其代入式(16)，将平衡矩阵A进行奇异值分解：

A=UVWT

(31)

(32)

(33)

μ1≥μ2≥…≥μb≥0

(34)

式中，U∈Rdn×dn和W∈Rb×b均为正交矩阵；V∈Rdn×b为由矩阵A的非负奇异值组成的对角矩阵，其奇异值大小按式(34)排列.

在整个迭代步骤中，s同样存在如下两种情况.

(1)平衡态(s=1)

在这种情况下，矩阵U和矩阵W的零空间为：

(35)

(36)

式中，m∈Rdn为m(=dn-rA)个无穷小机构所组成的矢量.

定义机构矩阵M∈Rdn×(dn-rA)如下：

(37)

考虑到s=b-rA=1，式(36)可重新写成：

(38)

若sign(q1)≡sign(q0)，即矢量q1与q0的符号关系相同，则可判定矢量q1为满足式(16)的唯一自应力模态.

(2)非平衡态(s=0)

此时结构内部不存在自应力模态，亦即式(31)中所定义的矩阵A不存在零空间，表明矩阵V中A的右奇异值μb不为零，由此式(16)不存在非零解.此时若将矩阵W的右奇异矢量基wb(对应于矩阵V中最小的右奇异值μb)作为近似的q，则矢量q1与q0不一定会有相同的符号关系.为了克服这个困难，采用本文所提出的方法时，设计人员必须逐一检查矩阵W中的所有列，直到寻找到某个wj(j=b,b-1,…,1)，使其满足该wj中的各元素与q0有着相同的符号关系，即sign(wj)≡sign(q0).通过上述步骤，整个找形过程可以寻找出合适的q，使得下式满足：

Aq≈0

(39)

(40)

(41)

需要指出的是，对于某个给定的张拉整体结构，根据本方法所求得的力密度系数矢量q不一定是唯一的.换言之，对于相同的关联矩阵C和同样符号的q0，可能存在其他的矢量q满足条件.

2.3 稳定性计算

张拉整体结构的切线刚度矩阵KT表达式如下[21]：

⊗E

(42)

式中，

KE为结构的线刚度矩阵；

KG为与自应力模态有关的几何刚度矩阵；

e为各单元的杨氏模量；

a为各单元的横截面面积；

l0为各单元的初始长度；

I∈R3×3为单位矩阵；

⊗为张量积.

如若切线刚度矩阵为正定，则可判定整个结构体系是稳定的(在刚体位移忽略的情况下).亦即对于任意一个无穷小位移d，矩阵KT的二次型是正定的[21]

dT(KT)d>0

(43)

或，

式中，rb为式(20)中定义的独立的刚体机构数目.

利用这个准则，任意自适应张拉整体结构的稳定性可通过结构切线刚度矩阵的特征值来判定.

定义的不平衡力矢量υf∈Rdn如下[11]

υf=Aq

(45)

或将不平衡力分别在x-，y-和z-方向上投影

υx=Ex

(46)

υy=Ey

(47)

υz=Ez

(48)

引入Euclidean范数,设计误差κ可以表述为

(49)

3 算例

基于第二节中的理论可知，根据事先给定的不同信息，存在两种不同的张拉整体结构找形方法.第一种方法是通过给定结构的几何形态找形，即通过给定结构的节点坐标求解各单元的力密度矢量.第二种方法是根据有限的节点拓扑关系(关联矩阵)和单元类型，同时寻找出合适的节点坐标和力密度矢量，从而实现找形的目的.在本节中，分别将这两种方法称为几何形态法和并行法.

3.1 截顶四面体

考虑一个由二杆六索组成的二维张拉整体结构(一类)，如图2(a)所示.通过引入虚拟杆单元将两个支撑转换为自由节点从而得到自平衡张拉整体结构.当实施了第二节中所述找形方法后，虚拟杆单元将重新转换为相应的两个支撑.加入虚拟杆单元后的自平衡二维张拉整体结构拓扑形式如图2(b)所示，其中粗、细和虚线分别代表压杆、拉索和虚拟杆单元.

图2 (a) 有支撑的二维二杆张拉整体结构 (b) 利用虚拟杆单元移除支撑后的自平衡二维张拉整体结构

几何形态法.根据本方法所需条件，将事先指定的结构节点坐标列于表3中.

表3 二维二杆张拉整体结构节点坐标

节点123456x2.005.005.002.000.007.00y-2.00-2.002.002.000.000.00

如表4列出了最终得到的归一化力密度系数(表示相对于杆1的力密度，下同)，与Zhang等[21]的结果相符.需要指出的是，Zhang等所采用的方法，将某些单元的方向指定为已包含在结构自平衡方程中的几何约束方向.随后，通过指定两组独立的力和节点坐标，从自平衡约束方程中按顺序求解出特定的力和节点坐标.这也是与本方法的主要区别所在.

表4 二维二杆张拉整体结构力密度系数

力密度q1q2q3q4q5q6q7q8q9Zhang等[64]1.00000.66671.00001.00000.66671.0000-0.5000-0.5000-0.5714本文1.00000.66671.00001.00000.66671.0000-0.5000-0.5000-0.5714

图3 二维二杆张拉整体结构的稳定形态

并行法.节点坐标、几何对称性、单元长度和力密度系数均未知，所需的条件仅为结构的关联矩阵和单元类型.根据本章所提出的找形方法，自动分配初始力密度矢量如下：

找形得到的归一化力密度系数矢量为：

如图3展示了找形结束并且移除虚拟杆单元后结构的最终稳定形态.找形过程仅经过一次迭代即达到收敛，且设计误差κ=3.4667×10-16

两种找形方法中力密度矩阵E均为半正定矩阵，这表明在不考虑材料特性和预应力的影响下，结构是超稳定的.进一步而言，自应力模态刚化了所有的无穷小机构位移模态，使得整个结构体系在除了三个刚体方向之外均能保持稳定.总而言之，利用本章提供的找形方法，在已知结构节点拓扑关系(关联矩阵)和单元类型的前提下，通过强制满足两个必要的秩亏条件，能够找到几何稳定的自平衡张拉整体结构.

3.2 三维五杆张拉整体结构

考虑一个由5杆20索组成的三维张拉整体结构(一类)，如图4(a)所示.加入虚拟杆单元后的自平衡三维张拉整体结构拓扑形式如图4(b)所示.

图4 (a) 有支撑的三维五杆张拉整体结构 (b) 利用虚拟杆单元移除支撑后的自平衡三维张拉整体结构

表5 三维五杆张拉整体结构节点坐标

Node123456789101112131415x-9.51-5.889.519.510.00-4.76-4.76-2.94-2.942.942.944.764.760.000.00y3.09-8.09-8.093.0910.001.541.54-4.05-4.05-4.05-4.051.551.555.005.00z0.000.000.000.000.002.50-2.502.50-2.502.50-2.502.50-2.502.50-2.50

几何形态法.同二维二杆张拉整体结构，将事先指定的结构节点坐标列于表5中.最终得到的归一化力密度系数矢量为：

并行法.同二维二杆张拉整体结构，已知信息仅为结构的关联矩阵和单元类型.根据并行找形法，自动分配初始力密度矢量如下：

找形得到的归一化力密度系数矢量为：

找形过程仅经过一次迭代即达到收敛，且设计误差为κ=4.0715×10-16满足条件.图5展示了找形结束并且移除虚拟杆单元后结构的最终稳定形态.在忽略六个刚体位移模态(rb=6)后,最终得到一个自应力模态(s=1)和十个无穷小机构位移模态(mim=10).本例中力密度矩阵E为半负定矩阵，表明结构不是超稳定的.通过研究该结构的切线刚度矩阵可以发现其为正定，从而证明该结构是几何稳定的.

图5 三维五杆张拉整体结构稳定形态(a)俯视图(b)轴测图

4 结论

本文主要研究了张拉整体结构初始找形问题，主要在力密度法的基础上，引入几何拓扑理论和矩阵分析理论，提出了考虑结构几何拓扑、稳定性、误差评估等诸多因素的张拉整体结构的初始找形方法.本方法的提出将张拉整体结构的力密度矩阵和平衡矩阵的关系紧密结合在一起，并通过结构的切线刚度矩阵特征进一步判定体系的稳定类型.本文的方法均由编程实现，不需要反复试算，大幅提升了找形的计算效率.最后结合两个张拉整体结构的具体实例进行分析比较，结果表明本文所提出的找形方法对于分析解决带有约束支撑的张拉整体结构是非常适用的.