一种空间六杆张拉整体结构找形分析与设计

2020-08-03冯晓东周倩倩章万鹏张佳丹赵容舟

西安建筑科技大学学报（自然科学版） 2020年3期

冯晓东，周倩倩，章万鹏，何溯，张佳丹，赵容舟

(1.绍兴文理学院土木工程学院，浙江绍兴 312000 2.浙江大学建筑工程学院，浙江杭州 310000)

“张拉整体”这一概念由Snelson和Fuller在20世纪60年代提出，该类结构体系一般是由若干离散的杆元和索元组成，其中杆元受压，索元受拉.结构中的受拉索元存在着自应力，在未受外力作用情况下，该结构可起到自适应平衡作用.由于这类智能结构存在着一系列的潜在优点，使其自出现以后，“张拉整体”概念便在众多工程和科学领域得到了广泛的应用.在生物力学领域中研究细胞骨架力学特性时通常运用张拉整体模型，并将其视为合适的研究模型.但现有的张拉整体模型并不能够较好地解决这些几何和拓扑形态极为复杂的聚合物纤维结构.因此复杂且非对称的张拉整体结构自由形态设计仍需进行深入的研究.

近年来，针对张拉整体结构研究所得的找形方法主要有:以多面体几何分析为基础发展的几何分析法、力密度法及其他系列关于张拉整体结构找形的衍生算法[4-5].此外，还有将结构找形问题转化为约束优化问题的非线性规划法；在原有找形基础上的改进力密度法；利用有限单元法确定结构初始形态[9]；张拉整体结构找形中蚁群算法的应用[4]；多自应力模态张拉整体结构的数值找形方法；基于平衡方程和几何相容方程的形式解决找形问题[13]等.

现有找形方法大多基于一定的条件假定，其中如以力密度系数作为符号变量并需假定各构件长度[14]；通过假定结构对称以达到缩减矩阵简化找形步骤的目的[15].由于事先确定结构的初始假定条件存在难度.因此，对于快速求解复杂的、非对称的张拉整体结构的找形问题仍存在一定的挑战.

本文为解决一类张拉整体结构的找形问题，提出了一种新型数值分析找形方法.分别采用谱分解和奇异值分解法处理结构体系的平衡矩阵和力密度矩阵.根据力密度矩阵和平衡矩阵的最小秩亏条件，采用双向循环迭代的方式，以相当小数量的迭代即可获得任意稳定态，求解节点坐标和力密度，即:满足于切线刚度矩阵为正定矩阵或超稳定态.同时，张拉整体结构自由形态的几何刚度矩阵必为正定矩阵.现有的找形方法往往需要假定张拉整体结构的单元长度、初始节点坐标、力密度矩阵的半正定性或几何对称性等条件，但本文所呈现的方法中，只根据已有结构的空间维数、单元类型以及节点间的几何拓扑关系，就可完成找形，本文较其他找形方法具有明显的优势.

1 力密度找形法及相应的秩亏条件

1.1 基本假定

张拉整体结构找形需要基于以下六条基本假定：

(1)已知结构的几何拓扑连接，其形态由节点坐标表示；

(2)节点连接均为铰接；

(3)不考虑任何外力作用；

(4)忽略结构自重.

(5)不考虑任何整体或局部的弯曲变形；

(6)不考虑耗散力对结构体系的影响.

1.2 张拉整体结构的自平衡方程

对任意张拉整体结构，关联矩阵CS∈Rb×(n+nf)可用于表示结构的连接关系.其中,b为单元数量，d为结构空间维度.根据节点的约束情况，可分为自由节点和固定节点，分别用n和nf来表示.

基于上述定义，假设连接单元k两端的节点编号分别为i和j(i

(1)

根据节点的连接关系，CS可分为C和Cf，分别表示自由节点间和固定节点间的连接关系.考虑这两者在数字标记顺序上的优先关系(固定节点可优于自由节点)，故关联矩阵CS可表示为

(2)

对任意单元可根据静力平衡关系建立如下平衡方程：

(3)

(4)

(5)

式中：x,y,z为节点在各方向上的坐标；t为对应单元内力；lk为对应单元长度；px,py,pz(∈Rn)分别为作用在自由节点各方向上的外力.

定义上式qk=t/lk(k=1,2…b)为力密度，且矩阵表示为q={q1,q2,…,qb}T∈Rb.故整个结构体系的力密度矩阵Q∈Rb×b可表示成

Q=diag(q)

(6)

且任意节点在x,y,z三个方向上需满足平衡方程组：

CTQCx+CTQCfxf=px

(7)

CTQCy+CTQCfyf=py

(8)

CTQCz+CTQCfzf=pz

(9)

因为张拉整体结构本身具备无固定节点的特征，故在忽略外力和自重的条件下，整个结构体系可视为一个空间自由体系，可根据对应的节点坐标来确定体系的几何形状，且有

px=0

(10)

py=0

(11)

pz=0

(12)

选用矩阵E∈Rn×n和Ef∈Rn×nf简化公式，

E=CTQC

(13)

Ef=CTQCf

(14)

根据式可得，矩阵任意列或行元素之和均为0，因此对任意一个张拉整体结构而言，矩阵E必然为一个对称的奇异方阵[16]，表示张拉整体结构的力密度矩阵.

方程组(7)～(9)和方程组(10)～(14)可表示为：

Ex=0

(15)

Ey=0

(16)

Ez=0

(17)

继续对方程组～合并与简化，得

E[x,y,z]=CTQC[x，y，z]=[0，0，0]

(18)

式中:[x，y，z]为张拉整体结构的节点坐标矩阵.

将式(13)代入式(15)～(17)，整合后可得张拉整体结构的自平衡方程：

(19)

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

1.3 两个必要的秩亏条件

由于d维张拉整体结构至少存在一种可行自应力且需满足以下秩亏条件[16]来保证结构自应力模态数s=b-rA≥1和独立机构位移模态数m=dn-rA[17]

秩亏条件一，考虑半正定矩阵E，即

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

(20)

式中：nE为矩阵E的秩亏.

秩亏条件二，考虑平衡矩阵A，即

s=rA=rank(A)

(21)

此外，式(12)为式(10)无解的必要条件.

其中独立机构位移模态数包含刚体位移模态数和无穷小机构位移模态数两部分[18].无穷小机构位移模态数mim可由以下式子得到.

mim=m-rb

(22)

(23)

式中：m为独立机构位移模态数；rb为刚体位移模态数.

2 基于矩阵分析的找形步骤

本文提出了一种避免事先假定(单元长度，几何对称性等)的找形方法，但需已知结构的空间维数、拉压单元布置情况以及关联矩阵，便能通过算法得到找形结果[19].

由已知的结构拉压关系，如下定义矢量：

(24)

图1 自由形态张拉整体结构找形流程图Fig.1 Process of the proposed form-finding procedure for free-form tensegrity structures

(25)

(26)

2.1 力密度矩阵的特征值分解

将力密度矩阵特征值分解，得

E=ΗΛΗT

(27)

HHT=In

(28)

式中：H∈Rn×n为正交矩阵，hi∈Rn为矩阵E第i列的特征向量基，

Λ∈Rn×n其对角线上的元素λi为hi对应的特征值；

In∈Rn×n为单位矩阵；

根据上述分解结果得，矩阵E的零特征值个数与零空间的维数一致[20].此时，定义Λ中对角线元素(特征值)不大于零的个数为k，存在以下两种情况

(29)

(30)

Chi=0

(31)

det|ld(ld)T|=0

(32)

(33)

式中：ld为d维空间中b个单元的长度.

当矩阵E不满足式子时，结构体系则处于非自平衡状态.在这种状态下，整个迭代过程将迫使矩阵的秩亏增大使其能够满足最小秩亏条件，且矩阵E并非定为半正定矩阵.因此需要引入张切线刚度矩阵KT去判别拉整体结构的稳定性，其表达式如下所示

(34)

式中：KE为弹性刚度矩阵；KG为几何刚度矩阵；e为各单元材料的杨氏模量(矢量)；a为各单元材料的横截面面积(矢量)；l0为各单元的初始长度(矢量)；I∈R3×3为单位矩阵；⊗为张量积.

根据式子确定忽略刚体位移数量.当切线刚度矩阵处于正定状态，则整个结构体系将是稳定的.

dT(KT)d>0

(35)

(36)

式中：rb为式中定义的刚体位移模态数.

结合上述两种情况可知，正交矩阵H前四组特征向量基(分别相应于前四个最小的特征值)中选择三组可能为张拉整体结构的节点坐标，再通过循环迭代的方式，使得特征值逐步递减为零，使得结构体系最终能够处于平衡状态.

(37)

2.2 平衡矩阵的奇异值分解

对平衡矩阵A奇异值分解可得：

A=UVWT

(38)

(39)

(40)

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

(41)

式中：U∈Rdn×dn和W∈Rb×b均为正交矩阵；

V∈Rdn×b为由矩阵A的非负奇异值组成的对角矩阵.

在算法计算流程中，对于自应力模态数s可分为两种情况.

(1)平衡态(s=1)

通过平衡矩阵的奇异值分解可得到力密度和机构的矢量空间基.因此，矩阵U和矩阵W应有如下形式的零空间：

(42)

(43)

式中：m∈Rdn为m个无穷小机构所组成的量.

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

(44)

由于s=b-rA=1，式(34)可写成如下形式：

(45)

若矢量q1与q0符号相同，即sign(q1)≡sign(q0)，则矢量q1为满足齐次方程的唯一自应力模态.

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

在结构自应力模态数为0的情况下，q不存在非零解.此时，若将矩阵W的右奇异矢量基wb(对应于矩阵V中最小的右奇异值μb)作为近似的q，则矢量q1与q0可能会存在不同符号.为了解决该问题，必须逐一尝试矩阵W所有列，找到某个wj(j=b,b-1,…,1)的所有元素都与q0符号相同，即sign(wj)≡sign(q0).此时，近似的q可用矢量wj直接使用.经过上述步骤，所得到的q，使得下式满足

Aq≈0

(46)

综上所述，为满足式(25)和式(26)所述的两个秩亏条件，通过式和式的循环迭代，得到最终的结构形态及自应力.值得注意的是，对于给定的张拉整体结构，可先判断结构的自应力模态数s，若为1，则本方法所求得的力密度系数矢量q是唯一的.若大于1，则可能存在满足条件的多种矢量q，而本方法所求得的力密度系数矢量q则是一种可行解.

2.3 设计误差评估

鉴于张拉整体结构为自平衡结构体系，故采用不平衡力矢量υf∈Rdn作为评估计算结果精度的指标，其定义如下.

υf=Aq

(47)

在各个方向上不平衡力矢量υf分别表示为：

υx=Ex

(48)

υy=Ey

(49)

υz=Ez

(50)

引入Euclidean范数，则设计误差κ可以表示为

(51)

文中，控制设计误差κ在10-10以内.

3 算例

本节主要是通过一类张拉整体结构数值算例验证该找形方法的有效程度和适用程度.根据算例结果表明，本文所研究的找形方法在应用于分析解决稳定形态和超稳定形态下的张拉整体结构时，非常适用.

3.1 6杆24索张拉整体结构

以空间6杆24索张拉整体结构为例，求解该结构稳定状态下的结构形态及自应力.该结构共有30个单元，12个节点，每个节点上均连接4根索和1根杆.结构的拓扑连接及构件拉压关系如图2所示.其中编号25～30为杆单元，其余为索单元.

图2 6杆24索张拉整体结构拓扑连接图Fig.2 Topological connection of tensegrity structure with 6 struts and 24 cables

图3 无指定节点坐标下得到的6杆24索张拉整体结构的几何形态 Fig.3 The obtained geometric configuration of tensegrity structure (6 struts and 24 cables) without specified node coordinates

根据本论文所研究的找形方法，自动分配初始力密度矢量如下：

此结构包含1个自应力模态和1个内部机构位移模态，通过乘积力的正定性，以验证此结构体系为静不定且动不定结构，几何稳定.如图3，展示了本结构在无指定节点坐标情况下的几何形态.

指定张拉整体部分节点坐标：(0，0，0)、(0，1，0)、(2，0，0)、(4，3，0)、(3，2，1).根据本文算法，最终结构经12次迭代后收敛，找到满足限定条件下的张拉整体结构如图4，节点坐标见表1，并在表2中列出此结构归一化后的可行自应力.

图4 找形得到的6杆24索张拉整体结构(状态1)Fig.4 The obtained tensegrity structure with 6 struts and 24 cables (configuration 1)

表1 图4所示结构的各节点坐标

在相同的初始条件下，仅仅改变部分指定节点，具体更改如下(1)( 0，0，0)、(3)( 5，0，-1)、(6)( 0，3，0)、(9)( 2，2，5)、(12)( 6，2，1).同样以图1中的计算流程对新节点坐标下的张拉整体结构进行找形计算.最终得到的张拉整体结构拓扑形态如图5所示，节点坐标与构件自应力如表3和表4所示.

表2 图4所示结构各构件的自应力

图5 找形得到的6杆24索张拉整体结构(状态2)Fig.5 The obtained tensegrity structure with 6 struts and 24 cables (configuration 2)

表3 图5所示结构的各节点坐标

三维6杆24索张拉整体结构在有无事先指定节点的两种情况条件下得到的找形结果表明，该算法在无指定节点情况下，找形结果可接近于经典扩展八面体张拉整体，构件受力也较为均匀.而在指定部分节点情况下，仍可找到一种稳定的张拉整体结构形态.两者找形结果不同是由于自应力模态数为1的张拉整体结构形态与结构的节点坐标为一一对应的关系.即指定节点坐标找形结果定会出现特定的张拉整体结构与其对应[23].

表4 图5所示结构各构件的自应力

4 张拉整体结构模型制作

在现有成果中，对结构找形的研究基本处于数值计算阶段，制作实际模型并进行验证的研究相对较少.本节将基于算例中的找形结果，制作三维六杆24索张拉整体结构，并将所得模型与理论计算结果进行比较.同时，尝试改变初始模型的节点位置，比较改变前后模型的状态.结合模型制作期间存在的技术问题，对结果进行初步分析总结.

4.1 模型制作过程

此模型含有6个杆单元，24个索单元，12个节点.根据算例的找形结果，其拉索单元的力密度、长度及内力均相等.由此确定结构单元尺寸：杆单元受压，选择长度为0.5 m，直径为10 mm的PVC管；索单元受拉，选择长度为0.7 m的棉绞线.对结构的计算模型进行力学模型简化，完成初始模型制作，具体流程如下：

(1)首先剪裁六段长度为0.5 m的PVC管，在每段PVC管两端，用铅笔画出深15 mm、宽3 mm且位于同一直径方向上的黑实线，再使用钢锯、刻刀、剪刀等工具分别在PVC管两端，沿黑实线刻出凹槽，并进行适当打磨使其光滑.

(2)裁剪六段长度为0.7 m的棉绞线，并在每段棉绞线中点处用黑色签字笔做上记号.以棉绞线中点为基准，卡入PVC管一端的凹槽中，将棉线两端拉直至另一端凹槽，每边预留长度为5 mm，并使用胶带进行固定，如图6所示.

(3)用上述方法制作完成六组杆件，完成准备工作，采用临时支撑将6根压杆进行临时固定.同时，为方便穿索过程,任意选取两根压杆采用弹性固定.

(4)将各个杆端按照x，y，z的坐标方向顺序依次进行连接，模型具备初步形态时，调整每个索段位置并拉紧棉线，以对6杆24索张拉整体结构模型施加一定的力.在结构模型具有一定刚度时，拆除临时支撑装置.

图6 安装杆件的前期准备Fig.6 Preparation for installation of members

图7 张拉成型后的六杆张拉整体结构Fig.7 The formed configuration of tensegrity structure with 6 struts

自应力的产生机理：通过对索单元施加预应力，使索拉伸改变节点坐标.同时导致其他构件单元的产生拉伸或者压缩，进而达到施加预应力的效果.在施加过程中，若索单元受压则会松弛.因此，在对索单元拉伸时，需要保证其他索单元仍处于受拉紧绷的状态.基于此原则及自应力的产生机理，按上述制作方法后调试不同的索单元，最终可得到6杆24索张拉整体结构模型，如图7所示.

4.2 不同约束条件下的结构形态分析

基于图7所示初始模型，对上述制作流程进行调整.即在第3步和第4步进行自动构型时，任意选取五个节点对其固定，通过调节其他剩余索杆单元，使其在调整后仍能保持为自稳定状态，并观察模型在此时的结构形态.

经过初步实验，总结出下述规律：

(1)当束的节点包含两根平行杆段时，稳定后各节点沿杆段方向的位移变化较大，整体形态较初始构型呈压缩或拉伸状态.

(2)当约束的节点在空间内相对分散时，稳定后大部分节点的移动相对情况(1)较小，但个别节点的移动较大，整体呈现不规则状.

4.3 原因分析与存在问题

对比理论解，结合制作过程中遇到的技术障碍，对实验结果进行初步的分析和总结，并提出下一阶段可深入进行研究的内容，具体如下：

(1)对于初始结构，其形态在杆件位置和棉线的张拉状态进行略微调整后，与理论结果较为接近.误差主要来源考虑是实验材料，工具以及制作流程存在一定缺陷.具体包括：剪裁处的光滑度和水平度不易保证；杆件两端凹槽的方向存在一定偏角；棉绞线的质量参差不一，力学性能指标无法得到准确保证；棉线的拉伸过程只能在多次尝试后凭借经验进行，没有合适的定量工具进行控制.

(2)对于存在约束的结构，不同的约束条件会产生不同结果，与理论结果也存在一定差异，主要差异考虑在对节点进行固定时，位置定位不够精准.对于第一种情况，由于固定节点包含两根同方向杆段，固定位置相对容易且准确，最终形态与理论结果较为接近.第二种情况由于涉及固定的杆段数量较多，节点位置的控制相对繁琐，易发生偏移，最终结果的成型率相对较低，规律性不够明显.