楊海峰， 吳子燕， 孫 蓬， 張 坤

(1.西北工業(yè)大學(xué)力學(xué)與土木建筑學(xué)院 西安，710129) (2.西北工業(yè)大學(xué)計(jì)算機(jī)學(xué)院 西安，710129)

基于逆有限元的應(yīng)變模態(tài)損傷檢測方法*

楊海峰1， 吳子燕1， 孫 蓬2， 張 坤1

(1.西北工業(yè)大學(xué)力學(xué)與土木建筑學(xué)院 西安，710129) (2.西北工業(yè)大學(xué)計(jì)算機(jī)學(xué)院 西安，710129)

基于光纖光柵傳感器(fiber Bragg grating,簡稱FBG)逆有限元方法(inverse finite element method,簡稱iFEM)，僅利用有限測點(diǎn)的應(yīng)變數(shù)據(jù)進(jìn)行全域應(yīng)變場重構(gòu)，得到近似完全測量應(yīng)變模態(tài)，提高了直接采用實(shí)測應(yīng)變數(shù)據(jù)來構(gòu)建應(yīng)變模態(tài)損傷指標(biāo)的實(shí)用性。利用基于損傷應(yīng)變模態(tài)差分原理的損傷指標(biāo)法，只需用損傷后應(yīng)變模態(tài)數(shù)據(jù)即能定位損傷，并給出了損傷指標(biāo)數(shù)學(xué)模型。計(jì)算結(jié)果表明，基于光纖光柵傳感器和逆有限元方法可以快速進(jìn)行全域應(yīng)變場重構(gòu)，為基于應(yīng)變模態(tài)的損傷檢測提供數(shù)據(jù)保障，而應(yīng)變模態(tài)差分曲線只在損傷處發(fā)生劇烈變化，損傷程度不同，曲線突變程度略有不同但規(guī)律一致。最后以某板的實(shí)驗(yàn)結(jié)果證明了該方法的有效性。

逆有限元法；損傷檢測；應(yīng)變模態(tài)；光纖光柵傳感器

引 言

結(jié)構(gòu)損傷識(shí)別作為結(jié)構(gòu)健康監(jiān)測系統(tǒng)的核心內(nèi)容，在航空航天、機(jī)械和土木工程等領(lǐng)域得到應(yīng)用。由于實(shí)際結(jié)構(gòu)的早期微小損傷不容易被檢測出來，會(huì)對結(jié)構(gòu)的安全性與可靠性產(chǎn)生影響，甚至引發(fā)重大事故。通常結(jié)構(gòu)損傷會(huì)引起模態(tài)參數(shù)的變化，國內(nèi)外學(xué)者提出了很多識(shí)別損傷的指標(biāo)，如固有頻率、振型、位移類模態(tài)參數(shù)(模態(tài)柔度、位移、位移模態(tài))和應(yīng)變類模態(tài)參數(shù)(曲率模態(tài)、應(yīng)變、應(yīng)變模態(tài))等[1-4]。

已有大量研究表明，固有頻率僅能識(shí)別損傷是否存在,而不能識(shí)別其位置和程度。模態(tài)柔度對單損傷有較好的識(shí)別能力，但對多損傷效果不明顯。位移模態(tài)僅能識(shí)別損傷位置和總體的損傷程度，但不能識(shí)別局部損傷程度。應(yīng)變模態(tài)對損傷識(shí)別優(yōu)于位移模態(tài)，特別是在結(jié)構(gòu)局部應(yīng)力集中或裂紋檢測時(shí)，它能獲得結(jié)構(gòu)的動(dòng)應(yīng)變，且應(yīng)變值只在損傷附近變化明顯，遠(yuǎn)離損傷處改變不大。

采用應(yīng)變模態(tài)至少有兩方面的好處，其一是可免去由位移到應(yīng)變計(jì)算過程中所帶來的誤差，而這種誤差往往是很難控制的。因?yàn)橛晌灰频綉?yīng)變是微分過程，位移的微小改變將被放大，從而引起應(yīng)變參數(shù)的顯著變化。因此，應(yīng)變類參數(shù)(應(yīng)變模態(tài)和曲率模態(tài))對結(jié)構(gòu)局部微小損傷反應(yīng)較敏感，而位移模態(tài)對此反應(yīng)甚微。其二是可直接研究某些關(guān)鍵部位的應(yīng)變，如應(yīng)力集中問題和局部結(jié)構(gòu)變動(dòng)對附近受力情況的影響等。

應(yīng)變類損傷指標(biāo)大多需要利用損傷前后的模態(tài)數(shù)據(jù)。實(shí)際上，由于設(shè)計(jì)資料不全、施工離散性和材料劣化等問題，要得到損傷前的精確狀態(tài)非常困難。因此應(yīng)變類損傷指標(biāo)法在實(shí)際工程中尚難以推廣應(yīng)用[5]。為此，筆者利用光纖光柵傳感器和逆有限元理論，提出基于逆有限元的應(yīng)變模態(tài)損傷檢測方法，僅利用有限測點(diǎn)的應(yīng)變數(shù)據(jù)且無需損傷前的模態(tài)數(shù)據(jù)即可進(jìn)行損傷檢測。

1 基本理論

1.1 應(yīng)變模態(tài)基本理論

應(yīng)變是位移的一階導(dǎo)數(shù),對應(yīng)于每一階位移模態(tài)，必有其對應(yīng)的固有應(yīng)變分布狀態(tài)，這種與位移模態(tài)相對應(yīng)的固有應(yīng)變分布狀態(tài)稱之為應(yīng)變模態(tài)，它能反映結(jié)構(gòu)的固有特征。設(shè)彈性結(jié)構(gòu)的變形位移為u，ν和w，運(yùn)用模態(tài)疊加進(jìn)行應(yīng)變?chǔ)舩推導(dǎo)

(1)

(2)

(3)

位移響應(yīng)表示為

(4)

根據(jù)應(yīng)力應(yīng)變關(guān)系，得到應(yīng)變響應(yīng)表達(dá)式為

(5)

其中：εx，εy，εz為正應(yīng)變；γxy，γyz，γzx為剪應(yīng)變。

通常剪切應(yīng)變無法直接測量，筆者將光纖光柵傳感器按照應(yīng)變花的形式布置，如圖1所示，再利用平面應(yīng)變理論推導(dǎo)出剪切應(yīng)變。

圖1 光纖光柵傳感器布置形式Fig. 1 Placement of fiber Bragg grating sensor

由平面應(yīng)變狀態(tài)幾何關(guān)系分析可知

其中：εα和γα分別為α截面上的線(正)應(yīng)變和切應(yīng)變。

將α=45°代入式(6)，可得

(8)

γxy=εx+εy-2ε45°

(9)

按筆者方式貼放光纖傳感器并結(jié)合式(9)，即可獲得待測點(diǎn)的應(yīng)變?chǔ)舩,εy和γxy。

1.2 損傷應(yīng)變模態(tài)差分?jǐn)?shù)學(xué)模型

(10)

其中：xi為節(jié)點(diǎn)i的坐標(biāo)值。

設(shè)xi+1-xi=h，并假定h充分小，不計(jì)三次冪項(xiàng)及更高次冪項(xiàng)，等間距差分式為

(11)

曲線的變化程度通常利用導(dǎo)數(shù)反映，但在未歸一化之前，各階應(yīng)變模態(tài)差分曲線數(shù)量級存在差異，沒有可比性，此外還存在試驗(yàn)或計(jì)算誤差等引起的突變假象問題，導(dǎo)致直接從差分曲線進(jìn)行突變程度定量判定較為困難。為此引入有效極值點(diǎn)概念[6]：若差分曲線任意相鄰兩極值符號相異，則它們之間有且僅有一個(gè)零值點(diǎn)，該極值點(diǎn)稱為有效極值點(diǎn)。

基于有效極值點(diǎn)定義的損傷位置指標(biāo)數(shù)學(xué)模型為

(12)

其中：ID(j)為某階應(yīng)變模態(tài)差分曲線第j個(gè)有效零值點(diǎn)損傷位置指標(biāo)值；xi,yi分別為差分曲線有效極值點(diǎn)與相應(yīng)的極值；l為總體跨度尺寸。

該指標(biāo)具有明確的位置坐標(biāo)，且ID(j)值越大，損傷可能性越大，再結(jié)合差分曲線變化規(guī)律，即可判定損傷位置。

1.3 基于光纖光柵傳感器和逆有限元的應(yīng)變模態(tài)損傷識(shí)別

將應(yīng)變模態(tài)指標(biāo)應(yīng)用于實(shí)際結(jié)構(gòu)損傷檢測，首先要解決測量手段問題。和傳統(tǒng)的應(yīng)變片相比，光纖光柵對應(yīng)力應(yīng)變變化非常敏感，測量精度也非常高，應(yīng)變測量分辨率已達(dá)1 με，利用光纖光柵傳感器復(fù)用技術(shù)還可同時(shí)測量多點(diǎn)應(yīng)變，甚至埋入被測構(gòu)件內(nèi)部，實(shí)現(xiàn)構(gòu)件內(nèi)部應(yīng)力應(yīng)變測量[7-8]。

但若想獲取未布置傳感器區(qū)域的應(yīng)力應(yīng)變信息，以及區(qū)域內(nèi)任一點(diǎn)的應(yīng)力應(yīng)變，要根據(jù)已知測點(diǎn)數(shù)據(jù)進(jìn)行結(jié)構(gòu)全域應(yīng)力應(yīng)變反演。在這個(gè)過程中，逆有限元法被證明是一種比較有效的方法。

離散應(yīng)變測量值重構(gòu)結(jié)構(gòu)中每個(gè)結(jié)點(diǎn)的應(yīng)變和位移屬于逆數(shù)學(xué)問題。在不必要滿足存在性、唯一性和穩(wěn)定性的意義上，逆問題是不適定的。對于這類問題，通常的解法是基于正則化算子構(gòu)建近似解，但此近似解只有在測量數(shù)據(jù)變化小的情況下是穩(wěn)定的。Tessler等[9-11]基于有限元思想，結(jié)合變分原理，提出逆有限元法，其本質(zhì)是首先根據(jù)傳感器布局，對結(jié)構(gòu)進(jìn)行合理的有限元網(wǎng)格劃分(劃分方式與傳感器布局有關(guān)，筆者采用均布方式)，然后利用實(shí)驗(yàn)測得的結(jié)構(gòu)表面應(yīng)變，選取合適的有限元形函數(shù)對單元進(jìn)行反演推導(dǎo)，再將單元方程擴(kuò)展至結(jié)構(gòu)整體方程，由此實(shí)現(xiàn)僅利用有限測點(diǎn)數(shù)據(jù)進(jìn)行結(jié)構(gòu)全域應(yīng)力應(yīng)變反演。其變分公式基于對最小二乘泛函求最小值，且該泛函僅使用符合一階線性剪切變形理論的應(yīng)變測量數(shù)據(jù)。由于該公式容許橫向剪切變形，因此適用于薄的及適中厚度的板殼結(jié)構(gòu)。變分原理的主要好處是適用于具有C0連續(xù)性的有限單元離散化，因此可以將魯棒性算法應(yīng)用到復(fù)雜的民用和航空結(jié)構(gòu)，特別適合針對下一代實(shí)時(shí)結(jié)構(gòu)健康監(jiān)測系統(tǒng)的使用。

運(yùn)用逆有限元法進(jìn)行結(jié)構(gòu)全域應(yīng)力應(yīng)變場反演的主要理論推導(dǎo)過程如下。

首先，使用由Mindlin理論應(yīng)變值和對應(yīng)的應(yīng)變測量值構(gòu)成的最小二乘差分項(xiàng)構(gòu)建誤差泛函，其中Mindlin應(yīng)變值以位移的形式表示。受懲罰參數(shù)控制的正則化項(xiàng)強(qiáng)迫其服從橫向剪切應(yīng)變的約束條件，由于這個(gè)設(shè)定，使所有應(yīng)變協(xié)調(diào)關(guān)系自然明確地得到滿足，所以該方法與材料的彈性和慣性特性無關(guān)[9]。

平板殼單元的變形按Mindlin理論中的位移向量定義為

(13)

其中：u=u(x,y),v=v(x,y)分別為中平面x和y方向上的位移；θx=θx(x,y)，θy=θy(x,y)分別為以x軸負(fù)方向和y軸正方向?yàn)榉ň€的轉(zhuǎn)角；w=w(x,y)為撓度變量，其沿厚度坐標(biāo)z∈[-t,t]恒定不變；2t為板殼的總厚度。

應(yīng)變位移關(guān)系式為

(14)

與中平面拉伸相關(guān)聯(lián)的膜應(yīng)變?yōu)?/p>

Leu

(15)

同理，彎曲曲率及橫向剪切可分別寫為k(u)≡Lku和g(u)≡Lgu。

圖2 使用應(yīng)變花測量板的上下表面應(yīng)變Fig.2 Measurement of strain of using FBG strain flower

結(jié)構(gòu)單元位移向量由最小二乘泛函取極值時(shí)獲得，可表示為

(16)

利用位移有限單元法，通過使用適當(dāng)?shù)男魏瘮?shù)插值式(16)中的位移場，可得

(17)

其中：N為具有C0連續(xù)性的非等參形函數(shù)矩陣。

Keue=fe

(18)

其中：Ke僅為位置坐標(biāo)的函數(shù)；fe為坐標(biāo)和實(shí)測應(yīng)變值的函數(shù)。

最后，通過標(biāo)準(zhǔn)有限單元集成操作，將離散的單元矩陣組裝到整體線性方程系統(tǒng)中，即可實(shí)現(xiàn)位移自由度從單元局部坐標(biāo)系到總體坐標(biāo)系的轉(zhuǎn)換。通過施加Mindlin位移邊界條件，位移自由度解存在且唯一。一旦得到結(jié)點(diǎn)位移，可通過單元級別的簡單計(jì)算得到平滑的單元應(yīng)變。這些計(jì)算都基于線性方程，因此計(jì)算效率很高。

運(yùn)用逆有限元方法進(jìn)行全域反演之后，可獲取近似完全測量的應(yīng)變模態(tài)，解決了傳統(tǒng)應(yīng)變模態(tài)損傷指標(biāo)需要完備的實(shí)測振型問題，大大提高直接采用測量數(shù)據(jù)構(gòu)建應(yīng)變模態(tài)損傷指標(biāo)進(jìn)行損傷檢測的實(shí)用性。

2 算例及分析

以圖3所示的懸臂板為例，對筆者提出的算法進(jìn)行驗(yàn)證。板尺寸為65 cm×30 cm×0.68 cm，左端采用三點(diǎn)固定，材料彈性模量E=210 GPa，泊松比ν=0.3，測試儀器包括光纖傳感器、光纖光柵波長解調(diào)儀等。損傷設(shè)置情況為在兩個(gè)單元(1157號、809號)處的單元?jiǎng)偠确謩e降低20%和50%(如圖4所示)，其位置(x,y)坐標(biāo)分別為(51 cm×17 cm，28 cm×12 cm)。在距右端5 cm處(距上邊緣10 cm)施加集中力P，采用懸掛標(biāo)準(zhǔn)砝碼的方式施加激勵(lì)，模擬集中作用力(實(shí)驗(yàn)所用的板右側(cè)有一微孔可以吊裝標(biāo)準(zhǔn)砝碼荷載)。

圖3 構(gòu)件尺寸圖(單位：cm)Fig.3 Size chart of plate(unit:cm)

圖4 實(shí)驗(yàn)現(xiàn)場與有限元模型圖Fig.4 Experimental field and finite element model (FEM)

首先，利用布置在板上的光纖傳感器(每個(gè)測點(diǎn)為一組應(yīng)變花)測得18個(gè)點(diǎn)的應(yīng)變值；然后，利用逆有限元法進(jìn)行結(jié)構(gòu)全域應(yīng)力應(yīng)變反演與重構(gòu)，以獲取近似完全測量的應(yīng)變模態(tài)，逆有限元法反演重構(gòu)的撓度圖以及損傷后的各階應(yīng)變模態(tài)如圖5，6所示。

從圖5可以看出，利用逆有限元法對懸臂板撓度的反演重構(gòu)結(jié)果與有限元仿真軟件計(jì)算的撓度結(jié)果非常接近(從數(shù)據(jù)上對比顯示平均相對誤差僅為1.73%，因數(shù)據(jù)量過多，此處未列出)，表明逆有限元法在全域反演方面具有非常好的效果。

通過逆有限元法反演可得到近似完全測量的應(yīng)變模態(tài)，這是普通應(yīng)變測量無法做到的。利用反演后的損傷應(yīng)變模態(tài)數(shù)據(jù)計(jì)算三維應(yīng)變差分結(jié)果及應(yīng)變模態(tài)差分曲線結(jié)果如圖7，8，9所示。應(yīng)變模態(tài)差分曲線及損傷位置指標(biāo)為曲線形式，因此只取經(jīng)過損傷點(diǎn)的二維應(yīng)變模態(tài)差分曲線上的數(shù)據(jù)進(jìn)行計(jì)算。

圖5 逆有限元法反演重構(gòu)撓度與有限元軟件數(shù)值分析結(jié)果Fig. 5 Deflection results of iFEM method and FEM software

圖6 損傷后應(yīng)變模態(tài)Fig. 6 Strain model shape before and after damage

圖7 五階應(yīng)變模態(tài)的三維應(yīng)變差分圖Fig. 7 3D Strain difference diagram of order 5

圖8 第1階應(yīng)變模態(tài)差分曲線Fig.8 Strain modal difference curve of first order

圖9 第3階應(yīng)變模態(tài)差分曲線Fig.9 Strain modal difference curve of third order

類別損傷程度20%50%1階ID(1)18.457.60ID(2)26.8232.21ID(3)18.567.3ID(4)3.525.253階ID(2)9.1110.81ID(3)30.1743.52ID(4)10.1211.33ID(5)22.866.23

黑體數(shù)字對應(yīng)損傷位置

從圖6可以看出，僅僅從應(yīng)變模態(tài)圖中很難發(fā)現(xiàn)損傷情況(應(yīng)變模態(tài)圖中左邊界處的3個(gè)峰值是由于板的左端采用三點(diǎn)固定夾持方式導(dǎo)致)，在假設(shè)能夠獲取損傷前模態(tài)數(shù)據(jù)的前提下或許還可以利用損傷前后的數(shù)據(jù)對比分析來判斷損傷情況，但在無法獲取損傷前數(shù)據(jù)時(shí)，僅根據(jù)模態(tài)數(shù)據(jù)來分析損傷將非常困難。應(yīng)變差分圖中，雖然一階模態(tài)差分曲線并不能明顯看出損傷，但在三階模態(tài)差分曲線中可以發(fā)現(xiàn)在損傷發(fā)生的位置曲線有明顯的突變(由于噪聲干擾及數(shù)據(jù)降噪處理后，20%損傷處的突變較小)。這也說明無論是單一損傷還是多個(gè)損傷，都可以在應(yīng)變模態(tài)差分曲線中得到體現(xiàn)，表明筆者提出方法的有效性。

由圖8，9還可知，在歸一化之前，各階應(yīng)變模態(tài)差分?jǐn)?shù)據(jù)的數(shù)量級相差懸殊，無法統(tǒng)一對比，只能借助歸一化之后的損傷指標(biāo)ID(j)進(jìn)行分析。若某點(diǎn)在某階或多階ID值越大，損傷可能性越大。實(shí)際計(jì)算結(jié)果如表1所示，可以看出，各階的損傷處ID值均為最大值(1階的ID(2)和3階的ID(3)，分別對應(yīng)著剛度降低20%和50%的兩處損傷位置)。節(jié)點(diǎn)處損傷可通過其余階指標(biāo)值ID(j)定位，再結(jié)合差分曲線變化規(guī)律判定出損傷位置。

3 結(jié)束語

基于光纖光柵傳感器和逆有限元方法，利用有限測點(diǎn)的應(yīng)變數(shù)據(jù)進(jìn)行全域應(yīng)變場重構(gòu)，得到近似完全測量應(yīng)變模態(tài)，解決了傳統(tǒng)應(yīng)變模態(tài)損傷指標(biāo)需要完備的實(shí)測振型問題。利用基于損傷應(yīng)變模態(tài)差分原理的損傷指標(biāo)法，僅用損傷后的數(shù)據(jù)進(jìn)行了損傷定位。計(jì)算結(jié)果表明，逆有限元法是根據(jù)已知測點(diǎn)數(shù)據(jù)進(jìn)行結(jié)構(gòu)全域應(yīng)力應(yīng)變反演比較有效的一種方法，通過逆有限元全域反演可得到近似完全測量的應(yīng)變模態(tài)，這是普通應(yīng)變測量無法做到的。還發(fā)現(xiàn)除支座附近外，應(yīng)變模態(tài)差分曲線只在損傷處發(fā)生劇烈變化，損傷程度不同，突變程度略有不同但規(guī)律一致。懸臂板的實(shí)驗(yàn)結(jié)果也證明所提方法的有效性。

[1] Yan Yunju, Cheng Li, Wu Ziyan, et al. Development in vibration-based structural damage detection technique [J]. Mechanical Systems and Signal Processing, 2007, 21(5): 2198-2211.

[2] Wu H. The challenges of structural health monitoring technologies in civil infrastructure [C]∥ 8th International Workshop on Structural Health Monitoring． Stanford, California：DEStech Publications, Inc., 2011.

[3] Pandey A K, Biswas M. Damage diagnosis of truss structures by estimation of flexibility change [J]. Modal Analysis, 1995, 10(2): 104-117.

[4] 周邵萍,郝占峰,韓紅飛,等.基于應(yīng)變模態(tài)差和神經(jīng)網(wǎng)絡(luò)的管道損傷識(shí)別[J]. 振動(dòng)、測試與診斷, 2015, 35(2): 334-338.

Zhou Shaoping, Hao Zhanfeng, Han Hongfei, et al. Damage identification in straight pipeline using strain modal difference and neural network [J]. Journal of Vibration, Measurement & Diagnosis，2015,35(2):334-338. (in Chinese)

[5] 顧培英, 陳厚群, 李同春,等. 用應(yīng)變模態(tài)技術(shù)診斷梁結(jié)構(gòu)的損傷[J].地震工程與工程振動(dòng),2005,25(4):50-53.

Gu Peiying, Chen Houqun, Li Tongchun， et al. Application of strain mode in diagnosing damage to beam [J]. Earthquake Engineering and Engineering Vibration, 2005,25(4): 50-53.(in Chinese)

[6] 顧培英. 基于應(yīng)變模態(tài)技術(shù)的結(jié)構(gòu)損傷診斷直接指標(biāo)法研究[D]．南京：河海大學(xué)，2006.

[7] 姜德生, 羅裴, 梁磊. 光纖布拉格光柵傳感器與基于應(yīng)變模態(tài)理論的結(jié)構(gòu)損傷識(shí)別[J]. 儀表技術(shù)與傳感器, 2003(2):17-19.

Jiang Desheng, Luo Pei, Liang Lei. Fiber Bragg grating sensor and recognition of structural damage based on strain modal theory [J]. Instrument Technique and Sensor, 2003(2): 17-19.(in Chinese)

[8] 柯孟龍. 基于光纖光柵傳感和應(yīng)變動(dòng)力參數(shù)的損傷識(shí)別研究[D]. 武漢：武漢理工大學(xué)，2011.

[9] Tessler A. Structural analysis methods for structural health management of future aerospace vehicles[R].2007,NASA/TM-2007-214871.https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20070018347.pdf.

[10] Tondreau G, Deraemaeker A. Automated data-based damage localization under ambient vibration using local modal filters and dynamic strain measurements: experimental applications [J]. Journal of Sound and Vibration, 2014, 333(26): 7364-7385.

[11] Gherlone M, Cerracchio P, Mattone M, et al. Shape sensing of 3D frame structures using an inverse finite element method [J]. International Journal of Solids and Structures, 2012 (49): 3100-3112.

10.16450/j.cnki.issn.1004-6801.2017.01.023

*國家自然科學(xué)基金資助項(xiàng)目(51278420);中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金資助項(xiàng)目(3102015BJ(II)MYZ13); 西北工業(yè)大學(xué)研究生創(chuàng)意創(chuàng)新種子基金資助項(xiàng)目(Z2016094);西北工業(yè)大學(xué)基礎(chǔ)研究基金資助項(xiàng)目(JCY20130125)

2015-11-23；

2016-05-31

O327;TH17

楊海峰，男，1976年9月生，博士。主要研究方向?yàn)榻Y(jié)構(gòu)健康監(jiān)測與結(jié)構(gòu)損傷檢測。曾發(fā)表《基于MIGA算法的結(jié)構(gòu)模型修正研究及其應(yīng)用》(《振動(dòng)、測試與診斷》2012年第32卷第2期)等論文。 E-mail: hfyang@nwpu.edu.cn