随着社会的进步和计算机的快速发展，建筑的造型越来越多样化。自由曲面因形式多样、线条柔美、造型流畅，越来越多地出现在了现代建筑中。但是自由曲面不同于传统的标准曲面，无法用解析函数来描述，传统标准曲面上的网格划分方法无法直接应用于自由曲面上，在自由曲面上布置杆件长度一致、分布均匀、线条流畅的单层网格是一个困难且复杂的问题。

　　 实际自由曲面工程的网格划分大都采用人工经验或者半经验的方法，如米兰新贸易会展中心^[1]的曲面网格采用人工划分的方法，网格质量严重依赖于设计师的水平;上海世博会阳光谷^[2]的曲面网格是采用在NURBS参数平面人工划分网格，然后将其映射至曲面，得到的网格较为扭曲，需要后期对网格进行优化。

　　 一般而言，高质量的自由曲面网格很难一步到位，因此需对网格进行优化调整，以改善网格的均匀性、流畅性和美观性。文献[2]中采用索网有限元的方法对阳光谷网格进行优化调整，得到了最终的结构网格;文献[3]采用改进的Laplacian光顺法对网格进行优化调整，改善了优化调整过程中网格的收缩情况。本文从几何形态上出发，提出了基于杆长均匀性的网格优化方法，并基于Rhino+Grasshopper平台编制了网格均匀性优化程序，该程序以网格杆长标准差最小为优化目标，通过调整网格节点在曲面上的位置，最终实现对自由曲面网格的优化。

　　 NURBS(Non-Uniform Rational B-Spline)曲线曲面，即非均匀有理B样条曲线曲面，是一种通用的描述自由曲线曲面的方法，通过控制点、权因子、基函数三个因素就可精确表示任意自由曲线曲面。NURBS为描述任意形式的曲线或曲面提供了统一的数学表达式，文献[4]对NURBS曲线曲面有较为系统的介绍。

　　 一条k次NURBS曲线采用有理分式表达形式如下:

　　 式中:u为曲线上的参数;P(u)为参数u在曲线上对应的三维点阵;ω_i(i=0,1，…，n)为权因子;d_i(i=0,1，…，n)为控制顶点，顺序连接成控制多边形;N_i,k(u)为k次规范B样条基函数，由节点矢量按德布尔-考克斯递推公式决定。

　　 类似于NURBS曲线，一张k×l次NURBS曲面的有理分式表达形式如下:

　　 式中:u,v代表曲面两个方向的参数;P(u,v)为曲面上参数u,v对应的三维点阵;ω_ij(i=0,1，…，m;j=0,1，…，n)为权因子;d_ij(i=0,1，…，m;j=0,1，…，n)为控制点网格;N_i,k(u)(i=0,1，…，m)和N_j,l(v)(j=0,1，…，n)分别为u向k次和v向l次的规范B样条基函数，分别由u向和v向节点矢量按德布尔递推公式决定。

　　 自由曲面网格质量主要取决于单元形状的好坏、杆件长度的均匀性和网格线条的流畅性，网格线条流畅性难以定量评价，故在自由曲面网格优化研究中，常将网格的单元形状指标和杆件长度指标作为网格质量评价指标^[5]。

　　 对于网格形状指标而言，网格形状越接近正多边形(如正三角形或正四边形)，则网格质量越好;对于杆件长度指标而言，杆件长度差异越小，则网格均匀性越好^[5]。其中网格单元形状指标和杆件长度指标并不是完全独立的，尤其是三角形网格，形状质量越好，越接近正三角形，其杆件的长度也更为均匀。杆件是构成网格的基本单元，考虑到经济性和美观性，杆件长度尽量均匀、规格尺寸尽量统一是结构优化的一个目标。

　　 本文以杆件长度指标为参考指标，对自由曲面网格进行优化研究。优化程序采用统计学的方式，以杆件长度为样本，利用杆长平均值l^—和杆长标准差δ_l来计算杆件长度指标。

　　 式中:l_i为网格各杆件长度;N为杆件总数。

　　 当δ_l值越小，则表示杆件长度越均匀，即网格越均匀。

　　 优化程序采用在物理空间直接优化的方法，以网格杆长标准差最小为优化目标，通过调整网格节点在曲面上的位置，实现对空间网格的均匀性优化，优化流程图如图1所示。优化时将网格节点分为固定节点、边界非固定节点和内部非固定节点三类，针对不同的网格节点类型采取相应的处理措施。

　　 对于开放型曲面，为使网格能更加贴合曲面轮廓，当网格节点位于曲面的角点时，将其定义为固定节点。对于其他不想在优化过程中发生位置改变的网格节点可通过人为输入的方式将其定义为固定节点。固定节点在优化过程中维持坐标不变。

　　 曲面角点是曲面的特征点，可通过下述方法确定:先提取曲面边界线，然后提取曲面边界线的端点作为曲面角点。

　　 当网格节点与曲面边界线的最短距离小于某一容许误差时可视该节点为边界节点，若边界节点为非固定节点，这些节点将沿该边界线进行调整以维持网格轮廓贴合曲面轮廓。边界线上节点的最优布置为等弦长布置，程序中采用中垂面的方法进行调整。

　　 图2所示曲线为目标曲面的边界线，P₁,P₂,P₃，…，P_n为该曲线上的网格节点，其中P₁和P_n为曲面的角点，也即优化过程的固定点，其他节点采用中垂面的方法进行调整，直到连接这些节点的杆件长度相等。以P₂节点的调整为例介绍该类节点的调整步骤。

　　 Step1:维持点P₁和P₃不动，计算出P₁和P₃的中点P_1,3，根据P_1,3和P₃坐标计算由点P_1,3指向P₃的向量，根据点P_1,3和向量计算出以为法向量，以点P_1,3为原点的参考平面。

　　 Step2:计算参考平面与边界曲线的交点P'_1,3。若参考平面与曲线只有一个交点，则该点为P'_1,3;若有多个交点，以位于节点P₁和P₃之间的交点为P'_1,3。

　　 Step3:令P₂=P'_1,3，即完成了点P₂的调整。

　　 Step4:重复步骤Step1～Step3调整其他节点。

　　 按照Step1～Step4对边界节点进行多次循环调整，并计算每次调整后节点P₁,P₂,P₃，…，P_n间连接杆件长度的标准差δ_i，并判断(ε为精度设置值)，若满足要求，则跳出循环;否则继续下一循环。

　　 网格节点大部分为曲面内部节点，对于曲面内部节点，采取逐点优化的方式进行优化调整。网格内部节点周围拓扑连接的节点往往有多个，对该节点的优化调整过程也就是在目标曲面上寻找目标点使得该点与周围各节点连接杆件长度渐趋相等的过程，即各杆件长度的标准差逐渐减小的过程。典型的内部节点优化单元如图3所示。

　　 图3所示优化单元中曲面内部节点P₀分别与节点P₁,P₂，…，P₆相连，优化过程中固定P₁,P₂，…，P₆位置不动，调整P₀在目标曲面上的位置，使得点P₀与周围各点的连接杆件长度标准差最小。当点P₁,P₂，…，P₆共球面时，最优P₀点与周围各点距离相等，此时杆长标准差等于0。曲面内部点采用网格法^[6]进行逐点调整，具体方法如下:

　　 Step1:计算节点P₁,P₂，…，P₆在目标曲面上所对应的曲面参数u,v值。

　　 Step2:选取上述u,v值中的最大值和最小值u_max,u_min,v_max,v_min，创建曲面搜索域Ω={(u,v)|u_min<u<u_max,v_min<v<v_max}。对于某方向闭合的曲面，如u向闭合的柱形曲面，若优化单元横跨曲面上u=0曲线，这时需将搜索域进行调整，具体调整方式为先将位于曲面终止侧的点(u_s,max/2<u<u_s,max点)的u值更新为u-u_s,max(u_s,max为曲面u向参数最大值)，然后再选取u,v值中的最大值和最小值并计算曲面搜索域Ω。

　　 根据曲面参数u_i,v_j和曲面参数方程可计算得到节点P_i,j在整体坐标系下的坐标。程序中采用Rhino中的Surface.PointAt()函数计算P_i,j坐标。

　　 Step4:计算点P_i,j与点P₁,P₂，…，P₆距离的标准差δ_i,j，然后计算同一个小格的四个角点标准差的平均值

　　 Step5:取上述计算的标准差平均值最小的小格作为新的搜索区域Ω，取该小格四个角点坐标的平均值得到点P_opt，采用函数Surface.ClosestPoint()在目标曲面上找到P_opt的最近点P'_opt，点P'_opt即为最优点，调整点P₀至该位置。

　　 Step6:重复步骤Step3～Step5直至点P₀与周围各点的杆长标准差的相对改变量小于精度要求。

　　 Step7:重复步骤Step1～Step6完成其余内部非固定节点的优化。

　　 对所有节点进行多次迭代优化，直至网格杆件均匀性满足要求。

　　 为验证本文网格优化方法的正确性和适用性，本节以三个实例进行说明。

　　 如图5所示，边长为9m×9m的平面四边形被划分成9个四边形网格，对该初始网格进行优化，最终得到图6中边长为3m×3m的9个四边形网格。由于本例有数值解，因此较易检验本优化方法的正确性。本文优化方法不改变网格的拓扑结构，因而只要初始网格的拓扑结构不变，不管其形状如何，最终都能优化得到图6所示的形状。

　　 采用等参线网格划分^[7]方法对图7所示曲面进行网格划分，得到图8所示三角形网格，该网格由540个三角形单元组成。等参线网格划分方法基于曲面的结构线进行划分，因而划分的网格质量受曲面结构线影响较大，曲面结构线的不均匀会直接导致网格质量较差。从图8可以看出网格杆长差异较大、均匀性较差，且网格较为扭曲、流畅性较差。

　　 将初始网格和目标曲面输入优化程序，优化最高迭代步设置为50步，对该网格进行优化调整。优化调整后的网格如图9所示，可以看出，优化后的网格均匀性和流畅性均得到了明显的改善。网格杆长最大值由优化前的4.40m减小为优化后的3.03m，杆长最小值由优化前的0.88m增大为1.76m，杆长平均值由优化前的2.50m变为优化后的2.46m。优化过程中网格杆长标准差的变化如图10所示。

　　 采用映射法^[2]对图11所示曲面布置网格，映射法基本思路是先在曲面的二维参数域上进行网格划分(图12)，然后将二维网格映射至曲面得到空间网格(图13)。该方法是在平面上进行网格划分，因此操作简单且易于构造网格拓扑形式，但由于是非线性映射，得到的空间网格会扭曲变形，网格质量较差，故需对其进行优化调整。图13中网格由1 323个三角形单元组成，从图中可以看出，网格杆长差异较大，均匀性较差。

　　 将初始网格和曲面输入优化程序，优化最高迭代步设置为30步，对该网格进行优化调整。优化调整后的网格如图14所示，可以看出，优化后的网格均匀性和流畅性均得到了明显的改善。网格杆长最大值由优化前的5.78m减小为优化后的3.87m，杆长最小值由优化前的1.44m增大为1.54m，杆长平均值由优化前的2.67m变为优化后的2.60m。优化过程中网格杆长标准差的变化如图15所示。

　　 本文提出了基于杆长均匀性的网格优化方法，以平面四边形网格优化算例验证了本方法的正确性，以两个自由曲面网格优化算例证明了本方法的适用性。结合这三个算例得出以下结论:

　　 (1)通过在Rhino+Grasshopper平台上编制网格均匀性优化程序，实现网格划分和网格优化在同一平台上相继进行，使网格设计与优化调整的衔接更加紧密、操作更加方便，同时也极大地提高了工作效率。

　　 (2)基于杆长均匀性的网格优化方法是对空间网格的直接优化，优化调整过程保证网格节点始终位于曲面上，且能智能识别曲面角点和边界点，避免了网格形状失真的情况出现，也保证了网格在整个优化过程中始终与曲面贴合，实现了结构网格与建筑曲面的完美融合。