索桁张拉结构是指由单榀索桁架通过一定布置原则形成的张拉整体结构，这类结构的共同点是整体结构可拆分为单榀索桁架。又由于索桁张拉结构具有独特的拓扑关系及自重轻、传力路径清晰等优点，因此在实际工程越来越被设计者所青睐。目前，索桁张拉结构大致可分为轮辐式索桁结构^[1]和环形交叉索桁结构^[2]两类，分别如图1和2所示。轮辐式索桁结构是由一系列沿径向分布的单榀索桁架组成，环形交叉索桁结构是由一系列单榀索桁架互相交叉编织构成。在没有预张力作用下，索桁张拉结构是一种机构，不具有特定形状且没有承载能力；对结构施加预张力后，形成具有特定形状和一定刚度的结构。由此引出本研究的两个问题：合理形状判断和独立自应力模态求解。

已有研究方法大都是先设计结构的初始形状；然后从结构拓扑关系出发，寻求结构的平衡矩阵；再对平衡矩阵使用奇异值分解法（SVD），判断结构的稳定性，并求解结构的独立自应力模态。如果结构的稳定性不符合要求，利用平衡矩阵修改节点坐标，最终可求出合理的结构位形^[3-4]。有一些求解张拉结构自应力模态的方法，如动力松弛法^[5-6]、力密度法^[5]等。以上方法一般均需要借助计算机编程求解，不易被工程技术人员掌握。此外，如张爱林^[7]、Wang^[8]、董石磷^[9]等从平面索桁结构的节点出发，通过求解每个节点的平衡方程组获得索桁结构的自应力模态，但该方法需要求解一系列方程组，应用起来不方便，同时不能直接利用有限元软件进行求解。田广宇等^[10]基于有限元提出逐点去约束法求解索桁结构的可行预应力，但该方法存在需要应用生死单元技术、坐标转换、反复迭代求解确定结构的合理位形等缺点，应用时较复杂。作者基于文献[10]中的逐点去约束法提出基于合理位形改进的逐点去约束法。本文方法基于索桁张拉结构的特点，从合理位形出发推导了形状合理性的简洁判断公式。然后，利用改进的逐点去约束法求解索桁张拉结构的自应力模态。

1 形状判定策略

由于索桁张拉结构是由一系列的单榀平面索桁架按照一定的规律集合而成，所以可先把整体结构分解为平面索桁架。只要把平面索桁架的合理形状设计出来，通过逆组装过程即可形成整体结构。

从平面索桁架的合理几何形状出发，推导了维持这种合理形状的内在关系。假设如图3所示的索桁架处于平衡状态（每个节点均受力平衡），其推导具体过程如下。

以整体结构为研究对象，对节点P₀取矩可得：

${T_4} \cdot {h_4} - {N_4} - {d_4} = 0$

(1)

以上部所有节点为研究对象，并求出每个节点水平方向上的平衡关系式，经过化简可得出一般通用公式为：

${T_4} = {T_i} \cdot \cos\; {\alpha _i}$

(2)

同理，可推导出下部节点在水平方向上的一般通用公式：

${N_4} = {N_i} \cdot \cos\; {\beta _i}$

(3)

以所有竖柱为研究对象，并求出每根竖柱在竖直方向上的平衡方程，经过化简可得出一般通用公式为：

${T_i} \cdot \sin\; {\alpha _i} = {N_i} \cdot \sin\; {\beta _i}$

(4)

将式（2）和（3）代入式（4）中，可得：

${T_4} \cdot \tan\; {\alpha _i} = {N_4} \cdot \tan\; {\beta _i}$

(5)

式（5）说明具有合理形状的索桁架，其上、下弦杆件的水平向力分量相等。进一步可知，节点处的竖向力也使竖柱处于平衡状态。水平方向的平衡力和竖直方向的平衡力使索桁架处于整体平衡状态。

根据图3，由几何关系可知：

$\left\{ {\begin{aligned} & {\tan\; {\alpha _0} = {h_1}/{l_0}}, \\ & {{\tan}\;{\beta _{0}} = {d_1}/{l_0}} \end{aligned}} \right.$

(6)

$\left\{ {\begin{aligned} & {\tan\; {\alpha _1} = ({h_2} - {h_1})/{l_1}}, \\ & {{\tan}\;{\beta _{1}} = ({d_2} - {d_1})/{l_1}} \end{aligned}} \right.$

(7)

由式（6）和（7），并基于式（5）可得：

$\frac{{\tan\; {\alpha _0}}}{{\tan\; {\beta _0}}} = \frac{{{h_1}}}{{{d_1}}} = \frac{{{T_4}}}{{{N_4}}}$

(8)

$\frac{{\tan\; {\alpha _1}}}{{\tan\; {\beta _1}}} = \frac{{{h_2} - {h_1}}}{{{d_2} - {d_1}}} = \frac{{{T_4}}}{{{N_4}}}$

(9)

基于差比定理，由式（8）和（9）可得：

$\frac{{\tan\; {\alpha _1}}}{{\tan\; {\beta _1}}} = \frac{{{h_2} - {h_1}}}{{{d_2} - {d_1}}} = \frac{{{h_1}}}{{{d_1}}} = \frac{{{h_2}}}{{{d_2}}} = \frac{{\tan\; {\alpha _0}}}{{\tan\; {\beta _0}}} = {\text{常数}}C$

(10)

同理，可得一般公式为：

$\frac{{{h_1}}}{{{d_1}}} = \frac{{{h_2}}}{{{d_2}}} = \frac{{{h_3}}}{{{d_3}}} = \frac{{{h_4}}}{{{d_4}}} = \cdot \cdot \cdot = \frac{{{h_i}}}{{{d_i}}} = {\text{常数}}C$

(11)

式（11）说明具有合理形状的索桁架所有竖柱的上半部分与下半部分的比为某一固定常数C。

因此，当索桁架所有竖柱的上半部分与下半部分的比为某一固定常数C时，上、下弦杆件的水平向力分量相等，该索桁架结构处于平衡状态且存在独立自应力模态。需要说明的是，固定常数C的取值可参考文献[11]。文献[11]中给出上弦索（稳定索）的经济矢跨比为1/25～1/20，下弦索（承重索）的经济矢跨比为1/20～1/15，根据上、下弦索矢跨比可获得合理的固定常数C。

2 求解自应力模态的有限元理论

从平衡矩阵理论出发推导出基于有限元软件求解索桁张拉结构的有限元方法。根据文献[12-13]可知，平衡矩阵理论的基本方程为：

${{A}} \cdot {{t}} = {{f}}$

(12)

式中， ${{A}} $ 为平衡矩阵，t为杆件内力向量，f为外荷载产生的等效节点力荷载向量。式（12）实质上是用矩阵表示所有节点处力的平衡方程的集合。平衡矩阵 ${{A}} $ 的数学性质可用线性代数中的秩、子空间等表示，反映了结构在几何和力学上的特性。式（12）对整体结构进行分析，平衡矩阵 ${{A}} $ 为每个节点所有单元方向向量的集成，但获得 ${{A}} $ 的过程比较复杂，而有限元软件又无法直接给出矩阵 ${{A}} $ ，一般需要自己编程。同时，用平衡矩阵求解时，忽略了结构的变形和材料属性，不能真实反映结构的自应力模态。

索桁张拉结构一般以平面索桁架出现在结构中，如轮辐式结构^[1]和环形交叉索桁结构^[2]，其中索桁架是平面受力单元。索桁架的每个节点均连接一根撑杆（或连接索）和两根弦索，称其为“基本单元”（图4（a））。在每个基本单元中，若杆件的内力在外力作用下能相互平衡，并使得节点处于合理位形，那么将所有构件的内力集成到式（12）中的失量t中时，式（12）也自然成立。此时，所求杆件的内力是满足结构初始位形的一组预应力值；结构若不受外力作用，所求预应力则为结构的自应力模态。

图4（b）中，节点1在撑杆、弦索1和2的内力F₁、F₂和F₃及外力P共同作用下处于平衡状态且位于设计坐标。以图4（b）中的坐标系建立节点的平衡方程，可写成式（13）的形式：

$\left[\!\!\!\!\!\!\!\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\cos\; {\theta _1}}&{\cos\; {\theta _2}}&0 \end{array}} \\ {\begin{array}{*{20}{c}} {\sin\; {\theta _1}}&{\;\sin\; {\theta _2}}&1 \end{array}} \end{array}}\!\!\!\!\!\!\!\! \right]\left[\!\!\!\! {\begin{array}{*{20}{c}} {{F_1}} \\ {{F_2}} \\ {{F_3}} \end{array}}\!\!\!\! \right] = \left[\!\!\!\! {\begin{array}{*{20}{c}} {{P_1}} \\ {{P_2}} \end{array}} \!\!\!\! \right]$

(13)

因为求解结构的自应力模态，所以不考虑结构在外力作用下的情况。令P=0，式（13）可写成式（14）：

$\left[\!\!\!\!\!\!\!\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\cos\; {\theta _1}}&{\cos\; {\theta _2}}&0 \end{array}} \\ {\begin{array}{*{20}{c}} {\,\sin\; {\theta _1}}&\;{\sin\; {\theta _2}}&1 \end{array}} \end{array}}\!\!\!\!\!\!\!\! \right]\left[ \!\!\!\!{\begin{array}{*{20}{c}} {{F_1}} \\ {{F_2}} \\ {{F_3}} \end{array}}\!\!\!\! \right] = \left[ \!\!\!\!{\begin{array}{*{20}{c}} 0 \\ 0 \end{array}}\!\!\!\! \right]$

(14)

式（14）的系数矩阵为式（15）， ${{A}} $ 为 $2\times3 $ 阶矩阵，其秩等于2。

${{A}} = \left[\!\!\!\!\!\!\!\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\cos\; {\theta _1}}&{\cos\; {\theta _2}}&0 \end{array}} \\ {\begin{array}{*{20}{c}} {\sin\; {\theta _1}}&{\;\sin\; {\theta _2}}&1 \end{array}} \end{array}}\!\!\!\!\!\!\!\! \right]$

(15)

因此，式（14）的基础解系只有1个自由向量，可以表达成式（16）：

${\left( { {{{k{F}}_1}} , {{{k{F}}_2}} , {{{k{F}}_3}} } \right)^{\rm{T}}} = {\left( {\frac{{{\rm{cos}}\;{\theta _2}}}{{\sin ({\theta _2} - {\theta _1})}},\frac{{{\rm{cos}}\;{\theta _1}}}{{\sin ({\theta _1} - {\theta _2})}},1} \right)^{\rm{T}}}$

(16)

式中：F₁、F₂为索的拉力；F₃为撑杆的压力；k为实数，θ₁和θ₂分别为两根弦索与x轴正方向的夹角，且cos θ₁、sin θ₁、cos θ₂和sin θ₂均不为0。只要F₁、F₂和F₃按照式（16）的比例取值，均能够满足式（14），并使节点1处于平衡状态，这一组比例关系被称为基本单元的自应力模态^[14]。

应用有限元软件对图4的基本单元求解时，无法直接求解式（13）和（14），而是求解式（17）：

$({{{K}}_{\rm{E}}} + {{{K}}_{\rm{G}}}) \cdot \Delta {{U}} = {{P}} + {{R}}$

(17)

式中，K_E为结构的弹性刚度矩阵，K_G为几何刚度矩阵， $ \Delta {{U}} $ 为节点位移向量，P为节点等效荷载向量，R为求解过程中由于位移U的高阶项所产生的不平衡力向量。

求解式（17）时可用力迭代法或位移迭代法^[13]，其迭代求得的解为式（13）与（14）的解。对于图1所示的基本单元而言，位移迭代法或力迭代法^[13]的收敛速度都非常快，一般只需要迭代2次即可达到较高的精度。根据有限元软件求解过程的简易程度，建议采用力迭代法。

本文改进的逐点去约束法正是通过求解这种基本单元从而获得整体结构的预应力分布，这也是本文改进方法的基本思想。其求解策略可形象地表述为“化整为零”，把所求的各单元的内力施加到整体结构中又可称为“化零为整”。

3 基于合理位形改进的逐点去约束法求解索桁张拉结构自应力模态的过程

通过对逐点去约束法不足之处的研究，作者提出改进的逐点去约束法，解决以下两个问题并给出具体办法。

问题1：原方法中采用生死单元技术^[13]把单元内力作为集中力施加给节点。首先，应用生死单元技术时，需要用到节点坐标转换，在实际操作中不太方便。其次，生死单元技术在有限元软件中又会引起求解不收敛，这是通用有限元软件的缺陷。最后，采用生死单元技术并没有删除该单元，而是对其刚度乘以一个很小的数值，导致计算的自应力模态和理论值会有误差。尤其是当索桁架跨度大、节点较多时，如用原方法，此时节点误差不断积累，会引起求解的结果与理论值偏差较大，甚至失真。

针对问题1，改进方法不采用生死单元技术，因此避免了使用坐标转换，简化了有限元求解流程，也确保了计算的收敛性，同时消除了因使用生死单元技术而导致的累计误差，提高了计算精度。

问题2：原方法先求解结构的自应力模态，再判断结构是否存在自应力模态。当结构不存在自应力模态时，先通过修改节点坐标，再利用力迭代法或位移迭代法进行迭代求解，直到修改后的位形使结构存在独立自应力模态。由于需要反复修正节点坐标并迭代求解，当设计位形与合理位形相差较大时，原方法的计算量将会过大。

针对问题2，根据第1节内容可知，只要满足索桁架所有竖柱的上半部分与下半部分之比为某一固定常数C时，结构处于平衡状态且存在自应力模态。因此，使用式（11）对结构形状先判定再求解其自应力模态，这样避免了原方法的麻烦，极大地简化了求解过程。

基于ANSYS有限元软件，以图5为例说明了基于合理位形改进的逐点去约束法（本文方法）求解索桁张拉结构自应力模态过程。

具体步骤如下：

1）利用式（11）判断结构是否存在自应力模态。如果结构不存在自应力模态，利用式（11）调整结构的几何位形，使调整后的结构存在自应力模态。

2）如图5所示，先约束所有节点，仅释放节点5处的约束。对单元E4、E5和E14施加初始内力，通过位移迭代或力迭代可求出当节点5处于设计坐标时单元E4、E5和E14的内力为 $ F_{{\rm{E}}4}^1 $ 、 $ F_{{\rm{E}}5}^1 $ 和 $ F_{{\rm{E}}14}^1 $ （上标表示第1次去掉约束并求解，下标表示单元编号，下同）。

3）恢复节点5处的约束，同时释放节点4处的约束。把 $ F_{{\rm{E}}4}^1 $ 、 $ F_{{\rm{E}}5}^1 $ 和 $ F_{{\rm{E}}14}^1 $ 作为已知内力分别赋予单元E4、E3和E13，通过计算可求出单元E4、E3和E13的内力为 $F_{{\rm{E}}4}^2$ 、 $F_{{\rm{E}}3}^2$ 和 $F_{{\rm{E}}13}^2$ 。此时 $F_{{\rm{E}}4}^2 \ne F_{\rm{E4}}^1$ ，把 $F_{{\rm{E}}4}^2$ 、 $F_{{\rm{E}}3}^2$ 和 $F_{{\rm{E}}13}^2$ 同时乘以关系系数 $ \;\beta = F_{{\rm{E}}4}^2/F_{{\rm{E}}4}^1 $ ，可获得在 $ F_{{\rm{E}}4}^1 $ 作用下式（17）的解。

4）依次释放其余约束，重复步骤3）可求出所有单元的内力。

5）校核所求的自应力模态是否正确。调整组合系数，可使上、下弦索所对应节点处撑杆的内力相同。若不相同，说明所求的自应力模态不正确，检查关系系数 $\;\beta $ 是否正确，修正后即可。

若不考虑外荷载作用，所有单元内力所构成的集合是式（14）的基础解系，即结构的自应力模态；若考虑外荷载作用，所有单元内力所构成的集合是式（13）的特解；二者叠加为式（13）的通解，即结构的可行预应力。计算流程如图6所示。

4 算　例 4.1 算例1

以文献[10]中的算例为例，如图7所示的平面索桁架，可视为一个轮辐式索桁结构位于同一平面内的两榀索桁架，其中，单元E7和E8为刚性杆，其余为拉索，索杆材料参数见表1。

根据式（11）判断可知，在预应力作用下该结构存在自应力模态。若对单元E1施加5 kN的内力，将E1～E4、E7、E8内力的计算结果和理论值进行对比，说明了本文方法的正确性，结果见表2。

由表2可知，本文改进方法和原方法所求得的结果与理论值基本一致，且前者精度高于后者。本文方法所求得的结果均是迭代2次所得的结果，收敛速度非常快。原方法和本文方法计算结果之间的误差主要是由于原方法采用生死单元技术时并没有删除该单元，而是赋予单元较小的弹性模量。由于该结构的节点很少，因此累计误差较小，但当索桁架的节点加多时，原方法导致的累计误差不容忽视。

4.2 算例2

以文献[10]中的算例为例，如图8所示的轮辐式屋盖共有16榀索桁架。通过任意中心轴的剖面投影可得到2榀对称位置上的索桁架，其节点编号、规格尺寸如图9所示。其中，单元3–4、5–6分别为受压的飞柱和撑杆，其他单元为拉索。单元3–5、5–7为上弦索，4–6、6–7为下弦索，1–3、2–4分别为上、下环索的投影。索截面面积为7.069 mm²，弹性模量为1.8×10⁵ MPa；压杆截面面积为143.728 mm²，弹性模量为2.06×10⁵ MPa。

根据本文方法的计算步骤，首先对结构形状进行判断。由图9可得，索桁架飞柱上、下部分长度比为 ${C_1} $ = 100/200=0.5，撑杆的上、下部分长度比为 ${C_2}$ =70/150=0.467。由于 ${C_1} \ne {C_2}$ ，所以结构形状设计不合理，即不存在独立的自应力模态，与文献[10]结论一致。因此需要对结构位形进行重新设计。修改原则：一般不改变上弦索节点的坐标，以保证屋盖的外形美观和排水方便；下弦索节点一般不能向下移动，以保证视线和采光的要求。

根据上述修改原则，不改变上弦索节点的坐标，使下弦索节点坐标向上移动。假设飞柱和撑杆上、下部分长度比值为 ${C_0} = 0.588$ ，可以得到飞柱下半部分长度为170 mm，撑杆下半部分长度为119 mm（图10）。即节点4向上移动200–170=30 mm，节点6向上移动150–119=31 mm，这与文献[10]中采用逐点去约束法所得结果基本一致。但在文献[10]方法中需要8次反复迭代求解才能得到上述结果，过程繁琐、工作量大。

本文方法求得的解与文献[10]对比结果如表3所示。由表3可知，单元6–7和3–4的误差相对稍大，分别为0.259%和0.189%，其余均在0.1%以内，精度较高。同时，本文方法所获得的值均是在进行2次迭代后所得到的结果，收敛速度很快。

单元6–7和3–4的误差相对稍大的原因是：文献[10]中节点4和6分别向上移动31和30 mm，本文方法对节点4和6分别向上移动30和31 mm，各节点相差1 mm。

4.3 算例3

设计一个跨度为160 m（半径为80 m）的环形交叉索桁架结构，并对其进行找形研究。根据体育场采光要求和座位需求量，取露天部分半径39.44 m，取环形交叉索桁结构的环向等分数为15，单榀索桁架跨越5个环向等分段。根据文献[10]所给出的经济矢跨比范围确定上、下弦索的矢高分别为5.459和7.582 m；然后设定索桁架外轮廓线为某一抛物线，通过数值计算可得出索桁架其余点的坐标值。所设计的环形交叉索桁结构透视图如图11所示，单榀索桁架几何尺寸如图12所示。索的弹性模量为1.80×10⁵ MPa，杆的弹性模量为2.06×10⁵ MPa。采用本文方法计算所得自应力模态结果与理论值对比如表4所示。

由表4可知，单元5–6和7–8的误差相对稍大，分别为0.105%和0.196%，其余均在0.1%以内，精度较高。同时，本文方法所获得的单元内力值均是在进行2次迭代后所得到的结果，如表5所示；如果继续迭代，误差会更小，但变化不大。算例验证了本文方法的稳定性和效率。

5 结　论

通过对索桁张拉结构形状判定和自应力模态求解方法的研究，得到以下结论：

1）基于索桁张拉结构所具有的拓扑关系，并从平面索桁架的合理位形出发推导了形状判定的简洁判断公式。通过该判定条件可快速设计出具有合理位形的索桁架，极大地简化了结构稳定性判断的流程。

2）本文方法基于先形状判断后自应力模态求解的思路，回避了原逐点去约束方法中先求解自应力模态后判断形状的反复求解的缺点；且改进的逐点去约束法未使用原方法中的生死单元技术，保证了求解过程的收敛性，避免了坐标转换，并简化了求解流程；消除了因使用生死单元技术而导致的累计误差，提高了计算精度。

3）本文方法能够在有限元软件上直接应用，并考虑结构的变形和材料属性，更能真实反映结构的内力情况。通过算例表明本文的改进方法具有系统的求解流程，且收敛快，精度高。