基于T-Spline 的全自动几何拓扑修复方法

2019-09-15池宝涛张见明鞠传明

自动化学报 2019年8期

池宝涛张见明鞠传明

自由曲线曲面造型技术[1]是计算机辅助几何设计(Computer aided geometric design,CAGD)的核心基础.非均匀有理B 样条(Non-uniform rational B-Splines,NURBS)[2]以其优良特性广泛应用于自由曲线曲面造型,而成为计算机辅助设计和制造(Computer aided design and manufacturing,CAD/CAM)的一个标准.目前,通用的商业计算机辅助工程(Computer aided engineering,CAE)软件计算仿真流程可分为前处理、中间求解计算、后处理三大模块.其中数值模拟的前处理模块主要用于导入专业计算机辅助设计(Computer aided design,CAD)软件构建的几何模型的标准数据文件或对已有的CAE 几何模型进行简单操作或修改、网格生成、施加边界条件及定义物理属性和求解参数等.复杂问题数值模拟难以实现自动化的主要性能瓶颈在于前处理涉及大量的人工交互,且其效率高低严重依赖于用户经验和知识水平.根据美国Sandia 国家实验室Michael Hardwick 和Robert Clay 提供的数据,如图1 所示,在数值模拟过程中,前处理模块的用时占据了整个数值模拟过程用时的绝大部分,其中对于CAD 模型的处理用时占60%(图1 中步骤1∼3),网格生成用时占20%(图1 中步骤4 和步骤5),而真正用于数值计算的时间仅占整个数值模拟过程的4% 左右(图1 中步骤8),因此,实现前处理的自动化是提高全自动CAE 分析效率和精度的关键.然而,由于高质量的网格生成及高效可靠的CAE 分析对几何模型通常有特殊的要求,而初始输入的几何模型由于其来源的多样化,在进行CAD几何模型设计时没有考虑或很难预见到这些特殊的要求,因此在导入到CAE 系统后的几何模型通常会引入几何“噪声”.

图1 美国Sandia 国家实验室数值模拟过程用时数据统计Fig.1 Statistics related to numerical computation from Sandia National Laboratories of USA

CAD 模型的几何“噪声”种类难以尽数,且其来源及其产生原因多种多样,主要可归纳为两类:一是CAD 模型在几何造型设计制造时存在缺陷.在传统的几何造型设计及产品研发过程中,几何形状设计与物理分析方法完全分属不同的工程领域,设计者主要关注CAD 模型的构建,而忽视不理想的几何特征对后续物理分析的影响.二是不同系统之间的数据传递造成的CAD 模型几何数据及拓扑信息缺失.在常见的传统商业软件中,CAD 系统与CAE系统两者相互独立,两者之间无法进行交互.在进行产品设计或性能分析时,需要将CAD 模型导入到CAE 系统中,在此数据传递过程中会造成原始几何数据及拓扑信息的丢失.另外,在传统的工程分析计算中,CAE 分析模型的模型描述往往对原始CAD几何模型进行修改或简化,无法体现几何设计模型中包含的复杂且精确的几何信息,CAE 分析模型与CAD 几何模型不统一,且两者之间存在巨大差异,导致CAE 分析自动化程度低.

CAD 模型几何拓扑修复[3]主要是指在保证原有CAD 模型不失真的前提下,对CAD 模型中存在的几何瑕疵(如短边、窄面、退化边、退化面、非连续光滑边界及尖锐特征等)进行处理,将存在的“错误”的“脏”几何转换为可满足CAE 分析的“干净”几何.目前,基于不同的应用需求,国内外学者对于CAD 模型修复方法的研究主要分为两类:一类基于连续曲面;另一类基于离散曲面.现在主流的CAD软件通常基于连续曲面造型,主要依赖NURBS 的建模功能实现复杂模型的建模,几何精度高,但受拓扑限制且存储数据量大.文献[4]提出了一种基于连续曲面的方法，利用少许连续曲面片进行拟合处理曲面间的非连续缺陷.该方法具有局限性,且在修复过程中存在不成功情况.文献[5]提出了一种基于连续曲面的方法,采用“实操作”消除细小几何特征,操作复杂,拟合数据量大,且在力学性能分析中,局部细小特征恰恰为应力集中的地方,是研究的关键对象.文献[6]提出了一种基于连续曲面的方法利用“实操作”修复实体模型,对CAD 模型中的几何瑕疵进行增删操作,利用“实操作”后的连续曲面替代原实体模型中的曲面,该方法存在以下缺陷:1)完成修复的主要步骤需要大量的人工交互,对于复杂实体模型的几何修复效率低且错误率高;2)修复后的模型与输入原模型几何数据及拓扑信息存在巨大差异,改变了原始CAD 模型的真实信息.基于离散曲面进行拓扑修复方法主要通过大量的三角形面片进行曲面表征,相关几何计算是线性的、效率高,但存储数据量大、几何精度差,无法满足高精度的数值模拟要求.文献[7]对离散模型的修复方法研究的最新进展进行了系统的归纳和总结,并概括总结了典型的离散模型错误,包括法向不一致、孤立点和边、非二边流形、退化单元、孔洞、细缝和单元重叠、相交、拓扑噪音、数据噪音、锐利几何特征丢失等.该方法分类复杂,且多数错误为在离散模型中存在,或在将连续参数曲面模型转化为离散模型的过程中引入的错误,而在基于参数曲面设计的CAD 模型中并不存在.文献[8]提出了一种采用离散曲面的方法,利用大量的三角形面片进行模型处理,该方法难以满足高精度的数值模拟要求且存储数据量大.此外,离散曲面表征仅涉及离散面片间等低层拓扑,无法直接用于需高层拓扑支持的操作,如若需要,需重构连续曲面模型中的边界表征数据(Boundary representation,B-Rep).参数曲面造型在数值模拟中应用最为广泛,如在等几何分析及计算流体力学等领域,模型的几何精度对保证数值模拟精度非常关键,因此高精度的数值模拟通常要求采用参数曲面造型方法,甚至要求采用高阶参数曲面.本文从高精度数值模拟的角度出发,采用第一类基于连续曲面方法进行全自动几何拓扑修复.

基于本文提出的一种基于T-Spline 的全自动几何拓扑修复方法,开发了一套CAD/CAE 集成系统,统一了几何模型与计算分析模型,实现了CAE与CAD 两者的无缝集成,所有拓扑修复操作及后续CAE 分析计算均在同一环境下进行,避免了几何模型在CAE 与CAD 系统间转换时造成的数据丢失.基于边界面法[9]的完整实体CAE 分析软件[10]是将边界积分方程与计算机图形学相结合,直接利用CAD 实体模型中的边界表征数据实现复杂结构CAE 分析自动化,将几何实体模型的曲面片作为物理量插值的基本单位,并且用于CAE 分析的网格单元可从原始输入CAD 模型中计算用于物理量插值和数值积分所需的精确几何信息,如切向量、法向量、雅克比等,从而避免了几何误差的引入并保证了计算结果的高精度.与主流商业CAE 软件不同的是,该软件是直接基于UG 平台,利用面向对象的程序设计语言Visual C++二次开发的,设计了用于读取原始输入CAD 模型B-Rep 数据的统一几何数据接口模块,降低了进程间频繁交互的通讯频度.该系统采用本文方法设计了全自动几何拓扑修复模块,且全自动几何拓扑修复与后续CAE 分析在同一环境下进行,避免了不同系统之间的数据传递造成的CAD 模型几何数据及拓扑信息缺失,统一了几何模型与分析模型,实现CAE 与CAD 两者的无缝集成.完整实体CAE 分析软件前处理模块框架如图2所示.

图2 完整实体CAE 分析软件前处理模块框架Fig.2 Preprocessing module of center for complete solid analysis software for engineering structures

完整实体CAE 分析软件的总体设计原则为:

1)软件既融于UG 又独立于UG.为了提高软件的独立性与可移植性,本软件的实现仅在几何造型方面依赖于UG 而其他子模块又尽可能地独立于UG,如用户界面模块、数据管理模块、计算仿真模块、后处理显示模块等.这种设计可使软件仅需对几何接口模块做部分修改就能移植到其他CAD 软件中.

2)减少与CAD 软件接口函数的交互.由于CAD 软件一般都提供二次开发的接口函数,因而可以直接获取到CAD 模型的几何数据,然而这些接口函数在调用时相当于一个黑匣子,频繁的交互会增加进程间的通讯频度,使得数据管理难度增加,影响计算效率.因此可以通过读取原始CAD 模型的B-Rep 数据来构建用于CAE 分析的B-Rep 数据,并且对参数曲线与参数曲面进行自定义解析,从而实现与CAD 系统的解耦.

3)采用微软基础库(Microsoft foundation classes,MFC)设计软件的用户交互界面,避免使用CAD 软件提供的用户界面设计工具.因为CAD软件提供的用户界面设计工具提供的界面设计功能没有MFC 的丰富全面且不利于软件移植到其他CAD 软件中.

4)软件中的各个模块应尽量相对独立,且每个模块要能在不修改源代码的情况下实现扩展.为各模块和子模块设计C++类并进行独立封装,使用继承机制与合成机制实现各模块类的功能扩展.如果一个类需要调用另一个类可以通过设计第三方类进行类与类之间的数据传递.完整实体CAE 分析软件整体框架如图3 所示.

本文提出了一种基于T-Spline 的全自动几何拓扑修复方法,与传统基于连续曲面的拓扑修复方法相比,基于T-Spline 比基于NURBS 进行拓扑修复最大优点在于允许T 型节点的存在,在保证同等几何精度的前提下,减少不必要曲面数据存储,而且不受拓扑限制.与传统基于离散曲面的拓扑修复方法相比,基于T-Spline 拓扑修复的模型表征几何精度高,可满足高精度的数值模拟,且存储数据量大幅降低.本文方法可对具有几何“噪声”和细小特征的复杂结构进行全自动拓扑修复操作,成功解决了CAD 几何模型中存在的短边、窄面、退化边、退化面、非连续光滑边界及尖锐特征等常见的不必要几何特征,利用新生成的“虚边”、“虚面”处理几何“噪声”,同时重构CAD 几何模型的B-Rep,在一定程度上可以明显改善网格生成的质量,降低网格的生成规模,进而有利于减小计算规模,缩短计算时间,提高CAE 分析计算的精确性和可靠性.分析软件界面如图4 所示.

1 NURBS 和T-Spline

样条理论和方法的建立为自由曲线曲面造型研究提供了重要的理论依据和技术工具.NURBS 以其优异特性而广泛应用于CAD 造型和CAE 分析中,如非负性、单位分解性、凸包性、紧支性、灵活的连续性以及保型的p 细分(升阶)或h 细分(插入节点)等.但是NURBS 曲面受到张量积的拓扑限制,使得NURBS 曲面在造型中控制点冗余.此外,在创建复杂几何模型时,裁剪曲面往往是造型中一个很重要的操作,在裁剪区域曲面交界处不可避免会产生缝隙,从而破坏几何模型曲面与曲面直接的连续性.为了克服NURBS 的固有缺陷,文献[11]提出一种新的样条T-Spline,即T-Spline 控制网格中允许T 型控制点存在的样条曲线曲面.T-Spline是在NURBS 的基础上推广和发展的,它突破了传统的B 样条曲面控制网格必须满足的拓扑限制的局限,可以有效消除冗余的控制点,而且在曲面拼接、曲面加细、曲面简化等方面具有独到之处.在一维空间,T-Spline 曲线与NURBS 曲线基本是一致的.在二维和三维空间,T-Spline 的控制网格中允许存在T 型控制点,且在保证几何造型准确度的前提下,可有效减少不必要的控制点,减少曲面解析存储的数据量.二者的区别如图5 所示.

图3 完整实体CAE 分析软件整体框架Fig.3 The software framework of center for complete solid analysis software for engineering structures

图4 完整实体CAE 分析软件界面Fig.4 The software interface of center for complete solid analysis software for engineering structures

图5 NURBS 与T-Spline 区别Fig.5 The difference between NURBS and T-Spline

一条p次NURBS 或T-Spline 曲线定义为

其中,Pi是控制点,ωi是权因子,Ni,p(u)是定义在非均匀节点矢量U上的p次样条基函数.

基于De Boor-Cox 递归公式,第i个p次样条基函数Ni,p(u)定义如下:

一张在u方向p次、v方向q次的NURBS 或T-Spline 曲面定义为

其中,Pi,j是控制点,ωi,j是权因子,Ni,p(u)和Nj,q(v)分别是定义在非均匀节点矢量U和V上的非有理样条基函数Ni,p(u)和Nj,q(v),定义同上.

2 非理想几何特征分类、识别及拓扑修复

全自动几何拓扑修复是CAE 分析前处理过程中一个非常复杂的关键问题.本文将几何模型的BRep 中存在的短边、窄面、退化边、退化面、非连续光滑边界及细小曲面、尖锐特征等常见几何噪声统称为非理想几何特征.图6 和图7 显示了钢架焊缝和汽车桥壳模型的几何噪声.在实体建模中,复杂工程模型的设计往往需要大量的布尔操作以及几何造型的系统误差,由此在模型结构中存在大量的裁剪曲面,且曲面之间在边界交界处不满足C2连续.由于在CAD 几何模型中存在大量的短边、窄面、退化边、退化面、非连续光滑边界及尖锐特征等几何特征,这些不理想的几何特征成为网格全自动化生成及CAE 全自动分析计算发展过程中的重要瓶颈和亟待解决的重要问题之一.本文提出的基于TSpline 全自动几何拓扑修复方法主要用于处理非连续光滑边界和一般非理想几何特征,也可适用于退化边、短边、退化面或窄面的处理.对于退化边、短边、退化面或窄面的处理,可采用T-Spline 对相应边或面进行拟合重构操作,但是鉴于退化边、短边、退化面和窄面的特殊性(长度比或面积比与邻边或邻面相差悬殊的孤立单边或孤立单面),从保证CAD模型不失真的前提和简化计算的角度出发,采用特征简化的方法更有优势且效率更高.

图6 钢架焊缝模型中常见的几何噪声Fig.6 Geometric noises of weld in the steel frame

图7 汽车桥壳模型中常见的几何噪声Fig.7 Geometric noises in the automobile axle housing

从高质量网格生成的角度来审视CAD 模型几何拓扑修复问题,本文归纳了4 类影响最终网格单元质量的曲面特征,分类如下:

1)退化边或短边.本文定义的退化边和短边为与邻近连接边长度比相差悬殊的孤立单边,如图8(a),模型中含有退化边,且退化边的长度远小于网格生成时指定的最小长度lmin,在退化边或短边附近生成的网格会局部加密或生成质量较差的狭长网格单元.为了保证几何拓扑的连续性,通过将其合并到与之相邻的曲线,或为了减少计算量,将其视为一个虚点可避免这一问题,如图9 所示.

图8 CAD 模型中常见的非理想几何特征Fig.8 Nonideal geometric features in CAD model

2)退化面或窄面.本文定义的退化面和窄面为与其邻近连接面面积比相差悬殊的孤立单面,如图8(b),模型中含有退化面,且该曲面的面积小于网格生成时指定的最小单元面积Smin,在退化面或窄面附近生成的网格会局部加密,甚至生成网格失败.为了保证几何拓扑的连续性,通过将其合并到与之相邻的曲面,或为了减少计算量,将其视为一条虚边可避免这一问题,如图9 所示.

3)非连续光滑边界.如图8(c),曲面1 和曲面2的边界交接处之间存在非连续的细小短边(高亮显示),网格生成首先需要对几何模型进行边界离散,在进行边离散的同时需要考虑边的曲率、长度等因素,由于这些细小短边的存在会导致边离散效果不理想,其间会产生狭长单元,生成网格质量不理想或者甚至生成网格失败.通过对这些细小短边重新参数化,并按照弦长比进行二叉树细分,在这些细小短边上进行型值点采样,从而达到自适应合理布置样点的目的.利用T-Spline 曲线重新拟合这些细小短边,生成一条新的虚边替代原有的非连续光滑边界,如图10 所示.

图9 退化边、退化面的处理Fig.9 Topology recovery for degenerated edges and degenerated faces

图10 非连续光滑边界的处理Fig.10 Topology recovery for fragmentary smooth boundary edges

4)一般非理想几何特征.如图8(d),在几何模型的B-Rep 中存在细小曲面、尖锐特征、细缝等曲面特征,这些几何特征的局部尺寸小于网格生成时指定的最小特征尺寸,在这些几何特征附近网格会局部加密,导致生成网格数量的量级增加.通过基于T-Spline 曲面的全自动几何拓扑修复方法可避免这一问题.由于一般非理想几何特征的拓扑修复较为复杂,对这一部分算法我们着重进行描述.

3 基于T-Spline 全自动几何拓扑修复

本文基于T-Spline 全自动几何拓扑修复算法,实现对复杂CAD 几何模型中非理想几何特征的自动识别、曲面探测及T-Spline 曲面重构的全自动几何拓扑修复,主要算法流程如下:

为了快速高效全自动进行几何拓扑修复,首先对目标特征进行自动识别,并将自动探测的一般非理想几何特征进行Delaunay 三角化和曲面重新参数化,从而为后序曲面拟合操作提供统一的参数空间.然后在重新参数化的T-Spline 曲面参数空间通过四叉树生成T-网格,并构建T-网格–散点–原几何曲面的映射机制.在T-Spline 曲面拟合操作过程中,所有的几何数据(如坐标、切向量等)都是直接基于原几何曲面进行计算,避免几何误差对拟合结果产生不良影响.

3.1 一般非理想几何特征的自动识别

对于一般非理想几何特征的识别,首先对需要进行拓扑修复的CAD 模型特征进行分块,每个块由它包含的曲面片和域边界组成.域边界是首尾相连的域顶点序列和边组成的环,域边界可表示任意复杂的块的拓扑,一个块可以由多个环构成.本文采用基于面的边界扩展生长法[12]将需要重整的曲面片连接成块状(如图11 所示).在自动识别过程中,将极小面片、狭长面片等小特征曲面作为“种子面片”,块的生长过程就是从边界的起始边出发,以当前面片与相邻面片在公共边界处的光滑度和面片几何中心处的法矢夹角为判断准则并分别给定一合适阈值(如图12 所示).若满足准则,则该边为其相邻域面的共边界,则将相邻面片加入块中,并对原有块的拓扑结构进行编辑(如域顶点的编辑、域边的合并等),更新域顶点和域边的邻面和连接边信息并修改上一步的原始边界;若不满足准则,则遍历下一条边,直到所有边界搜寻结束.最终连接成块的曲面集拥有一个外边界,允许存在多个内边界,而且边界可以是任意形状.该算法对任意拓扑的几何模型均可实施,适应性好,且具有良好的鲁棒性,如图13 所示.

图11 基于面的边界扩展生长法示意图Fig.11 The seeded region growing algorithm images

图12 基于相邻面片间公共边界处的光滑度和几何中心处的法矢夹角判断准则示意图Fig.12 The judgment criteria for seeded region growing algorithm images

图13 一般非理想几何特征的自动识别Fig.13 Automatic identification for nonideal geometric features

1)域顶点的编辑

对域顶点的编辑操作包括删除、合并等.删除算法指当域顶点只有两个相邻边时,可以删除该域顶点,将两个邻边合并成一个域边;当该域顶点有三条及三条以上的相邻边时,则该域顶点至少为三个域面的公共点,则不对这类域顶点进行删除操作.合并算法指域顶点可以与其共域边的其他域顶点进行合并,合并的结果是删除当前域顶点和共域边.

2)域边的合并和优化

在块边界交接处有时会存在非连续的细小短边,且相邻边连接光滑,则将块边界上的这些细小短边重新拟合生成一条新边,并更新块边界的拓扑信息,如前述对非连续光滑边界的处理.

3.2 一般非理想几何特征的Delaunay 三角化

本文在传统Delaunay 三角剖分的基础上改进了一种高效增量插点算法和三角形单元的快速定位搜索算法.Delaunay 三角化[13]方法具有良好的数学理论支撑,且具有算法效率高和生成网格质量好等优点,在不同领域得到广泛的研究和应用.由于平面Delaunay 三角化算法只适用于体网格的生成,不能直接应用于三维曲面面网格生成,本文采用Delaunay 三角化算法在曲面的参数空间生成三角形网格,然后引入黎曼度量将参数三角网映射到三维空间得到三维曲面面网格.在Delaunay 三角化方法中,点的插入占用了绝大部分的运行时间,因此,插点算法的优劣会对Delaunay 网格生成效率产生直接影响.Delaunay 三角剖分的重要特性及三角化过程如图14 和图15 所示.

图14 Delaunay 三角剖分的重要特性Fig.14 The most important characteristic of Delaunay triangulation

图15 Delaunay 三角化过程Fig.15 The process of Delaunay triangulation

本文改进的高效B-W 增量插点算法步骤如下:

1)建立初始网格.构建一个包含所有样本点的矩形,并将矩形初始划分为两个三角形.

2)插入新点,定位基单元.逐个插入新点,插入新点后使用点的快速定位算法(后面详述)确定该点所在的三角形作为基单元.

3)生成空腔.由空外接圆特性,利用基单元搜索找出所有包含P 点的全部三角形,将当前这些三角形删除,形成多边形空腔.

4)生成新单元.形成有效的空腔后,将插入点依次与多边形空腔边界连接,形成新单元,建立点–边–单元之间的拓扑关系,将新生成的单元加入三角网.

5)重复上述步骤,直到所有的节点全部插入为止.

3.2.1 基单元的快速定位

快速定位出任意给定点位于哪一个三角形的算法是Delaunay 三角化的关键.由于这一过程频繁地在网格生成中执行,因此当三角形数目很大时,点的搜索算法将直接影响程序的执行效率.

本文中的K-d 树采用二叉树的数据结构[14],对目标点或基单元进行快速搜索,如图16 所示.在二叉树中,树节点及其邻近节点区域搜索的时间复杂度为O(logn).以参数空间的离散三角网作为背景网格,借助三角形重心数目作为细分规则生成二叉树,建立网格–散点–参数三角网的映射机制,提高目标点或基单元的搜索效率.

图16 基于三角形重心构建的K-d 树Fig.16 Construct the K-d tree based on triangle center

3.2.2 临近单元的快速搜索

基单元确定后,需对所有被插入点打破空外接圆特性的三角形进行快速搜索.本文改进了一种根据点边关系的快速方向定位算法,该方法可以迅速地定位当前插入点位于三角网中的哪一个三角形内部,并在一定程度上提高程序的执行效率,如图17所示.主要算法步骤如下:

1)在给定点附近快速找到一个三角形,作为初始三角形,即基单元;

2)从基单元开始,根据三角形的面积坐标进行快速搜索,并对包含给定点的所有临近三角形进行快速定位,直至搜索到目标三角形.

图17 临近单元的搜索路径Fig.17 The search path for finding the target triangle

3.2.3 改进目标单元快速定位搜索算法效率

为说明本文改进高效增量插点算法和三角形单元的快速定位搜索算法效率明显优于传统搜索算法,本文通过如图18 所示几何模型进行验证.首先在几何模型的参数空间生成Delaunay 三角网,并基于参数空间的三角形重心构建K-d 树.然后,在其参数空间进行随机布点,通过调用C++系统库函数中的随机函数rand(),将随机函数的最大值randmax()设定为106,并将产生的随机点规则化到限定的坐标范围内,调用UG 二次开发函数UF MODL ask point containment(),将不在面上的点(圆孔内部的点)进行剔除,据统计共产生了628 611 个面上的点.表1 给出了分别采用传统方法和本文改进算法的搜索时间对比数据.

图18 验证改进目标单元快速搜索算法效率算例Fig.18 An example to verify the efficiency of target cell searching algorithm

表1 本文改进快速搜索方法与传统搜索方法数据对比Table 1 Comparison between traditional method and our method

3.3 Delaunay 三角化网格曲面的重新参数化

为了快速高效地重建T-Spline 曲面,对Delaunay 三角化网格曲面进行参数化处理,即对自动探测到要处理的非理想几何特征重新参数化得到一个统一的参数空间,且最大程度地保持曲面几何度量(如长度、面积、角度)在映射过程中不发生畸变,如图19 所示.

图19 Delaunay 三角网格曲面的重新参数化Fig.19 Reparameterization of Delaunay triangulation

本文在最小二乘法保形映射(Least squares conformal maps,LSCM)[15]方法理论基础上,通过利用Delaunay 三角化网格曲面的每个三角面片的线性函数去逼近Cauchy-Riemann 方程,通过建立二次逼近误差函数求得在最小二乘意义下的整体的线性拟共形映射,求得三角网格上各个顶点对应的参数值,达到非理想几何特征曲面数据的参数化,为实现T-Spline 曲面重建做准备.LSCM 方法是一种保形的参数化方法,既可以实现保角映射又可以实现保面积映射,如图20 所示.

本文对Delaunay 三角化网格曲面ST=[V,E,T]的每个三角面片建立局部坐标系,其顶点的局部坐标为Zi=(xi,yi),i=1,2,3.已知该三角网格有N个顶点、n个三角面片.设S是其对应参数曲面的方程,则Pi=S(ui,vi),则S的逆映射:U:(x,y,z)→(u,v,0).

图20 三角网格的映射示意图Fig.20 Reparameterization map for a triangle

在每个三角面片上满足柯西–黎曼方程,即

定义目标保角能量函数为

将所有三角形面片上的能量相加,所得的最终能量方程为

其中,i表示第i个三角形,1,2,3 表示该三角形的局部标号;Q表示R3空间三角形面片顶点局部坐标关系式,U表示R2空间三角形面片顶点局部坐标关系式;Mi是投影在Z=0 上的Z1,Z2,Z3组成对应三角形面积的2 倍.

以三角形①为例,Q1,1,Q1,2,Q1,3表示R3空间三角形①顶点(xj,yj)局部坐标关系式(复数形式),u1,1,u1,2,u1,3表示R2空间三角形②顶点(uj,vj)局部坐标关系式(复数形式).

在参数空间至少选取两个基准点作为已知点,即可求解上述Cauchy-Riemann 方程.为最大程度实现高质量的网格曲面重新参数化且保持几何度量在映射过程中不发生畸变,本文选取两个点作为已知点,令u1=x1,v1=y1;u2=x2,v2=y2;则在参数空间还有N −2 个顶点的参数值待求,令Q=(QN)=(QN−2,Q2)为n×N矩阵,n表示三角形个数,N表示所有三角形顶点个数.

其中,qi,j表示第i个三角形第j个顶点关系式值.

为了使目标函数最小,则

其中,上标r 代表复数的实部,i 代表复数的虚部.

对各个矩阵组装完成后对x进行求解,即可得到参数空间各个顶点对应的参数坐标(uj,vj).

3.4 自适应T-Spline 曲面重建算法

本文将T-Spline 和移动最小二乘法(Moving least square,MLS)或径向基函数(Radial basis function,RBF)相结合,提出了一种新的自适应TSpline 曲面重建算法,实现对复杂CAD 几何模型中一般非理想几何特征的全自动几何拓扑修复.与NURBS 曲面相比,T-Spline 曲面允许T 型节点出现,本身具备局部自适应剖分能力,重建曲面能满足任意规定的精度,具有自适应性,并且无需专门的曲面简化的过程,能处理任意拓扑网格,重建过程易操作,重建曲面品质高.由于双三次T-Spline 曲面在工程实际中广泛应用,本文采用双三次T-Spline 曲面.

一般来说,自适应T-Spline 曲面重建流程如图21 所示.首先由Delaunay 三角网格曲面的重新参数化将采样点映射到参数空间,然后利用四叉树或二叉树细分的方法在曲面参数空间自动生成有效的T-网格,最后求解相应控制点,同时在误差较大区域进行局部修正,使拟合误差满足要求.重建实例及拟合曲面误差评价数据如图22 及表2 所示.

表2 自适应T-Spline 拟合曲面误差及网格质量评价Table 2 Relevant experimental data of T-Spline fitting for automatic surface reconstruction

3.4.1 T-网格的构造

在进行T-Spline 曲面重建[16]时,核心问题在于如何构造合理高效的T-网格.一个完整的T-网格包括相应的控制顶点、T-网格子面、网格边和相应的边权值.本文在二维参数空间上利用四叉树细分构造T-网格,细分准则根据点的密集程度和曲面曲率大小决定相应细分次数来决定.T-Spline 网格边及对应的节点间距满足以下规则:

规则1.在T-网格的任意一个子面中,相对的边的节点间距之和必须相等.

图21 自适应T-Spline 曲面重建流程Fig.21 The flowchart of automatic T-Spline surface reconstruction algorithm

图22 自适应T-Spline 曲面重建实例Fig.22 An example of automatic T-Spline surface reconstruction algorithm

规则2.在T-网格中,位于一个子面一条边上的T-交点,如果可以与对边的T-交点相连而不违背规则1,那么该边就包含于T-网格中.

自适应T-网格的构造过程如下:

1)定义最大细分次数以控制细分四边块数量;

2)定义叶子采样点数目阈值或曲面曲率细分终止条件;

3)重复上述步骤,直到满足细分终止条件为止;

4)在T-网格最外层添加重节点和相应虚边,在重节点处添加相应边界条件(包括切矢和扭矢边界条件),虚边赋予一非负值(一般置0)即可;

5)T-网格创建后,T-网格边赋予边权值(节点间距),并对T-网格中各个顶点创建相应节点矢量,确定相应基函数.

3.4.2 T-Spline 曲面拟合算法

设P={p0,p1,p2,···,pn}为T-网格中待求的控制顶点集,Bi(u,v)是由T-网格的节点矢量确定的相应基函数.

定义T-Spline 曲面为

光顺项

重建曲面的目标就是在使E最小化时,令待求控制点Pi的坐标(xi,yi,zi)为X、原始曲面样点Sk的坐标(xk,yk,zk)为b,则上式等价于求解线性方程组AX=b.AX=b两边同乘AT后,等价于求解方程组ATAX=ATb,矩阵A=(aki)=(Bi(uk,vk)).

根据实验经验,σ一般取值在0.01∼0.1 之间,本文取值0.05.其中Auu,Auv,Avv的元素分别为基函数的二阶偏导解此方程组可得到相应T-Spline 曲面所有的控制顶点集P,进而得到T-Spline 曲面.

3.4.3 T-Spline 局部节点插入算法

在重建曲面误差检测过程中,对于某局部区域误差大于给定阈值的型值点,需要将其对应的局部区域进行节点插入[17],形成新的T 网格.此过程利用T-Spline 曲面拟合算法中求解得到的控制点,重新计算局部节点插入后相应的控制点,不再对曲面进行大规模反算,大大减少计算量,提高算法效率.算法的具体步骤如下:

1)误差检测.对于一给定的容许误差ε>0,若原始曲面上一点与T-样条曲面的距离d(pi,S)>ε,则找出pi对应的T-样条曲面区域作为待优化区域.

2)对步骤1)确定的T-网格区域插入新节点.

3)找出插入T-网格节点后基函数改变的控制顶点区域,利用节点插入公式重新计算插入节点后对应的相应控制点.

4)重新分配拟合数据.将求解所得控制点P 赋值到T-网格中.

5)单步优化结束后,需要再对误差进行评价,直到拟合精度满足要求,中止迭代优化.

T-Spline 局部节点插入算法的一个实例及拟合曲面误差评价数据如图23 及表3 所示.

3.4.4 任意自由边界或含空洞的一般非理想几何特征曲面拟合

对于任意自由边界或含空洞的一般非理想几何特征,传统方法通常将任意自由边界曲面在重新参数化过程中通过固定边界为正方形或其他形状作为T-Spline 拟合的参数空间,如图24 所示,但不可避免地会产生长度、面积、角度上的扭曲变形,拟合曲面的几何度量在映射过程中发生畸变.与传统方法相比,本文将T-Spline 与散乱点拟合算法MLS 或RBF 相结合[18],提出了一种新的自适应T-Spline曲面重建算法,在使拟合T-Spline 曲面最大程度地保持原始子曲面的几何度量的同时,提高全自动几何拓扑修复效率,实现对任意自由边界或含空洞的一般非理想几何特征的全自动曲面重构,重建实例及拟合曲面误差评价数据如图25 及表4 所示.

图23 T-Spline 局部节点插入算法实例Fig.23 An example of T-Spline local refinement algorithm

表3 T-Spline 局部节点插入算法拟合曲面误差及网格质量评价Table 3 Relevant experimental data of T-Spline fitting by the local refinement algorithm

表4 任意自由边界非理想几何特征的拟合曲面误差及网格质量评价Table 4 Relevant experimental data of T-Spline fitting for nonideal geometric features with arbitrary boundary

图24 任意自由边界或含空洞的非理想几何特征Fig.24 Nonideal geometric features with complex boundary or holes

图25 任意自由边界的非理想几何特征自动修复Fig.25 Topology recovery for nonideal geometric features with arbitrary boundary

对于任意自由边界或含空洞的一般非理想几何特征,本文仍采用自由边界的重新参数化,最大程度地保持曲面的几何度量在映射过程中不发生畸变.在曲面采样时,对于落在原始子曲面内的型值点,为保证求解精度,通过映射关系直接通过原始子曲面几何信息求得;对于落在原始子曲面外和空洞内的型值点,利用散乱点拟合算法计算相应的型值点,再进行T-Spline 曲面拟合即可.

3.4.5 拟合T-Spline 曲面的误差及网格质量评价

对一般非理想几何特征进行T-Spline 曲面拟合处理后,需要对曲面拟合效果进行误差评价.若曲面拟合精度满足要求,则中止迭代优化;反之,对于不满足拟合精度要求对应区域的控制网格则需要自适应剖分加密,实现对局部细小特征逼近和表征.本文采用最大拟合误差和均方差拟合误差作为误差评价标准.另外,网格的生成质量从侧面也可以反映拟合T-Spline 曲面的好坏,高质量曲面网格有利于提高CAE 分析的计算精度,本文采用网格质量因子对拟合T-Spline 曲面的网格生成质量进行评价.

最大拟合误差和均方差拟合误差的计算如下:

其中,pi为T-Spline 曲面上的点,Si为原始曲面上的点.

三角形的质量α因子[19]定义如下:

其中,α是三角形的面积与边长的平方和之比,0<α≤1,α越大表示质量越好.

四边形的质量因子[20]定义如下:

令四个三角形∆ABC,∆ABD,∆CDB,∆CDA的因子依次为α1,α2,α3,α4,假设α1≥α2≥α3≥α4,则四边形ABCD的质量因子β定义为

其中,0<β≤1,β越大表示四边形质量越好.

三角形及四边形网格的质量因子如图26 所示.

图26 三角形及四边形网格的质量因子Fig.26 The quality factor of triangular and quadrilateral mesh generation

4 全自动几何拓扑修复及网格生成实例

本文利用一种基于T-Spline 的全自动几何拓扑修复方法,对具有几何噪声和细小特征的复杂结构(如汽车桥壳模型、钢架模型等)进行全自动几何拓扑修复,成功处理了几何模型中存在的短边、窄面、退化边、退化面、非连续光滑边界及尖锐特征等常见的不必要几何特征,拟合T-Spline 曲面的精度满足要求,网格生成的质量良好,降低了网格的生成规模,并成功运用到实际应用中解决工程问题.表5∼8 以及图27 和图28 是两个基于T-Spline 的全自动几何拓扑修复及网格生成的实例,从中可以看出,本文算法的可行性及优越性,

图27 汽车桥壳模型全自动几何拓扑修复Fig.27 The results of automobile axle housing model after automatic topology recovery

5 结论

从上述实例可知,本文提出的基于T-Spline 的全自动几何拓扑修复方法在实际应用中具有可行性,相比其他几何拓扑修复方法[21−24],本文提出的方法具有以下特点:

图28 钢架焊缝模型全自动几何拓扑修复Fig.28 The results of steel frame model after automatic topology recovery

1)本文基于T-Spline 全自动几何拓扑修复算法,实现对复杂CAD 几何模型中非理想几何特征的自动识别、曲面探测及T-Spline 曲面重构的全自动几何拓扑修复.对于几何模型中存在的短边、窄面、退化边、退化面、非连续光滑边界及尖锐特征等常见的不必要几何特征,本文的全自动几何拓扑修复算法具有通用性和可行性.

2)开发了一套新的CAD/CAE 集成系统,统一了几何模型与分析模型,所有操作在同一环境下进行,避免了不同系统之间转换造成的几何数据及拓扑信息丢失,实现CAE 与CAD 两者的无缝集成.