孫振山，毛虎平，賈冰潮，張穌寧

(中北大學(xué) 能源與動力工程學(xué)院，山西 太原 030051)

0 引 言

有限元法是求解靜力學(xué)問題常用的方法，在平面結(jié)構(gòu)分析中也不例外. 平面靜力學(xué)分析以力學(xué)原理和數(shù)值分析算法為基礎(chǔ)，通過構(gòu)造整體剛度矩陣建立本構(gòu)方程，在指定的設(shè)計要求和約束條件下，計算其位移、應(yīng)力等，進而為實際的生產(chǎn)加工提供參考.

有限元分析法是機械、土木、汽車、船舶及航空航天等領(lǐng)域中平面靜力學(xué)分析普遍采用的基礎(chǔ)手段. 對于結(jié)構(gòu)復(fù)雜且需要精確表示的平面，目前依賴于傳統(tǒng)的建模方法，利用CAD與CAE建模，然后將模型導(dǎo)入有限元分析軟件進行計算. 隨著設(shè)計模型復(fù)雜程度的增加，對有限元網(wǎng)格劃分的質(zhì)量、計算速度、仿真精度有了更高的要求. 在有限元分析過程中網(wǎng)格劃分占據(jù)約80%的時間，而且有限元分析所建立的幾何模型是近似的而非精確的表示，所引入的幾何誤差對一些關(guān)鍵部件的影響較大. 通常，工程師會通過不斷細化網(wǎng)格來達到所需要的計算精度，但這樣不僅會產(chǎn)生大量的網(wǎng)格單元，也會降低計算速度，所需要的計算成本也會隨之增加. 有限元分析和設(shè)計的過程需要計算機不斷地修改和重新劃分網(wǎng)格，此過程需要與CAD模型頻繁交互，過程繁瑣. 此外，在有限元分析過程中以C0連續(xù)性為代表的常規(guī)單元無法滿足高精度高性能的分析要求.

研究者希望超越現(xiàn)有網(wǎng)格劃分的分析方法，直接基于精確的幾何模型進行建模、分析和優(yōu)化，將所需要的幾何模型通過控制點精確地表示出來，但此方法相較于成熟的有限元分析方法而言也存在一些缺點. 2005年，Hughes等[1]提出了基于NURBS(Non-Uniform Rational B-Splines)基函數(shù)的等幾何分析(Isogeometric Analysis, IGA)數(shù)值計算方法，該方法將所需要的幾何模型通過NURBS基函數(shù)建立精確的參數(shù)化模型，通過對節(jié)點、節(jié)點矢量、控制點及控制點權(quán)重的計算，選擇合適的分段多項式基函數(shù)，構(gòu)造出分段有理多項式曲線以及非有理的張量積曲面.

等幾何分析方法以樣條(如B樣條、T樣條等)函數(shù)為基礎(chǔ)，將CAD幾何模型與結(jié)構(gòu)分析緊密結(jié)合，從而實現(xiàn)精確的幾何分析[2-3]. 將等幾何分析應(yīng)用于平面靜力學(xué)分析時，直接建立起了IGA與平面邊界的聯(lián)系，這有利于模型參數(shù)化[4-5]. 大部分的CAD和CAE軟件在建模時采用布爾和運算，這將會產(chǎn)生大量的剪切曲面，使得設(shè)計難度增加，而等幾何分析直接基于物理模型建模，省去了大量的剪切操作. 與傳統(tǒng)的有限元分析相比，等幾何分析在優(yōu)化過程中有如下優(yōu)點：(1) 通過精確的NURBS基函數(shù)對模型進行建模分析，從源頭上消除了有限元分析所引入的幾何誤差，使得計算精度更高；(2) 避免了在有限元建模過程中與計算機進行繁瑣交互的問題；(3) 等幾何分析將計算所得的控制點作為設(shè)計變量，在保證分析精度要求的同時具有更好的連續(xù)性，從而有效克服了有限元分析方法中C0連續(xù)性所引起的計算困難[6-8]的問題. 目前，國內(nèi)外學(xué)者對等幾何分析也開展了大量相關(guān)研究，如袁沛等[9]采用了一種基于等幾何分析的非結(jié)構(gòu)化T樣條構(gòu)建復(fù)雜模型的方法；賈悅等[10]結(jié)合傳統(tǒng)等幾何分析配置法和傳統(tǒng)等幾何分析伽遼金法，提出了一種基于PHT-樣條函數(shù)的加強等幾何分析配置方法；張中林等[11]引入了基于推廣B樣條的等幾何分析方法對二維彈性問題進行邊界元分析；胡丹丹等[12]針對等幾何分析中復(fù)雜物理區(qū)域上的偏微分方程求解問題，提出了多片參數(shù)域上雙三次樣條投影映射的方法. 在等幾何法應(yīng)用方面，李明超等[13]基于等幾何分析開展了水利結(jié)構(gòu)數(shù)值仿真計算新方法的研究；柏硌等[14]為解決零虧格邊界網(wǎng)格模型的T樣條實體重建問題，提出了一種基于八叉樹細分和漸進迭代最小二乘擬合算法的T樣條實體構(gòu)建算法；陳龍等[15]將等幾何分析應(yīng)用于單齒嚙合的齒輪接觸問題上，此法在解決接觸仿真分析方面具有優(yōu)良的性能.

綜上所述，等幾何分析由于精確的參數(shù)化建模，使得其分析更加方便，求解效率更高，所需要的計算成本相對更低，因此已被廣泛應(yīng)用于許多領(lǐng)域. 但與成熟的有限元技術(shù)相比，等幾何法起步較晚，目前還沒有一款商用軟件將等幾何分析成功應(yīng)用于多種工程實例. 節(jié)點向量是等幾何分析的關(guān)鍵要素之一，節(jié)點向量的等距分布可能會產(chǎn)生奇異方程組，降低矩陣計算效率. 本文借助matlab平臺，將平均值法融入到節(jié)點向量的計算中，并將平均值等幾何分析方法應(yīng)用于經(jīng)典力學(xué)問題，建立了平面靜力學(xué)系統(tǒng)化程序. 采用平均值方法對節(jié)點向量進行計算，以NURBS基函數(shù)對平面結(jié)構(gòu)進行參數(shù)化建模，以二維平面靜力學(xué)問題為例展現(xiàn)此算法的可行性.

1 等幾何分析

1.1 NURBS基函數(shù)

非均勻有理B樣條(NURBS，Non-Uniform Rational B-Splines)以B樣條方法為基礎(chǔ)，具有非均勻性、有理性和B樣條三個特征，即控制點的影響范圍可以改變、有理多項式定義幾何結(jié)構(gòu)和控制多邊形或控制多面體來構(gòu)造曲線曲面.k次有理基函數(shù)

(1)

式中：ωi為對應(yīng)控制點的權(quán)值；Ni,k(ξ)為k次B樣條基函數(shù)，由de Boor-Cox[16-18]提出的遞推公式定義為

(2)

為了表示面，二維NURBS基函數(shù)可表示為

(3)

NURBS曲面是一個雙線性分段有理函數(shù)，表示為

(4)

式中：Pi,j為控制點.

1.2 節(jié)點向量計算方法

建立一個系數(shù)矩陣為(n+1)×(n+1)的線性方程組

(5)

式中：{Qk}為一組數(shù)據(jù)點.

(6)

(7)

u0=…=up=0,um-p=…=um=1,

(8)

式中：n為控制點的最大下標(控制點的個數(shù)為n+1)；p為曲線的次數(shù)；m+1為節(jié)點個數(shù)(控制點、節(jié)點從0開始計數(shù)).

1.3 等幾何分析

圖1 所示為NURBS曲面等幾何分析示意圖，可以看出，等幾何分析首先通過控制點建立物理空間，控制點部分在曲面上，部分不在曲面上，然后采用等參元思想將物理空間轉(zhuǎn)化到[-1，1]的父空間，再將父空間轉(zhuǎn)化到[0，1]參數(shù)空間，得到節(jié)點向量，最后建立等幾何分析模型.

圖1 NURBS曲面等幾何分析示意圖[2]Fig.1 Schematic diagram of geometric analysis of NURBS surface

采用NURBS基函數(shù)離散時，涉及到參數(shù)變換問題. NURBS基函數(shù)定義在標準[0,1]空間，而物理空間不在此空間，因此需要進行轉(zhuǎn)換. 數(shù)值積分一般采用Gauss積分，而其積分空間在標準[-1,1]空間，需要從[-1,1]空間轉(zhuǎn)換到[0,1]空間.

(9)

對應(yīng)的雅可比公式為

(10)

從物理空間[a,b]到標準[0,1]空間轉(zhuǎn)化的雅可比公式為

(11)

在優(yōu)化中，用等幾何分析法代替有限元法對物理場進行數(shù)值計算. 采用NURBS基函數(shù)可以直接表示幾何模型并將其作為結(jié)構(gòu)分析的形函數(shù). 參數(shù)坐標(ξ,η)中的物理量可以通過控制點的值獲得

(12)

對于一個線彈性問題，控制方程離散后獲得

Ku=F,

(13)

式中：K表示全局剛度矩陣;u表示位移向量;F表示與控制點相關(guān)的外力.等幾何分析中，同樣通過單元剛度矩陣Ke來構(gòu)造全局剛度矩陣K，Ke表示為

(14)

(15)

而

(16)

式中：Ni為NURBS單元的基函數(shù);J1為

(17)

雅可比矩陣J2表示為

(18)

對于完全彈性的各向同性材料，彈性矩陣D表示為

(19)

式中：E為彈性模量；μ為泊松比.

將單元剛度矩陣Ke逐個按索引編號組裝到整體矩陣中，組裝過程與傳統(tǒng)的有限元方法類似，得到關(guān)于所有控制點變量的整體方程

Ku=F.

(20)

整體剛度矩陣

(21)

整體等效控制點載荷向量

(22)

2 基于平均值等幾何法的靜力學(xué)分析

等幾何分析與有限元分析類似，采取一種類等參元的思想. 不同之處是等幾何分析一般采用NURBS基函數(shù)，而有限元法采用拉格朗日基函數(shù). NURBS基函數(shù)具有局部支撐性、規(guī)范性、可微性和非負性等性質(zhì). 其中，可微性使等幾何分析建立的曲線曲面更加光滑，連續(xù)性更好，同時它又可以出現(xiàn)尖點. NURBS基函數(shù)在節(jié)點區(qū)間內(nèi)部存在任意階導(dǎo)數(shù)，在節(jié)點處，它是p-k次連續(xù)可微的，k為該節(jié)點的重復(fù)度. 等幾何分析在建模過程中直接基于精確的幾何模型建模，可以避免有限元分析在劃分網(wǎng)格時可能引入的誤差. 圖2 為二維平面靜力學(xué)等幾何分析的基本流程圖.

圖2 等幾何分析平面靜力學(xué)分析流程圖Fig.2 Flow chart of plane statics analysis for isogeometric analysis

本文采用matlab編程實現(xiàn)二維平面問題的等幾何分析. 平面問題的等幾何分析屬于線性靜態(tài)分析，首先定義平面的材料屬性，然后固定該平面的一端使其x軸位移為0，對另一端施加集中力載荷. 之后將平面模型導(dǎo)入，將平面模型從物理域向幾何域轉(zhuǎn)化，并參數(shù)化幾何域模型. 根據(jù)算法計算平面模型的節(jié)點向量以及基函數(shù)，再根據(jù)要求對其模型進行節(jié)點插入、網(wǎng)格細分. 引入等參元的思想，計算每個等幾何單元的剛度矩陣，再將它們整合到整體剛度矩陣，把邊界條件施加到整體剛度矩陣中，最終求得平面的位移、應(yīng)力和應(yīng)變.

算法具體實現(xiàn)過程如下：

1) 采用式(8)計算節(jié)點向量；

2) 定義基本變量，包括控制點、基函數(shù)次數(shù)、節(jié)點向量和材料屬性,調(diào)用彈性矩陣算法；

3) 進行單元細化，調(diào)用單元細化算法(hRefinement2d)和生成等幾何二維網(wǎng)格算法(generateIGA2DMesh)；

4) 提取控制點數(shù)和單元數(shù)量，為節(jié)點編號；

5) 通過調(diào)用plotPlate算法繪制幾何模型；

6) 找尋相應(yīng)的邊界節(jié)點編號和受力單元節(jié)點編號對其進行初始化設(shè)置；

7) Gauss積分；

8) 編寫循環(huán)單元剛度矩陣，并將之組裝到整體剛度矩陣；

9) 對受力點施加節(jié)點力F；

10) 進行邊界條件處理，簡化剛度矩陣；

11) 通過u=KF求解位移；

12) 后處理，繪制位移云圖.

3 算例分析

為了驗證改進后的等幾何法的可行性，選擇一個普通的平面靜力學(xué)模型進行計算.

3.1 懸臂梁問題

懸臂梁結(jié)構(gòu)如圖3 所示，懸臂梁長L=30 m，梁的高度h=6 m，單位厚度. 平面左端固定，右端施加集中力載荷F=-10 000 N, 方向垂直向下. 材料為線彈性，彈性模量E=210 GPa，泊松比μ=0.3.

圖3 懸臂梁結(jié)構(gòu)Fig.3 Cantilever beam structure

3.2 懸臂梁解析解

圖4 懸臂梁模型Fig.4 Cantilever beam model

此問題的位移解析解和應(yīng)力解析解[19]為

u(x,y)=

(23)

v(x,y)=

(24)

(25)

由式(23)～式(25)計算得，此問題x方向的最大位移解析解 (ux)max=3.571×10-6m，y方向的最大位移解析解(uy)max=2.446 4×10-5m，法向應(yīng)力最大值解析解(σxx)max=5×104Pa.

3.3 有限元分析計算

采用ansys workbench對懸臂梁模型進行有限元分析，網(wǎng)格劃分結(jié)果如圖5 所示，計算結(jié)果如圖6 所示.

圖5 網(wǎng)格劃分結(jié)果Fig.5 Meshing result

(a) x方向位移變化云圖(b) y方向位移變化云圖(c) 法向應(yīng)力分布云圖圖6 有限元分析結(jié)果Fig.6 Analysis results of finite element

3.4 等幾何分析計算

計算u方向的節(jié)點向量，其中，n=5,p=3,m=9.

利用式(5)和式(6)計算d：

|Q1-Q0|=6,|Q2-Q1|=6,

|Q3-Q2|=6,|Q4-Q3|=6,

|Q5-Q4|=6,

d=30.

利用式(8)計算uj+p：

u0=…=u3=0,u6=…=u9=1，

計算v方向的節(jié)點向量，其中，n=2,p=2,m=5，直接根據(jù)式(8)可知

v0=…=v2=0,v3=…=v5=1，

所以，v=[0 0 0 1 1 1.

采用matlab編寫計算過程. 首先建立NURBS基礎(chǔ)數(shù)據(jù)庫，如表1 所示.

通過調(diào)用函數(shù)生成控制點網(wǎng)格和節(jié)點向量網(wǎng)格(其中經(jīng)過3次插值細分)，如圖7 和圖8 所示.

通過matlab編程初始化邊界，計算單元剛度矩陣并組裝，再通過Gauss積分計算節(jié)點力，邊界條件處理調(diào)用求解，最后的后處理結(jié)果如圖9 所示.

表1 NURBS基函數(shù)Tab.1 NURBS basis function

圖7 控制點網(wǎng)格

圖8 節(jié)點網(wǎng)格Fig.8 Node grid

(a) x方向位移變化云圖

(b) y方向位移變化云圖

(c) 法向應(yīng)力分布云圖圖9 等幾何分析結(jié)果Fig.9 Isogeometric analysis results

3.5 計算結(jié)果對比

對比有限元分析(FEA)和等幾何分析(IGA)計算的x方向和y方向的位移變化以及法向應(yīng)力(見圖6 和圖9)可知，兩種方法所得位移變化和法向應(yīng)力分布云圖基本一致，且最大值點位置相同，最大值與解析解誤差如表2 所示.

表2 FEA與IGA法的計算結(jié)果對比Tab.2 Comparison of calculation results in FEA and IGA

在相同的計算條件下，由matlab編程所產(chǎn)生的節(jié)點數(shù)量和網(wǎng)格數(shù)量等于或者大于由有限元分析生成的節(jié)點數(shù)量和網(wǎng)格數(shù)量. 為了保證兩種方法運算時間對比的公平性，采用控制變量，選擇與解析解精度相同的有限元解，保證了求解結(jié)果和解析解的一致性，從而使獲得的運算時間更具有代表性，同時，為了減小計算機硬件資源對運算效率分析的影響，特別是AVX指令集在CPU頻率不同時的運算效率差距很大，將頻率固定在3.0 GHz，內(nèi)存頻率固定在3 200 MHz，以減少自動睿頻對分析結(jié)果的影響. 分析時間定義為：有限元分析時間從網(wǎng)格劃分結(jié)束后開始到后處理結(jié)束(人為設(shè)置邊界條件和后處理方案不計入分析時間)；等幾何分析時間指全部分析、計算、后處理時間. 總時長定義為：網(wǎng)格劃分時間和分析時間的總和. 每個算例中，計算所耗費的時間都取10次的平均值. 運算時間統(tǒng)計結(jié)果如表3 所示.

表3 FEA與IGA法的運算時間對比Tab.3 Comparison of calculation time in FEA and IGA

由表2 可知，本文改進后的等幾何分析計算的x方向位移變化、y方向位移變化以及法向應(yīng)力變化比有限元法計算的精度更高. 由表3 可知，在保證計算精度的前提下，本文等幾何分析計算的總時長比有限元分析計算的總時長縮短了47.72%.

4 結(jié) 論

本文利用平均值法對等幾何分析的節(jié)點向量求解進行改進，建立了改進等幾何分析處理平面靜力學(xué)問題的結(jié)構(gòu)化流程. 該流程用平均值法代替原等距分布法，建立NURBS基礎(chǔ)信息庫，然后通過等幾何建模、單元細分、施加邊界條件、后處理計算四個步驟進行求解，并將其應(yīng)用于分析經(jīng)典靜力學(xué)問題. 通過等幾何分析法求解懸臂梁問題得到的x方向、y方向位移變化云圖與有限元分析一致，且最大值與解析解相比的誤差分別為0.84%和0.06%，法向應(yīng)力最大值的誤差比有限元的誤差減小了約5%，驗證了本文方法的可行性；在保持計算精度的情況下，改進后的等幾何分析的運算總時長比有限元分析縮短了47.72%，證明了本文方法的高效性.