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 ∗)
Dosen tetap pada Jurusan Matematika Fakultas MIPA Universitas Hasanuddin Makassar
PENDAHULUAN Pendekatan persamaan integral untuk masalah nilai batas elastostatik sudah 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 (lihat misalnya Brebbia dan Dominguez (1989)). Dalam pertentangan, applikasi dari metode elemen batas ke masalah untuk material elastis isotropik tak-homogen sangat terbatas jumlahnya dan hal ini disebabkan oleh karena kesulitan dalam memperoleh 1
fungsi Green sebagai kernel dari persamaan integral batas yang bersesuaian. Baru-baru ini Manolis and 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´e λ dan µ sama. Hal ini berarti nilai ratio Poisson’s adalah 0, 25 yang membatasi applikasi. Tetapi seperti yang Manolis dan Shaw (1996) katakan nilai ratio Poisson’s ini adalah biasa untuk material batuan (lihat Turcotte dan Schubert (1982)). Tulisan ini mengembangkan kerja Manolis dan Shaw (1996) untuk membangun suatu prosedur perturbasi untuk menentukan solusi masalah statik dua dimensi untuk media isotropik tak-homogen dengan parameter Lam´e λ(x) = λ(0) g(x) + λ(1) (x), µ(x) = µ(0) g(x) + µ(1) (x), dimana x = (x1 , x2 , x3 ) adalah suatu vektor di R3 , g(x) suatu fungsi yang memenuhi batasan tertentu, λ(0) = µ(0) adalah konstanta dan suatu parameter kecil. Dalam batasan ini, bentuk parameter Lam´e di atas menawarkan pilihan keberagaman parameter λ(x) dan µ(x) yang cukup luas.
PERSAMAAN DASAR Dengan merujuk pada kerangka Cartesian Ox1 x2 x3 persamaan kesetimbangan (the equilibrium equations) dalam suatu material elastis tanpa gaya badan (body force) dapat ditulis dalam bentuk σij,j = 0,
(1)
dimana σij untuk i, j = 1, 2, 3 melambangkan tensor stress, tanda koma terindeks mengindikasikan turunan parsial terhadap kordinat 2
spasial xj dan konvensi penjumlahan indeks berulang (jumlahan dari 1 sampai 3) diterapkan. Hubungan displacement-stress dapat dituliskan sebagai σij = λ(x) δij uk,k + µ(x) (ui,j + uj,i ) , (2) dimana uk untuk k = 1, 2, 3 menyatakan displacement dan δij adalah operator delta Kronecker. Juga dalam Pers. (2) λ(x) dan µ(x) dengan x = (x1 , x2 , x3 ) adalah parameter Lam´e yang diasumsikan fungsi yang punya turunan kedua (twice differentiable). Substitusi Pers. (2) ke Pers. (1) menghasilkan [λ(x) δij uk,k + µ(x) (ui,j + uj,i )],j = 0.
(3)
MASALAH NILAI BATAS Suatu material isotropik tak-homogen berada pada daerah Ω dalam ruang tiga dimensi R3 dengan batas ∂Ω yang terdiri atas sejumlah berhingga permukaan tertutup dan mulus (smooth) bagian demi bagian. Pada ∂Ω1 displacement ui ditentukan dan pada ∂Ω2 vektor stress Pi = σij nj diketahui dimana ∂Ω = ∂Ω1 ∪ ∂Ω2 dan n = (n1 , n2 , n3 ) melambangkan vektor normal satuan mengarah keluar pada permukaan ∂Ω. Dikehendaki untuk menentukan displacement dan stress dalam material. Dengan demikian solusi untuk Pers. (3) yang berlaku dalam domain Ω dan memenuhi syarat batas pada ∂Ω akan dicari.
PENURUNAN PERSAMAAN DENGAN KOFISIEN KONSTAN Dalam pasal ini prosedur yang dikembangkan dalam Manolis dan Shaw (1996) dipakai untuk menurunkan suatu metode persamaan integral batas untuk kelas kofisien λ(x) and µ(x) tertentu. Penurunan ini diperoleh dengan memperkenalkan suatu transformasi dari peubah 3
tak-bebas ui (x) untuk mentransformasikan Pers. (3) ke suatu persamaan dengan kofisien konstan. Kofisien λ(x) dan µ(x) disyaratkan mengikuti bentuk λ(x) = λ(0) g(x) ,
µ(x) = µ(0) g(x),
(4)
dimana λ(0) dan µ(0) konstan. Misalkan ψi (x) = g 1/2 (x) ui (x)
(5)
Substitusi Pers. (4) dan Pers. (5) ke dalam Pers. (3) menghasilkan h
i
g 1/2 λ(0) δij ψk,k + µ(0) (ψi,j + ψj,i ) h
− λ
(0)
1/2 ψk g,ki
1/2 ψi g,jj
(0)
,j
(0)
1/2
i
+µ + µ ψj g,ij 1 g −1/2 [g,k ψk,i − g,i ψk,k ] = 0. − λ(0) − µ(0) 2
(6)
Jika g(x) diasumsikan memiliki bentuk g(x) = (αx1 + βx2 + γx3 + δ)2 ,
(7)
dimana α, β, γ dan δ adalah konstanta dan λ(0) = µ(0) ,
(8)
sehingga λ(x) = µ(0) (αx1 + βx2 + γx3 + δ)2 = µ(x) maka Pers. (6) dapat ditulis sebagai h
i
λ(0) δij ψk,k + µ(0) (ψi,j + ψj,i )
,j
=0
(dengan λ(0) = µ(0) ).
(9)
Dengan kata lain, jika ψi adalah solusi displacement dari persamaan kesetimbangan untuk material isotropik homogen dengan konstanta 4
Lam´e λ(0) dan µ(0) maka solusi bersesuaian dari persamaan kesetimbangan untuk material elastis isotropik tak-homogen dengan parameter Lam´e λ(x) dan µ(x) yang diberikan dalam bentuk multiparameter Pers. (4) dapat ditulis, dari Pers. (5), dalam bentuk ui (x) = g −1/2 (x) ψi (x) = (αx1 + βx2 + γx3 + δ)−1 ψi (x). Tensor stress diperoleh dari Pers. (2) diberikan oleh [g]
[ψ]
σij = −ψk σijk + g 1/2 σij dimana h
[g]
i
σijk = λ(0) δij (g 1/2 ),k + µ(0) δki (g 1/2 ),j + δkj (g 1/2 ),i , [ψ]
σij
= λ(0) δij ψk,k + µ(0) (ψi,j + ψj,i )
dan vektor stress [g]
[ψ]
Pi = −ψk Pik + g 1/2 Pi ,
(10)
dimana [g]
[ψ]
[g]
Pik = σijk nj ,
Pi
[ψ]
= σij nj .
(11)
Persamaan integral batas untuk solusi Pers. (9) dengan ψi diberikan [ψ] pada ∂Ω1 dan Pi diberikan pada ∂Ω2 diberikan dalam Brebbia dan Dominguez (1989) dalam bentuk η ψj (x0 ) =
Z ∂Ω
h
[ψ]
Φij Pi
i
− Γij ψi ds.
(12)
dimana x0 adalah titik sumber, η = 0 jika x0 ∈ / Ω ∪ ∂Ω, η = 1 jika 1 x0 ∈ Ω dan η = 2 jika x0 ∈ ∂Ω dan ∂Ω memiliki tangen membelok
5
secara kontinyu pada x0 . Juga untuk kasus tiga dimensi 1 [(3 − 4ν)δij + d,i d,j ] , (13) 16πµ(0) (1 − ν)d 1 = − 8π(1 − ν)d2 # " ∂d {(1 − 2ν)δij + 3d,i d,j } + (1 − 2ν) (ni d,j − nj d,i ) (14) ∂n
Φij = Γij
dan untuk kasus dua dimensi 1 1 = (15) (3 − 4ν) log δij + d,i d,j , (0) 8πµ (1 − ν) d 1 = − 4π(1 − ν)d " # ∂d {(1 − 2ν)δij + 2d,i d,j } + (1 − 2ν) (ni d,j − nj d,i ) , (16) ∂n
Φij Γij
dimana d = kx − x0 k, ν = λ(0) /(2(µ(0) + λ(0) )) dan ∂d/∂n = d,k nk . Substitusi Pers. (5) dan Pers. (10) dalam Pers. (12) menghasilkan ηg
1/2
(x0 )uj (x0 ) =
Z ∂Ω
h
[g]
i
g −1/2 Φij Pi − g 1/2 Γkj − Pki Φij uk ds.
(17) Persamaan integral batas ini dapat dipakai untuk penentuan solusi ui dan σij pada semua titik dalam domain Ω.
METODE PERTURBASI Prosedur elemen batas yang disebutkan dalam pasal sebelumnya menawarkan metode numerik efektif untuk penentuan solusi ui (x) pada saat g(x) memiliki bentuk Pers. (7) dan parameter λ(0) dan µ(0) memenuhi Pers. (8). Pada pasal ini suatu prosedur baru ditemukan untuk kasus pada saat kofisien λ(x) dan µ(x) diganggu sedikit 6
(perturbed) di sekitar λ(0) g(x) dan µ(0) g(x) sementara Pers. (7) dan Pers. (8) tetap dipertahankan. Kofisien λ(x) dan µ(x) disyaratkan mengambil bentuk λ(x) = λ(0) g(x) + λ(1) (x), µ(x) = µ(0) g(x) + µ(1) (x),
(18) (19)
dengan λ(0) = µ(0)
1/2
and g,ij = 0
dimana λ(1) dan µ(1) fungsi yang memiliki turunan kedua. Sehingga dari Pers. (3) h n
oi
g λ(0) δij uk,k + µ(0) (ui,j + uj,i )
=
,j
h
i
− λ(1) δij uk,k + µ(1) (ui,j + uj,i )
,j
.
(20)
Penggunaan transformasi (5) dan mengikuti analisis yang digunakan untuk menurunkan Pers. (6) dari Pers. (3) menghasilkan h
i
λ(0) δij ψk,k + µ(0) (ψi,j + ψj,i ) h
=
,j
i
−g −1/2 λ(1) δij uk,k + µ(1) (ui,j + uj,i )
,j
.
(21)
Solusi dari Pers. (21) dicari dalam bentuk ψi (x) =
∞ X
(r)
r ψi (x).
(22)
r=0
Dari Pers. (5) dan Pers. (22) displacement uk dapat dituliskan dalam bentuk deret ∞ uk (x) =
X r=0
7
(r)
r uk (x),
(23)
(r)
(r)
dimana uk berkorespondensi dengan ψk menurut relasi (r)
(r)
ψk (x) = g 1/2 uk (x),
untuk r = 0, 1, . . . .
(24)
Substitusi Pers. (22) ke dalam Pers. (21) dan menyamakan kofisien dari suku-suku dalam berpangkat menghasilkan h
(r)
(r)
(r)
λ(0) δij ψk,k + µ(0) ψi,j + ψj,i
i ,j
= h(r)
untuk r = 0, 1, . . . , (25)
dimana h(0) (x) = 0,
(26) h
(r−1)
h(r) (x) = −g −1/2 λ(1) δij uk,k
(r−1)
+ µ(1) (ui,j
(r−1)
+ uj,i
i
)
(27)
,j
untuk r = 1, 2, . . .. Persamaan integral untuk Pers. (25) adalah η(x0 )
(r) ψj (x0 )
=
Z ∂Ω
+
Z Ω
Γij (x, x0 )
(r) ψi (x)
− Φij (x, x0 )
[ψ (r) ] Pi (x)
(r)
hi (x) Φij (x, x0 ) dS(x)
ds(x) (28)
untuk r = 0, 1, . . ., dimana [ψ (r) ]
Pi
h
(r)
(r)
(r)
= λ(0) δij ψk,k + µ(0) ψi,j + ψj,i
i
nj .
Juga [ψ (r) ]
Pi
(r)
= g 1/2 Pi
(r)
[g]
+ uk Pik
untuk r = 0, 1, . . . ,
(29)
(30)
dimana (r)
h
(r)
(r)
(r)
Pi (x) = λ(0) δij uk,k + µ(0) ui,j + uj,i
8
i
nj ,
[g]
dan Pik diberikan dalam Pers. (11). Sehingga persamaan integral (28) dapat ditulis dalam bentuk (r)
η(x0 ) g 1/2 (x0 ) uj (x0 ) = [g] Pki (x)
Z
i
Z
Φkj (x, x0 ) −
n
h
(r)
ui (x) g 1/2 (x) Γij (x, x0 )−
∂Ω (r) Pi (x)
h
io
g 1/2 (x) Φij (x, x0 )
(r)
Ω
hi (x) Φij (x, x0 ) dS(x).
ds(x) + (31)
Nilai Pi dapat ditulis sebagai (0)
Pi = gPi
+
∞ X
(r)
(r)
r (gPi
+ Gi ),
(32)
r=1
dimana (r)
h
(r−1)
Gi (x) = λ(1) δij uk,k
(r−1)
+ µ(1) (ui,j
(r−1)
+ uj,i
i
) nj .
Untuk memenuhi syarat batas seperti yang disebutkan dalam pasal sebelumnya disyaratkan bahwa (0)
u i = ui (0) Pi = g −1 Pi
pada ∂Ω1 , pada ∂Ω2 ,
dimana ui mengambil nilai tetapannya pada ∂Ω1 dan Pi mengambil nilai tetapannya pada ∂Ω2 . Akibatnya, bahwa untuk r = 1, 2, . . . dari Pers. (23) dan Pers. (32) diperoleh (r)
ui = 0 (r) (r) Pi = −g −1 Gi
pada ∂Ω1 , pada ∂Ω2
Sekarang, persamaan integral (31) dapat digunakan untuk menentukan nilai numerik dari anu (unknowns) pada batas ∂Ω dan nilai 9
(r)
numerik dari ui dan turunannya dalam domain Ω untuk r = 0, 1, . . .. Pers. (23) dan Pers. (32) kemudian memberikan nilai ui dan Pi di seluruh domain Ω.
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 memperoleh hasil numerik prosedur elemen batas baku digunakan (lihat misalnya Clements (1981)). Untuk variasi parameter elastisitas terpilih ruas kanan dari Pers. (27) dianggap sangat kecil sehingga hanya perlu untuk mempertahankan dua suku pertama dari ekspresi deret Pers. (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 Pers. (31) domain dibagi atas sejumlah M 2 sel bujur sangkar yang sama dan fungsi yang diintegralkan (integrand) diasumsikan konstan dan mengambil nilainya pada titik pusat dari setiap sel. Tetapi, (r−1) nilai dari solusi iterasi sebelumnya ui beserta turunannya dalam Pers. (27) dihitung hanya pada sejumlah 5 × 5 titik pusat dari bujur sangkar bagian bersisi l/5 dan diasumsikan konstan di dalam setiap (r−1) bujur sangkar bagian itu. Dengan kata lain nilai ui beserta turunannya dalam suatu bujur sangkar bagian tertentu adalah sama untuk sejumlah (M/5)2 sel yang termuat dalam bujur sangkar bagian itu. Problem 1 Perhatikan masalah nilai batas yang memenuhi solusi 0 analitik u01 = x01 /[4.4(1 + 0.1x01 )], u02 = x02 /[4.4(1 + 0.1x01 )] dan σ11 = 10
0 0 1 + 0.12x01 , σ12 = −0.1x02 /4.4, σ22 = 1 + 0.32x01 untuk suatu material dengan kofisien elastisitas
λ0 (x) = 1.2(1 + 0.1x01 )2 , µ0 (x) = (1 + 0.1x01 )2 ,
(33) (34)
dimana x0i = xi /l, u0i = ui /u, σij0 = σij /P (untuk i, j = 1, 2), λ0 = λ/λ, µ0 = µ/λ, l, λ dan u masing-masing adalah panjang, modulus elastis dan displacement rujukan, serta P = λu/l. Kofisien elastisitas Pers. (33) dan Pers. (34) memiliki bentuk Pers. (18) dan Pers. (19) dengan g(x) = (1 + 0.1x01 )2 , λ(0) = µ(0) = λ, λ(1) = λ(1 + 0.1x01 )2 , µ(1) = 0 dan = 0.2. Syarat batasnya (lihat Gambar 1) adalah P10 dan u02 P10 dan P20 P10 dan u02 u01 dan u02
diketahui diketahui diketahui diketahui
pada pada pada pada
AB, BC, CD, AD,
dimana Pi0 = Pi /P . Tabel 1 – 2 memperlihatkan solusi eksak dan MEB untuk beberapa titik dalam domain dan untuk kasus dimana batas domain dibagi atas 80, 160 dan 320 segments. 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 11
x02 6
(0, 1)
P10 , u02 diketahui
D
C
u01 , u02 diketahui
P10 , P20 diketahui
A P10 , u02
B - x0 1 (1, 0)
diketahui Gambar 1: Geometri untuk Problem 1
Tabel 1: Solusi MEB dan eksak untuk σij0 dari Problem 1 Posisi (x01 , x02 ) (.1,.5) (.3,.5) (.5,.5) (.7,.5) (.9,.5) Posisi (x01 , x02 ) (.1,.5) (.3,.5) (.5,.5) (.7,.5) (.9,.5)
80 segmen 0 0 σ12 σ22 .9967 -.0124 1.0063 1.0027 -.0122 1.0223 1.0093 -.0124 1.0377 1.0165 -.0124 1.0528 1.0241 -.0125 1.0698 320 segmen 0 0 0 σ11 σ12 σ22 .9997 -.0130 1.0060 1.0052 -.0121 1.0208 1.0112 -.0122 1.0359 1.0174 -.0122 1.0511 1.0213 -.0134 1.0673 0 σ11
12
160 segmen 0 0 σ12 σ22 .9997 -.0127 1.0055 1.0043 -.0122 1.0214 1.0106 -.0124 1.0365 1.0173 -.0123 1.0516 1.0209 -.0109 1.0696 Eksak 0 0 0 σ11 σ12 σ22 1.0027 -.0114 1.0073 1.0082 -.0114 1.0218 1.0136 -.0114 1.0364 1.0191 -.0114 1.0509 1.0245 -.0114 1.0655 0 σ11
Tabel 2: Solusi MEB dan eksak untuk u0i dari Problem 1 Posisi (x01 , x02 ) (.1,.5) (.3,.5) (.5,.5) (.7,.5) (.9,.5)
80 segmen u01 u02 .0217 .1124 .0650 .1102 .1067 .1081 .1469 .1062 .1857 .1044
160 segmen u01 u02 .0220 .1124 .0655 .1102 .1073 .1081 .1476 .1062 .1864 .1044
320 segmen u01 u02 .0222 .1124 .0657 .1102 .1076 .1081 .1479 .1062 .1867 .1044
Eksak u02 .0225 .1125 .0662 .1103 .1082 .1082 .1487 .1062 .1877 .1043 u01
terlihat pada Gambar 2. Material ini memiliki kofisien elastisitas λ0 (x) = 1.5(1 + αx01 + βx02 )2 , µ0 (x) = (1 + αx01 + βx02 )2 .
(35) (36)
Kofisien elastisitas pada Pers. (35) dan Pers. (36) memiliki bentuk Pers. (18) dan Pers. (19) dengan g(x) = (1+αx01 +βx02 )2 , λ(0) = µ(0) = λ, λ(1) = λ(1 + αx01 + βx02 )2 , µ(1) = 0 dan = 0.5. Bila λ = 2.49 × 103 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 u01 = 0, P10 = 0, P10 = 0, P10 = 0,
u02 = 0, P20 = 0, P20 = −1, P20 = 0,
pada pada pada pada
AB, BC, CD, AD.
Empat kasus yang diperhatikan menyangkut kofisien keelastikan material λ dan µ yaitu kasus material homogen (α = β = 0) dan tiga kasus material tak-homogen (α = 0, β = 0.1; α = 0.1, β = 0; α = 0.1, β = 0.1). 13
D(-.5,.5) ? ? ? ? ? ? ? ? ? C(.5,.5) x02 6 - x0
1
A(-.5,-.5) @ @ @ @ @ @ @ @ @ B(.5,-.5) @ @ @ @ @ @ @ @ @ Gambar 2: Geometri untuk Problem 2
Gambar 3 dan 4 memperlihatkan hasil deformasi sisi material dan sisi daerah −0.25 ≤ x01 ≤ 0.25, −0.25 ≤ x02 ≤ 0.25 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 kofisien keelastikan khususnya rigiditas µ. Sistem kordinat baru Xi0 dalam Gambar 3 dan 4 adalah sistem untuk material hasil deformasi. Dalam hal ini variabel kordinatnya didefinisikan sebagai Xi0 = x0i + u0i untuk i = 1, 2.
SUMMARY Metode elemen batas untuk solusi suatu kelas masalah nilai batas untuk material isotropik tak-homogen telah berhasil ditemukan. Metode ini cukup mudah digunakan untuk memperoleh solusi nu14
0.2 0.1 0
X20
-0.1 -0.2 -0.3
α = 0.0, β α = 0.0, β α = 0.1, β α = 0.1, β
= 0.0 = 0.1 = 0.0 = 0.1
-0.2
0
-0.4 -0.5 -0.6
-0.4
0.2
0.4
0.6
X10 Gambar 3: Hasil deformasi sisi material
0 -0.05 -0.1
X20 -0.15 -0.2
α = 0.0, β α = 0.0, β α = 0.1, β α = 0.1, β
= 0.0 = 0.1 = 0.0 = 0.1
-0.1
0
-0.25 -0.3 -0.35 -0.3
-0.2
0.1
0.2
0.3
X10 Gambar 4: Hasil deformasi sisi daerah −0.25 ≤ x01 ≤ 0.25, −0.25 ≤ x02 ≤ 0.25 dalam material
15
merik dari suatu masalah tertentu. Hasil numerik yang telah diperoleh mengindikasikan bahwa metode ini mampu memberikan solusi yang akurat.
DAFTAR RUJUKAN Brebbia, C. A. dan Dominguez, J., (1989), Boundary Elements An Introductory Course, Computational Mechanics Publications, Boston. Clements, D. L., (1981), Boundary Value Problems Governed by Second Order Elliptic Systems, Pitman. 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, Halaman 59–83. Rizzo, F. J., (1967), An Integral Equation Approach to Boundary Value Problems of Classical Elastostatics. Quarterly of Applied Mathematics, Volume 25, Halaman 83–95. Turcotte, D. L. dan Schubert, P., (1982), Geodynamic Applications of Continuum Physics to Geological Problems, Wiley, New York.
16
