Нэг хэмжээст тархалтын тэгшитгэлийг Латтис Больцманы Арга ашиглан бодох
 |
| Молекул динамикаас макро орчны физик хэмжигдэхүүнрүү. Эх сурвалж: Transport Phenomena - Delft University of Technology |
Латтис больцманы аргаар нэг хэмжээст
тархалтын тэгшитгэлийг бодох код бичихэд өгөх тэмдэглэл.
Бодох тэгшитгэл:
Кодыг үе бүрчлэн тайлбарлан бичье.
program diffusionLBM
parameter (m=100) ! m is the number of lattice
nodes
real f1(0:1), f2(0:m), rho(0:m),
feq(0:m), x(0:m)
integer i
open(2,file='resultLBM.dat')
dt=1.0 ! time step
dx=1.0 ! space step
x(0)=0.0
Энэ хэсэгт хэрэглэгдэх хэмжигдэхүүн,
параметрүүдийг зарлаж зарим өгөгдлүүдийг өгч байна. f1(0:1) ба f2(0:m) нь
зангилаанаас гарч буй түгэлтийн функцыг илтгэнэ. Эдгээр нь нэгэнт m =100 ширхэг
зангилаатай учир 0 – ээс 100 хүртэл дугаарлагдах бүрдэл (array) тоо байна.
Өөрөөр хэлбэл f1 нь бодолтын алхам бүрд 100 ширхэг тоон цувааг үүсгэнэ гэсэн
үг. rho(0:m) хувьсагч
бөгөөд энэ бодлогийн хувьд температур нь юм. Гэхдээ температурыг масштаблан
хэмжээсгүй хувьсагч болгон бодож байгаа юм. feq(0:m) нь тэнцвэрт түгэлтийн функц бөгөөд мөн
л зуун ширхэг утга агуулах бүрдэл (array) юм. x(0:m) нь бодлогийн геометрийн уртыг илэрхийлэг
бөгөөд зүүнээс баруун тийш тоологдоно.
do i=1,m
x(i)=x(i-1)+dx
end do
csq=dx*dx/(dt*dt)
alpha=0.25 ! diffusive coefficient
omega=1.0/(alpha/(dt*csq)+0.5) ! eq.3.23
mstep=200 ! total number of time
steps
twall=1.0 ! left hand wall
temperature
Дээрх хэсгийн эхний гогцоо нь бодлогийн
геометрийг бүтээж байна. Дараа нь тайвшралтын давтамжийг тодорхойлохын тулд
температур тархалтын коэф-оор дамжуулан тэгшитгэл 3.23 ашиглан дараах байдлаар
тодорхойлно.
Бодлогийг 200с хүртэл хийнэ. Зүүн хананаас
өгөгдөх хэмжээсгүй температур нь twall=1.0 байна. Энэ нь гол хувьсагч болох rho(0:m) –ын анхны нөхцөл
юм.
do i=0,m
rho(i)=0.0 ! initial value of the
domain temperature
f1(i)=0.5*rho(i)
f2(i)=0.5*rho(i)
end do
Дээрх хэсэгт өгөгдсөн геомерт дотор анхны
нөхцлийг бий болгож бүх зангилаан дээрх температурыг 0 гэсэн хэмжээсгүй
температурт тохируулж байна. Мөн бүх түгэлтийн функцын нийлбэр нь хамаарах
хэмжигдэхүүнтэй тэнцүү байна гэдгээс түгэлтийн функыг дараах байдлаар
тодорхойлно.
Гэх мэтээр тодорхойлох ба энэ нь бүх
зангилааны турш бодогдоно.
! main loop
do kk=1,mstep
! collision process
do i=0,m
rho(i)=f1(i)+f2(i) ! dependent variable d.function
feq(i)=0.5*rho(i) ! equilibrium distribution function
! since k1=k2=0.5, then feq1=feq2=feq
f1(i)=(1.-omega)*f1(i)+omega*feq(i)
f2(i)=(1.-omega)*f2(i)+omega*feq(i)
enddo
! streaming process
do i=1,m-1
f1(m-i)=f1(m-i-1) ! f1 streaming
f2(i-1)=f2(i) !
f2 streaming
enddo
! boundary condition
f1(0)=twall-f2(0) ! constant temperature boundary condition,
x=0
f1(m)=f1(m-1) ! adiabatic boundary condition, x=L
f2(m)=f2(m-1) ! adiabatic boundary condition, x=L
enddo
Дээрх нь гол бодолтын гогцоог харуулах ба
хэсэг бүрт нь тайлбарлавал мөргөлдөөн, урсгал, хязгаарын нөхцлийн шинэчлэл
гэсэн гурван хэсгээс тогтоно. Эхний гогцоонд өгөгдсөн түгэлтийн функцаар гогцоо
бүрийн хамаарах хэмжигдэхүүнийг (rho(i)=f1(i)+f2(i)) тодорхойлно. Үүний дараа
тэнцвэрт түгэлтийн функцыг тодорхойлохдоо уг хамаарах хэмжигдэхүүнийг ашиглана.
Бодлогод D1Q2 схем хэрэглэж байгаа учраас k индекс нь 1 ба 2 гэсэн утга авах ба
тэнцвэрт түгэлтийн функц нь 2 ширхэг байна. Гэвч ижил үр дүн гарах учир
бодлогод зөвхөн feq1=feq2=feq хэмээн feq – ээр авсан болно. Жинлэх хүчин зүйл нь w1
+ w2 = 1 байна гэдгээр тус бүр нь 0.5тай (k1=k2=0.5) тэнцэх ёстой. Ингээд
мөргөлдөх алхамыг f1(i)=(1.-omega)*f1(i)+omega*feq(i) байдлаар зангилаа бүрд,
түгэлтийн функц бүрд бодно. Мөргөлдөх үеийн түгэлтийн функц нь өмнөх урсгалаас
ирсэн түгэлтийн функц, тэнцвэрт түгэлтийн функцаар тодорхойлогдоно. Харин
мөргөлдсөний дараа чиглэлийн дагуу дараачийн зангилааруу шилжих болно.
Энэ зурганд тасархай зураасаар f2(i) түгэлтийн
функцыг, тод зураасаар f1(i) функыг харуулсан болно. Мөргөлдөөний дараа
чиглэлийн дагуу f1(i) нь зүүнээс баруунлуу шилжих ба f1(1)=f1(0) тэнцүү, харин
f1(2)=f1(1) тэй тэнцүү гэх мэтээр эцэстээ f1(100)=f1(99) байх ёстой. Харин f2(i)
түгэлтийн функцын хувьд чиглэлийн дагуу баруунаас зүүнлүү урсах ба илэрхийлбэл f2(0)=f2(1),
харин f2(1)=f2(2), эцэстээ f2(99)=f2(100) болно. Хязгаарын нөхцл баруун ба зүүн
захад байх ёстой ба зүүн захын 0 гэсэн зангилаан дээр хэмжээсгүй температур 1
гэсэн утгыг өгөх ёстой. Иймд 0 цэгийн хувьд хоёр дахь тэгшитгэлээс улбаалан f1(0)=Twall-f2(0)
байна. Бодолтын алхам бүрд 0 цэгийн нэгдүгээр түгэлтийн функц нь энэ утгыг
дамжуулах болно. Энэ бол зүүн зах дээрх Диричлетийн хязгаарын нөхцөл юм.
Сүүлийн хоёр нь температурын графиент тэг байх адиабат хязгаарын нөхцөл юм. Энэ
мэтээр дээрх бодолтыг 200 хүртэл давтан хийнэ. Энэ үед бодогдсон урсгалаар
тодорхойлогдсон түгэлтийн функцууд хамаарах хэмжигдэхүүнийг тодорхойлж улмаар
мөргөлдөөний үеийн түгэлтийн функцыг тодорхойлох болно.
do i=0,m
write(2,*) x(i), rho(i)
enddo
Ингээд бодолт дуусаж үр дүнг хэвлэнэ. Энэ
аргаар бодсон үр дүнг Төгсгөлөг ялгаварын аргаар бодсон үр дүнтэй харьцуулж
харуулбал
Тооцоонд бараг ялгаа гарсангүй. Учир нь тархалтын тэгшитгэлийн төгсгөлөг
ялгаварын тэгшитгэл болон ЛБА-ын мөргөлдөөний тэгшитгэлүүд нь бараг адилтгал
тэгшитгэлүүд болж таардаг. Нюфлотоор хэрхэн дээрх графикийг байгаалах кодыг авч
үзье.
gnuplot> plot [0:30] 'resultFDM.dat' w points pt 7, 'resultLBM.dat' w
linespoints
gnuplot> set term png
Terminal type set to 'png'
Options are 'nocrop font "arial,12" fontscale 1.0 size 640,480 '
gnuplot> set output 'result.png'
gnuplot> replot
gnuplot> set term win
Terminal type set to 'windows'
Options are '0 color solid noenhanced'
Дээрх цэнхэрээр ялгасан нь терминалийн
өөрөө хэвлэх мэссэж юм. Ингээд ЛБА ба ТЯА-ын фортран кодуудыг хавсаргавал.
Латтис Больцманы код
1: program diffusionLBM
2: parameter (m=100) ! m is the number of lattice nodes
3: real fo(0:m), f1(0:1), f2(0:m), rho(0:m), feq(0:m), x(0:m)
4: integer i
5: open(2,file='resultLBM.dat')
6: dt=1.0
7: dx=1.0
8: x(0)=0.0
9: ! lattice arrangement
10: do i=1,m
11: x(i)=x(i-1)+dx
12: end do
13: csq=dx*dx/(dt*dt)
14: alpha=0.25 ! diffusive coefficient
15: omega=1.0/(alpha/(dt*csq)+0.5) ! eq.3.23
16: mstep=200 ! total number of time steps
17: twall=1.0 ! left hand wall temperature
18: ! initial condition
19: do i=0,m
20: rho(i)=0.0 ! initial value of the domain temperature
21: f1(i)=0.5*rho(i)
22: f2(i)=0.5*rho(i)
23: end do
24: do kk=1,mstep
25: ! main loop
26: ! collision process
27: do i=0,m
28: rho(i)=f1(i)+f2(i) ! dependent variable d.function
29: feq(i)=0.5*rho(i) ! equilibrium distribution function
30: ! since k1=k2=0.5, then feq1=feq2=feq
31: f1(i)=(1.-omega)*f1(i)+omega*feq(i)
32: f2(i)=(1.-omega)*f2(i)+omega*feq(i)
33: enddo
34: ! streaming process
35: do i=1,m-1
36: f1(m-i)=f1(m-i-1) ! f1 streaming
37: f2(i-1)=f2(i) ! f2 streaming
38: enddo
39: ! boundary condition
40: f1(0)=twall-f2(0) ! constant temperature boundary condition, x=0
41: f1(m)=f1(m-1) ! adiabatic boundary condition, x=L
42: f2(m)=f2(m-1) ! adiabatic boundary condition, x=L
43: enddo
44: ! end of loop
45: do i=0,m
46: write(2,*) x(i), rho(i)
47: enddo
48: stop
49: close(2)
50: end program
Төгсгөлөг ялгаварын аргын код
1: program diffusionfdm
2: parameter (m=100)
3: real dens,fo(0:m),f(0:m)
4: integer i
5: open(2,file='resultFDM.dat')
6: dx=1.0
7: dt=0.500
8: alpha=0.25
9: mstep=400
10: do i=0,m
11: fo(i)=0.0
12: end do
13: fo(0)=1.0 ! initial condition for old value of f at x=0.
14: f(0)=1.0 ! initial condition for updated value of f at x=0.
15: fo(m)=fo(m-1) ! initial condition for old value of f at x=L
16: f(m)=f(m-1) ! initial condition for updated value of f at x=L
17: do kk=1,mstep
18: ! main loop
19: do i=1,m-1
20: f(i)=fo(i)+dt*alpha*(fo(i+1)-2.*fo(i)+fo(i-1))/(dx*dx)
21: end do
22: do i=1,m-1
23: fo(i)=f(i) ! updating
24: end do
25: fo(m)=f(m-1) ! updating the boundary condition at x=L
26: end do
27: ! end of the main loop
28: x=0.0
29: do i=0,m
30: write(2,*) x,f(i)
31: x=x+dx
32: end do
33: stop
34: end
Төгсөв.
No comments:
Post a Comment