Изотерм бус үл шахагдах шингэний урсгалыг загварчлах
Томъёоллоороо бол изо гэдэг нь ижил гэсэн
утгыг илэрхийлэх ба изотерм гэдэг нь бодлогод авч үзэж буй материалд агуулагдах
дулааны тоо хэмжээнд өөрчлөлт орохгүй гэдэг нөхцлийг илтгэнэ. Өөрөөр хэлбэл бодлогийн
туршид температур тогтмол байх эсхүл температуртай холбоотой хэсгийг бодлогод
тооцохгүй үлдээнэ гэсэн үг юм. Харин изотерм бус гэдэг нь шингэний урсгалтай
хамт дулааны тоо хэмжээнд өөрчлөлт орох урсгалыг хэлнэ. Дулааны өөрчлөлт
урсгалын бусад физик хэмжигдэхүүнд нөлөөлнө. Үүний нэг жишээ болгон гэрлийн
чийдэн доторх аргонд дулаан хэрхэн дамжих процессыг харъя. Гэрлийн ламфан дотор
вольфрамын утас цахилгаан гүйдэлийн улмаар гэрэлтэж Жоулын нэгдүгээр хуулиар
дулаан ялгаруулах ба уг дулаан аргонд шингэж гэрлийн чийдэн тодорхи аргоны
урсгалыг бий болгоно.
Figure
1. Гэрлийн шилэн доторхи температурын орон
Дулаан тархах болон аргоны урсгах энэхүү
холимог үзэгдлийг байгалийн конвекци гэж нэрлэдэг.
Figure
2. Гэрлийн чийдэн доторх аргоны урсгал
Ерөнхийдөө хоёр төрлийн дулаан тархалтын
процесс байх ба үүний нэг нь дээрх жишээн д дурьдсан байгалийн конвекцийн
үзэгдэл юм. Дулаан тархалт урсгалаар тархах ба ямар нэгэн хүчний улмаар
тархалтын процесс идэвхижиж болно. Үүнийг хүчилсэн концвекцийн үзэгдэл гэж
нэрлэнэ. Хүчилсэн конвекцийн үзэгдлийн эх үүсвэр хүч нь ихэвчлэн гадны хүч байх
ба энэ төрлийн бодлогод Прандтлын тоо Prandtl number /Pr/ чухал үүрэг
гүйцэтгэнэ. Энэ тоо нь хэмжээсгүй тоо бөгөөд зунгааралтын коэффициент ба дулаан
тархалтын коэффицентын ноогдвороор анх Германы физикч Лудвик Прандтл
илэрхийлсэн учир түүний нэрээр нэрлэжээ.
Альфа нь дулаан тархалтын коэффициент ба
алфа=k/(ро cp) [м2/c] гэж илэрхийлэгдэнэ.
К нь дулаан дамжуулах чадвар, [В/м К]
гэсэн нэгжтэй байна. Ватт метр келвин.
Ср нь дулаан багтаамж [Ж/кг К] юм.
Жоуль килограмм келвин
Прандтлын тоо нэгээс их бол гадны хүчний
нөлөө их байх ба момент [хөдөлгөөний тоо хэмжээний] зонхилсон, нэгээс бага бол
дулаан тархалт зонхилсон буюу байгалийн конвекцит урсгал гэж үзэж болох юм. Энд
хүчилсан конвекцийн үзэгдлийн жишээг ЛБА-аар хэрхэн загварчлахыг D2Q9 сүлжээн
дээр авч үзье.
Халуун хоолой доторхи ламинар урсгал
Гаднаас үйлчлэх хүч нь температураас
хамаарах функц биш учраас моментийн тэгшитгэл болон скаляр орны тэгшитгэлүүд
тусдаа бодогдох болно. Скаляр оронд бидний авч үзэх температур эсвэл ууссан
бодис, масс тархалт, энергийн тэгшитгэл гэх мэт хамаарах болно. Харин байгалийн
конвекцийн үзэгдэлд урсгал ба дулаан тархалтыг идэвхижүүлэгч нь дулаанаас
хамааралтай функц байх учир дээрх тэгшитгэлүүд нэгэн дор тооцогдож болно. Иймд
хоёр өөр түгэлтийн функц тус бүр моментийн тэгшитгэлд f, температурын орны
тэгшитгэлд g гээд дараах байдлаар тооцогдоно.
Дээрхи мөргөлдөөний тэгшитгэлд байх тайвшралтын
давтамж нь момент ба скаляр оронд мөн өөр байдлаар тодорхойлогдох ёстой.
Үүнд алфа ба ню нь томъёо 1-д ашигласантай
ижил. Тэнцвэрт түгэлтийн функцыг адвекцийн бодлогын хувьд дараах байдлаар
тодорхойлно.
Хүйтэн tin=20oC
шингэн халуун tw=80oC хоолой [тэгш өнцөгт] дундуур анхны хурд uin=0.001м/с хурдтай урсана гэж бодлогын өгөгдлийг авъя. Хоолойн хана тогтмол
температураа бодлогийн туршид барих ба хоолойн өндөр 6.0см, урт нь 120см байна.
Бодлогын хувьд дундаж температур болох 50oC градустай үе дахь
кинематик зунгааралт ба дулаан тархалтын коэффициент тус бүр ню=6.0x10-7
м2/c, алфа=1.58x10-7 м2/c байна. Рейнольдсын
тоо ба Прандтлын тоог тус бүр тодорхойлбол:
Хэмжээсгүй температурыг дараах байдлаар
тодорхойлно.
ЛБА-д геометрийн нөлөөллийг тусгахын тулд
уртын масштабыг өндөрөөр авъя. Бодлогийн геометрийг тайлбарлавал,
Хэмжээсгүй температурын хувьд x=0 үед
төттө=0, у=0 ба y/H=1.0 үед төттө=1.0 байна.
хоолойн гаралтын хэсэгт x=L үед хэмжээсгүй температурын орон зайн
өөрчлөлт 0 байна гэж үзнэ. Ингээд фортран кодыг эхнээс нь тайлбарлая.
parameter (n=1000,m=50)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8), cx(0:8), cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
real g(0:8,0:n,0:m), geq(0:8,0:n,0:m), th(0:n,0:m)
integer i
Босоо чиглэлд 50 сүлжээ 6см-ийн зайнд,
хэвтээ чиглэлд 1000 сүлжээ 120см зайнд тус тус багтана. Масштабны коэффициент 0.12см
байна. Момент ба температур тархалт тус бүрдээ тооцогдох учраас тус бүрд
түгэлтийн функц f(0:8,0:n,0:m),
g(0:8,0:n,0:m) гэж, тэнцвэрт түгэлт мөн л feq(0:8,0:n,0:m), geq(0:8,0:n,0:m) гэж зарлагдана.
real w(0:8), cx(0:8), cy(0:8) гэдэг
мөрөнд эхнийх нь сүлжээний жинлэх хүчин зүйл, сүүлчийн хоёр нь босоо ба хэвтээ
сүлжээний хурднууд сүлжээний чиглэл дүрт харгалзах болно. real u(0:n,0:m), v(0:n,0:m)
нь макро орчин дахь хэвтээ ба босоо хурдны байгуулагчууд rho(0:n,0:m), th(0:n,0:m) нар
тус бүр нягт ба хэмжээсгүй температурууд болно. Зарлах хэсэг дууссан учир
тогтмол параметрүүдийг дараах байдлаар өгнө.
cx(:)=(/0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0/)
cy(:)=(/0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0/)
w(:)=(/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./)
uo=0.1
rhoo=5.00
dx=1.0
dy=dx
dt=1.0
tw=1.0
th=0.0
g=0.0
visco=0.05
pr=3.8
alpha=visco/pr
Re=uo*m/alpha
omega=1.0/(3.*visco+0.5)
omegat=1.0/(3.*alpha+0.5)
Эхний гурван мөрөөр сүлжээний хэвтээ, босоо хурд, жинлэх хүчин зүйлийг өгч
байна. Сүлжээний хурд сүлжээ бүрт 9 ширхэг байх ба байгуулагчийг өмнөх жишээний
/хөндийн урсгал/ зурагтай харьцуулж харж болно. Анхны хурдыг uo=0.1, анхны нягтыг rhoo=5.00,
орон зайн ба хугацааны алхамыг dx=dy=dt=1.0,
ханын температурыг tw=1.0, анхны шингэний
температурыг th=0.0 өгч Рейнольдсын тоог (6)-аар,
тайвшралтын хэлбэлзлүүдийг тус бүр (4)-өөр omega-wm, omegat-ws гэж тодорхойлов. Ингээд
шингэний анхны нөхцлийг дараах байдлаар өгнө.
do j=0,m
do i=0,n
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
end do
do j=1,m-1
u(0,j)=uo
v(0,j)=0.0
end do
Эхний хоёр гогцоонд шингэний нягт ба хурднуудыг бүх сүлжээний хувьд анхны
байдлаар өгнө. Үүний дараа оролтын хэсэгт хэвтээ хурд оролтын хурдтай тэнцүү
байх учир дахин j=1,m-1 үед u(0,j)=uo гэж өгнө. Ингэснээр урсгалыг эхлүүлэх
боломж бий болно. Температурын орон нэгэнт 0 учир түүнийг өгөхгүйгээр бодлогийн
гол гогцоог эхлүүлж болно. Анхны нөхцлийг өгч дууссан учир бодлогийн гол
гогцоонд шилжье.
DO kk=1,mstep
call
collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
call streaming(f,n,m)
call sfbound(f,n,m,uo)
call rhouv(f,rho,u,v,cx,cy,n,m)
do j=0,m
do i=0,n
sum=0.0
do k=0,8
sum=sum+g(k,i,j)
th(i,j)=sum
end do
end do
end do
! collision for
scalar
call
collt(u,v,g,geq,th,omegat,w,cx,cy,n,m)
! streaming for
scalar
call streaming(g,n,m)
call gbound(g,tw,w,n,m)
print *,
th(n/2,m/2),v(0,m/2),rho(0,m/2),u(n,m/2)
print*, v(n,m/2),rho(n,m/2)
END DO
Эхний ногоон хэсэгт өмнөх бодлоготой ижилээр моментийн тэгшитгэлийг
шийдвэрлэж нягт ба хурдны бүрдүүлэгчийг тооцох дэд програмууд байна. Харин
дараагийн цэнхэр хэсэгт шинээр нэмэгдэж орж буй температурын тархалтыг бодохын
тулд бүх сүлжээнд температурыг түгэлтийн функаар тодорхойлох гогцоог оруулж өгч
байна. Энд байх g(k,i,j) нь температурын
орны түгэлтийн функц ба параметрийн өгөгдлийн хэсэгт g=0.0
гэж өгөгдсөн байгаа. Харин th(i,j) нь
хэмжээсгүй температурыг илэрхийлэх төттө юм. Сүүлчийн ногоон хэсэгт
температурын хувьд мөргөлдөх процессийг дэд програмаар call collt, урсгалыг call
streaming(g,n,m) дэд програмаар гэхдээ моменттой ямарч ялгаагүй, тус тус
үнэлж хязгаарын нөхцлийг call gbound(g,tw,w,n,m) –аар шийдвэрлэнэ. Ингээд
сүүлийн ногоонд дурьдагдсан дэд кодуудыг тайлбарлая.
Дулаан тархалтын мөргөлдөх процесс call collt
Моментийн тэгшитгэлээр үнэлэгдсэн хурд дулааны тархалтанд хэрэглэгдэнэ.
Учир нь дулаан урсгалын хурдаар тархах ёстой.
subroutine
collt(u,v,g,geq,th,omegat,w,cx,cy,n,m)
do i=0,n
do j=0,m
do k=0,8
geq(k,i,j)=th(i,j)*w(k)*(1.0+3.0*(u(i,j)*cx(k)+v(i,j)*cy(k)))
g(k,i,j)=omegat*geq(k,i,j)+(1.0-omegat)*g(k,i,j)
end do
end do
end do
return
end
Гол үйлдэл цэнхэрээр ялгасан хэсэгт бий болох ба эхний мөр нь тэгшитгэл (5)-
ийн дагуу температурын тэнцвэрт түгэлтийн функцийг geq(k,i,j)=th(i,j)*w(k)*(1.0+3.0*(u(i,j)*cx(k)+v(i,j)*cy(k)))
гэж тодорхойлно. Үүнд дээрх өгүүлбэр ёсоор макро орчны хурднууд нь моментийн
бодолтоос олдсон байна. Хоёр дахь мөр нь мөргөлдөх алхамыг тэгшитгэл (3) – ийн
дагуу бодно.
Моментийн бодолтын хязгаарын нөхцөл call sfbound
Тэг эрэмбийн момент тархалт буюу нягтын бодолтын хувьд хязгаарын нөхцөл нь
өмнөх жишээ болох хөндийн урсгалтай ижил боловч оролт ба гаралтын нөхцөл үл
ялиг ялгаатай байна. Үүнийг тасдаж аван
доор талйбарлая.
! flow in on west boundary
rhow=(f(0,0,j)+f(2,0,j)+f(4,0,j)+2.*(f(3,0,j)+f(6,0,j)+f(7,0,j)))
rhow=rhow/(1.-uo)
f(1,0,j)=f(3,0,j)+2.*rhow*uo/3.
f(5,0,j)=f(7,0,j)+rhow*uo/6.
f(8,0,j)=f(6,0,j)+rhow*uo/6.
end do
! account for open boundary condition at the outlet
do j=1,m
f(1,n,j)=2.*f(1,n-1,j)-f(1,n-2,j)
f(5,n,j)=2.*f(5,n-1,j)-f(5,n-2,j)
f(8,n,j)=2.*f(8,n-1,j)-f(8,n-2,j)
end do
return
end
Баруун хэсэгт хоолойн урсгалын оролтын хэсэг байх учир оролтын хязгаарын
нөхцлийг өгнө. Энэ хэсэгт f1, f5, f8 нар урсгалын процессоос олдохгүй учир
тэдгээрийг f3, f7, f6 –аар дамжуулан тодорхойлох болно. Ингэхдээ хязгаарын
нөхцөлтэй нормал тэнцвэрт нөхцлийг ашиглана.
Дээрхи хоёр тэгшитгэлээс хөөн хязгаарын нөхцлийг дараах байдлаар тодорхойлж
болно. Тод хараах бичигдсэн макро орчны хурд, сүлжээний хурд нь вектор
хэлбэрээр байгаа учир тод хараах тэмдэглэсэн болно.
Дээрх тэгшитгэлд
индекс w нь зөвхөн west гэдэг утгыг илэрхийлж байгаа юм. Нягтыг эхний
тэгшитгэл ёсоор rhow=(f(0,0,j)+f(2,0,j)+f(4,0,j)+2.*(f(3,0,j)+f(6,0,j)+f(7,0,j)))/(1.-uo) гэж тодорхойлох ба код урт учир таслан
хоёр мөр болгон бичсэн болно. Ингээд f1, f5, f8 тус бүрийг томъёо ёсоор
тодорхойлно. Хоёр дахь хэсэгт гаралтын хэсэгт зориулан нээлттэй хязгаарын
нөхцлийг өгнө. Бид гаралтын хурд ямар байхыг мэдэхгүй, мөн сүлжээний f3, f6, f7
гэсэн түгэлтийн функцуудыг мэдэхгүй учир буцаж ойх хязгаарын нөхцлөөс ялгаатайгаар
түгэлтийн функцад экстраполяци хийж тодорхойлно. Үүнийг нээлттэй
хязгаарын нөхцөл гэж нэрлэнэ.
Дулаан тархалтын хязгаарын нөхцөл call gbound
! West boundary condition, the temperature is 0
do j=0,m
g(1,0,j)=-g(3,0,j)
g(5,0,j)=-g(7,0,j)
g(8,0,j)=-g(6,0,j)
end do
! Right hand
boundary condition, open
do j=0,m
g(6,n,j)=2.*g(6,n-1,j)-g(6,n-2,j)
g(3,n,j)=2.*g(3,n-1,j)-g(3,n-2,j)
g(7,n,j)=2.*g(7,n-1,j)-g(7,n-2,j)
g(2,n,j)=2.*g(2,n-1,j)-g(2,n-2,j)
g(0,n,j)=2.*g(0,n-1,j)-g(0,n-2,j)
g(1,n,j)=2.*g(1,n-1,j)-g(1,n-2,j)
g(4,n,j)=2.*g(4,n-1,j)-g(4,n-2,j)
g(5,n,j)=2.*g(5,n-1,j)-g(5,n-2,j)
g(8,n,j)=2.*g(8,n-1,j)-g(8,n-2,j)
end do
! Top boundary conditions, T=1.0
do i=0,n
g(8,i,m)=tw*(w(8)+w(6))-g(6,i,m)
g(7,i,m)=tw*(w(7)+w(5))-g(5,i,m)
g(4,i,m)=tw*(w(4)+w(2))-g(2,i,m)
g(1,i,m)=tw*(w(1)+w(3))-g(3,i,m)
end do
!Bottom boundary conditions
do i=0,n
g(2,i,0)=tw*(w(2)+w(4))-g(4,i,0)
g(6,i,0)=tw*(w(6)+w(8))-g(8,i,0)
g(5,i,0)=tw*(w(5)+w(7))-g(7,i,0)
end do
return
end
Эхний хэсэгт
оролтын хэсгийн хязгаарын нөхцлийг буцаж ойх нөхцлөөр тодорхойлсон байна.
Ингэхдээ сөрөг тэмдэг авах нь тооцоолол сүлжээний хурдны байгуулагчийн сөрөг
тэмдэгийг арилгах үүрэгтэй юм. Баруун гар тал буюу гаралтын хэсэг дээр өмнөхтэй
адил байдлаар нээлттэй хязгаарын нөхцлийг хэрэглэнэ. Хамгийн сүүлчийн сүлжээний түгэлтийн функцууд нь
сүүлээсээ хоёр ба гуравдах сүлжээний түгэлтийн функцыг ашиглан арифметик дундажаар
тодорхойлогдоно. Энэ нь гаралтын хэсэгт
температурын өөрчлөлт байхгүй dтөттө/dt=0 гэдэгтэй утга нэг нөхцөл юм. Дээд ба
доод хана нь тогтмол температураар шингэнийг халаах учир Диричлет /нэг хэмжээст
бодлогийг харах/ хязгаарын нөхцлийг ашиглана. Энэ хязгаарын нөхцөлд хязгаар
дахь цүлхэлтийн утга 0 байх ёстой гэсэнд үндэслэн дараах тэгшитгэлийг жишээ
болгон бичиж болно.
Энд байх
тэнцвэрт түгэлтийн функцуудыг дараах байдлаар тодорхойлж болно.
Эндээс g6 нь урсгалын процессоос тооцогдох боломжтой учир
хязгаарын нөхцөл g8-ийн хувьд
дараах байдалтай болно.
Дээрх
тэгшитгэлээр бичсэн код нь g(8,i,m)=tw*(w(8)+w(6))-g(6,i,m)
бөгөөд бусад үл мэдэгдэх g4, g7,
g1, доод хязгаарт g2, g6, g5 –
ийн
хувьд ижил байдлаар бодогдоно.
Температурын тархалтыг үнэлэх дэд програм subroutine
tcalcu
Хязгаарын
нөхцлийг өгч дууссаны дараа бодлогийн хүрээнд тодорхойлогдоогүй түгэлтийн функц
байхгүй болно. Иймд хугацааны алхамд харгалзах температурын тархалтыг бодох
боломж бий болно.
do j=1,m-1
do i=1,n-1
ssumt=0.0
do k=0,8
ssumt=ssumt+g(k,i,j)
end do
th(i,j)=ssumt
end do
end do
return
end
Дулааны тэгшитгэлд
бодогдсон түгэлтийн функцуудын нийлбэр макро орчин дахь хэмжээсгүй температурыг
тодорхойлно. Харин моментийн хэсэгт шингэний нягт ба хурдыг тодорхойлно.
Үр дүн
Оролтын хэсэг
дах температурын тархалт
Оролтын хэсэг
дахь босоо хурдны байгуулагч
Оролтын хэсэг дахь хэвтээ хурдны байгуулагч
Ламинар урсгалын төлөвших үе
Ингээд бүрэн кодыг хавсаргая.
ЛБА-аар халуун хоолой дотуурх ламинар ургалыг тооцох фортран код
1: program heatedchannel
2: parameter (n=1000,m=50)
3: real f(0:8,0:n,0:m)
4: real feq(0:8,0:n,0:m),rho(0:n,0:m)
5: real w(0:8), cx(0:8),cy(0:8)
6: real u(0:n,0:m), v(0:n,0:m)
7: real g(0:8,0:n,0:m), geq(0:8,0:n,0:m),th(0:n,0:m)
8: integer i
9: open(2,file='uvfield.data')
10: open(3,file='uvely.data')
11: open(4,file='vvelx.data')
12: !
13: cx(:)=(/0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0/)
14: cy(:)=(/0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0/)
15: w(:)=(/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./)
16: uo=0.1
17: rhoo=5.00
18: dx=1.0
19: dy=dx
20: dt=1.0
21: tw=1.0
22: th=0.0
23: g=0.0
24: visco=0.05
25: pr=3.8
26: alpha=visco/pr
27: Re=uo*m/alpha
28: print *, 'Re=', Re
29: omega=1.0/(3.*visco+0.5)
30: omegat=1.0/(3.*alpha+0.5)
31: mstep=10000
32: do j=0,m
33: do i=0,n
34: rho(i,j)=rhoo
35: u(i,j)=0.0
36: v(i,j)=0.0
37: end do
38: end do
39: do j=1,m-1
40: u(0,j)=uo
41: v(0,j)=0.0
42: end do
43: ! main loop
44: do kk=1,mstep
45: call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
46: call streaming(f,n,m)
47: call sfbound(f,n,m,uo)
48: call rhouv(f,rho,u,v,cx,cy,n,m)
49: do j=0,m
50: do i=0,n
51: sum=0.0
52: do k=0,8
53: sum=sum+g(k,i,j)
54: th(i,j)=sum
55: end do
56: end do
57: end do
58: ! collision for scalar
59: call collt(u,v,g,geq,th,omegat,w,cx,cy,n,m)
60: ! streaming for scalar
61: call streaming(g,n,m)
62: call gbound(g,tw,w,n,m)
63: print *, th(n/2,m/2),v(0,m/2),rho(0,m/2),u(n,m/2)
64: print*, v(n,m/2),rho(n,m/2)
65: END DO
66: ! end of the main loop
67: call result(u,v,th,uo,n,m)
68: stop
69: end program
70: ! end of the main program
71: subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
72: real f(0:8,0:n,0:m)
73: real feq(0:8,0:n,0:m),rho(0:n,0:m)
74: real w(0:8), cx(0:8),cy(0:8)
75: real u(0:n,0:m), v(0:n,0:m)
76: DO i=0,n
77: DO j=0,m
78: t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
79: DO k=0,8
80: t2=u(i,j)*cx(k)+v(i,j)*cy(k)
81: feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
82: f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
83: END DO
84: END DO
85: END DO
86: return
87: end
88: subroutine collt(u,v,g,geq,th,omegat,w,cx,cy,n,m)
89: real g(0:8,0:n,0:m),geq(0:8,0:n,0:m),th(0:n,0:m)
90: real w(0:8),cx(0:8),cy(0:8)
91: real u(0:n,0:m),v(0:n,0:m)
92: do i=0,n
93: do j=0,m
94: do k=0,8
95: geq(k,i,j)=th(i,j)*w(k)*(1.0+3.0*(u(i,j)*cx(k)+v(i,j)*cy(k)))
96: g(k,i,j)=omegat*geq(k,i,j)+(1.0-omegat)*g(k,i,j)
97: end do
98: end do
99: end do
100: return
101: end
102: subroutine streaming(f,n,m)
103: real f(0:8,0:n,0:m)
104: ! streaming
105: DO j=0,m
106: DO i=n,1,-1 !RIGHT TO LEFT
107: f(1,i,j)=f(1,i-1,j)
108: END DO
109: DO i=0,n-1 !LEFT TO RIGHT
110: f(3,i,j)=f(3,i+1,j)
111: END DO
112: END DO
113: DO j=m,1,-1 !TOP TO BOTTOM
114: DO i=0,n
115: f(2,i,j)=f(2,i,j-1)
116: END DO
117: DO i=n,1,-1
118: f(5,i,j)=f(5,i-1,j-1)
119: END DO
120: DO i=0,n-1
121: f(6,i,j)=f(6,i+1,j-1)
122: END DO
123: END DO
124: DO j=0,m-1 !BOTTOM TO TOP
125: DO i=0,n
126: f(4,i,j)=f(4,i,j+1)
127: END DO
128: DO i=0,n-1
129: f(7,i,j)=f(7,i+1,j+1)
130: END DO
131: DO i=n,1,-1
132: f(8,i,j)=f(8,i-1,j+1)
133: END DO
134: END DO
135: return
136: end
137: subroutine sfbound(f,n,m,uo)
138: real f(0:8,0:n,0:m)
139: do j=0,m
140: ! flow in on west boundary
141: rhow=(f(0,0,j)+f(2,0,j)+f(4,0,j)+2.*(f(3,0,j)+f(6,0,j)+f(7,0,j)))
142: rhow=rhow/(1.-uo)
143: f(1,0,j)=f(3,0,j)+2.*rhow*uo/3.
144: f(5,0,j)=f(7,0,j)+rhow*uo/6.
145: f(8,0,j)=f(6,0,j)+rhow*uo/6.
146: end do
147: ! bounce back on south boundary
148: do i=0,n
149: f(2,i,0)=f(4,i,0)
150: f(5,i,0)=f(7,i,0)
151: f(6,i,0)=f(8,i,0)
152: end do
153: ! bounce back, north boundary
154: do i=0,n
155: f(4,i,m)=f(2,i,m)
156: f(8,i,m)=f(6,i,m)
157: f(7,i,m)=f(5,i,m)
158: end do
159: ! account for open boundary condition at the outlet
160: do j=1,m
161: f(1,n,j)=2.*f(1,n-1,j)-f(1,n-2,j)
162: f(5,n,j)=2.*f(5,n-1,j)-f(5,n-2,j)
163: f(8,n,j)=2.*f(8,n-1,j)-f(8,n-2,j)
164: end do
165: return
166: end
167: subroutine gbound(g,tw,w,n,m)
168: real g(0:8,0:n,0:m)
169: real w(0:8)
170: ! Boundary conditions
171: ! West boundary condition, the temperature is given, tw
172: do j=0,m
173: g(1,0,j)=-g(3,0,j)
174: g(5,0,j)=-g(7,0,j)
175: g(8,0,j)=-g(6,0,j)
176: end do
177: ! Right hand boundary condition, open
178: do j=0,m
179: g(6,n,j)=2.*g(6,n-1,j)-g(6,n-2,j)
180: g(3,n,j)=2.*g(3,n-1,j)-g(3,n-2,j)
181: g(7,n,j)=2.*g(7,n-1,j)-g(7,n-2,j)
182: g(2,n,j)=2.*g(2,n-1,j)-g(2,n-2,j)
183: g(0,n,j)=2.*g(0,n-1,j)-g(0,n-2,j)
184: g(1,n,j)=2.*g(1,n-1,j)-g(1,n-2,j)
185: g(4,n,j)=2.*g(4,n-1,j)-g(4,n-2,j)
186: g(5,n,j)=2.*g(5,n-1,j)-g(5,n-2,j)
187: g(8,n,j)=2.*g(8,n-1,j)-g(8,n-2,j)
188: end do
189: ! Top boundary conditions, T=1.0
190: do i=0,n
191: g(8,i,m)=tw*(w(8)+w(6))-g(6,i,m)
192: g(7,i,m)=tw*(w(7)+w(5))-g(5,i,m)
193: g(4,i,m)=tw*(w(4)+w(2))-g(2,i,m)
194: g(1,i,m)=tw*(w(1)+w(3))-g(3,i,m)
195: end do
196: !Bottom boundary conditions, Adiabatic
197: ! g(1,i,0)=g(1,i,1)
198: ! g(2,i,0)=g(2,i,1)
199: ! g(3,i,0)=g(3,i,1)
200: ! g(4,i,0)=g(4,i,1)
201: ! g(5,i,0)=g(5,i,1)
202: ! g(6,i,0)=g(6,i,1)
203: ! g(7,i,0)=g(7,i,1)
204: ! g(8,i,0)=g(8,i,1)
205: ! T=0.0
206: do i=0,n
207: g(2,i,0)=tw*(w(2)+w(4))-g(4,i,0)
208: g(6,i,0)=tw*(w(6)+w(8))-g(8,i,0)
209: g(5,i,0)=tw*(w(5)+w(7))-g(7,i,0)
210: end do
211: return
212: end
213: subroutine tcalcu(g,th,n,m)
214: real g(0:8,0:n,0:m),th(0:n,0:m)
215: do j=1,m-1
216: do i=1,n-1
217: ssumt=0.0
218: do k=0,8
219: ssumt=ssumt+g(k,i,j)
220: end do
221: th(i,j)=ssumt
222: end do
223: end do
224: return
225: end
226: subroutine rhouv(f,rho,u,v,cx,cy,n,m)
227: real f(0:8,0:n,0:m),rho(0:n,0:m),u(0:n,0:m),v(0:n,0:m)
228: real cx(0:8),cy(0:8)
229: do j=0,m
230: do i=0,n
231: ssum=0.0
232: do k=0,8
233: ssum=ssum+f(k,i,j)
234: end do
235: rho(i,j)=ssum
236: end do
237: end do
238: DO i=1,n
239: DO j=1,m-1
240: usum=0.0
241: vsum=0.0
242: DO k=0,8
243: usum=usum+f(k,i,j)*cx(k)
244: vsum=vsum+f(k,i,j)*cy(k)
245: END DO
246: u(i,j)=usum/rho(i,j)
247: v(i,j)=vsum/rho(i,j)
248: END DO
249: END DO
250: do j=0,m
251: v(n,j)=0.0
252: end do
253: return
254: end
255: subroutine result(u,v,th,uo,n,m)
256: real u(0:n,0:m),v(0:n,0:m),th(0:n,0:m)
257: write(2,*)'VARIABLES =X, Y, U, V, T'
258: write(2,*)'ZONE ','I=',n+1,'J=',m+1,',','F=BLOCK'
259: do j=0,m
260: write(2,*)(i,i=0,n)
261: end do
262: do j=0,m
263: write(2,*)(j,i=0,n)
264: end do
265: do j=0,m
266: write(2,*)(u(i,j),i=0,n)
267: end do
268: do j=0,m
269: write(2,*)(v(i,j),i=0,n)
270: end do
271: do j=0,m
272: write(2,*)(th(i,j),i=0,n)
273: end do
274: do j=0,m
275: write(3,*)j/float(m),u(5,j)/uo,u(n/2,j)/uo,u(n-10,j)/uo
276: end do
277: do i=0,n
278: write(4,*) i/float(n),v(i,m/2)/uo
279: end do
280: return
281: end
No comments:
Post a Comment