一种索杆张力结构的找形方法

摘要

本发明公开了一种索杆张力结构的找形方法，属于空间结构的建筑设计和结构设计领域。本发明方法步骤是：设定索杆张力结构的拓扑关系，并给出其约束条件和压杆长度；确定索杆张力结构找形分析中的目标函数；设定每根拉索和压杆的长度权重系数以及每个边界条件的权重系数，并给出梯度法优化中的初始步长；计算每个节点坐标变量在xk处的梯度向量。本发明从能量的角度出发，并将索杆张力结构的边界条件等约束写入梯度法优化方法的功能函数中，从而较容易地解决了具有约束条件的张力结构的找形问题。传统找形方法中加入了过多的约束条件以及几何限制，本发明初始限制条件较少，可以找到大量非规则体系。

说明书

技术领域

本发明公开了一种基于梯度法优化的索杆张力结构的找形方法，属于空间结构的建筑 设计和结构设计领域。

背景技术

索杆张力结构的刚度是由预应力提供的，它们需要一个过程达到自平衡状态，这个状 态一般称为预应力态。而结构的预应力态取决于结构的形状，所以达到预应力态的过程也 称之为找形。而索杆张力结构的找形同时也是一个找力的过程，在具体的找形方法中可以 以形状参数为变量，也可以以内力参数为变量。因此，可以将找形方法划为“找力”和“找 形”两大类，前者主要指搜索可行预应力或预应力优化，后者的代表性方法有力密度法、 非线性有限元法和动力松弛法。力密度法最早由Schek提出，用于索网等只有受拉构件结 构的找形分析中。随后，众多学者对力密度法进行了改进，使之能够运用于含有拉压杆的 张拉整体结构的找形。但在运用力密度法对索杆张力结构找形时，确定力密度后，结构的 几何构形并不唯一，而这一问题一直没有很好的解决。

在一个三维索杆结构中，假定自由节点i与节点j、节点k相连，如图1所示。索杆体 系的节点平衡方程可以表述为：

$\left(\begin{array}{c}(x_i-x_j)f_{\mathrm{ij}}/l_{\mathrm{ij}}+(x_i-x_k)f_{\mathrm{ik}}/l_{\mathrm{ik}}=p_{\mathrm{ix}}\\ (y_i-y_j)f_{\mathrm{ij}}/l_{\mathrm{ij}}+(y_i-y_k)f_{\mathrm{ik}}/l_{\mathrm{ik}}=p_{\mathrm{iy}}\\ (z_i-z_j)f_{\mathrm{ij}}/l_{\mathrm{ij}}+(z_i-z_k)f_{\mathrm{ik}}/l_{\mathrm{ik}}=p_{\mathrm{iz}}\end{array}\right)$式1

对于任意一个连接节点i、j的单元(i，j)，其内力为f_ij，长度为l_ij；作用于节点i的外 荷载为p_i，在三个方向的分量分别为：p_ix，p_iy，p_iz。力密度为单元内力和长度的比值，即 为q_ij＝f_ij/l_ij。式1可整理为：

$\left(\begin{array}{c}(q_{\mathrm{ij}}+q_{\mathrm{ik}})x_i-q_{\mathrm{ij}}{x}_j-q_{\mathrm{ik}}{x}_k=p_{\mathrm{ix}}\\ (q_{\mathrm{ij}}+q_{\mathrm{ik}})y_i-q_{\mathrm{ij}}{y}_j-q_{\mathrm{ik}}{y}_k=p_{\mathrm{iy}}\\ (q_{\mathrm{ij}}+q_{\mathrm{ik}})z_i-q_{\mathrm{ij}}{z}_j-q_{\mathrm{ik}}{z}_k=p_{\mathrm{iz}}\end{array}\right)$式2

假定索杆张力结构有b个单元，n_s个节点，其中n个自由节点，n_f个固定节点，引入 b×n的拓扑矩阵：

式3

在C_s中将自由节点排列在固定节点之前，C_s又可以拆分为自由节点拓扑矩阵和固定 节点拓扑矩阵：C_s＝[C C_f]，则自由节点的力平衡关系可表示为(以x方向为例)：

C^TQCx+C^TQC_fx_f＝P_x    式4

式中：x为自由节点坐标向量；x_f为固定节点坐标向量；P_x为节点外荷载x方向向量； C为自由节点拓扑矩阵；C_f为固定节点拓扑矩阵；Q为力密度对角矩阵。

在初始预应力确定的过程中，不考虑外荷载的作用，并且如果仅考虑结构是自平衡的 (没有固定点)，则结构的平衡可以简化为：

$\left(\begin{array}{c}\mathrm{Dx}=0\\ \mathrm{Dy}=0\\ \mathrm{Dz}=0\end{array}\right)$式5

其中D为力密度矩阵，可以表示为：

D＝C^TQC    式6

图2所示为平面张拉整体结构，由4根拉索和2根压杆组成，压杆通过拉索相连，其 本身互不相连。对于这种自平衡结构，需满足平衡条件式5。如果式5中力密度矩阵D是 满秩矩阵，则节点坐标需满足：

x＝y＝z＝0    式7

这样图2所示的结构就缩为1个点了。如果力密度矩阵D是非满秩矩阵，则可以通过奇异 值分解等数学手段得到其零向量空间。假定4根拉索和2根压杆的力密度比值为： 1∶1∶1∶1∶-1∶-1，则可以获得多组节点坐标解。如图3和图4所示的两种几何构形都满足式5。 该结构的力密度矩阵D和节点坐标可以表示为：

$D=\left(\begin{array}{cccc}-1& -1& 1& 1\\ -1& -1& 1& 1\\ 1& 1& -1& -1\\ 1& 1& -1& -1\end{array}\right),$$x=a\left(\begin{array}{c}1\\ 1\\ 1\\ 1\end{array}\right)+b\left(\begin{array}{c}1\\ -1\\ 0\\ 0\end{array}\right)+c\left(\begin{array}{c}0\\ 0\\ 1\\ -1\end{array}\right)$式8

式中a，b，c为系数，可取任意值。

如果4根拉索和2根压杆的力密度比值不是：1∶1∶1∶1∶-1∶-1，则节点坐标无解。结构的 力密度矩阵D和节点坐标可以表示为：

$D=\left(\begin{array}{cccc}-3& -1& 2& 2\\ -1& -3& 2& 2\\ 2& 2& -1& -3\\ 2& 2& -3& -1\end{array}\right),$$x=a\left(\begin{array}{c}1\\ 1\\ 1\\ 1\end{array}\right),$$y=b\left(\begin{array}{c}1\\ 1\\ 1\\ 1\end{array}\right),$$z=c\left(\begin{array}{c}1\\ 1\\ 1\\ 1\end{array}\right)$式9

结构的所有节点都集中于一点上。

发明内容

为了克服力密度法找形过程中几何构形不唯一的缺点，本发明从能量的角度出发，提 供了一种基于优化方法的索杆结构找形方法。另外，现有的找形方法大多针对自平衡的索 杆体系，而土木工程中应用的结构大多有约束条件(支座等边界条件)。本发明将这些约 束条件写入梯度法优化方法的功能函数中，从而较容易地解决了该问题。

一种索杆张力结构的找形方法，其特征在于，包括以下步骤：

第一步、设定索杆张力结构的拓扑关系，并给出其约束条件和压杆长度；

第二步、确定索杆张力结构找形分析中的目标函数，

$f\left(x\right)=\underset{j=1}{\overset{b-n_s}{\Sigma }}{w}_j{l_j}^4\left(x\right)+\underset{m=1}{\overset{n_s}{\Sigma }}w\_{\mathrm{str}}_m{(L_m\left(x\right)-l_m)}^4+\underset{i=1}{\overset{n_f}{\Sigma }}{w\_\mathrm{cor}}_i{(X_i-x_i)}^2$

其中，x为所有节点坐标向量，表示为：{x₁，x₂，x₃，…，x_p，…，x_d*n}^T；n为结构的节点数； d为结构的维数，当体系为二维结构时，d＝2；当体系为三位结构时，d＝3；l_j(x)为第j根拉 索的计算长度，是坐标向量x的函数；L_m为第m根压杆的计算长度，是坐标向量x的函数； l_m为第m根压杆的给定长度；X_i是步骤1给出第i个约束条件的节点坐标值；b为索杆张 力结构的杆件数；n_s是索杆张力结构的压杆数量；n_f是约束条件的数量；w_j为第j根拉索 的长度权重系数；w_str_m为第m根压杆的长度权重系数；w_cor_i为第i个边界条件的权重 系数；

第三步、设定每根拉索和压杆的长度权重系数以及每个边界条件的权重系数，并给出 梯度法优化中的初始步长step；

第四步、输入所有节点坐标向量的初始值x⁰，并带入目标函数，计算得到目标函数为 f⁰；

第五步、在第k步循环过程中，其对应的节点坐标向量为x^k，其目标函数值为f^k，计 算每个节点坐标变量在x^k处的梯度向量为：

$▿f={[\frac{\partial f\left(x^k\right)}{\partial x_1},\frac{\partial f\left(x^k\right)}{\partial x_2},...,\frac{\partial f\left(x^k\right)}{\partial x_{d*n}}]}^T$

判断梯度向量中第p个元素是否大于零，如果大于零，x_p^k+1＝x_p^k-step；否则，x_p^k+1＝ x_p^k+step。从而可以得到下一个循环步的节点坐标向量x^k+1；

第六步、计算x^k+1对应的目标函数值f^k+1，如果f^k+1＞f^k，则step＝step/4；否则，step不 变，重复第五步，直至完成所有循环步，输出最终节点坐标向量x。

本发明的有益效果是，本发明提出的索杆体系“找形”的方法适合于绝大多数索杆体 系，无论是简单的平面结构，还是复杂的空间结构。在给定的初始条件下，该发明可以找 到同时满足能量最小和一定形态要求的索杆张力体系。该发明的另外一个优点是可以在目 标函数中直接引入边界约束条件，从而可以很方便地求解考虑约束条件的索杆张拉体系的 几何形状。传统找形方法中加入了过多的约束条件以及几何限制，因此适合规则索杆体系 的找形，本发明初始限制条件较少，可以找到大量非规则体系。

附图说明

下面结合附图及实施例对本发明进一步详细说明：

图1索杆张力结构节点平衡示意图。

图2两杆四索平面张力结构示意图。

图3两杆四索平面张力结构构形一。

图4两杆四索平面张力结构构形二。

图5是本发明计算流程图。

图6是平面索杆张力结构示意图。

图7是平面索杆张力结构计算构形。

图8是四压杆索穹顶结构示意图。

图9是四压杆索穹顶结构计算构形。

图10是五压杆索穹顶结构示意图。

图11是五压杆索穹顶结构计算构形。

具体实施方式

下面结合具体实施案例对本发明进行更为详细的表述：

组成索杆张力结构的拉索和压杆，都只受到轴力的作用，其应变能V可以表述为：

$V=\frac{1}{2}\frac{F^{2}l}{\mathrm{EA}}---\left(1\right)$

其中F为杆件轴力，l为杆件长度，EA为杆件的轴向刚度。式(1)还可以表述为：

$V=\frac{1}{2}q\frac{F}{\mathrm{EA}}{l}^2={\mathrm{wl}}^2---\left(2\right)$

式中q为杆件的力密度，w为权重系数。

和力密度法相等价的能量函数П(x)可以表示为：

$\Pi \left(x\right)=\underset{j}{\Sigma }{q}_j{l_j}^2\left(x\right)---\left(3\right)$

其中，q_j为第j根杆件的力密度，l_j(x)为第j根杆件的计算长度，向量x为节点的x，y，z向 的坐标，可以表述为：{x₁，x₂，x₃，…，x_p，…，x_d*n}^T。

考虑式(3)的驻值条件为：

$\frac{\partial \Pi \left(x\right)}{\partial x}=▿\Pi =\underset{j=1}{\overset{b}{\Sigma }}{2q}_j{l}_j▿l_j=0---\left(4\right)$

式(4)和节点的平衡方程等价，所以式(3)为力密度法所对应的能量函数。

假定压杆的长度为定值，则可以利用Lagrange乘子法将能量函数П(x，λ)修改为：

$\Pi (x,\lambda )=\underset{j}{\Sigma }{q}_j{l_j}^2\left(x\right)+\underset{m}{\Sigma }{\lambda }_m(L_m\left(x\right)-l_m)---\left(5\right)$

其中，λ_m为第m根压杆所对应的Lagrange乘子；L_m(x)为第m根压杆的计算长度，是坐标 向量x的函数；l_m为第m根压杆的给定长度。对于图2所示的索杆结构，如果我们将所有 拉索的力密度比值定为：1∶1∶1∶1，两根压杆的长度也相等，则图3和图4所示的两种几何 构形都符合式(5)所对应的驻值条件。如果将式(5)的功能函数П(x，λ)假定为其它形式：

$\Pi (x,\lambda )=\underset{j}{\Sigma }{w}_j{l_j}^4\left(x\right)+\underset{m}{\Sigma }{w\_\mathrm{str}}_m{(L_m\left(x\right)-l_m)}^4---\left(6\right)$

其中w_j为第j根拉索的长度权重系数；w_str_m为第m根压杆的长度权重系数。式(6)给 出的函数没有实际的物理意义。这个时候将图2所示结构的拉索的力密度比值定为： 1∶1∶1∶1，两根压杆的长度也相等，根据式(3)的驻值条件，可以得到图3所示的几何构形 是满足条件的唯一解。这主要是因为，式(3)给出的功能函数和长度的平方相关，而根 据勾股定理可知，在直角三角形中，斜边长度的平方永远等于两直角边长度的平方和。所 以满足式(3)驻值条件的几何构形有无穷多组。而式(6)中没有考虑索杆张力结构的边 界条件，确定考虑边界条件约束的索杆张力结构初始平衡态包含如下步骤：

(1)设定索杆张力结构的拓扑关系，并给出其约束条件和压杆长度；

(2)确定索杆张力结构找形分析中的目标函数，

$f\left(x\right)=\underset{j=1}{\overset{b-n_s}{\Sigma }}{w}_j{l_j}^4\left(x\right)+\underset{m=1}{\overset{n_s}{\Sigma }}{w\_\mathrm{str}}_m{(L_m\left(x\right)-l_m)}^4+\underset{i=1}{\overset{n_f}{\Sigma }}{w\_\mathrm{cor}}_i{(X_i-x_i)}^2---\left(7\right)$

其中，x为所有节点坐标向量，表示为：{x₁，x₂，x₃，…，x_p，…，x_d*n}^T，n为结构的节点数， d为结构的维数(当体系为二维结构时，d＝2；当体系为三位结构时，d＝3)；l_j(x)为第j 根拉索的计算长度，是坐标向量x的函数；L_m(x)为第m根压杆的计算长度，是坐标向量x 的函数；l_m为第m根压杆的给定长度；X_i是步骤1给出第i个约束条件的节点坐标值；b 为索杆张力结构的杆件数；n_s是索杆张力结构的压杆数量；n_f是约束条件的数量；w_j为第 j根拉索的长度权重系数；w_str_m为第m根压杆的长度权重系数；w_cor_i为第i个边界条件 的权重系数。

(3)设定每根拉索和压杆的长度权重系数以及每个边界条件的权重系数，并给出梯 度法优化中的初始步长step。

(4)输入所有节点坐标向量的初始值x⁰，并带入目标函数，计算得到目标函数为f⁰。

(5)在第k步循环过程中，其对应的节点坐标向量为x^k，其目标函数值为f^k。计算每 个节点坐标变量在x^k处的梯度向量为：

$▿f={[\frac{\partial f\left(x^k\right)}{\partial x_1},\frac{\partial f\left(x^k\right)}{\partial x_2},...,\frac{\partial f\left(x^k\right)}{\partial x_{d*n}}]}^T---\left(8\right)$

判断梯度向量中第p个元素是否大于零，如果大于零，x_p^k+1＝x_p^k-step；否则，x_p^k+1＝ x_p^k+step。从而可以得到下一个循环步的节点坐标向量x^k+1。

(6)计算x^k+1对应的目标函数值f^k+1。如果f^k+1＞f^k，则step＝step/4；否则，step不变， 重复第5步，直至完成所有循环步，输出最终节点坐标向量x。

利用图5所示的流程编制相应程序，结构具体算例对本发明的具体实施方式作出更为 详细的说明：

算例1图6所示的二维索杆结构是由2根压杆和6根拉索组成。其已知条件为：压 杆(杆件编号分别为7，8)的长度均为4；已知节点5的坐标为(0，0)，节点6的坐标 为(0，7)；压杆的长度权重系数为1000，边界条件的权重系数为5000。将6根拉索的权重 系数矩阵w设定为[1 1 1 1 1 1]^T，则6个节点的坐标为：x＝[2.1692 4.8308 4.8308 2.1692 0.0000 7.0000]^T；y＝[-1.9995 -1.9995 1.9993 1.9993 -0.0000 -0.0000]^T。其几何形状如图7所示。

算例2四压杆索穹顶结构是由4根压杆和16根拉索组成，如图8所示。其已知条件 为：压杆的长度矩阵均为1；固定支座节点1～4坐标分别为(0，0，0)，(10，0，0)，(10，10，0) (0，10，0)；压杆的长度权重系数为1000，边界条件的权重系数为5000。

若将16根拉索的权重系数矩阵w设定为[1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2]^T，则8个自由 节点的坐标见表1。其几何形状如图9所示。

表1四压杆索穹顶节点坐标

算例3五压杆索穹顶结构是由5根压杆和20根拉索组成，如图10所示。其已知条件 为：压杆的长度矩阵均为2；固定支座节点1～5坐标分别为(-10.0000，0，0)，(-3.0902， -9.5106，0)，(8.0902，-5.8779，0)，(8.0902，5.8779，0)，(-3.0902，9.5106，0)；压杆 的长度权重系数为1000，边界条件的权重系数为5000。

若将16根拉索的权重系数矩阵w设定为[1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2]^T，则10个 自由节点的坐标如表2所示。其几何形状如图11所示。

表2五压杆索穹顶节点坐标