Бүлэг 11 - Конвектив диффузын үзэгдэл гэж ярьдаг юм байна. Гэхдээ конвектив гэдэг нь монгол хэлний нэвчилт гэдэг үгэтэй дүйдэг гэнэ. Тэгэхээр Нэвчилт-Тархалтын үзэгдэл гэж ярьж бас болох л юм.

Конвекци-Тархалт /Convection-Diffusion problem/

11.1 Ууссан бодисын тээвэрлэлт

Бүлэг 4 – т ууссан бодисын тээвэрлэлтийг тархалтыг авч үзэлгүйгээр хэлэлцэж байсаныг санаж байгаа байх. Ууссан бодис (давс эвсэл хатуу хаягдал гэх мэт, мөн температур ижил асуудал болно.) нь дундаж хурдаар тээвэрлэгдэн мөн орон зайд тархаж байсан гэх мэт хэд хэдэн шалтгааны улмаас дээрхи тээвэрлэлтийн тооцоолол нь бодит биш мэт санагдаж болно. Өөрөөр хэлбэл голын хөндлөн огтлолын дагуу хурдны хуваарилалт (энэ нь маш чухал үзүүлэлт) өөр өөр байна мөн турбулент холилдолт, молекулын тархалт гэх мэт зүйлсийг авч үзээгүй гэсэн үг юм. Илүү дэлгэрэнгүй судалгааг Fischer et al (1979) – ээс уншаарай. Ууссан бодисын тархалтыг физик процессын үүднээс төсөөлөн бодож үзвэл энэ нь концентраци ихтэй газраас багатай газарлуу дээрх процессуудын улмаас тээвэрлэгдэх ёстой. Үүнийг зураг 11.1 – т харуулав. Үүнийг бид энд тархалт (diffusion) гэж нэрлэж байгаа юм. Хурдны өөрчлөлтөөр тархаж буй үед, дисперсийн илэрхийлэл ихэвчлэн хэрэглэгдэнэ. Дундаж хурдаар тээвэрлэгдэж процессыг конвекци гэх ба нийлүүлбэл конвекци-тархалтын үзэгдэл нэгэн зэрэг явагддаг гэж ойлгож болно.

Зарим хөндлөн огтлолын хувьд тухайн огтлолд хөндлөн ууссан бодисын тээвэрлэгдэх ба энэ нь дараах хоёр хэсгийг агуулна.

Үүнд: А – хөндлөн огтлолын талбай, u – дундаж хурд, c – хөндлөн огтлол дахь дундаж концентраци болно.

Тархалтын коэффициент D нь нэг тохиолдлоос нөгөөд хүчтэй өөрчлөгдөж болно. Ердийн байгалийн голд үүнийг дараах байдлаар тодорхойлж болно.

Энд а – нь усны гүн ба k – нь туршилгын хүчин зүйл байх ба энэ нь 50 – аас 500 хооронд өөрчлөгдөж байдаг.

Ууссан бодисын массын балансын тэгшитгэл нь дараах байдалтай байна.

Тэгшитгэл 11.1 – д үүнийг орлуулбал:

Дифференциалчилж, усны балансын тэгшитгэл 4.2 – ийг тооцоонд авч, A болон D – г тогтмол гэж үзсэнээр конвекци-тархалтын тэгшитгэлийн оргиналь хэлбэрийг гарган авах болно.

Энэ нь энгийн долгионы тэгшитэл (4.3) болон тархалтын тэгшитгэл (7.4) – ийг нэгтгэн авч үзсэн тэгшитгэл мэт харагдаж байгаа боловч математик үүднээсээ өөрийн гэсэн шинж чанартай байдаг. Жишээлбэл бид энэ тэгшитгэлд z=x-ut гэсэн хувьсах хэмжигдэхүүнийг авч үзвэл тархалтын тэгшитгэлийг гарган авах болно.

Тийм болохоор, конвекци-тархалтын тэгшитгэл нь маш үр дүнтэйгээр харьцангуй хүрээн дотор дундаж урсгалтай хөдөлж байгаа тархалтыг илэрхийлдэг. Бид шаардлагатай адил тооны болон адил төрлийн анхны болон хязгаарын нөхцлүүдийг авч үзэх хэрэгтэй. Жишээлбэл

-          Бүсийн цэг бүрт нэг анхны нөхцөл

-          Хугацааны функ эсвэл хязгаар бүрт нэг хязгаарын нөхцөл гэх мэт

11.2 Тооцолох арга

Конвекци-тархалтын тэгшитгэлийн тоон шийдэл нь цэвэр тархалтын тэгшитгэлийн ойролцоолох аргатай маш ижил. Зургаан-цэгт “молекул” –тай адил торны цэгийг ашилан 7.3 дугаар хэсэгт үзсэн шиг Кранк-Ничолсоны аргаар тэгшитгэлийн ойролцооллыг тэнцвэрт байдлаар бичиж болно.

Хэрэв thetha=0 бол энэ арга нь FTCS төрлийн ил арга болно. Эсрэгээр энэ нь далд арга юм. Сүүлчийн тохиолдолд ялгаварын тэгшитгэлийн систем нь хязгаарын нөхцөлтэйгээ цуг 7.4 хэсэгт үзсэн шиг Томасын алгоритмоор бодогдоно.

Энэхүү жишээнд, энэ аргын нэмэгдэлтийн хүчин зүйл нь дараах байдлаар тодорхйологдоно.

Энд мэдэгдэж буй параметрүүд нь

Тархалтын параметр lamda=2Ddelta t/delta x²

Коурантын тоо sigma=u delta t/delta x байна.

Жишээнд, өмнөх шиг 1/2.le.thetha.le.1 байхад нөхцөлт бус тогтвортой арга болохыг харж болно. Түүгээрч барахгүй ил арга нь дараах хоёр нөхцлийн дор тогтвортой байна.

Ил аргын тухайд, тэгшитгэл 11.9-д байх хоёр нөхцөл нь илүү хязгаарлалттай байгааг харах нь сонирхолтой юм. Өгөгдсөн торны хэмжээ delta x үед баруун гар талын тэнцэтгэл нь дараах байдалтай болно.

Энд бид нэгжийн Рейнолдсын тоо R нь торны тэгш бус байдал дээр ямар нэгэн юм илэрхийлж байгааг харж болно.

11.3 Хэрэглээ

Практик талдаа тохиолд, тэгшитгэл 11.5 дахь авч үзэж буй коэффициент нь ихэвчлэн тогтмол биш байдаг. Тиймээс энэ нь хоёр үндсэн тэгшитгэл болох 11.1 болон 11.2 – ийг шууд гаргалгаа хийхэд зөвлүүштэй зүйл юм. Зураг 11.2 – т байгаа торыг хэрэглэн дараах байдлаар гүйцэтгэж болно.

Тээвэрлэлт s нь торны цэгийн хооронд тодорхойлогдоно. Ингээд массын балансын тэгшитгэлийг дараах байдлаар гаргалгааг хийж болно.

Тээвэрлэлт нь төвийн ялгаварын тэгшитгэлийг хэрэглэн дараах байдлаар тооцогдоно.

Тогтмол коэффициенттэй үед, үндсэн тэгшитгэлүүд нь Кранк-Ничолсоны аргаруу хөтлөнө.

Тогтмол дундаж урсац нь 0.002м/с байх голдирол дахь давсны тоо хэмжээний гаргалтын авч үзье. Тархалтын коэффициент нь D=1.6м²/с гэж тодорхойлогдов. Тодорхой хугацаанд 1500м – ийн урттай давсны “бүрхүүл” үүссэн болохыг харна. Хэрэв доод хашицруу 10км – т очвол энэ нь ямар харагдах бол? Энэ нь дараах өгөгдөлтэй үед Кранк-Ничолсоны арга дээр суурилсан тооцон бодох загварын симулци юм. Thetha=0.55, delta x=250м, delta t=5000s

Үр дагаварт нь, sigma=0.4 ба lamda=0.256 байна. Тэгэхээр ил арга нь сайн ажилласан гэсэн үг. Зураг 11.3 – д симуляцын үр дүнг харуулах ба эндээс хамгийн их концентраци нь өөрийнхөө оргиналь утгын хагас орчим хүртэл буурсаныг тодорхойлж болно. Хэрэв тархалтыг үл тооцсон бол оргил утга нь бууралтгүйгээр хэвээр үлдэх байсан.

Ингээд тооцоолон бодох Фортран кодийг сонирхоё.

!

!    陰解法（クランクニコルソン法）を用いた１次元移流拡散方程式の解法

!

      parameter(pi=3.141592, g=9.8              )

      parameter(nt=1000,dt=1.0 )

      parameter(peri=nt*dt/5.  )

      parameter(nx=100,dx=5.0  )

      parameter(theta=0.6) !クランクニコルソン {күранкүникорүсан}гэж хираганаагаар бичжээ.

      parameter(difc=2.0       )

      parameter(uamp=1.0       )

      dimension a(nx),b(nx),c(nx),d(nx),ans(nx)

      dimension c1(nx),u(nx)

      open(10,file='result.dat')

      do i=1,nx

           c1(i)=0.0

          if(i.ge.nx/2-5.and.i.le.nx/2+5) then

           c1(i)=1.0

          endif

      enddo

        do n=1,nt

        time=n*dt

        vel=uamp*sin(2.*pi*time/peri)

        do i=1,nx

        u(i)=vel

        enddo

            do i=1,nx

        if(i.eq.1) then

         b(i)= 1.0

         c(i)=-1.0

         d(i)= 0.0

        else if(i.eq.nx) then

         a(i)=-1.0

         b(i)= 1.0

         d(i)= 0.0

        else

          uu=(u(i)+u(i-1))/2.

        alfa=(uu+abs(uu))/(2.*dx)

        beta=(uu-abs(uu))/(2.*dx)

        a(i)=   (-dt*(alfa+difc/dx**2)           )*theta

        b(i)=1.+(((alfa-beta)+2.*difc/dx**2)*dt  )*theta

        c(i)=   ( dt*(beta-difc/dx**2)           )*theta

        xx=dt*(alfa*(c1(i)-c1(i-1))+beta*(c1(i+1)-c1(i)))*(1.-theta)

        yy=dt*(difc*(c1(i+1)-2.*c1(i)+c1(i-1))/dx**2    )*(1.-theta)

        d(i)= c1(i)-xx+yy

        endif

            enddo

            call tridag(a,b,c,d,ans,nx,1,1)

        c1=ans

        if(mod(n,2).eq.0) then

        write(10,'(100f8.3)') (c1(i),i=1,nx)

        endif

        enddo ! time

        close(10)

            stop

            end

Энд та бүхэн Call tridag гэсэн мөн байгааг харж байгаа байх. Энэ нь тусдаа байгаа дэд програмыг дуудаж оруулж ирэх комманд юм. Ихэнх програм нь бидний үзэж байгаагаас хэд дахин том, хэдэн зуун мөр агуулсан байдаг учир тэдгээрт хэрэглэгдэх математик илэрхийлэл, томъёолол зэргийг дэд програмын тусламжтайгаар өөр файл болгон хадгалж болдог. Энэ нь алдааг хайхад тусад нь шүүх боломж олгодог, өөр төрлийн програмд шууд ашиглах боломжтой гэх мэт олон талын ач холбогдолтой байдаг. Энэ удаагын дэд програмаар Гурван диагональт nxn хэмжээст матрицийг Гауссын зайлуулах аргаар бодох Томасын алгоритмын кодыг авч үзнэ.

      SUBROUTINE TRIDAG(A,B,C,R,U,N,II,JJ)

      PARAMETER (NMAX=500)

      DIMENSION GAM(NMAX),A(N),B(N),C(N),R(N),U(N)

      IF(B(1).EQ.0.)PAUSE

      BET=B(1)

      U(1)=R(1)/BET

      DO 11 J=2,N

        GAM(J)=C(J-1)/BET

        BET=B(J)-A(J)*GAM(J)

        IF(BET.EQ.0.) THEN

              write(*,*) II,JJ,A(J),B(J),C(J),GAM(J),BET,J,N

              PAUSE

              ENDIF

        U(J)=(R(J)-A(J)*U(J-1))/BET

11    CONTINUE

      DO 12 J=N-1,1,-1

        U(J)=U(J)-GAM(J+1)*U(J+1)

12    CONTINUE

      RETURN

      END

Хэрэв таны бодлого эсвэл тоон загварчлал Гурван диагональт матрицаар
шийдэгдэхээр бол (Кранк-Ничолсоны арга) энэ дэд програмын кодыг шууд авч ямар ч өөрчлөлт хийлгүйгээр ашиглах боломжтой гэсэн үг юм.

Ингээд Конвекци-Тархалтын тоон шийдлийн үр дүнг харъя.

Figure 1. Тоон шийдлийн үр дүн. Концентраци хугацааны туршид хэрхэн тархаж, конвекцлогдож байгааг харж болно.