一種SPH流體仿真邊界校正方法

劉旭，班曉娟，楊鳴遠(yuǎn)，賀亮亮

（北京科技大學(xué)計(jì)算機(jī)與通信工程學(xué)院，北京100083）

一種SPH流體仿真邊界校正方法

劉旭，班曉娟，楊鳴遠(yuǎn)，賀亮亮

（北京科技大學(xué)計(jì)算機(jī)與通信工程學(xué)院，北京100083）

使用光滑粒子流體動(dòng)力學(xué)方法進(jìn)行流體仿真，并提出一種邊界校正方法。使用快速泊松盤(pán)采樣算法對(duì)容器邊界進(jìn)行采樣，生成邊界粒子，對(duì)邊界粒子質(zhì)量進(jìn)行差值估算，計(jì)算邊界粒子對(duì)流體粒子的作用力，以此來(lái)仿真流體與邊界的相互作用。該方法可避免穿刺、滯留等現(xiàn)象的發(fā)生。通過(guò)實(shí)驗(yàn)驗(yàn)證了該算法的正確性。

流體仿真；光滑粒子流體動(dòng)力學(xué)；邊界校正

計(jì)算機(jī)圖形學(xué)的一個(gè)重要目標(biāo)就是逼真地再現(xiàn)自然界各物體的運(yùn)動(dòng)及形態(tài)。流體作為自然界的重要組成部分，其運(yùn)動(dòng)復(fù)雜而多變，流體仿真一直以來(lái)都是圖形學(xué)工作者研究的難點(diǎn)與熱點(diǎn)。光滑粒子流體動(dòng)力學(xué)((smoothed particle hydrodynamics, SPH)[1-2]方法最早用于解決三維開(kāi)放空間的天體物理學(xué)問(wèn)題，由于其去網(wǎng)格化對(duì)于仿真流體的潑濺、卷浪等現(xiàn)象具有天然優(yōu)越性，因此近年來(lái)越來(lái)越多的研究者開(kāi)始研究流體仿真SPH方法。Müller等[3]首次開(kāi)創(chuàng)性地提出了SPH流體仿真的算法框架。Becker和Teschner[4]對(duì)不可壓縮流體仿真進(jìn)行研究，提出了弱可壓縮SPH方法((weakly compressible SPH,WCSPH)。他使用Tait方程取代氣體狀態(tài)方程，利用對(duì)稱的壓力公式并提出了一種新的表面張力建模方法。Solenthaler和Pajarola[5]提出了一種預(yù)測(cè)矯正 SPH法(predictive-corrective incompressible SPH,PCISPH)，該方法可大大增加仿真的時(shí)間步長(zhǎng)，從而提高了仿真系統(tǒng)的整體效率。邊界條件是SPH流體仿真的一個(gè)重要研究點(diǎn)，許多研究人員將剛體邊界使用粒子來(lái)表示，Monaghan等[6]通過(guò)給與剛體邊界粒子相鄰的流體粒子施加邊界力來(lái)解決。Ihmsen等[7]將流體粒子密度外推至邊界粒子，然后對(duì)流體重新計(jì)算。Schechter和Bridson[8]對(duì)其改進(jìn)，使用一種ghostSPH方法，得到了令人滿意的結(jié)果。

SPH方法流體仿真計(jì)算量及大，許多研究人員對(duì)于SPH加速問(wèn)題進(jìn)行了研究。文獻(xiàn)[9]中使用自適應(yīng)大小粒子方法提高效率。Solenthaler和Gross[10]根據(jù)細(xì)節(jié)需要使用兩種大小的粒子對(duì)流體進(jìn)行模擬，減少了粒子數(shù)量從而提高了效率。文獻(xiàn)[11]提出了一種基于背景網(wǎng)格的并行化算法。文獻(xiàn)[12]針對(duì)大規(guī)模場(chǎng)景不宜使用背景網(wǎng)格的缺點(diǎn)，使用hash表來(lái)存儲(chǔ)粒子。文獻(xiàn)[13]考慮到Cache命中率問(wèn)題，針對(duì)背景網(wǎng)格使用粒子Z排序方法，針對(duì)大規(guī)模場(chǎng)景的hash存儲(chǔ)提出了緊湊hash策略，該方法可提升2～3倍的運(yùn)行效率。很多研究者也對(duì)使用GPU并行計(jì)算SPH進(jìn)行了研究[14-15]。更多的流體仿真知識(shí)可參考文獻(xiàn)[1]，更多SPH知識(shí)可參考文獻(xiàn)[2]。

流體的運(yùn)動(dòng)往往伴隨與其他物體的交互，因此使用SPH方法對(duì)流體和其他物體的交互進(jìn)行模擬也是一個(gè)研究重點(diǎn)。Mao和Yang[16]對(duì)不可混溶流體交互進(jìn)行研究，提出了一種不可混溶流體交互仿真方法。Du等[17]使用三角形網(wǎng)格表示可變形體，并通過(guò)持續(xù)碰撞檢測(cè)實(shí)現(xiàn)了一種流體與可變形體交互方法。Akinci等[18]提出了一種通用的流體剛體仿真方法，該方法使用一種矯正的密度評(píng)估方法來(lái)計(jì)算流體和剛體的相互作用力。SPH方法使用粒子對(duì)流體建模，在動(dòng)畫(huà)制作中需要對(duì)流體自由表面進(jìn)行表面重構(gòu)以進(jìn)行后期渲染。Mathieu和Cani-Gascuel[19]使用Blobby球方法生成的表面，并施加“表面張力”來(lái)平滑表面，從而克服了不自然的凹凸效果；Müller等[3]用SPH方法代替三維高斯函數(shù)對(duì)空間標(biāo)量場(chǎng)進(jìn)行疊加；Zhu和Bridson[20]在Blobby球方法的基礎(chǔ)上對(duì)粒子附近的密度變化進(jìn)行補(bǔ)償，并進(jìn)一步對(duì)空間標(biāo)量場(chǎng)進(jìn)行了一次平滑，從而得到更光滑的表面。Yu和Turk[21]提出了一種用各向異性核函數(shù)(anisotropickernel)構(gòu)造表面的方法。陳沸鑌等[22-23]則進(jìn)一步用SPH方法模擬了剛體的破裂等現(xiàn)象。

但是由于剛體建模與流體建模不同，剛體表面由網(wǎng)格構(gòu)成且粒子分布較為稀疏。而基于SPH的流體模擬，流體由粒子構(gòu)成。這樣剛體與流體的接觸面上兩種粒子的數(shù)量差極大，不但會(huì)影響雙方受力的計(jì)算，還會(huì)產(chǎn)生液體粒子穿進(jìn)或穿過(guò)剛體的現(xiàn)象。因此本文將針對(duì)復(fù)雜場(chǎng)景流體運(yùn)動(dòng)進(jìn)行SPH仿真，使用泊松盤(pán)對(duì)場(chǎng)景邊界進(jìn)行采樣，此處將剛體視為場(chǎng)景的一部分，構(gòu)建場(chǎng)景邊界粒子，并估算其質(zhì)量，使用校正的公式來(lái)計(jì)算流體粒子與邊界的交互作用。

1 SPH流體仿真

1.1 SPH方法

SPH方法的核心思想是以離散化粒子的形式來(lái)表征連續(xù)的場(chǎng)，并對(duì)場(chǎng)物理量使用積分近似的方式進(jìn)行計(jì)算。場(chǎng)量A（x）可用下式近似計(jì)算：

其中，mj、ρj分別表示粒子質(zhì)量和密度；W （x-xj,h）為積分光滑核函數(shù)，一般滿足3個(gè)條件：歸一化條件、緊支持性和Dirac函數(shù)性；h為其積分半徑。核函數(shù)可取滿足條件的多個(gè)函數(shù)，文獻(xiàn)[2]中使用B-樣條函數(shù)，文獻(xiàn)[3]使用Spiky核函數(shù)并且在計(jì)算粘性力時(shí)設(shè)計(jì)了專(zhuān)門(mén)的核函數(shù)，此外，高斯函數(shù)也經(jīng)常被用做光滑核函數(shù)。本文使用三次B-樣條函數(shù)作為核函數(shù)，取h=4r，r為粒子半徑。

密度：流體運(yùn)動(dòng)滿足質(zhì)量守恒，即滿足連續(xù)方程：

其中，v表示粒子速度。在SPH流體仿真中，粒子密度使用插值進(jìn)行計(jì)算，由式(1)可知粒子密度計(jì)算公式：

該公式可保證系統(tǒng)質(zhì)量守恒。由該式可知粒子密度僅與當(dāng)前狀態(tài)周?chē)Ｗ臃植加嘘P(guān)，密度的計(jì)算是SPH流體仿真的基礎(chǔ)。只有當(dāng)周?chē)Ｗ映尸F(xiàn)各向球狀均勻分布時(shí)，使用上式可得近似結(jié)果，而在流體邊界處粒子分布不滿足該條件，將導(dǎo)致粒子密度值偏小，導(dǎo)致邊界問(wèn)題，這將在后一章進(jìn)一步討論。

壓力：流體運(yùn)動(dòng)滿足動(dòng)量守恒，流體運(yùn)動(dòng)動(dòng)量方程如下：

其中，P表示壓強(qiáng)，g表示外力場(chǎng)。通過(guò)推導(dǎo)[2]，可使用下式計(jì)算粒子所受壓力：

由該式可知，相鄰粒子相互作用力是對(duì)稱的，且能保證線動(dòng)量和角動(dòng)量守恒。

狀態(tài)方程：通過(guò)式(3)求得粒子密度之后，需要對(duì)粒子壓強(qiáng)進(jìn)行計(jì)算，以求得粒子壓力。文獻(xiàn)[3]使用氣體狀態(tài)方程P=k(ρ-ρ0)來(lái)計(jì)算粒子壓強(qiáng)，其中，k為常數(shù)， ρ0為常量密度。該方法將導(dǎo)致比較大的可壓縮性。本文使用文獻(xiàn)[4]中的Tait方程進(jìn)行密度計(jì)算：

1.2 人工粘度

為了保證數(shù)值計(jì)算穩(wěn)定性，在SPH流體仿真中通常加入粘性力，文獻(xiàn)[3]中對(duì)流體運(yùn)動(dòng)控制NS方程中的粘力項(xiàng)進(jìn)行推導(dǎo)得到粘力計(jì)算公式：

該式計(jì)算的粘力并不對(duì)稱。文獻(xiàn)[2]中提出了一種人工粘性力：

其中，

式中，α是可調(diào)參數(shù)，ε=0.01，εh2的引入是為避免xij2=0。

本文使用式(8)計(jì)算粘性力，該方法產(chǎn)生的人工粘力是對(duì)稱的，即可以保證數(shù)字計(jì)算的穩(wěn)定性，同時(shí)還可以模擬流體運(yùn)動(dòng)的粘性效果。

1.3 流體表面張力

表面張力是流體的一大重要特性，文獻(xiàn)[3-4]都對(duì)表面張力進(jìn)行了研究，本文采用文獻(xiàn)[4]的方法。該方法基于這樣一個(gè)事實(shí)，表面張力的微觀機(jī)理是分子間引力，通過(guò)在相鄰粒子間加入引力來(lái)對(duì)表面張力進(jìn)行模擬，如下式：

其中，κ為可調(diào)參數(shù)，用來(lái)控制表面張力的強(qiáng)弱。在流體內(nèi)部，由于周?chē)Ｗ臃植季鶆?，各周?chē)Ｗ訉?duì)粒子i產(chǎn)生的表面張力合力項(xiàng)為0，在流體表層則不為0。

1.4 時(shí)間積分

仿真的每一步，先根據(jù)粒子分布，計(jì)算出各粒子密度ρi，然后依據(jù)式(6)求得粒子壓強(qiáng) Pi，最后分別計(jì)算粒子壓力FiP、粘性力Fiv和表面張力Fie，則每個(gè)粒子加速度為：

其中，g為外力場(chǎng)，通常為重力。之后需要對(duì)粒子進(jìn)行時(shí)間積分，常用的時(shí)間積分方法有：蛙跳法(leap-frog,LF)、預(yù)測(cè)校正法(predictor-corrector)、龍格-庫(kù)塔法(Runge-Kutta,RK)等，本文采用簡(jiǎn)單易用的蛙跳法進(jìn)行實(shí)現(xiàn)，可參考文獻(xiàn)[25]。積分時(shí)間步長(zhǎng)通過(guò)CFL(courant-friedrichs-lewy)條件[4]確定。

2 基于質(zhì)量校正的邊界處理

流體的邊界包括流體和空氣的交界面以及流體和容器壁的交界面，本文只討論流體和容器壁交界面的處理。為了簡(jiǎn)化處理過(guò)程，將容器看成剛體。本文結(jié)合文獻(xiàn)[8]和文獻(xiàn)[11]的方法，通過(guò)引入邊界粒子來(lái)實(shí)現(xiàn)容器壁和流體的相互作用。首先將邊界采樣為邊界粒子，進(jìn)而通過(guò)計(jì)算流體粒子和邊界粒子之間的相互作用來(lái)計(jì)算容器壁和流體的相互作用。

2.1 邊界泊松采樣

為得到邊界粒子，需要對(duì)以網(wǎng)格形式存在的容器模型進(jìn)行采樣，本文對(duì)容器表面按泊松盤(pán)分布（Poisson disk distribution，一種藍(lán)色噪聲分布）進(jìn)行采樣，具體參考文獻(xiàn)[26]。注意，在整個(gè)仿真過(guò)程中對(duì)每個(gè)容器或剛體只進(jìn)行一次初始采樣，該采樣結(jié)果將在后續(xù)整個(gè)仿真過(guò)程中使用，因此，對(duì)容器或剛體的采樣并不會(huì)帶來(lái)額外的開(kāi)銷(xiāo)。表面采樣結(jié)果如圖1所示。

圖1剛體表面采樣結(jié)果

2.2 邊界粒子校正計(jì)算

考慮到邊界粒子的影響，流體粒子的密度計(jì)算式(2)需要加入邊界粒子的加權(quán)求和結(jié)果：

式中，用fi來(lái)表示流體粒子i，bk表示邊界粒子k，用j迭代粒子if的所有流體鄰居粒子，用k迭代所有邊界鄰居粒子。

文獻(xiàn)[11]對(duì)上式進(jìn)行了改進(jìn)?？紤]到邊界粒子質(zhì)量設(shè)置不合理或分布不均勻，將導(dǎo)致計(jì)算錯(cuò)誤，使用下式估算邊界粒子質(zhì)量：

其中，ρ0表示流體的常量密度，為邊界粒子所代表邊界區(qū)域體積的估算值。使用Ψbi（ρ0）取代邊界粒子質(zhì)量進(jìn)行計(jì)算可提高穩(wěn)定性。

因此，式(12)可改為：

容器壁和流體最重要的相互作用就是壓力，用下式計(jì)算邊界粒子對(duì)流體粒子的壓力產(chǎn)生的加速度：

式中，pfi＞0時(shí)取k=2。pfi＜0時(shí)，邊界粒子和流體粒子相互吸引，可對(duì)k進(jìn)行調(diào)整(0≤k≤2)以實(shí)現(xiàn)不同吸附效果，本文中取k=1。

為模擬容器壁與流體之間的摩擦力或?qū)崿F(xiàn)流體和剛體的交互，需要計(jì)算邊界粒子和流體粒子摩擦力，摩擦力的計(jì)算借鑒了人工粘度公式(8)：

式中的Πik同式(9)。

2.3 算法流程

本文流體仿真算法流程如下：

3 實(shí)驗(yàn)結(jié)果

本文通過(guò)實(shí)驗(yàn)驗(yàn)證了算法的有效性，程序的運(yùn)行平臺(tái)為Inteli5-3470（四核，3.20GHz，6MB Cache）、8GB內(nèi)存。仿真算法和表面重構(gòu)算法用C++語(yǔ)言采用多線程技術(shù)實(shí)現(xiàn)，仿真算法中的臨近粒子搜索使用空間背景網(wǎng)格進(jìn)行哈希查找。表面重構(gòu)采用文獻(xiàn)[15]的方法，使用各向異性核函數(shù)構(gòu)造顏色場(chǎng)，然后使用步進(jìn)立方體算法重構(gòu)表面，其中矩陣奇異值分解采用的是NIST的JAMA和TNT開(kāi)源數(shù)學(xué)庫(kù)。使用OpenGL三維圖形庫(kù)實(shí)時(shí)顯示仿真和表面構(gòu)造的結(jié)果，并使用OpenCV庫(kù)錄制視頻。為了對(duì)復(fù)雜容器和場(chǎng)景進(jìn)行建模以及制作動(dòng)畫(huà)，本文使用Blender軟件，后期高質(zhì)量流體效果的渲染使用光線跟蹤引擎Mitsuba。

為了驗(yàn)證邊界處理算法的正確性，本文進(jìn)行了噴泉仿真實(shí)驗(yàn)與城市洪澇實(shí)驗(yàn)。噴泉實(shí)驗(yàn)中流體粒子數(shù)約為100～200k，邊界粒子數(shù)150k，粘度系數(shù)α=0.08，表面張力系數(shù)κ=0.01，時(shí)間步長(zhǎng)為0.0002s。圖2為噴泉實(shí)驗(yàn)結(jié)果。

圖2 噴泉實(shí)驗(yàn)結(jié)果

城市洪澇實(shí)驗(yàn)中流體粒子數(shù)約為50～200k，邊界粒子數(shù)80k，粘度系數(shù)α=0.08，表面張力系數(shù)κ=0.005，時(shí)間步長(zhǎng)為0.0005s。圖3為城市洪澇實(shí)驗(yàn)結(jié)果。

手?jǐn)噭?dòng)水的實(shí)驗(yàn)中流體粒子數(shù)約為1.17M，邊界粒子數(shù)63k，粘度系數(shù)α=0.05，表面張力系數(shù)κ=0.005，時(shí)間步長(zhǎng)為0.000134s，如圖4所示。

圖3 城市洪澇實(shí)驗(yàn)不同時(shí)刻渲染結(jié)果

圖4 手?jǐn)噭?dòng)水的模擬

4結(jié)論

本文通過(guò)實(shí)驗(yàn)驗(yàn)證了所使用的SPH流體仿真和邊界處理算法的有效性，并制作了噴泉?jiǎng)赢?huà)，實(shí)驗(yàn)結(jié)果逼真地反應(yīng)了潑濺、吸附、表面張力、粘性等流體特征。但是算法計(jì)算量較大，今后工作將對(duì)如何提高算法效率進(jìn)行研究，此外還將研究流體剛體交互計(jì)算方法。

[1] Bridson R, Müller-fischer M. Fluid simulation: SIGGRAPH 2007 course notes video filesassociated w ith this course are available from the citation page[C]//ACM SIGGRAPH 2007Courses.New York,USA,2007:1-81.

[2]Monaghan J J.Smoothed particle hydrodynam ics[J]. Reportson Progress in Physics,2005,68(8):1703-1759.

[3]Müller M,Charypar D,Gross M.Particle-based fluid simulation for interactive applications [C]//ACM SIGGRAPH/Eurographics Symposium on Computer Animation.San Diego,CA,USA,2003:154-159.

[4]Becker M,Teschner M.Weakly compressible SPH for free surface flows[C]//ACM SIGGRAPH/Eurographics symposium on Computer animation.Sw itzerland,2007: 209-217.

[5] Solenthaler B, Pajarola R. Predictive-corrective incompressible SPH[J].ACM Transactions on Graphics, 2009,28(3):40:1-40:6.

[6]Monaghan JJ,Kos A,Issa N.Fluidmotion generated by impact[J].Journal of Waterway Port Coastal and Ocean Engineering,2003,129(6):250-259.

[7]Ihmsen M,Akinci N,Gissler M.Boundary handling and adaptive time-stepping for PCISPH[C]//Workshop in Virtual Reality Interactions and Physical Simulation, VRIPHYS(2010).Darmstadt,Germany,2010:79-88.

[8] Schechter H,Bridson R.Ghost SPH for animating water[J]. ACM Transcationson Graphics,2012,31(4):61:1-61:8.

[9]Adams B,Pauly M,Keiser R,et al.Adaptively sampled particle fluids[J].ACM Transactions on Graphics,2007, 26:48:1-48:8

[10]Solenthaler B,GrossM.Two-scale particle simulation[J]. ACM Transactionson Graphics,2011,30(4):81:1-81:8.

[11]Kalojanov J,Slusallek P.A parallel algorithm for construction of uniform grids [C]//In HPG '09: Proceedings of the 1st ACM Conference on High PerformanceGraphics.New York,USA,2009:23-28.

[12]Bell N,Yu Yizhou,Mucha P J.Particle-based simulation of granularmaterials[C]//ACM SIGGRAPH/Eurographics Symposium on Computer Animation,New York,USA, 2005:77-86.

[13]Ihmsen M,Akinci N,Becker M,et al.A parallel SPH implementation on multicore CPUs [J].Computer Graphics Forum,2011,30(1):99-112.

[14]Harada T,Koshizuka S,Kawaguchi Y.Smoothed particle hydrodynamics on GPUs [C]//In Proceedings of Computer Graphics International.Petropolis,RJ,Brazil, 2007:63-70.

[15]Goswam i P,Schlegel P,Solenthaler B,et al.Interactive SPH simulation and rendering on the GPU [C]// Proceedings of the 2010 ACM SIGGRAPH/Eurographics Symposium on Computer Animation.Madrid,Spain, 2010:55-64.

[16]Mao Hai,Yang Y.Fluid-fluid collision for particle-based immiscible fluid-fluid collision[C]//GI'06 Proceedingsof Graphics Interface 2006.Toronto Canada,2006:49-55.

[17]Du Peng,Tang M in,Meng Chang,et al.A fluid/cloth coupling method for high velocity collision simulation[C]// Proceedings of the 11th ACM SIGGRAPH International Conference on Virtual-Reality Continuum and its Applications in Industry.Singapore,2012:309-314.

[18]AkinciN,Ihmsen M,AkinciG,etal.Versatile rigid-fluid coupling for incompressible SPH[J].ACM Transactions on Graphics,2012,31:62:1-62:8.

[19]Mathieu D,Cani-Gascuel M.Active implicit surface for animation [C]//Graphics Interface.Vancouver,BC, Canada,1998:143-150.

[20]Zhu Yongning,Bridson R.Animating sand as a fluid[J]. ACM Transactionson Graphics,2005,24(3):965-972.

[21]Yu Jihun, Turk G. Reconstructing surfaces of particle-based fluids using anisotropic kernels[J].ACM Transactionson Graphics,2013,32(5):1-5,12.

[22]陳沸鑌,王長(zhǎng)波,謝步瀛.可控性的雙層粒子剛體脆性破裂模擬動(dòng)畫(huà)[J].圖學(xué)學(xué)報(bào),2015,36(1):111-116.

[23]陳沸鑌,謝步瀛,冉修遠(yuǎn).基于細(xì)分粒子的剛體破裂動(dòng)畫(huà)[J].圖學(xué)學(xué)報(bào),2014,35(5):669-675.

[24]Liu G R,Liu M B.光滑粒子流體動(dòng)力學(xué)——一種無(wú)網(wǎng)格粒子法[M].韓 旭,楊 剛,強(qiáng)洪夫,譯.長(zhǎng)沙:湖南大學(xué)出版社,2005:1-433.

[25]Desbrun M,CaniM.Smoothed particles:a new paradigm for animating highly deformable bodies [C]//6th Eurographics Workshop on Computer Animation and Simulation'96.Poitiers,France,1996:61-76.

[26]Bridson R.Fast Poisson disk sampling in arbitrary dimensions [C]//In ACM SIGGRAPH Technical Sketches2007.San Diego,CA,USA,2007:1.

A Boundary Correction Method of SPH Fluid Simulation

Liu Xu, Ban Xiaojuan, Yang Mingyuan, He Liangliang

(Department of Computer Science and Technology, University of Science and Technology Beijing, Beijing 100083, China)

This article presents a boundary correction method in fluid simulation based on smoothed particle hydrodynamics. We sample boundaries as boundary particles with fast Poisson disk algorithm,interpolate mass of boundary particles, and calculate force between fluid particles and boundaryparticles as simulation of fluid-boundary interaction. The presented method can also avoid penetration and holdup. Experiments demonstrated conduct experiments to prove the validity of our method.

fluid simulation; smoothed particle hydrodynamics; boundary correction

TP391

A

2095-302X(2015)03-0462-06

2014-12-26；定稿日期：2015-02-08

國(guó)家自然科學(xué)基金資助項(xiàng)目(61272357，61300074)

劉旭(1987-)，男，北京人，博士研究生。主要研究方向?yàn)橛?jì)算機(jī)圖形學(xué)、流體模擬及表面重構(gòu)。E-mail：liuxu.ustb@gmail.com

班曉娟(1970-)，女，遼寧朝陽(yáng)人，教授，博士。主要研究方向?yàn)槿斯ぶ悄堋⒂?jì)算機(jī)動(dòng)畫(huà)。E-mail：Banxj@ustb.edu.cn