计算力学学报

参数	本文结果	文献[8]（Nitsche法）	文献[10]（光顺rTBS法）	ABA QUS
W_max/mm	－0.302	－0.301	－0.308	－0.386

参数	本文结果	文献[8]（Nitsche法）	文献[10]（光顺rTBS法）	ABA QUS
W_max/mm	－0.302	－0.301	－0.308	－0.386

基于Bézier三角形的曲面网格划分及等几何板壳分析

PDF下载

罗飞 ¹^,² , 郭玉杰 ² , 孙芳斌 ² , 张兴 ¹

计算力学学报 | 研究论文 2025,42(5): 825-830

收起

计算力学学报 | 研究论文 2025, 42(5): 825-830

基于Bézier三角形的曲面网格划分及等几何板壳分析

全屏

罗飞¹^,², 郭玉杰², 孙芳斌², 张兴¹

作者信息

^1.中国飞机强度研究所　强度与结构完整性重点实验室，西安　710065

^2.南京航空航天大学　飞行器先进设计技术国防重点学科实验室，南京　210016

郭玉杰^*（1986-），男，博士，副教授（E-mail：yujieguo@nuaa.edu.cn）.

Bézier triangular surface mesh generation and isogeometric plate and shell analysis

Fei LUO¹^,², Yujie GUO², Fangbin SUN², Xing ZHANG¹

Affiliations

^1.Key Laboratory of Strength and Structural Integrity, China Aircraft Strength Research Institute, Xi'an 710065, China

^2.National Defense Key Discipline Laboratory of Advanced Design Technology of Aircraft, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

出版时间: 2025-10-28 doi: 10.7511/jslx20231118001

文章导航

摘要

收起

针对单张NURBS曲面提出了一套考虑裁剪与非裁剪形式的Bézier三角形网格划分方法，并将其应用于等几何Kirchhoff-Love壳结构的分析中。首先，将NURBS曲面插值为Bézier曲面；其次，在参数域内计算裁剪曲线与各个参数节点的拓扑关系，并在裁剪曲线上取点，生成参数域内的Bézier轮廓曲线集合；基于该轮廓曲线集合，生成参数域上的三角形网格；最后，采用映射法生成物理域上的Bézier三角网格。通过三个模型验证了算法的适应性与鲁棒性，并计算了网格质量，结果表明整体网格质量良好；在此基础上，以Scordelis-Lo's Roof壳模型为例，采用罚函数法施加Kirchhoff-Love壳单元间的转动约束，验证了基于Bézier三角形的Kirchhoff-Love壳单元的精度。

关键词

等几何分析 / Bézier三角网格划分 / Bézier曲面映射 / 曲线裁剪 / Kirchhoff-love壳

Abstract

收起

This paper presents a Bezier triangle meshing method that considers both clipped and non-clipped forms for a single NURBS surface. The proposed method is applied to analyze isogeometric Kirchhoff-Love shell structures. The process begins by interpolating NURBS surfaces into Bezier surfaces. Subsequently, the topological relationship between the clipping curve and each parameter node is calculated within the parameter domain. A Bezier contour curve set is then generated in the parameter domain by selecting points along the clipping curve. Utilizing this contour curve set, a triangular mesh is generated in the parameter domain. Finally, the Bezier triangle mesh in the physical domain is created through a mapping method. The adaptability and robustness of the algorithm are verified through three models, and the mesh quality is assessed. The results demonstrate favorable overall mesh quality. Building upon this foundation, the paper illustrates the application of a rotation constraint between Kirchhoff-Love shell elements, using the penalty function method with Scordelis-Lo's Roof shell model as an example. The accuracy of Kirchhoff-Love shell elements based on Bezier triangles is subsequently validated.

Key words

isogeometric analysis / Bézier triangular mesh / Bézier surface mapping / Curves trim / Kirchhoff-love shell

引用本文

罗飞, 郭玉杰, 孙芳斌, 张兴. 基于Bézier三角形的曲面网格划分及等几何板壳分析. 计算力学学报, 2025 , 42 (5) : 825 -830 . DOI: 10.7511/jslx20231118001

Fei LUO, Yujie GUO, Fangbin SUN, Xing ZHANG. Bézier triangular surface mesh generation and isogeometric plate and shell analysis[J]. Chinese Journal of Computational Mechanics, 2025 , 42 (5) : 825 -830 . DOI: 10.7511/jslx20231118001

正文

收起

1　引言

收起

Hughes等^[1]基于等参单元的概念提出了等几何分析（IGA）方法，IGA不仅采用非均匀有理B样条（NURBS）描述几何模型，还以NURBS函数插值位移场^[2]，因此IGA消除了几何离散误差，并且IGA还具有高阶连续、几何精确等特性，使得其非常适于曲面等薄壁结构的分析。NURBS由于具有张量积的性质，因此只能描述较为简单的形状。但工程结构常较为复杂，难以采用单张NURBS张量积曲面进行描述，进而需要采用裁剪、布尔、拼接等操作方式才能构建所需曲面构型。针对裁剪曲面，采用四边形单元分析时，需要对裁剪单元进行额外处理；Kim等^[3]提出了一种处理二维剪裁曲面的NURBS增强单元，在积分时，根据裁剪单元的形状对其进行分类，并连接对应的节点以及裁剪交点，从而在曲面参数域中构建裁剪单元的三角网格。采用三角形单元进行离散比四边形单元更加灵活，因此国内外许多学者对其进行了研究。Farin^[4]首先在计算机辅助设计（CAGD）中提出伯恩斯坦Bézier三角形的概念，给出了两个三角面片连续的约束条件，但是该公式并不能保证适用于任意阶次的三角形单元。Lasser^[5]提出了从一个Bézier四边形面片（母片）中提取出任意三角形面片（子片）的方法，通过构造辅助控制点，可以将在Bézier张量积曲面参数域内的子三角形映射到物理域，通过调整参数域中子片以及母片的阶数，能够生成任意阶数的Bézier三角曲面。

曲壳是工程中常见的结构形式，在对曲壳进行力学分析时，常采用Kirchhoff-Love壳单元，由于不考虑横向剪切，Kirchhoff-Love壳单元能够以较少自由度模拟曲壳的位移场，节省计算时间。基于NURBS函数的高阶连续特点，Kiendl等^[6]首先发展了NURBS基的Kirchhoff-Love壳单元，采用将单元交界面上的自由度直接耦合在一起的方式，满足了Kirchhoff-Love壳单元在单元交界面所需的C¹连续性约束条件。Schuβ等^[7]基于拉格朗日乘子以及修改基函数的方式，提出了一种针对多片Kirchhoff-Love壳单元的耦合方法。Guo等^[8]采用具有变分一致特性的Nitsche方法实现了薄壳结构的本质边界条件施加，提升了薄壳结构分析的精度。Herrema等^[9]采用罚函数法施加多片NURBS的界面约束，误差在工程允许范围内。Zareh等^[10]在已有Bézier三角网格的基础上，通过Clough-Tocher和Powell-Sabin方法进一步将大网格划分为若干小网格以满足连续性要求。此外，等几何分析还应用到边界元求解^[11]等问题中。

本文针对单张NURBS曲面的三角网格划分进行研究，提出了一套考虑裁剪与非裁剪的NURBS曲面Bézier三角网格自动划分方法；进一步将划分的Bézier三角形应用到薄壳计算中，研究了基于Kirchhoff-Love薄壳理论的等几何分析方法。

2　理论及方法

收起

2.1　B样条和NURBS曲线曲面

B样条基函数以递归形式定义在一个单调不减的节点向量U={u₁，…，u_n+p+₁}上，对于每个节点值，总有u_i≤u_i+₁（i=1，…，n+p），n是控制点数量，p是曲线的阶数。

(1)

(2)

一条NURBS曲线可定义为

(3)

P_i是曲线的第i个控制点，w_i是权重。采用张量积形式可以定义一张NURBS曲面为

(4)

(5)

{P_ij}形成两个方向的控制点网络，{ω_ij}是权因子，{R_ij}是非均匀有理B样条函数，p和q是NURBS曲面在两个参数方向的阶数^[12]。

2.2　Bézier三角形

两变量的伯恩斯坦基定义为

(6)

采用面积坐标，其可定义为

(7)

在此基础上，Bézier三角形可定义为

(8)

T（ζ₁，ζ₂，ζ₃）是Bézier三角形上的一点，ζ₁、ζ₂、ζ3是参数三角形上的面积坐标^[2]，其中i=（i，j，k），i、j、k是对应控制点的下标，P_i是对应的控制点。

2.3　映射算法

选择文献[5]的映射公式将Bézier四边形面片映射为若干张Bézier三角形面片，相关公式为

(9)

(10)

S（ζ₁，ζ₂，ζ₃）是曲面三角形面片上的控制点，计算S（ζ₁，ζ₂，ζ₃）首先需要在Bézier四边形面片的参数域内选取一个参数三角形，记参数三角形的阶数为N，则最终映射得到的物理三角形的阶数为N×（l+m），其中l与m分别是Bézier四边形面片两个方向的阶次。本文参数三角形的阶数只考虑N=1的情况，

是辅助控制点的插值系数，

是辅助构造控制点，其由开花算法计算得到，

，

是中间迭代系数。

2.4　Bézier曲面三角形壳单元的应用

2.4.1　等几何离散

采用划分的Bézier三角网格对K-L壳进行离散，K-L壳^[13]需要计算二阶微分，因此在单元界面上要保持一阶连续，而Bézier三角形在单元界面上为C⁰连续，因此本文基于罚函数方法施加三角形界面上的转动约束。由等几何的概念可知，采用伯恩斯坦-Bézier基进行单元离散后，域上一点处的位移，可以由单元的控制点位移插值得到，其定义如下，其中，u是壳中面层上一点的位移，U_i是控制点位移。

(11)

在小变形线性假设下，单元的刚度矩阵为

(12)

耦合刚度矩阵为

(13)

其中f是罚因子，Γ是对应的单元交界面。{φ_T}是界面上两个单元沿切线方向转角的差值。

3　曲面三角形网格划分流程

收起

经典NURBS曲面三角形网格划分方法在处理裁剪区域的思路为在参数域中的裁剪曲线上进行取点，然后将所取的点向参数域四条边界映射，连接对应的映射点，形成若干矩形单元格，随后对矩形单元格一分为二，得到参数域三角形，最终经过映射得到物理域三角形^[2]，该思路简单，容易实现，但是对曲面进行三角形网格剖分时容易生成畸形网格。

本文提到的NURBS曲面三角形网格划分算法思路为先形成参数域中的贝塞尔曲线集，再结合二维网格划分算法进行参数域网格划分，因为参数域网格为自动生成，因此保证了模型的网格质量。本文算法共分为六个步骤。

（1）在犀牛软件中进行模型创建，保存模型为3dm格式并由程序读取该模型。

（2）对曲面在两个参数方向进行节点插值，将NURBS面片转换为Bézier面片；假设NURBS面片在两个方向的节点长度为a×b（不考虑重复度，后续提到的节点重复度均不考虑重复度），则经过节点插入转换后可以得到（a－1）×（b－1）张Bézier面片。

（3）对参数域的所有节点进行遍历，判断节点是否在裁剪曲线内，并进行分类，生成（a－1）×（b－1）个初始节点单元（每个初始节点单元包含四个节点，其存储在数组A₁中），与之对应物理域上（a－1）×（b－1）张Bézier面片。

（4）对所有初始节点单元遍历，每个节点单元的四条边分别与裁剪曲线求交；若有交点，则将交点保存到一个4×n_i的二维数组A₂中，n_i是对应第i边上的裁剪交点个数；再将A₂中所有交点在各条边上按逆时针的拓扑结构进行排序；图1是数组A₂中的交点位置变化情况。

（5）对于每一条边，将数组A₁与A₂中的点按照对应的拓扑关系（该点是否在模型的参数域内）轮流填入到一个1 m的数组A₃，其中，m为数组A₁中在单元上的点数量与数组A₂中点的数量之和。在对应裁剪曲线上取点，并按逆时针连接节点单元所有的点，使其成为一条闭合的一次Bézier曲线。对所有节点单元遍历，得到（a－1）×（b－1）条一次Bézier轮廓曲线。将这些曲线输出为.spline格式文件，并基于该轮廓，通过二维网格划分程序XMESH^[14]生成二维三角网格，最后返回模型参数域内的网格信息。

（6）采用式（9）将每一个二维网格映射到物理域中。若网格不满足要求，对模型进行升阶或节点插入，然后跳转到第（2）步，若满足要求，则网格划分结束。图2为裁剪单元边界处理流程。

4　网格划分算例

收起

将模型的形状质量系数作为衡量网格质量的一项指标；对单个三角形，其形状质量系数为

(14)

其中φ的取值在（0，1）之间，S_Δ表示该三角形的面积，L₁、L₂、L₃分别是三角形对应三边的长度；当φ值越大，网格越趋近于等边三角形，表明网格形状质量越好。对于整体模型，可将所有网格的形状质量系数的平均值以及标准差作为评价指标，具体可表示为

(15)

(16)

其中

是网格的整体质量系数，δ_G是整体质量系数的标准差，反映整个网格质量的差异性；M是模型的网格总数。

4.1　开孔半圆柱模型

开孔壳模型在U和V方向的基函数阶数分别为2阶和1阶，节点向量长度均为5；其NURBS几何描述如图3所示。划分的三角形单元如图4所示。由图4可知，在裁剪区域附近，网格较密，在模型非裁剪边界附近，网格较粗。网格质量系数为0.9137，标准差为0.0673。

4.2　多孔半圆柱壳模型

多孔半圆柱壳模型包含五条NURBS裁剪曲线，在U方向的阶数为2阶，节点向量长度为5，V方向为1阶，节点向量长度为7。利用本文算法可得最终网格划分模型（如图5所示）。该模型形状质量系数为0.9175，大于0.9，标准差为0.0723，低于0.1，有较高的网格质量。

4.3　锥形曲壳

本算例研究一锥形曲壳的网格划分问题。图6为采用本文算法得到的物理域网格，图7为对照网格^[16]。通过对比发现，本文模型所得网格规整性略好，且网格数量略少，数量为文献[16]的2/3。

5　Scordelis-Lo's Roof壳结构分析

收起

将上述三角形网格算法应用于壳结构的分析，采用算例为Scordelis-Lo's Roof壳，其几何描述如图8所示，该壳结构长度为50 mm，半径为25 mm，圆弧张开角度为80°，杨氏模量E=432 e+8 MPa，泊松比ν=0，该壳结构受Z方向的均布力，大小为p=－90 N/mm²，位移边界条件为曲边简支。本算例罚因子取为1.0e+8。图9是p=5时的结构Z向位移云图；图10为本文算法计算得到的Roof壳的弯曲力M₁₂，图11为ABAUQS计算得到的Roof壳的弯曲力M₁₂，可以看出，二者在分布形式上基本一致。图12为当单元阶次p从3升高为5时（单元数为337），该壳结构的Z向最大位移的收敛曲线，可知该收敛曲线波动不大，最终结果收敛于－0.302 mm，与其他文献和ABAUQS的对比结果列入表2；由表2可以看出，本文结果与参考文献中的结果相近，与ABAQUS有所差别。其原因为ABAQU目前采用的方法仍是经典有限元法，本文以及其余两篇文献均采用了IGA方法，其消除了几何离散误差。

图13为Scordelis-Lo's Roof壳Z方向的位移随单元数增加的收敛曲线，可知基于本文的Bezier三角形单元计算结果收敛更快，其最终收敛到－0.303 mm。通过本算例可知本文网格划分算法以及等几何壳单元的有效性及精确性。

6　结论

收起

本文针对单张NURBS曲面，通过将NURBS曲面以节点插值的方式转化为Bézier面片，并将曲面的裁剪问题转化为与Bézier面片对应的二维张量积参数域四条边的相交问题，从而简化了曲线裁剪的复杂性，同时结合成熟网格划分算法保证了网格划分结果的质量。通过三个算例展示了本文提出的网格算法的通用性。采用罚函数法将三角单元间的转角进行耦合，解决了单元间转角不连续的问题，耦合方式简单且能满足精度要求，通过计算Scordelis-Lo's Roof模型验证了本文壳单元算法的正确性。

本文还存在一些不足之处，如在网格划分过程中，为了保证相邻Bézier面片划分出的三角网格在相邻边界处具有相同的控制点分布，本文将所有Bézier面片（剪裁与非裁剪面片）的轮廓作为输入条件进行二维网格划分，导致二维网格划分过程中约束边界过多，从而造成一定的网格划分困难。

基金

收起

国家自然科学基金面上项目(11972187)
中央高校基本科研业务费(NT2022004)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

Hughes

T J R

, Cottrell

J A

, Bazilevs

. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39-41): 4135-4195.

[2]

Xia

S T

, Wang

X L

, Qian

X P

. Continuity and convergence in rational triangular Bézier spline based isogeometric analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 297: 292-324.

[3]

Kim

H J

, Seo

Y D

, Youn

S K

. Isogeometric analysis for trimmed CAD surfaces[J]. Computer Methods in Applied Mechanics and Engineering, 2009, 198(37-40): 2982-2995.

[4]

Farin

. Triangular Bernstein-Bézier patches[J]. Computer Aided Geometric Design, 1986, 3(2): 83-127.

[5]

Lasser

. Triangular subpatches of rectangular Bézier surfaces[J]. Computers & Mathematics with Applications, 2008, 55(8): 1706-1719.

[6]

Kiendl

, Bletzinger

K U

, Linhard

, et al. Isogeometric shell analysis with Kirchhoff-love elements[J]. Computer Methods in Applied Mechanics and Engineering, 2009, 198(49-52): 3902-3914.

[7]

Schuβ

, Dittmann

, Wohlmuth

, et al. Multi-patch isogeometric analysis for Kirchhoff-love shell elements[J]. Computer Methods in Applied Mechanics and Engineering, 2019, 349: 91-116.

[8]

Guo

Y J

, Ruess

. Nitsche's method for a coupling of isogeometric thin shells and blended shell structures[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 881-905.

[9]

Herrema

A J

, Johnson

E L

, Proserpio

, et al. Penalty coupling of non-matching isogeometric Kirchhoff-love shell patches with application to composite wind turbine blades[J]. Computer Methods in Applied Mechanics and Engineering, 2019, 346: 810-840.

[10]

Zareh

, Qian

X P

. Kirchhoff-love shell formulation based on triangular isogeometric analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2019, 347: 853-873.

[11]

韩少燕, 姜人伟, 高汝鑫, 等. 基于等几何分析的边界元法求解二维Laplace方程[J]. 计算力学学报, 2023, 40(1): 105-110.

(HAN

Shao-yan

, JIANG

Ren-wei

, GAO

Ru-xin

, et al. Isogeometric boundary element analysis for 2D Laplace equations[J]. Chinese Journal of Computational Mechanics, 2023, 40(1): 105-110.

[12]

Piegl

, Tiller

. The NURBS book[M]. Springer Science & Business Media, 1996.

[13]

Guo

Y J

, Ruess

. Weak Dirichlet boundary conditions for trimmed thin isogeometric shells[J]. Computers & Mathematics with Applications, 2015, 70(7): 1425-1440.

[14]

Engvall

, Evans

J A

. Isogeometric triangular bernstein-bézier discretizations: Automatic mesh generation and geometrically exact finite element analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2016, 304: 378-407.

[15]

Breitenberger

, Apostolatos

, Philipp

, et al. Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shell structures[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 401-457.

[16]

丁慧, 罗尧治. 自由形态网壳结构网格生成的等参线分割法[J]. 浙江大学学报(工学版), 2014, 48(10): 1795-1801, 1834.

(DING

Hui

, LUO

Yao-zhi

. Isoparametric line dividing grid generation method for free-form lattice shells[J]. Journal of Zhejiang University (Engineering Science), 2014, 48(10): 1795-1801, 1834. (in Chinese)

2025年第42卷第5期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.7511/jslx20231118001

接收时间：2023-11-18
首发时间：2026-03-24
出版时间：2025-10-28

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2023-11-18
修回日期：2024-02-07

基金

国家自然科学基金面上项目(11972187)

中央高校基本科研业务费(NT2022004)

作者信息

^1.中国飞机强度研究所　强度与结构完整性重点实验室，西安　710065

^2.南京航空航天大学　飞行器先进设计技术国防重点学科实验室，南京　210016

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/jslxxb/CN/10.7511/jslx20231118001

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT