1
Investigasi Numerik Dynamic Stall pada Vertical Axis Wind Turbine Airfoil Go Peter Christian Tejo Hutomo dan Herman Sasongko Teknik Mesin, Fakultas Teknologi Industri, Institut Teknologi Sepuluh Nopember (ITS) Jl. Arief Rahman Hakim, Surabaya 60111 Indonesia e-mail:
[email protected]
Abstrak—Pada operasinya, gerakan berputar dari VAWT mengakibatkan variasi angle of attack (α) yang besar pada blades selama rotor berevolusi. Hal ini mengakibatkan beban aerodinamis yang signifikan dan unsteady yang dikenal sebagai peristiwa dynamic stall. Simulasi numerik menggunakan model turbulensi SST k-ω digunakan untuk menginvestigasi fenomena fisik yang terjadi pada airfoil VAWT skala atap rumah yang beroperasi pada low tip speed ratio. Dalam penelitian ini digunakan airfoil dengan profil NACA 0012 yang beroperasi pada blade Reynolds number sebesar 1.35 x 105. Hasil kalkulasi CFD menunjukkan perbedaan yang mencolok pada gaya-gaya yang bekerja pada airfoil dibandingkan dengan kondisi statis dikarenakan adanya fenomena travelling vortices. Pembentukan leading edge dan trailing edge vortex yang menjadi karakteristik dari dynamic stall dapat diperlihatkan dengan jelas sehingga memberikan informasi yang jelas tentang perkembangan aliran fluida di sekitar blade dari VAWT. Adapun efek dynamic stall dalam peningkatan lift terlihat secara signifikan, posisi critical stall angle tertunda sekitar 10˚ dan lift yang diproduksi meningkat hingga dua kali lipat dari kasus statis. Hasil analisis medan aliran memperlihatkan bahwa stall terjadi akibat adanya insertion flow yang terjadi di dekat leading edge dari airfoil. Evolusi serta interaksi vortex pada fase upwind maupun downwind juga disajikan dalam makalah ini untuk memberikan pengetahuan lebih dalam tentang peristiwa dynamic stall pada VAWT. Kata Kunci— CFD, dynamic stall, unsteady aerodynamics, VAWT.
I. PENDAHULUAN
K
ebutuhan akan energi yang ramah dan terbarukan menjadi semakin besar seiring dengan meningkatnya kesadaran akan tingginya emisi karbon dan dampaknya terhadap lingkungan. Beberapa sumber alternatif untuk mengganti bahan bakar fosil telah diidentifikasi dan tenaga angin merupakan salah satu sumber energi yang paling menjanjikan [1]. Pertumbuhan ini disebabkan adanya ladang angin yang besar baik off-shore maupun on-shore serta minat yang kuat dalam memanen energi angin baik di daerah perkotaan maupun di pinggiran kota [2]. Saat ini terdapat dua kategori turbin angin modern, yaitu Horizontal Axis Wind Turbine (HAWT) dan Vertical Axis Wind Turbine (VAWT). Tipe HAWT lebih dominan secara komersial karena telah melalui berbagai riset dan pengembangan dalam beberapa dekade terakhir. Namun, hampir tidak mungkin untuk membangun pembangkit listrik tenaga angina skala besar dengan menggunakan HAWT karena kondisi lahan di wilayah perkotaan dan pesisir pantai. Vertical axis wind turbine (VAWT) telah dikenal sebagai teknologi yang tepat untuk pembangkitan energi dari tenaga
angin di perkotaan. Telah diketahui bahwa VAWT mengalami unsteady aerodynamics yang sangat kompleks [2]. Gerakan berputar dari turbin mengakibatkan variasi angle of attack (AoA) yang besar pada blades selama rotor berevolusi. Hal ini mengakibatkan beban aerodinamis yang signifikan dan unsteady yang dikenal sebagai peristiwa dynamic stall [3]. Dynamic stall adalah fenomena yang melibatkan serangkaian separasi dan reattachment aliran yang terjadi di sekitar lifting surfaces ketika dikenai gerakan unsteady. Fenomena ini belum benar-benar dipahami meskipun fenomena ini sangat penting untuk menentukan performa dan batas operasional dari helikopter, flapping wings, dan turbin angin. Dynamic stall diawali dengan terjadinya leading edge separation yang diikuti terbentuknya leading edge vortex (LEV). Leading edge vortex berkembang dan bergerak downstream sepanjang permukaan airfoil. Vortex mencapai trailing edge dan terurai diikuti dengan terbentuknya trailing edge vortex (TEV) yang mengindikasikan awal terjadinya stall. Trailing edge vortex kemudian akan terurai seiring dengan naiknya angle of attack diikuti dengan rusaknya leading edge vortex. Separasi leading edge vortex ini diikuti dengan penurunan lift secara dramatis dan peningkatan pitching moment secara mendadak. Sehingga peristiwa ini dapat menyebabkan terjadinya violent vibration dan beban aerodinamis yang besar dan sangat berbahaya bagi struktur airfoil. Pengaruh dynamic stall pada aerodinamika dan performa dari VAWT yang beroperasi pada low blade speed ratios telah banyak dipelajari. Beberapa simulasi CFD (Computational Fluid Dynamics) juga telah dilakukan untuk menangkap dan memahami medan aliran serta evolusi dari dynamic stall vortex. Sebagai contoh, Wang et al [4,5] menyimpulkan bahwa model turbulensi SST k-ω memberikan hasil prediksi yang lebih baik daripada model turbulensi Wilkox k-ω untuk kasus 2D. Dijelaskan juga bahwa tingkat turbulensi dari upstream mempengaruhi lokasi terjadinya lift stall. Meskipun demikian, penelitian-
penelitian ini dilakukan untuk menginvestigasi peristiwa dynamic stall pada airfoil yang berosilasi, dimana simulasi ini tidak mengakomodasi adanya plunging motion dari airfoil saat VAWT beroperasi [5]. Lebih lagi, VAWT beroperasi pada range angle of attack yang lebih luas dibandingkan dengan airfoil yang berosilasi maupun Harozontal axis wind turbine (HAWT) yang pada umumnya bekerja pada angle of attack yang positif dengan range yang lebih sempit sehingga menghasilkan pola dynamic stall yang berbeda. Oleh karena itu, studi yang dilakukan saat ini meliputi simulasi numerik berbasis CFD untuk menginfestigasi gaya-gaya dan medan aliran disekitar airfoil
2 dari single bladed Darrieus turbine skala atap rumah yang beroperasi pada Reynolds number yang rendah guna memahami fenomena fisik yang terjadi untuk pengembangan lebih lanjut dari turbin ini. II. METODOLOGI PENELITIAN Model CFD dengan 2D geometrical configuration digunakan untuk merepresentasikan domain analisa single bladed Darrieus turbine. Hal ini diambil berdasarkan literatur yang relevan [6,7] yang menunjukkan bahwa model 2D sudah cukup untuk mengungkap faktor-faktor yang mempengaruhi performa dan medan aliran di sekitar VAWT. Kontribusi dari blade end effects dan supporting arm junction effects dibaikan, tetapi hal ini dapat diterima karena dianggap sebagai efek sekunder. ANSYS Fluent 15.0 digunakan untuk semua simulasi yang dilakukan dalam penelitian ini. Kode yang digunakan merupakan metode finite volume untuk menyelesaikan persamaan-persamaan fluida. Model turbulensi k-ω SST dengan Kato-Launder Correction digunakan dalam kalkulasi ini untuk memprediksi turbulent production. Secara umum, kalkulasi numerik menggunakan model turbulensi k-ω SST menerapkan aliran fully turbulent untuk hasil yang lebih robust. Namun, kalkulasi fully turbulent akan menyebabkan gaya-gaya lebih rendah dibandingkan dengan kalkulasi dengan efek transisi dikarenakan adanya penebalan boundary layer. Hal ini diperkuat oleh Bangga et al [8] yang menyatakan bahwa dengan adanya surface roughness yang tinggi di dekat leading edge akan menyebabkan transisi aliran laminar ke turbulen yang lebih cepat atau prematur dimana hal ini akan menurunkan performa dari suatu airfoil. Maka dari itu, Intermittency crossflow transition model diaktifkan untuk meningkatkan akurasi simulasi turbin angin. Production Limiter digunakan untuk mengurangi turbulent kinetic energy yang berlebihan pada daerah stagnasi serta Low Reynolds Number Correction untuk meningkatkan akurasi pada daerah viscous sublayer.
konstan, oleh karena itu digunakan pressure-based solver. Persamaan incompressible, URANS diselesaikan untuk seluruh domain aliran. Second order implicit transient formulation dipilih untuk meningkatkan akurasi. Semua variabel solusi diselesaikan menggunakan second order upwind discretization karena sebagian besar aliran dapat diasumsikan tidak in line dengan mesh. Algoritma Coupled digunakan dalam pressure-velocity coupling algorithm dan semua persamaan untuk mendapatkan solusi yang diselesaikan secara simultan. Turbin yang diteliti memiliki diameter 1.2m, terdiri atas 1 blade dengan profil NACA 0012, panjang chord 0.15m dan blade pitch angle 0º (lihat Gambar 1). Meshing dibagi menjadi 2 domain, yaitu rotating domain dan stationary domain. Rotating domain terhubung dengan stationary domain via sliding interface boundary condition yang mengonservasi massa dan momentum. Untuk menangkap karakteristik secara akurat dari boundary layer dan gaya-gaya pada blade, ketinggian baris pertama pada elemen meshing di sekitar blade diatur agar memenuhi nilai non-dimensional wall distance (y+) kurang dari 1. Rotating domain dibuat 1c dari luar permukaan blade. Mesh dibuat dengan sangat halus dengan angle skewness dari tiap cell kurang dari 0.6. Untuk mengeliminasi efek wall, maka wall diletakkan sejauh 4D di atas dan 4D di bawah turbin. Untuk membuat kecepatan freestream yang seragam pada inlet maka posisi velocity inlet diletakkan sejauh 4D upstream. Posisi pressure outlet diletakkan sejauh 12D downstream dari turbin sehingga memungkinkan perkembangan wake secara penuh [5]. Adapun meshing yang digunakan dalam kalkulasi dan domain analisis ditunjukkan pada Gambar 2 dan 3.
Gambar 2. Meshing yang digunakan dalam kalkulasi
Gambar 3. Domain analisis R c θ
= turbine radius (0.6 m) = chord length (0.15 m) = azimuth angle
β = blade pitch angle (00) ω = kecepatan rotasi turbin u∞ = kecepatan freestream
Gambar 1. Geometri turbin Kasus yang ditinjau dalam penelitian ini merupakan aliran unsteady dengan kerapatan massa yang dianggap
Seperti yang telah dijelaskan sebelumnya, kecepatan relatif yang dialami blade saat blade berputar sepanjang azimuth menyebabkan perubahan angle of attack. Gambar 3 mengilustrasikan kecepatan aliran di sekitar blade VAWT pada posisi azimuth angle (θ). Azimuth angle diukur dari sumbu y negatif (titik referensi) berputar secara clockwise ke posisi dimana blade berada. AoA sendiri didefinisikan sebagai sudut yang dibentuk antara chord line dari airfoil
3 dengan undisturbed relative incoming flow. Secara vektor dan geometri (lihat Gambar 4) dapat ditentukan hubungan antara azimuth angle dengan AoA seperti Persamaan 1.
Rec = 1.35 x 105
Gambar 4. Ilustrasi vektor kecepatan pada Darrieus turbine blade tan α =
U∞ sinθ ωR+u∞ cosθ
α = arctan (
=
NACA 0012
sinθ λ+cosθ
sinθ λ+cosθ
)
(1)
Gambar 5. Validasi model numerik dengan data eksperimen dan hasil kalkulasi XFOIL untuk profil NACA 0012 pada Re = 1.35 x 105
Untuk memvalidasi simulasi numerik yang dilakukan, maka kalkulasi statis untuk NACA 0012. Setting kalkulasi yang digunakan untuk kasus ini juga sama dengan yang digunakan dalam simulasi turbin, hanya saja karena kasusnya statis maka intermittency crossflow transition model dan sliding mesh tidak digunakan. Kalkulasi yang dilakukan merupakan kasus fully turbulent flow sesuai dengan eksperimen yang dilakukan oleh Lee dan Gerontakos [13]. Adapun kasus yang ditinjau dalam penelitian ini menerapkan aliran steady untuk angle of attack yang rendah (0˚ ≤ α ≤ 10˚) dan aliran unsteady untuk angle of attack yang tinggi (α > 10˚) karena tingginya tingkat unsteadiness pada angle of attack yang tinggi. Sebagai pembanding, kalkulasi dengan XFOIL boundary layer code juga dilakukan. III. HASIL DAN DISKUSI Karakteristik dynamic stall dari single bladed Darrieus turbine telah dilakukan secara numerik. Meshing yang digunakan dalam studi ini merupakan hybrid mesh yang mengombinasikan structured grid di dekat airfoil dan unstructuredgrid untuk domain fluida. Grid convergence study telah dilakukan dan meshing yang dipilih terdari dari total cells sebanyak 304476. Jumlah grid pada suction side dan pressure side dari airfoil adalah sebanyak 250 nodes dengan ketinggian baris pertama dari elemen meshing sebesar 5x10-4 c untuk memprediksi fenomena pada laminar sublayer secara akurat. Beberapa time step size (Δt) yang ekivalen dengan specific rotational displacement sepanjang azimuth juga dievaluasi. Time step size yang dipilih adalah sebesar 1˚ ω-1 yang telah mampu menangkap vortex evolution dari VAWT dengan baik. Adapun periodic convergence tercepai setelah 4 kali revolusi turbin. Lift coefficient dimonitor sepanjang 7 revolusi dan tidak ditemukan perbedaan setelah periodic convergence tercapai. Oleh karena itu, data yang dipresentasikan dalam makalah ini adalah data yang diperoleh dari 3 revolusi terakhir. Gambar 5 menunjukkan perbandingan antara kalkulasi CFD yang dilakukan dalam studi ini dengan eksperimen dari Lee dan Gerontakos [9] dan kalkulasi XFOIL code. Terlihat bahwa model numerik yang dibuat menghasilkan hasil yang lebih baik daripada hasil kalkulasi XFOIL dan sangat sesuai dengan data eksperimen.
Gambar 6. Perbandingan lift coefficient pada kondisi statis dan dinamis
Gambar 7. Ilustrasi symmetrical airfoil dalam rectilinear dan curvelinear flowfield. Gambar 6 menunjukkan perbandingan antara lift coefficient (Cl) pada kondisi statis dan dinamis. Lift coefficient yang dihasilkan oleh kondisi dinamik lebih tinggi dibandingkan dengan kondisi statis, baik dalam kondisi upstroke maupun downstroke. Posisi terjadinya stall pada VAWT juga tertunda sekitar 10 derajat dari α ≈ 13° menjadi sekitar 23°. Diperkirakan curvelinear motion dari aliran berkontribusi terhadap peristiwa ini. Terjadinya curvelinear flow pada VAWT menyebabkan kelakuan aerodinamis di sekitar airfoil berubah mengakibatkan tekanan yang lebih besar pada radius terluar dari kelengkungan aliran (lihat
4 Gambar 7) sehingga symmetrical airfoil menghasilkan lift pada AoA = 0˚ (lihat Gambar 6). Curvelinear flow juga menyebabkan tertundanya separasi pada suction side dari airfoil dan juga tertundanya dynamic stall vortex sehingga lift yang dihasilkan lebih besar. Gambar 8 menunjukkan perbandingan medan aliran untuk kasus statis dan dinamis pada fase upstroke. Dapat dilihat bahwa separasi massif yang terjadi pada kondisi statis tidak terjadi pada kondisi dinamis.
Static
Dynamic
α = 0˚
TEV bersama LEV yang telah terseparasi bergerak bersama menuju downstream turbin membentuk vortex pair. Vortex pair ini menghasilkan tekanan yang tinggi pada suction side dari airfoil, sehingga menyebabkan turunnya lift lebih lanjut. Ketika vortex pair ini bergerak menuju downstream turbin, mereka akan berinteraksi dengan blade pada fase downwind (lihat Gambar 9). Terjadinya blade vortex interaction ini menyebabkan terjadinya fluktuasi lift coefficient pada fase downwind dengan rata-rata yang bernilai rendah [10,11]. Hal ini diperkuat oleh pernyataan Bangga [12] yang menjelaskan bahwa perbedaan range operasi angle of attack dan reduced frequency dapat menghasilkan force polar yang berbeda meskipun bekerja pada Reynolds number yang sama. θ = 64˚ α = 20.2˚
θ = 69˚ α = 21.6˚
θ = 72˚ α = 23.4˚
θ = 74˚ α = 22.9˚
θ = 81˚ α = 24.6˚
θ = 85˚ α = 25.5˚
α = 0˚
Dynamic Stall
α = 13.5˚
α = 13.5˚ θ = 88˚ α = 26.2˚
α = 30˚
α = 30˚
θ =94˚ α = 27.3˚
θ = 97˚ α = 27.9˚
Gambar 10. Proses terjadinya dynamic stall
Gambar 8. Vorticity field dan streamlines untuk upstrokeupwind phase. Clockwise dan counter clockwise vorticity ditunjukkan dengan warna biru dan merah. Seluruh vorticity contour yang ditampilkan diplot dengan contour levels yang sama.
θ = 150˚ α = 24˚
θ = 180˚ α = 0˚
θ = 210˚ α = -18.6˚
Material from the Outside of LEV Clockwise TEV pocket Counter Clockwise TEV
Gambar 11. Insertion dan counter clockwise TEV yang menyebabkan dynamic stall (α = 22.90°)
θ = 240˚ α = -30˚
θ = 270˚ α = -26.6˚
θ = 300˚ α = -19.1˚
Gambar 9. Blade vortex interaction Hasil yang sangat berbeda ditunjukkan pada fase downwind dimana negative lift yang dihasilkan oleh VAWT lebih rendah daripada kondisi statis. Hal ini terjadi karena saat blade memasuki fase downwind kecepatan aliran lokal yang dialami blade lebih rendah daripada saat fase upwind (karena angin bertiup dari downstream airfoil). Lebih lagi, Pada awal fase downstroke, terlihat bahwa counter clockwise
Gambar 10 menunjukkan proses terjadinya dynamic stall pada VAWT. Berbeda dengan proses terjadinya dynamic stall pada umumnya yang diawali dengan terbentuknya LEV, dynamic stall pada kasus ini diawali dengan terbentuknya clockwise TEV. TEV inilah yang berkontribusi menambah suction pada bagian suction side dari airfoil, sehingga lift yang dihasilkan terus meningkat. Clockwise TEV berkembang seiring naiknya angle of attack kemudian barulah terbentuk LEV. Perkembangan LEV menyebabkan naiknya lift lebih lanjut. Dynamic stall terjadi saat terbentuknya counter clockwise TEV yang diiringi dengan lepasnya LEV dari permukaan airfoil akibat adanya insertion flow yang berasal dari luar kantung LEV yang kemudian membentuk vortex yang berputar counter clockwise dan mendesak LEV hingga lepas dari permukaan airfoil (lihat Gambar 10). LEV yang telah lepas dari permukaan airfoil ini kemudian berkembang dan terkonveksi menuju trailing edge.
5 Akibat naiknya angle of attack, LEV mendesak clockwise TEV hingga akhirnya bertabrakan dengan counter clockwise TEV. Counter clockwsie TEV pecah dan ter-shedding, sedangkan clockwise TEV terhisap oleh LEV yang kemudian terus berkembang dan meninggalkan airfoil. Hal ini menyebabkan turunnya lift coefficient secara drastis. Fenomena perkembangan aliran fluida pada operasi VAWT tipe Darrieus telah didiskusikan dalam makalah ini. Hasilnya menunjukkan bahwa komputasi numerik mampu menangkap data kualitatif berupa fenomena travelling vortex yang merupakan karakter utama dari dynamic stall. Disimpulkan juga bahwa insertion flow di dekat leading edge memainkan peranan penting dalam karakteristik aliran di sekitar airfoil dari VAWT, sehingga desain bentuk leading edge dari airfoil sangat penting dalam menentukan performa dari VAWT. IV. KESIMPULAN/RINGKASAN Simulasi numerik menggunakan URANS telah dilakukan untuk menginvestigasi terjadinya peristiwa dynamic stall pada VAWT tipe Darrieus. Disimpulkan bahwa pola dynamic stall dipengaruhi oleh perbedaan range operasi angle of attack, Reynolds number serta reduced frequency. Pada kasus ini, dynamic stall diawali dengan terbentuknya clockwise TEV yang kemudian diikuti dengan terbentuknya LEV. Hal ini mengakibatkan kenaikan lift secara signifikan serta penundaan terjadinya stall sebesar 10˚ yaitu dari angle of attack 13˚ pada kondisi statis, menjadi 23˚ pada operasi VAWT. Stall terjadi ketika LEV lepas dari permukaan airfoil akibat desakan insertion flow seiring dengan terbentuknya counter clockwise TEV. Blade vortex interaction yang terjadi pada fase downwind mengakibatkan berfluktuasinya lift coefficient dengan rata-rata yang bernilai rendah. Namun, karena studi yang dilakukan adalah untuk kasus 2D maka diperlukan kalkulasi 3D maupun eksperimen untuk memvalidasi hasil yang telah didapatkan. UCAPAN TERIMA KASIH Penulis mengucapkan terimakasih kepada Prof. Dr.-ing Herman Sasongko dan Galih Senja Titah Aji Bangga, S.T, M.T yang telah membimbing penulis, serta Laboratorium Mekanika Fluida dan Pencampuran, Jurusan Teknik Kimia, ITS, yang memfasilitasi penulis untuk melakukan simulasi dengan software ANSYS Fluent 15.0. DAFTAR PUSTAKA [1]
[2]
[3]
[4]
[5]
Danao, L.A.M., Edward, J., Eboibi, O., Howell, R. 2013. The Performance of a Vertical Axis Wind Turbine in Fluctuating Wind-A Numerical Study. Proceedings of the World Congress on Engineering 2013 Vol III, WCE 2013, July 3 - 5, 2013, London, U.K. Qin, N., Howell, R., Durrani, N., Hamada, K., Smith, T. 2011. Unsteady Flow Simulation and Dynamic Stall Behaviour of Vertical Axis Wind Turbine Blades. Wind Engineering Vol. 35, no. 4, 2011 PP 511-510 Scheurich, Frank. 2011. Modelling the Aerodynamics of Vertical Axis Wind Turbine. Ph.D Thesis of the School of Engineering, University of Glasgow, Glasgow, U.K. Wang, S., Ingham, D.B., Ma, L., Pourkashanian, M., Tao, Z. 2012. Turbulence Modeling of Deep Dynamic Stall at Relatively Low Reynolds Number. Journal of Fluid and Structures 33 (2012) pp.191209. Wang, S., Ingham, D. B., Ma, L., Pourkashanian, M., Tao, Z. 2010. Numerical Investigations on Dynamic Stall of Low Reynolds Number Flow around Oscillating Airfoils. Computers & Fluids 39, 1529–1541.
[6]
Dai, Y. M., and Lam, W. 2009. Numerical Study of Straight-Bladed Darrieus-Type Tidal Turbine. Proceedings of the Institution of Civil Engineers - Energy, 162(2), pp. 67 -76. [7] Consul, C. A., Willden, R. H. J., Ferrer, E., and Mcculloch, M. D. 2009. Influence of Solidity on the Performance of a Cross-Flow Turbine. Proceedings of the 8th European Wave and Tidal Energy Conference, Uppsala, Sweden. [8] Bangga, G.S.T.A., Lutz, Th., Krämer, E. 2015. An Examination of Rotational Effects on Large Wind Turbine Blades. 11th EAWE PhD Seminar, Stuttgart, Germany. [9] Lee, T., and Gerontakos, P. 2004. Investigation of Flow over an Oscillating Airfoil. Journal of Fluid Mechanics, 512, pp. 313-341. [10] Simão Ferreira, C. J.. 2007. The Near Wake of the VAWT, 2D and 3D Views of the VAWT Aerodynamics. Master Thesis of Aerospace Engineering Department, TU Delft, The Netherlands. [11] Tsai, H. C., Colonius, T. 2014. Coriolis Effect on Dynamic Stall in a Vertical Axis Wind Turbine at Moderate Reynolds Number. 32nd AIAA Applied Aerodynamics Conference,Atlanta, USA. [12] Bangga, G.S.T.A. 2013. Assessment of Modified Two-Equations URANS Turbulence Model to Predict the Onset of Dynamic Stall. Master Thesis of Mechanical Engineering Department, Institut Teknologi Sepuluh Nopember.
