一种基于联合迭代算法树状结构形优化设计方法

本发明涉及建筑结构工程技术领域，特别涉及一种基于联合迭代算法树状结构形优化设计方法。

树状结构由于美观的造型以及高效的荷载传递形式而广泛的应用与机场和火车站等大型公共建筑中，如图1所示。关于树状结构的科学问题中，形分析一直是广大学者的关注点。逆吊法是树状结构的一种很典型的分析方法。在数值分析技术快速发展之前，上世纪末和本世纪初，树状结构的形主要是通过试验进行。Kolodziejczyk通过将丝线模型浸在水中，利用水的表面张力作用实现树状结构形；Buelow提出干丝线模型方法对树状结构进行形，结果表明利用丝线模型对简单的树状结构的形具有较高的精度。

伴随着数值分析技术的发展，数值计算分析方法也被广泛的用于树状结构的形分析。Hunt等提出利用虚拟支座对树状结构进行形的分析方法。国内学者武岳在一系列深入研究的基础上，提出了逆吊递推形法。张倩等通过只受拉连续折现索单元解决了树状结构的形问题。作者在将双单元与数值逆吊法进行结合，提出了基于构件欧拉临界力的长度优化算法。基于数值逆吊法的树状结构形方法，利用双单元模拟丝线模型，进而达到了对树状结构进行形的目的。但是双单元的抗弯刚度需要靠经验设置，抗弯刚度太大会引起收敛速度慢，太小会导致计算结果不收敛的问题。在以往研究的基础上，本发明将平衡矩阵法与欧拉稳定承载力计算理论结合，在利用平衡矩阵法对树状结构进行形的同时，利用欧拉稳定承载力理论对分枝内力与几何长度的比值进行修正，达到了对树状结构稳定承载力进行优化的目标。

平衡矩阵法已被广泛的应用于索网等张拉结构的形中。利用平衡矩阵法，根据每个构件的平衡矩阵以及边界节点的坐标可以求得在满足给定平衡矩阵条件下自由节点的坐标。平衡矩阵法主要针对受拉杆系结构的形提出，因此未考虑受压构件的稳定承载力问题。

为了解决大型复杂树状结构形以及传统方法未考虑受压构件的稳定承载力的问题，本发明采用如下技术方案：

一种基于联合迭代算法树状结构形优化设计方法，首先计算每一个分枝内力与欧拉稳定承载力的比值，得到矢量R：

步骤1：根据树状结构初始树状结构几何构型，使各节点受力平衡为前提，根据拓扑矩阵，荷载，固定点坐标，内力与分枝几何长度的比值对树状结构进行初次形；

步骤2：计算得到分枝的临界荷载向量：

步骤3：根据树状结构分枝的实际内力矢量F以及稳定承载力矢量Πcr可得矢量R：

R＝F./Γcr； (13)

式中F为所有分枝的内力构成的矢量F＝[F1，F2，···，Fn]；符号./为表示列向量各个对应元素分别相除；

步骤4：计算R中所有元素的均值ave(R)，当R向量中的第i个元素Ri小于ave(R)时，说明该分枝内力与其失稳临界荷载的比值偏小；

进一步，计算得到每一个构件的内力与长度比值的改变量Δq；将计算得到的Δq直接与原来的q叠加，如式(15)，当Δqi为负值时，减小qi，当Δqi为正值时，增大qi；然后根据更新后的q再次利用平衡矩阵法形；

Δq＝(R-ave(R)).*Γcr/l/ψ (14)

q＝q+Δq (15)

式中ψ为调节收敛速度的变量，为保证计算稳定，ψ为一个大于1的数，本发明中案例取值为20。

进一步，通过利用更新的矢量q再次进行形计算，直到Δq是10-3为止，实现分枝长度的优化，最终的优化目标是所有构件同时达到稳定极限承载力，使树状结构的稳定承载力最大化。

图1为本发明具体实施中树状结构分级示意；

图2为本发明具体实施中树状结构分枝类型示意；

图3为本发明具体实施中构件分枝内力与几何长度的比值调整示意图；

图4为本发明具体实施中树状结构形-优化流程图；

图5为本发明具体实施中树状结构初始形状及构件编号；

图6为本发明具体实施中形结果对比；

其中，(1)SI和μI；(2)SI和μIII；

图7为本发明具体实施中各分枝Δq数值大小；

图8为本发明具体实施中不同条件下形结果；

其中，(1)SI和μII；(2)SI和μIV；

图9为本发明具体实施中不同荷载幅值下形结果；

其中，(1)I；(2)II；

图10为本发明具体实施中收敛曲线；

图11为本发明具体实施中优化后树状结构内力图；

其中，(1)弯矩图(单位:N×m)；(2)轴力图(单位:N)；

图12为本发明具体实施中空间树状结构初始形态；

其中，(1)侧视图；(2)透视图；

图13为本发明具体实施中空间树状结构初始形态；

其中，(1)μI；(2)μII；(3)μIII；(4)μIV；

下面结合附图对本发明具体实施方式作进一步的详细说明。

为了解决大型复杂树状结构形以及传统平衡矩阵法中未考虑受压构件的稳定承载力的问题，如图1至图13所示，本发明提供了一种基于联合迭代算法树状结构形优化设计方法，其特征在于，包括以下步骤：

步骤1：根据树状结构建立几何模型，所述树状结构的拓扑是指不同节点之间通过分枝相互连接的网络。本发明对树状结构进行分析之前首先对分枝的级别进行定义，如图1所示。最底部的构件为树干，然后依次向上为一级分枝，二级分枝，······。

根据上一级每个分枝所分出的下一级分枝的数量，可以将树状结构定义为不同的类型，如所示分别为2分枝、3分枝、4分枝和5分枝树状结构，如图2所示。为简化分析，本发明假定同一个树状结构所包含的分枝类型相同。

步骤2：基于平衡矩阵法对树状结构形，所述平衡矩阵法利用索杆结构体系的拓扑矩阵，利用每一个构件的内力与其长度为比值q＝F/l进行求解。进而得到所有自由节点的空间坐标。

树状结构形的目标是合理安排节点的位置，以使各个分枝只受到轴力的作用。从力学角度讲，树状结构的形与基于平衡矩阵法的索网结构形完全一致。因此本发明利用平衡矩阵法对树状结构进行形分析。根据树状结构各个节点之间的连接关系可以得到树状结构的拓扑矩阵C＝[Cf Cr]。C为一m×n阶矩阵，其行对应于相应的树状结构分枝信息，其列对于相应的节点信息，其中Cf中仅包含自由节点信息，Cr中仅包含约束节点信息。对于树状结构，将承受荷载的节点以及树干底部节点作为约束节点。拓扑矩阵中各个元素的取值根据式(1)进行确定。

对于树状结构，假设不承受荷载的节点个数为nf，约束节点的个数为nr，第i个分枝所承受的压力为Fi。则形后，树状结构自由节点的坐标满足式(2)。

式中，为与第j号节点相连的分枝的个数，(xi,1 yi,1 zi,1)、(xi,2 yi,2 zi,2)分别为与j节点相连的第i个分枝的两个节点的坐标值。将式(2)写成矩阵的形式如式(3)所示，其中三个子方程分别代表树状结构x，y和z三个方向的平衡方程。

式中U，V和W如(4)所示，其中u、v、w为各个分枝两个节点间的坐标差矢量，可通过式(5)得到；q为一m维矢量，其元素为各个分枝内力与几何长度的比值。由q可得到m阶对角方阵Q＝diag(q)。

由于U，V，W以及Q均为对角矩阵，则式(6)成立。根据式(5)，(6)可将式(3)中x向的平衡方程改写为如式(7)所示。将式(7)展开可得式(8)。

Uq＝Qu；Vq＝Qv；Wq＝Qw (6)

式中，下标’f’表示是与自由节点相关，下标’r’表示与约束节点相关。

假定树状结构所有自由节点的外力矢量fxf＝0，且定义则式(8)-a，变为式(9)。

Dfxf＝-Drxr (9)

同理可得到所有自由节点与约束节点的坐标关系式，如式(10)所示。

则根据约束节点的坐标即可确定自由节点的坐标，完成树状机构的形。但是从式(8)可知，自由节点坐标的求解需要提前知道各个分枝构件的内力与几何长度的比值。而在实际工程中，树状结构各个分枝的内力与几何长度的比值相差较大，越接近树干内力越大，很难提前给出一组合理的内力与几何长度的比值矢量q。因此接下来就是解决如何寻一组合理的内力与几何长度的比值矢量。

步骤3：基于平衡矩阵法的稳定承载力优化

步骤3.1：根据步骤2形后，在竖向荷载的作用下，树状结构的各个分枝只受到轴向压力的作用。在轴向压力作用下，构件所能承受的稳定荷载通过式(11)确定。由于树状结构中各个分枝在荷载作用下都受到压力的作用，而上节所述平衡矩阵法是针对受拉的索网结构提出，未考虑构件在受压作用下的稳定问题，因此对于树状结构不能直接采用。

式中μi、li、Ei、Ii分别为构件计算长度系数、几何长度、材料弹性模量以及截面的轴惯性矩。

在利用平衡矩阵法对树状结构进行形时，最理想的目标是直接利用每个构件的欧拉稳定承载力除以其几何长度作为其力密度。形后树状结构的各个分枝在竖向荷载作用下会同时失稳，树状结构的稳定承载力会最大化。但是从式(11)可知，欧拉稳定承载力的确定需要分枝的几何长度，而形之前无法确定各个分枝的几何长度。且当分枝的截面惯性矩提前指定的前提下，在给定荷载下构件未必会达到失稳临界荷载。

因此，设构件在给定荷载作用下的实际内力与稳定承载力的比值为式中F为在给定荷载作用下，构件的实际轴力。在荷载作用下，各个构件的内力随荷载的增加成比例增加。则当各个构件的比值R相同时，各个构件会同时达到稳定临界荷载，以此保证树状结构整体稳定承载力最大。

由于提前并不知道各个构件的几何长度，因此需要利用迭代算法进行求解。首先，假定树状结构的初始形状，根据在给定荷载作用下分枝的内力以及初始形状的几何长度利用平衡矩阵法对树状结构进行初次形。在此基础上计算得到分枝的临界荷载向量，如式(12)所示。

式中μi、li、Ii分别为第i个分枝的计算长度系数、几何长度、截面的轴惯性矩；n为树状结构分枝的数量。

则根据树状结构分枝的实际内力矢量F以及稳定承载力矢量Πcr可得矢量R。

R＝F./Γcr (13)

式中F为所有分枝的内力构成的矢量F＝[F1，F2，···，Fn]；符号./为表示列向量各个对应元素分别相除。

当矢量R中各元素相等时树状结构的稳定承载力达到最大化。为了使矢量R中各元素趋于相等，计算R中所有元素的均值ave(R)。当R向量中的第i个元素Ri小于ave(R)时，说明该分枝内力与其失稳临界荷载的比值偏小。此时应该减小该分枝的内力与几何长度的比值，以优化分枝的内力和长度。如图3所示一级分枝平面树状结构。树干与分枝的长度和内力相等。此时各个构件的R值相等。当树干长度减小，而分枝内力和长度不变时，树干的R值会低于平均值，此时应降低树干的内力与几何长度的比值，降低内力与几何长度的比值后，节点位置上移，树干内力减小，长度增加，各个构件的R值趋于相同。

因此可利用式(14)计算得到每一个构件的内力与几何长度的比值的改变量Δq。将计算得到的Δq直接与原来的q叠加，如式(15)，当Δqi为负值时，减小qi，当Δqi为正值时，增大qi。然后根据更新后的q再次利用平衡矩阵法形。

Δq＝(R-ave(R)).*Γcr/l/ψ (14)

q＝q+Δq (15)

式中ψ为调节收敛速度的变量，为保证计算稳定，ψ为一个大于1的数，本发明中案例取值为20。

步骤3.2：在利用平衡矩阵法对树状结构形的过程中需要求解所有分枝在给定荷载作用下的内力大小。为了求解方便，同时考虑大变形效应，内力的求解部分在ANSYS中进行。在Matlab中调用ANSYS，根据节点的初始坐标建立模型并求解，然后提取所有分枝的内力。根据式(13)计算矢量R，随后根据式(15)计算所有分枝的考虑稳定承载力以后的力密度。

依据所得到的最新的力密度矢量，利用平衡矩阵法再次进行形分析，调节各个分枝的几何长度。以此循环，直到Δq足够小为止。

另外，树状结构的树干一般都是竖直方向，因此在建立树状结构数值模型时，将树干上下端节点的水平坐标设置为一致。虽然树干的下部节点为约束节点，但是在迭代过程中，下节点的水平坐标随着树干顶部节点的水平移动而移动。在迭代过程中树干底部节点的竖向坐标保持不变。顶部节点的所有坐标保持不变。

具体实施算例

步骤1：平面树状结构

步骤1.1：根据初始结构概况，所述本节利用所提出的优化算法首先对一平面树状结构进行了分析。为简化计算，将同一级所有分枝的计算长度系数设置为相同的数值。基于所提出的优化方法，对平面树状结构的形优化进行了系统分析。将各级分枝的截面特性及计算长度系数设置为不同的参数，树干及各级分支的截面参数如表1所示，SI代表树状结构截面参数组合的编号。其中，IG、AG和Ak、Ik分别表示树干和第k级分枝的轴惯性矩和截面面积；各级分支的计算长度系数如表2所示，μI、μII、μIII、μIV表示四种各分枝计算长度系数组合编号，μG和μk分别表示树干和第k级分枝的计算长度系数。所采用的树状结构的荷载及初始形状以及构件编号如图5所示，该初始形状为人为假设，共有4级分枝，受到不均匀荷载的作用，总高为1.4m，总宽为1.5m。

表1不同情况下平面树状结构截面参数

表2不同情况下平面树状结构计算长度系数

步骤1.2：根据分析结果对比验证，所述为了验证本发明所提出的树状结构形优化算法的可靠性，本节将分析结果与作者提出的虚拟温度法的分析结果进行对比分析，已验证本发明算法的可靠性。计算基于计算长度系数μI和μIII展开，计算结果对比如图6所示。从对比结果可以看出，基于本发明改进平衡矩阵法的计算结果与虚拟温度法计算结果总体上吻合较好。个别构件由于数值计算的误差会存在轻微差别。例如长度系数μIII的对比结果，树状结构边缘分枝的计算误差较大。引起误差的原因是无论是虚拟温度算法还是本发明中的改进平衡矩阵算法，都是从总体均值的角度对树状结构进行优化。根据树状结构的受力特点，无法做到使每一个分枝的R值相等。因为同一节点会与同一级分枝的多个构件相连，无法同时满足所有构件，只能使所有分枝的R值尽量接近。

图7所示为基于μIII组合所得到的最终各个构件的Δq值。从计算结果可以看出，树状结构边部分枝的Δq值基本都是正数，而且数值偏大，内部分枝的Δq值普遍是负数，数值偏小。Δq值大于0时说明该分枝的实际内力与欧拉临界荷载的比值偏小；相反，Δq值大于0时说明该分枝的实际内力与欧拉临界荷载的比值偏大。随着荷载的增加，Δq为正时的构件会首先失稳，且Δq值越大，该构件越早失稳。由于边缘构件内力较小，因此构件的稳定承载力由内部构件控制。从图6所示对比结果可以看出，对于内部构件来说，两种算法得到的结算结果高度吻合。另外，图7所示Δq的值与其内力相比普遍较小，基本可以忽略。

步骤1.3：为了验证本发明计算方法对不同荷载幅值计算结果的稳定性，同时也揭示荷载大小对形结果的影响，将图5所示荷载放大2倍，放大后的荷载用荷载I代替，放大之前的荷载以荷载I表示，将不同荷载幅值下的形结果进行对比，如图9所示。从高对比结果可以看出，荷载大小不会影响形结果，形结果与荷载位置和荷载分布相关。

将荷载II工况下每次迭代计算后所得到的Δq矢量的所有元素的绝对值和Sum(Δq)作为残差以判断形结果的收敛速度。残差随迭代次数的变化曲线如图10所示。从显示结果可以看出，残差迭代到100次左右以后就不在发生变化。所以迭代100次以后计算结果收敛，用时5分钟左右。最终残差减小为0。优化后树状结构的内力云图如图11所示。从图中所示结果可以看出，优化后所有分枝的弯矩基本减小为0，达到了形的目的。

步骤2：空间树状结构

步骤2.1：根据初始结构概况，所述在平面树状结构的基础上，本节分析了包含1365个分枝的空间树状结构，包含1366个节点。空间树状结构的初始形状如图12所示。该树状结构为4分枝，即下一级构件数量是上一级构件数量的四倍。树干与一级分枝、二级分枝、三级分枝、四级分值和五级分枝的惯性矩的比值分别为1:1:0.5:0.25:0.1:0.01。各级分支的计算长度系数如表3所示。荷载作用在顶部节点且方向竖直向下，荷载大小为5kN。约束顶部节点的水平位移以及树干底部节点的所有自由度。通过该算例的结果分析验证所提出的算法对空间受力树状结构的适用性。

表3不同情况下空间树状结构计算长度系数

步骤2.2：根据优化结果分析，所述基于所提出的优化方法对空间受力树状结构进行了形分析，优化计算迭代3000次，最终得出了不同计算长度系数组合下树状结构的优化结果，如图13所示。从计算结果可知，所提出的优化算法对于大型的空间树状结构的形优化同样适用。避免了双单元法中由于刚度系数设置不恰当而导致的计算不收敛的问题。由于树状结构的总高度不变，不同的计算长度系数的取值会改变不同级别分值的相对长度。通过调整不同级别分值的相对长度来提高树状结构的整体稳定承载力。

果可知，μ2和μ3的提高使2级分枝和3级分枝的长度减小，而其它级别分枝的长度增加。因此，本发明所提出的方法是在树状结构总高度不变，各级分枝截面特性指定的前提下进行优化。