簡易檢索 / 詳目顯示

研究生: 陳明揚
Chen, Ming-Yang
論文名稱: 多顆小型垂直軸風機流場之數值研究
Numerical Study of Flow Around Multiple Small Veritcal-Axis Wind Turbines
指導教授: 黃啟鐘
Hwang, Chii-Jong
學位類別: 碩士
Master
系所名稱: 工學院 - 航空太空工程學系
Department of Aeronautics & Astronautics
論文出版年: 2012
畢業學年度: 100
語文別: 中文
論文頁數: 109
中文關鍵詞: 垂直軸風機多顆垂直風機分離渦流模擬噪音
外文關鍵詞: VAWT, Multiple VAWT, DES, Aerodynamic Efficiency, Noise
相關次數: 點閱:106下載:5
分享至:
查詢本校圖書館目錄 查詢臺灣博碩士論文知識加值系統 勘誤回報
  • 近年來為解決能源危機問題,各種替代能源研究正蓬勃發展。現今全球利用風機發電明顯增加,本計算使用商業套裝軟體FLUENT探討葉片式垂直軸風機(Straight-Bladed Vertical Axis Wind Turbine, VAWT)流場及其相關之效率與噪音現象。利用SIMPLEC/QUICK之數值方法、k-ω紊流模型與分離渦流模擬法DES(Detached Eddy Simulation)在四邊形/三角形網格上求解非穩態不可壓縮納維史托克方程式(Navier-Stokes Equation),接著採用聲場模型FW-H(Ffowcs Williams and Hawkings Equation)方程式計算流場中之噪音分貝值。首先為了解紊流模型與DES特性,進行葉片翼型為NACA0018之風機流場計算,並與相關之數值與實驗值比較以評估本數值模擬。接著進行葉片翼型為NACA0015之單顆風機流場模擬,以了解網格建構、疏密及軟體內部參數設定之影響。最後,針對葉片翼型為NACA0015之雙顆風機流場,在不同位置以及旋轉方向下進行探討,並與相關之實驗值做比較。在噪音研究方面,首先進行機翼葉片為NACA0015之單顆風機聲場計算以探討自由流通過單顆風機時在各個位置造成之聲音分貝值。接著進行葉片翼型為NACA0015之雙顆風機流場與聲場模擬,以了解兩顆風機之距離所造成之聲音分貝值。總之,藉由本文之模擬解果可得知雙顆風機之流場與噪音之現象。

    In the recent years, variety of alternative energy is booming to solve the problem of energy crisis. The electric generation using wind power is increasing around the globe curently and we used software package (FLUENT) to investigate the power efficiency and noise of the flow field around the straight-bladed vertical axis wind turbine in this study. In the present computations, the SIMPLEC / QUICK method,k-ω turbulence model and Detached Eddy Simulation model are used to solve unsteady incompressible Navier-Stokes Equation on the quadrilateral / triangular meshes, then the Ffowcs Willams and Hawkings equation is introduced to calculate the decibel value of acoustic field. In order to understand the characteristics of turbulent model and DES, the flow field around wind turbine with NACA0018 blade is computed first, and the results are compared with numerical and experimental values in the related literatures. Then, the flow field around wind turbine with NACA0015 blade is simulated to understand the effects due to mesh construction, mesh intensity and parameter setting. Finally, the flows around wind turbine arrays with different location and direction of rotation between two NACA0015 wind turbine blades are studied, and the numerical result are compared with experimental data. For the noise study, the acoustic field around a single wind turbine with NACA0015 blade is computed to understand the relation between dB value and position of wind turbine.Finally, the flow field and acoustic field simulation around two wind turbines are processed to calculate the dB values for the different distance between two turbines. In the conclusion, the flow fields and noise phenomenon between two wind turbines can be investigated .

    目錄 摘要……………………………………………………………………………I Abstract……………………………………………………………………...III 誌謝…………………………………………………………………………..V 目錄………………………………………………………………………….VI 主要符號說明…………………………………………………………...XVIII 第一章 緒 論…………………………………………………………….1 1-1 前言………….…………………………………………………….1 1-2 動機與目的….…………………………………………………….1 1-3 文獻回顧…………………………………………………………..3 1-4 基礎理論…………………………………………………………..6 1-4-1 Betz極限………………………………………………….6 1-4-2 座標系統………………………………………………….7 1-4-3 自行啟動能力(Self Starting)……………………………..9 1-5 研究內容…………………………………………………………..9 第二章 數 值 方 法 與 數 學 模 式….……………………………10 2-1 數值方法…………………………………………………………10 2-1-1 統御方程式.…………………………………………….10 2-1-2 SIMPLEC演算法……………………………………….12 2-1-3 QUICK法……………………………………….………14 2-2 DES模型…………………………………………………………15 2-3 噪音模型…………………………………………………………16 2-3-1 Ffowcs Williams and Hawkings方程式………………...16 2-3-2 聲壓值……………………………………….………….18 第三章 前 處 理 器 與 求 解 器 之 設 置……………………….19 3-1 NACA0018 之風機流場模擬…………………………………..19 3-1-1 NACA0018 之幾何外型與網格建立………………….19 3-1-2 NACA0018 之流場、邊界條件與參數設定………….20 3-2 NACA0015 之單顆風機流場模擬……………………………..22 3-2-1 NACA0015 之幾何外型與網格建立………………….22 3-2-2 NACA0015 之流場、邊界條件與參數設定………….22 3-2-3 NACA0015 之幾何外型與疏密網格建立…………….24 3-2-4 NACA0015 之流場、邊界條件與參數設定………….25 3-3 NACA0015 之雙顆風機流場模擬……………………………..26 3-3-1 NACA0015 之幾何外型與網格建立………………….26 3-3-2 NACA0015 之流場、邊界條件與參數設定………….27 第四章 結 果 與 討 論……………………………………………….29 4-1 NACA0018之風機流場模擬探討………………………………29 4-2 NACA0015之單顆風機流場與聲場模擬探討…………………30 4-3 NACA0015之雙顆風機流場與聲場模擬探討…………………32 第五章 結 論 與 建 議……………………………………………….34 5-1 結論………………………………………………………………34 5-2 建議………………………………………………………………36 參考文獻…………………………………………………………………….38 表目錄 表4-1 不同滑移半徑下之效率係數比較(λ=2)………………………….89 表4-2 觀測點位置表……………………………………………………...90 表4-3 兩種網格滑移區域半徑為1.5m和2m上觀測點之噪音值……..90 表4-4 疏和密網格上不同觀測點上之噪音值(滑移區域半徑為3.75m) .91 表4-5 雙顆風機在不同配置下之效率係數…………………………….103 表4-6 觀測點之座標與噪音值(x/D,y/D)=(1.5,-1.5)……………………103 表4-7 觀測點之座標與噪音值(x/D,y/D)=(2.0,1.5)…………..…………104 表4-8 雙顆風機在不同配置下之效率係數(密網格)………………...…105 圖目錄 圖1-1 全球總裝機容量MW[27]…………………………………………42 圖1-2 經濟部能源局預期台灣總裝機容量[28]………………………….42 圖1-3 Savonius型 [7]…………………………………………………….43 圖1-4 Darrieus型 [7]……………………………………………………..43 圖1-5 直立葉片垂直軸風機(Straight-Bladed VAWT)[7]…………….44 圖1-6 動量理論之示意圖 [22]…………………………………………...44 圖1-7 風機葉片座標定義 [22]…………………………………………...45 圖1-8 相對速度向量圖 [22]……………………………………………...45 圖1-9 不同尖端速度比之攻角隨旋轉角度位置的變化 [22]…………...46 圖1-10 典型風機效率值隨尖端速度比變化圖 [10]…………………….46 圖1-11 相對速度與攻角關係圖 [10]…………………………………….47 圖1-12 效率值隨尖端速度比變化圖 [4]………………………………...47 圖1-13 各角度下NACA0015之風力機所提供之Reaction torque [22]..48 圖1-14 m/s,尖端速度比與扭力係數關係 [22]………………..48 圖1-15 m/s,尖端速度比與扭力係數關係 [22]………………..49 圖2-1 交錯格點系統 [23]………………………………………………...49 圖2-2 SIMPLEC計算流程圖 [23]……………………………………….50 圖2-3 QUICK-Scheme一維控制體積系統 [23]………………………...50 圖3-1 NACA0018 流場密網格…………………………………………..51 圖3-2 NACA0018 內部滑移區域密網格………………………………..51 圖3-3 NACA0018環繞葉片區域密網格………………………………...52 圖3-4 NACA0018 密網格偏斜度分布…………………………………..52 圖3-5 NACA0018 流場疏網格…………………………………………..53 圖3-6 NACA0018 內部滑移區域疏網格………………………………..53 圖3-7 NACA0018 環繞葉片區域疏網格………………………………..54 圖3-8 NACA0018 疏網格偏斜率分布…………………………………..54 圖3-9 NACA0015 流場網格內部滑移半徑為1.5m…………………….55 圖3-10 NACA0015 內部滑移區域網格半徑為1.5m…………………...55 圖3-11 NACA0015環繞葉片區域網格………………………………….56 圖3-12 NACA0015網格偏斜分布……………………………………….56 圖3-13 NACA0015 流場網格內部滑移半徑為2m……………………..57 圖3-14 NACA0015 內部滑移區域網格半徑為2m……………………..57 圖3-15 NACA0015環繞葉片區域網格………………………………….58 圖3-16 NACA0015網格偏斜分布……………………………………….58 圖3-17 NACA0015 流場疏網格…………………………………………59 圖3-18 NACA0015 內部滑移區域疏網格………………………………59 圖3-19 NACA0015環繞葉片區域疏網格……………………………….60 圖3-20 NACA0015疏網格偏斜分布…………………………………….60 圖3-21 NACA0015 流場密網格…………………………………………61 圖3-22 NACA0015 內部滑移區域密網格………………………………61 圖3-23 NACA0015環繞葉片區域密網格……………………………….62 圖3-24 NACA0015密網格偏斜分布…………………………………….62 圖3-25 NACA0015 (x/D,y/D)=(1.5,-1.5)流場網格………………………63 圖3-26 NACA0015 (x/D,y/D)=(1.5,-1.5)內部滑移區域網格……………63 圖3-27 NACA0015環繞葉片區域網格………………………………….64 圖3-28 NACA0015整體網格偏斜分布………………………………….64 圖3-29 NACA0015 (x/D,y/D)=(2.0,1.5)流場網格……………………….65 圖3-30 NACA0015 (x/D,y/D)=(2.0,1.5)內部滑移區域網格…………….65 圖3-31 NACA0015環繞葉片區域網格………………………………….66 圖3-32 NACA0015整體網格偏斜分布………………………………….66 圖4-1 文獻實驗與模擬之升力係數隨攻角變化[10]…………………….67 圖4-2 文獻實驗與模擬之阻力係數隨攻角變化[10]…………………….67 圖4-3 疏、密網格模擬之升力隨攻角之變化……………………………68 圖4-4 疏、密網格模擬之阻力隨攻角之變化……………………………68 圖4-5 密網格與文獻模擬之升力係數隨攻角之變化…………………...69 圖4-6 密網格與文獻模擬之阻力係數隨攻角之變化…………………...69 圖4-7 λ=2,θ=0° 渦度分布[10]………………………………………70 圖4-8 λ=2,θ=0°,疏網格搭配DES+SST k-ω model渦度分布(左)、 壓力係數分布與流線(右)…………………………………………70 圖4-9 λ=2,θ=0°,密網格搭配DES+SST k-ω model渦度分布(左)、 壓力係數分布與流線(右)…………………………………………71 圖4-10 λ=2,θ=90°渦度分布[10]……………………………………..71 圖4-11 λ=2,θ=90°,疏網格搭配DES+SST k-ω model渦度分布(左)、 壓力係數分布與流線(右)…………………………………………72 圖4-12 λ=2,θ=90°,密網格搭配DES+SST k-ω model渦度分布(左)、 壓力係數分布與流線(右)…………………………………………72 圖4-13 λ=2,θ=180°渦度分布[10]……………………………………73 圖4-14 λ=2,θ=180°,疏網格搭配DES+SST k-ω model渦度分布(左)、 壓力係數分布與流線(右)…………………………………………73 圖4-15 λ=2,θ=180°,密網格搭配DES+SST k-ω model渦度分布(左)、 壓力係數分布與流線(右)…………………………………………74 圖4-16 λ=2,θ=270°渦度分布[10]……………………………………74 圖4-17 λ=2,θ=270°,疏網格搭配DES+SST k-ω model渦度分布(左)、 壓力係數分布與流線(右)…………………………………………75 圖4-18 λ=2,θ=270°,密網格搭配DES+SST k-ω model渦度分布(左)、 壓力係數分布與流線(右)…………………………………………75 圖4-19 λ=2,θ=0°,文獻模擬之Isovalues of Ux/U∞ 値分布圖[10]..76 圖4-20 λ=2,θ=0°,疏網格搭配DES+SST k-ω model之Isovalues of Ux 分布…………………………………………………………...76 圖4-21 λ=2,θ=0°,密網格搭配DES+SST k-ω model之Isovalues of Ux 分布…………………………………………………………...77 圖4-22 θ=0°,內部滑移區半徑分別為1.5m、2m壓力係數分布……77 圖4-23 θ=90°,內部滑移區半徑分別為1.5m、2m壓力係數分布…..78 圖4-24 θ=180°,內部滑移區半徑分別為1.5m、2m壓力係數分布…78 圖4-25 θ=270°,內部滑移區半徑分別為1.5m、2m壓力係數分布…79 圖4-26 θ=0°,內部滑移區半徑分別為1.5m、2m速度大小分布……79 圖4-27 θ=90°,內部滑移區半徑分別為1.5m、2m速度大小分布…..80 圖4-28 θ=180°,內部滑移區半徑分別為1.5m、2m速度大小分布…80 圖4-29 θ=270°,內部滑移區半徑分別為1.5m、2m速度大小分布…81 圖4-30 θ=0°,內部滑移區半徑分別為1.5m、2m渦度大小分布……81 圖4-31 θ=90°,內部滑移區半徑分別為1.5m、2m渦度大小分布…..82 圖4-32 θ=180°,內部滑移區半徑分別為1.5m、2m渦度大小分布…82 圖4-33 θ=270°,內部滑移區半徑分別為1.5m、2m渦度大小分布…83 圖4-34 θ=0°,內部滑移區半徑為3.75m疏、密網格之壓力係數分布….83 圖4-35 θ=90°,內部滑移區半徑為3.75m疏、密網格之壓力係數分….84 圖4-36 θ=180°,內部滑移區半徑為3.75m疏、密網格之壓力係數布 ...........................................................................................................84 圖4-37 θ=180°,內部滑移區半徑為3.75m疏、密網格之壓力係數分布 ……………………………………………………………………...85 圖4-38 θ=0°,內部滑移區半徑為3.75m疏、密網格之速度大小分布 …………………………………………………………………..….85 圖4-39 θ=90°,內部滑移區半徑為3.75m疏、密網格之速度大小分布 ……………………………………………………………………...86 圖4-40 θ=180°,內部滑移區半徑為3.75m疏、密網格之速度大小分布 ……………………………………………………………………...86 圖4-41 θ=270°,內部滑移區半徑為3.75m疏、密網格之速度大小分布 ……………………………………………………………………...87 圖4-42 θ=0°,內部滑移區半徑為3.75m疏、密網格之渦度大小分布 …………………………………………………………………..….87 圖4-43 θ=90°,內部滑移區半徑為3.75m疏、密網格之渦度大小分布 ……………………………………………………………………...88 圖4-44 θ=180°,內部滑移區半徑為3.75m疏、密網格之渦度大小分布 ……………………………………………………………………...88 圖4-45 θ=270°,內部滑移區半徑為3.75m疏、密網格之渦度大小分布 ……………………………………………………………………...89 圖4-46 觀測點位子……………………………………………………….89 圖4-47 分貝與音量關係表[29]…………………………………………..91 圖4-48 分貝隨頻率變化於A位置:滑移區域為(a)1.5m(b)2m(c)3.75疏 網格(d)3.75密網格………………………………………………...92 圖4-49 分貝隨頻率變化於B位置:滑移區域為(a)1.5m(b)2m(c)3.75疏 網格(d)3.75密網格………………………………………………...92 圖4-50 分貝隨頻率變化於C位置:滑移區域為(a)1.5m(b)2m(c)3.75疏 網格(d)3.75密網格………………………………………………...93 圖4-51 分貝隨頻率變化於D位置:滑移區域為(a)1.5m(b)2m(c)3.75疏 網格(d)3.75密網格………………………………………………...93 圖4-52 分貝隨頻率變化於E位置:滑移區域為(a)1.5m(b)2m(c)3.75疏 網格(d)3.75密網格………………………………………………...94 圖4-53 分貝隨頻率變化於F位置:滑移區域為(a)1.5m(b)2m(c)3.75疏網 格(d)3.75密網格…………………………………………………...94 圖4-54 分貝隨頻率變化於G位置:滑移區域為(a)1.5m(b)2m(c)3.75疏 網格(d)3.75密網格………………………………………………...95 圖4-55 分貝隨頻率變化於H位置:滑移區域為(a)1.5m(b)2m(c)3.75疏 網格(d)3.75密網格………………………………………………...95 圖4-56 分貝隨頻率變化於I位置:滑移區域為(a)1.5m(b)2m(c)3.75疏網 格(d)3.75密網格…………………………………………………...96 圖4-57 分貝隨頻率變化於J位置:滑移區域為(a)1.5m(b)2m(c)3.75疏網 格(d)3.75密網格…………………………………………………...96 圖4-58 風機交互作用下之P/Po圖(x/D,y/D)=(1.5,-1.5)…………………97 圖4-59 雙顆風機交互作用下之速度大小分布(x/D,y/D)=(1.5,-1.5)…... 98 圖4-60 雙顆風機交互作用下之渦度大小分布(x/D,y/D)=(1.5,-1.5)…... 99 圖4-61 雙顆風機交互作用下之壓力(P/Po)分布(x/D,y/D)=(2.0, 1.5)…100 圖4-62 雙顆風機交互作用下之速度大小分布(x/D,y/D)=(2.0, 1.5)…...101 圖4-63 雙顆風機交互作用下之渦度大小分布(x/D,y/D)=(2.0, 1.5)…...102 圖4-64 觀測點位置(x/D,y/D)=(1.5,-1.5)………………………………..103 圖4-65 觀測點位置(x/D,y/D)=(2.0,1.5)…………………………………103 圖4-66 各觀測點之分貝隨頻率變化(x/D,y/D)=(1.5,-1.5)……………..107 圖4-67 各觀測點之分貝隨頻率變化(x/D,y/D)=(2.0,1.5)………………109

    [1] Iida, A., Mizuno, A., and Fukudome, K,. “Numerical Simulation of Aerodynamic Noise Radiated from Vertical Axis Wind Turbines,” Proceedings of the 18th International Congress on Acoustics, 2004, CD-ROM.
    [2] Carlos Simao Ferreira, Gerard van Bussel, and Gijs van Kuik, “2D CFD Simulation of Dynamic Stall on a Vertical Axis Wind Turbine: Verification and Validation with PIV Measurements,” AIAA, 2007.
    [3] Carlos Simao Ferreira, Ferard van Bussel and Gijs van Kuik, “Simulating Dynamic Stall in a 2D VAWT: Modeling Strategy, Verification and Validation with Particle Image Velocimetry Data,” Journal of Physicsl, 2007.
    [4] Mazharul Islam, and David S-K. Ting, and Amir Fartaj, “Desirable Airfoil Features for Smaller-Capacity Straight-Bladed VAWT,” Wind Engineering, Vol. 31, pp. 165-196, 2007.
    [5] Mazharul Islam, and David S-K. Ting, and Amir Fartaj, “Design of a Special-Purpose Airfoil for Smaller-Capacity Straight-Bladed VAWT,” Wind Engineering, Vol. 31, pp. 401-424, 2007.
    [6] Bravo R., Tullis S., and Ziada S., “Performance Testing of a Small Vertical-Axis Wind Turbine,” 21st Canadian Congress of a Applied Mechanics, Ryerson University Toronto, Ontario Canada , June 3-7, 2007.
    [7] Mazharul Islam, David S-K. Ting, and Amir Fartaj, “Aerodynamic Models for Darrieus-Type Straight-Bladed Vertical Axis Wind Turbines. ” Renewable and Sustainable Energy Reviews, 2008.
    [8] Mazharul Islam, and Amir Fartaj, and Amir Fartaj, “Analysis of the Design Parameters Related to a Fixed-Pitch Straight-Bladed Vertical Axis Wind Turbine,” Wind Engineering, Vol. 32 , pp. 491-507, 2008.
    [9] Sylvain Antheaume, Thierry Maître, and Jean-Luc Achard, “Hydraulic Darrieus Turbines Efficiency for Free Fluid Flow Conditions versus Power Farms Conditions,” Renewable Energy, Vol. 33 , pp. 2186-2198, 2008.
    [10] Amet E., Maître T., Pellone C., and Achard J.-L.,“2D Numerical Simulations of Blade-Vortex Interaction in a Darrieus Turbine,” ASME, Vol. 131, 2009.
    [11] Shengyi Wang, Lin Ma, Derek B Ingham, Mohamed Pourkashanian, and Zhi Tao, “Turbulence Modelling of Deep Dynamic Stall at Low Reynolds Number,” Proceedings of the World Congress on Engineering Vol. II, 2010.
    [12] Akinari Shigetomi, Yuichi Murai, Yuji Tasaka, Yasushi Takeda, “Interactive Flow Field around Two Savonius Turbines,” Renewable Energy, Vol. 36, pp. 536-545, 2011
    [13] Habtamu Beri, and Yingxue Yao, “Double Multiple Stream Tube Model and Numerical Analysis of Vertical Axis Wind Turbine,” Energy and Power Engineering, Vol. 3, pp. 262-270, 2011.
    [14] Habtamu Beri, and Yingxue Yao, “Effect of Camber Airfoil on Self Starting of Vertical Axis Wind Turbine,” Journal of Environmental Science and Technology, Vol. 4, pp. 302-312, 2011,
    [15] John O. Dabiri, “Potential Order-of-Magnitude Enhancement of Wind Farm Power Density via Counter-Rotating Vertical-Axis Wind Turbine Arrays,” Journal of Renewable and Sustainable Energy, Vol. 3, pp. 043104-1-043104-12, 2011.
    [16] Kjellin J., Bulow F., Eriksson S., Deglaire P., Leijon M. and Bernhoff H., “Power Coefficient Measurement on a 12 kW Straight Bladed Vertical Axis Wind Turbine, ” Renewable Energy, Vol. 36, pp.3050-3053, 2011.
    [17] Marco Raciti Castelli, Alessandro Englaro and Ernesto Benini, “The Darrieus Wind Turbine : Proposal for a New Performance Prediction Model Based on CFD,” Energy, Vol. 36, pp. 4919-4934, 2011.
    [18] Huimin Wang, Jianliang Wang, Ji Yao, Weibin Yuan and Liang Cao, “Analysis on the Aerodynamic Performance of Vertical Axis Wind Turbine Subjected to the Change of Wind Velocity,” SciVerse ScienceDirect, Vol. 31, pp. 213-219, 2012.
    [19] Liang Cao, Huimin Wang, Jianliang Wang, Ji Yao, and Weibin Yuan, “Analysis on the Influence of Rotational Speed to Aerodynamic Performance of Vertical Axis Wind Turbine,” SciVerse ScienceDirect, Vol. 31, pp. 245-250, 2012.
    [20] Ji Yao, Jianliang Wang, Weibin Yuan, Huimin Wang, and Liang Cao, “Analysis on the Influence of Turbulence Model Changes to Aerodynamic Performance of Vertical Axis Wind Turbine.” SciVerse ScienceDirect, Vol. 31, pp. 274-281, 2012.
    [21] 韓尚育, ”小型垂直軸風機氣動力效率與噪音之數值研究,” 成功大學碩士論文, 2011.
    [22] 謝承翰, “垂直軸風力機扭力與功率的檢測與模擬,” 成功大學碩士論文, 2009.
    [23] FLUENT6.3 User’s Guide. 2006.
    [24] 陳廷彰, “風扇噪音分析與有孔平板受風扇噪音激振的聲場分析,” 中央大學碩士論文, 2008.
    [25] 王俊傑, “抽油煙機之流場模擬與氣動噪音分析,” 國立交通大學碩士論文, 2008.
    [26] FLUENT3.1 Tutorial Guide., “Aero-Noise Prediction of Flow Across a Circular Cylinder,” Fluent Inc., 2002.
    [27] World Wind Energy Association, “World Wind Energy Report,” WWEA, 2010.
    [28] Bureau of Energy , Ministry of Economic Affairs.
    http://www.moeaboe.gov.tw/
    [29] Anthony L., Rogers, James F., Manwell, Sally Wright,M. S , “Wind Turbine Acoustic Noise,” Renewable Energy Research Laboratory, 2006.

    下載圖示 校內:2014-08-22公開
    校外:2015-08-22公開
    QR CODE