以攀达穹顶（Pantadome）^[1-2]为代表的“机构施工法”将网格结构划分为若干片区组成的机构，并在地面完成组装，而后顶升或提升到设计构型. 机构施工法可以大幅降低脚手架和高空施工成本，且快速的顶升或提升有效地缩短了工期. 针对机构施工法的研究已经取得了一些成果. 陈志华^[3]介绍了攀达穹顶体系及其施工特点. 王小盾等^[4-5]研究攀达穹顶的几何构成和运动特点. 罗尧治等^[6-7]提出柱面网壳的“折叠展开式”提升施工方法. 张其林等^[8]将攀达穹顶的运动分析看作弹性变形和刚体位移相耦合的问题，基于有限单元平衡方程给出运动路径跟踪的迭代求解列式. 蒋本卫等^[9]将三心圆柱面网壳简化为平面连杆机构，利用动力松弛法来跟踪机构的运动形态，以系统切线刚度矩阵的最小特征值为指标，对运动过程中的形态稳定性进行判别. 蒋旭东等^[10]将攀达穹顶简化为空间铰接板片机构，推导机构的刚体运动协调方程，给出系统刚体运动路径的计算策略. 祖义祯等^[11]建立以驱动杆件伸长量为控制变量的连杆机构运动控制方程，借助弧长法进行系统运动路径的求解. 蔡建国等^[12]基于Hamilton原理，对简化的攀达穹顶模型展开过程进行数值模拟. 除此之外，机构施工法还存在网格结构的合理划分以转化为机构、顶升或提升构件的合理布置等问题，但相关的研究工作报道很少.

本文讨论“机构施工法”中面临的网格结构划分问题，提出一种几何学的划分策略. 基本思路是首先预测撤除杆件后机构的刚体位移，并将该位移强制施加到网格结构上，让网格结构从设计构型反向运动到安装面的附近，此时那些变形较大的区域需要通过撤除杆件来释放变形的区域. 数学上，将结构中所有杆件伸长量的平方和作为描述结构变形程度的指标，发现该指标可以表示为位移的二次型，并尝试用该二次型矩阵的特征向量来近似构造刚体位移. 理论上，选择那些特征值较小且与运动趋势方向夹角较小的特征向量来构造位移，可以使结构在变形较小的情况下朝安装面运动. 从设计构型开始，利用该方法不断构造近似的刚体位移方向；然后强迫网格结构按该位移方向小幅度运动并逼近安装面，这样不仅可以控制运动过程中结构的总体变形始终处于较低的水平，而且那些变形相对较大的区域逐步凸显并成为可实施划分的区域. 通过一个立体拱架和一个空间网格结构算例，验证以上划分策略的有效性.

对于一个具有N个自由度、b根杆件的网格结构，协调方程^[13]可以表示为

(1) $ {Bd} ={e}. $

式中：e和d分别为相对于构型G^t的杆件伸长（缩短）向量和节点位移向量；B为结构的协调矩阵，具体形式^[13]为

(2) $ {B} {\rm{ = }}{{B}_{\rm{L}}}{\rm{ + }}{o}({d}). $

其中B_L是协调矩阵中不含d的线性部分，o（d）为所有包含d的非线性部分. B_L可以表示为向量集合的形式：

(3) $ {{B}}_{\rm{L}}^{}{\rm{ = }}{\left[ {{{b}}_1^{\rm{T}}, \cdots ,{{b}}_k^{\rm{T}}, \cdots ,{{b}}_b^{\rm{T}}} \right]^{\rm{T}}}. $

式中：向量b_k表示杆件k对B_L的贡献，具体表示为

(4) $ \begin{split} {{b}}_k^{}{\rm{ = }} & \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0,&{\cdots}, \end{array}}&{\displaystyle\frac{{{x_m} - {x_n}}}{{{l_k}}}},&{\displaystyle\frac{{{y_m} - {y_n}}}{{{l_k}}}},&{\displaystyle\frac{{{z_m} - {z_n}}}{{{l_k}}},\cdots,} \end{array}} \right.\\ & {\left. {\begin{array}{*{20}{c}} { - \displaystyle\frac{{{x_m} - {x_n}}}{{{l_k}}}},&{ - \displaystyle\frac{{{y_m} - {y_n}}}{{{l_k}}}},&{ - \displaystyle\frac{{{z_m} - {z_n}}}{{{l_k}}}},&{\begin{array}{*{20}{c}} {\cdots},&0 \end{array}} \end{array}} \right]_{1 \times N}}{\rm{.}} \end{split} $

其中l_k（k = 1, $ \cdots$, b）为杆k的长度，m、n为杆k两端节点的编号，（x_i, y_i, z_i）分别为节点i（i = m, n）在构型G^t的坐标.

假设网格结构是几何不变的. 由式（1）可知，任何非零位移d都会使网格结构产生非零伸长e. 此外，对结构施加不同的位移，杆件的伸长一般也不同. 将所有杆件伸长的平方和L作为描述结构变形程度的指标，根据式（1）可得

(5) $ L {\rm{ = }}{{{e}}^{\rm{T}}}{{e}} = {{{d}}^{\rm{T}}}{{B}}_{}^{\rm{T}}{{Bd}} = {{{d}}^{\rm{T}}}{{Qd}}. $

式中：Q为N×N维矩阵. 可见，L是位移d的二次型. 二次型矩阵Q是关于d的非线性矩阵，通常难以直接解析. 如果考虑d足够小，则可将式（2）中包含d的非线性部分o（d）忽略，

(6) $ {{Q}} \approx {{{Q}}_{\rm{L}}} = {{B}}_{\rm{L}}^{\rm{T}}{{B}}_{\rm{L}}^{}. $

对Q_L进行特征值分解^[14]，可得

(7) $ {{{Q}}_{\rm{L}}}{\rm{ = }}{{\varPhi \varLambda }}{{{\varPhi }}^{\rm{T}}}. $

式中：Φ为Q_L的特征向量ϕ_i（i = 1, $ \cdots $, N）构成的正交矩阵；Λ为Q_L的特征值λ_i（i = 1, $ \cdots $, N）组成的对角矩阵，且特征值从大到小排列. 由于网格结构的几何不变性，由式（6）可知Q_L为正定矩阵，故λ_i（i = 1, $ \cdots $, N）均大于0. 注意到Φ为正交矩阵，ϕ_i（i = 1, $ \cdots $, N）可以作为N维空间的一组基. 于是，N维空间的任意位移向量d可以表达为如下线性组合的形式：

(8) $ {{d}}{\rm{ = }}{{\varPhi \omega }}{\rm{ = }}\sum\limits_{i = 1}^N {{\omega _i}^{}{{\phi} _i}}. $

式中：ω为ω_i（i = 1, $ \cdots $, N）组成的向量，ω_i为ϕ_i（i= 1, $ \cdots $, N）的组合系数. 将式（8）代入式（5）并考虑式（7），整理后可得

(9) $ L = {{d}^{\rm{T}}}{{{Q}}_{\rm{L}}}{d} = {{{\omega }}^{\rm{T}}}{{{\varPhi }}^{\rm{T}}}{{{Q}}_{\rm{L}}}{{\varPhi \omega }} = \sum\limits_{i = 1}^N {{\omega _i}^2{\lambda _i}}. $

在B_L不变的前提下，λ_i是确定的，因此L与ω_i（i = 1, $ \cdots $, N）的分布及相对大小有关. d的指向表达了网格结构在N维空间的运动方向，d的长度表示为网格结构在该方向运动的距离. 为了评价结构在不同运动方向的变形大小，可以统一将d转化为单位向量的形式，即

(10) $ {{d}_{\rm{n}}} = {{d}}/{{\sqrt {{{d}^{\rm{T}}}{d}} }}. $

式中：d_n仅表示运动方向，可以表示为 $ {{{d}}_{\rm{n}}}{\rm{ = }}{{\varPhi }}{{{\omega }}_{\rm{n}}}$，但 ${{\omega }}_{\rm{n}}^{\rm{T}}{{{\omega }}_{\rm{n}}}{\rm{ = }}1$. 注意到d_n为单位长度的向量，若将式（9）中的d替换为d_n，则L仅与向量ω_n中的元素分布有关，且λ_N<L<λ₁^[14]. 这表明网格结构沿不同方向运动的变形程度是不同的. 当ω_n中较大的元素恰好对应于那些较大的λ_i时，L相对较大，反之则较小. 若网格结构本身为几何可变，则存在机构位移模态^[15]，即至少λ_N = 0. 当ω_n中仅对应λ_N（即机构位移模态）取非零值时，L= 0表明体系的变形为刚体位移， ${{{d}}_{\rm{n}}}{\rm{ = }}{{\varPhi }}{{{\omega }}_{\rm{n}}} $为网格结构的刚体位移方向.

L不仅是描述网格结构变形程度的指标，而且反映了网格结构自身拓扑对特定运动的限制程度. L较小表明网格结构的拓扑对特定运动的限制作用较弱，则撤除少量杆件可以释放体系的变形而产生机构位移. 当估计网格结构划分为机构后的刚体位移时，需要使该位移对结构产生的L尽量小.

在机构施工法中，网格结构的运动路径应该是可逆的. 在进行运动分析时，一般将问题描述为从已知的设计构型朝安装面运动. 于是，对于运动路径上的任意构型，可以将该构型到安装面的方向作为运动趋势方向，数学上可以表示为

(11) $ {{{d}}_{\rm{a}}} = \left[ {\cdots,{d_{xi}},{d_{yi}},{d_{zi}}, \cdots } \right]. $

式中：d_xi、d_yi、d_zi分别为该构型中节点i到安装面的最近距离在3个坐标轴方向的分量. 若安装面由多个折面构成，则可以采用到最近折面的距离作为节点i到安装面的距离. 若安装面为曲面，则可以将该曲面用多个折面来近似拟合该曲面. 当d_a仅用于表达构型到安装面的运动趋势方向时，可以利用式（10）将d_a转化为单位向量.

式（11）定义的运动趋势向量d_a一般不是期望的机构位移方向. 若强迫网格结构按d_a运动，则结构的变形指标L较大，即需要拆除较多的杆件来释放由此产生的结构变形. 对于任意运动路径上的构型，应期望其撤除杆件后的机构位移d_r能够尽量靠近d_a. 数学上，d_r既要保证所引起的L较小，也需要与d_a的夹角尽量小.

为了能够从Φ中挑选出来用于构造d_r的特征向量集合Φ_r，可以用Ritz向量法^[16-17]进行求解. Ritz向量法是求解矩阵少部分特征向量的方法，得到的特征向量按对应的特征值从大到小顺序排列. 若指定Ritz向量法的初始向量为d_a，则求解得到的特征向量组在N维空间中都与d_a保持较小夹角. 考虑到Q_L与其逆矩阵Q_L⁻¹的特征向量相同，它们的特征值互为倒数，因此将Q_L⁻¹作为Ritz向量法的求解矩阵，求得的特征向量组总是对应Q_L较小的特征值. 在Ritz向量法求解前，需要人为指定Φ_r中包含的特征向量数目（即Φ_r的维度），但Φ_r的维度不可以任意指定. 维度过高会导致d_r引发较大的结构变形，维度过低则使d_r无法与d_a保持足够小的夹角. 采用半经验性的方法来确定Φ_r的维度. 具体思路如下：先根据设计构型Q_L特征值的分布特点，初选一个合适的阈值λ_t；然后将小于λ_t的特征值数量作为设计构型及后续所有中间构型计算d_r时的Φ_r维度. 观察最终接近安装面构型的变形集中情况并进行划分. 若划分后的刚体板块数目过多，则可以减少Φ_r的维度；若刚体板块过少，则系统刚体位移运动方向较难与d_a保持较小的夹角，以致无法运动到安装面附近，此时应加大Φ_r的维度.

在求得Φ_r后，d_r可以表示为如下线性组合的形式：

(12) $ {{{d}}_{\rm{r}}} = {{{\varPhi }}_{\rm{r}}}{{{\omega }}_{\rm{r}}}. $

式中：ω_r为对应于Φ_r中里兹向量的组合系数向量. 为了使d_r与d_a保持最小的夹角，可以令两者残差长度最小，即

(13) $ {\left\| {{{{E}}_{\rm{d}}}} \right\|_2} = ({{{d}}_{\rm{r}}} - {{{d}}_{\rm{a}}}){({{{d}}_{\rm{r}}} - {{{d}}_{\rm{a}}})^{\rm{T}}}. $

根据最小二乘法原理^[14]可知，ω_r可以通过下式确定：

(14) $ {{{\omega }}_{\rm{r}}} = {{\varPhi }}_{\rm{r}}^{\rm{T}}{{{d}}_{\rm{a}}}. $

从设计构型开始，针对运动路径上的每一个构型，采用式（12）来不断估计刚体位移方向d_r，并强迫网格结构逐渐向安装面运动. 为了反映网格结构的运动形态，可以采用下式对d_r规定一个较小的步长l_s，以确保系统小幅度地朝安装面运动：

(15) $ {{{d}}_{\rm{r}}} = \frac{{{l_{\rm{s}}}{{{\varPhi }}_{\rm{r}}}{{{\omega }}_{\rm{r}}}}}{{{{\left\| {{{{\omega }}_{\rm{r}}}} \right\|}_2}}}. $

尽管网格结构被强迫运动，但采用以上方法构造的d_r始终将结构的变形维持在较低的水平. 从后面的算例可以看出，强迫位移引起的变形主要集中在网格结构的局部区域，其他大部分区域的变形都较小. 当结构最终逼近安装面时，这些变形较大的区域非常集中，表明必须在该区域进行划分，以释放这些累积的变形.

网格结构的划分区域可以通过对比设计构型和最终靠近安装面构型的相对变形来确定. 相对变形较小的区域可以考虑为同一个刚体块，相对变形较大的区域应作为块与块间的划分界线. 变形大小的描述可以针对每个节点进行，即通过分析与该节点相连的节点在变形前、后的相对坐标差，定量判别结构的局部变形大小. 对于双层或多层网格结构，可以仅考虑同一弦层（上弦层或下弦层）节点的相对坐标变化.

在设计构型下，令节点i的坐标向量为n_i=[x_i, y_i, z_i]^T. 在同一弦层，通过杆件与节点i相连的节点数为s. 于是，s个相连节点相对于节点i的位置可以通过如下向量来描述（见图1（a））：

(16) $ {{C}_j}{\rm{ = }}{\left[ {\Delta {{n}_1},{\rm{ }} \cdots {\rm{, }}\Delta {{n}_j}{\rm{, }} \cdots ,\Delta {{n}_s}} \right]^{\rm{T}}}. $

式中： $\Delta {{n}_j}$= n_j−n_i（j=1, $ \cdots $, s），n_j=[x_j, y_j, z_j]. 对于最终靠近安装面的构型，可以对节点i及其相连节点定义一个相对坐标差的向量，记为（见图1（b））

(17) $ {C}_{j}^{'}{\rm{=}}{{\left[ \Delta {n}_{1}^{'},\cdots , \Delta {n}_{j}^{'}, \cdots ,\Delta {n}_{s}^{'} \right]}^{\rm{T}}}. $

将最终构型的节点i平移到设计构型的节点i处 （见图1（c））. 如果在2个构型中，节点i及其相连节点的相对位置未发生变化，则必然存在一个正交的坐标转换R，满足：

(18) $ {{{C}}_{j}}={RC}_{j}^{'}. $

如果在2个构型中，节点i与其相连节点的相对位置不同，则C_j'不能旋转变换到C_j，即存在

(19) $ {E} = {\rm{ }}{{C}_j} - {R}{{C}_j}^\prime \ne 0. $

式中：E为C_j与经过旋转变换后的C_j'的差. 即使C_j'和C_j不重合，也总是存在一个坐标转换矩阵R，使向量E的长度||E||₂最小，此时C_j'与C_j贴合最紧密. 此时，R可以利用最小二乘法获得，具体原理和求解步骤可以参见文献[18, 19].

以网格结构的某个弦层作为考察对象，对该层内的每个节点采用式（19）计算相应的||E||₂. 对于变形相对较小的区域，其中节点的||E||₂较低，因而可以划分为同一刚体区块. ||E||₂较高的节点区域可以认为应通过拆除杆件并设置铰接连接，以释放变形.

如图2（a）所示为一圆柱面立体拱架的设计构型. 图2（b）中，n为特征值阶数. 拱架跨度为30 m，高度为15 m，厚度为1.2 m，纵向长度为4 m. 拱架跨度方向的同层弦杆长度相同，纵向2个网格长度相等. 2个拱脚节点铰接于z=0的水平面上，且将该水平面作为安装面.

采用式（11）可以计算结构朝安装面的运动趋势方向d_a. 对设计构型的Q_L进行特征值分析发现，特征值最大值为5.86，最小值为5.71×10⁻⁴，分布如图2（b）所示. 将λ_t定为最大特征值的百分之一，即0.058 6，则可以统计小于λ_t的特征值为15，即Φ_r的维度. 从设计构型开始运动，令步长l_s=1.0 m. 针对运动路径上的每个构型，采用式（15）来构造刚度位移d_r，并强迫拱架逐步靠近安装面. 经过25步的强迫变形，拱架到达最接近于安装面的位置，如图3（a）所示. 利用式（19），定量计算上、下弦节点的局部变形指标||E||₂，用云图的形式显示于图3（b）. 图中，浅色和深色分别对应较低和较高数值. 杆件中的颜色根据该杆件两端节点的||E||₂进行线性插值后确定. 对比图3（a） 、（b）可以看出，拱架结构最终构型变形集中的区域（弯折明显区域）与||E||₂较大（深色）的区域吻合. 这说明根据||E||₂来确定网格结构的划分区域是有效的. ||E||₂指明了网格结构中哪些区域应通过撤除杆件以实现大变形的转动，于是需要进一步根据该局部区域的几何可变性特点，确定应撤除哪些杆件以保证所需的刚体转动. 对应于图3（b）中||E||₂较大的区域撤除下弦杆件，拱架转化为如图4所示的机构. 将该拱架的机构位移模态作为Φ_r，采用以上方法可以求得划分后机构距离安装面最近的最终构型，如图4所示.

若不构造d_r，直接让拱架以同样的步长l_s沿d_a运动，则能够获得一系列朝安装面运动的构型. 针对按d_r和d_a运动所计算的构型，分别计算结构的变形指标L. 如图5所示为2种情况下前25个构型的变形指标值的变化情况. 可以看出，由d_a计算得到的构型很难维持较低的结构变形，而d_r能够保证.

如图6（a）所示为一折板型双层网架的设计构型及其安装面. 网架外形由5个平面组成，安装面由9个平面构成. 网架周圈的下弦节点三向铰接于安装面的四边. 网架及安装面的几何尺寸见图6（b）、（c）.

采用式（11）计算网架指向安装面的运动趋势方向d_a. 对设计构型的Q_L进行特征值分析，最终确定Φ_r维度为16. 从设计构型开始，取2个构型间的运动步长为30.0 m，不断迫使网架按式（15）构造的位移d_r向安装面运动. 经过27步的强迫变形，网架到达最接近于安装面的位置，如图7所示. 利用式（19），进一步定量计算上弦节点的局部变形指标||E||₂，采用云图的形式，如图8所示.

对比图7、8可以看出，侧面4个平面内均形成了2个弯折变形明显的区域. 除此之外，侧面4条脊线处反映出大变形. 上述区域均应通过撤除杆件转化为能够进行大变形转动的铰接区域，如图9的虚线所示. 图中，箭头指示板件的运动趋势. 脊线与其两侧平面内横向的铰接成相交关系. 如果虚线处仅形成以该线为轴转动的铰，由于脊线处节点的运动趋势无法协调（如图9的板件1、2所示），相邻2个侧面的板件在几何上不具备刚体运动的能力. 为了有效形成机构，将脊线处的杆件完全撤除以使侧面的刚体板件不再相邻. 网架划分后的最终形式如图10（a）所示. 将撤除杆件后的网架的机构位移模态作为Φ_r，采用同样的方法可以求得距离安装面最近的最终构型，如图10（b）~（d）所示.

本文将网格结构中所有杆件伸长量的平方和作为描述结构变形程度的指标，并将其表示成为位移的二次型. 利用Ritz向量法来选择该二次型矩阵的特征向量，以近似构造与运动趋势方向最接近的刚体位移. 进一步强迫网格结构按所构造的刚体位移运动至安装面附近，那些变形集中的区域应通过划分来释放变形的区域. 2个算例的计算结果反映了提出划分策略的有效性. 杆件伸长量的平方和是一个宏观的指标，通过将该指标的数值维持在较小的水平上，根据该方法所构造的近似刚体位移能够使网格结构中的大多数区域变形保持在很低的水平，这是该划分策略有效的重要原因. 此外，本文方法只是获得了网格结构按估计的刚体位移强迫运动后变形集中的区域. 如何撤除杆件以实现该区域的大转动是主观能动的. 譬如算例1中，可以撤除变形集中区域的上弦杆件来实现结构在该区域的大角度转动. 网格结构的最终划分除决定于安装面的形式外，还决定于Φ_r的维度. 随着Φ_r维度的增大，刚性板块数量会增多，划分后的网格结构会更贴合于安装面，但在工程中增加了机构运动控制的复杂性. 当在实际工程中采用本文的划分方法时，Φ_r的维度往往应进行多次试算和比较来合理确定.