Jurnal Matematika dan Sains Vol. 7 No. 1, April 2002, hal 35 - 42
Pemakaian Metode Numerik Pada Sirkulasi Udara di Sekitar Bangunan Tradisional Bali I Gusti Bagus Wijaya Kusuma Program Studi Teknik Mesin, Fakultas Teknik, Universitas Udayana Kampus Bukit Jimbaran, Badung, 80361, Bali Phone: 62 – 361 – 703321, Fax. : 62 – 361 – 701806 Email:
[email protected] Diterima tanggal 5 Juli 2001, disetujui untuk dipublikasikan 22 Maret 2002
Abstrak Pola aliran udara pada sekumpulan bangunan tradisional Bali tidak dapat diestimasi dengan menggunakan Metode percobaan ataupun dengan memakai bangunan-bangunan modern, karena letak dan fungsi gedung pada sekumpulan bangunan tradisional Bali tersebut berhubungan erat satu sama lainnya. Adalah tidak tepat apabila kajian dilakukan dengan memakai satu bangunan yang terisolasi (isolated building) serta Metode simetri untuk mengetahui pola aliran udara di seluruh bangunan. Metode numerik dengan berdasarkan fi nite volume digunakan pada simulasi ini yakni untuk mempercepat proses analisa serta memberikan data yang lebih detail dari pola aliran udara di sekitar bangunan 3 -dimensi. Hasil menunjukkan bahwa Metode numerik memberikan hasil yang akurat tentang pola aliran udara di sekitar bangunan bangunan tradisional tersebut. Kata kunci : Metode numerik, sirkulasi udara, bangunan tradisional Abstract The airflow patterns around traditional Balinese buildings cannot be estimated by experiment method nor using modern buildings, since the position and the function of a building on a cluster of traditional Balinese buildings are linked each other. It is not correct if the investigation uses an isolated building and symmetrical method in order to understand the airflow patterns of all buildings. Numerical method based on finite volume is used in this simulation since it gives quicker analyses and capable of delivering more detailed and comprehensive information about the flow structure around 3-dimension buildings. Results show that numerical method gives more accurate prediction of airflow patterns around traditional buildings. Keywords : Numerical Method, air circulation, traditional building fungsinya berhubungan satu dengan lainnya. Metode numerik merupakan satu Metode alternatif dalam upaya menjawab permasalahan ini. Selain itu, Metode ini dapat mempercepat proses analisa, jauh lebih murah bila dibandingkan dengan Metode terowongan angin, serta dapat menyajikan informasi yang lebih lengkap dan komprehensip dari struktur aliran udara. Pola aliran udara di sekitar bangun-an senantiasa turbulen, sehingga kajian numerik harus memperhitungkan pola aliran turbulen tersebut. Pola aliran turbulen disajikan pa da Gambar 1. Metode pendekatan numerik untuk aliran turbulen yang paling banyak dipakai adalah memakai model standar k-ε. Namun kondisi – kondisi batas yang spesial haruslah diaplikasikan karena adanya sedikit kelebihan estimasi, yang berakibat pada rendahnya harga distribusi tekanan pada ujung depan bangunan apabila dibandingkan dengan metode terowongan angin 1,2), seperti
1. Pendahuluan Pola aliran udara di sekitar bangunan bangunan tradisional Ba li belum pernah dikaji oleh para peneliti secara menyeluruh. Padahal, interaksi antara aliran udara dan perpindahan panas di sekitar bangunan tradisional Bali memegang peran yang sangat penting dalam upaya meningkatkan kenyamanan termis penghuni, selain untuk melestarikan arsitektur tradisional Bali itu sendiri. Aliran udara pada bangunan tradi-sional Bali sangatlah kompleks serta saat ini masih sulit untuk dilakukan dengan mempergunakan simulasi terowongan angin ataupun dengan menggunakan model berupa bangunan modern. Aliran udara pada sekumpulan bangunan tradisional Bali tidak dapat diestimasi dengan menggunakan sebuah bangunan yang terisolasi dan mempergunakan kondisi-kondisi batas yang simetris, karena masing-masing bangunan baik posisi serta
35
36
JMS Vol. 7 No. 1, April 2002
disajikan pada Gambar 2. Untuk mengatasi hal tersebut penulis melakukannya dengan cara mereduksi harga k (energi kinetik) dan ε (energi disipasi), dimana berakibat pada berkurangnya arus eddy (eddy vortex) dan memperkecil produksi turbulensi pada bagian depan selubung bangunan tersebut.
utuh) merupakan satu kesulitan untuk dikaji secara eksperimen di laboratorium. section 4
section 4
section 3 section 2 Gedung 1 section 1 Gedung 4 section 4 section 3 section 2 section 1
Separation bubble
Reattachment point
section 3 section 2 section 1
section 4 section 3 section 2 section 1
Gedung 2
section 1 section 2 section 3 section 4
Reattachment length
Gambar 1. Pola sirkulasi udara pada sebuah gedung terisolasi Pressuredistribution C
0,5 A
Experiment(Castro& Robins,1977)
D
0
Cp
Windtunnel(Okadaet al,1991)
-0,5
Numerical(k-e,Qasim et al, 1992)
-1
Computational (k-e, Selvametal,1996)
-1,5
A
B
C
Gedung 3
section 1 section 2 section 3 section 4 Gedung 5
(a) Gambar 3. (a) Susunan gedung pada arsitektur tradisional Bali.
1 B
Gedung 6
D
Buildingsurfacelength
Beberapa asumsi perlu disampaikan, yaitu: 1. Wilayah yang dinyatakan sebagai atmospheric boundary layer (hingga ketinggian 500 m) telah dipenuhi dalam simulasi ini. Lapis batas permukaan pada ketinggian 10 m adalah lapisan yang mana harga momentum ke arah vertikal, panas dan fluks massa adalah konstan, dengan demikian teori similaritas dari Monin-Obukhov dapat diaplikasikan dalam kajian ini3). 2. Atmosfer diasumsikan dalam kondisi kesetimbangan hidrostatik, sehingga persamaan Boussinesq dapat dipergunakan dalam pendekatan4). 2.1 Kondisi batas
Gambar 2. Perbandingan distribusi tekanan pada satu bangunan terisolasi Harga k dan ε merupakan variabel bebas yang berfluktuasi. Semakin kecil harga yang dipakai, semakin kecil pula arus eddy pada bagian depan bangunan. Kedua parameter ini memiliki hubungan yang signifikan serta belum pernah dinyatakan oleh para peneliti.
Kondisi batas (boundary condition) untuk simulasi lapis batas atmosfer (atmospheric boundary layer) adalah untuk lapis batas homogen dengan karakteristik sebagai berikut: (a) Inlet Kondisi batas untuk kecepatan masuk adalah dinyatakan sebagai:
uy
2. Deskripsi Model Bentuk dari sekumpulan bangunan tradisional Bali dapat dilihat pada Gambar 3(a). Susunan bangunan yang unik (dimana letak gedung dan fungsinya adalah satu kesatuan yang
u 10m dimana
=
uy
ketinggian y,
log( y / y0 ) log( 10 / y0 )
(1)
adalah kecepatan tengah pada u10 m
adalah
kecepatan
pada
JMS Vol. 7 No. 1, April 2002
37
ketinggian 10, y 0 adalah kekasaran permukaan yang sebesar 0.010 m untuk rumput yang pendek. Selain itu, harga kecepatan baik ke arah y dan z adalah nol telah dimasukkan juga sebagai kondisi batas 5). Harga intensitas turbulensi adalah sebesar 6.2%, mengacu pada klasifikasi wilayah permukaan untuk wilayah pemukiman yang agak jarang penduduknya 6). Temperatur udara adalah sebesar 301 K mengacu pada hasil penelitian udara di Bali7). Persamaan momentum adalah: ∂u ∂P ∂ ∂ui ∂υt ∂u j uj i = − + (υ+υt ) + ∂xj ∂xi ∂x j ∂x j ∂xj ∂xi
Harga tegangan geser fluida adalah konstan dan dinyatakan sebagai ∂u (3) (µl + µt ) = τ w = ρ uτ2 ∂y Persamaan (3) di atas selanjutnya digunakan untuk mencari harga u τ , sedangkan para -meter µl adalah kekentalan laminar, µt adalah kekentalan turbulen, τw adalah tegangan geser fluida di permukaan selubung bangunan, ρ adalah massa jenis udara dan uτ adalah kecepatan gesek. Persamaan energi kinetik turbulensi, k dan laju energi dissipasi, ε pada kondisi masuk dinyatakan sebagai berikut: (4)
2
∂ µ t ∂ε ∂u ε ε2 − C2ρ = 0 (4) + C1µt ∂z σε ∂y k ∂y k Kekentalan Eddy:
υ t = Cµ
k2 ε
p
menyatakan harga tengah dari tekanan
ε=
uτ2 Cµ
uτ3 K ( y + y0 )
(8)
(9)
(10)
dimana K adalah konstanta von Karman’s (≅ 0.41) dan C µ = 0.09. (b) Bidang solid Semua kecepatan ke arah x, y dan z (u, v dan w) adalah nol pada bidang solid, serta kondisi ini disebut sebagai no-slip boundary condition. Temperatur pada permukaan tanah adalah konstan sebesar 305 K mengacu pada penelitian yang dilakukan di Bali7). Harga fluks panas telah dimasukkan dalam perhitungan dan besarnya bervariasi, tergantung pada posisi bangunan yang dikaji 7). Implementasi dari kondisi batas untuk aliran turbulen diawali dengan mengkaji terhadap harga y +, yakni lapis batas minimum yang diperlukan antara permukaan solid terhadap fluida, dimana pada kondisi dekat sekali dengan dinding memiliki harga y+ ≤ 11.63. Sehingga, apabila harga y+ > 11.63, maka node pertama (dari dinding solid) dinyatakan telah mengikuti kaidah logaritmik dari lapisan batas turbulen. Hubungan tersebut dinyatakan sebagai: uτ =
1 ln( Ey + ) K
(11)
dan (7)
u menyatakan komponen kecepatan pada sumbu x, y dan z., k menyatakan energi kinetik untuk kondisi turbulen ε merupakan energi disipasi untuk kondisi turbulen υt merupakan kekentalan eddy υ merupakan kekentalan udara ρ merupakan kerapatan massa dari udara __
uτ y + y0 ln K y 0
K=
(6)
dengan harga tekanan:
p 2 P= + k ρ 3
u=
(2)
dimana i = 1, 2, dan 3 menunjukkan arah dan j adalah indeks yang menunjukkan penjumlahan dari 1 hingga 3.
∂u ∂ µ t ∂k = + µt − ρε = 0 ∂z σ k ∂y ∂y
yang mana σk , σε , C 1 and C2 adalah konstanta dari model dengan σk , = 1.0, σ ε = 1.3, C1 = 1.44 and C2 = 1.923). Persamaan vektor di atas dapat diselesaikan dengan menggunakan persamaan
0. 75 9.24 σT , l − 1 × σT , t T + = σT ,t u + + 1 + 0 .28 exp − 0.007 σT , l σ T ,t
y+ =
ρu τ δ µ
(12)
(13)
dimana K adalah konstanta von Karman (0.41), E adalah sebuah konstanta integrasi (pada permukaan dinding yang halus serta memiliki tegangan geser yang konstan) bernilai sebesar
38
JMS Vol. 7 No. 1, April 2002
9,793. σ T, l adalah bilangan Prandtl untuk kondisi laminar (0.707) dan σ T,t adalah bilangan Prandtl untuk kondisi turbulen (≈ 0.9). (c) Outlet Kondisi lapis batas yang homogen dari Neumann dipergunakan dalam simulasi ini, dimana dui = dk = dε = 0 . dx j dx j dx j 2.2. Prosedur perhitungan numerik Semua persamaan diferensial di atas diselesaikan secara iterasi dengan Metode Picard, namun disertai perlakuan secara relaksasi dengan parameter relaksasi sebesar 0,7 setiap siklusnya. Parameter relaksasi bertujuan untuk mempercepat proses konvergensi. Hal ini dilakukan mengingat Metode Picard memberikan kestabilan dalam perhitungan namun memerlukan waktu yang lebih panjang untuk mencapai kriteria konvergen. Kriteria konvergen akan tercapai apabila jumlah dari kesalahan residual dari seluruh domain adalah lebih kecil dari harga toleransi yang diberikan (dalam simulasi ini adalah sebesar 10-6).
persamaan diferensial di atas adalah dengan metode finite volume. Produksi adalah hasil yang disebabkan karena adanya persamaan (4) dan (5) dimana akan berpengaruh terhadap tebal lapisan batas antara bidang solid terhadap udara (y + ). Semakin tinggi produksi yang dihasilkan maka semakin kecil tebal lapis batas yang dihasilkan, yang berakibat pada semakin tingginya tingkat penyimpangan yang dihasilkan. Besarnya produksi yang dihasilkan dengan menggunakan persamaan standar dari model k-ε adalah Pk = C µ εS 12 , dimana 1 S1 = 2
1
∂u i ∂u j + ∂x j ∂x i
2
2 k ε
(13)
Sedangkan produksi yang dihasilkan dengan menggunakan persamaan Kato dan Launder adalah Pk = C µ ε S 1 S 2
(14) 1
2 2 1 ∂ui ∂u j k + dimana S1 = dan 2 ∂x j ∂xi ε 1
2 2 1 ∂ui ∂u j k S2 = − 2 ∂xj ∂xi ε
a b c a=
first block grid (uniform grid) diselesaikan dengan Metode k-ε b = second bloc k grid (uniform grid) diselesaikan dengan Metode k-ε c = third block grid (non-uniform grid) diselesaikan dengan Metode k-ε Gambar 3b. Pengaturan grid dalam kajian ini Grid berbentuk segiempat diper-gunakan dalam simulasi ini, akan tetapi dengan bentuk yang mengikuti bangunan itu sendiri (body fittedgrid), seperti disajikan pada Gambar 3(b). Konfigurasi dari bentuk body fitted-grid ini serta dikombinasikan dengan harga k dan ε yang kecil akan memperbaiki keakuratan hasil serta waktu untuk proses iteras i. Metode penyelesaian
Dari persamaan (13) dan (14) di atas, jelaslah bahwa model k -ε menghasilkan produksi yang lebih besar daripada model yang disampaikan oleh Kato dan Launder, sehingga akan menghasilkan penyimpangan yang lebih tinggi terhadap hasil yang diharapkan. Pada wilayah sepa rasi (di belakang gedung ataupun di antara gedung), harga parameter regangan S1 dari model standar k -ε menjadi sangat besar pada wilayah stagnasi (kemandegan terhadap aliran udara setelah bertumbukan dengan permukaan bangunan), sehingga akan meningkatkan harga produksi Pk di atas level yang seharusnya (dibandingkan dengan Metode eksperimen tetapi untuk satu bangunan yang terisolasi), seperti dalam Gambar 2. Dengan mempergunakan harga ε yang kecil dalam simulasi perhitungan, maka produksi yang dihasilkan bila menggunakan model standar k-ε akan dapat diperbaiki (diperkecil) sehingga keakuratan hasil dapat ditingkatkan, baik bila dibandingkan terhadap Metode eksperimen (untuk satu bangunan yang terisolasi) maupun terhadap Kato dan Launder.
JMS Vol. 7 No. 1, April 2002
39
Inner region, diselesaikan dengan oneequations (LES)
Outer region , diselesaikan dengan twoequations
adalah selain mamp u memperbaiki produksi di sekitar permukaan bangunan juga mampu mempercepat proses konvergensinya. Selain itu, harga kesalahan relatif dari residu juga dapat diperkecil, jauh lebih kecil dari yang telah diperoleh oleh Zhou dan Stathopoulus2), seperti disajikan pada Gambar 4(a) dan 4(b).
Pressure distribution
1,5 1 Present study (bodyfitted grid)
0,5 0
Experiment(Castroand Robin, 1977)
-0,5
Computational (Selvam, 1996)
Cp
Penyelesaian persamaan dengan model standar k-ε juga disebut sebagai penyelesaian dengan two-equation method , karena menggunakan dua buah parameter (k dan ε ) untuk menyelesaikan persamaan aliran fluida. Zhou dan Stathopoulus 2) mencoba mempergunakan Metode yang berbeda, dima na mereka memakai one-equation method (Large Eddy Simulation, LES) untuk menyelesaikan persamaan di sekitar permukaan bangunan (solid wall) namun menggunakan two -equation method untuk menyelesaikan persamaan di luar lapisan batasnya (boundary layer), seperti disajikan dalam Gambar 3(c). Gambar lengkap dari simulasi disampaikan pada Gambar 3(d).
-1
Computational (Zhou and Stathopoulus, 1997)
-1,5 -2
Gambar 3c. Metod e two-layers
A
B
C
D
Buildingsurfacelength
(a)
Fully developed flow at inlet
Relative residues of the continuity equation
2.25 H
q”
q” q” q” q”
q”
q”
q” H q”
q”q”
Ts is constant
z
l2
l3
x 1 2 3 1 = kelompok rumah pertama 2 = kelompok rumah kedua 3 = kelompok rumah ketiga l 1 = jarak antara kelompok rumah pertama dan kedua l 2 = jarak antara kelompok rumah kedua dan ketiga H = tinggi rumah q” = fluks panas di permukaan T s = temperatur permukaan Gambar 3d. Bentuk geometri dan kondisi batas yang dipergunakan Akan tetapi, hasil yang dicapai oleh penulis dengan menggunakan perpaduan dari Gambar 3(b) dan harga k dan ε yang rendah
0,09 0,08
Relative residue
q”
0,07 0,06
Zhou
0,05 0,04 0,03
Present study
0,02 0,01 0,00
0
10
20
30
40 50 60
70
Number of iterations
(b) Gambar 4(a) Perbandingan distribusi tekanan pada bangunan terisolasi (b) Relative residues dari iterasi yang dipergunakan
40
3. Hasil dan Pembahasan Pola aliran udara di sekitar bangunan tradisional Bali disajikan pada Gambar 5(a) dan 5(b). Melihat struktur aliran di sekitar bangunan, maka penggunaan Metode numerik dengan model standar k-ε dapat menghasilkan struktur yang baik dan detail, khususnya pada daerah dimana separasi aliran terjadi. Hal ini terjadi karena model standar k -ε telah dimodifikasi sedemikian rupa, baik dengan jalan memberikan kondisi batas yang tepat, koreksi terhadap produksi yang akan dihasilkan, koreksi terhadap Metode iterasinya, koreksi terhadap bentuk grid yang dipergunakan untuk komputasi dan dengan menggunakan Metode penyelesaian finite volume untuk bentuk body-fitted grid tersebut. Dalam Metode penyelesaian finite volume untuk bentuk body-fitted grid, maka solusi algoritma adalah menyelesaikan persamaan momentum untuk memperoleh harga kecepatan. Penyelesaian persamaan momentum tersebut dilakukan dengan menggunakan persamaan tekanan (7), serta selanjutnya harga tekanan dikoreksi dengan menggunakan persamaan massa. Dengan demikian, maka akan terjadi saling koreksi dalam menyelesaikan persamaan momentum, tekanan dan massa untuk mendapatkan solusi yang sebenarnya. Karena menggunakan persamaan tekanan untuk mengoreksi penyelesaian persamaan momentum, maka penggunaan body-fitted grid adalah untuk menghindari pemakaian kondisi batas pada persamaan tekanan, dengan demikian akan terjadi interaksi yang tepat antara bentuk grid dengan Metode penyelesaiannya. Pola aliran udara di sekitar bangunan tradisional Bali tidak dapat dikaji dengan mempergunakan percobaan biasa, karena belum adanya wind tunnel yang khusus untuk model aliran udara di sekitar bangunan. Meskipun hasil Metode numerik lebih detail, namun Metode ini masih tetap memerlukan kajian kembali, khususnya terhadap pemilihan Metode dan kondisi- kondisi batas yang diterapkan. Pola aliran di sekitar bangunan tradisional Bali menunjukkan adanya hubungan yang erat antara letak gedung dan fungsinya, khususnya dalam upaya untuk meningkatkan kenyamanan termis para penghuninya. Streamline pada Gambar 5(a) menunjukkan bahwa aliran mengalami pemisahan (separasi) pada bagian depan gedung searah dengan aliran (windward side) serta menghasilkan arus putar pada bagian belakang gedung (leeward side). Arus putar (eddy vortex) di antara dua bangunan juga menunju kkan
JMS Vol. 7 No. 1, April 2002
adanya turbulensi yang tinggi, dimana hal ini tidak dapat dipantau dengan menggunakan percobaan biasa. Arus putar ini yang akan menumbuk bangunan sehingga memerlukan konstruksi yang kokoh untuk mengantisipasinya. Semakin tinggi arus putar, semakin kokoh konstruksi bangunan yang diperlukan. Arus putar juga berakibat semakin tingginya udara yang mengalir di selubung bangunan. Untuk meningkatkan kenyamanan penghuni dan mengurangi akibat yang dihasilkan oleh arus putar terhadap selubung bangunan, maka bagian – bagian selubung bangunan yang menerima arus putar tinggi tersebut harus dalam kondisi terbuka (tanpa dinding). Reattachment length
1
Reattachment point
4
2
6
2 3
5
3
5 (a)
Separasi
Daerah ini memiliki energi kinetik turbulensi yang tinggi, sehingga untuk meningkatkan kenyamanan termal penghuni dan mengurangi efek tumbukan pada gedung, maka selubung bangunannya harus terbuka (tanpa dinding)
Gambar 5(a). Plot streamline pada dasar gedung, (b). Plot energi disipasi pada dasar gedung
JMS Vol. 7 No. 1, April 2002
41
Sesaat setelah membentur bagian depan gedung, maka secara serentak aliran udara akan mencapai kemandegan (stagnasi) pada jarak sekitar 1,5 kali tinggi bangunan, yang ditandai dengan tercapainya harga nol untuk tegangan geser di permukaan (ground). Hal ini terlihat dari semakin jarangnya garis – garis/alur streamline itu sendiri. Turbulensi juga menjadi semakin berkurang pada titik stagnasi. Hal ini menunjukkan bahwa pada titik stagnasi terjadi proses penurunan kecepatan yang berakibat pada berkurangnya energi kinetik dan dissipasi (turbulensi) karena tegangan geser fluida di permukaan adalah nol. Titik ini disebut sebagai titik reattachment. Sedangkan jarak dari barisan gedung pertama hingga titik reattachment disebut reattachment length. Setelah titik reattachment ini tercapai, maka aliran berangsur-angsur meningkat kembali hingga mencapai kondisi aliran kembang penuh (fully developed flow ), yakni profil aliran yang sama seperti pada kondisi inlet. Kondisi fully developed flow , secara analitis baru akan tercapai pada saat jarak gedung diperpanjang hingga mencapai 1,5 kali reattachment length itu sendiri, atau sekitar 2,25 kali tinggi bangunan. Hal ini ditunjukkan oleh distribusi tekananan, dimana distribusi pada ruas gedung paling belakang hampir sama dengan yang terjadi pada ruas gedung paling depan. Selanjutnya, dalam arsitektur tradisional Bali letak pekarangan yang berada pada bagian tengah dari kumpulan gedung akan berakibat pada terdistribusinya udara secara merata ke semua gedung secara alami (natural ve ntilation), sehingga secara tidak langsung adalah bertujuan untuk meningkatkan kenyamanan termal dari penghuninya.
produksi yang akan dihasilkan, koreksi terhadap Metode iterasinya, koreksi terhadap bentuk grid yang dipergunakan untuk komputasi dan menggunakan Metode penyelesaian dengan finite volume untuk bentuk body-fitted grid tersebut. Meskipun metode numerik memberikan hasil yang lebih detail, namun keakuratannya masih perlu dikaji kembali dengan menggunakan Metode observasi yang lain, sebagai misal dengan menggunakan laser doppler anemometer.
4. Kesimpulan
6.
Metode numerik dengan model standar k-ε dapat dipergunakan serta menghasilkan struktur yang baik dan detail meskipun dipakai untuk mengkaji adanya separasi aliran. Hal ini terjadi karena model standar k-ε telah dimodifikasi sedemikian rupa, baik dengan jalan memberikan kondisi batas yang tepat, koreksi terhadap
Daftar Pustaka 1.
2.
3.
4.
5.
7.
Murakami S, Mochida, A, Ooka, R, Kato, S dan Iizuka, S, J. “Numerical prediction of flow around a building with various turbulence models: comparison of k-ε EVM, ASM, DSM, and LES with wind tunnel test”, dalam ASHRAE Transactions 102, 741-753 (1996). Zhou, YS dan Stathopoulus, T, J. “Application of two-layer methods for the evaluation of wind effects on a cubic building, dalam ASHRAE Transaction , 102 , 754-764 (1996). Monin, A.S. dan Obukhov, A.M. J. “Dimensionless characteristics of turbulence in the atmospheric surface layer”, dalam Doklady AN SSSR , 93, 223-226 (1953). Yoshida, A. J. “Two-dimensional numerical simulation of thermal structure of urban polluted atmosphere (effects of aerosol characteristics)”, dalam Atmospheric Environment, Part B , 25, 17-23 (1991). Hoxey, R.P. dan Richards, P. J. J. “Flow patterns and pressure field around a fullscale building”, dalam J . Wind Engineering and Industrial Aerodynamics, 50, 203-212 (1993). Wieringa, J, J. ”Updating the Davenport roughness classification” dalam J. Wind Engineering and Industrial Aerodynamics , 41-44, 357-368 (1992). Wijaya Kusuma, I G B, J. “Numerical Investiga tion of Airflow and Heat Transfer in Traditional Balinese uildings”, Ph.D. Thesis, Brunel University, 1999.
