找形分析、荷载分析和裁剪样式生成是薄膜结构设计的主要阶段.薄膜结构经过找形分析, 获得了由离散点描述的膜曲面.膜曲面裁剪样式生成是膜结构设计中的关键[1,2,3], 其目的是将已经由找形分析得到的三维膜曲面转化成离散的平面膜片, 在市售的膜卷材上进行下料分析.定义膜片首先要定义膜片沿经向的裁剪线.当裁剪线确定以后, 可以沿着裁剪线将曲面裁成离散的膜片, 但每一膜片来说仍是三维曲面, 因此还需将空间膜片展开成平面膜片.

对于曲面, 已经证明的可展面只有3类[4]: 柱面、锥面, 以及曲线的切线曲面.展开后, 空间曲面和展开平面上曲线的弧长不变是曲面的可展性标准.薄膜结构找形获得的曲面绝大多数是不可展曲面, 即不可能无误差展开.针对膜面的展开, 存在着不同的方法, 如王启文等人[5]采用了无约束极值法, 此方法是一种比拟的方法, 对初始形状的依赖性较强;Xia等人[6]和卫东等人[7]分别提出了刚性展开的几何法和平面映射法, 即按照三角形边长相等原则将膜面展开, 这是一种基于约束的三角片翻转的方法, 对于可展曲面是可以采用的, 但对于不可展曲面, 由一个三角形单元确定的第三个节点位置和其他单元确定的此节点位置不会重合, 即使采用平均位置也可能导致层叠误差的发生;钱基宏等人[8]采用最小变形能原理推导了几何非线性有限元求解算法, 这种算法在节点较多时需要求解大型方程, 在一定程度上影响计算速度.

在飞机、汽车、服装和造船等行业, 复杂曲面的展开问题也是计算机辅助设计研究的热点和难点.Parida等人[9]提出了将离散三角化片逐一变换到指定平面上的算法, 但在展开时容易产生裂纹和较大误差;Shimada等人[10]使用有限元算法计算曲面展开, 基于曲面离散表达式, 运用分区域曲面的展开思想进行曲面展开, 但所使用的曲面分区域算法过于简单, 无法处理复杂曲面的展开问题.

本文针对薄膜结构的曲面展开, 提出一种基于弹簧-质点系统的动力求解方法.在膜片展开过程中.平面膜片的局部精度可以通过弹簧弹性变形系数进行控制, 易满足膜片边界相容的要求.

欲将空间的膜片展平到平面上, 获得二维膜片, 首先需要在选定的平面上定义一个和空间膜片在参数域与空间域具有相同拓扑结构的二维膜片.对于平面的选择, 如果曲率不太大, 可以直接向x-y平面上投影来获得膜片的平面拓扑形式, 但必须保证不产生自层叠现象;如果产生了自层叠现象, 则可选择另外一个平面进行投影, 以保证不产生自层叠的平面投影.采用这种做法尽管初始映射的结果非常粗糙, 但三角网格在参数域与空间域具有相同的拓扑结构.在获得了平面拓扑后, 为了加快收敛速度, 还可以修正初始平面膜片, 引入比例因子S将平面膜片进行比例变换:

S=max (Lq/Lp) . (1)

式中: Lq为描述空间膜片的三角形的边长, Lp为描述平面膜片的三角形边长.

由平面三角形网格可以建立弹簧-质点系统, 其中三角形的节点作为质点, 边长采用弹簧代替, 这样系统的物理量与几何量相对应, 如弹簧的力、弹性变形能可由平面网格节点间距与曲面三角形网格节点间距的差异来获得;质点质量可由包含此质点的平面膜片的单元面积来确定.空间膜片的网格形状和位置与展开后二维膜片的网格形状和位置的差异, 可视为一种存储于弹簧-质点系统中的弹性变形能.弹簧-质点系统见图1.其中, Pi为质点, 质点Pi和Pj之间采用弹簧连接.在变形过程中, 如果平面上描述二维膜片的三角形单元的节点Pi和Pj的间距大于对应的空间节点的间距, 则对质点Pi和Pj施加拉力, 否则施加压力.

定义弹簧-质点系统在质点Pi处的弹性力为

fi=∑j=1mC(|PiPj|−dj)NPiPj. (2)fi=∑j=1mC(|ΡiΡj|-dj)ΝΡiΡj.(2)

式中: C为弹簧弹性变形系数, |PiPj|为平面网格上质点Pi和Pj的间距, dj为空间膜片上Pi和Pj的间距, NPiPj为平面上Pi指向Pj的单位矢量, m为围绕质点Pi的弹簧总数.

在曲面的不同部位定义不同的弹簧弹性变形系数, 可对膜片展开的局部精度进行控制.为了满足施工要求和精度要求, 要求相邻膜片的边界线长度相等, 因此可以通过对平面膜片边界线上的弹簧赋予较大的弹簧弹性变形系数, 以使相邻膜片的裁剪线长度误差满足精度要求.

通过弹簧-质点系统的变形, 可由初始平面膜片得到对应的空间膜片展开形状.采用拉格朗日运动方程描述系统运动:

Mq⋅⋅+Dq?+Kq=gq+fq. (3)Μq⋅⋅+Dq?+Κq=gq+fq.(3)

式中: q为质点坐标, M、D和K分别为系统质量矩阵、阻尼矩阵和刚度矩阵, gq为系统内力向量, fq为系统外力向量.在进行膜片展开计算时, gq和fq为零, 忽略阻尼项, 得到简化的拉格朗日运动方程为

Mq⋅⋅+Kq=0. (4)Μq⋅⋅+Κq=0.(4)

式中: Kq为式 (2) 确定的弹性力, 即Kq=-f.

通常采用欧拉方法求解方程 (4) .当Δt很小时, 质点Pi的加速度为常量, 则整个弹簧-质点系统的平衡由每个质点的平衡组成, 质点质量矩阵采用集中质量矩阵.对于质点Pi, 在典型时刻t到t+Δt, 由方程 (4) 得到简化的拉格朗日运动方程为

mi=ρ3∑k=1nAk, (5)mi=ρ3∑k=1nAk,(5)

q⋅⋅ti=fti/mi, (6)q⋅⋅it=fit/mi,(6)

q?t+Δti=q?ti+Δtq⋅⋅ti, (7)q?it+Δt=q?it+Δtq⋅⋅it,(7)

qt+Δti=qti+Δtq?ti+Δt22q⋅⋅ti. (8)qit+Δt=qit+Δtq?it+Δt22q⋅⋅it.(8)

式中: mi为质点Pi的质量, ρ为膜材单位面积的质量, Ak为围绕质点Pi的第k个三角形单元的面积, n为围绕质点Pi的三角形的总数, q⋅⋅tiq⋅⋅it为质点Pi在时刻t的加速度, q?tiq?it和q?t+Δtiq?it+Δt分别为质点Pi在时刻t和t+Δt的速度, qtiit和qt+Δtiit+Δt分别为质点Pi在时刻t和t+Δt的坐标.

从式 (7) 可以看出, 前一时刻的速度对后一时刻的速度是有影响的, 即所谓的惯性作用.这样就会产生一种累计效应, 并在平衡位置附近达到最大速度.当加速度反向时, 质点不会立刻反向运动, 仍然继续按照原来方向运动, 这样会产生较大的振荡, 极易造成发散;同时在振荡过程中, 很难获得满足一定精度要求的节点平衡位置.

如图2所示, 一空间等边三角形ABC, 节点坐标为 A (1, 0, 0) , B (0, 1, 0) , C (0, 0, 1) , 将其投影到x-y平面上, 得到了初始平面三角形ABC′, 其中C′为坐标系原点.图3为当时间步长为0.1 s, 迭代次数为200时, 边线AC′长度L随迭代次数的振荡图.可以看出, 如果平面初始形状与最终平衡形状相差较大, 则会产生较大的振荡, 甚至造成单元反向和破坏, 平衡位置很难获得;同时由于没有考虑阻尼项, 结果极有可能导致发散.针对这个问题, 本文提出了一种忽略惯性的做法, 即认为在时刻t到t+Δt为匀速运动, 速度为时刻t+Δt的速度, 则由式 (7) 和 (8) 得

q?t+Δti=Δtq⋅⋅tiq?it+Δt=Δtq⋅⋅it, (9)

qt+Δti=qti+Δtq?t+Δtiqit+Δt=qit+Δtq?it+Δt. (10)

图4为当忽略惯性作用时, AC′计算长度L随迭代次数的变化曲线.可以看出, 这种做法可以保证质点即使穿过平衡位置后, 也会立刻反向, 再次向平衡位置逼近.当接近平衡位置时, 加速度和速度都会变小, 因此收敛性很好, 收敛速度非常快.

面积误差和形状误差可以作为展开精度的评定标准.在面积展开过程中, 二维膜片的面积会不停地改变, 最终的相对面积误差为

ES=|A-A′|/A. (11)

式中: A为膜片展开前的实际面积, A′为空间膜片展开后对应的面积.

形状误差可由展开前描述空间膜片的三角形边界总长度L与展开后描述二维膜片的三角形边界总长度L′之差来衡量.最终的相对边长误差为

EL=|L-L′|/L. (12)

计算中采用相对面积误差、形状误差或迭代次数来控制计算精度要求.

由于在计算过程中展开面积在不停地变化, 方程 (4) 中的质量矩阵和弹性力也是不停变化的, 方程表现出非线性特性, 计算步骤如下:

1) 向平面投影得到空间膜片的平面拓扑;

2) 由式 (5) 计算各节点质量mi;

3) 由式 (2) 计算各节点力;

4) 由式 (6) 、 (9) 和 (10) 计算各节点新位置;

5) 由式 (11) 计算系统的相对面积误差ES, 由式 (12) 计算系统的相对边长误差EL;

6) 当ES或EL达到收敛精度, 或迭代次数达到指定数目, 退出迭代, 否则按照2) ～5) 步继续迭代.

例 1 本算例为一个等张力的曲面, 由于对称性, 可以从中抽取四十分之一的膜片进行展开, 以获得各个膜片的裁剪样式.其中膜结构顶部直径为10 m, 底部直径为30 m, 预张力为30 N/cm, 高为9 m.图5为初始平衡形状.图6为当空间膜片采用本文提出的薄膜曲面展开算法展开时, 面积残余误差和长度残余误差随迭代次数的变化图.膜片的空间面积为68.199 m2, 单元各边界长度总和为204.853 m;展开后面积为68.199 m2, 单元各边界长度总和为204.847 m.由图6可以看出, 经过20000次迭代后面积残余误差ES=0.03%, 长度残余误差EL=0.001%, 收敛精度非常高;当迭代次数继续增加时, 精度可以继续提高.迭代次数较多主要是精度控制较小的缘故, 同时由于无大型方程组求解, 计算速度较快, 整个求解过程耗时不到1 min (CPU为Pentium Ⅳ, 2 G) .另外, 在迭代逼近过程中, 没有出现振荡和发散现象, 收敛情况很好.图7为膜片的展平图, 可以看出没有出现自层叠现象, 符合裁剪样式要求.

例 2 本算例对一等张力马鞍面, 采用本文的膜片展开算法来求解以网格作为裁剪线的膜片展开样式.图8是经、纬向为等张力的初始马鞍形曲面, 即在设计时固定两边, 并且另外两边的给定位移满足抛物线的方程, 膜面张力在经、纬向均为40 N/cm.整个膜面在平面上的投影边长分别为18和16 m.将整个马鞍面分为4块膜片.图9为马鞍面膜片的样式展开图, 面积精度和形状精度均为0.3%, 迭代次数控制在5000次.

本文提出的空间膜片展开算法具有以下特点:

(1) 可用来展开复杂的膜片, 对于可展曲面可以得到精确的展开解;对于不可展曲面可以得到有面积误差和形状误差控制的近似解.

(2) 在展开过程中通过改变弹簧弹性变形系数来控制膜片展开的局部精度, 易满足膜结构裁剪中边界相容的要求.

(3) 由于算法中不需要求解大型方程组, 且收敛速度快, 计算时间大大减少.

(4) 克服了以往刚性展开算法中的层叠现象和裂纹较多的缺点.