Vol. 14, No. 1,  77-86,  Juli  2017
Metode Elemen Batas (MEB) Untuk Solusi Numerik Masalah  Statik Dari Material Elastis Isotropik Tak-Homogen
Mohammad Ivan Azis †

## Abstract
A  boundary  element  method  is  derived  for the  solution  of  static  boundary  value problems for inhomogeneous isotropic elastic materials. Some particular problems are considered to illustrate the application of the method.
Keywords: Boundary  element  method,  static,  boundary  value  problem, inhomogeneous, isotropic, elastic.

## 1.  Pendahuluan
Pendekatan  persamaan  integral  untuk  masalah  nilai  batas  elastostatik  lama  telah diperkenalkan oleh Rizzo (1967). Sejak saat itu sudah banyak penulis yang ikut menggunakan metode  elemen  batas  untuk  menentukan  solusi  numerik  dari  berbagai  masalah  statik  untuk material elastis, isotropik dan homogen (Brebbia & Dominguez, 1989).
Dalam  pertentangan,  aplikasi  dari  metode  elemen  batas  ke  masalah  untuk  material elastis  isotropik  tak-homogen  sangat  terbatas  jumlahnya  dan  hal  ini  disebabkan  oleh  karena kesulitan dalam memperoleh fungsi Green sebagai kernel dari persamaan integral batas yang bersesuaian.  Baru-baru  ini  Manolis  &  Shaw  (1996)  menurunkan  suatu  fungsi  Green  dari persamaan gelombang vektor untuk media isotropik heterogen lemah ( mildly heterogeneous ). Fungsi Green ini diturunkan untuk material dengan parameter yang bervariasi secara tertentu dan  secara  khusus  terbatas  untuk  kasus  dimana  parameter  Lamé  dan  sama.  Hal  ini berarti nilai rasio Poison adalah 0, 25 yang membatasi aplikasi. Akan tetapi, menurut Manolis & Shaw (1996) nilai rasio Poison ini adalah sesuatu biasa untuk material batuan (Turcotte & Schubert, 1982).
Tulisan ini mengembangkan kerja Manolis & Shaw (1996) untuk membangun suatu prosedur  perturbasi  untuk  menentukan  solusi  masalah  statik  dua  dimensi  untuk  media isotropik tak-homogen dengan parameter Lamé
) ( ) ( ) ( ) ( ) ( ) ( ) 1 ( 0 ) 1 ( 0 x x g x x x g x           dimana ) , , ( 3 2 1 x x x x  adalah  suatu  vektor  di 3 R , ) ( x g suatu  fungsi  yang  memenuhi batasan tertentu, ) 0 ( ) 0 (    adalah konstanta dan  suatu parameter kecil. Dalam batasan ini,    bentuk  parameter  Lamé  di  atas  menawarkan  pilihan  keberagaman  parameter ) ( x  dan ) ( x  yang cukup luas.

## 2. Persamaan Dasar
Dengan  merujuk  pada  kerangka  Cartesian 3 2 1 x x Ox persamaan  kesetimbangan  ( the equilibrium  equations )  dalam  suatu  material  elastis  tanpa  gaya  badan  ( body  force )  dapat ditulis dalam bentuk
† Staf Pengajar pada Jurusan Matematika FMIPA Universitas Hasanuddin Makassar
Mohammad Ivan Azis
0 ,  j ij  (1)
dimana ij  untuk 3 , 2 , 1 ,  j i melambangkan  tensor  stress,  tanda  koma  terindeks mengindikasikan  turunan  parsial  terhadap  kordinat  spasial j x dan  konvensi  penjumlahan indeks berulang (jumlahan dari 1 sampai 3) diterapkan. Hubungan displacement-stress dapat dituliskan sebagai
) )( ( ) ( , , , i j j i k k ij ij u u x u x        (2)
dimana k u untuk 3 , 2 , 1  k menyatakan  displacement  dan ij  adalah  operator  delta Kronecker. ) ( x  dan ) ( x  dengan ) , , ( 3 2 1 x x x x  dalam  persamaan  (2)  adalah  parameter Lamé  yang  diasumsikan  fungsi  yang  punya  turunan  kedua  ( twice  differentiable ).  Subtitusi persamaan (2) ke (1) menghasilkan
    0 ) ( ) ( , , ,    j i j j i k k ij x u x      (3)

## 3. Masalah Nilai Batas
Suatu  material  isotropik  tak-homogen  berada  pada  daerah  dalam  ruang  tiga dimensi 3 R dengan  batas   yang  terdiri  atas  sejumlah  berhingga  permukaan  tertutup  dan mulus  ( smooth ),  pada  bagian  demi  bagian.  Pada 1   displacement i u ditentukan,  dan  pada
2   vektor  stress j ij i n P   diketahui,  dimana
2 1         dan ) , , ( 3 2 1 n n n n  melambangkan  vektor  normal  satuan  mengarah  keluar  pada  permukaan   .  Dikehendaki untuk  menentukan  displacement  dan  stress  dalam  material.  Dengan  demikian  solusi  untuk persamaan  (3)  yang  berlaku  dalam  domain  dan  memenuhi  syarat  batas  pada   akan dicari.

## 4. Penurunan Persamaan dengan Koefisien Konstan
Dalam  pasal  ini  prosedur  yang  dikembangkan  dalam  Manolis  &  Shaw  (1996) dipakai untuk menurunkan suatu metode persamaan integral batas untuk kelas kofisien ) ( x  dan ) ( x  tertentu.  Penurunan ini  diperoleh  dengan  memperkenalkan  suatu transformasi  dari peubah  tak-bebas ) ( x u i untuk  mentransformasikan  persamaan  (3)  ke  suatu  persamaan dengan koefisien konstan. Kofisien ) ( x  dan ) ( x  disyaratkan mengikuti bentuk
, ) ( ) ( , ) ( ) ( ) 0 ( ) 0 ( x g x x g x   
   (4)
dimana ) 0 (  dan ) 0 (  konstan. Misalkan
) ( ) ( ) ( 2 / 1 x u x g x i i   (5)
Subtitusi persamaan (4) dan (5)  ke dalam persamaan (3), menghasilkan
(6)         0 2 1 ) ( , , , , 2 / 1 ) 0 ( ) 0 ( 2 / 1 , ) 0 ( 2 / 1 , ) 0 ( 2 / 1 , ) 0 ( , , ) 0 ( , ) 0 ( 2 / 1                     k k i i k k ij j ji i ki k j i j j i k k ij
g g g g g g g        
Mohammad Ivan Azis
Jika ) ( x g diasumsikan memiliki bentuk
  2 3 2 1 ) (         x x x x g (7)
dimana    , , dan  adalah konstanta, dan
) 0 ( ) 0 (    (8) sehingga   ) ( ) ( 2 3 2 1 ) 0 ( x x x x x           
 ,  maka  persamaan  (6)  dapat  ditulis sebagai     0 , , ) 0 ( , ) 0 (       j i j j i k k ij   
(9)
dengan ) 0 ( ) 0 (    .  Dengan  kata  lain,    jika i  adalah  solusi  displacement  dari  persamaan kesetimbangan  untuk  material  isotropik  homogen  dengan  konstanta  Lamé ) 0 (  dan ) 0 (  , maka  solusi  bersesuaian  dari  persamaan  kesetimbangan  untuk  material  elastis  isotropik  tak- homogen  dengan  parameter  Lamé ) ( x  dan ) ( x  yang  diberikan  dalam  bentuk multiparameter pada persamaan (4) dapat ditulis dari (5), dalam bentuk
  ) ( ) ( ) ( ) ( 1 3 2 1 2 / 1 x x x x x x g x u i i i             
Tensor stress yang diperoleh dari  persamaan (2), dituliskan dalam bentuk
    
    ik g ijk k ij g    2 / 1
dimana
  ) ( , , ) 0 ( , ) 0 ( i j k k k k ij ik           
dan vektor stress
         i g ik k i P g P P 2 / 1 (10)
dimana
        n P n P ij i j g ijk g ik       , (11)
Persamaan  integral  batas  untuk  solusi  persamaan  (9)  dengan i  diberikan  pada 1   dan
] [  i P diberikan pada 2   (dalam Brebbia & Dominguez , 1989), dalam bentuk
    ds P x i ij i ij j          
 ) ( 0 (12)
dimana 0 x adalah titik sumber, 0   jika 1 , 1 0      x jika   0 x , dan 2 1   jika    0 x , dan   memiliki tangen membelok secara kontinyu pada 0 x .  Juga untuk kasus tiga dimensi
    , ) ( ) ( ) ( , 2 / 1 , 2 / 1 ) 0 ( , 2 / 1 ) 0 ( i kj j ki k ij g ijk g g g         
Mohammad Ivan Azis
    , , , 4 3 ) 1 ( 16 1 ) 0 ( j i ij ij d d v d v     
 
(13)                  i j j i j i ij ij d n d n v d d v n d d v , , )( 2 1 ( , , 3 ) 2 1 ( ) 1 ( 8 1 2  
(14)
Sedangkan untuk kasus dua dimensi,
           j i ij ij d d d v d v , , 1 log 4 3 ) 1 ( 8 1 ) 0 (   (15)                  i j j i j i ij ij d n d n v d d v n d d v , , )( 2 1 ( , , 2 ) 2 1 ( ) 1 ( 4 1 2  
(16) dimana   ) 0 ( ) 0 ( ) 0 ( 0 ( 2 / ,        v x x d dan k k d n d  , /    .  Subtitusi  persamaan
(5) dan (10) ke persamaan (12), menghasilkan
    ds u P g P g x u x g k ij g ki kj i ij j          ) ( ) ( ) ( ) ( 2 / 1 2 / 1 0 0 2 / 1    (17)
Persamaan integral batas ini dapat dipakai untuk penentuan solusi i u dan ij  pada semua titik dalam domain  .

## 5. Metode Perturbasi
Prosedur  elemen  batas  yang  disebutkan  dalam  pasal  sebelumnya  menawarkan metode  numerik  efektif  untuk  penentuan  solusi ) ( x u i pada  saat ) ( x g memiliki  bentuk  (7) dan  parameter ) 0 (  dan ) 0 (  yang memenuhi persamaan (8). Pada pasal ini suatu prosedur baru ditemukan untuk kasus pada saat koefisien ) ( x  dan ) ( x  diganggu sedikit ( perturbed )
di sekitar ) ( ) 0 ( x  dan ) ( ) 0 ( x  , sementara persamaan (7) dan (8) tetap dipertahankan.
Koefisien ) ( x  dan ) ( x  disyaratkan mengambil bentuk
) ( ) ( ) ( ) 1 ( ) 0 ( x x g x     
(18) ) ( ) ( ) ( ) 1 ( ) 0 ( x x g x      (19)
dengan
) 0 ( ) 0 (    dan 0 2 / 1  j i g
dimana ) 1 (  dan ) 1 (  fungsi yang memiliki turunan kedua. Sehingga dari persamaan (3)
 
       j i j j i k k ij u u u g , , ) 0 ( , ) 0 (        . , , ) 1 ( , ) 1 ( j i j j i k k ij u u u        (20)
Penggunaan  transformasi  (5)  dan  mengikuti  analisis  yang  digunakan  untuk  memenuhi persamaan (6) dari (3) menghasilkan
Mohammad Ivan Azis
          i j i i k k ij , , ) 0 ( , ) 0 (        j i j i i k k ij cg , , ) 1 ( , ) 1 ( 2 / 1           (21)
Solusi dari persamaan (21) dicari dalam bentuk
) ( ) ( ) ( 0 x x r i r r i        (22)
Dari persamaan (5) dan (22)  displacement k u dapat dituliskan dalam bentuk deret
) ( ) ( ) ( 0 x u x u r k r r k      (23)
dimana ) ( r k u berkorespondensi dengan ) ( r k  menurut relasi
), ( ) ( ) ( 2 / 1 ) ( x u g x r k r k   untuk ,... 1 , 0  r (24)
Subtitusi persamaan (22) ke dalam (21), dan menyamakan koefisien dari suku-suku dalam  berpangkat menghasilkan
    ) ( ) ( , ) ( , ) 0 ( ) ( , ) 0 ( r j r i j r j i r k k ij h          (25)
dimana
0 ) ( ) 0 (  x h (26)     j r i j r j i r k k ij r u u u g x h ) 1 ( , ) 1 ( , ) 1 ( ) 1 ( , ) 0 ( 2 / 1 ) ( ) (            (27)
untuk ... , 2 , 1  r .
Persamaan integral untuk persamaan (25) adalah
  ) ( ) ( ) , ( ) ( ) , ( ) ( ) ( ) ( 0 ) ( 0 0 ) ( 0 x ds x P x x x x x x x r i ij r i ij r j                 ) ( ) , ( ) ( 0 ) ( x ds x x x h ij r i (28)
untuk ... , 1 , 0  r , dimana       j r i j r j i r k k ij i n P r ) ( , ) ( , ) 0 ( ) ( , ) 0 ( ) (           Juga     g ik r i r P r k P g i u P ) ( ) ( 2 / 1 ) (    , untuk ... , 1 , 0  r (29)
dimana
      j r i j r j i r k k ij i n u u u x P r ) ( , ) ( , ) 0 ( ) ( , ) 0 ( ) ( ) (        (30)
dan   g ik P diberikan  dalam  persamaan  (11).  Sehingga  persamaan  integral  (28)  dapat  ditulis dalam bentuk
      ) , ( ) ( )[ ( { ) ( ) ( ) ( 0 2 / 1 ) ( 0 ) ( 0 2 / 1 0 x x x g x u x u x g x ij r i r j    )]
, ( ) ( 0 x x x P k j g k i 
Mohammad Ivan Azis
   ) ( )]} , ( ) ( )[ ( 0 2 / 1 ) ( x ds x x x g x P ij r i    ) ( ) , ( ) ( { 0 ) ( x ds x x x h ij r i (31)
Nilai i P dapat ditulis sebagai
         1 ) ( ) ( ) 0 ( r r i r i r i i G gP r gP P  (31) dimana
    j r ji r ij r k k ij r i n u u u x G ) 1 ( ) 1 ( ) 1 ( ) 1 ( , ) 1 ( ) ( ) (          .
Untuk  memenuhi  syarat  batas  seperti  yang  disebutkan  dalam  pasal  sebelumnya  disyaratkan bahwa
i i u u  ) 0 ( pada 1   , i i
P g P 1 ) 0 (   pada 2   ,
dimana i u mengambil  nilai tetapannya  pada 1   dan i P mengambil  nilai  tetapannya  pada
2   . Akibatnya, bahwa untuk ... , 2 , 1  r dari persamaan (23) dan (31) diperoleh
0 ) (  r i u pada 1   , ) ( 1 ) ( r i r i G g P    pada 2  
Sekarang, persamaan integral (31) dapat digunakan untuk menentukan nilai numerik dari  anu  ( unknowns )  pada  batas   ,  dan  nilai  numerik  dari ) ( r i u dan  turunannya  dalam domain  untuk r = 0, 1, .... Persamaan (23) dan (31) kemudian memberikan nilai i u dan i P di seluruh domain  .

## 6. Hasil Numerik
Dalam  pasal  ini  beberapa  masalah  renggangan  bidang  ( plane  strain )  akan diselesaikan  secara  numerik  dengan  menggunakan  persamaan  integral  yang  diperoleh  pada pasal sebelumnya. Dalam pemakaian metode ini, untuk mendapatkan hasil numerik prosedur elemen batas baku, digunakan misalnya metode yang terdapat dalam Clements (1981). Untuk variasi  parameter  elastisitas  terpilih  ruas  kanan  dari  persamaan  (27)  dianggap  sangat  kecil sehingga hanya perlu untuk mempertahankan dua suku pertama dari ekspresi deret persamaan (23).
Untuk  semua  masalah  yang  diperhatikan  domain  berupa  suatu  bujur  sangkar bersisi l dan masing-masing sisi dari bujur sangkar itu dibagi menjadi sejumlah M (kelipatan dari 5) segmen dengan panjang sama. Integral numerik Gaussian dipakai untuk penghitungan integral  garis  pada  setiap  segmen.  Sedangkan  untuk  integral  domain  dalam  persamaan  (31), domain dibagi atas sejumlah 2 M sel bujur sangkar yang sama dan fungsi yang diintegralkan ( integrand )  diasumsikan  konstan  dan  mengambil  nilainya  pada  titik  pusat  dari  setiap  sel.
Tetapi,  nilai  dari  solusi  iterasi  sebelumnya ) 1 (  r i u beserta  turunannya  dalam  persamaan  (27) dihitung  hanya  pada  sejumlah  5  ×  5  titik  pusat  dari  bujur  sangkar  bagian  bersisi l /5,  dan diasumsikan  konstan  di  dalam  setiap  bujur sangkar  bagian  itu.  Dengan  kata  lain nilai ) 1 (  r i u beserta  turunannya  dalam  suatu  bujur  sangkar  bagian  tertentu  adalah  sama  untuk  sejumlah
2 5 / M sel yang termuat dalam bujur sangkar bagian itu.

## Mohammad Ivan Azis
Problem 1. Perhatikan masalah nilai batas yang memenuhi solusi analitik
    ' 1 ' 2 ' 2 ' 1 ' 1 ' 1 1 . 0 1 ( 4 . 4 / , 1 . 0 1 ( 4 . 4 / x x u x x u     dan ' 1 ' 22 ' 1 ' 11 32 . 0 1 , 12 . 0 1 x x       untuk suatu material dengan kofisien elastisitas
2 ' 1 ' ) 1 . 0 1 ( 2 . 1 ) ( x x    , (32)
2
' 1 ' ) 1 . 0 1 ( ) ( x x    (33) dimana      P u u u l x x ij ij i i i i / , / , / ' ' '
  (untuk 2 . 1 ,  j i ),
            , , / , / ' ' l dan  u masing-masing  adalah  panjang,  modulus  elastis  dan
displacement  rujukan,  serta l u P /     .  Koefisien  elastisitas  persamaan  (32)  dan  (33) memiliki
bentuk persamaan (18) dan (19) dengan
0 , ) 1 .
0 ( , , ) 1 . 0 1 ( ) ( ) 1 ( 2 ' 1 ) 1 ( ) 0 ( ) 0 ( 2 ' 1                x x x g dan 2 . 0   . Syarat batasnya (lihat Gambar 1) adalah
' 1 P dan ' 2 u diketahui pada AB,
' 1 P dan ' 2 u diketahui pada BC,
' 1 P dan ' 2 u diketahui pada CD,
' 1 P dan ' 2 u diketahui pada AD dimana   P P P i i / ' .
Tabel  1  dan  2  memperlihatkan  solusi  eksak  dan  MEB  untuk  beberapa  titik  dalam domain  dan  untuk  kasus  dimana  batas  domain  dibagi  atas  80,  160  dan  320  segmen.  Hasil dalam tabel konvergen ke solusi eksak sejalan dengan bertambahnya jumlah segmen. Hal ini sesuai dengan apa yang diharapkan, karena semakin kecil ukuran segmen semakin akurat nilai integral  numerik.  Khusus  untuk  kasus  320  segmen  solusi  displacement  memperlihatkan keakuratan sampai pada desimal keempat, dan solusi stress pada desimal ketiga.
Problem 2 . Perhatikan masalah nilai batas untuk suatu material bujur sangkar tak homogen yang diberi beban dari atas dan terklem pada bagian bawah serta dibiarkan bebas pada sisi kiri dan kanan, seperti yang diberikan pada gambar berikut.

## Mohammad Ivan Azis
Tabel 1. Solusi MEB dan eksak untuk ' ij  dari Problem 1.
Tabel 2. Solusi MEB dan eksak untuk ' i u dari Problem 1.
Pada Gambar 2 di bawah ini, terlihat material ini memiliki koefisien elastisitas
2 ' 2 ' 1 ' ) 1 ( 5 . 1 ) ( x x x  
    (34) 2 ' 2 ' 1 ' ) 1 ( ) ( x x x  
    (35)
Koefisien elastisitas pada persamaan (34) dan (35) memiliki bentuk persamaan (18) dan (19) dengan 2 ' 2 ' 1 ) 1 ( ) ( x x x g      ,
0 , ) 1 ( , ) 1 ( 2 ' 1 ' 1 ) 1 ( ) 0 ( ) 0 (                 x x dan 5 . 0   .  Bila 3 10 49 . 2     ksi,  maka  material  yang  diperhatikan  adalah magnesium alloy .    Kordinat  titik  sudut  material  ini  adalah    A(-0.5,-0.5),    B(0.5,-0.5),    C(0.5,0.5)    dan D(0.5,0.5), dan syarat pada batas adalah
0 , 0 ' 2 ' 1   u u pada AB, 0 , 0 ' 2 ' 1   u u pada BC, 0 , 0 ' 2 ' 1   u u pada CD, 0 , 0 ' 2 ' 1   u u
pada AD
Empat  kasus  yang  diperhatikan  menyangkut  kofisien  keelastikan  material  dan  yaitu kasus material homogen ( 0     ) dan tiga kasus material tak-homogen
) 1 . 0 , 1 . 0 , 0 , 1 . 0 , 1 . 0 , 0 (             .

## Mohammad Ivan Azis
Gambar 2. Geometri untuk Problem 2.
Gambar  3  dan  4  memperlihatkan  hasil  deformasi  sisi  material  dan  sisi  daerah 25 . 0 25 . 0 , 25 . 0 25 . 0 ' 2 ' 1       x x dalam  material.  Dari  kedua  gambar  ini  dapat disimpulkan bahwa fungsi ketakhomogenan g sangat mempengaruhi hasil deformasi. Dalam hal  ini  peningkatan  nilai  fungsi  ketakhomogenan g akan  meningkatkan  nilai  koefisien keelastikan khususnya rigiditas  . Sistem  kordinat baru ' i X dalam Gambar 3 dan 4 adalah sistem  untuk  material  hasil  deformasi.  Dalam  hal  ini  variabel  kordinatnya  didefinisikan sebagai ' ' ' i i i u x X   untuk 2 , 1  i .

## 7. Penutup
Metode  elemen  batas  untuk  solusi  suatu  kelas  masalah  nilai  batas  untuk  material isotropik tak-homogen telah berhasil ditemukan.
Gambar 3. Hasil Deformasi Sisi Material.

## Mohammad Ivan Azis
Gambar 4. Hasil Deformasi Sisi Daerah
, 25 . 0 25 . 0 ' 1    x 25 . 0 25 . 0 ' 2    x Dalam Material.
Metode  ini  cukup  mudah  digunakan  untuk  memperoleh  solusi  numerik  dari  suatu masalah  tertentu.  Hasil  numerik  yang  telah  diperoleh  mengindikasikan  bahwa  metode  ini mampu memberikan solusi yang akurat.

## Daftar Pustaka
[1]   Brebbia,  C.A.  dan  Dominguez,  J.,  (1989),  “Boundary  Elements  An  Introductory Course”, Computational Mechanics Publications , Boston.
[2]   Clements, D.L., (1981),”Boundary Value Problems Governed By Second Order Elliptic Systems”, Pitman.
[3].  Manolis, R.P. dan Shaw, R.P., (1996), “Green’s function for the vector wave    equation in a mildly heterogeneous continuum”, Wave Motion , Volume 24, pp : 59–83.
[4].   Rizzo,    F.J.,  (1967),  “An  integral  equation  approach  to  boundary  value  problems  of classical elastostatics”, Quarterly of Applied Mathematics ”, Volume 25, pp : 83–95.
[5].   Turcotte, D.L. &  Schubert, P., (1982), Geodynamic Applications of Continuum Physics to Geological Problems , John Wiley & Sons Inc., New York.