Нэг хэмжээст дулаан дамжуулалт /1D heat conduction/

Нэг хэмжээст дулаан дамжуулалтын загвар

Нэг хэмжээст хатуу хавтангаар дулаан дамжих үзэгдлийг загварчлая. Хавтангийн анхны температурыг T=0.0 байна гэж үзье. Хугацаа t >0 үед хавтангийн зүүн хязгаараас T_w=1.0 гэсэн температур өгөгдөнө. Хавтангийн уртыг 100 нэгж гэж үзье. Хавтангийн дулаан тархалтын коэффициент гэвал хавтангаар дамжих температурыг t=200 нэгж хугацаанд СБА болон ТЯА-аар тооцож үр дүнгүүдийг харьцуулъя.

Бодох тэгшитгэл

Энд Ф тархаж буй эсвэл дамжуулагдаж буй дулаан, алфа нь орчны дамжуулалтын коэффициент болно. Үүнийг төгсгөлөг ялгаварын аргаар ойролцоолбол:

Сүлжээний Больцманны арга болон Төгсгөлөг ялгаварын аргаар бодсон үр дүнг харьцуулбал:

1 хэмжээст хавтгай орчинд дулаан дамжих үзэгдлийг СБА болон ТЯА-аар бодсон үр дүн

Хоёр аргын фортран эх кодуудыг тус бүрд нь өгөөд СБА-ын бодолтыг тайлбарлая. Эхлээд СБА-ын кодыг үзье.

Одоо төгсгөлөг ялгаварын аргын кодыг үзье.
Сүлжээний Больцманны аргын кодны ажиллагааг тайлбарлая. 

        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

Ингээд бодолт дуусаж үр дүнг хэвлэнэ. Энэ аргаар бодсон үр дүнг Төгсгөлөг ялгаварын аргаар бодсон үр дүнтэй харьцуулж харуулбал

Энэ зургийг хамгийн эхэнд харсан ш дээ. Зүүн талаас 1 гэсэн температур өгөхөд тархалт нь шугаман бус байгааг харж болно. 

Тооцоонд бараг ялгаа гарсангүй. Учир нь тархалтын тэгшитгэлийн төгсгөлөг ялгаварын тэгшитгэл болон ЛБА-ын мөргөлдөөний тэгшитгэлүүд нь бараг адилтгал тэгшитгэлүүд болж таардаг. Нюфлотоор хэрхэн дээрх графикийг байгаалах кодыг авч үзье.

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'