一种联合映射法和Delaunay法的拟桁架自由曲面网格划分方法与流程

本发明涉及空间结构网格划分方法，尤其涉及一种自由曲面网格结构的

网格划分方法。

背景技术：

建筑网格是建筑设计中常用的一种方法。建筑网格设计法以其有条理、有秩序的组织结构和有规律、有节奏的韵律变化，深受广大建筑师的喜爱。近些年来，我国出现了较多的大跨度空间网壳结构，均把结构网格外露形成建筑网格，如阳光谷、深圳湾体育中心、辰山植物园等等。如何使得外露的结构网格，在保证网格规则性、杆件流畅性、杆长一致性等条件下，更好的体现建筑师赋予建筑的灵魂，是当今结构师和建筑师不得不面对的一个重要问题。

网格划分起源于上世纪60年代，最初应用于有限元领域，在这个过程中，波前法、Delaunay法、映射法等网格划分方法相继被提出并发展，但很多时候，现有的网格划分方法所生成的网格是非结构化网格，并不能满足建筑美观的要求，即网格大小均匀、杆件线条流畅等美学指标。近年来，有不少学者在这一领域进行了研究探索，并取得了一定的成果，提出了基于主应力轨迹线的波前推进法生成曲面网格；基于映射思想提出了自定义单元法来生成网格；采用等参线分割法和改进的Delaunay点云网格生成法进行曲面网格划分；基于曲面展开和Delaunay优化来生成网格。

技术实现要素：

本发明的目的是克服现有技术的不足，提供一种自由曲面网格结构的网格划分方法。

联合映射法和Delaunay法的拟桁架自由曲面网格划分方法，具体步骤为：首先将NURBS曲面从空间域映射到平面域上，得到相应的平面网格划分区域A；其次在区域A的包络矩形内生成均匀的点阵{P₀}，求{P₀}到区域A边界的距离函数d(x,y)，判断点阵分布的{P₀}中各点处在区域A的内部还是外部，剔除在区域A外的点，得到在区域A内的点阵{P_in}，剔除点阵{P_in}中部分点，得到与期望的网格密度基本一致的点阵分布{P}；接着将点阵{P}用Delaunay法连成三角形网格T0，剔除区域A外的三角形，得到网格T；将网格T作为平面桁架结构，用调节函数h(x,y)代表点阵{P}中两点连线的中点(x,y)处的网格杆长，网格T的实际杆长与期望杆长之差，作为平面桁架的不平衡力，在不平衡力作用下移动各节点；进行循环迭代求解桁架结构的平衡位置，直到各点受到的最大不平衡力不大于收敛控制值；最后将平面网格映射回空间曲面，生成自由曲面网格。

上述方法中，各步骤可采用如下方式实现：

通过调节函数h(x,y)，控制不同位置的网格疏密，调节函数设为：

h(x,y)＝(1-1/k)/(g×h₀)×|d(x,y)|+1/k

其中|d(x,y)|表示点(x,y)到边界线的最小距离，h₀为期望网格杆长，k为期望网格杆长控制加密系数，k＝1表示不加密；g为期望网格杆长控制大小的过渡区间。

另外，可通过旋转初始点阵{P₀}，控制网格的整体走向；旋转前后的点分别为p(x,y)、p(x’,y’)，以点(x₀,y₀)为圆心旋转α度，按下列公式得到p(x’,y’)：

x’＝(x-x₀)×cos(a)-(y-y₀)×sin(a)+x₀

y’＝(x-x₀)×sin(a)+(y-y₀)×cos(a)+y₀

通过距离函数d(x,y)的正负性，判断点阵{P₀}处在区域A的内部还是外部点的方法为：先将曲面的边界线离散为小线段；再将边界线映射到平面后，得到由一个或多个多边形表示的区域A；然后采用射线法判断点在单个多边形的内外，即位于区域A内时，d(x,y)小于零，位于区域A外时，d(x,y)大于零。

对区域A内的点阵{P_in}进行调整的方法为：先求出{P_in}中各点的调节函数h(x,y)，再根据概率论中的蒙特·卡罗原则，按正比于1/h(x,y)²的概率，随机删除部分点，得到与期望网格密度基本一致的点阵分布{P}。

所述的区域A外的三角形为重心处于区域A外的三角形。

本发明与背景技术相比具有的有益效果是：

1)本发明能适应一般的自由曲面，生成的网格大小、网格走向可控、线条流畅，克服了目前网格划分方法的不足，既满足了建筑美学的审美要求，又达到了传力连续的效果。

2)本发明能够实现程序化，可快速实现网格的自动生成，大大提高了网格划分的质量和速度。

附图说明

图1是联合映射法和Delaunay法的拟桁架自由曲面网格划分过程示意图；

其中：

a)是待划分NURBS曲面；

b)是NURBS曲面在参数域上的边界线；

c)是参数域上的初始布点；

d)是参数域上的剔除边界之外的点；

e)是参数域上的筛选后的点

f)是用Delaunay法在参数域上生成三角形初始网格；

g)是在参数域上拟桁架法迭代41次时的网格；

h)映射回NURBS曲面的网格；

图2是α＝0°，h0＝4.0，k＝1，g＝0时的曲面网格划分结果；

图3是α＝30°，h0＝4.0，k＝1，g＝0时的曲面网格划分结果；

图4是α＝0°，h0＝2.0，k＝2，g＝10时的曲面网格划分结果；

图5是α＝0°，h0＝2.0，k＝3，g＝10时的曲面网格划分结果。

具体实施方式

下面结合附图和一个具体曲面作为实例对本发明做进一步阐述。

某NURBS曲面如图1.a)所示，利用联合映射法和Delaunay法的拟桁架自由曲面网格划分方法对其进行网格划分。首先将该NURBS曲面从空间域映射到平面域上，得到相应的平面网格划分区域A，如图1.b)所示；其次在区域A的包络矩形内生成均匀等间距分布的点阵{P₀}，如图1.c)所示。求{P₀}到区域A边界的距离函数d(x,y)，判断点阵分布于{P₀}中各点处在区域A的内部还是外部。此处的判断是通过通过距离函数d(x,y)的正负性判断来实现的，具体方法为：先将曲面的边界线离散为小线段；再将边界线映射到平面后，得到由一个或多个多边形表示的区域A；然后采用射线法判断点在单个多边形的内外，即位于区域A内时，d(x,y)小于零，位于区域A外时，d(x,y)大于零。

对每个点经过上述判断后，剔除在区域A外的点，得到在区域A内的点阵{P_in}，如图1.d)所示。剔除点阵{P_in}中部分点，得到与期望的网格密度基本一致的点阵分布{P}。对区域A内的点阵{P_in}进行调整的一种具体方法为：先求出{P_in}中各点的调节函数h(x,y)，再根据概率论中的蒙特·卡罗原则，按正比于1/h(x,y)²的概率，随机删除部分点，得到与期望网格密度基本一致的点阵分布{P}，如图1.e)所示。此处，h(x,y)越大，表示相应位置的期望网格密度越小，所需的节点越少。

接着将点阵{P}用Delaunay法连成三角形网格T0，如图1.f)所示。剔除区域A外的三角形(区域A外的三角形为重心处于区域A外的三角形)，得到网格T；将网格T视为平面桁架结构，用调节函数h(x,y)代表点阵{P}中两点连线的中点(x,y)处的网格杆长，网格T的实际杆长与期望杆长之差，作为平面桁架的不平衡力，在不平衡力作用下移动各节点。调节函数h(x,y)代表了某一网络杆长，可用于控制不同位置的网格疏密，调节函数设为：

h(x,y)＝(1-1/k)/(g×h₀)×|d(x,y)|+1/k

其中|d(x,y)|表示点(x,y)到边界线的最小距离，h₀为期望网格杆长，k为期望网格杆长控制加密系数，k＝1表示不加密；g为期望网格杆长控制大小的过渡区间(在划分网格之前，预先对各区域三角形的网格杆长设定一个期望值，过渡区间就是指两个具有不同期望值的三角网格的重心间的曲面距离)。

对上述过程进行循环迭代，求解桁架结构的平衡位置，即直到各点受到的最大不平衡力均不大于预先设置的收敛控制值。例如，本实例中通过循环迭代41次后达到桁架结构的平衡位置，其结果如图1.g)所示。最后将平面网格映射回空间曲面，生成自由曲面网格。最终映射回NURBS曲面的网格如图1.h)所示。

另外，针对不同的自由曲面，可根据其曲面的特征，通过旋转初始点阵{P₀}，控制网格的整体走向；旋转前后的点分别为p(x,y)、p(x’,y’)，以点(x₀,y₀)为圆心逆时针旋转α度，按下列公式得到p(x’,y’)：

x’＝(x-x₀)×cos(a)-(y-y₀)×sin(a)+x₀

y’＝(x-x₀)×sin(a)+(y-y₀)×cos(a)+y₀

下面将上述方法应用到两个具体算例中，以便更好地展示本发明的效果。两个算例中，由于中间过程的实施步骤与上一实例基本相同，因此不再展示中间过程的各附图。两个算例中具体参数和结果如下：

算例1，某NURBS曲面，高约15.6(均采用无量纲，下同)，在参数平面内的包络矩形的左下角为(-11.0,-11.7)，右上角为(68.5,64.6)，设定不同的角度α＝0°及α＝30°，其余参数均为：h₀＝4，k＝1，g＝0，点阵的旋转中心取包络矩形的中心，最大不平衡力的收敛控制值取f＝0.1kh₀，进行网格划分。在α＝0°及α＝30°的旋转角度下，得到如图2、图3的网格划分结果。两者的网格呈现明显不同的整体走向，说明通过本方法通过调整旋转角度，可调控整体网格走向，以符合不同的曲面特征。

算例2，某NURBS曲面，由算例1中的曲面经裁剪后得到，选取内边界线作为加密线，最大不平衡力的收敛控制值也取f＝0.1kh₀，分别采用α＝0°，h₀＝2.0，k＝2，g＝10及α＝0°，h₀＝2.0，k＝3，g＝10的设定参数，进行网格划分。k＝2和k＝3的参数下，得到如图4、图5的网格划分结果。两者的网格尺寸呈现明显的差异，说明本方法可控制不同区域网格尺寸。