薄膜结构是一种新型的张力结构, 它的基本组成单元为支承柱、张拉索和覆盖的膜材。在膜结构的设计过程中必须首先确定结构的初始形状, 即“找形”[1?2]。动力松弛法是一种常用的找形方法, 最早由Hotter在1964年分析混凝土核容器时提出, 而Barnes成功地将这一方法发展百富策略白菜网于索网及膜结构的找形。动力松弛法将结构体系离散化后, 假定质量集中于结点上, 并在结点上施加激振力, 使结构产生振动, 最后跟踪体系的振动并判断其平衡位置[3]。

目前对该方法的研究较为普遍, 文献[4]从网格畸变的角度通过调整参数和算法, 改良动力松弛法。文献[5]则构造了不同的迭代格式以优化该算法。文献[6]研究了动力松弛法的各参数对其收敛性和收敛速度的影响。本文考虑到建筑师对膜面造型的变化要求, 对动力松弛法进行了改进, 使改进后的动力松弛法可以根据建筑师的不同要求, 得到更加自由丰富的膜面造型。此外, 本文还研究了动力松弛法中阻尼参数的取值问题, 提出了简单可行的参数取值方法。

动力松弛法是百富策略白菜网动力学原理将结构的静力问题转化为动力问题, 从而求解膜结构平衡状态的数值方法。在求解过程中不需要形成结构的刚度矩阵, 因而可以节约计算机内存。在求解过程中容易改变结构的边界约束, 特别适合于高强非线性问题。

动力松弛法的基本思想是:首先将膜结构离散为网格, 在假定的初始形状中给定应力分布, 从而形成结构体系的不平衡力。在不平衡力的作用下结构会产生运动, 逐步跟踪结构各个离散点的运动, 最终结构的运动停止在其平衡位置, 即可得到找形结果 (平衡曲面) 。

索膜结构体系的总势能为:

式中:φ为结构总势能;U为结构应变能;V为外力势能。

当结构离散为网格后, 式 (1) 可表示为:

式中, {P}和{d}分别为结构的外荷载向量和节点位移向量。

式 (2) 中对节点位移dij (节点j的i方向) 取偏导数:

式中:Pij为作用于节点j在i方向的外力矢量;dij表示j节点在i方向的位移。

当结构在外力作用下处于平衡状态时, 结构的总势能最小, 即 ;而当结构没有达到平衡状态时, 则 , 即结构在各个节点存在不平衡力Ri:

结构在不平衡力的作用下即会产生振动。根据达朗贝尔原理, 结构振动中的任一时刻其运动平衡方程为:

式中:mi为结构节点虚拟集中质量;ci为结构节点阻尼系数;di为结构节点的加速度;di为结构节点的速度。

将式 (5) 用中心有限差分形式表示:

令 , 解式 (6) , 得到:

则结构在 时刻各个节点的坐标为:

通过式 (7) 和式 (8) 即可跟踪结构各个时刻的振动状态。

百富策略白菜网动力松弛法时, 时间步长∆t的选择对算法的收敛性有较大影响。如果∆t的取值大于应力波在结构节点间传递的最短时间, 则算法不收敛[7]。阻尼系数ci的取值不会影响到算法的收敛性, 但影响收敛速度。当ci为临界阻尼时, 算法的收敛速度最快。文献[2]给出了临界阻尼公式:

式中:f为结构的基频,  ;N为结构振动一周期的迭代次数, 一般假定ci=0, 对结构试算一次, 即可得到N;im为结构各个节点的虚拟集中质量。

文献[7]中对索网结构的动力松弛法参数做了如下改进:采用Barnes在文献[8]中给出的保证解收敛的不等式, 即∆t的上限值:

式中Si是结构节点i的刚度。

取时间步长∆t=1, 则:

根据式 (9) 和式 (11) , 可以得到:

传统的动力松弛法适用性不强, 对初始的假定形状要求比较高, 而且对找形结果很难进行再调整。文献[7]对其进行了一定程度的改进, 显著增强了对索网结构找形的适用性, 但需要试算结构的基频以得到结构的临界阻尼。本文通过人为调整膜单元各边内力和改进节点阻尼的计算方法, 不需要试算即可直接得到结构临界阻尼的近似值, 并使得改进后的动力松弛法更容易控制其找形结果。

本文对阻尼的改进以文献[7]为基础, 通过推导膜单元的节点刚度, 重新计算阻尼。

文献[7]给出短程线索网结构节点刚度Si的计算公式:

式中:m是和节点i相连的单元数;Tj是单元j的拉力;Lj是单元j的长度。

将式 (13) 推广到图1, 在三角形单元j中节点刚度Si可以改写为:

对三角形膜单元j的节点i建立力平衡方程[9], 则Sij可进一步写成:

式中σj为单元j的膜面应力。

考查单自由度体系的振动, 其临界阻尼 , k为刚度。将式 (11) 代入, 可得各个节点阻尼系数:

根据式 (15) , 则节点不平衡力{iR}相应为:

式中{X i}为节点i的坐标向量。

等应力膜面因为应力处处相等, 所以对材料的利用最充分, 是最为结构工程师接受的一种膜面形状。但是在某些情况下, 我们需要用到不等应力膜面。

例如, 考察悬链线的方程:

当r<b时, 悬链线是不存在的, 即这种边界条件的等应力膜面是不存在的, 这时, 我们只能进行不等应力膜面找形。

在文献[1]中, 通过引入膜面材料弹性模量的方法, 使结构具有一定刚度, 即设定初始应力为σ0, 按下式迭代计算各个时间点的应力:

式中:[D]为单元的材料矩阵;{∆ε}为单元应变。由于这种方法需要重新计算应力和不平衡力, 所以本文从另外的角度入手。

式 (15) 中, 每个单元中的节点刚度Si是与节点相连的两条边的刚度之和。在等应力膜面找形中, 每次迭代后, 需要按照节点的新坐标和设定的找形应力σ0重新计算刚度。在不等应力膜面找形中, 以初始设定应力σ0为已知条件, 按下式计算j单元中边l的初始刚度lkj (如图2) :

klj=σ0t2 tanαlj (19)

在后面的找形迭代过程中, 每条边的刚度lkj保持不变, 节点刚度Si和节点不平衡力{iR}分别按下式计算:

式中:n为和节点i相关联的边数;lk为第l条边的刚度,  或m′=2, 即与三角形膜单元l边相关联的单元数为1或2, lkj见式 (19) 。在实际工程中, 建筑师们通常会对膜结构的造型提出一定要求。图3是一个标准的马鞍形膜面及其轮廓示意。有些情况下, 建筑师会要求调整脊线1L和谷线2L的曲率, 即抬高P点, 如图4所示, 虚线为抬高P点后的结构轮廓。从力学上考虑, 为了要达到这种效果, 需要增加脊线1L的内力或者降低谷线L2的内力。在索网的找形中, 这种需要可以直接通过改变脊索和谷索的预应力值来满足建筑师的要求。然而在膜结构找形中, 需要进行更为复杂的算法处理, 这也是本文另一重要改进之处, 即通过人为调整各单元内力, 改变最终的找形形状, 以满足建筑师的要求。

分析图4, 不难看出P点抬高后, 脊线1L的长度缩短, 内力变大;谷线2L的长度变长, 内力变小。同时考虑到脊线和谷线的内力永远不会小于0, 即内力和长度之间的关系趋势, 如图5。

按照图4、图5, 本文引入函数 以调整三角形单元各个边l的内力, 从而得到建筑师要求的初始平衡曲面形状:

式中:lk′为各个边的等效刚度,  α>0为控制参数, 其取值根据建筑师要求的初始平衡曲面的形状而定, 带有试算性质;d为网格划分控制尺寸;lL为三角形单元边长。

则式 (20) 和式 (21) 变化为:

如图6, 已知三角形单元中k1、k2、k3, 需要求解三角形单元应力。定义单元局部坐标系中, x向应变为εx, y向应变为εy, 剪应变为γxy, 三角形各个边的应变 , 则有以下换算关系:

由上式可得:

式中:

式 (25) 可写成:

由虚功原理可得:

式中:A为三角形单元面积;t为单元厚度;{T}为各边内力向量。

将式 (27) 代入 可得:

即:

式中:

根据以上理论, 可以按照以下步骤进行膜结构初始形状的确定:

1) 假定结构初始形状, 用三角形单元离散结构, 设定膜面应力σ0和索拉力0T。

2) 进行不等应力膜面找形时, 分别按式 (16) 、式 (23) 和式 (24) 计算节点阻尼ic、刚度Si和不平衡力{iR};进行等应力膜面找形时, 分别按式 (15) 、式 (16) 和式 (17) 计算Si、ic和{iR}。

3) 按式 (8) 计算下一时刻节点的坐标, 重复步骤2) , 以不平衡力{iR}作为收敛控制条件, 控制迭代次数。需要注意的是, 对于不等应力膜面, 每次迭代中的节点刚度Si保持不变, 即重复步骤2) 时无需重新计算Si和ic, 迭代结束后利用式 (29) 计算膜面应力;对于等应力膜面, 则需要在每次迭代过程中按式 (15) 计算节点刚度Si, 按式 (16) 计算节点阻尼ic。

算例1.图7所示的菱形平面双曲抛物面索网, 周边刚性边界, 设定索初始拉力为800 k N。给定收敛精度为0.1 k N。文献[7]中将中心点向上抬高0.5 m (如图7) , 给出了找形中迭代次数与系数N (N可以调节阻尼大小, 详见式 (9) ) 的关系, 见表1。N=100时, 对应临界阻尼, 迭代次数最少。当采用改进后的阻尼参数 (式 (16) ) , 迭代次数为68次, 接近文献[7]中的临界阻尼取值时的迭代次数, 图8为找形后的曲面形状。

算例2.马鞍面膜结构平面尺寸10m×10m, 初始形状及网格划分如图9, 柔性边界, 进行等应力膜面找形。膜面应力取1k N/m, 边索拉力取30k N, 控制误差0.01k N, 最终迭代305次。该膜结构膜面的理论解析解为 , 数值解与理论解比较, 最终误差在0.8%以内 (如图10) 。

算例3.等应力悬链面膜结构的解析方程为[10]:

式中, b=10m, h=17.627m, r=10m?30mm, 边界条件为上下环边固定约束, 如图11。等应力膜面找形, 应力取1k N/m, 控制误差0.01k N, 最终迭代259次。数值解与理论解比较, 最终误差在1.8%以内。找形后平衡形状为图12。

算例4.假定初始形状如图13, 上下环周边约束, 进行不等应力膜面找形。设定初始迭代应力1k N/m, 控制误差0.01k N, 网格划分尺寸d=100mm。图14?图21分别为α=2、α=3 (见式 (22) ) 后结构找形结果。

从图14、图18及算例3的悬链面方程可以看出, 此时所找到的平衡曲面已经不再是悬链面, 而是一个光滑的不等应力曲面。用水平面去截图14、图18, 所得到的曲线有些是圆, 而有些则是椭圆 (图15、图19) 。本文即通过选取不同的α值, 得到不同形状的初始平衡曲面, 以满足建筑师的要求, 这是经典动力松弛法无法做到的。

(1) 通过对典型马鞍面的轮廓分析, 以改变膜曲面造型为目标, 通过引入幂函数 修正了动力松弛法中节点刚度的取值, 使动力松弛法可以根据建筑师的要求, 得到曲率变化各不相同的找形结果。算例通过不同的参数取值 (α=2, α=3) , 成功调整了膜结构最终的找形形状, 验证了改进后动力松弛法的正确性。

需要指出的是, 所引入的幂函数表达式只是一个概念函数。作者的主要目的是提出这样一个思路:在找形过程中可以通过人为调整各三角形膜单元内力, 而调整最终找形结果, 使所找到的平衡曲面更尽人意。

(2) 讨论了动力松弛法中阻尼系数的取值, 提出了相对简便的阻尼计算公式。通过和文献[6]中的算例对比, 证实了本文所提出的阻尼计算方法和结构的临界阻尼相接近, 却避免了文献[6]中的试算过程。