于麗潔 孫紅巖

摘 要：通過對(duì)Waterbomb鋪嵌模型的結(jié)構(gòu)分析，提出雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)并對(duì)其進(jìn)行模擬。根據(jù)與棱相鄰的三角形面片是否對(duì)稱，將Waterbomb鋪嵌模型的非邊界棱（邊）歸為兩類，建立兩類彎曲彈簧，設(shè)置不同的彎曲彈性系數(shù)，該質(zhì)點(diǎn)彈簧系統(tǒng)被稱為雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)。然后運(yùn)用隱式歐拉法和塊坐標(biāo)下降法求解最小化問題，得到下一時(shí)刻物體運(yùn)動(dòng)變形的速度和位置，實(shí)現(xiàn)對(duì)Waterbomb鋪嵌模型快速、穩(wěn)定和準(zhǔn)確的模擬。實(shí)驗(yàn)證明，在誤差允許范圍內(nèi)，該算法速度快于牛頓法。相比單彎曲（基本）系統(tǒng)，雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)提升了模擬穩(wěn)定性、真實(shí)性，并減小了誤差。

關(guān)鍵詞：質(zhì)點(diǎn)彈簧系統(tǒng);彈性變形模擬;薄殼;Waterbomb折紙鋪嵌;雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)

DOI：10. 11907/rjdk. 181788

中圖分類號(hào)：TP317.4文獻(xiàn)標(biāo)識(shí)碼：A文章編號(hào)：1672-7800（2019）001-0186-06

Abstract：According to the analysis of waterbomb tessellation model structure， bi-bending mass-spring system is proposed to simulate waterbomb tessellation. According to the symmetry of the triangular facets adjacent to the edge， the non boundary edges of the waterbomb tessellation model are classified into two classes， and two kinds of bending springs are established with different bending elasticity coefficients. This mass-spring system is called the bi-bending mass-spring system. Then， we use the implicit Euler method and solve the minimization problem by block coordinate descent method （also known as alternating optimization） to get the velocity and position of the movement and deformation of the object at the next time. This paper realizes the fast， stable and accurate simulation of waterbomb tessellation model. Within the allowable error range， the algorithm is faster than Newtons method. Compared with the single bending （basic） system， the bi-bending mass-spring system further improves the stability and realism of the simulation and reduces the error.

0 引言

基于物理的模擬方法符合客觀物體運(yùn)動(dòng)和變形規(guī)律，能夠很好地保持運(yùn)動(dòng)和變形的真實(shí)性。在基于物理的模擬中，模型的表示方法包括質(zhì)點(diǎn)彈簧模型、有限元模型、無網(wǎng)格表達(dá)等。質(zhì)點(diǎn)彈簧模型由一系列質(zhì)點(diǎn)和連接各個(gè)質(zhì)點(diǎn)的彈簧構(gòu)成，是一種最簡(jiǎn)單的彈性體模擬方法，在圖形學(xué)中得到了快速發(fā)展。不同種類的彈簧具有不同特性。

在圖形學(xué)的物理模擬中，薄殼容易彎曲、寬度與厚度之比大于100，在未變形狀態(tài)下是彎曲的結(jié)構(gòu)[1]。而與之相對(duì)的薄板在未受力狀態(tài)下是平整的。常見的薄殼結(jié)構(gòu)有轎車外殼、草帽、紙張、雞蛋殼、樹葉等。

紙張本身是剛性的薄殼結(jié)構(gòu)。當(dāng)紙張被折疊，形成具有負(fù)泊松比的鋪嵌模型時(shí)，紙張整體便具有了彈性，成為柔性曲面，可形成各種形狀。折紙（Origami）既是一種古老的藝術(shù)形式，又是一種很好的娛樂方式。人們可以通過折紙藝術(shù)將具有水密性的硬質(zhì)材料如紙張、金屬等通過彎曲或折疊形成復(fù)雜的三維形狀。

拼嵌模式包括waterbomb tessellation、Miura-ori、Yoshimura pattern、Resch's pattern、Curved folding、Huffman's Origami Tessellations等。Waterbomb樣模鋪嵌或水雷鋪嵌是應(yīng)用最廣泛的折紙圖案之一。它是基于三角形可平面折疊的圖案模式，具有周期對(duì)稱，泊松比為負(fù)，屬于具有各向異性的拉漲材料[2]。Waterbomb樣模鋪嵌的負(fù)泊松比特性通過剛性的完整紙張折疊得到，無需通過剪切實(shí)現(xiàn)。Waterbomb樣模分為六折痕樣模和八折痕樣模[3]，兩種樣?？尚纬刹煌匿伹赌Ｐ汀Ｓ捎谇婕扔邪键c(diǎn)也有凸點(diǎn)，因此Waterbomb鋪嵌具有足夠的柔性，可形成多種折紙圖案[4]。

1 相關(guān)工作

質(zhì)點(diǎn)彈簧模型原理簡(jiǎn)單，易于操作，在更加注重模擬速度和穩(wěn)定性、而不十分注重模擬精度的實(shí)時(shí)模擬中應(yīng)用廣泛。2009年，Lu等[5]用雙層質(zhì)點(diǎn)彈簧系統(tǒng)對(duì)植物葉片進(jìn)行模擬。2013年，劉雪梅等[6]基于粘彈性力學(xué)模型耦合無網(wǎng)格伽遼金與質(zhì)點(diǎn)彈簧對(duì)軟組織進(jìn)行形變仿真，在大變形和拓?fù)涓淖儏^(qū)域使用無網(wǎng)格伽遼金法，其它區(qū)域使用質(zhì)點(diǎn)彈簧法，兩個(gè)區(qū)域之間建立過渡單元。高精度及與拓?fù)錈o關(guān)的無網(wǎng)格伽遼金法與高效簡(jiǎn)單的質(zhì)點(diǎn)彈簧組合的方式有效解決了精度與效率平衡的問題，實(shí)現(xiàn)了高精度實(shí)時(shí)形變仿真。2013年，陳瑛等[7]根據(jù)織物拉伸形變過程中受力與形變程度的非線性關(guān)系特征，提出分段線性化彈性系數(shù)的改進(jìn)質(zhì)點(diǎn)彈簧模型，解決了織物模擬中超彈性的問題，提高了模擬真實(shí)感。2013年，Liu等[8]將隱式歐拉方法轉(zhuǎn)換為能量最小化問題，并引入彈簧方向作為輔助未知變量，用塊下降法對(duì)最小化問題進(jìn)行求解，實(shí)現(xiàn)了對(duì)質(zhì)點(diǎn)彈簧系統(tǒng)的快速實(shí)時(shí)模擬。2015年，王崴等[9]對(duì)繩索形變進(jìn)行模擬，將傳統(tǒng)正方形質(zhì)點(diǎn)彈簧網(wǎng)格進(jìn)行改進(jìn)，提出一種蜂窩狀的彈簧質(zhì)點(diǎn)模型，利用Verlet-梯形預(yù)測(cè)-校正法進(jìn)行數(shù)值積分以求得下一時(shí)刻速度和位置。實(shí)驗(yàn)結(jié)果表明，蜂窩狀質(zhì)點(diǎn)彈簧模型有效，算法提高了計(jì)算效率和精確度。2017年，王舒珍等[10]在傳統(tǒng)四邊形模型的基礎(chǔ)上對(duì)質(zhì)點(diǎn)彈簧模型進(jìn)行改進(jìn)，去除剪切彈簧，以結(jié)構(gòu)彈簧代替，對(duì)軟組織形變進(jìn)行模擬。

自1987年Terzopoulos等[11]對(duì)彈性變形物體進(jìn)行模擬后，基于物體的模擬方法被成功運(yùn)用到立體、薄殼、布料、流體的變形模擬中。在薄殼模擬中，2003年，Grinspun等[1]提出了一種簡(jiǎn)單的離散薄殼模型，該模型無需重復(fù)之前模型復(fù)雜的連續(xù)介質(zhì)力學(xué)公式和算法。該方法主要關(guān)注不可拉伸的薄殼，即薄殼的變形主要為等距變換。彎曲能量用相鄰面片二面角的變化進(jìn)行描述。在此基礎(chǔ)上，2007年Garg等[12]提出實(shí)現(xiàn)更簡(jiǎn)單、計(jì)算更高效的基于折葉（hinge-based）的物理模型，對(duì)易彎曲或接近剛性的、各向同性或各向異性的薄板或薄殼進(jìn)行模擬。一個(gè)彎曲折頁(yè)由兩個(gè)相鄰三角形面片形成。由于折頁(yè)的方法簡(jiǎn)單，計(jì)算復(fù)雜度低，在網(wǎng)格模型模擬中大受歡迎。2014年，Bouaziz等[13]提出投影動(dòng)力學(xué)，能夠支持一系列幾何約束，在同一框架中對(duì)立體、布料、薄殼進(jìn)行基于實(shí)例的模擬。其中，薄殼的彎曲能量用變形狀態(tài)下和未變形狀態(tài)下平均曲率的變化（差的平方）描述。對(duì)于等距變換，該彎曲能量可以用附加的（auxiliary）旋轉(zhuǎn)矩陣重寫為拉普拉斯-貝爾特粒米算子，然后用分段線性的帽子基函數(shù)（hat basis）將該連續(xù)彎曲能量進(jìn)行離散化。2016年，Zhang等[14]提出用交互工具對(duì)3D中空打印模型變形進(jìn)行設(shè)計(jì)。根據(jù)模型各部分彎曲彈性的不同，中空模型各部分被設(shè)計(jì)成不同厚度。2017年，Liu等[15]將投影動(dòng)力學(xué)解釋為擬-牛頓法，并用L-BFGS進(jìn)行加速，能夠支持多種不同超彈材料的模擬，包括新虎克彈性材料和基于樣條的材料。

在折紙領(lǐng)域方面，2003年Grinspun等[1]提出可通過修改標(biāo)準(zhǔn)布料模擬器快速實(shí)現(xiàn)離散薄殼模型。2006年，Burgoon等[16]以該薄殼模型為基礎(chǔ)，考慮了紙張彈性，提出將紙張折疊成為簡(jiǎn)單模型的交互界面。該方法沒有重點(diǎn)處理折痕多的紙張間碰撞問題，折出了簡(jiǎn)單的折紙模型，如5個(gè)折痕的小狗模型。2010年，Tomohiro Tachi在文獻(xiàn)[4]中提出設(shè)計(jì)三維自由曲面的方法和Freeform Origami 系統(tǒng)，可對(duì)給定的圖案進(jìn)行自由編輯;在文獻(xiàn)[17]中，Tomohiro Tachi提出一種折出任意三維形狀的方法Origamizer Method， 可將材料折疊成給定多面體的折痕圖。2013年，Tomohiro Tachi[18]研究了由給定多面體曲面設(shè)計(jì)Resch式折紙鋪嵌自由形狀的方法;在文獻(xiàn)[19]中，Tomohiro Tachi提出交互式模擬和設(shè)計(jì)折紙的計(jì)算方法，其中折紙形狀由每個(gè)樣模（圖元，基本型）彎曲情況決定。2013年，Davis等[20]對(duì)David Huffman的折紙鋪嵌進(jìn)行了重建。

2 Waterbomb鋪嵌上的雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)

2.1 模擬對(duì)象

本文對(duì)圖1（a）所示的折紙鋪嵌圖案類型進(jìn)行模擬，其模擬模型來自Tomohiro Tachi的Freefrom Origami系統(tǒng)[4]，見圖1（b）。

2.2 雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)

Waterbomb鋪嵌模型的質(zhì)點(diǎn)彈簧模型主要包括結(jié)構(gòu)彈簧和彎曲彈簧。

將紙張沿折線向上折疊，稱為峰折或山折（moutain fold），向下折疊稱為谷折（valley fold），向上向下折兩次稱為中性折（genderless fold），帶有折痕模式的紙張平面圖稱為折痕圖或CP圖（crease pattern），最基本的折痕圖稱為waterbomb 樣模（base）、圖元或基本型，將一組waterbomb 樣模在三維空間按照一定規(guī)律重復(fù)拼嵌便組成waterbomb 鋪嵌或拼嵌（tessellation）。Waterbomb樣模邊界邊、峰折折痕和谷折折痕統(tǒng)稱為邊（或棱）。

六折痕Waterbomb樣模解析如圖2所示。

2.2.1 結(jié)構(gòu)彈簧

在所有棱（包括Waterbomb樣模邊界邊、峰折折痕和谷折折痕）上設(shè)置結(jié)構(gòu)彈簧，如圖3中直線所示。

2.2.2 彎曲彈簧

（1）單彎曲彈簧。將所有彎曲彈簧均使用同一彈性系數(shù)的質(zhì)點(diǎn)彈簧系統(tǒng)稱為單彎曲質(zhì)點(diǎn)彈簧系統(tǒng)。以單個(gè)Waterbomb樣模為例，對(duì)如何建立其鋪嵌模型質(zhì)點(diǎn)彈簧系統(tǒng)中的彎曲彈簧進(jìn)行說明。單個(gè)Waterbomb樣模的質(zhì)點(diǎn)彈簧模型如圖3所示。

對(duì)于薄殼結(jié)構(gòu)來說，彎曲變換是一個(gè)重要變換。根據(jù)對(duì)圖2的分析，本文以棱為中心（面向棱）構(gòu)建Waterbomb彎曲彈簧，如圖3中的虛線點(diǎn)狀曲線所示。例如，與邊OB相鄰的三角形分別為三角形OBE和三角形OBD，其上與邊OB相鄰的頂點(diǎn)分別為點(diǎn)E和點(diǎn)D。連接ED，建立在邊OB上的彎曲彈簧。

（2）雙彎曲彈簧。由于Waterbomb鋪嵌模型變形時(shí)具有各向異性特性，本文在系統(tǒng)中設(shè)定具有不同彎曲彈性系數(shù)的彎曲彈簧，即雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)。仍然以單個(gè)Waterbomb樣模為例說明Waterbomb折紙鋪嵌模型雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)的構(gòu)建（見圖4）。通過對(duì)Waterbomb樣模鋪嵌的結(jié)構(gòu)進(jìn)行分析，發(fā)現(xiàn)根據(jù)棱兩邊三角形面片結(jié)構(gòu)是否對(duì)稱，可以將Waterbomb鋪嵌模型非邊界棱分為3類：①與棱相鄰的兩個(gè)三角形面片關(guān)于棱是軸對(duì)稱的，該類棱包括：峰折折痕（如邊OE）和Waterbomb樣模中邊界較長(zhǎng)的邊（如邊BD）;②與棱相鄰的兩個(gè)三角形面片關(guān)于棱是中心對(duì)稱的，該類棱包括Waterbomb樣模中邊界較短的邊（如邊EB）;③與棱相鄰的兩個(gè)三角形面片關(guān)于棱不具有對(duì)稱性，該類棱包括谷折折痕（如邊OB）。

以棱為中心，根據(jù)不同類型的棱（面向棱）構(gòu)建Waterbomb的彎曲彈簧（見圖4）。

根據(jù)以上分析，本文將Waterbomb折紙樣模鋪嵌模型的棱（或邊）分為3類。本文將與棱兩邊相鄰的三角形面片對(duì)稱（包括軸對(duì)稱和中心對(duì)稱）的第1、第2類棱統(tǒng)稱為Ⅰ類棱，Ⅰ類棱包括Waterbomb樣模所有邊界邊和“峰折”折痕;將與棱兩邊相鄰的三角形面片不具有對(duì)稱性的第3類棱統(tǒng)稱為Ⅱ類棱。

本文在Ⅰ類棱和Ⅱ類棱上分別建立彎曲彈性系數(shù)不同的彎曲彈簧。將Ⅰ類棱上的彎曲彈簧稱為峰彎曲彈簧（或者相同彎曲彈簧），Ⅱ類棱上的彎曲彈簧稱為谷彎曲彈簧（或者不同彎曲彈簧）。

圖4展示了一個(gè)Waterbomb樣模的雙彎曲質(zhì)點(diǎn)彈簧構(gòu)建。以棱OE為例：與棱OE相鄰的三角形面片為三角形OEB和三角形OEA，從而與棱相對(duì)頂點(diǎn)為點(diǎn)A和點(diǎn)B。因此對(duì)于棱OE，本文連接鄰接點(diǎn)點(diǎn)A和點(diǎn)B建立相同彎曲彈簧AB。同理，在棱OB上連接鄰接點(diǎn)點(diǎn)E和點(diǎn)D，建立谷彎曲彈簧ED或不同彎曲彈簧ED。

3 基本概念與受力分析

設(shè)所模擬物體質(zhì)點(diǎn)數(shù)為[m]，經(jīng)歷的離散時(shí)刻依次為[t1，t2，t3，？]，時(shí)間步長(zhǎng)用[h]表示。物體在時(shí)刻[tn]的位置坐標(biāo)為[qn∈R3m]。系統(tǒng)所受的力用非線性函數(shù)表示為[f：R3m][→R3m]，則[tn]時(shí)刻系統(tǒng)所受的力為[f（qn）]。[f=-？E]，其中[E]為包括內(nèi)部力和外部力的勢(shì)能函數(shù)，[E：R3m→R]。

設(shè)[M∈R3m×3m]為對(duì)角質(zhì)量矩陣，[vn]為[tn]時(shí)刻的速度。

物體所受的力包括內(nèi)力和外力。內(nèi)力是指質(zhì)點(diǎn)間的相互作用力，主要有彈性變形產(chǎn)生的力。本文中內(nèi)力主要是彈簧彈性力。根據(jù)胡克定律，彈簧彈性勢(shì)能為：[Es=12kp1-p2-r2]。其中，[k]是彈性系數(shù)，[r]為彈簧原長(zhǎng)（[r0]），[p1，p2∈R3]為彈簧的兩個(gè)端點(diǎn)。外力指物體受到的外部作用力，包括物體本身的重力、外界施加的力、碰撞反應(yīng)力、風(fēng)力等。本文暫時(shí)只考慮物體所受的重力。外力[fext=Mg]，其中[g]為重力常數(shù)。

4 數(shù)值方法

本文采用文獻(xiàn)[8]的隱式歐拉法和數(shù)值求解方法。

4.1 隱式歐拉法

4.2 數(shù)值求解

最小化問題可用塊坐標(biāo)下降法（又稱交替優(yōu)化法）求解，以下稱為local-global方法。首先，為[x]指定一個(gè)初始猜測(cè)值[y]。在局部過程中，固定[x]計(jì)算最優(yōu)的[d];在全局過程中，固定[d]計(jì)算最優(yōu)的[x]。一直重復(fù)該過程直到收斂或達(dá)到最大迭代次數(shù)。

5 實(shí)驗(yàn)結(jié)果與分析

本部分首先從誤差、穩(wěn)定性及模擬速度方面，對(duì)單彎曲（基本）質(zhì)點(diǎn)彈簧模型的實(shí)驗(yàn)結(jié)果進(jìn)行分析和討論，確定迭代次數(shù)和彎曲彈性系數(shù)的范圍，然后對(duì)雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)的實(shí)驗(yàn)結(jié)果進(jìn)行分析討論。

本文實(shí)驗(yàn)環(huán)境為64位Windows10 操作系統(tǒng)，Intel（R） Core（TM） i56300HQ CPU@ 2.30GHz 。本文采用的實(shí)驗(yàn)?zāi)Ｐ腿鐖D1所示，該模型質(zhì)點(diǎn)個(gè)數(shù)為67，邊數(shù)為170，三角形面片數(shù)為104。在實(shí)驗(yàn)中，設(shè)定時(shí)間步長(zhǎng)為[h=1/30s]，質(zhì)點(diǎn)總質(zhì)量為10g，結(jié)構(gòu)彈簧拉伸彈性系數(shù)為80，重力常量為10kg/N，阻尼系數(shù)為0.001。

5.1 參數(shù)范圍

單彎曲質(zhì)點(diǎn)彈簧模型即基本質(zhì)點(diǎn)彈簧模型。本文首先用單彎曲質(zhì)點(diǎn)彈簧模型確定模擬參數(shù)大致范圍。

5.1.1 實(shí)驗(yàn)效果

從圖5可以看出，與沒有彎曲彈簧相比，有彎曲彈簧的質(zhì)點(diǎn)彈簧系統(tǒng)模擬的Waterbomb鋪嵌模型折痕不會(huì)隨時(shí)間慢慢展平，很好地保持了模型變形的細(xì)節(jié)特征，由此說明彎曲彈簧的重要性。雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)對(duì)彎曲彈簧進(jìn)行細(xì)分，提高了模擬精確性。

5.1.2 誤差與穩(wěn)定性分析

其中，[x0]為初始的預(yù)測(cè)值（[x0=y]），[x*]為采用牛頓法迭代到收斂時(shí)的值，[xi]表示第[i]次迭代時(shí)的位置。從定義可以看出，將由牛頓法得到的值作為基準(zhǔn)（“準(zhǔn)確值”）。相對(duì)誤差主要描述了算法的收斂性，“絕對(duì)誤差”則表示由本文算法得到的位置與準(zhǔn)確值之間的差異。

在物理模擬中，模擬速度和精度往往是一對(duì)矛盾。本文綜合模擬速度和精確度對(duì)Waterbomb結(jié)構(gòu)的彎曲彈性系數(shù)進(jìn)行選擇，并設(shè)置結(jié)構(gòu)彈簧拉伸系數(shù)均為80。

在單彎曲質(zhì)點(diǎn)彈簧系統(tǒng)中，本文采用不同的彎曲彈性系數(shù)檢測(cè)算法收斂性。彎曲彈性系數(shù)分別為10、20、100、500時(shí)的相對(duì)誤差值，如圖6所示。其中，相對(duì)誤差值是綜合第30、60、90幀后得到的平均值。

從圖6可以看出，本文將Tiantian Liu等[8]提出的質(zhì)點(diǎn)彈簧快速模擬方法運(yùn)用在Waterbomb樣模鋪嵌模型的模擬時(shí)，對(duì)于不同的彎曲彈性系數(shù)，該方法最終均可收斂。整體而言，彎曲彈性系數(shù)越低，收斂速度越快。

此外，實(shí)驗(yàn)結(jié)果證實(shí)，相對(duì)誤差并不是迭代次數(shù)的遞減函數(shù)。當(dāng)?shù)螖?shù)為20次或30次時(shí)，相對(duì)誤差增大，有些幀的相對(duì)誤差遠(yuǎn)遠(yuǎn)大于1。相對(duì)誤差描述了快速模擬方法相對(duì)于牛頓法的收斂情況。對(duì)實(shí)驗(yàn)結(jié)果整理后發(fā)現(xiàn)，[g（x0）-g（x*）<0]即[g（x0）

本文分別在圖6和圖7中描繪了彎曲彈性系數(shù)為10、20、100、500時(shí)相對(duì)誤差和絕對(duì)誤差的曲線，經(jīng)綜合分析可知，相對(duì)誤差和絕對(duì)誤差最終均會(huì)收斂。

除去彎曲彈性系數(shù)[bs=20]的情況，可以發(fā)現(xiàn)彈性彎曲系數(shù)越大，相對(duì)誤差和絕對(duì)誤差越小，收斂速度越快。當(dāng)[bs=10]時(shí)，相對(duì)誤差和絕對(duì)誤差尤其?。合鄬?duì)誤差最大在1左右。

綜上分析，當(dāng)彎曲彈性系數(shù)[bs=10]時(shí)，實(shí)驗(yàn)中的相對(duì)誤差和絕對(duì)誤差均最小。由此，本文選擇彈性彎曲系數(shù)為10對(duì)Waterbomb樣模鋪嵌進(jìn)行模擬。

5.1.3 模擬速度分析

在物理模擬中，除了精確度外，模擬速度是一個(gè)重要因素，而模擬精確度和速度往往不能兼顧。實(shí)驗(yàn)記錄了迭代次數(shù)分別為10、15、20、30、40時(shí)各數(shù)值方法的計(jì)算效率數(shù)值，如表1所示。

在表1中，當(dāng)?shù)螖?shù)分別為10、15、20、30、40時(shí)，本文算法（local-global）的fps數(shù)值依次遞減，分別為19.5、14、11、7.5、5.6，即隨著迭代次數(shù)增多，模擬速度降低。牛頓迭代法的fps數(shù)值均為11，這是因?yàn)楸疚乃玫腤aterbomb模型的復(fù)雜度較低，當(dāng)?shù)螖?shù)分別為10、15、20、30、40時(shí)，牛頓法已經(jīng)收斂。

5.1.4 參數(shù)范圍確定

在單彎曲質(zhì)點(diǎn)彈簧系統(tǒng)中，從模擬效果、誤差及穩(wěn)定性、模擬速度方面綜合考慮，通過對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行討論分析，選擇迭代次數(shù)為15、彎曲彈性系數(shù)小于等于10。此時(shí)，模擬穩(wěn)定，誤差在允許范圍內(nèi)，并且模擬速度快于牛頓法。

5.2 雙彎曲質(zhì)點(diǎn)彈簧模型

5.2.1 實(shí)驗(yàn)效果

根據(jù)實(shí)驗(yàn)結(jié)果，選擇相同彎曲彈性系數(shù)[sb=5]，不同彎曲彈性系數(shù)[db=10]，迭代次數(shù)為15，對(duì)Waterbomb樣模鋪嵌模型進(jìn)行模擬。實(shí)驗(yàn)結(jié)果如圖8所示。

5.2.2 誤差與穩(wěn)定性分析

雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)的Waterbomb折紙模型模擬的相對(duì)誤差和絕對(duì)誤差曲線如圖9、圖10所示。

從圖9可以看出，雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)中的相對(duì)誤差比單彎曲質(zhì)點(diǎn)彈簧系統(tǒng)更穩(wěn)定、誤差更小（保持在1附近及以下）。從絕對(duì)誤差方面來看，當(dāng)[sb=5，db=10]時(shí)，雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)的絕對(duì)誤差數(shù)量級(jí)為10-4，在誤差允許范圍內(nèi)，如圖10所示。相對(duì)誤差和絕對(duì)誤差的計(jì)算綜合了第30、60、90幀誤差的平均值。圖10幀數(shù)較少，大致描繪了誤差走勢(shì)。

在物理模擬中，模擬速度和精度往往是一對(duì)矛盾。根據(jù)對(duì)Waterbomb的結(jié)構(gòu)分析，綜合模擬速度、精確度及穩(wěn)定性，本文選擇彎曲彈性系數(shù)為5、不同彎曲彈性系數(shù)為10、迭代次數(shù)為15時(shí)對(duì)Waterbomb折紙模型進(jìn)行物理模擬。此時(shí)，模擬穩(wěn)定性良好，誤差在允許范圍內(nèi)，并且速度快于牛頓法。

5.2.3 參數(shù)確定

首先，在單彎曲質(zhì)點(diǎn)彈簧模型上，對(duì)模擬Waterbomb鋪嵌模型的參數(shù)進(jìn)行初步確定;然后，根據(jù)Waterbomb鋪嵌模型結(jié)構(gòu)特征，建立雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)。實(shí)驗(yàn)結(jié)果證明，相對(duì)于單彎曲質(zhì)點(diǎn)彈簧系統(tǒng)，雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)的穩(wěn)定性更好、誤差進(jìn)一步降低;最終選擇模擬參數(shù)為：相同彎曲彈性系數(shù)等于5、不同彎曲彈性系數(shù)等于10、迭代次數(shù)等于15。此時(shí)，算法模擬速度快于牛頓法。

6 結(jié)語(yǔ)

本文根據(jù)基本質(zhì)點(diǎn)彈簧模型確定參數(shù)范圍，為Waterbomb鋪嵌的每條邊均建立結(jié)構(gòu)彈簧，并以每條非邊界邊為基礎(chǔ)建立彎曲彈簧。采用隱式歐拉法和塊坐標(biāo)下降法（或交替優(yōu)化法）求解方程，得到物體運(yùn)動(dòng)變形的下一位置坐標(biāo)。通過對(duì)實(shí)驗(yàn)效果的誤差及穩(wěn)定性分析發(fā)現(xiàn)，當(dāng)?shù)螖?shù)為15、彎曲彈性系數(shù)小于等于10時(shí)，模擬效果穩(wěn)定、誤差在允許范圍內(nèi)且效果真實(shí)。然后，根據(jù)與棱相鄰的三角形面片是否對(duì)稱，將非邊界棱歸為兩大類，并設(shè)定不同的彎曲彈性系數(shù)，建立雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)。實(shí)驗(yàn)結(jié)果表明，相對(duì)于單彎曲質(zhì)點(diǎn)彈簧系統(tǒng)，雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)進(jìn)一步提升了穩(wěn)定性、減小了誤差、提高了模擬真實(shí)性。本文實(shí)現(xiàn)了對(duì)Waterbomb鋪嵌模型的穩(wěn)定模擬，在誤差允許范圍內(nèi)，模擬速度比牛頓法更快。雙彎曲質(zhì)點(diǎn)彈簧系統(tǒng)進(jìn)一步提高了模擬穩(wěn)定性、降低了模擬誤差。

本文沒有建立對(duì)角相望彈簧。在Waterbomb樣模中，如圖4中的E點(diǎn)和F點(diǎn)之間會(huì)發(fā)生彎曲，但是在本文算法中并沒有建立該彎曲彈簧。因?yàn)樵搹澢鷱椈墒腔邳c(diǎn)的，而本文彎曲彈簧是基于棱建立的。因此，在以后的研究中，可嘗試基于點(diǎn)建立質(zhì)點(diǎn)彈簧系統(tǒng)。另外，本文將與棱相鄰的三角形面片是軸對(duì)稱還是中心對(duì)稱視為同種情況（對(duì)稱）處理的做法不夠精確。由于彈簧對(duì)二面角的符號(hào)不敏感，所以本文并不是嚴(yán)格意義上的薄殼，但是由于Waterbomb鋪嵌的結(jié)構(gòu)特性，運(yùn)用質(zhì)點(diǎn)彈簧對(duì)其進(jìn)行模擬的做法也得到了很好的效果。

參考文獻(xiàn)：

[1] GRINSPUN E，HIRANI A，DESBRUN M，et al. Discrete shells [C]. Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation， 2003： 62-67.

[2] RAFSANJANI A， PASINI D. Bistable auxetic mechanical metamaterials inspired by ancient geometric motifs[J]. Extreme Mechanics Letters，2016， 9： 291-296.

[3] CHEN Y，F(xiàn)ENG H，MA J，et al. Symmetric waterbomb origami [J]. Proceedings of the Royal Society of London，2016，472 （2190）：1-20.

[4] TACHI T. Freeform variations of origami[J]. Journal for Geometry and Graphics，2010， 14（2）： 203-215.

[5] LU S， ZHAO C，GUO X，et al. Bi-layered mass-spring model for leaf motions[J]. Journal of System Simulation，2009，21（14）：4343-4388.

[6] 劉雪梅，毛磊，李運(yùn)華，等. 耦合無網(wǎng)格伽遼金與質(zhì)點(diǎn)彈簧實(shí)現(xiàn)軟組織形變仿真[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2013，25（1）：1-6.

[7] 陳瑛，吳明珠，胡新榮. 一種改進(jìn)的質(zhì)點(diǎn)-彈簧織物建模與仿真[J]. 華南師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2013， 45（5）：43-46.

[8] LIU T，BARGTEIL A W，O'BRIEN JAMES F，et al. Fast simulation of mass-spring systems[J]. ACM Transactions on Graphics，2013，32（6）：803-812.

[9] 王崴，周誠(chéng)，楊云，等. 基于改進(jìn)彈簧-質(zhì)點(diǎn)模型的柔性繩索仿真[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2015，27（11）：2230-2236.

[10] 王舒珍，陳國(guó)棟，李琦. 基于改進(jìn)的質(zhì)點(diǎn)-彈簧模型的軟組織形變仿真研究[J]. 通化師范學(xué)院學(xué)報(bào)，2017（1）：64-69.

[11] TERZOPOULOS D，PLATT J，BARR A，et al. Elastically deformable models[C]. Conference on Computer Graphics & Interactive Tec，1987：205-214.

[12] GARG A，GRINSPUN E，WARDETZKY M，et al. Cubic shells[C]. Siggraph/Eurographics Symposium on Computer Animation， 2007： 91-98.

[13] BOUAZIZ S，MARTIN S，LIU T，et al. Projective dynamics： fusing constraint projections for fast simulation[J]. ACM Transactions on Graphics，2014，33（4）：1-11.

[14] ZHANG X，LE X，WU Z，et al. Data-driven bending elasticity design by shell thickness[J]. Computer Graphics Forum，2016，35（5）：157-166.

[15] LIU T，BOUAZIZ S，KAVAN L，et al. Quasi-Newton methods for real-time simulation of hyperelastic meterials[J]. ACM Transactions on Graphics，2017，36（3）：1-16.

[16] BURGOON R，WOOD Z J，GRINSPUN E. Discrete shells origami[C]. Proceedings of the 21st International Conference on Computers and Their Applications， 2006：180-187.

[17] TACHI T. Origanizing polyhedral surfaces[J]. IEEE Transactions on Visualization and Computer Graphics，2010， 16（2）： 298-311.

[18] TACHI T. Designing freeform origami tessellations by generalizing Resch's patterns[C]. Proceedings of ASME 2013 IDETC/CIE， 2013：1-10.

[19] TACHI T，OBREBSKI J B， TARCZEWSKI R. Interactive form-finding of elastic origami[C]. Proceedings of the International Association for Shell and Spatial Structures（IASS） Symposium，2013：23-27.

[20] DAVIS E， DEMAINE E D， DEMAINE M. L， et al. Reconstructing David Huffman's origami tessellations [C]. Proceedings of the ASME 2013 International Design Engineering Technical Conferences & Computersand Information in Engineering Conference， 2013：1-10.

（責(zé)任編輯：江 艷）