Pusat Penelitian Informatika - LIPI
PEMBUATAN PROGRAM SIMULASI PREDIKSI PERILAKU DINAMIK SISTEM DUALROTOR DENGAN MENGGUNAKAN SOFTWARE MATLAB Tri Admono Pusat Penelitian Telimek -LIPI
ABSTRACT This thesis presents the discussion of the making of the dynamic behavior prediction simulation of dualrotor system by using MATLAB software. Dualrotor system has two rotors wherethe axis of both rotors is in one line and rotate together. Both the rotors of the dualrotor system are connected by intershaftbearing. Dualrotor system are usually found in compressor and turbine of air craft engine. The operation of rotary machine at safe speed is a great importance, because when operated close to its natural frequency, the stability of the machine will be disturbed. The prediction of the dynamic behavior of dualrotor system will be done by using rotor dynamic knowledge. The behavior state in this thesis related to natural frequency, mode shape, and vibration respon.In analyzing the dynamic behavior of rotor system, it is assumed that rotor shaft is only influenced by bending, so that the effect of torsion is negligible. Bearing and seals aren’t have damping and stiffness in x and z direction. The axial force are neglected .
The result of simulation
program are displayed on graphics that show the influence of every component of rotor shaft to the dynamic behavior. It is also can be seen that the influence of the changing of stiffness and damping to the dynamic behavior
ABSTRAK Dalam tulisan ini akan membahas pembuatan program simulasi prediksi perilaku dinamik sistem dualrotor dengan menggunakan software MATLAB. Sistem dualrotor memiliki dua buah rotor dengan kedua sumbunya saling berimpit dan mengalami rotasi secara bersamaan. Antara kedua poros sistem dualrotor ini dihubungkan oleh intershaftbearing. Sistem dualrotor biasanya ditemui pada kompresor atau turbin pesawat terbang. Pengoperasian mesin-mesin rotasi jenis ini pada kecepatan putar yang aman menjadi suatu masalah yang sangat penting, mengingat apabila dioperasikan disekitar frekuensi pribadinya akan sangat mengganggu kestabilannya. Prediksi perilaku dinamik sistem
Kedeputian Ilmu Pengetahuan Teknik
1
Bandung, 29 – 30 Juli 2003
dualrotor tersebut akan dilakukan dengan menggunakan ilmu dinamika rotor (rotordynamics). Perilaku yang dimaksud adalah frekuensi pribadi, mode shape, dan respon getaran. Dalam menganalisis perilaku dinamik system rotor ini diasumsikan poros rotor hanya dipengaruhi oleh beban bending, sehingga momen torsi yang timbul akibat putaran poros diabaikan. Bantalan dan seals dianggap tidak mempunyai redaman dan kekakuan hanya dalam arah x dan z. Gaya aksial pada rotor juga diabaikan. Hasil dari program simulasi ini berupa grafik yang memperlihatkan pengaruh dari komponen poros rotor pada perilaku dinamik. Disini kita dapat melihat pengaruh perubahan redaman dan kekakuan pada perilaku dinamik tersebut.
LATAR BELAKANG Penggunaan mesin-mesin rotasi telah merambah ke berbagai bidang penerapan, antara lain: pada stasiun tenaga, sistem propulsi, mesin pesawat terbang, machine tools, automobil, peralatan medis, peralatan rumah tangga dan lain sebagainya, sehingga dapat dikatakan dari berbagai macam mesin
sangat sulit menemukan mesin yang tidak melibatkan
komponen rotasi. Sistem dual rotor adalah suatu sistem rotor yang memiliki dua buah rotor dimana kedua sumbu saling berimpit mengalami rotasi bersamaan dengan kecepatan putar sama atau berbeda dengan arah putaran sama atau berlawanan dan dihubungkan oleh bantalan antar poros (intershaft bearing) sebagaimana gambar 3.3. Sistem dual rotor biasanya ditemui pada kompresor atau turbin pesawat terbang, yaitu terdiri dari low pressure (LP) turbin atau kompresor dan high pressure (HP) turbin atau kompresor. Karakteristik rotor 1 dan rotor 2 saling berhibungan dan saling mempengaruhi (couple). Pada awal permulaan abad ke-20, perancangan mesin rotasi telah mencapai tahap dimana plant memerlukan mesin rotasi yang bisa berputar dengan kecepatan tinggi, dalam banyak kasus telah ditemui timbulnya satu atau lebih kecepatan kritis pada mesin rotasi [1]. Putaran poros rotor yang relatif tinggi pada mesin-mesin rotasi, disertai dengan beban torsi maupun bending berakibat timbulnya getaran (vibration) pada mesin tersebut, jika putaran mesin tersebut tidak berada dalam rentang kecepatan operasi yang dianjurkan atau katakanlah sampai masuk pada daerah putaran kritisnya maka akan terjadi resonansi. Hal
2
Pemaparan Hasil Litbang 2003
Pusat Penelitian Informatika - LIPI
ini akan mengakibatkan berbagai kerusakan (failure) antara lain pada bantalan, defleksi rotor yang berlebihan atau mesin bisa mengalami shutdown secara tiba-tiba.(2) Pada makalah ini analisis sistem poros rotor mengambil asumsi hanya dipengaruhi oleh beban bending, sehingga torsi yang timbul akibat putaran poros diabaikan [3]. Bantalan dan seals tidak mempunyai redaman dan sedangkan kekakuan hanya dalam arah x dan z saja. Gaya aksial pada rotor diabaikan. Radius poros pertama (dianggap pejal) sama dengan radius dalam disk. Rotor pertama dan kedua untuk dualrotor ini berada pada sumbu yang sama (sesumbu).
KARAKTERISTIK ELEMEN SISTEM POROS ROTOR [3] Elemen dasar dari sistem poros rotor adalah: disk, poros, bantalan dan seals. Persamaan energi kinetik diperlukan untuk mendapatkan karakteristik dari poros, rotor dan ketakseimbangan massa. Persamaan energi regangan juga diperlukan untuk mendapatkan karakteristik dari poros. Persamaan energi kinetik dan energi regangan untuk elemenelemen dasar dari sistem poros rotor tersebut dapat dilihat pada buku refernsi.(3) Persamaan Lagrange diterapkan untuk mendapatkan persamaan gerak dari sistem poros rotor. d ∂T dt ∂q io
dengan
i
∂T ∂U − ∂q + ∂q = Fq i i i
adalah jumlah derajad kebebasan (dof), qi adalah koordinat umum yang
independen, Fqi adalah gaya-gaya luar yang digeneralisasi.
MODEL ELEMEN HINGGA SISTEM DUALROTOR Karakteristik Model Elemen Hingga Sistem Dualrotor Untuk keperluan analisis sistem poros rotor dengan jumlah derajad kebebasan banyak (multi degree of freedom) maka perlu disusun karakteristik elemen-elemen rotor berdasarkan metode elemen hingga. Proses penyusunan elemen-elemen rotor dapat dilihat pada buku referensi.(3) Pada sistem dualrotor karakteristik rotor 1 dan rotor 2 saling berhubungan dan saling mempengaruhi (couple). Putaran untuk masing-masing rotor dapat dinyatakan :
Kedeputian Ilmu Pengetahuan Teknik
3
Bandung, 29 – 30 Juli 2003
Ω1 = aΩ 2 + b
dimana
a dan b adalah
variabel konstan, dalam aplikasi kerekayasaan nilai-ailai tersebut
ditentukan dari eksperimen. Hubungan karekteristik poros 1 dan poros 2 dinyatakan pada gambar 3.3, dimana w1 perpindahan
pada rotor satu dan
u2 dan w2 perpindahan
u1 dan
pada rotor dua yang berubungan
dengan titik yang sama pada sumbu rotor. Bantalan antar poros dianggap pendek, pengaruh slope diabaikan. Ambil harga R1 dan R1 bisa
reaksi rotor 1 dan 2 pada bantalan. Maka
ditulis dalam bentuk matrik:
R R1 = 1x R1z
maka:
R2 adalah
R1x = [ Fu1, Fu 2 ]t
k xz u1 − u2 cxx + k zz w1 − w2 czx
k R1 = xx k zx
R1z = [ Fw1, Fw2 ]t
cxz u&1 − u&2 czz w&1 − w& 2
Persamaan di atas dapat ditulis dalam bentuk ringkas sebagai: u −u u& − u& R1 = K 1 2 + C 1 2 − w w 2 1 w&1 − w& 2
k K = xx k zx
dengan :
k xz k zz
c C = xx czx
cxz czz
Karena massa bantalan diabaikan : R1 + R2 = 0
Akhirnya persamaan dapat disusun menjadi: R1 K = R2 − K
u1 − K w1 C + K u2 − C w2
u&1 − C w&1 & C u2 w& 2
Penerapan persamaan Lagrange pada rotor 1 dan rotor 2 menghasilkan persamaan : Rotor 1 :
M1δ&&1 + C1 (Ω1 )δ&1 + K1δ1 = F1 (t ) − R1
Rotor 2 :
M 2δ&&2 + C2 (Ω 2 )δ&2 + K 2δ 2 = F2 (t ) − R2
Karena persamaan-persamaan tersebut saling menghilangkan maka persamaan gerak sistem dual rotor dapat diidentikkan dengan persamaan monorotor sebagai berikut : Mδ&& + C (Ω1, Ω2 )δ& + Kδ = F (t ) ,
dimana M adalah matrik massa simetrik, C adalah matrik asimetrik dan kekakuan asimetrik dan K
δ
K
adalah matrik
adalah vektor perpindahan nodal. Sedangkan matrik M , C dan
secara berurutan adalah assembly dari matrik global M1 , M2 ; C1(Ω1) , C2(Ω2) serta redaman
viskos
4
C,
kemudian
K1 , K 2 dan
kekakuan bantalan
K
.
Pemaparan Hasil Litbang 2003
Pusat Penelitian Informatika - LIPI
Gambar 1 menunjukkan suatu hubungan antara poros pertama dengan poros kedua yang dilakukan oleh intershaftbearing.
intershaftbearing
Gambar 1 Model sistem dual rotor
Pengglobalan pada sistem dual rotor Pengglobalan matrik pada sistem dual rotor pada dasarnya hampir mirip dengan sistem monorotor berbeda hanya karena ada dua buah poros pada sistem dual rotor. Jadi masingmasing poros harus dibuat terlebih dahulu matrik globalnya baru kemudian dari dua matrik global tersebut disusun matrik global sistem dual rotor. Pemodelan bantalan antar poros (intershaft bearing) juga harus diperhatikan, sebab menyangkut pengisian matrik krakteristik bantalan antar poros pada matrik global dual rotor. Pengglobalan matrik untuk masing-masing poros rotor telah dijelaskan di atas, kemudian susunan matrik global untuk dual rotor (misalkan matrik massa global) ditunjukkan pada gambar 2. Kemudian jika massa bantalan diabaikan susunan matrik kekakuan global dan matrik redaman global untuk sistem dual rotor berbeda dengan matrik massa globalnya, karena harus memasukkan pemodelan bantalan antar poros pada matrik tersebut. Perbedaan tersebut ditunjukkan pada gambar 3.
Kedeputian Ilmu Pengetahuan Teknik
5
Bandung, 29 – 30 Juli 2003
Cara pengisian matrik karakteristik elemen-elemen rotor pada matrik global dual rotor juga berbeda dengan sistem monorotor khususnya pengisian untuk rotor luar. Hal ini disebabkan karena pemodelan metode elemen hingga tidak kontinyu untuk memodelkan interaksi antara dual rotor jadi ada satu elemen yang hilang (missing element) pada pemodelan tersebut. Salah satu caranya adalah dengan memisalkan rotor dalam (inner rotor) dibagi menjadi N nodal, sedangkan rotor luar (outer rotor) dibagi menjadi M nodal, maka model elemen hingga diberi penomoran nodal secara berurutan jadi seluruhnya terdapat
N + M = P nodal.
Dengan mengambil jumlah derajad kebebasan tiap nodal empat
maka ukuran matrik global dual rotor adalah 4 P . Sebagai contoh akan dimasukkan matrik karakteristik massa disk rotor luar ke dalam matrik global dual rotor, dari model elemen hingga disk diletakkan pada nodal d maka diambil suatu harga dum = ((d − N ) − 1) * 4
selanjutnya merujuk pada letak elemen-elemen matrik karakteristik disk, maka di dalam matrik global letaknya pada (baris, kolom) sebagai berikut (dum + 1, dum + 1)
adalah letak elemen
MD
(dum + 2, dum + 2)
adalah letak elemen
MD
(dum + 3, dum + 3)
adalah letak elemen
I Dx
(dum + 4, dum + 4)
adalah letak elemen
I Dx
Demikian juga untuk memasukkan matrik karakteristik bantalan ke dalam matrik global kekakuan dan redaman bisa dipakai cara yang sama. Keuntungan dari metode ini adalah sangat mudah diterapkan pada program komputer.
6
Pemaparan Hasil Litbang 2003
Pusat Penelitian Informatika - LIPI
1 1 4 5 8 9 12
4 5
8 9
12
4(M+N)
Matriks global Poros 1
1 2
4N N-2 N-1
Matriks global Poros 2
1 2
M-2 M-1
4(M+N)
Gambar 2 Matriks massa global system dualrotor 1
45
8 9
4(b-1)
4(M+N)-3
4(M+N)
12 1 4 5 8 9 12
Matriks global poros 1
1 2
Intershaft bearing
Intershaft bearing
Matriks lokal Poros 1 (b-1)4
Matriks lokal Poros 2 Matriks global Poros 2
Intershaft bearing
M-2
4(M+N)-3
M-1
4(M+N)
Gambar 3 Matriks kekakuan dan redaman global system dualrotor
Kedeputian Ilmu Pengetahuan Teknik
7
Bandung, 29 – 30 Juli 2003
SOLUSI PERSAMAAN (3) Solusi persamaan dicari untuk mendapatkan karakteristik perilaku dinamik diantaranya: • Frekuensi pribadi sebagai fungsi putaran dalam bentuk diagram Campbell dan daerah (zone) instabilitas. • Respon gaya eksitasi khususnya yang dibahas di sini adalah eksitasi ketakseimbangan massa terhadap sistem poros rotor. Pada umumnya sistem persamaan memiliki derajad tinggi, untuk pemecahannya bisa menggunakan metode langsung (direct method) atau bisa juga dengan metode pseudo modal (pseudo modal method). Dalam makalah ini digunakan metode pseudo modal (pseudo modal method).
PROGRAM SIMULASI PREDIKSI PERILAKU DINAMIK SISTEM DUALROTOR Pemecahan persamaan diatas diselesaikan dengan software MATLAB, sedangkan dalam program simulasinya menggunakan fasilitas yang dimiliki oleh MATLAB yaitu graphical user interface (GUI). Dibawah ini ditampilkan modul program simulasi dinamika rotor jenis dualrotor.
Gambar 4 Tampilan modul dualrotor
8
Pemaparan Hasil Litbang 2003
Pusat Penelitian Informatika - LIPI
Gambar 5 Model elemen hingga sistem dual rotor [3] Model elemen hingga sistem dual rotor Dari gambar 5 terlihat bahwa model dualrotor berdasarkan literature[3] terdiri dari 7 elemen beam untuk poros dalam (inner shaft) dan 4 elemen beam untuk poros luar (outer shaft). Untuk lebih jelasnya dapat dilihat pada gambar 4. Dimensi untuk penampang poros sistem dual rotor adalah sebagai berikut: poros dalam (1)
: jari-jari =
poros luar (2) : jari-jari dalam = jari-jari luar
=
1,524 cm 2.54 cm 3,048 cm
kedua poros berputar searah dengan perbandingan kecepatan angular antara poros 2 dan poros 1 sebesar 1,5. Data disk untuk sistem dual rotor diperlihatkan pada tabel 1. Tabel 1 Data disk dual rotor[3] Disk Massa (kg) -2
2
D1
D2
D3
D4
10,51
7,01
3,5
7,01
IDx 10 (kg m )
4,295 2,145 1,355
3,39
IDy 10-2 (kg m2)
8,59
6,78
4,29
2,71
Karena data dimensi disk-disk tidak diketahui maka penulis memasukkan datanya seperti pada isi gambar 4 (default-nya). Sistem poros rotor terdiri dari empat buah bantalan identik dengan karakteristik: Tabel 2 Karakteristik bantalan dual rotor[3] Bantalan kxx=kzz
1
2
3
4
2.63x107 1.75x107 0.875x107 1.75x107
kxz=kzx
Kedeputian Ilmu Pengetahuan Teknik
0
0
0
0
9
Bandung, 29 – 30 Juli 2003
Panjang tiap elemen poros dinyatakan dalam bentuk koordinat nodal untuk masing-masing poros sebagai berikut: Tabel 3a Koordinat nodal poros 1 Nodal
1
2
3
4
5
6
7
8
Absis(cm) 0 7,62 15,875 24,13 32,385 40,64 45,72 50,8
Tabel 3b Koordinat nodal poros 2 Nodal
9
10
11
12
13
Absis(cm) 15,24 20,32 27,94 35,56 40,64 Prediksi perilaku dinamik sistem poros rotor jenis dualrotor dilakukan dengan mengisi data input seperti tampak pada gambar 4. Data input default yang ada mengacu pada literatur [3]. Dari data yang ada apabila kemudian di-klik pada tombol “RUN” maka akan dapat ditampilkan prediksi perilaku dinamiknya. Diagram Campbell Dengan meng-klik tombol “Campbell” kita dapat melihat tampilan diagram Campbell, tampak seperti gambar 6. Disini ditampilkan sepuluh frekuensi pribadi terendah pertama sebagai fungsi dari putaran poros rotor. Respon ketidakseimbangan massa Sedangkan apabila di-klik tombol “Campbell” maka akan tampil diagram respon ketidakseimbangan massa, seperti pada gambar 7. Sebuah massa tak seimbang sebesar 70 g mm diasumsikan berada pada sistem poros rotor. Prediksi respon ketidakseimbangan massa dilakukan dengan menganggap ketidakseimbangan massa terletak pada disk 4 nodal 7. Lokasi respon yang dipilih yaitu pada nodal 7 dan nodal 12.
10
Pemaparan Hasil Litbang 2003
Pusat Penelitian Informatika - LIPI
Diagram Campbell
400
10
350
F=N/60 10
300
-3
6949,35
250
F=0,5N/60
10
Amplitudo (m)
F(Hz)
Respon Massa Unbalance
-2
200 150
10
100 10
11.980
-4
nodal12
-5
nodal 7
-6
50 0
0
2000
4000
6000 8000 N(rpm)
10000
12000
10
14000
-7
0
2000
Gambar 6 Diagram Campbell
4000 6000 8000 10000 Kecepatan putar, N (rpm)
12000
14000
Gambar.7. Respon Ketidak seimbangan massa
Modeshape Apabila di-klik pada tombol “Modeshape 1” maka akan tampil gambar modeshape 1, begitu juga dengan modeshape yang lain. Pada gambar 8 ditampilkan enam modeshape pertama pada kecepatan 14000 rpm. Modeshape 1
Modeshape 2
-4
-4
x 10
x 10
1.5
1.5
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
-1.5 2
-1.5 2
1 x 10
1
8 6
0
-4
-1 -2
2 0
4 F=63,3012 Hz (BW)
8 6
0
-4
x 10
-1 -2
Modeshape 3
F=139,4209 Hz (FW)
-5
x 10
x 10
6
6
4
4
2
2
0
0
-2
-2
-4
-4
-6 1
-6 1
0.5 x 10
2 0
Modeshape 4
-5
-4
4
0.5
8 6
0
4
-0.5 -1
2 0
F=173,5023 Hz (BW)
Kedeputian Ilmu Pengetahuan Teknik
x 10
-4
8 6
0
4
-0.5 -1
2 0
F=201,7329 Hz (FW)
11
Bandung, 29 – 30 Juli 2003
Modeshape 6
Modeshape 5 -5
-5
x 10
x 10 6
6
4
4
2
2
0
0
-2
-2
-4
-4
-6 1
-6 1
0.5
0.5
8 6
0
-4
x 10
x 10
4
-0.5 -1
2 0
8 6
0
-4
4
-0.5 -1
F=264,0566 Hz (BW)
2 0
F=296,0223 Hz (FW)
Gambar 8 Enam modeshape pertama pada kecepatan 14000 rpm Defleksi Apabila di-klik tombol “Defleksi rotor” maka akan tampil defleksi rotor yang terjadi pada poros saat melewati putaran kritis 6949,35 rpm akibat eksitasi massa takseimbang, seperti tampak pada gambar 9. Defleksi Rotor
Diagram Stabilitas
-3 x 10
-4
-4
-5
2 0 -2 -4 4 8
2 x 10
-4
Komponen real nilai eigen
4
-6 -7 -8 -9 -10
6
0
4 -2
2 -4
N=6949,35 rpm
0
Gambar 9 Defleksi rotor sebesar 1,9e-4 m pada kecepatan kritis, N=6949,35 rpm
-11 -12 0
2000
4000 6000 8000 10000 Kecepatan putar, N(rpm)
12000
14000
Gambar 10 Komponen real solusi persamaan karakteristik sistem dualrotor
Stabilitas Sedangkan apabila kita meng-klik tombol “Stabilitas” maka akan ditampilkan diagram stabilitas sistem poros rotor (gambar 10) Dari gambar tersebut terlihat bahwa bagian real dari akar-akar persamaan karakteristik (nilai eigen) sistem dualrotor tersebut untuk kecepatan putar sampai 14000 rpm masih
12
Pemaparan Hasil Litbang 2003
Pusat Penelitian Informatika - LIPI
berharga negatif. negatif. Oleh karena itu berdasar teori di atas sistem poros rotor adalah stabil.
KESIMPULAN Berdasarkan pemaparan pada bab-bab sebelumnya , maka dapat diambil beberapa kesimpulan antara lain: 1. Dengan program simulasi ini prediksi perilaku dinamik sistem poros rotor akan dapat diketahui dengan segera dan dengan biaya yang murah. 2. Dengan program simulasi ini kita dapat mengubah harga kekakuan sistem untuk mempengaruhi perilaku dinamik sistem. 3. Dengan program simulasi ini kita dapat mengubah harga redaman sistem untuk mempengaruhi perilaku dinamik sistem.
DAFTAR PUSTAKA Childs, Dara. (1993), Turbomachinery Rotordynamics, John Wiley and Sons Ltd, Totonto, Canada. Vance, John, M. (1987),” Rotordynamics of Turbomachinery”, John Wiley and Sons Ltd, Toronto, Canada. Lalanne, M., Ferraris, G. (1990), Rotor dynamics Prediction in Engineering, John Wiley and Sons Ltd, West Sussex, England. Widodo, Achmad. (2002),” Prediksi secara numerik perilaku dinamik sistem poros-rotor di industri”, Departemen Teknik Mesin, Pascasarjana, ITB. Lalanne, M., Berthier, P., Der Hagopian, J. (1983),” Mechanical Vibration for Engineers”, John Wiley and Sons Ltd. Thomson, W.T. (1981),”Theory of Vibration with Applications”, 2nd edition, Prentice Hall.
Kedeputian Ilmu Pengetahuan Teknik
13
Bandung, 29 – 30 Juli 2003
Lalanne, M., Ferraris, G., Maissonneuve, V. (1996),” Prediction of the Dynamics behavior of Non-Symmetric Coaxial Co-or Counter Rotating Rotors”, Journal of Sound and Vibration, 195(4), pp649-666. Lalanne, M., Ferraris, G. (1999),” Use of the Campbell Diagram in Rotordyanamics,” Proceeding of DETC’99 ASME, Nevada.
14
Pemaparan Hasil Litbang 2003