李運(yùn)平 周來水 王志國 陸友太

南京航空航天大學(xué)江蘇省精密與微細(xì)制造技術(shù)重點(diǎn)實(shí)驗(yàn)室，南京，210016

0 引言

B－spline曲線曲面作為 STEP工業(yè)標(biāo)準(zhǔn)NURBS曲線曲面的特例(其權(quán)因子全為1)，是曲面造型領(lǐng)域中描述工業(yè)產(chǎn)品幾何外形的重要方法。由B－spline曲面的定義可知，通過改變節(jié)點(diǎn)矢量與控制頂點(diǎn)都可以實(shí)現(xiàn)曲面的變形［1-2］。然而，修改幾何參數(shù)的方式對(duì)于復(fù)雜形狀的曲面變形設(shè)計(jì)非常困難，且難以達(dá)到工程上的幾何光順性要求和設(shè)計(jì)精度。

為了避免這些缺陷，簡化設(shè)計(jì)過程，產(chǎn)生了基于物理的造型方法。Terzopoulos等［3］將基于能量模型的可變形曲線曲面造型技術(shù)引入到計(jì)算機(jī)圖形學(xué)中，模擬出旗幟在風(fēng)中飄擺及地毯飄落等動(dòng)態(tài)過程，為基于物理的變形技術(shù)奠定了基礎(chǔ)。近年來，基于物理的變形方法受到研究者的廣泛關(guān)注［4-5］。Celniker等［6］將有限元法引入到曲線曲面變形中，將曲面劃分為三角形板單元，在載荷作用與幾何約束下計(jì)算曲面變形。之后，基于有限元法的曲面變形技術(shù)得到了廣泛研究［7-10］。文獻(xiàn)［7-8］采用有限元法求解曲線曲面變形的思想，運(yùn)用3次均勻B－spline描述變形曲線曲面，給出了相應(yīng)的有限元網(wǎng)格模型，以控制頂點(diǎn)為求解未知量，建立總剛度矩陣方程，計(jì)算滿足約束的曲線曲面。成思源等［9］推導(dǎo)了 B－spline曲線曲面的可變形模型運(yùn)動(dòng)方程，并對(duì)其有限元求解方法進(jìn)行了研究，用特征點(diǎn)作用下可變形B－spline曲線曲面的變形驗(yàn)證了有限元法的有效性。采用有限元網(wǎng)格近似參數(shù)曲面會(huì)導(dǎo)致曲面自身的拓?fù)浒l(fā)生改變，因此，難以精確描述曲面間的光滑拼接關(guān)系。文獻(xiàn)［11］將初始幾何體整體或局部地嵌入到橫截面為圓形的梁結(jié)構(gòu)中，通過施加載荷使梁結(jié)構(gòu)產(chǎn)生形變，幾何體作為梁結(jié)構(gòu)的一部分在載荷作用下隨梁結(jié)構(gòu)一起產(chǎn)生變形。文獻(xiàn)［12］基于控制網(wǎng)格建立了剛架模型，實(shí)現(xiàn)了單張B－spline曲面形狀修改。文獻(xiàn)［11-12］避免了變形過程中曲面拓?fù)涞母淖?，但是僅適用于單張B－spline曲面。

筆者借鑒剛架受力變形的思想，提出了一種剛架模型驅(qū)動(dòng)多張B－spline曲面變形算法。通過“合并”曲面片的控制網(wǎng)格的手段創(chuàng)建剛架模型，將曲面變形問題轉(zhuǎn)化為剛架模型的變形。根據(jù)幾何約束、曲面片間光滑拼接條件和選取的優(yōu)化目標(biāo)，計(jì)算曲面片控制頂點(diǎn)的變化，實(shí)現(xiàn)多張B－spline曲面片變形。

1 多張B－spline曲面剛架模型的創(chuàng)建

1.1 剛架模型

空間剛架模型是由梁單元相互連接而成的結(jié)構(gòu)，用于驅(qū)動(dòng)曲面變形，如圖1所示。梁單元的節(jié)點(diǎn)稱為剛架的節(jié)點(diǎn)，節(jié)點(diǎn)的自由度為6，根據(jù)自由度的不同，可將節(jié)點(diǎn)分為自由節(jié)點(diǎn)、固支節(jié)點(diǎn)與鉸支節(jié)點(diǎn)。

圖1 剛架模型

1.2 B-spline曲面的數(shù)學(xué)描述

B－spline曲面由2個(gè)方向的控制頂點(diǎn)網(wǎng)格、2個(gè)節(jié)點(diǎn)矢量和單變量B－spline基函數(shù)的乘積定義。假設(shè)給定由(m+1)×(n+1)個(gè)控制頂點(diǎn)Sij(i=0，1，…，m;j=0，1，…，n)構(gòu)成的陣列組成一張控制網(wǎng)格，分別給定參數(shù)u(次數(shù)為k)與v(次數(shù)為l)和節(jié)點(diǎn)矢量U與V，就可以定義一張張量積為k×l次的B－spline曲面，其方程為

式中，Ni，k(u)是由節(jié)點(diǎn)矢量U定義的k次B－spline基函數(shù);Nj，l(v)是由節(jié)點(diǎn)矢量V定義的l次B－spline基函數(shù);N(u，v)為B－spline基函數(shù)乘積矢量;S為控制頂點(diǎn)矢量。

1.3 單張B-spline曲面的剛架模型

為了更好地反映控制頂點(diǎn)位置的改變對(duì)曲面形狀的影響，用曲面的控制頂點(diǎn)定義剛架模型的節(jié)點(diǎn)，依據(jù)控制網(wǎng)格的幾何拓?fù)潢P(guān)系生成相應(yīng)的剛架模型。圖2所示為基于控制網(wǎng)格所創(chuàng)建的單張B－spline曲面的剛架模型，圓球?yàn)楣?jié)點(diǎn)。

1.4 多張B-spline曲面的剛架模型

圖2 單張B-spline曲面的剛架模型

工業(yè)產(chǎn)品的幾何外形通常由4邊域B－spline曲面片拼接而成，且曲面片之間一般維持G1連續(xù)，即曲面片的每個(gè)公共頂點(diǎn)與位于兩側(cè)的各一個(gè)相鄰非公共頂點(diǎn)共線，如圖3所示，其中，P、P'為公共邊界處的重合點(diǎn)，P1、P2分別為兩側(cè)對(duì)應(yīng)控制頂點(diǎn)。采取去掉重復(fù)控制頂點(diǎn)的方法合并控制網(wǎng)格，建立多張曲面的剛架模型，如圖4所示。圖4中，黑色球?yàn)楣讨Ч?jié)點(diǎn)，灰色球?yàn)樽杂晒?jié)點(diǎn)。由圖4可知，剛架模型的節(jié)點(diǎn)和合并后的控制網(wǎng)格頂點(diǎn)是一一對(duì)應(yīng)的，故節(jié)點(diǎn)的位移量等價(jià)于控制網(wǎng)格頂點(diǎn)的位移量。

圖3 G1連續(xù)條件下多張B-spline曲面和控制網(wǎng)格

圖4 多張B-spline曲面剛架模型

2 多張B－spline曲面約束變形計(jì)算

2.1 基于有限元法的剛架模型節(jié)點(diǎn)位移計(jì)算

剛架模型在承受載荷的情況下，節(jié)點(diǎn)會(huì)產(chǎn)生相應(yīng)的位移，根據(jù)線彈性有限元方法可知，載荷和節(jié)點(diǎn)位移滿足如下方程:

式中，K為整體剛度矩陣;Δ為節(jié)點(diǎn)位移列陣;P為節(jié)點(diǎn)載荷列陣。

在曲面造型過程中，為滿足工程的設(shè)計(jì)要求，將曲面設(shè)計(jì)要求抽象為幾何約束(點(diǎn)約束、點(diǎn)和法矢約束、等參線約束等)，在幾何約束下對(duì)曲面進(jìn)行變形。根據(jù)多曲面剛架模型節(jié)點(diǎn)和控制網(wǎng)格頂點(diǎn)的一一對(duì)應(yīng)關(guān)系，將曲面的約束變形問題轉(zhuǎn)化為剛架模型變形。以剛架模型節(jié)點(diǎn)位移為變量，構(gòu)建了載荷最小與節(jié)點(diǎn)位移量最小的優(yōu)化目標(biāo)，并建立統(tǒng)一的表達(dá)式，用于組合使用兩個(gè)最小化模型來滿足不同的變形要求。多曲面約束變形問題可描述為以剛架節(jié)點(diǎn)位移Δ為變量的函數(shù):

式中，‖KΔ‖2為節(jié)點(diǎn)載荷;‖Δ‖2為節(jié)點(diǎn)位移量;α、β為權(quán)因子，用于調(diào)整目標(biāo)函數(shù)中每個(gè)優(yōu)化目標(biāo)在變形中所占的比重;ArΔ=br為第r個(gè)約束的約束方程;Ar為約束方程的系數(shù)矩陣;br為列陣;cn為約束的個(gè)數(shù)。

由于優(yōu)化問題(式(3))中的目標(biāo)函數(shù)為Δ的二次函數(shù)，約束方程為線性方程組，因此上述約束優(yōu)化問題為標(biāo)準(zhǔn)的等式二次規(guī)劃問題，可采用Lagrange乘子法和罰函數(shù)法進(jìn)行求解。按照罰函數(shù)法將式(3)改寫成

式中，I為單位矩陣。

由式(4)可以看出，罰函數(shù)法將約束優(yōu)化問題轉(zhuǎn)化為線性方程組的求解問題。本文先采用超松弛迭代預(yù)處理共軛梯度法(SSOR－PCG)求解節(jié)點(diǎn)位移Δ，然后根據(jù)節(jié)點(diǎn)和控制頂點(diǎn)的對(duì)應(yīng)關(guān)系，計(jì)算曲面片控制頂點(diǎn)的新位置，完成多張曲面的約束變形。

由于總剛度矩陣K具有稀疏、對(duì)稱等特性，本文采用一維變帶寬壓縮存儲(chǔ)方法存儲(chǔ)矩陣，并選取SSOR－PCG迭代求解方程組。將式(4)簡記為FΔ=u，D為F的對(duì)角陣，L為F的嚴(yán)格下三角陣，求解方程組偽代碼如下:

2.2 曲面片光滑拼接約束

為保證變形過程中曲面片之間的光滑拼接，將曲面片之間的光滑拼接關(guān)系作為約束條件。在圖3 中，曲面S1、S2為G1連續(xù)曲面，即P1、P、P2位于一條直線上，則?0≤t≤1，滿足:

將式(5)改寫為以節(jié)點(diǎn)位移Δ為參數(shù)的線性方程組:

式中，Ai'j'(0≤ i'≤ m'，0≤ j'≤ n')為 3×6維矩陣，m'+1、n'+1分別為合并后的整合控制網(wǎng)格U、V向控制頂點(diǎn)數(shù)目。

記與 P1、P、P2對(duì)應(yīng)的節(jié)點(diǎn)索引分別為 i'j'、k'l'、r's'(0≤k'≤m';0≤r'≤m';0≤l'≤n';0≤ s'≤ n')，則矩陣 A 除 Ai'j'、Ak'l'、Ar's'外，其余元素全為零矩陣。Ai'j'、Ak'l'、Ar's'分別為

2.3 點(diǎn)約束變形

假定S0是曲面片Pw(u，v)=Nw(u，v)Sw上參數(shù)為(u0，v0)的點(diǎn)，T0為變形后曲面片P'w(u，v)上對(duì)應(yīng)的目標(biāo)點(diǎn)。 則變形后的曲面片 P'w(u，v)滿足:

式(7)可以改寫為

將上述方程轉(zhuǎn)化為以節(jié)點(diǎn)位移為變量的約束方程AΔ=b。其中，A的形式與式(6)中的系數(shù)矩陣一致，b= ΔP0=［Δx Δy Δz］T。

初始各曲面片的每一個(gè)控制頂點(diǎn)對(duì)應(yīng)唯一的節(jié)點(diǎn)，元素Ai'j'取值主要有兩種情況:

(1)若控制頂點(diǎn)Swij和節(jié)點(diǎn)nodei'j'對(duì)應(yīng)，記Nwij=Nwi(u)Nwj(v)，則元素Ai'j'可以表示為

(2)若節(jié)點(diǎn)nodei'j'和任一控制頂點(diǎn)Swij未重合，則 Ai'j'=0。

在多點(diǎn)約束的情況下，點(diǎn)約束可能會(huì)分布在不同的曲面片上。通過上述討論分別計(jì)算將每個(gè)點(diǎn)約束轉(zhuǎn)化為約束矩陣ArΔ=br。

圖5所示為G1連續(xù)條件下4張B－spline曲面點(diǎn)約束變形結(jié)果。圖5a所示為G1連續(xù)條件下4張B－spline曲面及點(diǎn)約束，圖5b、圖5c所示分別為給定不同權(quán)因子的曲面變形結(jié)果。由圖5b、圖5c可以看出，變形后的曲面片都滿足給定的點(diǎn)約束條件，但隨著β的逐漸增大(節(jié)點(diǎn)位移最小優(yōu)化目標(biāo)逐漸占較大比重)，曲面將逐漸由整體光順模型(圖5b)向局部尖銳模型(圖5c)過渡。因此，調(diào)節(jié)權(quán)因子可以獲得形狀各異的變形結(jié)果，為設(shè)計(jì)需求提供充足的素材。

圖5 G1連續(xù)條件下4張B-spline曲面點(diǎn)約束變形

2.4 點(diǎn)與法矢約束變形

假設(shè)變形后的曲面片P'w(u，v)在參數(shù)點(diǎn)(u0，v0)處能精確通過點(diǎn)T0，且曲面在(u0，v0)處的法矢 n=(nx，ny，nz)，即曲面片 P'w(u，v)滿足:

式(8)可以根據(jù)點(diǎn)約束下多張B－spline曲面變形中約束系數(shù)矩陣的求法轉(zhuǎn)換成AΔ=b的形式;式(9)、式(10)表達(dá)形式一致，其解決方法也類似。以式(9)為例，其約束方程式可以表示為

式中，Bi'j'為1×6維矩陣。

對(duì)約束矩陣元素 Bi'j'同樣需要進(jìn)行分情況討論:

(1)若控制頂點(diǎn)Swij和節(jié)點(diǎn)nodei'j'對(duì)應(yīng)，則元素Bi'j'為

(2)若節(jié)點(diǎn)nodei'j'和任一控制頂點(diǎn)Swij未重合，則 Bi'j'=0。

同理，式(10)可以表示成CΔ=d的形式。

圖6所示為G1連續(xù)條件下基于點(diǎn)和法矢約束的3張B－spline曲面變形。圖6a所示為初始曲面及點(diǎn)與法矢約束，圖6b、圖6c所示為不同權(quán)因子值下B－spline曲面的變形效果。

圖6 G1連續(xù)條件下3張B-spline曲面點(diǎn)與法矢約束變形

3 變形實(shí)例

基于剛架模型的多曲面變形方法已通過編程實(shí)現(xiàn)，實(shí)例中梁單元材料的彈性模量E=1kPa、泊松比μ=0.4，通過交互添加幾何約束使物體達(dá)到滿足特定需求的變形結(jié)果。圖7所示為自由曲面的點(diǎn)約束變形，圖7a中的自由曲面由3張B－spline曲面片通過G1拼接而成，黑色的點(diǎn)為初始曲面片上的點(diǎn)，黑色三角形為其相應(yīng)的目標(biāo)點(diǎn)，最終的變形結(jié)果如圖7b所示，實(shí)現(xiàn)了曲面精確變形到點(diǎn)約束位置的要求。

圖7 自由曲面點(diǎn)約束變形

和傳統(tǒng)的曲面變形方法相比，文中剛架模型的變形計(jì)算采用有限元的直接剛度法，由于總剛度矩陣是正定的，理論上一定可解。因此運(yùn)用剛架變形原理進(jìn)行復(fù)雜曲面的變形設(shè)計(jì)時(shí)，避免了拼接約束系統(tǒng)的可解性問題，本文方法可確保曲面片之間在變形過程中保持G1連續(xù)。表1所示為圖5～圖7曲面變形前后曲面片之間G1連續(xù)性條件分析數(shù)據(jù)，平均正弦為公共邊界點(diǎn)和兩側(cè)對(duì)應(yīng)的控制頂點(diǎn)連線的矢量夾角正弦的平均值。由表1可以看出，變形前后的平均正弦值基本保持一致，即曲面片之間的光滑拼接狀態(tài)在變形中得到了有效的保持。

表1 多曲面變形前后G1連續(xù)性分析

運(yùn)用本文算法計(jì)算多曲面變形時(shí)，要求每張曲面為4邊域B－spline曲面，且兩曲面片公共邊界處及兩側(cè)對(duì)應(yīng)的控制頂點(diǎn)在一條直線上。在約束變形計(jì)算中，約束條件過多會(huì)導(dǎo)致剛架模型處于過約束狀態(tài)，進(jìn)而無法生成滿足指定約束的變形結(jié)果。圖8所示為G1連續(xù)條件下基于點(diǎn)約束的2張B－spline曲面變形。其中，圖8a所示為初始曲面、曲面的剛架模型及點(diǎn)約束，圖8b所示為基于載荷最小模型的下B－spline曲面的變形結(jié)果。由圖8b可以看出，變形結(jié)果僅通過1個(gè)目標(biāo)點(diǎn)。

圖8 G1連續(xù)條件下2張B-spline曲面點(diǎn)約束變形

4 結(jié)語

結(jié)合剛架模型受力變形原理，提出了一種多張B－Spling曲面變形設(shè)計(jì)算法。給出了G1連續(xù)條件下多張B－spline曲面剛架模型的建立，以及曲面片之間光滑拼接條件的處理方法。算法實(shí)現(xiàn)了多張B－spline曲面在幾何約束下的變形協(xié)調(diào)設(shè)計(jì)，而且還能使曲面片之間保持變形前的光滑拼接狀態(tài)。

為增加變形算法的實(shí)用性可將本文方法擴(kuò)展應(yīng)用于NURBS表示的曲面變形設(shè)計(jì);文中只提供了節(jié)點(diǎn)載荷最小和節(jié)點(diǎn)位移最小兩個(gè)優(yōu)化目標(biāo)，后續(xù)工作可以進(jìn)一步組合更多不同意義的優(yōu)化目標(biāo)，以便達(dá)到理想的曲面變形效果。

［1］Piegl L.Modifying the Shape of Rational B－splines.Part 1:Surfaces［J］.Computer－Aided Design，1989，21(9):538-546.

［2］李亞娟，汪國昭.改變B樣條曲面與NURBS曲面的節(jié)點(diǎn)［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2005，17(5):986-989.Li Yajun，Wang Guozhao.On Knot Modifications of B－spline or NURBS Surface［J］.Journal of Computer－Aided Design ＆ Computer Graphics，2005，17(5):986-989.

［3］Terzopoulos D，Platt J，Barr A.Elastically Deformable Models［J］.Computer Graphics，1987，21(4):205-214.

［4］Wang Z G，Zhou L S，Wang X P.Direct Manipulation of B－spline Surfaces［J］.Chinese Journal of Mechanical Engineering，2005，18(1):103-108.

［5］王青，柯映林，李江雄.基于力密度方法的NURBS曲線和曲面變形框架［J］.機(jī)械工程學(xué)報(bào)，2007，43(3):135-142.Wang Qing，Ke Yinglin，Li Jiangxiong.Shape Modification of NURBS Curves and Surfaces Based on Force Density Method［J］.Journal of Mechanical Engineering，2007，43(3):135-142.

［6］Celniker G，Gossard D.Deformable Curve and Surface Finite Elments for Free－form Shape Design［J］.Computer Graphics，1991，25(4):257-266.

［7］Guan Z D，Jing L，Ning T，et al.Study and Application of Physics－based Deformable Curves and Surfaces［J］.Computers ＆ Graphics，1997，21(3):305-313.

［8］成思源，張湘?zhèn)?可變形B樣條曲線曲面模型及其有限元求解［J］.重慶大學(xué)學(xué)報(bào)，2003，26(2):75-77.Cheng Siyuan，Zhang Xiangwei.Deformable B－spline Curves and Surfaces Model and Its Finite Element Solution［J］.Chinese Journal of Chongqing University，2003，26(2):75-77.

［9］經(jīng)玲，席平，唐榮錫.有限元方法在變形曲線曲面造型中的應(yīng)用［J］.計(jì)算機(jī)學(xué)報(bào)，1998，21(3):245-251.Jing Ling，Xi Ping，Tang Rrongxi.Application of Finite Element Method in Deformable Curve and Surface［J］.Journal of Computers，1998，21(3):245-251.

［10］程仙國，劉偉軍，卞宏友.基于外載荷的B樣條曲面變形［J］.機(jī)械工程學(xué)報(bào)，2011，47(9):146-150.Cheng Xianguo，Liu Wweijun，Bian Hongyou.Shape Modification of B－spline Surface via External Loads［J］.Journal of Mechanical Engineering，2011，47(9):146-150.

［11］Wang Z G，Wang X P，Bao Y D，et al.Shape Modification by Beam Model in FEM［J］.Chinese Journal of Aeronautics，2010，23(2):246-251.

［12］Yu C S，Zhou L S，Wang Z G，et al.Shape Modification of B－spline Surfaces Based on Rigid Frame via Geometric Constraint Optimization［C］//Proceedings of the Conference on Intelligent Human－machine Systems and Cybernetices.Hangzhou，2011:120-124.