五軸機(jī)床直驅(qū)回轉(zhuǎn)工作臺(tái)熱特性仿真研究*

趙 繼 馬 強(qiáng) 王紫瑄 于天彪

(東北大學(xué)機(jī)械工程與自動(dòng)化學(xué)院，遼寧 沈陽(yáng) 110819)

?

五軸機(jī)床直驅(qū)回轉(zhuǎn)工作臺(tái)熱特性仿真研究*

趙 繼 馬 強(qiáng) 王紫瑄 于天彪

(東北大學(xué)機(jī)械工程與自動(dòng)化學(xué)院，遼寧 沈陽(yáng) 110819)

針對(duì)直驅(qū)回轉(zhuǎn)工作臺(tái)熱變形導(dǎo)致五軸機(jī)床產(chǎn)生加工誤差的問(wèn)題，以耳軸式回轉(zhuǎn)工作臺(tái)虛擬樣機(jī)為研究對(duì)象進(jìn)行仿真分析。綜合利用Fluent模塊與ANSYS Workbench建立工作臺(tái)的流-熱-固單向耦合仿真模型，對(duì)工作臺(tái)的熱特性進(jìn)行研究。在給定工況下，得到工作臺(tái)達(dá)到熱平衡狀態(tài)時(shí)的溫度分布與熱變形情況。通過(guò)簡(jiǎn)化模型計(jì)算冷卻液對(duì)流換熱系數(shù)的理論估算值，并與仿真結(jié)果比較，驗(yàn)證仿真分析有較高的可信度。研究結(jié)果可為工作臺(tái)的后續(xù)結(jié)構(gòu)改進(jìn)設(shè)計(jì)與冷卻液的使用提供參考。

回轉(zhuǎn)工作臺(tái)；熱特性；耦合仿真；對(duì)流換熱系數(shù)

高動(dòng)態(tài)性能的五軸機(jī)床回轉(zhuǎn)工作臺(tái)普遍采用力矩電動(dòng)機(jī)直接驅(qū)動(dòng)。力矩電動(dòng)機(jī)產(chǎn)熱量較高，且封裝在工作臺(tái)內(nèi)部，一般需要通過(guò)冷卻液散熱。在精密加工中，機(jī)床熱變形引起的加工誤差可達(dá)總誤差的75%[1]。可通過(guò)物理(如改變材料、結(jié)構(gòu)、散熱等)和軟件補(bǔ)償兩種方式來(lái)減小熱變形誤差[2]。但其前提是了解工作臺(tái)熱特性，找到其熱變形的規(guī)律。國(guó)內(nèi)外學(xué)者從多個(gè)方面對(duì)力矩電動(dòng)機(jī)及工作臺(tái)的熱特性進(jìn)行研究，如建立電動(dòng)機(jī)三維暫態(tài)溫度場(chǎng)計(jì)算模型[3]、實(shí)機(jī)測(cè)量溫升與熱變形[4]、研究切削液對(duì)機(jī)床熱變形的影響[5]等。在高端數(shù)控機(jī)床的設(shè)計(jì)階段即應(yīng)重視其熱特性[6]。因此，本文以耳軸式回轉(zhuǎn)工作臺(tái)虛擬樣機(jī)為研究對(duì)象，利用Fluent與ANSYS Workbench進(jìn)行流—熱—固耦合仿真，分析其溫度分布與熱變形情況，并驗(yàn)證了仿真結(jié)果的可信度。研究結(jié)果可為工作臺(tái)結(jié)構(gòu)改進(jìn)及冷卻液的使用提供重要參考。

1 工作臺(tái)有限元分析模型建立

1.1 工作臺(tái)結(jié)構(gòu)性能參數(shù)

本文所研究的工作臺(tái)采用B/C軸設(shè)計(jì)，YRT轉(zhuǎn)臺(tái)軸承連接回轉(zhuǎn)部件，直驅(qū)力矩電動(dòng)機(jī)提供回轉(zhuǎn)動(dòng)力。工作臺(tái)的三維模型如圖1所示，其部分結(jié)構(gòu)性能參數(shù)如表1所示。

1.2 仿真分析模型建立

為減小計(jì)算規(guī)模，提高效率，在導(dǎo)入ANSYS 前，應(yīng)對(duì)模型進(jìn)行適當(dāng)簡(jiǎn)化，省略對(duì)分析影響很小的孔、倒角和圓角等小特征[7]。簡(jiǎn)化后工作臺(tái)模型為對(duì)稱結(jié)構(gòu)，在仿真研究中以對(duì)稱模式進(jìn)行分析可大幅減小計(jì)算規(guī)模，縮短計(jì)算時(shí)間[8]，并可避免因分網(wǎng)不均造成的仿真誤差，提高計(jì)算精度。導(dǎo)入ANSYS的模型如圖2所示。

表1 工作臺(tái)性能參數(shù)表

主要參數(shù)數(shù)值回轉(zhuǎn)工作臺(tái)直徑/mm600工作臺(tái)最大承重/kg160工件最大回轉(zhuǎn)直徑/mm480擺動(dòng)定位精度/重復(fù)精度/(″)25/10回轉(zhuǎn)定位精度/重復(fù)精度/(″)15/10B軸轉(zhuǎn)角/(°)-90～+90C軸轉(zhuǎn)角/(°)0/360

仿真分析模型需要正確的材料參數(shù)。工作臺(tái)的主體材料為HT200，兩軸承的材料為軸承鋼GCr15SiMn，電動(dòng)機(jī)繞組為銅，牌號(hào)T1，電動(dòng)機(jī)定轉(zhuǎn)子結(jié)構(gòu)材料為35SiMn。常溫下材料的參數(shù)如表2所示。

表2 結(jié)構(gòu)材料參數(shù)表

材料參數(shù)HT200GCr15SiMnT135SiMn密度/(t/m3)7.37.828.937.8比熱容/(J/(kg·K))510443386436導(dǎo)熱系數(shù)/(W/(m·K))523739848線膨脹系數(shù)/(×10-6/K)10.111.217.110.2彈性模量/GPa122216137214泊松比0.30.30.310.27

2 熱變形分析

2.1 熱邊界條件確定

工作臺(tái)的熱分布由其產(chǎn)熱與散熱情況決定。熱源分為內(nèi)部熱源與外部熱源。工作臺(tái)的內(nèi)部熱源主要為電動(dòng)機(jī)與軸承，外部熱源包括工件切削熱，其他設(shè)備傳熱，人、陽(yáng)光等，其影響相對(duì)較小，但更為復(fù)雜。本文著重對(duì)內(nèi)部熱源進(jìn)行研究，忽略外部熱源的影響。工作臺(tái)產(chǎn)生的熱量主要通過(guò)3種方式散失：冷卻液、空氣與基座。力矩電動(dòng)機(jī)產(chǎn)生大量熱量，其中絕大部分被冷卻液帶走，故冷卻液對(duì)流換熱部分是研究的重點(diǎn)。

力矩電動(dòng)機(jī)的主要損耗為銅耗與鐵耗，其中銅耗占總損耗的絕大部分[9]。假定電動(dòng)機(jī)損耗全部轉(zhuǎn)換為熱量，忽略風(fēng)阻損耗等附加損耗的影響[10]。在線性工作區(qū)，電動(dòng)機(jī)的銅耗PCu可通過(guò)電動(dòng)機(jī)常數(shù)和轉(zhuǎn)矩求得，其計(jì)算公式如下：

(1)

軸承的產(chǎn)熱源于摩擦，若摩擦損耗全部轉(zhuǎn)化為熱量，則其生熱功率即為摩擦力矩與軸承旋轉(zhuǎn)角速度的乘積[11]。總摩擦力矩與轉(zhuǎn)速、潤(rùn)滑劑參數(shù)、載荷、滾動(dòng)體與滾道彈性變形量等因素有關(guān)。根據(jù)軸承廠商的產(chǎn)品手冊(cè)，YRT系列軸承摩擦力矩與轉(zhuǎn)速近似呈線性關(guān)系。由該關(guān)系曲線可估算軸承的生熱功率。

電動(dòng)機(jī)與軸承的生熱功率由工作臺(tái)設(shè)計(jì)要求的典型工況確定。經(jīng)計(jì)算，B軸電動(dòng)機(jī)的生熱功率為2 153.8 W，軸承生熱功率691 W；C軸電動(dòng)機(jī)的生熱功率為2 421.5 W，軸承生熱功率為282.7 W。

由于機(jī)床工作臺(tái)溫升較小，通過(guò)熱輻射散失的熱量忽略不計(jì)[12]。空氣對(duì)流換熱情況復(fù)雜：工作臺(tái)同時(shí)有運(yùn)動(dòng)與靜止部件，運(yùn)動(dòng)部件攪動(dòng)空氣，又會(huì)使空氣相對(duì)靜止部件運(yùn)動(dòng)。各部件的形狀尺寸及運(yùn)動(dòng)部件的運(yùn)動(dòng)速度、時(shí)間等因素都會(huì)影響空氣換熱。另一方面，臺(tái)面上的工件、切屑及切削液也會(huì)影響工作臺(tái)散熱，切削熱甚至?xí)魅牍ぷ髋_(tái)，成為熱源。因此，求解工作臺(tái)表面空氣對(duì)流換熱系數(shù)的精確值是沒(méi)有意義的。為簡(jiǎn)化仿真分析，這里假設(shè)工作臺(tái)臺(tái)面無(wú)工件，且其轉(zhuǎn)動(dòng)較慢，表面空氣按照自然對(duì)流換熱取為10 W/(m2·K)。

散熱邊界條件計(jì)算的重點(diǎn)是冷卻液對(duì)流換熱系數(shù)的確定。冷卻液在冷卻套內(nèi)的管路如圖3所示。液體從入口處流入，分為兩路繞到電動(dòng)機(jī)另一端匯合，進(jìn)入下一層散熱管路并再次分開(kāi)繞到電動(dòng)機(jī)另一端，以此往復(fù)，直至最終在出口管路處匯合流出。

根據(jù)電動(dòng)機(jī)廠商的要求，冷卻液的主要成分為添加防腐劑的水。合適的混合冷卻液為25%～30% 的乙二醇水溶液。冷卻液的參數(shù)隨溫度變化較為明顯。本文設(shè)定冷卻水入口溫度維持在20 ℃，間隔1 ℃獲得冷卻液從20 ℃到27 ℃的參數(shù)。根據(jù)生產(chǎn)商的數(shù)據(jù)，25%乙二醇水溶液的參數(shù)如表3所示。

表3 冷卻液參數(shù)表

溫度/℃密度/(kg/m3)比熱容/(kJ/(kg·K))導(dǎo)熱系數(shù)/(W/(m·K))運(yùn)動(dòng)粘度/(mm2/s)動(dòng)力粘度/Pa·s2010353.830.4921.950.0020182110343.830.4931.90.001962210343.830.4931.840.0019022310343.830.4941.790.0018502410333.830.4941.750.0018072510333.830.4951.70.0017562610323.830.4961.660.0017132710323.830.4961.610.001661

對(duì)流換熱系數(shù)為單位溫差作用下通過(guò)單位面積的熱流量。其物理意義根據(jù)如下牛頓冷卻公式得到：

Φ=AhΔt

(2)

式中：Φ為熱流量，W；A為面積，m2；Δt為流體和物體表面的溫差，約定永遠(yuǎn)取正值，℃；H為對(duì)流換熱系數(shù)，W/(m2·K)。

傳熱過(guò)程中的很多因素都會(huì)影響對(duì)流換熱系數(shù)的大小，如流體的物性(密度、粘度、導(dǎo)熱系數(shù)、比熱容等)、流動(dòng)的形態(tài)(層流、紊流)、流動(dòng)的成因(自然對(duì)流與強(qiáng)迫對(duì)流)、物體表面的形狀及尺寸。因此，管路中冷卻液各部分的對(duì)流換熱系數(shù)不同，即各部分的冷卻能力不同。本文通過(guò)Fluent計(jì)算冷卻液的對(duì)流換熱系數(shù)分布，并作為邊界條件導(dǎo)入后續(xù)的熱分析中。

工作臺(tái)共使用兩臺(tái)力矩電動(dòng)機(jī)，兩電動(dòng)機(jī)的冷卻原理、冷卻液都相同，只有型號(hào)不同，對(duì)應(yīng)冷卻套尺寸及冷卻液流量不同。下文只分析C軸電動(dòng)機(jī)，B軸電動(dòng)機(jī)計(jì)算方法相似，不再重復(fù)說(shuō)明。根據(jù)電機(jī)廠商產(chǎn)品手冊(cè)，C軸電動(dòng)機(jī)的冷卻液參考流量為15.59 L/min，參考溫升為5 ℃。為減小計(jì)算規(guī)模，此處假設(shè)電動(dòng)機(jī)產(chǎn)生的熱量都由冷卻液帶走。對(duì)流換熱系數(shù)的分布云圖如圖4所示。從入口處開(kāi)始沿直徑方向，各節(jié)點(diǎn)的對(duì)流換熱系數(shù)如圖5所示。

由圖4可知，冷卻液由單路分散為兩路的位置，對(duì)流換熱系數(shù)遠(yuǎn)大于其他部分。具體表現(xiàn)為入口附近，冷卻液沖擊管壁后分為兩路，流動(dòng)狀況復(fù)雜，流速快，對(duì)流換熱系數(shù)非常大；每一層冷卻液匯聚成單路進(jìn)入下一層并再次分開(kāi)的位置系數(shù)較大，但變化較為平緩。圖5描述了各節(jié)點(diǎn)對(duì)應(yīng)冷卻液的換熱能力，其橫軸代表入口方向沿直徑的距離，縱軸代表節(jié)點(diǎn)處的對(duì)流換熱系數(shù)。分為兩路后的冷卻液，其對(duì)流換熱系數(shù)比較穩(wěn)定，在2 000～2 380 W/(m2·K)。將Fluent的仿真結(jié)果作為載荷導(dǎo)入穩(wěn)態(tài)熱分析模塊可求解工作臺(tái)溫度分布。

2.2 工作臺(tái)穩(wěn)態(tài)熱分析

本文研究工作臺(tái)處于熱平衡狀態(tài)的熱變形。確定熱邊界條件后，用ANSYS的穩(wěn)態(tài)熱分析模塊進(jìn)行溫度分布的求解。熱穩(wěn)態(tài)即流入系統(tǒng)的熱量加上系統(tǒng)自身生成的熱量與流出系統(tǒng)的熱量相等，此時(shí)任意節(jié)點(diǎn)的溫度不隨時(shí)間變化。

假定初始溫度為20 ℃，環(huán)境溫度持續(xù)穩(wěn)定在20 ℃，則處于熱平衡狀態(tài)的工作臺(tái)溫度分布如圖6所示。如圖可知，工作臺(tái)最高溫度為30.672 ℃，位于C軸電動(dòng)機(jī)定子，最低溫度位于冷卻液入口附近，表明冷卻液帶走大量的熱。

2.3 工作臺(tái)穩(wěn)態(tài)結(jié)構(gòu)分析

將溫度分布作為載荷導(dǎo)入結(jié)構(gòu)靜態(tài)分析模塊，對(duì)工作臺(tái)支撐的外圍施加固定約束。工作臺(tái)穩(wěn)態(tài)熱變形如圖7所示，其最大變形量達(dá)到78 μm。參照?qǐng)D7仿真模型的坐標(biāo)系，臺(tái)面中心點(diǎn)在X、Y、Z三個(gè)方向的位移分別為0、-2.3 μm、54.1 μm。X方向位移為0，這是由于分析采用對(duì)稱模型，臺(tái)面中心點(diǎn)位于對(duì)稱面Y-Z平面，其X向位移永遠(yuǎn)為0。臺(tái)面所處位置靠近B軸軸線，即靠近擺動(dòng)部件的回轉(zhuǎn)中心，因此受熱時(shí)Y方向位移很小，這體現(xiàn)出工作臺(tái)結(jié)構(gòu)對(duì)稱設(shè)計(jì)的優(yōu)勢(shì)。而Z方向位移超過(guò)了50 μm，這需要后續(xù)改進(jìn)結(jié)構(gòu)，并通過(guò)熱誤差補(bǔ)償模塊進(jìn)行軟件補(bǔ)償。

2.4 仿真結(jié)果可信性驗(yàn)證

工作臺(tái)實(shí)際的溫度分布及熱變形難以進(jìn)行理論計(jì)算，然而冷卻液的對(duì)流換熱系數(shù)可以通過(guò)理論計(jì)算給出估算值。本節(jié)對(duì)Fluent仿真結(jié)果進(jìn)行驗(yàn)證，以檢驗(yàn)仿真研究的可信性。

在衛(wèi)浴空間加入天窗設(shè)計(jì)，可使原本狹小的空間一下子就變得空曠起來(lái)，圓弧形的天窗順應(yīng)房頂結(jié)構(gòu)，仿佛一塊LED屏幕，在保證私密性之余，不僅為沐浴間提供了室外的別樣美景，也最大化地提供了光照、通風(fēng)，使浴室獲得別樣的明亮溫暖。

為簡(jiǎn)化計(jì)算，研究單管路內(nèi)遠(yuǎn)離入口部分的冷卻液，并假設(shè)此處冷卻液的流動(dòng)不受管路兩端流體沖擊與匯合的影響，冷卻液溫度為20 ℃。此時(shí)問(wèn)題簡(jiǎn)化為一種簡(jiǎn)單的管內(nèi)強(qiáng)制對(duì)流換熱。

不同流動(dòng)狀態(tài)下，流體的對(duì)流傳熱規(guī)律不同。管內(nèi)強(qiáng)制對(duì)流，流體的流動(dòng)狀態(tài)分為層流和湍流兩種，其可通過(guò)雷諾數(shù)Re判別，雷諾數(shù)的計(jì)算公式如下：

(3)

式中：u為速度，m/s；d為直徑，m；υ為流體的運(yùn)動(dòng)粘度，m2/s。

一般認(rèn)為當(dāng)Re<2 300時(shí)，流體的流動(dòng)狀態(tài)為層流；當(dāng)Re>10 000時(shí)，流體的流動(dòng)狀態(tài)為旺盛湍流；當(dāng)2 300≤Re≤10 000時(shí)，為層、湍流的過(guò)渡狀態(tài)，流體在這3種流動(dòng)狀態(tài)下的對(duì)流傳熱規(guī)律不同[13]。

經(jīng)計(jì)算，雷諾數(shù)Re=5 131，因此冷卻液的流動(dòng)狀態(tài)為層、湍流過(guò)渡狀態(tài)。

(4)

式中：cp為定壓比熱容，J/(kg·K)；μ為流體的動(dòng)力粘度，Pa·s；λ為導(dǎo)熱系數(shù)，W/(m·s)。

定義努塞爾數(shù)Nu，對(duì)于過(guò)渡狀態(tài)，根據(jù)豪森公式確定對(duì)流換熱準(zhǔn)則，Nu的計(jì)算公式如下：

(5)

式中：de為當(dāng)量直徑，m；l為管道總長(zhǎng)，m；μw為以壁面溫度為定性溫度的流體動(dòng)力粘度，Pa·s。

對(duì)流換熱系數(shù)的計(jì)算公式如下：

(6)

經(jīng)計(jì)算，冷卻液對(duì)流換熱系數(shù)h的理論估算值為2 240 W/(m2·K) ，該值在Fluent仿真結(jié)果穩(wěn)定段的變化區(qū)間內(nèi)。因此，仿真結(jié)果的可信度較高。

3 結(jié)語(yǔ)

應(yīng)用Fluent與ANSYS Workbench對(duì)五軸機(jī)床耳軸式回轉(zhuǎn)工作臺(tái)虛擬樣機(jī)進(jìn)行流—熱—固單向耦合分析，得出如下結(jié)論：

(1) 冷卻液的冷卻能力與其流動(dòng)狀態(tài)密切相關(guān)。平穩(wěn)流動(dòng)時(shí)，管路內(nèi)液體的對(duì)流換熱系數(shù)穩(wěn)定在2 000～2 380 W/(m2·K)，在劇烈運(yùn)動(dòng)的位置，對(duì)流換熱系數(shù)有數(shù)倍的增大，這將增強(qiáng)冷卻液的整體換熱能力。

(2) 在給定的加工要求下運(yùn)行至熱平衡狀態(tài)，工作臺(tái)最大溫升為10 ℃，位于C軸電動(dòng)機(jī)定子附近，最低溫度位于冷卻液入口管路附近，冷卻液的冷卻效果顯著。

(3) 熱平衡時(shí)，工作臺(tái)的最大熱變形達(dá)到78 μm，臺(tái)面中心點(diǎn)沿機(jī)床B軸的位移量達(dá)到54 μm，大于精度的設(shè)計(jì)標(biāo)準(zhǔn)，在工作臺(tái)的后續(xù)設(shè)計(jì)中需要通過(guò)改進(jìn)結(jié)構(gòu)或采用軟件補(bǔ)償?shù)姆绞綔p小熱變形造成的加工誤差。

[1]Mayr J, Jedrzejewski J, Uhlmann E, et al. Thermal issues in machine tools[J]. CIRP Annals - Manufacturing Technology, 2012, 61(2): 771-791.

[2]Gebhardt Michael. Thermal behaviour and compensation of rotary axes in 5-axis machine tools[D]. ETH-Zürich, 2014.

[3]鄒繼斌,張洪亮,江善林,等. 電磁穩(wěn)態(tài)條件下的力矩電機(jī)三維暫態(tài)溫度場(chǎng)分析[J]. 中國(guó)電機(jī)工程學(xué)報(bào), 2007(21): 66-70.

[4]王秀山,楊建國(guó),余永昌,等. 雙轉(zhuǎn)臺(tái)五軸數(shù)控機(jī)床熱誤差建模、檢測(cè)及補(bǔ)償實(shí)驗(yàn)研究[J]. 中國(guó)機(jī)械工程, 2009,20(4): 405-408.

[5]Mayr J, Gebhardt M, Massow B B, et al. Cutting fluid influence on thermal behavior of 5-axis machine tools[J]. 6th CIRP International Conference on High Performance Cutting, HPC2014, 2014, 14: 395-400.

[6]張曙. 機(jī)床產(chǎn)品創(chuàng)新與設(shè)計(jì)[M]. 南京: 東南大學(xué)出版社, 2014.

[7]于天彪,王學(xué)智,關(guān)鵬,等. 超高速磨削機(jī)床主軸系統(tǒng)模態(tài)分析[J]. 機(jī)械工程學(xué)報(bào), 2012，48(17): 183-188.

[8]Mian N S, Fletcher S, Longstaff A P, et al. Efficient thermal error prediction in a machine tool using finite element analysis[J]. Measurement Science and Technology, 2011, 22(8): 085107.

[9]張亮,王艷萍,蔣正華. 機(jī)殼水冷無(wú)刷力矩電動(dòng)機(jī)溫升估算研究[J]. 微電機(jī), 2013(12): 27-29，38.

[10]孫振宏. 基于流—固—熱耦合的電主軸熱特性仿真分析與實(shí)驗(yàn)研究[D]. 天津: 天津大學(xué), 2014.

[11]鄧小雷. 數(shù)控機(jī)床主軸系統(tǒng)多物理場(chǎng)耦合熱態(tài)特性分析研究[D]. 杭州：浙江大學(xué), 2014.

[12]張耀滿,高冠濱,王旭東,等. 加工中心主軸部件及其主軸箱的熱特性有限元分析[J]. 組合機(jī)床與自動(dòng)化加工技術(shù), 2005(4): 43-44，52.

[13]戰(zhàn)洪仁. 工程傳熱學(xué)基礎(chǔ)[M]. 北京: 中國(guó)石化出版社, 2014.

如果您想發(fā)表對(duì)本文的看法，請(qǐng)將文章編號(hào)填入讀者意見(jiàn)調(diào)查表中的相應(yīng)位置。

Simulation study on thermal characteristic of direct drive rotary table in 5-axis machining tools

ZHAO Ji, MA Qiang, WANG Zixuan, YU Tianbiao

(School of Mechanical Engineering and Automation, Northeastern University, Shenyang 110819, CHN)

Thermal deformation of direct drive rotary table causes error of five axis machining tools. To solve the problem, a simulation analysis model is established based on a trunnion rotary table virtual prototype. With the help of Fluent and ANSYS Workbench, the fluid-thermal-solid unidirectional coupling simulation model can be used to study the thermal characteristic of table. In the given conditions, the temperature distribution and thermal deformation of the table are obtained, when the table reaches thermal equilibrium state. Simulation results are compared with the theoretical estimates of convective heat transfer coefficient through the simplified model, and it is verified that the simulation has high credibility. The results can provide a reference for subsequent works such as optimal design of structure and use of cooling liquid.

rotary table; thermal characteristic; coupling simulation; convective heat transfer coefficient

* 遼寧省科技攻關(guān)項(xiàng)目( 2012220031)

TG502.15

A

10.19287/j.cnki.1005-2402.2016.11.026

趙繼，男，1959年生，博士，教授，研究方向?yàn)橹悄芫苤圃臁?/p>

(編輯 汪 藝)

2016-04-19)

161133