CN109446611B - 一种强耦合树状结构找形优化设计方法 - Google Patents

Abstract

本发明提出一种强耦合树状结构找形优化设计方法，流程包括：设定初始参数；建立具有双单元数值模型的树状结构；若m为奇数，施加反向载荷；求解每个分枝内力，计算所有分枝的温度载荷；若m为偶数，施加温度载荷；提取所有节点竖向位移；根据提取点竖向位移，改变节点位置；直到最大节点位移小于阈值，结束找形过程。本发明采用双单元数值模型，杆单元的横截面积远大于梁单元，给梁单元赋予一个很小的抗弯刚度，以克服双单元树状模型由于刚度矩阵不满秩而导致的无法求解的问题；本发明具有很高的精度，很高的效率，避免了繁重的编程工作，将优化与找形两个工作耦合在一个迭代程序中，且收敛速度快，并保证树状结构整体稳定承载能力的最大化。

Description

一种强耦合树状结构找形优化设计方法

技术领域

本发明涉及建筑结构工程技术领域，特别涉及一种强耦合树状结构找形优化设计方法。

背景技术

树状结构由于美观的造型以及高效的荷载传递形式而逐渐得到了广泛的应用。树状结构可以将大范围内的竖向荷载高效的传递到地面上的一点，从而被广泛地应用于大型建筑结构当中。由于美观的造型，树状结构往往成为地标性建筑。由于其良好的力学性能，很多学者已对其进行了研究。

通过合理的安排布置树状结构各个分支的空间位置，可以是树状结构的各个分枝只受到轴向压力的作用。因此，目前关于树状结构的主要问题就是通过找形分析以寻求最合理的受力形式。目前有关树状结构找形的计算方法主要有Kolodziejczyk提出的通过将丝线模型浸在水中，利用水的表面张力作用实现树状结构找形；Buelow提出得干丝线模型方法对树状结构进行找形；Hunt提出将树状结构的全部节点看成是铰接，然后施加竖向滑动的虚拟支座，通过迭代减小虚拟支座的反力来进行找形；陈志华提出的通过连续折现索单元来进行树状结构的找形以及武岳提出的逆吊递推找形法，逆吊递推找形法生成的树状结构的所有分枝均只受轴力作用，完全满足找形的要求。

目前有关树状结构找形的方法已有很多，比较成熟。但是大多需要繁重的编程过程，很难被广大的工程技术人员所掌握。因此，积极探索更为简单实用的找形方法对于树状结构的推广应用是大有裨益的。另外，目前已有的找形方法没有对各个分枝的几何长度进行优化。众所周知，找形后树状结构的各个分枝只受到轴向压力的作用，在压力作用下，分枝会出现失稳。所有分枝的最小稳定承载能力确定了整个树状结构的稳定承载能力。因此，对杆件的稳定承载能力进行优化可以提高整体结构的稳定承载能力，而目前有关该项研究的成果很少。

发明内容

为了解决上述需要繁重的编程过程还是很难被广大的工程技术人员所掌握的问题，本发明提出了利用双单元数值模型的迭代过程对树状结构进行找形分析的数值分析方法，在此基础上，利用各个分枝在给定荷载作用下的内力对其几何长度进行优化，提出了集找形与优化为一体的高效迭代算法-树状结构强耦合找形优化算法。一种强耦合树状结构找形优化设计方法，包括如下流程：

步骤1：根据设计需要设定初始参数，包括：树状结构高度、截面尺寸、总分级数，根据设计或者实际需要设置施加荷载的位置，施加荷载的节点坐标在迭代过程中是固定的；根据节点刚度设置每个单元的计算长度系数μ，以下公式中的μ_i是指第i个分枝的计算长度系数；

步骤2：建立双单元数值模型，根据初始参数，组成具有双单元数值模型的树状结构，树干是竖直方向的，树状结构的每个分支为一个双单元数值模型，双单元数值模型包括梁单元与杆单元，设置梁单元的截面面积a与杆单元的截面积A，设置梁单元的轴惯性矩I；

所述建立双单元模型的方法为：在有限元软件中先建立普通的梁单元数值模型，并将梁单元复制，将复制后的梁单元类型改为杆单元模型，但节点号不变，将梁单元与杆单元一起称为双单元数值模型。

每个双单元数值模型中杆单元的截面面积为梁单元的100～10000倍。

步骤3：设置正整数m初始值为1；

步骤4：判断若m为奇数，向具有双单元数值模型树状结构施加反向载荷，即竖直向上的点载荷；

步骤5：对施加反向载荷的具有双单元数值模型树状结构进行非线性静力分析，求解每个分枝内力，其中，第i分枝的轴向内力为F_i；提取所有节点水平位移，其中，第j节点的水平位移为d_xj和d_yj；提取所有分枝的几何长度，第i个分枝的几何长度为l_i，根据各个分枝的轴向轴力按下列步骤计算所有分枝的温度载荷△T，包括步骤5.1～步骤5.5：

步骤5.1：根据压杆稳定及欧拉承载力公式得到优化后每个分枝的稳定极限承载能力F_cri：

其中，E_i是指第i个分枝所用材料的弹性模量，I_i为第i个分枝的梁单元的轴惯性矩，l_i为第i个分枝的几何长度，μ_i是指第i个分枝的计算长度系数；

步骤5.2：计算分枝在给定载荷作用下的实际内力与稳定承载力的比值R：

其中，R_i为第i个分枝在给定载荷作用下的实际内力与稳定承载力的比值，F_i为第i分枝的轴向内力；

步骤5.3：计算分枝在给定载荷作用下的实际内力与稳定承载力的比值的平均值R_ave：

其中，n为所有分枝的数量；

步骤5.4：计算R_i与R_ave的差值△R_i，当△R_i为负数时，即R_ave小于R_i，则分枝的内力大，则施加负温度载荷，以减小分枝几何长度，否则施加正温度载荷，以增加分枝几何长度：

△R_i＝R_ave-R_i (4)

步骤5.5：计算分枝的温度载荷△T：

其中，A_i为第i个分枝杆单元的截面积，α_i为第i个分枝的线膨胀系数；

步骤6：根据提取节点水平位移，改变节点位置，其中，第j个节点位置改变为：x_j＝x_j+d_xj,y_j＝y_j+d_yj，x_j为第j个节点水平x轴位置，y_j为第j个节点水平y轴位置；

步骤7：判断若m为偶数，以步骤5中计算的温度载荷数值，对具有双单元数值模型树状结构施加温度载荷；

步骤8：在有限元软件中通过施加位移约束的方式，约束树状结构顶部节点竖向及水平自由度，同时约束树干底部节点的所有自由度。

步骤9：进行非线性静力分析，提取所有节点竖向位移，其中第j个节点竖向位移为：d_zj；

步骤10：根据提取点竖向位移，改变节点位置，其中第j个节点位置改变为：z_j＝z_j+d_zj，z_j为第j个节点竖直z轴位置；

步骤11：判断最大节点位移是否小于阈值，即Max(d_xj,d_yj,d_zj)<error，其中，error为阈值，如果是，则结束找形过程，否则，m＝m+1，返回步骤4重新进行找形。

有益技术效果：

1、由于双单元数值模型的每个线单元由两个单元组成，即杆单元和梁单元。杆单元的横截面积远大于梁单元，给梁单元赋予一个很小的抗弯刚度，以克服双单元树状模型由于刚度矩阵不满秩而导致的无法求解的问题。

2、本发明基于双单元的找形迭代程序具有很高的精度，同时也具有很高的效率，避免了繁重的编程工作。

3、根据压杆的欧拉稳定承载能力对各个分枝的长度进行了优化，则可保证树状结构整体稳定承载能力的最大化。

4、由找形优化迭代程序，经过反复迭代以调整每个分枝的长度，从而使各个分枝的内力比趋于相同。

5、本方法可以将优化与找形两个工作耦合在一个迭代程序中，而且收敛速度快。

附图说明

图1为本发明具体实施中采用双单元示意图；

图2是本发明具体实施中一种强耦合树状结构找形优化设计方法流程图；

图3是本发明具体实施中算例1的树状结构初始形状及荷载示意图；

图4是本发明具体实施中算例1的找形后树状结构内力示意图；

其中，(a)为算例1的找形后树状结构弯矩示意图；(b)为算例1的找形后树状结构轴力示意图；

图5是本发明具体实施中算例1的找形后树状结构对比示意图；

其中，(a)为本发明所得结果；(b)为Hunt所得结果；

图6是本发明具体实施中树状结构分级及几何长度示意图；

图7是本发明具体实施中树状结构初始形状示意图；

图8是本发明具体实施方式中优化过程中树状结构形态变化示意图；

其中，(a)为第10次；(b)为第50次；(c)为第100次；(d)为第200次；(e)为第300次；(f)为第500次；(g)为第1000次；(h)为第2000次；(i)为第6000次；(j)为第10000次；

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

其中，(a)为树状结构弯矩云图；(b)为树状结构轴力云图；

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

图11是本发明具体实施中不同计算长度系数下树状结构优化结果示意图；

其中，(a)s_I和μ_I；(b)s_I和μ_II；(c)s_I和μ_III；(d)s_I和μ_IV

图12是本发明具体实施中空间树状结构示意图；

其中，(a)为空间树状结构前视图；(b)为空间树状结构俯视图；

图13是本发明具体实施中空间树状结构内力云图；

其中，(a)为空间树状结构轴力云图；(b)为空间树状结构弯矩云图；

图14是本发明具体实施中不同计算长度系数下树状结构优化结果示意图；

其中，(a)为s_I和μ_I的俯视图、前视图、透视图；(b)为s_I和μ_II的俯视图、前视图、透视图。

具体实施方式

下面结合附图和具体实施实例对发明做进一步说明：一种强耦合树状结构找形优化设计方法，如图2所示，包括如下流程：

步骤1：根据设计需要设定初始参数，包括：树状结构高度、截面尺寸、总分级数，本实施例分为4级，根据设计或者实际需要设置施加荷载的位置，施加荷载的节点坐标在迭代过程中是固定的；根据节点刚度设置每个单元的计算长度系数μ，以下公式中的μ_i是指第i个分枝的计算长度系数；

步骤2：建立双单元数值模型，根据初始参数，组成具有双单元数值模型的树状结构，树干是竖直方向的，树状结构的每个分支为一个双单元数值模型，双单元数值模型包括梁单元与杆单元，如图1所示，设置梁单元的截面面积a与杆单元的截面积A，a＝A/(100～10000)，设置梁单元的轴惯性矩I；基于上述所提出的迭代程序对图3所示树状结构进行了找形分析。该树状结构所用双单元数值模型中，杆单元的截面面积设为1×10^-3m²，梁单元的截面面积设为1×10^-6m²，梁单元的轴惯性矩为2×10^-10m⁴。

所述建立双单元模型的方法为：在有限元软件中先建立普通的梁单元数值模型，并将梁单元复制，将复制后的梁单元类型改为杆单元模型，但节点号不变，将梁单元与杆单元一起称为双单元数值模型。

每个双单元数值模型中杆单元的截面面积为梁单元的100～10000倍。杆单元的横截面积远大于梁单元，给梁单元赋予一个很小的抗弯刚度，以克服双单元树状模型由于刚度矩阵不满秩而导致的无法求解的问题。由于具有很小的抗弯刚度，双单元数值模型可以很好的模拟树形结构。

步骤3：设置正整数m初始值为1；

步骤4：判断若m为奇数，向具有双单元数值模型树状结构施加反向载荷，即竖直向上的点载荷，10kN；

步骤5：对施加反向载荷的具有双单元数值模型树状结构进行非线性静力分析，求解每个分枝内力，如图4所示，其中，第i分枝的轴向内力为F_i；提取所有节点水平位移，其中，第j节点的水平位移为d_xj和d_yj；提取所有分枝的几何长度，第i个分枝的几何长度为l_i，根据各个分枝的轴向轴力按下列步骤计算所有分枝的温度载荷△T，包括步骤5.1～步骤5.5：

步骤5.1：根据压杆稳定及欧拉承载力公式得到优化后每个分枝的稳定极限承载能力F_cri：

其中，E_i是指第i个分枝所用材料的弹性模量，I_i为第i个分枝的梁单元的轴惯性矩，μ_i为第i个分枝的计算长度系数，l_i为第i个分枝的几何长度；

由于树状结构往往都采用圆形截面，因此可不比区分强轴和弱轴。根据式(1)可知，分枝的稳定极限承载能力与分枝的截面特性，几何长度，计算长度系数以及材料属性密切相关。加入树状结构各个分枝同时失稳，则可保证树状结构整体稳定承载能力的最大化。在此，设分枝在给定荷载作用下的实际内力与稳定承载力的比值为R，如(2)所示，则当各个分枝的比值相同时，则保证树状结构整体稳定承载能力的最大化。由于分枝截面往往根据设计已经确定，因此E_i可以视为常数，计算长度系数与节点的抗弯刚度相关，对于特定的节点，计算长度系数可视为定值。因此，R值的大小可以通过分枝的几何长度来改变。

步骤5.2：计算分枝在给定载荷作用下的实际内力与稳定承载力的比值R：

其中，R_i为第i个分枝在给定载荷作用下的实际内力与稳定承载力的比值，F_i为第i分枝的轴向内力；

步骤5.3：计算分枝在给定载荷作用下的实际内力与稳定承载力的比值的平均值R_ave：

其中，n为所有分枝的数量；对于R值偏大的分枝，则属于是较危险分枝，可以通过减小分枝的几何长度l来减小R的大小。同样，对于R值较小的分枝，则可通过增加l的大小来增大R值。在双单元数值模型中，分枝几何长度的大小可以通过施加温度荷载的方式改变。由于温度在实际中并不存在，因此称为虚拟温度法。

步骤5.4：计算R_i与R_ave的差值△R_i，当△R_i为负数时，即R_ave小于R_i，则分枝的内力大，则施加负温度载荷，以减小分枝几何长度，否则施加正温度载荷，以增加分枝几何长度：

△R_i＝R_ave-R_i (4)

步骤5.5：计算分枝的温度载荷△T：

其中，A_i为第i个分枝杆单元的截面积，α_i为第i个分枝的线膨胀系数；

根据“木桶理论”，树状结构的整体稳定承载能力由稳定承载能力最小的分枝决定。因此，可以说当树状结构各个分枝同时失效时，树状结构的节点定位，分枝长度将是最为合理的。因此，根据压杆的欧拉稳定承载能力对各个分枝的长度进行了优化理论的研究。树状结构的结构分级及分枝几何长度如图6所示。

步骤6：根据提取节点水平位移，改变节点位置，其中，第j个节点位置改变为：x_j＝x_j+d_xj,y_j＝y_j+d_yj，x_i为第j个节点水平x轴位置，y_j为第j个节点水平y轴位置；

步骤7：判断若m为偶数，以步骤5中计算的温度载荷数值，对具有双单元数值模型树状结构施加温度载荷；

步骤8：在有限元软件中通过施加位移约束的方式，约束树状结构顶部节点竖向及水平自由度，同时约束树干底部节点的所有自由度。

步骤9：进行非线性静力分析，提取所有节点竖向位移，其中第j个节点竖向位移为：d_zj；

步骤10：根据提取点竖向位移，改变节点位置，其中第j个节点位置改变为：z_j＝z_j+d_zj，z_j为第j个节点竖直z轴位置；

步骤11：判断最大节点位移是否小于阈值，即Max(d_xj,d_yj,d_zj)<error，其中，error为阈值，如果是，则结束找形过程，否则，m＝m+1，返回步骤4重新进行找形。

为简化计算，本实施例将同一级所有分枝的计算长度系数设置为相同的数值。基于所提出的优化方法，本节对平面树状结构的找形优化进行了系统分析。将各级分枝的截面特性及计算长度系数设置为不同的参数，树干及各级分支的截面参数如表1所示。其中，I_G，A_G和A_i，I_i分别表示树干和第i级分枝的轴惯性矩和截面面积；各级分支的计算长度系数如表2所示。同理，μ_G和μ_i分别表示树干和第i级分枝的计算长度系数。所采用的树状结构的荷载及初始形状如图7所示。

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

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

假设树状结构的截面特性和计算长度系数分别采用S_I和μ_I，在给定荷载作用下进行了找形优化分析。图8所示为经过不同迭代次数后，树状结构的形状。从所示计算结果可以看出，通过采用所提出的方法，不合理的分枝长度可以被快速的优化。由于树干的截面面积和轴惯性矩较大，且计算长度系数小，因此经过迭代后树干的长度减小。由于树状结构边缘分枝的内力较小，因此优化后，边缘分枝的几何长度普遍大于树状结构内部分枝的几何长度。同时还可以看出，在经过2000次迭代后，树状结构的形状几乎不再改变，说明本方法同时作为树状结构的找形和优化方法是非常高效的。

图9所示为优化后树状结构的内力云图。从图中可以看出，在优化分枝长度的同时，本方法同时还可以为树状结构找形，找形后分枝的弯矩数量级约是轴力的1％，基本可以忽略。本方法可以将优化与找形两个工作耦合在一个迭代程序中，而且收敛速度快，图10所示为节点18，位置见图8，在收敛过程中x向和y向位移的变化过程。从图中可以看出，节点18的位移很快减小为0，在大约迭代2500次后，节点18的位置不再变化，整个树状结构的优化找形过程已经完成。

将表1和表2所示的不同参数进行了组合，分别得到了不同参数组合下树状结构的最优结果，如图11所示。图11a表明，当各级分枝的计算长度系数与树干相同时，树干的长度要明显小于其它情况。从图11c可以看出，随着第4级分枝计算长度系数的增大，该级分枝的几何长度明显减小。图11d也表明随着第1级分枝计算长度系数的相对增大，第1级分枝的几何长度也明显减小。因此，计算结果表明，计算长度系数的影响可以在本计算方法中得到精确的反应。

在实际工程中，树状结构都是空间的。为了验证本发明所提方法对于空间树状结构的适用性，本节对空间树状结构的优化找形进行了分析，所建柱状结构如图12所示。该树状结构的所有顶点不再一个平面内，而是位于一个近似球形的曲面内。所有节点施加10kN的竖直方向的集中力。不同情况下各个分枝的计算长度系数和截面特性如表3和表4所示。

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

表4不同情况下空间树状结构截面参数

优化后树状结构的内力云图如图13所示。从计算结果可以看出，树状结构轴力最大约是文具最大值的两万分之一。弯矩值相对于轴力很小，以致可以忽略，认为分枝之手轴力的作用。

不同情况下空间树状结构优化找形后的结果如图14所示。从分析结果可以看出，所提出的方法对于三维空间在树状结构同样可以快速高效的得到精确的结果，具有较强的工程适用性。基于S_I和μ_I所得结果的一级分枝长度要明显小于基于S_I和μ_II所得到的一级分枝长度，从不同的结果对比可以看出，分枝的计算长度系数同样可以在找形结果中得到精确的反应。

本发明对树状结构的找形优化进行了较为系统的研究。在利用双单元建立的树状结构数值模型的基础上，提出了树状结构的找形迭代程序，结果显示该迭代程序可以高效快速的得到满足要求的树状形态。在给定荷载下，该树状结构的各个分枝只受到轴力的作用，弯矩很小可以忽略。

在找形的基础上，本发明基于压杆欧拉稳定承载力对树状结构的分枝长度进行了优化分析。提出了基于虚拟温度法的优化程序，并将找形程序与优化程序结合提出了树状结构的强耦合找形优化算法。本方法不但可以保证树状结构的各个分支在给定荷载作用下只受到轴力的作用，同时能根据各个分支的截面特性和内力对其长度进行优化，进而提高整体树状结构的稳定承载能力。分析结果表明，该算法可以高效的应用于空间树状结构的找形分析。

实验结果：

将本发明找形后所得树状结构的结果与Hunt所得结果进行了对比研究，其中Hunt为论文：Hunt J,Haase W,Sobek W.A design tool for spatial tree structures[J].Journal of the International Association for Shell and Spatial Structures,2009,50(1):3-10，如图5所示。从计算结果可以看出，两个高度吻合，说明本发明基于双单元的找形迭代程序具有很高的精度，同时也具有很高的效率，避免了繁重的编程工作。

Claims (3)

1.一种强耦合树状结构找形优化设计方法，其特征在于，包括如下流程：

步骤1：根据设计需要设定初始参数，包括：树状结构高度、截面尺寸、总分级数，根据设计或者实际需要设置施加荷载的位置，施加荷载的节点坐标在迭代过程中是固定的；根据节点刚度设置每个单元的计算长度系数μ；

步骤2：建立双单元数值模型，根据初始参数，组成具有双单元数值模型的树状结构，树干是竖直方向的，树状结构的每个分支为一个双单元数值模型，双单元数值模型包括梁单元与杆单元，设置梁单元的截面面积a与杆单元的截面积A，设置梁单元的轴惯性矩I；

步骤3：设置正整数m初始值为1；

步骤4：判断若m为奇数，向具有双单元数值模型树状结构施加反向载荷，即竖直向上的点载荷；

步骤5：对施加反向载荷的具有双单元数值模型树状结构进行非线性静力分析，求解每个分枝内力，其中，第i分枝的轴向内力为F_i；提取所有节点水平位移，其中，第j节点的水平位移为d_xj和d_yj；提取所有分枝的几何长度，第i个分枝的几何长度为l_i，根据各个分枝的轴向轴力按下列步骤计算所有分枝的温度载荷△T，包括步骤5.1～步骤5.5：

步骤5.1：根据压杆稳定及欧拉承载力公式得到优化后每个分枝的稳定极限承载能力F_cri：

其中，E_i是指第i个分枝所用材料的弹性模量，I_i为第i个分枝的梁单元的轴惯性矩，l_i为第i个分枝的几何长度，μ_i是指第i个分枝的计算长度系数；

步骤5.2：计算分枝在给定载荷作用下的实际内力与稳定承载力的比值R：

其中，R_i为第i个分枝在给定载荷作用下的实际内力与稳定承载力的比值，F_i为第i分枝的轴向内力；

步骤5.3：计算分枝在给定载荷作用下的实际内力与稳定承载力的比值的平均值R_ave：

其中，n为所有分枝的数量；

步骤5.4：计算R_i与R_ave的差值△R_i，当△R_i为负数时，即R_ave小于R_i，则分枝的内力大，则施加负温度载荷，以减小分枝几何长度，否则施加正温度载荷，以增加分枝几何长度：

△R_i＝R_ave-R_i (4)

步骤5.5：计算分枝的温度载荷△T：

其中，A_i为第i个分枝杆单元的截面积，α_i为第i个分枝的线膨胀系数；

步骤6：根据提取节点水平位移，改变节点位置，其中，第j个节点位置改变为：x_j＝x_j+d_xj,y_j＝y_j+d_yj，x_j为第j个节点水平x轴位置，y_j为第j个节点水平y轴位置；

步骤7：判断若m为偶数，以步骤5中计算的温度载荷数值，对具有双单元数值模型树状结构施加温度载荷；

步骤8：在有限元软件中通过施加位移约束的方式，约束树状结构顶部节点竖向及水平自由度，同时约束树干底部节点的所有自由度；

步骤9：进行非线性静力分析，提取所有节点竖向位移，其中第j个节点竖向位移为：d_zj；

步骤10：根据提取点竖向位移，改变节点位置，其中第j个节点位置改变为：z_j＝z_j+d_zj，z_j为第j个节点竖直z轴位置；

步骤11：判断最大节点位移是否小于阈值，即Max(d_xj,d_yj,d_zj)<error，其中，error为阈值，如果是，则结束找形过程，否则，m＝m+1，返回步骤4重新进行找形。

2.根据权利要求1所述一种强耦合树状结构找形优化设计方法，其特征在于，所述建立双单元模型的方法为：在有限元软件中先建立普通的梁单元数值模型，并将梁单元复制，将复制后的梁单元类型改为杆单元模型，但节点号不变，将梁单元与杆单元一起称为双单元数值模型。

3.根据权利要求1所述一种强耦合树状结构找形优化设计方法，其特征在于，每个双单元数值模型中杆单元的截面面积为梁单元的100～10000倍。

Images