Saturday, February 25, 2017

Стефаны бодлого /Stefan problem solved/

Стефаны бодлого 

Математик болон физикийн шийд нь олдсон үндсэн, сонгодог бодлогууд ихэвчлэн шийдийг нь олсон эсвэл их хувь нэмэр оруулсан эрдэмтэнийхээ нэрээр нэрлэгддэг. Тухайлбал бидний өмнө шийдвэрлэсэн Стокесийн бодлогууд, одоо шийдэх гэж буй Стефаны бодлого, арван жил, их сургуульд таарч байсан Диричлетийн бодлого, Кошийн бодлого, физикийн алдарт Ньютон-Пеписийн бодлого, математикийн Васелийн бодлого, Гауссын дугуйн бодлого, геометрт Жанжин* Напелионы бодлого - хоолонд хүртэл нэрээ мөнхөлсөн Напелион шүү дээ, магадлалын онол дахь төрсөн өдрийн бодлого, бүр Черилийн төрсөн өдрийн бодлого гэх мэт хүн болгон бараг нэг нэг бодлоготой байгаа бол зарим хүмүүс олонтой, жишээ нь Хилбертийн 23 бодлого гэж байх юм. Амьдал дээр ч гэсэн бид бүгд асуудалтай учирч түүнийгээ яаж шийдвэрлэх вэ гэдгээ боддог. Тэр бодлого маань нийтлэг чанартай бол аятайхан математикаар илэрхийлчихвэл тусгай бодлоготой болох боломж байна аа. Цахим ертөнцөд хүний нэрээр нэрлэгдсэн тэгшитгэл, хууль, физик үзэгдэл, коэффициент, хэмжээсгүй тоонуудын жагсаалт гэж байхаас хүний нэрээр нэрлэгдсэн стандарт сонгодог бодлогуудын жагсаалт байдаггүй юм байна. За энэ ч яахав Стефаны бодлогоо тодорхойлъё.
Сонгодог Стефаны бодлого нь нэгэн төрлийн орчинд бий болж байгаа фазын өөрчлөлт, уг орчин дахь температурын тархалт зэргийг тодорхойлдог бодлого юм.
Стефаны бодлогын үр дүн. Дараа нь гоё хувилбарыг нь Нюфлот кодтой нь оруулна.
Температурын тархалт нь дулааны тэгшитгэлээр илэрхийлэгдэх ба фазийн шилжилт нь уг нэгэн төрлийн орчинд хатуу-шингэн-хийн аль нэг хоёр нь харилцсан харилцах үеийг бий болгоно. Тухайлбал бид мөс хайлах үзэгдлийг Стефаны бодлогоор шийдвэрлэх ба фазын өөрчлөлтийн үр дүнд бий болон харилцах үе нь шингэн-хатуу харилцах үе байна. Энэ харилцах үеийн байрлалыг хугацаанаас хамааруулан тодорхойлох, мөс болон усны орчинд температур ямар маягаар тархаж байгааг Стефаны бодлого бодож өгнө гэсэн үг юм. Бид тэгэхээр чөлөөт хязгаартай /шингэн-хатуу харилцах үе хугацаанаас хамаарч өөрчлөгдөнө/, дифференциал тэгшитгэлийг бодох хэрэгтэй гэсэн үг. Бодлогын математикаар илэрхийлбэл:


Стефаны бодлогын аналитик шийд
Энэ бодлогыг анх 1831 оны үед Ламе болон Клапэирон гэх хүмүүс дурдаж байсан бол Стефан 1890 онд томъёолж аналитик шийдийг нь олсон байна. Бодлогын аналитик шийд нь шингэн-хатуу харилцах үеийн байрлал, температурын тархалтыг
гэсэн хоёр томъёогоор олж болох ба томъёонд байгаа кси хэмээх параметрийг
гэх тодорхойгүй функцаас олно гэж байгаа. Хамгийн сонирхолтой нь erf гэх алдааны функц болон тодорхойгүй функцаас кси-г яаж олох гэдэг л байгаа юм. Алдааны функц нь алдааны тархалтыг харуулах сигмойд хэлбэрийн функц ба тархалтын тэгшитэглийг тайлбарлахад, эсвэл ститикстик магадлалд хэмжилтийн алдааг тодорхойлоход хэрэглэгддэг байна. Ингээд 10 см урттай мөсийг нэг талаас нь 25 градусын хэмтэй халуунаар халааж, нөгөө үзүүрийг нь хөлдөх 0 градусын температурт барих юм. Дээрх томъёонуудаар бодохын тулд 10 см мөсөө 100 торонд хувааж 1сек -ийн хугацааны алхамтайгаар бодохоор Фортран хэлээр код бичсэнийг доор харуулна. Тэгшитгэл 28-д байгаа альфа нь усны дулаан тархалтын коэффициент, t нь мэдээж ахиж буй хугацаа, харин тэгшитгэл 29-д т макс т мэлт энэ тэр нь 25 градус, 0 градус болж байгаа юм. хүртвэрт буй алдааны функцын ард хаалтанд х нь зайн алхам, хуваарьт байгаа алдааны функцын хаалтанд кси нь параметрүүд болно. Функц 30-д байгаа St -г Стефаны тоо гэх ба
гэж олно. Үүнд с нь усны дулаан багтаамж, L нь нуугдмал дулаан буюу хайлахын хувийн дулаан зэрэг болно. 
За тодорхойгүй функцыг анх хараад яаж бодноо гэж айдас хүрээгүй нэг шалтгаан бол Инженерийн гидравлик үзэж байхдаа Шезийн коэффициент, критик гүн гэх мэтийг дөхөх аргаар яг иймэрхүү тодорхойгүй функцаас олдог байсан юм. Харин энэ удаад арай өөр дөхөх арга ашиглана. Энэ бол Ньютон-Рапсоны арга юм. Энэ аргаар язгуур олоход их амархан. Ньютон-Рапсоны аргын шийд олох тэгшитгэл бол
юм. Энд функцыг тухайн функцын уламжлалд харьцуулж анхны гөрдөж өгсөн утгаас хасаж дараагийн тохирол нь сайжирсан шийдийг олох явдал юм. Хэдэн удаа итераци хийн тэр чинээгээр тохирол сайжирна гэсэн үг юм. Функц 30-ийн нэгдүгээр эрэмбийн уламлал бас сонирхолтой. Хувьсагч нь кси бөгөөд тэнцүүгийн тэмдэгийн өмнө гэхэд гурван функцын үржвэр, энэ нь гурван функцын үржвэрийн дүрмийг уламжлалдаа ашиглана гэсэн үг. Амьдралдаа дээд тал нь 2 функцын үржвэрийн уламжлал авч байсан бол энэ удаагых гайхалтай тохиолдол байлаа. Ингээд гурван функцын үржвэрийн уламжлалыг яаж авахыг санавал:
Харин мань алдааны функц өөрөө нэг
иймэрхүү уламжлал өгдөг юм байна. Ерөнхийдөө алдааны функцыг Эйлэрийн тооны зэрэгт функцаар төсөөлж болно. 
ингээд 30 функцын нэгдүгээр эрэмбийн уламжлал нь
болж байгаа юм. Та бас өөрийгөө шалгаад үзээрэй. 
За ингээд гол асуудал нь алдааны функцыг яаж бодох вэ гэдэг байгаа юм. Фортран хэл дээр бол янз бүрийн хэцүү функцуудыг бодох дэд програмууд маш их олдоцтой байдаг. Зарим өргөн дэлгэр ашиглагддаг функцууд нь үндсэн функц болоод хэлэнд програмчлагдсан байдаг. Алдааны функц бол тэдний нэг бөгөөд erf гэж бичээд хаалтанд тоогоо оруулангуут утгыг бодоод өгнө. Гэхдээ бусад хэлнүүдэд жишээ нь Пайтон, Руби эсвэл Жавад энэ үндсэн функцаар өгөгдөөгүй бол яаж вэ? Фортртаны 77-оос өмнөх хувилбарт ч гэсэн алдааны функцыг бодох үндсэн функц байхгүй байжээ. Тэр үед алдааны функцыг ингэж боджээ:
Алдааны функцыг бодох эхний Фортран хувилба - энийг Проф. Хосояамада бичиж өгсөн юм.
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
! test program for erfnc

      integer err
      real x,fx
      write(*,600) ' Test results from erfnc !'
      write(*,*  ) '  x          erf(x)'

      do 10 x=0.,5.,0.5
      call erfnc(x,fx,err)
      write(*,610) x,fx
   10 continue
  600 format(a/)
  610 format(f6.1,f15.6)
      end


      subroutine erfnc(x,fx,err)
      integer i,err
      real x,fx
      real*8 a(6),w,wx
      data a/1.00000000d0,  .0705230784d0, .0422820123d0,
     1        .92705272d-2, .1520143d-3  , .2765672d-3   /

      if(x.lt.0.) then
       err=999
      else
       err=0
       if(x.eq.0.) then
        fx=0.
       elseif(x.ge.10.0) then
        fx=1.
       else
        wx=x
        w=.430638d-4
        do 10 i=6,1,-1
         w=w*wx+a(i)
   10   continue
        fx=1.d0-1.d0/w**16
       endif
      endif

      return
      end

Харин хоёр дахь хувилбар нь Басик хэлнээс Фортранруу орчуулсан хуучин код юм.

 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
      program write

      write(*,*) erf(1.)

      stop
      end



      real function erf(x)

      data a1/0.0705230784/
      data a2/0.0422820123/
      data a3/9.27052e-3  /
      data a4/1.52014e-4  /
      data a5/2.76567e-4  /
      data a6/4.30638e-5  /

      y=abs(x)

      apollo=1.+y*(a1+y*(a2+y*(a3+y*(a4+y*(a5+y*a6)))))
      saturn=1.-(1./apollo)**((1./apollo)*16.)*16.

      if(x.eq.y) then
       erf=+saturn
      else
       erf=-saturn
      endif

      return
      end

Эдгээрийг  алдааны функцыг шууд бодох боломжгүй хэлнүүдрүү хөрвүүлээд авч болох нь байна. Харин миний хувьд фортраны erf() гэдэг үндсэн функцаа ашигласан юм. 

 Кси-г Ньютон-Рапсоны аргаар олох. 
Стефаны тоо мэдэгдэж байвал кси гуайг Ньютон-Рапсоны аргаар дараах байдлаар олж болно. 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
      subroutine liamda (St,lamda)
      parameter (pi=3.1415926)
      parameter (iter=20)
      real lamda

      lamda=0.1 !initial guess
      do i=1,iter
      funct=St/sqrt(pi)-lamda*exp(lamda**2)*erf(lamda)
      derivative=-exp(lamda**2)*erf(lamda)-2.*lamda**2
     &           *exp(lamda**2)*erf(lamda)-lamda*2./sqrt(pi)
      lamda=lamda-funct/derivative
      end do

      return
      end

Энд funct нь функц өөрөө, derivative гэдэг нь дээр хэдүүлээ гаргасан нэгдүгээр эрэмбийн уламжлал юм. За кси (би кодонд лиамда гэж бичсэн байгаа) -г олчихсон юм чинь Стефаны бодлогыг хэрхэн кодолж аналитик аргаар шийд олсоныг харъя. 

 Стефаны бодлогын аналитик шийд олох код

 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
! stefan problem analytical solution from Book
!"Mathematical modeling of Melting and freezing process"
! Code by Ayurzana.B
      program stefan

      parameter (m=100, t=200, dis=0.1 ) !grid size, time, slab (m)
      parameter (tm=0., th=25.0, tc=0.) !melt,hot,cold
      parameter (tdl=0.143e-6, tds=1.1819e-6) !thermal Dif (m2/s-1)
      parameter (La=333.4,  pi=3.1415926 ) !kJ/kg and pi
      parameter (chl=4.1868, chs=0.138) !specific heat (J/gCels)
      real tem(t,m), xs(t) ! temperature and melt front
      real dx,dt,time,lamda   ! spacing, time, parameter
      real St
      character*6 rr
   
      open(1,file='interface.dat')
      open(2,file='temp2.dat')
      open(7,file='flowan.dat')
   
      St=chl*(th-tm)/La
      !lamda=0.706*sqrt(St)*(1.-0.21*(0.5642*St)**(0.93-0.15*St))
      call liamda(St,lamda)
      t_t=dis**2/(4.*lamda**2*tdl)
      dt=t_t/t
      dx=dis/m
   
      t_p1=0.01**2/(4.*lamda**2*tdl) !1cm
      t_p3=0.03**2/(4.*lamda**2*tdl) !3cm
      t_p5=0.05**2/(4.*lamda**2*tdl) !5cm
      t_p9=0.09**2/(4.*lamda**2*tdl) !9cm
   
      print*, t_t/3600, 50*dt/3600, 100*dt/3600, 150*dt/3600, 200*dt/3600
      xs=0.
      tem=0.
   
      do i=1,t
         time=i*dt
         xs(i)=2.*lamda*sqrt(tdl*time) ! free boundary
         do j=1,m
              if(j*dx.le.xs(i)) then
            tem(i,j)=th+(tm-th)*erf((j*dx)/(sqrt(4.*tdl*time)))
     &                         /erf(lamda)
              else
!         if(j*dx.ge.xs(i)) tem(i,j)=tm
            tem(i,j)=tm*erf((j*dx)/(sqrt(4.*tdl*time)))
     &                         /erf(lamda)
              end if
  
         end do
   
         write(1,*) i, xs(i)/dx !time/3600.
         write(7,*) (i*dt)/3600., tem(i,90)    


      end do
   
      do i=1,t
      write(2,200) (tem(i,j), j=1,m)
      enddo
  200 format(1000e11.3)
   
      print*, St, lamda, t_p1, t_p3, t_p5, t_t
      
   
      close (1)
      close (2)
      close (7)
      stop
      end program stefan
Гол ажил 36-47-р мөрөнд явагдаж байгаа юм. 38-д тэгшитгэл 28-ыг, 41 ба 45-д 29-ээр тэгшитгэлийг бодож байгаа юм. 

 Сүлжээний Больцманы аргаар Стефаны бодлогыг бодох
Одоо сүлжээний Больцманы аргыг фазийн шилжилттэй бодлогыг хэрхэн бодож чадах эсэхийг нь аналитик шийдтэй харьцуулж шалгах үлдэж байна. 1 хэмжээст бодлого учир D1Q3 загвар ашиглана. Бодолтын алгоримт бол ерөнхийдөө Стокесын 2-р бодлоготой төстэй. Энэ дурьдмаар байгаа гол зүйл бол сүлжээний тэнцвэрт орох хугацааг
гэж тодорхойлно гэдгийг үзүүлэх юм. Энд хугацааны алхам, торны алхам энэ тэр оролцсон байгаа. 
10 см мөсний халж байгаа үзүүр дээр Диричлетийн хязгаарын нөхцөл, нөгөө үзүүр дээр нь хоёрдугаар эрэмбийн экстраполяцийн хязгаарын нөхцлийг өгнө. Интерполяци их боддог байсан бол экстраполяци амархаан. За тооцон бодох фортран кодыг сонжвол:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
! stefan problem lattice Boltzmann solution
! Code by Ayurzana.B
      program OneDLBM
      parameter (m=100, dis=0.1 ) !grid size, time, slab (m)
      parameter (tm=0., th=25.0, tc=0.)  !melt,hot,cold
      parameter (tdl=0.143e-6, tds=1.1819e-6) !thermal Dif (m2/s-1)
      parameter (La=333.4,  pi=3.1415926 )    !kJ/kg and pi
      parameter (chl=4.1868, chs=0.138) !specific heat (J/gCels)
      parameter (rh=1.0, twc=0.2)        !pi, temperature unit
      real f(0:2,m), w(0:2), cy(0:2)    !distribution functions, LB
      real tem(m),treal(m)     ! temp and melt front
      real dx,dt,time,lamda    ! spacing, time, parameter
      real St, lf(m), xs(200000), En(m)
      real L_La,L_tdl,L_tds
      character*6 rr
   
      open(1,file='result.dat')
      open(12,file='flowlb.dat') !
      open(13,file='track.dat')
      open(14,file='interface.dat')

      adj=0.2
      St=chl*(th-tm)/La
!      lamda=0.5*sqrt(St)*(1.-0.21*(0.5642*St)**(0.93-0.15*St))
      call liamda(St,lamda)
      t_t=dis**2/(4.*lamda**2*tdl)
      dt=1.0!0.2
      nt=int(t_t/dt)
      dx=dis/m
      L_La=La*dt**2/dx**2 !lattice latent heat
      cp_w=St*L_La        !1.2 lattice water specific heat
      cp_i=(cp_w*chs)/chl
      ome=1.0/(tdl*dt/dx**2+0.5)+adj    !1.9-1.0 relaxation parameter
      L_tdl=tdl*dt/dx**2   !lattice thermal diff of water
      L_tds=tds*dt/dx**2   !lattice thermal diff of ice

      print*, nt, L_La, cp_w, ome, St, L_tdl, L_tds
      npr=nt/200
      pause
      w(:)=(/4./6.,1./6.,1./6./)
      cy(:)=(/0.0,1.0,-1.0/)
!initial conditions 
        xs=0.     
        do j=1,m
          lf(j)=0.0
          tem(j)=0.0
          En(j)=0.0
          treal(j)=0.
          do k=0,2
             f(k,j)=w(k)*tem(j)
          end do   
        end do

! main loop
      do i=1,nt
         time=i*dt
! streaming
         do j=m,1,-1
           f(1,j)=f(1,j-1)
           f(2,m-j)=f(2,m-j+1)
         end do
! boundary condition
         f(2,m)=2,0*f(2,m-2)-f(2,m-1)!
         f(1,1)=rh*(w(1)+w(2))-f(2,1)
! macroscopic variables
         do j=1,m
           tem(j)=(f(1,j)+f(2,j)+f(0,j))
           if(tem(j).lt.0.) tem(j)=0.
           if(lf(j).eq.1.0) then
           treal(j)=(tem(j)-twc)/((rh-twc)/th)
           else
           treal(j)=0.0
           end if
         end do
! phase change
         do j=1,m
           cp=(1.-lf(j))*cp_i+lf(j)*cp_w
           En(j)=cp*tem(j)

      if (En(j).gt.cp*twc) lf(j)=1.0  !1.25e-4
      if (En(j).lt.cp*twc) lf(j)=0.0       !2.50e-4

      if(lf(j)+lf(j+1).eq.1.0) xs(i)=j*dx
         end do
! collision
         do j=1,m

           f(0,j)=(1.-ome)*f(0,j)+ome*w(0)*tem(j)
           f(1,j)=(1.-ome)*f(1,j)+ome*w(1)*tem(j)
           f(2,j)=(1.-ome)*f(2,j)+ome*w(2)*tem(j)
!                                     equalibrium DFs
         end do
         if(mod(i,npr).eq.0) write(1,'(2000f8.3)') (treal(j),j=1,m)
         if(mod(i,npr).eq.0) write(13,'(2000f8.3)') (lf(j),j=1,m)

      
      if(mod(i,npr).eq.0) write(14,*) i/npr, xs(i)/dx

      end do
   
      close (1)
      close (12)
      close (13)
      close (14)
      end
Жаахан урт ч гэсэн үндсэн үйлдлүүд нь 54-95-р мөрөнд л хийгдэж байгаа юм. Хэрэг болж магадгүй гээд кодуудыг бүтнээр ямар нэгэн хасалтгүйгээр оруулж байгаа юм. 

 За үр дүнгүүдээ харъя. 
Босоо тэнхлэгт хугацаа, хэвтээд 10 см мөсний урт байгаа юм. Улаан талаас температур 25 градусын хэм хүртэл халж мөс хайлах ба мөс ба усны харилцах үе хар ба улаан зураасаар байна гэж тооцоологдсон байна. 

Энэ графикт цэнхэрээр хэд хэдэн хугацаан дахь мөсний уртын дагуу температурын тархалтыг харуулж байгаа юм. Тэдгээр хугацаа нь тус бүр 1см, 3, 5, 9см - т мөс хайлж хүрэх үеийн хугацаанууд юм. Эдгээр байрлалуудад хугацааны турш температур хэрхэн өөрчлөгдөж вэ гэдгийг улаан графикаар харуулсан байна. 
Ньюфлотоор дээрх хоёр графикийг хэрхэн байгуулж зурахыг дараах нийтлэлээс олж мэдэх болно. Тооцон бодох математик, тооцон бодох шингэний динамиктай холбоотой сонгодог бодлогууд цаашид шийдэгдэн үргэлжлэх байх гэж бодож байна. 

 Аналитик шийдний бүрэн кодийг GitHub
Сүлжээний Болцманы аргын бүрэн кодыг GitHub
боловсруулсан: Усны барилгын инженер Б.Аюурзана
Нийтэлсэн at
Labels: , , ,

1 comment:

  1. EnkhtaivanMarch 26, 2017 at 3:22 AM

    Их л сайхан тайлбарлаж бичиж. Хэрэгтэй хүн нь уншиж хэрэглээсэй.

    ReplyDelete

Subscribe to: Post Comments (Atom)