Journal of Mechanical Design and Testing 2(1), (2020), 47-65 DOI: http://dx.doi.org/10.22146/jmdt.v2i1.53724 e-ISSN: 2685-8029, ISSN: 2716-1293

## Simulasi Numerik Prediksi Erosi Aliran Fasa Cair-Gas Berpasir Pada Belokan Pipa Sumur Minyak Mentah Menggunakan Model Finnie, McLaury dan Oka
H. J. Atmaja * , S. Kamal, Sugiyono
Departemen Teknik Mesin dan Industri, Fakultas Teknik, Universitas Gadjah Mada Jalan Grafika No. 2, Yogyakarta 55281 Indonesia E-mail: *jaya_77a@yahoo.com

## Abstrak
Keberadaan pasir dalam fluida produksi di industri perminyakan merupakan hal yang tak dapat dihindari.  Pada  kecepatan  fluida  tertentu,  pasir  berpotensi  menimbulkan  erosi  pada  dinding perpipaan yang apabila proses tersebut terjadi secara terus menerus maka dapat menyebabkan penipisan dinding pipa, mengurangi tekanan kerja maksimum fluida yang diijinkan, kegagalan pipa dan akhirnya fluida hidrokarbon keluar dari sistem tertutup pipa. Kondisi ini berpotensi menimbulkan kerugian dalam hal keselamatan kerja, kehilangan finansial dan isu lingkungan. Kompleksitas mekanisme erosi akibat fluida berpasir didalam perpipaan memerlukan bantuan perangkat  lunak  untuk  menghitung  besaran  dan  profil  erosi  didalam  pipa.  Makalah  ini memprediksi laju erosi aliran multifasa cair-gas-padat pada belokan pipa yang didominasi fasa cair dengan CFD menggunakan ANSYS Fluent 18.2 dengan validasi berdasarkan eksperimen yang dilakukan peneliti lain yang sudah dipublikasi. Dari hasil simulasi, laju erosi dengan model Oka  lebih  mendekati  hasil  eksperimen  dibandingkan  model  Finnie  dan  McLaury.  Hal  ini dimungkinkan karena kondisi eksperimen lebih sesuai dengan model erosi yang dibangun oleh Oka.  Hasil  simulasi  lanjutan  menunjukkan  bahwa  laju  erosi  pada  pipa  meningkat  dengan kenaikan laju alir dari minyak mentah dan gas, ukuran partikel dan laju alir partikel. Laju erosi menurun dengan meningkatnya diameter pipa.
Kata kunci : Komputasi dinamika fluida, laju erosi maksimum, laju alir minyak mentah, laju alir gas, laju alir partikel, ukuran partikel, diameter pipa.

## Abstract
The sand presence in production fluid in oil and gas industry is unavoidable. At a certain fluid velocity, sand potentially causing erosion in the piping wall. If the process occurs continuously it  can  reduce  pipe  wall  thickness,  reducing  the maximum allowable  working pressure  of  the fluid, piping failure and hydrocarbon loss of containment. This condition will result in losses in terms  of  worker  safety,  financial  and  environmental  issues.  The  complexity  of  the  erosion mechanism due to sandy fluid in the piping system requires software to calculate the magnitude and pattern of erosion in the pipe. This study predicted the erosion rate of liquid dominated multiphase  flow  at  pipe  elbow  by  using  CFD  ANSYS  Fluent  18.2  and  validated  based  on published experimental studies conducted by other researcher. The simulation results show that
erosion rate of Oka model has the lowest relative difference with experimental result compared with Finnie and McLaury model. This may be caused by experimental condition that was more suitable with erosion model developed by Oka. Further simulation revealed that erosion rate at the  pipe  wall  increased  with  increasing  crude  oil  and  gas  flowrate,  particle  size  and  particle flowrate. Erosion rate decreased by increasing pipe diameter.
Keywords : Computational Fluid Dynamics (CFD), maximum erosion rate, crude oil flowrate, particle mass flowrate, particle size, pipe diameter.

## 1. PENDAHULUAN
Pengangkatan  ( lifting ) minyak mentah  dari reservoir ke permukaan, baik dengan  aliran alami ( natural flow ) atau aliran dengan bantuan pompa atau gas yang diinjeksikan ( artificial lift ) memungkinkan material padatan selain minyak mentah (misalnya pasir) ikut terbawa. Material pasir atau padatan lain yang berasal dari reservoir ketika terbawa ke permukaan akan menimbulkan masalah  pada  proses  pengolahan  fluida  dan  peralatan  penunjang  produksi.  Pasir  dapat menimbulkan erosi, kebuntuan pipa dan korosi dibawah endapan pasir.
Pencegahan  pasir  untuk  naik  ke  permukaan  sebenarnya  dapat  dilakukan  dengan menggunakan  peralatan  pekerjaan  sumur  ( workover) atau  perawatan sumur  ( well  service )  seperti pemasangan saringan pasir ( sand screen ) dan gravel pack didalam tubing annulus . Saringan pasir tidak mampu  untuk  menyaring  partikel  kurang  dari  50  µm  untuk  tidak  terbawa  aliran  fluida  ke permukaan. Begitu pula gravel pack lebih efektif untuk mencegah material padat yang kasar atau berukuran  tertentu  saja,  sedangkan  material  yang  lebih  kecil  atau  halus  masih  bisa  melewati peralatan tersebut (Parsi et al, 2014).
Selain  itu,  pencegahan  material  padatan  dari reservoir ke  permukaan  berpotensi menimbulkan  masalah  lain  yang  lebih  besar  yaitu  terkumpulnya  material  padat  tersebut  di reservoir atau tubing produksi yang justru akan menghambat aliran fluida itu sendiri. Untuk itu penanganan  pasir  pada  fluida  produksi  bergeser  dari  upaya  untuk  pencegahan  ( prevention ) menjadi pengelolaan ( management ). Artinya pasir dengan ukuran tertentu diperkenankan untuk naik  ke  permukaan  bersama  fluida  produksi,  tetapi  dikelola  dengan  cara-cara  tertentu  di permukaan.
Keberadaan  pasir  dalam  fluida  produksi  sudah  menjadi  hal  yang  umum  ditemukan pada proses produksi minyak mentah. Properti pasir pada kecepatan aliran fluida yang rendah sangat  rentan  terhadap  pengendapan  (sedimentasi)  yang  dapat menyebabkan hambatan  pada aliran  fluida,  meningkatkan  jatuh  tekanan  ( pressure  drop )  dan  mempercepat  proses  terjadinya korosi  pada  dinding  dalam  pipa  terutama  pada  bagian  yang  tertutupi  oleh  pasir  tersebut. Sementara pada kecepatan fluida yang tinggi, adanya pasir berpotensi menimbulkan erosi yang apabila proses tersebut terjadi secara terus menerus maka dapat menyebabkan penipisan dinding pipa, mengurangi tekanan kerja maksimum fluida yang diijinkan, kegagalan pipa dan akhirnya fluida hidrokarbon keluar dari sistem tertutup pipa. Kondisi ini akan berpotensi menimbulkan kerugian dalam hal keselamatan pekerja, kehilangan finansial dan isu lingkungan.
Studi erosi pada sistem perpipaan sudah intensif dilakukan para peneliti sejak lebih dari lima dekade tetapi kebanyakan penelitian berfokus pada aliran multifasa gas-cair-padat dengan fasa dominan adalah gas (Edwards et al, 2001; Humphrey, 1990; Salama, 2000), sementara masih sangat sedikit yang melakukan penelitian baik berbasis eksperimen maupun komputasi dinamika fluida (CFD) pada aliran slurry yang didominasi oleh fasa cair.
Dalam studi ini dianalisa pengaruh beberapa parameter aliran gas-cair-padat ( slurry ) yang didominasi  fasa  cair  terhadap  potensi  terjadinya  erosi  pada  sistem  perpipaan  sumur  migas hususnya  pada  belokan  pipa.  Simulasi  dilakukan  dengan  menggunakan software komputasi dinamika  fluida  ANSYS  Fluent  18.2  untuk  memprediksi  erosi  dan  memahami  secara  lebih mendalam parameter-parameter yang mempengaruhi laju erosi.

## 2. SIMULASI NUMERIK
Simulasi  pada  studi  ini menggunakan software  ANSYS Fluent  18.2.  Pada bagian ini akan dibahas secara lebih rinci geometri komputasi yang digunakan, mesh dan model multifasa yang dipilih.

## 2.1. Domain Komputasi
Simulasi numerik pada studi ini dilakukan untuk mensimulasikan laju erosi akibat pasir pada aliran multifasa cair-gas dengan geometri seperti ditunjukkan pada Gambar 1. Geometri adalah pipa vertikal dengan belokan keluaran arah horizontal. Diameter pipa adalah 50,8 mm dengan radius curvature 1,5 kali diameter. Panjang pipa sebelum dan setelah belokan masing masing adalah 42 dan 3 kali diameter.
Untuk  mendapatkan  mesh  terstruktur  maka  prosedur  diterapkan  dengan  membagi geometri menjadi beberapa bagian seperti ditunjukkan pada Gambar 1. Penampang pipa dibagi menjadi lima luasan, dibentuk dari persegi dengan pusat yang sama dengan lingkaran dan garis terdekat yang menghubungkan sudut persegi dengan lingkaran ( Fluent 6.0 User’s Guide , 2001). Panjang persegi adalah X dan panjang garis hubung adalah Y. Rasio X dan Y adalah 2. Tabel 1 menunjukkan studi kualitas mesh dengan memvariasikan X dan Y. Sesuai dengan tujuan simulasi maka mesh divalidasi dengan melihat pengaruh kerapatan mesh terhadap perubahan laju erosi. Tabel 2 menunjukkan parameter yang digunakan untuk validasi mesh.
Gambar 1. Domain komputasi
Tabel 1. Pengaturan grid penampang pipa Section Mesh 1  Mesh 2  Mesh 3 Mesh 4 X 8 12 14 16 Y 4 6 7 8
Gambar 2. Metode meshing penampang pipa
Tabel 2. Parameter validasi mesh
Parameter Fluida Air Kecepatan (m/s) 6,3 Laju alir massa cairan (kg/s) 12,7 Laju alir massa partikel (kg/s) 0,1285 Konsentrasi volume partikel (%) 0,37

## 2.2. Kondisi Simulasi Aliran Multifasa Cair-Gas-Partikel
Studi prediksi erosi aliran cair-gas-padat dalam belokan pipa ini menggunakan kerangka eulerian-lagrangian . Model eulerian digunakan untuk memecahkan persamaan aliran multifasa cair- gas  dan  kerangka lagrangian dengan  model  fasa  diskrit  ( Discrete  Phase  Model/DPM )  digunakan untuk menentukan lintasan partikel abrasif dalam aliran dan tumbukannya dengan dinding pipa yang digunakan dalam perhitungan laju erosi.
Model yang digunakan untuk mensimulasikan aliran multifasa adalah eulerian mixture 2 fasa  dengan  cairan  sebagai  fasa  primer  dan  gas  sebagai  fasa  sekunder.  Interaksi  antar  fasa dinyatakan  dalam  tegangan  permukaan.  Pasir  sebagai  fasa  diskrit  diasumsikan  tidak mempengaruhi profil aliran fluida karena pembebanan partikel yang rendah (ANSYS Workbench Help ). Simulasi  ini  menggunakan  model steady  state dengan  pertimbangan  bahwa  lintasan partikel  didalam  aliran  dari  masukan  hingga  keluaran  sudah  mencukupi  untuk  menghasilkan informasi  paramater  yang  digunakan  dalam  perhitungan  erosi. Slip  velocity digunakan  untuk mengakomodasi  perbedaan  kecepatan  antara  fasa  kontinyu  dan  fasa  terdispersi. Implicit  body formulation diterapkan  untuk  memperbaiki  konvergensi  solusi  dengan  memperhitungkan kesetimbangan  gradien  tekanan  dan body  force (gravitasi  dan  tegangan  permukaan)  dalam persamaan momentum. Model viscous k-  realizable dengan enhanced wall treatment digunakan untuk memperhitungkan faktor turbulensi didalam aliran dan mengangani viskositas didaerah dekat dinding ( near-wall region ).
Metode solusi yang diterapkan dalam simulasi menggunakan scheme SIMPLE, gradient menggunakan Least Square Cell Based , tekanan menggunakan PRESTO!, momentum dan volume fraction menggunakan QUICK, turbulence kinetic energy dan dissipation rate keduanya menggunakan Second Order Upwind untuk memperkecil galat perhitungan (ANSYS Workbench Help ). Untuk menentukan besaran dan profil erosi pada material target, ANSYS Fluent 18.2 menyediakan beberapa model erosi yang bisa digunakan, yaitu model Generic , Finnie, McLaury dan Oka. Fitur ini dapat ditemukan dengan mengaktifkan Erosion/Accretion didalam tab Physical Models setelah Discrete Phase Model diaktifkan pada tahap setup model. Pada studi ini, model erosi yang digunakan adalah Finnie, McLaury dan Oka. Parameter fluida, partikel dan material target yang dimasukkan ke dalam Fluent untuk keperluan validasi diperoleh berdasarkan informasi yang ditunjukkan pada Tabel 3. Adapun pada simulasi lanjutan untuk melihat pengaruh perubahan parameter fluida, partikel dan ukuran pipa terhadap laju dan profil erosi menggunakan data yang ditunjukkan pada Tabel 2 dan 3.

## 2.3. Pemodelan Aliran Fluida
Persamaan kontinuitas untuk campuran adalah:
𝜕 𝜕𝑡 (𝜌 𝑚 ) + ∇ ∙ (𝜌 𝑚 𝑣⃑ 𝑚 ) = 𝑚̇ (1)
Dimana 𝑣⃑ 𝑚 adalah rata-rata kecepatan massa.
𝑣⃑ 𝑚 = ∑ 𝛼 𝑘 𝑛 𝑘=1 𝜌 𝑘 𝑣 ′ 𝑘 𝜌 𝑚 (2)
Dengan 𝜌 𝑚 adalah kerapatan campuran.
𝜌 𝑚 = ∑ 𝛼 𝑘 𝑛 𝑘=1 𝜌 𝑘 (3)
𝛼 𝑘 adalah fraksi volume fasa k dan 𝑚̇ mengGambarkan perpindahan massa.
Persamaan momentum untuk campuran diperoleh dengan menjumlahkan persamaan momentum tiap fasa.
𝜕 𝜕𝑡 (𝜌 𝑚 𝑣⃑ 𝑚 ) + ∇ ∙ (𝜌 𝑚 𝑣⃑ 𝑚 𝑣⃑ 𝑚 ) = −∇𝑝 + ∇ ∙ [𝜇 𝑚 (∇𝑣⃑ 𝑚 + ∇𝑣⃑ 𝑚 𝑇 )] + 𝜌 𝑚 𝑔⃑ + 𝐹⃑ + ∇ ∙ (∑ 𝛼 𝑘 𝜌 𝑘 𝑣⃑ 𝑑𝑟,𝑘 𝑛 𝑘=1 𝑣⃑ 𝑑𝑟,𝑘 ) (4)
Dengan n adalah banyaknya fasa, 𝐹⃑ adalah body force, 𝜇 𝑚 adalah viskositas campuran dengan formula:
𝜇 𝑚 = ∑ 𝛼 𝑘 𝜇 𝑘 𝑛 𝑘=1 (5)
𝑣⃑ 𝑑𝑟,𝑘 adalah kecepatan drift untuk fasa sekunder k dengan formula:
𝑣⃑ 𝑑𝑟,𝑘 = 𝑣⃑ 𝑘 − 𝑣⃑ 𝑚 (6)
Persamaan perpindahan untuk k dan  pada model k-  adalah:
𝜕 𝜕𝑡 (𝜌𝑘) + 𝜕 𝜕𝑥 𝑖 (𝜌𝑘𝑢 𝑗 ) = 𝜕 𝜕𝑥 𝑖 [(𝜇 + 𝜇 𝑡 𝜎 𝑘 ) 𝜕𝑘 𝜕𝑥 𝑗 ] + 𝐺 𝑘 + 𝐺 𝑏 − 𝜌𝜀 − 𝑌 𝑀 + 𝑆 𝑘 (7)
dan 𝜕 𝜕𝑡 (𝜌𝜀) + 𝜕 𝜕𝑥 𝑖 (𝜌𝜀𝑢 𝑗 ) = 𝜕 𝜕𝑥 𝑗 [(𝜇 + 𝜇 𝑡 𝜎 𝜀 ) 𝜕𝜀 𝜕𝑥 𝑗 ] + 𝜌𝐶 1 𝑆𝜀 − 𝜌𝐶 2 𝜀 2 𝑘 + √𝑣𝜀 + 𝐶 1𝜀 𝜀 𝑘 𝐶 3𝜀 𝐺 𝑏 + 𝑆 𝜀 (8)
dimana:
𝐶 1 = 𝑚𝑎𝑥 [0.43, 𝜂 𝜂 + 5 ] (9)
dan 𝜂 = 𝑆 𝑘 𝜀 (10)
dengan: 𝐺 𝑘 = energi kinetik turbulen karena gradien kecepatan rata-rata. 𝐺 𝑏 = energi kinetik turbulen karena buoyancy. 𝑌 𝑀
= energi akibat fluktuasi dilatasi dalam turbulen compressible pada keseluruhan laju disipasi.
𝐶 2 , 𝐶 1𝜀 = konstanta. 𝜎 𝑘 , 𝜎 𝜀 = bilangan Prandtl turbulen untuk k dan  . 𝑆 𝑘 , 𝑆 𝜀 = parameter yang bisa didefinisikan pengguna ( user-defined ).

## 2.4. Pemodelan Gerak Partikel
Prediksi  lintasan  partikel  fasa  diskrit  yang  digunakan  oleh  ANSYS  Fluent  diperoleh dengan  mengintegralkan  gaya  kesetimbangan  ( force  balance )  partikel  yang  dituliskan  dalam kerangka lagrangian .  Persamaan  gaya  kestimbangan  yang  menunjukkan  inersia  partikel  akibat gaya yang bekerja pada partikel yang jika dituliskan pada koordinat-x adalah:
𝑑𝑢 𝑝 𝑑𝑡 = 𝐹 𝐷 (𝑢 − 𝑢 𝑝 ) + 𝑔 𝑥 (𝜌 𝑝 − 𝜌) 𝜌 𝑝 + 𝐹 𝑥 (11)
dimana 𝐹 𝐷 (𝑢 − 𝑢 𝑝 ) adalah gaya drag per unit massa partikel.
𝐹 𝐷 = 18𝜇 𝜌 𝑝 𝑑 𝑝 2 𝐶 𝐷 𝑅𝑒 24 (12)
dengan 𝑢 adalah  kecepatan  fasa  fluida, 𝑢 𝑝 adalah  kecepatan  partikel, 𝜇 adalah  viskositas molekular fluida, 𝜌 adalah kerapatan fluida, 𝜌 𝑝 adalah kerapatan partikel dan 𝑑 𝑝 adalah diameter partikel. 𝑅𝑒 adalah bilangan Reynolds relatif yang didefinisikan sebagai:
𝑅𝑒 ≡ 𝜌𝑑 𝑝 |𝑢 𝑝 − 𝑢| 𝜇 (13)
Haider dan Levenspiel [6] mengajukan persamaan koefisien drag yang memperhitungkan faktor bentuk partikel (  ), melalui persamaan:
𝐶 𝐷 = 24 𝑅𝑒 (1 + 𝑏 1 𝑅𝑒 𝑏2 ) + 𝑏 3 𝑅𝑒 𝑏 4 + 𝑅𝑒 (14)
dimana: 𝑏 1 = exp⁡(2,3288 − 6,4581∅ + 2,4486∅ 2 ) 𝑏 2 = 0,0964 + 0,5565∅ 𝑏 3 = exp(4,905 − 13,8944∅ + 18,4222∅ 2 − 10,2599∅ 3 ) 𝑏 4 = exp(1,4681 + 12,2584∅ − 20,7322∅ 2 + 15,8855∅ 3 )
Faktor bentuk ( ∅ ) didefinisikan sebagai:
∅ = 𝑠 𝑆 (15)
dimana s adalah luas permukaan bulatan yang memiliki volume yang sama dengan partikel, dan S adalah luas permukaan partikel yang sebenarnya.
Persamaan  lintasan  dan  persamaan  tambahan  yang  mengGambarkan  perpindahan massa atau panas  dari  dan ke partikel  dipecahkan  mengintegralkan  secara  bertahap  terhadap waktu diskrit. Jika persamaan (11) diintegralkan terhadap waktu menghasilkan kecepatan partikel ditiap titik sepanjang lintasan dengan lintasan sendiri yang diprediksi sebagai:
𝑑𝑥 𝑑𝑡 = 𝑢 𝑝 (16)
Jika diasumsikan gaya yang bekerja pada partikel adalah konstan, maka persamaan lintasan dapat dituliskan dalam bentuk yang lebih sederhana:
𝑑𝑢 𝑝 𝑑𝑡 = 1 𝜏 𝑝 (𝑢 − 𝑢 𝑝 ) (17)
dimana 𝜏 𝑝 adalah waktu relaksasi partikel.
Fluent menggunakan trapezoidal scheme untuk mengintegralkan persamaan (17):
𝑢 𝑝 𝑛+1 − 𝑢 𝑝 𝑛 ∆𝑡 = 1 𝜏 (𝑢 ∗ − 𝑢 𝑝 𝑛+1 ) + ⋯ (18)
dimana n menunjukkan banyaknya iterasi, dan
𝑢 ∗ = 1 2 (𝑢 𝑛 + 𝑢 𝑛+1 ) (19)
𝑢 𝑛+1 = 𝑢 𝑛 + ∆𝑡𝑢 𝑝 𝑛 ∙ ∇𝑢 𝑛 (20)
Persamaan (16) dan (17) dipecahkan secara simultan untuk menentukan kecepatan dan posisi partikel terhadap waktu.

## 2.5. Model Erosi

## 2.5.1. Model Erosi Finnie
Berdasarkan  beberapa  asumsi,  Finnie  mengajukan  dua  persamaan  untuk memperkirakan volume material yang hilang oleh butiran abrasif tunggal untuk sudut tumbukan rendah dan tinggi sebagai berikut:
𝜀 𝑉𝑃 = 𝑚𝑝𝑉𝑝2 𝑝  (sin(2𝛼) − 6  sin 2 (𝛼)) untuk tan 𝛼 ≤  6 (21)
𝜀 𝑉𝑃 = 𝑚𝑝𝑉𝑝2 𝑝  (  cos2 𝛼 6 ) untuk tan 𝛼 ≥  6 (22)
Dimana 𝜀 𝑉𝑃 adalah volume material permukaan  yang  dipindahkan, 𝑚 𝑝 adalah  massa partikel, Vp adalah kecepatan partikel, p adalah konstanta tegangan aliran plastis,  adalah sudut tumbukan partikel,  adalah rasio komponen  gaya vertikal pada permukaan partikel terhadap komponen gaya horizontal.  adalah rasio kedalaman kontak l terhadap y t (sesuai Gambar 3) memiliki nilai konstan, Finnie menggunakan nilai 2 untuk parameter ini. Persamaan pertama berlaku  untuk  sudut  tumbukan  rendah  yang  berkaitan  dengan  keadaan  dimana  partikel meninggalkan permukaan sambil memotong. Persamaan kedua berlaku untuk sudut tumbukan tinggi pada keadaan dimana gerak horizontal ujung partikel terhenti saat memotong.
Gambar 3. Abrasif yang menumbuk permukaan dan memotong material. Vektor kecepatan awal membentuk sudut  terhadap pusat gravitasi partikel (Finnie, 1960).

## 2.5.2. Model Erosi McLaury
McLaury  mengajukan  model  prediksi  laju  erosi  partikel  dengan  fluida  pembawa  air. Model ini dikembangkan untuk mensimulasikan laju erosi pada aliran cair-padat ( slurry) . Model McLaury dituliskan sebagai berikut (Oka et al, 2005).
𝐸 = 𝐴⁡𝑉 𝑛 ⁡𝑓(𝛾) (23) 𝐴 = 𝐹⁡𝐵ℎ 𝑘 (24) Dimana: F = konstanta empiris V = kecepatan tumbukan partikel Bh
= angka kekerasan Brinell material dinding
k
= -0,59 untuk carbon steel (nilai k berbeda untuk material lain)
dengan:
𝑓(𝛾) = 𝑏𝛾 2 + 𝑐𝛾⁡⁡⁡⁡⁡⁡⁡⁡⁡𝑢𝑛𝑡𝑢𝑘⁡𝛾 ≤ 𝛾 𝑙𝑖𝑚 (25)
𝑓(𝛾) = 𝑥⁡ cos 2 𝛾⁡sin(𝑤𝛾) + 𝑦⁡ sin 2 (𝛾) + 𝑧⁡⁡⁡⁡⁡𝑢𝑛𝑡𝑢𝑘⁡𝛾 > 𝛾 𝑙𝑖𝑚 (26)
Nilai konstanta b, c, w, x dan y pada model diatas ditentukan berdasarkan eksperimen. Nilai z ditentukan sedemikian sehingga kedua persamaan fungsi sudut diatas berdekatan pada 𝛾 = 𝛾 𝑙𝑖𝑚 , dengan 𝛾 𝑙𝑖𝑚 adalah sudut transisi.

## 2.5.3. Model Erosi Oka
Model erosi Oka dikembangkan berdasarkan pengujian yang dilakukan pada beberapa jenis  material  target seperti  besi  murni,  aluminum,  baja karbon  dan baja  stainless.  Laju  erosi akibat  tumbukan  partikel  dengan  material  target pada sudut  tumbukan 𝛼 dinyatakan  sebagai (Oka et al, 2005;
Sedrez et al, 2019): 𝐸(𝛼) = 𝑔(𝛼)⁡𝐸 90 (27)
𝐸 90 adalah laju erosi pada sudut tumbukan 90 0 dan 𝑔(𝛼) adalah fungsi sudut tumbukan yang dinormalkan, dinyatakan oleh dua fungsi trigonometri dan angka kekerasan awal material 𝐻𝑣 (dalam GPa), sesuai persamaan berikut:
𝑔(𝛼) = (sin 𝛼) 𝑛1 (1 + 𝐻𝑣(1 − sin 𝛼)) 𝑛2
(28)
𝑛1 dan 𝑛2 adalah  eksponen  yang  ditentukan  oleh  kekerasan  material  dan  kondisi tumbukan  lainnya  seperti  properti  partikel  termasuk  bentuk  partikel. 𝑔(𝛼) ditentukan berdasarkan (Oka et al, 2005) Representasi grafis 𝑔(𝛼) ditunjukkan seperti Gambar 4.
Gambar 4. Konsep erosi akibat deformasi plastis berulang dan pemotongan [9].
Suku pertama  dari persamaan  (28)  mengGambarkan  deformasi  plastis  berulang  atau karakteristik getas dan dengan komponen vertikal dari energi tumbukan partikel jika 𝑛1 = 2. Suku kedua menunjukkan fenomena pemotongan, dimana erosi relatif dan lebih efektif terjadi pada sudut tumbukan yang lebih rendah. 𝐸 90 dinyatakan sebagai:
𝐸 90 = 𝐾⁡(𝐻𝑣) 𝑘1 (𝑣) 𝑘2 (𝐷) 𝑘3 (29)
dimana k1 , k2 dan k3 adalah faktor eksponen yang dipengaruhi oleh parameter lain. K adalah sembarang unit yang menunjukkan properti partikel seperti bentuk dan tingkat kekerasan, yang nilainya  berbeda-beda  untuk  jenis  partikel  berbeda.  ANSYS  Fluent  sudah  memasukkan  nilai konstanta-konstanta empiris tersebut ke dalam masing-masing model erosi, sehingga pengguna hanya memasukkan parameter fluida, partikel dan material target.

## 2.6. Prosedur Simulasi
Secara garis besar langkah-langkah simulasi dalam studi ini secara adalah:
a. Membuat geometri sebagai domain komputasi.
b. Melakukan studi kualitas mesh.
c. Pemodelan aliran multifasa dan model fasa diskrit sesuai dengan kondisi eksperimen.
d. Validasi simulasi berdasarkan hasil eksperimen.
e. Variasi parameter aliran fluida dan partikel.
Geometri  komputasi  yang  digunakan  untuk  validasi  simulasi  ini  berdasarkan eksperimen  yang  dilakukan  oleh  Sedrez seperti  ditunjukkan pada  Gambar  1,  dengan kondisi eksperimen  sesuai  data  yang  ditunjukkan  pada  Tabel  3.  Eksperimen  yang  dilakukan  Sedrez menghasilkan laju erosi 0,12, 0,18, 0,26 dan 0,40 mm/day untuk masing-masing kasus 1, 2, 3 dan 4.
Tabel 3. Kondisi eksperimen Sedrez (Sedrez, 2019) Diatemer partikel (μm) 300 Kerapatan partikel (kg/m 3 ) 2650 Kerapatan material pipa (kg/m 3 ) 7800 Parameter Kasus 1 Kasus 2 Kasus 3 Kasus 4 Fluida air air dan udara air dan udara air dan udara Kecepatan (m/s) 6,3 6 dan 2,4 5,5 dan 3,76 6,31 dan 5,46 Laju alir massa cairan dan gas  (kg/s) 12,7 12,1 dan 0,00505 11,2 dan 0,00933 12,8 dan 0,01356 Laju alir massa partikel (kg/s) 0,1285 0,1150 0,0914 0,1131 Konsentrasi volume partikel (%) 0,37 0,36 0,31 0,33 3. HASIL DAN PEMBAHASAN

## 3.1. Validasi Mesh dan Validasi Simulasi
Dengan menggunakan prosedur pengaturan mesh sesuai Tabel 1 dan parameter sesuai Tabel  2  diperoleh  hasil  seperti  ditunjukkan  pada  Gambar  5.  Dari  Gambar  5  terlihat  bahwa deviasi laju erosi maksimum (ER) antara mesh 4 dan 3 sudah sangat kecil, yaitu 2,2%, 0,4% dan 8,3% masing-masing untuk model Finnie, McLaury dan Oka. Deviasi antara mesh 4 dan 2 untuk model Finnie, McLaury dan Oka masing-masing adalah 20%, 15,9% dan 24,9%. Mesh 2 dapat digunakan untuk simulasi berikutnya dengan pertimbangan bahwa perbedaan 20% bisa diterima untuk simulasi prediksi erosi ( Fluent 6.0 User’s Guide , 2001).
Gambar 5. Grafik uji mesh independency
Hasil validasi simulasi ditunjukkan pada Tabel 4 dengan representasi grafis ditunjukkan pada  Gambar  5.  Terlihat  bahwa  bahwa  model  erosi  Oka  lebih  mendekati  hasil  eksperimen dengan perbedaan relatif 0,3% pada kasus 1 dan 0,8% pada kasus 4. Perbedaan relatif terbesar 24% pada kasus 2, diikuti perbedaan relatif 21,6% pada kasus 3. Model Finnie dan McLaury memberikan perbedaan relatif yang cukup besar pada semua kasus yaitu lebih dari 100% seperti ditunjukkan pada Tabel 4. Namun demikian terlihat bahwa pada dasarnya ketiga model erosi tersebut menunjukkan kecenderungan hasil yang hampir sama. Kecenderungan yang sama pada
model Finnie, McLaury dan Oka menunjukkan bahwa ketiga model memiliki pendekatan yang sama  dalam  menentukan  parameter  yang  mempengaruhi  erosi.  Ketiganya  menggunakan kecepatan partikel, fungsi sudut dan konstanta yang tergantung dari tingkat kekerasan material baja karbon dalam dasar pemodelannya.
Tabel 4. Perbandingan simulasi model erosi dengan hasil eksperimen
Kasus ER (mm/day) Finnie McLaury Oka Eksperimen Kasus 1 0,3910 0,9072 0,1196 0,12 Kasus 2 0,4187 0,9249 0,1362 0,18 Kasus 3 0,5782 1,2406 0,2038 0,26 Kasus 4 1,0889 2,2043 0,4032 0,40
Gambar 3.2. Perbandingan hasil simulasi dan eksperimen
Setelah  simulasi  sudah  divalidasi  maka  variasi  dilakukan  untuk  melihat  pengaruh perubahan parameter fluida dan partikel terhadap laju erosi maksimum dan profil erosi yang terjadi.  Simulasi  menggunakan  fluida  minyak  mentah,  gas  alam  dan  partikel  pasir  dengan parameter sebagaimana ditunjukkan pada Tabel 5. Skenario variasi yang dilakukan pada studi ini ditunjukkan pada Tabel 6.
Tabel 5. Parameter fluida dan partikel dalam simulasi
Fluida Kerapatan (kg/m3) Viskositas (kg/m-s) Tegangan Permukaan (N/m) Minyak hidrokarbon 879,6 3,839 x 10-3 - Air 998,2 1,003 x 10-3 - Gas Alam 0,762 1,300 x 10-5 0.072 Minyak mentah water cut 90% 986,34 1,287 x 10-3 Partikel Diameter 30 – 300 µm Kerapatan 2.224 kg/m 3
Tabel 6. Skenario variasi parameter dalam simulasi Parameter Data Masukan Setara Diameter pipa (mm) 50,8;  76,2; 101,16 ( 2, 3 dan 4 inch) Laju  alir  minyak  mentah (kg/s) 7,260 - 14,520 (4000 – 8000 bbls/day) Laju alir gas (kg/s) 0,0025 – 0,075 (0,1 – 1 MMSCFD) Laju alir partikel (kg/s) 0,000926 – 0,001852 2.60. – 100 kg/day)

## 3.2. Pengaruh Laju Alir Massa Minyak Mentah Terhadap Laju Erosi
Pada  diameter  pipa  tetap  dengan  laju  alir  gas  dan  diameter  partikel  tetap,  laju  erosi maksimum  meningkat  dengan  kenaikan  laju  alir  minyak  mentah  seperti  ditunjukkan  pada Gambar 6.
Gambar 6. Pengaruh laju alir minyak mentah terhadap laju erosi
Dari  hasil  simulasi  dapat  dilihat  bahwa  laju  erosi  maksimum  meningkat  dengan kenaikan  laju  alir minyak mentah.  Keempat  model  menunjukkan kecenderungan  yang  sama. Sebagaimana hasil yang diperoleh pada tahap validasi, model erosi Oka memberikan hasil yang paling kecil dibandingkan dua model lainnya, sementara model McLaury memberikan hasil laju erosi yang paling besar. Hasil ini sesuai dengan teori bahwa laju erosi berbanding lurus dengan kecepatan.  Dalam  hal  ini  peningkatan  laju  alir  fasa  cair  berarti  peningkatan  kecepatan  aliran fluida yang pada akhirnya meningkatkan kecepatan partikel yang menumbuk dinding material target.  Kondisi ini  sesuai  dengan hasil penelitian  Zhang  et al.  yang menunjukkan bahwa  laju erosi  meningkat  dengan  kenaikan  kecepatan  fluida  cair. Pada  Gambar 7  terlihat bahwa  erosi sebagian besar terdistribusi ke bagian luar belokan dengan kenaikan laju alir fasa cair. Hal ini menunjukkan bahwa fluida dengan kecepatan yang lebih tinggi bergerak dari inti (sumbu pipa vertikal) aliran ke arah punggung belokan.  Studi eksperimental aliran dua fasa air-udara yang dilakukan Yadav menunjukkan bahwa fluida di daerah belokan mengalami gaya sentripetal yang mengarahkan aliran dari pusat ke bagian luar belokan.
0.0000 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 3000 4000 5000 6000 7000 8000 E R (m m /da y) Laju alir minyak mentah (BPD) Finnie McLaury Oka
Gambar 7. Profil erosi pada laju alir fasa cair (a) 7,26 kg/s (b) 9,075 kg/s (c) 10,890 kg/s (d) 12,705 kg/s (e) 14,520 kg/s.

## 3.3. Pengaruh Laju Alir Massa Gas Terhadap Laju Erosi
Pada diameter pipa tetap dengan laju alir minyak mentah dan diameter partikel tetap, laju erosi maksimum meningkat dengan kenaikan laju alir gas seperti ditunjukkan pada Gambar 8.  Hal  ini  disebabkan  dengan  meningkatnya  laju  alir  gas  maka  kecepatan  fluida  semakin meningkat sesuai dengan penjumlahan kecepatan superfisial masing-masing fasa didalam fluida.
Gambar 8. Pengaruh laju alir fasa gas terhadap laju erosi
Dari Gambar 9 dapat dilihat bahwa pada laju alir gas rendah, erosi maksimum terjadi di belokan  bagian  samping  dan  semakin  bergeser  ke  tengah  seiring  peningkatan  laju  alir  gas. Disamping  itu  peningkatan  aliran  gas  cenderung  menyebabkan  lebih  banyak  partikel  yang menumbuk dinding pipa. Keadaan ini sesuai dengan penelitian yang dilakukan oleh Xu et al. bahwa  kenaikan  kecepatan  gas  menyebabkan  kapasitas  pengangkutan  partikel  oleh  gas meningkat dan meningkatkan energi kinetik partikel sehingga tekanan sentuh lebih besar pada dinding dalam pipa.
0.00 0.50 1.00 1.50 2.00 2.50 0 0.02 0.04 0.06 0.08 E R (m m /day ) Laju alir massa gas (kg/s) Finnie McLaury Oka
Gambar 9. Profil erosi pada laju alir gas (a) 0,0025 kg/s; (b) 0,0125 kg/s; (c) 0,025 kg/s;
(d) 0,05 kg/s; (e) 0,075 kg/s.
Perbedaan besaran laju erosi maksimum antara ketiga model dengan kenaikan laju alir fluida dapat dilihat dari beberapa hal. Pertama adalah bagaimana perbedaan ketiga model dalam menentukan  eksponen  kecepatan,  Finnie  menggunakan  nilai  2,  McLaury  menggunakan  nilai nilai 1,73, sementara Oka meggunakan eksponen 0,19. Hal ini yang memungkinkan laju Erosi yang dihasilkan oleh model Oka paling kecil dibandingkan model lainnya.  Kedua, penentuan konstanta  empiris  pada  ketiga  model  menggunakan  pendekatan  yang  berbeda  sesuai  dengan kondisi eksperimen masing-masing model. Ketiga, acuan sudut tumbukan kritis yang digunakan ketiga model juga berbeda, Finnie mengambil sudut kritis pada 18,5 0 , McLaury menggunakan sudut 15 0 , sementara Oka secara tidak langsung menunjukkan sudut sekitar 40 0 sebagai sudut tumbukan maksimum, seperti ditunjukkan pada Gambar 4.
Dari  grafik  pengaruh  laju  alir  fluida  terhadap  laju  erosi  maksimum,  terlihat  bahwa semakin tinggi laju alir (kecepatan) menghasilkan perbedaan relatif laju erosi maksimum yang semakin besar antar model erosi. Hal ini mengindikasikan bahwa penerapan model erosi dalam suatu simulasi sangat bergantung kepada rentang kecepatan aliran yang digunakan.

## 3.4. Pengaruh Diameter Pipa Terhadap Laju Erosi
Hasil simulasi dengan variasi diameter pipa ditunjukkan pada Gambar 10. Dari Gambar tersebut dapat dilihat bahwa semakin besar diameter pipa maka semakin kecil laju erosi yang terjadi. Hal ini disebabkan dengan bertambahnya diameter pipa maka kecepatan aliran fluida semakin rendah dan laju erosi semakin rendah sebanding dengan eksponen kecepatan masing- masing  model.  Xu  et  al.  mengobservasi  dari  penelitiannya  bahwa  kenaikan  diamater  pipa meningkatkan  dispersi  partikel  sehingga  tumbukan  partikel  pada  dinding  dalam  pipa  juga terdispersi, dengan demikian laju erosi maksimum menurun.
Gambar 10. Pengaruh diameter pipa terhadap laju erosi
Gambar 11. Profil erosi pada pipa diameter (a) 2 inch (b) 3 inch (c) 4 inch.
Gambar 11 menunjukkan bahwa distribusi erosi pada pipa 2 inch hampir sama dengan pipa 3 inch yaitu terdistribusi dominan pada punggung belokan, yang membedakan adalah pada besaran laju erosi. Pada pipa berdiameter 4 inch selain menurunnya laju erosi, profil erosi tidak lagi  pada  punggung  belokan  melainkan  bergeser  ke  bawah.  Hal  ini  dimungkinkan  dengan meningkatnya  diameter  pipa  kecepatan  aliran  menurun  sehingga  gaya  sentripetal  yang menyebabkan aliran dari daerah inti mengarah ke dinding belokan bagian luar juga menurun. Penurunan dua kali kecepatan setara dengan penurunan empat kali gaya sentripetal fluida. Selain itu penurunan kecepatan menyebabkan bilangan Stokes menurun, artinya dominasi aliran fluida terhadap gerakan partikel tidak terlalu besar, sehingga inersia partikel cenderung mengarahkan partikel ke bawah.

## 3.5. Pengaruh Laju Alir Massa Partikel Terhadap Laju Erosi
Hasil  simulasi  variasi  laju  alir  partikel  abrasif  ditunjukkan  pada  Gambar  12.  Dari Gambar terlihat bahwa laju erosi bertambah dengan kenaikan laju alir massa partikel. Kondisi ini disebabkan dengan meningkatnya laju alir massa partikel maka semakin banyak partikel yang
menumbuk dinding dalam pipa sehingga meningkatkan laju erosi. Fenomena peningkatan erosi akibat  meningkatnya  laju  alir  massa  partikel  memiliki  kecenderungan  yang  sama  dengan kenaikan laju alir fluida cair, sebagaimana ditunjukkan pada Gambar 6.
Gambar 12. Pengaruh laju alir partikel terhadap laju erosi
Gambar 13 memperlihatkan profil erosi dengan variasi laju alir massa fluida 60 – 100 kg/day. Pada laju alir partikel rendah, besaran laju erosi terlihat rendah dan distribusi erosi yang tidak  terlalu  banyak  pada  daerah  punggung  belokan.  Laju  erosi  meningkat  seiring  dengan meningkatnya distribusi erosi didaerah punggung belokan ketika laju alir massa partikel mulai dari 70 hingga 100 kg/day dengan pola yang sama dengan profil erosi pada variasi laju alir fasa cair.  Di  dalam  ANSYS  Fluent,  laju  alir  massa  partikel  digunakan  untuk  menentukan  massa partikel yang kemudian akan dimasukkan kedalam perhitungan model erosi. Karena itu laju erosi dengan peningkatan massa partikel menghasilkan kenaikan yang linier, sedangkan peningkatan kecepatan fluida menghasilkan kecenderungan laju erosi yang naik secara eksponensial sesuai dengan konstanta eksponen yang digunakan didalam masing-masing model erosi.
Gambar 13. Profil erosi pada laju alir partikel (a) 60kg/day (b) 70kg/day (c) 80kg/day (d) 90kg/day (e)

## 100kg/day.

## 3.6. Pengaruh Ukuran Partikel Terhadap Laju Erosi
Hasil  simulasi  variasi  ukuran  partikel  ditunjukkan  pada  Gambar  14.  Dari  Gambar tersebut  terlihat  bahwa  laju  erosi  maksimum  meningkat  dengan  semakin  besarnya  ukuran
0.00 0.01 0.01 0.02 0.02 0.03 0.03 0.04 0.04 0.05 50 60 70 80 90 100 E R (m m /day ) Laju alir massa partikel (kg/s) Finnie McLaury Oka
(a)
partikel.  Hal  ini  sesuai  dengan  penelitian  yang  dilakukan  Desale  yang  menyimpulkan  bahwa peningkatan ukuran partikel menyebabkan meningkatnya laju erosi. Kenaikan diameter partikel berarti  meningkatnya  massa  partikel,  yang  pada  kecepatan  tetap  memberikan  energi  kinetik tumbukan  terhadap  material  target  lebih  besar.  Hasil  simulasi  juga  didukung  oleh  hasil eksperimen  yang  dilakukan  Okita  yang  menunjukkan  bahwa  semakin  besar  ukuran  partikel cenderung meningkatkan rasio erosi pada sembarang nilai viskositas.
Gambar 14. Pengaruh ukuran partikel terhadap laju erosi.
Gambar 15. Profil erosi pada variasi ukuran partikel
(a) 30  m; (b) 100  m; (c) 200  m; (d) 300  m.
Gambar  15  menunjukkan  distribusi  erosi  yang  hampir  sama  pada  variasi  diameter partikel  30 – 300  m  dengan  dominasi  erosi  terbesar  pada  daerah  punggung  belokan,  yang membedakan adalah pada intensitas laju erosinya.  Pada diameter 200  m, terlihat konsentrasi laju erosi tinggi lebih padat dibandingkan pada diameter partikel 300  m.  Hal ini dimungkinkan bahwa  pada  partikel  berdiameter  lebih  besar  atau  lebih  berat,  kecenderungan  partikel  untuk turun  dari streamline fluida  pembawa  lebih  besar.  Peningkatan  diamater  partikel  juga  berarti meningkatkan  bilangan  Stokes  yang  berbanding  lurus  dengan  kuadrat  diameter  partikel, sehingga kecenderungan partikel untuk tidak didominasi oleh aliran fluida semakin besar. Dapat dikatakan bahwa pada  diameter  partikel 300  m,  sebagian partikel  terbawa  oleh aliran fluida hingga menumbuk dinding dan sebagian lain keluar dari streamline fluida yang  mengarah ke dinding pipa.
0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0 100 200 300 E R (m m /d ay ) Ukuran partikel (mikron) Finnie McLaury Oka a)

## 4. KESIMPULAN
Pada  studi  ini  prediksi  erosi  dengan  menggunakan  CFD  dilakukan  dengan menggunakan  tiga model  erosi  yaitu model Finnie,  McLaury  dan  Oka.  Dari  semua  skenario simulasi yang dibuat, model McLaury menghasilkan laju erosi yang paling besar dibandingkan dua model  lainnya.  Model Oka secara konsisten menghasilkan laju  erosi  yang paling  rendah. Perbedaan  ini  disebabkan  ketiga  model  erosi  dibangun  dari  pendekatan  mekanistik  maupun eksperimen yang berbeda, sehingga dimungkinkan masing-masing model erosi hanya berlaku pada kondisi tertentu saja. Dalam studi ini model Oka lebih unggul dibandingkan model Finnie dan McLaury dimana model Oka menghasilkan perbedaan relatif yang paling rendah terhadap eksperimen yang dijadikan acuan dalam validasi.
Hasil simulasi lanjutan menunjukkan bahwa laju erosi meningkat dengan peningkatan laju alir fasa cair, laju alir fasa gas, laju alir partikel, ukuran partikel dan laju erosi menurun dengan kenaikan diameter pipa.
Pada semua variasi laju alir fasa cair, profil erosi relatif sama yaitu kerapatan erosi terjadi pada punggung belokan akibat gaya sentripetal yang mengarahkan fluida dari inti pipa ke bagian luar belokan. Pada kenaikan laju alir gas, profil erosi terdistribusi di punggung belokan hingga bagian pipa keluaran yang disebabkan oleh peningkatan kapasitas pengangkutan partikel oleh fluida gas. Kenaikan ukuran pipa menyebabkan kecepatan aliran menurun dan gaya sentripetal menurun  sehingga  tumbukan  partikel  pada  punggung  belokan  berkurang.  Kenaikan  laju  alir partikel berarti peningkatan jumlah partikel yang menumbuk dinding pipa sehingga laju erosi meningkat.  Lintasan  partikel  berdiameter  kecil  cenderung  didominasi  aliran  fluida  sehingga menghasilkan profil erosi terdistribusi pada bagian samping maupun punggung belokan sesuai dengan  turbulensi  didaerah  belokan.  Sementara  semakin  besar  diameter  partikel  terdapat kecenderungan  sebagian  partikel  keluar  dari streamline fluida  pembawa  saat  masuk  ke  daerah belokan  karena  dominasi  fluida  terhadap  partikel  semakin  kecil  (bilangan  Stokes  meningkat dengan  kenaikan  diameter  partikel)  sehingga  semakin  sedikit  yang  menumbuk  punggung belokan.

## 5. DAFTAR PUSTAKA
ANSYS Fluent 18.2. ANSYS Workbench Help .
Arabnejad H., A.Mansouri, S.A. Shirazi, B.S. McLaury, 2015, Evaluation of solid particle erosion equations and models for oil and gas industry applications, Journal Paper of SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers , SPE-174987-MS. Chen  X.,  B.S.  McLaury,  S.A.  Shirazi,  2004,  Application  and  Experimental  Validation  of  a Computational Fluid Dynamics (CFD)-based Erosion Prediction Model in Elbows and Plugged-Tees, Journal of Computers & Fluids , Vol. 33, pp. 1251-1272.
Desale,  G.R.,  Gandhi,  B.K.,  Jain,  S.C,  2009,  Particle  Size  Effects  on  The  Slurry  Erosion  of Aluminum Alloy (AA 6063). Journal of Wear , Vol. 266, Page 1066-1071.
Edwards J.K. B.S McLaury, S.A. Shirazi, 2001, Modeling solid particle erosion in elbows and plugged tees. ASME Journal of Energy Resources Technology, Vol. 123, Page 277-284. Finnie I, 1960,  Erosion of Surface by Solid Particles, Journal of Wear , Vol. 3, Page 87-103. Fluent 6.0 User’s Guide, 2001.
Oka Y.I., K.Okamura, T.Yoshida, 2005, Practical estimation of erosion damage caused by solid particle impact Part 1: Effects of impact parameters on a predictive equation, Journal of Wear , Vol. 259, Page 95-101.
Oka Y.I., K.Okamura, T.Yoshida, 2005, Practical estimation of erosion damage caused by solid particle  impact  Part  2:  Mechanical  Properties  of  Materials  Directly  Associated  with Erosion Damage, Journal of Wear , Vol. 259, Page 102-109.
Okita,  R.,  Y.  Zhang,  B.S.  McLaury,  S.A.  Shirazi,  2012,  Experimental  and  Computational Investigation  to  Evaluate  the  Effect  of  Fluid  Viscosity  and  Particle  Size  on  Erosion Damage. Journal of Fluids Engineering , Vol. 134.
Parsi M., K. Najmi, F. Najafifard, S. Hassani, B.S. McLaury, S.A. Shirazi, 2014, A comprehensive review of solid particle erosion modeling for oil and gas wells and pipelines applications, Journal of Natural Gas Science and Engineering , Vol. 21, pp. 850-873.
Sedrez, T.A., S.A. Shirazi., Y.R. Rajkumar, K. Sambath, H.J. Subramani, 2019, Experiments and CFD  Simulations  of  Erosion  of  a  90 0 elbow  in  liquid-dominated  liquid-solid  and dispersed-bubble-solid flows. Journal of Wear , 426-427. Page 570-580.
Vieira, R.E., A. Mansouri, B.S. McLaury, S.A. Shirazi, 2016, Experimental and computational study of erosion in elbows due to sand particles in air flow, Journal of Powder Technology , Vol. 288, pp. 339-353.
Xu, J., Z. Lian, J. Hu, M. Luo, 2018, Prediction of Maximum Erosion Rate of Gas-Solid Two- Phase Flow Pipelines. Journal of Energies .
Yadav, M. S., 2013, Interfacial Area Transport Across Vertical Elbows in Air-Water Two Phase Flow, Dissertation ,  Department  of  Mechanical  and  Nuclear  Engineering,  The Pennsylvania State University.
Zhang,  Y.,  E.P.  Reuterfors,  B.S.  McLaury,  S.A.  Shirazi,  E.F.  Rybicki,  2007,  Comparison  of Computed and Measured Particle Velocities and Erosion in Water and Air Flows, Journal of Wear , Vol. 263, pp. 330-338.