1. PENDAHULUAN

Dalam teknologi perawatan prediktif, pengamatan dan analisis data getaran pada suatu sistem poros rotor merupakan hal yang penting. Pengukuran getaran pada suatu sistem poros rotor umumnya dilakukan dengan memasang akselerometer pada rumah bantalan. Dalam transmisinya menuju ke rumah bantalan, gaya yang dibangkitkan oleh elemen-elemen rotasi akan ditransmisikan melalui elemen rotasi pada bantalan. Oleh karena itu, agar analisis spektrum getaran suatu sistem poros rotor dapat dilakukan dengan baik maka informasi mengenai geometri dan karakteristik bantalan perlu diketahui dengan baik.

Salah satu informasi yang diperlukan untuk menganalisis data getaran adalah hubungan antara gaya radial dengan defleksi yang terjadi pada bantalan. Penentuan hubungan ini dapat dilakukan dengan menggunakan metode sederhana yang tidak memperhitungkan adanya kelonggaran pada bantalan. Pada metode sederhana ini, defleksi yang terjadi dapat langsung ditentukan berdasarkan gaya yang diberikan. Jadi dalam metode ini bila gaya yang diberikan ke bantalan sama dengan nol, maka defleksi yang terjadi juga sama dengan nol.

Disamping metode tersebut, penentuan hubungan antara gaya dan defleksi dapat pula dilakukan dengan memperhitungkan adanya kelonggaran pada bantalan. Dalam metode ini, besarnya defleksi yang terjadi ditentukan dengan cara iterasi karena defleksi maksimum yang terjadi dan besarnya sudut daerah beban (load zone) sangat dipengaruhi oleh besarnya gaya yang diberikan. Dalam proses iterasi tersebut diperlukan data mengenai kekakuan kontak dari elemen rotasi di daerah beban. Kekakuan kontak ini ditentukan dengan memanfaatkan teori Hertzian [1]. Metode yang dipakai untuk menentukan hubungan antara gaya dan defleksi yang dikemukakan di atas telah dilakukan oleh Harris [2], dan White [3]. Perbedaan antara metode Harris dan metode White terletak pada cara iterasinya. Namun karena hasil yang diperoleh dari kedua metode ini tidak jauh berbeda, maka proses analisis yang digunakan dalam penelitian ini akan didasarkan pada hasil yang diperoleh dari metode Harris. Dengan mengetahui hubungan antara gaya dan defleksi serta dengan memperhatikan faktor deformasi elemen tunggal pada titik kontak bantalan, kekakuan bantalan untuk berbagai gaya yang diberikan dapat ditentukan. Setelah analisis teoritik dengan menggunakan metode Harris selesai dilakukan dalam penelitian ini, selanjutnya dilakukan kaji eksperimental. Hasil eksperimental tersebut dimanfaatkan untuk validasi proses pemodelan yang dilakukan oleh Harris [2] dan White [3]. Selanjutnya, data hasil pengujian dan pemodelan tersebut dibandingkan untuk mengetahui metode pemodelan yang paling sesuai dengan data pengujian.

2. DISTRIBUSI BEBAN PADA ELEMEN ROTASI BANTALAN ROL

2.1 Pembebanan Bantalan dalam Arah Radial

Kekakuan statik dari suatu bantalan didefinisikan sebagai laju perubahan gaya terhadap defleksi. Secara matematik hubungan ini dapat dinyatakan sebagai berikut

\[K = \frac{dF}{d\delta} \tag{1}\]

dimana K menyatakan kekakuan bantalan, F menyatakan gaya pada bantalan, dan δ menyatakan defleksi yang terjadi.

Dalam perumusan di atas, kekakuan statik dari bantalan akan sangat dipengaruhi oleh distribusi gaya yang terjadi di antara dua elemen rotasi dan hubungan antara beban maksimum yang dapat diterima oleh masingmasing elemen rotasi dengan beban yang diberikan. Untuk jenis bantalan rol yang memiliki kontak garis, metode sederhana menyatakan bahwa hubungan antara gaya dan defleksi yang adalah sebagai berikut [2]

\[F = K_r \delta^{1,1} \tag{2}\]

Dalam analisis di atas perlu diketahui bagaimana sebenarnya distribusi beban yang terjadi pada seluruh elemen rol dari suatu bantalan serta defleksi yang terjadi pada masing-masing elemen. Besar defleksi radial dari setiap elemen rotasi pada suatu bantalan yang ditumpu pada rumah bantalan yang kaku dan mengalami beban dalam arah radial dinyatakan dengan persamaan (3) berikut

\[\delta_{\psi} = \delta_{r} \cos \psi - \frac{1}{2} C_{d} \tag{3}\]

dimana \(\delta_r\) menyatakan besar defleksi arah radial pada sudut \(\psi\)=0 dan \(C_d\) menyatakan besarnya kelonggaran.

Gambar 1 berikut ini memberikan ilustrasi mengenai notasi, dimensi dan kondisi defleksi pada bantalan rol serta hubungannya dengan kelonggaran yang dimilikinya

Gambar 1 Ilustrasi defleksi pada bantalan rol (a) Sebelum terjadi defleksi

(c)

(b) Setelah ada defleksi arah radial

(c) Bentuk distribusi beban

Selanjutnya, persamaan (3) dapat dinyatakan sebagai fungsi dari defleksi maksimum yaitu sebagai berikut:

\[\delta_{\psi} = \delta_{\text{max}} \left[ 1 - \frac{1}{2\epsilon} (1 - \cos \psi) \right]\] (4)

dimana

\[\varepsilon = \frac{1}{2} \left( 1 - \frac{C_d}{2\delta_r} \right) \tag{5}\]

Dari persamaan (4) dan (5) dapat diketahui bahwa sudut batas daerah beban sangat dipengaruhi oleh kelonggaran dari bantalan dan besarnya sudut ini dapat dinyatakan dalam persamaan berikut

\[\psi_1 = \cos^{-1}(\frac{C_d}{2\delta_r}) \tag{6}\]

Jadi untuk kelonggaran sama dengan nol maka ψ₁=90°.

Pengaruh kelonggaran dari suatu bantalan terhadap sudut batas dari daerah beban dapat ditunjukkan dengan jelas dalam Gambar 2 berikut.

Gambar 2 Pengaruh kelonggaran terhadap sudut batas daerah beban

(a) \(C_d=0\), \(\epsilon=0,5\), \(\psi_i=\pm90^\circ\)

(b) \(C_d > 0\), \(0 < \epsilon < 0.5\), \(0 < \psi_l < 90^\circ\), contoh: clearance

(c) \(C_d\)<0, 0,5<\(\epsilon\)<1,0, 90°<\(\psi_l\)<180°,contoh: beban awal

Dari gambar di atas tampak bahwa besarnya beban yang terjadi pada posisi w mengikuti hubungan berikut:

\[\frac{Q_{\psi}}{Q_{\text{max}}} = \left(\frac{\delta_{\psi}}{\delta_{\text{max}}}\right)^{n} \tag{7}\]

Dalam persamaan di atas \(Q_{\psi}\) menyatakan beban yang terjadi pada posisi sudut \(\psi\) sedangkan Qmax menyatakan beban maksimum pada posisi \(\psi\) =0.

Selanjutnya, dengan memakai persamaan (4), maka persamaan (7) dapat dinyatakan sebagai berikut:

\[Q_{\psi} = Q_{\text{max}} \left[ 1 - \frac{1}{2\varepsilon} (1 - \cos \psi) \right]^{n}\] (8)

Sesuai dengan kondisi keseimbangan statik, maka besar gaya radial (F_r) yang bekerja sama dengan jumlah gaya dalam arah vertikal, sehingga:

\[F_{r} = \sum_{\psi=0}^{\psi=\pm\psi_{I}} Q_{\psi} \cos \psi \tag{9}\] atau

\[F_{r} = Q_{\text{max}} \sum_{\psi=0}^{\psi=\pm\psi_{1}} \left[ 1 - \frac{1}{2\varepsilon} (1 - \cos\psi) \right]^{n} \cos\psi \qquad (10)\]

Selanjutnya, persamaan (10) dapat dinyatakan dalam bentuk integrasi sebagai berikut

\[F_{r} = ZQ_{max}x \frac{1}{2\pi} \int_{-\psi_{l}}^{+\psi_{l}} \left[1 - \frac{1}{2\varepsilon}(1 - \cos\psi)\right]^{n} \cos\psi \ d\psi \ (11)\] atau

\[F_{r} = ZQ_{max}J_{r}(\varepsilon) \tag{12}\]

dimana

\[J_{r}(\varepsilon) = \frac{1}{2\pi} \int_{-\psi_{1}}^{+\psi_{1}} \left[ 1 - \frac{1}{2\varepsilon} (1 - \cos \psi) \right]^{n} \cos \psi \ d\psi \quad (13)\]

Dalam persamaan (11) dan (12), Z menyatakan jumlah rol, sedangkan \(J_r(\varepsilon)\) menyatakan parameter sebagai fungsi dari jenis bantalan, sudut daerah beban dan kelonggaran dari bantalan.

Padahal, dari persamaan (2) dapat diperoleh hubungan berikut:

\[Q_{\text{max}} = K_n \delta_{\psi=0}^n = K_n (\delta_r - \frac{1}{2} C_d)^n\] (14)

sehingga dengan mensubstitusikan persamaan ini ke dalam persamaan (12) dapat diperoleh hubungan:

\[F_{r} = ZK_{n} \left(\delta_{r} - \frac{1}{2}C_{d}\right)^{n} J_{r}(\varepsilon)\] (15)

Apabila data teknis dan kelonggaran dari bantalan diketahui maka persamaan (15) dapat diselesaikan dengan metode iterasi. Dalam hal ini, harga ε mulamula diasumsikan, kemudian ψ₁ dihitung dengan menggunakan persamaan (5) dan (6). Setelah itu persamaan (13) digunakan untuk menentukan J. Jika data-data hasil perhitungan tersebut dimasukkan ke dalam persamaan (11) untuk gaya F, tertentu, maka akan diperoleh data \(\delta_r\). Apabila data \(\delta_r\) ini selanjutnya dimasukkan kembali kedalam persamaan (5), maka akan diperoleh data ε yang baru. Apabila harga ε yang diasumsikan mula-mula tidak sama dengan harga ε hasil perhitungan dari persamaan (5) maka proses iterasi ini harus diulangi. Proses iterasi ini harus dilakukan sampai harga ε yang diasumsikan sama dengan yang diperoleh dari persamaan (5) atau dicapai perbedaan harga sekitar ≈1x10^-6. Sebagai ilustrasi dari proses iterasi itu, dalam Gambar 3 berikut ini ditampilkan grafik antara Jr(ε) dengan ε untuk bantalan bola yang memiliki kontak titik dan bantalan rol yang memiliki kontak garis.

Gambar 3 Grafik antara Jr(ε) dengan ε untuk bantalan bola dan bantalan rol

Agar kedua kurva dalam Gambar 3 menjadi lebih jelas terutama untuk daerah dengan faktor beban \(\varepsilon\) <1, maka dalam Gambar 4 berikut ini kedua kurva tersebut ditampilkan dalam skala semi-log.

Gambar 4 Grafik antara Jr(ε) dengan ε dalam skala semi-log

Data \(J_r(\epsilon)\) diperoleh dengan melakukan proses integrasi dari \(-\psi_1\) sampai dengan \(\psi_1\). Jumlah elemen rotasi yang berada pada daerah integrasi tersebut (daerah beban) akan dijelaskan dalam uraian berikut.

2.2 Distribusi Gaya pada Elemen Rotasi di Daerah Beban

Defleksi sudut yang dialami oleh setiap elemen rotasi ditentukan oleh jumlah elemen rotasi yang berada di daerah beban. Besar daerah beban ini dibatasi oleh sudut \(\psi_l\). Apabila sudut antara masing-masing elemen

rotasi ditentukan dengan menggunakan persamaan (16) berikut

\[\theta = \frac{2\pi}{Z} \tag{16}\]

maka proses integrasi persamaan (15) harus dilakukan untuk seluruh sudut \(\psi_i\)=i\(\theta\) yang terletak di antara \(\pm \psi_i\), dimana i adalah bilangan asli =1,2,3,... . Selanjutnya, jika beban stedi dikenakan pada elemen rotasi dari bantalan, maka besar sudut yang membatasi daerah beban akan tetap konstan. Tetapi apabila elemen rotasi yang berada di daerah beban diputar sebesar sudut \(\phi\) dari kondisi awal \(\phi\)=0, maka akan terjadi perubahan jumlah elemen rotasi yang ada di daerah beban, seperti ditunjukkan dalam Gambar 5 berikut.

Gambar 5 Posisi angular elemen rotasi melalui daerah beban

• (a) Elemen rotasi memasuki daerah beban
• (b) Elemen rotasi berputar sebesar φ°.

Pada Gambar 5 di atas terlihat bahwa pada saat elemen rotasi tersebut berputar sebesar \(\phi\), jarak angular antara masing-masing elemen rotasi tidak berubah karena dibatasi oleh cage. Oleh karena itu, posisi angular dari setiap elemen rotasi yang berada pada sudut putar sebesar \(\phi\) antara 0 dan \(\theta\) dapat dinyatakan dalam persamaan berikut:

\[\psi_i = \phi + i\theta \tag{17}\] di mana \(\psi_i\) terletak antara -\(\psi_1\) dan \(\psi_i\)

Jadi, dapat disimpulkan bahwa penyelesaian persamaan (15) dapat dilakukan untuk berbagai posisi elemen rotasi dan hasil yang diperoleh akan tetap sama.

Dengan diperolenya data \(J_r(\epsilon)\) dan \(\delta_r\) hasil iterasi, serta data Z dan \(C_d\) dari spesifikasi teknik, maka langkah selanjutnya adalah menghitung kekakuan kontak, \(K_n\) pada elemen rotasi bantalan. Data ini diperlukan untuk menentukan gaya radial yang bekerja pada masingmasing elemen rotasi. Persamaan yang digunakan untuk menghitung kekakuan kontak dari elemen rotasi bantalan akan dijelaskan dalam uraian berikut.

2.3 Kekakuan Kontak Elemen Rotasi Tunggal

Penentuan kekakuan kontak dapat dilakukan dengan memanfaatkan teori kontak Hertzian pada kontak dua permukaan bola. Untuk bantalan bola, defleksi lintasan dalam dan lintasan luar ditentukan dengan persamaan berikut

\[\delta_i = 1,55 \left( \frac{F^2}{E^2} \frac{d_i + D}{d_i D} \right)^{1/3}\] (18)

\[\delta_o = 1.55 \left( \frac{F^2}{E^2} \frac{d_o - D}{d_o D} \right)^{1/3} \tag{19}\]

dimana F: gaya yang diberikan (N)

E:modulus elastisitas bahan rol (Nm^-2)

Persamaan (18) dan (19) dapat dituliskan dalam bentuk persamaan (2) dimana kekakuan kontak untuk lintasan dalam dan lintasan luar dihitung dari hubungan berikut

\[K = \frac{F}{\delta^{1,5}} \tag{20}\]

Sedangkan untuk bantalan rol, kekakuan kontak permukaannya dinyatakan dengan persamaan

\[K_n = 7.86 \times 10^4 L^{\frac{8}{3}}\] (21)
dimana L menyatakan panjang rol (mm)

Menurut White [3], kekakuan total dari suatu bantalan dapat dihitung dengan memasukkan faktor elemen geometri dan kelonggaran dari bantalan tersebut dan dinyatakan dalam persamaan berikut

\[\overline{K} = \frac{F}{x_m - \frac{1}{2}C_d} = \frac{F}{\delta}\] (22)

Hasil perhitungan dari persamaan (22) menghasilkan hubungan antara gaya yang bekerja dengan defleksi yang terjadi dan dinyatakan dalam persamaan berikut.

Bantalan bola \[\delta_b = 4.36 \times 10^{-8} \frac{F_m^{\frac{1}{13}}}{D^{\frac{1}{13}}}\] (23)

Bantalan rol \[\delta_1 = 3.06 \times 10^{-10} \frac{F_m^{0.9}}{L^{0.8}}\] (24)

3. KEKAKUAN TEORITIK DARI BANTALAN ROL

Data geometri dari bantalan rol yang akan dianalisis kekakuannya ditunjukkan dalam Tabel 1 berikut.

Tabel 1 Data geometri bantalan rol

Hubungan antara gaya dan defleksi yang dihasilkan dari perhitungan berdasarkan teori sederhana, teori Harris dan teori White ditunjukkan dalam Gambar 6 berikut.

Gambar 6 Hubungan antara gaya dengan defleksi pada bantalan rol

Pada Gambar 6 terlihat bahwa dengan memasukkan faktor kelonggaran dari bantalan pada teori White dan teori Harris, telah berhasil diperoleh data yang cukup Pada daerah dimana kondisi ketaklinierannya baik. cukup kecil, dapat ditunjukkan bahwa teori White menghasilkan harga kekakuan yang lebih besar daripada teori Harris. Perbedaan ini karena estimasi kekakuan elemen yang dikembangkan oleh Harris sedikit berbeda dengan White.

Setelah hubungan antara gaya dengan defleksi berhasil diperoleh, selanjutnya kekakuan bantalan sebagai fungsi dari defleksi yang terjadi dihitung dari gradien kurva dalam Gambar 6. Hubungan antara kekakuan bantalan dengan defleksi dinyatakan dalam Gambar 7 berikut.

Gambar 7 Kekakuan bantalan rol sebagai fungsi defleksi

Pada Gambar 7 di atas terlihat bahwa kekakuan bantalan tidak linier terhadap defleksi yang terjadi. Dalam Gambar 8 berikut ini akan disajikan hubungan antara kekakuan bantalan terhadap defleksi yang telah dinormalkan terhadap diameter rol dari bantalan sehingga menjadi bilangan tak berdimensi.

Gambar 8 Data non dimensional kekakuan bantalan

Maksud dari penentuan data nondimensional kekakuan bantalan terhadap defleksi pada hasil perhitungan di atas adalah untuk menunjukkan bahwa data-data tersebut untuk bantalan rol berada pada pada kurva yang sama dengan bantalan bola.

Dengan mengasumsikan bahwa massa efektif dari bantalan dan poros rotor adalah m, maka variasi frekuensi pribadi terhadap beban dapat dihitung dengan menggunakan persamaan berikut.

\[\overline{\omega_n} = \sqrt{\frac{\overline{K}}{m}}\] (25)

Hasil dari perhitungan frekuensi pribadi sebagai fungsi dari beban ditunjukkan dalam Gambar 9 berikut.

Gambar 9 Variasi frekuensi pribadi sebagai fungsi besar beban

Pada Gambar 9 di atas terlihat bahwa frekuensi pribadi dari bantalan menjadi tinggi sekali dengan cepat akibat pemberian gaya yang cukup kecil. Frekuensi pribadi dari bantalan ini proposional terhadap harga \(\sqrt{K}\) karena massa dinamik bantalan yang digunakan adalah konstan terhadap variasi gaya yang bekerja. Frekuensi pribadi dari bantalan tersebut menjadi relatif konstan pada beban sekitar 4000 N, penambahan beban selanjutnya hanya menyebabkan perubahan kecil dari harga frekuensi pribadi tersebut. Besarnya slope kekakuan bantalan pada beban yang kecil dapat terjadi karena adanya ketaklinieran akibat kelonggaran pada bantalan. Disamping itu, penentuan massa, m, juga merupakan kesulitan tersendiri karena tidak mudah untuk menentukan harga massa dinamik dari poros yang ditahan oleh bantalan.

4. KAJI EKSPERIMENTAL

Setelah karakteristik teoritik dari bantalan rol diuraikan , maka pada bahasan berikut ini akan dijelaskan mengenai penentuan kekakuan bantalan rol secara eksperimental. Hasil yang diperoleh dari pengujian ini selanjutnya dibandingkan dengan harga kekakuan teoritik.

Susunan perangkat uji untuk menentukan kekakuan bantalan rol ditunjukkan dalam Gambar 10 berikut.

Gambar 10 Susunan perangkat uji untuk menentukan kekakuan bantalan rol

Dalam Gambar 10 di atas tampak bahwa penentuan kekakuan bantalan rol dilakukan dengan menggunakan aktuator hidrolik sebagai penghasil gaya dan tessa modul untuk mengukur defleksi yang terjadi. Modus pengoperasian aktuator hidrolik yang digunakan dalam pengujian ini adalah modus gaya karena yang ingin diatur adalah besar beban yang terjadi. Jenis gaya yang dihasilkan diatur dengan menggunakan pembangkit sinyal dan pengendali servo. Besarnya gaya yang diningkan dan gaya yang terukur oleh load cell ditampilkan pada monitor pengatur hidrolik.

Dalam pengujian ini jenis gaya yang diberikan adalah gaya statik. Gaya yang diberikan memiliki rentang antara 0 – 4000N dengan perubahan naik dan turun setiap 250N. Besaran gaya yang dicatat adalah gaya yang diinginkan dan gaya yang terukur oleh load cell. Sedangkan besarnya defleksi yang terjadi dapat diketahui dari data yang ditunjukkan oleh tessa modul. Data gaya dan defleksi yang diperoleh kemudian digambarkan dalam suatu grafik yang menyatakan hubungan antara gaya yang diberikan dengan defleksi yang terjadi. Data hasil pengujian bantalan rol disajikan dalam Gambar 11 berikut.

Gambar 11 Kurva antara gaya radial dan defleksi hasil pengujian pada bantalan rol

Dalam Gambar 11 di atas ditunjukkan hubungan antara gaya dengan defleksi hasil pengujian untuk berbagai posisi sudut dari bantalan. Dari grafik di atas dapat diketahui bahwa data hasil pengujian untuk berbagai posisi sudut ternyata menunjukkan hasil yang hampir sama. Dengan demikian, hubungan antara gaya dan defleksi dapat diwakili dengan data dari salah satu posisi sudut. Agar kaji banding antara data hasil pemodelan dan data hasil pengujian dapat dilakukan dengan lebih baik, maka kedua data tersebut ditampilkan secara bersama-sama seperti terlihat dalam Gambar 12 berikut.

Dari Gambar 12 dapat diketahui bahwa data hasil pengujian ternyata tersebar di antara ketiga metode perhitungan yang digunakan dalam pemodelan. Untuk menentukan metode perhitungan yang paling mendekati hasil pengujian, perlu dihitung perbedaan antara hasil pengujian dengan hasil pemodelan untuk ketiga jenis metode yang dipakai. Dalam Tabel 2 berikut ini akan ditunjukkan perbedaan antara hasil pengujian dengan data hasil perhitungan yang didasarkan pada teori White, sederhana dan Harris.

Tabel 2 Perbedaan hasil pengujian dengan pemodelan

Gambar 12 Kaji banding antara data gaya dengan defleksi hasil pemodelan dan pengujian

Pada baris terakhir dari Tabel 2 di atas dinyatakan harga rata-rata dari perbandingan antara defleksi eksperimen dengan defleksi teoritik untuk ketiga teori. Apabila perbedaan antara harga perbandingan tersebut dengan harga l semakin kecil, berarti bahwa teori yang dipakai dalam pemodelan semakin mendekati data hasil pengujian.

Berdasarkan data dalam Tabel 2 di atas dapat diketahui bahwa teori White merupakan model teoritik terbaik yang memiliki perbedaan dengan harga 1 sebesar 0,002. Kemudian diikuti oleh teori Harris yang memiliki perbedaan sebesar 0,024 dan yang terakhir adalah teori sederhana dengan harga perbedaan sebesar 0,044.

5. KESIMPULAN

Dari analisis teoritik dan kaji eksperimental yang telah dilakukan dalam makalah ini, dapat ditarik beberapa kesimpulan berikut:

Efek ketaklinieran karena adanya bantalan dapat ditunjukkan dengan baik dengan menggunakan grafik antara gaya radial dengan defleksi yang terjadi.

2. Model teoritik yang telah dikembangkan untuk memperkirakan kekakuan bantalan ternyata sudah cukup baik dalam mendekati hasil pengujian.

Diantara ketiga model teoritik yang dikembangkan untuk menentukan kekakuan bantalan rol, ternyata teori White menghasilkan kesalahan terkecil terhadap data pengujian. Kemudian disusul oleh teori Harris dan yang terakhir adalah teori sederhana.

Dengan memperhitungkan faktor kelonggaran pada perhitungan kekakuan bantalan ternyata diperoleh hasil perhitungan defleksi yang lebih baik.

6. DAFTAR PUSTAKA

• Hertz, H.,: "On the Contact of Rigid Elastic Solids and on Hardness", MacMillan, U.K., 1896.
• Harris, T. A.,: "Rolling Bearing Analysis", 3^rd Ed., John Wiley & Sons, Inc., USA, 1991.
• White, M. F.,:"Rolling Element Bearing Vibratrion Transfer Characteristics: Effect of Stiffness", Journal of Applied Mechanics, Vol. 46, U.K., 1979.
• 4. ----,:"MATLAB User's Guide & Reference Guide for Microsoft Windows", The Math Works Inc., USA, 1994.