Hey, what's going on?

Archive for the ‘nuclear engineering’ Category

Installing SRAC Code on Mandrivalinux 2010.0 x86_64

Posted by Syeilendra Pramuditya on January 18, 2016

SRAC: Standard Reactor Analysis Code, developed by Japan Atomic Energy Agency

SRAC code was developed long time ago, so there are problems when you try to install it on modern linux distros. The biggest problem of all is I think the original SRAC works only with F77/G77, whereas modern linux distros usually only provide gfortran.

Some installation tips:

Tip 1

Using the control panel, make sure that the following packages have been installed: gcc-cpp, gcc-c++, gcc-gfortran, tcsh, make, automake, ncompress, mscompress

Tip 2

Edit these 2 files:

SRAC/tool/install/SysDpnd/linux-g77/PrepInst.sh

SRAC/tool/install/SysDpnd/linux-g77/lmmake/lmmk/lmmk.sh

Open the files, you will see “set F77 = G77”, add a new line after that line: “set F77 = gfortran”

Now you just test run the code by executing SRAC/smpl/shr/Test.sh

Sample of input and output files are available here.

This has been tested and works for me.

*Of course you can just simply install SRAC 2006 which supports gfortran compiler.

Posted in nuclear engineering, software & simulation | Leave a Comment »

Catatan Perjalanan: Trieste, Italia (2)

Posted by Syeilendra Pramuditya on February 17, 2015

Alhandulillah.. hari Minggu siang tanggal 15 Feb saya sudah sampai di ICTP.. 🙂

Dua minggu lalu saya sebenarnya sakit, hari Selasa 3 Feb saya kena panas sampai 39 C. Tiga hari kemudian, Kamis sore 5 Feb, saya periksa darah dan ternyata trombositnya sudah turun ke 166rb.. astaghfirullah.. deg2an juga takut kena demam berdarah atau typhus. Malam itu saya langsung diopname di RS Borromeus..

Kalau demam berdarah, seminggu sejak mulai panas insya Allah bisa pulang dari RS, yang gawat kalau typhus, di RS bisa 2 minggu dan juga harus istirahat di rumah 2 minggu lagi.. wah hawatir juga.. bisa2 ga jadi ke Italia nih..

Setelah hasil lab keluar, ternyata saya kena demam berdarah dan bukan typhus, alhamdulillah.. lho sakit kok alhamdulillah??? 😀

Hari Senin 9 Feb saya sudah diperbolehkan pulang.. alhamdulillah..

Hari Sabtu 14 Feb saya berangkat ke bandara Soetta menggunakan bis Primajasa, berangkat dari rumah jam 7 pagi dan naik bis yang jam 8.30 pagi. Jam 11 saya sudah sampai di terminal 2 Soetta (hanya 2.5 jam!). Makan siang di Hokben, lalu muter2 di bandara sambil menunggu waktu check in. Sekitar jam 15.30 saya check in, dan jam 16.30 saya boarding. Pesawat berangkat hampir tepat waktu, yaitu sekitar jam 17.30. ICTP sudah mengatur tiket pesawat saya, dan mereka memberi saya maskapai Etihad Airways. Saya transit dua kali dalam perjalanan ini, yaitu di Abu Dhabi dan Roma.

Penerbangan dari Soetta ke Abu Dhabi memakan waktu 8 jam. Pesawatnya jumbo jet besar, deretan kursinya ada 10 kolom (3-4-3). Saya dapat seat number 18J. Sekitar pukul 22.30 waktu setempat (15, Jan, 1.30 WIB) saya alhamdulillah mendarat di Abu Dhabi. Di pesawat yang saya naiki ada banyak sekali jamaah umrah dari Indonesia, demikian juga di bandara Abu Dhabi, saya melihat banyak sekali jamaah umrah dengan seragam yang berbeda2. Bandaranya terlihat bagus dan sangat modern, juga terlihat sangat besar. Seragam polisinya terlihat berbeda, mereka hanya memakai baju khas Arab warna putih yang mirip daster ibu2 itu, kemudian memakai rompi warna candy green bertuliskan Airport Police, tentu lengkap dengan sorban di kepala.

Rombongan Umrah asal Indonesia di Abu Dhabi Airport

Rombongan Umrah asal Indonesia di Abu Dhabi Airport

Setelah menunggu sekitar 4 jam, saya melanjutkan penerbangan jam 3.00 waktu setempat (6.00 WIB). Berangkatnya agak delayed, jadwalnya sebenarnya jam 2.40. Pesawatnya tidak sebesar yang sebelumnya, kursinya ada 8 kolom (2-4-2) dan saya mendapat seat number 26D. Perjalanan memakan waktu sekitar 6.5 jam. Sekitar pukul 6.30 waktu setempat (15 Jan, 12.30 WIB) saya alhamdulillah mendarat di Roma. Bandaranya bagus, bahkan mungkin lebih bagus dari Abu Dhabi. Terutama toiletnya yang terlihat bagus, bersih, rapih, dan modern.

Transit 3 jam di Rome Airport

Transit 3 jam di Rome Airport

Setelah menunggu hampir 3 jam di Roma, saya melanjutkan penerbangan jam 9.20 waktu setempat (15.20 WIB). Pesawatnya lebih kecil lagi, kursinya hanya 6 kolom (3-3) dan saya mendapat seat number 10B. Perjalanan hanya memakan waktu sekitar 1 jam saja. Pukul 10.30 waktu setempat saya alhamdulillah mendarat di Trieste. Bandaranya tidak terlalu besar, tapi bersih, rapih, dan nyaman. Suhu udara di Trieste sedang lumayan dingin, sekitar 10 C.

Information Counter di Trieste Airport

Information Counter di Trieste Airport

Kemudian saya membeli tiket bis di information counter yang terletak di dalam bandara, harganya 3.30 euro. Bis nomor 51 dan berangkat dari depan bandara jam 12.05. Harus agak hati2 jangan sampai salah naik bis, karena bis 51 ini ternyata ada 2 rute, satu ke Trieste dan satu lagi ke Udine, dan waktu datangnya bersamaan juga. Naik bis yang berhenti di tempat parkir bertuliskan “Trieste” (di sebelahnya persis adalah tempat parkir bertuliskan “Udine”). Bis nya datang tepat waktu dan segera berangkat. Perjalanan memakan waktu sekitar 50 menit, dan sekitar jam 1 siang saya sampai di bus stop Grignano-Miramare (tepat sebelum terowongan). Ternyata ada satu orang lagi yang turun di situ, saya tanya ternyata dia juga menuju ICTP. Orang Vietnam, namanya Hai Hong Vo dari Ho Chi Minh University. Kami ternyata akan mengikuti workshop yang sama, dan juga tinggal di asrama yang sama, Adriatico Guest House (AGH). Dari bus stop, jalan menuju AGH agak tidak normal, bukan berupa jalan besar, tapi jalan kecil menurun mirip gang, terdiri dari tangga2 yang lumayan terjal, jadinya agak repot juga harus angkat2 koper gede saya.

Gang ga jelas menuju ICTP Adriatico Guest House

Gang ga jelas menuju ICTP Adriatico Guest House

Kalau dihitung2, saya berangkat dari rumah tanggal 14 Feb jam 7.00 WIB, dan sampai di Adriatico tanggal 15 Feb jam 13.00 (19.00 WIB), total perjalanan memakan waktu 36 jam! Kalau flight+transit nya memakan waktu sekitar 22 jam.. Walah.. pantes rasanya capek bgt pas sampe asrama..

ICTP Adriatico Guest House

ICTP Adriatico Guest House

Kamarnya tipe shared room, saya kebagian sekamar dengan anak Argentina, dia sampai di AGH satu jam sebelum saya, namanya Agustin Baceyro Ferra, katanya pegawai Komisi Nuklir Argentina. Kamarnya lumayan juga, twin bed, twin desk, dan kamar mandi dengan fasilitas lengkap, dan juga air keran siap minum. Window view nya? Wah jangan ditanya dah, mantappp abisss! Pemandangan langsung ke laut Adriatik! ICTP punya dua asrama, Adriatico dan Galileo, kalau Galileo tempatnya agak ke atas, jadi ga bisa liat laut. Alhamdulillah saya dapat di Adriatico.. 🙂

Window view dari kamar saya: laut Adriatik yang mantap abis

Window view dari kamar saya: laut Adriatik yang mantap abis

Walaupun capek bgt, setelah menaruh koper di kamar, saya sempatkan untuk jalan2 sebentar di sekitar AGH. Tepat di sebelah AGH ada taman yang luas dan indah, namanya Parco di Miramare, dan di dalamnya ada sebuah kastil yang juga indah, namanya Castello di Miramare. Lumayan ramai juga di situ, banyak orang2 lokal dan juga turis yang sedang jalan2.

Castello di Miramare

Castello di Miramare

Setelah itu saya kembali ke asrama, capek dan ngantuknya udah bener2 ga ketahan, akhirnya langsung tidur deh sore itu.. hehe..

(Bersambung..)

Posted in just a story, nuclear engineering | 2 Comments »

Catatan Perjalanan: Trieste, Italia (1)

Posted by Syeilendra Pramuditya on February 16, 2015

Kali ini saya insya Allah mengunjungi kota Trieste di Italia. Alhamdulillah.. 🙂

Saya terpilih sebagai salah satu peserta untuk mengikuti sebuah workshop saintifik di Centro Internazionale di Fisica Teorica Abdus Salam (CIFT), atau juga dikenal sebagai The Abdus Salam International Centre for Theoretical Physics (ICTP).

Nama kegiatannya adalah “Joint ICTP-IAEA Training Course on Physics and Technology of Water Cooled Reactors through the Use of PC-based Simulators“. Yup, jadi ini adalah kerjasama antara ICTP dan International Atomic Energy Agency (IAEA). Acara akan berlangsung selama 2 minggu full, dari tanggal 16 (Senin) sampai 27 (Jumat) Februari 2015.

The Invitation

The Invitation

Kalau bepergian keluar negeri, apalagi yang harus segera diurus kalau bukan paspor dan visa. Paspor saya kebetulan masih lama masa berlakunya, jadi saya hanya tinggal mengurus visa saja. Setelah saya menerima berbagai e-documents dari ICTP, saya langsung mengurus visa. Visa yang saya ajukan adalah visa Schengen, berjenis Study visa. Sebagai catatan, ada berbagai macam visa Schengen yang bisa diajukan di Kedubes Italia. Untuk visa selain study, visa harus diurus di sebuah agen visa yang bernama VFS Global, dan bukan langsung di Kedubes. Biaya pembuatan visa juga ternyata lumayan mahal, antara 1,2 – 1,5 juta. Nah hanya kalau kita apply Study visa kita bisa datang langsung ke Kedubes Italia. Berapa biayanya? Gratis total! 🙂

Oh iya Kedubes Italia beralamat di Jalan Diponegoro 45, Jakarta (6°12’2″S 106°50’23″E). Dekat dengan taman Suropati (~1 km). Saya tidak foto gedungnya, karena takut ditangkap polisi yang berjaga di luar.. hehe

Submission time adalah jam 11 – 12 siang, jadi saya naik kereta yang jam 6.35 dari Bandung, dan sampai di Gambir sekitar jam 10. Dari Gambir tinggal naik taksi dan bayar sekitar 25rb. Waktu tempuhnya kurang dari 30mnt. Ketika masuk Kedubes, kita harus meninggalkan KTP dan HP di pos keamanan. Yang antri sedikit sekali, klo ga salah ga sampe 5 orang aja, ga banyak orang kita yang main ke Italia kali ya. Beda banget dengan Kedubes Jepang yang selalu crowded. Dokumen2 yang diperlukan bisa dilihat disini atau disini.

Setelah kita submit dokumen2nya, kita akan diberi sebuah kartu untuk nanti ketika pengambilan visa dan paspor kita. Saya datang ke Kedubes Italia pada hari Selasa (20 Jan), dan diminta datang lagi hari Jumat (23 Jan) untuk pengambilan. Waktu pengambilan adalah jam 12 – 13 siang. Katanya mbak nya ga perlu nelfon dulu, langsung datang aja, optimis bgt mbak nya.

Kartu Pengambilan Visa

Kartu Pengambilan Visa

Oh iya, harus diperhatikan juga detail perjalanan kita. Dari info yang diberikan ICTP, saya akan transit di Abu Dhabi dan Roma. Lalu apa perlu visa transit UAE? Untuk transit sampai dengan 8 jam, tidak perlu visa transit UAE. Jadi saya tidak perlu mengurus visa transit.

Hari Jumat (23 Jan) saya datang lagi ke Kedubes Italia. Saya sampai di Gambir jam 10, dan sampai di Kedubes jam 10.30. Kecepetan sih.. kan katanya loket buka jam 12. Tapi ya sudah saya tunggu aja di dalam. Di dalam cuma ada 1 orang (mba2) yang menunggu, beneran jarang banget orang Indonesia yang ke Italia berarti ya. Mba2 ini juga sama2 apply visa studi, jadi katanya dia udah 3 taun sekolah di Florence, tapi karena sekarang pindah sekolah, harus apply visa lagi katanya. Nah tapi mba2 ini dapet panggilan interview dulu, aneh juga padahal dia udah lama di Italia. Saya yang belum pernah ke Italia malah ga dipanggil interview.

Jam 11 loket dibuka, daaan.. nomor saya yang dipanggil pertama! haha..

Alhamdulillah.. proses permohonan visa saya bisa selesai hanya dalam 3 hari, dan dikabulkan. 🙂

I got my very first Schengen Visa! Terimakasih Pak Federico 🙂

Posted in just a story, nuclear engineering | 3 Comments »

Pemecahan persamaan difusi neutron 1D single-group secara numerik dalam c++

Posted by Syeilendra Pramuditya on June 11, 2014

//Pemecahan persamaan difusi neutron 1D single-group secara numerik dalam c++
//Created: 11 Juni 2014
//Last Modified: -
//Programmer: Syeilendra Pramuditya

#include <iostream>
#include <cmath>
using namespace std;

float findflux(float a, float b, float c, float d, float e, float f){
float out1;
//Diturunkan dari persamaan difusi neutron 1D single-group
//Metode finite difference
out1=(((a+b)/(c*c))+(d/e)) / ((f/e)+(2/(c*c)));
return out1;
}

int main(){

int nx,itermax,i,iter,konvergen;
float L,dx,sigma_a,D,errormax,error,dummy;
float FluxOld[50],FluxNew[50],S[50];

//=============input===================
nx=50;//jumlah partisi
L=3;//dimensi domain
dx=L/nx;//lebar partisi
sigma_a=1;//cross section makroskopik
D=1.0/6;//koefisien difusi
errormax=1e-3;//kriteria konvergensi
itermax=1000;//loop limit
//=====================================

//==========syarat batas===============
FluxOld[0]=0; //flux di x=0 harus nol
FluxOld[nx]=0; //flux di x=L harus nol
FluxNew[0]=0; //flux di x=0 harus nol
FluxNew[nx]=0; //flux di x=L harus nol
S[0]=0; //source di x=0 harus nol
S[nx]=0; //source di x=L harus nol
//=====================================

//==========tebakan awal===============
for (i=1; i<=nx-1; i++){
FluxOld[i]=300; //tebakan awal flux
S[i]=100; //tebakan awal source
};
//=====================================

iter=0;
do{
iter++;
konvergen=1;
	for(i=1;i<=nx-1;i++){
	FluxNew[i]=findflux(FluxOld[i+1],FluxOld[i-1],dx,S[i],D,sigma_a); //Jacobi	
	}
	
	for(i=0;i<=nx;i++){
	error=fabs(FluxNew[i]-FluxOld[i]);
	if (error > errormax) konvergen=0; 
	}
	
	for(i=1;i<=nx-1;i++){FluxOld[i]=FluxNew[i];}
cout<<"Iterasi "<<iter<<"\n";
}while(konvergen==0 && iter<itermax);

cout<<"\n";
cout<<"Distribusi Flux\n";
for (i=0; i<=nx; i++){
cout<<FluxNew[i]<<"\n";
}

}

Posted in nuclear engineering, programming | Leave a Comment »

Revised Cheng-Todreas Correlations – Improvements in the Transition Flow Region

Posted by Syeilendra Pramuditya on April 8, 2014

Recently Cheng et al. published an article describing a revision to their correlations to improve the prediction of pressure drop coefficient in the transition flow region.

A computer code implementing the original correlations was developed and described in my previous post.

In the current post, I provide a modified code, in which the revision in the transition region has been incorporated.

 

Posted in nuclear engineering, software & simulation | Leave a Comment »

Basic Thermal-Hydraulics Experiments for Nuclear Engineering Students

Posted by Syeilendra Pramuditya on April 14, 2013

Teaching Assistantship Final Report (2013.03.12).

Posted in nuclear engineering, study & live in japan | Leave a Comment »

Water thermodynamic properties (in polynomials) at pressure of 155 bar (typical PWR operating condition)

Posted by Syeilendra Pramuditya on January 23, 2013

In the previous article, water thermodynamic properties at (nearly) atmospheric condition (p = 1 bar) were provided in form of easy-to-use polynomials based on XSteam libray.

In this article, water thermodynamic properties in form of polynomials at pressure of 155 bar, which is typical operating condition of PWR nuclear reactor, are provided. Please refer to XSteam homepage for more information.

Please note the range of applicability of the polynomials: 283.15 K < T < 613.15 K or equally 10 C < T < 340 C

Density of water as a function of temperature at pressure of 155 bar

 

 Specific heat capacity of water as a function of temperature at pressure of 155 bar

Specific heat capacity of water as a function of temperature at pressure of 155 bar

 

Thermal conductivity of water as a function of temperature at pressure of 155 bar

 

Dynamic viscosity of water as a function of temperature at pressure of 155 bar

Dynamic viscosity of water as a function of temperature at pressure of 155 bar

 

Saturation temperature of water as a function of pressure

 

Volumetric thermal expansion coefficient can be obtained from density equation by using this method.

The polynomials shown in this article are practically applicable for any pressure, since liquid water can be considered as an incompressible fluid. The relative differences of water properties between p = 1 bar and 160 bar are shown in the following table:

Property
Relative difference between p = 1 bar and 160 bar
Density
0.7%
Specific heat capacity
0.8%
Thermal conductivity
1.3%
Dynamic viscosity
0.5%

The weak dependency of the properties on pressure is shown by the following figures (calculated using XSteam):

Density of water as a function of pressure at temperature of 323.15 K (50 C)

Specific heat capacity of water as a function of pressure at temperature of 323.15 K (50 C)

Thermal conductivity of water as a function of pressure at temperature of 323.15 K (50 C)

Dynamic viscosity of water as a function of pressure at temperature of 323.15 K (50 C)


Ready-to-use Fortran functions

C     ******************************************************************
C     Calculation of water density as a function of temperature
C     WTRRHO -> kg/m3
C     T -> Kelvin
C     Available from http://wp.me/p61TQ-Tw
      REAL*8 FUNCTION WTRRHO(T)
      IMPLICIT NONE
      INTEGER I
      REAL*8 VAR1,C(0:6),T
       C(6)=0.0D0
       C(5)=0.0D0
       C(4)=0.0D0
       C(3)=-6.989D-6
       C(2)=6.312D-3
       C(1)=-2.422D0
       C(0)=1.351D3
       VAR1=0.0D0
        DO I=0,6
        VAR1=VAR1+C(I)*T**DFLOAT(I)
        END DO
      WTRRHO=VAR1
      RETURN
      END
C     ******************************************************************
C     ******************************************************************
C     Calculation of water specific heat capacity as a function of temperature
C     WTRCP -> kJ/kg.K
C     T -> Kelvin
C     Available from http://wp.me/p61TQ-Tw
      REAL*8 FUNCTION WTRCP(T)
      IMPLICIT NONE
      INTEGER I
      REAL*8 VAR1,C(0:6),T
       C(6)=0.0D0
       C(5)=0.0D0
       C(4)=2.276D-9
       C(3)=-3.755D-6
       C(2)=2.299D-3
       C(1)=-6.166D-1
       C(0)=6.515D1
       VAR1=0.0D0
        DO I=0,6
        VAR1=VAR1+C(I)*T**DFLOAT(I)
        END DO
      WTRCP=VAR1
      RETURN
      END
C     ******************************************************************
C     ******************************************************************
C     Calculation of water thermal conductivity as a function of temperature
C     WTRTC -> W/m.K
C     T -> Kelvin
C     Available from http://wp.me/p61TQ-Tw
      REAL*8 FUNCTION WTRTC(T)
      IMPLICIT NONE
      INTEGER I
      REAL*8 VAR1,C(0:6),T
       C(6)=0.0D0
       C(5)=0.0D0
       C(4)=0.0D0
       C(3)=0.0D0
       C(2)=-5.712D-6
       C(1)=4.757D-3
       C(0)=-2.946D-1
       VAR1=0.0D0
        DO I=0,6
        VAR1=VAR1+C(I)*T**DFLOAT(I)
        END DO
      WTRTC=VAR1
      RETURN
      END
C     ******************************************************************
C     ******************************************************************
C     Calculation of water dynamic viscosity as a function of temperature
C     WTRMU -> Pa.s
C     T -> Kelvin
C     Available from http://wp.me/p61TQ-Tw
      REAL*8 FUNCTION WTRMU(T)
      IMPLICIT NONE
      INTEGER I
      REAL*8 VAR1,C(0:6),T
       C(6)=0.0D0
       C(5)=0.0D0
       C(4)=6.746D-13
       C(3)=-1.32D-9
       C(2)=9.609D-7
       C(1)=-3.093D-4
       C(0)=3.738D-2
       VAR1=0.0D0
        DO I=0,6
        VAR1=VAR1+C(I)*T**DFLOAT(I)
        END DO
      WTRMU=VAR1
      RETURN
      END
C     ******************************************************************
C     ******************************************************************
C     Calculation of water saturation temperature as a function of pressure
C     WTRTST -> Kelvin
C     P -> Pa
C     Available from http://wp.me/p61TQ-Tw
      REAL*8 FUNCTION WTRTST(P)
      IMPLICIT NONE
      INTEGER I
      REAL*8 VAR1,C(0:6),P,PBAR
      PBAR=P*1.0D-5
       C(6)=-4.795D-10
       C(5)=2.594D-7
       C(4)=-5.512D-5
       C(3)=5.847D-3
       C(2)=-3.288D-1
       C(1)=1.053D1
       C(0)=3.705D2
       VAR1=0.0D0
        DO I=0,6
        VAR1=VAR1+C(I)*PBAR**DFLOAT(I)
        END DO
      WTRTST=VAR1
      RETURN
      END
C     ******************************************************************

Related links:

Posted in nuclear engineering | Tagged: , , , , , , , , | 3 Comments »

Subchannel Analysis Method at A Glance

Posted by Syeilendra Pramuditya on December 9, 2012

Posted in nuclear engineering | Tagged: , , | Leave a Comment »

Nuclear engineering schools around the world

Posted by Syeilendra Pramuditya on November 28, 2012

The University of California Berkeley campus, it has a nuclear engineering program.

For those who are looking for nuclear engineering schools around the world, here is the list of some of them:

Nuclear Science & Engineering OpenCourseWare (OCW)

Posted in nuclear engineering | Tagged: | 1 Comment »

ANSYS Workbench 13 Tutorial: 2D Meshing with prism layers (or boundary layers) on wall boundaries

Posted by Syeilendra Pramuditya on October 15, 2012

This post shows step-by-step how to create the following 2D mesh:

Posted in nuclear engineering, software & simulation | Tagged: , , , , , | Leave a Comment »

Laminar velocity profile in square flow channel

Posted by Syeilendra Pramuditya on October 13, 2012

In my previous article, I discussed about laminar velocity profile for 2D problem, the numerical solution technique is easily extendible to 3D problem. For example to solve isothermal laminar flow in a square channel problem, as illustrated below:

For this case, linear momentum equation takes the form:

The program is basically very similar to the 2D problem, just in the current case the calculated x-velocity is stored in 2D (or two-index) array. The Fortran code is shown below:

!     Laminar velocity profile in square flow channel
!     By Syeilendra Pramuditya - https://syeilendrapramuditya.wordpress.com
!     October 2012
      program main
      implicit none
      integer j,jmin,jmax,k,kmin,kmax,iter,itermax,conv,prnt
      parameter(jmin=0, jmax=40, kmin=0, kmax=40)
      real*8 mu,rho,Ly,Lz,uy1,uy2,uz1,uz2,errmax,Re,Dh,uinit,f,mdpdx,
     &alpha,uold(jmin:jmax,kmin:kmax),unew(jmin:jmax,kmin:kmax),dy,dz,
     &y(jmin:jmax),z(kmin:kmax),dy2,dz2,beta,err

      !input
      Ly=1.0d-2
      Lz=Ly !square channel
      uy1=0.0d0
      uy2=0.0d0
      uz1=0.0d0
      uz2=0.0d0
      itermax=5000
      errmax=1.0d-6
      Re=5.0d2

      !props
      mu=1.0d-3
      rho=1.0d3

      !calc internal vars
      Dh=4.0d0*(Ly*Lz)/(2.0d0*Ly+2.0d0*Lz)
      uinit=Re*mu/(rho*Dh)
      f=14.0d0/Re !friction factor
      mdpdx=(f/Dh)*0.5d0*rho*uinit*uinit ! -dp/dx
      alpha=mdpdx/mu

      !init all arrays
      do j=jmin,jmax
      do k=kmin,kmax
      uold(j,k)=0.0d0
      unew(j,k)=0.0d0
      end do
      end do

      !mesh
      dy=Ly/(jmax-jmin)
      dy2=dy*dy
      y(jmin)=0.0d0
      do j=jmin+1,jmax
      y(j)=y(j-1)+dy
      end do
      dz=Lz/(kmax-kmin)
      dz2=dz*dz
      z(kmin)=0.0d0
      do k=kmin+1,kmax
      z(k)=z(k-1)+dz
      end do

      beta=1.0d0/(2.0d0/dy2+2.0d0/dz2)

      !init guess of x-velocity
      do j=jmin+1,jmax-1
      do k=kmin+1,kmax-1
      uold(j,k)=uinit
      unew(j,k)=uinit
      end do
      end do

      !init boundary condition
      do k=kmin,kmax
      uold(0,k)=uy1
      uold(jmax,k)=uy2
      unew(0,k)=uy1
      unew(jmax,k)=uy2
      end do
      do j=jmin,jmax
      uold(j,0)=uz1
      uold(j,kmax)=uz2
      unew(j,0)=uz1
      unew(j,kmax)=uz2
      end do

      !calculate
      conv=0
      iter=0
      do while((conv.eq.0).and.(iter.le.itermax))
      conv=1
       !calc new u
       do j=jmin+1,jmax-1
       do k=kmin+1,kmax-1
       unew(j,k)=beta*(uold(j+1,k)+uold(j-1,k))/dy2 +
     &           beta*(uold(j,k+1)+uold(j,k-1))/dz2 + beta*alpha
       err=dabs(unew(j,k)-uold(j,k))/unew(j,k)
       if(err.ge.errmax) conv=0 !check convergence
       end do
       end do
       !update u
       do j=jmin,jmax
       do k=kmin,kmax
       uold(j,k)=unew(j,k)
       end do
       end do
      iter=iter+1
      end do

      !output result
      prnt=6
      write(prnt,*) 'Solution Converged in',iter,'iterations'
      write(prnt,*) ' '
      write(prnt,*) 'The result has been saved to the following files:'
      write(prnt,*) 'xveloc.dat -> calculated x velocity [m/s]'
      write(prnt,*) 'ypos.dat   -> y coordinate [m]'
      write(prnt,*) 'zpos.dat   -> z coordinate [m]'

      open(10,file="xveloc.dat")
      do j=jmin,jmax
      write(10,1000) (unew(j,k),k=kmin,kmax)
      end do
      close(10)
      open(10,file="ypos.dat")
      do j=jmin,jmax
      write(10,1000) y(j)
      end do
      close(10)
      open(10,file="zpos.dat")
      do k=kmin,kmax
      write(10,1000) z(k)
      end do
      close(10)

 1000 format(1P50E10.3)

      stop
      end

The calculated result can be plotted using the following Matlab m-file:

y=dlmread('ypos.dat',' ');
z=dlmread('zpos.dat',' ');
[Y,Z]=meshgrid(y,z);
VX = dlmread('xveloc.dat',' ');
surf(Y,Z,VX)
shading interp
title('x-Velocity [m/s]')
xlabel('y (cm)'),ylabel('z (cm)'),zlabel('x-Velocity [m/s]')
axis tight;
colorbar;
view(0,90);

And here is the resulting plot:

x-Velocity profile at Re = 1000

The real power of numerical solution methods is their ability to approximate the solution of a given mathematical equation for complex geometry. So let’s make this problem a bit more interesting by making the geometry a bit more complex, and we do this by “putting” an immersed rectangular solid object into the flow channel. This is achieved simply by applying a zero-velocity boundary condition at the object’s location. The modified code is available here, and the result is shown below:

x-Velocity profile with immersed solid rectangular object inside the flow channel

Posted in nuclear engineering, software & simulation | Tagged: , | Leave a Comment »

Problem Solving Series – Laminar Velocity Profile Between Infinite Parallel Plates

Posted by Syeilendra Pramuditya on October 12, 2012

The Document

The Code

!     Laminar velocity profile between infinite parallel plates
!     By Syeilendra Pramuditya - https://syeilendrapramuditya.wordpress.com
!     October 2012
      program main
      implicit none
      integer j,jmin,jmax,iter,itermax,conv,prnt
      parameter(jmin=0, jmax=10)
      real*8 uold(jmin:jmax),unew(jmin:jmax),L,uinit,dy,errmax,
     &mu,rho,Re,f,mdpdx,alpha,err,ubot,utop,y(jmin:jmax),
     &umath(jmin:jmax)

      !input
      L=1.0d-2
      ubot=0.0d0
      utop=0.0d0
      itermax=1000
      errmax=1.0d-6
      Re=5.0d2

      !props
      mu=1.0d-3
      rho=1.0d3

      !calc internal vars
      uinit=Re*mu/(rho*L) !average x-veloc
      f=24.0d0/Re !friction factor
      mdpdx=(f/L)*0.5d0*rho*uinit*uinit ! -dp/dx
      alpha=0.5d0*mdpdx/mu

      !init all arrays
      do j=jmin,jmax
      uold(j)=0.0d0
      unew(j)=0.0d0
      end do

      !mesh
      dy=L/(jmax-jmin)
      y(jmin)=0.0d0
      do j=jmin+1,jmax
      y(j)=y(j-1)+dy
      end do

      !init guess of x-velocity
      do j=jmin+1,jmax-1
      uold(j)=uinit
      end do

      !init boundary condition
      uold(jmin)=ubot
      uold(jmax)=utop
      unew(jmin)=ubot
      unew(jmax)=utop

      !calculate
      conv=0
      iter=0
      do while((conv.eq.0).and.(iter.le.itermax))
      conv=1
       !calc new u
       do j=jmin+1,jmax-1
       unew(j)=0.5d0*(uold(j+1)+uold(j-1))+alpha*dy*dy
       err=dabs(unew(j)-uold(j))/unew(j)
       if(err.ge.errmax) conv=0 !check convergence
       end do
       !update u
       do j=jmin,jmax
       uold(j)=unew(j)
       end do
      iter=iter+1
      end do

      !exact solution for static walls
      do j=jmin,jmax
      umath(j)=0.5d0*(1.0d0/mu)*mdpdx*(L*y(j)-y(j)*y(j))
      end do

      !result
      prnt=6
      write(prnt,*) 'Re =',Re
      write(prnt,*) 'Solution Converged in',iter,'iterations'
      write(prnt,*) ' '
      write(prnt,*) '      y [m]      u[m/s]  umath[m/s]'
      do j=jmin,jmax
      write(prnt,1000) real(y(j)),real(unew(j)),real(umath(j))
      end do
      write(prnt,*) ' '
      write(prnt,*) 'Result has been saved to fort.10 file'
      prnt=10
      write(prnt,*) 'Re =',Re
      write(prnt,*) 'Solution Converged in',iter,'iterations'
      write(prnt,*) ' '
      write(prnt,*) '      y [m]      u[m/s]  umath[m/s]'
      do j=jmin,jmax
      write(prnt,1000) real(y(j)),real(unew(j)),real(umath(j))
      end do

 1000 format(1P20E12.3)

      stop
      end

Posted in nuclear engineering, software & simulation | Tagged: , , | Leave a Comment »

Basic Thermal-Hydraulics Experiments for Nuclear Engineering Students

Posted by Syeilendra Pramuditya on September 28, 2012

During last spring semester, I worked as a teaching assistant (TA) for my department, my job was to teach new Japanese master students about some thermal-hydraulics experiments. I thought them the procedures of doing experiments, how to run some apparatuses, how to use some data acquisition equipments , some practical techniques, and most importantly safety aspect related to each of the experiments.

The activity was conducted once a week, every Wednesday morning from 9.00 to 12.15 noon. So every Tuesday I was also responsible to ensure that all the equipments, tools, and apparatuses were ready for the activity the next day.

There was also a Japanese senior technician (Matsui-san) working with me to prepare some materials such as heater pin samples and orifice for flow measurement.

The monthly salary was not so bad, about 10% of that of an assistant professor here :). Since we the TAs received the money, the department asked us to report about what have we done and how was the activity going. It was just a 10-min presentation plus 10-min Q&A, attended by 7 nuclear engineering professors (Profs. Iio, Igashira, Obara, Yano, Oguri, Matsumoto, Matsuda). The professors asked some questions to me, not so difficult though. Alhamdulillah everything went well. So, below is my presentation material.

Posted in nuclear engineering, study & live in japan | Leave a Comment »

Liquid Sodium Thermodynamic and Transport Properties

Posted by Syeilendra Pramuditya on February 8, 2012

The program can be accessed from this link.

Source code can be downloaded from this link. (or from here)

Related link:

Posted in nuclear engineering, software & simulation | Tagged: , , , , , , , | Leave a Comment »

Water Thermodynamic Properties

Posted by Syeilendra Pramuditya on August 20, 2011

Update+Fortran code Click here

Please cite as:
Syeilendra Pramuditya, Water Thermodynamic Properties, ITB Physics Department – Technical Document, http://portal.fi.itb.ac.id/tecdoc/waterprop

This post has been cited in these scientific articles:

  1. Numerical investigation of heat transfer in extended surface microchannels, published in International Journal of Heat and Mass Transfer (2016)
  2. Dual-pulse nonlinear photoacoustic technique: a practical investigation, published in Biomedical Optics Express (2015)
  3. Quantification of cooling channel heat transfer in low pressure die casting, Master Thesis at the University of British Columbia (2014)
  4. Stanford University Mechanical Engineering Design Proposal
  5. Investigating a New Approach for Harvesting Low-grade Thermal Energy Using an Electrochemical System, Research Poster, The Ohio State University
  6. DESARROLLO DE UN PRODUCTO INNOVADOR PARA LA PRODUCCIÓN DE HIELO EN EL REFRIGERADOR DOMÉSTICO, THESIS, UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO

There are many empirical and analytical correlations have been developed for thermodynamic properties of water. These correlations generally can be found from scientific papers as well as text books.

Correlations presented in this article are polynomial fit to data obtained from XSteam library for Matlab. XSteam itself is a kind of digital library of water properties based on the International Association for Properties of Water and Steam Industrial Formulation 1997 (IAPWS IF-97). These polynomial correlations are meant to be used in computer codes, and they are valid under the following conditions:

p=\textrm{1 bar}

5\;C \leq T \leq 95\;C

or equivalently:

278.15\;K \leq T \leq 368.15\;K


Water density as a function of temperature (p = 1 bar)

Polynomial fit:

\rho(T)={1001.1-0.0867T-0.0035T^2} (T in Celcius, unit is kg/m3)

\rho(T)={765.33+1.8142T-0.0035T^2} (T in Kelvin, unit is kg/m3)


Water specific heat capacity as a function of temperature (p = 1 bar)

Polynomial fit:

C_p(T)={4.214-2.286\times10^{-3}T+4.991\times10^{-5}T^2-4.519\times10^{-7}T^3+1.857\times10^{-9}T^4} (T in Celcius, unit is kJ/kg.C)

C_p(T)={28.07-0.2817T+1.25\times10^{-3}T^2-2.48\times10^{-6}T^3+1.857\times10^{-9}T^4} (T in Kelvin, unit is kJ/kg.K)


Water thermal conductivity as a function of temperature (p = 1 bar)

Polynomial fit:

k(T)={0.5636+1.946\times10^{-3}T-8.151\times10^{-6}T^2} (T in Celcius, unit is W/m.C)

k(T)={-0.5752+6.397\times10^{-3}T-8.151\times10^{-6}T^2} (T in Kelvin, unit is W/m.K)


Water dynamic viscosity as a function of temperature (p = 1 bar)

Polynomial fit:

\mu(T)={1.684\times10^{-3}-4.264\times10^{-5}T+5.062\times10^{-7}T^2-2.244\times10^{-9}T^3} (T in Celcius, unit is Pa.s)

\mu(T)={9.67\times10^{-2}-8.207\times10^{-4}T+2.344\times10^{-6}T^2-2.244\times10^{-9}T^3} (T in Kelvin, unit is Pa.s)


Water volumetric thermal expansion coefficient as a function of temperature (p = 1 bar)

The volumetric thermal expansion coefficient is evaluated from the density equation, as described in this article.

Polynomial fit:

\beta(T)=7.957\times10^{-5}+7.315\times10^{-6}T (T in Celcius, unit is /C)

\beta(T)=-1.908\times10^{-3}+7.318\times10^{-6}T (T in Kelvin, unit is /K)


Related links:

Posted in nuclear engineering | Tagged: , , , , , , | 10 Comments »