葉片裂紋故障早期監(jiān)測(cè)預(yù)警方法研究?

張 婭 陳 康 王維民 張旭龍

（北京化工大學(xué)機(jī)電工程學(xué)院）

葉片裂紋故障早期監(jiān)測(cè)預(yù)警方法研究?

張 婭 陳 康 王維民 張旭龍

（北京化工大學(xué)機(jī)電工程學(xué)院）

振動(dòng)問題一直是限制風(fēng)機(jī)向大型化高參數(shù)化方向發(fā)展的主要障礙，而葉片振動(dòng)導(dǎo)致的疲勞斷裂故障和事故尤為突出。對(duì)風(fēng)機(jī)葉片同步共振進(jìn)行了研究，提出基于非線性最小二乘擬合和GARIV方法相結(jié)合的葉片同步振動(dòng)的振動(dòng)參數(shù)精確辨識(shí)方法，準(zhǔn)確捕捉葉片的動(dòng)態(tài)固有頻率。基于激光傳感器的葉尖定時(shí)方法，通過捕捉葉片裂紋早期固有頻率的微小變化，實(shí)現(xiàn)葉片裂紋故障的早期預(yù)警。仿真結(jié)果表明，所提出的方法可以達(dá)到萬分之七的測(cè)量精度，實(shí)驗(yàn)研究結(jié)果表明，本文所提出的方法可以進(jìn)行葉片裂紋的早期識(shí)別。

葉片振動(dòng)；參數(shù)識(shí)別；模擬仿真；實(shí)驗(yàn)測(cè)試；葉尖定時(shí)；固有頻率

0 引言

葉片是通風(fēng)機(jī)、鼓風(fēng)機(jī)和透平壓縮機(jī)等旋轉(zhuǎn)設(shè)備的核心部件之一，其工作安全與工作效率對(duì)風(fēng)機(jī)的平穩(wěn)運(yùn)行意義重大[1-5]。振動(dòng)問題一直是限制風(fēng)機(jī)向大型高參數(shù)化方向發(fā)展的主要障礙，而葉片振動(dòng)導(dǎo)致的疲勞斷裂故障和事故尤為突出。

葉片振動(dòng)測(cè)量主要分為接觸式和非接觸式兩種方法[6]，葉尖計(jì)時(shí)法是目前非接觸式葉片測(cè)振的主要方法。目前國際上公開且有效的基于葉尖定時(shí)法原理的葉片振動(dòng)參數(shù)識(shí)別算法主要有單參數(shù)法[6]（速矢端跡法）、雙參數(shù)法[7]（橢圓擬合法）、正弦擬合法[8]（包括三參數(shù)、四參數(shù)、六參數(shù)法）、自回歸法[9]（即AR法以及基于其發(fā)展而來的GAR法及GARIV法）、Number of Aliases法[8]、基于速矢端跡法的倍頻遍歷法[10]。

各個(gè)方法各有優(yōu)缺點(diǎn)，如雙參數(shù)法對(duì)信號(hào)質(zhì)量要求較高，倍頻遍歷法需要較多的傳感器，速矢端跡法無法求出倍頻值等。因此需要新的方法來更好地進(jìn)行葉片振動(dòng)參數(shù)的識(shí)別。

1 葉片振動(dòng)頻率辨識(shí)方法

本課題組提出的結(jié)合速矢端跡擬合法和輔助參數(shù)全局自回歸法（GARIV法）的葉片參數(shù)識(shí)別方法。該方法具體為先通過掃頻擬合確定同步振動(dòng)中諧共振點(diǎn)處的中心頻率、幅值等參數(shù)；然后對(duì)存在的諧共振中心進(jìn)行恒速運(yùn)轉(zhuǎn)，通過GARIV方法獲得對(duì)應(yīng)諧共振中心的倍頻值，進(jìn)而獲得較為完整的葉片同步振動(dòng)參數(shù)。本方法的提出可在傳感器使用的數(shù)量上減少到4支，安裝的位置較為緊密，傳感器夾角較小，且可以獲得較高的參數(shù)辨識(shí)精度。

1.1 葉片振動(dòng)倍頻值辨識(shí)方法

根據(jù)自回歸法原理[9]，求解的方程中含有兩個(gè)未知數(shù)，a1=2cos(ωnΔtp)及振動(dòng)恒偏直流分量xDC，求解需要至少4支葉頂傳感器，且傳感器的安裝時(shí)間隔夾角相等。假設(shè)4支00葉頂傳感器監(jiān)測(cè)得到的振動(dòng)位移分別為x1，x2，x3，x4?？捎蟹匠探M如下：

通過上式可以解出a1值的大小，結(jié)合a1=2cos(ωnΔtp)可以獲得葉片振動(dòng)頻率ωn值。葉片同步振動(dòng)過程中諧共振中心對(duì)應(yīng)的倍頻值可表示為nEO=ωn/Ω：其中Ω運(yùn)轉(zhuǎn)轉(zhuǎn)頻值（單位：Hz）。將 a1=2cos(ωnΔtp)帶入式nEO=ωn/Ω中可得：

當(dāng)傳感器等間距分布之間的夾角為α?xí)r（與下文中的傳感器安裝夾角相同），則Δtp表示為：

帶入式（2）可得：

推導(dǎo)中假設(shè)葉片組運(yùn)轉(zhuǎn)轉(zhuǎn)速為諧共振中心處轉(zhuǎn)速，但通過分析式（4）可知運(yùn)轉(zhuǎn)轉(zhuǎn)速對(duì)振動(dòng)倍頻的求解沒有影響。因此可以認(rèn)為只要a1值的求解準(zhǔn)確，則在任何運(yùn)轉(zhuǎn)轉(zhuǎn)速下都可以準(zhǔn)確求解出振動(dòng)倍頻值。

GARIV法的倍頻求取方法與以上基本相同，故不再詳細(xì)敘述。

1.2 固有頻率的精確識(shí)別方法

基于掃頻擬合和GARIV法的葉片同步振動(dòng)參數(shù)辨識(shí)流程如圖1所示。

圖1 葉片同步振動(dòng)參數(shù)辨識(shí)流程圖Fig.1 Flow chart of blade synchronous vibration parameter identification

由于葉片實(shí)際的振動(dòng)是多倍頻的，在掃頻過程中可能監(jiān)測(cè)到較多個(gè)明顯的同步振動(dòng)區(qū)間，對(duì)各傳感器采樣信號(hào)進(jìn)行同步振動(dòng)區(qū)間定位，截取不同傳感器采集到的同一編號(hào)葉片的同一諧共振區(qū)間分別進(jìn)行最小二乘法擬合，獲得葉片的大部分同步振動(dòng)參數(shù)，包括諧共振中心頻率、振動(dòng)幅值、振動(dòng)恒偏量等。

然后根據(jù)存在的諧共振中心頻率選取合適的測(cè)試轉(zhuǎn)速進(jìn)行恒速運(yùn)轉(zhuǎn)測(cè)試，4支周向均布的葉間計(jì)時(shí)傳感器采集的振動(dòng)位移信息經(jīng)過GARIV法分析得到振動(dòng)倍頻值，進(jìn)而確定存在的不同諧共振區(qū)對(duì)應(yīng)的振動(dòng)倍頻值。在實(shí)驗(yàn)過程中電機(jī)轉(zhuǎn)速會(huì)存在一定的轉(zhuǎn)速波動(dòng)，當(dāng)運(yùn)轉(zhuǎn)頻率與諧共振中心頻率一致時(shí)，葉片振動(dòng)位移可能會(huì)隨轉(zhuǎn)速波動(dòng)。故在實(shí)際測(cè)量中選擇在諧共振中心頻率值附近選取運(yùn)轉(zhuǎn)頻率。

最后可根據(jù)擬合獲得的諧共振中心頻率、振動(dòng)幅值等參數(shù)，結(jié)合諧共振對(duì)應(yīng)的振動(dòng)倍頻值得到葉片的固有頻率，做出葉片振動(dòng)坎貝爾圖。

固有頻率計(jì)算公式如下：

其中,fn為倍頻值nEO下對(duì)應(yīng)的共振中心頻率。

1.3 模擬仿真驗(yàn)證

為了進(jìn)行下一步的實(shí)驗(yàn)驗(yàn)證，搭建了相關(guān)實(shí)驗(yàn)臺(tái)，實(shí)驗(yàn)臺(tái)的介紹將在后文進(jìn)行介紹，實(shí)驗(yàn)臺(tái)葉輪的剖視圖如圖2（a）所示，該葉輪通過錐面定位安裝在電機(jī)軸上。為了降低葉片的固有頻率從而使得在更低的轉(zhuǎn)速下激起高階的固有頻率，在葉片根部進(jìn)行葉片局部減薄，結(jié)構(gòu)如圖2（b），最小厚度1.5mm。輪盤俯視圖如圖2（c）所示，整個(gè)輪盤由32個(gè)葉片組成，輪盤直徑138mm，葉尖厚度3mm。在輪盤上部開有鍵相槽和為了平衡鍵相而設(shè)計(jì)的平衡槽。根據(jù)傳感器安裝孔位置和運(yùn)轉(zhuǎn)方向，將32個(gè)葉片進(jìn)行了編號(hào)區(qū)分，具體可參見圖2（a）中紅色數(shù)字表示。

圖2 葉輪結(jié)構(gòu)圖及葉片編號(hào)圖Fig.2 Impeller structure and blade number

用ANSYS軟件對(duì)其模態(tài)進(jìn)行分析，得到前四階頻率及模態(tài)，及得到葉片部分振動(dòng)坎貝爾圖見圖3和圖4。

圖3 前4階模態(tài)振型圖Fig.3 Former fourth order mode

圖4 高速直葉片振動(dòng)坎貝爾圖Fig.4 Campbell diagram of high-speed straight blade test bench

綜合以上，利用文獻(xiàn)[11]所述Simulink模型建立過程進(jìn)行算法的模擬仿真。采用4葉片模型，葉片間無耦合，但葉輪為失諧葉輪，各個(gè)葉片的固有頻率有一定差異。

設(shè)置葉片參數(shù)見表1。

表1 葉片參數(shù)設(shè)置表Tab.1 Blade parameter settings

假設(shè)在運(yùn)轉(zhuǎn)過程中葉片的30倍頻共振被激起，即共振中心頻率在76.17Hz附近，每個(gè)葉片受到的激振力函數(shù)可以表示為：

式中,Fb表示激振力幅值常量;Ω表示轉(zhuǎn)子系統(tǒng)運(yùn)轉(zhuǎn)轉(zhuǎn)速值;b表示葉片編號(hào);φ表示葉片b的延遲相位，為方便計(jì)算通常取φ=0,設(shè)置Fb=0.1N。

1）掃頻模擬仿真

設(shè)置掃頻范圍為60到90Hz升速，升速時(shí)間為30s，升速速率為1Hz/s（60rpm/s）。采用4支傳感器，夾角6°均布。

通過4支傳感器對(duì)4葉片系統(tǒng)進(jìn)行葉尖振動(dòng)信號(hào)采集，并對(duì)模擬采出信號(hào)進(jìn)行分析，分為3種情況進(jìn)行對(duì)比，即無噪聲理想信號(hào)，信噪比為0.8的信號(hào)，信噪比為0.6的信號(hào)。所添加噪聲為高斯白噪聲，具體添加方式為：

其中，x表示理想仿真所得的傳感器信號(hào)數(shù)據(jù)；SNR表示信噪比，其取值范圍為SNR∈(0,1]，SNR越大，信號(hào)質(zhì)量越好；n為信號(hào)長度；f為添加高斯噪聲后的信號(hào)數(shù)據(jù)。

1號(hào)葉片在1號(hào)傳感器上的信噪比0.6的位移響應(yīng)信號(hào)及擬合情況如下圖5所示：

圖5 信噪比0.6的位移響應(yīng)信號(hào)及擬合曲線Fig.5 SNR=0.6 displacement response signal and fitting curve

以下以一號(hào)葉片為例進(jìn)行分析，速矢端跡擬合法擬合結(jié)果匯總見表2。

表2 不同信號(hào)質(zhì)量情況下的振動(dòng)中心頻率識(shí)別結(jié)果Tab.2 Vibration center frequency recognition results for different signal quality cases

如表2所示，在以上三種信號(hào)質(zhì)量下辨識(shí)的振動(dòng)中心頻率幾乎沒有差別，可選取76.2Hz（4 572r/min）為下一步恒速運(yùn)轉(zhuǎn)測(cè)試的運(yùn)轉(zhuǎn)頻率。

2）GARIV法恒速模擬仿真

選恒速運(yùn)轉(zhuǎn)頻率為76.2Hz，運(yùn)轉(zhuǎn)時(shí)間設(shè)置為30s，其他設(shè)置同前。1號(hào)傳感器采集的1號(hào)葉片信噪比0.6情況下的位移信號(hào)如圖6所示。

圖6 各傳感器信噪比0.6下恒速位移響應(yīng)信號(hào)Fig.6 Displacement response signal measured by each sensor（SNR=0.6）

GARIV法識(shí)別結(jié)果匯總見表3。

表3 不同信號(hào)質(zhì)量情況下的振動(dòng)倍頻識(shí)別情況表Tab.3 Identification of vibration frequency doubling indifferent signal quality cases

由表3可以看出在不同的信號(hào)質(zhì)量下，GARIV法識(shí)別的振動(dòng)倍頻值的精確度有一定的差異，但這個(gè)差異基本可以忽略不計(jì)，在76.2Hz情況下用GARIV法分析得到的倍頻值取整后均為30。

3）綜合評(píng)估識(shí)別結(jié)果

結(jié)合速矢端跡擬合法識(shí)別結(jié)果與GARIV法所得結(jié)果nEO=30，故可通過式（5）計(jì)算得到各個(gè)葉片的固有頻率ωn，結(jié)果匯總于表4，并與預(yù)設(shè)值進(jìn)行誤差計(jì)算。

表4 各個(gè)葉片固有頻率ωn識(shí)別結(jié)果（單位：Hz）Tab.4 The natural frequency recognition results of each blade（Unit:Hz）

由表4可知，該分析方法在不同的信號(hào)質(zhì)量下進(jìn)行分析得到的結(jié)果十分準(zhǔn)確，誤差均在萬分之7以內(nèi)。

3 葉片故障監(jiān)測(cè)實(shí)驗(yàn)研究

3.1 實(shí)驗(yàn)臺(tái)結(jié)構(gòu)與測(cè)試系統(tǒng)

為能通過實(shí)驗(yàn)驗(yàn)證前文所述葉片同步振動(dòng)辨識(shí)算法，構(gòu)建了高速直葉片振動(dòng)實(shí)驗(yàn)臺(tái)，具體細(xì)節(jié)和整體概貌如圖7所示。高速直葉片振動(dòng)監(jiān)測(cè)實(shí)驗(yàn)臺(tái)主要由底座、支撐體（電機(jī)）、操作平臺(tái)、測(cè)試輪盤、護(hù)罩、永磁激勵(lì)部件及其他輔助元部件組成。測(cè)試輪盤（具體結(jié)果見圖2）外安裝有金屬保護(hù)罩，護(hù)罩周向開有呈特定角度（最小6°）的若干個(gè)傳感器安裝孔，頂端設(shè)有鍵相傳感器（通過磁座固定在操作平臺(tái)上）安裝孔。護(hù)罩頂部開有12個(gè)永磁激勵(lì)部件安裝孔。每一個(gè)激勵(lì)部件都內(nèi)嵌一個(gè)鈷基永磁鐵，用以提供激振力，通過調(diào)整部件數(shù)量及分布情況，可開展若干不同激勵(lì)條件下的測(cè)試實(shí)驗(yàn)。

圖7 高速直葉片實(shí)驗(yàn)臺(tái)和測(cè)試系統(tǒng)圖Fig.7 High-speed straight blade test rig and testing system

測(cè)試系統(tǒng)中采用5通道激光源（1個(gè)通道用于監(jiān)測(cè)鍵相、其余通道用于葉尖計(jì)時(shí)），光功率可在0～1W之間連續(xù)調(diào)節(jié)。葉頂光纖傳感器的安裝要求為間隔夾角6°，為實(shí)現(xiàn)轉(zhuǎn)速控制，通過模擬輸出卡對(duì)變頻器進(jìn)行輸出控制，已達(dá)到按設(shè)定參數(shù)進(jìn)行升降速運(yùn)轉(zhuǎn)、恒速運(yùn)轉(zhuǎn)等。

3.2 實(shí)驗(yàn)測(cè)量

3.3.1 掃頻實(shí)驗(yàn)

在掃頻測(cè)量中，實(shí)驗(yàn)臺(tái)均布安裝6個(gè)磁鐵激勵(lì)，運(yùn)轉(zhuǎn)轉(zhuǎn)頻范圍為60～132Hz，通過4支葉頂傳感器監(jiān)測(cè)的葉片振動(dòng)位移隨轉(zhuǎn)頻的變化情況如圖8所示。注意不同傳感器對(duì)應(yīng)的葉片編號(hào)一定要修正為物理下的同一葉片。

圖8 4支傳感器監(jiān)測(cè)到的某一葉片振動(dòng)位移響應(yīng)情況Fig.8 Blade vibration displacement response by four sensors

由圖8可知葉片在該段運(yùn)轉(zhuǎn)區(qū)間內(nèi)被激起較大振幅的諧共振區(qū)間有4個(gè)，圖中用紅色圓圈示出。對(duì)各諧共振區(qū)間進(jìn)行截取共振區(qū)間片段擬合，下面將以1號(hào)葉頂傳感器采集的該葉片第2和第3個(gè)共振區(qū)間的擬合情況進(jìn)行說明。

如圖9截取范圍分別為73～78Hz，92～97Hz，用紅色曲線示出擬合曲線。

圖9 各段響應(yīng)曲線與擬合曲線之間的對(duì)比圖Fig.9 Comparison of the original and the fitting curve

通過LM法擬合獲得的目標(biāo)曲線的擬合參數(shù)，表5為4支傳感器擬合得到的各諧共振區(qū)內(nèi)的諧共振中心頻率的平均值。

表5 4支傳感器各段擬合參數(shù)結(jié)果平均值Tab.5 Average of the fitting parameters of each resonance interval by four sensors

3.3.2 恒速運(yùn)轉(zhuǎn)求取葉片振動(dòng)倍頻值

為了避免轉(zhuǎn)速波動(dòng)帶來的影響，以及運(yùn)行轉(zhuǎn)速下葉片被激起的振動(dòng)位移大小合適，實(shí)驗(yàn)中將恒速運(yùn)轉(zhuǎn)轉(zhuǎn)速設(shè)定在諧共振中心頻率附近，無需保證與中心頻率完全一致。

下面對(duì)上述的2個(gè)諧共振區(qū)分別恒速運(yùn)轉(zhuǎn)測(cè)量，圖10為運(yùn)轉(zhuǎn)轉(zhuǎn)頻分別為93.1Hz、74.35Hz時(shí)葉片振動(dòng)位移監(jiān)測(cè)情況。圖中對(duì)振動(dòng)位移數(shù)據(jù)進(jìn)行了7階線性平滑濾波處理。

圖10 不同運(yùn)轉(zhuǎn)頻率下的葉片振動(dòng)位移監(jiān)測(cè)圖Fig.10 Monitoring of blade vibration displacement under different rotating freguency

上述2個(gè)恒速運(yùn)轉(zhuǎn)測(cè)量采用GARIV法分析得到的葉片振動(dòng)倍頻值分別為24.45，30，取整后可得24，30。表6為各諧共振對(duì)應(yīng)的葉片振動(dòng)倍頻值和振動(dòng)頻率。

表6 各諧共振中心對(duì)應(yīng)的葉片動(dòng)頻值Tab.6 Blade dynamic frequency values corresponding to the harmonic resonance

3.3.3 整機(jī)分析

為研究整機(jī)葉片的振動(dòng)特性，實(shí)驗(yàn)葉輪上的32個(gè)葉片都可以按照上小節(jié)所述過程進(jìn)行分析。圖11為不同編號(hào)葉片的不同諧共振求得頻率變化曲線，從圖中分析可知整機(jī)32個(gè)葉片的一階振動(dòng)頻率不相同，為失諧葉片，這是由葉輪的制造誤差造成的。從單個(gè)葉片在不同倍頻下的共振角度分析，不同倍頻下求得的振動(dòng)頻率均相差不大，且隨著倍頻值的降低葉片振動(dòng)頻率值逐漸增大，這是由于隨著轉(zhuǎn)速的增大，轉(zhuǎn)子葉片的離心剛度增大，進(jìn)而葉片振動(dòng)頻率也有所增大。

圖11 各葉片在各倍頻下振動(dòng)頻率變化曲線圖Fig.11 The identified natural frequencies of the blades at each EO

4 結(jié)論

當(dāng)葉片上產(chǎn)生裂紋等故障時(shí)，其振動(dòng)頻率變化是一個(gè)明顯的特征。因此，通過識(shí)別葉片的固有頻率變化可對(duì)葉片裂紋故障進(jìn)行早期預(yù)警。提出了基于非線性最小二乘擬合和GARIV方法相結(jié)合的葉片同步振動(dòng)時(shí)的轉(zhuǎn)動(dòng)參數(shù)精確辨識(shí)方法，能夠準(zhǔn)確捕捉葉片的動(dòng)態(tài)固有頻率。仿真結(jié)果表明，本文所提出的方法可以達(dá)到萬分之7的測(cè)量精度，實(shí)驗(yàn)研究結(jié)果表明，本文所提出的方法對(duì)葉片固有頻率的識(shí)別是有效的。

[1]太興宇,楊樹華,馬輝,等.轉(zhuǎn)子-葉片系統(tǒng)固有特性分析及試驗(yàn)研究[J].風(fēng)機(jī)技術(shù),2017（2）:29-35.

[2]陶麗樺,王鍵,劉艷,等.葉片載荷分布對(duì)離心葉輪氣動(dòng)性能影響的數(shù)值研究[J].風(fēng)機(jī)技術(shù),2015（6）:13-20，92.

[3]董學(xué)智,張志鑫,許劍,等.串列風(fēng)機(jī)葉片間隙尺寸對(duì)風(fēng)機(jī)性能的影響研究[J].風(fēng)機(jī)技術(shù),2015（6）:44-48.

[4]賀長波,李宏坤,張曉雯.流體激振下葉片疲勞失效機(jī)理與檢測(cè)方法綜述[J].風(fēng)機(jī)技術(shù),2015（5）:69-77.

[5]李輝,王軍,周水清,等.葉片數(shù)對(duì)多翼離心風(fēng)機(jī)性能影響的分析[J].風(fēng)機(jī)技術(shù),2017（2）:19-22，47.

[6]薩勃洛斯基,吳士祥,鄭叔琛.渦輪機(jī)葉片振動(dòng)的非接觸測(cè)量[M].北京:國防工業(yè)出版社,1986.

[7]歐陽濤,段發(fā)階,閆明,等.雙參數(shù)法辨識(shí)葉片同步振動(dòng)的研究[J].傳感器與微系統(tǒng),2010,29（3）:42-45.

[8]Carrington,Ian B.Development of blade tip timing data analysis techniques[J].Journal of Dental Education,2002,47（2）:947-951.

[9]Gallegogarrido J,Dimitriadis G,Carrington I B,et al.A Class of Methods for the Analysis of Blade Tip Timing Data from Bladed Assemblies Undergoing Simultaneous Resonances—Part II: Experimental Validation[J].International Journal of Rotating Machinery,2007,2007（4）:981-1077.

[10]歐陽濤,郭文力,段發(fā)階,等.基于葉尖定時(shí)的旋轉(zhuǎn)葉片同步振動(dòng)辨識(shí)新方法[J].振動(dòng)與沖擊,2011,30（8）:249-252.

[11]歐陽濤.基于葉尖定時(shí)的旋轉(zhuǎn)葉片振動(dòng)檢測(cè)及參數(shù)辨識(shí)技術(shù)[D].天津：天津大學(xué)，2011.

Early Monitoring and Warning Method for Crack Faults in Blades

Ya ZhangKang ChenWei-ming WangXu-long Zhang
(Beijing University of Chemical Technology)

Vibration has always been a major obstacle in the development of enlarged and high-parametric rotating equipment.Blade vibration is the main reason of fatigue failure and rotating machine accidents.In this paper,the synchronous resonance of blades is studied with the combination of the nonlinear least squares fitting and the GARIV method to accurately capture the dynamic eigen frequency of the blade.An early warning for blade crack failures is realized by using the blade tip timing(BTT)method of laser sensors,which captures the small changes of the eigen frequency due to an appearing blade crack.The simulation results show the measurement accuracy of the proposed method can achieve 0.7‰.Thus,it can be used in the early identification of blade cracks,which is confirmed by experimental results.

blade vibration,parameter identification,simulation,experimental test,blade tip timing,eigen frequency

TB53；TK474

1006-8155-(2017)04-0045-07

A

10.16492/j.fjjs.2017.04.0008

國家自然科學(xué)基金(51775030)

2017-05-15 北京 100029