SIMULASI DAN PERANCANGAN SISTEM KONTROL SUSPENSI SEMI AKTIF MODEL SEPEREMPAT KENDARAAN

SIMULASI DAN PERANCANGAN SISTEM KONTROL SUSPENSI SEMI AKTIF MODEL SEPEREMPAT KENDARAAN

TUGAS AKHIR

Oleh: Aloysius Dodi Setyobudi 133 94 080

JURUSAN TEKNIK FISIKA FAKULTAS TEKNOLOGI INDUSTRI INSTITUT TEKNOLOGI BANDUNG 2001



Tugas Akhir untuk melengkapi syarat sebagai Sarjana Teknik Fisika Institut Teknologi Bandung

Dibimbing oleh: Dipl.-Ing.Ir. Nyoman Bangsing Dr. Ir. IGN. Wiratmaja Puja




Disetujui oleh : Pembimbing I

Pembimbing II

Dipl.-Ing.Ir. Nyoman Bangsing

Dr. Ir. IGN. Wiratmaja Puja

NIP: 131 41 7437

NIP: 131 83 5240


Kupersembahkan kepada Mama-Papa dan Adik-adik tercinta, Serta Robertha Eka Woro A., AMK.

iii

ABSTRAK Tugas akhir ini menawarkan satu alternatif untuk sistem suspensi aktif, dimana sistem suspensi yang akan diteliti menggunakan peredam kejut dengan fluida yang kekentalannya dapat dikontrol, yang dikenal dengan fluida elektrorheologi (FER). FER bersifat reversible, yaitu dapat berubah dari mengalir bebas sampai menjadi semi-padat bila padanya diberikan pengaruh medan listrik. Teknologi FER ini memiliki keunggulan dibandingkan dengan teknologi konvensional, diantaranya daya beban listrik maksimal yang diperlukan hanya kurang dari 4 Watt. Di samping itu konstruksi perangkatnya relatif sederhana sehingga dapat menekan biaya produksi. Dalam penelitian ini dilakukan pemodelan terhadap peredam kejut FER (jenis TXER-6) menggunakan pendekatan geometri pelat paralel dan cylindricalaxisymmetric 1-dimensi Newton. Dari jenis FER yang dipilih, dengan jarak antar elektroda 0,5 mm dan panjang celah 63,5 mm, bila dikenakan medan listrik sebesar 4 kV/mm akan dihasilkan gaya redam maksimal sebesar 546,37 N. Pada saat tidak ada medan listrik, perangkat ini akan menjadi seperti peredam kejut konvensional biasa. Pada simulasi yang dilakukan, sebagai masukan sistem dipilih gangguan jalan bump. Metoda yang dipilih dalam perancangan pengontrolnya adalah metoda kontrol optimal. Pengujian sistem suspensi pasif menggunakan masukan jalan bump dengan tinggi 0.04 m dan lebar bump 1.2 m ketika kecepatan kendaraan konstan 20 km/jam selama 1.5 detik menghasilkan percepatan vertikal (rms) badan kendaraan sebesar 3,2527 m/s2, defleksi suspensi (rms) 0,01334 m, dan defleksi ban sebesar 0,00698 m (rms). Pada kontrol optimal dengan umpan balik keadaan, sistem aktif dan FER, masing-masing, akan menghasilkan perbaikan pada percepatan (rms) badan kendaraan sebesar 74,34% (menjadi 0,8348 m/s2) dan 51,89% (1,5648 m/s2), defleksi (rms) ruang suspensi sebesar 5,17% dan 2,19%, dan pada defleksi (rms) ban sebesar 66,39% dan 45,70%. Pada kontrol optimal yang menggunakan indeks performansi yang melibatkan pengoptimalan perintah gaya kontrol, sistem aktif dan FER memberikan hasil perbaikan pada percepatan vertikal (rms) badan kendaraan sebesar 27,69% (2,3520 m/s2), defleksi ruang suspensi sebesar 6,35% dan defleksi ban sebesar 26,59%.

iv

KATA PENGANTAR Segala puji syukur penulis panjatkan kepada Allah Yang Maha Murah karena atas Kuasa-Nya dan Kemurahan HatiNya-lah penulis akhirnya sanggup menuntaskan laporan tugas akhir yang berjudul “Simulasi dan Perancangan Sistem kontrol Suspensi Semi Aktif Model Seperempat Kendaraan” untuk melengkapi syarat sebagai Sarjana Teknik Fisika – ITB. Penulisan laporan tugas akhir ini dapat tuntas karena peran serta besar dari mereka yang berada di lingkungan penulis mulai dari tahun awal pengerjaan hingga saat ini. Atas dorongan, dukungan, dan peranan mereka yang demikian besar bagi penulis, pada pengantar laporan ini penulis hendak menyampaikan penghargaan dan terima kasih sebesar-besarnya kepada: 1. Bapak Dipl.-Ing.Ir. Nyoman Bangsing, selaku dosen wali dan dosen pembimbing utama, atas kesabaran dan kebesaran hatinya tetap mendukung penulis untuk menuntaskan tugas akhir ini. 2. Bapak Dr. Ir. IGN Wiratmaja Puja, selaku dosen pembimbing kedua, yang telah banyak memberikan inspirasi awal pada tugas akhir ini, juga atas kesabaran dan kebesaran hatinya untuk tetap mendukung penulis hingga menuntaskan studinya. 3. Bapak Dr.-Ing. Parsaulian Siregar dan Bapak Dr.Ir. Nugraha, atas kesediaannya untuk menjadi anggota tim penguji tugas akhir ini. 4. Bapak Dr.-Ing.Ir. Yul Y. Nazaruddin, MSc., DIC., yang telah bersedia menerima kehadiran Bapak dan Ibu saya kala itu. 5. Segenap Bapak dan Ibu Dosen Teknik Fisika yang pernah menuangkan ilmu yang dimiliki kepada penulis selama menempuh pendidikan. 6. Segenap staf TU dan Perpustakaan Teknik Fisika yang bersedia dengan rela membantu kelancaran administrasi selama penulis membutuhkannya. 7. Teman-teman Teknik Fisika Angkatan 1994 yang telah berpencar menempuh perjuangan atas dukungan moralnya. 8. Mbak Anjar dan seluruh penghuni GEMA, atas dukungan doanya. 9. Teman-teman penghuni Burung Gereja 2, maupun alumninya, karena tetap rajin mengingatkan penulis untuk berjuang menuntaskan laporan ini dan bersedia menemani siang dan malam. 10. Teman-teman penghuni mailing-list [email protected]

v 11. Bapak Ir. Degus Rustianto, Ir. Bambang Sunarjo, Pak Purwadiyono dan Mbak Useu Sopiah atas kesediaannya menerima penulis untuk bergabung di lingkungan kerja di PT Solar Services Indonesia dan mengembangkan kemampuan dan wawasan penulis. 12. Special thanks for Mr. Charles Lozinger, for your appreciation. 13. Mr. Ryuji Aizawa, Chief Researcher of Nippon Shokubai Co. Ltd., for your kindness to sent us your TX-ER6 material sample. Thank you very much for the future of the controllable fluids. 14. Dan terlebih kepada Mama dan Papa, yang dengan penuh kasih, bersedia tanpa pernah lelah mengingatkan penulis untuk tetap berjuang dan berjuang tanpa pernah putus asa, karena Tuhan Maha Murah bagi mereka yang selalu berjuang. 15. Sekali lagi dan selalu kepada Papa-Mama, Dik Rini dan Dik Niken, serta tambatan hatiku – Eka Woro, yang dengan penuh kasih dan cinta menemani, mendampingi, membantu penuh pengorbanan demi terselesaikannya studi penulis di jenjang S1 ini.

Sepatah dua patah kata dalam kata pengantar ini tidak akan cukup untuk menggambarkan rasa syukur dan terima kasih penulis kepada semua pihak yang telah berperan dan ada di dalam kehidupan penulis. Juga kemampuan saya mengingat satu persatu pihak-pihak yang terlibat tidak akan dapat menghilangkan penghargaan saya kepada mereka yang tidak sempat tersebutkan pada pengantar yang singkat ini.

Akhirnya semoga semua yang telah penulis alami selama proses pengerjaan hingga selesainya pengerjaan laporan tugas akhir ini sedikitnya akan bermanfaat bagi para pembaca dan civitas akademika ITB, khususnya bagi Jurusan Teknik Fisika.

Bandung, 22 September 2001

Aloysius Dodi Setyobudi

vi

DAFTAR ISI Judul ...............................................................................................................................i Lembar Pengesahan .........................................................................................................ii Abstrak ............................................................................................................................iii Kata Pengantar ................................................................................................................iv Daftar Isi ..........................................................................................................................vi Daftar Gambar .................................................................................................................ix Daftar Tabel .....................................................................................................................xi

BAB I PENDAHULUAN 1.1

Latar Belakang ......................................................................................................1

1.2

Tujuan ...................................................................................................................3

1.3

Metodologi Pengerjaan Tugas Akhir ....................................................................4

1.4

Batasan Masalah ....................................................................................................4

1.5

Sistematika Pembahasan .......................................................................................5

BAB II DASAR TEORI 2.1

Konsep Mekanika Fluida .......................................................................................6

2.1.1 Definisi Fluida .......................................................................................................6 2.1.2 Fluida Newton .......................................................................................................7 2.1.3 Fluida Non-Newton ...............................................................................................8 2.1.4 Aliran Internal Laminar pada Pelat Parallel Tak Terbatas .....................................10 2.2

Fluida Elektrorheologi (FER) ...............................................................................13

2.2.1 Sejarah Perkembangan FER ..................................................................................13 2.2.2 Mekanisme Kerja FER ..........................................................................................13 2.2.3 Karakteristik FER ..................................................................................................15 2.3

Getaran Mekanis ....................................................................................................16

2.4

Sistem Suspensi Seperempat Kendaraan................................................................18

2.4.1 Model Ruang Keadaan Sistem Suspensi Aktif Seperempat Kendaraan ...............19 2.5

Optimasi Sistem ....................................................................................................20

2.6

Sistem Kontrol Optimal Quadratic .......................................................................21

2.6.1 Integral Performansi Secara Umum ......................................................................24

vii

BAB III PEMODELAN SISTEM 3.1

Karakteristik bahan Fluida Elektrorheologi ..........................................................26

3.2

Sketsa Model Peredam Kejut Elektrorheologi ......................................................27

3.3

Pemodelan Peredam Kejut Menggunakan Pendekatan Geometri 1D Pelat Paralel .....................................................................................29

3.3.1 FER sebagai Fluida Newton ..................................................................................30 3.3.2 FER sebagai Plastik Bingham ...............................................................................31 3.3.3 Parameter-Parameter Model Pelat Paralel .............................................................37 3.4

Pemodelan Peredam Kejut Menggunakan Pendekatan Geometri 1D Cylindrical Axisymmetric ...............................................................40

3.4.1 Aliran Geser Newton .............................................................................................41 3.4.2 Aliran Fluida Bingham ..........................................................................................43 3.4.3 Parameter-Parameter Peredam Kejut FER ............................................................48 3.5

Model Suspensi Semi-Aktif Seperempat Kendaraan ............................................51

BAB IV PERANCANGAN SISTEM KONTROL OPTIMAL 4.1

Indeks Performansi yang Digunakan .....................................................................55

4.2

Sistem Kontrol Optimal pada Suspensi Semi Aktif Seperempat Kendaraan dengan Peredam Kejut Variabel ......................................56

4.2.1 Perancangan Pengontrol Optimal ..........................................................................57 4.2.2 Sistem Pengontrol Optimal ...................................................................................59

BAB V SIMULASI MODEL FER DAN SISTEM SUSPENSI SEMI AKTIF MODEL SEPEREMPAT KENDARAAN 5.1

Simulasi Model Peredam Kejut FER ....................................................................61

5.2

Pemilihan Parameter Model Seperempat Kendaraan ............................................66

5.3

Pemilihan Parameter Pembobot ............................................................................73

5.4

Respon Frekuensi Sistem Pasif dan Aktif Penuh ...................................................75

5.5

Respon Sistem terhadap Masukan Jalan Bump .....................................................77

BAB VI USAHA PERBAIKAN SISTEM KONTROL 6.1

Indeks Performansi Baru .......................................................................................82

6.2

Respon Sistem terhadap Masukan Jalan Bump .....................................................83

viii

BAB VII KESIMPULAN DAN SARAN 7.1

Kesimpulan ...........................................................................................................92

7.2

Saran untuk Pengembangan Lebih Lanjut ............................................................93

DAFTAR PUSTAKA LAMPIRAN A

LISTING PROGRAMS

A.1 M-file: txer6.m ......................................................................................................A1 A.2 M-file: erdamp.m ..................................................................................................A4 A.3 M-file: ermodel.m .................................................................................................A5 A.4 M-file: ermodel1.m ...............................................................................................A7 A.5 M-file: ermodel2.m ...............................................................................................A7 A.6 M-file: ermodel3.m ...............................................................................................A8 A.7 M-file: cardata.m ...................................................................................................A9 A.8 M-file: oplores.m ...................................................................................................A9 A.9 M-file: ferlaw.m ................................................................................................. A11 A.10 M-file: carmodel.m ............................................................................................ A12 A.11 M-file: carplot.m ................................................................................................ A15 LAMPIRAN B

VARIASI PARAMETER PEREDAMAN

B.1 Perubahan Koefisien Redaman dengan Jarak Celah Konstan (= 2 mm)................B1 B.2 Perubahan Koefisien Redaman dengan Panjang Saluran Celah Konstan (= 63.5 mm) ............................................................................................B2

B.3 Listing Program Variasi Parameter Redaman .......................................................B3 LAMPIRAN C

SIMULINK MODEL

ix

DAFTAR GAMBAR No. Gambar

Deskripsi Gambar

Halaman

2.1

Perilaku (a) padatan dan (b) fluida, bila dikenakan gaya geser secara konstan

6

2.2

Deformasi elemen fluida

7

2.3

9

2.5

(a) Shear stress, τ, dan (b) viskositas nyata (apparent), η, sebagai fungsi dari laju deformasi untuk aliran satu dimensi dari berbagai jenis fluida non-Newton Volume kontrol untuk analisis aliran laminar antara pelat paralel tak terbatas stasioner Properti Stress Fluida ER

2.6

Mekanisme Efek ER

15

2.7

Model Suspensi Pasif Seperempat Kendaraan

18

2.8

Model Suspensi Aktif

19

2.9

Sistem kontrol optimal

22

3.1

Shear Stress terinduksi terhadap medan listrik pada laju geser 200/sec untuk TX-ER6 Gambar Konsep Peredam kejut yang berisi FER

26

3.3

Volume kontrol untuk analisis aliran antara pelat paralel tak terbatas stasioner

29

3.4

Profil kecepatan yang dapat terjadi pada pelat paralel

32

3.5

Volume kontrol untuk analisis aliran laminar penuh pada pipa

40

3.6

Profil kecepatan pada celah melingkar untuk bahan plastik Bingham dengan adanya gradien tekanan linier pada arah aksial (sumbu x) Model Bingham peredam kejut dengan FER

43

Model suspensi semi-aktif seperempat kendaraan dengan peredam kejut variabel Diagram blok model matematis sistem suspensi seperempat kendaran Model suspensi semi aktif seperempat kendaraan dengan peredam kejut model Bingham Profil ketebalan plug berkurang ketika gaya bertambah untuk medan listrik yang konstan: – – 1kV/mm; –. 2 kV/mm; -- 3 kV/mm; – 4 kV/mm.

51

2.4

3.2

3.7 3.8 3.9 3.10 5.1

5.2 5.3

Profil ketebalan plug terhadap perubahan medan listrik ketika diberikan gaya konstan: – – 100 N; –. 500 N; -- 1000 N; – 1500 N; .. 2000 N. Profil kecepatan pada celah elektroda dengan adanya gaya 500 N.

10 14

28

50

53 54 61

62 63

DAFTAR GAMBAR

x

5.4

Profil kecepatan pada celah elektroda dengan adanya gaya 600N.

63

5.5

Profil kecepatan pada celah elektroda dengan adanya gaya 1000 N.

64

5.6

Gaya peredam terhadap eksitasi sinusoida 5 Hz dengan Amplituda sebesar 2 cm Respon frekuensi model Raymon M.

65, 66

Respon frekuensi defleksi suspensi sistem dengan konfigurasi seperti pada tabel 3.2 Respon step percepatan badan kendaraan dengan ρ1 yang berbeda (1, 10, 100, 1.000,10.000, 1.000.000) Respon frekuensi sistem pasif (garis kontinu) dan sistem aktif (garis putus-putus) Masukan jalan “bump” dengan amplituda 0.1 m

70, 71

5.7 5.8 5.9 5.10 5.11 5.12 5.13 6.1 6.2 6.3 6.4 6.5 6.6

67

75 76, 77 78

Perbandingan tanggapan waktu sistem pasif (garis kontinu) dan sistem aktif penuh (garis putus-putus): Gaya dari Peredam FER (garis putus-putus) vs. Gaya Aktif (garis kontinu) Gaya dari peredam FER (garis kontinu) dan Gaya aktif (garis putus-putus) Perbandingan tanggapan waktu terhadap masukan jalan bump

79, 80

Gaya dari peredam FER (garis kontinu) dan gaya aktif umpan balik keadaan penuh Perbandingan tanggapan waktu sistem kontrol umpan balik keadan penuh terhadap masukan jalan bump Perbandingan tanggapan waktu sistem kontrol umpan balik keluaran terhadap masukan jalan bump Gaya dari peredam kejut FER dan gaya aktif umpan balik keluaran

86

81 83 84, 85

87, 88 89,90 90

DAFTAR GAMBAR

xi

DAFTAR TABEL

No. Tabel

Deskripsi Tabel

2.1

Parameter-parameter model sistem suspensi

19

3.1

Koefisien polinom Yield Stress sebagai fungsi dari Medan Listrik

27

5.1

Konfigurasi parameter model seperempat kendaraan

5.2

Parameter sistem suspensi seperempat kendaraan awal

72

5.3

“Step Response Data”untuk perubahan konstanta pembobot

74

5.4

Harga RMS dari tanggapan sistem terhadap masukan jalan bump

81

6.1

Harga RMS dari tanggapan sistem terhadap masukan jalan bump

85

6.2

Harga RMS dari tanggapan sistem dengan kecepatan kendaraan 20 km/jam terhadap masukan jalan bump dengan tinggi 4 cm dan lebar 1,2 m

91

Halaman

68, 69

DAFTAR TABEL

1

BAB I PENDAHULUAN 1.1 LATAR BELAKANG Profil permukaan jalan yang tidak rata akan menyebabkan kenyamanan berkendaraan

terganggu,

dimana

pengemudi/penumpang

kendaraan

akan

terguncang apabila kendaraan melewati jalan yang rusak dan berlubang di sanasini. Untuk mengurangi ataupun meredam pengaruh getaran/guncangan yang dirasakan pengguna kendaraan maka diperlukan suatu sistem suspensi kendaraan. Pada sistem suspensi konvensional atau suspensi pasif, komponen suspensi hanya terdiri dari pegas dan peredam kejut. Suspensi jenis ini belum dilengkapi kemampuan untuk selalu beradaptasi dengan perubahan profil permukaan jalan. Di lain pihak masyarakat pengguna kendaraan menuntut perbaikan dalam hal kenyamanan maupun keamanan dalam berkendaraan. Untuk memenuhi tuntutan tersebut, maka banyak pihak mencoba untuk mengembangkan sistem suspensi aktif. Sistem suspensi aktif tidak lain adalah sistem suspensi pasif, dimana padanya telah ditambahkan elemen aktif, yang dapat berdefleksi/bergerak sedemikian rupa, sehingga roda akan selalu mengikuti perubahan profil permukaan jalan dan meminimalkan getaran/guncangan yang dirasakan pengguna kendaraan dengan mempertahankan badan kendaraan pada ketinggian/level yang tetap. Sistem suspensi aktif biasanya terdiri dari aktuator linier (elemen aktif), sensor getaran dan pengontrol elemen aktif. Mengingat daya yang diperlukan dalam pengontrolan riil sistem suspensi aktif cukup besar maka elemen aktif yang banyak dipilih adalah aktuator hidrolik. Namun perlu diingat bahwa elemen ini harganya relatif mahal. Oleh karena itu, pada tugas akhir ini dicoba dibahas salah satu komponen alternatif pada sistem suspensi aktif, yaitu menggunakan peredam kejut dengan fluida yang kekentalannya dapat dikontrol, seperti fluida elektro-rheologi (FER). FER adalah fluida yang memiliki kemampuan dapat balik (reversible) untuk

2

berubah dari mengalir bebas menjadi fluida viskos linier sampai menjadi semipadat bila diberikan pengaruh medan listrik. Karakteristik bahan FER adalah kompleks dan nonlinier karena adalah fungsi dari medan listrik, beban terpasang, amplituda strain dan frekuensi eksitasi medan listrik. Pada awalnya, peredam kejut yang akan digunakan dalam simulasi adalah modifikasi dari peredam kejut konvensional jenis twin-tube atau floating-piston, dengan memberikan tegangan listrik yang berbeda pada kedua pelat tabung yang saling berhadapan agar resistansi aliran fluida dapat diatur. FER, yang diperoleh dari pihak manufaktur, diisikan ke dalam peredam kejut hasil modifikasi tadi. Peredam kejut ini akan terlebih dahulu diidentifikasi melalui model matematis dan disimulasikan bersama dengan sistem kontrolnya menggunakan perangkat lunak, kemudian barulah akan divalidasikan dengan mengambil data karakteristik melalui eksperimen. Peredam kejut berteknologi FER ini diharapkan memiliki keunggulankeunggulan yang berarti bila dibandingkan dengan penggunaan teknologi konvensional dalam merancang peredam kejut aktif, diantaranya: 1.

Efek elektrorheologi efisien karena hanya membutuhkan daya listrik yang rendah. Arus listrik yang diambil adalah kurang dari 1 mA dari aplikasi tegangan 4 kV sehingga daya beban listrik kurang dari 4 W.

2.

Konstruksi perangkat peredam kejutnya sederhana karena pengaruh sifat nonmekanik fluida elektro-rheologi. Hal ini potensial untuk mengurangi biaya produksi dan meningkatkan kemampuan perangkat peredam kejut tersebut.

3.

Efek elektrorheologi mampu merespons sinyal kontrol dengan cepat.

3

Penelitian lebih lanjut setelah tugas akhir ini dinilai sangat strategis, mengingat belum banyak instansi/perusahaan di dunia yang melakukan penelitian di bidang ini, dan diharapkan akan memberikan kontribusi pada pengembangan sistem suspensi aktif, terutama dalam industri otomotif di Indonesia. Untuk itu dibutuhkan usaha yang berkesinambungan setahap demi setahap untuk membangun pemahaman yang utuh dan menyeluruh tentang teknologi ini. Tugas akhir ini hanyalah salah satu langkah awal kecil yang diharapkan membawa khususnya Jurusan Teknik Fisika dan Institut Teknologi Bandung, pada umumnya, sebagai perintis dalam bidang pengontrolan sistem suspensi aktif FER di Indonesia.

1.2 TUJUAN Adapun hasil yang ingin dicapai pada tugas akhir ini adalah: 1.

Menurunkan model matematis suspensi semi aktif model seperempat kendaraan yang menggunakan peredam kejut dengan FER.

2.

Merancang pengontrol optimal pada simulasi terhadap model yang dibuat, yaitu dengan memperhatikan hal-hal berikut: •

Meningkatkan

faktor

kenyamanan,

hal

ini

dilakukan

dengan

meminimumkan percepatan vertikal kendaraan terhadap acuan pada saat kendaraan tidak bergerak. •

Meningkatkan faktor keamanan, yaitu dengan menjaga agar ban tetap mengikuti permukaan jalan. Hal ini dilakukan dengan meminimumkan defleksi ban, sehingga diharapkan ban akan tetap melekat pada permukaan jalan.

4

1.3 METODOLOGI PENGERJAAN TUGAS AKHIR Pada pengerjaan tugas akhir ini, seperti telah diungkapkan dalam sub bab latar belakang di atas, dilakukan tahapan-tahapan sebagai berikut: 1. Melakukan studi literatur tentang FER, untuk mendapatkan model matematika yang dapat merepresentasikan dinamika FER. Bertujuan untuk dapat memilih bahan FER yang sesuai sehingga mendapatkan parameter penting karakter bahan yang dapat diimplementasikan pada sistem kontrol. 2. Menurunkan model matematis peredam kejut yang menggunakan fluida elektrorheologi. 3. Merancang sistem pengontrol suspensi semi aktif model seperempat kendaraan berdasarkan teori sistem kontrol optimal. 4. Melakukan simulasi sistem hasil point 2 dan 3 di atas. 5. Mengevaluasi performansi sistem dan memperbaiki model peredam kejut beserta performansi sistem pengontrolnya.

1.4 BATASAN MASALAH Masalah yang dibahas dalam studi perancangan ini dibatasi pada hal-hal berikut : 1. Dalam tugas akhir ini akan digunakan model fluida Bingham untuk membantu memodelkan FER, seperti yang dianjurkan oleh Kamath, G.M., dkk.[2]. Pemodelan ini diperlukan agar dapat dirancang sistem kontrol dengan performansi yang diinginkan. 2. Karakter bahan elektrorheologi dianggap hanya merupakan fungsi dari salah satu peubah pada property bahan tersebut. 3. Sistem suspensi yang digunakan adalah model seperempat kendaraan. 4. Kursi penumpang diangap sebagai satu kesatuan dengan massa badan kendaraan (massa sprung). 5. Gangguan yang ditinjau hanya berasal dari gangguan akibat ketidakrataan permukaan jalan bump.

5

1.5 SISTEMATIKA PEMBAHASAN Pembahasan pada tugas akhir ini disusun dalam bab-bab sebagai berikut: I.

Pendahuluan, berisikan latar belakang masalah, tujuan penelitian, metodologi pengerjaan tugas akhir, batasan masalah, dan sistematika pembahasan.

II.

Dasar Teori, berisikan teori-teori yang mendasari tugas akhir ini, yaitu Elektro-rheologi, Teori Getaran Mekanis, dan Teori Kontrol.

III. Pemodelan Fluida Elektrorheologi dan Gerak Kendaraan, berisikan dasar

serta proses penurunan model suspensi yang menggunakan fluida elektrorheologi serta spesifikasi data teknis seperempat kendaraan yang akan disimulasikan dalam tugas akhir ini. IV. Perancangan Sistem, yang berisikan metoda perancangan sistem kontrol peredam kejut yang akan diterapkan pada sistem suspensi tersebut. V. Simulasi Hasil Perancangan dan Interpretasinya, menguraikan hasil simulasi sistem suspensi aktif maupun pasif terhadap model peredam kejut. VI. Usaha Perbaikan Sistem Kontrol, yang membahas dan memperbaiki hasil simulasi kontrol sistem suspensi. VII. Kesimpulan dan Saran, yang berisikan kesimpulan hasil penelitian dan saransaran penulis untuk pengembangan lebih lanjut dari tugas akhir ini.

6

BAB II DASAR TEORI 2.1 KONSEP MEKANIKA FLUIDA[2] Tugas Akhir ini tidak terlepas dari disiplin ilmu mekanika fluida. Beberapa pengertian dasar mekanika fluida digunakan sebagai landasan pengembangan dan pembahasan-pembahasan selanjutnya dalam pengerjaan Tugas Akhir ini.

2.1.1

DEFINISI FLUIDA

Fluida adalah substansi yang berubah bentuk (deform) secara kontinu bila padanya dikenakan tekanan geser (shear/tangential stress) sekecil apapun. Fluida meliputi fasa cair (liquid) dan gas. Perbedaan antara fluida dan padatan (solid) dapat jelas terlihat dari perilakunya (behavior). Padatan berdeformasi jika dikenakan gaya tekan permukaan tapi deformasinya tidak kontinu terhadap waktu. Pada Gambar 2.1 diperlihatkan perbedaan perilaku padatan (Gbr. 2.1.a) dan fluida (Gbr. 2.1.b) bila dikenakan gaya geser F secara konstan pada pelat bagian atas. Pada Gbr. 2.1.a diperlihatkan bahwa selama batas elastis bahan padatan tidak terlampaui, deformasi adalah sebanding dengan tekanan geser (shear stress) yang dialaminya, τ = F A , dengan A adalah luas permukaan kontak dengan pelat. Sedangkan pada Gbr. 2.1.b, ketika gaya F dikenakan pada pelat sebelah atas, elemen fluida terus berdeformasi selama dikenakan gaya. Fluida yang kontak langsung dengan pelat mempunyai kecepatan yang sama dengan pelat tersebut, dengan anggapan tidak ada slip pada perbatasan pelat dan fluida. Bentuk elemen fluida, untuk waktu berturutan t0
F

F t0 t 1

t2 t0 > t1 > t2

(a) Gambar 2.1

(b)

Perilaku (a) padatan dan (b) fluida, bila dikenakan gaya geser secara konstan

7

2.1.2

FLUIDA NEWTON

Tanpa ada tekanan geser (atau mulai dari saat ini penulis sebut langsung dengan shear stress), tidak akan terjadi deformasi fluida. Fluida dapat diklasifikasikan secara umum berdasarkan hubungan antara shear stress yang dikenakan dan laju δl M

M'

P

P'

Force, Fx Velocity, δu

δα

Fluid element y at time, t

Fluid element at time, t + δt

δy

x N

O δx

Gambar 2.2

Deformasi elemen fluida

deformasinya. Pada Gambar 2.2 ditunjukkan deformasi elemen fluida di antara dua pelat tak terbatas. Pelat atas bergerak pada kecepatan konstan, δu, karena pengaruh dikenakannya gaya F konstan seperti pada gambar. Shear stress yang dikenakan pada elemen fluida diberikan oleh persamaan berikut: δFx dFx = δA y → 0 δ A dAy y

τ yx = lim

(2.1)

dimana δAy adalah luas daerah kontak elemen fluida dengan pelat, dan δFx adalah gaya yang diteruskan pelat kepada elemen fluida. τyx menyatakan shear stress pada bidang y dengan arah x. Pada selang waktu δt, elemen fluida terdeformasi dari kedudukan MNOP ke M’NOP’. Laju deformasinya adalah: δα dα = δt →0 δt dt

Laju deformasi = lim

(2.2)

Untuk menghitung shear stress, τyx, pernyataan dα/dt harus diganti dengan besaran yang terukur. Jarak δl antara titik M dan M’ adalah δl = δu δt, atau untuk sudut kecil δl = δy δα. Sehingga diperoleh persamaan:

δα δu = , dan dengan δt δy

mengambil limit dari kedua sisi persamaan diperoleh: dα du = dt dy

(2.3)

8 Jadi, elemen fluida pada Gambar 2.2 jika dikenakan shear stress τyx akan mengalami laju deformasi (shear rate) sebesar du/dy. Fluida yang laju deformasinya sebanding proporsional dengan shear stress atau menurut persamaan berikut: τ yx ∝

du dy

(2.4)

digolongkan sebagai fluida Newton. Fluida Newton yang umum kita jumpai adalah air dan udara. Bila pada beberapa macam fluida Newton dikenakan shear stress yang sama besarnya, dan diperoleh laju deformasi yang berbeda maka fluida yang memiliki resistansi yang lebih besar untuk berdeformasi (lebih sulit berdeformasi) disebut lebih kental (viskositas lebih besar). Konstanta proporsionalitas pada persamaan (2.4) adalah viskositas absolut (dinamik), µ. Sehingga hukum viskositas Newton untuk aliran satu dimensi diberikan seperti persamaan berikut: τ yx = µ

du dy

(2.5)

Dimensi dan besaran untuk persamaan (2.5) seperti pada tabel berikut: τ (shear stress) du/dy (shear strain rate) µ (absolute viscosity)

Dimensi F/L2 1/t Ft/L2

simbol SI Pa (= N/m2) 1/s Pa ⋅ s (= N.s/m2)

Perbandingan antara viskositas absolut, µ, terhadap densitas, ρ [M/L3], disebut viskositas kinematik, ν [L2/t]. Keterangan dimensi: F = gaya; L = panjang; M = massa; t = waktu.

2.1.3

FLUIDA NON NEWTON

Fluida non-Newton adalah fluida yang shear-stress dan laju deformasinya tidak berbanding langsung secara proporsional. Cukup banyak fluida yang sering kita temui adalah fluida non-Newton, salah satunya adalah pasta gigi. Pasta gigi berlaku seperti fluida ketika ditekan keluar dari tabungnya. Tetapi pasta gigi tidak mengalir keluar sendiri ketika tutup tabung dibuka. Karena terdapat ambang tekanan (yield stress) dimana pasta gigi masih berlaku sebagai padatan. Dengan kata lain, definisi fluida diatas hanya berlaku untuk material yang mempunyai yield stress sama dengan nol. Fluida non-Newton umumnya mempunyai perilaku tergantung pada waktu (time-dependent) ataupun tidak tergantung waktu (time-independent). Salah satu

9

contoh sifat tidak tergantung pada waktu dapat dilihat pada diagram rheologi

Shear stress, τ

Bingham plastic Pseudoplastic Dilatant Newtonian

Apparent viscosity, η

(Gambar 2.3) berikut.

Pseudoplastic

Dilatant

Newtonian Deformation rate, du/dy (a)

Deformation rate, du/dy (b)

Gambar 2.3 (a) Shear stress, τ, dan (b) viskositas nyata (apparent), η, sebagai fungsi dari laju deformasi untuk aliran satu dimensi dari berbagai jenis fluida non-Newton Sedangkan hubungan antara τyx dan du/dy untuk fluida time-independent untuk aliran satu dimensi dinyatakan dalam persamaan: τ yx

 du  = k    dy 

n

(2.6)

dimana n disebut index perilaku aliran (flow behavior index), dan k adalah index konsistensi. Untuk n = 1 dan k = µ, persamaan ini menjadi persamaan viskositas Newton. Untuk memastikan τyx bertanda sama dengan du/dy, persamaan (2.6) dituliskan dalam bentuk: τ yx

du =k dy

n −1

Pernyataan η = k du dy

du du =η dy dy

n −1

(2.7)

ini disebut sebagai viskositas nyata (apparent). Sebagian

besar fluida non-Newton memiliki viskositas nyata yang relatif lebih besar dibandingkan dengan viskositas air. Fluida yang viskositas nyatanya mengecil seiring meningkatnya laju deformasi (n<1) disebut fluida pseudoplastic. Kebanyakan fluida non-Newton termasuk kelompok ini, misalnya larutan polymer, suspensi koloid dan pulp kertas dalam air. Jika viskositas nyatanya membesar terhadap peningkatan laju deformasi (n>1), fluida tersebut disebut dilatant. Suspensi tepung pati dan pasir termasuk fluida dilatant.

10

“Fluida” yang memiliki sifat padatan hingga batas minimum tekanan ambang (yield stress), τy, dilewati dan kemudian mempunyai hubungan linier antara stress dan laju deformasi, disebut sebagai Bingham plastic ideal. Persamaan shear stress untuk Bingham plastic ideal diberikan oleh persamaan (2.8) berikut:  du  du τ yx = τ y sgn   + µ p , dy  dy  du dy

= 0,

τ yx > τ y

(2.8a)

τ yx < τ y

(2.8b)

Suspensi tanah liat, lumpur sumur bor, dan pasta gigi adalah contoh fluida Bingham plastic.

2.1.4

ALIRAN INTERNAL LAMINAR PADA PELAT PARALEL TAK TERBATAS STASIONER

Aliran yang seluruhnya dibatasi permukaan solid disebut aliran internal. Aliran internal yang bergerak dengan karakter struktur alir yang berlapis-lapis (laminae) kemudian disebut aliran internal laminar. Pada aliran laminar ini tidak terjadi pencampuran secara makroskopik antara lapisan-lapisan alir yang saling berdekatan. Bila pada aliran laminar disuntikkan zat pewarna, lapisan tipis yang dibentuk zat pewarna tersebut akan tetap terlihat sebagai garis-garis tunggal; tidak terjadi dispersi zat pewarna selain difusi secara perlahan-lahan karena pergerakan molekul. Untuk aliran incompressible melalui pipa, keadaan aliran (laminar atau turbulent) ditentukan menggunakan nilai dari sebuah parameter tak berdimensi, yaitu bilangan Reynolds, Re = ρV D µ , dengan ρ adalah densitas fluida, V kecepatan aliran rata-rata, D diameter pipa dan µ adalah viskositas fluida. Pada aliran yang dibahas ini, aliran akan turbulent untuk bilangan Reynolds lebih besar dari 1400.

a 2

a

dy

τyx

p dx

Control volume y x

Gambar 2.4 Volume kontrol untuk analisis aliran laminar antara pelat paralel tak terbatas stasioner

11

Bilangan Reynolds digunakan untuk mendapatkan memeriksa apakah solusi yang didapat adalah valid. Pada Gambar 2.4, dua pelat paralel terpisah sejauh a. Pelat dianggap tak terbatas pada arah z, dan tidak memiliki perubahan properti pada arah ini. Aliran juga dianggap tunak dan incompressible. Kecepatan untuk komponen x pada pelat bagian atas dan bawah adalah nol karena keadaan tanpa slip pada dinding (Keadaan batas: pada y=0, u=0 dan pada y=a, u=0). Kecepatan tidak berubah terhadap x sehingga tergantung hanya pada y, u = u(y). Persamaan momentum pada volume kontrol adalah: Fx = 0 ; Kita lihat bahwa gaya normal (gaya tekanan/pressure) bekerja pada sisi kiri dan kanan dan gaya tangensial (gaya geser/shear) bekerja pada sisi atas dan bawah volume kontrol. Jika tekanan pada pusat elemen adalah p, maka gaya tekanan pada sisi kiri dan kanan adalah: ∂p dx   p− dydz ; dan ∂x 2  

∂p dx   − p+ dydz ∂x 2  

Jika shear stress pada pusat elemen adalah τyx maka gaya geser pada bagian sisi bawah dan atas adalah: dτ yx dy   dx dz ; dan − τ yx − dy 2  

dτ dy   τ yx + yx dx dz dy 2  

Dengan menjumlahkan gaya-gaya tersebut, diperoleh: −

∂p dτ yx + =0 ∂x dy

atau

dτ yx dy

=

∂p ∂x

(2.9)

Persamaan (2.9) harus berlaku untuk semua x dan y, untuk ini dibutuhkan dτ yx dy

=

∂p = konstan ∂x

 ∂p  dan bila diintegrasikan menghasilkan: τ yx =   y + c1  ∂x 

(2.10)

yang menunjukkan bahwa shear stress berubah secara linier terhadap y. Untuk fluida Newton, shear stress diberikan oleh persamaan (2.5), lalu: µ

du  ∂p  =   y + c1 dy  ∂x 

Diperoleh profil kecepatan: u =

1  ∂p  2 c1   y + y + c2 2µ  ∂x  µ

(2.11)

12

Konstanta c1 dan c2 diperoleh dari syarat-syarat batas. Apabila pada y=0 dan y=a , 1  ∂p  u=0, maka c2=0 dan c1 = −  a 2  ∂x  sehingga profil kecepatan: u=

2 1  ∂p  2 a 2  ∂p   y   y  − = y ay     −     2 µ  ∂x  2 µ  ∂x   a   a 

(

)

(2.12)

Distribusi shear stress, τyx, menjadi sebagai berikut: 1  ∂p   ∂p   y 1   ∂p   ∂p  τ yx =   y + c1 =   y −  a = a   −  2  ∂x   ∂x   a 2   ∂x   ∂x  r r Sedangkan laju aliran volume diberikan oleh: Q = ∫ V ⋅ dA

(2.13)

A

a

untuk kedalaman l pada arah z: Q = ∫ ul dy ; atau dengan menggunakan Persamaan 0

(2.12) diperoleh laju aliran volume untuk kedalaman l sebagai berikut: Q 1  ∂p  3 =−  a l 12µ  ∂x 

(2.14a)

Bila ∂p ∂x konstan, tekanan p berubah secara linier sepanjang x dan ∂p p 2 − p1 − ∆p = = ∂x L L disubstitusikan ke dalam (2.14) menjadikan laju aliran volume: Q a 3 ∆p = l 12 µ L

(2.14b)

Kecepatan rata-rata, V , diberikan oleh persamaan berikut: V=

Q a 2 ∆p =− A 12 µ L

(2.14c)

13

2.2 FLUIDA ELEKTRO-RHEOLOGI (FER)

2.2.1

SEJARAH PERKEMBANGAN ELEKTRORHEOLOGI*

Sejak abad 19, ilmuwan Duff (1896) dan Quinke (1897) sudah mulai mempelajari respons elektro-rheologi walaupun belum mendapatkan perhatian sebesar saat Winslow (1947) melakukan penelitian tentang fenomena elektro-viskos. Winslow memperkenalkan konsep untuk mengontrol kekentalan dari fluida elektroviskos dengan menggunakan medan listrik. Daya tahan (resistansi) alir fluida itu meningkat sebanding dengan medan listrik ketika dikenakan padanya medan listrik AC sebesar 4 kV/mm. Beliau mengamati adanya struktur berserat yang terdiri dari rantai-rantai partikel yang timbul teratur searah medan listrik yang dikenakan. Winslow menyatakan hipotesanya bahwa rantai-rantai partikel yang terinduksi medan listrik inilah yang meningkatkan kekentalan fluida. Peristiwa ini bahkan lebih dikenal dengan nama efek Winslow. Ketidakmampuan untuk mengontrol keadaan fluida dengan cepat dan tepat pada awal pengenalan teknologi ini menyebabkan fluida elektro-rheologi sulit memperoleh perhatian. Walaupun efek elektro-rheologi telah dikenal lebih dari 50 tahun, baru sekitar satu dekade terakhir para ilmuwan dan ahli rekayasa menoleh kembali kepada teknologi material ini karena kebutuhan akan peredam kejut dan suspensi aktif, serta semakin berkembangnya teknologi rekayasa smart material, menunjukkan potensi penggunaan material fluida elektro-rheologi yang besar. Bidang-bidang yang ’menjanjikan’ bagi pengembangan aplikasi teknologi ini antara lain adalah peredaman getaran mekanik, menggantikan katup-katup dan aktuator hidrolik pada pesawat terbang (aplikasi pada ruang angkasa), suspensi otomotif dan peredam getaran pada bangunan tahan gempa.

2.2.2

MEKANISME KERJA FER[5]

Ketika medan listrik eksternal dikenakan pada fluida elektro-rheologi (FER), kekentalan (viskositas) fluida meningkat dengan karakter Bingham plastic. Pada keadaan ini titik yield (batas minimum yield stress) dapat diamati sebagai fungsi dari besarnya medan listrik yang dikenakan. Dan ketika medan listrik dihilangkan, *

Dari berbagai sumber

14

Gambar 2.5 Properti stress Fluida ER viskositas fluida kembali normal pada keadaan awal dengan karakter seperti layaknya fluida Newton. Fluida elektro-rheologi (FER) dapat digolongkan secara umum ke dalam type dispersi dan type uniform. Type dispersi FER terdiri dari butiran-butiran padatan yang terdispersi pada cairan insulator (insulating oil). Butiran-butiran padat itu adalah bahan dielektrik yang mampu menghantar listrik ataupun dapat dipolarisasi. Medium pendispersi haruslah insulator (non-conducting), ataupun dapat berupa cairan organik. Sedangkan salah satu contoh type uniform FER adalah kristal cair (liquid crystal) yang bila ada medan listrik di sekitarnya koefisien viskositasnya akan meningkat. Konduktivitas FER ini biasanya sangat rendah sehingga daya listrik yang dibutuhkan untuk dapat mengubah dari Newton ke Bingham–plastic juga rendah walaupun membutuhkan tegangan listrik yang tinggi. FER adalah material yang sangat unik dan mempunyai fungsi yang fantastis. Jika sinyal listrik diberikan sebagai sinyal input pada FER, akan diperoleh keluaran berupa sinyal mekanik (viskositas fluida) dengan respons yang cepat dan reversible. Mekanisme kerja FER digambarkan seperti pada Gambar 2.6. Untuk selanjutnya butiran-butiran dielektrik disebut fasa terdispersi dan cairan insulator disebut medium pendispersi. Pada saat tidak ada medan listrik, dispersi ER berlaku seperti dispersi pada umumnya (Gambar 2.6a). Ketika medan listrik dikenakan pada dispersi ER, partikel-partikel fasa terdispersi terpolarisasi . Partikel-partikel yang terpolarisasi ini tarik-menarik karena adanya inter aksi elektrostatik. Tarik-menarik ini membentuk deretan rantai-rantai yang terjadi diantara elektroda positif dan negatif (Gambar 2.6b). Deretan rantai-rantai inilah yang menyebabkan daya melawan

15

terhadap aliran shear dan aliran tekanan, mengakibatkan terjadinya yield pada shear stress (efek ER).

Gambar 2.6 Mekanisme efek ER

2.2.3

KARAKTERISTIK FER

Karakteristik FER yang biasa dijumpai adalah sebagai berikut: 1. Shear stress terinduksi besar. 2. Kerapatan arus yang terjadi kecil, sehingga konsumsi daya juga kecil. 3. Viskositas pada medan nol adalah rendah, dan fluiditasnya baik. 4. Respons sangat cepat (dalam milidetik). 5. Durabilitas cukup baik. 6. Partikel fasa terdispersi tidak mengendap untuk waktu yang lama. 7. Dapat bereaksi terhadap medan listrik AC ataupun DC. 8. Shear stress terinduksi hampir tidak dipengaruhi shear rate. 9. Shear stress terinduksi yang besar dapat diperoleh walaupun mendapat shear rate yang tinggi. 10. Jangkauan daerah temperatur operasi lebar.

16

2.3 GETARAN MEKANIS Getaran mekanis berhubungan dengan gerak osilasi benda dan gaya yang berhubungan dengan gerak itu. Suatu gerakan benda disebut getaran jika gerakan benda itu tampak berulang kembali dalam interval waktu tertentu. Pada umumnya, getaran merupakan bentuk energi sisa, dan pada berbagai kasus getaran tersebut tidak diinginkan. Sebagai contoh getaran semacam ini adalah getaran yang dialami penumpang ataupun muatan yang berada di dalam kendaraan akibat ketidakrataan permukaan jalan. Secara umum, terdapat dua buah getaran yaitu getaran bebas dan getaran paksa. Getaran bebas terjadi bila suatu sistem berosilasi akibat gaya yang ada dalam sistem itu sendiri. Hal ini berarti tidak ada gaya luar yang bekerja pada sistem. Sistem yang mengalami getaran ini, akan bergetar pada frekuensi alaminya. Getaran paksa terjadi jika sistem mengalami eksitasi dari gaya luar. Jika gaya luar ini merupakan suatu osilasi, maka sistem akan dipaksa untuk bergetar pada frekuensi eksitasi tersebut. Jika frekuensi eksitasi sama dengan frekuensi alami dari sistem, maka akan terjadi resonansi. Resonansi akan mengakibatkan terjadinya osilasi yang besar. Osilasi yang besar biasanya tidak diinginkan, sehingga frekuensi alami menjadi hal yang sangat penting dalam masalah getaran. Dalam kondisi nyata, hampir seluruh sistem yang bergetar akan mengalami redaman. Redaman terjadi akibat energi yang terbuang menjadi bentuk energi panas yang ditimbulkan oleh gesekan. Jika redaman sistem kecil, maka pengaruh redaman sangat kecil terhadap frekuensi alami. Redaman merupakan hal yang penting dalam membatasi osilasi pada saat terjadi resonansi. Untuk menggambarkan gerak suatu sistem diperlukan sistem koordinat tertentu. Jumlah koordinat bebas yang diperlukan untuk menggambarkan gerak suatu sistem dikenal sebagai derajat kebebasan sistem. Suatu partikel bebas yang mengalami gerak umum dalam ruang akan memiliki tiga derajat kebebasan. Benda kaku akan memiliki enam derajat kebebasan, yaitu tiga komponen posisi dan tiga komponen sudut yang menyatakan orientasinya. Benda elastik kontinu akan memerlukan derajat kebebasan yang tak berhingga, mengingat benda demikian akan bergerak dalam jumlah koordinat yang tak berhingga. Dalam kebanyakan kasus, bagian dari benda dapat dianggap tegar (rigid) dan sistem secara dinamis dapat dianggap ekivalen dengan sistem yang mempunyai derajat kebebasan berhingga.

17

Untuk menurunkan persamaan gerak sistem terdapat berbagai cara seperti metoda energi, hukum gerakan Newton, persamaan Lagrange, dan lain sebagainya. Dalam perancangan ini akan digunakan hukum Newton II, yang dapat dituliskan sebagai berikut:

∑ (Gaya) = m&x& ∑ (Torsi) = Iα&& dimana :

(2.15a) (2.15b)

m = massa I = momen inersia x = perpindahan pada arah gerak massa α = sudut rotasi dari momen inersia

Apabila suatu sistem gerakan memiliki n buah derajat kebebasan, maka akan diperlukan n buah persamaan diferensial orde dua untuk menggambarkan sistem tersebut. Persamaan gerak sistem selanjutnya dapat dinyatakan dalam bentuk matriks sebagai berikut: M&x& + Cx& + Kx = F

dimana :

(2.16)

M = matriks massa

= [mij]

C = matriks redaman

= [cij]

K = matriks kekakuan = [kij] x = vektor perpindahan = [xi] F = vektor Gaya Luar

= [Fi]

i, j = 1, 2, 3, … Sistem yang memiliki n derajat kebebasan akan memiliki n buah frekuensi alami. Frekuensi alami merupakan sifat yang tidak akan berubah walaupun digunakan sistem koordinat yang berbeda untuk menggambarkan gerak sistem. Penggunaan sistem koordinat yang berbeda akan menghasilkan matriks M, K, dan C yang berbeda pula. Modus getar (mode of vibration) sistem akan berhubungan dengan frekuensi alami tersebut. Modus getar sistem itu sendiri akan ditentukan oleh kondisi awal sistem. Kondisi awal tertentu akan dapat menghasilkan sistem yang bergetar pada salah satu frekuensi alaminya saja. Modus getar ini disebut sebagai modus alami.

18

2.4 SISTEM SUSPENSI SEPEREMPAT KENDARAAN Sistem suspensi kendaraan terdiri atas susunan pegas, peredam kejut, ban dan konstruksi penopangnya. Secara kasar pembagian gaya adalah massa beban mobil ditopang oleh pegas dan peredam kejut yang disusun paralel. Pegas dan peredam kejut menumpu pada suatu sistem sumbu, dan total sistem sumbu ini bertumpu pada ban. Hubungan ini pada implementasinya dapat disusun dengan berbagai macam posisi dan ukuran. Sebuah model harus diturunkan dari sistem suspensi tersebut untuk keperluan analisis dan perhitungan. Model yang diinginkan adalah model yang dapat dengan mudah diperoleh persamaan matematisnya guna kepentingan analisis. Bila gerak yang ditinjau hanya pada arah vertikal maka model tersederhana adalah seperti terlihat pada Gambar 2.7. Model terdiri atas : 1. Massa sprung, m1, berupa massa benda tegar. 2. Massa unsprung, m2, juga berupa massa benda tegar. 3. Pegas suspensi, k1, adalah pegas murni yang memenuhi hukum Hooke: F = k .∆x 4. Pegas model, k2, representasi dari ban, adalah juga pegas murni. 5. Dan terakhir adalah peredam kejut murni yang memiliki model matematis: F = c.x& (hukum Stokes)

Mengingat semua komponen tersebut dapat dikarakterisasi oleh sebuah parameter, maka ada lima parameter penting yang harus ada. Lima parameter tersebut diperinci pada tabel 2.1. x1

m1 k1

c x2

m2 k2

x0

Gambar 2.7 Model suspensi pasif seperempat kendaraan

19

Tabel 2.1 Parameter-parameter model sistem suspensi Parameter

Simbol m1 m2 k1 k2 C

Massa sprung Massa unsprung Konstanta pegas suspensi Konstanta pegas model ban Konstanta redaman peredam kejut

Satuan Kg Kg N/m N/m Ns/m

Model sistem suspensi aktif dibangun dengan menambahkan sebuah aktuator yang dipasang paralel dengan pegas dan peredam kejut. Model konfigurasi ini adalah seperti tampak pada Gambar 2.8. Bagi sistem suspensi kehadiran aktuator memberikan tambahan gaya secara dinamis. Dengan demikian variabel penting yang menghubungkan sistem suspensi dengan aktuator adalah gaya (F).

2.4.1

DINAMIKA SISTEM SUSPENSI AKTIF MODEL SEPEREMPAT KENDARAAN

Penurunan model dinamik sistem dilakukan dengan menggunakan hukum Newton. Skema model diperlihatkan pada Gambar 2.8. x1

m1 k1

u

c x2

m2 k2

x0

Gambar 2.8 Model suspensi aktif seperempat kendaraan Untuk massa sprung, m1, berlaku : m1 &x&1 = −k1 ( x1 − x2 ) − c( x&1 − x& 2 ) + u dan untuk massa unsprung, m2, berlaku : m2 &x&2 = − k1 ( x 2 − x1 ) − c( x& 2 − x&1 ) − u + k 2 .( x0 − x2 )

Persamaan di atas diatur sehingga berbentuk seperti :

20

&x&1 = − &x&2 =

k1 k u c c x& 2 + x&1 + x1 + 1 x2 − m1 m1 m1 m1 m1

(2.17)

k k1 k + k2 u c c x& 2 − x&1 − x1 − 1 x2 + + 2 x0 m2 m2 m2 m2 m2 m2

(2.18)

Persamaan dinamika di atas dapat ditulis dalam bentuk vektor keadaan seperti berikut: − c m1   c v& =  m 2  1   0

di mana:

c

m1 − c m2 −1 1

−

k1

k1

m1

m2 0 0

   1 0   m1  0  k − 1  − 2 m 2 u +   w ; m2 v +  0  0  0      0   − 1   0  0

(2.19 a)

 v1   x&1  v   x&  2  , dan w = x& 0 v =  2 =   v3   x1 − x2      v 4   x 2 − x 0 

Bila percepatan massa sprung, defleksi ruang suspensi, dan defleksi ban adalah besaran-besaran yang dapat diukur dari sistem tersebut, maka didapatkan vektor keluaran sistem: − c m1  y= 0   0  di mana:

c

m1 0 0

− k1 1 0

m1

 0   0 v +    1  

 m1  0 u ; 0  

1

(2.19 b)

 &x&1  y =  x1 − x2  .    x2 − x0 

2.5 OPTIMASI SISTEM Dalam perancangan sistem kontrol, dikehendaki agar sistem memenuhi spesifikasi performansi tertentu. Spesifikasi performansi ini biasanya dinyatakan dengan index performansi (IP) atau kriteria optimisasi. IP adalah fungsi yang nilainya menandakan seberapa ‘baik’ performansi aktual sistem bila dibandingkan dengan performansi yang diinginkan. Pemilihan IP sangat penting, karena bentuk IP

21

menentukan hasil sistem kontrol menjadi linier, nonlinier, stasioner ataupun berubah terhadap waktu. Perekayasa kontrol merumuskan index ini berdasarkan kebutuhan. Suatu sistem kontrol dapat dinyatakan optimal bila nilai-nilai parameter yang dipilih dapat menghasilkan IP terpilih bernilai minimum atau maksimum, tergantung pada keadaan. Nilai-nilai parameter-parameter yang optimal tergantung langsung pada IP yang dipilih. IP harus mampu memberikan kejelasan antara pengaturan parameter yang optimal dan yang tidak optimal. IP juga harus menghasilkan sebuah bilangan positif atau nol. Dan agar dapat digunakan pada sistem, IP harus merupakan fungsi dari parameter-parameter sistem yang akan menghasilkan minimum atau maksimum.

2.6 SISTEM KONTROL OPTIMAL LINEAR QUADRATIC Pemilihan vektor kontrol u(t) dalam perancangan sistem pada persamaan (2.17) haruslah dapat menghasilkan nilai index performansi (IP) yang optimal. Aturanaturan linier kontrol, u(t) = –K x(t), dengan K adalah matriks r x n, dapat diperoleh dari IP, dengan batas-batas integrasi 0 dan ∞, seperti berikut: ∞

J = ∫ L(x, u) dt 0

dengan L(x,u) adalah fungsi quadratic atau fungsi Hermitian. Vektor kontrol linier u(t) = –K x(t) dapat dituliskan sebagai berikut:  k11  u1  k u   21  2  . .   =−   . .  . .     k r1 u r 

k12 k 22 . . . kr2

. .

. .

. .

.

.

.

k1n   x1  k 2 n   x 2  .  .    .  .  .  .    k rn   x n 

Perancangan sistem kontrol optimal dan sistem regulator optimal yang berdasar pada IP quadratic tersederhanakan dengan mencari elemen-elemen matriks K. Indeks performansi quadratic: ∞

J = ∫ (x * Qx + u * Ru ) dt 0

(2.20a)

dengan Q matriks Hermitian definit non-negatif atau matriks simetris riil, R matriks Hermitian definit positif atau matriks simetris riil. Perlu dicatat bahwa IP quadratic

22

riil adalah IP quadratic kompleks (IP Hermitian) untuk kasus khusus. IP untuk sistem yang memiliki vektor dan matriks riil dapat dituliskan dengan persamaan berikut: ∞

(

)

I = ∫ x T Qx + u T Ru dt 0

(2.20b)

Tinjau kembali persamaan sistem ruang keadaan baku seperti berikut: x& (t ) = Ax(t ) + Bu(t )

(2.21a)

y (t ) = Cx(t ) + Du(t )

(2.21b)

Matriks K dari vektor kontrol optimal: u

x

x& = Ax + Bu

-K

Gambar 2.9 Sistem kontrol optimal u(t ) = − Kx (t )

(2.22)

harus dapat meminimisasi IP quadratic diatas, persamaan (2.20a). Matriks Q dan R menyatakan hubungan kepentingan antara galat (error) dan pengeluaran energi sinyal kontrol. Jika elemen-elemen yang tidak diketahui dari matriks K dapat ditentukan dan juga berhasil meminimisasi nilai IP, maka u(t ) = − Kx (t ) adalah optimal untuk setiap keadaan awal x(0). Diagram blok pada Gambar 2.9 menunjukkan konfigurasi optimal tersebut. Jika persamaan (2.22) disubstistusikan ke persamaan (2.21a) akan diperoleh hubungan berikut: x& = Ax − BKx = (A − BK )x Dalam hal ini matriks (A – BK) diasumsikan stabil, atau memiliki nilai-nilai eigen yang riil negatif. Jika persamaan (2.22) disubstitusikan ke persamaan (2.20a), akan diperoleh Indeks Performansi yang dapat dituliskan sebagai berikut: J=

∞

∞

∫0 (x * Qx + x * K * RKx ) dt = ∫0

x * (Q + K*RK )x dt

(2.23)

Guna menyelesaikan masalah optimisasi parameter, digunakan persamaan (2.24) untuk menghasilkan persamaan (2.25).

23

x * (Q + K * RK)x = −

d (x * Px ) dt

(2.24)

x * (Q + K * RK)x = − x& * Px − x * Px& = − x * [(A − BK ) * P + P(A − BK )]x (2.25) Persamaan yang terakhir ini harus berlaku untuk setiap x dan untuk itu dibutuhkan persamaan (2.26).

(A − BK ) * P + P (A − BK )

= − (Q + K * RK )

(2.26)

Dengan menggunakan metoda Liapunov kedua, jika (A – BK) adalah matriks yang stabil maka ada sebuah matriks definit positif P yang memenuhi (2.26). Selanjutnya persamaan (2.24) diintegralkan menurut persamaan (2.23) dengan memperhatikan bahwa x(∞) = 0, maka IP dapat dituliskan sebagai berikut: J = x * (0) P x(0)

(2.27)

Persamaan-persamaan berikut ini akan memberikan solusi permasalahan kontrol optimal linear quadratic. Matriks optimal K diberikan oleh persamaan berikut: K = R −1 B * P

(2.28)

Sehingga vektor kontrol optimal diberikan oleh persamaan: u(t ) = − K x(t ) = − R −1B * P x(t )

(2.29)

Matriks P pada persamaan (2.28) harus dapat memenuhi persamaan (2.26) atau persamaan tereduksi ini: A * P + PA − PBR −1B * P + Q = 0

(2.30)

Persamaan (2.30) ini dikenal sebagai persamaan matriks–tereduksi Riccati. Langkah-langkah perancangan sistem kontrol optimal linear quadratic adalah sebagai berikut: 1. Pecahkan persamaan (2.30), persamaan Riccati, untuk mendapatkan matriks P. 2. Substitusikan matriks P tersebut ke dalam persamaan (2.28). Matriks K yang didapat adalah yang optimal. Persyaratan stabil dari (A – BK) ekivalen dengan memeriksa matriks berikut ini:

24

 Q2   12   Q A    . Rank   =n .     .  1 n −1  Q 2 A  1

(2.31)

Keadaan rank dari matriks pada persamaan (2.31) yang sama dengan orde system (n) adalah cukup untuk memeriksa kestabilan matriks (A – BK). Pendekatan lain untuk mendapatkan matriks penguatan umpan-balik optimal K: 1. Matriks P pada persamaan (2.26) dijadikan sebagai fungsi dari K. 2. Substitusi matriks P tersebut pada persamaan (2.27) sehingga IP menjadi fungsi dari K. 3. Hitung elemen-elemen K sehingga IP terminimisasi. Minimisasi J terhadap elemen-elemen kij dari K dapat diperoleh dengan mengeset ∂J ∂k ij sama dengan nol dan menyelesaikannya untuk mendapatkan nilai optimal kij.

2.6.1

INTEGRAL PERFORMANSI SECARA UMUM

Banyak kasus dapat diformulasikan dengan indeks performansi seperti pada persamaan (2.20), tetapi dalam beberapa kasus dijumpai term persilangan 2x * N * u = x * N * u + u * Nx

dalam integral performansi.

Penguatan optimum

dalam kasus seperti ini dapat dicari dengan metode yang sama seperti pada sub bab sebelum ini. Penguatan optimum akan didapatkan seperti pada persamaan (2.32). ˆ = R −1 (B * P + N) K

(2.32)

dimana matriks P memenuhi persamaan matriks Riccati (2.33) berikut ini: A * P + PA − PBR −1 B * P + Q = 0 dengan

(2.33)

A = A − B R −1 N

(2.34)

Q = Q − N* R −1 N

(2.35)

Dari sini dapat dicoba cara lain untuk membuktikan penguatan tersebut. Sinyal kontrol u dapat diambil seperti pada persamaan (2.36) u = v − R −1 N x

(2.36)

25

Lalu disubstitusikan ke dalam persamaan umum ruang keadaan sistem seperti di bawah ini: x& = ( A − BR -1 N)x + B v = A x + B v

(2.37)

Indeks/integral performansi yang akan diminimalkan adalah: ∞

J = ∫ (x* Q x + x* N* u + u* N x + u* R u ) dt

(2.38)

0

IP pada persamaan (2.38) bila menggunakan persamaan (2.36) akan menjadi seperti persamaan berikut ini:

(

) (

)

(

) (

x* Q x + x* N* v − R −1 N x + v * − x * N * R −1 N x + v * −x * N * R −1 R v − R −1 N x = x* (Q − N* R N ) x + v* Rv

(2.39)

Usaha meminimalkan nilai dari persamaan (2.38) untuk proses sebenarnya adalah ekivalen dengan minimisasi dari persamaan (2.40): J = ∫ (x * Q x + v * Rv ) dt ∞

(2.40)

0

x& = A x + B v . Nilai minimum dari J dapat

untuk persamaan matriks proses

diperoleh dari vektor kontrol yang menggunakan matriks K dari sub bab terdahulu. Persamaan (2.41) berikut ini adalah vektor kontrol untuk mengoptimalkan persamaan (2.40). v= −Kx

(2.41)

dimana matriks penguatan untuk v diberikan oleh: K = R −1 B * P . Matriks P harus memenuhi persamaan (2.33) agar menghasilkan persamaan berikut ini:

(

)

ˆ x u = − R −1 B * P − R −1 N x = − K ˆ diberikan pada persamaan (2.32). dengan K

(2.42)

)

26

BAB III PEMODELAN SISTEM

3.1 KARAKTERISTIK BAHAN FLUIDA ELEKTRORHEOLOGI Fluida elektrorheologi yang akan digunakan pada tugas akhir ini adalah fluida TX-ER6 (Nippon Shokubai Co., Ltd., 1999). Typical data yang diberikan pihak manufaktur antara lain adalah diagram shear stress terinduksi terhadap medan listrik, seperti pada gambar 3.1 berikut.

Gambar 3.1 Shear stress terinduksi terhadap medan listrik pada laju geser 200/sec untuk TX-ER6[5] Yield stress dinamik untuk fluida tersebut dapat didekati sebagai persamaan polinom (derajat banyak). Setelah membandingkan beberapa orde pendekatan kurva dengan memperhatikan besar residual norm, penulis mengambil pendekatan polinom orde lima sehingga diperoleh fungsi medan listrik seperti berikut: τ y ( E ) = a5 E 5 + a4 E 4 + a3 E 3 + a 2 E 2 + a1 E + a0

(3.1)

27

Koefisien polinom dapat diperoleh menggunakan dengan metode least-squares fit, yang cara mendapatkannya dapat dilihat pada lampiran A.1. Untuk medan listrik AC 50 Hz: a0 = 2,7645 Pa, a1 = 40,2403 Pa mm kV-1, a2 = – 27,5004 Pa mm2kV-2, dan a3 = 51,9837 Pa mm3 kV-3, a4 = – 13,3152 Pa mm4 kV-4, a5 = 1,2737 Pa mm5 kV-5 dengan residual norm = 5.2968 x 10-012. Bagi medan listrik DC: a0 = 11.0263 Pa, a1 = – 82,0472 Pa mm kV-1, a2 = 237,9117 Pa mm2 kV-2, dan a3 = – 126,6655 Pa mm3 kV-3, a4 = 38,7015 Pa mm4 kV-4, a5 = –3,6857 Pa mm5 kV-5 dengan residual norm = 1,5664 x 10-011. Dari pihak manufaktur disebutkan bahwa TXER6 dapat dianggap mempunyai viskositas absolut nominal, µ0, konstan sebesar 0,035 Poise. Tabel 3.1 berikut ini memudahkan penulis untuk meninjau koefisienkoefisien polinom tersebut.

Tabel 3.1 Koefisien polinom Yield Stress sebagai fungsi dari Medan Listrik

AC 50 Hz

a0 2,7645

a1 40,2403

a2 –27,5004

a3 51,9837

a4 –13,3152

a5 1,2737

DC

11,0263

–82,0472

237,9117

–126,6655

38,7015

–3,6857

3.2 SKETSA MODEL PEREDAM KEJUT ELEKTRORHEOLOGI Petugas akhir mencoba mengembangkan model matematis peredam kejut FER sebagai dasar untuk perancangan lebih jauh. Penurunan model matematis peredam ini menggunakan mekanisme aliran fluida Newton dan plastik-Bingham untuk memperoleh harga peredaman viskos ekivalennya. Asumsi-asumsi dan batasan-batasan yang digunakan pada pemodelan ini diantaranya adalah: 1. Analisis dilakukan dengan pendekatan aliran pelat paralel stasioner 1 dimensi dan dengan batasan bahwa celah elektroda relatif cukup kecil terhadap diameter piston. 2. Pada analisis aliran geser Newton, aliran adalah quasi-steady dan berkembang penuh. 3. Tidak ada aliran yang melewati kepala piston.

28

4. Rugi-rugi tahanan gesek fluida yang keluar dan masuk celah saluran diabaikan. 5. Yield stress dinamik dan viskositas plastik adalah hanya merupakan fungsi dari tegangan listrik yang diaplikasikan. Dan hanya keduanya-lah properti dari material yang diperhatikan. 6. Inersia fluida diabaikan. 7. Komponen-komponen peredam dianggap tidak bermassa. 8. Tekanan yang didapat dari gerakan piston berubah linier sepanjang celah elektroda. Sketsa gambar konsep peredam kejut FER[1] yang digunakan adalah seperti berikut: FR d1

sekat

batang piston fluida ER saluran ER

sumber V tegangan

silinder luar piston silinder dalam

P1 d2 P2 d3

L

a GAS

piston mengambang Gambar 3.2

Gambar konsep peredam kejut yang berisi FER

Pada penurunan model matematis, reservoir gas seperti pada gambar konsep tidaklah diperhitungkan. Reservoir gas, yang secara nalar diberi tekanan yang

29

lebih besar dari drop tekanan maksimum yang diperkirakan terjadi pada piston, akan menaikkan tekanan total di dalam peredam kejut. Hal ini akan menghindari terjadinya tekanan di sebelah atas piston lebih kecil dari tekanan atmosfer sehingga mencegah terjadinya kavitasi selama kompresi dan mencegah udara terambil masuk lewat sekat batang piston ke dalam. Pada rancangan ini dikembangkan gaya redaman yang bergantung pada kecepatan karena adanya drop tekanan sepanjang celah elektroda begitu piston peredam dikenakan gaya.

3.3 PEMODELAN PEREDAM KEJUT MENGGUNAKAN PENDEKATAN GEOMETRI 1D PELAT PARALLEL Kita lihat kembali Gambar 2.4 yang ditampilkan pada Gambar 3.3 berikut ini. a 2

a

dy

τyx

p

Control volume y

dx

Gambar 3.3

x

Volume kontrol untuk analisis aliran antara pelat

paralel tak terbatas stasioner Persamaan (2.9) yang merupakan persamaan kesetimbangan gaya, juga akan ditampilkan kembali pada persamaan berikut ini: dτ

yx ∂p = = konstan dy ∂x

Tekanan, p, diasumsikan berubah linier sepanjang celah dan jarak ∆x = L , maka persamaan diatas dapat dituliskan seperti pada persamaan (3.2). dτ ∆p = = konstan dy ∆x

⇔

dτ ∆p = = konstan dy L

(3.2)

Persamaan (3.2) bila diintegralkan akan diperoleh persamaan berikut: τ=

∆p y + c1 L

(3.3)

30

Persamaan (3.2) dan (3.3) ini adalah acuan awal analisis yang akan dilakukan pada sub bab selanjutnya.

3.3.1

FER SEBAGAI FLUIDA NEWTON

Konstanta peredaman viskos, C, akan dicari menggunakan mekanisme aliran Newton, yang adalah pendekatan model FER tanpa medan listrik. Pada kasus fluida Newton, shear stress τ proporsional terhadap gradien kecepatan pada celah, seperti berikut: τ = µ0

du dy

(3.4)

dimana µ0 adalah viskositas absolut (pada laju geser nol), yang dapat digunakan sebagai pendekatan sifat FER saat tanpa ada medan listrik. Substitusikan (3.4) ke (3.4) lalu diintegrasi untuk memperoleh profil kecepatan: µ0

du ∆p = y + c1 dy L

u( y ) =

∆p 2 c1 y + y + c2 µ0 2µ0 L

(3.5)

c1 dan c2 dapat dihitung dari syarat-syarat batas: u(y = 0) = 0 dan u(y = a) = 0. c1 = −

a∆p ; 2L

c2 = 0

Profil kecepatan pada celah elektroda akan menjadi tampak seperti pada persamaan (3.6). u N ( y) =

∆p (y 2 − ay ) 2µ0 L

(3.6)

Dari prinsip kontinuitas, fluks volume yang melalui celah elektroda QN, harus sama dengan fluks volume yang dipindahkan kepala piston, QP = A v0, dimana A adalah luas penampang kepala piston dan v0 kecepatan piston. l∆p ( Q N = ∫ u N ( y ) l dy = y 2 − ay ) ∫ 2µ0 L 0 0 a

a 3l QN = − ∆p 12µ 0 L

a

(3.7)

31

dimana:

a

: jarak antar dinding elektroda

l

: kedalaman saluran

L

: panjang saluran

∆p : beda tekanan pada piston µ0 : viskositas absolut fluida Karena QN = QP, dan mengingat bahwa hubungan beda tekanan dan gaya F adalah ∆p = − F A , maka diperoleh ∆p = −

12 µ 0 L Av0 a 3l

F = C N v0

(3.8) (3.9)

dimana konstanta peredaman (damping) viskos CN adalah sebagai berikut: C N = ΓN µ 0

(3.10)

dengan ΓN merupakan faktor geometri peredam kejut: 12 LA2 ΓN = a 3l

3.3.2

(3.11)

FER SEBAGAI PLASTIK BINGHAM Material plastik Bingham berlaku seperti benda padat di bawah batas shear

stress tertentu yang disebut yield stress. Keadaan pada tingkat stress ini biasanya disebut keadaan pre-yield. Pada keadaan pre-yield ini, material kaku (rigid) dan tidak mengalir. Ketika shear stress melampaui yield stress, material disebut berada pada keadaan post-yield dan berlaku seperti fluida viskos. Atau dapat dikatakan bahwa shear stress pada material harus melebihi yield stress dinamik sebelum material dapat mengalir. Perilaku material ini dapat dilihat seperti pada Persamaan (2.8). Yield stress dinamik, τy, diasumsikan semata-mata adalah fungsi polinom dari medan listrik. Dan µp adalah viskositas plastik dan diasumsikan tidak bergantung pada besar medan listrik untuk menyederhanakan analisis. Aliran Newton dapat dilihat sebagai kasus khusus dari aliran plastik Bingham dengan yield stress dinamik sama dengan nol.

32

a2

δp

Daerah 2

a

y

Daerah 3

a1

Daerah 1

x

Gambar 3.4

Profil kecepatan yang dapat terjadi pada pelat parallel

Profil kecepatan pada pelat elektroda parallel untuk plastik Bingham dengan adanya gradien tekanan linier sepanjang sumbu x dapat dilihat pada gambar (3.4). Profil aliran dibagi menjadi 3 daerah. Pada daerah pertama dan ketiga, shear stress lebih besar dari yield stress, sehingga kedua daerah ini berada pada keadaan post-yield. Daerah kedua adalah bersifat padat dan kaku karena shear stress padanya lebih kecil dari pada yield stress. Daerah 2 disebut pada pre-yield dan mengalami aliran padatan (plug). Masing-masing daerah akan dilihat tersendiri profil kecepatannya. Daerah 1 (post-yield): 0 < y < a1, τ > τ y , du dy > 0 , sehingga shear stress: τ =τy + µ

du dy

(3.12)

Persamaan ini dimasukkan ke dalam persamaan (3.3): τy + µ

du ∆p = y + c1 dy L

persamaan diatas dapat dituliskan sebagai: µ

du ∆p = y + c *1 , dy L

dengan

c *1 = c1 − τ y

(3.13)

Persamaan (3.13) bila diintegralkan untuk memperoleh profil kecepatan akan menjadi seperti berikut ini: ∆p 2 c *1 u= y + y + c2 2 µL µ dimana c*1 dan c2 dapat dihitung dari syarat-syarat batas di bawah ini: u ( y = 0)

=0

du dy

=0 y =a1

(3.14)

33

Aliran padatan berkecepatan konstan sehingga percepatan pada titik pertemuan antara fluida alir dan padatan haruslah nol, dan kecepatan pada dinding stasioner adalah nol. Syarat batas tersebut akan memperoleh persamaan berikut ini: c *1 = −

a1 ∆p L

⇒

c1 = τ y −

a1 ∆p L

c2 = 0 Profil kecepatan pada daerah 1 akan menjadi seperti yang tertera pada persamaan (3.15). a − 1 ∆p ∆ p y, u 1B ( y ) = y2 + L 2 µL µ ∆p 2 (y − 2a1 y ) u 1B ( y ) = 2 µL

atau

(3.15)

Fluks volume pada daerah 1 dapat diketahui seperti berikut:  π∆p 3  Q = ∫ u ( y ) 2πy dy = ∫  y − 2a1 y 2 )dy ( µL  0 0 a1

B 1

a1

B 1

a1

π∆p 1 4 2 5π ∆p 4 3 ( Q = a1 4 y − 3 a1 y ) = − µL 12 µL 0 B 1

(3.16)

Daerah 2 (pre-yield): a1 < y < a2, τ < τ y , du dy = 0 Ini adalah daerah aliran padat (plug), dan dari du dy = 0 diketahui bahwa u 2B = konstan = u p

(3.17)

dimana up adalah kecepatan plug, dan shear stress tampak pada persamaan (3.18). τ =

∆p y + c1 , L

(3.18)

Pada batas-batas plug, shear stressnya adalah sama dengan yield stress FER. τ ( y = a1 ) = − τ y

dan

τ ( y = a2 ) = τ y

Shear stress tersebut di atas ini bila masing-masing dimasukkan ke dalam persamaan (3.18) akan menghasilkan persamaan berikut: −τ y =

∆p a1 + c1 L

τy =

∆p a 2 + c1 L

34

Persamaan pertama dikurangi oleh persamaan kedua untuk menghilangkan c1 akan menghasilkan persamaan (3.19). τy =

∆p ( a 2 − a1 ) 2L

(3.19)

Memecahkan persamaan di atas untuk mencari tebal plug, δp = (a2 – a1), dan dengan mengingat bahwa ∆p = − F A pada kasus ini adalah negatif, maka akan diperoleh persamaan (3.20). a 2 − a1 =

2 Lτ y ∆p

=

2 LAτ y F

=δp

(3.20)

Fluks volume pada daerah 2 diperoleh sebagai berikut: Q = ∫ 2πy u2B ( y ) dy = u pπ (a 22 − a12 ) a2

B 2

(3.21)

a1

Daerah 3 (post-yield): a2 < y < a, τ > τ y , du dy < 0 , sehingga: τ = −τ y + µ

du dy

(3.22)

Persamaan (3.22) dimasukkan ke persamaan (3.3) menjadi: −τ y + µ

du ∆p y + c1 = dy L

⇔

µ

du ∆p y + c **1 = dy L

(3.23)

dimana c**1 = c1 + τy. Lalu diintegralkan untuk memperoleh profil kecepatan: u=

∆p 2 c **1 y + y + c3 2 µL µ

(3.24)

Dengan menggunakan syarat-syarat batas berikut: u( y = a )

=0

du dy

=0 y =a2

akan diperoleh solusi untuk konstanta-konstanta yang belum diketahui pada persamaan (3.24), seperti berikut: c **1

=−

∆p a 2 , dan L

35

=

c3

− ∆p 2 ∆pa 2 ∆p ( a= a + − a 2 + 2a 2 a ) µL 2 µL 2 µL

Profil kecepatan pada daerah 3 adalah: ∆p 2 ∆p ∆p ( a2 y + y − − a 2 + 2a 2 a ) 2 µL µL 2 µL ∆p 2 ( = y − 2a 2 y + 2a 2 a − a 2 ) 2 µL

u 3B ( y ) =

(3.25)

Fluks volume yang melalui daerah 3: Q3B = ∫ 2πy u3B ( y ) dy =

π ∆p ( y 3 − 2a 2 y 2 + (2a 2 a − a 2 )y ) dy ∫ µL a 2

]

a

a2

a

Q3B =

π ∆p µL

[

Q3B =

π ∆p µL

[(

Q3B =

π ∆p 1 4 1 (− 4 a + 3 a2 a 3 + 12 a22 a 2 − a23a + 125 a24 ) µL

y 4 − 23 a 2 y 3 + (a 2 a − 12 a 2 )y 2

1 4

1 4

a a2

a 4 − 23 a2 a 3 + a 2 a 3 − 12 a 4 ) − (14 a24 − 23 a2 a 23 + a 23a − 12 a 22 a 2 )

]

(3.26)

Pada persamaan-persamaan di atas dijumpai dua bilangan tak diketahui, yaitu a1 dan a2, dan kita hanya menjumpainya dalam satu persamaan (3.20): a 2 − a1 =

2 LAτ y F

⇔

a1 = a 2 −

2 LAτ y F

(3.27)

Persamaan lain dapat dibangun dari kecepatan plug, up, yang konstan: u p = u1B ( a1 ) = u 3B ( a 2 ) = konstan

(3.28)

atau lebih jelasnya dapat dilihat pada persamaan berikut: u 3B ( a 2 ) − u1B ( a1 ) = 0

(3.29)

Dapat dilihat lagi bahwa kecepatan pada a1 adalah seperti pada persamaan (3.30): u1B ( a1 ) =

∆p 2 ∆p 2 ( a1 − 2a12 ) = − a1 2 µL 2 µL

Dan kecepatan pada a2 tampak seperti berikut:

(3.30)

36

∆p 2 ∆p ∆p a2 − 2a 22 + − a 2 + 2a 2 a 2 µL 2 µL 2 µL ∆p 2 a 2 − 2a 2 a + a 2 =− 2 µL

(

u 3B (a 2 ) =

(

)

)

(3.31)

Menggunakan persamaan (3.30) dan (3.31), berarti (3.29) dapat dituliskan seperti berikut: −

∆p 2 ∆p 2 a 2 − 2a 2 a + a 2 + a1 = 0 2 µL 2 µL

(

)

∆p 2 a1 − a 22 + 2a 2 a − a 2 = 0 2 µL

(

⇔

)

(3.32)

Persamaan (3.27) bila disubstitusikan ke persamaan (3.32) akan diperoleh persamaan berikut ini: −

4 L2 A 2τ y2 F

2

+

4 LAτ y F

a 2 − 2a 2 a + a 2 = 0

Persamaan di atas harus diselesaikan untuk memperoleh a2 berikut ini: a2 =

4 L2 A2τ y2 − a 2 F

2

(3.33)

4 F LAτ y − 2 F a 2

sehingga dapat diperoleh a1 seperti dituliskan di bawah ini: a1 =

4 L2 A 2τ y2 − a 2 F

2

4 F LAτ y − 2 F a 2

−

2 LAτ y F

(3.34)

Begitu a1 dan a2 dapat dihitung maka profil kecepatan dapat diketahui dari persamaan (3.15), (3.25) dan (3.28) sehingga didapatkan distribusi kecepatan pada tiap daerah.

37

Kecepatan piston dapat dihitung dari fluks volume yang dipindahkan piston, Qp, yang harus sama dengan fluks volume yang melalui celah elektroda, QB. Qp = QB

(3.35)

Total fluks volume yang melalui celah elektroda ialah: QB = Q1B + Q2B + Q3B

(3.36)

dimana fluks volume masing-masing daerah dapat diketahui dari (3.16), (3.21) dan (3.26). Kecepatan piston dapat dituliskan pada persamaan (3.37): Q1B + Q2B + Q3B v0 = A

(3.37)

Jadi redaman viskos ekivalen, CB, sama dengan: C B =

3.3.3

F v0

(3.38)

PARAMETER-PARAMETER MODEL PELAT PARALEL

Metodologi pemecahan persamaan-persamaan hasil sub bab sebelum ini dapat dijelaskan sebagai berikut: Bila diberikan dimensi peredam kejut, properti fluida dan gaya, persamaan (3.27) dan (3.28) digunakan untuk mencari batas-batas plug (a1 dan a2). Begitu a1 dan a2 diketahui maka persamaan-persamaan (3.16), (3.21) dan (3.26) dapat dihitung sehingga fluks total volume yang melalui celah elektroda dapat diperoleh dari persamaan (3.36). Fluks total volume yang telah diketahui itu kemudian digunakan untuk mencari kecepatan piston dengan persamaan (3.37) dan barulah redaman viskos ekivalennya diketahui melalui persamaan (3.38). Bila kita lihat kembali model peredam kejut elektrorheologi pada awal bab ini, drop tekanan pada perangkat ini dipengaruhi oleh dua komponen, yaitu komponen viskos yang independent terhadap medan listrik ∆p N dan komponen yield stress terinduksi yang dependent pada medan listrik ∆pτ . Dari persamaan (3.8) dan (3.20) drop tekanan yang diberikan oleh fluida adalah: ∆P = ∆p N + ∆pτ = − dimana:

12 µ 0 L A 2L v0 − τ (E ) 3 (a2 − a1 ) y a l

a : jarak antar dinding elektroda

(3.39)

38

a2, a1 : jarak-jarak batas daerah alir plug l : kedalaman celah L : panjang saluran ∆p: beda tekanan pada piston µ0 : viskositas absolut fluida A : luas penampang kepala piston ν0 : kecepatan gerak piston τy : yield stress, yang merupakan fungsi polinom medan listrik dan diasumsikan bahwa pada saat FER sebagai fluida Bingham, fluks volume yang melalui daerah 1 dan 3 (perbatasan fluida terhadap dinding) jauh lebih kecil dari fluks volume alir plug. Hal ini dapat terjadi bila a1 sangat tipis ( a1 → 0 ) dan a2 mendekati a ( a 2 → a ) sehingga fluks volume alir plug-lah yang paling berpengaruh terhadap terjadinya drop tekanan pada perangkat ini. Dengan demikian persamaan (3.39) dapat dituliskan sebagai berikut: ∆P = ∆p N + ∆pτ = −

12 µ 0 L A 2L v0 − τ y (E ) 3 a l a

(3.40)

Adapun dimensi peredam kejut yang dipergunakan dalam pemodelan ini seperti yang terlihat pada gambar 3.2 mengikuti dimensi prototype peredam kejut Petek, K. Nicholas[2], yaitu (dalam mm): d1 = 12,7

d2 = 42,5

d3 = 46,6 L = 63,5

a = 0,5

Dengan kedalaman saluran l menggunakan keliling dari diameter (d3 + a = 47,1 mm), yaitu l = 147,1836 mm. Selanjutnya akan digunakan medan listrik DC, untuk menyederhanakan analisis, dan persamaan 3.40 yang menghasilkan persamaan berikut: ∆P = ∆p N + ∆pτ = −

12 µ 0 L A 2L v0 − τ y (E ) 3 a l a

dengan parameter-parameter persamaan sebagai berikut: µ0 = 0,035 Pa.s;

L

A

= ¼ . π . (42,5)2 mm2 = 1418,625 mm2; a3

l

= 147,1836 mm;

= 63,5 mm; = 0,125 mm3;

39

akan menjadikan persamaan tersebut seperti persamaan (3.41) berikut: ∆P = ∆p N + ∆pτ = − 2056,465 v0 − 254 τ y (E )

(3.41)

Atau, untuk menyatakannya dalam persamaan gaya, kita lihat kembali persamaan (3.9): FN = C N v0 ;

dimana

(

C N = 12 L A 2 µ 0

) (a l ) 3

dan persamaan (3.20), dimana gaya F terhadap drop tekanan ∆P dinyatakan sebagai dituliskan di bawah ini: F = − ∆P A sehingga Fτ = {(2 LA)τ y a} sgn(v0 ) .

Kemudian persamaan gaya sistem ini dapat dituliskan sebagai berikut: 2L A  F = C N v0 +  τ y  sgn(v0 )  a  dimana:

(3.42a)

C N = 2917353,24 N.s.(mm)-1 = 2917,35324 N.s.m-1;

(2 L A a ) = 360330,85988 mm-2 = 0,36033 m-2; dan 1,  sgn(v0 ) = 0, − 1, 

v0 > 0 v0 = 0 v0 < 0

Persamaan (3.42a) dapat dituliskan sebagai berikut: F = 2917,35324 v0 + 0,36033τ y sgn(v0 )

(3.42b)

F adalah gaya yang diberikan oleh sistem dengan masukan kecepatan piston v0 dan yield stress τy, yang merupakan fungsi polinom dari medan listrik DC. Lengkapnya persamaan gaya dapat dituliskan kembali seperti berikut ini: F = 2917,3532 v0 + 0,3603 (−3,6857 E 5 + 38,7015 E 4 − 126,6655 E 3 + 237,9117 E 2 − 82,0472 E + 11,0263) sgn(v0 ) (3.42c) dengan batasan pada E: 0 ≤ E ≤ 4 kV/mm . E = 0:

F = 2917,35 v0 + 3,97 sgn(v0)

(3.43a)

E maksimal = 4 kV/mm:

F = 2917,35 v0 + 546,37 sgn(v0)

(3.43b)

40

Koefisien redaman ekivalennya, CE, dapat dicari dari hubungan F = CE v0, sehingga diperoleh persamaan berikut ini: C E = 2917,3532 +

0,3603 (−3,6857 E 5 + 38,7015 E 4 − 126,6655 E 3 + 237,9117 E 2 v0 − 82,0472 E + 11,0263) sgn(v0 )

(3.42d)

untuk v0 tidak sama dengan nol. 3,9728 sgn(v0) v0

Bila E = 0 kV/mm:

CE=0 = 2917,3532 +

E maksimal = 4 kV/mm:

CEmaks = 2917,3532 +

(3.43c)

546,3539 sgn(v0) v0

(3.43d)

3.4 PEMODELAN PEREDAM KEJUT MENGGUNAKAN PENDEKATAN GEOMETRI 1D CYLINDRICAL AXISYMMETRIC Peredam kejut yang hendak dikembangkan memiliki bentuk pipa, maka selanjutnya dirasa lebih mengena bila disertakan juga pemodelan axisymmetric 1D dalam tugas akhir ini. Model itu dapat kita lihat pada gambar (3.5). Annular control volume

R

r x

p, τrx dr dx

dr Annular control volume

r

Gambar 3.5 Volume kontrol untuk analisis aliran laminar penuh pada pipa Persamaan kesetimbangan gaya pada kontrol volume tersebut adalah ρ

∂u ∂τ rx τ ∂p + + = ∂t ∂r r ∂x

(3.44)

Dengan u adalah kecepatan, τ shear stress, r koordinat radial, x koordinat longitudinal, dan p adalah tekanan yang didapat dari gerakan kepala piston. ρ adalah densitas fluida. Asumsi kita adalah bahwa aliran tunak dan p berubah

41

secara linier sepanjang celah elektroda sehingga persamaan (3.44) dapat ditulis sebagai berikut: dτ τ ∆p + = dr r L

(3.45)

Persamaan ini adalah acuan awal untuk mengembangkan persamaanpersamaan berikutnya.

3.4.1

ALIRAN GESER NEWTON

Mekanisme aliran geser Newton digunakan sebagai pendekatan terhadap model peredam kejut FER tanpa adanya medan listrik. Shear stress, τN, pada fluida Newton adalah proporsional terhadap gradien kecepatan pada celah. τ N = µ0

du = µ 0γ& dr

(3.46)

dimana γ& adalah laju shear strain. Bila persamaan (3.46) disubstitusikan pada persamaan (3.45): µ0

d 2 u µ 0 du ∆p + = r dr L dr 2

(3.47)

Profil kecepatan pada celah sebagai fungsi jari-jari silinder dapat diperoleh dengan mengintegrasi persamaan diatas: u ( R) =

∆P 2 r + D1 ln r + D0 4µ 0 L

(3.48)

D1 dan D0 adalah konstanta integrasi yang diperoleh dari syarat batas peredam kejut. Kemudian dari Gambar 3.1 sebagai bentuk sketsa dan menggunakan syarat batas: u(R1) = 0

u(R2) = 0

profil kecepatan pada celah adalah u NR (r ) = Fluks

volume

ln( R2 r )  ∆P  2 2 ln( R1 r ) − R12  r + R2  4µ 0 L  ln( R2 r ) ln( R2 R1 )  yang

melalui

celah

elektroda

mengintegrasikan profil kecepatan yang ada pada celah itu

(3.49) diketahui

dengan

42

R2

Q = ∫ 2π r u (r ) d r

(3.50)

R

lalu menggunakan persamaan (3.49) pada persamaan (3.50), kemudian mengintegrasikan persamaan (3.50) dan menyederhanakannya, diperolehlah fluks volume yang melalui celah elektroda, sebagaimana tertera pada persamaan (3.51). Q NR

∆P A22  R1 =  8π µ 0 L  R2 

[

2  1 − (R1 R2 )  + ln (R2 R1 )  4

]

2

 − 1 

(3.51)

dimana A2 = π R22 . Fluks volume yang melalui celah elektroda harus sama dengan fluks volume yang dipindahkan oleh kepala piston, QP = A v0, dimana A luas kepala piston. Lalu pemecahan gaya F peredam kejut adalah seperti berikut: F = CNR v0,

(3.52)

dimana konstanta peredaman CNR dinyatakan oleh persamaan (3.53). C NR = µ 0 ΓNR

(3.53)

ΓNR hanya tergantung pada geometri aliran peredam kejut, dituliskan pada persamaan (3.54). ΓNR

 A = 8π L   A2

  

2

  R 1 1 −    R2

[

2  1 − (R1 R2 )  − ln (R2 R1 )  4

]

2

  

−1

(3.54)

Terlihat pada model ini bahwa gaya fluida Newton adalah hasil perkalian dari viskositas absolut, µ0, geometri peredam kejut, ΓNR , dan kecepatan piston, v0. Konstanta peredaman viskos dapat dimaksimalkan dengan menerapkan strategi seperti berikut ini (1) memaksimalkan panjang elektroda, L; (2) mengurangi jarak celah, atau memaksimalkan perbandingan luas A/A2; dan (3) memperbesar viskositas, µ0, fluida.

43

3.4.2

ALIRAN FLUIDA BINGHAM

R2

Daerah 3 a

Daerah 2 Daerah 1

Rpo

δp Rpi

R1 r x

Gambar 3.6

Profil kecepatan pada celah melingkar untuk bahan plastik Bingham dengan adanya gradien tekanan linier pada arah aksial (sumbu x)

Aplikasi-aplikasi yang membutuhkan gaya besar akan menyebabkan FER cenderung berada pada keadaan post-yield sehingga diharapkan pemodelan fluida dalam Bingham-plastic dapat merepresentasikan karakter peredam kejut yang akurat. Profil kecepatan pada celah untuk fluida Bingham dapat digambarkan seperti pada gambar (3.6). Sedangkan sifat fluida Bingham dapat dinyatakan dalam persamaan: du  du  , τ = τ y sgn   + µ dr  dr 

τ > τy

(3.55a)

du =0, dr

τ < τy

(3.55b)

disini τy adalah yield stress dinamik dan diasumsikan sebagai fungsi polinom dari medan listrik. µ adalah viskositas plastis dan diasumsikan independent terhadap kuat medan listrik agar analisis kita menjadi lebih sederhana. Rpi dan Rpo pada gambar menyatakan jari-jari dalam dan luar plug. Selanjutnya akan diturunkan persamaan-persamaan pada ketiga daerah pada gambar diatas satu-persatu. Daerah 1 (post-yield): R1 < r < Rpi, τ > τ y , du dr > 0 , sehingga shear stress:

44

τ =τy + µ

du dr

(3.56)

yang jika diturunkan terhadap r menjadi seperti berikut: dτ d 2u =µ 2 dr dr

(3.57)

sehingga persamaan (3.45) menjadi seperti persamaan (3.58). µ

d 2 u τ y µ du ∆P + + = r r dr L dr 2

(3.58)

Persamaan (3.58) diintegralkan dan dipecahkan untuk mencari uBR1(r), yang dituliskan pada persamaan (3.59). u BR1 (r ) =

∆P 2 τ y r − r + C1 ln(r ) + D1 4µ L µ

(3.59)

dimana C1 dan D1 adalah konstanta yang dapat diperoleh dari syarat batas daerah 1, yaitu kecepatan alir pada dinding elektroda sama dengan nol, dan percepatan alir pada titik temu antara fluida cair dan padat juga sama dengan nol karena plug melaju dengan kecepatan konstan. Syarat batas tersebut yang digambarkan pada persamaan di bawah ini: u ( R1 ) = 0 , dan

du dr

=0 r = R pi

Bila syarat batas itu diterapkan untuk mencari profil kecepatan, maka profil kecepatan menjadi seperti pada persamaan (3.60) berikut ini: u BR1 (r ) =

 r  τ y   r  ∆P  2 2 2 r − R1 − 2 R pi ln  − r − R1 − R pi ln  4µ L   R1  µ   R1 

(3.60)

Daerah 2 (pre-yield): Rpi < r < Rpo, τ < τ y , du dr = 0 , Ini adalah daerah alir padatan (plug) dan karena du dr = 0 maka kecepatan pada daerah ini dituliskan oleh persamaan (3.61). uBR2 = konstan = up

(3.61)

Dari persamaan (3.45) yang diintegrasikan terhadap r dan dipecahkan untuk τ, diperolehlah fungsi terhadap r seperti berikut:

45

τ (r ) =

C ∆P r+ 2 2L r

(3.62)

dimana C2 adalah konstanta yang didapat dari syarat batas pada daerah ini. Pada batas-batas plug, shear stress sama dengan yield stress bahan fluida. Syarat-syarat batas tersebut dituliskan pada persamaan-persamaan di bawah ini: τ(Rpi) = τy,

τ(Rpo) = –τy,

dan

Bila syarat batas tersebut diterapkan pada persamaan (3.62) dan mengalikan masing-masing persamaan dengan Rpi dan Rpo diperolehlah persamaan berikut: τ y R pi =

∆P 2 R pi + C 2 2L

(3.63a)

τ y R po =

∆P 2 R po + C 2 2L

(3.63b)

kemudian {[persamaan (3.63a)] – [persamaan (3.63b)]} akan menghasilkan persamaan yield stress (3.64). τy =

∆P (R pi − R po ) 2L

(3.64)

Memecahkannya untuk tebal plug, δp = (Rpo – Rpi), maka diperolehlah persamaan di bawah ini: δ p = R po − R pi =

τy

(3.65a)

∆P 2 L

karena ∆P = − F A adalah negatif pada kasus kita ini, maka tebal plug menjadi: δp =

2 Lτ y

(3.65b)

∆P

Daerah 3 (post-yield): Rpo < r < R2 , τ > τ y , du dr < 0 , sehingga shear stress: τ = −τ y + µ Persamaan

(3.66)

du dr

dan

(3.66) (3.57)

dimasukkan

ke

persamaan

(3.45),

lalu

mengintegrasinya dan memecahkannya untuk mencari profil kecepatan, yang dituliskan pada persamaan (3.67).

46

u BR 3 (r ) =

∆P 2 τ y r + r + C 3 ln(r ) + D3 4µ L µ

(3.67)

C3 dan D3 juga adalah konstanta yang didapat dari syarat batas pada daerah ini. Pada perbatasan dengan padatan, dimana r = Rpo, fluida geser mengalir pada kecepatan plug konstan. Dan kecepatan pada batas dengan dinding luar harus sama dengan nol. Syarat-syarat batas tersebut dapat ditulis sebagai berikut: du dr

u ( R2 ) = 0 , dan

=0 r = R po

Bila syarat batas tersebut di atas diterapkan pada persamaan (3.67) akan didapatkan profil kecepatan: u BR 3 (r ) =

 r ∆P  2 2 2 r − R2 − 2 R po ln 4µ L   R2

 r  τ y   + r − R2 − R po ln  R2  µ 

  

(3.68)

Pada persamaan-persamaan yang didapat sejauh ini masih dijumpai dua peubah yang tak diketahui, Rpo dan Rpi, yang hanya dijumpai berhubungan dalam satu persamaan, yaitu persamaan (3.65). Persamaan kedua didapat dari plug yang mengalir dengan kecepatan konstan, up, yaitu seperti pada persamaan (3.69b).

atau:

u p = u BR1 ( R pi ) = u BR 3 ( R po )

(3.69a)

u BR1 ( R pi ) − u BR 3 ( R po ) = 0

(3.69b)

Persamaan (3.69b) dan (3.65a) digunakan untuk mendapatkan jari-jari dalam dan luar plug, Rpo dan Rpi. Begitu Rpo dan Rpi dihitung, profil kecepatan dapat diketahui dari persamaan (3.60), (3.68) dan (3.69a) sehingga didapatkan distribusi kecepatan pada tiap daerah. Fluks volume yang dipindahkan oleh piston, Qp, sama dengan fluks volume yang melalui celah elektroda, QBR, atau dalam persamaan: Qp = QBR

(3.70)

Fluks volume total yang melalui celah elektroda adalah sebagai berikut: QBR = QBR1 + QBR2 + QBR3

(3.71)

Dimana masing-masing fluks volume yang melalui tiap-tiap daerah dihitung dengan persamaan berikut:

47

QBR1 = ∫ u BR1 (r ) 2πr dr R pi

(3.72)

R1

QBR 2 = ∫

R po

u BR 2 (r ) 2πr dr

(3.73)

QBR 3 = ∫ u BR 3 (r ) 2πr dr

(3.74)

R pi R2

R po

Lalu memasukkan masing-masing kecepatannya, persamaan (3.60), (3.61) dan (3.68), ke dalam persamaan fluks volume (3.72), (3.73) dan (3.74), mengintegrasinya dan menyederhanakannya sehingga didapatkan persamaanpersamaan sebagai berikut: QBR1 =

 R pi π ∆P  2 2 2 2 4  3R pi − R1 R pi − R1 − 4 R pi ln 8µ L   R1

−

(

(

(

+

)

  

πτy   R pi 2 2 3 (R pi − R1 ) 7 R pi + R pi R1 − 2 R1 − 6 R pi ln 6 µ   R1

2 QBR 2 = u p π R po − R pi2

QBR 3 =

)(

)

  

(3.75)

)

(3.76)

 R po π ∆P  2 2 2 2 4  − 3R po + R2 R po − R2 + 4 R po ln 8µ L   R2

(

)(

)

  

πτy   R po 2 2 3 − (R po − R2 ) 7 R po + R po R2 − 2 R2 + 6 R po ln 6 µ   R2

(

)

  

(3.77)

Kemudian dengan fluks volume Qp = A v0, kecepatan piston dan peredaman viskos ekivalen diberikan oleh CBR didapat dengan: v0 =

QBR1 + QBR 2 + QBR 3 A

(3.78)

C BR = F v 0

(3.79)

Metodologi pemecahan persamaan-persamaan diatas adalah sebagai berikut: Bila diberikan dimensi peredam kejut, properti fluida, dan gaya, maka persamaan (3.69b) dan (3.65a) digunakan untuk mencari jari-jari dalam, Rpo, dan luar, Rpi, plug. Baru kemudian persamaan (3.75), (3.76) dan (3.77) digunakan untuk mengetahui fluks volume masing-masing daerah

dan

kemudian

48

dijumlahkan untuk mendapatkan fluks volume total yang melalui celah elektroda. Setelah itu kecepatan kepala piston dapat dihitung dari persamaan (3.78) dan peredaman viskos ekivalennya diketahui dari persamaan (3.79).

Catatan patut diberikan pada metodologi ini, yaitu bahwa metodologi ini hanya berlaku sesaat (instantaneous), pada saat dikenakan gaya dengan besar tertentu dan pada saat dikenakan medan listrik tertentu barulah kita dapat menghitung tebal plug. Sasaran akhirnya adalah untuk mengetahui besar redaman viskos ekivalennya. Metodologi pemecahan ini tidak dapat kita terapkan bila yang diinginkan adalah besar harga koefisien redaman viskos ekivalennya dengan cara memberikan tegangan tertentu (yield stress tertentu) terhadap input gaya acak. Untuk hal yang demikian ini haruslah dibuat suatu tabel koefisien redaman viskos sebagai hasil dari perhitungan pengenaan medan listrik pada saat diberikan masukan gaya tertentu. Pada input gaya tertentu bila dikehendaki peredam kejut yang memiliki koefisien redaman viskos ekivalen tertentu dan memiliki respon sistem seperti yang diinginkan dengan cara memberikan input medan listrik pada peredam kejut, maka kita akan bergantung pada cara look-up table untuk mendapatkan hasil.

3.4.3

PARAMETER-PARAMETER PEREDAM KEJUT FER

Seperti pada sub-bab 3.3.3, parameter-parameter geometri peredam kejut yang digunakan pada model ini adalah seperti berikut: R = ½ d2 = 21,25 mm (jari-jari penampang kepala piston) A = 1418,625 mm2 (luas penampang kepala piston) R1 = ½ d3 = 23,3 mm (jari-jari luar elektroda dalam) a

= 0,5 mm (jarak celah antar elektroda)

R2 = ½ d3 + 2a = 23,8 mm (jari-jari dalam elektroda luar) L

= 63,5 mm (panjang saluran)

µ0 = 0,035 Pa.s (viskositas absolut fluida)

49

Pada saat tidak dikenakan medan listrik, untuk selanjutnya disebut keadaan pasif, fluida berlaku sebagai fluida Newton dengan gaya F peredam kejut diberikan oleh persamaan (3.52), dengan konstanta peredaman CNR dinyatakan oleh persamaan (3.53) yang terdiri dari faktor viskositas absolut dan faktor geometri peredam kejut, seperti pada persamaan (3.54). Bila digunakan parameter-parameter tersebut, maka faktor geometri peredam kejut, ΓNR , dapat dihitung untuk kemudian memperoleh harga konstanta peredaman pasif CNR, yaitu 2901,8 N.s/m.

Pada saat aktif, dalam hal ini berarti saat dikenakan medan listrik, yield stress dinamik FER adalah fungsi polinom berorde lima dari medan listrik E sebagaimana diuraikan pada persamaan (3.40) dan tabel 3.1. Berikut dituliskan kembali fungsi polinom tersebut untuk medan listrik DC: τ y ( E ) = −3,686 E 5 + 38,702 E 4 − 126,666 E 3 + 237,912 E 2 − 82,047 E + 11,026 (3.80) dengan batasan E pada 0 – 4 kV/mm, di mana τy (0) = 11,026 Pa dan τy (4 kV/mm) = 1516,3 Pa.

Selanjutnya peredam kejut FER ini dapat pula digambarkan oleh elemen friksi yang di-parallel dengan peredam viskos seperti digambarkan persamaan berikut ini: F = µ 0 ΓNR x& + f b sgn( x& )

(3.81)

di mana ΓNR adalah koefisien geometri peredam kejut dan fb adalah gaya friksi yang disebabkan kekakuan fluida karena adanya medan listrik. Komponen gaya yang terakhir ini diperoleh dari keadaan pada persamaan (3.65) dengan ∆P = − F A yang menghasilkan persamaan berikut ini: fb =

2L A τ y (E) a

Persamaan (3.81) akan tampak seperti pada ilustrasi di bawah ini.

(3.82)

50

Bila disertakan batasan medan listrik seperti persamaan (3.80) di atas, dapat diketahui minimum dan maksimumnya, dapat dilihat pada persamaan (3.83) di bawah ini. x

F

CNR

Gambar 3.7 Model Bingham peredam kejut dengan FER f b ( E = 0) = f min = 3,97

N,

(3.83a)

f b ( E = 4 kV/mm) = f maks = 546,37

N.

(3.83b)

Dan µ 0 ΓNR dengan µ0 diasumsikan konstan adalah koefisien redaman viskos, CNR, seperti yang diuraikan pada (3.53). Persamaan (3.81) lebih jelasnya diuraikan pada persamaan dibawah ini: F = 2901,8 x& + fb sgn( x& )

(3.84)

Patut dicatat bahwa pada model ini bila mana pada suatu titik kecepatan piston sama dengan nol maka gaya pada elemen friksi sama dengan gaya yang dikenakan pada sistem model ini.

51

3.5 MODEL SUSPENSI SEMI AKTIF SEPEREMPAT KENDARAAN Sebelum dapat memecahkan masalah sintesis pengontrol sistem suspensi dibutuhkan perumusan masalah dalam bentuk persamaan matematis. Perumusan suspensi aktif telah diperkenalkan pada sub-bab 2.4.1. Model inilah yang menjadi acuan sintesis pengontrol sistem suspensi yang menggunakan peredam kejut elektrorheologi ini. Peredam kejut elektrorheologi dapat dimodelkan sebagai peredam kejut continuously variable (berubah secara kontinu) pada sistem suspensi semi-aktif.

x1

m1 k1

cmin + c(t) x2

m2 k2

Gambar 3.8

x0

Model suspensi semi aktif seperempat kendaraan dengan peredam kejut variabel

Pada gambar diatas, m1 dan m2 menyatakan massa kendaraan (sprung) dan massa assembly ban serta suspensi kendaraan (unsprung), k1 dan k2 adalah konstanta kekakuan pegas suspensi dan ban. Sedangkan cmin + c(t) merupakan konstanta redaman peredam kejut variabel, yang dapat bervariasi antara cmin dan cmin+cmax, dimana c max + c min > c min ≥ 0 ; sehingga c(t) memenuhi 0 ≤ c(t) ≤ c max . Persamaan gerak sistem menjadi sebagai berikut: m1 &x&1 = − k 1 ( x1 − x 2 ) − [c min + c(t )]( x&1 − x& 2 )

(3.85a)

m 2 &x&2 = −k 1 ( x 2 − x1 ) − [c min + c(t )]( x& 2 − x& 1 ) + k 2 .( x 0 − x 2 )

(3.85b)

Bila vektor keadaan dan gangguan dinyatakan seperti berikut:

52

 v1   x1 − x 2  v   x&  1  , v =  2 =  v 3   x 2 − x 0     &  v 4   x 2 

w = x& 0

dan

(3.86a)

dan karena tidak semua komponen keadaan dapat diukur, anggaplah dipunyai keluaran berupa percepatan massa sprung, defleksi ruang suspensi dan defleksi ban. Matriks vektor keluaran dapat dinyatakan sebagai berikut:  &x&1  y =  x1 − x 2   x 2 − x0 

(3.86b)

Sistem persamaan pada persamaan (3.85) dapat dituliskan sebagai persamaan ruang keadaan berikut: v& = Av + B[c(v 2 − v 4 )] + Gw

(3.87a)

y = Cv + D[c(v 2 − v 4 )]

(3.87b)

dengan matriks-matriksnya sebagai berikut:  0 − k m A= 1 1  0   k1 m 2

1 − c min m1 0 c min m 2

 0  − 1 m  1 B= ,  0     1 m2  − k1 m1 C =  1  0

0 0 0 − k 2 m2

−1  c min m1  ; 1   − c min m 2 

0 0 G =  , − 1   0 − c min m1 0 0

0 c min m1  0 0  ; 1 0 

− 1 m1  D =  0   0 

(3.87c)

53

Diagram blok sistem tersebut dapat dilihat pada gambar di bawah ini.

Gambar 3.9

Diagram blok model matematis sistem suspensi seperempat kendaraan

Bila

disertakan

koefisien

redaman

ekivalen

dari

peredam

kejut

elektrorheologi seperti yang tertera pada persamaan (3.42d), maka dapat diketahui bahwa: c min = 2917,3532 Ns/m, c(t ) =

(3.88a)

0,3603 (−3,6857 E 5 + 38,7015 E 4 − 126,6655 E 3 + 237,9117 E 2 v2 − v4 − 82,0472 E + 11,0263) sgn(v 2 − v 4 ) Ns/m (3.88b)

c (t) = 0,

jika (v2 – v4) = 0

(3.88c)

v0, kecepatan piston peredam kejut, dianggap sama dengan kecepatan defleksi ruang suspensi, (v2 – v4). Peredam kejut FER dapat pula digambarkan oleh model Bingham seperti pada gambar 3.7, sehingga sistem suspensi seperempat kendaraan akan tampak seperti gambar 3.10. Model tersebut menyerupai sistem suspensi semi aktif lainnya, di mana komponen aktuator gayanya dalam kasus ini digantikan oleh elemen gaya friksi Coulomb, yang merupakan fungsi dari yield stress fluida, yang adalah fungsi dari medan listrik. Persamaan gerak sistem pada gambar di atas adalah: m1 &x&1 = − k1 ( x1 − x 2 ) − c NR ( x&1 − x& 2 ) − f b sgn( x&1 − x& 2 )

(3.89a)

m2 &x&2 = − k1 ( x 2 − x1 ) − c NR ( x& 2 − x&1 ) − f b sgn( x& 2 − x&1 ) + k 2 ( x0 − x 2 )

(3.89b)

54

Model ini perlu diberi catatan penting, yaitu jika pada suatu titik kecepatan yang didapat piston sama dengan nol maka gaya yang diberikan oleh elemen friksi adalah sama dengan gaya yang dikenakan padanya. x1

m1 k1

cNR x2

m2 k2

Gambar 3.10

x0

Model suspensi semi aktif seperempat kendaraan dengan peredam kejut model Bingham

Vektor keadaan dapat diambil seperti yang tertera pada persamaan (3.86) dan bahwa cNR adalah koefisien redaman minimal, cmin, maka sistem ini dapat dituliskan dalam persamaan ruang keadaan seperti berikut: v& = A v + B [ f b sgn(v 2 − v 4 )] + G w

(3.90a)

dan karena tidak semua keadaan dapat diukur maka dapat diambil keluaran berupa percepatan massa sprung, defleksi ruang suspensi dan defleksi ban. Persamaan ruang keadaan keluaran yang demikian adalah y = C v + D [ f b sgn(v 2 − v 4 )]

(3.90b)

Di mana matriks-matriks ruang keadaan A, B, G, C dan D sama dengan seperti pada persamaan 3.87.

55

BAB IV PERANCANGAN SISTEM KONTROL OPTIMAL Penerapan kontrol optimal pada sistem suspensi semi aktif telah banyak dilakukan. Haæ dan Youn (1992)[2] meneliti penggunaan kontrol optimal dengan preview untuk suspensi semi aktif pada seperempat kendaraan. Sonny Setiadi (1996)[26] melakukan studi perancangan kontrol dengan metoda regulator optimal untuk sistem suspensi aktif setengah kendaraan. Penelitian lain dilakukan oleh Satriyo (1998) pada sistem suspensi semi aktif seperempat kendaraan. Tugas akhir ini hendak mempelajari dan merancang suatu sistem kontrol pada sistem suspensi yang menggunakan peredam kejut berfluida elektrorheologi. Teori kontrol optimal akan dipergunakan untuk menurunkan suatu aturan kontrol optimal (optimal control law). Aturan kontrol ini akan meminimalkan indeks performansi yang melibatkan kenyamanan berkendaraan (ride comfort), roadholding, dan ruang defleksi suspensi.

4.1 INDEKS PERFORMANSI YANG DIGUNAKAN Permasalahan utama dalam perancangan sistem kontrol optimal adalah bahwa sistem suspensi kendaraan hendak dikontrol sedemikian rupa sehingga memenuhi kriteria kenyamanan dan keamanan. Karena itu pemilihan fungsi indeks performansi menentukan karakter sistem kontrol. Sistem seperempat kendaraan yang hendak dikontrol telah diketahui pada bab terdahulu. Indeks Performansi pada sistem suspensi seperempat kendaraan ini haruslah melibatkan kenyamanan berkendara, keamanan sistem suspensi karena adanya keterbatasan gerak (ruang) defleksi suspensi, dan keamanan berkendara yang menghendaki ban tetap menempel pada permukaan jalan. Indeks Performansi (IP) terhadap model sistem suspensi semi aktif dipilih sebagai berikut.

56

T

I= Indeks

[

]

1 2 2 ρ1 (x1 − x 2 ) + ρ 2 (x 2 − x 0 ) + ρ 3 &x&12 dt ∫ 20

performansi

tersebut

melibatkan

(4.1)

variansi

dari

percepatan

kompartemen penumpang (badan kendaraan) &x&1 , ruang defleksi suspensi

(x1 − x2 ) dan

defleksi ban (x 2 − x0 ) . Konstanta pembobot, ρ1, ρ2 dan ρ3 akan

dipilih kemudian. Pemilihan konstanta-konstanta pembobot ini akan dilakukan dengan metoda trial and error.

Maksud dari pemilihan indeks performansi tersebut adalah bahwa dengan meminimalkan variansi percepatan badan kendaraan akan diperoleh peningkatan isolasi kompartemen penumpang terhadap gangguan permukaan jalan sehingga akan meningkatkan kenyamanan penumpang selama berkendaraan.

Dengan

mengurangi variansi defleksi ban, yang merupakan komponen waktu gaya ban terhadap jalan, akan menjaga ban tetap menempel dengan jalan sehingga akan meningkatkan stabilitas. Dan menjaga supaya ruang gerak suspensi kecil adalah penting untuk mengurangi kemungkinan suspensi terbentur bagian bawah badan kendaraan karena keterbatasan ruang gerak suspensi.

4.2 SISTEM KONTROL OPTIMAL PADA SUSPENSI SEMI AKTIF SEPEREMPAT

KENDARAAN

DENGAN

PEREDAM

KEJUT

VARIABEL Model sistem dapat dilihat pada persamaan gerak sistem suspensi, persamaan (3.85a,b) bersama matriks keadaan, v, dan gangguan, w, sedangkan bentuk ruang keadaan sistem telah dituliskan pada persamaan (3.87). Bila persamaan tersebut digunakan dalam IP (persamaan 4.1) maka akan diperoleh IP sebagai fungsi masukan dan variabel keadaan sistem seperti yang dituliskan pada persamaan (4.2). T

{

}

1 2 I = ∫ v T Qv + 2 v T N [(v 2 − v 4 )c ]+ R [(v 2 − v 4 )c ] dt 20 dengan matriks-matriks pembobot sebagai berikut:

(4.2)

57

 ρ 3 k12 + ρ1 m12  1  ρ 3 k1c min Q= 2 m1  0   − ρ 3 k1c min

ρ 3 k1c min 2 ρ 3 c min 0

0 0 ρ 2 m12

2 − ρ 3 c min

0

 k1    ρ 3  c min  N= 2 , m1  0    − c min 

R=

− ρ 3 k1c min   − ρ 3 c min2   0  2 ρ 3 c min 

ρ3 m12

Diasumsikan bahwa Q adalah matriks simetris dan definit non-negatif, sedang R adalah bilangan skalar positif. Mengingat ρ3 hadir dalam semua komponen matriks pembobot maka harga ρ3 diambil sama dengan satu sehingga dari sini hanya perlu dipilih ρ1 dan ρ2.

4.2.1

PERANCANGAN PENGONTROL OPTIMAL

Permasalahan kontrol yang perlu dipecahkan adalah untuk mencari aturan kontrol yang dapat mengatur perubahan faktor redaman c(t) yang adalah input kontrol. Faktor redaman tersebut sama dengan persamaan di bawah ini: c(t ) = f [ v(τ ), w(τ + t p ), t 0 ≤ τ ≤ t ] ,

t0 ≤ t ≤ T

(4.3)

Persamaan (4.3) harus diupayakan dapat meminimalkan kriteria I yang dituliskan

di

atas

terhadap

syarat

pertidaksamaan

koefisien

redaman

c min ≤ c ≤ c max . Perlu ditekankan bahwa persamaan syarat kontrol diatas adalah bersifat causal karena masukan kontrol pada sebarang waktu t hanya tergantung pada informasi yang tersedia pada saat itu, khususnya pada keadaan yang lalu, sekarang dan pada gangguan hingga saat t + tp. Selanjutnya untuk menyederhanakan perancangan disini tp = 0, tidak tersedia waktu preview. Diasumsikan lebih lanjut perubahan koefisien redaman terjadi secara seketika, independen terhadap waktu. Mengikuti persamaan (4.2) dan persamaan (4.3), gaya kontrol pada sistem suspensi ini diberikan oleh persamaan (4.4). u = c (v 2 − v 4 )

(4.4)

58

Sistem suspensi seperempat kendaraan dengan peredam kejut variabel ini dapat dianggap sebagai sistem suspensi aktif penuh dengan beberapa pembatasan. Solusi optimal sistem ini berkaitan dekat dengan solusi optimal sistem aktif penuh sehingga akan lebih mudah dimengerti bila dibahas dahulu pemecahan untuk sistem aktif. Sistem aktif penuh diberikan oleh persamaan berikut ini: v& = A v + B u + G w

(4.5)

dimana matriks-matriks A, B dan G sama dengan yang diberikan sub bab 3.5. Indeks performansi sistem suspensi aktif diberikan oleh persamaan di bawah ini: T

[

]

1 J = ∫ v T Q v + 2 v T N u + R u 2 dt 20

(4.6)

dengan Q, N dan R sama dengan yang tertera pada persamaan (4.2), yang mana cmin adalah koefisien redaman peredam kejut pasif. Jika pasangan (A, B) stabil yaitu bila matriks keterkontrolan n × n berikut memiliki rank sama dengan n, atau memiliki n vektor kolom yang independent. [B AB K A n −1B] dan jika matriks Q memenuhi keadaan sebagai berikut:  Q1 2   12  Q A  rank  =n   M  1 2 n −1  Q A  Pasangan matriks (A, Q) tersebut bila dapat teramati maka solusi kontrol, yang mampu meminimalkan J pada (4.6), diperoleh sama dengan persamaan sebagai berikut:

(

)

u 0 (t ) = − R −1 N T + B T P v (t )

(4.7)

dimana P adalah solusi definit positif dari persamaan aljabar Riccati P A n + A Tn P − P B R −1 B T P + Q n = 0

(4.8)

59

4.2.2

SISTEM KONTROL OPTIMAL

Persamaan hukum kontrol yang optimal telah didapat pada persamaan (4.7). Secara intuitive terlihat bahwa gaya kontrol pada persamaan (4.4) untuk sistem semi aktif haruslah dapat mengikuti gaya optimal untuk sistem aktif seperti yang tertera pada persamaan (4.7), tanpa syarat batas. Pemecahan yang mungkin dilakukan untuk memperoleh harga indeks performansi yang optimal adalah dengan melakukan pengaturan terhadap c agar dapat meminimalkan fungsi berikut: L(c) = [c (v 2 − v 4 ) − u 0 ] 2

(4.9)

dengan batasan c min ≤ c ≤ c max . Perubahan koefisien peredaman c kemudian diusulkan diatur menurut aturan-aturan berikut. (1) Jika c min <

u0 u0 ≤ c max , maka c = (v 2 − v 4 ) v2 − v4

(4.10a)

(2) Jika

u0 ≤ c min , v2 − v4

maka c = cmin

(4.10b)

(3) Jika

u0 > c max , v2 − v4

maka c = cmax

(4.10c)

Dan pada v2 – v4 = 0, harga c tidak berubah. Koefisien redaman pada tugas akhir ini telah dijabarkan oleh persamaan (3.89b), yaitu seperti berikut: c=

0,3603 (−3,6857 E 5 + 38,7015 E 4 − 126,6655 E 3 + 237,9117 E 2 v2 − v4 − 82,0472 E + 11,0263) sgn(v 2 − v 4 )

[ Ns/m ]

Persamaan terakhir ini akan membawa kita pada suatu pendekatan baru. Sistem suspensi seperempat kendaraan dengan peredam kejut model Bingham seperti pada Gambar 3.9 serupa dengan sistem suspensi aktif tetapi komponen aktuator gayanya digantikan oleh elemen friksi Coulomb, yang dalam kasus pada tugas akhir ini merupakan fungsi dari yield stress fluida. Yield stress fluida peredam kejut ini sendiri dianggap semata-mata adalah fungsi dari medan listrik.

60

Gaya pengontrol sistem suspensi semi aktif pada persamaan (4.4) dapat pula ditulis sebagai berikut: u = f b sgn(v 2 − v 4 )

(4.11)

dimana fb sama dengan yang diberikan pada persamaan (3.82). Persamaan (4.9) pun dapat pula dituliskan sebagai berikut: L( f b ) = [ f b sgn(v 2 − v 4 ) − u o ]

2

(4.12)

dengan fb dibatasi oleh karakter fluida peredam kejut f min ≤ f b ≤ f maks . Di mana f min = fb (Emin = 0 kV/mm) dan f max = fb (Emax = 4 kV/mm).

Sehingga input kontrol dapat diatur menurut aturan-aturan berikut: (1) Jika (v 2 − v 4 ) > 0 è sgn(v 2 − v 4 ) = 1 maka f b = u 0

dengan f min ≤ f b ≤ f max

(4.32)

(2) Jika (v 2 − v 4 ) < 0 è sgn(v 2 − v 4 ) = −1 maka f b = − u 0

dengan f min ≤ f b ≤ f max

(3) Dan jika (v 2 − v 4 ) = 0 , maka fb = 0 untuk setiap f b ∈ [ f min , f max ] sehingga persamaan dinamika menjadi seperti berikut: v& = A v + G w y =C v

(4.33) (4.34)

61

BAB V SIMULASI MODEL FER DAN SISTEM SUSPENSI SEMI AKTIF MODEL SEPEREMPAT KENDARAAN Listing program simulasi pada bab ini dapat dilihat pada Lampiran A. 5.1 SIMULASI MODEL PEREDAM KEJUT FER Pertama-tama penulis melakukan simulasi model peredam kejut yang dikembangkan pada Bab III Pemodelan Sistem, untuk melihat karakteristik FER pada peredam kejut tersebut. Metodologi yang digunakan pada simulasi ini adalah hasil yang dikembangkan pada sub bab 3.4 (pendekatan geometri 1 dimensi Cylindrical Axisymmetric). Gambar 5.1 berikut ini menampilkan profil ketebalan plug FER pada celah peredam kejut dengan beban gaya yang bertambah. Pada pemberian medan listrik konstan, penambahan beban gaya akan menyebabkan berkurangnya ketebalan

Gambar 5.1 Profil ketebalan plug berkurang ketika gaya bertambah untuk medan listrik yang konstan: – – 1kV/mm; –. 2 kV/mm; -- 3 kV/mm; – 4 kV/mm.

62

Gambar 5.2 Profil ketebalan plug terhadap perubahan medan listrik ketika diberikan gaya konstan: – – 100 N; –. 500 N; -- 1000 N; – 1500 N; .. 2000 N.

plug. Hal ini mengindikasikan bahwa koefisien redaman viskos ekivalen peredam kejut juga berkurang. Dapat pula dikatakan dengan cara lain, yaitu bahwa pada gaya konstan, penambahan besar medan listrik akan menyebabkan koefisien redaman bertambah pula. Usaha memperbesar koefisien redaman peredam kejut ini dibatasi oleh kekuatan dielektrik bahan FER yang digunakan, pada TXER6 medan listrik maksimal yang dapat diberikan adalah 4 kV/mm. Peredam kejut hasil pemodelan mampu memberikan gaya maksimum sebesar 546,37 N pada medan listrik 4 kV/mm (dapat dilihat kembali persamaan (3.83)). Gambar-gambar selanjutnya pada sub bab ini akan mengilustrasikan tentang profil kecepatan terhadap pemberian medan listrik yang bertambah dari 0 – 4 kV/mm pada celah elektroda peredam kejut untuk mengantisipasi beberapa tingkat beban gaya yang diberikan. Pada aplikasi gaya 500 N (gambar 5.3), kecepatan alir FER pada celah elektroda pada 4 kV/mm sama dengan nol, yang berarti fluida tidak mengalir dan berlaku masih seperti padatan. Pada aplikasi gaya 600 N (gambar 5.4), kecepatan fluida pada 4 kV/mm (kurva paling kiri, dekat sumbu y) sedikit lebih besar dari nol dan mayoritas fluida berada pada keadaan padat dengan kecepatan konstan.

63

Hal tersebut menandakan bahwa beda tekanan pada peredam yang disebabkan oleh gaya 600 N sudah melampaui (lebih besar dari) batas yield stress maksimal FER pada 4 kV/mm untuk peredam kejut ini.

Gambar 5.3 Profil kecepatan pada celah elektroda dengan adanya gaya 500 N. Kurva menunjukkan profil kecepatan pada saat dikenakan medan listrik 0 kV/mm; 1 kV/mm; 2 kV/mm; 3 kV/mm; 3,5 kV/mm, berturut-turut dari kanan ke kiri.

Gambar 5.4 Profil kecepatan pada celah elektroda dengan adanya gaya 600 N. Kurva menunjukkan profil kecepatan pada saat dikenakan medan listrik 0 kV/mm; 1 kV/mm; 2 kV/mm; 3 kV/mm; 3,5 kV/mm; 4 kV/mm berturut-turut dari kanan ke kiri.

64

Gambar 5.5 Profil kecepatan pada celah elektroda dengan adanya gaya 1000 N. Kurva menunjukkan profil kecepatan pada saat dikenakan medan listrik 4 kV/mm; 3,5 kV/mm; 3 kV/mm; 2 kV/mm; 1 kV/mm; 0 kV/mm berturut-turut dari kiri ke kanan.

Pada gambar 5.5, profil kecepatan pada 1 kV/mm dan 0 kV/mm didapatkan hampir sekurva, yang berarti yield stress pada 1 kV/mm mendekati stress karena gaya beban sebesar 1000 N itu. Radius sebesar 0,0233 m pada ketiga gambar di atas adalah jari-jari luar elektroda dalam, sedangkan 0,0238 m adalah jari-jari dalam elektroda luar peredam kejut.

Selanjutnya petugas akhir melakukan simulasi untuk sistem peredam kejut seperti pada persamaan (3.81). Simulasi sistem dilakukan terhadap eksitasi masukan sinusoidal 5 Hz dengan amplituda 2 cm (4 cm peak-to-peak). Hasil simulasi dapat diamati pada gambar 5.6.

65

a) Gaya terhadap waktu

b) Gaya terhadap kecepatan

4 3 0

Gambar 5.6

Gaya peredam terhadap eksitasi sinusoida 5 Hz dengan amplituda sebesar 2 cm

66

c) Gaya terhadap perubahan posisi

4 3 0

Gambar 5.6

Gaya peredam terhadap eksitasi sinusoida 5 Hz dengan amplituda sebesar 2 cm (...lanjutan)

5.2 PEMILIHAN PARAMETER MODEL SEPEREMPAT KENDARAAN Model yang hendak diamati perlu ditentukan secara definitif supaya dapat dilanjutkan ke simulasi sistem kontrol. Batasan pemilihan pada koefisien peredaman minimal c

min

= 2901,8 Ns/m. Maka penentuan perbandingan cmin/m1,

k1/m1 dan k2/m1 sangat menentukan karakteristik model karena akan menentukan massa sprung yang kemudian akan menentukan harga konstanta-konstanta pegas, yang berarti akan menentukan harga frekuensi pribadi model kendaraan. Pemilihan cmin/m1 dimulai dari konfigurasi model yang didapatkan Saudara Raymond Mudrig, yaitu menggunakan perbandingan parameter-parameter m1, m2, k1, k2 dan cmin. cmin/m1 = 1,64 sehingga m1 = 1779 kg; k1/m1 = 245 sehingga k1 = 436.237 N/m; k2/m1= 2016 sehingga k2 = 3.586.896 N/m, dan kemudian m1/m2= 0.12 sehingga m2=213 kg. Konfigurasi tersebut memiliki respon frekuensi seperti

67

terlihat pada gambar 5.7 berikut. Pada gambar itu terlihat dua buah puncak yang menandakan letak frekuensi pribadi sistem, yaitu masing-masing untuk massa sprung dan massa unsprung.

Gambar 5.7 Respon frekuensi model Raymon

Frekuensi pribadi massa sprung didapatkan pada 14,9 rad/sec, setara dengan 2,4 Hz, dan 127 rad/sec untuk massa unsprung, atau sama dengan 20,2 Hz. Namun karena harga-harga konstanta pegas yang digunakan pada model ini sangat besar, yang berarti pegas yang dipergunakan adalah sangat kaku, maka patut dicoba harga perbandingan lain yang memberikan harga konstanta pegas lebih kecil. Tetapi mengubah harga perbandingan konstanta pegas berarti menggeser posisi frekuensi pribadi yang terjadi, sedangkan mengubah perbandingan cmin/m1 akan mengubah bentuk puncak pada plot respon frekuensi. Pada tabel 5.1 dapat dilihat beberapa konfigurasi perbandingan parameter lup terbuka sistem. Penentuan frekuensi pribadi dilakukan dengan menerapkan pendekatan persamaan di bawah ini: k eq =

k1 k 2  m  k 2 + k1 1 + 2  m1  

(5.1a)

68

fs =

1 2π

k eq

fu =

1 2π

k1 + k 2 m2

m1

untuk massa sprung, dan

(5.1b)

untuk massa unsprung.

(5.1c)

Tanggapan frekuensi sistem suspensi pasif didapat bila input gaya u sama dengan nol. Tanggapan frekuensi masing-masing konfigurasi tersebut dapat diamati pada gambar 5.8. Pada konfigurasi yang mempunyai harga konstanta pegas suspensi k1 ≈ 10.000 N/m, atau lebih kecil dari itu, hanya dijumpai satu puncak pada frekuensi di sekitar 18 rad/sec ( ≈ 3 Hz) yang menandakan bahwa pegas suspensi tersebut relatif lunak terhadap massa sprung dan unsprung. Pemilihan konfigurasi ini menginginkan didapatkannya dua puncak yang masing-masing menandakan frekuensi pribadi untuk massa sprung dan massa unsprung beserta harga-harga parameter massa dan konstanta pegas yang mendekati harga-harga pada literatur referensi. Pada dasarnya konfigurasi model dapat dipilih sebarang, namun agar dapat mendekati model sebenarnya perlulah dicari harga parameter yang sesuai.

Tabel 5.1 Konfigurasi parameter model seperempat kendaraan No. Model

A

B

C

D

E

F

cmin/m1

2

10

20

2

10

20

M2/m1

0,125

0,125

0,125

0,125

0,125

0,125

k1/m1

200

200

200

30

30

30

k2/m1

2000

2000

2000

360

360

360

M1 (kg)

1451

290

145

1451

290

145

M2 (kg)

181

36

18

181

36

18

k1 (N/m)

290187

58037

29019

43528

8706

4352

k2 (N/m)

2901868

580374

290187

522336

104467

52234

ωn

13,4 [2,1]

13,4

13,4

5,3 [0,8]

5,3

5,3

132,7 [21,1]

132,7

132,7

55,9 [8,9]

55,9

55,9

sprung (rad/sec) [fs Hz]

ω n unsprung(rad/sec) [fu Hz]

69

Tabel 5.1 Konfigurasi parameter model seperempat kendaraan (…lanjutan) No.

G

H

I

J

K

L

cmin/m1

2

10

20

6

6

6

m2/m1

0,125

0,125

0,125

0,125

0,125

0,125

k1/m1

60

60

60

200

60

30

K2/m1

500

500

500

2000

500

360

m1

1451

290

145

484

484

484

m2

181

36

18

60,5

60,5

60,5

k1

87056

17411

8706

120911

29019

14509

k2

725467

145093

72547

967289

241822

174112

ω n [fs] sp

7,3 [1,2]

7,3 [1,2]

7,27 [1,2]

13,4 [2,1]

7,3 [1,2]

5,3 [0,8]

ω n [fu]unsp

66,9 [10,7]

66,9 [10,7]

66,9 [10,7]

132,7 [21,1]

66,9 [10,7]

55,9 [8,9]

Pada massa sprung, frekuensi pribadi yang lazimnya dijumpai adalah pada 1 – 3 Hz; dan frekuensi pribadi pada massa unsprung berada di antara 10 – 20 Hz. Sedangkan konstanta pegas suspensi konvensional umumnya berada di sekitar harga 20 kN/m dan konstanta pegas ban kurang lebih 200 kN/m. Dengan massa sprung seperempat kendaraan antara 200 hingga 500 kg dan massa unsprung antara 40 – 60 kg. Massa unsprung yang lebih tinggi membutuhkan konstanta peredaman yang lebih tinggi pula. Dari kedua belas konfigurasi parameter diatas, yang memiliki massa dan konstanta pegas mendekati model umum seperempat kendaraan diatas adalah konfigurasi H dan K. Frekuensi pribadi kedua model tersebut adalah 1,2 Hz (7,3 rad/sec) untuk massa sprung dan 10,7 Hz (66,9 rad/sec) untuk massa unsprung. Penulis memilih konfigurasi model H karena harga konstanta pegas dan massanya lebih mendekati referensi ( m1 ≈ 250 kg; m2 ≈ 50 kg; k 2 ≈ 120 000 N/m ) dari pada K.

70

Gambar 5.8

Respon frekuensi defleksi suspensi sistem dengan konfigurasi seperti pada tabel 5.1

71

Gambar 5.8

Respon frekuensi defleksi suspensi sistem dengan konfigurasi seperti pada tabel 5.1 (…lanjutan)

Pada tabel 5.2 berikut diuraikan sekali lagi parameter model awal untuk simulasi hasil rancangan sistem kontrol.

72

Tabel 5.2 Parameter sistem suspensi seperempat kendaraan awal No. Model

H

cmin/m1 (Perbandingan konstanta redaman terhadap massa sprung/ ¼ kendaraan)

10

m2/m1 (Perbandingan massa unsprung terhadap massa sprung)

0,125

k1/m1 (Perbandingan kekakuan pegas suspensi terhadap massa sprung)

60

k2/m1 (Perbandingan kekakuan pegas ban terhadap massa sprung)

500

Massa sprung ( ¼ kendaraan), m1 (kg)

290

Massa unsprung (ban, pelek, poros), m2 (kg)

36

Konstanta kekakuan pegas suspensi, k1 (N/m)

17411

Konstanta kekakuan pegas ban, k2 (N/m)

145093

Frekuensi pribadi massa sprung, ω n rad/sec [fn Hz]

7,3 [1,2]

Frekuensi pribadi massa unsprung,

ω u rad/sec [fu

Hz]

66,9 [10,7]

Konfigurasi awal H ini selanjutnya digunakan pada carmodel.m, kemudian tiap keadaan, masukan dan keluaran dapat dinyatakan dengan perintah sebagai berikut: states={'su_def' 'sp_vel' 'us_def' 'us_vel'}; inputs={'er_force' 'road_vel'}; outputs={'sp_acc' 'su_def' 'us_def'}; a = su_def

sp_vel

us_def

us_vel

su_def

0

1.00000

0

-1.00000

sp_vel

-60.00000

-10.00000

0

10.00000

us_def

0

0

0

1.00000

us_vel

480.00000

80.00000 -4000.00000

-80.00000

b = er_force

road_vel

su_def

0

0

sp_vel

-0.00345

0

us_def

0

-1.00000

us_vel

0.02757

0

su_def

sp_vel

us_def

us_vel

sp_acc

-60.00000

-10.00000

0

10.00000

su_def

1.00000

0

0

0

us_def

0

0

1.00000

0

c =

73

d = er_force

road_vel

sp_acc

0.00345

0

su_def

0

0

us_def

0

0

Continuous-time system.

Dua masukan pada sistem kontrol ini adalah gaya kontrol dari peredam kejut FER (er_force) dan masukan gangguan berupa kecepatan vertikal permukaan jalan (road_vel).

5.3 PEMILIHAN PARAMETER PEMBOBOT Matriks-matriks pembobot yang digunakan adalah hasil dari indeks performansi pada persamaan (4.2): 3600 + ρ 1  600 Q =  0   − 600

600 0 100 0 0 ρ2 − 100 0

− 600 − 100  , 0   100 

 0,2068   0,0345  , N=  0    − 0,0345

R = 1,188 × 10 −5

Q dan R adalah matriks-matriks simetris dan definit non-negatif. Aturanaturan kontrol yang diperoleh pada bab terdahulu, dan juga performansi lup tertutup sistem, sebenarnya tergantung pada pemilihan konstanta pembobot ρ1 dan ρ2 yang tergantung pada perancang sistem itu. Harga ρ1 yang terlalu kecil akan menyebabkan sistem over-damped dan mempunyai tanggapan sistem yang lambat. Sedangkan memperbesar ρ1 tampak akan menyebabkan sistem berosilasi (damping ratio menjadi lebih kecil), yang secara nalar berarti mengurangi kenyamanan berkendara, tetap hal ini terlihat meningkatkan (mempercepat) tanggapan sistem.

74

Tanggapan sistem terhadap masukan step dapat diamati pada gambar 5.9, ρ1 diambil terhadap ρ2 konstan (=10000). Berikut ini karakter dari tanggapan step sistem. POS: persentase overshoot; Tr: waktu naik (rise time); Ts: waktu penetapan (settling time) 2%; dan PeakTime: waktu puncak.

Tabel 5.3

“Step response data” untuk perubahan konstanta pembobot

ρ1

POS

Tr

Ts

PeakTime

1 10 100 1000 10 000 1000 000

0 0,14424 5.0098 7.6134 126.3723 872.0914

1,82 1.34 0.665 0.215 0.005 0

2,415 1.865 1.775 0.92 0.515 1.12

2,5 2.34 1.22 0.52 0.02 0.005

ρ2 1 10 100 1000 10 000 1000 000

POS 30.4315 30.3343 29.3943 22.4205 7.6134 5.198

Tr 0.02 0.02 0.02 0.025 0.215 0.53

Ts 1.77 1.77 1.67 1.27 0.92 1.45

PeakTime 0.42 0.42 0.42 0.42 0.52 0.99

Pada perubahan ρ2 dijumpai bahwa pada harga yang relatif kecil atau sama dengan besar ρ1, sistem berosilasi dan mempunyai jejak yang hampir sama. Namun pada harga ρ2 yang relatif lebih besar terhadap ρ1 jejak tanggapan menunjukkan peningkatan rasio redaman. Tabel di atas menunjukkan tanggapan sistem dengan perubahan ρ2 terhadap ρ1 konstan (=1000). Parameter ρ1 sendiri adalah pembobot pada defleksi suspensi, dalam simulasi sistem akan digunakan ρ1 sama dengan 1000, dan ρ2, pembobot untuk defleksi ban, sama dengan 10000. Sedangkan ρ3 pada persamaan 4.2, yang adalah pembobot percepatan massa sprung, diambil sama dengan satu.

75

1000

1

Gambar 5.9 Respon step percepatan badan kendaraan dengan ρ1 yang berbeda (1, 10, 100, 1.000,10.000, 1.000.000)

5.4 RESPON FREKUENSI SISTEM PASIF DAN AKTIF PENUH Pada sub bab ini, akan dibandingkan performansi sistem suspensi pasif terhadap sistem suspensi aktif. Respon sistem pasif yang dimaksud adalah respon sistem tanpa masukan kontrol u. Sedangkan sistem aktif penuh disini adalah model sistem yang menggunakan umpan balik keadaan dengan gaya kendali pada peredam semi-aktif aktif tanpa batasan tertentu. Respon frekuensi sistem pasif dan aktif penuh tersebut dapat dilihat pada gambar 5.10, masing-masing untuk percepatan vertikal massa unsprung, defleksi suspensi dan defleksi ban. Garis kontinu menandakan tanggapan untuk sistem pasif dan garis putus-putus mewakili sistem suspensi aktif.

76

Gambar 5.10 Respon frekuensi sistem pasif (garis kontinu) dan sistem aktif (garis putus-putus)

77

Gambar 5.10 Respon frekuensi sistem pasif (garis kontinu) dan sistem aktif (garis putus-putus) (…lanjutan) Terlihat dari gambar bahwa sistem lup tertutup suspensi aktif secara umum dapat menekan penguatan (gain) pada percepatan vertikal massa kendaraan, yang berarti terjadi perbaikan kriteria kenyamanan. Sistem aktif melakukan perbaikan tersebut dengan kompensasi penurunan pada kriteria keamanan seperti tampak pada respon defleksi suspensi dan defleksi ban.

5.5 RESPON SISTEM TERHADAP MASUKAN JALAN BUMP Masukan jalan bump yang diberikan menggunakan bentuk sinusoida setengah perioda seperti pada persamaan 5.2.   vπ  t ,  Y sin  y=  d   0,

untuk t0 < t < t0 + 0,5T

(5.2)

untuk nilai t yang lain

Pada gambar 5.11 digunakan amplituda Y sama dengan 0.04 m dengan frekuensi 13.89π (43.64) rad/sec (sama dengan 6.95 Hz), yang adalah representasi

78

dari jarak jalan bump d sama dengan 1.2 m untuk kecepatan kendaraan v sama dengan 16.67 m/s untuk durasi simulasi 1,5 detik.

Gambar 5.11 Masukan jalan “bump” dengan amplituda 0.04 m Tanggapan sistem terhadap masukan jalan tersebut selanjutnya ditunjukkan pada Gambar 5.12.

Garis putus-putus panjang menunjukkan tanggapan dari

sistem aktif, sedangkan sistem pasif oleh garis kontinu. Sebagaimana konsisten dengan hasil tanggapan frekuensi pada sub bab terdahulu, sistem aktif memberikan hasil yang lebih baik dari pada sistem pasif pada faktor kenyamanan. Dari jejak tanggapan waktu bump ini didapatkan harga rms percepatan vertikal massa kendaraan, defleksi ruang suspensi dan defleksi ban pada sistem pasif masing-masing sama dengan 3,9455 m/s2, 0,01197 m, dan 0,01470 m, sedangkan pada sistem aktif masing-masing keluaran sistem tersebut menjadi 1,6178 m/s2, 0,02257 m, dan 0,02078 m. Berarti sistem aktif memberikan perbaikan pada kenyamanan sebesar 58,99%, dengan kompensasi performansi ruang suspensi menurun sebesar 88,49%, dan keamanan berkendaraan berkurang hingga 41,42%.

79

(a)

(b)

Gambar 5.12 Perbandingan tanggapan waktu sistem pasif (garis kontinu), sistem aktif penuh (garis putus-putus) dan sistem suspensi ER (titik-titik): (a) Percepatan vertikal massa kendaraan; (b) Defleksi suspensi

80

(c)

Gambar 5.12 Perbandingan tanggapan waktu sistem pasif (garis kontinu), sistem aktif penuh (garis putus-putus) dan sistem suspensi ER (titik-titik): (c) Defleksi ban Tingkat perbaikan yang didapatkan adalah hasil perbaikan bila sistem aktif ini memiliki aktuator gaya yang dapat memberikan masukan gaya pengontrol sesuai dengan permintaan pengontrol. Dengan kata lain hal ini dapat dicapai bila tidak ada batasan kemampuan gaya yang diberikan oleh aktuator. Sejak awal, tugas akhir ini menggunakan “aktuator gaya” yang sebenarnya adalah peredam kejut sistem pasif itu sendiri. Aktuator gaya yang digunakan adalah peredam kejut yang menggunakan fluida elektrorheologi, yang dari hasil pemodelan akan memiliki gaya tambahan maksimal ketika diberikan medan listrik 4 kV/mm, yaitu gaya kontrol maksimal sebesar 546,37 N. Pada gambar 5.12 di atas, respon sistem suspensi FER dapat diamati dari jejak berupa garis titik-titik. Gambar 5.13 di halaman berikut menggambarkan kemampuan peredam kejut FER mengikuti perintah sistem kontrol. Sinyal kontrol aktif terpotong pada sistem FER karena adanya keterbatasan gaya maksimal seperti telah diuraikan pada paragraf sebelum ini.

81

Gambar 5.13 Gaya dari Peredam FER (garis kontinu) dan gaya kontrol aktif (garis putus-putus) Harga-harga rms untuk masing-masing sistem beserta perbaikan sistem terhadap sistem pasif dapat dilihat pada tabel berikut ini: Tabel 5.3

Harga RMS dari respon sistem terhadap masukan jalan bump Percepatan Vertical (m/s2)

Defleksi Suspensi (m)

Defleksi ban (m)

Perbaikan (%)

Sistem Pasif

3,9455

0,01197

0,01470

-

Sistem Aktif Penuh

1,6178

0,02257

0,02078

[58,99 -88,49 -41,42]

Sistem FER

3,6629

0,01364

0,01543

[7,16 -13,94 -4,98]

Gambar 5.12 dan gambar 5.13 memperlihatkan bahwa respon waktu defleksi suspensi dan defleksi ban sistem aktif dan sistem FER memiliki harga yang lebih tinggi dari pada sistem pasif.

Sistem FER memberikan respon terbaiknya

terhadap masukan sinyal kontrol yang harus “diikutinya”. Hal ini berdampak pada lebih kecilnya perbaikan pada kriteria kenyamanan dibandingkan dengan sistem aktif, namun kompensasi yang lebih baik diberikan sistem FER pada kriteria keamanan suspensi dan keamaan berkendaraan.

82

BAB VI USAHA PERBAIKAN SISTEM KONTROL 6.1 INDEKS PERFORMANSI BARU Dalam Bab IV Perancangan Sistem Kontrol, indeks performansi yang dipilih tidak menyertakan usaha mengoptimalkan masukan gaya pengendali u pada sistem suspensi. Gaya yang diperoleh oleh pengontrol sudah kita lihat pada sub bab terdahulu, dimana perintah masukan gaya maksimal berada jauh di atas kemampuan ‘aktuator’ yang ada pada permasalahan ini. Selanjutnya Indeks Performansi pada persamaan (6.1) akan melibatkan masukan gaya kendali selain dari percepatan massa sprung (badan kendaraan), defleksi ruang suspensi dan defleksi ban. T

[

]

1 2 2 I = ∫ ρ 1 (x1 − x 2 ) + ρ 2 (x 2 − x0 ) + ρ 3 &x&12 + ρ 4 u 2 dt 20

(6.1)

IP ini dapat dinyatakan sebagai fungsi dari vektor keluaran y dalam bentuk matriks sebagai berikut: T

I=

[

]

1 y T Q y y + R u 2 dt 2 ∫0

(6.2)

Matriks pembobot Qy dan R diuraikan seperti di bawah ini: ρ3 Q y =  0  0

0 ρ1 0

 ; R=ρ ; 4  ρ 2  0 0

Konstanta-konstanta pembobot, ρ1, ρ2 dan ρ3 pada Qy diambil sama dengan konstanta pembobot pada bab V, sedangkan ρ4 sama dengan besar (ρ3/m12). Bentuk IP di atas masih perlu distandarisasi menjadi seperti persamaan (4.6). Standarisasi ini dilakukan terhadap matriks-matriks pembobot baru Qy1, Ry, dan Ny sebagai berikut: Q y1 = C T Q y C ; R y = R + D T Q y D ; N y = C T Q y D

(6.3)

83

Hasil-hasil dari pengontrolan dengan IP baru ini dapat dilihat pada sub bab berikut ini. Selanjutnya sistem kontrol yang menggunakan IP seperti pada persamaan (6.2) akan disebut dengan sistem umpan balik keluaran, sedangkan sistem pada bab V akan disebut sistem dengan umpan balik keadaan penuh.

6.2 RESPON SISTEM TERHADAP MASUKAN JALAN BUMP Gambar berikut menjelaskan respon waktu terhadap masukan jalan bump seperti pada gambar 5.11 dari gaya aktif pada sistem FER dan gaya yang dikehendaki oleh pengontrolan aktif penuh yang menggunakan matriks-matriks pembobot baru di atas.

Gambar 6.1

Gaya dari peredam FER (garis kontinu) dan Gaya aktif umpan balik keluaran (garis putus-putus)

Tampak pada gambar di atas bahwa gaya maksimal yang diminta oleh sistem kontrol aktif dengan IP baru tersebut lebih kecil dari pada gaya aktif maksimal pada bab sebelumnya. Pada gambar 5.13, gaya maksimal sinyal kontrol adalah 2930,4 N, sedangkan pada gambar 6.1 harga maksimal sinyal kontrol yang diminta oleh sistem aktif sama dengan 1055,2 N. Dengan demikian persentase rasio kemampuan suspensi FER mengikuti sinyal kontrol aktif seperti pada

84

gambar 6.1 lebih besar, yaitu 51,78%, dibandingkan dengan sinyal kontrol pada bab V, yang hanya sebesar 18,64%. a.)

b.)

Gambar 6.2

Perbandingan tanggapan waktu terhadap masukan jalan bump: (a) Percepatan vertikal massa sprung; (b) Defleksi suspensi

85

c.)

Gambar 6.2

Perbandingan tanggapan waktu terhadap masukan jalan bump: (c) Defleksi ban

Pada gambar 6.2 di atas garis kontinu mewakili suspensi pasif, garis putusputus mewakili sistem suspensi aktif, dan garis titik-titik mewakili sistem suspensi FER. Tabel berikut menghitung harga rms tanggapan sistem yang diperoleh di atas.

Tabel 6.1

Harga RMS dari tanggapan sistem terhadap masukan jalan bump Percepatan Vertical (m/s2)


Defleksi ban (m)

Perbaikan (%)

Sistem Pasif

3,9455

0,01197

0,01470

-

Sistem Aktif Penuh

3,6250

0,01441

0,01589

[8,14 -20,33 -8,15]

Sistem FER

3,6538

0,01349

0,01539

[7,39 -12,68 -4,75]

86

Secara keseluruhan sistem suspensi FER akan mengadaptasi perintah kontrol sebaik mungkin sehingga akan memberikan hasil sistem yang performansinya tidak seperti sistem aktif penuh bila gaya yang diminta pengontrol melewati batas kemampuan gaya peredam kejut FER. Dalam usaha pencapaian sistem yang kenyamanan, keamanan suspensi dan keamanan berkendaraannya sebaik mungkin, lazimnya besar sinyal kontrol tidak menjadi suatu faktor yang penting. Pencapaian kriteria IP oleh sistem kontrol yang diperoleh pada bab V masih lebih baik dari pada sistem pada sub bab ini, karena pada bab tersebut tidak dipedulikan berapa besar sinyal kontrol yang harus diberikan oleh sistem aktif penuh. Kompromi harus dijalankan antara pengoptimalan sinyal kontrol, faktor kenyamanan, faktor keamanan suspensi, dan keamanan berkendaraan. Gambar-gambar berikut ini akan membandingkan hasil simulasi sistem kontrol dengan umpan balik keadaan dan sistem kontrol dengan umpan balik keluaran yang didapatkan pada gangguan jalan bump setinggi 4 cm dan lebar 1,2 m. Hasil ini didapatkan bila kecepatan kendaraan 20 km/jam.

Gambar 6.3

Gaya dari peredam FER (garis kontinu) dan gaya aktif umpan balik keadaan penuh (garis putus-putus)

87

a.)

b.)

Gambar 6.4

Perbandingan tanggapan waktu sistem kontrol umpan balik keadaan penuh terhadap masukan jalan bump: (a) Percepatan vertikal massa sprung; (b) Defleksi suspensi

88

c.)

Gambar 6.4

Perbandingan tanggapan waktu sistem kontrol umpan balik keadaan penuh terhadap masukan jalan bump: (c) Defleksi ban

Pada gambar 6.4, garis kontinu mewakili sistem pasif, sistem aktif penuh dengan umpan balik keadaan penuh digambarkan oleh garis putus-putus, dan garis titik-titik mewakili tanggapan sistem FER. Gambar 6.5 di bawah ini adalah tanggapan sistem pasif (garis kontinu), aktif dengan umpan balik keluaran (garis putus-putus), dan tanggapan sistem FER (garis titik-titik) terhadap sinyal aktif umpan balik keluaran. Pada gambar 6.5 tampak hanya dua jejak sistem. Hal ini terjadi karena sistem aktif dengan umpan balik keluaran memberikan permintaan sinyal kontrol pada peredam kejut FER yang dapat diberikan oleh peredam kejut FER ini seperti yang digambarkan pada gambar 6.6. Harga RMS dari tanggapan sistem-sistem tersebut dapat diamati pada tabel 6.2. Pada tabel ini terlihat bahwa harga RMS untuk sistem aktif dengan umpan balik keluaran sama dengan harga RMS untuk sistem FER yang menerima sinyal kontrol dari sistem aktif tersebut.

89

a.)

b.)

Gambar 6.5

Perbandingan tanggapan waktu sistem kontrol umpan balik keluaran terhadap masukan jalan bump: (a) Percepatan vertikal massa sprung; (b) Defleksi suspensi

90

c.)

Gambar 6.5

Gambar 6.6

Perbandingan tanggapan waktu sistem kontrol umpan balik keluaran terhadap masukan jalan bump: (c) Defleksi ban

Gaya dari peredam FER dan gaya aktif umpan balik keluaran

91

Tabel 6.2 Harga RMS dari tanggapan sistem dengan kecepatan kendaraan 20 km/jam terhadap masukan bump dengan tinggi 4 cm dan lebar 1,2 m

Sistem Pasif State Feedback FER State Feedback Output Feedback FER Output Feedback

Percepatan (m/s2)


Defleksi ban (m)

Perbaikan (%)

3,2527 0,8348 1,5648 2,3520 2,3520

0,01334 0,01265 0,01305 0,01250 0,01250

0,00698 0,00234 0,00379 0,00512 0,00512

5,17 2,19 6,35 6,35

[74,34 [51,89 [27,69 [27,69

66,39] 45,70] 26,59] 26,59]

92

BAB VII KESIMPULAN DAN SARAN 7.1 KESIMPULAN Dari hasil pemodelan peredam FER dan simulasi sistem suspensi semi-aktif yang menggunakan peredam FER dengan menggunakan metoda kontrol linearquadratic-regulator (LQR) optimal , dapat ditarik beberapa kesimpulan sebagai berikut: 1. Pemodelan peredam kejut yang menggunakan fluida ER dengan pendekatan geometri

pelat

paralel

1-dimensi

dapat

dianggap

memadai

untuk

menggambarkan sistem peredam kejut berfluida ER sebagai peredam kejut continuously variable. 2. Model matematis peredam kejut memiliki gaya maksimal ketika diberikan medan listrik 4 kV/mm sebesar 546,37 N, sehingga bila pengontrol meminta masukan gaya yang melebihi kapasitasnya maka masukan gaya sistem kontrol akan terpotong sebesar gaya maksimum peredam yang digunakan tersebut. 3. Pengujian sistem suspensi pasif menggunakan masukan jalan bump dengan tinggi 0.04 m dan lebar bump 1.2 m ketika kecepatan kendaraan konstan 60 km/jam selama 1.5 detik menghasilkan percepatan vertikal (rms) badan kendaraan sebesar 3,9455 m/s2, defleksi suspensi (rms) 0.01197 m, dan defleksi ban sebesar 0.0147 m (rms). Pada kontrol optimal dengan umpan balik keadaan, sistem aktif dan FER, masing-masing, akan menghasilkan perbaikan pada percepatan (rms) badan kendaraan sebesar 58,99% (menjadi 1,6178 m/s2) dan 7,16% (3,6629 m/s2), dengan kompensasi menambah defleksi (rms) ruang suspensi sebesar 88,49% dan 13,94%, dan juga menambah defleksi (rms) ban sebesar 41,42% dan 4,98%. Pada kontrol optimal

yang

menggunakan

indeks

performansi

yang

melibatkan

pengoptimalan perintah gaya kontrol, sistem aktif dan FER, masing-masing, memberikan hasil perbaikan pada percepatan vertikal (rms) badan kendaraan sebesar 8,14% (3,625 m/s2) dan 7,39% (3,6538 m/s2), dengan kompensasi

93

bertambahnya defleksi ruang suspensi sebesar 20,33% dan 12,68% dan bertambahnya defleksi ban sebesar 8,15% dan 4,75%. 4. Pengujian sistem suspensi pasif menggunakan masukan jalan bump dengan tinggi 0.04 m dan lebar bump 1.2 m ketika kecepatan kendaraan konstan 20 km/jam selama 1.5 detik menghasilkan percepatan vertikal (rms) badan kendaraan sebesar 3,2527 m/s2, defleksi suspensi (rms) 0,01334 m, dan defleksi ban sebesar 0,00698 m (rms). Pada kontrol optimal dengan umpan balik keadaan, sistem aktif dan FER, masing-masing, akan menghasilkan perbaikan pada percepatan (rms) badan kendaraan sebesar 74,34% (menjadi 0,8348 m/s2) dan 51,89% (1,5648 m/s2), defleksi (rms) ruang suspensi sebesar 5,17% dan 2,19%, dan pada defleksi (rms) ban sebesar 66,39% dan 45,70%. Pada kontrol optimal yang menggunakan indeks performansi yang melibatkan pengoptimalan perintah gaya kontrol, sistem aktif dan FER memberikan hasil perbaikan pada percepatan vertikal (rms) badan kendaraan sebesar 27,69% (2,3520 m/s2), defleksi ruang suspensi sebesar 6,35% dan defleksi ban sebesar 26,59%.

7.2 SARAN UNTUK PENGEMBANGAN LEBIH LANJUT 1. Perlunya penerapan model peredam FER pada konstruksi fisik/prototype peredam, untuk menguji secara nyata fisik parameter-parameter dalam pemodelan yang sudah dilakukan. 2. Perancangan sistem suspensi FER perlu dilakukan pada orde sistem yang lebih kompleks dan dilakukan dengan berbagai macam profil permukaan jalan. 3. Perlunya penelitian lebih lanjut terhadap penerapan model peredam FER ini pada berbagai perancangan sistem kontrol selain sebagai sistem suspensi kendaraan.

DAFTAR PUSTAKA

1. Artikel, SHOCK ABSORBER USES ELECTRORHEOLOGICAL FLUID, Automotive Engineering Magazine vol. 100(6), Juni 1992 2. Fox, R.W., McDonald, A.T., INTRODUCTION TO FLUID MECHANICS, 4th ed., John-Willey & Sons, 1994 3. Kamath, G.M., Hurt, M.K., and Wereley, N.M., ANALYSIS AND TESTING OF

BINGHAM

PLASTIC

BEHAVIOR

IN

SEMI-ACTIVE

ELECTRORHEOLOGICAL FLUID DAMPERS, Smart Materials and Structures vol.5, pp. 576-590, 1996 4. Kamath, G.M., and Wereley, N.M., A NONLINEAR VISCOELASTICPLASTIC MODEL FOR ELECTRORHEOLOGICAL FLUIDS. Maryland: Smart Mater. Struct. 6, pp 351-359, 1997 5. ---, INTRODUCTION TO ER FLUID: NEW ERF, TX-ER SERIES, Nippon Shokubai Co.Ltd., 1999 6. Kaster, Judith, ELECTRORHEOLOGICAL MATERIALS PAPER. 1994 7. Hac, A., Youn, I., OPTIMAL SEMI-ACTIVE SUSPENSION WITH PREVIEW BASED ON A QUARTER CAR MODEL, Journal of Vibration and Acoustics, vol.114, 1992 8. Spencer

Jr,

B.F.,

Dyke,

S.J.,

Sain,

M.K.,

Carlson.

J.D.,

PHENOMENOLOGICAL MODEL OF A MAGNETORHEOLOGICAL DAMPER, to appear in the ASCE Journal of Engineering Mechanics, 10 Maret 1996 9. Kelly, S.G., FUNDAMENTALS OF MECHANICAL VIBRATIONS, 2ND International Editions, McGraw-Hill, 2000 10. McCabe, W.L., Smith, J.C, Harriott, P., UNIT OPERATIONS OF CHEMICAL ENGINEERING. 5th International Ed., McGraw-Hill, 1993

11. ---, A DESIGNER'S GUIDE TO HYDRAULIC SHOCK ABSORBER SELECTION. New York: Taylor Devices, Inc., 1999 12. Ogata, Katsuhiko, MODERN CONTROL ENGINEERING, 2nd Ed. New Delhi: Prentice-Hall of India, 1996 13. Radcliffs,

C.J.,

STATE

FEEDBACK

CONTROL

OF

ELECTRORHEOLOGICAL FLUIDS, ASME Int. Congress and Exhibition, Nov 17-22, 1996 (with 1/6/97 corrections) 14. Sudrajat, A., Nazaruddin, Y.Y., Chandra, A.D., Amir, R., MENENTUKAN INDEKS KENYAMANAN (RIDE COMFORT INDEX) KENDARAAN MENGGUNAKAN SENSOR AKSELEROMETER PADA TIGA ARAH PENGUKURAN, PPI–KIM, 1997 15. Hanafi, D., Nazaruddin, Y.Y., IDENTIFIKASI DINAMIKA SUSPENSI MODEL SEPEREMPAT KENDARAAN DENGAN MENGGUNAKAN DATA HASIL PENGUKURAN WAKTU NYATA, PPI–KIM, 1997 16. Lin, J., Kanellakopoulos, I., NONLINEAR DESIGN OF ACTIVE SUSPENSIONS, IEEE Control Systems Magazine vol. 17(3), June 1997 17. Doyle, John C., Francis, Bruce A., Tannenbaum, Allen R., FEEDBACK CONTROL THEORY, Singapore: Macmillan Publishing Company, 1992 18. Feigel, H., Romano, N., NEW VALVE TECHNOLOGY FOR ACTIVE SUSPENSION, SAE-SP-96/1136 #960727 19. Soliman,

A.M.A.,

El-Tawwab,

A.M.A.,

Crolla,

D.A.,

ADAPTIVE

CONTROL STRATEGIES FOR A SWITCHABLE DAMPER SUSPENSION SYSTEM, SAE-SP-96/1136 #960939 20. Sharma, C.S, MECHANICAL VIBRATIONS ANALYSIS. Delhi: Khanna Publishers, 1983 21. Shahian, B., and M. Hassul, CONTROL SYSTEM DESIGN USING MATLAB. Englewood Cliffs, NJ: Prentice-Hall, Inc., 1993 22. Glisson, T.H., INTRODUCTION TO SYSTEM ANALYSIS. Singapore: McGraw-Hill, Inc., 1987

23. Ziemer, R.E., Tranter, W.H., and D.R. Fannin, SIGNALS AND SYSTEMS: CONTINUOUS AND DISCRETE, 3rd Ed. New York: Macmillan Publishing, Co., 1993 24. Ogata, Katsuhiko, SOLVING CONTROL ENGINEERING PROBLEMS WITH MATLAB. Engelewood Cliffs, NJ: Prentice-Hall, Inc., 1994 25. The University of California, A CONTROL ALGORITHM UTILIZING ELECTRO-RHEOLOGICAL

(ER)

FLUIDS

FOR

VIBRATION

SUPPRESSION, U.S. Patent No. 5,398,785, UC Case No.: B92-068-1. Berkeley, CA: http://socrates.berkeley.edu/ ~otl 26. Margolis, Donald L., THE RESPONSE OF ACTIVE AND SEMI-ACTIVE SUSPENSIONS TO REALISTIC FEEDBACK SIGNALS, Vehicle System Dynamics vol. 11, pp. 267-282, 1982 27. Martinus, Donny, PERANCANGAN SISTEM KONTROL SEMI AKTIF PADA

PEREDAM

GETAR

MODEL

KENDARAAN

SETENGAH

MENGGUNAKAN KONTROL OPTIMAL, Tugas Akhir, Jurusan Teknik Fisika ITB, 1996 28. Setiadi, Sonny, STUDI PERANCANGAN KONTROL DAN SIMULASI SISTEM

SUSPENSI

AKTIF

MODEL

KENDARAAN

SETENGAH

MENGGUNAKAN METODE REGULATOR OPTIMAL, Tugas Akhir, Jurusan Teknik Fisika ITB, 1996 29. Rudy, PEMILIHAN PARAMETER FUNGSI PEMBOBOT DALAM PERANCANGAN KONTROL TEGAR UNTUK SISTEM SUSPENSI KENDARAAN DENGAN ALGORITMA GENETIK, Tugas Akhir, Jurusan Teknik Fisika ITB, 1998 30. Mudrig, A. Raymon, PERANCANGAN DAN PEMBUATAN PROTOTIPE SISTEM SUSPENSI AKTIF UNTUK KENDARAAN, Tugas Akhir, Jurusan Teknik Fisika ITB, 1999

A1

LAMPIRAN A LISTING PROGRAMS Berikut ini adalah listing program simulasi yang dilakukan petugas akhir. A.1 M-FILE: txer6.m % The Properties of TX-ER6 (R) Nippon Shokubai Co., Ltd., contains: % Data Table of Induced Shear Stress under DC and AC 50Hz Electric Field % at 200/sec Shear Rate (assumed constant for all condition) % Description of Data Tables: % EF

= electric field (kV/mm)

% tyac = induced shear stress (Pa) under AC 50Hz Electric Field % tydc = induced shear stress (Pa) under DC Electric Field % and hereby, produces polynomial equations for both induced shear stress % as functions of AC and DC electric field. % tyaca = induced shear stress (Pa) as a polynomial function of AC 50Hz % electric field % tydca = induced shear stress (Pa) as a polynomial function of DC electric % field % compiled by: A.D. Setyobudi, 133 94 080 - TF-ITB

(c)1999

EF = [0.0000; 0.0392; 0.0784; 0.1176; 0.1569; 0.1961; 0.2353; 0.2745; 0.3137; 0.3529; 0.3922; 0.4314; 0.4706; 0.5098; 0.5490; 0.5882; 0.6275; 0.6667; 0.7059; 0.7451; 0.7843; 0.8235; 0.8627; 0.9020; 0.9412; 0.9804; 1.0196; 1.0588; 1.0980; 1.1373; 1.1765; 1.2157; 1.2549; 1.2941; 1.3333; 1.3725; 1.4118; 1.4510; 1.4902; 1.5294; 1.5686; 1.6078; 1.6471; 1.6863; 1.7255; 1.7647; 1.8039; 1.8431; 1.8824; 1.9216; 1.9608; 2.0000; 2.0392; 2.0784; 2.1176; 2.1569; 2.1961; 2.2353; 2.2745; 2.3137; 2.3529; 2.3922; 2.4314; 2.4706; 2.5098; 2.5490; 2.5882; 2.6275; 2.6667; 2.7059; 2.7451; 2.7843; 2.8235; 2.8627; 2.9020; 2.9412; 2.9804; 3.0196; 3.0588; 3.0980; 3.1373; 3.1765; 3.2157; 3.2549; 3.2941; 3.3333; 3.3725; 3.4118; 3.4510; 3.4902; 3.5294; 3.5686; 3.6078; 3.6471; 3.6863; 3.7255; 3.7647; 3.8039; 3.8431; 3.8824; 3.9216; 3.9608; 4.0000];

ty= [0.7

0.7; 1.296456353

1.944684529; 2.592912705 3.88936905;

3.889369058 5.834053587;

5.185825411 7.778738116;

6.482281763 9.723422645;

7.778738116 11.66810717;

9.075194468 13.6127917;

10.37165082 15.55747623;

LAMPIRAN A

LISTING PROGRAMS

A2

11.66810717 17.50216076;

12.96456353 19.44684529;

14.26101988 21.39152982;

15.55747623 23.33621435;

16.85393258 25.28089888;

28.65168539 27.22558341;

29.49438202 29.17026793;

30.33707865 31.11495246;

32.02247191 33.05963699;

33.70786517 35.00432152;

35.39325843 36.94900605;

42.13483146 38.89369058;

50.56179775 40.83837511;

55.61797753 42.78305964;

58.98876404 44.72774417;

63.20224719 46.67242869;

75.84269663 48.61711322;

84.26966292 50.56179775;

84.26966292 54.77528090;

88.48314607 58.98876404;

92.69662921 63.20224719;

101.1235955 67.41573034;

109.5505618 75.84269663;

117.9775281 84.26966292;

117.9775281 88.48314607;

126.4044944 92.69662921;

130.6179775 101.1235955;

134.8314607 109.5505618;

143.2584270 113.7640449;

151.6853933 117.9775281;

160.1123596 126.4044944;

168.5393258 130.6179775;

176.9662921 134.8314607;

193.8202247 143.2584270;

202.2471910 151.6853933;

210.6741573 160.1123596;

219.1011236 164.3258427;

235.9550562 172.7528090;

244.3820225 185.3932584;

252.8089888 193.8202247;

269.6629213 202.2471910;

282.3033708 210.6741573;

294.9438202 219.1011236;

303.3707865 227.5280899;

320.2247191 235.9550562;

337.0786517 248.5955056;

353.9325843 261.2359551;

370.7865169 269.6629213;

379.2134831 278.0898876;

387.6404494 286.5168539;

404.4943820 294.9438202;

412.9213483 307.5842697;

429.7752809 320.2247191;

446.6292135 328.6516854;

455.0561798 337.0786517;

471.9101124 345.5056180;

488.7640449 362.3595506;

505.6179775 370.7865169;

522.4719101 387.6404494;

539.3258427 404.4943820;

556.1797753 412.9213483;

573.0337079 421.3483146;

589.8876404 429.7752809;

606.741573

623.5955056 455.0561798;

648.8764045 471.9101124;

674.1573034 488.7640449;

691.0112360 497.1910112;

724.7191011 514.0449438;

741.5730337 526.6853933;

783.7078652 543.5393258;

800.5617978 564.6067416;

842.6966292 573.0337079;

872.1910112 589.8876404;

893.2584270 606.7415730;

926.9662921 623.5955056;

446.6292135;

LAMPIRAN A

LISTING PROGRAMS

A3

952.2471910 640.4494382;

977.5280899 657.3033708;

1019.662921 665.7303371;

1053.370787 682.5842697;

1087.078652 707.8651685;

1112.359551 724.7191011;

1146.067416 733.1460674;

1179.775281 758.4269663;

1213.483146 775.2808989;

1247.191011 792.1348315;

1272.471910 808.9887640;

1306.179775 825.8426966;

1339.887640 851.1235955;

1365.168539 867.9775281;

1398.876404 884.8314607;

1424.157303 901.6853933;

1491.573034 918.5393258;

1500

943.8202247];

tydc = ty(:,1); % induced shear stress under dc electric field tyac = ty(:,2); % induced shear stress under ac electric field

%% Add or Remove '%' sign to hide or display plot! %disp('Press any key to plot Induced Shear Stress vs. Electric Field... ') %pause % Press a key to continue %figure(1), %plot(EF, tydc, '-'), hold on %title('Induced Shear Stress for TX-ER6 under Electric Field at 200/sec') %xlabel('Electric Field (kV/mm)'), ylabel('Induced Shear Stress(Pa)'),grid on %plot(EF, tyac, '.'), hold on %xlabel('Electric Field (kV/mm)'), ylabel('Induced Shear Stress(Pa)'),grid on %text(0.25, 1400, '(-) DC field'), text(0.25,1300, '(.) AC 50Hz field') %disp('Press a key to continue...') %pause % Press a key to continue

%% Returns the coefficients of least-squares-fit polynomial, in descending %% powers: tydc = aa(1).*EF.^5 + aa(2).*EF.^4 + ... + aa(5).*EF + aa(6) => %% 5 order polynomial [aa, sa] = polyfit(EF, tydc, 5); % for TX-ER6 under DC Electric field [bb, sb] = polyfit(EF, tyac, 5); % for TX-ER6 under AC 50Hz Electric field

%% Produces new ty from polynomial [tydca, delta] = polyval(aa, EF, sa); [tyaca, deltb] = polyval(bb, EF, sb);

%% Returns the coefficients of least-squares-fit polynomial, in descending %% powers: tydc = aa(1).*EF.^4 + ... + aa(4).*EF + aa(5) =>

LAMPIRAN A

LISTING PROGRAMS

A4

%% 4 order polynomial %[aa1, sa1] = polyfit(EF, tydc, 4); %%for TX-ER6 under DC Electric field %[bb1, sb1] = polyfit(EF, tyac, 4); %%for TX-ER6 under AC 50Hz Electric field %aa1=abs(aa1); %%Makes polynomial coefficient positive %% aa1= [1.8449

4.0684

43.2146

27.0725

-2.5308]; before abs() command

%% Produces new ty from polynomial %[tydca1, delta1] = polyval(aa1, EF, sa1); %[tyaca1, deltb1] = polyval(bb1, EF, sb1);

%disp('Next polynomial equations is determined by least-squares-fit of the') %disp('dynamic yield stress data as a function of electric field') %disp('Press any key to continue...') %pause %plot(EF, tydca, 'x'), grid on %xlabel('Electric Field (kV/mm)'), ylabel('Induced Shear Stress (Pa)') %plot(EF, tyaca, 'o'), grid on %xlabel('Electric Field (kV/mm)'), ylabel('Induced Shear Stress (Pa)') %text(0.25, 1200, '(x) Poly (DC)'), text(0.25,1100, '(o) Poly (AC)')

%clo=input('Close the ER figure window? (y/n) ','s'); %if (clo=='y') close; end %clear clo ty %% end of file

A.2 M-FILE: erdamp.m % ER Damper model, 1D Cylindrical Axisymmetric Approximation % used for: Force = cnr * (speed)

+ ct * ty *sgn(speed)

% ER Damper geometric data R = 0.02125; % = 21.25 mm; piston head radius (m) R1= 0.0233; % outer wall radius of inner electrode (m) a = 0.0005; % = 0.5 mm; gap separation length (m) R2= R1+a; % inner wall radius of outer electrode (m) L = 0.0635; % = 63.5 mm; duct length (m) mu= 0.035; % zero-voltage viscosity (Pa.s) AA=pi.*R^2; % piston head area (m^2) A2=pi.*R2^2; % outer electrode area (m^2)

LAMPIRAN A

LISTING PROGRAMS

A5

% Geometry coefficient cylindrical approximation: tonr=8*pi*L*(AA/A2)^2 *inv(1-(R1/R2)^4-(((1-(R1/R2)^2)^2)/log(R2/R1))); % Passive damping coefficient: cnr=mu*tonr; % Geometrical component of active damping coefficient: ct = 2*L*AA/a; % end of file

A.3 M-FILE: ermodel.m % Gives Velocity profiles through the gap VS. Gap Radius % Before execute this file, Users must give inputs as described below: % FF : trial force applied to the damper % ty : yield stress as a function of electric field (Pa) erdamp; clear tonr R; % çload er damper geometric data % load txer6 data: txer6; clear EF ty bb sb tyac tydc tydca tyaca delta deltb; FF=[100 400 500 600 1000 1500 2000 2500]; % ç trial forces (N) EF=[0 1 2 3 3.5 4]; % ç electric field input (kV/mm) for i=3 % change i to the required number of FF column number for j=length(EF):-1:1 %change the first number to length of EF matrix clear ri rii riii ubr1 ubr2 ubr3; %clears old values dp=-(FF(i))/AA; %delta pressure is always negative in our case EEF=EF(j); [ty, delta] = polyval(aa, EEF, sa); % gives yield stress point dpa=dp/(4*mu*L); dpp=2.*L*ty/dp; %plug thickness, dpp, is a function of applied force % Equation: ubr1(rpi)-ubr3(rpo)=0 s2=symsub('rpi',dpp); %oke s3=s2^2-R2^2-2*s2^2*log(s2/R2); %oke s4=s2-R2-s2*log(s2/R2); %oke s5='rpi*log(rpi)'; %oke s6=symmul('rpi',log(R1)); %oke s7=symmul('2*rpi',s5-s6); %oke s8=symsub('rpi',s5-s6)-R1; %oke s9=symsub('rpi^2',s7)-R1^2; %oke

LAMPIRAN A

LISTING PROGRAMS

A6

% Solve rpi and rpo, step by step rpii=solve(dpa*(s9-s3)-(ty/mu)*(s8+s4)); rpi=numeric(rpii);% converts symbolic form, rpii, to numeric form, rpi. rpoo=rpii-dpp; rpo=numeric(rpoo);% we've got both rpi and rpo clear s2 s3 s4 s5 s6 s7 s8 s9 rpii rpoo % clears unnecessary constants % Prevents rpi and rpo exceed R1 and R2 limits: if (rpiR2) rpi=R1; rpo=R2; end % Velocity distribution & volume flux for region 1 ri=[R1:0.0000001:rpi]; %qbr1=0.5*pi*dpa*((3*rpi^2-R1^2)*(rpi^2-R1^2)-4*rpi^4*log(rpi/R1)) ... %

-(pi*ty/(6*mu))*((rpi-R1)*(7*rpi^2+rpi*R1-2*R1^2)-...

%

6*rpi^3*log(rpi/R1));

for hi=1:length(ri) r=ri(hi); ubr1(:,hi)=dpa*(r^2-R1^2-2*rpi^2*log(r/R1))-(ty/mu)*(r-R1-... rpi*log(r/R1)); end % Velocity distribution & volume flux for region 2 (PLUG REGION) rii=[rpi:0.0000001:rpo]; r=rpi; %qbr2=dpa*(r^2-R1^2-2*rpi^2*log(r/R1))-(ty/mu)*(r-R1-rpi*log(r/R1)) ... %

*pi*(rpo^2-rpi^2);

for hii=1:length(rii) r=rpo; %ubr2(:,hii)=dpa*(r^2-R1^2-2*rpi^2*log(r/R1))-(ty/mu)*(r-R1-... %

rpi*log(r/R1));

ubr2(:,hii)=dpa*(r^2-R2^2-2*rpo^2*log(r/R2))+(ty/mu)*(r-R2-... rpo*log(r/R2)); end % Velocity distribution & volume flux for region 3 riii=[rpo:0.0000001:R2]; %qbr1=0.5*pi*dpa*((-3*rpo^2+R2^2)*(rpo^2-R2^2)+4*rpo^4*log(rpo/R2)) ... % %

+(pi*ty/(6*mu))*(-(rpo-R2)*(7*rpo^2+rpo*R2-... 2*R2^2)+6*rpo^3*log(rpo/R2));

for hiii=1:length(riii)

LAMPIRAN A

LISTING PROGRAMS

A7

r=riii(hiii); ubr3(:,hiii)=dpa*(r^2-R2^2-2*rpo^2*log(r/R2))+(ty/mu)*(r-R2-... rpo*log(r/R2)); end

rr=[ri rii riii]; ubr=[ubr1 ubr2 ubr3]; %% Add/remove '%' in front of below lines to toggle the figure %figure(1); %plot(ubr, rr); hold on %xlabel('Velocity (m/s)'); ylabel('Radius (m)'); end end % end of file

A.4 M-FILE: ermodel1.m % This m file gives plug thickness versus electric field diagram erdamp; clear tonr R;% ç load er damper geometric data txer6; clear ty bb sb tyac tydc deltb; % ç load txer6 data FF=[100 500 1000 1500 2000]; % ç trial forces (N) for i=1:5 dp=abs(FF(i))/AA; EEF=0:0.05:4; % ç electric field increase steps [ty, delta] = polyval(aa, EEF, sa); dop=2000*L*ty/dp; %plug thickness,dop, as a function of applied force [hh,h]=size(dop); for hi=1:h if (dop(hi)>1000*a) dop(hi)=1000*a; end end plot(EEF,dop); hold on; xlabel('Electric field (kV/mm)'); ylabel('Plug thickness (mm)') end %end of file

LAMPIRAN A

LISTING PROGRAMS

A8

A.5 M-FILE: ermodel2.m % This m file gives plug thickness versus force diagram erdamp; clear tonr R; % ç load er damper geometric data txer6; clear EF ty bb sb tyac tydc deltb; % ç load txer6 data EE=[0.5 1 2 2.5 3 4]; % for constant electric field for i=1:6 EF=EE(i); [ty, delta] = polyval(aa, EF, sa); ff=[50:2:2000];% ç varying applied force dp=abs(ff)/AA; dop=2000*L.*ty./dp; %plug thickness,dop, as a function of applied force plot(ff,dop); hold on; axis([0 2000 0 5]); xlabel('force (N)'); ylabel('plug thickness (mm)'); end %end of file

A.6 M-FILE:ermodel3.m % This m file gives plug thickness versus applied force diagram % under some electric field constants % ------ inter active ----------------------------------------erdamp; clear tonr R;% ç load er damper geometric data txer6; clear EF ty bb sb tyac tydc deltb; % ç load txer6 data %EE=[0.5 1 2 2.5 3 4]; % for constant EE electric field flag=0; %sign to begin looping while(flag==0) ef=input('Type electric field magnitude (0-4 kV/mm) : '); ef=abs(ef); if isempty(ef) ef=0; else if ef>4 ef=4; disp('Maximum allowed electric field is 4 kV/mm, your ef now: 4 kV/mm'); end end

LAMPIRAN A

LISTING PROGRAMS

A9

[ty,delta]=polyval(aa,ef,sa); ff=[1:1:1000];% ç varying applied force dp=abs(ff)/AA; % plug thickness,dop, as a function of applied force dop=2000.*L.*ty./dp; [hh,h]=size(dop); for hi=1:h if (dop(hi)>1000*a) dop(hi)=1000*a; end end plot(ff,dop,'-'); hold on; xlabel('force (N)'); ylabel('plug thickness (mm)'); sat=input('Satisfy? (y/n)','s');if (sat=='y') flag=1; end; %End of looping end %end of file

A.7 M-FILE: cardata.m %Contains data for FER quarter car model erdamp; %ç load er-damper parameters % Minimum Damping constant (cmin) cmin = cnr; % = 2901.9; dalam N.s/m clear R L a mu R1 R2 AA A2 le cn ct tonr cnr %configuration no. H : c2m1=10; m22m1=0.125; k12m1=60; k22m1=500; m1= cmin/c2m1; m2= m1.*m22m1; k1= m1.*k12m1; k2= m1.*k22m1;

% State Space Equation Matrices % vdot = A.*v + [B1|B2].*[u; w] % y = C2.*v + [D1|D2].*[u; w] A=[0 1

0

-1; -k1/m1

-cmin/m1

0

0 0

0

1;

cmin/m2

-k2/m2 -cmin/m2];

k1/m2

B1=[0; -1/m1; 0;

1/m2];

B2=[0; 0;

0];

-1;

C2=[-k1/m1 -cmin/m1 1

0

0

cmin/m1;

cmin/m1; 0 0;

LAMPIRAN A

LISTING PROGRAMS

A10

0

0

1 0];

D1=[1/m1; 0; 0]; D2=[0;

0; 0];

%end of file

A.8 M-FILE: oplores.m % Semi-active suspension system on a Quarter-Car model using ERF-based damper % This file contains Open-loop response of the system cardata; % ç load model data % Obtain transfer function for y(1), y(2), y(3) UNDER INPUT W [num12,den2]=ss2tf(A,B2,C(1,:),D2(1,1));% t/f of body acceleration (y(1)) [num22,den2]=ss2tf(A,B2,C(2,:),D2(2,1));% t/f of suspension deflection (y(2)) [num32,den2]=ss2tf(A,B2,C(3,:),D2(3,1));% t/f of tire deflection (y(3)) w=logspace(0,3,300); [mag2,phase2,w]=bode(num12,den2,w); [mag3,phase3,w]=bode(num22,den2,w); [mag4,phase4,w]=bode(num32,den2,w); % Mencari frequensi (pribadi) peak: 1: under input w mak2=max(mag2);[mm2,nn2]=size(mag2); indek2=(mak2*ones(mm2,nn2)==mag2); indekw2=indek2.*w; wmax2=sum(indekw2); % Mencari frequensi (pribadi) peak: 2: under input w mak3=max(mag3);[mm3,nn3]=size(mag3); indek3=(mak3*ones(mm3,nn3)==mag3); indekw3=indek3.*w; wmax3=sum(indekw3); % Mencari frequensi (pribadi) peak: 3: under input w mak4=max(mag4);[mm4,nn4]=size(mag4); indek4=(mak4*ones(mm4,nn4)==mag4); indekw4=indek4.*w; wmax4=sum(indekw4); figure(1); subplot(221);ma2 = 20*log10(mag2); semilogx(w,ma2); grid on; text(wmax2-2,max(ma2)-3,['max='num2str(wmax2)'rad/sec,'num2str(max(ma2))... 'dB']); title('Body Acceleration');xlabel('Frequency (rad/sec)'); ylabel('Magnitude (dB)');

LAMPIRAN A

LISTING PROGRAMS

A11

subplot(222);ma3=20*log10(mag3); semilogx(w,ma3);grid on; %text(wmax3-2,max(ma3)-3,['max='num2str(wmax3)'rad/sec,'num2str(max(ma3))... 'dB']); title('Susp. Deflection');xlabel('Frequency (rad/sec)'); ylabel('Magnitude (dB)'); subplot(224);ma4=20*log10(mag4); semilogx(w,ma4); grid on; %text(wmax4-2,max(ma4)-3,['max='num2str(wmax4)'rad/sec,'num2str(max(ma4))... 'dB']); title('Tyre Deflection');xlabel('Frequency (rad/sec)'); ylabel('Magnitude (dB)'); %end of file

A.9 M-FILE: ferlaw.m function [fn]=ferlaw(u0,x2,x4) % [fn, ef]=ferlaw(u0, x2, x4) % Adaptive Control law for Semi-active Continuous Damper % % Gives a new force value fn and

electric field ef

% for a required control input force u0 and velocity (x2-x4) % x2: sprung mass velocity (m/s) % x4: unsprung mass velocity (m/s) clc ni = nargin; error(nargchk(1,3,ni)); clear ni txer6; %load txer6 polynomial curve data erdamp; %load erdamper data

umin=ct*polyval(aa,0); umax=ct*polyval(aa,4); %min force = 3.97313002642733; %max force = 5.463709527096947e+002;

% The following conditional structures determine new f and ef

values

vel=sign(x2-x4); switch vel

case 1 %umin
LAMPIRAN A

LISTING PROGRAMS

A12

if (abs(u0)>umin)&(abs(u0)0)&(abs(u0)<=umin) fn=sign(u0)*umin; ef=0; elseif (abs(u0)>=umax) fn=sign(u0)*umax; ef=4; end

case -1 %uminumin)&(abs(u0)0)&(abs(u0)<=umin) fn=sign(u0)*umin; ef=0; elseif (abs(u0)>=umax) fn=sign(u0)*umax; ef=-4; end

case 0 fn=0;

end % end of file

A.10 M-FILE: carmodel.m %% Builds an FER based quarter car MIMO model format short e erdamp; %load er-damper parameters %% Remove or add '%' to toggle the functions/syntaxes/graphics ! cardata; % load cardata

LAMPIRAN A

LISTING PROGRAMS

A13

%% specify the state space as an LTI object, and attach names to the %% states, inputs, and outputs: %states={'sudef' 'spvel' 'usdef' 'usvel'}; %inputs={'erforce' 'roadvel'}; %outputs={'spacc' 'sudef' 'usdef'}; %sys_er=ss(A,B,C,D,'statename',states,'inputname',inputs,'outputname',outputs ); %% you can display this continuous-time model by typing: sys_er %% our model has two inputs (er force, and road velocity), and three ouputs %% (sprung acceleration, suspension deflection, and unsprung deflection).

%% ---- Test Controllable & Observable ---------[mA,nA]=size(A); if (rank(ctrb(A,B1)) ~= mA) | (rank(obsv(A,C)) ~= nA) disp(' ') disp('

********************************')

disp('

Uncontrollable or Unobservable')

disp('

Check the plant matrix !!')

disp('

********************************')

break end

%% ----- Optimal Gain: State weighting ----%% State Weighting Inputs ro1=1000; ro2=10000; %% <-- you can change this values. They're not unique ! ro3=1; ro4=ro3/m1^2;%% <-- you can change this values. They're not unique !

%% State Weighting matrices Q=[(ro3*(k1.^2)/(m1.^2))+ro1 ro3*k1*cmin/(m1.^2) 0 -ro3*k1*cmin/m1.^2

ro3*k1*cmin/(m1.^2)

0 -ro3*k1*cmin/m1.^2;

(ro3*cmin.^2)/m1.^2

0 -(ro3*cmin.^2)/m1.^2;

0 -(ro3*cmin.^2)/m1.^2

ro2

0;

0 (ro3*cmin.^2)/m1.^2];

N=[ro3*k1/m1.^2; ro3*cmin/m1.^2; 0; -ro3*cmin/m1.^2]; R=ro3/m1.^2;%

%% Regulator, New States: AN = A - B1*inv(R)*N'; CN = C - D1*inv(R)*N';

LAMPIRAN A

LISTING PROGRAMS

A14

QN = Q - N*inv(R)*N';

QQ = QN.^0.5; if (rank(ctrb(AN,B)) ~= mA) | (rank(obsv(AN,QQ)) ~= nA) disp(' ') disp('

***************************************')

disp('

Uncontrollable or Unobservable')

disp(' disp('

Check the plant weighting matrix !!') ***************************************')

break end [K,Pa,clp]=lqr(A,B1,Q,R,N); %% Thus far, we've got open-loop system matrices: %% A, B1, C and D1 from carmodel.m, and weighting constants Q, N, R. %% The open-loop matrices: aol=A; bol=B2; col=C; dol=D2; %% The closed-loop regulator matrices are introduced below: acl=A-B1*K; bcl=B2; ccl=C-D1*K; dcl=D2;

%% --------- Optimal Gain: Output weighting --------Qy=diag([ro3 ro1 ro2]); Ry=[ro4];

%% New PI Equation QQy=C'*Qy*C; RRy=Ry+D1'*Qy*D1; Ny=C'*Qy*D1;

%% New State & New Weighting Matrices %An=A-B*inv(Ry)*N'; %Bn=B; %Cn=C-D*inv(Ry)*N';

LAMPIRAN A

LISTING PROGRAMS

A15

%Dn=D; %Qn=Qy-N*inv(Ry)*N'; %Rn=Ry; [Ky, Pay, clpy]=lqr(A,B1,QQy,RRy,Ny);

%% The new closed-loop matrices for output weighting are introduced below: acly=A-B1*Ky; bcly=B2; ccly=C-D1*Ky; dcly=D2; % end of file

A.11 M-FILE: carplot.m carmodel; %% Then we will see the transfer function of each inputs to each output: %tf(sys_er); %step(sys_er(:,2)); %h11=tf(num11,den); h12=tf(num12,den); %h21=tf(num21,den); h22=tf(num22,den); %h31=tf(num31,den); h32=tf(num32,den); %H=[h11 h12; h21 h22; h31 h32]; %% H is system's transfer functions concatenation

%%

------------------------------------

%% Specify sampling time & points t=(0:0.005:2.5)'; w=logspace(0,3,600); oup=step(aol,bol,col,dol,1,t); % Compute passive step response oup=oup/oup(length(oup)); out=step(acl,bcl,ccl,dcl,1,t); % Compute active step response out=out/out(length(out)); outy=step(acly,bcly,ccly,dcly,1,t); %Compute output weighting step response outy=outy/outy(length(outy)); [mag1,phase,w]=bode(aol,bol,col,dol,1,w);% mag & phase of open-loop [cmag1,cphase,w]=bode(acl,bcl,ccl,dcl,1,w);% mag & phase of closed-loop

LAMPIRAN A

LISTING PROGRAMS

A16

[cmagy,cphasey,w]=bode(acly,bcly,ccly,dcly,1,w);%mag&phase of output weight %clpq=eig(acl); % closed-loop poles array yp1=oup(:,1);% step response of passive body acc. %yp2=oup(:,2);% step response of passive susp. deflection %yp3=oup(:,3);% step response of passive tyre deflection yq1=out(:,1);% step response of body acc. %yq2=out(:,2);% step response of susp. deflection %yq3=out(:,3);% step response of tyre deflection yqy1=outy(:,1); % step response of output weighting moq1=20*log10(mag1(:,1)); % open-loop Bode magnitude of Body acc. moq2=20*log10(mag1(:,2)); % open-loop Bode magnitude of susp. defl. moq3=20*log10(mag1(:,3)); % open-loop Bode magnitude of tyre deflection %poq1=phase(:,1); % open-loop Bode phase of body acc. %poq2=phase(:,2); % open-loop Bode phase of susp. defl. %poq3=phase(:,3); % open-loop Bode phase of tyre defl. mcq1=20*log10(cmag1(:,1)); % closed-loop magnitude of body acc. mcq2=20*log10(cmag1(:,2)); % c.l. magnitude of susp. defl. mcq3=20*log10(cmag1(:,3)); % c.l. magnitude of tyre defl. %pcq1=cphase(:,1); % phase of body acc. %pcq2=cphase(:,2); % phase of susp. defl. %pcq3=cphase(:,3); % phase of tyre defl. [posp1,trp1,ts2p1,tpp1]=stepchar(t,yp1); [pos1, tr1, ts21, tp1]=stepchar(t,yq1); [posy1,try1,ts2y1,tpy1]=stepchar(t,yqy1); %[pos2(i,:), tr2(i,:), ts22(i,:), tp2(i,:)]=stepchar(t,yq2(:,i)); %[pos3(i,:), tr3(i,:), ts23(i,:), tp3(i,:)]=stepchar(t,yq3(:,i)); mcqy1=20*log10(cmagy(:,1)); % c/l magnitude of body acc. output weighting mcqy2=20*log10(cmagy(:,2)); % c/l magnitude of susp. defl. output weighting mcqy3=20*log10(cmagy(:,3)); % c/l magnitude of tyre defl. output weighting %pcqy1=cphasey(:,1); % phase of body acc. output weighting %pcqy2=cphasey(:,2); % phase of susp. defl. output weighting %pcqy3=cphasey(:,3); % phase of tyre defl. output weighting

%% ------------------------- WARNING !!! ----------------------------%% If your machine has a low memory resource, %% Copy these lines to your clipboard first one line to another, %% if your computer is having low memory resources. %% Otherwise your figures will not displayed properly.

LAMPIRAN A

LISTING PROGRAMS

A17

%% ------------------------------------------------------------------close all; %figure(1); title('Step response'); %subplot(221); %plot(t,yq1,'r--'); hold on; title('Step response of body acceleration'); grid on %plot(t,yp1,'r:'); hold on; %plot(t,yqy1,'r-.'); %subplot(222); %plot(t,yq2);hold on; title('Step response of suspension deflection');grid on %subplot(223); %plot(t,yq3); hold on;title('Step response of tyre deflection'); grid on figure(2); %subplot(221); semilogx(w, moq1,'b-');hold on ; semilogx(w,mcq1,'r--');hold on;title('Body acceleration'); grid on semilogx(w,mcqy1,'r-.'); figure(3); %subplot(222); semilogx(w, moq2,'b-');hold on ; semilogx(w,mcq2,'r--'); hold on; title('Suspension deflection'); grid on semilogx(w,mcqy2,'r-.'); figure(4);%subplot(223); semilogx(w, moq3,'b-'); hold on ;%title('Tyre deflection'); grid on semilogx(w,mcq3,'r--'); hold on; title('Tyre deflection'); grid on semilogx(w,mcqy3,'r-.'); %figure(5); %plot(t,yq2,'r--');hold on; title('Step response of suspension deflection'); %grid on %plot(t,yp2,':'); %subplot(221); %semilogx(w,poq1); title('Body acceleration'); hold on ;grid on %subplot(222); %semilogx(w,poq2); title('Suspension deflection'); hold on ;grid on %subplot(223); %semilogx(w,poq3); hold on ;title('Tyre deflection'); grid on %figure(6); %subplot(221);

LAMPIRAN A

LISTING PROGRAMS

A18

%semilogx(w,mcq1, lin(i)); hold on ;title('Body acceleration'); grid on %subplot(222); %semilogx(w,mcq2); hold on ;title('Suspension deflection'); grid on %subplot(223); %semilogx(w,mcq3); hold on ;title('Tyre deflection'); grid on clear c2m1 m22m1 k12m1 k22m1 cmin m1 m2 k1 k2 QQ mAN nAN %% end of file

LAMPIRAN A

LISTING PROGRAMS

B1

LAMPIRAN B VARIASI PARAMETER PEREDAMAN Berikut ini adalah koefisien redaman hasil dari perubahan pada jarak celah antar elektroda dan perubahan panjang saluran celah. Cnr: koefisien redaman dengan pendekatan geometri 1D Cylindrical Axisymmetric. Cnp: koefisien redaman dengan pendekatan geometri 1D Pelat Parallel. Tb: koefisien geometri pada gaya peredaman aktif (fb). B.1 PERUBAHAN KOEFISIEN REDAMAN DENGAN JARAK CELAH KONSTAN (= 2 mm) % Perubahan koefisien redaman terhadap perubahan panjang. Duct length (mm) | Cnr

|

Cnp

| Tb

0.0535

2444.8628

2444.8812

0.30359

0.0555

2536.2596

2536.2786

0.31493

0.0575

2627.6563

2627.6761

0.32628

0.0595

2719.0531

2719.0735

0.33763

0.0615

2810.4498

2810.4709

0.34898

0.0635

2901.8465

2901.8684

0.36033

0.0655

2993.2433

2993.2658

0.37168

0.0675

3084.64

3084.6632

0.38303

0.0695

3176.0368

3176.0606

0.39438

0.0715

3267.4335

3267.4581

0.40573

0.0735

3358.8303

3358.8555

0.41708

0.0755

3450.227

3450.2529

0.42842

0.0775

3541.6237

3541.6503

0.43977

0.0795

3633.0205

3633.0478

0.45112

0.0815

3724.4172

3724.4452

0.46247

0.0835

3815.814

3815.8426

0.47382

B.2 PERUBAHAN KOEFISIEN REDAMAN DENGAN PANJANG SALURAN CELAH KONSTAN (= 63.5 mm) % Perubahan koefisien redaman terhadap perubahan celah elektroda Gap (m)

|

Cnr

|

Cnp

| Tb

0.0002

45632.2893

45341.693

0.90083

0.00025

23338.7816

23214.9468

0.72066

LAMPIRAN B


B2

0.0003

13491.8292

13434.5757

0.60055

0.00035

8487.2557

8460.2576

0.51476

0.0004

5679.7432

5667.7116

0.45041

0.00045

3984.8209

3980.615

0.40037

0.0005

2901.8465

2901.8684

0.36033

0.00055

2177.8849

2180.2166

0.32757

0.0006

1675.746

1679.322

0.30028

0.00065

1316.6225

1320.8322

0.27718

0.0007

1053.0452

1057.5322

0.25738

0.00075

855.2589

859.8128

0.24022

0.0008

703.9666

708.464

0.22521

B.3 LISTING PROGRAM VARIASI PARAMETER REDAMAN

% ER Damper geometric data R =0.02125; % = 21.25 mm; piston head radius (m) L =0.0635; % = 63.5 mm; duct length (m) a =0.0005; % = 0.5 mm; gap separation length (m) mu =0.035; % zero-voltage viscosity (Pa.s) R1= 0.0233; % outer wall radius of inner electrode (m) R2=R1+a; % inner wall radius of outer electrode (m) AA=pi.*R^2; % piston head area (m^2) A2=pi.*R2^2; % outer electrode area (m^2) le=pi.*0.0471; % duct depth

% ------ varying duct length ------LL=[0.0535:0.0020:0.0835]; disp(' % Perubahan koefisien redaman terhadap perubahan panjang'); disp(' Duct length (mm) | Cnr

|

Cn

| Tb');

[junk, ii]=size(LL);

for i=1:ii cn=(12*LL(:,i)*AA^2*mu)/(a^3*le);%passive damping coef parallel plate approx. ct=2*LL(:,i)*AA/a; %geometrical component of active damping coefficient % Geometry coefficient for cylindrical approximation tonr=8*pi*LL(:,i)*(AA/A2)^2 *inv(1-(R1/R2)^4-(((1-(R1/R2)^2)^2)/log(R2/R1))); % Passive damping coefficient cnr=mu*tonr; cnri(i,:)=cnr;

LAMPIRAN B


B3

disp([' 'num2str(LL(:,i))' 'num2str(cnr)' ' num2str(cn) '

'num2str(ct)]);

disp(''); end clear ii

% ------- varying gap distance ------aa=[0.0002:0.00005:0.0008]; disp(' % Perubahan koefisien redaman terhadap perubahan celah elektroda'); disp(' Gap (m) | Cnr

|

Cn

|

Tb');

[junk, ii]=size(aa); for i=1:ii R2=R1+aa(:,i); % inner wall radius of outer electrode (m) A2=pi.*R2^2; % outer electrode area (m^2) cn=(12*AA^2*mu*L)/(aa(:,i)^3*le);%passive damping coef parallel plate approx. ct = 2*L*AA/aa(:,i); % geometrical component of active damping coefficient % Geometry coefficient for cylindrical approximation tonr=8*pi*L*(AA/A2)^2 *inv(1-(R1/R2)^4-(((1-(R1/R2)^2)^2)/log(R2/R1))); % Passive damping coefficient cnr=mu*tonr; cnri(i,:)=cnr; disp([''num2str(aa(:,i))' 'num2str(cnr)' ' num2str(cn) '

' num2str(ct)]);

end clear ii % end of file

LAMPIRAN B


C1 t Clock

Waktu

AN OPTIMAL ACTIVE SUSPENSION BASED ON QUARTER-CAR MODEL WITH LINEAR QUADRATIC CONTROLLER

Sine Wave1 road

du/dt Product

Derivative

yy1

Plant + Regulator (Output Feedback)

Demux

yy2 yy3

Demux1

Step yp1 Plant Suspensi Pasif

Demux

yp2

run first: carmodel.m

yp3 Demux

y1 Plant + Regulator

Demux

y2 y3

Demux2

y1a K*u Demux

y2a

D K*u

1 s

B Sum

y3a

Sum2

K*u

Demux4 C

x

A

K*u K*u G

uoper

uop

ER Damper MATLAB Function

-K

STATE WEIGHTED CONTROLLER

K*u Mux

Demux

yy1a K*u Demux

yy2a

D1 K*u

1 s

B1 Sum1

A1

K*u

yy3a

Sum3 Demux6

x

C1

K*u K*u G1

uopery

ER Damper1 MATLAB Function

uopy

-K1

OUTPUT WEIGHTED CONTROLLER

K*u Mux

Demux

LAMPIRAN C - DIAGRAM BLOK SIMULASI

SIMULASI DAN PERANCANGAN SISTEM KONTROL SUSPENSI SEMI AKTIF MODEL SEPEREMPAT KENDARAAN

Recommend Documents