椭圆鞍形单层正交索网结构解析计算方法与形态优化

王博厚，孙 逊，方立新

(1.东南大学土木工程学院，江苏，南京 210096；2.东南大学建筑设计研究院有限公司，江苏，南京 210096；3.东南大学建筑学院，江苏，南京 210096)

索网结构因其在充分利用材料强度实现结构功能的同时，能保持整体结构的轻巧简洁，能适应建筑功能及外观需求等优点，受到工程结构界的广泛关注。单层正交索网，由相互正交、曲率相反的两组钢索连接而形成，形成负高斯曲率的曲面。索网悬挂在强大的边缘构件(外环梁)上。单层正交索网可按其平面投影形状分类，如平面投影为菱形、椭圆形、圆形的单层正交索网结构。

1952 年建成的道顿竞技馆被公认为是第一个具有现代意义的大跨度索网屋盖[1]。随后在国内外的大型场馆、公共建筑中得到广泛应用，如伦敦奥运自行车馆[2]、苏州奥体中心游泳馆[3]、国家速滑馆[4]等。对于索网结构的研究，主要包括找形与荷载分析[5]、解析计算[6-7]、静力与动力性能研究[8]、形态优化[9-11]、断索分析[12]、施工模拟等方面。

用于索网结构初始形态确定和荷载状态分析的主要方法有：非线性有限元法[13]、动力松弛法[14]和力密度法[15]。而非线性有限元法是一种更为精确的数值计算方法[16]，且随着有限元计算软件的发展，该方法在找形及荷载分析中应用最广。

沈世钊等[6]、金问鲁[7]将鞍形索网连续化为薄膜，分别利用平衡与变形协调关系、能量原理建立以竖向位移为未知量的非线性方程，其所得解析解精度极大程度依赖于所假定的位移函数。

肖康丽等[8]对索网结构静力性能进行参数分析，研究矢跨比、边界形状、索初始预张力、索直径等参数对变形、内力的影响，其结论可用于指导结构优化设计。目前在索膜结构设计中，大多根据建筑造型要求，在所允许的范围内调整几何参数以确定最终形态，在此过程中虽然会进行一定的方案比选，但该过程大多依赖于设计者的经验，很难保证最终形态在给定约束条件下最优[17]。

对索网结构形态的优化，关键问题在于选取合适的优化目标。形态优化研究根据优化目标可分为单目标优化、多目标优化，目前学者所提出的优化目标可概括为材料用量最少、竖向刚度最大两大类[9-11]，对于材料用料的优化主要针对索和下部支撑结构。已有形态优化方法需借助有限元、优化算法自编程序实现，计算过程繁杂且不易收敛，不能快速判断结构形态是否合理。

索网屋盖结构中，外环梁材料用量大、受力复杂，其截面极大地影响着结构的经济性和安全性，在优化设计中未受到广泛关注。外环梁大多处压弯状态，重力荷载方向的弯矩由支承条件决定，而鞍面的形态对外环梁的面内弯矩有极大影响。外环梁越接近轴压构件，其材料利用越充分，则可相应地减小外环截面尺寸和变形[18]。

针对以上问题，本文首先利用平衡方程、变形协调方程，以竖向位移为未知量推导了正交单层索网连续化计算方法，提出一个新的位移函数，采用伽辽金法求解。其次以外环梁平面内弯矩最小为优化目标，以中央承重索垂度为优化变量，以索力、竖向位移和矢跨比为约束条件。利用外环梁上任一点弯矩计算公式，求解方程得到结构最优形态。最后对一个简单算例分别采用有限元法、解析法计算，验证位移、索力、外环弯矩计算公式和形态优化的准确性。在结构初步设计阶段，利用本方法仅需求解一元方程，便可快速得到结构的最优形态，避免了有限元编程、反复试算的过程。

1 索网连续化计算方法

目前，单层正交索网的计算理论主要有离散化理论和连续化理论[6]。连续化理论通过假设索网在受力变形后的位移函数并以之为未知量，进而求出索网面内的索力分布，避免了离散化理论中大量方程、未知量造成的求解困难。

本文利用连续化理论，通过平衡微分方程、变形协调方程建立以竖向位移为未知量的方程，以变分原理求解，推导正交索网的解析公式。

1.1 基本假定

索网计算时采用如下基本假设：① 索是理想柔性的，只能承受拉力；② 索材料符合胡克定律；③ 两个方向的索之间永远保持接触，能相互传递竖向力的作用；④ 只考虑小垂度问题，且仅考虑竖向均布荷载作用；⑤x方向索为承重索，y方向索为稳定索，相同方向索的初始预拉力、截面面积、间距均相同。

1.2 索网的三种形态

根据索受力状态和外荷载的不同，索网结构分为3 种状态：零状态、初始平衡态、荷载态。假设结构在初始态下的曲面形状为z0(x,y)，索自重折算为竖向均布荷载q0(x,y)，两个方向索折算到单位宽度内的薄膜预张力为Hx0、Hy0，折算到单位宽度内的拉伸刚度为EAx、EAy。

假设对结构施加的外荷载为Δq(x,y)，竖向荷载变为q=q0+Δq。索网在外荷载作用下产生三个方向的位移分量u、v、w，但在小垂度情况下水平位移分量u、v与竖向位移分量w相比是高阶微量，可以忽略其对曲面形状的影响，变形后索网形状变为z=z0+w。索水平张力获得增量ΔHx、ΔHy，达到Hx=Hx0+ΔHx、Hy=Hy0+ΔHy。

索网结构在初始态时的各项参数均为已知，在荷载态时两正交方向的水平拉力Hx(或ΔHx)、Hy(或ΔHy)及位形z(或w)是未知的，求解三个未知量需要三个方程，即x、y方向的变形协调方程和z方向的平衡微分方程。索在受到外荷载作用从初始态转换为荷载态的过程中，索中应变增量很小，但位移较大。这时的平衡条件应建立在变形后的位形上，应变表达式也应包含位移的二次项，即在建立方程时应考虑几何非线性的影响。

1.3 索网平衡方程

在荷载态，从连续化的索网中任取微分单元dxdy，受力情况如图1 所示。

图1 索网微分单元受力简图Fig.1 Stress diagram of cable-net differential element

考虑x、y方向的平衡条件，可得Hx=Hx(y)、Hy=Hy(x)，即Hx在x方向不变，Hy在y方向不变。再考虑z方向平衡条件并将以上结论代入，可得：

式(1)即为荷载态时索网的平衡微分方程。同理可得初始态时索网的平衡微分方程：

若将q=q0+Δq、z=z0+w、Hx=Hx0+ΔHx、Hy=Hy0+ΔHy代入式(1)，并减去式(2)，可得索网平衡方程的增量形式：

1.4 索网变形协调方程

以x方向任意位置的承重索MN为例，假设MN从初始态转换为荷载态时，两端产生的位移为(uM、vM、wM)、(uN、vN、wN)，索温度变化为Δtx。

1) 由几何关系推导索长变化

考察在初始态长度为dsx0的微元AB，在荷载态移动到A′B′，其长度伸长至dsx，如图2 所示。

图2 由几何关系计算索长简图Fig.2 Schematics of calculating cable length from geometric relations

由几何关系可知，微元AB的伸长量为：

在小垂度问题中，略去水平位移u、v的高阶微量，将根号利用级数展开并保留级数前两项，并将z=z0+w代入，可得：

承重索MN的总伸长量为：

式中：M、N为承重索的两个边界点；uN、uM表示这两点在x方向的位移。

2) 由物理关系推导索长变化

由物理关系可知，索的伸长由拉伸变形和温度变化引起，即：

式中：lx为索MN的水平投影长度；ΔHx为索MN索力增量的水平分量；α 为材料线膨胀系数；Δtx为索MN温度变化。在小垂度问题中，认为(dz0/dx)2与1 相比是微量。则式(7)可简化为：

文献[6]表明：当f/l≤0.1 时，这种近似引起的索力误差在5%左右。令式(6)与式(8)相等，可得x方向承重索MN的变形协调方程，即：

同理，y方向稳定索PQ的变形协调方程为：

1.5 索网近似解

从理论上讲，当给定索网初始状态、边界条件、荷载增量时，即已知z0、Hx0、Hy0、q0、Δq，由变形协调方程和平衡方程，可解出ΔHx、ΔHy和w三个未知量。但由于三个未知量是函数式，三个方程是积分或微分方程，直接求解并不容易。

实际应用中可求助于基于变分原理的各种近似解法，建立简便的实用计算公式。常用的近似解法有配点法、力矩法、最小二乘法和伽辽金法等，因伽辽金法求解精度相对较高，采用伽辽金法对结构进行近似求解。具体步骤如下：

1) 根据索网的边界条件和荷载分布情况，假设一个位移试函数w(x,y)，其中包含有限个待定常数，一般形式如下：

式中：fi(x,y) (i=0, 1, …,n)为选定的已知函数，f0在边界处满足给定的边界条件；ci(i=1, 2, …,n)为n个待定参数。

2) 将假定的位移函数w代入变形协调方程，可得索水平张力增量ΔHx、ΔHy的函数，可写为：

3) 将w、z0、ΔHx和ΔHy代入式(3)，得包含各待定参数ci的函数S。

4) 按平衡条件要求，函数S应恒等于0，但由于位移函数是事先假定的，不可能精确。此方法要求函数S与所选位移函数中的各已知函数fi正交，即：

可得n个方程，近似求解出n个待定参数ci。

2 椭圆鞍形正交索网解析解

平面投影为椭圆的鞍形索网由马鞍面和椭圆柱面相切而成。马鞍面与椭圆柱面相交而成的空间曲线即为结构的外环，其在xoy面上的投影即椭圆柱面在xoy面上的投影。索网结构的平面及侧面投影、几何参数、坐标体系如图3 所示，侧视图中实线代表索初始态位置，虚线代表索荷载态位置。初始态下马鞍面、椭圆柱面方程如下：

图3 索网结构平面及侧视图Fig.3 Plan and side view of cable-net structure

式中：f1为中央承重索的跨中垂度；f2为中央稳定索的跨中拱度，设H为结构最高点、最低点之间的高差，则H=f1+f2；a、b分别为椭圆的两个半轴长度。可见，f1、f2、a、b这4 个参数将决定索网结构的形态。

由于索网结构跨度大、自重轻，在初始态忽略索网自重，令q0=0，根据初始态平衡方程，可得：

即两个方向索的初始预张力与曲面在相应方向的曲率成反比。若已知一个方向的预张力，则可确定满足平衡方程的另一方向的预张力。

2.1 位移函数的选取

位移函数w需满足以下条件：① 索网承受对称荷载时，位移函数应同样具有对称的性质；② 索网的边缘构件支承于下部结构之上，本身刚度很大，可认为索网边界上的竖向位移为0。位移函数的选取将极大地影响位移解的精度，此前文献均采用如下形式的位移函数，区别在于所选取的待定参数和已知函数的个数不同。

式中：wmn为需求解的待定参数。对于承受竖向均布荷载的情况，沈世钊等[6]取1 个待定参数w00，其位移函数精度差，不能准确描述索网的变形，但仅需求解一元方程。金问鲁[7]取3 个待定参数w00、w20、w02，其位移函数可准确地描述索网变形，与有限元解十分接近；但需联立求解方程组，很难写出各待定参数的表达式。

因此，本文提出一个新的位移函数，该函数既满足对称性和边界条件，又在仅取一个待定参数的条件下，有较高的精度，以求解索网竖向位移及内力 。位移函数形式如下：

式中：w0为待定参数，表示索网中点的竖向位移。

2.2 单参数三次方程解

采用伽辽金法求解位移函数中的待定参数，对马鞍面方程式(15)、位移函数式(19)代入变形协调方程式(9)与式(10)，沿x、y方向积分，可得索力增量：

将马鞍面方程式(15)、竖向位移函数式(19)及式(20)、式(21)代入平衡方程式(3)，经整理后得：

式中：A0、A1、A2、均为x、y的函数。

由式(13)，伽辽金变分方程为：

将函数S代入伽辽金变分方程，利用椭圆积分运算并整理后，可得用于确定待定参数w0的一元三次方程，即：

一元三次方程可用卡丹求根公式求解，但表达式过于复杂，实际中采用迭代法求解。

2.3 单参数一次方程解

索网在外荷载作用下的中点竖向位移与跨度的比值w0/a、w0/b通常是微量。在不考虑支座位移和温度变化的条件下，为使一元三次方程便于求解与表达，忽略式(24)中w0/a及w0/b的二次项、三次项，将其简化为一次方程：

得待定参数w0的线性表达式：

与三次方程相比，此表达式大大简化，可快速计算出索网中点在荷载态下的竖向位移，并在此基础上计算各节点竖向位移及各索索力增量。索网中任一节点的位移可由式(19)计算，任一根索的索力增量计算公式便简化为：

3 索网形态优化

外环梁大多处压弯状态，重力荷载方向的弯矩由支承条件决定，而鞍面的形态对外环梁的面内弯矩有极大影响，外环梁越接近轴压构件，其材料利用越充分，可相应地减小外环截面尺寸和变形。因此在结构平面投影尺寸a和b、高差H、承重索预拉力Hx0不变的条件下，以外环梁平面内弯矩最小为优化目标，以跨中承重索垂度f1为优化变量，以索力、位移及矢跨比为约束条件，采用解析方法求解结构最优形态。

3.1 外环荷载

为简化计算做出以下假定：① 实际工程中马鞍面高差与跨度相比较小，将空间曲梁简化为平面构件；② 只考虑外环所受平面内弯矩；③ 外环横截面处处相等；④ 由于外环刚度较大，计算时不考虑几何非线性的影响；⑤ 将索对外环的集中力作用转换为分布荷载。

马少华[19]将索对外环的作用转换为双向均布荷载，根据平衡条件推导出外环任一点轴力、弯矩的表达式。但其未考虑索力的真实分布情况，所得弯矩表达式精度较差。以外侧受拉为正，其外环上任一点弯矩M表达式为：

本文考虑索力的真实分布情况，由于对称性，仅取1/4 外环进行分析，外环受力如图4 所示。

图4 1/4 外环受力简图Fig.4 Load condition of 1/4 outer-ring

承重索、稳定索初始索力Hx0、Hy0呈均匀分布，索力增量ΔHx、ΔHy呈曲线分布，索力增量按式(27)、式(28)计算。如何选取形态优化时采用的外荷载Δq，关系到形态优化的优劣，目前主要有两种思路：① 取结构全寿命周期内发生概率最大的荷载组合作为外荷载；② 取各组合中最大、最小荷载增量的平均值作为外荷载[20]。

3.2 外环弯矩及其分布

考虑y方向平衡条件，可得：

式中：NA表示A点轴力。再对B点取距并将NA代入，可得B点弯矩值：

以外侧受拉为正，结合式(29)和式(31)，假设外环上任一点弯矩M表达式为：

式(32)中的弯矩分布符合边界条件，即θ=0°时，M=MA；θ=90°时，M=MB。MA为A点弯矩，可由等截面闭合环的弯矩沿周长积分为零的条件确定，即解得：

将MA代回式(32)，便得到外环上任一点的弯矩M表达式，式中 γ 的含义同式(29)。

1/4 外环弯矩分布曲线如图5 所示，弯矩在A、B两点达到最大，即结构最高点与最低点。

图5 1/4 外环弯矩图Fig.5 Bending moment of 1/4 outer-ring

3.3 索网最优形态求解

若式(34)中括号内结果等于0，外环上任一点弯矩便为0。以f1为变量，固定a、b、H、Hx0保持不变，改变索网空间形态。将式(17)、式(26)代入，可得关于f1的一元四次方程式(35)。求解得到f1后，便得到外环平面内弯矩处处为0 时的索网结构形态。

3.4 约束条件

索网结构设计时还需遵循一定的约束条件，也就是根据安全和设计要求对变量的限制，主要以索力、竖向位移、矢跨比为约束条件。

1) 索力。在任意荷载组合下，拉索不能松弛，索力应始终为正，且最大索力应满足材料强度限值，以索水平力近似代替索力。索力约束条件为：

式中：αc为考虑材料安全储备的安全系数；[σ]为材料容许应力。

2) 竖向位移。结构的竖向位移应满足《索结构技术规程》[21]的要求，即最大挠度不宜大于跨度的1/250。位移约束条件为：

式中，[uz]为节点竖向位移限值。

3) 矢跨比。文献[21]给出了承重索垂度、稳定索拱度与其跨度比值的适宜取值区间，还应根据建筑功能、造型具体考虑其限值。矢跨比约束条件为：

式中，[f/l]max、[f/l]min分别为矢跨比上限、下限值，应分别考虑承重索与稳定索。

综上，用解析法对椭圆单层正交索网结构进行优化设计的完整步骤为：① 给定初始参数，包含椭圆长半轴长度a、短半轴长度b、高差H、承重索初始索力Hx0，确定索力、位移和矢跨比的限值；② 求解式(35)，计算使外环上任一点弯矩为0 时的承重索垂度f1，得到结构最优形态；③ 验算所得形态在各种荷载组合下是否满足约束条件；④ 若满足，则该形态为结构的最优形态；若不满足，可调整初始索力Hx0或索截面面积Ax、Ay，返回步骤② 再次求解。

4 有限元验证

本节将前文所推导的索网位移和索力计算公式、外环弯矩计算公式、形态优化结果，与有限元分析结果进行对比，以验证解析计算方法的准确性。为便于比较，取与文献[7]P212 算例相同的几何参数、荷载条件，具体参数如表1 所示。

表1 椭圆单层正交鞍形索平网结构参数Table 1 Parameters of elliptic saddle single-layer orthogonal cable-net structure

4.1 有限元找形、荷载分析方法

采用ANSYS APDL 参数化设计语言编写找形及荷载分析程序。由于索中预拉力远大于自重引起的张力，索可用只受拉的杆单元模拟，采用Link10 单元模拟悬索，初始预拉力通过初应变施加；采用Surf154 单元模拟膜面施加均布荷载；采用Beam189 单元模拟外环梁。程序分为以下3 部分：

1) 找形分析

对索网结构进行荷载分析之前，首先应通过找形分析，建立初始态时索网结构模型。找形分析将初始预张力分布、边界条件作为已知，寻求对应的平衡形态。本文采用一种改进的综合找形法[22]进行初始形态分析。该方法在支座位移法找形后，约束边界点，采用节点平衡法迭代求解，此时形态变化较大，控制位移达到要求精度；再恢复真实弹模、初应变，再度采用节点平衡法迭代求解，此时索力变化较大，控制索力达到要求精度。共3 次找形过程，具体分析流程见图6(a)。

图6 有限元分析流程图Fig.6 Finite element analysis process

对案例进行分析，找形后各节点z坐标值与理论双曲抛物面最大相差3.5 mm，索力值与预设初始预拉力最大相差0.47 kN。该方法找形结果准确，且找形过程共经历7 次计算，收敛快、计算效率高。

2) 索网荷载分析

在找形后的模型上，借助表面效应单元施加竖向均布荷载，并将其转换为节点荷载施加于索网自由节点上，具体分析流程见图6(b)。

3) 外环荷载分析

在荷载态模型基础上，记录边界节点上的节点力，删除索网并建立外环，将所记录的节点力施加给外环梁，具体流程见图6(c)。

4.2 索网位移、索力对比

将有限元所得位移、索力结果与上文及文献[6 - 7]的解析计算结果进行对比，近似求解方法均采用伽辽金法。在所有节点中，中点附近位移最大，其值决定了结构是否能满足挠度限值，因此仅对比中点竖向位移结果，如表2 所示。在所有拉索中，中央承重索、稳定索的索力变化量最大，其索力值决定了所有索是否松弛与是否满足强度要求，因此仅对比中央承重索、稳定索索力增量，如表3 所示。

表2 不同求解方式中点竖向位移解对比Table 2 Comparison of vertical displacement of midpoints in different solutions

表3 不同求解方式中央索索力增量对比Table 3 Comparison of central cable force increment in different solutions

应用本文位移函数的单参数一次方程解所得位移、索力结果，与本文单参数三次解相比，结果相差很小，说明将三次方程简化为一次方程是可行的。与有限元相比，其误差在可接受范围内；与文献[6]单参数解相比，位移解结果精度大大提高；与文献[7]三参数解相比，虽然精度稍有差距，但所需运算量大大减小且表达式简便，便于后续求解结构最优形态。因此可使用本文所提出的计算方法快速估算索网结构在各种荷载作用下的位移、内力。

4.3 外环弯矩对比

1) 不同位置弯矩值对比

取表1 算例，对1/4 外环计算结果进行对比，每15°取一测点，分别采用有限元法、式(29)和式(34)计算，外环弯矩值随角度θ 的变化如图7所示。相比于式(29)的计算方法，本文的弯矩计算公式所得结果与有限元解更接近，在结构最高点(θ=0°)、最低点(θ=90°)误差分别为1.1%、3.4%。

图7 不同角度θ 时外环弯矩值对比Fig.7 Comparison of outer-ring moment at different angles θ

2) 不同形态弯矩值对比

固定H=5 m，使f1在2.8 m～4.0 m 区间内变化，f2、Hy0随之改变，其余参数保持不变，分别采用有限元法、式(34)计算外环A、B点弯矩值。弯矩值随f1的变化如图8 所示，解析解与有限元结果基本吻合。A点弯矩随着f1的增加，从内侧受拉转变为外侧受拉，B点则相反；A、B点弯矩均接近0 时，外环越接近轴压状态。

图8 不同垂度f1 时外环弯矩值对比Fig.8 Comparison of outer-ring moment at different sags f1

综上，通过与有限元结果对比，验证了本文所提出的外环弯矩计算公式的准确性。

4.4 最优形态对比

采用表1 案例，固定H=5 m，改变f1、f2以求解结构最优形态。通过式(35)解得f1=3.52 m，则f2=1.48 m，即为解析方法求得的索网最优形态。

实际情况中并不存在一个使外环梁弯矩处处为0 的形态，在有限元模型中以外环弯矩绝对值的最大值最小为目标，以f1为变量，寻求最优形态。求得f1=3.53 m 时，为有限元法得到的结构最优形态。

分别取表1 案例、解析法与有限元法所得最优形态的几何参数，采用有限元法进行荷载分析，外环梁上最大弯矩、轴力值对比如表4 所示。

表4 不同形态下最大弯矩、轴力值对比Table 4 Comparison of maximum bending moment and axial force under different shapes

与表1 案例相比，解析法与有限元法所得最优形态下，外环梁上最大弯矩均大幅度减小并接近于0，最大轴力小幅度增加。外环中应力大大减小，可减小外环截面面积，使得材料强度充分利用。解析法与有限元法形态优化结果仅相差0.01 m，且外环上弯矩、轴力最大值相差很小，解析法所得形态优化结果具有相当高的精度。

5 结论

本文推导了椭圆鞍形正交单层索网连续化计算方法和外环梁上任一点弯矩计算公式。随后，提出一种快速求解此类结构最优形态的方法，该方法以外环面内弯矩最小为优化目标，以中央承重索垂度为优化变量，以索力、竖向位移和矢跨比为约束条件，利用外环弯矩计算公式，得到外环弯矩最小时的结构形态。本文得到的主要结论如下：

(1) 对椭圆鞍形单层正交索网结构的连续化计算方法中所假设的位移函数进行了改进，利用变形协调方程、平衡微分方程结合伽辽金法推导出结构荷载态时的位移、内力。改进后的位移函数既具有较高的精度又使得表达式简便，并采用有限元验证了计算公式的准确性，可利用本文公式快速计算出给定结构的位移、内力。

(2) 利用所得索力，对外环施加更符合真实情况的荷载作用，改进了外环上任一点弯矩的计算公式。通过与有限元对比分析，本文外环弯矩计算公式具有较高的精度。

(3) 提出一种快速求解索网结构最优形态的方法，该方法以外环梁面内弯矩最小为优化目标，以中央承重索垂度为优化变量，以索力、竖向位移和矢跨比为约束条件。利用外环弯矩计算公式，得到以优化变量为未知数的一元方程，求解方程便快速得到结构的最优形态。解析方法所得结构最优形态与有限元结果基本一致，优化结果具有相当高的精度。在结构初步设计阶段，利用本文解析计算方法，可快速、准确的获得结构最优形态，避免了有限元法编程、反复试算的繁杂过程。