程正南，黃璐璐，王閃

（安徽工業(yè)大學(xué)機(jī)械學(xué)院，馬鞍山　243002）

基于Bézier提取NURBS的二維線彈性等幾何分析

程正南，黃璐璐，王閃

（安徽工業(yè)大學(xué)機(jī)械學(xué)院，馬鞍山243002）

等幾何分析方法是用統(tǒng)一的語(yǔ)言描述幾何模型和計(jì)算模型，打通了CAD與CAE模型上的不一致問題。研究基于Bézier提取NURBS的二維線彈性等幾何分析方法，將Bézier提取操作符應(yīng)用于NURBS，使NURBS的基函數(shù)分解成Bernstein多項(xiàng)式，形成基于Berntein多項(xiàng)式的NURBS形式。相對(duì)于傳統(tǒng)有限元分析，此方法Bézier單元高階連續(xù)，求解精度更高；并以薄壁圓形梁模型算例驗(yàn)證了該等幾何分析的實(shí)用性。

等幾何分析；非均勻有理B樣條；Berntein多項(xiàng)式；有限元分析

傳統(tǒng)工程實(shí)踐中，計(jì)算機(jī)輔助設(shè)計(jì)（CAD）與計(jì)算機(jī)輔助工程（CAE）是基于兩個(gè)不同的平臺(tái)，兩者之間的信息是單向傳輸，即在CAD系統(tǒng)中進(jìn)行建模，通過CAD與CAE的接口導(dǎo)入到CAE分析系統(tǒng)內(nèi)，CAE軟件計(jì)算之前通過有限元軟件把CAD模型進(jìn)行前處理，也即網(wǎng)格劃分，這樣無(wú)法與CAD幾何建立起直接的聯(lián)系，這樣的單向信息傳輸使得在計(jì)算分析過程復(fù)雜耗時(shí)。

Hughes等［1］從統(tǒng)一CAD與CAE的角度，提出了一種新的以樣條理論為基礎(chǔ)的數(shù)值計(jì)算方法-等幾何分析［2］（Isogeometry Analysis，IGA）。其思想是將樣條函數(shù)構(gòu)建精確幾何模型［3］，用此樣條基函數(shù)作為形函數(shù)分析該幾何模型，使得設(shè)計(jì)、分析模型采用統(tǒng)一的描述，解決了傳統(tǒng)數(shù)值方法的求解域與幾何模型非一致問題，不僅提高了求解精度也節(jié)約了大量的時(shí)間。等幾何分析優(yōu)勢(shì)多，應(yīng)用廣，廣泛應(yīng)用于板殼分析［4］、流體力學(xué)［5］、電磁場(chǎng)［6］等。雖然等幾何分析方法應(yīng)用和優(yōu)勢(shì)眾多，但對(duì)于特殊結(jié)構(gòu)體中，存在局部不能進(jìn)行細(xì)化等問題，為此，Bazilevs［7］等人及王平［8］等人將T樣條以及PHT樣條相繼用于等幾何分析方法以彌補(bǔ)NURBS的不足。Bézier、B樣條與NURBS的意義眾多［9］，利用Bézier與NURBS優(yōu)勢(shì)互補(bǔ)的性質(zhì)在等幾何分析應(yīng)用中較少。

本文研究了基于Bézier提取NURBS等幾何分析方法，通過節(jié)點(diǎn)嵌入的方式把Bernstein多項(xiàng)式引入到NURBS基函數(shù)中，形成新的基于Bernstein多項(xiàng)式的NURBS形式。并將該方法應(yīng)用到經(jīng)典的二維線彈性實(shí)例中。最后將其計(jì)算結(jié)果與經(jīng)典有限元軟件計(jì)算的結(jié)果進(jìn)行了比較。

1　等幾何分析基礎(chǔ)

1.1B樣條基本概念

節(jié)點(diǎn)矢量Ξ=｛ξ1，ξ2，…ξi，ξi+1，…ξn+pmax+1｝是一個(gè)非減的實(shí)數(shù)序列。其中，ξi稱為節(jié)點(diǎn)，Ξ稱為節(jié)點(diǎn)矢量，p為階次，n為最高階的基函數(shù)個(gè)數(shù)。在等幾何分析中，參數(shù)空間中的單元是由節(jié)點(diǎn)劃分的。用Ni，p（ξ）表示第i個(gè)p次B樣條基函數(shù)，其定義為：

B樣條曲線定義為：

｛Pi｝是控制點(diǎn)，Ni，p（ξ）是定義在節(jié)點(diǎn)矢量上的p次B樣條基函數(shù)。

B樣條基函數(shù)的基本性質(zhì)：非負(fù)性、單位分解性、緊支性、高階連續(xù)性、遞歸求導(dǎo)性、線性無(wú)關(guān)性和非插值性等［9］。

1.2NURBS基本概念

NURBS理論包含了所有B樣條理論的性質(zhì)，同時(shí)它本身也有自己的優(yōu)點(diǎn)。NURBS在B樣條方法的基礎(chǔ)上引入了權(quán)因子與分母，使得在設(shè)計(jì)CAD模型時(shí)更加靈活方便。當(dāng)權(quán)重均為1時(shí)，此時(shí)的NURBS也即B樣條。

NURBS曲線，如圖1（a）所示，在ξ方向上p次的NURBS曲線定義為：

NURBS曲面，如圖1（b）所示，一張?jiān)讦畏较騪次、η方向q次的雙變量分段有理矢值函數(shù)的NURBS曲面為：

圖1　幾何圖形（黑色為控制點(diǎn)）

2　等幾何分析基本原理

2.1Bézier曲線［9］

一條p次Bézier曲線定義為

式中，Bernstein多項(xiàng)式為

其中p為階次，｛Pi｝為控制點(diǎn)，令B1，0（ξ）=1；如果i＜1或i〉p+1，則B1，0（ξ）=0。

2.2Bézier提取運(yùn)算符與應(yīng)用

2.2.1NURBS中節(jié)點(diǎn)嵌入

將待插入的節(jié)點(diǎn)ξ*∈［ξk，）ξk+1插入節(jié)點(diǎn)矢量Ξ形成新的節(jié)點(diǎn)矢量為Ξ=｛ξ1，ξ2，…，ξk，ξ*，ξk+1，…ξn+p+1｝，節(jié)點(diǎn)向量被更新后，以前與節(jié)點(diǎn)矢量Ξ相匹配的控制點(diǎn)和權(quán)重也需更新，得到控制點(diǎn)權(quán)重及其坐標(biāo)更新公式分別為：

2.2.2Bézier提取運(yùn)算符

在B樣條原有的節(jié)點(diǎn)矢量中加入其內(nèi)部已有的節(jié)點(diǎn)，并且其重復(fù)度等于曲線的次數(shù)時(shí)，則B樣條曲線與Bézier曲線形狀一樣，只不過相應(yīng)的控制點(diǎn)會(huì)增加，此時(shí)曲線連續(xù)性和單元之間連續(xù)性并未改變［10］。

Bézier分解是一個(gè)節(jié)點(diǎn)嵌入的操作，Bézier提取運(yùn)算符C是根據(jù)當(dāng)初始節(jié)點(diǎn)矢量嵌入新的節(jié)點(diǎn)矢量后，控制點(diǎn)做出相應(yīng)的變化規(guī)則。相當(dāng)于細(xì)化規(guī)

節(jié)點(diǎn)嵌入后的控制點(diǎn)變化

根據(jù)B樣條曲線公式（2），嵌入節(jié)點(diǎn)后的Bézier

曲線和原始B樣條曲線幾何參數(shù)和形狀未改變，得

由上式得B樣條基函數(shù)與Bernstein多項(xiàng)式的關(guān)系為

其中，wb=CTw（w為NURBS權(quán)重點(diǎn)）因此NURBS的基函數(shù)公式變?yōu)?/p>

得出Bézier控制點(diǎn)與NURBS控制點(diǎn)的關(guān)系為：

則NURBS曲線用Bernstein多項(xiàng)式表示為

根據(jù)單變量提取運(yùn)算符，同樣二變量提取運(yùn)算符可應(yīng)用于NURBS曲面中。二變量提取運(yùn)算符可以被定義為單變量提取運(yùn)算符的張量積形式

2.3線彈性

根據(jù)虛功方程［11］，二維線彈性等效積分弱形式為：

應(yīng)力與應(yīng)變的關(guān)系為

式中，Ri是NURBS基函數(shù)，aei是控制點(diǎn)P的位移，N為控制點(diǎn)總數(shù)。將式（25）、（26）帶入式（24），由于虛位移是任意的，經(jīng)過化簡(jiǎn)得

3　實(shí)例分析

以四分之一薄板圓形梁為例，內(nèi)外半徑分別為a和b。定義參數(shù)，彈性模量E=1MPa，泊松比μ= 0.25，a=5mm，b=10mm。

圖2　四分之一薄壁圓形梁

3.1右側(cè)底部1/4處受到豎直向下的集中力，大小為1N

邊界條件及受力如圖2（a）所示。假設(shè)曲面x方向與y方向節(jié)點(diǎn)矢量均為［0 0 0 1 1 1］。將圓形梁劃分為4x4單元數(shù)，取2階NURBS和Bézier單元進(jìn)行粗略的網(wǎng)格細(xì)分如圖3所示。

圖3　模型離散4x4單元（黑色五角星為控制點(diǎn)）

從圖3可看出Bézier控制點(diǎn)多于NURBS控制點(diǎn)。由于分析均是計(jì)算單元上控制點(diǎn)受力情況，因此這將對(duì)計(jì)算精度有重要影響。

圖4　有限元分析和等幾何分析的第一主應(yīng)力（單位：MPa）

傳統(tǒng)有限元分析與等幾何分析Mises應(yīng)力云圖如圖4所示。從圖4可看出，等幾何分析結(jié)果與經(jīng)典有限元分析的結(jié)果很相近，Mises應(yīng)力最大值均發(fā)生在最左上角處，其值為2.6MPa附近。

3.2底部受到豎直向上的均布線拉力，大小為0.2N/mm

受力及邊界條件如圖2（b）所示。當(dāng)模型的網(wǎng)格不斷細(xì)化，也即模型的單元數(shù)不斷增加時(shí)，有限元結(jié)點(diǎn)和等幾何分析控制點(diǎn)不斷增多，A點(diǎn)的位移計(jì)算也會(huì)趨于穩(wěn)定，其結(jié)果如表1所示。傳統(tǒng)有限元與等幾何計(jì)算方法的數(shù)值解隨單元數(shù)增加的變化趨勢(shì)如圖5所示。

圖5　有限元分析與等幾何分析計(jì)算結(jié)果比較

表1　等幾何分析與有限元分析A點(diǎn)縱向位移（單位：mm）

由圖5可知，在單元數(shù)相等的情況下，等幾何分析遠(yuǎn)比有限單元分析的精度高的多。從等幾何分析看出，階數(shù)等于3，初始劃分結(jié)果接近真實(shí)值，這在復(fù)雜的幾何體表現(xiàn)的更為明顯。并且圓梁模型在等幾何離散過程中始終保持著原有形狀，而有限元離散中模型簡(jiǎn)化為多邊形。

4　結(jié)語(yǔ)

通過節(jié)點(diǎn)嵌入方法，引入Bézier提取NURBS等幾何分析方法，用此方法對(duì)經(jīng)典的1/4圓形梁進(jìn)行分析。從上面可知，該方法對(duì)單元細(xì)分方面更方便，由NURBS基本性質(zhì)，此方法單元高階連續(xù)。最后把等幾何分析在求解的精確度方面與有限分析做了比較，結(jié)果表明等幾何分析的求解精度明顯優(yōu)于有限元分析，且計(jì)算效率和單元的間的連續(xù)性方面也優(yōu)于有限元分析。在未來(lái)工作中等幾何分析方法的優(yōu)勢(shì)將會(huì)更突出。

［1］Hughes T J R，Cottrell J A，Bazilevs Y.Isogeometric analysis：CAD，finite elements，NURBS，exact geometry and mesh refinement［J］.Computer Methods inAppliedMechanicsandEngineering，2005，194（39-41）：4135-4195.

［2］Cottrell J A，Hughes T J R，Bazilevs Y.Isogeometric analysis toward integration of CAD and FEA［M］.New York：John Wiley&Sons，Ltd，2009.

［3］Xu G，Mourrain B，Duvigneau R.Analysis-suitable volume parameterization of multi-block computational domain in isogeometric applications［J］.Computer-AidedDesign，2013，45（2）：395-404.

［4］Benson D J，Bazilevs Y，HSU M C，et al.Isogeometricshellanalysis：theReissner-Mindlinshell［J］.Computer Methods in Applied Mechanics and Engineering，2010，199（5-8）：276-289.

［5］Ming-ChenHSU，AkkermanI，BAZILEVSY. High-performance computing of wind turbine aerodynamics using isogeometric analysis［J］.Comput& Fluids，2011，49（1）：93-100.

［6］Buffa A，Sangalli G，Vazquez R.Isogeometric analysis in electromagnetics：B-splines approximation［J］. Computer Methods in Applied Mechanics and Engineering，2010，199（17-20）：1143-1152.

［7］Bazilevs Y，Calo V M，Cottrell J A，et al.Isogeometric analysis using T-splines［J］.Computer Methods in Applied Mechanics and Engineering，2010，199（5-8）：229-263.

［8］Wang Ping，Xu Jinlan，Deng Jiansong，et al.Adaptive isogeometric analysis using rational PHT-splines［J］.Computer-Aided Design，2011，43（11）：1438-1448.

［9］Piegl L，Tiller W.The NURBS book［M］.2nd ed. New York：Springer，1997：81-116.

［10］Borden，M.J.，M.A.Scott，J.A.Evans，and T.J. R.Hughes.Isogeometric finite element data structures based on Bézier extraction of NURBS［J］.International Journal for Numerical Methods in Engineering，2011，87：15-47.

［11］徐芝綸.彈性力學(xué)［M］.第三版.北京：2002.8：105-128.

Isogeometric Analysis of Dimensional Linear Elasticity Based on Bézier Extraction of NURBS

CHENG Zhengnan，HUANG Lulu，WANG Shan
（School of Mechanical Engineering，Anhui University of Technology，Maanshan，243002）

In isogeometric analysis（IGA），identical language is shared for geometric design and numerical simulation，which solves the inconsistencies between CAD and CAE model.IGA of dimensional linear elasticity is based on Bézier extraction of NURBS in the article，the Bézier extraction operator is used in description of NURBS，then the Bézier extraction operator decomposes the NURBS basis functions to Bernstein polynomials，which forms new NURBS based on Bernstein polynomials.With respect to the traditional finite element analysis（FEA），which generates high degree-continuous of Bézier elements，at the same time the solving accuracy is higher；The thin-walled circular beam model is chosen as an illustrative example to verify the practicability of IGA.

isogeometric analysis；NURBS，Bernstein polynomials；finite element analysis

TP122

A

1672-9870（2016）04-0130-05

2016-03-17

教育部人文社科研究項(xiàng)目（13YJAZH106）

程正南（1989-），男，碩士研究生，E-mail：604834726@qq.com