抗磁懸浮振動(dòng)能量采集器動(dòng)力學(xué)響應(yīng)的仿真分析?

秦立振 張振宇 張坤 丁建橋 段智勇 蘇宇鋒

1)(鄭州大學(xué)機(jī)械工程學(xué)院,鄭州 450001)

2)(鄭州大學(xué)物理工程學(xué)院,鄭州 450001)

1 引 言

近年來(lái),植入式醫(yī)療設(shè)備、無(wú)線傳感器網(wǎng)絡(luò)等得到了廣泛的應(yīng)用,但是這些設(shè)備的電能供應(yīng)仍然依賴(lài)于傳統(tǒng)電池[1].傳統(tǒng)電池具有使用方便、成本低廉的優(yōu)點(diǎn),但同時(shí)也具有持久性差、儲(chǔ)能有限的缺點(diǎn).在一些無(wú)法更換電池或更換電池的成本比較高的場(chǎng)所,傳統(tǒng)電池顯現(xiàn)了其固有的局限性.因此,開(kāi)發(fā)一種在用電設(shè)備有效壽命結(jié)束之前無(wú)需維護(hù)的永久電池成為各國(guó)研究者關(guān)注的課題之一.自然環(huán)境中存在著各種各樣的能量,利用環(huán)境中的能量為電子元器件供電是取代傳統(tǒng)電池的有效途徑[2].

能量采集技術(shù)的主要原理是利用光伏效應(yīng)、壓電效應(yīng)、熱電效應(yīng)、電磁感應(yīng)等,將周?chē)h(huán)境中的能量轉(zhuǎn)化為電能,然后利用儲(chǔ)能元件進(jìn)行電能存儲(chǔ),最后通過(guò)相應(yīng)的電源管理電路,將存儲(chǔ)的電能提供給外部負(fù)載,進(jìn)而為用電設(shè)備供電.能量采集不僅涉及能量轉(zhuǎn)化利用的問(wèn)題,同時(shí)也涉及優(yōu)質(zhì)能量來(lái)源的問(wèn)題.相對(duì)于光、風(fēng)、熱等,振動(dòng)能量具有分布范圍廣、能量密度高、易于采集的特點(diǎn),因此振動(dòng)能量采集是目前研究較多的一種能量采集方案[3].

振動(dòng)能量的采集存在不同的途徑,基于不同的原理、方法和材料.1939年,Braunbek和Schweben[4]首次利用非均勻電磁場(chǎng)實(shí)現(xiàn)了微小片狀石墨的穩(wěn)定懸浮,此后,抗磁懸浮得到了越來(lái)越多的關(guān)注和研究.意大利都靈理工大學(xué)的de Pasquale等[5]利用釹鐵硼(NdFeB)強(qiáng)磁鐵和高定向熱解石墨(HOPG)抗磁材料,在常溫下獲得了1—2 mm的懸浮高度.抗磁懸浮由于能夠在常溫下實(shí)現(xiàn)無(wú)摩擦穩(wěn)定懸浮,所以得到了越來(lái)越多的應(yīng)用.目前,抗磁懸浮原理已被應(yīng)用到生物醫(yī)學(xué)[6,7]、傳感器[8]、微執(zhí)行器[9]、微型電機(jī)[10]等領(lǐng)域.基于這一原理,我們課題組開(kāi)展了微型抗磁懸浮振動(dòng)能量采集器的研究[11,12],本文研究了該能量采集器的振動(dòng)響應(yīng)特性.

2 結(jié)構(gòu)模型和工作原理

抗磁懸浮振動(dòng)能量采集器的主要部件包括:提升永磁體、懸浮永磁體、上熱解石墨板、下熱解石墨板以及布置在上下熱解石墨板內(nèi)表面的平面螺旋線圈,其結(jié)構(gòu)如圖1所示.

圖1 抗磁懸浮振動(dòng)能量采集器的結(jié)構(gòu)模型Fig.1.The structure model of diamagnetic levitation vibration energy harvester.

由于提升永磁體和懸浮永磁體都是圓柱形且磁化方向都沿軸向,所以提升永磁體的磁場(chǎng)會(huì)在水平方向上形成一個(gè)磁勢(shì)阱.根據(jù)能量最小原理,懸浮永磁體在水平方向上能夠穩(wěn)定懸浮且平衡點(diǎn)位于兩永磁體的公共軸線上.高定向熱解石墨是常溫下抗磁性最強(qiáng)的材料,表現(xiàn)為對(duì)外界磁場(chǎng)的排斥作用.在豎直方向上,懸浮永磁體同時(shí)受到重力、提升磁體吸引力的豎直分量、熱解石墨板的抗磁力.在滿足特定的尺寸參數(shù)條件下,豎直方向上也能夠穩(wěn)定懸浮,所以懸浮永磁體處于一個(gè)三維勢(shì)阱之中能夠穩(wěn)定懸浮.當(dāng)系統(tǒng)受到外界干擾時(shí),懸浮磁體可以在平衡位置附近振動(dòng),從而將外界振動(dòng)轉(zhuǎn)化為懸浮磁體的振動(dòng).根據(jù)法拉第電磁感應(yīng)定律,懸浮磁體振動(dòng)的機(jī)械能轉(zhuǎn)化為線圈內(nèi)部的電能[11].

3 結(jié)構(gòu)參數(shù)和受力分析

抗磁懸浮振動(dòng)能量采集器的結(jié)構(gòu)參數(shù)和材料參數(shù)直接決定了懸浮磁體的受力狀態(tài),進(jìn)而影響能量采集器的振動(dòng)響應(yīng)特性和能量轉(zhuǎn)化效率,其結(jié)構(gòu)材料參數(shù)如表1所列.

當(dāng)懸浮永磁體發(fā)生振動(dòng)時(shí),其不但會(huì)受到重力、抗磁力、提升永磁體的吸引力,還會(huì)受到空氣和感應(yīng)線圈的阻尼力.文獻(xiàn)[13—16]給出了永磁體之間以及永磁體和熱解石墨板之間受力的理論關(guān)系式,但這些關(guān)系式并不是顯式且過(guò)于復(fù)雜,不便于進(jìn)行動(dòng)力學(xué)響應(yīng)分析,因此采用COMSOLTM仿真并結(jié)合數(shù)值擬合的方法獲得磁力和抗磁力的近似表達(dá)式.本研究中懸浮永磁體受力的仿真計(jì)算模型如圖2(a)和圖2(b)所示.

表1 抗磁懸浮能量采集器結(jié)構(gòu)參數(shù)Table 1.The structure parameters of diamagnetic levitation energy harvester.

圖2 磁力和抗磁力計(jì)算模型Fig.2.The calculational model of magnetic force and diamagnetic force.

永磁體和熱解石墨板周?chē)目諝庥蚍譃閮蓪?外層設(shè)為無(wú)限域模擬磁體的彌散場(chǎng),且網(wǎng)格劃分較粗,內(nèi)層是有限域且網(wǎng)格劃分較細(xì),以節(jié)省仿真時(shí)間并兼顧仿真精度,磁體的和熱解石墨板的網(wǎng)格最細(xì)以提高仿真精度,單元類(lèi)型采用自由剖分四面體網(wǎng)格,最大生長(zhǎng)率為1.6,曲率因子為0.6,狹窄區(qū)域解析度為0.5.當(dāng)網(wǎng)格的精度不斷增加時(shí),仿真結(jié)果越趨近于定值且變化越來(lái)越慢.為了兼顧時(shí)間和精度,分析中網(wǎng)格最大單元尺寸選定如下:無(wú)限空氣域5 mm,有限空氣域2 mm,磁體和熱解石墨板0.8 mm.圖2(a)和圖2(b)中模型單元數(shù)分別為661370和920624.

對(duì)仿真結(jié)果進(jìn)行數(shù)據(jù)擬合求得磁力F1的解析式為

式中,H為提升永磁體和懸浮永磁體之間的距離,單位為mm;F1單位為μN(yùn).

抗磁力F2的解析式為

式中,h為懸浮永磁體和熱解石墨板之間的距離,單位為mm;F2單位為μN(yùn).

感應(yīng)線圈工作時(shí)懸浮永磁體受到的電磁阻力Fεγ

式中,Bri為第i匝線圈所在位置磁感應(yīng)強(qiáng)度在徑向上的分量,Bz為磁感應(yīng)強(qiáng)度在軸向上的分量,ri為第i匝線圈的半徑.

對(duì)(3)式利用響應(yīng)面法求得懸浮磁體周?chē)艌?chǎng)在軸向和徑向上的分布,最終求得其電磁阻尼力的表達(dá)式為

懸浮磁體的運(yùn)動(dòng)屬于低幅低頻運(yùn)動(dòng),因此其受到的空氣阻力可以用(5)式近似計(jì)算:

式中,S為迎風(fēng)面積,k為阻力系數(shù).

以?xún)蔁峤馐鍖?duì)稱(chēng)面的中心為零點(diǎn),向上設(shè)為正方向,統(tǒng)一坐標(biāo)系后,在豎直方向上振動(dòng)的動(dòng)力學(xué)方程如(6)式所示:

其中,a1,a3由F1和F2決定,ω=2πf為外界激勵(lì)頻率,A為激勵(lì)振幅.

4 運(yùn)動(dòng)特性

懸浮磁體所受磁力、抗磁力和重力的合力隨熱解石墨板間距L的變化曲線如圖3所示.通過(guò)分析其受力曲線可以發(fā)現(xiàn),其受力曲線的非線性程度隨著兩熱解石墨板間距的增大而增強(qiáng).

圖3 懸浮磁體所受合力Fig.3.The resultant curves of the fl oating magnet.

圖4 懸浮永磁體勢(shì)能曲線Fig.4.The potential energy curve of the fl oating magnet.

圖4是懸浮永磁勢(shì)能隨熱解石墨板間距變化的曲線.發(fā)現(xiàn)當(dāng)兩熱解石墨板間距L＞7.7 mm時(shí),勢(shì)能曲線會(huì)出現(xiàn)兩個(gè)勢(shì)阱,而在L＜7.7 mm時(shí),只存在一個(gè)勢(shì)阱.由能量最小原理可知,懸浮永磁體的穩(wěn)定平衡位置就是磁勢(shì)阱的位置.

圖5是當(dāng)熱解石墨板間距L分別為7 mm和8 mm時(shí),懸浮磁體的靜止平衡現(xiàn)象.由實(shí)驗(yàn)現(xiàn)象可知熱解石墨板間距的變化會(huì)導(dǎo)致懸浮永磁體穩(wěn)態(tài)類(lèi)型發(fā)生改變,即單穩(wěn)態(tài)類(lèi)型和雙穩(wěn)態(tài)類(lèi)型,分別如圖5(a)和圖5(b)所示.

圖5 懸浮磁體的單穩(wěn)態(tài)(a)和雙穩(wěn)態(tài)(b)現(xiàn)象Fig.5. The monostable(a)and bistable(b)phenomenon of the fl oating magnet.

4.1 單穩(wěn)態(tài)類(lèi)型時(shí)的動(dòng)力學(xué)響應(yīng)

文獻(xiàn)[17]給出了當(dāng)能量采集器處于非工作的單穩(wěn)態(tài)狀態(tài)時(shí)懸浮磁體振動(dòng)的幅頻關(guān)系方程:

當(dāng)能量采集器處于單穩(wěn)態(tài)類(lèi)型時(shí),懸浮磁體只有一個(gè)平衡位置.當(dāng)受到外界激勵(lì)時(shí),懸浮磁體在平衡位置附近往復(fù)運(yùn)動(dòng).選擇熱解石墨板間距L=7,7.4 mm研究單穩(wěn)態(tài)類(lèi)型時(shí)懸浮磁體的振動(dòng)特性.圖6是外界激勵(lì)幅值為1 mm、感應(yīng)線圈處于開(kāi)路狀態(tài)(非工作狀態(tài))時(shí)懸浮磁體的幅頻響應(yīng)曲線.可以發(fā)現(xiàn)懸浮磁體的幅頻響應(yīng)曲線向右偏斜,呈現(xiàn)“硬彈簧”特質(zhì),并且隨著熱解石墨板間距的增加,曲線向右偏斜程度增強(qiáng).當(dāng)熱解石墨板間距L為7 mm和7.4 mm時(shí),懸浮磁體振動(dòng)的基頻分別為2.37 Hz和1.5 Hz.增大熱解石墨板間距,懸浮磁體的振動(dòng)響應(yīng)頻帶變寬并出現(xiàn)多解和跳躍現(xiàn)象.

圖6 振幅隨激勵(lì)頻率變化的曲線Fig.6.The curve of amplitude changes with excited frequency.

圖7是熱解石墨板間距分別為7 mm和7.4 mm,感應(yīng)線圈工作時(shí),懸浮磁體的相對(duì)位移隨外界激振幅值和頻率的變化關(guān)系.當(dāng)外界激勵(lì)相同,增大熱解石墨板間距能獲得更大的相對(duì)運(yùn)動(dòng)位移.

圖7 單穩(wěn)態(tài)類(lèi)型時(shí)位移峰值隨外界激勵(lì)的變化Fig.7.Displacement peaks change with the external excitation in monostable type.

4.2 雙穩(wěn)態(tài)類(lèi)型時(shí)的動(dòng)力學(xué)響應(yīng)

當(dāng)采集器處于雙穩(wěn)態(tài)類(lèi)型時(shí),此時(shí)的動(dòng)力學(xué)方程的解Z(t)可以看作一勢(shì)子在勢(shì)場(chǎng)U(z)中的運(yùn)動(dòng),當(dāng)沒(méi)有外界激勵(lì)時(shí),U(z)=a1z2/2+a3z4/4,是一種雙穩(wěn)勢(shì)場(chǎng)結(jié)構(gòu),由dU(z)/dz=0很容易得出系統(tǒng)平衡點(diǎn)Zm1,Zm2,Zb位置,分別如(8)和(9)式所示:

當(dāng)受到外界周期激勵(lì)時(shí),系統(tǒng)的勢(shì)函數(shù)受到周期調(diào)制,由U(z)變?yōu)閂(z),如(10)式所示:

由(10)式可以看出,勢(shì)函數(shù)的勢(shì)阱被周期性地抬高或加深.此時(shí),系統(tǒng)存在一臨界幅值A(chǔ)c,使得方程保持雙穩(wěn)結(jié)構(gòu).當(dāng)外界周期激勵(lì)的幅值A(chǔ)＜Ac時(shí),勢(shì)阱周期性的抬高或加深不足以使勢(shì)壘消失,勢(shì)子無(wú)法越過(guò)勢(shì)壘進(jìn)入另一勢(shì)阱,只能在單阱中做小幅振蕩;當(dāng)外界周期激勵(lì)幅值A(chǔ)＞Ac時(shí),由于勢(shì)壘可能隨著時(shí)間的變化而消失,勢(shì)子能夠越過(guò)勢(shì)壘在兩勢(shì)阱間進(jìn)行大位移周期性跳躍.臨界幅值A(chǔ)c表示系統(tǒng)勢(shì)函數(shù)V(z)在Acos(2πft)為最大值或最小值時(shí)恰好從雙穩(wěn)態(tài)結(jié)構(gòu)變成單穩(wěn)結(jié)構(gòu).

圖8 雙穩(wěn)態(tài)類(lèi)型時(shí)懸浮永磁體共振分岔圖Fig.8.The resonant bifurcation diagram of fl oating magnets in bistable type.

兩熱解石墨板的間距不僅影響采集器的穩(wěn)態(tài)類(lèi)型同時(shí)也影響懸浮磁體的平衡位置.當(dāng)L=7.8 mm時(shí),懸浮永磁體在x=±0.51 mm處穩(wěn)定懸浮.動(dòng)力學(xué)方程的初始條件可以任選其一,此處選擇x=0.51 mm.感應(yīng)線圈開(kāi)路時(shí),懸浮磁體振動(dòng)幅值隨外界激勵(lì)的分岔圖,如圖8所示.由圖8可以看出,感應(yīng)線圈開(kāi)路時(shí),當(dāng)外界激勵(lì)頻率為1 Hz,激勵(lì)振幅從0.2—4 mm的變化過(guò)程中,系統(tǒng)出現(xiàn)了2倍周期—4倍周期—混沌—雙倍周期-單倍周期變化的現(xiàn)象.

表2是不同外界激振參數(shù)下,當(dāng)熱解石墨板間距L=7.8 mm時(shí)懸浮磁體的振動(dòng)相圖及其對(duì)應(yīng)的Poincare映射圖.當(dāng)非線性系統(tǒng)和線性系統(tǒng)進(jìn)行耦合時(shí),兩系統(tǒng)之間的強(qiáng)弱關(guān)系會(huì)直接影響到系統(tǒng)的振動(dòng)響應(yīng)特性:當(dāng)其中一個(gè)系統(tǒng)弱到整個(gè)系統(tǒng)可以看做兩個(gè)獨(dú)立的系統(tǒng),或者其中一個(gè)系統(tǒng)強(qiáng)到占統(tǒng)治地位時(shí),整個(gè)系統(tǒng)都將處于周期運(yùn)動(dòng)狀態(tài);只有當(dāng)兩系統(tǒng)的強(qiáng)弱不相上下時(shí),兩系統(tǒng)相互影響強(qiáng)烈,運(yùn)動(dòng)才變得十分復(fù)雜,從而才可能形成混沌[18].

圖9是兩熱解石墨板間距為7.8 mm、感應(yīng)線圈處于工作狀態(tài)、外界激勵(lì)頻率為1 Hz、激勵(lì)振幅從0.5 mm增大到8.5 mm的過(guò)程中(每次增加1 mm),懸浮磁體穩(wěn)態(tài)振動(dòng)的運(yùn)動(dòng)相圖.從圖9可以看出隨著激振幅值A(chǔ)的增加,倍周期現(xiàn)象和混沌現(xiàn)象消失,懸浮磁體進(jìn)行周期振動(dòng).

圖9 外界激勵(lì)頻率為1 Hz時(shí)雙穩(wěn)態(tài)類(lèi)型的運(yùn)動(dòng)相圖Fig.9.The vibration phase diagram when the excitation is 1 Hz in bistable type.

圖10是外界激勵(lì)頻率在0—3.2 Hz,激勵(lì)振幅在0—10 mm之間變化、系統(tǒng)處于雙穩(wěn)態(tài)類(lèi)型時(shí),懸浮磁體穩(wěn)態(tài)振動(dòng)的峰值.當(dāng)激勵(lì)振幅較小時(shí),懸浮磁體僅能在鞍點(diǎn)一側(cè)振動(dòng);當(dāng)激勵(lì)振幅足夠大時(shí),懸浮磁體跨過(guò)鞍點(diǎn)做阱間周期振動(dòng).

表2 雙穩(wěn)態(tài)類(lèi)型時(shí)懸浮磁體的運(yùn)動(dòng)相圖和對(duì)應(yīng)的Poincare映射圖Table 2.The vibration phase diagram and the Poincare mapping diagram of the fl oating magnetic in bistable type.

圖10 雙穩(wěn)態(tài)類(lèi)型時(shí)位移隨外界激勵(lì)的變化Fig.10.The displacement changes with the external excitation in bistable type.

5 總 結(jié)

抗磁懸浮振動(dòng)能量采集器中懸浮永磁體的運(yùn)動(dòng)特性直接決定著其輸出特性,懸浮磁體受力的非線性使得能量采集器的動(dòng)力學(xué)響應(yīng)和能量輸出特性較為復(fù)雜.本文利用COMSOLTM建立仿真模型,仿真獲得懸浮磁體所受的磁力和抗磁力,并進(jìn)行數(shù)據(jù)擬合得到磁力和抗磁力的解析表達(dá)式.然后,建立了懸浮永磁體的動(dòng)態(tài)響應(yīng)方程以研究其動(dòng)態(tài)響應(yīng)特性.在感應(yīng)線圈開(kāi)路(非工作狀態(tài))時(shí),單穩(wěn)態(tài)類(lèi)型的幅頻響應(yīng)曲線呈現(xiàn)向右偏斜的現(xiàn)象,并隨著熱解石墨板間距和非線性擾動(dòng)的增大而增強(qiáng).同時(shí)發(fā)現(xiàn)在相同激勵(lì)下,石墨板間距越小,采集器的響應(yīng)幅值越大.當(dāng)感應(yīng)線圈開(kāi)路且采集器處于雙穩(wěn)態(tài)類(lèi)型時(shí),通過(guò)分析懸浮磁體的振動(dòng)相圖及其對(duì)應(yīng)的Poincare映射圖發(fā)現(xiàn)采集器出現(xiàn)了倍周期、4倍周期和混沌等非線性系統(tǒng)特有的現(xiàn)象.當(dāng)感應(yīng)線圈處于工作狀態(tài)且采集器處于雙穩(wěn)態(tài)類(lèi)型時(shí)由于受到電磁阻尼力的影響,倍周期和混沌現(xiàn)象消失,懸浮磁體的振動(dòng)頻率和外界激勵(lì)頻率保持一致.該研究是對(duì)微型抗磁懸浮振動(dòng)能量采集器動(dòng)態(tài)特性的深入探索,對(duì)該能量采集器的結(jié)構(gòu)設(shè)計(jì)具有重要的指導(dǎo)意義,加快了該能量采集器的實(shí)用化進(jìn)程.基于本文的研究,可以對(duì)能量采集器進(jìn)行優(yōu)化,以提高其響應(yīng)特性和輸出性能.此外,本文也是對(duì)能量采集器中非線性振動(dòng)特性研究的補(bǔ)充和豐富,對(duì)理解和利用非線性振動(dòng)具有重要的意義.下一步的工作是研究分析該能量采集器的輸出特性,開(kāi)展動(dòng)態(tài)研究實(shí)驗(yàn),并對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比分析,推動(dòng)該能量采集器向?qū)嵱没~進(jìn).

[1]Ma X C,Ye R F,Zhang T L,Zhang X Q 2016Acta Phys.Sin.65 177701(in Chinese)[馬星晨,葉瑞豐,張?zhí)順?lè),張曉青2016物理學(xué)報(bào)65 177701]

[2]Liu L,Yuan F G 2013J.Sound.Vib.332 455

[3]Dai X Z,Wen Y M,Li P,Yang J,Jiang X F 2010Acta Phys.Sin.59 2137(in Chinese)[代顯智,文玉梅,李平,楊進(jìn),江小芳2010物理學(xué)報(bào)59 2137]

[4]Braunbek W 1939Z.Phys.112 764

[5]de Pasquale G,Siyanbalapitiya C,Somà A,Wang J 2009International Conference on Electromagnetics in Advanced ApplicationsTorino,Italy,September 14–18,2009 p465

[6]Cao J P,Yin D C,Qian A R,Tian Z C,Xu H Y,Huang Y P,Shang P 2011Chin.J.Space Sci.31 63(in Chinese)[曹建平,尹大川,騫愛(ài)榮,田宗成,續(xù)惠云,黃勇平,商澎2011空間科學(xué)學(xué)報(bào)31 63]

[7]Sun Y L,Chen Z H,Chen X H,Yin C,Li D J,Ma X L,Zhao F,Zhang G,Shang P,Qian A R 2015IEEE Trans.Bio-Med.Eng.62 900

[8]Billot M,Piat E,Abadie J,Stemp fl é P 2016Sensor.Actuat.A:Phys.238 266

[9]Hiber W,Jakoby B 2013IEEE Sens.J.13 2786

[10]Liu W,Chen W Y,Zhang W P,Huang X G 2008Electron.Lett.44 681

[11]Su Y F,Zhang K,Ye Z T,Zhang K P 2016Instrum.Tech.Sensor10 28(in Chinese)[蘇宇鋒,張坤,葉志通,張?chǎng)H鵬2016儀表技術(shù)與傳感器10 28]

[12]Ye Z T,Duan Z Y,Takahata K,Su Y F 2015Appl.Phys.A118 91

[13]Qin L Z,Su Y F,Liu W F 2016J.Electron.Meas.Instrum.30 1438(in Chinese)[秦立振,蘇宇鋒,劉武發(fā)2016電子測(cè)量與儀器學(xué)報(bào)30 1438]

[14]Ravaud R,Lemarquand G,Babic S,Lemarquand V,Akyel C 2010IEEE Trans.Magn.46 3585

[15]Su Y F,Ye Z T,Xiao Z M,Takahata K 2015IEEE 10th International Conference on NEMSXi’an China,July 7–11,2015 p116

[16]Liu K,Zhang W P,Liu W,Chen W Y,Li K,Cui F,Li S P 2010Microsyst.Technol.16 431

[17]Xu W 2001Numerical Analysis Methods for Stochastic Dynamical System(Beijing:Science Press)p179(in Chinese)[徐偉 2001非線性隨機(jī)動(dòng)力學(xué)的若干數(shù)值方法和應(yīng)用(北京:科學(xué)出版社)第179頁(yè)]

[18]Lai Z H 2014Ph.D.Dissertation(Tianjin:Tianjin University)(in Chinese)[賴(lài)志慧 2014博士論文 (天津:天津大學(xué))]